Monte Carlo 积分方法

Monte Carlo 积分方法是一种求解定积分的数值方法．假定需要计算的是下面的积分

\[ I=\int_a^b f(x)dx. \]

先选取一个连续型随机变量 \(X\)，要求其概率密度函数 \(p(x)\) 满足

\[ \forall x\in[a,b],\quad p(x)>0. \]

然后令 \(Y=g(X)=\frac{f(X)}{p(X)}\)．注意到

\[ \mathrm{E}\left[Y\right]=\mathrm{E}\left[\frac{f(X)}{p(X)}\right]=\int_{a}^{b} \frac{f\left(x\right)}{p\left(x\right)} p\left(x\right) d x=I, \]

我们可以将求解定积分 \(I\) 的任务，直接转化为估计随机变量 \(Y\) 的期望 \(\mathrm{E}\left[Y\right]\).

估计 \(\mathrm{E}\left[Y\right]\) 的一个常见做法，就是用 \(Y\) 为总体进行简单随机抽样，然后以样本均值

\[ \overline{Y}=\frac{1}{N}\sum_{i=1}^N Y_i \]

作为 \(\mathrm{E}\left[Y\right]\) 的估计量．具体方案如下：

设 \(X_i\sim p(x)(i=1,2,\cdots,n)\) 是独立同分布的 \(n\) 个随机变量，则 \(Y_i=\frac{f(X_i)}{p(X_i)}\) 也是独立同分布的 \(n\) 个随机变量．令 \[ I_N=\frac{1}{N}\sum_{i=1}^N\frac{f(X_i)}{p(X_i)}=\frac{1}{N}\sum_{i=1}^N Y_i, \] 以 \(I_N\) 作为 \(\mathrm{E}\left[Y\right]=I\) 的估计量．

统计学告诉我们：样本均值是总体期望的一致无偏估计量．这里我们可以再次证明一下这个结论．首先，\(I_N\) 的期望为

\[ \begin{aligned} \mathrm{E}\left[I_{N}\right]&=\mathrm{E}\left[\frac{1}{N}\sum_{i=1}^N Y_i\right]=\frac{1}{N}\mathrm{E}\left[\sum_{i=1}^N Y\right] = \mathrm{E}\left[Y\right] =I, \end{aligned} \]

这说明 \(I_N\) 是无偏的．其次，当 \(N\to\infty\) 时，\(I_N\) 的方差

\[ \begin{aligned} \mathrm{Var}\left[I_{N}\right]&=\mathrm{Var}\left[\frac{1}{N} \sum_{i=1}^{N} Y_i\right]\\ &=\frac{1}{N^2} \sum_{i=1}^{N} \mathrm{Var}\left[Y_i\right] \\ &=\frac{1}{N} \mathrm{Var}\left[Y\right] \to 0, \end{aligned} \]

这说明 \(I_N\) 是一致的．证明完毕．

重要性采样

对于一个估计量，我们总希望减小它的方差．估计量 \(I_{N}\) 的方差表达式为

\[ \mathrm{Var}\left[I_{N}\right]=\frac{1}{N} \mathrm{Var}\left[\frac{f(X)}{p(X)}\right]. \]

我们清楚的是：当 \(p(x)=kf(x)\) 时，也即 \(p(x)=\frac{1}{I}f(x)\) 时，方差会取得最小值 \(0\)．但是，求出 \(I\) 是不可能的，因为这正是我们 Monte Carlo 积分方法的求解目标．我们所能做的，只是让 \(p(x)\) 的形态尽量与 \(f(x)\) 相接近．这种削减方差的手段叫作 重要性采样（importance sampling）．\(f(x)\) 取值越大的区域，对积分的贡献会越多，也就越重要，需要更多的采样．

如果我们简单地让 \(X\) 服从 \([a,b]\) 上的均匀分布．此时估计量为

\[ I_N=\frac{b-a}{N}\sum_{i=1}^Nf(X_i). \]

