量子アニーリング

ここでは、横磁場の無い「古典イジング模型」の基底状態をどうやって求めるのかを考えてみましょう。

(適当に書いているので間違いなどあると思うので鵜呑みにしないようにしてください)

横磁場の無い所謂 「古典イジング模型」 の基底状態を実時間発展で求めるにはどうすればよいでしょうか?

例えば4 qubit状態を考えて |0000> から実時間発展させて見ましょう。

実時間発展でこの状態が発展してエネルギーが最小になる状態に落ち着いてくれればよいのですが、実はそうはなりません。

これはこの状態自体がハミルトニアンの固有状態になっており時間発展してくれないためです。

(横磁場がある場合は、量子的揺らぎの効果があるので先のページの例の様に実時間発展してくれます)


そこで便利な手法が量子アニーリングです(この手法に特化した量子コンピュータもあります)。

簡単に言えば

H_initial = (1-dt/N)*h*Z_x(i); H_final = (dt/N)*J*Z(i)*Z(i+1)
H_total = H_initial + H_final

という様にハミルトニアンを用意します。

ここで dt は時間発展の幅で N*dt は最終的な発展した時間です(つまりdt=t/Nt は時間発展の最終時間)。

最初は横磁場のみのハミルトニアンで、これが時間発展することで少しづつ古典イジング模型のハミルトニアンになっていくという訳です。

この場合、始めに横磁場イジング模型での基底状態 qc.h(q) を持ってくると(断熱変化なので)時間発展中も基底状態となります。

つまり時間発展の最後において、古典イジング模型の基底状態が得られることになります。


ゲート方式での量子アニーリングの実際のコードを載せておくと、

# ライブラリをいろいろとインポート
import math
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import time
from tqdm import tqdm
from qiskit import QuantumRegister, ClassicalRegister, QuantumCircuit
from qiskit import BasicAer, execute

# shots:何回繰り返し観測するか. backend:シミュレータを設定
shots = 2000
backend = BasicAer.get_backend("qasm_simulator")

# 必要なイジングゲート
# Z方向の外場項
def unit_Z(qc, q0, h):
    qc.rz(2*h, q0)
# 相互作用項   
def unit_ZZ(qc, q0, q1, k):
    qc.cx(q0, q1)
    qc.rz(2*k, q1)
    qc.cx(q0, q1)
# X(横)方向の外場項
def unit_X(qc, q0, h):
    qc.rx(2*h, q0)

# スピンの期待値を計算するためにqubitをカウントする関数:例 |0101> -> 2
def digit_sum(n):
    digit_str = str(n)
    sum = 0
    for i in range(len(digit_str)):
        sum = sum + int(digit_str[i])
    return sum

# Ising model の quantum annealing
N = 4
# coupling constants for ZZ and coefficient of Z
J = 1.0; m = 0.1
# coefficient of X
h = 1.0
# step number
t  = 0.5; Nt = 100

x =[]; y=[]
# 時間発展を指定
for i in tqdm(range(Nt)):
    x.append(i)
# 回路を設定    
    q  = QuantumRegister(N, "q")
    c  = ClassicalRegister(N, "c")
    qc = QuantumCircuit(q, c)
# 初期状態を横磁場イジング模型の基底状態で準備
    qc.h(q)
# 時間発展を実行
    for j in range(i):
        for k in range(N):
            unit_X (qc, q[k], h*t*(1.0-j/Nt))
            unit_ZZ(qc, q[k], q[(k+1)%N], J*t*j/Nt)
            unit_Z (qc, q[k], m*t*j/Nt)        
# 測定
    for j in range(N):
        qc.measure(q[j], c[j])        
    job = execute(qc, backend=backend, shots=shots)
    result = job.result()
    counts = result.get_counts()
# スピン方向の期待値を計算
    M = 0
    for digit, count in counts.items():
        M = M +(abs(2*digit_sum(digit)/N-1))*count/shots
    y.append(M)
    time.sleep(0.1)
# 図をプロット
plt.plot(x,y, "-")
plt.show()

という感じになります。これを実行すると下の図が得られます。

横軸は時間、縦軸はスピンの方向の期待値です。

無事にスピンの方向が全て揃った「エネルギーの最も低い状態」が得られましたね(縦軸の ±1 がすべて揃った状態)。

ちなみに J の符号を変えるとスピン方向がランダムな方がエネルギーが下がるので 0 に漸近する結果になります。

z 方向の外場 m の大きさを J が負の場合でいろいろ変えてみるのも面白いので試してみましょう。


この手法を利用すると事で、いろいろな最適化問題を解くことも可能となります。

(解きたい問題をイジング模型で表現することが必要ですが)

例えば 記事 でアニーリングマシンを用いて経路の最適化を行っていることを、ゲート方式の量子コンピュータで計算できます。

まあゲート方式の場合、イジング模型に焼き直さずとも原理的に解けるわけなので、すこし冗長な感もありますね。