3D数学系列之——再谈特卡洛积分和重要性采样

一、前篇文章回顾

  在前一篇文章3D数学系列之——从“蒙的挺准”到“蒙的真准”解密蒙特卡洛积分！中，介绍了下面的这些朴素的公式：
$$\begin{array}{rl}& \frac{灰色部分中点数}{总的点数}\approx \frac{灰色部分面积}{矩形面积}\\ & \Rightarrow 灰色部分面积\approx 矩形面积\times \frac{灰色部分中点数}{总的点数}\end{array}$$

$$灰色部分面积\approx 矩形面积\times P（P为点落在灰色部分的概率）$$

  虽然最终通过一个简单的概率乘法就可以计算出积分值，但是如何得到这个精确的概率值从而得到比较精确的积分值却不是那么简单的。并且这个方法直接用于计算或者说模拟计算来说，不确定的东西太多。下面我们就继续从积分的基本公式开始，讨论一下看看有没有进一步更形式化和更精确一点的方法。

二、积分的黎曼和形式

  对于数值积分计算来说，积分的黎曼和形式是最基本的形式：
$$\underset{a}{\overset{b}{\int }}f(x)\mathrm{d}x\approx \sum _{n=1}^{N}f(x_n)\frac{b-a}{N}\text{\hspace{1em}(1)}$$
  这是一个常见和基本的可以用于积分计算程序的基本公式。

三、积分的概率形式（蒙特卡洛积分）

  接着我们对黎曼和使用著名的”陶哲轩瞪眼法“进行观察。

  首先对于求和来说，第一项是 $f(x_n)$ ，这里我们一般会比较容易知道 $f(x)$ 的解析形式或者可以知道它的一些离散值。

  在实际计算中需要生成的是 $x_n$ 的序列，这个序列是跟后面的 $\frac{(b-a)}{N}$ 紧密关联的，因为需要保证 $x_i-x_{i-1}=\frac{b-a}{N}$ 即 $x_n$ 序列是个等差数列，间隔是就是 $\frac{(b-a)}{N}$ 。

  这时，我们先打开思路，既然前一讲中已经可以通过随机变量的形式来计算积分，那么这时我们可不可以用一组随机数的序列来代替 $x_n$ 的序列呢？

  既然讲到随机，那么回忆一下概率论中关于“均匀分布随机变量”的概率密度函数（The Probability Distribution Function，缩写为 PDF）：
