モンテカルロ法（HPC・科学計算）

モンテカルロ法とは何を解いているのか

モンテカルロ法は、確率的サンプリングによって決定論的な問題（積分・期待値・偏微分方程式の解など）を近似する数値計算手法です。核心は単純で、求めたい量を「ある確率分布のもとでの期待値」として書き直し、その分布から乱数でサンプルを取り、標本平均で期待値を推定するという発想です。

例えば区間 [0,1] 上の積分 I = ∫f(x)dx は、一様分布からのサンプル x_1, ..., x_N を使い、次の標本平均で推定できます。

モンテカルロ積分の基本形

  I_hat = (1/N) * Σ f(x_i),  x_i ~ Uniform(0,1)

  大数の法則により N → 無限大 で I_hat → I に概収束

この定式化が強力なのは、被積分関数の次元が増えても推定量の構造が変わらない点です。台形則やシンプソン則のような格子求積は、次元 d が上がると格子点数が n^d で爆発します（次元の呪い）。モンテカルロ法はサンプル点を格子ではなく分布から無作為に取るため、次元に依存しない収束オーダーを維持できます。

収束速度O(1/√N)の意味と限界

モンテカルロ推定量の標準誤差は、中心極限定理から次のように導かれます。

標準誤差の導出

  Var(I_hat) = Var(f(x)) / N = σ²/N
  標準誤差   = σ/√N

  → 誤差を1桁減らすには N を100倍にする必要がある

この O(1/√N) という収束速度は、次元 d によらず一定です。対照的に、格子求積の誤差は次元 d に対して O(n^(-r/d))（r は求積公式の次数）のように悪化します。次元がある閾値を超えると、モンテカルロ法の O(1/√N) の方が格子求積より速く収束するようになります。これが高次元問題でモンテカルロ法が選ばれる根拠です。

一方で欠点も明確です。次元によらず遅い。桁数を1つ稼ぐのに100倍のサンプル数が要るため、高精度が必要な低次元問題では他の数値解析手法（求積・スペクトル法）に劣ります。したがって「モンテカルロ法を使うべきか」の判断は、次元の高さと要求精度のトレードオフで決まります。

分散低減技術 ── 標準誤差の定数項を削る

O(1/√N) というオーダー自体は変えられませんが、標準誤差の係数である分散 σ² を小さくすることで、同じサンプル数でも実用精度に到達できます。これが分散低減技術の狙いです。

重要度サンプリング

被積分関数 f(x) の値が大きい領域からサンプルを優先的に取り、サンプリング分布と真の分布のズレを重みで補正する手法です。

重要度サンプリング

  I = ∫ f(x) p(x) dx = ∫ [f(x)p(x)/q(x)] q(x) dx
  → q(x) からサンプリングし、重み w(x) = p(x)/q(x) を掛けて平均

  I_hat = (1/N) * Σ f(x_i) w(x_i),  x_i ~ q(x)

理想的な q(x) は f(x)p(x) に比例する分布で、この場合分散はゼロに近づきます。ただし q(x) の裾が真の分布より薄いと、ごく少数のサンプルが極端に大きな重みを持ち、分散が逆に爆発する危険があります（重みの分散が発散するケース）。輸送計算では、吸収断面積が大きく粒子が滅多に届かない領域の寄与を拾うために重要度サンプリングが多用されますが、重み関数の設計を誤ると少数サンプルに支配された不安定な推定になる点に注意が必要です。

制御変量法

真の期待値が既知の関数 g(x) を使い、f(x) との相関を利用して分散を削る手法です。

制御変量法

  期待値が既知: E[g(x)] = μ_g （解析的に分かっている）

  補正推定量: I_hat = (1/N)Σ f(x_i) - c * [(1/N)Σ g(x_i) - μ_g]

  最適な c = Cov(f,g) / Var(g)  のとき分散は最小化され
  Var(I_hat) = (1 - ρ²) * Var(f)/N   （ρ は f と g の相関係数）

f と g の相関 ρ が1に近いほど分散削減効果が大きく、ρ = 1 なら理論上分散はゼロになります。実務では、非線形問題を線形近似した解や、簡略化した解析解を制御変量として使うのが典型です。重要度サンプリングが「サンプリング分布を変える」のに対し、制御変量法は「サンプリング後に既知の情報で補正する」という違いがあり、両者は併用可能です。

並列疑似乱数生成のストリーム分割

モンテカルロ法はサンプル同士が独立であるため、サンプル生成を並列プロセス・スレッドに分割しやすいという特性があります。しかし、疑似乱数生成器（PRNG）は内部状態を持つ逐次的なアルゴリズムであるため、単純に各並列ワーカーが同じ生成器を異なるシードで初期化すると、周期の重複や系列間の相関という問題が生じ得ます。これを避けるための設計が並列PRNGのストリーム分割です。

代表的な方式は次の2つです。

方式1: ブロックスキップ（block-splitting）

  周期 P の生成器を、ワーカー数 M 個に分割
  ワーカー k の担当区間 = [k * (P/M), (k+1) * (P/M))
  → 各ワーカーは自分の開始点までジャンプ（skip-ahead）してから逐次生成

方式2: パラメータ化（異なるジャンプ多項式 / 異なるパラメータ集合）

  各ワーカーに、周期が重ならないことが数学的に保証された
  異なるパラメータのPRNG（例: 異なる素数を使う多重再帰生成器）を割り当てる

ブロックスキップ方式では、周期 P が非常に大きい生成器（例えばメルセンヌ・ツイスタは周期 2^19937 - 1）を使うことで、ワーカー数を大きくしても各区間が十分長く保たれます。skip-ahead操作は、生成器の状態遷移を行列べき乗のように表現し、対数時間で任意ステップ先へジャンプできるアルゴリズム（多くの生成器で実装されている）を用います。

パラメータ化方式は、各ワーカーが独立した数学的性質を持つ生成器インスタンスを使うため相関のリスクを構造的に排除できますが、パラメータ集合を安全に用意する事前検証コストがかかります。

まとめ

• モンテカルロ法は期待値としての積分を標本平均で推定する手法で、誤差は次元によらず O(1/√N)。次元が高い問題では格子求積より有利になる。
• 収束オーダー自体は変えられないため、分散 σ² を減らす重要度サンプリングと制御変量法が実務上の精度向上の主軸になる。
• 重要度サンプリングは提案分布の裾の設計を誤ると分散が悪化し、制御変量法は代理関数との相関の高さが効果を左右する。
• 並列実行では、乱数生成器のブロックスキップやパラメータ化によるストリーム分割で、ワーカー間の周期重複・相関を避ける設計が精度保証の前提になる。
• 並列計算パターン全般の基礎は/programming/、数値計算基盤としてのプロセッサ特性は/semiconductor/も参照。