View source for Metropolis-adjusted Langevin algorithm

{{Short description|Markov Chain Monte Carlo algorithm}}
In [[computational statistics]], the '''Metropolis-adjusted Langevin algorithm (MALA)''' or '''Langevin Monte Carlo (LMC)''' is a [[Markov chain Monte Carlo]] (MCMC) method for obtaining [[pseudo-random number sampling|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]]:
* new states are proposed using ([[overdamped]]) '''[[Langevin dynamics]]''', which use evaluations of the [[gradient]] of the target [[probability density function]];
* these proposals are accepted or rejected using the '''[[Metropolis–Hastings algorithm]]''', which uses evaluations of the target probability density (but not its gradient).
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,<ref name="Besag1994">{{cite journal
| author = J. Besag
| title = Comments on "Representations of knowledge in complex systems" by U. Grenander and MI Miller
| journal = Journal of the Royal Statistical Society, Series B
| year = 1994
| volume = 56
| pages = 591&ndash;592
}}</ref> (although the method [[Smart Monte Carlo]] was already introduced in 1978 <ref>{{cite journal |last1=Rossky |first1=P.J. |last2=Doll |first2=J.D. |last3=Friedman |first3=H.L. |title=Brownian Dynamics as smart Monte Carlo simulation |journal=Journal of Chemical Physics |date=1978 |volume=69 |issue=10 |page=4628|doi=10.1063/1.436415 |bibcode=1978JChPh..69.4628R }}</ref>) and its properties were examined in detail by  [[Gareth Roberts (statistician)|Gareth Roberts]] together with [[Richard Tweedie]]<ref name="RobertsTweedie1996">{{cite journal
| author = G. O. Roberts and R. L. Tweedie
| title = Exponential convergence of Langevin distributions and their discrete approximations
| journal = Bernoulli
| year = 1996
| volume = 2
| number = 4
| pages = 341&ndash;363
| doi = 10.2307/3318418
| jstor = 3318418
| url = http://projecteuclid.org/euclid.bj/1178291835
}}</ref> and [[Jeff Rosenthal]].<ref name="RobertsRosenthal1998">{{cite journal
| author = G. O. Roberts and J. S. Rosenthal
| title = Optimal scaling of discrete approximations to Langevin diffusions
| journal = Journal of the Royal Statistical Society, Series B
| year = 1998
| volume = 60
| number = 1
| pages = 255&ndash;268
| doi = 10.1111/1467-9868.00123
| s2cid = 5831882
}}</ref>  Many variations and refinements have been introduced since then, e.g. the [[manifold]] variant of Girolami and Calderhead (2011).<ref name="GirolamiCalderhead2011">{{cite journal
| author = M. Girolami and B. Calderhead
| title = Riemann manifold Langevin and Hamiltonian Monte Carlo methods
| journal = Journal of the Royal Statistical Society, Series B
| year = 2011
| volume = 73
| number = 2
| pages = 123&ndash;214
| doi = 10.1111/j.1467-9868.2010.00765.x
| citeseerx = 10.1.1.190.580
}}</ref> The method is equivalent to using the [[Hamiltonian Monte Carlo]] (hybrid Monte Carlo) algorithm with only a single discrete time step.<ref name="GirolamiCalderhead2011" />

==Further details==

Let <math>\pi</math> denote a probability density function on <math>\mathbb{R}^{d}</math>, one from which it is desired to draw an ensemble of [[independent and identically distributed]] samples.  We consider the overdamped Langevin  [[Itô diffusion]]

: <math>\dot{X} = \nabla \log \pi(X) + \sqrt{2} \dot{W}</math>

