哈密顿蒙特卡洛：高效探索高维空间

哈密顿蒙特卡洛（HMC）的核心思想

哈密顿蒙特卡洛（Hamiltonian Monte Carlo, HMC）是一种结合了物理力学原理的马尔可夫链蒙特卡洛（MCMC）方法。其核心在于通过模拟物理系统的哈密顿动力学来高效探索高维概率空间。HMC利用目标分布的梯度信息，避免了随机游走的低效性，特别适用于高维复杂分布的抽样。

HMC的数学基础

HMC将目标分布 $p(\mathbf{x})$ 类比为物理系统的势能 $U(\mathbf{x})$，并引入辅助动量变量 $\mathbf{p}$ 对应动能 $K(\mathbf{p})$。系统的哈密顿量为： $$ H(\mathbf{x}, \mathbf{p}) = U(\mathbf{x}) + K(\mathbf{p}) = -\log p(\mathbf{x}) + \frac{1}{2}\mathbf{p}^T\mathbf{M}^{-1}\mathbf{p} $$

哈密顿动力学方程（时间演化）为： $$ \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 \log p(\mathbf{x}) $$

HMC的实现步骤

初始化参数

选择质量矩阵 $\mathbf{M}$（通常设为单位矩阵）
设置步长 $\epsilon$ 和轨迹长度 $L$
定义跃迁核的接受概率： $$ \alpha = \min\left(1, \exp\left(H(\mathbf{x}^{(t)},\mathbf{p}^{(t)}) - H(\mathbf{x}^,\mathbf{p}^)\right)\right) $$

动量刷新 从正态分布采样新动量： $$ \mathbf{p} \sim \mathcal{N}(0, \mathbf{M}) $$

模拟动力学 使用蛙跳积分（Leapfrog Integrator）近似哈密顿演化： $$ \mathbf{p}(t+\epsilon/2) = \mathbf{p}(t) + (\epsilon/2)\nabla \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) + (\epsilon/2)\nabla \log p(\mathbf{x}(t+\epsilon)) $$

接受/拒绝 根据哈密顿量变化决定是否接受新状态 $(\mathbf{x}^, \mathbf{p}^)$

Python代码实现框架

import numpy as np

def hmc(target_log_prob, initial_x, n_samples, step_size, n_leapfrog, mass_matrix):
    dim = len(initial_x)
    samples = np.zeros((n_samples, dim))
    x = initial_x.copy()
    
    for i in range(n_samples):
        p = np.random.multivariate_normal(np.zeros(dim), mass_matrix)
        current_x = x.copy()
        current_p = p.copy()
        
        # Leapfrog integration
        p += 0.5 * step_size * grad_log_prob(current_x)
        for _ in range(n_leapfrog):
            x += step_size * np.linalg.solve(mass_matrix, p)
            p += step_size * grad_log_prob(x)
        p += 0.5 * step_size * grad_log_prob(x)
        
        # Negate momentum for reversibility
        p = -p
        
        # Metropolis acceptance
        current_H = -target_log_prob(current_x) + 0.5 * current_p.T @ np.linalg.solve(mass_matrix, current_p)
        proposed_H = -target_log_prob(x) + 0.5 * p.T @ np.linalg.solve(mass_matrix, p)
        if np.log(np.random.rand()) < current_H - proposed_H:
            samples[i] = x
        else:
            samples[i] = current_x
            x = current_x
    
    return samples

HMC的参数调优

步长选择 需要通过试验调整步长 $\epsilon$ 使接受率保持在60%-70%之间。可以使用自适应方法如No-U-Turn Sampler（NUTS）自动调节。

轨迹长度 短轨迹导致随机游走行为，长轨迹增加计算成本。经验法则是选择轨迹长度 $L\epsilon$ 使状态在参数空间移动足够距离。

哈密顿蒙特卡洛：高效探索高维空间

哈密顿蒙特卡洛（HMC）的核心思想

HMC的数学基础

HMC的实现步骤

Python代码实现框架

HMC的参数调优

HMC的优势与局限

进阶发展方向

应用案例

性能优化建议

全站热榜

创作者周榜