量子アニーリングシミュレータを作ろう - 技術リソース - Amplify - 量子アニーリングと共に進化するクラウド

量子アニーリングの原理では量子アニーリングにおける基本方程式を解説しました。量子アニーリングでは解きたいターゲット関数 \(H_B\) に量子効果 \(H_A\) を加えて、全体として下記のようなハミルトニアン (エネルギー関数) を構成します。

\[ H_{\mathbf{QA}} \left(t\right) = A\left(t\right) H_A + B\left(t\right) H_B \]

両者がアニーリング開始時刻と終時刻で切り替わるように \(A\left(t\right)\) , \(B\left(t\right)\) を制御します。下記とみなせるような時間変化をさせます。

\[ H_{\mathbf{QA}} \left(t = 0\right) = A_{t = 0} H_A \quad \rightarrow \quad H_{\mathbf{QA}} \left(t = \tau\right) = B_{t = \tau} H_B \]

初期状態として適切に基底状態 (最低エネルギー状態) を選べば、ターゲット関数 \(H_B\) の詳細はわからなくても、自動的に最適解が探索されます。

ここでは、具体的な問題に対して量子アニーリングを実行したときに、どのように最適解が探索されるのかを見ていきます。小規模な問題に対し、量子アニーリングシミュレータのコードを作成することで、量子アニーリングに対する理解を深めていきます。

今回は Python3 でコーディングします。数値計算ライブラリとして NumPy, SciPy, Matplotlib を使用しますが、Anaconda で環境構築しているなら動作すると思います。

問題設定: 最小値探索問題 #

量子アニーリングマシンでは、解きたい最適化問題をイジング模型に変換する必要がありました。なぜなら量子ビット間の相互作用はイジング模型で表現されるからです。しかし、基底状態を探索するという量子アニーリングの概念自体では \(H_B\) の詳細構造( \(J_{ij}\) や \(h_i\) の値やグラフ構造)に触れる必要はありません。今回作成するのはあくまでシミュレータなので、物理モデル制約やハードウェアの制約は気にしないことにします。

そこで大胆に、 \(N\) ビット \(2^N\) の通りの入力に対し、「何らかの計算」をした結果である \(2^N\) 通りの出力の中から最小値となる入力を検索する問題を考えます。「何らかの計算」という部分のイジングモデル表現は省略します¹。つまり入力値に対して必ず同じ値を返す未ソートのデータセットの最小値とその入力の探索です。これは当然 \(2^N\) の計算量がかかるので難しい問題といえます。

シミュレータでは、適当なシードによる疑似乱数から生成されたガウス分布とします。

import numpy as np
N = 5
E = np.random.normal(0, scale=N / 2, size=2**N)

ビット数を \(N\) とし、分布幅を \(N\) 程度に抑えるために標準偏差は \(N/2\) としました。サイズ \(2^N\) の配列 \(E\) から最小値を検索する問題です。ターゲット関数の値は本来入力値から計算されますが、今回は省略して \(E[i]\) で得られようにしました。もちろん、自分の解きたい問題に対応する関数 \(E\) を定義しても大丈夫です。

ハミルトニアンの作成 #

これで最適化したいターゲット関数ができました。次に全体のハミルトニアン \(H\) を作成します。今回は以下のようなアニーリングスケジュールに従うことにしました。

\[ H = \left( 1 - \frac{t}{\tau} \right) H_A + \frac{t}{\tau} H_B \]

アニーリング開始時刻 \(t = 0\) では \(H_A\) のみ、終了時刻 \(t=\tau\) では \(H_B\) のみになるような時間変化になっています。量子効果 \(H_A\) として横磁場を導入します。

\[ H_A = - \sum_{i}^N{\sigma_i^x} \]

ターゲット関数 \(H_B\) は前述の通りですが、形式的に \(N\) ビット( \(N\) スピン) を入力として持ちます。

\[ H_B = E(\sigma_0^z, \sigma_1^z, \cdots, \sigma_{2^N-1}^z) \]

