Stochastic Gradient Hamiltonian Monte Carlo with Variance
Reduction for Bayesian Inference

arXiv:1803.11159v1 [cs.LG] 29 Mar 2018

Zhize Li†
Tianyi Zhang†
Jian Li
IIIS, Tsinghua University
{zz-li14@mails, tianyi-z16@mails, lijian83@mail}.tsinghua.edu.cn
Abstract
Gradient-based Monte Carlo sampling algorithms, like Langevin dynamics and Hamiltonian Monte Carlo, are important methods for Bayesian inference. In large-scale settings, full-gradients are not affordable and thus stochastic
gradients evaluated on mini-batches are used as a replacement. In order to reduce the high variance of noisy stochastic gradients, [Dubey et al., 2016] applied the standard variance reduction technique on stochastic gradient Langevin
dynamics and obtained both theoretical and experimental improvements. In this paper, we apply the variance reduction tricks on Hamiltonian Monte Carlo and achieve better theoretical convergence results compared with the
variance-reduced Langevin dynamics. Moreover, we apply the symmetric splitting scheme in our variance-reduced
Hamiltonian Monte Carlo algorithms to further improve the theoretical results. The experimental results are also
consistent with the theoretical results. As our experiment shows, variance-reduced Hamiltonian Monte Carlo demonstrates better performance than variance-reduced Langevin dynamics in Bayesian regression and classification tasks
on real-world datasets.

1

Introduction

Gradient-based Monte Carlo algorithms are useful tools for sampling posterior distributions. Similar to gradient
descent algorithms, gradient-based Monte Carlo generates posterior samples iteratively using the gradient of loglikelihood.
Langevin dynamics (LD) and Hamiltonian Monte Carlo (HMC) [Duane et al., 1987, Neal et al., 2011] are two
important examples of gradient-based Monte Carlo sampling algorithms that are widely used in Bayesian inference.
Since calculating likelihood on large datasets is expensive, people use stochastic gradients [Robbins and Monro,
1951] in place of full gradient, and have, for both Langevin dynamics and Hamiltonian Monte Carlo, developed their
stochastic gradient counterparts[Welling and Teh, 2011, Chen et al., 2014]. Stochastic gradient Hamiltonian Monte
Carlo (SGHMC) usually converges faster than stochastic gradient Langevin dynamics (SGLD) in practical machine
learning tasks like covariance estimation of bivariate Gaussian and Bayesian neural networks for classification on
MNIST dataset, as demonstrated in [Welling and Teh, 2011]. Similar phenomenon was also observed in [Chen et al.,
2015] where SGHMC and SGLD were compared on both synthetic and real-world datasets. Intuitively speaking,
comparing against SGLD, SGHMC has a momentum term that enables it to explore the parameter space of posterior
distribution much faster when the gradient of log-likelihood becomes smaller.
Very recently, [Dubey et al., 2016] borrowed the standard variance reduction techniques from the stochastic optimization literature [Johnson and Zhang, 2013, Defazio et al., 2014] and applied them on SGLD to obtain two variancereduced SGLD algorithms (called SAGA-LD and SVRG-LD) with improved theoretical results and practical performance. Because of the superiority of SGHMC over SGLD in terms of convergence rate in a wide range of machine
learning tasks, it would be a natural question whether such variance reduction techniques can be applied on SGHMC
to achieve better results than variance-reduced SGLD.
† denotes

equal contribution

1

1.1

Our contribution

1. We propose two variance-reduced versions of Hamiltonian Monte Carlo algorithms (called SVRG-HMC and
SAGA-HMC) using the standard approaches from [Johnson and Zhang, 2013, Defazio et al., 2014]. Compared
with SVRG/SAGA-LD [Dubey et al., 2016], our algorithms guarantee improved theoretical convergence results
due to the extra momentum term in HMC (see Corollary 3).
2. Moreover, we combine the proposed SVRG/SAGA-HMC algorithms with the symmetric splitting scheme [Chen
et al., 2015, Leimkuhler and Shang, 2016] to extend them to 2nd-order integrators, which further improve the
dependency on step size (see Theorem 2 and 5). We denote these two algorithms as SVRG2nd-HMC and
SAGA2nd-HMC.
3. Finally, we evaluate our algorithms on real-world datasets and compare them with SVRG-LD [Dubey et al.,
2016]; as it turns out, our algorithms converge markedly faster than the benchmarks.

1.2

Related work

Langevin dynamics and Hamiltonian Monte Carlo are two important sampling algorithms that are widely used in
Bayesian inference. Many literatures studied how to develop the variants of them to achieve improved performance,
especially for scalability for large datasets. [Welling and Teh, 2011] started this direction with the notable work
stochastic gradient Langevin dynamics (SGLD). [Ahn et al., 2012] proposed a modification to SGLD reminiscent
of Fisher scoring to better estimate the gradient noise variance, with lower classification error rates on HHP dataset
and MNIST dataset. [Chen et al., 2014] developed the stochastic gradient version of HMC (SGHMC), with a quite
nontrivial approach (momentum) different from SGLD. ding2014bayesian further improved SGHMC by a new dynamics to better control the gradient noise, and the proposed stochastic gradient Nosé-Hoover thermostats (SGNHT)
outperforms SGHMC on MNIST dataset.
Variant settings of Markov Chain Monte Carlo (MCMC) are also considered. [Girolami and Calderhead, 2011]
enhanced LD and HMC by exploring the Riemannian structure of the target distribution, with Riemannian manifold
LD and HMC (RMLD and RMHMC, respectively). [Byrne and Girolami, 2013] developed geodesic Monte Carlo
(GMC) that is applicable to Riemannian manifolds with no global coordinate systems. Large-scale variants of RMLD,
RMHMC and GMC with stochastic gradient were developed by [Patterson and Teh, 2013], [Ma et al., 2015] and
[Liu et al., 2016], respectively. [Ahn et al., 2014] studied the behaviour of stochastic gradient MCMC algorithms for
distributed posterior inference.

2

Preliminary

Qn
Let X = {xi }ni=1 be a d-dimensional dataset that follows the distribution
Pr(X|θ) = i=1 Pr(xi |θ). Then, we are
Qn
interested in sampling the posterior distribution Pr(θ|X) ∝ Pr(θ)
Pn i=1 Pr(xi |θ) based on Hamiltonian Monte Carlo
algorithms. Let [n] denote the set {1, 2, . . . , n}. Define f (θ) = i=1 fi (θ)−log Pr(θ), where fi (θ) = − log Pr(xi |θ)
and i ∈ [n]. Similar to [Dubey et al., 2016], we assume that each fi is L-smooth and G-Lipschitz, for all i ∈ [n].
The general algorithmic framework maintains two sequences for t = 0, 1, . . . , T − 1 by the following discrete time
procedure:
√
˜ t + 2Dh · ξt
(1)
pt+1 = (1 − Dh)pt − h∇
θt+1 = θt + hpt+1

(2)

and then return the samples {θ1 , θ2 , . . . , θT } as an approximation to the stationary distribution Pr(θ|X). θt is the
parameter we wish to sample and pt is an auxiliary variable conventionally called the “momentum”. Here h is step
˜
size, D is a constant
Pcalled “friction”, ξt ∼ N(0, Id ) and ∇t is a mini-batch approximation of the full gradient ∇f (θt ).
n
˜
If we set ∇t = b i∈I ∇fi (θt ), I being a b-element index set uniformly randomly drawn (with replacement) from
{1, 2, . . . , n} as introduced in [Robbins and Monro, 1951], then the algorithm becomes SGHMC.

