こんにちは。Juliaです。今日は横磁場イジングモデルの量子アニーリングを量子モンテカルロ法で実装して速度を検証してみます。本格的に使いたいというよりはJuliaの練習で使いますが、思ったよりも使い勝手がいいので、普通にSDKにできそうな勢いです。
参考
こちらを参考にさせていただきました。
https://gist.github.com/shinaoka/ada5a6db8d8f3fe4ff81a0a90370b6dc
横磁場イジング量子アニーリング
ある、ハミルトニアンの定式化を解くのに量子アニーリングを利用します。今回の定式化は+1と-1のスピンの値を使って、ハミルトニアンを作りますが、上三角行列の形で実装をしました。
量子モンテカルロ
量子モンテカルロ法を使うと、２次元の量子系の計算を３次元の古典系に落とすことができます。鈴木トロッター展開を使って、定式化をしました。今回はトロッタ数を指定して実行できるようにしました。
コード
早速みてみましょう。
@show VERSION
import Pkg; Pkg.add("BenchmarkTools")
using BenchmarkTools, Random, LinearAlgebra
少し慣れてきました。ベンチマークツールをパッケージとして導入し、RandomとLinearAlgebraを読み込みました。
VERSION = v"1.5.1"
次に今回の定式化ですが、とりあえず任意のサイズの問題をランダムで作ります。eleに量子ビット数を入れると、ランダム値で上三角行列を作って返してくれます。これを今回は定式化として使います。もちろん任意の問題を入れてもOKです。
function rands(ele)
return triu(rand(rng,Float64,(ele,ele)))
end
次にエネルギーの値の計算ですが、
function Energy(qq,jj)
return qq'*jj*qq
end
ハミルトニアンと量子ビットを使ってこのように計算できます。早速メインの計算を見てみましょう。
function run_opt!(q, J, T::Float64, r::Float64, Gs::Float64, Gf::Float64, N, tro, ite, rng)
#初期の磁場の強さをここで設定
G = Gs
#磁場が弱まるまで計算する
while G > Gf
for i=1:ite
#フリップさせる量子ビット数の指定
x = rand(rng, 1:N)
#トロッター番号を指定
y = rand(rng, 1:tro)
dE = 0
#差だけを計算
for j=1:N
if j == x
dE += -2*q[y,x]*J[x,x] / tro
else
dE += -2*q[y,j]*q[y,x]*J[j,x] / tro
end
end
#トロッタ間の計算
dE += q[y,x]*(q[(y+tro-1)%tro+1,x]+q[(y+1)%tro+1,x])*log(1/tanh(G/T/tro))*T
#メトロポリス法で判定
if exp(-dE/T) > rand(rng)
q[y,x] *= -1
end
end
#一定の回数計算したら磁場を弱める
G *= r
end
return G
end
次に設定値です。
#量子ビット数
N = 5
#温度は固定にしました
T = 0.02
#初期の磁場
Gs = 10.0
#ターゲットの磁場
Gf = 0.02
#トロッタ数
tro = 8
#イテレーションの回数
ite = 10000
#磁場の減衰の係数
r = 0.95
#乱数生成はシードを指定してメルセンヌツイスターで
rng = MersenneTwister(4649)
#量子ビットの値を-1か+1で初期化それをトロッタ数だけ実行
q = rand(rng, (-1,1), tro, N)
#ランダムにハミルトニアンを作る
J = rands(N)
そうすると、このようにハミルトニアンができます。
5×5 Array{Float64,2}:
0.673513 0.673642 0.969669 0.647173 0.2787540.0
0.396891 0.83374 0.778207 0.3100230.0 0.0
0.377737 0.31789 0.8472050.0 0.0 0.0
0.278436 0.05608120.0 0.0 0.0 0.0 0.647878
では、早速実行してみましょう。
Gf_ = run_opt!(q, J, T, r, Gs, Gf, N, tro, ite, rng)
println()
println("Gfinal = ",Gf_)
println(Energy(q[1,:],J))
show(stdout, "text/plain", q)
これは、
Gfinal = 0.019154898067750878
1.1515795449183275
8×5 Array{Int64,2}:
-1 1 -1 -1 1
-1 1 -1 -1 1
-1 1 -1 -1 1
-1 1 -1 -1 1
-1 1 -1 -1 1
-1 1 -1 -1 1
-1 1 -1 -1 1
-1 1 -1 -1 1
こうなりました。各トロッタがきちんと値が揃っていて、エネルギー値が出ています。最後にベンチマークを見てみましょう。
q = rand(rng, (-1,1), tro, N)
rng = MersenneTwister(4649)
@benchmark run_opt!(q, J, T, r, Gs, Gf, N, tro, ite, rng)
こちらを実行したところ、
BenchmarkTools.Trial:
memory estimate: 16 bytes
allocs estimate: 1
minimum time: 188.530 ms (0.00% GC)
median time: 191.590 ms (0.00% GC)
mean time: 193.643 ms (0.00% GC)
maximum time: 214.312 ms (0.00% GC)
samples: 26
evals/sample: 1
いかがでしょうか。200ms未満で計算ができています。これは体感でもかなり早いと思います。
最後にPythonと比較
最後に同様のコードをpythonで作ったものを比較してみます。pythonコードは、
import numpy as np
import time
np.random.seed(10)
def rands(ele):
return np.triu(np.random.rand(ele,ele))
def Energy(qq,jj):
return qq @ jj @ qq
Tf = 0.02
Gs = 10
Gf = 0.02
tro = 8
ite = 10000
R = 0.95
J = rands(5)
def sqa():
G = Gs
N = len(J)
q = np.random.choice([-1,1], tro*N).reshape(tro, N)
while G > Gf:
for _ in range(ite):
x = np.random.randint(N)
y = np.random.randint(tro)
dE = 0
for j in range(N):
if j == x:
dE += -2*q[y][x]*J[x][x] / tro
else:
dE += -2*q[y][j]*q[y][x]*J[j][x] / tro
dE += q[y][x]*(q[(tro+y-1)%tro][x]+q[(y+1)%tro][x])*np.log(1/np.tanh(G/Tf/tro))*Tf
if np.exp(-dE/Tf) > np.random.rand():
q[y][x] *= -1
G *= R
print(Energy(q[0], J))
print(q)
start = time.time()
sqa()
print(time.time() - start)
ちょっと問題の種類が違ってしまいましたが、何回かやってみて、やはり45秒前後かかりました。200倍近く早くなっています。
0.8764328810175842
[[-1 -1 1 -1 1]
[-1 -1 1 -1 1]
[-1 -1 1 -1 1]
[-1 -1 1 -1 1]
[-1 -1 1 -1 1]
[-1 -1 1 -1 1]
[-1 -1 1 -1 1]
[-1 -1 1 -1 1]]
46.937604665756226
パラメータの設定などはほぼ同じで、書き方も同じですが、Juliaもかなり早く実行ができました。体感で速度がわかるので使っていて楽しいですね。以上です。今後は汎用計算や量子古典計算、NISQ計算にも利用してみたいと思います。
では、良いJuliaを！