# The shifted ODE method for underdamped Langevin MCMC

James Foster, Terry Lyons and Harald Oberhauser

###### Abstract

In this paper, we consider the underdamped Langevin diffusion (ULD) and
propose a numerical approximation using its associated ordinary differential
equation (ODE). When used as a Markov Chain Monte Carlo (MCMC) algorithm, we
show that the ODE approximation achieves a $2\hskip 0.7113pt$-Wasserstein
error of $\varepsilon$ in $\mathcal{O}\big{(}\sqrt[3]{d\hskip
0.7113pt}/\varepsilon^{\frac{2}{3}}\big{)}$ steps under the standard
smoothness and strong convexity assumptions on the target distribution. This
matches the complexity of the randomized midpoint method proposed by Shen and
Lee [NeurIPS 2019] which was shown to be order optimal by Cao, Lu and Wang.
However, the main feature of the proposed numerical method is that it can
utilize additional smoothness of the target log-density $f$. More concretely,
we show that the ODE approximation achieves a $2\hskip 0.7113pt$-Wasserstein
error of $\varepsilon$ in
$\mathcal{O}\big{(}d^{\frac{2}{5}}/\varepsilon^{\frac{2}{5}}\big{)}$ and
$\mathcal{O}\big{(}\sqrt{d\hskip 0.7113pt}/\varepsilon^{\frac{1}{3}}\big{)}$
steps when Lipschitz continuity is assumed for the Hessian and third
derivative of $f$. By discretizing this ODE using a third order Runge-Kutta
method, we can obtain a practical MCMC method that uses just two additional
gradient evaluations per step. In our experiment, where the target comes from
a logistic regression, this method shows faster convergence compared to other
unadjusted Langevin MCMC algorithms.

##  1 Introduction

Traditionally used to model coarse-grained particle systems [34, 38, 14], the
underdamped Langevin diffusion (ULD) is given by the following stochastic
differential equation (SDE):

$\displaystyle dx_{t}$ $\displaystyle=v_{t}\,dt,$ (1.1) $\displaystyle dv_{t}$
$\displaystyle=-\gamma v_{t}\,dt-u\nabla f(x_{t})\,dt+\sqrt{2\gamma
u}\,dW_{t}\hskip 0.7113pt,$

where $x,v\in\mathbb{R}^{d}$ represent the position and momentum of a
particle, $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ denotes a scalar potential,
$\gamma>0$ is a friction coefficient and $W$ is a $d$-dimensional Brownian
motion.

Under mild assumptions, it is known that (1.1) admits a unique strong solution
which is ergodic with respect to the Gibbs measure with Hamiltonian
$H(x,v)=f(x)+\frac{1}{2u}\|v\|^{2}$ (we refer the reader to Proposition 6.1 in
[50] for further details). That is, the diffusion process has a stationary
measure $\pi$ on $\mathbb{R}^{2d}$ with density proportional to
$\exp(-H(x,v))$. In particular, the (marginal) stationary distribution for $x$
is proportional to $\exp(-f(x))$ and thus solving (1.1) is a method of
sampling from unnormalized log-concave densities. When the SDE solution cannot
be obtained exactly, it must be appropriately discretized and so numerical
methods are required to apply ULD to log-concave sampling problems. In this
paper, we develop an ordinary differential equation (ODE) based method for
ULD.

Langevin-type equations are widely used in statistical physics as they are
natural extensions of Hamiltonian systems that can model randomness [14, 46].
More recently, they have been applied to sampling and optimization problems in
machine learning [56, 44, 43]. As a result, the study of (underdamped)
Langevin dynamics is an active area of research with an extensive literature.
Below we give a brief overview for some of this related work:

Numerical methods for underdamped Langevin dynamics

A variety of approximations have been proposed for ULD. Some noteworthy
examples are the unadjusted and Metropolis-adjusted OBABO schemes [4, 47, 54],
the randomized midpoint method [52, 31] (which was shown to be order optimal
in [8]), the left-point111In [12], the proposed approximation is referred to
as the “discrete underdamped Langevin diffusion”. We shall refer to this
approach as the “left-point method” (simply to distinguish it from other
methods). method [12, 16, 11], Strang splitting [3, 35] and quasi-symplectic
Runge-Kutta methods [45]. We shall compare our approach to several of these
methods in our numerical experiment.

Probabilistic couplings for underdamped Langevin dynamics

The analysis of ULD and its discretizations often makes use of probabilistic
couplings. In [22], the authors used reflection and synchronous couplings to
establish a (Wasserstein) contraction property for the diffusion. Similar
contraction theorems were later obtained in [12, 16] where only synchronous
couplings were used. The exponential contractivity of Langevin dynamics allows
one to obtain non-asymptotic estimates between the diffusion and approximation
processes directly from local mean and mean-square estimates [39]. Contraction
results have also been established for the OBABO discretization of (1.1) using
synchronous and Gaussian couplings [47]. However, unlike in [39], the author
of [47] obtains non-asymptotic estimates without comparing to the SDE solution
(instead the OBABO chain is compared to another Markov chain with the desired
stationary measure). In [19], the authors propose the use of couplings to
empirically estimate the convergence of a numerical solution to the stationary
measure of an ergodic SDE. For example, they use a combination of reflection,
synchronous and maximal couplings to estimate the convergence of an Euler-
Maruyama discretization of ULD with nonconvex potential
$f(x)=\big{(}\|x\|_{2}^{2}-1\big{)}^{2}$.

Applications of Langevin dynamics to MCMC algorithms and optimization

After the breakthrough papers of Neal [49] along with Girolami and Calderhead
[30], there has significant interest in Hamiltonian Monte Carlo (HMC) methods
for sampling and integral estimation. Since these MCMC methods are based on
Hamiltonian systems, they are intimately connected to Langevin dynamics [21,
43, 10, 54, 2]. More recently, the Metropolis-adjusted OBABO sampler (an MCMC
algorithm which converges to ULD [47]) was shown to empirically outperform the
standard HMC method [54]. Alongside MCMC [18, 54, 55], underdamped Langevin
diffusions have also been considered for optimization due to their natural
parallels with existing momentum-based algorithms [28, 29, 42, 36, 53].

### 1.1. Contributions

Our key contribution is an ODE-based numerical method for ULD, which we refer
to as the “shifted ODE method”. To the best of our knowledge, this is the
first derivative-free approximation of ULD that can achieve third order
convergence (a third order numerical method is given in [45] but it requires
further derivatives of $\nabla f$).

######  Definition 1.1 (Shifted ODE method).

Let $\\{t_{n}\\}_{n\hskip 0.7113pt\geq\hskip 0.7113pt0}$ be a sequence of
times with $t_{0}=0$, $t_{n+1}>t_{n}$ and step sizes
$h_{n}=t_{n+1}-t_{n}\hskip 0.7113pt$, for $n\geq 0$. We construct a numerical
solution $\\{(\widetilde{x}_{n}\hskip 0.7113pt,\widetilde{v}_{n})\\}_{n\hskip
0.7113pt\geq\hskip 0.7113pt0}$ by setting $(\widetilde{x}_{0}\hskip
0.7113pt,\widetilde{v}_{0}):=(x_{0}\hskip 0.7113pt,v_{0})$ and for each $n\geq
0$, defining $(\widetilde{x}_{n+1}\hskip 0.7113pt,\widetilde{v}_{n+1})$ as

$\displaystyle\Bigg{(}\begin{matrix}\,\widetilde{x}_{n+1}\\\\[-3.0pt]
\,\widetilde{v}_{n+1}\end{matrix}\Bigg{)}:=\Bigg{(}\begin{matrix}\,\overline{x}_{1}^{\hskip
0.7113ptn}\\\\[-3.0pt] \,\overline{v}_{1}^{\hskip
0.7113ptn}\end{matrix}\Bigg{)}-\big{(}H_{n}-6K_{n}\big{)}\Bigg{(}\begin{matrix}0\\\\[-3.0pt]
\sqrt{2\gamma u}\end{matrix}\,\Bigg{)},$

where $\\{(\overline{x}_{t}^{\hskip 0.7113ptn},\overline{v}_{t}^{\hskip
0.7113ptn})\\}_{t\hskip 0.7113pt\in\hskip 0.7113pt[0,1]}$ solves the following
(rescaled) Langevin-type ODE,

$\displaystyle\frac{d}{dt}\Bigg{(}\begin{matrix}\,\overline{x}^{\hskip
0.7113ptn}\\\\[-4.0pt] \,\overline{v}^{\hskip
0.7113ptn}\end{matrix}\Bigg{)}=\Bigg{(}\begin{matrix}\,\overline{v}^{\hskip
0.7113ptn}\\\\[-4.0pt] \,-\gamma\hskip 0.7113pt\overline{v}^{\hskip
0.7113ptn}-u\nabla f\big{(}\,\overline{x}^{\hskip
0.7113ptn}\big{)}\end{matrix}\,\Bigg{)}\hskip
0.7113pth_{n}+\big{(}W_{n}-12K_{n}\big{)}\Bigg{(}\begin{matrix}0\\\\[-4.0pt]
\sqrt{2\gamma u}\end{matrix}\,\Bigg{)},$ (1.2)

with initial condition

$\displaystyle\Bigg{(}\begin{matrix}\,\overline{x}_{0}^{\hskip
0.7113ptn}\\\\[-3.0pt] \,\overline{v}_{0}^{\hskip
0.7113ptn}\end{matrix}\Bigg{)}:=\Bigg{(}\begin{matrix}\,\widetilde{x}_{n}\\\\[-3.0pt]
\,\widetilde{v}_{n}\end{matrix}\Bigg{)}+\big{(}H_{n}+6K_{n}\big{)}\Bigg{(}\begin{matrix}0\\\\[-3.0pt]
\sqrt{2\gamma u}\end{matrix}\,\Bigg{)},$

and the random Gaussian vectors $\\{W_{n}\\}_{n\hskip 0.7113pt\geq\hskip
0.7113pt0}\hskip 0.7113pt,\,\\{H_{n}\\}_{n\hskip 0.7113pt\geq\hskip
0.7113pt0}$ and $\\{K_{n}\\}_{n\hskip 0.7113pt\geq\hskip 0.7113pt0}$ are
independent with $W_{n}\sim\mathcal{N}\big{(}0,h_{n}\hskip
0.7113ptI_{d}\big{)}\hskip
0.7113pt,\,H_{n}\sim\mathcal{N}\big{(}0,\frac{1}{12}h_{n}\hskip
0.7113ptI_{d}\big{)}$ and $K_{n}\sim\mathcal{N}\big{(}0,\frac{1}{720}\hskip
0.7113pth_{n}\hskip 0.7113ptI_{d}\big{)}$ for all $n\geq 0$. Moreover, we
define $\\{\hskip 0.7113ptW_{n}\hskip 0.7113pt,\hskip 0.7113ptH_{n}\hskip
0.7113pt,\hskip 0.7113ptK_{n}\hskip 0.7113pt\\}$ from the same Brownian motion
$W$ (see Theorem 3.5 for details).

$\displaystyle W_{n}$ $\displaystyle:=W_{t_{n+1}}-W_{t_{n}}\hskip 0.7113pt,$
$\displaystyle H_{n}$
$\displaystyle:=\frac{1}{h_{n}}\int_{t_{n}}^{t_{n+1}}\Big{(}\big{(}W_{t}-W_{t_{n}}\big{)}-\frac{t-t_{n}}{h_{n}}\hskip
0.7113ptW_{n}\Big{)}\,dt,$ $\displaystyle K_{n}$
$\displaystyle:=\frac{1}{h_{n}^{2}}\int_{t_{n}}^{t_{n+1}}\bigg{(}\frac{1}{2}h_{n}-(t-t_{n})\bigg{)}\Big{(}\big{(}W_{t}-W_{t_{n+1}}\big{)}-\frac{t-t_{n}}{h_{n}}\hskip
0.7113ptW_{n}\Big{)}\,dt.$

We discuss the intuition behind this method in Section 3 and perform an error
analysis in the appendices. In Table 1.1, we compare this ODE approach to some
related methods:

[h] Summary of complexities for ULD methods (with respect to $2$ -Wasserstein
error). Numerical method Smoothness assumptions Number of steps $n$ to achieve
on the strongly convex $f$ an error of $W_{2}\big{(}\hskip
0.7113pt\widetilde{x}_{n}\hskip 0.7113pt,e^{-f}\hskip
0.7113pt\big{)}\leq\varepsilon$ Shifted ODE method Lipschitz Gradient
$\mathcal{O}\big{(}\hskip 0.7113pt\sqrt[3]{d\hskip
0.7113pt}/\varepsilon^{\frac{2}{3}}\hskip 0.7113pt\big{)}$* (this paper) \+
Lipschitz Hessian $\mathcal{O}\big{(}\hskip
0.7113ptd^{\frac{2}{5}}/\varepsilon^{\frac{2}{5}}\hskip 0.7113pt\big{)}$* \+
Lipschitz third derivative $\mathcal{O}\big{(}\hskip 0.7113pt\sqrt{d\hskip
0.7113pt}/\varepsilon^{\frac{1}{3}}\hskip 0.7113pt\big{)}$ Strang splitting
[3] Lipschitz Gradient $\mathcal{O}\big{(}\sqrt{d\hskip
0.7113pt}/\varepsilon\hskip 0.7113pt\big{)}$$\dagger$ \+ Lipschitz Hessian
$\mathcal{O}\big{(}\sqrt{d\hskip 0.7113pt}/\sqrt{\varepsilon}\hskip
0.7113pt\big{)}$$\dagger$ OBABO scheme [47] Lipschitz Gradient
$\mathcal{O}\big{(}\sqrt{d\hskip 0.7113pt}/\varepsilon\hskip 0.7113pt\big{)}$
\+ Lipschitz Hessian $\mathcal{O}\big{(}\sqrt{d\hskip
0.7113pt}/\sqrt{\varepsilon}\hskip 0.7113pt\big{)}$ Randomized midpoint
Lipschitz Gradient $\mathcal{O}\big{(}\hskip 0.7113pt\sqrt[3]{d\hskip
0.7113pt}/\varepsilon^{\frac{2}{3}}\hskip 0.7113pt\big{)}$ method [52] Left-
point method [12] Lipschitz Gradient $\mathcal{O}\big{(}\hskip
0.7113pt\sqrt{d\hskip 0.7113pt}/\varepsilon\big{)}$

  * $\ast$

We expect these complexities to be reduced when the ODE is discretized using
standard high order ODE solvers. In this case, we expect the complexities to
be $\mathcal{O}\big{(}\sqrt{d\hskip 0.7113pt}/\varepsilon\hskip
0.7113pt\big{)}$ and $\mathcal{O}\big{(}\sqrt{d\hskip
0.7113pt}/\sqrt{\varepsilon}\hskip 0.7113pt\big{)}$ respectively. On the other
hand, if a randomized ODE solver is applied (such as a randomized midpoint
method), we conjecture that it is possible to achieve these complexities. This
may be a topic for future work.

  * $\dagger$

These orders are conjectured, but seem quite likely given the method’s
stochastic Taylor expansion.

In practice, the ODE (1.2) cannot be solved exactly, so must be discretized
using an appropriate ODE solver. For the resulting method to demonstrate third
order convergence, we expect (but have not proven) that it suffices to use a
third order numerical ODE solver. Therefore we consider two approaches which
are derived from the standard third order Runge-Kutta method [5] and the
fourth order splitting method of Forest and Ruth [23]. In order to distinguish
them, we refer to these discretizations of ULD as the SORT222Shifted ODE with
Runge-Kutta-Three. and SOFA333Shifted ODE with Fourth-order-splitting Applied.
methods. Due to the structure of the “shifted ODE”, the SORT and SOFA methods
only require two and three additional evaluations of $\nabla f$ per step
respectively. This means the SORT and SOFA methods have similar cost to the
methods in Table 1.1, which use one extra evaluation per step (except the
randomized midpoint which uses two). However by using more evaluations of
$\nabla f$ per step, the SORT and SOFA methods can exhibit faster convergence
when $\nabla f$ has more smoothness (see the numerical experiment).

######  Definition 1.2 (The SORT method).

Let $\\{t_{n}\\}_{n\hskip 0.7113pt\geq\hskip 0.7113pt0}$ be a sequence of
times with $t_{0}=0$, $t_{n+1}>t_{n}$ and step sizes
$h_{n}=t_{n+1}-t_{n}\hskip 0.7113pt$. We construct a numerical solution
$\\{(\hskip 0.7113pt\overrightarrow{x\,}_{n}\hskip
0.7113pt,\overrightarrow{v\,}_{n})\\}_{n\hskip 0.7113pt\geq\hskip 0.7113pt0}$
for (1.1) by setting $\hskip 0.7113pt(\hskip
0.7113pt\overrightarrow{x\,}_{0}\hskip
0.7113pt,\overrightarrow{v\,}_{0})\hskip 0.7113pt:=\hskip 0.7113pt(x_{0}\hskip
0.7113pt,v_{0})\hskip 0.7113pt$ and for $n\hskip 0.7113pt\geq\hskip
0.7113pt0$, defining $\hskip 0.7113pt(\hskip
0.7113pt\overrightarrow{x\,}_{n+1}\hskip
0.7113pt,\overrightarrow{v\,}_{n+1}\hskip 0.7113pt)\hskip 0.7113pt$ as
follows:

$\displaystyle\overrightarrow{v\,}_{n}^{\hskip 0.7113pt(1)}$
$\displaystyle:=\overrightarrow{v\,}_{n}+\sqrt{2\gamma u}\hskip
0.7113pt\big{(}H_{n}+6K_{n}\big{)},$
$\displaystyle\overrightarrow{x\,}_{n}^{\hskip 0.7113pt(1)}$
$\displaystyle:=\overrightarrow{x\,}_{n}+\bigg{(}\frac{1-e^{-\frac{1}{2}\gamma
h_{n}}}{\gamma}\bigg{)}\overrightarrow{v\,}_{n}^{\hskip
0.7113pt(1)}-\bigg{(}\frac{e^{-\frac{1}{2}\gamma h_{n}}+\frac{1}{2}\gamma
h_{n}-1}{\gamma^{2}}\bigg{)}u\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}\big{)}$ $\displaystyle\hskip
31.29802pt+\bigg{(}\frac{e^{-\frac{1}{2}\gamma h_{n}}+\frac{1}{2}\gamma
h_{n}-1}{\gamma^{2}h_{n}}\bigg{)}\sqrt{2\gamma u}\hskip
0.7113pt\big{(}W_{n}-12K_{n}\big{)}\hskip 0.7113pt,$
$\displaystyle\overrightarrow{x\,}_{n+1}$
$\displaystyle:=\overrightarrow{x\,}_{n}+\bigg{(}\frac{1-e^{-\gamma
h_{n}}}{\gamma}\bigg{)}\overrightarrow{v\,}_{n}^{\hskip
0.7113pt(1)}-\bigg{(}\frac{e^{-\gamma h_{n}}+\gamma
h_{n}-1}{\gamma^{2}}\bigg{)}\bigg{(}\frac{1}{3}u\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}\big{)}+\frac{2}{3}u\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}^{\hskip 0.7113pt(1)}\big{)}\bigg{)}$
$\displaystyle\hskip 31.29802pt+\bigg{(}\frac{e^{-\gamma h_{n}}+\gamma
h_{n}-1}{\gamma^{2}h_{n}}\bigg{)}\sqrt{2\gamma u}\hskip
0.7113pt\big{(}W_{n}-12K_{n}\big{)}\hskip 0.7113pt,$
$\displaystyle\overrightarrow{v\,}_{n}^{\hskip 0.7113pt(2)}$
$\displaystyle:=e^{-\gamma h_{n}}\overrightarrow{v\,}_{n}^{\hskip
0.7113pt(1)}-\frac{1}{6}\hskip 0.7113pte^{-\gamma h_{n}}u\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}\big{)}h_{n}-\frac{2}{3}\hskip
0.7113pte^{-\frac{1}{2}\gamma h_{n}}u\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}^{\hskip
0.7113pt(1)}\big{)}h_{n}-\frac{1}{6}u\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n+1}\big{)}h_{n}$ $\displaystyle\hskip
31.29802pt+\bigg{(}\frac{1-e^{-\gamma h_{n}}}{\gamma
h_{n}}\bigg{)}\sqrt{2\gamma u}\hskip 0.7113pt\big{(}W_{n}-12K_{n}\big{)},$
$\displaystyle\overrightarrow{v\,}_{n+1}$
$\displaystyle:=\overrightarrow{v\,}_{n}^{\hskip 0.7113pt(2)}-\sqrt{2\gamma
u}\hskip 0.7113pt\big{(}H_{n}-6K_{n}\big{)},$

where $\\{W_{n},H_{n},K_{n}\\}_{n\hskip 0.7113pt\geq\hskip 0.7113pt0}$ are the
Gaussian random vectors given in definition 1.1.

######  Definition 1.3 (The SOFA method).

Let $\\{t_{n}\\}_{n\hskip 0.7113pt\geq\hskip 0.7113pt0}$ be a sequence of
times with $t_{0}=0$, $t_{n+1}>t_{n}$ and step sizes
$h_{n}=t_{n+1}-t_{n}\hskip 0.7113pt$. We construct a numerical solution
$\\{(\hskip 0.7113pt\overrightarrow{x\,}_{n}\hskip
0.7113pt,\overrightarrow{v\,}_{n})\\}_{n\hskip 0.7113pt\geq\hskip 0.7113pt0}$
for (1.1) by setting $\hskip 0.7113pt(\hskip
0.7113pt\overrightarrow{x\,}_{0}\hskip
0.7113pt,\overrightarrow{v\,}_{0})\hskip 0.7113pt:=\hskip 0.7113pt(x_{0}\hskip
0.7113pt,v_{0})\hskip 0.7113pt$ and for $n\hskip 0.7113pt\geq\hskip
0.7113pt0$, defining $\hskip 0.7113pt(\hskip
0.7113pt\overrightarrow{x\,}_{n+1}\hskip
0.7113pt,\overrightarrow{v\,}_{n+1}\hskip 0.7113pt)\hskip 0.7113pt$ as
follows:

$\displaystyle\overrightarrow{v\,}_{n}^{\hskip 0.7113pt(0)}$
$\displaystyle:=\overrightarrow{v\,}_{n}+\sqrt{2\gamma u}\hskip
0.7113pt\big{(}H_{n}+6K_{n}\big{)},$
$\displaystyle\overrightarrow{v\,}_{n}^{\hskip 0.7113pt(1)}$
$\displaystyle:=e^{-\gamma(\frac{1}{2}+\phi)h_{n}}\overrightarrow{v\,}_{n}^{\hskip
0.7113pt(0)}+\Big{(}-u\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}\big{)}\hskip 0.7113pth_{n}+\sqrt{2\gamma
u}\hskip
0.7113pt\big{(}W_{n}-12K_{n}\big{)}\Big{)}\bigg{(}\frac{1-e^{-\gamma\left(\frac{1}{2}+\phi\right)h_{n}}}{\gamma
h_{n}}\bigg{)}\hskip 0.7113pt,$ $\displaystyle\overrightarrow{x\,}_{n}^{\hskip
0.7113pt(1)}$
$\displaystyle:=\overrightarrow{x\,}_{n}+\overrightarrow{v\,}_{n}^{\hskip
0.7113pt(1)}\big{(}1+2\phi\big{)}h_{n}\hskip 0.7113pt,$
$\displaystyle\overrightarrow{v\,}_{n}^{\hskip 0.7113pt(2)}$
$\displaystyle:=e^{\gamma\phi h_{n}}\overrightarrow{v\,}_{n}^{\hskip
0.7113pt(1)}+\Big{(}-u\nabla f\left(\overrightarrow{x\,}_{n}^{\hskip
0.7113pt(1)}\right)h_{n}+\sqrt{2\gamma u}\hskip
0.7113pt\big{(}W_{n}-12K_{n}\big{)}\Big{)}\bigg{(}\frac{1-e^{\gamma\phi
h_{n}}}{\gamma h_{n}}\bigg{)}\hskip 0.7113pt,$
$\displaystyle\overrightarrow{x\,}_{n}^{\hskip 0.7113pt(2)}$
$\displaystyle:=\overrightarrow{x\,}_{n}^{\hskip
0.7113pt(1)}-\overrightarrow{v\,}_{n}^{\hskip
0.7113pt(2)}\big{(}1+4\phi\big{)}h_{n}\hskip 0.7113pt,$
$\displaystyle\overrightarrow{v\,}_{n}^{\hskip 0.7113pt(3)}$
$\displaystyle:=e^{\gamma\phi h_{n}}\overrightarrow{v\,}_{n}^{\hskip
0.7113pt(2)}+\Big{(}-u\nabla f\left(\overrightarrow{x\,}_{n}^{\hskip
0.7113pt(2)}\right)h_{n}+\sqrt{2\gamma u}\hskip
0.7113pt\big{(}W_{n}-12K_{n}\big{)}\Big{)}\bigg{(}\frac{1-e^{\gamma\phi
h_{n}}}{\gamma h_{n}}\bigg{)}\hskip 0.7113pt,$
$\displaystyle\overrightarrow{x\,}_{n+1}$
$\displaystyle:=\overrightarrow{x\,}_{n}^{\hskip
0.7113pt(2)}+\overrightarrow{v\,}_{n}^{\hskip
0.7113pt(3)}\big{(}1+2\phi\big{)}h_{n}\hskip 0.7113pt,$
$\displaystyle\overrightarrow{v\,}_{n}^{\hskip 0.7113pt(4)}$
$\displaystyle:=e^{-\gamma(\frac{1}{2}+\phi)h_{n}}\overrightarrow{v\,}_{n}^{\hskip
0.7113pt(3)}+\Big{(}-u\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n+1}\big{)}\hskip 0.7113pth_{n}+\sqrt{2\gamma
u}\big{(}W_{n}-12K_{n}\big{)}\Big{)}\bigg{(}\frac{1-e^{-\gamma\left(\frac{1}{2}+\phi\right)h_{n}}}{\gamma
h_{n}}\bigg{)}\hskip 0.7113pt,$ $\displaystyle\overrightarrow{v\,}_{n+1}$
$\displaystyle:=\overrightarrow{v\,}_{n}^{\hskip 0.7113pt(4)}-\sqrt{2\gamma
u}\hskip 0.7113pt\big{(}H_{n}-6K_{n}\big{)},$

where $\\{W_{n},H_{n},K_{n}\\}_{n\hskip 0.7113pt\geq\hskip 0.7113pt0}$ are the
Gaussian vectors from definition 1.1 and $\phi$ is given by

$\displaystyle\phi:=\frac{-1+\sqrt[3]{2\hskip
0.7113pt}\,}{2\big{(}2-\sqrt[3]{2\hskip 0.7113pt}\,\big{)}}\hskip 0.7113pt.$

Despite being straightforward to implement, the above methods are difficult to
analyse. At the very least, it should be much easier to study the shifted ODE
method due to the similarities in structure between the SDE (1.1) and the ODE
(1.2). For example, both differential equations have continuous-time solutions
and admit similar Taylor expansions. As our analysis is based on these
continuous-time Taylor expansions, it does not apply if the ODE is
discretized. That said, we expect these methods to achieve a $2$-Wasserstein
error of $\varepsilon$ in $\mathcal{O}\big{(}\sqrt{d\hskip
0.7113pt}/\varepsilon^{\frac{1}{k}}\hskip 0.7113pt\big{)}$ steps when the
first $k$ derivatives of $f$ are Lipschitz continuous. In addition, we do not
study the numerical stability of the SORT and SOFA methods. Therefore it is
likely that further improvements can be made to the proposed ODE solvers
(particularly for large step sizes). For example, it may be possible to
discretize (1.2) using ideas from [6, 7], where the stability of splitting
methods is considered in the HMC setting.

Although this is not explored in the paper, it is also straightforward to
incorporate adaptive step sizes into these methodologies. That is, given
$(W_{n}\hskip 0.7113pt,H_{n}\hskip 0.7113pt,K_{n})$, we can generate the same
triple over the intervals $[t_{n},t_{n+\frac{1}{2}}]$ and
$[t_{n+\frac{1}{2}},t_{n+1}]$ where
$t_{n+\frac{1}{2}}:=\frac{1}{2}(t_{n}+t_{n+1})$. This procedure can be viewed
as an extension of Lévy’s construction of Brownian motion and is detailed in
[24]. We expect that adaptive step sizes are also possible for the other ULD
numerical methods (given in section 4). However this may be a topic for future
work.

### 1.2. Organization of the paper

Having introduced our proposed numerical methods, we now outline the paper. In
Section 2 we establish our key notation, definitions and assumptions. In
Section 3 we give an informal derivation of the shifted ODE method and present
our main results (which imply the $2$ -Wasserstein complexities given in Table
1.1). In Section 4, we discuss some related methods for ULD in preparation of
our numerical experiment. In Section 5, we empirically compare the SORT and
SOFA methods to these related methods in an example where the target density
comes from a logistic regression. The technical details for our analysis of
the shifted ODE method are given in the appendix.

##  2 Notation and definitions

In this section, we present some of the notation, definitions and assumptions
for the paper.

Throughout, we use $\|\cdot\|_{2}$ to denote the standard Euclidean norm on
$\mathbb{R}^{d}$ and $\mathbb{R}^{2d}$. Along with the Euclidean norm, we use
the usual inner product $\langle\hskip 0.7113pt\cdot\hskip 0.7113pt,\hskip
0.7113pt\cdot\hskip 0.7113pt\rangle$ with $\langle
v,v\rangle=\|v\|_{2}^{2}\hskip 0.7113pt$. For a function $g$, we use Big O
notation $O(g)$ and let $\mathcal{O}(g)$ be the class
$O(g)\cdot(\log(g))^{\hskip 0.7113ptO(1)}$.

### 2.1. Assumptions on $f$

The function $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is assumed to be
$m$-strongly convex

$\displaystyle f(y)\geq f(x)+\big{\langle}\hskip 0.7113pt\nabla
f(x),y-x\big{\rangle}+\frac{1}{2}m\|x-y\|_{2}^{2}\hskip 0.7113pt,$ (2.1)

and at least twice continuously differentiable with an $M$-Lipschitz
continuous gradient $\nabla f$,

$\displaystyle\|\nabla f(x)-\nabla f(y)\|_{2}\leq M\|x-y\|_{2}\hskip
0.7113pt,$ (2.2)

for all $x,y\in\mathbb{R}^{d}$. It is straightforward to show that the
conditions (2.1) and (2.2) are equivalent to the Hessian of $f$ being positive
definite and satisfying $mI_{d}\preccurlyeq\nabla^{2}f(x)\preccurlyeq
MI_{d}\hskip 0.7113pt$ (where $I_{d}$ is the $d\times d$ identity matrix and
$A\preccurlyeq B$ means that $x^{\mathsf{T}}A\hskip 0.7113ptx\leq
x^{\mathsf{T}}B\hskip 0.7113ptx$ for $x\in\mathbb{R}^{d}\hskip 0.7113pt$).

At certain points in our analysis, we make additional smoothness assumptions
on $f$. Before we give these assumptions, we first recall the notions of
tensors and operator norms. For $k\geq 1$, a $k$-tensor on $\mathbb{R}^{d}$ is
simply an element of the $d^{\hskip 0.7113ptk}$\- dimensional space
$\mathbb{R}^{d\times\cdots\times d}$. We can interpret a $k$-tensor on
$\mathbb{R}^{d}$ as a multilinear map into $\mathbb{R}$ with $k$ arguments
from $\mathbb{R}^{d}$. Hence for a $k$-tensor $T$ on $\mathbb{R}^{d}$ (with
$k\geq 2$) and a vector $v\in\mathbb{R}^{d}$, we define $Tv:=T(v,\cdots)$ as a
$(k-1)$-tensor on $\mathbb{R}^{d}$. Using this, we can define the operator
norm of $T$ recursively as

$\displaystyle\|T\|_{\text{op}}$ $\displaystyle:=\begin{cases}\,\,\,\hskip
0.7113pt\|T\|_{2}\hskip 0.7113pt,&\text{if}\,\,T\in\mathbb{R}^{d},\\\\[3.0pt]
\,\underset{\begin{subarray}{c}v\hskip 0.7113pt\in\hskip
0.7113pt\mathbb{R}^{d},\\\\[2.0pt] \|v\|_{2}\hskip 0.7113pt\leq\hskip
0.7113pt1\end{subarray}}{\sup}\|Tv\|_{\text{op}}\hskip
0.7113pt,&\text{if}\,\,T\,\,\text{is a }k\text{-tensor on
}\mathbb{R}^{d}\text{ with }k\geq 2\hskip 0.7113pt.\end{cases}$

With a slight abuse of notation, we will write $\|\cdot\|_{2}$ instead of
$\|\cdot\|_{\text{op}}$ for all $k$-tensors. In Section D.1, we additionally
assume that the Hessian of $f$ is $M_{2}$-Lipschitz continuous:

$\displaystyle\left\|\nabla^{2}f(x)-\nabla^{2}f(y)\right\|_{2}\leq
M_{2}\|x-y\|_{2}\hskip 0.7113pt,$ (2.3)

for all $x,y\in\mathbb{R}^{d\times d}$. In Section D.2, we shall additionally
assume that the function $f$ is three times continuously differentiable and
its third derivative is $M_{3}$-Lipschitz continuous:

$\displaystyle\left\|\nabla^{3}f(x)-\nabla^{3}f(y)\right\|_{2}\leq
M_{2}\|x-y\|_{2}\hskip 0.7113pt,$ (2.4)

for all $x,y\in\mathbb{R}^{d\times d\times d}$.

### 2.2. Probability notation

Suppose $\big{(}\Omega,\mathcal{F},\mathbb{P}\,;\\{\mathcal{F}_{t}\\}_{t\hskip
0.7113pt\geq\hskip 0.7113pt0}\big{)}$ is a filtered probability space carrying
a standard $d$-dimensional Brownian motion. The only SDE that we analyse in
this paper is (1.1), which admits a unique strong solution that is ergodic and
whose stationary measure $\pi$ has a density $\pi(x,v)\propto
e^{-f(x)+\frac{1}{2u}\|v\|_{2}^{2}}$ under our assumptions on $f$ (see
Proposition 6.1 in [50]). Numerical SDE solutions are obtained at times
$\\{t_{n}\\}_{n\hskip 0.7113pt\geq\hskip 0.7113pt0}$ with $t_{0}=0$,
$t_{n+1}>t_{n}$ and step sizes $h_{n}=t_{n+1}-t_{n}\hskip 0.7113pt$. The
results in Table 1.1 assume fixed steps. Given a random variable $X$, taking
its values in $\mathbb{R}^{d}$ or $\mathbb{R}^{2d}$, we define the
$\mathbb{L}_{p}$ norm of $X$ as

$\displaystyle\|X\|_{\mathbb{L}_{p}}:=\mathbb{E}\Big{[}\|X\|_{2}^{p}\Big{]}^{\frac{1}{p}},$

for $p\geq 1$. Similarly, for a stochastic process $\\{X_{t}\\}$, the
$\mathcal{F}_{t_{n}}$-conditional $\mathbb{L}_{p}$ norm of $X_{t}$ is

$\displaystyle\|X_{t}\|_{\mathbb{L}_{p}^{n}}:=\mathbb{E}_{n}\Big{[}\|X_{t}\|_{2}^{p}\Big{]}^{\frac{1}{p}}=\mathbb{E}\Big{[}\|X_{t}\|_{2}^{p}\,\big{|}\,\mathcal{F}_{t_{n}}\Big{]}^{\frac{1}{p}},$

for $t\geq t_{n}\hskip 0.7113pt$. For $0\leq s\leq t$, we use
$X_{s,t}:=X_{t}-X_{s}$ to denote the increment of $X$ over $[s,t]$.

### 2.3. Coupling and Wasserstein distance

Let $\mu$, $\nu$ be probability measures on $\mathbb{R}^{d}$. A coupling
between $\mu$ and $\nu$ is a random variable $Z=(X,Y)$ for which $X\sim\mu$
and $Y\sim\nu$. For $p\geq 1$, we define the $p$ -Wasserstein distance between
$\mu$ and $\nu$ as

$\displaystyle W_{p}\big{(}\mu,\nu\big{)}:=\inf_{(X,Y)\hskip
0.7113pt\sim\hskip 0.7113pt(\mu\times\nu)}\|X-Y\|_{\mathbb{L}_{p}}\hskip
0.7113pt,$ (2.5)

where the above infimum is taken over all couplings of random variables $X$
and $Y$ with distributions $X\sim\mu$ and $Y\sim\nu$. In this paper, we obtain
non-asymptotic bounds for the $2$ -Wasserstein distance between the target
distribution and shifted ODE approximation.

### 2.4. Diffusion and approximation processes

We shall use $W$ to denote a standard $d$-dimensional Brownian motion and
define the (piecewise) Brownian bridge process $B$ as

$\displaystyle
B_{t}=W_{t}-\bigg{(}W_{t_{n}}+\bigg{(}\frac{t-t_{n}}{h_{n}}\bigg{)}\hskip
0.7113ptW_{t_{n},t_{n+1}}\bigg{)},$ (2.6)

for $t\in[t_{n},t_{n+1}]$ and $n\geq 0$. Using $W$ and $B$, we can define the
following Gaussian vectors:

$\displaystyle W_{n}$ $\displaystyle:=W_{t_{n},t_{n+1}}\hskip 0.7113pt,$ (2.7)
$\displaystyle H_{n}$
$\displaystyle:=\frac{1}{h_{n}}\int_{t_{n}}^{t_{n+1}}B_{t_{n},t}\,dt,$ (2.8)
$\displaystyle K_{n}$
$\displaystyle:=\frac{1}{h_{n}^{2}}\int_{t_{n}}^{t_{n+1}}\bigg{(}\frac{1}{2}h_{n}-(t-t_{n})\bigg{)}B_{t_{n},t}\,dt.$
(2.9)

Since “$\hskip 0.7113pt12K_{n}$” appears in the shifted ODE, we also consider
the “shifted” Brownian bridge:

$\displaystyle\widetilde{B}_{t}:=B_{t}+12K_{n}\bigg{(}\frac{t-t_{n}}{h_{n}}\bigg{)},$

for $t\in[t_{n},t_{n+1})$ and $n\geq 0$. In the SDE (1.1), we use
$\sigma:=\sqrt{2\gamma u}$ to simply the notation. In our analysis of ULD, we
shall use the following diffusion and approximation processes:

1\. The underdamped Langevin diffusion $\\{(x_{t},v_{t})\\}_{t\hskip
0.7113pt\geq\hskip 0.7113pt0}$ is defined by the SDE,

$\displaystyle dx_{t}$ $\displaystyle=v_{t}\,dt,$ $\displaystyle dv_{t}$
$\displaystyle=-\gamma v_{t}\,dt-u\nabla f(x_{t})\,dt+\sigma\,dW_{t}\hskip
0.7113pt,$

with initial condition $(x_{0},v_{0})\sim p_{0}\hskip 0.7113pt$, where $p_{0}$
is a distribution on $\mathbb{R}^{2d}$. We usually set $p_{0}$ to the unique
invariant measure $\pi$ of ULD, which has density $\pi(x,v)\propto
e^{-f(x)+\frac{1}{2u}\|v\|_{2}^{2}}$. In this case, we have
$(x_{t},v_{t})\sim\pi$ for all $t\geq 0\hskip 0.7113pt$, and so the
$2$-Wasserstein distance between (the law of) an approximation
$\widetilde{x}_{n}$ at time $t_{n}$ and $\pi$ can be estimated using $\|\hskip
0.7113pt\widetilde{x}_{n}-x_{t_{n}}\|_{\mathbb{L}_{2}}$.

2\. Instead of using definition 1.1, it will be more convenient to define the
shifted ODE approximation
$\\{(\widetilde{x}_{n},\widetilde{v}_{n})\\}_{n\hskip 0.7113pt\geq\hskip
0.7113pt0}$ by setting $(\widetilde{x}_{0}\hskip
0.7113pt,\widetilde{v}_{0})\sim\widetilde{p}_{0}$ (for some distribution
$\widetilde{p}_{0}$ on $\mathbb{R}^{2d}\hskip 0.7113pt$) and defining
$(\widetilde{x}_{n+1},\widetilde{v}_{n+1})$ as

$\displaystyle\Bigg{(}\begin{matrix}\,\widetilde{x}_{n+1}\\\\[-3.0pt]
\,\widetilde{v}_{n+1}\end{matrix}\Bigg{)}:=\Bigg{(}\begin{matrix}\,\widehat{x}_{t_{n+1}}^{\hskip
0.7113ptn}\\\\[-3.0pt] \,\widehat{v}_{t_{n+1}}^{\hskip
0.7113ptn}\end{matrix}\Bigg{)}+12K_{n}\Bigg{(}\begin{matrix}\,0\\\\[-3.0pt]
\,\sigma\end{matrix}\,\Bigg{)},$

where $\big{\\{}\big{(}\widehat{x}_{t}^{\hskip
0.7113ptn},\widehat{v}_{t}^{\hskip 0.7113ptn}\big{)}\big{\\}}_{t\hskip
0.7113pt\in\hskip 0.7113pt[t_{n},t_{n+1}]}$ solves the following ODE,

$\displaystyle\frac{d}{dt}\Bigg{(}\begin{matrix}\,\widehat{x}^{\hskip
0.7113ptn}\\\\[-3.0pt] \,\widehat{v}^{\hskip
0.7113ptn}\end{matrix}\Bigg{)}=\Bigg{(}\begin{matrix}\,\widehat{v}^{\hskip
0.7113ptn}+\sigma\big{(}H_{n}+6K_{n}\big{)}\\\\[-3.0pt] \,-\gamma\big{(}\hskip
0.7113pt\widehat{v}^{\hskip
0.7113ptn}+\sigma\big{(}H_{n}+6K_{n}\big{)}\big{)}-u\nabla
f\big{(}\,\widehat{x}^{\hskip
0.7113ptn}\big{)}\end{matrix}\,\Bigg{)}+\frac{W_{n}-12K_{n}}{h_{n}}\,\Bigg{(}\begin{matrix}\,0\\\\[-3.0pt]
\,\sigma\end{matrix}\,\Bigg{)}\hskip 0.7113pt,$

with initial condition $\big{(}\hskip 0.7113pt\widehat{x}_{t_{n}}^{\hskip
0.7113ptn},\widehat{v}_{t_{n}}^{\hskip
0.7113ptn}\big{)}:=\big{(}\widetilde{x}_{n}\hskip
0.7113pt,\widetilde{v}_{n}\big{)}$.

3\. We also consider a sequence of shifted ODE approximations
$\\{(x_{n}^{\prime},v_{n}^{\prime})\\}_{n\hskip 0.7113pt\geq\hskip 0.7113pt0}$
given by

$\displaystyle\Bigg{(}\begin{matrix}\,x_{n+1}^{\prime}\\\\[-3.0pt]
\,v_{n+1}^{\prime}\end{matrix}\Bigg{)}:=\Bigg{(}\begin{matrix}\,\wideparen{x}^{\hskip
0.7113ptn}_{t_{n+1}}\\\\[-3.0pt] \,\wideparen{v}^{\hskip
0.7113ptn}_{t_{n+1}}\end{matrix}\Bigg{)}+12K_{n}\Bigg{(}\begin{matrix}0\\\\[-3.0pt]
\sigma\end{matrix}\,\Bigg{)},$

where $\big{\\{}\big{(}\wideparen{x}^{\hskip
0.7113ptn}_{t},\wideparen{v}^{\hskip
0.7113ptn}_{t}\big{)}\big{\\}}_{t\in[t_{n},t_{n+1}]}$ solves the following
ODE,

$\displaystyle\frac{d}{dt}\Bigg{(}\begin{matrix}\,\wideparen{x}^{\hskip
0.7113ptn}\\\\[-3.0pt] \,\wideparen{v}^{\hskip
0.7113ptn}\end{matrix}\Bigg{)}=\Bigg{(}\begin{matrix}\,\wideparen{v}^{\hskip
0.7113ptn}+\sigma\big{(}H_{n}+6K_{n}\big{)}\\\\[-3.0pt] \,-\gamma\big{(}\hskip
0.7113pt\wideparen{v}^{\hskip
0.7113ptn}+\sigma\big{(}H_{n}+6K_{n}\big{)}\big{)}-u\nabla
f\big{(}\,\wideparen{x}^{\hskip
0.7113ptn}\big{)}\end{matrix}\,\Bigg{)}+\frac{W_{n}-12K_{n}}{h_{n}}\,\Bigg{(}\begin{matrix}\,0\\\\[-3.0pt]
\,\sigma\end{matrix}\,\Bigg{)}\hskip 0.7113pt,$

