モンテカルロ法のはじまり

「案ずるより産むが易し」とは、事前に思い悩んでいたが実際に行動してみると意外に簡単であった、といった意味ですが、アメリカの数学者スタニスワフ・ウラムはまさにこの発想でモンテカルロ法を考案しました。

1946年のある日、病床でソリティアをして暇つぶししていた数学者ウラムは、「ランダムに配られたソリティアが最後まで完成する確率はどのくらいなのだろうか?」といういかにも数学者らしい疑問を抱きました。しかしいくら考えても解けず、しまいには面倒になり「こんな難しい問題を解くくらいなら100回トライして成功確率を見た方が早いんじゃないか…」と思うに至りました。ここで「何回もトライしてみる」というアイデアに対しウラムは思うところがあったようで、これを友人のジョン・フォン・ノイマンに話しました。同じくアイデアの重要性に気付いたノイマンはウラムと一緒になり、このアイデアを形にすべく活動を開始しました。形になったアイデアはノイマンにより「モンテカルロ法」と命名されました [1]。モンテカルロとはカジノで有名なモナコのモンテカルロ地区 (図1) のことですが、ギャンブルという確率的試行の世界で踊らされて借金を重ねるウラムの叔父を重ね合わせてコードネームとして使っていたものが次第に定着したとも言われています。

図1 モンテカルロの街並み
ref: Wikipedia「モンテカルロ」の項

当時はちょうどデジタル計算機がようやく実用化の軌道に乗り始めた頃で、モンテカルロ法は計算機上で数値計算されることを想定しての研究でした。繰り返し数値実験を行うモンテカルロ法がデジタル計算機と非常に相性が良かったのです。やがてモンテカルロ法の数学的基礎と計算機による実施手法を確立したノイマンやウラムらは、当時推進されていた次世代核兵器の開発に向け、中性子が物質中を動き回る様子を探る計算にモンテカルロ法を適用し、大きな成果を挙げることとなります。また、モンテカルロ法は計算機の発展との相乗効果により急速に普及し、自然科学分野に留まらず金融分野等の他分野にも広まってゆくことになります。

モンテカルロ法を支える擬似乱数

順調に開発されたかのようにさらりとモンテカルロ法のはじまりを述べましたが、モンテカルロ法の実用化には解決すべき致命的な問題点が1つありました。ランダム性です。「神はサイコロを振らない」かどうかはアインシュタイン以降今日に至るまで物理学界で議論され続けているようですが [2]、正確性が取り柄のデジタル計算機がサイコロを振れないことは誰の目にも明らかでした。モンテカルロ法はランダムに試行を重ねた結果を統計処理するため、入力データのランダム性が崩れれば出力された結果の統計も意味を持たなくなってしまいます。数値計算におけるランダム性とはランダムな数字、すなわち「乱数」に帰着しますが、乱数を上手く生成することはモンテカルロ法の根底を支える非常に重要な技術と言えます。

ノイマンはこの乱数の生成問題を「擬似乱数生成器」により解決しました。擬似乱数とは一見乱数列のように見えるが実際は確定的な計算によって求める数列 [3] により生成された乱数のことで、擬似乱数列を生成するアルゴリズムさえあれば、サイコロを振れない計算機でもランダムな数字を生成することができます。当時ノイマンが考案したアルゴリズムは「平方採中法」と呼ばれています。例えば4桁の擬似乱数の場合、ある数字を2乗して8桁の数字を生成し、その中央4桁を取って次の数字とします。最初の数字を6239とすると、 \[ 6239^2 = 38925121 \rightarrow 9251 \\ 9251^2 = 85581001 \rightarrow 5810 \\ 5810^2 = 33756100 \rightarrow 7561 \] といったように、4桁の数字が次々と生成されます。モンテカルロ法研究の傍らで考案された手法にしては非常によくできています。ここでは最初の数字として6239を与えましたが、この数字を変えると別の乱数列が得られますし、同じ数字を指定すると同じ乱数列がいつでも得られます。このような乱数列を与えるための初期情報を「乱数の種」と言います。

これにより、コンピュータもサイコロを振ることが可能になりました。擬似乱数もモンテカルロ法と同様、長年にわたって研究がなされ、古典的な定番アルゴリズムである線形合同法に加え、近年のより高品質なアルゴリズムとしてMersenne-Twister法やXorshiftなど、数多くの擬似乱数が開発されていくことになります。

整然とランダムに!? ～ 準乱数による準モンテカルロ法

モンテカルロ法にランダム性が必須であると述べましたが、対象とする問題 (例えば数値積分など) によってはただ単純に乱数を適用すればよいわけでもないことが後年の研究により分かってきました。これは使用している乱数の品質を高めればいいというものではなく、もっと根本的なところに問題があります。例えば図2の左側の図は擬似乱数をMersenne-Twister法により発生させて作った2次元分布ですが、場所によって密度の薄いところと濃いところがあります。ランダムである以上密度の濃淡が発生するのは当然で、逆にそれがランダム性の証左でもあるのですが、この濃淡が収束性を落とす可能性があるというわけです。

