N-QueensパズルのQUBOを高速生成し，Gurobiで解く

# ポイント
* N-QueensパズルのQUBOは，one-hot制約の数式を展開するといった方法を用いず，直接的に作ることができる
* Gurobiのコールバック関数を設定し，指定した値の解が得られたら解探索を終了する方法

# N-QueensパズルのQUBOを高速生成
[N-QueensパズルのQUBOをSymPyで生成し，Gurobiで解く](https://qiita.com/koji_nakano/items/860bd9fb46072af986af)では，one-hot制約およびzero-one-hot制約を満たす式を構築し，SymPyを用いて展開してQUBOを生成しました．数式処理を伴うため，ボードのサイズが大きくなると計算にかなりの時間がかかります．ここでは、より直観的かつ簡単にN-QueensパズルのQUBOを生成する方法を用います．
以前と同様に，バイナリ変数の$n\times n$行列$X=(x_{i,j})$を使用してボードを表現します．ここで，$x_{i,j}=1$は座標$(i,j)$にQueenが配置されていることを意味します．このバイナリ変数の総数を最大化しつつ，クイーンの位置が衝突する場合にペナルティを課すようなQUBOを次のように構築します．
* 全ての1次項$x_{i,j}$の係数を$-1$と設定することで，$x_{i,j}=1$となるバイナリ変数の個数を最大化します
* 座標$(i,j)$と$(i',j')$が同一行、同一列、または同一斜めライン上にある場合、2次項$x_{i,j}x_{i',j'}$の係数を$+1$と設定し，Queenが衝突することをペナルティとします 

これにより，$X$がN-Queensパズルの正しい解である場合に限り，最適解は$-n$となります．$n+1$個の$x_{i,j}=1$が存在する場合，1次項からの寄与は$-n-1$となりますが、少なくとも2箇所での衝突が発生するため、2次項からの寄与は$+2$以上になります．したがって，そのような解は$-n+1$以上の値を持ち，最適解とはなり得ません．$n+1$個を超える$x_{i,j}=1$が存在する場合，解はさらに悪化します．よって，$X$がN-Queensパズルの正しい解である場合に限り，かつその場合にのみ，最適解が$-n$になります．

こうして得られたQUBOは、実際には[N-QueensパズルのQUBOをSymPyで生成し，Gurobiで解く](https://qiita.com/koji_nakano/items/860bd9fb46072af986af)で導出されたものと本質的に同じです．具体的には，本稿の方法で得られるQUBOの全係数を2倍にしたものが，後者のQUBOに完全に一致します．

# Pythonプログラム
* `find_conflict()`で，`self.conflict`に衝突する座標のペア（非順序対）を求めます.
* `var[i,j]`はGurobiの変数で$x_{i,j}$に相当します.
* `obj`には，`self.conflict`を参照しつつ，N-Queensパズルの条件を満たすための`var[i,j]`の項を追加しています.
* `mycallback`はGurobiによって呼び出されるコールバック関数です．この関数は、最適解$-n$を得た時点でGurobiの実行を直ちに停止するように設定されています．

```Python
#!/usr/bin/env python3

import argparse
import gurobipy as gp


class NQUEEN:
    def __init__(self, n):
        self.n = n

    def find_conflict(self):
        self.conflict = set()
        for x in range(self.n):
            for y in range(self.n):
                for i in range(x + 1, self.n):
                    self.conflict.add(frozenset([(x, y), (i, y)]))
                for j in range(y + 1, self.n):
                    self.conflict.add(frozenset([(x, y), (x, j)]))
                for d in range(1, self.n - x):
                    if y - d >= 0:
                        self.conflict.add(frozenset([(x, y), (x + d, y - d)]))
                    if y + d < self.n:
                        self.conflict.add(frozenset([(x, y), (x + d, y + d)]))

    def mycallback(self, model, where):
        if where == gp.GRB.Callback.MIPSOL:
            if int(model.cbGet(gp.GRB.Callback.MIPSOL_OBJ)) == -self.n:
                model.terminate()

    def gurobi(self, timelimit):
        self.model = gp.Model()
        self.model.setParam("TimeLimit", timelimit)
        var = {}
        for i in range(self.n):
            for j in range(self.n):
                var[i, j] = self.model.addVar(vtype=gp.GRB.BINARY, name=f"x_{i}_{j}")
        obj = [-var[i, j] * var[i, j] for i in range(self.n) for j in range(self.n)]
        for u, v in self.conflict:
            obj.append(var[u] * var[v])
        self.model.setObjective(gp.quicksum(obj), sense=gp.GRB.MINIMIZE)

        def wrapper(model, where):
            self.mycallback(model, where)

        self.model.optimize(wrapper)

    def print(self):
        for i in range(self.n):
            for j in range(self.n):
                print(int(self.model.getVarByName(f"x_{i}_{j}").X), end="")
            print()
        print(f"optimal={-self.n} obtained={int(self.model.ObjVal)}")


def main():
    parser = argparse.ArgumentParser(description="Solve N-Queens puzzle using Gurobi")
    parser.add_argument(
        "-n", default=32, type=int, help="Value of N of N-QUEENS (default:32)"
    )
    parser.add_argument(
        "-t", default=100, type=int, help="Time limit of Gurobi optimizer (default:100)"
    )
    args = parser.parse_args()

    nqueen = NQUEEN(args.n)
    nqueen.find_conflict()
    nqueen.gurobi(args.t)
    nqueen.print()


if __name__ == "__main__":
    main()
```

# 実行例
$n=32$の場合，最適解として$-32$が得られ，即座に処理が完了しています．SymPyを利用した場合、$n$が増加するにつれて（例えば$n=100$など），QUBOの生成にはかなりの時間がかかります．一方、ここで示した手法では、QUBOの生成が瞬時に行えます．

```
Gurobi Optimizer version 10.0.0 build v10.0.0rc2 (linux64)

CPU model: Intel(R) Xeon(R) Gold 6338 CPU @ 2.00GHz, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 64 physical cores, 64 logical processors, using up to 32 threads

Optimize a model with 0 rows, 1024 columns and 0 nonzeros
Model fingerprint: 0xf7c7dec7
Model has 53600 quadratic objective terms
Variable types: 0 continuous, 1024 integer (1024 binary)
Coefficient statistics:
  Matrix range     [0e+00, 0e+00]
  Objective range  [0e+00, 0e+00]
  QObjective range [2e+00, 2e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [0e+00, 0e+00]
Found heuristic solution: objective 0.0000000
Found heuristic solution: objective -27.0000000
Found heuristic solution: objective -32.0000000
Presolve time: 0.00s

Explored 0 nodes (0 simplex iterations) in 0.28 seconds (0.38 work units)
Thread count was 1 (of 64 available processors)

Solution count 3: -32 -27 0
No other solutions better than 0

Solve interrupted
Best objective -3.200000000000e+01, best bound -, gap -

User-callback calls 212, time in user-callback 0.00 sec
10000000000000000000000000000000
00000000000000000100000000000000
00000000000000000000001000000000
00001000000000000000000000000000
00000000000000000000000100000000
00000000000000000000000000000001
00000000000000000000100000000000
00000000000000000000000000000010
00000000000000000000010000000000
00000000000010000000000000000000
00100000000000000000000000000000
00000010000000000000000000000000
00010000000000000000000000000000
01000000000000000000000000000000
00000000100000000000000000000000
00000000000000000000000001000000
00000000000000000000000000010000
00000000000001000000000000000000
00000000000000000000000000100000
00000000000000001000000000000000
00000100000000000000000000000000
00000000000100000000000000000000
00000000000000000000000000001000
00000000001000000000000000000000
00000001000000000000000000000000
00000000000000000000000010000000
00000000000000000001000000000000
00000000000000000000000000000100
00000000000000100000000000000000
00000000000000000010000000000000
00000000010000000000000000000000
00000000000000010000000000000000
optimal=-32 obtained=-32
```

# 参考文献
* Shunsuke Tsukiyama, Koji Nakano, Yasuaki Ito, Takashi Yazane, Junko Yano, Takumi Kato, Shiro Ozaki, Rie Mori, Ryota Katsuki: **Solving the N-Queens Puzzle by a QUBO Model with Quadratic Size**. Proc. of International Symposium on Computing and Networking, 59-67, 2023. ([DOI](https://doi.org/10.1109/CANDAR60563.2023.00015))