Unbiased approximation of posteriors via coupled particle Markov chain Monte Carlo

MCMC is a powerful methodology for the approximation of complex distributions. MCMC is intrinsically iterative and, while asymptotically unbiased, the size of the bias and the Monte Carlo error of generated estimates given a finite number of iterations are often difficult to quantify. Moreover, MCMC will typically not allow full exploitation of the computational potential of modern distributed-computing techniques. Recently, Jacob et al. (2020b) propose a method for unbiased MCMC estimation based on coupling of Markov chains, building on ideas by Glynn and Rhee (2014). The algorithm is embarrassingly parallel, and the unbiasedness provides immediate quantification of the Monte Carlo error. However, devising effective coupling for MCMC algorithms targeting a given posterior can be highly challenging. See e.g. the construction of coupled MCMC for a horseshoe regression model in Biswas et al. (2021) and the development for posteriors based on Hamiltonian Monte Carlo (HMC) in Heng and Jacob (2019).

Considering the scope for unbiased MCMC and its current limitations, we propose a coupled MCMC algorithm where the coupling mechanism is not in principle specific to the posterior at hand, resulting in a general recipe for unbiased MCMC. We do so by working on an appropriate augmented space. Specifically, we build on recent advances for unbiased estimation in state-space models (Middleton et al., 2019; Jacob et al., 2020a) which devise coupling strategies through particle filters independently of the shape of the target distribution. Our work extends these ideas by embedding a general posterior distribution in a state-space model using adaptive SMC (Del Moral et al., 2006). This results in a methodology that broadens the class of posteriors that can be treated via unbiased MCMC, by exploiting the coupling methods feasible within the SMC framework. As our methodology couples Markov chains that consist of particle MCMC methods, we refer for convenience to our method as ‘coupled particle MCMC’. To illustrate the potential of our methodology, we apply coupled particle MCMC to (i) Gaussian graphical models, which are substantially more complex than models considered previously in the literature on unbiased MCMC, mainly due to discreteness and dimension of graph space; (ii) horseshoe regression.

The MCMC step resulting from our state space embedding combines particle MCMC methods by Middleton et al. (2019) and Jacob et al. (2020a), namely coupled particle independent Metropolis-Hastings (PIMH) and conditional SMC, respectively. The focus of these works is unbiased approximation for the smoothing distribution of a state-space model; the high-dimensionality therein relates to the number of time steps, whereas the state space is implicitly assumed to be low-dimensional. We consider more complex state spaces of a much higher dimension than Middleton et al. (2019) and Jacob et al. (2020a), provide effective adaptation of the SMC, and consider both PIMH and conditional SMC to improve mixing (as compared to using only PIMH) of the MCMC. We empirically investigate the trade-off between mixing and coupling via a number of numerical experiments.

The structure of the paper is as follows. Section 2 reviews unbiased MCMC estimation based on coupled Markov chains. Section 3 introduces coupled particle MCMC for unbiased estimation, for the case of a general posterior. Section 4 contains theoretical results for the meeting time of the coupled Markov chains and the resulting unbiased estimator deduced from the literature on coupled conditional SMC for state-space models (Lee et al., 2020). Section 5 applies the proposed general methodology to simulated data, including a case with horseshoe regression, where a comparison with the MCMC coupling approach in Biswas et al. (2021) is carried out. Section 6 considers Gaussian graphical models as an example where effective coupling of MCMC is highly non-trivial. We finish with a discussion in Section 7.

