2003年08月14日
準モンテカルロ法
現在の lucille では、モンテカルロ法を用いていますが、そろそろ準モンテカルロ(Quasi Monte Carlo, QMC)法に移行しようと考えています。
lucille でモンテカルロ法を用いる主要な部分は以下になります。
1. アンチエイリアスのための、スーパーサンプリングのサブピクセル位置
2. ディフューズ面でのシェーディング時に、半球上にレイを飛ばすときのレイの方向
3. 面光源のサンプル
4. フォトントレーシング時のロシアンルーレット
乱数には、Mersenne Twister を使用しています。
Mersenne Twister: A random number generator
準モンテカルロ法では、乱数ではなく、低食い違い量列(Low Discrepancy Sequence, LDS)と呼ばれる、決定論的な数列を用います。つまりはすでに決まっている数列を乱数として用いるのです。
簡単に云えば、乱数を使用するのがモンテカルロ法で、乱数のかわりに低食い違い量列を用いるのが準モンテカルロ法、になります。
準モンテカルロ法は、低次元の積分においてモンテカルロ法よりも収束が早いという特徴があります。サンプル数をNとすると収束(平均自乗誤差(Root Mean Squared, RMS))は、モンテカルロ法では 1/Sqrt[N] であり、準モンテカルロ法では、(Log2[N])^d /N になるそうです(d は次元数)。
2次元の場合で考えると、モンテカルロ法を用いて精度を2倍にするときは、サンプル数を4倍にしなくてはならない(1/2のRMS = 1/2 * 1/Sqrt[N] = 1/Sqrt[4N])のに対して、準モンテカルロ法を用いて精度を2倍にしたい場合は、サンプル数も2倍(1/2のRMS = 1/2 * (Log2[N])^2/N = 1/2N)と線形で行えるということになります。
低次元の積分において、という制限ですが、CGで高次元の積分を扱うことはあまりありません。
上で述べた lucille でモンテカルロ法を使用している場合でも、
1. XY平面なので2次元
2. 半球面でのサンプリングは、一意に平面でのサンプリングに対応しているので2次元
3. 面光源は面であるので2次元
4. [0,1]の範囲なので1次元
となり、1次元と2次元の問題だけです。
これ以上の次元は、あるとしても任意の3次元ベクトルのサンプルということで3次元までぐらいではないでしょうか。
まとめると以下のようになります。
| モンテカルロ法 | 準モンテカルロ法 |
| 疑似乱数 | 低食い違い量列 |
| RMS 1/Sqrt[N] | RMS (Log2[N]^d)/N |
低食い違い量列を生成するコードは単純なので、モンテカルロ法から準モンテカルロ法に変更するのは簡単です。
少ないサンプル数で高品質な結果が得られるので、使わないわけにはいきません。
準モンテカルロ法を用いたレンダリングでは、Alexander Keller博士 が有名です。
http://www.uni-kl.de/AG-Heinrich/Alex.html
準モンテカルロ法を用いたレンダラ McRender を初めとして、準モンテカルロ法のコンピュータグラフィックスの応用で著名な業績を修めている人です。
また、同氏のサイトでは、各種準モンテカルロ乱数(低食い違い量列)を生成するコードライブラリ SamplePack が公開されていますので、低食い違い量列の実装にはこれをつかえばよいでしょう。
2003年08月16日
アンチエイリアシング
lucille に準モンテカルロ法を用いる第一歩として、まずは、スーパーサンプリングによるアンチエイリアシングのための、サブピクセルサンプル位置の生成に準モンテカルロ法(quasi-Monte Carlo, QMC)を適用しました。
今回準モンテカルロ法に用いた低食い違い量列(Low Discrepancy Sequence, LDS)は Hammersley points です。この数列は QMC で多くの場合に用いられている Halton 列(Halton Sequence)よりもより一様性があるようです。
Hammersley points およびそれを用いたアンチエイリアシングには、
Strictly Deterministic Sampling Methods in Computer Graphics
Alexander Keller, (mental images technical report, 2001) in "Monte Carlo Ray Tracing", SIGGRAPH'2003 Course #44, San Diego, July 2003.
を参照しました。てゆうかコードが全部載っているので、そのまま写しただけです。
通常のモンテカルロ法(層別サンプリング)と、準モンテカルロ法の結果を以下に示します。サンプル数は 1 ピクセルあたり 2x2 の 4 サンプルです。
| MC | QMC |
 |  |
 |  |
