3-PyQUBO with OpenJij¶
ここでは,PyQUBOによるコスト関数のQUBOへの変換方法とSimulated Annealing、そしてOpenJijへの変数受け渡しなどについてクリーク被覆問題を例にご紹介いたします。
pipを用いてpyqubo
をインストールしましょう。
[ ]:
!pip install pyqubo
PyQUBO を使った QUBOの定式化¶
PyQUBO
は、直感的にQUBOやIsing modelを定式化するライブラリです。今回はクリーク被覆問題を例にしてPyQUBOを使ってみます。
クリーク被覆問題に関しては、こちら (T-Wave:クリーク被覆問題) が参考になります。
クリーク被覆問題のQUBO表現による定式化を紹介します。 グラフ \(G=(V, E)\)を\(n\)個のクリークで被覆できるかという問題です。 QUBO表現は以下のようになります(ここ:T-Wave と同じ記法を用います)。
今回はグラフとクリークの数\(n\)を以下のように与えます。
[1]:
# 頂点の数を定義します。
N_VER = 8
# 色の数を定義します。
N_COLOR = 4
# グラフを定義します。0~7の頂点があったときにどの頂点同士が線で結ばれているかを定義しています。
graph = [(0,1), (0,2), (1,2), (5,6), (2,3), (2,5), (3,4), (5,7), (7, 6)]
PyQUBOによる定式化¶
PyQUBOから必要となるクラスをimport しておきます。
[2]:
from pyqubo import Array, Constraint, solve_qubo
最初にQUBOを表現する変数を用意しましょう。Arrayを使って変数列を作ります。 今回は(頂点の数)×(色の数)だけ変数が必要です。その分をshape
を用いて準備します。
[3]:
x = Array.create('x', shape=(N_VER,N_COLOR), vartype='BINARY')
[4]:
# 第一項 (制約項)を定義します。
H_A = Constraint(sum((1-sum(x[v,i] for i in range(1,N_COLOR)))**2 for v in range(N_VER)), label='HA')
# 第二項 (コスト、目的関数)を定義します。
H_B = sum((-1+sum(x[v,i] for v in range (N_VER)))/2*(sum(x[v,i] for v in range (N_VER))) - sum(x[u,i]*x[v,i] for (u,v) in graph) for i in range (1,N_COLOR))
# ハミルトニアン全体を定義します。
Q = H_A+H_B
Constraint
関数を用いることで、「この項は制約項である」とスクリプトに認識させることができます。Q.compile().to_qubo()
で簡単にQUBO(Pythonの辞書型)に変換することができます。.compile
すればあとはそのまま各ソルバーに投げることが可能です。[5]:
# モデルをコンパイルします。
model = Q.compile()
qubo, offset = model.to_qubo()
qubo
にはPythonの辞書型で格納されたQUBOが、そしてoffset
にはQUBO化した際に現れる定数(無視してよい)が代入されます。
PyQUBOに備わっていたSimulated Annealingソルバーのsolve_qubo
は使用が非推奨となりました。ここではD-Wave Ocean SDKのnealというライブラリを直接呼び出して実行しましょう。
[6]:
# nealのSAを用いて解きます。
import neal
sampler = neal.SimulatedAnnealingSampler()
raw_solution = sampler.sample_qubo(qubo)
print(raw_solution)
x[0][1] x[0][2] x[0][3] x[1][1] x[1][2] x[1][3] ... x[7][3] energy num_oc.
0 0 1 0 0 1 0 ... 0 -8.0 1
['BINARY', 1 rows, 1 samples, 24 variables]
得られた解の中で一番低いエネルギーの状態を抽出します。
[7]:
raw_solution.first.sample
[7]:
{'x[0][1]': 0,
'x[0][2]': 1,
'x[0][3]': 0,
'x[1][1]': 0,
'x[1][2]': 1,
'x[1][3]': 0,
'x[2][1]': 0,
'x[2][2]': 1,
'x[2][3]': 0,
'x[3][1]': 0,
'x[3][2]': 0,
'x[3][3]': 1,
'x[4][1]': 0,
'x[4][2]': 0,
'x[4][3]': 1,
'x[5][1]': 1,
'x[5][2]': 0,
'x[5][3]': 0,
'x[6][1]': 1,
'x[6][2]': 0,
'x[6][3]': 0,
'x[7][1]': 1,
'x[7][2]': 0,
'x[7][3]': 0}
得られた解を見てみると、‘x[0][0]’: 1のように文字列をキーにした辞書型で格納されていることがわかります。このままだと、今後の解析がしづらいです。PyQUBOにはそれを扱いやすい形に直すデコード機能.decode_sample()
があります。
[8]:
# 得られた結果をデコードします。
decoded_sample = model.decode_sample(raw_solution.first.sample, vartype="BINARY")
# さらに解を見やすくする処理を追加します。
# .array(変数名, 要素番号)で希望する要素の値を抽出することができます。
x_solution = {}
for i in range(N_VER):
x_solution[i] = {}
for j in range(1,N_COLOR):
x_solution[i][j] = decoded_sample.array('x', (i, j))
x_solution
[8]:
{0: {1: 0, 2: 1, 3: 0},
1: {1: 0, 2: 1, 3: 0},
2: {1: 0, 2: 1, 3: 0},
3: {1: 0, 2: 0, 3: 1},
4: {1: 0, 2: 0, 3: 1},
5: {1: 1, 2: 0, 3: 0},
6: {1: 1, 2: 0, 3: 0},
7: {1: 1, 2: 0, 3: 0}}
(0,1,2), (3, 4), (5,6,7)にそれぞれ分けられたようです。これは今回与えたグラフ上でそれぞれクリークになっています。
得られた解が制約を満たしているか確認する方法をご紹介します。それには.constraints(only_broken=True)
で結果を表示しましょう。これは制約項(今回は第一項)が破れているとき(0でないとき)、どのように破れたかを記録してくれます。
[9]:
print(decoded_sample.constraints(only_broken=True))
{}
OpenJij に投げる¶
[10]:
# OpenJijをインポートします。
import openjij as oj
# SQAを用いて問題を解きます。
sampler = oj.SQASampler()
# QUBOに先ほど.to_quboで作成したものを代入します。
response = sampler.sample_qubo(Q=qubo)
print(response)
x[0][1] x[0][2] x[0][3] x[1][1] x[1][2] x[1][3] ... x[7][3] energy num_oc.
0 0 1 0 0 1 0 ... 0 -8.0 1
['BINARY', 1 rows, 1 samples, 24 variables]
sampler
の部分を取り替えるだけで、他のアルゴリズムやマシンに投げることができます。興味のある方は試してみると良いでしょう。
ではOpenJijで返ってきた結果を、PyQUBOのデコーダーを使ってデコードしてみましょう。 具体的には以下のようにします。
[11]:
# エネルギーが一番低い状態を取り出します。
dict_solution = response.first.sample
# 得られた結果をデコードします。
decoded_sample = model.decode_sample(raw_solution.first.sample, vartype="BINARY")
# さらに解を見やすくする処理を追加します。
# .array(変数名, 要素番号)で希望する要素の値を抽出することができます。
x_solution = {}
for i in range(N_VER):
x_solution[i] = {}
for j in range(1,N_COLOR):
x_solution[i][j] = decoded_sample.array('x', (i, j))
x_solution
[11]:
{0: {1: 0, 2: 1, 3: 0},
1: {1: 0, 2: 1, 3: 0},
2: {1: 0, 2: 1, 3: 0},
3: {1: 0, 2: 0, 3: 1},
4: {1: 0, 2: 0, 3: 1},
5: {1: 1, 2: 0, 3: 0},
6: {1: 1, 2: 0, 3: 0},
7: {1: 1, 2: 0, 3: 0}}
まとめ¶
PyQUBOを使って定式化する方法とOpenJijとの連携を見ました。
手順としては
pyqubo.Arrayで変数を用意
QUBOを定式化
QUBOをコンパイルして辞書型に変換
辞書型QUBOを受け付けるOpenJijなどのソルバーに投げて最適化問題を解く
返ってきた解を添字をキーとした辞書型にしてdecodeする
参考:PyQUBO公式ドキュメント https://pyqubo.readthedocs.io/en/latest/reference/array.html?highlight=arry%20create