直觉上看，这就是将积分区间等分成 \(N\) 段，用 \(N\) 个面积为 \[\frac{b-a}{N}\cdot f(X_i)\] 的矩形来估计曲边梯形．如果我们希望在某个子区间采样数翻倍，那么这个子区间需要分成两半，将一个矩形换成宽度为 \(\frac{b-a}{2N}\) 的两个矩形．推而广之就是采样概率越大的区域，单个样本的权重就越小．这也是为什么 \(I_{N}\) 中要除以 \(p(X_i)\) 的原因．

采样余弦

给定法向量 \(\mathbf{n}\)，我们用 \(H^2_{\mathbf{n}}\) 表示底部大圆与 \(\mathbf{n}\) 垂直的上半球面．考虑对积分

\[ I=\int_{H^2_{\mathbf{n}}}f(\boldsymbol{\omega}_i)(\mathbf{n}\cdot\boldsymbol{\omega}_i)d\Omega_i. \]

进行重要性采样．

注意到

\[ \int_{H^2_{\mathbf{n}}}\frac{\mathbf{n}\cdot\boldsymbol{\omega}_i}{\pi}\:d\Omega_{i}=1, \]

可以设 \(\eta_{\mathbf{n}}\) 是半球 \(H^2\) 上的一个概率分布，其概率密度函数为

\[ pdf_{\mathbf{n}}(\boldsymbol{\omega}_i) =\frac{\mathbf{n}\cdot\boldsymbol{\omega}_i}{\pi}. \]

我们利用 \[\eta_{\mathbf{n}}\] 进行重要性采样，由此得到 \(I\) 的统计量

\[ \widehat{I}=\frac{\pi}{N}\sum_{k=1}^Nf\left(\boldsymbol{\omega}_i^{(k)}\right),\quad \boldsymbol{\omega}_i^{(k)}\sim\eta_{\mathbf{n}}, \]

其中 \(N\) 为采样数．下面我们将具体描述如何为 \(\eta_{\mathbf{n}}\) 生成样本．

我们首先讨论如何为 \(\eta_{(0,0,1)}\) 生成样本．

使用逆变换法可以证明：若 \(\xi_{1}^{(k)}, \xi_{2}^{(k)}\) 独立且服从均匀分布 \(U(0,1)\)，则 \[ \boldsymbol{\nu}_i^{(k)}=\left(\sqrt{\xi_{1}} \cos \left(2 \pi \xi_{2}\right), \sqrt{\xi_{1}} \sin \left(2 \pi \xi_{2}\right), \sqrt{1-\xi_{1}}\right) \] 服从分布 \(\eta_{(0,0,1)}\)．

对于一般的 \(\mathbf{n}\)，我们会提供两种采样 \(\eta_{\mathbf{n}}\) 的方案：坐标变换法和四元数旋转法．

坐标变换法：

设

\[ \mathbf{u}= \begin{cases} (0,0,1) &\mathrm{if}\; \mathbf{n}\ne(0,0,\pm1), \\ (1,0,0) & \mathrm{if}\; \mathbf{n}=(0,0,\pm1). \end{cases} \]

切空间的基向量在世界空间下的矩阵 \(T_{\mathbf{n}}\)可以表示成

\[ T_{\mathbf{n}}=\left(\mathbf{t},\mathbf{b},\mathbf{n}\right), \]

其中

\[ \begin{aligned} \mathbf{t}&=\frac{\mathbf{u}\times\mathbf{n}}{\Vert\mathbf{u}\times\mathbf{n}\Vert},\\ \mathbf{b}&=\mathbf{n}\times\mathbf{t}. \end{aligned} \]

于是我们得到了 \(\eta_{\mathbf{n}}\) 的采样方案．

采样方案 设 \(\xi_{1}^{(k)}, \xi_{2}^{(k)}\) 独立且服从均匀分布 \(U(0,1)\)，则

\[ \boldsymbol{\omega}_i^{(k)}=T_{\mathbf{n}}\boldsymbol{\nu}_i^{(k)}=T_{\mathbf{n}} \begin{bmatrix} \sqrt{\xi_{1}} \cos \left(2 \pi \xi_{2}\right)\\ \sqrt{\xi_{1}} \sin \left(2 \pi \xi_{2}\right)\\ \sqrt{1-\xi_{1}} \end{bmatrix} \]

服从概率密度为

\[ pdf_{\mathbf{n}}(\boldsymbol{\omega}_i) =\frac{\mathbf{n}\cdot\boldsymbol{\omega}_i}{\pi}. \]

