低差异序列 (low-discrepancy sequences)之Halton序列均匀产生多维随机数的介绍与实现

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

	f⬅1
	r⬅0

	while i > 0 do
		f⬅f/b
		r⬅r+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

参照

低差异序列的构造 Constructions of low-discrepancy sequences

补充原文

原文中陈述了很多具体的例子，而缺乏了一些Halton序列本身的说明，使用场景、以及与其他序列使用对比的差异，故在此处进行补充，更详细的介绍可参考https://zhuanlan.zhihu.com/p/20197323

Halton和Hammersley可以生成在无穷维度上分布均匀的点集，它们都基于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