Resampled Importance Sampling & Weighted Reservoir Sampling

このページは レイトレ Advent Calendar 2020 の20日目の記事です。

Spatiotemporal Reservoir Resampling for Real-Time Ray Tracing with Dynamic Direct Lighting を読んでみて Resampled Importance Sampling と Weighted Reservoir Sampling が面白いと思ったので少しまとめてみます。
Resampled Importance Sampling は通常の Importance Sampling では用いられない pdf に従ったサンプリングも行うことができます。また、 Resampled Importance Sampling は通常の Importance Sampling よりも計算に必要なメモリ容量が多くなるので、 Weighted Reservoir Sampling を用いて計算に必要なメモリ容量を削減します。

Importance Sampling

Direct Lighting では、ある点 $y$ から方向 $\omega$ への反射輝度 $L$ は次のような積分で表されます、

$L \left( y, \omega \right) = \int_A \rho \left( y, \overrightarrow{xy} \leftrightarrow \omega \right) L_e \left( x \rightarrow y \right) G \left( x \leftrightarrow y \right) V \left( x \leftrightarrow y \right) dA_x \tag{1}$

$\rho$ は BSDF 、 $L_e$ は光源の放射輝度、 $G$ は幾何項、 $V$ は可視度項になります。簡単のために $y$ と $\omega$ の表記を省略すると式 (1) は、

$L = \int_A f \left( x \right) dx, \qquad where \quad f(x) \equiv \rho \left( x \right ) L_e \left( x \right) G \left( x \right) V \left( x \right) \tag{2}$

と表せます。

式 (2) の積分は分析的に解くことは難しいため、通常 Monte Carlo 法を用いて解析的に解きます。 Monte Carlo 法は確率的手法であるため分散によるノイズが発生します。 分散を低減するために Importance Sampling (以降 IS と表記) がよく用いられます。
IS では $N$ 個のサンプル $x_i$ を pdf $p \left( x_i \right)$ に従ってサンプリングします。その推定量は以下のようになります、

$\left< L \right>_{is}^{N} = \frac{1}{N} \sum_{i=1}^{N} \frac{f \left( x_i \right)}{p \left( x_i \right)} \approx L \tag{3}$

IS は $f \left( x \right) \neq 0$ となる $x$ に対して常に $p \left( x \right) \gt 0$ となる場合に unbiased となります。 また、 $p \left( x \right)$ が $f \left( x \right)$ に対して相関が強いほど分散を低減することができます。
IS を用いるためには pdf で定義される分布に従ってサンプルを生成できる必要があります。一般的にサンプリングの方法として inverse cumulative distribution function (逆CDF) や rejection sampling が用いられます。

Resampled Importance Sampling

Direct Lighting を IS を用いて計算する場合は $\rho \cdot L_e \cdot G \cdot V$ に比例した pdf に従ってサンプリングするのが理想ですが、このような pdf は closed-form で表現できず複雑なため 逆CDF によるサンプリングは難しくなります。
Resampled Importance Sampling (以降 RIS と表記) はこのような複雑な pdf を用いて IS を行うことができます。 RIS では2つの pdf を用います。1つはソース pdf $p$ で $f$ との相関は強くないが 逆CDF ができる pdf です (例えば $p \propto L_e$ など)。もう1つは ターゲット pdf $\hat{p}$ で $f$ との相関は強いが正規化が難しく 逆CDF ができない pdf です (例えば $\hat{p} \propto \rho \cdot L_e \cdot G$ など)。

先ず、 $M$ 個 ($M \ge 1$) の候補サンプル $\boldsymbol{x} = \left\{ x_1, \ldots, x_M \right\}$ を ソース pdf $p$ に従ってサンプリングします。その候補サンプルの中から $z$ 番目、 $x_z$ ($z \in \left\{1, \ldots, M \right\}$) を重み $w \left( x \right) = \frac{\hat{p} \left( x \right)}{p \left( x \right)}$ に従ってランダムに選びます。候補の中から $x_z$ が選択される確率は、

$p \left( z | \boldsymbol{x} \right) = \frac{w \left( x_z \right)}{\sum_{i=1}^{M} w \left( x_i \right)} \tag{4}$

となります。上記の手順でサンプリングされた $x_z$ はターゲット pdf $\hat{p}$ に近似した分布になります。この $x_z$ を使った推定量は以下のようになります、

$\left< L \right>_{ris}^{N,M} = \frac{1}{N} \sum_{i=1}^{N} \left( \frac{f \left( x_{iz} \right)}{\hat{p} \left( x_{iz} \right)} \cdot \left( \frac{1}{M} \sum_{j=1}^{M} w \left( x_{ij} \right) \right) \right) \tag{5}$

式 (5) の直感的な理解は、この推定量は pdf $\hat{p} \left( x_z \right)$ に従って $x_z$ をサンプリングしているように振る舞いますが、 実際には $\hat{p}$ は正規化されておらず $x_z$ は $\hat{p}$ では無く $\hat{p}$ に近似した分布に従ってサンプリングされているため、 $\frac{1}{M} \sum_{j=1}^{M} w \left( x_j \right)$ の項でその補正をかけている形になります。