的分布．

四元数旋转法：

• 当 \(\mathbf{n}=(0,0,1)\) 时，我们可以使用逆变换采样．设 \(\xi_{1}^{(k)}, \xi_{2}^{(k)}\) 独立且服从均匀分布 \(U(0,1)\)，则 \[ \boldsymbol{\nu}_i^{(k)}=\left(\sqrt{\xi_{1}} \cos \left(2 \pi \xi_{2}\right), \sqrt{\xi_{1}} \sin \left(2 \pi \xi_{2}\right), \sqrt{1-\xi_{1}}\right) \] 服从分布 \(\eta_{(0,0,1)}\)．

• 当 \(\mathbf{n}=(0,0,-1)\) 时，\(-\boldsymbol{\nu}_i^{(k)}\) 服从分布 \(\eta_{(0,0,-1)}\)．

• 当 \(\mathbf{n}=(n_x,n_y,n_z)\ne(0,0,\pm 1)\) 时，将向量 \((0,0,1)\) 沿旋转轴 \[ (0,0,1)\times\mathbf{n}=(-n_y,n_x,0) \] 旋转成向量 \(\mathbf{n}\) 的四元数为 \[ q_\mathbf{n} = \cos\frac{\alpha}{2}+\sin\frac{\alpha}{2}\left(-\frac{n_y}{\sqrt{n_x^2+n_y^2}}i+\frac{n_x}{\sqrt{n_x^2+n_y^2}}j\right), \] 其中 \(\alpha\) 是旋转角．因为 \[ \cos\alpha=(0,0,1)\cdot \mathbf{n}=n_z. \] 我们有 \[ \begin{aligned} q_\mathbf{n}&=\sqrt{\frac{1+\cos\alpha}{2}}+\sqrt{\frac{1-\cos\alpha}{2}}\left(-\frac{n_y}{\sqrt{n_x^2+n_y^2}}i+\frac{n_x}{\sqrt{n_x^2+n_y^2}}j\right)\\ &=\sqrt{\frac{1+n_z}{2}}+\sqrt{\frac{1-n_z}{2}}\left(-\frac{n_y}{\sqrt{n_x^2+n_y^2}}i+\frac{n_x}{\sqrt{n_x^2+n_y^2}}j\right). \end{aligned} \] 由此可以计算出四元数旋转变换 \(\mathrm{rot}_q:p\mapsto qpq^{-1}\) 对应的旋转矩阵 \[ R_\mathbf{n}= \begin{pmatrix} \tfrac{n_x^2n_z+n_y^2}{n_x^2+n_y^2}&-n_xn_y(1+n_z)&n_x\\ -n_xn_y(1+n_z)&\tfrac{n_x^2+n_y^2n_z}{n_x^2+n_y^2}&n_y\\ -n_x&-n_y&n_z \end{pmatrix}. \] 不难看出，\(\boldsymbol{\omega}_i^{(k)}=R_\mathbf{n}\boldsymbol{\nu}_i^{(k)}\) 服从分布 \(\eta_{\mathbf{n}}\)．

采样 GGX 法线分布

统计量

给定法向量 \(\mathbf{n}\)，我们用 \(H^2_{\mathbf{n}}\) 表示底部大圆与 \(\mathbf{n}\) 垂直的上半球面．考虑对积分

\[ I=\int_{H^2_{\mathbf{n}}}f(\boldsymbol{\omega}_i)D(\mathbf{h})(\mathbf{n}\cdot\boldsymbol{\omega}_i)d\Omega_i, \]

进行重要性采样,其中

\[ D(\mathbf{h})=\dfrac{\alpha^{2}}{\pi\left((\mathbf{n} \cdot \mathbf{h})^{2}\left(\alpha^{2}-1\right)+1\right)^{2}},\quad \mathbf{h}=\frac{\boldsymbol{\omega}_i+\boldsymbol{\omega}_o}{\Vert\boldsymbol{\omega}_i+\boldsymbol{\omega}_o\Vert}. \]

根据面积元关系

\[ 4(\mathbf{h} \cdot \boldsymbol{\omega}_o)d\Omega_\mathbf{h}=d\Omega_i, \]

