準モンテカルロ法で使用される低不一致列の解説、次元数による有効性の違い

はじめに

モンテカルロ法で数値積分の誤差を半分にするには、サンプル数を4倍にしなければなりません。収束がこれほど遅いのは、乱数が偶然つくる点の偏り（点が固まる領域と空いてしまう領域）が原因です。準モンテカルロ法は、この偏りを最初から抑えるように設計された「低不一致列」という数列で乱数を置き換え、同じサンプル数でより小さい誤差を狙います。ただし、この優位性は次元が上がると失われていきます。低不一致列がなぜ効くのかを理論から追い、Halton列とSobol列を実装したうえで、何次元あたりで乱数に追いつかれるのかを実験で確かめます。

対象読者:

数値積分や数値計算の効率化に関心のある方
モンテカルロ法の基礎知識があり、さらに理解を深めたい方
準モンテカルロ法の実装に興味がある方

記事のポイント:

準モンテカルロ法で用いられる低不一致列の概念と、その利点
代表的な低不一致列（Halton列、Sobol列など）の特性と実装方法
不一致度（discrepancy）の概念と、それが数値積分の精度に与える影響
準モンテカルロ法の性能と次元の関係、高次元の問題への対処法

サンプル点の選び方：乱数vs格子点vs低不一致列

数値積分や数値計算では、計算に使う点（サンプル点）の選び方ひとつで結果の精度が変わります。代表的な3つの方法を比べます。

方法	特徴
乱数（モンテカルロ法）	- ランダムに点を選ぶ - 確率的な意味で、どの領域にも均等に点が分布 - 点と点の間に偏りが生じる可能性がある
格子点	- 等間隔に点を配置 - 配置が規則的 - 特定の方向に点が並んでしまう - 次元が増えると、必要な点の数が爆発的に増加する（次元の呪い）
低不一致列（準モンテカルロ法）	- 乱数と格子点の中間的な性質 - 規則的すぎず、ランダムすぎない - どの方向から見ても分布が均一に近い

格子点による方法は、低次元では効率的ですが、次元が増えると必要な点の数が指数関数的に増大する「次元の呪い」が問題となります。乱数を用いるモンテカルロ法は、次元の呪いには強いものの、点の偏りによって収束が遅いという欠点があります。

低不一致列は、これらの問題を解決するために考案されました。各点は、それまでに配置された点との位置関係を踏まえて置かれ、空間全体を均一に埋めるように設計されています。格子点法の規則性に起因する問題を避けつつ、乱数よりも効率よくサンプリングできます。

乱数の一様性と不一致度

「良い」サンプル点の条件を考えるために、一様性と不一致度という2つの見方を導入します。

一様性は、直感的には「点がどの領域にも同じくらいの密度で分布している」ことを指します。ただ、数値計算ではこの直感だけでは足りず、点列が空間をどれだけ均一に埋めているかを定量的に測る尺度として「不一致度（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列）の不一致度を比較したものです。

10次元程度までは準モンテカルロ法の不一致度が明確に低く、10〜15次元で両者の差が縮まり始めます。15次元を超えると、モンテカルロ法の方が低い誤差を示すことがあります。

なぜ高次元で逆転が起きるのか

逆転の主因は、低不一致列の均一性が次元の増加に耐えられなくなることです。高次元では点列の均一性を保つのが難しくなり、Halton列では基底の素数が大きくなることによる数値的な不安定性も生じます。次元間の相関から構造的な偏りが出る場合もあります。一方、純粋な乱数は次元の影響を受けにくく、何次元でも $O(1/\sqrt{N})$ の収束率を保ちます。低次元で効いていた設計上の工夫が、高次元では偏りの原因に変わるという構図です。

次元数を見て手法を選ぶ

実務で高次元の積分に当たったときは、まず次元数で当たりをつけます。10次元以下なら準モンテカルロ法を第一候補にし、10〜15次元では問題の特性を見て両手法を使い分け、15次元を超えたらモンテカルロ法を検討します。

次元そのものを減らせる場合もあります。主成分分析（PCA）で実効的な次元を落とす、あるいは問題の構造から寄与の大きい次元を特定し、その低次元部分にだけ準モンテカルロ法を適用する方法です。また、スクランブル準モンテカルロ法のようにランダム化を加える手法は、高次元での準モンテカルロ法の欠点を補い、モンテカルロ法よりも良い結果を出せることがあります。

まとめ

低不一致列は空間を均一に埋めるように設計された数列で、Koksma-Hlawka不等式が示すとおり、不一致度を下げることがそのまま数値積分の誤差を抑えることにつながります。2次元の積分実験では、Sobol列がモンテカルロ法の $O(1/\sqrt{N})$ を上回る $O(1/N)$ に近い収束を示しました。一方で不一致度の比較からは、15次元を超えるあたりでこの優位性が失われることも確認できました。準モンテカルロ法は無条件に速いわけではなく、効く範囲が次元数で決まる手法です。

導入するなら、まず自分の問題の実効的な次元数を見積もるところから始めてください。10次元以下であれば、SciPyのqmc.Sobolに乱数生成を差し替えるだけで収束の改善を試せます。それより高次元なら、次元削減やスクランブル化を組み合わせるか、モンテカルロ法のままにする判断も含めて、実際の誤差で比較するのが確実です。