跳转至

Hamiltonian dynamics for sampling on discrete spaces

作者: Raphaël Barboni, Sebastiano Grazzi, Giacomo Zanella
主题: 统计计算 / 算法
相关性: 5/10
机构绿灯: Bocconi University（US News 前 50，免分进入精读）
链接: https://arxiv.org/abs/2608.17961

一、领域脉络与小综述

这个方向是什么

这个子方向致力于为离散状态空间（如超立方体、排列空间、多元分类空间）设计高效的马尔可夫链蒙特卡洛（MCMC）采样算法。其根本挑战在于，离散空间缺乏连续空间中的微分几何结构（如梯度），使得基于梯度的算法（如哈密顿蒙特卡洛 HMC）无法直接应用。当前主流方法包括：基于连续嵌入的梯度近似、基于非可逆 lifting 的加速、以及基于分段确定性马尔可夫过程（PDMP）的连续时间采样。该方向正处于从“特例方法”向“通用框架”过渡的阶段。

发展脉络（history）

  • 奠基工作：可逆 MCMC 与 lifting 思想。Roberts & Rosenthal (2004) [62] 系统总结了 MCMC 的基础理论。Chen, Lovász & Pak (1999) [21] 提出了“lifting”概念，通过将每个状态分裂为多个副本来加速混合，并证明了混合时间最多可降至平方根量级。Neal (2004) [48] 和 Sun et al. (2010) [65] 进一步从渐近方差角度证明了非可逆链优于可逆链。
  • 主要进展：离散空间上的非可逆方法。Diaconis, Holmes & Neal (2000) [24] 分析了非可逆 Metropolis-Hastings 算法在超立方体上的加速。Zanella (2020) [67] 提出了“informed proposals”框架，通过局部平衡函数设计高效提议。Bierkens (2016) [11] 系统发展了非可逆 Metropolis-Hastings 的涡量矩阵方法。Andrieu & Livingstone (2021) [4] 建立了 skew-reversible 过程的 Peskun-Tierney 排序。Gagnon & Maire (2024) [30] 提出了渐近 Peskun 排序，用于比较 lifted 采样器。
  • 当前 frontier：PDMP 与连续时间方法。Bouchard-Côté et al. (2018) [16] 提出了 Bouncy Particle Sampler，Bierkens et al. (2019) [12] 提出了 Zig-Zag 过程，两者都是连续时间非可逆 PDMP，在连续空间上实现了弹道加速。Livingstone et al. (2025) [42] 建立了局部平衡马尔可夫跳过程的通用理论。Jansson et al. (2025) [36] 提出了通过 skew-detailed balance 构造非可逆连续时间采样器的通用机制。
  • 本文的位置：本文试图将 HMC 的加速优势（弹道 vs 扩散）推广到离散空间，通过增广连续动量变量，构造一类通用的非可逆 PDMP，并建立其不变性、遍历性、缩放极限和近似方案。它填补了“离散空间上通用非可逆 HMC 框架”的空白。

子线索聚类

  1. 基于 lifting 和 skew-reversibility 的方法：核心思想是通过增广状态空间（如引入离散动量）打破可逆性，实现加速。代表工作：Chen et al. (1999) [21], Diaconis et al. (2000) [24], Andrieu & Livingstone (2021) [4], Gagnon & Maire (2024) [30]。这些方法通常需要精心设计方向集，且当拒绝率较高时可能退化为扩散行为。
  2. 基于连续嵌入的梯度近似方法：通过将离散空间嵌入到连续空间，利用梯度信息构造提议。代表工作：Grathwohl et al. (2021) [35], Zhang et al. (2022) [68], Nishimura et al. (2020) [51], Pakman & Paninski (2013) [54]。这些方法的适用性受限于是否存在自然的连续嵌入，且梯度信息仅在间断面上有效。
  3. 基于 PDMP 的连续时间方法：利用分段确定性过程实现非可逆采样，通常需要精确模拟事件时间。代表工作：Bouchard-Côté et al. (2018) [16], Bierkens et al. (2019) [12], Andrieu et al. (2021) [3], Lu & Wang (2022) [44]。这些方法主要针对连续空间，本文将其思想反转，用于离散空间。

