量子アニーラーのシミュレーション

量子アニーラーは、量子ビットの間の結合が決まった形のグラフの上で、イジングモデルのエネルギーが小さくなるスピンの割り当てを探す装置です。同じイジングモデルを PyQBPP で式にして解けば、量子アニーラーが返すはずの解が得られます。つまり、実機を使わずに量子アニーラーの動作をシミュレーションできます。実機に問題を送る前にプログラムを確かめたり、実機が返した解を最適解と比べたりするのに使えます。

このページでは、D-Wave の量子アニーラー 2000Q の結合の形であるキメラグラフの上で、実機と同じ範囲の乱数の係数を与えたイジングモデルを解きます。 2000Q はすでに提供が終わっていて、現在の D-Wave の量子アニーラーは Advantage2 です。ここでは構造が分かりやすい 2000Q のキメラグラフを使いますが、Advantage2 にも同じ方法で対応できます（Advantage2 の場合）。

イジングモデル

量子ビット $i$ はスピン $s_i\in\lbrace -1,+1\rbrace$ を表します。量子ビットを頂点、結合を辺とするグラフの辺の集合を $E$ とすると、イジングモデルのエネルギーは次の式です：

\[E(s) = \sum_{i} h_i s_i + \sum_{(i,j)\in E} J_{ij}\, s_i s_j\]

$h_i$ は量子ビット $i$ にかかる外場、$J_{ij}$ は量子ビット $i$ と $j$ の結合の強さで、どちらも実数です。 $J_{ij}<0$ の結合は $s_i$ と $s_j$ が同じとき、$J_{ij}>0$ の結合は異なるときに、エネルギーを下げます。利用者は $h_i$ と $J_{ij}$ を与え、量子アニーラーは $E(s)$ の小さいスピンの割り当てを返します。 2000Q で与えられる値の範囲は、$-2\le h_i\le 2$、$-1\le J_{ij}\le 1$ です。

キメラグラフ

2000Q の量子ビットは、キメラグラフの形に結合しています。キメラグラフは、8 量子ビットからなるセルを格子状に並べたグラフです。セルの 8 量子ビットは 4 個ずつの 2 組（0 組目と 1 組目）に分かれていて、次のように結合します：

セルの中では、0 組目の各量子ビットが 1 組目の 4 量子ビットすべてと結合します。
0 組目の量子ビットは、下のセルの 0 組目の同じ位置の量子ビットと結合します。
1 組目の量子ビットは、右のセルの 1 組目の同じ位置の量子ビットと結合します。

2000Q は 16×16 個のセル（2048 量子ビット）でした。このページでは、2×2 個のセル（32 量子ビット、80 個の結合）を使います：

2×2 セルのキメラグラフ

薄い灰色の四角がセルで、黄色の 4 量子ビットが 0 組目、灰色の 4 量子ビットが 1 組目です。 0 組目は横に並べてあり、下のセルと縦の辺で結合します。 1 組目は縦に並べてあり、右のセルと横の辺で結合します。

$i$ 行 $j$ 列のセルの $u$ 組目の $k$ 番目（$k=0,1,2,3$）の量子ビットには、番号 $((i\cdot C+j)\cdot 2+u)\cdot 4+k$ を付けます（$C$ はセルの列数）。この番号の付け方は、D-Wave の Ocean SDK が作るキメラグラフと同じです。

QUBO への変換

PyQBPP のソルバーは、値が 0 か 1 のバイナリ変数の式を解きます。スピン $s_i$ とバイナリ変数 $x_i$ を $s_i=2x_i-1$ で対応させると（$x_i=1$ のとき $s_i=+1$、$x_i=0$ のとき $s_i=-1$）、$E(s)$ はバイナリ変数の式になります。 PyQBPP では、スピン変数で $E(s)$ を書いて spin_to_binary() を呼ぶと、この置き換えができます。置き換えた式の値は元の $E(s)$ と同じなので、ソルバーが返すエネルギーはそのままイジングモデルのエネルギーです。

PyQBPP プログラム

