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序列生成器。

## 相关

- 低差异序列
- 拟蒙特卡罗方法