Unbiased Markov chain Monte Carlo with couplings

Introduction

Markov chain Monte Carlo (MCMC) methods constitute a popular class of algorithms to approximate high-dimensional integrals arising in statistics and other fields (Liu 2008; Robert and Casella 2004; Brooks et al. 2011; Green et al. 2015). These iterative methods provide estimators that are consistent as the number of iterations grows large but potentially biased for any fixed number of iterations, which discourages the parallel execution of many short chains (Rosenthal 2000). Consequently, efforts have focused on exploiting parallel processors within each iteration (Tjelmeland 2004; Brockwell 2006; Lee et al. 2010; Jacob et al. 2011; Calderhead 2014; Goudie et al. 2017; Yang et al. 2017) and on the design of parallel chains targeting different distributions (Altekar et al. 2004; Wang et al. 2015; Srivastava et al. 2015). Still, MCMC estimators are ultimately justified by asymptotics in the number of iterations, which is discordant with current trends in computing hardware, characterized by increasing parallelism but stagnating clock speeds.

In this paper we propose a general construction to produce unbiased estimators of integrals with respect to a target probability distribution from MCMC kernels. The lack of bias means that these estimators can be implemented on parallel processors in the framework of Glynn and Heidelberger 1991, without communication between processors. Confidence intervals can be constructed with asymptotic guarantees in the number of processors, in contrast with standard MCMC confidence intervals that are justified asymptotically in the number of iterations (Flegal et al. 2008; Gong and Flegal 2016; Atchadé 2016; Vats et al. 2018, e.g.). The lack of bias has additional benefits, as discussed in Section 5.5 in which we make use of its interplay with the law of iterated expectations to perform modular inference; see also the discussion in Section 6.

Our contribution follows the path-breaking work of Glynn and Rhee 2014, which uses couplings to construct unbiased estimators of integrals with respect to an invariant distribution. They illustrate their construction on Markov chains represented by iterated random functions, leveraging the contraction properties of such functions. Glynn and Rhee 2014 also consider Harris recurrent chains for which an explicit minorization condition holds. Previously, McLeish 2011 employed similar debiasing techniques to obtain “nearly unbiased” estimators from a single MCMC chain. More recently Jacob et al. 2019 remove the bias from conditional particle filters (Andrieu et al. 2010) by coupling chains so that they meet in finite time. The present article brings this type of “Rhee–Glynn” construction to generic MCMC algorithms, with a novel analysis of estimator efficiency and a variety of examples. Our proposed construction involves couplings of MCMC algorithms, which we discuss for generic Metropolis–Hastings and Gibbs samplers.

Couplings have been used to study the convergence properties of MCMC algorithms from both theoretical and practical points of view (Reutter and Johnson 1995; Johnson 1996; Rosenthal 1997; Johnson 1998; Neal 1999; Roberts and Rosenthal 2004; Johnson 2013; Johndrow and Mattingly 2017, e.g.). Couplings also underpin perfect samplers (Propp and Wilson 1996; Murdoch and Green 1998; Casella et al. 2001; Flegal and Herbei 2012; Lee et al. 2014; Huber 2016). A notable aspect of the approach of Glynn and Rhee 2014 preserved in our method is that only two chains have to be coupled for the proposed estimator to be unbiased, without further assumptions on the state space or target distribution. Thus the approach applies more broadly than perfect samplers (Glynn 2016, see) while yielding unbiased estimators rather than exact samples. Coupling pairs of Markov chains also forms the basis of the approach of Neal 1999, with a similar motivation for parallel computation. The proposed estimation technique also shares aims with regeneration methods (Mykland et al. 1995; Brockwell and Kadane 2005, e.g.), and we propose a numerical comparison in Section 5.2.

In Section 2 we introduce our estimators and present a coupling of random walk Metropolis–Hastings chains as an illustration. In Section 3 we establish the efficiency properties of these estimators, discuss the verification of key assumptions, and describe the use of the proposed estimators on parallel processors in light of results from e.g. Glynn and Heidelberger 1991. In Section 4 we describe how to couple some important MCMC algorithms and illustrate the effect of dimension on algorithm performance with a multivariate Normal target. Section 5 contains more challenging examples including a multimodal target, a comparison with regeneration methods, sampling problems in large-dimensional discrete spaces arising in Bayesian variable selection and Ising models, and an application to modular inference. We discuss our findings in Section 6. Scripts in R (R Core Team 2015) are available at https://github.com/pierrejacob/unbiasedmcmc and supplementary materials are available online.

Unbiased estimation from coupled chains

The chains stay together after meeting, i.e. $X_{t}=Y_{t-1}$ for all $t\geq\tau$.

By construction, each of the marginal chains $(X_{t})_{t\geq 0}$ and $(Y_{t})_{t\geq 0}$ has initial distribution $\pi_{0}$ and transition kernel $P$. Assumption 2.1 requires these chains to result in a uniformly bounded $(2+\eta)$-moment of $h$; more discussion on moments of Markov chains can be found in Tweedie 1983. Since $X_{0}$ and $Y_{0}$ may be drawn from any coupling of $\pi_{0}$ with itself, it is possible to set $X_{0}=Y_{0}$. However, $X_{1}$ is then generated from $P(X_{0},\cdot)$, so that $X_{1}\neq Y_{0}$ in general. Thus one cannot force the meeting time to be small by setting $X_{0}=Y_{0}$. Assumption 2.2 puts a condition on the coupling operated by $\bar{P}$, and would not in general be satisfied for an independent coupling. Coupled kernels must be carefully designed, using e.g. common random numbers and maximal couplings, for Assumption 2.2 to be satisfied. We present a simple case in Section 2.2 and further examples in Section 4. We stress that the state space is not assumed to be discrete, and that the constants $D$ and $\eta$ of Assumption 2.1 and $C$ and $\delta$ of Assumption 2.2 do not need to be known to implement the proposed approach. Assumption 2.3 typically holds by design; coupled chains that stay identical after meeting are termed “faithful” in Rosenthal 1997.

Before presenting examples and enhancements to the estimator above, we discuss the relationship between our approach and existing work. There is a rich literature applying forward couplings to study Markov chains convergence (Johnson 1996; Johnson 1998; Thorisson 2000; Lindvall 2002; Rosenthal 2002; Johnson 2013; Douc et al. 2004; Nikooienejad et al. 2016), and to obtain new algorithms such as perfect samplers (Huber 2016) and the methods of Neal 1999 and Neal and Pinto 2001. Our approach is closely related to Glynn and Rhee 2014, who employ pairs of Markov chains to obtain unbiased estimators. The present work combines similar arguments with couplings of MCMC algorithms and proposes further improvements to remove bias at a reduced loss of efficiency.

Indeed Glynn and Rhee 2014 did not apply their methodology to the MCMC setting. They consider chains associated with contractive iterated random functions (Diaconis and Freedman 1999, see also), and Harris recurrent chains with an explicit minorization condition. A minorization condition refers to a small set $\mathcal{C}$, $\lambda>0$, an integer $m\geq 1$, and a probability measure $\nu$ such that for all $x\in\mathcal{C}$ and some measurable set $A$, $P^{m}(x,A)\geq\lambda\nu(A)$. Such a condition is said to be explicit if the set, constant and probability measure are known by the user. Finding explicit small sets that are useful in practice can present a technical challenge, even for MCMC experts (Cowles and Rosenthal 1998, see discussion and references in). When available, explicit minorization conditions can also be employed to identify regeneration times, yielding estimators amenable to parallel computation in the framework of Mykland et al. 1995 and Brockwell and Kadane 2005. By contrast Johnson 1996; Johnson 1998 and Neal 1999 address the question of coupling MCMC algorithms so that pairs of chains meet exactly, without analytical knowledge on the target distribution. The present article focuses on the use of couplings of this type in the framework of Glynn and Rhee 2014.

2 Coupled Metropolis–Hastings example

