Goals

• Brief reminder of Metropolis–Hastings

• Review some popular choices of candidate kernel

• Some known results (I will try to give a range of different approaches)

• Again, if time, I’ll mention some of my own contributions/interests

Metropolis–Hastings recap

• General `accept-reject’ \(\pi\)-reversible kernel \[ P(x,dx') = \alpha(x,x')Q(x,dx') + \left( \int (1-\alpha(x,x'))Q(x,dx') \right) \delta_x(dx') \]

• Last time we showed that you should set \[ \alpha(x,x') = \min\left(1, \frac{\pi(dx')Q(x',dx)}{\pi(dx)Q(x,dx')} \right) \]

• This time we will discuss \(Q\), and some connections (e.g. beyond MH)

Two standard approaches to choosing \(Q\)

• Choose a very simple \(Q\)
• Choose a \(Q\) based on some \(\pi\)-invariant stochastic process

• In both cases there are typically tuning parameters in \(Q\), the choice of which is an active area of study and crucial for performance

Random Walk Metropolis: a very simple \(Q\)

• Set \(Q(x,dx') := q(x'-x)m(dx')\), \(q\) symmetric (e.g. \(N(x, h^2I)\))

• With this choice \(q(x-x')/q(x'-x) = 1\), meaning \[ \alpha(x,x') = \min\left( 1, \frac{\pi(dx')}{\pi(dx)} \right) \]

• Uphill moves are accepted

Tuning the Random Walk Metropolis

• Typically \(Q\) is chosen to be \(N(x,h^2I)\)
  • \(h\) small \(\implies\) all moves are accepted, but exploration is slow
  • \(h\) large \(\implies\) all moves are rejected
  • How should we choose \(h\)?

• Roberts, Gelman & Gilks (1997) proposed an approach based on diffusion limits

• It’s now very established and similar approaches have been used many times

Diffusion limits for Random Walk Metropolis

• Take continuous-time jump process with generator \[ L_df := \lambda_d (P_d - I) f, \]

• \(P_d\) is MH kernel using \(Q_d\). Choose \(Q_d\) as \(N(x, h_d^2 I)\)
  • Set \(h_d = \ell \cdot d^{-\gamma}\) for \(\ell,\gamma \in (0,\infty)\)
  • Set \(\lambda_d = d^{\gamma}\)

• If \(\gamma\) and is chosen correctly then for a suitable class \(f\) (e.g. restricted to \(x_1\)) and some forms of \(\pi\) \[ L_d f(x) \to h(\ell)\mathcal{L}f(x) \] uniformly in \(x\), where \(\mathcal{L}f\) is the generator of a \(\pi\)-invariant diffusion

• Essentially this means:
  • The chain takes \(O(d^\gamma)\) steps to explore state space (\(\gamma = 1\) for random walk)
  • Maximising \(h(\ell) \implies\) acceptance rate guidelines (e.g. \(\mathbb{E}[\alpha(x,x')] = 0.234\) for RWM)

Diffusion limits picture

Overdamped Langevin diffusion

• Consider the diffusion process (with generator \(\mathcal{L}\)) \[ dX_t = \frac{1}{2}\nabla \log\pi(X_t)dt + dW_t \]

• Using the Fokker–Planck equation \[ \mathcal{L}^* \pi = 0 \implies \pi \text{ is invariant.} \]

• Satisfying this here is as simple as noting that for all \(i\) \[ \left[ \frac{\partial}{\partial x_i}\log\pi(x)\right]\pi(x) = \frac{\partial}{\partial x_i} \pi(x) \]

Using Langevin processes in MCMC

• We cannot exactly simulate the diffusion in general

• Instead we can use \[ Y = X_i + \frac{h^2}{2}\nabla\log\pi(X_i) + h\xi_{i+1}, \quad \xi_{i+1} \sim N(0,I). \] (Euler–Maruyama discretisation).

• Can use as a proposal within MH creating MALA (Metropolis-adjusted Langevin algorithm)

• Some prefer to use it without correction and show that under suitable assumptions the error is not too large, e.g. \[ \mathcal{D}(\pi_h,\pi) \leq \epsilon(h) \] for a suitable metric \(\mathcal{D}\), where \(\pi_h\) is the equilibrium law of the discretised process

Contraction for Langevin diffusions

• Write \(\pi(dx) \propto e^{-U(x)}dx\). We show exponential convergence in Wasserstein-2 distance if \(U\) is strongly convex \[ \mathcal{W}_2(\mu,\nu) := \inf_{\xi \in \mathcal{C}(\mu,\nu)} \left( \mathbb{E}_\xi|X-Y|^2 \right)^{1/2}, \] where \(\mathcal{C}(\mu,\nu):= \{ \xi \in \mathcal{P}(\mathcal{X} \times \mathcal{Y}) : \xi(A \times Y) = \mu(A), ~ \xi(X \times B) = \nu(B) \}\)

• Define two Langevin diffusions \((X_t)_{t\geq 0}\) and \((Y_t)_{t\geq 0}\) s.t.
  • \(X_0 \sim \mu_0\) and \(Y_0 \sim \nu_0\), with \(\xi_0 := \inf \mathcal{C}(\mu_0,\nu_0)\)
  • \((X_t)_{t\geq 0}\) and \((Y_t)_{t \geq 0}\) have the same Brownian part: synchronous coupling

• This implies \[ \frac{d}{dt} \|X_t - Y_t \|^2 = 2\langle X_t - Y_t, \frac{d}{dt}(X_t - Y_t) \rangle = -\langle X_t - Y_t, \nabla U(X_t) - \nabla U(Y_t) \rangle \]

Contraction for Langevin diffusions

• Now, assuming \(U\) is strongly convex \(\exists m>0\) s.t. \[ \frac{d}{dt} \|X_t - Y_t \|^2 = - \langle X_t - Y_t, \nabla U(X_t) - \nabla U(Y_t) \rangle \leq -m\|X_t - Y_t\|^2 \]

• Writing \(f(t) := \|X_t - Y_t\|^2\) then \(f'(t) \leq -mf(t) \implies f(t) \leq e^{-mt}f(0)\) (Gronwall). Taking expectations \[ \mathbb{E}\|X_t - Y_t \|^2 \leq e^{-mt} \mathbb{E}\|X_0 - Y_0 \|^2 \]

• Writing \(X_t \sim \mu_t\) and \(Y_t \sim \nu_t\) and noting \(W_2(\mu_t,\nu_t)^2 \leq \mathbb{E}_\xi\|X_t - Y_t\|^2\) gives \[ W_2(\mu_t, \nu_t) \leq e^{-\frac{m}{2}t}W_2(\mu_0,\nu_0). \]

• If we set \(\nu_0 = \pi\), then by \(\pi\)-invariance \(\nu_t = \pi\) \(\forall t\), giving the convergence bound \(W_2(\mu_t,\pi) \leq e^{-mt/2}W_2(\mu_0,\pi)\)

• See e.g. Bolley, Gentil & Guillin (2012)

There are many more things known about Langevin MCMC…

• Using diffusion limits it is known that:
  • The MALA chain takes \(O(d^{1/3})\) steps to explore state space
  • Optimal \(\mathbb{E}\alpha \approx 0.574\)

• General conditions under which the diffusion is non-explosive and converges to equilibrium at an exponential rate are given in Roberts & Tweedie (1996) (very broad)

• Lots of active work on studying these processes as they are very regular

Hamiltonian Monte Carlo

• The state-of-the-art sampler

• First introduced by physicists in 1987

• Introduced to Statistics by Radford Neal and others (since early 90s)

• Not popular until influential review paper (Neal, 2011) and Stan software

HMC idea

• Call \(x\) position, introduce momentum \(p\), and measure \[ \mu(dx,dp) = \pi(dx)\nu(dp) \propto e^{-H(x,p)}dxdp \]

• Use \(\mu\)-invariant deterministic dynamics (fast exploration) combined with samples from \(\nu\) (ensures ergodicity)

• Deterministic dynamics are Hamiltonian \(H = -\log\pi(x) + K(p)\) \[ \begin{aligned} \dot{x}_t &= \nabla K(p_t) \\ \dot{p}_t &= \nabla \log\pi(x_t) \end{aligned} \]
• If \((X_0,P_0) \sim \mu\), then \(\varphi_t(X_0, P_0) \sim \mu\) for all \(t\), where \(\varphi_t\) solves the above system

• Choice of kinetic energy \(K(p)\) left to the user (more later)

Hamiltonian Monte Carlo demo

Show simulation

HMC as a continuous-time process

• PDMP: Write \(z = (x,p)\) and \(J := \left( \begin{array}{cc} 0 & I \\ -I & 0 \end{array} \right)\), and set \[ Lf(z) = \langle J\nabla H(z), \nabla f(z) \rangle + \lambda \int [f(z') - f(z)] \delta_x(dx')\nu(dp') \]

• Diffusion: Hamiltonian dynamics + Ornstein–Uhlenbeck process on \(p\) \[ \begin{aligned} dX_t &= P_tdt \\ dP_t &= \nabla \log\pi(X_t)dt -\gamma P_tdt + dW_t \end{aligned} \]

• The processes are different, but in my experience convergence behaviour is similar…

My life in HMC: geometric ergodicity

• There are many more details to implementation that I won’t mention (e.g. numerical integrators, MH correction)

• Early on I worked on convergence to equilibrium for HMC, precisely when this will be geometric, meaning \[ \mathcal{D}(\delta_x P^n, \pi) \leq CV(x)\rho^n, \] for \(\rho < 1\), \(C < \infty\), \(V:\mathcal{X} \to [1,\infty)\), where \[ \mathcal{D}(\mu,\nu) := \sup_{|f| \leq V} |\mu(f) - \nu(f)| \]

• Such a bound induces CLT for ergodic averages of \(f \in L^2(\pi)\) (if \(P\) is \(\pi\)-reversible)

• Proof follows drift and minorisation approach (e.g. Meyn & Tweedie). Also called Harris’ theorem.

• A very neat and direct general proof of Harris’ theorem is given in Hairer & Mattingley (2011)

My life in HMC: kinetic energy

• I also worked on the best choice of \(K\). Typically \[ K(p) = \frac{1}{2}p^T M^{-1} p, \] meaning \(\nu(dp) \propto e^{-K(p)}dp\) is Gaussian.

• Using some theory and some heuristics, writing \(U = -\log\pi\) (potential), we recommended that \[ \nabla K = (\nabla U)^{-1} \] when \(\|p\|\) is large, but that \(K\) behaves quadratically for \(\|p\|\) small.

• This implies \(K\) should be a regularised Legendre transform of \(U\). Connections to Codina’s talks…

The future

• There is now interest in general piece-wise deterministic processes (these are not reversible but \((\mu,Q)\)-self-adjoint, from last week)

• I am also working with collaborators on a class of jump processes that can be used as proposals, with nice properties

• What is really missing are lower bounds (like Nesterov results in optimisation). Until these arrive we can keep tinkering…

References

• Bolley, F., Gentil, I., & Guillin, A. (2012). Convergence to equilibrium in Wasserstein distance for Fokker–Planck equations. Journal of Functional Analysis, 263(8), 2430-2457.
• Hairer, M., & Mattingly, J. C. (2011). Yet another look at Harris’ ergodic theorem for Markov chains. In Seminar on Stochastic Analysis, Random Fields and Applications VI (pp. 109-117). Springer, Basel.
• Livingstone, S., Betancourt, M., Byrne, S., & Girolami, M. (2019). On the geometric ergodicity of Hamiltonian Monte Carlo. Bernoulli, 25(4A), 3109-3138.
• Livingstone, S., Faulkner, M. F., & Roberts, G. O. (2019). Kinetic energy choice in Hamiltonian/hybrid Monte Carlo. Biometrika, 106(2), 303-319.
• Meyn, S. P., & Tweedie, R. L. (2012). Markov chains and stochastic stability. Springer Science & Business Media.
• Roberts, G. O., & Rosenthal, J. S. (2001). Optimal scaling for various Metropolis-Hastings algorithms. Statistical science, 16(4), 351-367.
• Roberts, G. O., & Tweedie, R. L. (1996). Exponential convergence of Langevin distributions and their discrete approximations. Bernoulli, 2(4), 341-363.