###### Abstract

In this work we consider the unbiased estimation of expectations w.r.t.
probability measures that have non-negative Lebesgue density, and which are
known point-wise up-to a normalizing constant. We focus upon developing an
unbiased method via the underdamped Langevin dynamics, which has proven to be
popular of late due to applications in statistics and machine learning.
Specifically in continuous-time, the dynamics can be constructed to admit the
probability of interest as a stationary measure. We develop a novel scheme
based upon doubly randomized estimation as in [17, 18], which requires access
only to time-discretized versions of the dynamics and are the ones that are
used in practical algorithms. We prove, under standard assumptions, that our
estimator is of finite variance and either has finite expected cost, or has
finite cost with a high probability. To illustrate our theoretical findings we
provide numerical experiments which verify our theory, which include
challenging examples from Bayesian statistics and statistical physics.

Key words: Underdamped Langevin dynamics, unbiased estimation, maximal
couplings, Markov chain simulation
AMS subject classifications: 60J22, 65C05, 65C40, 82C31, 62G08, 35Q56
Code available at: https://github.com/ruzayqat/unbiased_uld
Corresponding author: Hamza Ruzayqat. E-mail<EMAIL_ADDRESS>
Unbiased Estimation using the Underdamped
Langevin Dynamics

HAMZA RUZAYQAT, NEIL K. CHADA, & AJAY JASRA

Applied Mathematics and Computational Science Program,
Computer, Electrical and Mathematical Sciences and Engineering Division,
King Abdullah University of Science and Technology, Thuwal, 23955-6900, KSA.

E-Mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
##  1 Introduction

We consider a class of probability measures on the measurable space
$(\mathbb{R}^{d},\mathcal{B}(\mathbb{R}^{d}))$ with non-negative Lebesgue
densities known point-wise up-to a normalizing constant. The objective of this
article is to consider simulation-based methods, which can return stochastic
estimates of finite expectations of functions w.r.t. the afore-mentioned
probability measure, that are unbiased, that is, on average are equal to the
expectation of interest. This latter task is of interest in a wide variety of
areas such as applied mathematics, physics and statistics; see e.g. [34] for a
book length introduction.

The primary methodology that is used in the literature to approximate
expectations is that of Markov chain Monte Carlo (MCMC) methods. These are
schemes which generate ergodic Markov chains whose stationary distribution is
exactly the one of interest and there are numerous variants of MCMC, such as
random walk Metropolis-Hastings, Hamiltonian Monte Carlo and non-reversible
MCMC; see e.g. [34]. In addition to this, are methods based upon uncorrected
time discretizations of continuous-time processes, which also have the
appropriate distribution of interest as a stationary distribution, such as the
underdamped Langevin algorithm e.g. [13, 30]. In particular with the
underdamped Langevin algorithm, it has been empirically observed to converge
with a better rate, to an invariant distribution, than that of the overdamped
Langevin dynamics, which is much simpler in comparison [7, 9, 12]. As a
result, these latter methods are of interest often due to their relative ease
of simulation relative to the former and have gained significant popularity in
the statistics and machine learning literature [14, 38]. In both of the
classes of methods that we have mentioned, in general, without starting these
Markov chains from draws from the target distribution of interest, one seldom
returns unbiased estimates which is the main interest in this article.
Unbiasedness can be desirable in certain contexts, for instance when computing
sensitivities for stochastic gradient algorithms e.g. [1].

We focus upon contributing to the class of unbiased MCMC algorithms and trying
to enhance the applicability of such schemes. Unbiased estimation can be
achieved in at least two ways, by exact simulation [2, 3] (e.g. starting the
chain from the correct probability of interest) or unbiased approximation. The
former has been investigated many years ago in the guise of coupling from the
past MCMC algorithms (e.g. [31]), but due to the complications of doing so,
the application of such methods is rather rare. The latter is often based upon
the pioneering works on unbiased estimation found in [16] (see also [27, 33]).
The idea of these latter methods, in the Markov chain context, is to work with
a pair of Markov chains on a product space and simulate them until they are
equal (the _meeting time_); there is then a novel identity which ensures that
the estimate is unbiased. Creating a methodology for allowing the Markov
chains to meet was developed in the paper [20] and subsequent to this, several
modifications and improvements were given in [28, 40]. However there has been
recent interest of unbiased methodologies, related to that latter way, where
we provide some of these works [5, 6, 17, 36, 37].

One of the main issues of the methodology developed in [20] and the sequels,
is that one needs to consider sometimes quite complex coupling of pairs of
Markov chains. This can be quite non-trivial to achieve and, in some scenarios
such as when the target density is multimodal, rather inefficient, leading to
large variances in estimation. This was partially addressed in [35] which
considered an unbiased version of the Schrödinger-Föllmer sampler (SFS). The
latter is a diffusion process on a bounded time domain $[0,1]$, that
transports a degenerate distribution at 0, to the target of interest, assuming
that the latter is absolutely continuous w.r.t. a $d-$dimensional standard
Gaussian. Even under Euler time-discretization, the process cannot be
simulated as the drift term is complicated resulting it in being intractable,
but several mechanisms are available; see [35]. The authors in that paper show
that by using doubly randomized unbiased schemes (e.g. [21, 22]) an unbiased
version of the SFS can be developed which provides unbiased estimates without
having to resort to complex coupling mechanisms. One of the drawbacks of the
methodology, however, is the SFS method in the beginning; when approximating
the drift, the method can struggle to well-represent complex probability
measures.

###  1.1 Contributions

The contribution of this article is to develop a new version of the
underdamped Langevin algorithm which, even when only working with time-
discretizations of the process, can deliver unbiased estimates of expectations
w.r.t. the class of probability measures under considerations. Our motivation
is that, like similarly to the work in [35], it does not require very complex
coupling techniques and also relies on a double randomization scheme which was
developed in [18]. The method also retains the advantages of the unbiased MCMC
methods in [20], that being that the estimates can be computed in an
embarrassingly parallel manner, which provides computational speed-ups versus
traditional MCMC algorithms. In particular from our work our highlighting
contribution is that we establish that our estimate is unbiased and of finite
variance and, either has finite expected cost to compute, or has finite cost
with high probability. These results rely heavily upon the works in [18] and
[11], where we assume similar assumptions for the former, and the latter of
which provides $\mathscr{V}$-uniform ergodicity of the discretized Markov
kernel that we use, where $\mathscr{V}$ is a chosen Lypaunov function. To
highlight our theoretical findings we provide numerical experiments on a range
of interesting and challenging examples, which arise in statistics and
physics. These include a Bayesian logistic regression problem, a double well
potential model and finally a Ginzburg-Landau model. We demonstrate both the
estimator achieving finite variance and unbiasedness. We compare it to the SFS
[35] to demonstrate the gains achieved from our proposed estimator.

###  1.2 Outline

This article is structured as follows. In Section 2 we detail our proposed
methodology based on the underdamped Langevin dynamics. This will lead onto
Section 3 where we establish that our estimator is unbiased and of finite
variance. Numerical experiments are then conducted in Section 4, where we
illustrate our methods on several challenging examples. We conclude our
remarks and future areas of research in Section 5. Finally the appendix houses
a technical result used in Section 3, and the majority of the algorithms
introduced.

##  2 Approach

In this section we first provide a common notation which will be used
throughout the manuscript. We then will introduce our approach for unbiased
estimation. This will include a review and discussion on the underdamped
Langevin dynamics, and in our context. Finally we will discuss how this is
related to our new unbiased estimator in the context of maximal couplings and
provide our unbiased estimator through various algorithms.

###  2.1 Notations

Let $(\mathsf{X},\mathcal{X})$ be a measurable space. For
$\varphi:\mathsf{X}\rightarrow\mathbb{R}$ we write
$\mathcal{B}_{b}(\mathsf{X})$, to denote the collection of bounded measurable
functions. For $\varphi\in\mathcal{B}_{b}(\mathsf{X})$, we write the supremum
norm as $\|\varphi\|=\sup_{x\in\mathsf{X}}|\varphi(x)|$. We denote the Borel
sets on $\mathbb{R}^{d}$ as $B(\mathbb{R}^{d})$. The $d-$dimensional Lebesgue
measure is written as $dx$. For a metric
$\mathsf{d:\mathsf{X}\times\mathsf{X}}\rightarrow\mathbb{R}^{+}$ on
$\mathsf{X}$ and a function $\varphi:\mathsf{X}\rightarrow\mathbb{R}$,
$\textrm{Lip}_{\mathsf{d}}(\mathsf{X})$ are the Lipschitz functions (with
finite Lipschitz constants), that is for every
$(x,w)\in\mathsf{X}\times\mathsf{X}$,
$|\varphi(x)-\varphi(w)|\leq\|\varphi\|_{\textrm{Lip}}\mathsf{d}(x,w)$.
$\mathscr{P}(\mathsf{X})$ denotes the collection of probability measures on
$(\mathsf{X},\mathcal{X})$. For a finite measure $\mu$ on
$(\mathsf{X},\mathcal{X})$ and a $\varphi\in\mathcal{B}_{b}(\mathsf{X})$, the
notation $\mu(\varphi)=\int_{\mathsf{X}}\varphi(x)\mu(dx)$ is used. For
$(\mathsf{X}\times\mathsf{W},\mathcal{X}\vee\mathcal{W})$ a measurable space
and $\mu$ a non-negative finite measure on this space, we use the tensor-
product of functions notation for
$(\varphi,\psi)\in\mathcal{B}_{b}(\mathsf{X})\times\mathcal{B}_{b}(\mathsf{W})$,
$\mu(\varphi\otimes\psi)=\int_{\mathsf{X}\times\mathsf{W}}\varphi(x)\psi(w)\mu(d(x,w))$.
Given a Markov kernel $K:\mathsf{X}\rightarrow\mathscr{P}(\mathsf{X})$ and a
finite measure $\mu$, we use the notations $\mu
K(dx^{\prime})=\int_{\mathsf{X}}\mu(dx)K(x,dx^{\prime})$ and
$K(\varphi)(x)=\int_{\mathsf{X}}\varphi(x^{\prime})K(x,dx^{\prime}),$ for
$\varphi\in\mathcal{B}_{b}(\mathsf{X})$. The iterated kernel is
$K^{n}(x_{0},dx_{n})=\int_{\mathsf{X}^{n-1}}\prod_{i=1}^{n}K(x_{i-1},dx_{i})$.
For $A\in\mathcal{X}$, the indicator function is written as
$\mathbb{I}_{A}(x)$. $\mathbb{Z}^{+}$ is the set of non-negative integers.
$I_{d}$ denotes the $d\times d$ identity matrix. The transpose of a vector or
matrix $x$ is denoted as $x^{T}$. We denote min($a,b$) as $a\wedge b$. The
operator $\nabla$ denotes the gradient and $\Delta$ the Laplacian, while
$\Delta_{l}$ denotes a time step-size of $2^{-l}$, $l\in\mathbb{Z}^{+}$.

###  2.2 Framework

We consider the case of underdamped Langevin dynamics:

$\displaystyle dX_{t}$ $\displaystyle=$ $\displaystyle V_{t}dt,$ (2.1)
$\displaystyle dV_{t}$ $\displaystyle=$ $\displaystyle\left(b(X_{t})-\kappa
V_{t}\right)dt+\sigma dB_{t},$ (2.2)

where $\\{B_{t}\\}_{t\geq 0}$ is a standard $d-$dimensional Brownian motion,
$(\kappa,\sigma)\in(0,\infty)^{2}$ are given friction and diffusion
coefficients and $b:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ is a typically of
gradient form, $b=-\nabla U$. In the latter case there is, under fairly weak
assumptions an invariant measure of the process $\\{X_{t},V_{t}\\}_{t\geq 0}$,
with Lebesgue density, $\pi(x,v)$:

$\pi(x,v)\propto\exp\left\\{\frac{-\kappa\left(2U(x)+\|v\|^{2}\right)}{\sigma^{2}}\right\\}.$

In practice, one often resorts to time-discretization of the dynamics
(2.1)-(2.2), for which we focus upon the Euler discretization of step-size
$\Delta_{l}=2^{-l}$, $l\in\mathbb{N}_{0}$:

$\displaystyle X_{(k+1)\Delta_{l}}$ $\displaystyle=$ $\displaystyle
X_{k\Delta_{l}}+V_{k\Delta_{l}}\Delta_{l},$ (2.3) $\displaystyle
V_{(k+1)\Delta_{l}}$ $\displaystyle=$ $\displaystyle
V_{k\Delta_{l}}+\left(b(X_{k\Delta_{l}})-\kappa
V_{k\Delta_{l}}\right)\Delta_{l}+\sigma\left(B_{(k+1)\Delta_{l}}-B_{k\Delta_{l}}\right),$
(2.4)

with $k\in\mathbb{N}_{0}$. Under assumptions (see [11]), for $l$ large enough,
the discrete-time Markov chain expressed as
$\\{X_{k\Delta_{l}},V_{k\Delta_{l}}\\}_{k\in\mathbb{N}_{0}}$ admits a unique
invariant measure $\eta_{l}$ and moreover will converge (in an appropriate
sense) to it geometrically quickly. In addition, we should have that
$\eta_{l}$ will converge to $\pi$, by the convergence of Euler approximations.
As it will facilitate the approach we are to introduce, we shall modify (2.3)
to

$X_{(k+1)\Delta_{l}}=X_{k\Delta_{l}}+V_{k\Delta_{l}}\Delta_{l}+\sigma_{l}\Gamma_{k,l},$
(2.5)

where, $\\{\sigma_{l}\\}_{l\in\mathbb{N}_{0}}$ is any sequence of non-negative
and decreasing constants, that converge to zero and for each
$l\in\mathbb{N}_{0}$, $\\{\Gamma_{k,l}\\}_{k\in\mathbb{N}_{0}}$ is a sequence
of i.i.d. $d-$dimensional Gaussian random variables of zero mean and
covariance matrix the identity multiplied by $\Delta_{l}$ and that this family
of sequences are independent of all other random variables (and of each
other). As this modification falls under the framework of [11], under
assumptions, that for $l$ large enough, the discrete-time Markov chain
$\\{X_{k\Delta_{l}},V_{k\Delta_{l}}\\}_{k\in\mathbb{N}_{0}}$, following the
dynamics (2.4)-(2.5), admits a unique invariant measure $\pi_{l}$ (which is
likely to be different to $\eta_{l}$) and moreover will converge (in an
appropriate sense) to $\pi_{l}$ geometrically quickly. In addition, we should
have that $\pi_{l}$ will converge to $\pi$, by the convergence of Euler
approximations and that the noise in (2.5) will disappear.

The approach here is as follows: by using only the dynamics (2.4)-(2.5), we
will show that one can produce an estimator of $\pi(\varphi\otimes
1)=\int_{\mathbb{R}^{2d}}\varphi(x,v)\pi(x,v)dxdv$ that is unbiased, where
$\varphi:\mathbb{R}^{2d}\rightarrow\mathbb{R}$ is $\pi-$integrable. In
practice, of course, we will only be interested in the marginal on the $X-$co-
ordinate of $\pi$. In other words, we will be interested in unbiasedly
estimating the integral $\int_{\mathbb{R}^{d}}\Psi(x)\tilde{\pi}(x)dx$, with
$\Psi(x):\mathbb{R}^{d}\to\mathbb{R}$ is $\tilde{\pi}-$integrable and
$\tilde{\pi}(x)\propto\exp{\\{-(2\kappa/\sigma^{2})U(x)\\}}$ for some
differentiable potential function $U$.

###  2.3 A Conditional-Type Max-Coupling

In computational statistics a common procedure is to couple various
distributions $\mu_{1}$ and $\mu_{2}$. Particular examples of this can include
independent coupling, or optimal couplings based on the theory of optimal
transport. For this work we make use of maximum couplings [39].

To construct our approach, we will need the following simple idea. Suppose
that we have two positive Lebesgue densities $\mu_{1}(x_{1},y)$ and
$\mu_{2}(x_{2},y)$ on spaces $\mathsf{X}_{1}\times\mathbb{R}^{2d}$ and
$\mathsf{X}_{2}\times\mathbb{R}^{2d}$, where $\mathsf{X}_{1}$ and
$\mathsf{X}_{2}$ are two possibly different dimensional sub-spaces of powers
of the real line. Suppose also, that we know pointwise the conditional
densities, for $j\in\\{1,2\\}$, $x_{j}\in\mathsf{X}_{j}$ fixed:

$\mu_{j}(y|x_{j})=\frac{\mu_{j}(x_{j},y)}{\int_{\mathbb{R}^{d}}\mu_{j}(x_{j},y)dy}.$

Our objective is to sample from a coupling of $\mu_{1}$ and $\mu_{2}$, so that
there is a non-zero probability that the $Y-$co-ordinate can be equal. Let
$\check{\mu}$ be any coupling of the marginals $\mu_{1}(x_{1})$ and
$\mu_{2}(x_{2})$ and $\overline{\mu}(d(y,y^{\prime})|x_{1},x_{2})$ be the
maximal coupling of $\mu_{1}(\cdot|x_{1})$ and $\mu_{2}(\cdot|x_{2})$ then one
way to achieve our given objective is to sample from the distribution
associated to the probability:

$\check{\mu}\left(d(x_{1},x_{2})\right)\overline{\mu}\left(d(y,y^{\prime})|x_{1},x_{2}\right),$

which can be broadly thought of as a posterior. Note that it is straight-
forward to show that the marginal of $(y,x_{j})$ is $\mu_{j}$ for
$j\in\\{1,2\\}$. This concept can be straight-forwardly extended to the case
of 4 targets, as in [18, Section 3.2.2.], using much the same construction.

###  2.4 Method

The method that we pursue is an adaptation of the approach in [18] which was
originally developed for unbiased inference for Bayesian inverse problems. The
contribution here is that we provide an original coupling construction that
one might expect is easier to apply than to conventional Markov kernels such
as Metropolis-Hastings. In addition, to verify the unbiasedness, one needs to
merge the theory of [11] and [18]. The former will be discussed later, with
various assumptions stated before the main theorem is presented. We will now
discuss our unbiased strategy, based on the work in [18].

####  2.4.1 Overall Strategy

Let $l_{*}\in\mathbb{N}_{0}$ be given and $\mathbb{P}_{L}$ be any positive
probability mass function on $\mathbb{N}_{l_{*}}:=\\{l_{*},l_{*}+1,\dots\\}$.
Let $\\{\xi_{l}\\}_{l\in\mathbb{N}_{l_{*}}}$ be any sequence of independent
random variables, such that

$\displaystyle\mathbb{E}[\xi_{l_{*}}]$ $\displaystyle=$
$\displaystyle\pi_{l_{*}}(\varphi)$ $\displaystyle\mathbb{E}[\xi_{l}]$
$\displaystyle=$
$\displaystyle\pi_{l}(\varphi)-\pi_{l-1}(\varphi)=:[\pi_{l}-\pi_{l-1}](\varphi)\quad
l\in\\{l_{*}+1,l_{*}+2,\dots\\}.$

Now, let $L$ be a random variable with probability $\mathbb{P}_{L}$ that is
independent of the sequence $\\{\xi_{l}\\}_{l\in\mathbb{N}_{l_{*}}}$ then

$\displaystyle\widehat{\pi}(\varphi)=\frac{\xi_{l}}{\mathbb{P}_{L}(l)},$ (2.6)

is an unbiased estimator of $\pi(\varphi)$; see [27, 33] for the initial
statement and proof. Moreover, if

$\sum_{l\in\mathbb{N}_{l_{*}}}\frac{\mathbb{E}[\xi_{l}^{2}]}{\mathbb{P}_{L}(l)}<+\infty,$
(2.7)

the estimator $\widehat{\pi}(\varphi)$ has finite variance. There is also the
independent sum-estimator, which can work better than this estimator, but we
shall not discuss that for now. The main challenge is then to construct the
sequence $\\{\xi_{l}\\}_{l\in\mathbb{N}_{l_{*}}}$.

Typically, one will run $M\in\mathbb{N}$ independent replicates of (2.6),
where for each replicate $i$, $l_{i}\sim\mathbb{P}_{L}$, and then use the
average

$\left(\widehat{\pi}(\varphi)\right)_{\text{avg}}:=\frac{1}{M}\sum_{i=1}^{M}\left(\widehat{\pi}(\varphi)\right)^{(i)},$

where $\left(\widehat{\pi}(\varphi)\right)^{(i)}$ represents the $i$-th
independent replicate of the estimate.

####  2.4.2 Unbiased Approximation of $\pi_{l_{*}}(\varphi)$

Throughout the section $l\in\mathbb{N}_{l_{*}}$ is fixed; although we are
interested in the case $l=l_{*}$ the subsequent exposition holds for any fixed
$l$. We will use a strategy that was developed in [16] (see also [20]) for
obtaining the given estimator. The approach is to construct a Markov chain on
the product space $\mathsf{U}^{2}$, where $\mathsf{U}=\mathbb{R}^{2d}$. From
herein, we consider the Markov transition associated to (2.4)-(2.5) over a
unit time interval. We denote this kernel as
$K_{l}:\mathsf{U}\rightarrow\mathscr{P}(\mathsf{U})$, where $\mathscr{P}(U)$
are the collection of probability measures on measurable space
$(\mathsf{U},\mathcal{U})$, $\mathcal{U}$ is the Borel $\sigma-$field on
$\mathbb{R}^{2d}$. The reason for this, will become apparent as we continue in
our exposition.

In order to follow the construction in [16] we will need to construct a Markov
chain, $\\{U_{n,l,}\tilde{U}_{n,l}\\}_{n\in\mathbb{N}_{0}}$, on
$(\mathsf{U}^{2},\mathcal{U}\vee\mathcal{U})$ so that marginally, at a given
discrete time, new positions of the chain have the kernel $K_{l}$ as a
marginal. More precisely, we need to construct a Markov kernel,
$\check{K}_{l}:\mathsf{U}^{2}\rightarrow\mathscr{P}(\mathsf{U}^{2})$, so that
for any $(u_{l},\tilde{u}_{l},A)\in\mathsf{U}^{2}\times\mathcal{U}$:

$\int_{A\times\mathsf{U}}\check{K}_{l}\left((u_{l},\tilde{u}_{l}),d(u_{l}^{\prime},\tilde{u}_{l}^{\prime})\right)=\int_{A}K_{l}(u_{l},du_{l}^{\prime})\quad\textrm{and}\quad\int_{\mathsf{U}\times
A}\check{K}_{l}\left((u_{l},\tilde{u}_{l}),d(u_{l}^{\prime},\tilde{u}_{l}^{\prime})\right)=\int_{A}K_{l}(\tilde{u}_{l},d\tilde{u}_{l}^{\prime}).$
(2.8)

To build the coupling $\check{K}_{l}$ we shall decompose this into a mixture
of two kernels
$\check{Q}_{l}:\mathsf{U}^{2}\rightarrow\mathscr{P}(\mathsf{U}^{2})$ and
$\check{P}_{l}:\mathsf{U}^{2}\rightarrow\mathscr{P}(\mathsf{U}^{2})$. We will
explain why this is required, as the discussion progresses. The simulation of
the first kernel $\check{Q}_{l}$ is described in Algorithm 2. On inspection of
Algorithm 2, it is clear that this samples from a coupling of $K_{l}$ which is
such that if $u=\tilde{u}$ then $u^{\prime}=\tilde{u}^{\prime}$. We note
however, that if $u\neq\tilde{u}$, then there is zero probability that
$u^{\prime}=\tilde{u}^{\prime}$ and being able to achieve the latter is
critical in our construction. This motivates our second kernel $\check{P}_{l}$
whose simulation is given in Algorithm 3.

The kernel $\check{P}_{l}$ is essentially the same as $\check{Q}_{l}$ except
at the final time step, one has a non-zero probability that
$u^{\prime}=\tilde{u}^{\prime}$ _irrespective_ of $u,\tilde{u}$ and note that
if $u=\tilde{u}$ one will always simulate $u^{\prime}=\tilde{u}^{\prime}$.
Note also, that it easily follows that $\check{Q}_{l}$ samples from a coupling
of $K_{l}$; the kernel is a simple application of the method in Section 2.3.
In Algorithm 3, step 4. the maximal coupling can be sampled, for instance,
using the maximal coupling algorithm in [39] or the reflection maximal
coupling described in [20], where the later performs better for high-
dimensional models and therefore we adopt it. Then, for some $\alpha\in(0,1)$
fixed we set

$\check{K}_{l}=\alpha\check{Q}_{l}+(1-\alpha)\check{P}_{l}.$

The reason for including the kernel $\check{Q}_{l}$ is that it encourages the
positions of $u_{n,l},\tilde{u}_{n,l}$ to be close, but can never achieve
equality; this latter task can be achieved by $\check{P}_{l}$ and is why this
latter kernel is used.

To initialize the Markov chain, for some initial probability
$\nu_{l}\in\mathscr{P}(\mathsf{U})$ we use the following idea as in [16]. Let
$\check{\nu}_{l}$ be any coupling of $(\nu_{l},\nu_{l})$ then the initial
probability is taken as:

$\overline{\nu}_{l}\left(d(u^{\prime},\tilde{u}^{\prime})\right)=\int_{\mathsf{U}^{2}}\check{\nu}_{l}\left(d(u,\tilde{u})\right)K_{l}(u,du^{\prime})\delta_{\\{\tilde{u}\\}}(d\tilde{u}^{\prime}),$

which evolves from the Markov kernel $\check{K}_{l}$ and satisfies the
properties discussed through (2.8).

Recall we are interested in producing an unbiased estimator of the quantity

$\pi_{l}(\varphi):=\varphi(U_{k,l})+\sum_{n=k+1}^{\infty}\\{\varphi(U_{n,l})-\varphi({U}_{n-1,l})\\}.$
(2.9)

However doing so can be challenging, in particular because we have an infinite
sum. Therefore what can be done instead, as eluded too from the discussion
above, is to construct another chain $(\tilde{U}_{n,l})$ such that
$U_{n,l}=\tilde{U}_{n,l}$, for $n\geq\tau$, where when the two chains meet,
this is known as the meeting time $\tau_{l}$, which is defined as follows

$\tau_{l}:=\inf\\{n\geq 1:U_{n,l}=\tilde{U}_{n,l}\\},$

where we require that $\tau_{l}$ is almost surely finite and that for each
$n\geq\tau_{l}$, $U_{n,l}=\tilde{U}_{n,l}$ almost surely. Therefore what we
can do is replace the infinite limit with the meeting time $\tau_{l}$ with a
lag of 1, where now we have the following unbiased estimator of
$\pi_{l}(\varphi)$ for any $k\in\mathbb{N}_{0}$

$\widehat{\pi_{l}(\varphi)}_{k}:=\varphi(U_{k,l})+\sum_{n=k+1}^{\tau_{l}-1}\\{\varphi(U_{n,l})-\varphi(\tilde{U}_{n,l})\\},$
(2.10)

with
$\sum_{n=k+1}^{\tau_{l}-1}\\{\varphi(U_{n,l})-\varphi(\tilde{U}_{n,l})\\}=0$
if $\tau_{l}-1<k+1$. In order to prove that (2.10) is an unbiased estimator,
with finite variance and expect cost, we summarize these weak assumptions,
without going into exact details.

  * 1.

We require convergence of the marginal chains, i.e,

$\lim_{n\rightarrow\infty}\mathbb{E}[\varphi(U_{n,l})]=\mathbb{E}_{\pi}[\varphi(U)].$

  * 2.

The meeting time $\tau_{l}:=\inf\\{n\geq 1:U_{n,l}=\tilde{U}_{n,l}\\}$ has
geometric tails

$\mathbb{P}(\tau>n)\leq C\rho^{n},\quad\textrm{for}\quad C<\infty,\
\rho\in(0,1).$

  * 3.

Faithfulness: once both chains meet at the meeting time $\tau_{l}$, then
$U_{n,l}=\tilde{U}_{n,l}$ for $n\geq\tau_{l}$.

We remark that a time-averaged extension is also possible; let
$(m,k)\in\mathbb{Z}^{+}\times\mathbb{Z}^{+}$, with $m\geq k$, then one can
also use (see [20]) the time-averaged estimator

$\widehat{\pi_{l}(\varphi)}_{T,k,m}:=\frac{1}{m-k+1}\sum_{n=k}^{m}\varphi(U_{n,l})+\sum_{n=k+1}^{\tau_{l}-1}\Big{(}1\wedge\frac{n-k}{m-k+1}\Big{)}\\{\varphi(U_{n,l})-\varphi(\tilde{U}_{n,l})\\},$
(2.11)

which recovers (2.10) in the case $m=k$.

Figure 1: Example illustration of the time-averaged estimator
$\widehat{\pi_{l}(\varphi)}_{T,k,m}$ from (2.11), which includes the meeting
time $\tau=7$ with respect to the two chains $U$ and $\tilde{U}$. Each line
corresponds to a chain, for which we are aiming to couple.

####  2.4.3 Unbiased Approximation of $[\pi_{l}-\pi_{l-1}](\varphi)$

Our objective is now to provide, for $l\in\\{l_{*}+1,l*+2,\dots\\}$ fixed, an
unbiased estimator of $[\pi_{l}-\pi_{l-1}](\varphi)$, which will essentially
use the approach detailed in the previous section. Indeed, one could simply
use the method outlined above, independently, for $\pi_{l}$ and $\pi_{l-1}$
and independently for each $l\in\\{l_{*}+1,l*+2,\dots\\}$. However, this is
unlikely to provide an estimator that can achieve (2.7) and hence the variance
of such an approach is infinite and not useful in practice. We therefore
present an alternative method.

The idea we use is to generate a Markov chain on $(\mathsf{Z},\mathcal{Z})$,
where $\mathsf{Z}=\mathsf{U}^{4}$ and
$\mathcal{Z}=(\mathcal{U}\vee\mathcal{U})\vee(\mathcal{U}\vee\mathcal{U})$ and
we write

$Z_{n,l,l-1}=\left((U_{n,l},\tilde{U}_{n,l}),(U_{n,l-1},\tilde{U}_{n,l-1})\right).$

The associated Markov kernel
$\check{K}_{l,l-1}:\mathsf{Z}\rightarrow\mathscr{P}(\mathsf{Z})$ will be
developed below. The construction of our Markov chain, should be that,
marginally $\\{U_{n,l}\\}_{n\in\mathbb{N}_{0}}$ (resp.
$\\{\tilde{U}_{n,l}\\}_{n\in\mathbb{N}_{0}}$) is a Markov chain with kernel
$K_{l}$ and marginally $\\{U_{n,l-1}\\}_{n\in\mathbb{N}_{0}}$ (resp.
$\\{\tilde{U}_{n,l-1}\\}_{n\in\mathbb{N}_{0}}$) is a Markov chain with kernel
$K_{l-1}$. It is explicitly assumed that (at the very least)
$(Z_{n,l,l-1})_{n\in\mathbb{N}_{0}}$ is constructed so that the stopping time
$\check{\tau}_{l,l-1}:=\tau_{l}\vee\tau_{l-1}$ is almost surely finite. In
addition, the pair of chains on each level should be faithful, i.e. for
$s\in\\{l,l-1\\}$, we have

$\displaystyle U_{n,s}=\tilde{U}_{n,s},\mbox{ for all }n\geq\tau_{s}.$ (2.12)

Hence for time $n\geq\check{\tau}_{l,l-1}$, $Z_{n,l,l-1}$ only has a distinct
state on each level.

Our Markov kernel $\check{K}_{l,l-1}$ will be a type of mixture constituting
two distinct Markov kernels
$\check{Q}_{l,l-1}:\mathsf{Z}\rightarrow\mathscr{P}(\mathsf{Z})$ and
$\check{P}_{l,l-1}:\mathsf{Z}\rightarrow\mathscr{P}(\mathsf{Z})$. We describe
the simulation of $\check{Q}_{l,l-1}$ in Algorithm 5. This kernel samples, for
$s\in\\{l,l-1\\}$, from a coupling of $K_{s}$ which is such that if
$u_{s}=\tilde{u}_{s}$ then $u_{s}^{\prime}=\tilde{u}_{s}^{\prime}$ and if
$u_{s}\neq\tilde{u}_{s}$, then there is zero probability that
$u_{s}^{\prime}=\tilde{u}_{s}^{\prime}$. The simulation is also such that the
outputs across the levels should be quite dependent and, in some sense, close.
The additional kernel $\check{P}_{l,l-1}$ that we require is described in
Algorithm 6. The kernel $\check{P}_{l,l-1}$ is essentially the same as
$\check{Q}_{l,l-1}$ except at the final time step, one has a non-zero
probability that, for $s\in\\{l,l-1\\}$,
$u_{s}^{\prime}=\tilde{u}_{s}^{\prime}$ _irrespective_ of
$u_{s},\tilde{u}_{s}$ and note that if $u_{s}=\tilde{u}_{s}$ one will always
simulate $u_{s}^{\prime}=\tilde{u}_{s}^{\prime}$. Note also, that it easily
follows that $\check{Q}_{l,l-1}$ samples from a couplings of $K_{l}$ and
$K_{l-1}$. In Algorithm 6, step 4. the synchronous pairwise reflection maximal
coupling is described in [4, 18] and the simulation thereof is similar to the
reflection maximal coupling. Then for $\alpha\in(0,1)$ we set

$\check{K}_{l,l-1}=\mathbb{I}_{D^{2}}(u_{l},\tilde{u}_{l},u_{l-1},\tilde{u}_{l-1})\check{Q}_{l,l-1}+\mathbb{I}_{(D^{2})^{c}}(u_{l},\tilde{u}_{l},u_{l-1},\tilde{u}_{l-1})[\alpha\check{Q}_{l,l-1}+(1-\alpha)\check{P}_{l,l-1}],$

where $D=\\{(u,\tilde{u})\in\mathsf{U}^{2}:u=\tilde{u}\\}$.

To initialize the Markov chain, for some initial probabilities,
$s\in\\{l,l-1\\}$, $\nu_{s}\in\mathscr{P}(\mathsf{U})$ we use the following
idea. We need a kernel
$\overline{K}_{l,l-1}:\mathsf{U}^{2}\rightarrow\mathscr{P}(\mathsf{U}^{2})$
described in Algorithm 4 for the initialization. Let $\check{\nu}_{s}$ be any
coupling of $(\nu_{s},\nu_{s})$ and $\check{\nu}_{l,l-1}$ be a coupling of
$(\check{\nu}_{l},\check{\nu}_{l-1})$ then the initial probability is taken
as:

$\overline{\nu}_{l,l-1}\left(dz_{l,l-1}^{\prime}\right)=\int_{\mathsf{U}^{4}}\check{\nu}_{l,l-1}\left(dz_{l,l-1}\right)\overline{K}_{l,l-1}\left((u_{l},u_{l-1}),d(u_{l}^{\prime},u_{l-1}^{\prime})\right)\delta_{\\{\tilde{u}_{l},\tilde{u}_{l-1}\\}}\left(d(\tilde{u}_{l}^{\prime},\tilde{u}_{l-1})\right).$

The process is then sampled at subsequent time-points, using the Markov kernel
$\check{K}_{l,l-1}$.

Finally then one can estimate $[\pi_{l}-\pi_{l-1}](\varphi)$ as follows, for
any $k\in\mathbb{N}_{0}$:

$\widehat{[\pi_{l}-\pi_{l-1}](\varphi)}_{k}:=\widehat{\pi_{l}(\varphi)}_{k}-\widehat{\pi_{l-1}(\varphi)}_{k},$
(2.13)

where $\widehat{\pi_{s}(\varphi)}_{k}$ is computed using (2.10) based on the
pair of chains $\\{U_{n,s},\tilde{U}_{n,s}\\}_{n\in\mathbb{Z}^{+}}$ on level
$s\in\\{l,l-1\\}$. One can also employ time-averaging, for
$(m,k)\in\mathbb{N}_{0}\times\mathbb{N}_{0}$ satisfying $m\geq k$:

$\displaystyle\widehat{[\pi_{l}-\pi_{l-1}](\varphi)}_{T,k,m}$
$\displaystyle:=\widehat{\pi_{l}(\varphi)}_{T,k,m}-\widehat{\pi_{l-1}(\varphi)}_{T,k,m},$
(2.14)

where $\widehat{\pi_{s}(\varphi)}_{T,k,m}$ is computed using (2.11) based on
the pair of chains $\\{U_{n,s},\tilde{U}_{n,s}\\}_{n\in\mathbb{N}_{0}}$ on
level $s\in\\{l,l-1\\}$.

####  2.4.4 Final Methodology and Estimator

We now consolidate the above discussion by summarizing our proposed
methodology to unbiasedly estimate $\pi(\varphi)$.

Algorithm 1 Unbiased estimator $\widehat{\pi(\varphi)}$.

Input: $\mathbb{P}_{L}$,
$(u_{0},\tilde{u}_{0})=\left((x_{0},v_{0}),(\tilde{x}_{0},\tilde{v}_{0})\right)\in\mathsf{U}^{2}$
and
$(u_{0,s},\tilde{u}_{0,s})=\left((x_{0,s},v_{0,s}),(\tilde{x}_{0,s},\tilde{v}_{0,s})\right)\in\mathsf{U}^{2}$,
$s\in\\{F,C\\}$.

  1. 1.

Sample $L\sim\mathbb{P}_{L}$.

  2. 2.

If $L=l_{*}$, initialize the chains by setting
$\tilde{U}_{0,l_{*}}=\tilde{u}_{0}$ and sampling $U_{0,l_{*}}$ by running the
dynamics in (2.4) - (2.5) at level $l_{*}$ up to time 1 with initial values
$(x_{0},v_{0})$. Then generate a Markov chain
$\\{U_{n,l},\tilde{U}_{n,l}\\}_{n\in\mathbb{N}}$ according to
$\check{K}_{l_{*}}$, and compute the estimator
$\widehat{\pi_{0}(\varphi)}_{k}$ in (2.10) or
$\widehat{\pi_{0}(\varphi)}_{T,k,m}$ in (2.11).

  3. 3.

If $L>l_{*}$, initialize the chains by setting
$(\tilde{U}_{0,L},\tilde{U}_{0,L-1})=(\tilde{u}_{0,F},\tilde{u}_{0,C})$ and
sampling $(U_{0,L},U_{0,L-1})$ from $\overline{K}_{L,L-1}$ by running
Algorithm 4 with initial values $(u_{0,F},u_{0,C})$. Then generate a Markov
chain $\\{Z_{n,L,L-1}\\}_{n\in\mathbb{N}}$ according to $\check{K}_{L,L-1}$,
and compute the estimator $\widehat{[\pi_{L}-\pi_{L-1}](\varphi)}_{k}$ in
(2.13) or $\widehat{[\pi_{L}-\pi_{L-1}](\varphi)}_{T,k,m}$ in (2.14).

Return the single-term estimator

$\widehat{\pi(\varphi)}_{S,k}:=\frac{1}{\mathbb{P}_{L}(L)}\Big{(}\mathbb{I}_{\\{l_{*}\\}}(L)\leavevmode\nobreak\
\widehat{\pi_{l_{*}}(\varphi)}_{k}+\mathbb{I}_{\\{l_{*}+1,l_{*}+2,\dots\\}}(L)\leavevmode\nobreak\
\widehat{[\pi_{L}-\pi_{L-1}](\varphi)}_{k}\Big{)},$ (2.15)

or the time-averaged estimator

$\widehat{\pi(\varphi)}_{S,T,k,m}:=\frac{1}{\mathbb{P}_{L}(L)}\Big{(}\mathbb{I}_{\\{l_{*}\\}}(L)\leavevmode\nobreak\
\widehat{\pi_{l_{*}}(\varphi)}_{T,k,m}+\mathbb{I}_{\\{l_{*}+1,l_{*}+2,\dots\\}}(L)\leavevmode\nobreak\
\widehat{[\pi_{L}-\pi_{L-1}](\varphi)}_{T,k,m}\Big{)}.$ (2.16)

##  3 Theoretical Properties

In this section we provide theory related to the discussed methodology in
Section 2. Specifically we provide theoretical justification, where we provide
our main result which states that our estimator has finite variance, and
either has finite expected cost, or has finite cost with a high probability.
In order to prove this we require a number of standard assumptions, which have
been used in previous works. We omit the proof of theorem in this section,
where we defer it to the Appendix.

###  3.1 Assumptions

Before discussing our required assumptions, we provide some new common
notation in the setting of Lyapunov functions. Let $0<C<+\infty$ be a constant
and $\mathsf{d}$ a metric on $\mathsf{U}$; we define the set

$\displaystyle\mathsf{B}(C,\Delta_{l},\mathsf{d})$ $\displaystyle:=$
$\displaystyle\\{u_{1:4}\in\mathsf{U}^{4}:\forall(i,j)\in\\{1,\dots,4\\},\mathsf{d}(u_{i},u_{j})\leq
C\Delta_{l}\\}.$

In what follows, for any $l\in\mathbb{N}_{l_{*}}$, under assumptions (see
[11]) the Markov kernel $K_{l}$ admits a unique invariant measure $\pi_{l}$.
We introduce a Lyapunov function $\mathscr{V}:\mathsf{U}\rightarrow[1,\infty)$
which will be used in our assumptions below. For
$f:\mathsf{U}\rightarrow\mathbb{R}$ and a given Lyapunov function
$\mathscr{V}$, we define the collection of functions
$\mathcal{L}_{\mathscr{V}}:=\\{f:\sup_{u\in\mathsf{U}}|f(u)|/\mathscr{V}(u)<+\infty\\}$
and we write $\|f\|_{\mathscr{V}}=|f(u)|/\mathscr{V}(u)$.

  * (A1)

    1. 1.

There exist a $\mathscr{V}:\mathsf{U}\rightarrow[1,\infty)$,
$(\lambda,b,\mathsf{C})\in(0,1)\times(0,\infty)\times\mathcal{U}$ such that
for any $(l,u)\in\mathbb{N}_{l_{*}}\times\mathsf{U}$

$K_{l}(\mathscr{V})(u)\leq\lambda\mathscr{V}(u)+b\mathbb{I}_{\mathsf{C}}(u).$

    2. 2.

There exists a $(\rho,C)\in(0,1)\times(0,\infty)$ such that for any
$n\in\mathbb{N}$

$\sup_{l\in\mathbb{N}_{l_{*}}}\sup_{u\in\mathsf{U}}\sup_{\varphi\in\mathcal{L}_{\mathscr{V}}}\frac{|K_{l}^{n}(\varphi)(u)-\pi_{l}(\varphi)|}{\mathscr{V}(u)}\leq
C\|\varphi\|_{\mathscr{V}}\rho^{n},$

with $\mathscr{V}$ as in 1..

  * (A2)

There exist $(C,\rho)\in(0,\infty)\times(0,1)$ such that:

    1. 1.

For any $n\in\mathbb{N}$

$\mathbb{E}[\mathbb{I}_{\\{\tau_{0}>n\\}}]\leq C\rho^{n}.$

    2. 2.

For any $(l,n,z)\in\mathbb{N}\times\mathbb{N}\times\mathsf{X}^{4}$

$\mathbb{E}[\mathbb{I}_{\\{\tau_{l}>n\\}}|Z_{0,l,l-1}=z]\leq C\rho^{n}.$

  * (A3)

There exist a $C<\infty$ and metric on $\mathsf{X}$,
$\tilde{\mathsf{d}}:\mathsf{X}\times\mathsf{X}\rightarrow\mathbb{R}^{+}$ such
that for any
$(l,\varphi,(x,y))\in\mathbb{Z}^{+}\times\mathcal{B}_{b}(\mathsf{X})\cap\textrm{Lip}_{\tilde{\mathsf{d}}}(\mathsf{X})\times\mathsf{X}\times\mathsf{X}$

$|K_{l}(\varphi)(x)-K_{l}(\varphi)(y)|\leq
C(\|\varphi\|\vee\|\varphi\|_{\textrm{Lip}})\tilde{\mathsf{d}}(x,y).$

  * (A4)

There exist $(C,\beta_{1})\in(0,\infty)\times(0,\infty)$ such that for any
$(l,\varphi)\in\mathbb{N}\times\mathcal{B}_{b}(\mathsf{X})$

    1. 1.

$|[\pi_{l}-\pi](\varphi)|\leq C\|\varphi\|\Delta_{l}^{\beta_{1}}.$

    2. 2.

$\sup_{x\in\mathsf{X}}|K_{l}(\varphi)(x)-K_{l-1}(\varphi)(x)|\leq
C\|\varphi\|\Delta_{l}^{\beta_{1}}.$

  * (A5)

    1. 1.

There exist $(C,\epsilon,\beta_{2})\in(0,\infty)^{3}$, such that for the
metric $\tilde{\mathsf{d}}$ as in (A(A3)) and any
$(l,n)\in\mathbb{N}\times\mathbb{N}$

$\displaystyle\mathbb{E}[\mathbb{I}_{\mathsf{B}(C,\Delta_{l}^{\beta_{2}},\tilde{\mathsf{d}})^{c}}(Z_{0,l,l-1})]$
$\displaystyle\leq$ $\displaystyle C\Delta_{l}^{\beta_{2}(2+\epsilon)},$
$\displaystyle\mathbb{E}[\mathbb{I}_{\mathsf{B}(C,\Delta_{l}^{\beta_{2}},\tilde{\mathsf{d}})^{c}\times
B(C,\Delta_{l}^{\beta_{2}},\tilde{\mathsf{d}})}(Z_{n,l,l-1},Z_{n-1,l,l-1})]$
$\displaystyle\leq$ $\displaystyle C\Delta_{l}^{\beta_{2}(2+\epsilon)}.$

    2. 2.

$\tilde{d}^{4(2+\epsilon)}\in\mathcal{L}_{\mathscr{V}\otimes\mathscr{V}}$,
with $\tilde{\mathsf{d}},\epsilon$ as in 1. and $\mathscr{V}$ as in (A(A1)).

###  3.2 Discussion of Assumptions

Let us briefly discuss each assumption in-turn. Assumption (A(A1)) is related
to the ergodicity of the underdamped Langevin sampler. The specific result we
have used, arises in the work [11], where the authors show that one can attain
$\mathcal{V}$-uniform ergodicity, for the discretized set of equations. These
conditions hold so long as $l_{*}$ is large enough, which we shall also
assume. Assumptions (A(A3))-(A(A4)) are the same as in [18], which are
discussed in detail there, except for (A(A5)) 2\. which is a simple
integrability assumption on the metric $\tilde{\mathsf{d}}$ which allows one
to remove the need for $\mathsf{U}$ to be compact in [18] and is not the case
in our context.

###  3.3 Main Result and Implications

We are now in a position to state our main result, which is that our estimator
is unbiased with finite variance.

######  Proposition 3.1.

Assume (A(A1)-(A4)). Then there exists a choice of positive probability mass
function $\mathbb{P}_{L}$, such that for the metric $\tilde{\mathsf{d}}$ in
(A(A3)) and any
$\varphi\in\mathcal{B}_{b}(\mathsf{X})\cap\textrm{\emph{Lip}}_{\tilde{\mathsf{d}}}(\mathsf{X})$,
(2.15) is an unbiased and finite variance estimator of $\pi(\varphi)$.

######  Remark 3.1.

As in [18, Theorem 2.1.], one can easily extend this result to the case of
(2.16).

Proposition 3.1 establishes that there exists a choice of positive probability
mass function $\mathbb{P}_{L}$ which ensures that (2.15) is unbiased and
finite variance, but not how to find one. That is the topic of the current
discussion and follows very much that given in [18, Section 2.4.3.]. Under our
assumptions, noting (2.7), using the argument contained in the Appendix and in
[18], the variance of (2.15) is upper-bounded by (with $C<\infty$ a generic
constant that does not depend upon $l$, but whose value may change on each
appearance)

$\displaystyle
C\sum_{l\in\mathbb{N}_{l_{*}}}\frac{\mathbb{E}\left[\xi_{l}^{2}\right]}{\mathbb{P}_{L}(l)}\leavevmode\nobreak\
\leq\leavevmode\nobreak\
C\sum_{l\in\mathbb{N}_{l_{*}}}\frac{\Delta_{l}^{2\beta}}{\mathbb{P}_{L}(l)},$
(3.1)

