2021

[1]James A. Brofos

1]Department of Statistics and Data Science, Yale University, 24 Hillhouse
Ave., New Haven, 06510, Connecticut, USA

2]Department of Statistics, Iowa State University, 2438 Osborn Drive, Ames,
50011, Iowa, USA

# Geometric Ergodicity in Modified Variations of Riemannian Manifold and
Lagrangian Monte Carlo

<EMAIL_ADDRESS>Vivekananda Roy<EMAIL_ADDRESS>Roy R. Lederman
<EMAIL_ADDRESS>[  [

###### Abstract

Riemannian manifold Hamiltonian (RMHMC) and Lagrangian Monte Carlo (LMC) have
emerged as powerful methods of Bayesian inference. Unlike Euclidean
Hamiltonian Monte Carlo (EHMC) and the Metropolis-adjusted Langevin algorithm
(MALA), the geometric ergodicity of these Riemannian algorithms has not been
extensively studied. On the other hand, the manifold Metropolis-adjusted
Langevin algorithm (MMALA) has recently been shown to exhibit geometric
ergodicity under certain conditions. This work investigates the mixture of the
LMC and RMHMC transition kernels with MMALA in order to equip the resulting
method with an “inherited” geometric ergodicity theory. We motivate this
mixture kernel based on an analogy between single-step HMC and MALA. We then
proceed to evaluate the original and modified transition kernels on several
benchmark Bayesian inference tasks.

###### keywords:

Markov chain Monte Carlo, Riemannian Manifold Hamiltonian Monte Carlo,
Statistical Computing, Lagrangian Monte Carlo

##  1 Introduction

Bayesian inference seeks to combine subjective sources of information with
observational data. By specifying one’s prior beliefs and correctly capturing
sources of uncertainty in a stochastic system, one may employ the Bayesian
approach in order to capture and reason about uncertainty under the posterior
distribution. Formally, the posterior distribution is a probability
distribution $\Pi:\mathfrak{B}(\mathbb{R}^{m})\to[0,1]$ with density
$\pi:\mathbb{R}^{m}\to\mathbb{R}_{+}$. As a practical matter in Bayesian
inference, one is concerned with the computation of expectations with respect
to $\Pi$ of integrable functions $h:\mathbb{R}^{m}\to\mathbb{R}$; that is, we
wish to compute $\underset{q\sim\Pi}{\mathbb{E}}[h(q)]$. Since $\pi$ is
generally intractable, these expectations are not available in closed form and
Monte Carlo methods based on samples from $\pi$ are often used to approximate
these means. For instance, one may compute the mean and variance of the
posterior by employing appropriate choices of expectation. Therefore, a
central problem in Bayesian inference is the generation of samples from the
posterior distribution.

One popular Monte Carlo method is based on approximate samples from the
posterior distribution generated via the technique of Markov chain Monte Carlo
(MCMC). In MCMC, given some initial state of the chain, $q^{0}$, a sequence of
random variables $(q^{1},q^{2},\ldots)$ is generated inductively: the
conditional distribution $q^{n}|q^{n-1},\ldots,q^{0}$ obeys the Markov
property. Under certain desirable regularity conditions, the random variables
$q^{n}$ will have a distribution which, asymptotically, converges to $\Pi$.
Establishing the rate of convergence, or at least an upper bound on that rate,
serves to quantify how quickly the chain $(q^{1},q^{2},\ldots)$ mixes toward
the target distribution. In the case where the rate of convergence is
geometrically fast, the Markov chain enjoys a strong stability property in
that the estimator $n^{-1}\sum_{i=1}^{n}h(q^{i})$ of
$\underset{q\sim\Pi}{\mathbb{E}}[h(q)]$ can be equipped with a central limit
theorem under additional technical assumptions. This central limit theorem has
important practical implications in that it allows the use of asymptotically
valid standard errors of MCMC estimates, which, in turn, can be used to decide
how long to run the Markov chains (Roy, 2020).

One of the most popular MCMC methods is Euclidean Hamiltonian Monte Carlo
(EHMC) (Neal, 2010b). In EHMC, one computes approximate solutions to
Hamilton’s equations of motion using numerical integrators, taking care to
ensure that the resulting Markov chain satisfies detailed balance in a phase-
space consisting of the position variables $q\in\mathbb{R}^{m}$ and auxiliary
momentum variables $p\in\mathbb{R}^{m}$. The rate of convergence of HMC was
established by Livingstone et al. (2016) in the presence of conditions on the
log-density (see also Durmus et al., 2020). Traditional EHMC, however,
struggles in distributions that exhibit multiple spatial scales (Pourzanjani
and Petzold, 2019; Betancourt, 2012). This observation led to the development
of geometric methods of HMC, specifically the Riemannian manifold Hamiltonian
Monte Carlo (RMHMC) method of Girolami and Calderhead (2011).

Unlike EHMC, RMHMC adapts to the second-order structure of the posterior,
which allows it to align its proposals in the direction of the posterior that
exhibits the greatest local variation. However, the sophisticated form of the
Hamiltonian employed in RMHMC necessitates the use of complex numerical
integrators that are significantly more expensive than the numerical
integrator employed in EHMC. This concern was partially alleviated in the
introduction of Lagrangian Monte Carlo (LMC) in Lan et al. (2015), which was
able to construct a MCMC method with a more efficient numerical integration
procedure, while simultaneously continuing to take advantage of second-order
geometric knowledge in order to efficiently traverse the posterior. Both RMHMC
and LMC can be viewed as generalizations of the Euclidean Hamiltonian Monte
Carlo algorithm, which incorporate second-order geometric information about
the posterior into the Markov chain transition kernel. Incorporating second-
order information allows these geometric methods of MCMC to explore the
typical region of the target distribution more efficiently.

However, neither RMHMC nor LMC have, as of yet, been equipped with a geometric
ergodicity theory. The focus in this work is in establishing the geometric
ergodicity of modified versions of RMHMC and LMC. In implementations of RMHMC
and LMC it is common to randomize the number of integration steps, and to take
one-step with positive probability. In the case of Euclidean HMC, the
equivalence between single-step HMC and the Metropolis-adjusted Langevin
algorithm (MALA) is critical for establishing geometric ergodicity of HMC in
Livingstone et al. (2016). Our approach is two-fold. First, we propose a
simple modification to both RMHMC and LMC which is motivated by a unique
correspondence between single-step Euclidean HMC and MALA. In particular, we
propose that instead of applying the transition kernel corresponding to RMHMC
or LMC with a single integration step, one instead applies the transition
kernel of the manifold Metropolis-adjusted Langevin algorithm (MMALA). Recent
work in Roy and Zhang (2023) gave conditions under which MMALA is
geometrically ergodic. We may imbue these modified variations of RMHMC and LMC
with a geometric ergodicity theory via inheritance once one establishes that
the RMHMC and LMC transition kernels (with a fixed, non-random number of
integration steps) are reversible in the required sense. This construction is
described in section 3. In section 4 we proceed to a numerical evaluation of
the proposed modification of RMHMC and LMC, with special attention given to
the probability of applying the MMALA transition kernel. We begin, however, in
section 2 with an overview of required mathematical concepts from numerical
integration and Markov chains.

##  2 Preliminaries

In section 2.1 we review the generalized leapfrog and Lagrangian leapfrog
employed in geometric MCMC methods for performing Bayesian inference. Section
2.2 covers Markov chains based on involutions as well as those based on
discretizations of Langevin diffusions. We also discuss geometric ergodicity,
mixture transition kernels, and the equivalence between single-step HMC and
MALA in the Euclidean regime.

###  2.1 Numerical Integrators

Hamiltonian and Lagrangian mechanics are equivalent, with the momentum and
velocity being related by the Legendre transform (Marsden and Ratiu, 2010). In
this section, we will consider numerical integrators of these equations of
motion. Throughout our discussion, we will consider a Hamiltonian function of
the following form.

######  Definition 1.

Let $\pi:\mathbb{R}^{m}\to\mathbb{R}_{+}$ be a smooth probability density and
let $\mathbf{G}:\mathbb{R}^{m}\to\mathrm{PD}(m)$. The Riemannian Hamiltonian
is the function $H:\mathbb{R}^{m}\times\mathbb{R}^{m}\to\mathbb{R}$ defined by

$\displaystyle
H(q,p)=-\log\pi(q)+\frac{1}{2}\log\mathrm{det}(\mathbf{G}(q))+\frac{1}{2}p^{\top}\mathbf{G}^{-1}(q)p.$
(1)

Moreover, assuming
$\int_{\mathbb{R}^{m}}\int_{\mathbb{R}^{m}}\exp(-H(q,p))~{}\mathrm{d}q~{}\mathrm{d}p<+\infty$,
the Riemannian density is the probability density
$\pi(q,p)\propto\exp(-H(q,p))$.

######  Remark 1.

The term
$\frac{1}{2}\log\mathrm{det}(\mathbf{G}(q))+\frac{1}{2}p^{\top}\mathbf{G}^{-1}(q)p$
appearing in the definition of the function $H$ in eq. 1 is chosen such that
the conditional density $\pi(q,p)$ satisfies

$\displaystyle\pi(p|q)=\mathrm{Normal}(0,\mathbf{G}(q)).$ (2)

It is easily seen that the $q$-marginal density of $\pi(q,p)$ obtained by
marginalizing out $p$ is $\pi(q)$.

Hamiltonian functions $H:\mathbb{R}^{m}\times\mathbb{R}^{m}\to\mathbb{R}$
produce Hamilton’s equations of motion which, by definition, are solutions to
the coupled differential equations

$\displaystyle\dot{q}_{t}$ $\displaystyle=\nabla_{p}H(q_{t},p_{t})$ (3)
$\displaystyle\dot{p}_{t}$ $\displaystyle=-\nabla_{q}H(q_{t},p_{t}).$ (4)

Except in special cases, a closed-form solution for the map
$t\mapsto(q_{t},p_{t})\in\mathbb{R}^{m}\times\mathbb{R}^{m}$ does not exist;
this necessitates the use of numerical integrators in order to approximate the
equations of motion obeying eqs. 3 and 4. One such example is the generalized
leapfrog method, which we now define.

######  Definition 2.

The generalized leapfrog integrator with step-size $\epsilon\in\mathbb{R}$
applied to the Riemannian Hamiltonian in eq. 1 is the map
$(q,p)\mapsto\hat{\Phi}_{\epsilon}(q,p)$ defined by,

$\displaystyle\breve{p}$
$\displaystyle=p-\frac{\epsilon}{2}\nabla_{q}H(q,\breve{p})$ (5)
$\displaystyle\tilde{q}$
$\displaystyle=q+\frac{\epsilon}{2}\left({\nabla_{p}H(q,\breve{p})+\nabla_{p}H(\tilde{q},\breve{p})}\right)$
(6)
$\displaystyle=q+\frac{\epsilon}{2}\left({\mathbf{G}^{-1}(q)+\mathbf{G}^{-1}(\tilde{q})}\right)\breve{p}$
(7) $\displaystyle\tilde{p}$
$\displaystyle=\breve{p}-\frac{\epsilon}{2}\nabla_{q}H(\tilde{q},\breve{p})$
(8)
$\displaystyle=p-\frac{\epsilon}{2}\nabla_{q}H(q,\breve{p})-\frac{\epsilon}{2}\nabla_{q}H(\tilde{q},\breve{p}),$
(9)

where $\hat{\Phi}_{\epsilon}(q,p)=(\tilde{q},\tilde{p})$.

######  Remark 2.

One can prove (Leimkuhler and Reich, 2005; Hairer et al., 2006) that the
generalized leapfrog integrator is a symplectic transformation and that
therefore $|{\mathrm{det}(\nabla\Phi(q,p))}|=1$.

######  Remark 3.

We must now give a remark about notation. Numerical integration plays a
central role in our analysis, and we give particular attention to multiple
steps of integration. At the same time, we will also discuss Markov chains,
which also consist of multiple steps. As an attempt to differentiate these two
notions of step, we will use a lower index to refer to steps of numerical
integration, whereas we will use upper indices to denote Markov chain steps.
Given initial data $(q_{0},p_{0})\in\mathbb{R}^{m}\times\mathbb{R}^{m}$, we
denote $k$-steps of the generalized leapfrog integrator by
$(q_{k},p_{k})=\hat{\Phi}^{k}_{\epsilon}(q_{0},p_{0})=\hat{\Phi}_{\epsilon}(q_{k-1},p_{k-1})$.
Similarly, we let
$\breve{p}_{k}=p_{k-1}-\frac{\epsilon}{2}\nabla_{q}H(q_{k-1},\breve{p}_{k})$,
which we call the intermediate momentum at the $k$-th step.

A Hamiltonian of the form in eq. 1 may be transformed into a Lagrangian as

$\displaystyle
L(q,\dot{q})=-\log\pi(q)+\frac{1}{2}\log\mathrm{det}(\mathbf{G}(q))-\frac{1}{2}\dot{q}^{\top}\mathbf{G}(q)\dot{q}.$
(10)

As noted at the beginning of this subsection, Lagrangian and Hamiltonian
mechanics are formally equivalent, with the momentum $p\in\mathbb{R}^{m}$
being related to the velocity $\dot{q}\in\mathbb{R}^{m}$ according to the
Legendre transformation $p=\mathbf{G}(q)\dot{q}$. The Lagrangian produces
equivalent equations of motion called the Euler-Lagrange equations as

$\displaystyle\nabla_{q}L(q_{t},\dot{q}_{t})=\frac{\mathrm{d}}{\mathrm{d}t}\nabla_{\dot{q}}L(q_{t},\dot{q}_{t}).$
(11)

One advantage possessed by the Lagrangian formalism over the Hamiltonian
approach is that one may identify explicit numerical integrators (Lan et al.,
2015) of the equations of motion, such as the Lagrangian leapfrog.

######  Definition 3.

The Lagrangian leapfrog integrator with step-size $\epsilon\in\mathbb{R}$
applied to the Riemannian Hamiltonian in eq. 1 is the map
$(q,p)\mapsto\tilde{\Phi}_{\epsilon}(q,p)$ defined by,

$\displaystyle v$ $\displaystyle=\mathbf{G}^{-1}(q)p$ (12)
$\displaystyle\breve{v}$
$\displaystyle=v-\frac{\epsilon}{2}\mathbf{G}^{-1}(q)\nabla_{q}U(q)-\frac{\epsilon}{2}\Omega(q,v)\breve{v}$
(13) $\displaystyle\tilde{q}$ $\displaystyle=q+\epsilon\breve{v}$ (14)
$\displaystyle\tilde{v}$
$\displaystyle=\breve{v}-\frac{\epsilon}{2}\mathbf{G}^{-1}(\tilde{q})\nabla_{q}U(\tilde{q})-\frac{\epsilon}{2}\Omega(\tilde{q},\breve{v})\tilde{v}$
(15)
$\displaystyle=v-\frac{\epsilon}{2}\mathbf{G}^{-1}(q)\nabla_{q}U(q)-\frac{\epsilon}{2}\Omega(q,v)\breve{v}-\frac{\epsilon}{2}\mathbf{G}^{-1}(\tilde{q})\nabla_{q}U(\tilde{q})-\frac{\epsilon}{2}\Omega(\tilde{q},\breve{v})\tilde{v}$
(16) $\displaystyle\tilde{p}$
$\displaystyle=\mathbf{G}^{-1}(\tilde{q})\tilde{v}$ (17)

where $(\tilde{q},\tilde{p})=\tilde{\Phi}_{\epsilon}(q,p)$ and

$\displaystyle U(q)$
$\displaystyle=-\log\pi(q)+\frac{1}{2}\log\mathrm{det}(\mathbf{G}(q))$ (18)
$\displaystyle\Omega_{ij}(q,v)$
$\displaystyle=\sum_{k=1}^{m}\Gamma^{i}_{kj}(q)v_{k}$ (19)
$\displaystyle\Gamma^{i}_{kj}(q)$
$\displaystyle=\frac{1}{2}\sum_{l=1}^{m}\mathbf{G}^{-1}_{kl}(q)\left({\frac{\partial}{\partial
q_{i}}\mathbf{G}_{lj}(q)+\frac{\partial}{\partial
q_{j}}\mathbf{G}_{il}(q)-\frac{\partial}{\partial
q_{l}}\mathbf{G}_{ij}(q)}\right).$ (20)

The Jacobian determinant of the transformation $\tilde{\Phi}_{\epsilon}$
(defined in definition 3) is,

$\displaystyle|{\mathrm{det}(\nabla\tilde{\Phi}_{\epsilon}(q,p))}|=\mathrm{det}(\mathbf{G}^{-1}(q))\mathrm{det}(\mathbf{G}(\tilde{q}))\frac{|{\mathrm{det}(\mathrm{Id}+\epsilon\Omega(q,v)/2)\mathrm{det}(\mathrm{Id}-\epsilon\Omega(\tilde{q},\tilde{v})/2)}|}{|{\mathrm{det}(\mathrm{Id}+\epsilon\Omega(\tilde{q},\bar{v})/2)\mathrm{det}(\mathrm{Id}-\epsilon\Omega(q,\bar{v})/2)}|}.$
(21)

######  Remark 4.

Given initial data $(q_{0},p_{0})\in\mathbb{R}^{m}\times\mathbb{R}^{m}$, we
denote $k$-steps of the generalized (or Lagrangian) leapfrog integrator by
$(q_{k},p_{k})=\tilde{\Phi}^{k}_{\epsilon}(q_{0},p_{0})=\tilde{\Phi}_{\epsilon}(q_{k-1},p_{k-1})$.

In the case where $\mathbf{G}$ is a constant function of $q$, the generalized
leapfrog integrator reduces to the standard leapfrog integrator, which is
defined as follows.

######  Definition 4.

Let $\pi:\mathbb{R}^{m}\to\mathbb{R}_{+}$ be a smooth probability density. A
Euclidean Hamiltonian is a smooth map
$H:\mathbb{R}^{m}\times\mathbb{R}^{m}\to\mathbb{R}$ of the form,

$\displaystyle H(q,p)=-\log\pi(q)+\frac{1}{2}p^{\top}\mathbf{G}^{-1}p,$ (22)

where $\mathbf{G}\in\mathrm{PD}(m)$.

######  Definition 5.

Consider a Euclidean Hamiltonian as defined in definition 4. The Euclidean
leapfrog integrator is the map $(q,p)\mapsto\check{\Phi}_{\epsilon}(q,p)$
defined by

$\displaystyle\breve{p}$ $\displaystyle=p+\frac{\epsilon}{2}\nabla\log\pi(q)$
(23) $\displaystyle\tilde{q}$
$\displaystyle=q+\epsilon\mathbf{G}^{-1}\breve{p}$ (24)
$\displaystyle\tilde{p}$
$\displaystyle=\breve{p}+\frac{\epsilon}{2}\nabla\log\pi(\tilde{q}),$ (25)

where $\check{\Phi}_{\epsilon}(q,p)=(\tilde{q},\tilde{p})$.

The Euclidean leapfrog integrator is the de-facto standard numerical
integrator employed in EHMC. This is because of its accuracy and computational
efficiency; in contrast to the generalized leapfrog integrator in definition
2, the Euclidean leapfrog is fully explicit. Nevertheless, the generalized
leapfrog and Euclidean leapfrog are exactly equivalent when applied to a
Hamiltonian in the form of definition 4; this is made precise in the following
result.

######  Proposition 1.

When the generalized leapfrog (definition 2) or Lagrangian leapfrog
(definition 3) is applied to a Euclidean Hamiltonian in the form of definition
4 (in the sense that we take $\mathbf{G}(q)\equiv\mathbf{G}$ for every
$q\in\mathbb{R}^{m}$), the resulting map is equivalent to the Euclidean
leapfrog method in definition 5.

Proof: In the case of the generalized leapfrog integrator, this is immediate
from the fact that $\nabla_{q}H(q,p)=-\nabla\log\pi(q)$ which does not depend
on $p$ and $\nabla_{p}H(q,p)=\mathbf{G}p$ which does not depend on $q$.

In the case of the Lagrangian leapfrog, the form of the Euclidean Hamiltonian
means that $\Omega(q,v)$ is uniformly zero (because all of the derivatives of
the metric vanish). The Lagrangian leapfrog integrator may be re-expressed as

$\displaystyle v$ $\displaystyle=\mathbf{G}^{-1}p$ (26)
$\displaystyle\breve{v}$
$\displaystyle=v+\frac{\epsilon}{2}\mathbf{G}^{-1}\nabla\log\pi(q)$ (27)
$\displaystyle=\mathbf{G}^{-1}\breve{p}$ (28) $\displaystyle\tilde{q}$
$\displaystyle=q+\epsilon\breve{v}$ (29)
$\displaystyle=q+\epsilon\mathbf{G}^{-1}\breve{p}$ (30)
$\displaystyle\tilde{v}$
$\displaystyle=\breve{v}+\frac{\epsilon}{2}\mathbf{G}^{-1}\nabla\log\pi(\tilde{q})$
(31) $\displaystyle=\mathbf{G}^{-1}\tilde{p},$ (32)

where $\breve{p}$ and $\tilde{p}$ are as in eqs. 23 and 25.

###  2.2 Markov Chains

The purpose of this section is to review fundamentals of Markov chains,
including their construction, convergence properties, and central limit
theorems. Markov chains are stochastic processes that may be inductively
defined by their transition kernel. Since the transition kernel yields a
probability measure depending only on the current state, one sees by
inspection that the Markov chain satisfies the Markov property.

######  Definition 6.

A Markov chain transition kernel on $\mathbb{R}^{m}$ is a function
$K:\mathbb{R}^{m}\times\mathfrak{B}(\mathbb{R}^{m})\to[0,1]$ such that (i) for
$q\in\mathbb{R}^{m}$, $K(q,\cdot)$ is a probability measure and (ii) for fixed
$A\in\mathfrak{B}(\mathbb{R}^{m})$ the function $q\mapsto K(q,A)$ is
measurable.

Within the context of simulation-based inference and mainly Bayesian
inference, one is interested in constructing a Markov chain which converges to
a specified target probability distribution; as noted in the introduction,
this target distribution is typically known by its density up to a normalizing
constant. We now define two notions of convergence.

######  Definition 7.

Given a Markov chain transition kernel $K$, a Markov chain is a sequence of
random variables defined inductively by $q^{n+1}|q^{n}\sim K(q^{n},\cdot)$.

######  Definition 8.

The total variation distance between two probability measures
$\Pi:\mathfrak{B}(\mathbb{R}^{m})\to[0,1]$ and
$\Xi:\mathfrak{B}(\mathbb{R}^{m})\to[0,1]$ is defined by

$\displaystyle\|\Pi-\Xi\|_{\mathrm{TV}}=2\sup_{A\in\mathfrak{B}(\mathbb{R}^{m})}|{\Pi(A)-\Xi(A)}|.$
(33)

######  Definition 9.

A Markov chain with transition kernel $K$ is said to be ergodic if

$\displaystyle\lim_{n\to\infty}\|K^{n}(q,\cdot)-\Pi(\cdot)\|_{\mathrm{TV}}=0,$
(34)

where for $A\in\mathfrak{B}(\mathbb{R}^{m})$, $K^{n}(q,A)=\mathrm{Pr}(q^{n}\in
A|q^{0}=q)$ is the $n$-step Markov transition probability.

It will turn out that one can establish ergodicity of a Markov chain if one
can establish three separate properties: irreducibility, aperiodicity, and
stationarity. We now define each of these concepts.

######  Definition 10.

A Markov chain transition kernel $K$ is said to be $\Pi$-irreducible if for
every set $A\in\mathfrak{B}(\mathbb{R}^{m})$ with $\Pi(A)>0$ there exists
$n\in\mathbb{N}$ for which $K^{n}(q,A)>0$ for all $q\in\mathbb{R}^{m}$.

######  Definition 11.

Given a Markov chain transition kernel
$K:\mathbb{R}^{m}\times\mathfrak{B}(\mathbb{R}^{m})\to[0,1]$, a set
$C\subseteq\mathbb{R}^{m}$ is called small if there exists a $n\in\mathbb{N}$,
a $\delta>0$, and a probability measure
$\nu:\mathfrak{B}(\mathbb{R}^{m})\to[0,1]$ such that for any $q\in C$ and
$A\in\mathfrak{B}(\mathbb{R}^{m})$ we have $K^{n}(q,A)\geq\delta\cdot\nu(A)$.

We may therefore say that the set $C$ is $(n,\delta,\nu)$-small.

######  Definition 12.

A Markov chain transition kernel $K$ is called aperiodic if there exists a
small set $C$ for which the greatest common divisor of the set

$\displaystyle\left\\{{n^{\prime}\in\mathbb{N}:\exists~{}\delta>0~{}\text{and
probability measure}~{}\nu~{}\text{such
that}~{}C~{}\text{is}~{}(n^{\prime},\delta,\nu)\text{-small}}\right\\}$ (35)

is one.

######  Definition 13.

A Markov chain with transition kernel $K$ is said to be stationary for the
probability measure $\Pi$ if for any Borel set
$A\in\mathfrak{B}(\mathbb{R}^{m})$ we have
$\underset{q\sim\Pi}{\mathbb{E}}K(q,A)=\Pi(A)$.

######  Definition 14.

A Markov chain with transition kernel $K$ is called geometrically ergodic if
there exists $\rho\in(0,1)$ and function $V:\mathbb{R}^{m}\to\mathbb{R}_{+}$
such that

$\displaystyle\|K^{n}(q,\cdot)-\Pi(\cdot)\|_{\mathrm{TV}}\leq V(q)\rho^{n}.$
(36)

It is clear from the definitions that geometric ergodicity is a stronger form
of convergence than mere ergodicity. In the latter case, however, simple
conditions under which a Markov chain is ergodic can be provided.

######  Theorem 2 (Tierney (1994)).

Suppose that $K$ is a Markov chain transition kernel that is
$\Pi$-irreducible, aperiodic, and for which $\Pi$ is the stationary
probability measure. Then $K$ produces an ergodic Markov chain.

However, in addition to giving an upper bound on the mixing of the Markov
chain into the target probability measure, geometric ergodicity is also useful
for establishing a central limit theorem for expectations computed from the
sequence $(q^{1},q^{2},\ldots)$. We make this precise as follows.

######  Theorem 3 (Meyn and Tweedie (1993)).

Let $h:\mathbb{R}^{m}\to\mathbb{R}$ be a function for which
$\underset{q\sim\Pi}{\mathbb{E}}|{h(q)}|^{2+\gamma}<\infty$ for some
$\gamma>0$. Let $K$ be a Markov chain transition kernel that is aperiodic,
$\Pi$-irreducible, and for which $\Pi$ is the unique stationary distribution.
Assume further that $K$ converges geometrically to $\Pi$. Define
$\mathcal{S}_{K}[h]=n^{-1}\sum_{i=1}^{n}h(q^{i})$ where $q^{i+1}|q^{i}\sim
K(q^{i},\cdot)$. Then, as $n\to\infty$,

$\displaystyle\sqrt{n}\left({\mathcal{S}_{K}[h]-\underset{q\sim\Pi}{\mathbb{E}}\left[h(q)\right]}\right)\overset{d}{\to}\mathrm{Normal}(0,\tau^{2}_{h})$
(37)

where $\tau^{2}_{h}$ is a constant, less than infinity, that depends on $h$
(and $K$). The quantity $\tau^{2}_{h}$ has a closed-form given by

$\displaystyle\tau^{2}_{h}=\underset{q\sim\Pi}{\mathrm{Var}}\left(h(q)\right)+2\sum_{k=1}^{\infty}\underset{q\sim\Pi}{\mathrm{Cov}}\left(h(q),h(q^{k})\right).$
(38)

Establishing that a Markov chain has $\Pi$ as a stationary distribution is
often easily achieved by showing that the chain satisfies detailed balance
with respect to $\Pi$. As we require detailed balance when discussing the
marginal transition kernels of RMHMC and LMC in section 3, we define this
notion now.

######  Definition 15.

A Markov chain transition kernel is said to satisfy detailed balance with
respect to the probability distribution $\Pi$ (equivalently, the Markov chain
is called reversible with respect to $\Pi$) if for any sets
$Q,Q^{\prime}\in\mathfrak{B}(\mathbb{R}^{m})$,

$\displaystyle\int_{Q}K(q,Q^{\prime})~{}\Pi(\mathrm{d}q)=\int_{Q^{\prime}}K(q,Q)~{}\Pi(\mathrm{d}q).$
(39)

######  Remark 5.

When detailed balance holds for a Markov chain transition kernel $K$, it
follows that $\Pi$ is the stationary probability measure (definition 13) of
the Markov transition kernel $K$.

####  2.2.1 Metropolis-Hastings Kernels and Generalized Langevin Algorithms

In this section we review Metropolis-Hastings Markov chains. We begin by
formally defining these objects before proceeding to a special case of
Metropolis-Hastings method based on Langevin diffusion, which recalls the
constructions considered in Roy and Zhang (2023).

######  Definition 16 (Metropolis-Hastings Algorithm).

Let $\pi:\mathbb{R}^{m}\to\mathbb{R}_{+}$ be a probability density and, for
each $q\in\mathbb{R}^{m}$, let
$\tilde{\pi}(\cdot|q):\mathbb{R}^{m}\to\mathbb{R}_{+}$ be a probability
density depending on $q$. The Metropolis-Hastings transition kernel with
proposal density $\tilde{\pi}(\cdot|q)$ is,

$\displaystyle\begin{split}K(q,A)&=\int_{A}\min\left\\{{1,\frac{\pi(\tilde{q})\tilde{\pi}(q|\tilde{q})}{\pi(q)\tilde{\pi}(\tilde{q}|q)}}\right\\}\tilde{\pi}(\tilde{q}|q)~{}\mathrm{d}\tilde{q}\\\
&\qquad+~{}\left({1-\int_{\mathbb{R}^{m}}\min\left\\{{1,\frac{\pi(\tilde{q})\tilde{\pi}(q|\tilde{q})}{\pi(q)\tilde{\pi}(\tilde{q}|q)}}\right\\}\tilde{\pi}(\tilde{q}|q)~{}\mathrm{d}\tilde{q}}\right)\mathrm{1}\left\\{{q\in
A}\right\\},\end{split}$ (40)

for $q\in\mathbb{R}^{m}$ and $A\in\mathfrak{B}(\mathbb{R}^{m})$.

######  Definition 17 (Generalized Metropolis-Adjusted Langevin Algorithm).

Fix $\epsilon\in\mathbb{R}\setminus\left\\{{0}\right\\}$. Let
$c_{\epsilon}:\mathbb{R}^{m}\to\mathbb{R}^{m}$ and let
$\mathbf{A}:\mathbb{R}^{m}\to\mathrm{PD}(m)$. The generalized Metropolis-
adjusted Langevin algorithm is an instance of the Metropolis-Hastings
algorithm with proposal density

$\displaystyle\tilde{\pi}(\tilde{q}|q)=\mathrm{Normal}(\tilde{q}\mathchar
24635\relax\;c_{\epsilon}(q),\epsilon^{2}\mathbf{A}(q)).$ (41)

We denote the Markov chain transition kernel of the generalized Metropolis-
adjusted Langevin algorithm with step-size $\epsilon$ by
$J_{\epsilon}:\mathbb{R}^{m}\times\mathfrak{B}(\mathbb{R}^{m})\to[0,1]$.

######  Definition 18 (Riemannian Manifold Metropolis-Adjusted Langevin
Algorithm).

In the special case where

$\displaystyle c_{\epsilon}(q)$
$\displaystyle=\frac{\epsilon^{2}}{2}\mathbf{A}(q)\nabla\log\pi(q)+\frac{\epsilon^{2}}{2}\Gamma(q)$
(42) $\displaystyle\Gamma_{i}(q)$
$\displaystyle=\sum_{j=1}^{m}\frac{\partial}{\partial
q_{j}}\mathbf{A}_{ij}(q)$ (43)

we call the resulting method the Riemannian manifold Metropolis-adjusted
Langevin algorithm (MMALA).

######  Remark 6.

The form of the proposal in eq. 42 is based on an Euler-Maruyama
discretization of a Langevin diffusion on a manifold with Riemannian metric
$\mathbf{G}(q)=\mathbf{A}^{-1}(q)$; see Xifara et al. (2014) for details. In
summary, eq. 42 is the drift component of the Euler-Maruyama discretization
(with step-size $\epsilon^{2}$) applied to the following stochastic
differential equation:

$\displaystyle\mathrm{d}X_{t}=\frac{1}{2}\mathbf{G}^{-1}(X_{t})\nabla\log\pi(X_{t})~{}\mathrm{d}t+\Gamma(X_{t})~{}\mathrm{d}t+\mathbf{G}^{-1/2}(X_{t})~{}\mathrm{d}B_{t},$
(44)

where $B_{t}$ is Euclidean Brownian motion at time $t\in\mathbb{R}$. When
$\mathbf{G}(q)$ is the Fisher information matrix, the term
$\mathbf{G}^{-1}(X_{t})\nabla\log\pi(X_{t})$ can be identified as the natural
gradient of the function $\log\pi$ under the geometry generated by the Fisher
metric (Amari and Nagaoka, 2000) at position $X_{t}$. On the other hand, the
term $\Gamma(X_{t})~{}\mathrm{d}t+\mathbf{G}^{-1/2}(X_{t})~{}\mathrm{d}B_{t}$
corresponds to manifold Brownian motion, since its infinitesimal generator is
the Laplace-Beltrami operator on the manifold (Hsu, 2002). Therefore, as the
stochastic differential equation in eq. 44 consists of a term comprising the
manifold gradient of a log-density and another term comprising manifold
Brownian motion, it is called the Riemannian Langevin equation in
correspondence with the Euclidean case.

######  Definition 19 (Simplified Riemannian Manifold Metropolis-Adjusted
Langevin Algorithm).

