ポイント
- SymPyでスピン変数($-1/+1$)を扱う方法
- Gurobiでスピン変数を扱う方法
N-QueensパズルのIsingモデル化
N-QueensパズルのQUBOをSymPyで生成し，Gurobiで解くのIsingモデル版です．
QUBOからIsingモデルへは自動的に変換可能ですが，ここでは敢えてその変換を使用せずに，Isingモデルを直接構築します．
Isingモデルは，スピン値（$-1$および$+1$）を取る変数の二次式です．$n \times n$のスピン変数$S = (s_{i,j})$を用いて，$n \times n$のボードを表現します．ここで，$s_{i,j} = +1$は，座標$(i,j)$にクイーンが配置されていることを意味します．各行および各列に$+1$が一つだけ存在する（one-hot制約）場合，さらに，各対角線上で$+1$が0個または1個だけ存在する（zero-one-hot制約）場合，変数$S$はN-Queensパズルの解を表しています．
Ising modelのone-hot制約とzero-one-hot制約
$n$個のスピン変数$s_0, s_1, \ldots, s_{n-1}$に対して、one-hot制約が成立するのは以下を満たした時です．
$$
s_0+s_1+\cdots +s_{n-1}=1\times (+1)+ (n-1)\times(-1)=-n+2
$$
これは，式
$$
(s_0+s_1+\cdots +s_{n-1}+n-2)^2
$$
が最小値0を取る場合に相当します．また、zero-one-hot制約を満たすのは，さらに以下を満たす時です．
$$
s_0+s_1+\cdots +s_{n-1}=0\times (+1)+ n\times(-1)=-n
$$
つまり，式
$$
(s_0+s_1+\cdots +s_{n-1}+n)(s_0+s_1+\cdots+ s_{n-1}+n-2)
$$
が最小値0をとる場合に該当します．したがって、全ての行と列がone-hot制約を満たし、さらに全ての斜めラインがzero-one-hot制約を満たす二次式を合計することで，N-Queensパズルを解くためのIsingモデルを構築できます．
Gurobiでスピン変数を扱う方法
Gurobi Optimizerでは，スピン変数を直接扱う機能は（たぶん）提供されていません．そのため、$-1$と$+1$の値のみを取る整数変数を定義するために，2値変数を介する必要があります．具体的には，2値変数$x_{i,j}$と整数変数$s_{i,j}$の間に，次の制約を設定します。
$$
2x_{i,j}-s_{i,j}=1
$$
この制約により，$x_{i,j}$が0の場合$s_{i,j}$は$-1$となり，$x_{i,j}$が1の場合は$s_{i,j}$が$+1$になります．
Pythonプログラム
SymPyを利用してスピン変数を扱うために，sp_Spinクラスを定義しています．このクラスでは，スピン変数が指数演算により偶数乗されると1になるという性質を利用しています．self.Sを$S$を保持するのに用います．HとVはそれぞれself.Sの行と列にあるスピン変数を表し，AとBはself.Sの斜めライン上のスピン変数を表しています．HとVの各要素がone-hot制約を，AとBの各要素がzero-one-hot制約を満たす二次式の合計がformulaとなります．このformulaを展開して得られる項のリストがself.isingです．
Gurobiを使用してこのIsingモデルから解を求めるために，varXとvarSをGurobiの変数として用います．ここで，varX[i,j]は2値変数で$x_{i,j}$に対応し，varS[self.S[i][j]]は最小値が$-1$，最大値が$+1$の整数変数で$s_{i,j}$に対応します．最小値の指定は重要で，指定しない場合Gurobiはデフォルトで$0$が最小値となるため，明示的に$-1$を設定します。$2x_{i,j}-s_{i,j}=1$の制約をaddConstrsメソッドを使用して追加します．最後に，self.isingに含まれる各一次項と二次項を，objにvarSの項として加えていきます．
#!/usr/bin/env python3
import argparse
import sympy as sp
import gurobipy as gp
class sp_Spin(sp.Symbol):
def _eval_power(self, exp):
if exp % 2 == 0:
return sp.Integer(1)
else:
return self
class NQUEEN:
def __init__(self, n):
self.n = n
self.S = [[sp_Spin(f"s_{i}_{j}") for j in range(self.n)] for i in range(self.n)]
def generate(self):
H = [[self.S[i][j] for j in range(self.n)] for i in range(self.n)]
V = [[self.S[i][j] for i in range(self.n)] for j in range(self.n)]
A = [
[self.S[i + k][k] for k in range(self.n) if 0 <= i + k and i + k < self.n]
for i in range(-self.n, 2 * self.n)
]
B = [
[self.S[i - k][k] for k in range(self.n) if 0 <= i - k and i - k < self.n]
for i in range(-self.n, 2 * self.n)
]
formula = (
sum([(sum(x) + len(x) - 2) ** 2 for x in H])
+ sum([(sum(x) + len(x) - 2) ** 2 for x in V])
+ sum([(sum(x) + len(x)) * (sum(x) + len(x) - 2) for x in A if len(x) >= 2])
+ sum([(sum(x) + len(x)) * (sum(x) + len(x) - 2) for x in B if len(x) >= 2])
)
self.ising = formula.expand().as_coefficients_dict()
def gurobi(self, timelimit):
self.model = gp.Model()
self.model.setParam("TimeLimit", timelimit)
varX = {}
varS = {}
for i in range(self.n):
for j in range(self.n):
varX[i, j] = self.model.addVar(vtype=gp.GRB.BINARY, name=f"x_{i}_{j}")
varS[self.S[i][j]] = self.model.addVar(
lb=-1.0, ub=1.0, vtype=gp.GRB.INTEGER, name=f"s_{i}_{j}"
)
self.optimal = -self.ising.pop(1, 0)
self.model.addConstrs(
(
2 * varX[i, j] - varS[self.S[i][j]] == 1
for i in range(self.n)
for j in range(self.n)
)
)
obj = []
for key, val in self.ising.items():
if isinstance(key, sp.Mul):
quad = key.expand().args
obj.append(val * varS[quad[0]] * varS[quad[1]])
else:
obj.append(val * varS[key])
self.model.setObjective(gp.quicksum(obj), sense=gp.GRB.MINIMIZE)
self.model.optimize()
def print(self):
for i in range(self.n):
for j in range(self.n):
print(
"-" if self.model.getVarByName(f"s_{i}_{j}").X < 0 else "+", end=""
)
print()
print(f"optimal={self.optimal} obtained={int(self.model.ObjVal)}")
def main():
parser = argparse.ArgumentParser(
description="Solve N-Queens puzzle through SymPy using Gurobi"
)
parser.add_argument(
"-n", default=8, type=int, help="Value of N of N-QUEENS (default:8)"
)
parser.add_argument(
"-t", default=10, type=int, help="Time limit of Gurobi optimizer"
)
args = parser.parse_args()
nqueen = NQUEEN(args.n)
nqueen.generate()
nqueen.gurobi(args.t)
nqueen.print()
if __name__ == "__main__":
main()
実行例
8-QueensパズルをSymPyでIsingモデル化してGurobiで解くと，値$-1264$の最適解が得られ，8-Queensパズルの正しい解になっています．
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (linux64)
Thread count: 16 physical cores, 32 logical processors, using up to 32 threads
Optimize a model with 64 rows, 128 columns and 128 nonzeros
Model fingerprint: 0x984ca412
Model has 728 quadratic objective terms
Variable types: 0 continuous, 128 integer (64 binary)
Coefficient statistics:
Matrix range [1e+00, 2e+00]
Objective range [4e+01, 5e+01]
QObjective range [4e+00, 4e+00]
Bounds range [1e+00, 1e+00]
RHS range [1e+00, 1e+00]
Found heuristic solution: objective -1200.000000
Presolve removed 64 rows and 64 columns
Presolve time: 0.00s
Presolved: 728 rows, 792 columns, 2184 nonzeros
Variable types: 0 continuous, 792 integer (792 binary)
Root relaxation: objective -1.456000e+03, 64 iterations, 0.00 seconds
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 -1456.0000 0 64 -1200.0000 -1456.0000 21.3% - 0s
H 0 0 -1240.000000 -1456.0000 17.4% - 0s
H 0 0 -1248.000000 -1456.0000 16.7% - 0s
H 0 0 -1256.000000 -1456.0000 15.9% - 0s
H 0 0 -1264.000000 -1436.0000 13.6% - 0s
0 0 -1370.6667 0 64 -1264.0000 -1370.6667 8.44% - 0s
0 0 -1304.0000 0 64 -1264.0000 -1304.0000 3.16% - 0s
0 2 -1304.0000 0 64 -1264.0000 -1304.0000 3.16% - 0s
Cutting planes:
Gomory: 1
MIR: 354
Zero half: 149
RLT: 289
BQP: 227
Explored 111 nodes (2489 simplex iterations) in 0.69 seconds
Thread count was 32 (of 32 available processors)
Solution count 5: -1264 -1256 -1248 ... -1200
No other solutions better than -1264
Optimal solution found (tolerance 1.00e-04)
Best objective -1.264000000000e+03, best bound -1.264000000000e+03, gap 0.0000%
----+---
-+------
---+----
------+-
--+-----
-------+
-----+--
+-------
optimal=-1264 obtained=-1264