Information-geometric Markov Chain Monte Carlo methods using Diffusions

Introduction

There are three objectives to this article. The first is to introduce geometric concepts that have recently been employed in Monte Carlo methods based on Markov chains to a wider audience. The second is to clarify what a ‘diffusion on a manifold’ is, and how this relates to a diffusion defined on Euclidean space. Finally we review the state of the art in the field, and suggest avenues for further research.

The connections between some Monte Carlo methods commonly used in Statistics, Physics and application domains such as Econometrics, and ideas from both Riemannian and Information geometry were highlighted by Girolami & Calderhead , and the potential benefits demonstrated empirically. Two Markov chain Monte Carlo methods were introduced, the manifold Metropolis-adjusted Langevin algorithm and Riemannian manifold Hamiltonian Monte Carlo. Here we focus on the former for two reasons. First, the intuition for why geometric ideas can improve standard algorithms is the same in both cases. Second, the foundations of the methods are quite different, and since the focus of the article is on using geometric ideas to improve performance, we considered a detailed description of both to be unnecessary. It should be noted, however, that impressive empirical evidence exists for using Hamiltonian methods in some scenarios (e.g. ). We refer interested readers to .

We take an expository approach, providing a review of some necessary preliminaries from Markov chain Monte Carlo, diffusion processes and Riemannian geometry. We assume only a minimal familiarity with measure-theoretic probability. More informed readers may prefer to skip these sections. We then provide a full derivation of the Langevin diffusion on a Riemannian manifold, and offer some intuition for how to think about such a process. We conclude Section 4 by presenting the Metropolis-adjusted Langevin algorithm on a Riemannian manifold.

A key challenge in the geometric approach is which manifold to choose. We discuss this in Subsection 4.4, and review some candidates that have been suggested in the literature, along with the reasoning for each. Rather than provide a simulation study here, we instead reference studies where the methods we describe have been applied in Section 5. In Section 6 we discuss several open questions which we feel could be interesting areas of further research, and of interest to both theorists and practitioners.

Throughout $\pi(\cdot)$ will refer to an $n$-dimensional probability distribution, and $\pi(x)$ its density with respect to Lebesgue measure.

Markov Chain Monte Carlo

Markov chain Monte Carlo (MCMC) is a set of methods for drawing samples from a distribution $\pi(\cdot)$ defined on a measurable space $(\mathcal{X},\mathcal{B})$ whose density is only known up to some proportionality constant. Although the $i$th sample is dependent on the $(i-1)$th, the Ergodic Theorem ensures that for an appropriately constructed Markov chain with invariant distribution $\pi(\cdot)$, long-run averages are consistent estimators for expectations under $\pi(\cdot)$. As a result, MCMC methods have proven useful in Bayesian Statistical inference, where often the posterior density $\pi(x|y)\propto f(y|x)\pi_{0}(x)$ for some parameter $x$ is only known up to a constant . Here we briefly introduce some concepts from general state space Markov chain theory together with a short overview of MCMC methods. The exposition follows .

for all $A\in\mathcal{B}$. If $P(x,\cdot)$ admits a density $p(x^{\prime}|x)$, this can be equivalently written:

for any $A\in\mathcal{B}$. Certain conditions are required for (5) to hold, but for all Markov chains presented here these are satisfied (though see ).

A useful condition which is sufficient (though not necessary) for $\pi(\cdot)$ to be an invariant distribution is reversibility, which can be shown by the relation

Integrating over both sides with respect to $x$ we recover (4). In words, a chain is reversible if at stationarity the probability that $x_{i}\in A$ and $x_{i+1}\in B$ is equal to the probability that $x_{i+1}\in A$ and $x_{i}\in B$. The relation (6) will be the primary tool used to construct Markov chains with a desired invariant distribution in the next section.

with probability one . This is a Markov chain analogue to the Law of Large numbers.

It follows directly from the Kipnis-Varadhan Theorem that an estimator $\hat{t}_{m}$ from a reversible Markov chain for which $X_{0}\sim\pi(\cdot)$ satisfies:

provided that $\sum_{i=1}^{\infty}i|\rho^{(0,i)}|<\infty$, where $\rho^{(0,i)}=\text{Corr}_{\pi}[g(X_{0}),g(X_{i})]$. We will refer to the constant $\tau$ as the autocorrelation time for the chain .

Equation (8) implies that for large enough $m$, $\text{Var}[\hat{t}_{m}]\approx\tau\text{Var}[\bar{t}_{m}]$. In practical applications, the sum in (8) is truncated to the first $p-1$ realisations of the chain, where $p$ is the first instance at which $|\rho^{(0,p)}|<\epsilon$ for some $\epsilon>0$ .

The measures arising from (8) give some intuition for what sort of Markov chain gives rise to efficient estimators. However, in practice the chain will never be at stationarity. So we also assess Markov chains according to how far away they are from this point. For this, we need to measure how close $P^{m}(x_{0},\cdot)$ is from $\pi(\cdot)$, which requires a notion of distance between probability distributions.

Although there are several appropriate choices , a common option in the Markov chain literature is the Total variation distance

