哈密顿蒙特卡洛(HMC)方法

哈密顿蒙特卡洛（ HMC ）采样方法

〖摘要〗快速给出下一个状态的提议值是 MCMC 方法的关键环节。对于状态有限的离散概率质量函数而言，可以采用随机游走的方式选择下一个状态的提议值，然后使用 Metropolis 更新步骤；但对于连续的概率密度函数而言， 随机游走方式显然不利于快速遍历状态空间。哈密顿蒙特卡洛方法利用 Hamilton 动力学的可逆性、能量守恒、体积保持等特性，为构造马氏链提供了一种快速生成提议状态的方法，该方法与 MCMC 中的 Metropolis 更新（或其他更新方法）步骤结合，可以快速生成给定概率分布的样本。

〖原文〗 Radford M. Neal (2011), MCMC Using Hamiltonian Dynamics, Handbook of Markov Chain Monte Carlo.

1 概述

马尔可夫链蒙特卡罗 (MCMC) 起源于 Metropolis 等人 的经典论文 (1953)。它被用于模拟理想化状态下分子系统的状态分布。不久之后，引入了另一种分子模拟方法( Alder 和 Wainwright，1959 年)，其中分子运动是确定性的，遵循牛顿运动定律，可被形式化为 哈密顿动力学 。为了找到气态材料的特性，这些方法是渐近等效的，因为即使在确定性模拟中，材料的每个局部区域都会受到来自其他区域的随机影响。尽管 MCMC 和 分子动力学 的应用领域有很大重叠，但这两种方法在接下来的几十年中保持着各自独立的存在 ( 参见 Frenkel 和 Smit，1996 )。 1987 年，Duane, Kennedy, Pendleton 和 Roweth 的一篇里程碑意义的论文 将 MCMC 和分子动力学方法统一在了一起，并称其方法为 “混合蒙特卡罗”，缩写为 “HMC”。该命名中，HMC 的缩写形式中所包含的 “哈密顿蒙特卡罗” 似乎更具体、更具描述性，因此本文将主要使用该命名。Duane 等人 并未将 HMC 应用于分子模拟，而是应用于量子色动力学的晶格场理论模拟。 HMC 的统计学应用始于本文作者将其应用于神经网络模型 ( Neal，1996a )。在对 MCMC 方法的综述中，作者还提供了一个关于 HMC 的统计教程 ( Neal，1993，第 5 章)。

HMC 在统计问题文献 （如 Ishwaran，1999；Schmidt，2009） 和面向统计的综述文献（例如 Liu，2001，第 9 章） 中还存在着其他应用场景，但至今（注：论文发表时间为 2011 年）仍然未被统计学家们重视 （或许也包括格场理论以外的物理学家）。本文的主体内容包括：

• 首先描述 哈密顿动力学。

  • 尽管该术语在物理学以外的领域可能会比较陌生，但 HMC 所需的哈密顿动力学特征非常基础，需要对哈密顿动力学的微分方程进行离散化才能适应计算机实现，而通常使用的 “蛙跳式” 方案容易理解。

• 在对哈密顿动力学的介绍之后，描述了使用它来构造一个 MCMC 的方法。

  • 第一步是根据待采样的概率分布来定义哈密顿函数。除了我们感兴趣的变量 （位置变量） 之外，还需要引入辅助的 “动量” 变量，这些辅助变量通常具有独立的高斯分布。
  • 将这些动量变量的简单更新与位置变量的 Metropolis 更新 交替进行，在 Metropolis 更新中，根据哈密顿动力学计算路径，并提议新状态 （使用 jumpfrog 方法实现）。
  • 以此方式提议的新状态，比简单随机游走提议方法对状态空间的探索更快 （避免随机游走的另一种方法是使用短路径，但只部分替换路径间的动量变量，后续的路径往往会朝着相同方向移动）。

• 在介绍了基本 HMC 方法之后，作者讨论了调整 “jumpfrog” 跳跃步长和步数的实际问题，以及 HMC 随维度扩展的一些理论结果。

• 然后作者介绍了 HMC 的一些变体。

• 通过查看路径开始和结束处的状态 “窗口”，可以提高许多问题的 HMC 接受概率。

• 对于许多统计问题，路径的近似计算（例如，使用数据的子集）可能是有益的。使用 “剪切路径” 可以使 HMC 的调整更容易，当然使用错误步长计算的路径也只需很少的计算时间。

• 最后，当存在多个孤立峰值时，“退火” 方法可能会很有用。

2 哈密顿动力学

哈密顿动力学 （也译为 哈密顿动力学） 可以提供有用的、具有直觉的物理解释。在二维空间中，动力学可以被可视化为冰球在一个高度随位置起伏的表面（ surface ）上做无摩擦滑动的过程。该系统的状态由冰球位置（ 向量 $q$ ）和冰球动量 （ 向量 $p$，质量乘以其速度）　组成。根据物理学原理，冰球的势能 $U(q)$ 与其当前高度成正比，其动能 $K(p)$ 则等于 $|p|^2 /(2m)$ （$m$ 为冰球质量）。冰球在该表面的水平部分以 $p/m$ 的恒定速度移动；如果遇到上升坡度，冰球动量使其保持运动，但动能减小、势能增大，一直到动能减小到零时停止；此后冰球将向低处滑动，此时动能增加而势能减少。

在哈密顿动力学的非物理 MCMC 应用中，上述案例中的 位置变量 将与我们感兴趣的随机变量对应，其势能则可以表示为负对数概率密度；而动量变量　（每个位置变量对应一个）　将被人工引入作为辅助变量。

上述解释可能有助于推动下面的讨论，但如果您发现其他情况，那么该动力学也可以被简单理解为一组微分方程组造成的结果。

2.1 哈密顿方程

哈密顿动力学在一个 $d$ 维的位置向量 $q$ 和一个 $d$ 维动量向量 $p$ 上实施，也就是说整个状态空间具有 $2d$ 个维度。整个系统可以用 $q$ 和 $p$ 的哈密顿函数 $H(q,p)$ 进行建模描述。

2.1.1 基础方程

哈密顿方程：位置 $q$ 和动量 $p$ 随时间 $t$ 的变化情况，由哈密顿函数关于两者的偏导数决定，即有：

$$\begin{align*} &\frac{d q_i}{d t}=\frac{\partial H}{\partial p_i}　\tag{2.1} \\ &\frac{d p_i}{d t}=-\frac{\partial H}{\partial q_i} \tag{2.2} \end{align*}$$

其中，$i=1, \ldots, d$，表示维度。

对于任意的持续时间段 $s$，上述公式可以用于定义从时刻 $t$ 到时刻 $t + s$ 的状态迁移映射 $T_s$ ( $H$ 和 $T_s$ 均不依赖于 $t$ )。

我们也可以换一种形式，将向量 $q$ 和 $p$ 合并为维度为 $2d$ 的新向量 $z =(q,p)$ ，并将哈密顿方程改写为：

$$\frac{d z}{d t}=J \nabla H(z) \tag{2.3}$$

其中 $\nabla H$ 是 $H$ 的梯度 （即　$[\nabla H]_k=\partial H / \partial z_k$ ），并且

$$J=\left[\begin{array}{cc} 0_{d \times d} & I_{d \times d} \\ -I_{d \times d} & 0_{d \times d} \tag{2.4} \end{array}\right]$$

是一个 $2d \times 2d$ 的矩阵，其象限用单位矩阵和零矩阵定义。

2.1.2 势能和动能的建模

在哈密顿蒙特卡洛方法（ HMC ）中，我们通常使用 势能 和 动能 来表达哈密顿函数，表示如下：

$$H(q, p)=U(q)+K(p) \tag{2.5}$$

其中　$U(q)$ 为势能，定义为待采样分布 $q$ 的负对数概率密度（加上任何方便的常数）。$K(p)$ 称为动能，定义为

$$K(p)=\frac{p^T M^{-1} p}{2} \tag{2.6}$$

这里的 $M$ 被称为质量矩阵（Mass Matrix），是一个对称的正定矩阵。实际使用中，该矩阵通常是对角矩阵，并且通常是单位矩阵的标量倍数。上述 $K(p)$ 的这种形式对应于均值为 $0$、协方差矩阵为 $M$ 的高斯分布的负对数概率密度 （加上某个常数）。

使用 $H$ 和 $K$ 的上述形式，哈密顿方程（公式 2.1 和 公式 2.2 ）可以改写成如下形式：

$$\begin{align*} \frac{d q_i}{d t} &=\left[M^{-1} p\right]_i \tag{2.7}\\ \frac{d p_i}{d t} &=-\frac{\partial U}{\partial q_i} \tag{2.8} \end{align*}$$

2.1.3 一维的情况

现在在一个维度上考虑一个简单示例　（此时 $q$ 和 $p$ 是标量并且不带下标），其中哈密顿函数定义如下：

$$\begin{align*} H(q, p)&=U(q)+K(p) \tag{2.9}\\ U(q)&=\frac{q^2}{2} \tag{2.10}\\ K(p)&=\frac{p^2}{2} \tag{2.11} \end{align*}$$

正如我们即将在 第 3.1 节 中看到的那样，这对应于 $q$ 服从均值为 $0$、方差为 $1$ 的高斯分布。由该哈密顿函数 ( 遵循 公式 2.7 和 公式 2.8 ) 得出的动力学方程为：

$$\begin{align*} \frac{d q}{d t} &=p \tag{2.12}\\ \frac{d p}{d t} &=-q \tag{2.13} \end{align*}$$

对于某两个常数 $r$ 和 $a$，解具有以下形式：

$$q(t)=r \cos (a+t), \quad p(t)=-r \sin (a+t) \tag{2.14}$$

因此，$T_s$ 是在 $(q,p)$ 平面上，围绕原点顺时针旋转 $s$ 弧度的一种映射。

注意：在更高维度上，哈密顿动力学通常不具有这种简单的周期形式，但此示例确实说明了一些重要性质，接下来我们将讨论它们。

2.2 哈密顿动力学的性质

哈密顿动力学的一些性质对于利用它来构造马尔可夫链蒙特卡罗更新至关重要。

性质 1 ：状态迁移的可逆性

哈密顿动力学具备可逆性，即从时刻 $t$ 的状态 $(q(t)，p(t))$ 到时刻 $t + s$ 的状态 $( q(t + s)，p(t + s) )$ 之间是一对一的，映射 $T_s$ 具有唯一的逆映射 $T_{-s}$。该逆可以通过对 公式 2.1 和 公式 2.2 中的时间偏导数取负号获得。

当哈密顿函数具有 公式 2.5 的可分解形式、且 $K(p)= K(-p)$ （ 例如：采用 公式 2.6 中的动能公式）时，逆映射 $T_{-s}$ 可以通过对 $p$ 取反、应用 $T_s$ 映射、对 $p$ 再次取反来获得。

对于 公式 2.9 的简单一维示例，$T_{-s}$ 则变成了按照 $s$ 弧度做的逆时针旋转，抵消了 $T_s$ 的顺时针旋转。

哈密顿动力学的可逆性意味着： 采用哈密顿动力学的 Metropolis 更新能够保持分布的不变性。马尔可夫链状态迁移的可逆性需要其中用于提出提议状态的动力学基础具备可逆性。

性质 2 ：哈密顿量守恒

哈密顿动力学的第二个特性是保持哈密顿函数值的不变性 ( 即能量守恒 )。

从 公式 2.1 和 公式 2.2 可以很容易地看出这一点，如下所示：

$$\frac{d H}{d t}=\sum_{i=1}^d\left[\frac{d q_i}{d t} \frac{\partial H}{\partial q_i}+\frac{d p_i}{d t} \frac{\partial H}{\partial p_i}\right]=\sum_{i=1}^d\left[\frac{\partial H}{\partial p_i} \frac{\partial H}{\partial q_i}-\frac{\partial H}{\partial q_i} \frac{\partial H}{\partial p_i}\right]=0 \tag{2.15}$$

如果使用 公式 2.9 的一维示例，其哈密顿函数值正好是 $(q,p)$ 距原点的平方距离的一半，而且 $(q,p)$ 的解 ( 公式 2.14 ) 与原点保持恒定距离，这说明一维情况下哈密顿量 $H$ 守恒。

哈密顿量 $H$ 守恒意味着：对于指定提议的 Metropolis 更新，理论上的接受概率应当为 $1$。但实际上我们只能使 $H$ 近似守恒，因此无法完全实现接受概率为 $1$ 的理想目标。

性质 3： 体积不变性

哈密顿动力学的第三个基本特性是：它能够保持 $T_s$ 映射前后 $(q，p)$ 空间体积的不变性 （被称为 Liouville 定理）。如果通过 $T_s$ 对 $(p,q)$ 空间中具有体积 $V$ 的某个区域 $R$ 中的点做映射，则该映射构成的镜像区域 $R^\prime$ 将具有相同的体积 $V$。

如果使用 公式 2.9 所示的一维情况，其状态解是一种旋转（公式 2.14 ），显然不会改变体积，这样的旋转同样也不会改变区域的形状。

但多维情况会比这复杂很多，大部分映射会改变形状但维持体积不变性。在多维度情况下，哈密顿动力学允许在一个维度上进行拉升，只要能够保证在其他方向上会被挤压以保持体积的不变性即可。

哈密顿动力学的体积不变性意味着：我们无需再考虑 Metropolis 更新的接受概率因体积变化而产生的任何变化，即其雅可比矩阵的行列式值为 $1$。这一点对于计算效率非常重要，因为当我们使用其他非哈密顿动力学提议新状态时，通常需要计算其对应的雅可比矩阵行列式，在大多数情况下这是不可行的。如果所采用的动力学映射具备体积不变性，则可以有效省去这些复杂计算。

---

哈密顿动力学体积不变性的证明

方法 1：证明矢量场散度为零

哈密顿动力学的体积保持能力可以采用多种方式得到证明，其中值得注意的第一种方法是： 公式 2.1 和 公式 2.2 定义的矢量场 散度（Divergence） 为零，通过下式可以看出：

$$\sum_{i=1}^d\left[\frac{\partial}{\partial q_i} \frac{d q_i}{d t}+\frac{\partial}{\partial p_i} \frac{d p_i}{d t}\right]=\sum_{i=1}^d\left[\frac{\partial}{\partial q_i} \frac{\partial H}{\partial p_i}-\frac{\partial}{\partial p_i} \frac{\partial H}{\partial q_i}\right]=\sum_{i=1}^d\left[\frac{\partial^2 H}{\partial q_i \partial p_i}-\frac{\partial^2 H}{\partial p_i \partial q_i}\right]=0 \tag{2.16}$$

Arnold （1989) 已经证明：具有零散度的矢量场具有体积保持能力。

方法 2：

在这里，我没有利用散度特性，而是非正式地证明，哈密顿动力学能够保持体积。我采用的方法是：假设体积保持等效于 $T_s$ 映射拥有值为 $1$ 的雅可比矩阵行列式。这个假设仅与该行列式的作用有关，与变换对积分和概率密度函数的影响无关。

$T_s$ ( 被视为 $z =(q，p)$ 的映射 ) 的雅可比矩阵记为 $B_s$，其形状为 $2d×2d$。通常 $B_s$ 取决于映射之前 $q$ 和 $p$ 的值。当 $B_s$ 为对角矩阵时，很容易看到：其对角线元素的绝对值是 $T_s$ 在各维度上拉伸或压缩区域形状的因子，这些因子的乘积（对于对角矩阵，就是 $\mid \operatorname{det}(B_s) \mid$ ）则是该区域映射前后的体积变化因子。我不会在这里证明一般性的结果，但请注意，如果要旋转所使用的坐标系，$B_s$ 不再是对角矩阵，但 $B_s$ 的行列式仍然具有不变性，因此仍将得到映射前后的体积变化因子。

让我们首先考虑一维哈密顿动力学的体积保持性 （即 $d = 1$ ），此时可以删除 $p$ 和 $q$ 的下标。当 $\delta$ 接近于 $0$ 时，可以将 $T_\delta$ 表示为：

$$T_\delta(q, p)=\left[\begin{array}{l} q \\ p \end{array}\right]+\delta\left[\begin{array}{l} d q / d t \\ d p / d t \end{array}\right]+\text { terms of order } \delta^2 \text { or higher. }\tag{2.17}$$

从 公式 2.1 和 公式 2.2 中可以获得时间导数，而雅可比矩阵可以写成：

$$B_\delta=\left[\begin{array}{cc} 1+\delta \frac{\partial^2 H}{\partial q \partial p} & \delta \frac{\partial^2 H}{\partial p^2} \\ -\delta \frac{\partial^2 H}{\partial q^2} & 1-\delta \frac{\partial^2 H}{\partial p \partial q} \end{array}\right]+\text { terms of order } \delta^2 \text { or higher. } \tag{2.18}$$

该雅可比矩阵的行列式为：

$$\operatorname{det}\left(B_\delta\right)=1+\delta \frac{\partial^2 H}{\partial q \partial p}-\delta \frac{\partial^2 H}{\partial p \partial q}+\text { terms of order } \delta^2 \text { or higher } \tag{2.19}$$

由于当 $x$ 接近于 $0$ 时，$\log(1 + x) \approx x$ ，因此除去 $\delta^2$ 阶或更高阶的项外 （我们稍后会看到它也正好为零），$\log \operatorname{det}(B_δ)=0$ 。现在针对一个不接近于 $0$ 的时间间隔 $s$ 来考虑 $\log \operatorname{det}(B_s)$ 。此时设 $δ= s/n$ （对于某个整数 $n$），我们可以将 $T_s$ 表示为 $n$ 次 $T_\delta$ 的组合（来自路径中的 $n$ 个点）， 因此 $\operatorname{det}(B_s)$ 是在这些点上做 $n$ 折 $\operatorname{det}(B_\delta)$ 的乘积。可以发现：

$$\begin{align*} \log \operatorname{det}\left(B_s\right) &=\sum_{i=1}^n \log \operatorname{det}\left(B_\delta\right) \\ &=\sum_{i=1}^n\left\{\text { terms of order } 1 / n^2 \text { or smaller }\right\} \\ &=\text { terms of order } 1 / n \text { or smaller. } \tag{2.20} \end{align*}$$

请注意， 公式 2.16 中的求和中，$B_\delta$ 值可能随 $i$ 变化，因为 $q$ 和 $p$ 会沿着产生 $T_s$ 的路径发生变化。但如果假设路径不是奇异的，$B_δ$ 的这种变化沿任何特定路径都是有界的。取极限 $n \rightarrow \infin$ 时，可以得出结论 $\log \operatorname{det}(B_s)= 0$ ，所以雅可比矩阵行列式 $\operatorname{det}(B_s)= 1$，因此 $T_s$ 保持了体积。

当维度 $d> 1$ 时，则采用相同的参数。雅可比矩阵将具有以下形式 （比较 公式 2.14），其中显示的每个块元素都是一个 $d×d$ 的子矩阵，行列索引为 $i$ 和 $j$:

$$B_\delta=\left[\begin{array}{cc} I+\delta\left[\frac{\partial^2 H}{\partial q_j \partial p_i}\right] & \delta\left[\frac{\partial^2 H}{\partial p_j \partial p_i}\right] \\ -\delta\left[\frac{\partial^2 H}{\partial q_j \partial q_i}\right] & I-\delta\left[\frac{\partial^2 H}{\partial p_j \partial q_i}\right] \end{array}\right]+\text { terms of order } \delta^2 \text { or higher. } \tag{2.21}$$

当 $d = 1$ 时，此矩阵的行列式将是一个 $\delta^2$ 阶或更高阶的加项，因为所有 $δ$ 阶的项都被消去了。上面讨论的其余部分无需更改即可适用。

---

性质 4：满足辛条件

体积保持也是哈密顿动力学符合辛条件的结果。令 $z =(q，p)$，并且像 公式 2.4 中一样定义 $J$ ，则辛条件为：

映射 $T_s$ 的雅可比矩阵 $B_s$ 满足：

$$B_s^T J^{-1} B_s=J^{-1} \tag{2.22}$$

上式隐含了体积保持，因为其左右两边的行列式相等： $\operatorname{det}\left(B_s^T\right) \operatorname{det}\left(J^{-1}\right) \operatorname{det}\left(B_s\right)=\operatorname{det}\left(J^{-1}\right)$ ，意味着 $\operatorname{det}\left(B_s\right)^2$ 为 $1$ 。

当 $d> 1$ 时，辛条件甚至强于体积保持。

可以将哈密顿动力学和辛条件概括为：$J$ 是任意一个 $J^T=-J$ 且 $\operatorname{det}(J) \neq 0$ 的矩阵。

哈密顿动力学最可贵之处在于：可逆性、体积不变性和辛条件都可以被精确维护，即便哈密顿动力学有时因为实践需要是近似的。

2.3 哈密顿方程的求解

对于计算机实现而言，哈密顿方程中 $p$ 和 $q$ 的求解必须通过时间离散化 ( 使用一些小步长的 $\varepsilon$ ) 来近似。从时间为零时的状态开始，迭代地（近似）计算 $\{1\varepsilon, 2 \varepsilon, 3 \varepsilon\}$ 等时刻的状态。

在讨论如何执行此操作时，我们将先假设哈密顿函数的形式为 $H(q, p)=U(q)+K(p)$，如 公式 2.4 所示。同时，尽管动能公式可以采取很多种形式，但此处仅假设 $K(p)=\frac{p^T M^{-1} p}{2}$ （ 如 公式 2.6 所示），同时假定 $M$ 为对角矩阵，对角线元素分别为 $m_1，...m_d$ ，此时动能可分解为：

$$K(p)=\frac{1}{2}\sum_{i=1}^d \frac{p_i^2}{ m_i} \tag{2.23}$$

2.3.1 欧拉方法

解微分方程组最著名的近似方法是欧拉方法。对于哈密顿方程，此方法对位置和动量的每个维度上的分量执行以下步骤 （索引为 $i = 1,...,d$ ）：

$$\begin{align*} &p_i(t+\varepsilon)=p_i(t)+\varepsilon \frac{d p_i}{d t}(t)=p_i(t)-\varepsilon \frac{\partial U}{\partial q_i}(q(t)) \tag{2.24} \\ &q_i(t+\varepsilon)=q_i(t)+\varepsilon \frac{d q_i}{d t}(t)=q_i(t)+\varepsilon \frac{p_i(t)}{m_i} \tag{2.25} \end{align*}$$

上面关于时间的导数是由 公式 2.7 和 公式 2.8 的哈密顿方程形式导出的。如果从给定的 $t = 0$ 时刻的值 $q_i(0)$ 和 $p_i(0)$ 开始 ，我们可以迭代上述步骤以获得路径时刻 $1\varepsilon ，2\varepsilon，3\varepsilon, \ldots$ 的位置和动量值。因此在 $τ/\varepsilon$ 步之后 （假设 $τ/\varepsilon$ 为整数），可以得到 $q(τ)$ 和 $p(τ)$ 的近似值。

图 1(a) 显示了用欧拉方法近似的 公式 2.9 的哈密顿动力学结果，从 $q(0)= 0$ 和 $p(0)= 1$ 开始，使用 $\varepsilon = 0.3$ 的步长执行了 $20$ 步 ( 即持续时间为 $τ= 0.3×20 = 6$ )。结果不是很好：欧拉方法产生的路径发散到了无穷大，但真实路径应当是一个圆。使用较小 $\varepsilon$ 值，并使用更多的步长，应该能够在 $τ= 6$ 时产生更为准确的结果，但虽然发散较慢，却并没有完全消除。

图 1： 当 $H(q,p) = q^2/2 +p^2/2$ 时，使用三种方法逼近哈密顿动力学的结果。初始状态为 $q = 0$，$p = 1$。对于 (a)、(b) 和（c），步长为 $\varepsilon = 0.3$，对于 (d)，步长为 $\varepsilon = 1.2$ 。每种方法都显示了 $20$ 步的模拟路径以及真实路径（灰色）。

2.3.2 欧拉方法的改进

通过稍微修改欧拉方法，可以获得更好的结果，如下所示：

$$\begin{align*} &p_i(t+\varepsilon)=p_i(t)-\varepsilon \frac{\partial U}{\partial q_i}(q(t)) \tag{2.26}\\ &q_i(t+\varepsilon)=q_i(t)+\varepsilon \frac{p_i(t+\varepsilon)}{m_i} \tag{2.27} \end{align*}$$

改进方法在为位置变量 $q_i$ 计算新值时，采用本次迭代得到的动量新值 $p_i$ 。当然，我们也可以先更新位置变量 $q_i$ ，然后用其新值来更新动量变量 $p_i$ 。

图 1(b) 显示了修改后步长为 $\varepsilon = 0.3$ 的新欧拉方法结果。尽管并不完美，但其产生的路径比之前的欧拉方法更接近真实路径，并且没有趋于无穷大的趋势。这种更好的性能与新方法对体积的精确保持有关，新方法有助于避免发散到无穷大或螺旋成原点（ 因为这些现象通常涉及到体积膨胀到无穷大或收缩到零 ）。

要看到：尽管时间的离散化非常有限，新欧拉方法依然可以精确地保持体积，请注意，利用 公式 2.26 从 $(q(t)，p(t))$ 到 $(q(t)，p(t +\varepsilon ))$ 的转换，和通过 公式 2.27 从 $(q(t)，p(t +\varepsilon))$ 到 $(q(t +\varepsilon)，p(t +\varepsilon))$ 的两种单步变换都是 “斜切” 变换，也就是说，其中只有部分变量（ 要么是 $p_i$ 要么是 $q_i$ ）会发生变化，而其他变量不发生变化。任何斜切变换都具备体积保持能力，因为其雅可比矩阵的行列式为 $1$ ( 行列式中唯一的非零项是对角元素，而斜切变换的对角线元素全为 $1$ )。

2.3.3 蛙跳法（ Leapfrog）

使用蛙跳法可以获得更好的效果，其工作方式如下：

$$\begin{align*} p_i(t+\varepsilon / 2) &=p_i(t)-(\varepsilon / 2) \frac{\partial U}{\partial q_i}(q(t)) \tag{2.28}\\ q_i(t+\varepsilon) &=q_i(t)+\varepsilon \frac{p_i(t+\varepsilon / 2)}{m_i} \tag{2.29}\\ p_i(t+\varepsilon) &=p_i(t+\varepsilon / 2)-(\varepsilon / 2) \frac{\partial U}{\partial q_i}(q(t+\varepsilon)) \tag{2.30} \end{align*}$$

我们从动量变量的半步 $p(t+\varepsilon/2)$ 开始，然后用动量的更新值对位置 $q(t+\varepsilon)$ 进行整步更新，最后使用位置变量的新值对动量进行剩余半步 $p(t+\varepsilon)$ 的更新。理论上，可以使用任何动能函数实现类似方案，只需用 $\partial K / \partial p_i$ 代替上面的 $p_i / m_i$ 即可。

当我们第二次将 公式 2.28 到 公式 2.30 应用到从 $t + 1\varepsilon$ 到 $t + 2\varepsilon$ 的状态迁移时，可以将第一次更新的后半步（ 从 $p_i(t+\varepsilon / 2)$ 到 $p_i(t+\varepsilon)$ ）与第二次更新的前半步 （ 从 $p_i(t+\varepsilon)$ 到 $p_i(t+\varepsilon+\varepsilon / 2)$ ）进行结合。蛙跳方法看上去与 公式 2.26 和 公式 2.27 中的欧拉改进法非常相似，不同之处在于，蛙跳在路径开始和末尾仅对动量执行了半步操作，并且计算出的动量值的时间标签偏移了 $\varepsilon / 2$。

因为 公式 2.28 到 公式 2.30 都是斜切变换，所以蛙跳法也可以精确地保持体积。由于其对称性，也可以通过取 $-p$、应用相同数量的蛙跳、再次取 $-p$ 的方式实现。

图 1(c) 显示了使用蛙跳方法的结果，步长为 $\varepsilon = 0.3$，在此图范围内与真实路径几乎没有区别。在 图 1(d) 中，展示了使用 $\varepsilon= 1.2$ 的蛙跳法结果（仍然是 $20$ 个步骤，可以看到四个周期，而不是一个周期）。使用较大的步长，可以清晰地看到逼近误差，但是路径仍然基本保持稳定（并且将无限期保持稳定）。只有当步长大小接近 $\varepsilon= 2$ 时，路径才会开始变得不稳定。

2.3.4 误差特性

我将简要讨论离散化动力学方程解法产生的误差在步长 $\varepsilon$ 趋近于零时的极限表现。 Leimkuhler and Reich(2004) 提供了更详细的讨论。

对于可用的方法，当 $\varepsilon$ 变为零时误差也应当变为零，因此，状态自身的误差上界将同样适用于状态的可微函数。例如，当状态 $(q，p)$ 的误差不超过阶数 $\varepsilon^2$ 时，其可微函数 $H(q，p)$ 的误差也不应大于阶数 $\varepsilon^2$ 。

在本文中，局部误差是执行一步之后（从时间 $t$ 到 $t +\varepsilon$）的状态误差；全局误差是模拟持续某个时间间隔 $s$ 后的累积误差，因此需要执行 $s/\varepsilon$ 步。如果局部误差的阶为 $\varepsilon^p$，则全局误差的阶将是 $\varepsilon^{p-1}$ ，也就是说：阶数为 $\varepsilon^p$ 的局部误差通过 $s/\varepsilon$ 步累积，会得到阶数为 $\varepsilon^{p-1}$ 的误差。如果固定 $\varepsilon$ 而增加路径累积时间 $s$ 的话，则误差通常会随着 $s$ 呈指数增加。有趣的是，用辛方法模拟哈密顿动力学时通常不会发生这种情况，如 图 1 所示。

欧拉方法及其改进方法的局部误差阶数为 $\varepsilon ^2$ 、全局误差为 $\varepsilon$ 。蛙跳法的局部误差阶为 $\varepsilon^3$ 、全局误差阶为 $\varepsilon^2$ （如 Leimkuhler 和 Reich，2004，第 4.3.3 节 所示）。这种区别来自于蛙跳法的可逆性：任何可逆方法的全局误差都必须具有 $\varepsilon$ 的偶数阶。

3 哈密顿动力学用于 MCMC（HMC）

使用哈密顿动力学从分布中采样需要将目标分布的概率密度函数转换为势能函数 $U(q)$ ，并引入新的 “动量” 变量，使之与感兴趣的目标随机变量一起构成哈密顿系统。在完成如此的哈密顿系统之后，就可以在生成马尔可夫链的每次迭代中利用哈密顿动力学实现 Metropolis 更新了（注：根据 MCMC 算法，在马尔可夫链的每次迭代由 “提议” 和 “接受/决绝” 两个环节组成，第二个环节带有一定的接受概率，因此也被称为 Metropolis 更新 ）：

• 首先对动量变量 $p$ 做随机采样
• 然后利用哈密顿动力学给出位置变量 $q$ 的新提议
• 最后进行传统的 Metropolis 更新

3.1 正则分布

待采样的目标分布可以通过统计力学中的 正则分布 ( Canonical Distributions ) 概念与势能函数建立关系。

3.1.1 什么是正则分布？

对于某个物理系统，其在状态 $\mathbf{x}$ 上的正则分布具有如下概率密度函数形式：

$$P(\mathbf{x} = x)=\frac{1}{Z} \exp \left(\frac{-E(\mathbf{x} = x)}{T}\right) \tag{2.31}$$

其中 $E(x)$ 为状态 $x$ 对应的能量函数，$T$ 为系统温度 （ 温度越高，系统能量越大，此处假设温度以玻尔兹曼常数为单位 ）， $Z$ 是为了保证函数求和或积分结果为 $1$ 而设置的归一化常数。当我们对概率密度函数 $P(x)$ 的某些特定分布形式感兴趣时，可以通过一个温度固定为常数（如 $T=1$）的正则分布得到，其方法也很简单，只需定义能量函数即可：令能量函数 $E(x)= -\log P(x)-\log Z$ ，其中 $Z$ 是任意合理的正常数。

3.1.2 哈密顿系统的建模

哈密顿系统可以被建模为： 满足 $H(q,p) = U(q) + K(p)$ 条件约束的正则分布。哈密顿动力学的能量（ $H$ ）守恒性质，意味着：理论上的哈密顿路径将在具有恒定概率密度的一个超曲面内移动。

据前所述，哈密顿函数（ $H(q,p)$ ）是 “位置变量” $q$ 和“动量变量” $p$ 所构成联合状态的能量函数，根据 公式 2.31，联合状态的正则分布可以定义如下：

$$P(q, p)=\frac{1}{Z} \exp \left(\frac{-H(q, p)}{T}\right) \tag{2.32}$$

当能量函数可分解为动能和势能时，即 $H(q，p)= U(q)+ K(p)$ 时，上式的联合概率密度可以进一步改写成分解形式：

$$P(q, p)=\frac{1}{Z} \exp \left(\frac{-U(q)}{T}\right) \exp \left(\frac{-K(p)}{T}\right) \tag{2.33}$$

可以看出，$q$ 和 $p$ 所对应的能量相互独立，两者也都服从正则分布，其其能量函数分别为 $U(q)$ 和 $K(p)$。我们用 $q$ 表示感兴趣的变量（位置变量），并引入同样数量的动量变量 $p$ 以使哈密顿动力学能够运行。

3.1.3 目标分布与哈密顿系统

在贝叶斯统计中，模型参数的后验分布通常是关注焦点，因此，这些参数在哈密顿动力学系统中将扮演位置变量 $q$ 的角色。根据 第 3.1.1 节 末，对于我们感兴趣的某种概率密度形式，可以通过构造势能函数的方式，将其转换为哈密顿系统所需的正则分布。根据该方法，后验分布可以表示为( 假设 $T = 1$ )：

$$U(q)=-\log [\pi(q) L(q \mid D)] \tag{2.34}$$

上式中 $\pi(q|D) = \pi(q) L(q|D)$ 代表后验分布的分子项，$\pi (q)$ 为先验分布，$L(q | D)$ 是给定数据集 $D$ 时的似然函数。

3.2 HMC 算法

现在，我们有了关于哈密顿动力学的基础知识。 HMC 只能用于从 $\mathbb{R}^d$ 上的连续分布进行采样，可以对其进行密度函数求值。目前，我暂时假设密度在任何地方都不为零（但这在 5.1 节中有所放松），同时还假设能够计算对数密度函数的偏导数。

3.2.1 基本思路

HMC 算法总体上是在从 公式 2.33 定义的正则分布中采样，只是其中 $q$ 对应感兴趣的目标分布，由势能函数 $U(q)$ 指定，同时需要扩展动量变量，并由动量函数 $K(p)$ 指定。由于直接从 $q$ 的正则分布中采样比较困难，所以可以考虑设计一个可处理的 $p$ 的分布并从中采样，然后利用哈密顿动力学的能量守恒、可逆、体积不变性等特性，间接得到 $p$ 的样本。

我们可以根据需要选择动量变量 $p$ 的分布，它应当独立于 $q$，并由动能函数 $K(p)$ 指定。 在实践中，HMC 的动能函数大多采用二次形式，如 公式 2.6 所示，这意味着 $p$ 是一个具有零均值的多元高斯分布。通常，$p$ 的各组分被指定为相互独立，且各自具有方差 $m_i$ 。产生该分布（设置 $T = 1$ ）的动能函数形式化为：

$$K(p)=\frac{1}{2}\sum_{i=1}^d \frac{p_i^2}{m_i} \tag{2.35}$$

我们将在 第 4 节 看到 $m_i$ 的选择对性能产生的影响。

3.2.2 HMC 算法的两个阶段

根据上一节中的基本思路，HMC 算法的每次迭代都应当有两个步骤：第一个只改变动量；第二个可能会改变位置和动量。这两个步骤中的任何一个都能够保持 $(q,p)$ 联合正则分布的不变性，因此其组合也能够保证联合正则分布的不变性（即分布不变、哈密尔顿量不变、体积不变）。

（1）第一步：随机抽取动量变量值

动量变量的新值随机抽取自多元高斯分布（也可以是其他形式），与位置变量的当前值无关。公式 2.35 所定义的动能函数中，$d$ 个动量分量之间相互独立，并且各自服从均值为 $0$ 、方差为 $m_i$ 的高斯分布。由于这一步中，并没有改变 $q$ 的值，并且 $p_i$ 是从给定 $q$ 时的条件分布中得到的（这里的巧妙之处在于：由于 $p$ 和 $q$ 相互独立，所以其实条件分布 $P(p|q)$ 等于边缘分布 $P(p)$ ），因此本步骤显然能够保持联合正则分布的不变性。

（2）第二步：哈密顿提议及 Metropolis 更新

本步骤执行 Metropolis 更新，首先利用哈密顿动力学给出新状态的提议值。具体做法是：

从当前状态 $(q,p)$ 开始，使用蛙跳方法（或其他一些能够保持体积不变性的可逆方法）对 $L$ 步的哈密顿动力学进行模拟，步长为 $\varepsilon$ 。$L$ 和 $\varepsilon$ 是算法参数，需要对其进行调优以获得良好的性能（如 第 4.2 节 所述）。对这个 $L$ 步路径末端的动量变量求反（即 $-p$ ），就能给出一个新状态的提议值 $(q^\star,p^\star)$ 。这个提议的新状态值按照以下接受概率被用于马尔可夫链的下一个状态：

$$\min \left[1, \exp \left(-H\left(q^*, p^*\right)+H(q, p)\right)\right]=\min \left[1, \exp \left(-U\left(q^*\right)+U(q)-K\left(p^*\right)+K(p)\right)\right] \tag{2.36}$$

如果提议的状态未被接受（即被拒绝），则下一个状态与当前状态相同（并在通过马尔可夫链的平均状态估计某个状态函数的期望时再次计数）。在路径末端对动量变量的求反使 Metropolis 提议对称，这是使上述接受概率有效所需要的。在实践中不需要进行这种求反，因为动能函数通常是偶函数 $K(p) = K(-p)$，并且动量将在下一次迭代的第一步使用之前会被替换。 （假设 HMC 更新是唯一执行的更新操作。）

如果我们将 HMC 视为从 $q$ 和 $p$ 的联合分布中采样，则第二步使用哈密顿动力学发现的新提议值，使 $(q,p)$ 的概率密度保持不变或几乎不变。也就是说，移动到具有不同概率密度的 $(q,p)$ 点仅通过 HMC 迭代中的第一步完成，其中动量 $p$ 被新值替换。幸运的是，$p$ 的这种替换可以大大改变 $(q,p)$ 的概率密度，因此移动到具有不同概率密度的点并不是问题。仅从位置变量 $q$ 的角度来看，$(q,p)$ 的哈密顿动力学可以产生具有较大概率密度差异的 $q$ 值（等效地，差异较大的势能值 $U(q)$ ）。

显然，动量变量的采样对于获得 $q$ 的正确分布至关重要。如果不重新进行采样，$H(q,p) = U(q) + K(p)$ 将保持不变；由于 $U(q)$ 和 $K(p)$ 是非负的，如果不对 $p$ 重新做采样，则势能 $U(q)$ 永远不会超过 $H(q,p)$ 的初始值。

用 R 语言编写的 HMC 算法实现单次迭代的函数见 图 2。其前两个参数是函数：U 和 grad_U，U返回给定 $q$ 值时的势能，grad_U 返回给定 $q$ 值时 $U$ 的偏导数向量。其他参数还包括：蛙跳的步长参数 epsilon、蛙跳的步数 L、路径开始的当前位置 current_q。动量变量在此函数中首先被采样，并在最后被丢弃，最终只返回下一个位置的提议值。动能被设定为最简单的形式，$K(p) = \sum p^2_i/2$（ 即所有 $m_i$ 都是 $1$ ）。在这段代码中，$p$ 和 $q$ 使用向量运算同时更新所有分量。该 HMC 实现可以从 我的网页 获得。网页中还由其他程序，提供了一些有助于实际使用的额外功能，以及 第 5 节 中即将说明的一些 HMC 变体。

图 2: 哈密顿蒙特卡洛算法

3.2.3 HMC 的 正则分布不变性

上面的 Metropolis 更新对于 $q$ 和 $p$ 的正则分布（ $T = 1$ ）是可逆的，这种条件也被称为 “细致平衡条件”，可以用如下非形式化的形式表述。

假设将 $(q，p)$ 空间划分为区域集合 $A_k$，每个区域都具有相同的体积 $V$ 。让 $A_k$ 相对于 $L$ 步蛙跳操作并附加一个负动量操作后的镜像为 $B_k$。由于蛙跳的可逆性，$B_k$ 也将对空间作出划分，并且由于蛙跳步骤能够保持了体积（取负值后），每一个 $B_k$ 也将具有体积 $V$ 。对于所有的 $i$ 和 $j$，如果能够满足下式，则细致平衡条件可以得到保证。

$$P\left(A_i\right) T\left(B_j \mid A_i\right)=P\left(B_j\right) T\left(A_i \mid B_j\right) \tag{2.37}$$

其中，$P$ 是正则分布下的概率，而 $T(X \mid Y)$ 是接受由当前状态 $Y$ 向状态 $X$ 转移的条件概率。显然，当 $i \neq j$，并且 $T\left(A_i \mid B_j\right)=T\left(B_j \mid A_i\right)=0$ 时，公式 2.37 将得到满足。

由于哈密顿函数几乎在任何地方都是连续的，因此随着区域 $A_k$ 和 $B_k$ 趋近于无限小，每个小区域的哈密顿函数值变得越来越恒定，其值在区域 $X$ 上为 $H_x$，且区域内的正则概率密度和转移概率也变得趋于恒定。此时，对于 $i=j=k$ 可以将 公式 2.37 进一步改写为：

$$\frac{V}{Z} \exp \left(-H_{A_k}\right) \min \left[1, \exp \left(-H_{B_k}+H_{A_k}\right)\right]=\frac{V}{Z} \exp \left(-H_{B_k}\right) \min \left[1, \exp \left(-H_{A_k}+H_{B_k}\right)\right] \tag{2.38}$$

很容易看出此式成立。

细致平衡条件意味着：

该 Metropolis 更新能够保持 $q$ 和 $p$ 变量正则分布的不变性。这可以从以下推导看出来。令 $R(X)$ 为小区域 $X$ 中 Metropolis 更新拒绝提议状态的拒绝概率，假设当前状态服从正则分布，则下一个状态在小区域 $B_k$ 中的概率等于如下两者之和：一是当前状态在 $B_k$ 且提议状态被拒绝的概率；二是当前状态在其他区域，其迁移至 $B_k$ 的提议被接受的概率。进而 下一个状态在区域 $B_k$ 中的概率可以被记为：

$$\begin{align*} P\left(B_k\right) R\left(B_k\right)+\sum_i P\left(A_i\right) T\left(B_k \mid A_i\right) &=P\left(B_k\right) R\left(B_k\right)+\sum_i P\left(B_k\right) T\left(A_i \mid B_k\right) \\ &=P\left(B_k\right) R\left(B_k\right)+P\left(B_k\right) \sum_i T\left(A_i \mid B_k\right) \\ &=P\left(B_k\right) R\left(B_k\right)+P\left(B_k\right)\left(1-R\left(B_k\right)\right) \\ &=P\left(B_k\right) \tag{2.35} \end{align*}$$

因此，HMC 中的 Metropolis 更新使正则分布保持了不变性。

由于动量变量采样、Metropolis 更新（通过哈密顿动力学提议）两个步骤都能够保持正则分布的不变性，因此 HMC 算法整体上也能够保持此特性。

3.2.4 HMC 的分布可遍历性

通常，HMC 算法具有 “可遍历性”。它不会被困在状态空间的某些子集中，并将渐近地收敛于其（唯一）不变分布。在 HMC 迭代中，可以为动量变量采样任何值，然后以任意方式影响位置变量。但如果路径中的蛙跳步骤导致某种状态函数的精确周期性时，可能会失去可遍历性。

例如，使用 公式 2.9、公式 2.10、公式 2.11 的简单哈密顿函数，精确解 ( 公式 2.12 ) 的周期为 $2\pi$。使用步数为 $L$ 、步长为 $\varepsilon$ 、且 $\varepsilon L = 2\pi$ 的蛙跳找到的近似路径，可能会重返相同位置坐标。具有这样 $L$ 和 $\varepsilon$ 值的 HMC 不具备遍历性。对于 $L$ 和 $\varepsilon$ 的邻近值，HMC 在理论上可能是遍历的，但要花很长时间才能移动到整个状态空间。

这个潜在的不可遍历性问题可以通过随机选择较小间隔的 $\varepsilon$ 、 $L$ 或同时两者来解决( Mackenzie，1989) )。每次都这样做可能是明智的选择。尽管在实际问题中，变量之间的交互通常可以防止发生任何确切的周期，但近似周期可能仍会大大降低 HMC 的速度。

3.3 HMC 演示及其优点

在这里，我将说明 HMC 的一些实际问题，并展示其与简单方法 （如 随机游走 Metropolis ） 相比具有更高效的采样潜力。我将用简单的高斯分布来做演示，以便能够与已知结果做比较。不过应当清楚，HMC 通常用于更复杂的分布。

3.3.1 二维问题中的路径

考虑从两个变量构成的二维高斯分布中采样，两者的均值为 $0$ 、标准差为 $1$ 、两者相关系数为 $0.95$，我们将其视为 “位置变量”。同时引入另外两个对应的呈高斯分布的 “动量变量”，其均值均为 $0$、标准差为 $1$、两者的相关系数为 $0$。此时的哈密顿函数被定义为：

$$H(q, p)=q^T \Sigma^{-1} q / 2+p^T p / 2, \quad \text { with } \Sigma=\left[\begin{array}{cc} 1 & 0.95 \\ 0.95 & 1 \end{array}\right] \tag{2.36}$$

图 3 显示了基于该哈密顿函数的路径，例如：可以用来在 HMC 中提议一个新状态，该路径使用 $L = 25$ 步数和 $\varepsilon = 0.25$ 步长来计算。由于整个状态空间是四维的，因此该图在左边两张子图中分别展示了位置变量和动量变量的坐标，而第三个子图则显示了每个蛙跳步骤之后的哈密顿函数。

图 3：二维高斯分布的路径，采用步数为 $25$ 、步长为 $0.25$ 的蛙跳模拟生成。绘制的椭圆表示均值的一倍标准差。初始状态为 $q = [− 1.50,−1.55]^T$ 和 $p = [− 1,1]^T$

注意：

此路径与随机游走并不相似，位置变量从左下角开始，系统地向右上方向移动，直到到达右上角为止，在到达该点后，开始朝相反方向运动。该运动的一致性是由动量变量引起的（想象一下冰球的运动）。而动量变量 $p$ 沿对角矩阵方向的投影只能缓慢变化，因为该方向上坡度较小，导致很多蛙跳步骤在对角线上的运动方向保持不变。当大尺度对角运动发生时，将产生小尺度震荡，会在与“峡谷”（由变量之间的高度相关性造成）交叉的方向上来回运动。

控制这些小的振荡会限制可用步长的范围。从 图 3 的右图可以看出，哈密顿函数的值也有波动 （如果模拟足够精确的话应当是常数） 。使用更大的步长，会导致振荡加剧。在临界步长时 ( 此例中为 $\varepsilon = 0.45$ )，路径开始变得不稳定，并且哈密顿函数的值会无限制地增长。但是只要步长小于此临界值，无论蛙跳步数有多少，哈密顿函数的误差总是有界的。我们并不能保证所有哈密顿量的误差都不增长，但对于许多比高斯族更复杂的分布来说确实可行。

可以看出，沿路径的哈密顿函数误差往往是正的，而不是负的。在此示例中，路径末端的误差为 $+0.41$，因此，如果将此路径用于 HMC 提议，则接受该端点作为下一个状态的概率将是 $exp(-0.41)= 0.66$。

3.3.2 从二维分布中进行采样

图 4 和 图 5 显示了使用 HMC 和简单随机游走 Metropolis 方法分别从双变量高斯中采样的结果，本例情况与上例相似，但位置变量之间的相关性更强，为 $0.98$。

在前一个节的例子中，HMC 使用的动能（定义了动量分布）为 $K(p)= p^Tp / 2$ 。HMC 方法的 $20$ 次迭代结果见 图 4 右图 （每次迭代步数为 $L=20$, 步长为 $\varepsilon=0.18$ ）。选择这些值是为了使路径长度 $\varepsilon L$ 足以移动到分布中较远的点，同时也不至于过大，以至于路径因自身加倍而浪费计算时间。这些路径的拒绝率为 $0.09$ 。

图 4: 随机游走 Metropolis 方法的 $20$ 次迭代（每次迭代更新 $20$ 次）和 Hamiltonian Monte Carlo 方法的 $20$ 次迭代（每条路径 $20$ 步蛙跳），作用于二维高斯分布，其边缘分布的标准差为 $1$，两个变量之间的相关性为 $0.98$。注意图中仅绘制了位置变量的坐标，椭圆表示均值的一倍标准差范围。

图 4 还显示了随机游走 Metropolis 方法的 $400$ 次迭代中每个 $20$ 次迭代的状态，采用以当前状态为均值、相关性为 $0$ 、标准差为 $0.181$ 的二维高斯提议分布。此示例提议的标准差为 $0.18$，与 HMC 提议的步长相等，因此，这些随机游走提议中的状态变化与 HMC 的单跳步骤具有可比性。这些随机游走提议的拒绝率为 $0.37$。

在 图 4 中可以看到，在 HMC 路径 ( 图 3 中所示 ) 期间，系统运动产生了比相应数量随机游走 Metropolis 迭代更大的状态变化。

图 5 说明了 $20×200$ 次随机游走 Metropolis 迭代和 $200$ 次 HMC 迭代的长运行时差异。

图 5: 从上面显示的 $20$ 次迭代开始的 $200$ 次迭代，仅绘制了其中第一个位置坐标。可见 HMC 的遍历性更好、效率更高。

3.3.3 避免随机游走的好处 – 快

如上所述，避免了随机游走行为是哈密顿蒙地卡罗的主要好处之一。

在此示例中：

（1）随机游走

• 迭代步长的限制：由于两个位置变量之间存在强相关，因此要使随机游走 Metropolis 的接受概率提高，所提议更新的幅度就应与最受约束的方向上存在的标准差相当 ( 在此示例中为 $0.14$，协方差矩阵的最小特征值的平方根 )。使用一次吉布斯采样扫描产生的变化将具有相似的幅度大小。

• 迭代次数的限制：达到几乎不依赖于当前状态的新状态所需的迭代次数，主要由较少约束方向上所需的探索时间决定。对于本例而言，该方向的标准差为 $1.41$，大约是约束最大方向上标准差的十倍。 因此，我们可能期望需要大约 $10$ 次随机游走的 Metropolis 迭代，这时所接受的提议将迁移至几乎独立的状态。但实际上所需的数量大约是其平方 – 大约 $100$ 次迭代接受一次提议 – 因为随机游走的 Metropolis 提议无法保持朝着同一方向移动。为证明这一点，请注意从某个初始状态开始 $n$ 次迭代的随机游走 Metropolis ，其位置方差将正比于迭代次数 $n$ ( 直到该方差与该状态的总体方差相当 )，因为位置是每次迭代当中大部分独立运动的总和。因此，移动量的标准差（ 给出了典型的移动量 ）正比于 $\sqrt{n}$。

（2）HMC

HMC 的蛙跳步长同样受最大约束方向的限制，但许多蛙跳步骤的运动方向基本一致。因此，在 $n$ 步之后移动的距离将倾向于与 $n$ 成比例，直到移动距离变得与分布的整体宽度相当。与随机游走运动相比， HMC 的优势将大致等于最小限制方向和最大限制方向上标准差之比 ( 此处约为 $10$ )。

因为避免随机游走非常有益，所以本示例中随机游走 Metropolis 提议的最佳标准差实际上远大于此处使用的值 $0.18$。 一个 $2.0$ 的提议标准差将给出非常低的接受概率 ( $0.06$)，但是在提议接受概率很低的情况下，这被远距离移动（至一个几乎独立的点）所补偿，从而产生一种与 HMC 一样高效的方法。但是，这种以小接受概率进行大变更的策略，只适用于分布仅在一个方向上受到严格限制的情况，

3.3.4 从 $100$ 维分布中进行采样

可以通过 $100$ 维多元高斯分布来更典型地说明 HMC 和随机游走 Metropolis 的行为，其中变量之间相互独立，且均值为 $0$，标准差为 $0.01,0.02, \ldots, 0.99, 1.00$。假设我们并不知道此分布的详细信息，因此采用与上述相同的简单的、旋转对称的动能函数 $K(p)= p^Tp/2$ 实现 HMC，同时另外使用随机游走 Metropolis 提议，其中对每个变量的更新是独立的，均具有相同的标准差。如以下在 第 4.1 节 中所述，这两种采样方法的性能对于旋转都是不变的，因此本示例说明了它们如何在协方差矩阵特征值的平方根为 $0.01, 0.02, \ldots ,0.99,1.00$ 的任何多元高斯分布上执行。

对于这个问题，位置坐标 $q_i$ 和相应的动量坐标 $p_i$ 都是独立的，因此用于模拟路径的蛙跳对于每个 $(q_i,p_i)$ 对都是独立运行的。但是，路径是否被接受取决于蛙跳离散化导致的哈密顿函数总误差，该总误差是每个 $(q_i,p_i)$ 对引起的误差之和。保持较小误差需要将蛙跳步长限制为一个近似等于标准差最小值 $0.01$ 的值，这意味着需要许多步数，蛙跳才能移动到与最大标准差 $1.00$ 相当的距离。

与此保持一致，我使用步数为 $L = 150$ 、每次迭代的步长 $\varepsilon$ 随机选择（ 均匀地抽取自 $(0.0104,0.0156)$， 即 $0.013±20％$ ）的 HMC 来模拟路径；同时使用的随机游走 Metropolis ，其提议标准差从 $(0.0176,0.0264)$ 均匀得出，即 $0.022±20％$。这两种方法都接近最佳设置。 HMC 的拒绝率为 $0.13$，随机游走 Metropolis 的拒绝率为 $0.75$。

图 6 显示了 $1000$ 次 HMC （右）和随机游走 Metropolis（左） 迭代的结果，随机游走方法的每次迭代使用 $150$ 次更新，一边每次迭代所需时间与 HMC 方法之间具备可比性。该图显示了具有最大标准差的最后一个变量。随机游走 Metropolis 的这些值的自相关性明显高于 HMC。图 7 显示了使用 HMC 和随机游走 Metropolis 运行获得的 $100$ 个变量中，每个变量的均值和标准差的估计值（ 估计值只是 $1000$ 次迭代的值的样本均值和标准差）。除里前面几个新变量（具有最小标准偏差）外，HMC 的平均估计误差大约是随机游走 Metropolis 平均估计误差的 $1/10$ ，HMC 的标准差估计也更好。

图 6：对于 $100$ 维示例，具有最大标准差的变量的值，来自随机游走 Metropolis 运行和 $L = 150$ 时的 HMC 运行。为了匹配计算时间，随机游走 Metropolis 每次迭代进行了 $150$ 次更新。

图 7：使用随机游走 Metropolis（左）和 HMC（右）对 $100$ 维变量的均值（上）和标准差（下）进行估计。 $100$ 个变量在水平轴上用该变量的真实标准差做了标注，而估计值体现在垂直轴上。

本例中的跳步步长随机化遵循 第 3.2 节 末尾讨论的建议。在这个例子中，不采用随机化步长（ 如固定为 $\varepsilon = 0.013$ ）实际上会导致问题：标准差接近 $0.31$ 或 $0.62$ 的变量变化很慢，因为 $150$ 次步长为 $0.013$ 的蛙跳，会产生近一个完整或半个周期，因此一个可接受的路径不会对变量的绝对值产生太大影响。

4 HMC 的实践理论

要从上一节中说明的 HMC 获得好处（包括避免随机游走），就需要对 $L$ 和 $\varepsilon$ 进行适当调整。我将在下面讨论 HMC 的调优，并展示如何通过利用变量规模及其相关性知识来提升性能。在简要讨论了单独使用 HMC 不足以满足需要时该怎么做之后，我将讨论 HMC 的另一个好处（与简单的 Metropolis 方法相比）：维数可扩展。

4.1 线性变换的效果

像所有 MCMC 方法一样，如果将要采样的变量通过乘以一些非奇异矩阵 $A$ 来转换，则 HMC 的性能可能会发生变化。但如果同时将相应的动量变量乘以 $(A^T)^{-1}$，则性能保持不变 （每次迭代的计算成本除外）。这些事实为 HMC 的运行提供了思路，特别是当我们对变量规模和相关性有所了解时，能够极大地帮助我们提高性能。

令新变量 $q \prime = Aq$ 。 $q \prime $ 的概率密度将由 $P \prime (q \prime )= P( A^{-1} q \prime )/ | \operatorname{det}(A)|$ 给出，其中 $P(q)$ 是 $q$ 的密度。如果该分布满足势能函数 $U(q)$ 的正则分布（请参阅 第 3.1 节 ），则可以得到 $q \prime $ 的分布为 $U \prime (q \prime )= U(A^{-1} q \prime )$ 的正则分布。 （因为 $|\operatorname{det}(A)|$ 是一个常数，我们不需要在势能中包含一个 $\log | \operatorname{det}(A)|$ 项）。

我们可以为相应的动量变量选择所需的分布，可以使用和以前一样动能函数。我们也可以选择将动量变量转换为 $p \prime = A^{-T} p$，并使用新的动能 $K \prime (p \prime )= K( A^T p \prime )$。如果我们使用二次动能 $K(p)= p^TM^{-1}p / 2$ （见 公式 2.6 ），则新动能为：

$$K^{\prime}\left(p^{\prime}\right)=\left(A^T p^{\prime}\right)^T M^{-1}\left(A^T p^{\prime}\right) / 2=\left(p^{\prime}\right)^T\left(A M^{-1} A^T\right) p^{\prime} / 2=\left(p^{\prime}\right)^T\left(M^{\prime}\right)^{-1} p^{\prime} / 2 \tag{2.37}$$

其中， $M^\prime = (AM^{-1}A^T)^{-1} = (A^{-1})^TMA^{-1}$ 。

如果我们使用以这种方式转换的动量变量，则新变量 $(q \prime ,p \prime )$ 的动力学本质上复制了 $(q,p)$ 的原始动力学，因此 HMC 的性能也将是一致的。要看到这一点，请注意，如果我们遵循 $(q \prime ,p \prime )$ 的哈密顿动力学，则原始变量的结果将如下所示（参见 公式 2.7 和 公式 2.8）：

$$\frac{d q}{d t}=A^{-1} \frac{d q^{\prime}}{d t}=A^{-1}\left(M^{\prime}\right)^{-1} p^{\prime}=A^{-1}\left(A M^{-1} A^T\right)\left(A^T\right)^{-1} p=M^{-1} p \tag{2.38}$$

