　この記事はNEEとMISの実装 (MIS編) - MochiMochi3DのおまけでMISのバランスヒューリスティックの導出はどうやっているのかの流れを書いたものです。

MISの理論的な話

　MISによるレンダリング結果を見れば大幅に分散（ノイズ）が減っていることが見てわかると思います。実用的にはバランスヒューリスティックはとても良い手法であることはわかりましたが、どうやってバランスヒューリスティックなんてものが出てきたのかということが疑問に浮かびます。また、これより最適な方法はないのかということも頭に浮かびます。

　MISの製作者であるVeachという方の論文[Multiple Importance Sampling]を見るとちゃんとその辺のことが書いてありました。この記事ではこの論文にあった証明を書いていこうと思います。ただし自分の解釈で書いていくため、間違いがある可能性があることに留意してください。

multi-sapmle Model

　今、ある関数 $f$ についてある積分範囲 $\Omega$ で積分値を調べるとする。その積分は

$\displaystyle{\int_{\Omega} f(x) d\mu(x)} \tag{1}$

と表します。これに対してモンテカルロ積分をして数値積分することを考えます。

　 $f(x)$ に対してn個のサンプリング戦略があるとします。通常のモンテカルロ積分であれば、これらから１つだけ選んでサンプリングを行いますが一般的には１つのサンプリング戦略がすべての場合で効率的であるとは限りません。場合によっては他のサンプリング戦略がいいこともありますし、同等の値を出すような戦略も存在します。

　なので、１つのサンプリング方法を使い続けるのではなく、うまいこと効率のいいサンプリング戦略をｎ個の中から必ずとれるようにしたいことになります。そこで、「ｎ個のサンプリング戦略を別々に行って、その結果をうまく合わせれば弱点を補えるのではないか？」というアイディアによって生まれたのがMultiple Importance Samplingとなります。

　Veachの論文ではいくつかサンプリングした結果をウェイト $\omega$ をかけて総和を取るという形でサンプリング戦略の合成考えています。このモデルはmulti-sample Modelと呼ばれています。

　 $f(x)$ に対してサンプリング戦略が $n$ 個存在しているものとして、multi-sample Modelは

$\displaystyle{ F = \sum_{i = 1}^{n} \frac{1}{n_i} \sum_{j=1}^{n_i} \frac{\omega_i(X_{i,j}) f(X_{i,j})}{p_i(X_{i,j})}} \tag{2}$

と表されます。各文字について

$n_i$ : i番目サンプリング戦略でのサンプル数,

$X_{i,j}$ : i番目サンプリング戦略でj番目のサンプル（乱数）

$\omega_i$ : i番目サンプリング戦略におけるウェイト関数

$p_i$ : i番目サンプリング戦略におけるPDF

これが意味する所はそれぞれのサンプリング戦略で別々に $n_i$ 回サンプリングを行い、その平均を合わせるという形で求めるということである。この際、この合成がうまいこと行くように調整項としてウェイト関数 $\omega_i$ をサンプリングごとに入ってきます。

　ウェイト関数 $\omega_i$ は当然何でもいいというわけではなく、条件として

$\displaystyle{ \sum_{i=1}^{n} \omega_i = 1}$

を満たす必要があります。こうなる理由は次に説明します。

　こうしてアイディアであった複数のサンプリングを使ったモンテカルロ積分を考えることができました。これでうまいこと行くようなウェイト関数 $\omega_i$ を作れば複数のサンプリング戦略を混ぜたMultiple Importance Samplingができるというわけです。

multi-sample Modelの収束値

　こうした複数サンプリングのモデルとしてmulti-sample Modelが作られましたが、これが実際に通常のモンテカルロ積分のように積分値にちゃんと収束しなければ意味がありません。このモデルの期待値について求めると

$\displaystyle{E(F) = \int_{\Omega} \sum_{i = 1}^{n} \frac{1}{n_i} \sum_{j=1}^{n_i} \frac{\omega_i(x) f(x)}{p_i(x)} p_i(x) d\mu(x) }$

$\displaystyle{ = \sum_{i = 1}^{n} \frac{1}{n_i} \sum_{j=1}^{n_i}\int_{\Omega} \frac{\omega_i(x) f(x)}{p_i(x)} p_i(x) d\mu(x) }$

$\displaystyle{ = \sum_{i = 1}^{n} \frac{1}{n_i} n_i\int_{\Omega} \omega_i(x) f(x) d\mu(x) }$

