インテル® Distribution for Python* におけるより高速な乱数生成

インテル® oneMKLインテル® ディストリビューションの Python*

この記事は、インテル® デベロッパー・ゾーンに公開されている「Faster random number generation in Intel® Distribution for Python*」(https://software.intel.com/en-us/blogs/2016/06/15/faster-random-number-generation-in-intel-distribution-for-python) の日本語参考訳です。


インテル® Distribution for Python* 2017 Beta Update 1 では、Numpy* の拡張として、numpy.random の設計を反映し、パフォーマンスを向上するためインテル® MKL のベクトル統計ライブラリーを使用した numpy.random_intel が追加されました。

numpy.randomnumpy.random_intel に置換するだけで、パフォーマンスを向上できます。

In [1]: import numpy as np

In [2]: from numpy import random, random_intel

In [3]: %timeit np.random.rand(10**5)
1000 loops, best of 3: 1.05 ms per loop

In [4]: %timeit np.random_intel.rand(10**5)
10000 loops, best of 3: 146 µs per loop

出力および時間測定に使用したシステム構成:

ソフトウェア: インテル® Distribution for Python* 3.5.1 2017 Beta Update 1 (Jun. 8, 2016)、Numpy* 1.11.0、インテル® MKL 11.3.3。
ハードウェア: インテル® Xeon® プロセッサー E5-2698 v3 @ 2.30GHz (2 ソケット、16 コア、ハイパースレッディング無効)、64GB RAM。
オペレーティング・システム: Ubuntu* 14.04 LTS。

次の表は、一般的な分布におけるパフォーマンスの向上を示しています。

100,000 の確率変数のサンプリングを 256 回繰り返した時間測定の合計
分布 timing(random) (秒) timing(random_intel) (秒) スピードアップ (倍)
uniform(-1, 1) 0.357 0.034 10.52
normal(0, 1) 0.834 0.081 10.35
gamma(5.2, 1) 1.399 0.267 5.25
beta(0.7, 2.5) 3.677 0.556 6.61
       
randint(0, 100) 0.228 0.053 4.33
poisson(7.6) 2.990 0.052 57.44
hypergeometric(214, 97, 83) 11.353 0.517 21.96

numpy.random_intel モジュールを使用すると、インテル® MKL でサポートされている、(RandomState クラス・コンストラクターの brng 引数、または seed 初期化メソッドを使用して指定できる) 異なる基本擬似乱数ジェネレーターを利用することができます。デフォルトの基本擬似乱数ジェネレーターは、numpy.random と同じです。

次のコードは、異なる基本擬似乱数ジェネレーターで 100,000 の標準一様確率変数のサンプリングのパフォーマンスを測定し、メルセンヌツイスター基本乱数ジェネレーター ‘MT19937’ と比較します。

import numpy as np
from timeit import timeit, Timer
from numpy import random_intel
from operator import itemgetter

def timer_brng(brng_name):
    return Timer('numpy.random_intel.uniform(0,1,size=10**5)', 
        setup='import numpy.random_intel; numpy.random_intel.seed(77777, brng="{}")'.format(brng_name))

_brngs = ['WH', 'PHILOX4X32X10', 'MT2203', 'MCG59', 'MCG31', 'MT19937', 'MRG32K3A', 'SFMT19937', 'R250']
tdata = sorted([(brng, timer_brng(brng).timeit(number=1000)) for brng in _brngs], key=itemgetter(1))

def relative_perf(tdata):
	base = dict(tdata).get('MT19937')
	return [(name, t/base) for name, t in tdata]

relative_perf(tdata)

出力は次のとおりです。

In [9]: relative_perf(tdata)
Out[9]:
[('MCG31', 0.607988362606216),
 ('SFMT19937', 0.6520599397085953),
 ('MCG59', 0.6522106463148241),
 ('MT2203', 0.8041550837721791),
 ('MT19937', 1.0),
 ('PHILOX4X32X10', 1.021985386332785),
 ('R250', 1.2991300341128584),
 ('WH', 1.4285221807604567),
 ('MRG32K3A', 2.148446327408357)]

numpy.random_intel は、leapfrog および skipahead メソッドを使用して RandomState クラスを初期化し、異なるスレッドで使用されるランダムストリームが相関していないことを保証します。

次に例を示します。

rs1 = np.random_intel.RandomState(77777, brng='SFMT19937')
rs2 = np.random_intel.RandomState(77777, brng='SFMT19937')
# skip the first stream 2**60 steps ahead
rs1.skipahead(2**60)

s1 = rs1.randint(0,2**31,size=10**5)
s2 = rs2.randint(0,2**31,size=10**5)

ピアソンの相関係数を使用して、2 つのサンプルが相関していないという帰無仮説をテストします。

from scipy.stats.stats import pearsonr
pearsonr(s1, s2)

出力は次のとおりです。

In[12]: pearsonr(s1, s2)
Out[12]:
(0.0034069712002612325, 0.28131564831676692)

データが仮説と矛盾しないことを示しています。

詳細についてはモジュールのドキュメントを、ベクトル統計についてはインテル® MKL のドキュメントをそれぞれ参照してください。

皆様からのフィードバックをお待ちしております。

コンパイラーの最適化に関する詳細は、最適化に関する注意事項を参照してください。

タイトルとURLをコピーしました