N-QueensパズルのIsingモデルをSymPyで生成し，量子アニーラーで解く

ポイント

• D-Wave量子アニーラーを用いてIsingモデルを与えて解探索する方法
• D-Wave量子アニーラーのchainsgrangthとchain_break_fraction
• D-Wave量子アニーラーに与えるIsingモデルの要求解像度の問題

D-Wave量子アニーラーとIsingモデル

N-QueensパズルのIsingモデルをSymPyで生成し，Gurobiで解くでは、Gurobiを用いて解いていましたが，ここではD-Waveの量子アニーラーに変更します．スピン変数の集合$S$に対して，そのIsingモデルは以下のように表すことができます：
$$
H(S)= \sum_{\binom{S}{2}}^{\lbrace u,v\rbrace} J_{\lbrace u,v\rbrace}uv+\sum_{S}^{u} h_u u
$$
ここで，$\binom{S}{2}=\lbrace\lbrace u,v\rbrace| u,v\in S, u\neq v\rbrace$は2つ選ぶすべての組み合わせ（非順序対）を示し，$J_{\lbrace u,v\rbrace}$は2次項$uv$の係数，$h_u$は1次項$u$の係数を意味します．
N-QueensパズルのIsingモデルをSymPyで生成し，Gurobiで解くでは，N-Queensパズルのクイーンの配置を2次元行列$S=(s_{i,j})$を用いて表現し，これを基にIsingモデルを構築しました．量子アニーラーは，このモデルの$J$と$h$を入力として受け取り，量子アニーリングを繰り返し実行することで解を探索します。

Isingモデルではなく，QUBOモデルを用いることも可能です．しかし，量子アニーラーはQUBOを直接扱うことができないので，D-Waveのツール内では，等価なIsingモデルに自動的に変換後，量子アニーラーで処理されます．モデルの要求解像度を把握するには，Isingモデルを確認する必要があります．そのため，直接Isingモデルを構築することにします。

D-Wave量子アニーラーへの二次項係数Jと一次項係数hの指定方法

D-Wave量子アニーラーにIsingモデルを適用する際，モデルのパラメータである$J$と$h$をPythonの辞書形式で指定する必要があります．

• 二次項係数$J$: スピン変数$u$と$v$から成る集合$\lbrace u,v\rbrace$をキーとし，$J_{\lbrace u,v\rbrace}$を値とする辞書Jを以下のように定義します：
  $$
  \lbrace \lbrace u,v\rbrace :J_{\lbrace u,v\rbrace}\mid J_{\lbrace u,v\rbrace}\neq 0 \rbrace
  $$
• 一次項係数$h$: 変数$u$をキーとし、その係数$h_u$を値とする辞書hを以下のように定義します：
  $$
  \lbrace u :h_u\mid h_u\neq 0 \rbrace
  $$

これらの辞書Jとhを引数にして、量子アニーリングを実行する関数を呼び出します.

Pythonプログラム

self.isingのIsingモデルの各項をもとに，辞書Jとhを作成します．これらの辞書を引数として，sample_ising関数を呼び出して，量子アニーリングを実行します．さらに，chainstrengthとnum_readsというパラメータも指定します．

• chainstrengthは，1つのスピン変数を複数のqubit cellに割り当てる際に，それらの結合（chain）の強さを指定します．この値の適切な設定は、最適解を得る上で非常に重要であり、設定が強すぎるても弱すぎても，望ましい結果が得られません．
• num_readsは，量子アニーリングを行う回数を指定します．このパラメータで指定された回数だけ量子アニーリングを繰り返し実行し，最も良い解を出力します．

#!/usr/bin/env python3