2

The above discrete time procedure provides an approximation to the continuous Hamiltonian Monte Carlo diffusion
process (θ, p):
dθ = pdt
dp = −∇θ f (θ)dt − Dpdt +

(3)

√
2DdW

(4)

Here W is a Wiener process. According to [Chen et al., 2015], the stationary joint distribution of (θ, p) is π(θ, p) ∝
p| p
e−f (θ)− 2 .
How do we evaluate the quality of the samples {θ1 , θ2 , . . . , θT }? Assuming φ : Rd → R is a smooth test function,
PT
we wish to upper bound the Mean-Squared Error (MSE) E(φ̂− φ̄)2 , where φ̂ = T1 t=1 φ(θt ) is the empirical average,
˜ t to
and φ̄ = Eθ∼Pr(θ|X) φ(θ) is the population average. So, the objective of our algorithm is to carefully design ∇
2
˜
minimize E(φ̂ − φ̄) in a faster way, where ∇t is a stochastic approximation of ∇f (θt ).
˜ t influences the value of E(φ̂ − φ̄)2 , define ψ(θ, p) to be the solution to the Poisson
To study how the choice of ∇
equation Lψ = φ(θ) − φ̄, L being the generator of Hamiltonian Monte Carlo diffusion process. In order to analyze the
theoretical convergence results related to the MSE E(φ̂ − φ̄)2 , we inherit the following assumption from [Chen et al.,
2015].
Assumption 1 ([Chen et al., 2015]) Function ψ is bounded up to 3rd-order derivatives by some real-valued function
Γ(θ, p), i.e. kDk ψk ≤ Ck Γqk where Dk is the kth order derivative for k = 0, 1, 2, 3, and Ck , qk > 0. Furthermore, the expectation of Γ on {(θt , pt )} is bounded, i.e. supt E[Γq (θt , pt )] < ∞ and that Γ is smooth such that
sups∈(0,1) Γq (sθ + (1 − s)θ0 , sp + (1 − s)p0 ) ≤ C(Γq (θ, p) + Γq (θ0 , p0 )), ∀θ, p, θ0 , p0 , q ≤ max 2qk for some constant
C > 0.
˜ t − ∇f (θt )) · ∇ for all t = 0, 1, 2, . . . , T − 1. When the above assumption holds, we
Define operator ∆Vt = (∇
have the following theorem by [Chen et al., 2015].
For the rest of this paper, for any two values A, B > 0, we say A . B if A = O(B), where the notation O(·) only
hides a constant factor independent of algorithm parameters T, n, D, h, G, b.
˜ t be an unbiased estimate of ∇f (θt ) for all t. Then under assumption 1, for a
Theorem 1 ([Chen et al., 2015]) Let ∇
smooth test function φ, the MSE of SGHMC is bounded in the following way:
PT −1
1
2
1
t=0 E(∆Vt ψ(θt , pt ))
2
T
+
+ h2
(5)
E(φ̂ − φ̄) .
T
Th
Similar to the [A2] assumption in [Dubey et al., 2016], we also need to make the following assumption which
˜ t − ∇f (θ)k2 .
relates ∆Vt φ(θ, p) to the difference k∇
˜ t − ∇f (θt )k2 for all 0 ≤ t < T .
Assumption 2 (∆Vt ψ(θt , pt ))2 . k∇
Combined with Theorem 1, Assumption 2 immediately yields the following corollary.
Corollary 1 Under Assumptions 1 and 2, we have:
PT −1
1
2
˜
1
t=0 Ek∇t − ∇f (θt )k
2
T
E(φ̂ − φ̄) .
+
+ h2
T
Th
˜ t to be the Robbins & Monro approximation of ∇f (θt ) [Robbins and Monro,
As we mentioned before, if we take ∇
1951], then it becomes SGHMC and the following corollary holds since all fi ’s are G-Lipschitz and E(X − EX)2 ≤
EX 2 for any random variable X.
Corollary 2 Under Assumptions 1 and 2, the MSE of SGHMC is bounded as:
PT −1
1
2
˜
1
t=0 Ek∇t − ∇f (θt )k
2
T
E(φ̂ − φ̄) .
+
+ h2
T
Th
n2 G2
1
≤
+
+ h2
bT
Th
3

(6)

3

Variance Reduction for Hamiltonian Monte Carlo

In this section, we introduce two versions of variance-reduced Hamiltonian Monte Carlo based on SVRG [Johnson
and Zhang, 2013] and SAGA [Defazio et al., 2014] respectively.

3.1

SVRG-HMC

In this subsection, we propose the SVRG-HMC algorithm (see Algorithm 1) which is
based on the SVRG algorithm.

˜ tK+k = −∇ log Pr(θtK+k )+ n P
As can be seenP
from Line 8 of Algorithm 1, we use ∇
i∈I ∇fi (θtK+k )−∇fi (w) +
b
n
g , where g is i=1 ∇fi (θtK ), as the stochastic estimation for the full gradient ∇f (θtK+k ).
Note that we initialize θ0 , p0 to be zero vectors in the algorithm only to simplify the theoretical analysis. It would
still work with an arbitrary initialization.
Algorithm 1: SVRG-HMC

10

parameters T, K, b, h > 0, Dh < 1, D ≥ 1;
initialize θ0 = p0 = 0;
for t = 0, 1, . . . , TP
/K − 1 do
n
compute g = i=1 ∇fi (θtK );
w = θtK ;
for k = 0, 1, . . . , K − 1 do
uniformly sample an index subset I ⊆ [n], |I| = b;


˜ tK+k = −∇ log Pr(θtK+k ) + n P
∇f
(θ
)
−
∇f
(w)
+ g;
∇
i
tK+k
i
i∈I
b
√
˜
ptK+k+1 = (1 − Dh)ptK+k − h∇tK+k + 2DhξtK+k ;
θtK+k+1 = θtK+k + hptK+k+1 ;

11

return {θt }1≤t≤T ;

1
2
3
4
5
6
7
8
9

The following theorem shows the convergence result for MSE of SVRG-HMC (Algorithm 1). We defer all the
proofs to Appendix A.
Theorem 2 Under Assumptions 1 and 2, the MSE of SVRG-HMC is bounded as:
n n2 G2 L2 n2 K 2 h2  √n2 G2 + D2 d 2 o
1
2
E[(φ̂ − φ̄) ] . min
,
+
+ h2
bT
bT
D − L2 n2 K 2 h3 b−1
Th

(7)

To see how the SVRG-HMC (Algorithm 1) is compared with SVRG-LD [Dubey et al., 2016], we restate their
results as following.
Theorem 3 ([Dubey et al., 2016]) Under Assumptions 1 and 2, the MSE of SVRG-LD is bounded as:
E[(φ̂ − φ̄)2 ] .

1
min{n2 G2 , n2 K 2 (n2 L2 h2 G2 + hd)}
+
+ h2
bT
Th

(8)

We assume n2 G2 > n2 K 2 (n2 L2 h2 G2 + hd). Otherwise the MSE upper bound of SVRG-LD would be equal to
SGHMC (see (6) and (8)). We then omit the same terms (i.e., second and third terms) in the RHS of (7) and (8). Then
we have the following lemma.
√
2
2 2
2 2
2 2
2d
Lemma 1 Let RHMC = L nbTK h D−Ln2 nG2 K+D
(i.e., the first term in the RHS of (7)) and RLD =
2 h3 b−1
(i.e., the first term in the RHS of (8)). Then the following inequality holds.

RHMC ≤ max

n 1
L o LD
,
R
D2 nK
4

n2 K 2 (n2 L2 h2 G2 +hd)
bT

(9)

√
In particular, if D ≥ 1/L h, then (9) becomes:
RHMC ≤

L LD
R
nK

(10)

Note that K is suggested to be 2n by [Johnson and Zhang, 2013] or n/b by [Dubey et al., 2016]. we obtain the
following corollary from Lemma 1.
Corollary 3 If K is n/b as suggested by [Dubey et al., 2016], then (10) becomes:
RHMC ≤

bL LD
R
n2

(11)

2

n
In other words, the SVRG-HMC is O( bL
) times faster than SVRG-LD, in terms of the convergence bound related to
the variance reduction (i.e., the first terms in RHS of (7) and (8)). Note that n is the size of the dataset which can be
very large, b is the mini-batch size which is usually a small constant and L is the Lipschitz smooth parameter for f (θ).

We also want to mention that the convergence proof for SVRG-HMC (i.e., Theorem 2) is a bit more difficult than that
for SVRG-LD [Dubey et al., 2016] due to the momentum variable p (see Line 9 of Algorithm 1).

3.2

SAGA-HMC

In this subsection, we propose the SAGA-HMC algorithm by applying the SAGA framework [Defazio et al., 2014]
to the Hamiltonian Monte Carlo. The details are described in Algorithm 2. Similar to the SVRG-HMC, we initialize θ0 , p0 to be zero vectors in the algorithm only to simplify the analysis; it would still work with an arbitrary
initialization.
Algorithm 2: SAGA-HMC

11

parameters T, b, h > 0, Dh < 1, D ≥ 1;
initialize θ0 = 0, p0 = 0;
i
initialize an array
Pn α0 = θ0 , ∀i ∈ [n];
compute g = i=1 ∇fi (α0i );
for t = 0, 1, . . . , T − 1 do
uniformly randomly pick a set I ⊆ [n] such that |I| = b;
˜ t = −∇ log Pr(θt ) + n P (∇fi (θt ) − ∇fi (αti )) + g;
∇
i∈I
b
√
˜ t + 2Dhξt ;
pt+1 = (1 − Dh)pt − h∇
θt+1 = θt + hpt+1 ;
i
update αt+1
P = θt , ∀i ∈i I;
g ← g + i∈I (∇fi (αt+1
) − ∇fi (αti ));

12

return {θt }Tt=1 ;

1
2
3
4
5
6
7
8
9
10

The following theorem shows the convergence result for MSE of SAGA-HMC. The proof is deferred to Appendix
A.
Theorem 4 Under Assumptions 1 and 2, the MSE of SAGA-HMC is bounded as:
n n2 G2 L2 n4 h2  √n2 G2 + D2 d 2 o
1
2
,
+
+ h2
E[(φ̂ − φ̄) ] . min
bT
T b3
D − L2 n4 h3 b−3
Th
Note that SAGA-HMC can be compared with SAGA-LD [Dubey et al., 2016] in a very similar manner to SVRGHMC (i.e., Lemma 1 and Corollary 3). Thus we omit such a repetition.

5

4

Variance-reduced SGHMC with Symmetric Splitting

Symmetric splitting is a numerically efficient method introduced in [Leimkuhler and Shang, 2016] to accelerate the
gradient-based algorithms. We note that one additional advantage of SGHMC over SGLD is that SGHMC can be
combined with symmetric splitting while SGLD cannot [Chen et al., 2015]. So it is quite natural to combine symmetric
splitting with the proposed SVRG-HMC and SAGA-HMC respectively to see if any further improvements can be
obtained.
The symmetric splitting scheme breaks the original recursion into 5 steps:
h
pt
2
h
(1)
pt = pt − Dpt
2
√
(2)
(1)
˜ t + 2Dhξt
pt = pt − h∇
h (2)
(2)
pt+1 = pt − Dpt
2
h
(1)
θt+1 = θt + pt+1
2
(1)

θt

= θt +

(12)
(13)
(14)
(15)
(16)

If we eliminate the intermediate variables, then

√
Dh
Dh 
˜ t + 2Dhξt
) (1 −
)pt − h∇
2
2
h
h
= θt + pt+1 + pt
2
2

