Information-geometric Markov Chain Monte Carlo methods using Diffusions

Recent work incorporating geometric ideas in Markov chain Monte Carlo is reviewed in order to highlight these advances and their possible application in a range of domains beyond Statistics. A full exposition of Markov chains and their use in Monte Carlo simulation for Statistical inference and molecular dynamics is provided, with particular emphasis on methods based on Langevin diffusions. After this geometric concepts in Markov chain Monte Carlo are introduced. A full derivation of the Langevin diffusion on a Riemannian manifold is given, together with a discussion of appropriate Riemannian metric choice for different problems. A survey of applications is provided, and some open questions are discussed.


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 [1] 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 [2,3] were highlighted by Girolami & Calderhead [1], 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. [4]). We refer interested readers to [5,6].
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 Metropolisadjusted 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 π(·) will refer to an n-dimensional probability distribution, and π(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 π(·) defined on a measurable space (X , B) whose density is only known up to some proportionality constant. Although the ith sample is dependent on the (i − 1)th, the Ergodic Theorem ensures that for an appropriately constructed Markov chain with invariant distribution π(·), long-run averages are consistent estimators for expectations under π(·). As a result, MCMC methods have proven useful in Bayesian Statistical inference, where often the posterior density π(x|y) ∝ f (y|x)π 0 (x) for some parameter x is only known up to a constant [7]. Here we briefly introduce some concepts from general state space Markov chain theory together with a short overview of MCMC methods. The exposition follows [8].

Markov Chain Preliminaries
A time-homogeneous Markov chain {X m } m∈N is a collection of random variables X m each of which is defined on a measurable space (X , B), such that: for any A ∈ B. We define the transition kernel P (x m−1 , A) = P[X m ∈ A|X m−1 = x m−1 ] for the chain to be a map for which P (x, ·) defines a distribution over (X , B) for any x ∈ X and P (·, A) is measurable for any A ∈ B. Intuitively, P defines a map from points to distributions in X . Similarly we define the m-step transition kernel to be We call a distribution π(·) invariant for {X m } m∈N if for all A ∈ B. If P (x, ·) admits a density p(x |x), this can be equivalently written: The connotation of (3) and (4) is that if X m ∼ π(·), then X m+s ∼ π(·) for any s ∈ N. In this instance we say the chain is 'at stationarity'. Of interest to us will be Markov chains for which there is a unique invariant distribution which is also the limiting distribution for the chain, meaning that for any x 0 ∈ X for which π(x 0 ) > 0 lim m→∞ P m (x 0 , A) = π(A) (5) for any A ∈ B. Certain conditions are required for (5) to hold, but for all Markov chains presented here these are satisfied (though see [8]). A useful condition which is sufficient (though not necessary) for π(·) to be an invariant distribution is reversibility, which can be shown by the relation π(x)p(x |x) = π(x )p(x|x ).
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 ∈ A and x i+1 ∈ B is equal to the probability that x i+1 ∈ A and x i ∈ B. The relation (6) will be the primary tool used to construct Markov chains with a desired invariant distribution in the next section.

Monte Carlo estimates from Markov chains
Of most interest here are estimators constructed from a Markov chain. The Ergodic Theorem states that for any chain {X m } m∈N satisfying (5) and any g ∈ L 1 (π) we have that with probability one [7]. This is a Markov chain analogue to the Law of Large numbers. The efficiency of estimators of the formt m = i g(X i )/m can be assessed through the autocorrelation between elements in the chain. We will assess the efficiency oft m relative to estimators t m = i g(Z i )/m, where {Z i } m∈N is a sequence of independent random variables each having distribution π(·). Provided Var π [g(Z i )] < ∞ then Var[t m ] = Var π [g(Z i )]/m. We now seek a similar result for estimators of the formt m .
It follows directly from the Kipnis-Varadhan Theorem [9] that an estimatort m from a reversible Markov chain for which X 0 ∼ π(·) satisfies: We will refer to the constant τ as the autocorrelation time for the chain .
Another commonly used measure of efficiency is the effective sample size m ef f = m/τ , which gives the number of independent samples from π(·) needed to give an equally efficient estimate for E π [g(X)]. Clearly minimising τ is equivalent to maximising m ef f . 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 , ·) is from π(·), which requires a notion of distance between probability distributions.
Although there are several appropriate choices [11], 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 B according to µ(·) and ν(·). If both distributions admit densities, (9) can be written which is proportional to the L 1 distance between µ(x) and ν(x). Our metric · T V ∈ [0, 1], with · T V = 1 for distributions with disjoint supports and µ(·) − ν(·) T V = 0 implying µ(·) ≡ ν(·).
Typically for an unbounded X the distance P m (x 0 , ·) − π(·) T V will depend on x 0 for any finite m. So bounds on the distance are often sought via some inequality of the form where V : X → [1, ∞) depends on x 0 and is called a drift function, and f : N → [0, ∞) depends on the number of iterations m (and is often defined such that f (0) = 1). 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 π(·) 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 formt m , so that the shape of the distribution is asymptotically Gaussian. For more detail on this see [12,13].
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. [14]). 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.

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 π(·). 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 [7]. At iteration i, a sample is drawn from some candidate transition kernel Q(x i−1 , ·), 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 , ·) admits a density q(x |x i−1 ) for all x i−1 ∈ X (though see [8]). In this case a single step is shown below (the wedge notation a ∧ b denotes the minimum of a and b).

