準モンアロ法カノレ

伏見正則111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111
記の性質が次元 s によらないところにある.

1. はじめにところで，/が解析的性質が良い関数である場合には

J.vonNeumann と s. M. Ulam は乱

1945年頃，台形公式などの数値積分法を多次元に拡張したものを使数を使って決定論的な種々の数学的問題を解く方法を考って I の近似値を求めることもできる.しかしこの場合

え，このような方法の総称としてモンテカルロ法(以下にはつの次元方向を m 等分するとすれば，全部でMC 法と略す)とし、う名称を用いた.現在では，対象が(m+l)' 回 f の値を計算しなければならず，本来確率的な要素を含む場合も含めて，乱数を用いて実次元 s とともに指数関数的に増大すると L 、う欠点があ験する方法のことを MC 法と呼ぶことが多い.

る.

MC 法では，何回か実験を繰り返し，各国に得られた計算量がそこで，虫のよい要望が出てくる:上の 2 つの方法の測定値の算術平均をもって，知りたい特性値の推定値と“中間的な"方法で，するのが普通である.

計算の回数 N を急激に大きくしなくてもよく，かつ，1 回の実験では区間 [0 ，lJ 上の有効数字の桁数を増すために関数s 個の一様乱数を用いて測定値を得るものとすれば，知計算量が次元とともにあまり増えない方法はないものだ

りたい特性値は形式的に s 重積分

ろうか?一一このような要望にある意味で応えるのが準

I=~>~:六 x.. ...， xs )仇仏(1 )の形に書き表わされる.われわれは，これをNA N = 三-;-L:ん( 2).
J.Vn=lによって推定する.

ただし，んは n 回目の実験で使用モンテカルロ法 (Quasi-Monte C
arlo method ，以下QMC 法と略す)と呼ばれるヰよりな(耳なれない?)方法である.これにはいくつかの方法があり，関数 f の解析性がきわめて良い場合に有効な優良格子点 (goodlattice points) 法もそのひとつであるが，本稿では解した乱数の組に対応する f の値である.確率変数として析性が良くなくても使える“差異の小さい点列 (low­の AN の期待値 E(A N ) は I に等しく

discrepancy sequence)

は I の不偏推定量であり)，(すなわち AN

各回の実験で使用する乱数の組が互いに独立で・あれば ， AN の分散 Var(AN) はn2V叫 A N )=ヤ2. 点列の差異と積分の観差(1)式の積分領繊，q2=~> -S:{f(X.. ..., xs) ー I}2 dx l...dxs、

MC 法の欠点は，誤差 (A N ー 1) の目安としての ANの標準偏差 .;V示目玉了が ， Nを増していったときにN"~ という遅い速度でしか小さくならず，たとえば有効数字を十進で 1 桁増やすためには N を 100 倍に増やさなければならないというような効率の惑さである.一方，
MC 法の長所は，関数 f が微分できないとか不連続であるなど，解析的には性質が悪くても適用でき，また上

ふしみまさのり東京大学工学部計数工学科〒 113 文京区本郷 7-3ー l590(22)すなわち s 次元の単位超立方体

G 内の N 個の点の集合 SN={X.. … ， XN} の( 5 次元)
差異 (discrepancy)となることはよく知られている.

を用いる方法"について述べる.というのは，|ν (J;N)̲"n I

N3=sup| 一τ7一 -vol(J) !
J I1Vで定義される量である.ここで上限 sup は
J=[O， t!l x … X[0， t8)~CBという形のすべての(超)直方体にわたってとるものとする.また， ν (J;N) は J に含まれる SN の点の個数を表わし，vol(J) は J の体積である.

1 次元 (5=1) の場合には ， Ð'J' は SN が一様分布からのランダムサンプルと見なせるかどうかの検定に使われるコルモゴロ

フ・スミルノフ統計量に一致する .Djf》は，© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

その多次オベレーションズ・リサーチ

元への一種の拡張でありSN の G 内における一様性が成り立つならば，その点列の最後の座標成分を(上記差異の小さい点列が定積分(1)の評価に有効である置き換えて得られる N 個の点の集合の差異に対してはとされる根拠は，次の Koksma-Hlawka の不等式で(8)式の形の上界評価ができることが知られている.

を表わす 1 つの尺度である.のある.

(4 )から(7)を作ったのと同様に) (n-1)/N で係数 B.(6) を最小にするためには ， b , =2 , b z =3， ・II‑ANI~V(f)D'J>(3)

ここに V (f) は関数 f の Hardy-Krause の意味での変動と呼ばれ，C. における f の変化の激しさを表わというふうに ，b k として小さい方から h 番目の素数をogs)とるのがよい.しかし，そうしても log B.=O(slであり ， B. は s とともに急激に大きくなってしまう.

す量である.その厳密な定義は省略するが，連続性や微3) Sobol'の方法