$$\frac{d p}{d t}=A^T \frac{d p^{\prime}}{d t}=-A^T \nabla U^{\prime}\left(q^{\prime}\right)=-A^T\left(A^{-1}\right)^T \nabla U\left(A^{-1} q^{\prime}\right)=-\nabla U(q) \tag{2.39}$$

这与 $(q,p)$ 的哈密顿动力学相匹配。

如果 $A$ 是一个正交矩阵 （例如旋转矩阵），即 $A^{-1}=A^T$ ，则将 $q$ 和 $p$ 同时乘以 $A$ 并不能改变 HMC 的性能（因为 $(A^{T})^{-1}=A$ )。如果我们选择动量为旋转对称分布，且 $M = mI$ （ 即动量变量之间相互独立，每个变量都具有各自的方差），则这种正交变换将不会更改动能函数 ( 因此不会更改动量变量的分布 ) ，因为我们将有 $M^\prime =(A(mI)^{-1}A^T)^{-1}=mI$ 。

这样的旋转不变性也适用于随机游走 Metropolis 方法，在该方法中提议分布是旋转对称的（例如，高斯与协方差矩阵 $mI$ ）。相反，Gibbs 采样不是旋转不变的，也不是使用 Metropolis 算法依次更新每个变量的可选方案 ( Gibbs 提议每次仅更改一个变量 )。但 Gibbs 采样对于变量的扩展 ( 通过对角矩阵进行转换 ) 是不变的，这对于 HMC 或随机游走 Metropolis 是不正确的，除非以相应方式转换动能或提议分布。

假设我们有一个关于 $q$ 协方差矩阵的估计 $\Sigma$，并且还假设 $q$ 至少具有大致高斯分布。我们如何使用这些信息来改善 HMC 的性能？

一种方法是对变量进行变换，以使其协方差矩阵接近于单位矩阵，方法是查找 Cholesky 分解 $\sum = LL^T$，$L$ 为下三角，让 $q \prime = L^{-1} q$。然后，让动能函数为 $K(p)= p^Tp/2$。由于动量变量是独立的，并且位置变量接近于独立、方差接近于一（如果我们的估计 $\sum$ 和 $q$ 接近高斯的假设是好的），那么 HMC 应该在使用少量蛙跳步数的路径上表现良好，该路径会将所有变量移动导一个邻近的独立点。实际上，估计值 $\sum$ 可能不是很好，但是与使用原始 $q$ 变量使用相同的动能相比，此变换仍可以提高性能。使用估计协方差 $\sum$ 的等效方法是保持原始 $q$ 变量，但使用动能函数 $K(p)= p^T \sum p/ 2$ 。即让动量变量具有协方差 $\sum - 1$。可以通过将动能转化为与 $q \prime = L^{-1} q$ （参见 公式 4.1 ）的转换相对应的等价关系，得出 $K\left(p^{\prime}\right)=\left(p^{\prime}\right)^T M^{\prime-1}$ ，且 $M^{\prime}=\left(L^{-1}\left(L L^T\right)\left(L^{-1}\right)^T\right)^{-1}=I$。

使用这种动能函数来补偿位置变量之间的相关性在分子动力学中具有悠久的历史(Bennett，1975)。当维数高时，此技术的实用性受到矩阵运算的计算成本的限制。

即使在高维问题中，使用对角矩阵 $\sum$ 也是可行的。当然，这仅提供有关变量的不同尺度的信息，而不是它们的相关性。此外，当实际相关性为非零时，不清楚使用什么尺度。做出最佳选择可能是不可行的。在给定所有其他变量的情况下，对每个变量的条件标准差进行近似估计是可能的—正如我对贝叶斯神经网络模型所做的那样(Neal，1996a)。如果这也不可行，则对变量的边际标准差进行近似估计可能比对它们使用相同的比例更好。

4.2 HMC 的调优

使用 HMC 的一个实际障碍是需要为蛙跳步长 $\varepsilon$ 选择合适的值，以及蛙跳步数 $L$，它们共同决定模拟时间中路径的步长 $\varepsilon$ 和步数 $L$。大多数 MCMC 方法都需要调整参数，但当条件分布适合直接采样时，Gibbs 采样是一个例外。但是，在某些方面调整 HMC 比调整简单 Metropolis 方法更困难。

4.2.1 初步运行和路径图

调整 HMC 通常需要使用 $\varepsilon$ 和 $L$ 的早期试验值进行初步运行。在判断这些运行的效果时，应检查能够指示整体收敛情况的 路径图(Trace plot)。对于贝叶斯推断问题，高层次的超参数通常移动最慢。势能函数 $U(q)$ 的值通常也具有关键意义。这些量的自相关程度能够作为马尔可夫链在状态空间中探索程度的指示器。理想情况下，我们希望某一次 HMC 迭代后的状态几乎独立于先前状态。

不幸的是，如果初步运行的时间不足以达到平衡点，则可能会产生误导。达到平衡点前的最优 $\varepsilon$ 、 $L$ 值可能与到达平衡点后有所不同。甚至在达到平衡后，不同位置的最佳选择也不同。因此，如果有必要，在 HMC 的每次迭代中，可以从适合于状态空间不同区域的值中随机选择 $\varepsilon$ 和 $L$ 值（或按顺序使用这些值）。

建议使用不同的随机启动状态进行多次运行（初步运行和最终运行），以便检测出可能孤立的峰值。请注意，与其他对状态进行局部更新的 MCMC 方法相比，在孤立峰值问题上 HMC 并不会做的更好。当发现存在孤立峰值时，需要采取一些措施来解决它，仅仅将限制在单一峰值上的多次运行组合起来，通常是无效的。沿着路径（ 第 5.7 节 ）对 HMC 进行 “退火” 有时可以帮助处理多峰值问题。

4.2.2 步长的选择

选择合适的蛙跳步长 $\varepsilon$ 是至关重要的。步长太大，将导致通过模拟路径提出的状态的接受概率非常低。步长过小或与步长过小相同，都会浪费计算时间，或者如果路径长度 \varepsilon L 太短（则 $L$ 不大），则（更差）会导致随机游走导致探索缓慢。足够，请参阅下文。幸运的是，如图 3 所示，步长的选择几乎与完成多少跳越步长无关。哈密顿函数值（将决定拒绝率）的误差通常不会随着蛙跳步数的增加而增加，可以在一个简单的一维问题中看到稳定性问题，其中使用以下哈密顿函数：

$$H(q, p)=\frac{q^2}{2 \sigma^2}+\frac{p^2}{2} 、\tag{2.40}$$

此定义的 $q$ 分布是标准差为 $\sigma$ 的高斯分布。这个系统的一个蛙跳步骤（对于任何二次哈密顿量）将是一个从 $(q(t), p(t))$ 到 $(q(t+\varepsilon), p(t+\varepsilon))$ 的线性映射。参考 公式 2.28 至 公式 2.30，我们看到这个映射可以用矩阵乘法表示，如下所示：

$$\left[\begin{array}{c} q(t+\varepsilon) \\ p(t+\varepsilon) \end{array}\right]=\left[\begin{array}{cc} 1-\varepsilon^2 / 2 \sigma^2 & \varepsilon \\ -\varepsilon / \sigma^2+\varepsilon^3 / 4 \sigma^4 & 1-\varepsilon^2 / 2 \sigma^2 \end{array}\right]\left[\begin{array}{c} q(t) \\ p(t) \tag{2.41} \end{array}\right]$$

迭代此映射是产生稳定路径还是发散到无穷大的路径，取决于上述矩阵的特征值的大小，即

$$\left(1-\frac{\varepsilon^2}{2 \sigma^2}\right) \pm\left(\frac{\varepsilon}{\sigma}\right) \sqrt{\varepsilon^2 / 4 \sigma^2-1} \tag{2.42}$$

当 $\varepsilon /\varsigma> 2$ 时，这些特征值是实数，并且至少一个特征值的绝对值大于 $1$。因此，使用蛙跳方法计算出的路径将是不稳定的。当 $\varepsilon/\varsigma<2$ 时，特征值是复数，并且都具有下式的平方幅度

$$\left(1-\frac{\varepsilon^2}{2 \sigma^2}\right)^2+\left(\frac{\varepsilon^2}{\sigma^2}\right)\left(1-\frac{\varepsilon^2}{4 \sigma^2}\right)=1 \tag{2.43}$$

因此，用 $\varepsilon <2\varsigma$ 计算的路径是稳定的。

对于动能使用 $K(p) = p^Tp/2$（如上例）的多维问题，$\varepsilon$ 的稳定性极限将（大致）由最受约束方向上的分布宽度确定。对于高斯分布，这将是$q$ 的协方差矩阵最小特征值的平方根。 $K(p) = p^TM^{−1}p/2$ 的二次哈密顿量的稳定性可以通过应用使 $K(p \prime ) = (p \prime )^Tp \prime /2$ 的线性变换来确定，如上文 第 4.1 节 所述。

使用产生不稳定路径的，H 的增长幅度与 $L$ 呈指数关系，因此接受概率将非常小。对于低维问题，使用稍低于稳定性极限的 $\varepsilon$ 值足以产生良好的接受概率。但是，对于高维问题，可能需要进一步减小步长，以将误差保持在产生良好接受概率的水平。这将在 第 4.4 节 中进一步讨论。

当使用产生不稳定路径的步长 $\varepsilon$ 时，$H$ 的值随 $L$ 呈指数增长，因此接受概率将非常小。对于低维问题，使用略低于稳定性极限的 $\varepsilon$ 值足以产生良好的接受概率。然而，对于高维问题，步长可能需要进一步减小，以将 $H$ 中的误差保持在产生良好接受概率的水平。这将在 第 4.4 节 中进一步讨论。

选择过大的 $\varepsilon$ 值可能会对 HMC 的性能产生非常不利的影响。在这方面，HMC 比随机游走的 Metropolis 对调整更敏感。对于随机游走的 Metropolis ，需要选择一倍的提议标准差，性能会平稳下降；但是对于 HMC 而言，当 $\varepsilon$ 超过稳定性极限时，HMC 会产生急剧退化 （但在高维问题中，提议标准差过大的随机游走 Metropolis 的退化也可能非常明显，使得这种区别变得不太明显）。

当稳定性极限恒定时，过大步长导致的 HMC 性能下降不会是一个严重问题。该问题在初步运行中会很明显，因此大部分可以被修复。真正的危险是：状态空间中若干区域可能具有不同的稳定性极限，并且都具有较大概率。如果在稳定性极限较大的区域开始初步运行，则选择比该极限稍小的 $\varepsilon$ 似乎是合适的。然而，如果这个 $\varepsilon$ 高于某个其他区域的稳定性极限，则运行可能永远不会访问该区域，即使其概率很大，这也会产生非常错误的结果。要了解为什么会发生这种情况，请注意，如果运行确实访问了所选 $\varepsilon$ 会产生不稳定性的区域，它将在那里停留很长时间，因为该 $\varepsilon$ 的接受概率将非常小。由于该方法仍然保留正确的分布不变性，因此运行很少从所选 $\varepsilon$ 导致稳定路径的区域移动到该区域。可能出现此问题的一个简单情况是从具有非常轻的尾部（比高斯分布更轻）的分布中采样时，密度的对数将比二次方下降得更快。在尾部，对数密度的梯度会很大，为了稳定需要小步长。参见 Roberts 和 Tweedie (1996) 在 Langevin 方法的背景下对此的讨论（参见 第 5.2 节 ）。

当步长太大时，HMC 的性能急剧下降不会如果稳定极限是恒定的，那么这将是一个严重的问题-从初步运行来看，该问题将很明显，因此可以解决。真正的危险是稳定性极限可能会有所不同对于状态空间的几个区域，它们都有很大的概率。如果在稳定性极限较大的区域开始初步运行，则选择比该极限小的一点似乎是合适的。但是，如果这超出了其他区域的稳定性极限，则运行可能永远不会访问该区域，即使它具有很大的可能性，也会产生严重错误的结果。要了解为什么会发生这种情况，请注意，如果跑步者确实访问了选定的将产生不稳定状态的区域，它将在该区域停留很长时间，因为接受该区域的可能性很小。由于该方法仍然使正确的分布不变，因此只能从选定的 $\varepsilon$ 指向稳定路径的区域很少运行到该区域。可能会出现此问题的一个简单情况是，从具有很多轻尾巴的分布（比高斯亮）进行采样分布），其徽标的密度下降速度将比平方下降快。在尾部，对数密度的梯度会很大，并且需要较小的步长才能保持稳定性。参见 Roberts 和 Tweedie(1996)中关于 Langevin 方法的讨论（参见 5.2 节）。

通过从某些分布中随机选择，可以缓解此问题。即使此分布的主成分太大，也可能偶尔为 $\varepsilon$ 选择较小的值。（请参见 第 3.2 节 ，以随机改变步长）。$\varepsilon$ 的随机选择应在路径开始时进行一次，而不是每次蛙跳步骤，因为即使所有选择都低于稳定性极限，每个步骤的随机变化也会导致 $H$ 的随机游走，而不是 图 3 中所示的有界误差。 第 5.6 节 中描述的 “捷径” 过程当随机选择的步长不合适时，可以将其视为节省计算时间的方法。

4.2.3 路径步数的选择

如果 HMC 要系统地探索状态空间而不是通过随机游走，则选择合适的路径长度至关重要。许多分布难以采样，因为它们在某些方向上受到严格约束，而在其他方向上没有约束。探索约束较少的方向最好使用足够长的路径来完成，该路径的长度足以到达该方向上与当前点相距较远的点。但是路径可能会太长，如 图 3 所示。该图左侧所示的路径有点长，因为它会反转方向，然后在可能达到其长度一半左右的路径处结束如果路径长一点，结果可能会更糟，因为路径不仅要花费更长的时间计算，而且可能会在起点附近终止。

对于更复杂的问题，人们不能指望通过观察来选择合适的路径长度如 图 3 所示。找到最小约束变量的线性组合将很困难，而且当典型情况下，最小约束“方向”实际上是非线性曲线或曲面时，将是不可能的。

因此，错误似乎是必要的。对于认为非常困难的问题，$L = 100$ 的路径可能是合适的起点。如果初步运行（具有合适的 $\varepsilon$ ，请参见上文）显示 HMC 仅经过一次迭代即可达到近独立的点，则可以尝试使用较小的 $L$ 值 （除非这些初步运行实际上是足够的，在这种情况下，当然不需要重复运行）。如果相反，在 $L = 100$ 的运行中存在较高的自相关性，那么接下来可以尝试在 $L = 1000$ 的运行中进行。

就像在 第 3.2 节 和 第 3.33 节 末讨论的，我们可能希望随机改变路径的长度 ( 在相当小的间隔内 ) ，以避免选择对于某些变量或变量组合恰好产生近似周期性的路径长度。

4.2.4 采用多种步长

使用 第 4.1 节 中的结果，我们可以利用变量相对尺度的信息来提高 HMC 的性能。这可以通过两种等效的方式来完成。如果 $s_i$ 是适合 $q_i$ 的尺度，我们可以通过设置 $q \prime _i = q_i / s_i$ 来变换 $q$，或者我们可以改用 $K(p)= p^TM{-1}p$ 的动能函数，其中 $M$ 是对角元素 $m_i=1/s_i^2$ 的对角矩阵。

通常利用此信息最方便的方法是对不同位置和动量变量使用不同步长。若要了解其工作原理，请考虑使用 $m_i=1/s_i^2$ 进行蛙跳更新（参照 公式 2.28 至 2.30）：

$$\begin{align*} p_i(t+\varepsilon / 2) &=p_i(t)-(\varepsilon / 2) \frac{\partial U}{\partial q_i}(q(t)) \tag{2.44}\\ q_i(t+\varepsilon) &=q_i(t)+\varepsilon s_i^2 p_i(t+\varepsilon / 2), \tag{2.45}\\ p_i(t+\varepsilon) &=p_i(t+\varepsilon / 2)-(\varepsilon / 2) \frac{\partial U}{\partial q_i}(q(t+\varepsilon)) . \tag{2.46} \end{align*}$$

将 $\left(q^{(0)}, p^{(0)}\right)$ 定义为蛙跳步骤开始时的状态（即 $\left.(q(t), p(t))\right)$，定义 $\left(q^{(1)}, p^{(1)}\right)$ 为最终状态 （即 $(q(t+\varepsilon), p(t+\varepsilon))$，并将 $p^{(1/2)}$ 定义为中途动量（即 $p(t+\varepsilon / 2)$ 。我们现在可以将上面的跳过步骤重写为：

$$\begin{align*} &p_i^{(1 / 2)}=p_i^{(0)}-(\varepsilon / 2) \frac{\partial U}{\partial q_i}\left(q^{(0)}\right), \tag{2.47}\\ &q_i^{(1)}=q_i^{(0)}+\varepsilon s_i^2 p_i^{(1 / 2)} \\ &p_i^{(1)}=p_i^{(1 / 2)}-(\varepsilon / 2) \frac{\partial U}{\partial q_i}\left(q^{(1)}\right) . \tag{2.48} \end{align*}$$

如果现在定义重新定标的动量变量 $\tilde{p}_i=s_i p_i$ ，并逐步调整 $\varepsilon_i=s_i \varepsilon$，则可以将蛙跳更新写为：

$$\begin{align*} \tilde{p}_i^{(1 / 2)} &=\tilde{p}_i^{(0)}-\left(\varepsilon_i / 2\right) \frac{\partial U}{\partial q_i}\left(q^{(0)}\right), \tag{2.49}\\ q_i^{(1)} &=q_i^{(0)}+\varepsilon_i \tilde{p}_i^{(1 / 2)} \tag{2.50}\\ \tilde{p}_i^{(1)} &=\tilde{p}_i^{(1 / 2)}-\left(\varepsilon_i / 2\right) \frac{\partial U}{\partial q_i}\left(q^{(1)}\right) . \tag{2.51} \end{align*}$$

这就像所有 $m_i=1$ 的蛙跳式更新，但对于不同的 $(q_i，p_i)$ 对具有不同的步长。当然，$(q,\tilde{p})$ 的连续值不再可以解释为在一致的时间点遵循哈密顿动力学，但这对于在 HMC 中使用这些路径没有任何意义。请注意，当我们对每个路径之前的动量进行采样时，每个 $p_i$ 均独立于均值为零且方差为 $1$ 的高斯分布而不论 $s_i$ 的值如何。这种多步长方法通常更方便，尤其是在估计比例尺 $s_i$ 不固定的情况下，如 第 4.5 节 所述，动量仅部分刷新( 第 5.3 节 )。

4.3 将 HMC 与其他 MCMC 更新组合在一起

对于某些问题，MCMC 仅使用哈密顿蒙特卡洛法将是不可能的或不可取的。当某些变量是离散变量并且对数概率的对数概率密度的导数相对于某些变量而言昂贵或无法计算时，将有必要进行非 HMC 更新的两种情况。然后可以将 HMC 仅适用于其他变量。另一个例子是，当特殊的 MCMC 更新被设计出来时，它可能以 HMC 无法实现的方式 （例如，通过在其他隔离的峰值之间移动） 但不会完全替代 HMC 的方式进行收敛。如下面 第 4.5 节 所述，贝叶斯分层模型也可能最好结合使用 HMC 和其他方法（例如 Gibbs 采样）来处理。

在这种情况下，可以将所有或部分变量的一个或多个 HMC 更新与一个或多个其他更新交替使用，从而使所有变量的期望联合分布保持不变。可以将 HMC 更新视为既保持该相同的联合分布不变性，又可以将 HMC 更改的变量的条件分布保持为不变，前提是给定 HMC 更新期间固定变量的当前值。这些是等效视图，因为可以将联合密度计算为该条件密度乘以固定变量的边际密度，从单个 HMC 更新的角度来看，这只是一个常数，因此可以将其排除在势能之外。

