3-PyQUBO with OpenJij

Open in Colab

ここでは,PyQUBOによるコスト関数のQUBOへの変換方法とSimulated Annealing、そしてOpenJijへの変数受け渡しなどについてクリーク被覆問題を例にご紹介いたします。

pipを用いてpyquboをインストールしましょう。

[ ]:
!pip install pyqubo

PyQUBO を使った QUBOの定式化

PyQUBOは、直感的にQUBOやIsing modelを定式化するライブラリです。
これまでの章ではPyQUBOを用いない場合をご紹介してきました。そこではQUBOなどを定式化したのち、自分で式を展開してpythonのコードに落とし込んでいました。しかしPyQUBOを用いることで、その手間をなくすことができます。PyQUBOはQUBO化, Ising model化を施すときの計算ミスや実装ミスを減らすことができる便利なライブラリです。

今回はクリーク被覆問題を例にしてPyQUBOを使ってみます。

クリーク被覆問題に関しては、こちら (T-Wave:クリーク被覆問題) が参考になります。

クリーク被覆問題のQUBO表現による定式化を紹介します。 グラフ \(G=(V, E)\)\(n\)個のクリークで被覆できるかという問題です。 QUBO表現は以下のようになります(ここ:T-Wave と同じ記法を用います)。

\[H = A\sum_v \left(1-\sum^n_{i=1} x_{v, i}\right)^2 + B \sum^n_{i=1}\left[ \frac{1}{2}\left(-1+\sum_{v \in V} x_{v,i}\right)\sum_{v \in V} x_{v, i} - \sum_{(u, v)\in E} x_{u,i} x_{v, i}\right]\]
第一項は各頂点に1色だけ塗るという制約項です。第二項は分割した部分グラフがどれだけクリーク(完全グラフ)に近いかを示しています。どちらの項も0にならなければいけませんが、ここでは第一項をペナルティ項、第二項をコスト(目的関数)として扱うことにします。
このQUBOをPyQUBOを使って表現してみましょう。

今回はグラフとクリークの数\(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')
これで(N_VER)行(N_COLOR)列のバイナリ変数’x’の作成ができました。
次にQUBOを作ります。比較的数式通りに記述することができます。
[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
第一項でPyQUBOのConstraint関数を用いることで、「この項は制約項である」とスクリプトに認識させることができます。
作成したコスト関数は,以下のようにQ.compile().to_qubo() で簡単にQUBO(Pythonの辞書型)に変換することができます。
OpenJijやD-Wave Oceanでは、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))
{}
今回は制約が満たされているため、空の辞書となっています。
decodeの機能はこのように制約が満たされているかを自動的にチェックできる、非常に便利な関数です。

OpenJij に投げる

先ほどはPyQUBOのSAで問題を解いてみました。次はOpenJijに投げてみましょう。
OpenJijでもSAを実行できますが、せっかくなのでPyQUBOには実装されていないSQA(Simulated quantum annealing)を用いたいと思います。
[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との連携を見ました。

手順としては

  1. pyqubo.Arrayで変数を用意

  2. QUBOを定式化

  3. QUBOをコンパイルして辞書型に変換

  4. 辞書型QUBOを受け付けるOpenJijなどのソルバーに投げて最適化問題を解く

  5. 返ってきた解を添字をキーとした辞書型にしてdecodeする

という流れになります。
PyQUBOは定式化、制約の評価を行う非常に便利で協力なツールです。
様々なソルバーを提供するOpenJijとうまく組み合わせて使うことで、より快適な開発を行うことができます。

参考:PyQUBO公式ドキュメント https://pyqubo.readthedocs.io/en/latest/reference/array.html?highlight=arry%20create