最適化されたロシアンルーレット

ロシアンルーレットのルーレット確率$q_i$を根拠を持って決定したい。


ロシアンルーレット

ロシアンルーレットとはモンテカルロ積分を行う際に、寄与の評価を間引くことで実行効率を上げるテクニックだ。推定量$F$が次のような寄与$C_i$からなるとする。

$ F = C_1 + C_2 + C_3 + .... + C_n $

寄与$C_i$の評価にそこそこコストがかかる場合、全ての$C_i$の評価を行わず、$C_i$をいくつか間引くことで全体の実行効率を上げる事が出来る。単に間引くだけでは推定量$F$の結果にバイアスが生じるので、ルーレット確率$q_i$とルーレット確率で割られる前の寄与$t_i$を用いて次のようにする。

$ C_i = \begin{cases} t_i/q_i & \text{probability }q_i \nonumber \\\ 0 & \text{otherwise} \end{cases} \\\ $

これは次に示すように不偏性が保たれている。

$ \begin{align} E[C_i] & = q_i (ti/qi) + (1-q_i) 0 \nonumber \\\ & = ti \nonumber \\\ \end{align} $

簡単な数列の和の推定の例を示す。

|2|6|4|6|5|8|1|8|6|3|9|2| = 60|

この数列の総和は60となる。次にロシアンルーレットを行いこの数列の総和を推定する。ルーレット確率を75%にし、残った結果を1.0/0.75=4/3倍にする。

|-|6|4|6|-|8|1|8|-|3|9|2|x4/3 = 62.3| |2|-|4|6|5|-|1|8|6|-|9|2|x4/3 = 56.0| |2|6|-|6|5|8|-|8|6|3|-|2|x4/3 = 61.3| |2|6|4|-|5|8|1|-|6|3|9|-|x4/3 = 58.7|

いくつかの場合をピックアップして並べたが元の総和60にほぼ近いものになることがわかる。また不偏性は保たれているので$E[F]$は60になる。数列のそれぞれの値の評価にかかるコストが同じだとすると、25%間引かれたので実行にかかるコストも25%割り引かれる。一方でこの操作は分散を大きくしてしまっていることもわかる。仮に極端な例として11/12を間引いた場合を示す。

|2|-|-|-|-|-|-|-|-|-|-|-|x12 = 24| |-|6|-|-|-|-|-|-|-|-|-|-|x12 = 72| |-|-|4|-|-|-|-|-|-|-|-|-|x12 = 48| |-|-|-|6|-|-|-|-|-|-|-|-|x12 = 72| |-|-|-|-|5|-|-|-|-|-|-|-|x12 = 60| |-|-|-|-|-|8|-|-|-|-|-|-|x12 = 96| |-|-|-|-|-|-|1|-|-|-|-|-|x12 = 12| |-|-|-|-|-|-|-|8|-|-|-|-|x12 = 96| |-|-|-|-|-|-|-|-|6|-|-|-|x12 = 72| |-|-|-|-|-|-|-|-|-|3|-|-|x12 = 36| |-|-|-|-|-|-|-|-|-|-|9|-|x12 = 108| |-|-|-|-|-|-|-|-|-|-|-|2|x12 = 24|

依然として不偏性は保たれているが、分散がとても大きくなってしまったことが分かる。これらの例ではルーレット確率を全体で固定しているが、寄与毎にルーレット確率を変更したとしても全体の推定量は不偏のままになることに注意が必要だ。つまりルーレット確率を何らかの根拠に基づき、寄与毎に変動させることで生じる分散を減らすことができるのだ。

実行効率(Efficiency)

数値積分において推定量$F$の分散$\sigma^2$は寄与$C_i$の評価の数に反比例して下がる。つまり分散$\sigma^2$を半分にすることと、推定量$F$を得るための計算コスト$T$を半分にすることは実行効率の観点から見ると同じになる。ロシアンルーレットをしたことによる分散の増加を間引きによる効率化が上回れば結果的に同じ時間で分散が少ない結果が得られるということになる。それを実行効率(Efficiency)という指標にして一つに式にまとめると次のようになる。

$ \epsilon = \cfrac{1}{\sigma^2 T} $

実行効率を最大化する事を考えた場合、ロシアンルーレットの確率$q_i$はこの$\epsilon$が最大になる(と予想される)ように決定しなくてはいけない。素朴なロシアンルーレットでは前記の例のように$q_i$を固定値にしていたり、何らかの方法で寄与を推定したのちにそれをスケールした値を用いていることが多いが実行効率を最大化しているとは言い難い。ロシアンルーレットが実行効率にどのような影響を及ぼすかを考えるために、ロシアンルーレットを行う前の実行効率をまず考える。推定量$F$の分散の平均値$\sigma_0^2$と、$F$を得るために使われた計算コストの平均値を$T_0$として次のように表される。

