On Markov chain Monte Carlo methods for tall data

Markov chain Monte Carlo methods are often deemed too computationally intensive to be of any practical use for big data applications, and in particular for inference on datasets containing a large number n of individual data points, also known as tall datasets. In scenarios where data are assumed independent, various approaches to scale up the Metropolis-Hastings algorithm in a Bayesian inference context have been recently proposed in machine learning and computational statistics. These approaches can be grouped into two categories: divide-and-conquer approaches and, subsampling-based algorithms. The aims of this article are as follows. First, we present a comprehensive review of the existing literature, commenting on the underlying assumptions and theoretical guarantees of each method. Second, by leveraging our understanding of these limitations, we propose an original subsampling-based approach which samples from a distribution provably close to the posterior distribution of interest, yet can require less than O(n) data point likelihood evaluations at each iteration for certain statistical models in favourable scenarios. Finally, we have only been able so far to propose subsampling-based methods which display good performance in scenarios where the Bernstein-von Mises approximation of the target posterior distribution is excellent. It remains an open challenge to develop such methods in scenarios where the Bernstein-von Mises approximation is poor.

Authors

• 15 publications
• 74 publications
• 30 publications
08/30/2021

Lagged couplings diagnose Markov chain Monte Carlo phylogenetic inference

Phylogenetic inference is an intractable statistical problem on a comple...
01/14/2015

Unbiased Bayes for Big Data: Paths of Partial Posteriors

A key quantity of interest in Bayesian inference are expectations of fun...
07/14/2017

Big Data vs. complex physical models: a scalable inference algorithm

The data torrent unleashed by current and upcoming instruments requires ...
09/10/2021

Diagnostics for Monte Carlo Algorithms for Models with Intractable Normalizing Functions

Models with intractable normalizing functions have numerous applications...
10/14/2021

Divide-and-Conquer Monte Carlo Fusion

Combining several (sample approximations of) distributions, which we ter...
11/07/2017

Bayesian Inference of Selection in the Wright-Fisher Diffusion Model

The increasing availability of population-level allele frequency data ac...
01/19/2010

Bayesian inference for queueing networks and modeling of internet services

Modern Internet services, such as those at Google, Yahoo!, and Amazon, h...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Performing inference on tall datasets, that is datasets containing a large number of individual data points, is a major aspect of the big data challenge. Statistical models, and Bayesian methods in particular, commonly demand Markov chain Monte Carlo (MCMC) algorithms to make inference, yet running MCMC on such tall datasets is often far too computationally intensive to be of any practical use. Indeed, MCMC algorithms such as the Metropolis-Hastings (MH) algorithm require at each iteration to sweep over the whole dataset. Frequentist or variational Bayes approaches are thus usually preferred to a fully Bayesian analysis in the tall data context on computational grounds. However, they might be difficult to put in practice or justify in scenarios where the likelihood function is complex; e.g. non-differentiable (Chernozhukov & Hong, 2003). Moreover, some applications require precise quantification of uncertainties and a full Bayesian approach might be preferable in those instances. This is the case for example for applications from experimental sciences, such as cosmology (Trotta, 2006) or genomics (Wright, 2014)

, were such big data problems abound. Consequently, much efforts have been devoted over recent years to develop scalable MCMC algorithms. These approaches can be broadly classified into two groups: divide-and-conquer approaches and subsampling-based algorithms. Divide-and-conquer approaches divide the initial dataset into batches, run MCMC on each batch separately, and then combine these results to obtain an approximation of the posterior: Subsampling approaches aim at reducing the number of individual data point likelihood evaluations necessary at each iteration of the MH algorithm.

After briefly reviewing the limitations of MCMC for tall data, introducing our notation and two running examples in Section 2, we first review the divide-and-conquer literature in Section 3. The rest of the paper is devoted to subsampling approaches. In Section 4, we discuss pseudo-marginal MH algorithms. These approaches are exact in the sense that they target the right posterior distribution. In Section 5, we review other exact approaches, before relaxing exactness in Section 6. Throughout, we focus on the assumptions and guarantees of each method. We also illustrate key methods on two running examples. Finally, in Section 7, we improve over our so-called confidence sampler in (Bardenet et al. , 2014), which samples from a controlled approximation of the target. We demonstrate these improvements yield significant reductions in computational complexity at each iteration in Section 8. In particular, our improved confidence sampler can break the barrier of number of individual data point likelihood evaluations per iteration in favourable cases. Its main limitation is the requirement for cheap-to-evaluate proxies for the log-likelihood, with a known error. We provide examples of such proxies relying on Taylor expansions.

All examples can be rerun or modified using the companion IPython notebook222The IPython notebook and a static html render of it can both be found at http://www.2020science.net/research/scaling-mcmc-methods. to the paper, available as supplementary material.

2 Bayesian inference, MCMC, and tall data

In this section, we describe the inference problem of interest and the associated MH algorithm. We also detail the two running examples on which we benchmark key methods in Section 4, 5 and 6.

2.1 Bayesian inference

Consider a dataset

 \cX={x1,...,xn}⊂\Xset⊂Rd, (1)

and a parameter space . We assume the data are conditionally independent with associated likelihood given a parameter value and we denote the associated average log-likelihood

 ℓ(θ)=1nn∑i=1logp(xi|θ)=1nn∑i=1ℓi(θ). (2)

We follow a Bayesian approach where one assigns a prior to the unknown parameter, so that inference relies on the posterior distribution

 π(θ)=p(θ|x)∝γ(θ)≜p(θ)enℓ(θ), (3)

where denotes an unnormalized version of . In most applications, is intractable and we will focus here on Markov chain Monte Carlo methods (MCMC; Robert & Casella, 2004) and, in particular, on the Metropolis-Hastings (MH) algorithm to approximate it.

2.2 The Metropolis-Hastings algorithm

A standard approach to sample approximately from is to use MCMC algorithms. To illustrate the limitation of MCMC in the tall data context, we focus here on the MH algorithm (Robert & Casella, 2004, Chapter 7.3). The MH algorithm simulates a Markov chain of invariant distribution . Then, under weak assumptions, see e.g. (Douc et al. , 2014, Theorem 7.32)

, the following central limit theorem holds for suitable test functions

 √Niter⎡⎢⎣1NiterNiter% ∑k=0h(θk)−∫h(θ)π(θ)dθ⎤⎥⎦→\cN(0,σ2lim(h)), (4)

where convergence is in distribution.

The pseudocode of MH targeting a generic distribution is given in Figure 1. In the case of Bayesian inference with independent data (3), Step 1 is equivalent to setting

 logα(θ,θ′) = log[p(θ′)p(θ)q(θ|θ′)q(θ′|θ)]+n[ℓ(θ′)−ℓ(θ)]. (5)

When the dataset is tall (), evaluating the log likelihood ratio in (5) is too costly an operation and rules out the applicability of such a method. As we shall see, two possible options are to either divide the dataset into tractable batches, or approximate the acceptance ratio in (5) using only part of the dataset.

2.3 Running examples

We will evaluate some of the described approaches on two illustrative simple running examples. We fit a one-dimensional normal distribution

to i.i.d. points drawn according to and lognormal observations , respectively. The latter example illustrates a misspecification of the model. We assign a flat prior

. For all algorithms, we start the chain at the maximum a posteriori (MAP) estimate. The MH proposal is an isotropic Gaussian random walk, whose stepsize is first set proportional to

and then adapted during the first iterations so as to reach acceptance. When applicable, we also display the number of likelihood evaluations per iteration, and compare it to the evaluations required at each iteration by the MH algorithm.

In Figure 2, we illustrate the results of iterations of vanilla MH on each of the two datasets. MH does well, as the posterior coincides with that of a longer reference run of iterations in each case, and the autocorrelations show a fast exponential decrease. The Bernstein-von Mises approximation (van der Vaart, 2000, Chapter 10.2), a Gaussian centered at the true value, with covariance minus the scaled inverse Fisher information, is a very good approximation to the posterior in both cases. We are thus in simple cases of heavy concentration of the posterior, where subsampling should help a lot if it is to be of any help in tackling tall data problems.

3 Divide-and-conquer approaches

A natural way to tackle tall data problems is to divide the data into batches, run MH on each batch separately, and then combine the results.

3.1 Representing the posterior as a combination of batch posteriors

Assume data are divided in batches . Relying on the equality

 p(θ|\cX)∝B∏i=1p(θ)1/Bp(xi|θ), (6)

Huang & Gelman (2005) propose to combine the batch posterior approximations using Gaussian approximations or importance sampling. Scott et al.  (2013) propose to average samples across batches, noting this is exact under Gaussian assumptions. Neiswanger et al.  (2014) propose to run an MCMC chain on each batch targeting an artificial batch posterior

 πi(θ)∝p(θ)1/Bp(xi|θ),

fit a smooth approximation to each batch posterior, and multiply them. These methods are however theoretically justified only when batch posteriors are Gaussian, or when the size of each batch goes to infinity, to guarantee that the used smooth approximation of each batch posterior is accurate.

There are few results available on how the properties of combined estimators scale with the number of batches . Neiswanger et al.  (2014) fit a kernel density estimator to the samples of each batchwise chain, and multiply the resulting kernel density estimators. A sample from this mixture approximation to is then obtained through an additional MCMC step. Under simplifying assumptions (all MCMC chains are assumed being independent draws from their targets, for example), a bound on the MSE of the final estimator is obtained. However, this bound explodes as the kernel bandwidth goes to zero, and more importantly, it is exponential in the number of batches . In a tall data context, the number of batches is expected to grow with to ensure that the size of each batch is less than . Thus, the proposed bound is currently not informative for tall data.

As pointed out by Wang & Dunson (2013), if the supports of the are almost disjoint, then the product of their approximations will be a poor approximation to . To improve the overlap between the approximations of the ’s, Wang & Dunson (2013) propose to replace the posterior in (6) by the product of the Weierstrass transforms of each batch posterior. When the approximation of is an empirical measure, its Weierstrass transform corresponds to a kernel density estimator. The product of the Weierstrass transforms can be interpreted as the marginal distribution of an extended distribution on , where the first copies of are associated to the batches, and the remaining copy is conditionally Gaussian around a weighted mean of the first copies. Unfortunately, sampling from the posterior of this artificial model is difficult when one only has access to approximate samples of each .

Although it is not strictly speaking a Monte Carlo method, we note that Xu et al.  (2014) and Gelman et al.  (2014) propose an expectation-propagation-like algorithm that similarly tackles the issue of disjoint approximate batch posterior supports. Each batch of data points is represented by its individual likelihood times a cavity distribution. The cavity distribution is itself the product of the prior and a number of terms that represent the contributions of other batches to the likelihood. The algorithm iterates between 1) simulating from each batchwise likelihood times a batch-specific cavity distribution, and 2) fitting each batch-specific cavity component. Again, while these approaches are computationally feasible and appear to perform well experimentally, it is difficult to assert the characteristics of the proposed approximation of the posterior and there is no convergence guarantee available for this iterative algorithm.

3.2 Replacing the posterior by a geometric combination of batch posteriors

Another avenue of research consists in avoiding multiplying the batch posteriors by replacing the target by a different combination of the latter.

By introducing a suitable metric on the space of probability measures such as the Wasserstein metric, it is for example possible to define the barycenter or the median of a set of probability measures.

Minsker et al.  (2014) propose to output the median of the batch posteriors, while Srivastava et al.  (2014) use the Wasserstein barycenter, which can be computed efficiently in practice using the techniques developed by Cuturi & Doucet (2014). While this idea has some appeal, the statistical meaning of these median or mean measures is unclear, and the robustness of the median estimate advocated in (Minsker et al. , 2014) may also be a drawback, as in some circumstances valuable information contained in some batches may be lost.

