蒙特卡洛应用：RTX 光线追踪算法 ReSTIR 原理

接上文：蒙特卡洛方法
相关论文：《Spatiotemporal reservoir resampling for real-time ray tracing with dynamic direct lighting: 2020》

ReSTIR 算法全称为 GIalgorithm based on Spatiotemporal Reservoir Resampling，是一种在实时路径追踪渲染算法，用于提高间接光照的采样率，来达到降低路径追踪噪点的目的，它非常适用于计算 GI，目前 UE5.4 的 Lumen 就已实现了基于 ReSTIR 的 Gather……

同时，它也是目前 NVIDIA RTX 应用的最新光线追踪算法之一，未来可期

在阅读这一章之前，请确保你已经完全理解了蒙特卡洛方法中的所有内容

四、蒙特卡洛与 ReSTIR

4.1 重采样重要性采样(RIS) —— 蒙特卡洛积分求解器

ReSTIR，我们先从其中的 "IR" 开始介绍

将前面 3.4 节的 SIR 方法真正用于重要性采样，即作为采样分布映射获得服从目标 PDF 分布的样本，这种重要性采样就被称为 RIS —— 重采样重要性采样(Resampled Importance Sampling)。RIS 也是首次把重采样引入渲染的开山之作

RIS 的思路也很明确：使用 SIR 的方式获得 n 个样本，并将这 n 个样本代入蒙特卡洛积分方程 $(x)=\int_{\Omega} f(x) d x\approx\frac{1}{n} \sum_{i=1}^{n} \frac{f\left(x_{i}\right)}{p\left(x_{i}\right)}$ 进行求解

代入之后是什么样的呢？这里我们回到前面 SIR 的部分：

我们目标是从概率密度函数为 $g(x)$ 的分布中进行采样， $g(x)$ 为 Target PDF
从一个非常简单的 PDF： $p(x)$ 中取出了 m 个样本，记为 $X=x_1, x_2, \cdots, x_m$
而通过 SIR 方法进一步获得的样本 Y 中有 n 个样本，分别为 $_1, \cdots, y_n$

那么最终的蒙特卡洛积分结果就应如下

$\int_{\Omega} f(x) d x\approx I_{R I S}^{n, m}=\frac{1}{n} \sum_{i=1}^n \frac{f\left(y_i\right)}{g\left(y_i\right)}$

而 $y_1, \cdots, y_n$ 服从的，应该是 SIR PDF： $p^{\prime}(z)=\frac{g\left(x_z\right)}{\frac{1}{m} \sum_{i=1}^m w\left(x_i\right)}$ ：因此可以将 SIR PDF 作为重要性采样的 PDF，即

$I_{R I S}^{n, m}=\frac{1}{n} \sum_{i=1}^n \frac{f\left(y_i\right)}{p^{\prime}\left(y_i\right)}= \frac{1}{n} \sum_{i=1}^n\left(\frac{f\left(y_i\right)}{g\left(y_i\right)} \cdot\left(\frac{1}{m} \sum_{j=1}^m w\left(x_{ j}\right)\right)\right)$

到此，我们就得到了一个简易版本的蒙特卡洛积分求解器

事实上，对于上面的蒙特卡洛积分求解器其实是一个弱化版的求解器。可以发现，我们在执行 SIR 步骤的时候：Y 样本中每一个 $y_i$ ，其实都是从同一个 X 样本中选取的

理论上，可以对要选取的每一个 $y_i$ ，都重新生成一份 X 样本，并且每次生成的 X 样本的总数量，也不再是固定的 m，而是不同的 $m_i$ ，这样做的一个好处就是：更不容易出现相同的两个样本 $y_{i1}, y_{i2}$ ，且最终方差也可能会更小，此时 $I_{R I S}^{n, m}=\frac{1}{n} \sum_{i=1}^n \frac{f\left(y_i\right)}{p^{\prime}\left(y_i\right)}= \frac{1}{n} \sum_{i=1}^n\left(\frac{f\left(y_i\right)}{g\left(y_i\right)} \cdot\left(\frac{1}{m_i} \sum_{j=1}^{m_i} w\left(x_{i j}\right)\right)\right)$ ，其中 $w\left(x_{i j}\right)$ 就是指当获得 n 个最终样本的第 i 个样本时，SIR 流程中确定的 m 个样本组成的 X 中的第 j 个元素的权重