以下の PyQBPP プログラムは、2×2 個のセルのキメラグラフで、乱数の係数を与えたイジングモデルを解きます：

import random
import pyqbpp.d as qbpp


def chimera(rows, cols):
    L = 4  # セルの各組の量子ビット数

    def q(i, j, u, k):
        return ((i * cols + j) * 2 + u) * L + k

    nodes = list(range(rows * cols * 2 * L))
    edges = []
    for i in range(rows):
        for j in range(cols):
            for k in range(L):
                for l in range(L):
                    edges.append((q(i, j, 0, k), q(i, j, 1, l)))
                if i + 1 < rows:
                    edges.append((q(i, j, 0, k), q(i + 1, j, 0, k)))
                if j + 1 < cols:
                    edges.append((q(i, j, 1, k), q(i, j + 1, 1, k)))
    return nodes, edges


ROWS, COLS = 2, 2  # セルの行数と列数
nodes, edges = chimera(ROWS, COLS)

random.seed(1)
h = {v: random.uniform(-2, 2) for v in nodes}
J = {e: random.uniform(-1, 1) for e in edges}

s = qbpp.var("s", shape=len(nodes))
f = 0
for v in nodes:
    f += h[v] * s[v]
for (a, b), w in J.items():
    f += w * s[a] * s[b]
f.spin_to_binary()
f.simplify_as_binary()

solver = qbpp.EasySolver(f)
sol = solver.search(time_limit=1)

spin = {v: 2 * sol(s[v]) - 1 for v in nodes}
E = 0
for v in nodes:
    E += h[v] * spin[v]
for (a, b), w in J.items():
    E += w * spin[a] * spin[b]
print(f"energy = {sol.energy:.4f}")
print(f"E(s)   = {E:.4f}")
for i in range(ROWS):
    for j in range(COLS):
        v = (i * COLS + j) * 8
        left = " ".join(f"{spin[v + k]:+d}" for k in range(4))
        right = " ".join(f"{spin[v + k]:+d}" for k in range(4, 8))
        print(f"cell ({i},{j}): {left} | {right}")

このプログラムでは、次のことを行っています：

import pyqbpp.d で、式の係数を実数（float）にしています（実数（double）係数）。
関数 chimera(rows, cols) は、rows 行 cols 列のセルからなるキメラグラフを作り、ノード（量子ビットの番号）のリスト nodes と辺のリスト edges を返します。中の q(i, j, u, k) が量子ビットの番号を計算し、セルの中の結合、下のセルとの結合、右のセルとの結合の 3 種類の辺を edges に入れます。
外場 h[v] を $-2$ 以上 $2$ 以下、結合 J[e] を $-1$ 以上 $1$ 以下の一様乱数で決めます。どちらも、ノードと辺ごとの値を辞書で持ちます。
スピン変数 s で $E(s)$ を式 f にし、f.spin_to_binary() でバイナリ変数の式に変えてから、f.simplify_as_binary() で簡約します。式を作る部分は、nodes と係数 h、J しか使っていません。
EasySolver で、制限時間を 1 秒として解きます。
sol(s[v]) はバイナリ変数の値（0 か 1）なので、2 * sol(s[v]) - 1 でスピン（$+1$ か $-1$）に戻します。そのスピンで $E(s)$ を計算し直し、ソルバーのエネルギー sol.energy と並べて表示します。
最後に、セルごとのスピンを表示します。$i$ 行 $j$ 列のセルの量子ビットは、$8(i\cdot C+j)$ 番からの 8 個です。

このプログラムの出力は以下のとおりです：

energy = -46.2007
E(s)   = -46.2007
cell (0,0): +1 -1 -1 -1 | +1 +1 +1 +1
cell (0,1): +1 +1 -1 +1 | +1 +1 -1 -1
cell (1,0): +1 -1 -1 +1 | +1 +1 -1 -1
cell (1,1): +1 +1 +1 +1 | -1 -1 +1 +1

energy と E(s) が一致していて、ソルバーのエネルギーがイジングモデルのエネルギーそのものであることが分かります。各行は 1 つのセルのスピンで、| の左が 0 組目、右が 1 組目です。量子アニーラーと同じように、各量子ビットのスピンと、そのエネルギーが得られました。