To conclude, divide-and-conquer techniques appear as a natural approach to handle tall data. However, the crux is how to efficiently combine the batch posterior approximations. The main issues are that the batch posterior approximations potentially have disjoint supports, that the multiplicative structure of the posterior (3) leads to poor scaling with the number of batches, that theoretical guarantees are often asymptotic in the batch size, and that cheap-to-sample combinations of batch posteriors are difficult to interpret.

4 Exact subsampling approaches: Pseudo-marginal MH

Pseudo-marginal MH (Lin et al. , 2000; Beaumont, 2003; Andrieu & Roberts, 2009)

is a variant of MH, which relies on unbiased estimators of an unnormalized version of the target. Pseudo-marginal MH is useful to help understand several potential approaches to scale up MCMC. We start by describing pseudo-marginal MH in Section

4.1. Then, we present two pseudo-marginal approaches to tall data in Section 4.2 and Section 4.3.

4.1 Pseudo-marginal Metropolis-Hastings

Assume that instead of being able to evaluate , we have access to an unbiased, almost-surely non-negative estimator of the unnormalized target . Pseudo-marginal MH substitutes a realization of to in Step 1. Similarly, it replaces in Step 1 by the realization of that was computed when the parameter value was last proposed. Pseudo-marginal MH is of considerable practical importance, with applications such as particle marginal MH (Andrieu et al. , 2010) and MCMC versions of the approximate Bayesian computation paradigm (Marin et al. , 2012). It is thus worth investigating its use in the context of tall data problems.

The possibility to use an unbiased estimator of

comes at a price: first, the asymptotic variance

in (4) of an MCMC estimator based on a pseudo-marginal chain will always be larger than that of an estimator based on the underlying “marginal” MH (Andrieu & Vihola, 2015). Second, the qualitative properties of the underlying MH may not be preserved, meaning that the rate of convergence to the invariant distribution may go from geometric to subgeometric, for instance; see Andrieu & Roberts (2009) and Andrieu & Vihola (2015) for a detailed discussion. In practice, if the variance of is large for some value , then an MH move to might be accepted while largely overestimates . In that case, it is difficult for the chain to leave , and pseudo-marginal MH chains thus tend to get stuck if the variance of the involved estimators is not controlled. When some tunable parameter allows to control this variance, Doucet et al.  (2015) show that, in order to minimize the variance of MCMC estimates for a fixed computational complexity, the variance of the log-likelihood estimator should be kept around when the ideal MH having access to the exact likelihood generates quasi-i.i.d samples from ; or set to around when it exhibits very large integrated autocorrelation times. In practice, the integrated autocorrelation times of averages under the ideal MH are unknown as this algorithm cannot be implemented. In this common scenario, Doucet et al.  (2015) recommend keeping the variance around 1.5 as this is a value which ensures a small penalty in performance even in scenarios where 1.0 or 3.0 are actually optimal. They also show that the penalty incurred for having a variance too small (i.e. inferior to 0.2) or too large (i.e. superior to 10) is very large. When mentioning pseudo-marginal MH algorithms, we will thus comment on the variance of the logarithm of the involved estimators , or, if not available, of their relative variance.

4.2 Unbiased estimation of the likelihood using unbiased estimates of the log-likelihood

As described in Section 4.1, pseudo-marginal MH requires an almost-surely nonnegative unbiased estimator of the unnormalized posterior at , for any in . It is easy to check that, by sampling from the dataset with or without replacement, we obtain the following unbiased estimator of the log-likelihood

 n^ℓ(θ)=ntt∑i=1logp(x∗i|θ). (7)

We denote by the subsampling estimate of the average log-likelihood and denote by its variance. Obviously, exponentiating (7) does not provide an unbiased estimate of the likelihood However, an interesting question is whether one can design a procedure which outputs an unbiased, almost-surely nonnegative estimate of using unbiased estimates of such as . Without making any further assumption about , it was recently shown by Jacob & Thiery (2013) that it is not possible. However, this can be done if one further assumes, for instance, that there exists such that for all , see (Jacob & Thiery, 2013, Section 3.1) who rely on a technique by Rhee & Glynn (2013) generalizing (Bhanot & Kennedy, 1985). Unfortunately, as we shall see, the resulting estimator typically has a very large relative variance, resulting in very poor performance of the associated pseudo-marginal chain.

We apply (Rhee & Glynn, 2013, Theorem 1) to build an unbiased non-negative estimator of , which is equivalent to defining . For , let

 D∗j=ntt∑i=1logp(x∗i,j|θ)−na(θ), (8)

be an unbiased estimator of , where the ’s are drawn with replacement from for each , and are further independent across . In (Jacob & Thiery, 2013, Section 3.1),

is an integer-valued random variable whose tails do not decrease too fast, in the sense that

. To ease computations, we take to be geometric with parameter . This corresponds to the lightest tails allowed by (Jacob & Thiery, 2013, Section 3.1), since . Finally, let

 Y\defeqena(θ)[1+N∑k=11P(N≥k)1k!k∏j=1D∗j]. (9)

By (Rhee & Glynn, 2013, Theorem 1), is a non-negative unbiased estimator of the likelihood . As mentioned in Section 4.1, it is crucial, if we want to plug in a pseudo-marginal algorithm, to control the variance of its logarithm. The variance of is difficult to compute, so we use here the relative variance of as a proxy.

Proposition 4.1

Let and be the almost surely non-negative estimator of defined in (9). Then its relative variance satisfies

 \VarYe2nℓ(θ)≥e−2n(ℓ(θ)−a(θ))+2n√(1+\eps)[σt(θ)2+(ℓ(θ)−a(θ))2]n√(1+\eps)[σt(θ)2+(ℓ(θ)−a(θ))2]+\cO(1). (10)

The proof of Proposition 4.1 can be found in Appendix A. We can interpret (10) as follows: in order for the relative variance of not to increase exponentially with , it is necessary that is of order . But if of order , so that the batchsize would have to be of order , which is impractical. It is also necessary that is of order to control the term in . This means that should be taken of order , but then the mean of the geometric variable will be of order . This entails that the number of terms in the randomly truncated series (9) should be of order , which defeats the purpose of using this estimator.