Before further examination of our estimator and its properties, we present a coupling of Metropolis–Hastings (MH) chains that will typically satisfy Assumptions 2.1-2.3 in realistic settings; this coupling was proposed in Johnson 1998 as part of a method to diagnose convergence. We postpone discussion of other couplings of MCMC algorithms to Section 4. We recall that each iteration $t$ of the MH algorithm (Hastings 1970) begins by drawing a proposal $X^{\star}$ from a Markov kernel $q(X_{t},\cdot)$, where $X_{t}$ is the current state. The next state is set to $X_{t+1}=X^{\star}$ if ${U\leq\pi(X^{\star})q(X^{\star},X_{t})}/({\pi(X_{t})q(X_{t},X^{\star})})$, where $U$ denotes a uniform random variable on $$, and$X_{t+1}=X_{t}$ otherwise.

We define a pair of chains so that each proceeds marginally according to the MH algorithm and jointly so that the chains will meet exactly after a random number of steps. We suppose that the pair of chains are in states $X_{t}$ and $Y_{t-1}$, and consider how to generate $X_{t+1}$ and $Y_{t}$ so that $\{X_{t+1}=Y_{t}\}$ might occur.

If $X_{t}\neq Y_{t-1}$, the event $\{X_{t+1}=Y_{t}\}$ cannot occur if both chains reject their respective proposals, $X^{\star}$ and $Y^{\star}$. Meeting will occur if these proposals are identical and if both are accepted. Marginally, the proposals follow ${X^{\star}|X_{t}\sim q(X_{t},\cdot)}$ and $Y^{\star}|Y_{t-1}\sim q(Y_{t-1},\cdot)$. If $q(x,x^{\star})$ can be evaluated for all $x,x^{\star}$, then one can sample from a maximal coupling between the two proposal distributions, which is a coupling of $q(X_{t},\cdot)$ and $q(Y_{t-1},\cdot)$ maximizing the probability of the event $\{X^{\star}=Y^{\star}\}$. How to sample from maximal couplings of continuous distributions is described in Thorisson 2000 and in Section 4.1. One can accept or reject the two proposals using a common uniform random variable $U$. The chains will stay together after they meet: at each step after meeting, the proposals will be identical with probability one, and jointly accepted or rejected with a common uniform variable. This coupling requires neither explicit minorization conditions nor contractive properties of a random function representation of the chain.

3 Time-averaged estimator

The estimator $H_{k:m}(X,Y)$ requires $\tau-1$ calls to $\bar{P}$ and $\max(1,m+1-\tau)$ calls to $P$, which is overall comparable to $m$ calls to $P$ when $m$ is large. Indeed, for the proposed couplings, calls to $\bar{P}$ are approximately twice as expensive as calls to $P$. Therefore, the cost of $H_{k:m}(X,Y)$ is comparable to $2(\tau-1)+\max(1,m+1-\tau)$ iterations of the underlying MCMC algorithm. Thus both the variance and the cost of $H_{k:m}(X,Y)$ will approach those of MCMC estimators for large values of $k$ and $m$. This motivates the use of the estimator $H_{k:m}(X,Y)$ with $m>k$, which allows us to control the loss of efficiency associated with the removal of burn-in bias in contrast with the basic estimator $H_{k}(X,Y)$ of Section 2.1. We discuss the choice of $k$ and $m$ in further detail in Section 3 and in the subsequent experiments. A variant of (2.1) can be obtained by considering a time lag greater than one between the two chains $(X_{t})_{t\geq 0}$ and $(Y_{t})_{t\geq 0}$, with the meeting time defined as the first time $t$ for which $\{X_{t}=Y_{t-\text{lag}}\}$ occurs. This introduces another tuning parameter but is found to be fruitful in Biswas and Jacob 2019.

We conclude this section with a few remarks on practical implementations. First, the test function $h$ does not have to be specified at run-time in Algorithm 1. One can store the coupled chains and choose the test function later. Also, one typically resorts to thinning the output of an MCMC sampler if the memory cost of storing chains is prohibitive, or if the cost of evaluating the test function of interest is significant compared to the cost of each MCMC iteration (Owen 2017, e.g.). This is feasible in the proposed framework: one could consider a variation of Algorithm 1 where each call to the Markov kernels $P$ and $\bar{P}$ would be replaced by multiple calls to them. We also observe that the proposed estimators can take values outside of the range of the test function $h$; for instance they can take negative values even if the range of the test function contains only non-negative values.

Finally, we stress the difficulty inherent in choosing an initial distribution $\pi_{0}$. The estimators are unbiased for any choice of $\pi_{0}$, including point masses, but this choice has an impact on both the computing cost and the variance. There is also a choice about whether to draw $X_{0}$ and $Y_{0}$ independently from $\pi_{0}$ or not; in our experiments we use independent draws. We will see in Section 5.1 that unfortunate choices of initial distributions can severely affect the performance of the proposed estimators. This suggests trying more than one choice of initialization, especially in the setting of multimodal targets. Overall the choice of $\pi_{0}$ and its relative importance compared to standard MCMC are open questions.

4 Signed measure estimator

We can formulate the proposed estimation procedure in terms of a signed measure $\hat{\pi}$ defined by

Properties and parallel implementation

The proofs of the results of this section are in the supplementary materials. Our first result establishes the basic validity of the proposed estimators.

Section 3.1 studies the variance and efficiency of $H_{k:m}(X,Y)$, Section 3.2 concerns the verification of Assumption 2.2 using drift conditions, and Section 3.3 discusses estimation on parallel processors in the presence of a budget constraint.

We consider the impact of $k$ and $m$ on the efficiency of the proposed estimators, which will then suggest guidelines for the choice of these tuning parameters. Estimators $H^{(r)}_{k:m}(X,Y)$, for $r=1,\ldots,R$, can be generated independently and averaged. More estimators can be produced in a given computing budget if each estimator is cheaper to produce. The trade-off can be understood in the framework of Glynn and Whitt 1992, see also Rhee and Glynn 2012; Glynn and Rhee 2014, by defining the asymptotic inefficiency as the product of the variance and expected cost of the estimator. That product is the asymptotic variance of $R^{-1}\sum_{r=1}^{R}H^{(r)}_{k:m}(X,Y)$ as the computational budget, as opposed to the number of estimators $R$, goes to infinity (Glynn and Whitt 1992). Of primary interest is the comparison of this asymptotic inefficiency with the asymptotic variance of standard MCMC estimators. We start by writing the time-averaged estimator of (2.1) as

The Markov kernel $P$ is $\pi$-invariant, $\varphi$-irreducible and aperiodic, and there exists a measurable function $V:\;\mathcal{X}\to[1,\infty)$, $\lambda\in(0,1)$, $b<\infty$ and a small set $\mathcal{C}$ such that for all $x\in\mathcal{X}$,

Suppose that Assumptions 2.2-2.3 and 3.1 hold, with a function $V$ for which the integral $\int V(x)\pi_{0}(dx)$ is finite. If the function $h$ is such that $\sup_{x\in\mathcal{X}}|h(x)|/V(x)^{\beta}<\infty$ for some $\beta\in[0,1/2)$, then for all $m\geq k\geq 0$ we have

for some constants $C_{\delta,\beta}<+\infty$, and $\delta_{\beta}=\delta^{1-2\beta}\in(0,1)$, with $\delta\in(0,1)$ as in Assumption 2.2.

Using Proposition 3.3, equation (3.1) becomes

The variance of $H_{k:m}(X,Y)$ is thus bounded by the mean squared error of an MCMC estimator plus additive terms that vanish geometrically in $k$ and polynomially in $m-k$.

Dropping the third term on the right-hand side of (3.2), which is of smaller magnitude than the second term, assuming that $\text{MSE}_{k:m}>0$ and that $m>\tau$ with large probability, we obtain the approximate inequality

This informal series of approximations suggests that we can retrieve an asymptotic efficiency comparable to the underlying MCMC estimators with appropriate choices of $k$ and $m$ that depend on the distribution of the meeting time $\tau$. These choices are thus sensitive to the coupling of the chains, and not only to the performance of the underlying MCMC algorithm. Choosing $m$ as a multiple of $k$, such as $5k$ or $10k$, makes intuitive sense when considering that $k/m$ is the proportion of iterations that are simply discarded in the event that $\tau<k$. In other words, the bias of MCMC can be removed at the cost of an increased variance, which can in turn be reduced by choosing large enough values of $k$ and $m$. This results in a tradeoff with the desired level of parallelism: one might prefer to keep $k$ and $m$ small, yielding a suboptimal efficiency for $H_{k:m}(X,Y)$, but enabling more independent copies to be generated in a given computing time.

