# Implicit Langevin Algorithms
for Sampling From Log-concave Densities

Liam Hodgkinson School of Mathematics and Physics, University of Queensland,
Australia. Email<EMAIL_ADDRESS>Robert Salomone Australian School
of Business, and School of Mathematics and Statistics, University of New South
Wales, Australia. Email<EMAIL_ADDRESS>Fred Roosta School of
Mathematics and Physics, University of Queensland, Australia, and
International Computer Science Institute, Berkeley, USA. Email:
<EMAIL_ADDRESS>
###### Abstract

For sampling from a log-concave density, we study implicit integrators
resulting from $\theta$-method discretization of the overdamped Langevin
diffusion stochastic differential equation. Theoretical and algorithmic
properties of the resulting sampling methods for $\theta\in[0,1]$ and a range
of step sizes are established. Our results generalize and extend prior works
in several directions. In particular, for $\theta\geq 1/2$, we prove geometric
ergodicity and stability of the resulting methods for all step sizes. We show
that obtaining subsequent samples amounts to solving a strongly-convex
optimization problem, which is readily achievable using one of numerous
existing methods. Numerical examples supporting our theoretical analysis are
also presented.

##  1 Introduction

Effectively sampling from arbitrary unnormalized probability distributions is
a fundamental aspect of the Monte Carlo method, and is central in Bayesian
inference. The most common cases involve probability densities $\pi$ with
support on all of $\mathbb{R}^{d}$, which can be written in the unnormalized
form as

$\displaystyle\pi(\bm{x})\propto\exp(-f(\bm{x})),\qquad\bm{x}\in\mathbb{R}^{d}.$

The sampling problem concerns the construction of a set of points
$\\{\bm{X}_{k}\\}$ whose empirical distribution approaches $\pi$ in some
appropriate sense. A standard approach is _Markov Chain Monte Carlo_ (MCMC),
in which approximate sampling from $\pi$ is accomplished by simulating a
$\pi$-ergodic Markov chain. By the ergodic theorem, this provides consistent
Monte Carlo estimators for expectations involving the density $\pi$. The most
popular approach to generate such a set of points is the _Metropolis-Hastings_
algorithm (Hastings, 1970), which constructs a $\pi$-ergodic Markov chain by
generating a proposal from a given transition kernel and implements an
acceptance criterion for these proposals (see Robert and Casella (1999) for an
overview of such methods). While geometric rates of convergence (geometric
ergodicity) can be guaranteed in a wide variety of settings, performance is
highly susceptible to the underlying proposal. However, the effectiveness of
Metropolis-Hastings methods diminish in higher dimensions, as step sizes must
be scaled inversely with dimension, making rapid exploration of the space
unlikely (see for example, Roberts and Rosenthal (2001)).

Many of these issues lie with the steadfast requirement of consistency: that
the sample empirical distribution should _asymptotically_ be the same as
$\pi$. Ensuring this requirement in turn can result in incurring serious
penalty to the mixing rate of the chain. However, when seeking a fixed
(finite) number of samples, which is almost always the case in practice,
consistency is not necessarily a decisive property. Therefore, it has recently
become popular to consider rapidly converging Markov chains whose stationary
distributions are only _approximations_ to $\pi$ with a bias of adjustable
size; see for example (Dalalyan, 2017a, b; Wibisono, 2018; Cheng et al., 2018;
Cheng and Bartlett, 2018). While the resulting Monte Carlo estimator is no
longer consistent, it will often have dramatically smaller variance. This is
an example of a _bias-variance trade-off_ , where a biased method can require
significantly less computational effort to reach the same mean-squared error
as an asymptotically unbiased Metropolis-Hastings chain (Korattikara et al.,
2014).

The most studied of these methods is the _unadjusted Langevin algorithm_
(ULA), seen in Roberts and Tweedie (1996), which is constructed by considering
the overdamped Langevin diffusion equation, given by the stochastic
differential equation (SDE)

$\displaystyle\bm{L}_{0}\sim\pi_{0},\qquad\mathrm{d}\bm{L}_{t}=-\frac{1}{2}\nabla
f(\bm{L}_{t})\mathrm{d}t+\mathrm{d}\bm{W}_{t},$ (1)

and employing the forward Euler integrator, also known as Euler–Maruyama
approximation (Kloeden and Platen, 2013), to obtain iterates of the form

$\displaystyle\bm{X}_{k+1}=\bm{X}_{k}-\frac{h}{2}\nabla
f(\bm{X}_{k})+\sqrt{h}\bm{Z}_{k}.$ (2)

Here, $\bm{W}_{t}$ denotes $d$-dimensional standard Brownian motion, $\pi_{0}$
is some arbitrary (possibly deterministic) initial distribution, each
$\bm{Z}_{k}$ is an independent standard $d$-dimensional normal random vector,
and $h$ is the _step size_ parameter representing the temporal mesh size of
the Euler method. Since 2 is explicitly defined, it is often referred to as
_explicit_ Euler scheme. This main interest in 1 lies in the well known fact
that, under certain mild conditions and regardless of $\pi_{0}$, for any $t>0$
the distribution of $\bm{L}_{t}$ is absolutely continuous (so we may consider
its density $\pi_{t}$ on $\mathbb{R}^{d}$), and $\bm{L}_{t}$ is an ergodic
Markov process with limiting distribution $\pi$, that is,
$\pi_{t}(\bm{x})\to\pi(\bm{x})$ as $t\to\infty$ for all
$\bm{x}\in\mathbb{R}^{d}$ (Kolmogorov, 1937).

However, unlike the Langevin SDE, the distribution of samples obtained from
ULA 2 will, generally speaking, not converge to $\pi$ as $t\to\infty$. More
precisely, ULA is an asymptotically _biased_ sampling algorithm, with
corresponding bias proportional to step size (temporal mesh size). Despite
this, in situations where MCMC fails to perform well, for example, high-
dimensional problems, ULA can provide approximate samples from the target
density with acceptable accuracy (Durmus and Moulines, 2016).

The theoretical properties of ULA, including geometric ergodicity (Hansen,
2003; Roberts and Tweedie, 1996), and performance in high dimensions (Durmus
and Moulines, 2016) are well understood. Of particular relevance to us is the
recent work of Dalalyan (2017a), Dalalyan (2017b), and Durmus and Moulines
(2017), concerning the stability of ULA. Although it does not possess a single
technical definition, stability of stochastic processes is often well-
understood conceptually — some common characterizations include non-
evanescence and Harris recurrence (Meyn and Tweedie, 2012, p. 15). To
establish stability, the aforementioned works develop theoretical guarantees
in the form of error bounds on the 2-Wasserstein metric between iterates of
ULA and the target distribution. Doing so gives conditions under which ULA is
bounded in probability, which in turn implies non-evanescence (Meyn and
Tweedie, 2012, Proposition 12.1.1), and Harris recurrence, of the
corresponding Markov chain (Meyn and Tweedie, 2012, Theorem 9.2.2). The
inexact case, where $\nabla f$ is approximated to within an absolute
tolerance, is also considered (Dalalyan and Karagulyan, 2017). Some
alternative unadjusted explicit methods have also been considered; these are
usually derived using other diffusions whose stationary distributions can also
be prescribed (Cheng et al., 2017).

As a direct result of the explicit nature of the underlying discretization
scheme, the main issue with ULA-type algorithms is that they are stable only
up to a fixed step size, beyond which the chain is no longer ergodic. In fact,
Roberts and Tweedie (1996) actively discourage the use of ULA for this reason,
and show that ULA may be transient for large step sizes. Stability is an
essential concept when designing and analyzing methods for the numerical
integration of continuous-time differential equations (Ascher, 2008). In some
cases, this step size must be taken extremely small to remain stable. This
becomes a major hindrance to the performance of the method in practice.
Drawing comparisons to the theory of ordinary differential equations (ODEs) by
dropping the stochastic term, in these cases, the Langevin diffusion is said
to be _stiff_ (Ascher and Petzold, 1998). Ill-conditioned problems, such as
sampling from any multivariate normal distribution with a covariance matrix
possessing a large condition number (Golub and Van Loan, 2012, §2.6.2), are
likely to induce a stiff Langevin diffusion (Lambert, 1991, §6.2). The
negative side-effects associated with ill-conditioning as well as the
restrictions on step size are often only exacerbated in high-dimensional
problems.

In this light, a natural alternative to using explicit schemes with careful
choice of step size is to consider _implicit_ variants. From the established
theory of numerical solutions of ODEs, it is well-known that implicit
integrators have larger regions of stability than explicit alternatives, that
is, one can take larger steps without unboundedly magnifying the underlying
discretization errors (Ascher, 2008). Motivated by this, we can instead
consider the $\theta$-method scheme (Ascher, 2008, p. 84), which when applied
to Langevin dynamics 1, yields general iterations of form

$\displaystyle\bm{X}_{k+1}=\bm{X}_{k}-\frac{h}{2}\big{[}\theta\nabla
f(\bm{X}_{k+1})+(1-\theta)\nabla f(\bm{X}_{k})\big{]}+\sqrt{h}\bm{Z}_{k},$ (3)

for some $\theta=[0,1]$. The special cases of $\theta=0,1$ and $1/2$
correspond to forward, backward, and trapezoidal integrators, respectively. Of
course, for $\theta=0$, 3 reduces to the explicit Euler scheme 2. As the
choice of $\theta\in(0,1]$ define the endpoint in an implicit way, such
integrators are often referred to as implicit.

To our knowledge, there have only been a handful of efforts to study the
properties of sampling algorithms obtained from such implicit schemes. A
universal analysis of sampling schemes based on general Langevin diffusion was
conducted in Mattingly et al. (2002). There, it was shown that the implicit
Euler scheme, and other numerical methods satisfying a certain minorization
condition are geometrically ergodic for sufficiently small step sizes under
the assumption that $f$ is ‘essentially quadratic’. In a more focused
analysis, Casella et al. (2011) investigated the ergodic properties of a few
implicit schemes (including the $\theta$-method scheme) for a restricted
family of one-dimensional super-Gaussian target densities. They found that, in
this setting, the $\theta$-method results in a geometrically ergodic chain for
_any_ step size $h>0$, provided that $\theta\geq 1/2$, and suggested the same
might be true in higher dimensions. Under slightly weaker assumptions than
strong convexity, Kopec (2014) conducted a weak backward error analysis
providing error bounds on the expectation of the fully implicit Euler scheme
($\theta=1$) with respect to suitable test functions. More recently, Wibisono
(2018) considered the $\theta=1/2$ case and provided a rate of convergence of
the scheme towards its biased stationary distribution under the 2-Wasserstein
metric, assuming strong convexity and small step sizes. Despite these efforts,
it is still unclear how implicit schemes compare with explicit schemes more
generally for large step sizes, and what the effect of $\theta$ is on the bias
of the method.

The aim of this work is to study the $\theta$-method sampling scheme 3 for all
$\theta\in(0,1]$, as it applies to the relevant case of strongly log-concave
distributions (that is, where $f$ is a strongly convex function). Such
distributions arise frequently in Bayesian regression problems (Bishop and
Tipping, 2003), for example generalized linear models (GLMs) with a Gaussian
prior (Chatfield et al., 2010).

##### Contributions.

To those ends, the contributions of this work are as follows:

  1. 1.

We show that the transition density associated with 3 has a closed form
solution. Then, using this, we establish conditions for geometric ergodicity,
in terms of $\theta$, the step size $h$, Lipschitz continuity, tail behavior,
and semi-convexity of $f$ (Theorem 1). By doing so, we show stability of the
$\theta$-method scheme in multivariate settings for _any_ step size when
$\theta\geq 1/2$ and $f$ is strongly convex, proving the conjecture by Casella
et al. (2011).

  2. 2.

We provide non-asymptotic theoretical guarantees for long-time behavior of 3,
which extend those of Dalalyan and Karagulyan (2017, Theorem 1) to the general
implicit case (Theorem 2).

  3. 3.

As for $\theta>0$, iterations of 3 involve solving a non-linear equation, we
study the effect of inexact solutions of the underlying sub-problems. We
propose practically computable termination criteria and quantify the effect of
approximating each iterate on the convergence rate and long-term bias of the
chain with this criteria.

  4. 4.

We establish large step size asymptotics for $\theta>0$ via a central limit
theorem as $h\to\infty$ (Theorem 3). As a consequence, we develop an effective
default heuristic choice of step size.

  5. 5.

Finally, we demonstrate the empirical performance of the implicit
$\theta$-method scheme in a series of numerical experiments; namely sampling
from high-dimensional Gaussian distributions, and the posterior density of a
Bayesian logistic regression problem involving a real data set.

