论文：并行化 Metropolis-Hastings 算法的一般结构

0摘要

马尔可夫链蒙特卡罗方法 (MCMC) 是解决许多现代统计和计算问题的重要工具；然而，一个主要的限制是这些算法固有的顺序性质。在本文中，我们提出了 Metropolis-Hastings 算法的自然推广，该算法允许使用现有的 MCMC 方法并行化单个链。我们通过并行提出多个点来做到这一点，然后在提出的点上构建和采样有限状态马尔可夫链，使得整个过程具有正确的目标密度作为其平稳分布。我们的方法通常适用且易于实施。我们演示了如何使用这种结构来大大提高各种现有 MCMC 方法的计算速度和统计效率，包括 Metropolis-Adjusted Langevin 算法和自适应 MCMC。此外，我们展示了它如何允许在哈密顿蒙特卡罗方法中使用每个积分步骤的原则方法；我们的方法增加了算法参数选择的鲁棒性，并在几乎没有额外计算成本的情况下提高了蒙特卡罗估计的准确性。

马尔可夫链蒙特卡洛| 贝叶斯推理 | 并行计算 | 哈密顿动力学

1 意义

现代统计中的许多计算问题在很大程度上依赖于马尔可夫链蒙特卡罗 (MCMC) 方法。 这些算法允许我们评估任意概率分布； 然而，由于马尔可夫特性，它们本质上是连续的，这严重限制了它们的计算速度。 我们提出了一种通用方法，允许现有 MCMC 方法的可扩展并行化。 我们通过在多个提案上定义一个有限状态马尔可夫链来做到这一点，以确保渐近收敛到正确的平稳分布。 在示例模拟中，我们展示了整体计算性能提高了两个数量级。

2 介绍

自 1970 年代推出以来，Metropolis-Hastings 算法彻底改变了计算统计 (1)。通过构建收敛到正确的平稳分布的马尔可夫链，从仅已知一个常数的任意概率分布 $\pi (X)$ 中抽取样本的能力使得贝叶斯推理能够实际应用于对大量科学现象进行建模，以及导致 Metropolis-Hastings 被认为是 20 世纪最重要的 10 种算法之一 (2)。尽管可用计算能力定期增加，但马尔可夫链蒙特卡罗 (MCMC) 算法的计算成本仍然很高；可能需要数千次迭代才能获得所需数量的低方差估计，并且通常需要为每组建议的参数评估复杂的统计模型。此外，许多 Metropolis-Hastings 算法受到其固有的顺序性质的严重限制。

已经提出了许多方法来提高 MCMC 的统计效率，尽管这些算法可以保证渐近收敛到平稳分布，但它们在有限次数的迭代中的性能可能会有很大差异。 因此，许多研究工作都集中在开发转移核，以使移动能够在远离当前点的地方提出并随后以高概率被接受，例如考虑到参数空间的相关结构 (3, 4 )，或通过使用哈密顿动力学 (5) 或扩散过程 (6)。 最近对提案内核的一项调查表明，更奇异的分布，例如 Bactrian 内核，也可能用于提高 MCMC 算法的统计效率 (7)。 高维和多模态分布的有效探索极具挑战性，并为许多进一步的方法提供了动力 (8)。

MCMC 算法的计算效率是另一个重要问题。已经提出了一些方法来利用现代计算中越来越低成本的并行性，最直接的方法是同时运行多个马尔可夫链（9）。这些可能会探索相同的分布，或相关分布的某些产品，如并行回火 (10)。此外，其他链的位置还可用于指导提议机制（11）。最近的工作将多个链的使用与自适应 MCMC 结合起来，试图使用这些多个信息源来学习适当的提案分布 (12, 13)。有时，特定的 MCMC 算法可以直接进行并行化，例如独立的 Metropolis-Hastings (14) 或切片采样 (15)，通过仔细重新参数化 (16) 或在专业硬件（例如图形处理单元）上实现的某些统计模型确实如此(GPU) (17, 18);然而，这些方法通常是针对特定问题的，并不普遍适用。对于涉及大量数据的问题，在某些情况下，并行化也可以通过在多台机器上同时使用标准 MCMC 方法对数据进行分区和分析每个子集来实现 (19)。然而，这些方法中的各个马尔可夫链都基于标准的顺序 Metropolis-Hastings 算法。