2 Verifying Assumption 2.2

We discuss how Assumption 3.1 on the Markov kernel $P$ can be used to verify Assumption 2.2, on the shape of the meeting time distribution. Informally, Assumption 3.1 guarantees that the bivariate chain $\{(X_{t},Y_{t-1}),\;t\geq 1\}$ visits $\mathcal{C}\times\mathcal{C}$ infinitely often, where $\mathcal{C}$ is a small set. If there is a positive probability of the event $\{X_{t+1}=Y_{t}\}$ for every $t$ such that $(X_{t},Y_{t-1})\in\mathcal{C}\times\mathcal{C}$, then we expect Assumption 2.2 to hold. The next result formalizes that intuition. The proof is based on a modification of an argument by Douc et al. 2004. We introduce $\mathcal{D}=\{(x,y)\in\mathcal{X}\times\mathcal{X}:\;x=y\}$. Then Assumption 2.3 reads $\bar{P}((x,x),\mathcal{D})=1$ for all $x\in\mathcal{X}$.

Suppose that $P$ satisfies Assumption 3.1 with a small set $\mathcal{C}$ of the form $\mathcal{C}=\{x:V(x)\leq L\}$ where $\lambda+b/(1+L)<1$. Suppose also that there exists $\epsilon\in(0,1)$ such that

Then there exists a finite constant $C^{\prime}$ and a $\kappa\in(0,1)$, such that for all $n\geq 1$,

where $\pi_{0}(V)=\int V(x)\pi_{0}(dx)$. Hence Assumption 2.2 holds as long as $\pi_{0}(V)<\infty$.

Note that if Assumption 3.1 holds with a small set of the form $\mathcal{C}=\{x:\;V(x)\leq L\}$ for some $L>0$, then it also holds for $\mathcal{C}=\{x:\;V(x)\leq L^{\prime}\}$ for all $L^{\prime}\geq L$. In that case one can always choose $L$ large enough so that $\lambda+b/(1+L)<1$. Hence the main restriction in Proposition 3.4 is the assumption that the small sets in Assumption 3.1 are of the form $\{x:V(x)\leq L\}$, i.e. level sets of $V$. This is known to be true in some cases. For instance it is known from Theorem 2.2 of Roberts and Tweedie 1996b that for a large class of Metropolis-Hastings algorithms, any non-empty compact set is a small set, and therefore for these algorithms it suffices to check that the level sets of the drift function $V$ are compact. Common examples of drift functions include $V(x)=c/\sqrt{\pi(x)}$ (Roberts and Tweedie 1996b; Jarner and Hansen 2000; Atchade 2006), $V(x)=ce^{b|x|}$ (Roberts and Tweedie 1996a) or the example in Pal and Khare 2014, which all have compact level sets under mild regularity conditions.

The work of Middleton et al. 2018 contains results that generalize Propositions 3.3 and 3.4 to Markov chains satisfying polynomial drift conditions (Andrieu and Vihola 2015, e.g), leading to polynomial tails for the associated meeting times.

3 Parallel implementation under budget constraints

Our main motivation for unbiased estimators comes from parallel processing; see Sections 5.5 and 6 for other motivations. Independent unbiased estimators with finite variance can be generated on separate machines, and combined into consistent and asymptotically Normal estimators. If the number of estimators is pre-specified, this follows from the central limit theorem for i.i.d. variables. We might prefer to specify a time budget, and generate as many estimators as possible within the budget. The lack of bias allows the application of a variety of results on budget-constrained parallel simulations, which we briefly review here, following Glynn and Heidelberger 1990; Glynn and Heidelberger 1991.

Couplings of MCMC algorithms

We consider couplings of various MCMC algorithms that satisfy Assumptions 2.2-2.3. These couplings are widely applicable and do not require extensive analytical knowledge of the target distribution. We stress that they are not optimal in general, and we expect that other constructions would yield more efficient estimators. We begin in Section 4.1 by reviewing maximal couplings.

where $d_{\text{TV}}(p,q)=\nicefrac{{1}}{{2}}\int_{\mathcal{X}}|p(x)-q(x)|dx$ is the total variation distance. By the coupling inequality (Lindvall 2002), this proves that the algorithm implements a maximal coupling.

Let $z=\Sigma^{-1/2}(\mu_{1}-\mu_{2})$ and $e=z/\|z\|$. We independently draw $\dot{X}\sim s$ and $U\sim\mathcal{U}()$ and let

The above procedure outputs a pair $(\dot{X},\dot{Y})$ that follows a coupling of $s$ with itself. We then define $(X,Y)=(\mu_{1}+\Sigma^{1/2}\dot{X},\mu_{2}+\Sigma^{1/2}\dot{Y})$. On the event $\{\dot{Y}=\dot{X}+z\}$, we have $X=Y$. On the event $\{\dot{Y}\neq\dot{X}+z\}$, the vector $\dot{X}-2(e^{\prime}\dot{X})e$ is the reflection of $\dot{X}$ through the hyperplane orthogonal to $e$ that passes through the origin. We show that the output $(X,Y)$ follows a maximal coupling of $p$ and $q$, which we refer to as a maximal coupling with reflection on the residuals, or a “reflection-maximal coupling”. First we show that $\dot{Y}$ follows $s$, closely following the argument in Bou-Rabee et al. 2018. For a measurable set $B$, we compute

The first integral above becomes $\int\mathds{1}_{B}(w)\min\left(s(w-z),s(w)\right)dw$, after a change of variables $w:=x+z$. To simplify the second integral we make the change of variables $w:=x-2(e^{\prime}x)e$. Since this corresponds to a reflection with respect to a plane orthogonal to $e$, we have $dw=dx$, and $x=w-2(e^{\prime}w)e$, thus

To verify that the procedure corresponds to a maximal coupling of $p$ and $q$, we observe that

Finally, for discrete distributions with common finite support, a procedure for sampling from a maximal coupling is described in Section 5.4, with a cost that is also deterministic.

2 Metropolis–Hastings

In Section 2.2 we described a coupling of MH chains due to Johnson 1998; we summarize the coupled kernel $\bar{P}((X_{t},Y_{t-1}),\cdot)$ in the following procedure.

Sample $(X^{\star},Y^{\star})|(X_{t},Y_{t-1})$ from a maximal coupling of $q(X_{t},\cdot)$ and $q(Y_{t-1},\cdot)$.

If $U\leq\min(1,\pi(X^{\star})q(X^{\star},X_{t})/\pi(X_{t})q(X_{t},X^{\star}))$, then $X_{t+1}=X^{\star}$, otherwise $X_{t+1}=X_{t}$.

If $U\leq\min(1,\pi(Y^{\star})q(Y^{\star},Y_{t-1})/\pi(Y_{t-1})q(Y_{t-1},Y^{\star}))$, then $Y_{t}=Y^{\star}$, otherwise $Y_{t}=Y_{t-1}$.

Here we address the verification of Assumptions 2.1-2.3 for this algorithm. Assumption 2.1 can be verified for MH chains under conditions on the target and the proposal (Nummelin 2002; Roberts and Rosenthal 2004). In some settings the explicit drift function given in Theorem 3.2 of Roberts and Tweedie 1996b may be used to verify Assumption 2.2 as in Section 3.2. The probability of coupling at the next step given that the chains are in $X_{t}$ and $Y_{t-1}$ can be controlled as follows. First, the probability of proposing the same value $X^{\star}$ depends on the total variation distance between $q(X_{t},\cdot)$ and $q(Y_{t-1},\cdot)$, which is typically strictly positive if $X_{t}$ and $Y_{t-1}$ are in bounded subsets of $\mathcal{X}$. Furthermore, the probability of accepting $X^{\star}$ is often strictly positive on bounded subsets of $\mathcal{X}$, for instance when $\pi(x)>0$ for all $x\in\mathcal{X}$. Assumption 2.3 is satisfied by design thanks to the use of maximal couplings and common uniform variable $U$ in the above procedure.

Different considerations drive the choice of proposal distribution in standard MCMC and in our proposed estimators. In the case of random walk proposals with variance $\Sigma$, larger variances lead to smaller total variation distances between $q(X_{t},\cdot)$ and $q(Y_{t-1},\cdot)$ and thus larger probabilities of proposing identical values. However meeting events only occur if proposals are accepted, which is unlikely if $\Sigma$ is too large. This trade-off could lead to a different choice of $\Sigma$ than the optima known for the marginal chains (Roberts et al. 1997), and deserves further investigation.