当同时使用 HMC 和其他更新时，与仅执行 HMC 时使用的路径相比，使用 HMC 最好使用更短的路径。这样可以更频繁地执行其他更新，这大概有助于采样。但是，找到最佳折衷方案可能很困难。 HMC 的一种变体减少了这种折衷的需要，将在下面的 第 5.3 节 中进行介绍。

4.4 维度扩展

在 第 3.3 节 中，说明了 HMC 的主要好处之一：HMC 避免了随机游走方法在探索状态空间时存在的低效率问题。对于大多数实际问题而言，这种优势都在某种程度上存在。对于维数较大（中等到高）的问题，HMC 相对于简单随机游走 Metropolis 方法的另一个好处：对于给定的精度水平，随着维数增加，所需计算时间的增长速度会变慢。 注意这里我将只考虑达到平衡后的采样性能，而不是从某种初始状态（如非典型的分布状态）达到平衡所需的时间，因为这很难分析。

4.4.1 通过复制创建维数增加的分布

为了讨论性能随维度变化的情况，需要作出分布随维度 $d$ 变化的某种假设。

这里我将假设通过添加变量的独立副本方式来增加维数，也就是说， $q = (q_1,...,q_d)$ 的势能函数具有 $U(q) = \sum u_i(q_i)$ 的形式，而 $u_i$ 是独立抽取自某些分布的函数。当然，这不是实际问题的真实情况，但是增加维数的一个合理模型。例如，在统计物理学中，大型系统的偏远区域通常几乎是独立的。请注意，独立性假设本身并不重要，因为如 第 4.1 节 所述，当动力学能量函数（或随机游走提议分布）是旋转对称的，则通过旋转坐标系改变独立性，并不会改变 HMC（以及简单随机游走 Metropolis）的性能。

对于这种变量之间相互独立的分布形式，吉布斯采样可以执行得非常好（假设它是可行的），在每次扫描所有变量后产生一个独立点。当单一变量的更新时间不随 $d$ 增长时，对每个变量单独应用 Metropolis 更新也可以很好的工作。不过，这些方法对旋转无法保持不变性，因此其良好性能无法推广到更复杂的分布，我们希望通过下面的分析获得洞察力。

4.4.2 HMC 和随机游走 Metropolis 的维度扩展性能

在这里，我将非正式地讨论了 HMC 和随机游走在维度扩展的表现，部分遵循了克鲁特兹（Cruetz，1988，第三节）的步伐。

克鲁埃茨指出，当使用任何 Metropolis 风格的算法对密度 $P(x)=(1 / Z)exp(-E( X))$ 进行采样时，都遵循下式：

$$1=\mathrm{E}\left[P\left(x^*\right) / P(x)\right]=\mathrm{E}\left[\exp \left(-\left(E\left(x^*\right)-E(x)\right)\right)\right]=\mathrm{E}[\exp (-\Delta)] \tag{2.52}$$

其中 $x$ 是当前状态，假设按照 $P(x)$ 分布，$x^*$ 是提议状态，$\Delta=E\left(x^*\right)-E(x)$ 。 杰森（Jensen）不等式意味着能量差的期望是非负的：

$$\mathrm{E}[\Delta] \geq 0 . \tag{2.53}$$

该不等式通常是严格的。

当 $U(q)=\Sigma u_i\left(q_i\right)$，并且为每个 $i$ 独立生成提议时，我们可以将这些关系应用于单个变量（或变量对）或整个状态。对于单个变量（或对），我将 $E\left(x^*\right)-E(x)$ 记为 $\Delta_1$，其中 $x=q_i$ 、 $E(x)=u_i\left(q_i\right)$，或者 $x=$$\left(q_i, p_i\right)$ 、 $E(x)=u_i\left(q_i\right)+p_i^2 / 2$ 。对于整个状态，我将 $E\left(x^*\right)-E(x)$ 记为 $\Delta_d$，其中 $x=q$ 、$E(x)=U(q)$，或 $x=(q, p)$ 、$E(x)=U(q)+K(p)$ 。对于随机游走 Metropolis 和 HMC，通过复制变量来增加维度将导致能量差异的增加，因为 $\Delta_d$ 是所有变量的 $\Delta_1$ 之和，每个变量都有正均值。这会导致接受概率的降低：等于 $\min \left(1, \exp \left(-\Delta_d\right)\right)$，除非提议分布的宽度或蛙跳步长能够减小以得到补偿。

更具体地说，对于能够独立改变每个变量提议值的随机游走 Metropolis，提议状态和当前状态之间的势能差将是每个变量独立差的总和。如果为每个提议的变化固定一个标准差 $\varsigma$ ，则该势能差的均值和方差都将随 $d$ 线性增长，这将导致接受概率逐渐降低。为了保持合理性能，$\varsigma$ 必须随着 $d$ 的增加而减少。此外，到达一个近独立的点所需的迭代次数将与 $\mathrm{\varsigma}^{-2}$ 成正比，因为探索是通过随机游走进行的。

类似地，当使用 HMC 从各分量独立的 $q$ 分布中采样时（ 使用动能 $K(p)=\Sigma p_i^2 / 2$ ），不同的 $\left(q_i, p_i \right)$ 在路径模拟过程中不存在相互作用：根据哈密顿动力学，每个 $\left(q_i, p_i\right)$ 仅涉及势能方程和动能方程中的一项。因此，路径的模拟时间长度没有必要随着维度增加而增加。然而，路径端点的接受概率与蛙跳近似最终导致的 $H$ 误差有关，这是一个与每个 $\left(q_i, p_i\right)$ 都相关的误差总和。对于固定步长 $\varepsilon$ 和固定路径长度 $\varepsilon L$，$H$ 误差的均值和方差都会随 $d$ 线性增长。如果不被 $\varepsilon$ 的减少抵消的话，这将导致接受概率随着维度的增加逐渐降低。到达一个独立点所需的蛙跳步数将与 $\varepsilon^{-1}$ 成正比。

为了查看哪种方法扩展性更好，需要确定随着 $d$ 增加，必须以多快速度减少 $\varsigma$ 和 $\varepsilon$，以保持合理的接受概率。随着 $d$ 增加和 $\varsigma$ 或 $\varepsilon$ 变为零，$\Delta_1$ 也将变为零。使用 $\exp \left(-\Delta_1\right)$ 作为 $1-\Delta_1+\Delta_1^2 / 2$ 的二阶近似值，连同方程 $2.53$，我们发现

$$\mathrm{E}\left[\Delta_1\right] \approx \frac{\mathrm{E}\left[\Delta_1^2\right]}{2} . \tag{2.54}$$

由此可见，$\Delta_1$ 的方差是 $\Delta_1$ 均值的两倍（当 $\Delta_1$ 较小时），这意味着 $\Delta_d$ 的方差是 $\Delta_d$ 均值的两倍（即使 $\Delta_d$ 不小）。因此，为了获得好的接受概率，我们必须保持 $\Delta_d$ 的平均值接近 $1$，因为同样大的标准差情况下，不会保存大的均值（这会产生相当频繁的接受概率，因为 $\Delta_d$ 偶尔会出现负值）。

对于对称的提议分布，可以观察如何通过直接平均提议和其逆的 $\Delta_1$ 来扩展 $\varsigma$ 。假设某个变量的提议为 $x^∗ = x + c$ ，并假设 $c = a$ 和 $c = -a$ 的可能性相同。将势能 $U(x^∗)$ 近似为二阶的 $U(x)+ cU^\prime(x)+ c^2U^{\prime\prime}(x)/ 2$ ，我们发现 $\Delta_1 = U(x^∗)-U(x)$ 的平均值为 $2U^{\prime\prime}(x)$ ( $c = a$ 和 $c = -a$ ）。在具有标准差 $a$ 的分布上平均和在 $x$ 的分布上的平均分布上，我们看到 $E \left[\Delta_1\right]$ 与 $\varsigma^2$ 成正比。因此，$E \left[\Delta_d\right]$ 与 $d\varsigma^2$ 成正比，因此我们可以通过使 $\varsigma$ 与 $d^{−1 / 2}$ 成比例来保持合理的接受概率。达到一个几乎独立的点所需的迭代次数将与 $\varsigma^{-2}$ 成正比，与之成正比。所需的计算时间通常与 $2$ 成正比。

如 第 2.3 节 末尾所述，使用蛙跳离散化来模拟固定长度的路径时，$H$ 中的误差与 $\varepsilon^2$ 成正比（对于 $\varepsilon$ 足够小）。单个$(q_i，p_i)$ 的 $H$ 误差与 $\Delta_1$ 相同，因此我们看到 $\Delta^2_1$ 与 $\varepsilon^4$ 成正比。然后， 公式 4.19 表示 $E \left[\Delta_1\right]$ 也与 $\varepsilon^4$ 成比例。 $H$ 所有变量的平均总误差 $E \left[\Delta_d\right]$ 将与 $d\varepsilon^4$ 成比例，因此我们必须使 $\varepsilon$ 与 $d^{−1 / 4}$ 成比例，以保持合理的接受概率。因此，达到几乎独立的点的蛙跳更新次数将增加 $d^{1/4}$，并且计算时间通常将增加 $d^{5/4}$，这比 $d^2$ 增长的随机游走 Metropolis 要好得多。

4.4.3 最佳接受概率

通过扩展上面的分析，我们可以确定当使用 $\varsigma$ 或 $\varepsilon$ 的最优选择时，提议的接受概率是多少。这在调整算法时很有帮助，当然，前提是待采样的分布是高维的，并且能够采用复制变量方式对该分布进行充分建模。

为了找到这个接受概率，首先注意到，由于 Metropolis 方法满足细致平衡，因此 $\Delta_d$ 为负数的被接受提议的概率必须等于 $\Delta_d$ 为正数的被接受提议的概率。由于所有具有负 $\Delta_d$ 的提议都被接受，因此接受概率只是提议具有负 $\Delta_d$ 概率的两倍。对于较大的 $d$，中心极限定理意味着 $\Delta_d$ 的分布是高斯分布，因为它是 $d$ 个独立的 $\Delta_1$ 值的总和 （假设每个 $\Delta_1$ 的方差是有限的）。我们在上面看到 $\Delta_d$ 的方差是其均值的两倍，$\mathrm{E}\left[\Delta_d\right]=\mu$。因此，对于大的 $d$，接受概率可以写成如下（Gupta 等，1990）：

$$P(\text { accept })=2 \Phi\left(\frac{(0-\mu)}{\sqrt{2 \mu}}\right)=2 \Phi(-\sqrt{\mu / 2})=a(\mu) \tag{2.55}$$

其中 $\Phi(z)$ 是均值为 $0$、 方差为 $1$ 的高斯变量的累积分布函数。

对于随机游走 Metropolis，获得独立点的成本将与 $1 /\left(a^2\right)$ 成正比，其中 $a$ 是接受概率。我们在上面看到 $\mu=\mathrm{E}\left[\Delta_d\right]$ 与 $\varsigma^2$ 成正比，因此成本遵循成比例：

$$C_{\mathrm{rw}} \propto \frac{1}{(a(\mu) \mu)} \tag{2.56}$$

数值计算表明，当 $\mu=2.8$ 和 $a(\mu)=0.23$ 时，这种情况最小化。

对于 HMC，获得独立点的成本将与 $1 /(a \varepsilon)$ 成正比，正如我们在上面看到的，$\mu$ 与 $\varepsilon^4$ 成正比。由此得到：

$$C_{\mathrm{HMC}} \propto \frac{1}{\left(a(\mu) \mu^{1 / 4}\right)} \tag{2.57}$$

数值计算表明，最小值是在 $\mu=0.41$ 和 $a(\mu)=0.65$ 时得到的。

Roberts 等人（1997）曾经采用更正式的分析获得了随机游走 Metropolis 的最佳接受概率 $23$ 。我在上面推导出的 HMC 最佳接受概率 $65$ 与先前关于此模型的研究（Neal，1994，图 2）以及实际高维问题（Creutz，1988，图 2 和 3；Sexton 和 Weingarten，1992 年，表 1）的经验结果一致。 Kennedy 和 Pendleton (1991) 曾经获得了将 HMC 应用于多元高斯分布及其逆分布的明确和严格的结果。

4.4.4 探索势能分布

上述 HMC 更好的扩展行为主要取决于动量变量的重采样。通过考虑方法探索势能 $U(q)=\sum u_i(q_i)$ 的程度，我们可以看到这一点。由于 $U(q)$ 是独立项的和，因此其标准差将与 $d^{1 / 2}$ 成比例增加。

继 Caracciolo 等人之后（1994），我们注意到单次 Metropolis 更新的预期势能变化将不超过 $1$ 阶。直觉上，大的向上变化不太可能被接受，并且由于 Metropolis 更新满足细致平衡条件，大的向下变化也很少（处于平衡状态）。因为 $U$ 的变化会
遵循随机游走（再次由于细致平衡），至少需要 $(d^{1/2}/ 1)^2 = d$ 阶 Metropolis 更新来探索 $U$ 的分布。

在 HMC 迭代的第一步中，对动量进行重采样变量通常会将动能改变与 $d^{1 / 2}$ 成比例的量，因为动能也是独立项的和，因此具有标准差 $d^{1 / 2}$ 增长（更确切地说，其标准差为 $d^{1 / 2}/2^{1/2}$ ） 。如果 HMC 的第二步提议有一个远点，则到路径，已经在动能和势能之间平均分配。如果接受此路径的终点，则单个 HMC 迭代产生的势能变化将与 $d^{1 / 2}$ 成比例，与其总体变化范围相当。因此，与随机游走的 Metropolis 相反，我们可能希望即使对于高维分布，也仅需进行少量 HMC 迭代就可以移动到几乎独立的点。

分析方法对 $U$ 分布的探索程度如何也可以洞悉其在分布上的性能变量复制不能很好地建模它们，我们将在下一部分中看到。

4.5 分层模型的 HMC

许多贝叶斯模型是分层定义的。大量低级参数的先验分布由较少的高级 “超参数” 确定，而这些“超参数”又可能具有由更高级别的超参数确定的先验分布。例如，在具有许多预测变量的回归模型中，回归系数可能具有高斯先验分布，均值为零，方差为超参数。这个超参数可以被赋予一个广泛的先验分布，因此它的后验分布主要由数据决定。

可以以一种很明显的方式将 HMC 应用于这些模型（在获取方差超参数的对数之后，因此它们将不受约束）。但是，将 HMC 仅应用于较低级别的参数可能会更好，原因我现在将讨论。 （有关将 HMC 应用于变量子集的一般性讨论，请参见 第 4.3 节。）

我将使用我在贝叶斯神经网络模型（Neal，1996a）方面的工作作为示例。这样的模型通常有几组低级参数，每组都有一个相关的方差超参数。这些超参数的后验分布反映了数据的重要方面，例如哪些预测变量与任务最相关。从后验分布中采样这些超参数值的效率通常可以决定 MCMC 方法的整体效率

我仅将 HMC 用于贝叶斯神经网络模型中的低级参数，超参数在 HMC 更新期间是固定的。这些 HMC 更新与超参数的 Gibbs 采样更新交替，它们（在模型的更简单版本中）在给定低级参数的情况下是独立的，并且具有标准形式的条件分布。通过仅将 HMC 用于低级参数，可以使用基于当前超参数值的启发式方法设置所使用的蛙跳步长。 （我使用 第 4.2 节 末尾描述的多步长方法，相当于对不同的参数使用不同的质量值 $m_i$。）例如，网络“权重”的大小在“隐藏单元”之外的连接上” 确定似然函数对连接到隐藏单元的权重变化的敏感程度；因此，这些传出连接上的权重方差对于设置传入连接上的权重的步长很有用。如果超参数被相同的 HMC 更新更改为更改较低级别的参数，则使用它们来设置步长将无效，因为反向路径将使用不同的步长，因此不会回溯原始路径。如果没有设置步长的好方法，低级参数的 HMC 可能效率会低得多。

Choo (2000) 通过使用 HMC 的修改绕过了这个问题，其中通过仅更新超参数的跳跃式步骤与仅更新低级参数的跳跃式步骤交替来模拟路径。这个过程既保持可逆性又保持体积（尽管不一定是辛性），同时允许使用超参数的当前值设置低级参数的步长（反之亦然）。然而，由于分层模型的第二个问题，性能并没有像希望的那样提高。

在这些贝叶斯神经网络模型和许多其他层次模型中，低级参数和超参数的联合分布高度偏斜，概率密度从一个高后验概率区域到另一个区域变化很大。除非控制低级参数方差的超参数具有非常窄的后验分布，否则超参数和低级参数的联合后验密度从方差低到高时变化很大

例如，假设在其高后验概率区域中，方差超参数变化 $4$ 倍。如果这个超参数控制 $1000$ 个低级参数，它们的典型先验概率密度将变化 $21000 = 1.07 × 10301$，对应于 $\log(21000) = 693$ 的势能范围，标准差为 $693/121/2 = 200$（因为均匀分布的方差是其范围的十二分之一）。在 第 4.4 节 结束时，一次 HMC 迭代仅通过对动量变量的重新采样来改变能量，这充其量会导致势能的变化，标准偏差约为 $d^{1/2}/2^{3/2}$。对于此示例，如果有 $1000$ 个低级参数，则为 $11.2$，因此大约需要 $(200 / 11.2)2 = 319$ 个 HMC 迭代才能达到独立点。

对于此示例，使用 Gibbs 抽样可能会获得类似的性能。然而，对于神经网络模型，没有可行的方法使用 Gibbs 采样来进行低级参数的后验分布，但 HMC 可以应用于给定超参数的低级参数的条件分布。然后可以使用 Gibbs 采样来更新超参数。正如我们所见，尝试使用 HMC 更新超参数也不会提高性能，而通过 Gibbs 采样更新它们更容易。

Choo (2000) 尝试了另一种可能改进这一点的方法——通过让 $θ_i = φ_i exp(κ/2)$ 重新参数化低级参数 $θ_i$，所有参数均具有方差 $exp(κ)$，然后使用采样对 $κ$ 和 $φ_i$ HMC。重新参数化消除了 HMC 无法有效采样的概率密度的极端变化。然而，他发现很难为 $κ$ 设置合适的步长，并且 $H$ 中的误差倾向于随着路径长度而增长，这与 HMC 仅用于低级参数的典型情况不同。使用 “退火” 技术（见 第 5.7 节 ）是另一种可能性。

尽管 HMC 不能消除所有困难，但 HMC 对于贝叶斯神经网络模型非常有用——事实上，没有它，对于大多数应用程序来说可能不可行。将 HMC 至少用于低级参数可以为其他层次模型产生类似的好处（例如 Ishwaran，1999），尤其是当这些低级参数的后验相关性很高时。正如在 HMC 的任何应用中一样，通常需要仔细调整步长和路径长度。

5 HMC 的变种

可以以多种方式修改 图 2 的基本 HMC 算法，以提高其效率，或使其可用于更广泛的分布范围。在本节中，我将首先讨论哈密顿方程的蛙跳式离散化的替代方法，并说明 HMC 如何处理对变量（例如，必须为正数的变量）进行约束的分布。然后，我将讨论 HMC 的一种特殊情况（仅完成一个跳过步骤时），并说明如何扩展它以产生另一种避免随机游走的方法，当 HMC 并非所有变量都更新时，这可能会很有用。 HMC 的大多数应用都可以受益于使用一种变体，其中使用状态的“窗口”来增加接受概率。另一种广泛适用的技术是使用哈密顿近似值来计算路径，同时在确定是否接受路径终点时仍通过使用精确的哈密顿函数获得正确的结果。可以通过使用“捷径”方法来辅助 HMC 的调整，该方法可以避免在步长不合适的情况下计算整个路径。退火方法有可能帮助解决具有多种峰值或高度偏态的分布。这里有许多其他变体，我将无法对其进行回顾，例如使用“影子哈密顿函数”，这是不精确模拟所精确保持的。真正的哈密顿语（Izaggu_irre 和汉普顿，2004 年），并且辛普积分方法的使用比跳越方法更为复杂（例如 Creutz 和 Gocksch，1989 年），其中包括 Girolami 等人 （2009 年）最近的一项提议，采用了一个不可分离的哈密顿函数的辛积分子，其最小动能 公式 2.6 取决于 $q$，允许基于局部信息进行 “适应”。

5.1 拆分方法

跳蛙法不是唯一的可逆和保体积的 Hamilton 方程离散化，因此可用于 HMC。已经设计了许多“辛积分方法”，主要用于 HMC 以外的应用（例如，模拟太阳系数百万年以测试其稳定性）。可以设计出比蛙跳法具有更高精度的方法（例如，参见 McLachlan 和 Atela，1992）。随着维数的增加，对 HMC 使用这种方法将产生比跳蛙方法渐近更好的性能。然而，经验表明，越级方法在实践中很难被击败。
尽管如此，值得更全面地了解如何模拟哈密顿动力学，因为这也指出了如何处理对变量的约束，以及可能的改进，例如利用部分解析解。

5.1.1 拆分哈密顿量

哈密顿动力学的许多辛离散化可以通过将哈密顿量“拆分”成几个项来导出，然后对于每个项连续地模拟由该项定义的动力学一个小的时间步长，然后重复这个过程直到所需的总模拟时间到了。如果可以分析地完成每个项的模拟，我们将获得一个可行的动力学辛近似。

Leimkuhler 和 Reich（2004 年，第 4.2 节）以及 Sexton 和 Weingarten（1992 年）描述了这种通用方案。假设哈密顿量可以写成 $k$ 项的总和，如下所示：

$$H(q, p)=H_1(q, p)+H_2(q, p)+\cdots+H_{k-1}(q, p)+H_k(q, p) 。$$

还假设我们可以精确地实现基于每个 $H_i$ 的哈密顿动力学，对于 $i=1,\ldots,k$，其中 $T_{i,\varepsilon}$ 是通过应用基于 $H_i$ 的动力学定义的映射时间$\varepsilon$。如 Leimkuhler 和 Reich 所示，如果 $H_i$ 是两次可微的，则这些映射的组成，$T_{1, \varepsilon} \circ T_{2, \varepsilon} \circ \cdots \circ T_{k-1 , \varepsilon} \circ T_{k, \varepsilon}$ 是基于 $H$ 的哈密顿动力学的有效离散化，它将在 $\varepsilon$ 变为零时重现极限内的精确动力学，全局误差为订购 $\varepsilon$ 或更少。

此外，这种离散化将保留体积，并且将是辛的，因为每个 $T_{i, \varepsilon}$ 映射都满足这些属性。如果$H_i$ 的序列是对称的，那么离散化也将是可逆的——即$H_i(q, p)=H_{k-i+1}(q, p)$。如第 5.2.3 节末尾所述，任何可逆方法都必须在 $\varepsilon$ 中具有偶数阶的全局误差（Leimkuhler 和 Reich，2004，第 4.3.3 节），这意味着全局误差必须是阶数 $ \varepsilon^2$ 或更好。

我们可以从哈密顿量的对称分裂中推导出越级法。如果$H(q, p)=U(q)+K(p)$，我们可以将哈密顿量写为

$$H(q, p)=\frac{U(q)}{2}+K(p)+\frac{U(q)}{2}，$$

这对应于 $H_1(q, p)=H_3(q, p)=U(q) / 2$ 和 $H_2(q, p)=K(p)$ 的拆分。基于 $H_1$ 的哈密顿动力学是（ 公式 2.1 和 公式 2.2 ）：

$$\begin{align*} \frac{d q_i}{d t} &=\frac{\partial H_1}{\partial p_i}=0, \\ \frac{d p_i}{d t} &=-\frac{\partial H_1}{\partial q_i}=-\frac{1}{2} \frac{\partial U}{\partial q_i} 。 \end{align*}$$

对时间 $\varepsilon$ 应用这个动态，只需将 $-(\varepsilon / 2) \partial U / \partial q_i$ 添加到每个 $p_i$，这是跳跃步骤的第一部分（公式 5.18）。基于$H_2$的动态如下：

$$\begin{align*} &\frac{d q_i}{d t}=\frac{\partial H_2}{\partial p_i}=\frac{\partial K}{\partial p_i}, \\ &\frac{d p_i}{d t}=-\frac{\partial H_2}{\partial q_i}=0 。 \end{align*}$$

如果 $K(p)=\frac{1}{2} \sum p_i^2 / m_i$，对时间 $\epsilon$ 应用这个动态会导致将 $\varepsilon p_i / m_i$ 添加到每个 $q_i$，这是跳跃式步长公式 5.19 的第二部分。最后，$H_3$ 产生了跳跃步的第三部分（公式 5.20），这与第一部分相同，因为 $H_3=H_1$。

5.1.2 拆分以利用部分解析解

拆分可以提供帮助的一种情况是，当势能包含一个可以单独进行分析处理的项时。例如，贝叶斯后验分布的势能将是减去参数的对数先验密度和减去对数似然的总和。如果先验是高斯的，则对数先验密度项将是二次的，并且可以进行解析处理（参见第 5.2.1 节末尾的一维示例）。

我们可以通过使用修改后的拆分来针对这种情况修改越级方法。假设 $U(q)=U_0(q)+U_1(q)$，其中 $U_0$ 与动能 $K(p)$ 一起在分析上是易处理的。我们使用拆分

$$H(q, p)=\frac{U_1(q)}{2}+\left[U_0(q)+K(p)\right]+\frac{U_1(q)}{2}，$$

即 $H_1(q, p)=H_3(q, p)=U_1(q) / 2$ 且 $H_2(q, p)=U_0(q)+K(p)$。 $p$ 的前半步和后半步与普通蛙跳法一致，仅基于 $U_1$。 但中间 $q$ 的完整步骤被替换为基于时刻 $\varepsilon$ 时的哈密顿量 $U_0(q)+K(p)$ 的解析解，该步骤在普通跳蛙法中通常只是将 $\varepsilon p$ 添加到 $q$ 而已。

通过这个过程，应该可以使用更大的步长（因此在路径中使用更少的步长），因为部分势能已经被分离出来并被精确处理。然而，精确处理先验的好处可能是有限的，因为先验通常受可能性支配。

5.1.3 以可变的计算成本拆分势能

如果势能函数可以拆分为两项，拆分也有帮助，其中一项需要比另一项更少的计算时间来评估（Sexton 和 Weingarten，1992)。假设 $U(q)=U_0(q)+U_1(q)$，其中 $U_0$ 的计算成本比 $U_1$ 便宜，并且让动能为 $K(p)$。对于某些 $M>1$ ，我们可以使用以下拆分：

