Monthly Archives: August 2014

Multiple importance samplingの罠

はじめに

レイトレ合宿のアドベントカレンダー8週目の記事として, Multiple importance sampling (MIS) で躓きやすい個所についていくつか実例を挙げて紹介したいと思います.

Multiple importance sampling (MIS) とは, Veachら [1], [2] によって開発されたモンテカルロ積分のための手法で, ひとつの分布から重点的サンプリングすると効率の悪い分布を, 複数の分布からのサンプルを組み合わせて用いることで効率よくサンプリングできます.

ですが, MISがかかわると非常に実装がバグりやすいことで有名です.
Bidirectional path tracing (BPT) [4], [2] がそのひとつで, 正しく実装するためにはMISに対する正しい理解と慎重な実装が求められます.
本記事では実例をあげながらなぜ躓きやすいのか, どういう部分に気を付けて実装すればいいのかについて解説していこうと思います.

MISの概要

モンテカルロ積分によって推定したい積分を

    \[ I=\int_{\Omega} f(x)d\mu(x) \]

とします. このとき, n個の異なる確率分布を考え, その密度関数をp_1,\cdots,p_nとします.
p_iに対する仮定として, {\rm supp}(f)\subseteq \cup_i^n {\rm supp}(p_i) (すべてのf(x)>0となるxで, 少なくともひとつの分布でp_i(x)>0となる) を考えます.
このとき,

    \[ F = \sum_{i=1}^n \frac{1}{n_i} \sum_{j=1}^{n_i} w_i(X_{i,j})\frac{f(X_{i,j})}{p_i(X_{i,j})} \]

multi-sample estimatorと呼びます.
ここで, w_i(x)x\in \cup_i^n {\rm supp}(p_i)に対して定義される重み関数です.

MISを正しく用いるには, 上のFがどのような場合にunbiasedになるのかを把握する必要があります.
重要な特徴として, 重み関数w_iが以下の条件を満たす場合にFはunbiasedとなります.

  • 条件1 : すべてのx\in \cup_i^n {\rm supp}(p_i)に対して, \sum_i^n w_i(x)=1.
  • 条件2 : p_i(x)=0ならばw_i(x)=0.

簡単のため, n=2のときで証明します. 一般の場合もほぼ同じだと思います.
\Omega_1 = {\rm supp}(p_1), \Omega_2 = {\rm supp}(p_2)とします.

    \begin{align*} \mathbb{E}[F] &= \int_{\Omega_1} w_1(x)\frac{f(x)}{p_1(x)}p_1(x)d\mu(x)+ \int_{\Omega_2} w_2(x)\frac{f(x)}{p_2(x)}p_2(x)d\mu(x) \\ &= \int_{\Omega_1} w_1(x)f(x)d\mu(x) + \int_{\Omega_2} w_2(x)f(x)d\mu(x) \\ &= \int_{\Omega_1\setminus\Omega_2} w_1(x)f(x)d\mu(x) \\ &+ \int_{\Omega_1\cap\Omega_2}\left( w_1(x) + w_2(x) \right) f(x)d\mu(x) \\ &+ \int_{\Omega_2\setminus\Omega_1} w_2(x)f(x)d\mu(x) \\ \end{align*}

ここで, x\in\Omega_1\cap\Omega_2に対して条件1より

    \[ \int_{\Omega_1\cap\Omega_2}\left( w_1(x) + w_2(x) \right) f(x)d\mu(x) = \int_{\Omega_1\cap\Omega_2} f(x)d\mu(x). \]

また, x\in\Omega_1\setminus\Omega_2に対して条件2よりw_2(x)=0.
条件1よりw_1(x)=1-w_2(x)=1.
ゆえに,

    \[ \int_{\Omega_1\setminus\Omega_2} w_1(x)f(x)d\mu(x) = \int_{\Omega_1\setminus\Omega_2}f(x)d\mu(x). \]

x\in\Omega_2\setminus\Omega_1のときも同様.
以上により,

    \begin{align*} \mathbb{E}[F] &= \int_{\Omega_1\setminus\Omega_2} f(x)d\mu(x) + \int_{\Omega_1\cap\Omega_2} f(x)d\mu(x) + \int_{\Omega_2\setminus\Omega_1} f(x)d\mu(x) \\ &= \int_{\Omega} f(x)d\mu(x) = I \end{align*}

で, Fはunbinasedです.