\(I\) 可以改写为关于 \(\mathbf{h}\) 的积分

\[ \begin{align*} I&=\int_{S^2}f(\boldsymbol{\omega}_i)D(\mathbf{h})(\mathbf{n}\cdot\boldsymbol{\omega}_i)^+d\Omega_i\\ &=\int_{H^2_{\mathbf{n}}}4f(\boldsymbol{\omega}_i)D(\mathbf{h})(\mathbf{h} \cdot \boldsymbol{\omega}_o)(\mathbf{n}\cdot\boldsymbol{\omega}_i)^+d\Omega_\mathbf{h}. \end{align*} \]

注意到

\[ \int_{H^2_{\mathbf{n}}}D(\mathbf{h})(\mathbf{n}\cdot\mathbf{h}) \:d\Omega_{\mathbf{h}}=1, \]

可以设 \(\eta_{\mathbf{n}}\) 是半球 \(H_{\mathbf{n}}^2\) 上的一个概率分布，其概率密度函数为

\[ pdf_{\mathbf{n}}(\mathbf{h}) =D(\mathbf{h})(\mathbf{n}\cdot\mathbf{h}). \]

我们利用 \(\eta_{\mathbf{n}}\) 进行重要性采样，由此得到 \(I\) 的统计量

\[ \widehat{I}=\frac{4}{N}\sum_{k=1}^N\frac{\left(\mathbf{h}^{(k)} \cdot \boldsymbol{\omega}_o\right)\left(\mathbf{n}\cdot\boldsymbol{\omega}_i^{(k)}\right)^+}{\mathbf{n}\cdot\mathbf{h}^{(k)}}f\left(\boldsymbol{\omega}_i^{(k)}\right),\quad \mathbf{h}^{(k)}\sim\eta_{\mathbf{n}}, \]

其中 \(N\) 为采样数，

\[\boldsymbol{\omega}_i^{(k)}=\mathbf{r}(\boldsymbol{\omega}_o,\mathbf{h}^{(k)})=2(\boldsymbol{\omega}_o\cdot\mathbf{h}^{(k)})\mathbf{h}^{(k)}-\boldsymbol{\omega}_o.\]

下面我们将具体描述如何为 \(\eta_{\mathbf{n}}\) 生成样本 \(\mathbf{h}^{(k)}\)．

采样方案

为应用逆变换采样法，先将 \(\mathbf{h}^{(k)}\) 在 \(\mathbf{n}\) 处的切空间上参数化为球坐标

\[ (\theta_{\mathbf{h}^{(k)}},\varphi_{\mathbf{h}^{(k)}})\in[0,\pi/2]\times[0,2\pi]. \]

设 \(\xi_{1}^{(k)}, \xi_{2}^{(k)}\) 独立且服从均匀分布 \(U(0,1)\)．由

\[ cdf_{\theta,\varphi}\left(\theta_{\mathbf{h}^{(k)}},\varphi_{\mathbf{h}^{(k)}}\right)=\int_{0}^{\varphi_{\mathbf{h}^{(k)}}} \mathrm{d} \varphi_{\mathbf{h}} \int_{0}^{\theta_{\mathbf{h}^{(k)}}} \frac{\alpha^{2} \cos \theta_{\mathbf{h}} \sin \theta_{\mathbf{h}} }{\pi\left(\left(\alpha^{2}-1\right) \cos ^{2} \theta_{\mathbf{h}} +1\right)^{2}} d\theta_{\mathbf{h}} \]

可知 \(\theta_{\mathbf{h}^{(k)}},\varphi_{\mathbf{h}^{(k)}}\) 是独立的．于是，我们可以分别求出 \(\theta_{\mathbf{h}^{(k)}},\varphi_{\mathbf{h}^{(k)}}\) 的边缘分布

