关于蒙特卡洛重要性采样无偏性、有效性的条件

---

前言

前几周看到一个面试问题，给哥们整玉玉了。

蒙特卡洛采样中的最优 PDF ，它满足什么条件时采样具有无偏性，什么条件下方差最小？

今天整理了一些相关的数学证明。

蒙特卡洛积分与重要性采样

求$f(x)$的期望，其中 x 服从 p(x)分布：

$$E_{x\sim p(x)}(f(x))=\int f(x)p(x)dx$$

分子分母同乘$q(x)$得到：

$$E_{x\sim p(x)}(f(x))=\int f(x)\frac{p(x)}{q(x)}q(x)dx$$

实际上是求$f(x)\frac{p(x)}{q(x)}$的期望，其中 x 服从$q(x)$分布：

$$E_{x\sim p(x)}(f(x))=E_{x\sim q(x)}(f(x)\frac{p(x)}{q(x)})$$

离散化

$$E_{x\sim p(x)}(f(x))=\frac{1}{N}\sum _{x_i\sim p(x)}f(x_i)=\frac{1}{N}\sum _{x_i\sim q(x)}f(x_i)\frac{p(x_i)}{q(x_i)}$$

蒙特卡洛重要性采样无偏性条件

根据$E_{x\sim p(x)}(f(x))=\frac{1}{N}\sum _{x_i\sim q(x)}f(x_i)\frac{p(x_i)}{q(x_i)}$，我们对$p(x)$和$q(x)$进行分类谈论：

Situation	$p(x)>0$	$p(x)=0$
$q(x)>0$	应该被计算，实际上也确实被计算了	对最终答案没有贡献，实际上也被计算了，但是无伤大雅	实际上能被采样到
$q(x)=0$	应该被计算，但是因为$q(x)=0$，我们采样不到这样的点（而且$q(x)=0$作分母无意义）	对最终答案没有贡献，实际上没有被计算（正确的）	实际上采样不到
	对最终结果有贡献，应该被计算	无效计算

更直观的，我们可以通过交并图来理解：

。。。。。。

蒙特卡洛重要性采样有效性条件

方差等于平方的期望减去期望的平方：

$$D_{x\sim q(x)}(f(x)\frac{p(x)}{q(x)})=E_{x\sim q(x)}(f(x{)}^2\frac{p(x{)}^2}{q(x{)}^2})-E_{x\sim q(x)}(f(x)\frac{p(x)}{q(x)})\cdot E_{x\sim q(x)}(f(x)\frac{p(x)}{q(x)})$$

$q(x)$分布下的期望可以转换为$p(x)$分布下的期望：

$$D_{x\sim q(x)}(f(x)\frac{p(x)}{q(x)})=E_{x\sim q(x)}(f(x{)}^2\frac{p(x{)}^2}{q(x{)}^2})-E_{x\sim p(x)}(f(x))\cdot E_{x\sim p(x)}(f(x))$$

配方：配出一个在$p(x)$分布下的方差$D_{x\sim p(x)}(f(x))$，同时，我们称后面的差值为$\mathrm{△}D$。

$$D_{x\sim q(x)}(f(x)\frac{p(x)}{q(x)})=D_{x\sim p(x)}(f(x))+(E_{x\sim q(x)}(f(x{)}^2\frac{p(x{)}^2}{q(x{)}^2})-E_{x\sim p(x)}(f(x{)}^2))=D_{x\sim p(x)}(f(x))+\mathrm{△}D$$

要使得在新分布$q(x)$中采样的方差比原分布$p(x)$中采样的方差更小，则需要$\mathrm{△}D<0$。

$$\mathrm{△}D=\int f(x{)}^{2}p(x)\cdot \frac{p(x)}{q(x)}dx-\int f(x{)}^{2}p(x)dx$$

由于$f(x{)}^2>=0$并且$p(x)$是概率密度函数，$p(x)>=0$，所以$f(x{)}^{2}p(x)>=0$。 使得$\mathrm{△}D<0$的条件，需要$\frac{p(x)}{q(x)}<1$，即$p(x)<q(x)$。

然而，$p(x)$和$q(x)$是概率密度函数，它们满足在积分域上和为 1，所以$p(x)<q(x)$不可能在积分域上处处成立。