pt+1 = (1 −

(17)

θt+1

(18)

˜ t is computed at θt(1) (which is θt + h pt ) instead
Same as before ξt ∼ N(0, Id ). Note that the stochastic gradient ∇
2
of θt (see (1), (12) and (14)). As shown in [Chen et al., 2015], this symmetric splitting scheme is a 2nd-order local
integrator. Then it improves the dependency of MSE on step size h, i.e., the third term in the RHS of (5) changes to
be h4 , which is a higher order term than the original h2 . It means that we can allow larger step size h by using this
symmetric splitting scheme (note that h < 1).
Similarly, we can further improve the convergence results for SVRG/SAGA-HMC by combining the symmetric
splitting scheme. We give the details of the algorithms and theoretical results for SVRG2nd-HMC and SAGA2ndHMC in the following subsections.

Algorithm 3: SVRG2nd-HMC

10

parameters T, K, b, h > 0, Dh < 1, D ≥ 1;
initialize θ0 = p0 = 0;
for t = 0, 1, . . . , TP
/K − 1 do
n
compute g = i=1 ∇fi (θtK + h2 ptK );
w = θtK + h2 ptK ;
for k = 0, 1, . . . , K − 1 do
uniformly sample an index subset I ⊆ [n], |I| = b;
˜ tK+k = −∇ log Pr(θtK+k + h ptK+k ) + n P (∇fi (θtK+k + h ptK+k ) − ∇fi (w)) + g;
∇
i∈I
2
b
2
√
Dh
˜
ptK+k+1 = (1 − Dh
2DhξtK+k );
2 )((1 − 2 )pt − h∇tK+k +
θtK+k+1 = θtK+k + h2 ptK+k+1 + h2 ptK+k ;

11

return {θt }Tt=1 ;

1
2
3
4
5
6
7
8
9

6

4.1

SVRG2nd-HMC

In this subsection, we propose the SVRG2nd-HMC algorithm (see Algorithm 3) by combining our SVRG-HMC
(Algorithm 1) with the symmetric splitting scheme.
The convergence result for SVRG2nd-HMC is provided in Theorem 5. It shows that the dependency of MSE on
step size h can be improved from h2 to h4 (see (7) and (19)).
Theorem 5 Under Assumptions 1 and 2, the MSE of SVRG2nd-HMC is bounded as:
n n2 G2 L2 n2 K 2 h2  √n2 G2 + D2 d 2 o
1
2
+
,
+ h4
E(φ̂ − φ̄) . min
2
2
2
3
−1
bT
bT
D−L n K h b
Th

4.2

(19)

SAGA2nd-HMC

In this subsection, we propose the SAGA2nd-HMC algorithm (see Algorithm 4) by combining our SAGA-HMC
(Algorithm 2) with the symmetric splitting scheme. The convergence result and algorithm details are described below.
Algorithm 4: SAGA2nd-HMC

11

parameters T, b, h > 0, Dh < 1, D ≥ 1;
initialize θ0 = 0, p0 = 0;
i
initialize an array
∈ [n];
Pn α0 = 0, ∀i
compute g = i=1 ∇fi (α0i );
for t = 0, 1, . . . , T − 1 do
uniformly randomly pick a set I ⊆ [n] such that |I| = b;
˜ t = −∇ log Pr(θt + h pt ) + n P (∇fi (θt + h pt ) − ∇fi (αti )) + g;
∇
i∈I √
2
b
2
Dh
˜ t + 2Dhξt );
pt+1 = (1 − Dh
)((1
−
)p
−
h
∇
t
2
2
θt+1 = θt + h2 pt+1 + h2 pt ;
h
i
update αt+1
P = θt + 2 pit , ∀i ∈ I;
g ← g + i∈I (∇fi (αt+1 ) − ∇fi (αti ));

12

return {θt }Tt=1 ;

1
2
3
4
5
6
7
8
9
10

Theorem 6 Under Assumptions 1 and 2, the MSE of SAGA2nd-HMC is bounded as:
n n2 G2 L2 n4 h2  √n2 G2 + D2 d 2 o
1
2
E(φ̂ − φ̄) . min
,
+
+ h4
bT
T b3
D − L2 n4 h3 b−3
Th

5

Experiment

We present experimental results in this section. We compare the proposed SVRG-HMC (Algorithm 1), as well as as
its symmetric splitting variant SVRG2nd-HMC (Algorithm 5), against SVRG-LD [Dubey et al., 2016] on Bayesian
linear and logistic regression. The experimental results of SAGA variants are almost same as the SVRG variants. In
accordance with the theoretical analysis, all algorithms have fixed step size h, and all HMC-based algorithms have
fixed friction; a grid search is performed to select the best step size and friction for each algorithm. The minibatch size
b is chosen to be 10 for all algorithms, and K is set to be n/b.
The experiments are tested on the real-world UCI datasets1 . The information of the standard datasets used in our
experiments are described in the following Table 1. For each dataset (regression or classification), we partition the
dataset into training (70%), validation (10%) and test (10%) sets. The validation set is used to select step size as well
as friction for HMC-based algorithms in an 8-fold manner.
1 The

UCI datasets can be downloaded from https://archive.ics.uci.edu/ml/datasets.html

7

Table 1: Summary of standard UCI datasets for Bayesian regression and classification
datasets concrete noise parkinson
bike
pima diabetic
eeg
size
1030
1503
5875
17379 768
1151
14980
features
8
5
21
12
8
20
15
The Bayesian regression experiments were conducted on the first four UCI regression datasets and the Bayesian
classification experiments were conducted on the last three UCI classification datasets.

5.1

Bayesian Regression

In this subsection we study the performance of those aforementioned algorithms on Bayesian linear regression. Say
we are provided with inputs Z = {(xi , yi )}ni=1 where xi ∈ Rd and yi ∈ R. The distribution of yi given xi is modelled
as Pr(yi |xi ) = N (β | xi , σ), where the unknown parameter β follows a prior distribution of N(0, Id ). The gradients of
log-likelihood can thus be calculated as ∇β log Pr(yi |xi , β) = (yi − β | xi )xi and ∇β log Pr(β) = −β. The average
test Mean-Squared Error (MSE) is reported in Figure 1.

noise

concrete

Test MSE
0.0

0.5

1.0

1.5

2.0

2.5

Number of pass through data

102

0

parkinson

6 × 101

1

2

3

4

5

Number of pass through data
bike

105

SVRG-LD
SVRG-HMC
SVRG2nd-HMC

104

Test MSE

Test MSE

103

3.0

SVRG-LD
SVRG-HMC
SVRG2nd-HMC

SVRG-LD
SVRG-HMC
SVRG2nd-HMC

104

Test MSE

SVRG-LD
SVRG-HMC
SVRG2nd-HMC

103

103

102

101
0.0

0.5

1.0

1.5

2.0

2.5

Number of pass through data

3.0

0.0

0.2

0.4

0.6

0.8

Number of pass through data

1.0

Figure 1: The x-axis and y-axis represent number of passes through the entire training dataset and average test MSE
respectively. For the bike dataset, we have omitted the first 10 MSE values from the diagram because otherwise the
diagram would scale badly as MSE values are very large in the first several iterations.

8

As can be observed from all five Bayesian regression experiments above, SVRG-HMC as well its symmetric
splitting counterpart SVRG2nd-HMC, converge markedly faster than SVRG-LD in the first pass through the whole
dataset.

5.2

Bayesian Classification

0.650
0.675
0.700
0.725

SVRG2nd-HMC
SVRG-HMC
SVRG-LD

0.750
0.775
0.800
0.0

0.5

1.0

1.5

2.0

2.5

Number of pass through data

3.0

eeg

diabetic

0.60

0.68750

0.62
0.64
0.66
0.68

SVRG2nd-HMC
SVRG-HMC
SVRG-LD

0.70
0.72
0.0

0.5

1.0

1.5

2.0

2.5

Number of pass through data

3.0

Average Test log-likelihood

pima

0.625

Average Test log-likelihood

Average Test log-likelihood

In this subsection we study classification tasks using Bayesian logistic regression. Suppose input data Z = {(xi , yi )}
where xi ∈ Rd , yi ∈ {0, 1}. The distribution of the output yi is modelled as Pr(yi = 1) = 1/(1 + exp(−β | xi )),
where the model parameter follows a prior distribution of N(0, Id ). Then
 the gradient of log-likelihood and log-prior
can be written as ∇β log Pr(yi |xi , β) = yi − 1/(1 + exp(−β | xi )) xi and ∇β log Pr(β) = −β. The average test
log-likelihood is reported in the Figure 2.

0.68775
0.68800
0.68825

SVRG2nd-HMC
SVRG-HMC
SVRG-LD

0.68850
0.68875
0.68900
0.0

0.2

0.4

0.6

0.8

Number of pass through data

1.0

Figure 2: Performance comparison of SVRG-HMC, SVRG2nd-HMC, SVRG-LD on classification tasks. The x-axis
and y-axis represent number of passes through the entire training dataset and average test log-likelihood respectively.
Similar to the Bayesian regression, SVRG-HMC as well its symmetric splitting counterpart SVRG2nd-HMC,
converge markedly faster than SVRG-LD for the Bayesian classification tasks. Also, the experimental results suggest
that SVRG2nd-HMC converges more quickly than SVRG-HMC, which is consistent with Theorem 2 and 5.

6

Conclusion

In this paper, we propose four variance-reduced Hamiltonian Monte Carlo algorithms, i.e., SVRG-HMC, SAGAHMC, SVRG2nd-HMC and SAGA2nd-HMC for Bayesian Inference. These proposed algorithms guarantee improved
theoretical convergence results and converge markedly faster than the benchmarks in practice. We would like to
note that, our variance-reduced Hamiltonian Monte Carlo samplers are not Markovian procedures, but fortunately our
theoretical analysis does not rely on any properties of Markov processes, and so it does not affect the correctness of
Theorem 2, 4, 5 and 6.
For future work, one thing to do with the comparison between variance-reduced Hamiltonian Monte Carlo and
variance-reduced Langevin dynamics is to try some more sophisticated tasks like Bayesian neural networks. Moreover,
we note that stochastic thermostat [Ding et al., 2014] could outperform both SGLD and SGHMC. It might be interesting to study if a variance-reduced variant of stochastic thermostat could also beat SVRG-LD and SVRG/SVRG2ndHMC both theoretically and experimentally.

Acknowledgments
We would like to thank Chang Liu for helpful discussions.

9

A

Missed Proofs

In this appendix, we provide the detailed proofs for Lemma 1, and Theorem 2, 4, 5 and 6.

A.1

Proof of Theorem 2

According to Corollary 1, we have:
E[(φ̂ − φ̄)2 ] .

T −1
1 X
1
E[k∆t k2 ] +
+ h2
2
T t=0
Th

(20)

˜ t − ∇f (θt ) is the additive error in estimating the full gradient ∇f (θt ). By applying the variance
where ∆t = ∇
PT
reduction technique, we need to upper bound the summation t=1 E[k∆t k2 ].
˜ t , we have:
Unpacking the definition of ∆t and ∇
T
−1
X

E[k∆t k2 ]

t=0

=

T
−1
X

E[k − ∇ log Pr(θt ) +

t=0

=

T
−1
X

i∈I

n

n2 E[k

t=0

≤

−1
2 T
X

n
b

nX
(∇fi (θt ) − ∇fi (θb Kt cK )) + g − ∇f (θt )k2 ]
b

1X
1X
(∇fi (θt ) − ∇fi (θb Kt cK )) −
(∇fj (θt ) − ∇fj (θb Kt cK ))k2 ]
b
n j=1
i∈I

Ei∈[n] k∇fi (θt ) − ∇fi (θb Kt cK )k2

(21)

t=0

The Inequality (21) is due to E(X − EX)2 ≤ EX 2 for any random variable X. In the rightmost summation, index i
is picked uniformly random from [n] = {1, 2, . . . , n}.
Then, we bound the RHS of (21) as follows:
T −1
T −1
n2 X
n2 X 2
2
t
Ei∈[n] k∇fi (θt ) − ∇fi (θb K cK )k ≤
L Ekθt − θb Kt cK k2
b t=0
b t=0

≤

≤

T −1
L2 n2 X
K
b t=0

t−1
X

Ekθj+1 − θj k2

t
j=b K
cK

T −1
L2 n2 K 2 X
Ekθt+1 − θt k2
b
t=0

The first inequality is by L-smoothness of all fi ’s, and the second one is by Cauchy’s inequality.
By our algorithm, kθt+1 − θt k2 = h2 kpt+1 k2 , so we need to upper bound Ekpt+1 k2 for each 0 ≤ t < T .
By the recursion of pt+1 ,
Ekpt+1 k2
˜t +
=Ek(1 − Dh)pt − h∇

√

2Dhξt k2

=Ek(1 − Dh)pt − h∇f (θt ) − h∆t +

√

2Dhξt k2

=Ek(1 − Dh)pt − h∇f (θt )k2 + h2 E[k∆t k2 ] + 2Dhd
p
≤(1 − Dh)2 Ekpt k2 + 2Gnh(1 − Dh) Ekpt k2 + h2 n2 G2
+ h2 E[k∆t k2 ] + 2Dhd
10

(22)

The third equality holds because E[∆
pt ] = E[ξt ] = 0 and ∆t and ξt are independent. The first inequality takes
advantage of k∇f k ≤ nG and Ekpt k ≤ Ekpt k2 .
PT
Define S = t=1 E[kpt k2 ]. Then, taking a grand summation over t = 0, 1, . . . , T − 1,
S ≤ (1 − Dh)2 S + 2nGh(1 − Dh)

T
−1 p
X

Ekpt k2

t=0

+ T (h2 n2 G2 + 2Dhd) + h2

T
−1
X

E[k∆t k2 ]

t=0

√ √
≤ (1 − Dh) S + 2nGh(1 − Dh) T S
2

L2 n2 K 2 h4
S
b
The second inequality again contains an implicit Cauchy’s inequality.
Rearranging the terms we have:
r
S
L2 n2 K 2 h4 S
2
(1 − (1 − Dh) −
) − 2nGh(1 − Dh)
− (h2 n2 G2 + 2Dhd) ≤ 0
b
T
T
p
Solving a quadratic equation with respect to S/T and ignoring constant factors, we have:
r
√
√
nG + n2 G2 + D2 d
S
n2 G2 + D2 d
.
.
T
D − L2 n2 K 2 h3 b−1
D − L2 n2 K 2 h3 b−1
+ T (h2 n2 G2 + 2Dhd) +

(23)

From (21) and (22), we have:
T
−1
X

E[k∆t k2 ] ≤

t=0

Recall that kθt+1 − θt k2 = h2 kpt+1 k2 and S =

T −1
L2 n2 K 2 X
Ekθt+1 − θt k2
b
t=0

PT

t=1

E[kpt k2 ], we have:

T
L2 n2 K 2 h2
1X
E[k∆t k2 ] .
T t=1
b

√

n2 G2 + D2 d
D − L2 n2 K 2 h3 b−1

!2
(24)

On the other hand, we can bound (21) as follows:
T
−1
X

E[k∆t k2 ] ≤

t=0

≤

T −1
n2 X
Ei∈[n] k∇fi (θt ) − ∇fi (θb Kt cK )k2
b t=0

4T n2 G2
b

(25)

where (25) holds due to all fi ’s are G-Lipschitz and Cauchy’s inequality ka + bk2 ≤ 2(kak2 + kbk2 ).
Now, the proof of Theorem 2 is finished by combining (20), (24) and (25).

A.2

Proof of Lemma 1

Since G2 > K 2 (n2 L2 h2 G2 + hd) > K 2 n2 L2 h2 G2 , we have h <
D−

1
nKL .

Therefore, we have

h3 L2 n2 K 2
1
L2 n2 K 2
>D− 3 3 3
b
n K L
b
1
=D−
nKbL
 D − 0.1 ≥ 0.9D
11



where the last line holds since nKbL  10 and D > 1 (see Line 1 of Algorithm 1). Thus, the proof is reduced to
2 2 2 2
comparing h2 L2 n2 G2 + hd and h LDn2 G + h2 L2 d asymptotically.
Clearly,
h2 L2 n2 G2
+ h2 L2 d
D2
1
h2 L2 n2 G2
+
≤
hL2 d
D2
nKL
L
h2 L2 n2 G2
+ hd
=
D2
nK
L
1
≤ max{ 2 ,
}(h2 L2 n2 G2 + hd)
D nK
√
L
L
} turns to be nK
which is very small. The reason is that:
Note that if D ≥ 1/L h, then max{ D12 , nK
1
L2 h
1
≥ 2 1
L nKL
nK
=
L

D2 ≥

where the second inequality holds due to h <

A.3

1
nKL

which is mentioned above.



Proof of Theorem 4

˜ t − ∇f (θt ), it suffices to upper bound PT −1 E[k∆t k2 ] according to (20). Unpacking the definition
Defining ∆t = ∇
t=0
˜ t , we have:
of ∆t and ∇
T
−1
X

E[k∆t k2 ]

t=0

=

T
−1
X

˜ t − ∇f (θt )k2 ]
E[k∇

t=0

=

T
−1
X

E[k

t=0

=

i∈I

T
−1
X

n2 E[k

t=0

≤

≤

n
X
1X
(∇fi (θt ) − ∇fi (αti )) −
(∇fj (θt ) − ∇fj (αtj ))k2 ]
b
j=1
i∈I

−1
2 T
X

n
b

n
X
nX
(∇fi (θt ) − ∇fi (αti )) −
(∇fj (θt ) − ∇fj (αtj ))k2 ]
b
j=1

Ei∈[n] [k∇fi (θt ) − ∇fi (αti )k2 ]

(26)

t=0

T −1
L2 n2 X
Ei∈[n] [kθt − αti k2 ]
b t=0

The first inequality is because E(X − EX)2 ≤ EX 2 for any random variable X; the second inequality holds due to
L-smoothness of fi ’s.

12

Let γ = 1 − (1 − 1/n)b . Next we upper bound each E[kθt − αti k2 ] in the following manner.
E[kθt − αti k2 ] =

t−1
X

E[kθt − θj k2 ] Pr(αti = θj )

j=0

=

t−1
X

E[kθt − θj k2 ](1 − γ)t−j−1 γ

j=0

= h2

t−1
X

E[kpt + pt−1 + · · · + pj+1 k2 ](1 − γ)t−j−1 γ

j=0
t−1
X
≤ h2
(t − j)(1 − γ)t−j−1 γ(E[kpt k2 ] + E[kpt−1 k2 ] + · · · + E[kpj+1 k2 ])
j=0

= h2

t
X

E[kpj k2 ]

j=1

≤ h2

t
X

< h2

(t − k)(1 − γ)t−k−1 γ

k=0

E[kpj k2 ]

j=1
t
X

j−1
X

∞
X

(t − j + k + 1)(1 − γ)t−j+k γ

k=0

E[kpj k2 ](

j=1

1
+ t − j)(1 − γ)t−j
γ

The second equality is by direct calculation Pr(αti = θj ) = (1 − γ)t−j−1 γ; the first inequality is a direct application
of Cauchy’s inequality; the last inequality is a weighted summation of geometric series (1 − γ)t−j+k , k ≥ 0.
Summing over all t and i, we then have:
T
−1
X

Ei∈[n] [kθt − αti k2 ] ≤ h2

t=0

t
T
−1 X
X

E[kpj k2 ](

1
+ t − j)(1 − γ)t−j
γ

E[kpj k2 ](

1
+ t − j)(1 − γ)t−j
γ

t=0 j=1

= h2

t
T
−1 X
X
t=0 j=1

≤ h2

T
−1
X

E[kpt k2 ]

t=1

≤ h2

T
−1
X

≤ h2

(

j=t

1
+ j − t)(1 − γ)j−t
γ

E[kpt k2 ]

∞
X
1
( + j − t)(1 − γ)j−t
γ
j=t

E[kpt k2 ]

2
γ2

t=1
T
−1
X

T
−1
X

t=1
T −1

The last inequality is because (1 − 1/n)b ≤

=

2 2X
h
E[kpt k2 ]
γ2
t=1

≤

T
8h2 n2 X
E[kpt k2 ]
b2 t=1

1
b
1+ n−1

, and thus γ = 1 − (1 − 1/n)n ≥

holds as mini-batch size b is smaller than dataset size n.

13

b
n−1
b
1+ n−1

>

b
2n ;

the last inequality

PT
Similar to the previous subsection, we derive upper an upper bound on t=1 E[kpt k2 ]. By recursion of pt+1 ’s, we
have:
√
˜ t + 2Dhξt k2
Ekpt+1 k2 = Ek(1 − Dh)pt − h∇
√
= Ek(1 − Dh)pt − h∇f (θt ) − h∆t + 2Dhξt k2
= Ek(1 − Dh)pt − h∇f (θt )k2 + h2 E[k∆t k2 ] + 2Dhd
p
≤ (1 − Dh)2 Ekpt k2 + 2Gnh(1 − Dh) Ekpt k2 + h2 n2 G2
+ h2 E[k∆t k2 ] + 2Dhd
The third equality holds because E[∆t ]p= E[ξt ] = 0 and ∆t and ξt are independent. The first inequality takes
advantage of k∇f k ≤ nG and Ekpt k ≤ Ekpt k2 .
PT
Define S = t=1 E[kpt k2 ]. Then, taking a grand summation over t = 0, 1, . . . , T − 1,
2

S ≤ (1 − Dh) S + 2nGh(1 − Dh)

T
−1 p
X

Ekpt k2

t=0

+ T (h2 n2 G2 + 2Dhd) + h2

T
−1
X

E[k∆t k2 ]

t=0

√ √
≤ (1 − Dh) S + 2nGh(1 − Dh) T S
2

+ T (h2 n2 G2 + 2Dhd) +

8h4 L2 n4
S
b3

Rearranging the terms we have:
r
S
8h4 L2 n4 S
(1 − (1 − Dh) −
) + 2nGh(1 − Dh)
+ (h2 n2 G2 + 2Dhd) ≤ 0
3
b
T
T
p
Solving a quadratic equation with respect to S/T and ignoring constant factors, we have:
r
√
√
S
nG
n2 G2 + D2 d
n2 G2 + D2 d
.
+
.
T
D − h3 L2 n4 b−3
D − h3 L2 n4 b−3
D − h3 L2 n4 b−3
2

Plugging it in
T
−1
X

E[k∆t k2 ] ≤

t=0

T −1
8h2 L2 n4 X
Ekpt k2
b3
t=0

we have:
T
h2 L2 n4
1X
E[k∆t k2 ] .
T t=1
b3

√

n2 G2 + D 2 d
D − h3 L2 n4 b−3

!2
(27)

Similar to (25), we can bound (26) as follows:
T
−1
X

E[k∆t k2 ]

t=0

≤

T −1
n2 X
Ei∈[n] [k∇fi (θt ) − ∇fi (αti )k2 ]
b t=0

≤

4T n2 G2
b

(28)

where (28) holds due to all fi ’s are G-Lipschitz and Cauchy’s inequality ka + bk2 ≤ 2(kak2 + kbk2 ).
Now, the proof of Theorem 4 is finished by combining (20), (27) and (28).

14

A.4

Proof of Theorem 5

To prove this theorem, we need the following theorem from [Chen et al., 2015]. The theorem shows that SGHMC
with symmetric splitting can improve the
˜ t be an unbiased estimate of ∇f (θt ) for all t. Then under Assumption 1, for
Theorem 7 ([Chen et al., 2015]) Let ∇
a smooth test function φ, the MSE of SGHMC with symmetric splitting is bounded in the following way:
PT −1
1
2
1
t=0 E(∆Vt ψ(θt , pt ))
2
T
+
+ h4
E(φ̂ − φ̄) .
(29)
T
Th
˜ t − ∇f (θt + h pt ). According to Assumption 2, we have:
Now, we define ∆t = ∇
2
E(φ̂ − φ̄)2 ≤

T −1
1
1 X
Ek∆t k2 +
+ h4
T 2 t=0
Th

According to (30), we mainly need to bound the term
˜ t,
we unfold the definition of ∇

PT −1
t=0

(30)

Ek∆t k2 for our SVRG2nd-HMC algorithm. First,

X
h
h
˜ t = −∇ log Pr(θt + h pt ) + n
(∇fi (θt + pt ) − ∇fi (θb Kt cK + pb Kt cK ))
∇
2
b
2
2
i∈I

+

n
X

∇fi (θb Kt cK

i=1

h
+ pb Kt cK )
2

Then,
Ek∆t k2 = Ek

h
nX
h
(∇fi (θt + pt ) − ∇fi (θb Kt cK + pb Kt cK ))
b
2
2
i∈I
n
X

−

(∇fi (θt +

i=1

h
h
pt ) − ∇fi (θb Kt cK + pb Kt cK ))k2
2
2

h
nX
h
≤ Ek
(∇fi (θt + pt ) − ∇fi (θb Kt cK + pb Kt cK ))k2
b
2
2
i∈I

n2
h
h
Ei∈[n] k∇fi (θt + pt ) − ∇fi (θb Kt cK + pb Kt cK )k2
≤
b
2
2
n2 L2
h
h
2
≤
Ekθt + pt − θb Kt cK − pb Kt cK k
b
2
2
t−1
2 2
X
n L K
h
h
≤
Ekθj+1 + pj+1 − θj − pj k2
b
2
2
t
j=b K cK

Taking a summation we have:
T
−1
X

Ek∆t k2 ≤

t=0

T −1
n2 L2 K X
b
t=0

t−1
X

Ekθj+1 +

t
j=b K
cK

h
h
pj+1 − θj − pj k2
2
2

≤

T −1
n2 L2 K 2 X
h
h
Ekθt+1 + pt+1 − θt − pt k2
b
2
2
t=0

=

T −1
n2 L2 K 2 h2 X
Ekpt+1 k2
b
t=0

15

(31)

The last equality follows by the recursion θt+1 = θt + h2 pt+1 + h2 pt .
By definition of pt+1 we have:
Ekpt+1 k2

√
Dh
Dh 
˜ t + 2Dhξt k2
) (1 −
)pt − h∇
=Ek(1 −
2
2
√
Dh 2
h
Dh
=(1 −
) Ek(1 −
)pt − h∆t − h∇f (θt + pt ) + 2Dhξt k2
2
2
2

√
Dh 2 
Dh
h
=(1 −
) Ek(1 −
)pt − h∇f (θt + pt )k2 + Ek 2Dhξt k2 + Ekh∆t k2
2
2
2

Dh 2 
Dh 2
Dh
2
=(1 −
) (1 −
) Ekpt k + 2(1 −
)nGEkpt k + n2 G2 h2 + 2Dhd + h2 Ek∆t k2
2
2
2

p
Dh 2 
Dh 2
Dh
2
≤(1 −
) (1 −
) Ekpt k + 2(1 −
)nG Ekpt k2 + n2 G2 h2 + 2Dhd + h2 Ek∆t k2
2
2
2
PT
2
2 2 2
Define S = t=1 Ekpt k and M = n G h + 2Dhd. Taking a grand summation of the above inequality for
t = 0, 1, 2, . . . , T − 1, we have:


T
−1 p
T
−1
X
X
Dh 2
Dh
Dh 2
S ≤ (1 −
Ekpt k2 + T M + h2
Ek∆t k2
) (1 −
) S + 2(1 −
)nG
2
2
2
t=0
t=0


T
−1
X
√
Dh 2
Dh 2
Dh
≤ (1 −
) (1 −
) S + 2(1 −
)nG T S + T M + h2
Ek∆t k2
2
2
2
t=0


T
−1
X
√
Dh 2
Dh 2
Dh
2
2
≤ (1 −
) (1 −
) S + 2(1 −
)nG T S + T M + h
Ek∆t k
2
2
2
t=0


√
Dh 2
Dh
n2 L2 K 2 h4
Dh 2
≤ (1 −
) (1 −
) S + 2(1 −
)nG T S + T M +
S
2
2
2
b
q
Rewriting it as a quadratic inequality with respect to TS , we have:
r

Dh 2 n2 L2 K 2 h4  S
Dh 3
S
Dh 2
Dh 4
) − (1 −
)
− 2(1 −
) nG
− (1 −
) M ≤0
1 − (1 −
2
2
b
T
2
T
2
Solve the inequality and ignore constant factors:
r
√
nG + n2 G2 + D2 d
S
.
T
D − L2 n2 K 2 h3 b−1
√
n2 G2 + D2 d
.
D − L2 n2 K 2 h3 b−1
Similar to the proof of Theorem 2, it easily follows that:
!2
√
T −1
1 X
L2 n2 K 2 h2
n2 G2 + D 2 d
2
Ek∆t k .
(32)
T t=0
b
D − L2 n2 K 2 h3 b−1
Similar to (25), we can bound (31) as follows:
h
h
n2
Ei∈[n] k∇fi (θt + pt ) − ∇fi (θb Kt cK + pb Kt cK )k2
b
2
2
4n2 G2
≤
b
where (33) holds due to all fi ’s are G-Lipschitz and Cauchy’s inequality ka + bk2 ≤ 2(kak2 + kbk2 ).
Now, the proof of Theorem 5 is finished by combining (30), (32) and (33).
E[k∆t k2 ] ≤

16

(33)


A.5

Proof of Theorem 6

˜ t − ∇f (θt + h pt ). By Assumption 2 and Inequality (29), we
Similar to the proof of Theorem 5, we define ∆t = ∇
2
have:
PT −1
1
2
1
t=0 E(∆Vt ψ(θt , pt ))
2
T
+
+ h4
E(φ̂ − φ̄) .
T
Th
T −1
1 X
1
≤ 2
Ek∆t k2 +
+ h4
T t=0
Th
˜ t , we have:
Unpacking the definition of ∆t and ∇
T
−1
X

E[k∆t k2 ] =

t=0

T
−1
X

˜ t − ∇f (θt )k2 ]
E[k∇

t=0

=

T
−1
X

E[k

t=0

=

i∈I

T
−1
X

n2 E[k

t=0

≤

n 
 X

h
h
1 X
∇fi (θt + pt ) − ∇fi (αti ) −
∇fj (θt + pt ) − ∇fj (αtj ) k2 ]
b
2
2
j=1
i∈I

−1
2 T
X

n
b

n 
 X

h
h
n X
∇fi (θt + pt ) − ∇fi (αti ) −
∇fj (θt + pt ) − ∇fj (αtj ) k2 ]
b
2
2
j=1

Ei∈[n] [k∇fi (θt +

t=0

h
pt ) − ∇fi (αti )k2 ]
2

T −1
L2 n2 X
h
≤
Ei∈[n] [kθt + pt − αti k2 ]
b t=0
2

Let γ = 1 − (1 − 1/n)b . Then,
E[kθt +
=

t−1
X

h
pt − αti k2 ]
2

E[kθt +

h
h
h
pt − θj − pj k2 ] Pr(αti = θj + pj )
2
2
2

E[kθt +

h
h
pt − θj − pj k2 ](1 − γ)t−j−1 γ
2
2

j=0

=

t−1
X
j=0
2

=h

t−1
X

E[kpt + pt−1 + · · · + pj+1 k2 ](1 − γ)t−j−1 γ

j=0
t−1
X
≤h
(t − j)(1 − γ)t−j−1 γ(E[kpt k2 ] + E[kpt−1 k2 ] + · · · + E[kpj+1 k2 ])
2

j=0
2

=h

t
X

2

E[kpj k ]

j=1

≤ h2

t
X
j=1

< h2

t
X
j=1

j−1
X

(t − k)(1 − γ)t−k−1 γ

k=0

E[kpj k2 ]

∞
X

(t − j + k + 1)(1 − γ)t−j+k γ

k=0

E[kpj k2 ](

1
+ t − j)(1 − γ)t−j
γ
17

(34)

(35)

Summing over all t and i, we then have:
T
−1
X

T −1

Ei∈[n] [kθt +

t=0

t

XX
h
1
pt − αti k2 ] ≤ h2
E[kpj k2 ]( + t − j)(1 − γ)t−j
2
γ
t=0 j=1
= h2

T
−1 X
t
X

E[kpj k2 ](

t=0 j=1

≤ h2

T
−1
X

E[kpt k2 ]

t=1

≤ h2

T
−1
X

≤ h2

t=1

(

j=t

1
+ j − t)(1 − γ)j−t
γ

E[kpt k2 ]

∞
X
1
( + j − t)(1 − γ)j−t
γ
j=t

E[kpt k2 ]

X
2
2
= 2 h2
E[kpt k2 ]
2
γ
γ
t=1

t=1
T
−1
X

T
−1
X

1
+ t − j)(1 − γ)t−j
γ

T −1

T
8h n X
E[kpt k2 ]
b2 t=1
2 2

≤
The last inequality is because (1 − 1/n)b ≤

1
b
1+ n−1

, and thus γ = 1 − (1 − 1/n)n ≥

(36)
b
n−1
b
1+ n−1

>

b
2n ;

the last

inequality holds as mini-batch size b is P
smaller than dataset size n.
T
Now, we derive an upper bound on t=1 E[kpt k2 ]. By recursion of pt+1 ’s, we have:
Ekpt+1 k2
√
Dh
Dh
˜ t + 2Dhξt )k2
)((1 −
)pt − h∇
2
2
√
Dh 2
Dh
h
= (1 −
) Ek(1 −
)pt − h∆t − h∇f (θt + pt ) + 2Dhξt k2
2
2
2

√
Dh 2 
h
Dh
= (1 −
) Ek(1 −
)pt − h∇f (θt + pt )k2 + Ek 2Dhξt k2 + Ekh∆t k2
2
2
2

Dh 2 
Dh 2
Dh
2
= (1 −
) (1 −
) Ekpt k + 2(1 −
)nGEkpt k + n2 G2 h2 + 2Dhd + h2 Ek∆t k2
2
2
2

p
Dh 2 
Dh 2
Dh
2
≤ (1 −
) (1 −
) Ekpt k + 2(1 −
)nG Ekpt k2 + n2 G2 h2 + 2Dhd + h2 Ek∆t k2
2
2
2
PT
2
2 2 2
Define S = t=1 Ekpt k and M = n G h + 2Dhd. Taking a grand summation of the above inequality for
t = 0, 1, 2, . . . , T − 1, we have:

T
−1 p
X
Dh 2
Dh 2
Dh
S ≤ (1 −
) (1 −
) S + 2(1 −
)nG
Ekpt k2
2
2
2
t=0

T
−1
X
+ T M + h2
Ek∆t k2
= Ek(1 −

t=0
T
−1
X
√
Dh 2
Dh
(1 −
) S + 2(1 −
)nG T S + T M + h2
Ek∆t k2
2
2
t=0

!

T
−1
X
√
Dh 2
Dh
(1 −
) S + 2(1 −
)nG T S + T M + h2
Ek∆t k2
2
2
t=0


√
Dh 2
Dh 2
Dh
8n4 L2 h4
≤ (1 −
) (1 −
) S + 2(1 −
)nG T S + T M +
S
2
2
2
b3

!

Dh 2
≤ (1 −
)
2
Dh 2
≤ (1 −
)
2

18

Similar to the proof of theorem 5, we solve a quadratic inequality with respect to
r

q

S
T,

and then,

√
S
n2 G2 + D2 d
.
T
D − h3 L2 n4 b−3

(37)

From (35), (36), (37) and the definition of S, we have:
T
−1
X

E[k∆t k2 ]

t=0
T −1
8h2 L2 n4 X
Ekpt k2
b3
t=0
!2
√
8h2 L2 n4 T
n2 G2 + D2 d
≤
b3
D − h3 L2 n4 b−3

≤

(38)

Similar to (25), we can bound (34) as follows:
T
−1
X

E[k∆t k2 ]

t=0

≤

T −1
h
n2 X
Ei∈[n] [k∇fi (θt + pt ) − ∇fi (αti )k2 ]
b t=0
2

≤

4T n2 G2
b

(39)

where (39) holds due to all fi ’s are G-Lipschitz and Cauchy’s inequality ka + bk2 ≤ 2(kak2 + kbk2 ).
Now, the proof of Theorem 6 is finished by combining (30), (38) and (39).



References
Sungjin Ahn, Anoop Korattikara, and Max Welling. Bayesian posterior sampling via stochastic gradient fisher scoring.
In Proceedings of the 29th International Conference on Machine Learning, pages 1771–1778, 2012.
Sungjin Ahn, Babak Shahbaba, and Max Welling. Distributed stochastic gradient mcmc. In International Conference
on Machine Learning, pages 1044–1052, 2014.
Simon Byrne and Mark Girolami. Geodesic monte carlo on embedded manifolds. Scandinavian Journal of Statistics,
40(4):825–845, 2013.
Changyou Chen, Nan Ding, and Lawrence Carin. On the convergence of stochastic gradient mcmc algorithms with
high-order integrators. In Advances in Neural Information Processing Systems, pages 2278–2286, 2015.
Tianqi Chen, Emily Fox, and Carlos Guestrin. Stochastic gradient hamiltonian monte carlo. In International Conference on Machine Learning, pages 1683–1691, 2014.
Aaron Defazio, Francis Bach, and Simon Lacoste-Julien. Saga: A fast incremental gradient method with support for
non-strongly convex composite objectives. In Advances in Neural Information Processing Systems, pages 1646–
1654, 2014.
Nan Ding, Youhan Fang, Ryan Babbush, Changyou Chen, Robert D Skeel, and Hartmut Neven. Bayesian sampling
using stochastic gradient thermostats. In Advances in Neural Information Processing Systems, pages 3203–3211,
2014.
19

Simon Duane, Anthony D Kennedy, Brian J Pendleton, and Duncan Roweth. Hybrid monte carlo. Physics letters B,
195(2):216–222, 1987.
Kumar Avinava Dubey, Sashank J Reddi, Sinead A Williamson, Barnabas Poczos, Alexander J Smola, and Eric P
Xing. Variance reduction in stochastic gradient langevin dynamics. In Advances in Neural Information Processing
Systems, pages 1154–1162, 2016.
Mark Girolami and Ben Calderhead. Riemann manifold langevin and hamiltonian monte carlo methods. Journal of
the Royal Statistical Society: Series B (Statistical Methodology), 73(2):123–214, 2011.
Rie Johnson and Tong Zhang. Accelerating stochastic gradient descent using predictive variance reduction. In Advances in Neural Information Processing Systems, pages 315–323, 2013.
Benedict Leimkuhler and Xiaocheng Shang. Adaptive thermostats for noisy gradient systems. SIAM Journal on
Scientific Computing, 38(2):A712–A736, 2016.
Chang Liu, Jun Zhu, and Yang Song. Stochastic gradient geodesic mcmc methods. In Advances in Neural Information
Processing Systems, pages 3009–3017, 2016.
Yi-An Ma, Tianqi Chen, and Emily Fox. A complete recipe for stochastic gradient mcmc. In Advances in Neural
Information Processing Systems, pages 2917–2925, 2015.
Radford M Neal et al. Mcmc using hamiltonian dynamics. Handbook of Markov Chain Monte Carlo, 2(11), 2011.
Sam Patterson and Yee Whye Teh. Stochastic gradient riemannian langevin dynamics on the probability simplex. In
Advances in Neural Information Processing Systems, pages 3102–3110, 2013.
Herbert Robbins and Sutton Monro. A stochastic approximation method. The annals of mathematical statistics, pages
400–407, 1951.
Max Welling and Yee W Teh. Bayesian learning via stochastic gradient langevin dynamics. In Proceedings of the 28th
International Conference on Machine Learning, pages 681–688, 2011.

20