ハミルトニアン \(H\) は \(2^N \times 2^N\) の行列です。これの作り方をおさらいすると、横磁場の場合には行 \(i\) 列 \(j\) に対し

対角要素:

\(H[i, i]=B(t) E(i)\)

非対角要素:

\(H[i, j]= - A \left( t \right)\) ( \(j\) は \(i\) の1ビット反転)
\(H[i, j]=0\) (上記以外)

と設定します。非対角要素を埋めるときの \(j\) は \(i\) のどのビットを反転するかで \(N\) 通りあります。

時刻 \(t\) における横磁場イジング模型のハミルトニアンを生成する関数 create_tfim は次のように定義されます。

Tau = 1

def scheduleE(time):
return time / Tau

def scheduleG(time):
return (Tau - time) / Tau

def create_tfim(time=0, hamiltonian=None):
# Set diagonal part
v = scheduleE(time)
if hamiltonian is None:
hamiltonian = np.diag(v * E)
else:
for i in range(2**N):
hamiltonian[i, i] = v * E[i]

# Set off-diagonal part
g = - 1 * scheduleG(time)
for i in range(2**N):
for n in range(N):
j = i ^ (1 << n)
hamiltonian[i, j] = g
return hamiltonian

scheduleG と scheduleE はそれぞれ \(A\left(t \right)\) , \(B\left( t \right)\) に対応します。例えば \(t / \tau = 0.5\) のハミルトニアンを print すると下記のように出力されます。

[[ 0.83337305 -0.5        -0.5        ...,  0.          0.          0.        ]
 [-0.5         1.47158252  0.         ...,  0.          0.          0.        ]
 [-0.5         0.          1.13847563 ...,  0.          0.          0.        ]
 ...,
 [ 0.          0.          0.         ..., -1.52027369  0.         -0.5       ]
 [ 0.          0.          0.         ...,  0.          0.92076888 -0.5       ]
 [ 0.          0.          0.         ..., -0.5        -0.5        -0.08976278]]

以上でシミュレータの入力に対する準備が整いました。

断熱的に時間変化させる場合 #

断熱的、すなわち極めてゆっくりと時間変化させる場合には、各時刻 \(t\) においてハミルトニアンのスナップショットを取ることで、時間に依存していないものとみなせます。この場合、各時刻の状態とエネルギーは時間に依存しないシュレディンガー方程式 \( H \left| \psi \right\rangle = E \left| \psi \right\rangle \) を解くことで求められます。

それでは実際に数値的にこの固有方程式を解いてみましょう。

from scipy.linalg import eigh

H = None
step = 0.01
time_steps = [step * i for i in range(int(Tau / step) + 1)]
for i, t in enumerate(time_steps):
H = create_tfim(t, H)
evals_all, evecs_all = eigh(H)
# ADDITIONAL CODE HERE

時刻を step ずつ動かしながら各時刻の固有値 evals_all と固有ベクトル evecs_all を求めていきます。今回は SciPy の eigh を使用しました。各時刻において状態は必ず evec_all の線形結合になっていることをシュレディンガー方程式は意味します。

各時刻 \(t/\tau\) を横軸に、固有値(エネルギー)をプロットすると下図(上半分)のようになります。

エネルギー図

量子アニーリングでは初期状態 \(t = 0\) を基底状態(最低エネルギー状態)にセットします。実は、この方程式の \(t = 0\) の解、つまり \(H = -\sum_{i}\sigma_i^x\) の場合は解析的に解くことができて \(E = -N\) , \(\left| \psi \right\rangle = \left( 1/2^{N/2}, \cdots \right)\) です。この状態から量子アニーリングを開始すると、図の一番下の青線を各時刻で辿り、最終的に \(t/\tau = 1\) に到達します。 \(t/\tau = 1\) はつまり解きたいターゲット関数なので最適解に到達することになります。その後、スピンの状態を観測することで最適解が得られます。

