WELL: Mersenne Twisterよりも優れた擬似乱数生成器

現在、コンピュータ上のシミュレーション用の擬似乱数生成器としては、Mersenne Twisterがよく用いられています。たとえばMersenne Twisterは、RubyやRの組み込み擬似乱数生成器として採用されています。

ですが、Mersenne Twisterには結果となる数値のバラつき方などに欠点があるのだそうです。次の論文では、よりよい擬似乱数生成器として「WELL」を提案しています。

以下、自分が理解したことを書きます。

循環する間隔の長さ

Mersenne Twister、WELLともに、生成される乱数列は2k回で循環するという特徴を持っています。kが小さいと、同一の乱数列が繰り返し現れてしまうため、嬉しくありません。つまり、kは大きければ大きいほど優れています。

Mersenne Twisterは、k=19937であり、充分すぎるほどに優れています。

バラつき (equidistribution)

生成される乱数を[0,1)の区間の実数に正規化した時、t個の連続する乱数からなるベクトルはt次元の単位超立方体上の点としてマッピングできます。乱数列は2k回で循環する、つまり2k種類の初期状態を持ち得るので、擬似乱数生成器が発生し得る点は2k個存在します。

ここで、[0,1)の区間を2l個に均等分割すると、t次元の単位超立方体は2tl個の超立方体に分割されます。擬似乱数生成器が発生する点は、この分割された超立方体それぞれに、同じ数だけ、つまり2k-tl個ずつ配分されることが望ましい。バラついていることになるから。あるt, lについて、各超立方体に同じ数の点が配分される場合、(t,l)-equidistributedと呼びます。

擬似乱数生成器がwビットの乱数を生成する場合、考える必要があるのはl≦wの場合だけです。たとえばw=32, l=64について考えると、乱数が32ビット値なのに、[0,1)の区間を264個にぶったぎったとしても、半分の区間には一個も乱数が入らない232個の区間にしか乱数が入らない*1。つまり同数だけ配分されることはありえないので、考えても意味が無い。

同様に、l≦floor(k/t)の場合だけ考えればよいことも分かります。l>floor(k/t)だと、分割された超立方体の数が2k個を超えてしまうので、すべての点を配分することができない。またlを固定すると、tについて考える必要があるのはt≦floor(k/l)の場合だけです。

ここで、あるlについて、(t,l)-equidistributedとなる最大のtをtlと呼びます。tlはfloor(k/l)であることが望ましい。また、δlとΔ1を次のように定義します。

  • δl=floor(k/l)-tl
  • \Delta_1=\sum_{l=1}^{w} \delta_l

Δ1が0である場合、連続いくつの乱数を取ったとしてもk/w個くらいまでの乱数をとったときに*2、それによって構成されるベクトルは、一様にバラつくことになります。このような場合、擬似乱数生成器はmaximally-equidistributed (ME)であると言います。

Mersenne TwisterはΔ1=6750であるため、MEではありません。つまり、連続する乱数によって構成されるベクトルは、超立方体上で偏って分布します。

状態遷移行列の固有多項式の非0係数

Mersenne Twisterや、本論文が提案するWELLなどの擬似乱数生成器は、いずれも次のような形式に一般化できます。

まず、擬似乱数生成器はkビットの状態を持ちます。初期状態x0からi回乱数を生成した後の状態xiを次の列ベクトルで表します。

  • {\bf x_i}=\left(\begin{array}{c} x_{i,0} \\ x_{i,1} \\ \vdots \\ x_{i,k-1}\end{array}\right)

ただし、xi,jは0か1、つまり1ビットを表します。

この時、状態xiからwビットの乱数のビット列yiを生成する操作はyi=Bxiと表せます。ただしBはw×k行列です。

また、状態xi-1からxiに移す操作はxi=Axi-1と表せます。ただしAはk×k行列です。

Aの固有多項式の0でない係数の数をN1とします。で、N1はk/2に近い方が良いのだそうです。

Mersenne Twisterはk=19937に対してN1=135なので、あまり良くない。

比較

WELLでは以上の指標を改善しています。

アルゴリズム w k Δ1 N1
Mersenne Twister (MT19937) 32 19937 6750 135
WELL512a 32 512 0 225
WELL1024a 32 1024 0 407
WELL19937c 32 19937 0 8585

これら指標以外のランダム性検定については触れられていません。

文中にあるWELL1024aのコード例を見ても分かる通り、簡潔に実装できることもよい特性だと思います。*3

性能については、Mersenne Twisterとそれほど変わらないとのことです。

ゼロからの脱出

Mersenne Twisterには、初期状態のビット列に0が多いと、しばらくの間、乱数のビット列に0が多い状態が続くという弱点があります。WELLにも同様の弱点がありますが、比較的早めに0の割合が小さくなります。

これは、初期状態をSHA-2で設定するなどの工夫をすれば、実用上問題にならない気がします。

*1:2014-01-01T00:54+09:00 編集。

*2:2014-01-01T01:05+09:00 編集。

*3:2014-01-01T01:25+09:00 編集。