Max-Cut問題を解こう
長い冬が終わり、梅の花が咲き始め、春の陽気漂うこの時期にふと解きたくなる問題がある。
そう、「Max-Cut問題」だ。
（※ 個人差がありますし、別に自分は解きたいとは思いません）
Max-Cut問題はグラフの頂点を２つのグループに分ける時に、グループ間の辺の本数が最大になるような分け方を考える問題。
量子コンピュータ用のフレームワークであるQiskitでは各辺に重みをつけたグラフを描き、Max-Cut問題を作ることができる。
import matplotlib.pyplot as plt
import numpy as np
import networkx as nx
n=5 # ノードの数
G=nx.Graph()
G.add_nodes_from(np.arange(0,n,1))
elist=[(0,1,1.0),(0,3,1.0),(1,2,1.0),(2,3,1.0),(2,4,1.0),(3,4,1.0)]
# tuple is (i,j,weight) where (i,j) is the edge
G.add_weighted_edges_from(elist)
colors = ['r' for node in G.nodes()]
pos = nx.spring_layout(G)
def draw_graph(G, colors, pos):
default_axes = plt.axes(frameon=True)
nx.draw_networkx(G, node_color=colors, node_size=600, alpha=.8, ax=default_axes, pos=pos)
edge_labels = nx.get_edge_attributes(G, 'weight')
nx.draw_networkx_edge_labels(G, pos=pos, edge_labels=edge_labels)
draw_graph(G, colors, pos)
さらに上記のMax-Cut問題は各辺の重みを使って以下の行列で表現することができる。
w = np.zeros([n,n])
for i in range(n):
for j in range(n):
temp = G.get_edge_data(i,j,default=0)
if temp != 0:
w[i,j] = temp['weight']
print(w)
# 結果
[[0. 1. 0. 1. 0.]
[1. 0. 1. 0. 0.]
[0. 1. 0. 1. 1.]
[1. 0. 1. 0. 1.]
[0. 0. 1. 1. 0.]]
行列のi行j列の値は、i番目のノードとj番目へのノードをつなぐ辺の重みであり、辺で結ばれていないものや同じノード同士は０となる。
この行列がいわゆるQUBO（二次非制約二項最適化問題）であり、イジングモデルの問題を解く際に必要な情報になる。
以下ではこのQUBOとQiskitのツールを用いてMax-Cut問題を解いていく。
Max-Cut用のモジュール 「max_cut」
Qiskitにはその用途に合わせてTerra、Ignis、Aer、Aquaという4つのモジュール群が用意されているが（FF4のゴルベーザ四天王みたい！...FFやったことないけど）、今回は量子コンピュータ用のアルゴリズムを提供するAquaを利用する。
イジングモデルの問題はAquaが提供するモジュールの１つであるqiskit.optimization.applications.isingからmax_cutをインポートして解くことができる。
この他にもqiskit.optimization.applications.isingには巡回セールスマン問題やナップサップ問題など数理最適化問題ごとのイジングモデルが用意されているため、自身が解きたい問題テーマに合わせてイジングモデルを選択する。
from qiskit.optimization.applications.ising import max_cut
qubitOp, offset = max_cut.get_operator(w) # QUBOを引数として渡す
print('Offset:', offset)
print('Ising Hamiltonian:')
print(qubitOp.print_details())
max_cut.get_operator(w)は引数としてQUBOを渡すことで対象のMax-Cut問題のハミルトニアンを作る関数。
返り値はWeightedPauliOperator型の以下の結果を返す。
フムフム………ナニソレ？？？
# 結果
Offset: -3.0
Ising Hamiltonian:
IIIZZ (0.5+0j)
IIZZI (0.5+0j)
IZIIZ (0.5+0j)
IZZII (0.5+0j)
ZIZII (0.5+0j)
ZZIII (0.5+0j)
急によくわからない数字と昇竜拳のコマンドみたいなものと0.5の羅列が出てきた。
WeightedPauliOperator型・・・？
0と1のみを変数として使うQUBOではMax-Cut問題は以下の数式の最小化問題として記述される。…①式
\begin{align}
E(x)&=\sum (-x_i-x_j+2x_ix_j)\\
\end{align}
0の代わりに-1を用いるイジングモデルでは、
x_i = \frac{s_i+1}{2}
として式を再記述する必要がある。
（QUBOの1がイジングの1、QUBOの0がイジングの-1に対応）
| QUBO ($x_i$) | イジング ($s_i$) |
|---|---|
| 1 | 1 |
| 0 | -1 |
したがって、①式は
\begin{align}
E(s)&=\sum(-\frac{s_i+1}{2}-\frac{s_j+1}{2}+2\frac{s_i+1}{2}\frac{s_j+1}{2})\\
&=\sum(-\frac{1}{2}s_i-\frac{1}{2}s_j-1+\frac{1}{2}s_is_j+\frac{1}{2}s_i+\frac{1}{2}s_j+\frac{1}{2})\\
&=-\frac{1}{2}\sum(1-s_is_j)
\end{align}
となる。…②式
これを初めに作ったMax-Cut問題に沿って展開すると以下のようになる。…③式
\begin{align}
E(s)&=-\frac{1}{2}\sum(1-s_is_j)\\
&=0.5s_0s_1+0.5s_0s_3+0.5s_1s_2+0.5s_2s_3+0.5s_2s_4+0.5s_3s_4-3.0
\end{align}
ん？なんかそれっぽいものが出てきた。
# 結果
Offset: -3.0
Ising Hamiltonian:
IIIZZ (0.5+0j)
IIZZI (0.5+0j)
IZIIZ (0.5+0j)
IZZII (0.5+0j)
ZIZII (0.5+0j)
ZZIII (0.5+0j)
max_cut.get_operator(w)によって、WeightedPauliOperator型として返却された結果を再掲する。
ここで示されてるOffset: -3.0は③式の定数部分-3.0に該当する箇所を意味しており、WeightedPauliOperatorのパラメータのうち、atolに該当するものだ。
(0.5+0j)は③式の$s_is_j$にかかる係数である。（+0jは虚数部を示すため今回は虚数部は０、つまり実数のみ）
そしてIやZはパウリゲート（IゲートおよびZゲート）であり、上記の係数と合わせてWeightedPauliOperatorのpaulisというパラメータに該当する。
このパラメータは「重み付きのパウリゲート」のリストであり、重みとパウリゲートのペアを配列として持っている。
すなわちIIIZZ (0.5+0j)は0.5の重みをもつIゲート３つとZゲート２つから成る量子回路であり、Zゲートが各$s_i$に対応する形となっている。
（例えばIIIZZは$s_0s_1$、IIZZIは$s_1s_2$。なんで右からになるんだろうとかいうのはまだよく分かってない...）
細かい話をすればもっと深くなるのだろうが、とりあえず急に現れたよくわからない結果に対して、それっぽい意味付けができた。
このあとも少し見慣れない関数が出てくるが以下のようにすればMax-Cut問題の解を得ることができる。
qp = QuadraticProgram()
qp.from_ising(qubitOp, offset)
qp.to_docplex().prettyprint()
exact = MinimumEigenOptimizer(NumPyMinimumEigensolver())
result = exact.solve(qp)
print(result)
# 結果
optimal function value: -5.0
optimal value: [0. 1. 0. 1. 0.]
status: SUCCESS
Max-Cut問題の答えは[0. 1. 0. 1. 0.]、すなわち１番目と３番目のノードが同じグループになるようにグラフをカットしてあげるのが最適解。めでたしめでたし。
参考