LoginSignup
2
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@koji_nakano(中野 浩嗣)

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

Posted at

ポイント

  • 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の量子アニーラーが扱うことができる解像度の限界に近いと考えられます．

2
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
2
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?