こんにちは。Juliaです。今日は横磁場イジングモデルの量子アニーリングを量子モンテカルロ法で実装して速度を検証してみます。本格的に使いたいというよりはJuliaの練習で使いますが、思ったよりも使い勝手がいいので、普通にSDKにできそうな勢いです。 <h2>参考</h2>こちらを参考にさせていただきました。<a href="https://gist.github.com/shinaoka/ada5a6db8d8f3fe4ff81a0a90370b6dc" rel="noopener noreferrer" target="_blank">https://gist.github.com/shinaoka/ada5a6db8d8f3fe4ff81a0a90370b6dc</a> <h2>横磁場イジング量子アニーリング</h2>ある、ハミルトニアンの定式化を解くのに量子アニーリングを利用します。今回の定式化は+1と-1のスピンの値を使って、ハミルトニアンを作りますが、上三角行列の形で実装をしました。 <h2>量子モンテカルロ</h2>量子モンテカルロ法を使うと、２次元の量子系の計算を３次元の古典系に落とすことができます。鈴木トロッター展開を使って、定式化をしました。今回はトロッタ数を指定して実行できるようにしました。 <h2>コード</h2>早速みてみましょう。 <pre class="ql-syntax" spellcheck="false">@show VERSION
import Pkg; Pkg.add(&#34;BenchmarkTools&#34;)
using BenchmarkTools, Random, LinearAlgebra
</pre> 少し慣れてきました。ベンチマークツールをパッケージとして導入し、RandomとLinearAlgebraを読み込みました。 <pre class="ql-syntax" spellcheck="false">VERSION = v&#34;1.5.1&#34;
</pre> 次に今回の定式化ですが、とりあえず任意のサイズの問題をランダムで作ります。eleに量子ビット数を入れると、ランダム値で上三角行列を作って返してくれます。これを今回は定式化として使います。もちろん任意の問題を入れてもOKです。 <pre class="ql-syntax" spellcheck="false">function rands(ele)
  return triu(rand(rng,Float64,(ele,ele)))
end
</pre> 次にエネルギーの値の計算ですが、 <pre class="ql-syntax" spellcheck="false">function Energy(qq,jj)
  return qq&#39;*jj*qq
end
</pre> ハミルトニアンと量子ビットを使ってこのように計算できます。早速メインの計算を見てみましょう。 <pre class="ql-syntax" spellcheck="false">function run_opt!(q, J, T::Float64, r::Float64, Gs::Float64, Gf::Float64, N, tro, ite, rng)
 #初期の磁場の強さをここで設定
  G = Gs
 #磁場が弱まるまで計算する
  while G &gt; 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) &gt; rand(rng)
        q[y,x] *= -1
      end
    end
 #一定の回数計算したら磁場を弱める
    G *= r
  end
  return G
end
</pre> 次に設定値です。 <pre class="ql-syntax" spellcheck="false">#量子ビット数
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)
</pre> そうすると、このようにハミルトニアンができます。 <pre class="ql-syntax" spellcheck="false">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
</pre> では、早速実行してみましょう。 <pre class="ql-syntax" spellcheck="false">Gf_ = run_opt!(q, J, T, r, Gs, Gf, N, tro, ite, rng)

println()
println(&#34;Gfinal = &#34;,Gf_)
println(Energy(q[1,:],J))
show(stdout, &#34;text/plain&#34;, q) 
</pre> これは、 <pre class="ql-syntax" spellcheck="false">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
</pre> こうなりました。各トロッタがきちんと値が揃っていて、エネルギー値が出ています。最後にベンチマークを見てみましょう。 <pre class="ql-syntax" spellcheck="false">q = rand(rng, (-1,1), tro, N)
rng = MersenneTwister(4649)
@benchmark run_opt!(q, J, T, r, Gs, Gf, N, tro, ite, rng)
</pre> こちらを実行したところ、 <pre class="ql-syntax" spellcheck="false">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
</pre> いかがでしょうか。200ms未満で計算ができています。これは体感でもかなり早いと思います。 <h2>最後にPythonと比較</h2>最後に同様のコードをpythonで作ったものを比較してみます。pythonコードは、 <pre class="ql-syntax" spellcheck="false">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 &gt; 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) &gt; np.random.rand():
        q[y][x] *= -1
    G *= R
  print(Energy(q[0], J))
  print(q)

start = time.time()
sqa()
print(time.time() - start)
</pre> ちょっと問題の種類が違ってしまいましたが、何回かやってみて、やはり45秒前後かかりました。200倍近く早くなっています。 <pre class="ql-syntax" spellcheck="false">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
</pre> パラメータの設定などはほぼ同じで、書き方も同じですが、Juliaもかなり早く実行ができました。体感で速度がわかるので使っていて楽しいですね。以上です。今後は汎用計算や量子古典計算、NISQ計算にも利用してみたいと思います。 では、良いJuliaを！

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

Yuichiro Minato