Before we introduce coupled particle MCMC, we review coupled MCMC and introduce notation. Denote the parameter space by $\mathcal{X}\subseteq\mathbb{R}^{d_{x}}$, $d_{x}\geq 1$. Let $x\in\mathcal{X}$ denote the parameter, $y\in\mathcal{Y}\subseteq\mathbb{R}^{d_{y}}$, $d_{y}\geq 1$, the data, and $\pi(x)$ the density of the posterior of interest w.r.t. some dominating measure. The unbiased construction requires a pair of coupled ergodic Markov chains $\{x(t-1),\bar{x}(t)\}$, $t\geq 1$, on $\mathcal{X}\times\mathcal{X}$ with both chains having $\pi(x)$ as equilibrium distribution. The coupling is such that $x(t-1)$ and $\bar{x}(t)$ have the same distribution for any $t\geq 1$, and the two chains meet at some random time $\tau<\infty$ a.s. so that $x(t)=\bar{x}(t)$ for $t\geq\tau$. Under standard conditions, the posterior expectation $\pi(h)=\int_{\mathcal{X}}h(x)\,\pi(x)\,dx$ of a statistic $h:\mathcal{X}\mapsto\mathbb{R}$ can be obtained as $\pi(h)=\lim_{t\to\infty}E\left[h\{x(t)\}\right]$, with all expectations assumed finite. Writing the limit as a telescoping sum and using that fact that $x(t-1)$, $\bar{x}(t)$ admit the same law, gives for any fixed $k\geq 0$ (Glynn and Rhee, 2014)

We assume that the technical conditions that permit the exchange of summation and expectation in the last step hold in our setting. Thus, the quantity

is an unbiased estimator of $\pi(h)$. Process $\{\bar{x}(t)\}$ is typically initialised at $\bar{x}(1)=x(0)\sim p_{0}$, for a law $p_{0}$. Both chains will evolve marginally according to the same MCMC chain with the posterior $\pi(x)$ as invariant distribution, with a coupling applied for the joint transition $[\,x(t),\bar{x}(t)\mid x(t-1),\bar{x}(t-1)\,]$, $t\geq 2$. We adopt this setting for the rest of the paper.

Coupled particle MCMC requires some further notation. Denote the prior density on the parameter by $p(x)$. Let $p(y\mid x)$ be the density of the data $y$ given $x$. The posterior density follows from Bayes’ rule as $\pi(x)\propto p(x)\,p(y\mid x)$. For any $\alpha\in[0,1]$, we denote by $\pi_{\alpha}(x)$ the tempered posterior proportional to $p(x)\,p(y\mid x)^{\alpha}$. Thus, $\pi_{0}(x)=p(x)$ and $\pi_{1}(x)=\pi(x)$. All densities are assumed to be determined w.r.t. appropriate reference measures. As our method involves adaptive SMC, we assume we can sample $x$ from its prior $p(x)$ and evaluate the likelihood $p(y\mid x)$. The method requires, for any $\alpha\in(0,1]$, the construction of an MCMC step with $\pi_{\alpha}(x)$ as its invariant distribution. We refer to this step as the ‘inner’ MCMC step as it is a part of a more involved ‘outer’ particle MCMC step. For integers $i\leq j$, denote the range $\{i,\dots,j\}$ by ${i\!:\!j}$. We use the colon notation for collections of variables, i.e., $x_{i:j}=\{x_{i},\dots,x_{j}\}$ and $x^{i:j}=\{x^{i},\dots,x^{j}\}$.

Existing analysis of coupled conditional SMC applies to our context which aims to approximate the posterior $\pi(x)$. Building on Lee et al. (2020), we derive results for Algorithm 6 with $\rho=0$, such that only conditional SMC is used. Appendix B contains the proofs which mostly consist of a mapping from the smoothing problem considered in Lee et al. (2020) to our context of posterior approximation. That mapping results in the following assumptions.

Many models satisfy Assumption 1, including the Gaussian graphical model in Section 6. The likelihood from linear regression in Section 5.2 violates it if the number of predictors is greater than the number of observations. The assumption relates to the boundedness of potential functions in the SMC literature, which is a common assumption (Del Moral, 2004, Section 2.4.1). Andrieu et al. (2018, Theorem 1) show that Assumption 1 is essentially equivalent to uniform ergodicity of the Markov chains produced by conditional SMC in Algorithm 6 with $\rho=0$.

In the context of particle filtering, Jacob et al. (2020a, Section 3) establish unbiasedness and finite variance of the estimator $\hat{h}_{k}$ in (1), like we do, but without the boundedness in Assumption 2 and instead use an assumption jointly on $h(x)$ and the Markov chain generated by conditional SMC. The simpler and more restrictive Assumption 2 provides a bound on the variance of $\hat{h}_{k}$ in terms of the number of particles $N$ in Proposition 2.

Assumptions 1 and 2 imply Assumptions 6 and 4 of Middleton et al. (2019), respectively. Therefore, the estimator $\hat{h}_{k}$ from Algorithm 8 with $\rho=1$, such that it uses only PIMH, is unbiased and has finite variance, and $\tau<\infty$ almost surely per Proposition 8 of Middleton et al. (2019) and its proof under Assumptions 1 and 2.

Proposition 1 contains no conditions on the quality of the coupling of the inner MCMC steps in Step 4c of Algorithm 4. This confirms that the SMC machinery by itself enables coupling and thus unbiased estimation for the posterior $\pi(x)$, though coupling of the inner MCMC steps can lead to substantially smaller meeting times as explored empirically in Appendix D, an aspect of the method that the above theoretical results fail to capture. Moreover, the meeting time $\tau$ can be made equal to 2 with arbitrarily high probability by increasing the number of particles $N$.

Whether a sufficient increase of $N$ is practicable depends on the SMC sampler in Algorithm 2 which we adapt to $\pi(x)$. For instance, under additional assumptions, Theorem 9 from Lee et al. (2020) states that the meeting probability $\mathrm{Pr}(x_{0:S}^{\prime}=\bar{x}_{0:S}^{\prime})$ does not vanish if the number of particles scales as $N=\mathcal{O}(2^{S}S)$ where $S$ is the number of temperatures.

Proposition 2 confirms the unbiasedness of $\hat{h}_{k}$. Moreover, the variance penalty resulting from the bias correction sum in (1) vanishes as $N\to\infty$ or $k\to\infty$. The latter confirms the role of $k$ as the number of burn-in iterations in the MCMC estimation. The constant $c$ in Proposition 2 is the same as in Proposition 1. The rates of convergence implied by the upper bounds in Propositions 1 and 2 are probably conservative as noted by Lee et al. (2020).

In all our applications, Algorithm 1 sets up a Feynman-Kac model via a preliminary run that uses $N_{0}=10^{4}$ particles, an ESS threshold of $\gamma_{0}=0.8$ to determine $\alpha_{0:S}$ and a correlation threshold of $\zeta_{0}=0.95$ for $m_{0:S}$. Subsequent executions of SMC algorithms use adaptive resampling with ESS threshold $\gamma=0.5$, i.e. we resample, (Step 2b of Algorithm 2) only if ESS drops below $\gamma N$. All empirical results use Rao-Blackwellization when computing $\bar{h}_{k}^{l}$.

To evaluate algorithmic performance for a test function $h$, we consider the product ‘$\hat{\operatorname{var}}(\bar{h}_{k}^{l})\times\text{time}$’ of the empirical variance of $\bar{h}_{k}^{l}$ and the average computation time to obtain one $\bar{h}_{k}^{l}$ across the number, $R$, of independent runs. This product captures the trade-off between quantity and quality of unbiased estimates: the variance of the average of independently and identically distributed estimates $\bar{h}_{k}^{l}$ is proportional to $\hat{\operatorname{var}}(\bar{h}_{k}^{l})$ and inversely proportional to the number of estimates. The number of estimates is in turn inversely proportional to the average computation time when working with a fixed computational budget. Thus, ‘$\hat{\operatorname{var}}(\bar{h}_{k}^{l})\times\text{time}$’ is proportional to the variance of the final estimator for fixed computation time. We compute ‘$\hat{\operatorname{var}}(\bar{h}_{k}^{l})\times\text{time}$’ for $l=1,\dots,{l_{\max}}$, for some $l_{\max}$, where $k\leq l$ is set for each $l$ such that $\hat{\operatorname{var}}(\bar{h}_{k}^{l})$ is minimized.

The mixing of the outer MCMC and the meeting times both influence $\operatorname{var}(\bar{h}_{k}^{l})$ with larger meeting times and worse mixing corresponding to higher variance, though the manner in which such effects combine is non-trivial. Thus, we present mixing and meeting times in our empirical results for further insight. Mixing is monitored via calculation of the integrated autocorrelation time for each run based on the chain $\sum_{i=1}^{N}w^{i}(t)\,h\{x_{S}^{i}(t)\}/\sum_{i=1}^{N}w^{i}(t)$, $t=l/2,\dots,l$, computed via the R package LaplacesDemon (Statisticat, LLC., 2020).