之前已经研究过使用多个提议并行化 Metropolis-Hastings 的想法；然而，这种尝试的主要缺点是缺乏计算效率。诸如 Multiple Try Metropolis (20)、Ensemble MCMC (21) 和“预取”方法 (22,23) 等算法都允许并行计算多个提议或未来提议的路径，尽管随后只有一个提议或路径被接受马尔可夫链，导致剩余点的计算浪费。另一类方法涉及将拒绝状态合并到具有适当权重的蒙特卡洛估计器中，这样得到的估计仍然是无偏的。 Frenkel (24) 对这些方法进行了概述。特别是，Tjelmeland (25) 提出的工作在每次迭代中利用了多个被拒绝的状态，并且与这里提出的想法有关，尽管最近的调查表明，包括所有被拒绝的提议实际上有时可能会导致蒙特卡罗估计与仅使用可接受的状态相比，具有渐近较大的方差（26）。

在本文中，我们提出了 Metropolis-Hastings 的概括，可用于并行化各种现有的 MCMC 算法，包括前面提到的大多数算法，从简单的随机游走 Metropolis 到最近的基于 Langevin 的算法，在欧几里德或黎曼空间上定义，以及自适应 MCMC。 我们提出的方法具有高度可扩展性，并且可以比现有算法在时间归一化效率方面提供几个数量级的改进。 该方法为哈密顿蒙特卡洛（HMC）方法提供了额外的实际好处，通过提高调整参数选择的鲁棒性并提供利用每次迭代计算的中间积分步骤的原则方法。

单个提议点的 Metropolis-Hastings

最初的 Metropolis 算法很容易理解，因为它是通过对称提议分布 (27) 实现的。 如果建议点的概率大于当前点的概率，则始终接受移动。 否则，它只被接受的概率与新点和旧点的比率成正比； 直观地说，它的访问频率较低，以确保访问的长期频率与其潜在概率成正比。 Hastings (28) 的扩展确保在使用非对称提议时收敛到正确的平稳分布； 新的接受率有效地消除了由于现在一些反向移动与正向移动的提议概率不同这一事实所带来的偏差。

令 $K(x_i,x_j)$ 是状态空间 $X$ 上的马尔可夫链，由当前点 $x_i$ 为条件的非对称提议概率分布定义。Metropolis-Hastings 算法 (29) 的一个有用解释是我们希望将马尔可夫链 $K$ 变成另一个具有平稳分布 $\pi (X)$ 的马尔可夫链。根据 Metropolis−Hastings 算法，我们提出了从 $x_i$ 到 $x_j$ 的移动，概率为 $K(x_i,x_j)$，并且以接受率为 $A(x_i,x_j)\in [0,1]$ 接受此移动。我们观察到 $A(.,.)$ 也能被视为马尔科夫链，如果它是可逆的，则收敛到所需的平稳分布，
$$\pi (x_i)K(x_i,x_j)A(x_i,x_j)=\pi (x_j)K(x_j,x_i)A(x_j,x_i)\text{\hspace{1em}(1)}$$
这是必要平衡条件成立的一个容易满足的充分条件，
$$\pi (x_i)=\int \pi (x_j))K(x_j,x_i)A(x_j,x_i)d{x}_j\text{\hspace{1em}(2)}$$
对于任意两个状态 $x_i$ 和 $x_j$,由（1）可知
$$A(x_j,x_i)=\frac{1}{R(x_i,x_j)}A(x_i,x_j)\le 1，\text{\hspace{1em}(3)}$$
其中 $R(x_i,x_j)=\frac{\pi (x_j)K(x_j,x_i)}{\pi (x_i)K(x_i,x_j)}$，由此，可逆 Metropolis−Hastings 算法 (29) 的接受概率的一般形式如下：$0\le A(x_i,x_j)\le min(1,R(x_i,x_j))$,由于 $A(x_i,x_j)\le R(x_i,x_j)$,从不等式3中可知,因此$A(x_i,x_j)$ 是我们恢复最佳 Metropolis-Hastings 接受概率（30）的最大值。

