common.title

Juliaで横磁場イジング量子アニーリングを量子モンテカルロで実装

Yuichiro Minato 2 years ago

#量子アニーリング #Julia

2

こんにちは。Juliaです。今日は横磁場イジングモデルの量子アニーリングを量子モンテカルロ法で実装して速度を検証してみます。本格的に使いたいというよりはJuliaの練習で使いますが、思ったよりも使い勝手がいいので、普通にSDKにできそうな勢いです。


参考

こちらを参考にさせていただきました。

https://gist.github.com/shinaoka/ada5a6db8d8f3fe4ff81a0a90370b6dc


横磁場イジング量子アニーリング

ある、ハミルトニアンの定式化を解くのに量子アニーリングを利用します。今回の定式化は+1と-1のスピンの値を使って、ハミルトニアンを作りますが、上三角行列の形で実装をしました。


量子モンテカルロ

量子モンテカルロ法を使うと、2次元の量子系の計算を3次元の古典系に落とすことができます。鈴木トロッター展開を使って、定式化をしました。今回はトロッタ数を指定して実行できるようにしました。


コード

早速みてみましょう。


@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
  <span class="hljs-comment">#トロッタ間の計算       </span>

      dE += q[y,x](q[(y+tro-1)%tro+1,x]+q[(y+1)%tro+1,x])log(1/tanh(G/T/tro))*T

  <span class="hljs-comment">#メトロポリス法で判定       </span>

      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 += -2q[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を!

Yuichiro Minato

@yuichiro_minato

blueqat CEO

© 2022, blueqat Inc. All rights reserved