MCMC加速——Gaussian Processes to Speed up Hybrid Monte Carlo for Expensive Bayesian Integrals高斯过程加速HMC

高斯过程为昂贵的贝叶斯积分加速混合蒙特卡罗

0. 总结

混合蒙特卡洛 (HMC) 通常是计算难以解析的贝叶斯积分的首选方法。然而，这种方法的成功可能需要对（未归一化的）后验及其偏导数进行大量评估。在后验计算成本高昂的情况下，这可能导致 HMC 的计算负载不可接受。我建议对 HMC 所需的大多数计算使用后验（对数）的高斯过程模型。在这个方案中，即使 GP 模型有些不准确，也只需要偶尔评估实际的后验，以确保生成的样本具有完全期望的分布。该方法在 10 维问题上进行了演示，其中 200 次评估足以从后验生成 100 个大致独立的点。因此，所提出的方案允许对具有后验计算要求的模型进行贝叶斯处理，例如涉及计算机模拟的模型。
关键词：高斯过程； 马尔科夫链蒙特卡罗； 混合蒙特卡罗；非参数模型； 衍生过程； 正交; 设计。

1. 介绍

积分评估是贝叶斯推理的核心：wrt 积分。后验可以用于预测和对似然wrt进行积分。先验用于获得模型比较的边际似然。后者通常比前者更困难，因为先验通常相当广泛并且可能性非常高。在本文中，我们考虑前者，并将参考后验 wrt。我们希望将其积分为目标密度。

对于许多有趣的统计模型，直接分析积分是不可能的，所以经常使用马尔可夫链蒙特卡罗 (MCMC) 方法。简单形式的 MCMC 可能需要对目标密度进行大量评估才能产生可​​靠的结果。主要原因之一是目标密度往往表现出变量之间的强依赖关系，使得目标密度高的区域形成狭窄的扩展流形；这种流形被分解成几个不同的区域（多模态）的情况更加困难，这里不再具体讨论。

目标密度中变量之间的强依赖关系可能是大多数贝叶斯推理问题的典型特征，流形内的扩展方向对应于观测数据的可能解释族。积分wrt。与基于点估计的范式相比，这些不确定性是贝叶斯处理的主要优点之一。

变量之间的强相关性使得难以探索目标密度，同时保持在高概率区域内。在后验本身在计算上需要评估的情况下，这使得简单的 MCMC 方案不切实际。计算要求苛刻的后验出现在几种情况下：1）分层模型，其中较低层次的层次已被整合以减少参数数量，通常会导致其余参数的复杂后验，以及 2）使用计算机模拟作为整体的模型其规范的一部分，例如对于许多地球物理或气象应用来说是典型的。在这些情况下，对单个点的目标密度的评估本身就需要计算。在这种情况下，显然我们需要尽可能有效地使用从每次评估中获得的信息。

本文提出的新方法是对 Duane 等人的 Hybrid Monte Carlo (HMC) 方法的详细说明。 (1987)，使用高斯过程 (GP) 模型来指导搜索并避免对目标密度函数的过多引用。 O’Hagan (1991) 以前曾建议使用高斯过程模型来评估积分，称为 Bayes-Hermite 求积，另见优秀的 O’Hagan (1992)。在该方法中，引入了“重要性”密度，并使用高斯过程对被积函数乘以后验密度和重要性密度的比率进行建模。对于分析的易处理性，重要性密度必须是正常的，尽管肯尼迪 (1998) 放宽了这一要求。正如我们将看到的，当前的方法试图解决完全相同的问题，但方式却截然不同。 Duane 等人的混合蒙特卡罗 (HMC) 方法。 （1987）是一种相当普遍适用的方法，通过抑制（低效）随机游走行为，从高度相关的密度中更有效地采样。尽管 HMC 通常比依赖随机游走的采样方案快得多，但它仍可能需要对目标密度进行大量评估。在本文中，我提出了对 HMC 的扩展，它允许进一步减少目标评估的数量。关键思想是在大多数计算中使用目标密度的近似模型，并且仅偶尔需要评估真实目标；这是以这样一种方式完成的，即生成的样本将具有完全期望的分布；相关思想参见 Neal (1996) 的第 3.5 节和 Liu (2001) 的第 9 章。本文的其余部分结构如下。在第 2 节中讨论了寻求抑制随机游走的动机，并概述了混合蒙特卡罗方法。在第 3 节中，回顾了高斯过程模型的一些结果，包括对函数导数的观察值进行调节。新算法在第 4 节中给出，并在几个简单问题上展示了它的性能。结论和进一步讨论在第 5 节中给出。

2. 随机游走方法与混合蒙特卡罗

在高度相关的分布中，依赖随机游走的采样方法在探索目标密度方面效率低下。使用大小为 $s$ 的步长，使用随机游走移动一段距离 $L$ 需要预期的 $L^2/s^2$ 步，因为轨迹通常会自行折叠。这在图 1 中的一个简单的二元正态情况下进行了说明，用于非常常用的吉布斯抽样和 Metropolis 抽样方法。

由于典型的步长被限制在最窄方向的数量级以便获得良好的接受概率(严格来说，这仅适用于 $D>2$ 维； 在 $D=2$ 维中，接受率仅随提案宽度线性下降，因此最佳宽度是两个方向中较长的一个。 所以图 1 在这方面有点误导。)，随着相关强度的增加，随机游走的二次探索时间将导致性能下降。 Gibbs 采样算法的情况类似。这个问题最直接的解决方案是进行变量转换以消除强依赖关系，但对于贝叶斯推理中遇到的复杂分布，这可能非常困难。其他方法包括 Neal (1998b) 的有序过松弛方案和混合蒙特卡洛算法，可以看作是生成 Metropolis 提议的更精细的方法。

Neal (1993) 在统计设置中回顾了 HMC 算法。在 HMC 中引入了一个虚构的动力系统，其中变量 $\theta$ 被解释为位置，并由一组“动量”变量 $\varphi$ 增强。 Metropolis 提案是通过虚拟时间模拟动力系统提出的。由于物理系统中的动量，随机游走被抑制。物理系统的能量通过通常的二元性与概率相关：
因此，势能 $E_{pot}$ 与我们希望从中采样的目标分布 $p(\theta )$ 相关。 动能与动量的平方有关。 注意 $\theta$ 和$\varphi$ 因素的联合概率成为独立部分； 因此，我们可以通过从关节采样并忽略动量从 $p(\theta )$ 中获取样本。

在恒定总能量 $H$ 下，系统通过虚构时间的演化由 Hamilton 方程控制：

不幸的是，对于大多数目标密度，方程 (2) 是耦合的非线性微分方程（涉及对数后验参数的偏导数），我们无法以封闭形式求解。 相反，我们通过迭代一阶步骤来近似解决方案：

在实践中，我们使用稍微修改过的“跨越式”方案，将前半步作为动量，然后对位置和动量进行全步迭代，最后再对动量进行半步； 这样做是为了使近似可逆（Metropolis 步骤的详细平衡需要）。

Hybrid Monte Carlo 的完整算法变为以下步骤的迭代：

1. 动量变量的 Gibbs 样本，其密度为标准normal（等式 (1)），
2. 模拟物理系统 $L$ 步以生成提案,
3. 在 Metropolis 步骤中使用联合密度 $p(\theta ,\varphi )$ 接受或拒绝提议的状态。此过程从联合密度 $p(\theta ,\varphi )$生成样本，我们通过忽略动量从边缘 $p(\theta )$ 中获得所需样本。
   

• 图 1. 从相关双变量正态密度采样的示例； 绘制密度的一个标准偏差轮廓。 显示了三种不同采样算法的样本轨迹：左：Gibbs 采样的 20 次（一半）迭代，中：20 个 Metropolis 提案（被拒绝的提案用虚线表示），右：20 组 Hybrid Monte Carlo 的越级迭代。 这三个轨迹代表了相当的计算负担，尽管 HMC 除了密度值之外还需要评估导数。

图 1 右侧显示了 HMC 在抑制随机游走方面的成功。请注意，所谓的 Langevin Monte Carlo 方法等价于 $L=1$ 的 HMC，但这当然消除了对随机游走的理想抑制。

如果我们已经准确地求解了 Hamilton 方程，那么 Metropolis 步骤总是会接受，因为总能量是由动力学守恒的。然而，由于我们只是大致遵循动态，因此能量可能会发生变化，这可能会导致拒绝。但是请注意，尽管动态是近似的，但获得的样本具有完全所需的分布。

2.1 使用混合蒙特卡罗
在实践中使用 HMC 过程时，需要解决两个主要问题以确保算法运行良好。首先，我们需要为越级迭代选择步长 $\epsilon$。如果使用等式 (3) 中的单步大小（适当的步长可能在 θ 的不同坐标之间变化很大； 接近最优的单个步长可能不容易估计。），则通过监控结果接受率很容易做到这一点。其次，一旦确定了良好的步长，我们需要确定合适的模拟长度 L。这通常可以通过一些实验来完成，例如，通过监控参数的自协方差函数并增加 L 直到获得大致独立的样本。太短的轨迹会导致无法抑制随机游走，太长的轨迹会浪费计算。

一旦确定了 $\epsilon$ 和 $L$ 的合理值，就可以运行 HMC 算法以产生所需的样本。在图 1 中，二元正态有一个方向，其长度是垂直方向的 10 倍。使用简单的 Metropolis 方案，我们预计步长会受到最窄方向的限制，因此预计大约需要 102 步才能生成大致独立的样本。在此示例中，我们看到 HMC 能够使用 20 次评估（后验及其偏导数）遍历分布，从而提供改进的性能。然而，我们可能仍然需要大量的评估来获得足够的样本。本文的关键思想是我们可以使用目标模型而不是目标本身进行动力学模拟，从而大大节省了必要的评估次数。事实证明，高斯过程是一个非常方便的后验建模工具。

3. 高斯过程回归

我们可以使用无噪声 (在实践中，我添加了一点点噪声（通过向协方差矩阵添加少量的恒等式）来改善数值条件。) 高斯过程来指定平滑函数上的分布。 考虑一组 $n$ 对 $D$ 维协变量 $x$ 和因变量 $y$： { x ( c ) , y ( c ) ∣ c = 1 , . . . , n } \{x(c), y(c) | c = 1, . . . , n\} {x(c),y(c)∣c=1,...,n}。 我们指定因变量具有联合正态分布，其协方差矩阵取决于协变量。 例如：

其中 $x_d$ 是第 $d$ 个协变量，粗体字用于表示整个 $n$ 个观测值集的聚合。 注意，这里 ${\Sigma }_{pq}=Cov(y(p),y(q))=C(x(p),x(q))$，因变量之间的协方差，是协变量的函数 $C$。 这定义了具有平稳（高斯）协方差函数的零均值高斯过程，有关 GP 模型的更多背景信息，请参见 O’Hagan (1978) 和 Neal (1998a)。 协方差由超参数参数化：长度尺度 $w_d$ 和整体尺度 $w_0$。 显然，协方差函数还有许多其他可能的选择，但在这里我们将坚持等式 (4)。

观察到一些 y 变量的值后，我们有兴趣预测对应于协变量的新值 $x^{\ast }$ 的 $y^{\ast }$，它是通过根据观察值对关节进行调节而获得的，产生一个正态：

先验地，我们可能对超参数的特定值没有任何明确的偏好。 因此，原则上，我们应该对这些设置一个先验并整合出 Neal (1998a) 中的超参数。 为简单起见，我将简单地优化可能性 wrt的 超参数。 对数似然 $L$ 及其偏导数 Williams 和 Rasmussen (1996) 由下式给出：

可以使用标准数值工具（例如共轭梯度方法）方便地优化似然性（受局部最大值警告）。

3.1. 高斯过程导数

由于微分是线性算子，所以高斯过程的导数又是高斯过程。 例如，导数的均值等于预测均值的导数：

我们可以类似地找到导数分布的方差。 对于目前的目的同样重要的是，我们希望观察偏导数的值并以此为条件来推断目标函数。 函数值和导数之间以及导数之间的协方差由下式给出：

所以我们可以简单地写出一个扩展的协方差矩阵（在所有值和偏导数之间），并将等式（5）中的观察向量扩展为包括观察到的导数，从而做出如图 2 所示的预测。

• 图 2. 从具有高斯协方差函数（$w_0=1$ 和 $w_1=1$）的高斯过程中抽取的随机函数示例：左：随机抽取四个函数，中间：随机抽取四个函数，条件是 3 个用 + 号表示的观测值， 右：与中间图相同的观察结果，但现在四个随机函数也以这些点的导数为条件，由小切线段表示。 可以看出关于函数的更多信息如何限制函数的分布，以及函数的不确定性如何随 x 变化。

4. GPHMC 算法

所提出算法的关键思想是将 HMC 算法与势能的 GP 模型相结合，或者等效地是负对数目标密度； 形式上，我们用 $\theta$ 识别 $x$，用减去对数（非标准化）目标密度的依赖变量 $y$识别。 零均值 GP 用于对减去经验均值的观测值进行建模。 GPHMC 算法是通过在一个数字上评估目标函数来初始化的，比如 D 个随机位置（可能从先验采样）。 该算法现在分两个阶段进行：

探索阶段:
此阶段的目标是收集有关目标密度的信息。我们使用具有长轨迹的 HMC，L = 1000。我们将 HMC 的势能设置为 $E_{pot}=\mu -\sigma$，这导致动态模拟寻找高目标密度的区域和 GP 模型不确定形状的区域目标密度。沿着轨迹监测不确定性的值，如果不确定性达到阈值 $\sigma =3$，进一步的模拟似乎不太可能产生有意义的信息，因此停止模拟，并评估目标密度（及其偏导数） .通过这种方式，算法会寻找设计点的良好位置，如果少量评估就足够了，这一特征至关重要。然后根据 Metropolis 标准模拟下一个轨迹，或者从新点开始，或者从前一个轨迹的开始（以避免在非常低概率的点重新开始）。在每个轨迹之间，协方差函数的超参数通过 20 次似然共轭梯度优化迭代进行细化，等式 (6)包括新收集的信息4(当收集到合理数量的数据点后，会逐渐移除初始（随机）数据点，因为它们往往具有非常低的目标密度值，这使得它们对于推理不重要，并且可能导致数值不稳定 GP模型。)

采样阶段
此阶段的目标是使用前一阶段生成的 GP 模型从目标分布中生成所需的样本。 我们使用 $E_{pot}=\mu$ 和 $L=1000$ 的 HMC。轨迹的终点受 Metropolis 接受步骤的影响，需要对目标密度进行单一评估； 这确保了生成的样本将具有完全期望的分布，尽管目标密度的 GP 模型可能仍然不完善。 拒绝率表示 GP 模型捕获目标密度的程度； 太多的拒绝要求延长探索阶段。

应该注意的是，在新算法中，步长 $\epsilon$ 和轨迹长度 $L$ 的设置远没有原始 HMC 算法那么重要。步长设置得足够低，以至于长轨迹在探索阶段很早就被接受。更难的是 $L$ 参数，必须仔细调整以获得 HMC 中的最佳性能，现在可以简单地设置为一些相当大的值（例如 L = 1000），因为模拟不再参考实际目标密度，而是使用GP 模型，假设评估成本要低得多。该属性可能使新方案在实践中更易于使用。

所提出方法的计算开销（不包括目标函数及其导数的评估）由两个主要部分组成。首先，我们需要遵循 L 次迭代的动力学，它最初需要 $C(x,x)$ 的单次求逆，它在（保留）约束 $c=i(D+1)$ 的数量上是三次的，其中 $i$ 是迭代。在采样阶段，现在可以在时间 $O(Lc)$ 中跟踪整个动态；在探索阶段，也需要 GP 的不确定性，复杂度为 $O(L{c}^2)$。其次，优化（6）中的超参数 需要 $O(c^3)$。正如第 5 节中所讨论的，这种相当大的开销使得该方法最适合目标函数的计算成本非常高的情况。

4.1 Demo

在本节中，我将通过将其应用于两个已知密度的问题来证明所提出方法的可行性。
第一个demo是与 $\exp (-8(x_1^2+x_2^2-1{)}^2)$ 成比例的 $D=2$ 维目标密度，见图 3。

GPHMC 方法在两个阶段的每个阶段中运行 100 次迭代。 在采样阶段没有拒绝，证实了目标密度的模型相当准确，实际上在这个例子中 100 次探索性迭代似乎是过大的。 通过自协方差函数判断，生成的样本似乎基本上是独立的。 总之，我们从所需分布中生成了 100 个独立样本，评估了目标函数 200 次及其偏导数 100 次。