通过从转换核 $K(x_i,.)$ 中提出一个点 $x_j$ ,并以概率 $A(x_i,x_j)$ 接受，很容易看出其满足详细平衡并且链保持正确的平稳分布。 在这个阶段重要的是要认识到，当我们接受或拒绝提议的步骤时，我们实际上是从有限状态马尔可夫链 模拟。 传统的 Metropolis-Hastings 方法导致基于图 1 中所示的 两个状态马尔可夫链在每次迭代中采取一个步骤，我们注意到，由于在这个有限状态马尔可夫链中只进行了一次转换，剩余的 转移概率 $A(x_j,.)$ 不立即使用。

Metropolis-Hastings 的多个提议点

我们现在考虑在多个提议状态上定义马尔可夫链 $A(.,.)$ 的情况，并表明这也可以被构造为获得具有正确平稳分布的采样算法。推导和证明这种方法有效性的一个有用的表示是定义在产品空间上的 Metropolis-Hastings 算法。我们首先注意到联合分布 $p(x_{1:N+1})$ 可能被分解为 $N+1$ 个不同的方式 $p(x_{1:N+1})=p(x_i)p(x_{-i}|x_i)\equiv \pi (x_i)K(x_i,x_{-i})$ 其中，使用符号 $K(x_i,x_{-i})=p(x_{[1:i-1,i+1:N+1]}|x_i)$，我们观察到第 $i$ 个分解被定义为 $x_i$ 根据感兴趣的目标密度分布，而其他点 $x_{-i}$ 根据以第 $i$ 个点为条件的提议核分布。 在 Tjelmeland (25) 之后，我们可以引入我在整数 $[1:N+1]$ 上定义的离散均匀辅助随机变量，指示应该使用该联合概率分布的哪个因式分解，即 $p(x_{1:N+1},I=i)=\frac{1}{N+1}p(x_i)p(x_{-i}|x_i)$。

我们简短描述的广义 Metropolis−Hastings 算法等价于探索乘积空间 $p(x_{1:N+1},I)$ 的单个马尔可夫链； 使用两个不同的转换核的组合，每个核都保留了潜在的联合平稳分布。首先，我们更新以 $x_i$ 和 $I=i$ 为条件的变量 $x_{-i}$，这清楚地保留了正确的平稳分布，因为我们可以直接从提议内核 $K(x_i,.)$ 中采样。我们注意到我们可以自由选择这个提议内核的形式，并且在本文后面，我们考虑使用 基于朗之万扩散和哈密顿动力学的内核。其次，我们使用转移矩阵 $A(.,.)$ 对以当前状态 $x_{1:N+1}$ 为条件的辅助变量 $I$ 进行采样，其中$A(i,j)$ 定义了从 $I=i$ 到 $I=j$ 的转换概率。我们通过构造看到，当 $I=i$ 时，随机变量 $x_i$ 具有正确的目标密度 $\pi (x_i)$，这就是我们在每次迭代时收集的样本。此外，我们可以计算转移矩阵的平稳分布，记为 $A^{\infty }$，并用它直接对 $I$ 进行采样。广义 Metropolis-Hastings 算法有两个提高效率的来源；第一个是可以并行计算多个提议点的似然性，第二个是使用多个提议时提高接受率。在推导适当的转移概率并证明满足平衡条件之前，我们现在给出算法概述。广义 Metropolis-Hastings 算法进行如下。

1: 初始化开始点 ${\tilde{x}}_1$，辅助变量 $I=1$，计数 $n=0$
2: 对于每一次MCMC迭代，进行如下操作：
3：          ~~~~~~~~         以 $I$ 为条件更新 ${\tilde{x}}_{-I}$，比如从proposal kernel中抽取 $N$ 个新点 $p({\tilde{x}}_{-I}|{\tilde{x}}_I)=K({\tilde{x}}_I,.)$
4：          ~~~~~~~~         以 ${\tilde{x}}_{1:N+1}$ 为条件计算 $I$ 的平稳分布，比如 $\forall j\in [1,...,N+1],p(I=j|{\tilde{x}}_{1:N+1})=A^{\infty }(.,j)\propto \pi ({\tilde{x}}_j)K({\tilde{x}}_j,{\tilde{x}}_{-j})$,其在等式（5）中并行求得
5：          ~~~~~~~~         for m=1:N do
6:                       ~~~~~~~~~~~~~~~~~~~~~                      直接从辅助变量 $I$ 的平稳分布中采样，得到样本 $x_{n+m}={\tilde{x}}_{i.}$ 。
7:            ~~~~~~~~~~            end for
8:            ~~~~~~~~~~           更新：$n=n+N$
9: end for

