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.


Introduction
Performing inference on tall datasets, that is datasets containing a large number n 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 O(n) 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 notebook 1 to the paper, available as supplementary material.

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.

Bayesian inference
and a parameter space Θ. We assume the data are conditionally independent with associated likelihood n i=1 p(x i |θ) given a parameter value θ and we denote (θ) the associated average log-likelihood We follow a Bayesian approach where one assigns a prior p(θ) to the unknown parameter, so that inference relies on the posterior distribution π(θ) = p(θ|x) ∝ γ(θ) p(θ)e n (θ) , 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.

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 (θ k ) k≥0 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 h where convergence is in distribution.
When the dataset is tall (n 1), 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.

Running examples
We will evaluate some of the described approaches on two illustrative simple running examples. We fit a one-dimensional normal distribution p(·|µ, σ) = N (·|µ, σ 2 ) to 10 5 i.i.d. points drawn according to X i ∼ N (0, 1) and lognormal observations X i ∼ log N (0, 1), respectively. The latter example illustrates a misspecification of the model. We assign a flat prior p(µ, log σ) ∝ 1. 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 1/ √ n and then adapted during the first 1 000 iterations so as to reach 50% acceptance. When applicable, we also display the number of likelihood evaluations per iteration, and compare it to the n evaluations required at each iteration by the MH algorithm.
In Figure 2, we illustrate the results of 10 000 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 50 000 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.

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.

