使用并行计算改进基于独立 Metropolis-Hastings 的估计

0.摘要

在本文中，我们考虑了如果提议的值独立于马尔可夫链的当前值，则可以通过通用 Metropolis-Hastings 算法利用并行原始功率这一事实的含义。特别是，我们提出了对独立 Metropolis-Hastings 算法的改进，该算法显着降低了从 MCMC 输出派生的任何估计量的方差，计算成本为零，因为这些改进是基于可以并行产生的固定数量的目标密度评估。本文开发的技术不会危及算法的马尔可夫收敛特性，因为它们基于 Gelfand 和 Smith (1990) 的 Rao-Blackwell 原理，已在 Casella 和 Robert (1996)、Atchadé 和Perron (2005)，以及 Douc 和 Robert (2011)。我们在玩具正态示例和经典概率回归模型上说明了这些改进，但强调它们适用于独立 Metropolis-Hastings 适用的任何情况。本文中使用的代码可作为补充材料。
关键词：加速； MCMC算法；并行化；排列;随机发生器；Rao-Blackwell。

1.介绍

Metropolis-Hastings (MH) 算法提供了一种迭代和收敛方案来从复杂的目标密度 $\pi$ 中进行采样。 该算法的每次迭代都会生成一个马尔可夫链的新值，该值依赖于前一次迭代的结果。 基本的马尔可夫原理很容易理解，并导致了一个通用的收敛原理，例如 Robert 和 Casella (2004) 所描述的。 然而，由于其马尔可夫性质，该算法并不直接并行化，这在 R 等较慢的语言中造成困难（R Development Core Team 2006）。 然而，以非常低的成本提供的并行内核数量的增加推动了对“并行友好”算法的越来越多的兴趣，即对可以从标准计算机上可用的并行处理单元中受益的算法（参见 ，例如，Lee 等人 2009；Suchard 等人 2010；Holmes 等人 2011）。

除了并行独立运行 $p$ 个 MCMC 算法并合并结果的基本方案之外，已经使用不同的技术来增强通用 Metropolis-Hastings (MH) 算法的某种程度的并行性。例如，一个自然的入口是依靠马尔可夫链的更新属性（Mykland, Tierney, and Yu 1995; Robert 1995; Hobert et al. 2002），等待所有 $p$ 链显示更新事件，然后使用块同 iid，但无法消除马尔可夫的约束。

Rosenthal (2000) 还指出了考虑老化时间的难题：虽然对于单次 MCMC 运行，老化时间基本上可以忽略不计，但在运行并行链时确实会产生明显的偏差（除非完美采样可以实施）。 Craiu 和 Meng (2005) 将对偶耦合和分层与完美抽样相结合。使用不同的方法，Craiu、Rosenthal 和 Yang (2009) 依靠 p 个平行链来构建自适应 MCMC 算法，本质上考虑到目标密度在链上的乘积是它们的目标，这种观点显然会影响收敛多链的属性。 Corander、Gyllenberg 和 Koski (2006) 利用并行化构建了一种不可逆算法，可以避免特定邻域结构的缩放效应，因此专注于一种非常特殊的问题。

MH 算法的一个特殊家族是独立的 Metropolis-Hastings (IMH) 算法，其中提议分布（以及因此提议的值）不依赖于马尔可夫链的当前状态。 由于这个特性，这个特定的算法更容易并行化，因此可以被认为是高效并行马尔可夫链蒙特卡罗算法的一个很好的构建块，如第 2 节所述。我们将专注于计算似然函数构成 MH 算法执行时间的主要部分。 Wraith 等人的文章中提供了这种设置的最现实示例(2009)。该模型基于一个非常复杂的 Fortran 程序，该程序转换了几个宇宙学实验的结果，因此对计算时间的要求很高。 在这个模型中，Wraith 等人 (2009) 使用自适应重要性采样和大规模并行化，而不是 MCMC。

文章的计划如下：在第 2 节中回顾标准 IMH 算法，然后描述我们的改进方案，这里称为“block Independent Metropolis-Hastings”（block IMH）。 这种改进取决于第 3 节中详细描述的 { 1 , . . . , p } \{1,..., p\} {1,...,p} 上的排列选择。我们在第 4 节中演示了块 IMH 和 Rao-Blackwellization 之间的联系。一个玩具示例的结果贯穿始终文章和一个实际的概率回归示例在第 5 节中描述，作为该方法的说明。

