Consistency and fluctuations for stochastic gradient Langevin dynamics

Introduction

We are entering the age of Big Data, where significant advances across a range of scientific, engineering and societal pursuits hinge upon the gain in understanding derived from the analyses of large scale data sets. Examples include recent advances in genome-wide association studies (Hirschhorn and Daly, 2005; McCarthy et al., 2008; Wang et al., 2005), speech recognition (Hinton et al., 2012), object recognition (Krizhevsky et al., 2012), and self-driving cars (Thrun, 2010). As the quantity of data available has been outpacing the computational resources available in recent years, there is an increasing demand for new scalable learning methods, for example methods based on stochastic optimization (Robbins and Monro, 1951b; Srebro and Tewari, 2010; Sato, 2001; Hoffman et al., 2010), distributed computational architectures (Ahmed et al., 2012; Neiswanger et al., 2013; Minsker et al., 2014), greedy optimization (Harchaoui and Jaggi, 2014), as well as the development of specialized computing systems supporting large scale machine learning applications (Gonzalez, 2014).

Recently, there has also been increasing interest in methods for Bayesian inference scalable to Big Data settings. Rather than attempting a single point estimate of parameters typical in optimization-based or maximum likelihood settings, Bayesian methods attempt to obtain characterizations of the full posterior distribution over the unknown parameters and latent variables in the model, hence providing better characterizations of the uncertainties inherent in the learning process, as well as providing protection against overfitting. Scalable Bayesian methods proposed in the recent literature include stochastic variational inference (Sato, 2001; Hoffman et al., 2010), which applies stochastic approximation techniques to optimizing a variational approximation to the posterior, parallelized Monte Carlo (Neiswanger et al., 2013; Minsker et al., 2014), which distributes the computations needed for Monte Carlo sampling across a large compute cluster, as well as subsampling-based Monte Carlo (Welling and Teh, 2011; Ahn et al., 2012; Korattikara et al., 2014), which attempt to reduce the computational complexity of Markov chain Monte Carlo (MCMC) methods by applying updates to small subsets of data.

In this paper we study the asymptotic properties of the stochastic gradient Langevin dynamics (SGLD) algorithm first proposed by Welling and Teh (2011). SGLD is a subsampling-based MCMC algorithm based on combining ideas from stochastic optimization, specifically using small subsets of data to estimate gradients, with Langevin dynamics, a MCMC method making use of gradient information to produce better parameter updates. Welling and Teh (2011) demonstrated that SGLD works well on a variety of models and this has since been extended by Ahn et al. (2012, 2014) and Patterson and Teh (2013b).

The stochastic gradients in SGLD introduce approximations into the Markov chain, whose effect has to be controlled by using a slowly decreasing sequence of step sizes. Welling and Teh (2011) provided an intuitive argument that as the step-size decreases the variations introduced by the stochastic gradients gets dominated by the natural stochasticity of Langevin dynamics, the result being that the stochastic gradient approximation should wash out asymptotically and that the Markov chain should converge to the true posterior distribution.

In this paper, we make this intuitive argument more precise by providing conditions under which SGLD converges to the targeted posterior distribution; we describe a number of characterizations of this convergence. Specifically, we show that estimators derived from SGLD are consistent (Theorem 7) and satisfy a central limit theorem (CLT) (Theorem 8); the bias-variance trade-off of the algorithm is discussed in details in Section 5. In Section 6 we prove that, when observed on the right (inhomogeneous) time scale, the sample path of the algorithm converges to a Langevin diffusion (Theorem 9).

Our analysis reveals that for a sequence of step-sizes with algebraic decay $\delta_{m}\asymp m^{-\alpha}$ the optimal choice, when measured in terms of rate of decay of the mean squared error (MSE), is given for $\alpha_{\star}=1/3$; the choice $\delta_{m}\asymp m^{-\alpha_{\star}}$ leads to an algorithm that converges at rate $\mathcal{O}(m^{-1/3})$. This rate of convergence is worse than the standard Monte-Carlo $m^{-1/2}$-rate of convergence. This is not due to the stochastic gradients used in SGLD, but rather to the decreasing step-sizes.

These results are asymptotic in the sense that they characterise the behaviour of the algorithm as the number of steps approaches infinity. Therefore they do not necessarily translate into any insight into the behaviour for finite computational budgets which is the regime in which the SGLD might provide computational gains over alternatives. The mathematical framework described in this article show that the SGLD is a sound algorithm, an important result that has been missing in the literature.

SJV and YWT acknowledge EPSRC for research funding through grant EP/K009850/1 and EP/K009362/1. AHT is grateful for financial support in carrying out this research from a Singaporean MoE grant.

Stochastic Gradient Langevin Dynamics

where $\Delta\varphi\overset{{\mbox{\tiny def}}}{=}\sum_{i=1}^{d}\nabla^{2}_{i}\varphi$ denotes the standard Laplacian operator. The motivation behind the choice of Langevin diffusions is that, under certain conditions, they are ergodic with respect to the distribution $\pi$; for example, (Roberts and Tweedie, 1996; Stramer and Tweedie, 1999a, b; Mattingly et al., 2002) describe drift conditions of the type described in Section 3.2 that ensure that the total variation distance from stationarity of the law at time $t$ of the Langevin diffusion (1) decreases to zero exponentially quickly as $t\to\infty$.

Given a time-step $\delta>0$ and a current position $\theta_{t}$, it is often straightforward to simulate a random variable $\theta_{\star}$ that is approximately distributed as the law of $\theta_{t+\delta}$ given $\theta_{t}$. For stochastic differential equations, the Euler-Maruyama scheme (Maruyama, 1955) might be the simplest approach for approximating the law of $\theta_{t+\delta}$. For a Langevin diffusion this reads

for a standard $d$-dimensional centred Gaussian random variable $\eta$. To fully correct the discretization error, one can adopt a Metropolis-Hastings accept-reject mechanism. The resulting algorithm is usually referred to as the Metropolis-Adjusted-Langevin algorithm (MALA) (Roberts and Tweedie, 1996). Other discretizations can be used as proposals. For example, the random walk Metropolis-Hastings algorithm uses the discretization of a standard Brownian motion as the proposal, while the Hamiltonian Monte Carlo (HMC) algorithm (Duane et al., 1987) is based on discretizations of an Hamiltonian system of differential equations. See the excellent review of Neal (2010) for further information.

In this paper, we shall consider the situation where the target $\pi$ is the density of the posterior distribution under a Bayesian model where there are $N\gg 1$ i.i.d. observations, the so called Big Data regime,

Here, both computing the gradient term $\nabla\log\pi(\theta_{t})$ and evaluating the Metropolis-Hastings acceptance ratio require a computational budget that scales unfeasibly as $\mathcal{O}(N)$. One approach is to use a standard random walk proposal instead of Langevin dynamics, and to efficiently approximating the Metropolis-Hastings accept-reject mechanism using only a subset of the data (Korattikara et al., 2014; Bardenet et al., 2014).