4.2 蓄水池采样(RS, Reservoir Sampling)

接下来就是介绍 ReSTIR，其中的 "Re"

在介绍蓄水池采样算法前：再举一个完全不需要复杂的数学推导，生活中的例子

在一个军训基地里，开学的第一天上午，来了 N = 2087 名来军训的学生，现在你作为教官，要从中选出 K = 100 名作为最终主方队的队员，并且由于你在鹅城当过县长，所以你主打的就是一个公平……，思来想去，你最终采用了最质朴的抽签的方式，即让这 2087 名学生依次抽签，抽中“幸运签”的学生被选中作为最终主方队的队员

但是紧接着问题就来了，这2087名学生并非全部要来这个军训基地里训练的学生，甚至连一半不到，这可不，第一天下午2点，又来了387名学生，下午4点，又来了566名学生……甚至第二天上午，还在源源不断的来新生

第二天上午的你，发现这个问题变得棘手了起来，到现在为止，学生总数 N = 3040，可是你昨天就已经从前2087名学生中选出了100名队员，事实上，你应该在确保不会再有新学生进来之后再抽取，现在你要保证公平只能有一种方案了：即前面兴师动众，让2000名学生同时抽签的结果直接废弃，等到后面不会再有新生了之后，再统一举办一次更浩大的“抽签仪式“

可是，真的只有这一种方案了嘛，确定要重新抽签嘛？前面的努力白费不说，后面再举行一次更浩大的全员抽签是否过于浪费场地和时间？

其实，以上的问题用数学描述的话，一句话就可以概括：从 n 个数据的池子中，等概率抽选 m 个样本，但是，n 的大小是动态的，可能后面会非常大且不可预知。请问你如何确保在 n 每次扩容后，都能从当前的数据池中等概率选出 m 个样本，要求总复杂度 O(n)

此题的一个非常完美的方案，即是：蓄水池采样(Reservoir Sampling)，它的流程如下（总样本数量 n 未知，你要等概率抽取 m 个样本）

遍历所有的 n 个数据，对于前 m 个数据，无脑全部加入“蓄水池”，此时蓄水池已满
从第 m+1 个数据开始，如果当前遍历到样本池中的第 i 个数据，则在 [1, i] 范围生成一个随机整数 $p$ ，如果 $p \leq m$ ，则将该数据与蓄水池中的第 $p$ 个数据替换，否则该数据丢弃
重复步骤②，直到遍历完所有的 n 个数据，如果 n 此时扩容，则继续执行步骤②……

这个算法即可确保在任意时刻处理完前 k 个数据后，蓄水池中的 m 个数据都是以 $\frac{m}{k}$ 的概率获得的

到此，如果你是前面题目中的教官，在了解了这个蓄水池采样算法后，大可不必放弃之前的抽样结果，更不需要在人全部到齐后再进行一次完全抽样，只需要对剩下新来的学生进行如上步骤②的操作即可：即假设当前学生是第 k 个到达军训基地的学生，你只需在 [1, k] 范围内生成一个随机数 $p$ ，如果 $p \leq k$ ，就让这个学生和主方队中的第 $p$ 个队员交换，原先的队员出队他进队，这样在任意时刻，当前的这最终主方队的100队员，一定都是当前总样本下最公平的选择

蓄水池采样策略，本质上就是一个“流式”采样策略

4.2.1 加权蓄水池采样(WRS, Weighted Reservoir Sampling)

可以发现，前面举得例子中，样本的均匀的，如果样本分布不均匀，就需要对其进行扩展，加权蓄水池采样应运而生，不用害怕公式有多复杂，只需要知道它的本质思路是和蓄水池采样一样的，就可以

扩展后的问题：从 n 个数据的池子 X 中，抽选 m 个样本，其中数据池 X 中的第 i 个样本 $x_i$ 的权重为 $w\left(x_i\right)$ ，除此之外，n 的大小是动态的，可能后面会非常大且不可预知。请问你如何确保在 n 每次扩容后，都能从当前的数据池中选出 m 个样本，且确保其满足重要性采样规则

加权蓄水池采样(WRS, Weighted Reservoir Sampling) 的步骤如下：当前蓄水池 Y 中已有 $t-1$ 个样本，那么蓄水池的下一个要被选择的样本 $x_t$ 按照如下规则选出，这个 $x_t$ 被称为 Reservoir

重新遍历总样本池 X 的所有数据，假设已经遍历了前 $k-1$ 个样本，当前备选样本为 $x_t$ ，当前权重和 $w_{\text {sum }}=\sum_{i=1}^{k-1} w\left(x_i\right)$
紧接着遍历总样本池 X 中的第 k 个样本，并在 [0, 1] 范围生成一个随机整数 $p$
如果 $p>\frac{w\left(x_k\right)}{w_{s u m}+w\left(x_k\right)}$ ，则无视当前样本
$p\leq\frac{w\left(x_k\right)}{w_{s u m}+w\left(x_k\right)}$ ，则将当前备选样本丢弃，选择当前样本替换它，即 $x_t=x_k$ ！
但无论替换以否，都更新权重和 $w_{\text {sum }}=w_{\text {sum }}+w\left(x_k\right)$
重复步骤②，直到遍历完所有的 n 个数据，最终的 Reservoir $x_t$ 加入池 y，继续重复步骤①……

用代码描述：对于当前的 Reservoir $x_t$ ，其定义如下

struct Reservoir
{
    float y;    // 当前备选样本 x_t，即 Reservoir
    float wsum; // 当前权重和 w_sum
    float M;    // 当前总样本池大小，也可以理解为当前遍历到了第几个样本，即上面的 k
    float W;    // 这个先不管，后面会给出介绍，注意，它并不是 w(x_k)！后者应用小写 w 表示
};

计算伪代码如下：

Spatiotemporal Reservoir Resampling (ReSTIR) - Theory and Basic Implementation | Shubham Sachdeva Blog

试求出当前考察的第 k 个样本，被作为最终 Reservoir 的概率：

很明显，当前第 k 个样本被最终选中的概率为第 k 个样本在遍历到它时直接被选中的概率乘以后面所有的样本都不被选中的概率

第 k 个样本在第一次遍历到它当时被选中的概率为 $p_1=\frac{w\left(x_k\right)}{\sum_{i=1}^k w\left(x_i\right)}$
第 m 个样本 $(k<m\leq n)$ 不被选中的概率 $1-p_2(m)=\frac{\sum_{i=1}^{m-1} w\left(x_i\right)}{\sum_{i=1}^m w\left(x_i\right)}$
因此后面所有的样本都不被选中的概率为

所以第 k 个样本被最终选中的概率

$\begin{aligned} p\left(x_t\right) & =p_1 \cdot \prod_{m=k+1}^{m \leq n} \frac{\sum_{i=1}^{m-1} w\left(x_i\right)}{\sum_{i=1}^m w\left(x_i\right)} \\ & =\frac{w\left(x_k\right)}{\sum_{i=1}^k w\left(x_i\right)} \cdot \frac{\sum_{i=1}^k w\left(x_i\right)}{\sum_{i=1}^{k+1} w\left(x_i\right)} \cdot \frac{\sum_{i=1}^{k+1} w\left(x_i\right)}{\sum_{i=1}^{k+2} w\left(x_i\right)} \cdots \frac{\sum_{i=1}^{n-1} w\left(x_i\right)}{\sum_{i=1}^nw\left(x_i\right)} \\ & = \frac{w\left(x_k\right)}{\sum_{i=1}^n w\left(x_i\right)} \end{aligned}$

这个结果是不是看上去很熟悉，没错，它正是前面 3.4 节的 SIR 中每一个 $x_z$ 在 X 中应被选中的概率！

4.2.2 蓄水池 Reservoir 合并

假设当前有两个总样本池，分别为 A 和 B，分别有 n 和 m 个总样本，都按照 WRS 的方式采样

其中 A 的某一次 Reservoir 选择

已经考察了全部 n 个样本，每个样本都是以 $p_A\left(x_i\right)=\frac{w\left(x_i\right)}{w_{\text {sum}_A}}$ 的概率被考察的（前面已证明）
最终选择的样本 $x_t$ 已经确定，积累的权重和 $w_{\text {sum}_A}=\sum_{i=1}^n w\left(x_i\right)$

其中 B 的某一次 Reservoir 选择：

已经考察了全部 m 个样本，每个样本都是以 $p_B\left(y_i\right)=\frac{w\left(y_i\right)}{w_{\text {sum}_B}}$ 的概率被考察的
最终选择的样本 $y_t$ 已经确定，积累的权重和 $w_{\text {sum}_B}=\sum_{i=1}^m w\left(y_i\right)$

现在如果想要把总样本池 AB 合并，得到一个有 n + m 个总样本的新样本池，是否可以直接将上述算法合并，即是否可以把 B 的最终选择看作是 A 的一个额外样本，以 $w_{\text {sum}_B}$ 为权重被样本池 A 考察

答案是可以的！对于 A 中的每一个考察过的样本 $x_i$ ，新样本的引入降低了它被选中的概率，只有当同时满足如下条件时， $x_i$ 才会被最终选中

$x_i$ 在 A 中被最终选中。这个概率 $p_A\left(x_i\right)=\frac{w\left(x_i\right)}{w_{\text {sum}_A}}$
B 考察的最终样本 $y_t$ 没有被选中，这个概率为 $1-p(B)=\frac{w_{\text {sum}_A}}{w_{\text {sum}_A}+w_{\text {sum}_B}}$

那么原来 A 中的每个样本 $x_i$ 被选中的概率就变成了

$\begin{aligned} p^{\prime}\left(x_i\right) & =p_A\left(x_i\right) \cdot\left(1-p_B\left(x_i\right)\right) \\ & =\frac{w\left(x_i\right)}{w_{\text {sum}_A}+w_{\text {sum}_B}} \end{aligned}$

类似的可以证明 B 中每个样本 $y_i$ 被选中的概率就变成了

$p^{\prime}\left(y_i\right)=\frac{w\left(y_i\right)}{w_{\text {sum}_A}+w_{\text {sum}_B}}$

也就是说，对于合并后的有 n + m 个总样本的新样本池，每个样本被选中的概率就为

$p\left(z_i\right)=\frac{w\left(z_i\right)}{\sum_{j=1}^{n+m} w\left(z_j\right)}$ ，和前面 4.2.1 推导一致

证毕：两个 Reservoir 可以按照如上思路进行合并操作

4.2.3 总结：WRS 的优点

其不用关心样本总数，样本池可以是动态的
样本选择操作相互独立，可以多线程并行计算最终要选出的样本
样本池之间可以合并，支持分治处理

4.3 结合 WRS 及 RIS

前面通过 4.1 重采样重要性采样(RIS) 得到了蒙特卡洛积分求解器：

$\int_{\Omega} f(x) d x\approx I_{R I S}^{n, m}=\frac{1}{n} \sum_{i=1}^n \frac{f\left(y_i\right)}{p^{\prime}\left(y_i\right)}= \frac{1}{n} \sum_{i=1}^n\left(\frac{f\left(y_i\right)}{g\left(y_i\right)} \cdot\left(\frac{1}{m} \sum_{j=1}^m w\left(x_{ j}\right)\right)\right)$

其中 SIR PDF： $p^{\prime}(z)=\frac{g\left(x_z\right)}{\frac{1}{m} \sum_{i=1}^m w\left(x_i\right)}$

对于其中的任意一个样本，就有

$\begin{aligned} I_{R I S}^{1, M} & =\frac{f(y)}{p^{\prime}(y)} =f(y) \cdot\left(\frac{1}{{g}(y)} \cdot \frac{1}{m} \sum_{j=1}^m w\left(x_j\right)\right) \end{aligned}$

再看 3.4 重要性重采样(SIR) 为了解决原分布 $g(x)$ 没有解析解或者太过于复杂的问题，先构造了一个简单的 proposal PDF： $p(x)$ ，并为其每一个样本分配了一个权重 $w\left(x_i\right)=\frac{g\left(x_i\right)}{p\left(x_i\right)}$

现结合 WRS，我们可以对这个 proposal PDF 进行如上 4.2.1 的操作，也就是从这个 proposal PDF 中获取 Reservoir $x_t$ ，以及对应的 $p(x_t)$

这个样本 $x_t$ 代入到前面的 $I_{R I S}^{1, M}$ 中，就有：

$\begin{aligned} I_{R I S}^{1, M} & =\frac{f(x_t)}{p^{\prime}(x_t)} =f(x_t) \cdot\left(\frac{1}{{g}(x_t)} \cdot \frac{1}{m} \sum_{j=1}^m w\left(x_j\right)\right) \end{aligned}$

对应伪代码如下，function shadowPixel 亦是在计算 $I_{R I S}^{1, M}$ ，伪代码中的 $\hat{p}_q\left(x\right)$ 就是原分布 $g(x)$

回到前面 WRS 对 Reservoir 结构的定义

struct Reservoir
{
    float y;    // 当前备选样本 x_t，即 Reservoir
    float wsum; // 当前权重和 w_sum
    float M;    // 当前总样本池大小，也可以理解为当前遍历到了第几个样本，即上面的 k
    float W;    // 这个先不管，后面会给出介绍，注意，它并不是 w(x_k)！后者应用小写 w 表示
};

就也能一一对上了

struct 中的 float y：对应样本 $x_t$
struct 中的 float wsum：对应 $\sum_{j=1}^m w\left(x_j\right)$
struct 中的 float M：对应 proposal PDF 中的总样本数 m

到此，也可理解 struct 中红色注释的部分到底是什么，W 即 $\frac{1}{g\left(x_t\right)} \cdot \frac{1}{m} \sum_{j=1}^m w\left(x_j\right)$

最后就能得到 $\begin{aligned} I_{R I S}^{1, M} & =\frac{f(x_t)}{p^{\prime}(x_t)} =f(x_t) \cdot W \end{aligned}$

是的没错，这里举得例子，正是通过 WRS 方法，从前面从 X 中继续选出一组新的样本 Y 的方式：即 3.4 重要性重采样(SIR) 中标红的那一句话 “然后从 X 中继续选出一组样本 Y，每个 X 中的某个样本被选中”，收回伏笔，X 中每个被选中的样本，正是每次的 Reservoir $x_t$

五、完备的实时光线追踪算法 ReSTIR 理论

根据前面 1~4 章尤其是第四章的内容，可以进行一个阶段性的总结

前面做的所有的一切，都是为了最终求出 $F(x)=\int_{\Omega} f(x) d x$ 这样的积分结果

但是可惜，计算机不会计算积分，因为我们只能尽可能暴力更多的样本 $x_1, \cdots, x_n$ ，代入到 $f(x)$ ，并加权求和，也就是事实上，我们只能计算离散型积分 $F(x) \approx \frac{1}{N} \sum_{i=0}^{i<N} \frac{f\left(x_i\right)}{p\left(x_i\right)}$
结束了嘛？不，这个积分并不是我们平时上课做习题时遇到的那种，可以用二维坐标系表示的简单积分，而是一种甚至无法用图像完全表示出来的复杂积分，它可能是路径的积分，也可能是三维球体的积分，分布也并不均匀，甚至会有多个参数，因此，在保证计算量不大的情况下为了得到更加精确的积分结果，样本的选择就至关重要 —— 在选择样本时，我们需要考虑它的概率分布 $p\left(x_i\right)=P D F$
到底该如何选出服从任意概率分布的样本呢？SIR 采样分布映射方式给出了我们思路：那就是先从一个简单的 proposal PDF 中抽出一组样本，并给予每个样本被选中的概率 $p(z\mid X)=\frac{w\left(x_z\right)}{\sum_{i=1}^m w\left(x_i\right)}, \ w(x)=\frac{{g}(x)}{p(x)}$
实际应用上，我们会使用 WPS 流式策略进行选择，每一次流式操作可以选出一个样本，它也被称为 Reservoir
不同的两个 Reservoir 可以进行合并，同时也可以非常灵活地组织考察样本的时机，把样本分布在空间上(Spartial)、时间上(Temporal)，甚至是多线程上进行考察，这也是 ReSTIR 中的 S(spartial) 和 T(temporal) 两个字母的由来

将以上的顺序反过来，就是最高效的求解积分结果的一个完整流程，它基于蒙特卡洛方法

5.1 结合 RIS 及 MIS

恭喜你成功通关到了这儿，但很可惜，这只是个开始……，类比到现在火热的游戏《黑神话·悟空》，你只不过是通关了铁头娃，仅此而已。好戏还在后头

在前面的 RIS 中，全部的 m 个样本都来源于同一个 proposal PDF，结合 MIS，我们其实可以让这 m 个样本来自于不同的 proposal PDF，当然，这些 proposal PDF 需要有相同的域，每个 proposal PDF 的权重为 $w_i(x)=w_{i, m i s}(x) \frac{\hat{p}(x)}{p_i(x)}$