同じサンプル数ですが、明らかに QMC の方が画質が上です。
すごいですね、QMC。これで少しは Arnold に近づいたかな。
2003年09月17日
QMC: van der Corput 列
準モンテカルロ法(quasi-Monte Carlo, QMC)では、乱数の代わりに低食い違い量列(Low Discrepancy Sequences, LDS) と呼ばれる数列を用います。
この低食い違い量列は、"一様"に散らばり、かつ"規則的なランダム"性があることが特徴です。
低食い違い量列による準モンテカルロ法は、乱数によるモンテカルロ法よりも真の値への収束が早くなるので、CGや金融などのスピードが求められる分野で特に用いられています。
乱数にもいくつかの生成方法があるように、LDS にもいくつかの種類があります。
LDS の数学的な理論、食い違い度(discrepancy)の測り方はとりあえずすっとばします。
なぜかというと、LDS について調べようとして LDS に関する文献を当たると、そのほとんど全てが、
"...。詳しい数学的な証明は [Niderreiter92] で"
なんて感じで、Harald Niderreiter 博士の以下の文献を参照しているからです。
Harald Niederreiter
"Random Number Generation and Quasi-Monte Carlo Methods"
Cbms-Nsf Regional Conference Series in Applied Mathematics, No 63, 1992.
というわけで、LDS の数学的な詳細についての解説はこの本が入手できてからにします。
とりあえず今回は、実際に LDS の数列とはどんなものかについて述べます。
van der Corput 列(van der Corput sequence)
まずは van der Corput 列です。
van der Corput 列は、最も単純な LDS になります。
van der Corput 列をはじめ、すべての低食い違い量列は、根基逆関数(radical inverse function) を元にして生成されます。
根基逆関数
10進法の整数 n を、 b 進法(b >= 2) で表し、このときの j+1 桁目の数字を a_{j}(n) とします(LDSの文献では以下の式の定義が j = 0 から始まるのですが、日本語で 0 桁目というのは変なので j+1 桁目という表現を使いました)。
すると、10進法での n は、
で表すことができます。たとえば 10 進数での 6 は 2 進数で表すと 110 になるので、a(0) = 0, a(1) = 1, a(2) = 1 になります。
このとき、基数 b の根基逆関数は、以下で定義されます。
つまり、b 進法での \Phi_{b}(n) は、 n を小数点で対称に折り返したものになります。また、 n >= 0 の全ての整数 n で、 \Phi_{b}(n) は区間 [0, 1) に入ると云う有益な性質があります。
以下は基数 2 のときの van der Corput 列です。
van der Corput列は、1 次元の LDS になります。
基数 b の取り方ですが、なるべく最小の素数(2, 3, 5, 7, ...) にするのがよいそうです。 これは後の多次元での LDS でのほうが重要になります。
プログラム
van der Corput 列は、a_{b}(n) を求めればよく、これは剰余で簡単に処理できます。
以下は RenderPark のソースコードから。
double vdC(int i, int base)
{
double h=0., f, factor;
f=factor=1./(double)base;
while (i>0) {
h += (double)(i%base) * factor;
i /= base;
factor *= f;
}
return h;
}
これは任意の基数に対応していますが、基数が 2 のときなどはシフト演算などで高速に処理できます。また、van der Corput 列に限らず、 LDS 生成のより高速なバージョンは、Ilja Friedel 氏および Alexander Keller 氏による、
Caltech Multi-Res Modeling Group - low discrepancy sequence library
や、Thomas Kollig 氏および Alexander Keller 氏による、
が参考になると思います。
LDS でのサンプル点の生成
van der Corput 列で、サンプル点を取りたい時は、i=0 から始めて vdC() を呼び出します。N個のサンプル点が取りたいときは、単純に N-1 で打ち切るだけでよく、どんな N にたいしても、 van der Corput 列は常に、[0, 1) で一様かつランダムに散らばるという有益な性質があります。
グラフ
層別サンプリングによる乱数。N=10
van der Corput 列。N=10, 基数=2。
2003年09月18日
QMC: Halton 列
Halton 列(Halton sequence) は、 van der Corput 列を次元数分並べたベクトルになります。
たとえば、2次元の Halton 列は、(\Phi_{b1}(n), \Phi_{b2}(n)) と、2 つの van der Corput 列のベクトルになります。
2 次元以上では、サンプルがベクトル(点)になるので、Halton 点群(Halton point set) とも呼ばれるようです。
s 次元の Halton 列 は以下で定義されます。
x_{n}^{s} は、n 番目のサンプル点です。ここで、基数 b1, b2, ... には、重ならない素数で、なるべく最小のものを割り当てます。実際には、 2, 3, 5, 7, ... と順に素数を割り当てていけばよいでしょう。
例として、b1 = 2, b2 = 3 のときの 2 次元の Halton 列を示します。
2 次元の Halton 列 (b1=2, b2=3)。
1 次元の Halton 列というのは、 van der Corput 列と同等です。
プログラム
Halton 列を求めるプログラムです。次元数は 10 以下に限定しています。
#include <stdio.h>
double vdC(int i, int base)
{
double h=0., f, factor; f=factor=1./(double)base;
while (i>0) {
h += (double)(i%base) * factor;
i /= base;
factor *= f;
}
return h;
}
void halton(double *x, int i, int dim)
{
int primes[10] = {2, 3, 5, 7, 9, 11, 13, 17, 19, 23};
int j;
if (dim >= 10) {
fprintf(stderr, "too big dimension\n");
return;
}
for (j = 0; j < dim; j++) {
x[j] = vdC(i, primes[j]);
}
}
グラフ
乱数(層別サンプリング, N=16)
Halton 列(2 次元, N=16, b1=2, b2=3)
QMC: Hammersley 点群
Hammersley 点群(Hammersley point set)も多次元の LDS で、Halton 列とほとんど同じです。 Hammersley 点群も Hammersley 列とも呼ばれるようですが、 Hammersley 点群の方がメジャーのようです。
Hammersley 点群は、最初の要素が n / N で、あとは Halton 列の構成と同じです。N は生成するサンプル数です。
Hammersley 点群は、サンプル数 N に依存するので、前もってサンプルの数が決定していなくてはならない点が、Halton 列との大きな違いです。そのため、 van der Courput 列や Halton 列では、n を 1 つずつ増やしてインクリメンタルにサンプル点を求めることができますが、 Hammersley 点群はできません。
その代わり、 Hammersley 点群は、 Halton 列よりも少し低い食い違い度(Discrepancy, 正確には Star Discrepancy) を持ちます。食い違い度が低いということは、より一様にサンプル点が散らばるということです。
プログラム
Hammersley 点群を求めるコードです。 x には n 個分のサンプル点が一度に返されます。
#includedouble vdC(int i, int base) { double h=0., f, factor; f=factor=1./(double)base; while (i>0) { h += (double)(i%base) * factor; i /= base; factor *= f; } return h; } void hammersley(double *x, int n, int dim) { int primes[10] = {2, 3, 5, 7, 9, 11, 13, 17, 19, 23}; int i, j; if (dim >= 10) { fprintf(stderr, "too big dimension\n"); return; } for (j = 0; j < n; j++) { x[j * dim + 0] = (double)j / (double)n; for (i = 1; i < dim; i++) { x[j * dim + i] = vdC(j, primes[i - 1]); } } }
グラフ
2 次元の Hammersley 点群(N=16, b1=2)
Hammersley 点群には、 2 次元のときかつ、サンプル数 N が (2^n)^2 のとき、つまり一辺のサンプル数が 2^n で表すことができる場合、面白い性質がります。
それは、区間 [0, 1)^2 を N * N のグリッドで区切ったときに、上記の条件の場合、サンプル点は、必ずグリッドの交点上にあり、またタテヨコのグリッドの区切り線に必ず1つのみのサンプル点が存在するようになるのです。
上記の N=16, b1=2 の 2 次元の Hammersley 点群のグラフがその例で、実際にグリッドを引いてみると、
となって確認することができます。
2003年10月01日
Interleaved Sampling
コンピュータグラフィックスのレンダリング問題の多くはサンプリング問題になるかと思います。特にグローバルイルミネーション問題を解く時にはいかにして少ないサンプル数で期待する答えになるようにするかというのが重要になります。
サンプリング点の生成には、低食い違い量列などを用いた規則的なサンプリングパターンと、乱数を用いた不規則なサンプリングパターンの2種類があります。
低食い違い量列などによる決定論的で規則的なサンプルパターンだとノイズが出ない代わりに、サンプル数が少ないとエイリアシングが生じてしまう・・・
乱数による不規則なサンプルパターンだとエイリアシングをノイズに置き換えることができるが、ノイズを減らして品質を上げるにはより多くのサンプル数が必要になる・・・
と、両者一長一短があります。
これのちょうど中間らへんにあたるのが、 インターリーブサンプリング(Interleaved Sampling)です。
Interleaved Sampling
Alexander Keller and Wolfgang Heidrich, 2001
日本語に訳すと、間引きサンプリング、引き伸ばしサンプリングという感じでしょうか。
以下、論文と同様に、アンチエイリアシングのコンテキストで考えます。
通常では、ピクセルのサンプリング点は、1ピクセルの領域内(1ピクセル周期)で生成します。そのため規則的なサンプリングでは、どのピクセルも同じサンプリングパターンを持ってしまうことになります。
結果、そのためパターン的なエイリアシングが生じ易くなります。
で、Keller 氏と Heidrich 氏は提案しました。 1 x 1 ピクセル周期でやるからダメなんだ。2 x 2 ピクセルの周期でやればいいんじゃないか、と。
こうすることで、各ピクセルの隣同士ではサンプリングパターンが異なるのでパターン的なエイリアシングを避けることができます。
1ピクセルにとらわれるな、さもなくば全体を見失ってしまう。という精神ですね。
これがインターリーブサンプリングです。
至極単純な方法ですが、論文の結果の画像を見ると結構品質はいいようです。
2 x 2 ピクセルの周期でサンプリングするときでも、各ピクセルのサンプル数は 1 x 1 ピクセルの周期の時と同じにするようにします。
図の例だと、2 x 2 ピクセルの範囲では 16 個のサンプル点を生成しますが、1ピクセルには必ず 4 個のサンプル点というのは変わらないということです。
確かに直感的には品質があがりそうだということは感じとれますが、なぜ 2 x 2 周期でやると品質があがるのか、またどれくらいあがるのかということについての数学的な証明は述べられてはいないようです。考えるな、感じるんだ。って云ってくれてるのかもしれません。
なんかこの考えでいくと、 3 x 3 ピクセル周期、 4 x 4 ピクセル周期と周期を上げるほど品質があがるんでしょうか。
2003年10月08日
QMC: 列と点群
低食い違い量列(Low-Discrepancy Sequence)について調べていると、Halton は Sequence(列) と表記しているのに Hammersley は point set(点群、もしくは点集合) と表記していたりで、列と点群はどちらも同じ意味なのかそれともなにか違いがあるのかなーと疑問に思っていました。今までは文献の文脈から、１次元の場合が Sequence で 2次元以上を point set と云うのだろうと思っていたのですが、
QMC 本によると、
o サンプル数 N に依存するのが point set
o 依存しないのが sequence
であるとのこと。
つまり Hammersley 点群のように生成するサンプル数 N に依存する、そのため N を変えた場合はすべてを計算しなおさなければならないような(コンピュータでの計算に不向きな)ものが point set であり、van der Corput や Halton のように N に依存しない、そのため以前に計算した値には依存せずにインクリメンタルに新しい値を計算できる(コンピュータの計算に向いている)のが sequence と云うそうです。
そのため、LDS も、Low-Discrepancy Point Sets と Low-Discrepancy Sequence の 2 種類があります。またこれらはおのおの quasirandom point sets(準乱数点群)と quasirandom sequence(準乱数列) とも呼ばれるそうです。
2003年10月11日
QMC: 食い違い度
準モンテカルロ法(QMC) で用いられる低食い違い量列(LDS)が、どれくらい一様かつランダム的に"良く"散らばるかを測るものさしが、食い違い度(Discrepancy)になります。
食い違い度の定義には各種ありますが、最も一般的に使用される食い違い度の定義はスター食い違い度(star discrepancy)と、(極値的)食い違い度( (extreme) discrepancy )になります。
スター食い違い度
N 個のサンプル点(x1,x2,..., xn)の集合を P とします。
P = (x1, x2, ..., xn)
次元数を s とすると、おのおののサンプル点は s 次元のベクトルになります。
また、簡単のために、[0,1)^{s} の単位超立方体内の区間で考えます。つまり 1 次元であれば数直線上で [0, 1) の区間で、 2 次元であれば各辺が 1 の単位正方形でものごとを考えます。この単位超立方体の区間を I^{s} で表します。
このとき、 P のスター食い違い度 D_{N}^{*}は、
と定義されます(「スター」は"*"をスターと読むことから)。
ここで、 J* は、
になります。つまり、 J* はI^{s} 内で原点を含むあらゆる組み合わせの部分空間になります(一方のカドは原点固定で、もう一方のカドをいろいろと動かしてできる領域)。
たとえば 2 次元の場合、 [0,1)^{2} の単位正方形内に含まれ、[0, u1) x [0, u2) で表されるいろいろな組み合わせの長方形(ヨコ u1, タテ u2)になります。
ここで注意しておきたいのは、I^{s} も J* も半開区間であることです(右側が閉じていない)。このことはこの後に J* の面積を考えるときに重要になります。
sup は、上限(supremum)を表し、sup f(x) は f(x) が取りうる値のうちで最大(上限)のものを定義します。 sup はほとんど最大値を表す max と似ていますが、max が定義できないものにも、sup が定義できるのがあります。
たとえば、f(x)=-1/x^2 の関数では、f(x)<0 となりますが、f(x) はゼロにはならないので、max f(x)=0 は定義できませんが、取り得る値の上限が0ということで、sup f(x)=0 は定義することができます。
A(J*, P)は、空間 I^s に散らばっている点群 P のうち、部分空間 J* に含まれる点の個数を表します。 これを 2 次元の場合で図示すると以下のようになります。
短い破線の部分空間では点が 3 つ含まれているので A = 3、長い破線の部分空間では点が 2 つ含まれれいるので A = 2 になります。
\lambda(J*) は J* のルベーグ測度(ルベーグの意味での面積)になります。今回ではこれはつまりJ*の面積と同じです。
なぜここで普通に面積 |J*|としないでルベーグ測度で lambda(J*) とするのでしょうか。それは、J* が半開区間であるからです。
結論としては今回の場合では、この半開区間の範囲のルベーグの意味での面積は今までの考えでいう面積と同じであるのですが、より厳密にこの半開区間の面積(ルベーグ測度)を考えておくことにします。
我々が小学校で習った意味での面積は、どちらも閉区間の閉じた範囲において与えられます。
図は閉じている
たとえば2次元で考えると上図のようになります。
[0,u1] x [0,u2] で囲まれる面積は、ヨコの長さ u1 、タテの長さ u2 の積、u1 u2 になります。しかし、今回の場合 J* の張る区間の片方が開区間で閉じていません。
図は上と右が閉じていない
そのため、今までの考えでの面積を厳密には与えることができません( [0, u1) x [0,u2) の面積を定義することができない)。
そこで、ルベーグ測度という概念で面積というものを考えると、このような閉じていない空間においても面積を与えることができるようになるのです。そして今回の場合はこのルベーグ測度はいつも通りの面積と同じになります。
よって、J* のルベーグ測度というのは、J* が閉区間の場合の面積 |J*| と同等になります。
最終的にまとめると、スター食い違い度とは、
I^{s} に含まれるすべての取り得る組み合わせの部分空間 J* のうちで、A を N で割った値と J* の面積 |J*| の差が最も大きくなる(上限)のときの値、
ということになります。
極値的食い違い度
スター食い違い度では、J* 区間には必ず原点が含まれていましたが、この制約をなくし、一般に区間を自由にとった場合、これを(極値的)食い違い度( (extreme) discrepancy)と呼びます。
(極値的)食い違い度 D_{N} は以下のように定義されます。
ここで、 J は [u11, u21) x [u12, u22) x ... x [u1s, u2s) \in I^{s} になります。
この J での 2 次元での A を図示すると以下のようになります。
2003年10月07日
QMC本
準モンテカルロ法(QMC)をやる上でたぶんバイブルであろう、
Harald Niederreiter
"Random Number Generation and Quasi-Monte Carlo Methods"
Cbms-Nsf Regional Conference Series in Applied Mathematics, No 63, 1992.
がやっと届きました。なのでこれからはもうちょっと詳しく QMC について説明できるようになると思います。証明とかちゃんと載っていますし。とりあえず次回の QMC ネタでは食い違い度の定義と Koksma の不等式の証明をやります。
2003年11月19日
準モンテカルロ法とインポータンスサンプリング
QMC とインポータンスサンプリングとの組み合わせについての論文をメトロポリス光輸送とヘミスフィアラジオシティを調べているときに見つけました。
まだちょこっと読んだだけですが、準モンテカルロ法(QMC)での LDS (準乱数) は、基本的にモンテカルロ法での乱数と置き換えることで、そのままインポータンスサンプリング(importance sampling)にも利用できると思っていたのですが、どうもそう簡単にはいかないようです。
Laszlo Szirmay-Kalos, Balazs Csebfalvi, Werner Purgathofer
Importance Driven Quasi-Random Walk Solution of the Rendering Equation
WSCG '98 (Sixth European Conference in Central Europe on Computer Graphics and Visualization), pp. 379--385, 1998.
pdf download(ResearchIndex)
以下も同著者による似たような論文でレンダリング方程式と準モンテカルロ法について。
Laszlo Szirmay-Kalos, Tibor Foris and Werner Purgathofer,
Quasi-Monte Carlo Global Light Tracing with Infinite Number of Rays
WSCG '98 (Sixth European Conference in Central Europe on Computer Graphics and Visualization), pp. 386--393, 1998
pdf download(ResearchIndex)
2003年12月07日
QMC: Koksma の不等式の準備
準モンテカルロ法で用いられる低食い違い量列(LDS) は、誤差は食い違い度で測られ、これはモンテカルロ法で用いられる乱数とは異なり、確定的な誤差の上限(deterministic error bound)を持つことが知られています。
乱数では、そのモンテカルロ法への利用による誤差は分散(variance)で与えられ、これは 1 / Sqrt(N) になることが知られています。しかしこれは確率的な誤差の上限(probabilistic error bound)を持ちます。
つまり、 LDS では誤差上限の値は決定しているのでそれを超える状況は起こりませんが、乱数では誤差上限の値は平均としての値であるので、必ずしもすべての状況で誤差上限内に収まることは断言できないということです。
たとえば乱数でのこの例を挙げると、サイコロの目は平均して(期待値は) 3.5 になります。乱数を用いる場合では、 N 回サイコロを振れば 1 / Sqrt(N) のオーダーですべての試行の平均が 3.5 に近づいていくことになりますが、もし N 回振ったときのサイコロの目がすべてゼロだった場合の状況を考えてみると判るかと思います。いくら N を大きくしていこうとも、このような状況に遭遇した場合(実際には起きにくいですが)、その時の誤差は最も最悪の値になります。
LDS の定理 1
LDS を始めとして、 QMC で用いられる数列では、乱数の時とはことなり、確定的な誤差上限が与えられます。これは 1 次元では Koksma の不等式(Koksma's inequality)、多次元では Koksma-Hlawka の不等式(Koksma-Hlawka inequality) により証明されます。
( Koksma はコクスマ、 Hlawka はラウカって読むのかなぁ...)
まずは Koksma の不等式を示したいのですが、それを証明するのに LDS の定理を 1 つ証明しておく必要がありますので、まずそちらを先に示したいと思います。
基本的には QMC 本の受け売りなのですが、 QMC 本では証明は簡潔に述べられているので、個人的に理解する上で必要と感じた分の加筆とかを加えています。
で、さらにその前に、理解しやすくするために、まず食い違い度のグラフを示しておきます。
例として、 Halton 列 を取り上げます。 N = 8 のときの Halton 列をグラフで表示すると以下になります。
食い違い度 (スター食い違い度)は以下で与えられます。
この点列に対して、D*の定義の絶対値の中の A/N のグラフは以下のようになります。
これはつまり x を 0 から 1 まで動かしたときに、 [0, x) に含まれる点の個数を表します。
絶対値の中のもう一つの λ は、これは 1 次元では y = x のグラフと同じことです。
これらグラフを 1 つにまとめると、
となります。 D* はこの 2 つの線の差が最も大きくなる部分を持って与えられます。つまり両方の線が一致するようになるほど食い違い度が低いということになります。
例えば、以下のグラフ、
では、差が大きくなっています。これは点群が一様に散らばらずに偏っているためです。この場合食い違い度は大きなります。
定理 1
N 個のサンプル点からなる低食い違い量列 P={x_{1}, x_{2}, ..., x_{N}} があり、
0 <= x_{1} <= x_{2} <= ... <= x_{N} <= 1
であるならば、
が成り立つ。
証明
D* は連続関数であるので(この証明はホントはまた必要なのですが、今回は自明であるとします)、
0 < x_{1} < x_{2} < ... < x_{N} < 1
と仮定できる。 x_{0} = 0, x_{N+1} = 1 とする。
ここで D* の定義は、連続関数であるため [x_{n}, x_{n+1}] の区間へ分割して考えても同じなので(上記のグラフで各点で y 軸方向に線を引いて短冊に区切って考える)、
ここで x_{n} < u <= x_{n+1} の区間では、 定義から u 個の点が含まれるので、 A = n になります。
ここで n/N - u は線形なので、上限は u = x_{n} か u = x_{n+1} のどちらかになります(グラフでいうと短冊の角のどちらかで絶対値が最大を取る)。
ここで max の中身を並べて考え、 x_{0}=0, x_{N+1}=1 を考慮すると、
n=0 : (0/N - 0, 0/N - x_{1})
n=1 : (1/N - x_{1}, 1/N - x_{2})
n=2 : (2/N - x_{2}, 2/N - x_{3})
...
n=N : (N/N - x_{N}, N/N - 1)
ここで n=0 のときの最初と n= N のときの最後にはゼロができるので、2 番目の項をそれぞれ下に1 つづつずらすと n = 0 の行が消え、
n=1 : (1/N - x_{1}, 0/N - x_{1})
n=2 : (2/N - x_{2}, 1/N - x_{2})
...
n=N : (N/N - x_{N}, N-1/N - x_{N})
になります。よって、
。
ここで、以下は N=3, n=2 のときの |n/N - x|, |n-1/N - x| (0<= x < 1) のグラフです。
これを見ると判るように、|n/N - x| と |n-1/N - x| のグラフは、 n/N と n-1/N の中間 2n-1/(2N)で交わり、ここで両式の max は最小値を取り、その時の値は 1/(2N) になります(上記のグラフでは N=3 なので 1/6)。これはどの区間の n についても当てはり、また丁度最小値の部分で左右対称となっているので、 max(|n/N - x_{n}|, |n-1/N - x_{n}|) = 1/(2N) + max(|2n-1/(2N) - x_{n}|) と書き直す事ができます。よって最終的に、
を示すことができます■
またこの定理が示すことは、 D* は 1/(2N) よりは小さくはならないという確定的な下限を持つということです。
この定理の多次元の証明は、 QMC 本によると簡単にはいかないようです。詳細は、
H.Niederreiter,
Discrepancy and convex programming,
Ann. Mat. Pura Appl., 93(1973), pp. 89-97.
L. de Clerck,
A method for exact caculation of stardiscrepancy of plane sets applied to the sequences of Hammersley,
Monatsh. Math., 101(1986), pp. 261-278.
を参照して欲しいとのこと。後者の de Clerck の論文では 2 次元の場合での明示的な公式を与えているそうです。
2003年12月10日
Latin Supercube Sampling
ラテン超方格法(Latin Hypercube Sampling, LHS, LHC 法, ラテン・ハイパーキューブ法)を調べているときに、そういえば libseq には、 さらに Latin Supercube Sampling という手法のコードがあったなぁということで、 ラテンスーパーキューブ法の論文を調べてみました。
Art B. Owen
Latin Supercube Sampling for Very High Dimensional Simulations
ACM Transaction on Modeling and Computer Simulation, vol 8, num 1, pp. 71-102, 1998.
これは、 LHS と QMC を組み合わせたもので、特にファイナンスや光輸送などの高次元の積分が必要なコンテキストで使えるとのこと。
以下 Abstract 日本語訳
本論文では、光輸送やファイナンス、待ち行列(queueing)などで生じる高次元のシミューレションのための、ラテン・スーパキューブサンプリング(LSS)を提案します。
LSS は広く利用されている 2 つの手法である、ラテン・ハイパーキューブサンプリング(LHS)と準モンテカルロ法(QMC)の組み合わせとして開発されました。
LSS では、入力の変数は部分集合へとグループ化し、そのおのおのの
部分集合に低次元の QMC 法を適用します。
QMC (サンプル)点は、この部分集合内でランダムな順番で現れます。
QMC 法は高次元の問題では効率が失われることが観測されてきました。
本論文では、もし入力の変数をうまくグループ化することができれば、LSS は QMC の利点を拡張することができることを示します。
例題を有用に見せるため、変数をグループ化するいくつかの提案を行います。
悪いグループ化でも、それでも LHS と同じくらいの結果を期待できます。
本論文ではまた、 QMC 法、その乱列化版(RQMC)、および QMC を高次元へと拡張する既存の手法のサーベイ(調査)を行います。
さらに、 RQMC に LSS を適用したものは、 LSS と QMC よりもより信頼度が高くなることを示します。
なんか LSS 結構凄そうだ... 乱列化 QMC(Randomized QMC) ももっとよく調べておきたいです。
ラテン超方格法は、これは層別サンプリングが多次元になるとサンプル数が増えてしまうのでこれを減らすという手法なのですが、これについてはまたいずれ良く調べてから書く事にしたいと思います。
2003年12月11日
Efficient Bidirectional Path Tracing by Randomized Quasi-Monte Carlo Integration
乱列化 QMC を用いて双方向パストレーシングを行う論文。
Tomas Kolling and Alexander Keller,
Efficient Bidirectional Path Tracing by Randomized Quasi-Monte Carlo Integration
K.-T. Fang, F.J. Hickernell, and H. Niederreiter (eds.), Monte Carlo and Quasi-Monte Carlo Methods 2000, Springer-Verlag, Berlin, pp. 290-305, 2000
QMC をそのまま多次元の問題である双方向パストレーシングに適用するとエイリアシングやノイズが生じることがあるので、 RQMC を用いるとより良い結果が得られることが示されています。
以下 abstract の日本語訳
コンピュータグラフィックスの用途では、通常非積分関数が非有界な変動(unbounded variation)を持つため、モンテカルロ積分に対して準モンテカルロ法(quasi-Monte Carlo Method)はサンプルからの誤差の推定を行う事ができず、また決定的な誤差の上限も得る事ができません。
我々は乱列化準モンテカルロ(randomized quasi-Monte Carlo)積分を応用し、低食い違い量サンプリング(low-discrepancy sampling)の利点を生かし、また同時に分散の推定を行うこともできるより効率的なアルゴリズムによる双方向パストレーシングを提案します。
QMC では Koksma-Hlawka の不等式により、関数が(Hardy-Klause の意味での)有界変動を持てば誤差の確定的な上限が得られるのですが、 パストレーシングでは測度寄与関数(measurement contribution function)が未知の不連続関数であったりするので非有界(unbounded)になるのでよくないね、ということが述べられています。それでも実際に既存研究(といってもこれも Keller 氏によるもの)では、 QMC でレンダリングしたものは MC よりも良い結果が得られることが示されているとのこと。
ここで非有界になるような状況があるのかと不思議に思いました。そのような状況になるのはたぶん無限大のある関数が存在するときぐらいだと思うからです。測度寄与関数に含まれるのは、 BRDF f, 可視関数 V, インポータンス関数 W が主なものになります。
そのなかで無限大を取る関数があるのかと思ったのですが、鏡などの理想鏡面反射や理想透過・屈折がそれにあたりますね。その場合 BRDF はディラックのデルタ関数(ある一点では無限大をとり、それ以外ではゼロ、積分すると 1 もしくは定数になるという積分するときに役立つ関数。普通の関数の形式としては表現できないのでこれは超関数と呼ばれる)を含むので、非有界になります。
また有界であったとしても、可視関数はゼロか 1 かなので、これにより変動が大きくりますし、またインポータンス関数 W もほとんど至る所でゼロになる関数なのでさらに変動が大きくなる可能性があります。
LSS(ラテンスーパーキューブ) についても言及されており、 LSS ではエイリアシングが生じるような場合があるとのこと。
本論文での RQMC 法は、通常の MC 法に比べて 2 から 5 倍収束が早くなることが示されています。
それにしても QMC の分野の発展もレンダリング分野並みに速い感じがします。金融が引っ張ってるからかな。
QMC で使われている部分にしか足をつっこんでいませんが数論も奥が深いです。
それにつられて mathworld で数学のコトを調べることが多くなりましたが、ホントこのサイト凄いです。レンダリングアルゴリズムに関するこのようなサイトをつくってみたいなぁ。
PARI/GP
PARI/GP は数論にフォーカスを当てた数学ソフトです。
http://pari.math.u-bordeaux.fr/
Mathematica に比べてあまり知られていないものかと思っていましたが、調べると結構日本語での解説サイトとかあります。
楕円曲線が高速に計算できるとのことですので、暗号分野とかでも使われているのかな。
LDS の検算や検証につかえそうです。
2003年12月12日
Towards real-time procing of complex financial derivatives
金融分野における、高次元(1500次元程度)での QMC の有用性を示す有名な論文。
NINOMYA, S. AND TEZUKA, S.
Towards real-time procing of complex financial derivatives
Applied Mathematical Finance 3 (1996), pp. 1-20, 1996
QMC を利用する事で、MC よりも 100 倍から 1000 倍のオーダで収束が早くなり、テストの例では 1 ベースポイント(=0.0001)精度の評価を行うのに MC では 2 日半かかる計算が QMC では 1 分で計算できることが示されています。
ただしこれは、 LDS として一般化 Niederreiter 列(generalized Niederreiter sequence)の一部である一般化 Sobol' 列(generalized Sobol' sequence)および一般化 Faure 列(generalized Faure sequence)を用いた場合です。
クラシックな LDS、たとえば Halton, Sobol', Faure 列などを用いた場合では逆に MC よりも悪くなることが示されています。
金融の分野では、フロントオフィス(金融派生商品を売る窓口?)で、顧客の要望に応じて膨大な要素からなる価格シミュレーション(多次元積分)をあれこれ設定を変えてすぐさま結果を出す必要があるようなので、QMC の収束の早さに昔から注目が集まっているみたいです。リアルタイムグラフィックスよりも計算量が求められそうです。
2003年12月16日
超一様分布列, ディスクレパンシー法
手塚集先生の日本語の書物を読んでみたところ、
discrepancy はディスクレパンシー、
low-discrepancy sequence は超一様分布列、
LDS を用いた QMC はディスクレパンシー法
と訳されていました。
超一様分布列と云うのは、LDS が、乱数による一様分布列よりもさらに一様に分布するからだそうです。
うーむ、こっちに合わせたほうがよいのだろうか...
2003年12月17日
LDS 特許
当時(1995) MBS(モーゲート証券) のプライシングを LDS を用いることで 100 倍以上高速に計算できるソフトウェア(FinDer)が金融の世界で大きな話題になったそうです。
で、なにが問題かというと、金融における LDS の利用に対して、この FinDer の開発を行ったブラウン大学の J. F. Traub 氏らにより特許が日米で申請されており、メリケンではすでに特許が成立されていることです。
LDS を用いて金融のプライシングを行うことのすべてが特許の対象になっているようです。
http://www1.cs.columbia.edu/~traub/html/body_patent_information.html
というわけで金融分野では LDS を利用したソフトはうかつに公開できない模様(とくに 一般化 Faure 列を使うものとか)。
CG レンダリングへの適用の場合はどうなんでしょうか?
2003年12月20日
日本語での QMC レンダリング論文
日本の学会での、QMC レンダリングに関する日本語での論文で有名そうなもの。
土井 章男 , 山口 美春 , 千葉 則茂
Low Discrepancy数列のコンピュータグラフィックスへの適用
情報処理学会論文誌, Vol. 38, No.7, pp. 1328-1335
大渕 竜太郎，青野 雅樹
Low-discrepancy sequenceを用いた準モンテカルロレンダリング
情報処理学会 グラフィクスとCAD研究会，1996年8月，Vol. 96, No. 77, pp. 91-96.
前者は 2 次元平面のサンプリングが主な話題、後者は一次レイがヒットしたシェーディング点での直接光のサンプリングが主な話題です。
どちらも多次元(複数反射)での場合は考慮されてはいません。
2004年01月17日
QMC テストレンダリング1
まず手始めに直接光のサンプリングを QMC で行うようにしました。
QMC の結果。各シェーディング点あたり 32 サンプル。
さすがに 32 サンプルだとサンプル数が足りないと思うので、床の影はエイリアシングのようになっています。しかしそれ以外はノイズがなくていい感じです。
MC の結果。各シェーディング点あたり 48 サンプル。
サンプル数が QMC と MC では違います。 QMC では 一般化スクランブル Hammerlsey 点集合(generalized scrambled Hammersley point set)を用いているので、サンプル数が 2^n 個の時が一番層別化される(stratified)ためです。 MC では、仰角のサンプル数は方位角の 3 倍を割り当てて層別サンプリングを行っています。そのためサンプル数は 3 * n^2 になります。
MC では床はノイズが出てソフトな感じがするので床だけ見れば QMC よりもいい感じですが、全体的にノイズが出てしまっています。
今回はディフューズ面のシーンで定数色のドームライトからのライティングであり、またそれほど遮蔽が多いシーンではありません。そのため QMC は MC よりもサンプル数が少ないにもかかわらずなかなかいい結果を出しています。
もうすこし照明分布や BRDF の変化が大きかったり、遮蔽が多いシーンでも QMC が有効かどうかをテストする必要があります。
コード
サンプル点を求めるコードは以下のようになります。
extern double generalized_scrambled_halton(int i, int offset,
int dim, int **p);
extern double generalized_scrambled_hammersley(int i, int offset,
int n, int dim, int **p);
...
for (i = 0; i < nsamples; i++) {
samples[2 * i + 0] = generalized_scrambled_hammersley(
i, 0, nsamples, 1,
perm);
samples[2 * i + 1] = generalized_scrambled_hammersley(
i, 0, nsamples, 2,
perm);
u = generalized_scrambled_halton(
ray->i, 0, 3,
perm);
v = generalized_scrambled_halton(
ray->i, 0, 4,
perm);
samples[2 * i + 0] = mod_1(u + samples[2 * i + 0]);
samples[2 * i + 1] = mod_1(v + samples[2 * i + 1]);
}
/* x = x mod 1 */
static double
mod_1(double x)
{
return (x - floor(x));
}
まず、1 次元と 2 次元の一般化スクランブル Hammersley 点集合で局所的なサンプル点を求めます。 perm は Faure の置換列 ( generalized van der Corput 列 参照)です。 generalized_scrambled_hammersley() と generalized_scrambled_halton() については実装のコード qmc.c を参照してください(こういうときに Viewcvs が便利ですね)。
次の一般化スクランブル Halton 列は、大域的な積分で考えたときのサンプル点（u,v)の生成です。1, 2 次元はスクリーンのサブピクセル位置にすでに用いていますので、 3, 4 次元を用います。 ray->i はレイのインスタンス番号です。これは低食い違い量列による関連度(correlation)を解決するために用いられます。(u,v)は"オフセット"として考えることができます。
最終的なサンプル点は、この (u, v) を足して、 1 でモジュロを取る ( [0, 1) の区間に納めるため)ことで生成されます。
あとはこのサンプル点を球面座標に変換してレイの方向を求めれば OK です。
2004年01月15日
100 次元までの素数
落ちつけ…………心を平静にして考えるんだ…こんな時どうするか……
２…３、５…７…１１…１３…１７……１９。
落ちつくんだ…『素数』を数えて落ちつくんだ…
『素数』は１と自分の数でしか割ることのできない孤独な数字……
わたしに勇気を与えてくれる。
Enrico Pucci. 1972-2011(?)
QMC の論文をいくつか読み、 QMC をどのように多次元のレンダリング問題に適用すればよいかが分ってきました。
これから lucille を全面的に QMC に移行しようと思います。
19 日のリリース( 0.1.3 ) には間に合いませんが、その次のリリース ( 0.2 ) に含めようと思います。
0.2 では以下の内容を予定しています。
o QMC への全面的移行
o RenderMan シェーダ言語の搭載(前回やった シェーダ言語の実装 の lucille への統合)
o モデラーへの対応(候補: Maya, Blender など)
次元数
レンダリングの積分問題では、実際どれくらいの次元が必要になるか考えてみました。
次元というのは、つまりは MC でいうところの乱数が必要なイベントがいくつあるかということになります。
レイが各面で反射するときに、まず光源のサンプルとレイが跳ね返る方向のサンプルが必要になります。
光源が面光源だとすると、平面のサンプルであるので 2 次元、方向のサンプルには仰角と方位角で２次元、計 4 次元が必要です。
ここにロシアンルーレットの利用を含めると 1 次元の追加になります。
他にも BRDF のサンプルとかいろいろ考慮しても、だいたい１回の反射で 5 - 8 次元程度でしょう。
現実的なシーンでは、レイの反射は多くても 10 回もあれば十分だと思います。
(鏡が向かい合わせになっていて無限に反射があるシーンなどの極端な場合を除く)
これに加えて、最初のスクリーンピクセルのサンプルに 2 次元、モーションブラーがある場合は時間のサンプルに 1 次元、被写界深度がある場合はレンズのサンプルに 2 次元が追加になります。
これらをまとめると、だいたいレンダリングの積分問題は 100 次元もあれば十分ではないでしょうか。
こうして見ると、金融が数千次元の問題であることに比べれば、レンダリングの問題においては次元はそれほど憂慮することではない気がします。
さらに、Keller 博士の Strictly Deterministic Sampling によれば、軌道分割(trajectory splitting) と依存サンプリング(dependent sampling)を用いれば、実質的にはレイの各反射における面光源のサンプルや反射方向のサンプルは低次元(2, 3 次元)の問題とできるので、次元が高いことによるエイリアシングや収束の問題というよく云われている QMC の欠点というのはほぼ解決されると思います。
どちらかというと、レンダリングの非積分関数は可視関数やBRDF関数により非連続であったり非有界な関数になるので、これをどううまく対処するかが重要なようです。
いずれにせよ、実際にプログラムで組んでみてこれら理論を確かめてみないとダメですね。
Strictly Deterministic Sampling methods in computer graphics で述べられている QMC 分散レイトレーシングを参考にして lucille の分散レイトレーサを QMC 分散レイトレーサにしようと思います(グロスなどの blurry 反射をする BRDF が無ければグローバルフォトンマップ+分散レイトレーサも悪くない気がします)。Strictly Deterministic Sampling の内容を理解するにはいくつかの数学的予備知識(数論、モンテカルロ積分論)が必要になりますが、理解してしまえば実装はとても簡単にできます。
100 次元までの素数
QMC では、各次元に互いに素な(重ならない)素数を 1 つづつ割り当てていくのがよいそうです。つまり 100 次元までの QMC 積分を行いたい場合、 100 番目までの素数を求めておく必要があります。 100 次元までの素数のリストです。
int primes[] = {
2, 3, 5, 7, 11, 13, 17, 19, 23, 29, /* 10 */
31, 37, 41, 43, 47, 53, 59, 61, 67, 71, /* 20 */
73, 79, 83, 89, 97, 101, 103, 107, 109, 113, /* 30 */
127, 131, 137, 139, 149, 151, 157, 163, 167, 173, /* 40 */
179, 181, 191, 193, 197, 199, 211, 223, 227, 229, /* 50 */
233, 239, 241, 251, 257, 263, 269, 271, 277, 281, /* 60 */
283, 293, 307, 311, 313, 317, 331, 337, 347, 349, /* 70 */
353, 359, 367, 373, 379, 383, 389, 397, 401, 409, /* 80 */
419, 421, 431, 433, 439, 443, 449, 457, 461, 463, /* 90 */
467, 479, 487, 491, 499, 503, 509, 521, 523, 541 /* 100 */
};
2004年01月22日
Latin hypercube sampling
モンテカルロ法で層別サンプリング(stratified sampling)を行う場合、次元が高くなると層別化が難しくなるという問題があります。
例えば、各次元を N で分割するとすると、次元数を d とすると全部で N^d 個のサンプル点が必要になってしまいます。
逆に全体でのサンプル数を N で固定するとすると、各次元に割り当てるサンプル数をこの N を分解して求めなければなりませんが、この N の分解というのが非常に困難になります。例として、 N が素数で全部で N^d 個のサンプル点を生成するとします。各次元に割り当てるサンプル数を (N1, N2, ..., Nd) とすると、
N^d = N1 * N2 * ... * Nd
となりますが、左辺の N は素数なので、そのべきの数もまた素数 N でしか割り切ることができません。そのため右辺が N1 = N2 = .... Nd = N の時しか成り立ちません(右辺に N 以外があるとすると、その数で割り切れてしまうことになるから)。
N が素数でなくとも、N を d 個の数に分解する問題はなかなか困難です。そこでこの高次元での層別化を容易にするのが、ラテン超方格法になります。
ラテン超方格法
ラテン超方格法(Latin hypercube, ラテンハイパーキューブ法) は McKay ら [1] により提案され、 d 次元の層別化においてちょうど N 個のサンプル点を生成する手法です。
ラテン超方格法はまた、N-ルークサンプリング(N-rook sampling)、N-クイーンサンプリング(N-queen sampling)とも呼ばれます。ルークとクイーンはどちらもチェスの駒に由来します。ラテン超方格法は、 Stein [2] によりサンプル数を多く取れば必ず機能することが示されており、 また Owen [3] によりラテン超方格法の中心極限定理が示されています。
まず、 2 次元における通常の層別サンプリングの例を示します。
この例では、分割数を N=4 にしています。そのため全部で 16 個のサンプル点が生成されています。
ラテン超方格法では、各次元(各軸)の各分割において 1 つのセルのみにサンプル点を生成します。ラテン超方格法の例を示します。
ここでグレーのセルがラテン超方格法でサンプルの生成されてセルを表します。各次元の各分割において 1 つのセルのみにサンプル点を生成するので、この 2 次元の場合では各分割のタテヨコにサンプルされるセルが重複して存在することはありません。そのためどんな高次元に対しても、丁度分割数と同じ N 個のサンプル点を生成することができます。これにより無駄なサンプル点を生成することなく自由に N をコントロールすることができるという利点があります。
ランダム順列
ラテン超方格法を実装するのには、ランダム順列(random permutation)を用います。これはトランプのシャッフル操作と同じです。
まず、 N = 4 を例として、普通に数字を順にならべた列
(0, 1, 2, 3)
から初めます。これを乱数を用いて入れ替え操作を行います。これは C++ STL の random_shuffle() で行うこともできます。コードとしては以下のようになります。
void
random_permute(int n, double *x)
{
int i, j;
double tmp;
for (i = 1; i < n; i++) {
j = (int)(random() * (double)(i + 1));
/* swap index i and j. */
tmp = x[j];
x[j] = x[i];
x[i] = tmp;
}
}
random() は [0,1) の区間で乱数を返すルーチンです。これにより、どの要素にも 1 / N の確率で各数字が現れます。ランダム順列を適用すると、
(1, 3, 2, 0)
のような並びが生成されます。
サンプルセルの生成
上記のランダム順列を次元数分繰り返して、 d 個の組を生成します。今回の 2 次元の例だと、
(1, 3, 2, 0)
(2, 0, 3, 1)
のようになります。あとはこの組の各要素を縦に取っていってサンプルセルの位置とします。この例だと (1,2), (3,0), (2,3), (0,1) になります。この位置は上図のラテン超方格法でサンプルされるセルの位置とちょうど一致しているのを確認してください。
コード
ラテン超方格法のコードは以下のようになります。
void
latin_hypercube(double *x, int n, int d)
{
int i, j;
double **replications; /* d replications of n points */
double r;
replications = (double **)malloc(sizeof(double) * d);
if (!replications) {
fprintf(stderr, "mudah!\n"); exit(-1);
}
for (i = 0; i < d; i++) {
replications[i] = (double *)malloc(sizeof(double) * n);
if (!replications[i]) {
fprintf(stderr, "mudah!\n"); exit(-1);
}
/* fill with identity mapping (0, 1, ..., n) */
for (j = 0; j < n; j++) {
replications[i][j] = j;
}
/* randomly permute array */
random_permute(n, replications[i]);
}
/* Now replications become like this(here for exapmple n = 5, d = 3).
*
* replications[0] (1, 3, 4, 0, 2)
* replications[1] (0, 1, 3, 2, 4)
* replications[2] (4, 2, 0, 1, 3)
*/
/* Then, perform stratified random sampling.
* Grid space is choosen from each column of replications.
*/
for (i = 0; i < n; i++) {
for (j = 0; j < d; j++) {
r = random() / n; /* [0, 1/n) */
x[j] = (double)replications[j][i] / n + r;
}
}
for (i = 0; i < d; i++) {
free(replications[i]); replications[i] = NULL;
}
free(replications); replications = NULL;
}
[1] McKay, M. D., Beckman, R. J., and Conover, W. J., "A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code," Technometrics, Vol. 21, No. 2, 1979, pp. 239-24.
[2] M.L. Stein. Large sample properties of simulations using latin hypercube sampling. Technometrics, 29(2):143-151, 1987.
[3] A. B. Owen, A central limit theorem for Latin hypercube sampling, Journal of the Royal Statistical Society Ser. B 54(1992), 541-551.
[4] 湯前祥二、鈴木輝好. モンテカルロ法の金融工学への応用. 朝倉書店 2000. amazon
2004年01月28日
フィボナッチ格子点
Efficient Bidirectional Path Tracing by Randomized Quasi-Monte Carlo Integration でちょこっと言及されていたフィボナッチ格子点(Fibonacci lattice points)というのが気になっていたので調べてみました。
フィボナッチ格子点とは、基本的に 2 次元で生成される点で、以下のように定義されます [1] (原書が入手できていないので以下はネットで見つけたものです)。
ここで Fk は k 番目のフィボナッチ数、 { } は値の少数部を取ることを示しています。生成される点の数は Fk - 1 個になります。
フィボナッチ格子点のプロット図。
[0, 1)^2 に収まり、規則的に並んでいることが見て取れます。
例
たとえば、 k = 7 のとき、 F_{k-1} = 8, F_{k} = 13 となり、全部で 12 個のサンプル点が生成されます。
x1 = (0.000000, 0.000000)
x2 = (0.076923, 0.615385)
x3 = (0.153846, 0.230769)
x4 = (0.230769, 0.846154)
x5 = (0.307692, 0.461539)
x6 = (0.384615, 0.076923)
x7 = (0.461538, 0.692308)
x8 = (0.538462, 0.307693)
x9 = (0.615385, 0.923077)
x10 = (0.692308, 0.538462)
x11 = (0.769231, 0.153846)
x12 = (0.846154, 0.769231)
--
このフィボナッチ格子点は、積分においてのサンプル点の生成(lattice rule, 訳は格子則?)としてとして用いられたのが元々のようです。フィボナッチ数は、とても奥が深くいろいろな研究がされています。たとえば今回のと関連がありそうなものを挙げるとパッキング問題 [2] など。
しかし、このフィボナッチ格子点が LDS であるかどうかは分りませんでした。 ただ、今回とは別の、フィボナッチ多項式を用いて LDS を構成するというのは手塚先生ら [3] が開発しています。
今回のものとは異なるかもしれませんが、フィボナッチ数によるサンプル点の生成は、redqueen や Arnold でもつかわれているようです。
フィボナッチ格子点のコード
以下もしくは cvs を参照。
void
fibonacci_lattice_2D(double *x, int k)
{
int i;
int f1 = fibonacci(k - 1);
int f2 = fibonacci(k);
if (k < 3) {
fprintf(stderr, "k must be >= 3.\n");
return;
}
for (i = 0; i < f2 - 1; i++) {
x[2 * i + 0] = (double)i / (double)f2;
x[2 * i + 1] = mod_1((double)(i * f1) / (double)f2);
}
}
double
mod_1(double x)
{
double v = x - floor(x);
assert(v >= 0.0);
assert(v < 1.0);
return v;
}
int
fibonacci(int k)
{
if (k < 3) return 1;
return fibonacci(k - 1) + fibonacci(k - 2);
}
[1] I.H. Sloan and S. Joe. Lattice methods for multiple integration. Oxford Science Publications. Oxford: Clarendon Press., 1994.
[2] Andrew V. Reztsov; Ian H. Sloan. On 2D packings of cubes in the torus Proc. Amer. Math. Soc. 125 (1997), 17-26.
[3] Fast generation of low discrepancy points based on Fibonacci polynomials, Proceedings of the 24th conference on Winter simulation table of contents. Pages: 433 - 437, 1992
2004年02月09日
Bit reversal
LDS を生成する SamplePack や、論文 Efficient Multidimensional Sampling では、基数 2 の van der Corput 列の生成に以下のコードが用いられています。
inline double RI_vdC(uint bits, uint r = 0)
{
bits = ( bits << 16) | ( bits >> 16);
bits = ((bits & 0x00ff00ff) << 8) | ((bits & 0xff00ff00) >> 8);
bits = ((bits & 0x0f0f0f0f) << 4) | ((bits & 0xf0f0f0f0) >> 4);
bits = ((bits & 0x33333333) << 2) | ((bits & 0xcccccccc) >> 2);
bits = ((bits & 0x55555555) << 1) | ((bits & 0xaaaaaaaa) >> 1);
bits ^= r;
return (double)bits / (double)0x100000000LL;
}
とまあ摩訶不思議なビット操作をするコードなのですが、 Hacker's Delight をぱらぱらと見ていたら、これとビンゴなコードが。
Hacker's Delight によれば、上記コードの最初の 6 行のビット演算の部分は、ビットの逆順を求めるコードです。
たとえば、
00000000000000000000000001000101 (2 進)
というビット列があった場合、この逆順(左右ひっくり返したもの)は以下のようになります。
10100010000000000000000000000000 (2 進)
今回は 32 bit の操作ですが、Hacker's Delight には 32bit 以外のビット幅での操作のコードもあります。
このコードの応用としては、エンディアンの入れ替えや FFT があるそうです。このビット逆順のテクニックは、REXX という言語のサブルーチンコードで見つけられるとのこと。やっていることは、まず最初に隣り合う 1 ビット同士を入れ替え、次に範囲を広げて隣り合う 2 ビット同士を入れ替え、・・・ということです。上記の RI_vdC では Hacker's Delight のコードとは少しことなります(RI_vdC ではシフト量の多い16 ビットシフトが先)が、出てくる結果は同じです。
ビット逆順のコードは約 30 インストラクションで処理でき、また分岐・除算なしです。Hacker's Delight ではさらに細かく、回転シフト(rotate shift, PowerPC や Java の <<< など) を使えば 1 インストラクション減らすことができるとか。Hacker's Delight すごいなぁ... 深く読めばいろいろ有益なものが見つかりそうです。
がんばればレイ・三角形の交差判定もすべてビット演算でできるかなぁ..
vdC 例
van der Corput 列の定義 の性質から、このビット逆順が非常に基数 2 の van der Corput 列の生成に向いており、効率的に van der Corput 列を求める事ができることが納得できると思います。
たとえば、 i = 1 のときを考えてみると、vdC(1) = 0.5 になります(基数 2 のとき)。これを上記の RI_vdC() で見ていくと、まず
00000000000000000000000000000001 (2 進) = 1
のビット逆順は、
10000000000000000000000000000000 (2 進) = 2147483648(unsigned int)
になります。
r はランダムを加える要素なので、今回は r=0 として考えて無視します。
0x100000000LL (16 進) = 最初のビットが 1、他がすべてゼロの 33 bit 値 = 4294967296
なので、
(double)2147483648 / (double)4294967296 = 0.5
となるのが確かめられます。
最初にこのコードを SamplePack や Efficient Multidimensional Sampling で見た時は、Kollig 氏か Keller 博士が独自に考えたものでそのままソースコードに組み込むのはまずいのかなぁ...と思っていたのですが、これで lucille に組み込んでも問題なさそうですね。しかし基数 2 の場合限定なのでそれほどこのビット演算版が効率的になるかはわかりません。ビット数は 32bit なので精度が問題になのでは? と思うかもしれませんが、レンダリングにおいてインスタンス変数 i が 42 億になるような状況は今のところ思い浮かばないのでたぶん問題ないのではないかと思います。しかし金融のようにサンプル点を多く生成することが必要な状況だと、もしかしたら問題になるかもしれません。
Mersenne Twister でもビット演算が多く用いられています。結構数論関係はビット演算で効率的に演算できる部分というのが多いのかな。
2004年02月22日
有限体
多くの低食い違い量列は、 Faure 列( Faure sequence) を元に生成されてきましたが、 Niederreiter 氏と Xing 氏により、有限体(finite field)に基づいた全く新しい低食い違い列である Niederreiter-Xing 列(Niederreiter-Xing sequence, NX 列) が考案されました。
C. Xing and H. Niederreiter,
A construction of low-discrepancy sequences using global function fields,
Acta Arith., 73 (1995), pp. 87-102.
NX 列は、既存の LDS の中でもかなり優秀なもののようで、それなりの高次元でも使えるようです。
NX 列の実装は [1] で入手することができます。
有限体
この NX 列の理論の元になっている有限体とは、足したり掛けたりの演算を行っても有限の位数(要素数)しか取りえない世界(体)のことです。有限で閉じているので、コンピュータでの演算に適しており、実際有限体は楕円曲線、暗号、符号理論、低食い違い量列などと幅広く応用されているようです。
Niederreiter 氏と Xing 氏は有限体についても多くの出版物を出しており、低食い違い量列への有限体の応用の解説も含んでいるものとしては、いずれ理解して役に立つだろうということで購入した、
Rational Points on Curves over Finite Fields: Theory and Applications (London Mathematical Society Lecture Note Series, N285/000)
Harald Niederreiter , Chaoping Xing. Cambridge Univ Pr , ISBN: 0521665434 , 2001.
が挙げられます。
現在は読んでもちんぷんかんぷんですが、有限体は理解しているとなんかいろいろ後々役に立ちそうです。
それにしても Niederreiter 博士はいったいいくつの論文と書物を書いているんだ!!!
[1] G. Pirsic, `A Software Implementation of Niederreiter-Xing Sequences', a revised/corrected version appeared in the Proceedings of MCQMC 2000.
2004年08月21日
計算統計 I -- 確率計算の新しい手法
汪 金芳 (著), 手塚 集 (著), 上田 修功 (著), 田栗 正章 (著), 樺島 祥介 (著), 甘利 俊一 (編集), 竹村 彰通 (編集), 竹内 啓 (編集), 伊庭 幸人(編集)
計算統計 I—確率計算の新しい手法 統計科学のフロンティア 11
岩波書店, ISBN: 4000068512, 2003.
amazon.co.jp
伊庭幸人先生 + ベイズ統計がらみで検索していたところ発見した書籍。
本のタイトルからは想像できかったので知らなかったのですが、実はこの本の第二章は手塚先生による超一様分布列(低食い違い量列, low-discrepancy sequences, LDS)の解説です。
日本語による LDS の解説がほとんど無いなかで、本書は日本語による LDS の貴重な文献であると思います。
(まあ 湯前本 を除くと、書籍になっている LDS に関する日本語の解説のほとんどは手塚先生によるものなのですが...)
とくに Sobol' のネット理論(Sobol's net theory)、Owen's randomization(Randomized QMC)や一般化 Niderreiter 列についても解説されています。
これは日本語では初めて解説されるものではないでしょうか?
(欲を言えば有限体理論による新しい LDS である NX 列(Niderreiter-Xing sequence) も取り上げてほしかった...)
また、なぜ手塚先生は LDS を超一様分布列と訳しているのか など、LDS に関する用語の日本語訳の問題についても取り上げており、これも LDS の歴史を知る上でとても貴重で参考となるものです。
超一様分布列は、かつて「準乱数」とよばれていたものとほとんど同じものである。 なぜ呼び方が変わったのかについて説明することが、この分野の歴史と背景を述べ ることにもつながるため、この話から入ろう。 ... 今日のようにコンピュータがこれだけ普及し、「乱数」もさまざまな分野で応用さ れるようになると、準乱数と擬似乱数の違いを明確にするために、欧米では、 low-discrepancy sequences という用語が今日広く使われている。そして、 それに呼応して日本でも、その訳語として「超一様分布列」が用いられるようにな ってきている。 ではなぜ、「low-discrepancy」の訳が「超一様分布」なのか?これには少しこの 分野の歴史から説明することが必要となる。 ... そこでは、分布の不規則性、すなわち一様分布からのズレを測る尺度としてディス クレパンシーという量が定義され、ディスクレパンシーができるだけ小さい、つま り「一様性のできるだけ高い」有限点列を構成することが主要な研究テーマとなっ た。こうした理由から、「low-discrepancy」と「超一様分布」が結びつくので ある。(本書より一部引用)
うーむ、しかしそのような背景説明なくして、超一様分布列というのは low-discrepancy sequences である、というのを想像するのは難しいので、私のような初学者などにはこの訳は少々取っ付きにくい感じがします(まだ低食い違い量列(低い食い違い量を持つ列)や低齟齬列のほうがなんとなく想像できる)。
さて、本書を読んで目から鱗であったのは、LDS(準乱数)のデザインには「乱数らしい」性質を持つことは一切要求されていないということでした。
てっきりランダムに散らばった感じがしつつも一様に分布するものが LDS だと思っていましたが、本書によれば LDS は数値積分の収束を早くなることを目的としてデザインされた点列であるだけであるとのこと。
擬似乱数も実際にはランダムではなくてコンピュータで生成される「乱数らしさ」を持つ決定的に生成される点列です。
そういう意味ではコンピュータ上では擬似乱数も LDS も決定的な点列であり、ただその目的が異なっているだけことなっているということです。
(ただ、では「(擬似)乱数」を用いるモンテカルロ法と、「LDS」を用いる準モンテカルロ法も同様に同じ理論なのかと言うとそうではなく、前者は確率論で論ぜられ、後者は数論で論ぜられるのでそこは区別が必要だと思います。)
本書の第一章、第三章もベイズ統計がらみで参考になりそうです。
このシリーズは全 12 巻構成。全部買ったら大変だ...