明細書
準乱数生成装置及び方法並びに関数 f の多重積分計算装置及び方法 [技術分野]
この発明は、 一般化されたニーダーライタ ' シーケンス
(General ized Niederrei ter Sequence) 及びフォーリレ · ンーケンス (Faure sequences) と含むローデイ スクレノ ンシ · シーケンス V low-discrepancy sequences, 準 BL ともいつ。 ) を筒速 生成する 方法及びシステムに関するものである。 この準乱数とは、 極めて一様に 分布した点列である。
[背景技術]
最近の金利自由化の傾向にともない、 銀行や証券会社などのさまざま な金融機関では、 金融派生商品の価格づけなどの目的で、 高速なコンビ ュ一タを駆使するようになってきている。 高速なコンピュータが使用さ れる理由は、 様々な顧客每のパラメータや経済指標などの極めて多数の 変動する要因があるにも拘わらず、 顧客の要望に迅速に応える必要があ るからである。 このような金融向けの数学的モデルを提供したり、 その数学的モデル における方程式を合理的な手法で解く アルゴリズムを与えるための分野 と して、 金融工学 (Financial Engineering) が提唱され、 近年、 盛ん に研究されるようになってきている。 金融工学、 特に金融派生商品の価格づけで頻繁に必要となる計算に、 多重積分計算がある。 一般的に、 コンピュータを用いて積分計算を行う
には、 Si分区問を微小な領域に等分し、 その個々の区 (HJにおける関数の 値と区間の幅の檳を累算する、 という技法が用いられる。 しかし、 多重 檳分計算の場合、 もし k重楱分計算を行うとすると、 1 の区間の分割数 が Nであるとき、 N kもの反復計算が必要になり、 k = 5程度でも、 実用 的な分割数 Nに対して、 通常のコンピュータでは、 実際的な計箅時間で 処理できる範囲を超えてしまう。 このため、 一種の乱数を用いたモンテカルロ (Monte Car l o ) 法と呼 ばれる多重楱分計算のための方法が以前から提案されている。 この方法 は、 多重積分計算の k次元の積分区間の範囲内のー铩乱数の k次元座標 を多数生成し、 生成された座標を積分される関数に頓次代入して、 その 代入した値が所定の条件を満たす確率を求める、 というものである。 れを数式で表すと、 次のとおりである。
この場合、 N個の点柒合 x . ( i =1 N ) を与えることがアルゴリ ズムそのものである。 金融派生商品の価格づけにおける価格の期待 iH F; I" ^には、 この式において、 それぞれの確率モデルから決まる Γ ( X ) が 用いられることになる。 また、 この嵇の問题では、 i分 15差を小さ くするようなアルゴリ ズム が 「良いアルゴリズム j ということになる。 一方、 誤差の方を一定 (例 えば 0. 001 ) にして考えれば、 できるだけ少ない点集合 (すなわち、 小
さい N) でその誤差を達成できるのが
1"良いアルゴリズム j である。 こ のようなアルゴリズムでは、 一般的に点の数に比例して計算時間がかか るので、
Γ良いアルゴリズム」 はすなわち
Γ速いアルゴリ ズム」 でもあ る。 線形合同法などの方法によつて発生されたランダムな点が使用される モンテカルロ法では、 誤差は、 ほぼ Νの平方根の逆数に比例することが 知られている。 これに従えば、 簡単に言うと、 0.001の誤差を達成する には、 10
6個の点が必要である。 一方、 準乱数を Ν個の点集合として用いるアルゴリズムもある。 この 場合は、 理論的には、 誤差がほぽ 1 /Νに比例することが知られている, これだと、 0.001の誤差を達成するには、 Ν =10
3個の点で十分というこ とになる。 従って、 準乱数を使用すると、 モンテカルロ法よ りも、 計算 がおおよそ 10
3倍速い、 ということになる。 準乱数と して、 幾つかの方法が提案されているが、 その中でも、 フォ 一ル · シーケンスとよばれる数列が、 この分野の専門家によって、 高い 次元の多重積分計算に対して有望視されている準乱数である。 これの詳 細については、 以下の論文を参照されたい。
P.P. Boyle, New Variance Reduction Methods for efficient Monte Carlo Simulations, Conference on Advanced Mathematics on
Derivatives, RISK Magazine, (September, 1994). また、 この数列のアルゴリズムとして、 下記の論文に記述されている 技法が知られており、 また、 現在使用されている。
P. Brat ley, B. L. Fox, and H. Niederreiter, Implementation and Tests of Low-Discrepancy Sequences, ACM Trans. Modeling and Computer Simulation, 2, (July, 1992), 195-213. しかし、 この論文に記述されている方法でフオール · シーケンスを発 生させると、 本願発明者の実験によれば、 線形合同法の乱数を発生させ るよ りも実に 60倍もの時間がかかってしまう。 特に金融派生商品の価格づけの計算では、 準乱数による多重積分計算 は、 多くの場合、 理論上、 通常のモンテカルロ法を用いた多重積分計算 よ りも 1 00倍程度の高速化を達成するが、 それは、 準乱数の 1 つの点 を発生するのに要する時間が、 モンテカルロ法における線形合同法によ る 1 つ点の発生に要する時間とほぼ等しいと仮定してのことであって、 もしフ ォール · シーケンスを発生させるのに、 線形合同法の乱数を発生 させるよ りも 60倍もの時間がかかってしまうと、 準乱数を使用するこ との利点が大幅に割り引かれてしまう。
[発明の開示]
本発明の目的は、 準乱数を高速で生成するための方法及びシステムを 提供することにある。 本発明の他の目的は、 そのようにして発生される準乱数を用いて、 金 融派生商品の価格づけなどに使用される多重積分の計算を高速化するこ とにある。 本発明によれば、 従来の数列 X = T V が、 Y = TG v で置き換え
られる。 但し、 この Gは、 その対角成分が全て 1 であり、 その
j 一 i = l を満たすような 成分が一 1 であり、 その他の成分が全て
0であるような正方行列である。 すると、 本発明の非自明な知見によれば、 π+l番目の準乱数は、 nの b進展開 (n , n n,.... η,) をベク トル成分と した ν = ( η,,....
Di m-1 1 n 1
. n ,, η ) において、 n .≠ b— 1 であるような最小の《i を求め、 生
ID - J Π J
成行列 Tの j列と、 n番目の準乱数を、 mo d bで加算するこ とによ つて、 求められる。 これは、 Y Al = Y ADD. T e .と表すことがで
nTJ n b j
きる。 このとき、 実質的に加算しか必要でないので、 実際の行列どう し の掛け算を要する方法よりは、 はるかに少ない計算量で、 高速に準乱数 を求めることができる。 また、 それらの計算を複数のプロセッサに行わせることによ り、 よ り 高速なシステム及び方法を提供する。 これは、 前述の行列 Gの性質を用 いるものであり、 Y =TG V よ り、 数 nを与えると任意の順番の準乱
n n
数 Y を生成することができる。 このことを用いて、 複数のプロセッサ π
に各プロセッサにて処理するための数 Sを与えれば、 その数 Sに対する
Υ を各々のプロセッサは生成する。 Υ が生成されれば、 上述した掛け
S S
算を行わない加算方式にて Υ +1を生成するこ とができるのでその処理 を繰り返し、 各プロセッサがその処理割り当ての終了の数 tに対する
Ytまで生成する。 そうすれば、 システム全体としてはよ り高速に準乱 数を生成することができるようになる。 但し、 被精分関数によっては準 乱数の計算よ り被積分関数の計算の方が非常に手間がかかる場合もある。 よって、 Y l = Y AD Dk T e .による掛け算を用いない加算方式で
n+] n b j
行っても、 Y = Τ G ν によ り掛け算を用いる方式を用いても、 処理時
n n
間に大差ない場合もあるので、 その場合には各プロセッサにおける Y から Y,の生成は掛け箅を用いる方式を用いて計算してもよい。
S+J t よ り具体的に本発明の構成を示す。
1次元の準乱数を生成するには、 以下のような構成が必要となる。 す なわち、 (a)n番目 (nは、 1以上の整数) の準乱数 Υπの值を b進数
(bは、 2以上の素数) 展開し、 生成された m個の成分
y y ) を順番に第 1記憶部に記憶する手段と、 (b)数 nの b進数展 in- J m
開 (η , η , η,) を計算し、 第 2記憶部に順番に記憶する手段 m m一 1 1
と、 (c)第 2記憶部を順番に走査して、 n .≠ b— 1 となる最小の順番 j を検出する検出手段と、 (d)検出手段に応答し、 予め記饯された、 準乱 数の生成行列 Tの第 j列の成分 ( ·ΐ , 1 j B β~ r 1 j ) 及び第
ID 1記憶部 に記憶された準乱数 Ynの成分を順番に取り出し、 同一の順番で取り出 された生成行列 Τの第 j列の成分 j .及び準乱数 Ynの成分 y iを法 bで 加算し、 加算結果を頫番に第 3記憶部に記憶する手段とを必要とする。 これによ り、 行列の掛け算処理が遅いという技術的問題を解決するとと もに、 コンピュータにおいて用いられる準乱数の生成を高速に行う装置 を提供することができる。 また、 k次元の準乱数を生成するには、 (a) n番目 (nは、 1以上の 整数) の準乱数 (u (1) u (k) ) の個々の座標値 u )
π n n
( 1 i ≤ k ) を、 b進数 ( bは、 kよ り大きいか又は等しい素数) 展 開し、 各座標値につき、 生成された m個の成分を、 順番に第 1記憶部に 記憶するステップと、 (b) 数 nの b進数展開 (η , n η,) を
Dl ID— J J 計算し、 第 2記億部に順番に記憶するステップと、 (c) 第 2記憶部を順 番に走査して、 n.≠ b— l なる最小の順番 jを検出する検出ステツ
プと、 (d) 検出ステップに応答して、 予め記憶された、 準乱数の第 i番 目 ( 1≤ i ≤k) の座標の生成行列 T(i)の第 j列の成分及び n番目の 準乱数の第 i番目の座標値の成分を順番に取り出し、 同一の順番で取り 出された生成行列 Τ(ί)の第 j列の成分及び n番目の準乱数の第 i番目 の座標値の成分とを法 bで加算し、 加箅結果を順番に第 3記慷部に記憶 するステップとを含むようにする。 このような k次元の準乱数を用いて、 関数 f の多重楱分を行うには、
(a) k個の要素に関連する関数 f ( xr x„ xk) を規定するステ ップと、 (b) n番目 (nは、 1以上の整数) の準乱数
(u u (k)) の個々の座標値 u "ノ ( 1≤ i ≤k) を、 b進 n n n
数 <bは、 kよ り大きいか又は等しい素数) 展開し、 各座標値につき、 生成された m個の成分を順番に第 1記憶部に記億するステップと、 (c) 数 nの b進数展開 (η , η , η η > を計算し、 第 2記憶部に順番 m m-1 1
に記憶するステップと、 (d) 第 2記憶部を順番に走査して、
n .≠ b— 1 となる最小の順番 jを検出する検出ステップと、 (e) 検出 ステップに応答して、 予め記憶された、 準乱数の第 i番目 ( 1≤ i ≤k) の座標のための生成行列 T ")の第 j列の成分及び n番目の準乱数の第 i番目の座標値の成分を順番に取り出し、 同一の順番で取り出された生 成行列 Τ )の第《5列の成分及び η番目の準乱数の第 i番目の座標値の 成分とを法 bで加算し、 加算結果を順番に第 3記憶部に記憶する ステツ プと、 (f) m個の加算結果が生成されるごとに b進数逆展開し、 生成さ れた u +1 u)を第 4記憶部に記憶するステップと、 (g) ステップ(f)が 座標値 u +1 ) ( 1≤ i ≤ k ) を第 4記憶部に記憶した場合には、 f (u un+1(k)) の値を計算するステップと、 (h) ステツ プ(b)乃至ステップ(g)を繰り返し、 f ( u (1) u (k)) の値を
n n
n = l から n = N (Nは、 1 よ り大きい予定の整数) までアキュムレー タにて累積するステップと、 U) アキュムレータに応答して、 累棲値を Nで除算するステップとを含む。 このようにすると、 積分値が高速に求 められる。 さ らに上述したように、 先の準乱数生成装置を並列処理に応用した場 合には、 複数のプロセッサと、 複数のプロセッサの各々に、 処理開始番 号 sを含む処理割当情報を送信する手段とを有する準乱数生成システム が構成される。
そして、 各プロセッサが、 (a) 処理割当情報を受信し、 処理開始番号 sに対応する準乱数 ( u (1) u„(k)) を生成し、 第 1記憶部に記
s s
悌する手段と、 (b) η番目の準乱数 ( u (1) un (k)) の個々の座 標值 u (i) ( 1 ≤ i ≤k) を、 b進数 < bは、 kよ り大きいか又は等し n
い素数) 展開し、 各座標值にっき、 生成された m個の成分を順番に第 2 記憶部に記憶する手段と、 (c) 数 nの b進数展開 (nm, η , η,) を計算し、 第 3記憶部に順番に記憶する手段と、 (d) 第 3記憶部を頃番 に走査して、 n .≠ b— 1 となる最小の頼番 j を検出する検出手段と、
(e) 検出手段に応答して、 予め記憶された、 準乱数の第 i 番目
( 1 ≤ i ≤ k ) の座標のための生成行列 T( の第 j列の成分及び n番 目の準乱数の第 i番目の座標値の成分を順番に取り出し、 同一の順番で 取り出された生成行列 T( の第 j列の成分及び n番目の準乱数の第 i 番目の座標値の成分とを法 bで加算し、 加算結果を順番に第 4記憶部に 記憶する手段と、 (Π m個の加算結果が生成されるごとに b進数逆展開 し、 生成された u を第 1記億部に記憶する手段と、 (g) 手段(b) 乃至手段(f)を制御して、 s + 1番目の準乱数
( u n (1) u ^ (k)) から当該プロセッサにおける処理終了番号 s+1 s+1
tに対応する準乱数 ( ut U) ut (k)) までを第 1記憶部に記恡さ せる制御手段とを有する。 これによ り、 個々のプロセッサの計算量は、 当然プロセッサの数だけ减少したことになるため、 よ り高速に計算が行 われること となる。 の場合、 各プロセッサに存在する手段 (a) は、 数 nの b進数展開 n n rij) を計算し、 その計算結果によ り生成されるべク
•r
トル V = ( n, n n ) と生成行列 T (i) ( 1 ≤ i ≤ k ) と行 列 G
(i)
とを用いて、 n番目の举乱数の i番目の座標值 u を生成する生成手 段を含むようにすることも考えられる。 当然、 この並列処 ¾を行う準乱数生成システムは、 並列処? J1を行う関 ¾ f の多 31 ¾i分を計^するシステムにも応用可能である。 さ らに、 並列処¾における準乱数生成システムの他の態搽においては 複数のプロセッサと、 複数のプロセッサの各々に、 処 ¾[¾始番号 sを含 む処理割当情報を送信する手段とを有し、 各プロセッサが、 (a) 処理割 当情報を受信する手段と、 (b) 数 nの b進数展^ ( n . η , η , )
m m_ 1 1 を計算し、 その計算結果によ り生成されるベク トル V ( n
(i)
n n ) と生成行列 T ( 1 ≤ i ≤ k ) と行列 G
m一 1 m
とを用いて、 n番目の準乱数の i番目の座標値 u
(1)を生成する生成手
n
段と、 (c) 生成手段が処理開始番号 sに対応する举乱数 (u (1) u (k)) から当該プロセッサにおける処理終了番号 t に対応する準乱数 S
( ut (1),..., ut (k)) を生成するように、 数 n ( s≤ n≤ t ) と生成 行列 T(1ノ ( 1≤ i ≤ k) と行列 Gを、 生成手段に入力する制御手段と を有するようにすることも考えられる。 本形態は、 準乱数が前述の準乱 数生成システムにて計算するのよ りも遅く計算されるこ と となるが、 例 えば後の積分のための計算において被積分関数 f が複雑で評価に時間が かかるようであれば、 並列に処理していることもあり、 十分な計算速度 を得ることができる。 なお、 複数のプロセッサにて同一の個数の準乱数 を計箅される場合には、 処¾ ^始番号 sの送信を省略するこ とも可能と なる場合もある。 当然、 この 95の準乱数生成システムも、 l ¾ f の多 31¾分計^シス テムに応用可能である。
[図面の f¾ '単な説明]
第 1図は、 1次元の^乱数を生成する処¾を示すフローチヤ一卜であ る
第 2図は、 本発明を実施するための装置構成であるコンピュータを示 すブロック図である。
第 3図は、 本発頃を実施するための装置構成を示す一例である。
第 4図は、 k次元の準乱数を生成する処理を示すフローチヤ一卜であ る。
第 5図は、 関数 f の多重積分計算の処理を示すフローチヤ一卜である c 第 6図は、 関数 f の多重積分計算を行う装置の一例を示すブロック図 である。
第 7図は、 複数のプロセッサにおいて並列処理する際のシステム構成 の一例を示すブロック図である。
第 8図は、 閧数 f の多重積分計算を並列処理する際のフローチヤ一ト を示す図である。 これらの図面にて用いられる符合の主なものを説明しておく。
1 0 制御部 1 2, 1 6, 1 8, 20, 32 記憧部
1 6 ポインタ 22 後処理ブロック 30 逆展開部
34 演算部 36 アキュムレータ 40 システム制御部
42, 44. 46 プロセッサ
[発明を実施するための最良の形態]
A :本発明の原理
以下では便宜上、 1次元で説明を行う。 しかし、 後で説明するように これは、 容易に k次元に拡張可能である。 数列を、 X ( n = 0. 1 , 2 bm'] ) で表す。 こ こで、 bは、
n
次元 kよ りも大きいかまたは等しい素数、 また、 mは通常、
1 o gb230に近い整数とする。 定 ¾に り、 数列 X は、 生成行列
T = ( c ..) を用いて、 X. T V と書く ことができる。 このような行 列 Tの選び方については、 上記 P. Bratley他の論文を参照されたい。 こ
m-1.
こで、 ベク トル V は、 整数 nの b進表現 n =n b n2b+n】【こ対
n m
して、 n )を対応させたものである
上記 P. Bratley他の論文では、 X = T v という行列の掛け箅をその
n n
まま行ったので、 計算に非常に時間がかかった。
一方、 本発明によれば、 X = T V で直接行列の掛け算を行う代わり
n n
に、 Y = T G V によって数列 Y が発生される。 数列 Y は、 Y
n n n n n
( n = 0 , 1 , 2, ···, bm- 1 ) であり、 bはやはり次元 kよ りも大 きいかまたは等しい素数である。 こ こで、 行列 Gは、 以下に示すように
(法 bで) 対角線上の要素 i i ( i = 1 m) が全て 1 で、 その 1 の要素の右隣の要素 i j ( j - 1 = i . 1 ≤ i ≤m, 1 ≤ j ≤ m ) 力 一 1 で、 他の要素が全て 0であるような疎行列 (sparse matrix) であ る。 このような行列 Gは正則なので、 Y 自体も^乱数となることはす n
ぐ分かる。
本発明者は、 このと き、 極めて非自明にも、 次の うな定 が成立す るこ とに 1到した。
定理
e . ( j = 1 , 2 , m) を、 長さ mの j番目の単位ベク トルとす る。 そこで、 g = G V とおく と、 g = g ADD. e .が成立する
n n n+l n b j ここで、 ADDbは、 座標毎の bを法とする加算を表し、 j は、 べク 卜 レ V = ( η, , n )において、 n .≠ b— 1 を ¾たす最小の添字と n 1 ra j
する。 検証
1 カ +1に変化するとき、 その b進表現がどのように変化するかを調 ベてみる。 例えば、 b = 3と して、 nが 8から 9に変わるとする。
8 = 2::: 3 + 2であるから、 8の 3進表現は、 (022) である。 よ つて、 ベク トル v8は、 以下のように表される。
2
0 また、 9 3ムによ り、 9の 3進表 ¾は、 ( 1 00) である つて ベク トル v9は、 以下のようにあらわされる
0
V = 0
1 よって、 3 X 3の場合の行列 Gを使用すると、 g。と、 g。とはそれぞ れ、 次のよ うに計算される
方、 8の 3進表現 (022) において、 3— 1 =2でない最小の添 字は、 j = 3である。 そこで、 g
8 ADD
3 e
3を計算してみると
^8 ADD3 S3 9 9
となって、 g
8 A D D
3 e
3が、 g
9に一致するこ とが分かる。 この一 般的な証明は、 任意の bにおいて、 nについての数学的帰納法を使用す るこ とによって行う こ とができる。 このよ うな式 g n= s ADD, e .が成立すると き、 その両辺に左
】
惻から行列 Tを掛ける と、 = G 且つ Y = T G であるので、
Y Al = Y ADD, T e .
】 b
が成立する。 と ころ力'、 T e .はま さに、 行列 Tの第 j列なので、
Y ,は、 丫 に、 bを法と して行列 Tの B« j列を加算するだけで求まる +1
このと き の -も決定する必¾があるが、 ベク トル
において、 n ·≠ t>— 1 を ¾たす最小の添字 j を求める こ とは、 単にべ ク トル の成分を走杏するだけでよ く 、 単な処理である。
因みに、 1 00000個の点の生成を、 I BM社の
R i s c S y s t e m/6000上で Cコンパイラを使って作成したプ ログラムによ り、 上記 P. Bratleyの論文に記述されている方法と、 本発 明の方法の両方で行ってみたところ、 前者は 60秒かかったのに対し、 後者では 7秒しかかからず、 約 8. 5倍の計算速度の向上が見て取れた c 以上述べたアルゴリズムをまとめると、 第 1図のようになる。 すなわ ち、 ( 1 ) 第 n番目の準乱数 Y を入力する (ステップ 1 0 1 0 ) 。
n
{ 2 ) nの b進数展開 ( η η ) を求める (ステップ 1 020 ) 。
( 3 ) n .≠ b— 1 となる最小の添字 j を求める (ステップ 1 0 30 ) 。
( 4 ) 生成行列 Tの第 j列ベク トルと Y を m o d bで加算し、 それ を Y +1とする (ステップ 1 040 ) 。 では、 これらの動作を行うハー ドウェアの説明を行う。
上述の動作をコンピュータ · プログラムで行う場合には、 第 2図のよ うになる。 この構成は、 バス 1 0 1 に、 演算及び入出力制御機能をもつ中央処理 装置 (C P U ) 1 02、 プログラムをロー ド し、 また、 C P U 1 02の ための作業領域を与える主記憶 (R AM) 1 04、 コマン ドや文字列な どをキー入力するためのキーボード 1 06と、 中央処理装置を制御する ためのオペレーティング ' システム、 行列 Tを初期化するための内容を 含むファイル、 本発明に係る金融派生商品の価格づけ計算プログラムな どを格納したハードディスク 1 08と、 データベースの検索結果を表示 するためのディ スプレイ装置 1 1 0と、 ディ スプレイ装置 1 1 0の画面 上の任意の位置をポイント してその位置情報を中央処理装置に伝えるた
めのマウス 1 12を接統した通常の構成である。 この実施例では、 I B M社から販売されている R I SC S y s t e m/6 OOO ( I BMの 商標) が使用される。 しかし、 やはり I BM社から提供されるている PS/2 ( I BMの商標) 、 P SZ55 ( I BMの商標) などのパーソ ナル · コンピュータのシステムで本発明を実施することも可能であって, コンピュータ · プログラムにおいて実施された本発明は、 特定のハード ゥエア構成に限定されるものではない。 オペレーティング - システムと してこの実施例では、 A I X ( I BM の商標) が使用される。 しかし、 P SZ2または P SZ55などのハー ドゥエア構成を使用する場合は、 W i n d ows (マイクロソフ トの商 標) 、 OSノ 2 ( I BMの商標) 、 PC - DOS、 MS-DOS (マイ クロソフ 卜の登録商標) などのオペレーティング - システムを使用する ことも可能であり、 コンピュータ · プログラムにおいて本発明を実施す るために、 特定のオペレーティング ' システム環境に限定されるもので はない。 では、 専用ハードウエアにて実施をする場合を第 2図に示す。
制御部 1 0は、 Y nと数 nと数 bを人力と して受け取る。 この人力 Y n は、 結采記憶部 20からの出力を用いる場合がある。 制御部 1 0は、 準 乱数展開記憶部 12、 数 n展開記憶部 14、 生成行列記憶部 18、 結果 記憶部 20に接続されており、 各記憶部には必要なデータを分けて記憶 する。 但し、 これらの記憶部は、 分けてあっても、 1つの記憶部を分割 して用いてもよい。 最終の出力は結果記憶部に格納され、 出力は後の処 理 (例えば、 関数 f の多重積分) のための後処理ブロック 22に接続さ れている。 1次元における各記憶部の記憶位置数は、 準乱数展開記憶部
1 2、 数 n展開記憧部 1 4、 結果記憶部 2 0が、 数 n及び n番目の準乱 数 Y が展開されるだけの数 mであり、 生成行列記憶部 1 8は、 m X m 個の記憶位置がある。 制御部 2 2は算術演箅を行うことができ、 次のような動作を行う。 最初に、 生成行列記慷部 1 8には、 予め生成行列 Tが記憶されており、 必要に応じて任意の列の要素を取り出せるようになつていることを述べ ておく。 まず、 数 n及び数 nに対する準乱数 Y が入力されると、 制御部 2 2 n
は n及び Y をそれぞれ b進数展開 (bは 2以上の素数) する。 すなわ
n
ち、 Y は y bm— ' · · y9b+y iと展開され、 この y . · * y9 , が準乱 数展開記憶部 1 2に順番に記憶される。 順番は からでも からでもよ いが、 後の処理のため取り出し方は記憶する順番に対応させる。 同様に 数 nも n=n bB 】+ · · · i^b+rhと展開され、 この η , · '!^,^が数!!展
m 2 J m z 1 開記憶部 1 4に順番に記憶される。 そして制御部〗 0は、 数 n展開記億部 1 4を が格納された記憶位置 から走査する。 この走査にはボインタ 1 6が用いられ、 から走査して いって、 ポインタ 1 6が移動するごとにカウンタが 1 インク リ メン ト さ れる。 制御部 1 0は、 各記憶位置において n .≠ b— 1 であるかを検査 し、 それに該当する最小のカウンタ値 j を検出する。 それに応答して、 制御部 1 0は生成行列記憶部 1 8から第 j列目の成分 (j , · · , j9. j , ) を取り出す。 この取り出し方は、 j' からでも j ¾からでもよいが、 j から
m J m であれば、 Y の成分を取り出す際には y から取り出さなければならな
n m
い。 同様に、 から取り出すのであれば、 から取り出さなければなら
ない。 どちらにしても、 準乱数展開記使部 1 2の成分と、 生成行列記 «部 1
8の成分を取り出し、 各々法 bにおいて ^ + , ·
1 1 · · j ra + の計算を m
する。 この計算結果 <y ' , · · y ' . y ' ) を、 結果記馑部
m Δ 1 2 0に順番に 記慷することとなる。 これにて、 γ +1の計算が完了したこととなる。 実際には、 これは b進数展開された Y であるから元に戻す必要もあ るが、 後の計算にも用いるので結果記憶部 2 0では b進数展開された形 式で記慷する。 さらに Y .2を計算したい場合には、 結果記憶部 2 0から今度は Y +1 を Y nと して制御部 1 0に入力することになる。 しかし、 前回は Y nが 進数展開されていないことを前提と していたが、 今回は既に展開されて いるので、 その処理は不要となる。 よって、 最初に制御部 1 0は数 n + 1 を b進数展開するこ とになる。 このこ とを鑑みれば、 準乱数展開 記憶部 1 2と結果記憶部 2 0は同一と し、 使い回すこともできる。 以下, 同様の処理であるから省略する。 この第 2図に示した装置において、 準乱数展開記悌部 1 2に記憶され る Y nの b進数展開されたものは、 生成行列記憶部 1 8に記憶された生 成行列 Tの j列目の成分との加算を行う際に、 適切な成分との加算が必 要であるから、 記憶及び取り出しの順番を整えるため、 F I F 0
( Fi rst I n F i rst Out ) バッファを用いてもよい。 また、 数 n展開記憶部 1 4を走査する際には から走査し、 カウン タをインクリメント していたが、 n から走査して、 カウンタの初期値
Π1
を mとし、 1つづつデクリメントする方法もある。 しかし、 Γ^-≠1)一 1 となる最小の jを求めるわけであるから、 デクリメントする方が時問 がかかる。 さらに、 数 nが人力されるようになっているが、 制御部 10にてカウ ンタを保持し、 1つの準乱数が計算されるごとにカウン卜アップして、 その値を用いるようにすることもできる。 以上の説明は、 Y から Y +1の計算を説明したものであるが、 実際に は が人力され、 そこから処理が開始される。
B :高次元の場合の処理
上記説明をふまえて、 次に k次元の準乱数 (特に、 フ ォール · シーケ ンス) を生成する場合の処理を、 第 4図を参照して説明する。 この第 4 図は、 第 2図に示した装置において、 本発明をコンピュータ · プログラ ムにおいて実施した時の処理を示すものである。 こ こでは、 k次元空間の第 n番目の点 (u n (1) u (k)) が与え
n
られたときに、 第 n + 1番目の点 ( u u +1 (k)) を計算す る処理が説明されるが、 実際的には、 初期的に第 1番目の点
( Uj(1) u k)) の座標値を与え、 その後第 4図の処理を反復的 に適用することによって、 任意の数の点が生成されることになる。 また、 k次元空間での準乱数では、 上述した生成行列 Tが、 個別の座 標毎に個別に用意されるこ とに留意されたい。 こ こでは、 i番目の生成 行列を T u)とあらわすものとする。
第 4図のステップ 202では、 k次元空間の第 n番目の点
( u u (k)) が与えられ、 また、 生成すべき点の座標の番号 n n
をあらわす値 i が 1 にセッ トされる。 ステップ 204では、 所定の素数 b (但し、 bは、 次元の数 kよ りも 大きいかまたは等しい) に従い、 nが b進展開され、 数字 nに対応する m元のベク トル (n ny n )を求める計算が行われる。 ステップ 206では、 ステップ 204で生成された m元のべク 卜ル
(η,,η,, π )の要素を迪つて、 η .≠ b— 1 を満たす最小の j が見出 される。 ステップ 208では、 本発明の原理に従い、 第 ί 座標の生成行列 T ) の第 j列ベク トルと、 u ( 1)を mo d bで加算し、 それによつて
n
u ^(i)が求められる。
n+1 ステップ 2 1 0では、 i = kかどうかが判断される。 もしそうなら、 第 n + 1番目の点の全ての座標が計算されたことになる。 i = kでなければ、 第 n + 1番目の全ての点の座標の計算がまだ完了 していないので、 ステップ 2 1 2で i を 1增分してから、 ステップ 20 8に戻り、 第 n + 1番目の点の次の座標の点の計算を行う。 では、 専用ハー ドウェアにおいて実施する場合には、 どのようにする のかを説明する。
第 3図では 1次元の場合のみを考えたが、 容易に k次元に拡張できる, すなわち、 第 3図における準乱数展開記憶部 1 2が 1 X m個の記 1*位置 を有するものであつたのに対し、 これを m X mにすればよい。 また、 生 成行列記億部 1 8には生成行列丁が 1 つ記億されていたが、 この場合に は k個必要になるので、 m x m x k個の記憶位置を有する記憶部とすれ ぱよい。 また結果記憶部 2 0においても、 生成された準乱数の座標値の 個々の成分を記憶するために、 m X m個の記憶位置を有するように拡張 すればよい。 そして、 制御部 1 0には、 Y の代わりに u ( 1 )
π n
u (k)を入力する。 生成行列制御部 1 8においては、 生成行列 Tを 1 つ 格納する m X mのブロックごとに管理するようにし、 そのブロックを選 択する選択信号を制御部 1 0から生成行列記憶部 1 8に出力するように することもできる。 但し、 制御部 1 0には u u )を 1度に一つだけ入力するようにするな
n
らば、 準乱数展開記憶部 1 2は 1次元の場合と同様であってもよい。 準 乱数展開記悽部 1 2とは異なり、 生成行列記憧部 1 8は、 ί によって生 成行列 Τが異なり、 すべての i に対し生成行列 Τが必要となるので、 前 述の大きさの記憶位置を必要とする。 また、 結果記憶部 2 0は、 結果を 1 つ分のみ記憶する場合には 1次元の場合と同様であってもよいが、 結 果記憶部 2 0から制御部 1 0へのフィードバックのループを用いるので、 k次元分の記憶位置を有することが望ま しい。
C : 多重積分計算への応用
次に、 k次元の準乱数を用して多重積分計算を行う処理について、 第 5図を参照して説明する。 第 5図は、 第 2図に示すコンピュータ · シス テムにおいて、 本発明をコンピュータ · プログラムにて実施した場合の
処理を示したものである。 また、 ここでは f (χΓ x2 ,xk) という k 個の引数をもつ関数を、 その個々の引数が張る k次元空間内で多重積分 するものとする。 第 5図のステップ 302では、 結果の積分値が格納される変数 Sと、 積分計算用に生成された k次元の準乱数の点の数を格納する変数 nとに 0が格納される。 ステップ 304では、 nが 1増分され、 ステップ 306では、 第 4図 に示したステップを利用して、 k次元空間の第 n番目の点
( u u (k)) が生成される。
η π ステップ 308では、 f ( u (1) u (k)) が計算され、 その結
n n
果の値が Sに加算される。 尚、 ステップ 306では、 u (,)
( l ≤i≤ k> の各々は、 b進展開された個々の值を要素とするべク ト ルと して扱われたが、 ステップ 308で f に代人されるとき、
u (i) ( 1 ≤i≤ k) はスカラー値と して扱われることに留意されたい ( n ステップ 3 1 0では、 nが予定の数 Nに達したかどうかが判断され、 もしそうでないなら、 nを 1 だけ增分するステップ 304の処理に戻る ( 尚、 Nの値は、 必要な計算の精度に応じて決定され、 Nの値を大き く と る程精度は向上するが、 その分計算時問が長く かかるこ とになる。 通常 のモンテカルロ法のよ うに、 線形合同法によつて乱数を生成した場合、 生成される点のばらつきの度合があま り大きく なく、 Nの値を相当に大 きく と らないと十分な精度が得られない。 準乱数はこの点で改善された 性質を示し、 生成される点は、 線形合同法によって生成された点よ りも
相当によいばらつきの度合を与え、 結果的に、 よ り少ない数の Nで十分 な精度を達成することを可能ならしめる。 ステップ 3 1 0で nが予定の数 Nに達したと判断されたなら、 ステツ プ 3 1 2で Sの値を Nで割って正規化した値が Sに格納され、 結局、 こ の Sが積分值をあらわすことになる。 では、 第 3図に戻って多重積分を行うための後処理ブロック 22の説 明を行う。 この後処理ブロック 22の詳細を第 6図に示す。 第 6図にお いて、 結果記憶部 20 'は k次元分の記憶位置を有するように拡張され ている。 この結采記憶部 20'は、 数 bを他の入力と して有する逆展開 部 30に接続されている。 この逆展開部 30は、 出力をスカラー値記憶 部 32に記憶するようになっている。 スカラー値記憧部 32には 1の準 乱数の k個の座標值を記億できるよう、 k個の記憶位置を有する。 そし て演算部 34は、 スカラー値記憶部 32に接続され、 関数 f と準乱数の 生成個数 Nも入力とされる。 演算部 34は、 アキュムレータ 36に接続 され、 アキュムレータ 36の出力は演算部 34に戻される。 次に第 6図の動作を示す。 結果記憶部 20 '以前の部分は第 3図にお いて説明したので、 結果記憶部 20 'には、 既に un+1 (1) u +] (k) のそれぞれにっき b進数展開されたものが記億されているものと して以 下の動作を説明する。 逆展開部 30は、 un+1 U),..., u ./k)のそれ ぞれにっき b進数展開されたものを結果記憶部 20 'から読み出し、 そ れぞれ数 bにて逆展開することによ り un+1 (1) un+1(k)を生成す る。 k次元であれば、 k個の un+1 (1)が生成されたと ころで、 演算部 3
4は u u Al (k)を取り出す。
n+1 n+1
その前に、 演算部 34には被積分関数 f (xr x2 xk) が入力され. 演算部 34にセッ トされている。 よって、 この変数 Xiに un+1 (1)を入力 すると、 関数値 f ( u un+1 (k)) が計算される。 この関数 值は、 アキュムレータ 36に出力される。 演算部 34は、 次々に生成さ れる準乱数に対し関数値を計算し、 アキュムレータ 36はその関数値を 加算、 蓄檳されていく。 アキュムレータ 36は、 カウンタを有しており 演算部 34から N回入力されたところで、 蓄積された値を演算部 34に 戻す。 演算部 34でも Nは人力されているので、 この Nによ りアキュム レー夕 36からの人力を除算する。 これによ り、 関数
f (x n (1 ) X n (k)) の k重檳分が行われたこととなる。 上記部分においては説明を省略したか、 アキュムレータ 36のカウン タは、 Nを示したところでリセッ トされ、 また累積していた値について も リセッ トする。
D : 金融派生商品の価格づけ計算への応用
多 稻分計箅を必 ¾とする典型的な応用分野と して、 金融派生商品の 価格づけ計 の例について説明する。 以下に示すのは、 抵当 された 証券の ¾在 E ( P V) の計算式である。
E(P V) = J P V( i , w, a ) f ( i , w, a ) d i dwda ここで、 P Vは次のように与えられる。
ί - v w7 ak Π - w Π— νゥ)
P V = C— ÷ C- ÷ C——— 7— -斗 …
上記式において、 Cは、 抵当上での定期的な支払いであり、 i tは、 月 tでの利子率であり、 w, tは、 月 tでの予納率であり、 a .は、 残り
χ-ι
の年金である。
実際的な計算では、 かなり長期間に渡って計算が行われ、 すると、 i .の数が大きく なるので、 一般的には 360次元程度の多重橫分が必 要となり、 これは、 従来の方法では、 何日もの計算を要するが、 本発明 の技法を用いることによ り、 この計算を数時間以内で完了し得るように なる。 以上説明してきた方法及び装置を用いると、 高速に準乱数が生成でき るので、 関数 f の多重積分の計算も高速になった。 では、 この計算を複 数のプロセッサにおいて実行することを考える。
E :並列処理について
上述では、 Y = Y ADD, T e .と Y = T G v とが等価である
n+1 n b j π π ことを示した。 すなわち、 任意の ηに対して Υ ηを求めるには、 生成行 列 Τと行列 G (先に示した) と数 ηを b進数 (bは 2以上の素数) した ものをベク トルと した V の行列の掛け算にて計算できる。 任意の数 n π
に対応できるのは、 行列 Gの性質による。 そして、 Y がー度生成され れば、 Y は上述した Y =丫 ADD. T e .にて行列の掛け算な
n+1 n+1 n b j
しの高速な方法にて計算される。 任意の nに対して Y を計算できるの で、 Y Al= Y ADD. T e .にて必要な nまで 1から空回し計算を行 n+】 n b j
う必要がない。 このことが重要である。 このこ とを念頭に置きながら第 7図及び第 8図を参照する。
第 7図は、 複数のプロセッサ (42, 44, 46) と、 システム制御
部 4 0がバス 4 8にて接統されているシステムを表すものである。 この システム制御部 4 0は、 第 1 プロセッサ 4 2等のプロセッサとは他に設 けてもよいし、 複数のプロセッサ 《4 2 , 4 4 , 4 6 ) のうち 1 つに含 まれるものでもよい。 また、 システム制御部 4 0は、 第 1 プロセッサ等 と同様のものでも、 専用のものでもよい。 このシステムの動作を第 8図を用いて説明する。
システム制御部 4 0は、 このシステムにおいて処理に用いることので きるプロセッサの数及びそのプロセッサの I Dを認識し、 どのプロセッ ザには全準乱数のどの部分を割り当てるかを計算する。 システム全体の ジョブに偏りがなければ、 各プロセッサに均等に分割すればよい。 但し. 一般的には処理開始番号 s及び処理終了番号 t を処理割当情報と して与 える。 処 S終了番号 tはそのプロセッサにおける準乱数生成個数に Sき 換えることができる。 また、 常に各プロセッサに均等に処理が割り当て られるのであれば、 準乱数生成個数又は処理終了番号 tは各プロセッサ にいちいち知らせることもないので、 それらの送信は省略することもあ る(ステップ 6 0 ) 。 各プロセッサでは、 各々異なる処理開始番号 sを受信し、
Y = T G V から Y をそれぞれ生成する(ステップ 6 2 ) 。 このよ うに
Π Π S
任意の sに対応する Y を生成できるのは、 先ほど述べたように行列 G の性質である。 一旦 Y sが生成できると、 Y n ++11 = Y n A D D . b T e j.に従い、 受信し た処理終了番号 tに対応する Y tまでの準乱数を、 生成する(ステップ 6 4 ) 。 この処理の詳細は、 上述で行っているので省略する。 これは、 第
4図に示した手頤や、 第 3図 (拡張したものも含む) で示した装置にて 実施することができる。 さ らに、 本発明において多重穰分を行う場合に は、 被積分関数 f も評価しなければならないが、 この被積分鬨数 f の評 価に時間がかかる場合には、 Y A1 = Y ADD. T e .にて計箅する方
n+] n b j
法を用いなく とも、 いちいち行列の掛け算を行う Y n = T G V nに従って 計算しても全体の実行時間は実質的に同じである場合もあるので、 ステ ップ 64は、 どちらの方法を用いてもよい。 但し、 準乱数の生成のみの 観点では Y n ++1l = Y n ADD. b T e j .に従って計算する方が高速であるこ とは明らかである。 各プロセッサにて割り当てられた準乱数 Y から Y,までが生成された
s t
のち、 又は各準乱数が生成されるたびに、 システム制御部 40から送ら れてきた関数 f に代入する。 そして、 計算された関数値 f ( Y n )
( s≤ n≤ t ) を全て加算し、 S . ( iはプロセッサの番号、 プロセッ サは h個ある。 ) を計算する (ステップ 66) 。 この計算値 S. ( 1≤ i≤ h ) は、 すべてのプロセッサからシステム 制御部 40に送信される (ステップ 68) 。 システム制御部 40では、 送付された S.を全て加算し、 全準乱数生成個数 Nで除算する (ステツ プ 70) 。 これにて、 関数 f の多重積分計算は完了した。 上述の説明では、 生成した準乱数に対する関数値 f ( Y ) を計算し た後、 その加算を各プロセッサにて行ったが、 加算する前にシステム制 御部 40に送信してもよい。 但し、 パス 48のトラフィ ックや各プロセ ッサにおける送信処理等のため他の処理が遅く なる場合もある。 また、 各プロセッサにて準乱数を生成したと ころで、 システム制御部 40に送
信することも可能である。 この場合には、 関数 f を各プロセッサに送信 する必要は無く なる。 このようにした時の問題は、 関数値の加箅前に送 信するのと同様である。 第 3図及び第 6図で示した装置を各プロセッサに含ませる場合には、 第 6図におけるアキュムレータ 36からの出力を演算部 34に戻さず、 アキュムレータ 36の出力をそのプロセッサの出力とする変更を必要と する。 また、 アキュムレータ 36には準乱数の生成個数 Nの入力を必要 としているが、 この Nは、 1にしておきシステム制御部 40においてま とめて除算する。 各プロセッサにおいて計算する準乱数の数が同じなら ば、 Nをそのプロセッサにおいて処理する準乱数の個数に置き換えても よい。 以上、 並列処理の説明においては簡単のため Y をもって説^したが、 これは (u n u n Al (k)) とするよ うに拡張することは前述の
+1 +1
ようにできるので、 それによ り多重橫分処理を行うとの意味である。 なお、 第 8図のステップ 64においては Y +1 = Y ADD, T e .又 は Y =TG ν に従って計箅すると した力1、 各プロセッサはどちらに従 n n
つても計算できるようにしておき、 システム制御部 40が関数 f の性質 を検知して、 どちらの方法を各プロセッサが実施すべきかを示すスィッ チ情報を送信するようにしてもよい。 また、 LANなどで接続された分散処理環境で複数台の処理装置を接 統し、 各々の処理装置のローカル ' メモリを分散して使用することによ り、 巨大な次元の処理に対処したりすることも上述の簡単な変形にて実
施できる。 以上述べた、 コンピュータ ' プログラムによ り実行する形態をとる場 合には、 各ステップに対応するプログラム · コードにて本発明を実現し, それを含む、 フロッピーディスクや C D— ROM、 MOディスク等の記 憶媒体にて流通させることは、 当業者が容易に実施可能な事項である。
[産業上の利用可能性]
以上説明したように、 この発明は、 準乱数の生成において、 行列の 1 つの列とのモジュロの和を行うという技法を提供することによ り、 準乱 数の生成速度を極めて向上させる、 という顕著な効果を与えるものであ る。 これによつて、 準乱数を利用するさまざまな数値計算 (典型的には 多重積分) の計算速度が実質的に改善される。 さ らに、 並列処理を行うことによ りよ り高速に準乱数が計算され、 そ れによって多重積分の計算等も高速に行えるようになつた。 並列処理の 場合には、 全体の積分計算時間を考えた場合に、 行列の掛け算でも十分 な計算速度を得るようにすることもできた。 尚、 上記実施例では、 金融派生商品の価格づけ計算への応用例につい て説明したが、 計算物理、 流体力学、 コンピュータ · グラフィ ックスな どの広範な分野に本発明を適用することが可能である。 このような準乱 数を使用する計算技法の背景については、 例えば、 丄 F. Traub,
H. ozniakowski, "Breaking Intractabi 1 ity", SCIENTIFIC AMERICAN
January 1994などを参照されたい。