一个点的代价 / 读 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 图而已。