Representing the posterior as a combination of batch posteriors
Assume data X are divided in B batches x 1 , . . . , x B . Relying on the equality Huang & Gelman (2005) propose to combine the batch posterior approximations using Gaussian approximations or importance sampling. Scott et al. (2013)  across batches, noting this is exact under Gaussian assumptions. Neiswanger et al. (2014) propose to run an MCMC chain on each batch x i targeting an artificial batch posterior 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 B. 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 B. In a tall data context, the number of batches is expected to grow with n to ensure that the size of each batch is less than O(n). Thus, the proposed bound is currently not informative for tall data.
As pointed out by Wang & Dunson (2013), if the supports of the π i are almost disjoint, then the product of their approximations will be a poor approximation to π. To improve the overlap between the approximations of the π i '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 π i 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 Θ B+1 , where the first B copies of θ are associated to the B batches, and the remaining copy is conditionally Gaussian around a weighted mean of the first B copies. Unfortunately, sampling from the posterior of this artificial model is difficult when one only has access to approximate samples of each π i .
Although it is not strictly speaking a Monte Carlo method, we note that Xu et al. (2014) andGelman 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.

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.

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

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 5. Similarly, it replaces γ(θ) in Step 5 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 σ 2 lim 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 com-putational complexity, the variance of the log-likelihood estimator should be kept around 1.0 when the ideal MH having access to the exact likelihood generates quasi-i.i.d samples from π; or set to around 3.0 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 loglikelihood 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 x * 1 , . . . , x * t from the dataset X with or without replacement, we obtain the following unbiased estimator of the log-likelihood n (θ) We denote byˆ (θ) the subsampling estimate of the average log-likelihood and denote by σ t (θ) 2 its variance. Obviously, exponentiating (7) does not provide an unbiased estimate of the likelihood e n (θ) . However, an interesting question is whether one can design a procedure which outputs an unbiased, almost-surely nonnegative estimate of e n (θ) using unbiased estimates of n (θ) such as nˆ (θ). Without making any further assumption about nˆ (θ), 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 a(θ) such that i (θ) > a(θ) for all i, 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 γ(θ)/p(θ), which is equivalent to definingγ(θ). For j ≥ 1, let be an unbiased estimator of n ( (θ) − a(θ)), where the x * i,j 's are drawn with replacement from X for each i, and are further independent across j. In (Jacob & Thiery, 2013, Section 3.1), N is an integer-valued random variable whose tails do not decrease too fast, in the sense that P(N ≥ k) ≥ C(1 + ) −k . To ease computations, we take N to be geometric with parameter /(1 + ). This corresponds to the lightest tails allowed by (Jacob & Thiery, 2013, Section 3 By (Rhee & Glynn, 2013, Theorem 1), Y is a non-negative unbiased estimator of the likelihood e n (θ) . As mentioned in Section 4.1, it is crucial, if we want to plugγ(θ) = Y × p(θ) in a pseudo-marginal algorithm, to control the variance of its logarithm. The variance of log Y is difficult to compute, so we use here the relative variance of Y as a proxy.
Proposition 4.1 Let θ ∈ Θ and Y be the almost surely non-negative estimator of e n (θ) defined in (9). Then its relative variance satisfies 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 Y not to increase exponentially with n, it is necessary that nσ t (θ) is of order 1. But σ t (θ) if of order t −1/2 , so that the batchsize t would have to be of order n 2 , which is impractical. It is also necessary that √ 1 + is of order 1 + n −1 to control the term in ( (θ) − a(θ)). This means that should be taken of order n −1 , but then the mean (1 + ) −1 of the geometric variable N will be of order n. This entails that the number of terms in the randomly truncated series (9) should be of order n, 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 Y 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 t, 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.

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 i (θ) ≥ b i (θ) is available for each i and θ. For simplicity, we further assume that b i (θ) = b(θ, x i ) only depends on i through x i . 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 i; 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 This distribution satisfies two important features: it admits π(θ) as a marginal distribution, and its pointwise evaluation only requires to evaluate i (θ) for those i's for which z i = 1. Note that evaluatingπ(θ, z) however requires to evaluate n i=1 exp(b i (θ)), and the bounds b i (θ) 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) 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π(z|θ). The expected cost of one conditional MH iteration on θ at equilibrium, that is the average number of indices i such that z i = 1, is O(n), 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 z 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 1.
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 Hence, the marginal distribution of z i under this extended model is given by where I θ exp(b(θ, x))dx. Using Bayes' theorem, we obtain accordingly An obvious unbiased estimator of the unnormalized posterior is thus given bŷ where each z i is drawn independently given θ from (12). Note that in the case of logistic regression, if b i (θ) is chosen to be the quadratic lower bound given in (MacLaurin & Adams, 2014), its integral I θ 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 nI θ , loosely speaking. Although the pseudo-marginal variant of Firefly we propose has the disadvantage of requiring the integrals I θ to be tractable, it comes with two advantages. First, the sampling of z 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, 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 1 will only be possible if most bounds b i (θ) are very tight, and the bigger n, the tighter the bounds have to be. These conditions will typically not be met when a fixed fraction of "outlier" x i '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 i (θ) 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 z are 1 at each iteration, and the number of likelihood evaluations per full joint iteration is thus roughly the fraction of points for which z i has been resampled. We chose the fraction of resampled points to be 10% here, and initialized z to have 10% 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 z i '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 10%, 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 z can easily take θ out of the posterior mode. To be fair, we thus discarded the first 1 000 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 z i draws the µ-component of the chain towards the few large values of (x i − µ) 2 which are bright. Since the bright points are rarely updated, the chain mixes very slowly.

Other exact approaches
Other exact approaches have been proposed, which do not rely on pseudo-marginal MH. Welling & Teh (2011) proposed an algorithm based on stochastic gradient Langevin dynamics (SGLD). This is an iterative algorithm which at iteration k + 1 uses the following update rule

Forgetting about acceptance: stochastic approximation approaches
is an unbiased estimate of the score computed at each iteration using a random subsample x * i,k 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 (2011) is that if one suppresses the MH acceptance step, computes an unbiased estimate of the score but introduces a sequence of stepsizes ( k ) that decreases to zero at the right rate, then is an approximation to π. The algorithm has been analyzed recently in , where it has been established that it provides indeed a consistent estimate of the target. Additionally,  a central limit theorem holds with convergence rate N −1/3 iter , which is slower than the traditional Monte Carlo rate N −1/2 iter . 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 5.1. The stepsize k is chosen proportional to k −1/3 , following the recommendations of Teh et al. (2014). We show the results of two choices for the subsample size t: 10% and 1% of the data, with respectively 10 000 and 100 000 iterations, so that both runs amount to the same 10% 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. Banterle et al. (2015) remarked that if we decompose the acceptance ratio in a product of positive functions

Delayed acceptance
then the MH-like algorithm that accepts the move from θ to θ with probability B i=1 min ρ i (θ, θ ), 1 still admits π as invariant distribution. In practice, in the case of tall data, we can divide the dataset into B batches and use for example This allows us to reject candidate θ without having to compute the full likelihoods and the calculations of ρ i (θ, θ ) can be done in parallel. However, as remarked by Banterle et al. (2015), the resulting Markov chain has a larger asymptotic variance σ 2 lim 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 25%. 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.

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 n (θ) based on random samples x * 1 , . . . , x * t from the dataset X ; see (7). Similarly, one can also easily obtain an unbiased estimator of the average log likelihood ratio [ (θ ) − (θ)] 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 π.

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 n independent Bernoulli variables and let be an unbiased estimator of the average log likelihood (θ), where z i ∼ B(1, λ) i.i.d. Now one could think of plugging estimatesγ(θ) = p(θ)e nˆ (θ) in Steps 2 and 3 of MH in Figure 1. However, asγ(θ) is not an unbiased estimator of γ(θ), this algorithm samples from a target distribution which is proportional to p(θ)Ee nˆ (θ) = γ(θ), where the expectation is w.r.t the distributions of the Bernoulli random variables {z i }. Now where I r denotes a set of r distinct indices in {1, . . . , n}, x Ir = {x i ; i ∈ I r }, and with the convention p(x ∅ |θ) = 1. Each subsampled likelihood contributes to the target, exponentiated to the power 1/λ, resulting in a nontrivial mixture of rescaled data likelihood terms. To further simplify, assume p(x Ir |θ) ≈ p r (θ) for each set of indices I r , that is, the variance of the likelihood under subsampling is small, then where B(n, λ) denotes the binomial distribution with parameters n and λ. Noticing that p r (θ) is roughly exponentially decreasing in r, the values or r 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 nλ, and the power 1/λ makes this term of the same scale as p(x 1 , . . . , x n |θ). 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 instead of (17) one ends up with Again, all subsampled likelihoods contribute, but without being further exponentiated. Still, the result is a much broadened target, as values of r that are larger than nλ 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, 1 or 3 depending on hypotheses, in order for pseudomarginal MH to be efficient. This means that λ should be such that is of order 1 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), can be of order 1 if λ ∼ n −1 . 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.

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 i (θ ) − i (θ) are heavy-tailed. Missing log likelihoods in the tails will lead to erroneous decisions, and uncontrolled results. 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 p(θ)Ee nˆ (θ)−b(θ) , whereb(θ) is an estimate of the bias b(θ) satisfying Ee nˆ (θ) = e n (θ)+b(θ) . 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 i being proportional to | i (θ)|. 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 i (θ). 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 10% 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 x → p(x|θ) 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.