We perform experiments with a $d$-dimensional Normal target distribution $\mathcal{N}(0,V)$, where $V$ is the inverse of a matrix drawn from a Wishart distribution with identity scale matrix and $d$ degrees of freedom. This setting, borrowed from Hoffman and Gelman 2014, yields Normal targets with strong correlations and a dense precision matrix. Below, each independent run is performed with an independent draw of $V$. We consider Normal random walk proposals with variance $\Sigma$ set to $V/d$. The division by $d$ heuristically follows from the scaling results of Roberts et al. 1997. We initialize the chains either from the target distribution, or from a Normal centered at $(1,\ldots,1)$ with identity covariance matrix. We first couple the proposals with a maximal coupling given by Algorithm 2. The resulting average meeting times, based on $1,000$ independent runs, are given in Figure 1a. The plot indicates an exponential increase of the average meeting times with the dimension, under both initialization strategies. In passing, this illustrates that meeting times can be large even if the chains marginally start at stationarity, i.e. in a setting where there is no burn-in bias.

Next we perform the same experiments with the reflection-maximum coupling described in the previous section. The results are shown in Figure 1b. The average meeting times now increase at a rate that appears closer to linear in the dimension. This is to be compared with established theoretical results on the linear performance of standard MH estimators with respect to the dimension (Roberts et al. 1997). A formal justification of the scaling observed in Figure 1b is an open question, and so is the design of more effective coupling strategies.

3 Gibbs sampling

Gibbs sampling is another popular class of MCMC algorithms, in which components of a Markov chain are updated alternately by sampling from the target conditional distributions (Robert and Casella 2004, Chapter 10 of), implemented e.g. in the software packages JAGS (Plummer et al. 2003). In Bayesian statistics, these conditional distributions sometimes belong to a standard family such as Normal, Gamma, or Inverse Gamma. Otherwise, the conditional updates might require MH steps. We can introduce couplings in each conditional update, using either maximal couplings of the target conditionals, if these are standard distributions, or maximal couplings of the proposal distributions in MH steps targeting the target conditionals. Controlling the probability of meeting at the next step over a set, as required for the application of Proposition 3.4, can be done on a case-by-case basis. Drift conditions for Gibbs samplers also tend to rely on case-by-case arguments (Rosenthal 1996, see e.g.).

Gibbs samplers tend to perform well for targets with weak correlations between the components being updated; otherwise Gibbs chains are expected to mix poorly. We perform numerical experiments on Normal target distributions in varying dimensions to observe the effect of correlations on the meeting times of coupled Gibbs chains. For each target $\mathcal{N}(0,V)$, we introduce an MH-within-Gibbs sampler, where each univariate component $i$ is updated with a single Metropolis step, using Normal proposals with variance $V_{i,i}$. Here an iteration of the sampler refers to a complete scan of the components. Figure 2a presents the median meeting times as a function of the dimension, when $V$ is the inverse of a Wishart draw as in the previous section. In this highly correlated setting, the meeting times scale poorly with the dimension. The plot presents the median instead of the average, because we have stopped the runs after $500,000$ iterations; the median is robust to this truncation, but not the average. We remark that shorter meeting times are obtained when initializing the chains away from the target distribution.

Next we consider a Normal target with covariance matrix $V$ defined by $V_{i,j}=0.5^{-|i-j|}$, which induces weak correlations among components; the inverse of $V$ is tridiagonal. In that case, the same Gibbs sampler performs much more favorably, as we can see from Figure 2b. The average meeting times seem to scale sub-linearly with the dimension, under both choices of initializations $\pi_{0}$. Couplings of other Gibbs samplers will be encountered in the numerical experiments of Section 5.

4 Coupling of other MCMC algorithms

Among extensions of the MH algorithm, Metropolis-adjusted Langevin algorithms (Roberts and Tweedie 1996a, e.g.) are characterized by the use of a proposal distribution given current state $X_{t}$ that is Normal with mean $X_{t}+h\nabla\log\pi(X_{t})/2$ and variance $h\Sigma$, with tuning parameter $h>0$ and covariance matrix $\Sigma$. Maximal couplings or reflection-maximal couplings of the proposals could be readily implemented to obtain faithful chains. Going further in the use of gradient information, Hamiltonian or Hybrid Monte Carlo (Duane et al. 1987; Neal 1993; Neal 2011, HMC,) is a popular MCMC algorithm for large-dimensional targets. In Heng and Jacob 2019, the framework of the present article is applied to pairs of Hamiltonian Monte Carlo chains, with a focus of the verification of Assumptions 2.1-2.3 in that context. Such couplings are analyzed in detail in Mangoubi and Smith 2017; Bou-Rabee et al. 2018 to obtain convergence rates for the underlying chains. We refer to Heng and Jacob 2019 for more details, and provide for completeness some experiments on the Normal target described above in the supplementary materials.

The present article generalizes unbiased estimators obtained by coupling conditional particle filters in Jacob et al. 2019. These algorithms, introduced in Andrieu et al. 2010, target the distribution of latent processes given observations and fixed parameters for nonlinear state space models. The couplings of conditional particle filters in Jacob et al. 2019 involve a combination of common random numbers and maximal couplings. Couplings of particle independent Metropolis–Hastings, which is a particular case of Metropolis–Hastings with an independent proposal distribution, are simpler to design and considered in Middleton et al. 2019.

The design of generic and efficient MCMC kernels is a topic of active ongoing research (see e.g. Murray et al. 2010; Goodman et al. 2010; Pollock et al. 2016; Vanetti et al. 2017; Titsias and Yau 2017, and references therein). Any new kernel could lead to unbiased estimators with the proposed framework, as long as appropriate couplings can be implemented.

Illustrations

Section 5.1 illustrates the impact of $k$, $m$, and the initial distribution $\pi_{0}$, identifying a situation where some care is required. Section 5.2 considers the removal of the bias from a Gibbs sampler previously considered for perfect sampling and regeneration methods. Section 5.3 introduces an Ising model and a coupling of a replica exchange algorithm, and we present experiments performed on parallel processors. Section 5.4 considers a high-dimensional variable selection example, with an MH algorithm previously shown to scale linearly with the number of variables. Finally, Section 5.5 focuses on the problem of approximating the cut distribution arising in modular inference, which illustrates the appeal of unbiased estimators beyond parallel computing.

We use a bimodal target distribution and a random walk MH algorithm to illustrate our method and highlight some of its limitations. In particular, we consider a mixture of univariate Normal distributions with density $\pi(x)=0.5\cdot\mathcal{N}(x;-4,1)+0.5\cdot\mathcal{N}(x;+4,1)$, which we sample from using random walk MH with Normal proposal distributions of variance $\sigma_{q}^{2}=9$. This enables regular jumps between the modes of $\pi$. We set the initial distribution $\pi_{0}$ to $\mathcal{N}(10,10^{2})$, so that chains are likely to start closer to the mode at $+4$ than the mode at $-4$. Over $1,000$ independent runs, we find that the meeting time $\tau$ has an average of $20$ and a $99\%$ quantile of $105$.

We present the results in Table 1. First, we see that the inefficiency is sensitive to the choice of $k$ and $m$. Second, we see that when $k$ and $m$ are sufficiently large we can retrieve an inefficiency comparable to that of the underlying MCMC algorithm. The ideal choice of $k$ and $m$ will depend on tradeoffs between inefficiency, the desired level of parallelism, and the number of processors available. We present a histogram of the target distribution, obtained using $k=200$, $m=2,000$, in Figure 3a. These histograms are produced by averaging unbiased estimators of expectations of indicator functions, corresponding to consecutive intervals. Confidence intervals at level $95\%$ are obtained from the central limit theorem and are represented as grey boxes, with vertical bars showing the point estimates.

