Control Variates for Stochastic Gradient MCMC

TL;DR

This paper introduces control variates to stochastic gradient MCMC methods, significantly reducing variance and computational cost, making the algorithms more scalable for large datasets under certain conditions.

Contribution

It proposes a novel variance reduction technique for stochastic gradient MCMC using control variates, improving efficiency and scalability of the methods.

Findings

Variance reduction leads to dataset-size independent computational cost.

Zero variance control variates can be applied for free to improve inference.

Analysis under log-concavity assumptions shows theoretical efficiency gains.

Abstract

It is well known that Markov chain Monte Carlo (MCMC) methods scale poorly with dataset size. A popular class of methods for solving this issue is stochastic gradient MCMC. These methods use a noisy estimate of the gradient of the log posterior, which reduces the per iteration computational cost of the algorithm. Despite this, there are a number of results suggesting that stochastic gradient Langevin dynamics (SGLD), probably the most popular of these methods, still has computational cost proportional to the dataset size. We suggest an alternative log posterior gradient estimate for stochastic gradient MCMC, which uses control variates to reduce the variance. We analyse SGLD using this gradient estimate, and show that, under log-concavity assumptions on the target distribution, the computational cost required for a given level of accuracy is independent of the dataset size. Next we show that a different control variate technique, known as zero variance control variates can be applied to SGMCMC algorithms for free. This post-processing step improves the inference of the algorithm by reducing the variance of the MCMC output. Zero variance control variates rely on the gradient of the log posterior; we explore how the variance reduction is affected by replacing this with the noisy gradient estimate calculated by SGMCMC.

Equations

$${\theta }_t={\theta }_0-{\int }_0^t\nabla f({\theta }_s)ds+\sqrt{2}d{B}_t,$$

$$\nabla \hat{f}(\theta ):=\nabla f_0(\theta )+\frac{N}{n}\sum _{i\in S_k}\nabla f_i(\theta )$$

$${\theta }_{k+1}={\theta }_k-\frac{h_k}{2}\nabla \hat{f}({\theta }_k)+{\zeta }_k,$$

$$\nabla f(\theta )=\nabla f(\hat{\theta })+[\nabla f(\theta )-\nabla f(\hat{\theta })],$$

$$\left[\nabla \hat{f}(\theta )-\nabla \hat{f}(\hat{\theta })\right]=\nabla f_0(\theta )-\nabla f_0(\hat{\theta })+\frac{1}{n}\sum _{i\in S}\frac{1}{p_i}\left[\nabla f_i(\theta )-\nabla f_i(\hat{\theta })\right]$$

$$\nabla f_i(\theta )\approx \nabla f_i(\hat{\theta })\text{mbox}if\theta \approx \hat{\theta }.$$

$$\nabla \tilde{f}(\theta )=\nabla f(\hat{\theta })+[\nabla \hat{f}(\theta )-\nabla \hat{f}(\hat{\theta })],$$

$${\theta }_{k+1}={\theta }_k-h_k\nabla \hat{f}(\theta ),$$

$$f(\theta )-f({\theta }^{\prime })-\nabla f({\theta }^{\prime }{)}^{\top }(\theta -{\theta }^{\prime })$$

$$\Vert \nabla f(\theta )-\nabla f({\theta }^{\prime })\Vert$$

$$\left\lVert\nabla f_{i}(\theta)-\nabla f_{i}(\theta^{\prime})\right\rVert\leq L_{i}\left\lVert\theta-\theta^{\prime}\right\rVert,\mbox{ for $i=0,\ldots,N$}.$$

$$\mathbb{E}{\Vert {\xi }_k\Vert }^2\le \frac{{\left({\sum }_{i=1}^N{L}_i\right)}^2}{n}\mathbb{E}{\Vert {\theta }_k-\hat{\theta }\Vert }^2.$$

$$\mathbb{E}{\Vert {\xi }_k\Vert }^2\le \frac{M^2}{n}\mathbb{E}{\Vert {\theta }_k-\hat{\theta }\Vert }^2.$$

$$f_i(\theta )-f_i({\theta }^{\prime })-\nabla f_i({\theta }^{\prime }{)}^{\top }(\theta -{\theta }^{\prime })$$

$$\Vert \nabla f_i(\theta )-\nabla f_i({\theta }^{\prime })\Vert$$

$$f(\theta )-f({\theta }^{\prime })-\nabla f({\theta }^{\prime }{)}^{\top }(\theta -{\theta }^{\prime })$$

$$\Vert \nabla f(\theta )-\nabla f({\theta }^{\prime })\Vert$$

$$\mathbb{E}{\Vert \nabla \hat{f}(\theta )-\nabla f(\theta )\Vert }^2\le \frac{2N^2{\sigma }^2}{n},$$

$$W_q(\mu ,\nu )={\left[\underset{\gamma \in \Gamma (\mu ,\nu )}{\inf }{\int }_{{\mathbb{R}}^d\times {\mathbb{R}}^d}{\Vert \theta -{\theta }^{\prime }\Vert }^{q}d\gamma (\theta ,{\theta }^{\prime })\right]}^{\frac{1}{q}},$$

$$W_2({\nu }_k,\pi )\le (1-A{)}^K{W}_2({\nu }_0,\pi )+\frac{C}{A}+\frac{B^2}{C+\sqrt{A}B},$$

$$A$$

$$B$$

$$C$$

$$h$$

$$Kh$$

$$n$$

$$\beta =\max \left\{\frac{1}{2R^2+1},\frac{{ϵ}_0^2}{64R^2{\alpha }^{2}d}\right\},$$

$$C_1=\min \left\{2R^2+1,\frac{64R^2{\alpha }^{2}d}{{ϵ}_0^2}\right\}.$$

$$Kn\le \left[C_1\log \left[m\mathbb{E}{\Vert {\theta }_0-\bar{\theta }\Vert }^2+d\right]+C_1\log \frac{4}{{ϵ}_0^2}+1\right]\left[C_{2}m\mathbb{E}{\Vert \hat{\theta }-\bar{\theta }\Vert }^2+C_{2}d+1\right].$$

$$\mathbb{E}{\Vert \hat{\theta }-{\theta }^{\ast }\Vert }^2\le \frac{4D^2}{m^{2}K}.$$

$$\mathbb{E}[g(\theta )]\approx \frac{1}{M}\sum _{i=1}^{M}g({\theta }^{(i)}).$$

$$\tilde{g}(\theta )=g(\theta )+h(\theta ),$$

$$h(\theta )=\Delta Q(\theta )+\nabla Q(\theta )\cdot \mathbf{z},$$

$$\hat{h}(\theta )=\Delta Q(\theta )+\nabla Q(\theta )\cdot \hat{\mathbf{z}},$$

Full text

Keywords: Stochastic gradient MCMC; Langevin dynamics; scalable MCMC; control variates; computational cost; big data

