Reconciling Langevin MC methods as one-step HMC versus as diffusion or brownian motion

Question

I have a basic understanding of Hamiltonian monte carlo and why it works. I've read that Langevin MC is basically a special case of HMC when you only step the dynamics forward a single timestep before resampling the momentum -- and I see how this results in the "simplified" equations where you don't even need to explicitly write down the momentum -- you can just update by adding random noise to the gradient of the log unnormalized posterior.

I've also read about the Langevin process as being some sort of diffusion or brownian motion model. While this makes sense to me on an intuitive level, and I see how one might derive a Langevin update from that physical process, I don't understand how it fits into the framework of MCMC.

In particular: I trust HMC because one can prove (ignoring discretization error) the detailed balance criterion. By extension, LMC, framed as a special case of HMC, can also be trusted. However I don't understand how the concept of brownian motion or diffusion can be used to prove that LMC converges to the posterior.

To summarize: is the idea of Langevin dynamics as brownian motion / diffusion merely an intuition pump, or can it be mathematically formalized to prove that LMC converges to the posterior? And if so, how?

Answer 1

The easiest way to understand why Langevin dynamics targets the "correct distribution" is to look at the corresponding Fokker-Planck equation.

Let me be more precise. Let us assume that our target distribution has the following density:

$\pi(x) = \frac1{Z} \exp(-U(x))$,

where $x \in \mathbb{R}^d$, $U$ is often called the potential energy, and $Z$ is the normalizing constant.

The Langevin algorithm, either the Metropolis-adjusted version (MALA) or the unadjusted version (ULA), is based on the Langevin diffusion, that is described by the following stochastic differential equation (SDE):

$dX_t = -\nabla U(X_t)dt + \sqrt{2}dB_t$,

where $B_t$ denotes the standard Brownian motion. To understand how the probability density function of $X_t$ evolves in time, we can use a very useful tool, which is called the Fokker-Planck equation (FPE). For notational simplicity let us assume we are in the scalar case, i.e. $d=1$ (for $d>1$ the idea is the same). In this case the FPE reads:

$\partial_t p(x,t) = \partial_x (\partial_x U(x) p(x,t)) + \partial_x^2 p(x,t)$,

where $p(x,t)$ is the probability density function of $X_t$ at time $t$. Now the main trick is this:

Let us assume the process $(X_t)_{t \geq 0}$ (which is a Markov process) is ergodic with its invariant measure. Then, when $p(x,t)$ reaches to this invariant measure, it cannot deviate from it anymore (since it's the invariant measure). Hence, when $p(x,t)$ converges to the invariant measure, then $\partial_t p(x,t)$ must be equal to zero (since it won't change over the time $t$).

Given this observation, in order to verify that the invariant measure of the Langevin equation is indeed $\pi$, we only need to check if $\partial_t p(x,t)$ becomes $0$ or not, when we replace $p(x,t)$ by $\pi(x)$. Now, let's see if this really happens:

\begin{align} \partial_t p(x,t) &= \partial_x (\partial_x U(x) \pi(x)) + \partial_x^2 \pi(x) \\ &= \partial_x (\partial_x U(x) \pi(x) + \partial_x \pi(x)) \end{align} By using the fact that $\partial_x U(x) = -\partial_x \log \pi(x)= - \frac1{\pi(x)} \partial_x \pi(x)$, we obtain: \begin{align} \partial_t p(x,t) &= \partial_x (- \frac1{\pi(x)} \pi(x) \partial_x \pi(x) + \partial_x \pi(x)) \\ &= \partial_x (- \partial_x \pi(x) + \partial_x \pi(x)) \\ &= 0. \end{align}

Then, this shows that $\pi$ is an invariant distribution of the process $(X_t)_t$, and if $(X_t)_t$ has a unique invariant distribution (for instance if $\nabla U$ is Lipschitz), then $\pi$ is the unique invariant distribution.

Now, if we go back to ULA (unadjusted Langevin algorithm) or MALA (Metropolis adjusted Langevin algorithm), they are both based on the Euler discretization of the Langevin SDE:

$Y_{n+1} = Y_{n} - \eta \nabla U(Y_n) + \sqrt{2\eta} Z_{n+1}$,

where $Z_n$ is a standard Gaussian random variable. This scheme is called ULA. The process $(Y_n)_n$ is still a Markov process, but it doesn't target $\pi$ anymore due to the discretization error. Hence, one can couple it with a Metropolis acceptance step to get rid of this error. The resulting algorithm is called MALA.

I hope it's more clear now.

Answer 2

Here is an elementary and not entirely rigorous "proof" that Langevin dynamics satisfies detailed balance. I only bothered considering the univariate case.

Let $\pi$ denote the distribution we are concerned with (and the distribution for which we'll prove detailed balance).

A single step of LMC is: \begin{align} x \leftarrow x + \epsilon \eta - \frac{\epsilon^2}{2} \nabla \log \pi(x) \quad \text{where $\eta \sim \mathcal{N}(0, 1)$} \end{align}

Consider any two points $x$ and $x'$, with $\delta x = (x'-x)/\epsilon$. If we assume $\log \pi(x)$ is locally linear, then \begin{align} \nabla \log \pi(x) = \nabla \log \pi(x') = \frac{\log \pi(x') - \log \pi(x)}{\epsilon \delta x} \end{align} With that in mind, we can write: \begin{align} P(x \rightarrow x') &= P(x' -x = \epsilon \eta - \frac{\epsilon^2}{2} \nabla \log \pi(x)) \\ &= P(\epsilon \delta x = \epsilon \eta - \frac{\epsilon}{2 \delta x} (\log \pi(x') - \log \pi(x))) \\ &= P(\eta = \delta x + \frac{1}{2 \delta x} (\log \pi(x') - \log \pi(x))) \end{align} For conciseness, write $\phi = \frac{1}{2} (\log \pi(x') - \log \pi(x))$

\begin{align} \log P(x \rightarrow x') &= -\frac{1}{2} (\delta x + \frac{\phi}{\delta x})^2 \\ &= -\frac{1}{2} ((\delta x)^2 + 2 \phi + (\frac{\phi}{\delta x})^2) \end{align} And analogously (although I don't bother to write out the steps):

\begin{align} \log P(x' \rightarrow x) &= -\frac{1}{2} ((\delta x)^2 - 2 \phi + (\frac{\phi}{\delta x})^2) \end{align} Thus the difference is

\begin{align} \log P(x' \rightarrow x) - \log P(x \rightarrow x') &= -2 \phi \end{align} Therefore we can complete the proof with: \begin{align} \log \pi(x) - \log \pi(x') &= \log P(x' \rightarrow x) - \log P(x \rightarrow x') \\ \log \pi(x) + \log P(x \rightarrow x') &= \log \pi(x') + \log P(x' \rightarrow x) \end{align}

---

Edit: Actually my assumption about $\log \pi$ being roughly linear doesn't seem to be reasonable, unless I can guarantee that $x$ and $x'$ are close enough to each other. I'm sure there's a way to patch the proof, but I'm not quite sure how exactly.