$\displaystyle{ = \sum_{i = 1}^{n} \int_{\Omega} \omega_i(x) f(x) d\mu(x) }$

$\displaystyle{ = \int_{\Omega}\sum_{i = 1}^{n} \omega_i(x) f(x) d\mu(x) }$

となります。この値は(1)式に等しくならなくてはならないため、ウェイトの条件として

$\displaystyle{ \sum_{i=1}^{n} \omega_i = 1} \tag{3}$

得られるわけです。なのでこの条件下のウェイトであれば、multi-sample Modelはちゃんと積分値へと収束することが求められます。

$\displaystyle{E(F) = \int_{\Omega} f(x) d\mu(x) }$

バランスヒューリスティックの導出

　MISのウェイト関数としてバランスヒューリスティックというものがあります。

$\displaystyle{\omega_i(x) = \frac{n_i p_i(x)}{\sum_{k} n_k p_k(x)}} \tag{4}$

なんとなく納得できる形ではありますがこれがどこから来たかというと、分散が最も小さくなるようなウェイトの取り方を導出するとこれが得られます（厳密にはさらに分散を小さくするウェイトの取り方がある可能性はあります）。

multi-sample Modelの分散を求める

　まず、簡単のためにiサンプリング戦略におけるj番目のサンプリングを $F_{i,j}$ と表すことにします。

$\displaystyle{ F_{i,j} = \frac{\omega_i(X_{i,j}) f(X_{i,j})}{p_i(X_{i,j})} } \tag{5}$

この時、iサンプリング戦略の収束値 $\mu_i$ は

$\displaystyle{ \mu_i = E(F_{i,j}) }$

$\displaystyle{ = \int_{\Omega} \omega_i(x) f(x) d\mu(x) }$

となります。

　次にmulti-sample Modelの分散について求めます。multi-sample Modelは

$\displaystyle{F = \sum_{i = 1}^{n} \frac{1}{n_i} \sum_{j=1}^{n_i} F_{i,j}}$

となっています。これの分散は

$\displaystyle{V(F) = V(\sum_{i = 1}^{n} \frac{1}{n_i} \sum_{j=1}^{n_i} F_{i,j})}$