不过这种情况下最后仍然会得到一个唯一的 SIR PDF，所以作为重要性采样的 PDF 还是只有一个，它在形式上并不是真正意义上的 MIS，其本质上只是 MIS for Proposals

其流程示意图如右下所示，如果只考虑同一个同一个 proposal PDF，也就是前面最基础的 RIS，就是左下图

Understanding The Math Behind ReStir DI · A Graphics Guy's Note

当然，也可以对于每个不同的 proposal PDF 都求出对应的 SIR PDF，这种纯正的 MIS 方法为 MIS for Samples

5.2 条件期望与重期望

条件概率公式：事件 A 和事件 B 同时发生的概率，等于在发生事件 A 的条件下事件 B 发生的概率乘以事件 A 发生的概率： $(A B)=P(A) \cdot P(B \mid A)$ ，其中

P(AB) 称为联合概率(joint probability)
P(A) 称为边缘概率(marginal probability)
P(B|A) 称为条件概率(conditional probabilty)

其中某个随机变量的边缘概率可以通过对联合概率在另一个随机变量域上积分得到，这相当于穷举了另一个随机变量所有取值的情况下该随机变量的概率并求和，设有两个相关的随机变量 $A=f(x), \ B=g(y)$ ，则有

$\begin{aligned} & P(A)=\int_{\Omega} P(A B) \mathrm{d} y \\ & P(B)=\int_{\Omega} P(A B) \mathrm{d} x \end{aligned}$

所以在已知联合概率的情况下，我们也可以求出单个随机变量的期望

$\begin{aligned} E(A) & =\int_{\Omega} A \cdot P(A) \mathrm{d} x \\ & =\int_{\Omega} A\left(\int_{\Omega} P(A B) \mathrm{d} y\right) \mathrm{d} x \\ & =\int_{\Omega} \int_{\Omega} A \cdot P(A B) \mathrm{d} y \mathrm{~d} x \end{aligned}$

条件期望：设 X 和 Y 都是离散随机变量，则在 Y = y 的情况下，X 的期望 $\mathrm{E}(X \mid Y=y)=\sum_{x \in \mathcal{X}} x \mathrm{P}(X=x \mid Y=y)$

但如果 X 是一个连续随机变量，满足 $X=g(x)$ ，Y 还是离散随机变量不变，条件期望就为 $E(g(x) \mid Y=y)=\int_{\Omega} g(x) p(x \mid y) \mathrm{d} x$

对于 Y，每确定一个 Y = y，都能得到一个条件期望值，在 y 不固定的情况下，就有 $h(y)=E(X \mid Y)=\sum_{x \in X} x P(X=x \mid Y)$ ，X 为连续随机变量同理

此时再对 $h(y)$ 求期望，就是在求任意 Y 取值的情况下， $E(X \mid Y)$ 的期望，就可得重期望公式

$\begin{aligned} E[E(X \mid Y)] & =\int_{\Omega}\left(\int_{\Omega} g(x) p(x \mid y) \mathrm{d} x\right) p(y) \mathrm{d} y \\ & =\int_{\Omega} \int_{\Omega} g(x) \underbrace{p(x \mid y) p(y)}_{p(x, y)} \mathrm{d} x \mathrm{~d} y \\ & =\int_{\Omega} \int_{\Omega} g(x) p(x, y) \mathrm{d} x \mathrm{~d} y=E(X) \end{aligned}$

通过这个重期望公式 $E(X)=E[E(X \mid Y)]$ ，可以得到一个更有用的推论

$E(X Y)=E[E(X Y \mid X)]=E[X \cdot E(Y \mid X)]$ 马上就能用得上

5.2.1 无偏性 RIS 证明

考虑 MIS 和 RIS 结合后的情况，对于 target PDF： ${g}(x)$ ，现使用了 m 个不同的 proposal PDF $p_i(x)$ 生成了 m 个样本 $X=x_1, x_2, \cdots, x_M$ 。之后从这组样本中最终选定的样本为 y，对应的 SIR PDF 为 $p(y)$