$ \epsilon = \cfrac{1}{\sigma_0^2 T_0} $

ここで今後の計算を簡単にするため$T_0$は、現在注目している寄与の評価コストで正規化しているとする(レイトレの文脈ではシーンとの交差判定の回数を単位として実行コストを評価することが多い)。ここからロシアンルーレットを行ったことによる分散の増加と計算コストの減少を加味すれば、ロシアンルーレットをした場合の実行効率がわかる。注目している寄与は$(1-q_i)$の確率で計算されず、$T_0$は現在注目している寄与の評価コストで正規化しているので、ロシアンルーレットを行ったことによる計算コストは単純に$T_0$から$(1-q_i)$を減算すればよい。また分散の増加は次のようになる。

$ \begin{align} V[C_i] & = E[C_i^2] - E[C_i]^2 \nonumber \\\ & = (q_i (t_i/q_i)^2 + (1-q_i)0^2) - t_i^2 \nonumber \\\ & = t_i^2((1/q_i) - 1) \nonumber \end{align} $

ロシアンルーレットを行うことによる分散と評価コストの変化量が分かったので、ロシアンルーレットを行ったことによる実行効率は次のようになる。

$ \epsilon_{RR} = \cfrac{1}{(\sigma_0^2 + t_i^2((1/q_i) - 1) )(T_0 - (1 - q_i) )} \nonumber \\\ $

$\epsilon_{RR}$を最大化する$q_i$を知るために、$qi$で微分する(計算を簡単にするため、$1/\epsilon{RR}$を微分している)。

$ \begin{align} \cfrac{d(1/\epsilon_{RR})}{dq_i} & = \cfrac{d(\sigma_0^2 + t_i^2((1/q_i) - 1) )(T_0 - (1 - q_i) )}{dq_i} \nonumber \\\ & = (-t_i^2 q^{-2})(T_0 - (1-q_i)) + (\sigma_0^2 + t_i^2(1/q_i -1 ) ) \nonumber \\\ & = \sigma_0^2 + t_i^2(q_i^{-2}(1-T_0) - 1 ) \nonumber \\\ \end{align} $

$d(1/\epsilon_{RR})/dqi$が0になるところが$\epsilon{RR}$が最大値になるところなので、そのときの$q_i$を出す。

$ \sigma_0^2 + t_i^2(q_i^{-2}(1-T_0) - 1 ) = 0\\\ q_i^{-2}(1-T_0) = -\cfrac{\sigma_0^2}{t_i^2} + 1 \\\ q_i^{-2} = \cfrac{-\sigma_0^2/t_i^2 + 1}{1-T_0} \\\ q_i^{-2} = \cfrac{1}{t_i^2} \cfrac{t_i^2 - \sigma_0^2}{1-T_0} \\\ q_i = t_i \sqrt{\cfrac{T_0-1}{\sigma_0^2-t_i^2}} $

この最適なルーレット確率$q_i$を$\delta$として書きなおすと次のようになる。

$ C_i = \begin{cases} t_i/\delta & \text{probability }\delta \nonumber \\\ 0 & \text{otherwise} \end{cases} \\\ \delta = \sqrt{\cfrac{\sigma_0^2-t_i^2}{T_0-1}} $

となり、最も実行効率が良くなるルーレット確率$\delta$を得られた。

変形

最適なルーレット確率が求まったのでこれを使ってロシアンルーレットを行えばいいということになるが、考えなくてはいけないことがある。まず算出された$q_i$は1を超えうる。ルーレット確率$q_i$は1を下回らなければ使えないので単純に1でクランプする。

$ C_i = \begin{cases} min(1,t_i / \delta) & \text{probability }\delta \nonumber \\\ 0 & \text{otherwise} \end{cases} \\\ \delta = \sqrt{\cfrac{\sigma_0^2-t_i^2}{T_0-1}} $

さらにこれよりも大きな問題がある。ロシアンルーレットは確率$\delta$に基づいて寄与を評価するか否かを決めるが、その確率にすでに評価後の$t_i$が含まれているのだ。評価を間引きしたいのに、それを決定する確率に寄与の評価が必要なのでは本末転倒になってしまう。これらの問題を解決するために$\delta$を$q_i=1$のものに固定してしまうという案がある。

$ \text{where } q_i=1 \\\ t_i/\delta = 1 \\\ \delta = t_i \\\ \cfrac{\sigma_0^2-t_i^2}{T_0-1} = t_i^2 \\\ \sigma_0^2-t_i^2 = t_i^2(T_0-1) \\\ t_i^2 = \sigma_0^2 / T_0 $

この$\delta$を$\delta^{\ast}$とすると

$ \begin{align} \delta^{\ast} & = \sqrt{\cfrac{\sigma_0^2-(\sigma_0^2 / T_0)}{T_0-1}} \nonumber \\\ & = \sqrt{\sigma_0^2/T_0} \nonumber \\\ \end{align} $

