Halton序列

在統計學中,Halton序列是用於生成空間中的點的序列,如Monte Carlo模擬的數值方法,雖然這些序列是確定性的,但它們的差異性很低,也就是說,在許多方面看起來是隨機的。它們在1960年首次提出,是準隨機數列的一個例子。它們概括了一維Van der Corput序列

用於生成\(R^2\)中(0,1)x(0,1)點的Halton序列的例子

Halton數列是以質數為基的確定性方法構造的。舉個簡單的例子,讓我們把Halton序列的一個維度基於2,另一個基於3。為了生成2的序列,我們首先將區間\((0,1)\)分成兩半,然後分成四分之一、八分之一等,這就產生了

\[\frac{1}{2},\frac{1}{4},\frac{3}{4},\frac{1}{8},\frac{5}{8},\frac{3}{8},\frac{7}{8},\frac{1}{16},\frac{9}{16}...
\]

等價的,這個序列的第n個數字是用二進位制表示的數字n,倒過來,並寫在小數點之後。這對任何基數都是如此。舉個例子,要找到上述序列的第六個元素,我們要寫\(6=1*2^2+1*2^1+0*2^0=110_2\),可以倒置並放在小數點之後,得到\(0.011_2=0*2^{-1}+1*2^{-2}+1*2^{-3}=\frac{3}{8}\)。所以上面的序列與\(0.1_2,0.01_2,0.11_2,0.001_2,0.101_2,0.011_2,0.111_2,0.0001_2,0.1001_2\)相同

為了生成3的序列,我們把區間\((0,1)\)分成三份,然後是九份,二十七份,等等...這就產生了(同理表示成三進位制的數,然後進行相應操作)

\[\frac{1}{3},\frac{2}{3},\frac{1}{9},\frac{4}{9},\frac{7}{9},\frac{2}{9},\frac{5}{9},\frac{8}{9},\frac{1}{27},...
\]

當我們把它們配對起來時,我們會得到一個單位方格中的點的序列。

\[(\frac{1}{2},\frac{1}{3}),(\frac{1}{4},\frac{2}{3}),(\frac{3}{4},\frac{1}{9}),(\frac{1}{8},\frac{4}{9}),(\frac{5}{8},\frac{7}{9}),(\frac{3}{8},\frac{2}{9}),(\frac{7}{8},\frac{5}{9}),(\frac{1}{16},\frac{8}{9}),(\frac{9}{16},\frac{1}{27}).
\]

儘管標準的Halton序列在低維情況下表現的很好,但由高質數生成的序列之間存在相關問題。例如,如果我們從質數17和19開始,前16個對點數:\((\frac{1}{17},\frac{1}{19}),(\frac{2}{17},\frac{2}{19}),(\frac{3}{17},\frac{3}{19})...(\frac{16}{17},\frac{16}{19})\)具有完美的線性相關性。為了避免這種情況,通常會刪除前20個條目,或者根據選擇的指數來刪除其他預定的數量。還提出了其他幾種辦法,最突出的解決方案之一是scrambled Halton序列,它使用在構建標準序列中使用的係數的排列組合。另一個解決方案是leaped Halton,它會在標準序列中跳過點(例如,只有每409個點(也可以是其他沒有在Halton核心序列中使用的質數),才能取得顯著的改進)。

實現

虛擬碼

algorithm Halton-Sequence is
inputs: index i
base b
output: result r f1
r0 while i > 0 do
ff/b
rr+f * (i mod b)
i[i/b] return r

下面的生成器函式 generator function (Python)中給出了另一種實現方式,它可以產生以\(b\)為基數的Halton序列的後續數字。這種演算法在內部只使用整數,這使得它對四捨五入的錯誤具有很強的健壯性。

def halton(b):
"""Generator function for Halton sequence."""
n, d = 0, 1
while True:
x = d - n
if x == 1:
n = 1
d *= b
else:
y = d // b
while x <= y:
y //= b
n = (b + 1) * y - x
yield n / d

參照

補充原文

原文中陳述了很多具體的例子,而缺乏了一些Halton序列本身的說明,使用場景、以及與其他序列使用對比的差異,故在此處進行補充,更詳細的介紹可參考https://zhuanlan.zhihu.com/p/20197323

HaltonHammersley可以生成在無窮維度上分佈均勻的點集,它們都基於Van der Corput序列

Halton序列的定義很簡單:

\[X_i:=(\varPhi_{b_1}(i),...,\varPhi_{b_n}(i))
\]

既是每一個維度都是一個基於不同底數\(b_n\)的Van der Corput序列,其中\(b_1...b_n\)互為質數(例如第\(1\)到第\(n\)個質數)

Hammersley點集的定義和Halton非常相似

以下是Hammersley點集的定義

\[X_i:=(\frac{i}{N},\varPhi_{b_1}(i),...,\varPhi_{b_{n-1}}(i))
\]

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

上面左邊的圖為第1-100個Halton序列中的二維的樣本點,\((\varPhi_2(i),\varPhi_3(i))^{99}_{i=0}\),右邊則為數量為100的二維Hammersley樣本點集,\((\frac{i}{100},\varPhi_2(i))^{99}_{i=0}\)。可以看出來它們的分佈都遠比一般的偽隨機數更加均勻。Hammersley的差異性比Halton更稍低一些,但是代價是必須預先知道點的數量,並且一旦固定無法更改虛幻引擎4中對環境貼圖的Filter取樣就是用的點集大小固定為1024的Hammersley點集。Halton雖然差異性稍高,但可以不受限制的生成無窮多個點,更適合於沒有固定樣本個數的應用,例如任何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為底的序列,一開始的點會分別是

\(\frac{1}{29},\frac{2}{29},\frac{3}{29}...\)所以造成了點都集中落在了一條直線上面。解決這個問題的方法也很簡單,Scrambling。Scrambling的方法有許多種,例如最簡單的XOR Scrambling適用於以2為底數的序列。對於Halton來說,一個比較常用的方法是Faure Scrambling

\[\varPhi _b\left( i \right) =\sum_{l=0}^{M-1}{\sigma _b\left( a_l\left( i \right) \right) 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的性質。

實現

實現虛擬碼

double make_halton_sequence(int index, int base) {
double f = 1, r = 0;
while (index > 0) {
f = f / base;
r = r + f * (index % base);
index = index / base;
}
return r * randpoint_scale;
} void halton_random()
{
// 此處討論生成二維隨機數,如要產生多維,base需要是不重複的質數即可
// 三維:base 2 3 5
for (uint i = 0u; i < max_num; ++i)
{
draw_points.emplace_back(make_halton_sequence(i, 2), make_halton_sequence(i, 3));
// 三維
draw_points.emplace_back(make_halton_sequence(i, 2), make_halton_sequence(i, 3)
, make_halton_sequence(i, 5));
}
}
  • Halton序列在底數較大的時候,樣本個數會有很嚴重的correlation。所以需要採用Scrambling解決這個問題
  • RadicalInverse的實現的效率依賴於一個迴圈,將索引Index的數字左右顛倒。這一步驟可以通過一次將多個連續數字的左右顛倒連同Faure Scrambling預計算出來,存在一個查詢表裡。執行的時候直接將索引的多個數字提取出來,然後直接查表得到結果。下面給出一段以5作為底數的Halton序列實現

    詳細做法可參考:HALTON The Halton Quasi Monte Carlo (QMC) Sequence

與Hammersley序列的對比

Hammersley序列的介紹與實現可參考這篇:

Halton序列無需在生成隨機數之前,知道需要生成隨機點的個數,但是在用一些比較大的質數作為底數時,Halton序列的分佈在點的數量不那麼多的時候並不會均勻的分佈,只有當點的數量接近底數的冪的時候分佈才會逐漸均勻

效果對比

Halton序列比一般的偽隨機數更加地分佈均勻,因為此處是沒有對Halton進行優化的,即沒有Scrambling,可從另一幅圖看到,Hammersley序列比未優化的Halton序列相對來說更加地均勻,但未優化的效果也可以說是比較不錯的了

引用

翻譯自:https://en.wikipedia.org/wiki/Halton_sequence

引用部落格:

可拓展閱讀: