Bayesian Learning via Stochastic Gradient Langevin Dynamics

Max Welling
welling@ics.uci.edu
D. Bren School of Information and Computer Science, University of California, Irvine, CA 92697-3425, USA
Yee Whye Teh
ywteh@gatsby.ucl.ac.uk
Gatsby Computational Neuroscience Unit, UCL, 17 Queen Square, London WC1N 3AR, UK

Abstract
In this paper we propose a new framework
for learning from large scale datasets basedon iterative learning from small mini-batches.
By adding the right amount of noise to a
standard stochastic gradient optimization algorithm we show that the iterates will converge to samples from the true posterior distribution as we anneal the stepsize. Thisseamless transition between optimization andBayesian posterior sampling provides an inbuilt protection against overﬁtting. We also
propose a practical method for Monte Carloestimates of posterior statistics which monitors a “sampling threshold” and collects samples after it has been surpassed. We apply
the method to three models: a mixture ofGaussians, logistic regression and ICA withnatural gradients.

1. Introduction
In recent years there has been an increasing amountof very large scale machine learning datasets, ranging
from internet traﬃc and network data, computer vision, natural language processing, to bioinformatics.
More and more advances in machine learning are now
driven by these large scale data, which oﬀers the opportunity to learn large and complex models for solving many useful applied problems. Recent successes
in large scale machine learning have mostly been optimization based approaches. While there are sophisticated algorithms designed speciﬁcally for certain typesof models, one of the most successful class of algorithms are stochastic optimization, or Robbins-Monro,
algorithms. These algorithms process small (miniAppearing in Proceedings of the 28 th International Conference on Machine Learning, Bellevue, WA, USA, 2011.
Copyright 2011 by the author(s)/owner(s).

)batches of data at each iteration, updating model
parameters by taking small gradient steps in a costfunction. Often these algorithms are run in an online setting, where the data batches are discarded after processing and only one pass through the data is
performed, reducing memory requirements drastically.
One class of methods “left-behind” by the recent advances in large scale machine learning are the Bayesian
methods. This has partially to do with the negative
results in Bayesian online parameter estimation (Andrieu et al., 1999), but also the fact that each iterationof typical Markov chain Monte Carlo (MCMC) algorithms requires computations over the whole dataset.
Nevertheless, Bayesian methods are appealing in theirability to capture uncertainty in learned parameters
and avoid overﬁtting. Arguably with large datasets
there will be little overﬁtting. Alternatively, as wehave access to larger datasets and more computational
resources, we become interested in building more complex models, so that there will always be a need toquantify the amount of parameter uncertainty.
In this paper, we propose a method for Bayesian learning from large scale datasets. Our method combines
Robbins-Monro type algorithms which stochastically
optimize a likelihood, with Langevin dynamics which
injects noise into the parameter updates in such a waythat the trajectory of the parameters will converge
to the full posterior distribution rather than just the
maximum a posteriori mode. The resulting algorithm
starts oﬀ being similar to stochastic optimization, thenautomatically transitions to one that simulates samples from the posterior using Langevin dynamics.
In Section 2 we introduce the two ingredients of our
method: stochastic optimization and Langevin dynamics. Section 3 describes our algorithm and how
it converges to the posterior distribution. Section 4
describes a practical method of estimating when ouralgorithm will transition from stochastic optimization
to Langevin dynamics. Section 5 demonstrates our al-

Stochastic Gradient Langevin Dynamics

gorithm on a few models and Section 6 concludes.

not collapse to just the MAP solution:

2. Preliminaries∆θt =Let θ denote a parameter vector, with p(θ) a priordistribution, and p(x|θ) the probability of data item
x given our model parameterized by θ. The posterior
distribution of a set of N data items X = {xi }N
i=1∏N
is: p(θ|X) ∝ p(θ) i=1 p(xi |θ). In the optimization
literature the prior regularizes the parameters while
the likelihood terms constitute the cost function to
be optimized, and the task is to ﬁnd the maximum
a posteriori (MAP) parameters θ∗ . A popular classof methods called stochastic optimization (Robbins &Monro, 1951) operates as follows. At each iteration t,
a subset of n data items Xt = {xt1 , . . . , xtn } is given,and the parameters are updated as follows:
(
)
n
ϵt
N∑
∆θt =
∇ log p(θt ) +∇ log p(xti |θt )
(1)
2
n i=1
where ϵt is a sequence of step sizes. The generalidea is that the gradient computed on the subset isused to approximate the true gradient over the wholedataset. Over multiple iterations the whole dataset
is used and the noise in the gradient caused by using
subsets rather than the whole dataset averages out.
For large datasets where the subset gradient approximation is accurate enough, this can give signiﬁcantcomputational savings over using the whole dataset tocompute gradients at each iteration.
To ensure convergence to a local maximum, in addition
to other technical assumptions, a major requirement
is for the step sizes to satisfy the property∞∑t=1

ϵt = ∞∞∑

ϵ2t < ∞

(2)

t=1

Intuitively, the ﬁrst constraint ensures that parameters
will reach the high probability regions no matter how
far away it was initialized to, while the second ensuresthat the parameters will converge to the mode insteadof just bouncing around it. Typically, step sizes ϵt =a(b + t)−γ are decayed polynomially with γ ∈ (0.5, 1].
The issue with ML or MAP estimation, as stochastic optimization aims to do, is that they do not capture parameter uncertainty and can potentially overﬁtdata. The typical way in which Bayesian approaches
capture parameter uncertainty is via Markov chain
Monte Carlo (MCMC) techniques (Robert & Casella,2004). In this paper we will consider a class of MCMC
techniques called Langevin dynamics (Neal, 2010). As
before, these take gradient steps, but also injects Gaussian noise into the parameter updates so that they do

)
(
N
∑
ϵ
∇ log p(θt ) +∇ log p(xi |θt ) + ηt2
i=1

ηt ∼ N (0, ϵ)(3)

The gradient step sizes and the variances of the injected noise are balanced so that the variance of thesamples matches that of the posterior. Langevin dynamics is motivated and originally derived as a discretization of a stochastic diﬀerential equation whose
equilibrium distribution is the posterior distribution.
To correct for discretization error, one can take (3)to just be a proposal distribution and correct using
Metropolis-Hastings. Interestingly, as we decrease ϵ
the discretization error decreases as well so that the rejection rate approaches zero. However typical MCMC
practice is to allow an initial adaptation phase where
the step sizes are adjusted, followed by ﬁxing the step
sizes to ensure a stationary Markov chain thereafter.
More sophisticated techniques use Hamiltonian dynamics with momentum variables to allow parameters
to move over larger distances without the ineﬃcient
random walk behaviour of Langevin dynamics (Neal,2010). However, to the extent of our knowledge all
MCMC methods proposed thus far require computations over the whole dataset at every iteration, resulting in very high computational costs for large datasets.

3. Stochastic Gradient Langevin
Dynamics
Given the similarities between stochastic gradient algorithms (1) and Langevin dynamics (3), it is natural to consider combining ideas from the two approaches. This allows eﬃcient use of large datasets
while allowing for parameter uncertainty to be captured in a Bayesian manner. The approach is straightforward: use Robbins-Monro stochastic gradients, add
an amount of Gaussian noise balanced with the step
size used, and allow step sizes to go to zero. The proposed update is simply:
(
)
n
ϵt
N∑
∆θt =
∇ log p(θt ) +∇ log p(xti |θt ) + ηt2
n i=1
ηt ∼ N (0, ϵt )(4)where the step sizes decrease towards zero at rates satisfying (2). This allows averaging out of the stochasticity in the gradients, as well as MH rejection rates that
go to zero asymptotically, so that we can simply ignore
the MH acceptance steps, which require evaluation ofprobabilities over the whole dataset, all together.

Stochastic Gradient Langevin Dynamics

In the rest of this section we will give an intuitive argument for why θt will approach samples from the posterior distribution as t → ∞. In particular, we willshow that for large t, the updates (4) will approach
Langevin dynamics (3), which converges to the posterior distribution. Letg(θ) = ∇ log p(θ) +N∑∇ log p(xi |θ)(5)

be the true gradient of the log probability at θ andN∑
ht (θ) = ∇ log p(θ) +∇ log p(xti |θ) − g(θ) (6)n i=1n

The stochastic gradient is then g(θ)+ht (θ), with ht (θ)
a zero mean random variable (due to the stochasticityof the data items chosen at step t) with ﬁnite variance
V (θ), and (4) is,ϵt(g(θt ) + ht (θt )) + ηt ,2

ηt ∼ N (0, ϵt )∑ ϵt(g(θt ) + ht (θt ))2
t=t +1
ts+1

(8)

s

i=1

∆θt =

∥θt − θts ∥2 ≪ 1 for t between ts and ts+1 . Making
the assumption that the gradient g(·) vary smoothly(e.g. they are Lipschitz continuous in the models in
Section 5), the total stochastic gradient is:(7)

There are two sources of stochasticity in (7): the injected Gaussian noise with variance ϵt , and the noise in
the stochastic gradient, which has variance ( ϵ2t )2 V (θt ).
The ﬁrst observation is that for large t, ϵt → 0, andthe injected noise will dominate the stochastic gradientnoise, so that (7) will be eﬀectively Langevin dynamics (3). The second observation is that as ϵt → 0,
the discretization error of Langevin dynamics will be
negligible so that the MH rejection probability will approach 0 and we may simply ignore this step.
In other words, (4), (7) eﬀectively deﬁne a nonstationary Markov chain such that the tth step transition operator, for all large t, will have as its equilibrium distribution the posterior over θ. The next question we address is whether the sequence of parametersθ1 , θ2 , . . . will converge to the posterior distribution.
Because the Markov chain is not stationary and thestep sizes reduce to 0, it is not immediately clear that
this is the case. To see that this is indeed true, wewill show that a subsequence θt1 , θt2 , . . . will converge
to the posterior as intended so the whole sequence will
also converge.
First ﬁx an ϵ0 such that 0 < ϵ0 ≪ 1. Since {ϵt } satisfy
the step size property (2),can ﬁnd a subsequence∑we
ts+1
t1 < t2 < · · · such that t=tϵ → ϵ0 as s → ∞.
s +1 tSince the injected noise at each step is independent,
for
∑ts+1
large enough s the total injected noise, ∥ t=t
η
∥
t 2,
s +1
√
between steps ts and ts+1 will be O( ϵ0 ). We now
show that the total noise due to the stochasticity of
the gradients among these steps will be dominated by
the total injected noise. Since ϵ0 ≪ 1, we may take

ts+1
∑ ϵt
ϵ0
= g(θts ) + O(ϵ0 ) +ht (θt )22
t=t +1s

Since the parameters did not vary much between ts
and ts+1 , the stochasticity in ht (θt ) will be dominatedby the randomness in the choice of the mini-batches.
Assuming that these are chosen randomly and independently, ht (θt ) for each t will be basically iid (ifmini-batches were chosen by random partitioning of
the whole dataset, ht (θt ) will be negatively correlated
instead and will not change the results here). Thus
∑ts+1 ϵt∑ ϵ2the variance of t=t
ht (θt ) is O( t 4t ) ands +1 2)(√∑ts+1 ϵ2t
ϵ0
= g(θts ) + O(ϵ0 ) + O
t=ts +1 4
2
ϵ0= g(θts ) + O(ϵ0 )2The last equation says that the total stochastic gradient step is approximately the exact gradient step at
θts with a step size of ϵ0 , with a deviation dominatedby O(ϵ0 ). Since this is in turn dominated by the total√
injected noise which is O( ϵ0 ), this means that the sequence θt1 , θt2 , . . . will approach a sequence generatedby Langevin dynamics with a ﬁxed step size ϵ0 , so it
will converge to the posterior distribution. Note alsothat it will have inﬁnite eﬀective sample size.
The implication of this argument is that we can use
stochastic gradient Langevin dynamics as an “anytime” and general-purpose algorithm. In the initial
phase the stochastic gradient noise will dominate andthe algorithm will imitate an eﬃcient stochastic gradient ascent algorithm. In the later phase the injected
noise will dominate, so the algorithm will imitate a
Langevin dynamics MH algorithm, and the algorithm
will transition smoothly between the two. Howevera disadvantage is that to guarantee the algorithm towork it is important for the step sizes to decrease tozero, so that the mixing rate of the algorithm willslow down with increasing number of iterations. To
address this, we can keep the step size constant once
it has decreased below a critical level where the MH
rejection rate is considered negligible, or use this algorithm for burn-in, but switch to a diﬀerent MCMC
algorithm that makes more eﬃcient use of the wholedataset later. These alternatives can perform better

Stochastic Gradient Langevin Dynamics

but will require further hand-tuning and are beyond
the scope of this paper. The point of this paper isto demonstrate a practical algorithm that can achieve
proper Bayesian learning using only mini-batch data.

4. Posterior Sampling
In this section we consider the use of our stochastic
gradient Langevin dynamics algorithm as one which
produces samples from the posterior distribution. We
ﬁrst derive an estimate of when the algorithm willtransition from stochastic optimization to Langevindynamics. The idea is that we should only start collecting samples after it has entered its posterior sampling phase, which will not happen until after it becomes Langevin dynamics. Then we discuss how thealgorithm scales with the dataset size N and give arough estimate of the number of iterations required forthe algorithm to traverse the whole posterior. Finally
we discuss how the obtained samples can be used toform Monte Carlo estimates of posterior expectations.
4.1. Transition into Langevin dynamics phase
We ﬁrst generalize our method to allow for preconditioning, which can lead to signiﬁcant speed ups by
better adapting the step sizes to the local structure of
the posterior (Roberts & Stramer, 2002; Girolami &Calderhead, 2011). For instance, certain dimensions
may have a vastly larger curvature leading to muchbigger gradients. In this case a symmetric preconditioning matrix M can transform all dimensions to thesame scale. The preconditioned stochastic gradientLangevin dynamics is simply,()ϵt∆θt = M g(θt ) + ht (θt ) + ηt , ηt ∼ N (0, ϵt M )2
As noted previously, whether the algorithm is in thestochastic optimization phase or Langevin dynamics
phase depends on the variance of the injected noise,
which is simply ϵt M , versus that of the stochastic gradient. Since the stochastic gradient is a sum over thecurrent mini-batch, if its size n is large enough thecentral limit theorem will kick in and the variations
ht (θt ) around the true gradient g(θt ) will become normally distributed. Its covariance matrix can then be
estimated from the empirical covariance:N2 ∑(sti − st )(sti − st )⊤ (9)n2 i=1n

V (θt ) ≡ V [ht (θt )] ≈where sti = ∇ log p(xti |θt ) + N1 ∇ log p(θt ) is
score
∑the
n
of data item i at iteration t and st = n1 i=1 sti is2
the empirical mean. Note that V (θt ) = Nn Vs , where

Vs is the empirical covariance of the scores {sti }, so2scales as Nn . From this we see that the variance of
ϵ2 N 2the stochastic gradient step is t4n M Vs M , so that to
get the injected noise to dominate in all directions, weneed the condition
11
ϵt N 2λmax (M 2 Vs M 2 ) = α ≪ 1
4n

(10)

where λmax (A) is the largest eigenvalue of A. In otherwords, if we choose a stepsize such that the sample
threshold α ≪ 1, the algorithm will be in its Langevindynamics phase and will be sampling approximatelyfrom the posterior.
We can now relate the step size at the sampling threshold to the posterior variance via the Fisher information, which is related to Vs as IF ≈ N Vs , and to theposterior variance Σθ ≈ IF−1 . Using these relationships
as well as (10), we see that the step size at the sampling threshold is ϵt ≈ 4αn
N λmin (Σθ ). Since Langevindynamics explores the posterior via a random walk,using this step size implies that we need on the orderof N/n steps to traverse the posterior, i.e. we process
the whole dataset. So we see this method is not asilver bullet. However, the advantage of the method
is its convenience: stochastic optimization smoothly
and automatically transitions into posterior sampling
without changing the update equation. Even without
measuring the sampling threshold one will enjoy the
beneﬁt of protection against overﬁtting and the ability
to perform Bayesian learning. Measuring the sampling
threshold will only be important if one needs to faithfully represent the posterior distribution with a ﬁnite
collection of samples.
4.2. Estimating Posterior Expectations
Since θ1 , θ2 , . . . converges to the posterior distribution,
we can estimate the posterior expectation E[f (θ)] of
some function f (θ) by simply taking the sample av∑Terage T1 t=1 f (θt ) (as typically in MCMC, we may
remove the initial burn-in phase, say estimated using
the sampling threshold). Since f (θt ) is an asymptotically unbiased estimator for E[f (θ)], this sample average will be consistent. Observe however that because
the step size decreases, the mixing rate of the Markov
chain decreases as well, and the simple sample average will over-emphasize the tail end of the sequence
where there is higher correlation among the samples,
resulting in higher variance in the estimator. Instead
we propose to use the step sizes to weight the samples:∑T
ϵt f (θt )(11)E[f (θ)] ≈ t=1∑T
t=1 ϵt

Stochastic Gradient Langevin Dynamics
3

3

Accuracy on test data

2

1

0.85
0.85

-2−2

Log joint probability per datum

2

00
-1−1

-3−3

1

-4−4

0

0

−1

−1

0

−3

−3
−1

0

1

2

−1

0

1

2

Figure 1. True and estimated posterior distribution.
0

0

10

average rejection rate

noise variance

10

−2

10

−4

∇θ1 noise10∇θ2 noise
0

10

2

4

10
10
iteration

6

10

−2

10

−3

10

−8

10

−6

−4

10

10

−2

10

step size

Figure 2. Left: variances of stochastic gradient noise andinjected noise. Right: rejection probability versus step size.
We report the average rejection probability per iteration ineach sweep through the dataset.

∑∞
Since t=1 ϵt = ∞, this estimator will be consistent
as well. The intuition is that the rate at which theMarkov chain mixes is proportional to the step size, sothat we expect the eﬀective sample size of {θ1 , . . . , θT }∑T
to be proportional to t=1 ϵt , and that each θt willcontribute an eﬀective sample size proportional to ϵt .

5. Experiments
5.1. Simple Demonstration
We ﬁrst demonstrate the workings of our stochastic
gradient Langevin algorithm on a simple example involving only two parameters. To make the posterior
multimodal and a little more interesting, we use a mixture of Gaussians with tied means:θ1 ∼ N (0, σ12 ) ;

22
44
66
88
Number of iterations through whole dataset

10

10

0.65
0.65
00

Accuracy after 10 iterations
Accuracy
0.5
1
1.5
0.5
1.5
Numberof iterations 1
through wholedataset22

Figure 3. Average log joint probability per data item (left)
and accuracy on test set (right) as functions of the number of sweeps through the whole dataset. Red dashed line
represents accuracy after 10 iterations. Results are averaged over 50 runs; blue dotted lines indicate 1 standard
deviation.

−1

10

injected noise

−6

10

0.7
0.7

-6−6

−2

−2

0.75
0.75

-5−5
-7−70

0.8
0.8

θ2 ∼ N (0, σ22 )xi ∼ 12 N (θ1 , σx2 ) + 21 N (θ1 + θ2 , σx2 )where σ12 = 10, σ22 = 1 and σx2 = 2. 100 data points
are drawn from the model with θ1 = 0 and θ2 = 1.
There is a mode at this parameter setting, but also asecondary mode at θ1 = 1, θ2 = −1, with strong negative correlation between the parameters. We ran thestochastic gradient Langevin algorithm with a batch-size of 1 and using 10000 sweeps through the wholedataset. The step sizes are ϵt = a(b + t)−γ whereγ = .55 and a and b are set such that ϵt decreases
from .01 to .0001 over the duration of the run. We see
from Figure 1 that the estimated posterior distribution is very accurate. In Figure 2 we see that there areindeed two phases to the stochastic gradient Langevinalgorithm: a ﬁrst phase where the stochastic gradientnoise dominates the injected noise, and a second phase
where the converse occurs. To explore the scaling of
the rejection rate as a function of step sizes, we reran
the experiment with step sizes exponentially decreasing from 10−2 to 10−8 . In the original experiment the
dynamic range of the step sizes is not wide enough forvisual inspection. Figure 2(right) shows the rejection
probability decreasing to zero as step size decreases.
5.2. Logistic Regression
We applied our stochastic gradient Langevin algorithm
to a Bayesian logistic regression model. The probability of the ith output yi ∈ {−1, +1} given the corresponding input vector xi is modelled as:p(yi |xi ) = σ(yi β ⊤ xi )(12)1where β are the parameters, and σ(z) = 1+exp(−z).
The bias parameter is absorbed into β by including 1
as an entry in xi . We use a Laplace prior for β with ascale of 1. The gradient of the log likelihood is:∂log p(yi |xi ) = σ(−yi β ⊤ xi )yi xi
∂β

(13)

while the gradient of the prior is simply −sign(β),
which is applied elementwise.
We applied our inference algorithm to the a9a dataset
derived by (Lin et al., 2008) from the UCI adultdataset. It consists of 32561 observations and 123 features, and we used batch sizes of 10. Results from 50

Stochastic Gradient Langevin Dynamics

runs are shown in Figure 3, with the model trainedon a random 80% of the dataset and tested on theother 20% in each run. We see that both the joint
probability and the accuracy increase rapidly, with thejoint probability converging after at most 10 iterations,while the accuracy converging after less than 1 iteration through the dataset, demonstrating the eﬃciencyof the stochastic gradient Langevin dynamics.
5.3. Independent Components Analysis
In the following we will brieﬂy review a popular ICA
algorithm based on stochastic (natural) gradient optimization (Amari et al., 1996). We start from a probabilistic model that assumes independent, heavy tailed
marginal distributions,
[
]
∏
∏
T
p(x, W ) = | det(W )|pi (wi x)
N (Wij ; 0, λ)
i

ij

(14)
where we have used a Gaussian prior over the weights.
It has been found that the eﬃciency of gradient descentcan be signiﬁcantly improved if we use a natural gradient. This is implemented by post-multiplication of thegradient with the term W T W (Amari et al., 1996). If
we choose pi (yi ) = 4 cosh12 ( 1 y ) with yi = wiT x, we get
2i.
DW = ∇W log[p(X, W )] W T W =
(
)
N
∑
1
T
NI −
tanh( yn )yn W − λW W T W
2
n=1

This sampler, proposes a new state W ∗ , as in Eqn.16.
Note however that we sum over all data-cases and that
we do not anneal the stepsize. Secondly, we need toaccept or reject the proposed step based on all thedata-cases in order to guarantee detailed balance. The
proposal distribution is given by (suppressing dependence on t),[]
1∗
q(W → W ∗) = N W ; W + εDW ; εM(17)2
where the quadratic function in the exponent is conveniently computed as,−111
tr[(δW − εDW )(W T W )−1 (δW − εDW )T ]
2ε
2
2
(18)
with δW = W ∗ − W and the normalization constantrequires the quantity det M = det(W T W )D . The accept/reject step is then given by the usual Metropolis
Hastings rule:[]p(W ∗ )q(W ∗ → W )p(accept) = min 1,(19)p(W )q(W → W ∗ )Finally, to compute the sampling threshold of Eqn.10,
we can use
1(15)

The term W T W acts like a preconditioning matrix (seesection 4.1), Mij,kl = δik (W T W )jl which is symmetric
under the exchange (i ↔ k, j ↔ l). It can be shown−1
T−1
that the inverse of M is given by√ M =√δ(W W ) ,T
and√ the matrix1 square root as M = δ W W with
W T W = U Λ 2 U T if W T W = U ΛU T .

1

M 2 V(s)M 2 =
(20)
[(
)
]
1
1
T
covn
∇ log p(W ) + ∇ log p(xi |W ) (W W ) 2
Nwith covn the sample covariance over the mini-batchof n data-cases.
To show the utility of our Bayesian variant of ICA wedeﬁne the following “instability” metric for independent components:∑
Ii =var(Wij )var(xj )(21)

The update equation for Langevin dynamics thus becomes,√1
Wt+1 = Wt + εt DW + η t W T W
2(16)

where every element of η t is normally distributed withvariance εt : ηij,t ∼ N [0, εt ]. Our stochastic version
simply approximates the part of the gradient thatsums over data-cases with a sum over a small minibatch of size n and multiplies the result with N/n tobring it back to the correct scale. We also anneal thestepsizes according to εt ∝ a(b + t)−γ .
To assess how well our stochastic Langevin approachcompares against a standard Bayesian method we implemented the ”corrected Langevin” MCMC sampler.

jwhere var(Wij ) is computed over the posterior samples and var(xj ) is computed over the data-cases.
The reason that we scale the variance of the weight
entry Wij with the variance∑ of xj is that the variance of the sources yi = j Wij xj is approximately
equal for all i because they are ﬁt to the distributionpi (yi ) = 4 cosh12 ( 1 y ) .
2i

5.3.1. Artificial Data
In the ﬁrst experiment we generated 1000 data-cases
IID in six channels. Three channels had high kurtosis
distributions while three others where normally distributed. We ran stochastic Langevin dynamics with

Stochastic Gradient Langevin Dynamics
Amari Distance Stoc. Lan.

Amari Distance Corr. Lan.

6

4

2
2

4

6

6

4

2
0

8

iteration

8

0.5

1

1.5

iteration

4

x 10

100

Instability Metric

8

Instability Metric Corr. Lan.

200

Instability Metric

Amari distance

Amari distance

Instability Metric Stoc. Lan.

10

10

150

100

50

0

2

1

4

2

3

4

5

80
60
40
20
0

6

1

Sorted Component ID

x 10

2

3

4

5

6

Sorted Component ID

Figure 4. Left two figures: Amari distance over time for stochastic Langevin dynamics and corrected Langevin dynamics.
Thick line represents the online average. First few hundred iterations were removed to show the scale of the fluctuations.
Right two figures: Instability index for the 6 independent components computed in section 5.3.1 for stochastic Langevindynamics and corrected Langevin dynamics.
PDF W(1,1) vs W(1,2) Stoc. Lan.

PDF W(1,1) vs W(1,2) Corr. Lan.

0
−2

5
4

0
−2

W(2,1)

2

W(2,1)

4

2

W(1,2)

W(1,2)

4

0

W(1,1)

5

2

0

−6
0

3

1

−4

−4
−6
−5

PDF W(1,1) vs W(2,1) Corr. Lan.

PDF W(1,1) vs W(2,1) Stoc. Lan.
5

6

−4

−2

0

2

4

−5
−5

W(1,1)

0

5

−4

−2

0

2

4

W(1,1)

W(1,1)

Figure 5. Posterior density estimates for artificial dataset for stochastic Langevin and corrected Langevin dynamics measured across the W11 − W12 and W1,1 − W2,1 axes.

a batch-size of 100 for a total of 500,000 iterations
and a polynomial annealing schedule εt = N4 t−0.55 .
After around 10,000 iterations the sampling thresholdat α = 0.1 was met. At that point we recorded the“mixing distance” as∑D0 = εt and collected samples
only when the sum t εt from the last sample timeexceeded D0 (in other words, as εt decreases we collect fewer samples per unit time). We note that simply
collecting all samples had no noticeable impact on theﬁnal results. The last estimate of W was used to initialize corrected Langevin dynamics (this was done toforce the samplers into the same local maximum) after
which we also collected 500, 000 samples. For corrected
Langevin we used a constant stepsize of ε = 0.1
N .
The two left ﬁgures of Figure 4 show the Amari distance (Amari et al., 1996) over time for stochastic
and corrected Langevin dynamics respectively. The
right two ﬁgures show the sorted values of our proposed instability index. Figures 5 show two dimensional marginal density estimates of the posterior distribution of W . ICA cannot determine the Gaussiancomponents and this fact is veriﬁed by looking at theposterior distribution. In fact, the stochastic Langevinalgorithm has mixed over a number of modes that presumably correspond to diﬀerent linear combinations of
the Gaussian components. To a lesser degree the corrected Langevin has also explored two modes. Dueto the complicated structure of the posterior distri-bution the stability index varies strongly between thetwo sampling algorithms for the Gaussian components(and in fact also varies across diﬀerent runs). We veriﬁed that the last three components correspond to stable, high kurtosis components.
5.3.2. MEG DataWe

downloaded

the

MEG

datasetfrom

http://www.cis.hut.fi/projects/ica/eegmeg/MEG data.html.

There are 122 channels and 17730 time-points, fromwhich we extracted the ﬁrst 10 channels for our
experiment. To initialize the sampling algorithms,
we ﬁrst ran fastICA (Hyvarinen, 1999) to ﬁnd an
initial estimate of the de-mixing matrix W . We
then ran stochastic Langevin and corrected Langevindynamics to sample from the posterior. The settings
were very similar to the previous experiment with a−0.55
schedule of εt = 0.1
for stochastic Langevin andN ta constant stepsize of 1/N for corrected Langevin.
We obtained 500, 000 samples for stochastic Langevinin 800 seconds and 100,000 samples for corrected
Langevin in 9000 seconds.
We visually veriﬁedthat the two dimensional marginal distributions of
stochastic Langevin and corrected Langevin dynamics
were very similar. The instability values are shown inﬁgure 6. Due to the absence of Gaussian components
we see that the stability indices are very similar acrossthe two sampling algorithms. It was veriﬁed that

Stochastic Gradient Langevin Dynamics
Intability Index Corr. Lan.
0.07

0.06

0.06

Instability Index

Instability Index

Instability Index Stoc. Lan.
0.07

0.05
0.04
0.03
0.02

0.04
0.03
0.02
0.01

0.01
0

0.05

1

2

3

4

5

6

7

8

Sorted Components

9

10

0

1

2

3

4

5

6

7

8

9

10

Sorted Components

Figure 6. Instability indices of 10 components for MEG
dataset for stochastic Langevin (left) and corrected
Langevin (right) respectively.

the most stable component corresponded to a highly
kurtotic source (kurtosus = 15.4), while the most
unstable component was closer to Gaussian noisewith a kurtosis of 3.4 (2 corresponds to Gaussian).
These ﬁndings verify that the stochastic Langevinprocedure produces accurate posterior distributions
that are in full agreement with a well established
MCMC procedure.

6. Discussion
Stochastic gradient optimization is among the most effective algorithms if we measure “predictive accuracy
obtained per unit of computation” (Bottou & Bousquet, 2008). Due to subsampling noise, the parameter estimates ﬂuctuate around their MAP values. The
common wisdom is that one must anneal these stepsizes to zero to reach the ﬁxed point. However, we
argue that one should not optimize beyond the scaleof the posterior distribution. The posterior represents
the intrinsic statistical scale of precision and trying todetermine parameter values with more precision runs
the risk of overﬁtting at additional computational cost.
MCMC sampling from the posterior distribution doesof course address the overﬁtting issue. However, general MCMC algorithms need to see all the data at every iteration, and thus lose the beneﬁts of the stochastic approximation approaches. This paper oﬀers forthe ﬁrst time a surprisingly simple solution that represents the best of both worlds: stick with stochastic
gradients but sample from the posterior nevertheless.
But perhaps the biggest advantage of stochastic gradient Langevin dynamics is the fact that stochastic
optimization seamlessly transitions into posterior sampling. By simply adding Gaussian noise with the correct variance our method performs “early stopping”
automatically without ever having to worry about it.
In fact, we have shown that with a polynomial annealing schedule the obtained samples will asymptotically
represent the posterior distribution faithfully.

We believe that this work represents only a tentativeﬁrst step to further work on eﬃcient MCMC samplingbased on stochastic gradients. Interesting directions ofresearch include stronger theory providing a solid proof
of convergence, deriving a MH rejection step basedon mini-batch data, extending the algorithm to theonline estimation of dynamical systems, and deriving
algorithms based on more sophisticated Hamiltonian
Monte Carlo approaches which do not suﬀer from random walk behaviour.

Acknowledgements
This material is based upon work supported by the National Science Foundation under Grant No. 0447903,
1018433 (MW) and the Gatsby Charitable Foundation(YWT).

References
Amari, S., Cichocki, A., and Yang, H.H. A new algorithm
for blind signal separation. In Neural Information Processing Systems, volume 8, pp. 757–763, 1996.
Andrieu, C., de Freitas, N., and Doucet, A. Sequential
MCMC for Bayesian model selection. In Proceedings ofthe IEEE Signal Processing Workshop on Higher-Order
Statistics, pp. 130–134, 1999.
Bottou, L. and Bousquet, O. The tradeoffs of large scalelearning. In Advances in Neural Information Processing
Systems, volume 20, pp. 161–168, 2008.
Girolami, M. and Calderhead, B. Riemann manifold
Langevin and Hamiltonian Monte Carlo methods. Journal of the Royal Statistical Society B, 73:1–37, 2011.
Hyvarinen, A. Fast and robust fixed-point algorithms forindependent component analysis. IEEE Transactions on
Neural Networks, 10(3):626–634, 1999.
Lin, C.-J., Weng, R. C., and Keerthi, S. S. Trust region
Newton method for large-scale logistic regression. Journal of Machine Learning Research, 9:627–650, 2008.
Neal, R. M. MCMC using Hamiltonian dynamics. In
Brooks, S., Gelman, A., Jones, G., and Meng, X.-L.
(eds.), Handbook of Markov Chain Monte Carlo. Chapman & Hall / CRC Press, 2010.
Robbins, H. and Monro, S. A stochastic approximation
method. Annals of Mathematical Statistics, 22(3):400–
407, 1951.
Robert, C. P. and Casella, G. Monte Carlo statistical methods. Springer Verlag, 2004.
Roberts, G. O. and Stramer, O. Langevin diffusions andmetropolis-hastings algorithms. Methodology and Computing in Applied Probability, 4:337–357, 2002.