分可能性といった条件が成り立たなくても有界になりう

無限列を生成する方法で，差異に関する評価式は( 5)と同じ形であるが，係数 B， が log B , =O(s loglogs)る量である.

となり3. 差異の小さい点列の例

s:<::5 では (6 )式の B. より小さくなる.系列の構成法は次のとおりである.

差異の小さ L 、点ヂu を構成する方法は L 、ろいろ提案されSη=(X~D ，… ， X~>) の各座標成分 zF 〉をそれぞれ独立に次の漸化式を用いて計算する.ているが，そのうちのいくつかを紹介しよう.

X~k>=O1) Hal旬n の方法
X n =( 仇 1(n-1) ， ・ "，'þÞ .(n ー 1)) ， n=I , 2 , …

ここに b" … ， b， は，(4 )どの 2 つをとっても互いに素な自然数の組であり，仇 (n) は radical-inverse functionと呼ばれる関数でありn の b 進展開を n= L:. j=oajbiとすると ， ，þÞ (n)= L:. j=oa j b-j-1

である . (4) は無限点列であり，最初の N 点からなる集合の差異は次の不等式

X~")， =x~k>EÐv~Nn>' 1 豆島孟 s， n ミ(9 )ここに， EÐ は 2 進法での繰り上りなしの足し算を表わす.また m(n) は ， n の 2 進展開を… dsd z んとすると

き ， dm=O となる最小の m を表わす.uLf} は次の漸化式によって定まる 2 進小数である.

(記号を簡単にするため，上っき添字 h は省略する)
世間 =C ， vm- ， EÐczvm-z ';B' ・・ EÐc p - ， v m - p + 1 EÐv明 -pEÐを満たす.
1一九[v m _p /2 P] , m:<=:p+1
ただし係数 Cj は ，一一gh(6)法

方yの・4sHb 咽
一戸 E

B一h 畑
九一
一円以剛

DW~B.N-l (l og N)'+O(N‑l(logN)8‑l) (5)(10)

zP+C， ZP-l 十… +CP_ ， Z 十 1 がガロア体 GF(2) 上の原始多項式となるように定め，また各次元長ごとに異なる原始多項式を使用するものとする.
漸化式( 10) を使うための初期値世間 (1ζmζ p) は ，Xn=(州一 1) ，"'， Øb._l(n-1) ， 竺示)，n=l , …, N

<1 で，かっ小数点以下第 m ピットが 1 ，それより下位(7)b" …, b._ ， に対する条件は 1) と同様であり，差異はD'
J>~B，_， N-l (l og N)H+O(N‑l(logN)B‑2)(8)を満たす.のピットはすべて 0 となるように定める.

4) Faure の方法
以上の 3 つの方法で構成される点列の差異の上界評価
式の主要項の係数 B. は，いずれも s→∞のとき無限大になる.これに対して，本方法によって得られる系列にこれらの点列を使うと，誤差の上界評価(3)がo(N‑l(logN)') あるいは O(N-l( logN) &-1) となり
MC 法における誤差の確率的評価(標準偏差)がO(N ちであるのといちじるしい対照をなす.また，対する B. は

3=一一一一一一一 ，

1 (れー 1\1B，=~( ー」一二一 )16(log2)2' -,- s!¥2logr.J(s 二三 3)

rs:S 以上で最小の素数1)と 2) を比べると， 1) はオーダの意味で劣っており，であり，すべての s:<=:2 について他より小さく使い途がないと思われるかもしれないが ， N を順次増しのとき B.→O となる.ていき，任意のところで止められるというメリットがある.なお，一般に，ある無限点列の差異に対して(5 )の形の上界評価(係数 B， は上記のものと違ってよい)1991 年 12 月号