\[ \begin{aligned} cdf_{\theta}\left(\theta_{\mathbf{h}^{(k)}}\right)&=cdf_{\theta,\varphi}\left(\theta_{\mathbf{h}^{(k)}},2\pi\right)\\ &=\int_{0}^{2\pi} \mathrm{d} \varphi_{\mathbf{h}} \int_{0}^{\theta_{\mathbf{h}^{(k)}}} \frac{\alpha^{2} \cos \theta_{\mathbf{h}} \sin \theta_{\mathbf{h}} }{\pi\left(\left(\alpha^{2}-1\right) \cos ^{2} \theta_{\mathbf{h}} +1\right)^{2}} d\theta_{\mathbf{h}}\\ &=\int_{0}^{\theta_{\mathbf{h}^{(k)}}} \frac{2\alpha^{2} \cos \theta_{\mathbf{h}} \sin \theta_{\mathbf{h}} }{\left(\left(\alpha^{2}-1\right) \cos ^{2} \theta_{\mathbf{h}} +1\right)^{2}} d\theta_{\mathbf{h}}\\ &=\left.\frac{\alpha^{2}}{\alpha^{2}-1} \frac{1}{\left(\alpha^{2}-1\right) \cos^2\theta_{\mathbf{h}}+1}\right|_{0} ^{\theta_{\mathbf{h}^{(k)}}}\\ &=\frac{\alpha^{2}}{\alpha^{2}-1} \frac{1}{\left(\alpha^{2}-1\right)\cos^2\theta_{\mathbf{h}^{(k)}}+1}-\frac{1}{\alpha^{2}-1} \\ &=\frac{1}{\alpha^{2}-1} \frac{\alpha^2-1-\left(\alpha^{2}-1\right)\cos^2\theta_{\mathbf{h}^{(k)}}}{\left(\alpha^{2}-1\right)\cos^2\theta_{\mathbf{h}^{(k)}}+1} \\ &= \frac{1-\cos^2\theta_{\mathbf{h}^{(k)}}}{\left(\alpha^{2}-1\right)\cos^2\theta_{\mathbf{h}^{(k)}}+1}, \\ \end{aligned} \]

\[ \begin{aligned} cdf_{\varphi}\left(\varphi_{\mathbf{h}^{(k)}}\right)&=cdf_{\theta,\varphi}\left(\frac{\pi}{2},\varphi_{\mathbf{h}^{(k)}}\right)\\ &=\int_{0}^{\varphi_{\mathbf{h}^{(k)}}} \mathrm{d} \varphi_{\mathbf{h}} \int_{0}^{\frac{\pi}{2}} \frac{\alpha^{2} \cos \theta_{\mathbf{h}} \sin \theta_{\mathbf{h}} }{\pi\left(\left(\alpha^{2}-1\right) \cos ^{2} \theta_{\mathbf{h}} +1\right)^{2}} d\theta_{\mathbf{h}}\\ &=\frac{1}{2\pi}\varphi_{\mathbf{h}^{(k)}}. \end{aligned} \]

使用 cdf 的逆变换，就可得到球坐标下 \(\mathbf{h}^{(k)}\) 的采样方案

\[ \begin{aligned} \theta_{\mathbf{h}^{(k)}}&=\arccos\sqrt{\frac{1-\xi_{1}^{(k)}}{1-(1-\alpha^2)\xi_{1}^{(k)}}},\\ \varphi_{\mathbf{h}^{(k)}}&=2\pi\xi_{2}^{(k)}. \end{aligned} \]

因为上式的球坐标是在切空间中生成的，我们需要先将它们转换到笛卡尔坐标系下，再用坐标变换法转成成世界空间下的坐标．沿用之前定义的 \(T_{\mathbf{n}}\)，我们得到了世界空间下 \(\mathbf{h}^{(k)}\) 的采样方案．

采样方案 设 \(\xi_{1}^{(k)}, \xi_{2}^{(k)}\) 独立且服从均匀分布 \(U(0,1)\)，则

\[ \mathbf{h}^{(k)}=T_{\mathbf{n}} \begin{bmatrix} \sin\theta_{\mathbf{h}^{(k)}}\cos\varphi_{\mathbf{h}^{(k)}}\\\\ \sin\theta_{\mathbf{h}^{(k)}}\sin\varphi_{\mathbf{h}^{(k)}}\\\\ \cos\theta_{\mathbf{h}^{(k)}} \end{bmatrix}=T_{\mathbf{n}} \begin{bmatrix} \alpha\cos2\pi\xi_{2}^{(k)}\sqrt{\dfrac{\xi_{1}^{(k)}}{1-(1-\alpha^2)\xi_{1}^{(k)}}}\;\\ \\ \alpha\sin2\pi\xi_{2}^{(k)}\sqrt{\dfrac{\xi_{1}^{(k)}}{1-(1-\alpha^2)\xi_{1}^{(k)}}}\\ \\ \sqrt{\dfrac{1-\xi_{1}^{(k)}}{1-(1-\alpha^2)\xi_{1}^{(k)}}} \end{bmatrix}. \]

