马尔可夫链蒙特卡洛（MCMC）中的哈密顿蒙特卡洛（HMC）算法原理与采样过程
字数 2576 2025-11-07 12:33:00

马尔可夫链蒙特卡洛（MCMC）中的哈密顿蒙特卡洛（HMC）算法原理与采样过程

题目描述
哈密顿蒙特卡洛（HMC）是一种结合了物理系统动力学思想的MCMC采样算法，用于从复杂的高维概率分布中高效采样。它通过引入动量变量，模拟哈密顿动力学轨迹，显著降低了随机游走行为，比传统Metropolis-Hastings算法收敛更快。本题要求理解HMC的物理类比、动力学方程构建、蛙跳积分方法，以及接受概率的计算过程。

1. 物理类比与概率映射

  • 目标分布：假设需要从目标分布 \(p(\mathbf{x})\) 采样，其中 \(\mathbf{x} \in \mathbb{R}^d\) 是位置变量。
  • 动量变量：引入人工动量变量 \(\mathbf{p} \in \mathbb{R}^d\)，服从高斯分布 \(p(\mathbf{p}) = \mathcal{N}(\mathbf{0}, \mathbf{M})\)，其中 \(\mathbf{M}\) 为质量矩阵（通常设为单位矩阵）。
  • 联合分布：定义哈密顿函数 \(H(\mathbf{x}, \mathbf{p}) = -\log p(\mathbf{x}) + \frac{1}{2} \mathbf{p}^\top \mathbf{M}^{-1} \mathbf{p}\)，其中 \(-\log p(\mathbf{x})\) 对应势能 \(U(\mathbf{x})\)，动能 \(K(\mathbf{p}) = \frac{1}{2} \mathbf{p}^\top \mathbf{M}^{-1} \mathbf{p}\)。
  • 核心思想：采样联合分布 \(p(\mathbf{x}, \mathbf{p}) \propto \exp(-H(\mathbf{x}, \mathbf{p}))\)，边际分布即目标 \(p(\mathbf{x})\)。

2. 哈密顿动力学方程
哈密顿系统由以下微分方程描述：

\[\frac{d\mathbf{x}}{dt} = \frac{\partial H}{\partial \mathbf{p}} = \mathbf{M}^{-1} \mathbf{p}, \quad \frac{d\mathbf{p}}{dt} = -\frac{\partial H}{\partial \mathbf{x}} = \nabla_{\mathbf{x}} \log p(\mathbf{x}). \]

  • 动量更新依赖目标分布梯度，梯度信息引导采样方向，避免随机游走。
  • 动力学轨迹保持哈密顿量 \(H(\mathbf{x}, \mathbf{p})\) 守恒（理想情况下）。

3. 离散化与蛙跳积分
实际求解需离散化时间，常用蛙跳法（Leapfrog Integrator）保持数值稳定性：
设步长为 \(\epsilon\)，迭代 \(L\) 步：

  1. 半步动量更新：

\[ \mathbf{p}(t + \epsilon/2) = \mathbf{p}(t) + \frac{\epsilon}{2} \nabla_{\mathbf{x}} \log p(\mathbf{x}(t)). \]

  1. 全步位置更新：

\[ \mathbf{x}(t + \epsilon) = \mathbf{x}(t) + \epsilon \mathbf{M}^{-1} \mathbf{p}(t + \epsilon/2). \]

  1. 剩余半步动量更新：

\[ \mathbf{p}(t + \epsilon) = \mathbf{p}(t + \epsilon/2) + \frac{\epsilon}{2} \nabla_{\mathbf{x}} \log p(\mathbf{x}(t + \epsilon)). \]

  • 重复 \(L\) 次得到候选状态 \((\mathbf{x}^*, \mathbf{p}^*)\)。

4. 接受-拒绝步骤
由于离散化误差，哈密顿量不严格守恒，需按Metropolis准则接受候选状态：

\[\alpha = \min\left(1, \frac{\exp(-H(\mathbf{x}^*, \mathbf{p}^*))}{\exp(-H(\mathbf{x}, \mathbf{p}))}\right) = \min\left(1, \exp(H(\mathbf{x}, \mathbf{p}) - H(\mathbf{x}^*, \mathbf{p}^*))\right). \]

  • 若接受，新状态为 \((\mathbf{x}^*, \mathbf{p}^*)\)；否则保留原状态。
  • 动量在每次采样前重新随机化，确保遍历性。

5. 算法流程总结

  1. 初始化 \(\mathbf{x}_0\)，设置步长 \(\epsilon\)、步数 \(L\)、样本数 \(N\)。
  2. 对每个样本 \(i\)：
    • 从 \(\mathcal{N}(\mathbf{0}, \mathbf{M})\) 采样新动量 \(\mathbf{p}\)。
    • 通过 \(L\) 步蛙跳积分从 \((\mathbf{x}_i, \mathbf{p})\) 计算候选 \((\mathbf{x}^*, \mathbf{p}^*)\)。
    • 以概率 \(\alpha\) 接受 \(\mathbf{x}_{i+1} = \mathbf{x}^*\)，否则 \(\mathbf{x}_{i+1} = \mathbf{x}_i\)。
  3. 返回样本序列 \(\{\mathbf{x}_1, \dots, \mathbf{x}_N\}\)。

关键参数影响

  • 步长 \(\epsilon\)：过大导致拒绝率高，过小增加计算成本。
  • 步数 \(L\)：需平衡探索效率与路径周期性（避免回归起点）。
  • 调参方法：可基于平均接受率（如目标0.65）自适应调整 \(\epsilon\)。