次に、時間変化とともにどのように最適解が強調されていくかを見ていきます。各時刻の固有ベクトル evecs_all[:, i]² はスピン状態 \(\left| \psi \right\rangle\) そのものを表します。一般に固有ベクトルは複素数ですが、ベクトル成分を確率に変換するには、絶対値の2乗を取るという約束でした。

そこで、下記のような関数を用意して各時刻の全状態の確率を記録してみます。

def amp2prob(vec):
p = [np.abs(z)**2 for z in vec]
return np.array(p)

amp2prob の引数に evecs_all[:, 0]³ を与えると、スピンの取りうる全 \(2^N\) 状態それぞれについて観測確率が得られます。これをプロットしたのが図の下半分パートです。初期状態では全状態が等確率になっていますが、時間変化とともに偏りが現れます。図の場合では、 \(E[23]\) が最適解だったようですが、最終的に \(23\) に対応する状態の確率が1に収束する様子が見られます。

このように極めてゆっくりと量子アニーリングを実行すると、必ず最適解が得られることになります。

アニーリングスケジュールに従う場合 #

次に限られた時間でのアニーリング、つまり現実的な性能について見ていきます。時間変化するハミルトニアンに対し、ある初期状態からの時間発展は、時間に依存するシュレディンガー方程式

\[ \imath \hbar \frac{\mathrm{d} }{\mathrm{d} t}\left| \psi \left (t \right ) \right\rangle = H \left (t \right ) \left| \psi \left (t \right ) \right\rangle \] に従います。 \(H\) は行列、 \(\left| \psi \right\rangle\) はベクトル、左辺に虚数単位 \(\imath\) があるので、この方程式は複素数の連立微分方程式です。 \(\hbar\) はプランク定数と呼ばれますが、以後 \(1\) とします。

数値積分ライブラリとして SciPy の complex_ode を使用します。complex_ode は微分方程式の右辺を配列で返す関数を引数に取ります。そのため、まずはシュレディンガー方程式を少しだけ変形します。

\[ \frac{\mathrm{d} }{\mathrm{d} t}\left| \psi \left (t \right ) \right\rangle = - \imath H \left (t \right ) \left| \psi \left (t \right ) \right\rangle \]

ハミルトニアン \(H\) はスパースなので、非ゼロの要素だけに着目すれば十分です。右辺の \(H \left (t \right ) \left| \psi \left (t \right ) \right\rangle\) について、行 \(i\) に対する計算を行う関数を下記のように定義しました。

def genrate_diffeq(idx, offdiag_indices, time, vec):
v = E[idx] * scheduleE(time)
g = -1 * scheduleG(time)
return v * vec[idx] + g * np.sum(vec[offdiag_indices])

この関数は非ゼロ要素のインデックス配列を与えることで行列要素全てに対する演算を回避しています。ハミルトニアンの非ゼロ要素は \(N\) から自動的に決まるため、事前に計算しておきます。

diffeq_array = [lambda vec, t, _i=i, _ij=np.array([i ^ (1 << n) for n in range(N)]):
genrate_diffeq(_i, _ij, t, vec) for i in range(2**N)]

以上より、complex_ode に与える微分方程式の右辺は次の関数で定義できます。

def simdiffeq_rhs(t, vec):
return np.array([- 1j * f(vec, t) for f in diffeq_array],
dtype=np.complex_)

後は初期状態をセットして数値積分を計算することで量子アニーリングのシミュレーションが実行されます。

from scipy.integrate import complex_ode

vec0 = np.array([2**(-N / 2)] * 2**N, dtype=np.complex_)
r = complex_ode(simdiffeq_rhs)
r.set_initial_value(vec0)

steps = 100
i = 0
t = 0
while r.successful() and i <= steps:
v = r.integrate(t) if t != 0 else vec0
# ADDITIONAL CODE HERE
t += Tau / steps

上では一度に積分せずに、Tau/steps ごとに状態の解析が行えるようにしています。以上で実装はほぼ完了です。あとは何を測定するかというコードを埋め込んでいきます。

アニーリングのプロセスと性能評価 #

