Hessian-Hamiltonian MC Rendering
Hessian-Hamiltonian MC Rendering
Introduction
本论文提出一种Hessian-Hamiltonian MC Rendering算法，简称H2MC，该算法基于Metropolis Light Transport，引入了Hamiltonian力学的思路，将光路贡献和转移概率类比为重力和势能，很好的提高了MLT中的accept rate，意味着有更高的收敛效率，但本身因为需要计算光路的一阶导，以及二阶导（Hessian Matrix），计算量比较大，因此，适用于渲染复杂场景，比如caustics，多次反弹的glossy材质以及运动效果（时间维度的求导）。
这类场景可以认为是在一个高维场景下，多数能量分布在局部微小区域内，因此增大了采样的难度。如上图(a)在distant光源下的24K金环，caustics是由表面极高的glossy材质产生，这样，我们可得(b),红星周围两次反射的光路形成的一个切片，这里，x轴表示屏幕空间下的一个维度，
表示BRDF采样方向，如(c)所示。可见，大多数区域的贡献都是0。图(d)是常用的MCMC算法下，采用isotropic的mutation，采样的接受率不高。图(e)是用Manifold Exploration算法，只是低维空间下的一个切线。图(f)是论文中的算法，在某一个采样点构建一个anisotropic的分布，有效的在目标区域采样，尽管有些采样被拒，但接受率还是很高。(g)和(h)体现了本算法采样更好的遍历性(ergodicity)。(i)是高斯分布图，考虑了宽度，采样可以得到该高斯分布的解析解。
要达到这个结果，存在两个关键问题，第一要基于本地信息获取这种anisotropy，第二就是基于该信息采样。首先，通过automatic differentiation(AD)获取光路的导数，但导数(梯度)只能告诉我们增长最明显的方向，我们还需要二阶导数，Hessian matrix来告诉我们坐标之间的关联，这样，我们通过Hessian可以捕捉anisotropy。然后，我们采用了Hamiltonian力学，模拟一个粒子的运动，在当前位置，随机指定一个初始速度，同时该粒子因为重力，会受到高光（贡献）区域的吸引，在该场景下给定时间T，计算粒子的新位置，实现了采样的能力。这里并没有直接采用Hamiltonian MC，因为该算法需要数值积分来获取一个样本，这个是光追中无法接受的。本论文中，给定一个初始的动量，在一个高斯分布的二阶函数上采用Hamiltonian力学，得到符号高斯分布的最终位置，这样产生了标准的Metropolis-Hasting采样。传统的MLT是基于当前采样点的isotropic分布，而本论文中的高斯分布是anisotropic，根据当前位置的梯度和Hessian二阶导数，中心朝向高贡献方向。
这里不介绍MLT的相关内容，直接进入Hamiltonian MC部分。
Hessian Hamiltonian MC
Hamiltonian MC如上图所示，首先，从之前的path space slice，我们可以得到一个target function以及一个随机获取的当前的采样点
；接着，我们把这个target function翻转，这样，原本较高的值则变为较小的值，满足了重力的规则，这时给
一个初始的动量，通常符号高斯分布；这样，在给定时间T时，我们可以获取新的样本
，然后执行MLT算法。因为重力的作用，新的样本往往具有较高的贡献值，则意味着较高的接受率。
哈密顿力学需要满足的微分方程如下：
这里，x是位置，p是初始随机的动量，E是energy，是势能（重力）和动能的和，A是质量的倒数
：
这里，把公式
代入公式
，可得：
这样，我们定义了一个轨迹，描述位置
和动量
随时间
变化。为了讲哈密顿力学引入到MCMC，我们定义一个指数函数如下：
这样，
可以看作一个高斯分布，期望值为0，协方差矩阵为
。如此，通过
随机获取一个动量
，在固定时间T下，模拟哈密顿力学获取一个新的位置
，我们可以定义一个接受函数，实现MLT和哈密顿力学的结合：
这里，哈密顿力学满足三个重要特性：
- 时间可逆，我们可以从
作为输入，其结果会是
。
- 映射保证体积不变，从
映射到
时，position-momentum空间保持不变，这样，我们不需要进行Jacobian determinant。
- 能量守恒，势能和动能之间的互转，总能量不变，理论上接受率始终为1，完美了。
然而，还有一个问题，该微分方程没有解析解，需要迭代求解，这就不适用于光追了。于是，论文中提出了H2MC的方式来解决这个问题。其核心思想就是泰勒公式的二阶导的近似，话说，我有点不确定这样的结果还能保证无偏吗？
这里，对任意offset x进行如下的近似解：
如此，则
，将公式(3)合并为一个：
记
，该方程的解析解：
这里，
：
这里，
是AH的特征向量，而
中，
用AH的特征值
，
和
则变为AG和
在
上的投影，则一个光路x对应的N维空间的解：
尽管现在我们可以求得x(t)的解析解，避免了迭代求解的问题，但产生了另一个问题，MLT中，从
转移到
时，新的位置的梯度和Hessian通常会不同，因此需要计算
对应的
值。这里，在固定时间T下，对于
，
到
是一个线性变换，而
是已知的高斯分布，我们可以记作：
因为
的期望是0，协方差矩阵为
，我们可以很容易的计算在
点对应的期望和协方差：
这样，我们就可以用该pdf来计算从
到
的转移概率，从而通过公式(3)计算接受率。
另一个问题则是二阶泰勒公式导致的，这个约等于会在一些情况下带来不准确，于是，这里引入了一个期望为0，用户指定的方差
(isotropic)的高斯分布，这样，最终的高斯分布是：
我们可以根据该高斯分布计算对应的pdf值，这里，我们记
，表示从位置x到位置y：
这样，完成了论文中对Hamiltonian MC的升级 Hessian-Hamiltonian MC ，主要的变化是根据梯度和Hessian，获取了一个anisotropic的高斯分布（黄色部分），同时，为了准确性，过滤掉偏移过大的采样点，增加了一个isotropic的高斯分布（紫色区域），最终形成了我们最终采用的高斯分布（红色区域）。
这里，还有两个小问题，如何选择A和T，由公式(5)，我们发现，当我们让A等于协方差矩阵时，同时设置T等于
时，公式(6)可得，在新的位置
仍然保持相同的协方差A。同时，如果target function
符合高斯分布，则
对应的Hessian H的逆的负
应该就是
的协方差A，同时，A作为质量的倒数，所有的特征值必须为正，我们可得：
这样，A和H共享相同的特征向量，更方便应用。
上图是金环场景下传统的HMC和H2MC之间的对比，HMC具有较高的接受率，并且更好的遍历了整个轨迹，但需要迭代计算且计算量过大。相比而言，H2MC只需要较小的采样量更好的覆盖整个轨迹。
应用
该算法在应用中需要和现有的MLT结合，同样分为large step和small step两种，前者采用MMLT和PSSMLT，后者提供了两种mutation的方式：一个是
的small step，这里，用户设置一个roughness的阈值来区分specular和non-specular，对于specular，用random number来采样，而对于non-specular，则基于球面坐标来采样，这样的好处是避免后者在随机数生成时出现维度上的相关性，提供discrepancy。同时，对于NEE算法中的implicit ray做了优化，以及将时间作为一个维度来支持。最后，也支持lens perturbation。
上图是在
中优化随机数采样方式的改善(c)，以及在H2MC时的效果(d)。
Result
上图是渲染10分钟后的效果，明显，H2MC的效果更好，噪点更少。
上图体现了在不同区域高斯分布的anisotropic的能力
上图体现了在复杂场景，特别是物体弯曲比较多的场景下，采用Hessian的优势。
个人总结
为了理解这篇论文，花了一天的时间来推导哈密顿的那个微分方程，了解了欧拉求解最速降线，拉格朗日方程的过程以及两者之间的相似性，这是启发哈密顿力学的起源。这个过程收获很大。
其次，非常佩服Hamiltonian MC这个想法，巧妙的借助哈密顿力学，将动能和势能之间的转换类比到MC采样之间的状态转移，进而有效的提高采样的接受率。
最后，当然是本论文通过泰勒二阶导，以及解析解的线性关系，提出了H2MC来解决渲染问题，让人不仅感叹，竟然还可以这样来搞渲染。
尽管该论文对数学物理要求很高，但整个过程逻辑连续，一气呵成，再看看该论文的作者，除了一作，其他四位都至少是海贼王中大将级别的角色，阵容太强大了。