Adaptive subsampling with T-tests
Still assuming the noise of the log likelihood is Gaussian, given a drawn θ ∈ Θ, one can try to adaptively choose the size of the subsample {x * 1 , . . . , x * t } 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 log α(θ, θ ) > u, or equivalently with u ∼ U [0,1] 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 10 000 iterations of Austerity MH on our two running examples in Figure 5. The parameters are = 0.05, corresponding to the p-value threshold in the aforementioned T-test, and an initial subsample size of 100 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 100, 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 X . Finally, the reductions in the number of samples needed per iteration are quite interesting: half of the iterations require less than 4% 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 n 1. 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.

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 O(n) 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 x * 1 , . . . , x * t be drawn independently with replacement from X . Let Λ * t (θ, θ ) 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 u ∼ U[0, 1]. Why, then, not rely on the subsampling noise to guarantee exploration, and accept a move if and only if 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 We do not assume that the Λ * t (θ, θ )'s are Gaussianly distributed, but we make the parametric assumption that the second and third absolute moments σ 2 and ρ of − log p(x * i |θ ) + log p(x * i |θ) are known and independent of θ, θ . Applying the Berry-Esseen inequality (van der Vaart & Wellner, 1996) for any u ∈ R, where Φ is the cdf of a N (0, 1) variable, and is the average log likelihood ratio. When u = 0, (23) yields Now let C, λ > 0 be such that for any x ∈ R, Hence, the acceptance probability of an algorithm that would accept the move from θ to θ by checking whether Λ * (θ, θ ) > 0 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 β = λ √ t nσ . 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 O p (n −1/2 ) if the MCMC proposal is a Gaussian random walk with covariance similar to that of π, see Bardenet et al. (2014). This means that t has to be of order n for the temperature β to be of order 1, and this approach is thus bound to use O(n) 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.

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 x → log[p(x|θ )/p(x|θ)] was Lipschitz with known constant. Then, having observed the log likelihood ratio at some points {x * i , i = 1, . . . , t} ⊂ X , one could build a lower and an upper bound for the complete log likelihood ratio 1 n simply by associating each x i with the nearest point among {x * 1 , . . . , x * t }. 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 then concentration inequalities such as Hoeffding's or Bernstein's, yield confidence bounds c t (δ) such that where the probability is taken over x * 1 , . . . , x * t drawn uniformly from X , 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 T such that the right MH acceptance decision is taken with a user-specified probability 1 − δ. 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 O(n). 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 O(δ) of that of the ideal MH. Importantly, we showed that our sampler then admits a limiting distribution, which is also within O(δ) 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 O(n) 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 c t (δ) being of order t −1/2 . In practice, the results of the vanilla confidence sampler on our running examples are shown in Figure 7. We set δ = 0.1 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 T for the subsample size in general, but only once when the previous iteration required computing all n 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 O(n) 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.

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.

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 c t (δ). In practice, one can think of the empirical Bernstein bound of Audibert et al. (2009) whereσ t,θ,θ is the sample standard deviation of the log likelihood ratio , i = 1, . . . , t , and C θ,θ 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.
The bottleneck for the performance of the confidence sampler was identified in (Bardenet et al. , 2014) as the expectation w.r.t. π(θ)q(θ |θ) of the variance of the log likelihood ratio log p(x|θ )/p(x|θ) w.r.t. to the empirical distribution of the observations. We now propose a control variate technique inspired from the Firefly MH of MacLaurin & Adams (2014) to lower this variance down when an accurate and cheap proxy of the log likelihood is known.
We require a proxy for the log likelihood ratio that may decrease the variance of the log likelihood ratio or its range. More precisely, let ℘ i (θ, θ ) be such that for any θ, θ ∈ Θ, can be bounded uniformly in 1 ≤ i ≤ n, and the bound is cheap to compute.
We now simply remark that the acceptance decision (21) in MH is equivalent to checking whether Building the confidence sampler on (28) leads to the same pseudocode as in Figure 9, except that Step 14 is replaced by the condition in Step 19 is replaced by

and the condition in
Step 21 becomes For completeness, we restate here in Proposition 7.1 that the vanilla confidence sampler inherits the uniform ergodicity of the underlying MH sampler, that its target is within O(δ) of π, and that the difference in speed of convergence is also controlled by δ. Let P be the underlying MH kernel, andP the kernel of the confidence sampler described in this section.
Even in the presence of proxies, the proofs of (Bardenet et al. , 2014, Lemma 3.1, Proposition 3.2) apply with straightforward modifications, so that we can extend Proposition 7.1 to the proxy case. The major advantage of this new algorithm is that the sample standard deviation σ t,θ,θ and range C θ,θ in the concentration inequality (27) are replaced by those of If ℘ i (θ, θ ) is a good proxy for the log likelihood ratio, one can thus expect significantly more accurate confidence bounds, leading in turn to reduction in the number of samples used.