有限時間のアニーリングでは、一般にアニーリング時間 \(\tau\) が長いほど最適解が得られる確率 \(p\) が高くなります。ここでは量子アニーリングの性能評価についてどのような要素が影響をあたえるのか、数値実験を行いながら見ていきます。

ところで、量子アニーリングで最適解を要求する場合でも、必ずしも確率 \(p=1\) で最適解に到達する必要はありません。なぜなら、試行回数を \(1/p\) とすることで1回程度は最適解に到達することを期待できるからです。また \(p=1\) という制約は、一般的には断熱的と見なせる程度のアニーリング時間が必要となり、実用性からも現実的ではありません。

また、必ずしも最適解を必要としないという場面もあります。実用的には最適解に近い解であれば十分という場面も多く、常に厳密解が必要という状況は特殊だと思われます⁴。

以上の二つの視点に基づき、

どの程度の確率で最適解が得られるか
最適解からの平均的な距離はどの程度か

について調査します。

最適解の得られる確率 #

まずは \(\tau\) を変えながら最適解の得られる確率がどのように変化するかを調査します。上記の数値積分を実行した後に、最適解に対応するインデックス i に対して、v[i] の絶対値の二乗をとることで、最適解の得られる確率が得られます。下のコードのようにターゲット関数 E を固定した状況で Tau でループさせた時の到達確率は次の図のようになりました。

i = np.argmin(E)
vec0 = np.array([2**(-N / 2)] * 2**N, dtype=np.complex_)
for Tau in [1 / 4, 1 / 2, 1, 2, 3, 4, 6, 8, 12, 16, 24, 32, 48, 64, 96,
128, 192, 256, 384, 512, 768, 1024, 1536, 2048, 3072, 4096]:
#
# 数値積分コード
#
print(Tau, abs(v[i])**2)

アニーリング時間と性能

\(\tau\) の値そのものについて評価を行うのは難しいですが、ある程度長時間のアニーリングを実行すると、急速に \(p\) が \(1\) に近づく様子が見られます。どの程度時間をかけるか、試行回数を増やすかというのはトレードオフの関係にありますが、問題とマシンの特性に応じて選ぶ必要があると思われます。

今のところ実際の量子コンピュータにはコヒーレンス時間 (量子ビットの寿命) の制約があるため、長時間のアニーリングは難しいようです。一方で試行回数は例えば 1000 回などに設定できるようなので、このグラフで言うと左側の特性が対応するように思います。

基底状態とのオーバラップとエネルギーギャップ #

次にアニーリングの途中のプロセスを見ていきます。常に基底状態をたどるような過程であれば最適解に到達しますが、短い時間のアニーリングでは励起状態 (最低エネルギー状態以外) に遷移していると思われます。

現在の状態がどの程度基底状態と重なり合っているかは、各時刻でのハミルトニアンの固有状態との内積を取ることで測定できます。先程の数値積分コードを少し変更し下記のようにします。

vec0 = np.array([2**(-N / 2)] * 2**N, dtype=np.complex_)
r = complex_ode(simdiffeq_rhs)
r.set_initial_value(vec0)

steps = 100
i = 0
t = 0
while r.successful() and i <= steps:
v = r.integrate(t) if t != 0 else vec0
H = create_tfim(t)
evals, evecs = eigh(H, eigvals=(0, 4))          # 固有値を下から5つ取得
overlap = np.abs(np.vdot(evecs[:, 0], v))**2    # 基底状態と v の内積の絶対値2乗
probs = amp2prob(v)                             # 全スピン状態各々の観測確率
t += Tau / steps

複素ベクトルなので、内積では片方の複素共役を取る必要があります。その後確率に変換するために絶対値の2乗を行います。これでアニーリング途中で時刻 Tau / steps 毎にどの程度基底状態の留まることが出来たかを評価できます。

アニーリングの途中において、時間とともにどのように状態確率が変化するかを表したのが下の図です。上から \(\tau = 1\) , \(8\) , \(64\) , \(512\) , \(2048\) の結果を並べています。