Next, we consider a more challenging case by setting $\sigma_{q}^{2}=1$, again with $\pi_{0}=\mathcal{N}(10,10^{2})$. These values make it difficult for the chains to jump between the modes of $\pi$. Over $R=1,000$ runs we find an average meeting time of $769$, with a $99\%$ quantile of $9,186$. When the chains start in different modes, the meeting times are often dramatically larger than when the chains start by the same mode. One can still recover accurate estimates of the target distribution, but $k$ and $m$ have to be set to larger values. With $k=20,000$ and $m=30,000$, we obtain the $95\%$ confidence interval $[0.397,0.430]$ for $\int\mathds{1}(x>3)\pi(dx)\approx 0.421$. We show a histogram of $\pi$ in Figure 3b.

Finally we consider a third case, with $\sigma_{q}^{2}=1$ as before but now with $\pi_{0}$ set to $\mathcal{N}(10,1)$. This initialization makes it unlikely for a chain to start near the mode at $-4$. The pair of chains typically converge around the mode at $+4$ and meet in a small number of iterations. Over $R=1,000$ replications, we find an average meeting time of $9$ and a $99\%$ quantile of $35$. A $95\%$ confidence interval on ${\int\mathds{1}(x>3)\pi(dx)}$ obtained from the estimators with $k=50$, $m=500$ is $[0.799,0.816]$, far from the true value of $0.421$. The associated histogram of $\pi$ is shown in Figure 3c.

Sampling $9,000$ additional estimators yields a $95\%$ confidence interval $[-0.353,1.595]$, again using $k=50$, $m=500$. Among these extra $9,000$ values, a few correspond to cases where one chain jumped to the left-most mode before meeting the other. This resulted in large meeting times and thus a large empirical variance for $H_{k:m}$. Upon noticing a large empirical variance one can then decide to use larger values of $k$ and $m$. We conclude that although our estimators are unbiased and are consistent in the limit as $R\to\infty$, poor performance of the underlying Markov chains combined with ill-chosen initializations can still produce misleading results for any finite $R$, such as $1,000$ in this example.

2 Gibbs sampler for nuclear pump failure data

Next we consider a classic Gibbs sampler for a model of pump failure counts, used e.g. in Murdoch and Green 1998 to illustrate perfect samplers for continuous distributions, and in Mykland et al. 1995 to illustrate their regeneration approach. Here we focus on a comparison with the regeneration approach, which was motivated by similar practical concerns as this paper, in particular to avoid an arbitrary choice of burn-in, construct confidence intervals on the expectations of interest, and make principled use of parallel processors. In that paper the authors show how to construct regeneration times – random times between which the chain forms independent and identically distributed “tours”. The authors define a consistent estimator for arbitrary test functions, whose asymptotic variance takes a simple form. The estimator is then obtained by aggregating over these independent tours.

The data consist of operating times $(t_{n})_{n=1}^{K}$ and failure counts $(s_{n})_{n=1}^{K}$ for $K=10$ pumps at the Farley-1 nuclear power station, as first described in Gaver and O’Muircheartaigh 1987. The model specifies $s_{n}\sim\text{Poisson}(\lambda_{n}t_{n})$ and $\lambda_{n}\sim\text{Gamma}(\alpha,\beta)$, where $\alpha=1.802$, $\beta\sim\text{Gamma}\left(\gamma,\delta\right)$, $\gamma=0.01$, and $\delta=1$. The Gibbs sampler for this model consists of the following update steps:

Here $\text{Gamma}(\alpha,\beta)$ refers to the distribution with density $x\mapsto\Gamma(\alpha)^{-1}\beta^{\alpha}x^{\alpha-1}\exp(-\beta x)$. We initialize all parameter values to 1 (the initialization is not specified in Mykland et al. 1995). To form our estimator we apply maximal couplings at each conditional update of the Gibbs sampler, as described in Section 4.3.

We begin by drawing $1,000$ meeting times independently. Following the guidelines of Section 3.1, we set $k=7$, corresponding to the $99\%$ quantile of $\tau$ and $m=10\cdot k=70$. For the regeneration approach, Mykland et al. 1995 gives a set of tuning parameters which we adopt below. Applying the regeneration approach to 1,000 Gibbs sampler runs of 5,000 iterations each, we observe on average 1,996 complete tours per run with an average length of 2.50 iterations per tour. These values agree with the count of 1,967 tours of average length 2.56 reported in Mykland et al. 1995. We observe a posterior mean estimate for $\beta$ of 2.47 with a variance of $1.89\times 10^{-4}$ over the 1,000 independent runs, which implies an efficiency value of $(5,000\cdot 1.89\times 10^{-4})^{-1}=1.06$. This exceeds the efficiency of $0.94$ achieved by our estimator with the choice of $k=7$ and $m=70$. On the other hand, the regeneration approach often requires more extensive analytical work with the underlying Markov chain; we refer to Mykland et al. 1995 for a detailed description. For reference, the underlying Gibbs sampler achieves an efficiency of $1.08$, based on a long run of $5\times 10^{5}$ iterations and a burn-in of $10^{3}$ iterations. More extensive comparisons with other regeneration approaches such as that of Brockwell and Kadane 2005 would deserve investigation.

3 Ising model

We consider an Ising model on a $32\times 32$ square lattice with periodic boundaries. This provides a setting where a basic MCMC sampler can mix slowly depending on an inverse temperature parameter $\theta$, and where a replica exchange strategy as in Geyer 1991 can be helpful. We also use this example to illustrate the use of our estimators on a large computing cluster, with the considerations reviewed in Section 3.3. For $i$ and $j$ in $\{1,\ldots,32\}^{2}$ we write $i\sim j$ if $i$ and $j$ are neighbors in the square lattice with periodic boundaries. We write $x_{i}\in\{-1,+1\}$ for the spin at location $i$, and $x=\{x_{i}\}$ for the full grid. We write $t(x)$ for the “natural statistic” ${t(x)=0.5\sum_{i\in\{1,\ldots,32\}^{2}}\sum_{j\sim i}x_{i}x_{j}}$ summing the products of pairs of neighbors. The $0.5$ multiplier here results in each pair of neighboring sites only being counted once. Under the model, the probability associated with a grid $x$ is $\pi_{\theta}(x)\propto\exp(\theta t(x))$, where $\theta>0$ denotes an inverse temperature parameter that calibrates the degree of correlation between neighboring sites.

We consider a single-site Gibbs sampler, called a heat bath algorithm in this context, to approximate the distribution $\pi_{\theta}$ given a value of $\theta$. One iteration of the algorithm consists of a sweep through all the locations $i\in\{1,\ldots,32\}^{2}$. For each $i$ we draw $x_{i}$ from its conditional distribution under $\pi_{\theta}$ given all the other spins. It can be checked that the conditional probability of $\{x_{i}=+1\}$ given the other spins equals $\exp(\theta s_{i})/(\exp(\theta s_{i})+\exp(-\theta s_{i}))$, where $s_{i}$ denotes the sum of spins over the four neighbors of $i$. We initialize the chains by drawing spins uniformly in $\{-1,+1\}$ at each site, independently across sites.

A simple strategy to couple heat bath chains consists of sampling from the maximal coupling of each conditional distribution. For a grid of $\theta$ values from $0.3$ and $0.48$, we run 100 pairs of chains until they meet. We then plot the average meeting time as a function of $\theta$ in Figure 4a, noting that the average meeting time increases sharply to values above $10^{6}$ as $\theta$ approaches its critical value (see the related discussion in Propp and Wilson 1996). We conclude that it would be expensive to produce unbiased estimators based on the heat bath algorithm for values of $\theta$ above $0.48$, for reasons related to the behavior of the underlying algorithm.

There are several ways to address the degeneracy of the heat bath algorithm as $\theta$ increases. Specialized algorithms have been proposed to jointly update groups of spins (Swendsen and Wang 1987; Wolff 1989). Here, we consider an approach based on an ensemble of $N$ chains that regularly exchange their states, a technique often termed replica exchange or parallel tempering. Following e.g. Geyer 1991, we introduce $N$ chains, $x^{(1)}$, …, $x^{(N)}$, with each $x^{(n)}$ targeting $\pi_{\theta^{(n)}}$ with different values of $\theta^{(n)}$ ordered as $\theta^{(1)}<\ldots<\theta^{(N)}$. Each iteration of the algorithm proceeds as follows. With probability $p_{\text{swap}}\in(0,1)$, for $n\in\{1,\ldots,N-1\}$ (sequentially), we propose exchanging the states $x^{(n)}$ and $x^{(n+1)}$ corresponding to $\theta^{(n)}$ and $\theta^{(n+1)}$. We accept this swap with probability $\min(1,\pi_{\theta^{(n)}}(x^{(n+1)})\pi_{\theta^{(n+1)}}(x^{(n)})/(\pi_{\theta^{(n)}}(x^{(n)})\pi_{\theta^{(n+1)}}(x^{(n+1)})))$, which simplifies to ${\min(1,\exp((\theta^{(n)}-\theta^{(n+1)})(t(x^{(n+1)})-t(x^{(n)}))))}$. Otherwise we perform a full sweep of single-site Gibbs updates, independently across chains.

