Sobol序列

256個點：來自前256個點的2,3 Sobol序列（上）與偽隨機數源（下）的比較。Sobol序列更均勻覆蓋空間。（紅色=1,..,10，藍色=11,..,100，綠色=101,..,256）

Sobol序列（亦稱LPτ序列或基2的(t, s)序列）是一類准隨機低差異序列。它們由俄國數學家Ilya M. Sobol'（Илья Меерович Соболь）於1967年首次提出。

這些序列使用以二為基底將單位區間依次劃分為更細的均勻子區間，並在每一維度中對坐標重排。

## 在s維單位超立方體中的良好分布

設1=Is = 0,1s為s維單位超立方體，且f是定義在Is上的可積實函數。Sobol'最初的動機是構造序列xn（位於Is中），使得$$\lim_{n\to\infty} \frac{1}{n} \sum_{i=1}^n f(x_i) = \int_{I^s} f$$且收斂儘可能快。

使和收斂到積分，較為明顯的是點xn應儘量填滿Is以最小化空洞。另一個好的性質是xn在低維面的投影也應儘量少留空洞。因此簡單的均勻填充Is並不合適，因為在低維投影中許多點會重合，因而對積分估計無用。

這些良好分布稱為基$b$的(t,m,s)-網與(t,s)-序列。先定義基$b$的基本s-區間為形式$$\prod_{j=1}^s \left[ \frac{a_j}{b^{d_j}}, \frac{a_j+1}{b^{d_j}} \right]$$其中$a_j$與$d_j$為非負整數，且對所有$j\in${1,…,s}有$a_j < b^{d_j}$。

給定整數$0\leq t\leq m$，基$b$的(t,m,s)-網是位於$I^s$的$b^m$個點的序列$x_n$，使得對任意基$b$的基本區間 `''P''`，當其超體積為$\lambda(P)=b^{t-m}$時，有 $\operatorname{Card} P \cap \{x_1, ..., x_{b^m}\} = b^t$。

給定非負整數$t$，基$b$的(t,s)-序列是一個無限序列$x_n$，滿足對任意整數$k \geq 0, m \geq t$序列$\{x_{kb^m}, ..., x_{(k+1)b^m-1}\}$是基$b$的一個(t,m,s)-網。

Sobol'在他的文章中描述了Πτ-網格與LPτ-序列，它們分別是基2的(t,m,s)-網與(t,s)-序列。基$b$的(t,m,s)-網與(t,s)-序列（也稱為Niederreiter序列）這一術語由Harald Niederreiter在1988年提出。「Sobol序列」這一術語在晚期的英語文獻中被用來與Halton序列、Faure及其它低差異序列相區分。

## 快速算法

Antonov與Saleev提出了一種更高效的格雷碼實現。

在生成Sobol數時，使用格雷碼$G(n)=n \oplus \lfloor n/2 \rfloor$代替$n$來構造第n個點是很有幫助的。

設已生成直到$n-1$的所有Sobol序列值，並在記憶體中保存了所有維度所需的值$x_{n-1,j}$。由于格雷碼$G(n)$與前一項$G(n-1)$只在單一位上不同（該位為n−1的最右側零位），因此對每一維度只需做一次異或（XOR）運算即可將所有$x_{n-1}$推導為$x_n$，即

$$x_{n,i} = x_{n-1,i} \oplus v_{k,i}.$$

## 額外的均勻性性質

Sobol'引入了稱為性質A與A'的附加均勻性條件。

定義：若對任意二進段（而非任意子集）長度為$2^d$的d維序列段，單位超立方體沿每條邊各自二等分所得到的$2^d$個超立方體中恰有且只有一個抽樣點，則稱該低差異序列滿足性質A。

定義：若對任意二進段（而非任意子集）長度為$4^d$的d維序列段，單位超立方體沿每條邊各自四等分所得到的$4^d$個超立方體中恰有且只有一個抽樣點，則稱該低差異序列滿足性質A。

