一个点的代价 / 读 Owen(2020)

Sobol’ 序列的第一个点永远是原点——$(0, 0, \dots, 0)$。

很多人觉得这个点不对劲。它刚好落在单位超立方体的角落里。如果你要用 Sobol’ 点做 Gaussian 分布的变换，逆正态 CDF 会把原点映射到 $-\infty$。这显然没法用。于是大家很自然地把第一个点丢掉，从第二个点开始取 $n$ 个。

这个操作叫 burn-in，在 MCMC 里是标准动作。问题是——Sobol’ 序列不是 Markov chain。

Art B. Owen 在 2020 年写了一篇短文，标题干脆利落：On dropping the first Sobol’ point。结论也干脆利落：别丢。

丢掉第一个点，你的 QMC 估计的均方根误差（RMSE）可能从 $O(n^{-3/2})$ 退化到 $O(n^{-1})$——在 $n$ 很大的时候，差别大约是 $\sqrt{n}$ 倍。

QMC 为什么比 MC 好

先回顾一下基本设定。

Monte Carlo（MC）的 RMSE 是 $O(n^{-1/2})$，对任何 $f \in L^2$ 都成立。这个收敛速度稳如老狗，但慢。

Quasi-Monte Carlo（QMC）试图做得更好。它不随机抽样，而是用精心构造的低差异序列（low-discrepancy sequence）来填满 $[0,1]^d$。Sobol’ 序列是最常用的一种。

Sobol’ 序列在 base 2 下是一个 $(t,d)$-sequence。这意味着它的前 $2^m$ 个点构成一个 $(t,m,d)$-net——对于某些被称为 elementary interval 的特殊长方体，里面的点数恰好和体积成正比。

举个二维的例子：前 16 个 Sobol’ 点会把 $[0,1]^2$ 均匀地填满所有 $1/4 \times 1/4$ 的格子，每个格子刚好一个点。$1 \times 1/16$ 的竖条，每个也刚好一个点。$1/16 \times 1$ 的横条，同理。16 个点平衡了 80 个不同形状的 elementary interval。

这就是 QMC 比 MC 好的根源——不是随机均匀，而是结构性的均匀。MC 的均匀是靠大数定律逼出来的，QMC 的均匀是设计出来的。

在这个条件下，scrambled Sobol’ 点的 RMSE 可以达到 $O(n^{-3/2}(\log n)^{(d-1)/2})$，远好于 MC 的 $O(n^{-1/2})$。

丢掉原点，丢掉结构

现在回到那个操作：丢掉第一个点。

论文的 Figure 1 画得很清楚。前 16 个 Sobol’ 二维点构成一个完美的 $(0,4,2)$-net。原点 $(0,0)$ 是其中之一（画成了同心圆）。

把它丢掉，换上第 17 个点——$(1/32, 17/32)$。

结果呢？左下角的 $1/4 \times 1/4$ 格子空了，它正上方的格子里有了两个点。$(0,4,2)$-net 的性质——没了。数字网（digital net）结构的分析基础——没了。

Owen 把这个问题用一句话说透了：

丢掉第一个点之后，新的估计等于原估计加上一个 $O(1/n)$ 的修正项。当原估计的误差是 $O(n^{-3/2})$ 时，这个修正项反而成了主导——丢掉第一个点不仅没有帮到你，还拖慢了整个收敛速度。

对 plain MC 来说，丢掉一个点无关痛痒。对 QMC 来说，丢掉一个点等于丢掉整个结构的数学保证。你以为只是在做 burn-in，实际上是把一台精密仪器砸了一个零件。

不是推测，是实测

Owen 在论文里跑了四个例子，三个合成的加一个真实的。

第一个例子：$g_0(\boldsymbol{x}) = \sum_{j=1}^d (e^{x_j} - e + 1)$，一个光滑的加性函数，$\mu=0$。用 scrambled Sobol’ 保留第一个点，RMSE 紧贴 $n^{-3/2}$ 的参考线。去掉第一个点，RMSE 紧贴 $n^{-1}$。

第二个例子：$g_1(\boldsymbol{x}) = (\sum_{j=1}^d x_j)^2$，有二阶交互。结果一模一样——保留：$n^{-3/2}$，丢掉：$n^{-1}$。

第三个例子：$g_2(\boldsymbol{x}) = \prod_{j=1}^d (e^{x_j} - e + 1)$，纯 $d$ 维交互，更难搞。丢掉第一个点的劣势没那么夸张了，但保留仍然明显更好。

