Fixstars Amplifyで，制約充足問題(N-Queen問題)に挑戦してみた

はじめに

こんにちは，@khrive　です．
再帰を使ったバックトラック等でパズル問題を解く方法は従来から数多く例示されています．
近年では，量子アニーリングを適用して身近な最適化問題，制約充足問題を求解するノウハウが広がっていると思います．FixstarsのAmplifyは，とてもわかり易く使いやすいプラットフォームと開発フレームワークで大変ありがたく使わせもらっています．今回,Amplifyを使って制約充足問題を解くこと，具体的にはN-Queen問題の求解に挑戦します．

N-Queen問題の求解例

• 動作環境

  1. MacOS sonoma 14.1.2, Python 3.8.13, pip 23.3.1, VSCode,JupyterLab
  2. WindowsではWSL2　UBUNTUなどLinux環境が必要
  3. Google colaboraotory (https://colab.research.google.com/)

• 事前準備
  Fixstars Amplifyにサインアップ，Amplify Access Tokenを取得する：
  URL: (https://amplify.fixstars.com/ja/) の「無料でアクセストークンを入手」をクリックして
  進める

!python3 -m pip install update pip
!python3 -m pip install --upgrade pip
!python3 -m pip install amplify

• 作成したコード

import numpy as np
from amplify import BinaryPoly,sum_poly,Solver,SymbolGenerator
from amplify.constraint import equal_to,clamp,less_equal,one_hot
from amplify.client import FixstarsClient

# Amplify Access Tokenを取得して上書き　（https://amplify.fixstars.com/ja/user/token）
MyAmplifyToken = 'xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx'

#  ソルバの取得
def get_solver():
    client = FixstarsClient()
    client.token = MyAmplifyToken
    client.parameters.timeout = 1000
    return Solver(client)

#  結果のデコード
def print_Q(mat):
    def elm(m):
        return 'Q ' if m>0 else '_ '
    for i in range(len(mat)):
        line = ''
        for j in range(len(mat[i])):
            line += elm(mat[i][j])
        print(line)
    print()    

#  ななめ方向の制約条件
def enum_diags(N):
    diags = []
    for i in range(N):
        diag = [(i-j,j) for j in range(N) if i-j>=0]
        diags += [diag]
    for i in range(N):
        diag = [(i+j,j) for j in range(N) if i+j<N]
        diags += [diag] 
    return diags

#  求解
def solve(N=8):
    rN = range(N)
    gen = SymbolGenerator(BinaryPoly)
    q = gen.array(N,N)
    diags = enum_diags(N)

    # 行に対する制約
    # row_const = [one_hot(q[n]) for n in rN]
    row_const = [equal_to(sum_poly([q[i,j] for i in rN]),1) for j in rN]
 
    # 列に対する制約
    # col_const = [one_hot(q[:, i]) for i in rN]
    col_const = [equal_to(sum_poly([q[i,j] for j in rN]),1) for i in rN]
    
    diag_const =  [ clamp(sum_poly([q[i,j] for i,j in d]),0,1) for d in diags ]
    # diag_const =  [ less_equal(sum_poly([q[i,j] for i,j in d]),1) for d in diags ]   
    
    # 行，列，斜めの制約条件を合成，問題設定
    constraints = sum(row_const) + sum(col_const) + sum(diag_const)
    
    # Amplify　ソルバ取得，求解
    solver = get_solver()
    results = solver.solve(constraints)

    # 求解結果のデコード
    print(results)
    if results[0]:
        vs = q.decode(results[0].values)
        print_Q(vs)

# 求解実行
if __name__ == '__main__':
    solve(N=8)
    solve(N=16)
    solve(N=25)

• 実行結果　（実行の都度，変わります）

<amplify.amplify.SolverResult object at 0x7fa4188c6e70>
_ Q _ _ _ _ _ _ 
_ _ _ _ Q _ _ _ 
_ _ _ _ _ Q _ _ 
_ _ _ _ _ _ _ Q 
Q _ _ _ _ _ _ _ 
_ _ _ Q _ _ _ _ 
_ _ _ _ _ _ Q _ 
_ _ Q _ _ _ _ _ 

<amplify.amplify.SolverResult object at 0x7fa3e9628730>
Q _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ Q _ _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ Q _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ Q _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ Q _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ Q _ _ _ 
_ Q _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ _ Q _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ Q 
_ _ _ _ _ _ _ Q _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ Q _ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ Q _ 
_ _ _ _ _ _ _ _ Q _ _ _ _ _ _ _ 
_ _ _ _ _ _ Q _ _ _ _ _ _ _ _ _ 
_ _ Q _ _ _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ Q _ _ _ _ _ _ _ _ _ _ 

<amplify.amplify.SolverResult object at 0x7fa418a7e570>
_ _ _ _ _ _ _ _ _ _ _ _ _ Q _ _ _ _ _ _ _ _ _ _ _ 
_ Q _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ Q _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ _ Q _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ Q _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ Q _ _ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ Q _ _ _ 
_ _ _ _ _ _ _ _ Q _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ Q _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ Q _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ Q _ 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ Q _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ Q 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ Q _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ Q _ _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ Q _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ Q _ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ Q _ _ _ _ _ 
Q _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ Q _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ Q _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ Q _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ Q _ _ _ _ _ _ _ _ _ _ _ _ _ 
_ _ Q _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ Q _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

制約条件違反に対するペナルティ項の重み係数＝ハイパーパラメタ調整

上記の求解結果をざっと検証すると，行と列方向の制約条件は充足する一方で，斜め方向制約を充足しない傾向が散見されます．チェスボードと問題の左右対称性，点対称性を考えると，行／列／斜め方向の制約違反のペナルティ重みは同程度と考えられますが，斜め方向の制約違反を抑制する方法として，

    constraints = sum(row_const) + sum(col_const) + sum(diag_const)

　上の右辺３項目のsum(diag_const)に適切な重み$w$　を乗せて試す価値はありそうです．

    w = 3
    constraints = sum(row_const) + sum(col_const) + w * sum(diag_const)

として，改めて求解した結果が以下です：

<amplify.amplify.SolverResult object at 0x7fa419419bf0>
_ _ Q _ _ _ _ _ 
Q _ _ _ _ _ _ _ 
_ _ _ _ _ _ Q _ 
_ _ _ _ Q _ _ _ 
_ _ _ _ _ _ _ Q 
_ _ _ _ _ Q _ _ 
_ _ _ Q _ _ _ _ 
_ Q _ _ _ _ _ _ 

<amplify.amplify.SolverResult object at 0x7fa3e962aa30>
Q _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ Q _ _ _ _ _ _ _ _ _ 
_ _ _ Q _ _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ _ Q _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ Q _ _ _ _ 
_ Q _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ _ Q _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ Q _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ Q _ 
_ _ _ _ _ _ _ _ _ _ Q _ _ _ _ _ 
_ _ _ _ _ _ _ _ Q _ _ _ _ _ _ _ 
_ _ _ _ _ Q _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ Q _ _ _ _ _ _ 
_ _ _ _ Q _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ Q 
_ _ Q _ _ _ _ _ _ _ _ _ _ _ _ _ 

<amplify.amplify.SolverResult object at 0x7fa418a5bcb0>
_ _ _ _ _ _ _ _ Q _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ Q _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
_ _ Q _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
Q _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ Q 
_ _ _ _ _ _ _ _ _ _ _ Q _ _ _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ Q _ 
_ _ _ Q _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ Q _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ Q _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ Q _ _ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ Q _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ Q _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ Q _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ Q _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ Q _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ Q _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ Q _ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ Q _ _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ Q _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
_ Q _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ Q _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ _ Q _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ _ Q _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ Q _ _ _ _ _ _ _ _ 

斜め方向の制約違反は若干改善できたでしょうか．
Amplifyドキュメント(https://amplify.fixstars.com/ja/docs/solver.html) によれば実行可能解か？どの制約条件が満たされなかったか？を調べるis_feasible, check_constraints()を使う手法が書かれています．

QUBO変数を使って条件式を表す定石例　

　量子アニーリングにおける量子ビットは，古典ビットと似た感覚で{0,1}で表現しますが，制約条件には下記のように頻出するパターンが存在します．

条件式頻出パタン

• 同時にはどれか一つの変数だけ、必ず1である（他は全て0）
  Amplifyではone_hot() 関数を利用可能

\sum_{i=0}^{n-1}x_i = 1  \iff (\sum_{i=0}^{n-1}x_i -1)^2 = 0

• どれか1つの変数だけ同時に1であってもよい（１ or 0でも良い）
  Amplifyではclamp(), less_equall(),sum_ploy() 関数を利用可能

\sum_{i=0}^{n-1}x_i = 1  \space  or \space  \sum_{i=0}^{n-1}x_i = 0

\iff

(\sum_{i=0}^{n-1}x_i-1)*(\sum_{i=0}^{n-1}x_i) = 0

おわりに

　FixstarsのAmplifyを活用して，N-Queen問題（制約充足問題）の求解に挑戦しました．とても使いやすいプラットフォームとフレームワークであり，習得すれば何かと役に立つと思います．量子コンピューティング応用や実用化までにまだ時間的猶予があるので，社会課題や身近な最適化問題に適用してノウハウを磨き，スキル人材を育成には良いタイミングだと思います．