$$H(q, p)=\frac{U_1(q)}{2}+\sum_{m=1}^M\left[\frac{U_0(q)}{2 M}+\frac{K(p )}{M}+\frac{U_0(q)}{2 M}\right]+\frac{U_1(q)}{2}。$$

我们将 $k=3 M+2$ 项标记为 $H_1(q, p)=H_k(q, p)=U_1(q) / 2$ 并且对于 $i=1，\ldots, M$, $ H_{3 i-1}(q, p)=H_{3 i+1}(q, p)=U_0(q) / 2 M$ 和 $H_{3 i}(q, p)=K(p )M$ 。由此产生的离散化可以看作是一种嵌套的跳跃式方法。 $M$ 内部越级步仅涉及 $U_0$，并使用有效步长 $\varepsilon / M$。对于$p$，外部跳蛙步仅使用$U_1$ 进行半步，并用$M$ 内部跳蛙步替换中间$q$ 的更新。
如果 $U_0$ 的计算成本比 $U_1$ 便宜得多，我们可以为 $M$ 使用较大的值，而不会增加太多计算时间。我们可以使用的步长 $\varepsilon$ 将主要受 $U_1$ 的属性限制，因为 $U_0$ 的有效步长要小得多，$\varepsilon / M$。通常可以使用比标准跳跃方法更大的 $\varepsilon$，因此我们将需要更少的路径步数，更少的 $U_1$ 计算。

5.1.4 根据数据子集进行拆分

当从独立数据点的贝叶斯模型的后验分布中采样时，可以通过将势能拆分为数据子集的项来节省计算时间。

假设我们将数据划分为子集 $S_m$，对于 $m=1，\ldots，M$，通常大小大致相等。然后我们可以将对数似然函数写为 $\ell(q)=\sum_{m=1}^M \ell_m(q)$，其中 $\ell_m$ 是基于 $S_m$ 中数据点的对数似然函数。如果 $\pi(q)$ 是参数的先验密度，我们可以令$U_m(q)=-\log (\pi(q)) / M-\ell_m(q)$，并将哈密顿量拆分为如下：

$$H(q, p)=\sum_{m=1}^M\left[\frac{U_m(q)}{2}+K(p) / M+\frac{U_m(q)}{2}\right ]$$

也就是说，我们让 $k=3 M$ 项为 $H_{3 m-2}(q, p)=H_{3 m}(q, p)=U_m(q) / 2$ 和 $H_{ 3 m-1}(q, p)=$$K(p)/m$。由此产生的步长为$\varepsilon$ 的离散化有效地执行了步长为$\varepsilon / M$ 的$M$ 跨越步，其中第$m$ 步使用$M U_m$ 作为势能函数。

如果数据集是冗余的，具有许多相似的数据点，则此方案可能是有益的。然后我们期望 $M U_m(q)$ 与 $U(q)$ 大致相同，并且我们可能希望我们可以将 $\varepsilon$ 设置为比标准跳跃方法大 $M$ 倍，获得$M$ 次计算的类似结果。然而，在实践中，路径末端的 $H$ 误差将大于标准的跳跃式，因此增益将小于此值。我发现（Neal，1996a，第 3.5.1 和 3.5.2 节）该方法对神经网络模型很有用，尤其是与下面第 5.5.4 节中描述的窗口 HMC 程序结合使用时。

请注意，与上面的其他示例不同，这种拆分不是对称的，因此产生的离散化是不可逆的。但是，它仍然可以用于生成 $\mathrm{HMC}$ 的提议，只要子集的标记在每次迭代中都是随机的，因此反向路径与正向路径具有相同的生成概率。 （但是，这种拆分的一些对称变化可能会产生更好的结果。）

5.1.5 处理约束

基于拆分的参数显示了如何处理对被采样变量的约束。在这里，我将只考虑对变量的某些子集的单独约束，对 $q_i$ 的约束采用 $q_i \leq u_i$ 或 $q_i \geq l_i$ 或两者的形式。对于任何可微函数 $G$，类似的方案可以处理形式为 $G(q) \geq 0$ 的约束。

我们可以通过让违反任何约束的 $q$ 值的势能是无限的来对变量施加约束，这将使这些点的概率为零。要了解如何处理这种无限的势能，我们可以查看接近无穷大的势能函数的极限，以及相应的动力学极限。
为了说明，假设$U_*(q)$ 是忽略约束的势能，并且$q_i$ 被约束为小于$u_i$。我们可以将极限作为以下势能函数的 $r \rightarrow \infty$（这是可以使用的众多函数之一）：

$$U(q)=U_*(q)+C_r\left(q_i, u_i\right), \quad \text { 其中 } C_r\left(q_i, u_i\right)= \begin{cases}0, & \text { if } q_i \leq u_i, \\ r^{r+1}\left(q_i-u_i\right)^r, & \text { if } q_i>u_i .\end{cases}$$

很容易看出 $\lim _{r \rightarrow \infty} C_r\left(q_i, u_i\right)$ 对于任何 $q_i \leq u_i$ 都是零，对于任何 $q_i>u_i$ 都是无穷大。对于任何有限$r>1，U(q)$ 是可微的，所以我们可以用它来定义哈密顿动力学。

为了模拟基于这个 $U(q)$ 的动力学，在动能 $K(p)=\frac{1}{2} \sum p_i^2 / m_i$ 的情况下，我们可以使用公式 5.29 的拆分，与 $U_1(q)=U_*(q)$ 和 $U_0(q)=C_r\left(q_i, u_i\right)$ ：

$$H(q, p)=\frac{U_*(q)}{2} + \left[C_r\left(q_i, u_i\right)+K(p)\right] + \frac{U_*(q)} {2}$$

这产生了一种跳跃方法的变体，其中 $p$（方程 $5.18$ 和 5.19）的半步保持不变，但修改了 $q$（方程 5.19）的全步以说明对 $q_i 的约束。在计算 $q_i^{\prime}=q_i(t)+\varepsilon p_i(t+\varepsilon / 2) / m_i$ 后，我们检查 $q_i^{\prime}>u_i$。如果不是，$C_r\left(q_i, u_i\right)$ 的值必须在从 $q_i$ 到 $q_i^{\prime}$ 的整个路径中为零，我们可以设置 $q(t+\varepsilon) $ 到 $q_i^{\prime}$。但是如果$q_i^{\prime}>u_i$，则基于哈密顿量$C_r\left(q_i, u_i\right)+K(p)$ 的动力学将受到$C_r$ 项的影响。这个术语可以看作是一座陡峭的山丘，当 $q_i$ 越过 $u_i$ 时，它将被攀登，直到到达 $C_r$ 等于 $\frac{1}{2} p_i^2 / m_i$，此时 $p_i$ 将为零。 （如果 $r$ 足够大，因为它会在 $r \rightarrow \infty$ 的极限内，这点将在整个步骤结束之前到达。）然后我们将下山，$p_i $ 越来越负值，直到我们再次达到 $q_i=u_i$，此时 $p_i$ 将只是 $p_i$ 原始值的负数。然后我们继续，现在朝着相反的方向移动，远离上限。

如果几个变量都有约束，我们必须对每个变量都遵循这个过程，如果一个变量同时有上限和下限，我们必须重复这个过程，直到没有一个约束被违反。最终结果是方程 $5.19$ 的 $q$ 的完整步骤被图 5.8 所示的过程替换。直观地说，路径只是从约束给出的“墙壁”反弹。如果 $U_*(q)$ 是常数，这些反弹是动力学中唯一有趣的方面，并且该过程有时被称为“台球”（例如，参见 Ruján，1997）。

5.2 朗之万方法

哈密尔顿蒙地卡罗的一种特殊情况是：用于提出新状态的路径仅包含一个跳跃步。

假设我们使用动能 $K(p)=(1/2) \sum p^2_i$。 HMC 的一个跳跃步骤的迭代可以用以下方式表示。我们从均值为 $0$ 、方差为 $1$ 的高斯分布中对动量变量 $p$ 采样，然后提出新值 $q^*$ 和 $p^*$，如下所示：

$$\begin{align*} &q_i^*=q_i-\frac{\varepsilon^2}{2} \frac{\partial U}{\partial q_i}(q)+\varepsilon p_i \\ &p_i^*=p_i-\frac{\varepsilon}{2} \frac{\partial U}{\partial q_i}(q)-\frac{\varepsilon}{2} \frac{\partial U}{\partial q_i}\left(q^*\right) . \end{align*}$$

我们接受 $q^*$ 作为新状态的概率为：

$$\min \left[1, \exp \left(-\left(U\left(q^*\right)-U(q)\right)-\frac{1}{2} \sum_i\left(\left(p_i^*\right)^2-p_i^2\right)\right)\right]$$

否则保持 $q$ 作为新状态。方程 $5.30$ 在物理学中被称为“朗之万方程”的一种，因此这种方法在晶格场理论文献中被称为朗之万蒙特卡罗 (LMC)（例如，Kennedy，1990）。

也可以删除任何明确提及的动量变量，并将此方法视为执行 Metropolis-Hastings 更新，其中 $q^*$ 从高斯分布提出，其中 $q_i^*$ 是独立的，具有 $q_i 的平均值-\left(\varepsilon^2 / 2\right)\left[\partial U / \partial q_i\right](q)$ 和 $\varepsilon^2$ 的方差。由于这个提议不是对称的，它必须根据 $q^*$ 和 $q$ 的概率密度的比率以及从 $q^* 提出 $q$ 的概率密度的比率来接受或拒绝。 $，反之亦然（Hastings，1970）。要使用一个跳跃步查看与 $\mathrm{HMC}$ 的等价性，我们可以将 Metropolis-Hastings 接受概率写成如下：

$$\min \left[1, \frac{\exp \left(-U\left(q^*\right)\right)}{\exp (-U(q))} \prod_{i=1}^d \frac{\exp \left(-\left(q_i-q_i^*+\left(\varepsilon^2 / 2\right)\left[\partial U / \partial q_i\right]\left(q^*\right)\right)^2 / 2 \varepsilon^2\right)}{\exp \left(-\left(q_i^*-q_i+\left(\varepsilon^2 / 2\right)\left[\partial U / \partial q_i\right](q)\right)^2 / 2 \varepsilon^2\right)}\right]$$

要看到与 公式 5.14 相同，请注意，使用 公式 5.12 和 公式 5.13 ，我们可以编写

$$\begin{align*} p &=\frac{1}{\varepsilon}\left[q_i^*-q_i+\frac{\varepsilon^2}{2} \frac{\partial U}{\partial q_i}(q)\right], \\ p^* &=-\frac{1}{\varepsilon}\left[q_i-q_i^*+\frac{\varepsilon^2}{2} \frac{\partial U}{\partial q_i}\left(q^*\right)\right] . \end{align*}$$

将这些代入公式 5.32 后，很容易看出与公式 5.33 的等价性。

在这个 Metropolis-Hastings 形式中，LMC 方法首先由 Rossky 等人提出。 （1978 年）用于物理模拟。也可以使用没有接受/拒绝步骤的近似 Langevin 方法（有关此问题的讨论，请参见 Neal，1993 年，第 5.3 节） - 例如，在 Grenander 和 Miller (1990) 的一篇关于复杂模型的统计推断的论文中，其中 J. Besag (p. 591) 的讨论中还提出了接受/拒绝步骤。
虽然 LMC 可以看作是 HMC 的一个特例，但它的性质却大相径庭。由于 LMC 更新是可逆的，并且通常只对状态进行很小的更改（因为 & 通常不会很大），因此 LMC 将通过低效率的随机游走来探索分布，就像随机游走 Metropolis 更新一样。

然而，随着维度的增加，LMC 比随机游走 Metropolis 具有更好的缩放行为，这可以从与第 5.4.4 节中的分析平行的分析中看出（Creutz，1988；Kennedy，1990）。跳跃步的局部误差为 $\varepsilon^3$，因此 $\mathrm{E}\left[\Delta_1^2\right]$，一个变量在 $H$ 中的平均平方误差为阶 $\varepsilon^6$。从方程 5.27 中，E [ $\left.\Delta\right]$ 也将是 $\varepsilon^6$ 阶，并且具有 $d$ 自变量，$\mathrm{E}\left[\Delta_d\right] $ 的阶数为 $d\varepsilon^6$，因此 $\varepsilon$ 必须按 $d^{-1 / 6}$ 缩放以保持合理的接受率。由于 LMC 通过随机游走探索分布，因此达到几乎独立点所需的迭代次数将与 $\varepsilon^{-2}$ 成正比，$d^{1 / 3}$达到几乎独立点的计算时间随着 $d^{4 / 3}$ 的增长而增长。这比随机游走 Metropolis 计算时间的 $d^2$ 增长要好，但比 HMC 与足够长以达到几乎独立点的路径一起使用时的 $d^{5 / 4}$ 增长要差。