核心问题与瓶颈

  • 核心问题 1：如何在不依赖连续嵌入的情况下，为离散空间设计通用的非可逆动力学？
  • 核心问题 2：如何实现从扩散到弹道的加速，尤其是在目标分布具有异质性（如非对称、有约束）时？
  • 核心问题 3：如何平衡计算成本（每次迭代需要评估所有邻居的目标密度）与采样效率？
  • 当前瓶颈：现有非可逆方法（如 skew-reversible MH）在目标分布不对称时，由于拒绝率不为零，容易退化为扩散行为（O(n²) 混合时间）。而基于连续嵌入的方法受限于嵌入的自然性。

⚠️ 作者的 framing

  • 作者把缺口 frame 成：离散空间上缺乏一种通用的、不依赖连续嵌入的、能实现弹道加速的非可逆 HMC 框架。他们声称本文的离散 HMC 即使在异质性目标下也能保持弹道行为，而标准非可逆方法（如 skew-reversible MH）会退化为扩散。
  • 被淡化或回避的竞争路线：
    • 基于连续嵌入的方法（如 [51, 54]）被描述为“概念上不同”，且“适用性受限于嵌入是否自然”。作者强调自己的方法“内在定义在离散空间上”，但未深入比较在超立方体特例下，两种方法的计算复杂度（例如，[51] 的算法每次迭代只需一次目标评估，而本文的 Algorithm 1 需要 O(D) 次）。
    • skew-reversible MH 被指出在拒绝率不为零时会失去弹道行为，但作者未讨论其变体（如带自适应拒绝的方案）是否能缓解此问题。
  • 什么明显该被引 / 该存在、却没出现在 intro 里？
    • 高阶 U-统计量的计算复杂度：本文的离散 HMC 算法（Algorithm 1）需要评估所有邻居的目标密度，这本质上是一个“局部”计算。对于某些离散空间（如排列空间），邻居数量可能很大（O(n²)）。这与研究者熟悉的高阶 U-统计量的张量收缩复杂度问题有潜在联系——如何利用图结构或组合恒等式来加速这种“全邻居”评估？本文未提及。
    • 统计-计算权衡：本文的缩放极限（Theorem 6.3）暗示了 O(n) 的混合时间，但 Algorithm 1 每次迭代的计算成本是 O(D)（D 为图度数）。对于超立方体，D = n，因此总计算复杂度为 O(n²)。这与 reversible MH 的 O(n²) 混合时间相比，并无优势。作者在 Section 7.2 中提出了近似方案来降低每次迭代成本，但未从“信息-计算差距”角度分析这种权衡。这是一个值得研究者去查的张力点。

张力

未见明显对立引用。各工作主要在设定和适用场景上不同，而非结论矛盾。

二、最核心、最简单的例子 / 数学问题

第一步：符号、模型、可观测数据交代清楚

  • 符号：

    • X：离散状态空间，有限或可数。例如 X = {0,1}^n（超立方体）。
    • π：定义在 X 上的目标概率分布，我们要从中采样。
    • p ∈ R^d：连续动量变量，维度 d ≥ 1。
    • ρ：动量 p 的边际分布，通常取标准高斯分布 ρ(p) ∝ exp(-||p||²/2)。
    • ˆπ(x, p) = π(x) ρ(p)：增广空间 X × R^d 上的乘积目标分布。
    • Q(p) ∈ R^{X×X}：给定动量 p 时，状态 x 的连续时间马尔可夫链生成元矩阵。Q_{x,y}(p) ≥ 0（x≠y）是跳转率，Q_{x,x}(p) = -∑_{y≠x} Q_{x,y}(p)。
    • μ_x(p) ∈ R^d：给定状态 x 时，动量 p 的漂移函数。
    • γ ≥ 0：动量刷新率。γ = 0 表示无刷新。
    • ˆL：增广过程 (X_t, P_t) 的无穷小生成元。
    • ¯Q：一个 π-可逆的马尔可夫链生成元，作为构造 Q(p) 的基础。
    • σ_{x,y} ∈ R^d：从 x 到 y 的“方向”向量，满足 σ_{x,y} = -σ_{y,x}。
    • H_{x,y}(p)：权重函数，用于构造 Q_{x,y}(p) = ¯Q_{x,y} H_{x,y}(p)。
    • N_x：状态 x 的邻居集合（¯Q_{x,y} > 0 的 y）。
    • n_x：N_x 的大小。

  • 模型：

    • 数据生成机制：我们有一个目标分布 π 定义在离散空间 X 上。我们无法直接采样，只能通过构造一个马尔可夫链来近似。
    • 统计模型：我们构造一个连续时间马尔可夫过程 (X_t, P_t) 在 X × R^d 上，其不变分布为 ˆπ(x, p) = π(x) ρ(p)。该过程由两部分组成：
      1. 状态跳转：X_t 以速率 Q_{X_t, y}(P_t) 跳转到邻居 y。
      2. 动量漂移：P_t 在两次跳转之间按照常微分方程 dP_t = μ_{X_t}(P_t) dt 演化。
    • 已知量：π（可计算到归一化常数）、ρ（已知，如标准高斯）、¯Q（用户选择，如随机游走生成元）、σ_{x,y}（用户定义的方向）。
    • 要估的对象：π 的期望 E_π[g(X)]。

  • 可观测数据：

    • 可观测：(X_t, P_t) 的轨迹。我们可以模拟这个连续时间过程（或离散近似），得到一系列状态和动量值。
    • 想要但观测不到：π 本身。我们只能通过链的遍历平均来估计其期望。π 的归一化常数通常是未知的。

