3-PyQUBO with OpenJij

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

pip で openjij と pyqubo をインストールしましょう。

[1]:
#!pip install openjij
# !pip install pyqubo
[2]:
!pip show openjij  # openjij == 0.0.5 を仮定する
Name: openjij
Version: 0.0.5
Summary: Framework for ising model and QUBO
Home-page: https://openjij.github.io/OpenJij/
Author: Jij Inc.
Author-email: openjij@j-ij.com
License: Apache License 2.0
Location: /usr/local/miniconda3/lib/python3.6/site-packages
Requires: requests, numpy
Required-by:
[3]:
!pip show pyqubo  # pyqubo == 0.3.0 を仮定する
Name: pyqubo
Version: 0.3.0
Summary: PyQUBO allows you to create QUBOs or Ising models from mathematical expressions.
Home-page: https://github.com/recruit-communications/pyqubo
Author: Recruit Communications Co., Ltd.
Author-email: rco_pyqubo@ml.cocorou.jp
License: Apache 2.0
Location: /usr/local/miniconda3/lib/python3.6/site-packages
Requires: six, numpy, dimod, dwave-neal
Required-by:

PyQUBO を使った QUBOの定式化

PyQUBOは直感的にQUBOやIsingモデルを定式化するライブラリです。
PyQUBOがない場合は QUBOなどを定式化したあと自分で式を展開してpythonのコードに落とし込んでいましたが、PyQUBOのおかげでその手間をなくすことができます。 PyQUBOで計算ミスや、実装ミスを減らすことができる便利なライブラリです。

今回はクリーク被覆問題という問題を例にして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\)を以下のように与えます。

[4]:
#頂点の数
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 しておきます。

[5]:
from pyqubo import Array, Constraint, solve_qubo

まずQUBOを表現する変数を用意しましょう。Arrayを使って変数列を作ります。 今回は 頂点の数 × 色の数が必要なので その分の shape で用意しましょう。

[6]:
x = Array.create('x', shape=(N_VER,N_COLOR), vartype='BINARY')

次にQUBOを作ります。比較的数式通りに記述することができます。

[7]:
# 第一項
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

作成したコスト関数は,以下のようにQ.compile().to_qubo() で簡単にQUBO(Pythonの辞書型)に変換することができます。OpenJijやD-Wave Ocean ではQUBOはPythonの辞書型で表現されることが前提になっているのでcompile すればあとはそのまま各ソルバーに投げることができます。

[8]:
# モデルをコンパイル
model = Q.compile()
qubo, offset = model.to_qubo()

PyQUBOにもSimulated Annealingのソルバーが付いています。D-Wave Ocean SDKの中のnealというライブラリを中で呼んでいます。 solve_qubo(qubo)を用いることでSAを実行できます。

また出てきた解をデコードして読みやすい形に整えてくれるメソッドがあるのでそれを使ってデコードします。

[9]:
#PyQUBOによるSA
raw_solution = solve_qubo(qubo)
# 得られた結果をデコードする
decoded_solution, broken, energy = model.decode_solution(raw_solution, vartype="BINARY")
[10]:
decoded_solution['x']
[10]:
{0: {1: 1, 2: 0, 3: 0},
 1: {1: 1, 2: 0, 3: 0},
 2: {1: 1, 2: 0, 3: 0},
 3: {1: 0, 2: 0, 3: 1},
 4: {1: 0, 2: 0, 3: 1},
 5: {1: 0, 2: 1, 3: 0},
 6: {1: 0, 2: 1, 3: 0},
 7: {1: 0, 2: 1, 3: 0}}

(0,1,2) と (3, 4), (5,6,7) のグループに分けられたようです。これは今回与えたグラフ上でちゃんとそれぞれクリークになっています。

また broken はペナルティ項(今回は第一項)が破れていると(0じゃないと)どのように破れたかを記録してくれます。
decodeの機能はこのように制約が満たされているかを自動的にチェックでき、非常に便利な関数です。

OpenJij に投げる

ではOpenJijに投げてみましょう。
OpenJijでもSAを実行できますが、せっかくなのでPyQUBOには実装されていないSQA(Simulated quantum annealing)を使って解いてみます。
[11]:
import openjij as oj

sampler = oj.SQASampler()
response = sampler.sample_qubo(Q=qubo)

sampler の部分を取り替えれば他のアルゴリズムやマシンに投げることができるので試してみてください。

ではOpenJijで返ってきた結果をPyQUBOのデコーダーを使ってデコードしてみましょう。 PyQUBOのデコーダーは引数としてsolverの結果を、コンパイルして作った辞書型QUBOの持つindexをキーとし、値を2値変数の値とした辞書型であると期待しています。

なので辞書型で渡してやれば良いです。具体的には以下のようにします。

[12]:
dict_solution = response.samples[0]
decoded_solution, broken, energy = model.decode_solution(dict_solution, vartype="BINARY")

# openjijの内部の処理で順番が崩れてしまうので sort して見やすくしておく
x_solution = dict(sorted(decoded_solution['x'].items()))
{key:dict(sorted(value.items())) for key, value in x_solution.items()}
[12]:
{0: {1: 0, 2: 0, 3: 1},
 1: {1: 0, 2: 0, 3: 1},
 2: {1: 0, 2: 0, 3: 1},
 3: {1: 1, 2: 0, 3: 0},
 4: {1: 1, 2: 0, 3: 0},
 5: {1: 0, 2: 1, 3: 0},
 6: {1: 0, 2: 1, 3: 0},
 7: {1: 0, 2: 1, 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