which informally gives the largest possible difference between the probabilities of a single event in $\mathcal{B}$ according to $\mu(\cdot)$ and $\nu(\cdot)$. If both distributions admit densities, (9) can be written

which is proportional to the $L_{1}$ distance between $\mu(x)$ and $\nu(x)$. Our metric $\|\cdot\|_{TV}\in$, with $\|\cdot\|_{TV}=1$ for distributions with disjoint supports and $\|\mu(\cdot)-\nu(\cdot)\|_{TV}=0$ implying $\mu(\cdot)\equiv\nu(\cdot)$.

Typically for an unbounded $\mathcal{X}$ the distance $\|P^{m}(x_{0},\cdot)-\pi(\cdot)\|_{TV}$ will depend on $x_{0}$ for any finite $m$. So bounds on the distance are often sought via some inequality of the form

A Markov chain is called geometrically ergodic if $f(m)=r^{m}$ in (11) for some $0<r<1$. If in addition to this $V$ is bounded above, the chain is called uniformly ergodic. Intuitively, if either condition holds then the distribution of $X_{m}$ will converge to $\pi(\cdot)$ geometrically quickly as $m$ grows, and in the uniform case this rate is independent of $x_{0}$. As well as providing some (often qualitative if $r$ is unknown) bounds on the convergence rate of a Markov chain, geometric ergodicity implies that a central limit theorem exists for estimators of the form $\hat{t}_{m}$, so that the shape of the distribution is asymptotically Gaussian. For more detail on this see .

In practice several approximate methods also exist to assess whether a chain is close enough to stationarity for long-run averages to provide suitable estimators (e.g. ). The MCMC practitioner also uses a variety of visual aids to judge whether an estimate from the chain will be appropriate for his or her needs.

2 Markov Chain Monte Carlo

Now that we have introduced Markov chains we turn to simulating them. The objective here is to devise a method for generating a Markov chain which has a desired limiting distribution $\pi(\cdot)$. In addition we would strive for the convergence rate to be as fast as possible, and the effective sample size to be suitably large relative to the number of iterations. Of course, the computational cost of performing an iteration is also an important practical consideration. Ideally any method would also require limited problem-specific alterations, so that practitioners are able to use it with as little knowledge of the inner workings as is practical.

Although other methods exist for constructing chains with a desired limiting distribution, a popular choice is the Metropolis–Hastings algorithm . At iteration $i$, a sample is drawn from some candidate transition kernel $Q(x_{i-1},\cdot)$, and then either accepted or rejected (in which case the state of the chain remains $x_{i-1}$). We focus here on the case where $Q(x_{i-1},\cdot)$ admits a density $q(x^{\prime}|x_{i-1})$ for all $x_{i-1}\in\mathcal{X}$ (though see ). In this case a single step is shown below (the wedge notation $a\wedge b$ denotes the minimum of $a$ and $b$).

The ‘acceptance rate’ $\alpha(x_{i-1},x^{\prime})$ governs the behaviour of the chain. If it is typically close to one then many proposed moves are accepted and so the current value in the chain is constantly changing. If it is on average close to zero then many proposals are rejected so the chain will remain in the same place for many iterations. However, $\alpha\approx 1$ is typically not ideal, often resulting in a large autocorrelation time (see below). The challenge in practice is to find the right acceptance rate to balance these two extremes.

Combining the ‘proposal’ and ‘acceptance’ steps, the transition kernel for the resulting Markov chain is

is the average probability that a draw from $Q(x,\cdot)$ will be rejected, and $\delta_{x}(A)=1$ if $x\in A$ and zero otherwise. A Markov chain defined in this way will have $\pi(\cdot)$ as an invariant distribution, since the chain is reversible for $\pi(\cdot)$. We note here that

in the case that the proposed move is accepted, and that if the proposed move is rejected then $x_{i}=x_{i-1}$ so the chain is reversible for $\pi(\cdot)$. It can be shown that $\pi(\cdot)$ is also the limiting distribution for the chain .

The convergence rate and autocorrelation time of a chain produced by the algorithm are dependent on both the choice of proposal $Q(x_{i-1},\cdot)$ and the target distribution $\pi(\cdot)$. For simple forms of the latter, less consideration is required when choosing the former. A broad objective among researchers in the field is to find classes of proposal kernels that produce chains which converge and mix quickly for a large class of target distributions. We first review a simple choice before discussing one which is more sophisticated, and will be the focus of the rest of the article.

3 Random walk proposals

An extremely simple choice for $Q(x,\cdot)$ is one for which:

where $\|\cdot\|$ denotes some appropriate norm on $\mathcal{X}$, meaning the proposal is symmetric. In this case, the acceptance rate reduces to:

In addition to simplifying calculations, (14) strengthens the intuition for the method, since proposed moves with higher density under $\pi(\cdot)$ will always be accepted. A typical choice for $Q(x,\cdot)$ is $\mathcal{N}(x,\lambda^{2}\Sigma)$, where the matrix $\Sigma$ is often chosen in an attempt to match the correlation structure of $\pi(\cdot)$, or simply taken as the identity . The tuning parameter $\lambda$ is the only other user-specific input required.

