o1:从蒙特卡洛到MCMC与光线追踪:采样方法的本质与应用(长篇推导版)
引言:为什么我们需要“聪明”的采样?
在计算机图形学(尤其是光线追踪渲染)中,我们往往希望估计一些对高维空间进行积分的物理公式。然而,用均匀采样来直接做蒙特卡洛估计,往往需要极大的采样量才能保证结果不出现明显噪点。这是由于在现实渲染场景中,能量分布并非均匀——有些方向或区域贡献极大,有些几乎可以忽略。简单均匀撒点显然会浪费大量算力。
为解决这一难题:
- 重要性采样 (Importance Sampling) 通过显式设计“提议分布”\(q(x)\),让我们更频繁地在高贡献区域采样,以此降低估计的方差。
- 马尔可夫链蒙特卡洛 (MCMC) 通常在我们对分布形状没有清晰把握或维度很高时,通过构造马尔可夫链,在目标分布上“漫步”,间接地逼近该分布,从而自动在高贡献区域聚焦更多采样。
本博客将详细阐述这两条路线完整的公式推导,并结合光线追踪的常见案例,说明如何落地应用。请做好长篇数学推导的准备,一起来深入探究采样方法的魅力吧。
第一章:蒙特卡洛方法基础
1.1 蒙特卡洛积分:用随机数计算定积分
(1) 从定积分到期望
我们先来看经典的定积分问题:
对于一维情况,如果 \(f\) 没有解析积分的闭式解,我们可以采用数值方法(如梯形法、辛普森法等)来近似。但一旦问题扩展到高维甚至超高维,传统数值方法的误差会十分明显,甚至无法满足需求。
蒙特卡洛方法的思路是将定积分转换为期望。若我们令 \(X \sim \mathrm{Uniform}(a,b)\),那么
原因在于对于均匀分布\(\mathrm{Unif}[a,b]\),我们有
从而得到
一旦能够将积分写成期望形式,就可以用样本均值来估计期望:
当 \(N \to \infty\) 时,根据大数定律 \(\hat{I} \to I\) 几乎必然收敛;按中心极限定理,其误差大约在 \(O\bigl(\frac{1}{\sqrt{N}}\bigr)\) 的尺度上。
(2) 高维情况与均匀采样的痛点
在高维情形下,如果我们令 \(X \sim \mathrm{Uniform}(\mathcal{D})\) 其中 \(\mathcal{D}\) 是一个高维域,理论上同样可以用
但问题在于:\(f(\mathbf{x})\) 在高维域通常分布极为不均,只有极少数区域贡献巨大,其它大部分区域贡献很小。均匀采样有可能很难“撞到”高贡献区域,导致方差极高、收敛极慢,这即是所谓的“维数灾难 (Curse of Dimensionality)”。
(3) 估计圆面积的经典例子
- 在二维情形,我们可以通过在\([-1,1]\times[-1,1]\)区域均匀撒点,然后数多少点落在单位圆内,再乘以4就得到估计的圆面积。
- 随着采样点\(N\)的增大,圆面积估计会逐渐逼近\(\pi\)。
这虽是蒙特卡洛方法最常见的启蒙示例,却也揭示了它的通性:只要能对目标区域或目标分布采样,就能通过样本均值近似所需的积分或期望。
1.2 蒙特卡洛的方差与收敛
任何蒙特卡洛估计都存在方差的问题。以上一节的一维情形为例:
其期望显然是 \(I\)。若定义 \(\sigma^2 = \mathrm{Var}(f(X))\),则可以得到:
所以我们最直观的手段是:增大 \(N\) 以降低方差;实际上往往需要指数级增长的 \(\!N\) 才能在高维空间中获得较好精度。
这就引出了重要性采样和MCMC:我们无法无限地增加 \(N\),那就必须靠更好地利用对 \(f\) 或其分布的了解,从而在有限样本下有效降低方差。
第二章:重要性采样(Importance Sampling)
2.1 从均匀采样到重要性采样
统一一下符号,考虑我们需要计算的积分形式
假设我们找到一个更“贴近”\(f(x)\)的分布 \(q(x)\) ,只要 \(q(x) >0\) 覆盖了 \(f(x)\) 的支撑集(non-zero 区域),就可以写成:
通过从 \(q(x)\) 独立采样 \(x_1, x_2,\dots, x_N\),我们定义估计量
就能保证 \(\mathbb{E}[\hat{I}_{\mathrm{IS}}] = I\),实现无偏估计。
为什么有用?
与原本的均匀采样相比,如果 \(q(x)\) 在 \(f(x)\) 大值区域取值概率更高,这时对那些区域会采更多样本、而对不重要区域则少采或不采,可以显著降低方差。在计算机图形中,最常见的例子就是:光源很亮,或者BRDF方向比较集中的地方,其贡献显著,适当“重点采样”可将噪点大幅减小。
2.2 重要性采样估计的方差及最优分布
要理解重要性采样的好处,让我们看看其方差。定义:
方差为
而
根据无偏性,有
因此
若要使 \(\mathrm{Var}[w(x)]\) 尽可能小,我们需要最小化
通过拉格朗日乘子法可证明,当
方差可达到最小值;并且当 \(q^*(x) = \frac{|f(x)|}{\int |f(x)|\,dx}\) 时,\(\mathrm{Var}\bigl[w(x)\bigr]\) 达到理论下界。这也就说明了最理想的分布就是与 \(f\) 自身正比。
2.3 实际实现中的难点
- 我们往往不知道 $ f(x) $ 的确切分布形态,或者算出 $ \int |f(x)| ,dx $ 并不容易;
- 即使知道了,可能也难以从这种分布 \(q^*(x) \propto |f(x)|\) 实际采到样本;
- 在高维场景下,构造高维分布和高效采样算法本身就是不小的挑战。
这些难点正是我们在具体应用中选择近似或启发式重要性分布、甚至选择放弃显式重要性分布,转用MCMC来做高维探索的原因。
2.4 光线追踪中的例子:对光源进行重要性采样
例如,一个场景下有一个很亮的小光源,若只进行均匀采样,可能多数方向都无法“看到”光源,导致大量样本贡献几乎为零。通过对光源方向进行采样,就能在每次弹射光线时都大概率命中光源,从而减少噪点,提升渲染质量和效率。
第三章:马尔可夫链蒙特卡洛(MCMC)
3.1 为什么需要马尔可夫链?
如果我们能够给出一个比较好的、贴近\(f(x)\)的分布 \(q(x)\) 并能高效采样,那么重要性采样其实是非常好的方法。然而,在大多数高维概率推断或复杂光传输场景中,我们往往无法直接设计出合适的 \(q(x)\)。此时就轮到马尔可夫链蒙特卡洛 (MCMC)。
MCMC 的基本思路是:构造一个马尔可夫链 \(\{x_t\}_{t=1}^\infty\),让它在目标分布\(\pi(x)\)上“游走”,当它“收敛/混合”后,链中的样本就能近似视作来自 \(\pi\) 的采样,从而可做蒙特卡洛估计:
其中 \(B\) 是“燃烧期”(Burn-in)的长度,表示需要丢弃前面一段还没有收敛的样本。
3.2 马尔可夫链与平稳分布
一个马尔可夫链 \(\{x_t\}\) 由转移核 \(T(x'|x)\) 定义:当前状态是 \(x\) 时,下一个状态 \(x'\) 的分布即 \(T(\cdot|x)\)。若这个转移核满足一些条件(如不可约 irreducible、非周期 aperiodic、正再生 positive recurrent 等),就存在且只存在一个平稳分布 \(\pi(x)\) 使得:
当马尔可夫链运行到平稳后,分布就不会再变,随后的样本就可以视为采样自 \(\pi\)。MCMC恰恰是要精心设计这样的 \(\pi\) 和 \(T\),使之在满足理论收敛条件的情况下,最终分布即是我们想要逼近的目标分布。
3.3 Metropolis–Hastings算法:细致平衡的关键
MCMC 的实现方法很多,其中Metropolis–Hastings (MH) 算法最常用。它的核心在于细致平衡 (Detailed Balance) 约束:
如果满足这个约束,就能保证 \(\pi\) 是一个平稳分布。MH算法如何强制这个约束呢?过程大致如下:
- 提议分布 \(q(x'|x)\):先随意选一个我们能采样的分布,比如某个对称分布 \(q(x'|x) = q(x|x')\) 或者简单的高斯跳步。
- 接受率:从 \(q(x'|x_t)\) 采到一个候选 \(x'\) 后,不一定全部接受,而是用一个概率 \(\alpha(x'|x_t)\) 来决定是否“跳”过去。
- Metropolis–Hastings 给出的接受率为:\[\alpha(x'|x_t) = \min\!\Bigl(1,\;\frac{\pi(x')\,q(x_t|x')}{\pi(x_t)\,q(x'|x_t)}\Bigr), \]其中 \(\pi(x)\) 是(未归一化的)目标分布。
- Metropolis–Hastings 给出的接受率为:
- 细致平衡证明:
- 若接受,则 \(x_{t+1}=x'\);否则 \(x_{t+1}=x_t\)。
- 可以证明这样设计的转移核 \(T\) 就会满足细致平衡,从而\(\pi\) 成为链的平稳分布。
(1) MH算法的具体步骤
- 选定一个初始状态 \(x_0\)。
- 第 \(t\) 步时:
- 从提议分布 \(q(\cdot|x_t)\) 采一个候选样本 \(x'\)。
- 计算接受率\[\alpha = \min\!\Bigl(1,\;\frac{\pi(x')\,q(x_t|x')}{\pi(x_t)\,q(x'|x_t)}\Bigr). \]
- 生成一个 \(\mathrm{Uniform}(0,1)\) 数 \(u\)。如果 \(u<\alpha\),则接受 \(x'\);否则拒绝并保持 \(x_{t+1} = x_t\)。
- 当 \(t\) 足够大时,样本 \(\{x_t\}\) “大致”分布于 \(\pi\)。可丢弃前面一些“混合期”的样本,再用后续样本做蒙特卡洛估计。
(2) 细致平衡推导要点
我们想要
其中
如果我们设
便能保证细致平衡成立。主要推理步骤如下:
- 当 \(\frac{\pi(x')\,q(x|x')}{\pi(x)\,q(x'|x)} \le 1\) 时,接受率 \(\alpha(x'|x) = \frac{\pi(x')\,q(x|x')}{\pi(x)\,q(x'|x)}\)。可带入验证平衡条件。
- 当 \(\frac{\pi(x')\,q(x|x')}{\pi(x)\,q(x'|x)} > 1\) 时,\(\alpha=1\),也同样能满足细致平衡。
从而得出经过足够长时间后,马尔可夫链的平稳分布正是 \(\pi\)。
3.4 诊断与调参:Burn-in、自相关、步长
- Burn-in:一般在链开始时会有一段不稳定阶段(尚未收敛到目标分布),我们需要丢弃此段样本。
- 自相关:链的样本并非独立,尤其提议分布\(q\)选得不好时,这些样本的相关性较高,导致估计效率下降。可以借助自相关时间(autocorrelation time)等指标来判断,必要时只保留每隔几步的样本。
- 步长 (proposal scale):若选择过小,链移动缓慢,相关性高;选择过大,又容易导致接受率过低,频繁拒绝,也降低效率。MH算法往往需要反复调参或使用自适应MCMC技术来改善。
第四章:重要性采样 vs MCMC:在什么场景该选用?
4.1 比较二者的优缺点
-
重要性采样
- 优点:样本是独立的;如果能很好地匹配 \(f(x)\) 的峰值区域,方差通常很小。
- 缺点:需要事先设计可采样的 \(q(x)\);在高维或复杂分布下,往往难以满足。
-
MCMC
- 优点:不需要显式知道归一化常数,只需能计算 \(\pi(x')/\pi(x)\) 即可;适合极其复杂或高维度的分布。
- 缺点:样本具相关性;需要设置提议分布、Burn-in、步长等等;链收敛速度也可能很慢。
在渲染中,如果能明确光源、材质峰值等重要区域,重要性采样往往更直接有效;若场景中有非常复杂的光路,如焦散等,或分布带有极稀疏但很高贡献的部分,MCMC(尤其是 Metropolis 光传输)就成为可行手段。
4.2 融合应用举例
- 粒子滤波:通过重要性采样更新粒子权重,然后再进行重采样,从而结合了“基于观测更新分布”与“去除低权重粒子”的优点,本质上是一个顺序化的重要性采样过程。
- 自适应MCMC:在马尔可夫链迭代过程中,根据历史样本适时更新提议分布 \(q\),让它逐渐逼近最优分布,从而降低自相关、提升接受率,这种融合了重要性采样思路与马尔可夫链迭代自适应方式。
第五章:光线追踪中的蒙特卡洛革命
5.1 路径追踪:蒙特卡洛积分的经典落地
所谓路径追踪 (Path Tracing),就是通过蒙特卡洛方法来解渲染方程
其中 \(\mathbf{x}\) 是表面点,\(\boldsymbol{\omega}_o\)和\(\boldsymbol{\omega}_i\) 分别是出射和入射方向,\(f_r\) 是 BRDF,\(L_e\) 表示自发光部分,\(\Omega\) 是半球积分域。
(1) 基本蒙特卡洛表达
如果我们在半球 \(\Omega\) 上做均匀采样(或余弦加权采样),可以得到一个近似公式:
这里 \(\boldsymbol{\omega}_i^*\) 为采样方向。
然而,如果场景中灯光数量多,或者某些方向贡献突出,简单的均匀采样可能噪点很大。
(2) 重要性采样在光线追踪中的应用
- 基于 BRDF 的重要性采样:若材质具有高光(glossy reflection),反射分布 \(f_r\) 可能尖锐集中,这种情况下按照 BRDF 分布来采样入射方向效果更佳;
- 基于光源的采样:如果光源有位置或形状优势,针对光源表面重点采样,可以急速减少阴影或强光直射处的噪点。
(3) 多重重要性采样 (Multiple Importance Sampling)
有时我们既想对 BRDF 做重要性采样,又想对光源做重要性采样。可以考虑将二者结合成“多重重要性采样 (MIS)”,并且用任意多种扩展(如 Power Heuristic)来组合权重,确保得到无偏且方差更低的结果。
5.2 MCMC渲染:Metropolis Light Transport (MLT)
当场景存在焦散或者很难通过一次次随机反弹就恰好打到光源或某些特定路径时,我们可以利用MCMC来进行路径采样。典型方案是Metropolis Light Transport (MLT):
- 把一条完整的光路(由顶点及其连接关系组成)视作马尔可夫链中的一个状态;
- 在当前路径的基础上,随机做局部突变或全局突变(如随机修改某个反射点、或替换整段子路径);
- 通过 Metropolis–Hastings 接受率,决定是否更新到新路径;
- 当有某次采样驶入一个特别“高能量”区域(比如镜面聚焦+折射焦散等),由于局部随机扰动会保持在这附近,能显著减少噪点。
MLT 在复杂的光传输场景中(尤其是有大量焦散现象)相对于简单路径追踪会更有效地集中采样到关键的路径,从而在相同样本数下得到更佳的图像质量。
5.3 代码示例(简化版)
以下为一个极简化的路径追踪示例,突出蒙特卡洛递归思想;实际生产渲染器要复杂得多:
Color pathTraceRay(const Ray& ray, Scene& scene, int depth) {
if (depth > MAX_DEPTH) return Color(0,0,0);
Intersection hit = scene.intersect(ray);
if (!hit.hitObject) {
// 没撞到场景对象,则返回天空或背景
return scene.getBackgroundColor();
}
// 1. 存在自发光
Color Le = hit.material.emission;
// 2. 重要性采样:根据 BRDF 获得新方向
// or 根据光源分布来采样
Vector3 newDir = sampleDirectionByBRDF(hit.normal, hit.material);
float pdf = brdfPDF(hit.normal, newDir, hit.material);
// 3. 计算间接光: 递归
Ray newRay(hit.position, newDir);
Color L_indirect = pathTraceRay(newRay, scene, depth+1);
// 4. 蒙特卡洛估计结合BRDF和pdf
Color brdfVal = evaluateBRDF(hit.normal, newDir, ray.dir, hit.material);
// cosTerm = dot(newDir, hit.normal)
float cosTerm = std::max(0.f, dot(newDir, hit.normal));
Color result = Le + brdfVal * L_indirect * cosTerm / pdf;
return result;
}
这里的 sampleDirectionByBRDF 就是一个重要性采样的核心:它让采样更倾向于材质反射高峰区域,而不是在半球上随便乱挑方向。
第六章:总结与进阶
6.1 核心结论
- 蒙特卡洛方法的本质:利用随机采样近似积分或期望,其精度随样本数\(\sqrt{N}\)比例收敛。
- 重要性采样:如果能有效设计贴近目标函数的分布,就可显著降低方差;但在高维或不熟悉目标分布时困难较大。
- 马尔可夫链蒙特卡洛:通过在目标分布上构造马尔可夫链实现采样,几乎只需能计算相对比例 \(\pi(x')/\pi(x)\),不需显式归一化常数,因而广泛用于复杂概率推断、渲染中的焦散等场景。
6.2 延伸阅读
- 哈密尔顿蒙特卡洛 (HMC):使用目标函数的梯度信息,让马尔可夫链更“有方向”地移动,在统计推断中常见。
- 随机梯度MCMC:结合大数据场景、随机梯度下降与 MCMC 的想法,用于处理巨量数据的贝叶斯学习问题。
- Metropolis Light Transport 原始论文:Eric Veach & Leonidas J. Guibas, Metropolis Light Transport, SIGGRAPH 1997。
- 多重重要性采样 (MIS):Eric Veach, Robust Monte Carlo Methods for Light Transport Simulation, PhD thesis (Stanford University, 1997)。
6.3 思考题
- 若你已经通过 MCMC 得到一批相关样本,是否还能结合重要性采样做二次修正?思考由 MCMC 后验分布和一些先验启发式的逆向结合。
- 在实时光线追踪中(如游戏引擎),如何将重要性采样与低差异序列 (low-discrepancy sequences) 结合以减少噪点,又能维持较快帧率?
- 当目标分布是高斯混合这种“多峰”结构时,传统 MH 算法如果初始点落在一个峰值附近,可能很难跳到其他峰值区域,如何对提议分布或算法进行改进?
附录:更完整的数学推导
A.1 从大数定律到中心极限定理
- 大数定律 (Law of Large Numbers):若 \(\{X_i\}\) 是同分布且独立的随机变量,\(\mathbb{E}[X_i]=\mu\) 有限,则\[\frac{1}{N}\sum_{i=1}^N X_i \xrightarrow[]{a.s.} \mu,\quad N\to\infty. \]这保证了蒙特卡洛估计的无偏和收敛。
- 中心极限定理 (Central Limit Theorem):在同样条件下,\[\sqrt{N}\Bigl(\frac{1}{N}\sum_{i=1}^N X_i - \mu\Bigr) \xrightarrow{d} \mathcal{N}\bigl(0,\sigma^2\bigr), \]这给出了估计量方差收敛随 \(\frac{1}{\sqrt{N}}\) 的速度。
A.2 重要性采样方差的推导
让 \(X\sim q(x)\),定义权重 \(w(X) = f(X)/q(X)\)。
其中
并且
故
最优 \(q^*(x)\) 使得 \(\mathrm{Var}[w(X)]\) 最小,满足 \(q^*(x) \propto |f(x)|\)。
A.3 Metropolis–Hastings细致平衡证明
令
我们需要
即
当我们取
可以分别讨论“比率 \(\le 1\)” 和“比率 \(\ge 1\)” 两种情况,一一验证能满足详细平衡,故 \(\pi\) 为平稳分布。
A.4 MCMC收敛依据 (Perron-Frobenius定理简述)
在可数或连续状态空间上,如果马尔可夫链满足不可约、非周期、\(\pi\)-遍历等条件,便存在唯一的平稳分布 \(\pi\),并且从任意初始分布出发都可以收敛到 \(\pi\)。这类定理为MCMC 的合理性提供了坚实理论依据。
过渡与Callback设计
- 章节过渡:我们先看了“蒙特卡洛”的通用方法,发现高维/不均匀时方差暴涨,引到“重要性采样”能大量减小方差。然而,当无法显式设计分布时,再转向“MCMC”通过马尔可夫链的方式实现目标分布的采样。
- 在光线追踪中呼应:蒙特卡洛路径追踪不可或缺,但要更高效地对焦散等场景采样,可借助 MCMC(MLT)。
- 方差与计算成本的综合权衡:所有采样方法的目标,都围绕“如何又快又好地逼近积分”这一点展开。
结语
至此,我们已经从最朴素的蒙特卡洛积分讲到重要性采样,再到马尔可夫链蒙特卡洛 (MCMC),并在光线追踪中看到这些采样技术的集大成应用:
- 当可以显式构造贴近目标函数的分布时,重要性采样能很直接地降低方差。
- 如若目标分布极其复杂或不足以显示设计提议分布,MCMC 能够通过构造马尔可夫链逐步逼近目标分布。
- 在实际图形渲染中,两者常可结合,如多重重要性采样 (MIS) 或 Metropolis Light Transport (MLT) 等,用于处理不同类型的光路和材质特性。
希望这篇较为详尽的博客能让你对采样方法在高维积分、统计推断及光线追踪渲染中的原理与细节有更深刻认识。如果你对更深层次的理论(如HMC、随机梯度MCMC等)或更多渲染细节感兴趣,欢迎继续拓展阅读或留言讨论。祝学习愉快!

浙公网安备 33010602011771号