对 X 中的每一个样本，计算权重 $w_i(x)=\frac{g(x)}{p_i(x)}$ ，那么前面 4.3 的蒙特卡洛求解公式在只考虑唯一样本 y 的情况下，形式就是 $I_{R I S}^{1, M}=f(y) \cdot\left(\frac{1}{{g}(y)} \cdot \frac{1}{m} \sum_{j=1}^m w_j\left(x_j\right)\right)=f(y) \cdot W(y)$

不过为了引出后面的内容，现在还需要多一个考量 z：也就是 y 在 X 中的下标，通俗易懂的讲就是 y 是 X 样本组中的第几个（ $y=x_z$ ），同时使用一个更一般的随样本变化的权重值 $m(x)$ 代替上面的 $\frac{1}{m}$ ，此时

$I_{R I S}^{1, M}=f(x_z) \cdot\left(\frac{1}{g(x_z)} \cdot m(x_z) \sum_{j=1}^m w_j\left(x_j\right)\right)=f(y) \cdot W(X, z)$

现在尝试求出 SIR PDF $p(y)$ :

对于唯一样本 $X=x_1, x_2, \cdots, x_M$ ，它被选出的概率为 $p(X)=\prod_{i=1}^m p_i\left(x_i\right)$ ，很好理解，就是每个 $x_i$ 都被选中的概率全部乘起来

确定了唯一样本集 X 后，从 X 中选出第 z 个样本的概率为

$p(z \mid X)=\frac{w_z\left(x_z\right)}{\sum_{i=1}^m w_i\left(x_i\right)}, w_i(x)=\frac{g(x)}{p_i(x)}$

因此，同时确定 X 分布并刚好选中其第 z 个样本的概率就为 $\left[\prod_{i=1}^m p_i\left(x_i\right)\right] \frac{w_z\left(x_z\right)}{\sum_{i=1}^m w_i\left(x_i\right)}$

在样本 X 确认后，选中我们想要的样本 y 的概率就为 $p(y)=\sum_{i=1}^M p(y, z=i)$ ，这个公式可以理解为 y 可能出现在任意一个位置，既然这样出现在每个位置后被选中的概率加在一起就是最终概率

根据边缘概率的算法，可得

其中 $X^{i \rightarrow y}$ 表示 y 位于 X 中的第 i 个位置，即 $X^{i \rightarrow y}=x_1, \cdots, x_{i-1}, y, x_{i+1}, \cdots, x_m$ ，这个公式可以理解为当 y 出现在第 k 个位置时，其它位置随便是哪个样本都行，所以它其实是关于其它 m−1 个位置的边缘概率

回到前面的公式： $W(X, z)=\frac{1}{g(x_z)} \cdot m(x_z) \sum_{i=1}^m w_i\left(x_i\right)$ ，现尝试求出 $y=x_z$ 时它的期望

这个公式是离散求和，且项的数量只有 m 个，大可全部求出来直接求期望

选中 y 的情况下，z = i 的概率 $p((z=i) \mid y)=\frac{p(y, z=i)}{p(y)}$ ，由此可得 $W(X, z)$ 的期望

$\begin{aligned} E\left[W(X, z) \mid x_z=y\right] & =\sum_{i=1}^m W\left(X^{i \rightarrow y}, i\right) p((z=i) \mid y) \\ & =\frac{1}{p(y)} \sum_{i=1}^m W\left(X^{i \rightarrow y}, i\right) p((z=i) \mid y) \\ & =\frac{1}{p(y)} \sum_{i=1}^m \int \cdots \int W\left(X^{i \rightarrow y}, i\right) p\left(X^{i \rightarrow y}, i\right) \mathrm{d} x_1 \cdots \mathrm{d} x_m \end{aligned}$

代入前面确定 X 分布并刚好选中其第 z 个样本的概率 $p(X,z)=\left[\prod_{i=1}^m p_i\left(x_i\right)\right] \frac{w_z\left(x_z\right)}{\sum_{i=1}^m w_i\left(x_i\right)}$

可得：

$E\left[W(X, z) \mid x_z=y\right]=\frac{1}{p(y)} \sum_{i=1}^m \int \cdots \int W\left(X^{i \rightarrow y}, i\right) \frac{w_i\left(x_i\right)}{\sum_{i=j}^m w_j\left(x_j\right)}\left[\prod_{j=1}^m p_j\left(x_i\right)\right] \mathrm{d} x_1 \cdots \mathrm{d} x_m$

