# On the Convergence of Stochastic Gradient MCMC Algorithms with High-Order Integrators

Recent advances in Bayesian learning with large-scale data have witnessed emergence of stochastic gradient MCMC algorithms (SG-MCMC), such as stochastic gradient Langevin dynamics (SGLD), stochastic gradient Hamiltonian MCMC (SGHMC), and the stochastic gradient thermostat. While finite-time convergence properties of the SGLD with a 1st-order Euler integrator have recently been studied, corresponding theory for general SG-MCMCs has not been explored. In this paper we consider general SG-MCMCs with high-order integrators, and develop theory to analyze finite-time convergence properties and their asymptotic invariant measures. Our theoretical results show faster convergence rates and more accurate invariant measures for SG-MCMCs with higher-order integrators. For example, with the proposed efficient 2nd-order symmetric splitting integrator, the mean square error (MSE) of the posterior average for the SGHMC achieves an optimal convergence rate of L^-4/5 at L iterations, compared to L^-2/3 for the SGHMC and SGLD with 1st-order Euler integrators. Furthermore, convergence results of decreasing-step-size SG-MCMCs are also developed, with the same convergence rates as their fixed-step-size counterparts for a specific decreasing sequence. Experiments on both synthetic and real datasets verify our theory, and show advantages of the proposed method in two large-scale real applications.

## Authors

• 69 publications
• 17 publications
• 169 publications
• ### (Non-) asymptotic properties of Stochastic Gradient Langevin Dynamics

Applying standard Markov chain Monte Carlo (MCMC) algorithms to large da...
01/02/2015 ∙ by Sebastian J. Vollmer, et al. ∙ 0

• ### High-Order Stochastic Gradient Thermostats for Bayesian Learning of Deep Models

Learning in deep models using Bayesian methods has generated significant...
12/23/2015 ∙ by Chunyuan Li, et al. ∙ 0

Stochastic gradient Hamiltonian Monte Carlo (SGHMC) is an efficient meth...
02/29/2020 ∙ by Ruqi Zhang, et al. ∙ 0

Stochastic gradient MCMC (SG-MCMC) has played an important role in large...
10/21/2016 ∙ by Changyou Chen, et al. ∙ 0

• ### A Convergence Analysis for A Class of Practical Variance-Reduction Stochastic Gradient MCMC

Stochastic gradient Markov Chain Monte Carlo (SG-MCMC) has been develope...
09/04/2017 ∙ by Changyou Chen, et al. ∙ 0

• ### On Convergence-Diagnostic based Step Sizes for Stochastic Gradient Descent

Constant step-size Stochastic Gradient Descent exhibits two phases: a tr...
07/01/2020 ∙ by Scott Pesme, et al. ∙ 0

• ### Bayesian Sparse learning with preconditioned stochastic gradient MCMC and its applications

In this work, we propose a Bayesian type sparse deep learning algorithm....
06/29/2020 ∙ by Yating Wang, et al. ∙ 18

##### This week in AI

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

## 1 Introduction

In large-scale Bayesian learning, diffusion based sampling methods have become increasingly popular. Most of these methods are based on Itô diffusions, defined as:

 dXt =F(Xt)dt+σ(Xt)dWt . (1)

Here represents model states, the time index, is Brownian motion, functions and ( not necessarily equal to ) are assumed to satisfy the usual Lipschitz continuity condition.

In a Bayesian setting, the goal is to design appropriate functions and , so that the stationary distribution, , of the Itô diffusion has a marginal distribution that is equal to the posterior distribution of interest. For example, 1st-order Langevin dynamics (LD) correspond to , and , with being the identity matrix; 2nd-order Langevin dynamics correspond to , , and for some . Here is the unnormalized negative log-posterior, and is known as the momentum [1, 2]. Based on the Fokker-Planck equation [3], the stationary distributions of these dynamics exist and their marginals over are equal to , the posterior distribution we are interested in.

Since Itô diffusions are continuous-time Markov processes, exact sampling is in general infeasible. As a result, the following two approximations have been introduced in the machine learning literature

[4, 1, 2], to make the sampling numerically feasible and practically scalable: 1) Instead of analytically integrating infinitesimal increments , numerical integration over small step is used to approximate the integration of the true dynamics. Although many numerical schemes have been studied in the SDE literature, in machine learning only the 1st-order Euler scheme is widely applied. 2) During every integration, instead of working with the gradient of the full negative log-posterior , a stochastic-gradient version of it, , is calculated from the -th minibatch of data, important when considering problems with massive data. In this paper, we call algorithms based on 1) and 2) SG-MCMC algorithms. To be complete, some recently proposed SG-MCMC algorithms are briefly reviewed in Appendix A. SG-MCMC algorithms often work well in practice, however some theoretical concerns about the convergence properties have been raised [5, 6, 7].