我们还可以找到当使用最优 $\varepsilon$ 时，当对具有重复 $d$ 次的独立变量的分布进行抽样时，LMC 的接受率将是多少。至于随机游走 Metropolis 和 $\mathrm{HMC}$，接受率由 $\mu=\mathrm{E}\left[\Delta_d\right]$ 由方程 $5.28$ 给出。使用 $\mathrm{LMC}$ 获得几乎独立点的成本与 $1 /\left(a(\mu) \varepsilon^2\right)$ 成正比，因为 $\mu$ 与 $\varepsilon 成正比^6$，我们可以将成本写为

$$C_{\mathrm{LMC}} \propto \frac{1}{\left(a(\mu) \mu^{1 / 3}\right)}$$

数值计算表明，当 a(μ)为 0.57 时，这是最小化的，这是 Roberts and Rosenthal(1998)更正式获得的结果。如果对要采样的分布的 LMC 行为类似于对复制依赖变量进行采样时的行为，则这可能很有用。

5.3 动量的部分更新

LMC 算法中使用的单步蛙跳通常不足以移动到近独立的点，因此 LMC 将通过低效的随机游走来探索分布。这就是为什么 HMC 通常与许多跨越式步骤的路径一起使用。正如 Horowitz (1991) 所提出的，即使路径仅包含一个越级步，也可以抑制随机游走行为的另一种方法是仅部分刷新路径之间的动量。

假设动能具有典型形式$K(p)=p^T M^{-1} p / 2$。 $p$ 的以下更新将使动量分布保持不变（均值为 0 且协方差为 $M$ 的高斯分布）：

$$p^{\prime}=\alpha p+\left(1-\alpha^2\right)^{1 / 2} n .$$

这里，$\alpha$ 是区间 $[-1,+1]$ 中的任意常数，$n$ 是均值为零且协方差矩阵 $M$ 的高斯随机向量。要看到这一点，请注意如果 $p$ 具有所需的高斯分布，则 $p^{\prime}$ 的分布也将是高斯分布（因为它是独立高斯的线性组合），均值为 0，协方差为 $\ alpha^2 M+\left(1-\alpha^2\right) M=M$。

如果 $\alpha$ 仅略小于 $1$，则 $p^{\prime}$ 将类似于 $p$，但这种重复更新最终会为动量变量生成一个几乎独立于初始值的值。当 $\alpha=0$ 时，$p^{\prime}$ 只是设置为从其高斯分布中抽取的随机值，与之前的值无关。请注意，当 $M$ 是对角线时，每个动量变量 $p_i$ 的更新独立于其他动量变量的更新。

方程 $5.34$ 的部分动量更新可以替代标准 HMC 算法中动量的完全替换。这给出了一个广义的 HMC 算法，其中一个迭代由三个步骤组成：

1. 使用公式 5.34 更新动量变量。设新动量为$p^{\prime}$。
2. 提出一个新状态 $\left(q^*, p^*\right)$，通过应用 $L$ 步长为 $\varepsilon$ 的跳跃步，从 $\left(q, p^{\prime }\right)$，然后否定动量。以概率接受 $\left(q^*, p^*\right)$

$$\min \left[1, \exp \left(-U\left(q^*\right)+U(q)-K\left(p^*\right)+K\left(p^{\prime}\right)\right)\right] .$$

如果 $\left(q^*, p^*\right)$ 被接受，令 $\left(q^{\prime \prime}, p^{\prime \prime}\right)=\left(q^ *, p^*\right)$;否则，令 $\left(q^{\prime \prime}, p^{\prime \prime}\right)=\left(q, p^{\prime}\right)$。

3. 对动量求反，使得新状态为$\left(q^{\prime \prime},-p^{\prime \prime}\right)$。

每个步骤中的转换- $(q, p) \rightarrow\left(q, p^{\prime}\right),\left(q, p^{\prime}\right) \rightarrow\left( q^{\prime \prime}、p^{\prime \prime}\right)$ 和 $\left(q^{\prime \prime}、p^{\prime \prime}\right) \rightarrow$ $\left(q^{\prime \prime},-p^{\prime \prime}\right)$-保持$(q, p)$ 的规范分布不变。因此，整个更新也使规范分布保持不变。三个过渡也都满足细化平衡，但三者的顺序组合不满足细化平衡（$\alpha=0$ 时除外）。这是至关重要的，因为如果组合是可逆的，当 $L$ 很小时它仍然会导致随机游走行为。

请注意，省略上面的步骤 (3) 将产生一个有效的算法，但是，远离抑制随机游走，该方法（$\alpha$ 接近 1）将产生几乎来回运动，因为运动方向会与步骤（2）中接受的每个路径反向。随着步骤 (3) 中的反转，只要步骤 (2) 中的路径被接受，运动就会沿相同方向继续，因为 $p$ 的两个否定将取消。每当路径被拒绝时，运动就会反转，因此如果要抑制随机游走行为，则拒绝率必须保持较小。

如果 $\alpha = 0$，则上述算法与标准 HMC 相同，因为步骤 (1) 将完全替换动量变量，步骤 (2) 与标准 HMC 相同，步骤 (3) 将不起作用，因为无论如何，在下一次迭代的步骤（1）中，动量将立即被替换。

由于该算法可以被视为标准 HMC 的泛化，带有一个额外的 $α$ 参数，人们可能会认为它会提供改进，前提是调整 $α$ 以获得最佳性能。然而，Kennedy 和 Pendleton (2001) 表明，当该方法应用于高维多元高斯分布时，只能获得很小的常数因子改进，没有更好的维度缩放。使用长路径（ $L$ 大）和不是非常接近 $1$（但不是零，因此最佳选择不是标准 HMC）的 $α$ 值可以获得最佳性能。如果 $L$ 很小，则需要保持非常低的拒绝率（通过使用小的 $\varepsilon$ ）来抑制随机游走，这使得该方法不如标准 HMC 有利。

令人失望的是，在对多元高斯进行采样时，这种泛化只获得了很小的改进，这是由于可能也适用于其他分布的限制。但是，该方法可能比人们想象的更有用。出于 第 4.3 节 和 4.5 节 中讨论的原因，我们经常将 HMC 更新与其他 MCMC 更新结合起来（也许对于 HMC 未更改的变量）。然后可能会在使用长路径以提高 HMC 效率和使用较短路径以便可以更频繁地完成其他 MCMC 更新之间进行权衡。如果出于这个原因要使用比最优路径更短的路径，则将 $α$ 设置为大于零可以减少否则会导致的随机游走行为。

此外，可以使用下面描述的“窗口”方法来降低拒绝率。结合窗口方法对部分动量刷新的分析可能会发现，使用中等长度的路径以及大于零的 $α$ 值会产生更大的改进。

5.4 接受概率窗口

图 3（右图）显示了 $H$ 中的误差如何沿用蛙跳法计算的典型路径变化。由于在最受限的方向（或多个方向，对于高维分布）上模拟运动时的误差，会发生快速振荡，周期为 $2$ 到 $3$ 步。当使用长路径为 HMC 提出状态时，本质上是随机的，该路径是否结束于 $H$ 中的误差为负或接近于零的状态，因此将以接近于 $1$ 的概率被接受，或者它是否恰好以 $H$ 中具有较大正误差的状态结束，并且相应地较低的接受概率。如果我们能以某种方式消除这些振荡，我们可能会获得所有路径的高接受概率。我介绍了一种实现这一结果的方法，该方法在路径的开始和结束处使用状态 “窗口”（Neal，1994）。

在这里，我将把该方法作为一种通用技术的应用来介绍，在该技术中，我们概率地映射到不同空间中的状态，在这个新空间中执行马尔可夫链转换，然后概率地映射回原始状态空间（Neal， 2006）。

我们的原始状态空间由位置和动量变量对 $(q,p)$ 组成。我们将映射到 $W$ 对序列 $[(q0,p0), . . . ,(qW−1,pW−1)]$，其中 $i> 0$ 时的每个 $(q_i,p_i)$ 是对 $(q_{i−1},p_{i−1})$ 应用一个蛙跳步（具有一些固定步长 $\varepsilon$）的结果.请注意，即使新空间中的一个点似乎包含的数字是原始空间中一个点的 $W$ 倍，但新空间的实际维数与旧空间相同，因为 $W$ 对的整个序列已由 $(q_0,p_0)$ 确定。

为了从 $(q，p)$ 概率映射到序列 $[\left[(q_0，p_0), \ldots ,(q_{W-1}，p_{W-1})\right]$，我们从$\{0, \ldots ,W-1\}$，并将新状态下的 $(q_s，p_s)$ 设置为当前状态 $(q，p)$。新状态中的其他 $\left(q_i, p_i\right)$ 是使用来自 $\left(q_s, p_s\right)$ 的蛙跳步骤获得的，对于 $i>s$，或向后蛙跳步骤（即完成步长 $-\varepsilon$ ) 为 $i<s$。很容易看出，利蛙跳步保留体积这一事实，如果我们的原始状态以概率密度 $P(q, p)$ 分布，那么获得序列的概率密度为 $\left[\left(q_0 , p_0\right), \ldots,\left(q_{W-1}, p_{W-1}\right)\right]$ 由这个过程是：

$$P\left(\left[\left(q_0, p_0\right), \ldots,\left(q_{W-1}, p_{W-1}\right)\right]\right)=\frac{1}{W} \sum_{i=0}^{W-1} P\left(q_i, p_i\right),$$

因为我们可以从匹配序列中任何对的$(q, p)$ 对中获得这个序列，并且我们将从这些对中的每一个开始生成序列的概率是$1 / W$（仅当$s$ 的随机选择将这对放在序列中的正确位置）。

在映射到 $W$ 对的序列后，我们现在执行 Metropolis 更新，以保持方程 $5.35$ 定义的序列分布不变，然后再映射回原始状态空间。为了获得 Metropolis 提议，我们执行 $L-W+1$ 步跳蛙（对于某些 $L \geq W-1$，从 $\left(q_{W-1}, p_{W-1}\right)$，产生 $\left(q_W, p_W\right)$ 到 $\left(q_L, p_L\right)$。然后我们提出序列 $\left[\left(q_L,-p_L\right), \ldots,\left(q_{L-W+1},-p_{L-W+1}\right)\right]$。我们通过通常的 Metropolis 标准接受或拒绝这个提议序列，接受概率为：

$$\min \left[1, \frac{\sum_{i=L-W+1}^L P\left(q_i, p_i\right)}{\sum_{i=0}^{W-1} P\left(q_i, p_i\right)}\right],$$

与 $P(q, p) \propto \exp (-H(q, p))$。 （注意这里$H(q, p)=H(q,-p)$，并且从提议的序列开始将对称地导致提议的原始序列。）
这次 Metropolis 更新给我们留下了序列 $\left[\left(q_L, p_L\right), \ldots,\left(q_{L-W+1}, p_{L-W+1}\right)\right]$，称为“接受窗口”或序列 $\left[\left(q_0, p_0\right), \ldots,\left(q_{W-1}, p_{W-1}\right) \right]$，称为“拒绝窗口”。 （请注意，如果 $L+1<2 W$，这些窗口将重叠。）我们将所选窗口中的对标记为 $\left[\left(q_0^{+}, p_0^{+}\right), \ ldots,\left(q_{W-1}^{+}, p_{W-1}^{+}\right)\right]$。我们现在通过从这个序列到单个对的概率映射来生成窗口化 $\mathrm{HMC}$ 更新的最终状态，选择 $\left(q_e^{+}, p_e^{+}\right)$

如果所选窗口中的序列根据公式 5.35 分布，则 $\left(q_e^{+}, p_e^{+}\right)$chosen 对将根据 $P(q, p) \propto 分布、exp (-H(q, p))$，根据需要。要看到这一点，让 $\left(q_{e+n}^{+}, p_{e+n}^{+}\right)$ 成为应用 $n$ 跨越步骤的结果（如果 $n <0$ ) 从 $\left(q_e^{+}, p_e^{+}\right)$ 开始。 $\left(q_e^{+}, p_e^{+}\right)$ 将从序列映射到单个对的概率密度可以写成如下，考虑到所有可以包含 $\left(q_e^{+}, p_e^{+}\right)$ 及其概率：

$$\frac{P\left(q_e^{+}, p_e^{+}\right)}{\sum_{i=0}^{W-1} P\left(q_i^{+}, p_i^{+}\right)} \text {. }$$

如果所选窗口中的序列是根据 公式 5.21 分布的，则选择的对(q + e，p + e)将根据 P(q，p)∝exp(-H(q，p))进行分布，如预期的。要看到这一点，令(q + e + n，p + e + n)是从(q + e，p + e)开始应用 nleapfrog 步骤（如果 nn <0 则为后向步骤）的结果。考虑到所有可能包含(q + e，p + e)的序列及其概率，可以将从对等映射到单对映射而产生的(q + e，p + e)概率密度写成如下：

$$\sum_{k=e-W+1}^e\left[\frac{1}{W} \sum_{i=k}^{k+W-1} P\left(q_i^{+}, p_i^{+}\right)\right] \frac{P\left(q_e^{+}, p_e^{+}\right)}{\sum_{i=k}^{k+W-1} P\left(q_i^{+}, p_i^{+}\right)}=P\left(q_e^{+}, p_e^{+}\right)$$

因此，整个过程将保持正确的分布不变性。

当 W> 1 时，不会出现在 第 3.2 节 末尾讨论的遍历性的潜在问题，因为存在非零概率转移到仅一步步跃迁的状态，其中 q 可能不同从当前状态的值任意地选择它。

似乎 HMC 窗口过程要求将所有 2W 状态保存在接受和拒绝窗口中，因为从接受窗口或拒绝窗口中选择一个状态时，这些状态中的任何一个都可能变为新状态。 。但是实际上，实际上最多需要保存三个状态-起始状态，以便可以在初始反向仿真之后恢复正向仿真，加上拒绝窗口的一种状态和接受窗口的一种状态，其中一种将成为新状态在选择了这些窗口之一之后。由于每个窗口中的状态是按顺序生成的，因此决定刚生成的状态是否应替换当前为该窗口保存的状态。假设到目前为止所看到的状态的概率密度之和为 si = p1 +··+ p_i。如果刚产生的状态具有概率密度 p_i + 1，则用概率 p_i + 1 /(si + p_i + 1)替换从该窗口保存的先前状态。

我证明了(Neal，1994)与标准 HMC 相比，使用窗口可以提高 HMC 在多元高斯分布上的系数是两个或多个，其中在某些方向上的标准差要比在其他方向上的标准差大得多，这是因为 公式 5.22 中的接受概率使用了窗口中状态的概率密度的平均值，从而进行了平滑处理路径的不精确模拟会产生 H 振荡。根据经验，发现窗口方法的优势随维数的增加而增加。对于高维分布，使用最佳阶跃大小时的接受概率约为 0.85，大于 HMC 的理论值 0.65（请参见 第 4.4 节 ）。

这些多变量高斯分布的结果是在窗口大小为 W 的情况下获得的，远小于路径长度，L。对于不太规则的分布，使用更大的窗口可能是有利的。当 W = L / 2 时，验收测试确定新状态是从路径的上半部分（包括当前状态）还是后半部分；然后从一半或另一半中选择新状态，其概率与该一半中状态的概率密度成正比。选择 Wguard 可以防止路径的最后几个状态具有较低的概率密度(highH)。

HMC 的窗口化变体可能会使 HMC 的其他变体更具吸引力。一种这样的变体( 第 5.1 节 )将哈密顿函数分成与数据子集相对应的许多项，这往往会导致更高的错误（同时节省计算量）。在窗口上平均时，错误的影响较小。如 第 5.3 节 所述，使用部分动量刷新时，需要非常低的抑制率。使用窗口更容易获得低拒绝概率（即，不需要大幅度减少压力），这使得部分动量更新更加有吸引力。

Qin 和 Liu(2001)引入了窗口 HMC 的一种变体。在它们的版本中，跳蛙步骤从开始状态开始，接受窗口由这些步骤的最后一个之后的状态组成。然后从接受窗口中选择一个与其概率密度成正比的概率状态。如果选择的状态在结束之前是 kkstates，则从开始状态开始执行 k 向后蛙跳跳跃，并且由这些阶梯所找到的状态以及从开始状态开始直至 w-k-1 的那些状态将形成拒绝窗口。接受窗然后由 公式 5.22 的类似物赋予概率成为下一个状态；秦和刘的过程与原始的窗口 HMC 过程非常相似。

Qin 和 Liu 的过程的一个缺点是，当拒绝接受窗口时，状态保持不变，而在原始过程中，从拒绝窗口中选择状态（可能是当前状态，但通常不会是）。唯一的不同是，根据秦和刘的过程，从当前状态到接受状态的步数范围从 L-W + 1 到 L（平均 L-(W + 1)/ 2)，而从 L-2W + 2 到 L（平均 L-W + 1）对于原始的窗口 HMC 程序，在 Qin 和 Liu 的过程中，所计算的蛙跳步数从 LtoL + W-1 变化，并且在原始过程中固定为 L。如果 W≪L，则这些差异是微小的。秦和林声称他们的程序在高维多元高斯分布上比原始方法有更好的表现，但是他们的实验是有缺陷的。

秦和刘(2001)也引入了更有用的想法，即在接受窗口和拒绝窗口中对状态进行加权加权，这是非常有用的。也可以合并到原始过程中。当从当前状态映射到一系列加权状态时，选择当前状态的位置的概率等于权重，当计算接受概率或从接受或拒绝窗口中选择状态时，状态的概率密度乘以状态的概率重量。秦和刘使用权重来支持距当前状态更远的状态，这通常可以通过导致移动到某个远点而有用，同时如果远点的概率密度较低，则可以选择一个更近的点。或者，如果将窗口视为消除 H 误差的一种方法，则采用权重更好的“低通滤波器”的对称权重将是有意义的。

5.5 近似路径方法

哈密顿蒙特卡罗的有效性并不取决于在模拟路径时是否使用正确的哈密顿函数。事实上，我们可以改用一些近似的哈密顿函数，只要在其动力学路径仿真过程中，能够可逆且保持体积，不过在计算最终的接受概率时必须使用精确的哈密顿函数。当动能方程为简单形式（ 如 公式 2.23 ） 时，无需寻找动能的近似值，但势能通常更为复杂且计算成本较高 — 例如，势能可能涉及许多数据点的对数似然之和，而该数据又无法被简单的充分统计量汇总。使用许多蛙跳步骤的路径时，如果可以快速、准确地近似势能，就可以节省很多计算时间，同时仍然获得准确结果（区别于继承自 MCMC 的固有采样变化）。

存在许多近似势能的有用方法。例如，如果势能的计算需要迭代的数值方法，则迭代次数可能少于获得机器精度结果所需的次数。在贝叶斯统计应用中，可以通过查看数据子集来以低成本获得对非归一化后验密度（其对数给出势能）的近似值。通常这可能不是一个好的策略，但我发现它对高斯过程模型很有用（Neal，1998；Rasmussen 和 Williams，2006），因为其计算时间是数据点数量的立方，所以即使点数的小幅减少也会产生有用的加速。

Rasmussen(2003) 提出了将势能建模为高斯过程的近似方法，该过程根据初始探索阶段所选位置的势能值来推断。该方法仅假设势能函数的平滑度，因此可以广泛使用。但它受到高斯过程推断成本的限制，因此对于中等维数的问题最有用，因为其精确势能估计代价非常高昂。