となります。定数である $n_i$ を外に出し（２乗になる）、分散の線形性(サンプル[tex:F_{i,j}がそれぞれ独立なので)から

$\displaystyle{= \sum_{i = 1}^{n} \frac{1}{n_i^2} \sum_{j=1}^{n_i} V(F_{i,j})}$

分散と期待値の関係性から、期待値に書き換えると

$\displaystyle{= \sum_{i = 1}^{n} \frac{1}{n_i^2} \sum_{j=1}^{n_i} E(F_{i,j}^2) - \sum_{i = 1}^{n} \frac{1}{n_i^2} \sum_{j=1}^{n_i} E(F_{i,j})^2}$

$F_{i,j}$ を展開して表すと

$\displaystyle{= \sum_{i = 1}^{n} \frac{1}{n_i^2} \sum_{j=1}^{n_i}\int_{\Omega} \frac{\omega_i(x)^2 f(x)^2}{p_i(x)^2} p_i(x)s d\mu(x)- \sum_{i = 1}^{n} \frac{1}{n_i^2} n_i \mu_i^2}$

第一項目について総和を取って、積分の中に入れると

$\displaystyle{= \int_{ \Omega} \sum_{i=1}^{n} \frac{\omega_i(x)^2 f(x)^2}{n_ip_i(x)} d\mu(x)- \sum_{i = 1}^{n} \frac{1}{n_i^2} n_i \mu_i^2} \tag{6}$

と求めることができます。これがmulti-sample Modelの分散になり、これをいかに小さくするかが問題となっていきます。この時、色々と変数はありますが結局のところ私たちが無理やり値を変えていいのはウェイト関数 $\omega_i(x)$ だけです。なので、分散を最小化するようなウェイト関数について考えていけばいいのではないかということになります。

第一項目の最小化

　まず、第一項目について最小化を考えていきます。

$\displaystyle{ \int_{ \Omega} \sum_{i=1}^{n} \frac{\omega_i(x) f(x)}{n_ip_i(x)} d\mu(x)} \tag{7}$

　この値は積分値で一見最小値をどうすればいいかがわからないと思います。この式は積分であるため、積分値は結局のところ各点 $x$ で

$\displaystyle{\sum_{i=1}^{n} \frac{\omega_i(x)^2 f(x)^2}{n_ip_i(x)}}$

という値を求めて、その総和を取っていくという話になります。仕方がって、各点 $x$ において上の値が各々最小になることができれば積分値(7)というのも最小になるはずです。　ということで、任意の点 $x$ において被積分関数が最小の値を取るようなウェイト関数を求めていく必要があります。これについては $x$ を固定して、その時のウェイト関数がどのようになるかを調べていけばいいというわけです。 $x$ を固定したものとして考えると、 $f$ は定数となり、 $x$ を除くと、

$\displaystyle{\sum_{i=1}^{n}\frac{\omega_i^2}{n_ip_i} } \tag{8}$

という式になります。これについて最小化していけば(7)式が最小化されることになります。　

　これを最小にするウェイト関数 $\omega_i$ について求めることになるのですが、これには式(3)という条件が付きます。こうした束縛条件が存在する時の最小化はラグランジュの未定乗数法を用いることで行えます。 ja.wikipedia.org

　束縛条件(3)と変数 $\lambda$ を用いて、

$\displaystyle{\sum_{i=1}^{n}\frac{\omega_i^2}{n_ip_i} + \lambda (\sum_{i=1}^{n} \omega_i - 1)} \tag{9}$

という式を考えます。ラグランジュの未定乗数法によって、これらは各偏微分が０になります。そのため、偏微分すると

$\displaystyle{\frac{\partial}{\partial \omega_i}(\sum_{i=1}^{n}\frac{\omega_i^2}{n_ip_i} + \lambda (\sum_{i=1}^{n} \omega_i - 1) )= 0}$

$\displaystyle{ 2 \frac{\omega_i}{n_i p_i} + \lambda= 0}$

$\displaystyle{ \omega_i= \frac{-1}{2}n_i p_i\lambda} \tag{10}$

と各ウェイト関数 $\omega_i$ と $\lambda$ の関係式が求められます。よって束縛条件を含めるとn+1個の方程式

$\displaystyle{ \omega_1= \frac{-1}{2}n_1 p_1\lambda , \omega_2= \frac{-1}{2}n_2 p_2\lambda,...,\omega_n= \frac{-1}{2}n_n p_n\lambda}$

$\displaystyle{ \sum_{i=1}^{n} \omega_i = 1}$

が求められます。この時、ウェイト関数の方程式を束縛条件に代入すると

$\displaystyle{ \frac{-1}{2}n_1 p_1\lambda +\frac{-1}{2}n_2 p_2\lambda+...+\frac{-1}{2}n_n p_n\lambda = 1}$

$\displaystyle{ \frac{-1}{2} \lambda \sum_{k=1}^{n}n_k p_k = 1}$

$\displaystyle{\lambda = - \frac{2}{\sum_{k=1}^{n}n_k p_k} }\tag{11}$

$\lambda$ が得られます。これを代入すると

$\displaystyle{ \omega_i= \frac{-1}{2}n_i p_i (- \frac{2}{\sum_{k=1}^{n}n_k p_k}) }$

$\displaystyle{ \omega_i= \frac{n_i p_i}{\sum_{k=1}^{n}n_k p_k} } \tag{12}$

バランスーヒューリスティックが求められます。すなわち、バランスーヒューリスティックとは式(7)を最小化するウェイト関数なのであります。

第二項目の最小化

　第二項目について

$\displaystyle{\sum_{i = 1}^{n} \frac{1}{n_i} \mu_i^2} \tag{13}$

こちらも最小化を考えてみます。こちらの場合、直接ウェイト関数が出てこないため $\mu_i$ について求めていきます（ウェイト関数はすでに求めているじゃんと思われるかもしれませんが一旦それは置いておきます）。　 $\mu_i$ については束縛条件として収束値 $\mu$ との関係式

$\displaystyle{ \sum_{i=1}^{n} \mu_i = \mu} \tag{14}$

があります。なので同様にラグランジュの未定乗数で最適化することができます。

$\displaystyle{\sum_{i = 1}^{n} \frac{1}{n_i} \mu_i^2 + \lambda(\sum_{i=1}^{n} \mu_i - \mu)}$

より、

$\displaystyle{\mu_i = \frac{n_i}{\sum_{k} n_k} \mu}$

が求められます。従って最小値は

$\displaystyle{\frac{1}{\sum_{k} n_k} \mu^2} \tag{15}$

となります。

バランスヒューリスティックは最適か？

　バランスヒューリスティックはあくまで第一項目の最適化を目的としており、第二項目に関しては最小化することは考慮されていません。そのため、バランスヒューリスティックでは第二項目が最小化されておらず、完全に最適ではない可能性が存在しています。

　実際に最小化されていないことを数式として求めることは $\mu_i$ と $\omega_i$ の関係は積分の関係にあることから非常に難しいです。第二項目の最小化は現実的には特殊な状態が必要なため、バランスヒューリスティックでは第二項目が最小化されないものとして考えても問題はないと考えられます。

　そうなると第二項目の非最小化による分散の影響を考えなくてはなりません。もしかしたら、第一項目より第二項目のほうが分散の影響が強く、こっちを最適化しなくてはならないのではないかということもあり得ます。

第二項目の範囲

　こうしたことを考えるため、第二項目の上限を調べます。 $n_i$ の内、最も小さい値 $\min{n_i}$ だけでn個だけ総和を取った値を

$\displaystyle{\min_i n_i = \sum_{i=1}^{n} \min{n_i}}$

と表すこととします。この時、

$\displaystyle{\sum_{i=1}^{n} n_i \geq \min_i n_i}$

という不等式が成り立ちます。これを(13)式に代入すると

$\displaystyle{\sum_{i=1}^{n} \frac{1}{n_i} \mu_i^2 \leq \frac{1}{\min_i n_i} \sum_{i = 1}^{n}\mu_i^2}$

この時、さらに

$\displaystyle{\sum_{i=1}^{n} \frac{1}{n_i} \mu_i^2 \leq \frac{1}{\min_i n_i} \sum_{i = 1}^{n}\mu_i^2 \leq \frac{1}{\min_i n_i} (\sum_{i = 1}^{n}\mu_i)^2}$

と求められるので、

$\displaystyle{\sum_{i=1}^{n} \frac{1}{n_i} \mu_i^2 \leq \frac{1}{\min_i n_i} \mu^2} \tag{16}$

つまり、第二項目については

$\displaystyle{\frac{1}{\sum_{k} n_k} \mu^2 \leq \sum_{i=1}^{n} \frac{1}{n_i} \mu_i^2 \leq \frac{1}{\min_i n_i} \mu^2} \tag{17}$

という範囲にあることが求められます。

理想との分散の差

　バランスヒューリスティックは完璧な最適化をすることができません。もし、仮に完璧な最適化ができるEstimater $F$ が存在するとして、パワーヒューリスティックによるEstimater $\tilde{F}$ と分散で差を取るとします。

$\displaystyle{V(\tilde{F}) - V(F)}$

この値はどれだけバランスヒューリスティックが最適化から離れているかを示すことになります。

　バランスヒューリスティックは第一項目を最適化していることは求めた通りであり、これ以上の最適化は存在しません。従って、この二つの分散での第一項目は等しい値となります。なので、この差は第二項目の差となります。

　当然、完璧な最適化をしている $F$ は第二項目も最適化をしています。すなわち第二項目は最小値(15)を取っていることになります。一方で、 $\tilde{F}$ は第二項目はどの値を取っているかはわかりません。しかしながら、第二項目の範囲は式(17)に書かれている通りであり、最大でも(16)の上限が存在しています。

　なのでこの分散の差について不等式として

$\displaystyle{V(\tilde{F}) - V(F) \leq (\frac{1}{\min_i n_i} - \frac{1}{\sum_{i}n_i} ) \mu^2} \tag{18}$

という形で得ることができます。 $n_i$ は各サンプルの回数であり、実用上ではそこまでサンプリングしなくても各分母が100とか1000といったぐらいの大きさになります。また、それぞれ極端に違うということはなかなかありません。そうなるとカッコ内の値というのは結構小さいものとなります（大体0.001とかそういうレベル）。従って実用上ではこの値というのはそこまで問題になるほどの値ではないと思われます。

　こうしたことから「バランスヒューリスティックはある程度最適解から離れているかもしれないけど、その差は精々(18)程度であり、かなり最適解に近い分散を持っている」ことが式(18)の意味となります。

one-sample Modelでの分散

　Veachの論文では確率的にサンプリング戦略を選択して１サンプルを行うone-sample Modelについても記述しています。こちらについても同様にバランスヒューリスティックを考えることができます。

　実はone-sample Modelにおいてバランスヒューリスティックはmulti-sample Modelとは異なり、最も分散が小さくなるウェイト関数であることが証明されているらしいです。（証明についてはちょっと理解できてないです）

$\displaystyle{V(\tilde{F}) \leq V(F)}$

MISはどんなシーンでも早くなるのか？

(ちょっとここら辺は自信がないです　飛ばしても大丈夫です）

　MISの画像とPathTracer単体、NEE単体のものと比べるとノイズは全然なく、結構強力な手法なんだとわかりました。なので何でもかんでもMIS使えばいいんじゃないみたいなことを考えたのですが、ほんとにMISってどんなシーンでも分散を減らしてくれるような万能なんだろうかっていう疑問がありました。今回の式で折角分散がわかったので１つのサンプリング戦略しか使わないサンプリング(PathtracerとかNEEとか)との分散の比較ができるんじゃないかってことで一応考えてみたことを書いてみます。

　ある１つのサンプリング戦略iのみを使うEstimaterを $F'$ とします。サンプル数を $N$ として、 $F'$ の分散は

$\displaystyle{V(F_{pt}) = \frac{1}{N} \int_{\Omega} \frac{f(x)^2}{p(x)} d\mu(x) - \frac{1}{N} \mu^2}\tag{19}$

という形になります。（PathTracerとかの分散でも見られるやつですね）

　multi-sample Modelについて、 $\omega_i = 1$ としてほかのウェイト関数 $\omega_k = 0$ というウェイト関数を考えた際、この時のmulti-sample Modelは $F'$ になるはずです。なので $F'$ は一つのmulti-sample Modelとしても考えられるはずです。従って式(19)も今まで考えていた第一項目、第二項目として分離できるはずです。　

バランスヒューリスティックとの比較

　まず、第二項目についてですがこれは(15)の分母部分がサンプル数の総和、つまり全サンプル数を意味しているため、(19)式の第二項目は(15)式と一致しており最適化されていることがわかります。従って第二項目を注目すればバランスヒューリスティックに勝っていることになりますが、その差というのは式(18)で表したように微小な値であるのであまり問題はありません。　

　結局のところ、大きな比較となる部分が第一項目となります。しかしながら、バランスヒューリスティックはこの第一項目を最適化したもののため、Estimater $F'$ は第一項目ではバランスヒューリスティックに同等あるいはそれより大きくなってしまうほかありません。その差については $f$ が未知の関数である以上求められませんが、サンプリング戦略が変わると分散も大きく変わるため、おそらく大きな差となるはずです。

　従って、１つのみのサンプリング戦略を使ったEstimater $F'$ は数量的に求めることはできませんがバランスヒューリスティックに打ち勝つことはほとんどないのだと考えられます。

$\displaystyle{V(\tilde{F}) \leq V(F')}$

しかしながら、すごい特殊な状況下ではバランスヒューリスティックのウェイト関数が丁度Estimater $F'$ と等しくなる場合があります。

　例えば、理想鏡面しか存在しないシーンがあった場合、PathTracerは完璧なサンプリング戦略となります。一方で、NEEなどのほかのサンプリング戦略をしたとしてもそもそもPDFが $p_i =０$ になるため、ウェイト関数がPathTracerのやつだけが１になり他の戦略ではすべて０となります。これはすなわちPathTracerのみを取ったEstimater $F'$ に相当することとなります。なので、結果としてEstimater $F'$ はバランスヒューリスティックを超えた、このシーンにおける最適解であることになります。

結論

　結論としては、「パワーヒューリスティックEstimater $\tilde{F}$ は現実的には $F'$ より分散が少ない。（ただし、コーナーケースとして $F'$ のほうが良い場合もありうる）」と考えられます。長々と話しましたが結局MISって全部のサンプリング戦略やるんだから絶対ノイズ減るじゃんっていう話です。

　ただし、これはあくまでサンプル数が同じ場合における話です。現実的にはMISはサンプリング回数が多いため、実行時間が長く、同じ実行時間で見ればほかの方法がサンプル数を回せてノイズを少なく出せるということもあり得ます。でもまあ、同じ時間でもMISのほうが明らかにノイズが小さくなり、また色々とサンプリング戦略の弱点をカバーできるので不自然な点もなくなりやすいです。事実上はMISを使うことが最も良いことになるのだと思います。