使用绝热演化/量子退火算法求解矩阵本征态

问题定义

定义一个\(N\times N\)大小的矩阵\(H\)，找到该矩阵的本征态。已知：若态矢量\(\left|\psi\right>\)为哈密顿矩阵\(H\)的本征矢，则有：

\[H\left|\psi(t)\right>=E\left|\psi(t)\right> \]

此处\(E\)为哈密顿矩阵\(H\)的本征能量，或称为本征值。

绝热演化与量子退火

绝热演化过程可以这么理解，在求解一个已知哈密顿矩阵\(H_1\)的本征态时，先制备一个容易计算出本征态的哈密顿矩阵\(H_0\)所对应的物理系统，并使得该物理系统出于对应的本征态\(\left|\psi(0)\right>\)。根据绝热近似，如果我们设计一条准静态的演化路径，使得系统哈密顿矩阵从\(H_0\)逐渐演化到\(H_1\)，此时可以测量的系统状态正对应一个\(H_1\)的本征态。这就相当于，我们利用一个物理系统的绝热演化过程，完成了一个矩阵本征态问题的求解。

基于绝热演化的原理，我们可以假想这样的一个场景：假如将我们所常见的一些问题如物流规划、频谱分配等组合优化问题，转换成QUBO模型对应成一个哈密顿矩阵，这样一来我们就可以通过对物理系统的操作，来实现实际问题的求解。D-wave这个公司就以此为出发点，发明了量子退火机，并且已经初步实现了其商业价值。

量子退火，实际上就是利用了绝热演化的原理：通过调制超导比特之间的耦合关系和对每个比特的控制，先制备一个本征能量较高的超导物理系统，然后精准控制物理温度缓慢降温，就可以实现到目标哈密顿矩阵的绝热演化。由于组合优化的求解过程是自洽的，因此可以根据前后两次不同温度所对应的能量值来判断是否需要继续演化，这使得量子退火机可以在既定的时间内找到一个极优解。注意这里不一定是最优解，但是根据量子退火理论，如果演化的时间足够长，理论上退火过程必然会演化到最优解，但是对于大部分的实际问题而言，极优解已经足够了。

绝热演化/量子退火算法Python模拟实现

首先我们定义一些常规的泡利矩阵：

import numpy as np
sigmai = np.array([[1, 0], [0, 1]], dtype = complex)
sigmax = np.array([[0, 1], [1, 0]], dtype = complex)
sigmay = np.array([[0, -1j], [1j, 0]], dtype = complex)
sigmaz = np.array([[1, 0], [0, -1]], dtype = complex)

这些泡利矩阵具有非常好的物理性质，并且可以构造任意的\(2\times 2\)矩阵，这里不作展开介绍。上面提到我们需要定义一个容易得到本征态的初始哈密顿矩阵\(H_0\)：

H0 = sigmax
print (H0)
print (np.linalg.eig(H0)[0])
print (np.linalg.eig(H0)[1])

这段python代码的执行结果如下：

[[0.+0.j 1.+0.j]
 [1.+0.j 0.+0.j]]
[ 1.+0.j -1.+0.j]
[[ 0.70710678-0.j  0.70710678+0.j]
 [ 0.70710678+0.j -0.70710678+0.j]]

这里打印的结果分别表示：\(H_0\)的矩阵形式、\(H_0\)的两个本征能量以及\(H_0\)的两个本征态矢量。这里我们可以以同样的方式来定义我们的目标哈密顿矩阵\(H_1\)：

H1 = (sigmax + sigmaz) / np.sqrt(2)
print (H1)
print (np.linalg.eig(H1)[0])
print (np.linalg.eig(H1)[1])

输出结果如下：

[[ 0.70710678+0.j  0.70710678+0.j]
 [ 0.70710678+0.j -0.70710678+0.j]]
[ 1.+0.j -1.+0.j]
[[ 0.92387953+0.j -0.38268343+0.j]
 [ 0.38268343-0.j  0.92387953+0.j]]

需要注意的是，这里计算出来的本征能量和本征态矢量是通过numpy中的库函数来实现的，不是通过绝热演化算法来实现的，这里的数据我们将作为一个对标的对象，测试最终绝热演化结果的准确度。除了上述定义的哈密顿矩阵之外，我们还需要定义一个常用的量子力学操作：归一化。

def uniform(state):
    s = 0
    for i in state:
        s += abs(i) ** 2
    for i in range(len(state)):
        state[i] /= np.sqrt(s)
    return state

定义好这些内容之后，我们开始正式定义一个绝热演化的函数，并将绝热演化过程的步骤数量steps作为唯一的参数：

def annealing_solver(steps):
    t = 0
    eg_vector1 = np.array([1, 1]) / np.sqrt(2)
    eg_value1 = 1
    energy1 = [eg_value1]
    Ht = H0
    h = 6.62607015e-34
    for s in range(steps):
        t += 1 / steps
        eg_vector_tmp1 = uniform(np.dot(Ht, eg_vector1) * (-1j) * (math.pi * 2 / steps) / h + eg_vector1)
        Ht = (1 - t) * H0 + t * H1
        eg_value1 = np.abs(eg_vector_tmp1[0]) * eg_value1 / np.abs(eg_vector1[0])
        eg_vector1 = eg_vector_tmp1
        energy1.append(eg_value1)
    print (np.abs(uniform(eg_vector1)))
    return energy1, uniform(eg_vector1)