我们注意到，对于 $N=1$，该算法简化为原始 Metropolis-Hastings 算法，在每次迭代中只提取一个建议点和一个样本。 在上面的算法中，我们选择了 $N$ 次采样，尽管这不一定等于提案的数量。 我们现在将推导出矩阵 $A(.,.)$的所有转移概率 使得在乘积空间上满足平衡条件。 我们首先注意到更新这个产品空间中的变量 $I$ 的细致平衡条件是
$$\frac{1}{N+1}\pi (x_i)K(x_i,x_{-i})A(i,j)=\frac{1}{N+1}\pi (x_j)K(x_j,x_{-j})A(j,i)\text{\hspace{1em}(4)}$$
对于所有 $i$ 和 $j$，平衡条件如下
$$\frac{1}{N+1}\pi (x_i)K(x_i,x_{-i}）=\sum _{j=1}^{N+1}\frac{1}{N+1}\pi (x_j)K(x_j,x_{-j})A(j,i)\text{\hspace{1em}(5)}$$
对于所有 $i$。 我们现在可以为转移概率推导出一个类似的表达式，就像我们在单个提案案例中所做的那样。

Proposition 1:

我们可以构造一个有限状态马尔可夫链 $A(.,.)$ 在给定当前状态集 $x_{1:N+1}$ 的情况下定义 $I$ 的转移概率
$$A(i|j)=\{\begin{array}{ll}\frac{1}{N}\min (1,R(i,j)),& ifj\ne j\\ 1-{\sum }_{j\ne i}A(i,j),& otherwise\end{array}\text{\hspace{1em}(6)}$$
其中
$$R(i,j)=\frac{\pi (x_j)K(x_j,x_{-j})}{\pi (x_i)K(x_i,x_{-i})}\text{\hspace{1em}(7)}$$

马尔可夫链的每一个转变，$A(i,.)$，满足细致平衡和平衡条件（方程 4 和 5）。;当使用此类更新对变量 $I$ 进行采样时， 联合分布 $p(X_{1:N+1},I)$是不变的。

从每对值的可逆性条件 $[I=i,I=j]$，我们有 $0\le A(i,j)\le min(1,R(i,j)),(j\ne i)$，其中 $R$ 如等式（7）所示，为了$A(.,.)$ 成为马尔可夫链，我们也要求 ${\sum }_{j=1}^{N+1}A(i,j)=1$ 并且 $A(i,i)=1-{\sum }_{j\ne i}A(i,j)\ge 0$。当 $N>1$ 时， 我们可能会尝试满足 $A(.,.)$ 的这些要求，通过考虑权重 $c_j\in [0,1]$ 使得 $0\le A(i,j)\le c_j\min (1,R(i,j))(j\ne i)$。如果我们选择 $A(i,j)$ 来最大化这个不等式，那么等式 4 中的可逆性条件意味着 $c_j=c_i$，即我们可以为所有提议的点选择一个常数 $c$。 因此我们有 ${\sum }_{j\ne i}c\min (1，R(i,j))\le 1$。取这个量的最大值意味着 ${\sum }_{j\ne i}c\le 1$，因此让 $c=1/N$ 满足这个不等式。 因此，现在很清楚 $A(.,.)$，在方程式(6)中定义， 它是一个马尔可夫链，满足等式（4）中的细致平衡条件。 因此也是等式（5）中的平衡条件， 从中直接得出 $A(.,.)$的平稳分布。

我们注意到接受率取决于当前点和所有提议的点，因此从 $I=i$ 转换到 $I=j$ 的概率可能会变得非常小。 出于这个原因，我们可以使用参考文献中介绍的类似方法 25。由此我们引入一个辅助变量 $z$ 并提出形式为 $\tilde{K}(x_i,x_{-i})=K(x_i,z)K(z,x_{-i})$的建议。 然后接受率简化为
$$R(i,j)=\frac{\pi (x_j)\tilde{K}(x_j,x_{-j})}{\pi (x_i)\tilde{K}(x_i,x_{-i})}=\frac{\pi (x_j)K(x_j,z)K(z,x_i)}{\pi (x_i)K(x_i,z)K(z,x_j)}\text{\hspace{1em}(8)}$$
并且，对于对称建议，接受率进一步简化为 $R(i,j)=\pi (x_j)/\pi (x_i)$。

3.材料和方法

我们现在详细介绍我们用于研究 Metropolis-Hastings 的这种广义构造的数值性能的 MCMC 算法、统计模型和效率度量。 对于模拟，我们采用贝叶斯方法，使得我们的目标密度是后验分布，$\pi (x)=p(x|y)\propto p(y|x)p(x)$。 因此，我们可以使用 Metropolis-Hastings 从后验分布中抽取样本，而无需明确计算由边际似然 $p(y)$ 给出的归一化常数。

Metropolis-Adjusted Langevin 算法
收敛到正确平稳分布的连续时间随机微分方程可用于在 MCMC 算法中进行有效移动；但是，仍然需要 Metropolis-Hastings 校正步骤，因为离散解不再收敛到正确的分布。这种方法效果很好，因为它考虑了目标密度的局部几何形状。 Metropolis-Adjusted Langevin 算法 (MALA) 的原始版本基于在欧几里得空间 (6) 中定义的扩散，尽管最近表明它们也可以在统计模型 (4) 诱导的黎曼流形上定义。在黎曼情况下，生成的提议机制具有由预期的 Fisher 信息定义的特定于位置的协方差矩阵，它满足度量张量 (31) 的属性；因此，生成的马尔可夫链通过在每一步考虑统计模型对其参数的微小变化的平均局部敏感性来提出转换。可以使用以下简化的流形 MALA (SmMALA) 转换核 (4) 提出计算有效的建议，我们将其用于本文中的后续数值模拟，$K(x_i,.)=N(x_i+{ϵ}^2/2)G(x_i{)}^{-1}{\nabla }_x\log p(x_i|y),{ϵ}^{2}G(x_i{)}^{-1})$，其中 $G(x_i)$ 是由预期的 Fisher 信息 (4) 给出的特定于位置的协方差矩阵。

自适应 MCMC
另一种方法是根据先前接受的 MCMC 移动自适应地学习适当的提议协方差结构。对于此类算法，必须单独证明收敛到平稳分布，因为链不再是马尔可夫（3）。出于我们模拟研究的目的，我们使用参考文献中详述的基于朗之万的自适应 MCMC 算法。 32. 转换核为
其中 D是有界漂移函数，δ 是固定常数，在每次迭代中计算的协方差矩阵 Λ 基于先前接受的具有递减适应的移动 (32)，使得链渐近收敛到正确的平稳分布。这种自适应 MCMC 算法可以用所提出的广义 Metropolis-Hastings 方法直接实现，方法是在每次迭代后简单地更新协方差矩阵，即在 N 个点被提出和采样之后。

HMC

略

计算采样效率

由于马尔可夫属性，使用 MCMC 生成的样本表现出一定程度的自相关。 评估此类方法的统计效率的一种方法是计算有效样本量（ESS），即从收集的后验样本总数中有效独立样本的数量。 对于每个协变量，我们可以使用以下标准度量，$ESS=N[1+2{\sum }_{k}r(k){]}^{-1}$，其中 $N$ 是后验样本的数量，${\sum }_{k}r(k)$ 是 $K$ 个单调样本自相关的总和，由初始单调序列估计估计器（34）。 另一种方法是考虑均方跳跃距离，$MSJD=(1/N){\sum }_{n=1}^{N-1}||x_{n+1}-x_n|{|}^2$，这与测量滞后 1 自相关有关，还可用于优化提案的参数 内核（35）。 我们利用这两种方法来研究广义 Metropolis-Hastings 的性能。

普通微分方程模型
基于常微分方程 (ODE) 系统的统计模型具有广泛的科学用途，从描述细胞调控网络到大规模化学过程，并且在此类机械模型中的推断提出了许多挑战，这都是由于涉及的计算工作量评估每组建议参数的解决方案，以及从潜在的强相关后验分布中采样的挑战，这是由模型方程中存在的任何非线性引起的。我们考虑 $\dot{u}=f(u,\theta )$ 形式的 ODE，我们对隐式定义的解状态 $u$ 进行数值求解。给定每个解状态 $y_i=u_i+{\gamma }_i$ 的不确定测量，根据所选误差模型定义 $\gamma$，我们希望推断模型参数$\theta$ 的后验分布。我们使用 Fitzhugh-Nagumo 模型作为这种系统 (4) 的说明性示例，该系统由两个状态组成，$\dot{V}=c(V-V^3/3+R)$ 和 $\dot{R}=-(V-a+bR)/c$。以下参考。如图 4 所示，我们使用 Fitzhugh−Nagumo ODE 模型在 $t=0$ 和 $t=20$ 之间生成的 $200$ 个数据点，模型参数 $a=0.2、b=0.2$ 和 $c=3$，初始条件 $V_0=-1$和 $R_0=1$。然后将 SD 等于 0.5 的高斯分布误差添加到数据中，随后用于推断模型参数。

逻辑回归模型
我们还考虑了一个由 $N\times D$ 设计矩阵 $X$ 定义的贝叶斯逻辑回归示例，其中 $N$ 是样本数，每个样本都使用 $D$ 协变量进行描述，使得 $p(y=1|X,\theta )=\sigma (X\theta )$ 和 $p(y=0|X,\theta )=1-\sigma (X\theta )$，其中 $\sigma$ 表示逻辑函数。 我们通过定义高斯先验 $p(\theta )=N(0,\alpha I)$，$\alpha =100$ 对回归系数 $\theta \in R^D$ 进行推断。 为了实现 SmMALA，该模型的度量张量直接如下： $G(\theta )=X^T\Lambda X+\alpha I$，其中 $\Lambda$ 是一个对角矩阵，元素为 $\Lambda (n,n)=\sigma ({\theta }^T{X}_{(n,\cdot )}^T)(1-\sigma ({\theta }^T{X}_{(n,\cdot )}^T)$; 见参考文献4 了解更多详情。

结果
统计效率随着提案数量的增加而增加。 我们使用越来越多的提案来研究广义 Metropolis-Hastings 的统计效率。 一旦选择了 $N$ 个提议，就可以在单个处理器上的 $N$ 个内核上或在 $N$ 个单独的计算机处理器上并行计算似然性和转移概率。 这导致每个样本的计算速度大约增加 N 倍，尽管我们牢记，由于通信开销会损失一些计算效率，这当然会受到算法的硬件和软件特定细节的影响执行。 最大的改进将来自复杂的数学模型，其解决方案的计算成本相对于通信费用而言是昂贵的。

4.结论

我们研究了 Metropolis-Hastings 算法的自然扩展，它允许通过多个提议在现有 MCMC 算法中直接并行化单个链。这是一个有点令人惊讶的结果，因为 Metropolis-Hastings 的主要限制之一是其固有的顺序性。我们研究了它与黎曼、自适应和哈密顿建议的使用，并展示了计算和统计效率的改进。

在计算资源有限的情况下，经常会出现这样一个问题：单条链运行时间更长，还是多条链运行时间更短。文献中给出的理论论据表明，马尔可夫链的单个较长运行是可取的（34），并且在这种情况下，所提出的方法将具有明显的价值。然而，我们注意到即使在多链设置中，Generalized Metropolis-Hastings 应该仍然有用，因为由于似然计算的并行化和降低的拒绝率，该方法能够为所有单独的链提供额外的效率改进。

广义 Metropolis-Hastings 直接适用于广泛的现有 MCMC 算法，并有望在利用 HMC 中的整个集成路径、降低所得蒙特卡罗估计器的整体方差和提高调整的鲁棒性方面特别有价值参数，以及加速贝叶斯推理对计算成本高昂的更复杂的数学模型。使用在多个提议点上定义的有限状态马尔可夫链来满足详细平衡的想法为算法设计提供了更大的灵活性。我们预计这种非常通用的方法将导致进一步的方法学发展，以及未来更有效和可扩展的并行 MCMC 方法。