〖摘要〗马尔可夫链蒙特卡罗算法通过对分布的局部性探索来模拟复杂的统计分布。这种局部特征虽然不要求使用者了解目标分布性质，但也同时会导致对目标分布更长时间的探索，并且随着问题维度和数据复杂性的增加，对模拟样本数量的要求会也会增加。有几种技术可用于加速蒙特卡罗算法的收敛，无论是在探索层面（如回火、哈密顿蒙特卡罗和部分确定性方法）还是在开发层面（使用 Rao-Blackwellisation 和可扩展方法）。本文是对这些方法进行的一个综述。

〖原文〗 Robert, C.P. et al. (2018) ‘Accelerating MCMC algorithms’, Wiley Interdisciplinary Reviews: Computational Statistics, 10(5), p. e1435. Available at: https://doi.org/10.1002/wics.1435.

1 概述

马尔可夫链蒙特卡罗（MCMC）算法已经使用了近 60 年，在 1990 年代初成为分析贝叶斯复杂模型的参考方法（Gelfand 和 Smith，1990 [41]）。这种方法的优势在于：它能够保证收敛到感兴趣的量，并且对这些量背后的目标分布要求最低。从这个意义上说，MCMC 算法是稳健或通用的，这与需要从目标分布直接模拟最标准的蒙特卡罗方法（例如，参见 Rubinstein，1981 [86]；Robert 和 Casella，2004 [79]）相反。

然而，这种稳健性可能会导致缓慢的收敛行为，因为对相关空间的探索可能需要很长时间，因为模拟通常在当前位置附近跳跃推进。换句话说，MCMC （尤其是其现成版本，如 Gibbs 采样和 Metropolis-Hastings 算法）通常是短视的，因为它提供了对局部区域的良好探索，但同时并不知道分布的全局支持。

与大多数其他模拟方法一样，总是存在通过利用目标分布的结构信息来创建高收敛 MCMC 算法的一些方法。在此，我们把范围主要限制在目标分布来自于计算机代码输出的现实情况，或者信息内容同样有限的类似情况。

加速 MCMC 算法的方法可以分为几类，从 “提高对目标分布的了解” 到 “修改算法中的提议方案”，也包括更好地利用原始 MCMC 算法结果。以下部分提供了有关这些方向的更多详细信息以及文献中提出的解决方案。

传统 MCMC 的局限性

MCMC 方法的历史（参见，例如 Cappé 和 Robert，2000 年 [22]）与蒙特卡洛方法大致同时开始，并与第一台计算机的概念相结合。它们被设计用于处理复杂目标分布的模拟，当复杂性源于目标密度的形状、相关数据的大小、要模拟的对象的尺寸或时间要求时。例如，目标密度 π(θ)\pi(\theta) 可能碰巧用无法解析求解的多重积分表示：

π(θ)=∫ω(θ,ξ)dξ\pi(\theta)=\int \omega(\theta, \xi) \mathrm{d} \xi

这需要模拟整个向量 (θ,ξ)(\theta, \xi)。在 ξ\xi 与数据具有相同维度的情况下，例如在隐变量模型中，要模拟的对象维度的这种显著增加会给标准蒙特卡洛方法带来计算困难，无法管理新目标 ω(θ,ξ)\omega(\theta, \xi)，以设计一种新的高效仿真算法。马尔可夫链蒙特卡罗 (MCMC) 算法允许通过模拟探索感兴趣空间（以及可能的辅助变量的补充空间）的马尔可夫链来解决这一计算挑战，而无需深入了解密度 π\pi，除了计算给定参数值 θ0\theta_0（如果达到归一化常数）和可能的梯度 ∇log⁡π(θ0)\nabla \log \pi\left( \theta_0\right)。该方法的验证（例如，Robert 和 Casella，2004 [79]）是马尔可夫链是遍历的（例如，Meyn 和 Tweedie，1993 [63]），即它在分布上收敛于密度为 π\pi 的分布，无论马尔可夫链在时间 t=0t=0 从哪里开始。

Metropolis-Hastings 算法是该原理的一般说明。基本算法是通过选择一个提议来构建的，即条件密度 K(θ′∣θ)K\left(\theta^{\prime} \mid \theta\right)（也称为马尔可夫核），马尔可夫链 {θt}t=1∞\left \{ \theta_t \right \}_{t=1}^{\infty} 是通过连续模拟转移推导出来的