由于此处涉及到的物理内容较多，这里简单展开说明一下。首先根据薛定谔方程有：

\[i\hbar\frac{d}{dt}\left|\psi(t)\right>=H_t\left|\psi(t)\right> \]

而绝热演化过程中的哈密顿量有：

\[H_t=\left(1-\frac{t}{T}\right)H_0+\frac{t}{T}H_1 \]

这个哈密顿量形式的意义在于，当\(t=0\)时，\(H_t=H_0\)，当\(t=1\)时，\(H_t=H_1\)。只要steps足够大，也就是绝热演化的时间足够长时，就可以演化到最终需要求解的本征态上。此时再代入薛定谔方程并将其写成可数值求解的差分形式：

\[i\hbar\left(\frac{\left|\psi(s+\frac{t}{T})\right>-\left|\psi(s)\right>}{1/steps}\right)=((1-s)H_0+sH_1)\left|\psi(s)\right> \]

则不断的计算\(\left|\psi(s+\frac{t}{T})\right>\)最后即可得到\(\left|\psi_1\right>=\left|\psi(\frac{T}{T})\right>\)为最终的目标本征态。我们可以用上面的这个例子来测试一下，首先我们尝试将steps设置为100：

import matplotlib.pyplot as plt
plt.figure()
energy1 = annealing_solver(100)[0]
print (np.abs(np.linalg.eig(H1)[1][0]))
plt.plot(energy1)
plt.show()

执行结果如下：

[0.92387953 0.38268343]
[0.92387953 0.38268343]

由于目标本征态所对应的本征能量比初始的本征能量高，因此随着迭代次数的增加，中间能量值也在逐步上升，并最终达到期望的本征值。同时对比最终的本征态矢量我们可以发现其实是一致的，只是通过绝热演化找到的是另外一组正交基，但同样也是一组合法的本征态矢量。

迭代步长对结果精度的影响

我们再考察一下，设置不同的steps会对结果的计算精度产生什么样的影响：

import matplotlib.pyplot as plt
plt.figure()
error1 = []
phase = 0.38268343 / 0.92387953
for i in range(49, 1000, 20):
    vector1 = annealing_solver(i)[1]
    error1.append(np.abs(vector1[1] / vector1[0]) - phase)
plt.plot(range(49, 1000, 20), error1)
plt.plot(range(49, 1000, 20), 1e-03 * np.ones(len(range(49, 1000, 20))), '-.')
plt.show()

其对应的结果图如下所示：

在组合优化常规问题中，并未声明对求解精度的要求，在其他领域中一般的精度要求在\(1\times 10^{-3}\)，所以我们这里也标识了要达到这个期望精度所需要的演化要求。

基于本征能量特点的另一种实现方案

在最前面我们提到过一个公式：\(H\left|\psi\right>=E\left|\psi\right>\)，这个公式在\(\left|\psi\right>\)是\(H\)的本征态时成立，\(E\)是一个常量。而我们前面考虑的绝热演化过程中，每一步得到的状态都是一组对应于当下哈密顿量的本征态，因此我们可以得到另外一种形式的迭代方程：

\[uniform(\left|\psi(s)\right>)=uniform\left(((1-s)H_0+sH_1)\left|\psi(s-\frac{t}{T})\right>\right) \]

让我们来用另外一组的本征态来测试该迭代方程的求解效果：

def annealing_solver(steps):
    t = 0
    eg_vector2 = np.array([1, -1]) / np.sqrt(2)
    eg_value2 = -1
    energy2 = [eg_value2]
    for s in range(steps):
        t += 1 / steps
        Ht = (1 - t) * H0 + t * H1
        eg_vector_tmp2 = uniform(np.dot(Ht, eg_vector2))
        eg_value2 = np.abs(eg_vector_tmp2[0]) * eg_value2 / np.abs(eg_vector2[0])
        eg_vector2 = eg_vector_tmp2
        energy2.append(eg_value2)
    print (np.abs(uniform(eg_vector2)))
    return energy2, uniform(eg_vector2)
if __name__ == '__main__':
    import matplotlib.pyplot as plt
    plt.figure()
    energy2 = annealing_solver(100)[0]
    print (np.abs(np.linalg.eig(H1)[1][1]))
    plt.plot(energy2)
    plt.show()

得到的结果迭代过程图如下所示，我们可以看到同样是可以演化到目标本征态的：

[0.38268343 0.92387953]
[0.38268343 0.92387953]

总结概要

根据绝热近似我们可以数值模拟绝热演化的过程，进而求解得到目标哈密顿矩阵的本征态。而量子退火也是基于绝热演化的原理，在实际的物理体系上通过控制物理温度和比特之间的耦合使得最终系统演化到一个能量最低的本征态，结果的验证可以自洽。

版权声明

本文首发链接为：https://www.cnblogs.com/dechinphy/p/annealer.html
作者ID：DechinPhy
更多原著文章请参考：https://www.cnblogs.com/dechinphy/

参考资料

1. Quantum Computation by Adiabatic Evolution, Edward Farhi and his collaborators, 28 Jan 2000, arXiv:quantum-ph/0001106.