重要性采样理论与多重重要性采样

重要性采样理论

译自:Physicall Based Rendering,Third Edition, Importance Sampling

额外补充了部分内容。

重要性采样在光线追踪中可以有效地加速噪点方差的降低进程。我们知道,在蒙托卡洛积分中,有:
$$
F_N=\frac{1}{N}\sum_{i=1}^{N}\frac{f(X_i)}{p(X_i)}
$$
当采样函数$p(x)$的变化趋势与被积函数$f(x)$越相近,蒙托卡洛积分的收敛速度就会越快。基于此,我们可以部分舍弃被积函数较小的部分,而将采样点集中在被积函数数值较高的部分,以此来得到更准确且更有效率的采样估算。

左图:均匀采样		右图:对光源的重要性采样

左图为均匀采样,右图为对光源的重要性采样

举个例子。假设现在有BSDF散射方程:
$$
L_o(p,\omega_o)=\int_{S^2}f(p,\omega_o,\omega_i)L_i(p,\omega_i)|cos\theta_i|d\omega_i
$$
当进行积分积分近似时,如果我们随机选取到了一个几乎和表面法线垂直的方向,它的余弦项$cos\theta$将接近于0。此时,这一轮BSDF的辐射度函数将基本上被浪费了,因此它对这一点辐射度的整体贡献十分微不足道。一般来说,如果我们能够通过控制采样函数,减少这种特殊方向的采样率,收敛速度将会明显提升。更一般地来讲,如果从与积分函数本身十分相近的分布函数中进行采样,同样可以提高收敛速度。

只要从与被积函数相近的概率分布中采样,就能显著降低采样噪点方差。证明如下:我们假设目前要计算一个积分,这个积分满足:
$$
I=\int_a^bf(x)dx
$$
其中,$f(x)$是一个单峰函数,且在区间$[a,b]$上存在概率分布函数$p(x)$,满足$p(x)\ge 0$。此时,我们使用蒙特卡洛积分积分来估计$I$,从$p(x)$中采样$n$个随机变量$x_1,x_2,…,x_n$,计算可得:
$$
\hat I=\frac{1}{n}\sum_{i=1}^n\frac{f(x_i)}{p(x_i)}
$$
$\hat I$为$I$的估计量。此时,计算$\hat I$的方差,可以得到:
$$
Var(\hat I)=\frac{1}{n^2}\sum_{i=1}^nVar(\frac{f(x_i)}{p(x_i)})
$$
根据方差的定义,可以得出:
$$
Var(\hat I)=\frac{1}{n^2}\sum_{i=1}^n(E[\frac{f^2(x_i)}{p^2(x_i)}]-E[\frac{f(x_i)}{p(x_i)}]^2)
$$
其中,$E[x]$表示$x$的期望。此时,使用概率密度函数,可以进一步化简:
$$
Var(\hat I)=\frac{1}{n^2}\sum_{i=1}^n(\int_a^bp(x)dx-\hat I^2)
$$
整理得到:
$$
Var(\hat I)=\frac{1}{n}\int_a^b(\frac{f^2(x)}{p(x)}-\hat I^2)p(x)dx
$$
此时很容易就能看出,当$p(x)$与$f(x)$的形状越接近时,产生的方差值将会越小。应用到光线追踪中就是路径追踪的收敛速度将会越快。

由于在采样时我们可以自由地选择采样函数,我们可以选择一个$p(x)$满足$p(x)\propto f(x)$,或者说$p(x)=cf(x)$,不难发现,
$$
c=\frac{1}{\int f(x)dx}
$$
要得到“形状类似于积分函数”的概率密度函数,就需要先得到这个积分函数;然而,我们最初的目标就是通过近似的方法来得到这个积分函数。不过,对于这个例子,如果我们能够从这个分布中进行抽样,我们就可以继续下去。假设$f(x)$和$p(x)$的形状完全一致,那末就能得到:
$$
\frac{f(X_i)}{p(X_i)}=\frac{1}{c}=\int f(x)dx
$$
由于$c$是一个常量,那么在我们的假设中,每一次采样得到的价值都相等——此时我们可以惊喜地发现,采样带来的方差是0。当然,这肯定是不合理的,因为我们甚至都还没有进行蒙特卡洛积分。然而,这也从另一方面说明了,如果我们能找到形状与$f(x)$相似的概率密度函数$p(x)$,采样方差就会减少。反之,如果采用了不当的分布,重要性采样则会使得方差增大。

多重重要性采样

蒙特卡洛积分提供了估算$\int f(x)dx$的方法,然而,我们会频繁地需要计算类似于$\int f(x)g(x)$形式的函数。如果我们对于$f(x)$和$g(x)$这两个函数分别各有一个重要性采样策略,如果我们不能很轻松地合并两个策略,应该采取其中的哪一个呢?正如上文所讨论的,一个糟糕的采样策略可能比均匀采样效率更低。

例如,考虑对直接光照的积分估量:
$$
L_o(p,\omega_o)=\int_{S^2}f(p,\omega_o,\omega_i)L_d(p,\omega_i)|cos\theta_i|d\omega_i
$$
如果我们使用对其中任一函数的重要性采样,另一函数的采样效率就会降低。