with initial condition $\big{(}\hskip 0.7113pt\wideparen{x}^{\hskip
0.7113ptn}_{t_{n}}\hskip 0.7113pt,\wideparen{v}^{\hskip
0.7113ptn}_{t_{n}}\big{)}\hskip 0.7113pt:=\hskip 0.7113pt\big{(}x_{n}\hskip
0.7113pt,v_{n}\big{)}$. In addition, we set $(\hskip
0.7113ptx_{0}^{\prime}\hskip 0.7113pt,v_{0}^{\prime}\hskip 0.7113pt)\hskip
0.7113pt:=\hskip 0.7113pt(x_{0}\hskip 0.7113pt,v_{0})\hskip 0.7113pt\sim\hskip
0.7113ptp_{0}\hskip 0.7113pt$.

4\. Since we use the same $2$-Wasserstein contractivity arguments as those
given in [16], we define certain linear combinations of
$\\{(x_{t_{n}},v_{t_{n}})\\}_{n\hskip 0.7113pt\geq\hskip 0.7113pt0}\hskip
0.7113pt$, $\\{(\widetilde{x}_{n},\widetilde{v}_{n})\\}_{n\hskip
0.7113pt\geq\hskip 0.7113pt0}$ and
$\\{(x_{n}^{\prime},v_{n}^{\prime})\\}_{n\hskip 0.7113pt\geq\hskip
0.7113pt0}\hskip 0.7113pt$:

$\displaystyle y_{n}$
$\displaystyle:=\Bigg{(}\begin{matrix}\,\big{(}\lambda\hskip
0.7113pt\widetilde{x}_{n}+\widetilde{v}_{n}\big{)}-\big{(}\lambda\hskip
0.7113ptx_{t_{n}}+v_{t_{n}}\big{)}\\\\[-3.0pt] \,\big{(}\eta\hskip
0.7113pt\widetilde{x}_{n}+\widetilde{v}_{n}\big{)}-\big{(}\eta\hskip
0.7113ptx_{t_{n}}+v_{t_{n}}\big{)}\end{matrix}\Bigg{)},$
$\displaystyle\widetilde{y}_{n}$
$\displaystyle:=\Bigg{(}\begin{matrix}\,\big{(}\lambda\hskip
0.7113pt\widetilde{x}_{n}+\widetilde{v}_{n}\big{)}-\big{(}\lambda\hskip
0.7113ptx_{n}^{\prime}+v_{n}^{\prime}\big{)}\\\\[-3.0pt] \,\big{(}\eta\hskip
0.7113pt\widetilde{x}_{n}+\widetilde{v}_{n}\big{)}-\big{(}\eta\hskip
0.7113ptx_{n}^{\prime}+v_{n}^{\prime}\big{)}\end{matrix}\Bigg{)},$
$\displaystyle y_{n}^{\prime}$
$\displaystyle:=\Bigg{(}\begin{matrix}\,\big{(}\lambda\hskip
0.7113ptx_{n}^{\prime}+v_{n}^{\prime}\big{)}-\big{(}\lambda\hskip
0.7113ptx_{t_{n}}+v_{t_{n}}\big{)}\\\\[-3.0pt] \,\big{(}\eta\hskip
0.7113ptx_{n}^{\prime}+v_{n}^{\prime}\big{)}-\big{(}\eta\hskip
0.7113ptx_{t_{n}}+v_{t_{n}}\big{)}\end{matrix}\Bigg{)},$

where $\lambda\in[0,\frac{1}{2}\gamma)$ and $\eta:=\gamma-\lambda\hskip
0.7113pt$.

We use the notation $\big{\\{}\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}\hskip
0.7113pt,\overrightarrow{v\,}_{n}\big{)}\big{\\}}_{n\hskip 0.7113pt\geq\hskip
0.7113pt0}$ to denote the other numerical solutions of ULD (SORT, SOFA, OBABO,
Strang splitting, randomized midpoint and left-point methods). ULD and its
various approximations will be defined from the same Brownian motion $W$.

##  3 Derivation and error analysis of the shifted ODE method

In this section, we shall discuss the theoretical aspects of the shifted ODE
method and present our main results regarding the $2$-Wasserstein convergence
of (the law of) $\widetilde{x}_{n}$ to $\pi$. For general SDEs, it is well
known that the best $\mathbb{L}_{2}$ approximations converge with rate
$O\big{(}N^{-\frac{1}{2}}\hskip 0.7113pt\big{)}$ where $N$ is the number of
Gaussian random variables used by the method [13, 17]. Hence we begin the
section by explaining why high order methods are possible for ULD.

### 3.1. The stochastic Taylor expansion of ULD

Perhaps the most important tool when developing and analysing numerical
methods for SDEs is the stochastic Taylor expansion. Just like for ODEs, this
allows one to express the solution as a sum of “polynomial” and “remainder”
terms. The challenge for SDEs is that these “polynomial” terms can now include
certain iterated integrals of Brownian motion that are difficult to generate
[27, 25]. By applying the Itô-Taylor expansion (Theorem 5.5.1 in [32]) to the
Langevin SDE (1.1), we see the significant terms are either monomials in
$h=t-s$ or Gaussian integrals of $W$.

######  Theorem 3.1 (High order Taylor expansion of ULD).

Consider the SDE (1.1) and suppose that the potential $f$ is three times
continuously differentiable. Then for $0\leq s\leq t$,

$\displaystyle x_{t}=x_{s}$ $\displaystyle+v_{s}\,h-\frac{1}{2}\big{(}\gamma
v_{s}+u\nabla
f(x_{s})\big{)}h^{2}+\sigma\int_{s}^{t}W_{s,r}\,dr-\sigma\gamma\int_{s}^{t}\int_{s}^{r_{1}}W_{s,{r_{2}}}\,dr_{2}\,dr_{1}$
(3.1) $\displaystyle+\frac{1}{6}\big{(}\gamma^{2}v_{s}+\gamma u\nabla
f(x_{s})-u\nabla^{2}f(x_{s})v_{s}\big{)}h^{3}$
$\displaystyle+\sigma\big{(}\gamma^{2}-u\nabla^{2}f(x_{s})\big{)}\int_{s}^{t}\int_{s}^{r_{1}}\int_{s}^{r_{2}}W_{s,r_{3}}\,dr_{3}\,dr_{2}\,dr_{1}$
$\displaystyle-\frac{1}{24}\big{(}\big{(}\gamma^{2}-u\nabla^{2}f(x_{s})\big{)}\big{(}\gamma
v_{s}+u\nabla
f(x_{s})\big{)}-u\gamma\nabla^{2}f(x_{s})v_{s}+u\nabla^{3}f(x_{s})(v_{s},v_{s})\big{)}h^{4}$
$\displaystyle+R^{\hskip 0.7113ptx}(h,x_{s},v_{s})\hskip 0.7113pt,$
$\displaystyle v_{t}=v_{s}$ $\displaystyle-\big{(}\gamma v_{s}+u\nabla
f(x_{s})\big{)}h+\sigma W_{s,t}-\sigma\gamma\int_{s}^{t}W_{s,r}\,dr$ (3.2)
$\displaystyle+\frac{1}{2}\big{(}\gamma^{2}v_{s}+u\gamma\nabla
f(x_{s})-u\nabla^{2}f(x_{s})v_{s}\big{)}h^{2}+\sigma\big{(}\gamma^{2}-u\nabla^{2}f(x_{s})\big{)}\int_{s}^{t}\int_{s}^{r_{1}}W_{s,r_{2}}\,dr_{2}\,dr_{1}$
$\displaystyle-\frac{1}{6}\big{(}\big{(}\gamma^{2}-u\nabla^{2}f(x_{s})\big{)}\big{(}\gamma
v_{s}+u\nabla
f(x_{s})\big{)}-u\gamma\nabla^{2}f(x_{s})v_{s}+u\nabla^{3}f(x_{s})(v_{s},v_{s})\big{)}h^{3}$
$\displaystyle-\sigma\gamma\big{(}\gamma^{2}-u\nabla^{2}f(x_{s})\big{)}\int_{s}^{t}\int_{s}^{r_{1}}\int_{s}^{r_{2}}W_{s,r_{3}}\,dr_{3}\,dr_{2}\,dr_{1}$
$\displaystyle-\sigma
u\int_{s}^{t}\int_{s}^{r_{1}}\nabla^{3}f(x_{s})\big{(}v_{s},(r_{2}-s)W_{s,r_{2}}\big{)}\,dr_{2}\,dr_{1}+R^{\hskip
0.7113ptv}(h,x_{s},v_{s}).$

where $h=t-s$ and the remainder terms $R^{\hskip 0.7113ptx}(h,x_{s},v_{s})$,
$R^{\hskip 0.7113ptv}(h,x_{s},v_{s})$ are given by

$\displaystyle R^{\hskip 0.7113ptx}(h,x_{s},v_{s})$
$\displaystyle=\gamma\big{(}u\nabla^{2}f(x_{s})-\gamma^{2}\big{)}\int_{s}^{t}\int_{s}^{r_{1}}\int_{s}^{r_{2}}\int_{s}^{r_{3}}(v_{r_{4}}-v_{s})\,dr_{4}\,dr_{3}\,dr_{2}\,dr_{1}$
(3.3) $\displaystyle\hskip
14.22636pt+u\big{(}u\nabla^{2}f(x_{s})-\gamma^{2}\big{)}\int_{s}^{t}\int_{s}^{r_{1}}\int_{s}^{r_{2}}\int_{s}^{r_{3}}\big{(}\hskip
0.7113pt\nabla f(x_{r_{4}})-\nabla
f(x_{s})\big{)}\,dr_{4}\,dr_{3}\,dr_{2}\,dr_{1}$ $\displaystyle\hskip
14.22636pt+u\gamma\int_{s}^{t}\int_{s}^{r_{1}}\int_{s}^{r_{2}}\int_{s}^{r_{3}}\big{(}\hskip
0.7113pt\nabla^{2}f(x_{r_{4}})v_{r_{4}}-\nabla^{2}f(x_{s})v_{s}\big{)}\,dr_{4}\,dr_{3}\,dr_{2}\,dr_{1}$
$\displaystyle\hskip
14.22636pt-u\int_{s}^{t}\int_{s}^{r_{1}}\int_{s}^{r_{2}}\int_{s}^{r_{3}}\big{(}\hskip
0.7113pt\nabla^{3}f(x_{r_{4}})(v_{r_{4}},v_{r_{3}})-\nabla^{3}f(x_{s})(v_{s},v_{s})\big{)}\,dr_{4}\,dr_{3}\,dr_{2}\,dr_{1}\hskip
0.7113pt,$ $\displaystyle R^{\hskip 0.7113ptv}(h,x_{s},v_{s})$
$\displaystyle=\gamma\big{(}u\nabla^{2}f(x_{s})-\gamma^{2}\big{)}\int_{s}^{t}\int_{s}^{r_{1}}\int_{s}^{r_{2}}\big{(}v_{r_{3}}-v_{s}-\sigma
W_{s,r_{3}}\big{)}\,dr_{3}\,dr_{2}\,dr_{1}$ (3.4) $\displaystyle\hskip
14.22636pt+u\big{(}u\nabla^{2}f(x_{s})-\gamma^{2}\big{)}\int_{s}^{t}\int_{s}^{r_{1}}\int_{s}^{r_{2}}\big{(}\nabla
f(x_{r_{3}})-\nabla f(x_{s})\big{)}\,dr_{3}\,dr_{2}\,dr_{1}$
$\displaystyle\hskip
14.22636pt+u\gamma\int_{s}^{t}\int_{s}^{r_{1}}\int_{s}^{r_{2}}\big{(}\nabla^{2}f(x_{r_{3}})v_{r_{3}}-\nabla^{2}f(x_{s})v_{s}\big{)}\,dr_{3}\,dr_{2}\,dr_{1}$
$\displaystyle\hskip
14.22636pt-u\int_{s}^{t}\int_{s}^{r_{1}}\int_{s}^{r_{2}}\big{(}\hskip
0.7113pt\nabla^{3}f(x_{r_{3}})(v_{r_{3}},v_{r_{2}})-\nabla^{3}f(x_{s})(v_{s},v_{s})\big{)}\,dr_{3}\,dr_{2}\,dr_{1}$
$\displaystyle\hskip 14.22636pt-\sigma
u\int_{s}^{t}\int_{s}^{r_{1}}\int_{s}^{r_{2}}\big{(}\hskip
0.7113pt\nabla^{3}f(x_{r_{3}})\big{(}v_{r_{3}},W_{s,r_{2}}\big{)}-\nabla^{3}f(x_{s})\big{(}v_{s},W_{s,r_{2}}\big{)}\big{)}\,dr_{3}\,dr_{2}\,dr_{1}\hskip
0.7113pt.$

###### Proof.

For $\tau\in[s,t]$, the SDE (1.1) can be written in integral form as

$\displaystyle x_{\tau}$ $\displaystyle=x_{s}+\int_{s}^{\tau}v_{r}\,dr,$ (3.5)
$\displaystyle v_{\tau}$
$\displaystyle=v_{s}-\gamma\int_{s}^{\tau}v_{r}\,dr-u\int_{s}^{\tau}\nabla
f(x_{r})\,dr+\sigma W_{s,\tau}\hskip 0.7113pt.$ (3.6)

Moreover, it then follows by Itô’s lemma that

$\displaystyle\nabla f(x_{r_{1}})$ $\displaystyle=\nabla
f(x_{s})+\int_{s}^{r_{1}}\nabla^{2}f(x_{r_{2}})\,dx_{r_{2}}$
$\displaystyle=\nabla
f(x_{s})+\int_{s}^{r_{1}}\nabla^{2}f(x_{r_{2}})v_{r_{2}}\,dr_{2}\hskip
0.7113pt,$ (3.7) $\displaystyle\nabla^{2}f(x_{r_{2}})v_{\tau}$
$\displaystyle=\nabla^{2}f(x_{s})v_{\tau}+\int_{s}^{r_{2}}\nabla^{3}f(x_{r_{3}})\big{(}dx_{r_{3}},v_{\tau}\big{)}$
$\displaystyle=\nabla^{2}f(x_{s})v_{\tau}+\int_{s}^{r_{2}}\nabla^{3}f(x_{r_{3}})(v_{r_{3}},v_{\tau})\,dr_{3}\hskip
0.7113pt.$ (3.8)

for $s\leq r_{1}\leq r_{2}\leq t$. The result now follows by repeatedly
substituting the four identities (3.5) – (3.8) into the integral equations
(3.5) – (3.6) and rearranging the various terms. ∎

######  Remark 3.2.

As a quick sanity check, we see that by integrating the lower order terms in
the expansion of $v$ over $[s,t]$, we obtain the corresponding terms in the
expansion of $x$.

By ignoring the remainder terms from the expansions (3.2) and (3.1), we can
derive a stochastic Taylor approximation for the diffusion $(x,v)$ at time $t$
given its value at time $s$. However this would require us to compute or
estimate the derivatives $\nabla^{2}f(x_{s}),\nabla^{3}f(x_{s})$ and in
applications, this may become computationally expensive (especially if $d$ is
high). That said, this does indicate that “third order” methods can be
implemented for ULD. Hence we seek a numerical method that captures the Taylor
expansions in Theorem 3.1 but only requires evaluations of $\nabla f$ (and is
therefore “derivative-free” in a certain sense).

### 3.2. Derivation of the shifted ODE method

Our key insight (which is really an insight from rough path theory [41]) is
that SDEs are a particular example of so-called “controlled differential
equations” (CDEs). In our setting, we consider the following CDE:

$\displaystyle d\hskip 0.7113pt\overline{x}_{t}$
$\displaystyle=\overline{v}_{t}\,dt,$ (3.9) $\displaystyle d\hskip
0.7113pt\overline{v}_{t}$ $\displaystyle=-\gamma\hskip
0.7113pt\overline{v}_{t}\,dt-u\nabla f(\hskip
0.7113pt\overline{x}_{t})\,dt+\sigma\,dX_{t}\hskip 0.7113pt,$

where $X:[0,T]\rightarrow\mathbb{R}^{d}$ is continuous and has bounded
variation so that (3.9) is well defined. Moreover, by constructing the path
$X$ using information about the Brownian motion $W$, we can obtain a CDE
solution that approximates the SDE solution. To understand how $X$ should be
constructed, we first note that (3.9) admits a Taylor expansion analogues to
(3.1) – (3.2). The only difference is that the chain rule is employed instead
of Itô’s lemma.

######  Theorem 3.3 (Taylor expansion of underdamped Langevin CDE).

Consider the CDE (3.9) and suppose $f$ is three times continuously
differentiable. Then for $0\leq s\leq t$,

$\displaystyle\overline{x}_{t}=\overline{x}_{s}$
$\displaystyle+\overline{x}_{s}\,h-\frac{1}{2}\big{(}\cdots\big{)}h^{2}+\sigma\int_{s}^{t}X_{s,r}\,dr-\sigma\gamma\int_{s}^{t}\int_{s}^{r_{1}}X_{s,{r_{2}}}\,dr_{2}\,dr_{1}$
(3.10)
$\displaystyle+\frac{1}{6}\big{(}\cdots\big{)}h^{3}+\sigma\big{(}\cdots\big{)}\int_{s}^{t}\int_{s}^{r_{1}}\int_{s}^{r_{2}}X_{s,r_{3}}\,dr_{3}\,dr_{2}\,dr_{1}-\frac{1}{24}\big{(}\cdots\big{)}h^{4}+R^{\hskip
0.7113pt\overline{x}}(h,\overline{x}_{s},\overline{v}_{s})\hskip 0.7113pt,$
$\displaystyle\overline{v}_{t}=\overline{v}_{s}$
$\displaystyle-\big{(}\gamma\hskip 0.7113pt\overline{v}_{s}+u\nabla f(\hskip
0.7113pt\overline{x}_{s})\big{)}h+\sigma
X_{s,t}-\sigma\gamma\int_{s}^{t}X_{s,r}\,dr+\frac{1}{2}\big{(}\cdots\big{)}h^{2}$
(3.11)
$\displaystyle+\sigma\big{(}\cdots\big{)}\int_{s}^{t}\int_{s}^{r_{1}}X_{s,r_{2}}\,dr_{2}\,dr_{1}-\frac{1}{6}\big{(}\cdots\big{)}h^{3}-\sigma\gamma\big{(}\cdots\big{)}\int_{s}^{t}\int_{s}^{r_{1}}\int_{s}^{r_{2}}X_{s,r_{3}}\,dr_{3}\,dr_{2}\,dr_{1}$
$\displaystyle-\sigma u\int_{s}^{t}\int_{s}^{r_{1}}\nabla^{3}f(\hskip
0.7113pt\overline{x}_{s})\big{(}\hskip
0.7113pt\overline{v}_{s},(r_{2}-s)X_{s,r_{2}}\big{)}\,dr_{2}\,dr_{1}+R^{\hskip
0.7113pt\overline{v}}(h,\overline{x}_{s},\overline{v}_{s}),$

where $h=t-s$, the coefficients $(\hskip 0.7113pt\cdots)$ are given in Theorem
3.1 and the remainder terms $R^{\hskip
0.7113pt\overline{x}}(h,\overline{x}_{s},\overline{v}_{s}),R^{\hskip
0.7113pt\overline{v}}(h,\overline{x}_{s},\overline{v}_{s})$ have the same form
as their stochastic counterparts (3.3), (3.4).

Thus, to approximate the SDE solution, we construct $X$ over the interval
$[s,t]$ so that

$\displaystyle 1.$ $\displaystyle\,\,\,\,X_{s,t}=W_{s,t}\hskip 0.7113pt,$
$\displaystyle 2.$
$\displaystyle\,\int_{s}^{t}X_{s,r}\,dr=\int_{s}^{t}W_{s,r}\,dr\hskip
0.7113pt,$ $\displaystyle 3.$
$\displaystyle\,\int_{s}^{t}\int_{s}^{r_{1}}X_{s,r_{2}}\,dr_{2}\,dr_{1}=\int_{s}^{t}\int_{s}^{r_{1}}W_{s,r_{2}}\,dr_{2}\,dr_{1}\hskip
0.7113pt.$

At this point, it is worth noting that the above three properties are not
sufficient for accurately approximating general SDEs. This is because a
general SDE will contain additional non-Gaussian iterated integrals of
Brownian motion in its Taylor expansion. The two most notable examples of such
integrals (neither of which appears for ULD) are

$\displaystyle\int_{s}^{t}W_{s,r}^{(i)}\,dW_{r}^{(j)},\hskip
14.22636pt\int_{s}^{t}\int_{s}^{r_{1}}W_{s,r_{2}}\,dW_{r_{2}}\,dr_{1}\hskip
0.7113pt,$

for $i,j\in\\{1,\cdots,d\\}$ with $i\neq j$. Due to the specific structure of
the Langevin SDE/CDE, we do not need to take these integrals into
consideration when constructing the path $X$.

To make the CDE (3.9) easier to discretize, we will assume that $X$ is
piecewise linear. Moreover, the path should have as few pieces as possible (to
reduce computational cost). The following theorem shows that, in general, it
is not enough for $X$ to have two pieces.

######  Theorem 3.4.

Let $a\in(s,t)$ and $b,c\in\mathbb{R}^{d}$. Suppose that the piecewise linear
path $X$ is defined on $[s,t]$ by $X_{r}:=b\big{(}\frac{r-s}{a-s}\big{)}$ for
$r\in[s,a]$ and $X_{r}:=b+(c-b)\big{(}\frac{r-a}{t-a}\big{)}$ for $r\in[a,t]$.
Then $X_{s,t}=c$, $\int_{s}^{t}X_{s,r}\,dr=\frac{1}{2}(b+c)(t-s)$ and
$\int_{s}^{t}\int_{s}^{r_{1}}X_{s,r_{2}}\,dr_{2}\,dr_{1}$ depends linearly on
$b,c$ and quadratically on $\hskip 0.7113pta$. In particular, for every
$b,\hskip 0.7113ptc\in\mathbb{R}^{d}$, it is bounded as a function of $a$.

###### Proof.

Using integration by parts, we can express the iterated integral of $X$ as

$\displaystyle\int_{s}^{t}\int_{s}^{r_{1}}X_{s,r_{2}}\,dr_{2}\,dr_{1}=(t-s)\int_{s}^{t}X_{s,r}\,dr-\int_{s}^{t}(r-s)\hskip
0.7113ptX_{s,r}\,dr\hskip 0.7113pt.$

By evaluating the above, we see that the integral is linear in $b,c$ and
quadratic in $a$. ∎

Therefore, the path $X$ should be a piecewise linear path on $[s,t]$ with
three pieces and thus defined by five variables: the discretization values
$b,c,d\in\mathbb{R}^{d}$ and times $a_{1},a_{2}\in(s,t)$. Our key observation
is that we can set $a_{1}=s$, $a_{2}=t$. That is, we give $X$ “vertical
pieces”. Since the Langevin equation has additive noise, it is trivial to
solve the CDE along these vertical pieces and so (3.9) reduces to a single ODE
(corresponding to the middle piece).