重み関数としては, balance heuristics

    \[ w_i(x)=\frac{n_i p_i(x)}{\sum_k n_k p_k(x)} \]

や, power heuristics

    \[ w_i(x)=\frac{(n_i  p_i(x))^\beta}{\sum_k (n_k p_k(x))^\beta} \]

がよく用いられます.

デルタ関数の扱い

本記事のメインパートです.

MISの扱いに関して注意する話題のひとつとして, 組み合わせる確率分布がデルタ関数のケースがあります. 光輸送シミュレーションにおいて, デルタ関数はさまざまな場所にあらわれますが, MISと同時に用いる場合いくつか注意が必要です.

サンプルできない分布は考慮しない

前節ではmulti-sample estimatorがunbiasedになる条件に付いて議論しました.
ここではその条件を拡張して, デルタ関数にも対応させることを考えます.

条件2′ : すべてのx\in \cup_i^n {\rm supp}(p_i)に対して, xp_iから生成される確率が0のとき, w_i(x)=0.

分布関数がゼロなのは「タテがゼロ」のケースで,
デルタ関数の場合は「ヨコがゼロ」のケースです.
いずれもサンプリングできません.

例として, 3つの確率密度関数p_1, p_2, p_2を考えて, そのうちのp_1x_0で定義されるデルタ関数だとしてp_1=\delta(x-x_0)とします. (実はこの書き方は数学的にあやしい : デルタ関数を測度として扱った際にルベーグ測度に対してRadom-Nikodym微分が存在しないので, 密度関数が存在しない. が, 便宜的に広く使われているのでここでは気にしないことにします^^).

ここでp_2(x_0)>0, p_3(x_0)>0ですが, \Omegaは連続なのでp_2, p_3ではX=x_0はサンプリングできません. このような場合, 条件2′を用いてw_2(x_0)=0, w_3(x_0)=0としてp_1, p_2を重みから外す必要があります. (結果として, w_1(x_0) = 1となります.)

デルタ関数の評価値について

実際のデルタ関数は1点で無限大となりますが, サンプリングの際に確率分布の評価値も同時に計算して, デルタ関数の場合には便宜的に有限の値を評価値としておくテクニックがよく用いられています. (私の実装ではcos thetaで割ったものを入れています. PBRTやmitsubaも同じようなことをやっています.)

このデルタ関数の評価値を有限にするテクニックですが, あくまでも便宜的なものなのでそのままの値を用いるとうまくいかないケースがあります. たとえばbalance heuristicsやpower heuristicsを用いている場合に, デルタ関数の評価値として有限の値を用いていると, 他の(非デルタ関数の)評価値の影響により正しい重みの値が計算できません. この場合, [2] p.314 にあるようにフラグを駆使して条件1,2′に従うような重みを計算します.

BPT用のシンプルな重み関数

最後に実際の例として, BPT用のシンプルな重み関数w_{Simple}の実装を示します.
実用上はあまり役に立たない重みですが, 「サンプルできない分布は考慮しない」ことがよくわかると思います.

BPTで用いる光路生成のための確率密度関数をp_{s,t}として, 頂点数nの光路\bar{x}_nに対して,

    \[ w_{Simple}(\bar{x}_n)=\frac{1}{\#\{ (s,t) : s+t=n \wedge P\{  p_{s,t} \text{ generates } \bar{x}_n \} > 0 \}}. \]

与えられたfullPathを生成できる確率分布の数を数えています.
デルタ関数を含む可能性のある場所として

  • s=0, t=nで, ポイントライトが用いられているとき (16行目)
  • s=n, t=0で, ピンホールカメラが用いられているとき (24行目)
  • その他のs,tで, specular BSDFとかdirect lightとかが用いられているとき (35行目)

をチェックしています.

 References

  1. E. Veach and L. Guibas, Optimally combining sampling techniques for Monte Carlo rendering, Procs. of SIGGRAPH 95, pp.419-428, 1995.
  2. E. Veach, Robust Monte Carlo methods for Light Transport Simulation, PhD Thesis, Stanford University, 1997.
  3. J. T. Kajiya, The Rendering Equation, Procs. of SIGGRAPH 86, pp.143-150, 1986.
  4. E. Veach and L. Guibas, Bidirectional estimators for light transport, Procs. of the Fifth Eurographics Workshop on Rendering, pp.147-162, 1994.