式 (5) は $M=1$ の時は $p$ を使って IS しているのと同じになり、 $M=\infty$ の時は $\hat{p}$ を使って IS しているのと同じになります。 $M$ が大きくなるほど $\hat{p}$ に近い分布で IS できますが計算量は多くなります。
以下の図は[1]からの引用ですが、ソース pdf $p$ を一様分布、ターゲット pdf $\hat{p}$ を $\hat{p} \propto cos \left( \theta \right) + sin^4 \left( 6 \theta \right)$ として $M$ を変化させた時の $x_z$ の分布を表しています、

RIS は $M,N \ge 1$ かつ $f \left( x \right) \neq 0$ となる $x$ に対して常に $p \left( x \right) \gt 0, \hat{p} \left( x \right) \gt 0$ となる場合に unbiased となります。

Weighted Reservoir Sampling

RIS を行う際に候補サンプルを一度に全て保持しようとすると $M$ 個分の候補サンプルのメモリ容量が必要となります。 Weighted Reservoir Sampling (以降 WRS と表記) を RIS と組み合わせることで $M$ 個の候補の中から $n$ 個を選択する場合 $n (+1)$ 個分のメモリ消費のみで行うことができます。ここでは $n=1$ とした場合、つまり候補サンプルの中から1個を選択する場合について説明します。
WRS は事前に候補数 $M$ を知っておく必要は無く、 $\left\{ x_1, x_2, \ldots, x_M \right\}$ を全て保持しておく必要もありません。例えば候補サンプルが $x_1$ から順にストリーミング入力されて総数 $M$ がわからない場合でも 1 個を選択することができます。

WRS では候補サンプルを $x_1$ から順に処理していき、 resevoir と呼ぶ入れ物 (ここでは候補サンプル 1 個分のメモリ容量) に確率的に挿入 (既に他の候補が入っている場合は交換)していきます。 $x_1$ から処理を始め $m$ 番目の $x_m$ までの処理を行っている時、 $x_i$ が reservoir に入っている確率は $\frac{w \left( x_i \right)}{\sum_{j=1}^{m} w \left( x_j \right)}$ です。ここで用いる $w \left( x \right)$ は RIS で用いる重みと同じものです。次の $x_{m+1}$ の処理では、reservoir に入っているものと以下の確率で交換します、

$\frac{w \left( x_{m+1} \right)}{\sum_{j=1}^{m+1} w \left( x_j \right)} \tag{6}$

式 (6) から逆に $x_{m+1}$ が reservoir に入らず $x_i$ が 残っている確率は、

$\frac{w \left( x_i \right)}{\sum_{j=1}^{m} w \left( x_j \right)} \left( 1 - \frac{w \left( x_{m+1} \right)}{\sum_{j=1}^{m+1} w \left( x_j \right)} \right) = \frac{w \left( x_i \right)}{\sum_{j=1}^{m+1} w \left( x_j \right)} \tag{7}$

となります。$m+1=M$ まで処理が終了した段階で $x_i$ が reservoir に入っている確率は、

$\frac{w \left( x_i \right)}{\sum_{j=1}^{M} w \left( x_j \right)} \tag{8}$

となり式 (4) と同じ形となります。このことから RIS と WRS を組み合わせても $x_z$ の分布を変えずにサンプリングできることがわかります。また、必要なメモリ容量を reservoir 分と重みの和 $\sum_{j=1}^{m} w \left( x_j \right)$ を保持する分のみに削減できます。

終わりに

RIS と WRS について調べましたが実はまだ試していません。以下の事についてまだ疑問が解決していないのでもう少し調べてみようと思っています。

1. RIS + WRS では乱数が多く必要になるが必要な個数を削減する方法はあるか
2. RIS と他のサンプリング戦略の Multiple Importance Sampling (MIS) を行う際の RIS 戦略の pdf を決定論的に計算する方法はあるか

1 については、通常の IS が 逆CDF を用いてサンプリングする際に乱数を 1 個使うのに対し、 RIS + WRS では候補サンプルのサンプリングや reservoir への確率的交換で多くの乱数を必要とします。何か必要な乱数の数を減らす方法はないでしょうか。
2 については、式 (5) から RIS で用いる pdf を $p_{ris} \left( x_z \right) = \frac{\hat{p} \left( x_{z} \right)}{\frac{1}{M} \sum_{j=1}^{M} w \left( x_{j} \right)}$ とすると、例えば MIS 時に $x_z$ から $p_{ris} \left( x_z \right)$ を評価する際に $\sum_{j=1}^{M} w \left( x_{j} \right)$ を計算するために乱数を用いて他の候補をサンプルをサンプリングする事になると思います。もっと決定論的に $p_{ris}$ を評価する方法はないでしょうか。

参考文献

1. Talbot J, Cline D, Egbert P: Importance Resampling for Global Illumination (2005)
2. Efraimidis PS, Spirakis PG: Weighted Random Sampling with a Reservoir (2006)
3. Efraimidis PS: Weighted Random Sampling over Data Streams (2015)
4. Bitterli B, Wyman C, Pharr M, Shirley P, Lefohn A, Jarosz W: Spatiotemporal Reservoir Resampling for Real-Time Ray Tracing with Dynamic Direct Lighting (2020)