Figure 3.1: Brownian motion approximated using piecewise linear paths with
vertical pieces.

Finally to derive the shifted ODE method, we determine the unknown variables
$b,c,d$. To simply calculations, we express the integrals of $W$ using
independent random vectors.

######  Theorem 3.5.

Let $W$ denote a Brownian motion. With $h=t-s$, we define the integrals:

$\displaystyle H_{s,t}$
$\displaystyle:=\frac{1}{h}\int_{s}^{t}\bigg{(}W_{s,r}-\frac{r-s}{h}W_{s,t}\bigg{)}\,dr,$
(3.12) $\displaystyle K_{s,t}$
$\displaystyle:=\frac{1}{h^{2}}\int_{s}^{t}\bigg{(}\frac{1}{2}h-(r-s)\bigg{)}\bigg{(}W_{s,r}-\frac{r-s}{h}W_{s,t}\bigg{)}\,dr.$
(3.13)

Then $W_{s,t}\sim\mathcal{N}\big{(}0,hI_{d}\big{)}$,
$H_{s,t}\sim\mathcal{N}\big{(}0,\frac{1}{12}hI_{d}\big{)}$,
$K_{s,t}\sim\mathcal{N}\big{(}0,\frac{1}{720}hI_{d}\big{)}$ are independent
and

$\displaystyle\int_{s}^{t}W_{s,r}\,dr$
$\displaystyle=\frac{1}{2}hW_{s,t}+hH_{s,t}\hskip 0.7113pt,$ (3.14)
$\displaystyle\int_{s}^{t}\int_{s}^{r_{1}}W_{s,r_{2}}\,dr_{2}\,dr_{1}$
$\displaystyle=\frac{1}{6}h^{2}W_{s,t}+\frac{1}{2}h^{2}H_{s,t}+h^{2}K_{s,t}\hskip
0.7113pt.$ (3.15)

###### Proof.

The first identity is trivial and the second identity follows by a direct
calculation (involving integration by parts). For each $i\in\\{1,\cdots,d\\}$,
we can identify $W_{s,t}^{(i)},H_{s,t}^{(i)},K_{s,t}^{(i)}$ as rescaled
coefficients from the polynomial expansion of $W^{(i)}$ (see Theorem 2.3 of
[25]). In particular, this means that $W_{s,t},H_{s,t},K_{s,t}$ are
independent centered Gaussian vectors. The variances of both $H_{0,1}$ and
$K_{0,1}$ can be directly computed using Fubini’s theorem and the Brownian
bridge covariance function, $\mathbb{E}[B_{r_{1}}B_{r_{2}}]=\min(r_{1}\hskip
0.7113pt,r_{2})-r_{1}r_{2}$ for $r_{1}\hskip 0.7113pt,r_{2}\in[0,1]$.

$\displaystyle\operatorname{Var}\big{(}H_{0,1}\big{)}$
$\displaystyle=\int_{0}^{1}\int_{0}^{1}\mathbb{E}\big{[}B_{r_{1}}B_{r_{2}}\big{]}\,dr_{2}\,dr_{1}=\frac{1}{12}\hskip
0.7113pt,$ $\displaystyle\operatorname{Var}\big{(}K_{0,1}\big{)}$
$\displaystyle=\int_{0}^{1}\int_{0}^{1}\bigg{(}\frac{1}{2}-r_{1}\bigg{)}\bigg{(}\frac{1}{2}-r_{2}\bigg{)}\mathbb{E}\big{[}B_{r_{1}}B_{r_{2}}\big{]}\,dr_{2}\,dr_{1}=\frac{1}{720}\hskip
0.7113pt.$

The result follows by the natural Brownian scaling
$(H_{s,t}\sim\sqrt{h}H_{0,1}\hskip 0.7113pt,K_{s,t}\sim\sqrt{h}K_{0,1})$. ∎

Using the identities from Theorem 3.5, we can derive a simple formula for the
path $X$.

######  Theorem 3.6.

Let $0\leq s<t$ and $h=t-s$. We define a path
$X:[s,t]\rightarrow\mathbb{R}^{d}$ as follows:

$\displaystyle
X_{r}:=\begin{cases}\big{(}H_{s,t}+6K_{s,t}\big{)}+\big{(}W_{s,t}-12K_{s,t}\big{)}\displaystyle\frac{r-s}{h}\hskip
0.7113pt,&\text{for}\,\,\,r\in(s,t)\\\ \,\,W_{r}\hskip
0.7113pt,&\text{for}\,\,r\in\\{s,t\\}\end{cases}.$

Then $\displaystyle X_{s,t}=W_{s,t}\hskip
0.7113pt,\int_{s}^{t}X_{s,r}\,dr=\int_{s}^{t}W_{s,r}\,dr\,$ and
$\displaystyle\int_{s}^{t}\int_{s}^{r_{1}}X_{s,r_{2}}\,dr_{2}\,dr_{1}=\int_{s}^{t}\int_{s}^{r_{1}}W_{s,r_{2}}\,dr_{2}\,dr_{1}\hskip
0.7113pt.$

###### Proof.

The result follows by a direct calculation with the identities (3.14) and
(3.15). ∎

By driving the CDE (3.9) with the above path, we obtain the shifted ODE method
(definition 1.1). The ODE (1.2) defined on $[0,1]$ is then obtained by a
change of variables.

### 3.3. Derivation of the SORT method

For simplicity, we consider the following ODE:

$\displaystyle\frac{d}{dr}\Bigg{(}\begin{matrix}\hskip
0.7113ptx_{r}^{n}\\\\[-4.0pt] \hskip
0.7113ptv_{r}^{n}\end{matrix}\Bigg{)}=\Bigg{(}\begin{matrix}\,v_{r}^{n}\\\\[-4.0pt]
\,-u\nabla f\big{(}x_{r}^{n}\big{)}\end{matrix}\Bigg{)},$ (3.16)

on $[t_{n},t_{n+1}]$. Then by applying the standard third order Runge-Kutta
method [5], we have

$\displaystyle\overrightarrow{x\,}_{n}^{\hskip 0.7113pt(1)}$
$\displaystyle:=\overrightarrow{x\,}_{n}+\frac{1}{2}\overrightarrow{v\,}_{n}\hskip
0.7113pth_{n}\hskip 0.7113pt,$ $\displaystyle\overrightarrow{v\,}_{n}^{\hskip
0.7113pt(1)}$ $\displaystyle:=\overrightarrow{v\,}_{n}-\frac{1}{2}u\nabla
f\big{(}\overrightarrow{x\,}_{n}\big{)}\hskip 0.7113pth_{n}\hskip 0.7113pt,$
$\displaystyle\overrightarrow{x\,}_{n}^{\hskip 0.7113pt(2)}$
$\displaystyle:=\overrightarrow{x\,}_{n}-\overrightarrow{v\,}_{n}\hskip
0.7113pth_{n}+2\overrightarrow{v\,}_{n}^{\hskip 0.7113pt(1)}\hskip
0.7113pth_{n}\hskip 0.7113pt,$ $\displaystyle\overrightarrow{v\,}_{n}^{\hskip
0.7113pt(2)}$ $\displaystyle:=\overrightarrow{v\,}_{n}+u\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}\big{)}h_{n}-2u\nabla
f\big{(}\overrightarrow{x\,}_{n}^{\hskip 0.7113pt(1)}\big{)}\hskip
0.7113pth_{n}\hskip 0.7113pt,$ $\displaystyle\overrightarrow{x\,}_{n+1}$
$\displaystyle:=\overrightarrow{x\,}_{n}+\frac{1}{6}\overrightarrow{v\,}_{n}\hskip
0.7113pth_{n}+\frac{2}{3}\overrightarrow{v\,}_{n}^{\hskip 0.7113pt(1)}\hskip
0.7113pth_{n}+\frac{1}{6}\overrightarrow{v\,}_{n}^{\hskip 0.7113pt(2)}\hskip
0.7113pth_{n}\hskip 0.7113pt,$ $\displaystyle\overrightarrow{v\,}_{n+1}$
$\displaystyle:=\overrightarrow{v\,}_{n}-\frac{1}{6}u\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}\big{)}h_{n}-\frac{2}{3}u\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}^{\hskip
0.7113pt(1)}\big{)}h_{n}-\frac{1}{6}u\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}^{\hskip 0.7113pt(2)}\big{)}h_{n}\hskip
0.7113pt.$

Our first observation is that $\overrightarrow{x\,}_{n}^{\hskip 0.7113pt(1)}$
approximates the solution at the midpoint of $[t_{n},t_{n+1}]$. and thus can
be improved by including the $O(h^{2})$ terms from the Taylor expansion (3.1).
Another improvement would be to replace $\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}^{\hskip 0.7113pt(2)}\big{)}$ with $\nabla
f\big{(}\hskip 0.7113pt\overrightarrow{x\,}_{n+1}\big{)}$ in the final line.
Since $\overrightarrow{x\,}_{n+1}$ does not depend on
$\overrightarrow{x\,}_{n}^{(2)}$, this makes $\overrightarrow{x\,}_{n}^{\hskip
0.7113pt(2)}$ obsolete and thus can be removed. Finally, we substitute for
$\overrightarrow{v\,}_{n}^{\hskip 0.7113pt(1)}$ and
$\overrightarrow{v\,}_{n}^{\hskip 0.7113pt(2)}$ in the formula of
$\overrightarrow{x\,}_{n+1}$ and rearrange the terms. After performing these
improvements, the Runge-Kutta method simplifies to the following:

$\displaystyle\overrightarrow{x\,}_{n+\frac{1}{2}}$
$\displaystyle:=\overrightarrow{x\,}_{n}+\frac{1}{2}\overrightarrow{v\,}_{n}\hskip
0.7113pth_{n}-\frac{1}{8}u\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}\big{)}h_{n}^{2}\hskip 0.7113pt,$
$\displaystyle\overrightarrow{x\,}_{n+1}$
$\displaystyle:=\overrightarrow{x\,}_{n}+\overrightarrow{v\,}_{n}\hskip
0.7113pth_{n}-\frac{1}{6}u\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}\big{)}h_{n}^{2}-\frac{1}{3}u\nabla
f\big{(}\overrightarrow{x\,}_{n+\frac{1}{2}}\big{)}\hskip
0.7113pth_{n}^{2}\hskip 0.7113pt,$ $\displaystyle\overrightarrow{v\,}_{n+1}$
$\displaystyle:=\overrightarrow{v\,}_{n}-\frac{1}{6}u\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}\big{)}h_{n}-\frac{2}{3}u\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n+\frac{1}{2}}\big{)}h_{n}-\frac{1}{6}u\nabla
f\big{(}\hskip 0.7113pt\overrightarrow{x\,}_{n+1}\big{)}h_{n}\hskip 0.7113pt,$

which coincides with the SORT method when $\gamma=\sigma=0$. If friction is
included in (3.16), and we thus write $\frac{d}{dr}v_{r}^{\hskip
0.7113ptn}=-\gamma v_{r}^{\hskip 0.7113ptn}-u\nabla
f\big{(}x_{r}^{n}\big{)}\hskip 0.7113pt$, we may apply a change of variable to
give

$\displaystyle\frac{d}{dr}\Bigg{(}\begin{matrix}\hskip
0.7113ptx_{r}^{n}\\\\[-4.0pt] \hskip
0.7113ptw_{r}^{n}\end{matrix}\Bigg{)}=\Bigg{(}\begin{matrix}e^{-\gamma
r}w_{r}^{\hskip 0.7113ptn}\\\\[-3.0pt] \,-u\hskip 0.7113pte^{\gamma r}\hskip
0.7113pt\nabla f\big{(}x_{r}^{n}\big{)}\end{matrix}\Bigg{)},$ (3.17)

where $w_{r}^{\hskip 0.7113ptn}:=e^{\gamma r}v_{r}^{\hskip 0.7113ptn}\hskip
0.7113pt$. By applying the third order Runge-Kutta method to this ODE and
applying similar arguments as before, we arrive at the SORT method without
noise terms. Finally, by replacing $\nabla f(\,\cdot\,)$ with $\nabla
f(\,\cdot\,)+\sigma\Big{(}\frac{W_{n}-12K_{n}}{h_{n}}\Big{)}$, we can obtain
the desired method.

### 3.4. Derivation of the SOFA method

To apply a splitting method, we express the ODE’s vector field as a sum of two
operators. Due to the structure of the ODE, we write

$\frac{d}{dr}\Bigg{(}\begin{matrix}\,\overline{x}_{r}^{\hskip
0.7113ptn}\\\\[-4.0pt] \,\overline{v}_{r}^{\hskip
0.7113ptn}\end{matrix}\Bigg{)}=\underbrace{\Bigg{(}\begin{matrix}\,\overline{v}_{r}^{\hskip
0.7113ptn}h_{n}\\\\[-4.0pt] \,0\end{matrix}\Bigg{)}}_{=:\hskip
0.7113ptA}+\underbrace{\Bigg{(}\begin{matrix}\,0\\\\[-4.0pt] \,-\gamma\hskip
0.7113pt\overline{v}_{r}^{\hskip 0.7113ptn}\hskip 0.7113pth_{n}-u\nabla
f\big{(}\hskip 0.7113pt\overline{x}_{r}^{\hskip 0.7113ptn}\big{)}\hskip
0.7113pth_{n}+\sigma\big{(}W_{n}-12K_{n}\big{)}\end{matrix}\Bigg{)}}_{=:\hskip
0.7113ptB}.$

In this case, the ODEs governed by $A$ and $B$ admit the following closed-form
solutions:

$\displaystyle\varphi_{t}^{A}\Bigg{(}\begin{matrix}x\\\\[-4.0pt]
v\end{matrix}\Bigg{)}$ $\displaystyle:=\Bigg{(}\begin{matrix}x+v\hskip
0.7113pth_{n}\hskip 0.7113ptt\\\\[-4.0pt] v\end{matrix}\Bigg{)},$
$\displaystyle\varphi_{t}^{B}\Bigg{(}\begin{matrix}x\\\\[-4.0pt]
v\end{matrix}\Bigg{)}$ $\displaystyle:=\Bigg{(}\begin{matrix}\,x\\\\[-4.0pt]
\,e^{-\gamma\hskip 0.7113pth_{n}t}v+\Big{(}-u\nabla f(x)\hskip
0.7113pth_{n}+\sigma\big{(}W_{n}-12K_{n}\big{)}\Big{)}\Big{(}\frac{1-e^{-\gamma
h_{n}t}}{\gamma\hskip 0.7113pth_{n}}\Big{)}\end{matrix}\Bigg{)}.$

It is important to note that these closed-form solutions are also valid when
$t$ is negative. Thus, it is straightforward to apply the fourth order
splitting due to Forest and Ruth [23]:

$\Bigg{(}\begin{matrix}\,\overline{x}_{1}^{\hskip 0.7113ptn}\\\\[-4.0pt]
\,\overline{v}_{1}^{\hskip
0.7113ptn}\end{matrix}\Bigg{)}\approx\varphi_{\frac{1}{2}+\phi}^{B}\circ\varphi_{1+2\phi}^{A}\circ\varphi_{-\phi}^{B}\circ\varphi_{-1-4\phi}^{A}\circ\varphi_{-\phi}^{B}\circ\varphi_{1+2\phi}^{A}\circ\varphi_{\frac{1}{2}+\phi}^{B}\Bigg{(}\begin{matrix}\,\overline{x}_{0}^{\hskip
0.7113ptn}\\\\[-4.0pt] \,\overline{v}_{0}^{\hskip
0.7113ptn}\end{matrix}\Bigg{)},$

where

$\phi:=\frac{-1+\sqrt[3]{2\hskip 0.7113pt}\,}{2\big{(}2-\sqrt[3]{2\hskip
0.7113pt}\,\big{)}}\hskip 0.7113pt.$

By applying this splitting method to the shifted ODE, we obtain the SOFA
method (definition 1.3). In addition to its high order of convergence, the
above splitting method is symplectic when used on Hamiltonian systems (which
the SDE (1.1) becomes with $\gamma=0$). Hence all that is required to show the
SOFA method is quasi-symplectic in the sense of [45] is to check that the
change in phase volume for every step does not depend on $\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n},\overrightarrow{v\,}_{n}\big{)}$:

$\displaystyle\frac{D\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n+1},\overrightarrow{v\,}_{n+1}\big{)}}{D\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n},\overrightarrow{v\,}_{n}\big{)}}$
$\displaystyle=\frac{D\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n+1},\overrightarrow{v\,}_{n+1}\big{)}}{D\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}^{\hskip
0.7113pt(2)},\overrightarrow{v\,}_{n}^{\hskip
0.7113pt(4)}\big{)}}\frac{D\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}^{\hskip
0.7113pt(2)},\overrightarrow{v\,}_{n}^{\hskip
0.7113pt(4)}\big{)}}{D\big{(}\hskip 0.7113pt\overrightarrow{x\,}_{n}^{\hskip
0.7113pt(2)},\overrightarrow{v\,}_{n}^{\hskip
0.7113pt(3)}\big{)}}\frac{D\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}^{\hskip
0.7113pt(2)},\overrightarrow{v\,}_{n}^{\hskip
0.7113pt(3)}\big{)}}{D\big{(}\hskip 0.7113pt\overrightarrow{x\,}_{n}^{\hskip
0.7113pt(2)},\overrightarrow{v\,}_{n}^{\hskip 0.7113pt(2)}\big{)}}$
$\displaystyle\hskip 42.67912pt\frac{D\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}^{\hskip
0.7113pt(2)},\overrightarrow{v\,}_{n}^{\hskip
0.7113pt(2)}\big{)}}{D\big{(}\hskip 0.7113pt\overrightarrow{x\,}_{n}^{\hskip
0.7113pt(1)},\overrightarrow{v\,}_{n}^{\hskip
0.7113pt(1)}\big{)}}\frac{D\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}^{\hskip
0.7113pt(1)},\overrightarrow{v\,}_{n}^{\hskip
0.7113pt(1)}\big{)}}{D\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n},\overrightarrow{v\,}_{n}^{\hskip
0.7113pt(0)}\big{)}}\frac{D\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n},\overrightarrow{v\,}_{n}^{\hskip
0.7113pt(0)}\big{)}}{D\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n},\overrightarrow{v\,}_{n}\big{)}}$
$\displaystyle=1\cdot e^{-\gamma(\frac{1}{2}+\phi)h_{n}}\cdot e^{\gamma\phi
h_{n}}\cdot e^{\gamma\phi h_{n}}\cdot e^{-\gamma(\frac{1}{2}+\phi)h_{n}}\cdot
1$ $\displaystyle=e^{-\gamma h_{n}},$

where $\big{\\{}\big{(}\hskip 0.7113pt\overrightarrow{x\,}_{n}^{\hskip
0.7113pt(i)},\overrightarrow{v\,}_{n}^{\hskip 0.7113pt(i)}\big{)}\big{\\}}$
are given by definition 1.3 and we use the same notation as in [45].
Similarly, any numerical method for ULD obtained by applying a symplectic
splitting method to the shifted ODE will be quasi-symplectic. As well as the
SOFA method, the Strang splitting method and OBABO scheme discussed in Section
4 are quasi-symplectic.

### 3.5. Error analysis of the shifted ODE method

In this section, we present our results concerning the convergence of the
shifted ODE method and discuss the techniques used in our error analysis.
Recall that we use $\\{(\hskip
0.7113pt\widetilde{x}_{n},\widetilde{v}_{n})\\}_{n\hskip 0.7113pt\geq\hskip
0.7113pt0}$ to denote the Markov chain generated by the shifted ODE method (at
times $\\{t_{n}\\}$). Our main theorem is given below

######  Theorem 3.7 (Convergence of the shifted ODE method).

Let $\big{\\{}(\widetilde{x}_{n}\hskip
0.7113pt,\widetilde{v}_{n})\big{\\}}_{n\hskip 0.7113pt\geq\hskip 0.7113pt0}$
denote the approximation of the underdamped Langevin diffusion
$\big{\\{}(x_{t},v_{t})\big{\\}}_{t\hskip 0.7113pt\geq\hskip 0.7113pt0}$
obtained using the shifted ODE method with a fixed step size $h>0$ and the
same underlying Brownian motion. We assume that the potential $f$ is
$m$-strongly convex and has an $M$-Lipschitz continuous gradient so that the
diffusion process has a unique stationary measure $\pi$. Suppose further that
$(x_{0},v_{0})\sim\pi$ and both processes have the same initial velocity
$\widetilde{v}_{0}=v_{0}\sim\mathcal{N}\big{(}0,u\hskip 0.7113ptI_{d}\big{)}$.
Let $h_{\max}>0$ be fixed. Then for $\lambda\in[0,\frac{\gamma}{2})$ there
exists a constant $c_{1}>0$, independent of $h$ and $d$ such that for $h\leq
h_{\max}$ and $n\geq 0$,

$\displaystyle\|\widetilde{x}_{n}-x_{t_{n}}\|_{\mathbb{L}_{2}}\leq\frac{\sqrt{2\lambda^{2}+2(\gamma-\lambda)^{2}}}{\gamma-2\lambda}\,e^{-n\alpha
h}\big{\|}\widetilde{x}_{0}-x_{0}\big{\|}_{\mathbb{L}_{2}}+c_{1}\sqrt{d\hskip
0.7113pt}\hskip 0.7113pth^{\frac{3}{2}}\hskip 0.7113pt,$ (3.18)

where the contraction rate $\alpha$ is given by

$\displaystyle\alpha$
$\displaystyle=\frac{\big{(}(\gamma-\lambda)^{2}-uM\big{)}\vee\big{(}um-\lambda^{2}\big{)}}{\gamma-2\lambda}\,.$
(3.19)

If we assume further that the Hessian of $f$ is $M_{2}$-Lipschitz continuous,
then there exists a constant $c_{2}>0$, independent of $h$ and $d$ such that
for $h\leq h_{\max}$ and $n\geq 0$,

$\displaystyle\|\widetilde{x}_{n}-x_{t_{n}}\|_{\mathbb{L}_{2}}\leq\frac{\sqrt{2\lambda^{2}+2(\gamma-\lambda)^{2}}}{\gamma-2\lambda}\,e^{-n\alpha
h}\big{\|}\widetilde{x}_{0}-x_{0}\big{\|}_{\mathbb{L}_{2}}+c_{2}\hskip
0.7113ptd\hskip 0.7113pth^{\frac{5}{2}}\hskip 0.7113pt.$ (3.20)

If we additionally assume that $f$ is three times continuously differentiable
and the third derivative is $M_{3}$-Lipschitz continuous, then there exists a
constant $c_{3}>0$, independent of $h$ and $d$ such that for $h\leq h_{\max}$
and $n\geq 0$,

