展览会建设网站平台的作用,阿里巴巴做网站,做网站怎么写代码,二级网站排名做不上去目录一、前篇文章回顾二、积分的黎曼和形式三、积分的概率形式#xff08;蒙特卡洛积分#xff09;四、误差五、蒙特卡洛积分计算与收敛速度六、重要性采样七、重要性采样方法和过程八、重要性采样的优缺点一、前篇文章回顾 在前一篇文章3D数学系列之——从“蒙的挺准”到“蒙…
目录一、前篇文章回顾二、积分的黎曼和形式三、积分的概率形式蒙特卡洛积分四、误差五、蒙特卡洛积分计算与收敛速度六、重要性采样七、重要性采样方法和过程八、重要性采样的优缺点一、前篇文章回顾 在前一篇文章3D数学系列之——从“蒙的挺准”到“蒙的真准”解密蒙特卡洛积分中介绍了下面的这些朴素的公式 灰色部分中点数总的点数≈灰色部分面积矩形面积⇒灰色部分面积≈矩形面积×灰色部分中点数总的点数\begin{align} \cfrac{ 灰色部分中点数 } {总的点数} \approx \cfrac{灰色部分面积} {矩形面积} \\[2ex] \Rightarrow 灰色部分面积 \approx 矩形面积 \times \cfrac{ 灰色部分中点数 } {总的点数} \end{align} 总的点数灰色部分中点数≈矩形面积灰色部分面积⇒灰色部分面积≈矩形面积×总的点数灰色部分中点数
灰色部分面积≈矩形面积×PP为点落在灰色部分的概率灰色部分面积 \approx 矩形面积 \times P \quad P为点落在灰色部分的概率 灰色部分面积≈矩形面积×PP为点落在灰色部分的概率 虽然最终通过一个简单的概率乘法就可以计算出积分值但是如何得到这个精确的概率值从而得到比较精确的积分值却不是那么简单的。并且这个方法直接用于计算或者说模拟计算来说不确定的东西太多。下面我们就继续从积分的基本公式开始讨论一下看看有没有进一步更形式化和更精确一点的方法。
二、积分的黎曼和形式 对于数值积分计算来说积分的黎曼和形式是最基本的形式 ∫abf(x)dx≈∑n1Nf(xn)b−aN(1)\int \limits_a^b {f}(x) \mathrm{d} x \approx \sum \limits_{n1}^N {f}(x_n) \cfrac{b-a}{N} \tag{1} a∫bf(x)dx≈n1∑Nf(xn)Nb−a(1) 这是一个常见和基本的可以用于积分计算程序的基本公式。
三、积分的概率形式蒙特卡洛积分 接着我们对黎曼和使用著名的”陶哲轩瞪眼法“进行观察。 首先对于求和来说第一项是 f(xn){f}(x_n)f(xn) 这里我们一般会比较容易知道 f(x){f}(x)f(x) 的解析形式或者可以知道它的一些离散值。 在实际计算中需要生成的是 xnx_nxn 的序列这个序列是跟后面的 (b−a)N\cfrac{(b-a)}{N}N(b−a) 紧密关联的因为需要保证 xi−xi−1b−aNx_{i}-x_{i-1} \cfrac{b-a}{N}xi−xi−1Nb−a 即 xnx_nxn 序列是个等差数列间隔是就是 (b−a)N\cfrac{(b-a)}{N}N(b−a) 。 这时我们先打开思路既然前一讲中已经可以通过随机变量的形式来计算积分那么这时我们可不可以用一组随机数的序列来代替 xnx_nxn 的序列呢 既然讲到随机那么回忆一下概率论中关于“均匀分布随机变量”的概率密度函数The Probability Distribution Function缩写为 PDF p(x){1b−a,axb0,其他(2)\mathrm{p}(x) \begin{cases} \cfrac{1}{b-a}, \quad a \lt x \lt b \\[2ex] 0, \quad 其他 \end{cases} \tag{2} p(x)⎩⎨⎧b−a1,axb0,其他(2) 接着我们用陶哲轩瞪眼法观察黎曼和公式与上面这个公式此时我们应灵光乍现想到下面这样的形式 ∫abf(x)dx≈∑n1Nf(xn)b−aN1N∑n1Nf(xn)×(b−a)根据乘以一个数等于除以一个数的倒数这个法则有1N∑n1Nf(xn)p(xn)(3)\int \limits_a^b {f}(x) \mathrm{d} x \approx \sum \limits_{n1}^N {f}(x_n) \cfrac{b-a}{N} \\[2ex] \cfrac{1}{N} \sum \limits_{n1}^{N} {f}(x_n) \times (b-a) \\[2ex] 根据乘以一个数等于除以一个数的倒数这个法则有 \\[2ex] \cfrac{1}{N} \sum \limits_{n1}^{N} \cfrac{{f}(x_n)}{\mathrm{p}(x_n)} \tag{3} a∫bf(x)dx≈n1∑Nf(xn)Nb−aN1n1∑Nf(xn)×(b−a)根据乘以一个数等于除以一个数的倒数这个法则有N1n1∑Np(xn)f(xn)(3) 啊哈太好了这一切看似是公式间的一种巧合或者说技巧性的接着我们思考下它的意义。因为刚才我们已经说过原来的xnx_nxn 序列是个等差序列本身就是均匀分布的所以其概率密度函数PDF就是均匀分布变量的 PDF。如果我们取一个随机的 xnx_nxn 序列并假设知道它的 PDF p(xn)\mathrm{p}(x_n)p(xn) 那么就可以根据上面这个式子的结果来计算积分值了。并且当 xnx_nxn 不断增多时和 1N∑n1Nf(xn)p(xn)\cfrac{1}{N} \sum \limits_{n1}^{N} \cfrac{{f}(x_n)}{\mathrm{p}(x_n)}N1n1∑Np(xn)f(xn) 也就越接近原始的积分值这也就是传说中的蒙特卡洛积分的基本形式。它是合理且有意义的起码我们不用只靠前一讲中的“蒙的真准”的方法了。需要注意的就是其中 f(xn)p(xn)\cfrac{{f}(x_n)}{p(x_n)}p(xn)f(xn) 项看上去虽然很奇怪但他其实是有实际含义的一般被称作概率平均。 当然使用随机序列 xnx_nxn 是有条件的即要求 ∫−∞∞p(x)1(4)\int \limits_{-\infty}^{\infty} \mathrm{p}(x) 1 \tag{4} −∞∫∞p(x)1(4) 接着我们继续观察还是使用 “陶哲轩瞪眼法”那个带有 PDF 的蒙特卡洛积分公式发现其中还有个 1N\cfrac{1}{N}N1 , 这时我们应该突然明白其实这个式子还是个平均值的公式而讲到平均值学过概率论的同学就应该立刻想到“数学期望”。一般的数学期望的定义 E[f(X)]∑n1∞f(xn)p(xn)E[{f}(X)] \mathop{\sum}\limits_{n 1}^{\infty} {f}(x_n) \mathrm{p}(x_n) E[f(X)]n1∑∞f(xn)p(xn) 也就是说如果我们知道了随机变量 xnx_nxn 的随机 PDF p(xn)\mathrm{p}(x_n)p(xn)那么随机变量函数的数学期望就可以用上式计算。当然如果 xnx_nxn 是连续的随机变量那么就可以用积分形式来表示上面的表达式 E[f(X)]∫−∞∞f(x)p(x)dxE[{f}(X)] \mathop{\int}\limits_{-\infty}^{\infty} {f}(x) \mathrm{p}(x) dx E[f(X)]−∞∫∞f(x)p(x)dx 这时为了验证蒙特卡洛积分公式的正确性我们来看看它的数学期望 E(F(X))E[1N∑n1Nf(xn)p(xn)]1N∑n1N∫abf(x)p(x)×p(x)dx1N∑n1N∫abf(x)dx∫abf(x)dxE(F(X)) E \left[ \cfrac{1}{N} \sum \limits_{n1}^{N} \cfrac{{f}(x_n)}{\mathrm{p}(x_n)} \right] \\[2ex] \cfrac{1}{N} \sum \limits_{n1}^{N} \int \limits_{a}^{b} \cfrac{{f}(x)}{\mathrm{p}(x)} \times \mathrm{p}(x) \mathrm{d}x \\[2ex] \cfrac{1}{N} \sum \limits_{n1}^{N} \int \limits_{a}^{b} {f}(x) \mathrm{d} x \\[2ex] \int \limits_{a}^{b} {f}(x) \mathrm{d}x E(F(X))E[N1n1∑Np(xn)f(xn)]N1n1∑Na∫bp(x)f(x)×p(x)dxN1n1∑Na∫bf(x)dxa∫bf(x)dx 这个结果简直太棒了它说明蒙特卡洛积分公式的数学期望就是我们要求的积分值。这也间接的证明我们前面靠瞪眼法得到的公式本质上是正确的。或者说蒙特卡洛积分的结果对于积分真值是无偏差的也就是说是原积分的无偏估计。这个可以由“大数定理”来保证。只是不同的随机变量序列 xnx_nxn 会以不同的“速度”靠近积分原值。
四、误差 按照数学家们的“尿性”一定会严格的分析一下整个“蒙的挺准”的过程中的误差。既然刚才已经分析了蒙特卡洛积分的数学期望那么继续根据概率论的知识我们来分析下它的误差。一般对于随机序列来说我们取其方差来评估其靠近真实值的程度。那么就让我们来计算下蒙特卡洛积分的方差 σ2[Fn(X)]σ2[1N∑n1Nf(xn)p(xn)]1N2∑n1N∫ab(f(x)p(x)−E(Fn(X)))2×p(x)dx1N[∫ab(f(x)p(x))2×p(x)dx−E(Fx(X))2]1N[∫abf(x)2p(x)dx−E(Fn(X))2](5)\sigma^2[F_n(X)] \\[2ex] \sigma^2 \left[ \cfrac{1}{N} \sum \limits_{n1}^{N} \cfrac{f(x_n)}{\mathrm{p}(x_n)} \right] \\[2ex] \cfrac{1}{N^2} \sum \limits_{n1}^{N} \int \limits_a^b \left( \cfrac{f(x)}{\mathrm{p}(x)} - E(F_n(X)) \right)^2 \times \mathrm{p}(x) \mathrm{d}x \\[2ex] \cfrac{1}{N}\left[ \int \limits_a^b \left( \cfrac{{f}(x)}{\mathrm{p}(x)} \right)^2 \times \mathrm{p}(x) \mathrm{d}x - E(F_x(X))^2 \right] \\[2ex] \cfrac{1}{N} \left[ \int \limits_a^b \cfrac{{f}(x)^2}{\mathrm{p}(x)} \mathrm{d}x - E(F_n(X))^2 \right] \tag{5} σ2[Fn(X)]σ2[N1n1∑Np(xn)f(xn)]N21n1∑Na∫b(p(x)f(x)−E(Fn(X)))2×p(x)dxN1a∫b(p(x)f(x))2×p(x)dx−E(Fx(X))2N1a∫bp(x)f(x)2dx−E(Fn(X))2(5) 上面最终包含了 1N\cfrac{1}{N}N1 项的最终的一堆公式就是蒙特卡洛积分的方差根据概率论的定义其标准差就是方差的算术平方根。 按照最终的证明结果蒙特卡洛积分值与积分真实值之间的标准差误差有如下的形式
误差∝1N(6)误差 \varpropto \cfrac{1}{\sqrt{N}} \tag{6} 误差∝N1(6) 即最终误差正比于N的平方根的倒数。
五、蒙特卡洛积分计算与收敛速度 有了前面这些数学知识的加持那么我们就可以来看看用蒙特卡洛积分来计算积分时到底“速度”如何 其实看到这里大家应该有个疑惑我们费劲巴拉的推导了半天到底在图个啥 那么如果你看了前一篇文章3D数学系列之——从“蒙的挺准”到“蒙的真准”解密蒙特卡洛积分首先应该想到有了蒙特卡洛积分公式我们是不是可以开始靠“蒙”来计算积分了答案是肯定的要计算蒙特卡洛积分我们就需要产生一个随机变量的序列 xnx_nxn 并且需要知道其概率密度函数PDFp(xn)\mathrm{p}(x_n)p(xn) 有时候不需要知道 PDF 的表达式只需要知道每个随机变量对应的概率值即可 然后代入公式的右侧进行计算 ∫abf(x)dx≈1N∑n1Nf(xn)p(xn)(3)\int \limits_a^b {f}(x) \mathrm{d} x \approx \cfrac{1}{N} \sum \limits_{n1}^{N} \cfrac{{f}(x_n)}{\mathrm{p}(x_n)} \tag{3} a∫bf(x)dx≈N1n1∑Np(xn)f(xn)(3) 这样计算的好处是我们可以用大一些计算能力充足时或者小一些计算能力不足时的随机变量序列来进行计算。 当然这样计算的限制是根据对蒙特卡洛方法标准差的计算最终可以知道如果我们想提高计算精度一倍那么差不多就需要将N提高4倍 1412\cfrac{1}{\sqrt{4}}\cfrac{1}{2}4121 即误差减小一半。所以这就是平常所说的蒙特卡洛积分法收敛慢的根本原因。 例如在我们关注的图形的光照积分估算过程中这是个麻烦的问题比如原先需要每个点采样1000条光线近似计算最终光照效果要提高一倍精确度从而改进画质效果的话那么每个点就需要采集4000条光线进行积分。因此基本上现实中很少直接用蒙特卡洛积分的基本形式来求积分值的。 所以最终这个结论是有点令人灰心丧气的不过请继续往后看这么好的方法一定已经有很多人帮我们想了很多办法来解决问题的。
六、重要性采样 OK看到这里你还没有晕过去的话请为自己点个赞先 鉴于前面说的直接运用蒙特卡洛方法进行积分计算时因为固有偏差的问题而导致算法过程收敛过慢N的平方根倒数的量级所以在计算效率上是不尽人意的。那么既然是因为偏差导致的问题我们就继续从偏差来分析看看 σ2[Fn(X)]1N[∫abf(x)2p(x)dx−E(Fn(X))2]1N[∫abf(x)2p(x)dx−(∫abf(x)dx)2]\sigma^2[F_n(X)] \\[2ex] \cfrac{1}{N} \left[ \int \limits_a^b \cfrac{{f}(x)^2}{\mathrm{p}(x)} \mathrm{d}x - E(F_n(X))^2 \right] \\[2ex] \cfrac{1}{N} \left[ \int \limits_a^b \cfrac{{f}(x)^2}{\mathrm{p}(x)} \mathrm{d}x - \left( \int \limits_a^b {f}(x) \mathrm{d}x \right)^2 \right] σ2[Fn(X)]N1a∫bp(x)f(x)2dx−E(Fn(X))2N1a∫bp(x)f(x)2dx−a∫bf(x)dx2 这时如果我们想要标准差方差的开方最小那么就需要让上式结果等于0此时我们继续使用瞪眼法可以发现这需要下面的等式成立 p(x)∣f(x)∣∫abf(x)dx(7)\mathrm{p}(x) \cfrac{|{f}(x)|}{\int \limits_a^b {f}(x) \mathrm{d}x} \tag{7} p(x)a∫bf(x)dx∣f(x)∣(7) 在很多其他的资料中推得这个式子之后忽略了讲解它的意义。乍看上去貌似我们要知道我们需要使用的随机变量的概率密度就需要先知道最终积分的结果使得这个式子变成了无用的“鸡肋”。 其实不然这个式子恰恰给我们指明了一个非常明确的方向。 它说明如果我们要使得方差、标准差最小那么就要找到一个随机变量序列使得它的概率密度函数与被积函数 f(x){f}(x)f(x) 在“形状”上要高度一致这是个很重要的信息。为什么这样说呢其实仔细观察这个式子可以发现其中的定积分其实是个常数 ∫abf(x)dx\int \limits_a^b {f}(x) \mathrm{d}xa∫bf(x)dx 所以整个式子表达的意思是说 p(x)∣f(x)∣SA∣f(x)∣(8)\mathrm{p}(x) \cfrac{|{f}(x)|}{S} A|{f}(x)| \tag{8} p(x)S∣f(x)∣A∣f(x)∣(8) 即最终当概率密度函数其实就是被积函数的一个常数倍的时候整个方差最小也即标准差会非常小。 这是整个重要性采样的出发点也就是说我们可以依据被积函数的特征构造一个与之匹配的概率分布函数的形状然后生成特意构造的随机序列再返回到蒙特卡洛积分 ∫abf(x)dx≈1N∑n1Nf(xn)p(xn)\int \limits_a^b {f}(x) \mathrm{d} x \approx \cfrac{1}{N} \sum \limits_{n1}^{N} \cfrac{{f}(x_n)}{\mathrm{p}(x_n)}a∫bf(x)dx≈N1n1∑Np(xn)f(xn) 进行计算那么只需要非常少的 xnx_nxn 就可以得到非常精确的积分结果因为这样的序列计算后使得蒙特卡洛积分最终的方差、标准差都变得非常小甚至为0。 也就是说这样的随机序列使得蒙特卡洛积分能够快速收敛也就是说咱们之前说的普通的蒙特卡洛积分过程误差与 1N\cfrac{1}{\sqrt{N}}N1 成正比导致我们必须付出平方倍的计算量才能使的误差变小的规律可以被打破了因为这个因子的后面那堆复杂的表达式在我们特意的设计下几乎 0所以我们完全不用浪费过多的计算量来获取较精确的积分结果。 这个过程就被称为重要性采样
七、重要性采样方法和过程 那么最终如何进行上面所说的重要性采样呢其实核心就是构造使的蒙特卡洛积分方差和标准差最小的随机数序列 xnx_nxn 。 1、根据定义 ∫abp(x)dx\int \limits_a^b \mathrm{p}(x) \mathrm{d}xa∫bp(x)dx 就是概率密度函数 p(x)p(x)p(x) 对应的累积分布函数Cumulative Distribution Function简称做 CDF通常其正式定义如下 CDF(x)∫−∞xp(x~)dx~(9)\mathrm{CDF}(x) \int \limits_{-\infty}^{x} p( \tilde{x} ) \mathrm{d} \tilde{x} \tag{9} CDF(x)−∞∫xp(x~)dx~(9) 然后选择一个跟被积分的函数“相类似”的概率密度函数 PDF 来计算其 CDF 2、接着推导出 CDF 的反函数 CDF−1(x)\mathrm{CDF}^{-1}(x)CDF−1(x) 3、然后生成一个均匀分布的随机变量序列 un∈[0,1]u_n \in [0,1]un∈[0,1] ; 4、使用 unu_nun 代入 CDF 的反函数 CDF−1(x)\mathrm{CDF}^{-1} (x)CDF−1(x) 计算出随机序列 xnx_nxn ; 5、最后将 xnx_nxn 代入蒙特卡洛积分中计算函数的积分值 通常在很小规模的 xnx_nxn 序列下就可以计算出很高精度的积分值。 举例来说假设我们需要求一个形似 ∫0πAsin(x)dx\int \limits_0^{\pi} A \sin(x) \mathrm{d}x0∫πAsin(x)dx 的积分此时根据被积函数为正弦函数的特征构造 p(x)csin(x)x∈[0,π]p(x) c \sin(x) \quad x \in [0,\pi]p(x)csin(x)x∈[0,π] 然后我们计算它的 CDF ∫0πcsin(x)dxc[−cos(x)]0π2c1∵∫0πp(x)dx1⇒c12用x替换积分上限得到CDF的表达式有CDF(x)∫0x12sin(t)dt12[−cos(t)]0x12(1−cos(x))\int \limits_{0}^{\pi} c \sin(x) \mathrm{d}x c \Biggl[-\cos(x)\Biggr]_0^{\pi} 2 c 1 \\[2ex] \because \quad \int_0^{\pi} p(x) \mathrm{d}x 1 \\[2ex] \Rightarrow c \cfrac{1}{2} \\[2ex] 用 x 替换积分上限得到CDF 的表达式有 \\[2ex] CDF(x) \int\limits_0^x \cfrac{1}{2} \sin(t)\mathrm{d}t \cfrac{1}{2} \Biggl[ -\cos(t) \Biggr]_0^{x} \cfrac{1}{2}(1-\cos(x)) 0∫πcsin(x)dxc[−cos(x)]0π2c1∵∫0πp(x)dx1⇒c21用x替换积分上限得到CDF的表达式有CDF(x)0∫x21sin(t)dt21[−cos(t)]0x21(1−cos(x)) 有了 CDF 那么我们来计算其反函数 y12(1−cos(x))⇒xarccos(1−2y)y \cfrac{1}{2}(1-\cos(x)) \\[2ex] \Rightarrow \quad x \arccos(1-2y) y21(1−cos(x))⇒xarccos(1−2y) 接着我们生成一个均匀分布的随机数序列 unun∼[0,1]u_n \quad u_n \sim [0,1]unun∼[0,1] 代入上式中计算得到指定分布的随机序列 xnarccos(1−2un)x_n \arccos(1-2u_n)xnarccos(1−2un)最终用这个序列按照蒙特卡洛积分公式计算即可快速得到 ∫0πAsin(x)dx\int \limits_0^{\pi} A \sin(x) \mathrm{d}x0∫πAsin(x)dx 的积分值。 那么详细的计算过程和程序就先略过了在后续的 IBL 教程中将有更详细的应用讲解和实际工程中的代码。有兴趣的同学可以自己编写程序比较一下即先用黎曼和及库函数sin计算得到积分结果然后与库函数cos的直接计算结果进行比较最好能求出方差然后再用刚才推导出的 CDF−1CDF^{-1}CDF−1 函数计算的序列代入蒙特卡洛积分公式中计算积分结果同样与cos函数直接计算的结果求差或者求方差看看有什么区别
八、重要性采样的优缺点 根据前面的叙述重要性采样的全部目的就是为了缩小 xnx_nxn 的规模同时还能保证一定的精度其实就是之前标题所说的从 “蒙的挺准” 到 “蒙的更准”而且计算量还大大变小了这对于很多需要数值积分的计算来说是非常非常好的优化方法。另外对于高维的积分来说重要性采样的方法依然有效这样带来的好处就是原本复杂度可能是 O(nα)α∈Z,α⩾1O(n^{\alpha}) \quad \alpha \in Z^,\alpha \geqslant 1O(nα)α∈Z,α⩾1 的问题有可能直接变成了 O(n)O(n)O(n) 的问题这对于我们关注的图形光照计算领域来说简直就是福音因为我们需要处理的图形往往都是 n3n^3n3 维度的并且在复杂光照渲染算法中需要大量的 3D 积分所以有了蒙特卡洛积分和重要性采样方法就使得该过程的计算量几何级数的被减少这最终使得实时 PBR 渲染成为可能。 当然重要性采样也有缺点根据刚才描述的过程大家应该能想到第一如果被积函数 f(x){f}(x)f(x) 过于复杂时我们基本无法简单的模拟出一个形似的 p(x)p(x)p(x) 函数从而整个过程有可能就前功尽弃第二个方面就是说通过 CDF−1CDF^{-1}CDF−1 函数计算出的序列 xnx_nxn 有可能使得 f(x){f}(x)f(x) 的值过于集中于一些较大函数值得地方从而导致最终积分的结果与真值之间产生正偏差。 当然幸运的是这两条缺点在我们的 PBR 渲染中几乎碰不到也就没什么太大的影响。