第二个示例在 D = 10 维中更具挑战性，其协方差矩阵的特征向量为 $(1,1,...,1{)}^T$，对应的特征值为 $1.0$，其他所有特征值为 $0.01$。因此，目标密度有 $9$ 个受限方向和一个宽 10 倍的方向，并且所有变量都是（正）相关的。使用标准密度可以确定生成的样本是否足够，并且不会使非参数过程变得更容易。还要注意，协方差函数的参数化不允许 GP 直接发现有一个长方向和 9 个受限方向（通过设置 $w$ 超参数），因为我选择了对角线而不是轴对齐的方向。此外，GP 使用平稳协方差函数来模拟目标密度的对数（它是二次的），这并不适合这种情况。当然，在任何实际情况下，我们都应该尝试使用反映我们对手头实际目标函数的信念的协方差函数；但是，我们通常不会对此有深入的了解，而且在实践中方便处理的协方差函数也是有限的。总而言之，我试图创建一个具有现实难度的示例，同时易于评估。

• 图 4. D = 10 维法线示例。 左上：程序生成的设计点（前10个是随机的N(0,1)）； 所有 10 个维度都绘制在彼此的顶部； 可以看到目标密度的相关性是如何逐渐恢复的，右上：随着发现高概率区域，负对数密度逐渐下降，左下：每个维度的 10 自相关函数表明生成的样本接近独立 ，右下：经验协方差矩阵的经验均值和特征值。

模拟结果如图 4 所示。两个阶段都允许进行 100 次迭代。在采样阶段，100 个建议样本中的 5 个被拒绝，这表明对数目标密度的 GP 近似值并不完美，但仍然足以满足高接受率。同样在此示例中，我们仅使用 200 次密度评估和 100 次偏导数评估从目标密度生成了 100 个大致独立的样本。经验均值在基于 100 个独立样本的预期容差范围内，并且 100 个样本的经验协方差的特征值似乎也很典型（当非正式地与相应的 Wishart 分布进行比较时）。所以有理由相信这个过程在这个例子中产生了正确的结果。

我们可以将此性能与更简单的 MCMC 方案进行比较。使用简单的 Metropolis 生成 100 个大致独立的样本将进行大约 10000 次评估（100 个样本乘以最宽与最受限方向的平方比）。具有导致大约 20% 拒绝的步长的 HMC 程序需要大约 L = 20（通过实验找到），从而导致 100 个大致独立的样本需要大约 2000 次评估。总之，这个例子表明 GPHMC 比 HMC 快大约一个数量级，而 HMC 本身比简单的 Metropolis 快大约一个数量级。任何任务获得的速度提升当然取决于目标密度下的依赖程度。

5. 结论

讨论和未来方向 已经提出了一个新版本的 HMC 算法，它适用于对计算上难以评估的分布进行积分。在这里介绍算法的方式，我们还需要对偏导数进行评估，这在某些情况下可能比较容易计算；例如，对于 GP 模型本身等式 (6) 就是这种情况，其中对似然性的评估是约束数量的三次方，但是一旦完成，每个偏导数只需要二次方时间。如果导数不容易计算，当然可以只使用函数值来运行算法。

一个 D = 10 维的简单示例显示了新方法的可行性。很自然地想知道该方法可以用于解决多大的问题。这取决于确定目标密度需要多少约束。即使是正态分布也需要 O(D2) 约束；由于我们正在进行非参数拟合，并且由于设计点选择的顺序性，我们可能无法接近全局最优设计，并且预计至少需要小倍数的 D2 约束。在 D = 10 的例子中，在探索阶段结束时，我们保留了 75 个（在 100 个中，因为我们逐渐删除了低密度的初始点），给出了 75×11 = 825 个约束。在现代笔记本电脑上，相应的 825×825 协方差矩阵的求逆大约需要一分钟。由于计算时间在约束的数量上成立方缩放，因此对于直接实现，总体计算时间预计会类似于 D6，不幸的是，这排除了远大于 D = 15 左右的大小。

有趣的未来工作应该将这种新提出的方法与 O’Hagan (1991) 的相关 Bayes-Hermite 求积法进行比较。特别是似乎值得研究是否可以处理更合适的对数密度或密度比的协方差函数。目前使用具有固定均值的平稳协方差似乎是不合适的，但在高密度区域（对于良好建模很重要的区域），这可能不是一个糟糕的模型。对这些想法的进一步理解和发展应该能够对计算复杂的模型进行贝叶斯处理，尤其是那些涉及计算机模拟的模型。