第二步：讲最小内核

最简特例：考虑最简单的离散空间——一维超立方体 X = {0, 1}，即只有两个状态。目标分布 π 任意，例如 π(0) = 0.9, π(1) = 0.1。动量是一维的 p ∈ R，取标准高斯分布 ρ(p) ∝ exp(-p²/2)。方向定义为 σ_{0,1} = 1, σ_{1,0} = -1。基础可逆生成元 ¯Q 取最简单的随机游走：¯Q_{0,1} = ¯Q_{1,0} = 1。

在这个特例下，核心思路是什么？

  1. 构造非可逆跳转率：根据 Corollary 2.2 和 Example 1（高斯动量），我们定义跳转率：

    • Q_{0,1}(p) = ¯Q_{0,1} max(0, σ_{0,1} p) = max(0, p)。即，只有当动量 p > 0 时，才能从 0 跳到 1。
    • Q_{1,0}(p) = ¯Q_{1,0} max(0, σ_{1,0} p) = max(0, -p)。即，只有当动量 p < 0 时，才能从 1 跳到 0。
    • 这实现了“动量方向决定跳转方向”的直觉。

  2. 定义动量漂移：根据 Corollary 2.2，漂移函数为常数：

    • μ_0(p) = ∑_{y≠0} ¯Q_{0,y} σ_{0,y} = 1 * 1 = 1。
    • μ_1(p) = ∑_{y≠1} ¯Q_{1,y} σ_{1,y} = 1 * (-1) = -1。
    • 这意味着，当处于状态 0 时，动量以速率 1 增加；当处于状态 1 时，动量以速率 1 减少。

  3. 过程演化：

    • 假设初始状态为 (X=0, P=0.5)。由于 P > 0，Q_{0,1}(p) > 0，存在一个正速率跳转到状态 1。同时，动量 P 以速率 μ_0 = 1 增加（dP/dt = 1）。
    • 在跳转发生前，动量线性增长：P(t) = 0.5 + t。跳转率 Q_{0,1}(P(t)) = max(0, 0.5 + t) 随时间增加。
    • 跳转时间 τ 由非齐次泊松过程决定。一旦跳转到状态 1，动量变为 P(τ)，然后开始以速率 μ_1 = -1 减少（dP/dt = -1）。此时 P(τ) > 0，所以 Q_{1,0}(P(τ)) = max(0, -P(τ)) = 0，无法立即跳回。动量必须减少到负值才能触发跳回。
    • 这个过程形成了一个确定性漂移 + 随机跳转的循环，类似于一个“反弹”运动。动量的大小决定了“惯性”，使得链倾向于沿同一方向连续移动，而不是来回震荡。

