# 問題定義 定義一個$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模擬實現 首先我們定義一些常規的泡利矩陣： ```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$： ```python H0 = sigmax print (H0) print (np.linalg.eig(H0)[0]) print (np.linalg.eig(H0)[1]) ``` 這段python程式碼的執行結果如下： ```bash [[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$： ```python H1 = (sigmax + sigmaz) / np.sqrt(2) print (H1) print (np.linalg.eig(H1)[0]) print (np.linalg.eig(H1)[1]) ``` 輸出結果如下： ```bash [[ 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`中的庫函式來實現的，不是通過絕熱演化演算法來實現的，這裡的資料我們將作為一個對標的物件，測試最終絕熱演化結果的準確度。除了上述定義的哈密頓矩陣之外，我們還需要定義一個常用的量子力學操作：歸一化。 ```python 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`作為唯一的引數： ```python def annealing_solver(steps): t = 0 eg_vector1 = np.array([1, 1]) / np.sqrt(2) eg_value1 = 1 energy1 = [eg_value1] Ht = H0 hbar = 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) / hbar + 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=(1-\frac{t}{T})H_0+\frac{t}{T}H_1$$ 這個哈密頓量形式的意義在於，當$t=0$時，$H_t=H_0$，當$t=1$時，$H_t=H_1$。只要`steps`=$\frac{T}{t}$足夠大，也就是絕熱演化的時間足夠長時，就可以演化到最終需要求解的本徵態上。此時再代入薛定諤方程並將其寫成可數值求解的差分形式： $$i\hbar(\left|\psi(s+\frac{t}{T})\right>-\left|\psi(s)\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： ```python 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() ``` 執行結果如下： ```bash [0.92387953 0.38268343] [0.92387953 0.38268343] ```  由於目標本徵態所對應的本徵能量比初始的本徵能量高，因此隨著迭代次數的增加，中間能量值也在逐步上升，並最終達到期望的本徵值。同時對比最終的本徵態向量我們可以發現其實是一致的，只是通過絕熱演化找到的是另外一組正交基，但同樣也是一組合法的本徵態向量。我們再考察一下，設定不同的`steps`會對結果的計算精度產生什麼樣的影響： ```python 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)$$ 讓我們來用另外一組的本徵態來測試該迭代方程的求解效果： ```python 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() ``` 得到的結果迭代過程圖如下所示，我們可以看到同樣是可以演化到目標本徵態的： ```bash [0.38268343 0.92387953] [0.38268343 0.92387953] ```  # 總結概要 根據絕熱近似我們可以數值模擬絕熱演化的過程，進而求解得到目標哈密頓矩陣的本徵態。而量子退火也是基於絕熱演化的原理，在實際的物理體系上通過控制物理溫度和位元之間的耦合使得最終系統演化到一個能量最低的本徵態，結果的驗證可以自洽。 # 參考資料 1. *Quantum Computation by Adiabatic Evolution*, Edward Farhi and his collaborators, 28 Jan 2000, arXiv:quantum-ph/