随机采样的动态方法起源于模拟哈密顿动力学下进行变化的物理系统的行为。在马尔科夫链 蒙特卡罗模拟中，目标是从一个给定的概率分布$p(z)$中采样。通过将概率仿真转化为哈密顿系统的形式，我们可以利用哈密顿动力学（Hamiltonian dynamics）的框架。为了与这个领域的文献保持一致，我们在必要的时候会使用相关动态系统的术语，这些术语会随着我们内容的推进而给出定义。 我们考虑的动力学对应于在连续时刻（记作$\tau$）下的状态变量$Z = \{z_i\}$的演化。经典的动力学由牛顿第二定律描述，即物体的加速度正比于施加的力，对应于关于时间的二阶微分方程。我们可以将一个二阶微分方程分解为两个相互偶合的一阶方程，方法是引入中间的动量（momentum）变量$r$，对应于状态变量$z$的变化率，元素为

$$r_i = \frac{dz_i}{d\tau} \tag{11.53}$$

从动力学的角度，$z_i$可以被看做位置（position）变量。因此对于每个位置变量，都存在一个对应的动量变量，位置和动量组成的联合空间被称为相空间（phase space）。

不失一般性，我们可以将概率分布$p(z)$写成下面的形式

$$p(z) = \frac{1}{Z_p}exp(-E(z)) \tag{11.54}$$

其中$E(z)$可以看做状态$z$处的势能（potential energy）。系统的加速度是动量的变化率，通过施加力（force）的方式确定，它本身是势能的负梯度，即

$$\frac{dr_i}{d\tau} = -\frac{\partial E(z)}{\partial z_i} \tag{11.55}$$

使用哈密顿框架重新写出这个动态系统的公式是比较方便的。为了完成这一点，我们首先将动能（kinetic energy）定义为

$$K(r) = \frac{1}{2}\Vert r \Vert^2 = \frac{1}{2}\sum\limits_i r_i^2 \tag{11.56}$$

系统的总能量是势能和动能之和，即

$$H(z,r) = E(z) + K(r) \tag{11.57}$$

其中$H$是哈密顿函数（Hamiltonian function）。使用式（11.53）、（11.55）、（11.56）和（11.57），我们现在可以将系统的动力学用哈密顿方程的形式表示出来，形式为

$$\begin{eqnarray} \frac{dz_i}{d\tau} = \frac{\partial H}{\partial r_i} \tag{11.58} \\ \frac{dr_i}{d\tau} = \frac{\partial H}{\partial z_i} \tag{11.59} \end{eqnarray}$$

在动态系统的变化过程中，哈密顿函数$H$的值是一个常数，这一点通过求微分的方式很容易看出来。

$$\begin{eqnarray} \frac{dH}{d\tau} &=& \sum\limits_i\left\{\frac{\partial H}{\partial z_i}\frac{dz_i}{d\tau} + \frac{\partial H}{\partial r_i}\frac{dr_i}{d\tau}\right\} \\ &=& \sum\limits_i\left\{\frac{\partial H}{\partial z_i}\frac{\partial H}{\partial r_i} - \frac{\partial H}{\partial r_i}\frac{\partial H}{\partial z_i}\right\} = 0 \tag{11.60} \end{eqnarray}$$

哈密顿动态系统的第二个重要性质是动态系统在相空间中体积不变，这被称为Liouville定理（Liouville's Theorem）。也就是说，如果我们考虑变量$(z, r)$空间中的一个区域，那么当这个区域在哈密顿动态方程下的变化时，它的形状可能会改变，但是它的体积不会改变。可以这样证明：我们注意到流场（位置在相空间的变化率）为

$$V = \left(\frac{dz}{d\tau}, \frac{dr}{d\tau}\right) \tag{11.61}$$

这个场的散度为0，即

$$\begin{eqnarray} div V &=& \sum\limits_i\left\{\frac{\partial}{\partial z_i}\frac{dz_i}{d\tau} + \frac{\partial}{\partial r_i}\frac{dr_i}{d\tau}\right\} \\ &=& \sum\limits_i\left\{+\frac{\partial}{\partial z_i}\frac{\partial H}{\partial r_i} - \frac{\partial}{\partial r_i}\frac{\partial H}{\partial z_i}\right\} = 0 \tag{11.62} \end{eqnarray}$$

现在考虑相空间上的联合概率分布，它的总能量是哈密顿函数，即概率分布的形式为

$$p(z,r) = \frac{1}{Z_H}exp(-H(z,r)) \tag{11.63}$$

使用体系的不变性和$H$的守恒性，可以看到哈密顿动态系统会使得$p(z, r)$保持不变。可以这样证明：考虑相空间的一个小区域，区域中$H$近似为常数。如果我们跟踪一段有限时间内的哈密顿方程的变化，那么这个区域的体积不会发生改变，从而这个区域的$H$的值不会发生改变，因此概率密度（只是$H$的函数）也不会改变。

虽然$H$是不变的，但是$z$和$r$会发生变换，因此通过在一个有限的时间间隔上对哈密顿动态系统积分，我们就可以让$z$以一种系统化的方式发生较大的变化，避免了随机游走的行为。

然而，哈密顿动态系统的变化对$p(z, r)$的采样不具有各态历经性，因为$H$的值是一个常数。为了得到一个具有各态历经性的采样方法，我们可以在相空间中引入额外的移动，这些移动会改变$H$的值，同时也保持了概率分布$p(z, r)$的不变性。达到这个目标的最简单的方式是将$r$的值替换为一个从以$z$为条件的概率分布中抽取的样本。这可以被看成Gibbs采样的步骤，因此根据11.3节，我们看到这也使得所求的概率分布保持了不变性。注意，$z$和$r$在概率分布$p(z, r)$中是独立的，我们看到条件概率分布$p(r|z)$是高斯分布，从中我们可以很容易地进行采样。

在这种方法的一个实际应用中，我们必须解决计算哈密顿方程的数值积分的问题。这会引入一些数值的误差，因此我们要设计一种方法来最小化这些误差产生的影响。事实上，可以证明，能够在Liouville定理仍然精确成立的条件下，对积分方法进行修改。这个性质在11.5.2节讨论混合蒙特卡罗算法时很重要。完成这件事的一种方法是蛙跳（leapfrog）离散化。这种方法使用下面的公式对位置变量和动量变量的离散时间近似$z$和$r$进行交替地更新。

$$\begin{eqnarray} \hat{r}_i\left(\tau + \frac{\epsilon}{2}\right) &=& \hat{r}_i(\tau) - \frac{\epsilon}{2}\frac{\partial E}{\partial z_i}(\hat{z}(r)) \tag{11.64} \\ \hat{z}_i(\tau + \epsilon) &=& \hat{z}_i(\tau) + \epsilon\hat{r}_i\left(\tau + \frac{\epsilon}{2}\right) \tag{11.65} \\ \hat{r}_i(\tau + \epsilon) &=& \hat{r}_i\left(\tau + \frac{\epsilon}{2}\right) - \frac{\epsilon}{2}\frac{\partial E}{\partial z_i}(\hat{z}(\tau + \epsilon)) \tag{11.66} \end{eqnarray}$$

我们看到，这种方法对动量变量的更新形式是半步更新，步长为$\epsilon$，接着是对位置变量的整步更新，步长为$\epsilon$，然后是对动量变量的第二个半步更新。如果我们连续地使用几次蛙跳，那么可以看到，对动量变量的半步更新可以结合到步长为$\epsilon$的整步更新中。于是，位置变量的更新和动量变量的更新互相之间以蛙跳的形式结合。为了将动态系统跪进一个时间间隔$\tau$，我们需要进行$\tau / \epsilon$个步骤。对连续时间动态系统的离散化近似引入的误差会在极限$\epsilon \to 0$的情况下趋于0，假设函数$E(Z)$是光滑的。然而，对于实际应用中使用的一个非零的$\epsilon$，一些保留的误差仍然会存在。我们会在11.5.2节看到在混合蒙特卡罗算法中，这些误差的影响如何被消除。

总结一下，哈密顿动力学方法涉及到交替地进行一系列蛙跳更新以及根据动量变量的边缘分布进行重新采样。

注意，与基本的Metropolis方法不同，哈密顿动力学方法能够利用对数概率分布的梯度信息以及概率分布本身的信息。在函数最优化领域有一个类似的情形。大多数可以得到梯度信息的情况下，使用哈密顿动力学方法是很有优势的。非形式化地说，这种现象是由于下面的事实造成的：在$D$维空间中，与计算函数本身的代价相比，计算梯度所带来的额外的计算代价通常是一个与$D$无关的固定因子。而与函数本身只能传递一条信息相比，$D$维梯度向量可以传递$D$条信 息。