$\displaystyle\|\widetilde{x}_{n}-x_{t_{n}}\|_{\mathbb{L}_{2}}\leq\frac{\sqrt{2\lambda^{2}+2(\gamma-\lambda)^{2}}}{\gamma-2\lambda}\,e^{-\frac{1}{2}n\alpha
h}\big{\|}\widetilde{x}_{0}-x_{0}\big{\|}_{\mathbb{L}_{2}}+c_{3}\hskip
0.7113ptd^{\hskip 0.7113pt\frac{3}{2}}\hskip 0.7113pth^{3}.$ (3.21)

######  Remark 3.8.

These error estimates follow directly using Theorems C.8, D.7 and D.13. In
these theorems and throughout the appendix, the various constants are given
explicitly. As $W_{2}\big{(}\hskip 0.7113pt\widetilde{\nu}_{n}\hskip
0.7113pt,e^{-f}\hskip
0.7113pt\big{)}\leq\|\widetilde{x}_{n}-x_{t_{n}}\|_{\mathbb{L}_{2}}$ where
$\widetilde{x}_{n}\sim\widetilde{\nu}_{n}$, this theorem gives the results in
Table 1.1.

###### Proof.

From the definition of $\\{y_{n}\\}_{n\hskip 0.7113pt\geq\hskip 0.7113pt0}$
(see the end of Section 2 or Theorem B.1), we have

$\displaystyle(\gamma-2\lambda)^{2}\|\hskip
0.7113pt\widetilde{x}_{n}-x_{t_{n}}\|_{\mathbb{L}_{2}}^{2}$
$\displaystyle=\|(\eta-\lambda)(\hskip
0.7113pt\widetilde{x}_{n}-x_{t_{n}})\|_{\mathbb{L}_{2}}^{2}$
$\displaystyle=\big{\|}\big{(}(\eta\hskip
0.7113pt\widetilde{x}_{n}+\widetilde{v}_{n})-(\eta\hskip
0.7113ptx_{t_{n}}+v_{t_{n}})\big{)}-\big{(}(\lambda\hskip
0.7113pt\widetilde{x}_{n}+\widetilde{v}_{n})-(\lambda\hskip
0.7113ptx_{t_{n}}+v_{t_{n}})\big{)}\big{\|}_{\mathbb{L}_{2}}^{2}$
$\displaystyle\leq 2\big{\|}(\lambda\hskip
0.7113pt\widetilde{x}_{n}+\widetilde{v}_{n})-(\lambda\hskip
0.7113ptx_{n}+v_{n})\big{\|}_{\mathbb{L}_{2}}^{2}+2\big{\|}(\eta\hskip
0.7113pt\widetilde{x}_{n}+\widetilde{v}_{n})-(\eta\hskip
0.7113ptx_{n}+v_{n})\big{\|}_{\mathbb{L}_{2}}^{2}$
$\displaystyle=2\|y_{n}\|_{\mathbb{L}_{2}}^{2}\hskip 0.7113pt,$

where we obtained the third line using Minkoski’s inequality. Since we can
estimate each $\|\hskip
0.7113pt\widetilde{x}_{n}-x_{t_{n}}\|_{\mathbb{L}_{2}}$ using
$\|y_{n}\|_{\mathbb{L}_{2}}$, the inequalities (3.18), (3.20) and (3.21) now
follow immediately from Theorems C.8, D.7 and D.13 respectively. ∎

In this paper, we are primarily interested in the dependence of the
$\mathbb{L}_{2}$ error on both the dimension $d$ and step size $h$. However,
under the standard assumptions on the potential ($f$ is $m$-strongly convex
and $\nabla f$ is $M$-Lipschitz continuous), one is additionally interested in
the dependence on the conditional number $\kappa:=\frac{M}{m}$ (see, for
example, Table 1 in [52]). Using the explicit formulae derived in the
appendices, it follows from Corollary C.9 that

$\displaystyle c_{1}$ $\displaystyle=\frac{2\hskip 0.7113pt\sigma
uM\Big{(}\frac{\sqrt{2}}{3}+\frac{\sqrt{10}+\sqrt{15}}{15}\Big{)}}{\big{(}(\gamma-\lambda)^{2}-uM\big{)}\vee\big{(}um-\lambda^{2}\big{)}}+O\big{(}h_{\max}\big{)}.$
(3.22)

Thus the leading term in the error estimate depends linearly on the condition
number $\kappa$, In particular, this means the shifted ODE method matches the
$\mathcal{O}\big{(}\kappa\sqrt{d}\hskip 0.7113pth^{1.5}\big{)}$ convergence
rate of the randomized midpoint method [52] (which was shown to be order
optimal in [8]).

For the rest of the section, we shall outline our approach to deriving the
estimates in Theorem 3.7. Our full error analysis of ULD (with calculations)
is given in the appendix.

Step 1. (Change of variable) For our analysis, we rewrite the shifted ODE
method:

$\displaystyle\Bigg{(}\begin{matrix}\,\widetilde{x}_{n+1}\\\\[-3.0pt]
\,\widetilde{v}_{n+1}\end{matrix}\Bigg{)}:=\Bigg{(}\begin{matrix}\,\widehat{x}_{t_{n+1}}^{\hskip
0.7113ptn}\\\\[-3.0pt] \,\widehat{v}_{t_{n+1}}^{\hskip
0.7113ptn}\end{matrix}\Bigg{)}+12K_{n}\Bigg{(}\begin{matrix}\,0\\\\[-3.0pt]
\,\sigma\end{matrix}\,\Bigg{)},$

where $\big{\\{}\big{(}\widehat{x}_{t}^{\hskip
0.7113ptn},\widehat{v}_{t}^{\hskip 0.7113ptn}\big{)}\big{\\}}_{t\hskip
0.7113pt\in\hskip 0.7113pt[t_{n},t_{n+1}]}$ solves the following ODE,

$\displaystyle\frac{d}{dt}\Bigg{(}\begin{matrix}\,\widehat{x}^{\hskip
0.7113ptn}\\\\[-3.0pt] \,\widehat{v}^{\hskip
0.7113ptn}\end{matrix}\Bigg{)}=\Bigg{(}\begin{matrix}\,\widehat{v}^{\hskip
0.7113ptn}+\sigma\big{(}H_{n}+6K_{n}\big{)}\\\\[-3.0pt] \,-\gamma\big{(}\hskip
0.7113pt\widehat{v}^{\hskip
0.7113ptn}+\sigma\big{(}H_{n}+6K_{n}\big{)}\big{)}-u\nabla
f\big{(}\,\widehat{x}^{\hskip
0.7113ptn}\big{)}\end{matrix}\,\Bigg{)}+\frac{W_{n}-12K_{n}}{h_{n}}\,\Bigg{(}\begin{matrix}\,0\\\\[-3.0pt]
\,\sigma\end{matrix}\,\Bigg{)}\hskip 0.7113pt,$ (3.23)

with initial condition $\big{(}\hskip 0.7113pt\widehat{x}_{t_{n}}^{\hskip
0.7113ptn},\widehat{v}_{t_{n}}^{\hskip
0.7113ptn}\big{)}:=\big{(}\widetilde{x}_{n}\hskip
0.7113pt,\widetilde{v}_{n}\big{)}$. This form of the shifted ODE method is
convenient to work with as we can easily Taylor expand (3.23) from the point
$\big{(}\widetilde{x}_{n}\hskip 0.7113pt,\widetilde{v}_{n}\big{)}$. Moreover,
from the above, we see that the shifted ODE method reduces to another ODE
approximation known as the “log-ODE” method (which is studied in [25]) when
$K_{n}=0$.

Step 2. (2-Wasserstein contractivity) Typically local errors propagate along
the trajectory of the SDE solution so that we would not expect
$\underset{n\hskip 0.7113pt\geq\hskip 0.7113pt0}{\sup}\hskip
0.7113pt\|\widetilde{x}_{n}-x_{t_{n}}\|_{\mathbb{L}_{2}}$ to exist. However
such error estimates (that hold for all $n\geq 0$) can be established when the
SDE exhibits 2-Wasserstein contractivity [39]. For ULD, contractivity does not
hold directly for the solution process $\\{(x_{t},v_{t})\\}_{t\hskip
0.7113pt\geq\hskip 0.7113pt0}$ but instead holds after introducing new
coordinates (this was first observed in [22]). We use the recent approach
given by [16], which considers

$\displaystyle\Bigg{\\{}\Bigg{(}\begin{matrix}\,\lambda\hskip
0.7113ptx_{t}+v_{t}\\\\[-3.0pt] \,\eta\hskip
0.7113ptx_{t}+v_{t}\end{matrix}\Bigg{)}\Bigg{\\}}_{t\hskip 0.7113pt\geq\hskip
0.7113pt0},$ (3.24)

where $\lambda\in[0,\frac{1}{2}\gamma)$ and $\eta:=\gamma-\lambda\hskip
0.7113pt$, to obtain $2$-Wasserstein contractivity. However, unlike in
previous works, we utilized the $2$-Wasserstein contractivity of the
approximation process. That is, using a synchronous coupling and essentially
the same proof as in [16], we establish a $2$-Wasserstein contraction property
for the shifted ODE approximation (Theorem B.1). Alternatively we could have
used the straightforward coupling approach proposed by [12], however this
would have restricted our choice of parameters for ULD to $\gamma=2$ and
$u=\frac{1}{M}$.

We then follow the arguments of [39] with the slight modification that the
contraction is applied to the approximation process instead of the diffusion.
This is illustrated below:

Figure 3.2: Graph outlining the general strategy for our error analysis.

By introducing a third process (obtained by applying the shifted ODE method to
the diffusion at each time $t_{n}$), we can estimate the error between the
coordinate changed approximation and diffusion processes at time $t_{n+1}$ in
terms of the same error at time $t_{n}$.

Step 3. (Local error estimates without global bounds) From figure 3.2, it is
clear that we need to derive error estimates between
$(x_{t_{n+1}},v_{t_{n+1}})$ and $(x_{t_{n+1}}^{\prime},v_{t_{n+1}}^{\prime})$.
Using the Lipschitz continuity of $\nabla f$ along with the fact that both
processes satisfy similar differential equations, it will be straightforward
to derive $\mathbb{L}_{p}$ estimates that are $O\big{(}\sqrt{d\hskip
0.7113pt}\hskip 0.7113pth^{2.5}\big{)}$ (Theorem C.5). Using the
$2$-Wasserstein contractivity of the shifted ODE approximation, this gives a
global $\mathbb{L}_{2}$ error estimate between $\\{(x_{t_{n}},v_{t_{n}})\\}$
and $\\{(\widetilde{x}_{n},\widetilde{v}_{n})\\}$ that is
$\mathcal{O}\big{(}\sqrt{d\hskip 0.7113pt}\hskip 0.7113pth^{1.5}\big{)}$.

Step 4. (Global $\mathbb{L}_{p}$ bounds) To derive the higher order error
estimates, we require that the diffusion and approximation processes are
globally bounded (in an $\mathbb{L}_{p}$ sense). Since the underdamped
Langevin diffusion is well studied, it has already been shown that

$\displaystyle\sup_{t\hskip 0.7113pt\geq\hskip
0.7113pt0}\big{\|}v_{t}\big{\|}_{\mathbb{L}_{p}}\leq C(p)\sqrt{d\hskip
0.7113pt}\hskip 14.22636pt\text{and}\hskip 14.22636pt\sup_{t\hskip
0.7113pt\geq\hskip 0.7113pt0}\big{\|}\nabla
f(x_{t})\big{\|}_{\mathbb{L}_{2}}\leq\sqrt{Md\hskip 0.7113pt},$ (3.25)

for $p\geq 1$ where $(x_{0}\hskip 0.7113pt,v_{0})\sim\pi$ and $C(p)$ is a
positive constant depending only on $p$ [15]. However, unlike in [16], these
bounds were not sufficient to derive high order estimates. Therefore the main
result in this part of our error analysis is the following global
$\mathbb{L}_{4}$ bound:

$\displaystyle\sup_{t\hskip 0.7113pt\geq\hskip 0.7113pt0}\big{\|}\nabla
f(x_{t})\big{\|}_{\mathbb{L}_{4}}\leq\sqrt[4]{3}\,\sqrt{Md\hskip
0.7113pt}\hskip 0.7113pt.$ (3.26)

Step 5. (Local error estimates with global bounds) Using these bounds and the
Lipschitz continuity of $\nabla f$ and $\nabla^{2}f$, we can estimate the
local $\mathbb{L}_{2}$ error as $O\big{(}d\hskip 0.7113pth^{3.5}\big{)}$. We
obtain this estimate by Taylor expanding the “$\nabla f$ terms” in both the
SDE and ODE. By the same arguments as before, we arrive at a global
$\mathbb{L}_{2}$ error estimate that is $\mathcal{O}\big{(}d\hskip
0.7113pth^{2.5}\big{)}$. However, as stated in Theorem 3.7, it is possible to
derive a third order $\mathbb{L}_{2}$ error estimate. To do this, we use the
same strategy of proof as [39] and consider the local mean deviation:

$\displaystyle\Bigg{\|}\mathbb{E}_{n}\Bigg{[}\Bigg{(}\begin{matrix}\,x_{n+1}^{\prime}-x_{t_{n+1}}\\\\[-3.0pt]
\,v_{n+1}^{\prime}-v_{t_{n+1}}\end{matrix}\Bigg{)}\Bigg{]}\Bigg{\|}_{\mathbb{L}_{2}}\hskip
0.7113pt.$

By extending our Taylor expansion of the “$\nabla f$ terms” using the third
derivative $\nabla^{3}f$, we can show that that the local mean deviation is
$O\big{(}d^{1.5}h^{4}\big{)}$. Intuitively, this means the
$O\big{(}h^{3.5}\big{)}$ component of the local error is unbiased and thus
propagates like a random walk. Since random walks grow like $O(\sqrt{N}\hskip
0.7113pt)$, this leads to a global $\mathbb{L}_{2}$ error estimate of
$\mathcal{O}\big{(}h^{3}\big{)}$. To apply this argument rigorously, we follow
the proof of Theorem 1 in [39] and obtain the desired
$\mathcal{O}\big{(}d^{1.5}h^{3}\big{)}$ estimate. Due to some technical
details, this error estimate only has half of the $2$-Wasserstein contraction
rate as the $\mathcal{O}\big{(}\sqrt{d\hskip 0.7113pt}\hskip
0.7113pth^{1.5}\big{)}$ and $\mathcal{O}\big{(}d\hskip 0.7113pth^{2.5}\big{)}$
estimates.

##  4 Related derivative-free numerical methods for ULD

In this section, we will survey some of the numerical methods recently
proposed for ULD. Before we discuss methods specifically designed for Langevin
dynamics, we first note a connection between the shifted ODE method and the
“log-ODE” method studied in [25].

######  Definition 4.1 (Log-ODE method for ULD).

When applied to SDE (1.1), each step of the log-ODE method (Definition 3.12 in
[25]) can be written as $\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n+1},\overrightarrow{v\,}_{n+1}\big{)}:=\big{(}Q_{1}^{n},P_{1}^{n}\big{)}$
where $\big{\\{}\big{(}Q_{t}^{n}\hskip
0.7113pt,P_{t}^{n}\big{)}\big{\\}}_{t\in[0,1]}$ solves the following ODE:

$\displaystyle\frac{d}{dt}\Bigg{(}\begin{matrix}\hskip
0.7113ptQ^{n}\\\\[-3.0pt] \hskip
0.7113ptP^{n}\end{matrix}\Bigg{)}=\Bigg{(}\begin{matrix}\,P^{n}\\\\[-3.0pt]
\,-\gamma P^{n}-u\nabla
f\big{(}Q^{n}\big{)}\end{matrix}\,\Bigg{)}h_{n}+\Bigg{(}\begin{matrix}\,0\\\\[-3.0pt]
\,\sigma\end{matrix}\,\Bigg{)}W_{n}+\Bigg{(}\begin{matrix}\,\sigma\hskip
0.7113pt\\\\[-3.0pt] \,-\gamma\sigma\end{matrix}\,\Bigg{)}h_{n}H_{n}\hskip
0.7113pt,$ (4.1)

with initial condition $\big{(}Q_{0}^{n}\hskip
0.7113pt,P_{0}^{n}\big{)}:=\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}\hskip
0.7113pt,\overrightarrow{v\,}_{n}\big{)}$.

Using the change of variables discussed in the previous section, it follows
that the log-ODE (4.1) can be obtained from the shifted ODE (3.23) simply by
setting $K_{n}=0$. Since the log-ODE method does not use the random vectors
$\\{K_{n}\\}_{n\hskip 0.7113pt\geq\hskip 0.7113pt0}$, it cannot match as many
terms in the stochastic Taylor expansion of ULD as the shifted ODE method. So
although we expect the log-ODE method to achieve a second order convergence
rate when applied to ULD, it is likely to be inferior to the higher order
shifted ODE method.

We shall now discuss four different methods for the underdamped Langevin
diffusion.

######  Definition 4.2 (Strang splitting [3]).

We construct a numerical solution
$\big{\\{}\big{(}\overrightarrow{x\,}_{n},\overrightarrow{v\,}_{n}\big{)}\big{\\}}_{n\hskip
0.7113pt\geq\hskip 0.7113pt0}$ for the SDE (1.1) by setting
$\big{(}\overrightarrow{x\,}_{0},\overrightarrow{v\,}_{0}\big{)}:=(x_{0}\hskip
0.7113pt,v_{0})$ and for $n\geq 0$, defining
$\big{(}\overrightarrow{x\,}_{n+1},\overrightarrow{v\,}_{n+1}\big{)}$ as

$\displaystyle\overrightarrow{v\,}_{n}^{\hskip 0.7113pt(1)}$
$\displaystyle:=\overrightarrow{v\,}_{n}-\frac{1}{2}u\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}\big{)}\hskip 0.7113pth_{n}\hskip 0.7113pt,$
$\displaystyle\Bigg{(}\begin{matrix}\,\overrightarrow{x\,}_{n+1}\\\\[-1.0pt]
\,\overrightarrow{v\,}_{n}^{\hskip 0.7113pt(2)}\end{matrix}\Bigg{)}$
$\displaystyle=\Bigg{(}\begin{matrix}\,\overrightarrow{x\,}_{n}\\\\[-1.0pt]
\,e^{-\gamma h_{n}}\overrightarrow{v\,}_{n}^{\hskip
0.7113pt(1)}\end{matrix}\,\Bigg{)}+\Bigg{(}\begin{matrix}\,\Big{(}\frac{1-e^{-\gamma
h_{n}}}{\gamma}\Big{)}\overrightarrow{v\,}_{n}^{\hskip
0.7113pt(1)}\\\\[-1.0pt]
\,0\end{matrix}\,\Bigg{)}+\sigma\Bigg{(}\begin{matrix}\int_{t_{n}}^{t_{n+1}}\int_{t_{n}}^{t}e^{-\gamma(t-s)}dW_{s}\hskip
0.7113ptdt\\\\[-1.0pt]
\int_{t_{n}}^{t_{n+1}}e^{-\gamma(t_{n+1}-t)}dW_{t}\end{matrix}\,\Bigg{)}\hskip
0.7113pt,$ $\displaystyle\overrightarrow{v\,}_{n+1}$
$\displaystyle:=\overrightarrow{v\,}_{n}^{\hskip
0.7113pt(2)}-\frac{1}{2}u\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n+1}\big{)}\hskip 0.7113pth_{n}\hskip 0.7113pt.$

The above splitting method is obtained by expressing the SDE (1.1) as the sum
of two subsystems but unlike in the derivation of the SOFA method, we split
the SDE as follows:

$d\Bigg{(}\begin{matrix}x_{t}\\\\[-3.0pt]
v_{t}\end{matrix}\Bigg{)}=\underbrace{\Bigg{(}\begin{matrix}\,v_{t}\,dt\\\\[-3.0pt]
-\gamma\hskip 0.7113ptv_{t}\,dt+\sigma\,dW_{t}\end{matrix}\Bigg{)}}_{=:\hskip
0.7113ptA}+\underbrace{\Bigg{(}\begin{matrix}\,0\\\\[-3.0pt] -u\nabla
f(x_{t})\,dt\end{matrix}\Bigg{)}}_{=:\hskip 0.7113ptB}.$

The subsystem $A$ corresponds to a physical Brownian motion [26] and thus can
be sampled exactly. The subsystem $B$ can be solved simply by evaluating the
gradient $\nabla f$. The stochastic integrals which appear in definition 4.2
follow a (joint) normal distribution.

######  Theorem 4.3 (Distribution of integrals required to simulate physical
Brownian motion).

Let $0\leq s\leq t$. Then

$\displaystyle\begin{pmatrix}\int_{s}^{t}e^{-\gamma(t-\tau)}dW_{\tau}\\\
\int_{s}^{t}\int_{s}^{\tau}e^{-\gamma(\tau-r)}dW_{r}\,d\tau\end{pmatrix}$
$\displaystyle\sim\mathcal{N}\left(\Bigg{(}\begin{matrix}0\\\\[-3.0pt]
0\end{matrix}\Bigg{)}\,,\begin{pmatrix}\frac{1-e^{-2\gamma h}}{2\gamma}\hskip
0.7113ptI_{d}&\frac{(1-e^{-\gamma h})^{2}}{2\gamma^{2}}\hskip
0.7113ptI_{d}\\\\[3.0pt] \frac{(1-e^{-\gamma h})^{2}}{2\gamma^{2}}\hskip
0.7113ptI_{d}&\frac{4e^{-\gamma h}-e^{-2\gamma h}+2\gamma
h-3}{2\gamma^{3}}\hskip 0.7113ptI_{d}\end{pmatrix}\right),$

where $h:=t-s$.

We refer the reader to Appendix E for additional details on the stochastic
integrals used by numerical methods for ULD. We note that the second integral
is often written as

$\int_{s}^{t}\int_{s}^{\tau}e^{-\gamma(\tau-r)}dW_{r}\,d\tau=\int_{s}^{t}\bigg{(}\frac{1-e^{-\gamma(t-r)}}{\gamma}\bigg{)}\,dW_{r}\hskip
0.7113pt,$

by Fubini’s theorem. We prefer the first form as it is clearer how it comes
from SDE $A$.

By considering the stochastic Taylor expansion of the above Strang splitting,
we would expect the method to achieve an $\mathbb{L}_{2}$ convergence rate of
$\mathcal{O}\big{(}\sqrt{d\hskip 0.7113pt}h\hskip 0.7113pt\big{)}$ when
$\nabla f$ is $M$-Lipschitz continuous and $\mathcal{O}\big{(}d\hskip
0.7113pth^{2}\hskip 0.7113pt\big{)}$ if, in addition, $\nabla^{2}f$ is also
assumed to be Lipschitz continuous. This conjecture is supported by our
numerical experiment, where the Strang splitting exhibited a second order rate
of convergence. Another particularly attractive feature of this numerical
method is that it only requires one additional evaluation of $\nabla f$ per
step (since, just as for our methods, $\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}\big{)}$ can be computed during the previous
step). Whilst this Strang splitting does not appear to be extensively studied
in the literature, we believe that it offers an attractive compromise between
accuracy and computational cost.

As previously stated, the splitting used to derive the SOFA method differs
from the one used to derive the above Strang splitting. In some sense, the
second splitting is more desirable as the associated numerical methods reduce
to the exact solution when $\nabla f=0$. We derived the SOFA method using the
first operator splitting simply because it results in less computations for
the position component. However, this may be negligible in practice since both
splitting approaches require the same number of gradient evaluations per step.
In any case, the choice of method for the shifted ODE is likely to be problem
dependent.

The following numerical scheme can be viewed as a variant of the Strang
splitting. Although we expect it only have a first order $\mathbb{L}_{2}$
convergence rate, it has the additional advantage that it can be Metropolis-
adjusted (in a similar fashion to standard HMC). Moreover, it was shown in
[47] to achieve a second order $2$ -Wasserstein convergence rate.

######  Definition 4.4 ((Unadjusted) OBABO scheme [4, 47, 54]).

Similar to Strang splitting, we construct a numerical solution
$\big{\\{}\big{(}\overrightarrow{x\,}_{n},\overrightarrow{v\,}_{n}\big{)}\big{\\}}_{n\hskip
0.7113pt\geq\hskip 0.7113pt0}$ for (1.1) by setting
$\big{(}\overrightarrow{x\,}_{0},\overrightarrow{v\,}_{0}\big{)}:=(x_{0}\hskip
0.7113pt,v_{0})$ and for $n\geq 0$, defining
$\big{(}\overrightarrow{x\,}_{n+1},\overrightarrow{v\,}_{n+1}\big{)}$ as

$\displaystyle\hskip 14.22636pt\overrightarrow{v\,}_{n}^{\hskip 0.7113pt(0)}$
$\displaystyle:=e^{-\frac{1}{2}\gamma
h_{n}}\overrightarrow{v\,}_{n}+\sigma\int_{t_{n}}^{t_{n}+\frac{1}{2}h_{n}}e^{-\gamma(t_{n}+\frac{1}{2}h_{n}-t)}dW_{t}\hskip
0.7113pt,\hskip 14.93768pt(O)$ $\displaystyle\hskip
14.22636pt\overrightarrow{v\,}_{n}^{\hskip 0.7113pt(1)}$
$\displaystyle:=\overrightarrow{v\,}_{n}^{\hskip
0.7113pt(0)}-\frac{1}{2}u\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}\big{)}\hskip 0.7113pth_{n}\hskip
0.7113pt,\hskip 109.5431pt(B)$ $\displaystyle\hskip
14.22636pt\overrightarrow{x\,}_{n+1}$
$\displaystyle:=\overrightarrow{x\,}_{n}+\overrightarrow{v\,}_{n}^{\hskip
0.7113pt(1)}h_{n}\hskip 0.7113pt,\hskip 151.51091pt(A)$ $\displaystyle\hskip
14.22636pt\overrightarrow{v\,}_{n}^{\hskip 0.7113pt(2)}$
$\displaystyle:=\overrightarrow{v\,}_{n}^{\hskip
0.7113pt(1)}-\frac{1}{2}u\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n+1}\big{)}\hskip 0.7113pth_{n}\hskip
0.7113pt,\hskip 99.58464pt(B)$ $\displaystyle\hskip
14.22636pt\overrightarrow{v\,}_{n+1}$ $\displaystyle:=e^{-\frac{1}{2}\gamma
h_{n}}\overrightarrow{v\,}_{n}^{\hskip
0.7113pt(2)}+\sigma\int_{t_{n}+\frac{1}{2}h_{n}}^{t_{n+1}}e^{-\gamma(t_{n+1}-t)}dW_{t}\hskip
0.7113pt,\hskip 26.31879pt(O)$ $\displaystyle\hskip
14.22636pt\overrightarrow{x\,}_{n+\frac{1}{2}}$
$\displaystyle:=\overrightarrow{x\,}_{n}+\frac{1}{4}\big{(}\overrightarrow{v\,}_{n}^{\hskip
0.7113pt(0)}+\overrightarrow{v\,}_{n}^{\hskip 0.7113pt(1)}\big{)}h_{n}\hskip
0.7113pt.$

Just as before, the OBABO scheme requires one additional evaluation of $\nabla
f$ per step. Its name is derived from the order in which the damping/sampling
(O), acceleration (A) and free transport parts (B) of the dynamics are
applied. For example, there are also BAOAB and ABOBA methods [33] (which we
view as being similar to Strang splitting). We note the scheme is
straightforward to implement as the stochastic integrals used in the (O)
operations are independent (since Brownian motion has independent increments).

We refer the reader to [47, 54] for details on the Metropolis-adjusted OBABO
scheme. In particular, the adjusted OBABO scheme was shown to empirically
outperform the standard HMC method in [54] (though it was only the third best
performing approach). Based on their stochastic Taylor expansions, the Strang
splitting should outperform the unadjusted OBABO scheme for $\mathbb{L}_{2}$
convergence and this is supported by our experiment. We also note that
$\big{\\{}\overrightarrow{x\,}_{n+\frac{1}{2}}\big{\\}}$ resembles a chain
produced by an ABOBA splitting, which has second order $\mathbb{L}_{2}$
convergence [33]. Since the OBABO scheme and Strang splitting have the same
$2$ -Wasserstein convergence rates, it is not clear how they would compare in
practice.

It was shown in [8] that the following method achieves an optimal
$\mathbb{L}_{2}$ convergence rate of $\mathcal{O}\big{(}\sqrt{d\hskip
0.7113pt}\hskip 0.7113pth^{1.5}\big{)}$ under minimal regularity assumptions
(i.e. $\nabla f$ is $M$-Lipschitz continuous). Unlike the previously discussed
methods, here the gradient $\nabla f$ is evaluated at a uniformly sampled
point within each interval $[t_{n},t_{n+1}]$ (referred to as the “randomized
midpoint”).

######  Definition 4.5 (Randomized midpoint method [52, 31]).

Let $\\{\alpha_{n}\\}_{n\hskip 0.7113pt\geq\hskip 0.7113pt0}$ be a sequence of
independent random variables with $\alpha_{n}\sim U[0,1]$. We construct a
numerical solution
$\big{\\{}\big{(}\overrightarrow{x\,}_{n},\overrightarrow{v\,}_{n}\big{)}\big{\\}}_{n\geq
0}$ by setting
$\big{(}\overrightarrow{x\,}_{0},\overrightarrow{v\,}_{0}\big{)}:=(x_{0}\hskip
0.7113pt,v_{0})$ and for $n\geq 0$, defining
$\big{(}\overrightarrow{x\,}_{n+1},\overrightarrow{v\,}_{n+1}\big{)}$ as

$\displaystyle\overrightarrow{x\,}_{n}^{(1)}$
$\displaystyle:=\overrightarrow{x\,}_{n}+\bigg{(}\frac{1-e^{-\gamma\alpha_{n}h_{n}}}{\gamma}\bigg{)}\overrightarrow{v\,}_{n}-u\bigg{(}\frac{e^{-\gamma\alpha_{n}h_{n}}+\gamma\hskip
0.7113pt\alpha_{n}h_{n}-1}{\gamma^{2}}\bigg{)}\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}\big{)}$ $\displaystyle\hskip
31.44029pt+\sigma\int_{t_{n}}^{t_{n}+\alpha_{n}h_{n}}\int_{t_{n}}^{t}e^{-\gamma(t-s)}dW_{s}\hskip
0.7113ptdt,$ $\displaystyle\overrightarrow{x\,}_{n+1}$
$\displaystyle:=\overrightarrow{x\,}_{n}+\bigg{(}\frac{1-e^{-\gamma
h_{n}}}{\gamma}\bigg{)}\overrightarrow{v\,}_{n}-uh_{n}\bigg{(}\frac{1-e^{-\gamma(1-\alpha_{n})h_{n}}}{\gamma}\bigg{)}\nabla
f\big{(}\hskip 0.7113pt\overrightarrow{x\,}_{n}^{(1)}\big{)}$
$\displaystyle\hskip
31.29802pt+\sigma\int_{t_{n}}^{t_{n+1}}\int_{t_{n}}^{t}e^{-\gamma(t-s)}dW_{s}\hskip
0.7113ptdt,$ $\displaystyle\overrightarrow{v\,}_{n+1}$
$\displaystyle:=e^{-\gamma
h_{n}}\overrightarrow{v\,}_{n}-u\big{(}e^{-\gamma(1-\alpha_{n})h_{n}}h_{n}\big{)}\nabla
f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}^{(1)}\big{)}+\sigma\int_{t_{n}}^{t_{n+1}}e^{-\gamma(t_{n+1}-t)}dW_{t}\hskip
0.7113pt.$

Unlike the previously discussed approaches, the randomized midpoint method
does not exploit the smoothness of $f$ and requires two additional evaluations
of $\nabla f$ per step. However, if the gradient $\nabla f$ is only Lipschitz
continuous, the randomized midpoint method should outperform the other
unadjusted Langevin MCMC methods discussed in this paper. In our numerical
experiment, the randomized midpoint method demonstrates $O(h^{1.5})$
convergence, but is outperformed by the SORT, SOFA and Strang splitting
methods since the target log-density comes from a logistic regression and is
thus infinitely differentiable. Generating the stochastic integrals is also
straightforward and is detailed in Appendix E.

The following numerical method for ULD (and the error analysis used to analyse
it) inspired much of our research. Since it does not yet appear to have a
widely used name, we shall refer to this approach as the “left-point method”.
The left-point method only requires one evaluation of $\nabla f$ per step and
achieves an $\mathbb{L}_{2}$ convergence rate of
$\mathcal{O}\big{(}\sqrt{d\hskip 0.7113pt}\hskip 0.7113pth\big{)}$. It is
obtained by solving (1.1) exactly with $\nabla f(x_{t})$ assumed to be
constant on $[t_{n}\hskip 0.7113pt,t_{n+1}]$. Thus, it differs from the
standard Euler -Maruyama method as only the $\nabla f$ term is fixed.

######  Definition 4.6 (Left-point method [12, 16, 11]).

We define a numerical solution
$\big{\\{}\big{(}\overrightarrow{x\,}_{n},\overrightarrow{v\,}_{n}\big{)}\big{\\}}$
for the SDE (1.1) by setting
$\big{(}\overrightarrow{x\,}_{0},\overrightarrow{v\,}_{0}\big{)}:=(x_{0}\hskip
0.7113pt,v_{0})$ and for $n\geq 0$, defining
$\big{(}\overrightarrow{x\,}_{n+1},\overrightarrow{v\,}_{n+1}\big{)}$ as

$\displaystyle\overrightarrow{x\,}_{n+1}$
$\displaystyle:=\overrightarrow{x\,}_{n}+\bigg{(}\frac{1-e^{-\gamma
h_{n}}}{\gamma}\bigg{)}\overrightarrow{v\,}_{n}-\bigg{(}\frac{e^{-\gamma
h_{n}}+\gamma h_{n}-1}{\gamma^{2}}\bigg{)}u\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}\big{)}+\sigma\int_{t_{n}}^{t_{n+1}}\int_{t_{n}}^{t}e^{-\gamma(t-s)}dW_{s}\hskip
0.7113ptdt,$ $\displaystyle\overrightarrow{v\,}_{n+1}$
$\displaystyle:=e^{-\gamma
h_{n}}\overrightarrow{v\,}_{n}-\bigg{(}\frac{1-e^{-\gamma
h_{n}}}{\gamma}\bigg{)}u\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}\big{)}+\sigma\int_{t_{n}}^{t_{n+1}}e^{-\gamma(t_{n+1}-t)}dW_{t}\hskip
0.7113pt.$

##  5 Numerical experiment

In this section, we shall empirically compare the proposed SOFA method with
the other unadjusted Langevin MCMC algorithms previously discussed. In our
experiment, the target density $\pi(\theta)\propto\exp(-f(\theta))$ comes from
a logistic regression on a real-world dataset (German credit data in [40])
where there are $m=1000$ individuals, each with $d=49$ features
$x_{i}\in\mathbb{R}^{d}$ and a label $y_{i}\in\\{-1,1\\}$ indicating whether
they are creditworthy or not. The model states that
$\mathbb{P}(Y_{i}=y_{i}\hskip 0.7113pt|\hskip
0.7113ptx_{i})=\big{(}1+e^{-y_{i}\hskip
0.7113ptx_{i}^{\mathsf{T}}\theta}\hskip 0.7113pt\big{)}^{-1}$ and so the
function $f$ is defined as

$f(\theta)=\frac{\delta}{2}\|\theta\|_{2}^{2}+\sum_{i=1}^{m}\log\big{(}1+\exp(-y_{i}\hskip
0.7113ptx_{i}^{\mathsf{T}}\theta)\big{)}\hskip 0.7113pt,$

where $\delta$ is a regularization parameter that we set to $\delta=0.1$. We
also set $\gamma=2$ and $u=1$ so that the SDE governing the underdamped
Langevin diffusion (equation (1.1)) becomes

$\displaystyle d\theta_{t}=v_{t}\,dt,\hskip 21.33955ptdv_{t}$
$\displaystyle=-2\hskip 0.7113ptv_{t}\,dt-\nabla
f(\theta_{t})\,dt+2\,dW_{t}\hskip 0.7113pt.$ (5.1)

In the regression, the weight vector $\theta$ takes values in
$\mathbb{R}^{49}$. For example, this means that $W$ is a $d$-dimensional
Brownian motion with $d=49$. We generate $\theta_{0}$ using a normal prior:

$\theta_{0}\sim\mathcal{N}\Big{(}0,10I_{d}\Big{)}\hskip 0.7113pt.$

We use the following Monte Carlo estimator as a proxy for the $\mathbb{L}_{2}$
approximation error.

######  Definition 5.1 (Strong error estimator).

Let $T>0$ be a fixed time and
$\big{\\{}\overrightarrow{\theta\,}_{k}^{h\begin{matrix}\\\\[-12.0pt]
\end{matrix}}\big{\\}}_{k\hskip 0.7113pt\geq\hskip 0.7113pt0}$ denote a
numerical solution of (5.1) computed at times $t_{k}:=kh$ using step size
$h=\frac{T}{N}$ where $N\geq 1$. Similarly, let
$\big{\\{}\overrightarrow{\theta\,}_{k}^{\frac{1}{2}h}\big{\\}}_{k\hskip
0.7113pt\geq\hskip 0.7113pt0}$ be the approximation obtained using step size
$\frac{1}{2}h$. We generate $n$ samples of these numerical solutions (at time
$T$) and define the estimator

$S_{N,n}:=\sqrt{\,\frac{1}{n}\sum_{i=1}^{n}\big{\|}\overrightarrow{\theta\,}_{N,i}^{h\begin{matrix}\\\\[-12.0pt]
\end{matrix}}-\overrightarrow{\theta\,}_{2N,i}^{\frac{1}{2}h}\big{\|}_{2}^{2}\,}\,,$
(5.2)

where each pair
$\big{(}\overrightarrow{\theta\,}_{N,i}^{h\begin{matrix}\\\\[-12.0pt]
\end{matrix}}\hskip
0.7113pt,\overrightarrow{\theta\,}_{2N,i}^{\frac{1}{2}h}\hskip
0.7113pt\big{)}$ is computed from the same sample path of Brownian motion and
each initial condition
$\overrightarrow{\theta\,}_{0,i}^{h\begin{matrix}\\\\[-12.0pt]
\end{matrix}}=\overrightarrow{\theta\,}_{0,i}^{\frac{1}{2}h}$ is sampled from
the normal distribution $\mathcal{N}\big{(}0,10\hskip 0.7113ptI_{d}\big{)}$.

######  Remark 5.2.

By the law of large numbers, the estimator (5.2) converges as
$n\rightarrow\infty$ to

$S_{N}:=\big{\|}\overrightarrow{\theta\,}_{N}^{h\begin{matrix}\\\\[-12.0pt]
\end{matrix}}-\overrightarrow{\theta\,}_{2N}^{\frac{1}{2}h}\big{\|}_{\mathbb{L}_{2}}\,,$

almost surely. For large step sizes, $S_{N}$ may be not close to the true
$\mathbb{L}_{2}$ error
$\big{\|}\overrightarrow{\theta\,}_{N}^{h\begin{matrix}\\\\[-12.0pt]
\end{matrix}}-\theta_{\hskip 0.7113ptT}\big{\|}_{\mathbb{L}_{2}}$.

For the experiment, we use a time horizon of $T=1000$, which is greater than
the mixing time of Hamiltonian Monte Carlo for this specific problem (see
example 3.2 in [1]). We compute the error estimator (5.2) using $n=100$
independent samples of
$\big{(}\overrightarrow{\theta\,}_{N}^{h\begin{matrix}\\\\[-12.0pt]
\end{matrix}}\hskip
0.7113pt,\overrightarrow{\theta\,}_{2N}^{\frac{1}{2}h}\hskip 0.7113pt\big{)}$.

We will now present our results for the numerical experiment that is described
above. Code for this experiment can be found at github.com/james-m-
foster/high-order-langevin.

Figure 5.1: Graph showing $S_{N,n}$ computed for various numerical methods and
step sizes $h=\frac{T}{N}\hskip 0.7113pt$.

From the above graph, we see that the SORT and SOFA methods are the most
accurate when a sufficiently small step size is used. For example, with a step
size of $h=0.005$, the “shifted ODE” methods are roughly 50 times more
accurate than the Strang splitting. Apart from the SOFA method, all of the
numerical methods exhibit the convergence rates that one would expect (see
Section 4 for details). In particular, we note that the OBABO scheme can
achieve second order $2$-Wasserstein convergence to the target distribution
$\pi$. However, since the estimator (5.2) corresponds to the $\mathbb{L}_{2}$
convergence to the SDE solution, the OBABO scheme only achieves first order
convergence in our numerical experiment. Since it is difficult to estimate the
$2$-Wasserstein distance between approximations and solutions of high-
dimensional underdamped Langevin dynamics (see example 4.3 in [19]), we were
unable to empirically show the second order convergence of the OBABO scheme.

Perhaps the most surprising aspect of Figure 5.1 is that the SOFA method
appears to exhibit a fourth order convergence rate (whereas the theory only
suggests a third order). One possible explanation for this is that the
constant in front of the “$\hskip 0.35565ptO(h^{3})$ component” of $S_{N}$ is
negligible for the range of steps sizes used in the numerical experiment. That
is, the noise term in the SDE (5.1) has a small impact on the accuracy of the
approximation (recall that the SOFA method will reduce to a fourth order
splitting method when $\sigma=0$).

Overall the shifted ODE and Strang splitting methods exhibit the fastest
convergence. Furthermore, the SORT and SOFA methods with a step size of
$h=0.01$ are comparable to the Strang splitting with smaller step size of
$h=0.025$ (and are thus more efficient). For MCMC-type applications, the OBABO
scheme may also perform well since it can be Metropolis-adjusted and achieves
a second order $2$-Wasserstein convergence rate when $\nabla f$, $\nabla^{2}f$
are Lipschitz continuous. If the target log-density $f$ only has Lipschitz
gradient, then the randomized midpoint method has the optimal order of strong
convergence [8]. However, for applications where approximate Langevin MCMC is
suitable and the target log-density $f$ has some degree of smoothness, the
theory suggests that the shifted ODE is state-of-the-art and this is supported
by our experiment for the SORT and SOFA methods.

##  6 Conclusion

In this paper, we present an ODE approximation for underdamped Langevin
dynamics and establish non-asymptotic error estimates (in $L^{2}(\mathbb{P})$)
under different assumptions on $f$ which then give bounds on the
$2$-Wasserstein convergence to the stationary distribution. Most notably, we
show that the “shifted ODE” method exhibits third order convergence when
$\nabla f,\nabla^{2}f$ and $\nabla^{3}f$ are all Lipschitz continuous (that
is, it achieves $W_{2}\big{(}\hskip 0.7113pt\widetilde{x}_{n}\hskip
0.7113pt,e^{-f}\hskip 0.7113pt\big{)}\leq\varepsilon$ in
$\mathcal{O}\big{(}\hskip 0.7113pt\sqrt{d\hskip
0.7113pt}/\varepsilon^{\frac{1}{3}}\hskip 0.7113pt\big{)}$ steps). To the best
of our knowledge, this is the first approximation of ULD to achieve
$1/\varepsilon^{\frac{1}{3}}$ dependence without requiring one to evaluate the
higher derivatives of $f$. By discretizing the shifted ODE using standard high
order ODE solvers, we derive practical MCMC algorithms that exhibited third
order convergence in our numerical experiment. In addition to the proposed
methods, we found that the Strang splitting (definition 4.2) is a conceptually
simple method that does not appear to be extensively studied in the literature
yet offers an attractive compromise between accuracy and computational cost.

##  7 Future work

We believe the SORT method requires further analysis and can be potentially
improved. For example, by modifying the noise terms, we arrive at the
following numerical method:

######  Definition 7.1 (A SORT-like method with midpoint computation).

Let $\\{t_{n}\\}_{n\hskip 0.7113pt\geq\hskip 0.7113pt0}$ be a sequence of
times with $t_{0}=0$, $t_{n+1}>t_{n}$ and step sizes
$h_{n}=t_{n+1}-t_{n}\hskip 0.7113pt$, for $n\geq 0$. We construct a numerical
solution $\\{(\overrightarrow{x\,}_{n}\hskip
0.7113pt,\overrightarrow{x\,}_{n+\frac{1}{2}}\hskip
0.7113pt,\overrightarrow{v\,}_{n})\\}_{n\hskip 0.7113pt\geq\hskip 0.7113pt0}$
for the SDE (1.1) by setting $(\overrightarrow{x\,}_{0}\hskip
0.7113pt,\overrightarrow{v\,}_{0}):=(x_{0}\hskip 0.7113pt,v_{0})$ where
$v_{0}\sim\mathcal{N}\big{(}0,u\hskip 0.7113ptI_{d}\big{)}$ and for $n\geq 0$,
defining $(\overrightarrow{x\,}_{n+1}\hskip
0.7113pt,\overrightarrow{x\,}_{n+\frac{1}{2}}\hskip
0.7113pt,\overrightarrow{v\,}_{n+1})$ by

$\displaystyle\overrightarrow{x\,}_{n}^{\hskip 0.7113pt(1)}$
$\displaystyle:=\overrightarrow{x\,}_{n}+\bigg{(}\frac{1-e^{-\frac{1}{2}\gamma
h_{n}}}{\gamma}\bigg{)}\overrightarrow{v\,}_{n}-\bigg{(}\frac{e^{-\frac{1}{2}\gamma
h_{n}}+\frac{1}{2}\gamma h_{n}-1}{\gamma^{2}}\bigg{)}u\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}\big{)}$ $\displaystyle\hskip
31.29802pt+\sqrt{2\gamma u}\hskip 0.7113pt\left(\hskip
0.7113pt\frac{3}{2}\int_{t_{n}}^{t_{n+1}}\bigg{(}\frac{e^{-\gamma(t_{n+1}-t)}+\gamma(t_{n+1}-t)-1}{\gamma^{2}h_{n}}\bigg{)}dW_{t}\right.$
$\displaystyle\hskip
99.58464pt-\left.\frac{1}{4}\int_{t_{n}}^{t_{n+1}}\bigg{(}\frac{1-e^{-\gamma(t_{n+1}-t)}}{\gamma}\bigg{)}\,dW_{t}\hskip
0.7113pt\right),$ $\displaystyle\overrightarrow{x\,}_{n+\frac{1}{2}}$
$\displaystyle:=\overrightarrow{x\,}_{n}+\bigg{(}\frac{1-e^{-\frac{1}{2}\gamma
h_{n}}}{\gamma}\bigg{)}\overrightarrow{v\,}_{n}-\bigg{(}\frac{e^{-\frac{1}{2}\gamma
h_{n}}+\frac{1}{2}\gamma
h_{n}-1}{\gamma^{2}}\bigg{)}\bigg{(}\frac{2}{3}u\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}\big{)}+\frac{1}{3}u\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}^{\hskip 0.7113pt(1)}\big{)}\bigg{)}$
$\displaystyle\hskip 31.29802pt+\sqrt{2\gamma
u}\int_{t_{n}}^{t_{n}+\frac{1}{2}h_{n}}\bigg{(}\frac{1-e^{-\gamma(t_{n}+\frac{1}{2}h_{n}-t)}}{\gamma}\bigg{)}\,dW_{t}\hskip
0.7113pt,$ $\displaystyle\overrightarrow{x\,}_{n+1}$
$\displaystyle:=\overrightarrow{x\,}_{n}+\bigg{(}\frac{1-e^{-\gamma
h_{n}}}{\gamma}\bigg{)}\overrightarrow{v\,}_{n}-\bigg{(}\frac{e^{-\gamma
h_{n}}+\gamma h_{n}-1}{\gamma^{2}}\bigg{)}\bigg{(}\frac{1}{3}\hskip
0.7113ptu\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}\big{)}+\frac{2}{3}\hskip 0.7113ptu\nabla
f\big{(}\hskip 0.7113pt\overrightarrow{x\,}_{n}^{\hskip
0.7113pt(1)}\big{)}\bigg{)}$ $\displaystyle\hskip 31.29802pt+\sqrt{2\gamma
u}\int_{t_{n}}^{t_{n+1}}\bigg{(}\frac{1-e^{-\gamma(t_{n+1}-t)}}{\gamma}\bigg{)}\,dW_{t}\hskip
0.7113pt,$ $\displaystyle\overrightarrow{v\,}_{n+1}$
$\displaystyle:=e^{-\gamma h_{n}}\overrightarrow{v\,}_{n}-\frac{1}{6}\hskip
0.7113pte^{-\gamma h_{n}}u\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}\big{)}h_{n}-\frac{2}{3}\hskip
0.7113pte^{-\frac{1}{2}\gamma h_{n}}u\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n}^{\hskip
0.7113pt(1)}\big{)}h_{n}-\frac{1}{6}u\nabla f\big{(}\hskip
0.7113pt\overrightarrow{x\,}_{n+1}\big{)}h_{n}$ $\displaystyle\hskip
31.29802pt+\sqrt{2\gamma
u}\int_{t_{n}}^{t_{n+1}}e^{-\gamma(t_{n+1}-t)}dW_{t}\hskip 0.7113pt.$

Here the stochastic integrals are jointly normal and thus straightforward to
generate. Similar to the SORT method, we conjecture that this method achieves
a $2$-Wasserstein error of $\varepsilon$ in $\mathcal{O}\big{(}\sqrt{d\hskip
0.7113pt}/\varepsilon^{\frac{1}{k}}\hskip 0.7113pt\big{)}$ steps when the
first $k$ derivatives of $f$ are Lipschitz continuous ($k=1,2,3$). However,
the key advantage of the above unadjsuted Langevin MCMC algorithm is that it
produces $2N$ samples $\\{\overrightarrow{x\,}_{\frac{1}{2}n}\\}_{1\hskip
0.7113pt\leq\hskip 0.7113ptn\hskip 0.7113pt\leq\hskip 0.7113pt2N}$ using only
$2N$ evaluations of $\nabla f$.

Another avenue for future research would be to explore if the “shifted ODE”
approach can be applied to related SDEs, such as high-order [48] or adaptive
[37] Langevin dynamics.

Acknowledgements All the authors were supported by the DataSig programme under
ESPRC grant EP/S026347/1 and by the Alan Turing Institute under EPSRC grant
EP/N510129/1. The first author would like to thank Mufan (Bill) Li and Niloy
Biswas for several interesting discussions on Langevin dynamics and the
coupling of Markov chains. We are also grateful to Sam Power for his valuable
comments regarding ODE integrators.

J. Foster, T. Lyons and H. Oberhauser
Mathematical Institute, University of Oxford, Oxford, OX2 6GG, UK.
The Alan Turing Institute, British Library, London, NW1 2DB, UK.
Email: {james.foster, terry.lyons, harald.oberhauser}@ maths.ox.ac.uk.

## References

  * [1] N. Biswas, P. E. Jacob, and P. Vanetti.  Estimating Convergence of Markov chains with $L$-Lag Couplings.  Advances in Neural Information Processing Systems, 2019\.
  * [2] N. Bou-Rabee and J. M. Sanz-Serna.  Randomized Hamiltonian Monte Carlo.  Annals of Applied Probability, 27(4):2159–2194, 2017\.
  * [3] E. Buckwar, M. Tamborrino, and I. Tubikanec.  Spectral density-based and measure-preserving ABC for partially observed diffusion processes. An illustration on Hamiltonian SDEs.  Statistics and Computing, 30(3):627–648, 2020.
  * [4] G. Bussi and M. Parrinello.  Accurate sampling using Langevin dynamics.  Physical Review E, 75(5):056707, 2007.
  * [5] J. C. Butcher.  Numerical Methods for Ordinary Differential Equations.  John Wiley and Sons, 2016.
  * [6] M. P. Calvo, Sanz-Alonso D., and Sanz-Serna J. M.  HMC: avoiding rejections by not using leapfrog and some results on the acceptance rate.  https://arxiv.org/pdf/1912.03253.pdf, 2020\.
  * [7] C. M. Campos and Sanz-Serna J. M.  Palindromic 3-stage splitting integrators, a roadmap.  Journal of Computational Physics, 346:340–355, 2017\.
  * [8] Y. Cao, J. Lu, and L. Wang.  Complexity of randomized algorithms for underdamped Langevin dynamics.  https://arxiv.org/abs/2003.09906, 2020\.
  * [9] F. Castell and J. G. Gaines.  The ordinary differential equation approach to asymptotically efficient schemes for solution of stochastic differential equations.  Annales de l’Institut Henri Poincaré, 32(2):231–250, 1996\.
  * [10] T. Chen, E. Fox, and C. Guestrin.  Stochastic gradient Hamiltonian Monte Carlo.  In International Conference on Machine Learning, pages 1683–1691, 2014.
  * [11] X. Cheng, N. S. Chatterji, Y. Abbasi-Yadkori, P. L. Bartlett, and M. I. Jordan.  Sharp convergence rates for Langevin dynamics in the nonconvex setting.  https://arxiv.org/abs/1805.01648, 2018\.
  * [12] X. Cheng, N. S. Chatterji, P. L. Bartlett, and M. I. Jordan.  Underdamped Langevin MCMC: A non-asymptotic analysis.  Proceedings of the 31st Conference On Learning Theory, Volume 75 of Proceedings of Machine Learning Research, 2018.
  * [13] J. M. C. Clark and R. J. Cameron.  The maximum rate of convergence of discrete approximations for stochastic differential equations.  in Stochastic Differential Systems Filtering and Control, ed. by Grigelionis (Springer Berlin), 1980.
  * [14] W. T. Coffey, Y. P. Kalmykov, and J. T. Waldron.  The Langevin Equation: With Applications to Stochastic Problems in Physics, Chemistry and Electrical Engineering.  World Scientifc, 2012.
  * [15] A. S. Dalalyan.  Further and stronger analogy between sampling and optimization: Langevin monte carlo and gradient descent.  Proceedings of the 2017 Conference on Learning Theory, volume 65 of Proceedings of Machine Learning Research:678–689, 2017.
  * [16] A. S. Dalalyan and L. Riou-Durand.  On sampling from a log-concave density using kinetic Langevin diffusions.  Bernoulli, 26(3):1956–1988, 2020.
  * [17] A. S. Dickinson.  Optimal Approximation of the Second Iterated Integral of Brownian Motion.  Stochastic Analysis and Applications, 25(5):1109–1128, 2007\.
  * [18] Z. Ding and Q. Li.  Random Coordinate Underdamped Langevin Monte Carlo.  https://arxiv.org/abs/2010.11366, 2020\.
  * [19] M. Dobson, Y. Li, and J. Zhai.  Using coupling methods to estimate sample quality for stochastic differential equations.  https://arxiv.org/abs/1912.10339, 2019\.
  * [20] S. S. Dragomir.  Some Gronwall Type Inequalities and Applications.  Nova Science Publishers, 2003.
  * [21] S. Duane, A. D. Kennedy, B. J. Pendleton, and D. Roweth.  Hybrid Monte Carlo.  Physics Letters B, 195(2), 1987.
  * [22] A. Eberle, A. Guillin, and R. Zimmer.  Couplings and quantitative contraction rates for Langevin dynamics.  Annals of Probability, 47:1982–2010, 2019.
  * [23] E. Forest and R. D. Ruth.  Fourth-order symplectic integration.  Physica D: Nonlinear Phenomena, 43(1):105–117, 1990\.
  * [24] J. Foster.  Numerical approximations for stochastic differential equations.  PhD thesis, University of Oxford, 2020.
  * [25] J. Foster, T. Lyons, and H. Oberhauser.  An optimal polynomial approximation of Brownian motion.  SIAM Journal on Numerical Analysis, 58(3):1393–1421, 2020\.
  * [26] P. Fritz, Gassiat. P, and T. Lyons.  Physical Brownian motion in magnetic field as a rough path.  Transaction or the American Mathematical Society, 367:7939–7955, 2015.
  * [27] J. G. Gaines and T. Lyons.  Random Generation of Stochastic Area Integrals.  SIAM Journal on Applied Mathematics, 54(4):1132–1146, 1994\.
  * [28] X. Gao, M. Gürbüzbalaban, and L. Zhu.  Global convergence of Stochastic Gradient Hamiltonian Monte Carlo for non-convex stochastic optimization: Non-asymptotic performance bounds and momentum-based acceleration.  https://arxiv.org/abs/1809.04618, 2018\.
  * [29] X. Gao, M. Gürbüzbalaban, and L. Zhu.  Breaking Reversibility Accelerates Langevin Dynamics for Non-Convex Optimization.  Advances in Neural Information Processing Systems, 2020\.
  * [30] M. Girolami and B. Calderhead.  Riemann manifold Langevin and Hamiltonian Monte Carlo methods.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73(2):123–214, 2011.
  * [31] Y. He, K. Balasubramanian, and M. A. Erdogdu.  On the Ergodicity, Bias and Asymptotic Normality of Randomized Midpoint Sampling Method.  Advances in Neural Information Processing Systems, 2020\.
  * [32] P. E. Kloeden and E. Platen.  Numerical Solution of Stochastic Differential Equations.  Springer, Berlin, 1992.
  * [33] B. Leimkuhler and Matthews. C.  Rational Construction of Stochastic Numerical Methods for Molecular Sampling.  Applied Mathematics Research eXpress, 2013(1):34–56, 2013\.
  * [34] B. Leimkuhler and C. Matthews.  Molecular Dynamics: With Deterministic and Stochastic Numerical Methods.  Interdisciplinary Applied Mathematics, Springer, 2015.
  * [35] B. Leimkuhler, C. Matthews, and G. Stoltz.  The computation of averages from equilibrium and nonequilibrium Langevin molecular dynamics.  IMA Journal of Numerical Analysis, 36(1):16–79, 2016\.
  * [36] B. Leimkuhler, T. Pouchon, T. Vlaar, and A. Storkey.  Constraint-Based Regularization of Neural Networks.  https://arxiv.org/abs/2006.10114, 2020\.
  * [37] B. Leimkuhler, M. Sachs, and G. Stoltz.  Hypocoercivity Properties of Adaptive Langevin Dynamics.  SIAM Journal on Applied Mathematics, 80(3):1197–1222, 2020\.
  * [38] X. Lelièvre and G. Stoltz.  Partial differential equations and stochastic methods in molecular dynamics.  Acta Numerica, 25:681–880, 2016.
  * [39] X. Li, D. Wu, L. Mackey, and M. A. Erdogdu.  Stochastic Runge-Kutta accelerates Langevin Monte Carlo and beyond.  Advances in Neural Information Processing Systems, 2019\.
  * [40] M. Lichman.  UCI machine learning repository.  https://archive.ics.uci.edu/ml, 2013\.
  * [41] T. Lyons, M. Caruana, and T. Lévy.  Differential equations driven by rough paths.  In École d’été de probabilités de Saint-Flour XXXIV-2004, edited by J. Picard in Volume 1908 of Lecture Notes in Mathematics, Springer, Berlin, 2007.
  * [42] Y. Ma, N. Chatterji, X. Cheng, N. Flammarion, and P. Bartlett.  Is There an Analog of Nesterov Acceleration for MCMC?  https://arxiv.org/abs/1902.00996, 2019\.
  * [43] Y. Ma, T. Chen, and E. Fox.  A complete recipe for stochastic gradient MCMC.  Advances in Neural Information Processing Systems, 2015\.
  * [44] S. Mandt, M. D. Hoffman, and D. M. Blei.  Stochastic gradient descent as approximate Bayesian inference.  Journal of Machine Learning Research, 18(1):4873–4907, 2017\.
  * [45] G. Milstein and M. Tretyakov.  Quasi-symplectic methods for Langevin-type equations.  IMA Journal of Numerical Analysis, 23(4):593–626, 2003\.
  * [46] G. N. Milstein and M. V. Tretyakov.  Stochastic Numerics for Mathematical Physics.  Springer, Berlin, 2004.
  * [47] P. Monmarché.  High-dimensional MCMC with a standard splitting scheme for the underdamped Langevin diffusion.  https://arxiv.org/pdf/2007.05455, 2020\.
  * [48] W. Mou, Y. Ma, M. J. Wainwright, P. L. Bartlett, and M. I. Jordan.  High-Order Langevin Diffusion Yields an Accelerated MCMC Algorithm.  https://arxiv.org/abs/1908.10859, 2019\.
  * [49] R. M. Neal.  MCMC using Hamiltonian dynamics.  in Handbook of Markov Chain Monte Carlo, ed. by Brooks, Gelman, Jones and Meng (Chapman and Hall/CRC Handbooks of Modern Statistical Methods), 2011\.
  * [50] G. A. Pavliotis.  Stochastic Processes and Applications.  Springer, New York, 2014.
  * [51] W. Rudin.  Real and complex analysis.  Third edition, McGraw-Hill Book Co., New York, 1987.
  * [52] R. Shen and Y. T. Lee.  The Randomized Midpoint Method for Log-Concave Sampling.  Advances in Neural Information Processing Systems, 2019\.
  * [53] U. Simsekli, L. Zhu, Y. W. Teh, and M. Gürbüzbalaban.  Fractional Underdamped Langevin Dynamics: Retargeting SGD with Momentum under Heavy-Tailed Gradient Noise.  Proceedings of the 36th International Conference on Machine Learning (ICML), 2020.
  * [54] Z. Song and Z. Tan.  Hamiltonian Assisted Metropolis Sampling.  https://arxiv.org/pdf/2005.08159, 2020\.
  * [55] Z. Wang and H. Delingette.  Quasi-symplectic Langevin Variational Autoencoder.  https://arxiv.org/abs/2009.01675, 2020\.
  * [56] M. Welling and Y. W. Teh.  Bayesian Learning via Stochastic Gradient Langevin Dynamics.  Proceedings of the 28th International Conference on Machine Learning (ICML), 2011.

##  Appendix A Preliminary theorems

We first recall the following standard results, which will be used within our
error analysis.

######  Theorem A.1.

For $p\geq 1$ and $a_{1},\cdots,a_{k}\in\mathbb{R}^{n}$, we have

$\displaystyle\Bigg{\|}\sum_{i=1}^{k}a_{i}\hskip 0.7113pt\Bigg{\|}_{2}^{p}\leq
k^{p-1}\sum_{i=1}^{k}\|a_{i}\|_{2}^{p}\hskip 0.7113pt.$

###### Proof.

For any $p\geq 1$, the function $g(x)=x^{p}$ is convex and strictly increasing
on $[0,\infty)$. Applying the triangle inequality followed by Jensen’s
inequality gives

$\displaystyle g\Bigg{(}\hskip 0.7113pt\frac{1}{k}\hskip
0.7113pt\Bigg{\|}\sum_{i=1}^{k}a_{i}\hskip 0.7113pt\Bigg{\|}_{2}\hskip
0.7113pt\Bigg{)}\leq
g\Bigg{(}\frac{1}{k}\sum_{i=1}^{k}\|a_{i}\|_{2}\Bigg{)}\leq\frac{1}{k}\sum_{i=1}^{k}g\big{(}\|a_{i}\|_{2}\big{)}.$

The result immediately follows from the definition of $g(x)$. ∎

######  Theorem A.2.

Let $z:[s,t]\rightarrow\mathbb{R}^{n}$ be such that $\|z\|_{2}^{p}$ is
integrable on $[s,t]$ for some $p\geq 1$. Then

$\displaystyle\bigg{\|}\int_{s}^{t}z_{r}\,dr\,\bigg{\|}_{2}^{p}\leq(t-s)^{p-1}\int_{s}^{t}\|z_{r}\|_{2}^{p}\,dr.$

###### Proof.

For any $p\geq 1$, the function $g(x)=x^{p}$ is convex and strictly increasing
on $[0,\infty)$. By viewing $\hskip
0.7113pt\frac{1}{t-s}\int_{s}^{t}\,\cdot\,\,d\tau\hskip 0.7113pt$ as an
expectation with respect to the uniform measure on $[s,t]$, we can apply
Jensen’s inequality to give

$\displaystyle\frac{1}{t-s}\int_{s}^{t}g\big{(}\|z_{r}\|_{2}\big{)}\,dr\,\leq\,g\bigg{(}\frac{1}{t-s}\int_{s}^{t}\|z_{r}\|_{2}\,dr\bigg{)}.$

The result immediately follows from the definition of $g(x)$. ∎

######  Theorem A.3 (Minkowski’s inequality).

Let $X$ and $Y$ be $\mathbb{R}^{n}$-valued random variables that are
$\mathbb{L}_{p}$-integrable for some $p\geq 1$. Then

$\displaystyle\big{\|}X+Y\big{\|}_{\mathbb{L}_{p}}\leq\big{\|}X\big{\|}_{\mathbb{L}_{p}}+\big{\|}Y\big{\|}_{\mathbb{L}_{p}}.$
(A.1)

###### Proof.

Setting $z:=X\cdot\mathbbm{1}_{[0,1]}+Y\cdot\mathbbm{1}_{[1,2]}$ and
$[s,t]=[0,2]$ in Theorem A.4 gives the result. ∎

######  Theorem A.4 (Minkowski’s inequality for integrals).

Let $z$ denote a $\mathbb{R}^{n}$-valued stochastic process which is
$\mathbb{L}_{p}$-integrable over an interval $[s,t]$ for some $p\geq 1$. Then

$\displaystyle\bigg{\|}\int_{s}^{t}z_{r}\,dr\,\bigg{\|}_{\mathbb{L}_{p}}\leq\int_{s}^{t}\|z_{r}\|_{\mathbb{L}_{p}}\,dr$
(A.2)

###### Proof.

By Theorem A.2, it is clear that the left hand side of (A.2) is finite.
Moreover, when $p=1$, the result directly follows by the linearity of
expectation. For $p>1$, we have

$\displaystyle\bigg{\|}\int_{s}^{t}z_{r}\,dr\hskip
0.7113pt\bigg{\|}_{\mathbb{L}_{p}}^{p}$
$\displaystyle=\mathbb{E}\Bigg{[}\hskip
0.7113pt\bigg{\|}\int_{s}^{t}z_{\tau}\,d\tau\hskip 0.7113pt\bigg{\|}_{2}\hskip
0.7113pt\bigg{\|}\int_{s}^{t}z_{r}\,dr\hskip 0.7113pt\bigg{\|}_{2}^{p-1}\hskip
0.7113pt\Bigg{]}$ $\displaystyle\leq\mathbb{E}\Bigg{[}\hskip
0.7113pt\int_{s}^{t}\|z_{\tau}\|_{2}\,d\tau\hskip
0.7113pt\bigg{\|}\int_{s}^{t}z_{r}\,dr\hskip 0.7113pt\bigg{\|}_{2}^{p-1}\hskip
0.7113pt\Bigg{]}$ $\displaystyle=\int_{s}^{t}\mathbb{E}\Bigg{[}\hskip
0.7113pt\|z_{\tau}\|_{2}\bigg{\|}\int_{s}^{t}z_{r}\,dr\hskip
0.7113pt\bigg{\|}_{2}^{p-1}\hskip 0.7113pt\Bigg{]}\hskip 0.7113ptd\tau$
$\displaystyle\leq\int_{s}^{t}\mathbb{E}\Big{[}\|z_{\tau}\|_{2}^{p}\Big{]}^{\frac{1}{p}}\mathbb{E}\Bigg{[}\hskip
0.7113pt\bigg{\|}\int_{s}^{t}z_{r}\,dr\hskip
0.7113pt\bigg{\|}_{2}^{(p-1)\cdot\frac{p}{p-1}}\hskip
0.7113pt\Bigg{]}^{1-\frac{1}{p}}\hskip 0.7113ptd\tau$
$\displaystyle=\int_{s}^{t}\|z_{\tau}\|_{\mathbb{L}_{p}}\,d\tau\,\bigg{\|}\int_{s}^{t}z_{r}\,dr\hskip
0.7113pt\bigg{\|}_{\mathbb{L}_{p}}^{1-\frac{1}{p}},$

where the penultimate line comes from Hölder’s inequality. The result now
follows. ∎

######  Remark A.5.

We also use the Minkowski’s inequalities for the $\mathcal{F}_{n}$-conditional
$\mathbb{L}_{p}$ norm:

$\displaystyle\|X\|_{\mathbb{L}_{p}^{n}}:=\mathbb{E}\big{[}\|X\|_{2}^{p}\,\big{|}\hskip
0.7113pt\mathcal{F}_{n}\big{]}^{\frac{1}{p}}.$

######  Theorem A.6.

Let $Z$ be a random variable in $\mathbb{R}^{n}$ and $\mathcal{F}$ be a
$\sigma$-algebra. Then for $p\geq 1$,

$\displaystyle\big{\|}\mathbb{E}[Z|\hskip
0.7113pt\mathcal{F}\,]\big{\|}_{\mathbb{L}_{p}}\leq\big{\|}Z\big{\|}_{\mathbb{L}_{p}}.$

###### Proof.

By Jensen’s inequality and the Tower law, we have

$\displaystyle\big{\|}\mathbb{E}[Z|\hskip
0.7113pt\mathcal{F}\,]\big{\|}_{\mathbb{L}_{p}}=\mathbb{E}\Big{[}\big{\|}\mathbb{E}[Z|\hskip
0.7113pt\mathcal{F}\,]\big{\|}_{2}^{p}\Big{]}^{\frac{1}{p}}\leq\mathbb{E}\Big{[}\hskip
0.7113pt\mathbb{E}\big{[}\|Z\|_{2}^{p}\,\big{|}\hskip
0.7113pt\mathcal{F}\,\big{]}\Big{]}^{\frac{1}{p}}=\mathbb{E}\Big{[}\|Z\|_{2}^{p}\Big{]}^{\frac{1}{p}}=\big{\|}Z\big{\|}_{\mathbb{L}_{p}},$

which gives the desired result. ∎

Just as in [39], we shall apply the following theorem to establish a global
error estimate.

######  Theorem A.7.

Suppose $X,Y,Z\in\mathbb{R}^{n}$ are square-integrable random variables and
that $Z$ is $\mathcal{F}$-measurable (where $\mathcal{F}$ is some
$\sigma$-algebra). Then for constants $c,h>0$, we have

$\displaystyle\|X+Y\|_{\mathbb{L}_{2}}^{2}\hskip 0.7113pt\leq\hskip
0.7113pt\|X\|_{\mathbb{L}_{2}}^{2}+2\hskip
0.7113pt\big{\|}X-\mathbb{E}[X|\hskip
0.7113pt\mathcal{F}\,]\big{\|}_{\mathbb{L}_{2}}\|Y-Z\|_{\mathbb{L}_{2}}+\frac{c}{h}\big{\|}\hskip
0.7113pt\mathbb{E}[X|\hskip
0.7113pt\mathcal{F}\,]\big{\|}_{\mathbb{L}_{2}}^{2}+\Big{(}1+\frac{h}{c}\Big{)}\|Y\|_{\mathbb{L}_{2}}^{2}\hskip
0.7113pt.$

###### Proof.

By the Cauchy-Schwarz and Young’s inequalities, we have

$\displaystyle\Big{|}\Big{\langle}\frac{c}{h}x,y\Big{\rangle}\Big{|}\leq\Big{\|}\frac{c}{h}x\Big{\|}_{2}\|y\|_{2}\leq\frac{1}{2}\Big{\|}\frac{c}{h}x\Big{\|}_{2}^{2}+\frac{1}{2}\|y\|_{2}^{2}\,,$

for vectors $x,y\in\mathbb{R}^{n}$. Rearranging the above produces the below
inequality:

$\displaystyle\big{|}\langle
x,y\rangle\big{|}\leq\frac{c}{2h}\|x\|_{2}^{2}+\frac{h}{2c}\|y\|_{2}^{2}\,.$

Applying this to the random variables $\mathbb{E}[X|\hskip
0.7113pt\mathcal{F}\,]$ and $Y$ gives

$\displaystyle\big{\|}\mathbb{E}\big{[}\langle X,Y\hskip
0.7113pt\rangle\big{|}\hskip
0.7113pt\mathcal{F}\,\big{]}\big{\|}_{\mathbb{L}_{1}}$
$\displaystyle\leq\big{\|}\mathbb{E}\big{[}\langle X-\mathbb{E}[X|\hskip
0.7113pt\mathcal{F}\,],Y\hskip 0.7113pt\rangle|\hskip
0.7113pt\mathcal{F}\,\big{]}\big{\|}_{\mathbb{L}_{1}}+\big{\|}\mathbb{E}\big{[}\langle\hskip
0.7113pt\mathbb{E}[X|\hskip 0.7113pt\mathcal{F}\,],Y\hskip
0.7113pt\rangle\big{]}|\hskip
0.7113pt\mathcal{F}\,\big{]}\big{\|}_{\mathbb{L}_{1}}$
$\displaystyle\leq\big{\|}\mathbb{E}\big{[}\langle X-\mathbb{E}[X|\hskip
0.7113pt\mathcal{F}\,],Y-Z\rangle|\hskip
0.7113pt\mathcal{F}\,\big{]}\big{\|}_{\mathbb{L}_{1}}+\big{\|}\langle\hskip
0.7113pt\mathbb{E}[X|\hskip 0.7113pt\mathcal{F}\,],Y\hskip
0.7113pt\rangle\big{\|}_{\mathbb{L}_{1}}$
$\displaystyle\leq\big{\|}X-\mathbb{E}[X|\hskip
0.7113pt\mathcal{F}\,]\big{\|}_{\mathbb{L}_{2}}\big{\|}Y-Z\big{\|}_{\mathbb{L}_{2}}+\frac{c}{2h}\|\mathbb{E}[X|\hskip
0.7113pt\mathcal{F}\,]\|_{\mathbb{L}_{2}}^{2}+\frac{h}{2c}\|Y\|_{\mathbb{L}_{2}}^{2}\hskip
0.7113pt,$

by the Cauchy-Schwarz inequality. The second line follows by employing the
Tower law (for the first term) and Theorem A.6 (for the second term).
Expanding $\|X+Y\|_{\mathbb{L}_{2}}^{2}$ gives

$\displaystyle\|X+Y\|_{\mathbb{L}_{2}}^{2}$
$\displaystyle=\|X\|_{\mathbb{L}_{2}}^{2}+2\hskip
0.7113pt\mathbb{E}\big{[}\langle X,Y\rangle\big{]}+\|Y\|_{\mathbb{L}_{2}}^{2}$
$\displaystyle=\|X\|_{\mathbb{L}_{2}}^{2}+2\hskip
0.7113pt\mathbb{E}\big{[}\mathbb{E}\big{[}\langle X,Y\rangle\big{|}\hskip
0.7113pt\mathcal{F}\,\big{]}\big{]}+\|Y\|_{\mathbb{L}_{2}}^{2}$
$\displaystyle\leq\|X\|_{\mathbb{L}_{2}}^{2}+2\hskip
0.7113pt\big{\|}\mathbb{E}\big{[}\langle X,Y\hskip
0.7113pt\rangle\big{|}\hskip
0.7113pt\mathcal{F}\,\big{]}\big{\|}_{\mathbb{L}_{1}}+\|Y\|_{\mathbb{L}_{2}}^{2}\hskip
0.7113pt,$

and the result follows. ∎

Along with Brownian motion, we also use Brownian bridge processes in our
analysis.

######  Definition A.8.

Let $W$ be a standard $d$-dimensional Brownian motion. Then for $n\geq 0$, we
consider the Brownian bridge

$\displaystyle
B_{t}=W_{t}-\bigg{(}W_{t_{n}}+\bigg{(}\frac{t-t_{n}}{h_{n}}\bigg{)}\hskip
0.7113ptW_{t_{n},t_{n+1}}\bigg{)},$

and define the “shifted” Brownian bridge

$\displaystyle\widetilde{B}_{t}:=B_{t}+12K_{n}\bigg{(}\frac{t-t_{n}}{h_{n}}\bigg{)},$

for $t\in[t_{n},t_{n+1}]$ where $h_{n}:=t_{n+1}-t_{n}\hskip 0.7113pt$.

We now give some basic $\mathbb{L}_{p}$ estimates for these processes as well
as the random variable $H_{n}+6K_{n}$ which appears in the shifted ODE (and
thus appears throughout our analysis).

######  Theorem A.9.

Let $W,B$ and $\widetilde{B}$ be the processes given by definition A.8. Then
for all $n\geq 0$ and $t\in[t_{n},t_{n+1}]$, we have

$\displaystyle\big{\|}\hskip 0.7113ptB_{t_{n},t}\big{\|}_{\mathbb{L}_{p}}$
$\displaystyle\leq C(p)\sqrt{d\hskip 0.7113pt}\hskip
0.7113pt(t-t_{n})^{\frac{1}{2}},$ $\displaystyle\big{\|}\hskip
0.7113pt\widetilde{B}_{t_{n},t}\big{\|}_{\mathbb{L}_{p}}$ $\displaystyle\leq
C_{B}(p)\sqrt{d\hskip 0.7113pt}\hskip 0.7113pt(t-t_{n})^{\frac{1}{2}},$
$\displaystyle\big{\|}H_{n}+6K_{n}\big{\|}_{\mathbb{L}_{p}}$
$\displaystyle=C_{K}(p)\sqrt{d\hskip 0.7113pt}\hskip
0.7113pt(h_{n})^{\frac{1}{2}},$

for $p\in\\{2,4,8\\}$ where the constants $C_{B}(p)$ and $C_{K}(p)$ are
defined as

$\displaystyle C(p)$ $\displaystyle:=\begin{cases}\hskip
0.7113pt1&\text{if}\,\,\,p=2\\\ \hskip 0.7113pt\big{(}1+2d^{\hskip
0.7113pt-1}\big{)}^{\frac{1}{4}}&\text{if}\,\,\,p=4\\\ \hskip
0.7113pt\big{(}1+12d^{\hskip 0.7113pt-1}+44d^{\hskip 0.7113pt-2}+48d^{\hskip
0.7113pt-3}\big{)}^{\frac{1}{8}}&\text{if}\,\,\,p=8\end{cases}\,,$
$\displaystyle C_{B}(p)$
$\displaystyle:=\left(1+\frac{\sqrt{5}}{5}\right)C(p),$ $\displaystyle
C_{K}(p)$ $\displaystyle:=\frac{\sqrt{30}}{15}C(p).$

(Note that $\hskip 0.7113ptC(p)\sqrt{d\,}$ is the $\mathbb{L}_{p}$ norm of a
standard $d$-dimensional normal random variable).

###### Proof.

To begin, we will consider a standard $d$-dimensional normal random variable
$Z$. The $\mathbb{L}_{2}$ and $\mathbb{L}_{4}$ norms of
$Z=\big{(}Z^{(1)},\cdots,Z^{(d)}\big{)}$ are straightforward to explicitly
compute.

$\displaystyle\big{\|}Z\hskip 0.7113pt\big{\|}_{\mathbb{L}_{2}}^{2}$
$\displaystyle=\mathbb{E}\Bigg{[}\sum_{k\hskip 0.7113pt=\hskip
0.7113pt0}^{d}\big{(}Z^{(k)}\big{)}^{2}\Bigg{]}$ $\displaystyle=d,$
$\displaystyle\big{\|}Z\hskip 0.7113pt\big{\|}_{\mathbb{L}_{4}}^{4}$
$\displaystyle=\mathbb{E}\left[\left(\,\sum_{k\hskip 0.7113pt=\hskip
0.7113pt0}^{d}\big{(}Z^{(k)}\big{)}^{2}\right)^{2}\,\right]$
$\displaystyle=\sum_{\begin{subarray}{c}0\hskip 0.7113pt\leq\hskip
0.7113pti,\hskip 0.7113ptj\hskip 0.7113pt\leq\hskip 0.7113ptd\\\\[1.0pt] i\neq
j\end{subarray}}\mathbb{E}\Big{[}\big{(}Z^{(i)}\big{)}^{2}\big{(}Z^{(j)}\big{)}^{2}\Big{]}+\sum_{k\hskip
0.7113pt=\hskip
0.7113pt0}^{d}\mathbb{E}\Big{[}\big{(}Z^{(k)}\big{)}^{4}\Big{]}$
$\displaystyle=d(d-1)+3d$ $\displaystyle=d^{2}+2d.$

Note that we used the fact that centred normal distributions have trivial odd
moments. Similarly, one can compute the $\mathbb{L}_{8}$ norm of $Z$ (though
the calculation is more involved).

$\displaystyle\big{\|}Z\hskip 0.7113pt\big{\|}_{\mathbb{L}_{8}}^{8}$
$\displaystyle=\mathbb{E}\left[\left(\,\sum_{k\hskip 0.7113pt=\hskip
0.7113pt0}^{d}\big{(}Z^{(k)}\big{)}^{2}\right)^{4}\,\right]$
$\displaystyle=\sum_{\begin{subarray}{c}0\hskip 0.7113pt\leq\hskip
0.7113pti,\hskip 0.7113ptj,\hskip 0.7113ptk,\hskip 0.7113ptl\hskip
0.7113pt\leq\hskip 0.7113ptd\\\\[1.0pt]
|\\{i,j,k,l\\}|=4\end{subarray}}\mathbb{E}\Big{[}\big{(}Z^{(i)}\big{)}^{2}\big{(}Z^{(j)}\big{)}^{2}\big{(}Z^{(k)}\big{)}^{2}\big{(}Z^{(l)}\big{)}^{2}\Big{]}+6\sum_{\begin{subarray}{c}0\hskip
0.7113pt\leq\hskip 0.7113pti,\hskip 0.7113ptj,\hskip 0.7113ptk\hskip
0.7113pt\leq\hskip 0.7113ptd\\\\[1.0pt]
|\\{i,j,k\\}|=3\end{subarray}}\mathbb{E}\Big{[}\big{(}Z^{(i)}\big{)}^{4}\big{(}Z^{(j)}\big{)}^{2}\big{(}Z^{(k)}\big{)}^{2}\Big{]}$
$\displaystyle\hskip 28.45274pt+3\sum_{\begin{subarray}{c}0\hskip
0.7113pt\leq\hskip 0.7113pti,\hskip 0.7113ptj\hskip 0.7113pt\leq\hskip
0.7113ptd\\\\[1.0pt] i\neq
j\end{subarray}}\mathbb{E}\Big{[}\big{(}Z^{(i)}\big{)}^{4}\big{(}Z^{(j)}\big{)}^{4}\Big{]}+4\sum_{\begin{subarray}{c}0\hskip
0.7113pt\leq\hskip 0.7113pti,\hskip 0.7113ptj\hskip 0.7113pt\leq\hskip
0.7113ptd\\\\[1.0pt] i\neq
j\end{subarray}}\mathbb{E}\Big{[}\big{(}Z^{(i)}\big{)}^{6}\big{(}Z^{(j)}\big{)}^{2}\Big{]}+\sum_{k\hskip
0.7113pt=\hskip
0.7113pt0}^{d}\mathbb{E}\Big{[}\big{(}Z^{(k)}\big{)}^{8}\Big{]}$
$\displaystyle=d(d-1)(d-2)(d-3)+18d(d-1)(d-2)+27d(d-1)+60d(d-1)+105d$
$\displaystyle=d^{4}+12d^{3}+44d^{2}+48d.$

Since the Brownian bridge on $[0,1]$ is a Gaussian process with covariance
function $K_{B}(s,t)=\min(s,t)-st$, it is clear that a rescaled bridge process
on $[t_{n},t_{n+1}]$ has variance

$\displaystyle\operatorname{Var}\big{(}B_{t_{n},t}\big{)}=\frac{(t_{n+1}-t)(t-t_{n})}{h_{n}}\hskip
0.7113pt,$

for $t\in[t_{n},t_{n+1}]$. Hence for $p\in\\{2,4,8\\}$, we have

$\displaystyle\big{\|}\hskip 0.7113ptB_{t_{n},t}\hskip
0.7113pt\big{\|}_{\mathbb{L}_{p}}=\bigg{(}\frac{(t_{n+1}-t)(t-t_{n})}{h_{n}}\bigg{)}^{\frac{1}{2}}\big{\|}Z\hskip
0.7113pt\big{\|}_{\mathbb{L}_{p}}\leq C(p)\sqrt{d\hskip 0.7113pt}\hskip
0.7113pt(t-t_{n})^{\frac{1}{2}}.$

for $t\in[t_{n},t_{n+1}]$. Finally, since $H_{n}$ and $K_{n}$ are independent
normal random variables with

$\displaystyle H_{n}\sim\mathcal{N}\bigg{(}0,\frac{1}{12}h_{n}\bigg{)},\hskip
14.22636ptK_{n}\sim\mathcal{N}\bigg{(}0,\frac{1}{720}h_{n}\bigg{)},$

we have
$\|H_{n}+6K_{n}\|_{\mathbb{L}_{p}}=\frac{\sqrt{30}}{15}(h_{n})^{\frac{1}{2}}\|Z\hskip
0.7113pt\|_{\mathbb{L}_{p}}$ for $n\geq 0$ and the result follows. ∎

######  Theorem A.10.

For $n\geq 0$, the triple iterated integral of the shifted Brownian bridge on
$[t_{n},t_{n+1}]$ can be expressed as

$\displaystyle\frac{1}{2}h_{n}^{2}\big{(}H_{n}+6K_{n}\big{)}=\int_{t_{n}}^{t_{n+1}}\int_{t_{n}}^{t}\widetilde{B}_{t_{k},s}\,ds\,dt.$

###### Proof.

The result follows by using integration by parts with the definitions of
$H_{n}$ and $K_{n}$.

$\displaystyle\frac{1}{2}h_{n}^{2}H_{n}+3h_{n}^{2}K_{n}$
$\displaystyle=\frac{1}{2}h_{n}\int_{t_{n}}^{t_{n+1}}B_{t_{n},t}\,dt+\int_{t_{n}}^{t_{n+1}}\bigg{(}\frac{1}{2}h_{n}-(t-t_{n})\bigg{)}B_{t_{n},t}\,dt+2h_{n}^{2}K_{n}$
$\displaystyle=\int_{t_{n}}^{t_{n+1}}\big{(}h_{n}-(t-t_{n})\big{)}B_{t_{n},t}\,dt+\int_{t_{n}}^{t_{n+1}}\int_{t_{n}}^{t}12K_{n}\bigg{(}\frac{s-t_{n}}{h_{n}}\bigg{)}\,ds\,dt$
$\displaystyle=\Bigg{[}\,\big{(}h_{n}-(t-t_{n})\big{)}\int_{t_{n}}^{t}B_{t_{n},s}\,ds\,\Bigg{]}_{t\hskip
0.7113pt=\hskip 0.7113ptt_{n}}^{t\hskip 0.7113pt=\hskip
0.7113ptt_{n+1}}-\int_{t_{n}}^{t_{n+1}}-\int_{t_{n}}^{t}B_{t_{n},s}\,ds\,dt$
$\displaystyle\hskip
42.67912pt+\int_{t_{n}}^{t_{n+1}}\int_{t_{n}}^{t}12K_{n}\bigg{(}\frac{s-t_{n}}{h_{n}}\bigg{)}\,ds\,dt$

$\displaystyle=\int_{t_{n}}^{t_{n+1}}\int_{t_{n}}^{t}\bigg{(}B_{t_{k},s}+12K_{n}\bigg{(}\frac{s-t_{n}}{h_{n}}\bigg{)}\bigg{)}\,ds\,dt$.
∎

Recall that each step of the shifted ODE approximation (over
$[t_{n},t_{n+1}]$) is given by

$\displaystyle\Bigg{(}\begin{matrix}\,x_{n+1}^{\prime}\\\\[-3.0pt]
\,v_{n+1}^{\prime}\end{matrix}\Bigg{)}:=\Bigg{(}\begin{matrix}\,\wideparen{x}^{\hskip
0.7113ptn}_{t_{n+1}}\\\\[-3.0pt] \,\wideparen{v}^{\hskip
0.7113ptn}_{t_{n+1}}\end{matrix}\Bigg{)}+12K_{n}\Bigg{(}\begin{matrix}0\\\\[-3.0pt]
\sigma\end{matrix}\,\Bigg{)},$

where $\big{\\{}\big{(}\wideparen{x}^{\hskip
0.7113ptn}_{t},\wideparen{v}^{\hskip
0.7113ptn}_{t}\big{)}\big{\\}}_{t\in[t_{n},t_{n+1}]}$ solves the following
ODE,

$\displaystyle\frac{d}{dt}\Bigg{(}\begin{matrix}\,\wideparen{x}^{\hskip
0.7113ptn}\\\\[-3.0pt] \,\wideparen{v}^{\hskip
0.7113ptn}\end{matrix}\Bigg{)}=\Bigg{(}\begin{matrix}\,\wideparen{v}^{\hskip
0.7113ptn}+\sigma\big{(}H_{n}+6K_{n}\big{)}\\\\[-3.0pt] \,-\gamma\big{(}\hskip
0.7113pt\wideparen{v}^{\hskip
0.7113ptn}+\sigma\big{(}H_{n}+6K_{n}\big{)}\big{)}-u\nabla
f\big{(}\,\wideparen{x}^{\hskip
0.7113ptn}\big{)}\end{matrix}\,\Bigg{)}+\frac{W_{n}-12K_{n}}{h_{n}}\,\Bigg{(}\begin{matrix}\,0\\\\[-3.0pt]
\,\sigma\end{matrix}\,\Bigg{)}\hskip 0.7113pt,$

with initial condition $\big{(}\hskip 0.7113pt\wideparen{x}^{\hskip
0.7113ptn}_{t_{n}}\hskip 0.7113pt,\wideparen{v}^{\hskip
0.7113ptn}_{t_{n}}\big{)}\hskip 0.7113pt:=\hskip 0.7113pt\big{(}x_{n}\hskip
0.7113pt,v_{n}\big{)}$. In addition, we will set $(\hskip
0.7113ptx_{0}^{\prime}\hskip 0.7113pt,v_{0}^{\prime}\hskip 0.7113pt)\hskip
0.7113pt:=\hskip 0.7113pt(x_{0}\hskip 0.7113pt,v_{0})$. We shall now present