存在保證性質A與A'的數學條件。

d維Sobol序列若且唯若滿足$$\det(\mathbf{V}_d) \equiv 1 \pmod{2},$$ 其中 Vd 為下述 d × d 的二值矩陣：$$\mathbf{V}_d := \begin{pmatrix} {v_{1,1,1}}&{v_{2,1,1}}&{\dots}&{v_{d,1,1}}\\ {v_{1,2,1}}&{v_{2,2,1}}&{\dots}&{v_{d,2,1}}\\ {\vdots}&{\vdots}&{\ddots}&{\vdots}\\ {v_{1,d,1}}&{v_{2,d,1}}&{\dots}&{v_{d,d,1}} \end{pmatrix},$$其中 vk,j,m 表示方向數 vk,j = (0.vk,j,1vk,j,2...)2在二進小數點之後的第 m 位數字。d維Sobol序列若且唯若滿足：$$\det(\mathbf{U}_d) \equiv 1 \pmod{2},$$其中 Ud 為下述 2d × 2d 的二值矩陣：$$\mathbf{U}_d := \begin{pmatrix} {v_{1,1,1}}&{v_{1,1,2}}&{v_{2,1,1}}&{v_{2,1,2}}&{\dots}&{v_{d,1,1}}&{v_{d,1,2}}\\ {v_{1,2,1}}&{v_{1,2,2}}&{v_{2,2,1}}&{v_{2,2,2}}&{\dots}&{v_{d,2,1}}&{v_{d,2,2}}\\ {\vdots}&{\vdots}&{\vdots}&{\vdots}&{\ddots}&{\vdots}&{\vdots}\\ {v_{1,2d,1}}&{v_{1,2d,2}}&{v_{2,2d,1}}&{v_{2,2d,2}}&{\dots}&{v_{d,2d,1}}&{v_{d,2d,2}} \end{pmatrix},$$其中 vk,j,m 表示方向數 vk,j = (0.vk,j,1vk,j,2...)2在二進小數點之後的第 m 位數字。

對性質A與A'的檢驗是相互獨立的。因此可以構造同時滿足A和A'的Sobol序列，或僅滿足其中之一的序列。

## Sobol數的初始化

要構造Sobol序列，需要選擇一組方向數$v_{i,j}$。初始化方向數的選擇存在一定自由度。因此針對不同維數可能得到不同的Sobol序列實現。若初始化數選擇不佳，會顯著降低Sobol序列在計算中的效率。

一種最容易的初始化選擇是令第l位最左側位為1，其餘位為0，即對所有$k$和$j$取$m_{k,j}=1$。此初始化通常稱為「單位初始化」。然而，該初始化即使在低維度下也不能通過性質A與A'的檢驗，因此屬於較差的初始化。

## 實現與可獲得性

針對不同維數的良好初始化數已由若干作者提供。例如，Sobol'提供了最高到51維的初始化數。同一組初始化數也被Bratley和Fox使用。

高維度的初始化數可在Joe和Kuo處獲得。Peter Jäckel在其著作《金融領域的蒙特卡羅方法》中提供了到32維的初始化數。

其它實現作為C、Fortran 77或Fortran 90例程出現在《Numerical Recipes》軟體集。基於Joe和Kuo的初始化數，有一個可支持最多1111維的自由/開源C實現（見NLopt庫），以及一個可支持最多21201維的Python實現和Julia實現。另一套可支持最多1111維的自由/開源實現可用於C++、Fortran 90、MATLAB與Python。

商業的Sobol序列生成器可在例如NAG Library中獲得。BRODA Ltd.提供具有附加均勻性性質A與A'的Sobol與擾動Sobol序列生成器，最大維度可達131072；這些生成器由BRODA與Ilya M. Sobol協同開發。MATLAB的統計工具箱亦包含最多到1111維的Sobol序列生成器。

## 相關

- 低差異序列
- 擬蒙特卡羅方法