本日は

Twitterのハッシュタグ Julia言語で話題となっている（というかしている）Ising ModelのJuliaのコードが速いというおはなしです．

Ising Model

Ising Model の簡単なイントロは Ising Modelを平易に解説してみるを読むといいでしょう．ここでは，”スピン”と呼ばれるプラス１またはマイナス１をとるものから構成される二次元領域の変化の記述を計算しているものだと思ってください．

Juliaでの実装

Juliaでの実装が Gen Kurokiさんのノートブック
にありますので抜粋します.リンク先にはコードの解説も書いてありますので参考にしてください．

#Julia言語最近では2次元イジング模型のメトロポリス法が「Hello world」扱いらしいので私が書いたテストコードも貼り付けておきます。https://t.co/AmpKeZg2eq
2D Ising Metropolis (Gen Kuroki, 2017-01-09)

10^9点の更新に8秒。最適化はもっとできます。自由な複製、改変、配布可。
— 黒木玄 Gen Kuroki (@genkuroki) 2018年1月16日

ising.jl


# 2D Ising Metropolis (Gen Kuroki, 2017-01-09)

const β_crit = log(1+sqrt(2))/2

rand_ising2d(n) = rand([-1, 1], n, n)

function ising2d_sum_of_adjacent_spins(s, m, n, i, j)
    s[ifelse(i+1 ≤ m, i+1, 1), j] + s[ifelse(i-1 ≥ 1, i-1, m), j] +
    s[i, ifelse(j+1 ≤ n, j+1, 1)] + s[i, ifelse(j-1 ≥ 1, j-1, n)]
end

function ising2d_sweep!(s, β, niters)
    m, n = size(s)
    prob = [exp(-2*β*k) for k in -4:4]
    for iter in 1:div(niters, m*n)
        for j in 1:n, i in 1:m
            s₁ = s[i,j]
            k = s₁ * ising2d_sum_of_adjacent_spins(s, m, n, i, j)
            s[i,j] = ifelse(rand() < prob[k+5], -s₁, s₁)
        end
    end
end

s = rand_ising2d(100)
@time ising2d_sweep!(s, β_crit, 10^9)
s

このコードを実行すると大体１０秒ほどかかります．この１０秒がどれぐらいの速さを表しているのかPythonのコードで比較してみます．

Python＋Numbaでの実装

JuliaがわからなくてもPythonならわかる方はいらっしゃると思いますのでJuliaの実装を移植したものを下に示します．PureなPythonはすごく遅いので，NumbaでJITコンパイルできるようにしています．

ising.py

# Reference: 
from math import log, sqrt, exp
from random import choice, random
from copy import deepcopy
import time
from numba import jit
import numpy as np


@jit('i8(i8[:,:],i8,i8,i8,i8)')
def ising2d_sum_of_adjacent_spins(s, m, n, i, j):
    i_bottom = i+1 if i+1 < m else 0
    i_top = i-1 if i-1 >= 0 else m-1
    j_right = j+1 if j+1 < n else 0
    j_left = j-1 if j-1 >= 0 else n-1
    return s[i_bottom, j]+s[i_top, j]+s[i, j_right]+s[i, j_left]


@jit('i8[:,:](i8[:,:],f8,i8)')
def ising2d_sweep(s, beta, niters):
    m, n = s.shape
    prob = np.exp(-2*beta*np.array([-4, -3, -2, -1, 0, 1, 2, 3, 4]))
    for _ in range(int(niters/(m*n))):
        for i in range(m):
            for j in range(n):
                s1 = s[i, j]
                k = s1*ising2d_sum_of_adjacent_spins(s, m, n, i, j)
                s[i, j] = -s1 if np.random.random() < prob[k+4] else s1
                # s[i, j] = -s1 if random() < prob[k+4] else s1 <- very slow
    return s

n = 100
beta_critical = log(1+sqrt(2))/2

def main():
    rand_ising2d = np.array([[choice([-1, 1])
                              for j in range(n)] for i in range(n)]).astype(np.int64)
    s_begin = deepcopy(rand_ising2d)
    begin = time.time()
    s_end = ising2d_sweep(rand_ising2d, beta_critical, 1e9)
    end = time.time()
    print("Elapsed=", end-begin)

if __name__ == '__main__':
    main()

