Metropolis-adjusted Langevin algorithm

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

${\displaystyle {\dot {X}}=\nabla \log \pi (X)+{\sqrt {2}}{\dot {W}}}$

driven by the time derivative of a standard Brownian motion ${\displaystyle W}$. (Note that another commonly-used normalization for this diffusion is

${\displaystyle {\dot {X}}={\frac {1}{2}}\nabla \log \pi (X)+{\dot {W}},}$

which generates the same dynamics.) In the limit as ${\displaystyle t\to \infty }$, this probability distribution ${\displaystyle \rho (t)}$ of ${\displaystyle X(t)}$ approaches a stationary distribution, which is also invariant under the diffusion, which we denote ${\displaystyle \rho _{\infty }}$. It turns out that, in fact, ${\displaystyle \rho _{\infty }=\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 ${\displaystyle \tau >0}$. We set ${\displaystyle X_{0}:=x_{0}}$ and then recursively define an approximation ${\displaystyle X_{k}}$ to the true solution ${\displaystyle X(k\tau )}$ by

${\displaystyle X_{k+1}:=X_{k}+\tau \nabla \log \pi (X_{k})+{\sqrt {2\tau }}\xi _{k},}$

where each ${\displaystyle \xi _{k}}$ is an independent draw from a multivariate normal distribution on ${\displaystyle \mathbb {R} ^{d}}$ with mean 0 and covariance matrix equal to the ${\displaystyle d\times d}$ identity matrix. Note that ${\displaystyle X_{k+1}}$ is normally distributed with mean ${\displaystyle X_{k}+\tau \nabla \log \pi (X_{k})}$ and covariance equal to ${\displaystyle 2\tau }$ times the ${\displaystyle d\times d}$ identity matrix.

In contrast to the Euler–Maruyama method for simulating the Langevin diffusion, which always updates ${\displaystyle X_{k}}$ according to the update rule

${\displaystyle X_{k+1}:=X_{k}+\tau \nabla \log \pi (X_{k})+{\sqrt {2\tau }}\xi _{k},}$

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

${\displaystyle {\tilde {X}}_{k+1}:=X_{k}+\tau \nabla \log \pi (X_{k})+{\sqrt {2\tau }}\xi _{k}.}$

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

${\displaystyle \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\},}

where

${\displaystyle q(x'\mid x)\propto \exp \left(-{\frac {1}{4\tau }}\|x'-x-\tau \nabla \log \pi (x)\|_{2}^{2}\right)}$

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

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

${\displaystyle {\tilde {X}}_{k+1}:=X_{k}+\tau A\nabla \log \pi (X_{k})+{\sqrt {2\tau A}}\xi _{k},}$

so that ${\displaystyle {\tilde {X}}_{k+1}}$ has mean ${\displaystyle X_{k}+\tau A\nabla \log \pi (X_{k})}$ and covariance ${\displaystyle 2\tau A}$.

For limited classes of target distributions, the optimal acceptance rate for this algorithm can be shown to be ${\displaystyle 0.574}$; if it is discovered to be substantially different in practice, ${\displaystyle \tau }$ should be modified accordingly.[3]