This paper is concerned with stochastic gradient Langevin dynamics (SGLD), an alternative approach proposed by Welling and Teh (2011). This follows the opposite route and chooses to completely avoid the computation of the Metropolis-Hastings ratio. By choosing a discretization of the Langevin diffusion (1) with a sufficiently small step-size $\delta\ll 1$, because the Langevin diffusion is ergodic with respect to $\pi$, the hope is that even if the Metropolis-Hastings accept-reject mechanism is completely avoided, the resulting Markov chain still has an invariant distribution that is close to $\pi$. Choosing a decreasing sequence of step-sizes $\delta_{m}\to 0$ should even allow us to converge to the exact posterior distribution. To further make this approach viable in large $N$ settings, the gradient term $\nabla\log\pi(\theta)$ can be further approximated using a subsampling strategy. For an integer $1\leq n\leq N$ and a random subset $\tau\overset{{\mbox{\tiny def}}}{=}(\tau_{1},\ldots,\tau_{n})$ of $[N]\equiv\{1,\ldots,N\}$ generated by sampling with or without replacement from $[N]$, the quantity

is an unbiased estimator of $\nabla\log\pi(\theta)$. Most importantly, this stochastic estimate can be computed with a computational budget that scales as $\mathcal{O}(n)$ with $n$ potentially much smaller than $N$. Indeed, the larger the quotient $n/N$, the smaller the variance of this estimate.

Stochastic gradient methods have a long history in optimisation and machine learning and are especially relevant in the large dataset regime considered in this article (Robbins and Monro, 1951a; Bottou, 2010; Hoffman et al., 2013). In this paper we will adopt a slightly more general framework and assume that one can compute an unbiased estimate $\widehat{\nabla\log\pi}(\theta,\mathcal{U})$ to the gradient $\nabla\log\pi(\theta)$, where $\mathcal{U}$ is an auxiliary random variable which contains all the randomness involved in constructing the estimate. Without loss of generality we may assume (although this is unnecessary) that $\mathcal{U}$ is uniform on $(0,1)$. The unbiasedness of the estimator $\widehat{\nabla\log\pi}(\theta,\mathcal{U})$ means that

for an i.i.d. sequence $\eta_{m}\sim\operatorname{N}(0,I_{d})$, and an independent and i.i.d. sequence $\mathcal{U}_{m}$ of auxiliary random variables. This is the equivalent of the Euler-Maruyama discretization (3) of the Langevin diffusion (1) with a decreasing sequence of step-sizes and a stochastic estimate to the gradient term. The analysis presented in this article assumes for simplicity that the initial position $\theta_{0}$ of the algorithm is deterministic; in the simulation study of Section 7, the algorithms are started at the MAP estimator. Indeed, more general situations could be analysed with similar arguments at the cost of slightly less transparent proofs. Note that the process $(\theta_{m})_{m\geq 0}$ is a non-homogeneous Markov chain, and many standard analysis techniques for homogeneous Markov chains do not apply.

with $T_{m}=\delta_{1}+\ldots+\delta_{m}$. The quantity $\pi_{m}(\varphi)$ thus approximates the ergodic average $T_{m}^{-1}\int_{0}^{T_{m}}\varphi(\theta_{t})\,dt$ between time zero and $t=T_{m}$. During the course of the proof of our fluctuation Theorem 8, we will need to consider more general averaging schemes than the one above. Instead, for a general positive sequence of weights $\omega=(\omega_{m})_{m\geq 1}$, we define the $\omega$-weighted sum

with $\Omega_{m}\overset{{\mbox{\tiny def}}}{=}\omega_{1}+\ldots+\omega_{m}$. Indeed, $\pi_{m}^{\omega}(\varphi)=\pi_{m}(\varphi)$ in the particular case $(\omega_{m})_{m\geq 1}=(\delta_{m})_{m\geq 1}$; we will consider the weight sequence $\omega=\{\delta_{m}^{2}\}_{m\geq 1}$ in the proof of Theorem 8.

In the rest of this paper, we will build a rigorous framework for understanding the properties of this SGLD algorithm, demonstrating that the heuristics and numerical evidences presented in Welling and Teh (2011) were indeed correct.

Assumptions and Stability Analysis

This section starts with the basics assumptions we will need for the asymptotic results to follow, and illustrates some of the potential stability issues that may occur, would the SGLD algorithm be applied without care.

Throughout this text, we assume that the sequence of step-sizes $\delta=(\delta_{m})_{m\geq 1}$ satisfies the following usual assumption.

The step-sizes $\delta=(\delta_{m})_{m\geq 1}$ form a decreasing sequence with

Indeed, this assumption is easily seen to also be necessary for the Law of Large Numbers of Section 4 to hold. Furthermore, we will need at several occasions to assume the following assumption on the oscillations of a sequence of step-sizes $(\omega_{m})_{m\geq 1}$.

The step-sizes sequence $(\omega_{m})_{m\geq 1}$ is such that $\omega_{m}\to 0$ and $\Omega_{m}\to\infty$ and

where $\Delta(\omega_{m}/\delta_{m})\overset{{\mbox{\tiny def}}}{=}\omega_{m+1}/\delta_{m+1}-\omega_{m}/\delta_{m}$.

Assumption 2 holds if $\delta=(\delta_{m})_{m\geq 1}$ satisfies Assumption (1) and the weights are defined as $\omega_{m}=\delta^{p}_{m}$, for some some exponent $p\geq 1$ small enough for $\Omega_{m}\to\infty$. This is because the first sum is less than \sum_{m\geq 1}\big{|}\Delta(\omega_{m}/\delta_{m})\big{|}/\Omega_{1}=\delta_{1}^{p-1}/\Omega_{1}, while the finiteness of the second sum can be seen as follows:

For any exponents $0<\alpha<1$ and $0<p<1/\alpha$ the sequences $\delta_{m}=(m_{0}+m)^{-\alpha}$ and $\omega_{m}=\delta_{m}^{p}$ satisfy both Assumption 1 and Assumption 2.

2 Stability

Unfortunately, stability of the continuous time Langevin diffusion does not always translate into good behaviour for its Euler-Maruyama discretization. For example, even if the drift term points towards the right direction in the sense that ${\left<\nabla\log\pi(\theta),\theta\right>}<0$ for every parameter $\theta$, it might happen that the magnitude of the drift term is too large so that the Euler-Maruyama discretization overshoots and becomes unstable. In a one dimensional setting, this would lead to a Markov chain that diverges in the sense that the sequence $(\theta_{m})_{m\geq 0}$ alternates between taking arbitrarily large positive and negative values. Lemma $6.3$ of (Mattingly et al., 2002) gives such an example with a target density $\pi(\theta)\propto\exp\{-\theta^{4}\}$. See also Theorem $3.2$ of (Roberts and Tweedie, 1996) for examples of the same flavours.

Guaranteeing stability of the Euler-Maruyama discretization requires stronger Lypanunov type conditions. At a heuristic level, one must ensure that the drift term $\nabla\log\pi(\theta)$ points towards the centre of the state space. In addition, the previous discussion indicates that one must also ensure that the magnitude of this drift term is not too large. The following assumptions satisfy both heuristics, and we will show are enough to guarantee that the SGLD algorithm is consistent, with asymptotically Gaussian fluctuations.

There exists an exponent $p_{H}\geq 2$ such that

This implies that $\operatorname{\mathbf{E}}\left[\|H(\theta,\mathcal{U})\|^{2p}\right]\lesssim V^{p}(\theta)$ for any exponent $0\leq p\leq p_{H}$.