Recently, [8, 5, 6] showed that the SGLD [4] converges weakly to the true posterior. In [7], the author studied the sample-path inconsistency of the Hamiltonian PDE with stochastic gradients (but not the SGHMC), and pointed out its incompatibility with data subsampling. However, real applications only require convergence in the weak sense, i.e., instead of requiring sample-wise convergence as in [7], only laws of sample paths are of concern***For completeness, we provide mean sample-path properties of the SGHMC (similar to [7]) in Appendix J.. Very recently, the invariance measure of an SG-MCMC with a specific stochastic gradient noise was studied in [9]. However, the technique is not readily applicable to our general setting.

In this paper we focus on general SG-MCMCs, and study the role of their numerical integrators. Our main contributions include: i) From a theoretical viewpoint, we prove weak convergence results for general SG-MCMCs, which are of practical interest. Specifically, for a th-order numerical integrator, the bias of the expected sample average of an SG-MCMC at iteration is upper bounded by with optimal step size , and the MSE by with optimal . This generalizes the results of the SGLD with an Euler integrator () in [8, 5, 6], and is better when ; ii) From a practical perspective, we introduce a numerically efficient 2nd-order integrator, based on symmetric splitting schemes [9]. When applied to the SGHMC, it outperforms existing algorithms, including the SGLD and SGHMC with Euler integrators, considering both synthetic and large real datasets.

## 2 Preliminaries & Two Approximation Errors in SG-MCMCs

In weak convergence analysis, instead of working directly with sample-paths in (1), we study how the expected value of any suitably smooth statistic of evolves in time. This motivates the introduction of an (infinitesimal) generator. Formally, the generator of the diffusion (1) is defined for any compactly supported twice differentiable function , such that,

 Lf(Xt) ≜limh→0+E[f(Xt+h)]−f(Xt)h=(F(Xt)⋅∇+12(σ(Xt)σ(Xt)T):∇∇T)f(Xt) ,

where , , means approaches zero along the positive real axis. is associated with an integrated form via Kolmogorov’s backward equationMore details of the equation are provided in Appendix B. Specifically, under mild conditions on , we can expand the operator up to the th-order () such that the remainder terms are bounded by . Refer to [10] for more details. We will assume these conditions to hold for the ’s in this paper. : , where denotes the exact solution of the diffusion (1) at time . The operator is called the Kolmogorov operator for the diffusion (1). Since diffusion (1) is continuous, it is generally infeasible to solve analytically (so is ). In practice, a local numerical integrator is used for every small step , with the corresponding Kolmogorov operator approximating . Let denote the approximate sample path from such a numerical integrator; similarly, we have . Let denote the composition of two operators and , i.e., is evaluated on the output of . For time , we have the following approximation

with compositions, where is obtained by decomposing into sub-operators, each for a minibatch of data, while approximation is manifested by approximating the infeasible with from a feasible integrator, e.g., the symmetric splitting integrator proposed later, such that is close to the exact expectation . The latter is the first approximation error introduced in SG-MCMCs. Formally, to characterize the degree of approximation accuracies for different numerical methods, we use the following definition.

###### Definition 1.

An integrator is said to be a th-order local integrator if for any smooth and bounded function , the corresponding Kolmogorov operator satisfies the following relation:

 Phf(x)=ehLf(x)+O(hK+1) . (2)

The second approximation error is manifested when handling large data. Specifically, the SGLD and SGHMC use stochastic gradients in the 1st and 2nd-order LDs, respectively, by replacing in and the full negative log-posterior with a scaled log-posterior, , from the -th minibatch. We denote the corresponding generators with stochastic gradients as , e.g., the generator in the -th minibatch for the SGHMC becomes , where . As a result, in SG-MCMC algorithms, we use the noisy operator to approximate such that , where denotes the numerical sample-path with stochastic gradient noise, i.e.,

 E[f(XeT)]\mathclapB1≃eh~LL∘…∘eh~L1f(X0)\mathclapB2≃~PLh∘…∘~P1hf(X0)=E[f(Xn,sT)]. (3)

Approximations and in (3) are from the stochastic gradient and numerical integrator approximations, respectively. Similarly, we say corresponds to a th-order local integrator of if . In the following sections, we focus on SG-MCMCs which use numerical integrators with stochastic gradients, and for the first time analyze how the two introduced errors affect their convergence behaviors. For notational simplicity, we henceforth use to represent the approximate sample-path .

## 3 Convergence Analysis

This section develops theory to analyze finite-time convergence properties of general SG-MCMCs with both fixed and decreasing step sizes, as well as their asymptotic invariant measures.

### 3.1 Finite-time error analysis