An example proxy: Taylor expansions
In general, the choice of proxy ℘ will be problem-dependent, and the availability of a good proxy at all is a strong assumption, although not as strong as our previous requirement in Bardenet et al. (2014) that the range (25) can be computed cheaply, which basically corresponds to ℘ i (θ, θ ) being identically zero for all i. Indeed, we show in this section that all models that possess up to third derivatives can typically be tackled using Taylor expansions as proxies. In Section 8, we detail the case of logistic regression and gamma linear regression.

Taylor expansions
We expand i around some reference value θ to obtain an estimatê where g i, and H i, are respectively the gradient and the Hessian of i at θ . The choice of θ is deferred to Section 7.2.2. Let us now define Indeed, the following holds Finally, assuming ∂ i ∂θ (j) ∂θ (k) ∂θ (l) θ can be bounded uniformly in i, j, k, l, the absolute difference i (θ ) − i (θ) − ℘ i (θ, θ ) can be bounded using the Taylor-Lagrange inequality. To conclude, all conditions of Section 7.1 are satisfied by the proxy ℘ i .

Drop proxies along the way
When the mass of the posterior is concentrated around the maximum likelihood estimator θ MLE , a single proxy -say a Taylor proxy centered at θ = θ MLE -can represent the target quite accurately. This is the proxy we used in the running examples of Section 6, see Figure 8. When the posterior does not concentrate, or the proposal is not local enough, such a proxy will be inaccurate, potentially resulting in insufficient subsampling gains. There are various tricks that can be applied. One can either precompute proxies across Θ if one has an idea where the modes of π are, and then use the closest proxy to the current state of the chain at each iteration. Alternately, if one agrees to look at the whole dataset every α iterations, we can drop proxies along the way, i.e. set θ to the current state of the chain every α MH iterations. The whole dataset needs to be browsed at each change of the reference point θ , since there is typically some preprocessing to do in order to compute later bounds. In the case of 2nd order Taylor expansions, for example, one has to compute the full gradient, Hessian, and any other quantity needed to bound the third derivatives. What the user should aim at is to sacrifice a proportion α of the budget of the ideal MH to make the remaining iterations cheaper. The proof of Proposition 7.1 easily generalizes to the case of proxies dropped every constant number of iterations. We demonstrate the empirical performance of such an approach in Sections 8.1.3 and 8.2.2.

A heuristic on the subsampling gain
In (Bardenet et al. , 2014), we presented a heuristic that showed the original confidence sampler required O(n) likelihood evaluations per iteration. At the time, it seemed every attempt at marrying subsampling and MH was fundamentally O(n). We first repeat here the heuristic from (Bardenet et al. , 2014), before arguing that the contributions of this paper can lower this budget to o(n), even O(1) up to polylogarithmic factors in very favourable conditions. Assuming a symmetric proposal and a flat prior, the stopping rule of the while loop in the original confidence sampler in Figure 9 is met whenever is of the same order as the confidence bound c t (δ), that is, when c t (δ) is of order 1/n. We consider the Bernstein bound in (27) and we assume the range C t,θ,θ grows with n strictly slower than √ n. The latter assumption is realistic: C t,θ,θ is often dominated by some power of max x i ∞ , and if x 1 , . . . , x n are drawn i.i.d. from a subgaussian distribution, then -Bianchi & Lugosi, 2006, Lemma A.13). The leading term of the Bernstein bound is proportional toσ t,θ,θ / √ t. In simple models such as logistic regression,σ t,θ,θ is proportional to θ − θ .
Assuming n is large enough that standard asymptotics apply and the target is approximately Gaussian, the results of Roberts & Rosenthal (2001) lead to choose the covariance matrix of the proposal such that θ − θ is of order n −1/2 . Summing up, we exit the while loop when which leads to t ∼ n. Now consider the confidence sampler with second-order Taylor proxies introduced in Section 7.2.1.σ t,θ,θ and C t,θ,θ now correspond to the standard deviation and range of Now let us assume the third-order derivatives at the reference point θ can be bounded, say by some constant times max i X i 3 ∞ as will be the case for the exponential family models of Section 8. Thenσ t,θ,θ and C t,θ,θ are dominated by But θ − θ is of order n −1/2 if standard asymptotics (van der Vaart, 2000) yield good approximations and θ is set to the maximum of the posterior. Alternatively, if one has implemented the strategy of dropping proxies regularly, then θ − θ should be of order n −1/2 since we assume the covariance matrix of the proposal distribution is of order 1/n. Again assuming that max i X i 3 ∞ grows, say, like ρ(n) = o(n 1/3 ), we now exit the while loop when Thus, when the target is approximately Gaussian and the chain is in the mode, the cost in likelihood evaluations per iteration of the confidence sampler with proxy is likely to be o(n).
The actual order of convergence depends on the rate of growth of the bounds on the third derivatives. For example, in the case of independent Gaussian data and still assuming (32), we have t = O(1) up to polylogarithmic factors.