A coupling of this algorithm involves a pair of ensembles with $N$ chains each; the two ensembles are identical if chain $n$ in the first ensemble equals chain $n$ in the second ensemble, for all $n\in\{1,\ldots,N\}$. We use common random numbers to decide whether to perform swap moves or single-site Gibbs moves, and whether to accept the proposed states in the event of a swap move. In the event of a single-site Gibbs move, we maximally couple each conditional update.

Throughout the following experiments we use $p_{\text{swap}}=0.01$, and introduce an equally spaced grid of $\theta$ values from $\theta^{(1)}=0.3$ to $\theta^{(N)}=0.55$ for several different choices of $N$. We note that these grids includes $\theta$ values at which we have seen that the single-site Gibbs sampler mixes poorly. Figure 4b shows the resulting average meeting times over 100 independent runs, as a function of the number of chains $N$. The average meeting time first decreases with the number of chains, but then increases again. A possible explanation is that the mixing of the chains first improves as $N$ increases, and then stabilizes; on the other hand it becomes harder for the ensembles to meet when $N$ increases since all chains in the ensembles have to meet. The minimum average meeting time is here attained for $N=16$ chains per ensemble.

4 Variable selection

We now consider an experiment like those of Yang et al. 2016. We define

and generate $Y$ given $X$ and $\beta^{\star}$ from the model with $\sigma^{2}=1$, $\sigma_{0}^{2}=1$, $n\in\{500,1000\}$, $p\in\{1000,5000\}$, and signal-to-noise parameter $\text{SNR}\in\{0.5,1,2\}$. We also set $s_{0}=100$, $g=p^{3}$, and $\kappa=2$ (exactly as in Yang et al. 2016; the value of $\kappa$ was obtained by personal communication) and generate the covariates $X$ using a multivariate normal distribution with covariance matrix $\Sigma$ either equal to a unit diagonal matrix or with entries $\Sigma_{ij}=\exp(-|i-j|)$. We refer to these two cases as the independent design and correlated design cases, respectively. We draw from the initial distribution $\pi_{0}$ by creating a vector of $p$ zeros, sampling $s_{0}$ coordinates uniformly from $\{1,\ldots,p\}$ without replacement, and setting the corresponding entries to $1$ with probability $0.5$.

For different values of $n$, $p$ and SNR, and the two types of design, we run coupled chains 100 times independently until they meet. We report the average meeting times in Tables 2 and 3. The average meeting times are of the order of $10^{4}$ to $10^{5}$, depending on the problem; the maximum is attained in the correlated design at $n=500,p=1000,\text{SNR}=2$. In contrast with this, the experiments in Yang et al. 2016 identify the scenario $n=500,p=5000,\text{SNR}=1$ as the most challenging one. This discrepancy deserves further study; it could be due to variations from a synthetic data set to another, or to differences in the criteria being reported.

To illustrate the impact of dimension, we focus on the independent design setting with $n=500$ and $\text{SNR}=1$, and consider values of $p$ between $100$ and $1000$. For each value of $p$, we run coupled chains $1,000$ times independently until they meet. We present violin plots representing the distributions of meeting times divided by $p$ in Figure 6a. The distribution of scaled meeting times appears to be approximately constant as a function of $p$, suggesting that meeting times increase linearly in $p$. This is consistent with the findings of Yang et al. 2016, where mixing times are shown to increase linearly in $p$.

Figure 6b shows the results in the form of $95\%$ confidence intervals shown as error bars, using (3.4), the CLT relevant when the time budget is fixed and the number of processors grows large. We observe that $\kappa$ has a strong impact on the probability of including the first 10 variables in this setting, and that the most satisfactory results are obtained for $\kappa=0.1$ rather than for $\kappa=2$, recalling that $\beta^{\star}$ has non-zero entries in its first 10 components. Note that the error bars are narrow but still noticeable, particularly for $\kappa=0.1$. On the same figure, the solid lines represent estimates obtained with 10 independent MCMC runs with $10^{6}$ iterations each, discarding the first $10^{5}$ iterations as burn-in. These MCMC estimates present noticeable variability in spite of the large number of iterations. In a standard MCMC setting, we might run chains for more iterations until the estimates agree across independent runs. In the proposed framework, we increase the precision by generating more independent unbiased estimators without necessarily modifying $k$ or $m$.

Figure 6b suggests that the variable selection procedure considered here is sensitive to the prior hyperparameter $\kappa$; we refer to Yang et al. 2016, and to Johnson 2013; Nikooienejad et al. 2016 for related discussions on Bayesian variable selection in high dimension and convergence of MCMC.

5 Cut distribution

Finally, our proposed estimator can be used to approximate the cut distribution, which poses a significant challenge for existing MCMC methods (Plummer 2014; Jacob et al. 2017). This illustrates another appeal of the unbiasedness property, beyond the motivation for parallel computation.

Consider two models, one with parameters $\theta_{1}$ and data $Y_{1}$ and another with parameters $\theta_{2}$ and data $Y_{2}$, where the likelihood of $Y_{2}$ might depend on both $\theta_{1}$ and $\theta_{2}$. For instance the first model could be a regression with data $Y_{1}$ and coefficients $\theta_{1}$, and the second model could be another regression whose covariates are the residuals, coefficients, or fitted values of the first regression (Pagan 1984; Murphy and Topel 2002). In principle one could introduce an encompassing model and conduct joint inference on $\theta_{1}$ and $\theta_{2}$ via the posterior distribution. In that case, misspecification of either model would lead to misspecification of the ensemble and thus to a misleading quantification of uncertainty, as noted in several studies (Liu et al. 2009; Plummer 2014; Lunn et al. 2009; McCandless et al. 2010; Zigler 2016; Blangiardo et al. 2011, e.g.).

The cut distribution (Spiegelhalter et al. 2003; Plummer 2014) allows the propagation of uncertainty about $\theta_{1}$ to inference on $\theta_{2}$ while preventing misspecification in the second model from affecting estimation in the first. The cut distribution is defined as

Here $\pi_{1}(\theta_{1})$ refers to the distribution of $\theta_{1}$ given $Y_{1}$ in the first model alone, and $\pi_{2}(\theta_{2}|\theta_{1})$ refers to the distribution of $\theta_{2}$ given $Y_{2}$ and $\theta_{1}$ in the second model. Often, the density $\pi_{2}(\theta_{2}|\theta_{1})$ can only evaluated up to a constant in $\theta_{2}$, which may vary with $\theta_{1}$. This makes the cut distribution difficult to approximate with MCMC algorithms (Plummer 2014).

A naive approach consists of first running an MCMC algorithm targeting $\pi_{1}(\theta_{1})$ to obtain a sample $(\theta_{1}^{n})_{n=1}^{N_{1}}$, perhaps after discarding a burn-in period and thinning the chain. Then for each $\theta_{1}^{n}$, one can run an MCMC algorithm targeting $\pi_{2}(\theta_{2}|\theta_{1}^{n})$, yielding $N_{2}$ samples. One might again discard some burn-in and thin the chains, or just keep the final state of each chain. The resulting joint samples approximate the cut distribution. However, the validity of this approach relies on a double limit in $N_{1}$ and $N_{2}$. Diagnosing convergence may also be difficult given the number of chains in the second stage, each of which targets a different distribution $\pi_{2}(\theta_{2}|\theta_{1}^{n})$.

We consider the example described in Plummer 2014, inspired by an investigation of the international correlation between human papillomavirus (HPV) prevalence and cervical cancer incidence (Maucort-Boulch et al. 2008). The first module concerns HPV prevalence, with data independently collected in $13$ countries. The parameter $\theta_{1}=(\theta_{1,1},\ldots,\theta_{1,13})$ receives a Beta$(1,1)$ prior distribution independently for each component. The data $(Y_{1},\ldots,Y_{13})$ consist of $13$ pairs of integers. The first represents the number of women infected with high-risk HPV, and the second represents population sizes. The likelihood specifies a Binomial model for $Y_{i}$, independently for each component $i$. The posterior for this model is given by a product of Beta distributions.

