Unbiased Markov chain Monte Carlo with couplings

Given a target probability distribution $\pi$ on a Polish space $\mathcal{X}$ and a measurable real-valued test function $h$ integrable with respect to $\pi$, we want to estimate the expectation $\mathbb{E}_{\pi}[h(X)]=\int h(x)\pi(dx)$. Let $P$ denote a Markov transition kernel on $\mathcal{X}$ that leaves $\pi$ invariant, and let $\pi_{0}$ be some initial probability distribution on $\mathcal{X}$. Our estimators are based on a coupled pair of Markov chains $(X_{t})_{t\geq 0}$ and $(Y_{t})_{t\geq 0}$, which marginally start from $\pi_{0}$ and evolve according to $P$. In particular, we suppose that $\bar{P}$ is a transition kernel on the joint space $\mathcal{X}\times\mathcal{X}$ such that $\bar{P}((x,y),A\times\mathcal{X})=P(x,A)$ and $\bar{P}((x,y),\mathcal{X}\times A)=P(y,A)$ for any $x,y\in\mathcal{X}$ and any measurable set $A$. We then construct the coupled Markov chain $(X_{t},Y_{t})_{t\geq 0}$ as follows. We draw $(X_{0},Y_{0})$ such that $X_{0}\sim\pi_{0}$ and $Y_{0}\sim\pi_{0}$. Given $(X_{0},Y_{0})$ we draw $X_{1}\sim P(X_{0},\cdot)$, and then for any $t\geq 1$, given $X_{0},(X_{1},Y_{0}),\ldots,(X_{t},Y_{t-1})$, we draw $(X_{t+1},Y_{t})\sim\bar{P}((X_{t},Y_{t-1}),\cdot)$. We consider the following assumptions.

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).

Under these assumptions we introduce the following motivation for an unbiased estimator of $\mathbb{E}_{\pi}[h(X)]$, following Glynn and Rhee (2014). We begin by writing $\mathbb{E}_{\pi}[h(X)]$ as $\lim_{t\to\infty}\mathbb{E}[h(X_{t})]$. Then for any fixed $k\geq 0$,

	$\displaystyle\mathbb{E}_{\pi}[h(X)]$	$\displaystyle=\mathbb{E}[h(X_{k})]+\sum_{t=k+1}^{\infty}(\mathbb{E}[h(X_{t})]-\mathbb{E}[h(X_{t-1})])$	expanding the limit as a telescoping sum,
		$\displaystyle=\mathbb{E}[h(X_{k})]+\sum_{t=k+1}^{\infty}(\mathbb{E}[h(X_{t})]-\mathbb{E}[h(Y_{t-1})])$	since the chains have the same marginals,
		$\displaystyle=\mathbb{E}[h(X_{k})+\sum_{t=k+1}^{\infty}(h(X_{t})-h(Y_{t-1}))\text{]}$	swapping the expectations and limit,
		$\displaystyle=\mathbb{E}[h(X_{k})+\sum_{t=k+1}^{\tau-1}(h(X_{t})-h(Y_{t-1}))]$	by Assumption 2.3.

We note that the sum in the last equation is zero if $k+1>\tau-1$. The heuristic argument above suggests that the estimator ${H_{k}(X,Y)=h(X_{k})+\sum_{t=k+1}^{\tau-1}(h(X_{t})-h(Y_{t-1}))}$ should have expectation $\mathbb{E}_{\pi}[h(X)]$. We observe that this estimator requires $\tau-1$ calls to $\bar{P}$ and $\max(1,k+1-\tau)$ calls to $P$; thus under Assumption 2.2 its cost has a finite expectation.

In Section 3 we establish the validity of the estimator under the three conditions above; this formally justifies the swap of expectation and limit. The estimator can be viewed as a debiased version of $h(X_{k})$. Unbiasedness is guaranteed for any choice of $k\geq 0$, but both the cost and variance of $H_{k}(X,Y)$ are sensitive to $k$; see Section 3.1. Thanks to this unbiasedness property, we can sample $R\in\mathbb{N}$ independent copies of $H_{k}(X,Y)$ in parallel and average the results to estimate $\mathbb{E}_{\pi}[h(X)]$ consistently as $R\to\infty$; we defer further considerations on the use of unbiased estimators on parallel processors to Section 3.3.

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, 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 (see also Diaconis and Freedman, 1999), 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 (see discussion and references in Cowles and Rosenthal, 1998). 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, 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).

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 $[0,1]$, 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.

To motivate our next estimator, we note that we can compute $H_{k}(X,Y)$ for several values of $k$ from the same realization of the coupled chains, and that the average of these is unbiased as well. For any fixed integer $m$ with $m\geq k$, we can run coupled chains for $\max(m,\tau)$ iterations, compute the estimator $H_{\ell}(X,Y)$ for each $\ell\in\{k,\ldots,m\}$, and take the average ${H_{k:m}(X,Y)=(m-k+1)^{-1}\sum_{\ell=k}^{m}H_{\ell}(X,Y)}$, as we summarize in Algorithm 1. We refer to $H_{k:m}(X,Y)$ as the time-averaged estimator; the estimator $H_{k}(X,Y)$ is retrieved when $m=k$. Alternatively we could average the estimators $H_{\ell}(X,Y)$ using weights $w_{\ell}\in\mathbb{R}$ for $\ell\in\{k,\ldots,m\}$, to obtain $\sum_{\ell=k}^{m}w_{\ell}H_{\ell}(X,Y)$. This will be unbiased if $\sum_{\ell=k}^{m}w_{\ell}=1$.

Rearranging terms in $(m-k+1)^{-1}\sum_{\ell=k}^{m}H_{\ell}(X,Y)$, we can write the time-averaged estimator as

$\displaystyle H_{k:m}(X,Y)$

$\displaystyle=\frac{1}{m-k+1}\sum_{\ell=k}^{m}h(X_{\ell})+\sum_{\ell=k+1}^{\tau-1}\min\left(1,\frac{\ell-k}{m-k+1}\right)(h(X_{\ell})-h(Y_{\ell-1})).$

(2.1)

The term $(m-k+1)^{-1}\sum_{\ell=k}^{m}h(X_{\ell})$ corresponds to a standard MCMC average with $m$ total iterations and a burn-in period of $k-1$ iterations. We can interpret the other term as a bias correction. If $\tau\leq k+1$, then the correction term equals zero. This provides some intuition for the choice of $k$ and $m$: large $k$ values lead to the bias correction being equal to zero with large probability, and large values of $m$ result in $H_{k:m}(X,Y)$ being similar to an estimator obtained from a long MCMC run. Thus we expect the variance of $H_{k:m}(X,Y)$ to be similar to that of MCMC estimators for appropriate choices of $k$ and $m$.

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 (e.g. Owen, 2017). 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.

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

$\displaystyle\hat{\pi}$

$\displaystyle=\frac{1}{m-k+1}\sum_{\ell=k}^{m}\delta_{X_{\ell}}+\sum_{\ell=k+1}^{\tau-1}\min\left(1,\frac{\ell-k}{m-k+1}\right)(\delta_{X_{\ell}}-\delta_{Y_{\ell-1}}),$

(2.2)

obtained by replacing test function evaluations by delta masses in (2.1), as in Section 4 of Glynn and Rhee (2014). The measure $\hat{\pi}$ is of the form $\hat{\pi}=\sum_{\ell=1}^{N}\omega_{\ell}\delta_{Z_{\ell}}$ where the weights satisfy $\sum_{\ell=1}^{N}\omega_{\ell}=1$ and where the atoms $(Z_{\ell})$ are values among the history of the coupled chains. Some of the weights $(\omega_{\ell})$ may be negative, making $\hat{\pi}$ a signed empirical measure. In this view the unbiasedness property states $\mathbb{E}[\sum_{\ell=1}^{N}\omega_{\ell}h(Z_{\ell})]=\mathbb{E}_{\pi}[h(X)]$ for a test function $h$.