Experiments
As a proof of concept, all experiments in this section avoid loading the dataset or proxy-related quantities into memory by building, maintaining and querying from a disk-based database using SQLite 2 .

Logistic regression 8.1.1 A Taylor proxy for logistic regression
In logistic regression, the likelihood is defined by and the label t i is in {−1, +1}. We can use the Taylor expansion proxy of Section 7.2.1, using

A toy example that requires O(1) likelihood evaluations
In this section, we consider the simple two-dimensional logistic regression dataset in (Bardenet et al. , 2014, Section 4.2.2), where the features within each class are drawn from a Gaussian. The dataset is depicted in Figure 10(a). We consider subsets of the dataset with increasing size log 10 n ∈ {3, 4, 5, 6, 7}, run a confidence MH chain for each n, started at the MAP, with δ = 0.1 and a single proxy around the MAP. We report the numbers of likelihood evaluations L at each iteration in Figure 10(b). The fraction of likelihood evaluations compared to MH roughly decreases by a factor 10 when the size of the dataset is multiplied by 10: the number of likelihood evaluations is constant for n large enough. In other words, 1 000 random data points at each iteration are enough to get within O(δ) of the actual posterior, the rest of the dataset appears to be superfluous. There is a saturation phenomenon. By relaxing the goal of sampling from π into sampling from a controlled approximation, we can break the O(n) barrier and in this particular example reach a cost per iteration of O(1).

The covtype dataset
We consider the dataset covtype.binary 3 described in Collobert et al. (2002). The dataset consists of 581,012 points, of which we pick n = 400, 000 as a training set, following the maximum training size in Collobert et al. (2002). The original dimension of the problem is 54, with the first 10 attributes being quantitative. To illustrate our point without requiring a more complex sampler than MH, we only consider the 10 quantitative attributes. We use the preprocessing and Cauchy prior recommended by Gelman et al. (2008).
We run 5 independent chains for 10 000 iterations, dropping proxies every 10 iterations as explained in Section 7.2.2. We obtain a Gelman-Rubin statistic of 1.01 (Robert & Casella, 2004, Section 12.3.4), which suggests the between-chain variance is low enough that we can stop sampling.
We estimate the number of likelihood evaluations L k at MH iteration k as follows. First, note that -dropping proxies or not-on a regular iteration where the proxy is not necessarily recomputed, L k can take values up to 2n, unlike MH, which can store the evaluation of the likelihood at the current state of the chain from one iteration to the next, and thus only requires n likelihood evaluations per iteration. Second, at an iteration where the proxy is recomputed, the whole data has to be read anyway, so that we choose here to perform a normal MH iteration. This requires the maximum 2n likelihood evaluations, Assuming the cost of the likelihood evaluation is the bottleneck, we neglect here the additional cost of computing the proxy itself, and only report L k = 2n when the proxy is recomputed. Third, whenever we compute the full likelihood at a state of the chain, we store it until the chain leaves that state, similarly to any implementation of MH. Thus, at an iteration that follows a full read of the data, i.e. L k−1 = 2n, we only count the likelihood evaluations of the proposed state.
We summarize the results in Figure 11. All runs use on average 27 to 42% of n likelihood evaluations per iteration. Since we compute the proxy every α = 10 iterations, there is a necessary 2 × 10 = 20% of n that is due to recomputing the proxy. We manually assessed the value of α, and recomputing the proxy less often increases the average number of likelihood evaluations (not shown). Thanks to these forced 20%, the rest of the iterations are considerably cheaper than n, since 50% of the iterations require less than 5% of the dataset, as shown in Figure 11(b). Relatedly, although subsampling implies a forced 2n likelihood evaluations to start and thus shows an initial delay in Figure 11(a), it quickly catches up and converges faster. The gains are two-or threefold, which is of limited overall practical interest, but we know from Section 7.2.3 and Figure 10(b) that increasing n will also improve the gain.