vm

与えられた次元 s に対してs→∞s 以上で最小の素数(上記の η) を r と書くことにする • Xn(n 二三 1) の各座標成分 zLK 〉はr 進小数として次のように定める.

© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

(23)591

XN-tl は fo(x) の選び方によらないので，いちど "1 ，...，zr〉 =Eou;K3(n) 「j-1Up

ここに叫kl(n) は ， n-1 の r 進展開を座に点列の生成を開始できるのである.

4.
2 理鵠的背景n‑1=:
Eaj(n)rij=O

表 1 に示す p， q および M(x) ， g(x) は次のようにとすると，して選ばれたものである.まず， GF(2) 上の p 次の任

y~ρ (n)= ゑ(])(k ー 1 )1叫 n)で与えられる.

以後はこれを使って即を計算して記録しておけば，modr(便宜上 ， k=1 の場合は0 0 =1

意の原始多項式 M(x) と ，と解釈す

p-1 次以下の任意の多項式g(x) を用いて上記のようにして生成される点列 SN のs 次元“メリット数"を次弐により定義する.

るものとする)pω=min:E {deg(hk(x))+1}

4. 高速に発生できる系列(14)k=l

前節で述べた方法には，いずれも点列の生成に時間がここで min は方程式:
Egk‑l(x)hk(x)=O modM(x)k=lかかると L 、う欠点がある.差異の小さい点列の生成が遅(15)いのであれば， QMC 法を用いるよりも通常の乱数によのすべての非零多項式解 (h 1 (x) ， … ， h， (x)) についてとる MC 法を使った方が同一時間内にもっと精度の良いるものとし，また deg(O)= ー 1解が得られると L 、う場合もありうる.そこで本節では，と，メリット数と差異の聞には次の関係がある.

最近筆者ら [9J が考案した高速に生成できる差異の小さ L 、点列を紹介する.生成に要する時聞は，一番速い通常の乱数の生成とほとんど同じである.
4.1)‑1(Un=Un_pEÐun_p+q , n~ ρ+1

したがって差異を小さくするためにはメリット数を大+1/Ar(x)))

ただしパラメタの組(ム q) は表 1 から

n=I ， 2 ， … ， 2 P ー 1

これに XO=O を加えた N=2P 個の点を使用する.

漸化式 (11) を使用するための初期値 U 1.… ， u p の設定法の原理は次のとおりである.まず GF(2) 上の p 次未満の多項式 fo(x) を任意に選ぶ(たとえば fo(x)=1). 次に表 1 の多項式の組 M(x) ， g(x) を使って，多項式の列 fl(X) ， … .fp(x) を次の漸化式により生成する.

mod M(x)そして，各 n についてん (x) 川l(x)の部分商 AJ(x) ， … ， Ar(x) と次の関係があることが知られている.

(12)

