4
Help us understand the problem. What are the problem?

More than 1 year has passed since last update.

posted at

updated at

Nクイーン問題をPyQUBOのSAで解いてみる

年末年始の時間を利用してPyQUBOの理解を深めるためにNクイーン問題をやってみました。

PyQUBOは、アニーリングマシンを使う上で、とても便利なQUBO作成ツールだと思います。
「n個のうちmだけを1にする問題をPyQUBOを使ってアニーリングマシンで解く」
https://qiita.com/morimorijap/items/196e7fc86ecff927bf40
などでも使っていますが、改めて、Nクイーン問題をPyQUBOを使って解いてみようと思います。

Nクイーン問題とは、Wikipediaの8クイーンは、マスが8x8個で8クイーン 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
も参考になります。

4x4の場合は、

Screen Shot 2020-01-03 at 17.09.22.png

のようなマスの表現を考え、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、1の変数を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を4x4配列の形式に作成することができます。
行の部分は、数式で以下のようになっていますので、

$$\sum_{row}\left(\sum_{i \in row} x_i - 1\right)^2 $$

これをPyQUBOで表すと、行ごとスライスしてそのまま足し2乗し、最後に合計するように表せます。

#行
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の変換を利用し、「\」「/」の2つに分けて記述することができます。(もっといい方法をご存知の方は教えてください)

#斜め 
# \の方
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}

Screen Shot 2020-01-03 at 21.16.39.png

となり、一応条件は満たしていますが、13にQueenが入ってもいいので、パラメータの調整などが必要な感じがします。
追記:斜めの一番角が抜けていましたので、最適解が出てなかったようです。
今回はちゃんと解が出るようになりました。

番外:
レゴのハリーポッターのアドベントカレンダー2019についていたチェス盤で実験(笑)
IMG_6770.JPG

PyQUBOについては、公式ドキュメント参考になります。
https://pyqubo.readthedocs.io/

以上。

追記:
2020/1/3 9:20pm
斜めの条件が間違っていたので編集しました。

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Sign upLogin
4
Help us understand the problem. What are the problem?