2. 改进 IMH 算法

2.1 标准 IMH 算法

我们在这里回顾了基本的 IMH 算法，假设我们可以采样的提议分布的可用性，并且在归一化常数之前，概率密度 $\mu$ 是已知的。 算法 1 中描述的独立 Metropolis-Hastings 算法生成与目标分布相对应的具有不变密度 $\pi$ 的马尔可夫链。

在 Monte Carlo 和 MCMC 算法的大图中，IMH 算法具有相当特殊的地位,与其他 MCMC 方案（Robert and Casella 2004）相比，它确实得到了更多的研究，但与更通用的随机游走 Metropolis-Hastings 算法相比，它无疑是一个不太实用的解决方案。例如，它本身很少使用，因为它需要推导对真实目标的相当好的近似，而这种近似通常是不可用的。另一方面，一阶近似和Metropolis-within-Gibbs 方案对于要求基于目标的高斯表示的IMH 局部移动并不陌生。 IMH算法的理论研究比比皆是，因为它与重要性采样等非马尔可夫模拟方法有很强的联系。与随机游走 Metropolis-Hastings 方案相反，IMH 算法可能具有非常好的收敛特性，并且还可能达到接近 1 的接受概率。此外，本文开发的并行化方案提供的方差减少潜在的巨大收益与随机游走 Metropolis-Hastings 算法相比，可以抵消原始 IMH 的较低效率。

在处理并行性时，IMH 算法的一个重要特征是它只能以迭代的方式工作，因为需要步骤 $t$ 的结果，即值 $x_t$ 来计算步骤 $t+1$ 的接受率。这种顺序构造对于算法的验证是强制性的，因为它的核心是马尔可夫属性（Robert and Casella 2004）。尽管如此，鉴于在 IMH 算法中，提议的值 $(y_t)$ 是独立于马尔可夫链的当前状态 $x_t$ 生成的，完全可以设想首先生成 $T$ 个提议的值 $y_t$，以及计算相关比率 ${\omega }_t=\pi (y_t)/\mu (y_t)$ 的变化。一旦这个计算要求完成，只需要迭代地考虑验收步骤。当 $y_t$ 的模拟和 ${\omega }_t$ 的推导可以并行实现时，这种两步观点极大地节省了计算时间，因为比率 $\rho (x_{t-1},y_t)$ 的剩余计算鉴于 ${\omega }_t$ 及其随后与均匀绘制的比较通常要快几个数量级。

在这方面，IMH 算法与随机游走 Metropolis-Hastings (RWMH) 算法有很大不同，后者无法预先处理接受率，因为建议的模拟值取决于马尔可夫链的当前值。因此，并行处理方案的普遍可用性可能会导致 IMH 算法再次流行起来。事实上，当利用 $p$ 个并行处理单元时，一个 IMH 可以运行 $p$ 倍于 RWMH 的迭代，而计算成本几乎相同，因为 RWMH 不能直接并行化。

为了更好地描述这种增加的计算能力，我们首先注意到，一旦产生了马尔可夫链的 $T$ 个连续值，该序列通常被处理为常规蒙特卡罗样本，以获得目标分布下的期望近似值， $E\pi [h(X)]$ 说，对于一些任意函数 $h$。我们在本文中提出了一种技术，该技术通过利用并行处理单元而不危及马尔可夫特性来提高这种期望的估计精度。

在介绍我们的改进方案之前，我们为表示 IMH 算法的单个步骤的运算符引入符号 $\vee$（读作“或”）。使用这种表示法，给定当前值 $x_t$ 和一系列 $p$ 独立建议值 $y_1,...,y_p\sim \mu$，IMH 算法根据图 1 中的图表从步骤 $t$ 进入步骤 $t+p$。

2.2 块 IMH 算法

