MAMHMC：微正则哈密顿蒙特卡洛如何用Metropolis修复能量壳偏差

最近在方法列表里看到一个名字：Metropolis Adjusted Microcanonical Hamiltonian Monte Carlo。长，拗口，读起来像在报口令。把它拆开其实是一个很自然的组合：一边是微正则系综里的确定性哈密顿流，一边是经典Metropolis接受规则，两者拼起来解决一个长期被忽视的问题——当微正则哈密顿蒙特卡洛遇上有限精度数值积分时，如何保证采样仍然精确。这篇博文想把这套设计的动机、数学结构、易错点以及实际调参时会遇到的问题讲透。适合已经接触过HMC、想在采样器层面做点探索的人，也适合被名字吓退、但好奇“微正则”到底在采样里意味着什么的朋友。

先说结论：MAMHMC不是普通HMC的竞争者，而是它的一个变体加固版。普通HMC靠每步刷新动量来保证遍历性，MAMHMC则把总能量固定到某个随机抽取的壳上，跑一段确定性可逆的哈密顿流，再通过Metropolis步骤把数值误差导致的壳偏移“打回去”。理解这套东西，关键不在于记住接受率公式，而在于想清楚能量壳、壳上测度和数值误差三者之间的真实关系。

1. 从一个“绕口的名字”说起：它到底保了哪两样保险

1.1 普通HMC的日常：刷新动量、跑轨迹、接受或拒绝

所有用过HMC的人应该都有这个印象：每次迭代要做三件事，重新从高斯分布里抽一组动量p，带着这组p跑L步leapfrog，最后按能量差接受或拒绝终点。其中“刷新动量”这一步最不起眼，却承担了最重要的职责——它让链有机会改变动能方向、跳出当前势能井，是随机性的主要来源。

代价也很直接：每步的起始点p都从零开始，轨迹之间几乎没有衔接，链在相空间里表现出明显的“跳跃感”。如果你把连续两次迭代的轨迹画出来，会发现它们像一堆随机方向的短线，而不是一条平滑延伸的路径。这在目标分布比较崎岖的高维空间里会造成较大的方差，也让每一步携带的几何信息被冲掉不少。

1.2 微正则思路：把“随机冲撞”换成“能量分层”

微正则HMC（MCHMC）想换一种方式提供随机性。它不再每一次都重新采一个动量，而是先随机抽一个总能量E，再把动量限制在满足动能+势能=E的状态集合上。这个集合在相空间里是一个高维曲面，通常叫能量壳（energy shell）。然后沿着哈密顿方程跑一段确定性轨迹，由于哈密顿流本身保持能量，这个轨迹会一直贴在能量壳上，不会跑出去。

这样一来，随机性只剩两个来源：每步开始前抽一个E，以及在能量壳上随机选一个动量方向。链的路径变得光滑、可逆、确定性强，理论上每一步都无需接受拒绝。名字里加上的Metropolis Adjusted，意思很直白：因为现实中的数值积分器并不能让能量真正守恒，轨迹端点会从能量壳上漏出去，所以需要补一道Metropolis阀门，把漏出去的端点概率性地拦下来。

2. 微正则HMC怎么把“能量”变成采样引擎

2.1 给定能量E，如何从当前q构造出合格的p

这是MCHMC和普通HMC最明显的分岔点。普通HMC里p直接按标准高斯采样，MCHMC则必须先解一个约束：势能U(q)和动能K(p)加起来要等于目标总能量E。给定当前q，可用动能是E-U(q)，于是动量p的方向随机，幅值必须满足一个关于质量矩阵的二次型约束。

最常见的做法是：先抽取一个各向同性的单位方向，再按2(E-U(q))开根号去缩放。这样构造出来的点(q,p)必然落在能量壳Σ_E上。这一步很朴素，但也是整个方法里最容易出bug的地方——如果E比U(q)还小，开根号直接无意义，状态就是一个空壳，后面会专门聊这个坑。

2.2 确定性流为什么能完成壳上的混合

状态到了壳上之后，接下来的L步leapfrog是一个标准的辛映射。辛映射有两个重要性质：可逆和时间方向对称。同一个轨迹反着跑，能精确回到起点；同时它保持相空间体积不变。再加上哈密顿流本身与能量梯度正交，轨迹会一直停留在同一能量壳上。

这听起来像“只在一个曲面上转圈，怎么可能混合得好”。但注意，壳上的运动并不是规整的环绕，而是受势能梯度影响、在壳内蜿蜒前进。只要轨迹足够长、势能形态不是完全可积到锁死状态，它就能在壳上近似遍历。也就是说，给定E之后，理论上能通过确定性路径在壳上“均匀”地逛一遍。这个“均匀”是指相对于壳上诱导的Liouville测度而言。

2.3 能量层必须正确配重，否则拼不回正则分布

链不能永远停留在同一个E上，否则只能在某个能量壳的切片里采样，得到的是微正则条件分布，和目标的正则分布π(q)∝exp(-U(q))不是一回事。MCHMC的解法是每次迭代抽取一个新的E，问题就变成：E应该怎么抽，才能让各个壳的贡献叠加起来正好等于目标分布？

这背后是一个经典的coarea公式思想。目标联合分布密度可以理解成一摞能量壳的加权平均，每个壳的权重既包含指数衰减因子exp(-E)，也包含该壳在相空间中的“面积”增长率。如果直接把E抽成均匀分布，大E的壳因相空间体积巨大而天然占优，q的边际会严重偏向高能量区域。所以实际中E的边缘分布通常取Gamma分布这一类，形状参数与系统自由度相关，用来对冲壳面积随E高次增长的趋势。具体实现时，很多代码会取动量维数的一半作为形状参数，但一定要结合质量矩阵和哈密顿量的标度重新核一遍，不能照抄。

2.4 一个常被误解的点：固定能量不等于冻结温度

我见过不少人把微正则HMC理解成“让系统在低温下跑，所以只探索低势能区域”。这是错的。固定E只是固定了动能与势能之和，并不限制两者的分配比例。同一个E壳上，既可能动能很大、势能很低，也可能势能很大、动能接近于零。一条长为L的轨迹会在这两种模式之间切换，所以它遍历的是同一个总能量下各种“冷热配置”的混合，而不是局限于势能极小点。

换句话说，普通HMC里的温度通过动量高斯分布的方差来控制，而MCHMC里没有单独的温度参数，只有一个总能量门槛。这个区别会在多峰问题里带来一个致命后果：如果E太低，轨迹确实翻不过高势垒。对应关系是“能不能翻山”取决于当前壳的E是否大于山高，而不是取决于随机抽到的动量够不够猛。

3. 能量壳上的数值裂缝：真正需要维修的对象

3.1 leapfrog的守恒是“近似”而不是“恒等”

任何辛积分器都不可能让连续哈密顿量精确不变。leapfrog在每步只做到二阶局部近似，单步能量误差量级通常在O(ε³)，这里ε是步长。长期运行时，轨道作为辛映射不会出现单调的能量漂移，但会围绕真实等能面做有界振荡。这个“有界”依赖步长足够小，一旦步长放大，振荡幅度就以多项式级别膨胀，轨迹末端离理论壳的距离肉眼可见。

MCHMC在理想世界里可以百分之百接受每个轨迹端点，因为理想动力学把点保持在壳上。但数值世界里，每一步的误差都让端点相对理论壳产生偏移，链的平稳分布不再是一个干净的能量壳约束，而是壳附近一个被模糊过的带状区域。更麻烦的是，这个带状区域的“厚度”在不同位置不一样，壳面积大的地方偏差暴露得更充分，最终污染到q的边际分布。

3.2 为什么传统HMC几乎没有这个烦恼

标准HMC其实也面对同样的leapfrog误差，但它没有把能量当作必须守住的硬约束。每步结束时，它会计算H_new和H_old的差值，并用exp(-(H_new-H_old))做接受概率。也就是说，误差被显式地转换成了接受/拒绝的权重，变成了算法的一部分。尽管大步长下接受率会下降，目标分布始终是精确的，不会因为积分器不完美而悄悄改变形状。

MCHMC本来为了效率去掉接受步骤，靠的是“动力学绝对可逆且保测”这个理想条件。一旦回到浮点数世界，等于拆掉了安全网，任何数值误差都直接进入样本分布。Metropolis Adjustment要做的，就是把这张安全网重新装上，同时尽量保留微正则流带来的平滑路径优势。

3.3 壳脱离如何把偏差注入样本：一个直观图像

可以把能量壳想象成一层层很薄的球壳，理论链只允许在某个壳上移动。数值积分让每一步像“脚滑了一下”，状态从壳上掉到了壳与壳之间的空隙里。如果径直接受这些掉出壳的点，就等价于在空隙区域也赋予了概率质量，等于给目标分布做了平滑卷积。光滑的偏差在低维可能感受不到，但维度越高，壳的面积增长越疯狂，同样的能量误差在切向上被放大得越明显，偏差就不再是“一小点噪声”，而是实实在在的分布变形。

4. Metropolis校正怎么在能量壳上正确落地

4.1 不能直接套exp(-ΔH)：三个原因

很多人看到“Metropolis Adjusted”，第一反应是把终点和起点的H差值丢进exp里，像标准HMC一样算接受率。这在MCHMC语境里是不对的。原因有三：第一，MCHMC里的p不是按高斯先验独立采样的，它的幅值被E约束住，直接套exp(-H)等于对同一个信息重复计费；第二，目标分布支撑在余维一的能量壳上，不是完整相空间里的正则分布，Lebesgue密度意义已经改变；第三，轨迹端点的能量误差主要来自积分器，不是来自提议核的性质，用完整能量差会过度惩罚，导致接受率低到没有实用价值。

4.2 壳上测度、切向体积与接受率的真实含义

正确做法是把目标分布写成约束在壳上的微正则分布，并考虑壳上诱导的测度。哈密顿流保持壳上的辛体积元，这意味着在精确动力学下，壳上的均匀分布是严格平稳的。所谓Metropolis调整，实际是在问：数值轨迹这一端相对理论壳偏离了这么多，若把它当作一个从旧点到新点的可逆提议，那么按壳上目标测度来看，新点是否该比旧点更不受欢迎？

如果我们坚持让目标始终是“固定E下的壳上均匀分布”，那么任何有限精度下产生的壳外状态都必须经过“壳外惩罚”或“投影回壳”等额外操作。接受概率里会涉及壳法向方向的局部体积因子，也就是|∇H|在终点和起点的比值。这个因子在标准HMC里被完整高斯动量的密度比悄悄抵消了，在MCHMC里却必须显式写出来，否则接受率会系统性偏向能量误差为正或为负的方向。

4.3 反向轨迹的角色：确定性流自带可逆提议

Metropolis调整的巧妙之处在于，确定性哈密顿流本身是可逆的。前向从(q,p)沿轨迹走到(q′,p′)，反向从那一点沿同一轨迹走回来，概率完全确定。这和随机提议有本质区别：不需要为反向核的密度比多操心，Jacobian由辛性质给出为1。因此唯一需要计算的就是壳上目标密度的比值，以及脱离壳所带来的权重修正。

这也就解释了为什么这个组合能work：理想情况下，壳上密度是常数，接受率为1，Metropolis阀门完全不干扰确定性路径；有数值误差时，阀门只负责矫正那些明显“出格”的状态，保留其余大部分光滑轨迹。

4.4 理想状态无限精度下的退化景观

在无限精度的思维实验里，MAMHMC退化为普通MCHMC，接受概率恒等于1，链完全确定性地在壳上跑。这种退化并不丢人，反而揭示出方法的核心哲学：Metropolis调整不是给微正则采样新增的负担，而是对现实计算环境的一种补偿。

从这个角度去读论文里的算法框图，很多细节都会豁然开朗。你会看到他们先构造能量壳上的状态，跑一段确定性流，再按壳上概率体积比做一个简单判断。它既不像标准HMC那样把每一步都变成随机跳跃，也不像原版MCHMC那样对误差装作看不见，而是把握在一个中间位置。

5. 调参视角：能量随机化、步长与多峰陷阱

5.1 E分布的形状参数不是随意定的

这是工程上最值得花时间的旋钮。前面说过，E的边缘分布需要和壳面积增长匹配。理想公式涉及势能函数的整体几何，不能只靠维度数一拍脑袋。我给一个可操作的启动方式：先用目标分布的一个粗略样本估计势能的典型范围，再让E的分布覆盖这个范围的若干倍标准差；运行一小段之后，检查接受率随E的分布曲线，如果高能区的接受率明显下降，就要收窄E的右尾或调低形状参数，反之则加宽。

另一个容易被忽视的点：E的抽样不能和当前状态完全解耦。如果当前q恰好落在一个势能很低的深井里，一个很低的E也能构造动量，轨迹就被锁在井旁边的壳上；如果E抽得太高，轨迹又会跑到势能很大的区域烧掉大量动能。这就需要E分布在不同势能区间都有足够的质量，否则链会长期盘踞在某个“舒适壳”附近。

5.2 步长和轨迹长度的作用方式不同于HMC

标准HMC里，步长主要影响接受率，轨迹长度主要影响去相关速度。在MAMHMC里，步长同时影响两层效果：一是leapfrog带来的壳偏离大小，二是Metropolis护门员的工作负荷。步长一放大，壳偏离方差成多项式增长，接受率会以指数衰减。所以实用的起步方式是把接受率目标放在0.6到0.8之间，然后用双误差条或简单二分法搜索步长。

轨迹长度的设置反而可以比HMC更保守。由于轨迹沿壳是连贯的，长度太长只是在同一个壳上绕很多圈，对跨壳混合没有帮助。我倾向于让轨迹长度适应壳的形状复杂度：目标分布是强相关高斯时，果壳拉得长一些有助于降低相邻样本自相关；遇到多峰或强烈非线性势能时，过长轨迹只会增加壳偏离的累积，拉低接受率。

5.3 能量屏障和壳连通性：Metropolis不背的锅

多峰分布是MAMHMC最需要小心的场景。假设两个模式之间隔着一个势垒ΔU，只有当抽样到的E大于当前所在的势垒高度时才可能穿越。如果E分布右尾太薄，链会在一个壳上徘徊很长时间，看起来基线已经收敛，实际上从未见过另一个峰。Metropolis调整只能保证单个壳内部的分布正确，无法帮助连通不同的壳段。

应对办法有两种方向。一种是刻意加宽E分布的尾，让高E壳被抽中的概率不至于太低；另一种是把E更新改成带惰性的随机游走，在相邻迭代中通过小步改变E，让它有机会缓慢爬过势垒对应的能量区间。前者简单粗暴但减小有效样本量，后者更细腻但需要额外引入E的转移核设计。无论哪种，都要在收敛诊断里专门看E的轨迹是否覆盖了足够宽的范围。

5.4 和标准HMC/NUTS放在一起选型时的直觉

我不建议在任何场景下无脑替换NUTS。NUTS的自适应步长和树状轨迹在强相关高斯类问题上极其可靠，而且绝大多数概率编程框架里开箱即用。MAMHMC的优势场景更偏向需要保留连续轨迹信息、目标较为光滑、能量壳连通性较好、且你能接受手动调E分布的情况。

6. 从零搭一个概念验证工程骨架（含三个实测坑）

6.1 核心循环骨架

如果只是想验证思路，不需要复杂的自动微分框架。下面这段Python风格骨架反映了MAMHMC的核心循环，注意它省略了严格的壳上接受率推导，只用于展示状态结构和检查点。

import numpy as np def leapfrog_flow(q, p, grad_U, M_inv, eps, L): p_half = p - 0.5 * eps * grad_U(q) for _ in range(L - 1): q = q + eps * M_inv @ p_half p_half = p_half - eps * grad_U(q) q = q + eps * M_inv @ p_half p = p_half - 0.5 * eps * grad_U(q) return q, p def energy_on_shell(q, p, M_inv): return 0.5 * p @ M_inv @ p + U(q) # U自行实现 def sample_momentum_on_shell(q, E, M_inv, U): available = E - U(q) if available <= 0: return None direction = np.random.randn(len(q)) norm_dir = np.sqrt(direction @ M_inv @ direction) direction = direction / norm_dir p = direction * np.sqrt(2.0 * available) return p def mamhmc_step(q, grad_U, M_inv, energy_dist, eps, L): E = energy_dist() p = sample_momentum_on_shell(q, E, M_inv, U) if p is None: return q, 0.0 # 实际实现必须按空壳概率加权 q_old, p_old = q, p q_new, p_new = leapfrog_flow(q, p, grad_U, M_inv, eps, L) # 这里需要按壳上测度设计接受率，下面只是占位示意 shell_dev = abs(energy_on_shell(q_new, p_new, M_inv) - E) accept = min(1.0, np.exp(-shell_dev * scale)) if np.random.rand() < accept: return q_new, shell_dev return q_old, shell_dev

注意这段代码里我故意没写正规的接受率。实际复现时，你需要根据原论文的壳上概率体积推导替换那几行占位逻辑，否则就是一个玩具。

6.2 工程坑一：负动能与空能量壳

最容易撞到的错误是available动能变成负数。原因往往是步长太大，轨迹钻进了势能极高的狭窄区域，导致下次抽样E时构造不出动量。如果你在循环里简单跳过这一步，等价于对E分布做条件截断，样本不再忠于目标分布。更隐蔽的是，有些实现为了省事会在E过低时重新抽E，直到available为正，这也会制造偏差。

应对方式是在算法框架内显式处理空壳事件：要么把E的下界与当前q的势能挂钩，并把这种截断写进概率密度；要么在采样q的同时保证不会进入势能高于E支撑的区域。调试时，遇到大量负动能，先怀疑步长而非E分布。

6.3 工程坑二：E序列的自相关才是真正的瓶颈

只看q的自相关时间会给你虚假的安全感。我曾经跑过一个模型，q的ESS看起来很高，但E的轨迹几乎卡在一个窄带上，很长时间不更新。原因是轨迹长度太短，每次独立抽E之后，链还没来得及跨过壳间的距离就被下一次E更新拽到新壳上。结果链只是在很多能量壳的“门口”打转，深层结构没有被探索。

解决方案是把E的时间序列也加入诊断，画一下自相关图。如果E的滞后相关系数长期高于0.5，就加大L或把E更新改成惯性更强的随机游走，让链在同一壳上多停留一会儿再换层。

6.4 工程坑三：能量误差报告是调参的第一信号

MAMHMC调试时最该打印的不是接受率，而是每步的壳偏离量shell_dev。这个指标告诉你积分器是否有能力把轨迹保持在该壳上。如果平均壳偏离已经超过了0.1而接受率还在硬撑，大概率是框架中的接受规则在强行“拽住”大量状态，导致有效步长比名义值小得多。这时候直接减小ε比调任何接受率系数都管用。

反过来，如果壳偏离很小但接受率仍然低，问题一定出在能量壳测度的体积因子处理上，回到理论推导逐项检查，不要继续盲目调参。

我自己的体会是，MAMHMC这类方法最迷人的地方不在“要不要接受步骤”，而在它逼着你去思考能量壳这个抽象对象在数值计算里的真实样貌。标准HMC让你轻松地忽略积分误差，MCHMC又让你天真地信任积分精度，Metropolis Adjusted版本恰好提供了一个反思的支点：我们既想要确定性轨迹的平滑优势，又不想丢掉概率上的严谨性。如果你有条件在自定义分布上做实验，我很建议拿它和NUTS做一轮逐项对比——先看E序列，再看q的边际，最后再比ESS/s，你对采样器的理解会整个换一层。