Computing the function $\Gamma$ in definition 18 may be inconvenient to
evaluate. It can be ignored while still yielding a Markov chain transition
kernel, giving the special case,

$\displaystyle
c_{\epsilon}(q)=\frac{\epsilon^{2}}{2}\mathbf{A}(q)\nabla\log\pi(q).$ (45)

The resulting method is called the simplified Riemannian manifold Metropolis-
adjusted Langevin algorithm, which was considered by Girolami and Calderhead
(2011).

We adopt the abbreviation SMALA to refer to the simplified Metropolis-adjusted
Langevin algorithm.

######  Remark 7.

There are three common choices of the function $\mathbf{A}$. The first is that
$\mathbf{A}$ is a constant function, in which case we simply write
$\mathbf{A}\in\mathrm{PD}(m)$. A second option is that $\mathbf{A}$ is chosen
as the inverse of the sum of the Fisher information matrix and the negative
Hessian of the log-prior, which can capture second-order geometry of both the
likelihood and the prior; the use of the inverse of the sum of the Fisher
information and the negative Hessian of the log-prior as a preconditioner is
the approach advocated by Girolami and Calderhead (2011). A third option is
that $\mathbf{A}$ is the inverse of the SoftAbs metric, which is a smooth
transformation of the Hessian of the log-density of the target distribution;
for details see Betancourt (2012).

####  2.2.2 Geometric Methods of Bayesian Inference

######  Definition 20 (Involutive Monte Carlo (Neklyudov et al., 2020)).

Let $\Phi:\mathbb{R}^{m}\to\mathbb{R}^{m}$ be a smooth involution and let
$\pi:\mathbb{R}^{m}\to\mathbb{R}_{+}$ be a probability density on
$\mathbb{R}^{m}$. The Markov chain transition kernel of involutive Monte Carlo
with target density $\pi$ is

$\displaystyle\begin{split}K(q,A)&=\min\left\\{{1,\frac{\pi(\Phi(q))}{\pi(q)}|{\mathrm{det}(\nabla\Phi(q))}|}\right\\}\mathbf{1}\left\\{{\Phi(q)\in
A}\right\\}\\\
&\qquad+~{}\left({1-\min\left\\{{1,\frac{\pi(\Phi(q))}{\pi(q)}|{\mathrm{det}(\nabla\Phi(q))}|}\right\\}}\right)\mathbf{1}\left\\{{q\in
A}\right\\},\end{split}$ (46)

for $q\in\mathbb{R}^{m}$ and $A\in\mathfrak{B}(\mathbb{R}^{m})$.

It is easily verified that the transition kernel of involutive Monte Carlo
satisfies detailed balance with respect to the density $\pi$. Involutive Monte
Carlo gives rise to two special transition kernels corresponding to RMHMC and
LMC.

######  Definition 21.

Let $\hat{\Phi}_{\epsilon}$ be as in definition 2 and let $H$ be as in
definition 1. The involution of Riemannian manifold Hamiltonian Monte Carlo
with step-size $\epsilon$ and $k$ integration steps is,

$\displaystyle\Phi=\mathbf{F}\circ\underbrace{\hat{\Phi}_{\epsilon}\circ\cdots\circ\hat{\Phi}_{\epsilon}}_{k~{}\mathrm{times}},$
(47)

where $\mathbf{F}(q,p)=(q,-p)$. The target density of Riemannian manifold
Hamiltonian Monte Carlo is $\pi(q,p)\propto\exp(-H(q,p))$. The transition
kernel of involutive Monte Carlo (definition 20) with involution $\Phi$ given
in eq. 47 is called the transition kernel of Riemannian manifold Hamiltonian
Monte Carlo with step-size $\epsilon$ and $k$ integration steps.

######  Remark 8.

It can be shown that if $\Phi$ is an invertible function and if
$\mathbf{F}\circ\Phi$ is an involution that
$\mathbf{F}\circ\Phi\circ\cdots\cdot\Phi$ is also an involution. This explains
why eq. 47 is also an involution. For further details, see proposition 15.

In the special case where $\mathbf{G}$ is a constant function of $q$, the
resulting Markov chain is called Euclidean Hamiltonian Monte Carlo.

######  Definition 22.

Let $\check{\Phi}_{\epsilon}$ be as in definition 5 and let $H$ be as in
definition 4. The involution of Euclidean Hamiltonian Monte Carlo (EHMC) with
step-size $\epsilon$ and $k$ integration steps is,

$\displaystyle\Phi=\mathbf{F}\circ\underbrace{\check{\Phi}_{\epsilon}\circ\cdots\circ\check{\Phi}_{\epsilon}}_{k~{}\mathrm{times}},$
(48)

where $\mathbf{F}(q,p)=(q,-p)$. The target density of Euclidean Hamiltonian
Monte Carlo is $\pi(q,p)\propto\exp(-H(q,p))$. The transition kernel of
involutive Monte Carlo (definition 20) with involution $\Phi$ given in eq. 48
is called the transition kernel of Euclidean Hamiltonian Monte Carlo with
step-size $\epsilon$ and $k$ integration steps.

######  Definition 23.

Let $\tilde{\Phi}_{\epsilon}$ be as in definition 3 and let $H$ be as in
definition 1. The involution of Lagrangian Monte Carlo with step-size
$\epsilon$ and $k$ integration steps is,

$\displaystyle\Phi=\mathbf{F}\circ\underbrace{\tilde{\Phi}_{\epsilon}\circ\cdots\circ\tilde{\Phi}_{\epsilon}}_{k~{}\mathrm{times}},$
(49)

with $\mathbf{F}(q,p)=(q,-p)$. The target density of Lagrangian Hamiltonian
Monte Carlo is $\pi(q,p)\propto\exp(-H(q,p))$. We call this the transition
kernel of Lagrangian Monte Carlo with step-size $\epsilon$ and $k$ integration
steps.

######  Remark 9.

Unlike the transition kernel employed in Riemannian manifold Hamiltonian Monte
Carlo, the Jacobian determinant of the involution appearing in Lagrangian
Monte Carlo must be computed through $k$ applications of eq. 21, where, as
before, $k$ is the number of integration steps.

In the special case where $\mathbf{G}$ is a constant function of $q$,
Riemannian manifold Hamiltonian Monte Carlo and Lagrangian Monte Carlo
produces a transition kernel that is equivalent to Euclidean Hamiltonian Monte
Carlo with constant mass matrix $\mathbf{G}$.

######  Proposition 4.

The transition kernels of RMHMC (or LMC) with step-size $\epsilon$ and $k$
integration steps when applied to the Euclidean Hamiltonian in definition 4
(in the sense that we take $\mathbf{G}(q)\equiv\mathbf{G}$ for every
$q\in\mathbb{R}^{m}$) are both equivalent to the EHMC transition kernel.

Proof: It follows from proposition 1 that the proposals produced by RMHMC and
LMC are identical to EHMC. It remains to be verified that the acceptance
decisions are identical, too. This can be verified by observing that the
Riemannian Hamiltonian in eq. 1 is the Euclidean Hamiltonian in eq. 22 up to
an additive constant (since $\mathbf{G}$ does not depend on position). Hence,
the Euclidean and Riemannian densities are proportional to one another, and
the acceptance probabilities will be identical. We now elaborate on the
connection between EHMC with a single-step and the Metropolis-adjusted
Langevin algorithm: in particular, these two methods can be constructed to be
exactly equivalent in $q$-space, producing equal proposals and identical
acceptance probabilities.

######  Proposition 5 (Neal (2010a)).

Let $\mathbf{A}\in\mathbb{R}^{m\times m}$ be a fixed positive definite matrix.
Consider EHMC with the Hamiltonian
$H(q,p)=-\log\pi(q)+\frac{1}{2}p^{\top}\mathbf{A}^{-1}p$ and MALA with
proposal distribution
$\tilde{q}|q\sim\mathrm{Normal}\left({q+\frac{\epsilon^{2}}{2}\mathbf{A}^{-1}\nabla\log\pi(q),\epsilon^{2}\mathbf{A}^{-1}}\right)$.
The marginal transition kernel of EHMC with step-size $\epsilon$ and a single
integration step can be constructed to be exactly equivalent to the transition
kernel of MALA.

A proof is provided in section A.1. The analysis of mixture transition kernels
will be central to our analysis. In standard implementations of Riemannian
manifold Hamiltonian and Lagrangian Monte Carlo, it is typical to randomize
the number of integration steps. This is done to avert any Markov chain
pathologies (such as irreducibility failures) that may result from using a
fixed number of integration steps. We therefore consider Markov chain
transition kernels which are mixtures.

######  Definition 24.

Let
$K_{\epsilon,k}:(\mathbb{R}^{m}\times\mathbb{R}^{m})\times\mathfrak{B}(\mathbb{R}^{m}\times\mathbb{R}^{m})\to\mathbb{R}_{+}$
be the Markov chain transition kernel of RMHMC (or LMC) with step-size
$\epsilon$ and $k$ integration steps. Let $\tilde{K}_{\epsilon,k}$ be the
marginal transition kernel (defined in lemma 6) of $K_{\epsilon,k}$. Define
the mixture transition kernel of RMHMC (or LMC) to be

$\displaystyle\tilde{K}_{\epsilon}=\sum_{k=1}^{\infty}\alpha_{k}\tilde{K}_{\epsilon,k},$
(50)

where $(\alpha_{1},\alpha_{2},\ldots)$ is a probability vector.

The fact that EHMC, RMHMC, and LMC all satisfy detailed balance in
$(q,p)$-space causes us to examine the progression of $q$-states alone in
between Gibbs resampling steps of the momentum. This leads to marginal Markov
chain transition kernels of the following form:

######  Lemma 6.

Let
$K:(\mathbb{R}^{m}\times\mathbb{R}^{m})\times\mathfrak{B}(\mathbb{R}^{m}\times\mathbb{R}^{m})\to\mathbb{R}$
be a Markov chain transition kernel. Suppose that $K$ satisfies detailed
balance with respect to the density $\pi(q,p)$ with $q$-marginal distribution
$\pi(q)$ and conditional density $\pi(p|q)$; i.e. $\pi(q,p)=\pi(q)\pi(p|q)$.
Consider the marginal chain constructed as follows. Given $q^{n}=q$, sample
$p^{n}\sim\pi(p|q)$, sample $(q^{n+1},p^{n+1})\sim K((q^{n},p^{n}),\cdot)$ and
discard both momenta. The transition kernel of the marginal chain satisfies

$\displaystyle\tilde{K}(q,Q)=\int_{\mathbb{R}^{m}}K((q,p),(Q,\mathbb{R}^{m}))~{}\pi(p|q)~{}\mathrm{d}p$
(51)

where $Q\in\mathfrak{B}(\mathbb{R}^{m})$.

A proof is given in section A.2

######  Proposition 7.

Let $\tilde{K}$ be the marginal transition kernel described in lemma 6. The
marginal chain satisfies detailed balance with respect to the distribution
whose density is $\pi(q)$.

A proof is provided in section A.3.

####  2.2.3 Ergodicity Theorems

The geometric ergodicity of MALA, MMALA, and SMALA have been examined in the
Markov chain literature. In section 3 we will see how the geometric ergodicity
of a single component of a mixture Markov chain transition kernel implies the
geometric ergodicity of the mixture itself. With this destination in mind, we
now recall two key results in this direction.

######  Theorem 8 (Roberts and Tweedie (1996)).

Consider the transition kernel $K$ of the Metropolis-adjusted Langevin
algorithm with $\mathbf{A}=\mathrm{Id}$ and
$c_{\epsilon}(q)=x+\frac{\epsilon^{2}}{2}\nabla\log\pi(q)$. Then $K$ is
geometrically ergodic under the following conditions:

  1. 1.

We have

$\displaystyle\liminf_{\|q\|\to\infty}\left({\|q\|-\|q+\frac{\epsilon^{2}}{2}\nabla\log\pi(q)\|}\right)>0.$
(52)

  2. 2.

We have

$\displaystyle\lim_{\|q\|\to\infty}\int_{(A(q)\cup I(q))\cap(A(q)\cap
I(q))}\tilde{\pi}(\tilde{q}|q)~{}\mathrm{d}\tilde{q}=0,$ (53)

where

$\displaystyle A(q)$
$\displaystyle=\left\\{{\tilde{q}\in\mathcal{X}:\pi(q)\tilde{\pi}(\tilde{q}|q)\leq\pi(y)\tilde{\pi}(q|\tilde{q})}\right\\}$
(54) $\displaystyle I(q)$
$\displaystyle=\left\\{{\tilde{q}\in\mathcal{X}:\|\tilde{q}\|\leq\|q\|}\right\\}.$
(55)

######  Theorem 9 (Roy and Zhang (2023)).

Consider the transition kernel $K$ of the generalized Metropolis-adjusted
Langevin algorithm. Then $K$ is geometrically ergodic under the following
conditions:

  1. 1.

There exist matrices $\mathbf{A}_{l}\in\mathrm{PD}(m)$ and
$\mathbf{A}_{u}\in\mathrm{PD}(m)$ such that
$\mathbf{A}_{l}\leq\mathbf{A}(q)\leq\mathbf{A}_{u}$.

  2. 2.

When $A\subset\mathbb{R}^{m}$ is bounded, the function
$c_{\epsilon}:\mathbb{R}^{m}\to\mathbb{R}^{m}$ is bounded on $A$.

  3. 3.

There is a quantity

$\displaystyle
C=\limsup_{\|q\|\to\infty}\int_{R(q)}\left({1-\tilde{\pi}(\tilde{q}|q)\min\left\\{{1,\frac{\pi(\tilde{q})\tilde{\pi}(q|\tilde{q})}{\pi(q)\tilde{\pi}(\tilde{q}|q)}}\right\\}}\right)~{}\mathrm{d}y$
(56)

which is strictly less than one, where

$\displaystyle
R(q)=\left\\{{\tilde{q}\in\mathbb{R}^{m}:\pi(q)\tilde{\pi}(\tilde{q}|x)>\pi(\tilde{q})\tilde{\pi}(q|\tilde{q})}\right\\}.$
(57)

  4. 4.

There exists $s>0$ such that,

$\displaystyle\liminf_{\|q\|\to\infty}\left({\|\mathbf{A}^{-1/2}_{u}(q)x\|-\|\mathbf{A}^{-1/2}_{u}(q)c_{\epsilon}(q)\|}\right)>\frac{\log(D(s))-\log(1-C)}{s},$
(58)

where,

$\displaystyle
D(s)=\epsilon^{-m/2}\left({\frac{\pi}{2}}\right)^{m/2-1}\left({\frac{\mathrm{det}(\mathbf{A}_{u})}{\mathrm{det}(\mathbf{A}_{l})}}\right)^{1/2}\exp(\epsilon
s^{2}/2)\int_{0}^{\infty}\exp\left({-(r-\epsilon
s)^{2}/2\epsilon}\right)r^{m-1}~{}\mathrm{d}r$ (59)

In section 3 we will examine mixture transition kernels and their geometric
ergodicity. In the case that one transition kernel is geometrically ergodic
and all transition kernels are reversible, the mixture transition kernel is
geometrically ergodic, too, from the following result.

######  Theorem 10 (Lee and Łatuszyński (2014)).

Let $(\alpha_{1},\alpha_{2},\ldots)$ be a sequence satisfying $\alpha_{i}\geq
0$ and $\sum_{i=1}^{\infty}\alpha_{i}=1$. Let
$\tilde{K}=\sum_{i=1}^{\infty}\alpha_{i}K_{i}$ be a mixture of reversible
transition kernels with invariant distribution $\Pi$. Suppose that $K_{1}$
satisfies the properties: (i) $\Pi$ is the unique invariant distribution, and
(ii) $K_{1}$ is geometrically ergodic. Then $\tilde{K}$ is geometrically
ergodic.

We term the derived geometric ergodicity of $\tilde{K}$ from $K_{1}$ as
“inherited” geometric ergodicity.

##  3 Inherited Geometric Ergodicity

In this section we discuss the notion of “inherited” geometric ergodicity,
wherein we modify RMHMC and LMC to be geometrically ergodic when MMALA is. In
the Euclidean case, by proposition 5, in cases where the Metropolis-adjusted
Langevin algorithm is geometrically ergodic, so is HMC by invoking theorem 10.
The situation in the geometric setting is more complicated, since there is no
analogue of proposition 5 to show that the marginal transition kernel of a
single-step of RMHMC or LMC is exactly equivalent to the Markov chain
transition kernel of the generalized Metropolis-adjusted Langevin algorithm
for a suitable choice of mean function. This leads us to propose the following
Markov chain transition kernel.

######  Definition 25.

Let
$K_{\epsilon,k}:(\mathbb{R}^{m}\times\mathbb{R}^{m})\times\mathfrak{B}(\mathbb{R}^{m}\times\mathbb{R}^{m})\to\mathbb{R}_{+}$
be the Markov chain transition kernel of RMHMC (or LMC) with step-size
$\epsilon$ and $k$ integration steps. Let $\tilde{K}_{\epsilon,k}$ be the
marginal transition kernel (defined in lemma 6). Let
$(\alpha_{1},\alpha_{2},\ldots)$ be a sequence satisfying $\alpha_{i}\geq 0$
and $\sum_{i=1}^{\infty}\alpha_{i}=1$. Let $J_{\epsilon}$ be the Markov chain
transition kernel of the generalized Metropolis-adjusted Langevin algorithm
(defined in definition 17). The Langevin mixture transition kernel of RMHMC
(or LMC), which we abbreviate by LMRMHMC (or LMLMC), is defined by

$\displaystyle\tilde{K}_{\epsilon}=\alpha_{1}J_{\epsilon}+\sum_{k=2}^{\infty}\alpha_{k}\tilde{K}_{\epsilon,k}.$
(60)

We call $\alpha_{1}$ the MMALA mixture weight.

The modified transition kernel simply replaces a single-step of RMHMC (or LMC)
by the transition kernel of the generalized Metropolis-adjusted Langevin
algorithm. In order to apply theorem 10, it is necessary to verify that the
marginal transition kernels of RMHMC and LMC are reversible. Fortunately, this
is readily shown as follows.

######  Corollary 11.

For fixed $\epsilon\in\mathbb{R}$ and $k\in\mathbb{N}$, the RMHMC and LMC
transition kernels with step-size $\epsilon$ and $k$ integration steps satisfy
detailed balance in $(q,p)$-space, it follows that their marginal chains are
reversible with respect to the distribution $\Pi$ by proposition 7.

######  Proposition 12.

Suppose $\alpha_{1}>0$. Under the conditions of theorem 9, the Markov chain
transition kernel of LMRMHMC (or LMLMC) in eq. 60 is geometrically ergodic.

Proof: This follows as an immediate consequence of theorem 10 using the fact
that the marginal transition kernels or RMHMC (or LMC) are reversible by
corollary 11.

######  Remark 10.

As a practical matter, we choose the mixture probabilities
$(\alpha_{1},\alpha_{2},\ldots)$ in the following way. For a particular target
distribution, we will choose a maximal number of integration steps
$k_{\mathrm{max}}$ for which $\alpha_{k}=0$ for $k>k_{\mathrm{max}}$. Then,
given a particular selection of $\alpha_{1}\in(0,1]$, we split the remaining
probability mass equally for each
$k\in\left\\{{2,\ldots,k_{\mathrm{max}}}\right\\}$; that is, the fraction can
be expressed as $\alpha_{k}=\frac{1-\alpha_{1}}{k_{\mathrm{max}}-1}$ for
$k=2,\ldots,k_{\mathrm{max}}$. In our experiments, we will consider variable
choices for the mixing parameter $\alpha_{1}$. The special case
$\alpha_{1}=0$, will correspond by convention to the unmodified RMHMC and LMC
transition kernels with a single-step computed using the prescribed
involution, as described in definition 24.

Mixing with the MMALA transition kernel also immediately establishes
irreducibility, aperiodicity, and the smallness of all compact sets, as the
following result reveals.

######  Lemma 13.

Let $\pi(q)\propto\exp(\mathcal{L}(q))$ be continuous and bounded on compact
sets and denote by $\Pi$ the probability measure with density $\pi$. Suppose
$\alpha_{1}>0$. The marginal transition kernel of LMRMHMC (or LMLMC) is
$\Pi$-irreducible, aperiodic, and all non-negligible compact sets are small.

Proof: With probability $\alpha_{1}$, the Markov chain transitions according
to a Metropolis-Hastings accept-reject decision with a normal proposal
distribution. Hence, for any set $A\subset\mathfrak{B}(\mathbb{R}^{m})$ for
which $\int_{A}\pi(q)~{}\mathrm{d}q>0$, we have,

$\displaystyle\tilde{K}_{\epsilon}(q,A)$
$\displaystyle\geq\alpha_{1}J_{\epsilon}(q,A)$ (61) $\displaystyle>0$ (62)

since a normal proposal distribution is non-vanishing everywhere on
$\mathbb{R}^{m}$. The fact that $\tilde{K}_{\epsilon}$ is aperiodic and that
all non-negligible compact sets are small follows as an immediate consequence
of Lemma 1.2 from Mengersen and Tweedie (1996).

######  Corollary 14.

From lemma 13, we have that the modified Markov chain of RMHMC (or LMC) is
$\Pi$-irreducible, aperiodic, and from proposition 7 $\Pi$ is the stationary
distribution. Hence, it follows from theorem 2 that the modified Markov chain
of RMHMC (or LMC) produces an ergodic Markov chain.

##  4 Experimentation

We turn now to the investigation of the proposed modified variations of RMHMC
and LMC. These examples are chosen to represent a wide class of posterior
distributions. In computing the convergence of the Markov chain under the
maximum mean discrepancy metric, we measure similarity between 10,000 i.i.d.
samples of the target distribution, and an equal number of independent Markov
chains; in this case, we compute the unbiased estimator of the maximum mean
discrepancy over the course of the first one-hundred sampling steps. For the
expected squared jump distance (ESJD) and the effective sample size (ESS)
metrics, we consider a ten replicates of a long Markov chain consisting of
1,000,000 samples, except in the case of the Fitzhugh-Nagumo posterior where
we sample only 100,000 times. In each experiment, we consider mixing LMC and
RMHMC with the transition kernel described in definition 18, except in the
case of Neal’s funnel distribution, wherein we consider a mixture with
definition 19. We begin in section 4.1 to describe measures and metrics by
which we may assess the performance of the Markov chains corresponding to
LMRHMC and LMLMC.

Code for reproducing these experimental results may be found at
https://tinyurl.com/29kz7krz.

###  4.1 Measures and Metrics

We now give details of the evaluation metrics by which we compare Markov
chains. We first recall the method of maximum mean discrepancy due to Gretton
et al. (2012).

######  Definition 26.

Let $k:\mathbb{R}^{m}\times\mathbb{R}^{m}\to\mathbb{R}$ be a positive definite
function that is symmetric in its arguments. Let $\Pi$ and $\Pi^{\prime}$ be
two probability distributions on $\mathbb{R}^{m}$. The squared maximum mean
discrepancy between $\Pi$ and $\Pi^{\prime}$ is defined by,

$\displaystyle\mathrm{MMD}^{2}(k,\Pi,\Pi^{\prime})=\underset{q,q^{\prime}\sim\Pi}{\mathbb{E}}k(q,q^{\prime})+\underset{q,q^{\prime}\sim\Pi^{\prime}}{\mathbb{E}}k(q,q^{\prime})-2\underset{q\sim\Pi,q^{\prime}\sim\Pi^{\prime}}{\mathbb{E}}k(q,q^{\prime}).$
(63)

Let $(q^{1}_{*},\ldots,q^{r}_{*})\sim\Pi$ and
$(q^{1},\ldots,q^{s})\sim\Pi^{\prime}$. An unbiased estimator of the squared
maximum mean discrepancy is,

$\displaystyle\begin{split}\mathrm{MMD}^{2}_{\mathrm{u}}(k,\left\\{{q^{i}_{*}}\right\\}_{i=1}^{r},\left\\{{q^{i}}\right\\}_{i=1}^{s})&=\frac{1}{r(r-1)}\sum_{i=1}^{r}\sum_{j\neq
i}^{r}k(q^{i}_{*},q^{j}_{*})+\frac{1}{s(s-1)}\sum_{i=1}^{s}\sum_{j\neq
i}^{s}k(q^{i},q^{j})\\\
&\qquad-~{}\frac{2}{rs}\sum_{i=1}^{r}\sum_{j=1}^{s}k(q^{i}_{*},q^{j}).\end{split}$
(64)

In our evaluations we adopt a squared exponential positive definite kernel
$k(q,q^{\prime})=\exp(-\|q-q^{\prime}\|^{2}/2h)$, where $h\in\mathbb{R}_{+}$
is a parameter called the kernel bandwidth. In our experiments we set $h$ to
be the median distance between i.i.d. samples from the target distribution.
Recall the definition of geometric ergodicity given in definition 14. The
measure of convergence is the total variation norm, which we can discuss
theoretically but cannot evaluate in a computational setting. Instead, we can
measure convergence to the target distribution $\Pi$ as a function of $n$, the
number of steps, by means of eq. 63, for which we can obtain an unbiased
estimate via eq. 64 if we have samples from $K^{n}(q^{0},\cdot)$ and samples
from $\Pi$. Indeed, as in definition 26, let
$(q^{1}_{*},\ldots,q^{r}_{*})\sim\Pi$ and, for a fixed initial position
$q^{0}$ and number of steps $n$ let $(q^{1},\ldots,q^{s})\sim
K^{n}(q^{0},\cdot)$ and we can compute an unbiased estimate of the squared
maximum mean discrepancy. In the latter case, $(q^{1},\ldots,q^{s})$ can be
obtained by running $s$ independent Markov chains for $n$ steps; independent
samples from the target distribution $\Pi$ may be available in certain
benchmark cases. In our experiments, i.i.d. samples from the target
distribution may be generated from the banana-shaped posterior, Neal’s funnel
distribution, the Fitzhugh-Nagumo posterior, and the multi-scale Student-$t$
distribution.

######  Remark 11.

Recall the definition of geometric ergodicity given in definition 14. Taking
logarithms reveals

$\displaystyle\log\|K^{n}(q,\cdot)-\Pi(\cdot)\|_{\mathrm{TV}}\leq
n\log\rho+\log V(q).$ (65)

Therefore, one may claim to see evidence of geometric ergodicity if, as a
function of $n$, there is a linear decrease in the total variation distance on
a logarithmic scale. Of course, we cannot directly compute the total variation
distance, but we may look for a similar negative linear trend when
$|{\mathrm{MMD}_{u}^{2}}|$ is plotted on a logarithmic scale.

As an additional measure of ergodicity, we consider comparing Markov chain
samples against i.i.d. samples via random projection onto one-dimensional sub-
spaces. Let $(q^{1}_{*},\ldots,q^{r}_{*})\sim\Pi$ and let $q^{n}\sim
K^{n}(q^{0},\cdot)$ for $n=1,\ldots,s$. Let $u$ be a random unit vector. We
compute the Kolmogorov-Smirnov (KS) statistic for the projections
$(u^{\top}q^{1}_{*},\ldots,u^{\top}q^{r}_{*})$ and
$(u^{\top}q^{1},\ldots,u^{\top}q^{s})$. Repeating this process for one-hundred
randomly generated unit vectors yields a distribution over Kolmogorov-Smirnov
statistics. The more tightly concentrated this distribution is near zero, the
closer the distribution of Markov chain iterates $(q^{1},\ldots,q^{s})$ is to
the collection of i.i.d. samples $(q^{1}_{*},\ldots,q^{r}_{*})$ from the
target distribution.

We also consider the expected squared jump distance (ESJD) (Gelman and
Pasarica, 2007), which measures the dissimilarity between subsequent states of
the Markov chain. Intuitively, a Markov chain that moves more efficiently
through the sample space (higher ESJD) will exhibit smaller sample auto-
correlation. Formally, the ESJD is
$\underset{q\sim\Pi}{\mathbb{E}}\|q^{\prime}-q\|^{2}$ where $q^{\prime}|q\sim
K(q,\cdot)$. Note that the choice of norm $\|\cdot\|$ is left to the
practitioner and may be selected to capture geometric properties of the target
distribution. This expectation is typically approximated as follows: let
$q^{k}$ be the state of the Markov chain at step $k$ and at step $k+1$ a
candidate state is generated, denoted $q^{k+1}_{\mathrm{prop}}$; the proposal
state is accepted with probability $\alpha_{k}$; the following empirical mean
is then taken as our approximation to the ESJD

$\displaystyle\mathrm{ESJD}(\left\\{{(q^{k},q^{k+1}_{\mathrm{prop}},\alpha_{k})}\right\\}_{k=1}^{n})=\frac{1}{n}\sum_{k=1}^{n}\alpha_{k}\|q^{k+1}_{\mathrm{prop}}-q^{k}\|^{2}.$
(66)

In our experiments in section 4.8, we observe that the ESJD can be misleading
as a measure. Therefore, we also introduce the median squared jump distance
(MSJD) which we define as

$\displaystyle\mathrm{MSJD}(\left\\{{(q^{k},q^{k+1}_{\mathrm{prop}},\alpha_{k})}\right\\}_{k=1}^{n})=\mathrm{Median}\left({\left\\{{\alpha_{k}\|q^{k+1}_{\mathrm{prop}}-q^{k}\|^{2}}\right\\}_{k=1}^{n}}\right).$
(67)

As the median, we expect the MSJD to exhibit less sensitivity to outliers than
the ESJD.

We additionally consider the effective sample size (ESS) as a metric for our
Markov chain procedures. We compute ESS using the technique of Kumar et al.
(2019), who describe the method succinctly as follows. For $i=1,\ldots,p$, let
$(q^{i,1},q^{i,2},\ldots,q^{i,n})$ be a sequence of $\mathbb{R}^{m}$-valued
parameters. The integer $p$ is called the number of chains. The effective
sample size of the $j$-th parameter is computed according to,

$\displaystyle\mathrm{ESS}\left({\left\\{{q^{i,1}_{j},\ldots,q^{i,n}_{j}}\right\\}_{i=1}^{p}}\right)$
$\displaystyle=\frac{pn}{\hat{\tau}}$ (68) $\displaystyle\hat{\tau}$
$\displaystyle=\left({2\sum_{k=1}^{r}\hat{\rho}_{2k}+\hat{\rho}_{2k+1}}\right)-1,$
(69)

where $\hat{\rho}_{k}$ is an estimate, based on all of the $p$ sequences, of
the autocorrelation of the $j$-th parameter with a $k$-step lag (for details
see Vehtari et al. (2021)) and $r$ is the smallest integer for which
$\hat{\rho}_{2(r+1)}+\hat{\rho}_{2(r+1)+1}\leq 0$. In our experiments, we set
$p=2$ by taking a single long Markov chain and splitting it in half at the
middle. As a practical matter, we will also report the minimum, over all
parameters of the posterior, ESS per second in order to represent the
computational efficiency of the method.

###  4.2 Banana-Shaped Distribution

(a) EHMC

(b) LMRMHMC

(c) LMLMC

Figure 1: We examine the ergodicity of the Markov chains as a function of the
sampling step in the banana-shaped posterior distribution. We observe that the
use of a smaller step-size degrades the ergodicity of the RMHMC, particularly
when compared against the LMC algorithm. In both RMHMC and LMC, we observe
that aggressively mixing with MMALA causes the Markov chain to mix more
slowly; on the other hand, a relatively small mixing probability leads to
ergodicity that is nearly indistinguishable from the unmodified
implementations of RMHMC and LMC.

(a) ESJD

(b) Min. ESS

(c) Min. ESS / Sec.

Figure 2: We show the ESJD, the minimum ESS, and the time-normalized minimum
ESS for inference in the banana-shaped posterior distribution. We observe that
for both RMHMC and LMC, aggressively mixing with MMALA causes the ESJD
distance to decrease; this correspondingly produces a decrease in the minimum
ESS. For LMC, the time-normalized minimum ESS is a decreasing function of the
mixing probability.

Our first example considers the following generative model, which represents
an example of non-identifiable parameters.

$\displaystyle\theta_{1},\theta_{2}$
$\displaystyle\overset{\mathrm{i.i.d.}}{\sim}\mathrm{Normal}(0,\sigma_{\theta}^{2})$
(70) $\displaystyle y_{i}|\theta_{1},\theta_{2}$
$\displaystyle\overset{\mathrm{i.i.d.}}{\sim}\mathrm{Normal}(\theta_{1}+\theta_{2}^{2},\sigma_{y}^{2})~{}~{}\mathrm{for}~{}i=1,\ldots,N.$
(71)

Given observations $\left\\{{y_{i}}\right\\}_{i=1}^{N}$, we wish to sample the
posterior distribution of $(\theta_{1},\theta_{2})$, a two-dimensional
posterior. In our experiments, we consider $n=100$. We employ Euclidean HMC
with a step-size of $\epsilon=0.1$; RMHMC is implemented with a step-size of
0.04; in LMC we use an integration step-size of $\epsilon=0.1$. In each case,
we set $k_{\mathrm{max}}=10$. For the geometric methods, the sum of the Fisher
information and the negative Hessian of the log-prior is used as a Riemannian
metric. The reason that a smaller integration step-size is employed in RMHMC
is that eq. 5 in the defining involution will fail to have a solution for
$\breve{p}$ for large step-sizes; this necessitates the use of a smaller
integration step.

In fig. 1, we visualize the statistic $|{\mathrm{MMD}_{u}^{2}}|$ over one-
hundred steps of the Markov chain for euclidean HMC (EHMC), RMHMC and LMC; we
color the RMHMC and LMC Markov chains according to how aggressively they mix
with MMALA. Neither EHMC nor RMHMC exhibit clear evidence of geometric
ergodicity in this target distribution; on the other hand, there exist a broad
range of mixing probabilities for which LMC exhibits a linear decrease, as a
function of $n$, in $|{\mathrm{MMD}_{u}^{2}}|$ on a logarithmic scale. In fig.
2 we visualize the ESJD, the minimum ESS, and the minimum ESS per second for
the three MCMC algorithms. We observe that LMC exhibits by far the strongest
performance on these metrics, while RMHMC languishes due to its small step-
size and the pathologies of the generalized leapfrog integrator applied to
this posterior.

###  4.3 Hierarchical Bayesian Logistic Regression

(a) ESJD

(b) Min. ESS

(c) Min. ESS / Sec.

Figure 3: We show the ESJD, the minimum ESS, and the time-normalized minimum
ESS for inference in the hierarchical Bayesian logistic regression posterior
distribution. We observe that both of the geometric methods do not perform as
well as standard HMC in this inference task. For each mixing weight with
MMALA, we observe that RMHMC enjoys a greater ESJD and minimum ESS when
compared against LMC; however, as LMC is a more computationally expedient
procedure, LMC has a greater time-normalized ESS. In all of these metrics,
however, EHMC dominates the geometric inference algorithms.

We consider sampling from the hierarchical Bayesian logistic regression model

$\displaystyle\alpha$ $\displaystyle\sim\mathrm{Gamma}(\omega,\theta)$ (72)
$\displaystyle\beta_{i}|\alpha$
$\displaystyle\overset{\mathrm{i.i.d.}}{\sim}\mathrm{Normal}(0,\alpha^{-1})~{}~{}\mathrm{for}~{}i=1,\ldots,m$
(73) $\displaystyle y_{i}|x_{i},\beta$
$\displaystyle\overset{\mathrm{indep.}}{\sim}\mathrm{Bernoulli}\left({\frac{1}{1+\exp(-x_{i}^{\top}\beta)}}\right)~{}~{}\mathrm{for}~{}i=1,\ldots,N.$
(74)

We consider a logistic regression dataset with $N=270$ observations and $m=14$
covariates. We set $\omega=10$ and $\theta=2$ in our experiments. We employ a
Metropolis-within-Gibbs sampling procedure wherein we alternate between
sampling the posterior distributions $\alpha|\beta,k,\theta$ and
$\beta|\alpha,\left\\{{(x_{i},y_{i})}\right\\}_{i=1}^{n}$; sampling the former
can be performed analytically, whereas we employ EHMC, RMHMC, and LMC to
sample the latter. In Euclidean HMC, we set $\epsilon=0.1$; in RMHMC and LMC
we set $\epsilon=0.8$. We set $k_{\mathrm{max}}=10$ as the upper bound on the
number of integration steps in each case. For implementing both LMRMHMC and
LMLMC, we use the sum of the Fisher information and the negative Hessian of
the log-prior as a metric.

In fig. 3 we show the ESJD, the minimum ESS, and the minimum ESS per second
for the hierarchical Bayesian logistic regression posterior. Neither of the
geometric methods perform well in this posterior, consistently under-
performing EHMC. A criticism of LMC is that its performance degrades
significantly in higher dimensions (Betancourt et al., 2014); one sees
evidence of this phenomenon in the smaller ESS generated by LMC; however, the
computational savings due to eliminating the fixed point iterations still
allow LMC to edge out a stronger time-normalized ESS compared to RMHMC.

###  4.4 Neal’s Funnel Distribution

(a) EHMC

(b) LMRMHMC

Figure 4: We examine the ergodicity of the Markov chains as a function of the
sampling step in Neal’s funnel distribution. We observe that EHMC struggles in
this distribution due to the multiple spatial scales inherent in the funnel-
shaped distribution. Moreover, we observe that aggressively mixing RMHMC with
SMALA causes convergence in MMD to become increasingly slow, with convergence
of the pure SMALA algorithm being slower than EHMC. On the other hand, a
modest amount of mixing with SMALA produces ergodicity that is nearly
indistinguishable from the unmodified implementation of RMHMC.

(a) ESJD

(b) KS

(c) Min. ESS / Sec.

Figure 5: We show the ESJD, the distribution of Kolmogorov-Smirnov statistics,
and the time-normalized minimum ESS for inference in Neal’s funnel
distribution. Although RMHMC and EHMC have comparable time-normalized ESS, it
is clear from the distribution of KS statistics that RMHMC produces samples
that are closer to the target distribution in EHMC.

Neal’s funnel distribution (Neal, 2003) is a density defined in the following
hierarchical manner.

$\displaystyle v$ $\displaystyle\sim\mathrm{Normal}(0,9)$ (75) $\displaystyle
x_{i}|v$
$\displaystyle\overset{\mathrm{i.i.d.}}{\sim}\mathrm{Normal}(0,\exp(-v))~{}~{}\mathrm{for}~{}i=1,\ldots,N.$
(76)

This distribution is shaped like a funnel, in which the thickness of the
“neck” is being controlled by the random variable $v$. This model is
reflective of posteriors encountered in hierarchical models with sparse data.
The objective in this task is to jointly sample $(v,x_{1},\ldots,x_{N})$,
producing a $(N+1)$-dimensional target distribution. In Euclidean HMC we
employ an integration step-size of $\epsilon=0.1$ and $k_{\mathrm{max}}=10$.
Our implementation of RMHMC uses $k_{\mathrm{max}}=20$ integration steps with
a step-size of $\epsilon=0.1$ with the SoftAbs metric. We do not consider LMC
in this task since we found it non-obvious how the SoftAbs structure could be
extended into the LMC framework while preserving the cubic computational cost
at each step. In this experiment, we use the SoftAbs metric.

In fig. 4, we visualize $|{\mathrm{MMD}_{u}^{2}}|$ for both EHMC and RMHMC. As
expected, EHMC struggles to sample from Neal’s funnel distribution due to the
multiscale phenomena. On the other hand, RMHMC exhibits much stronger
convergence properties, having a linear decrease over several possible mixing
probabilities with SMALA. We find that aggressively mixing with SMALA can be
counter-productive, however, due to the less efficient traversal of the target
distribution by single-step methods. In fig. 5 we show the ESJD, the
Kolmogorov-Smirnov statistics, and the minimum ESS per second. On all of these
metrics, RMHMC clearly outperforms EHMC.

###  4.5 Stochastic Volatility Model

(a) ESJD

(b) Min. ESS

(c) Min. ESS / Sec.

Figure 6: We show the ESJD, the minimum ESS, and the time-normalized minimum
ESS for inference in the posterior distribution of the stochastic volatility
model. We observe that both of the geometric methods with a modest mixing
probability with MMALA produce time-normalized ESS that are competitive with
EHMC. For each mixing weight with MMALA, we observe that LMC enjoys a greater
ESJD when compared against RMHMC.

We consider a stochastic volatility model with the following generative model.

$\displaystyle\frac{\phi+1}{2}$ $\displaystyle\sim\mathrm{Beta}(20,3/2)$ (77)
$\displaystyle\sigma^{2}$
$\displaystyle\sim\mathrm{InverseChiSquared}(10,1/20)$ (78) $\displaystyle
x_{1}|\phi,\sigma^{2}$
$\displaystyle\sim\mathrm{Normal}(0,\sigma^{2}/(1-\phi^{2}))$ (79)
$\displaystyle x_{t+1}|x_{t},\phi,\sigma^{2}$
$\displaystyle\sim\mathrm{Normal}(\phi
x_{t},\sigma^{2})~{}~{}\mathrm{for}~{}t=2,\ldots,T-1$ (80) $\displaystyle
y_{t}|\beta,x_{t}$
$\displaystyle\sim\mathrm{Normal}(0,\beta^{2}\exp(x_{t}))~{}~{}\mathrm{for}~{}t=1,\ldots,T.$
(81)

Additionally, the parameter $\beta$ is equipped with an improper prior
proportional to $\beta^{-1}$. In this example, we seek to generate samples
from the joint distribution $x_{1},\ldots,x_{T},\beta,\phi,\sigma^{2}$ given
observations $y_{1},\ldots,y_{T}$. We employ Metropolis-within-Gibbs-like
strategy wherein we alternate between sampling
$(x_{1},\ldots,x_{T})|(y_{1},\ldots,y_{T})\phi,\sigma^{2},\beta$ and
$(\phi,\sigma^{2},\beta)|\left\\{{(x_{i},y_{i})}\right\\}_{i=1}^{T}$; in the
former case we employ Euclidean HMC whereas in the latter case we compare
Euclidean HMC, RMHMC, and LMC. In sampling either distribution, the Riemannian
metric is chosen as the sum of the Fisher information and the negative Hessian
of the log-prior. In our experiments we set $T=1,000$. When using Euclidean
HMC to sample $(\phi,\sigma^{2},\beta)$ we use $k_{\mathrm{max}}=50$
integration steps and a step-size of $\epsilon=0.01$; in the case of RMHMC and
LMC we use $k_{\mathrm{max}}=6$ integration steps and a step-size of
$\epsilon=0.5$.

In fig. 6 we show the ESJD, the minimum ESS, and the minimum ESS per second.
We observe that employing a modest mixture probability with MMALA produces a
Markov chain that is marginally better than EHMC.

###  4.6 Log-Gaussian Cox-Poisson Process

(a) ESJD

(b) Min. ESS

(c) Min. ESS / Sec.

Figure 7: We show the ESJD, the minimum ESS, and the time-normalized minimum
ESS for inference in the posterior distribution of the log-Gaussian Cox-
Poisson model. For any mixing probability, the geometric methods enjoy larger
time-normalized effective sample sizes than EHMC, with LMC outperforming RMHMC
on these metrics. For nearly every mixing probability, we additionally observe
that the ESJD is greater for LMC and RMHMC than for EHMC.

We consider inference in a log-Gaussian Cox-Poisson model with the following
generative model:

$\displaystyle\Sigma_{(i,j),(i^{\prime},j^{\prime})}|\sigma^{2},\beta$
$\displaystyle=\sigma^{2}\exp\left({-\sqrt{(i-i^{\prime})^{2}+(j-j^{\prime})^{2}}/(N\beta)}\right)$
(82) $\displaystyle\mathrm{vec}(\mathbf{x})|\Sigma$
$\displaystyle\sim\mathrm{MultivariateNormal}(\mu\mathbf{1},\Sigma)$ (83)
$\displaystyle y_{ij}|x_{ij}$
$\displaystyle\sim\mathrm{Poisson}(\exp(x_{ij})/N^{2}),$ (84)

with priors $\beta\sim\mathrm{Gamma}(2,1/2)$ and
$\sigma^{2}\sim\mathrm{Gamma}(2,1/2)$. In this example, the objective is to
sample the joint distribution
$(\beta,\sigma^{2},\left\\{{x_{ij}}\right\\}_{i,j=1}^{N})$ given observations
$\left\\{{y_{ij}}\right\\}_{i,j=1}^{N}$. As in the case of stochastic
volatility model, we alternatively sample between
$\left\\{{x_{ij}}\right\\}_{i,j=1}^{N}$ given
$\left\\{{y_{ij}}\right\\}_{i,j=1}^{N}$, $\sigma^{2}$, and $\beta$, and
$(\beta,\sigma^{2})$ given $\left\\{{x_{ij}}\right\\}_{i,j=1}^{N}$. In each
case, the metric is given by the sum of the Fisher information and the
negative Hessian of the log-prior. In our experiments we set $N=16$. In the
former case we employ Euclidean HMC, whereas in the latter case we compare
Euclidean HMC, RMHMC and LMC. When implementing Euclidean HMC we employ
$k_{\mathrm{max}}=50$ integration steps and a step-size of $\epsilon=0.01$; in
RMHMC and LMC we use $k_{\mathrm{max}}=6$ integration steps and a step-size of
$\epsilon=0.5$.

In fig. 7 we show the ESJD, the minimum ESS, and the minimum ESS per second.
We observe that there is a range of mixture probabilities for which mixing
with the MMALA in LMRMHMC and LMLMC exhibit superior performance compared to
EHMC. Both RMHMC and LMC exhibit similar minimum ESS metrics, but due to its
computational advantage, LMC produces a larger minimum ESS per second.

###  4.7 Fitzhugh-Nagumo Model

(a) EHMC

(b) LMRMHMC

(c) LMLMC

Figure 8: We examine the ergodicity of the Markov chains as a function of the
sampling step in the Fitzhugh-Nagumo posterior distribution. We observe that
when a modest mixing probability is employed, the modified RMHMC and LMC
transition kernels mix more efficiently than the EHMC transition kernel, being
nearly indistinguishable from the unmodified RMHMC and LMC transition kernels.
When one uses a large mixing probability, we see that the mixing rate is
decreased.

(a) ESJD

(b) Min. ESS

(c) Min. ESS / Sec.

Figure 9: We show the ESJD, the minimum ESS, and the time-normalized minimum
ESS for inference in the Fitzhugh-Nagumo posterior distribution. LMC dominates
both RMHMC and EHMC in this inference task under the time-normalized ESS
metrics with its greater computational expediency. When timing is not
accounted for, we observe that RMHMC and LMC enjoy similar ESJD and minimum
ESS metrics. We observe that for both RMHMC and LMC, aggressively mixing with
MMALA causes the ESJD distance to decrease; this correspondingly produces a
decrease in the minimum ESS.

Given $\mathbb{R}$-valued parameters $a$, $b$, and $c$, the Fitzhugh-Nagumo
differential equations are defined by,

$\displaystyle\dot{v}_{t}$
$\displaystyle=c\left({v_{t}-\frac{v_{t}^{3}}{3}+r_{t}}\right)$ (85)
$\displaystyle\dot{r}_{t}$
$\displaystyle=-\left({\frac{v_{t}-a+br_{t}}{c}}\right).$ (86)

Given initial conditions $v_{0}$ and $r_{0}$, we consider the following
generative model:

$\displaystyle(a,b,c)$
$\displaystyle\overset{\mathrm{i.i.d.}}{\sim}\mathrm{Normal}(0,1)$ (87)
$\displaystyle\hat{r}_{t_{k}}|a,b,c,r_{t_{k}},\sigma^{2}$
$\displaystyle\overset{\mathrm{indep.}}{\sim}\mathrm{Normal}(r_{t_{k}},\sigma^{2})~{}~{}\mathrm{for}~{}k=1,\ldots,n$
(88) $\displaystyle\hat{v}_{t_{k}}|a,b,c,v_{t_{k}},\sigma^{2}$
$\displaystyle\overset{\mathrm{indep.}}{\sim}\mathrm{Normal}(v_{t_{k}},\sigma^{2})~{}~{}\mathrm{for}~{}k=1,\ldots,n,$
(89)

where $t_{1},\ldots,t_{n}$ are equally spaced points between $[0,T]$. In our
experiments, we set $v_{0}=1$, $r_{0}=-1$, $\sigma^{2}=1/4$, $T=10$, and
$n=200$. Our metric is given by the sum of the Fisher information and negative
Hessian of the log-prior. In our implementations, we use Euclidean HMC with a
step-size of $\epsilon=0.01$ and $k_{\mathrm{max}}=10$ integration steps. In
RMHMC and LMC, we employ an integration step-size of $\epsilon=0.5$ and
$k_{\mathrm{max}}=6$ integration steps.

Figure 8 shows the ergodicity measures for EHMC, RMHMC, and LMC. We observe
that each of these MCMC algorithms exhibit a linear decrease in
$|{\mathrm{MMD}_{u}^{2}}|$ on a logarithmic scale. In fig. 9 we visualize the
ESJD, the minimum ESS, and the minimum ESS per second for the three MCMC
algorithms. We observe that LMC exhibits the strongest performance in terms of
the minimum ESS per second, whereas RMHMC only approaches the time-normalized
performance of EHMC due to its complexity.

###  4.8 Multi-Scale Student Distribution

(a) EHMC

(b) LMRMHMC

(c) LMLMC

Figure 10: We examine the ergodicity of the Markov chains as a function of the
sampling step in the multiscale Student-$t$ distribution. Here we see that, in
the case of LMC, it is beneficial from an ergodicity perspective to mix
heavily with MMALA, whereas for RMHMC, only modest mixing probabilities can
produce a modified Markov chain that is competitive with the unmodified
version. The Euclidean HMC struggles in this posterior distribution due to the
multiple spatial scales.

(a) ESJD

(b) MSJD

Figure 11: We show the ESJD and the MSJD for inference in the multiscale
Student-$t$ distribution. Curiously, EHMC enjoys the largest ESJD, but this
does not correspond to more efficient sampling, as shown in the inferior
distribution of KS statistics. To examine this property further, we also
measured the MSJD, which reveals that EHMC is less effective in traversing the
parameter space compared to the RMHMC geometric method.

(a) KS

(b) Min. ESS / Sec.

Figure 12: We show the distribution of Kolmogorov-Smirnov statistics and the
time-normalized minimum ESS for inference in the multiscale Student-$t$
distribution. For the geometric methods, we observe that RMHMC has much larger
ESJD than LMC; however, in terms of time-normalized performance, it is best to
fully mix either method with MMALA. Indeed, we find that MMALA on its own
produces a distribution of KS statistics that is competitive with any of the
mixture kernels.

Fix $m\in\mathbb{N}$ and let $\Sigma\in\mathrm{PD}(m)$ and $\nu>2$. The
density function of the multivariate Student-$t$ distribution is,

$\displaystyle\pi(q)\propto\left({1+\frac{1}{\nu}x^{\top}\Sigma^{-1}x}\right)^{-(m+\nu)/2}.$
(90)

In our experiments we consider $m=20$, $\nu=5$, and
$\Sigma=\mathrm{diag}(1,\ldots,1,10^{4})\in\mathbb{R}^{m\times m}$. The
presence of severely differing spatial scales in the multivariate Student-$t$
distribution will cause Euclidean HMC with identity mass matrix to exhibit
highly oscillatory behavior (Pourzanjani and Petzold, 2019), which will limit
the efficiency of the method both in terms of ergodicity and the effective
sample size. We employ Euclidean HMC with a step-size of $\epsilon=0.8$, and
RMHMC and LMC with a step-size of $\epsilon=0.7$; in each case we employ
$k_{\mathrm{max}}=20$ integration steps. In the case of RMHMC and LMC, we
consider a metric given by the positive definite term in the Hessian of the
log-density of the multivariate Student-$t$ distribution.

In fig. 10 we show $|{\mathrm{MMD}_{u}^{2}}|$ as a function of the number of
Markov chain steps. All of the MCMC algorithms produce a linear decrease in
the estimate of the maximum mean discrepancy on a logarithmic scale; however,
the methods differ drastically in terms of the slope of this linear
relationship. We see that RMHMC exhibits by far the fastest convergence.
Notably, LMC exhibits slower convergence on this target distribution than
MMALA; this is due to the dimensionality of the posterior, in which LMC
struggles to maintain the Hamiltonian energy required to accept proposals. In
fig. 12 we show the ESJD, KS, and the minimum ESS per second. EHMC produces
the largest ESJD but this does not translate into a large ESS due to the
oscillatory behavior of the EHMC proposal mechanism, with RMHMC producing the
largest minimum ESS per second despite its computational complexity.

##  5 Conclusion

This work has considered methods by which to equip RMHMC and LMC with a
geometric ergodicity theory. The fundamental technique we adopt is to replace
the RMHMC (or LMC) transition kernel consisting of a single integration step
with the MMALA transition kernel. This modification is inspired by the
Euclidean case, in which single-step HMC and MALA can be constructed to be
exactly equivalent. By establishing reversibility of the marginal transition
kernels, geometric ergodicity can be inherited from MMALA. We evaluated the
modified variations of RMHMC and LMC, called LMRMHMC and LMLMC, respectively,
on a suite of Bayesian inference tasks. We found that aggressively mixing with
MMALA transition kernel can be detrimental for the performance of the Markov
chain on a variety of metrics, but that more modest mixing can produce
behaviors competitive with, or exceeding, the original RMHMC or LMC methods
while still imbuing the methods with a supporting theory of geometric
ergodicity.

Acknowledgments

We thank the Yale Center for Research Computing for use of the research
computing infrastructure. This material is based upon work supported by the
National Science Foundation Graduate Research Fellowship under Grant No.
1752134. Any opinion, findings, and conclusions or recommendations expressed
in this material are those of the authors(s) and do not necessarily reflect
the views of the National Science Foundation. The work is also supported in
part by NIH/NIGMS R01GM136780 and AFOSR FA9550-21-1-0317.

## Statements and Declarations

No competing interests to declare.

##  Appendix A Proofs

###  A.1 Proof of Proposition 5

Proof: The proposal distribution of MALA is

$\displaystyle\tilde{q}|q$
$\displaystyle\sim\mathrm{Normal}\left({q+\frac{\epsilon^{2}}{2}\mathbf{A}^{-1}\nabla\log\pi(q),\epsilon^{2}\mathbf{A}^{-1}}\right)$
(91) $\displaystyle\tilde{\pi}(\tilde{q}|q)$
$\displaystyle\propto\exp\left({-\frac{1}{2\epsilon^{2}}\left({\tilde{q}-q-\frac{\epsilon^{2}}{2}\mathbf{A}^{-1}\nabla\log\pi(q)}\right)^{\top}\mathbf{A}\left({\tilde{q}-q-\frac{\epsilon^{2}}{2}\mathbf{A}^{-1}\nabla\log\pi(q)}\right)}\right)$
(92)
$\displaystyle=\exp\left({-\frac{1}{2\epsilon^{2}}\left({\mathbf{A}\tilde{q}-\mathbf{A}q-\frac{\epsilon^{2}}{2}\nabla\log\pi(q)}\right)^{\top}\mathbf{A}^{-1}\left({\mathbf{A}\tilde{q}-\mathbf{A}q-\frac{\epsilon^{2}}{2}\nabla\log\pi(q)}\right)}\right).$
(93)

Therefore, the acceptance probability of MALA is,

$\displaystyle\begin{split}&\min\left\\{{1,\frac{\pi(\tilde{q})\tilde{\pi}(q|\tilde{q})}{\pi(q)\tilde{\pi}(\tilde{q}|q)}}\right\\}\\\
=&\min\left\\{{1,\frac{\pi(\tilde{q})}{\pi(q)}\cdot\frac{\exp\left({-\frac{1}{2\epsilon^{2}}\left({\mathbf{A}q-\mathbf{A}\tilde{q}-\frac{\epsilon^{2}}{2}\nabla\log\pi(\tilde{q})}\right)^{\top}\mathbf{A}^{-1}\left({\mathbf{A}q-\mathbf{A}\tilde{q}-\frac{\epsilon^{2}}{2}\nabla\log\pi(\tilde{q})}\right)}\right)}{\exp\left({-\frac{1}{2\epsilon^{2}}\left({\mathbf{A}\tilde{q}-\mathbf{A}q-\frac{\epsilon^{2}}{2}\nabla\log\pi(q)}\right)^{\top}\mathbf{A}^{-1}\left({\mathbf{A}\tilde{q}-\mathbf{A}q-\frac{\epsilon^{2}}{2}\nabla\log\pi(q)}\right)}\right)}}\right\\}.\end{split}$
(94)

The proposal $\tilde{q}$ given $q$ can be sampled by generating
$z\sim\mathrm{Normal}(0,\mathrm{Id})$ and setting

$\displaystyle\tilde{q}=q+\frac{\epsilon}{2}\mathbf{A}^{-1}\nabla\log\pi(q)+\epsilon\sqrt{\mathbf{A}^{-1}}z.$
(95)

In HMC, the leapfrog integrator is applied to the Hamiltonian
$H(q,p)=-\log\pi(q)+\frac{1}{2}p^{\top}\mathbf{A}^{-1}p$. The sequence of
updates is,

$\displaystyle\bar{p}$ $\displaystyle=p+\frac{\epsilon}{2}\nabla\log\pi(q)$
(96) $\displaystyle\tilde{q}^{\prime}$
$\displaystyle=q+\epsilon\mathbf{A}^{-1}\bar{p}$ (97)
$\displaystyle=q+\frac{\epsilon^{2}}{2}\mathbf{A}^{-1}\nabla\log\pi(q)+\epsilon\mathbf{A}^{-1}p$
(98) $\displaystyle\tilde{p}^{\prime}$
$\displaystyle=\bar{p}+\frac{\epsilon}{2}\nabla\log\pi(\tilde{q}^{\prime})$
(99)
$\displaystyle=p+\frac{\epsilon}{2}\nabla\log\pi(q)+\frac{\epsilon}{2}\nabla\log\pi(\tilde{q}^{\prime}).$
(100)

When sampling from the distribution with density proportional to
$\exp(-H(q,p))$, one sees by inspection that $p$ is independent of $q$ and
that $p\sim\mathrm{Normal}(0,\mathbf{A})$. Therefore, using the same $z$ as in
eq. 95, we can sample $p$ by setting $p=\sqrt{\mathbf{A}}z$ so that
$\epsilon\mathbf{A}^{-1}p=\epsilon\sqrt{\mathbf{A}^{-1}}z$; we therefore see
that, in this case, HMC and MALA produce exactly the same proposal (i.e.
$\tilde{q}^{\prime}=\tilde{q}$). Recall that the HMC acceptance probability
is,

$\displaystyle\min\left\\{{1,\frac{\pi(\tilde{q})}{\pi(q)}\cdot\frac{\exp\left({-\frac{1}{2}(\tilde{p}^{\prime})^{\top}\mathbf{A}^{-1}(\tilde{p}^{\prime})}\right)}{\exp\left({-\frac{1}{2}p^{\top}\mathbf{A}^{-1}p}\right)}}\right\\}.$
(101)

By rearranging eq. 97 we find,

$\displaystyle p$
$\displaystyle=\frac{1}{\epsilon}\left({\mathbf{A}\tilde{q}-\mathbf{A}q-\frac{\epsilon^{2}}{2}\nabla\log\pi(q)}\right)$
(102) $\displaystyle\tilde{p}^{\prime}$
$\displaystyle=\frac{1}{\epsilon}\left({\mathbf{A}\tilde{q}-\mathbf{A}q+\frac{\epsilon^{2}}{2}\nabla\log\pi(\tilde{q})}\right)$
(103)
$\displaystyle=-\frac{1}{\epsilon}\left({\mathbf{A}q-\mathbf{A}\tilde{q}-\frac{\epsilon^{2}}{2}\nabla\log\pi(\tilde{q})}\right).$
(104)

Substituting these into eq. 101 and comparing to LABEL:geometric-
ergodicity:eq:mala-acceptance shows that not only are $\tilde{q}$ and
$\tilde{q}^{\prime}$ identical but that the acceptance probabilities are also
identical. Therefore, the marginal chain of single-step HMC is exactly
equivalent to the MALA chain.

###  A.2 Proof of Lemma 6

Proof:

$\displaystyle\mathrm{Pr}(q^{n+1}\in Q|q^{n}=q)$
$\displaystyle=\mathrm{Pr}(q^{n+1}\in
Q~{}\mathrm{and}~{}p^{n+1}\in\mathbb{R}^{m}|q^{n}=q)$ (105)
$\displaystyle=\int_{\mathbb{R}^{m}}\mathrm{Pr}(q^{n+1}\in
Q~{}\mathrm{and}~{}p^{n+1}\in\mathbb{R}^{m}|q^{n}=q,p^{n}=p)~{}\pi(p|q)~{}\mathrm{d}p$
(106)
$\displaystyle=\int_{\mathbb{R}^{m}}K((q,p),(Q,\mathbb{R}^{m}))~{}\pi(p|q)~{}\mathrm{d}p.$
(107)

###  A.3 Proof of Proposition 7

Proof: Given $Q,Q^{\prime}\in\mathfrak{B}(\mathbb{R}^{m})$ we have

$\displaystyle\int_{Q}\tilde{K}(q,Q^{\prime})~{}\pi(q)~{}\mathrm{d}q$
$\displaystyle=\int_{Q}\int_{\mathbb{R}^{m}}K((q,p),(Q^{\prime},\mathbb{R}^{m}))~{}\pi(p|q)~{}\mathrm{d}p~{}\pi(q)~{}\mathrm{d}q$
(108)
$\displaystyle=\int_{Q}\int_{\mathbb{R}^{m}}K((q,p),(Q^{\prime},\mathbb{R}^{m}))~{}\pi(q,p)~{}\mathrm{d}p~{}\mathrm{d}q$
(109)
$\displaystyle=\int_{Q^{\prime}}\int_{\mathbb{R}^{m}}K((q,p),(Q,\mathbb{R}^{m}))~{}\pi(q,p)~{}\mathrm{d}p~{}\mathrm{d}q$
(110)
$\displaystyle=\int_{Q^{\prime}}\int_{\mathbb{R}^{m}}K((q,p),(Q,\mathbb{R}^{m}))~{}\pi(p|q)~{}\mathrm{d}p~{}\pi(q)~{}\mathrm{d}q$
(111) $\displaystyle=\int_{Q^{\prime}}\tilde{K}(q,Q)~{}\pi(q)~{}\mathrm{d}q,$
(112)

where in eq. 110 we have used the fact that the phase space chain satisfies
detailed balance with respect to $\pi(q,p)$. This verifies that the marginal
chain satisfies detailed balance.

###  A.4 Proof of Involution Composition

######  Proposition 15.

Suppose that
$\Phi:\mathbb{R}^{m}\times\mathbb{R}^{m}\to\mathbb{R}^{m}\times\mathbb{R}^{m}$
is an invertible function and let $\mathbf{F}$ be the momentum flip function
given in definition 21. Suppose that $\mathbf{F}\circ\Phi$ is an involution,
then $\mathbf{F}\circ\Phi^{k}$ is also an involution.

Proof: This will be proved by induction with the base case established by
assumption. As the inductive hypothesis, assume that $\mathbf{F}\circ\Phi^{k}$
is an involution. Using the fact that $\mathbf{F}\circ\Phi$ is an involution,
we immediately obtain that $\Phi^{-1}\circ\mathbf{F}=\mathbf{F}\circ\Phi$.
Using the inductive hypothesis, one also has
$\Phi^{-k}\circ\mathbf{F}=\mathbf{F}\circ\Phi^{k}$. Therefore, we obtain,

$\displaystyle\Phi^{-1}\circ\mathbf{F}=\mathbf{F}\circ\Phi$ (113)
$\displaystyle\implies$
$\displaystyle\Phi^{-1}\circ\mathbf{F}\circ\Phi^{k}=\mathbf{F}\circ\Phi^{k+1}$
(114) $\displaystyle\implies$
$\displaystyle\Phi^{-1}\circ\Phi^{-k}\circ\mathbf{F}=\mathbf{F}\circ\Phi^{k+1}$
(115) $\displaystyle\implies$
$\displaystyle\Phi^{-k-1}\circ\mathbf{F}=\mathbf{F}\circ\Phi^{k+1}$ (116)
$\displaystyle\implies$
$\displaystyle\mathbf{F}\circ\Phi^{k+1}\circ\mathbf{F}\circ\Phi^{k+1}=\mathrm{Id}.$
(117)

This verifies that $\mathbf{F}\circ\Phi^{k+1}$ is also an involution.

## Data Availability Statement

The datasets generated during and/or analysed during the current study are
available in the GitHub repository, https://tinyurl.com/29kz7krz.

## References

  * Amari and Nagaoka (2000) S. Amari and H. Nagaoka.  _Methods of Information Geometry_.  Translations of mathematical monographs. American Mathematical Society, 2000.  ISBN 9780821843024.  URL https://books.google.com/books?id=vc2FWSo7wLUC.
  * Betancourt (2012) M. Betancourt.  A general metric for Riemannian manifold Hamiltonian Monte Carlo, 12 2012.
  * Betancourt et al. (2014) M. Betancourt, Simon Byrne, Samuel Livingstone, and Mark Girolami.  The geometric foundations of Hamiltonian Monte Carlo.  _Stat. Sci._ , 23, 10 2014.
  * Durmus et al. (2020) Alain Durmus, Éric Moulines, and Eero Saksman.  Irreducibility and geometric ergodicity of Hamiltonian Monte Carlo.  _The Annals of Statistics_ , 48(6):3545–3564, 2020.
  * Gelman and Pasarica (2007) Andrew Gelman and Cristian Pasarica.  Adaptively scaling the Metropolis algorithm using expected squared jumped distance.  _Statistica Sinica_ , 2007.  URL http://dx.doi.org/10.2139/ssrn.1010403.
  * Girolami and Calderhead (2011) Mark Girolami and Ben Calderhead.  Riemann manifold langevin and Hamiltonian Monte Carlo methods.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 73(2):123–214, 2011.  URL https://rss.onlinelibrary.wiley.com/doi/abs/10.1111/j.1467-9868.2010.00765.x.
  * Gretton et al. (2012) Arthur Gretton, Karsten M. Borgwardt, Malte J. Rasch, Bernhard Schölkopf, and Alexander Smola.  A kernel two-sample test.  _J. Mach. Learn. Res._ , 13(null):723–773, March 2012.  ISSN 1532-4435.
  * Hairer et al. (2006) Ernst Hairer, Christian Lubich, and Gerhard Wanner.  _Geometric numerical integration_ , volume 31 of _Springer Series in Computational Mathematics_.  Springer-Verlag, Berlin, second edition, 2006.  ISBN 3-540-30663-3; 978-3-540-30663-4.  Structure-preserving algorithms for ordinary differential equations.
  * Hsu (2002) E.P. Hsu.  _Stochastic Analysis on Manifolds_.  Contemporary Mathematics. American Mathematical Society, 2002.  ISBN 9780821808023.  URL https://books.google.com/books?id=GDEPCgAAQBAJ.
  * Kumar et al. (2019) Ravin Kumar, Colin Carroll, Ari Hartikainen, and Osvaldo A. Martin.  ArviZ a unified library for exploratory analysis of Bayesian models in Python.  _The Journal of Open Source Software_ , 2019.  URL http://joss.theoj.org/papers/10.21105/joss.01143.
  * Lan et al. (2015) Shiwei Lan, Vasileios Stathopoulos, Babak Shahbaba, and Mark Girolami.  Markov chain Monte Carlo from lagrangian dynamics.  _Journal of Computational and Graphical Statistics_ , 24(2):357–378, Apr 2015.  ISSN 1537-2715.  URL http://dx.doi.org/10.1080/10618600.2014.902764.
  * Lee and Łatuszyński (2014) Anthony Lee and Krzysztof Łatuszyński.  Variance bounding and geometric ergodicity of markov chain Monte Carlo kernels for approximate Bayesian computation.  _Biometrika_ , 101(3):655–671, 2014.  ISSN 00063444.  URL http://www.jstor.org/stable/43304674.
  * Leimkuhler and Reich (2005) Benedict Leimkuhler and Sebastian Reich.  _Simulating Hamiltonian Dynamics_.  Cambridge Monographs on Applied and Computational Mathematics. Cambridge University Press, 2005.
  * Livingstone et al. (2016) Samuel Livingstone, Michael Betancourt, Simon Byrne, and Mark Girolami.  On the geometric ergodicity of Hamiltonian Monte Carlo.  _Bernoulli_ , 25, 01 2016.  10.3150/18-BEJ1083.
  * Marsden and Ratiu (2010) Jerrold E. Marsden and Tudor S. Ratiu.  _Introduction to Mechanics and Symmetry: A Basic Exposition of Classical Mechanical Systems_.  Springer Publishing Company, Incorporated, 2010.  ISBN 1441931430.
  * Mengersen and Tweedie (1996) K. L. Mengersen and R. L. Tweedie.  Rates of convergence of the Hastings and Metropolis algorithms.  _The Annals of Statistics_ , 24(1):101 – 121, 1996.  10.1214/aos/1033066201.  URL https://doi.org/10.1214/aos/1033066201.
  * Meyn and Tweedie (1993) S.P. Meyn and R.L. Tweedie.  _Markov Chains and Stochastic Stability_.  Springer-Verlag, London, 1993.
  * Neal (2003) Radford M. Neal.  Slice sampling.  _The Annals of Statistics_ , 31(3):705 – 767, 2003.  URL https://doi.org/10.1214/aos/1056562461.
  * Neal (2010a) Radford M. Neal.  MCMC using Hamiltonian dynamics.  _Handbook of Markov Chain Monte Carlo_ , 54:113–162, 2010a.
  * Neal (2010b) Radford M. Neal.  MCMC using Hamiltonian dynamics.  _Handbook of Markov Chain Monte Carlo_ , 54:113–162, 2010b.
  * Neklyudov et al. (2020) Kirill Neklyudov, Max Welling, Evgenii Egorov, and Dmitry P. Vetrov.  Involutive MCMC: a unifying framework.  In _Proceedings of the 37th International Conference on Machine Learning, ICML 2020, 13-18 July 2020, Virtual Event_ , volume 119 of _Proceedings of Machine Learning Research_ , pages 7273–7282. PMLR, 2020\.  URL http://proceedings.mlr.press/v119/neklyudov20a.html.
  * Pourzanjani and Petzold (2019) Arya A. Pourzanjani and Linda R. Petzold.  Implicit Hamiltonian Monte Carlo for sampling multiscale distributions, 2019.
  * Roberts and Tweedie (1996) Gareth O. Roberts and Richard L. Tweedie.  Exponential convergence of Langevin distributions and their discrete approximations.  _Bernoulli_ , 2(4):341 – 363, 1996.  URL https://doi.org/.
  * Roy (2020) Vivekananda Roy.  Convergence diagnostics for Markov chain Monte Carlo.  _Annual Review of Statistics and Its Application_ , 7:387–412, 2020.
  * Roy and Zhang (2023) Vivekananda Roy and Lijin Zhang.  Convergence of position-dependent MALA with application to conditional simulation in GLMMs.  _Journal of Computational and Graphical Statistics_ , 2023.  to appear.
  * Tierney (1994) Luke Tierney.  Markov Chains for Exploring Posterior Distributions.  _The Annals of Statistics_ , 22(4):1701 – 1728, 1994.  10.1214/aos/1176325750.  URL https://doi.org/10.1214/aos/1176325750.
  * Vehtari et al. (2021) Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and Paul-Christian Bürkner.  Rank-normalization, folding, and localization: An improved rˆ for assessing convergence of mcmc (with discussion).  _Bayesian Analysis_ , 16(2), Jun 2021.  ISSN 1936-0975.  10.1214/20-ba1221.  URL http://dx.doi.org/10.1214/20-BA1221.
  * Xifara et al. (2014) T. Xifara, C. Sherlock, S. Livingstone, S. Byrne, and M. Girolami.  Langevin diffusions and the Metropolis-adjusted langevin algorithm.  _Statistics & Probability Letters_, 91:14 – 19, 2014\.  ISSN 0167-7152.  URL http://www.sciencedirect.com/science/article/pii/S0167715214001333.