使用絕熱演化/量子退火算法求解矩陣本征態

問題定義

定義一個\(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.