5.3 分层与低差异序列
5.2 量完了那条 $N^{-1/2}$ 的定律，本节换一个问法：能不能不堆样本，而是把样本放得更聪明。$N^{-1/2}$ 里的指数是随机采样自带的代价，真正被被积函数与 pdf 的比值决定的是那个 $\sigma^2$——5.2 结尾说过，分层与低差异序列改的不是指数，是分母。实时 GI 那篇靠时间上的重投影把噪声摊到几十帧里（见 实时全局光照 ），离线光追没有前一帧可借，只能在单帧内把点的分布做对。
分层采样：把立方体切成格子
分层采样（stratified sampling）的动机很直接：如果 16 个点恰好挤在单位正方形的一角，这一次估计就废了。那就先把 $[0,1]^2$ 切成 $4 \times 4$ 的 16 个格子，每个格子强制放一个样本，再让样本在格子内部随机抖动。stratify_grid 就是两层循环：外层走格子编号，内层算 ((g << 16) + v) // GRID——$g$ 是格子号，$v$ 是格子内部的 Q16.16 抖动，除以 4 把「格子 + 抖动」压回 $[0,1)$。四个格子的值域恰好是 $[0, 16384)$、$[16384, 32768)$、$[32768, 49152)$、$[49152, 65536)$，一格一个，谁也不越界。
分层换来的是方差下降，不是偏差消除：格子内部仍然是均匀抖动，估计量依旧无偏，只是 $\sigma^2$ 被「每格必有一个点」这条约束削掉一截。代价是一条硬约束——样本数必须能被层数整除。脚本里 16 个格子、5 档 $N$ 从 16 起都是 16 的倍数，n // (GRID * GRID) 恰好把每层填满；如果写成 $N = 20$，20 // 16 = 1，实际只用上 16 个样本，剩下 4 个被悄悄丢掉，误差曲线会比 5.2 那张诊断表里的 $-0.5$ 提前跑偏。
四种采样器的误差对比
被积函数取 $d$ 维单位立方体里 $\sum x_i^2 \le 1$ 的占比，闭式解是单位 $d$ 球体积除以 $2^d$：$d = 2$ 时就是 $\pi/4 = 0.785398$。它既不像 $x^2$ 那样平滑（指示函数在球面上有跳变），又有真值可查，正好当四种采样器的共同靶子。每种采样器跑 16 个实现（换子流、换抖动量），报 16 个相对误差的 RMS：
#!/usr/bin/env python3
# file: lantern/sim/cg_trace_05_3_strata.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_trace_05_3_strata.py
"""5.3 分层与低差异序列：同一场景上四种采样器的误差对比，以及维度一升就失效。
场景是 d 维单位立方体里「Σ x_i² ≤ 1」的占比，闭式解 = 单位 d 球体积 / 2^d，d = 2 时
真值 π/4 = 0.785398。四种采样器：朴素随机、4×4 分层（stratify_grid）、Halton（整数版
radical inverse）、Sobol（基 2 方向数，整数版）。每种跑 16 个实现（换子流 / 换抖动 /
旋转 / 数字移位），报 RMS 相对误差并拟合 log-log 斜率。Sobol 用 (0,2q,2)-net 的对齐网格
性质自检：前 4^q 个点落在 2^q×2^q 网格里每格恰好一个，q = 1…6。方向数只覆盖 2 维，高维
退化改用 Halton 演示，固定 N = 1024 逐维升到 8 维。
产物：config/cg_trace_05_3_strata.json
"""
import json
import math
import os
import sys
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import cglib # noqa: E402
import cgtrace # noqa: E402
ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
ONE = cglib.ONE
SEED0 = 0x00C0FFEE_12345678
STREAM_STRIDE = 1 << 40
N_LIST = (16, 64, 256, 1024, 4096)
REPS = 16
GRID = 4 # 4×4 分层
DIMS = (2, 4, 6, 8)
DIM_N = 1024
PRIMES = (2, 3, 5, 7, 11, 13, 17, 19)
EXACT = {d: math.pi ** (d / 2) / math.gamma(d / 2 + 1) / (1 << d) for d in DIMS}
# 第 1 维写 None 走 van der Corput 特例（方向数 2^(32-i)）；第 2 维是 (a, m) = (0, (1,3))。
SOBOL = (None, (0, (1, 3)))
def stream(index):
"""第 index 条独立子流的起点（index 从 1 开始）。"""
return (SEED0 + index * STREAM_STRIDE) & cgtrace.RNG_MASK
def shifts(base, d):
"""从一条子流取 d 个 16 位移位量，供旋转 / 数字移位使用。"""
out, s = [], base
for _ in range(d):
s, v = cgtrace.rng_q16(s)
out.append(v)
return tuple(out)
def radical_inverse(i, base):
"""整数版 radical inverse：i 的 base 进制数字反序。Q32 中间量，输出 Q16.16。"""
f, r = 1 << 32, 0
while i:
f //= base
r += f * (i % base)
i //= base
return r >> 16
def halton_pts(n, d, sh):
"""前 n 个 d 维 Halton 点 + Cranley–Patterson 旋转（逐维加常数再模 1）。"""
return [tuple((radical_inverse(i, PRIMES[k]) + sh[k]) & (ONE - 1) for k in range(d))
for i in range(1, n + 1)]
def sobol_dirs(spec, bits=32):
"""方向数 V[1..bits]：spec 为 None 时是 van der Corput，否则按 (a, m) 递推。"""
if spec is None:
return [0] + [1 << (bits - i) for i in range(1, bits + 1)]
a, m = spec
s, v = len(m), [0] * (bits + 1)
for k in range(1, s + 1):
v[k] = m[k - 1] << (bits - k)
for k in range(s + 1, bits + 1):
v[k] = v[k - s] ^ (v[k - s] >> s)
for i in range(1, s):
v[k] ^= ((a >> (s - 1 - i)) & 1) * v[k - i]
return v
def sobol_pts(n, d, sh):
"""前 n 个 d 维 Sobol 点（0 基，整数版）+ 数字移位：每步翻转最低零位对应的方向数。"""
vs = [sobol_dirs(s) for s in SOBOL[:d]]
xs, out = [0] * d, [tuple(sh)]
for i in range(1, n):
c, t = 1, i - 1
while t & 1:
t >>= 1
c += 1
for k in range(d):
xs[k] ^= vs[k][c]
out.append(tuple((x >> 16) ^ sh[k] for k, x in enumerate(xs)))
return out
def naive(state, n, d):
"""朴素随机：n 个独立均匀 d 维点。"""
pts, s = [], state
for _ in range(n):
row = []
for _ in range(d):
s, v = cgtrace.rng_q16(s)
row.append(v)
pts.append(tuple(row))
return pts
def stratify_grid(state, n):
"""4×4 分层：16 个格子各放 n/16 个抖动样本，抖动落在格子内部。"""
pts, s = [], state
for gy in range(GRID):
for gx in range(GRID):
for _ in range(n // (GRID * GRID)):
row = []
for g in (gx, gy):
s, v = cgtrace.rng_q16(s)
row.append(((g << 16) + v) // GRID)
pts.append(tuple(row))
return pts
def score(pts, t):
"""相对误差 |Σ x_i² ≤ 1 的占比 − t| / t；Q16.16 判据是 Σ x_i² ≤ ONE。"""
return abs(sum(1 for p in pts if sum(cglib.fmul(v, v) for v in p) <= ONE) / len(pts) - t) / t
def rms(vals):
"""均方根。"""
return math.sqrt(sum(v * v for v in vals) / len(vals))
def fit(rows, key):
"""最小二乘斜率：log10(误差) ~ log10(N)。"""
xs = [math.log10(r["n"]) for r in rows]
ys = [math.log10(r[key]) for r in rows]
mx, my = sum(xs) / len(xs), sum(ys) / len(ys)
return sum((x - mx) * (y - my) for x, y in zip(xs, ys)) / sum((x - mx) ** 2 for x in xs)
def part1():
"""d = 2 场景：四种采样器 × 五档 N × 16 个实现的 RMS 相对误差。"""
t2, rows = EXACT[2], []
for n in N_LIST:
e = {k: [] for k in ("random", "strat", "halton", "sobol")}
for rep in range(REPS):
base = stream(rep + 1)
sh = shifts(base, 2)
e["random"].append(score(naive(base, n, 2), t2))
e["strat"].append(score(stratify_grid(base, n), t2))
e["halton"].append(score(halton_pts(n, 2, sh), t2))
e["sobol"].append(score(sobol_pts(n, 2, sh), t2))
rows.append({"n": n, **{k: round(rms(v), 6) for k, v in e.items()}})
return rows
def part2():
"""固定 N = 1024、逐维升：低差异序列相对随机的误差优势怎么衰减。"""
rows = []
for d in DIMS:
td, e = EXACT[d], {"random": [], "halton": [], "sobol": []}
for rep in range(REPS):
base = stream(200 + rep)
sh = shifts(base, d)
e["random"].append(score(naive(base, DIM_N, d), td))
e["halton"].append(score(halton_pts(DIM_N, d, sh), td))
if d <= len(SOBOL):
e["sobol"].append(score(sobol_pts(DIM_N, d, sh), td))
row = {"d": d, "truth": round(td, 6), "sobol": None}
for k in ("random", "halton", "sobol"):
if e[k]:
row[k] = round(rms(e[k]), 6)
if k != "random":
row[k + "Gain"] = round(row[k] / row["random"], 4)
rows.append(row)
return rows
def net_check(pts, q):
"""(0,2q,2)-net 的必要条件：前 4^q 个点在 2^q×2^q 网格里每格恰好一个。"""
cells = {}
for x, y in pts:
cells[(x >> (16 - q), y >> (16 - q))] = cells.get((x >> (16 - q), y >> (16 - q)), 0) + 1
return len(cells) == 4 ** q and max(cells.values()) == 1
def main():
os.makedirs(os.path.join(ROOT, "out"), exist_ok=True)
rows2, rowsd = part1(), part2()
fits = {k: round(fit(rows2, k), 4) for k in ("random", "strat", "halton", "sobol")}
nets = {str(q): net_check(sobol_pts(4 ** q, 2, (0, 0)), q) for q in range(1, 7)}
report = {"seed": str(SEED0), "streamStride": STREAM_STRIDE, "reps": REPS, "grid": GRID,
"nList": list(N_LIST), "truth2": EXACT[2], "fits": fits, "nets": nets,
"dims": list(DIMS), "dimN": DIM_N, "rows2": rows2, "rowsDim": rowsd,
"exact": {str(d): round(EXACT[d], 6) for d in DIMS},
"sobol": [list(s) if s else "vdC" for s in SOBOL]}
with open(os.path.join(ROOT, "config/cg_trace_05_3_strata.json"), "w", encoding="utf-8") as fh:
json.dump(report, fh, ensure_ascii=False, indent=2)
print(f"== 1. d = 2 场景（真值 π/4 = {EXACT[2]:.6f}），16 个实现的 RMS 相对误差")
for r in rows2:
print(f" N = {r['n']:>4}：随机 {r['random']:.6f}｜分层 {r['strat']:.6f}"
f"｜Halton {r['halton']:.6f}｜Sobol {r['sobol']:.6f}")
print("== 2. log-log 斜率（-0.5 即朴素 MC 定律）")
print(" " + "｜".join(f"{k} {fits[k]:+.4f}" for k in ("random", "strat", "halton", "sobol")))
print("== 3. Sobol 自检：(0,2q,2)-net 对齐网格每格恰好一个点")
print(" " + "｜".join(f"q = {q}: {nets[q]}" for q in map(str, range(1, 7))))
print(f"== 4. 固定 N = {DIM_N}，维度上升时低差异序列相对随机的误差比")
for r in rowsd:
s = "" if r["sobol"] is None else f"｜Sobol {r['sobol']:.6f}（×{r['sobolGain']:.3f}）"
print(f" d = {r['d']}：真值 {r['truth']:.6f}｜随机 {r['random']:.6f}"
f"｜Halton {r['halton']:.6f}（×{r['haltonGain']:.3f}）{s}")
return 0
if __name__ == "__main__":
sys.exit(main())
| $N$ | 朴素随机 | 4×4 分层 | Halton | Sobol |
|---|---|---|---|---|
| 16 | 0.119703 | 0.075326 | 0.087940 | 0.061446 |
| 64 | 0.049543 | 0.025500 | 0.026232 | 0.032969 |
| 256 | 0.021641 | 0.016915 | 0.007283 | 0.006916 |
| 1024 | 0.016068 | 0.009170 | 0.004357 | 0.002602 |
| 4096 | 0.006658 | 0.003355 | 0.001570 | 0.000862 |
低 $N$ 时四家几乎打平——$N = 16$ 那一行，随机 0.119703 与 Sobol 0.061446 只差一倍；$N$ 一大就分层。$N = 4096$ 时 Halton 的误差 0.001570 是随机 0.006658 的 $1/4.24$，把它换算成随机样本，要 7.4 万个才能追平 4096 个 Halton 点（$4.24^2 \approx 18$ 倍）。同样的计算预算，点放得好就等于把样本数乘了 18。脚本同时把每条曲线拟合成一个斜率：
== 2. log-log 斜率（-0.5 即朴素 MC 定律）
random -0.4980｜strat -0.5227｜halton -0.7103｜sobol -0.7987
== 3. Sobol 自检：(0,2q,2)-net 对齐网格每格恰好一个点
q = 1: True｜q = 2: True｜q = 3: True｜q = 4: True｜q = 5: True｜q = 6: True
斜率比 $-0.5$ 更陡，意思是误差不只被压小，压小的速度还更快。分层只陡了 0.02：它固定的是「一层一个点」这一件事，层内的抖动仍然是随机的，$4 \times 4$ 把 $\sigma^2$ 除掉一个常数因子，指数没动。Halton 与 Sobol 是低差异序列（low-discrepancy sequence, LDS）：它们不抖动，而是用一个确定的整数递推直接吐出「尽量不扎堆」的点列，误差里混进一项按 $(\log N)^d / N$ 走的准蒙特卡洛项，$N$ 一大就压过 $N^{-1/2}$ 项，斜率于是明显更陡。
Sobol：方向数与一次能挂掉的自检
Halton 的实现是整数版 radical inverse：把下标 $i$ 写成 base 进制、数字反序、小数点放最前，全程用 Q32 整数移位与整除，不碰浮点。Sobol 走另一条路——先准备一组方向数 $V_1, V_2, \dots$，第 $i$ 个点由 Gray 码决定这一步异或哪几个方向数：$X_i = X_{i-1} \oplus V_c$，其中 $c = 1 + (i-1)$ 的尾 1 个数。方向数本身由一条线性递推生成，$V_k = V_{k-s} \oplus (V_{k-s} \gg s)$，初值取 $V_1 = 2^{31}$、$V_2 = 3 \cdot 2^{30}$。
序列好不好不能靠肉眼看散点图。一条 base 2 的 Sobol 满足 $(0, 2q, 2)$-net：前 $4^q$ 个点落在 $2^q \times 2^q$ 的对齐网格里每格恰好一个。脚本对 $q = 1 \dots 6$ 逐个验算，六档全 True。这个自检很能咬人，三个都是眼睛看不出来的实现细节：把第 2 维初值从 $(1, 3)$ 改成 $(1, 1)$，$X_0$ 到 $X_3$ 全部落在对角线上，$q = 1$ 的四格只命中两格；把第 1 维也交给那条递推（$V_1 = 2^{31}$ 起、步长 $s = 1$），二维投影从 $q = 1$ 起就全 False；把点的生成改成从 $X_1$ 开始（跳过原点），$q = 1$ 还能蒙混过关，$q = 2$ 起立刻现形。
维度一升，优势就掉
低差异序列的误差上界里坐着一个 $d$：准蒙特卡洛是 $O((\log N)^d / N)$，维度在指数上。同一场景升到 4、6、8 维，固定 $N = 1024$：
| 维度 $d$ | 真值 | 朴素随机 | Halton | 随机 / Halton |
|---|---|---|---|---|
| 2 | 0.785398 | 0.013722 | 0.003970 | 3.46× |
| 4 | 0.308425 | 0.045469 | 0.022492 | 2.02× |
| 6 | 0.080746 | 0.107205 | 0.059948 | 1.79× |
| 8 | 0.015854 | 0.280642 | 0.194379 | 1.44× |
两件事同时发生。一是随机采样自己变差了：$d = 2$ 的相对误差 0.013722，$d = 8$ 涨到 0.280642。这不是采样器退化，是被积函数变了——命中率 $p$ 从 $\pi/4 = 0.785$ 掉到 0.0159（49 倍），而二项分布的相对标准差是 $\sqrt{(1-p)/(pN)}$，单看 $1/\sqrt{p}$ 这一项就预测误差要涨 15 倍；表里实测 20 倍，量级对上了，多出来的那截是 16 次实现的统计噪声。二是Halton 的优势收窄：相对随机的误差比从 3.46 倍掉到 1.44 倍，换算成「等效随机样本数」就是从 12 倍掉到 2 倍。8 维里那点红利，还不够抵消一次实现失误——这正是把低维调参经验外推到路径追踪（维数动辄几十上百）时会踩的坑。
P 层：BigInt 里的方向数与 Gray 码
跨语言断言在这里有个坑：基 2 的方向数会用到第 31 位（$V_1 = 2^{31}$），而 JavaScript 的位运算先把操作数压成 32 位有符号整数，$2^{31}$ 右移会变成负数，异或结果全错。所以 TypeScript 侧把方向数、点坐标、异或、右移一律走 BigInt，只在最后取 Q16.16 时转回 Number：
// snippet: 取自 lantern/core/src/lantern/cg_trace_05_3_strata.ts 的方向数递推与 Gray 码生成（逐字摘录，完整文件另含 radicalInverse、haltonPts、stratifyGrid 与 netCheck）
/** 方向数 V[1..bits]：spec 为 null 时是 van der Corput，否则按 (a, m) 递推。 */
export function sobolDirs(spec: null | [bigint, bigint[]], bits = 32): bigint[] {
if (spec === null) {
const v = [0n];
for (let i = 1; i <= bits; i += 1) v.push(1n << BigInt(bits - i));
return v;
}
const [a, m] = spec;
const s = m.length;
const v: bigint[] = new Array(bits + 1).fill(0n);
for (let k = 1; k <= s; k += 1) v[k] = m[k - 1] << BigInt(bits - k);
for (let k = s + 1; k <= bits; k += 1) {
v[k] = v[k - s] ^ (v[k - s] >> BigInt(s));
for (let i = 1; i < s; i += 1) v[k] ^= ((a >> BigInt(s - 1 - i)) & 1n) * v[k - i];
}
return v;
}
/** 前 n 个 d 维 Sobol 点（0 基，整数版）+ 数字移位：每步翻转最低零位对应的方向数。 */
export function sobolPts(n: number, d: number, sh: bigint[]): bigint[][] {
const vs = SOBOL.slice(0, d).map((s) => sobolDirs(s));
const xs: bigint[] = new Array(d).fill(0n);
const out: bigint[][] = [sh.slice()];
for (let i = 1; i < n; i += 1) {
let c = 1;
let t = i - 1;
while (t & 1) {
t >>= 1;
c += 1;
}
for (let k = 0; k < d; k += 1) xs[k] ^= vs[k][c];
out.push(xs.map((x, k) => (x >> 16n) ^ sh[k]));
}
return out;
}
八条断言全绿。三条硬口径：四种采样器 × 5 档的 RMS 相对误差复现到 $10^{-6}$ 内（子流起点、抖动、点位全部逐点重算）、log-log 斜率复现到 $10^{-4}$ 内（随机 $-0.4980$ 落在 $-0.5 \pm 0.03$，Halton 与 Sobol 必须分别比 $-0.6$、$-0.7$ 更陡）、Sobol 的 $q = 1 \dots 6$ 自检在两侧都为真。另有三条反向判据：四种采样器在每一档都必须赢过随机、$d = 8$ 的 Halton 增益必须比 $d = 2$ 大一倍以上、单位 $d$ 球占比的闭式解必须与 Python 的 math.gamma 口径对上。到这里，$\sigma^2$ 已经被压下去一截；下一节把这一节的估计量装进真正的渲染方程。
阅读导航：上一节：5.2 方差与收敛 · 下一节：6.1 渲染方程 。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。