次の 2 つの図は、このプログラムで作った問題とその解です。 1 つ目の図は問題、つまり乱数で決めた $h_i$ と $J_{ij}$ を表しています：

キメラグラフ上の乱数の係数 h と J

量子ビットの色は $h_i$ を表します。赤は $h_i>0$、青は $h_i<0$ で、色が濃いほど $|h_i|$ が大きくなります。項 $h_is_i$ は $s_i$ が $h_i$ と逆の符号のときに負になるので、赤い量子ビットは $-1$ に、青い量子ビットは $+1$ になろうとします。
緑の実線は $J_{ij}<0$ の結合で、両端のスピンをそろえようとします（引き合う）。紫の破線は $J_{ij}>0$ の結合で、両端のスピンを逆向きにしようとします（反発する）。線の太さは $|J_{ij}|$ の大きさを表します。

2 つ目の図は、この問題の解です：

キメラグラフ上のイジングモデルの解

黒がスピン $+1$、白がスピン $-1$ の量子ビットです。 32 量子ビットのうち 23 個は $h_i$ が望む向き（赤なら白、青なら黒）になっていますが、量子ビット 7 や 12 のように、結合の影響で逆向きになっているものもあります。橙の太線は満たされない結合、つまり $J_{ij}s_is_j>0$ となっていて、その項だけを見るとエネルギーを上げている結合です。すべての結合を同時に満たす割り当ては、一般には存在しません。たとえば閉路をたどりながら、結合ごとに望む関係（同じか異なるか）になるようにスピンを決めていくと、1 周して戻ったときに向きが食い違うことがあります。外場 $h_i$ が望む向きと、結合が望む向きがぶつかることもあります。最適解は、こうした損の合計が最小になる割り当てで、この例では 80 個の結合のうち 22 個が満たされていません。

全探索で確かめる

32 量子ビットのスピンの割り当ては $2^{32}$（約 43 億）通りなので、Exhaustive Solver ですべて調べられます。ソルバーを作って解く 2 行を次のように変えます：

solver = qbpp.ExhaustiveSolver(f)
sol = solver.search()

出力は EasySolver のときと同じで、EasySolver の解がエネルギー最小の割り当て（最適解）であることが確かめられます。理想的な量子アニーラーが返すのも、この最適解です。全探索にかかる時間は、GPU（A100）1 枚で約 2 秒、CPU 8 コアだけで約 30 秒でした。

大きなキメラグラフ

ROWS と COLS を 16 にすると、2000Q と同じ 2048 量子ビット、6016 個の結合になります。この大きさでは全探索はできませんが、EasySolver でそのまま解けます。ただし、制限時間 1 秒では最適解に届かないことがあり、制限時間を延ばすとよりエネルギーの小さい解が得られることがあります。

Advantage2 の場合

Advantage2 の量子ビットは、キメラグラフより結合の多い Zephyr グラフの形に結合しています（1 量子ビットあたりの結合は、キメラグラフでは最大 6 個、Zephyr グラフでは最大 20 個）。 $h_i$ と $J_{ij}$ の範囲も 2000Q とは異なります。

このプログラムで式を作って解く部分は、chimera() が返す nodes、edges と係数 h、J しか使っていません。そのため、chimera() の代わりに Zephyr グラフの nodes と edges を使い、乱数の範囲を Advantage2 の範囲に変えれば、同じプログラムで Advantage2 のシミュレーションになります。 D-Wave の Ocean SDK には、キメラグラフや Zephyr グラフのノードと辺のリストを作る関数があります。実機には使えない量子ビットや結合があることもあるので、実機の量子ビットと結合の一覧と $h_i$、$J_{ij}$ の範囲を Ocean SDK で取得して使えば、その実機と同じ形の問題を PyQBPP で解けます。このとき量子ビットの番号は飛び飛びになるので、変数の配列 s は「最大の番号 + 1」の大きさで作ります（式に現れない変数は解に影響しません）。