Langevin diffusions and the Metropolis-adjusted Langevin algorithm

Introduction

The Metropolis-adjusted Langevin algorithm (MALA) (e.g. Roberts and Rosenthal, 1998) and its manifold variant (MMALA) Girolami and Calderhead (2011) are Markov chain Monte Carlo methods based on diffusions. While theoretical properties of the former are better understood (e.g. Roberts and Rosenthal, 1998), the latter has been shown to be more effective in practice, producing more efficient estimates for the same computational budget in many experiments Girolami and Calderhead (2011). In this article we highlight two properties of the diffusion on which MMALA is based. First, we point out an unfortunate transcription error which has propagated through the literature, whereby a factor of a 1/2 has been missed from one of the terms Roberts and Stramer (2002); Girolami and Calderhead (2011). Second, we show that the corrected diffusion does not have the intended invariant density with respect to Lebesgue measure. It would seem logical that a similar diffusion which does preserve the intended probability density may prove a better basis for a Metropolis–Hastings algorithm. We therefore describe such a diffusion and the resulting sampling method, which we call PMALA (position-dependent MALA). We show that the incorrectly transcribed diffusion and that on which PMALA is based are equivalent in some cases, although the former leads to a more computationally costly algorithm; this equivalence explains to some extent why the error has been missed previously. Finally we describe simulation studies based on those in Girolami and Calderhead (2011) comparing PMALA and MMALA. In terms of effective sample size (ESS) PMALA outperforms MMALA when the two are not equivalent. PMALA always outperforms MMALA in terms of effective sample size (ESS) per second, since even when the two algorithms are equivalent, each step of PMALA involves fewer CPU operations.

Langevin diffusions

The law of the diffusion is described by the Fokker–Planck equation (e.g. Oksendal, 1998), which relates the evolution of the probability density function $u(\bm{x},t)$ for $\bm{X}_{t}$ to the drift and volatility \bm{b},\mbox{\boldmath\sigma},

where V(\bm{x})=\mbox{\boldmath\sigma}(\bm{x})\mbox{\boldmath\sigma}(\bm{x})^{T}. If $u(\bm{x},t)=\pi(\bm{x})$ for all $t$, then the process is stationary, and $\pi$ is the density of the invariant or stationary distribution of the diffusion, meaning that if $\bm{X}_{t}\sim\pi(\cdot)$ then $\bm{X}_{t+\tau}\sim\pi(\cdot)$ for all $\tau>0$ (e.g. Oksendal, 1998). One such is the Langevin diffusion, the solution to:

The Metropolis–Hastings algorithm simulates from a Markov chain which has a desired invariant density, $\pi(\bm{x})$. Expectations from this distribution can be approximated by averaging values across the chain (e.g. Gilks et al., 1996). At each iteration some proposal $\bm{x}^{\prime}$ is drawn from a distribution with density $q(\cdot|\bm{x})$ (where $\bm{x}$ represents the current value in the chain). The next value in the chain is set to be $\bm{x}^{\prime}$ with probability $\alpha(\bm{x}^{\prime}|\bm{x})$, or else $\bm{x}$, where:

for a chosen step size $h$, and then accepted with probability $\alpha(\bm{x}^{\prime}|\bm{x})$. Scaling properties of $h$ with $d$ and asymptotic optimal acceptance rates for the method are discussed in Roberts and Rosenthal (1998). A slight generalisation of (2) is the diffusion:

where $A$ is a positive-definite matrix, and $\surd{A}\surd{A}^{T}=A$. As with the diffusion (2), substitution of the drift and volatility terms from (4) into the Fokker–Planck equation leads to $\partial u/\partial t=0$, so that $\pi$ is the invariant density of (4). The Metropolis–Hastings scheme derived from (4) is known as ‘pre-conditioned MALA’ Roberts and Stramer (2002) and is well-suited to scenarios in which the components of $\pi$ are highly correlated or have very different marginal variances, but where these relationships vary little over the main posterior mass.

The MMALA algorithm Girolami and Calderhead (2011) is based on the discretisation of a diffusion with a position-dependent volatility matrix:

where $G(\bm{X}_{t})$ is some positive definite $d\times d$ matrix. The choice of $G$ is arbitrary, but some natural candidates arise by noting that the above process can be thought of as a diffusion defined on a Riemannian manifold, specified in local coordinates (Girolami and Calderhead (2011)). In the resulting algorithm, proposals are generated according to:

and then accepted or rejected according to (3). A similar scheme is proposed in Roberts and Stramer (2002), based on the same diffusion. For a suitable choice of $G(\bm{x})$, the position-dependent covariance matrix for proposals in (6) allows adaptation to the local curvature of the target density $\pi$, which has been shown to increase algorithm efficiency in a number of examples Girolami and Calderhead (2011).

A new position-dependent diffusion and MALA

In general, a diffusion with invariant density $\pi$ can be constructed by starting from (1) and selecting a drift and volatility such that

If the intention is to derive a Metropolis–Hastings proposal mechanism with a position-dependent covariance matrix, a natural starting point would be to simply set $A=A(\bm{X}_{t})$ in (4), giving $d\bm{X}_{t}=\frac{1}{2}A(\bm{X}_{t})\nabla\log\pi(\bm{X}_{t})dt+\surd{A(\bm{X}_{t})}~{}d\bm{W}_{t}$, and this diffusion forms the basis of the simplified MMALA algorithm of Girolami and Calderhead (2011). However, substituting the drift and volatility terms into (7) gives the requirement that:

for each $i$. Since $\partial/\partial x_{j}[\log\pi(\bm{x})]\pi(\bm{x})=\partial\pi(\bm{x})/\partial x_{j}$, (8) is only satisfied in general when $A$ is a constant matrix. A simple modification to the drift term, however, leads to a new diffusion which satisfies (7):

This diffusion has invariant density $\pi$ with respect to Lebesgue measure, and the additional drift term $\Gamma$ is of a simpler form than $\Omega$ in (5). The resulting Metropolis–Hastings proposal mechanism is:

We refer to the resulting Metropolis–Hastings method as ‘position-dependent MALA’ or, more succinctly, ‘PMALA’.

The remainder of this section details two connections between the diffusions (5) and (9) when $A(\bm{X}_{t})=G^{-1}(\bm{X}_{t})$. In describing these connections the following equivalent forms for the $i$th components of $\Omega(\bm{X}_{t})$ and $\Gamma(\bm{X}_{t})$ will be helpful. For clarity of exposition we suppress explicit dependence on $X_{t}$ of all four of these quantities.

The first connection arises because the diffusion (5) on which both MMALA and the algorithm of Roberts and Stramer (2002) are based contains a transcription error. The term $\Omega$ should be multiplied by a factor of $1/2$, giving the diffusion

This can be viewed as a deterministic mapping of (2) onto a Riemannian manifold with metric tensor $G$, with the first term being the covariant drift, and the second and third corresponding to a Brownian motion on the manifold Kent (1978). However, the density $\pi^{*}$ is not given with respect to the Lebesgue measure, but instead respect to the dimensional volume or Hausdorff measure of the manifold, which is coordinate invariant. We refrain from discussing this in detail, but note that this is related to the density $\pi$ with respect to the Lebesgue measure via the area formula (Federer, 1969, Theorem 3.2.5),

The diffusions defined by (12) and (9) are equal.

The volatilities of the two diffusions are the same, so we need only compare the drift terms. Substituting (13) into (9) gives a diffusion where the $i$th component of the drift term is

which, using (10), is the $i$th component of the drift in (12).

Thus the diffusion (5) arises as a result of both an error in transcription and omitting the determinant factor when changing reference measures. Interestingly, in certain circumstances these two mistakes appear to cancel, so that (5) does, in fact, have the correct invariant distribution.

If $G(\bm{x})$ is chosen such that for any combination of $1\leq j,k,m\leq d$:

for all $\bm{x}$, then (5) and (9) represent the same diffusion.

Since the volatitilites and the multipliers of $\nabla\log\pi$ in the drift are identical for the two diffusions, we need only show that $\Omega_{i}=\Gamma_{i}$ for all $i$. From (14) the second term in (11) can be rewritten as

on relabelling $j\leftrightarrow k$. The result follows since $G^{-1}_{jm}=G^{-1}_{mj}$.

This property arises in certain simple cases, which suggests, perhaps, how this mistake has thus far remained undetected. If the process is univariate ($d=1$), then (14) holds trivially. More generally, it also holds if $G$ is the (continuous) Hessian matrix of some real-valued function: in particular, in the case of a natural exponential family, such as a generalised linear model (GLM) with canonical link, the Fisher information matrix used by Girolami and Calderhead (2011), is equal to the Hessian of the negative log-likelihood function.