Algorithm 1 Metropolis-Hastings, single iteration
Require: The 'acceptance rate' α(x i−1 , x ) 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, α ≈ 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 for any A ∈ B, where is the average probability that a draw from Q(x, ·) will be rejected, and δ x (A) = 1 if x ∈ A and zero otherwise. A Markov chain defined in this way will have π(·) as an invariant distribution, since the chain is reversible for π(·). 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 π(·). It can be shown that π(·) is also the limiting distribution for the chain [7]. 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 , ·) and the target distribution π(·). 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.

Random walk proposals
An extremely simple choice for Q(x, ·) is one for which: where · denotes some appropriate norm on 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 π(·) will always be accepted. A typical choice for Q(x, ·) is N (x, λ 2 Σ), where the matrix Σ is often chosen in an attempt to match the correlation structure of π(·), or simply taken as the identity [15]. The tuning parameter λ 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 X tends to ∞ for a wide class of targets (e.g. [16,17]). 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. [18]), as no information about π(·) 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 π(·) is a one dimensional N (0, 1 2 ) distribution.  Several authors have also shown that for certain classes of π(·) the tuning parameter λ should be chosen such that λ 2 ∝ n −1 so that α 0 as n → ∞ [19]. Because of this we say that algorithm efficiency 'scales' O(n −1 ) as the dimension n of π(·) increases.
Ergodicity results for a Markov chain constructed using the RWM algorithm also exist [20,21,22]. At least exponentially light tails are a necessity for π(x) for geometric ergodicity, which means that π(x)/e − x → c as x → ∞, for some constant c. For super-exponential tails (where π(x) → 0 at a faster than exponential rate), additional conditions are required [20,22]. We demonstrate with a simple example why heavy-tailed forms of π(x) pose difficulties here (where π(x) → 0 at a rate slower then e − x ).
In the remainder of the article we will primarily discuss another approach to choosing Q which has been shown empirically [1] and in some cases theoretically [19] 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 [15].

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 continuoustime n-dimensional Markov processes (X t ) t≥0 for which any sample path t → X t (ω) 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 (X , B) as before. In the next section we provide some preliminaries, followed by an introduction to our main object of study, the Langevin diffusion.

Preliminaries
We focus on the class of time-homogeneous Itô diffusions, whose dynamics are governed by a stochastic differential equation of the form where (B t ) t≥0 is a standard Brownian motion and the drift vector b and volatility matrix σ are Lipschitz continuous [23].
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≥0 through the relation then we have giving the stochastic part of the relationship between X t+ t and X t for small enough t, see e.g. [24]. While (15) is often a suitable description of an Itô diffusion, it can also be characterised through an infinitessimal generator 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 ∈ C 0 (X ) as though A can be associated with the drift and volatility of (X t ) t≥0 by the relation where V ij (x) denotes the component in row i and column j of σ(x)σ(x) T [23]. 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 , ·). In the case of an Itô diffusion, P t (x 0 , ·) 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 , ·) is unknown the expectation and variance of X t ∼ P t (x 0 , ·) are given by the integral equations: where the second of these is a result of the Itô isometry [23]. Continuing the analogy, a natural question is whether a diffusion process has an invariant distribution π(·), and whether for any A ∈ B and any x 0 ∈ X , in some sense. For a large class of diffusions (which we confine ourselves to) this is in fact the case [25], and in addition (21) provides a means of finding π(·) given b and σ.
Setting the left-hand side of (21) to zero gives which can be solved to find π(·).

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 π(·). 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 [26], 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 π(x) is suitably regular so that ∇ log π(x) is well-defined and the derivatives in (23) exist, we can use (24) to construct a diffusion which has invariant distribution π(·).
Roberts & Tweedie [27] give sufficient conditions on π(·) under which a diffusion (X t ) t≥0 with dynamics given by (24) will be ergodic, meaning as t → ∞, for any x 0 ∈ X . They remark that these conditions 'should be appropriate for virtually all commonly encountered target densities' [27].

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, ·) 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 λ 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 π(x), combined with some random additive Gaussian noise, with variance λ 2 for each component. The relative weights of the deterministic and random parts are fixed, given as they are by the parameter λ. Typically if λ 1/2 λ then the random part of the proposal will dominate, and vice versa in the opposite case, though this also depends on the form of ∇ log π(x) [27].
Again since this is a Metropolis-Hastings method, choosing λ 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 → ∞ the optimal acceptance rate for the algorithm is 0.574 [19] for forms of π(·) which either have independent and identically distributed components or whose components only differ by some scaling factor [19]. In these cases, as n → ∞ the parameter λ must be ∝ 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 [27] highlight some cases in which MALA is either geometrically ergodic or not. Typically results are based on the tail behaviour of π(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 π(x) ∝ 1/(1 + x 2 ) as in the previous example. Then ∇ log π(x) = −2x/(1 + x 2 ) 2 → 0 as |x| → ∞. 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.
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 π(x) [28]. In this case, proposals take the form where λ is again a tuning parameter. It can be shown that provided Σ is a constant matrix, π(x) is still the invariant distribution for the diffusion on which (28) is based [29].

Geometric concepts in Markov Chain Monte Carlo
Ideas from Information geometry have been successfully applied to Statistics from as early as [30]. More widely, other geometric ideas have also been applied, offering new insight into common problems (e.g. [31,32]). A survey is given in [33]. 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.

Manifolds and Markov chains
We often make assumptions in MCMC about the properties of the space X in which our Markov chains evolve. Often X = R n , and in many cases where it is not then a simple re-parametrisation would make it so. But here by R n we mean the set of ordered n-tuples of real numbers, R n = {(a 1 , ..., a n ) : a i ∈ (−∞, ∞) ∀i}. The additional assumption that is often made is that R n is Euclidean, an inner product space with the induced distance metric 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 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 X which is 'likely to occur' under π(·), 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.
While the idea is attractive, it is far from a constructive definition. We only have the pre-requisite that (X , d) must be a metric space. However, as Langevin dynamics use gradient information, we will require (X , d) to be a space on which we can do differential calculus. Riemannian manifolds are an appropriate choice therefore, as the rules of differentiation are well understand for functions defined on them [34,35], while we are still free to define a more local notion of distance than Euclidean. In this section we write R n to denote the Euclidean vector space.

Preliminaries
We do not provide a full overview of Riemannian geometry here [34,35,36]. We simply note that for our purposes we can consider an n-dimensional Riemannian manifold (henceforth manifold) to be an n-dimensional metric space, in which distances are defined in a specific way. We also only consider manifolds for which a global coordinate chart exists, meaning that a mapping r : R n → M exists which is both differentiable and invertible, and for which the inverse is also differentiable (a diffeomorphism). Although this restricts the class of manifolds available (the sphere, for example, is not in this class), it is again suitable for our needs, and avoids the practical challenges of switching between coordinate patches. The connection with R n defined through r is crucial for making sense of differentiability in M . We say a function f : As has been stated, (29) can be induced via a Euclidean inner product, which we denote ·, · . However, it will aid intuition to think of distances in R n via curves We could think of the distance between two points in x, y ∈ R n as the minimum length among all curves which pass through x and y. If γ(0) = x and γ(1) = y the length is defined as giving the metric In R n the curve with minimum length will be a straight line, so that (32) agrees with (29). More generally we call a solution to (32) a geodesic [34].
In a vector space, metric properties can always be induced through an inner product (which also gives a notion of orthogonality). Such a space can be thought of as 'flat', since for any two points y and z, the straight line ay + (1 − a)z, a ∈ [0, 1] is also contained in the space. In general, manifolds do not have vector space structure globally, but do so at the infinitessimal level. As such we can think of them as 'curved'. We cannot always define an inner product, but we can still define distances through (32). We define a curve on a manifold M as γ M : [0, 1] → M . At each point γ M (t) = p ∈ M the velocity vector γ M (t) lies in an n-dimensional vector space which touches M at p. These are known as tangent spaces, denoted T p M , which can be thought as local linear approximations to M . We can define an inner product on each as g p : T p M → R, which allows us to define a generalisation of (31) as and provides a means to define a metric on the manifold as d(

Embeddings and local coordinates
Much of Riemannian geometry is written in a 'coordinate free' manner, and coordinate invariance is a pivotal concept in the subject [35]. For our purposes though, it suits calculations to work with a defined set of coordinates. We use the term 'local' coordinates here, since manifolds need only be locally mappable to R n . Although we restrict study to manifolds which are globally mappable to R n , we still use this convention. So far we have introduced manifolds as abstract objects. In fact, they can also be considered as objects which are embedded in some higher-dimensional Euclidean space. A simple example is any twodimensional surface, such as the unit sphere, lying in R 3 . If a manifold is embedded in this way, then metric properties can be induced from the ambient Euclidean space. The Nash embedding theorem [37] states that any manifold can be embedded in such a way. Figure 2: A two-dimensional manifold (surface) embedded in R 3 through r(x 1 , x 2 ) = (x 1 , x 2 , sin(x 1 )+1), parametrised by the local coordinates x 1 and x 2 . The distance between points A and B is given by the length of the curve γ(t) = (t, t, sin(t) + 1)).

A B
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 We can see that M is embedded in R 3 , but that any point can be identified using only two coordinates x 1 and x 2 . In this case, each T p M is a plane, and therefore a two-dimensional subspace of R 3 , so (i) it inherits the Euclidean inner product ·, · and (ii) any vector v ∈ T p M can be expressed as a linear combination of any two linearly independent basis vectors: a canonical choice is the partial derivatives ∂r/∂x 1 := r 1 and r 2 evaluated at x = r −1 (p) ∈ R 2 . The resulting inner product g p (v, w) between two vectors v, w ∈ T p M can be induced from the Euclidean inner product, as v, w = v 1 r 1 (x) + v 2 r 2 (x), w 1 r 1 (x) + w 2 r 2 (x) , and we use v i , w i to denote the components of v and w. To write (31) using this notation, we define the curve (31) can then be written which can be used in (32) as before.
The key point is that although we have started with an object embedded in R 3 , we can compute the Riemannian metric g p (v, w) (and hence distances in M ) using only the two-dimensional 'local' coordinates (x 1 , x 2 ). We also need not have explicit knowledge of the mapping r, only the components of the positive definite matrix G(x). The Nash embedding theorem in essence enables us to define manifolds by the reverse process: simply choose the matrix G(x) so that we define a metric space with suitable distance properties, and some object embedded in some higher-dimensional Euclidean space will exist for which these metric properties can be induced as above. So to define our new space, we simply choose an appropriate matrix-valued map G(x) (we discuss this choice in Subsection 4.4). If G(x) does not depend on x, then M has a vector space structure, and can be thought of as 'flat'. Trivially, G(x) = I gives Euclidean n-space.
We can also define volumes on a Riemannian manifold in local coordinates. Following standard coordinate transformation rules, we can see that for the above example the area element dx in R 2 will change according to a jacobian J = |(Dr) T (Dr)| 1/2 , where Dr = ∂(p 1 , p 2 , p 3 )/∂(x 1 , x 2 ). This reduces to J = |G(x)| 1/2 , which is also the case for more general manifolds [34]. We therefore define the Riemannian volume measure on a manifold M in local coordinates as If G(x) = I then this reduces to Lebesgue measure.

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 [38]. 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.
We turn first to (B t ) t≥0 , which we use to denote Brownian motion on a manifold. Intuitively, we may think of a construction based on embedded manifolds, by settingB 0 = p ∈ M , and for each increment sampling some random vector in the tangent space T p M , and then moving along the manifold in the prescribed direction for an infinitessimal period of time before re-sampling another velocity vector from the next tangent space [38]. In fact we can define such a construction using Stratonovitch calculus, and show that the infinitessimal generator can be written using only local coordinates [24]. Here we instead take the approach of generalising the generator directly from Euclidean space to the local coordinates of a manifold, arriving at the same result. We then deduce the stochastic differential equation describing (B t ) t≥0 in Itô form using (20).
For a standard Brownian motion on R n , A = /2, where denotes the Laplace operator Substituting A = /2 into (20) trivially gives b i (x) = 0 ∀i, V ij (x) = 1 {i=j} as required. The Laplacian f (x) is the divergence of the gradient vector field of some function f ∈ C 2 (R n ), and its value at x ∈ R n can be thought of as the average value of f in some neighbourhood of x [39].
To define a Brownian motion on any manifold, the gradient and divergence operators must be generalised. The gradient of a function on R n is the unique vector field such that for any unit vector u the directional derivative of f along u at x ∈ R n . On a manifold, the gradient operator ∇ M can still be defined such that the inner product g which is equal to the directional derivative along u as required.
The divergence of some vector field v at a point x ∈ R n is the net outward flow generated by v through some small neighbourhood of x. Mathematically, the divergence of v(x) ∈ R 3 is given by i ∂v i /∂x i . On a more general manifold the divergence is also a sum of derivatives, but here they are covariant derivatives. A short introduction is provided in Appendix A. Here we simply state that the covariant derivative of a vector field v at a point p ∈ M is the orthogonal projection of the directional derivative onto the tangent space T p M . Intuitively, a vector field on a manifold is a field of vectors each of which lie in the tangent space to a point p ∈ M . It only makes sense therefore to discuss how vector fields change along the manifold, or in the direction of vectors which also lie in the tangent space. Although the idea seems simple, the covariant derivative has some attractive geometric properties, notably it can be completely written in local coordinates, and so does not depend on knowledge of an embedding in some ambient space.
The divergence of a vector field v defined on a manifold M at the point p ∈ M is defined as: where e i denotes the ith basis vector for the tangent space T p M at p ∈ M , and v i denotes the ith 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. [40,41]), as: for some f ∈ C 2 0 (M ).
The generator of a Brownian motion on M is LB /2 [40]. Using (20) the resulting diffusion has dynamics given by Those familiar with the Itô formula will not be surprised by the additional drift term Ω(X t ). As Itô integrals do not follow the chain rule of ordinary calculus, non-linear mappings of martingales like (B t ) t≥0 typically result in drift terms being added to the dynamics (e.g. [23]). To define∇ we simply note that this is again the gradient operator on a general manifold, sõ ∇ = G −1 (x)∇. For the densityπ(x), we note that this density will now implicitly be defined with respect to the volume measure on the manifold. So to ensure the diffusion (39) has the correct invariant density with respect to Lebesgue measure, we definẽ Putting these three elements together, (39) becomes which, upon simplification, becomes It can be shown that this diffusion has invariant Lebesgue density π(x) as required [29]. Intuitively when a set is mapped onto the manifold, distances are changed by a factor G(x). So to end up with the initial distances, they must first be changed by a factor of 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 [29], and is given by where λ 2 is a tuning parameter. The nonlinear drift term here is slightly different to that reported in [1,28], for reasons discussed in [29].

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 π(x|y), as some of the discussion is directed towards Bayesian inference, where π(·) 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 [42] and others [2] for many years. Although generic measures of distance between distributions (such as total variation) are often appropriate, based on Informationtheoretic principles one can deduce that for a given parametric family {p x (y) : x ∈ 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 α = 0 connection employed, see [2] for details).
Because of this, Girolami & Calderhead [1] proposed a variant of the Fisher metric for geometric Markov chain Monte Carlo, as where π(x|y) ∝ f (y|x)π 0 (x) is the target density, f denotes the likelihood and π 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 [28], the objective is to choose G −1 (x) to match the covariance structure of π(x|y) locally. If the target density were Gaussian with covariance matrix Σ, then In the non-Gaussian case, the negative Hessian is no longer constant, but we can imagine that it matches the correlation structure of π(x|y) locally at least. Such ideas have been discussed in the geostatistics literature previously [43]. One problem with simply using (50) to define a metric is that unless π(x|y) is log-concave the negative Hessian will not be globally positive-definite, although Petra et al. [44] conjecture that it may be appropriate for use in some realistic scenarios, and suggest some computationally efficient approximation procedures [44].
, which is negative if x 2 > 1, so unusable as a proposal variance.
Girolami & Calderhead [1] 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. [45]), so provided a log-concave prior is chosen then (49) should be a suitable choice for G(x). Indeed, Girolami & Calderhead [1] provide several examples in which geometric MCMC methods using this Fisher metric perform better than their 'non-geometric' counterparts.
Betancourt [46] also starts from the viewpoint that the Hessian (50) is an appropriate choice for G(x), and defines a mapping from the set of n × n matrices to the set of positive-definite n × n matrices by taking a 'smooth' absolute value of the eigenvalues of the Hessian. This is done in a way such that derivatives of G(x) are still computable, inspiring the author to the name SoftAbs metric. For a fixed value of x, the negative Hessian H(x) is first computed, and then decomposed into U T DU , where D is the diagonal matrix of eigenvalues. Each diagonal element of D is then altered by the mapping t α : R → R, given by where α is a tuning parameter (typically chosen to be as large as possible for which eigenvalues remain non-zero numerically). The function t α acts as an absolute value function, but also uplifts eigenvalues which are close to zero to ≈ 1/α. 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 π(·). Many authors (e.g. [1,44]) 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 [1] 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 = |λ 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 × n matrices. The problem has in fact been well-studied in mathematical finance, in the context of finding correlations using incomplete data sets [47], and tackled using distances induced by the Frobenius norm. Approximate solution algorithms are discussed in Higham [47]. 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 π(x) ∝ 1/(1 + x 2 ), and set G(x) = | − ∂ 2 log π(x)/∂x 2 -. Then G −1 (x)∇ log π(x) = −x(1 + x 2 )/|1 − x 2 |, which no longer tends to 0 as |x| → ∞, 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| ≈ 1, but since this event occurs with probability 0 we do not see it as a major cause for concern.
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. [48] 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, ·) ≡ N (x + λ 2 Ω(x), λ 2 G −1 (x)). If we assume (like in the simplified manifold MALA) that Ω(x) ≈ 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 |x) = q(x|x ) 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 π(x|y) need not be known exactly. Examples from the article demonstrate the power of the approach [48].

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. [49] 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 [49].
Calderhead & Girolami [50] 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 [50]. 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 [51] 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 [51].
Konukoglu et al. [52] 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 [52]. 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 [52]. 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 radiofrequency ablation [52].

Discussion
The geometric viewpoint in not necessary to understand manifold variants of the MALA. Indeed, several authors [28,29] 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 [1]. 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) = Σ 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 [53]. 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 Λ(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 π(x) ∝ e −|x| , the Hessian is everywhere equal to zero! Much work in Information geometry has centred on the geometry of Hessian structures [54], 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 π(·) [18,55]. 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 [56]. The key challenge for MCMC is to define transition kernels for which proposed moves are inside the support for π(·). 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 [18]. However, the posterior may still look very different from the prior, as noted in [57], 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 π(·), 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 [58,59]. A natural extension to the work of Girolami & Calderhead [1] and Xifara et al. [29] 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 [59], though the most appropriate discretisation scheme for a second order process to facilitate sampling is unclear, and perhaps a question worthy of further exploration.
We have focused primarily here on the sample space X = R n , and on defining an appropriate manifold on which to construct Markov chains. In some inference problems, however, the sample space is a pre-defined manifold, for example the set of n × n rotation matrices, commonly found in the field of Directional Statistics [60]. Such manifolds are often not globally mappable to Euclidean n-space. Methods have been devised for sampling from such spaces [61,62]. In order to use the methods described here for such problems, an appropriate approach for switching between coordinate patches at the relevant time would need to be devised, which could be an interesting area of further study.
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. [63] have shown that in one dimension, the manifold MALA is geometrically ergodic for a class of targets of the form π(x) ∝ exp(−|x| β ) for any choice of β = 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 λ 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 π(·). 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 [17] notes that a requirement for optimal acceptance rate results for the RWM to be appropriate is that the curvature of π(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.

A Vector fields and the Covariant Derivative
Here we provide a short introduction to vector fields and differentiation on a smooth manifold, see [34,35]. The following geometric notation is used here: (i) vector components are indexed with a superscript, e.g. v = (v 1 , ..., v n ), and (ii) repeated subscript and superscripts are summed over, e.g. v i e i = i v i e i (known as the Einstein summation convention).
For any smooth manifold M , the set of all tangent vectors to points on M is known as the tangent bundle, and denoted T M .
A C r vector field defined on M is a mapping which assigns to each point p ∈ M a tangent vector v(p) ∈ T p M . In addition, the components of v(p) in any basis for T p M must also be C r [34]. We will denote the set of all vector fields on M as Γ(T M ). For some vector field v ∈ Γ(T M ), at any point p ∈ M the vector v(p) ∈ T p M can be written as a linear combination of some n basis vectors {e 1 , ..., e n } as v = v i e i . To understand how v will change in a particular direction along M , it only makes sense therefore to consider derivatives along vectors in T p M . Two other things must be considered when defining a derivative along a manifold: (i) how the components v i of each basis vector will change, and (ii) how each basis vector e i itself will change. For the usual directional derivative on R n the basis vectors do not change, as the tangent space is the same at each point, but for a more general manifold this is no longer the case: the e i 's are referred to as a 'local' basis for each T p M .
The covariant derivative D c is defined so as to account for these shortcomings. When considering differentiation along a vector u * / ∈ T p M , u * is simply projected onto the tangent space. The derivative with respect to any u ∈ T p M can now be decomposed into a linear combination of derivatives of basis vectors and vector components D c 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 [34], so (53) can be decomposed into The operator D c need therefore only be defined along the direction of basis vectors e i and for vector component v i and basis vector e i arguments. For components v i , D c e j [v i ] is defined as simply the partial derivative ∂ j v i := ∂v i /∂x j . The directional derivative of some basis vector e i along some e j is best understood through the example of a regular surface Σ ⊂ R 3 . Here D ej [e i ] will be a vector w ∈ R 3 . Taking the basis for this space at the point p as {e 1 , e 2 ,n}, wheren denotes the unit normal to T p Σ, we can write w = αe 1 + βe 2 + κn. The covariant derivative D c ej [e i ] is simply the projection of w onto T p Σ, given by w * = αe 1 + βe 2 . More generally at some point p in a smooth manifold M the covariant derivative D c ej [e i ] = Γ k ji e k (with upper and lower indices summed over). The coefficients Γ k ji are known as the Christoffel symbols: Γ k ji denotes the coefficient of the kth basis vector when taking the derivative of the ith with respect to the jth. If a Riemannian metric g is chosen for M , then they can be expressed completely as a function of g (or in local coordinates as a function of the matrix G). Using these definitions, (54) can be re-written as The divergence of a vector field v ∈ Γ(T M ) at the point p ∈ M is given by where again repeated indices are summed over. If M = R n , this reduces to the usual sum of partial derivatives ∂ i v i . On a more general manifold M , the equivalent expression is 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 [64] Γ j ij = |G(x)| − 1 2 ∂ i |G(x)| 1 2 , so (57) becomes where v = v(x).