Equation (12) ensures that on average the drift term $\widehat{\nabla\log\pi}(\theta)$ points towards the centre of the state space, while equations (10) and (11) provide control on the magnitude of the (stochastic) drift term. The drift condition (12) implies in particular that the Langevin diffusion (1) converges exponentially quickly towards the equilibrium distribution $\pi$ (Mattingly et al., 2002; Roberts and Tweedie, 1996). The proof of the Law of Large Numbers (LLN) and the Central Limit Theorem (CLT) both exploit the following Lemma.

Let the step-sizes $(\delta_{m})_{m\geq 1}$ satisfy Assumption 1 and suppose that the stability Assumptions 4 hold. For any exponent $0\leq p\leq p_{H}$ the following bounds hold almost surely,

Moreover, for any exponent $0\leq p\leq p_{H}$ we have $\pi(V^{p})<\infty$. If the sequence of weights $(\omega_{m})_{m\geq 1}$ satisfies Assumption 2 the following holds almost surely,

The technical proof can be found in Section B. The idea is to leverage condition (12) in order to establish that the function $V^{p}$ satisfies both discrete and continuous drift conditions.

3 Scope of the analysis

For a posterior density $\pi$ of the form (4) and the usual unbiased estimate to $\nabla\log\pi$ described in Equation (5), to establish that Equations (10) and (11) hold it suffices to verify that the prior density $p_{0}$ is such that $\left\|\nabla\log\operatorname{p}_{0}(\theta)\right\|^{2}\lesssim V(\theta)$ and that for any index $1\leq i\leq N$ the likelihood term $\operatorname{p}\left(\left.y_{i}\,\right|\theta\right)$ is such that

Indeed, in these circumstances, we have $\left\|H(\theta,\mathcal{U})\right\|^{2p_{H}}\lesssim\sum_{i=1}^{N}\,\left\|\nabla\log\operatorname{p}(y_{i}\mid\theta)\right\|^{2p_{H}}$. Several such examples are described in Section 7.

It is important to note that the drift Condition (12) typically does not hold for distributions with heavy tails such that $\nabla\log\pi(x)\to 0$ as $\left\|x\right\|\to\infty$ (Roberts and Tweedie, 1996). For example, the standard MALA algorithm is not geometrically ergodic when $\nabla\log\pi(x)$ converges to zero as $\left\|x\right\|\to\infty$ (Theorem $4.3$ of (Roberts and Tweedie, 1996)); indeed, the analysis of standard local-move MCMC algorithms when applied to target densities with heavy tails is delicate and typically necessitate other tools Stramer and Tweedie (1999b); Jarner and Roberts (2007); Kamatani (2014) than the approach based on drift conditions of the type (12). The analysis of the properties of the SGLD algorithm when applied to such heavy tail densities is out of the scope of this article. It is important to note that many more complex scenarios involving high-dimensionality, multi-modality, non-parametric settings where the complexity of the target distribution increases with the size of the data, or combination thereof, are examples of interesting and relevant situations where our analysis typically does not apply; analysing the SGLD algorithm when applied to these challenging target distributions is well out of the scope of this article.

Consistency

with a similar result for $\omega$-weighted empirical averages, under assumptions on the weight sequence $\omega$. The proofs of several results of this paper make use of the following elementary lemma.

Let ${\left(\Delta M_{k}\right)}_{k\geq 0}$ and ${\left(R_{k}\right)}_{k\geq 0}$ be two sequences of random variables adapted to a filtration ${\left(\mathcal{F}_{k}\right)}_{k\geq 0}$ and let ${\left(\Gamma_{k}\right)}_{k\geq 0}$ be an increasing sequence of positive real numbers. The limit

holds almost surely if the following two conditions are satisfied.

The process $M_{m}=\sum_{k\leq m}\Delta M_{k}$ is a martingale, i.e. $\operatorname{\mathbf{E}}\left[\left.\Delta M_{k}\,\right|\mathcal{F}_{k}\right]=0$ and

The sequence ${\left(R_{k}\right)}_{k\geq 0}$ is such that

The above lemma, whose proof can be found in the appendix A, is standard; Lamberton and Pages (2002) also follows this route to prove several of their results.

If in addition the sequence of weights $\{\omega_{m}\}_{m\geq 1}$ satisfies Assumption (2), a similar result holds almost surely for the $\omega$-weighted ergodic average:

Proof In the following, we write $\operatorname{\mathbf{E}}_{k}{\left[\,\cdot\,\right]}$ and $\operatorname{\mathbf{P}}_{k}{\left(\,\cdot\,\right)}$ to denote the conditional expectation $\operatorname{\mathbf{E}}\left[\left.\,\cdot\,\,\right|\theta_{k}\right]$ and conditional probability $\operatorname{\mathbf{P}}\left(\left.\,\cdot\,\,\right|\theta_{k}\right)$ respectively. We use the notation $\Delta\theta_{k}\overset{{\mbox{\tiny def}}}{=}(\theta_{k+1}-\theta_{k})$. Finally, for notational convenience, we only present the proof in the scalar case $d=1$, the multidimensional case being entirely similar. We will give a detailed proof of Equation (18) and then briefly describe how the more general Equation (19) can be proven using similar arguments. To prove Equation (18), we first show that the sequence $(\pi_{m})_{m\geq 1}$ almost surely converges weakly to $\pi$. Equation (18) is then proved in a second stage.

Weak convergence of $(\pi_{m})_{m\geq 1}$. To prove that almost surely the sequence $(\pi_{m})_{m\geq 1}$ converges weakly towards $\pi$ it suffices to prove that the sequence is almost surely weakly pre-compact and that any weakly convergent subsequence of $(\pi_{m})_{m\geq 0}$ necessarily (weakly) converges towards $\pi$. By Prokhorov’s Theorem (Billingsley, 1995) and Equation (13), because the Lyapunov function $V$ goes to infinity as $\|\theta\|\to\infty$, the sequence $(\pi_{m})_{m\geq 1}$ is almost surely weakly pre-compact. It thus remains to show that if a subsequence converges weakly to a probability measure $\pi_{\infty}$ then $\pi_{\infty}=\pi$.

To prove Equation (20) we use the following decomposition of $\pi_{m}(\mathcal{A}\varphi)$,

Let us prove that the first term of (21) converges almost surely to zero. The numerator is equal to the sum of $\sum_{k=1}^{m}\operatorname{\mathbf{E}}_{k-1}[\varphi(\theta_{k})]-\varphi(\theta_{k})$ and $\varphi(\theta_{m})-\varphi(\theta_{0})$. By boundedness of $\varphi$, the term ${\left\{\varphi(\theta_{m})-\varphi(\theta_{0})\right\}}/T_{m}$ converges almost surely to zero. By Lemma 6, to conclude is suffices to show that the martingale difference terms $\operatorname{\mathbf{E}}_{k-1}{\left[\varphi(\theta_{k})\right]}-\varphi(\theta_{k})$ are such that

Because $\varphi$ is Lipschitz, it suffices to prove that$\sum_{k\geq 1}\operatorname{\mathbf{E}}{\left(\left\|\theta_{k+1}-\theta_{k}\right\|^{2}\right)}/T_{k}^{2}$ is finite. The stability Assumption 4 and Lemma 5 imply that the supremum $\sup_{m}\,\operatorname{\mathbf{E}}{\left[V(\theta_{m})\right]}$ is finite. Since $\operatorname{\mathbf{E}}_{k}{\left[\left\|\theta_{k+1}-\theta_{k}\right\|^{2}\right]}\lesssim\delta^{2}_{k+1}\,V(\theta)+\delta_{k+1}$, it follows that $\operatorname{\mathbf{E}}{\left(\left\|\theta_{k+1}-\theta_{k}\right\|^{2}\right)}$ is less than a constant multiple of $\delta_{k+1}$. Under Assumption 1, because the telescoping sum $\sum_{k\geq 1}T^{-1}(k)-T^{-1}(k+1)$ is finite, the sum $\sum_{k\geq 1}\delta_{k}/T_{k}^{2}$ is finite. This concludes the proof that the first term in (21) converges almost surely to zero.

The second term of (21) equals \big{(}R_{0}+\ldots+R_{m-1}\big{)}/T_{m} with

We now show that there exists a constant $C$ such that the bound $|R_{k}|\leq C\,\delta_{k+1}^{3/2}$ holds for any $k\geq 0$. To do so, let $K>0$ be such that the support of the test function $\varphi$ is included in the compact set $\Omega=[-K,K]$. We examine two cases separately.

If $|\theta_{k}|>K+1$ then $\varphi(\theta_{k})=\mathcal{A}\varphi(\theta_{k})=0$ so that $|R_{k}|\leq\|\varphi\|_{\infty}\times\operatorname{\mathbf{P}}_{k}(\theta_{k+1}\in\Omega)$. Since $\theta_{k+1}-\theta_{k}={\left\{\frac{1}{2}\nabla\log\pi(\theta_{k})+H(\theta_{k},\mathcal{U})\right\}}\,\delta_{k+1}+\sqrt{\delta_{k+1}}\,\eta$ we have

If $|\theta_{k}|\leq K+1$, we decompose $R_{k}$ into two terms. A second order Taylor formula yields

Under Assumption 4, the quantities $[\nabla\log\pi(\theta_{k})]^{2}$ and $\operatorname{\mathbf{E}}[H^{2}(\theta_{k},\mathcal{U})]$ are upper bounded by a constant multiple of $V(\theta_{k})$. Since the function $\theta\mapsto\varphi^{{}^{\prime\prime}}(\theta)\,V(\theta)$ is globally bounded (because continuous with compact support) this shows that $R_{k,1}$ is less than a constant multiple of $\delta_{k+1}^{2}$. Since $|\theta_{k}|\leq K+1$, the bounds $\operatorname{\mathbf{E}}[H^{3}(\theta,\mathcal{U})]\lesssim V^{3/2}(\theta)$ and $\sup_{k\geq 0}\operatorname{\mathbf{E}}[V^{3/2}(\theta_{k})]<\infty$ (see Lemma 5) yield that $\operatorname{\mathbf{E}}_{k}|\Delta\theta_{k}|^{3}\leq 9\,\overline{C}\,(\delta_{k+1}^{3}+\delta_{k+1}^{3/2})\lesssim\delta_{k+1}^{3/2}$ with

Note that $\overline{C}$ is finite by Assumption 4 and Lemma 5.

We have thus proved that there is a constant $C$ such $\left|R_{k}\right|\leq C\,\delta_{k+1}^{3/2}$ for $k\geq 0$; it follows that the sum ${\left(R_{0}+\ldots+R_{m-1}\right)}/T_{m}$ is less than a constant multiple of ${\left(\delta_{1}^{3/2}+\ldots+\delta_{m}^{3/2}\right)}/T_{m}$. Under Assumption 1, this upper bound converges to zero as $m\to\infty$, hence the conclusion.

This ends the proof of the almost sure weak convergence of $\pi_{m}$ towards $\pi$.

Proof of Equation (18). By assumption we have $|\varphi(\theta)|\leq C_{p}\,V^{p}(\theta)$ for some constant $C_{p}>0$ and exponent $p<p_{H}/2$. To show that $\pi_{m}(\varphi)\to\pi(\varphi)$ almost surely, we will use Lemma 5 and the almost sure weak convergence, which guarantees that $\pi_{m}(\widetilde{\varphi})\to\pi(\widetilde{\varphi})$ for a continuous and bounded test function $\widetilde{\varphi}$.

For any $t>0$, the set $\Omega_{t}\overset{{\mbox{\tiny def}}}{=}\{\theta:V(\theta)\leq t\}$ is compact and Tietze’s extension theorem (Rudin, 1986, Theorem $20.4$) yields that there exists a continuous function $\widetilde{\varphi}_{t}$ with compact support that agrees with $\varphi$ on $\Omega_{t}$ and such that $\|\widetilde{\varphi}_{t}\|_{\infty}=\sup\{|\varphi(\theta)|:\theta\in\Omega_{t}\}$. We can indeed also assume that $|\widetilde{\varphi}_{t}(\theta)|\leq C_{p}\,V^{p}(\theta)$. Since Lemma 5 states that $\sup_{m}\pi_{m}(V^{p_{H}/2})$ is almost surely finite, it follows that

By the triangle inequality, we thus have,

On the right-hand-side, the term in the middle can be made arbitrarily small as $m\to\infty$ since $\pi_{m}$ converges weakly towards $\pi$, while the other two terms converges to zero as $t\to\infty$. This concludes the proof of Equation (18).

Proof of Equation (19). The approach is very similar to the proof of Equation (18) and for this reason we only highlight the main differences. The same argument shows that the sequence $\pi^{\omega}_{m}$ is tight and it suffices to show that $\pi^{\omega}_{\infty}(\mathcal{A}\varphi)=0$ for any weak limit $\pi^{\omega}_{\infty}$ of the sequence $(\pi^{\omega}_{m})_{m\geq 0}$ for obtaining the almost sure weak convergences of $(\pi^{\omega}_{m})_{m\geq 0}$ towards $\pi$. One can then upgrade this almost sure weak convergence to a Law of Large Numbers. To prove (19), we thus concentrate on proving that $\pi^{\omega}_{\infty}(\mathcal{A}\varphi)=0$. For a smooth and compactly supported test function $\varphi$ we use the decomposition $\pi^{\omega}_{m}(\mathcal{A}\varphi)=S_{1}(m)+S_{2}(m)+S_{3}(m)$ with

and prove that each term converges to zero almost surely. For $S_{1}(m)$, by Lemma 6 it suffices to show that $\sum_{k\geq 1}(\omega_{k}/\delta_{k})^{2}\,\operatorname{\mathbf{E}}{\left[{\left\{\operatorname{\mathbf{E}}_{k-1}[\varphi(\theta_{k})]-\varphi(\theta_{k})\right\}}^{2}\right]}/\Omega_{k}^{2}$ is finite. This follows from the bound \operatorname{\mathbf{E}}{\left[\big{(}\operatorname{\mathbf{E}}_{k-1}[\varphi(\theta_{k})]-\varphi(\theta_{k})\big{)}^{2}\right]}\lesssim\delta_{k} and the fact that $\sum_{m\geq 0}\omega_{m}^{2}/(\Omega_{m}^{2}\,\delta_{m})$ is finite. For $S_{2}(m)$, we can write it as

Because $\Omega_{m}\to\infty$, $(\omega_{m+1}/\delta_{m+1})/\Omega_{m}\to 0$ and $\varphi$ is bounded, one can concentrate on proving that $\Omega_{m}^{-1}\sum_{k=1}^{m}\varphi(\theta_{k})\,\Delta(\omega_{k}/\delta_{k})$ converges almost surely to zero. By Lemma 6, it suffices to verify that $\sum_{k\geq 1}\operatorname{\mathbf{E}}{\left[\left|\varphi(\theta_{k})\,\Delta(\omega_{k}/\delta_{k})\right|\right]}/\Omega_{k}$ is finite; this directly follows from the boundedness of $\varphi$ and Assumption 2. Finally, algebra shows that $S_{3}(m)=\Omega_{m}^{-1}\,\sum_{1}^{m}(\omega_{k}/\delta_{k})\,R_{k-1}$ with the quantity $R_{k}$ defined in Equation (22). It has been proved that there is a constant $C$ such that, almost surely, $|R_{k}|\leq C\,\delta_{k+1}^{3/2}$ for all $k\geq 0$. Since $\delta_{m}\to 0$, the rescaled sum $\Omega_{m}^{-1}\,\sum_{k\leq m}\omega_{k}\delta_{k}^{1/2}$ converges to zero as $m\to\infty$. It follows that $S_{3}(m)$ converges almost surely to zero.

Fluctuations, Bias-Variance Analysis, and Central Limit Theorem

is introduced so that the additive functional $\pi_{m}(\varphi)$ of the trajectory of the Markov process $\{\theta_{k}\}_{k\geq 0}$ can be expressed as the sum of a martingale and a remainder term. A central limit for martingales can then be invoked to describe the asymptotic behaviour of the fluctuations

where the random variables $\Theta\overset{\mathcal{D}}{\sim}\pi$ and $\mathcal{U}$ are independent.

Proof The proof follows the strategy described in Lamberton and Pages (2002), with the additional difficulty that only unbiased estimates of the drift term of the Langevin diffusion are available. We use the decomposition

A fifth order Taylor expansion and Equation (7) yields that

In the above, we have defined $\mathcal{C}^{(k)}_{n,i}\equiv{\left(2^{i}\,i!\,(n-i)!\right)}^{-1}\,\nabla^{n}h(\theta_{k})\widehat{\nabla\log\pi}(\theta_{k},\mathcal{U}_{k+1})^{i}\eta_{k+1}^{n-i}$; the quantity $\xi_{k}$ lies between $\theta_{k}$ and $\theta_{k+1}$. It follows from the expression (2) of the generator of the $\mathcal{A}$ of the Langevin diffusion (1) and decomposition (28) that $\pi_{m}(\varphi)-\pi(\varphi)=\mathscr{F}_{m}+\mathscr{B}_{m}+\mathscr{R}_{m}$ where the fluctuation and bias terms are given by

Remainder term: we start by proving that the term $\mathscr{R}_{m}$ is negligible. The term ${\left\{h(\theta_{m})-h(\theta_{0})\right\}}/T_{m}^{1/2}$ converges to zero in probability because $\left|h(\theta)\right|\lesssim V^{p_{H}}(\theta)$ and Lemma 5 shows that $\sup_{m\geq 0}\,\operatorname{\mathbf{E}}[V^{p_{H}}(\theta_{m})]$ is almost surely finite. Similarly, Assumptions 1 and 4 and Lemma 5 yield that

from which it follows that ${\left\{\sum_{k=0}^{m-1}\nabla^{5}h(\xi_{k})\,{\left(\theta_{k+1}-\theta_{k}\right)}^{5}\right\}}/{\left\{\sum_{k=0}^{m-1}\delta_{k+1}^{2}\right\}}$ converges to zero in probability; we have exploited the fact that $\nabla^{5}h$ is assumed to be globally bounded. Essentially the same argument yield that the high-order terms are asymptotically negligible: for $(n,i)\in\mathcal{I}_{\mathscr{R},p}$ and $p\in\{5,6,7,8\}$ the limit

holds in probability because the coefficients $\mathcal{C}^{(k)}_{n,i}$ are uniformly bounded in expectation and the quantity ${\left(\sum_{k=0}^{m-1}\delta_{k+1}^{(n+i)/2}\right)}/{\left(\sum_{k=0}^{m-1}\delta_{k+1}^{2}\right)}$ converges to zero since $(n+i)/2\geq 5/2$ and $\delta_{k}\to 0$. To conclude, one needs to verify that the low order terms are also negligible in the sense that the limit

holds in probability with $X^{(k)}_{1,1}=\nabla h(\theta_{k})\,H(\theta_{k},\mathcal{U}_{k+1})$ and $X^{(k)}_{2,0}=\nabla^{2}h(\theta_{k})(\eta_{k+1}^{2}-1)$ and $X^{(k)}_{2,1}=-\mathcal{C}^{(k)}_{2,1}$ and $X^{(k)}_{3,0}=-\mathcal{C}^{(k)}_{3,0}$. Since $\operatorname{\mathbf{E}}\left[\left.X^{(k)}_{n,i}\,\right|\mathcal{F}_{k}\right]=0$ where $\mathcal{F}_{k}=\sigma{\left(\theta_{0},\ldots,\theta_{k}\right)}$ is the natural filtration associated to the process ${\left(\theta_{k}\right)}_{k\geq 0}$ it follows that

We made use of the fact that the expectations $\operatorname{\mathbf{E}}{\left[(X^{(k)}_{n,i})^{2}\right]}$ are uniformly bounded for all $k\geq 0$ by the same arguments as above, and that the final expression converges to 0 since $n+i\geq 2$, $\delta_{m}\to 0$ and $T_{m}\to\infty$. This concludes the proof that the remainder term $\mathscr{R}_{m}$ is asymptotically negligible.

Fluctuation term: we now prove that the fluctuations term converges in distribution at Monte-Carlo rate towards a Gaussian distribution,

Using the standard martingale central limit theorem (e.g. Theorem $3.2$, Chapter $3$ of (Hall and Heyde, 1980)), it suffices to verify that for any $\varepsilon>0$ the following limits hold in probability,

Since $T_{m}^{-2}\,\sum_{k=0}^{m-1}\delta_{k+1}^{2}\to 0$, the conclusion follows.

Bias term: we conclude by proving that the bias term is such that the limit

holds in probability. The quantity $\mathscr{B}_{m}/{\left(\sum_{k=1}^{m}\delta_{k}^{2}/T_{m}\right)}^{-1}$ can also be expressed as

for a martingale difference term $\Delta M_{k}\equiv{\left(\mathcal{C}^{(k)}_{2,2}+\mathcal{C}^{(k)}_{3,1}+\mathcal{C}^{(k)}_{4,0}\right)}-\Psi(\theta_{k})$ where $\Psi(\theta_{k})\equiv\operatorname{\mathbf{E}}\left[\left.\mathcal{C}^{(k)}_{2,2}+\mathcal{C}^{(k)}_{3,1}+\mathcal{C}^{(k)}_{4,0}\,\right|\mathcal{F}_{k}\right]$ and ${\left(\mathcal{C}^{(k)}_{2,2}+\mathcal{C}^{(k)}_{3,1}+\mathcal{C}^{(k)}_{4,0}\right)}$ equals

