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

年末年始の時間を利用して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４の場合は、

<img width="268" alt="Screen Shot 2020-01-03 at 17.09.22.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/260261/9b32abd3-12ad-7419-7fa9-a8c916c9619a.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、１の変数を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}

<img width="268" alt="Screen Shot 2020-01-03 at 21.16.39.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/260261/20240fb4-19b2-2e69-a595-170b1fce2d9a.png">

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

番外：
レゴのハリーポッターのアドベントカレンダー2019についていたチェス盤で実験（笑）
![IMG_6770.JPG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/260261/86caad03-9289-319c-75dc-6532f0b2dae8.jpeg)


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

以上。

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