12.2 参考解与误差
站内 /graphics-gpu-profiling-debugging/ 那篇做的是 GPU 剖析与调试：帧时间、事件时间轴、绘制调用，量的是实时管线的开销。本章量的是离线渲染的误差——而误差这个词有个前提常被跳过：你得先有一个「对的」东西。本节就专门盯住这个前提：参考解自己有多少不确定度？
参考解也是一个估计量
上一节用「8 批 × 64 = 512 样本/像素」的结果当参考。可它同样是蒙特卡洛估计量，同样带方差。用一个带方差的东西去量另一个带方差的东西，误差里就混着两份不确定性。所以先把参考解做大：16 条独立批 × 64 样本 = 1024 样本/像素。
「批」是造不确定度的工具：16 条批各自独立（起点 seed + batch·batchStride + index·streamStride，batchStride = 2^52 远大于 streamStride = 2^40，流之间永不重叠），于是得到 16 个独立的批均值，它们的标准差除以 sqrt(16) 就是参考解自己的标准误：
# snippet: sim/cg_trace_12_2_reference.py（节选：批均值与对半拆）
# verify: /opt/homebrew/bin/python3 sim/cg_trace_12_2_reference.py
perA, totA = C.batch_sums(pixels, idx, light, tris, bvh, nb, ns, rf["seed"],
rf["batchStride"], rf["streamStride"])
refA = totA // (nb * ns)
bmA = perA / ns
semA = bmA.std(axis=0, ddof=1) / math.sqrt(nb)
sem_rms = float(np.sqrt((semA ** 2).mean())) # 图像级口径：RMS
# 对半拆：A₁ 与 A₂ 的实测差
sub1 = perA[:half].sum(axis=0) // (half * ns)
sub2 = perA[half:].sum(axis=0) // (half * ns)
d12_rms = float(np.sqrt(((sub1 - sub2).astype(np.float64) ** 2).mean()))
注意最后一行：图像级比较取的是标准误的 RMS（np.sqrt((semA ** 2).mean())），不是标准误的算术平均。这不是风格问题——逐像素的标准误相差几个数量级（最大 303 Q16.16，最小接近 0），取算术平均会把高方差像素的贡献摊平，算出来的「几倍标准误」就不再是理论值 2.000 对应的那个量。这也正是下面四条量法都要跟理论值对一遍的原因：口径对了，2.045 与 2.219 才会分别落在 2.000 与 2.236 附近。
四条互相独立的量法
参考解的不确定度不能只靠一条推理链，本节用了四条彼此独立的路：
| 量法 | 做法 | 理论比值 | 实测 |
|---|---|---|---|
| 批均值标准误 | 16 批均值的 s/√16 | — | 相对 0.005729 |
| 对半拆 | A₁(8 批) 与 A₂(8 批) 相减 | 2.000 | 2.045 |
| 换种子族 | 参考解 B（4 批，另一套种子）与 A 相减 | 2.236 | 2.219 |
| 确定性求积 | 全程无遮挡像素上的 33×33 中点网格 | — | 与 A 最大差 63 Q16.16 |
对半拆的理论值 2.000 来自「两个互不重叠的 8 批估计之差，方差是各自的 2 倍」；换种子族的 2.236 = sqrt(1 + 1/4) 来自「A 是 16 批、B 是 4 批」。两条都实测吻合到 3% 以内——这说明「标准误」这个口径本身是可信的，不是自说自话。
#!/usr/bin/env python3
# file: lantern/sim/cg_trace_12_2_reference.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_trace_12_2_reference.py
"""12.2 参考解与误差：参考解自己有多少不确定度。
「误差」只在**有一个参考解**的时候才有意义，而参考解本身也是估计量。四条量法：
1. **参考解 A**：16 条独立批 × 64 样本 = 1024 样本/像素，批均值的标准误 s/√16 就是
参考解自己的不确定度。图像级比较取 **RMS**——逐像素标准误相差几个数量级，取算术
平均会把高方差像素的贡献摊平，比值就不再是「几倍标准误」该有的量级。
2. **对半拆**：16 批拆成 8 + 8 得两条互不重叠的子参考 A₁、A₂。A₁ − A₂ 的实测 RMS 差
应当接近 √2 倍标准误——这是对「标准误」这个口径本身的检验。
3. **另一族种子**：换一套种子起点造参考解 B，与 A 相减。两次独立构造的「真值」差
多少，就是「参考解不是绝对真值」的直接证据（理论 √(1+1/4) = 2.236 倍标准误）。
4. **确定性求积**：在**全程无遮挡**的像素上，面光源的期望可以用矩形上的确定性网格
积分算出来，这条路与蒙特卡洛完全无关。判据是**可见性恒定**而不是贡献恒定——贡献
会随光源点位置连续变化（n·l、距离衰减、归一化都在变），恒定这个条件永远不成立。
产物：out/cg_trace_12_2_reference.json、out/cg_trace_12_2_error.csv、
out/cg_trace_12_2_ref.png、out/cg_trace_12_2_refsem.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")
def main():
os.makedirs(OUT, exist_ok=True)
spec = cglib.load_spec()
cfg = C.load_area_cfg()
rf = cfg["reference"]
w, h = rf["width"], rf["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 = len(idx)
nb, ns = rf["batches"], rf["batchSamples"]
# ---- 1. 参考解 A + 批均值给出的标准误 ----
t0 = time.time()
perA, totA = C.batch_sums(pixels, idx, light, tris, bvh, nb, ns, rf["seed"],
rf["batchStride"], rf["streamStride"])
refA_ms = int((time.time() - t0) * 1000)
refA = totA // (nb * ns)
bmA = perA / ns
semA = bmA.std(axis=0, ddof=1) / math.sqrt(nb)
lumA = float(refA.mean())
sem_rms = float(np.sqrt((semA ** 2).mean())) # 图像级口径：RMS
bad = int((semA.mean(axis=1) > refA.mean(axis=1) * 0.01).sum())
# ---- 2. 对半拆：A₁ 与 A₂ 的实测差 ----
half = nb // 2
sub1 = perA[:half].sum(axis=0) // (half * ns)
sub2 = perA[half:].sum(axis=0) // (half * ns)
d12_rms = float(np.sqrt(((sub1 - sub2).astype(np.float64) ** 2).mean()))
# ---- 3. 另一族种子：参考解 B ----
t0 = time.time()
_, totB = C.batch_sums(pixels, idx, light, tris, bvh, rf["altBatches"], ns,
rf["altSeed"], rf["batchStride"], rf["streamStride"])
refB_ms = int((time.time() - t0) * 1000)
refB = totB // (rf["altBatches"] * ns)
dAB_rms = float(np.sqrt(((refA - refB).astype(np.float64) ** 2).mean()))
# ---- 4. 误差指标表：独立扫掠 vs 参考解 A ----
levels = cfg["experiment"]["levels"]
t0 = time.time()
S, _ = C.sweep(pixels, idx, light, tris, bvh, levels, max(levels),
cfg["experiment"]["seed"], rf["streamStride"])
sweep_ms = int((time.time() - t0) * 1000)
rows = []
for n in levels:
err = S[n] / n - refA
q = C.rel_stats(np.abs(err), refA)
rmse_q16 = float(np.sqrt((err.astype(np.float64) ** 2).mean()))
rmse8 = float(np.sqrt((((err * 255.0 / 65536.0) ** 2).mean())))
rows.append({"n": n, "rmseQ16": round(rmse_q16, 4),
"relRmse": round(rmse_q16 / lumA, 8), "rmse8bit": round(rmse8, 6),
"psnrDb": round(float("inf") if rmse8 == 0 else 20 * math.log10(255 / rmse8), 4),
"p50": q["p50"], "p90": q["p90"], "p99": q["p99"], "max": q["max"],
"overRefSem": round(rmse_q16 / sem_rms, 4)})
# ---- 5. 确定性求积：只取「全程无遮挡」的像素 ----
cand = []
for k in range(m):
px = pixels[idx[k]]
if all(C.light_visible(px, lp, light, tris, bvh) for lp in C.light_grid(light, 9)):
cand.append(k)
n_cand = len(cand)
cand = cand[::max(1, n_cand // 64)][:64]
quad_rows = []
for k in cand:
q17 = C.quadrature(pixels[idx[k]], light, tris, bvh, 17)
q33 = C.quadrature(pixels[idx[k]], light, tris, bvh, 33)
mc = [int(refA[k, ch]) for ch in range(3)]
quad_rows.append({"pixel": idx[k], "refA": mc, "quad17": q17, "quad33": q33,
"quadDrift": int(max(abs(q33[c] - q17[c]) for c in range(3))),
"mcVsQuad": int(max(abs(mc[c] - q33[c]) for c in range(3)))})
quad_drift = max((r["quadDrift"] for r in quad_rows), default=0)
quad_gap = max((r["mcVsQuad"] for r in quad_rows), default=0)
ref_img = C.digest_light_only(pixels, idx, w, h, tuple(spec["backgroundQ8"]), refA)
sem_map = np.zeros((h, w), dtype=np.int64)
for k, i in enumerate(idx):
sem_map[divmod(i, w)] = int(semA[k].max())
cglib.save_gray_png(os.path.join(OUT, "cg_trace_12_2_refsem.png"), sem_map)
cglib.write_csv(os.path.join(OUT, "cg_trace_12_2_error.csv"),
["n", "rmseQ16", "relRmse", "rmse8bit", "psnrDb", "p50", "p90",
"p99", "max", "overRefSem"],
[[r["n"], r["rmseQ16"], r["relRmse"], r["rmse8bit"], r["psnrDb"],
r["p50"], r["p90"], r["p99"], r["max"], r["overRefSem"]] for r in rows])
report = {
"width": w, "height": h, "hitPixels": m,
"referenceA": {"batches": nb, "samplesPerBatch": ns, "totalSamples": nb * ns,
"ms": refA_ms, "lumaQ16": round(lumA, 4),
"relSem": round(sem_rms / lumA, 8), "maxSemQ16": int(semA.max()),
"pixelsSemOver1pct": bad, "sha256": ref_img},
"split": {"halfBatches": half, "rmsDiffQ16": round(d12_rms, 4),
"relDiff": round(d12_rms / lumA, 8),
"overSem": round(d12_rms / sem_rms, 4), "theory": 2.0},
"referenceB": {"batches": rf["altBatches"], "ms": refB_ms,
"rmsDiffQ16": round(dAB_rms, 4),
"relDiff": round(dAB_rms / lumA, 8),
"overSem": round(dAB_rms / sem_rms, 4), "theory": 2.2360679},
"sweepMs": sweep_ms, "rows": rows,
"quadrature": {"gridA": 17, "gridB": 33, "candidates": n_cand,
"pixels": len(quad_rows), "maxDriftQ16": quad_drift,
"maxMcVsQuadQ16": quad_gap,
"relGap": round(quad_gap / max(1.0, lumA), 8),
"rows": quad_rows[:8]},
}
with open(os.path.join(OUT, "cg_trace_12_2_reference.json"), "w", encoding="utf-8") as fh:
json.dump(report, fh, ensure_ascii=False, indent=2)
fh.write("\n")
print("== 12.2 参考解与误差 ==")
print(f"参考解 A：{nb} 批 × {ns} = {nb * ns} 样本/像素，{refA_ms} ms，"
f"平均亮度 {lumA:.2f} Q16.16，sha256 {ref_img[:16]}…")
print(f"参考解自身不确定度：相对标准误 {sem_rms / lumA:.6f}；"
f"相对标准误 > 1% 的像素 {bad} / {m}；逐像素最大 {int(semA.max())} Q16.16")
print(f"对半拆 A₁ vs A₂：RMS 差 {d12_rms:.2f} Q16.16 = {d12_rms / sem_rms:.3f} × "
f"标准误（理论 2.000）")
print(f"另一族参考 B（{rf['altBatches']} 批）：RMS 差 {dAB_rms:.2f} Q16.16 = "
f"{dAB_rms / sem_rms:.3f} × 标准误（理论 2.236），{refB_ms} ms")
print(f"{'N':>5} {'RMSE(Q16)':>11} {'相对RMSE':>11} {'RMSE(8位)':>11} {'PSNR dB':>9} "
f"{'p99':>8} {'max':>8} {'/参考标准误':>11}")
for r in rows:
print(f"{r['n']:>5} {r['rmseQ16']:>11.2f} {r['relRmse']:>11.6f} {r['rmse8bit']:>11.6f} "
f"{r['psnrDb']:>9.4f} {r['p99']:>8.4f} {r['max']:>8.4f} {r['overRefSem']:>11.3f}")
print(f"确定性求积（全程无遮挡像素 {n_cand} 个，取 {len(quad_rows)} 个）："
f"17×17 vs 33×33 最大漂移 {quad_drift} Q16.16；参考解 A vs 求积 最大差 "
f"{quad_gap} Q16.16（相对 {quad_gap / lumA:.6f}）")
return 0
if __name__ == "__main__":
sys.exit(main())
跑一遍：
# snippet: 本机实跑（工作区 lantern/）
cd lantern
/opt/homebrew/bin/python3 sim/cg_trace_12_2_reference.py
== 12.2 参考解与误差 ==
参考解 A：16 批 × 64 = 1024 样本/像素，108083 ms，平均亮度 4966.30 Q16.16，sha256 0ac7328aaeb5c669…
参考解自身不确定度：相对标准误 0.005729；相对标准误 > 1% 的像素 540 / 2734；逐像素最大 303 Q16.16
对半拆 A₁ vs A₂：RMS 差 58.18 Q16.16 = 2.045 × 标准误（理论 2.000）
另一族参考 B（4 批）：RMS 差 63.14 Q16.16 = 2.219 × 标准误（理论 2.236），27390 ms
N RMSE(Q16) 相对RMSE RMSE(8位) PSNR dB p99 max /参考标准误
1 906.70 0.182570 3.527960 37.1803 0.2569 1.5384 31.868
2 662.65 0.133428 2.578351 39.9040 0.2057 1.1029 23.290
4 469.39 0.094515 1.826391 42.8989 0.1321 1.1082 16.498
8 342.88 0.069040 1.334127 45.6269 0.1139 0.4478 12.051
16 239.94 0.048313 0.933586 48.7277 0.0843 0.2881 8.433
32 171.74 0.034582 0.668250 51.6320 0.0583 0.2923 6.036
64 118.28 0.023816 0.460218 54.8715 0.0439 0.2632 4.157
128 85.18 0.017152 0.331445 57.7226 0.0286 0.2245 2.994
256 63.86 0.012859 0.248478 60.2250 0.0219 0.1349 2.244
确定性求积（全程无遮挡像素 1754 个，取 64 个）：17×17 vs 33×33 最大漂移 8 Q16.16；参考解 A vs 求积 最大差 63 Q16.16（相对 0.012685）
误差指标表：RMSE、相对误差、分位数
独立扫掠（同一 seed/stride，与参考解不共享随机数）与参考解 A 相减，得到九个档位的误差指标：
| N | RMSE (Q16.16) | 相对 RMSE | RMSE (8 位) | PSNR (dB) | p99 | max | / 参考标准误 |
|---|---|---|---|---|---|---|---|
| 1 | 906.70 | 0.182570 | 3.527960 | 37.1803 | 0.2569 | 1.5384 | 31.868 |
| 8 | 342.88 | 0.069040 | 1.334127 | 45.6269 | 0.1139 | 0.4478 | 12.051 |
| 64 | 118.28 | 0.023816 | 0.460218 | 54.8715 | 0.0439 | 0.2632 | 4.157 |
| 256 | 63.86 | 0.012859 | 0.248478 | 60.2250 | 0.0219 | 0.1349 | 2.244 |
最后一列是本节的核心结论。N = 256 时，扫描解与参考解的差是参考解自身标准误的 2.244 倍——也就是说，再往上加样本，误差也不会低于这个量级，因为分母（参考解）本身就在抖。1024 样本/像素的参考解，只够支撑到「误差约为它自身标准误的两倍」这个结论，想再往下压，必须先把参考解做大。
还要注意 p99 与 max 的关系：N = 256 时 p99 是 0.0219，max 是 0.1349，差 6 倍。用 RMSE 说「误差 1.3%」和用 max 说「最坏的像素差 13.5%」，两个都对，但传达的信息完全不同。
为什么参考解不是绝对真值
「参考解」这个词很容易让人以为它是真值。它不是，理由有三条：
- 它自己带标准误。本节实测：相对标准误 0.005729，而且逐像素差别极大——2734 个命中像素里有 540 个的相对标准误超过 1%，最大的像素标准误是 303 Q16.16。
- 换一套种子，它就变了。参考解 B（4 批，另一套种子）与 A 的 RMS 差是 63.14 Q16.16，相当于图像平均亮度的 1.27%。两次独立构造的「真值」自己就差了 1.27%，这就是「参考解不是绝对真值」最直接的证据。
- 它对场景量程敏感。参考解是 Q16.16 定点数，每批求和后整除取均值，
//的截断是确定性的、但也是偏的。这不是主要误差来源（远小于统计误差），但它意味着参考解与真值之间永远差着一点系统性的东西。
所以本章的所有说法都写成「与参考解 A 的差」，而不是「误差」。这个措辞不是谦虚，是口径：换一个参考解，数字就会变。
确定性求积：一条与蒙特卡洛无关的路
前三条量法都是「用蒙特卡洛检验蒙特卡洛」。有没有一条完全不同的路？有——在全程无遮挡的像素上，面光源的照明就是一个矩形上的确定性积分，用中点网格求积算出来，与随机数无关。
判据必须是可见性恒定，不是「贡献恒定」：贡献会随光源点位置连续变化（n·l、距离衰减、归一化都在变），恒定这个条件永远不成立。脚本在 9×9 网格上筛出全程无遮挡的像素 1754 个，从里面取 64 个，各算 17×17 与 33×33 两种网格：
- 17×17 与 33×33 的最大漂移 8 Q16.16——网格加密一倍，求积结果几乎不动，说明这两档都已经收敛到积分值附近；
- 参考解 A 与 33×33 求积的最大差 63 Q16.16（相对 0.012685）。
63 这个数与参考解自己的标准误同量级（像素级最大 303，图像级 RMS 28.45）。换句话说：一条与蒙特卡洛毫无关系的确定性路径，给出的结果落在蒙特卡洛参考解自身的不确定度之内。这是本节最强的一条证据——不是「两个估计量碰巧接近」，而是两条独立的方法在同一量级上互相印证。
P 层：整数矩还原方差
跨语言这一层，一致性判据落在整数矩上。TypeScript 侧同样从 Σx 与 Σx² 还原样本方差和标准误，浮点只在「同一串运算顺序」这一层参与：
# snippet: core/src/lantern/cg_trace_12_metric.ts（节选：从整数矩还原方差）
# verify: node --test 'core/test/cg_trace_12_1_noise.test.ts'
/** 样本方差 s² = (n·Σx² − (Σx)²) / (n(n−1))，由整数矩还原。 */
export function sampleVariance(s: bigint, sq: bigint, n: number): number {
const nn = n;
return Math.max((Number(sq) * nn - Number(s) * Number(s)) / (nn * (nn - 1)), 0);
}
/** 均值的标准差 s/√n，由整数矩还原。 */
export function semQ16(s: bigint, sq: bigint, n: number): number {
return Math.sqrt(sampleVariance(s, sq, n) / n);
}
Math.max(..., 0) 那一层保护是必要的：n·Σx² − (Σx)² 在整数运算下恒非负，但转成 Number 之后可能因浮点舍入变成极小的负数，开方会得到 NaN。
小结
- 参考解也是估计量，它的不确定度必须被量出来，否则「误差」这个词没有分母。
- 四条量法互相印证：批均值标准误、对半拆 2.045、换种子族 2.219、确定性求积最大差 63 Q16.16。
- 1024 样本/像素的参考解，相对标准误 0.005729，但540 / 2734 个像素超过 1%。
- N = 256 的扫描解与参考解的差是参考解标准误的 2.244 倍——参考解是当前的地板。
下一节把误差反过来用：既然知道每个像素的误差量级，什么时候可以停止采样？
阅读导航：上一节：12.1 噪声与收敛的度量 · 下一节：12.3 采样的停止准则 。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。