One can consider the convergence behavior of $\hat{\pi}^{R}=R^{-1}\sum_{r=1}^{R}\hat{\pi}^{(r)}$ towards $\pi$, where $(\hat{\pi}^{(r)})$ for $r\in\{1,\ldots,R\}$ are independent replications of $\hat{\pi}$. Glynn and Rhee (2014) obtain a Glivenko–Cantelli result for a similar measure related to their estimator. In the current setting, assume for simplicity that $\pi$ is univariate or else consider only one of its marginals. To emphasize the importance of the number of replications $R$, we rewrite the weights and atoms as $\hat{\pi}^{R}=\sum_{\ell=1}^{N_{R}}\omega_{\ell}\delta_{Z_{\ell}}$. Introduce the function $s\mapsto\hat{F}^{R}(s)=\sum_{\ell=1}^{N_{R}}\omega_{\ell}\mathds{1}(Z_{\ell}\leq s)$ on $\mathbb{R}$. Proposition 3.2 below states that $\hat{F}^{R}$ converges to $F$ as $R\to\infty$ uniformly with probability one, where $F$ is the cumulative distribution function of $\pi$.

The function $s\mapsto\hat{F}^{R}(s)$ is not monotonically increasing because of negative weights among $(\omega_{\ell})$, which motivates the following comments regarding the estimation of quantiles of $\pi$. Assume from now on that the pairs $(\omega_{\ell},Z_{\ell})$ are ordered such that $Z_{\ell}\leq Z_{\ell+1}$. For any $q\in(0,1)$ there might more than one index $\ell$ such that $\sum_{i=1}^{\ell-1}\omega_{i}\leq q$ and $\sum_{i=1}^{\ell}\omega_{i}>q$; the quantile estimate might be defined as $Z_{\ell}$ for any such $\ell$. The convergence of $\hat{F}^{R}$ to $F$ indicates that all such estimates are expected to converge to the $q$-th quantile of $\pi$. Therefore the signed measure representation leads to a way of estimating quantiles of the target distribution in a consistent way as $R\to\infty$. The construction of confidence intervals for these quantiles, perhaps by bootstrapping the $R$ independent copies, stands as an interesting area for future research. Another route to estimate quantiles of $\pi$ would be to project marginals of $\hat{\pi}^{R}$ onto the space of probability measures, for instance using a generalization of the Wasserstein metric to signed measures (Mainini, 2012). One could also estimate $F$ using isotonic regression (Chatterjee et al., 2015), considering $\hat{F}^{R}(s)$ for various values $s$ as noisy measurements of $F(s)$.

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

$\displaystyle H_{k:m}(X,Y)$

$\displaystyle=\text{MCMC}_{k:m}+\text{BC}_{k:m},$

where $\text{MCMC}_{k:m}$ is the MCMC average $(m-k+1)^{-1}\sum_{\ell=k}^{m}h(X_{\ell})$ and $\text{BC}_{k:m}$ is the bias correction term. The variance of $H_{k:m}(X,Y)$ can be written

$\displaystyle\mathbb{V}[H_{k:m}(X,Y)]$

$\displaystyle=\mathbb{E}\left[(\text{MCMC}_{k:m}-\mathbb{E}_{\pi}[h(X)])^{2}\right]+2\mathbb{E}\left[(\text{MCMC}_{k:m}-\mathbb{E}_{\pi}[h(X)])\text{BC}_{k:m}\right]+\mathbb{E}\left[\text{BC}_{k:m}^{2}\right].$