Given an ergodicSee [11, 6] for conditions to ensure (1) is ergodic. Itô diffusion (1) with an invariant measure , the posterior average is defined as: for some test function of interest. For a given numerical method with generated samples , we use the sample average defined as to approximate . In the analysis, we define a functional that solves the following Poisson Equation:

 Lψ(Xlh)=ϕ(Xlh)−¯ϕ,or equivalently,1LL∑l=1Lψ(Xlh)=^ϕ−¯ϕ. (4)

The solution functional characterizes the difference between and the posterior average for every , thus would typically possess a unique solution, which is at least as smooth as under the elliptic or hypoelliptic settings [12]. In the unbounded domain of , to make the presentation simple, we follow [6] and make certain assumptions on the solution functional, , of the Poisson equation (4), which are used in the detailed proofs. Extensive empirical results have indicated the assumptions to hold in many real applications, though extra work is needed for theoretical verifications for different models, which is beyond the scope of this paper.

###### Assumption 1.

and its up to 3rd-order derivatives, , are bounded by a function§§§The existence of such function can be translated into finding a Lyapunov function for the corresponding SDEs, an important topic in PDE literatures [13]. See Assumption 4.1 in [6] and Appendix C for more details. , i.e., for , . Furthermore, the expectation of on is bounded: , and is smooth such that , for some .

We emphasize that our proof techniques are related to those of the SGLD [12, 6], but with significant distinctions in that, instead of expanding the function [6], whose parameter does not endow an explicit form in general SG-MCMCs, we start from expanding the Kolmogorov’s backward equation for each minibatch. Moreover, our techniques apply for general SG-MCMCs, instead of for one specific algorithm. More specifically, given a th-order local integrator with the corresponding Kolmogorov operator , according to Definition 1 and (3), the Kolmogorov’s backward equation for the -th minibatch can be expanded as:

 E[ψ(Xlh)] =~Plhψ(X(l−1)h)=eh~Llψ(X(l−1)h)+O(hK+1) =(I+h~Ll)ψ(X(l−1)h)+K∑k=2hkk!~Lklψ(X(l−1)h)+O(hK+1) , (5)

where is the identity map. Recall that , e.g., in SGHMC. By further using the Poisson equation (4) to simplify related terms associated with , after some algebra shown in Appendix D, the bias can be derived from (3.1) as:

 ∣∣ ∣∣E[ψ(Xlh)]−ψ(X0)Lh−1L∑lE[ΔVlψ(X(l−1)h)]−K∑k=2hk−1k!LL∑l=1E[~Lklψ(X(l−1)h)]∣∣ ∣∣+O(hK) .

All terms in the above equation can be bounded, with details provided in Appendix D. This gives us a bound for the bias of an SG-MCMC algorithm in Theorem 2.

###### Theorem 2.

Under Assumption 1, let be the operator norm. The bias of an SG-MCMC with a th-order integrator at time can be bounded as:

 ∣∣E^ϕ−¯ϕ∣∣=O(1Lh+∑l∥EΔVl∥L+hK) .

Note the bound above includes the term

, measuring the difference between the expectation of stochastic gradients and the true gradient. It vanishes when the stochastic gradient is an unbiased estimation of the exact gradient, an assumption made in the SGLD. This on the other hand indicates that if the stochastic gradient is biased,

might diverge when the growth of is faster than . We point this out to show our result to be more informative than that of the SGLD [6], though this case might not happen in real applications. By expanding the proof for the bias, we are also able to bound the MSE of SG-MCMC algorithms, given in Theorem 3.

###### Theorem 3.

Under Assumption 1, and assume is an unbiased estimate of . For a smooth test function , the MSE of an SG-MCMC with a th-order integrator at time is bounded, for some independent of , as

 E(^ϕ−¯ϕ)2≤C⎛⎝1L∑lE∥ΔVl∥2L+1Lh+h2K⎞⎠ .

Compared to the SGLD [6], the extra term

relates to the variance of noisy gradients. As long as the variance is bounded, the MSE still converges with the same rate. Specifically, when optimizing bounds for the bias and MSE, the optimal bias decreases at a rate of

with step size ; while this is with step size for the MSETo compare with the standard MCMC convergence rate of , the rate needs to be taken a square root.. These rates decrease faster than those of the SGLD [6] when . The case of for the SGHMC with our proposed symmetric splitting integrator is discussed in Section 4.

### 3.2 Stationary invariant measures

The asymptotic invariant measures of SG-MCMCs correspond to approaching infinity in the above analysis. According to the bias and MSE above, asymptotically () the sample average

is a random variable with mean

, and variance , close to the true . This section defines distance between measures, and studies more formally how the approximation errors affect the invariant measures of SG-MCMC algorithms.

First we note that under mild conditions, the existence of a stationary invariant measure for an SG-MCMC can be guaranteed by application of the Krylov–Bogolyubov Theorem [14]. Examining the conditions is beyond the scope of this paper. For simplicity, we follow [12] and assume stationary invariant measures do exist for SG-MCMCs. We denote the corresponding invariant measure as , and the true posterior of a model as . Similar to [12], we assume our numerical solver is geometric ergodic, meaning that for a test function , we have for any from the ergodic theorem, where denotes the expectation conditional on . The geometric ergodicity implies that the integration is independent of the starting point of an algorithm. Given this, we have the following theorem on invariant measures of SG-MCMCs.

###### Theorem 4.

Assume that a th-order integrator is geometric ergodic and its invariance measures exist. Define the distance between the invariant measures and as: . Then any invariant measure of an SG-MCMC is close to with an error up to an order of , i.e., there exists some such that: .

For a th-order integrator with full gradients, the corresponding invariant measure has been shown to be bounded by an order of [12, 9]. As a result, Theorem 4 suggests only orders of numerical approximations but not the stochastic gradient approximation affect the asymptotic invariant measure of an SG-MCMC algorithm. This is also reflected by experiments presented in Section 5.

### 3.3 SG-MCMCs with decreasing step sizes

The original SGLD was first proposed with a decreasing-step-size sequence [4], instead of fixing step sizes, as analyzed in [6]. In [5], the authors provide theoretical foundations on its asymptotic convergence properties. We demonstrate in this section that for general SG-MCMC algorithms, decreasing step sizes for each minibatch are also feasible. Note our techniques here are different from those used for the decreasing-step-size SGLD [5], which interestingly result in similar convergence patterns. Specifically, by adapting the same techniques used in the previous sections, we establish conditions on the step size sequence to ensure asymptotic convergence, and develop theory on their finite-time ergodic error as well. To guarantee asymptotic consistency, the following conditions on decreasing step size sequences are required.

###### Assumption 2.

The step sizes are decreasingActually the sequence need not be decreasing; we assume it is decreasing for simplicity., i.e., , and satisfy that 1) ; and 2) .

Denote the finite sum of step sizes as . Under Assumption 2, we need to modify the sample average defined in Section 3.1 as a weighted summation of : . For simplicity, we assume to be an unbiased estimate of such that . Extending techniques in previous sections, we develop the following bounds for the bias and MSE.

###### Theorem 5.

Under Assumptions 1 and 2, for a smooth test function , the bias and MSE of a decreasing-step-size SG-MCMC with a th-order integrator at time are bounded as:

 BIAS: ∣∣E~ϕ−¯ϕ∣∣=O(1SL+∑Ll=1hK+1lSL) (6) MSE: E(~ϕ−¯ϕ)2≤C(∑lh2lS2LE∥ΔVl∥2+1SL+(∑Ll=1hK+1l)2S2L) . (7)

As a result, the asymptotic bias approaches 0 according to the assumptions. If further assuming******The assumption of satisfies this requirement, but is weaker than the original assumption. , the MSE also goes to 0. In words, the decreasing-step-size SG-MCMCs are consistent.

Among the kinds of decreasing step size sequences, a commonly recognized one is for . We show in the following corollary that such a sequence leads to a valid sequence.

###### Corollary 6.

Using the step size sequences for , all the step size assumptions in Theorem 5 are satisfied. As a result, the bias and MSE approach zero asymptotically, i.e., the sample average is asymptotically consistent with the posterior average .

###### Remark 7.

Theorem 5 indicates the sample average asymptotically converges to the true posterior average . It is possible to find out the optimal decreasing rates for the specific decreasing sequence . Specifically, using the bounds for (see the proof of Corollary 6), for the two terms in the bias (6) in Theorem 5, decreases at a rate of , whereas decreases as . The balance between these two terms is achieved when , which agrees with Theorem 2 on the optimal rate of fixed-step-size SG-MCMCs. Similarly, for the MSE (7), the first term decreases as , independent of , while the second and third terms decrease as and , respectively, thus the balance is achieved when , which also agrees with the optimal rate for the fixed-step-size MSE in Theorem 3.

According to Theorem 5, one theoretical advantage of decreasing-step-size SG-MCMCs over fixed-step-size variants is the asymptotically unbiased estimation of posterior averages, though the benefit might not be significant in large-scale real applications where the asymptotic regime is not reached.

## 4 Practical Numerical Integrators

Given the theory for SG-MCMCs with high-order integrators, we here propose a 2nd-order symmetric splitting integrator for practical use. The Euler integrator is known as a 1st-order integrator; the proof and its detailed applications on the SGLD and SGHMC are given in Appendix I.

