# 最大切割問題介紹 最大切割問題（Max-Cut），也常作為最小切割問題（Min-Cut）出現，這兩個問題可以等價，只需要對權重值取負號即可。給定一個無向加權圖$G(V,E)$，找到一個方案將所有的節點$\{V\}$劃分為兩組$\{V_1\}$和$\{V_2\}$，使得這兩組點之間所連線的邊的權重之和最大（如果是最小切割問題就是權重和最小）。讓我們看一個實際的問題圖，該問題中一共包含4個節點：$v_1, v_2, v_3, v_4$，以及4條邊：$e_{1,2}, e_{2,3}, e_{3,4}, e_{2,4}$，為了簡化問題，我們令這四條邊的權重值都為1，即：$w_{1,2}=w_{2,3}=w_{3,4}=w_{2,4}=1$，該問題所對應的問題圖如下所示： ```python import networkx as nx import matplotlib.pyplot as plt G = nx.Graph() G.add_nodes_from([1, 2, 3, 4]) G.add_edges_from([(1, 2), (2, 3), (3, 4), (2, 4)]) plt.subplot(121) nx.draw(G, with_labels=True) ```  這個給定的示例問題較為簡單，我們容易看出，將給定的節點$\{v_1,v_2,v_3,v_4\}$化為$\{v_2\}\cup\{v_1,v_3,v_4\}$的時候，切割的邊集合為$\{e_{1,2},e_{2,3},e_{2,4}\}$，即切割的邊的權重和為3。因為在該問題中找不到能夠切割4條邊及以上的節點集合劃分方法，因此我們認為$\{e_{1,2},e_{2,3},e_{2,4}\}$的切割方案是問題的一個解（有時候最優解不一定是唯一解）。 # 量子絕熱演算法與量子退火 在上一篇文章中介紹了[使用絕熱演化/量子退火演算法求解矩陣本徵態](https://www.cnblogs.com/dechinphy/p/annealer.html)，這裡僅作簡單總結。量子絕熱演算法是基於準靜態物理過程，利用演化過程中保持本徵態的特性對目標哈密頓量的本徵態進行求解，其對應的薛定諤方程為： $$i\hbar\frac{d}{dt}\left|\psi\right>=H\left|\psi\right>$$ 以及： $$H\left|\psi\right>=E\left|\psi\right>$$ 假定初始哈密頓量為$H_0$以及目標哈密頓量為$H_1$，則過程中對應的哈密頓量為： $$H=(1-s)H_0+sH_1$$ 此處$s$表示演化時長。那麼我們可以瞭解到，在最終系統演化到$H_1$時，對應的系統狀態就是$H_1$的本徵態。 量子退火演算法是基於絕熱演化的原理而構造的一套基於退火機的組合優化問題求解方案，可以將初始態設定為一個本徵能量較高的狀態，最終通過緩慢降溫使得系統演化到另一個目標哈密頓量的基態。在數學上來說，只要演化的時間足夠長，該演算法可以保障一定能夠給出目標哈密頓量的最優解。 # 最大切割問題的Ising建模 最大切割問題的一個特點是，僅需要考慮任意兩點之間的連線關係，因此我們可以採用`Ising模型`對最大切割問題進行建模。首先我們考慮Ising模型的一般形式： $$H_{Ising}=\sum^{N}_{i=1}h_is_i+\sum^{N}_{i=1}\sum^{N}_{j=i+1}J_{i,j}s_is_j$$ 其中$s\in\{-1,1\}$用於表示物理體系的自旋向上與自旋向下，在處理最大切割問題時，可以作為處在兩個不同節點集合$\{V_1\}$和$\{V_2\}$的標記。 由於最大切割問題中不涉及到節點本身的一些屬性（物理學中可以稱之為偏移量），所以最大切割問題中的$E_{Ising}$中的第一項為`0`。這樣一來，對應到最大切割問題的無向加權圖，我們可以重新寫出Max-Cut問題的哈密頓量： $$H_{MaxCut}= -(\sigma^z_1\sigma^z_2+\sigma^z_2\sigma^z_3+\sigma^z_3\sigma^z_4+\sigma^z_2\sigma^z_4)$$ 由於這裡取的是最大切割，我們可以想象，如果是$s_i=s_j$，則結果一定是正的，這表示$i,j$兩個節點被分配到了同一個集合$V_1$或者$V_2$裡面。那麼如果切割的邊越多，我們希望它的能量值是越大的，因此需要取負號才能使得趨勢是正確的。*注：由於作者一直專注在程式碼實現層面，對於最大切割問題為什麼選取$\sigma^z$作為對自旋的操作量還沒有一個獨立的思考，只是從文章中直接摘錄，後續再補充理論解釋。* # 絕熱演化求解哈密頓量基態 首先，我們還是需要定義好泡利矩陣的形式： ```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) ``` 然後定義一個容易製備和求解本徵態的哈密頓量$H_0$及其本徵態向量，這裡用$\sigma^x$來構造： ```python H0 = np.kron(np.kron(np.kron(sigmax,sigmax),sigmax),sigmax) eigen_vector1 = np.array([1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]).T ``` 接下來根據最大切割問題定義一個目標哈密頓量$H_1=H_{MaxCut}$： ```python H1 = - np.kron(sigmaz, np.kron(sigmaz, np.kron(sigmai, sigmai))) - np.kron(sigmai, np.kron(sigmaz, np.kron(sigmaz, sigmai))) - np.kron(sigmai, np.kron(sigmaz, np.kron(sigmai, sigmaz))) - np.kron(sigmai, np.kron(sigmai, np.kron(sigmaz, sigmaz))) ``` 我們可以先用其他的計算方法計算該哈密頓量所對應的本徵態作為一個對比： ```python print (np.linalg.eig(H1)) ``` 得到的結果如下： ```bash (array([-4.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 2.-0.j, 2.-0.j, 2.-0.j, -2.-0.j, -2.-0.j, 2.-0.j, 2.-0.j, 2.-0.j, 0.+0.j, 0.+0.j, 0.+0.j, -4.+0.j]), array([[1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j], [0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j], [0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j], [0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j], [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j], [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j], [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j], [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j], [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j], [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j], [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j], [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j], [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j], [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j], [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j], [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j]])) ``` 由於中間計算過程的需要，我們定義一個量子力學常用操作`歸一化`函式： ```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 ``` 接下來我們還是按照薛定諤方程定義一個絕熱演化的函式，將演化的步長作為一個入參： ```python def annealing_solver(steps): t = 0 eg_vector1 = np.array([0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]).T eg_value1 = 1 energy1 = [eg_value1] Ht = H0 hbar = 6.62607015e-34 for s in range(steps): t += 1 / steps eg_vector1 = uniform(np.dot(Ht, eg_vector1) * (-1j) * (np.pi * 2 / steps) / hbar + eg_vector1) Ht = (1 - t) * H0 + t * H1 print (np.abs(uniform(eg_vector1))) return np.abs(uniform(eg_vector1)) ``` 在該演化函式中，我們已經定義了一個初始哈密頓量的本徵態用於演化，需要注意的是，並不是所有的初始態都能夠演化到目標哈密頓量的基態。用組合優化的語言來說就是，迭代的結果不一定能夠得到最優解，但是一般都可以得到一個較優解。讓我們先來看一下這個絕熱演化所得到的結果： ```python from matplotlib import pyplot as plt eg_vector = annealing_solver(100)[1] fig,ax=plt.subplots() ax.bar(range(16), eg_vector) plt.show() ``` 這裡我們可以得到哈密頓量演化的最終量子態分佈： ```bash [0. 0. 0. 0. 0.5 0.5 0. 0. 0. 0. 0.5 0.5 0. 0. 0. 0. ] ```  需要注意的是，這裡得到的結果並不是最大切割問題的最終解，我們只是得到了包含最終解的一個分佈。我們先逐一看下這個分佈中僅有的4個量子態所對應的切割結果。第五個位置的量子態的二進位制表示為：`0100`，對應的集合劃分為：$\{v_1,v_3,v_4\}\cup\{v_2\}$，這個解對應的邊的切割數為3,是理論的最優解。第六個位置對應的二進位制表示為：`0101`，對應的集合劃分為：$\{v_1,v_3\}\cup\{v_2,v_4\}$，這個解對應的邊的切割數為3,也是理論的最優解。同理，第11和第12個位置是對稱的結構，都是理論最優解。因此，我們到這裡就完整的利用量子絕熱演算法/量子退火演算法解決了一個最大切割問題，並得到了兩組不同的最優解。 # 技術彩蛋 雖然，上述章節的結果已經找到了兩組最優解，但是，這還並不是所有的理論最優解。這裡我們做了一個嘗試，把初始的本徵態設定為： ```python eg_vector1 = np.array([0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0]).T ``` 這樣以來我們得到的最終量子態就變成了： ```bash [0. 0. 0. 0. 0. 0. 0.5 0.5 0.5 0.5 0. 0. 0. 0. 0. 0. ] ```  我們從這個結果中可以發現，$\{v_1,v_4\}\cup\{v_2,v_3\}$也是一組理論最優解。從這裡我們可以得到一個結論，絕熱演化的結果具有侷限性，不能保障找到所有的最優解。但是在實際問題中，往往我們能夠找到一組最優解就已經足夠了。 # 參考連結 1. Edward Farhi, Jeffrey Goldstone, and Sam Gutmann, “*A quantum approximate optimization algorithm*” (2014), arXiv:141