Hamiltonian Monte Carlo几何理论

1 引言

哈密顿蒙特卡洛（Hamiltonian Monte Carlo, HMC）是现代MCMC方法中最重要的高效采样算法之一。¹ 与简单的随机游走采样相比，HMC利用几何信息（梯度）来指导采样，在高维复杂分布上表现出色。

本章从几何视角深入解析HMC的理论基础，包括哈密顿动力学、辛几何、以及黎曼流形上的采样。

2 哈密顿动力学基础

2.1 物理背景

HMC的灵感来自统计力学。在哈密顿力学中，粒子在势能场中运动：

势能 $U (q)$ ：定义粒子位置 $q$ 的能量
动能 $K (p)$ ：定义动量 $p$ 的能量
总能量 $H (q, p) = U (q) + K (p)$

2.2 哈密顿方程

哈密顿方程描述系统的演化：

\frac{d q}{d t} = \frac{\partial H}{\partial p} = \frac{\partial K ( p )}{\partial p}

\frac{d p}{d t} = - \frac{\partial H}{\partial q} = - \frac{\partial U ( q )}{\partial q}

物理意义：

位置的变化率 = 动量对时间的导数
动量的变化率 = 受力的负梯度（来自势能）

2.3 与MCMC的联系

将哈密顿动力学与贝叶斯推断联系起来：

位置 $q$ = 待采样参数
势能 $U (q) = - lo g π (q)$ （负对数后验）
动能 $K (p) = \frac{1}{2} p^{T} M^{- 1} p$ （通常取二次型）
总能量 $H (q, p) = - lo g π (q) + \frac{1}{2} p^{T} M^{- 1} p$

关键洞察：哈密顿动力学的轨迹在参数空间上探索得更远，避免随机游走的低效性。

3 HMC算法详解

3.1 基本算法

HMC算法流程：

给定当前状态 (q, p):
1. 采样初始动量 p ~ N(0, M)
2. 运行L步Leapfrog积分器模拟哈密顿动力学
3. 以概率 min(1, exp(ΔH)) 接受最终状态

3.2 Leapfrog积分器

Leapfrog（跳蛙）积分器是辛几何的数值方法，比标准欧拉法更稳定。

算法：

def leapfrog(q, p, epsilon, L, M, grad_U):
    """
    Leapfrog积分器
    q: 位置
    p: 动量
    epsilon: 步长
    L: 步数
    M: 质量矩阵
    grad_U: 势能梯度
    """
    p = p - 0.5 * epsilon * grad_U(q)  # 半步更新动量
    
    for i in range(L):
        q = q + epsilon * np.linalg.solve(M, p)  # 全步更新位置
        if i < L - 1:
            p = p - epsilon * grad_U(q)  # 全步更新动量
    
    p = p - 0.5 * epsilon * grad_U(q)  # 半步更新动量
    
    return q, p

为什么是Leapfrog：

交替更新位置和动量
每步只需要一次梯度计算
保持辛几何的”体积保持”性质

3.3 接受率分析

HMC的接受率取决于积分误差：

α = min (1, exp (- Δ H))

其中 $Δ H = H (q^{*}, p^{*}) - H (q, p)$ 。

理想情况：辛积分器精确保持 $H$ ，接受率为1。

实际情况：数值误差导致 $H$ 略有增加，接受率略低于1。

3.4 参数选择

参数	选择原则
步长 $ϵ$	目标接受率 65-90%
步数 $L$	使轨迹长约1-2个单位长度
质量矩阵 $M$	通常取单位阵或协方差估计

自适应调整：NUTS算法自动选择 $L$ 。

4 辛几何视角

4.1 辛流形基础

哈密顿动力学的几何框架是辛流形。²

辛形式：在 $(q, p)$ 空间中，辛形式定义为：

ω = i \sum d q^{i} \land d p_{i} = d q \land d p

基本性质：

$ω$ 是非退化的： $ω^{n} \neq = 0$
$ω$ 是闭合的： $d ω = 0$

4.2 辛映射的体积保持

定理：哈密顿动力学产生的映射是辛映射，自动保持辛形式 $ω$ 。

推论：辛映射保持相空间体积（Liouville定理）。

几何意义：

辛积分器（Leapfrog）保持 $ω$
标准欧拉法不保持 $ω$ ，导致相位误差累积

4.3 保持性质的优势

辛几何积分器相比普通数值方法的优势：

方面	普通积分器	辛积分器
能量守恒	短期振荡	长期稳定
体积保持	否	是
相位误差	累积	有界振荡
长期行为	漂移	稳定

4.4 可视化理解

考虑二维相位空间 $(q, p)$ ：

                    p
                    ↑
                    |     *  (最终位置)
                    |   *
                    | *
                    |/
            --------+------------→ q
                   /
                  /
                 *
                *

辛积分器保持相空间的几何结构，因此即使能量有微小波动，轨迹也不会”螺旋”到奇点。

5 黎曼几何视角

5.1 黎曼流形上的HMC

当目标分布的协方差结构随位置变化时，使用黎曼度量可以提高效率。³

黎曼HMC（RMHMC）：

H (q, p) = U (q) + \frac{1}{2} lo g ∣Σ (q) ∣ + \frac{1}{2} p^{T} Σ (q)^{- 1} p

其中 $Σ (q)$ 是位置依赖的协方差矩阵。

5.2 Fisher信息度量

最常用的黎曼度量是Fisher信息矩阵：

G (q) = I (q) = E [\nabla_{q} lo g p (x ∣ q) \nabla_{q} lo g p (x ∣ q)^{T}]

几何意义：

度量随位置自适应变化
在高曲率区域移动较慢
在低曲率区域移动较快

5.3 梯度信息利用

普通HMC：只用梯度 $\nabla U (q)$

黎曼HMC：用梯度 $\nabla U (q)$ 和 Fisher信息 $I (q)$

效率提升：在高度非均匀的分布上，RMHMC比普通HMC效率高数倍。

5.4 近似黎曼HMC

精确计算Fisher信息在高维时很昂贵。近似方法：

对角近似： $G (q) \approx diag (I_{1}, \dots, I_{d})$
块对角近似：分组变量共享度量
随机逼近：使用随机梯度估计

6 NUTS算法

6.1 No-U-Turn准则

NUTS（No-U-Turn Sampler）自动选择Leapfrog步数 $L$ 。⁴

U-Turn条件：

(q_{0, L} - q_{0, - L})^{T} p_{0, L} < 0 或 (q_{0, L} - q_{0, - L})^{T} p_{0, - L} < 0

物理意义：当轨迹开始”折返”时停止，避免重复采样。

6.2 NUTS算法

def nuts sampler(U, grad_U, epsilon, max_treedepth=10):
    """
    NUTS采样器
    """
    # 初始化
    theta = current_state
    epsilon = find_reasonable_epsilon(theta, U, grad_U)
    
    # 构建树
    theta_minus = theta_plus = theta
    rho_minus = rho_plus = momentum
    
    j = 0
    n = 1
    s = 1
    
    while s == 1 and j < max_treedepth:
        # 随机选择方向
        v = random.choice([-1, 1])
        
        if v == -1:
            theta_minus, rho_minus, _, _, = build_tree(
                theta_minus, rho_minus, v, epsilon, j)
        else:
            _, _, theta_plus, rho_plus, = build_tree(
                theta_plus, rho_plus, v, epsilon, j)
        
        # 接受条件
        if n > 1:
            alpha = min(1, N_subtree / n)
            if random.uniform(0, 1) < alpha:
                theta = choose_uniform_node()
        
        # U-Turn检查
        s = check_u_turn(theta_minus, theta_plus, rho_minus, rho_plus)
        
        n += N_subtree
        j += 1
    
    return theta

6.3 自适应参数调整

NUTS通常配合自适应步长调整：

预热期：运行若干迭代估计接受率
步长调整：目标接受率约65%
质量矩阵估计：使用样本协方差

6.4 收敛性保证

NUTS的理论保证：

无偏采样（渐近）
自动选择合适的轨迹长度
在高维问题上表现稳健

7 实用算法实现

7.1 完整HMC实现

import numpy as np
from scipy import stats
 