Much research has been conducted into properties of the random walk Metropolis algorithm (RWM). It has been shown that the optimal acceptance rate for proposals tends to $0.234$ as the dimension $n$ of the state space $\mathcal{X}$ tends to $\infty$ for a wide class of targets (e.g. ). The intuition for an optimal acceptance rate is to find the right balance between the distance of proposed moves and the chances of acceptance. Increasing the former will reduce the autocorrelation in the chain if the proposal is accepted, but if it is rejected the chain will not move at all, so autocorrelation will be high. Random walk proposals are sometimes referred to as blind (e.g. ), as no information about $\pi(\cdot)$ is used when generating proposals, so typically very large moves will result in a very low chance of acceptance, while small moves will be accepted but result in very high autocorrelation for the chain. Figure 1 demonstrates this in the simple case where $\pi(\cdot)$ is a one dimensional $\mathcal{N}(0,1^{2})$ distribution.

Several authors have also shown that for certain classes of $\pi(\cdot)$ the tuning parameter $\lambda$ should be chosen such that $\lambda^{2}\propto n^{-1}$ so that $\alpha\nrightarrow 0$ as $n\rightarrow\infty$ . Because of this we say that algorithm efficiency ‘scales’ $O(n^{-1})$ as the dimension $n$ of $\pi(\cdot)$ increases.

Ergodicity results for a Markov chain constructed using the RWM algorithm also exist . At least exponentially light tails are a necessity for $\pi(x)$ for geometric ergodicity, which means that $\pi(x)/e^{-\|x\|}\to c$ as $\|x\|\to\infty$, for some constant $c$. For super-exponential tails (where $\pi(x)\to 0$ at a faster than exponential rate), additional conditions are required . We demonstrate with a simple example why heavy-tailed forms of $\pi(x)$ pose difficulties here (where $\pi(x)\to 0$ at a rate slower then $e^{-\|x\|}$).

Example: Take $\pi(x)\propto 1/(1+x^{2})$, so that $\pi(\cdot)$ is a Cauchy distribution. Then if $X^{\prime}\sim\mathcal{N}(x,\lambda^{2})$, the ratio $\pi(x^{\prime})/\pi(x)=(1+(x^{\prime})^{2})/(1+x^{2})\to 1$ as $|x|\to\infty$. So if $x_{0}$ is far away from , the Markov chain will dissolve into a random walk, with almost every proposal being accepted.

In the remainder of the article we will primarily discuss another approach to choosing $Q$ which has been shown empirically and in some cases theoretically to be superior to the RWM algorithm, though it should be noted that random walk proposals are still widely used in practice and are often sufficient for more straightforward problems .

Diffusions

In MCMC we are concerned with discrete time processes. However, often there are benefits to first considering a continuous time process with properties we desire. One such is that some continuous time processes can be specified via a form of differential equation. In this section we derive a choice for a Metropolis–Hastings proposal kernel based on approximations to diffusions, those continuous-time $n$-dimensional Markov processes $(X_{t})_{t\geq 0}$ for which any sample path $t\mapsto X_{t}(\omega)$ is a continuous function with probability one. For any fixed $t$, we assume $X_{t}$ is a random variable taking values on the measurable space $(\mathcal{X},\mathcal{B})$ as before. In the next section we provide some preliminaries, followed by an introduction to our main object of study, the Langevin diffusion.

We focus on the class of time-homogeneous Itô diffusions, whose dynamics are governed by a stochastic differential equation of the form

implying that the drift dictates how the mean of the process changes over a small time interval, and if we define the process $(M_{t})_{t\geq 0}$ through the relation

giving the stochastic part of the relationship between $X_{t+\triangle t}$ and $X_{t}$ for small enough $\triangle t$, see e.g. .

While (15) is often a suitable description of an Itô diffusion, it can also be characterised through an infinitessimal generator $\mathcal{A}$, which describes how functions of the process are expected to evolve. We define this partial differential operator through its action on a function $f\in C_{0}(\mathcal{X})$ as

though $\mathcal{A}$ can be associated with the drift and volatility of $(X_{t})_{t\geq 0}$ by the relation

where $V_{ij}(x)$ denotes the component in row $i$ and column $j$ of $\sigma(x)\sigma(x)^{T}$ . Later on we shall use the generator characterisation of a diffusion to generalise it in some sense.

As in the discrete case, we can describe the transition kernel of a continuous time Markov process $P^{t}(x_{0},\cdot)$. In the case of an Itô diffusion, $P^{t}(x_{0},\cdot)$ admits a density $p_{t}(x|x_{0})$, which in fact varies smoothly as a function of $t$. The Fokker–Planck equation describes this variation in terms of the drift and volatility, and is given by

Although typically the form of $P^{t}(x_{0},\cdot)$ is unknown the expectation and variance of $X_{t}\sim P^{t}(x_{0},\cdot)$ are given by the integral equations:

where the second of these is a result of the Itô isometry . Continuing the analogy, a natural question is whether a diffusion process has an invariant distribution $\pi(\cdot)$, and whether

for any $A\in\mathcal{B}$ and any $x_{0}\in\mathcal{X}$, in some sense. For a large class of diffusions (which we confine ourselves to) this is in fact the case , and in addition (21) provides a means of finding $\pi(\cdot)$ given $b$ and $\sigma$. Setting the left-hand side of (21) to zero gives