我们建议充分利用模拟的建议值及其相应 $\omega$ 比率的计算。为此，我们介绍了块 IMH 算法，它由大小为 $p\times p$ 的块的连续模拟组成。在这个替代方案中，块的数量 $b$ 使得期望的迭代次数 $T$ 等于 $b\ast p$，以保持与标准 IMH 输出的比较公平。通常 $p$ 不需要校准，因为它表示代码可以利用的物理并行处理单元的数量。然而，原则上，这个数字 $p$ 可以任意设置，并且基于虚拟并行处理单元，缺点是计算成本增加。 （注意，块 IMH 算法也可以在没有并行能力的情况下实现；但它仍然提供了方差增益，可以抵消时间的增加。）在下面的例子中，我们取 $p$ 从 $4$ 到 $100$ 变化。我们首先解释如何模拟一个块，然后如何从一个块移动到下一个块。

图 2. 从步骤 $t_1$ 到步骤 $t_p$ 的模块仿真。 此处，建议值的循环排列用于说明目的。

一个 $p\times p$ 块包括$p$ 个马尔可夫链值并行产生的$p$个子代，所有这些都从当前状态 $x_t$ 的时间 $t$ 开始，并且都基于相同的建议模拟值 $y_1,...,y_p$。 $p$ 流之间的差异在于包含这些 $y_i$ 的顺序。 例如，这些顺序可能是 $y_1,...,y_p$ 的 $p$ 个循环排列，或者它们可能是随机排列，如第 3 节中详细讨论（和比较）的那样。块 IMH 算法如图 2 所示循环排列的集合。

应该清楚的是，在这个块中找到的 $p$ 个平行链中的每一个在单独取时都是长度为 $p$ 的有效 MCMC 序列。 因此，它可以作为常规 MCMC 输出进行处理。 特别是，如果 $x_t$ 是从平稳分布模拟的，则任何后续 $x_{t+j}^{(j)}$ 也都是从平稳分布模拟的。 但是，$p$ 个并行流的点是双倍的：

• 它旨在整合由选择 $y_k$ 的辅助顺序产生的部分随机性，接近 Perron (1999) 提倡的 $y_k$ 顺序统计的条件；
• 它还旨在部分整合在选择过程中生成统一变量所产生的随机性，因为块实现导致绘制 $p^2$ 个统一实现而不是标准IMH 设置的 $p$ 个统一实现。

这两点本质上相当于实现了一种新的 Rao-Blackwellization 技术（在第 4 节中绘制了更精确的联系）。 在独立设置中，每个 $y_k$ 在 $p$ 个平行链的 $p$ 个步骤中出现 $n_k\ge 0$ 次，也就是说，对于 $p^2$ 个实现。 因此，当考虑 $E\pi [h(X)]$ 的标准估计量 ${\hat{\tau }}_1$ 时，基于单个 MCMC 链，
$${\hat{\tau }}_1(x_t,y_{1:p})=\frac{1}{p}\sum _{k=1}^{p}h(x_{t+k})$$

该估计量必然具有比双平均更大的方差，
$${\hat{\tau }}_2(x_t,y_{1:p})=\frac{1}{p}\sum _{j=1}^p\sum _{k=1}^{p}h(x_{t+k}^{(j)})=\frac{1}{p^2}\sum _{k=0}^p{n}_{k}h(y_k)$$

其中 $y_0:=x_t$ 和 $n_0$ 是 $x_t$ 重复的次数。 （从 ${\hat{\tau }}^1$ 到 ${\hat{\tau }}^2$ 的方差减小的证明很容易从双重积分论点得出。）我们再次坚持使用 $p$ 个并行处理单元计算 ${\hat{\tau }}^2$ 不会比使用单个处理单元计算 ${\hat{\tau }}^1$ 花费更多时间这一引人注目的特征。

为了保持其马尔可夫验证，算法必须在时间 $t+p$ 正确地继续。一个明显的选择是随机选择 $p$ 个序列中的一个，并将对应的 $x_{t+p}^{(j)}$ 作为下一个并行块的起点 $x_{t+p}$ 的值。这种机制如图 3 所示。虽然从马尔可夫的角度来看是有效的，但由于序列是由常规 IMH 算法产生的，这意味着从块 IMH 算法推导出的链与原始 IMH 算法的收敛速度完全相同.块起点的另一种选择是利用通过块结构计算的 $y_k$ 上的权重 $n_k$。实际上，这些权重本质上充当了重要性权重，它们允许选择任何 $p^2{x}_{t+i}^{(j)}$ 作为传入块的起点，这对应于选择建议的 $y_k$ 之一概率与 $n_k$ 成正比。虽然这个提议确实减少了结果链的长度，但它不会影响估计方面（仍然涉及所有 $p^2$ 值）并且它可以提高收敛性，因为加权 $y_k$ 的行为类似于样本的离散版本从目标密度 $\pi$。我们将不再介绍此替代方案。