where the data $(Z_{1,i},Z_{2,i})_{i=1}^{13}$ are pairs of integers. The first component represents numbers of cancer cases, while the second is a number of woman-years of follow-up. The Poisson regression model might be misspecified, motivating departures from inference based on the joint model (Plummer 2014).

Here we can draw directly from the first posterior, denoted by $\pi_{1}(\theta_{1})$, and obtain a sample $(\theta_{1}^{n})_{n=1}^{N}$. For each $\theta_{1}^{n}$ we consider an MH algorithm targeting $\pi_{2}(\theta_{2}|\theta_{1}^{n})$, using a Normal random walk proposal with variance $\Sigma$. We couple this algorithm using reflection-maximal couplings of the proposals as in Section 4.1. In preliminary runs, starting with a standard bivariate Normal as an initial distribution and a proposal covariance matrix set to identity, we estimate the first two moments of the cut distribution, and we use them to refine the initial distribution $\pi_{0}$ and the proposal covariance matrix $\Sigma$. With these settings we obtain a distribution of meeting times shown in Figure 7a. We then set $k=100$, $m=10k$, and obtain approximations of the cut distribution represented by histograms in Figures 7b and 7c, using $N=10,000$ unbiased estimators. The overlaid curves correspond to a kernel density estimate obtained by running $m=1,000$ steps of MCMC targeting $\pi_{2}(\theta_{2}|\theta_{1}^{n})$ with $\theta_{1}^{n}$ drawn from $\pi_{1}(\theta_{1})$, for $n\in\{1,\ldots,N\}$, and keeping the final $m$-th state of each chain. The proposed estimators can be refined by increasing the number $N$ of independent replications, whereas the MCMC estimators would converge only in the double limit of $N$ and $m$ going to infinity.

Discussion

By combining the powerful technique of Glynn and Rhee 2014 with couplings of MCMC algorithms, unbiased estimators of integrals with respect to the target distribution can be constructed. Their efficiency can be controlled with tuning parameters $k$ and $m$, for which we have proposed guidelines: $k$ can be chosen as a large quantile of the meeting time $\tau$, and $m$ as a multiple of $k$. Improving on these simple guidelines stands as a subject for future research. In numerical experiments we have argued that the proposed estimators yield a practical way of parallelizing MCMC computations in a range of settings. We stress that coupling pairs of Markov chains does not improve their marginal mixing properties, and that poor mixing of the underlying chains can lead to poor performance of the resulting estimator. The choice of initial distribution $\pi_{0}$ can have undesirable effects on the estimators, as in the multimodal example of Section 5.1. Unreliable estimators would also result from stopping the chains before their meeting time.

Couplings of MCMC algorithms can be devised using maximal couplings, reflection couplings, and common random numbers. We have focused on couplings that can be implemented without further analytical knowledge about the target distribution or about the MCMC kernels. However, these couplings might result in prohibitively large meeting times, either because the marginal chains mix slowly, as in 5.1, or because the coupling strategy is ineffective, as in Section 4.2.

Regarding convergence diagnostics, the proposed framework yields the following representation for the total variation between $\pi_{k}$ and $\pi$, where $\pi_{k}$ denotes the marginal distribution of $X_{k}$:

Thanks to its potential for parallelization, the proposed framework can facilitate consideration of MCMC kernels that might be too expensive for serial implementation. For instance, one can improve MH-within-Gibbs samplers by performing more MH steps per component update, HMC by using smaller step-sizes in the numerical integrator (Heng and Jacob 2019), and particle MCMC by using more particles in the particle filters (Andrieu et al. 2010; Jacob et al. 2019). We expect the optimal tuning of MCMC kernels to be different in the proposed framework than when used marginally.

On top of enabling the application of the results of Glynn and Heidelberger 1991 to accomodate budget constraints, the lack of bias of the proposed estimators can be beneficial in combination with the law of total expectation, to implement modular inference procedures as in Section 5.5. In Rischard et al. 2018 the lack of bias is exploited in new estimators of Bayesian cross-validation criteria. In Chen et al. 2018 similar unbiased estimators are used in the expectation step of an expectation-maximization algorithm. There may be other settings where the lack of bias is appealing, for instance in gradient estimation for stochastic gradient descents (Tadić et al. 2017).

The authors are grateful to Jeremy Heng and Luc Vincent-Genod for useful discussions. The authors gratefully acknowledge support by the National Science Foundation through grants DMS-1712872 and DMS-1844695 (Pierre E. Jacob), and DMS-1513040 (Yves F. Atchadé).

References

Appendix A Proofs

A.2 Proof of Proposition 3.2

The above implies that for any finite $i\geq k$ we have

almost surely by the strong law of large numbers. Assumption 2.2 implies that this quantity goes to 0 as $i\to\infty$.

A.3 Proof of Proposition 3.3

for some arbitrary bounded sequence $(b_{t})_{t\geq 0}$. Fix an integer $N\geq k$, and set

The same argument as in the proof of Proposition 3.1 can be applied here and shows that $S_{k}^{(N)}$ is a Cauchy sequence in $L_{2}$ that converges to $S_{k}$, as $N\to\infty$, so that

where for a function $W:\;\mathsf{X}\to[1,\infty)$, the $W$-norm between two probability measures $\mu,\nu$ is defined as

and $|f|_{W}:=\sup_{x}|f(x)|/W(x)$. This result can be found in Theorem 15.0.1 of Meyn and Tweedie 2009. It follows from (A.1) that $\sum_{j\geq 0}|P^{j}(h-\pi(h))(x)|\leq|h|_{V^{\beta}}\sum_{j\geq 0}\|P^{j}(x,\cdot)-\pi\|_{V^{\beta}}\leq\frac{C_{\beta}|h|_{V^{\beta}}V^{\beta}(x)}{1-\rho_{\beta}}<\infty$. Hence the function

is well-defined and measurable (as a limit of a sequence of measurable functions) and satisfies $|g(x)|\leq C_{\beta}|h|_{V^{\beta}}V^{\beta}(x)/(1-\rho_{\beta})$. And since $PV^{\beta}$ is finite everywhere, by Lebesgue’s dominated convergence we deduce that $Pg$ is finite everywhere as well and

Hence, with $Z_{t}:=(X_{t},Y_{t-1})$, and $\bar{g}(x,y)=g(x)-g(y)$, we have

Using this and a telescoping sum argument, we write

Since $\bar{g}(Z_{t+1})=0$ on $\{\tau=t+1\}$, the last term in the above display reduces to $\sum_{t=k}^{N-1}(b_{t+1}-b_{t})\bar{g}(Z_{t+1})\mathds{1}(\tau>t+1)$, and we obtain

Let $\mathcal{F}_{t}$ denote the sigma-algebra generated by the variables $X_{0},(X_{1},Y_{0}),\ldots,(X_{t},Y_{t-1})$. Note that $\{\tau>t\}$ belongs to $\mathcal{F}_{t}$. Hence

In other words, $\left\{\left(\sum_{t=k}^{j}b_{t}\left(\bar{g}(Z_{t+1})-\bar{P}[\bar{g}](Z_{t})\right)\mathds{1}(\tau>t),\mathcal{F}_{j}\right),\;k\leq j\leq N-1\right\}$ is a martingale. The orthogonality of the martingale increments gives

We use this together with (A.2), the convexity of the squared norm, and Minkowski’s inequality to conclude that

Assumption 3.1 together with $\pi_{0}(V)<\infty$, implies that

as seen above. In conclusion, all the expectations appearing in (A.3) are upper bounded by some constant times terms of the form $\delta_{\beta}^{t}$. We conclude that

In the particular case of $\eta_{k:m}$, we have $\eta_{k:m}=\mathds{1}(\tau>k)\sum_{t=k}^{\tau-1}\min\left(1,\frac{t+1-k}{m+1-k}\right)\left(h(X_{t+1})-h(Y_{t})\right)$. Hence $b_{k}=0$, $b_{t}=(t-k)/(m-k+1)$ if $k<t\leq m+1$, $b_{t}=1$ if $t>m+1$. We then obtain the bound of Proposition 3.3.