HMC通过物理动力学有效探索高维空间，尤其适用于具有连续密度函数的复杂分布采样。

马尔可夫链蒙特卡洛（MCMC）中的哈密顿蒙特卡洛（HMC）算法原理与采样过程 题目描述 哈密顿蒙特卡洛（HMC）是一种结合了物理系统动力学思想的MCMC采样算法，用于从复杂的高维概率分布中高效采样。它通过引入动量变量，模拟哈密顿动力学轨迹，显著降低了随机游走行为，比传统Metropolis-Hastings算法收敛更快。本题要求理解HMC的物理类比、动力学方程构建、蛙跳积分方法，以及接受概率的计算过程。 1. 物理类比与概率映射 目标分布 ：假设需要从目标分布 \( p(\mathbf{x}) \) 采样，其中 \(\mathbf{x} \in \mathbb{R}^d\) 是位置变量。 动量变量 ：引入人工动量变量 \(\mathbf{p} \in \mathbb{R}^d\)，服从高斯分布 \( p(\mathbf{p}) = \mathcal{N}(\mathbf{0}, \mathbf{M}) \)，其中 \(\mathbf{M}\) 为质量矩阵（通常设为单位矩阵）。 联合分布 ：定义哈密顿函数 \( H(\mathbf{x}, \mathbf{p}) = -\log p(\mathbf{x}) + \frac{1}{2} \mathbf{p}^\top \mathbf{M}^{-1} \mathbf{p} \)，其中 \( -\log p(\mathbf{x}) \) 对应势能 \( U(\mathbf{x}) \)，动能 \( K(\mathbf{p}) = \frac{1}{2} \mathbf{p}^\top \mathbf{M}^{-1} \mathbf{p} \)。 核心思想 ：采样联合分布 \( p(\mathbf{x}, \mathbf{p}) \propto \exp(-H(\mathbf{x}, \mathbf{p})) \)，边际分布即目标 \( p(\mathbf{x}) \)。 2. 哈密顿动力学方程 哈密顿系统由以下微分方程描述： \[ \frac{d\mathbf{x}}{dt} = \frac{\partial H}{\partial \mathbf{p}} = \mathbf{M}^{-1} \mathbf{p}, \quad \frac{d\mathbf{p}}{dt} = -\frac{\partial H}{\partial \mathbf{x}} = \nabla_ {\mathbf{x}} \log p(\mathbf{x}). \] 动量更新依赖目标分布梯度，梯度信息引导采样方向，避免随机游走。 动力学轨迹保持哈密顿量 \( H(\mathbf{x}, \mathbf{p}) \) 守恒（理想情况下）。 3. 离散化与蛙跳积分 实际求解需离散化时间，常用 蛙跳法 （Leapfrog Integrator）保持数值稳定性： 设步长为 \(\epsilon\)，迭代 \( L \) 步： 半步动量更新 ： \[ \mathbf{p}(t + \epsilon/2) = \mathbf{p}(t) + \frac{\epsilon}{2} \nabla_ {\mathbf{x}} \log p(\mathbf{x}(t)). \] 全步位置更新 ： \[ \mathbf{x}(t + \epsilon) = \mathbf{x}(t) + \epsilon \mathbf{M}^{-1} \mathbf{p}(t + \epsilon/2). \] 剩余半步动量更新 ： \[ \mathbf{p}(t + \epsilon) = \mathbf{p}(t + \epsilon/2) + \frac{\epsilon}{2} \nabla_ {\mathbf{x}} \log p(\mathbf{x}(t + \epsilon)). \] 重复 \( L \) 次得到候选状态 \( (\mathbf{x}^ , \mathbf{p}^ ) \)。 4. 接受-拒绝步骤 由于离散化误差，哈密顿量不严格守恒，需按Metropolis准则接受候选状态： \[ \alpha = \min\left(1, \frac{\exp(-H(\mathbf{x}^ , \mathbf{p}^ ))}{\exp(-H(\mathbf{x}, \mathbf{p}))}\right) = \min\left(1, \exp(H(\mathbf{x}, \mathbf{p}) - H(\mathbf{x}^ , \mathbf{p}^ ))\right). \] 若接受，新状态为 \( (\mathbf{x}^ , \mathbf{p}^ ) \)；否则保留原状态。 动量在每次采样前重新随机化，确保遍历性。 5. 算法流程总结 初始化 \( \mathbf{x}_ 0 \)，设置步长 \(\epsilon\)、步数 \( L \)、样本数 \( N \)。 对每个样本 \( i \)： 从 \( \mathcal{N}(\mathbf{0}, \mathbf{M}) \) 采样新动量 \( \mathbf{p} \)。 通过 \( L \) 步蛙跳积分从 \( (\mathbf{x}_ i, \mathbf{p}) \) 计算候选 \( (\mathbf{x}^ , \mathbf{p}^ ) \)。 以概率 \( \alpha \) 接受 \( \mathbf{x} {i+1} = \mathbf{x}^* \)，否则 \( \mathbf{x} {i+1} = \mathbf{x}_ i \)。 返回样本序列 \( \{\mathbf{x}_ 1, \dots, \mathbf{x}_ N\} \)。 关键参数影响 步长 \(\epsilon\) ：过大导致拒绝率高，过小增加计算成本。 步数 \( L \) ：需平衡探索效率与路径周期性（避免回归起点）。 调参方法 ：可基于平均接受率（如目标0.65）自适应调整 \(\epsilon\)。 HMC通过物理动力学有效探索高维空间，尤其适用于具有连续密度函数的复杂分布采样。