In general, however, the diffusion (5) will not have the desired invariant density.

For some positive-valued, differentiable function $f$, set

It is then straightforward to show that $\Gamma=$ and $\Omega=[0,f^{\prime}(x_{2})/f(x_{2})]$, and hence the diffusions (5) and (9) have different drift coefficients. Moreover, the diffusion (5) can be written in the same form as that of (9); by matching the drift terms, it can be seen that the invariant density of (5) is actually proportional to $\pi(\bm{x})f(x_{2})$.

Experiments

We compared the performance of the MALA schemes across three of the scenarios considered in Girolami and Calderhead (2011): logistic regression on each of five different datasets; a stochastic volatility model; and a non-linear ODE model. As in Girolami and Calderhead (2011) we base the metric tensor, $G({\bf x})=A({\bf x})^{-1}$, on the expected Fisher information.

Initial tuning runs provided the optimal scaling parameter(s) ($\sqrt{h}$ in this article) in terms of ESS for each algorithm (on each dataset, where relevant). The initialisation, burn-in, and length of each Markov chain was exactly as in Girolami and Calderhead (2011), however we performed $100$ (rather than $10$) replicated runs for each chain.

Bayesian logistic regression and the non-linear ODE model are of most interest since in Girolami and Calderhead (2011) MMALA was found to outperform Riemann Manifold Hamiltonian Monte Carlo for these scenarios. Due to space considerations we therefore present detailed results for these scenarios; results for the stochastic volatility model showed the same pattern as for the non-linear ODE model. Where especially pertinent we provide brief details on the models themselves and the priors; for further details the reader is referred to Girolami and Calderhead (2011).

We perform Bayesian inference for a logistic regression model on each of five different datasets containing between $7$ and $25$ covariates. We choose a Gaussian prior for the parameter vector $\bm{\beta}\sim\mathcal{N}(\mathbf{0},\alpha I)$, so that with a design matrix is $X$ and link function $s(\cdot)$ the metric tensor is given by $G(\bm{\beta})=X^{T}\Lambda X+\alpha^{-1}I$, where $\Lambda$ is a diagonal matrix with elements $\Lambda_{i,i}=s(\bm{\beta}^{T}X^{T}_{i,\cdot})(1-s(\bm{\beta}^{T}X^{T}_{i,\cdot}))$. As noted above, this satisfies (14), so the diffusions on which PMALA and MMALA are based have the same law and we should expect the ESSs for these two algorithms to be the same up to Monte Carlo error.

For each Markov chain the ESS was computed for each parameter and the minimum, median and maximum of these was noted. Table 1 shows, for each algorithm and dataset, the means and their corresponding standard errors using the $100$ replicates. The CPU time and the mean (over replicates) minimum (over parameters) effective number of independent samples per second are also provided.

As expected, the ESSs for PMALA and MMALA are very similar. Since $\Gamma$ is computationally less costly to calculate than $\Omega$, PMALA is quicker and so obtains the larger ESS per second.

2 Non-linear differential equation model

We now consider the FitzHugh–Nagumo differential equations in Ramsay et al. (2007): $\dot{W}=c\left(W-{W^{3}}/{3}+R\right)$ and $\dot{R}=-(W-a+bT)/{c}.$ The simulated dataset and our independent priors for the parameter vector $(a,b,c)$ and the variance of the Gaussian observation noise are the same as those used in Girolami and Calderhead (2011). To be consistent with the appendix of Girolami and Calderhead (2011) and the associated Matlab code we assume $\beta\sim\mbox{Exp}(1)$.

Table 2 presents the mean ESS for each parameter with its standard error and shows that PMALA outperforms MMALA using this measure. CPU time and ESS/sec are also provided in the table; since each iteration of PMALA is also quicker, its advantage is even clearer when CPU time is accounted for.

Acknowledgements

T. Xifara was part-funded by North West Development Agency project N0003235 and the Greek State Scholarships Foundation. S. Livingstone is funded by a PhD Scholarship from Xerox Research Centre Europe. S. Byrne is funded by an EPSRC Postdoctoral Research Fellowship, EP/K005723/1. M. Girolami is funded by an EPSRC Established Career Research Fellowship, EP/J016934/1 and a Royal Society Wolfson Research Merit Award and is grateful to Prof. Jesus Sanz Serna for drawing his attention to the underlying Hausdorff measure of the diffusion (5) in a personal communication.

References