はじめに
本記事では、モンテカルロ法の改良版である準モンテカルロ法について解説します。モンテカルロ法は乱数を用いて数値計算や積分を行う手法ですが、準モンテカルロ法は「低不一致列」と呼ばれる特殊な数列を用いることで、より高速な収束を実現します。本記事では、準モンテカルロ法で使用される数列（疑似乱数）の性質について、理論的背景から実装、性能評価まで、具体例を交えながら詳しく解説します。
対象読者:
- 数値積分や数値計算の効率化に関心のある方
- モンテカルロ法の基礎知識があり、さらに理解を深めたい方
- 準モンテカルロ法の実装に興味がある方
記事のポイント:
- 準モンテカルロ法で用いられる低不一致列の概念と、その利点
- 代表的な低不一致列（Halton列、Sobol列など）の特性と実装方法
- 不一致度（Discrepancy）の概念と、それが数値積分の精度に与える影響
- 準モンテカルロ法の性能と次元の関係性、および高次元問題への対処法
サンプル点の選び方：乱数vs格子点vs低不一致列
数値積分や数値計算において、計算に使う点（サンプル点）をどのように選ぶかは、結果の精度を左右する重要な要素です。ここでは、代表的な3つの方法を比較します。
| 方法 | 特徴 |
|---|---|
| 乱数（モンテカルロ法） | - ランダムに点を選ぶ - 確率的な意味で、どの領域にも均等に点が分布 - 点と点の間に偏りが生じる可能性がある |
| 格子点 | - 等間隔に点を配置 - 規則的で美しい配置 - 特定の方向に点が並んでしまう - 次元が増えると、必要な点の数が爆発的に増加する（次元の呪い） |
| 低不一致列（準モンテカルロ法） | - 乱数と格子点の良いとこ取り - 規則的すぎず、ランダムすぎない、絶妙なバランス - あらゆる方向で均一な分布を実現 |
格子点による方法は、低次元では効率的ですが、次元が増えると必要な点の数が指数関数的に増大する「次元の呪い」が問題となります。乱数を用いるモンテカルロ法は、次元の呪いには強いものの、点の偏りによって収束が遅いという欠点があります。
低不一致列は、これらの問題を解決するために考案されました。各点は、それまでに配置された点との位置関係を考慮して配置され、空間全体を均一に埋めるように工夫されています。これにより、格子点法の規則性に起因する問題を避けつつ、乱数よりも効率的なサンプリングが可能になります。
乱数の一様性と不一致性
「良い」サンプル点とは何でしょうか。ここで重要になるのが「一様性」と「不一致性」です。
一様性は、直感的には「点がどの領域にも同じくらいの密度で分布している」ことを意味します。しかし、数値計算においては、単なる一様性だけでは不十分で、「不一致度（Discrepancy）」がより重要な指標となります。不一致度は、点列が空間をどれだけ均一に埋めているかを定量的に評価する尺度です。
不一致度は以下の式で定義されます：
D_N^*(x_1,\ldots,x_N) = \sup_{B \in \mathcal{B}} \left|\frac{1}{N}\sum_{n=1}^N \mathbf{1}_B(x_n) - \lambda(B)\right|この式は、以下のような意味を持ちます：
- 空間内の任意の部分領域
Bを考えます。 - その領域に含まれる点の割合（
\frac{1}{N}\sum_{n=1}^N \mathbf{1}_B(x_n)）を計算します。 - その領域の体積の割合（
\lambda(B)）と比較します。 - すべての可能な領域
Bについて、この差の最大値を求めます。
つまり、不一致度が小さいほど、点列は空間をより均一に埋めていることを意味します。「局所的に見ても大域的に見ても、点の密度が一様である」ことを保証する指標と言えます。
低不一致列の種類と特徴
準モンテカルロ法で用いられる代表的な低不一致列には、以下のようなものがあります。
- Halton列
- Sobol列
- Faure列
- Niederreiter列
これらの数列は、それぞれ異なる方法で生成されますが、共通して低い不一致度を実現します。ここでは、実装が比較的簡単でよく利用されるHalton列とSobol列について詳しく見ていきます。
Halton列
Halton列は、最も基本的な低不一致列の一つです。異なる素数を基底として用い、各次元で van der Corput 列を構成します。n番目のHalton点の第k次元成分は以下のように計算されます
\phi_{p_k}(n) = \sum_{j=0}^{\infty} a_j(n) p_k^{-j-1}ここで：
p_kはk番目の素数a_j(n)はnのp進展開の係数
以下に、2次元Halton列を生成するPythonコードを示します。
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
def vdc(n, base):
"""van der Corput列の生成"""
vdc, denom = 0, 1
while n:
denom *= base
n, remainder = divmod(n, base)
vdc += remainder / denom
return vdc
def halton_sequence(n_points, dim):
"""Halton列の生成"""
bases = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29] # 最初の10個の素数
result = np.zeros((n_points, dim))
for d in range(dim):
for i in range(n_points):
result[i, d] = vdc(i + 1, bases[d])
return result
# 2次元Halton列の生成と可視化
n_points = 1000
points = halton_sequence(n_points, 2)Halton列は、
- 構成が単純で実装が容易
- 低次元（10次元程度まで）で特に効果的
- 高次元になると、基底の素数が大きくなり、品質が低下する可能性がある
といった特徴があります。
Sobol列
Sobol列は、2を基底とする数列で、特に高次元での性能が良いとされています。各次元は以下の漸化式で生成されます。
x_n^{(i)} = x_{n-1}^{(i)} \oplus (v_c^{(i)} \cdot 2^{-c})ここで：
\oplusはビット単位のXOR演算v_c^{(i)}は方向数（事前に定められた値）cはnの最下位ビットの位置
以下にSobol列を生成するPythonコードを示します。
from scipy.stats import qmc
# Sobol列の生成
sampler = qmc.Sobol(d=2, scramble=False)
points = sampler.random(n=1000)Sobol列は、
- 高次元でも良好な性能を維持
- 方向数の選択が重要
- 2のべき乗個の点で特に良い性質を示す
といった特徴があります。
ここで、Halton列、Sobol列を含む、2次元上でのサンプリングの様子を可視化してみます。
不一致度の理論的解析
低不一致列の性能を理論的に裏付けるのが、Koksma-Hlawka不等式です。
\left|\int_{[0,1]^s} f(\mathbf{x})d\mathbf{x} - \frac{1}{N}\sum_{n=1}^N f(\mathbf{x}_n)\right| \leq V(f)D_N^*(\mathbf{x}_1,\ldots,\mathbf{x}_N)ここで：
V(f)は関数fの全変動（関数の滑らかさを表す指標）D_N^*は数列の不一致度
この不等式は、数値積分の誤差が、関数の変動と点列の不一致度の積で上から抑えられることを示しています。つまり、関数が滑らかであるほど、そして点列の不一致度が小さいほど、数値積分の結果は正確になります。
純粋な乱数列の不一致度は、確率的に以下のオーダーになります。
D_N^* = O\left(\sqrt{\frac{\log\log N}{N}}\right)一方、低不一致列では：
D_N^* = O\left(\frac{(\log N)^s}{N}\right)となり、次元sに依存するものの、N（サンプル数）に関して、より速い収束を示します。
実験：積分計算での収束性比較
理論的な結果を確かめるため、以下の2次元積分を例に、モンテカルロ法と準モンテカルロ法で数値計算を行い、収束性を比較します。
I = \int_0^1\int_0^1 \sin(\pi x)\cos(2\pi y)dxdyこの積分の厳密解は0です。
def integrand(x, y):
return np.sin(np.pi * x) * np.cos(2 * np.pi * y)
# サンプル数のリスト
ns = [100, 1000, 10000, 100000]
mc_errors = []
qmc_errors = []
for n in ns:
# モンテカルロ法
mc_points = np.random.uniform(0, 1, (n, 2))
mc_result = np.mean([integrand(x, y) for x, y in mc_points])
mc_errors.append(abs(mc_result))
# 準モンテカルロ法（Sobol列）
sampler = qmc.Sobol(d=2, scramble=False)
qmc_points = sampler.random(n=n)
qmc_result = np.mean([integrand(x, y) for x, y in qmc_points])
qmc_errors.append(abs(qmc_result))結果をグラフにプロットすると、以下のようになります。
このグラフから、以下のことがわかります。
- モンテカルロ法は
O(1/\sqrt{N})の収束率 - 準モンテカルロ法は
O(1/N)に近い収束率 - 特にサンプル数が多い場合、準モンテカルロ法の優位性が顕著
高次元での振る舞い
準モンテカルロ法は、低次元では優れた性能を発揮しますが、次元が高くなるにつれてその優位性が失われることがあります。
次元による性能比較
以下のグラフは、異なる次元でモンテカルロ法と準モンテカルロ法（Sobol列）の不一致度を比較したものです。
このグラフから、以下のことが観察できます。
- 低次元（1-10次元）: 準モンテカルロ法が明確に低い不一致度を示します。
- 中次元（10-15次元）: 両手法の性能差が縮小し始めます。
- 高次元（15次元以上）: モンテカルロ法の方が低い誤差を示すことがあります。
理論的解釈
高次元での性能逆転には、以下のような理由が考えられます。
- 準モンテカルロ法の限界:
- 高次元では、点列の均一性を維持することが困難になります。
- Halton列では、基底の素数が大きくなることによる数値的な不安定性が生じることがあります。
- 次元間の相関によって、構造的な偏りが生じる可能性があります。
- モンテカルロ法の安定性:
- 純粋な乱数性は、次元の影響を受けにくい性質があります。
- モンテカルロ法は、一貫して
O(1/\sqrt{N})の収束率を維持します。
実用的な対応策
高次元の問題に対処するためには、以下のようなアプローチが考えられます。
- 次元数に応じた手法の選択:
- 10次元以下: 準モンテカルロ法を積極的に活用します。
- 10-15次元: 問題の特性を考慮して、モンテカルロ法と準モンテカルロ法を使い分けます。
- 15次元以上: モンテカルロ法の使用を検討します。
- 次元削減:
- 主成分分析（PCA）などの手法を用いて、次元を削減します。
- 問題の構造を利用して、重要な次元を特定し、選択的に処理します。
- ハイブリッドアプローチ:
- 重要な低次元部分には準モンテカルロ法を、高次元部分にはモンテカルロ法を適用します。
- 両手法の利点を組み合わせた、適応的な戦略を検討します。スクランブル準モンテカルロ法と呼ばれる手法では、ランダム化を加えることで、準モンテカルロ法の欠点を補い、モンテカルロ法よりも良い結果を出せることがあります。
まとめ
本記事では、準モンテカルロ法における疑似乱数の性質について、理論と実践の両面から詳しく解説しました。
- 低不一致列: 準モンテカルロ法では、乱数の代わりに低不一致列（Halton列、Sobol列など）を用います。これらは、空間をより均一に埋めるように設計された数列です。
- 不一致度: 点列の均一性を測る指標が不一致度です。Koksma-Hlawka不等式により、数値積分の誤差は不一致度と関数の変動の積で抑えられます。
- 収束性: 低不一致列を用いることで、モンテカルロ法の
O(1/\sqrt{N})よりも速い、O(1/N)に近い収束が期待できます。 - 次元の影響: 準モンテカルロ法は低次元で特に有効ですが、高次元では性能が低下することがあります。問題に応じて、モンテカルロ法との使い分けや、次元削減などの対策が必要です。
準モンテカルロ法は、適切に利用すれば、モンテカルロ法を上回る性能を発揮します。しかし、次元の呪いの影響は完全には避けられないため、問題の特性を理解した上で適用することが重要です。
お気軽にご相談ください！
弊社のデジタル技術を活用し、貴社のDXをサポートします。
基本的な設計や実装支援はもちろん、機械学習・データ分析・3D計算などの高度な技術が求められる場面でも、最適なソリューションをご提案します。
初回相談は無料で受け付けております。困っていること、是非お聞かせください。
