年末年始の時間を利用してPyQUBOの理解を深めるためにNクイーン問題をやってみました。
PyQUBOは、アニーリングマシンを使う上で、とても便利なQUBO作成ツールだと思います。
「n個のうちmだけを１にする問題をPyQUBOを使ってアニーリングマシンで解く」
https://qiita.com/morimorijap/items/196e7fc86ecff927bf40
などでも使っていますが、改めて、Nクイーン問題をPyQUBOを使って解いてみようと思います。
Nクイーン問題とは、Wikipediaの８クイーンは、マスが８x８個で８クイーン https://ja.wikipedia.org/wiki/エイト・クイーン
に対して、マスがNxN個の問題を言います。
クイーンの動きは、上下左右斜めの8方向に、遮る物がない限り進めるので、マスの列、行、斜めにクイーンがいないように配置する問題となります。
これをアニーリング マシンで解くためにコスト関数の形式にすると、
N-Queen問題のハミルトニアン、 列column、行row、斜めdiagonal
$$H = \sum_{col}\left(\sum_{i \in col} x_i - 1\right)^2
+ \sum_{row}\left(\sum_{i \in row} x_i - 1\right)^2
+ \sum_{diag} \left((\sum_{i \in diag}
x_i ) (\sum_{i \in diag} x_i - 1 )\right)$$
となります。
斜めのコスト関数の表現が異なりますが、東北大学のT-Waveの記事になる「N-クイーン問題 をD-Waveマシンで解く」https://qard.is.tohoku.ac.jp/T-Wave/?p=884
も参考になります。
４x４の場合は、
のようなマスの表現を考え、PyQUBOで0,1のQUBO形式 Binary(q[0]) のようなArrayを表現します。
from pyqubo import Array, Placeholder, solve_qubo, Sum
from pprint import pprint
import numpy as np
まずは、PyQUBO など必要なものをインポートします。
4x4のマスを作るとします。
# N-Queenのマスの数
N = 4
#pyQUBOで0、１の変数をNxN正方マス分つくる
q = Array.create('q', (N*N), 'BINARY')
q_shape = np.reshape(q,(N,N))
print(q_shape)
そうすると、
[[Binary(q[0]) Binary(q[1]) Binary(q[2]) Binary(q[3])]
[Binary(q[4]) Binary(q[5]) Binary(q[6]) Binary(q[7])]
[Binary(q[8]) Binary(q[9]) Binary(q[10]) Binary(q[11])]
[Binary(q[12]) Binary(q[13]) Binary(q[14]) Binary(q[15])]]
このようにPyQUBOのBinaryを４x４配列の形式に作成することができます。
行の部分は、数式で以下のようになっていますので、
$$\sum_{row}\left(\sum_{i \in row} x_i - 1\right)^2 $$
これをPyQUBOで表すと、行ごとスライスしてそのまま足し２乗し、最後に合計するように表せます。
#行
row_const = 0.0
for row in range(N):
row_const += (sum(i for i in q_shape[row, :]) -1)**2
また、列もどうように以下のように表せます。
$$ \sum_{col}\left(\sum_{i \in col} x_i - 1\right)^2 $$
#列
col_const = 0.0
for col in range(N):
col_const += (sum(i for i in q_shape[:, col]) -1)**2
となります。
そして、問題の斜めの表現ですが、
$$\sum_{diag} \left((\sum_{i \in diag}
x_i ) (\sum_{i \in diag} x_i - 1 )\right)$$
をPyQUBOで表すことになりますが、Numpyの変換を利用し、「＼」「/」の２つに分けて記述することができます。（もっといい方法をご存知の方は教えてください）
#斜め
# ＼の方
diag_const = 0.0
xi = 0.0
xi_1 = 0.0
for k in range(-N+1,N):
xi += (sum(i for i in np.diag(q_shape,k=k)))
xi_1 += (sum(i for i in np.diag(q_shape,k=k))) -1
diag_const += xi * xi_1
と
#斜め
# /の方のために入れ替えて, ＼ をする。
diag_const_f = 0.0
xi = 0.0
xi_1 = 0.0
for k in range(-N+1,N):
xi += (sum(i for i in np.diag(np.fliplr(q_shape),k=k)))
xi_1 += (sum(i for i in np.diag(np.fliplr(q_shape),k=k))) -1
diag_const_f += xi * xi_1
となります。
これを全て足すとことで、エネルギー関数のハミルトニアンを表現することができます。
また、パラメータ調整のためにPyQUBOのPlaceholderを使いalpha,beta,gammaをつけておきます。
# エネルギー (ハミルトニアン) を構築
alpha = Placeholder("alpha")
beta = Placeholder("beta")
gamma = Placeholder("gamma")
#パラメータ
feed_dict = {'alpha': 1.0, 'beta': 1.0, 'gamma': 1.0}
H = alpha * col_const + beta * row_const + gamma *(diag_const + diag_const_f)
次にQUBOへ変換をします。
#QUBOのコンパイル
model = H.compile()
qubo, offset = model.to_qubo(feed_dict=feed_dict)
pprint(qubo)
だけで、QUBOを作ってくれます。
{('q[0]', 'q[0]'): -11.0,
('q[0]', 'q[10]'): 6.0,
('q[0]', 'q[11]'): 4.0,
('q[0]', 'q[12]'): 2.0,
('q[0]', 'q[13]'): 6.0,
('q[0]', 'q[14]'): 6.0,
('q[0]', 'q[15]'): 6.0,
('q[0]', 'q[1]'): 6.0,
('q[0]', 'q[2]'): 4.0,
・・・長いので、途中省略・・・
('q[8]', 'q[8]'): -19.0,
('q[8]', 'q[9]'): 14.0,
('q[9]', 'q[9]'): -21.0}
です。 これをそのまま、PyQUBOのSAで計算することができます。
# PyQUBOのSAで計算
solution = solve_qubo(qubo)
# pyQUBOのSAで計算結果のデコード
decoded_solution, broken, energy = model.decode_solution(solution, vartype="BINARY",
feed_dict=feed_dict)
pprint(solution)
結果は、
{'q[0]': 0,
'q[10]': 0,
'q[11]': 0,
'q[12]': 0,
'q[13]': 0,
'q[14]': 1,
'q[15]': 0,
'q[1]': 1,
'q[2]': 0,
'q[3]': 0,
'q[4]': 0,
'q[5]': 0,
'q[6]': 0,
'q[7]': 1,
'q[8]': 1,
'q[9]': 0}
となり、一応条件は満たしていますが、13にQueenが入ってもいいので、パラメータの調整などが必要な感じがします。
追記：斜めの一番角が抜けていましたので、最適解が出てなかったようです。
今回はちゃんと解が出るようになりました。
番外：
レゴのハリーポッターのアドベントカレンダー2019についていたチェス盤で実験（笑）
PyQUBOについては、公式ドキュメント参考になります。
https://pyqubo.readthedocs.io/
以上。
追記：
2020/1/3 9:20pm
斜めの条件が間違っていたので編集しました。