$$\mathrm{p}(x)=\{\begin{array}{l}\frac{1}{b-a},a<x<b\\ 0,其他\end{array}\text{\hspace{1em}(2)}$$
  接着我们用"陶哲轩瞪眼法"观察黎曼和公式与上面这个公式，此时我们应灵光乍现想到下面这样的形式：
$$\begin{gathered} \underset{a}{\overset{b}{\int }}f(x)\mathrm{d}x\approx \sum _{n=1}^{N}f(x_n)\frac{b-a}{N} \\ =\frac{1}{N}\sum _{n=1}^{N}f(x_n)\times (b-a) \\ 根据乘以一个数等于除以一个数的倒数这个法则，有： \\ =\frac{1}{N}\sum _{n=1}^N\frac{f(x_n)}{\mathrm{p}(x_n)}\text{\hspace{1em}(3)} \end{gathered}$$
  啊哈！太好了，这一切看似是公式间的一种巧合或者说技巧性的，接着我们思考下它的意义。因为刚才我们已经说过原来的$x_n$ 序列是个等差序列，本身就是均匀分布的，所以其概率密度函数（PDF）就是均匀分布变量的 PDF。如果我们取一个随机的 $x_n$ 序列，并假设知道它的 PDF $\mathrm{p}(x_n)$ 那么就可以根据上面这个式子的结果来计算积分值了。并且当 $x_n$ 不断增多时，和 $\frac{1}{N}\sum _{n=1}^N\frac{f(x_n)}{\mathrm{p}(x_n)}$ 也就越接近原始的积分值！这也就是传说中的蒙特卡洛积分的基本形式。它是合理且有意义的，起码我们不用只靠前一讲中的“蒙的真准”的方法了。需要注意的就是其中 $\frac{f(x_n)}{p(x_n)}$ 项看上去虽然很奇怪，但他其实是有实际含义的，一般被称作概率平均。

  当然使用随机序列 $x_n$ 是有条件的，即要求：
$$\underset{-\infty }{\overset{+\infty }{\int }}\mathrm{p}(x)=1\text{\hspace{1em}(4)}$$
  接着我们继续观察（还是使用 “陶哲轩瞪眼法”）那个带有 PDF 的蒙特卡洛积分公式，发现其中还有个 $\frac{1}{N}$ , 这时我们应该突然明白其实这个式子还是个平均值的公式，而讲到平均值，学过概率论的同学就应该立刻想到“数学期望”。一般的数学期望的定义：
$$E[f(X)]=\sum _{n=1}^{\infty }f(x_n)\mathrm{p}(x_n)$$
  也就是说如果我们知道了随机变量 $x_n$ 的随机 PDF： $\mathrm{p}(x_n)$，那么随机变量函数的数学期望就可以用上式计算。当然如果 $x_n$ 是连续的随机变量，那么就可以用积分形式来表示上面的表达式：
$$E[f(X)]=\underset{-\infty }{\overset{+\infty }{\int }}f(x)\mathrm{p}(x)dx$$
  这时，为了验证蒙特卡洛积分公式的正确性，我们来看看它的数学期望：
$$\begin{gathered} E(F(X))=E\left[\frac{1}{N}\sum _{n=1}^N\frac{f(x_n)}{\mathrm{p}(x_n)}\right] \\ =\frac{1}{N}\sum _{n=1}^N\underset{a}{\overset{b}{\int }}\frac{f(x)}{\mathrm{p}(x)}\times \mathrm{p}(x)\mathrm{d}x \\ =\frac{1}{N}\sum _{n=1}^N\underset{a}{\overset{b}{\int }}f(x)\mathrm{d}x \\ =\underset{a}{\overset{b}{\int }}f(x)\mathrm{d}x \end{gathered}$$
  这个结果简直太棒了！它说明蒙特卡洛积分公式的数学期望就是我们要求的积分值。这也间接的证明我们前面靠瞪眼法得到的公式本质上是正确的。或者说，蒙特卡洛积分的结果对于积分真值是无偏差的，也就是说是原积分的无偏估计。这个可以由“大数定理”来保证。只是不同的随机变量序列 $x_n$ 会以不同的“速度”靠近积分原值。

四、误差

  按照数学家们的“尿性”，一定会严格的分析一下整个“蒙的挺准”的过程中的误差。既然刚才已经分析了蒙特卡洛积分的数学期望，那么继续根据概率论的知识，我们来分析下它的误差。一般对于随机序列来说，我们取其方差来评估其靠近真实值的程度。那么就让我们来计算下蒙特卡洛积分的方差：
$$\begin{gathered} {\sigma }^2[F_n(X)] \\ ={\sigma }^2\left[\frac{1}{N}\sum _{n=1}^N\frac{f(x_n)}{\mathrm{p}(x_n)}\right] \\ =\frac{1}{N^2}\sum _{n=1}^N\underset{a}{\overset{b}{\int }}{\left(\frac{f(x)}{\mathrm{p}(x)}-E(F_n(X))\right)}^2\times \mathrm{p}(x)\mathrm{d}x \\ =\frac{1}{N}\left[\underset{a}{\overset{b}{\int }}{\left(\frac{f(x)}{\mathrm{p}(x)}\right)}^2\times \mathrm{p}(x)\mathrm{d}x-E(F_x(X){)}^2\right] \\ =\frac{1}{N}\left[\underset{a}{\overset{b}{\int }}\frac{f(x{)}^2}{\mathrm{p}(x)}\mathrm{d}x-E(F_n(X){)}^2\right]\text{\hspace{1em}(5)} \end{gathered}$$
  上面最终包含了 $\frac{1}{N}$ 项的最终的一堆公式，就是蒙特卡洛积分的方差，根据概率论的定义，其标准差就是方差的算术平方根。

  按照最终的证明结果，蒙特卡洛积分值与积分真实值之间的标准差（误差）有如下的形式：

$$误差\propto \frac{1}{\sqrt{N}}\text{\hspace{1em}(6)}$$
  即最终误差正比于N的平方根的倒数。

五、蒙特卡洛积分计算与收敛速度

  有了前面这些数学知识的加持，那么我们就可以来看看用蒙特卡洛积分来计算积分时到底“速度”如何？

  其实看到这里，大家应该有个疑惑，我们费劲巴拉的推导了半天，到底在图个啥？

  那么如果你看了前一篇文章（3D数学系列之——从“蒙的挺准”到“蒙的真准”解密蒙特卡洛积分！），首先应该想到有了蒙特卡洛积分公式，我们是不是可以开始靠“蒙”来计算积分了？答案是肯定的，要计算蒙特卡洛积分，我们就需要产生一个随机变量的序列 $x_n$ 并且需要知道其概率密度函数（PDF）$\mathrm{p}(x_n)$ （有时候不需要知道 PDF 的表达式，只需要知道每个随机变量对应的概率值即可） 然后代入公式的右侧进行计算：
$$\underset{a}{\overset{b}{\int }}f(x)\mathrm{d}x\approx \frac{1}{N}\sum _{n=1}^N\frac{f(x_n)}{\mathrm{p}(x_n)}\text{\hspace{1em}(3)}$$
这样计算的好处是，我们可以用大一些（计算能力充足时）或者小一些（计算能力不足时）的随机变量序列来进行计算。

  当然这样计算的限制是，根据对蒙特卡洛方法标准差的计算，最终可以知道如果我们想提高计算精度一倍，那么差不多就需要将N提高4倍（ $\frac{1}{\sqrt{4}}=\frac{1}{2}$ ，即误差减小一半）。所以这就是平常所说的蒙特卡洛积分法收敛慢的根本原因。

  例如在我们关注的图形的光照积分估算过程中，这是个麻烦的问题，比如原先需要每个点采样1000条光线近似计算最终光照效果，要提高一倍精确度从而改进画质效果的话，那么每个点就需要采集4000条光线进行积分。因此，基本上现实中很少直接用蒙特卡洛积分的基本形式来求积分值的。

  所以最终这个结论是有点令人灰心丧气的，不过请继续往后看，这么好的方法一定已经有很多人帮我们想了很多办法来解决问题的。

六、重要性采样

  OK，看到这里你还没有晕过去的话，请为自己点个赞先！

  鉴于前面说的直接运用蒙特卡洛方法进行积分计算时，因为固有偏差的问题，而导致算法过程收敛过慢（N的平方根倒数的量级）所以在计算效率上是不尽人意的。那么既然是因为偏差导致的问题，我们就继续从偏差来分析看看：
$$\begin{gathered} {\sigma }^2[F_n(X)] \\ =\frac{1}{N}\left[\underset{a}{\overset{b}{\int }}\frac{f(x{)}^2}{\mathrm{p}(x)}\mathrm{d}x-E(F_n(X){)}^2\right] \\ =\frac{1}{N}\left[\underset{a}{\overset{b}{\int }}\frac{f(x{)}^2}{\mathrm{p}(x)}\mathrm{d}x-{\left(\underset{a}{\overset{b}{\int }}f(x)\mathrm{d}x\right)}^2\right] \end{gathered}$$
  这时如果我们想要标准差（方差的开方）最小，那么就需要让上式结果等于0，此时我们继续使用瞪眼法可以发现，这需要下面的等式成立：
$$\mathrm{p}(x)=\frac{|f(x)|}{\underset{a}{\overset{b}{\int }}f(x)\mathrm{d}x}\text{\hspace{1em}(7)}$$
  在很多其他的资料中，推得这个式子之后，忽略了讲解它的意义。乍看上去，貌似我们要知道我们需要使用的随机变量的概率密度，就需要先知道最终积分的结果，使得这个式子变成了无用的“鸡肋”。

  其实不然，这个式子恰恰给我们指明了一个非常明确的方向。

  它说明，如果我们要使得方差、标准差最小，那么就要找到一个随机变量序列使得它的概率密度函数与被积函数 $f(x)$ 在“形状”上要高度一致！这是个很重要的信息。为什么这样说呢？其实仔细观察这个式子可以发现其中的定积分其实是个常数 $\underset{a}{\overset{b}{\int }}f(x)\mathrm{d}x$ ，所以整个式子表达的意思是说 ：
$$\mathrm{p}(x)=\frac{|f(x)|}{S}=A|f(x)|\text{\hspace{1em}(8)}$$
  即最终当概率密度函数其实就是被积函数的一个常数倍的时候，整个方差最小，也即标准差会非常小。

  这是整个重要性采样的出发点，也就是说我们可以依据被积函数的特征，构造一个与之匹配的概率分布函数的形状，然后生成特意构造的随机序列，再返回到蒙特卡洛积分 $\underset{a}{\overset{b}{\int }}f(x)\mathrm{d}x\approx \frac{1}{N}\sum _{n=1}^N\frac{f(x_n)}{\mathrm{p}(x_n)}$ 进行计算，那么只需要非常少的 $x_n$ 就可以得到非常精确的积分结果，因为这样的序列计算后使得蒙特卡洛积分最终的方差、标准差都变得非常小，甚至为0。

  也就是说这样的随机序列，使得蒙特卡洛积分能够快速收敛！也就是说，咱们之前说的普通的蒙特卡洛积分过程误差与 $\frac{1}{\sqrt{N}}$ 成正比导致我们必须付出平方倍的计算量，才能使的误差变小的规律可以被打破了，因为这个因子的后面那堆复杂的表达式在我们特意的设计下几乎 = 0，所以我们完全不用浪费过多的计算量来获取较精确的积分结果。 这个过程就被称为重要性采样！

七、重要性采样方法和过程

  那么最终如何进行上面所说的重要性采样呢？其实核心就是构造使的蒙特卡洛积分方差和标准差最小的随机数序列 $x_n$ 。

  1、根据定义 $\underset{a}{\overset{b}{\int }}\mathrm{p}(x)\mathrm{d}x$ 就是概率密度函数 $p(x)$ 对应的累积分布函数（Cumulative Distribution Function），简称做 CDF，通常其正式定义如下：
$$\mathrm{CDF}(x)=\underset{-\infty }{\overset{x}{\int }}p(\tilde{x})\mathrm{d}\tilde{x}\text{\hspace{1em}(9)}$$
  然后选择一个跟被积分的函数“相类似”的概率密度函数 PDF 来计算其 CDF；

  2、接着推导出 CDF 的反函数 ${\mathrm{CDF}}^{-1}(x)$ ；

  3、然后生成一个均匀分布的随机变量序列 $u_n\in [0,1]$ ;

  4、使用 $u_n$ 代入 CDF 的反函数 ${\mathrm{CDF}}^{-1}(x)$ 计算出随机序列 $x_n$ ;

  5、最后将 $x_n$ 代入蒙特卡洛积分中计算函数的积分值；

  通常在很小规模的 $x_n$ 序列下，就可以计算出很高精度的积分值。

  举例来说，假设我们需要求一个形似 $\underset{0}{\overset{\pi }{\int }}A\sin (x)\mathrm{d}x$ 的积分，此时根据被积函数为正弦函数的特征，构造 $p(x)=c\sin (x)x\in [0,\pi ]$ ，然后我们计算它的 CDF：
$$\begin{gathered} \underset{0}{\overset{\pi }{\int }}c\sin (x)\mathrm{d}x=c[-\cos (x){]}_0^{\pi }=2c=1 \\ \because {\int }_0^{\pi }p(x)\mathrm{d}x=1 \\ \Rightarrow c=\frac{1}{2} \\ 用x替换积分上限，得到CDF的表达式，有： \\ CDF(x)=\underset{0}{\overset{x}{\int }}\frac{1}{2}\sin (t)\mathrm{d}t=\frac{1}{2}[-\cos (t){]}_0^x=\frac{1}{2}(1-\cos (x)) \end{gathered}$$
  有了 CDF 那么我们来计算其反函数：
$$\begin{gathered} y=\frac{1}{2}(1-\cos (x)) \\ \Rightarrow x=\arccos (1-2y) \end{gathered}$$
  接着我们生成一个均匀分布的随机数序列 $u_n{u}_n\sim [0,1]$ 代入上式中计算得到指定分布的随机序列 $x_n=\arccos (1-2u_n)$，最终用这个序列按照蒙特卡洛积分公式计算即可快速得到 $\underset{0}{\overset{\pi }{\int }}A\sin (x)\mathrm{d}x$ 的积分值。

  那么详细的计算过程和程序就先略过了，在后续的 IBL 教程中将有更详细的应用讲解和实际工程中的代码。有兴趣的同学可以自己编写程序比较一下，即先用黎曼和及库函数sin计算得到积分结果，然后与库函数cos的直接计算结果进行比较，最好能求出方差，然后再用刚才推导出的 $CD{F}^{-1}$ 函数计算的序列代入蒙特卡洛积分公式中计算积分结果，同样与cos函数直接计算的结果求差，或者求方差，看看有什么区别？

八、重要性采样的优缺点

  根据前面的叙述，重要性采样的全部目的就是为了缩小 $x_n$ 的规模，同时还能保证一定的精度，其实就是之前标题所说的从 “蒙的挺准” 到 “蒙的更准”，而且计算量还大大变小了，这对于很多需要数值积分的计算来说是非常非常好的优化方法。另外对于高维的积分来说，重要性采样的方法依然有效！这样带来的好处就是，原本复杂度可能是 $O(n^{\alpha })\alpha \in Z^{+},\alpha ⩾1$ 的问题，有可能直接变成了 $O(n)$ 的问题，这对于我们关注的图形光照计算领域来说简直就是福音！因为我们需要处理的图形往往都是 $n^3$ 维度的，并且在复杂光照渲染算法中需要大量的 3D 积分，所以有了蒙特卡洛积分和重要性采样方法，就使得该过程的计算量几何级数的被减少，这最终使得实时 PBR 渲染成为可能。

  当然重要性采样也有缺点，根据刚才描述的过程，大家应该能想到，第一，如果被积函数 $f(x)$ 过于复杂时，我们基本无法简单的模拟出一个形似的 $p(x)$ 函数，从而整个过程有可能就前功尽弃；第二个方面就是说，通过 $CD{F}^{-1}$ 函数计算出的序列 $x_n$ 有可能使得 $f(x)$ 的值过于集中于一些较大函数值得地方，从而导致最终积分的结果与真值之间产生正偏差。

  当然幸运的是这两条缺点在我们的 PBR 渲染中几乎碰不到，也就没什么太大的影响。