with $\xi_{l}:=\widehat{[\pi_{l}-\pi_{l-1}](\varphi)}_{k}$ defined in (2.13)
and $\beta=\min\left\\{\beta_{1},\tfrac{\beta_{2}}{2}\right\\}$. The expected
cost is upper-bounded by

$C\sum_{l\in\mathbb{N}_{l_{*}}}\mathbb{E}\left[\check{\tau}_{l,l-1}\right]\Delta_{l}^{-1}\mathbb{P}_{L}(l)\leavevmode\nobreak\
\leq\leavevmode\nobreak\
C\sum_{l\in\mathbb{N}_{l_{*}}}\Delta_{l}^{-1}\mathbb{P}_{L}(l),$

where it is assumed that the cost to sample from the kernel $K_{l}$ is bounded
by $C\Delta_{l}^{-1}$ and that, by A(A2), one can upper-bound the expected
stopping time $\mathbb{E}\left[\check{\tau}_{l,l-1}\right]$. If one took
$\mathbb{P}_{L}(l)\propto\Delta_{l}^{\omega}$, then if $\omega\in(1,2\beta)$
both the variance and expected cost are finite. From standard results on Euler
discretization of diffusions one might expect that $\beta_{1}=1$ as confirmed
in Figure 2, where we estimate the strong error rate for Euler–Maruyama method
considering the models in Section 4. We also plot the second moment of the
increments $\xi_{l}$ in Figure 3. If we assume that
$\mathbb{E}[\xi_{l}^{2}]=\mathcal{O}(\Delta_{l}^{2\beta})$, then Figure 3
gives an estimate of $\beta$ for each model, which is approximately $1$, and
as a result, $\beta_{2}>2$.

Figure 2: Euler–Maruyama (E–M) weak error $|[\pi_{l}-\pi](\varphi)|$ against
the time step-size $\Delta_{l}$ for the models in Section 4, with $l$ the
discretization level, $\varphi(u)=u_{1}$ (the first co-ordinate),
$\pi(\varphi)$ the reference expectation which we take to be the average of
$5200$ samples of Euler-discretized $\pi_{L}(\varphi)$ at level $L=18$, a very
high level of discretization. These plots verify that the weak error rate
$\beta_{1}$ of E–M is approximately 1 as expected.

Figure 3: Second moment of the time-averaged estimator increments
$\xi_{l}=\widehat{[\pi_{l}-\pi_{l-1}](\varphi)}_{S,T,k,m}$ against
discretization level $l$ for the different models in Section 4. The function
considered here is $\varphi(u)=u_{1}$, $u\in\mathbb{R}^{2d}$. These plots
provide a numerical estimate of the parameter $\beta$ in (3.1) of value 1.

##  4 Numerical Results

In this section we seek to illustrate and verify our theoretical results using
three numerical examples arising in statistics and physics. Our experiments
will be based on demonstrating both finite variance and unbiasedness of the
proposed estimator. We will test the methodology on Bayesian logistic
regression problem, a double-well potential and Ginzburg-Landau model.
Finally, we will compare our methodology to the unbiased Schrodinger–Föllmer
sampler for one particular example.

In all examples below we choose $\kappa$ and $\sigma$ such that
$2\kappa=\sigma^{2}$ and set $b$ in (2.4) to $b(x)=-\nabla U(x)$, where
$U=-\log\pi$ and $\pi$ is the density of interest that we aim to sample from.
The objective is to estimate, unbiasedly, the expectation $\pi(\varphi)$. We
find that the algorithm works quite well when $\sigma$ is assigned a large
value, but it may collapse for small values. In our simulations, we set
$\sigma=3$ and take $\varphi(x)=x$.

###  4.1 Bayesian Logistic Regression

We consider the binary logistic regression model in which the binary
observations $\\{Y_{i}\\}_{i=1}^{n}$ are conditionally independent Bernoulli
random variables such that $Y_{i}\in\\{0,1\\}$ and

$\mathbb{P}(Y_{i}=1|X_{i}=x_{i},\alpha)=\rho(\alpha^{T}x_{i}),$

where $\rho:\mathbb{R}\to(0,1)$ defined by $\rho(w)=e^{w}/(1+e^{w})$ is the
logistic function and $X_{i},\leavevmode\nobreak\ \alpha\in\mathbb{R}^{d}$ are
the covariates and the unknown regression coefficients, respectively. The
prior density for the parameter $\alpha$ is a multivariate normal Gaussian
given by

$pr(\alpha)=\phi(\alpha;0,\Sigma_{\alpha}),$

where $\Sigma_{\alpha}$ is defined through its inverse
$\Sigma_{\alpha}^{-1}=\frac{1}{n}(\sum_{i=1}^{n}X_{i}X_{i}^{T})$. The
covariate vectors $\\{X_{i}\\}_{i=1}^{n}$ are sampled independently from
$\mathcal{U}\\{-1,1\\}^{d}$ which then are standardized. The density of the
posterior distribution of $\alpha$ is then given by

$\pi\left(\alpha\leavevmode\nobreak\ |\leavevmode\nobreak\
\\{X_{i}=x_{i},Y_{i}=y_{i}\\}_{i=1}^{n}\right)\propto\exp\left(-U(\alpha)\right),$

where

$U(\alpha)=-\sum_{i=1}^{n}\left[y_{i}\alpha^{T}x_{i}-\log(1+\exp(\alpha^{T}x_{i}))\right]+\frac{1}{2}\alpha^{T}\Sigma_{\alpha}^{-1}\alpha,$

and

$\nabla
U(\alpha)=-\sum_{i=1}^{n}y_{i}x_{i}+\sum_{i=1}^{n}\frac{\exp(\alpha^{T}x_{i})x_{i}}{1+\exp(\alpha^{T}x_{i})}+\alpha^{T}\Sigma_{\alpha}^{-1}.$

We take $d=5$ and $n=100$. The reference expectation is the mean of $2.6\times
10^{8}$ samples computed through the unbiased MCMC method proposed in [20]
with a random-walk Metropolis-Hasting kernel.

###  4.2 Double-Well Model

For our second example we consider sampling from $\pi(x)=\exp{(-U(x))}$, where
$U$ is a double-well potential given by

$U(x)=\frac{1}{4}\|x\|_{2}^{4}-\frac{1}{2}\|x\|_{2}^{2},\qquad
x\in\mathbb{R}^{d}.$

The double-well potential is one of several quartic potentials of substantial
importance in quantum mechanics and quantum field theory [23] for the
investigation of different physical phenomena or mathematical features. The
gradient of the potential is given by

$\nabla U(x)=(\|x\|_{2}^{2}-1)\leavevmode\nobreak\ x.$

In our simulations, we take $d=100$. The true reference expectation is
$\pi(\varphi)=0$.

###  4.3 Ginzburg-Landau Model

In the final example we consider the Ginzburg-Landau (GL) model [15, 19, 25],
a model that describes a thermodynamic system that undergoes continuous phase
transitions at a temperature $T=T_{c}$, from a high-temperature symmetric
phase to a low-temperature ordered phase in which some symmetry is broken. We
denote by the random variable $\psi\in\mathbb{R}$ the order parameter, which
is assumed to be spatially dependent, i.e., $\psi=\psi(x)$,
$x\in\Omega\subseteq\mathbb{R}^{3}$, and $\nabla_{x}\psi\cdot n=0$, where $n$
is a unit vector normal to the boundary $\partial\Omega$. In the absence of
external fields, the probability of a fluctuation $\psi(x)$ is given by

$\pi(\psi)\propto\exp\left\\{-U(\psi(x))\right\\},$

with $U$ is the GL free energy functional defined via

$U(\psi(x))=\int_{\Omega}\left[\frac{1-\overline{T}}{2}\psi^{2}(x)+\frac{\gamma\overline{T}}{2}\|\nabla_{x}\psi(x)\|^{2}+\frac{\zeta\overline{T}}{4}\psi^{4}(x)\right]dx,$

$\overline{T}=T_{c}/T$ and $\gamma,\zeta>0$. We now consider the discretized
GL model on a 3d-lattice with $d=d_{0}^{3}$ sites,
$(x_{i,j,k})_{i,j,k=1}^{d_{0}}$, where $d_{0}\in\mathbb{N}$, with spacing
equal to one. Periodic boundary conditions are applied to the system where we
have $x_{11,j,k}=x_{1,j,k}$ and $x_{0,j,k}=x_{10,j,k}$ and similar situation
for the second and the third coordinates. Each site represents a random
variable $\psi_{i,j,k}\in\mathbb{R}$. Set
$\bm{\psi}:=(\psi_{i,j,k})_{i,j,k=1}^{d_{0}}\in\mathbb{R}^{d}$, then an
approximation to $\pi$ is given by

$\widetilde{\pi}(\bm{\psi})\propto\exp\left\\{-\widetilde{U}(\bm{\psi})\right\\},$

with

$\displaystyle\widetilde{U}(\bm{\psi})$
$\displaystyle=\sum_{i,j,k=1}^{d_{0}}\Big{[}\frac{1-\overline{T}}{2}\psi_{i,j,k}^{2}+\frac{\gamma\overline{T}}{2}\left((\psi_{i+1,j,k}-\psi_{i,j,k})^{2}+(\psi_{i,j+1,k}-\psi_{i,j,k})^{2}+(\psi_{i,j,k+1}-\psi_{i,j,k})^{2}\right)+\frac{\zeta\overline{T}}{4}\psi_{i,j,k}^{4}\Big{]},$

such that a forward finite difference is used to approximate $\nabla_{x}\psi$.
The functional derivative $\delta U(\psi(x^{\prime}))/\delta\psi(x)$
is given by

$\displaystyle\frac{\delta U(\psi(x^{\prime}))}{\delta\psi(x)}$
$\displaystyle=\int_{\Omega}\Big{[}(1-\overline{T})\psi(x^{\prime})\delta(x-x^{\prime})+\gamma\overline{T}\nabla_{x^{\prime}}\psi(x^{\prime})\cdot\nabla_{x^{\prime}}\delta(x-x^{\prime})+\zeta\overline{T}\psi^{3}(x^{\prime})\delta(x-x^{\prime})\Big{]}dx^{\prime}$
$\displaystyle=(1-\overline{T})\psi(x)-\gamma\overline{T}\Delta_{x}\psi(x)+\zeta\overline{T}\psi^{3}(x),$
(4.1)

where we used the following identities

$\displaystyle\frac{\delta}{\delta\psi(x)}\psi(x^{\prime})$
$\displaystyle=\delta(x-x^{\prime}),$
$\displaystyle\int_{\Omega}\nabla_{x^{\prime}}\psi(x^{\prime})\cdot\nabla_{x^{\prime}}\delta(x-x^{\prime})dx^{\prime}$
$\displaystyle=-\int_{\Omega}\Delta_{x^{\prime}}\psi(x^{\prime})\delta(x-x^{\prime})dx^{\prime}=-\Delta_{x}\psi(x).$

Here $\delta(x)$ denotes the Dirac delta function. In the last identity we
used integration by parts with the fact that $\nabla_{x}\psi\cdot n=0$ for any
normal vector $n$ on the boundary. On the lattice, an approximation to the
functional derivative in (4.1), which we denote by
$\widetilde{\nabla_{\bm{\psi}}U}$, is given by

$\displaystyle\widetilde{\nabla_{\bm{\psi}}U}$
$\displaystyle=\Big{(}(1-\overline{T})\psi_{i,j,k}+\gamma\overline{T}\vartheta_{i,j,k}+\zeta\overline{T}\psi_{i,j,k}^{3}\Big{)}_{i,j,k=1}^{d_{0}},$

and

$\vartheta_{i,j,k}=\big{\\{}6\psi_{i,j,k}-(\psi_{i+1,j,k}+\psi_{i-1,j,k}+\psi_{i,j+1,k}+\psi_{i,j-1,k}+\psi_{i,j,k+1}+\psi_{i,j,k-1})\big{\\}},$

where a second-order central finite difference is used to approximate
$\Delta_{x}\psi$. Then $b-\widetilde{\nabla_{\bm{\psi}}U}$ in (2.4). The true
reference expectation is $\pi(\varphi)=0\in\mathbb{R}^{d}$, where $d=10^{3}$.
We also set $\overline{T}=2$, $\gamma=0.1$ and $\zeta=0.5$.

###  4.4 Simulation Settings

The time-averaged, unbiased estimator

$\displaystyle\left(\widehat{\pi}(\varphi)_{S,T,k,m}\right)_{\text{avg}}:=\frac{1}{M}\sum_{i=1}^{M}\left(\widehat{\pi}(\varphi)_{S,T,k,m}\right)^{(i)},$
(4.2)

where $\left(\widehat{\pi}(\varphi)_{S,T,k,m}\right)^{(i)}$ is the estimator
presented in (2.16) in Algorithm 1 and the cost of the single-level, time-
averaged estimator

$\displaystyle\left(\widehat{\pi_{L}}(\varphi)_{T,k,m}\right)_{\text{avg}}:=\frac{1}{M}\sum_{i=1}^{M}\left(\widehat{\pi_{L}}(\varphi)_{T,k,m}\right)^{(i)},$
(4.3)

where $\left(\widehat{\pi_{L}}(\varphi)_{T,k,m}\right)^{(i)}$ is presented in
(2.11) and $L$ is the discretization level. We will compare the cost of both
estimators versus the mean-squared errors (MSEs) that are obtained by running
50 independent simulations of each estimator and are given by

$\text{MSE}_{ub}=\frac{1}{50}\sum_{j=1}^{50}\left[\left(\widehat{\pi}(\varphi)_{S,T,k,m}\right)_{\text{avg}}^{(j)}-\pi(\varphi)\right]^{2},\quad\text{and}\quad\text{MSE}_{s}=\frac{1}{50}\sum_{i=1}^{50}\left[\left(\widehat{\pi_{L}}(\varphi)_{T,k,m}\right)_{\text{avg}}^{(j)}-\pi(\varphi)\right]^{2},$

where $\pi(\varphi)$ is the reference expectation. Let $\tau$ be the stopping
time, we set $k=100$ and $m=\min\\{2k,\tau-1\\}$, however, if $\tau<k-1$, we
set $k=0.5\leavevmode\nobreak\ \tau$.

For a given $\epsilon>0$, the goal is to obtain an MSE of order
$\epsilon^{2}$. Therefore, for the single-level estimator, we set
$L=\mathcal{O}(-\log_{2}(\epsilon))$ and $M=\mathcal{O}(\epsilon^{-2})$. For
the unbiased estimator, in practice, one has to truncate the values of $L$. As
a result, we set $\mathbb{P}_{L}(l)\propto
2^{-1.5l}\mathbb{I}_{\\{l_{*},\cdots,l_{\text{max}}\\}}(l)$, where
$l_{*},l_{\text{max}}\in\mathbb{N}_{0}$ and $l_{*}<l_{\text{max}}$. With this
choice, the probability of simulating (2.4)-(2.5) at a high discretization
level is very small. In all examples, we take $l_{*}=5$, $l_{\text{max}}=12$
and $M=\mathcal{O}(\epsilon^{-2})$ as in the single-level estimator. The cost
to compute the single-level estimator is
$\mathcal{C}_{\textrm{s}}:=2^{L+1}\sum_{j=1}^{M}\tau_{L}^{(j)}=2^{L+1}\leavevmode\nobreak\
M\leavevmode\nobreak\
\max_{j}\\{\tau_{L}^{(j)}\\}\leq\mathsf{C}2^{L}\leavevmode\nobreak\
M=\mathcal{O}(\epsilon^{-3})$ where $\tau_{L}^{(j)}$ is the stopping time for
the $j$-th replicate that is assumed to be bounded by A(A2). The cost of the
unbiased estimator is
$\mathcal{C}_{\textrm{ub}}:=\sum_{j=1}^{M}\mathcal{C}^{(j)}$ where

$\displaystyle\mathcal{C}^{(j)}=\left\\{\begin{array}[]{lr}\tau_{L_{j}}\leavevmode\nobreak\
2^{L_{j}+1}&\text{ if }L_{j}=l_{*},\\\
\check{\tau}_{L_{j},L_{j-1}}\leavevmode\nobreak\
(2^{L_{j}+1}+2^{L_{j}})&\text{ if }L_{j}>l_{*},\\\ \end{array}\right.\qquad
L_{j}\sim\mathbb{P}_{L}.$

######  Remark 4.1.

We note that even though the rate for the single-level estimator should be
$\mathcal{O}(\epsilon^{-3})$, what we observe numerically is different. As we
will discuss in the next subsection, the rates obtained are worse, where it
seems that we attain a rate of $\mathcal{O}(\epsilon^{-(3+\delta)})$, for
$\delta>0$. This is for each individual numerical example and we expect that
this is associated to the drift coefficient in the dynamics, which are in-turn
affected by the target. The regularity of the drift will influence the
practical rate of convergence of the discretization scheme and the time
convergence of the dynamics.

By Remark 3.1, the expected cost of each replicate is bounded by a constant
$\mathsf{C}$ and hence the overall expected cost is
$\overline{\mathcal{C}_{\textrm{ub}}}\leq\mathsf{C}M=\mathcal{O}(\epsilon^{-2})$.
Notice that in fact the cost of sampling from $\check{Q}_{L}$ is $2^{L}$ and
the cost of sampling from $\check{P}_{L}$ is $2^{L}-1$ plus the cost of
sampling from a maximal coupling, which we take it to be one, and therefore
the cost of sampling from
$\check{K}_{L}=\alpha\check{Q}_{L}+(1-\alpha)\check{P}_{L}$ is $2^{L+1}$.

We will also compare the unbiased SFS algorithm [35] with our unbiased ULD
sampler Algorithm 1. For a fair comparison, we consider comparing the average
machine time (in seconds) that one independent realization of the 50
realizations takes on a workstation of 52 cores. The unbiased SFS estimator is
based on the SF diffusion process,

$dX_{t}=b(X_{t},t)dt+dW_{t},\quad X_{0}=0,\leavevmode\nobreak\ t\in[0,1],$

where the drift term has the specific form

$b(x,t)=\nabla\log\mathbb{E}[f(x+W_{1-t})]=\frac{\mathbb{E}_{\phi}[\nabla
f(x+\sqrt{1-t}Z)]}{\mathbb{E}_{\phi}[f(x+\sqrt{1-t}Z)]}=\frac{1}{\sqrt{1-t}}\frac{\mathbb{E}_{\phi}[Zf(x+\sqrt{1-t}Z)]}{\mathbb{E}_{\phi}[f(x+\sqrt{1-t}Z)]},$

and $f$ corresponds to the analogous of a likelihood function for a standard
Gaussian prior, i.e. for $z\in\mathbb{R}^{d}$, $f(z)=\pi(z)/\phi(z),$ with
$\phi(z)$ the standard $d-$dimensional Gaussian density. For
$l\in\\{0,1,\dots\\}$ and $k\in\\{0,1,\dots,\Delta_{l}^{-1}-1\\}$, we discrete
the above SDE as

$\widetilde{X}_{(k+1)\Delta_{l}}^{l,N}=\widetilde{X}_{k\Delta_{l}}^{l,N}+\hat{b}(\widetilde{X}_{k\Delta_{l}}^{l,N},k\Delta_{l})\Delta_{l}+W_{(k+1)\Delta_{l}}-W_{k\Delta_{l}},$

where, for $x\in\mathbb{R}^{d}$

$\hat{b}(x,k\Delta_{l})=\frac{\frac{1}{N}\sum_{i=1}^{N}\nabla
f(x+\sqrt{1-k\Delta_{l}}Z^{i})}{\frac{1}{N}\sum_{i=1}^{N}f(x+\sqrt{1-k\Delta_{l}}Z^{i})},$

and, for $i\in\\{1,\dots,N\\}$,
$Z^{i}\stackrel{{\scriptstyle\textrm{i.i.d.}}}{{\sim}}\mathcal{N}_{d}(0,I)$.
The idea of the unbiased SFS is based upon double-randomization techniques
(e.g. [21, 22]) where both $N$ and $l$ are chosen at random in a specific
manner.

###  4.5 Results and Discussion

In Figure 4 (a)–(c), we plot the cost of the estimators in (4.2)-(4.3) against
their MSE for the models above. The MSE of the proposed unbiased estimator
presented in Algorithm 1 decays at the optimal rate of
$1/\mathcal{C}_{\text{ub}}=\mathcal{O}(\epsilon^{2})$ as shown in the plots,
which is as expected, and outperforms that of the single-level estimator given
in (4.3). As eluded to, from Remark 4.1 the actual rate which we obtain for
each model is around $\mathcal{O}(\epsilon^{-(3+\delta)})$, for $\delta>0$,
which differs for each model. This is especially the case for the GL model
which, compared to the other models, has a more complicated drift term.
Therefore we believe this constitutes to a modified rate for the single-level
estimator. Further investigation is needed to verify this, but this is beyond
the scope of this work and we leave it for future work. Nonetheless, as
stated, the unbiased estimator outperforms the single estimator for each model
example.

In Figure 4 (d), the cost of our unbiased estimator and that of the SFS,
measured in seconds, is plotted against the MSE of each method. As we can see,
the new methodology beats the unbiased SFS and in fact the MSE decays faster
than expected by the theory. This suggesting to attain similar values of MSE,
our unbiased estimator takes considerably less time than that using the SFS,
as presented in [35].

(a)  (b)

(c)  (d)

Figure 4: In (a)–(c), we plot the cost versus MSE for the single-level, time-
averaged estimator in (2.11) and the time-averaged unbiased estimator
presented in (2.16) in Algorithm 1. In (d), we plot the machine time versus
MSE for the unbiased estimator in Algorithm 1 and the unbiased SFS estimator
in [35] for the Bayesian logistic model.

##  5 Conclusion

Our motivation from this work was to provide a new unbiased estimator, which
is based on the discretized underdamped Langevin dynamics (ULD) (2.3)-(2.4).
The ULD has sparked recent interest in both the statistics and machine
learning community, and as a result we wanted to see if such a dynamics could
be used in the context of unbiased estimation. We introduced a new methodology
based on ULD, which was based on the double-randomization schemes used for
unbiased estimation. Subsequently we proved that our new estimator is unbiased
with finite variance under suitable assumptions. To verify our theory, we
implemented our methodology on a range of interesting model problems, such as
the stochastic Ginzburg-Landau model, the double-well model and a Bayesian
logistic regression problem. We also justified such an estimator by comparing
it to our known estimators, developed in a similar manner, such as the
unbiased Schrödinger–Föllmer sampler presented in [35].

From this work, there are a number of research directions one can consider. A
rather obvious one would be to use the current methodology aimed at unbiased
estimation for the both the score function and the Hessian [6, 17]. This has
already been considered previously, but in the context of our methodology
which may prove to be more useful. Another direction is to extend the
ergodicity result presented in [11], where one does not have the requirement
that $l_{*}$ has to be large enough to ensure $\mathcal{V}$-uniform
ergodicity. This may prove to be challenging, but from a theoretical
perspective would be of interest, especially if one can attain geometric
ergodicity. One could also consider the mean-field ULD, which has shown to be
promising for deep neural architectures [24], in terms of the trainability of
two-layer neural networks [8, 26]. Related to the above point, another
extension could be the perturbed ULD [10], which has demonstrated improvements
for sampling even more complex probability measures.

#### Acknowledgements

All three authors were supported by KAUST baseline funding.

##  Appendix A Proof

The proof of Proposition Proposition 3.1 is virtually identical to that of
[18, Theorem 2.1]. There is only place where the proof has to be modified and
we give the result below. Below $C$ is a generic constant which does not
depend upon $n,l$ and whose value changes upon each appearance. The
expectation operator $\mathbb{E}$ relates to law which generates the simulated
process used to compute (2.15).

######  Lemma A.1.

Assume (A(A1),(A4)). Then there exists a $C\in(0,\infty)$ such that for any
$(l,n)\in\mathbb{N}_{l_{*}}\times\mathbb{Z}^{+}$ we have:

$\displaystyle\mathbb{E}[\mathbb{I}_{\mathsf{B}(C,\Delta_{l}^{\beta_{2}},\tilde{\mathsf{d}})^{c}}(Z_{n,l,l-1})]$
$\displaystyle\leq$ $\displaystyle C(n+1)\Delta_{l}^{\beta_{2}(2+\epsilon)},$
$\displaystyle\max\left\\{\mathbb{E}[\tilde{\mathsf{d}}(U_{n,l},U_{n,l-1})^{2+\epsilon}],\mathbb{E}[\tilde{\mathsf{d}}(\tilde{U}_{n,l},\tilde{U}_{n,l-1})^{2+\epsilon}]\right\\}$
$\displaystyle\leq$ $\displaystyle
C(n+1)\Delta_{l}^{\beta_{2}(1+\tfrac{\epsilon}{2})},$

where $\beta_{2}$ and $\epsilon$ are as (A(A4)).

###### Proof.

The first inequality is proved in an identical manner to [18, Lemma A.3.], so
we consider only the second inequality. This latter inequality is shown for
$\mathbb{E}[\tilde{\mathsf{d}}(U_{n,l},U_{n,l-1})^{2+\epsilon}]$, only, as the
argument for the other term is the same up to changes in notation. Let
$n\in\mathbb{N}$, then we have

$\mathbb{E}[\tilde{\mathsf{d}}(U_{n,l},U_{n,l-1})^{2+\epsilon}]=\mathbb{E}[\tilde{\mathsf{d}}(U_{n,l},U_{n,l-1})^{2+\epsilon}\mathbb{I}_{\mathsf{B}(C,\Delta_{l}^{\beta_{2}},\tilde{\mathsf{d}})^{c}}(Z_{n,l,l-1})]+\mathbb{E}[\tilde{\mathsf{d}}(U_{n,l},U_{n,l-1})^{2+\epsilon}\mathbb{I}_{\mathsf{B}(C,\Delta_{l}^{\beta_{2}},\tilde{\mathsf{d}})}(Z_{n,l,l-1})].$

For the second term on the R.H.S. as
$Z_{n,l,l-1}\in\mathsf{B}(C,\Delta_{l}^{\beta_{2}},\tilde{\mathsf{d}})$, we
have $\tilde{\mathsf{d}}(U_{n,l},U_{n,l-1})^{2+\epsilon}\leq
C\Delta_{l}^{\beta_{2}(2+\epsilon)}$, so we focus on the first term on the
R.H.S.. Applying Cauchy-Schwarz and the first-statement of the Lemma, we have

$\mathbb{E}[\tilde{\mathsf{d}}(U_{n,l},U_{n,l-1})^{2+\epsilon}\mathbb{I}_{\mathsf{B}(C,\Delta_{l}^{\beta_{2}},\tilde{\mathsf{d}})^{c}}(Z_{n,l,l-1})]\leq\mathbb{E}[\tilde{\mathsf{d}}(U_{n,l},U_{n,l-1})^{2(2+\epsilon)}]^{1/2}C(n+1)\Delta_{l}^{\beta_{2}(1+\tfrac{\epsilon}{2})}.$

Then using
$\tilde{\mathsf{d}}(U_{n,l},U_{n,l-1})^{4(2+\epsilon)}\in\mathcal{L}_{\mathscr{V}\otimes\mathscr{V}}$,
we get

$\mathbb{E}[\tilde{\mathsf{d}}(U_{n,l},U_{n,l-1})^{2+\epsilon}\mathbb{I}_{\mathsf{B}(C,\Delta_{l}^{\beta_{2}},\tilde{\mathsf{d}})^{c}}(Z_{n,l,l-1})]\leq\mathbb{E}[\\{\mathscr{V}(U_{n,l})\mathscr{V}(U_{n,l-1})\\}^{\frac{1}{2}}]^{1/2}C(n+1)\Delta_{l}^{\beta_{2}(1+\tfrac{\epsilon}{2})}.$

Then applying Cauchy-Schwarz and using (A(A1)) 1. we have

$\mathbb{E}[\tilde{\mathsf{d}}(U_{n,l},U_{n,l-1})^{2+\epsilon}\mathbb{I}_{\mathsf{B}(C,\Delta_{l}^{\beta_{2}},\tilde{\mathsf{d}})^{c}}(Z_{n,l,l-1})]\leq
C(n+1)\Delta_{l}^{\beta_{2}(1+\tfrac{\epsilon}{2})},$

from which the proof can be completed. ∎

##  Appendix B Algorithms

In this Appendix we provide each algorithm required for our unbiased
estimator. This is related to the simulation from the associated kernels
$\check{Q}_{l}$, $\check{Q}_{l,l-1}$, $\check{P}_{l}$, $\check{P}_{l,l-1}$,
which are discussed in Algorithms 2 to 6. For the earlier algorithms, fuller
details can be found in [39]. Recall that our objective is to construct the
coupling, or coupled kernels, $\check{K}_{l}$ and $\check{K}_{l,l-1}$, which
can be decomposed with the abovely stated kernels where

$\displaystyle\check{K}_{l}$ $\displaystyle=$
$\displaystyle\alpha\check{Q}_{l}+(1-\alpha)\check{P}_{l},$
$\displaystyle\check{K}_{l,l-1}$ $\displaystyle=$
$\displaystyle\mathbb{I}_{D^{2}}(u_{l},\tilde{u}_{l},u_{l-1},\tilde{u}_{l-1})\check{Q}_{l,l-1}+\mathbb{I}_{(D^{2})^{c}}(u_{l},\tilde{u}_{l},u_{l-1},\tilde{u}_{l-1})[\alpha\check{Q}_{l,l-1}+(1-\alpha)\check{P}_{l,l-1}],$

where $\alpha\in(0,1)$. Finally we require the coupling
$\overline{K}_{l,l-1}$, which is required for the initialization before
sampling from $\check{K}_{l,l-1}$, and is presented in Algorithm 4.

Algorithm 2 Sampling from $\check{Q}_{l}$.

  1. 1.

Input: $l\in\mathbb{N}_{l_{*}}$,
$(u,\tilde{u})=\left((x_{0},v_{0}),(\tilde{x}_{0},\tilde{v}_{0})\right)\in\mathsf{U}^{2}$.

  2. 2.

Sample one sequence $(\Gamma_{0,l},\dots,\Gamma_{\Delta_{l}^{-1}-1,l})$ and
one sequence
$B_{\Delta_{l}},B_{2\Delta_{l}}-B_{\Delta_{l}},\dots,B_{1}-B_{1-\Delta_{l}}$.

  3. 3.

Run the recursion (2.4)-(2.5) with
$\\{\Gamma_{k,l}\\}_{k\in\\{0.\dots,\Delta_{l}^{-1}-1\\}}$ and
$\\{B_{(k+1)\Delta_{l}}-B_{k\Delta_{l}}\\}_{k\in\\{0,\dots,\Delta_{l}^{-1}-1\\}}$
simulated at the previous step, up-to time 1.

  4. 4.

Output: $u^{\prime}=(x_{1},v_{1})$ and
$\tilde{u}^{\prime}=(\tilde{x}_{1},\tilde{v}_{1})$ as generated in step 3..

Algorithm 3 Sampling from $\check{P}_{l}$.

  1. 1.

Input: $l\in\mathbb{N}_{l_{*}}$,
$(u,\tilde{u})=\left((x_{0},v_{0}),(\tilde{x}_{0},\tilde{v}_{0})\right)\in\mathsf{U}^{2}$.

  2. 2.

Sample one sequence $(\Gamma_{0,l},\dots,\Gamma_{\Delta_{l}^{-1}-2,l})$ and
one sequence
$B_{\Delta_{l}},B_{2\Delta_{l}}-B_{\Delta_{l}},\dots,B_{1-\Delta_{l}}-B_{1-2\Delta_{l}}$.

  3. 3.

Run the recursion (2.4)-(2.5) with the sequence
$\\{\Gamma_{k,l}\\}_{k\in\\{0.\dots,\Delta_{l}^{-1}-2\\}}$
and
$\\{B_{(k+1)\Delta_{l}}-B_{k\Delta_{l}}\\}_{k\in\\{0,\dots,\Delta_{l}^{-1}-2\\}}$
simulated at the previous step, up-to time $1-\Delta_{l}$.

  4. 4.

Sample
$\left((X_{1},V_{1}),(\tilde{X}_{1},\tilde{V}_{1})\right)|\left((x_{1-\Delta_{l}},v_{1-\Delta_{l}}),(\tilde{x}_{1-\Delta_{l}},\tilde{v}_{1-\Delta_{l}})\right)$
from a maximal coupling (e.g. the reflection maximal coupling described in
[20]) of $\displaystyle
p_{l}(x_{1},v_{1}|x_{1-\Delta_{l}},v_{1-\Delta_{l}})\quad\textrm{and}\quad
p_{l}(\tilde{x}_{1},\tilde{v}_{1}|\tilde{x}_{1-\Delta_{l}},\tilde{v}_{1-\Delta_{l}}),$
where the density $p_{l}$ is that of a $2d-$dimensional Gaussian distribution
whose moments can be determined from one step of (2.4)-(2.5).

  5. 5.

Output: $u^{\prime}=(x_{1},v_{1})$ and
$\tilde{u}^{\prime}=(\tilde{x}_{1},\tilde{v}_{1})$ as generated in step 4..

Algorithm 4 Sampling from $\overline{K}_{l,l-1}$.

  1. 1.

Input: $l\in\\{l_{*}+1,l_{*}+2,\dots\\}$, and
$(u_{l},u_{l-1})=\left((x_{0,l},v_{0,l}),(x_{0,l-1},v_{0,l-1})\right)\in\mathsf{U}^{2}$.

  2. 2.

Sample one sequence $(\Gamma_{0,l},\dots,\Gamma_{\Delta_{l}^{-1}-1,l})$ and
one sequence
$B_{\Delta_{l}},B_{2\Delta_{l}}-B_{\Delta_{l}},\dots,B_{1}-B_{1-\Delta_{l}}$.
Concatenate the Brownian motion to obtain
$B_{\Delta_{l-1}},B_{2\Delta_{l-1}}-B_{\Delta_{l-1}l},\dots,B_{1}-B_{1-\Delta_{l-1}}$
and similarly for $(\Gamma_{0,l},\dots,\Gamma_{\Delta_{l}^{-1}-1,l})$ to
obtain $(\Gamma_{0,l-1},\dots,\Gamma_{\Delta_{l-1}^{-1}-1,l-1})$.

  3. 3.

For $s\in\\{l,l-1\\}$: run the recursion (2.4)-(2.5) with the sequence
$\\{\Gamma_{k,s}\\}_{k\in\\{0.\dots,\Delta_{s}^{-1}-1\\}}$
and
$\\{B_{(k+1)\Delta_{s}}-B_{k\Delta_{s}}\\}_{k\in\\{0,\dots,\Delta_{s}^{-1}-1\\}}$
simulated at the previous step, up-to time 1.

  4. 4.

Output: $(u_{l},u_{l-1})=\left((x_{1,l},v_{1,l}),(x_{1,l-1},v_{1,l-1})\right)$
as generated in step 3..

Algorithm 5 Sampling from $\check{Q}_{l,l-1}$.

  1. 1.

Input: $l\in\\{l_{*}+1,l_{*}+2,\dots\\}$, and for $s\in\\{l,l-1\\}$,
$(u_{s},\tilde{u}_{s})=\left((x_{0,s},v_{0,s}),(\tilde{x}_{0,s},\tilde{v}_{0,s})\right)\in\mathsf{U}^{2}$.

  2. 2.

Sample one sequence $(\Gamma_{0,l},\dots,\Gamma_{\Delta_{l}^{-1}-1,l})$ and
one sequence
$B_{\Delta_{l}},B_{2\Delta_{l}}-B_{\Delta_{l}},\dots,B_{1}-B_{1-\Delta_{l}}$.
Concatenate the Brownian motion to obtain
$B_{\Delta_{l-1}},B_{2\Delta_{l-1}}-B_{\Delta_{l-1}l},\dots,B_{1}-B_{1-\Delta_{l-1}}$
and similarly for $(\Gamma_{0,l},\dots,\Gamma_{\Delta_{l}^{-1}-1,l})$ to
obtain $(\Gamma_{0,l-1},\dots,\Gamma_{\Delta_{l-1}^{-1}-1,l-1})$.

  3. 3.

For $s\in\\{l,l-1\\}$: run the recursion (2.4)-(2.5) with the sequence
$\\{\Gamma_{k,s}\\}_{k\in\\{0.\dots,\Delta_{s}^{-1}-1\\}}$
and
$\\{B_{(k+1)\Delta_{s}}-B_{k\Delta_{s}}\\}_{k\in\\{0,\dots,\Delta_{s}^{-1}-1\\}}$
simulated at the previous step, up-to time 1.

  4. 4.

Output: for $s\in\\{l,l-1\\}$,
$(u_{s},\tilde{u}_{s})=\left((x_{1,s},v_{1,s}),(\tilde{x}_{1,s},\tilde{v}_{1,s})\right)$
as generated in step 3..

Algorithm 6 Sampling from $\check{P}_{l,l-1}$.

  1. 1.

Input: $l\in\\{l_{*}+1,l_{*}+2,\dots\\}$, and for $s\in\\{l,l-1\\}$,
$(u_{s},\tilde{u}_{s})=\left((x_{0,s},v_{0,s}),(\tilde{x}_{0,s},\tilde{v}_{0,s})\right)\in\mathsf{U}^{2}$.

  2. 2.

Sample one sequence $(\Gamma_{0,l},\dots,\Gamma_{\Delta_{l}^{-1}-2,l})$ and
one sequence
$B_{\Delta_{l}},B_{2\Delta_{l}}-B_{\Delta_{l}},\dots,B_{1-\Delta_{l}}-B_{1-2\Delta_{l}}$.
Concatenate the Brownian motion to obtain
$B_{\Delta_{l-1}},B_{2\Delta_{l-1}}-B_{\Delta_{l-1}l},\dots,B_{1-\Delta_{l-1}}-B_{1-2\Delta_{l-1}}$
and similarly for $(\Gamma_{0,l},\dots,\Gamma_{\Delta_{l}^{-1}-2,l})$ to
obtain $(\Gamma_{0,l-1},\dots,\Gamma_{\Delta_{l-1}^{-1}-2,l-1})$.

  3. 3.

For $s\in\\{l,l-1\\}$: run the recursion (2.4)-(2.5) with the sequence
$\\{\Gamma_{k,s}\\}_{k\in\\{0.\dots,\Delta_{s}^{-1}-1\\}}$
and
$\\{B_{(k+1)\Delta_{s}}-B_{k\Delta_{s}}\\}_{k\in\\{0,\dots,\Delta_{s}^{-1}-1\\}}$
simulated at the previous step, up-to time $1-\Delta_{s}$.

  4. 4.

Sample

$\displaystyle\left((X_{1,l},V_{1,l}),(\tilde{X}_{1,l},\tilde{V}_{1,l})\right),\left((X_{1,l-1},V_{1,l-1}),(\tilde{X}_{1,l-1},\tilde{V}_{1,l-1})\right)\Big{|}\left((x_{1-\Delta_{l},l},v_{1-\Delta_{l},l}),(\tilde{x}_{1-\Delta_{l},l},\tilde{v}_{1-\Delta_{l},l})\right),$
$\displaystyle((x_{1-\Delta_{l-1},l-1},v_{1-\Delta_{l-1},l-1}),$
$\displaystyle(\tilde{x}_{1-\Delta_{l-1},l-1},\tilde{v}_{1-\Delta_{l-1},l-1}))$

from the synchronous pairwise reflection maximal coupling [4] of

$p_{l}(x_{1,l},v_{1,l}|x_{1-\Delta_{l},l},v_{1-\Delta_{l},l}),p_{l}(\tilde{x}_{1,l},\tilde{v}_{1,l}|\tilde{x}_{1-\Delta_{l},l},\tilde{v}_{1-\Delta_{l},l}),\quad\textrm{and}$
$p_{l-1}(x_{1,l-1},v_{1,l-1}|x_{1-1\Delta_{l-1},l-1},v_{1-1\Delta_{l-1},l-1}),p_{l-1}(\tilde{x}_{1,l-1},\tilde{v}_{1,l-1}|\tilde{x}_{1-\Delta_{l-1},l-1},\tilde{v}_{1-\Delta_{l-1},l-1}),$

  5. 5.

Output: for $s\in\\{l,l-1\\}$,
$(u_{s},\tilde{u}_{s})=\left((x_{1,s},v_{1,s}),(\tilde{x}_{1,s},\tilde{v}_{1,s})\right)$
as generated in step 4..

## References

  * [1] Benveniste, A., Metivier, N. & Priouret, P. (1990). _Adaptive Algorithms and Stochastic Approximations. Applications of Mathematics_. Spinger: New York.
  * [2] Beskos, A., Papaspiliopoulos, O., Roberts, G. O. & Fearnhead, P. (2006) Exact and computationally efficient likelihood-based estimation for discretely observed diffusion processes. J. R. Stat. Soc. Ser. B Stat. Methodol., 68, pp. 333–382.
  * [3] Beskos, A. & Roberts, G. (2005). Exact simulation of diffusions, Ann. Appl. Probab., 15, 422–2444.
  * [4] Bou-Rabee, N., Eberle, A., & Zimmer, R. (2020). Coupling and convergence for Hamiltonian Monte Carlo. Ann. Appl. Probab., 30, 1209–1250.
  * [5] Chada, N. K., Franks, J., Jasra, A., Law, K. J. H. & Vihola, M. (2021) Unbiased estimation of discretely observed hidden Markov models. _SIAM/ASA Journal on Uncertainty Quantification_ , 9(2), 763-787.
  * [6] Chada, N. K., Jasra, A. & Yu, F. (2022). Unbiased estimation of the Hessian for partially observed diffusions. _Proceedings of the Royal Society A.,_ (Accepted).
  * [7] Cheng, X., Chatterji, N. S. Bartlett, P. L., & Jordan, M. I. (2018). Underdamped Langevin MCMC: A non-asymptotic analysis. _Proceedings of Machine Learning Research_ , 75, 1–24.
  * [8] Chizat, L. & Bach F. (2018). On the global convergence of gradient descent for overparameterized models using optimal transport. Advances in neural information processing systems, 3040–3050.
  * [9] Dalalyan, A. S. (2017). Theoretical guarantees for approximate sampling from smooth and logconcave densities. J. R. Statist. Soc. Ser. B, 79 651–676.
  * [10] Duncan, A., Nusken, N. & Pavliotis, G. A. (2017). Using perturbed underdamped Langevin dynamics to efficiently sample from probability distributions. J Stat Phys., 169, 1098–1131.
  * [11] Durmus, A., Enfroy, A. Moulines, E., & Stoltz, G. (2021). Uniform minorization condition and convergence bounds for discretizations of kinetic Langevin dynamics. arXiv preprint.
  * [12] Durmus, A. & Moulines, E. (2016) Sampling from strongly log-concave distributions with the Unadjusted Langevin Algorithm. arxiv preprint.
  * [13] Frenkel, D. & Smit, B. (2002). _Understanding Molecular Simulation, From Algorithms to Applications_. Academic Press: New York.
  * [14] Gao, X., Gurbuzbalaban, M. & Lingjiong, Z. (2020). Breaking reversibility accelerates langevin dynamics for non-convex optimization. NIPS’20: Proceedings of the 34th International Conference on Neural Information Processing Systems, 498, 17850–17862.
  * [15] Ginzburg, V.L. & Landau, L.D. (1950). Zh. Eksp. Teor. Fiz. 20, 1064. English translation in: L. D. Landau, Collected papers (Oxford: Pergamon Press, 1965) p. 546
  * [16] Glynn, P. W. & Rhee, C. H. (2014). Exact estimation for Markov chain equilibrium expectations. _J. Appl. Probab._ , 51, 377–389.
  * [17] Heng, J., Houssineau, J. & Jasra, A. (2021). On unbiased score estimation for partially observed diffusions. arXiv preprint.
  * [18] Heng, J., Jasra, A. Law, K. J. H. & Tarakanov, A. (2021). On unbiased estimation for discretized models. arXiv preprint.
  * [19] Hohenberg, P.C. & Krekhov, A.P. (2015). An introduction to the Ginzburg-Landau theory of phase transitions and nonequilibrium patterns. _Physics Reports_ , 572, 1–42.
  * [20] Jacob, P., O’ Leary, J. & Atchadé, Y. (2020). Unbiased Markov chain Monte Carlo with couplings (with discussion). _J. R. Statist. Soc. Ser. B_ , 82, 543–600.
  * [21] Jasra, A., Law, K. J. H. & Lu, D. (2021). Unbiased estimation of the gradient of the log-likelihood in inverse problems. _Stat. Comp._ 31, 21.
  * [22] Jasra, A., Law, K. J. H. & Yu, F. (2022). Unbiased filtering of a class of partially observed diffusions. _Adv. Appl. Probab._ (to appear).
  * [23] Jelic, V. & and Marsiglio, F. (2012) The double-well potential in quantum mechanics: a simple, numerically exact formulation. _Eur. J. Phys_. 33 1651.
  * [24] Kazeykina, A., Ren, Z. Tan, X. & Yang, J. (2020) Ergodicity of the underdamped mean-field Langevin dynamics.arXiv preprint.
  * [25] Livingstone, S., Faulkner, M. F. & Roberts, G. O. (2019). Kinetic energy choice in Hamiltonian/hybrid Monte Carlo. _Biometrika_ , 106, 2, 303–319.
  * [26] Mei, S., A. Montanari, A, & P.-M. Nguyen P.-M. (2018), A mean field view of the landscape of two-layer neural networks, Proceedings of the National Academy of Sciences, 115, 7665–7671.
  * [27] McLeish, D. (2011). A general method for debiasing a Monte Carlo estimator. _Monte Carlo Meth. Appl._ , 17, 301–315.
  * [28] Middleton, L., Deligiannidis, G., Doucet, A. &Jacob, P. (2019). Unbiased smoothing using particle independent Metropolis-Hastings. In _Proc. Mach. Learn. Res._ , 89.
  * [29] Jacob, R. M. (2011). MCMC using Hamiltonian dynamics. _Handbook of Markov Chain Monte Carlo._ CRC Press.
  * [30] Pavliotis, G. A. (2014). Stochastic Processes and Applications. Springer.
  * [31] Propp, J. G., & Wilson, D. B. (1996). Exact sampling with coupled Markov chains and applications to statistical mechanics. _Rand. Struct. Algs._ , 9, 223–252.
  * [32] Rhee, C. H. (2013). Unbiased estimation with biased samples. Ph.D. thesis, Stanford University.
  * [33] Rhee, C. H. & Glynn, P. (2015). Unbiased estimation with square root convergence for SDE models. _Op. Res._ , 63, 1026–1043.
  * [34] Robert, C. & Casella, G. (2004). _Monte Carlo Statistical Methods_. Springer: New York.
  * [35] Ruzayqat, H., Beskos, A., Crisan, D., Jasra, A. & Kantas, N. (2022). Unbiased Estimation using a Class of Diffusion Processes. arXiv preprint.
  * [36] Ruzayqat, H. & Jasra A. (2022). Unbiased Parameter Inference for a Class of Partially Observed Lévy-Process Models. _Found. Data Sci._ , 4(2), 299–322.
  * [37] Ruzayqat, H. & Jasra A. (2020) Unbiased estimation of the solution to Zakai’s equation. _Monte Carlo Meth. Appl_., 26(2), 113–129.
  * [38] simsekli, U., Lingjiong, Z., Why Teh, Y. & Gurbuzbalaban, M. (2020) Fractional underdamped Langevin dynamics: retargeting SGD with momentum under heavy-tailed gradient noise. Proceedings of the 37th International Conference on Machine Learning, 832, 8970–8980.
  * [39] Thorisson, H. (2000). _Coupling, Stationarity and Regeneration_. Springer: New York.
  * [40] van den Boom, W., Jasra, A., De Iorio, M., & Beskos, A. (2022). Unbiased approximation of posteriors via coupled particle Markov chain Monte Carlo. _Stat. Comp._ (to appear).