準モンテカルロ法 (Quasi Monte Carlo Method; QMC)

準モンテカルロおよび、それに付随する研究・テクニックについてまとめたページです。 このページより圧倒的に詳しく厳密に書いてある日本語の解説¹ があるのでそちらを読むとより深い理解が得られると思います。

準モンテカルロ法の背景

ある関数 に対する数値積分を考えてみましょう。 素朴に区間を等間隔に切る数値積分の式は、

という形で求めることができます。 を極限に持っていくともちろん元の積分と一致しますが、この場合は が大きくなると計算量が指数関数的に増加します。

モンテカルロ積分は一様乱数を用いて積分区間から 個の点列をサンプリングして積分値を評価する手法であり、上記の問題を解決します。 式で書くと、

という形で近似計算を行います。 モンテカルロ積分の利点は、積分変数の次元 によらずに積分の評価ができる点にあります。 ただし、近似精度を上げるためには、 をより大きく取る必要があります。 具体的なオーダーの評価をしてみましょう。そのために、一様分布の点列による積分が元の積分と同じ値を取ることを使います。つまり、

で、 は一様分布です。この一様分布による積分を という形で書くとしましょう。 中心極限定理より、サンプル数 が大きい極限では、モンテカルロ積分が次の平均 と分散 を持つガウス分布に収束します。

よって誤差のオーダーは となることが分かりました。 (TODO: 別に一様分布による期待値計算をしなくても出せるので後で書き直す)

入力次元 の増加による指数関数オーダーの計算量の増加問題はモンテカルロ法で解決できたわけですが、先程の議論から誤差を小さくするために点列の取る数 をかなり大きくする必要があります。これを解決するための手法が準モンテカルロ法です。

準モンテカルロ法

準モンテカルロ法は、low-discrepancy sequence (LDS) を用いた数値積分方法です。

これが使われている応用場面は多岐にわたり、

• 物理ベースレンダリング ²
• ベイズ最適化 ³
• etc. (後で追記する)

最初に Discrepancy という不一致さの量を定義し、Discrepancy と積分誤差の評価を関係付ける Koksma-Hlawka inequality を説明し、最後に LDS を紹介します。

Discrepancy

次元ベクトルの 点集合 に対する discrepancy を次のように定義する。

ただし、 は 次元ルベーグ測度、 は超直方体

の集合で、 は超直方体 に入る の点の数を表す。

Star-discrepancy も先程の定義と殆ど同じだが、 の取る集合が、

である。

Koksma-Hlawka inequality

Van der Corput sequence ⁴

区間で一次元の LDS の一例が van der Corput sequence です。 この数列は、自然数の 進表記を逆順にしたものを小数点以下に並べたものです。

具体的に数式で表してみましょう。 正の整数 を、 進数で表現すると

となります。ここで全ての に対して を満たすとします。 このとき、 van der Corput sequence の 番目の値 は、

で表されます。

この式の定義があまり直感的でないと思う人は次のように考えても良いと思います。 進数の場合に限定しても一般性を失わないので基数が の場合で考えます。 から までを 進数で表したとします。 この 進数表記を右から左にかけて小数点第 位に設定した2進数 (これを逆 進数表記と暫定的に命名します) とそれを 進数に変換した値を下の表に示します。

このテーブルの一番右の列の値が van der Corput sequence に対応します。

van der Corput sequence の最初の 点の discrepancy は のオーダーで収束します。

数値シミュレーションをしてみましょう。

Copyimport jax; jax.config.update("jax_enable_x64", True)
import jax.numpy as jnp
import numpy as np
import matplotlib.pyplot as plt
from sympy.ntheory.factor_ import totient
 
def van_der_corput(n: int, b: int):
    x = 0
    inv_b = 1.0 / b
    while n > 0:
        x += (n % b) * inv_b
        n //= b
        inv_b /= b
 
    return x
 
def van_der_corput_seq(n: int, b: int) -> jnp.ndarray:
    return jnp.array([van_der_corput(_n, b) for _n in range(n)])
 
base = 10
N = int(3e2)
n_list = range(N)
seq = van_der_corput_seq(N, base)
 
fig, ax = plt.subplots(figsize=(6, 6), tight_layout=True)
for n in n_list:
    # Note: 縦軸の n に応じて今までに得られた数列 seq[:n] をプロットしている
    s = seq[:n]
    ax.scatter(s, [n] * len(s), s=1, c='royalblue', alpha=0.2)
ax.set_ylabel(r'$n$')
ax.set_xlabel(r'$x_n$')
ax.set_xlim([0, 1])
ax.set_ylim([0, N])
ax.invert_yaxis()
ax.set_title('Van der Corput sequence (base 10)')
plt.show()

上図より、Van der Corput sequence が稠密に を埋め尽くす様子が確認できました。

Halton Sequence ⁵

ハルトン数列は van der Corput sequence の多次元版です。 次元のハルトン数列を 番目の要素を で表すとします。 このとき任意の添字 に対する数列 は van der Corput sequence を満たします。ただし、各添字 に対して基数 は互いに素な数を使います。

例) 基数が (2, 3) の Halton sequence

CopyN = 100
x1 = van_der_corput_seq(N, 2)
x2 = van_der_corput_seq(N, 3)
cs = jnp.linspace(0, N, N)
 
fig, ax = plt.subplots(figsize=(7, 6))
sc = ax.scatter(x1, x2, c=cs, s=10)
cb = fig.colorbar(sc)
cb.set_label(r'$n$')
ax.set_xlabel(r'$x_1$')
ax.set_ylabel(r'$x_2$')
ax.set_xlim([0, 1])
ax.set_ylim([0, 1])
ax.set_title('(2, 3) Halton sequence')
plt.show()

例) (17, 19) の Halton sequence

CopyN = 100
x1 = van_der_corput_seq(N, 17)
x2 = van_der_corput_seq(N, 19)
cs = jnp.linspace(0, N, N)
 
fig, ax = plt.subplots(figsize=(7, 6))
sc = ax.scatter(x1, x2, c=cs, s=10)
cb = fig.colorbar(sc)
cb.set_label(r'$n$')
ax.set_xlabel(r'$x_1$')
ax.set_ylabel(r'$x_2$')
ax.set_xlim([0, 1])
ax.set_ylim([0, 1])
ax.set_title('(17, 19) Halton sequence')
plt.show()

Sobol Sequence ⁶

ある奇数 に対して、 と定義したとき、 の 進数を

で表すことにする。 に対して を満たす。

Original Sobol sequence は

で表される。ここで は の2進数表現 の各 bit である。 , , のとり方の一例を表にすると、

となる。ここで に対して のとり方には任意性が多少あるため、どのように選択すれば良いかという問題が残る。

を構成するために、次の原子多項式

を利用する。 はオイラー関数であり、係数 は 0 か 1 の値を取る。係数は任意に選んで良い。 この 次元原始多項式 は 個存在する。 この の係数を用いて を

で計算する。

https://i05nagai.github.io/memorandum/math/sobol_sequence.html

Antonov と Saleev ら ⁷ によると、Sobol 列の 番目の要素 を

で表しても asymptotic discrepancy には影響しない。ここで は のグレイコード表現 である。

Sobol sequence の 番目の要素を

で定義する。ここで は の2進表現の一番右のゼロビットを表す。

Kuronecker sequence

が与えられたとき、kuronecker sequence の 番目の要素を

で定義する。