据我所知，一个有趣的可能性是将 “精确势能” 表达为 “近似势能” 和 “近似误差” 的和，然后应用 第 5.1 节 中描述的分裂技术 — 利用近似的解析可处理性（例如，对于高斯近似，具有二次势能），或低的计算成本，因此可以使用许多小步骤以很少的成本准确地模拟其动力学。如果被近似项消除的势能变化量对于误差项而言允许一个较大的步长，则这将减少对精确势能的梯度进行评估的次数。

5.6 裁剪路径方法

哈密尔顿蒙地卡罗的一个主要缺点是，如 第 4.2 节 所述，其性能在很大程度上取决于其调整参数的设置，该参数至少包括蛙跳步长 $\epsilon$ 和蛙跳步数 $L$ 以及变化的窗口等。HMC 也有其他可调整的参数。路径长度($\epsilon L$)的最佳选择取决于目标分布的整体范围，因此要找到合适的路径长度，可能需要检查大量的 HMC 更新。相比之下，少量的蛙跳步数有助于揭示步长大小的好坏，这使得人们努力尝试在 HMC 运行过程中 “自适应” 设置步长的可能性。

Andrieu 和 Thoms(2008) 回顾了自适应 MCMC 方法的最新工作。正如他们所解释的那样，天真地根据以前的更新结果为每个 HMC 更新选择了一个步长（例如，如果前 10 条路径全部被拒绝，则将步长减小 20％；如果前 10 条路径中的少于两条被拒绝，则将步长提高 20％的扩展和变型）破坏了正确性的证明（特别是该过程不再是马尔可夫链），并且通常会从错误分布中生成点。但如果适应程度能够随时间下降，则可以获得正确的结果。此类自适应方法应当也可用于 HMC，这与其他任何可调 MCMC 方法几乎相同。

另一种方法 (Neal，2005，2007) 按时间设计不同节点的调整参数值，然后按照计划执行 MCMC 更新或随机选择参数而不考虑已实现状态，以便可以使用 MCMC 收敛和错误分析的常规证明。但该方法为了在参数不适合于分布时减少计算时间，需要使用经过微调的 MCMC 更新来实现。因此，大部分计算时间被用于为使用合适的参数值而做的更新。实际上，参数调整是从计算角度自适应地设置的，而不是从数学角度自适应地设置。

例如，步长过大的模拟路径可以在仅跳过几步之后就被拒绝，每当单个步骤（或一系列短步骤）导致哈密顿量的变化大于某个阈值时，都会拒绝 — 即如果 $| H(q(t +\epsilon)，p(t +\epsilon))-H(q( t)，p(t))|$ 大于阈值即拒绝。如果这在路径早期发生，那么在这种不合适的步长上只会浪费很少的计算时间。这样的路径提前终止是有效的，因为任何满足细致平衡的 MCMC 更新，如果被修改以消除某些状态之间的任何一种方式的转换，仍将满足细致平衡。

通过这个简单的修改，我们可以从一些分布中随机选择步长，而不会在那些变得太大的步长上浪费太多时间。但如果将阈值设置得足够小，以便在步长稍微太大时拒绝时，也许会在已经完成大量计算之后，终止本应被接受的路径。当步长小于最佳值时，尝试尽早终止路径会带来类似风险。

当步长似乎不合适时，终止路径的一个不太激进的替代方案是反转路径。假设我们以 $k$ 步为一组执行蛙跳。根据这 $k$ 个步骤上 $H$ 的变化，我们可以测试步长是否不合适 —— 例如，如果 $H$ 在 $k +1$ 个状态上的标准差大于某个上限阈值或小于某个下限阈值时，则无法通过测试（任何对反向序列产生相同决策的准则都是有效的）。当一组 $k$ 步蛙跳测试失败时，路径停留在这组开始的状态，而不是向前移动 $k$ 步，并且动量变量被求反。路径现在将精确地回溯到先前计算的状态（因此不需要重新计算），直到达到初始状态为止，此时将再次计算新状态。如果另一组 $k$ 步测试失败，路径将再次反转，之后路径的整个剩余部分将遍历已计算的状态，允许立即找到其端点而无需进一步计算。

如果没有任何一组 $k$ 步蛙跳测试失败，则此方案与标准 HMC 结果相同。如果在路径早期有两次失败，那么很少有计算时间会浪费在这个（很可能）不适当的步长上。在这两个极端之间，可能会发生一两次逆转，但不会在路径早期发生；路径端点通常不会接近初始状态，因此不会浪费已经执行的、不可忽略的计算（就像路径已经终止时一样）。

这种裁剪方案可以有效地为少量可调参数找到较好的值，当随机抽取这些参数时，经常会合理地选择较好的值。设置大量可调参数是不可行的，例如当维数很高时 公式 2.5 的 “质量矩阵” 中的条目。因为即使在早期发生了两次反转，当很少选中恰当值的情况下，不恰当参数值造成的成本将占主导地位。

5.7 退火方法

到目前为止，标准 HMC 及其变体均与其他局部 MCMC 方法（例如随机游走 Metropolis 法、吉布斯采样法）一样，很难在被低概率区域分割的峰值之间迁移。不过已经设计出几种通用方案来解决孤立峰值问题，这些方案大多考虑从一系列比目标分布更弥散的分布中进行采样，例如、平行退火方法(Geyer，1991; Earl and Deem，2005)、模拟退火方法 (Marinari 和 Parisi，1992)、退火迁移方法(Neal，1996b)、退火重要性采样方法(Neal，2001)等。这些方法大多通过调整 “温度” 参数 $T$ 获得：当 $T = 1$ 时，对应于目标分布； $T$ 值越大，分布的弥散程度越大。这些 “退火” 方法都可以与 HMC 结合使用。类似退火的行为也可以直接整合进 HMC 的新状态提议生成路径之中，被称为 “退火路径”。

在 “退火路径” 方法的最简单版本（Neal，1996b，第 6 节）中，前一半路径的每步蛙跳中，动量变量均乘以某个略大于 $1$ 的因子 $\alpha$；而后一半路径的每步蛙跳中，动量均除以相同的因子 $\alpha$ 。只要结果可逆并且除法/乘法完全配对，就可以用各种方式使用这些乘法和除法。最对称的一种方案是：在前一半路径中，分别在动量更新的前半步之前( 公式 2.28 ) 和后半步之后( 公式 2.30 ) 乘以因子 $\sqrt{\alpha}$。在后一半路径中，动量更新的前半步之前和后半步之后，相应地将动量除以 $\sqrt{\alpha}$。如果路径的跳蛙步数为奇数，则在路径正中间的蛙跳步骤中，将动量更新的上半步之前乘以 $\sqrt{\alpha}$ ，动量更新的后半步之后除以 $\sqrt{\alpha}$ 以保持对称。请注意，大多数 $\sqrt{\alpha}$ 的乘法（或除法）之前（或之后）会存在另外一个乘法（或除法），因此通常可以将其组合成一个完整的 $\alpha$ 因子乘法或除法。

很容易看出，对于这种退火路径，雅可比矩阵的行列式依然是 $1$，与标准 HMC 一样，因此其末端仍然可以用作新状态提议，在接受概率公式中也仍然不用包含雅可比因子。

将动量乘以稍大于 $1$ 的 $\alpha$ 会稍微增加 $H(q，p)$ 的值。如果 $H$ 的初始值是正则分布在 $T = 1$ 处的典型值，则在此乘法之后，$H$ 将会是比 $T$ 值稍高时的典型值。在起步阶段，$H(q，p)= K(p)+ U(q)$ 的变化完全是因为动能 $K(p)$ 发生了变化（因为 $p$ 变大了），但随后的动力学步骤中，$H$ 的增长会在 $K$ 和 $U$ 之间进行分配，并为 $q$ 产生一个比 $T = 1$ 时更为弥散的分布。在经过多次 $p$ 乘以 $\alpha$ 之后，$q$ 的分布就与 $T = 1$ 时完全不同了，这种更为弥散的分布，使得状态更容易在被低概率区域分割的峰值之间进行迁移。后一半路径中除以 $\alpha$ 导致 $H$ 又返回到 $T = 1$ 时的典型值，不过现在可能已经处于不同的峰值了。

如果 $\alpha$ 太大，则退火路径末端的接受概率会变小，因为末端处的 $H$ 值可能会比初始值大得多。要观察到这一点，请考虑仅包含一个蛙跳步骤的路径。如果 $\epsilon = 0$，则该步骤不执行任何操作，在动量的前半步更新之前乘以 $\sqrt{\alpha}$, 在动量的后半步更新之后除以 $\sqrt{\alpha}$ 来精确抵消，因此总体将保持不变，并且路径将会被接受。但是，由于我们希望发生某些事情，因此使用非零的 $\epsilon$，这将导致蛙跳步骤完成时动能下降（平均情况下），因为乘以 $\sqrt{\alpha}$ 后增加的 $H$ 值不再只分配给动能 $K$，而是同时分配给动能 $K$ 和势能 $U$。现在，$p$ 除以 $\sqrt{\alpha}$ 将无法再抵消与 $\sqrt{\alpha}$ 的乘法；与之相反，它会使 $H$ 的减少量小于之前的增加量（平均而言）。此时 $H$ 在路径末端处比初始值增大的趋势，可以被增加的蛙跳步骤缓解（ 如通过某个因子 $R$，将 $\alpha$ 降至 $\alpha^{1/R}$ ），这可以（大致）在达到路径中点时保持足够有效的温度。

图 9 展示了从均值为 $[0,0]$ 和 $[10,10]$ 、方差为 $I$ 和 $2I$ 的双变量双峰值高斯分布中采样的两条退火路径。每条路径由 $200$ 个蛙跳步数组成，步长为 $\epsilon= 0.3$ ，退火温度比例采用 $\alpha= 1.04$。左图显示路径的水平变化；右图显示路径的位置坐标。上面的图是从 $q = \left[-0.4 \quad -0.9\right]$ 和 $p = \left[0.7 \quad -0.9\right]$ 开始的路径，路径终点最终出现在另一种峰值$\left[10 10\right] 附近。底部图是以 $q = \left[0.1 \quad 1.0\right]$ 和 $p = \left[0.5 \quad 0.8\right]$ 开始的路径，路径终点出现在与起步时相同的峰值附近。顶部路径的 $H$ 变化量为 $0.69$，因此接受概率为 $\exp(-0.69) = 0.50$；底部路径的 $H$ 变化量为 $-0.15$，因此接受概率为 $1$。

通过使用这样的退火路径，HMC 可以对这两种完全分离的峰值进行采样：$11\%$ 的路径能够转移到另一种峰值并被接受，而相对来说，标准 HMC 的表现非常差，会在其中一种峰值中被困很长时间。图 9 中所选择的退火路径参数只是为了便于解释图像，并非最佳选择。使用更少的蛙跳步骤 $L$、更大的步长 $\varepsilon$ 和更大的 $\alpha$ 可以获得更为有效的采样，例如 $L = 20$ 、$\epsilon= 0.6$、 $\alpha= 1.5$ 时，能够实现约 $6\%$ 的峰值迁移概率。

上述退火方法的一个基本限制是（对于标准 HMC）：如果 $H$ 值远高于初始状态值，则不可能接受退火路径的终点。这对应于终点处的概率密度远低于当前状态的概率密度。因此，如果一种峰值既高又窄，另一种峰值又低又宽，尤其是在维数较高时，退火方法将无法以相等的总概率在两种峰值之间很好地移动。本人提出了一种针对此问题的改进 (Neal，1999)，其中，计划迁移到的点可以来自于退火路径的任何位置，而不仅仅是终点。不过，需要根据该点的 $H$ 值和累积雅可比因子来进行选择。这种选择计算比较简单，因为 $p$ 和 $\alpha$ 之间乘法的雅可比矩阵行列式为 $\alpha^d$ ，其中 $d$ 是维数。这种改进退火过程不仅可以在不同宽度的峰值之间移动，还可以用于在重尾分布的中心区域和尾部区域之间来回移动。

关于 HMC 这些变种的详细信息，可以在 本人网页 的 R 语言实现部分找到。

参考文献

Alder, B. J. and Wainwright, T. E. (1959). Studies in molecular dynamics. I. General method.Journal of Chemical Physics, 31:459–466.

Andrieu, C. and Thoms, J. (2008). A tutorial on adaptive MCMC.Statistics and Computing,18:343–373.

Arnold, V. I. (1989).Mathematical Methods of Classical Mechanics, 2nd edition. Springer.Bennett, C. H. (1975). Mass tensor molecular dynamics.Journal of Computational Physics,19:267–279.

Caracciolo, S., Pelissetto, A., and Sokal, A. D. (1994). A general limitation on Monte Carloalgorithms of Metropolis type.Physical Review Letters, 72:179–182. Also available athttp://arxiv.org/abs/hep-lat/9307021.

Choo, K. (2000).Learning Hyperparameters for Neural Network Models Using HamiltonianDynamics. MSc. Thesis, Dept. of Computer Science, University of Toronto,available athttp://www.cs.toronto.edu/∼radford/ftp/kiam-thesis.pdf.

Creutz, M. (1988). Global Monte Carlo algorithms for many-fermion systems.PhysicalReview D, 38:1228–1238.

Creutz, M. and Gocksch, A. (1989). Higher order hybrid Monte Carlo algorithms.PhysicalReview Letters, 63:9–12.

Duane, S., Kennedy, A. D., Pendleton, B. J., and Roweth, D. (1987). Hybrid Monte Carlo.Physics Letters B, 195:216–222.

Earl, D. J. and Deem, M. W. (2005). Parallel tempering: Theory, applications, and newperspectives.Physical Chemistry Chemical Physics, 7:3910–3916.

Frenkel, D. and Smit, B. (1996).Understanding Molecular Simulation: From Algorithms toApplications. Academic Press.

Geyer, C. J. (1991). Markov chain Monte Carlo maximum likelihood. InKeramidas, E. M.,editor,Computing Science and Statistics: Proceedings of the 23rd Symposium on the In-terface, pages 156–163. Interface Foundation.

Girolami, M., Calderhead, B., and Chin, S. A. (2009). Riemannian manifold HamiltonianMonte Carlo.http://arxiv.org/abs/arxiv:0907.1100, 35 pages.

Grenander, U. and Miller, M. I. (1990). Representations of knowledge in complex systems.Physics Letters B, 242:437–443.

Gupta, S., Irb ̈ack, A., Karsch, F., and Petersson, B. (1990). The acceptance probability inthe hybrid Monte Carlo method.Physics Letters B, 242:437–443.

Hastings, W. K. (1970). Monte Carlo sampling methods using Markovchains and theirapplications.Biometrika, 57:97–109.

Horowitz, A. M. (1991). A generalized guided Monte Carlo algorithm.Physics Letters B,268:247–252.

Ishwaran, H. (1999). Applications of hybrid Monte Carlo to generalized linear models:quasicomplete separation and neural networks.Journal of Computational and GraphicalStatistics, 8:779–799.

Izagguirre, J. A. and Hampton, S. S. (2004). Shadow hybrid Monte Carlo: an efficientpropagator in phase space of macromolecules.Journal of Computational Physics, 200:581–604.

Kennedy, A. D. (1990). The theory of hybrid stochastic algorithms. InProbabilistic Methodsin Quantum Field Theory and Quantum Gravity, pages 209–223. Plenum Press.

Kennedy, A. D. and Pendleton, B. (1991). Acceptances and autocorrelations in hybrid MonteCarlo.Nuclear Physics B (Proc. Suppl.), 20:118–121.

Kennedy, A. D. and Pendleton, B. (2001). Cost of the generalizedhybrid Monte Carloalgorithm for free field theory.Nuclear Physics B, 607:456–510. Also available athttp://arxiv.org/abs/hep-lat/0008020.

Leimkuhler, B. and Reich, S. (2004).Simulating Hamiltonian Dynamics. Cambridge Uni-versity Press.

Liu, J. S. (2001).Monte Carlo Strategies in Scientific Computing. Springer.

Mackenzie, P. B. (1989). An improved hybrid Monte Carlo method.Physics Letters B,226:369–371.

Marinari, E. and Parisi, G. (1992). Simulated tempering: A new MonteCarlo scheme.Europhysics Letters, 19:451–458.

McLachlan, R. I. and Atela, P. (1992). The accuracy of symplecticintegrators.Nonlinearity,5:541–562.

Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., andTeller, E. (1953).Equation of state calculations by fast computing machines.Journal of Chemical Physics,21:1087–1092.

Neal, R. M. (1993).Probabilistic Inference Using Markov Chain Monte Carlo Methods.Technical Report CRG-TR-93-1. Dept. of Computer Science, University of Toronto.

Neal, R. M. (1994). An improved acceptance procedure for the hybrid Monte Carlo algorithm.Journal of Computational Physics, 111:194–203.

Neal, R. M. (1996a).Bayesian Learning for Neural Networks. Lecture Notes in Statistics,No. 118. Springer-Verlag.

Neal, R. M. (1996b). Sampling from multimodal distributions using tempered transitions.Statistics and Computing, 6:353–366.

Neal, R. M. (1998). Regression and classification using Gaussian process priors (with dis-cussion). In J. M. Bernardo, et al., editor,Bayesian Statistics 6, pages 475–501. OxfordUniversity Press.

Neal, R. M. (1999). Markov chain sampling using Hamiltonian dynamics.Talkat the Joint Statistical Meetings, Baltimore, August 1999, slides are available athttp://www.cs.utoronto.ca/∼radford/ftp/jsm99.pdf.

Neal, R. M. (2001). Annealed importance sampling.Statistics and Computing, 11:125–139.

Neal, R. M. (2005). The short-cut Metropolis method. Technical Report No. 0506,Department of Statistics, University of Toronto (August 2005),28 pages, available athttp://arxiv.org/abs/math.ST/0508060.

Neal, R. M. (2006). Constructing efficient MCMC methods using temporary mapping andcaching. Talk at Columbia University, Dept. of Statistics, 11 December 2006, slides areavailable athttp://www.cs.utoronto.ca/∼radford/ftp/cache-map.pdf.

Neal, R. M. (2007). Short-cut MCMC: An alternative to adaptation. Talk at theThird Workshop on Monte Carlo Methods, Harvard, May 2007, slidesare available athttp://www.cs.utoronto.ca/∼radford/ftp/short-mcmc-talk.pdf.

Qin, Z. S. and Liu, J. S. (2001). Multipoint Metropolis method with application to hybridMonte Carlo.Journal of Computational Physics, 172:827–840.

Rasmussen, C. E. (2003). Gaussian processes to speed up hybridMonte Carlo for expensiveBayesian integrals. In J. M. Bernardo, et al., editor,Bayesian Statistics 7, pages 651–659.Oxford University Press.

Rasmussen, C. E. and Williams, C. K. I. (2006).Gaussian Processes for Machine Learning.MIT Press.

Roberts, G. O., Gelman, A., and Gilks, W. R. (1997). Weak convergence and optimal scalingof random walk Metropolis algorithms.Annals of Applied Probability, 7:110–120.

Roberts, G. O. and Rosenthal, J. S. (1998). Optimal scaling of discrete approximations toLangevin diffusions.Journal of the Royal Statistical Society, Series B, 60:255–268.

Roberts, G. O. and Tweedie, R. L. (1996). Exponential convergence of Langevin distributionsand their discrete approximations.Bernoulli, 2:341–363.

Rossky, P. J., Doll, J. D., and Friedman, H. L. (1978). Brownian dynamics as smart MonteCarlo simulation.Journal of Chemical Physics, 69:4628–4633.

Ruj ́an, P. (1997). Playing billiards in version space.Neural Computation, 9:99–122.

Schmidt, M. N. (2009). Function factorization using warped Gaussian processes. InProceed-ings, Twenty-Sixth International Conference on Machine Learning.

Sexton, J. C. and Weingarten, D. H. (1992). Hamiltonian evolution for the hybrid MonteCarlo method.Nuclear Physics B, 380:665–677