ん (x)=g( 叫ん -l(X)

メリット数は g(x) 川I/ (x) の連分数

展開
g(x)/M(x)=I/(A 1 (x)+!/(A 2 (x)+ ・ .
.

選ぶ.これを用いて点列 {X n } を次のように定める.

(13)
を形式的に

Ln(x)=:El~n) x-j という形のローラン級数に展開し，
j=lD)
J'=O((1ogN) I-1 /2P <Bl), N=2P2 次元の場合，2 進小数の列 {U n } を漸化式

Sη =(Un ， Uη+ 1.… ， U n +ト tl ，

そうするきくすればよいことになる.

生成法

によって作る.と定義する.

・p(2l=p+2‑ max deg(At(x))
l:>:i:>:r
したがって，(16)ρ を固定した場合，部分商がすべて 1 次式 (x または x+1) となる多項式の対 (g， M) がメリット数を最小にする.このような多項式対をフィボナッチ多項式対と呼ぶ.

(この命名は，フィボナ ':1 チ数列(1，1 ， 2 ， 3， 5 ， 8 ，…)の連続する 2 項の比の正則連分数展開の部分商がすべて l に等しいと L 、う事実との類似による).
フィボナッチ多項式対は次の漸化式により生成できる.

Fo(x)=I , FJ(x)=AJ(x)Ft(x)=At(x)Fト J(x)+F，← 2(X) ，
ここに ，

注2

(17)

At(x)(i~l) は z または x+1 のいずれでこれから得られる 2 進小数 L n (2) を適当なピット数もよい.

(1 6 ピット， 32 ピットなど)で‘打ち切ったものを h とす

Ft_J+Fト 2 )と違 L 、 ， At(x) として 2 つのものを使える

フィボナッチ数列の生成に使う漸化式(Fj=る.このように初期値の設定にはやや時間がかかるが，ので，これが完了すると，あとは 1 つの X n を生成するのに lであり.フィボナッチ多項式対 (Fト l(x) ， Ft(x)) は同回の@演算をすればよいので，きわめて高速に点列が生じ次数のものが一般に多数あることに注意されたい.

成できる.

しかも

592 (24)

N 個の点の集合 SN={XO ， X 1. …，(1 7)式によって生成されるのは多項式の二分木

一方 ，p 次の原始多項式 M(x) と ，© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

p- 1 次以下の多オペレーションズ・リサーチ

表 1(18) 式を満たすフィボナッチ多項式対 G(p， q) :(g, M)(数字は非零項のベき指数を表わす)

4 5 6 11 14 15 16 17
9 12 15 169・一
onzo--onu

守sq3'ir

一一

haaヲ

Gnu-o''

Ma 一M'lF2Rh

GGGGGGGGG

5 7 9 1
1 12 14 15
3 5 10 1
1 12 13 14

2 3 4 5 8 10 13 14 18

3 4 6 8 12 14 15 17

g一

2 4 6 10 12 13 14 15 16 20
3 4 5 6 7 9 1
0 16 17 19AU'a
一コ。

3581
49Meg 一9 創刊

乙5 6 7 9 10 12 1
3 14 15 16 18 19 22
3 6 8 14 16 19 21
1 13 14 16 17 18 19 20 21 23
4 5 7 8 9 1onυ
一・3e3r一一一

3 6 7 8 9 11 18 226 9 1
1 14 16 18 19 23 259・一Gnu-o'A

1 12 13 14 21 24
3 7 9 13 4 5 8 9 10 1
1 12 15 20 21 22 23 24 26 27 28
1 12 13 14 15 16 1
7 18 19 21 22 24 26 27
2 3 4 5 6 7 9 10 1Mg。

3 5 6 7 8 12 13 15 16 18 19 20 23 24 25 26 27 29 30 31
1 12 13 14 22 23 30
3 8 9 10 1項式 g(x) を用いて，漸化式(1 3) によって生成され=Fp(x) が (18) 式を満たすものを探せばる多項式の系列{ん (x)} から，前記のようにん (x)/リット数が最小の数列を高速に生成できることになる.

M(x) のローラン級数展開を基にして 2 進小数の系列

表 1 は，{U n } を作ったとき，った結果得られたものである.ただし，表に示しである

これが漸化式(11)を満たすための条件は2 次元メ10::;;p::;;31 の範囲でこのような全数探索を行な

各 P について複数の対が見つかったので，ある基準によgP(x)+g<l(.
x)+I=O modM(x)(18)が成り立つことである.ただし x P +x q +1 は GF(2)り選んだ 1 対ずつを挙げてある.
3 次元以上のメリット数については(1 6) 式のような上の原始 3 項式である.したがってフィボナッチ多項式簡単な表現は知られていないので，対 (Fp_l(x) ， Fp(x)) の中で ，式にもとづいて計算する必要がある.表 2 は，このようg(x)=Fp_I( .x), M(x)(14) および(1 5)にして得られたメリット数を示したものである.
G(p， q) により生成される点列の

表 25A
r

< >

ρ・ρ・< >

4a1 p<Sl̲
1￨

メリット数

< >

p<6'

1612

1
177G (I 7 , 5)

18

14

1
2

1
17G(18, 7)

19

14

13

12

1
1G(20, 3)

2
1

14

14

1212G(22 , 1)

23

17

1
715

13G(23, 5)

24

16

15

1515G(25 , 3)

26

20

19

1715G(28, 13)

29

24

23

18

18G(31 , 6)

32

24

24

2220

1991 年 12 月号

QMC 法と MC 法の誤差の比較(数値は相対誤差 x 10・を示す)

ヨ瓦旦竺!

?日斗
線形合同法

G(15 , 1)表 3①②③④⑤2140293710373 19816 16503

1432

4084

2612 34696

3698

583

3101

2763 34342

3610

5002393

853

8887

2510

115

369365275

1577

201463© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

(25)593

pseudorandomnumbers.MathematicsofCom.

4.
3 数値計算例

putation , Vol. 54 (1990) , 7
37‑
748.

本節で述べた方法の効果を確かめるために，表 1 の.L.:Algorithm659:[2J Bratley , P. , andFox , B

G(17 ， 5) を使って次の 5 つの 5 変数関数について計算

Implementing Sobol's quasirandom sequenceを行なってみた.

(exp( ‑ :i=1Xk)(XI X2 Xax • x s exp(:E ~=1 xi)(exp( ‑2: ~=1 X k)sin(2:%=1Xk)

generator. ACM Trans. Mathematical Soft.
即are ，

主 un syst色me denum駻ation(endimensions).

④行平2: i=三E

Acta Arithmetica , Vol.41 (1982) , 3
37
‑
351.

1+xl+2xIX~+3xIX~X;+4xIX~X;X:(Vol.14 (1988) , 8
8
‑
100.

uitesassoci馥s[3J Faure , H.:Discr駱ance de s

[4J Mullen , G. L. , andNiederreiter , H.:Opti
ｭmal characteristic polynomials fordigital+ 5XIX~X:X:x~
表 3 は，相対誤差 I( A N ー 1)/11 のmultistep pseudorandom

10 6 倍を示したもの

numbers.

白m­puting , Vol.39 (1987), 1
55‑
163.

である.比較のために，使用した計算機 (Sun-4) に内[5J Niederreiter, H.:Quasi.Monte Carlo meｭ

蔵されている線形合同法

thods and pseudo‑random numbers. Bulletin

Xn=2736731631558Xト 1 十 138 mod2岨.
ofthe American Mathematical Society , Volによる乱数生成関数 drand 48( )を使って得られた結84 (1978) , 957‑1041.

果も示してある.後者については，乱数の初期値 X o を[6J Niederreiter , H.:Pointsetsandsequences

5 通り変えて実験を行なっている.

withsmalldiscrepancy.Monatsheftef

表 3 に見られるとおり，本節で述べた方法による誤差は，通常の MC 法による誤差よりはるかに小さい.

MCMaｭthematik , Vol.104 (1987), 273‑337.

法では，誤差を 1 桁小さくするのに 100 倍のサンプル数. M.:The distributionofpoints[7J Sobol' , 1を必要とすることを思い起こせば，われわれの方法の有inacubeandtheapproximateevaluationof

効性がさらにはっきりする.

integrals. USSR Computational Mathematics

5.

andM athematical Physics , Vol.7(1967) , 86‑

おわりに

112.

差異の小さい点列を用いる準モンテカルロ法を紹介[8J Tezuka , S.:On the discrepancy ofGFSR

し，特にわれわれが最近考察した点列の高速発生法とそ

pseudorandom numbers. Journal of ACM ,の効果について述べた.これは，Vol.34 (1987), 939‑
949.

}]IJ 稿の著者手塚集民と[9J Tezuka , S. , and Fushimi , M.:Calculationの共同研究によるものであり，同氏の貢献に感謝する.

ofFibonaccipolynomialsf
orGFSRsequences

参宏文献

with lowd
iscrepancies. Submitted.

[1J Andre , D.A.
" Mullen , G.L. , andN
iederreiｭ
ter, H.:Figuresofmeritf
ordigitalmultistep[IOJ 水田成仁:差異の小さい点列の高速発生法に関する研究.東京大学工学部計数工学科卒業論文，

1991.

考多594(26)© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

オベレーションズ・リサーチ