which can be solved to find $\pi(\cdot)$.

2 Langevin diffusions

Given (23) our goal becomes clearer: find drift and volatility terms so that the resulting dynamics describe a diffusion which converges to some user-defined invariant distribution $\pi(\cdot)$. This process can then be used as a basis for choosing $Q$ in a Metropolis–Hastings algorithm. The Langevin diffusion, first used to describe the dynamics of molecular systems , is such a process, given by the solution to the stochastic differential equation

which is a sufficient condition for (23) to hold. So for any case in which $\pi(x)$ is suitably regular so that $\nabla\log\pi(x)$ is well-defined and the derivatives in (23) exist, we can use (24) to construct a diffusion which has invariant distribution $\pi(\cdot)$.

Roberts & Tweedie give sufficient conditions on $\pi(\cdot)$ under which a diffusion $(X_{t})_{t\geq 0}$ with dynamics given by (24) will be ergodic, meaning

as $t\to\infty$, for any $x_{0}\in\mathcal{X}$. They remark that these conditions ‘should be appropriate for virtually all commonly encountered target densities’ .

3 Metropolis-adjusted Langevin algorithm

We can use Langevin diffusions as a basis for MCMC in many ways, but a popular variant is known as the Metropolis-adjusted Langevin algorithm (MALA), whereby $Q(x,\cdot)$ is constructed through an Euler–Maruyama discretisation of (24) and used as a candidate kernel in a Metropolis–Hastings algorithm. The resulting $Q$ is

where $\lambda$ is again a tuning parameter.

Before we discuss the theoretical properties of the approach, we first offer intuition for the dynamics. From (27) it can be seen that Langevin-type proposals comprise a deterministic shift towards a local mode of $\pi(x)$, combined with some random additive Gaussian noise, with variance $\lambda^{2}$ for each component. The relative weights of the deterministic and random parts are fixed, given as they are by the parameter $\lambda$. Typically if $\lambda^{1/2}\gg\lambda$ then the random part of the proposal will dominate, and vice versa in the opposite case, though this also depends on the form of $\nabla\log\pi(x)$ .

Again since this is a Metropolis–Hastings method, choosing $\lambda$ is a balance between proposing large enough jumps and ensuring that a reasonable proportion are accepted. It has been shown that in the limit as $n\to\infty$ the optimal acceptance rate for the algorithm is $0.574$ for forms of $\pi(\cdot)$ which either have independent and identically distributed components or whose components only differ by some scaling factor . In these cases, as $n\to\infty$ the parameter $\lambda$ must be $\propto n^{-1/3}$, so we say algorithm efficiency scales $O(n^{-1/3})$. Note that these results compare favourably with the $O(n^{-1})$ scaling of the random walk algorithm.

Convergence properties of the method have also been established. Roberts & Tweedie highlight some cases in which MALA is either geometrically ergodic or not. Typically results are based on the tail behaviour of $\pi(x)$. If these tails are heavier than exponential, then the method is typically not geometrically ergodic, and similarly if the tails are lighter than Gaussian. However, in the in between case the converse is true. We again offer two simple examples for intuition here.

Example: Take $\pi(x)\propto 1/(1+x^{2})$ as in the previous example. Then $\nabla\log\pi(x)=-2x/(1+x^{2})^{2}\to 0$ as $|x|\to\infty$. So if $x_{0}$ is far away from 0, then the MALA will be approximately equal to the RWM algorithm, and so will also dissolve into a random walk.

Example: Take $\pi(x)\propto e^{-x^{4}}$. Then $\nabla\log\pi(x)=-4x^{3}$ and $X^{\prime}\sim\mathcal{N}(x-4\lambda^{2}x^{3},\lambda^{2})$. So for any fixed $\lambda$, there exists $c>0$ such that for $|x_{0}|>c$ we have $|4\lambda^{2}x^{3}|>>x$ and $|x-4\lambda^{2}x^{3}|>>\lambda$, suggesting that MALA proposals will quickly spiral further and further away from any neighbourhood of , and hence nearly all will be rejected.

For cases where there is strong correlation between elements of $x$ or each element has a different marginal variance, the MALA can also be ‘pre-conditioned’ in a similar way to the RWM, so that the covariance structure of proposals more accurately reflects that of $\pi(x)$ . In this case, proposals take the form

where $\lambda$ is again a tuning parameter. It can be shown that provided $\Sigma$ is a constant matrix, $\pi(x)$ is still the invariant distribution for the diffusion on which (28) is based .

Geometric concepts in Markov Chain Monte Carlo

Ideas from Information geometry have been successfully applied to Statistics from as early as . More widely, other geometric ideas have also been applied, offering new insight into common problems (e.g. ). A survey is given in . In this section we suggest why some ideas from differential geometry may be beneficial for sampling methods based on Markov chains. We then review what is meant by a ‘diffusion on a manifold’, before turning to the specific case of (24). After this we discuss what can be learned from work in Information geometry in this context.