块 IMH 算法的原始版本在算法 2 中描述。该算法由 $b$ 个块上的循环和每个块的 $p$ 个平行链上的内循环组成。这个内部循环的 $p$ 个步骤实际上是要并行实现的。算法 2 的输出是双倍的：

• 一条长度为 $T$ 的标准马尔可夫链，它由 $b$ 条长度为 $p$ 的链组成，每条链都是从算法 2 的第 12 行的 $p$ 条链中选择的，
• 一个 $p\times T$ 数组 $(x_t^k{)}_{t=1:T}^{k=1:p}$ ，估计器 ${\hat{\tau }}_2$ 是基于它的。
  

2.3 节省计算时间

由于目标密度 $\pi (y_k)$ 的逐点评估通常是算法中计算机密集程度最高的部分，因此采样额外的均匀变量在这里的影响可以忽略不计，与存储比原始 IMH 更大的向量相关的进一步成本也是如此.这尤其引人注目，因为多个链不需要比单个块执行时间存储得更远。这就是为什么，尽管我们在块 IMH 算法中采样了 $p$ 倍以上的 Uniform，但我们仍然认为它的运行成本与 IMH 算法大致相同。两者的目标密度评估数量确实相同，并且通常代表 Metropolis-Hastings 算法中计算时间的绝大部分。此外，均匀的的伪随机生成也可以受益于并行处理；例如，参见 L’Ecuyer 等人的工作。 （2001 年）。在接下来的蒙特卡罗实验中，将块 IMH 算法的不同版本以及标准 IMH 和重要性采样进行比较。我们强调，不与基于 $p$ 个独立并行链的普通并行算法进行比较的一个直接原因是它在成本方面没有多大意义。实际上，在目标密度评估方面，运行 $p$ 个相同长度 $T$ 的并行 MCMC 链确实要花费 $p$ 倍。显然，如果坚持运行 $p$ 个独立的链，例如从几个分散良好的起点初始化一个 MCMC 算法，那么这些链中的每一个都可以从我们的稳定方法中受益，这将改善最终的估计。此处针对尺寸为 $(p,p)$ 的方形块介绍了该方法，但块也可以是矩形的：该算法在使用 $r=p$ 排列时同样有效，导致 $(r,p)$ 块。我们在这里关注方块，因为当手边的机器提供 $p$ 个并行处理单元时，模拟建议值和均匀分布，并并行计算 $p$ 个建议值的目标密度和接受率是最有效的。再一次，具有 $p\times p$ 个正方形块的块 IMH 算法与原始 IMH 算法的成本大致相同，因为计算目标密度和接受率不仅仅补偿了在每个块末尾随机选择索引的成本。这等于说算法 1 的第 4 行和算法 2 的第 5 行是（迄今为止）各自算法中计算要求最高的。

2.4 玩具示例

我们现在介绍一个我们将在整篇文章中遵循的玩具示例。 目标 $\pi$ 是标准 $N(0,1)$ 正态分布的密度，提议 $\mu$ 是 $C(0,1)$ 柯西分布的密度。 因此，密度比为

我们只考虑积分 $\int x\pi (dx)$，在这种情况下 $\pi$ 的平均值为零。本例中 IMH 算法的接受率约为 70%。 （请注意，与其他 Metropolis-Hastings 算法相反，具有较高接受率的 IMH 被认为更有效；参见 Robert 和 Casella (2004)。）

在与随后提出的玩具示例相关的所有结果中，使用 10,000 次独立运行计算估计的方差。 $p$ 的值表示可用的并行处理单元的数量，范围从台式计算机的 4 个到集群或图形处理单元 (GPU) 的 100 个（对于配备多个 GPU 的计算机，这个值甚至可能更大）。 （这两个示例的代码都可作为补充材料获得。）