Under the assumptions of Theorem 8, the function $\Psi$ satisfies the hypothesis of Theorem 7 applied to the weight sequence $\{\delta_{k}^{2}\}_{k\geq 0}$; it follows that the first term in Equation (31) converge almost surely to $\mu(\varphi)$. It remains to prove that the second term in Equation (31) also converges almost surely to zero. By Lemma 6, it suffices to prove that the martingale

is bounded in $L^{2}$. Under the Assumption of Theorem 8, Lemma 5 yields that the martingale difference term $\Delta M_{k}$ is uniformly bounded in $L^{2}$ from which the conclusion readily follows.

For the standard choice of step-sizes $\delta_{m}=(m_{0}+m)^{-\alpha}$ the statistical fluctuations dominate in the range $1/3<\alpha\leq 1$, there is an exact balance between bias and fluctuations for $\alpha=1/3$, and the bias dominates for $0<\alpha<1/3$. The optimal rate of convergence is obtained for $\alpha=1/3$ and leads to an algorithm that converges at rate $m^{-1/3}$.

Diffusion limit

In this section we show that, when observed on the right (inhomogeneous) time scale, the sample path of the SGLD algorithm converges to the continuous time Langevin diffusion of Equation (1), confirming the heuristic discussion in Welling and Teh (2011).

If the drift function $\theta\mapsto\frac{1}{2}\nabla\log\pi(\theta)$ is globally Lipschitz, then the Itô’s map $\mathcal{I}$ is well defined and continuous. Further, the image $\mathcal{I}(W)$ under the Itô map of a standard Brownian motion $W$ on $[0,T]$ can be seen to be described by Langevin diffusion (1).

Define $T^{(r)}_{0}=0$ and $T^{(r)}_{k}=\delta^{(r)}_{1}+\cdots+\delta^{(r)}_{k}$ for each $k\geq 1$. The Markov chains are coupled to $W$ as follows:

for an i.i.d. collection of auxiliary random variables $(\mathcal{U}^{(r)}_{k})_{r\geq 1,k\geq 1}$. Note that $(\eta^{(r)}_{k})_{k\geq 1}$ form an i.i.d. sequence of $\operatorname{N}(0,1)$ variables for each $r$. We can construct piecewise affine continuous time sample paths $(S^{(r)})_{r\geq 1}$ by linearly interpolating the Markov chains,

for $x\in$. The approach then amounts to showing that each $S^{(r)}$ can be expressed as $\mathcal{I}(\widetilde{W}^{(r)})+e^{(r)}$, where $\widetilde{W}^{(r)}$ is a sequence of stochastic processes converging to $W$ and $e^{(r)}$ is asymptotically negligible, and making use of the continuity properties of the Itô map $\mathcal{I}$.

For convenience, we define $\widetilde{W}^{(r)}$ as the continuous piecewise affine processes that satisfies $\widetilde{W}^{(r)}(T^{(r)}_{k})=W(T^{(r)}_{k})$ for all $0\leq k\leq m(r)$ and that is affine in between. It follows that for any time $T^{(r)}_{k-1}\leq t\leq T^{(r)}_{k}$ we have

are finite, it readily follows that $\|e^{(r)}_{1}\|_{\infty}$ converges to zero in expectation.

Numerical Illustrations

In this section we illustrate the use of the SGLD method to a simple Gaussian toy model and to a Bayesian logistic regression problem. We verify that both models satisfy Assumption 4, the main assumption needed for our asymptotic results to hold. Simulations are then performed to empirically confirm our theory; for step-sizes sequences of the type $\delta_{m}=(m_{0}+m)^{-\alpha}$, both the rate of decay of the MSE and the impact of the sub-sampling scheme are investigated. The main purpose of this article is to establish the missing theoretical foundation of stochastic gradient methods for the approximation of expectations. For more exhaustive simulation studies we refer to Welling and Teh (2011); S. Ahn and Welling (2012); Patterson and Teh (2013a); Chen et al. (2014). By considering a logistic regression model, we demonstrate that the SGLD can be advantageous over the Metropolis-Adjusted-Langevin (MALA) algorithm if the available computational budget only allows a few iterations through the whole data set, see Section 7.2.2.

Consider $N$ independent and identically distributed observations $(x_{i})_{i=1}^{N}$ from the two parameters location model given by

We use a Gaussian prior $\theta\sim\operatorname{N}(0,\sigma_{\theta}^{2})$ and assume that the variance hyper-parameters $\sigma_{\theta}^{2}$ and $\sigma_{x}^{2}$ are both known. The posterior density $\pi(\theta)$ is normally distributed with mean $\mu_{p}$ and variance $\sigma_{p}^{2}$ given by

where $\bar{x}=(x_{1}+\ldots+x_{N})/N$ is the sample average of the observations. In this case, we have

for a random subset $\mathcal{I}_{n}(\mathcal{U})\subset[N]$ of cardinal $n$.

We verify in this section that Assumption (4) is satisfied for the following choice of Lyapunov function,

Since the error term $H(\theta,\mathcal{U})$ is globally bounded, the drift $(1/2)\nabla\log\pi$ and the Lyapunov function $V$ are linear, Assumptions (4).1 and (4).2 are satisfied. Finally, to verify Assumption (4).3, it suffices to note that since $\nabla\log\pi(\theta)=-(\theta-\mu_{p})/\sigma_{p}^{2}$ we have

In other words, Assumption (4).3 holds with $\alpha=\beta=1/\sigma_{p}^{2}$.

1.2 Simulations

We chose $\sigma_{\theta}=1$, $\sigma_{x}=5$ and created a data set consisting of $N=100$ data points simulated from the model. We used $n=10$ as the size of subsets used to estimate the gradients. We evaluated the convergence behaviour of SGLD using the test function $\mathcal{A}\varphi$ where $\varphi=\sin\left(x-\mu_{p}-0.5\sigma_{p}\right)$.

We are interested in confirming the asymptotic convergence regimes of Theorem 8 by running SGLD with a range of step sizes, and plotting the mean squared error (MSE) achieved by the estimate $\pi_{m}(\mathcal{A}\varphi)$ against the number of steps $m$ of the algorithm to determine the rates of convergence. We used step sizes $\delta_{m}=(m+m_{0}(\alpha))^{-\alpha}$, for $\alpha\in\{0.1,0.2,0.3,0.33,0.4,0.5\}$ where $m_{0}(\alpha)$ is chosen such that $\delta_{1}$ is less than the posterior standard deviation. According to the Theorem, the MSE should scale as $T_{m}^{-1}$ for $\alpha>1/3$, and $\sum_{k=1}^{m}\delta_{k}^{2}/T_{m}$ for $\alpha\leq 1/3$.

The observed MSE is plotted against $m$ on a log-log plot in Figure 1. As predicted by the theory, the optimal rate of decay is around $\alpha_{\star}=1/3$. To be more precise, we estimate the rates of decay by estimating the slopes on the log-log plots. This is plotted in Figure 2, which also shows a good match to the theoretical rates given in Theorem 8, where the best rate of decay is $2/3$ achieved at $\alpha=1/3$. Finally, to demonstrate that there are indeed two distinct regimes of convergence, in Figure 3 we have plotted the MSE multiplied by $T_{m}$. For $\alpha>1/3$, the plots remain flat, showing that the MSE does indeed decay as $T_{m}^{-1}$. For $\alpha<1/3$, the plots diverge, showing that the MSE decays at a slower rate than $T_{m}^{-1}$.

For $\alpha=0.33$, Figure 4 depicts how the MSE decreases as a function of the number of likelihood evaluations for subsample sizes $n=1,5,10,50,100$.

2 Logistic Regression

We verify in this section that Assumption (4) is satisfied for the following logistic regression model. Consider $N$ independent and identically observations $(y_{i})_{i=1}^{N}$ distributed as

for a random subset $\mathcal{I}_{n}(\mathcal{U})\subset[N]$ of cardinal $n$.

We verify in this section that Assumption (4) is satisfied for the Lyapunov function $V(\theta)=1+\|\theta\|^{2}$. Since $H(\theta,\mathcal{U})$ is globally bounded and $\|\nabla V(\theta)\|^{2}=\|\theta\|^{2}$ and

it is straightforward to see that Assumption (4).1 and (4).2 are satisfied. Finally,

2.2 Comparison of the SGLD and the MALA for logistic regression

We consider a simulated dataset where $d=3$ and $N=1000$. We set the input covariates $x_{i}=(x_{i,1},x_{i,2},1)$ with $x_{i,1},x_{i,2}\overset{\text{i.i.d.}}{\sim}\operatorname{N}(0,1)$ for $i=1\dots N$, and use a Gaussian prior $\theta\sim\operatorname{N}(0,I)$. We draw a $\theta_{0}\sim\operatorname{N}(0,I)$ and based on it we generate $y_{i}$ according to the model probabilities (36). In the following we compare MALA in SGLD by comparing their estimate for the variance of the first component.

The findings of this article show that SGLD-based expectation estimates converge at a slower rate of at most $n^{-\frac{1}{3}}$ compared to the standard rate of $n^{-\frac{1}{2}}$ for standard MCMC algorithms such as the MALA algorithm. In the following we demonstrate that in the non-asymptotic regime (allowing only a few passes through the data set) the SGLD can be advantageous. We start both algorithms at the MAP estimator and we ensure that this study is not biased due to different speeds in finding the mode of the posterior. For a fair comparison we tune the MALA to an acceptance rate of approximately $0.564$ following the findings of Roberts and Rosenthal (1998). For the SGLD-based variance estimate of the first component for $n=30$ we choose $\delta_{m}=(a\cdot m+b)^{-0.38}$ as step sizes and optimise over the choices of $a$ and $b$. This is achieved by estimating the MSE for choices of $a$ and $b$ on a log-scale grid based on $512$ independent runs. The estimates based on $20$ and $1000$ effective iterations through the data set the averages are visualised in the heat maps in Figure 5. That means we limit the algorithm to $200$ and $1000000$ likelihood evaluations, respectively. The figures indicate that the range of the good parameter choices seems to be the same in both cases. Using the heat map for the estimated MSE after 20 iterations through the data set, we pick $a=5.89\cdot 10^{7}$ and $b=7.90\cdot 10^{8}$ and compare the time behaviour of the SGLD and the MALA algorithm in Figure 6. The figure is a simulation evidence that the SGLD algorithm can be advantageous in the initial phase for the first few iterations through the data set. This recommends further investigation as the initial phase can be quite different from the asymptotic phase.

Conclusion

The CLT and bias-variance decomposition can be leveraged to show that it is optimal to choose a step-sizes sequences $(\delta_{m})_{m\geq 0}$ that scales as $\delta_{m}\asymp m^{-1/3}$; the resulting algorithm converges at rate $m^{-1/3}$. Note that this recommendation is different from the previously suggested Welling and Teh (2011) choice of $\delta_{m}\asymp m^{-1/2}$.

Our theory suggests that an optimally tuned SGLD method converges at rate $\mathcal{O}(m^{-1/3})$, and is thus asymptotically less efficient than a standard MCMC procedure. We believe that this result does not necessarily preclude SGLD to be more efficient in the initial transient phase, a result hinted at in Figure 4; the detailed study of this (non-asymptotic) phenomenon is an interesting venue of research. The asymptotic convergence rate of SGLD depends crucially on the decreasing step sizes, which is required to reduce the effect of the discretization bias due to the lack of a Metropolis-Hastings correction. Another avenue of exploration is to determine more precisely the bias resulting from the discretization of the Langevin diffusion, and to study the effect of the choice of step sizes in terms of the trade-off between bias, variance, and computation.

A Proof of Lemma 6

Recall Kronecker’s Lemma (Shiryaev, 1996, Lemma IV.3.2) that states that for a non-decreasing and positive sequence $b_{m}\to\infty$ and another real valued sequence $(a_{m})_{m\geq 0}$ such that the series $\sum_{m\geq 0}a_{m}/b_{m}$ converges the following limit holds,

For proving Equation (15) it thus suffices to show that the sums $\sum_{k\geq 0}\left|\Delta M_{k}\right|/T_{k}$ and $\sum_{k\geq 0}\left|X_{k}\right|/T_{k}$ are almost surely finite. This follows from Condition (16) ($L^{2}$ martingale convergence theorem) and Condition (17).

B Proof of Lemma 5

For clarity, the proof is only presented in the scalar case $d=1$; the multidimensional setting is entirely similar. Before embarking on the proof, let us first mention some consequences of Assumptions 4 that will be repeatedly used in the sequel. Since the second derivative $V^{{}^{\prime\prime}}$ is globally bounded and $(V^{\prime})^{2}$ is upper bounded by a multiple of $V$, we have that

and that the function $V^{1/2}$ is globally Lipschitz. By expressing the quantity $V^{p}(\theta+\varepsilon)$ as \big{(}V^{1/2}(\theta)+[V^{1/2}(\theta+\varepsilon)-V^{1/2}(\theta)]\big{)}^{2p}, it then follows that

Similarly, Definition (7), the bound $\|\nabla\log\operatorname{p}(\theta)\|^{2}\lesssim V(\theta)$ and Equation (10) yield that for any exponent $0\leq p\leq p_{H}$ the following holds,

For clarity, the proof of Lemma (5) is separated into several steps. First, we establish that the process $m\mapsto V^{p}(\theta_{m})$ satisfies a Lyapunov type condition; see Equation (40) below. We then describe how Equation (13) follows from this Lyapunov condition. The fact that $\pi(V^{p})$ is finite can be seen as a consequence of Theorem $2.2$ of (Roberts and Tweedie, 1996).

Discrete Lyapunov condition. Let us prove that there exists an index $m_{0}\geq 0$ and constants $\alpha_{p},\beta_{p}>0$ such that for any $m\geq m_{0}$ we have

Since for any $\varepsilon$ there exists $C_{\varepsilon}$ such that $V^{p-1}(\theta)\leq C_{\varepsilon}+\,\varepsilon V^{p}(\theta)$, for proving (40) it actually suffices to verify that we have

for some constants $\widetilde{\alpha_{p}},\widetilde{\beta_{p}}>0$ and index $m\geq 1$ large enough. A second order Taylor expansion yields that the left hand side of (41) is less than

for a random quantity $\xi$ lying between $\theta_{m}$ and $\theta_{m+1}$. Since $\operatorname{\mathbf{E}}_{m}[\theta_{m+1}-\theta_{m}]=\frac{1}{2}\,\nabla\log\operatorname{p}(\theta_{m})$, the drift condition (12) yields that the first term of (42) is less than

for $\alpha,\beta>0$ given by Equation (12). Consequently, for proving Equation (40), it remains to bound the second term of (42). Equation (37) shows that $|(V^{p})^{{}^{\prime\prime}}(\xi)|$ is upper bounded by a multiple of $|V^{p-1}(\xi)|$; the bound (38) then yields that $|V^{p-1}(\xi)|$ is less than a constant multiple of $|V^{p-1}(\theta_{m})|+|\theta_{m+1}-\theta_{m}|^{2(p-1)}$. It follows from the bound (39) on the difference $(\theta_{m+1}-\theta_{m})$ and the assumption $\operatorname{\mathbf{E}}[\,\|H(\theta,\mathcal{U})\|^{2p_{H}}\,]\lesssim V^{p_{H}}(\theta)$ that for any $\varepsilon>0$ one can find an index $m_{0}\geq 1$ large enough such that for any index $m\geq m_{0}$ the second term of (41) is less than a constant multiple of

for a constant $\beta_{p,\varepsilon}>0$. Equations (43) and (44) directly yield to Equation (41), which in turn implies to Equation (40).

Proof that $\sup_{m\geq 1}\;\operatorname{\mathbf{E}}[V^{p}(\theta_{m})]\;<\;\infty$ for any $p\leq p_{H}$. Equations (38) and (39) show that if $\operatorname{\mathbf{E}}[V^{p}(\theta_{m})]$ is finite then so is $\operatorname{\mathbf{E}}[V^{p}(\theta_{m+1})]$. Under the conditions of Lemma 5, this shows that $\operatorname{\mathbf{E}}[V^{p}(\theta_{m})]$ is finite for any $m\geq 0$. An inductive argument based on the discrete Lyapunov Equation (40) then yields that for any index $m\geq m_{0}$ the expectation $\operatorname{\mathbf{E}}[V^{p}(\theta_{m})]$ is less than

It follows that $\sup_{m\geq 1}\;\operatorname{\mathbf{E}}[V^{p}(\theta_{m})]$ is finite.

Proof that $\sup_{m\geq 1}\;\pi_{m}(V^{p})\;<\;\infty$ for any $p\leq p_{H}/2$. One needs to prove that the sequence $(1/T_{m})\sum_{k=m_{0}}^{m}\delta_{k+1}V^{p}(\theta_{k})$ is almost surely bounded. The discrete Lyapunov Equation (40) yields that $\delta_{k+1}V^{p}(\theta_{k})$ is less than $\delta_{k+1}\,\beta_{p}/\alpha_{p}-\operatorname{\mathbf{E}}_{k}[V^{p}(\theta_{k+1})-V^{p}(\theta_{k})]/\alpha_{p}$; this yields that $(1/T_{m})\sum_{k=m_{0}}^{m}\delta_{k+1}V^{p}(\theta_{k})$ is less than a constant multiple of

To conclude the proof, we prove that the last term in the above displayed Equation almost surely converges to zero; by Lemma 6, it suffices to prove that the quantity

is almost surely finite. We have $\operatorname{\mathbf{E}}\left[|V^{p}(\theta_{k+1})-\operatorname{\mathbf{E}}_{k}[V^{p}(\theta_{k+1})]|^{2}\right]\leq 2\times\operatorname{\mathbf{E}}[|\,V^{p}(\theta_{k+1})-V^{p}(\theta_{k})\,|^{2}]$ and the mean value theorem yields that $|V^{p}(\theta_{k+1})-V^{p}(\theta_{k})|\lesssim V^{p-1}(\xi)\,V^{\prime}(\xi)\,(\theta_{k+1}-\theta_{k})$ for some $\xi$ lying between $\theta_{k}$ and $\theta_{k+1}$. The bound $|V^{\prime}(\theta)|\lesssim V^{1/2}(\theta)$ and Equation (38) then yield that |V^{p}(\theta_{k+1})-V^{p}(\theta_{k})|\lesssim V^{p-1/2}(\theta_{k})\,\big{|}\theta_{k+1}-\theta_{k}\big{|}+\big{|}\theta_{k+1}-\theta_{k}\big{|}^{2p}. From the bound (39) and the assumption that $\operatorname{\mathbf{E}}[H(\theta,\mathcal{U})^{2p_{H}}]\lesssim V^{p_{H}}(\theta)$ it follows that the quantity in Equation (46) is less than a constant multiple of

Since $\operatorname{\mathbf{E}}\left[V^{2p}(\theta_{k})\right]$ is uniformly bounded for any $p\leq p_{H}/2$ and $\sum_{m\geq m_{0}}\delta_{m}/T^{2}(m)<\infty$ (because the sum $\sum_{m}T^{-1}(m+1)-T^{-1}(m)$ is finite), the conclusion follows.

Proof of $\pi(V^{p})\;<\;\infty$ for any $p\geq 0$. Since $V(\theta)\lesssim 1+\|\theta\|^{2}$, the drift condition (12) yields that Theorem $2.1$ of (Roberts and Tweedie, 1996) holds. Moreover, the bound $V^{p-1}(\theta)\leq C_{\varepsilon}+\varepsilon\,V^{p}(\theta)$ implies that there are constants $\alpha_{p,*}.\beta_{p,*}>0$ such that

where $\mathcal{A}$ is the generator of the Langevin diffusion (1). Theorem $2.2$ of (Roberts and Tweedie, 1996) gives the conclusion.

Proof that $\sup_{m\geq 1}\;\pi^{\omega}_{m}(V^{p})\;<\;\infty$ for any $p\leq p_{H}/2$. One needs to prove that the sequence $[1/\Omega_{m}]\times\sum_{k=m_{0}}^{m}\omega_{k+1}V^{p}(\theta_{k})$ is almost surely bounded. The bound $\delta_{k+1}V^{p}(\theta_{k})\lesssim\delta_{k+1}\,\beta_{p}/\alpha_{p}-\operatorname{\mathbf{E}}_{k}[V^{p}(\theta_{k+1})-V^{p}(\theta_{k})]/\alpha_{p}$ yields that $\pi^{\omega}_{m}(V^{p})$ is less than a constant multiple of

To conclude the proof, we establish that the following limits hold almost surely,

To prove Equation (48) it suffices to use the assumption that $\sum_{m\geq 0}\omega^{2}_{m}/[\delta_{m}\Omega^{2}_{m}]<\infty$ and then follow the same approach used to establish that the quantity (46) is finite. Lemma 6 shows that to prove Equation (49) it suffices to verify that

This directly follows from the assumption that \sum_{m\geq 0}\,\big{|}\Delta(\omega_{m}/\delta_{m})\big{|}/\Omega_{m}<\infty and the fact that $\sup_{m\geq 0}\,\operatorname{\mathbf{E}}[V^{p}(\theta_{m})]$ is finite.

References