7.1 概率密度与权重
站内 /graphics-real-time-gi/ 那篇做的是实时 GI:屏幕空间探针、光照贴图、辐照度缓存,目标是把间接光塞进 16 毫秒的预算里。本章换一条路——把采样本身当成研究对象,问的是方差从哪里来、怎么把它压下去。前一章已经把渲染方程写成了可估计的形式,这一章开始动它的旋钮。第一个旋钮就是概率密度。
估计量只有一个形状
蒙特卡洛把积分变成平均值:
$$\langle I \rangle = \frac{1}{N}\sum_{i=1}^{N}\frac{f(x_i)}{p(x_i)}$$
x_i 从概率密度 p 抽出,f 是待积分的量。这个式子只有一个形状——f/p。p 换成什么,估计量的期望不变(只要 p 在 f 非零处也非零),变的只有方差。
于是 p 要满足两件事,而它们经常被混为一谈:
- 归一:
∫ p = 1。这是「它是不是一个概率密度」的问题。漏掉归一常数,估计量会整体偏一个固定倍数——不随机、不收敛,就是错的。 - 形状:
p在f大的地方也要大。这是「重要性采样」这个名字的来源。当p ∝ f时f/p恒为常数,方差恰好为 0。
本节把这两件事分开量。归一在半球上验(三个 pdf,两个应得 1、一个应得 π);形状在一维区间上验(同一个积分 I = ∫₀¹ 2x dx = 1,三条 pdf,方差差两个量级)。最后把 fdiv 在定点下的截断误差和除零边界摆出来——这两件事都会在后面的 MIS 里咬人。
归一:在半球上估 ∫ p dω
余弦加权采样是路径追踪里最常用的 pdf:p_c(ω) = cosθ/π,θ 是方向与法线的夹角。1/π 这个常数很容易被漏掉——漏掉之后形状看起来「更简单」,代码照跑不误,只是整体偏了 π 倍。
要验归一,得用一个已知的采样器去估另一个 pdf 的积分。这里用半球上的均匀采样:p_u(ω) = 1/(2π),cosθ 均匀、方位角在单位圆上均匀。于是
$$\int_\Omega p(\omega),\mathrm{d}\omega ;=; \mathbb{E}_{\omega \sim p_u}\bigl[p(\omega)\cdot 2\pi\bigr]$$
三个被估的 pdf:cosθ/π(应 1)、1/(2π)(应 1)、以及漏掉常数的 cosθ(应 π)。1048576 个均匀样本,实跑结果:
被估的 p(ω) | 估得 ∫ p dω | 理论值 | 相对误差 |
|---|---|---|---|
cosθ/π | 1.000334 | 1 | +0.033% |
1/(2π) | 0.999954 | 1 | −0.005% |
cosθ(漏常数) | 3.142802 | π = 3.141593 | +0.038% |
三行里最值得看的是第三行。它不是「错了」——cosθ 确实是一个处处非负的函数,∫ cosθ dω = π 也是对的;错的是把它当概率密度用。只要它出现在分母上,估计量就被放大 π 倍。
第一行与第三行的比值是 3.142802 / 1.000334 = 3.14175,与 π 只差 5×10⁻⁵。这不是巧合:同一条采样流下,两个估计量的被积函数只差一个常数 1/π,所以比值精确等于 π(扣除定点截断)。归一常数漏没漏,用一次比值就能判。
第二行的 0.999954 不是采样噪声,是纯定点截断:fdiv(ONE, 2π) 把 1/(2π) 截断到 65533/65536,再乘 2π 就回不到 1。它给出了这一类测量的精度地板——大约 5×10⁻⁵,比蒙特卡洛噪声(第一、三行的 3×10⁻⁴)还小一个量级。
均匀半球采样本身也要用整数写出来。cosθ = u 直接取一个 rng_q16 就是均匀的,麻烦的是方位角——要在单位圆上取均匀点,又不能用 sin/cos。做法是单位圆上的拒绝采样:在 [−1,1]² 里抽点,落在单位圆外就丢掉,落进来再归一化。返回的方向是
# snippet: sim/cg_trace_07_1_pdf_weight.py(节选:均匀半球采样)
c = u # cosθ 均匀
s = fsqrt(ONE - fmul(c, c)) # sinθ = sqrt(1 − cos²θ)
# 方位角:单位圆上的拒绝采样,再归一化
cx, cy = fdiv(x, r), fdiv(y, r)
return (fmul(s, cx), fmul(s, cy), c)
拒绝采样有一个容易被忽视的后果:被拒绝的次数由随机流决定。同一条 SplitMix64 流上,Python 与 TypeScript 必须拒绝同样多的样本、在同样的位置拒绝,否则后面所有样本全部错位。这正是本节 P 层对照缓冲只取 16 个样本、却仍然要求逐字节一致的原因——采样的正确性不只在于分布对不对,还在于消耗随机数的顺序。
形状:匹配 pdf 让方差归零
归一只是及格线。真正决定方差的是形状。取一个能被手算的积分:
$$I = \int_0^1 2x,\mathrm{d}x = 1, \qquad f(x) = 2x$$
三条 pdf 依次逼近 f 的形状:
| 逆 CDF | 权重 f/p | 权重取值范围 | |
|---|---|---|---|
均匀 p = 1 | x = u | 2x | [0, 2] |
匹配 p = 2x | x = √u | 2x / 2x = 1 | 恒为 1 |
反匹配 p = 2(1−x) | x = 1 − √(1−u) | x/(1−x) | [0, 255] |
三条 pdf 的逆 CDF 都只用平方根——fsqrt(x) = isqrt(x << 16),整数、无三角函数,和底座里余弦采样用的是同一个原语。匹配那条的权重化简后恒等于 1,所以方差恰好为 0:每个样本都给出精确的积分值,不需要平均。
256 条独立子流 × 4096 样本,实跑结果:
| 均值 | 相对偏差 | 相对标准差 | 最大权重 | 除零次数 | |
|---|---|---|---|---|---|
均匀 1 | 1.000166 | +0.017% | 0.009193 | 2.0000 | 0 |
匹配 2x | 0.999987 | −0.001% | 0.000056 | 1.0000 | 14 |
反匹配 2(1−x) | 0.999981 | −0.002% | 0.042460 | 255.0000 | 0 |
三个均值都在 1 附近(相对偏差 ≤ 0.02%),这是无偏的直接证据:三条 pdf 形状差得很远,收敛到的却是同一个数。变的是方差——匹配 pdf 的相对标准差 5.6×10⁻⁵,比均匀的 9.2×10⁻³ 小两个量级。
反匹配那行更值得停一下。它的相对标准差只有 4.62 倍于均匀(0.0425 对 0.0092),但最大权重是 255。原因是 x/(1−x) 的分布在 [0,1] 上极度偏左:绝大多数样本落在 x 很小的地方,权重接近 0;只有极少数样本抽到 x → 1,权重暴涨。标准差被这些稀有样本抬高了 4.6 倍,而最大值暴露了真正的危险——一个样本就能贡献 255 倍的能量,这就是「反匹配」在画面上的表现:稀疏的亮斑。
匹配那行的「除零 14」也不是 bug。p = 2x 在 x < 1/65536 时被定点截断到 0,此时 f = 2x 同样是 0,0/0 无定义。1048576 个样本里出现 14 次,与 1/65536 的概率吻合。脚本的做法是跳过这些样本,这会给均值引入约 1.3×10⁻⁵ 的偏差——正好是实测的 −0.001%(0.999987 对 1)。
定点:fdiv 的截断与除零边界
权重 f/p 落到定点上就是一次 fdiv。底座里的实现是向零截断的整数除法:
# snippet: sim/cglib.py(节选:fdiv)
def fdiv(a, b):
if b == 0:
raise ZeroDivisionError("fdiv by zero")
q = (abs(a) << 16) // abs(b)
return -q if (a < 0) != (b < 0) else q
精确的 Q16.16 值是 f·65536/p(一个有理数),fdiv 给出的是它的整数部分。于是绝对误差恒小于 1,但相对误差会被权重的大小放大:
f | p | fdiv 得 q16 | 精确值 f·65536/p | 相对误差 |
|---|---|---|---|---|
| 65536 | 65536 | 65536 | 65536.000000 | 0 |
| 3 | 131072 | 1 | 1.500000 | −33.3% |
| 65536 | 196608 | 21845 | 21845.333333 | −1.5×10⁻⁵ |
| 5 | 458752 | 0 | 0.714286 | −100% |
| 17 | 327680 | 3 | 3.400000 | −11.8% |
| 65536 | 7 | 613566756 | 613566756.571429 | −9.3×10⁻¹⁰ |
权重在 1 量级时相对误差是 10⁻⁵ 量级——这就是上一节那个 5×10⁻⁵ 精度地板的来源。权重掉到 0.7 时,fdiv 直接给出 0:整个样本的贡献消失。所以不能拿「权重很小」当无害——在定点下它等于「这个样本不存在」。
另一半是 b == 0 时抛异常。这一条在本节只是被跳过,但到 7.3 会变成真正的分支:光源 pdf 在掠射角处会被截断到 0,fdiv 一旦没有 guard 就整条链路崩掉。定点代码里每个 fdiv 的分母都要先问一句「它可能是 0 吗」。
跑一次
# snippet: 本机实跑(工作区 lantern/)
cd lantern
/opt/homebrew/bin/python3 sim/cg_trace_07_1_parity.py
驱动脚本本身只有四十来行,四块计算都在 cg_trace_07_1_pdf_weight 里:
#!/usr/bin/env python3
# file: lantern/sim/cg_trace_07_1_parity.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_trace_07_1_parity.py
"""7.1 的紧凑驱动:调用 cg_trace_07_1_pdf_weight 的四个函数,打印正文引用的全部数字。
四块输出对应正文四件事:半球归一估计、三条 1D pdf 的均值/方差/最大权重、
fdiv 的定点截断误差、以及 P 层逐字节对照的整数缓冲。
"""
import hashlib
import os
import sys
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import cg_trace_07_1_pdf_weight as m # noqa: E402
ONE = 65536
def main():
norm = m.norm_estimates()
shape = m.shape_study()
wtab = m.weight_table()
buf, vals = m.parity_values()
print("== 7.1 概率密度与权重 ==")
print(f"归一(半球上估 ∫ p dω,{m.N_NORM} 个均匀样本,采样 pdf = 1/2π):")
for k, want in (("cosine", "1"), ("uniform", "1"), ("unnormalized", "π")):
print(f" p = {k:<12}:{norm[k] / ONE:.6f}(应 {want})")
print(f"形状(f(x) = 2x,I = 1;{m.N_TRIALS} 子流 × {m.N_SAMPLES} 样本):")
for k in ("uniform", "matched", "anti"):
d = shape[k]
print(f" p = {k:<8}:均值 {d['meanFloat']:.6f}|相对标准差 {d['relStd']:.6f}"
f"|最大权重 {d['maxWeightFloat']:.4f}|除零 {d['zeroPdfCount']}")
print("定点权重表(fdiv vs 精确 Q16.16 值 f·65536/p):")
for r in wtab:
print(f" f={r['f']:>7} p={r['p']:>8} -> q16={r['q16']:>9}"
f"|精确 {r['exactQ16']:>14.6f}|相对误差 {r['relErr']:+.9f}")
print(f"P 层对照缓冲 sha256 : {hashlib.sha256(buf).hexdigest()}")
print(f"缓冲值(6 个整数) : {vals}")
return 0
if __name__ == "__main__":
sys.exit(main())
== 7.1 概率密度与权重 ==
归一(半球上估 ∫ p dω,1048576 个均匀样本,采样 pdf = 1/2π):
p = cosine :1.000334(应 1)
p = uniform :0.999954(应 1)
p = unnormalized:3.142802(应 π)
形状(f(x) = 2x,I = 1;256 子流 × 4096 样本):
p = uniform :均值 1.000166|相对标准差 0.009193|最大权重 2.0000|除零 0
p = matched :均值 0.999987|相对标准差 0.000056|最大权重 1.0000|除零 14
p = anti :均值 0.999981|相对标准差 0.042460|最大权重 255.0000|除零 0
定点权重表(fdiv vs 精确 Q16.16 值 f·65536/p):
f= 65536 p= 65536 -> q16= 65536|精确 65536.000000|相对误差 +0.000000000
f= 131072 p= 65536 -> q16= 131072|精确 131072.000000|相对误差 +0.000000000
f= 65536 p= 131072 -> q16= 32768|精确 32768.000000|相对误差 +0.000000000
f= 3 p= 131072 -> q16= 1|精确 1.500000|相对误差 -0.333333333
f= 65536 p= 196608 -> q16= 21845|精确 21845.333333|相对误差 -0.000015259
f= 5 p= 458752 -> q16= 0|精确 0.714286|相对误差 -1.000000000
f= 1 p= 65536 -> q16= 1|精确 1.000000|相对误差 +0.000000000
f= 65536 p=65536000 -> q16= 65|精确 65.536000|相对误差 -0.008178711
f= 17 p= 327680 -> q16= 3|精确 3.400000|相对误差 -0.117647059
f= 65536 p= 7 -> q16=613566756|精确 613566756.571429|相对误差 -0.000000001
P 层对照缓冲 sha256 : dc20a72b1562c0de2740d98ce2d2b91598e03d79de3b48ff96c7b47e0fd6cea6
缓冲值(6 个整数) : [1234936, 1048528, 3879846, 1147366, 1048576, 1032008]
P 层把估计量钉死
P 层用 TypeScript 重写同一条采样流,全程 BigInt,断言对象是那六个整数的 sha256:
// file: lantern/core/test/cg_trace_07_1_parity.test.ts
// tier: P
// verify: node --test 'core/test/cg_trace_07_1_parity.test.ts'
// 7.1 的最小 P 复现:读 Python 侧 cg_trace_07_1_pdf.json,把归一化常数、
// 三种权重形状的极值、以及逐字节对照缓冲的 sha256 钉住。
import assert from "node:assert/strict";
import { createHash } from "node:crypto";
import { readFileSync } from "node:fs";
import { dirname, join } from "node:path";
import { test } from "node:test";
import { fileURLToPath } from "node:url";
import { parityBuffer, pdfWeight } from "../src/lantern/cg_trace_07_1_pdf_weight.ts";
const root = join(dirname(fileURLToPath(import.meta.url)), "..", "..");
const py = JSON.parse(
readFileSync(join(root, "out", "cg_trace_07_1_pdf.json"), "utf8"),
) as {
norm: { cosine: number; uniform: number; unnormalized: number };
shape: Record<string, { relStd: number; maxWeightFloat: number; zeroPdfCount: number }>;
parity: { values: number[]; bufferSha256: string; bufferLen: number };
};
test("【核心】归一化常数:匹配 pdf 的均值是 1.000334", () => {
assert.equal((py.norm.cosine / 65536).toFixed(6), "1.000334");
assert.equal((py.norm.uniform / 65536).toFixed(6), "0.999954");
assert.equal((py.norm.unnormalized / 65536).toFixed(6), "3.142802");
});
test("【核心】匹配 pdf 的权重恒为 1,相对标准差接近 0", () => {
assert.equal(py.shape.matched.maxWeightFloat, 1);
assert.ok(py.shape.matched.relStd < 1e-4, `匹配 pdf 的相对标准差 ${py.shape.matched.relStd}`);
assert.equal(py.shape.matched.zeroPdfCount, 14);
});
test("【核心】逐字节对照缓冲 sha256 与 Python 侧相同", () => {
const { buf, vals } = parityBuffer();
assert.equal(buf.length, py.parity.bufferLen);
assert.deepEqual(vals.map(Number), py.parity.values);
assert.equal(createHash("sha256").update(buf).digest("hex"), py.parity.bufferSha256);
});
test("【负例】反向 pdf 把最大权重放大到 255 倍", () => {
assert.equal(py.shape.anti.maxWeightFloat, 255);
assert.equal(py.shape.uniform.maxWeightFloat, 2);
assert.equal(pdfWeight(65536n, 65536n), 65536n);
assert.equal(pdfWeight(1n, 1000n * 65536n), 0n); // 1/1000 在 Q16.16 下截断为 0
});
P 层缓冲里塞的是同一条流上的整数(半球均匀采样三个累加值 + 三条 pdf 的权重和,各 16 个样本),不是浮点结果——浮点两侧不可能逐位对齐,整数可以。
✔ 【核心】归一化常数:匹配 pdf 的均值是 1.000334 (0.463208ms)
✔ 【核心】匹配 pdf 的权重恒为 1,相对标准差接近 0 (0.105334ms)
✔ 【核心】逐字节对照缓冲 sha256 与 Python 侧相同 (1.102709ms)
✔ 【负例】反向 pdf 把最大权重放大到 255 倍 (0.068959ms)
ℹ tests 4 / pass 4 / fail 0
小结
- 估计量只有一个形状
f/p;p换形状不改变期望,只改变方差。 p的第一条要求是归一。半球上实测:cosθ/π得 1.000334、1/(2π)得 0.999954、漏掉常数的cosθ得 3.142802(π)。比值 3.14175,一次就能判出漏没漏常数。p的第二条要求是形状。I = ∫₀¹ 2x dx上:匹配 pdf2x的相对标准差 5.6×10⁻⁵(权重恒为 1,方差恰为 0),均匀 pdf 9.2×10⁻³,反匹配 pdf2(1−x)4.2×10⁻² 且最大权重 255。- 三条 pdf 的均值都收敛到 1,相对偏差 ≤ 0.02%——形状不影响无偏性,只影响方差。
- 定点上
fdiv向零截断:绝对误差恒 < 1,权重小时相对误差放大(0.714 被截成 0),分母为 0 时抛异常。这两条在 7.3 的 MIS 权重里都会再次出现。
下一节把这条估计量用到真实场景上:同一个直接光照积分,用两条完全不同的采样路径去估。
阅读导航:上一节:6.3 无偏性与偏差 · 下一节:7.2 余弦加权与光源采样 。
继续阅读
探索更多技术文章
浏览归档,发现更多关于系统设计、工具链和工程实践的内容。