模拟实验的结果在下面的图 4-8 中以条形图的形式呈现，这表明与所比较的估计量相关的方差减少百分比，参考估计量是标准IMH 输出 ${\hat{\tau }}_1$。与块抽样的观点一致，相同的建议值和统一的绘制用于绘制在同一图表上的所有估计量（即对于给定的 $p$ 值），因此比较不会受到相关的附加噪声的干扰与模拟。

3. 排列

虽然算法 2 第 7 行中的排列选择与并行化验证无关，但它对方差改进具有重要影响，我们现在讨论几种自然选择。 Atchadé 和 Perron (2005) 的文章中出现了在 IMH 算法中测试建议值的各种阶数的想法，其中将排列选择为循环的。 我们首先列出排列的自然类型以及一些理由，然后我们凭经验比较它们对玩具示例的估计性能的影响。

3.1 五种自然排列

令 $S$ 为 { 1 , . . . , p } \{1, . . . , p\} {1,...,p}。 它的大小是 $p!$，因此太大而不能对所有排列进行平均，尽管这个解决方案是理想的。我们考虑在 $S$ 中找到 $p$ 个有效排列的更简单的选择，用 $({\sigma }_1,....,{\sigma }_p)$ 表示， 目标是有利于在第 2 节中定义的估计量 ${\hat{\tau }}_2$ 的方差最大可能减小的选择。

3.1.1 相同顺序

最基本的选择是在每个 $p$ 链上选择相同的排列：${\sigma }_1={\sigma }_2=\cdot \cdot \cdot ={\sigma }_p$。 这种选择听起来可能适得其反，但与 ${\hat{\tau }}_1$ 相比，使用这组排列我们仍然可以显着降低 ${\hat{\tau }}_2$ 的方差。 改进的原因是在 ${\hat{\tau }}_2$ 中使用的均匀性是 ${\hat{\tau }}_1$ 的 $p$ 倍，这导致了自然的 Rao-Blackwellization 现象，该现象将在第 4 节中详细研究。尽管如此，这组简单化的排列肯定不是最佳选择 因为它没有整合由建议值的任意排序产生的辅助随机性。

3.1.2 循环排列

3.1.3 随机排列

3.1.4 半随机半反转排列

3.1.5 分层随机排列

3.2 蒙特卡罗比较

我们比较了玩具示例中描述的五种排列类型。 图 4 显示了 $p$ 的各种值的结果，显示了与每个置换阶相关的 ${\hat{\tau }}_2$ 的方差减少，与原始 IMH 估计器 ${\hat{\tau }}_1$ 的方差相比。 对于 10,000 次独立复制中的每一次，块 IMH 算法在单个 $p\times p$ 块上启动，例如，使用第 2 节的符号 $b=1$，因为 $b$ 在此比较中没有任何作用。

如前所述，对 $p$ 条平行链中的每条链使用相同的 $y_k$ 顺序已经使估计量的方差显着降低了约 20%。 该模拟实验表明，三种随机排列（随机、半随机半反转和分层）在方差改进方面相当等效，并且它们明显优于循环排列建议，后者仅比“相同排列”方案。 因此，在接下来的蒙特卡洛实验中，我们将只使用随机阶解，最简单的随机方案。 考虑到对于具有并行能力的计算机而言，它基本上是免费获得的，当 $p\ge 32$ 时，像 35% 这样的改进量是相当可观的（Holmes et al. 2011）。

4.RAO-BLACKWELLIZATION

另一个可以超越经典 MH 算法的通用改进是 Rao–Blackwellization（Gelfand 和 Smith 1990；Casella 和 Robert 1996）。 在本节中，介绍了两种 Rao-Blackwellization 方法：一种计算免费，另一种计算量大。 然后我们在块 IMH 算法中实现这两种解决方案，并解释为什么“相同顺序”方案已经改进了 IMH 算法。

4.1 初级RAO-LACKWELLIZATION

