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

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

---

前言

之前看到一个问题：蒙特卡洛采样中的最优 PDF ，它满足什么条件时采样具有无偏性，什么条件下方差最小？

今天整理了一些相关的数学证明，和 VNDF 相关知识。

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

求$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)}$$

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

首先一个函数在某个分布上的期望可以表示为：

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

它可以通过过量的采样来近似计算：

$$E_{x\sim p(x)}(f(x))=\frac{1}{N}\sum_{x_i\in p(x)}^N f(x)$$

把两个方程联立，也就是说

$$\int f(x)p(x)dx=\frac{1}{N}\sum_{x_i\in p(x)}^N f(x)$$

变量代换：

$$\int g(x)dx=\frac{1}{N}\sum_{x_i\in p(x)}^N \frac{g(x)}{p(x)}$$

用渲染方程展开，这下看懂了：

$$L_o(v)=\int BRDF\cdot L_i(l) (n\cdot l)^+dl=\frac{1}{N}\sum_{l_i\in p(l)}^N \frac{BRDF\cdot L_i(l) (n\cdot l)^+}{p(l)}$$

如果我们的采样策略是在半球面上的均匀采样，那么概率密度函数是$P(l)=\frac{1}{2\pi}$，代入上式得到：

$$L_o(v)=\frac{1}{N}\sum_{l_i\in p(l)}^N BRDF\cdot L_i(l) (n\cdot l)^+\cdot 2\pi$$

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

根据$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))$，同时，我们称后面的差值为$\triangle D$。

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

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

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

由于$f(x)^2>=0$并且$p(x)$是概率密度函数，$p(x)>=0$，所以$f(x)^2p(x)>=0$。 使得$\triangle 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 \quad \text{and} \quad \int q(x)dx=1$$

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

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

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

关于 VNDF 的部分参考自Importance Sampling PDFs (VNDF, Spherical Caps)

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

回顾一下渲染方程：

$$L_o(v)=\int_{l\in\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) \vert \omega_i\cdot \omega_m \vert D(\omega_m)}{\vert \omega_i \cdot \omega_g \vert}$$

它们将 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这篇文章中，作者给出了另一种采样方法。

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

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

采样 VNDF 的技巧 - 其 3

Bounded VNDF Sampling for Smith–GGX Reflections

技巧其 3 相比于技巧其 2 的改进在于：增加了一个新的约束条件，排除了微表面法向$\omega_m$朝向摄像机但是反射光线$\omega_o$来自几何法平面以下的情况。

下表显示了 VNDF 和 Bounded VNDF 的条件。 普通的 VNDF 只要满足微表面法向$\omega_m$朝向视向量（$(\omega_i \cdot \omega_m) > 0$）即可，而 Bounded VNDF 则需要额外满足反射光线$\omega_o$不得低于几何法平面$(\omega_i \cdot \omega_g )> 0$。

	${(\omega_i \cdot \omega_m)}>0$	${(\omega_i \cdot \omega_m)}<0$
${(\omega_o \cdot \omega_g)}>0$	$\text{ Bounded VNDF} $
${(\omega_o \cdot \omega_g)}<0$
	VNDF

新的采样方法带来新的 Bounded VNDF 分布。他和原有的 VNDF 分布稍有不同，文中也推导了 Bounded VNDF 的概率密度函数。

Note

仔细地读了一下原文，跟着作者推导了一遍。 文中新加的限制是：反射光线$o$不得低于法平面以下$(o\cdot n)>0$。这个新的边界在 Unstretched 空间中是围绕法平面的圆圈。而在 Stretched 空间中会在球面上上下变化。 推导出新的边界在 Stretched 空间中的表达式： $$\acute{o}_z(\theta=\frac{\pi}{2},\phi)=\Big( \frac{(i_x+\cos{\phi})^2+(i_y+\sin{\phi})^2+i_z^2}{\frac{(i_x+\cos{\phi})^2}{\alpha_x^2}+\frac{(i_x+\cos{\phi})^2}{\alpha_y^2}+i_z^2}-1 \Big)\acute{i}_z$$ 这个新的编辑如下图红线所示。他有一个最小值，正是文中提到的$-k\acute{i}_z$。 我们采样最小值$-k\acute{i}_z$到球盖上顶点（1）作为反射向量的 z 分量$z\in(-k\acute{i}_z, 1)$。 所以最终采样出来的反射向量还是有一部分在法平面以下的，因为并非所有的采样点都位于新的边界以上（红线以上）。