Defining the mean squared error of the MCMC estimator as $\text{MSE}_{k:m}=\mathbb{E}\left[(\text{MCMC}_{k:m}-\mathbb{E}_{\pi}[h(X)])^{2}\right]$, the Cauchy–Schwarz inequality yields

$$\mathbb{V}[H_{k:m}(X,Y)]\leq\text{MSE}_{k:m}+2\sqrt{\text{MSE}_{k:m}}\sqrt{\mathbb{E}\left[\text{BC}_{k:m}^{2}\right]}+\mathbb{E}\left[\text{BC}_{k:m}^{2}\right].$$

(3.1)

To bound $\mathbb{E}[\text{BC}_{k:m}^{2}]$, we introduce a geometric drift condition on the Markov kernel $P$.

We refer the reader to Meyn and Tweedie (2009) for the definitions and core theoretical tools for working with Markov chains on a general state space, in particular Chapter 5 for aperiodicity, $\varphi$-irreducibility and small sets, and Chapter 15 for geometric drift conditions; see also Roberts and Rosenthal (2004). Geometric drift conditions are known to hold for various MCMC algorithms (e.g. Roberts and Tweedie, 1996a, b; Jarner and Hansen, 2000; Atchade, 2006; Khare and Hobert, 2013; Choi and Hobert, 2013; Pal and Khare, 2014). Assumption 3.1 often plays a central role in establishing geometric ergodicity (e.g. Theorem 9 in Roberts and Rosenthal, 2004). We show next that this assumption enables an informative bound on $\mathbb{E}[\text{BC}_{k:m}^{2}]$.

Using Proposition 3.3, equation (3.1) becomes

$$\mathbb{V}[H_{k:m}(X,Y)]\leq\text{MSE}_{k:m}+2\sqrt{\text{MSE}_{k:m}}\frac{\sqrt{C_{\delta,\beta}\delta_{\beta}^{k}}}{m-k+1}+\frac{C_{\delta,\beta}\delta_{\beta}^{k}}{(m-k+1)^{2}}.$$

(3.2)

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$.

In order to facilitate the comparison between the efficiency of $H_{k:m}(X,Y)$ and that of MCMC estimators, we add simplifying assumptions. First, the right-most terms of (3.2) decrease geometrically with $k$, at a rate driven by $\delta_{\beta}=\delta^{1-2\beta}$ where $\delta$ is as in Assumption 2.2. This motivates a choice of $k$ depending on the distribution of the meeting time $\tau$. In practice, we can sample independent realizations of the meeting time and choose $k$ such that $\mathbb{P}(\tau>k)$ is small; i.e. we choose $k$ as a large quantile of the meeting times.

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

$$\mathbb{E}[2(\tau-1)+\max(1,m+1-\tau)]\mathbb{V}[H_{k:m}(X,Y)]\lessapprox\left(m+\mathbb{E}(\tau)\right)\mathbb{V}[H_{k:m}(X,Y)]\\ \lessapprox(m-k+1)\text{MSE}_{k:m}\left(1+\frac{k+\mathbb{E}(\tau)}{m-k+1}\right)\left(1+\frac{2}{\sqrt{m-k+1}}\sqrt{\frac{C_{\delta,\beta}\delta_{\beta}^{k}}{(m-k+1)\text{MSE}_{k:m}}}\right).$$

As $k$ increases we expect $(m-k+1)\text{MSE}_{k:m}$ to converge to $\mathbb{V}[(m-k+1)^{-1/2}\sum_{t=k}^{m}h(X_{t})]$, where $X_{k}$ would be distributed according to $\pi$. Denote this variance by $V_{k,m}$. The limit of $V_{k,m}$ as $m\to\infty$ is the asymptotic variance of the MCMC estimator, denoted by $V_{\infty}$. Hence, for $k$ and $m-k$ both large, the loss of efficiency of the method compared to standard MCMC is approximately $1+(k+\mathbb{E}(\tau))/(m-k)$.

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.

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}$.

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 (e.g Andrieu and Vihola, 2015), leading to polynomial tails for the associated meeting times.

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, 1991).