class HamiltonianMonteCarlo:
    """哈密顿蒙特卡洛采样器"""
    
    def __init__(self, log_posterior, gradient, n_samples, warmup=500,
                 epsilon=None, L=None, M=None):
        self.log_posterior = log_posterior
        self.gradient = gradient
        self.n_samples = n_samples
        self.warmup = warmup
        self.epsilon = epsilon
        self.L = L
        self.M = M if M is not None else np.eye(d)
        
    def find_reasonable_epsilon(self, theta0):
        """寻找合理的初始步长"""
        epsilon = 1.0
        p = np.random.randn(len(theta0))
        theta = theta0.copy()
        
        # 尝试不同的epsilon值
        for _ in range(20):
            theta_new, p_new = self.leapfrog(theta, p, epsilon)
            
            # 检查接受率
            alpha = min(1, np.exp(
                self.log_posterior(theta_new) - self.log_posterior(theta)
                - 0.5 * p_new @ np.linalg.solve(self.M, p_new)
                + 0.5 * p @ np.linalg.solve(self.M, p)
            ))
            
            if 0.25 < alpha < 0.75:
                break
            elif alpha > 0.75:
                epsilon *= 2
            else:
                epsilon /= 2
        
        return epsilon
    
    def leapfrog(self, theta, p, epsilon):
        """Leapfrog积分器"""
        d = len(theta)
        
        # 半步动量更新
        p = p - 0.5 * epsilon * self.gradient(theta)
        
        # 全步位置更新
        theta = theta + epsilon * np.linalg.solve(self.M, p)
        
        # 半步动量更新
        p = p - 0.5 * epsilon * self.gradient(theta)
        
        return theta, p
    
    def sample(self, theta0):
        """执行采样"""
        d = len(theta0)
        samples = []
        
        # 初始化
        theta = theta0.copy()
        
        # 预热期
        if self.epsilon is None:
            self.epsilon = self.find_reasonable_epsilon(theta)
        if self.L is None:
            self.L = int(1.0 / self.epsilon)
        
        for i in range(self.warmup + self.n_samples):
            # 采样动量
            p = np.random.multivariate_normal(
                np.zeros(d), self.M)
            
            # 存储旧状态
            theta_old = theta.copy()
            p_old = p.copy()
            
            # Leapfrog积分
            theta_star, p_star = self.leapfrog(theta, p, 
                                               self.epsilon * np.ones(self.L))
            
            # Metropolis接受
            alpha = min(1, np.exp(
                self.log_posterior(theta_star) - self.log_posterior(theta)
                - 0.5 * p_star @ np.linalg.solve(self.M, p_star)
                + 0.5 * p @ np.linalg.solve(self.M, p)
            ))
            
            if np.random.uniform() < alpha:
                theta = theta_star
            else:
                theta = theta_old
            
            # 存储样本
            if i >= self.warmup:
                samples.append(theta)
        
        return np.array(samples)

7.2 应用案例：高斯混合模型

def demo_hmc_gmm():
    """HMC在高斯混合模型上的应用"""
    
    # 定义双峰分布
    def log_posterior(x):
        mean1, mean2 = -3, 3
        std = 1.0
        p1 = stats.norm.pdf(x, mean1, std)
        p2 = stats.norm.pdf(x, mean2, std)
        return np.log(0.5 * p1 + 0.5 * p2)
    
    def gradient(x):
        mean1, mean2 = -3, 3
        std = 1.0
        p1 = stats.norm.pdf(x, mean1, std)
        p2 = stats.norm.pdf(x, mean2, std)
        grad = (mean1 - x) / std**2 * p1 + (mean2 - x) / std**2 * p2
        return grad / (p1 + p2)
    
    # 运行HMC
    hmc = HamiltonianMonteCarlo(log_posterior, gradient, 
                                n_samples=5000, warmup=500)
    samples = hmc.sample(np.array([0.0]))
    
    return samples

8 与其他采样方法的对比

8.1 与Metropolis-Hastings对比

方面	MH算法	HMC
提议分布	任意对称分布	哈密顿动力学
接受率	变化大	通常较高
步长选择	敏感	较鲁棒
梯度利用	否	是
高维效率	低	高

8.2 与Gibbs采样对比

方面	Gibbs	HMC
条件分布	需要解析可采样	不需要
梯度信息	否	是
变量选择	顺序更新	同时更新
混合速度	依赖条件独立	依赖梯度

8.3 选择指南

问题特征
   │
   ├─ 维度低 (<10) ──> Gibbs / MH
   │
   ├─ 梯度可计算
   │    ├─ 分布高度非均匀 ──> RMHMC
   │    └─ 其他 ──> HMC/NUTS
   │
   └─ 梯度不可计算 ──> 切片采样 / 变分推断

9 理论基础

9.1 细致平衡条件

HMC满足MCMC的细致平衡条件：

π (x) K (x \to x^{'}) = π (x^{'}) K (x^{'} \to x)

其中 $K$ 是转移核。

证明思路：

Leapfrog积分精确保持体积
接受步骤满足可逆性
动量翻转确保细致平衡

9.2 遍历定理

定理：对于满足条件的分布 $π$ ，HMC生成的链是：

不可约：能访问所有状态
非周期：包含连续部分
正常返：期望返回时间有限

9.3 混合时间分析

HMC的混合时间（mixing time）取决于：

分布曲率：平坦区域混合快
步长选择：太大不稳定，太小效率低
质量矩阵：应接近协方差

10 与相关内容的联系

10.1 概率图模型

HMC是贝叶斯网络推理的重要工具：

bayesian-networks — 贝叶斯网络基础
junction-tree-algorithm — 精确推断（对比）

10.2 贝叶斯深度学习

HMC可用于训练贝叶斯神经网络：

bayesian-neural-networks — BNN基础
mcmc-methods — MCMC方法回顾

10.3 变分推断

当HMC太慢时，变分推断是替代方案：

variational-inference-deep-dive — 变分推断深度
variational-inference-advanced — 进阶方法

11 总结

本章从几何视角深入解析了HMC算法：

哈密顿动力学基础：物理对应、数值实现
辛几何视角：体积保持、能量守恒的数学本质
黎曼几何视角：自适应度量、Fisher信息
NUTS算法：自动参数调整
实现细节：完整代码和应用案例

HMC代表了几何信息在采样中的应用，是现代MCMC方法的重要里程碑。

参考文献

Neal, R.M. (2011). MCMC using Hamiltonian dynamics. Handbook of Markov Chain Monte Carlo, 54(2), 113-162. ↩
Betancourt, M. (2017). A conceptual introduction to Hamiltonian Monte Carlo. arXiv:1701.02434. ↩
Girolami, M. & Calderhead, B. (2011). Riemann manifold Langevin and Hamiltonian Monte Carlo methods. Journal of the Royal Statistical Society, 73(2), 123-214. ↩
Hoffman, M.D. & Gelman, A. (2014). The No-U-Turn Sampler: Adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research, 15(1), 1593-1623. ↩

Metaphor

探索

Hamiltonian Monte Carlo几何理论

Hamiltonian Monte Carlo几何理论

1 引言

2 哈密顿动力学基础

2.1 物理背景

2.2 哈密顿方程

2.3 与MCMC的联系

3 HMC算法详解

3.1 基本算法

3.2 Leapfrog积分器

3.3 接受率分析

3.4 参数选择

4 辛几何视角

4.1 辛流形基础

4.2 辛映射的体积保持

4.3 保持性质的优势

4.4 可视化理解

5 黎曼几何视角

5.1 黎曼流形上的HMC

5.2 Fisher信息度量

5.3 梯度信息利用

5.4 近似黎曼HMC

6 NUTS算法

6.1 No-U-Turn准则

6.2 NUTS算法

6.3 自适应参数调整

6.4 收敛性保证

7 实用算法实现

7.1 完整HMC实现

7.2 应用案例：高斯混合模型

8 与其他采样方法的对比

8.1 与Metropolis-Hastings对比

8.2 与Gibbs采样对比

8.3 选择指南

9 理论基础

9.1 细致平衡条件

9.2 遍历定理

9.3 混合时间分析

10 与相关内容的联系

10.1 概率图模型

10.2 贝叶斯深度学习

10.3 变分推断

11 总结

参考文献

Footnotes

关系图谱

目录

反向链接