For sampling methods based on Markov chains which explore the space locally, like the RWM and MALA, it may be advantageous to instead impose a different metric structure on the space $\mathcal{X}$, so that some points are drawn closer together and others pushed further apart. Intuitively, one can picture distances in the space being defined such that if the current position in the chain is far from an area of $\mathcal{X}$ which is ‘likely to occur’ under $\pi(\cdot)$, then the distance to such a typical set could be reduced. Similarly, once this region is reached, the space could be ‘stretched‘ or ‘warped’ so that it is explored as efficiently as possible.

2 Preliminaries

and provides a means to define a metric on the manifold as $d(x,y)=\inf\left\{L(\gamma_{M}):\gamma_{M}(0)=x,\gamma_{M}(1)=y\right\}$.

We seek to make these ideas more concrete through an example, the graph of a function $f(x_{1},x_{2})$ of two variables $x_{1}$ and $x_{2}$. The resulting map $r$ is

If $G(x)=I$ then this reduces to Lebesgue measure.

3 Diffusions on manifolds

By a ‘diffusion on a manifold’ in local coordinates, we actually mean a diffusion on Euclidean space which, when mapped onto $M$ through $r$, becomes the desired diffusion along the manifold. For example, a ‘Brownian motion on a sphere’ means the diffusion on Euclidean space that, when mapped onto the sphere, produces Brownian motion along the surface. If a realisation of Brownian motion was drawn in wet ink on the surface of the sphere, and the sphere was then rolled over a piece of flat paper, the resulting imprint on the paper would be a realisation of Brownian motion on Euclidean space . However, the mapping $r$ does not correspond to this rolling procedure, so the pre-image of the Brownian motion on the sphere under $r$ will not be Brownian motion on Euclidean space.

Our goal, therefore, is to define a diffusion on Euclidean space which, when mapped onto a manifold through $r$, becomes the Langevin diffusion described in (24) by the above procedure. Such a diffusion takes the form

where those objects marked with a tilde must be defined appropriately. The next few paragraphs are technical, and readers aiming to simply grasp the key points may wish to skip to the end of this Subsection.

On a manifold, the gradient operator $\nabla_{M}$ can still be defined such that the inner product $g_{p}(\nabla_{M}f(x),u)=D_{u}[f(x)]$. Setting $\nabla_{M}=G(x)^{-1}\nabla$ gives

which is equal to the directional derivative along $u$ as required.

The divergence of a vector field $v$ defined on a manifold $M$ at the point $p\in M$ is defined as:

where $e_{i}$ denotes the $i$th basis vector for the tangent space $T_{p}M$ at $p\in M$, and $v_{i}$ denotes the $i$th coefficient. This can be written in local coordinates (see Appendix A) as

Combining these two operators, we can define a generalisation of the Laplace operator, known as the Laplace–Beltrami operator (e.g. ), as:

The generator of a Brownian motion on $M$ is $\triangle_{LB}/2$ . Using (20) the resulting diffusion has dynamics given by

Those familiar with the Itô formula will not be surprised by the additional drift term $\Omega(X_{t})$. As Itô integrals do not follow the chain rule of ordinary calculus, non-linear mappings of martingales like $(B_{t})_{t\geq 0}$ typically result in drift terms being added to the dynamics (e.g. ).

Putting these three elements together, (39) becomes

It can be shown that this diffusion has invariant Lebesgue density $\pi(x)$ as required . Intuitively when a set is mapped onto the manifold, distances are changed by a factor $\sqrt{G(x)}$. So to end up with the initial distances, they must first be changed by a factor of $\sqrt{G^{-1}(x)}$ before the mapping, which explains the volatility term in (46).

The resulting Metropolis–Hastings proposal kernel for this ‘MALA on a manifold’ was clarified in , and is given by

where $\lambda^{2}$ is a tuning parameter. The nonlinear drift term here is slightly different to that reported in , for reasons discussed in .

4 Choosing a Metric

We now turn to the question of which manifold to choose, or equivalently how to choose $G(x)$. In this section we sometimes switch notation slightly, denoting the target density $\pi(x|y)$, as some of the discussion is directed towards Bayesian inference, where $\pi(\cdot)$ is the posterior distribution for some parameter $x$ after observing some data $y$. The problem statement is: what is an appropriate choice of distance between points in the sample space of a given probability distribution?

A related (but distinct) question is how to define a distance between two probability distributions from the same parametric family but with different parameters. This has been a key theme in Information geometry, explored by Rao and others for many years. Although generic measures of distance between distributions (such as total variation) are often appropriate, based on Information-theoretic principles one can deduce that for a given parametric family $\{p_{x}(y):x\in\mathcal{X}\}$, it is in some sense natural to consider this ‘space of distributions’ to be a manifold, where the Fisher information is the matrix $G(x)$ (with the $\alpha=0$ connection employed, see for details).

Because of this, Girolami & Calderhead proposed a variant of the Fisher metric for geometric Markov chain Monte Carlo, as

where $\pi(x|y)\propto f(y|x)\pi_{0}(x)$ is the target density, $f$ denotes the likelihood and $\pi_{0}$ the prior. The metric is tailored to Bayesian problems, which are a common use for MCMC, so the Fisher information is combined with the negative Hessian of the log-prior. One can also view this metric as the expected negative Hessian of the log target, since this naturally reduces to (49).

The motivation for a Hessian-style metric can also be understood from studying MCMC proposals. From (48) and by the same logic as for general pre-conditioning methods , the objective is to choose $G^{-1}(x)$ to match the covariance structure of $\pi(x|y)$ locally. If the target density were Gaussian with covariance matrix $\Sigma$, then

In the non-Gaussian case, the negative Hessian is no longer constant, but we can imagine that it matches the correlation structure of $\pi(x|y)$ locally at least. Such ideas have been discussed in the geostatistics literature previously . One problem with simply using (50) to define a metric is that unless $\pi(x|y)$ is log-concave the negative Hessian will not be globally positive-definite, although Petra et al. conjecture that it may be appropriate for use in some realistic scenarios, and suggest some computationally efficient approximation procedures .

Example: Take $\pi(x)\propto 1/(1+x^{2})$, and set $G(x)=-\partial^{2}\log\pi(x)/\partial x^{2}$. Then $G^{-1}(x)=(1+x^{2})^{2}/(2-2x^{2})$, which is negative if $x^{2}>1$, so unusable as a proposal variance.

Girolami & Calderhead use the Fisher metric in part to counteract this problem. Taking expectations over the data ensures that the likelihood contribution to $G(x)$ in (49) will be positive (semi-)definite globally (e.g. ), so provided a log-concave prior is chosen then (49) should be a suitable choice for $G(x)$. Indeed, Girolami & Calderhead provide several examples in which geometric MCMC methods using this Fisher metric perform better than their ‘non-geometric’ counterparts.

where $\alpha$ is a tuning parameter (typically chosen to be as large as possible for which eigenvalues remain non-zero numerically). The function $t_{\alpha}$ acts as an absolute value function, but also uplifts eigenvalues which are close to zero to $\approx 1/\alpha$. It should be noted that while the Fisher metric is only defined for models in which a likelihood is present, and for which the expectation is tractable, the SoftAbs metric can be found for any target distribution $\pi(\cdot)$.

Many authors (e.g. ) have noted that for many problems the terms involving derivatives of $G(x)$ are often small, and so it is not always worth the computational effort of evaluating them. Girolami & Calderhead propose the simplified manifold MALA, in which proposals are of the form

Using this method means derivatives of $G(x)$ are no longer needed, so more pragmatic ways of regularising the Hessian are possible. One simple approach would be to take the absolute values of each eigenvalue, giving $G(x)=U^{T}|D|U$, where $H(x)=U^{T}DU$ is the negative Hessian and $|D|$ is a diagonal matrix with $\{|D|\}_{ii}=|\lambda_{i}|$ (this approach may fall into difficulties if eigenvalues are numerically zero). Another would be choose $G(x)$ as the ‘nearest’ positive-definite matrix to the negative Hessian, according to some distance metric on the set of $n\times n$ matrices. The problem has in fact been well-studied in mathematical finance, in the context of finding correlations using incomplete data sets , and tackled using distances induced by the Frobenius norm. Approximate solution algorithms are discussed in Higham . We again provide a simple example suggesting that a ‘Hessian-style metric’ can alleviate some of the difficulties associated with heavy-tailed target densities.

Example: Take $\pi(x)\propto 1/(1+x^{2})$, and set $G(x)=|-\partial^{2}\log\pi(x)/\partial x^{2}$—. Then $G^{-1}(x)\nabla\log\pi(x)=-x(1+x^{2})/|1-x^{2}|$, which no longer tends to as $|x|\to\infty$, suggesting a manifold variant of MALA with a Hessian-style metric may avoid some of the pitfalls of the standard algorithm. Note that the drift may become very large if $|x|\approx 1$, but since this event occurs with probability we do not see it as a major cause for concern.

Example: Take $\pi(x)\propto e^{-x^{4}}$, and set $G(x)=|-\partial^{2}\log\pi(x)/\partial x^{2}|$. Then $G^{-1}(x)\nabla\log\pi(x)=-x/3$, which is $O(x)$, so alleviates the problem of spiralling proposals for light-tailed targets demonstrated by MALA in an earlier example.

Other choices for $G(x)$ have been proposed which are not based on the Hessian. These have the advantage that gradients need not be computed (either analytically or using computational methods). Sejdinovic et al. propose a Metropolis–Hastings method which can be viewed as a geometric variant of the RWM, where the choice for $G(x)$ is based on mapping samples to an appropriate feature space, and performing principal component analysis on the resulting features to choose a local covariance structure for proposals.

If we consider the RWM with Gaussian proposals to be an Euler–Maruyama discretisation of Brownian motion on a manifold, then proposals will take the form $Q(x,\cdot)\equiv\mathcal{N}(x+\lambda^{2}\Omega(x),\lambda^{2}G^{-1}(x))$. If we assume (like in the simplified manifold MALA) that $\Omega(x)\approx 0$, then we have proposals centred at the current point in the Markov chain with a local covariance structure (the full Hastings acceptance rate must now be used as $q(x^{\prime}|x)\neq q(x|x^{\prime})$ in general).