import argparse
import sympy as sp
from dwave.system import DWaveSampler, EmbeddingComposite


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 dwave(self, chainstrength, num_reads):
        self.optimal = -self.ising.pop(1, 0)
        self.chainstrength = chainstrength
        self.num_reads = num_reads
        h = {}
        J = {}
        for key, val in self.ising.items():
            if isinstance(key, sp.Mul):
                quad = key.expand().args
                J[frozenset([quad[0], quad[1]])] = val
            else:
                h[key] = val
        sampler = EmbeddingComposite(DWaveSampler())
        self.properties = sampler.properties["child_properties"]
        self.sampleset = sampler.sample_ising(
            h,
            J,
            chain_strength=chainstrength,
            num_reads=num_reads,
            label=f"{self.n}-Queens puzzle",
        )

    def print(self):
        solution = self.sampleset.first.sample
        timing = self.sampleset.info["timing"]
        for i in range(self.n):
            for j in range(self.n):
                print("-" if solution[self.S[i][j]] < 0 else "+", end="")
            print()
        print(f"optimal={self.optimal} obtained={int(self.sampleset.first.energy)}")
        print(
            f"annealer={self.properties['chip_id']} qubits={len(self.properties['qubits'])} couplers={len(self.properties['couplers'])}"
        )
        print(f"chainstrength={self.chainstrength} num_reads={self.num_reads}")
        print(
            f"qpu_sampling_time={timing['qpu_sampling_time']} qpu_access_time={timing['qpu_access_time']}"
        )
        print(f"chain_break_fraction={self.sampleset.first.chain_break_fraction}")


def main():
    parser = argparse.ArgumentParser(
        description="Solve N-Queens puzzle through SymPy using Quantum annealer"
    )
    parser.add_argument("-s", "--size", default=8, type=int, help="size (default:8)")
    parser.add_argument(
        "-c",
        "--chainstrength",
        default=10,
        type=int,
        help="chainstrength (default:10)",
    )
    parser.add_argument(
        "-n", "--num_reads", default=1000, type=int, help="num_reads (default:1000)"
    )
    args = parser.parse_args()

    nqueen = NQUEEN(args.size)
    nqueen.generate()
    nqueen.dwave(args.chainstrength, args.num_reads)
    nqueen.print()


if __name__ == "__main__":
    main()

実行例

8-Queensパズルを解いた結果，得られた解は$-1264$でした。これは、8-Queensパズルにおける最適解$-1264$に一致しており，正しい解となっています。ただし，繰り返し実行すると，この最適解を得られることもあれば、そうでない場合もあります。

chainstrengthを十分に大きく設定することで，すべてのqubit cellが同じスピン値を取るようになり，結果としてchain_break_fractionが0に近づきます．しかし，これによってqubit cellが取りうる値の自由度が制限され，最適解を見つけることが難しくなる可能性があります．そのため，chain_break_fractionが適度な値（例えば、約$0.1$程度）になるようにchainstrengthの値を調整することが重要です。

--+-----
-----+--
-+------
----+---
-------+
+-------
------+-
---+----
optimal=-1264 obtained=-1264
annealer=Advantage_system4.1 qubits=5627 couplers=40279
h's [38, 42, 46, 50]
J's [2]
chainstrength=10 num_reads=1000
qpu_sampling_time=128800.0 qpu_access_time=144561.97
chain_break_fraction=0.046875

10-Queensパズルの場合，最適解$-2620$に対して，得られた解は$-2588$であり，10-Queensパズルの正しい解ではありません．

--+-------
-------+--
----+-----
-+--------
-----+----
---------+
-----+----
----------
--------+-
----------
optimal=-2620 obtained=-2588
annealer=Advantage_system4.1 qubits=5627 couplers=40279
h's [50, 54, 58, 62, 66]
J's [2]
chainstrength=10 num_reads=1000
qpu_sampling_time=213580.0 qpu_access_time=229339.57
chain_break_fraction=0.27

10-Queensパズルで，chainstrength=15に増やしてみると，解$-2604$となり少し良くなりましたが，それでも最適解ではありません．

--+-----+-
+---------
---+-+----
-------+--
---------+
----------
------+---
----+-----
-+--------
-------+--
optimal=-2620 obtained=-2604
annealer=Advantage_system4.1 qubits=5627 couplers=40279
h's [50, 54, 58, 62, 66]
J's [2]
chainstrength=15 num_reads=1000
qpu_sampling_time=202400.0 qpu_access_time=218163.17
chain_break_fraction=0.06

Isingモデルの要求解像度の問題

8-Queensの場合，hとJの値の最大値は50，最小値は2です．この値の比率$50/2=25$をこのIsingモデルの要求解像度と考えることができます．量子アニーラーは基本的にアナログ処理を行うため，小さな値の差を区別することが容易ではありません．したがって，要求解像度はなるべく小さい方が望ましいとされ，この場合の25はD-Waveの量子アニーラーが扱うことができる解像度の限界に近いと考えられます．