1 Introduction

Markov chain Monte Carlo (MCMC), one of the most popular methods for Bayesian inference, scales poorly with dataset size. This is because standard methods require the whole dataset to be evaluated at each iteration. Stochastic gradient MCMC (SGMCMC) are a class of MCMC algorithms that aim to account for this issue. The algorithms have recently gained popularity in the machine learning literature, though they were originally proposed in Lamberton and Pagès (2002). These methods use efficient MCMC proposals based on discretised dynamics that use gradients of the log posterior. They reduce the computational cost by replacing the gradients with an unbiased estimate which uses only a subset of the data, referred to as a minibatch. They also bypass the acceptance step by making small discretisation steps (Welling and Teh, 2011; Chen et al., 2014; Ding et al., 2014; Dubey et al., 2016).

These new algorithms have been successfully applied to a range of state of the art machine learning problems (e.g. Patterson and Teh, 2013; Wenzhe Li, 2016). There is a variety of software available implementing these methods (Tran et al., 2016; Baker et al., 2016). In particular Baker et al. (2016) implements the control variate methodology we discuss in this article.

This paper investigates stochastic gradient Langevin dynamics (SGLD), a popular SGMCMC algorithm that discretises the Langevin diffusion. There are a number of results suggesting that while SGLD has a lower per iteration computational cost compared with MCMC, its overall computational cost is proportional to the dataset size (Welling and Teh, 2011; Nagapetyan et al., 2017). This motivates improving the computational cost of SGLD, which can be done by using control variates (Ripley, 2009). Control variates can be applied to reduce the Monte Carlo variance of the gradient estimate in stochastic gradient MCMC algorithms. We refer to SGLD using this new control variate gradient estimate as SGLD-CV.

We analyse the algorithm using the Wasserstein distance between the distribution defined by SGLD-CV and the true posterior distribution, by adapting recent results by Dalalyan and Karagulyan (2017). Our results are based on assuming the target distribution is strongly log-concave. We get bounds on the Wasserstein distance between the target distribution and the distribution we sample from at a given step of SGLD-CV. These bounds are in terms of the tuning constants chosen when implementing SGLD-CV. By making assumptions on how the posterior distribution changes with the number of data points, we are able to show that the computational cost required for a given level of accuracy does not grow with the dataset size. Though this is providing we can obtain an estimate, $\hat{\theta}$, of the posterior mode that is sufficiently close. Our results also show the impact on the computational cost of using a poor $\hat{\theta}$ value.

The algorithm requires some additional preprocessing steps before the computational cost benefits come into effect. These preprocessing steps include finding $\hat{\theta}$ and calculating the full log posterior gradient at $\hat{\theta}$. Both of these steps have a cost that is linear in the amount of data. However, the cost of finding $\hat{\theta}$ essentially replaces the burn-in of the chain, and we find empirically that this is often more efficient. The cost of these steps is analysed in more detail in the article.

The use of control variates has also been shown to be important for other Monte Carlo algorithms for simulating from a posterior with a cost that is sub-linear in the number of data points (Bardenet et al., 2016; Bierkens et al., 2016; Pollock et al., 2016; Nagapetyan et al., 2017). For previous work that suggests using control variates within SGLD see Dubey et al. (2016) and Chen et al. (2017). These latter papers, whilst showing benefits of using control variates, do not show that the resulting algorithm can have sub-linear cost in the number of data points. A recent paper, Nagapetyan et al. (2017), does investigate how SGLD-CV peforms in the limit as you have more data under similar log-concavity assumptions on the posterior distribution. They have results that are qualitatively similar to ours, including the sub-linear computational cost of SGLD-CV. Though they measure accuracy of the algorithm through the mean square error of Monte Carlo averages rather than through the Wasserstein distance.

Not only can control variates be used to speed up stochastic gradient MCMC by enabling smaller minibatches to be used; we show that they can be used to improve the inferences made from the MCMC output. In particular, we can use post-processing control variates (Mira et al., 2013; Papamarkou et al., 2014; Friel et al., 2016) to produce MCMC samples with a reduced variance. The post-processing methods rely on the MCMC output as well as gradient information. Since stochastic gradient MCMC methods already compute estimates of the gradient, we explore replacing the true gradient in the post-processing step with these free estimates. We also show theoretically how this affects the variance reduction factor; and empirically demonstrate the variance reduction that can be achieved from using these post-processing methods.

2 Stochastic gradient MCMC

Throughout this paper we aim to make inference on a vector of parameters $\theta\in\mathbb{R}^{d}$, with data $\mathbf{x}=\{x_{i}\}_{i=1}^{N}$. We denote the probability density of $x_{i}$ as $p(x_{i}|\theta)$ and assign a prior density $p(\theta)$. The resulting posterior density is then $p(\theta|\mathbf{x})\propto p(\theta)\prod_{i=1}^{N}p(x_{i}|\theta)$, which defines the posterior distribution $\pi$. For brevity we write $f_{i}(\theta)=-\log p(x_{i}|\theta)$ for $i=1,\dots N$, $f_{0}(\theta)=-\log p(\theta)$ and $f(\theta)=-\log p(\theta|\mathbf{x})$.

Many MCMC algorithms are based upon discrete-time approximations to continuous-time dynamics, such as the Langevin diffusion, that are known to have the posterior as their invariant distribution. The approximate discrete-time dynamics are then used as a proposal distribution within a Metropolis-Hastings algorithm. The accept-reject step within such an algorithm corrects for any errors in the discrete-time dynamics. Examples of such an approach include the Metropolis-adjusted Langevin algorithm (MALA; see e.g. Roberts and Rosenthal (1998)) and Hamiltonian Monte Carlo (HMC; see Neal (2010)).

2.1 Stochastic gradient Langevin dynamics

SGLD, first introduced by Lamberton and Pagès (2002), and popularised more recently by Welling and Teh (2011), is a minibatch version of the Metropolis-adjusted Langevin algorithm. At each iteration it creates an approximation of the true gradient of the log-posterior by using a small sample of data.

The SGLD algorithm is based upon the discretisation of a stochastic differential equation known as the Langevin diffusion. A Langevin diffusion for a parameter vector $\theta$ with posterior $p(\theta|\mathbf{x})\propto\exp(-f(\theta))$ is given by

[TABLE]

where $B_{t}$ is a $d$-dimensional Wiener process. The stationary distribution of this diffusion is $\pi$. This means that it will target the posterior exactly, but in practice we need to discretize the dynamics to simulate from it, which introduces error. A bottleneck for this simulation is that calculating $\nabla f(\theta)$ is an $O(N)$ operation. So to get around this, Welling and Teh (2011) replace the log posterior gradient with the following unbiased estimate

[TABLE]

for some subsample $S_{k}$ of $\{1,\dots,N\}$, with $|S_{k}|=n$. A single update of SGLD is then

[TABLE]

where $\zeta_{k}\sim N(0,h_{k})$.

MALA uses a Metropolis-Hastings accept-reject step to correct for the discretisation of the Langevin process. Welling and Teh (2011) bypass this acceptance step, as it requires calculating $p(\theta|\mathbf{x})$ using the full dataset, and instead use an adaptive rather than fixed stepsize, where $h_{k}\rightarrow 0$ as $k\rightarrow\infty$. The motivation is that the noise in the gradient estimate disappears faster than the process noise, so eventually, the algorithm will sample the posterior approximately. In practice, we found the algorithm does not mix well when the stepsize is decreased to zero, so in general a fixed small stepsize $h$ is used in practice, as suggested by Vollmer et al. (2016).

3 Control variates for SGLD efficiency

The SGLD algorithm has a reduced per iteration computational cost compared to traditional MCMC algorithms. However, there have been a number of results suggesting that the overall computational cost of SGLD is still $O(N)$ (Welling and Teh, 2011; Nagapetyan et al., 2017). The main reason for this result is that in order to control the variance in the gradient estimate, $\nabla\hat{f}(\theta)$, we need $n$ to increase linearly with $N$. So intuition would suggest trying to reduce this variance would reduce the computational cost of the algorithm.

A natural choice is to reduce this variance through control variates (Ripley, 2009). Control variates applied to SGLD have also been investigated by Dubey et al. (2016) and Chen et al. (2017), who show that the convergence bound of SGLD is reduced when they are used. Theoretical results, similar to ours below, on how the use of control variates can improve how the computational cost of SGLD scales with $N$ are given in Nagapetyan et al. (2017).

In Section 3.1, we show how control variates can be used to reduce the variance in the gradient estimate in SGLD, leading to the algorithm SGLD-CV. Then in Section 3.3 we analyse the Wasserstein distance between the distribution defined by SGLD-CV and the true posterior. There are a number of quantities that affect the performance of SGLD-CV, including the stepsize $h$, the number of iterations $K$ and the minibatch size $n$. We provide sufficient conditions on $h$, $K$ and $n$ in order to bound the Wasserstein distance. We show under certain assumptions, the computational cost, measured as $Kn$, required to bound the Wasserstein distance is independent of $N$.

3.1 Control variates for SGMCMC

Let $\hat{\theta}$ be a fixed value of the parameter, chosen to be close to the mode of the posterior $p(\theta|\mathbf{x})$. The log posterior gradient can then be re-written as

[TABLE]

where the first term on the right-hand side is a constant and the bracketed term on the right-hand side can be unbiasedly estimated by

[TABLE]

where $p_{1},\ldots,p_{N}$ are user-chosen, strictly positive probabilities, $S$ is a random sample from $\{1,\ldots,N\}$ such that $|S|=n$ and the expected number of times $i$ is sampled is $np_{i}$. The standard implementation of control variates would set $p_{i}=1/N$ for all $i$. Yet we show below that there can be advantages in having these probabilities vary with $i$; for example to give higher probabilities to sampling data points for which $\nabla f_{i}(\theta)-\nabla f_{i}(\hat{\theta})$ has higher variability.

If the gradient of the likelihood for a single observation is smooth in $\theta$ then we will have

[TABLE]

Hence for $\theta\approx\hat{\theta}$ we would expect the unbiased estimator

[TABLE]

to have a lower variance than the simpler unbiased estimator (2). This is because when $\theta$ is close to $\hat{\theta}$ we would expect the terms $\nabla\hat{f}(\theta)$ and $\nabla\hat{f}(\hat{\theta})$ to be correlated. This reduction in variance is shown formally in Lemma 1, stated in Section 3.2.

The gradient estimate (4) can be substituted into any stochastic gradient MCMC algorithm in place of $\nabla\hat{f}(\theta)$. We refer to SGLD using this alternative gradient estimate as SGLD-CV. The full procedure is outlined in Algorithm 1.

Implementing this in practice means finding a suitable $\hat{\theta}$, which we refer to as the centering value. We show below that for the computational cost of SGLD-CV to be $O(1)$, we require both $\hat{\theta}$ and the starting point of SGLD-CV, $\theta_{0}$, to be a $O(N^{-\frac{1}{2}})$ distance from the posterior mean.

In practice, we find $\hat{\theta}$ using stochastic optimisation (Robbins and Monro, 1951), and then calculate the full log posterior gradient at this point $\nabla f(\hat{\theta})$. We then start the algorithm from $\hat{\theta}$. In our implementations we use a simple stochastic optimisation method, known as stochastic gradient descent (SGD, see e.g. Bottou, 2010). The method works similarly to the standard optimisation method gradient descent, but at each iteration replaces the true gradient of the function with an unbiased estimate. A single update of the algorithm is as follows

[TABLE]

where $\nabla\hat{f}(\theta)$ is as defined in (2) and $h_{k}>0$ is a small tuning constant referred to as the stepsize. Provided the stepsizes $h_{k}$ satisfy the following conditions $\sum_{k}h_{k}^{2}<\infty$ and $\sum_{k}h_{k}=\infty$ then this algorithm will converge to a local maximum.

We show in Section 3.4, under our assumptions of log-concavity of the posterior, that finding $\hat{\theta}$ using SGD has a computational cost that is linear in $N$, and we can achieve the required accuracy with just a single pass through the data. As we then start SGLD-CV with this value for $\theta$, we can view finding the centering value as a replacement for the burn-in phase of the algorithm, and we find, in practice, that the time to find a good $\hat{\theta}$ is often quicker than the time it takes for SGLD to burn-in. One downside of this procedure is that the SGD algorithm, as well as the SGLD-CV algorithm itself needs to be tuned, which adds to the tuning burden.

In comparison to SGLD-CV, the SAGA algorithm by Dubey et al. (2016) also uses control variates to reduce the variance in the gradient estimate of SGLD. They show that this reduces the MSE of SGLD. The main difference is that their algorithm uses a previous state in the chain as the control variate, rather than an estimate of the mode. This means that SAGA does not require the additional optimisation step, so tuning should be easier. However we show in the experiments of Section 5, that the algorithm gets more easily stuck in local stationary points, especially during burn-in. For more complex examples, the algorithm was prohibitively slow to burn-in because of this tendency to get stuck. Dubey et al. (2016) also do not show that SAGA has favourable computational cost results.

3.2 Variance reduction

The improvements of using the control variate gradient estimate (4) over the standard (2) become apparent when we calculate the variances of each. For our analysis, we make the assumption that the posterior is strongly log-concave, formally defined in Assumption 1. This has become a common assumption when analysing gradient based samplers that do not have an acceptance step (Durmus and Moulines, 2016; Dalalyan and Karagulyan, 2017). In all the following analysis we use $\left\lVert\cdot\right\rVert$ to denote the Euclidean norm.

Assumption 1**.**

Strongly log-concave posterior: there exists positive constants $m$ and $M$, such that the following conditions hold for the negative log posterior

[TABLE]

for all $\theta,\theta^{\prime}\in\mathbb{R}^{d}$.

We further need a Lipschitz condition for each of the likelihood terms in order to bound the variance of our control-variate estimator of the gradient.

Assumption 2**.**

Lipschitz: there exists constants $L_{0},\ldots,L_{N}$ such that

[TABLE]

Using Assumption 2 we are able to derive a bound on the variance of the gradient estimate of SGLD-CV. This bound is formally stated in Lemma 1.

Lemma 1**.**

Under Assumption 2. Let $\theta_{k}$ be the state of SGLD-CV at the $k^{th}$ iteration, with stepsize $h$ and centering value $\hat{\theta}$. Assume we estimate the gradient using the control variate estimator with $p_{i}=L_{i}/\sum_{j=1}^{N}L_{j}$ for $i=1,\ldots,N$. Define $\xi_{k}:=\nabla\tilde{f}(\theta_{k})-\nabla f(\theta_{k})$, so that $\xi_{k}$ measures the noise in the gradient estimate $\nabla\tilde{f}$ and has mean 0. Then for all $\theta_{k},\hat{\theta}\in\mathbb{R}^{d}$, for all $k=1,\dots,K$ we have

[TABLE]

All proofs are relegated to the Appendix. If Assumption 1 also holds, then we can choose an $M=\sum_{i=0}^{N}L_{i}$, and our bound (8) implies

[TABLE]

We will use this form of the bound for the rest of the analysis. While it looks like picking $p_{i}$ will require estimates of the Lipschitz constants $L_{i}$; in practice, under Assumption 3 stated below, we can just use the standard $p_{i}=1/N$ for all $i$. We use $p_{i}=1/N$ in all our implementations in the experiments of Section 5.

In order to consider how SGLD-CV scales with $N$ we need to make assumptions on the properties of the posterior and how these change with $N$. To make discussions concrete we will focus on the following, strong, assumption that each likelihood-term in the posterior is strongly log-concave. As we discuss later, our results apply under weaker conditions.

Assumption 3**.**

Assume there exists positive constants $L$ and $l$ such that $f_{i}$ satisfies the following conditions

[TABLE]

for all $i\in 0,\dots,N$ and $\theta,\theta^{\prime}\in\mathbb{R}^{d}$.

Under this assumption the log-concavity constants, $m$ and $M$, of the posterior both increase linearly with $N$, as shown by the following Lemma.

Lemma 2**.**

Suppose Assumption 3 holds. Then the log-posterior, $f$, satisfies the following

[TABLE]

Thus the posterior is strongly log-concave with parameters $M=(N+1)L$ and $m=(N+1)l$.

To see the potential benefit of using control variates to estimate the gradient in situations where $N$ is large, we can now compare the variance bound from Lemma 1, with a bound on the variance of the simple estimator, $\nabla\hat{f}(\theta)$. If we assume that $\left\lVert\nabla f_{i}(\theta)\right\rVert$ is bounded by some constant $\sigma$, for all $i=0,\dots,N$ and for all $\theta\in\mathbb{R}^{d}$, then Dubey et al. (2016) show that for SGLD

[TABLE]

for all $\theta\in\mathbb{R}^{d}$.

We can see that the bound on the gradient estimate variance in (8) depends on the distance between $\theta_{k}$ and $\hat{\theta}$. Appealing to the Bernstein-von Mises theorem (Le Cam, 2012), under standard asymptotics we would expect the distance $\mathbb{E}\left\lVert\theta_{k}-\hat{\theta}\right\rVert^{2}$ to be $O(1/N)$, if $\hat{\theta}$ is within $O(N^{-1/2})$ of the posterior mean, once the MCMC algorithm has burnt in. As $M$ is $O(N)$, this suggests that using control variates could give an $O(N)$ reduction in variance, and this plays a key part in the computational cost improvements we show in the next section.

3.3 Computational cost of SGLD-CV

In this section, we investigate how applying control variates to the gradient estimate of SGLD reduces the computational cost of the algorithm.

In order to show this, we investigate the Wasserstein-Monge-Kantorovich (Wasserstein) distance $W_{2}$ between the distribution defined by the SGLD-CV algorithm at each iteration and the true posterior as $N$ is changed. For two measures $\mu$ and $\nu$ defined on the probability space $(\mathbb{R}^{d},B(\mathbb{R}^{d}))$, and for a real number $q>0$, the distance $W_{q}$ is defined by

[TABLE]

where the infimum is with respect to all joint distributions $\Gamma$ having $\mu$ and $\nu$ as marginals. The Wasserstein distance is a natural distance measure to work with for Monte Carlo algorithms, as discussed in Dalalyan and Karagulyan (2017); Durmus and Moulines (2016).

One issue when working with the Wasserstein distance is that it is not invariant to transformations. For example scaling all entries of $\theta$ by a constant will scale the Wasserstein distance by the same constant. A linear transformation of the parameters will result in a posterior that is still strongly log-concave, but with different constants $m$ and $M$. To account for this we suggest measuring error by the quantity $\sqrt{m}W_{2}$, which is invariant to scaling $\theta$ by a constant. Theorem 1 of Durmus and Moulines (2016) bounds the standard deviation of any component of $\theta$ by a constant times $1/\sqrt{m}$, so we can view the quantity $\sqrt{m}W_{2}$ as measuring the error on a scale that is relative to the variability of $\theta$ under the posterior distribution.

There are a number of quantities that will affect the performance of SGLD and SGLD-CV. These include the step size $h$, the minibatch size $n$ and the total number of iterations $K$. In the analysis to follow we find conditions on $h$, $n$ and $K$ that ensure the Wasserstein distance between the distribution defined by SGLD-CV and the true posterior distribution $\pi$ are less than some $\epsilon>0$. We use these conditions to calculate the computational cost, measured as $Kn$, required for this algorithm to reach the satisfactory error $\epsilon$.

The first step is to find an upper bound on the Wasserstein distance between SGLD-CV and the posterior distribution in terms of $h$, $n$, $K$ and the constants $m$ and $M$ declared in Assumption 1.

Proposition 1**.**

Under Assumptions 1 and 2, let $\theta_{k}$ be the state of SGLD-CV at the $k^{th}$ iteration of the algorithm with stepsize $h$, initial value $\theta_{0}$, centering value $\hat{\theta}$. Let the distribution of $\theta_{k}$ be $\nu_{k}$. Denote the expectation of $\theta$ under the posterior distribution $\pi$ by $\bar{\theta}$. If $h<\frac{2m}{2M^{2}+m^{2}}$, then for all integers $k\geq 0$,

[TABLE]

where

[TABLE]

$\alpha=7\sqrt{2}/6$* and $d$ is the dimension of $\theta_{k}$.*

The proof of this proposition is closely related to the proof of Proposition 2 of Dalalyan and Karagulyan (2017). The extra complication comes from our bound on the variance of the estimator of the gradient; which depends on the current state of the SGLD-CV algorithm, rather than being bounded by a global constant.

We can now use Proposition 1 to find conditions on $K$, $h$ and $n$ in terms of the constants $M$ and $m$ such that the Wasserstein distance is bounded by some positive constant $\epsilon_{0}/\sqrt{m}$ at its final iteration $K$.

Theorem 1**.**

Under Assumptions 1 and 2, let $\theta_{K}$ be the state of SGLD-CV at the $K^{th}$ iteration of the algorithm with stepsize $h$, initial value $\theta_{0}$, centering value $\hat{\theta}$. Let the distribution of $\theta_{k}$ be $\nu_{K}$. Denote the expectation of $\theta$ under the posterior distribution $\pi$ by $\bar{\theta}$. Define $R:=M/m$. Then for any $\epsilon_{0}>0$, if the following conditions hold:

[TABLE]

where

[TABLE]

$\alpha=7\sqrt{2}/6$, and $d$ is the dimension of $\theta_{k}$, then $W_{2}(\nu_{K},\pi)\leq\epsilon_{0}/\sqrt{m}$.

As a corollary of this result, we have the following, which gives a bound on the computational cost of SGLD, as measured by $Kn$, to achieve a required bound on the Wasserstein distance.

Corollary 1**.**

Assume that Assumptions 1 and 3 and the conditions of Theorem 1 hold. Fix $\epsilon_{0}$ and define

[TABLE]

and $C_{2}:=64R^{2}\beta/\epsilon_{0}^{2}$. We can implement an SGLD-CV algorithm with $W_{2}(\nu_{K},\pi)<\epsilon_{0}/\sqrt{m}$ such that

[TABLE]

The constants, $C_{1}$ and $C_{2}$, in the bound on $Kn$, depend on $\epsilon_{0}$ and $R=M/m$. It is simple to show that both constants are increasing in $R$. Under Assumption 3 we have that $R$ is a constant as $N$ increases. Corollary 1 suggests that provided $\left\lVert\theta_{0}-\bar{\theta}\right\rVert<c/\sqrt{m}$ and $\left\lVert\hat{\theta}-\bar{\theta}\right\rVert<c/\sqrt{m}$, for some constant $c$; then the computational cost of SGLD-CV will be bounded by a constant. Since we suggest starting SGLD-CV at $\hat{\theta}$, then technically we just need this to hold for $\left\lVert\hat{\theta}-\bar{\theta}\right\rVert$. Under Assumption 3 we have that $m$ increases linearly with $N$, so this corresponds to needing $\left\lVert\hat{\theta}-\bar{\theta}\right\rVert<c_{1}/\sqrt{N}$ as $N$ increases. Additionally, by Theorem 1 of Durmus and Moulines (2016) we have that the variance of the posterior scales like $1/m=1/N$ as $N$ increases, so we can interpret the $1/\sqrt{N}$ factor as being a measure of the spread of the posterior as $N$ increases. The form of the corollary makes it clear that a similar argument would apply under weaker assumptions than Assumption 3. We only need that the ratio of the log-concavity constants, $M/m$, of the posterior remains bounded as $N$ increases.

This corollary also gives insight into the computational penalty you pay for a poor choice of $\theta_{0}$ or $\hat{\theta}$. As $\left\lVert\theta_{0}-\bar{\theta}\right\rVert$ increases, the bound on the computational cost will increase logarithmically with this distance. By comparison the bound increases linearly with $\left\lVert\hat{\theta}-\bar{\theta}\right\rVert$.

3.4 Setup Costs

There are a number of known results on the convergence of SGD under the strongly log-concave conditions of Assumption 1. These will allow us to quantify the setup cost of finding the point $\hat{\theta}$ in this setting. More complex cases are explored empirically in the experiments in Section 5. Lemma 3 due to Nemirovski et al. (2009) quantifies the convergence of the final point of SGD.

Lemma 3**.**

Under Assumption 1, let $\hat{\theta}$ denote the final state of SGD with stepsizes $h_{k}=1/(mk)$ after $K$ iterations. Suppose $\mathbb{E}\left\lVert\nabla\hat{f}(\theta)\right\rVert^{2}\leq D^{2}$ and denote the true mode of $f$ by $\theta^{*}$. Then it holds that

[TABLE]

If we again assume that, as in Dubey et al. (2016), $\left\lVert\nabla f_{i}(\theta)\right\rVert$ is bounded by some constant $\sigma$, for all $i=0,\dots,N$ and $\theta\in\mathbb{R}^{d}$ then $D^{2}$ will be $O(N^{2}/n)$. This means that under Assumption 3, we will need to process the full dataset once before the SGD algorithm has converged to a mode $\hat{\theta}$ within $O(N^{-\frac{1}{2}})$ of the posterior mean. It follows that, for these cases there are two one off $O(N)$ setup costs, one to find an acceptable mode $\hat{\theta}$ and one to find the full log posterior gradient at this point $\nabla f(\hat{\theta})$.

4 Post-processing control variates

Control variates can also been used to improve the inferences made from MCMC by reducing the variance of the output directly. The general aim of MCMC is to estimate expectations of functions, $g(\theta)$, with respect to the posterior $\pi$. Given an MCMC sample $\theta^{(1)},\dots,\theta^{(M)}$, from the posterior $\pi$, we can estimate $\mathbb{E}[g(\theta)]$ unbiasedly as

[TABLE]

Suppose there exists a function $h(\theta)$, which has expectation 0 under the posterior. We can then introduce an alternative function,

[TABLE]

where $\mathbb{E}[\tilde{g}(\theta)]=\mathbb{E}[g(\theta)]$. If $h(\cdot)$ is chosen so that it is negatively correlated with $g(\theta)$, then the variance of $\tilde{g}(\theta)$ will be reduced considerably.

Mira et al. (2013) introduce a way of choosing $h(\theta)$ almost automatically by using the gradient of the log-posterior. Choosing $h(\cdot)$ in this manner is referred to as a zero variance (ZV) control variate. Friel et al. (2016) showed that, under mild conditions, we can replace the log-posterior gradient with an unbiased estimate and still have a valid control variate. SGMCMC methods produce unbiased estimates of the log-posterior gradient, and so it follows that these gradient estimates can be applied as ZV control variates. For the rest of this section, we focus our attention on SGLD, but these ideas are easily extendable to other stochastic gradient MCMC algorithms. We refer to SGLD with these post-processing control variates as SGLD-ZV.

Given the setup outlined above, Mira et al. (2013) propose the following form for $h(\theta)$,

[TABLE]

here $Q(\theta)$ is a polynomial of $\theta$ to be chosen and $\mathbf{z}=f(\theta)/2$. $\Delta$ refers to the Laplace operator $\frac{\partial^{2}}{\partial\theta_{1}^{2}}+\dots+\frac{\partial^{2}}{\partial\theta_{d}^{2}}$. In order to get the best variance reduction, we simply have to optimize the coefficients of the polynomial $Q(.)$. In practice, first or second degree polynomials $Q(\theta)$ often provide good variance reduction (Mira et al., 2013). For the rest of this section we focus on first degree polynomials, so $Q(\theta)=\mathbf{a}^{T}\theta$, but the ideas are easily extendable to higher orders (Papamarkou et al., 2014).

The SGLD algorithm only calculates an unbiased estimate of $\nabla f(\theta)$, so we propose replacing $h(\theta)$ with the unbiased estimate

[TABLE]

where $\mathbf{\hat{z}}=\nabla\hat{f}(\theta)/2$. By identical reasoning to Friel et al. (2016), $\hat{h}(\theta)$ is a valid control variate. Note that $\hat{\mathbf{z}}$ can use any unbiased estimate, and as we will show later, the better the gradient estimate, the better this control variate performs.

We set $Q(\theta)$ to be a linear polynomial $\mathbf{a}^{T}\theta$, so our SGLD-ZV estimate will take the following form

[TABLE]

Similar to standard control variates (Ripley, 2009), we need to find optimal coefficients $\mathbf{\hat{a}}$ in order to minimize the variance of $\tilde{g}(\cdot)$, defined in (14). In our case, the optimal coefficients take the following form (Friel et al., 2016)

[TABLE]

This means that SGLD already calculates all the necessary terms for these control variates to be applied for free. So the post-processing step can simply be applied once when the SGLD algorithm has finished, provided the full output plus gradient estimates are stored. With this in place, we can write down the full algorithm in the linear case, which is given in Algorithm 2. For higher order polynomials, the calculations are much the same, but more coefficients need to be estimated (Papamarkou et al., 2014).

The efficiency of ZV control variates in reducing the variance of our MCMC sample is directly affected by using an estimate of the gradient rather than the truth. For the remainder of this section, we investigate how the choice of the gradient estimate, and the minibatch size $n$, affects the variance reduction.

Assumption 4**.**

$\mathbb{V}\mathrm{ar}[\phi(\theta)]<\infty$* and $\mathbb{V}\mathrm{ar}[\hat{\psi}(\theta)]<\infty$. $\mathbb{E}_{\theta|\mathbf{x}}\left\lVert\nabla f_{i}(\theta)\right\rVert^{2}$ is bounded by some constant $\sigma$ for all $i=0,\dots N$, $\theta\in\mathbb{R}^{d}$.*

Theorem 2**.**

Under Assumption 4, define the optimal variance reduction for ZV control variates using the full gradient estimate to be $R$, and the optimal variance reduction using SGLD gradient estimates to be $\hat{R}$. Then we have that

[TABLE]

where $\xi_{S}(\theta)$ is the noise in the log-posterior gradient estimate.

The proof of this result is given in the Appendix. An important consequence of Theorem 2 is that if we use the standard SGLD gradient estimate, then the denominator of (15) is $O(n/N)$, so our variance reduction diminishes as $N$ gets large. However, if we use the SGLD-CV estimate instead (the same probably holds for other control variate algorithms such as SAGA), then under standard asymptotics, the denominator of (15) is $O(n)$, so the variance reduction does not diminish with increasing dataset size. It follows that for best results, we recommend using the ZV post-processing step after running the SGLD-CV algorithm, especially for large $N$. The ZV post-processing step can be immediately applied in exactly the same way to other stochastic gradient MCMC algorithms, such as SGHMC and SGNHT (Chen et al., 2014; Ding et al., 2014).

It is worth noting that there are some storage constraints for SGLD-ZV. This algorithm requires storing the full MCMC chain, as well as the gradient estimates at each iteration. So the storage cost is twice the storage cost of a standard SGMCMC run. However, in some high dimensional cases, the required SGMCMC test statistic is estimated on the fly using the most recent output of the chain and thus reducing the storage costs. We suggest that if the dimensionality is not too high, then the additional storage cost of recording the gradients to apply the ZV post-processing step can offer significant variance reduction for free. However, for very high dimensional parameters, the cost associated with storing the gradients may preclude the use of the ZV step.

5 Experiments

5.1 Logistic regression

We examine our approaches on a Bayesian logistic regression problem. The probability of the $i^{th}$ output $y_{i}\in\{-1,+1\}$ is given by

[TABLE]

We use a Laplace prior for $\beta$ with scale 1.

We used the cover type dataset (Blackard and Dean, 1999), which has 581012 observations, which we split into a training and test set. First we run SGLD, SGLD-CV and SAGA on the dataset, all with minibatch size 500. The method SAGA was discussed at the end of Section 3.1. To empirically support the scalability results of Theorem 1, we fit the model 3 times. In each fit, the dataset size is varied, from about 1% of the full dataset to the full dataset size $N$. The performance is measured by calculating the log predictive density on a held-out test set every 10 iterations. Some of our examples are high dimensional, so our performance measure aims to reduce dimensionality while still capturing important quantities such as the variance of the chain. We include the burn-in of SGLD and SAGA, to contrast with the optimisation step required for SGLD-CV which is included in the total computational time.

The results are plotted against time in Figure 1. The results illustrate the efficiency gains of SGLD-CV over SGLD as the dataset size increases, as expected from Theorem 1. SAGA outperforms SGLD-CV in this example because SGLD converges quickly in this simple setting. In the more complicated examples to follow, we show that SAGA can be slow to converge.

We also compare the log predictive density over a test set for SGLD-CV with and without ZV post-processing, averaged over 5 runs at different seeds. We apply the method to SGLD-CV rather than SGLD due to the favourable scaling results as discussed after Theorem 2. Results are given in Figure 2. The plot shows box-plots of the log predictive density of the SGLD sample before and after post-processing using ZV control variates. The plots show excellent variance reduction of the chain.

5.2 Probabilistic matrix factorization

A common recommendation system task is to predict a user’s rating of a set of items, given previous ratings and the ratings of other users. The end goal is to recommend new items that the user will rate highly. Probabilistic matrix factorization (PMF) is a popular method to train these models (Mnih and Salakhutdinov, 2008). As the matrix of ratings is sparse, over-fitting is a common issue in these systems, and Bayesian approaches are a way to account for this (Ahn et al., 2015).

In this experiment, we apply SGLD, SGLD-CV and SAGA to a Bayesian PMF problem, using a model similar to Ahn et al. (2015) and Chen et al. (2014). We use the Movielens dataset ml-100k222https://grouplens.org/datasets/movielens/100k/, which contains 100,000 ratings from almost 1,000 users and 1,700 movies. We use batch sizes of 5,000, with a larger minibatch size chosen due to the high-dimensional parameter space. As before, we compare performance by calculating the predictive distribution on a held out dataset every 10 iterations.

We investigate the scaling results of SGLD-CV and SAGA versus SGLD by varying the dataset size. We do this by limiting the number of users in the dataset, ranging from 100 users to the full 943. The results are given in Figure 3. Once again the scaling improvements of SGLD-CV as the dataset size increases are clear.

In this example SAGA converges slowly in comparison even to SGLD. In fact the algorithm converges slowly in all our more complex experiments. The problem is particularly bad for large $N$. This is likely a result of the starting point for SAGA being far from the posterior mode. Empirically, we found that the gradient direction and magnitude can update very slowly in these cases. This is not an issue for simpler examples such as logistic regression, but for more complex examples we believe it could be a sign that the algorithm is getting stuck in, or moving slowly through, local modes where the gradient is comparatively flatter. The problem appears to be made worse for large $N$ when it takes longer to update $g_{\alpha}$. This is an example where the optimisation step of SGLD-CV is an advantage, as the algorithm is immediately started close to the posterior mode and so the efficiency gains are quickly noted. This issue with SAGA could be related to the starting point condition for SGLD-CV as detailed in Corollary 1. Due to the form of the Wasserstein bound, it is likely that SAGA would have a similar starting point condition.

Once again we compare the log predictive density over a test set for SGLD-CV with and without ZV post-processing when applied to the Bayesian PMF problem, averaged over 5 runs at different seeds. Results are given in Figure 2. The plot shows box-plots of the log predictive density of the SGLD sample before and after post-processing using ZV control variates. The plots show excellent variance reduction of the chain.

5.3 Latent Dirichlet allocation

Latent Dirichlet allocation (LDA) is an example of a topic model used to describe collections of documents by sets of discovered topics (Blei et al., 2003). The input consists of a matrix of word frequencies in each document, which is very sparse, motivating the use of a Bayesian approach to avoid over-fitting.

Due to storage constraints, it was not feasible to apply SGLD-ZV to this problem, so we focus on SGLD-CV. We scraped approximately 80,000 documents from Wikipedia, and used the 1,000 most common words to form our document-word matrix input. We used a similar formulation to Patterson and Teh (2013), though we did not use a Riemannian sampler.

Once again in our comparison of SGLD, SGLD-CV and SAGA, we vary the dataset size, this time by changing the number of documents used in fitting the model, from 10,000 to the full 81,538. We use batch sizes of 50 documents. We measure the performance of LDA using the perplexity on held out words from each document, a standard performance measure for this model. The results are given in Figure 5. Here the scalability improvements of using SGLD-CV over SGLD are clear as the dataset size increases. This time the batch size is small compared to the dataset size, which probably makes the scalability improvements more obvious. The sudden drop in perplexity for the SGLD-CV plot occurs at the switch from the stochastic optimization step to SGLD-CV. This is likely a result of the algorithm making efficient use of the Gibbs step to simulate the latent topics.

An interesting aspect of this problem is that it appears to have a pronounced local mode where each of the methods become trapped (this can be seen by the blip in the plot at a perplexity of around 850). SGLD-CV is the first to escape followed by SGLD, but SAGA takes a long time to escape. This is probably due to a similar aspect as the one discussed in the previous experiment (Section 5.2). Similar to the previous experiment, we find that while SAGA seems trapped, its gradient estimate changes very little, which could be a sign that the algorithm is moving very slowly through an area with a relatively flat gradient, such as a local mode. A simple solution would be to start SAGA closer to the mode using a stochastic optimisation scheme.

6 Discussion

We have used control variates for stochastic gradient MCMC to reduce the variance in the gradient estimate. We have shown that in the strongly log-concave setting, and under standard asymptotics, this proposed SGLD-CV algorithm reduces the computational cost of stochastic gradient Langevin dynamics to $O(1)$. Our theoretical results give results on the computational cost under non-standard asymptotics also, and show there should be some benefit provided distance between the centering value $\hat{\theta}$ and the posterior mean inversely depends on $N$. The algorithm relies on a setup cost that estimates the posterior mode which replaces the burn-in of SGLD. We have explored the cost of this step both theoretically and empirically. We have empirically supported these scalability results on a variety of interesting and challenging problems from the statistics and machine learning literature using real world datasets. The simulation study also revealed that SGLD-CV was less susceptible to getting stuck in local stationary points than an alternative method that performs variance reduction using control variates, SAGA (Dubey et al., 2016). An interesting future extension would be to reduce the startup cost of SGLD-CV, along with introducing automatic step-size tuning.

We showed that stochastic gradient MCMC methods calculate all the information needed to apply zero variance post-processing control variates. This improves the inference of the output by reducing its variance. We explored how the variance reduction is affected by the minibatch size and the gradient estimate, and show using SGLD-CV or SAGA rather than SGLD can achieve a better variance reduction. We demonstrated this variance reduction empirically. A limitation of these post-processing control variates is they require the whole chain, which can lead to high storage costs if the dimensionality of the sample space is high. Future work could explore ways to reduce the storage costs of stochastic gradient MCMC.

7 Acknowledgements

The first author gratefully acknowledges the support of the EPSRC funded EP/L015692/1 STOR-i Centre for Doctoral Training. This work was supported by EPSRC grant EP/K014463/1, ONR Grant N00014-15-1-2380 and NSF CAREER Award IIS-1350133.

Appendix A Computational cost proofs

A.1 Proof of Proposition 1

Proof.

Let $\pi$ be the invariant distribution of the underlying dynamics, so that it has density $e^{-f(\theta)}=p(\theta|\mathbf{x})$, and define $W_{2}(\nu_{k},\pi)$ to be the Wasserstein distance between $\nu_{k}$ and $\pi$. Define $\xi_{k}$ to be the SGLD-CV gradient noise term. Then we can write a single step of SGLD-CV as

[TABLE]

We have that $\theta_{k}\sim\nu_{k}$, and follow similarly to the proof of Dalalyan and Karagulyan (2017, Proposition 2). First define $Y_{0}$ to be a draw from the invariant distribution $\pi$, such that the joint distribution of $Y_{0}$ and $\theta_{k}$ minimizes $\mathbb{E}\left\lVert Y_{0}-\theta_{k}\right\rVert^{2}$. Here $\left\lVert.\right\rVert$ denotes the Euclidean distance for $\mathbb{R}^{d}$. It follows that $\mathbb{E}\left\lVert Y_{0}-\theta_{k}\right\rVert^{2}=W_{2}^{2}(\nu_{k},\pi)$.

Let $B_{t}$ be a $d$-dimensional Wiener process, independent of $\theta_{k}$, $Y_{0}$ and $\xi_{k}$ but which we couple to the injected noise $\zeta_{k}$ so that $B_{h}=\sqrt{h}\zeta_{k}$. Now let $Y_{t}$, $t>0$, follow the diffusion

[TABLE]

Let $\Delta_{k}=Y_{0}-\theta_{k}$ and $\Delta_{k+1}=Y_{h}-\theta_{k+1}$. Since we started the process $Y_{t}$ from $Y_{0}\sim\pi$, then it follows that $Y_{t}\sim\pi$ for all $t>0$. Also since $W_{2}^{2}(\nu_{k+1},\pi)$ minimizes the expected squared distance between two random variables with marginals $\nu_{k+1}$ and $\pi$ then it follows that $W_{2}^{2}(\nu_{k+1},\pi)\leq\mathbb{E}\left\lVert\Delta_{k+1}\right\rVert^{2}$.

Let us define

[TABLE]

Then by the unbiasedness of the gradient estimation, $\xi_{k}$ has mean 0 regardless of the value of $\theta_{k}$. Thus

[TABLE]

We can then apply Lemmas 2 and 4 in Dalalyan and Karagulyan (2017), stated below in Lemmas 4 and 5, as well as applying the gradient noise bound in Lemma 1, to obtain a bound on $W_{2}^{2}(\nu_{k+1},\pi)$ given $W_{2}^{2}(\nu_{k},\pi)$.

Lemma 4**.**

With $U$ as defined in (17), if $h<2m/(2M^{2}+m^{2})$, then $\left\lVert\Delta_{k}-hU\right\rVert\leq(1-mh)\left\lVert\Delta_{k}\right\rVert$.

The original lemma by Dalalyan and Karagulyan (2017) assumed $h<2/(m+M)$, but this holds when $h<2m/(2M^{2}+m^{2})$ as $m\leq M$.

Lemma 5**.**

Under Assumption 1. Let $V$ be as defined in (18), then

[TABLE]

Finally we can apply Lemma 1, as stated in the main body, to get

[TABLE]

Using Theorem 1 of Durmus and Moulines (2016)

[TABLE]

It follows that

[TABLE]

Now using that $W_{2}^{2}(\nu_{k+1},\pi)\leq\mathbb{E}\left\lVert\Delta_{k+1}\right\rVert^{2}$ we get the following

[TABLE]

where $\alpha=7\sqrt{2}/6$. Gathering like terms we can further bound $W_{2}^{2}(\nu_{k+1},\pi)$ to get the following recursive formula

[TABLE]

where

[TABLE]

We can now apply Lemma 1 of Dalalyan and Karagulyan (2017), as stated below to solve this recurrence relation.

Lemma 6**.**

Let $A$, $B$ and $C$ be non-negative numbers such that $A\in(0,1)$. Assume that the sequence of non-negative numbers $x_{k}$, $k=0,1,\dots$, satisfies the recursive inequality

[TABLE]

for every integer $k>0$. Then for all integers $k\geq 0$

[TABLE]

To complete the proof all that remains is to check $A\in(0,1)$ so that Lemma 6 can be applied. Clearly $A<1$, since $n\geq 1$ we have

[TABLE]

and the RHS is positive when $h\in(0,2m/(2M^{2}+m^{2}))$. ∎∎

A.2 Proof of Theorem 1

Proof.

Starting from Proposition 1, we have that

[TABLE]

where

[TABLE]

Suppose we stop the algorithm at iteration $K$. Using (21), the following are sufficient conditions that ensure $W_{2}^{2}(\nu_{K},\pi)<\epsilon_{0}/\sqrt{m}$,

[TABLE]

The starting point $\theta_{0}$ is deterministic, so from Theorem 1 of Durmus and Moulines (2016)

[TABLE]

If we rewrite

[TABLE]

where $\gamma\in(0,1)$ is some constant and $R:=M/m$ as defined in the theorem statement, then it follows that we can write

[TABLE]

Since we have the condition

[TABLE]

then $\gamma\leq\frac{1}{2}$.

Now suppose, using (27), we set

[TABLE]

Then using the result for the deterministic starting point $\theta_{0}$ (25), we find that (28) implies that

[TABLE]

Using (27) and that our conditions imply $\gamma<1/2$. Hence (22) holds.

Using that for some real number $y\in[0,1]$, $\sqrt{1-y}\leq 1-y/2$, we can bound $A$ by

[TABLE]

As $\gamma\leq 1/2$, for (23) to hold it is sufficient that

[TABLE]

where $C/A_{0}\geq 2\alpha M\sqrt{hd}/m$. This leads to the following sufficient condition on $h$,

[TABLE]

Similarly for (24) it is sufficient that

[TABLE]

Now

[TABLE]

Leading to the following sufficient condition on $n$

[TABLE]

Now due to the conditions on $h$, define

[TABLE]

Then (24) will hold when

[TABLE]

∎∎

A.3 Proof of Lemma 1

Proof.

Our proof follows similarly to Dubey et al. (2016),

[TABLE]

Where the third line follows due to independence. For any random variable $R$, we have that $\mathbb{E}\left\lVert R-\mathbb{E}R\right\rVert^{2}\leq\mathbb{E}\left\lVert R\right\rVert^{2}$. Using this, the Lipschitz results of Assumption 2 and our choice of $p_{i}$, gives the following, where $\mathbb{E}_{I}$ refers to expectation with respect to the sampled datum index, $I$,

[TABLE]

from which the required bound follows trivially. ∎∎

A.4 Proof of Lemma 2

Proof.

Lipschitz condition: By the triangle inequality

[TABLE]

Strong convexity: We have that

[TABLE]

∎∎

Appendix B Postprocessing proofs

B.1 Proof of Theorem 2

Proof.

We start from the bound in Theorem $6.1$ of Mira et al. (2013), stating for some control variate $h$, the optimal variance reduction $R$ is given by

[TABLE]

so that in our case we have

[TABLE]

Then we can apply Lemmas 7, 8, defined in Section B.2, to get the desired result

[TABLE]

∎∎

B.2 Lemmas

Lemma 7**.**

Define $A=\sum_{i=1}^{d}a_{i}^{2}$, and let $\xi_{S}(\theta)=\widehat{\nabla\log p(\theta|\mathbf{x})}-\nabla\log p(\theta|\mathbf{x})$ be the noise in the gradient estimate. Then

[TABLE]

Proof.

We can condition on the gradient noise, and then immediately apply the Cauchy-Schwarz inequality to get

[TABLE]

∎∎

Lemma 8**.**

Under Assumption 4, define $A=\sum_{i=1}^{d}a_{i}^{2}$. Then $\mathbb{E}_{\theta|x}\left[h(\theta)\right]^{2}\leq A\sigma(N+1)/4$.

Proof.

Applying the Cauchy-Schwarz inequality

[TABLE]

∎∎