マシーンにもよりますが１５～１７秒ほどかかると思います．この結果を見てもJuliaが速いというがわかると思います．

Numbaは遅いのか（いいえ．）？

ここで，s[i, j] = -s1 if np.random.random() < prob[k+4] else s1の箇所を
s[i, j] = -s1 if 1 < prob[k+4] else s1 と置き換えてみてください．そうすると（同様に変更した）Juliaのコードの実行速度と同等のスピードを得ます．

　どこがボトルネックになっているかというと乱数生成の部分です．

Juliaが速い理由の一つとして，Juliaが持っている乱数生成のスピードがNumpyがもっているそれに勝っているからというのが挙げられます．

ということで，Ising Modelの実装には乱数生成アルゴリズムの選択が大事になりそうです．

Julia は dSFMT だそうです．
Numpyは (https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.random.set_state.html) などを見ると Mersenne Twister っぽいですね．
https://heavywatal.github.io/cxx/random.html のサイトも参考になりそう．

他の言語ではどうか

気になりますね．Juliaの実装を移植する形でソースを書いてみました．全部紹介すると記事自体が長くなりそうなので私のGistを参照してください．

Cython

Numbaときたら次はCythonでも試したくなるでしょう．

https://gist.github.com/terasakisatoshi/fc8496f7fb497e63ac0c276944c68d5d
libc.stdlib cimport rand を使って最適化を行っていますがNumbaには及ばず．

Go

Webアプリケーションを作りたかったモチベーションとCの代用として勉強中

https://gist.github.com/terasakisatoshi/22abc51cf6b2e156f632765f2dcecd6e
私の実装では根本的に遅いので誰かトライしてください・・・．

Nim

https://gist.github.com/terasakisatoshi/9c110ce080069305d85c9d4e093656eb
パフォーマンスに定評のあるNimですが，Nimでの random はXorshift128を使っているようですね．速度はNumba並みです．　個人的にはもっと速くなって欲しかった．

Rust

周りがRust推しなのでC++の代用として勉強中．

https://gist.github.com/terasakisatoshi/d54f0758e84592fc88573bedf5573cb6 Rustで乱数生成するときの処理速度の違いを参考にXorShiftRngを使うことにしました．結果は速いのですがJuliaに及ばす．

結果

https://gist.github.com/terasakisatoshi/39398cc3127e7a49e736bb58cf7ff3ea
にも記録しましたが，下記のとおりです：

niter=1e9の場合の計測時間表．

PC環境（Macbook12inch2015）

julia>versioninfo()
Platform Info:
OS: macOS (x86_64-apple-darwin14.5.0)
CPU: Intel(R) Core(TM) M-5Y51 CPU @ 1.10GHz
WORD_SIZE: 64
BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
LAPACK: libopenblas64_
LIBM: libopenlibm
LLVM: libLLVM-3.9.1 (ORCJIT, broadwell)

時間計測結果

単位は秒．小数第3位で四捨五入
Julia 10.58
Rust 12.05
Nim 17.37
Python(with Numba) 17.50
Cython 56.01
Go 1m23.67

わかったこと

地球は青かった．そして，Juliaは速かった．

最後に

皆さんも得意な言語でトライしてみてください．
タイトルは速度の云々に絞ってしまいましたが，　

実際にトライする際に上で抜粋したJuliaと比べて

導入しやすい（環境構築が容易か）
書きやすいか（”実装”の速さ，ライブラリを導入する場合は使うための学習コスト）
読みやすいか（可読性が良い）
乱数生成のアルゴリズムが選べるのであれば，質（生成スピードが優れている，　乱数の周期が長いかどうか）はどうか（IsingModeｌの実装においてはココが大事になってきます．）　なども吟味してみると良いかなと思いました．

https://t.co/AmpKeZg2eq

にある2Dイジング模型のテストコードは、型指定を行ったり、最適化のおまじない(inboundsの類)を使ったり、アルゴリズムでさらに工夫したりすることを、意図的に一切しないようにしています。

特に擬似乱数を32ビット浮動小数で生成するのは御法度。 #Julia言語
— 黒木玄 Gen Kuroki (@genkuroki) 2018年1月23日

速度でJuliaに勝つのが難しかった（私では無理だった）話．

本日は