となる。この$\delta^{\ast}$は固定値で、最適な$\delta$よりもすこし大きくなる(それでも$\sqrt{T_0/(T_0-1)}$倍を超える事はない)が先の問題を全て回避することができる。最終的なルーレット確率は次のようになる。

$ q_i = min(1,t_i / \delta^{\ast}) \\\ \delta^{\ast} = \sqrt{\sigma_0^2/T_0} \nonumber \\\ $

デモ

以下はR言語で記述されたデモだ。要素を評価するにはコストがかかる配列の総和を次の三つの戦略を用いて求め、標準誤差が目標値を下回った時点での総評価コストをそれぞれ求める。

  • 評価コストを無視し全て評価する(戦略1)
  • ルーレット確率固定のロシアンルーレットを行う(戦略2)
  • 平均評価コストと結果の分散から最適なルーレット確率を求めたロシアンルーレットを行う(戦略3)

※デモ、かなり怪しいので近日中に修正します。 {% highlight ruby %}

—————————

目標標準誤差

stdErrLimit <- 0.0001

—————————

徐々に期待値が大きくなる乱数

x<-0.0 pdist <- function(){ x<<-x+0.05; x<<-x%%1.0; return (runif(1,0,x)); }

—————————

戦略1: 頭から全て評価していく

totalEval<- 0 est1 <- c(1) while(TRUE) { est1[length(est1)+1] <- pdist() totalEval= totalEval+ 1 #print(var(est1)/totalEval) if( var(est1)/totalEval< stdErrLimit ) break; } print(“Strategy1”) print(totalEval)

—————————

戦略2: ルーレット確率が決め打ちのロシアンルーレット

totalEval <- 0 est2 <- c(1) rouletteProb <- 0.5 while(TRUE) { if( runif(1) < rouletteProb ) { est2[length(est2)+1] <- pdist() / rouletteProb totalEval= totalEval+ 1 # print(var(est2)/totalEval) if( var(est2)/totalEval < stdErrLimit ) break; } } print(“Strategy2”) print(totalEval)

—————————

戦略3: ルーレット確率を、評価コスト平均と推定量の分散から算出

numRR <- 1 totalEval <- 1 est3 <- c(1,0)

while(TRUE) { # 今回は毎回一回しか評価せずそれをロシアンルーレットにかけているので、 # 平均評価コストは過去のロシアンルーレット確率になる t0 <- numRR/totalEval # 推定量の分散 sigmaAst <- var(est3) # rouletteProb <- 1.0/((sigmaAst/t0)^0.5) # print(rouletteProb) if( runif(1) < rouletteProb ) { est3[length(est3)+1] <- pdist() / rouletteProb totalEval= totalEval+ 1 # print(var(est3)/totalEval) if( var(est3)/totalEval < stdErrLimit ) break; } numRR = numRR + 1 }

print(“Strategy3”) print(totalEval)

[1] “Strategy1” [1] 474 [1] “Strategy2” [1] 1759 [1] “Strategy3” [1] 114

{% endhighlight %}

全てを評価するよりも数倍早く収束した。また配列のデモよりもずっと面白いBPT(PT?)での比較はこちら

参照

[1] E. Veach, “Robust Monte Carlo Methods For Light Transport Simulation,” p.318-320, 1997.


メモ

  • 実際のレイトレでの実装は、平均評価コストと推定の分散は対象のピクセルによって大幅に変わるので過去何回のかの推定を保存して決める。今回のデモのように全ての結果を保存してルーレット確率を求めるようなことはしないようです。今回のデモのようにやると算出されたルーレット確率は一定値に収束するので平均評価コストと推定の分散がピクセルによってまちまちになる場合に適しません。

  • BPT比較画像にもありますが、固定値ルーレット確率はとんでもないピクセルをたまに生み出してしまうのでレイが変な所にハマる対処以外ではあまり使えないかもしれません。変な所にハマる対処は最大深度を設定するだけでもいいでしょうし。ちなみにレイトレ合宿1&2で提出されたUnbiasedなアルゴリズムを使っていると考えられるレイトレーサーのロシアンルーレットの実装の内訳は次のようになっていました。

  • ルーレット確率を何らかの簡易な推定から算出したロシアンルーレット: 1620(80%)

  • ロシアンルーレットなし: 320(15%)

  • ルーレット確率固定のロシアンルーレット: 120(5%)

  • 実行効率を最適化させたロシアンルーレット: 0/20(0%)

大きくまとめた1の実装も、マテリアルから固定値を引き出す実装、パススループットから出す実装、パスの深さも考慮する実装まで様々でしたが、このポストの方法を使っていると思われる実装は見つかりませんでした。(だれかがBPTを実装していたからそれには入っている?)