Pythonでのモンテカルロ法の実装を手抜きで速くする．

本日は

モンテカルロ法を実装したコードを早くしてみましょう．

import random
NUM=100000000
def monte():
    counter = 0
    for i in range(NUM):
        x = random.random()
        y = random.random()
        if x*x+y*y < 1.0:
            counter += 1
    pi = 4.0*counter/NUM
    print(pi)


def main():
    monte()

if __name__ == '__main__':
    main()

$ time python monte.py
real    1m37.223s
user    1m35.556s
sys 0m0.636s

自分の環境ですと100秒ほどかかるのでやってられないです．

Cython　を使いましょう．

明らかに型はこれだろ！ってのをコンピュータに教えます．そして型推定する時間を省くことにします．

#monte.pyx
import random

cdef int NUM = 100000000

cdef cmonte():
    cdef :
        int counter = 0
        int i=0
        double x
        double y
    for i in range(NUM):
        x = random.random()
        y = random.random()
        if x*x + y*y < 1.0:
            counter += 1

    cdef double pi = 4.0*counter/NUM
    return pi

def monte():
    pi=cmonte()
    print(pi)

pyximport　を使って通常のPythonのモジュールのようにインポートできるようにします．

#main.py
import pyximport
pyximport.install()
import monte
import time

print("start")
s=time.time()
monte.monte()
e=time.time()
print("elapsed",e-s)

$ time python main.py
start
3.14175328
elapsed 17.078665018081665

real    0m17.540s
user    0m17.211s
sys     0m0.158s

型つけのヒントがあるかないかで全然違ってきます．　
補足記事にてもっと速くしました↓
Cython手抜き実装のモンテカルロ法のスピードを速くする．

Numpy を使いましょう．

for 使ってるのがいけないんや．．． Numpy の配列を一気に生成して
各成分の二乗和を求めて　1．0 未満の部分だけをカウントすればOKという方式

import numpy as np
NUM=100000000

def no_for():
    xs = np.random.rand(NUM)
    ys = np.random.rand(NUM)
    square = xs*xs+ys*ys
    counter = np.sum(square < 1.0)
    pi = 4.0*counter/NUM
    print(pi)

def main():
    no_for()

if __name__ == '__main__':
    main()

$ time python monte.py
real    0m10.838s
user    0m5.213s
sys     0m3.200s

かなりましになりました． 書き方の違いはあれど結果だけ見るとCythonの呼び出しよりも速くできますね．
forを書かずにできるような方法があればこちらの方が良いということでしょうか．

Numba　の力を借りる

JITコンパイルさせれば．．．

import numpy as np
from numba import jit

NUM=100000000

@jit
def jit_no_for():
    xs = np.random.rand(NUM)
    ys = np.random.rand(NUM)
    square = xs*xs+ys*ys
    counter = np.sum(square < 1.0)
    pi = 4.0*counter/NUM
    print(pi)

def main():
    jit_no_for()

if __name__ == '__main__':
    main()

$ time python monte.py
real    0m5.836s
user    0m3.922s
sys     0m1.497s

よし，早くなりました．

実は．．．

元々のloopのコードに対してJITコンパイルをかけます．

import random
from numba import jit

NUM = 100000000

@jit
def jit_for():
    counter = 0
    for i in range(NUM):
        x = random.random()
        y = random.random()
        if x*x + y*y < 1.0:
            counter += 1

    pi = 4.0*counter/NUM
    print(pi)

def main():
    jit_for()

if __name__ == '__main__':
    main()

time python monte.py
real    0m3.122s
user    0m2.855s
sys     0m0.162s

ええ．．．結果的にはこの例では@jitのデコレータをつければよかったというオチです．

このケースではNumbaが勝ちましたが，ケースによってはNumpyで書いた方が勝つという場合もありますので，Numbaが最強というわけではないです．最適化できない場合は逆にJITのコストで遅くなることもあるので，@jit つければいいや．．．というものでもないです．　オールマイティで見ればNumpyが使いやすいですね．　Cythonはまだ勉強不足なので高速化の手段の一つとして勉強の一歩とできればと思います．

以上手を抜いて結構速くすることができました．