低差異序列（一）- 常見序列的定義及性質

01-27

高效的生成在高維空間分布均勻的隨機數是在計算機程序中非常常見的組成部分。對於一切需要採樣的演算法來說，分布均勻的隨機數就意味著更加優秀的樣本分布。光線傳遞的模擬（渲染）基於蒙特卡洛積分（Monte Carlo Integration），這個過程中採樣無處不在，所以好的樣本分布直接影響積分過程的收斂速度。與常見的偽隨機數對比，低差異序列（Low Discrepancy Sequence）非常廣泛的被用在圖形，甚至於金融領域。它們除了在高維空間中的分布更加均勻以外還有許多其他的性質更利於渲染程序的執行。本文標題圖片的左右兩邊分別用32個Sobol序列和偽隨機數作為樣本分布渲染，可以看出左邊的噪點比右邊少許多。這篇專欄會介紹幾種常見的低差異序列的定義（Van der Corput， Halton，Hammersley，Sobol，Rank-1 lattice），下一篇專欄（低差異序列（二）- 高效實現以及應用）則會專註於它們高效的實現，以及在渲染程序應用時候的一些問題。

什麼是Discrepancy

首先說說這裡均勻分布里的「均勻」指的是什麼。一個直觀的理解可以看下面的圖片，左邊為偽隨機數組成的二維點集，右邊則是由低差異序列點集的對整個空間的覆蓋更加完整。

更加嚴謹的定義則要引入Discrepancy（Definition of Discrepancy）的概念

$D_N(P) = sup_{Bin J}n left| frac{A(B)}{N} - lambda_s(B) right|$

公式看的眼花繚亂，用普通話描述一下就是，對於一個在 $[0,1]^n$ 空間中的點集，任意選取一個空間中的區域 $B$ ，此區域內點的數量 $A$ 和點集個數的總量 $N$ 的比值和此區域的體積 $lambda_s$ 的差的絕對值的最大值，就是這個點集的Discrepancy。分布越均勻的點集，任意區域內的點集數量占點總數量的比例也會越接近於這個區域的體積。

Radical Inversion與Van der Corput序列

接下來在介紹這些序列定義之前，先介紹一個基本的運算，Radical Inversion。這篇專欄將會介紹的所有序列都會用到這個運算過程。

$i=sum_{l=0}^{M-1}{a_l(i)b^l}$

$Phi_{b, C}(i)=(b^{-1}...b^{-M})left[ C left( a_0(i)...a_{M-1}(i) right)^T right]$

這個操作非常直觀， $b$ 是一個正整數，則任何一個整數 $i$ 如果先將 $i$ 表示成 $b$ 進位的數，然後把得到的數中的每一個位上的數字 $a_l(i)$ 排成一個向量，和一個生成矩陣（generator matrix） $C$ 相乘得到一個新的向量，最後再把這個新向量中鏡像到小數點右邊去就能得到這個數以 $b$ 為底數，以 $C$ 為生成矩陣的radical inversion $Phi_{b, C}(i)$ 。

上面的描述可能略微有些複雜，但如果矩陣 $C$ 是單位矩陣（Identity Matrix）的時候，這個過程會簡化很多，既是直接把這個 $b$ 進位的數鏡像到小數點右邊去即可。同時這也是Van der Corput序列的定義。

$Phi_{b}(i)=(b^{-1}...b^{-M})left( a_0(i)...a_{M-1}(i) right)^T=sum_{l=0}^{M-1}{a_l(i)b^{-l-1}}$

舉個例子，正整數8以2為底數的radical inverse的計算過程如下。首先算出8的2進位表示，1000。此處假設 $C$ 為單位矩陣，所以直接將1000鏡像到小數點右邊，0.0001。這個二進位數的值就是最終結果，把它轉換回10進位就得到1/16, 既 $Phi_2 (8)=1/16=0.0625$ 。下面的表給出了更多的以2為底數的Van der Corput序列的例子。

Van der Corput序列有幾個屬性：1）每一個樣本點都會落在當前已經有的點裡「最沒有被覆蓋」的區域。例如 $Phi_2 (3)=3/4$ 是剛好落在了 $[0,1)$ 區間中被覆蓋最少的區域（ $Phi_2 (0)=0,Phi_2 (1)=1/2,Phi_2 (2)=1/4$ ）。2）樣本個數到達 $b^m$ 個點時對 $[0,1)$ 會形成uniform的劃分。例如 $Phi_2 (0)=0,Phi_2 (1)=1/2,Phi_2 (2)=1/4,Phi_2(3)=3/4$ 。3）很多時候並不能夠代替偽隨機數，因為點的位置和索引有很強的關係，例如在以2為底的Van der Corput序列中，索引為基數時候序列的值大於等於0.5，偶數時則小於0.5

Halton序列與Hammersley點集

介紹完Van der Corput之後，Halton和Hammersley就非常簡單了。Halton和Hammersley可以生成在無窮維度上分布均勻的點集。它們都基於Van der Corput序列。

Halton序列的定義很簡單：

$X_i:=(Phi_{b_1}(i),...,Phi_{b_n}(i))$

既是每一個維度都是一個基於不同底數 $b_n$ 的Van der Corput序列，其中 $b_1...b_n$ 互為質數（例如第 $1$ 到第 $n$ 個質數）。

Hammersley點集的定義和Halton非常類似。

$X_i:=(frac{i}{N}, Phi_{b_1}(i),...,Phi_{b_{n-1}}(i))$

唯一不同的就是將第一個維度變成 $frac{i}{N}$ ，其中 $i$ 為樣本點的索引， $N$ 為樣本點集中點的個數。根據定義，Hammersley點集只能生成固定數目個樣本，而Halton序列則可以生成無窮個樣本（當然在計算機里我們只有有限的bit去表示有限個樣本點）。

上面左邊的圖為第1-100個Halton序列中的二維的樣本點， $(Phi_2(i),Phi_3(i))_{i=0}^{99}$ ，右邊則為數量為100的二維Hammersley樣本點集， $(frac{i}{100}, Phi_2(i))_{i=0}^{99}$ 。可以看出來他們的分布都遠比一般的偽隨機數更加均勻。Hammersley的Discrepancy比Halton更稍低一些，但是代價是必須預先知道點的數量，並且一旦固定沒法更改。虛幻引擎4中對環境貼圖的Filter採樣就是用的點集大小固定為1024的Hammersley點集。Halton雖然Discrepancy稍高但可以不受限制的生成無窮多個點，更適合於沒有固定樣本個數的應用，例如任何progressive或者adaptive的過程。

基於radical inversion的序列還都具有Stratified樣本的性質。因為每一個維度都是一個radical inversion，所以每一維度都具有所有之前提到的radical inversion的性質。其中之一就是點集個數到達 $b^m$ 個點時對 $[0,1)$ 會形成uniform的劃分。下圖是第1-12個Halton序列的二維點集，可以看出點0-7在X軸的投影和0-8在Y軸的投影都是均勻覆蓋。這也意味著在樣本數量等於每個維度底數的公倍數的時候，樣本會自然在每個維度上底數的倍數的strata中自然的形成stratified採樣。例如下圖中的第0-5個點，剛好在圖中落在2x3的strata中。

Halton序列的一個缺點是，在用一些比較大的質數作為底數時，序列的分布在點的數量不那麼多的時候並不會均勻的分布，只有當點的數量接近底數的冪的時候分布才會逐漸均勻。例如下圖中以29和31為底的序列，一開始的點會分別是1/29,2/29,3/29...所以造成了點都集中落在了一條直線上面。解決這個問題的方法也很簡單，Scrambling。Scrambling的方法有許多種，例如最簡單的XOR Scrambling適用於以2為底數的序列。對於Halton來說，一個比較常用的方法是Faure Scrambling。

$Phi_b(i)=sum_{l=0}^{M-1}{sigma _b(a_l(i))b^{-l-1}}$

如上面的公式所寫，Faure Scrambling的做法就是在做radical inverse的時候不直接將數字鏡像到小數點右邊，而在鏡像前先把每個數字通過一個permutation $sigma _b$ 轉換成另一個數字。不同的底數 $b$ 有不同的permutation $sigma$ 。例如 $sigma_4=[0,2,1,3]$ 。至於 $sigma _b$ 如何具體計算這裡不再展開，下一篇專欄在講實現時會給出參考鏈接。這裡值得一提的是Scrambling完全不會影響radical inversion序列分布的隨機性，因為radical inversion會自然的將空間均等劃分成底數 $b$ 的整數次冪個部分，scrambling本質上就是在交換這些均等劃分的部分，所以Scrambled後的序列依然具有radical inversion的性質。

Sobol序列

與Halton和Hammersley不同，Sobol序列的每一個維度都是由底數為2的radical inversion組成，但每一個維度的radical inversion都有各自不同的矩陣 $C$ 。 $X_i:=(Phi_{2, C_1}(i),...,Phi_{2, C_n}(i))$

因為完全以2為底數，所以Sobol序列的生成可以直接使用bit位操作實現radical inversion，非常高效。Sobol序列的分布具有不僅均勻，而且當樣本的個數為2的整數次冪時，在 $[0,1)^s$ 區間中以2為底的每個Elementary Interval中都有且只會有一個點，這意味著它可以生成和Stratified Sampling和Latin Hypercube同樣高質量分布的樣本（見下圖），同時又不需要預先確定樣本的數量或者將樣本儲存起來，並可以根據需要生成無限個樣本，非常適合progressive的採樣。這些性質也使得Sobol在需要一切對高維空間採樣的應用中，例如圖形，渲染以及金融領域，都非常流行，

因為Sobol序列需要一個生成矩陣 $C$ ，而且所有維度都以2為底，所以沒有Halton那樣在以比較大的數為底時需要用Scrambling來消除分布間的correlation這個問題。那麼如何計算出能生成如此高質量分布的矩陣 $C$ 呢？Quasi Monte Carlo的學者們已經花了數10年的時間搜索這種矩陣，現在我們可以在這個網頁（Sobol sequence generator ）找到可以生成21201維度的Sobol序列的矩陣。

Rank-1 Lattices

最後再簡單提一下Rank-1 Lattices。類似於Sobol序列，Rank-1 Lattices依賴於一個生成向量（Generator Vector） $G=(g_0,...g_{n-1})$ 。生成向量的質量直接影響到最終樣本的分布。給定一個生產向量後，產生一個點集的做法非常簡單：

$X_i:=frac{i}{N}(g_0,...,g_{n-1}) mod [0,1)^n$

既每個樣本都是用 $frac{i}{N}$ 乘以生成向量， $i$ 為樣本的索引， $N$ 為樣本總數，得到的乘積如果大於1，則對1求餘數（modulate）映射回 $[0,1)$ 範圍。

這種做法類似Hammersley也只能用於預先知道樣本數量並且樣本數量固定的應用，如果需要可持續的生成無線個樣本，也很簡單，只需將 $frac{i}{N}$ 換成一個radical inverse的序列即可：

$X_i:=Phi_b(i)(g_0,...,g_{n-1}) mod [0,1)^n$

Next...

這篇專欄文章介紹了幾個常見低差異序列的定義以及性質，下一篇文章（低差異序列（二）- 高效實現以及應用）會拋開理論著重談它如何高效的實現它們，以及在圖形程序中使用它們時會遇到的一些實際問題。例如如何將有限的序列用於所有的像素，如何多線程並行，如何利用這些性質並結合到一些光線傳遞的演算法中去。

Reference：http://web.maths.unsw.edu.au/~josefdick/MCQMC_Proceedings/MCQMC_Proceedings_2012_Preprints/100_Keller_tutorial.pdf