論文arxiv1302.5843 実装その２

1. 本記事の概要

NPな色々な問題に対してIsingで定式化した論文arxiv1302.5843があります。
D-waveの勉強がてら、それを実装してみようというもの。

2. 本記事の想定読者

• 最適化がなんとなくわかっている
• 理屈はいいから使えればいいや
• Pythonなんとなくわかる

レベルの人

3. 対象の問題

「2.2. Graph Partitioning」を対象とする。

頂点の数 $N=|V|$が偶数である無向グラフ $ G=(V,E) $ において、頂点数が同数の$N/2$ となるようグラフを分割した際に、エッジ数（辺）が最小となる分割方法は？

4. Ising式

論文を見ると2つのパーツに分かれていますね。頂点が同数となるような制限事項$H_A$と辺のCut方法$H_B$からなるようです。
すなわち、
$$ H_A = A(\sum_{n=1}^N s_i)^2 　$$

$$ H_B = B(\sum_{(uv) \in E} \frac{1-s_us_v}{2}) 　$$

だそうです。（Qiitaはきれいな数式が書けるので気持ちいですね！）

さて、グダグダ言わず、実装に移りましょう。（ほんと、Pyqubo超絶便利です。というか私、Pyquboなくしては何もできないんですが、いいのか！？）

問題生成

%matplotlib inline
import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import random
from pyqubo import Placeholder,Array,Constraint,Sum,solve_qubo

SIZE = 8
G = nx.Graph()

# サイズ分のVertixを作成し、適当にEdgeをつなぐ
for i in range(SIZE):
    G.add_node(i)
    if i != 0:
        for j in range(random.randint(1,6)):
            New_line = random.randint(0,i-1)
            if (G.get_edge_data(i,New_line) == None):
                G.add_edge(i,New_line,W=1)

# グラフの出来具合を確認する
pos = nx.spring_layout(G, k=5.,seed=0)
nx.draw_networkx(G,pos)
plt.rcParams['figure.figsize'] = (7.0, 7.0)
plt.show()

今回生成された問題は以下のようなもの。

SAで解を確認

# ハミルトニアンを構成
A = Placeholder("A")
B = Placeholder("B")

x = Array.create('x',SIZE,'SPIN')

H_A = Constraint(Sum(0,SIZE, lambda i:x[i])**2,label="H_A")

H_B = 0.0
for i in range(SIZE):
    for j in range(SIZE):
        if (G.get_edge_data(i,j) != None):
            H_B = H_B + (1.0 - x[i] * x[j]) / 2.0
            
H = A * H_A + B * H_B


# モデルをコンパイル
model = H.compile()
# QUBOを作成
feed_dict = {'A': 0.2,'B':0.1}
qubo, offset = model.to_qubo(feed_dict=feed_dict)

# PyquboのSAで解いてみる
sol = solve_qubo(qubo)
# 結果表示
pos = nx.spring_layout(G, k=5.,seed=0)
nx.draw_networkx(G,pos,node_color=list(sol.values()))
plt.rcParams['figure.figsize'] = (7.0, 7.0)
plt.show()

次にD-waveで試します。

D-WAVEで実行

from dwave.system.samplers import DWaveSampler
from dwave.system.composites import EmbeddingComposite
import dimod

sampler = EmbeddingComposite(DWaveSampler(solver='DW_2000Q_VFYC_5'))
response = sampler.sample_qubo(qubo, num_reads=100)

result = []
result = model.decode_dimod_response(response,feed_dict=feed_dict)

pos = nx.spring_layout(G, k=5.,seed=0)
nx.draw_networkx(G,pos,node_color=list(result[0][0]['x'].values()))
plt.rcParams['figure.figsize'] = (7.0, 7.0)
plt.show()

お、無事いい感じの答えが出ました！

めでたしめでたし