Gamma linear regression 8.2.1 A Taylor proxy for gamma regression
In gamma regression, the nonnegative response y is assumed to be gamma-distributed where Γ(κ, s) is the gamma distribution with shape parameter κ and scale parameters s. Assuming κ is known, the log likelihood is thus given by log p(y|x, θ) ∝ −κye −θ T x − κθ T x up to an additive constant, so that ∇ log p(y|x, θ) = κ ye −x T θ − 1 x, Hess(log p(y|x, θ)) = −κye −x T θ xx T and ∂ ∂θ (j) ∂θ (k) ∂θ (l) log p(y|x, θ) = κye −x T θ x (j) x (k) x (l) . Furthermore, we can bound The Taylor proxies of Section 7.2.1 can thus be applied.

The covtype dataset
As an application, we consider the covtype dataset again and regress the nonnegative feature "horizontal distance to nearest wildfire ignition" onto the other quantitative features. We run 5 independent chains for 10 000 iterations, dropping proxies every 10 iterations as explained in Section 7.2.2. We obtain a Gelman-Rubin statistic of 1.001, which again suggests we can stop sampling. We estimate the evaluation budget as in Section 8.1.3. We summarize the results in Figure 12.
All runs use on average 33 to 54% of n likelihood evaluations per iteration, from which 2 × 10 = 20% are due to recomputing the proxy every 10 iterations. Recomputing the proxy less often increases the average number of likelihood evaluations (not shown). Thanks to these forced 20% the rest of the iterations are considerably cheaper than n, since, as in Section 8.1.3, 50% of the iterations require less than 10% of the dataset. Relatedly, and similarly to the logistic regression task in Section 8.1.3 subsampling converges two or three times faster in this example. Again, this is a proof of concept that subsampling works, and we know from Section 7.2.3 and Figure 10(b) that increasing n will also improve the gain.

Discussion
We have reviewed recent advances in applying MCMC to tall datasets. Divide-and-conquer approaches have yet to solve the recombination problem, i.e. how to obtain a meaningful distribution in a stable manner from the output of individual chains on a growing number of smaller datasets. Subsampling approaches face different issues, namely that of approaching the right target at a known speed, and of keeping the overall budget in terms of likelihood evaluations per iteration low.
In this paper, we have proposed an original subsampling approach. We have showed that under strong ergodicity assumptions on the original MH sampler, our algorithm samples from a controlled approximation of the posterior target. While these strong assumptions are rarely satisfied in practice, our experiments suggest that our results extend to more general scenarios.
In terms of scaling, the introduced methodology is even able to lower the natural cost of O(n) subsamples per iteration to as low as O(1) in favourable scenarios. However, we have yet only observed these dramatic gains in contexts where the Bernstein-von Mises approximation is already excellent. On the positive side, our algorithm improves on other proposed subsampling approaches in this context. On the negative side, computing the Bernstein-von Mises approximation for regular models can be typically achieved in only a couple of passes over the data, using for example stochastic gradient to compute the maximum likelihood estimator, and the observed information matrix at this point to estimate the Hessian.

Appendix B: proof of Proposition 4.2
We write