QA tau = 1 QA tau = 8 QA tau = 64 QA tau = 512 QA tau = 2048

アニーリング時間 \(\tau\) が \(1\) 程度に小さい場合には、基底状態から多くの状態に遷移してしまってます。しかし \(\tau = 8\) 程度まで増やすと、最適解から数えて下から4状態のいずれかに収束するようになります。 \(\tau = 64\) では下から2状態、 \(\tau = 512\) 以上では最適解が最も高い到達確率になります。

これら様々な \(\tau\) について共通する性質は、大体 \(0.6 \le t/\tau \le 0.7\) の間で、基底状態とのオーバーラップに大きな変化があることです。実は、基底状態から励起状態への遷移確率はハミルトニアンの固有値の差が関係しています。

下の図は今回のハミルトニアンで \(t/\tau\) を変えながら固有値を低い順に5つを計算した結果です。

エネルギーギャップ

最低エネルギー状態と一つ上の状態の差(エネルギーギャップ)は \(0.6 \le t/\tau \le 0.7\) 付近で最も小さくなっています。エネルギーギャップ \(\Delta E\) と必要なアニーリング時間の間には、 \(T \propto \Delta E^{-2}\) の関係があることが知られています。もし事前にギャップ位置がわかっていればそこでじっくりと時間をかけてアニーリングすることで、精度を向上することもできます。

残留エネルギー #

先ほどの \(\tau = 1\) , \(8\) , \(64\) , \(512\) , \(2048\) の結果において、最適解が得られる確率を高めるには \(\tau = 512\) 程度以上の時間をかける必要がありました。しかし、誤差が許容できる場合はどうでしょうか。例えば最適解から数えて下から4状態で問題なければ、長時間かけなくとも \(\tau = 8\) で良いわけです。厳密性に目をつぶり、近似ソルバーとして見なした場合には、どの程度最適解に近づくことができたかという指標が重要になります。

そのような指標として、到達した状態に対するエネルギーの期待値と最適解のエネルギーの差 (残留エネルギー) で評価してみます。

i = np.argmin(E)
vec0 = np.array([2**(-N / 2)] * 2**N, dtype=np.complex_)
for Tau in [1 / 4, 1 / 2, 1, 2, 3, 4, 6, 8, 12, 16, 24, 32, 48, 64, 96,
128, 192, 256, 384, 512, 768, 1024, 1536, 2048]:
#
# 数値積分コード
#
print(Tau, np.dot(E, [np.abs(z)**2 for z in v]) - E[i])

残留エネルギー

予想通り \(\tau = 8\) 程度で残留エネルギーが十分小さくなる様子が見られます。もちろんどの程度まで最適解に近づける必要があるのかは状況次第ですが、これ以上 \(\tau\) を大きくしても割に合わないという意味では、量子アニーリングの限界ともいえます。繰り返しになりますが、試行回数とアニーリング時間、そして許容誤差のバランスに応じて適切なパラメータ設定を考察する必要があると考えられます。

まとめ #

本ページでは量子アニーリングの基本方程式に従うシミュレータを作成し、理想的な量子アニーリングマシンがどのように最適解を探索しているのかを見てきました。極めてゆっくりとアニーリングを実行した場合は必ず最適解を得ることも可能ですが、現実的な時間で実行を完了するためには、

アニーリング時間
試行回数
最適解からの許容誤差

を設定する必要があります。これらは目的や問題設定、ハードウェアの制約や特性に応じて決めるべき値です。また、今回は深掘りしませんでしたが、どのように量子効果を導入し制御するか、そしてスケジュールについても調整することが可能です。

これらは理論・数値実験の双方、そして実際のマシンを使って活発に議論されており、量子アニーリングの性能向上につながる研究として今後の発展が期待されています。

本来はどのような論理式でもイジング模型で表現できるはずです ↩︎
eigh から返却される固有ベクトルは転置されています ↩︎
固有値順にソートされているものとします ↩︎
決定問題をアニーリングで解くという状況は厳密解が必要になります ↩︎