The main idea of the symmetric splitting scheme is to split the local generator into several sub-generators that can be solved analytically††††††This is different from the traditional splitting in SDE literatures[15, 9], where instead of is split.. Unfortunately, one cannot easily apply a splitting scheme with the SGLD. However, for the SGHMC, it can be readily split into: , where

 LA=p⋅∇\mathchar28946\relax,LB=−Dp⋅∇p,LOl=−∇\mathchar28946\relax~U(\mathchar28946\relax)⋅∇p+2DIn:∇p∇Tp . (8)

These sub-generators correspond to the following SDEs, which are all analytically solvable:

 (9)

Based on these sub-SDEs, the local Kolmogorov operator is defined as:

 E[f(Xlh)]=~Plhf(X(l−1)h),where, ~Plh≜eh2LA∘eh2LB∘ehLOl∘eh2LB∘eh2LA,

so that the corresponding updates for consist of the following 5 steps:

 \mathchar28946\relax(1)lh=\mathchar28946\relax(l−1)h+p(l−1)hh/2⇒p(1)lh=e−Dh/2p(l−1)h⇒p(2)lh=p(1)lh−∇\mathchar28946\relax~Ul(\mathchar28946\relax(1)lh)h+√2Dhζl

where are intermediate variables. We denote such a splitting method as the ABOBA scheme. From the Markovian property of a Kolmogorov operator, it is readily seen that all such symmetric splitting schemes (with different orders of ‘A’, ‘B’ and ‘O’) are equivalent [15]. Lemma 8 below shows the symmetric splitting scheme is a 2nd-order local integrator.

###### Lemma 8.

The symmetric splitting scheme is a 2nd-order local integrator, i.e., the corresponding Kolmogorov operator satisfies: .

When this integrator is applied to the SGHMC, the following properties can be obtained.

###### Remark 9.

Applying Theorem 2 to the SGHMC with the symmetric splitting scheme (), the bias is bounded as: . The optimal bias decreasing rate is , compared to for the SGLD [6]. Similarly, the MSE is bounded by: , decreasing optimally as with step size , compared to the MSE of for the SGLD [6]. This indicates that the SGHMC with the splitting integrator converges faster than the SGLD and SGHMC with 1st-order Euler integrators.

###### Remark 10.

For a decreasing-step-size SGHMC, based on Remark 7, the optimal step size decreasing rate for the bias is , and for the MSE. These agree with their fixed-step-size counterparts in Remark 9, thus are faster than the SGLD/SGHMC with 1st-order Euler integrators.

## 5 Experiments

We here verify our theory and compare with related algorithms on both synthetic data and large-scale machine learning applications.

#### Synthetic data

We consider a standard Gaussian model where . 1000 data samples are generated, and every minibatch in the stochastic gradient is of size 10. The test function is defined as , with explicit expression for the posterior average. To evaluate the expectations in the bias and MSE, we average over 200 runs with random initializations.

First we compare the invariant measures (with ) of the proposed splitting integrator and Euler integrator for the SGHMC. Results of the SGLD are omitted since they are not as competitive. Figure 1 plots the biases with different step sizes. It is clear that the Euler integrator has larger biases in the invariant measure, and quickly explodes when the step size becomes large, which does not happen for the splitting integrator. In real applications we also find this happen frequently (shown in the next section), making the Euler scheme an unstable integrator.

Next we examine the asymptotically optimal step size rates for the SGHMC. From the theory we know these are for the bias and for the MSE, in both fixed-step-size SGHMC (SGHMC-F) and decreasing-step-size SGHMC (SGHMC-D). For the step sizes, we did a grid search to select the best prefactors, which resulted in for the SGHMC-F and for the SGHMC-D, with different values. We plot the traces of the bias for the SGHMC-D and the MSE for the SGHMC-F in Figure 2. Similar results for the bias of the SGHMC-F and the MSE of the SGHMC-D are plotted in Appendix K. We find that when rates are smaller than the theoretically optimal rates, i.e., (bias) and (MSE), the bias and MSE tend to decrease faster than the optimal rates at the beginning (especially for the SGHMC-F), but eventually they slow down and are surpassed by the optimal rates, consistent with the asymptotic theory. This also suggests that if only a small number of iterations were feasible, setting a larger step size than the theoretically optimal one might be beneficial in practice.

Finally, we study the relative convergence speed of the SGHMC and SGLD. We test both fixed-step-size and decreasing-step-size versions. For fixed-step-size experiments, the step sizes are set to , with chosen according to the theory for SGLD and SGHMC. To provide a fair comparison, the constants are selected via a grid search from to 0.5 with an interval of for , it is then fixed in the other runs with different values. The parameter in the SGHMC is selected within as well. For decreasing-step-size experiments, an initial step size is chosen within with an interval of for different algorithms‡‡‡‡‡‡Using the same initial step size is not fair because the SGLD requires much smaller step sizes., and then it decreases according to their theoretical optimal rates. Figure 3 shows a comparison of the biases for the SGHMC and SGLD. As indicated by both theory and experiments, the SGHMC with the splitting integrator yields a faster convergence speed than the SGLD with an Euler integrator.

