Stratification by Rank-1 Lattices

by syoyo

Alexander Keller,
Stratification by Rank-1 Lattices
Monte Carlo and Quasi-Monte Carlo Methods 2002.

I implemented “stratification by rank-1 lattices” sampling method,
which generates sample points in 2D domain.

http://lucille.atso-net.jp/svn/angelina/sampling/rank1lattice/

Alexander Keller 博士の “Stratification by Rank-1 Lattices”
を実装してみました。

ランク-1 格子(Rank-1 Lattice)とは、

lattice_eq

のような形式により生成される格子点のことです[1]。
ここで、N は生成するサンプルの数、z は s 次元のベクトルであり、
その要素はどれも N と互いに素(N との最大公約数が 1)であるような
ベクトルになります。{}は値の小数部のみを取り出すことを意味します
(つまり値が [0, 1)^s の単位超立方体に収まる)。

特に、この格子点が低食い違い量列の性質を満たすときは、
良い格子(good lattice)と呼ばれています。
ランク-1 格子の代表的なものに、フィボナッチ格子や
コロボフ型(Korobov form)の格子があります。


フィボナッチ格子

フィボナッチ格子点は、以下のような式で定義されます。

fibo_lattice
こで、Fk は k 番目のフィボナッチ数になります。
ここで、Fk は k 番目のフィボナッチ数になります。
k=9(Fk=34)のときのフィボナッチ格子をプロットすると以下のようになります。

rank1_fib1.png


コロボフ型の格子点

コロボフ型の格子点は、以下のように定義されます。

z = (1, a, a*a, a*a*a, … , a^s)

このベクトル z を使って s 次元空間に n 個の格子点を生成します。
省略して (n, a) と書かれたりします。もちろん、任意の n と a が
ランク-1 格子を満たすわけではありません。
ランク-1 格子を満たすような n と a の取り方にはいろいろあるようですが、
有名なものに (2^m, 17797), (2^m, 203) があるようです。 

Cranley-Patterson rotation

ここで、 Alexander Keller 博士お得意の、Cranley-Patterson rotation
を格子点に追加することを考えます。格子点は deterministic な点列であるため、
そのままでは相関(correlation)によりレンダリング結果にエイリアシングなど
が出たりしてしまいます。Cranley-Patterson rotation は、
サンプル点に一様乱数を加えることで、この相関を取り除きます。
式としては、ランク-1 格子に単純に一様乱数 ζ ( [0,1) ) によりオフセットを追加するだけです。

{j/N z + ζ}

Cranley-Patterson rotation を各フィボナッチ格子点に 15 点づつ適用
(上記と同じ k=9 のとき。全部で 34 * 15 点)したもののプロット図は以下のようになります。

rank1_rand

Reduced Cranley-Patterson rotation

さて、Cranely-Patterson rotation(CP 回転)を追加してみましたが、
あまりランダム性がありつつも一様に点群が分布していませんので、
少し分布の性質が悪いですね。
これはつまり、乱数を [0,1) の区間で生成しているので、うまく層別化が行われていない
というのがあります。ここで、フィボナッチ格子を良く見ると、格子点という名前の通り、
サンプル点を結ぶと格子模様が現れることがわかります。

rank1_fib_basis

ここで、格子の線が囲む各ひし形の領域内で一様乱数を生成すると、
一様乱数のオフセットが sweep する領域が狭くなり、
よりよく層別化が行われるのがわかると思います。

格子が囲むひし形の領域を定義する基底ベクトル(青色、ただし必ずしも直交にはならない)は、
以下のような 2×2 の基底行列で求まります。

bmat_eq

ここで、
j1_eq
j2_eq
となります。

基底行列を B とすると、

{j/N z + B ζ}

として CP rotation を行うのを、論文では reduced Cranley-Patterson rotation と
呼んでいます。

相関サンプリング

ここで、いままでは各フィボナッチ点に対して、毎度ランダムな一様乱数を CP rotation として与えてサンプル点を生成してきましたが、ループを逆にしても数学的には積分の結果には影響がありません。つまり CP rotation によるオフセットは一度だけ計算し、それをすべてのフィボナッチ格子点に適用するということです。
これは相関サンプリングと呼ばれています。

こうすることで、かなり相関性(フィボナッチ格子のかたち)を保ったまま一様なサンプル点を生成することができます(ただ、CG の場合はこれはエイリアシングを引き起こしやくなるので、このような相関サンプリングを使う場合は注意が必要です)。

rank1_reduced_random
相関サンプリングによる reduced RP rotation。(34 * 15 samples)

rank1_reduced_stratified
相関サンプリングによる reduced RP rotation。
ランダム値は 3×5 の層別サンプリングで生成(34 * (3*5) samples)。
一様度が上がっていますが、論文ほどきれいに層別化はされていないですね

ここで、さらに論文では、RP rotation のオフセットに、(8, 3) のコロボフ格子の
次元 3 と 4 を用いた場合も示しています。

rank1_correlated
34 fibonacchi lattice points * 8 samples( (8, 3) Korobov lattices)

なんでここでいきなりコロボフ格子(のしかも次元 3 と 4)をオフセットに使う
のかどうかはよくわかりませんが、
見た目的によい結果が得られたからでしょうか。

まとめ

次は実際にこの手法をレンダラに組み込んでみて、テストしてみたいと思います。
ところで、この手法(ランク-1 格子)もまだサンプル数 N はあらかじめ既知で
すから、適応サンプリングに使うには向かないかもしれないですね。

次は実際にこの手法をレンダラに組み込んでみて、テストしてみたいと思います。
ところで、この手法(ランク-1 格子)もまだサンプル数 N はあらかじめ既知ですから、
適応サンプリングに使うには向かないかもしれないですね。

適応サンプリングに適した、
N -1 の時のサンプル位置を変えずに、N 番目の点を
プログレッシブかつ常に statified に高速に生成する
サンプリング法を求める旅は続く…

参考文献

[1] I.H. Sloan and S. Joe. Lattice methods for multiple integration. Oxford Science Publications. Oxford: Clarendon Press., 1994.
[2] http://lucille.sourceforge.net/blog/archives/000199.html 

Advertisements