On Markov chain Monte Carlo methods for tall data

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 ${\cal 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 The 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.

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.

We follow a Bayesian approach where one assigns a prior $p(\theta)$ to the unknown parameter, so that inference relies on the posterior distribution

where $\gamma$ denotes an unnormalized version of $\pi$. In most applications, $\pi$ 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 The Metropolis-Hastings algorithm

A standard approach to sample approximately from $\pi(\theta)$ 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 $(\theta_{k})_{k\geq 0}$ of invariant distribution $\pi$. 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$

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

When the dataset is tall ($n\gg 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.

3 Running examples

We will evaluate some of the described approaches on two illustrative simple running examples. We fit a one-dimensional normal distribution $p(\cdot|\mu,\sigma)={\cal N}(\cdot|\mu,\sigma^{2})$ to $10^{5}$ i.i.d. points drawn according to $X_{i}\sim{\cal N}(0,1)$ and lognormal observations $X_{i}\sim\log{\cal N}(0,1)$, respectively. The latter example illustrates a misspecification of the model. We assign a flat prior $p(\mu,\log\sigma)\propto 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/\sqrt{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.

Assume data ${\cal X}$ are divided in $B$ batches ${\bf x}_{1},\dots,{\bf 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 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 ${\bf 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 $\pi$ 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 ${\cal 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 $\pi_{i}$ are almost disjoint, then the product of their approximations will be a poor approximation to $\pi$. To improve the overlap between the approximations of the $\pi_{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 $\pi_{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 $\Theta^{B+1}$, where the first $B$ copies of $\theta$ 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 $\pi_{i}$.

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.

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.

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.

Assume that instead of being able to evaluate $\gamma(\theta)$, we have access to an unbiased, almost-surely non-negative estimator $\hat{\gamma}(\theta)$ of the unnormalized target $\gamma(\theta)$. Pseudo-marginal MH substitutes a realization of $\hat{\gamma}(\theta^{\prime})$ to $\gamma(\theta^{\prime})$ in Step 1. Similarly, it replaces $\gamma(\theta)$ in Step 1 by the realization of $\hat{\gamma}(\theta)$ that was computed when the parameter value $\theta$ 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 $\gamma$ comes at a price: first, the asymptotic variance $\sigma_{\text{lim}}^{2}$ 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 $\hat{\gamma}(\vartheta)$ is large for some value $\vartheta\in\Theta$, then an MH move to $\vartheta$ might be accepted while $\hat{\gamma}(\vartheta)$ largely overestimates $\gamma(\vartheta)$. In that case, it is difficult for the chain to leave $\vartheta$, 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 $1.0$ when the ideal MH having access to the exact likelihood generates quasi-i.i.d samples from $\pi$; 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 $\hat{\gamma}(\theta)$, or, if not available, of their relative variance.

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

We apply (Rhee & Glynn, 2013, Theorem 1) to build an unbiased non-negative estimator of $\gamma(\theta)/p(\theta)$, which is equivalent to defining $\hat{\gamma}(\theta)$. For $j\geq 1$, let

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 ${\epsilon}$, 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.

3 Building γ^​(θ)\hat{\gamma}(\theta) with auxiliary variables

where $I_{\theta}{\;\triangleq\;}\int\exp(b(\theta,x))dx$. Using Bayes’ theorem, we obtain accordingly

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

where each $z_{i}$ is drawn independently given $\theta$ from (12). Note that in the case of logistic regression, if $b_{i}(\theta)$ is chosen to be the quadratic lower bound given in (MacLaurin & Adams, 2014), its integral $I_{\theta}$ 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_{\theta}$, loosely speaking.

Although the pseudo-marginal variant of Firefly we propose has the disadvantage of requiring the integrals $I_{\theta}$ 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.

Let $\theta\in\Theta$. 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}(\theta)$ 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.

As expected, the algorithm behaves erratically in the lognormal case, as failure to attempt a flip of each $z_{i}$ draws the $\mu$-component of the chain towards the few large values of $(x_{i}-\mu)^{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

$({\epsilon}_{k})$ is a sequence of time steps, $(\eta_{k})$ are independent ${\cal N}(0,I_{d})$ vectors and

is an unbiased estimate of the score computed at each iteration using a random subsample $\left\{x_{i,k}^{*}\right\}$ of the observations. This approach is reminiscent of the Metropolis-adjusted Langevin algorithm (MALA; Robert & Casella, 2004), where the proposal given by

is used in an MH acceptance step, where ${\epsilon}\sim{\cal N}(0,I_{d})$. 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 $({\epsilon}_{k})$ that decreases to zero at the right rate, then

is an approximation to $\pi$. 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 $N_{\text{iter}}^{-1/3}$, which is slower than the traditional Monte Carlo rate $N_{\text{iter}}^{-1/2}$. 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 ${\epsilon}_{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.

2 Delayed acceptance

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

then the MH-like algorithm that accepts the move from $\theta$ to $\theta^{\prime}$ with probability

still admits $\pi$ 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 $\theta^{\prime}$ without having to compute the full likelihoods and the calculations of $\rho_{i}(\theta,\theta^{\prime})$ can be done in parallel. However, as remarked by Banterle et al. , 2015, the resulting Markov chain has a larger asymptotic variance $\sigma_{\text{lim}}^{2}$ 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

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

The first approach one could try is to only use a random fixed proportion of data points to estimate $\pi$ at any newly drawn $\theta$. 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 $\theta$ drawn in an MH run, we draw $n$ independent Bernoulli variables and let

where $I_{r}$ denotes a set of $r$ distinct indices in $\{1,\dots,n\}$, $x_{I_{r}}=\{x_{i};i\in I_{r}\}$, and with the convention $p(x_{\emptyset}|\theta)=1$. Each subsampled likelihood contributes to the target, exponentiated to the power $1/\lambda$, resulting in a nontrivial mixture of rescaled data likelihood terms. To further simplify, assume $p(x_{I_{r}}|\theta)\approx p_{r}(\theta)$ for each set of indices $I_{r}$, that is, the variance of the likelihood under subsampling is small, then

where $B(n,\lambda)$ denotes the binomial distribution with parameters $n$ and $\lambda$. Noticing that $p_{r}(\theta)$ 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\lambda$, and the power $1/\lambda$ makes this term of the same scale as $p(x_{1},\dots,x_{n}|\theta)$. 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

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\lambda$ 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 pseudo-marginal MH to be efficient. This means that $\lambda$ should be such that

is of order $1$ in the case of (17). This entails that $\lambda$ 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 $\lambda\sim 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 $\pi$.

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.

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.

2.2 Adaptive subsampling with T-tests

Still assuming the noise of the log likelihood is Gaussian, given a drawn $\theta\in\Theta$, one can try to adaptively choose the size of the subsample $\{x_{1}^{*},\dots,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\alpha(\theta,\theta^{\prime})>u$, or equivalently

with $u\sim{\cal U}_{}$ 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 $(\theta,\theta^{\prime})$, 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 ${\epsilon}=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 ${\cal 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\gg 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 $\theta,\theta^{\prime}\in\Theta$, 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 $\Theta\times\Theta$ to make sure the CLT assumptions are realistic. Note that even then, there is no guarantee the above algorithms have $\pi$ for target, if any.

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 ${\cal 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 $\theta,\theta^{\prime}\in\Theta$, and let $x_{1}^{*},\dots,x_{t}^{*}$ be drawn independently with replacement from ${\cal X}$. Let $\Lambda_{t}^{*}(\theta,\theta^{\prime})$ 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\sim{\cal U}$. 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 $\Lambda_{t}^{*}(\theta,\theta^{\prime})$’s are Gaussianly distributed, but we make the parametric assumption that the second and third absolute moments $\sigma^{2}$ and $\rho$ of $-\log p(x_{i}^{*}|\theta^{\prime})+\log p(x_{i}^{*}|\theta)$ are known and independent of $\theta,\theta^{\prime}$. Applying the Berry-Esseen inequality (van der Vaart & Wellner, 1996) to the variables $-\log p(x_{i}^{*}|\theta^{\prime})+\log p(x_{i}^{*}|\theta)$ yields

is the average log likelihood ratio. When $u=0$, (23) yields

Bowling et al. , 2009 for instance, empirically found $C=0.0095$ and $\lambda=1.702$. Combining this bound with (24), we obtain

Hence, the acceptance probability of an algorithm that would accept the move from $\theta$ to $\theta^{\prime}$ by checking whether $\Lambda^{*}(\theta,\theta^{\prime})>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 $\pi^{\beta}$ with temperature $\beta=\frac{\lambda\sqrt{t}}{n\sigma}$. 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, $\sigma$ is of the order of $\|\theta-\theta^{\prime}\|$, which in turn should be of order ${\cal O}_{p}(n^{-1/2})$ if the MCMC proposal is a Gaussian random walk with covariance similar to that of $\pi$, see Bardenet et al. , 2014. This means that $t$ has to be of order $n$ for the temperature $\beta$ to be of order $1$, and this approach is thus bound to use ${\cal 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.

4 Confidence samplers

In (Bardenet et al. , 2014), we proposed a controlled approximation of the acceptance decision (21). Indeed, let us fix $\theta,\theta^{\prime}$ and momentarily assume that $x\mapsto\log[p(x|\theta^{\prime})/p(x|\theta)]$ was Lipschitz with known constant. Then, having observed the log likelihood ratio at some points $\{x_{i}^{*},i=1,\dots,t\}\subset{\cal X}$, one could build a lower and an upper bound for the complete log likelihood ratio

simply by associating each $x_{i}$ with the nearest point among $\{x_{1}^{*},\dots,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}(\delta)$ such that

where the probability is taken over $x_{1}^{*},\dots,x_{t}^{*}$ drawn uniformly from ${\cal 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-\delta$. 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 ${\cal 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 ${\cal O}(\delta)$ of that of the ideal MH. Importantly, we showed that our sampler then admits a limiting distribution, which is also within ${\cal O}(\delta)$ of $\pi$. 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 ${\cal 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}(\delta)$ 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 $\delta=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 ${\cal 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.

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}(\delta)$. In practice, one can think of the empirical Bernstein bound of Audibert et al. , 2009

where $\hat{\sigma}_{t,\theta,\theta^{\prime}}$ is the sample standard deviation of the log likelihood ratio

and $C_{\theta,\theta^{\prime}}$ 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. $\pi(\theta)q(\theta^{\prime}|\theta)$ of the variance of the log likelihood ratio $\log p(x|\theta^{\prime})/p(x|\theta)$ 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 $\wp_{i}(\theta,\theta^{\prime})$ be such that for any $\theta,\theta^{\prime}\in\Theta$,

$\sum_{i=1}^{n}\wp_{i}(\theta,\theta^{\prime})$ can be computed cheaply.

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 9 is replaced by

Let $P$ be uniformly geometrically ergodic, i.e., there exists an integer $m$, a probability measure $\nu$ on $\left(\Theta,\mathcal{B}\left(\Theta\right)\right)$ and $0\leq\rho<1$ such that for all $\theta\in\Theta$, $P^{m}(\theta,\cdot)\geq\left(1-\rho\right)\nu(\cdot)$ . Hence there exists $A<\infty$ such that

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 $\hat{\sigma}_{t,\theta,\theta^{\prime}}$ and range $C_{\theta,\theta^{\prime}}$ in the concentration inequality (27) are replaced by those of

If $\wp_{i}(\theta,\theta^{\prime})$ 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.

2 An example proxy: Taylor expansions

In general, the choice of proxy $\wp$ 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 $\wp_{i}(\theta,\theta^{\prime})$ 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.

2.2 Drop proxies along the way

When the mass of the posterior is concentrated around the maximum likelihood estimator $\theta_{\text{MLE}}$, a single proxy – say a Taylor proxy centered at $\theta_{\star}=\theta_{\text{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 $\Theta$ if one has an idea where the modes of $\pi$ 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 $\alpha$ iterations, we can drop proxies along the way, i.e. set $\theta_{\star}$ to the current state of the chain every $\alpha$ MH iterations. The whole dataset needs to be browsed at each change of the reference point $\theta_{\star}$, 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 $\alpha$ 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.

2.3 A heuristic on the subsampling gain

In (Bardenet et al. , 2014), we presented a heuristic that showed the original confidence sampler required ${\cal O}(n)$ likelihood evaluations per iteration. At the time, it seemed every attempt at marrying subsampling and MH was fundamentally ${\cal 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 ${\cal 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

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 $\|\theta-\theta^{\prime}\|$ is of order $n^{-1/2}$. Summing up, we exit the $\While$ loop when

Now consider the confidence sampler with second-order Taylor proxies introduced in Section 7.2.1. $\hat{\sigma}_{t,\theta,\theta^{\prime}}$ and $C_{t,\theta,\theta^{\prime}}$ now correspond to the standard deviation and range of

Now let us assume the third-order derivatives at the reference point $\theta_{\star}$ can be bounded, say by some constant times $\max_{i}\|X_{i}\|_{\infty}^{3}$ as will be the case for the exponential family models of Section 8. Then $\hat{\sigma}_{t,\theta,\theta^{\prime}}$ and $C_{t,\theta,\theta^{\prime}}$ are dominated by

But $\|\theta-\theta_{\star}\|$ is of order $n^{-1/2}$ if standard asymptotics (van der Vaart, 2000) yield good approximations and $\theta_{\star}$ is set to the maximum of the posterior. Alternatively, if one has implemented the strategy of dropping proxies regularly, then $\|\theta-\theta_{\star}\|$ 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}\|_{\infty}^{3}$ grows, say, like $\rho(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={\cal 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 http://www.sqlite.org/.

and the label $t_{i}$ is in $\{-1,+1\}$. We can use the Taylor expansion proxy of Section 7.2.1, using

1.2 A toy example that requires 𝒪⁡(1){\cal 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\in\{3,4,5,6,7\}$, run a confidence MH chain for each $n$, started at the MAP, with $\delta=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 ${\cal O}(\delta)$ 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 $\pi$ into sampling from a controlled approximation, we can break the ${\cal O}(n)$ barrier and in this particular example reach a cost per iteration of ${\cal O}(1)$.

1.3 The covtype dataset

We consider the dataset covtype.binary available at http://www.csie.ntu.edu.tw/ cjlin/libsvmtools/datasets/binary.html 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 $\alpha=10$ iterations, there is a necessary $2\times 10=20\%$ of $n$ that is due to recomputing the proxy. We manually assessed the value of $\alpha$, 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.

2 Gamma linear regression

In gamma regression, the nonnegative response $y$ is assumed to be gamma-distributed

where $\Gamma(\kappa,s)$ is the gamma distribution with shape parameter $\kappa$ and scale parameters $s$. Assuming $\kappa$ is known, the log likelihood is thus given by

The Taylor proxies of Section 7.2.1 can thus be applied.

2.2 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\times 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 ${\cal O}(n)$ subsamples per iteration to as low as ${\cal 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.

Further work should thus now focus on demonstrating the applicability of subsampling approaches to cases where it is either difficult to compute Bernstein-von Mises even if it is a good approximation (Chernozhukov & Hong, 2003), or – more importantly – cases where $n$ is not big enough that Bernstein-von Mises yields a good approximation.

The authors acknowledge Louis Aslett, Nando de Freitas, Pierre Jacob, François Septier, Matti Vihola, and Sebastian Vollmer for their comments and discussions on this paper and topic.

Appendix A: proof of Proposition 4.1

From (Rhee & Glynn, 2013, Theorem 1), the second moment of $Y$ is

with the convention $S_{-1}=0$ and $S_{0}=e^{na(\theta)}$. We note that

Now, since $k!k!\leq 4^{-k}(2k+1)!$, letting

Appendix B: proof of Proposition 4.2

References