Proofs of all results can be found in Appendix A.

##### Notation.

In the sequel, vectors and matrices are denoted by bold lowercase and
Romanized bold uppercase letters, for example, $\bm{v}$ and $\mathbf{V}$,
respectively. We denote the identity matrix by $\mathbf{I}$. Regular lower-
case and upper-case letters, such as s $m$ and $M$, are used to denote scalar
constants. Random vectors are denoted by italicized bold uppercase letters,
such as $\bm{X}$. For two symmetric matrices $\mathbf{A}$ and $\mathbf{B}$,
$\mathbf{A}\succeq\mathbf{B}$ indicates that $\mathbf{A}-\mathbf{B}$ is
symmetric positive semi-definite. For vectors, we let $\|\cdot\|$ denote the
Euclidean norm of appropriate dimension, and $\|\cdot\|_{L^{2}}$ denote the
$L^{2}$ norm acting on random vectors, that is,
$\|\bm{X}\|_{L^{2}}^{2}\coloneqq\mathbb{E}\|\bm{X}\|^{2}$. For matrices,
$\|\,\cdot\,\|_{2}$ denotes the spectral norm.

##  2 Implicit Langevin Algorithm (ILA)

In this section, we establish conditions under which the sequence of
$\theta$-method iterates 3 form a Markov chain that is geometrically ergodic.
For this, we impose the following assumption on the smoothness of $f$, which
ensures the existence of a unique solution to 1; see Ikeda and Watanabe (2014,
Theorem 2.4–3.1).

######  Assumption 1.

The function $f\in\mathcal{C}^{2}$ (it is twice continuously differentiable),
and $\nabla f$ is $M$-Lipschitz, that is,

$\displaystyle\|\nabla f(\bm{x})-\nabla f(\bm{y})\|\leq
M\|\bm{x}-\bm{y}\|,\qquad\mbox{for any }\bm{x},\bm{y}\in\mathbb{R}^{d}.$

Under 1, Dalalyan (2017b) shows that if $h<4/M$, iterations of ULA 2 are
stable. However, this restriction is a fundamental disadvantage of ULA. If $M$
is particularly large, as might be the case in ill-conditioned problems and in
high dimensions where the ULA is commonly applied, then the step size must be
taken very small, which results in slow mixing time of the chain and high
autocorrelation of the samples. In sharp contrast, we now show that for
appropriate choice of $\theta$, 3 does not suffer from this restriction.

The process of obtaining samples by iterating 3 is outlined in Algorithm 1.
For brevity, we henceforth refer to this procedure as the implicit Langevin
algorithm (ILA). Note that (3) can be rewritten as

$({\cal I}+\tfrac{1}{2}h\theta\nabla
f)(\bm{X}_{k+1})=\bm{X}_{k}-\tfrac{1}{2}h(1-\theta)\nabla
f(\bm{X}_{k})+\sqrt{h}\bm{Z}_{k},$

where ${\cal I}$ denotes the identity mapping. Assuming that solutions to 3
exist and are unique, that is ${\cal I}+\frac{1}{2}h\theta\nabla f$ is
globally invertible, this implies

$\bm{X}_{k+1}=({\cal I}+\tfrac{1}{2}h\theta\nabla
f)^{-1}\left[\bm{X}_{k}-\tfrac{1}{2}h(1-\theta)\nabla
f(\bm{X}_{k})+\sqrt{h}\bm{Z}_{k}\right].$ (4)

Conditions under which the procedure 4 is guaranteed to be well-defined are
discussed in §2.1. For the time being, it is also assumed that in Algorithm 1
that (3) can be solved _exactly_. The discussion of inexact solutions is
relegated to §3.

Input : \- Initial value $\bm{X}_{0}=\bm{x}_{0}\in\mathbb{R}^{d}$

\- Number of samples $n$

\- Step size $h>0$

\- $\theta$-method parameter $\theta\in(0,1]$

for  _$k=0,1,\ldots,n$_ do

Draw $\bm{Z}_{k}\sim\mathcal{N}(0,\mathbf{I})$

Solve 3 to obtain $\bm{X}_{k+1}$

end for

Algorithm 1 Implicit Langevin Algorithm (ILA)

###  2.1 Theoretical analysis of Algorithm 1

In this section, we establish sufficient conditions for the geometric
ergodicity of the sequence of iterates generated from Algorithm 1. To conduct
such an analysis, we require the transition kernel density
$p(\bm{y}\,|\,\bm{x})$ induced from (3). In general, this is only implicitly
defined, however, assuming ${\cal I}+\frac{1}{2}h\theta\nabla f$ is globally
invertible, $p(\bm{y}\,|\,\bm{x})$ is nevertheless available in closed form.
Assuming that $f\in\mathcal{C}^{2}$, this is true whenever $h$ is chosen so
that

$\mathbf{I}+\frac{h\theta}{2}\nabla^{2}f(\bm{x})\succ 0,\qquad\mbox{for all
}\bm{x}\in\mathbb{R}^{d}.$ (5)

Therefore, at the very least, for (3) to be well-defined as a sampling method,
we require $f$ to be semi-convex, that is, there exists some $\gamma>0$ such
that $\nabla^{2}f(\bm{x})+\gamma\mathbf{I}$ is positive-semidefinite for all
$\bm{x}\in\mathbb{R}^{d}$. For example, under Assumption 1, $f$ is $M$-semi-
convex and (5) holds if $h<\frac{2}{\theta M}$. This restriction on step size
can be removed entirely if $f$ is assumed to be convex.

From (4), for fixed $\bm{x}\in\mathbb{R}^{d}$, note that
$p(\,\cdot\,|\,\bm{x})$ is the probability density function of the random
variable

$\bm{Y}=(\mathcal{I}+\tfrac{1}{2}h\theta\nabla
f)^{-1}(\bm{x}-\tfrac{1}{2}h(1-\theta)\nabla f(\bm{x})+\sqrt{h}\bm{Z}),$

where $\bm{Z}\sim\mathcal{N}(\bm{0},\mathbf{I})$. As this is an invertible
transformation of a standard Gaussian random vector, by the change of
variables theorem (see for example, (Shao, 2008, Proposition 1.8)), we have

$p(\bm{y}\mid\bm{x})=\left|\det\left(\mathbf{I}+\frac{h\theta}{2}\nabla^{2}f(\bm{y})\right)\right|\phi\left(\bm{y}+\frac{h\theta}{2}\nabla
f(\bm{y})\,;\,\bm{x}-\frac{h(1-\theta)}{2}\nabla
f(\bm{x}),\;h\mathbf{I}\right),$ (6)

where $\phi(\,\cdot\,;\bm{\mu},\Sigma)$ is the density of a multivariate
Gaussian distribution with mean $\bm{\mu}$ and covariance matrix $\Sigma$, and
‘$\det$’ denotes the determinant of a matrix.

It can be seen from (6) that increasing $\theta$ (or $h$ when $\theta>0$)
alters the landscape of the transition density, and the shape of its level-
sets. To illustrate this, Figure 1 depicts the contour plots of the transition
kernel 6 for an anisotropic example problem with differing $\theta$ and
initial state for the same step size. It can be seen that the case with
$\theta=0$ (ULA) results in an isotropic proposal in all situations, whereas
other choices of $\theta$ (implicit methods) yield proposal densities that can
better adapt to the anisotropic target density.

Figure 1: Transition kernel for an example density and large step size on a
two–dimensional Bayesian Logistic Regression example for large $h$. While the
traditional $\theta=0$ case (ULA) imposes isotropic proposals, for other
choices of $\theta$, the proposal density adapts to the anisotropic target
density.

With an explicit expression for the transition density, we can investigate the
stability of the iterates given by Algorithm 1. The most convenient way of
doing this is by demonstrating geometric ergodicity of the chain induced by
the transition kernel (6). Recall that a Markov chain with $n$-fold transition
kernel $p_{n}(\bm{y}|\bm{x})$ is said to be _geometrically ergodic_ toward an
invariant density $\nu(\cdot)$ if there exist constants $C>0$ and $0<\rho<1$
such that

$\displaystyle\sup_{\bm{x}\in\mathbb{R}^{d}}\int_{\mathbb{R}^{d}}|p_{n}(\bm{y}\,|\,\bm{x})-\nu(\bm{y})|\mathrm{d}\bm{y}\leq
C\rho^{n},\qquad\mbox{for all }n=1,2,\dots.$

Similarly, a diffusion process $\bm{X}_{t}$ with transition kernel density
$p_{t}(\cdot\,|\,\cdot)$ is said to be _exponentially ergodic_ towards an
invariant density $\nu(\bm{x})$ if there exist constants $C,\lambda>0$ such
that

$\displaystyle\sup_{\bm{x}\in\mathbb{R}^{d}}\int_{\mathbb{R}^{d}}|p_{t}(\bm{y}\,|\,\bm{x})-\nu(\bm{y})|\mathrm{d}\bm{y}\leq
Ce^{-\lambda t},\qquad\mbox{for all }t>0.$

It was shown in (Hansen, 2003, Eqn. (12)) that the overdamped Langevin
diffusion (1) is exponentially ergodic provided the following assumption on
$f$ holds:

######  Assumption 2.

There exists a constant $0<m<\infty$ such that

$\displaystyle
m\coloneqq\liminf_{\|\bm{x}\|\to\infty}\displaystyle\frac{\langle\nabla
f(\bm{x}),\bm{x}\rangle}{\|\bm{x}\|^{2}}.$

Intuitively, 2 imposes super-Gaussian tails of the target distribution. Under
Assumptions 1 and 2, for any $\bm{x}$ and $\bm{y}$, there exists a constant
$c(\bm{y})\geq 0$ depending on $\bm{y}$, such that (Lemma 1)

$\langle\nabla f(\bm{x})-\nabla f(\bm{y}),\,\bm{x}-\bm{y}\rangle\geq
m\|\bm{x}-\bm{y}\|^{2}-c(\bm{y}),\qquad\mbox{for every
}\bm{x}\in\mathbb{R}^{d}.$ (7)

2 is not new, having also appeared in Kopec (2014), and appears to be among
the weakest assumptions one can make to effectively study these implicit
schemes. Clearly, 7 is a significantly weaker condition than strong convexity:

######  Assumption 3.

The function $f\in\mathcal{C}^{2}(\mathbb{R}^{d})$ is $m$-strongly convex,
that is, there exists $0<m<\infty$ such that

$\langle\nabla f(\bm{x})-\nabla f(\bm{y}),\bm{x}-\bm{y}\rangle\geq
m\|\bm{x}-\bm{y}\|^{2},\qquad\mbox{for any }\bm{x},\bm{y}\in\mathbb{R}^{d}.$

In fact, it is straightforward to show that the values of $m$ in Assumptions 2
and 3 must coincide. Furthermore, we remark that under Assumptions 1 and 3,
the spectrum of every Hessian matrix $\nabla^{2}f(\bm{x})$ is controlled to be
within $[m,M]\subset(0,\infty)$. Strong convexity is quite a natural
assumption in Bayesian regression problems, as it can be guaranteed for the
class of GLMs with Gaussian priors (DasGupta, 2011).

Under Assumptions 1 and 2, we can now establish geometric ergodicity of the
$\theta$-method scheme under certain conditions on $\theta$ and $h$ (Theorem
1).

######  Theorem 1.

For $f$ satisfying Assumptions 1 and 2 that is $\gamma$-semi-convex, the
iterates of the $\theta$-method scheme with associated transition kernel (6)
form a geometrically ergodic chain when $h<\frac{2}{\theta\gamma}$ provided we
also have either $\theta\geq 1/2$, or both $\theta<1/2$ and
$h<4m/[M^{2}(1-2\theta)]$.

While Theorem 1 establishes the geometric ergodicity of the chain towards
_some_ stationary distribution, in general, that distribution need not
necessarily be $\pi$. Nevertheless, under Assumptions 1 and 3, we have
established that the $\theta$-method discretization of the overdamped Langevin
equation is stable for _any_ step size, provided $\theta\geq 1/2$. For these
choices of $\theta$, this implies that ILA is less strict about step size
tuning than ULA. As will be seen in §5, this will prove to have a profound
effect on the performance of ILA relative to ULA on high-dimensional problems.

###  2.2 Asymptotic exactness for the normal distribution

Among all values for $\theta$, we draw special attention to the choice
$\theta=1/2$. This resulting integrator, also known as the trapezoidal scheme,
is known to be second-order accurate when applied to ODEs; that is, for a
quadratic function $F$, iterates of the trapezoidal scheme for solving
$y^{\prime}=F(y)$ yield points of the exact solution (Süli and Mayers, 2003,
§12.4). An important consequence of this is that the global error incurred in
the trapezoidal scheme is $\mathcal{O}(h^{2})$ as $h\to 0$. Unfortunately, as
a consequence of the Itô calculus, this property does not hold for numerical
solutions of stochastic differential equations. The construction of schemes
with $\mathcal{O}(h^{2})$ global error generally requires careful treatment of
the stochastic term — see for example Anderson and Mattingly (2009). However,
the notion of second-order accuracy itself carries over in a rather remarkable
way.

The case where $f$ is a quadratic form corresponds to sampling from a
(multivariate) normal distribution $\mathcal{N}(\bm{\mu},\Sigma)$, where

$f(\bm{x})=\tfrac{1}{2}(\bm{x}-\bm{\mu})^{\top}\Sigma^{-1}(\bm{x}-\bm{\mu}).$
(8)

It is easy to see that, in this particular setting, 3 becomes explicitly
solvable. Indeed, letting $\mathbf{Q}=\Sigma^{-1}$, we see that $\nabla
f(\bm{x})=\mathbf{Q}(\bm{x}-\bm{\mu})$, and so

$\bm{X}_{k+1}=\left(\mathbf{I}+\frac{h\theta}{2}\mathbf{Q}\right)^{-1}\left[\left(\mathbf{I}-\frac{h(1-\theta)}{2}\mathbf{Q}\right)(\bm{X}_{k}-\bm{\mu})+\sqrt{h}\bm{Z}_{k}\right]+\bm{\mu}.$
(9)

Observe that if $\bm{X}_{0}$ is chosen to be a fixed value, all of the
iterates $\bm{X}_{k}$ are normally distributed. As a consequence of Lévy
continuity, the stationary distribution of the ILA, if it exists, must also be
normally distributed. In particular, due to (9), it must have mean $\bm{m}$
and covariance $\mathbf{V}$ satisfying

$\displaystyle\bm{m}-\bm{\mu}$
$\displaystyle=(\mathbf{I}+\tfrac{1}{2}h\theta\mathbf{Q})^{-1}(\mathbf{I}-\tfrac{1}{2}h(1-\theta)\mathbf{Q})(\bm{m}-\bm{\mu})$
$\displaystyle\mathbf{V}$
$\displaystyle=(\mathbf{I}+\tfrac{1}{2}h\theta\mathbf{Q})^{-2}[(\mathbf{I}-\tfrac{1}{2}h(1-\theta)\mathbf{Q})^{2}\mathbf{V}+h\mathbf{I}].$

Since $\mathbf{Q}\neq\mathbf{0}$, it must be the case that $\bm{m}=\bm{\mu}$.
Solving for $\mathbf{V}$, the stationary distribution of the ILA is found to
be

$\mathcal{N}\left(\bm{\mu},\Sigma\left(\mathbf{I}+\tfrac{1}{2}h(\theta-\tfrac{1}{2})\mathbf{Q}\right)^{-1}\right).$

Here we encounter the remarkable fact that when $f$ is quadratic and
$\theta=1/2$, regardless of the step size chosen, ILA is _asymptotically
unbiased!_ To our knowledge, this was first observed in Wibisono (2018),
however, as a consequence of our analysis, we can now deduce that $\theta=1/2$
is the only choice of $\theta$ that yields this property. While asymptotic
exactness is unlikely to hold for other sampling problems, it suggests that
cases involving approximately quadratic $f$ should see near optimal
performance when $\theta=1/2$.

##  3 Inexact Implicit Langevin Algorithm (i-ILA)

It is clear that the utility of ILA is dependent on the solvability of 3.
Fortunately, this is made feasible by a useful reinterpretation of solutions
to (3) as those of a corresponding optimization problem. Indeed, the inverse
operator $({\cal I}+\tfrac{1}{2}h\theta\nabla f)^{-1}$ is quite commonly
considered in convex optimization, as it is equivalent to the proximal
operator $\bm{{\rm prox}}_{\frac{1}{2}h\theta f}$ defined by

$\bm{{\rm
prox}}_{f}(\bm{v})=\operatorname*{arg\,min}_{\bm{x}\in\mathbb{R}^{d}}\left\\{f(\bm{x})+\frac{1}{2}\|\bm{x}-\bm{v}\|^{2}\right\\},\qquad\bm{v}\in\mathbb{R}^{d}.$

This equivalence follows from that of an optimization problem and the root-
finding problem for its critical values (Parikh and Boyd, 2014, Eqn. (3.4)).
Therefore, (3) can be formulated as the following optimization problem (after
rescaling by $2/h$):

$\displaystyle\bm{X}_{k+1}$
$\displaystyle=\operatorname*{arg\,min}_{\bm{x}\in\mathbb{R}^{d}}\;F(\bm{x};\bm{X}_{k},\bm{Z}_{k}),$
(10a) where $\displaystyle F(\bm{x};\bm{y},\bm{z})\coloneqq\theta
f(\bm{x})+\frac{1}{h}\left\lVert\bm{x}-\left(\bm{y}-\frac{h(1-\theta)}{2}\nabla
f(\bm{y})+\sqrt{h}\bm{z}\right)\right\rVert^{2}.$ (10b)

This reinterpretation of (3) was also noted in Wibisono (2018), although only
the $\theta=1/2$ case was considered. Iterations of the form (10) are often
referred to as _proximal-point methods_ in the optimization literature
(Combettes and Pesquet, 2011; Parikh and Boyd, 2014). We remark that proximal
operators were used in Pereyra (2016) in the construction of a _proximal
unadjusted Langevin algorithm_ (P-ULA). In fact, iterates of their P-ULA
algorithm would correspond with (10) when $\theta=1$ and if the Gaussian term
$\sqrt{h}\bm{Z}_{k}$ were to be moved outside the proximal operator. As one
might expect, this discrepancy has a significant impact on the covariance of
each proposal as $h\to\infty$; it will be shown in Theorem 3 that the
asymptotic covariance of these proposals is generally anisotropic.

The implementation of Algorithm 1 now hinges entirely on our ability to solve
the subproblem 10. For the unadjusted Langevin algorithm where $\theta=0$,
this can be done trivially through a closed form solution. However, for
$\theta>0$, we generally have to resort to an iterative optimization scheme to
solve 10. Thus far, we have assumed that the optimization problem in 10 can be
solved exactly. However, more often than not this is infeasible, and one must
instead consider the effects of _approximate_ solutions of 10 in the overall
convergence of the chain. This results in a sampling variant, which is
henceforth referred to as i-ILA (for inexact ILA).

The most natural way of doing this is by measuring the error in the
corresponding root-finding problem (3) via the norm of the gradient of the
subproblem 10b, $\|\nabla F\|$. This is ideal because not only can it be
readily computed in practice, but also the termination criterion of many
iterative optimization algorithms is based on this norm falling below a given
tolerance; for example, see Nocedal and Wright (2006). Furthermore, efficient
algorithms for directly minimizing $\|\nabla F\|$, as a surrogate function for
optimization of $F$, have been recently proposed, which enjoy linear, that is,
geometric, convergence rates, even in the absence of smoothness or convexity
of $F$ (Roosta et al., 2018). In addition, for sampling in distributed
computational environments, such as when large-scale data cannot be stored on
a single machine, distributed variants of these surrogate optimization
algorithms have also been recently considered (Crane and Roosta, 2019). These
algorithms are particularly suitable as part of i-ILA since they are
guaranteed to rapidly and monotonically decrease $\|\nabla F\|$; recall that
$\|\nabla F\|$ need not be monotonically decreasing in optimization algorithms
that optimize $F$ directly.

With this in mind, we consider an inexact modification of Algorithm 1, shown
in Algorithm 2, for approximate sampling from $\pi$.

Input : \- Initial value $\bm{X}_{0}=\bm{x}_{0}\in\mathbb{R}^{d}$

\- Number of samples $n$

\- Step size $h>0$

\- $\theta$-method parameter $\theta\in[0,1]$

\- Sub-problem inexactness tolerance $\epsilon\geq 0$

for  _$k=0,1,\ldots,n$_ do

Draw $\bm{Z}_{k}\sim\mathcal{N}(0,\mathbf{I})$

Find $\bm{X}_{k+1}$ satisfying $\|\nabla
F(\bm{X}_{k+1};\bm{X}_{k},\bm{Z}_{k})\|\leq\epsilon$, where $F$ is defined in
(10b)

end for

Algorithm 2 Inexact Implicit Langevin Algorithm (i-ILA)

###  3.1 Theoretical analysis of Algorithm 2

The increased stability offered by Algorithm 1 has been established in Theorem
1. However, while Theorem 1 guarantees rapid convergence towards _some_
stationary distribution, closeness of the $\theta$-method iterates to the
target distribution $\pi$ and the effect of increasing $h$ on its bias as a
sampling method, has yet to be established. Furthermore, Algorithm 1 and its
guarantees given by Theorem 1 require exact solutions of the root-finding
problem 3, whereas Algorithm 2 allows for such problems to be solved only
inexactly. To address both of these problems, we devote this section to the
development of theoretical guarantees of Algorithm 2, inspired by the
techniques of Dalalyan (2017b). These guarantees come in the form of rate of
convergence estimates under the 2-Wasserstein metric, defined between two
probability measures $\nu$ and $\pi$ by

$\displaystyle
W_{2}(\nu,\pi)=\inf_{\bm{X}\sim\nu,\bm{Y}\sim\pi}\|\bm{X}-\bm{Y}\|_{2},$

where the infimum is taken over all couplings $(\bm{X},\bm{Y})$ of $\nu$ and
$\pi$, and is attained by some optimal coupling (Villani, 2008, Thm. 4.1). The
2-Wasserstein metric can be readily linked to other quantities of interest.
For example, from the Kantorovich-Rubinstein formula (Villani, 2008, Eqn.
(5.11)), for any $M$-Lipschitz function $\varphi$, we have that

$\displaystyle|\nu(\varphi)-\pi(\varphi)|\coloneqq\left|\int\varphi\,\mathrm{d}(\nu-\pi)\right|\leq
MW_{2}(\nu,\pi).$

Our guarantees will require the same assumptions on $f$ seen in Dalalyan
(2017b), that are Assumptions 1 and 3. Under these assumptions, the condition
number of $F$ in 10b can be written as

$\displaystyle\kappa_{h}\coloneqq\frac{1+\frac{1}{2}\theta
hM}{1+\frac{1}{2}\theta hm}.$ (11)

Recall that the condition number 11 encodes and summarizes the curvature (the
degree of relative flatness and steepness), of the graph of $F$. In
optimization, it is well-known that a large condition number typically amounts
to a more difficult problem to solve, and hence algorithms that can take such
contorted curvature into account (Newton-type methods, for example), are more
appropriate (Roosta-Khorasani and Mahoney, 2018; Xu et al., 2017). It is only
natural to anticipate that challenges corresponding to problem ill-
conditioning similarly carry over to sampling procedures as well. Indeed,
large ratios of $M/m$, which imply increasingly anisotropic level-sets for
$f$, can hint at more difficult sampling problems. For example, this
difficulty directly manifest itself in ill-conditioning of $F$, which in turn
results in more challenging sub-problems. Furthermore, in such situations,
taking a larger step size can only exacerbate the ill-condition of $F$. As a
result, similar to the role played by second-order methods in optimization,
one can naturally expect to see implicit methods to be more appropriate for
ill-conditioned sampling problems.

Under Assumptions 1 and 3, the discrepancy between the inexact variant of the
$\theta$-method given in Algorithm 2 and the target density $\pi$ under the
2-Wasserstein metric is described in Theorem 2.

######  Theorem 2.

Suppose $f$ satisfies Assumptions 1 and 3. Let $\theta\in(0,1]$ and let
$\nu_{t}$ denote the distribution of the iterate $\bm{X}_{t}$ obtained by
Algorithm 2, for each $t\geq 1$, starting from $\bm{X}_{0}\sim\nu_{0}$. Let
$\kappa_{h}$ be as in (11), and if $\theta<1$, let

$h^{\ast}=\frac{(\theta-\frac{1}{2})(M+m)+\sqrt{(\theta-\frac{1}{2})^{2}(M+m)^{2}+4\theta(1-\theta)mM}}{\theta(1-\theta)mM}.$
(12)

Furthermore,

  1. (i)

if $h\leq h^{\ast}$ or $\theta=1$, then let

$\rho=\frac{1-\frac{1}{2}h(1-\theta)m}{1+\frac{1}{2}h\theta
m},\quad\mbox{and}\quad C=\frac{\kappa_{h}}{\theta m};$ (13)

  2. (ii)

alternatively, if $\theta<1/2$ and $h^{\ast}<h<\frac{4}{M(1-2\theta)}$, or if
$1/2\leq\theta<1$ and $h>h^{\ast}$, then let

$\rho=\frac{\frac{1}{2}h(1-\theta)M-1}{\frac{1}{2}h\theta
M+1},\quad\mbox{and}\quad
C=\frac{\frac{1}{2}\kappa_{h}^{2}h}{2+\frac{1}{2}h(2\theta-1)M}.$ (14)

Then, for any $t\in\mathbb{N}$,

$W_{2}(\nu_{t},\pi)\leq\kappa_{h}\rho^{t}W_{2}(\nu_{0},\pi)+C\left(\epsilon+\min\left\\{2M\sqrt{hd}(2+\sqrt{hM}),4\sqrt{Md}\right\\}\right).$

Remark 1.  As $\theta\to 0$, for the transition point $h^{\ast}$ we have
$h^{\ast}\to\frac{4}{M+m}$. Moreover, at $\theta=0$ and $\epsilon=0$,
$\kappa_{h}=1$, and the bias term coincides with that in Dalalyan and
Karagulyan (2017, Theorem 1). Theorem 2 is thus a generalization of Dalalyan
and Karagulyan (2017, Theorem 1) to arbitrary $\theta\in[0,1]$ and error
$\epsilon$.

As Theorem 1 did for Algorithm 1, Theorem 2 shows that for $\theta\geq 1/2$,
Algorithm 2 is stable for all $h>0$. Theorem 2 does suggest that smaller
values of $\theta$ will achieve faster convergence rates and smaller biases
for small step sizes, although this does not appear to be the case in practice
(for example, refer to §5). Observe that, for $\theta>1/2$ and fixed $h$, the
bias term is in the order of $\mathcal{O}(M^{-1/2})$. This implies that
increasing the condition number when $m$ is bounded below (for example the
spherical Gaussian prior in Bayesian regression) results in smaller bias and
faster convergence. This is in sharp contrast to ULA whose performance
significantly degrades with increasing condition number in such settings.

Also, we would like to reiterate that, in stark contrast to what is observed
in Roberts and Rosenthal (2001) for Metropolis-Hastings algorithms, the rate
of convergence in Theorem 2 for Algorithm 2 is not dependent on the dimension
$d$ in any form other than through the appearances of $m$ and $M$. The
dimension appears in the bias term simply due to the natural expansion of the
Euclidean distance with dimension. In particular, following Durmus and
Moulines (2016), as the dependence on dimension is at most polynomial, this
lends credence to the claim that implicit Langevin methods are well-equipped
to handle high-dimensional sampling problems.

##  4 Asymptotics for large step size

While Theorem 2 provides an essential description of the behavior of Algorithm
2, the bounds presented there are tightest for smaller step sizes on the order
of $1/M$, and are less effective when $h$ is larger. Unfortunately, the most
useful applications of ILA will occur when $M$ is large, and so the small step
size ($h\to 0$) regime will not be all that relevant. Enabled by the increased
stability of ILA, we present a novel analysis of Algorithm 1 by establishing a
central limit-type theorem regarding asymptotic behavior of the iterates in
the $h\to\infty$ regime.

Before we begin with a formal analysis, we are able to obtain insight by
considering the behavior of the subproblem 10 as $h\to\infty$. For $h\gg 1$,
we have

$\frac{1}{h}\left\lVert\bm{x}-\bm{x}_{t}+\frac{h(1-\theta)}{2}\nabla
f(\bm{x}_{t})+\sqrt{h}\bm{z}_{t}\right\rVert^{2}=\,(1-\theta)\bm{x}\cdot\nabla
f(\bm{x}_{t})+\mathcal{O}(h^{-1/2})+C,$

where $C$ does not depend on $\bm{x}$, and so does not contribute to solving
(10). As a result, the iterates of the $\theta$ method in the $h\to\infty$
regime will satisfy the relations $\bm{x}_{t+1}=\bm{x}_{t}^{\theta}$, where we
let

$\nabla f(\bm{x}_{t}^{\theta})=\left(1-\frac{1}{\theta}\right)\nabla
f(\bm{x}_{t}),\qquad\bm{x}\in\mathbb{R}^{d}.$ (15)

Letting $\bm{x}^{\ast}$ denote the unique mode of $\pi$, iterating (15) gives

$\|\nabla f(\bm{x}_{t})-\nabla f(\bm{x}^{\ast})\|=\rho^{t}\|\nabla
f(\bm{x}_{0})-\nabla f(\bm{x}^{\ast})\|,$

where $\rho=\frac{1}{\theta}-1$. Under Assumption 3, we obtain

$\frac{m\rho^{t}}{M}\|\bm{x}_{0}-\bm{x}^{\ast}\|\leq\|\bm{x}_{t}-\bm{x}^{\ast}\|\leq\frac{M\rho^{t}}{m}\|\bm{x}_{0}-\bm{x}^{\ast}\|.$

Therefore, the behavior of the $\theta$-method for large $h$ is determined
according to the three regimes depicted in Table 1.

Table 1: Asymptotic behavior of iterates of (10) as $h\to\infty$ $0\leq\theta<1/2$ | $\rho>1$ |  $\|\bm{X}_{t}\|\to\infty$ (unbounded in probability)
---|---|---
$\theta=1/2$ | $\rho=1$ | iterates oscillate about the mode
$1/2<\theta\leq 1$ | $\rho<1$ |  $\bm{X}_{t}\to\bm{x}^{\ast}$ (collapse to the mode)

The $\theta<1/2$ case is clearly undesirable from a practical standpoint.
Moreover, the collapse towards the mode seen when $\theta$ is close to one
suggests enormous potential bias for large step sizes. On the other hand, the
$\theta=1/2$ case provides no damping effect whatsoever (a fact also supported
by Theorem 2), making it susceptible to rare large proposals. Based on this
preliminary analysis, for some small $\epsilon>0$, a choice of
$\theta=1/2+\epsilon$ appears to provide the safest, and potentially the most
accurate of our $\theta$-method samplers. This aligns with the rule-of-thumb
used for $\theta$-method discretization of ODEs (Ascher, 2008, p. 85). To
formally extend these characterizations to the implicit $\theta$-method scheme
3, in Theorem 3, a central limit theorem as $h\to\infty$ is obtained for a
single step of the scheme about the deterministic map
$\bm{x}\mapsto\bm{x}^{\theta}$.

######  Theorem 3.

Given any $f\in\mathcal{C}^{2}(\mathbb{R}^{d})$, consider iterations given by
10, where $\theta\in(0,1]$. Conditioned on $\bm{X}_{k}$, we have

$\displaystyle\sqrt{h}(\bm{X}_{k+1}-\bm{X}_{k}^{\theta})\underset{h\to\infty}{\overset{\mathcal{D}}{\longrightarrow}}\mathcal{N}\left(\bm{0},\frac{4}{\theta^{2}}\nabla^{2}f(\bm{X}_{k}^{\theta})^{-2}\right).$

Theorem 3 implies that as $h\to\infty$, the implicit $\theta$-method scheme
behaves similarly to a Markov chain $\\{\bm{W}_{k}\\}$ with transitions

$\bm{W}_{k+1}=\bm{W}_{k}^{\theta}+\frac{2}{\theta\sqrt{h}}\nabla^{2}f(\bm{W}_{k}^{\theta})^{-1}\bm{Z}_{k},$

whose dynamics mimic those of the map $\bm{x}\mapsto\bm{x}^{\theta}$, but with
an additional normally-distributed noise term at each step. Furthermore, the
variance of this noise term increases as the implicit component of the scheme
diminishes (taking $\theta\to 0$).

###  4.1 A heuristic choice for step size

A consequence of the proof of Theorem 3 is that the covariance
$\Sigma_{h}(\bm{x})$ of the proposal density from the transition kernel
$p(\bm{y}\,|\,\bm{x})$ behaves asymptotically as

$\Sigma_{h}(\bm{x})\approx
h\left(\mathbf{I}+\frac{h\theta}{2}\nabla^{2}f(\bm{x}^{\theta})\right)^{-2},\qquad\mbox{as
}h\to\infty.$

Conversely, it is relatively straightforward to show that

$\Sigma_{h}(\bm{x})\approx
h\left(\mathbf{I}+\frac{h\theta}{2}\nabla^{2}f(\bm{x})\right)^{-2},\qquad\mbox{as
}h\to 0.$

These two expressions coincide when $\bm{x}=\bm{x}^{\theta}=\bm{x}^{\ast}$,
where $\bm{x}^{\ast}$ denotes the mode. At this point, one might expect a
‘good’ transition kernel to resemble the Laplace approximation of the
distribution about $\bm{x}^{\ast}$, which has covariance
$\nabla^{2}f(\bm{x}^{\ast})^{-1}$. This suggests a heuristic for choosing a
good step size in practice, by taking $h$ as a solution to the one-dimensional
optimization problem

$\hat{h}_{\theta}\coloneqq\operatorname*{arg\,min}_{h\geq 0}\left\lVert
h\left(\mathbf{I}+\frac{h\theta}{2}\nabla^{2}f(\bm{x}^{\ast})\right)^{-2}-\nabla^{2}f(\bm{x}^{\ast})^{-1}\right\rVert_{E},$
(16)

where the norm $\|\cdot\|_{E}$ can be any matrix norm of choice. Solutions to
(16) can be obtained using off-the-shelf methods in univariate optimization,
such as golden section search (Cottle and Thapa, 2017, §9.5). We will show in
the next section that, for several examples, the step size obtained from 16
with Frobenius norm tends to be an effective choice in practice, especially
for $\theta=1/2$, where it reveals itself to be near optimal in all of our
experiments. For this choice of norm, 16 can be replaced by the equivalent
problem

$\displaystyle\hat{h}_{\theta}=\operatorname*{arg\,min}_{h\geq
0}\sum_{k=1}^{d}\left[h\left(1+\frac{h\theta}{2}\lambda_{k}\right)^{-2}-\frac{1}{\lambda_{k}}\right]^{2},$
(17)

where $\lambda_{1},\ldots,\lambda_{d}$ are the eigenvalues of
$\nabla^{2}f(\bm{x}^{\ast})$. One drawback is that solving (16) or (17)
require either inversion of $\nabla^{2}f(\bm{x}^{\ast}$), or knowledge of its
spectrum, respectively, both of which may be prohibitively expensive in high
dimensions. In many problems, however, it is reasonable to assume a certain
distribution of its spectrum; for example, that the eigenvalues
$\lambda_{1}\geq\cdots\geq\lambda_{d}$ of $\nabla^{2}f(\bm{x}^{\ast})$ decay
exponentially:

$\log\lambda_{k}\approx\left(1-\frac{k-1}{d-1}\right)\log
M+\frac{k-1}{d-1}\log m,\qquad k=1,\dots,d,$ (18)

where $m$ and $M$ take the place of the smallest and largest eigenvalues of
$\nabla^{2}f(\bm{x}^{\ast})$, respectively. Simplifying assumptions such as
these can be justified in problems where Hessian of $f$ is approximately low
rank, in the sense that it has a small stable rank (Roosta-Khorasani and
Ascher, 2015), and hence its spectrum decays fast. Under these assumptions,
solving 17 becomes more tractable; we shall make use of this for Figure 3 in
§5.

##  5 Numerical experiments

In this section, we evaluate the empirical performance of Algorithms 1 and 2
in high-dimensions as measured by a few discrepancy measures. Recall that the
total variation distance between any two absolutely continuous distributions
with densities $p$ and $q$ over $\mathbb{R}^{d}$ respectively is given by

$\displaystyle
d_{\mathrm{TV}}(p,q)\coloneqq\frac{1}{2}\int_{\mathbb{R}^{d}}|p(\bm{x})-q(\bm{x})|\mathrm{d}\bm{x}.$

Since the total variation metric is too difficult to directly estimate in
higher dimensions, we follow the standard approach in the literature (see for
example Durmus and Moulines (2016) and Maire et al. (2018)) and consider
instead the _mean marginal total variation_ (MMTV),

$\mathrm{MMTV}(p,q)\coloneqq\frac{1}{2d}\sum_{i=1}^{d}\int_{\mathbb{R}}|p_{i}(x)-q_{i}(x)|\mathrm{d}x,$

which we estimate as follows. First, kernel smoothing is applied to samples
for each marginal from an extended MCMC run, as well as samples obtained from
a single run of each method. The total variation between these estimated
univariate densities is then computed with high accuracy via Gauss-Kronrod
quadrature (Kahaner et al., 1989).

As a weakness of MMTV is its inability to adequately compare the covariances
within coordinates between the two sample sets, we also compare with a second
discrepancy measure; _maximum mean discrepancy_ (MMD) (Gretton et al., 2012;
Muandet et al., 2017). Letting $\mathcal{H}$ denote a reproducing kernel
Hilbert space over $\mathbb{R}^{d}$ with reproducing kernel $k$, MMD is
defined as the integral probability metric

$\displaystyle\mathrm{MMD}^{2}(p,q)$
$\displaystyle\coloneqq\left[\sup_{\|h\|_{\mathcal{H}}\leq
1}\int_{\mathbb{R}^{d}}h(x)[p(x)-q(x)]\mathrm{d}x\right]^{2}$
$\displaystyle=\mathbb{E}_{p,p}k(X,\tilde{X})-2\mathbb{E}_{p,q}k(X,Y)+\mathbb{E}_{q,q}k(Y,\tilde{Y}),$

where $X,\tilde{X}$ are independent random variables with distribution $p$,
and $Y,\tilde{Y}$ are independent random variables with distribution $q$.
These expectations can be estimated using samples from $p$ and $q$. In our
experiments, we use the Gaussian kernel:

$k(\bm{x},\bm{y})=\exp\left(-\frac{1}{2\sigma^{2}}\|\bm{x}-\bm{y}\|^{2}\right),$

where the kernel bandwidth parameter $\sigma$ is chosen so that $2\sigma^{2}$
is the median of $(\|\bm{x}_{i}~{}-~{}\bm{x}_{j}\|)_{i,j=1}^{n}$, where
$\bm{x}_{i}$ denotes the $i$-th sample taken from $q$.

###  5.1 High-dimensional Gaussian distributions

To highlight the effects of problem ill-conditioning, we once again consider
sampling from a multivariate Gaussian distribution, as in (8), with explicitly
computable iterates given by (9). It is easy to see that $f$ satisfies
Assumptions 1 and 3. To show efficacy in higher dimensions, we will consider
$d=1000$. Furthermore, to test the effects of ill-conditioning, we focus on
three choices of $\Sigma$ with condition numbers
$\kappa\in\\{1,100,10^{8}\\}$. Each $\Sigma$ is generated using the method of
Bendel and Mickey (1978) to uniformly sample a correlation matrix with
eigenvalues given by (18) for $m=1$ and $M=\kappa$. For simplicity, we take
$\bm{\mu}=\bm{0}$. MMTV and MMD discrepancies were computed between $\pi$ and
samples of $N=5000$ points generated by Algorithm 1 with
$\theta\in\\{0,1/2,1\\}$ and a variety of step sizes $h$ (encompassing $4/M$
and the step size heuristics in §4). Common random numbers were used, and no
burn-in period was applied. The results are shown in Figure 2.

| Mean Marginal Total Variation |  Maximum Mean Discrepancy
---|---|---

$\kappa=1\,$: |
---|---
$\kappa=100\,$: |
$\kappa=10^{8}\,$: |

Figure 2: MMTV estimates and MMD discrepancies for 5000 samples generated by
Algorithm 1 with $\theta\in\\{0,1/2,1\\}$,
$h\in[\frac{4}{100M},100\hat{h}_{1/2}]$, and target distribution $\pi$ given
by (8) with $\bm{\mu}=\bm{0}$ and $\Sigma$ a correlation matrix with condition
number $\kappa\in\\{1,10^{2},10^{8}\\}$.

Due to the rapid explosion in magnitude of samples generated by ULA when
$h\geq 4/M$, we only display discrepancies for ULA for $h<4/M$. This critical
value of is highlighted as a black solid vertical line.

In light of the large step size asymptotics, the existence of an “optimal”
step size for $\theta\geq 1/2$ as evidenced in these plots is perhaps not too
surprising. However, it is surprising to see that, especially for large
$\kappa$, this optimum is much greater than the maximum allowed step size of
$4/M$ for ULA. Most notable here is the greatly improved performance of the
implicit method ($\theta=1/2$) at this optimum over ULA for any allowable step
size. Moreover, in all cases, the optimal performance of ILA for $\theta=1/2$
exceeds that of the purely implicit method ($\theta=1$). These two facts are
not suggested by Theorem 2, implying that the large step size asymptotics
should indeed play a significant role in the analysis of implicit methods
moving forward.

In all cases, the step size heuristic for $\theta=1$ performs poorly,
suggesting the fully implicit case operates by a different mechanism that is
currently unknown to us. For $\kappa=1$,
$\hat{h}_{1/2}~{}=~{}4/M~{}\equiv~{}4$, which is clearly the optimal step
size, as it yields exact samples ($\bm{X}_{k+1}=\bm{Z}_{k}$). In fact,
$\theta=1/2$, $h=4$ is the only choice of $\theta$ and $h$ which results in
exact samples in this scenario. According to both estimated MMTV and MMD, the
step size heuristic is an almost optimal choice of $h$ for all $\kappa$
considered, even in high dimensions.

###  5.2 Logistic regression

We now consider sampling problems involving the Bayesian posterior densities
of generalized linear models (GLM), which have log-concave likelihood
functions, with Gaussian priors. For simplicity, and without loss of
generality, we consider radially symmetric Gaussians. For a GLM with this
choice of prior, posterior densities are of the form

$\displaystyle
f(\bm{x})=\sum_{i=1}^{n}\left(\Phi(\bm{a}_{i}^{T}\bm{x})-b_{i}\bm{a}_{i}^{T}\bm{x}\right)+\frac{\lambda}{2}\|\bm{x}\|^{2},$

where $(\bm{a}_{i},b_{i}),\;i=1,2,\cdots,n,$ are the response and covariate
pairs, $\bm{a}_{i}\in\mathbb{R}^{p}$, and the domain of $b_{i}$ depends on the
GLM. The cumulant generating function, $\Phi$, determines the type of GLM. For
example, in the case of logistic regression, $\Phi(t)=\log(1+\exp(t))$; see
McCullagh and Nelder (1989) for further details and applications. It is easy
to see that

$\displaystyle\nabla^{2}f(\bm{x})=\sum_{i=1}^{n}\bm{a}_{i}\bm{a}^{\top}_{i}\Phi^{{}^{\prime\prime}}(\bm{a}_{i}^{\top}\bm{x})+\lambda\mathbf{I}=\mathbf{A}^{\top}\mathbf{D}\mathbf{A}+\lambda\mathbf{I},$

where $\mathbf{A}\in\mathbb{R}^{n\times d}$ is a matrix whose $i$-th row is
$\bm{a}_{i}$, $\mathbf{D}$ is a diagonal matrix whose $i$-th diagonal element
is $\Phi^{{}^{\prime\prime}}(\bm{a}^{\top}\bm{x})$, and $\lambda$ is the
precision parameter of the prior. As a result, for Assumption 3, we have

$\lambda\leq m\leq
M\leq\|\mathbf{A}\|^{2}\sup_{t\in\mathbb{R}}\Phi^{{}^{\prime\prime}}(t)+\lambda.$
(19)

For our example, we consider Bayesian logistic regression in this setting,
yielding

$f(\bm{x})\varpropto\sum_{i=1}^{n}\left(\log\left(1+\exp(\bm{a}_{i}^{T}\bm{x})\right)-b_{i}\bm{a}_{i}^{T}\bm{x}\right)+\frac{\lambda}{2}||\bm{x}||^{2},$
(20)

and $\sup_{t\in\mathbb{R}}\Phi^{{}^{\prime\prime}}(t)\leq 1/4$. We use the
musk (version 1) dataset from the UCI repository (Dua and Graff, 2019), with
the prior precision parameter $\lambda=1$. These choices yield a target
distribution which is relatively ill-conditioned, whose Hessian
$\nabla^{2}f(\bm{x}^{\ast})$ about its mode $\bm{x}^{\ast}$ possesses a
condition number of $\kappa>2\times 10^{3}$. We estimate the values $m$ and
$M$ according to their lower and upper bounds in (19). For the target density
given according to (20), MMTV and MMD discrepancies were computed between
samples of $N=10000$ points generated using Algorithm 2 (with
$\theta\in\\{0,1/2,1\\}$, $\epsilon=10^{-9}$ and a variety of step sizes $h$
encompassing $4/M$ and the step size heuristics in §4 under the assumption
that eigenvalues are distributed according to (18)) and a gold standard run
comprised of 50,000 samples obtained from hand-tuned SMMALA (Girolami and
Calderhead, 2011). Due to the large difference in computation time between ULA
and ILA, for the sake of comparison, we also computed MMTV and MMD
discrepancies for the ULA algorithm using the same step sizes, now with a
thinning factor of 50. This factor was chosen so that the computation time of
ULA became roughly equivalent to the other ILA methods. Once again, common
random numbers were used, and no burn-in period was applied. The results are
shown in Figure 3, and follow a similar pattern to those found in the previous
example. The step size heuristic for $\theta=1/2$ performs admirably in this
case, yielding samples with smaller discrepancies to the gold standard run
than ULA for any reasonable step size, even when thinned to account for the
difference in computation time.

Mean Marginal TV | Maximum Mean Discrepancy
---|---
|

Figure 3: MMTV and MMD discrepancies between 10000 samples generated by
Algorithm 2 with $\theta\in\\{0,1/2,1\\}$ and ULA with a thinning factor of
50, over $h\in[\frac{4}{10M},100\hat{h}_{1/2}]$, and gold standard run, for
target distribution specified according to (20).

##  6 Conclusions

In the context of sampling from an unnormalized probability distribution, we
considered a general class of unadjusted sampling algorithms that are based on
implicit discretization of the Langevin dynamics. Unlike the traditional
Metropolis-adjusted sampling algorithms, these unadjusted methods relax the
requirement of consistency that the sample empirical distribution should
asymptotically be the same as the target distribution, and hence avoid
incurring serious penalty to the mixing rate of the chain. As a result, these
variants generate rapidly converging Markov chains whose stationary
distributions are only approximations of the target distribution with a bias
that is of adjustable size. When one seeks a fixed (finite) number of samples,
which is almost always the case in practice, this latter unadjusted view point
can offer greatly many advantages.

In this context, we focused on the class of discretization schemes generated
using $\theta$-method in the context of smooth and strongly log-concave
densities, explicitly deriving the transition kernel of the chain and
establishing the corresponding sub-problems that are formulated as
optimization problems. For smooth densities, the resulting implicit Langevin
algorithms (ILA) have been shown to be geometrically ergodic for $\theta\geq
1/2$, irrespective of the step size. We also considered inexact variants
(i-ILA) where the optimization sub-problems are solved only approximately. For
this, we established non-asymptotic convergence of the sample empirical
distribution to the target as measured by 2-Wasserstein metric, finding again
that for $\theta>1/2$, the resulting scheme is unconditionally stable for all
step sizes. Furthermore, the growth rate in the bias term, that is shown to
depend on problem’s condition number, is greatly diminished for $\theta>1/2$.
Together with our numerical experiments, this suggests that the implicit
methods are a more appropriate choice for ill-conditioned problems than
explicit schemes. Furthermore, the case $\theta=1/2$ appears to perform best
in practice, especially when paired with our default heuristic choice of step
size. The underlying reason for this is likely related to its asymptotic
exactness for the normal distribution. It was suggested in Wibisono (2018)
that the asymptotic bias of the $\theta=1/2$ case could be second-order
accurate, which would imply the increased performance we have observed.
Unfortunately, proving this claim remains an open problem.

Although enticing, extensions of these results to non-convex cases may prove
challenging due to the potential lack of unique solutions for the implicit
scheme, and the relative difficulty of non-convex optimization in general.
Nevertheless, one could find success in considering $f$ that is only strongly
convex outside of a compact region, as in Cheng et al. (2018). Furthermore,
although it has not been treated explicitly, we believe that implicit methods
should prove effective in big data problems, that is with
$f(\bm{x})=\sum_{i=1}^{n}f_{i}(\bm{x})$ and $n\gg 1$, where it might be
computationally prohibitive to evaluate $f$ or its gradient exactly. In this
regard, one can use optimization algorithms that can employ inexact oracle
information; see Roosta-Khorasani and Mahoney (2018) for example. The efficacy
of this approach would prove interesting for future research.

#### Acknowledgments

All authors have been supported by the Australian Centre of Excellence for
Mathematical and Statistical Frontiers (ACEMS) under grant number CE140100049.
Liam Hodgkinson acknowledges the support of an Australian Research Training
Program (RTP) Scholarship. Fred Roosta was partially supported by the
Australian Research Council through a Discovery Early Career Researcher Award
(DE180100923). Part of this work was done while Fred Roosta was visiting the
Simons Institute for the Theory of Computing.

## References

  * Anderson and Mattingly [2009] David F. Anderson and Jonathan C. Mattingly.  A weak trapezoidal method for a class of stochastic differential equations.  _arXiv preprint arXiv:0906.3475_ , 2009.
  * Ascher [2008] Uri M. Ascher.  _Numerical methods for evolutionary differential equations_.  SIAM, 2008.
  * Ascher and Petzold [1998] Uri M. Ascher and Linda Petzold.  _Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations_.  Other Titles in Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM, 3600 Market Street, Floor 6, Philadelphia, PA 19104), 1998.  ISBN 9781611971392.
  * Bendel and Mickey [1978] Robert B. Bendel and M. Ray Mickey.  Population correlation matrices for sampling experiments.  _Communications in Statistics-Simulation and Computation_ , 7(2):163–182, 1978.
  * Bishop and Tipping [2003] Christopher M. Bishop and Michael E. Tipping.  Bayesian regression and classification.  _Nato Science Series sub Series III Computer And Systems Sciences_ , 190:267–288, 2003.
  * Boucheron et al. [2013] Stéphane Boucheron, Gábor Lugosi, and Pascal Massart.  _Concentration inequalities: A nonasymptotic theory of independence_.  Oxford university press, 2013.
  * Casella et al. [2011] Bruno Casella, Gareth Roberts, and Osnat Stramer.  Stability of partially implicit Langevin schemes and their MCMC variants.  _Methodology and Computing in Applied Probability_ , 13(4):835–854, December 2011.
  * Chatfield et al. [2010] Chris Chatfield, Jim Zidek, and Jim Lindsey.  _An introduction to generalized linear models_.  Chapman and Hall/CRC, 2010.
  * Cheng and Bartlett [2018] Xiang Cheng and Peter Bartlett.  Convergence of Langevin MCMC in KL-divergence.  In _Proceedings of Algorithmic Learning Theory_ , volume 83, pages 186–211, 2018.
  * Cheng et al. [2017] Xiang Cheng, Niladri S. Chatterji, Peter L. Bartlett, and Michael I. Jordan.  Underdamped Langevin MCMC: A non-asymptotic analysis.  _arXiv preprint arXiv:1707.03663_ , 2017.
  * Cheng et al. [2018] Xiang Cheng, Niladri S. Chatterji, Yasin Abbasi-Yadkori, Peter L. Bartlett, and Michael I. Jordan.  Sharp Convergence Rates for Langevin Dynamics in the Nonconvex Setting.  _arXiv preprint arXiv:1805.01648_ , 2018.
  * Combettes and Pesquet [2011] Patrick L. Combettes and Jean-Christophe Pesquet.  Proximal splitting methods in signal processing.  In _Fixed-point algorithms for inverse problems in science and engineering_ , pages 185–212. Springer, 2011.
  * Cottle and Thapa [2017] Richard W. Cottle and Mukund N. Thapa.  _Linear and Nonlinear Optimization_.  International Series in Operations Research & Management Science. Springer New York, 2017.  ISBN 9781493970537.
  * Crane and Roosta [2019] Rixon Crane and Fred Roosta.  DINGO: Distributed Newton-Type Method for Gradient-Norm Optimization.  _arXiv preprint arXiv:1901.05134_ , 2019.
  * Dalalyan [2017a] Arnak S. Dalalyan.  Further and stronger analogy between sampling and optimization: Langevin Monte Carlo and gradient descent.  _arXiv preprint arXiv:1704.04752_ , 2017a.
  * Dalalyan [2017b] Arnak S. Dalalyan.  Theoretical guarantees for approximate sampling from smooth and log-concave densities.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 79(3):651–676, 2017b.
  * Dalalyan and Karagulyan [2017] Arnak S. Dalalyan and Avetik G. Karagulyan.  User-friendly guarantees for the Langevin Monte Carlo with inaccurate gradient.  _arXiv preprint arXiv:1710.00095_ , 2017.
  * DasGupta [2011] Anirban DasGupta.  _Probability for Statistics and Machine Learning: Fundamentals and Advanced Topics_.  Springer Texts in Statistics. Springer New York, 2011.  ISBN 9781441996336.
  * Dua and Graff [2019] Dheeru Dua and Casey Graff.  UCI machine learning repository, 2019.
  * Durmus and Moulines [2016] Alain Durmus and Eric Moulines.  High-dimensional Bayesian inference via the Unadjusted Langevin Algorithm.  _arXiv:1605.01559_ , May 2016.
  * Durmus and Moulines [2017] Alain Durmus and Eric Moulines.  Nonasymptotic convergence analysis for the unadjusted Langevin algorithm.  _The Annals of Applied Probability_ , 27(3):1551–1587, 2017.
  * 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.
  * Golub and Van Loan [2012] Gene H. Golub and Charles F. Van Loan.  _Matrix Computations_ , volume 3.  JHU Press, 4 edition, 2012.
  * Gretton et al. [2012] Arthur Gretton, Karsten M. Borgwardt, Malte J. Rasch, Bernhard Schölkopf, and Alexander Smola.  A kernel two-sample test.  _Journal of Machine Learning Research_ , 13(Mar):723–773, 2012.
  * Hansen [2003] Niels Richard Hansen.  Geometric ergodicity of discrete-time approximations to multivariate diffusions.  _Bernoulli_ , 9(4):725–743, 2003.
  * Hastings [1970] Wilfred K. Hastings.  Monte Carlo sampling methods using Markov chains and their applications.  _Biometrika_ , 57(1):97–109, 1970.  ISSN 00063444.
  * Ikeda and Watanabe [2014] Nobuyuki Ikeda and Shinzo Watanabe.  _Stochastic differential equations and diffusion processes_ , volume 24.  Elsevier, 2014.
  * Kahaner et al. [1989] David Kahaner, Cleve B. Moler, Stephen Nash, and George E. Forsythe.  _Numerical Methods and Software_.  Prentice-Hall series in computational mathematics. Prentice Hall, 1989\.
  * Kloeden and Platen [2013] Peter E. Kloeden and Eckhard Platen.  _Numerical solution of stochastic differential equations_ , volume 23.  Springer Science & Business Media, 2013.
  * Kolmogorov [1937] Andrei N. Kolmogorov.  Zur umkehrbarkeit der statistischen naturgesetze.  _Mathematische Annalen_ , 113(1):766–772, 1937\.
  * Kopec [2014] Marie Kopec.  Weak backward error analysis for overdamped Langevin processes.  _IMA Journal of Numerical Analysis_ , 35(2):583–614, 2014.
  * Korattikara et al. [2014] Anoop Korattikara, Yutian Chen, and Max Welling.  Austerity in MCMC land: Cutting the Metropolis-Hastings budget.  In _International Conference on Machine Learning_ , pages 181–189, 2014.
  * Lambert [1991] John Denholm Lambert.  _Numerical methods for ordinary differential systems: the initial value problem_.  John Wiley & Sons, Inc., 1991.
  * Maire et al. [2018] Florian Maire, Nial Friel, and Pierre Alquier.  Informed sub-sampling MCMC: approximate Bayesian inference for large datasets.  _Statistics and Computing_ , pages 1–34, Jun 2018.
  * Mattingly et al. [2002] Jonathan C. Mattingly, Andrew M. Stuart, and Desmond J. Higham.  Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise.  _Stochastic processes and their applications_ , 101(2):185–232, 2002.
  * McCullagh and Nelder [1989] Peter McCullagh and John A. Nelder.  _Generalized linear models_ , volume 37.  CRC press, 1989.
  * Meyn and Tweedie [2012] Sean P. Meyn and Richard L. Tweedie.  _Markov chains and stochastic stability_.  Springer Science & Business Media, 2012.
  * Muandet et al. [2017] Krikamol Muandet, Kenji Fukumizu, Bharath Sriperumbudur, and Bernhard Schölkopf.  Kernel mean embedding of distributions: A review and beyond.  _Foundations and Trends® in Machine Learning_ , 10(1-2):1–141, 2017.
  * Nocedal and Wright [2006] Jorge Nocedal and Stephen Wright.  _Numerical Optimization_.  Springer Science & Business Media, 2006.
  * Parikh and Boyd [2014] Neal Parikh and Stephen Boyd.  Proximal Algorithms.  _Foundations and Trends® in Optimization_ , 1(3):127–239, 2014.
  * Pereyra [2016] Marcelo Pereyra.  Proximal Markov Chain Monte Carlo algorithms.  _Statistics and Computing_ , 26(4):745–760, Jul 2016.  ISSN 1573-1375.
  * Robert and Casella [1999] Christian P. Robert and George Casella.  _Monte Carlo Statistical Methods_.  Springer texts in statistics. Springer, 1999.  ISBN 9780387987071.
  * Roberts and Rosenthal [2001] Gareth O. Roberts and Jeffrey S. Rosenthal.  Optimal scaling for various Metropolis-Hastings algorithms.  _Statistical science_ , 16(4):351–367, 2001.
  * 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, 12 1996.
  * Roosta et al. [2018] Fred Roosta, Yang Liu, Peng Xu, and Michael W. Mahoney.  Newton-MR: Newton’s Method Without Smoothness or Convexity.  _arXiv preprint arXiv:1810.00303_ , 2018.
  * Roosta-Khorasani and Ascher [2015] Farbod Roosta-Khorasani and Uri M. Ascher.  Improved bounds on sample size for implicit matrix trace estimators.  _Foundations of Computational Mathematics_ , 15(5):1187–1212, 2015.
  * Roosta-Khorasani and Mahoney [2018] Farbod Roosta-Khorasani and Michael W. Mahoney.  Sub-sampled Newton methods.  _Mathematical Programming_ , Nov 2018.
  * Shao [2008] Jun Shao.  _Mathematical Statistics_.  Springer Texts in Statistics. Springer New York, 2008.  ISBN 9780387217185.
  * Süli and Mayers [2003] Endre Süli and David F. Mayers.  _An introduction to numerical analysis_.  Cambridge University Press, 2003.
  * Villani [2008] Cédric Villani.  _Optimal Transport: Old and New_ , volume 338.  Springer Science & Business Media, 2008.
  * Wibisono [2018] Andre Wibisono.  Sampling as optimization in the space of measures: The Langevin dynamics as a composite optimization problem.  _arXiv preprint arXiv:1802.08089_ , 2018.
  * Xu et al. [2017] Peng Xu, Farbod Roosta-Khorasani, and Michael W. Mahoney.  Newton-Type Methods for Non-Convex Optimization Under Inexact Hessian Information.  _arXiv preprint arXiv:1708.07164_ , 2017.

##  Appendix A Proofs

In this section, we gather the detailed proofs of all the main results of this
work along with some additional technical lemmas. In particular, Sections A.1,
A.2 and A.3 give, respectively, the proofs of Theorems 1, 2 and 3.

###  A.1 Proof of Theorem 1 (Geometric Ergodicity)

To establish geometric ergodicity, we prove the stronger Proposition 1 below.
First, we connect Assumptions 1 and 2 to the lower bound (7).

######  Lemma 1.

The condition (7) holds under Assumptions 1 and 2.

###### Proof.

By 1, observe that for any $\bm{x},\bm{y}\in\mathbb{R}^{d}$, we have

$\displaystyle\frac{\langle\nabla f(\bm{x}+\bm{y})-\nabla f(\bm{x})-\nabla
f(\bm{y}),\bm{x}\rangle}{\|\bm{x}\|^{2}}$ $\displaystyle=\frac{\langle\nabla
f(\bm{x}+\bm{y})-\nabla
f(\bm{x}),\bm{x}\rangle}{\|\bm{x}\|^{2}}-\frac{\langle\nabla
f(\bm{y}),\bm{x}\rangle}{\|\bm{x}\|^{2}}$ $\displaystyle\leq\frac{\|\nabla
f(\bm{x}+\bm{y})-\nabla f(\bm{x})\|\|\bm{x}\|}{\|\bm{x}\|^{2}}+\frac{\|\nabla
f(\bm{y})\|\|\bm{x}\|}{\|\bm{x}\|^{2}}$
$\displaystyle=\frac{M\|\bm{y}\|+\|\nabla f(\bm{y})\|}{\|\bm{x}\|}.$

Hence, we have

$\displaystyle\liminf_{\|\bm{x}\|\to\infty}\frac{\langle\nabla
f(\bm{x})-\nabla f(\bm{y}),\bm{x}-\bm{y}\rangle}{\|\bm{x}-\bm{y}\|^{2}}$
$\displaystyle=\liminf_{\|\bm{x}\|\to\infty}\frac{\langle\nabla
f(\bm{x}),\bm{x}\rangle}{\|\bm{x}\|^{2}}$
$\displaystyle\quad+\liminf_{\|\bm{x}\|\to\infty}\frac{\langle\nabla
f(\bm{x}+\bm{y})-\nabla f(\bm{x})-\nabla
f(\bm{y}),\bm{x}\rangle}{\|\bm{x}\|^{2}}$
$\displaystyle=\liminf_{\|\bm{x}\|\to\infty}\frac{\langle\nabla
f(\bm{x}),\bm{x}\rangle}{\|\bm{x}\|^{2}}>0.$

The result is implied by the definition of limit infimum.

To state the result, we recall the definition of $V$-uniform ergodicity, as
seen in Meyn and Tweedie [2012].

######  Definition 1.

A $\nu$-ergodic Markov chain $\\{\bm{X}_{n}\\}_{n=0}^{\infty}$ with Markov
transition operator $\mathcal{P}$ on $\mathbb{R}^{d}$ is $V$-uniformly ergodic
for a measurable function $V:\mathbb{R}^{d}\to[1,\infty)$ if

$\sup_{\bm{x}\in\mathbb{R}^{d}}\sup_{|\phi|\leq
V}\frac{|\mathcal{P}^{n}\phi(\bm{x})-\pi(\phi)|}{V(\bm{x})}\to
0,\qquad\mbox{as }n\to\infty.$

By [Meyn and Tweedie, 2012, Theorem 16.0.1], any $V$-uniformly ergodic Markov
chain is also geometrically ergodic.

######  Proposition 1.

For any $s>0$ and $f$ satisfying Assumptions 1 and 2, let $V_{s}(\bm{x})$
denote the Lyapunov drift function

$V_{s}(\bm{x})=\exp\left(s\|\bm{x}-\bm{x}^{\ast}+\tfrac{1}{2}h\theta\nabla
f(\bm{x})\|\right),$ (21)

where $\bm{x}^{\ast}$ is a critical point of $f$. Supposing that (5) holds,
the $\theta$-method scheme with transition kernel (6) is $V_{s}$-uniformly
ergodic provided $\theta\geq 1/2$, or $\theta<1/2$ and

$h<\frac{4m}{M^{2}(1-2\theta)}.$

###### Proof.

It is immediately apparent from the positivity of (6) due to (5) that the
iterates of the $\theta$-method scheme are aperiodic and irreducible with
respect to Lebesgue measure. Furthermore, it follows from [Meyn and Tweedie,
2012, Proposition 6.2.8] that all compact sets are small. Therefore, by [Meyn
and Tweedie, 2012, Theorem 15.0.1] and [Meyn and Tweedie, 2012, Lemma 15.2.8],
it suffices to show that

$\limsup_{\|\bm{x}\|\to\infty}\frac{\mathcal{P}V_{s}(\bm{x})}{V_{s}(\bm{x})}=0,$

where $\mathcal{P}$ is the Markov transition operator of the $\theta$-method
scheme. Indeed, by the definition of $\limsup$, for a given $0<\lambda<1$,
there exists a $K>0$ such that

$\sup_{\|\bm{x}\|\geq
K}\frac{\mathcal{P}V_{s}(\bm{x})}{V_{s}(\bm{x})}\leq\lambda,$

and so $\mathcal{P}V_{s}(\bm{x})\leq\lambda V_{s}(\bm{x})+\sup_{\|\bm{x}\|\leq
K}\mathcal{P}V_{s}(\bm{x})$ for any $\bm{x}\in\mathbb{R}^{d}$. Letting
$\bm{X}_{1}$ denote the first step of the $\theta$-method scheme starting from
$\bm{X}_{0}=\bm{x}$, (4) implies

$\bm{X}_{1}+\tfrac{1}{2}h\theta\nabla
f(\bm{X}_{1})\sim\mathcal{N}\left(\bm{x}-\frac{h(1-\theta)}{2}\nabla
f(\bm{x}),\ h\mathbf{I}\right).$

Thus, by letting $\bm{Z}\sim\mathcal{N}(\bm{0},\mathbf{I})$, we obtain
$\mathcal{P}V_{s}(\bm{x})=\mathbb{E}\exp(sg(\bm{Z}))$, where

$g(\bm{z})=\|\bm{x}-\bm{x}^{\ast}-\tfrac{1}{2}h(1-\theta)\nabla
f(\bm{x})+\sqrt{h}\bm{z}\|.$

By the reverse triangle inequality,
$|g(\bm{z}_{1})-g(\bm{z}_{2})|\leq\sqrt{h}\|\bm{z}_{1}-\bm{z}_{2}\|$ for any
$\bm{z}_{1},\bm{z}_{2}\in\mathbb{R}^{d}$, and hence, $g$ is
$\sqrt{h}$-Lipschitz in $\bm{z}$. Consequently, we can apply the Gaussian
concentration inequality [Boucheron et al., 2013, Theorem 5.5] to reveal

$\mathbb{E}e^{sg(\bm{Z})}\leq\exp\left(s\mathbb{E}g(\bm{Z})+\frac{hs^{2}}{2}\right).$

Since by Jensen’s inequality,

$\mathbb{E}g(\bm{Z})\leq\sqrt{hd}+\|(\bm{x}-\bm{x}^{\ast})-\tfrac{1}{2}h(1-\theta)\nabla
f(\bm{x})\|.$

It follows that

$\frac{\mathcal{P}V_{s}(\bm{x})}{V_{s}(\bm{x})}\leq\exp(s\sqrt{hd}+\tfrac{1}{2}hs^{2}+s[T_{1}(\bm{x})-T_{2}(\bm{x})]),$

where

$T_{1}(\bm{x})=\|\bm{x}-\bm{x}^{\ast}-\tfrac{1}{2}h(1-\theta)\nabla
f(\bm{x})\|\quad\mbox{and}\quad
T_{2}(\bm{x})=\|\bm{x}-\bm{x}^{\ast}+\tfrac{1}{2}h\theta\nabla f(\bm{x})\|.$

Therefore, if we can show that $T_{1}(\bm{x})-T_{2}(\bm{x})\to-\infty$ as
$\|\bm{x}\|\to\infty$, then the result follows. Since
$T_{1}(\bm{x})-T_{2}(\bm{x})=(T_{1}(\bm{x})^{2}-T_{2}(\bm{x})^{2})/(T_{1}(\bm{x})+T_{2}(\bm{x}))$,
we may focus on the difference of the squares:

$\displaystyle T_{1}(\bm{x})^{2}-T_{2}(\bm{x})^{2}$
$\displaystyle=\|(\bm{x}-\bm{x}^{\ast})-\tfrac{1}{2}h(1-\theta)\nabla
f(\bm{x})\|^{2}-\|(\bm{x}-\bm{x}^{\ast})+\tfrac{1}{2}h\theta\nabla
f(\bm{x})\|^{2}$
$\displaystyle=\|\bm{x}-\bm{x}^{\ast}\|^{2}-h(1-\theta)\langle\nabla
f(\bm{x}),\bm{x}-\bm{x}^{\ast}\rangle+\tfrac{1}{4}h^{2}(1-\theta)^{2}\|\nabla
f(\bm{x})\|^{2}$
$\displaystyle\qquad-\|\bm{x}-\bm{x}^{\ast}\|^{2}-h\theta\langle\nabla
f(\bm{x}),\bm{x}-\bm{x}^{\ast}\rangle-\tfrac{1}{4}h^{2}\theta^{2}\|\nabla
f(\bm{x})\|^{2}$ $\displaystyle=-h\langle\nabla
f(\bm{x}),\bm{x}-\bm{x}^{\ast}\rangle+\tfrac{1}{4}h^{2}(1-2\theta)\|\nabla
f(\bm{x})\|^{2}$ $\displaystyle\leq-
hm\|\bm{x}-\bm{x}^{\ast}\|^{2}+c(\bm{x}^{\ast})+\tfrac{1}{4}h^{2}\max\\{0,1-2\theta\\}M^{2}\|\bm{x}-\bm{x}^{\ast}\|^{2},$
(22)

where the last inequality follows from (7). Provided that
$\frac{1}{2}h(1-2\theta)M^{2}<m$ or $\theta\geq 1/2$, (22) will be negative
for sufficiently large $\bm{x}$. Since also

$T_{1}(\bm{x})+T_{2}(\bm{x})\leq
2\|\bm{x}-\bm{x}^{\ast}\|+\tfrac{1}{2}h\|\nabla
f(\bm{x})\|\leq(2+\tfrac{1}{2}hM)\|\bm{x}-\bm{x}^{\ast}\|,$

for any $\epsilon>0$ and sufficiently large $\bm{x}$,

$T_{1}(\bm{x})-T_{2}(\bm{x})\leq\frac{-hm+\frac{1}{4}h^{2}\max\\{0,1-2\theta\\}M^{2}}{2+\frac{1}{2}hM}\|\bm{x}-\bm{x}^{\ast}\|+\epsilon,$

which implies the difference $T_{1}(\bm{x})-T_{2}(\bm{x})\to-\infty$ as
$\|\bm{x}\|\to\infty$, as required.

###  A.2 Proof of Theorem 2 ($W_{2}$ bounds)

Next, using techniques analogous to those of Dalalyan Dalalyan and Karagulyan
[2017], we prove Theorem 2. For the sake of brevity, we let $a\wedge b$ denote
the minimum of any two quantities $a$ and $b$. The following estimate is
fundamental to the argument.

######  Lemma 2.

Let $\bm{L}_{t}$ be the solution to the (overdamped) Langevin equation

$\mathrm{d}\bm{L}_{t}=-\tfrac{1}{2}\nabla
f\left(\bm{L}_{t}\right)\mathrm{d}t+\mathrm{d}\bm{W}_{t}$

for $f\in\mathcal{C}^{1}(\mathbb{R}^{d})$ such that $\nabla f$ is
$M$-Lipschitz continuous. Then for any $h>0$, if $\bm{L}_{0}\sim\pi$,

$\left\lVert\int_{0}^{h}\nabla f(\bm{L}_{t})-\nabla
f(\bm{L}_{0})\mathrm{d}t\right\rVert_{L^{2}}\leq\frac{1}{2}h[M\sqrt{hd}(2+\sqrt{hM})\wedge
4\sqrt{Md}].$

###### Proof.

Since $\bm{L}_{t}$ is stationary, for any $t\geq 0$, $\|\nabla
f(\bm{L}_{t})\|_{L^{2}}\leq\sqrt{Md}$ by [Dalalyan, 2017a, Lemma 2].
Therefore, $\left\lVert\int_{0}^{h}\nabla f(\bm{L}_{t})-\nabla
f(\bm{L}_{0})\mathrm{d}t\right\rVert_{L^{2}}\leq 2h\sqrt{Md}$. Furthermore,
following the same procedure as in [Dalalyan and Karagulyan, 2017, Lemma 4]

$\left\lVert\int_{0}^{h}\nabla f(\bm{L}_{t})-\nabla
f(\bm{L}_{0})\mathrm{d}t\right\rVert_{L^{2}}\leq\frac{1}{4}h^{2}M^{3/2}d^{1/2}+\frac{2}{3}h^{3/2}Md^{1/2}.\@qedbox{}$

With Lemma 2 in hand, we may proceed with the proof of the main result.

###### Theorem 2.

Letting $\bm{W}_{t}$ denote a $d$-dimensional standard Brownian motion
independent of $\bm{X}_{k}$ and $\bm{L}_{0}\sim\pi$, define the stochastic
process $\bm{L}$ by

$\bm{L}_{t}=\bm{L}_{0}-\frac{1}{2}\int_{0}^{t}\nabla
f(\bm{L}_{s})\mathrm{d}s+\bm{W}_{t},\qquad\qquad t\geq 0.$

Evidently, $\bm{L}_{t}$ is a realization of (1) and so is a reversible Markov
process with $\bm{L}_{t}\sim\pi$ for every $t\geq 0$. Now, couple the inexact
$\theta$-method scheme $\bm{X}_{k}$ satisfying

$\bm{X}_{k+1}=\bm{X}_{k}-\frac{h}{2}[\theta\nabla
f(\bm{X}_{k+1})+(1-\theta)\nabla
f(\bm{X}_{k})]+\sqrt{h}\bm{Z}_{k}+\bm{E}_{k},$

for an appropriate error term $\bm{E}_{k}$, to $\bm{L}_{t}$, by letting
$\bm{Z}_{k}=h^{-1/2}[\bm{W}_{(k+1)h}-\bm{W}_{kh}]$ for each $k\geq 1$, and
choosing $\bm{L}_{0}$ such that
$W_{2}(\pi_{0},\pi)=\|\bm{X}_{0}-\bm{L}_{0}\|_{L^{2}}$. Observing that, for
any $k\geq 1$,

$\displaystyle\bm{E}_{k}$ $\displaystyle=\frac{h}{2}\theta\nabla
f(\bm{X}_{k+1})+\bm{X}_{k+1}-\bm{X}_{k}+\frac{h}{2}(1-\theta)\nabla
f(\bm{X}_{k})-\sqrt{h}\bm{Z}_{k}$ $\displaystyle=\frac{h}{2}\nabla
F(\bm{X}_{k+1};\bm{X}_{k},\bm{Z}_{k}),$

by construction, $\|\bm{E}_{k}\|_{L^{2}}\leq\frac{1}{2}h\epsilon$. For each
$k$, let $\bm{D}_{k}=\bm{L}_{kh}-\bm{X}_{k}$, observing that
$W_{2}(\pi_{0},\pi)=\|\bm{D}_{0}\|_{2}$ and
$W_{2}(\pi_{k},\pi)\leq\|\bm{D}_{k}\|_{2}$. Choosing some $k\geq 1$, for the
sake of brevity, we denote $\bm{L}_{t}^{(k)}=\bm{L}_{kh+t}$, which now
satisfies

$\displaystyle\bm{L}_{h}^{(k)}=\bm{L}_{kh+h}$
$\displaystyle=\bm{L}_{kh}-\frac{1}{2}\int_{0}^{h}\nabla
f(\bm{L}_{kh+s})\mathrm{d}s+\bm{W}_{kh+h}-\bm{W}_{kh}$
$\displaystyle=\bm{L}_{0}^{(k)}-\frac{1}{2}\int_{0}^{h}\nabla
f(\bm{L}_{s}^{(k)})\mathrm{d}s+\sqrt{h}\bm{Z}_{k}.$

Altogether, we have

$\bm{D}_{k+1}=\bm{D}_{k}-\tfrac{1}{2}h[(1-\theta)\bm{U}_{k}+\theta\tilde{\bm{U}}_{k}]-[(1-\theta)\bm{V}_{k}+\theta\tilde{\bm{V}}_{k}]+\bm{E}_{k},$

where

$\displaystyle\bm{U}_{k}$ $\displaystyle=\nabla
f(\bm{X}_{k}+\bm{D}_{k})-\nabla f(\bm{X}_{k})$ $\displaystyle\bm{V}_{k}$
$\displaystyle=\frac{1}{2}\int_{0}^{h}\nabla f(\bm{L}_{s}^{(k)})-\nabla
f(\bm{L}_{0}^{(k)})\mathrm{d}s$ $\displaystyle\bm{\tilde{}}{\bm{U}}_{k}$
$\displaystyle=\nabla f(\bm{X}_{k+1}+\bm{D}_{k+1})-\nabla f(\bm{X}_{k+1})$
$\displaystyle\bm{\tilde{}}{\bm{V}}_{k}$
$\displaystyle=\frac{1}{2}\int_{0}^{h}\nabla f(\bm{L}_{h-s}^{(k)})-\nabla
f(\bm{L}_{h}^{(k)})\mathrm{d}s.$

An application of the fundamental theorem of calculus implies
$\bm{U}_{k}=\mathbf{F}_{k}\bm{D}_{k}$ and
$\bm{\tilde{}}{\bm{U}}_{k}=\mathbf{F}_{k+1}\bm{D}_{k+1}$, where

$\mathbf{F}_{k}=\int_{0}^{1}\nabla^{2}f(\bm{X}_{k}+t\bm{D}_{k})\mathrm{d}t.$

Altogether, $\bm{D}_{k+1}=\mathbf{S}_{k}\bm{D}_{k}+\bm{T}_{k}$ where
$\bm{T}_{k}=-(\mathbf{I}+\frac{h\theta}{2}\mathbf{F}_{k+1})^{-1}[(1-\theta)\bm{V}_{k}+\theta\tilde{\bm{V}}_{k}-\bm{E}_{k}]$
and

$\mathbf{S}_{k}=\left(\mathbf{I}+\frac{h\theta}{2}\mathbf{F}_{k+1}\right)^{-1}\left(\mathbf{I}-\frac{h(1-\theta)}{2}\mathbf{F}_{k}\right).$

It can be verified using induction that the solution to this first-order non-
homogeneous recurrence relation is given by

$\bm{D}_{k}=\mathbf{S}_{k-1}\cdots\mathbf{S}_{0}\bm{D}_{0}+\sum_{l=0}^{k-1}\mathbf{S}_{k-1}\cdots\mathbf{S}_{l+1}\bm{T}_{l}.$

Now, observe that by denoting
$G(\mathbf{X})=(\mathbf{I}-\frac{h(1-\theta)}{2}\mathbf{X})(\mathbf{I}+\frac{h\theta}{2}\mathbf{X})^{-1}$,
for any $l<k$,

$\mathbf{S}_{k-1}\cdots\mathbf{S}_{l}=\left(\mathbf{I}+\frac{h\theta}{2}\mathbf{F}_{k}\right)^{-1}G(\mathbf{F}_{k-1})\cdots
G(\mathbf{F}_{l+1})\left(\mathbf{I}-\frac{h(1-\theta)}{2}\mathbf{F}_{l}\right).$
(23)

Since the eigenvalues of $\nabla^{2}f$ are bounded above by $M$ and below by
$m$, so too are the eigenvalues of $\mathbf{F}_{k}$ for each $k$. Therefore,

$\displaystyle\|G(\mathbf{F}_{k})\|_{2}$
$\displaystyle=\max_{z\in[m,M]}\left|\frac{1-\frac{1}{2}h(1-\theta)z}{1+\frac{1}{2}h\theta
z}\right|$
$\displaystyle=\max\left\\{\frac{1-\frac{1}{2}h(1-\theta)m}{1+\frac{1}{2}h\theta
m},\frac{\frac{1}{2}h(1-\theta)M-1}{\frac{1}{2}h\theta
M+1}\right\\}\eqqcolon\rho.$ (24)

The transition between these regimes occurs at the point $h^{\ast}$ which is
the solution to

$\frac{1-\frac{1}{2}h(1-\theta)m}{1+\frac{1}{2}h\theta
m}=\frac{\frac{1}{2}h(1-\theta)M-1}{\frac{1}{2}h\theta M+1}$

over $h>0$. Equivalently, it is the solution to

$\tfrac{1}{2}h\theta
M+1-\tfrac{1}{4}h^{2}\theta(1-\theta)mM-\tfrac{1}{2}h(1-\theta)m\\\
=\tfrac{1}{2}h(1-\theta)M+\tfrac{1}{4}h^{2}\theta(1-\theta)mM-1-\tfrac{1}{2}h\theta
m,$

and therefore to the quadratic equation

$\tfrac{1}{2}h(1-2\theta)(m+M)+\tfrac{1}{2}h^{2}\theta(1-\theta)mM-2=0.$

It may be readily verified that $h^{\ast}$ as defined in (12) is the only
positive solution. Furthermore, $\rho<1$ provided that $\theta\geq 1/2$ or
$h<4/[M(1-2\theta)]$. Also, for any $j,k\geq 1$,

$\|(\mathbf{I}+\tfrac{1}{2}h\theta\mathbf{F}_{j})^{-1}\|_{2}\|\mathbf{I}-\tfrac{1}{2}h(1-\theta)\mathbf{F}_{k}\|_{2}\\\
\leq\frac{\max\\{1-\tfrac{1}{2}hm(1-\theta),\tfrac{1}{2}hM(1-\theta)-1\\}}{1+\frac{1}{2}hm\theta}\leq\kappa_{h}\rho,$
(25)

which further implies $\|\mathbf{S}_{j}\|_{2}\leq\kappa_{h}\rho$ for any $j$.
Now combining (23), (24), and (25), for $k>l$,
$\|\prod_{j=l}^{k-1}\mathbf{S}_{j}\|_{2}\leq\kappa_{h}\rho^{k-l}$, and hence,
altogether,

$\|\bm{D}_{k}\|_{2}\leq\kappa_{h}\rho^{k}\|\bm{D}_{0}\|_{2}+\sum_{l=0}^{k-1}\kappa_{h}\rho^{k-l}\|\bm{T}_{l}\|_{2}.$

Since $\bm{L}_{t}$ is reversible and stationary,
$\|\bm{V}_{k}^{\ast}\|_{2}=\|\bm{V}_{k}\|_{2}$, and by Lemma 2,
$\|\bm{V}_{k}\|_{2}\leq h[M\sqrt{hd}(2+\sqrt{hM})\wedge 4\sqrt{Md}]$.
Therefore

$\|\bm{T}_{k}\|_{2}\leq\frac{h[\frac{1}{2}\epsilon+M\sqrt{hd}(2+\sqrt{hM})\wedge
4\sqrt{Md}]}{1+\frac{1}{2}hm\theta}.$

Since

$\displaystyle 1-\frac{1-\frac{1}{2}h(1-\theta)m}{1+\frac{1}{2}h\theta m}$
$\displaystyle=\frac{1+\frac{1}{2}h\theta
m-1+\frac{1}{2}h(1-\theta)m}{1+\frac{1}{2}h\theta m}=\frac{\frac{1}{2}h\theta
m}{1+\frac{1}{2}h\theta m},$ $\displaystyle
1-\frac{\frac{1}{2}h(1-\theta)M-1}{\frac{1}{2}h\theta M+1}$
$\displaystyle=\frac{1+\frac{1}{2}h\theta
M+1-\frac{1}{2}h(1-\theta)M}{\frac{1}{2}h\theta
M+1}=\frac{2+\frac{1}{2}h(2\theta-1)M}{\frac{1}{2}h\theta M+1},$

applying the closed-form expression for the geometric series,

$\sum_{l=0}^{k-1}\rho^{k-l}\leq\frac{1}{1-\rho}=\max\left\\{\frac{1+\frac{1}{2}h\theta
m}{\frac{1}{2}h\theta m},\frac{\frac{1}{2}h\theta
M+1}{2+\frac{1}{2}h(2\theta-1)M}\right\\},$

and the result follows.

###  A.3 Proof of Theorem 3 (Central Limit Theorem)

###### Proof.

The proof makes use of Laplace’s method. From 6 and the change of variables
theorem, the Markov kernel $\tilde{p}_{h}(\bm{y}|\bm{x})$ for the transition
$\bm{X}_{k}\mapsto\sqrt{h}(\bm{X}_{k+1}-\bm{X}_{k}^{\theta})$ is given by

$\tilde{p}_{h}(\bm{y}|\bm{x})=(2\pi)^{-d/2}\det\left(\frac{1}{h}\mathbf{I}+\frac{\theta}{2}\nabla^{2}f(\bm{x}^{\theta}+h^{-1/2}\bm{y})\right)\times\\\
\exp\left(-\frac{1}{2h}\left\lVert\bm{x}^{\theta}-\bm{x}+h^{-1/2}\bm{y}+\frac{h\theta}{2}\nabla
f(\bm{x}^{\theta}+h^{-1/2}\bm{y})+\frac{h(1-\theta)}{2}\nabla
f(\bm{x})\right\rVert^{2}\right).$ (26)

Letting $q(\bm{y}|\bm{x})=\phi(\bm{y};\bm{0},\Sigma(\bm{x}))$ where
$\Sigma(\bm{x})=(4/\theta^{2})\nabla^{2}f(\bm{x}^{\theta})^{-2}$, it suffices
to show that $\tilde{p}_{h}(\bm{y}|\bm{x})\to q(\bm{y}|\bm{x})$ as
$h\to\infty$, pointwise in $\bm{y}$. Denoting
$H_{h}(\bm{y})=h^{-1}\mathbf{I}+\frac{1}{2}\theta\int_{0}^{1}\nabla^{2}f(\bm{x}^{\theta}+th^{-1/2}\bm{y})dt$,
since $\theta\nabla f(\bm{x}^{\theta})=-(1-\theta)\nabla f(\bm{x})$, the
exponent of (26) becomes

$-\frac{1}{2h}\|\bm{x}^{\theta}-\bm{x}+\sqrt{h}H_{h}(\bm{y})\bm{y}\|^{2}=-\frac{\|\bm{x}^{\theta}-\bm{x}\|^{2}}{2h}-\frac{(\bm{x}^{\theta}-\bm{x})^{\top}H_{h}(\bm{y})\bm{y}}{\sqrt{h}}-\frac{1}{2}\bm{y}^{\top}H_{h}(\bm{y})^{2}\bm{y}.$

Now, since $H_{h}(\bm{y})^{2}\to\Sigma(\bm{x})^{-1}$ and the determinant term
converges to $\det(\Sigma^{-1/2})$, the result follows.