そこでこの濃淡を限りなく小さくしながらもランダムに見える数列として準乱数 (LDS: Low-discrepancy sequence、低食い違い数列とも呼ばれる) が開発されました。下右図は準乱数の1つであるHalton列を用いた2次元分布ですが、分布の一様性を保持しながら全体の濃度が濃くなっていく様子が見て取れるかと思います。擬似乱数列の代わりに準乱数を利用するモンテカルロ法のことを特に「準モンテカルロ法」と呼びます。例えば数値積分をモンテカルロ法で計算する場合、擬似乱数を使ったモンテカルロ法ではいくら高品質な乱数発生器を使用したとしても精度を1桁改善するのに100倍程度のサンプリング点数が必要となりますが、準モンテカルロ法では10倍程度のサンプリング点数で済むことが知られており [4]、準モンテカルロ法の威力が分かります。

図2 擬似乱数Mersenne-Twisterと準乱数Halton列を使った2次元分布の比較

準乱数としてはHalton列、Sobol'列 (Sobolの後ろのアポストロフィー「'」はミスタイプではなく「Sobol'列」というのが正式な名前です) などが有名ですが、ここでは最も単純な準乱数の1つであるvan der Corput列について具体的に紹介します [5]。van der Corput列は \(1, 2, 3, \cdots\) という数列のそれぞれについて、「\(n\) 進数に変換し、小数点で折り返す」という処理を行うことによって得られます。例えば \(n=2\) のとき、数字の \(6\) は2進数で \(110.0_{(2)}\) ですが、小数点で折り返すと \(0.011_{(2)}\) となり、これを10進数にもどして \(0.375\) と変換できます。この処理を \(1\) から順番に実行すると、 \[ \begin{array}{ccrclcl} 1 &=& 1.0_{(2)} &\rightarrow& 0.1_{(2)} &=& 0.5 \\ 2 &=& 10.0_{(2)} &\rightarrow& 0.01_{(2)} &=& 0.25 \\ 3 &=& 11.0_{(2)} &\rightarrow& 0.11_{(2)} &=& 0.75 \\ 4 &=& 100.0_{(2)} &\rightarrow& 0.001_{(2)} &=& 0.125 \\ 5 &=& 101.0_{(2)} &\rightarrow& 0.101_{(2)} &=& 0.625 \\ 6 &=& 110.0_{(2)} &\rightarrow& 0.011_{(2)} &=& 0.375 \\ 7 &=& 111.0_{(2)} &\rightarrow& 0.111_{(2)} &=& 0.875 \\ \end{array} \] と変換していくことができ、結果として \[ 0.5, \ 0.25, \ 0.75, \ 0.125, \ 0.625, \ 0.375, \ 0.875, \ \cdots \] という数列が得られます。これがvan der Corput列となります。0から1の間をランダムに数字が発生して行っているように見えて、実はとても規則的です。上記は敢えて7つで止めましたが、この7つを分数として表記すると \[ \frac{4}{8}, \ \frac{2}{8}, \ \frac{6}{8}, \ \frac{1}{8}, \ \frac{5}{8}, \ \frac{3}{8}, \ \frac{7}{8}, \ \cdots \] となります。\(1/8\) ～ \(7/8\) が1回ずつ現れていることとなり、0～1までの間を理想的に等間隔で埋めていることが分かります。このようにして密度の濃淡なくしかもランダムに数を発生させるという相反する要求を実現しています。なお、van der Corput列の \(n\) は素数であれば何でもよく、そのため複数の \(n\) に対して多次元的に準乱数を発生させることも可能ですが、こうして発生させた多次元準乱数がHalton列になります。なお、上右図は横軸が \(n=2\)、縦軸が \(n=3\) としたvan der Corput列を組み合わせたHalton列です。

IsightにおけるMonte Carloコンポーネント

CTCで取扱いをしている自動化・統合化・最適化システムソフトウェアIsightでは、汎用的にモンテカルロ法の計算ができるようになっているMonte Carloコンポーネントが搭載されております。モンテカルロ法の擬似乱数発生器として線形合同法を、準モンテカルロ法の準乱数としてSobol'列をサポートしており、ドロップリストから簡単に選択して実行することができます。また一様乱数以外の確率分布も発生させられるほか、出力結果の整理のための統計的な処理機能も付属しています。

どうしても解けず袋小路に入ってしまった難しい問題は、ウラムのように発想を転換させ、Isightにサイコロを振らせてみると意外とあっさり解決するかもしれませんよ。