看到这么长的公式要不要慌，代入一开始的 $W(X, z)=\frac{1}{g(x_z)} \cdot m(x_z) \sum_{i=1}^m w_i\left(x_i\right)$ ，~~这样公式就更长了，不是~~，你就会发现这一切都可以化简，这里直接给出化简得最后结果：

$\begin{aligned} E\left[W(X, z) \mid x_z=y\right] =\frac{1}{p(y)} \sum_{i=1}^M m\left(x_i\right) \end{aligned}$

这个化简过程非常的“爽”，所以看到这里的你或许也可以试试自己尝试化简一下

还没结束！再回到前面 RIS 估计量： $I_{R I S}^{1, M}=f(x_z) \cdot\left(\frac{1}{g(x_z)} \cdot m(x_z) \sum_{j=1}^m w_j\left(x_j\right)\right)=f(y) \cdot W(X, z)$

求其期望

所以，当 $\sum_{i=1}^m m\left(x_i\right)=1$ 时，该估计量就是无偏的，对于不考虑 MIS 的情况， $m\left(x_i\right)=\frac{1}{M},\ p_i(x)=p(x)$ 也必然无偏，同时由于： $\begin{aligned} E\left[W(X, z) \mid x_z=y\right] =\frac{1}{p(y)} \sum_{i=1}^M m\left(x_i\right) \end{aligned}$ ，如果 $W(X, z)$ 收敛于 $\frac{1}{p(y)}$ ，结果也必然无偏，这样就可以用 $W(X, z)$ 来替代复杂的 SIR PDF 来参与计算

5.3 考虑了 MIS 的不同样本池 Reservoir 合并

前面介绍 Reservoir 合并时，假设的情况其实是理想的情况，而在实际 ReSTIR 计算时，还并没有那么简单

结合 WRS 及 RIS 的 Reservoir 合并实际情景如下：

两个样本池，分别为 A 和 B

样本 A 的 proposal PDF 为 $p_a(x)$ ，对应 target PDF 为 $g_a(x)$ ，已经考察了 $m_a$ 个样本，最终选择的样本为 $y_a$
邻近样本 B proposal PDF 为 $p_b(x)$ ，target PDF 为 $g_b(x)$ ，已经考察了 $m_b$ 个样本，最终选择的样本为 $y_b$ ，对应的结果权重为 $W\left(X_b, z_b\right)$ ，SIR PDF 为 $p'_b(x)$

一样我们想要把样本池 AB 合并，即把 B 的最终选择看作是 A 的一个额外样本，但若样本池 A 和 B 的 target PDF 大相径庭，甚至域也不相同的话（例如不同的光源），前面的合并方式 —— 直接使用 Reservoir 的 $w_{\text {sum}_B}$ 参与合并就可能会引入偏差

论文中使用了一种新的合并公式：在样本 A 处考察来自近邻的 $y_b$ 样本时使用的 target PDF 还是样本 A 的 target PDF： $g_a(x)$ ，而使用的 proposal PDF 则是样本 B 的 SIR PDF： $p'_b(x)$

此时考察 $y_b$ 的权重就变成了 $w_{y_b}(x) = \frac{g_a(x)}{p'_b(x)}$ ，由于在无偏的情况下， $W(X, z)$ 收敛于 SIR PDF 的倒数 $\frac{1}{p(y)}$ ，公式就可转变成 $w_{y_b}(x) = g_a(x)\cdot W\left(X_b, z_b\right)$

但是还没有结束！对于待考查样本 $y_b$ 而言，由于考察它使用的 proposal PDF 已经是经过 $m_b$ 轮逼近之后的，其更接近于样本 B 的 target PDF，因此它的质量会比单独的样本 $y_a$ 好上 $m_b$ 倍。所以权衡下来，它被选择的概率也应该提高 $m_b$ 倍，也就是 $w_{y_b}(x) = g_a(x)\cdot W\left(X_b, z_b\right)\cdot m_b$ ，其实本质上， $W\left(X_b, z_b\right)\cdot m_b$ 的操作，就对应着前面对 $w\left(x_i\right)=\frac{g\left(x_i\right)}{p\left(x_i\right)}$ 求和得到 $w_{\text {sum}_B}$ 的操作，只不过就是把原本 B 的 target PDF $g_b$ 替换成了 $g_a(x)$ 而已