再考虑使用一个近似于镜面的BRDF区域光照$L_d$的分布来进行近似采样。因为BRDF近似于一个镜面,积分的值在啊除了镜面反射方向外都接近于0。这意味着对所有方向的采样$L_d$都会给出几乎为0的贡献,同时方差将会变得比较大。更糟糕的是,当光源变得更大、更多的采样点集中到光源的方向时,概率密度函数的值继续下降,此时,对于BRDF在采样方向上非零的罕见方向而言,我们将不得不以一个很大的积分数值去和一个很小的PDF值相除。虽然在BRDF的分布中进行采样会是处理光源的一个好方法,但是对于漫反射、光泽度较高的BRDF以及小型的光源而言,从BRDF的分布中进行采样同样可能导致比从光源的分布中进行采样导致更高的方差。

不幸的是,从每个分布中抽取一个样本,看似最显而易见的解决方案——并对两个估计策略进行平均——并不能解决问题,因为在这种情况下方差是可加的。一旦估计策略已经产生了方差,即使它本身的方差较低,我们也不能通过将其和另一个估计策略做平均来消除它。

多重重要性采样(Multiple importance sampling,MIS)正是为了解决这一问题而生的,这是一种简单且容易实现的方法。最基本的思想是,当估值积分时,我们需要从多个采样分布中进行采样,并期待其中至少有一个分布可以合理地匹配被积函数的形状,即使我们并不知道哪一个分布是最合适的。MIS提供了一种方法来对每个采样策略的样本进行加权,从而消除由于积分和采样密度之间不匹配而导致的大方差峰值。这一采样方法甚至鼓励专门针对特殊情况构造特化的采样策略,因为它们可以在发生这些情况时降低方差,并通常来说成本开销不大。

在MIS中,如果有两个采样策略$p_f,p_g$被用于对$\int f(x)g(x)dx$进行采样,那么MIS给出的新蒙特卡洛积分式为:
$$
\frac{1}{n_f}\sum_{i=1}^{n_f}\frac{f(X_i)g(X_i)\omega_f(X_i)}{p_f(X_i)}+\frac{1}{n_g}\sum_{i=j}^{n_g}\frac{f(Y_j)g(Y_j)\omega_g(Y_j)}{p_g(Y_j)}
$$
其中,$n_f$是在$p_f$分布中的采样次数,$n_g$是在$p_g$中的采样次数,$\omega_f$和$\omega_g$是特殊的加权函数,使得估算策略的期望值等于被估计的被积函数$f(x)g(x)$。

权重函数考虑了样本$X_i,Y_j$的所有不同生成方式,而不是使用其中特定的一个。一个不错的权重函数是平衡启发式:
$$
\omega_s(x)=\frac{(n_sp_s(x))}{\sum_{i}(n_ip_i(x))}
$$
平衡启发式是一种可以验证的减少方差的优秀权重采样方法。

考虑这种方式对从概率分布$p_f$中得到的一个样本$X$,此时假设$p_f(X)$较低,如果$p_f$可以很好的匹配$f(x)$的形状,那么$f(x)$的价值也会较低;但如果$g(x)$的价值偏高,此时标准的重要性采样的估计结果:
$$
\frac{f(X)g(X)}{p_f(X)}
$$
将会偏高,因为$p_f(X)$的值较低,这也导致了较高的方差。而通过平衡启发式,$X$点的贡献将会是:
$$
\frac{f(X)g(X)\omega_f(X)}{p_f(X)}=\frac{f(X)g(X)n_f}{n_fp_f(X)+n_gp_g(X)}
$$
只要$p_g$的分布可以较好地匹配$g(x)$,那么由于$n_gp_g(X)$的存在,分母就不会过小,巨大的方差峰值也得以消除,即使$X$实际上是在以一个较差的匹配进行估计。事实上,另一个分布也会被用于采样,而这一新的分布将可以找到$X$在两种权重的平衡下产生的最大积分价值,以此来减少方差问题。

平衡启发式的一个简单实现是维护两个采样点和两个采样函数;在pbrt中,我们不需要考虑更一般的分布情况。

1
2
3
@ti.func
def BalanceHeuristic(nf: int, fpdf: float, ng: int, gpdf: float)->float:
return (nf * fpdf) / (nf * fpdf + ng * gpdf)

在实际情况中,我们往往会采用幂平衡启发式来进一步减少方差。对于指数$\beta$,构造它的平衡启发式:
$$
\omega_s(x)=\frac{(n_sp_s(x))^\beta}{\sum_{i}(n_ip_i(x))^\beta}
$$
而指数$\beta = 2$通常来说是一个不错的选择。

1
2
3
4
5
@ti.func
def BalanceHeuristic(nf: int, fpdf: float, ng: int, gpdf: float)->float:
f = nf * fpdf
g = ng * gpdf
return (f * f) / (f * f + g * g)

重要性采样理论与多重重要性采样
https://tridiamond.tech/post/ImportanceSampling.html
作者
想变成一只白泽
发布于
2023年3月31日
许可协议