We denote the proposed estimator by $H$ and its expectation, which is the object of interest here, by $\pi(h)$. Generating $H$ takes a random time $C$. We write $N(t)$ for the number of independent copies of $H$ that can be produced by time $t$. The sequence $(H_{n},C_{n})_{n\in\mathbb{N}}$ refers to i.i.d. copies of $(H,C)$, so that we can write $N(t)=\sup\{n\geq 0:C_{1}+\ldots+C_{n}\leq t\}$, with $N(t)=0$ if $t<C_{1}$. We add the subscript $\mathrm{p}$ to refer to objects associated with processor $\mathrm{p}\in\{1,\ldots,\mathrm{P}\}$.

A first result is that the estimator $\bar{H}_{\mathrm{p}}(t)$, defined for all $1\leq\mathrm{p}\leq\mathrm{P}$ as $0$ if $N_{\mathrm{p}}(t)=0$, and by the sample average of $H_{\mathrm{p}1},\ldots,H_{\mathrm{p}N_{\mathrm{p}}(t)}$ otherwise, is biased: $\mathbb{E}[\bar{H}_{\mathrm{p}}(t)]=\mathbb{E}[H]-\mathbb{E}[H\mathds{1}(C>t)]$. Corollary 6 of Glynn and Heidelberger (1990) states that if $\mathbb{E}[|H\exp(\alpha C)|]<\infty$ for some $\alpha>0$, then the bias is negligible compared to $\exp(-\alpha t)$ as $t\to\infty$. By Cauchy–Schwarz, $\mathbb{E}[|H\exp(\alpha C)|]^{2}$ is less than the product of $\mathbb{E}[H^{2}]$ and $\mathbb{E}[\exp(2\alpha C)|]$. In our context, $\mathbb{E}[H^{2}]$ is finite under Proposition 3.1, and $\mathbb{E}[\exp(2\alpha C)|]$ is finite for a range of values of $\alpha$ that depends on the value of $\delta$ in Assumption 2.2.

We can define an unbiased estimator of $\pi(h)$ with a slight modification of $\bar{H}_{\mathrm{p}}(t)$. For all ${\mathrm{p}\in\{1,\dots,\mathrm{P}\}}$, set $\tilde{H}_{\mathrm{p}}(t)=\bar{H}_{\mathrm{p}}(t)$ if $N_{\mathrm{p}}(t)>0$ and $\tilde{H}_{\mathrm{p}}(t)=H_{\mathrm{p}1}$ if $N_{\mathrm{p}}(t)=0$. With $\tilde{N}_{\mathrm{p}}(t)=\max(1,N_{\mathrm{p}}(t))$, then $\tilde{H}_{\mathrm{p}}(t)$ is the sample average of $H_{\mathrm{p}1},\ldots,H_{\mathrm{p}\tilde{N}_{\mathrm{p}}(t)}$. The computation of $\tilde{N}_{\mathrm{p}}(t)$ requires the completion of $H_{\mathrm{p}1}$, and thus we cannot necessarily return $\tilde{H}_{\mathrm{p}}(t)$ at time $t$, in contrast with $\bar{H}_{\mathrm{p}}(t)$. On the other hand, we have $\mathbb{E}[\tilde{H}_{\mathrm{p}}(t)]=\mathbb{E}[H]=\pi(h)$, i.e. the estimator is unbiased, provided that $\mathbb{E}[|H|]<\infty$ (Corollary 7 of Glynn and Heidelberger (1990)). We denote the average of $\tilde{H}_{\mathrm{p}}(t)$ over $\mathrm{P}$ processors by $\tilde{H}(\mathrm{P},t)=\mathrm{P}^{-1}\sum_{\mathrm{p}=1}^{\mathrm{P}}\tilde{H}_{\mathrm{p}}(t)$, which is unbiased for $\pi(h)$.