$$\int p(x)dx=1\text{and}\int q(x)dx=1$$

我们放宽条件：在积分域上重要的一部分区域上，$p(x)<q(x)$，而在另一部分区域上，$p(x)\ge q(x)$。 注意到公式中非负的权重$f(x{)}^{2}p(x)$，很像$f(x)$的二阶矩$E_{x\sim p(x)}(f(x{)}^2)$。

因此在权重$f(x{)}^{2}p(x)$大的位置上，我们调高$q(x)$，使得在该区域上满足$p(x)<q(x)$。而在权重$f(x{)}^{2}p(x)$小的地方，我们调低$q(x)$，使得在该区域上满足$p(x)\ge q(x)$。最终使得$\mathrm{△}D<0$。

以 VNDF 为例讨论快速蒙特卡洛收敛

上一节中讲到：选择一个合适的采样分布，“在权重$f(x{)}^{2}p(x)$大的位置上，调高$q(x)$，以降低方差”，能够使我们在较少的采样成本下快速接近准确值。

回顾一下渲染方程：

$$L_o(v)={\int }_{l\in \mathrm{\Omega }}BRDF\cdot L_i(l)(n\cdot l{)}^{+}dl$$

其中$f=BRDF\cdot L_i(l)(n\cdot l{)}^{+}$正是我们的采样函数（Always not Negative），而$p$是入射光线$l$的原始分布，在渲染方程中，我们对入射光线$l$进行均匀采样，所以$p$实际上是均匀分布。 我们的目标是找到一个分布$q$，使得在$f=BRDF\cdot L_i(l)(n\cdot l{)}^{+}$的位置上有更大的概率密度。

接下来，让我们看看 VNDF 是怎么做的。

VNDF（可见几何法向分布）

VNDF（Visible Normal Distribution Function）和通常的 NDF（Normal Distribution Function）最大的差异在于：在半球面上，只考虑了面朝摄像机一侧的微表面法向，背朝摄像机的法向占比为 0。

图中，NDF 中 15%面朝摄像机的微表面法向（$f>0$），在 VNDF 中占比为 30%。而 NDF 中 18%背朝摄像机的微表面法向（$f=0$），在 VNDF 中占比为 0%（不纳入 VNDF 考虑）。 因此 NVDF 相比于通常的 NDF 是更优的采样分布。

采样 VNDF 的技巧 - 其 1

在Understanding the Masking-Shadowing Function in Microfacet-Based BRDFs, Heitz这篇文章中，作者给出了很有巧思的采样方法。 文中使用的 VNDF:

$$D_{{\omega }_i}({\omega }_m)=\frac{G_1({\omega }_i,{\omega }_m)|{\omega }_i\cdot {\omega }_m|D({\omega }_m)}{|{\omega }_i\cdot {\omega }_g|}$$

它们将 VNDF 建模为椭球体（包含粗糙度和各向异性特性），通过对椭球体在观察向量$v$方向的投影区域进行采样，可以精确地得到半程向量$h$（或者叫做微表面法向量${\omega }_m$）。

具体步骤：

1. 图 a：椭球体和它的投影区域
2. 图 b：记正变换 T，它将椭球体变换到标准球面，投影区域也随之改变。
3. 图 c：对投影区域进行均匀采样，得到相对应于标准球面上的法向量$N_h$
4. 图 d：对法向量$N_h$作逆变换$T^{-1}$，得到椭球体上的法向量$N_e$。

怎么这么像 LTC，这种“转换进 cardinal 标准空间处理问题”的思路很流行呢。

采样 VNDF 的技巧 - 其 2

在Sampling Visible GGX Normals with Spherical Caps, Jonathan Dupuy, Anis Benyoub这篇文章中，作者给出了另一种采样方法。

这种方法来自一个很巧妙的观察：当我们看向一个完全镜面球体时，根据仰角的不同，反射光线均匀的分布在一个球冠表面上。记球坐标仰角$\varphi$，球冠表面的高度是$-\cos \varphi$。

因此我们在高度为$-\cos \varphi$的球冠表面上均匀采样出反射向量$l$，结合观察向量$v$计算得到服从 VNDF 分布的微表面法向${\omega }_m$。

采样 VNDF 的技巧 - 其 3