当前位置 博文首页 > 陆言君：量子绝热算法求解最大切割问题
最大切割问题（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}\}\)的切割方案是问题的一个解（有时候最优解不一定是唯一解）。
在上一篇文章中介绍了使用绝热演化/量子退火算法求解矩阵本征态，这里仅作简单总结。量子绝热算法是基于准静态物理过程，利用演化过程中保持本征态的特性对目标哈密顿量的本征态进行求解，其对应的薛定谔方程为：
以及：
假定初始哈密顿量为\(H_0\)以及目标哈密顿量为\(H_1\)，则过程中对应的哈密顿量为：
此处\(s\)表示演化时长。那么我们可以了解到，在最终系统演化到\(H_1\)时，对应的系统状态就是\(H_1\)的本征态。
量子退火算法是基于绝热演化的原理而构造的一套基于退火机的组合优化问题求解方案，可以将初始态设置为一个本征能量较高的状态，最终通过缓慢降温使得系统演化到另一个目标哈密顿量的基态。在数学上来说，只要演化的时间足够长，该算法可以保障一定能够给出目标哈密顿量的最优解。
最大切割问题的一个特点是，仅需要考虑任意两点之间的连接关系，因此我们可以采用Ising模型对最大切割问题进行建模。首先我们考虑Ising模型的一般形式：
其中\(s\in\{-1,1\}\)用于表示物理体系的自旋向上与自旋向下，在处理最大切割问题时，可以作为处在两个不同节点集合\(\{V_1\}\)和\(\{V_2\}\)的标记。
由于最大切割问题中不涉及到节点本身的一些属性（物理学中可以称之为偏移量），所以最大切割问题中的\(E_{Ising}\)中的第一项为0。这样一来，对应到最大切割问题的无向加权图，我们可以重新写出Max-Cut问题的哈密顿量：
由于这里取的是最大切割，我们可以想象，如果是\(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\}\)也是一组理论最优解。从这里我们可以得到一个结论，绝热演化的结果具有局限性，不能保障找到所有的最优解。但是在实际问题中，往往我们能够找到一组最优解就已经足够了。