driven by the time derivative of a standard [[Brownian motion]] <math>W</math>.  (Note that another commonly-used normalization for this diffusion is

: <math>\dot{X} = \frac{1}{2} \nabla \log \pi(X) + \dot{W},</math>

which generates the same dynamics.)  In the limit as <math>t \to \infty</math>, this probability distribution <math>\rho(t)</math> of <math>X(t)</math> approaches a stationary distribution, which is also invariant under the diffusion, which we denote <math>\rho_\infty</math>.  It turns out that, in fact, <math>\rho_\infty = \pi</math>.

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 <math>\tau > 0</math>.  We set <math>X_0 := x_0</math> and then recursively define an approximation <math>X_k</math> to the true solution <math>X(k \tau)</math> by

:<math>X_{k + 1} := X_k + \tau \nabla \log \pi(X_k) + \sqrt{2 \tau} \xi_k,</math>

where each <math>\xi_{k}</math> is an independent draw from a [[multivariate normal distribution]] on <math>\mathbb{R}^{d}</math> with [[expected value|mean]] 0 and [[covariance matrix]] equal to the <math>d \times d</math> [[identity matrix]].  Note that <math>X_{k + 1}</math> is normally distributed with mean <math>X_k + \tau \nabla \log \pi(X_k)</math> and covariance equal to <math>2 \tau</math> times the <math>d \times d</math> identity matrix.

In contrast to the Euler–Maruyama method for simulating the Langevin diffusion, which always updates <math>X_k</math> according to the update rule

:<math>X_{k + 1} := X_k + \tau \nabla \log \pi(X_k) + \sqrt{2 \tau} \xi_k,</math>

MALA incorporates an additional step.  We consider the above update rule as defining a ''proposal'' <math>\tilde{X}_{k + 1}</math> for a new state,

:<math>\tilde{X}_{k + 1} := X_k + \tau \nabla \log \pi(X_k) + \sqrt{2 \tau} \xi_k.</math>

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

:<math>\alpha := \min \left\{ 1 , \frac{\pi(\tilde{X}_{k + 1}) q(X_{k}\mid\tilde{X}_{k + 1})}{\pi({X}_{k}) q(\tilde{X}_{k + 1}\mid X_k)} \right\},</math>

where

:<math>q(x'\mid x) \propto \exp \left( - \frac{1}{4 \tau} \| x' - x - \tau \nabla \log \pi(x) \|_2^2 \right)</math>

is the transition probability density from <math>x</math> to <math>x'</math> (note that, in general <math>q(x'\mid x) \neq q(x\mid x')</math>).  Let <math>u</math> be drawn from the [[continuous uniform distribution]] on the interval <math>[0, 1]</math>.  If <math>u \leq \alpha</math>, then the proposal is accepted, and we set <math>X_{k + 1} := \tilde{X}_{k + 1}</math>;  otherwise, the proposal is rejected, and we set <math>X_{k + 1} := X_k</math>.

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 <math>\rho_{\infty} = \pi</math>.  Compared to naive Metropolis–Hastings, MALA has the advantage that it usually proposes moves into regions of higher <math>\pi</math> probability, which are then more likely to be accepted.  On the other hand, when <math>\pi</math> is strongly [[isotropy|anisotropic]] (i.e. it varies much more quickly in some directions than others), it is necessary to take <math>0 < \tau \ll 1</math> in order to properly capture the Langevin dynamics;  the use of a positive-definite [[preconditioning]] matrix <math>A \in \mathbb{R}^{d \times d}</math> can help to alleviate this problem, by generating proposals according to
:<math>\tilde{X}_{k + 1} := X_k + \tau A \nabla \log \pi(X_k) + \sqrt{2 \tau A} \xi_k,</math>
so that <math>\tilde{X}_{k + 1}</math> has mean <math>X_k + \tau A \nabla \log \pi(X_k)</math> and covariance <math>2 \tau A</math>.

For limited classes of target distributions, the optimal acceptance rate for this algorithm can be shown to be <math>0.574</math>; if it is discovered to be substantially different in practice, <math>\tau</math> should be modified accordingly.<ref name="RobertsRosenthal1998" />

==References==
<references />

[[Category:Monte Carlo methods]]
[[Category:Markov chain Monte Carlo]]
[[Category:Sampling techniques]]