#### Large-scale machine learning applications

For real applications, we test the SGLD with an Euler integrator, the SGHMC with the splitting integrator (SGHMC-S), and the SGHMC with an Euler integrator (SGHMC-E). First we test them on the latent Dirichlet allocation model (LDA) [16]. The data used consists of 10M randomly downloaded documents from Wikipedia, using scripts provided in [17]. We randomly select 1K documents for testing and validation, respectively. As in [17, 18], the vocabulary size is 7,702. We use the Expanded-Natural reparametrization trick to sample from the probabilistic simplex [19]. The step sizes are chosen from , and parameter from . The minibatch size is set to 100, with one pass of the whole data in the experiments (and therefore ). We collect 300 posterior samples to calculate test perplexities, with a standard holdout technique as described in [18].

Next a recently studied sigmoid belief network model (SBN) [20]

is tested, which is a directed counterpart of the popular RBM model. We use a one layer model where the bottom layer corresponds to binary observed data, which is generated from the hidden layer (also binary) via a sigmoid function. As shown in

[18], the SBN is readily learned by SG-MCMCs. We test the model on the MNIST dataset, which consists of 60K hand written digits of size for training, and 10K for testing. Again the step sizes are chosen from , from . The minibatch is set to 200, with 5000 iterations for training. Like applied for the RBM [21], an advance technique called anneal importance sampler (AIS) is adopted for calculating test likelihoods.

We briefly describe the results here, more details are provided in Appendix K. For LDA with 200 topics, the best test perplexities for the SGHMC-S, SGHMC-E and SGLD are 1168, 1180 and 2496, respectively; while these are 1157, 1187 and 2511, respectively, for 500 topics. Similar to the synthetic experiments, we also observed SGHMC-E crashed when using large step sizes. This is illustrated more clearly in Figure 4. For the SBN with 100 hidden units, we obtain negative test log-likelihoods of 103, 105 and 126 for the SGHMC-S, SGHMC-E and SGLD, respectively; and these are 98, 100, and 110 for 200 hidden units. Note the SGHMC-S on SBN yields state-of-the-art results on test likelihoods compared to [22], which was 113 for 200 hidden units. A decrease of 2 units in the neg-log-likelihood with AIS is considered to be a reasonable gain [20], which is approximately equal to the gain from a shallow to a deep model [22]. SGHMC-S is more accuracy and robust than SGHMC-E due to its 2nd-order splitting integrator.

## 6 Conclusion

For the first time, we develop theory to analyze finite-time ergodic errors, as well as asymptotic invariant measures, of general SG-MCMCs with high-order integrators. Our theory applies for both fixed and decreasing step size SG-MCMCs, which are shown to be equivalent in terms of convergence rates, and are faster with our proposed 2nd-order integrator than previous SG-MCMCs with 1st-order Euler integrators. Experiments on both synthetic and large real datasets validate our theory. The theory also indicates that with increasing order of numerical integrators, the convergence rate of an SG-MCMC is able to theoretically approach the standard MCMC convergence rate. Given the theoretical convergence results, SG-MCMCs can be used effectively in real applications.

#### Acknowledgments

Supported in part by ARO, DARPA, DOE, NGA and ONR. We acknowledge Jonathan C. Mattingly and Chunyuan Li for inspiring discussions; David Carlson for the AIS codes.

## References