A.4 Proof of Proposition 3.4

Here $Z_{n}$ is defined as $(X_{n},Y_{n-1})$ for all $n\geq 1$. The assumption in (3.3), within the statement of Proposition 3.4, implies that for $(x,y)\in\mathcal{C}\times\mathcal{C}$, $\bar{P}$ can be written as a mixture

where $\epsilon_{x,y}\geq\epsilon$, $\nu_{x,y}(dz)$ is a restriction of $\bar{P}\left((x,y),dz\right)$ on $\mathcal{D}$ (that is for any measurable subset $A$ of $\mathcal{D}$, $\bar{P}\left((x,y),A\right)=P\left((x,y),A\cap\mathcal{D}\right)/P\left((x,y),\mathcal{D}\right)$), and $R\left((x,y),dz\right)$ is the restriction of $\bar{P}\left((x,y),dz\right)$ on $(\mathcal{X}\times\mathcal{X})\setminus\mathcal{D}$. This means that whenever $(x,y)\in\mathcal{C}\times\mathcal{C}$ one can sample from $\bar{P}((x,y),\cdot)$ by drawing independently a Bernoulli random variable $J$, with probability of success $\epsilon_{x,y}$. Then if $J=1$, we draw from $\nu_{x,y}$, if $J=0$, we draw from $R\left((x,y),\cdot\right)$. From this decomposition, the proof of the proposition follows the same lines as in Douc et al. 2004, and we give the details only for completeness. We cannot directly invoke their result since their assumptions do not seem to apply to our setting.

Set $\bar{V}(x,y)=\frac{1}{2}(V(x)+V(y))$. First we show that the bivariate kernel satisfies a geometric drift towards $\mathcal{C}\times\mathcal{C}$. That is, there exists $\alpha\in(0,1)$ such that

Indeed for $(x,y)\notin\mathcal{C}\times\mathcal{C}$, since $V\geq 1$, and $\mathcal{C}=\{V\leq L\}$, $\bar{V}(x,y)\geq(1+L)/2$. In other words, $\frac{1}{2}\leq\bar{V}(x,y)/(1+L)$. Therefore,

with $\alpha=\lambda+\frac{b}{1+L}<1$. We set

In this section $\mathds{1}_{\mathcal{S}}(\cdot)$ refers to the indicator function on the set $\mathcal{S}$. Let $N_{n}$ denote the number of visits to $\mathcal{C}\times\mathcal{C}$ by time $n$. Then

The event $\{\tau>n,\;N_{n-1}\geq j\}$ implies that no success occurred within at least $j$ independent Bernoulli random variables each with probability of success at least $\epsilon$. Hence

For the second term, we have (since $B\geq 1$, and the chains stay together after meeting via Assumption 2.3),

Then use Markov’s inequality to conclude that

Since $\alpha<1$, there exists an integer $k_{0}\geq 1$ such that $\alpha B^{\frac{1}{k_{0}}}<1$. In that case for $n\geq k_{0}$ one can take $j=\lceil n/k_{0}\rceil$, to get

Suppose now that $Z_{n}\in\mathcal{C}\times\mathcal{C}$. Then $N_{n}=N_{n-1}+1$. Hence

Appendix B Hamiltonian Monte Carlo on multivariate Normals

As Section 4 in the main document, we perform experiments with a $d$-dimensional Normal target distribution $\pi=\mathcal{N}(0,V)$, where $V$ is the inverse of a matrix drawn from a Wishart distribution, with identity scale matrix and $d$ degrees of freedom. We provide average meeting times obtained in varying dimensions, when using the coupled HMC algorithm described in Heng and Jacob 2019. The latter article presents similar experiments but for a different choice of matrix $V$, thus we provide the present section for completeness.

Let us introduce a Markov kernel $P$ as a mixture of two kernels, an MH kernel $P_{\text{MH}}$ and an HMC kernel $P_{\text{HMC}}$. We first describe $P_{\text{MH}}$ and a coupling of it. The kernel is an MH kernel with Normal random walk proposals, with a covariance matrix equal to $10^{-8}$ times the identity matrix. The coupled version of $P_{\text{MH}}$ uses a maximal coupling of the proposals (as in Algorithm 2 of the main document).

The kernel $P_{\text{HMC}}$ corresponds to an HMC algorithm, with mass matrix given by the inverse of the target variance $V$. This preconditioning mechanism is motivated by considerations similar to those in Girolami and Calderhead 2011. In the present case of Normal distributions, this is particularly advantageous as it leads to a complete decoupling of the $d$ components of the target; see Proposition 3.1 in Bou-Rabee and Sanz-Serna 2018. We discretize Hamiltonian equations with a leap-frog integrator, using a stepsize of $\varepsilon=0.1\times d^{-1/4}$, and a number of steps of $L=1+\lfloor\varepsilon^{-1}\rfloor$, which corresponds to a trajectory length $\varepsilon L$ of approximately one. The coupling of such Hamiltonian kernels is done by using common random numbers for the momentum variables, i.e. a synchronous coupling [Bou-Rabee et al. 2018]. The kernels $P_{\text{MH}}$ and $P_{\text{HMC}}$, and their coupled counterparts, are combined into mixtures $P$ and $\bar{P}$, by assigning weights respectively of $0.05$ and $0.95$, i.e. an MH step is performed with probability $0.05$.

We consider two types of initialization $\pi_{0}$: either the target distribution $\pi$, or a Normal distribution $\mathcal{N}(1_{d},I_{d})$, with $1_{d}$ a vector of ones and $I_{d}$ the identity matrix. With the latter initialization, we observed a very low acceptance rate when using the stepsize $\varepsilon$ given above. We did not observe any issue when the chains were started from $\pi$. We also did not observe such issues when $V$ was replaced by the identity matrix. In principle, smaller stepsizes could be chosen, in order to increase the acceptance rate. However, this would result in more expensive iterations as $L$ is defined as $1+\lfloor\varepsilon^{-1}\rfloor$ above. Instead, we resort to the following heuristic strategy, which appears to solve the issue in the present example. We draw an initial position from $\mathcal{N}(1_{d},I_{d})$, then we perform $10$ steps of “unadjusted HMC”, with $\varepsilon=0.1\times d^{-1/4}$, and $L=1+\lfloor\varepsilon^{-1}\rfloor$ as described above. By “unadjusted HMC”, we refer to a scheme where the final point of the Hamiltonian trajectory is accepted with probability one, i.e. no MH correction is applied. This initialization procedure is simply a redefinition of $\pi_{0}$, and thus would not jeopardize the validity of the proposed estimators.

The results under both initializations are shown in Figure 8, where we observe average meeting times that increase very slowly with the dimension of the target distribution. Thanks to the initialization strategy described above, we obtain similar meeting times when starting from the chains from $\pi$ or from $\mathcal{N}(1_{d},I_{d})$.

Appendix C Baseball batting averages

We consider a classic Gibbs sampler discussed in the context of parallel computing in Rosenthal 2000. In Rosenthal 1996, it was proved that the chain produced by this Gibbs sampler converges in total variation to within $1\%$ of its stationary distribution after at most $140$ iterations. This type of derivation is technically challenging and has been done only in specific cases. Using this result, Rosenthal 2000 recommends to run parallel chains with a burn-in of $140$ iterations. We will compare this value with choices of $k$ and $m$ for the proposed unbiased estimators.

The data $(Z_{n})_{n=1}^{K}$ are baseball players’ batting averages taken from Table 1 of Morris 1983, where $K=18$. In the model, each $Z_{n}$ is assumed to follow $\mathcal{N}(\theta_{n},V)$ where $V$ is fixed to $0.00434$. Then, $\theta_{n}$ is assumed to follow $\mathcal{N}(\mu,A)$, where $\mu$ is given a flat prior and $A$ an Inverse Gamma $(a,b)$, with $a=-1$ and $b=2$. Here the Inverse Gamma $(\alpha,\beta)$ distribution has pdf $p(x;\alpha,\beta)=\Gamma(\alpha)^{-1}\beta^{\alpha}x^{-\alpha-1}\exp(-\beta/x)$. The Gibbs updates are as follows: