Ergodicity of Approximate MCMC Chains with Applications to Large Data Sets

Introduction

Markov chain Monte Carlo (MCMC) sampling is an indispensable tool for Bayesian computation. Most Metropolis-Hastings samplers require full evaluation of the posterior at two points at every step; some MCMC samplers require even more information, such as gradients of the likelihood function. In many modern applications, the computational cost of this full evaluation of the likelihood function can be prohibitively large. For instance, a full likelihood evaluation might involve processing a massive amount data, or computing the solution of a partial differential equation representing an underlying physical phenomenon. The result is that the inferential performance of naive implementations of MCMC algorithms can actually deteriorate as the amount of data grows, unless available computational resources grow even more quickly. The same problem has been noted for other techniques that are popular in computational statistics (see e.g., [Jord13, SVA15] for broader discussions of this problem in a similar spirit, outside of the context of MCMC). As mentioned in those papers, it might be possible to avoid actual loss of performance if we are aware that it is a possibility. However, it is much harder to ensure that data is being used efficiently under severe computational constraints.

Simulating a Markov chain that merely approximates the desired MCMC dynamics is an increasingly common approach to circumvent this difficulty. These approximations often rely on estimating, rather than evaluating, the posterior distribution of interest - for example, doing so based on a subsample of the available data [Beau03, OBBEM00, WeTe11, KCW13, BDC14, quiroz2014speeding, AFEB14, maire2015light]. While many approximate MCMC methods seem to be very successful in practice, they do not have the same convergence guarantees as standard MCMC samplers. In this paper, we present some general convergence results for such ‘approximate’ MCMC algorithms and discuss their applications to some recently proposed algorithms. Our results give quantitative bounds on the convergence in distribution of the Markov chain as well as convergence of finite samples drawn from the Markov chain. These bounds allow us to provide advice on how to choose parameters for various approximation schemes. As a consequence, they also give modest conditions under which certain approximation schemes are more efficient than their underlying MCMC algorithms. This result is unsurprising, but it is quite hard to prove for many existing approximation schemes for reasons discussed at the start of Section 4. While our applications focus on the problems posed by large data sets, our theoretical bounds are also relevant to MCMC samplers targetting intractable likelihoods; see [MCPS13] for applications of related ideas in that context.

Our paper has three main contributions. The first is providing a variety of robust ergodicity results for perturbed Markov chains, comparing the mixing properties and Monte Carlo errors of a perturbed chain to its base chain. The second is an application of these bounds to certain ‘interpolating’ chains in order to obtain a bias-variance tradeoff inequality. This inequality provides some advice on how to choose the best approximate sampler for a given computational budget. In particular, it gives sufficient conditions under which particular approximate MCMC samplers give smaller Monte Carlo errors than the ‘correct’ Metropolis-Hastings dynamics. The third is a discussion of the limits of our perturbation-based analysis for studying the tradeoff between computational and statistical efficiency for Markov chains, with a focus on some ‘obvious’ facts that cannot be proved in this framework.

2. Guide to the Paper

We begin by setting up some notation in Section 2, and prove our general perturbation bounds in Section 3. We apply our bounds to the austerity framework of [KCW13] in Section 4. This application begins with an introduction to the austerity framework and some definitions related to the idea of ‘computational complexity’ of MCMC. We continue by providing a general bias-variance tradeoff result in this framework. This bound provides sufficient conditions under which an approximate MCMC kernel is better than the true Metropolis-Hastings kernel. In Sections LABEL:SecCompChoice and LABEL:SecJust we note that the perturbation results in Section 3 give very poor bounds if applied directly to the true Metropolis-Hastings dynamics, and introduce interpolating chains and new perturbation bounds that allow us to avoid this problem. In Section LABEL:SecBigMisc, we discuss the sharpness of our results and point out some interesting phenomena that seem difficult to study via perturbation analysis.

3. Related Literature

The earlier papers [KCW13, Mitr05, FHL13] also give abstract convergence results for perturbations of Markov chains, though the first two primarily restrict their attention to uniformly ergodic chains and the latter does not give the required quantitative bounds for the types of comparisons we do in this paper.

Since writing our first version of this paper, a number of other preprints and articles studying the perturbation theory of Markov chains and applications to ‘big data’ have been released. This includes [BDC14, AFEB14, quiroz2014speeding, zhu2014big, friel2015exploiting, rudolf2015perturbation, maire2015light, bardenet2015markov, chen2014sublinear, chen2015subsampling, green2015bayesian, johndrow2015approx], the latter 9 of which refer to our earlier draft. We briefly discuss some relationships between these papers and our work.

In [BDC14, AFEB14], the authors discuss or introduce various approximate MCMC algorithms and prove that they are in fact small perturbations of the ‘true’ MCMC dynamics under various conditions. This justifies the use of perturbation analysis to bound the bias of these approximate MCMC samplers. However, when applied to subsampling MCMC algorithms, the bounds in both papers typically apply only when the average subsample size $n$ is on the same order as the total amount of data $N$. [quiroz2014speeding] introduces another subsampling MCMC algorithm and gives the first bounds on the decay rate of the bias of their algorithm that depends only on the subsample size $n$, not the total amount of data $N$. However, these bounds do not give an explicit decay rate in terms of $n$. The excellent survey article on approximate MCMC [bardenet2015markov] also provides arguments on the convergence of subsampling MCMC algorithms when $n\ll N$, among many other results.

The recent article [johndrow2015approx] introduces several new approximate MCMC algorithms, and also has some focus on testing when their approximate MCMC algorithms are ‘better’ than the default MCMC methods. This is an important aspect that we do not address in our paper. Understanding this issue for a wide class of algorithms in the key next step. The authors in [johndrow2015approx] provide a criterion that is very similar to the criterion introduced in Section 4.2 of this paper; but their criteria makes sense even in the regime $n\ll N$.

[rudolf2015perturbation] gives bounds on convergence of perturbed Markov chains under very general conditions, the main subject of our Section 3. Some of their bounds, like some of ours, are based on the ‘curvature’ of the original Markov chain. However, they consider many more technical conditions than we do.The other papers cited are also about approximate MCMC, but have core messages that do not greatly overlap with ours. We emphasize that some of these subsampling algorithms (e.g., [maire2015light, chen2014sublinear]) do not obviously require $n\approx N$ in order to be useful.

Of these papers, we feel that results in [quiroz2014speeding] are closest in spirit to our conclusions. Our main purpose in writing this note was to show that subsampling MCMC schemes can give much smaller Monte Carlo errors than ‘correct’ Metropolis-Hastings schemes, even when the amount of data $N$ is very large to the total computational resources available. Among other contributions, [quiroz2014speeding] showed that this was plausible by bounding the bias of subsampling MCMC in a way that depended only on the subsample size $n$, not the total amount of data available $N$. They also discussed notions of computational complexity that are similar to ours. Our analysis differs from theirs in a number of ways. Most obviously, the bias bounds in [quiroz2014speeding] were based on the pseudomarginal algorithm; our bounds apply to approximate chains that are not pseudomarginal. Technically, our bounds are obtained via perturbation theory, and we introduce new ‘interpolating’ chains that allow the powerful perturbation techniques to be applied to the small-$n$, large-$N$ regime. Finally, we bound both the bias and the mixing properties of approximate chains, while [quiroz2014speeding] focuses on the bias.

Finally, we compare and contrast our results to that of [AFEB14]. Most of the perturbation bounds in [AFEB14] are closely related to those in [Mitr05]. Although our general bounds are similar in spirit to [AFEB14], our arguments are quite different and there are situations where our bounds apply and theirs do not (and vice versa). See Remark 3.4 for a longer discussion on the relationship between our results and that of [Mitr05].

Notation

Our convergence results, like many in the Markov chain literature, are described in terms of the Wasserstein distance. For a pair of measures $\mu,\nu$ on a Polish space $(\Omega,d)$, let $\Pi(\mu,\nu)$ be the set of all couplings of $\mu$ and $\nu$. Then the Wasserstein distance between $\mu,\nu$ is given by

W_d(μ,ν) = inf_ζ∈Π(μ,ν) ∫_x,y ∈Ω d(x,y) ζ(dx,dy).

We study the convergence of Markov chains through the notions of mixing times, curvature and drift functions. Recall that the mixing time of a Markov chain $\{X_{t}\}_{t\geq 0}$ with kernel $K$ and stationary distribution $\pi$ on state space $\Omega$ is given by {equs}τ_mix= sup_X_0 = x ∈Ω inf{ t ≥0 : ∥ L(X_t) - π∥_TV ¡ 14 }.

We follow [Olli09] in defining the Ricci curvature of a transition kernel $K$ on $(\Omega,d)$ by {equs}κ(x,y) = 1 - Wd(K(x, ⋅), K(y, ⋅))d(x,y), and the curvature of the entire chain is defined to be {equs}κ= inf_x,y ∈Ω κ(x,y).

A positive curvature for a power $K^{s}$ of $K$ implies that $K$ has good mixing properties. For example, if the curvature of $K^{s}$ with respect to the metric $d(x,y)=\textbf{1}_{x\neq y}$ is $\kappa$, it is straightforward to check that the mixing time of $K$ is at most $\lceil\frac{s\,\log(4)}{\log((1-\kappa)^{-1})}\rceil\approx\frac{s}{\kappa}\log(4)$. It is worth noting that in many cases of interest, it is sufficient to calculate $\kappa(x,y)$ for $d(x,y)$ ‘small’; see, e.g., Prop 19 of [Olli09].

Following [JoOl10], we also define the eccentricity of a point $x\in\Omega$ by: {equs}E(x) = ∫_Ω d(x,y) π(dy).

Technical Results

We begin with a general mixing estimate that most of our bounds will follow from. First, an assumption:

This assumption is generally easy to check using the same calculations that one uses to establish a drift condition in the usual sense (see inequality (2)). Indeed, as long as the Lyapunov function is not pathological, this assumption will be implied by the usual drift condition.

In the common situation that a contraction bound holds uniformly, this gives:

Assume that $K$ has eccentricity $E(x)<\infty$ with respect to $\pi$. Assume that {equs} sup_x,y ∈Ω Wd( K(x,⋅), K(y,⋅))d(x,y) ≤(1 - α) and that inequality (3.3) is satisfied for some $0<\delta<\infty$.

Taking powers of $K$ and noting that all kernels are contractive in Total Variation distance, Lemma 3.3 has the immediate corollary:

We mention that a result very similar to Corollaries 1, 2 is also immediately implied by Corollary 3.1 of [Mitr05]; the constants are also very similar in our regime of interest. Although our Lemma 3.3 and the results in [Mitr05] both imply the same result in this restricted setting, they have different emphases. Their result is based on linear algebra; ours is purely probabilistic. Our results also apply to chains that are not uniformly ergodic, as chains can have Wasserstein contraction (or Total Variation contraction on small sets) without being uniformly ergodic. Finally, their result only applies for convergence in Total Variation, while our results explicitly allow the use of many other Wasserstein metrics; this flexibility can lead to bounds that are effectively much sharper if the metric is chosen carefully.

This last difference is most easily seen when the Markov chain $K$ satisfies inequality (3.3) for some fixed $\alpha>0$ throughout a non-compact state space, and for which the eccentricity satisfies $E(x)<\infty$ for each $x\in\Omega$ but $\sup_{x\in\Omega}E(x)=\infty$. These chains are generally geometrically ergodic but not uniformly ergodic. Lemma 3.3 can provide direct bounds on their finite-time bias while Corollary 3.1 of [Mitr05] does not apply. This class of Markov chains contains many Markov chains of interest. For example, it includes the Metropolis-Hastings chain with proposal distribution $L(x,\cdot)=\mathcal{N}(x,\sigma_{1})$ and target stationary distribution $\pi=\mathcal{N}(\mu,\sigma_{2})$ as long as $\sigma_{1}<\sigma_{2}$.

2. Convergence of Monte Carlo Estimates

We begin by showing that the bias is small. To state our result, we recall the definition of the trace $\{X_{t}^{(S)}\}$ of a Markov chain $\{X_{t}\}_{t\geq 0}$ onto a set $S$. Let $T(S)=\{t\geq 0\,:\,X_{t}\in S\}$ be an ordered set. When $|T(S)|=\infty$, we define {equs}{X_t^(S) }_t ≥0 = { X_t }_t ∈T(S), again viewed as an ordered set. We have:

For further study of this question, see e.g. [rudolf2015perturbation, FHL13].

Bias-Variance Tradeoff and Applications to Austerity Framework

In this section, we recall various ‘approximate’ MCMC chains, introduce measures of computational and statistical efficiency, give a general tradeoff result for these measures, and apply the tradeoff result to certain examples.

We consider a data set $\mathcal{X}=\{x_{1},\ldots,x_{N}\}$ and are interested in sampling from a posterior distribution {equs} π(θ) ≡π(θ—{x_i }_i=1^N) = p(θ) ∏_i=1^N p(θ—x_i).

One standard tool is the Metropolis-Hastings algorithm. Fix a reversible transition kernel $L$ on state space $\Omega$. We recall the following form of the Metropolis-Hastings algorithm from [BDC14]:

Algorithm 2, first suggested in [BDC14], is a representative approximate MCMC algorithm. The associated constants are {equs}f_t^∗ &= t-1N C_θ,θ’ = max_1 ≤i ≤N — log( p(x_i — θ’)) - log( p(x_i — θ)) — c_t = C_θ,θ’ 2(1-ft∗) log(2 / δt)t and any number $0<\gamma<\infty$.

All of the algorithms stated in this section, and in later sections, are defined for any reversible kernel $L$. In practice, we are most interested in proposal kernels of the form {equs} L(x,y) ≡L_σ(x,y) = f_σ(x-y), where $f_{\sigma}$ is a density with variance $\sigma$. Both in practice and all our theorems, the choice of $\sigma$ will be adapted to the variance of the stationary distribution of the underlying Markov chain - the variance of the proposal distribution should increase with the variance of the stationary distribution. In particular, a subsampling MCMC algorithm should generally have a proposal distribution $f_{\sigma}$ with a larger variance than the optimal variance of the original Metropolis-Hastings dynamics. See e.g. [rosenthal2011optimal] for a survey on the optimization of proposal distributions.

2. Measures of Computational and Statistical Efficiency

We introduce measures of computational efficiency for approximate MCMC, in the same tradition as earlier work such as [sherlock2014efficiency, doucet2015efficient, bornn2014use, quiroz2014speeding].