Metropolis-adjusted Langevin algorithm explained

In computational statistics, the Metropolis-adjusted Langevin algorithm (MALA) or Langevin Monte Carlo (LMC) is a Markov chain Monte Carlo (MCMC) method for obtaining random samples – sequences of random observations – from a probability distribution for which direct sampling is difficult. As the name suggests, MALA uses a combination of two mechanisms to generate the states of a random walk that has the target probability distribution as an invariant measure:

Informally, the Langevin dynamics drive the random walk towards regions of high probability in the manner of a gradient flow, while the Metropolis–Hastings accept/reject mechanism improves the mixing and convergence properties of this random walk. MALA was originally proposed by Julian Besag in 1994,[1] (although the method Smart Monte Carlo was already introduced in 1978 [2]) and its properties were examined in detail by Gareth Roberts together with Richard Tweedie[3] and Jeff Rosenthal.[4] Many variations and refinements have been introduced since then, e.g. the manifold variant of Girolami and Calderhead (2011).[5] The method is equivalent to using the Hamiltonian Monte Carlo (hybrid Monte Carlo) algorithm with only a single discrete time step.

Further details

Let

\pi

denote a probability density function on

Rd

, one from which it is desired to draw an ensemble of independent and identically distributed samples. We consider the overdamped Langevin Itô diffusion
X

=\nablalog\pi(X)+\sqrt{2}

W

W

. (Note that another commonly-used normalization for this diffusion is
X

=

1
2

\nablalog\pi(X)+

W

,

which generates the same dynamics.) In the limit as

t\toinfty

, this probability distribution

\rho(t)

of

X(t)

approaches a stationary distribution, which is also invariant under the diffusion, which we denote

\rhoinfty

. It turns out that, in fact,

\rhoinfty=\pi

.

Approximate sample paths of the Langevin diffusion can be generated by many discrete-time methods. One of the simplest is the Euler–Maruyama method with a fixed time step

\tau>0

. We set

X0:=x0

and then recursively define an approximation

Xk

to the true solution

X(k\tau)

by

Xk:=Xk+\tau\nablalog\pi(Xk)+\sqrt{2\tau}\xik,

where each

\xik

is an independent draw from a multivariate normal distribution on

Rd

with mean 0 and covariance matrix equal to the

d x d

identity matrix. Note that

Xk

is normally distributed with mean

Xk+\tau\nablalog\pi(Xk)

and covariance equal to

2\tau

times the

d x d

identity matrix.

In contrast to the Euler–Maruyama method for simulating the Langevin diffusion, which always updates

Xk

according to the update rule

Xk:=Xk+\tau\nablalog\pi(Xk)+\sqrt{2\tau}\xik,

MALA incorporates an additional step. We consider the above update rule as defining a proposal

\tilde{X}k

for a new state,

\tilde{X}k:=Xk+\tau\nablalog\pi(Xk)+\sqrt{2\tau}\xik.

This proposal is accepted or rejected according to the Metropolis-Hastings algorithm: set

\alpha:=min\left\{1,

\pi(\tilde{X
k

)q(Xk\mid\tilde{X}k)}{\pi({X}k)q(\tilde{X}k\midXk)}\right\},

where

q(x'\midx)\propto\exp\left(-

1
4\tau

\|x'-x-\tau\nablalog\pi(x)

2
\|
2

\right)

is the transition probability density from

x

to

x'

(note that, in general

q(x'\midx)q(x\midx')

). Let

u

be drawn from the continuous uniform distribution on the interval

[0,1]

. If

u\leq\alpha

, then the proposal is accepted, and we set

Xk:=\tilde{X}k

; otherwise, the proposal is rejected, and we set

Xk:=Xk

.

The combined dynamics of the Langevin diffusion and the Metropolis–Hastings algorithm satisfy the detailed balance conditions necessary for the existence of a unique, invariant, stationary distribution

\rhoinfty=\pi

. Compared to naive Metropolis–Hastings, MALA has the advantage that it usually proposes moves into regions of higher

\pi

probability, which are then more likely to be accepted. On the other hand, when

\pi

is strongly anisotropic (i.e. it varies much more quickly in some directions than others), it is necessary to take

0<\tau\ll1

in order to properly capture the Langevin dynamics; the use of a positive-definite preconditioning matrix

A\inRd

can help to alleviate this problem, by generating proposals according to

\tilde{X}k:=Xk+\tauA\nablalog\pi(Xk)+\sqrt{2\tauA}\xik,

so that

\tilde{X}k

has mean

Xk+\tauA\nablalog\pi(Xk)

and covariance

2\tauA

.

For limited classes of target distributions, the optimal acceptance rate for this algorithm can be shown to be

0.574

; if it is discovered to be substantially different in practice,

\tau

should be modified accordingly.

References

  1. J. Besag. Comments on "Representations of knowledge in complex systems" by U. Grenander and MI Miller. Journal of the Royal Statistical Society, Series B. 1994. 56. 591 - 592.
  2. Rossky . P.J. . Doll . J.D. . Friedman . H.L. . Brownian Dynamics as smart Monte Carlo simulation . Journal of Chemical Physics . 1978 . 69 . 10 . 4628. 10.1063/1.436415 . 1978JChPh..69.4628R .
  3. G. O. Roberts and R. L. Tweedie. Exponential convergence of Langevin distributions and their discrete approximations. Bernoulli. 1996. 2. 4. 341 - 363. 10.2307/3318418. 3318418.
  4. G. O. Roberts and J. S. Rosenthal. Optimal scaling of discrete approximations to Langevin diffusions. Journal of the Royal Statistical Society, Series B. 1998. 60. 1. 255 - 268. 10.1111/1467-9868.00123. 5831882.
  5. M. Girolami and B. Calderhead. Riemann manifold Langevin and Hamiltonian Monte Carlo methods. Journal of the Royal Statistical Society, Series B. 2011. 73. 2. 123 - 214. 10.1111/j.1467-9868.2010.00765.x. 10.1.1.190.580.