Hence for the reasons outlined in Section 4.1, we expect the pseudo-marginal MH relying on to be highly inefficient. Indeed, we have not been able to obtain reasonably mixing chains even on our Gaussian running example. We have experimented with various choices of , and with various values of , but none yielded satisfactory results. We conclude that this approach is not a viable solution to MH for tall data.

We note that Strathmann et al.  (2015) have recently proposed a different way to exploit the methodology of Rhee & Glynn (2013) in the context of tall data. However, their methodology does not provide unbiased estimates of the posterior expectations of interest. It only provides unbiased estimates of some biased MCMC estimates of these expectations, these MCMC estimates corresponding indeed to running an MCMC kernel on the whole dataset for a finite number of iterations. Strathmann et al.  (2015) suggest that it might be possible to combine their algorithm with the recent scheme of Glynn & Rhee (2014) to obtain unbiased estimates of the posterior expectations. It is yet unclear whether this could be achieved under realistic assumptions on the MCMC kernel.

4.3 Building ^γ(θ) with auxiliary variables

In (MacLaurin & Adams, 2014), the authors propose an alternative MCMC to sample from which, similarly to the method described previously, only requires evaluating the likelihood of a subset of the data at each iteration. Assume a bound is available for each and . For simplicity, we further assume that only depends on through . This is the case in the experiments of (MacLaurin & Adams, 2014), as well as ours. Note also that in Section 4.2, we used a bound that was uniform in the data index ; we could have used similarly a non-uniform bound, but this would have made the derivation of Proposition 4.1 unnecessarily heavy.