这个特例揭示了论文的核心数学思想：通过引入连续动量，将离散空间上的随机游走“提升”为一个具有惯性的非可逆过程。动量漂移 μ 和方向依赖的跳转率 Q(p) 共同作用，使得链在动量方向改变之前倾向于沿同一方向移动，从而减少回溯，实现加速。论文的一般情形（任意图、多维动量、非均匀目标）只是这个“动量驱动跳转”思想的推广，需要处理更复杂的图结构和目标分布带来的技术细节。

三、这篇论文做了什么

三句话

  1. 研究了什么问题：为离散状态空间设计一类通用的、非可逆的哈密顿蒙特卡洛（HMC）动力学，无需连续嵌入，并分析其理论性质和计算实现。
  2. 核心工具 / 方法：通过增广连续动量变量，构造分段确定性马尔可夫过程（PDMP），其生成元由方向依赖的跳转率和动量漂移组成。利用 lifting 理论和 hypocoercivity 分析进行理论刻画。
  3. 主要结论：建立了目标分布不变性的充要条件（Theorem 2.1）；证明了遍历性（Proposition 4.2）和指数收缩性（Proposition 4.3）；证明了非可逆过程在渐近方差上优于其可逆投影（Proposition 4.4）；在细粒度格点和超立方体两种高维设定下，推导了缩放极限，收敛到连续空间中的随机化 HMC 动力学，揭示了弹道加速（Theorems 5.2 & 6.3）。

关键设定与假设

  • 设定：状态空间 X 是离散的（有限或可数），配备一个 π-可逆的基础生成元 ¯Q 和一个方向集 {σ_{x,y}}。动量空间为 R^d，边际分布 ρ 具有正密度。过程 (X_t, P_t) 的生成元为 ˆL = ˆL_H + γ ˆL^D_p，其中 ˆL_H 由跳转和漂移组成，ˆL^D_p 是动量刷新（如完全刷新）。
  • 关键假设：
    • Assumption A：C_c^∞(X × R^d) 是生成元 ˆL 的一个核（core）。这是一个技术性假设，确保可以通过检验函数来验证不变性。
    • Assumption B：ˆL_rev（ˆL 的对称部分）是自伴的，且与一个闭对称 Dirichlet 型相关联。这保证了谱分析的有效性。
    • Assumption C（格点情形）：势函数 U 和 log ρ 是 C^2 且有界二阶导数；平衡函数 G 有界且 C^2；半径序列 R_n → ∞ 且 n^{-1}R_n → 0。这些保证了泰勒展开和边界效应的可控性。
    • Assumption D（超立方体情形）：log ρ 是 C^2 且有界二阶导数；平衡函数 G 有界且 C^2。这些保证了缩放极限推导中余项的可控性。
  • 相比已有文献的放宽或强化：
    • 放宽：不要求离散空间有自然的连续嵌入（对比 [51, 54]），也不要求目标分布具有对称性（对比 [24] 中分析的对称目标）。
    • 强化：要求基础生成元 ¯Q 是 π-可逆的，且方向集 σ_{x,y} 需满足反对称性。这为构造不变性提供了简洁的充分条件。

主要结果

  • Theorem 2.1（不变性条件）：给出了 ˆπ-不变性的充要条件：∇_p · (ρ μ_x) = ρ q_x（分布意义下），其中 q_x 由 Q 和 π 定义。这个条件将不变性验证转化为一个关于漂移 μ 和跳转率 Q 的偏微分方程。
  • Proposition 4.4（渐近方差排序）：对于任何 f ∈ L^2_0(ˆπ)，非可逆过程的渐近方差 Var(f, ˆL) 不大于其可逆投影 Var(f, ˆL_rev)。此外，对于仅依赖于 x 的函数 ˜f，当 γ → ∞ 时，Var(˜f ∘ τ, ˆL_rev) 趋近于基础可逆过程 ¯Q 的渐近方差，而 Var(˜f ∘ τ, ˆL) 则严格更小。这从理论上证明了非可逆性的优势。
  • Theorem 6.3（超立方体缩放极限）：对于超立方体 {0,1}^{3m} 上的约束目标 π^{(3m)}（式 36），一维动量离散 HMC 过程在时间缩放 m 下，其低维投影 (m^{-1}|X_{mt}|_2, V_{mt}) 弱收敛到连续空间 [0,1] × R 上的随机化 HMC 过程。这个极限过程具有弹道行为，其混合时间推测为 O(m)，而可逆 MH 为 O(m^2)。这个定理是本文的核心理论贡献，它严格证明了在异质性目标下，离散 HMC 仍能保持弹道加速，而标准非可逆方法（如 skew-reversible MH）会退化为扩散。

证明路线与技术技巧（理论型）

Theorem 6.3 的证明路线（3-5 步逻辑主干）：

  1. 降维与投影：利用目标 π^{(3m)} 的结构（前 m 个比特固定为 0，后 m 个比特固定为 1，中间 m 个比特均匀分布），证明过程 (X_t, P_t) 的动力学可以投影到一维变量 z_t = m^{-1}|X_t|_2 上。这是通过 Lemma A.1（lumpability）实现的，该引理给出了 Markov 链可被粗粒化的条件。
  2. 生成元计算：写出投影过程 (z_t, P_t) 的生成元 ˆL^{(m)}（式 37）。关键步骤是计算从 z 到 z ± 1/m 的总跳转率 h_m^±(z)，以及漂移项 h_m^+(z) - h_m^-(z)。这些量依赖于平衡函数 G 和目标分布的结构。
  3. 泰勒展开与缩放：对 h_m^±(z) 进行泰勒展开，利用 G 的平衡性质 G(t) = t G(1/t) 和光滑性，得到：
    • h_m^+(z) = h(z) + O(m^{-1} + m^{-2}(1-z)^{-1})
    • h_m^-(z) = h(z) + O(m^{-1} + m^{-2}z^{-1})
    • m (h_m^+(z) - h_m^-(z)) = h'(z) + O(m^{-1} z^{-1} (1-z)^{-1}) 其中 h(z) 是极限生成元中的关键函数。
  4. 应用标准缩放极限定理（Theorem 5.1）：将时间缩放 m 倍，即考虑过程 (Z_t^{(m)}, P_t^{(m)}) = (z_{mt}, V_{mt})。其生成元为 m ˆL^{(m)}。将泰勒展开结果代入，并利用测试函数 f ∈ C_c^∞([0,1] × R) 的紧支撑性，证明 m ˆL^{(m)} f 一致收敛到极限生成元 ˆL^{(∞)} f（式 38）。
  5. 边界控制：证明过程在有限时间内几乎不会触及边界 z=0 或 z=1。这通过控制边界附近的跳转次数（事件 E_m^2, E_m^3）和初始位置（事件 E_m^1）来实现，利用 Markov 不等式和期望的界。

关键跳跃点： - 从离散跳转到连续极限：最吃功夫的是证明 m (h_m^+(z) - h_m^-(z)) → h'(z) 的收敛性。这需要精确的泰勒展开，并处理 z 接近 0 或 1 时的奇异性。作者通过引入边界区域 (ε_m, m-ε_m) 并证明过程几乎不离开该区域来绕过这个困难。 - lumpability 的验证：证明投影过程 z_t 本身是 Markov 的（Lemma A.1），需要验证条件 (50) 和 (51)。这依赖于目标 π^{(3m)} 的对称性（在中间 m 个比特上可交换）和跳转率的特殊结构。

技术技巧点名： - PDMP 理论：用于定义和模拟连续时间过程。 - Lifting 理论：将离散过程提升到增广空间，实现非可逆性。 - Hypocoercivity 分析：用于证明指数收缩性（Proposition 4.3），但作者承认其给出的界可能不是最优的。 - Lumpability / 粗粒化：用于将高维过程投影到低维，简化缩放极限分析。 - 泰勒展开与余项估计：用于处理平衡函数 G 和目标比 π_y/π_x 的展开。 - 泊松过程与事件时间模拟：Algorithm 1 中利用独立非齐次泊松过程模拟跳转时间。

真实例子与应用

本文包含丰富的数值实验（Section 8），用于验证理论并展示性能： - 数据 / 场景： 1. 约束超立方体（Section 8.1, 图 2 上）：X = {0,1}^{3m}，目标 π^{(3m)} 如式 (36)，m=200。这是一个异质性目标，用于展示弹道加速。 2. 贝叶斯高斯混合模型（Section 8.1, 图 2 下）：X = [K]^n，K=2, n=500。后验分布 π 由式 (48) 定义。这是一个多元分类空间上的实际应用。 3. Plackett-Luce 排序模型（Section 8.1, 图 3）：X = Σ_n，n=200。目标 π 如式 (49)。这是一个排列空间上的应用。 - 怎么把本文方法用上去：对每个场景，定义状态空间 X、基础可逆生成元 ¯Q（如式 14 或 15）、方向 σ_{x,y}（如式 17, 19, 20）和动量分布 ρ（高斯或拉普拉斯）。然后运行 Algorithm 1（精确模拟）或其近似方案（Algorithm 2）。 - 得到什么结果： - 图 2 & 3：离散 HMC 的自相关函数衰减速度远快于基础可逆过程，轨迹探索更充分，直观展示了加速效果。 - 图 4：展示了在无刷新（γ=0）时，离散 HMC 的轨迹近似于周期性的确定性流（与 Theorem 6.3 的极限一致）。引入刷新后，周期被打破。同时比较了 Euler 和 leapfrog 离散化的误差，leapfrog 更好地保持了哈密顿量。 - 图 5：量化了 Euler 和 leapfrog + τ-leaping 近似方案的 TV 误差、动量标准差和加速比。结果显示 leapfrog 通常误差更小，且随着步长 h 增大，加速比增加，但误差也增大。 - 这个例子想说明什么：数值实验旨在验证理论预测（弹道加速），展示方法在多种离散空间上的通用性，并评估近似方案的计算-精度权衡。

🔎 结论是否比证明窄

  • 是。Theorem 6.3 严格证明的是一维动量（p ∈ R）在特定约束目标 π^{(3m)} 上的缩放极限。作者在 Section 9.2 中明确承认：“We conjecture that a scaling limit result similar to Theorem 6.3 should also hold for the discrete Hamiltonian dynamics with a multi-dimensional velocity when accelerating the process by an additional factor of O(√m). However, we were not able to come up with a proof.” 这表明论文的核心理论结论（弹道加速）在更一般的多维动量设定下仍是一个猜想，而非严格证明。数值实验（图 1）支持了 O(m^{3/2}) 的混合时间猜想，但这并非定理。
  • 此外，Proposition 4.3 给出的谱间隙下界 λ* 被作者认为“we do not believe the rate to be optimal”（Section 2.4），且与 Proposition 4.5 的下界不匹配。这表明对离散 HMC 的收敛速度的定量刻画仍不完整。

四、开放问题（点到为止，扎根具体语句）

  1. 多维动量的缩放极限：对于超立方体上的多维动量（p ∈ R^n），能否证明一个类似于 Theorem 6.3 的缩放极限，并确定加速因子（推测为 O(√m)）？扎根：Section 9.2 第一点：“We conjecture that a scaling limit result similar to Theorem 6.3 should also hold for the discrete Hamiltonian dynamics with a multi-dimensional velocity when accelerating the process by an additional factor of O(√m). However, we were not able to come up with a proof.”
  2. 最优收敛速率：能否为离散 HMC 建立最优的 L² 收敛速率（如通过 hypocoercivity 分析），使其与 Proposition 4.5 的下界匹配？扎根：Section 9.2 第二点：“Proposition 4.3 provides a contraction rate in L² for the discrete Hamiltonian dynamics which translates into an upper bound on the relaxation time. This upper bound does not match the lower bound in Proposition 4.5.”
  3. 无偏的“uninformed”近似方案：能否设计一种离散 HMC 的近似方案，每次迭代只需 O(1) 次目标评估（如通过梯度近似），同时保持目标分布不变（即无偏）并保留弹道行为？扎根：Section 9.2 第三点：“It would thus be interesting to investigate different discretization schemes of discrete Hamiltonian dynamics which are ‘uninformed’, that is, which only require O(1) target evaluations per iteration, while preserving the invariant measure and retaining a ballistic behaviour.”
  4. 方向自适应：能否自适应地调整方向 σ_{x,y} 的大小（类似于 HMC 中的质量矩阵），以减少各向异性带来的效率损失？扎根：Section 9.2 第四点：“Finally, we did not investigate adapting the magnitudes of the directions σ_{x,y}. Analogously to the mass matrix in classical HMC, such adaptation could reduce anisotropy across directions.”

Maintained by 陈星宇 · Homepage · Source on GitHub

评论