As no gradient information is needed, the Sejdinovic et al. metric can be used in conjunction with the pseudo-marginal MCMC algorithm, so that $\pi(x|y)$ need not be known exactly. Examples from the article demonstrate the power of the approach .

Survey of applications

Rather than conduct our own simulation study, we instead highlight some cases in the literature where geometric MCMC methods have been used with success.

Martin et al. consider Bayesian inference for a Statistical inverse problem, in which a surface explosion causes seismic waves to travel down into the ground (the subsurface medium). Often the properties of the subsurface vary with distance from ground level or because of obstacles in the medium, in which case a fraction of the waves will scatter off these boundaries and be reflected back up to ground level at later times. The observations here are the initial explosion and the waves which return to the surface, together with return times. The challenge is to infer the properties of the subsurface medium from this data. The authors construct a likelihood based on the wave equation for the data, and perform Bayesian inference using a variant of the manifold MALA. Figures are provided showing the local correlations present in the posterior, and therefore highlighting the need for an algorithm which can navigate the high density region efficiently. Several methods are compared in the paper, but the variant of MALA which incorporates a local correlation structure is shown to be the most efficient, particularly as the dimension of the problem increases .

Calderhead & Girolami dealt with two models for biological phenomena based on nonlinear dynamical systems. A model of circadian control in the Arabidopsis thaliana plant comprised a system of six nonlinear differential equations, with twenty two parameters to be inferred. Another model for cell signalling consisted of a system of six nonlinear differential equations with eight parameters, with inference complicated by the fact that observations of the model are not recorded directly . The resulting inference was performed using RWM, MALA and geometric methods, with the results highlighting the benefits of taking the latter approach. The simplified variant of MALA on a manifold is reported to have produced the most efficient inferences overall, in terms of effective sample size per unit of computational time.

Stathopoulos & Girolami considered the problem of inferring parameters in Markov jump processes. In the paper a linear noise approximation is shown, which can make inference in such models more straightforward, enabling an approximate likelihood to be computed. Models based on chemical reaction dynamics are considered; one such from chemical kinetics contained four unknown parameters, another from gene expression consisting of seven. Inference was performed using the RWM, the simplified manifold MALA and Hamiltonian methods, with the MALA reported as most efficient according to the chosen diagnostics. The authors note that the simplified manifold method is both conceptually simple and able to account for local correlations, making it an attractive choice for inference .

Konukoglu et al. designed a method for personalising a generic model for a physiological process to a specific patient, using clinical data. The personalisation took the form of patient-specific parameter inference. The authors highlight some of the difficulties of this task in general, including the complexity of the models and relative sparsity of the datasets, which often result in a parameter identifiability issue . The example discussed in the paper is the Eikonal-Diffusion model describing electrical activity in cardiac tissue, which results in a likelihood for the data based on a nonlinear partial differential equation, combined with observation noise . A method for inference was developed by first approximating the likelihood using a spectral representation, and then using geometric MCMC methods on the resulting approximate posterior. The method was first evaluated on synthetic data, and the authors note that inference took ‘less than five minutes’ to provide accurate estimates. The method was then repeated on clinical data taken from a study for ventricular tachycardia radio-frequency ablation .

Discussion

The geometric viewpoint in not necessary to understand manifold variants of the MALA. Indeed, several authors have discussed these algorithms without considering them to be ‘geometric’, rather simply Metropolis–Hastings methods in which proposal kernels have a position-dependent covariance structure. We do not claim that the geometric view is the only one that should be taken. Our goal is merely to point out that such position-dependent methods can often be viewed as methods defined on a manifold, and that studying the structure of the manifold itself may lead to new insights on the methods. For example, taking the geometric viewpoint and noting the connection with Information geometry enabled Girolami & Calderhead to adopt the Fisher metric for calculations . We list here a few open questions that the geometric viewpoint may help shed some insight on.

Computationally minded readers will have noted that using position-dependent covariance matrices adds a significant computational overhead in practice, with additional $O(n^{3})$ matrix inversions required at each step of the corresponding Metropolis–Hastings algorithms. Clearly there will be many problems for which the matrix $G(x)$ does not change very much, and therefore choosing a constant covariance $G^{-1}(x)=\Sigma$ may result in a more efficient algorithm overall. Geometrically, this would correspond to a manifold with low curvature. It may be that geometric ideas could be used to understand whether the manifold is flat enough that a constant choice of $G(x)$ is sufficient. To make sense of this truly would require a relationship between curvature, an inherently local property, and more global statements about the manifold. Many results in differential geometry, beginning with the celebrated Gauss-Bonnet theorem, have previously related global and local properties in this way . It is unknown to the authors whether results exist relating the curvature of a manifold to some global property, but this is an interesting avenue for further research.

A related question is when to choose the simplified manifold MALA over the full method. Problems in which the term $\|\Lambda(x)\|$ is sufficient large to warrant calculation correspond to those for which the manifold has very high curvature in many places, so again making some global statement related to curvature could help here.

Although there is a reasonable intuitive argument for why the Hessian is an appropriate starting point for $G(x)$, the lack of positive-definiteness may be seen as a cause for concern by some. After all, it could be argued that if the curvature is not positive-definite in a region, then how can it be a reasonable approximation to the local covariance structure. Indeed, for target densities of the form $\pi(x)\propto e^{-|x|}$, the Hessian is everywhere equal to zero! Much work in Information geometry has centred on the geometry of Hessian structures , and some insights from this field may help better understand the question of what appropriate metric to choose is.

Some recent work in high-dimensional inference has centred on defining MCMC methods for which efficiency scales $O(1)$ with respect to the dimension $n$ of $\pi(\cdot)$ . In the case where $X$ takes values in some infinite-dimensional function space, this can be done provided a Gaussian prior measure is defined for $X$. A striking result from infinite-dimensional probability spaces is that two different probability measures defined over some infinite dimensional space have a striking tendency to have disjoint supports . The key challenge for MCMC is to define transition kernels for which proposed moves are inside the support for $\pi(\cdot)$. A straight-forward approach is to define proposals for which the prior is invariant, since the likelihood contribution to the posterior typically will not alter its support from that of the prior . However, the posterior may still look very different from the prior, as noted in , so this proposal mechanism, though $O(1)$, can still result in slow exploration. Understanding the geometry of the support, and defining methods which incorporate the likelihood term but also respect this geometry so as to ensure proposals remain in the support of $\pi(\cdot)$, is an intriguing research proposition.

The methods reviewed in this paper are based on first order Langevin diffusions. Algorithms have also been developed which are based on second order Langevin diffusions, in which a stochastic differential equation governs the behaviour of the velocity of a process . A natural extension to the work of Girolami & Calderhead and Xifara et al. would be to map such diffusions onto a manifold and derive Metropolis–Hastings proposal kernels based on the resulting dynamics. The resulting scheme would be a generalisation of , though the most appropriate discretisation scheme for a second order process to facilitate sampling is unclear, and perhaps a question worthy of further exploration.

Alongside these geometric problems, we can also discuss geometric MCMC methods from a statistical perspective. The last example given in the previous section hinted that the manifold MALA may cope better with target distributions with heavy tails. In fact, Latuszynski et al. have shown that in one dimension, the manifold MALA is geometrically ergodic for a class of targets of the form $\pi(x)\propto\exp(-|x|^{\beta})$ for any choice of $\beta\neq 1$. This incorporates cases where tails are heavier than exponential and lighter than Gaussian, two scenarios under which geometric ergodicity fails for the MALA.

Finding optimal acceptance rates and scaling of $\lambda$ with dimension are two other related challenges. In this case the picture is more complex. Traditional results have been shown for Metropolis–Hastings methods in the case where target distributions are independent and identically-distributed, or some other suitable symmetry and regularity in the shape of $\pi(\cdot)$. Manifold methods are, however, specifically tailored to scenarios in which this is not the case, scenarios in which there is high correlation between components of $x$ which changes depending on the value of $x$. It is less clear how to proceed with finding relevant results which can serve as guidelines to practitioners here. Indeed, Sherlock notes that a requirement for optimal acceptance rate results for the RWM to be appropriate is that the curvature of $\pi(x)$ does not change too much, yet this is the very scenario in which we would want to use a manifold method.

Conclusions

We have discussed the merits of viewing the sample space of a statistical model as a Riemannian manifold for Markov chain Monte Carlo. We focused on the Metropolis-adjusted Langevin algorithm, and provided a full exposition of the relevant Markov chain and Riemannian geometry background in order to understand the method. After this we provided a full geometric derivation of a Langevin diffusion on a Riemannian manifold, and resulting Metropolis–Hastings algorithms. In the previous section, we have highlighted several open questions related to the emerging field of geometric Markov chain Monte Carlo, which we hope will inspire innovative new research in the field.

Acknowledgements

S. Livingstone is funded by a PhD Scholarship from Xerox Research Centre Europe. M. Girolami is funded by an EPSRC Established Career Research Fellowship, EP/J016934/1 and a Royal Society Wolfson Research Merit Award. The authors thank Michael Epstein and Simon Byrne for proofreading the article and giving useful suggestions.

Appendix A Vector fields and the Covariant Derivative

For any smooth manifold $M$, the set of all tangent vectors to points on $M$ is known as the tangent bundle, and denoted $TM$.

The covariant derivative $D^{c}$ is defined so as to account for these shortcomings. When considering differentiation along a vector $u^{*}\notin T_{p}M$, $u^{*}$ is simply projected onto the tangent space. The derivative with respect to any $u\in T_{p}M$ can now be decomposed into a linear combination of derivatives of basis vectors and vector components

where the argument $p$ has been dropped but is implied for both components and local basis vectors. The operator $D^{c}_{u}[v]$ is defined to be linear in both $u$ and $v$ and satisfy the product rule , so (53) can be decomposed into

The operator $D^{c}$ need therefore only be defined along the direction of basis vectors $\mathbf{e}_{i}$ and for vector component $v^{i}$ and basis vector $\mathbf{e}_{i}$ arguments.

The divergence of a vector field $v\in\Gamma(TM)$ at the point $p\in M$ is given by

where again repeated indices are summed. As has been previously stated, if a metric $g$ and coordinate chart is chosen for $M$, the Christoffel symbols can be written in terms of the matrix $G(x)$. In this case

References