第四个例子是真实世界的：一个十维的飞机机翼重量函数（wing weight function），来自实际的物理制造模型。保留第一个点的标准差明显低于丢掉后的 $O(n^{-1})$ 趋势线。

没有反例。没有一个例子中去掉第一个点会更好。

scrambling 才是正解

如果你担心原点是 $(0,0,\dots,0)$ 带来的问题——Gaussian 变换、Bayesian optimization 的 surrogate model、各种需要非零输入的场合——正确的做法不是丢掉它，是 scramble 它。

Scrambling 的做法是：对 Sobol’ 序列的每个点施加一个嵌套均匀扰动（nested uniform scramble）。扰动之后：

• 每个点 $\boldsymbol{x}_i$ 仍然服从 $\text{U}[0,1]^d$
• 整个点集以概率 1 仍然是一个 $(t,d)$-sequence
• 原点不再在边界上，而是变成单位方块内的一个随机均匀点

一口气解决了两个问题：原点不再是问题，数字网结构还在。

这也是 RQMC（Randomized QMC）的标准做法——scrambling 不仅消除了原点问题，还提供了误差估计的框架。你可以跑 $R$ 个独立的 scrambled 副本，用它们的方差来估计 RMSE。这在 plain QMC 里是做不到的，因为 QMC 本身是确定性的。

QMC 不是 MC

Owen 在 Discussion 里写了这样一段话：

MC and QMC and RQMC points all come as an $n \times d$ matrix of numbers in $[0,1]$ that we can then pipe through several functions … Despite that similarity, there are sharp differences in the properties of QMC and RQMC points that affect how we should use them.

这句话轻描淡写，但分量很重。QMC 和 MC 的输入看起来一模一样——都是 $n \times d$ 的矩阵。你把它们塞进同一个函数，MC 用 np.random.rand，QMC 用 scipy.stats.qmc.Sobol。表面上看只是换了一个随机数生成器。

但它们的"使用说明书"完全不同。

MC 不挑样本量。1000 和 1024 没区别。QMC 挑——Sobol’ 序列的最佳样本量是 2 的幂。Owen 的原文写得很直接：Using 1000 points of a Sobol’ sequence may well be less accurate than using 512。少跑 488 个点反而更准——这在 MC 的世界里是不可想象的。

MC 可以 thinning。MCMC 里每隔 $k$ 步取一个点，省存储空间。QMC 不能——把 van der Corput 序列每隔一个点取出来，你会得到所有点都在 $[0,1/2)$ 或者都在 $[1/2,1)$。一半的积分区间空了。

MC 可以 burn-in。Sobol’ 序列的第一个点是原点，不好看，丢掉——然后你就把 $(t,m,d)$-net 丢掉了。

这些 MC 视角下的"无害操作"，在 QMC 的框架里每一个都是性能杀手。

论文还点名了 MATLAB R2020a 的 sobolset 函数，里面居然提供了 Skip 和 Leap 两个参数，分别对应 burn-in 和 thinning。（好在默认值是关掉的。）Owen 对此的评价耐人寻味：It is not clear how one should use them safely. 翻译一下：我们也不知道怎么安全地用，所以你最好别用。

细思恐极

这篇论文的核心论点其实非常简单——别丢第一个点。一个七页的 proceedings paper，论证了一个单一的建议。

但它的说服力不在于复杂，而在于精准。几个合成函数加一个真实案例，log-log 图上的两条参考线，干干净净地把"丢掉第一个点会怎样"讲清楚了。没有模棱两可的"可能"、“或许”——RMSE 从 $n^{-3/2}$ 退化到 $n^{-1}$，硬数字。

我读这篇论文的时候想到一个更普遍的困境：很多工具在被"移植"到新领域时，使用者会带着旧领域的习惯。MCMC 出身的统计学家用 Sobol’ 序列，自然就想做 burn-in 和 thinning。写 sobolset 的 MATLAB 工程师，顺手就加了 Skip 和 Leap 参数。

习惯是一件很可怕的事。它让你觉得一切都在掌控中，直到有人画了一张 log-log 图告诉你，你一直在用一个 $\sqrt{n}$ 倍的减速器。

话说回来，我自己写 QMC 代码的时候，干的蠢事一定比丢掉第一个点更离谱。只是没人帮我画 log-log 图而已。