服从概率密度函数为

\[ pdf_{\mathbf{n}}(\mathbf{h}) =D(\mathbf{h})(\mathbf{n}\cdot\mathbf{h}) \]

的分布．

多重重要性采样

假定需要计算的是下面的积分

\[ I=\int_a^b f(x)g(x)dx, \]

并且已知概率密度函数 \(p_X(x)\) 和 \(p_Y(x)\) 分别能较好地吻合 \(f(x)\) 和 \(g(x)\) 的形状，即它们分别能对

\[ \int_a^b f(x)dx,\,\int_a^b g(x)dx \]

进行较好地重要性采样，那么该如何对 \(I\) 进行重要性采样呢？

很容易想到的一种方法是用 \(p_X(x)p_Y(x)\) 进行重要性采样，但是，有些时候为总体 \(Z\sim p_X(x)p_Y(x)\) 生成样本 \(Z_i\) 比较困难．这时，多重重要性采样（multiple importance sampling）不失为一种好的选择．

设 \(X_i\sim p_X(x)\)，\(Y_i\sim p_Y(x)\)，多重重要性采样使用了如下估计量

\[ I_{n_{f},n_g}=\frac{1}{n_{f}} \sum_{i=1}^{n_{f}} \frac{f\left(X_{i}\right) g\left(X_{i}\right) w_{f}\left(X_{i}\right)}{p_{X}\left(X_{i}\right)}+\frac{1}{n_{g}} \sum_{j=1}^{n_{g}} \frac{f\left(Y_{j}\right) g\left(Y_{j}\right) w_{g}\left(Y_{j}\right)}{p_{Y}\left(Y_{j}\right)}, \]

其中 \(w_f(x),w_g(x)\) 是满足 \(w_f(x)+w_g(x)=1\) 的权重函数．可以验证

\[ \begin{aligned} \mathrm{E}\left[I_{n_{f},n_g}\right]&= \mathrm{E}\left[\frac{f\left(X\right) g\left(X\right) w_{f}\left(X\right)}{p_{X}\left(X\right)}\right]+\mathrm{E}\left[\frac{f\left(Y\right) g\left(Y\right) w_{g}\left(Y\right)}{p_{Y}\left(Y\right)}\right]\\ &=\int_a^b\frac{f\left(x\right) g\left(x\right) w_{f}\left(x\right)}{p_X(x)}p_X(x)dx+\int_a^b\frac{f\left(x\right) g\left(x\right) w_{g}\left(x\right)}{p_Y(x)}p_Y(x)dx\\ &=\int_a^bf\left(x\right) g\left(x\right) w_{f}\left(x\right)dx+\int_a^bf\left(x\right) g\left(x\right) w_{g}\left(x\right)dx\\ &=\int_a^bf\left(x\right) g\left(x\right) (w_{f}\left(x\right)+w_{g}\left(x\right))\:dx\\ &=\int_a^bf\left(x\right) g\left(x\right)dx, \end{aligned} \]

即 \(I_{n_{f},n_{g}}\) 确实是 \(I\) 的无偏估计量．实践中，常使用如下形式的启发式权重函数

\[ w_{f}(x)=\frac{\left(n_{f} p_{X}(x)\right)^{\beta}}{\left(n_{f} p_{X}(x)\right)^{\beta}+\left(n_{g} p_{Y}(x)\right)^{\beta}},\quad w_{g}(x)=\frac{\left(n_{g} p_{Y}(x)\right)^{\beta}}{\left(n_{f} p_{X}(x)\right)^{\beta}+\left(n_{g} p_{Y}(x)\right)^{\beta}}. \]

当 \(\beta=2\) 时，往往效果会较好．