Asymptotic results on $\tilde{H}(\mathrm{P},t)$ can be found in Glynn and Heidelberger (1991) and are summarized below. We first have the consistency results: $\lim_{t\to\infty}\tilde{H}(\mathrm{P},t)=\lim_{\mathrm{P}\to\infty}\tilde{H}(\mathrm{P},t)=\pi(h)$ almost surely for all $t$ and $\mathrm{P}$, and if $\mathbb{E}[|H|^{1+\delta}]$ for some $\delta>0$ and if $\{t_{\mathrm{P}}\}$ is a sequence such that $\lim_{\mathrm{P}\to\infty}t_{\mathrm{P}}=\infty$, then $\tilde{H}(\mathrm{P},t_{\mathrm{P}})$ converges to $\pi(h)$ in probability as $\mathrm{P}\to\infty$. Next, we can construct confidence intervals for $\pi(h)$ based on $\tilde{H}(\mathrm{P},t)$, following the end of Section 3 in Glynn and Heidelberger (1991). Indeed, define

	$\displaystyle\hat{\sigma}_{1}^{2}(\mathrm{P},t)$	$\displaystyle=\frac{1}{\mathrm{P}-1}\sum_{\mathrm{p}=1}^{\mathrm{P}}\left(\tilde{H}_{\mathrm{p}}(t)-\tilde{H}(\mathrm{P},t)\right)^{2},\quad\tilde{\tau}(\mathrm{P},t)=\frac{1}{\mathrm{P}}\sum_{\mathrm{p}=1}^{\mathrm{P}}\frac{1}{\tilde{N}_{\mathrm{p}}(t)}\sum_{n=1}^{\tilde{N}_{\mathrm{p}}(t)}C_{\mathrm{p}n},$
	$\displaystyle\hat{\sigma}_{2}^{2}(\mathrm{P},t)$	$\displaystyle=\tilde{\tau}(\mathrm{P},t)\left[\left(\frac{1}{\mathrm{P}}\sum_{\mathrm{p}=1}^{\mathrm{P}}\frac{1}{\tilde{N}_{\mathrm{p}}(t)}\sum_{n=1}^{\tilde{N}_{\mathrm{p}}(t)}H_{\mathrm{p}n}^{2}\right)-\tilde{H}(\mathrm{P},t)^{2}\right],$

where $\tilde{N}_{\mathrm{p}}(t)=\max(1,N_{\mathrm{p}}(t))$. Then we have the three following central limit theorems,

	$\displaystyle\text{for fixed $t$ and $\mathrm{P}\to\infty$},\quad\frac{\sqrt{\mathrm{P}}}{\hat{\sigma}_{1}(\mathrm{P},t)}\left(\tilde{H}(\mathrm{P},t)-\pi(h)\right)$	$\displaystyle\to$	$\displaystyle\mathcal{N}(0,1),$		(3.4)
	$\displaystyle\text{for fixed $\mathrm{P}$ and $t\to\infty$},\quad\frac{\sqrt{\mathrm{P}t}}{\hat{\sigma}_{2}(\mathrm{P},t)}\left(\tilde{H}(\mathrm{P},t)-\pi(h)\right)$	$\displaystyle\to$	$\displaystyle\mathcal{N}(0,1),$		(3.5)
	$\displaystyle\text{if $t_{\mathrm{P}}\to\infty$ as $\mathrm{P}\to\infty$},\quad\frac{\sqrt{\mathrm{P}t_{\mathrm{P}}}}{\hat{\sigma}_{2}(\mathrm{P},t_{\mathrm{P}})}\left(\tilde{H}(\mathrm{P},t_{\mathrm{P}})-\pi(h)\right)$	$\displaystyle\to$	$\displaystyle\mathcal{N}(0,1).$		(3.6)