As noted in (MacLaurin & Adams, 2014), we can then define the following extended target distribution on

 \tpi(θ,z) ∝ p(θ)n∏i=1[exp(ℓi(θ)−exp(bi(θ))]ziexp(bi(θ))1−zi (11) = p(θ)n∏i=1exp(bi(θ))n∏i=1[exp(ℓi(θ)−bi(θ))−1]zi.

This distribution satisfies two important features: it admits as a marginal distribution, and its pointwise evaluation only requires to evaluate for those ’s for which . Note that evaluating however requires to evaluate , and the bounds

thus must be chosen so that this computation is cheap. This is the case for the lower bound of the logistic regression log-likelihood model discussed in

(MacLaurin & Adams, 2014), which is a quadratic form in , where is the label of datum . The idea of replacing the evaluation of the target by a Bernoulli draw and the evaluation of a lower bound has been exploited previously; see e.g. (Mak, 2005).

Any MCMC sampler could be used to sample from . MacLaurin & Adams (2014) propose an MH-within-Gibbs sampler that leverages the known conditional . The expected cost of one conditional MH iteration on at equilibrium, that is the average number of indices such that , is , and the constant is related to the expected relative tightness of the bound, see (MacLaurin & Adams, 2014, Section 3.1). The number of likelihood evaluations for an update of conditional on is explicitly controlled in (MacLaurin & Adams, 2014) by either specifying a maximum number of attempted flips, or implicitly specifying the fraction of flips to .

The authors of MacLaurin & Adams (2014) remarked that their methodology is related to pseudo-marginal techniques but did not elaborate. We show here how it is indeed possible to exploit the extended target distribution in (11) to obtain an unbiased estimate of an unnormalized version of . More precisely, we have

 p(xi|θ)=∑zi∈{0,1}p(xi,zi|θ)

where . Hence, the marginal distribution of under this extended model is given by

 p(zi=1|θ) = ∫p(zi=1,xi|θ)dxi (12) = ∫[exp(ℓi(θ))−exp(bi(θ))]dxi = 1−Iθ,

where

. Using Bayes’ theorem, we obtain accordingly

 p(xi|θ,zi=1)=exp(ℓi(θ))−exp(bi(θ))1−Iθ, p(xi|θ,zi=0)=exp(bi(θ))Iθ.

An obvious unbiased estimator of the unnormalized posterior is thus given by

 ^γ(θ)=p(θ)n∏i=1p(xi|θ,zi) (13)

where each is drawn independently given from (12). Note that in the case of logistic regression, if is chosen to be the quadratic lower bound given in (MacLaurin & Adams, 2014), its integral is a Gaussian integral and can thus be computed. Finally, similarly to the Firefly algorithm of MacLaurin & Adams (2014), the number of evaluations of the likelihood per iteration is , loosely speaking.

Although the pseudo-marginal variant of Firefly we propose has the disadvantage of requiring the integrals to be tractable, it comes with two advantages. First, the sampling of does not require to evaluate the likelihood at all. If computing all bounds does not become a bottleneck, this avoids the need to explicitly state a resampling fraction at the risk of augmenting the variance of the likelihood estimator. Second, the properties of this variant are easier to understand, as it is a ‘standard’ pseudo-marginal MH and hence the results from Section 4.1 apply. In particular, although it has the correct target distribution, the asymptotic variance of ergodic averages is inflated compared to the ideal algorithm.

As explained in Section 4.1, we consider the variance of the log likelihood estimator.

Proposition 4.2

Let . With the notations introduced in Section 4.3,

 \Varz[n∑i=1logp(xi|θ,zi)]=Iθ(1−Iθ)n∑i=1log2[Iθ1−Iθ(eℓi(θ)−bi(θ)−1)] (14)

The proof of Proposition 4.2 can be found in Appendix B. Proposition 4.2 can be interpreted as follows: the variance is related to how tight the bound is. In general, obtaining a variance of order will only be possible if most bounds are very tight, and the bigger

, the tighter the bounds have to be. These conditions will typically not be met when a fixed fraction of “outlier

’s correspond to untight bounds.

We give the results of the original Firefly MH on our running Gaussian and log normal examples in Figure 3. We bound each using a 2nd order Taylor expansion at the MLE and the Taylor-Lagrange inequality, see Section 7.2.1 for further details. This bound is very tight in both cases, so that we are in the favourable case where only a few components of are at each iteration, and the number of likelihood evaluations per full joint iteration is thus roughly the fraction of points for which has been resampled. We chose the fraction of resampled points to be here, and initialized to have of ones. Trying smaller fractions led to very slowly mixing chains even for the Gaussian example. Estimating the number of likelihood evaluations per full joint iteration as the sum of the number of resampled ’s and the number of “bright” points, we obtained in both the Gaussian and lognormal case an almost constant number of likelihood evaluations close to , so that only a few points are bright. This can be explained by the tightness of the Taylor bound, which leads Firefly MH to almost exclusively replace the evaluation of the likelihood by that of the Taylor bound. Finally, unlike the other algorithms we applied, we observed that a bad choice of the initial value of can easily take out of the posterior mode. To be fair, we thus discarded the first iterations as a burn-in before plotting.

As expected, the algorithm behaves erratically in the lognormal case, as failure to attempt a flip of each draws the -component of the chain towards the few large values of which are bright. Since the bright points are rarely updated, the chain mixes very slowly.

5 Other exact approaches

Other exact approaches have been proposed, which do not rely on pseudo-marginal MH.

5.1 Forgetting about acceptance: stochastic approximation approaches

Welling & Teh (2011) proposed an algorithm based on stochastic gradient Langevin dynamics (SGLD). This is an iterative algorithm which at iteration uses the following update rule

 θk+1=θk+\epsk+12[∇logp(θ)+ntt∑i=1∇logp(x∗i,k|θ)]+√\epsk+1ηk+1, (15)

is a sequence of time steps, are independent vectors and

 ntt∑i=1∇logp(x∗i,k|θ)

is an unbiased estimate of the score computed at each iteration using a random subsample of the observations. This approach is reminiscent of the Metropolis-adjusted Langevin algorithm (MALA; Robert & Casella 2004, Section 7.8.5), where the proposal given by

 θ′=θ+\eps2[∇logp(θ)+n∑i=1∇logp(xi|θ)]+√\epsη,

is used in an MH acceptance step, where . The point of Welling & Teh (2011) is that if one suppresses the MH acceptance step, computes an unbiased estimate of the score but introduces a sequence of stepsizes that decreases to zero at the right rate, then

 ⎛⎝Niter∑k=0\epsk⎞⎠−1N% iter∑k=0\epskδθk

is an approximation to . The algorithm has been analyzed recently in (Teh et al. , 2014), where it has been established that it provides indeed a consistent estimate of the target. Additionally, a central limit theorem holds with convergence rate , which is slower than the traditional Monte Carlo rate . It is yet unclear how SGLD compares to other subsampling schemes in theory: it may require a smaller fraction of the dataset per iteration, but more iterations are needed to reach the same accuracy.

In practice, we show the results of SGLD on our two running examples in Figure 4. The stepsize is chosen proportional to , following the recommendations of Teh et al.  (2014). We show the results of two choices for the subsample size : and of the data, with respectively and iterations, so that both runs amount to the same fraction of the budget of the vanilla MH in Figure 2. Both runs are still far from convergence on the lognormal example: subsampling draws the chain away from the support of the posterior, and one has to wait for smaller stepsizes to avoid overconfident moves. But then, the variance of the final estimate gets bigger. Constant stepsizes lead to comparable results (not shown).

Finally, we note that subsampling for Hamiltonian Monte Carlo (HMC; Duane et al. , 1987) has also been recently considered. Chen et al.  (2014) propose a modification of HMC that is inspired by the SGLD with decreasing stepsize of Welling & Teh (2011), while Betancourt (2014) explores why naive approaches suffer from unacceptable biases. The algorithm of Chen et al.  (2014)

is however a heuristic, which further relies on the subsampling noise being Gaussian. As demonstrated in

(Bardenet et al. , 2014) and in Section 6.2, relying on a Gaussian noise assumption can yield arbitrarily poor performance when this assumption is violated.

5.2 Delayed acceptance

Banterle et al.  (2015) remarked that if we decompose the acceptance ratio in a product of positive functions

 α(θ,θ′)=B∏i=1ρi(θ,θ′)

then the MH-like algorithm that accepts the move from to with probability

 B∏i=1min[ρi(θ,θ′),1]

still admits as invariant distribution. In practice, in the case of tall data, we can divide the dataset into batches and use for example

 ρi(θ,θ′)=p(θ′)1/Bp(xi|θ′)q(θ|θ′)p(θ)1/Bp(xi|θ)q(θ′|θ′).

This allows us to reject candidate without having to compute the full likelihoods and the calculations of can be done in parallel. However, as remarked by Banterle et al.  (2015), the resulting Markov chain has a larger asymptotic variance in (4) than the original MH, and it does not necessarily inherit the ergodicity of the original MH. Furthermore, by construction, every accepted point has to be evaluated on the whole dataset, and the average proportion of data points used is thus lower bounded by the acceptance rate of the algorithm, which in practice is often around . Overall, it is an easy-to-implement feature that does not add any bias, but its benefits are inherently limited, and speed of convergence might be affected.

6 Approximate subsampling approaches

In this Section, we consider again subsampling approaches where, at each MH iteration, a subset of data points is used to approximate the acceptance ratio (5). As mentioned in Section 4.2, it is very simple to obtain an unbiased estimator of the log-likelihood based on random samples from the dataset ; see (7). Similarly, one can also easily obtain an unbiased estimator of the average log likelihood ratio

 Λ∗t(θ,θ′)\defeq1tt∑i=1logp(x∗i|θ′)p(x∗i|θ). (16)

Note that unlike the exact approaches of Sections 4 and 5, the methods reviewed in Section 6 do not attempt to sample exactly from , but just from an approximation to .

6.1 Naive subsampling

The first approach one could try is to only use a random fixed proportion of data points to estimate at any newly drawn . We highlight that this leads to a nontrivial mixture target that is hard to interpret, where all subsampled posteriors appear, suitably rescaled. Assume that at each new drawn in an MH run, we draw independent Bernoulli variables and let

 ^ℓ(θ)=1nn∑i=1ziλℓi(θ) (17)

be an unbiased estimator of the average log likelihood , where i.i.d. Now one could think of plugging estimates in Steps 1 and 1 of MH in Figure 1. However, as is not an unbiased estimator of , this algorithm samples from a target distribution which is proportional to , where the expectation is w.r.t the distributions of the Bernoulli random variables . Now

 Een^ℓ(θ) = n∏i=1[λℓi(θ)1/λ+(1−λ)] = n∑r=0∑#Ir=rλr(1−λ)n−rp(xIr|θ)1/λ,

where denotes a set of distinct indices in , , and with the convention . Each subsampled likelihood contributes to the target, exponentiated to the power , resulting in a nontrivial mixture of rescaled data likelihood terms. To further simplify, assume for each set of indices , that is, the variance of the likelihood under subsampling is small, then

 Een^ℓ(θ)≈n∑r=0Crnλr(1−λ)n−rpr(θ)1/λ=ER∼B(n,λ)p1/λR(θ), (18)

where

denotes the binomial distribution with parameters

and . Noticing that is roughly exponentially decreasing in , the values or that are larger than the mode of the binomial probability mass function are unlikely to contribute a lot to (18). The largest subsample size contributing to (18) is thus roughly , and the power makes this term of the same scale as . Broadly speaking, subsampling has a “broadening” effect due to the contribution of the likelihoods of small subsamples.

Alternately, if one starts with the biased estimator of the average log likelihood

 ~ℓ(θ)=1nn∑i=0ziℓi(θ), (19)

instead of (17) one ends up with

 Een~ℓ(θ) = n∏i=1[λℓi(θ)+(1−λ)] (20) = n∑r=0∑#Ir=rλr(1−λ)n−rp(xIr|θ) ≈ ER∼B(n,λ)pR(θ).

Again, all subsampled likelihoods contribute, but without being further exponentiated. Still, the result is a much broadened target, as values of that are larger than are unlikely to contribute a lot. In this case, the broadening effect of subsampling is not only due to the contribution of small subsamples, but also to the absence of bias correction in (19).

We have thus seen that naive subsampling is nontrivial tempering, so that the target is not preserved. Additionally, as mentioned in Section 4.1, the variance of the log likelihood estimator needs to be kept around a constant, or depending on hypotheses, in order for pseudo-marginal MH to be efficient. This means that should be such that

 (1−λ)λn∑i=1ℓi(θ)2

is of order in the case of (17). This entails that should be close to 1, so that there is no substantial gain in terms of number of likelihood evaluations. In the case of (19),

 λ(1−λ)n∑i=1ℓi(θ)2

can be of order if . But then the leading terms in the mixture target (20) will be the subsampled likelihoods corresponding to small subsamples, so that the target will be very different from the actual target .

Overall, naive subsampling is a very poor approach. However, it allows us to identify the main issues a good subsampling approach should tackle: guaranteeing its target, not loosing too much convergence speed compared to MH, and cutting the likelihood evaluation budget. As shown in (Bardenet et al. , 2014), the first point is an algorithmic design issue, while the last two points are related to controlling the variance of the log likelihood ratios.

6.2 Relying on the CLT

Several authors have appealed to the central limit theorem to justify their assumption that the average subsampled log likelihoods and log likelihood ratios in (7) and (16

) are Gaussianly distributed.

If the noise of the log likelihood ratio estimate is normal with known variance and mean equal to the true log-likelihood ratio, Ceperley & Dewing (1999) have proposed an MH with a corrected acceptance ratio that is exact, i.e., that still targets . When the variance of the noise is not known, and is rather estimated, the method becomes inexact. Nicholls et al.  (2012) propose a heuristic argument to show that this inexact chain gives reasonable approximate results, but the Gaussian assumption remains crucial there. As shown in (Bardenet et al. , 2014) and in this paper in Figures 5 and 6, this assumption can be arbitrarily violated when subsampling tall data if the log likelihood ratios are heavy-tailed. Missing log likelihoods in the tails will lead to erroneous decisions, and uncontrolled results.

6.2.1 A pseudo-marginal approach with variance reduction under Gaussian assumption

Quiroz et al.  (2014) propose a methodology to use MH for tall data which also relies on the assumption that the log-likelihood estimator is Gaussian with mean , for every . By introducing a bias correction providing an approximately unbiased estimate of the likelihood, this corresponds to a pseudo-marginal MH algorithm whose target distribution is proportional to , where is an estimate of the bias satisfying . They rightly notice that, ideally, if one wants to keep the variance of average subsampled log likelihoods small, one should not subsample data points with or without replacement, but one should rather perform importance sampling with the weight of data point being proportional to . While this variance reduction approach obviously defeats the purpose of subsampling, Quiroz et al.  (2014) propose to use as weights an approximation of the log-likelihood, based e.g. on a Gaussian process or splines trained on a small subset of computed likelihoods . Finally, a heuristic to adaptively choose the size of the total subsample so as to keep the variance of the log likelihood controlled is proposed. The method is demonstrated to work on a bivariate probit model using only of the full dataset. However, as a general purpose method, it suffers from two limitations. First, it is based on Gaussian assumptions, which can be unreasonable as noted above and it is unclear whether it will be robust to these CLT approximations not being valid. Second, the proposed importance sampling step requires to learn a good proxy for for each drawn during the MCMC run. The fitted proxies should thus be cheap to train and evaluate, but at the same time accurate if any variance reduction is to be obtained.

Still assuming the noise of the log likelihood is Gaussian, given a drawn , one can try to adaptively choose the size of the subsample to be used in the unbiased estimators (7) or (16), so as to take the correct acceptance decision with high probability. Upon noting that the MH acceptance decision is equivalent to deciding whether , or equivalently

 1nn∑i=1logp(xi|θ′)p(xi|θ)>1nlogu−1nlog[p(θ′)q(θ|θ′)p(θ)q(θ′|θ)] (21)

with drawn beforehand, statistical tests can be used to assert whether (21) holds with a given level of “confidence”. As far as we are aware, Bulgak & Sanders (1988)

were the first to consider such a procedure. They used it in a simulated annealing algorithm maximizing a function defined as an expectation w.r.t a probability distribution, and approximated using Monte Carlo. Simulated annealing is a simple non-homogeneous variant of the MH algorithm where the the target distribution is annealed over the iterations. The same application received more attention later

(Alkhamis et al. , 1999; Wang & Zhang, 2006). Applied to the standard MH, the method has been considered by Singh et al.  (2012), and more recently by Korattikara et al.  (2014) specifically for tall data. Korattikara et al.  (2014) propose an MH-like algorithm called Austerity MH

that incorporates a sequential T-test to check (

21) for each pair , thus relying on several CLTs. They demonstrate dramatic reductions in the average number of subsamples used per MCMC iteration on particular applications. However, as noted in (Korattikara et al. , 2014; Bardenet et al. , 2014), the results can be arbitrarily far from the original MH when the CLT approximations are not valid.

We show the results of iterations of Austerity MH on our two running examples in Figure 5. The parameters are , corresponding to the p-value threshold in the aforementioned T-test, and an initial subsample size of at each iteration. In the Gaussian case, the posterior is rightly centered, but is slightly too wide. This is a tempering effect due to too small subsamples, while the CLT-based Student approximation seems reasonable, as shown in Figure 6. In the lognormal case, the departure of the chain from the actual posterior is more remarkable, and relatedly the CLT approximations of Austerity MH are inaccurate for the chosen initial subsample size of , as we demonstrate in Figure 6. This explains the strong mismatch of the chain and the posterior in Figure 5(b). The standard deviation of the fitted Gaussian is largely underestimated, due to small subsamples which do not include enough of the tails of the log likelihood ratios, which coincide with the tails of . Finally, the reductions in the number of samples needed per iteration are quite interesting: half of the iterations require less than

of the dataset for the lognormal case, but this is at the price of a large error in the posterior approximation. Augmenting the initial size of the subsample will likely make the CLT approximations tighter, but there is no generic answer as to which size to choose: any fixed choice will fail on an example where the log likelihood ratios have heavy enough tails. In both the Gaussian and the lognormal example, it is actually safer to go with the Bernstein-von Mises approximation, which costs little more than a run of stochastic gradient descent, and only requires one CLT approximation, for a sample of size

. This illustrates the danger of using CLT-based approximations for small sample sizes, which is related to asymptotic arguments on small batches in Section 3.

Overall, CLT-based approaches to MH with tall data lead to heuristics with interesting reductions in the number of samples used, but they have little theoretical backing so far and they are not robust to the involved CLT approximations being inaccurate. We note also that the CLT is assumed to provide a good approximation for the log likelihood or log likelihood ratio for any , which amounts to more than one Gaussian assumption. The approaches in this section should thus be applied with care. As a minimal sanity-check, we recommend using tests of Gaussianity across to make sure the CLT assumptions are realistic. Note that even then, there is no guarantee the above algorithms have for target, if any.

6.3 Exchanging acceptance noise for subsampling noise

This section is an original contribution, which illustrates a way to obtain subsampling algorithms with guarantees under weaker assumptions than Gaussianity. This approach is impractical, but it is of methodological and illustrative interest. First it illustrates a potentially useful technique to take advantage of subsampling noise. Second, it is our first illustration of the seemingly inevitable average number of subsamples required per MCMC iteration as soon as we do not use any CLT-based approximation and require theoretical guarantees.

Let , and let be drawn independently with replacement from . Let be the average subsampled log likelihood ratio defined in (16). Now, we remark that MH has some inherent noise in its acceptance decision (21), encapsulated by the uniform variable . Why, then, not rely on the subsampling noise to guarantee exploration, and accept a move if and only if

 Λ∗t(θ,θ′)+1nlog[p(θ′)q(θ|θ′)p(θ)q(θ′|θ)]>0 (22)

instead of (21)? This idea has been first used by Branke et al.  (2008)

to develop heuristics for simulated annealing in the presence of noise. We formalize this argument here in the context of subsampling. For the sake of simplicity, assume for a moment we have a flat prior and a symmetric proposal, so that (

22) becomes

 Λ∗t(θ,θ′)>0.

We do not assume that the ’s are Gaussianly distributed, but we make the parametric assumption that the second and third absolute moments and of are known and independent of . Applying the Berry-Esseen inequality (van der Vaart & Wellner, 1996) to the variables yields

 ∣∣ ∣∣P(−Λ∗t(θ,θ′)≤u)−Φ(u+Λn(θ,θ′)σ/√t)∣∣ ∣∣≤K(σ,ρ)√t (23)

for any , where is the cdf of a variable, and

 Λn(θ,θ′)\defeq1nn∑i=1logp(xi|θ′)p(xi|θ)

is the average log likelihood ratio. When , (23) yields

 ∣∣ ∣∣P(Λ∗t(θ,θ′)≥0)−Φ(Λn(θ,θ′)σ/√t)∣∣ ∣∣≤K(σ,ρ)√t. (24)

Now let be such that for any ,

 ∣∣∣Φ(x)−11+e−λx∣∣∣≤C.

Bowling et al.  (2009) for instance, empirically found and . Combining this bound with (24), we obtain

 ∣∣ ∣ ∣∣P(Λ∗t(θ,θ′)≥0)−11+e−λΛn(θ,θ′)σ/√t∣∣ ∣ ∣∣≤C+K(σ,ρ)√t.

Hence, the acceptance probability of an algorithm that would accept the move from to by checking whether is close to the acceptance probability of an MCMC algorithm with a Baker acceptance criterion (Robert & Casella, 2004, Section 7.8.1) that targets with temperature . Arguments such as (Bardenet et al. , 2014, Lemma 3.1, Proposition 3.2) could then help concluding that the distance between the kernels of both Markov chains is controlled, which would yield positive ergodicity results, in the line of (Bardenet et al. , 2014, Proposition 3.2). This reasoning shows again a close relation between subsampling and tempering, as in Section 6.1, with a clear link between the variance of the subsampled log likelihood ratios and the temperature.

Now, from a practical point of view, in simple applications such as logistic regression, is of the order of , which in turn should be of order if the MCMC proposal is a Gaussian random walk with covariance similar to that of , see Bardenet et al.  (2014). This means that has to be of order for the temperature to be of order , and this approach is thus bound to use subsamples per iteration! In conclusion, robustness to non-Gaussianity leads to requiring a fixed proportion of the whole dataset on average, even in the favourable case when one controls the first three moments of the subsampling noise and one swaps subsampling noise for the inherent MCMC acceptance noise.

6.4 Confidence samplers

In (Bardenet et al. , 2014), we proposed a controlled approximation of the acceptance decision (21). Indeed, let us fix and momentarily assume that was Lipschitz with known constant. Then, having observed the log likelihood ratio at some points , one could build a lower and an upper bound for the complete log likelihood ratio

 1nn∑i=1log[p(xi|θ′)p(xi|θ)],

simply by associating each with the nearest point among . These bounds could be refined by augmenting the set of observed log likelihoods ratios, until eventually one knows for sure whether (21) holds.

Now, concentration inequalities allow softer bounds and require less than this Lipschitz assumption. If one knows a bound for the range

 Cθ,θ′\defeqnmaxi=1∣∣ ∣∣log[p(xi|θ′)p(xi|θ)]∣∣ ∣∣, (25)

then concentration inequalities such as Hoeffding’s or Bernstein’s, yield confidence bounds such that

 P(∣∣ ∣∣1tt∑i=1log[p(x∗i|θ′)p(x∗i|θ)]−1nn∑i=1log[p(xi|θ′)p(xi|θ)]∣∣ ∣∣>ct(δ))≥1−δ, (26)

where the probability is taken over drawn uniformly from , with or without replacement. Borrowing from the bandit literature, we explain in (Bardenet et al. , 2014) how to leverage such confidence bounds to automatically select a subsample size such that the right MH acceptance decision is taken with a user-specified probability . Note that for our algorithm to bring any improvement over the ideal MH, the range (25) must be cheap to compute, i.e. cheaper than . This is the case for logistic regression, for example, but it is the major limitation of the approach in Bardenet et al.  (2014). We showed in (Bardenet et al. , 2014, Proposition 3.2) that if the ideal MH sampler is uniformly ergodic then the resulting algorithm inherits the uniform ergodicity of the ideal MH sampler, with a convergence speed that is within of that of the ideal MH. Importantly, we showed that our sampler then admits a limiting distribution, which is also within of . Uniform ergodicity is a very strong assumption and it would be worth extending these results to the geometrically ergodic scenario. There has recently been work in this direction (Alquier et al. , 2014; Pillai & Smith, 2014; Rudolf & Schweizer, 2015).

On the negative side, we demonstrated in (Bardenet et al. , 2014) that vanilla confidence samplers still require samples at each iteration at equilibrium, where the proportionality constant is the variance of the log likelihood ratio under subsampling. This statement relies on the leading term in being of order . In practice, the results of the vanilla confidence sampler on our running examples are shown in Figure 7. We set and we place ourselves in the favourable scenario where the algorithm has access to the actual range of each log likelihood ratio. The number of likelihood evaluations is estimated as follows: we take by default twice the detected value for the subsample size in general, but only once when the previous iteration required computing all likelihoods at the current state of the chain. Still, even in these favourable conditions, the algorithm basically requires essentially the whole dataset at each iteration.

Concentration inequalities are “worst-case” guarantees, and the theoretical results come at the price of a smaller reduction in the number of samples required. When the target is locally Gaussian, e.g. when Bernstein-von Mises yields a good approximation, there is potentially a lot to be gained, as empirically demonstrated by Korattikara et al.  (2014), for example. In the current paper, we propose in Section 7 a modified confidence sampler that can leverage concentration of the target to yield dramatic empirical gains while not sacrificing any theoretical guarantee of the confidence sampler. The basic tool is a cheap proxy for the log likelihood ratio that acts as a control variate in the concentration inequality (26). Using a 2nd order Taylor expansion centered at the maximum of the likelihood – obtained with a stochastic gradient descent for example – allows to replace many likelihood evaluations by the evaluation of this Taylor expansion. Figure 8 shows the results of this new confidence sampler with proxy on our running Gaussian and lognormal examples. Our algorithm outperforms all preceding methods, using almost no sample in the Gaussian case where it automatically detects that a quadratic form is enough to represent the log likelihood ratio. Finally, we demonstrate in Sections 7.2.3 and 8 that this new algorithm can require less than likelihood evaluations per iteration. Combined with the statements in Bardenet et al.  (2014) that each iteration is almost as efficient as the ideal MH, which is further supported by the match of the autocorrelation functions in Figures 8(c) and 8(d), this opens up big data horizons. We give full details on the confidence algorithm with proxy in Section 7.

7 An improved confidence sampler

In this section, we build upon the confidence sampler in (Bardenet et al. , 2014) by introducing likelihood proxies, which act as control variates for the individual likelihoods.

7.1 Introducing proxies in the confidence sampler

We start by recalling the pseudocode of the confidence sampler in (Bardenet et al. , 2014) in Figure 9, using sampling with replacement and a generic empirical concentration bound . In practice, one can think of the empirical Bernstein bound of Audibert et al.  (2009)

 ct(δ)=^σt,θ,θ′√2log(3/δ)t+6Cθ,θ′log(3/δ)t, (27)

where is the sample standard deviation of the log likelihood ratio

 {log[p(x∗i|θ′)p(x∗i|θ)],i=1,…,t},

and is their range, defined in (25). We emphasize that other choices of sampling procedure and concentration inequalities are valid, as long as they guarantee a concentration like (26). We refer the reader to (Bardenet et al. , 2014) for a proof of the correctness of the confidence sampler and implementation details.