在第 2.1 节的标准 IMH 算法中，可以通过简单的 Rao-Blackwellization 论证获得无成本的改进。 假设在步骤 $t+i$，$y_i$ 以概率 $\rho (x_{t+i-1},y_i)$ 被接受并以概率 $1-\rho (x_{t+i-1},y_i)$ 被拒绝，$y_i$ 的权重可以通过 $\rho (x_{t+i-1},y_i)$ 和 $x_{t+i-1}$ 对应的模拟值 $y_j$ 的权重可以类似地通过概率 $1-\rho (x_{t+i-1},y_i)$ 更新。 接下来考虑块 IMH 算法，在每个块的开头，我们可以定义 $p$ 个权重，用 $(w_k{)}_{k=1}^p$ 表示，初始化为 0，然后，对于 $p$ 个平行链中的第一个，用 $j$ 表示索引，使得 $x_{t+i-1}^{(1)}=y_j$ ，我们在每次 $t+i$ 更新这些权重为

5. 概率回归说明

6.结论

本文进行的蒙特卡罗实验表明，与标准 IMH 算法相比，所提出的方法显着提高了估计的精度。除了这些示例之外，我们还看到了多种情况，其中块 IMH 算法可用于改进具有挑战性问题的估计。首先，如前所述，IMH 算法可用于 Metropolis-within-Gibbs 算法（Gilks、Best 和 Tan 1995）。显然，如果对状态的每个组成部分执行单个 IMH 步骤，则无法应用块 IMH 技术而不产生额外成本。但是，多次而不是一次更新每个组件也是正确的。此外，均匀的 Gibbs 扫描很少是更新组件的最佳方式，并且已经研究了更复杂的方案，导致随机扫描 Gibbs 采样器和自适应 Gibbs 采样器（Latuszynski 和 Rosenthal 2010），其中更新给定组件的概率取决于组件，并且是沿着算法的迭代学习的。因此，如果一个组件的更新频率比另一个组件高 $n$ 倍，则可以连续执行 $n$ 个 IMH，这允许使用 $p=n$ 的块 IMH 技术。

IMH 步骤也用于顺序蒙特卡洛 (SMC) 采样器 (Chopin 2002;Del Moral 和 Doucet 2003)，以在重采样步骤后使粒子多样化。在这种情况下，可以通过在粒子上拟合（通常是高斯）分布来设计独立的提议。如果移动步骤连续重复多次，例如为了确保令人满意的粒子多样化，则可以使用块 IMH 算法。

与 SMC 相关，块 IMH 提供的方差减少可能有价值的另一个背景是粒子马尔可夫链蒙特卡罗方法类（Andrieu、Doucet 和 Holenstein 2010）。对于这些方法，在 MH 算法的每次迭代中都会计算一个粒子滤波器来估计目标密度，因此最重要的是充分利用这些估计所涉及的昂贵计算。因此，这是一个自然的并行化框架。

作为最后的消息，块 IMH 方法接近 100% 并行（除了在每个块的末尾随机抽取索引）。由于并行计算越来越容易使用，${\hat{\tau }}_3$带来的免费改进适用于IMH算法的所有实现。此外，即使不考虑并行计算，由于该方法使用了每个目标密度评估的大部分，因此在计算目标密度非常昂贵时带来了显着的改进。在这样的设置中，绘制 $p^2$ 而不是 $p$ 个统一变量的成本可以忽略不计，因此块 IMH 算法的运行时间与标准 IMH 算法大致相同。我们注意到，在算法中完成一个块所需的时间本质上是计算密度比 $w_i$ 所需的 $p$ 次的最大值。因此，如果这些时间变化很大，则随着标准 IMH 和块 IMH 算法的 $p$ 增加，计算时间的节省可能会减少。尽管如此，即使在这种极端情况下，在块 IMH 算法中使用 ${\hat{\tau }}_3$ 也会带来显着的方差改进，而基本上没有额外的成本。

补充材料 Python/R 代码：
补充文件包括一个 python 程序，可用于复制本文中提供的结果。 该程序名为 py-block-imh，也可在 Google Code 上找到。 调用 python 程序后，存档中包含的 R 程序用于创建图形。 请阅读存档中包含的 README.txt 文件以获取更多详细信息。 补充文件包含在单个存档中（并且可以通过单个下载获得）。 (py-block-imh.zip)