These results require moment conditions such as $\mathbb{E}[\tilde{H}_{\mathrm{p}}(t)^{2}]<\infty$. The central limit theorem in (3.4) will be used to construct confidence intervals in Sections 5.3 and 5.4.

We conclude this section with a remark on the setting where $t$ is fixed and the number of processors $\mathrm{P}$ goes to infinity. There, the time to obtain $\tilde{H}(\mathrm{P},t)$ would typically increase with $\mathrm{P}$. Indeed at least one estimator needs to be completed on each processor for $\tilde{H}(\mathrm{P},t)$ to be available. The completion time behaves as the maximum of independent copies of the cost $C$. Under Assumption 2.2, the completion time for $\tilde{H}(\mathrm{P},t)$ has expectation behaving as $\log\mathrm{P}$ when $\mathrm{P}\to\infty$, for fixed $t$. Other tail assumptions (Middleton et al., 2018) would lead to different behavior for the completion time associated with $\tilde{H}(\mathrm{P},t)$.

A maximal coupling between two distributions $p$ and $q$ on a space $\mathcal{X}$ is a distribution of a pair of random variables $(X,Y)$ that maximizes $\mathbb{P}(X=Y)$, subject to the marginal constraints $X\sim p$ and $Y\sim q$. We write $p$ and $q$ both for these distributions and for their probability density functions with respect to a common dominating measure, and refer to the uniform distribution on the interval $[a,b]$ by $\mathcal{U}([a,b])$. A procedure to sample from a maximal coupling is described in Algorithm 2; see e.g. Section 4.5 of Chapter 1 of Thorisson (2000), and Johnson (1998) where it is termed $\gamma$-coupling.

We justify Algorithm 2 and compute its cost. Denote by $(X,Y)$ the output of the algorithm. First, $X$ follows $p$ from step 1. To prove that $Y$ follows $q$, introduce a measurable set $A$. We write $\mathbb{P}(Y\in A)=\mathbb{P}(Y\in A,\text{step 1})+\mathbb{P}(Y\in A,\text{step 2})$, where the events $\{\text{step 1}\}$ and $\{\text{step 2}\}$ refer to the algorithm terminating at step 1 or 2. We compute

$\displaystyle\mathbb{P}\left(Y\in A,\text{step 1}\right)$

$\displaystyle=\int_{A}\int_{0}^{+\infty}\mathds{1}\left(w\leq q(x)\right)\frac{\mathds{1}\left(0\leq w\leq p(x)\right)}{p(x)}p(x)dwdx=\int_{A}\min(p(x),q(x))dx.$

We can deduce from this that $\mathbb{P}\left(\text{step 1}\right)=\int_{\mathcal{X}}\min(p(x),q(x))dx$. For $\mathbb{P}\left(Y\in A,\text{step 2}\right)$ to be equal to $\int_{A}(q(x)-\min(p(x),q(x)))dx$, we need

$\displaystyle\int_{A}(q(x)-\min(p(x),q(x)))dx$

$\displaystyle=\mathbb{P}\left(Y\in A|\text{step 2}\right)\left(1-\int_{\mathcal{X}}\min(p(x),q(x))dx\right),$

and we conclude that the distribution of $Y$ given $\{\text{step 2}\}$ should for all $x$ have a density $\tilde{q}(x)$ equal to ${(q(x)-\min(p(x),q(x)))/(1-\int\min(p(x^{\prime}),q(x^{\prime}))dx^{\prime})}$. Step 2 is a standard rejection sampler using $q$ as a proposal distribution to target $\tilde{q}$, which concludes the proof that $Y\sim q$. We also confirm that Algorithm 2 maximizes the probability of $\{X=Y\}$. Under the algorithm,

$$\mathbb{P}(X=Y)=\mathbb{P}(\text{step 1})=\int_{\mathcal{X}}\min(p(x),q(x))dx=1-d_{\text{TV}}(p,q),$$

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.