A3 - 古典イジング模型

Open in Colab

この章ではOpenJijのcore interface (core python interface)の使い方を説明し、簡単な計算のデモンストレーションを行います。

core interfaceは前回までのチュートリアルよりも下部のレイヤーのAPIです。対象読者としては前回までのOpenJijチュートリアルを一通り終えて、イジングモデルやモンテカルロ法などの用語を知っている方を想定します。次のような目的を持った読者を想定しています。

  • 最適化問題だけでなくサンプリングや研究用途などより専門的な用途にOpenJijを用いたい

  • アニーリングスケジュールの設定や使用するアルゴリズム等を直接触りたい

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

graph = G.Dense(N)
#sparseの場合は以下を使用します。
#graph = G.Sparse(N)

\(J_{ij}, h_i\)を設定します。今回は平均0、標準偏差1のGauss分布から生成される値を設定します。

[ ]:
!pip install numpy #乱数生成にnumpyを使います。
[2]:
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

縦磁場に関しては、graph[i]でも、graph[i,i]でもどちらでもアクセスできます。また、イジングモデルの定義上、\(J_{ij}\)\(J_{ji}\)は自動で同じ値となります。試しに以下のように出力を行なってみましょう。

[3]:
graph[20] = 0.5
print(graph[20,20])
print(graph[20])
graph[12,34] = -0.6
print(graph[12,34])
print(graph[34,12])
0.5
0.5
-0.6
-0.6

システムの設定 - System -

続いて計算を行うためのシステムを定義します。ここでは古典イジング模型か横磁場イジング模型か、また別の模型にするか等を選べます。

まずは古典イジング模型のシステムを作成してみます。system.make_classical_isingで作成できます。

[4]:
import cxxjij.system as S

mysystem = S.make_classical_ising(graph.gen_spin(), graph)

ここで、1つ目の引数にはランダムに生成したスピン、2つめにはGraphそのものを代入します。これにより初期のスピン配位がgraph.gen_spin()となる古典イジング模型のシステムの作成ができます。

システムに直接アクセスして、値を読むことも可能です。

[5]:
print(mysystem.spin)
[ 1. -1.  1. -1. -1. -1.  1. -1.  1.  1. -1.  1. -1.  1.  1. -1.  1. -1.
  1.  1.  1.  1. -1.  1. -1. -1.  1.  1.  1.  1. -1.  1. -1.  1. -1.  1.
 -1. -1.  1. -1. -1.  1. -1.  1.  1. -1. -1.  1. -1.  1.  1.  1. -1. -1.
  1.  1.  1. -1. -1. -1.  1.  1.  1.  1. -1. -1. -1. -1.  1. -1.  1. -1.
  1. -1.  1. -1.  1. -1. -1.  1.  1. -1.  1. -1.  1.  1. -1.  1. -1.  1.
 -1. -1.  1. -1.  1.  1. -1.  1. -1. -1.  1.]

古典イジング模型以外にも様々なシステムが用意されており、これを用途別に使うことが出来ます。用いるSystemによって初期化の方法は多少異なります。これは後ほど少しずつ紹介していきます。

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

Systemを定義した後はUpdaterを選択し、Algorithmを実行していきます。

Updater

Systemに対して使用できるUpdaterは決められており、古典イジング模型に対するUpdaterは主に

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

  • SwendsenWang (SwendsenWang法によるクラスターアップデート)

の2つが用意されています。 Algorithmを実行するには、スケジュールリストが必要になります。まずスケジュールリストを作成するところから始めましょう。

Algorithm

スケジュールリスト

スケジュールリストは(パラメータ, モンテカルロステップ数)のリストで与えられるものです。パラメータ部分に入力する値はSystemによって異なります。例えば、古典イジング模型ならばパラメータとして温度の逆数である逆温度\(\beta\)を設定します。例として以下のように設定してみましょう。

[6]:
schedule_list = [(0.01, 10),(10, 80),(0.1, 30)]
この場合、逆温度\(\beta=0.01\)で10モンテカルロステップ、\(\beta=10\)で80ステップ、\(\beta=0.1\)で30ステップの計120モンテカルロステップを実行することを意味します。
アニーリングを実行するにあたっては、逆温度は等比級数で増加させていくことが多いため、以下のようにutilityにあるmake_classical_schedule_listを使うとより便利です。
[7]:
import cxxjij.utility as U
schedule_list = U.make_classical_schedule_list(0.1, 50, 20, 10)
print(schedule_list)
[((beta: 0.100000) mcs: 20), ((beta: 0.199474) mcs: 20), ((beta: 0.397897) mcs: 20), ((beta: 0.793701) mcs: 20), ((beta: 1.583223) mcs: 20), ((beta: 3.158114) mcs: 20), ((beta: 6.299605) mcs: 20), ((beta: 12.566053) mcs: 20), ((beta: 25.065966) mcs: 20), ((beta: 50.000000) mcs: 20)]

上の例では\(\beta=0.1\)から\(\beta=50\)まで、各温度で20モンテカルロステップ計算しながら10段階で温度を変えていく設定例です。計200モンテカルロステップの計算を行います。

Algorithmの実行

続いて、Algorithmを実行します。Algorithm_[Updater]_runのように書くことで、指定したUpdaterで計算を行います。次の例ではSingleSpinFlipを実行しています。

[8]:
import cxxjij.algorithm as A
A.Algorithm_SingleSpinFlip_run(mysystem, schedule_list)

一瞬で処理が終わりましたが、この間に計200モンテカルロステップの計算が行われています。 > A.Algorithm_SingleSpinFlip_run(mysystem, seed, schedule_list)とすることで、seedを固定したまま計算を行うことができます。結果に再現性をもたせたい際に使うことができます。

callbackを使用することで、Algorithmの実行中に1モンテカルロステップごとのシステムを取得することができます。古典イジング模型の場合は、システムとパラメータ (逆温度)を引数を持つ関数を作成すれば良いです。
例として、以下ではシステムのエネルギーの値を記録するcallbackを作っています。
[9]:
energies = []

def callback_log_energy(system, beta):
    #graphは以前にGraphモジュールにて定義したオブジェクトです
    energies.append(graph.calc_energy(system.spin))

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

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

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

[ ]:
!pip install matplotlib
[11]:
%matplotlib inline
import matplotlib.pyplot as plt
plt.plot(range(len(energies)), energies)
plt.xlabel('Monte Carlo step')
plt.ylabel('energy')
plt.show()
../_images/ja_A003-LargeScaleMC_26_0.png

アニーリングが進むに連れ徐々にエネルギーが低くなっているのが分かります。Algorithmの動作中にシステムの様子を知りたい時に有用です。

結果の取得 -Result-

result.get_solutionで計算結果であるスピン列を取得できます。古典イジング模型の場合は直接mysystem.spinを参照することで、スピン列を取得も可能です。しかし、result.get_solutionはそれ以外のシステムについてもスピン列を取得できる便利なメソッドです。

[12]:
import 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]

このスピン列がアニーリングによって得られた答えです。ハミルトニアンの基底状態 (に近い状態)であることが期待されます。

C++ core interface

多少の違いはありますが、C++ core interfaceでも上記とほぼ同じことが可能です。以下の点に留意しましょう。

  • seed値を入れる引数には、乱数生成器 (C++11 random)を代入する必要があります。

  • Graphクラスで、\(J_{ij}, h_i\)へのアクセス方法が多少異なります。

今までの内容をC++ core interfaceで記述すると次のようになります。

#include <graph/all.hpp>
#include <system/all.hpp>
#include <updater/all.hpp>
#include <algorithm/all.hpp>
#include <result/all.hpp>
#include <utility/schedule_list.hpp>
#include <utility/random.hpp>
#include <random>

#include <iostream>

using namespace openjij;

int main(void){

    //generate dense graph with size N=100
    constexpr std::size_t N = 100;
    auto dense = graph::Dense<double>(N);

    //generate random engine
    auto rand_engine = std::mt19937(0x1234);
    //of course you can specify another random generator that is compatible with C++ random generator, say utility::Xorshift,
    //auto rand_engine = utility::Xorshift(0x1234);

    //Gaussian distribution
    auto gauss = std::normal_distribution<>{0, 1};

    //set interactions
    for(std::size_t i=0; i<N; i++){
        for(std::size_t j=0; j<N; j++){
            dense.J(i, j) = (i == j) ? 0 : gauss(rand_engine)/N;
        }
    }

    //set local fields
    for(std::size_t i=0; i<N; i++){
        dense.h(i) = gauss(rand_engine);
    }

    //create classical Ising system
    auto system = system::make_classical_ising(dense.gen_spin(rand_engine), dense);

    //generate schedule list
    //from beta=0.1 to beta=50, 20 Monte Carlo step for each tempearature, 10 samples
    auto schedule_list = utility::make_classical_schedule_list(0.1, 50, 20, 10);

    //or in the case of classical system, you can specify arbitrary schedules by using a std::vector of utility::Schedule<system::classical_system>
    //std::vector<utility::Schedule<system::classical_system>> schedule_list;
    //utility::Schedule<system::classical_system> schedule;
    //schedule.updater_parameter = {0.01};
    //schedule.one_mc_step = 10; //number of Monte Carlo step per temperature
    //schedule_list.push_back(schedule);
    //
    //schedule.updater_parameter = {10};
    //schedule.one_mc_step = 80; //number of Monte Carlo step per temperature
    //schedule_list.push_back(schedule);
    //
    //schedule.updater_parameter = {0.1};
    //schedule.one_mc_step = 30; //number of Monte Carlo step per temperature
    //schedule_list.push_back(schedule); //schedule_list -> [(0.01, 10), (10, 80), (0.1, 30)]


    //do annealing (updater: SingleSpinFlip)
    algorithm::Algorithm<updater::SingleSpinFlip>::run(system, rand_engine, schedule_list);

    //show spins
    std::cout << "The result spins are [";
    for(auto&& elem : result::get_solution(system)){
        std::cout << elem << " ";
    }

    std::cout << "]" << std::endl;
}