θt+1={θ′∼K(θ′∣θt) with probability {π(θ′)π(θt)×K(θt∣θ′)K(θ′∣θt)}∧1,θt otherwise. \theta_{t+1}= \begin{cases}\theta^{\prime} \sim K\left(\theta^{\prime} \mid \theta_t\right) & \text { with probability }\left\{\frac{\pi\left(\theta^{\prime}\right)}{\pi\left(\theta_t\right)} \times \frac{K\left(\theta_t \mid \theta^{\prime}\right)}{K\left(\theta^{\prime} \mid \theta_t\right)}\right\} \wedge 1, \\ \theta_t & \text { otherwise. }\end{cases}

如果得到的马尔可夫链 {θt}t=1∞\left\{\theta_t\right\}_{t=1}^{\infty} 是不可约的（即，在有限次迭代中访问 π\pi 支持的任何区域的概率为正），则该算法的这种接受-拒绝特性使其适合将 π\pi 作为其平稳分布。（平稳性可以很容易地显示出来，例如，通过使用使链时间可逆的所谓的细致平衡属性，参见 Robert 和 Casella，2004 [79]）。

考虑到从目标分布 π\pi 模拟样本的初始目标，像上面的 Metropolis-Hastings 算法这样的 MCMC 方法的性能通常会有很大差异，主要取决于提议 KK 和目标 π\pi 之间的对应性。例如，如果 K(θ∣θt)=π(θ)K(\theta|\theta_t) = \pi(\theta)，Metropolis-Hastings 算法将简化为 i.i.d。从目标采样，这当然是 i.i.d. 时的正式选项。从 π\pi 中采样被证明是不可能实现的。尽管存在马尔可夫链 θtt=1∞{\theta_t}^\infty_{t=1} 导致链的连续项之间存在负相关的罕见情况，使其比常规独立同分布的采样（Liu 等, 1995 [55]）更有效。，最常见的情况是模拟值之间的正相关之一（有时是一致的，参见 Liu 等, 1994 [54]）。此功能意味着算法效率降低，因此需要大量模拟才能达到与基于独立同分布的近似值相同的精度。模拟（不考虑计算时间的差异）。更一般地说，MCMC 算法可能需要大量迭代来摆脱其起点 θ0\theta_0 的吸引力并达到平稳性，以至于此类算法的某些版本无法在可用时间内收敛（即，在实践中，如果不是理论上）

因此，寻求以下加速是非常有意义的：

(1) 加速指定 MCMC 算法收敛到其平稳分布

(2) 加速指定 MCMC 算法收敛以计算其期望

(3) 加速对指定目标分布的 MCMC 探索

这些目标是相关的，但仍然有不同。例如，通过从目标分布模拟初始化的链可能仍然无法在可接受的迭代次数内探索整个支持区域。虽然这个问题没有最佳和通用的解决方案，但我们将在下面讨论尽可能通用的方法，而不是利用特定目标分布的数学结构的人工方法。

理想情况下，我们仅覆盖目标分布是计算机代码输出的现实性情况已知的现实情况。务实地，我们在这里还介绍了在应用于足够广泛的问题时需要更多努力和校准步骤的解决方案。

2 利用目标分布的几何形态

虽然尝试构建更有效和更快的 MCMC 算法没有尽头，而且这个（无尽的）目标需要考虑在有限的资源预算下设计此类替代方案的成本，但存在几种通用解决方案，使得给定的目标在构造算法之前，可以首先根据密度的几何（或拓扑）进行探索。尽管这种类型的方法以某种方式使我们偏离了改进现有算法的最初目的，但它们在本次调查中仍然有意义，因为它们允许几乎自动化的实现。

2.1 哈密顿蒙特卡洛方法

从这篇综述的角度来看，哈密顿（或混合）蒙特卡罗（HMC）是一种辅助变量技术，它利用连续时间马尔可夫过程从目标 ππ 中采样。这种方法来自物理学（Duane、Kennedy、Pendleton 和 Roweth，1987 年 [37]），并由 Neal（1999 年 [70]、2011 年[71]）和 MacKay（2002 年[58]）在统计学中推广。给定目标 π(θ)π(θ)，其中 θ∈Rdθ \in \mathbb{R}^d，人工辅助变量 ϑ∈Rdϑ \in \mathbb{R}^d 与密度 ϑ(ϑ∣θ)ϑ(ϑ|θ) 一起被引入，使得 (θ,ϑ)(θ, ϑ) 的联合分布以 π(θ)π(θ) 作为其边缘。虽然此表示具有完全的自由度，但 HMC 文献通常将 ϑϑ 称为位于 θθ 处的粒子的动量，以类比物理学。基于联合分布的表示：

ω(θ,θ)=π(θ)ϖ(θ∣θ)∝exp⁡{−H(θ,θ)}ω(θ,θ) = π(θ)ϖ(θ|θ) \propto \exp\{−H(θ,θ)\}

其中 H(⋅)H(·) 称为哈密顿量，哈密顿量蒙特卡洛 (HMC) 与所谓的哈密顿量方程生成的连续时间过程 (θt,ϑt)(θ_t, ϑ_t) 相关联。

dθtdt=∂H∂θ(θt,θt)dθtdt=−∂H∂θ(θt,θt)d\theta_t dt=\partial H\partial θ(θ_t,θ_t) d\theta_t dt=−\partial H\partial θ(θ_t,θ_t)

随着时间的推移保持哈密顿目标稳定，因为

dH(θt,θt)dt=∂H∂θ(θt,θt)dθtdt+∂H∂θ(θt,θt)dθtdt=0dH(\theta_t ,\theta_t ) dt=\partial H\partial θ(\theta_t ,\theta_t )d\theta_t dt+\partial H\partial θ(\theta_t ,\theta_t )d\theta_t dt = 0

显然，上述连续时间马尔可夫过程是确定性的，只探索给定的水平集，

{(θ,θ):H(θ,θ)=H(θ0,θ0)}\{(θ,θ):H(θ,θ) = H(\theta_0,\theta_0)\}

而不是整个增强状态空间 R2d\mathbb{R}2d，这会引发不可约性问题。这个问题的一个可接受的解决方案是在随机时间 {τn}n=1∞\{τ_n\}^{\infty}_{n=1} 刷新动量 ϑt∼ϖ(ϑ∣θt−)ϑ t \sim ϖ(ϑ|θ_{t^{-}})，其中 θt−θ_{t^{-}} 表示紧接在时间 tt 之前的 θθ 的位置，并且随机持续时间 {τn−τn−1}n=2∞\{τ_n − τ_{n−1}\}^{\infty}_{n=2} 服从指数分布。通过构造，连续时间哈密顿马尔可夫链可以被视为使用哈密顿动力学的特定分段确定性马尔可夫过程（Davis，1984 年 [28]，1993 年[29]；Bou-Rabee 等，2017 年 [17]），我们的目标 ππ 是其相关的边缘不变分布。

在进入概念的实际实施之前，让我们指出机器中的自由齿轮是条件密度 ϖ(θ∣θ)ϖ(θ|θ)，它通常被选为高斯密度，具有对应于目标协方差或局部曲率取决于黎曼 HMC 中的 θθ (Girolami & Calderhead, 2011 [44])。 Betancourt (2017 [12]) 支持这两个案例反对非高斯替代方案，Livingstone、Faulkner 和 Roberts (2017 [57]) 分析了 HMC 中不同的动能选择如何影响算法性能。对于固定的协方差矩阵，哈密顿方程变为：

dθtdt=M−1θtdθtdt=∇L(θt)d\theta_t dt=M−1\theta_t d\theta_t dt=\nabla L(\theta_t )

这是得分函数。因此，过程的速度（或动量）由这个得分函数驱动，即对数目标的梯度。

上面的描述仍然很概念化，因为没有通用的方法来产生这个连续时间过程，因为哈密顿方程在大多数情况下无法精确求解。此外，像欧拉方法这样的标准数值求解器会创建一个不稳定的近似值，当过程偏离其真实轨迹时会产生偏差。然而，存在一种离散化模拟技术，它可以产生马尔可夫链并且非常适合哈密顿方程，因为它保留了平稳分布（Betancourt，2017 [12]）。它被称为辛积分器，在具有恒定协方差的独立情况下的一个版本包括以下（所谓的蛙跳）步骤

ϑt + ε/2\theta_t +εεt+ε=εt+ε\nablaL(\theta_t )/2,=\theta_t +εM−1εt+ε/2,=εt+ε/2+εεL(\theta_t +ε)/2

其中 εε 是时间离散化步骤。使用从高斯辅助目标得出的 ϑ0ϑ0 的建议，并通过 Metropolis–Hastings 步骤决定接受 (θ T arepsilon , ϑ T arepsilon ) 的值可以限制错过目标的危险。请注意，前两个 leapfrog 步骤会导致 θtθ t 上的 Langevin 移动：

\theta_t + arepsilon = \theta_t + arepsilon 2M−1\nablaℒ(\theta_t )/2 + arepsilon M−1ϑt

因此与下文讨论的 Metropolis-adjusted Langevin 算法 (MALA) 相关联（有关 εε 最佳选择的理论讨论，请参阅 Durmus 和 Moulines，2017 [38]）。请注意，leapfrog 积分器是准确度（因为它是二阶准确度）和计算效率之间非常有吸引力的中间地带。

在实践中，重要的是要注意哈密顿动力学的离散化引入了两个自由参数，步长 $ arepsilon $ 和轨迹长度 $T arepsilon $，两者都需要校准。作为 HMC 的一个经验成功且流行的变体，Hoffman 和 Gelman（2014 [47]）的“不掉头采样器”（NUTS）根据原始对偶平均调整了 εε 的值。它还消除了通过递归算法选择轨迹长度 TT 的需要，该算法为许多向前和向后的蛙跳步骤构建一组候选建议，并在模拟路径后退时自动停止。

Rasmussen (2003 [76]) 提出了该领域的进一步加速步骤（另见 Fielding, Nott, & Liong, 2011 [40]），即用近似值 π(⋅)π^(⋅) 替换精确目标密度 π(⋅)π(·) 在 HMC 算法的多次迭代中计算速度更快。构建这种近似的一种通用方法是依赖高斯过程，当解释为目标密度 π(⋅)π(·) 的先验分布时，它仅在 θ，π(θ1)，…、π(θn)θ，π(θ 1)，\ldots 、π(θ n) 的某些值处观察到）（拉斯穆森和威廉姆斯，2005 年）。该解决方案正在加速算法，可能提高了几个数量级，但它引入了对蒙特卡洛方法的进一步近似，即使在 Fielding 等 (2011) 的 leapfrog 离散化结束时使用真实目标时也是如此。

Stan（以 Stanislas Ullam 命名，参见 Carpenter 等，2017 年 [23]）是一种用于贝叶斯推理的计算机语言，除其他近似技术外，它还实现了 NUTS 算法以消除手动调整。更准确地说，Stan 是一种概率编程语言，因为输入处于统计模型级别，连同数据，而不是 MCMC 算法的细节。算法部分在某种程度上是自动化的，这意味着当可以通过这种语言方便地定义模型时，它提供了一种替代产生原始链的采样器的方法。作为 HMC 带来的加速的说明，图 Figure1,1 转载自 Hoffman 和 Gelman (2014 [47])，显示了 NUTS 与随机游走 MH 和 Gibbs 采样器相比的性能。

3 通过问题分解加速 MCMC

近年来，“大” 数据集收集和分析的爆炸式增长给用于贝叶斯推理的 MCMC 算法带来了新的挑战。当在接受-拒绝步骤检查一个新的提议样本是否被接受时，MCMC 算法（例如 Metropolis-Hastings 版本）需要在每次迭代时扫描整个数据集，以评估似然函数。 MCMC 算法很难扩大规模，这极大地阻碍了它们在大数据环境中的应用。在某些情况下，数据集可能太大而无法容纳在单台机器上。也可能是保密措施将不同的数据库强加在不同的网络上，这可能会增加加密数据的负担（Aslett 等，2015 年 [5]）。在涉及数千或数十万次迭代的 MCMC 规模上，不同机器之间的通信可能被证明是不可能的。

3.1 可扩展的 MCMC 方法

近年来，人们努力设计可扩展的算法，即通过将问题分解为可管理或可扩展的部分来设法处理大规模目标的解决方案。粗略地说，这些方法可以分为两类（Bardenet、Doucet 和 Holmes，2015）：分而治之的方法和子采样方法。

分而治之的方法将整个数据集（表示为 ）分成批次 {1,…,k}\{1,\ldots ,k\}，并对每个数据批次独立运行单独的 MCMC 算法，就好像它们是独立的贝叶斯推理问题一样。1 这些方法然后将模拟参数结果组合在一起以近似原始后验分布。根据在 MCMC 阶段选择的批次的处理，这些方法可以进一步细分为两个更精细的组：后验方法和增强后验方法。次后验方法受独立乘积方程的启发：

π(θ)∝∣i=1k(π0(θ)1/k∣ℓ∈Xip(xℓ∣∣∣∣θ))=∣i=1kπi(θ)(1)π(θ) \propto |i=1k(π0(θ)1/k|ℓ\in Xip(xℓ||||θ))=|i=1kπi(θ) (1)

并且他们在各自的 MCMC 步骤中以密度 πi(θ)π i(θ)（最大为常数）为目标。因此，他们通过在每个批次上独立运行 MCMC 采样器来绕过通信成本（Scott 等，2016 年[88]），并且他们通常会增加 MCMC 混合率（在第二个产生的有效样本量中），因为后验分布 πi(θ)π i(θ ) 基于较小的数据集。例如，Scott 等 (2016 [88]) 通过高斯重新加权组合来自子后验的样本 πi(θ)π i(θ)。 Neiswanger、Wang 和 Xing (2013 [74]) 通过非参数和半参数方法估计子后验 πi(θ)π i(θ)，并且他们在这些估计量的乘积上运行额外的 MCMC 采样器以逼近真实的后验 π(θ)π(θ)。 Wang 和 Dunson (2013 [97]) 使用额外的 Weierstrass 采样器改进了这个乘积估计器，而 Wang、Guo、Heller 和 Dunson (2015 [98]) 通过使用阶跃函数划分样本空间来估计后验。

作为从 subposteriors 采样的替代方法，boosted subposteriors 方法的目标是组件

π∼i(θ)∝π0(θ)(∠ℓ∈Xip(xℓ∣θ))k(2)π \sim i(θ) \propto π0(θ)(∠ ℓ \in Xip(xℓ|θ))k (2)

在单独的 MCMC 运行中。由于它们在形式上相当于将每批重复 kk 次以生成与真实数据集大小相同的伪数据集，因此得到的增强子后验 π∼1(θ),…,π∼k(θ)π \sim 1(θ),\ldots ,π\sim k(θ) 在参数 θθ 的每个分量的方差作为真实后验，因此可以被视为一组真实后验的估计量。在随后的组合阶段，这些子后验被合并在一起，以构建一个更好的目标分布近似值。例如，Minsker、Srivastava、Lin 和 Dunson（2014 年[65]）使用增强子后验的几何中值来近似后验，将它们嵌入到相关的再生核 Hilbert 空间中，而 Srivastava、Cevher、Dinh 和 Dunson（2015 年 [89]）实现了这一目标使用 π∼1,…π \sim 1,\ldots ,π∼kπ \sim k 的质心，这些质心是根据 Wasserstein 距离计算的。

从不同于上述分而治之并行方案的角度来看，子采样方法旨在减少每次迭代时操作的单个数据点似然评估的数量，以加速 MCMC 算法。从一般的角度来看，这些方法可以进一步分为两个更细的类别：精确子采样方法和近似子采样方法，具体取决于它们的结果输出。精确子采样方法通常需要在每次迭代时随机大小的数据子集。解决此问题的一种方法是利用伪边缘 MCMC，通过构建对数据子集评估的目标密度的无偏估计量 (Andrieu & Roberts, 2009 [3])。 Quiroz、Villani 和 Kohn（2016 年 [75]）通过结合 Rhee 和 Glynn（2015 年 [78]）强大的去偏技术和 Deligiannidis、Doucet 和 Pitt（2015 年 [31]）的相关伪边缘 MCMC 方法遵循这个方向。另一个方向是使用分段确定性马尔可夫过程 (PDMP) (Davis, 1984 [28], 1993 [29])，该过程将目标分布作为其不变分布的边缘。这个 PDMP 版本需要对数似然函数梯度的无偏估计，而不是似然本身。通过对相关泊松过程的事件率函数使用足够严格的界限，PDMP 可以产生超高效的可扩展 MCMC 算法。弹力粒子采样器 (Bouchard-Côté, Vollmer, & Doucet, 2017 [18]) 和之字形采样器 (Bierkens, Fearnhead, & Roberts, 2016 [14]) 是两个相互竞争的 PDMP 算法，而 Bierkens 等 (2017 [15]) 统一和扩展这两种方法。此外，应该注意到 PDMP 产生不可逆马尔可夫链，这意味着与可逆 MCMC 算法（如 MH、HMC 和 MALA）相比，该算法在混合率和渐近方差方面应该更有效，因为在一些理论和实验著作中观察到 (Bierkens, 2016 [14]; Chen & Hwang, 2013 [27]; Hwang, Hwang-Ma, & Sheu, 1993 [48]; Sun, Gomez, & Schmidhuber, 2010 [91])。

近似子采样方法旨在构建目标分布的近似值。除了前面提到的 Rasmussen (2003 [76]) 和 Fielding 等 (2011 [40]) 的尝试。一个方向是通过使用数据子集来高精度地近似接受概率 (Bardenet 等, 2015 [10]; Bardenet, Doucet, & Holmes, 2014 [9])。另一种解决方案是基于对精确方法的直接修改。 Welling 和 Teh（2011 年 [99]）的开创性工作，随机梯度朗之万动力学 (SGLD)，是利用朗之万扩散

dθt=12Λ∇log⁡π(θt)dt+Λ1/2dBt,θ0∈Rd,t[0,∞)(3)d\theta_t = \frac{1}{2} Λ \nabla \log π(\theta_t )dt + Λ1/2dBt, \theta_0 \in \mathbb{R}^d, t [0, \infty )(3)

其中 ΛΛ 是用户指定的矩阵，ππ 是目标分布，4B t4 是 dd 维布朗过程。凭借 Euler-Maruyama 离散化和使用对数目标密度梯度的无偏估计，SGLD 及其变体 (Chen, Fox, & Guestrin, 2014 [26]; Ding 等, 2014 [32]) 通常在以下方面产生快速准确的结果与使用 MH 步骤的 MCMC 算法进行比较时的实践。

图 2 显示了共识蒙特卡罗算法（Scott 等，2016 [88]）与使用整个数据集的 Metropolis-Hastings 算法相比的时间要求，而图 3 显示了 Bardenet 等(2015 [10])的置信度采样器中似然评估的节省。

3.2 并行化和分布式方案

现代计算体系结构由多个计算单元构建而成，这些计算单元允许完全独立或通过某些通信进行并行处理。尽管 MCMC 的马尔可夫性质本质上是顺序的并且与并行化的概念有些不同，但是在文献中已经提出了几个部分解决方案来利用这些并行架构。最简单的方法是并行运行多个 MCMC 链，对所有其他链都视而不见，直到分配的计算时间耗尽。最后，对所有链的估计量进行平均。然而，这种天真的实现可能会受到以下事实的影响，即其中一些链在计算时间结束时尚未达到其静止状态，这会导致结果估计出现偏差。尽管在文献中可以找到几种方法（Guihenneuc-Jouyaux & Robert, 1998 [45]；Jacob, O’Leary, & Atchadé, 2017 [50]；Mykland 和 Tierney 等，1995 [69]）。在相反的极端，复杂的目标可能表示为涉及许多必须评估的项的乘积，每个项都可以在相乘之前归因于不同的线程。该策略需要在每个 MCMC 步骤中处理器之间进行通信。中间版本 (Jacob, Robert, & Smith, 2011 [51]) 包括并行运行多个马尔可夫链，并定期选择参考链，所有模拟都通过 Rao-Blackwell 方案循环使用（另请参阅 Calderhead，2014 年的类似方案）。 Martino、Elvira、Luengo、Corander 和 Louzada (2016 [61]) 提出了交互正交 MCMC 方法族 (O-MCMC)，旨在促进对状态的更好探索空间，特别是在高维和多峰值目标中。多个 MCMC 链并行运行，使用随机游走提议探索空间。平行链定期共享信息，也通过联合 MCMC 步骤，从而允许全局（协调）探索和局部近似的有效组合。 O-MCMC 方法还允许并行实施 Multiple Try Metropolis。在 Calderhead (2014 [19]) 中，Metropolis-Hastings 算法的推广允许直接并行化。在每次 MCMC 迭代中，每个建议的点都可以在不同的处理器中进行评估。最后，请注意关于可扩展 MCMC 的部分还包含可并行化的方法，例如 Angelino、Kohler、Waterland、Seltzer 和 Adams (2014 [4]) 的预取方法（有关相关方法，另请参见 Banterle、Grazian、Lee 和 Robert，2015 年 [8]） ，主要基于目标的近似值）。最近一项称为异步 MCMC（Terenin、Simpson 和 Draper，2015 年 [92]）的努力旨在通过减少并行线程之间的交换量来提高并行化收益，但这一概念在现阶段仍然保密。

4 通过改进来加速 MCMC

本着与上一节相同的精神，本节通过考虑 MCMC 算法本身的可能修改来扩展本文的目的，而不仅仅是利用给定 MCMC 算法的输出。例如，设计一个 HMC 算法是这个问题的答案，即使 “改进” 并没有得到保证。尽管如此，我们在这里的论点是，一旦提供了这个输出，就有可能以半自主的方式推导出新的提案。

4.1 模拟退火

dd 维状态空间 ΘΘ 上的目标分布 π(θ)π(θ) 可以表现出多峰值，概率质量位于状态空间的不同区域。大多数 MCMC 算法使用局部提议机制，该机制针对局部近似最优性进行调整，例如，参见 Roberts、Gelman 和 Gilks (1997 [81]) 以及 Roberts 和 Rosenthal (2001 [82])。通过构造，这些局部提议导致马尔可夫链“陷入”状态空间的一个子集中，这意味着在有限的运行时间内，链可能完全无法探索状态空间中的其他模式，从而导致样本偏差。加速 MCMC 的策略通常使用局部梯度信息，这会将链拉回到模式的中心，这与多模式设置中所需的相反。

有一系列方法可用于克服 MCMC 中的多峰值问题，其中大部分使用状态空间增强。允许马尔可夫链探索整个状态空间的辅助分布是有针对性的，然后传递它们的混合信息以帮助在真正的目标中混合。虽然上一节的子后验可以看作是下一节的特例，但这些方法最成功和最方便的实现是使用功率调节目标分布。对于 β∈(0,1]β \in (0, 1]，在逆温度水平 ββ 下的目标分布定义为

πβ(θ)=(β)[π(θ)]β其中(β)=[∫[π(θ)]βdθ]−1πβ(θ) = (β)[π(θ)]β 其中 (β) = [∫[π(θ)]βdθ]−1

因此，π1(θ)=π(θ)π 1(θ) = π(θ)。温度 β<1β < 1 使目标分布变平，允许链探索整个状态空间，前提是 ββ 值足够小。模拟回火 (ST) 和平行回火 (PT) 算法 (Geyer, 1991 [43]; Marinari & Parisi, 1992 [59]) 通常使用功率回火目标来克服多峰值问题。 ST 方法在增强状态空间 {B,Θ}\{B, Θ\} 上运行单个马尔可夫链，其中 B={β0,β1,…,βn}B = \{β 0, β 1, \ldots , β n\} 是 nn 个逆温度水平的离散集合，其中 1=β0>β1>…⟩>βn>01 = β 0 > β 1 > \ldots ⟩> β n > 0。该算法通过在空间的 ΘΘ 和 BB 分量的更新之间循环使用 Metropolis-within-Gibbs 策略。例如，提议的温度交换移动 βi⟩→⟩βjβi⟩→⟩βj 被概率接受

min⁡{1,πβj(θ)πβi(θ)}\min \{1,πβj(θ)πβi(θ)\}

为了保持详细的平衡。请注意，此接受率取决于通常未知的归一化常数 (ββ)，尽管有时可以估计它们，例如 Wang 和 Landau (2001 [96]) 以及 Atchadé 和 Liu (2004 [6])。如果边缘归一化常数的估计不切实际，则采用 PT 算法。该方法同时在 n+1n + 1 个温度水平中的每个温度水平运行马尔可夫链，目标是 ∣i=0n[π(θi)]βi|i=0n[π(θi)]βi 给出的联合分布。根据不再取决于边缘归一化常数的比率，接受相邻温度水平下链之间的交换移动。事实上，这种权力调节方法已成功应用于许多环境并被广泛使用，例如，Neal (1996 [72])、Earl 和 Deem (2005 [39])、Xie，Zhou 和 Jiang (2010 [102])、Mohamed、Calderhead、Filippone 等 (2012 [68]) 以及 Carter 和 White (2013 [24])。

在这两种方法中，都有一个“金发姑娘”原则来设置逆向温度计划。 “太大”的温度水平之间的间距会导致很少被接受的交换移动，从而延迟热态混合信息向冷态的传输。另一方面，太小的间距需要大量的中间温度水平，再次导致通过温度空间的缓慢混合。随着 ΘΘ 维数的增加，这个问题变得更加困难。

许多历史文献表明几何间距是最优的，即存在 c∈(0,1)c \in (0, 1) 使得 βi+1=cβii=0,1,…,nβ i + 1 = cβ i i = 0, 1, \ldots , n。然而，对于 ST 版本，Atchadé、Roberts 和 Rosenthal（2011 年 [7]）通过最大化 BB 空间中温度交换移动的（维度渐近）预期平方跳跃距离，将该问题视为最佳缩放问题。在限制性假设下，他们表明连续的逆温度水平之间的间距应与 O(d−1/2)O(d−1/2) 的维度成比例，以防止交换移动接受率的退化。根据 Gelman、Gilks 和 Roberts（1996 年[42]）的说法，对于从业者而言，该结果为最佳设置提供了指导，因为它建议在连续的逆温度水平之间相应的最佳交换移动接受率为 0.2340.234。最后，与历史上推荐的几何时间表相反，作者建议温度时间表设置应该连续构建，以便在连续级别之间产生大约 0.2340.234 的交换接受率；这是在 Miasojedow、Moulines 和 Vihola (2013 [64]) 中自适应实现的。 Roberts 和 Rosenthal（2014 年 [85]）证明使用预期平方跳跃距离作为混合速度的度量是合理的，在与 Atchadé 等（2011 [7]）相同的条件下，表明 ST 链的温度分量具有相关的扩散过程。

0.2340.234 接受率的目标为在某些设置中设置 ST/PT 算法提供了很好的指导，但对于遵循此规则以获得最佳设置的从业者来说，这是一个重要的警告。 Atchadé 等 (2011) 所做的假设以及 Roberts 和 Rosenthal (2014) 忽略了温度水平内混合的限制，而是假设相对于温度空间内的混合，这可以无限快地完成。 Woodard、Schmidler 和 Huber (2009a , 2009b[100][101]) 以及 Bhatnagalr 和 Randall (2016 [13]) 对 ST/PT 链的光谱差距进行了全面分析，他们的结论相当否定使用功率回火目标的 ST/PT 方法.本质上，在模式具有不同结构的情况下，ST/PT 算法达到给定收敛水平所需的时间可以在维度上呈指数增长。一个主要原因是基于功率的回火不保持不同温度水平区域之间的相对重量/质量，见图 4.4。这个问题可以在维度上呈指数级扩展。从实际的角度来看，在这些有限运行的高维非相同峰值结构设置中，交换接受率可能非常具有误导性，这意味着它们在诊断联运混合质量方面的用途有限。

4.2 自适应 MCMC

改进和校准 MCMC 算法以更好地对应预期目标是使算法更高效的自然步骤，前提是有足够的有关此目标分布的信息可用。例如，当与该目标关联的 MCMC 样本可用时，即使它没有完全探索目标的范围，它也包含一些信息，然后可以利用这些信息来构建新的 MCMC 算法。文献中可用的一些解决方案（例如，Liang、Liu 和 Carroll，2007 年 [53]）通过重复 MCMC 迭代块并在每个块后更新建议 KK 来进行，以达到特定的最优目标，如特定的接受率，如 0.2340.234 Metropolis-Hastings 台阶（Gelman 等，1996 年 [42]）。此方法的大多数版本都会根据先前的实现 (Robert & Casella, 2009 [80]) 或整个样本 (Douc, Guillin, Marin, & Robert, 2007 [33]) 更新随机游走建议的尺度结构，从而将该方法转变为迭代具有马尔可夫依赖性的重要性采样。 （也可以将其视为粒子过滤的静态版本，Doucet、Godsill 和 Andrieu，2000 年 [36]；Andrieu 和 Doucet，2002 年 [1]；Storvik，2002 年 [90]。）

其他自适应解决方案绕过这种初步的和有点特别的构造，而是旨在算法内的永久更新，其动机是持续适应不断提高与目标的对应性。为了保持方法的有效性 (Gelman 等, 1996 [42]; Haario, Saksman, & Tamminen, 1999 [46]; Roberts & Rosenthal, 2007 [83]; Saksman & Vihola, 2010 [87])，即算法产生的链收敛到预期目标，需要建立特定的收敛结果，因为标准 MCMC 算法背后的遍历定理不适用。如果不小心（见图 5），自适应 MCMC 算法可能会由于过度拟合而无法收敛。自适应性的一个缺点是建议分布的更新过于依赖早期的模拟，从而加强了对尚未探索的空间部分的排除。

因此，为了验证自适应 MCMC 方法，必须对算法施加更严格的约束。一个描述良好的解决方案 (Roberts & Rosenthal, 2009 [84]) 被称为递减适应。非正式地，它包括在两个连续的提案内核之间施加一个距离以均匀地减小到零。实际上，这意味着像 Haario 等 (1999 [46]) 的早期提议一样，通过类似山脊的因素来稳定提议的变化。该决议的一个缺点是，减少本身必须进行校准，并且很可能无法对原始提案带来重大改进。

4.3 多次尝试 MCMC

改进 MCMC 算法中使用的原始建议的一种完全不同的方法是考虑建立在不同原理和实验基础上的建议集合。多次尝试 MCMC 算法（Liu, Liang, & Wong, 2000 [56]; Bédard, Douc, & Moulines, 2012 [11]; Martino, 2018 [60]）遵循这一观点。顾名思义，多次尝试 MCMC 算法的出发点是同时提出马尔可夫链的 NN 个潜在移动 θ1t,…,θNtθ1t,\ldots ,θNt，而不是单个值。提议值θit可以根据以马尔可夫链的当前值 θt\theta_t 为条件的N个不同的提议密度 Ki(⋅∣θt)Ki(·|\theta_t ) 独立生成。其中一个 θitθit 是根据重要性采样权重 wit∝π(θit)/Ki(⋅∣θt)wit\propto π(θit)/Ki(⋅|\theta_t ) 选择的。然后，所选值被进一步的 Metropolis-Hastings 步骤接受，该步骤涉及重要性阶段的标准化常数比率，一个对应于先前所做的选择，另一个对应于为此目的创建的选择。事实上，除了计算重要性权重之和和生成不同变量的额外成本外，这种方法还面临着不可忽略的缺点，即需要 N−1N - 1 次补充模拟，这些模拟仅用于实现详细平衡和计算重要性的向后求和权重。当考虑一组独立的 Metropolis-Hastings 提议 q(θ)q(θ) 时，这种约束可能会消失，但这种设置很少是现实的，因为它需要一定量的先验知识或实验来构建相关分布。

文献中发现的另一种方法是集成蒙特卡洛（Iba，2000 年[49]；Cappé、Douc、Guillin、Marin 和 Robert，2008 年 [20]；Neal，2011 年 [73]；Martino，2018 年 [60]），如图 6 所示，它在每次迭代时生成一个完整的样本，以初始目标的产品为目标，更接近于粒子方法（Cappé、Guillin、Marin 和 Robert，2004 年 [21]；Mengersen 和 Robert，2003 年 [62]）。

这一原则的另一种实现称为延迟拒绝（Tierney & Mira, 1998 [93]; Mira, 2001 [66]; Mira & Sargent, 2003 [67]），一旦先前提议的值被拒绝，就会按顺序考虑提议，通过考虑加速 MCMC 几种可能性，如果顺序的话。这种方法的计算困难在于，随着延迟次数的增加，相关的接受概率会变得越来越复杂，这可能会削弱它相对于同时进行多次尝试的吸引力。另一个困难是以足够多样化的方式设计提案的顺序。

5 方差减少方法

由于 MCMC 的主要目标是生成形式感兴趣的数量的近似值

Ih=∫Θh(θ)π(θ)dθ\mathfrak{I}_h=\int_{\Theta} h(\theta) \pi(\theta) \mathrm{d} \theta

加速这些算法的另一种（和累积的）方法是提高从 MCMC 输出得出的近似值的质量。也就是说，给定一个 MCMC 序列 θ1,…,θT\theta_1, \ldots, \theta_T，收敛到 π(⋅)\pi(\cdot)，我们可以超越基本的蒙特卡罗近似

I^hT=1/T∑t=1Th(θt)\hat{\mathfrak{I}}_h^T=1 / T \sum_{t=1}^T h\left(\theta_t\right)

减少 I^hT\hat{\mathfrak{I}}_h^T 到 Ih\mathfrak{I}_h 的方差（如果不是收敛速度的话）。

在考虑 Ih\mathfrak{I}_h 的蒙特卡洛近似时，一个常见的综述是积分作为期望的表示不是唯一的（例如 Robert 和 Casella，2004 [79]）。这导致了重要性采样技术，其中可能以自适应方式（Douc 等，2007b [34]）利用替代分布代替 π(θ)\pi(\theta)，或在粒子滤波器中按顺序使用（Del Moral 等，2006 [30]； Andrieu 等，2011 年 [2]）。在本文框架内，给定 MCMC 采样器的结果也可以通过多种方式加以利用，从而改进 \mathgrak{I} h 的近似值。

5.1 Rao-Blackwellization 和其他平均技术

“Rao–Blackwellisation” 这个名称是由 Gelfand 和 Smith (1990 [41]) 在他们的基础 Gibbs 采样论文中创造的，此后它已成为减少积分近似方差的标准方法。虽然它本质上是从基本概率恒等式出发的。

π[h(θ)]=π1[π2h(θ)∣ξ],π[h(θ)] = π1[π2{h(θ)|ξ}],

当 ππ 可以表示为如下的边缘密度

π(θ)=∫Ξπ1(ξ)π2(θ∣ξ)dξ,π(θ) = ∫Ξπ1(ξ)π2(θ|ξ)dξ,

虽然充分性与蒙特卡洛近似没有明确的等价性，但该名称源于 Rao-Blackwell 定理（Lehmann 和 Casella，1998 年 [52]），该定理通过以充分统计为条件来改进给定的估计量。在蒙特卡洛设置中，这意味着方程式 (4) 可以通过部分集成版本进行改进

I Th=1T∑t=1TEπ2[h(θ)∣∣ξt](5)I~Th=1T∑t=1TEπ2[h(θ)||ξt] (5)

假设第二个连接的模拟序列 (ξtξ t) 可用并且条件期望很容易构建。例如，吉布斯采样 (Gelfand & Smith, 1990 [41]) 通常对这种 Rao–Blackwell 分解持开放态度，因为它依赖于来自多个条件分布的连续模拟，可能包括辅助变量和有害参数。特别是，一种称为切片采样器 (Robert & Casella, 2004 [79]) 的通用形式的 Gibbs 采样在每次迭代时产生一个或多个均匀变量。

然而，对于涉及拒绝的所有 MCMC 方法，首先是 Metropolis-Hastings 算法，可以使用更通用的 Rao-Blackwellization 类型（Casella 和 Robert，1996 [25]）。事实上，首先，可以导出或近似拒绝变量的分布，这导致对原始估计量的重要性校正。此外，接受-拒绝步骤依赖于一个统一变量，但这个统一变量可以被整合出来。即，给定一个由 Metropolis–Hastings 算法生成的样本 θ(1),…,θ(T)θ (1), \ldots , θ (T)，可以利用两个基础样本，建议值 ϑ1,…,ϑTϑ 1, \ldots , ϑ T 和统一的 u1,…,uTu 1, \ldots , u T，因此遍历均值可以重写为

ITh=1T∑t=1Th(θ(t))=1T∑t=1Th(θt)∑i=tTIθ(i)=θt。I^Th=1T∑t=1Th(θ(t))=1T∑t=1Th(\theta_t )∑i=tTIθ(i)=\theta_t 。

条件期望

I∼Th=1T∑t=1Th(θt)E[∑i=tTIθ(i)=θt∣θ1,…,θT]=1T∑t=1Th(θt)∑i=tTP(θ(i)=θt∣∣θ1,…,θT)I\sim Th=1T∑t=1Th(\theta_t )E[∑i=tTIθ(i)=\theta_t |θ1,\ldots ,θT]=1T∑t=1Th(\theta_t ){∑i=tTP(θ(i)= \theta_t ||θ1,\ldots ,θT)}

然后享有较小的方差。另请参阅 Tjelmeland (2004 [94]) 以及 Douc 和 Robert (2010 [35]) 基于多次尝试的相关改进。通过在每个 Metropolis-Hastings 迭代中整合决策步骤，可以考虑一个更基本（和更便宜）的版本：如果 θtθ t 是马尔可夫链的当前值并且 ϑtϑ t 是建议值，则被接受（作为 θt+1θ t + 1) 概率为 αtα t，版本

1T∑t=1Tαth(θt)+(1−αt)h(θt)1T∑t=1T{αth(\theta_t )+(1−αt)h(\theta_t )}

应该最常 2 带来对基本估计的改进（Liu 等，1995 年 [55]；Robert 和 Casella，2004 年 [79]）。

6 结论

加速 MCMC 算法听起来像是一个新的悖论，因为有很多方法可以加速给定的算法。然而，这种无 w 限倒退的停止规则是，实现这种加速所带来的额外痛苦可能会在某些时候克服额外的增益。尽管我们在本次调查中仅涵盖了一些可能的方向，而且大多是肤浅的，但我们因此鼓励最热情的读者保持对各种几乎免费的加速解决方案所带来的潜力的认识，并继续尝试设计更多在每个新的 MCMC 实施中进行微调。例如，对于我们中的至少一个人来说，在这个阶段总是考虑 Rao-Blackwellisation。因此，强烈建议至少保留一些这样的技巧供自己使用

参考文献

  • [1] Andrieu, C. and Doucet, A. (2002). Particle filtering for partially observed Gaussian state space models. J. Royal Statist. Society Series B, 64 827–836.
  • [2] Andrieu, C., Doucet, A. and Holenstein, R. (2011). Particle Markov chain Monte Carlo (with discussion). J. Royal Statist. Society Series B, 72 (2) 269–342.
  • [3] Andrieu, C. and Roberts, G. (2009). The pseudo-marginal approach for efficient Monte Carlo computations. The Annals of Statistics 697–725.
  • [4] Angelino, E., Kohler, E., Waterland, A., Seltzer, M. and Adams, R. (2014). Accelerating MCMC via parallel predictive prefetching. arXiv preprint arXiv:1403.7265.
  • [5] Aslett, L., Esperanc ̧a, P. and Holmes, C. (2015). A review of homomorphic encryption and software tools for encrypted statistical machine learning. arXiv:1508.06574.
  • [6] Atchad ́e, Y. F. and Liu, J. S. (2004). The Wang-Landau algorithm for Monte Carlo computation in general state spaces. Statistica Sinica, 20 209–33.
  • [7] Atchad ́e, Y. F., Roberts, G. and Rosenthal, J. (2011). Towards optimal scaling of Metropolis-coupled Markov chain Monte Carlo. Statistics and Computing, 21 555–568.
  • [8] Banterle, M., Grazian, C., Lee, A. and Robert, C. P. (2015). Accelerating Metropolis–Hastings algorithms by delayed acceptance. arXiv preprint arXiv:1503.00996.
  • [9] Bardenet, R., Doucet, A. and Holmes, C. (2014). Towards scaling up Markov chain Monte Carlo: an adaptive subsampling approach 405–413.
  • [10] Bardenet, R., Doucet, A. and Holmes, C. (2015). On Markov chain Monte Carlo methods for tall data. J. Machine Learning Res., 18 1515–1557.
  • [11] B ́edard, M., Douc, R. and Moulines, E. (2012). Scaling analysis of multiple-try MCMC methods. Stochastic Processes and their Applications, 122 758–786.
  • [12] Betancourt, M. (2017). A conceptual introduction to Hamiltonian Monte Carlo. ArXiv e-prints. 1701.02434.
  • [13] Bhatnagar, N. and Randall, D. (2016). Simulated tempering and swapping on mean-field models. Journal of Statistical Physics, 164 495–530.
  • [14] Bierkens, J. (2016). Non-reversible Metropolis–Hastings. Statistics and Computing, 26 1213–1228.
  • [15] Bierkens, J., Bouchard-Cˆot ́e, A., Doucet, A., Duncan, A. B., Fearnhead, P., Roberts,G. and Vollmer, S. J. (2017). Piecewise deterministic Markov processes for scalable Monte Carlo on restricted domains. arXiv preprint arXiv:1701.04244.
  • [16] Bierkens, J., Fearnhead, P. and Roberts, G. (2016). The zig-zag process and super-efficient sampling for Bayesian analysis of big data. arXiv preprint arXiv:1607.03188.
  • [17] Bou-Rabee, N., Sanz-Serna, J. M. et al. (2017). Randomized hamiltonian monte carlo. The Annals of Applied Probability, 27 2159–2194.
  • [18] Bouchard-Cˆot ́e, A., Vollmer, S. J. and Doucet, A. (2017). The bouncy particle sampler: A non-reversible rejection-free Markov chain Monte Carlo method. Journal of the American Statistical Association.
  • [19] Calderhead, B. (2014). A general construction for parallelizing Metropolis–Hastings algorithms. Proceedings of the National Academy of Sciences, 111 17408–17413.
  • [20] Capp ́e, O., Douc, R., Guillin, A., Marin, J.-M. and Robert, C. (2008). Adaptive importance sampling in general mixture classes. Statist. Comput., 18 447–459.
  • [21] Capp ́e, O., Guillin, A., Marin, J.-M. and Robert, C. (2004). Population Monte Carlo. J. Comput. Graph. Statist., 13 907–929.
  • [22] Capp ́e, O. and Robert, C. (2000). Markov Chain Monte Carlo: Ten years and still running! J. American Statist. Assoc., 95 1282–1286.
  • [23] Carpenter, B., Gelman, A., Hoffman, M., Lee, D., Goodrich, B., Betancourt, M., Brubaker, M., Guo, J., Li, P. and Riddell, A. (2017). Stan: A probabilistic programming language. Journal of Statistical Software, Articles, 76.
  • [24] Carter, J. and White, D. (2013). History matching on the imperial college fault model using parallel tempering. Computational Geosciences, 17 43–65.
  • [25] Casella, G. and Robert, C. (1996). Rao-Blackwellization of sampling schemes. Biometrika, 83 81–94.
  • [26] Chen, T., Fox, E. and Guestrin, C. (2014). Stochastic gradient hamiltonian monte carlo 1683–1691.
  • [27] Chen, T. and Hwang, C. (2013). Accelerating reversible Markov chains. Statistics & Probability Letters, 83 1956–1962.
  • [28] Davis, M. H. (1984). Piecewise-deterministic Markov processes: A general class of non-diffusion stochastic models. Journal of the Royal Statistical Society. Series B (Methodological) 353–388.
  • [29] Davis, M. H. (1993). Markov Models & Optimization, vol. 49. CRC Press.
  • [30] Del Moral, P., Doucet, A. and Jasra, A. (2006). Sequential Monte Carlo samplers. J. Royal Statist. Society Series B, 68 411–436.
  • [31] Deligiannidis, G., Doucet, A. and Pitt, M. K. (2015). The correlated pseudo-marginal method. arXiv preprint arXiv:1511.04992.
  • [32] Ding, N., Fang, Y., Babbush, R., Chen, C., Skeel, R. D. and Neven, H. (2014). Bayesian sampling using stochastic gradient thermostats 3203–3211.
  • [33] Douc, R., Guillin, A., Marin, J.-M. and Robert, C. (2007a). Convergence of adaptive mixtures of importance sampling schemes. Ann. Statist., 35(1) 420–448.
  • [34] Douc, R., Guillin, A., Marin, J.-M. and Robert, C. (2007b). Convergence of adaptive mixtures of importance sampling schemes. Ann. Statist., 35(1) 420–448. ArXiv:0708.0711.
  • [35] Douc, R. and Robert, C. (2010). A vanilla variance importance sampling via population Monte Carlo. Ann. Statist. To appear.
  • [36] Doucet, A., Godsill, S. and Andrieu, C. (2000). On sequential Monte-Carlo sampling methods for Bayesian filtering. Statist. Comp., 10 197–208.
  • [37] Duane, S., Kennedy, A. D., Pendleton, B. J. and Roweth, D. (1987). Hybrid Monte Carlo. Phys. Lett. B, 195 216–222.
  • [38] Durmus, A. and Moulines, E. (2017). Nonasymptotic convergence analysis for the unadjusted Langevin algorithm. Ann. Applied Probability, 27 1551–1587.
  • [39] Earl, D. J. and Deem, M. W. (2005). Parallel tempering: Theory, Applications, and New Perspectives. Physical Chemistry Chemical Physics, 7 3910–3916.
  • [40] Fielding, M., Nott, D. J. and Liong, S.-Y. (2011). Efficient MCMC schemes for computationally expensive posterior distributions. Technometrics, 53 16–28. https://doi.org/10.1198/TECH.2010.09195, URL https://doi.org/10.1198/TECH.2010.09195.
  • [41] Gelfand, A. and Smith, A. (1990). Sampling based approaches to calculating marginal densities. J. American Statist. Assoc., 85 398–409.
  • [42] Gelman, A., Gilks, W. and Roberts, G. (1996). Efficient Metropolis jumping rules. In Bayesian Statistics 5 (J. Berger, J. Bernardo, A. Dawid, D. Lindley and A. Smith, eds.). Oxford University Press, Oxford, 599–608.
  • [43] Geyer, C. J. (1991). Markov chain Monte Carlo maximum likelihood. Computing Science and Statistics, 23 156–163.
  • [44] Girolami, M. and Calderhead, B. (2011). Riemann manifold Langevin and Hamiltonian Monte Carlo methods. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73 123–214.
  • [45] Guihenneuc-Jouyaux, C. and Robert, C. P. (1998). Discretization of continuous Markov chains and Markov chain Monte Carlo convergence assessment. Journal of the American Statistical Association, 93 1055–1067.
  • [46] Haario, H., Saksman, E. and Tamminen, J. (1999). Adaptive proposal distribution for random walk Metropolis algorithm. Computational Statistics, 14(3) 375–395.
  • [47] Hoffman, M. D. and Gelman, A. (2014). The No-U-turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. J. Machine Learning Res., 15 1593–1623.
  • [48] Hwang, C.-R., Hwang-Ma, S.-Y. and Sheu, S.-J. (1993). Accelerating gaussian diffusions. The Annals of Applied Probability 897–913.
  • [49] Iba, Y. (2000). Population-based Monte Carlo algorithms. Trans. Japanese Soc. Artificial Intell., 16 279–286.
  • [50] Jacob, P. E., O’Leary, J. and Atchad ́e, Y. F. (2017). Unbiased Markov chain Monte Carlo with couplings. ArXiv e-prints. 1708.03625.
  • [51] Jacob, P., Robert, C. P. and Smith, M. H. (2011). Using parallel computation to improve independent Metropolis–Hastings based estimation. Journal of Computational and Graphical Statistics, 20 616–635.
  • [52] Lehmann, E. and Casella, G. (1998). Theory of Point Estimation (revised edition). Springer-Verlag, New York.
  • [53] Liang, F., Liu, C. and Carroll, R. (2007). Stochastic approximation in Monte Carlo computation. JASA, 102 305–320.
  • [54] Liu, J., Wong, W. and Kong, A. (1994). Covariance structure of the Gibbs sampler with application to the comparison of estimators and augmentation schemes. Biometrika, 81 27–40.
  • [55] Liu, J., Wong, W. and Kong, A. (1995). Covariance structure and convergence rates of the Gibbs sampler with various scans. Journal of Royal Statistical Society, B 57 157–169.
  • [56] Liu, J. S., Liang, F. and Wong, W. H. (2000). The multiple-try method and local optimization in Metropolis sampling. J. American Statist. Assoc., 95 121–134.
  • [57] Livingstone, S., Faulkner, M. F. and Roberts, G. O. (2017). Kinetic energy choice in hamiltonian/hybrid monte carlo. arXiv preprint arXiv:1706.02649.
  • [58] MacKay, D. J. C. (2002). Information Theory, Inference & Learning Algorithms. Cambridge University Press, Cambridge, UK.
  • [59] Marinari, E. and Parisi, G. (1992). Simulated tempering: a new Monte Carlo scheme. EPL (Europhysics Letters), 19 451.
  • [60] Martino, L. (2018). A Review of Multiple Try MCMC algorithms for Signal Processing. ArXiv e-prints. 1801.09065.
  • [61] Martino, L., Elvira, V., Luengo, D., Corander, J. and Louzada, F. (2016). Orthogonal parallel MCMC methods for sampling and optimization. Digital Signal Processing, 58 64–84.
  • [62] Mengersen, K. and Robert, C. (2003). Iid sampling with self-avoiding particle filters: the pinball sampler. In Bayesian Statistics 7 (J. Bernardo, M. Bayarri, J. Berger, A. Dawid, D. Heckerman, A. Smith and M. West, eds.). Oxford University Press, Oxford.
  • [63] Meyn, S. and Tweedie, R. (1993). Markov Chains and Stochastic Stability. Springer-Verlag, New York.
  • [64] Miasojedow, B., Moulines, E. and Vihola, M. (2013). An adaptive parallel tempering algorithm. Journal of Computational and Graphical Statistics, 22 649–664.
  • [65] Minsker, S., Srivastava, S., Lin, L. and Dunson, D. (2014). Scalable and robust Bayesian inference via the median posterior 1656–1664.
  • [66] Mira, A. (2001). ”on metropolis-hastings algorithms with delayed rejection”. Metron, 59 (3-4) 231–241.
  • [67] Mira, A. and Sargent, D. J. (2003). A new strategy for speeding Markov chain Monte Carlo algorithms. Stat. Methods Appl., 12 49–60.
  • [68] Mohamed, L., Calderhead, B., Filippone, M., Christie, M. and Girolami, M. (2012). Population mcmc methods for history matching and uncertainty quantification. Computational Geosciences, 16 423–436.
  • [69] Mykland, P., Tierney, L. and Yu, B. (1995). Regeneration in Markov chain samplers. Journal of the American Statistical Association, 90 233–241.
  • [70] Neal, R. (1999). Bayesian Learning for Neural Networks, vol. 118. Springer–Verlag, New York. Lecture Notes.
  • [71] Neal, R. (2011a). MCMC using Hamiltonian dynamics. In In Handbook of Markov Chain Monte Carlo (S. Brooks, A. Gelman, G. L. Jones and X.-L. Meng, eds.). CRC Press, New York.
  • [72] Neal, R. M. (1996). Sampling from Multimodal Distributions using Tempered Transitions. Statistics and Computing, 6 353–366.
  • [73] Neal, R. M. (2011b). MCMC Using Ensembles of States for Problems with Fast and Slow Variables such as Gaussian Process Regression. ArXiv e-prints. 1101.0387.
  • [74] Neiswanger, W., Wang, C. and Xing, E. (2013). Asymptotically exact, embarrassingly parallel MCMC. arXiv preprint arXiv:1311.4780.
  • [75] Quiroz, M., Villani, M. and Kohn, R. (2016). Exact subsampling MCMC. arXiv preprint arXiv:1603.08232.
  • [76] Rasmussen, C. E. (2003). Gaussian processes to speed up hybrid monte carlo for expensive bayesian integrals. In Bayesian Statistics 7 (J. Bernardo, M. Bayarri, J. Berger, A. Dawid, D. Heckerman, A. Smith and M. West, eds.). Oxford University Press, 651–659.
  • [77] Rasmussen, C. E. and Williams, C. K. I. (2005). Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning). The MIT Press.
  • [78] Rhee, C.-h. and Glynn, P. W. (2015). Unbiased estimation with square root convergence for sde models. Operations Research, 63 1026–1043.
  • [79] Robert, C. and Casella, G. (2004). Monte Carlo Statistical Methods. 2nd ed. Springer-Verlag, New York.
  • [80] Robert, C. and Casella, G. (2009). Introducing Monte Carlo Methods with R. Springer-Verlag, New York.
  • [81] Roberts, G., Gelman, A. and Gilks, W. R. (1997). Weak convergence and optimal scaling of random walk Metropolis algorithms. The Annals of Applied Probability, 7 110–120.
  • [82] Roberts, G. and Rosenthal, J. (2001). Optimal scaling for various Metropolis-Hastings algorithms. Statistical Science, 16 351–367.
  • [83] Roberts, G. and Rosenthal, J. (2007). Coupling and ergodicity of adaptive Markov Chain Monte carlo algorithms. J. Applied Proba., 44(2) 458–475.
  • [84] Roberts, G. and Rosenthal, J. (2009). Examples of adaptive MCMC. J. Comp. Graph. Stat., 18 349–367.
  • [85] Roberts, G. and Rosenthal, J. (2014). Minimising MCMC variance via diffusion limits, with an application to simulated tempering. The Annals of Applied Probability, 24 131–149.
  • [86] Rubinstein, R. Y. (1981). Simulation and the Monte Carlo Method. J. Wiley, New York.
  • [87] Saksman, E. and Vihola, M. (2010). On the ergodicity of the adaptive Metropolis algorithm on unbounded domains. Ann. Applied Probability, 20(6) 2178–2203.
  • [88] Scott, S. L., Blocker, A. W., Bonassi, F. V., Chipman, H. A., George, E. I. and McCulloch, R. E. (2016). Bayes and big data: The consensus Monte Carlo algorithm. International Journal of Management Science and Engineering Management, 11 78–88.
  • [89] Srivastava, S., Cevher, V., Dinh, Q. and Dunson, D. (2015). Wasp: Scalable Bayes via barycenters of subset posteriors 912–920.
  • [90] Storvik, G. (2002). Particle filters for state space models with the presence of static parameters. IEEE Trans. Signal Process., 50 281–289.
  • [91] Sun, Y., Schmidhuber, J. and Gomez, F. J. (2010). Improving the asymptotic performance of Markov chain Monte-carlo by inserting vortices 2235–2243.
  • [92] Terenin, A., Simpson, D. and Draper, D. (2015). Asynchronous Gibbs Sampling. ArXiv e-prints. 1509.08999.
  • [93] Tierney, L. and Mira, A. (1998). Some adaptive Monte Carlo methods for Bayesian inference. Statistics in Medicine, 18 2507–2515.
  • [94] Tjelmeland, H. (2004). Using all Metropolis–Hastings proposals to estimate mean values. Tech. Rep. 4, Norwegian University of Science and Technology, Trondheim, Norway.
  • [95] Vehtari, A., Gelman, A., Sivula, T., Jyl ̈anki, P., Tran, D., Sahai, S., Blomstedt, P., Cunningham, J. P., Schiminovich, D. and Robert, C. P. (2014). Expectation propagation as a way of life: A framework for Bayesian inference on partitioned data. ArXiv e-prints. 1412.4869.
  • [96] Wang, F. and Landau, D. (2001). Determining the density of states for classical statistical models: A random walk algorithm to produce a flat histogram. Physical Review E, 64 056101.
  • [97] Wang, X. and Dunson, D. (2013). Parallelizing MCMC via weierstrass sampler. arXiv preprint arXiv:1312.4605.
  • [98] Wang, X., Guo, F., Heller, K. and Dunson, D. (2015). Parallelizing MCMC with random partition trees. Advances in Neural Information Processing Systems 451–459.
  • [99] Welling, M. and Teh, Y. (2011). Bayesian learning via stochastic gradient Langevin dynamics 681–688.
  • [100] Woodard, D. B., Schmidler, S. C. and Huber, M. (2009a). Conditions for rapid mixing of parallel and simulated tempering on multimodal distributions. The Annals of Applied Probability 617–640.
  • [101] Woodard, D. B., Schmidler, S. C. and Huber, M. (2009b). Sufficient conditions for torpid mixing of parallel and simulated tempering. Electronic Journal of Probability, 14 780–804.
  • [102] Xie, Y., Zhou, J. and Jiang, S. (2010). Parallel tempering monte carlo simulations of lysozyme orientation on charged surfaces. The Journal of chemical physics, 132 02B602.