12.3 采样的停止准则
站内 /graphics-gpu-profiling-debugging/ 那篇做的是 GPU 剖析与调试：量一帧里每个 pass 花了多少时间，然后决定优化谁。本章量的是离线渲染的误差，而它的用途正好相反——不是事后分析，而是当场决定还要不要继续采样。前两节把噪声和误差都量成了数，这一节拿这些数去做决定。
三条准则，逐像素各自判停
均匀采样把每个像素都跑到同一个 N。可像素之间的方差差着几个数量级：平坦地面上的像素几十个样本就稳了，阴影边缘和面光源掠射区可能要上千。自适应采样的想法很直白——让每个像素自己决定什么时候停。
本节实现三条准则，逐像素判停。V（逐像素方差） 比的是流内 s/√N 与固定阈值 64 Q16.16（≈ 1/4 个 8 位色阶）；R（相对误差） 比的是 s/√N ÷ 当前均值 与千分之二十五（2.5%）；B（预算上限） 不设阈值，N 到 256 无条件停，是兜底，也是高方差区域的实际行为。
V 与 R 各跑两个最小样本数版本（N ≥ 4 与 N ≥ 32）。加最小样本数是必须的：方差估计本身是随机变量，上一节已经看到 N = 1 时它按定义为 0，而四个样本估出来的方差几乎全是噪声。这一档对照就是为了把「最小样本数该取多少」量出来。
同一条流，回放判停
三条准则共用同一条流：每像素跑满 256 个样本，在 4, 8, …, 256 共 64 个检查点记下整数前缀和 Σx、Σx²，再回放判断各自停在哪一档。共用一条流的好处是——三条准则看到的是同一批随机数，差异只来自判停规则本身，不掺入采样路径的不同。
#!/usr/bin/env python3
# file: lantern/sim/cg_trace_12_3_stopping.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_trace_12_3_stopping.py
"""12.3 采样的停止准则：什么时候可以不再采样。
三条准则逐像素判停，共用**同一条流**——每像素跑满 256 个样本，在 4, 8, …, 256 共
64 个检查点记下整数前缀和，再回放判断各自停在哪一档。三条准则看到的是同一批随机数，
差异只来自判停规则本身：
* **V（逐像素方差）**：流内 `s/√N` ≤ `varianceThresholdQ16`（64 Q16.16 ≈ 1/4 个
8 位色阶）。
* **R（相对误差）**：`s/√N ÷ 当前均值` ≤ `relErrorThresholdPermille`（千分之二十五）。
* **B（预算上限）**：N 到 256 无条件停，兜底，也是高方差区域的实际行为。
V 与 R 各跑两个最小样本数版本（N ≥ 4 与 N ≥ 32）：**方差估计本身是随机变量**，四样本
估出来的方差几乎全是噪声，这一档对照就是为了把「最小样本数」量出来。
四个必须量出来的东西：停在哪（直方图 + 总样本 + 实际 RMS）、自估误差诚不诚实（自估
`s/√N` 对实际误差的比值、超 3σ 计数）、两种失效（最小档就停 / 撞预算上限）、以及与
**同预算均匀采样**的对照（自适应到底是省了还是亏了）。
误差的参照是 12.2 的参考解 A（16 批 × 64 = 1024 样本/像素，同一 seed/stride），脚本会
断言它的 sha256 与 12.2 一致——两个脚本量的是同一个「真值」。
产物：out/cg_trace_12_3_stopping.json、out/cg_trace_12_3_stopping.csv、
out/cg_trace_12_3_{variance,relative,uniform,nmap}.png
"""
import json
import math
import os
import sys
import time
import numpy as np
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import cglib # noqa: E402
import cgtrace # noqa: E402
import cg_trace_12_common as C # noqa: E402
ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
OUT = os.path.join(ROOT, "out")
REF_SHA = "0ac7328aaeb5c669666292b8d3e6750d16163a95c572c1791f99e913646507d7"
def main():
os.makedirs(OUT, exist_ok=True)
spec = cglib.load_spec()
cfg = C.load_area_cfg()
st, rf = cfg["stopping"], cfg["reference"]
budget, step = st["budget"], st["checkEvery"]
var_thr = float(st["varianceThresholdQ16"])
rel_thr = float(st["relErrorThresholdPermille"])
w, h = st["width"], st["height"]
tris, _ = cgtrace.build_scene_tris(spec)
bvh = cgtrace.build_bvh(tris)
light = C.make_light(spec, cfg)
pixels = C.prepare_pixels(spec, tris, bvh, w, h, light)
idx = [i for i, p in enumerate(pixels) if p is not None]
m, ar = len(idx), np.arange(len(idx))
# ---- 参考解 A：与 12.2 同 seed/stride，断言 sha256 一致 ----
t0 = time.time()
perA, totA = C.batch_sums(pixels, idx, light, tris, bvh, rf["batches"],
rf["batchSamples"], rf["seed"], rf["batchStride"],
rf["streamStride"])
ref_ms = int((time.time() - t0) * 1000)
ref = (totA // (rf["batches"] * rf["batchSamples"])).astype(np.float64)
bmA = perA.astype(np.float64) / rf["batchSamples"]
sem_ref = bmA.std(axis=0, ddof=1) / math.sqrt(rf["batches"])
sem_ref_rms = float(np.sqrt((sem_ref ** 2).mean()))
ref_sha = C.digest_light_only(pixels, idx, w, h, tuple(spec["backgroundQ8"]), ref)
assert ref_sha == REF_SHA, f"参考解与 12.2 不一致：{ref_sha}"
# ---- 一条流跑满预算，取 64 个检查点 ----
t0 = time.time()
stats = {"rays": 0, "occluded": 0}
lv, S, Q = C.ladder_sums(pixels, idx, light, tris, bvh, budget, step,
st["seed"], st["streamStride"], stats)
lad_ms = int((time.time() - t0) * 1000)
n3 = np.array(lv, dtype=np.float64)[:, None, None]
Sf, Qf = S.astype(np.float64), Q.astype(np.float64)
mean = Sf / n3
sem = np.sqrt(np.maximum((Qf * n3 - Sf * Sf) / (n3 * (n3 - 1.0)), 0.0) / n3)
sem_px, rel_px = sem.max(axis=2), (sem / np.maximum(mean, 1e-9)).max(axis=2)
# 准则 B 的 mask 只在最后一档为真——「N 到预算就停」等价于「停在最后一档」。
L = len(lv)
base = {"V": sem_px <= var_thr, "R": rel_px * 1000.0 <= rel_thr,
"B": np.zeros_like(sem_px, dtype=bool)}
base["B"][L - 1] = True
i32 = lv.index(32)
rows, res = [], {}
for tag, kind, mi in [("V4", "V", 0), ("V32", "V", i32), ("R4", "R", 0),
("R32", "R", i32), ("B", "B", 0)]:
mask = base[kind].copy()
if mi > 0:
mask[:mi] = False
ns, i_eff, found, est = C.apply_criterion(mask, lv, budget, S)
err = est - ref
sem_at = sem[i_eff, ar] # (m, 3) 判停时自估的标准误
rmse = float(np.sqrt((err ** 2).mean()))
sem_rms = float(np.sqrt((sem_at ** 2).mean()))
at_min = ns == lv[mi]
tot = int(ns.sum())
nu = max(lv[0], min(budget, (tot // m // step) * step))
rmse_u = float(np.sqrt(((S[lv.index(nu)].astype(np.float64) / nu - ref) ** 2).mean()))
rows.append({
"tag": tag, "totalSamples": tot, "samplesPerPixelMean": round(float(ns.mean()), 3),
"stopAtMin": int(at_min.sum()),
"zeroVarAtMin": int((at_min & (sem_at.max(axis=1) == 0)).sum()),
"cappedAtBudget": int((~found).sum()),
"rmseQ16": round(rmse, 4), "maxErrQ16": round(float(np.abs(err).max()), 4),
"selfSemQ16": round(sem_rms, 4), "honestyRatio": round(rmse / sem_rms, 4),
"over3Sigma": int((np.abs(err).max(axis=1) >
3.0 * np.maximum(sem_at.max(axis=1), 1e-9)).sum()),
"uniformN": nu, "uniformRmseQ16": round(rmse_u, 4),
"uniformRatio": round(rmse_u / max(rmse, 1e-9), 4),
"hist": C.histogram(ns, lv, budget),
})
res[tag] = {"ns": ns, "i": i_eff, "found": found, "est": est, "rmse": rmse}
# ---- 失效模式细看（以准则 V4 为例）----
r4 = res["V4"]
err_v = r4["est"] - ref
sem_v = sem[r4["i"], ar]
early, capped = r4["ns"] == lv[0], ~r4["found"]
mid = ~early & ~capped
e_rmse = float(np.sqrt((err_v[early] ** 2).mean())) if early.any() else 0.0
c_rmse = float(np.sqrt((err_v[capped] ** 2).mean())) if capped.any() else 0.0
m_rmse = float(np.sqrt((err_v[mid] ** 2).mean())) if mid.any() else 0.0
e_over = int((np.abs(err_v[early]).max(axis=1) >
3.0 * np.maximum(sem_v[early].max(axis=1), 1e-9)).sum())
tot_mse = m * r4["rmse"] ** 2
share = {"early": round(int(early.sum()) * e_rmse ** 2 / tot_mse, 4),
"capped": round(int(capped.sum()) * c_rmse ** 2 / tot_mse, 4),
"middle": round(int(mid.sum()) * m_rmse ** 2 / tot_mse, 4)}
# ---- 落盘 ----
bg = tuple(spec["backgroundQ8"])
for tag, name in (("V4", "variance"), ("R4", "relative")):
cglib.save_png(os.path.join(OUT, f"cg_trace_12_3_{name}.png"),
C.build_image(pixels, idx, w, h, bg, res[tag]["est"]))
cglib.save_png(os.path.join(OUT, "cg_trace_12_3_uniform.png"),
C.build_image(pixels, idx, w, h, bg, S[lv.index(rows[0]["uniformN"])]
.astype(np.float64) / rows[0]["uniformN"]))
nmap = np.zeros((h, w), dtype=np.int64)
for k, i in enumerate(idx):
y, x = divmod(i, w)
nmap[y, x] = int(res["V4"]["ns"][k])
cglib.save_gray_png(os.path.join(OUT, "cg_trace_12_3_nmap.png"), nmap)
cols = ["tag", "totalSamples", "samplesPerPixelMean", "stopAtMin", "cappedAtBudget",
"rmseQ16", "maxErrQ16", "selfSemQ16", "honestyRatio", "over3Sigma",
"uniformN", "uniformRmseQ16", "uniformRatio"]
cglib.write_csv(os.path.join(OUT, "cg_trace_12_3_stopping.csv"), cols,
[[r[c] for c in cols] for r in rows])
report = {
"width": w, "height": h, "hitPixels": m, "budget": budget, "checkEvery": step,
"varianceThresholdQ16": var_thr, "relErrorThresholdPermille": rel_thr,
"reference": {"batches": rf["batches"], "samplesPerBatch": rf["batchSamples"],
"ms": ref_ms, "sha256": ref_sha,
"selfSemRmsQ16": round(sem_ref_rms, 4), "matches12_2": ref_sha == REF_SHA},
"ladder": {"ms": lad_ms, "samples": m * budget, "checkpoints": len(lv),
"shadowRays": stats["rays"], "occluded": stats["occluded"]},
"criteria": rows,
"failure": {"earlyStopPixels": int(early.sum()),
"zeroVariancePixels": rows[0]["zeroVarAtMin"],
"earlyStopRmseQ16": round(e_rmse, 4), "earlyStopOver3Sigma": e_over,
"budgetCappedPixels": int(capped.sum()),
"budgetCappedRmseQ16": round(c_rmse, 4),
"middlePixels": int(mid.sum()), "middleRmseQ16": round(m_rmse, 4),
"adaptiveRmseQ16": round(r4["rmse"], 4), "mseShare": share},
}
with open(os.path.join(OUT, "cg_trace_12_3_stopping.json"), "w", encoding="utf-8") as fh:
json.dump(report, fh, ensure_ascii=False, indent=2)
fh.write("\n")
print("== 12.3 采样的停止准则 ==")
print(f"参考解 A（与 12.2 同式）：{rf['batches']} 批 × {rf['batchSamples']} = "
f"{rf['batches'] * rf['batchSamples']} 样本/像素，{ref_ms} ms，"
f"sha256 {ref_sha[:16]}… 与 12.2 一致 {ref_sha == REF_SHA}，"
f"自身 RMS 标准误 {sem_ref_rms:.2f} Q16.16")
print(f"梯子：{m} 像素 × {budget} 样本 = {m * budget} 样本，{len(lv)} 个检查点，"
f"{lad_ms} ms，阴影光线 {stats['rays']}（遮挡 {stats['occluded']}）")
print(f"{'准则':<5} {'总样本':>9} {'样本/像素':>9} {'最小档停':>8} {'全零停':>7} "
f"{'撞上限':>7} {'RMS误差':>9} {'自估RMS':>9} {'诚实比':>7} {'超3σ':>6} "
f"{'均匀N':>6} {'均匀RMS':>9} {'均匀/自适':>9}")
for r in rows:
print(f"{r['tag']:<5} {r['totalSamples']:>9} {r['samplesPerPixelMean']:>9.2f} "
f"{r['stopAtMin']:>8} {r['zeroVarAtMin']:>7} {r['cappedAtBudget']:>7} "
f"{r['rmseQ16']:>9.2f} {r['selfSemQ16']:>9.2f} {r['honestyRatio']:>7.3f} "
f"{r['over3Sigma']:>6} {r['uniformN']:>6} {r['uniformRmseQ16']:>9.2f} "
f"{r['uniformRatio']:>9.3f}")
print(f"失效 (a) 最小档（N = 4）就停：{int(early.sum())} 像素（其中 "
f"{rows[0]['zeroVarAtMin']} 个四样本全零、方差估计为 0），"
f"RMS 误差 {e_rmse:.2f}，其中 {e_over} 个超过自估 3σ")
print(f"失效 (b) 撞预算上限：{int(capped.sum())} 像素，RMS 误差 {c_rmse:.2f}")
print(f"V4 误差三段分解（全体 RMS {r4['rmse']:.2f}）："
f"最小档停 {int(early.sum())} 像素 RMS {e_rmse:.2f} 占 MSE {share['early']:.3f}；"
f"中间档 {int(mid.sum())} 像素 RMS {m_rmse:.2f} 占 {share['middle']:.3f}；"
f"撞上限 {int(capped.sum())} 像素 RMS {c_rmse:.2f} 占 {share['capped']:.3f}")
return 0
if __name__ == "__main__":
sys.exit(main())
准则 B 的 mask 只在最后一档为真，这不是偷懒：np.argmax 在「全为假」时返回 0，会让「预算上限」错误地停在第一个检查点。写成「最后一档为真」正好等价于「N 到预算就停」。
跑一遍：
# snippet: 本机实跑（工作区 lantern/）
cd lantern
/opt/homebrew/bin/python3 sim/cg_trace_12_3_stopping.py
== 12.3 采样的停止准则 ==
参考解 A（与 12.2 同式）：16 批 × 64 = 1024 样本/像素，109813 ms，sha256 0ac7328aaeb5c669… 与 12.2 一致 True，自身 RMS 标准误 28.45 Q16.16
梯子：2734 像素 × 256 样本 = 699904 样本，64 个检查点，28795 ms，阴影光线 540826（遮挡 7882）
准则 总样本 样本/像素 最小档停 全零停 撞上限 RMS误差 自估RMS 诚实比 超3σ 均匀N 均匀RMS 均匀/自适
V4 328084 120.00 526 498 554 187.39 59.97 3.125 220 120 87.34 0.466
V32 369976 135.32 406 384 663 72.95 62.75 1.163 87 132 84.55 1.159
R4 107424 39.29 1119 498 314 307.82 140.11 2.197 361 36 158.62 0.515
R32 200844 73.46 2142 384 459 98.89 91.98 1.075 76 72 111.38 1.126
B 699904 256.00 0 0 0 64.13 57.18 1.122 34 256 64.13 1.000
失效 (a) 最小档（N = 4）就停：526 像素（其中 498 个四样本全零、方差估计为 0），RMS 误差 397.23，其中 187 个超过自估 3σ
失效 (b) 撞预算上限：554 像素，RMS 误差 120.20
V4 误差三段分解（全体 RMS 187.39）：最小档停 526 像素 RMS 397.23 占 MSE 0.865；中间档 1654 像素 RMS 54.99 占 0.052；撞上限 554 像素 RMS 120.20 占 0.083
实测：停在哪、花多少、错多少
| 准则 | 总样本 | 样本/像素 | 最小档停 | 全零停 | 撞上限 | RMS 误差 | 自估 RMS | 诚实比 | 超 3σ | 均匀 N | 均匀 RMS | 均匀/自适 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| V4 | 328084 | 120.00 | 526 | 498 | 554 | 187.39 | 59.97 | 3.125 | 220 | 120 | 87.34 | 0.466 |
| V32 | 369976 | 135.32 | 406 | 384 | 663 | 72.95 | 62.75 | 1.163 | 87 | 132 | 84.55 | 1.159 |
| R4 | 107424 | 39.29 | 1119 | 498 | 314 | 307.82 | 140.11 | 2.197 | 361 | 36 | 158.62 | 0.515 |
| R32 | 200844 | 73.46 | 2142 | 384 | 459 | 98.89 | 91.98 | 1.075 | 76 | 72 | 111.38 | 1.126 |
| B | 699904 | 256.00 | 0 | 0 | 0 | 64.13 | 57.18 | 1.122 | 34 | 256 | 64.13 | 1.000 |
（单位：RMS 误差与自估 RMS 为 Q16.16；「均匀/自适」= 同预算均匀采样的 RMS ÷ 自适应的 RMS，小于 1 表示自适应更差。）
这张表里有两件事值得停一下。
第一，「诚实比」列。判停用的是流内自估的 s/√N，与对参考解的实际误差比一比，比值接近 1 才说明「自己报的不确定度」可信。预算上限这一条（B）的比值是 1.122——在固定 N 下自估是诚实的。而 V4 的比值是 3.125：它自己以为误差 59.97，实际是 187.39，差了 3 倍。自适应会让自估变得不诚实，因为判停本身与样本方差相关——方差估小了的像素会最早停下来，而它们正是错得最狠的那些。「超 3σ」列就是这个现象的直接计数：V4 有 220 个像素的实际误差超过自估的 3 倍。
第二，「均匀/自适」列。把每条准则的总样本数平摊成均匀 N 再跑一遍：V4 是 0.466、R4 是 0.515，都小于 1——在这份配置下，用 4 个样本做最小档的自适应采样比均匀采样更差，而且是差一倍。把最小档提到 32（V32 / R32），比值变成 1.159 和 1.126，自适应才开始赢。自适应不是免费的，它的收益完全取决于最小样本数这个参数选得对不对。
两种失效模式
失效 (a)：最小档就判停。 526 个像素在 N = 4 就停了，其中 498 个是四样本全零——四个样本的贡献都是 0，于是 Σx² = 0、s² = 0、s/√N = 0，判据「≤ 64」立刻成立。四个样本为什么全零？脚本没有区分两种可能：四次采样点都被几何挡住，或者这个像素本身就处在背光面（n·l ≤ 0）——无论哪一种，这一轮的四次贡献都是 0。而零方差的观测是信息量最低的观测，朴素的自适应准则却把它当成了「已经收敛」。
这 526 个像素的 RMS 误差是 397.23，其中 187 个超过自估的 3σ。把 V4 全体 187.39 的 RMS 按 MSE 分解：
| 分组 | 像素数 | 该组 RMS | 占全体 MSE |
|---|---|---|---|
| 最小档停 | 526 | 397.23 | 0.865 |
| 中间档 | 1654 | 54.99 | 0.052 |
| 撞上限 | 554 | 120.20 | 0.083 |
全体误差的 86.5% 只来自那 526 个过早停止的像素。 中间那 1654 个像素的 RMS 只有 54.99，比预算上限那一条（64.13）还好。也就是说：V4 这条准则大部分时候工作得很好，坏就坏在那一小撮「四样本全零」的像素上。修法也直白——提高最小样本数，这正是 V32 把诚实比从 3.125 拉回 1.163 的原因。
失效 (b)：撞预算上限。 554 个像素跑满 256 个样本仍没达到阈值，它们的 RMS 误差是 120.20。这些是真正的高方差区域（光源掠射、部分遮挡的软阴影边缘）。对它们，s/√N 的收敛速度不足以在 256 个样本内降到 64 Q16.16——自适应准则对这类区域完全失效，只能靠预算兜底。这也是为什么准则 B 必须存在：没有它，程序会在这里无限循环。
P 层：判停规则也要跨语言一致
判停发生在浮点里（比较 s/√N 与阈值），但它的输入是整数矩。P 层用 BigInt 重算整数矩、在 TS 侧从矩还原 s/√N 再判停，断言四条规则的停止档位与 Python 侧逐个相同：
# snippet: core/test/cg_trace_12_1_noise.test.ts（节选：四条规则的判停断言）
# verify: node --test 'core/test/cg_trace_12_1_noise.test.ts'
const rules: [string, "V" | "R", number, number][] = [
["V4", "V", cfg.varianceThresholdQ16, 0],
["V32", "V", cfg.varianceThresholdQ16, I32],
["R4", "R", cfg.relErrorThresholdPermille, 0],
["R32", "R", cfg.relErrorThresholdPermille, I32],
];
for (const [tag, kind, thr, minIdx] of rules) {
const stop = firstStop(got.lv, got.sum, got.sumSq, kind, thr, minIdx);
assert.equal(stop, want[tag], `像素 ${i} 规则 ${tag}`);
}
注意最后一行的约定：全程不达标时返回最后一档，而不是 0 或 -1。这样「预算上限」这条兜底规则在两侧的语义完全一致——Python 侧的 apply_criterion 也是这个约定。
小结
- 三条准则共用一条流，64 个检查点回放判停，三条看到的是同一批随机数。
- 固定 N 下自估的标准误是诚实的（B 的诚实比 1.122），但朴素自适应会把它变成 3.125。
- 最小档 N = 4 时，526 个像素过早停止（498 个四样本全零），独占全体 MSE 的 86.5%。
- 同预算下 V4 / R4 比均匀采样更差（0.466 / 0.515）；最小档提到 32 才反超（1.159 / 1.126）。
- 高方差区域到不了阈值，554 个像素撞上预算上限——兜底规则必须存在。
下一章换一个角度：前面一直在跟随机误差打交道，而定点数还会带来另一类误差——求交的容差边界。
阅读导航：上一节：12.2 参考解与误差 · 下一节：13.1 求交的容差边界 。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。