core interfaceを用いた量子イジング模型の数値シミュレーション

core interfaceを用いた量子イジング模型の数値シミュレーション#

このチュートリアルでは、OpenJijのcore interfaceを用いた、量子 (主に横磁場)効果の入ったイジング模型の数値シミュレーションを行います。まずはGraphを定義し、数値シミュレーションを行う系の $J_{ij}, h_i$ を定義します。今回は、古典イジング模型の場合にも扱った、ランダム相互作用およびランダムな縦磁場をもつ系のシミュレーションを行なっていきます。

import openjij.cxxjij.graph as G
#問題サイズを100とします。
N = 100

graph = G.Dense(N)

import numpy as np
mu, sigma = 0, 1

for i in range(N):
    for j in range(N):
        #Jijの値が大きくなりすぎてしまうので、全体の係数を1/Nしています。
        graph[i,j] = 0 if i == j else np.random.normal()/N

for i in range(N):
    graph[i] = np.random.normal()/N

グラフの設定方法は前章の古典イジング模型の場合と同じです。

システムの設定横磁場イジング模型#

今回はシステムに横磁場イジング模型

H = s \left(\sum_{i<j}J_{ij}\sigma_i^z \sigma_j^z + \sum_{i=1}^{N}h_i \sigma_i^z \right) - \sum_{i=1}\Gamma (1-s) \sigma_i^x

を用います。

$\Gamma$ は固定されたまま、 $\beta$ 、 $s$ を変化させて量子モンテカルロ法を行います。デフォルトでは鈴木・トロッター分解による量子モンテカルロ法が実装されています。

連続虚時間量子モンテカルロ法も用意してはいますが、現在試験的実装となっています。

まずはシステムを生成してみます。system.make_transverse_isingで生成できます。

import openjij.cxxjij.system as S

mysystem = S.make_transverse_ising(graph.gen_spin(), graph, 1.0, 4)

ここで、1つ目の引数にはスピン列を、2つ目にはグラフ、3つ目には $\Gamma$ の値、4つ目にはtrotterスライスの数を入力します。これで、全てのtrotterスライスが graph.gen_spin()で初期化されたシステムができます。

mysystem.trotter_spinsで全てのtrotterスピンを表示します。縦方向が空間方向、横方向がtrotter方向です。全てのtrotterスライスが同じスピンで初期化されていることがわかります。

print(mysystem.trotter_spins)

[[-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.]
 [ 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.]
 [ 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.]
 [ 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.]
 [-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.]
 [-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.]
 [-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.]
 [ 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.]
 [-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.]
 [-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.]
 [ 1.  1.  1.  1.]]

graph.gen_spin()の代わりに上の二重Listを入れて直接trotterスピンを初期化することができます。

アルゴリズムの実行 -Updater, Algorithm-#

Updater#

量子モンテカルロ法に対しては、現状

SingleSpinFlip (メトロポリス・ヘイスティング法によるスピン1つずつのアップデート)

が使用可能です。

Algorithm#

スケジュールリスト#

スケジュールリストは(パラメータ, モンテカルロステップ数)のリストで与えられ、横磁場イジングモデルに対しては(( $\beta$ , $s$ ), モンテカルロステップ数)で与えます。例として以下のように設定してみましょう。

schedule_list = [((10, 0.1), 10),((12, 0.3), 80),((10, 0.8), 30)]

この場合、逆温度 $\beta=10, s=0.1$ で10モンテカルロステップ、 $\beta=12, s=0.3$ で80ステップ、 $\beta=0.1, s=0.8$ で30ステップの計120モンテカルロステップを実行することを意味します。
アニーリングを実行するにあたっては、以下のようにutilityにあるmake_transeverse_field_schedule_listを使うとより便利です。

import openjij.cxxjij.utility as U
schedule_list = U.make_transverse_field_schedule_list(10, 20, 10)
print(schedule_list)

[((beta: 10.000000, s: 0.000000) mcs: 20), ((beta: 10.000000, s: 0.111111) mcs: 20), ((beta: 10.000000, s: 0.222222) mcs: 20), ((beta: 10.000000, s: 0.333333) mcs: 20), ((beta: 10.000000, s: 0.444444) mcs: 20), ((beta: 10.000000, s: 0.555556) mcs: 20), ((beta: 10.000000, s: 0.666667) mcs: 20), ((beta: 10.000000, s: 0.777778) mcs: 20), ((beta: 10.000000, s: 0.888889) mcs: 20), ((beta: 10.000000, s: 1.000000) mcs: 20)]

上の例では $\beta=10$ で固定しながら $s=0$ から $s=1$ まで、各パラメータで20モンテカルロステップ計算しながら10段階で $s$ を変えていく設定例です。計200モンテカルロステップの計算を行います。 $s$ の変化についてはMorita, Nishimori (2008)の手法を適用しています。

Algorithmの実行#

続いて、Algorithmを実行します。前章と全く同じように書けます。

import openjij.cxxjij.algorithm as A
A.Algorithm_SingleSpinFlip_run(mysystem, schedule_list)

前章と同じようにcallbackを使ってみましょう。横磁場イジング模型の場合は、システムとパラメータ (逆温度 $\beta$ 、 $s$ )を引数を持つ関数を作成すれば良いです。
例として、以下ではシステムのエネルギーの値を記録するcallbackを作っています。

energies = []

def callback_log_energy(system, t):
    #graphは以前にGraphモジュールにて定義したオブジェクトです
    #各trotterスライスの平均値から、古典スピンの0、1を決めます。
    classical_spin = [-1 if np.mean(s)<0 else 1 for s in system.trotter_spins[:-1]] #最後のスピンは補助スピンのため、除く
    energies.append(graph.calc_energy(classical_spin))

このcallbackを用いて同じAlgorithmを実行します。

#スケジュールをもっと長く取ります (計20000モンテカルロステップ)
schedule_list = U.make_transverse_field_schedule_list(10, 200, 100)
A.Algorithm_SingleSpinFlip_run(mysystem, schedule_list, callback_log_energy)

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[9], line 3
      1 #スケジュールをもっと長く取ります (計20000モンテカルロステップ)
      2 schedule_list = U.make_transverse_field_schedule_list(10, 200, 100)
----> 3 A.Algorithm_SingleSpinFlip_run(mysystem, schedule_list, callback_log_energy)

Cell In[8], line 3, in callback_log_energy(system, t)
      1 energies = []
----> 3 def callback_log_energy(system, t):
      4     #graphは以前にGraphモジュールにて定義したオブジェクトです
      5     #各trotterスライスの平均値から、古典スピンの0、1を決めます。
      6     classical_spin = [-1 if np.mean(s)<0 else 1 for s in system.trotter_spins[:-1]] #最後のスピンは補助スピンのため、除く
      7     energies.append(graph.calc_energy(classical_spin))

KeyboardInterrupt: 

記録したシステムのエネルギーを、横軸をモンテカルロステップ、縦軸をエネルギーでプロットすると次のようになります。

!pip install matplotlib

Requirement already satisfied: matplotlib in /opt/conda/lib/python3.9/site-packages (3.5.2)
Requirement already satisfied: cycler>=0.10 in /opt/conda/lib/python3.9/site-packages (from matplotlib) (0.11.0)
Requirement already satisfied: packaging>=20.0 in /opt/conda/lib/python3.9/site-packages (from matplotlib) (21.3)
Requirement already satisfied: python-dateutil>=2.7 in /opt/conda/lib/python3.9/site-packages (from matplotlib) (2.8.2)
Requirement already satisfied: pyparsing>=2.2.1 in /opt/conda/lib/python3.9/site-packages (from matplotlib) (2.4.7)
Requirement already satisfied: fonttools>=4.22.0 in /opt/conda/lib/python3.9/site-packages (from matplotlib) (4.37.1)
Requirement already satisfied: pillow>=6.2.0 in /opt/conda/lib/python3.9/site-packages (from matplotlib) (9.2.0)
Requirement already satisfied: kiwisolver>=1.0.1 in /opt/conda/lib/python3.9/site-packages (from matplotlib) (1.4.4)
Requirement already satisfied: numpy>=1.17 in /opt/conda/lib/python3.9/site-packages (from matplotlib) (1.23.2)
Requirement already satisfied: six>=1.5 in /opt/conda/lib/python3.9/site-packages (from python-dateutil>=2.7->matplotlib) (1.16.0)

%matplotlib inline
import matplotlib.pyplot as plt
plt.plot(range(len(energies)), energies)
plt.xlabel('Monte Carlo step')
plt.ylabel('energy')
plt.show()

../../../_images/87bfa7fbe46a61bebf5ef5d400c2bbe22d753c15e927498fb0f784e373e1f303.png

結果の取得 -Result-#

result.get_solutionsで計算結果である古典スピンを取得します。この関数は最適化問題を解く観点にフォーカスを当てているため、trotterスライスの中でもっともエネルギーが低いスピン列を返します。

import openjij.cxxjij.result as R
print(R.get_solution(mysystem))

[-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, 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, 1, 1, -1, 1, -1, 1, -1, -1, -1, -1, -1, -1, 1, 1, -1, -1, 1, -1, 1, -1]