モンテカルロ法による円周率の求解によるベンチマーク with Python(Numpy, Numexpr, Numba)

モンテカルロ法による円周率の求解によるベンチマーク with Python(Numpy, Numexpr, Numba)

参考:

GoとPythonとGrumpyの速度ベンチマーク ~Googleのトランスパイラはどれくらい速い?~-Qiita

参考記事のGoの0.050 [sec]を目標に.

100回の平均.経過時間の計測に,平均処理の部分は除外する.

・普通に

import random
import statistics
import time

elpst = []
res = []
for _ in range(100):
    t0 = time.process_time()
    c = 0
    num = 300000
    for _ in range(num):
        x, y = random.random(), random.random()
        if x*x + y*y <= 1.0:
            c = c + 1.0
    t1 = time.process_time()
    res.append(c)
    elpst.append(t1-t0)

print(f'result: {4.0*statistics.mean(res)/num}, elapsed time: {statistics.mean(elpst)} [sec]')

result: 3.1415201333333336, elapsed time: 0.2465625 [sec]

・内包表記で

elpst = []
res = []
for _ in range(100):
    t0 = time.process_time()
    c = 0
    num = 300000
    c = sum([1 for _ in range(num) if random.random()**2 + random.random()**2 <= 1.0])
    t1 = time.process_time()
    res.append(c)
    elpst.append(t1-t0)

print(f'result: {4.0*statistics.mean(res)/num}, elapsed time: {statistics.mean(elpst)} [sec]')

result: 3.141735333333333, elapsed time: 0.25171875 [sec]

ifが律速な上,「sum」が遅いので逆に遅くなる結果に.

・numpyを使った場合

import numpy as np

num = 300000
elpst = []
res = []
for _ in range(100):
    t0 = time.process_time()
    x, y = np.random.rand(num), np.random.rand(num)
    c = ((x**2 + y**2) <= 1).sum()
    t1 = time.process_time()
    res.append(c)
    elpst.append(t1-t0)

print(f'result: {4.0*statistics.mean(res)/num}, elapsed time: {statistics.mean(elpst)} [sec]')

result: 3.14156, elapsed time: 0.0178125 [sec]

Goの約3倍速い.numpy.sum()よりも, numpy.ndarray.sum()の方が20μs程速いっぽい.

・Numbaで高速化(余り使った事がないので合っているか分からない)

import numpy as np
from numba import jit


@jit
def calc_pi(num = 300000):
    x, y = np.random.rand(num), np.random.rand(num)
    c = ((x**2 + y**2) <= 1).sum()
    return c


num = 300000
elpst = []
res = []
for _ in range(100):
    t0 = time.process_time()
    c = calc_pi()
    t1 = time.process_time()
    res.append(c)
    elpst.append(t1-t0)

print(f'result: {4.0*statistics.mean(res)/num}, elapsed time: {statistics.mean(elpst)} [sec]')

result: 3.141974, elapsed time: 0.01328125 [sec]

Goの約4倍速くなった.

numpyには確か「import numpy* as ne」してne.~~~ってやると高速化出来るテクニックがあった気がする(ドキュメントをぼーっと読んでいた時にたまたまみた).完全にうろ覚えで情報も出てこないので試せないけど,分かったら試して追記する.

追記:
Numexprか.「Numexpr is a fast numerical expression evaluator for NumPy」という事で,例えば「(x2 + y2) <= 1」の部分なんかは置き換えれば速くなるのかもしれない.

・Numpy

import numpy as np
import numexpr as ne


def calc_pi(num = 300000):
    x, y = np.random.rand(num), np.random.rand(num)
    c = (x**2 + y**2 <= 1).sum()
    return c


num = 300000
elpst = []
res = []
for _ in range(1000):
    t0 = time.process_time()
    c = calc_pi()
    t1 = time.process_time()
    res.append(c)
    elpst.append(t1-t0)

print(f'result: {4.0*np.mean(res)/num}, elapsed time: {np.mean(elpst)} [sec]')

result: 3.14133168, elapsed time: 0.01640625 [sec]

・Numpy + Numexpr

import numpy as np
import numexpr as ne


def calc_pi(num = 300000):
    x, y = np.random.rand(num), np.random.rand(num)
    c = (ne.evaluate('x**2 + y**2 <= 1')).sum()
    return c


num = 300000
elpst = []
res = []
for _ in range(1000):
    t0 = time.process_time()
    c = calc_pi()
    t1 = time.process_time()
    res.append(c)
    elpst.append(t1-t0)

print(f'result: {4.0*np.mean(res)/num}, elapsed time: {np.mean(elpst)} [sec]')

result: 3.1417125466666667, elapsed time: 0.015234375 [sec]

・Numpy + Numexpr + Numba

import numpy as np
import numexpr as ne
from numba import jit


@jit
def calc_pi(num = 300000):
    x, y = np.random.rand(num), np.random.rand(num)
    c = (ne.evaluate('x**2 + y**2 <= 1')).sum()
    return c


num = 300000
elpst = []
res = []
for _ in range(1000):
    t0 = time.process_time()
    c = calc_pi()
    t1 = time.process_time()
    res.append(c)
    elpst.append(t1-t0)

print(f'result: {4.0*np.mean(res)/num}, elapsed time: {np.mean(elpst)} [sec]')

result: 3.1423733333333335, elapsed time: 0.016234375 [sec]

・Numpy + Numba

import numpy as np
from numba import jit


@jit
def calc_pi(num = 300000):
    x, y = np.random.rand(num), np.random.rand(num)
    c = ((x**2 + y**2) <= 1).sum()
    return c


num = 300000
elpst = []
res = []
for _ in range(1000):
    t0 = time.process_time()
    c = calc_pi()
    t1 = time.process_time()
    res.append(c)
    elpst.append(t1-t0)

print(f'result: {4.0*np.mean(res)/num}, elapsed time: {np.mean(elpst)} [sec]')

result: 3.141626733333333, elapsed time: 0.010953125 [sec]

Numexprを使うとNumbaによる恩恵が得られない様だ.

まとめ(Goの0.050 [sec]を1倍として):
Python標準
 0.2465625 [sec](約5倍遅い) 
Numpy
 0.01640625 [sec](約3倍高速)
Numpy + Numexpr + Numba
 0.016234375 [sec](約3倍高速)
Numpy + Numexpr
 0.015234375 [sec](約3倍高速)
Numpy + Numba
 0.010953125 [sec](約5倍高速)

追記:
random.random()にしても,np.random.randにしても,ハーフオープンなので(),結果が微妙に確からしくないのでは,と思い,100万点を1万試行して結果をみてみた.

calc_pi: 3.1415980128000003
calc_pi2: 3.1415956308
円周率:3.141592653589793

誤差の範囲内.

カテゴリー: 未分類 パーマリンク

モンテカルロ法による円周率の求解によるベンチマーク with Python(Numpy, Numexpr, Numba) への2件のフィードバック

  1. ピンバック: pandasのstack/unstackはあまり効率的な方法ではない? | 粉末@それは風のように (日記)

  2. ピンバック: ColaboratoryでNumbaを使う | 粉末@それは風のように (日記)

コメントは受け付けていません。