• Chen et al. [2014] T. Chen, E. B. Fox, and C. Guestrin. Stochastic gradient Hamiltonian Monte Carlo. In ICML, 2014.
• Ding et al. [2014] N. Ding, Y. Fang, R. Babbush, C. Chen, R. D. Skeel, and H. Neven. Bayesian sampling using stochastic gradient thermostats. In NIPS, 2014.
• Risken [1989] H. Risken. The Fokker-Planck equation. Springer-Verlag, New York, 1989.
• Welling and Teh [2011] M. Welling and Y. W. Teh. Bayesian learning via stochastic gradient Langevin dynamics. In ICML, 2011.
• Teh et al. [2014] Y. W. Teh, A. H. Thiery, and S. J. Vollmer. Consistency and fluctuations for stochastic gradient Langevin dynamics. Technical Report arXiv:1409.0578, University of Oxford, UK, Sep. 2014. URL http://arxiv.org/abs/1409.0578.
• Vollmer et al. [2015] S. J. Vollmer, K. C. Zygalakis, and Y. W. Teh. (Non-)asymptotic properties of stochastic gradient Langevin dynamics. Technical Report arXiv:1501.00438, University of Oxford, UK, January 2015. URL http://arxiv.org/abs/1501.00438.
• Betancourt [2015] M. Betancourt. The fundamental incompatibility of Hamiltonian Monte Carlo and data subsampling. In ICML, 2015.
• Sato and Nakagawa [2014] I. Sato and H. Nakagawa. Approximation analysis of stochastic gradient Langevin dynamics by using Fokker-Planck equation and Itô process. In ICML, 2014.
• Leimkuhler and Shang [2015] B. Leimkuhler and X. Shang. Adaptive thermostats for noisy gradient systems. Technical Report arXiv:1505.06889v1, University of Edinburgh, UK, May 2015. URL http://arxiv.org/abs/1505.06889.
• Abdulle et al. [2015] A. Abdulle, G. Vilmart, and K. C. Zygalakis. Long time accuracy of Lie–Trotter splitting methods for Langevin dynamics. SIAM J. NUMER. ANAL., 53(1):1–16, 2015.
• Hasminskii [2012] R. Hasminskii. Stochastic stability of differential equations. Springer-Verlag Berlin Heidelberg, 2012.
• Mattingly et al. [2010] J. C. Mattingly, A. M. Stuart, and M. V. Tretyakov. Construction of numerical time-average and stationary measures via Poisson equations. SIAM J. NUMER. ANAL., 48(2):552–577, 2010.
• Giesl [2007] P. Giesl.

Construction of global Lyapunov functions using radial basis functions

.
Springer Berlin Heidelberg, 2007.
• Bogoliubov and Krylov [1937] N. N. Bogoliubov and N. M. Krylov. La theorie generalie de la mesure dans son application a l’etude de systemes dynamiques de la mecanique non-lineaire. Ann. Math. II (in French), 38(1):65–113, 1937.
• Leimkuhler and Matthews [2013] B. Leimkuhler and C. Matthews. Rational construction of stochastic numerical methods for molecular sampling. AMRX, 2013(1):34–56, 2013.
• Blei et al. [2003] D. M. Blei, A. Y. Ng, and M. I. Jordan. Latent Dirichlet allocation. JMLR, 2003.
• Hoffman et al. [2010] M. D. Hoffman, D. M. Blei, and F. Bach. Online learning for latent Dirichlet allocation. In NIPS, 2010.
• Gan et al. [2015a] Z. Gan, C. Chen, R. Henao, D. Carlson, and L. Carin. Scalable deep Poisson factor analysis for topic modeling. In ICML, 2015a.
• Patterson and Teh [2013] S. Patterson and Y. W. Teh.

Stochastic gradient Riemannian Langevin dynamics on the probability simplex.

In NIPS, 2013.
• Gan et al. [2015b] Z. Gan, R. Henao, D. Carlson, and L. Carin. Learning deep sigmoid belief networks with data augmentation. In AISTATS, 2015b.
• Salakhutdinov and Murray [2008] R. Salakhutdinov and I. Murray.

On the quantitative analysis of deep belief networks.

In ICML, 2008.
• Mnih and Gregor [2014] A. Mnih and K. Gregor. Neural variational inference and learning in belief networks. In ICML, 2014.
• Debussche and Faou [2012] A. Debussche and E. Faou. Weak backward error analysis for SDEs. SIAM J. NUMER. ANAL., 50(3):1734–1752, 2012.
• Kopec [2014] M. Kopec. Weak backward error analysis for overdamped Langevin processes. IMA J. NUMER. ANAL., 2014.

## Appendix A Representative Stochastic Gradient MCMC Algorithms

This section briefly introduces three recently proposed stochastic gradient MCMC algorithms, including the stochastic gradient Langevin dynamic (SGLD) [4], the stochastic gradient Hamiltonian MCMC (SGHMC) [1], and the stochastic gradient Nosé-Hoover thermostat [2] (SGNHT).

Given data , a generative model with model parameter , and prior , we want to compute the posterior:

 π(\mathchar28946\relax)≜p(\mathchar28946\relax|X)∝p(X|\mathchar28946\relax)p(\mathchar28946\relax)≜e−U(\mathchar28946\relax) .

### a.1 Stochastic gradient Langevin dynamics

The SGLD [4] is based on the following 1st-order Langevin dynamic defined as:

 d\mathchar28946\relax=−12∇\mathchar28946\relaxU(\mathchar28946\relax)dt+dW , (10)

where is the standard Brownian motion. We can show via the Fokker–Planck equation that the equilibrium distribution of (10) is:

 p(\mathchar28946\relax)=π(\mathchar28946\relax) .

As described in the main text, when sampling from this continuous-time diffusion, two approximations are adopted, e.g., a numerical integrator and a stochastic gradient version of the log-likelihood from the -th minibatch. This results in the following SGLD algorithm.

