量子絕熱算法求解最大切割問題

最大切割問題介紹

最大切割問題（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\)，該問題所對應的問題圖如下所示：

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}\}\)的切割方案是問題的一個解（有時候最優解不一定是唯一解）。

量子絕熱算法與量子退火

在上一篇文章中介紹了使用絕熱演化/量子退火算法求解矩陣本征態，這里僅作簡單總結。量子絕熱算法是基於准靜態物理過程，利用演化過程中保持本征態的特性對目標哈密頓量的本征態進行求解，其對應的薛定諤方程為：

\[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\)作為對自旋的操作量還沒有一個獨立的思考，只是從文章中直接摘錄，后續再補充理論解釋。

絕熱演化求解哈密頓量基態

首先，我們還是需要定義好泡利矩陣的形式：

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\)來構造：

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}\)：

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)))

我們可以先用其他的計算方法計算該哈密頓量所對應的本征態作為一個對比：

print (np.linalg.eig(H1))

得到的結果如下：

(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]]))

由於中間計算過程的需要，我們定義一個量子力學常用操作歸一化函數：

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

接下來我們還是按照薛定諤方程定義一個絕熱演化的函數，將演化的步長作為一個入參：

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))

在該演化函數中，我們已經定義了一個初始哈密頓量的本征態用於演化，需要注意的是，並不是所有的初始態都能夠演化到目標哈密頓量的基態。用組合優化的語言來說就是，迭代的結果不一定能夠得到最優解，但是一般都可以得到一個較優解。讓我們先來看一下這個絕熱演化所得到的結果：

from matplotlib import pyplot as plt
eg_vector = annealing_solver(100)[1]
fig,ax=plt.subplots()
ax.bar(range(16), eg_vector)
plt.show()

這里我們可以得到哈密頓量演化的最終量子態分布：

[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個位置是對稱的結構，都是理論最優解。因此，我們到這里就完整的利用量子絕熱算法/量子退火算法解決了一個最大切割問題，並得到了兩組不同的最優解。

技術彩蛋

雖然，上述章節的結果已經找到了兩組最優解，但是，這還並不是所有的理論最優解。這里我們做了一個嘗試，把初始的本征態設置為：

eg_vector1 = np.array([0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0]).T

這樣以來我們得到的最終量子態就變成了：

[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\}\)也是一組理論最優解。從這里我們可以得到一個結論，絕熱演化的結果具有局限性，不能保障找到所有的最優解。但是在實際問題中，往往我們能夠找到一組最優解就已經足夠了。

版權聲明

本文首發鏈接為：https://www.cnblogs.com/dechinphy/p/qaa-maxcut.html
作者ID：DechinPhy
更多原著文章請參考：https://www.cnblogs.com/dechinphy/

參考鏈接

1. Edward Farhi, Jeffrey Goldstone, and Sam Gutmann, “A quantum approximate optimization algorithm” (2014), arXiv:1411.4028.