### a.2 Stochastic gradient Hamiltonian MCMCs

The SGHMC [1] is based on the 2nd-order Langevin dynamic defined as:

 {d\mathchar28946\relax=pdtdp=−∇\mathchar28946\relaxU(\mathchar28946\relax)dt−Dpdt+√2DdW , (11)

where is a constant independent of and . Again we can show that the equilibrium distribution of (11) is:

 P(\mathchar28946\relax,p)∝e−U(\mathchar28946\relax)+pTp2 .

Similar to the SGLD, we use the Euler scheme to simulate the dynamic (11), shown in Algorithm 2.

### a.3 Stochastic gradient Nośe-Hoover thermostats

The SGNHT [2] is based on the Nośe-Hoover thermostat defined as:

 ⎧⎪ ⎪⎨⎪ ⎪⎩d\mathchar28946\relax=pdtdp=−∇\mathchar28946\relaxU(\mathchar28946\relax)dt−ξpdt+√2DdWdξ=(pTp/n−1) , (12)

If is independent of and , it can also be shown that the equilibrium distribution of (12) is [2]:

 P(\mathchar28946\relax,p,ξ)∝e−U(\mathchar28946\relax)−12pTp+12(ξ−D)2 .

The SGNHT is much more interesting than the SGHMC when considering subsampling data in each iteration, as the covariance in SGHMC is hard to estimate, a thermostat is used to adaptively control the system temperature, thus automatically estimate the unknown . The whole algorithm is shown in Algorithm 3.

## Appendix B More Details on Kolmogorov’s Backward Equation

The generator is used in the formulation of Kolmogorov’s backward equation, which intuitively tells us how the expected value of any suitably smooth statistic of evolves in time. More precisely:

###### Definition 11 (Kolmogorov’s Backward Equation).

Let , then

satisfies the following partial differential equation, known as

Kolmogorov’s backward equation:

 {∂u∂t(t,x)=Lu(t,x) ,t>0,x∈Rnu(0,x)=ϕ(x),x∈Rn (13)

Based on the definition, we can write so that is the transition semigroup associated with the Markov process [23] (also called the Kolmogorov operator). Note that the Kolmogorov’s backward equation can be written in another form as:

 u(t,x)=E[ϕ(Xt)]=etLϕ(x) , (14)

where is the exponential map operator associated with the generator defined as:

 etL≜I+∞∑i=1(tL)ii! ,

with being the identity map. This is obtained by expanding in time by using Taylor expansion [23]:

 u(t,x) =u(0,x)+∞∑i=1tii!didtiu(t,x)|t=0 =u(0,x)+∞∑i=1tii!di−1dti−1ddtu(t,x)|t=0 =u(0,x)+∞∑i=1tii!Ldi−1dti−1u(t,x)|t=0 =ϕ(x)+∞∑i=1tii!Liϕ(x)=etLϕ(x) . (15)

The form (14) instead of the original form (13) of the Kolmogorov’s backward equation is used in our analysis. To be able to expand the form (14) to some particular order such that remainder terms are bounded, the following assumption is required [24].

###### Assumption 3.

Assume 1) is with bounded derivatives of any order, furthermore, and 2) for some positive integer . Under these assumptions, series of the generator expansion can be bounded, thus (B) can be written in the following form [24, 6]:

 u(t,x)=ϕ(x)+ℓ∑i=1tii!Liϕ(x)+tℓ+1rℓ(F,ϕ)(x) , (16)

with for some constant .

## Appendix C More Comments on Assumption 1

Assumption 1 assumes that the solution functional of the Poisson equation (4) satisfies: and its up to 3-rd order derivatives, , are bounded by a function , i.e., for , . Furthermore, is smooth such that , for some . Finally, for . This is summarized as:

 suplEVp(Xlh) <∞ (17) sups∈(0,1)Vp(sX+(1−s)Y) ≤C(Vp(X)+Vp(Y)) (18) ∥Dkψ∥ ≤CkVpk (19)

Compared to the SGLD case [6], in our proofs, we only need be up to 3 in (19) instead of 4. More specifically, the proof for the bias only needs be up to 0 given other assumptions in this paper, and the proof for the MSE needs be up to 3.

As long as the corresponding SDE is hypoelliptic, meaning that the Brownian motion is able to propagate to the other variables of the dynamics [12], e.g., the model parameter in SGHMC, we can extend Assumption 4.1 of [6] to our setting. Thus we have that (17) is equivalent to finding a function ( is the dimension of , e.g., including the momentum in SGHMC), which tends to infinity as , and is twice differentiable with bounded second derivatives and satisfies the following conditions:

• is a Lyapunov function of the SDE, i.e., there exists constants , such that for , we have .

• There exists an exponent such that