$\displaystyle+\int_{0}^{h}(h-s)M\|\boldsymbol{V}_{s}\|_{\mathbb{L}_{2}}^{\boldsymbol{z}_{0}}\mathrm{d}s+\gamma\int_{0}^{h}(h-s)M(\|\boldsymbol{X}_{s}\|_{\mathbb{L}_{2}}^{\boldsymbol{z}_{0}}+\mathcal{D})\mathrm{d}s$
$\displaystyle\leq(h^{2}M/2)\bigg{[}(\gamma+\sqrt{T}/2)(|\boldsymbol{x}_{0}|+\mathcal{D})+|\boldsymbol{v}_{0}|+\sqrt{d}+(1+\gamma)\underset{0\leq
s\leq
h}{\sup}\|\boldsymbol{Z}_{s}\|_{\mathbb{L}_{2}}^{\boldsymbol{z}_{0}}\bigg{]}$
$\displaystyle\leq(M/2)\left[(2\gamma+\sqrt{T})(1+\mathcal{D})+2(1+\sqrt{d})+(1+\gamma)\sqrt{A}\right](1+|\boldsymbol{z}_{0}|^{2})^{1/2}h^{2}$
$\displaystyle\|\boldsymbol{X}_{h}-\boldsymbol{x}_{1}-\boldsymbol{G}_{h}^{\prime}\|_{\mathbb{L}_{2}}^{\boldsymbol{z}_{0}}$
$\displaystyle\leq(\gamma
h^{2}/2)|\boldsymbol{v}_{0}|+\int_{0}^{h}(h-s)\left(M|\boldsymbol{x}_{0}|+M\|\boldsymbol{X}_{s}\|_{\mathbb{L}_{2}}^{\boldsymbol{z}_{0}}+\gamma\|\boldsymbol{V}_{s}\|_{\mathbb{L}_{2}}^{\boldsymbol{z}_{0}}\right)\mathrm{d}s$
$\displaystyle\leq(h^{2}/2)\left[\gamma|\boldsymbol{v}_{0}|+M|\boldsymbol{x}_{0}|+(M+\gamma)\underset{0\leq
s\leq
h}{\sup}\|\boldsymbol{Z}_{s}\|_{\mathbb{L}_{2}}^{\boldsymbol{z}_{0}}\right]$
$\displaystyle\leq(1/2)\left[(\gamma+M)(1+A^{1/2})\right](1+|\boldsymbol{z}_{0}|^{2})^{1/2}h^{2}$

As a result, for $B\triangleq
2M(\gamma+\sqrt{T}/2)(1+\mathcal{D})+2M(1+\sqrt{d})+(\gamma+\gamma
M+M)A^{1/2}$ we obtain

$\|\boldsymbol{V}_{h}-\boldsymbol{v}_{1}\|_{\mathbb{L}_{2}}^{\boldsymbol{z}_{0}}+\|\boldsymbol{X}_{h}-\boldsymbol{x}_{1}-\boldsymbol{G}_{h}^{\prime}\|_{\mathbb{L}_{2}}^{\boldsymbol{z}_{0}}\leq
B(1+|\boldsymbol{z}_{0}|^{2})^{1/2}h^{2}.$

Plugging this bound into the decompositions (41) and (42), we get for any
$h\in(0,T]$

$\displaystyle\left|\mathbb{E}^{\boldsymbol{z}_{0}}[\boldsymbol{Z}_{h}-\boldsymbol{z}_{1}]\right|$
$\displaystyle\leq B(1+|\boldsymbol{z}_{0}|^{2})^{1/2}h^{2}$
$\displaystyle\left(\mathbb{E}^{\boldsymbol{z}_{0}}[|\boldsymbol{Z}_{h}-\boldsymbol{z}_{1}|^{2}]\right)^{1/2}$
$\displaystyle\leq(\sqrt{2T}B+\sqrt{8d\gamma})(1+|\boldsymbol{z}_{0}|^{2})^{1/2}h^{3/2}$

This yields the claim of Lemma 3 for $C=B\vee(\sqrt{2T}B+\sqrt{8d\gamma})$.

###  B.2 Proof of Proposition 4

We introduce $\Gamma(B)\triangleq\mathbb{P}(\boldsymbol{\xi}\in B)$, where
$\boldsymbol{\xi}\sim\mathcal{N}_{d}(\boldsymbol{0}_{d},\mathbf{I}_{d})$,
defined for any Borel set $B$ of $\mathbb{R}^{d}$. We also define the Gibbs
update, corresponding to the conditional distribution of $\mu$ given
$\boldsymbol{x}_{0}$.

$\mathbf{G}(\boldsymbol{z}_{0:L},\mathrm{d}\boldsymbol{z}_{0:L}^{\prime})\triangleq\delta_{\boldsymbol{x}_{0}}(\mathrm{d}\boldsymbol{x}_{0}^{\prime})\Gamma(\mathrm{d}\boldsymbol{v}_{0}^{\prime})\prod_{i=1}^{L}\mathbf{Q}_{h,\gamma}(\boldsymbol{z}_{i-1},\mathrm{d}\boldsymbol{z}_{i})$
(43)

The Gibbs kernel $\mathbf{G}$ is reversible with respect to $\mu$ by
construction. Built upon a deterministic proposal of the backward trajectory
$\beta(\boldsymbol{z}_{0:L})\triangleq(\varphi(\boldsymbol{z}_{L}),\varphi(\boldsymbol{z}_{L-1}),\cdots,\varphi(\boldsymbol{z}_{0}))$,
we introduce a Metropolis update defined for any Borel set $A$ of
$\mathbb{R}^{2d(L+1)}$

$\mathbf{M}(\boldsymbol{z}_{0:L},A)\triangleq(1\wedge
e^{-\Delta(x_{0:L})})\delta_{\beta(\boldsymbol{z}_{0:L})}(A)+(1-1\wedge
e^{-\Delta(x_{0:L})})\delta_{\boldsymbol{z}_{0:L}}(A).$

The Metropolis kernel $\mathbf{M}$ is also reversible with respect to $\mu$.
For $\gamma=0$, this follows from (11), (12) and [2, Theorem 3]. For
$\gamma>0$ the distribution $\mu$ admits a density with respect to Lebesgue’s
measure. From (11) and [13, Eq 5.13] we obtain that the corresponding density
$\mu$ is such that

$-\Delta(\boldsymbol{x}_{0:L})=\sum_{i=1}^{L}\log\left(\frac{\Pi_{*}(\boldsymbol{z}_{i})q_{h,\gamma}(\varphi(\boldsymbol{z}_{i}),\varphi(\boldsymbol{z}_{i-1}))}{\Pi_{*}(\boldsymbol{z}_{i-1})q_{h,\gamma}(\boldsymbol{z}_{i-1},\boldsymbol{z}_{i})}\right)=\log\left(\frac{\mu(\beta(\boldsymbol{z}_{0:L}))}{\mu(\boldsymbol{z}_{0:L})}\right).$
(44)

Finally, we consider the cycle $\mathbf{K}=\mathbf{G}\mathbf{M}\mathbf{G}$
defined on Borel sets $A$ of $\mathbb{R}^{2d(L+1)}$

$\mathbf{K}(\boldsymbol{z}_{0:L},A)\triangleq\int\mathbf{G}(\boldsymbol{z}_{0:L},\mathrm{d}\boldsymbol{z}_{0:L}^{\prime})\mathbf{M}(\boldsymbol{z}_{0:L}^{\prime},\mathrm{d}\boldsymbol{z}_{0:L}^{\prime\prime})\mathbf{G}(\boldsymbol{z}_{0:L}^{\prime\prime},A).$

The palindromic structure of $\mathbf{K}=\mathbf{G}\mathbf{M}\mathbf{G}$
ensures reversibility with respect to $\mu$. Since the transition
$\mathbf{G}(\boldsymbol{z}_{0:L},.)$ only depends on the starting position
$\boldsymbol{x}_{0}\in\mathbb{R}^{d}$ and $\Pi$ is the marginal of $\mu$, we
obtain that
$\mathbf{P}(\boldsymbol{x}_{0},A)\triangleq\mathbf{K}(\boldsymbol{z}_{0:L},A\times\mathbb{R}^{d(2L+1)})$
defines marginally a Markov kernel on $\mathbb{R}^{d}$, reversible with
respect to $\Pi$. In particular, the distribution of
$(\boldsymbol{X}_{n})_{n\geq 0}$ in Algorithm 2 coincides with the
distribution of a Markov chain generated by $\mathbf{P}$.

###  B.3 Proof of Proposition 5

The claim of the proposition follows from establishing that $\mathbf{P}$ is
both $\Pi$-irreducible and aperiodic; see [66, Theorem 4]. We prove these two
results here. For any $\gamma>0$ and any position
$\boldsymbol{x}_{0}\in\mathbb{R}^{d}$, the kernel $\mathbf{G}$ in (43) defines
a density with respect to Lebesgue’s measure on $\mathbb{R}^{d(2L+1)}$

$g(\boldsymbol{x}_{0},(\boldsymbol{v}_{0},\boldsymbol{z}_{1:L}))\triangleq\psi(\boldsymbol{v}_{0})\prod_{i=1}^{L}q_{h,\gamma}(\boldsymbol{z}_{i-1},\boldsymbol{z}_{i}).$

By definition, this conditional density is positive everywhere, and
$\Delta(\boldsymbol{x}_{0:L})<\infty$ for any
$\boldsymbol{x}_{0:L}\in\mathbb{R}^{d(L+1)}$. Subsequently, for any Borel set
$A$ of $\mathbb{R}^{d}$ such that $\Pi(A)>0$, we have

$\mathbf{P}(\boldsymbol{x}_{0},A)\geq\int
g(\boldsymbol{x}_{0},(\boldsymbol{v}_{0},\boldsymbol{z}_{1:L}))(1\wedge
e^{-\Delta(x_{0:L})})\mathrm{d}(\boldsymbol{v}_{0},\boldsymbol{z}_{1:L})>0$
(45)

Since $\mathbf{P}(\boldsymbol{x}_{0},A)>0$ for any
$\boldsymbol{x}_{0}\in\mathbb{R}^{d}$ and any $A$ such that $\Pi(A)>0$, we
conclude that $\mathbf{P}$ is $\Pi$-irreducible. Suppose now that $\mathbf{P}$
is periodic. Then there exist two disjoint sets $A_{1},A_{2}$ with
$\Pi(A_{1})>0$ and $\Pi(A_{2})>0$ such that
$\mathbf{P}(\boldsymbol{x}_{0},A_{2})=1$ for any $\boldsymbol{x}_{0}\in
A_{1}$. However (45) implies that $\mathbf{P}(\boldsymbol{x}_{0},A_{1})>0$ for
any such $\boldsymbol{x}_{0}\in A_{1}$, yielding a contradiction. We conclude
that $\mathbf{P}$ is aperiodic.

##  Appendix C Proofs of optimal scaling

###  C.1 Proof of Theorem 2

Analogously to (10) and (11) denote for $x,y\in\mathbb{R}$

$\varepsilon_{h}(x,y)~{}=~{}\phi(y)-\phi(x)-\frac{y-x}{2}\left(\phi^{\prime}(y)+\phi^{\prime}(x\right))-\frac{h^{2}}{8}\left(\phi^{\prime}(y)^{2}-\phi^{\prime}(x)^{2}\right)$
(46)

and for $x_{0:L}\in\mathbb{R}^{L+1}$

$\Delta_{h}(x_{0:L})~{}=~{}\sum_{i=1}^{L}\varepsilon_{h}(x_{i-1},x_{i})\,.$
(47)

The product structure of the potential assumed in Assumption 3 enforces a
product structure also on the associated total energy differences in the
following sense

$\Delta(\boldsymbol{x}_{0:L})~{}=~{}\sum_{j=1}^{d}\Delta_{h,j}(\boldsymbol{x}_{0:L}(j))\,,$
(48)

where $\Delta_{h,j}(\boldsymbol{x}_{0:L}(j))$ are IID copies of a random
variable $\Delta_{h}(\boldsymbol{x}_{0:L})$ given by a single component of the
Langevin trajectory $x_{0:L}$, with respect to the one-dimensional potential
$\phi$ (initiated in stationarity).

To prove that total energy difference random variables
$\Delta(\boldsymbol{x}_{0:L})$ satisfy a form of CLT, we use an extension of
the framework introduced in [73, Section 3]. More precisely, Theorem 2 is as a
direct application of [73, Theorem 8]. To invoke it, we need to understand the
asymptotic behavior of a single component of the Langevin trajectory as the
dimension increases. Specifically, we need to verify the following two
conditions:

######  Proposition 9.

Let $T>0$, $\gamma\geq 0$ and $\ell>0$ and take a sequence of time-steps $h\to
0$ and $L=\lfloor T/h\rfloor$. Assume the one-dimensional potential $\phi$
satisfies Assumption 3 and let $x_{0}\sim\Pi$. Then

  1. (i)

$\frac{1}{h^{4}}\mathbb{E}\left[\Delta_{h}^{2}(x_{0:L})\right]\quad\xrightarrow{h\to
0}\quad\Sigma~{}=~{}\mathbb{E}\left[\left(\int_{0}^{T}S(X_{t},V_{t})\mathrm{d}t\right)^{2}\right]$

  2. (ii)

$\frac{1}{h^{4}}\mathbb{E}\left[\Delta_{h}^{2}(x_{0:L})\mathds{1}_{\Delta_{h}>h}\right]\quad\xrightarrow{h\to
0}\quad 0\,.$

A proof is given is Section C.2.

Some care is required to justify that the results of [73, Section 3] do
actually carry over from the classical Metropolis-Hastings setting to MALT.
The results of [73, Section 3] concern objects called _log Metropolis-Hastings
random variables_ that are related to the Kullback-Leibler divergence between
the forward $\Pi(\mathrm{d}x)q(x,\mathrm{d}y)$ and reverse
$\Pi(\mathrm{d}y)q(y,\mathrm{d}x)$ transition kernels. The total energy
difference random variables defined here are a generalisation, they can also
be related to the Kullback-Leibler divergence between the forward distribution
$\mu(z_{0:L})$ and the skew-backward distribution $\mu(\beta(z_{0:L}))$.

The proofs of [73, Section 3] depend entirely on the symmetry property of log-
MH random variables featured in the following proposition. They do not rely on
the exact definition of log-MH random variables, only on their symmetric
property. Therefore, the only thing we need to do in order to extend the
results to the trajectory setting is to verify this symmetry is still
satisfied, the proofs from that point on remain literally the same.

######  Proposition 10.

Let $\Delta(\boldsymbol{x}_{0:L})$ be a total energy difference function of
MALT and abbreviate
$\boldsymbol{x}_{L:0}=(\boldsymbol{x}_{L},\dots,\boldsymbol{x}_{0})$. Then

  1. (i)

$\Delta(\boldsymbol{x}_{0:L})~{}=~{}-\Delta(\boldsymbol{x}_{L:0})$.

  2. (ii)

Let $f:\mathbb{R}\to\mathbb{R}$ be a measurable functions such that
$f\circ\Delta$ is integrable with respect to $\mu$. Then
$\mathbb{E}_{\mu}[f(-\Delta)e^{-\Delta}]=\mathbb{E}_{\mu}[f(\Delta)]$.

###### Proof.

Part (i) of the proposition is immediate from (10) and (11). To establish Part
(ii) we split the cases $\gamma>0$ and $\gamma=0$. In the positive friction
case, by (44) the measure $\mu$ has a positive density with respect to
Lebesgue measure on $\mathbb{R}^{2d(L+1)}$ such that

$\mu(\boldsymbol{z}_{0:L})e^{-\Delta(\boldsymbol{x}_{0:L})}~{}=~{}\mu(\beta(\boldsymbol{z}_{0:L}))\,.$

Combining this equality to the change of variable
$\boldsymbol{z}^{\prime}_{0:L}=\beta(\boldsymbol{z}_{0:L})$ with unit
Jacobian, we obtain

$\displaystyle\mathbb{E}_{\mu}[|f(-\Delta(\boldsymbol{x}_{0:L}))|e^{-\Delta(\boldsymbol{x}_{0:L})}]~{}$
$\displaystyle=~{}{\int}_{\mathbb{R}^{2d(L+1)}}|f(-\Delta(\boldsymbol{x}_{0:L}))|e^{-\Delta(\boldsymbol{x}_{0:L})}\mu(\boldsymbol{z}_{0:L})\mathrm{d}\boldsymbol{z}_{0:L}$
$\displaystyle=~{}{\int}_{\mathbb{R}^{2d(L+1)}}|f(\Delta(\boldsymbol{x}_{L:0}))|\mu(\beta(\boldsymbol{z}_{0:L}))\mathrm{d}\boldsymbol{z}_{0:L}$
$\displaystyle=~{}{\int}_{\mathbb{R}^{2d(L+1)}}|f(\Delta(\boldsymbol{x}^{\prime}_{0:L}))|\mu(\boldsymbol{z}^{\prime}_{0:L})\mathrm{d}\boldsymbol{z}^{\prime}_{0:L}~{}<~{}\infty\,.$

This ensures that the quantity
$\mathbb{E}[f(-\Delta(\boldsymbol{x}_{0:L}))e^{-\Delta(\boldsymbol{x}_{0:L})}]$
is integrable. Repeating the calculation without the absolute values yields
the desired identity. Now, when $\gamma=0$ the measure $\mu$ is characterized
by the density $\Pi_{*}$ on $\mathbb{R}^{2d}$. Each trajectory is a
deterministic function of its starting point
$\boldsymbol{z}_{0:L}=(\boldsymbol{z}_{0},\cdots,\boldsymbol{\theta}_{h}^{L}(\boldsymbol{z}_{0}))$.
We apply the change of variable
$\boldsymbol{z}^{\prime}_{0}=\varphi\circ\boldsymbol{\theta}^{L}_{h}(\boldsymbol{z}_{0})$
with unit Jacobian. Denote
$\boldsymbol{z}^{\prime}_{0:L}=(\boldsymbol{z}^{\prime}_{0},\cdots,\boldsymbol{\theta}_{h}^{L}(\boldsymbol{z}^{\prime}_{0}))$
and remark that time-reversibility of the Leapfrog update yields
$\boldsymbol{z}^{\prime}_{0:L}=(\varphi\circ\boldsymbol{\theta}_{h}^{L}(\boldsymbol{z}_{0}),\cdots,\varphi(\boldsymbol{z}_{0}))$.
Since $\varphi$ leaves $\Pi_{*}$ invariant, we obtain

$\displaystyle\mathbb{E}_{\mu}[|f(-\Delta(\boldsymbol{x}_{0:L}))|e^{-\Delta(\boldsymbol{x}_{0:L})}]$
$\displaystyle\quad=~{}{\int}_{\mathbb{R}^{2d}}\left|f\left(-\Delta\left(\boldsymbol{x}_{0:L}\right)\right)\right|e^{-\Delta\left(\boldsymbol{x}_{0:L}\right)}\Pi_{*}(\boldsymbol{z}_{0})\mathrm{d}\boldsymbol{z}_{0}$
$\displaystyle\quad=~{}{\int}_{\mathbb{R}^{2d}}\left|f\left(\Delta\left(\boldsymbol{x}_{L:0}\right)\right)\right|\Pi_{*}(\boldsymbol{\theta}^{L}_{h}(z_{0}))\mathrm{d}\boldsymbol{z}_{0}$
$\displaystyle\quad=~{}{\int}_{\mathbb{R}^{2d}}\left|f\left(\Delta\left(\boldsymbol{x}^{\prime}_{0:L}\right)\right)\right|\Pi_{*}(\boldsymbol{z}^{\prime}_{0})\mathrm{d}\boldsymbol{z}^{\prime}_{0}~{}<~{}\infty\,.$

Again, the same calculation without absolute values establishes the desired
identity. ∎

###  C.2 Proof of Proposition 9

The center of the proof for both (i) and (ii) is showing that

$\frac{1}{h^{2}}\Delta_{h}(x_{0:L})~{}=~{}\int_{0}^{T}S(X_{t},V_{t})\mathrm{d}t~{}+~{}R_{h}\,,$
(49)

where the reminder $R_{h}$ satisfies $\mathbb{E}[R_{h}^{2}]\to 0$ as
$h\to\infty$. This is done by approximating the leading terms of a Taylor’s
expansion of $\Delta_{h}(x_{0:L})$ over the entire Langevin trajectory and
controlling the square expectations of the error terms. For this task we need
to extend strong accuracy results (see Proposition 3) beyond Lipschitz
functions. This is enabled by the following lemma proven in Section C.3.

######  Lemma 4.

Let $T>0$, $\gamma\geq 0$ and $h\in(0,1]$ and take $L$ so that $hL\leq T$.
Assume the one-dimensional potential $\phi$ satisfies Assumption 3 and let
$x_{0}\sim\Pi$. There exists constants $A_{1},A_{2},A_{3},A_{4}>0$ independent
of the choices of $h$ and $L$ such that following statements hold:

  1. (i)

$\max_{0\leq i\leq
L}\max(\mathbb{E}[x_{i}^{8}],\mathbb{E}[v_{i}^{8}])~{}\leq~{}A_{1}\,,$

  2. (ii)

$\mathbb{E}\left[\left(\phi^{\prime}(x_{L})^{2}-\phi^{\prime}(X_{hL})^{2}\right)^{2}\right]~{}\leq~{}A_{2}h\,,$

  3. (iii)

$\sup_{1\leq i\leq
L}\mathbb{E}\left[\left(\left(\frac{x_{i}-x_{i-1}}{h}\right)^{3}-\hat{V}_{i-1}^{3}\right)^{2}\right]~{}\leq~{}A_{3}h\,,$

  4. (iv)

$\sup_{1\leq i\leq
L}\mathbb{E}\left[\left(V_{(i-1)h}^{3}-v_{i-1}^{3}\right)^{2}\right]~{}\leq~{}A_{4}h^{2/3}\,.$

To show (49) we transform the expression in (47) into the desired form. The
sum of third terms of (46) terms appearing in (47) can be dealt with using
Itô’s lemma

$\int_{0}^{T}V_{t}\phi^{\prime\prime}(X_{t})\phi^{\prime}(X_{t})\mathrm{d}t~{}=~{}\int_{0}^{T}\frac{\partial}{\partial
x}\left(\frac{1}{2}\phi^{\prime}(X_{t})^{2}\right)\mathrm{d}X_{t}~{}=~{}\frac{1}{2}\left(\phi^{\prime}(X_{T})^{2}-\phi^{\prime}(X_{0})^{2}\right).$

Since $x_{0}=X_{0}$, this implies that

$\displaystyle\begin{split}\mathbb{E}\bigg{[}\Big{(}&\sum_{i=1}^{L}\big{(}\phi^{\prime}(x_{i})^{2}-\phi^{\prime}(x_{i-1})^{2}\big{)}~{}-~{}2\int_{0}^{T}V_{t}\phi^{\prime\prime}(X_{t})\phi^{\prime}(X_{t})\mathrm{d}t\Big{)}^{2}\bigg{]}\\\
&=~{}\mathbb{E}\left[\left(\phi^{\prime}(x_{L})^{2}-\phi^{\prime}(x_{0})^{2}-\phi^{\prime}(X_{T})^{2}+\phi^{\prime}(X_{0})^{2}\right)^{2}\right]\\\
&\leq~{}2\mathbb{E}\left[\left(\phi^{\prime}(x_{L})^{2}-\phi^{\prime}(X_{hL})^{2}\right)^{2}\right]~{}+~{}2\mathbb{E}\left[\left(\phi^{\prime}(X_{T})^{2}-\phi^{\prime}(X_{hL})^{2}\right)^{2}\right]\quad\xrightarrow{h\to
0}\quad 0\,.\end{split}$ (50)

The first term vanishes by Lemma 4 (ii). The second vanishes because $hL\to
T$, the function $\phi^{\prime}$ is $M$-Lipschitz (grows at most linearly) and
$X_{t}$ is a stationary process with continuous paths and finite fourth
moments.

To control the rest of (47) we take advantage of [73, Lemma A.2.] and the
fundamental theorem of calculus. For each $1\leq i\leq L$ we rewrite using
$\int_{0}^{1}(1-2u)\mathrm{d}u=0$

$\displaystyle\phi(x_{i})-\phi(x_{i-1})$
$\displaystyle-\frac{1}{2}(x_{i}-x_{i-1})\left(\phi^{\prime}(x_{i})+\phi^{\prime}(x_{i-1}\right)$
$\displaystyle=~{}\frac{1}{2}(x_{i}-x_{i-1})^{2}\int_{0}^{1}(1-2u)\phi^{\prime\prime}(x_{i-1}+u(x_{i}-x_{i-1}))\mathrm{d}u$
$\displaystyle=~{}\frac{1}{2}(x_{i}-x_{i-1})^{3}\int_{0}^{1}(1-2u)u\int_{0}^{1}\phi^{(3)}(x_{i-1}+su(x_{i}-x_{i-1}))\mathrm{d}s\mathrm{d}u\,.$

Noting that $\int_{0}^{1}(1-2u)u\mathrm{d}u=-\tfrac{1}{6}$, we can relate this
to a term expressed with Langevin dynamics:

$\displaystyle\begin{split}\phi(x_{i})-\phi(x_{i-1})&-\frac{1}{2}(x_{i}-x_{i-1})\left(\phi^{\prime}(x_{i})+\phi^{\prime}(x_{i-1}\right)+\frac{1}{12}h^{3}V_{(i-1)h}^{3}\phi^{(3)}(X_{(i-1)h})\\\
&=~{}R_{i,1}~{}+~{}R_{i,2}~{}+~{}R_{i,3}\,,\end{split}$ (51)

where

$\displaystyle R_{i,1}$
$\displaystyle=\frac{1}{2}\left((x_{i}-x_{i-1})^{3}-h^{3}v_{i-1}^{3}\right)\int_{0}^{1}(1-2u)u\int_{0}^{1}\phi^{(3)}(x_{i-1}+su(x_{i}-x_{i-1}))\mathrm{d}s\mathrm{d}u$
$\displaystyle R_{i,2}$
$\displaystyle=\frac{1}{2}h^{3}\left(v_{i-1}^{3}-V_{(i-1)h}^{3}\right)\int_{0}^{1}(1-2u)u\int_{0}^{1}\phi^{(3)}(x_{i-1}+su(x_{i}-x_{i-1}))\mathrm{d}s\mathrm{d}u$
$\displaystyle R_{i,3}$
$\displaystyle=\frac{1}{2}h^{3}V_{(i-1)h}^{3}\int_{0}^{1}(1-2u)u\int_{0}^{1}\left(\phi^{(3)}(x_{i-1}+su(x_{i}-x_{i-1}))-\phi^{(3)}(X_{(i-1)h})\right)\mathrm{d}s\mathrm{d}u\,.$

We will show that of the $k=1,2,3$ we have

$\frac{1}{h^{6}}\sup_{1\leq i\leq
L}\mathbb{E}[R_{i,k}^{2}]\quad\xrightarrow{h\to 0}\quad 0\,.$ (52)

For the first two term $R_{i,1}$ and $R_{i,2}$ we use boundedness of
$\phi^{(3)}$ together with respectively points (iii) and (iv) of Lemma 4. For
the last term $R_{i,3}$ we first use Cauchy’s inequality (recall
$\mathbb{E}[V_{(i-1)h}^{12}]$ is finite since $V_{t}$ is a stationary standard
Gaussian process) together with Jensen’s inequality, Fubini’s theorem and the
fact that both $\phi^{(3)}$ and $\phi^{(4)}$ are bounded.

$\mathbb{E}\left[\left(\phi^{(3)}(x_{i-1}+su(x_{i}-x_{i-1}))-\phi^{(3)}(X_{(i-1)h})\right)^{4}\right]\\\
\leq~{}8\|\phi^{(4)}\|_{\infty}\|\phi^{(3)}\|_{\infty}^{3}\left(\mathbb{E}\left[\left|x_{i-1}-X_{(i-1)h}\right|\right]+\mathbb{E}\left[\left|x_{i}-x_{i-1}\right|\right]\right)\,.$

Then use Proposition 3 and Lemma 4 (iii) and (iv) to control this bound as
follows:

$\displaystyle\mathbb{E}[|x_{i}-$ $\displaystyle
x_{i-1}|]\leq\mathbb{E}\left[\left(x_{i}-x_{i-1}\right)^{6}\right]^{1/6}$
$\displaystyle=\mathbb{E}\left[\left(\left(x_{i}-x_{i-1}\right)^{3}-h^{3}v_{i-1}^{3}+h^{3}\left(v_{i-1}^{3}-V_{(i-1)h}^{3}\right)+h^{3}V_{(i-1)h}^{3}\right)^{2}\right]^{1/6}$
$\displaystyle\leq
3h\left(\mathbb{E}\Big{[}\big{(}\left((x_{i}-x_{i-1})/h\right)^{3}-v_{i-1}^{3}\big{)}^{2}\Big{]}+\mathbb{E}\Big{[}\big{(}v_{i-1}^{3}-V_{(i-1)h}^{3}\big{)}^{2}\Big{]}+\mathbb{E}\left[V_{(i-1)h}^{6}\right]\right)^{1/6}$

Finally, paths of $(X_{t},V_{t})$ are almost surely continuous and therefore
Riemann integrable. Hence almost surely

$h\sum_{i=1}^{L}V_{(i-1)h}^{3}\phi^{(3)}(X_{(i-1)h})\quad\xrightarrow{h\to
0}\quad\int_{0}^{T}V^{3}_{t}\phi^{(3)}(X_{t})\mathrm{d}t\,.$

Since $hL\to T$, $\phi^{(3)}$ is bounded and $V_{t}$ is standard Gaussian for
all $t$, the convergence also holds in $\mathbb{L}_{2}$-sense

$\mathbb{E}\left[\left(h\sum_{i=1}^{L}V_{(i-1)h}^{3}\phi^{(3)}(X_{(i-1)h})-\int_{0}^{T}V^{3}_{t}\phi^{(3)}(X_{t})\mathrm{d}t\right)^{2}\right]\quad\xrightarrow{h\to
0}\quad 0\,.$

Combining this convergence with (51) and (52), we obtain as $h\to 0$

$\frac{1}{h^{4}}\mathbb{E}\bigg{[}\Big{(}\sum_{i=1}^{L}\big{(}\phi(x_{i})-\phi(x_{i-1})-\frac{1}{2}(x_{i}-x_{i-1})(\phi^{\prime}(x_{i})+\phi^{\prime}(x_{i-1}))\big{)}-\int_{0}^{T}V^{3}_{t}\phi^{(3)}(X_{t})\mathrm{d}t\Big{)}^{2}\bigg{]}\to
0$

which, combined with (50), yields (49). Point (i) follows by integrating

$\frac{1}{h^{4}}\mathbb{E}\left[\Delta^{2}_{h}(x_{0:L})\right]~{}=~{}\mathbb{E}\left[\left(\int_{0}^{T}S(X_{t},V_{t})\mathrm{d}t\right)^{2}\right]-2\mathbb{E}\left[R_{h}\int_{0}^{T}S(X_{t},V_{t})\mathrm{d}t\right]+\mathbb{E}\left[R_{h}^{2}\right]$

where the last two terms vanish as $h\to 0$ (use Cauchy’s inequality for the
middle term).

Proof of point Proposition 9(ii) uses similar arguments. Recall that $V_{t}$
is a stationary Gaussian and that by Assumption 3 there exists constants $M,D$
such that $\|\phi^{\prime\prime}\|_{\infty}\leq M$ and $|\phi^{\prime}(x)|\leq
M(|x|+D)$. Note that by Assumption 3, Jensen’s inequality and Fubini’s theorem

$\displaystyle\mathbb{E}\left[\Big{(}\int_{0}^{T}S(X_{t},V_{t})\mathrm{d}t\Big{)}^{4}\right]$
$\displaystyle\leq\int_{0}^{T}\mathbb{E}\left[S(X_{t},V_{t})^{4}\right]\mathrm{d}t$
$\displaystyle=4\int_{0}^{T}\mathbb{E}\left[\frac{1}{12^{4}}V_{t}^{12}\phi^{(3)}(X_{t})^{4}+\frac{1}{4^{4}}V_{t}^{4}\phi^{\prime\prime}(X_{t})^{4}\phi^{\prime}(X_{t})^{4}\right]\mathrm{d}t$
$\displaystyle\leq\frac{T}{3^{4}4^{3}}\|\phi^{(3)}\|_{\infty}^{4}\mathbb{E}[V_{t}^{12}]+\frac{T}{4^{3}}M^{4}\mathbb{E}[V_{t}^{8}]^{1/2}\mathbb{E}\left[(M(|X_{t}|+D))^{8}\right]^{1/2}$
$\displaystyle<\infty.$

This and (49) show

$\displaystyle\frac{1}{h^{4}}\mathbb{E}\left[\Delta^{2}_{h}(x_{0:L})\mathds{1}_{\Delta_{h}(x_{0:L})>h}\right]$
$\displaystyle=\mathbb{E}\left[\Big{(}\int_{0}^{T}S(X_{t},V_{t})\mathrm{d}t\Big{)}^{2}\mathds{1}_{\Delta_{h}(x_{0:L})>h}\right]+\mathbb{E}\left[R_{h}^{2}\mathds{1}_{\Delta_{h}(x_{0:L})>h}\right]$
$\displaystyle+2\mathbb{E}\left[R_{h}\left(\int_{0}^{T}S(X_{t},V_{t})\mathrm{d}t\right)\mathds{1}_{\Delta_{h}(x_{0:L})>h}\right].$

The last two terms vanish since $\mathbb{E}[R_{h}^{2}]\to 0$. The first term
vanishes by a combination of Cauchy and Markov inequalities together with part
(i)

$\displaystyle\mathbb{E}\bigg{[}\Big{(}\int_{0}^{T}S(X_{t},V_{t})\mathrm{d}t\Big{)}^{2}$
$\displaystyle\mathds{1}_{\Delta_{h}(x_{0:L})>h}\bigg{]}\leq\bigg{(}\mathbb{E}\bigg{[}\Big{(}\int_{0}^{T}S(X_{t},V_{t})\mathrm{d}t\Big{)}^{4}\bigg{]}\mathbb{P}\left(\Delta^{2}_{h}(x_{0:L})>h^{2}\right)\bigg{)}^{1/2}$
$\displaystyle~{}\leq
h\,\mathbb{E}\left[\Big{(}\int_{0}^{T}S(X_{t},V_{t})\mathrm{d}t\Big{)}^{4}\right]^{1/2}\left(\frac{1}{h^{4}}\mathbb{E}\left[\Delta_{h}^{2}(x_{0:L})\right]\right)^{1/2}\xrightarrow{h\to
0}~{}0\,.$

###  C.3 Proof of Lemma 4

A uniform control of the distance between the numerical Langevin trajectory
and the Langevin diffusion is given by Proposition 3, which extends trivially
to functional values along the trajectories for all Lipschitz functions. Point
(i) allows us to extend these findings for functions that are not globally
Lipschitz but do not grow to rapidly. This is then done for specific functions
required for the analysis of MALT in points (ii), (iii) and (iv).

The proof of (i) we use a ”discrete Grönwall’s inequality”. We will show that
for $0\leq i\leq L-1$ and some constants $A,B>0$

$\max\left(\mathbb{E}[x_{i+1}^{8}],\mathbb{E}[v_{i+1}^{8}],B\right)~{}\leq~{}(1+Ah)\times\max\left(\mathbb{E}[x_{i}^{8}],\mathbb{E}[v_{i}^{8}],B\right)\,.$
(53)

Iterating this inequality and using $h\leq T/L$ and the inequality $(1+t)\leq
e^{t}$ gives us the following uniform bound, valid for all $0\leq i\leq L$:

$\displaystyle\max\left(\mathbb{E}\left[x_{i}^{8}\right],\mathbb{E}\left[v_{i}^{8}\right],B\right)~{}$
$\displaystyle\leq~{}(1+Ah)^{i}\times\max\left(\mathbb{E}\left[x_{0}^{8}\right],\mathbb{E}\left[v_{0}^{8}\right],B\right)$
$\displaystyle\leq~{}(1+AT/L)^{L}\times\max\left(\mathbb{E}\left[x_{0}^{8}\right],\mathbb{E}\left[v_{0}^{8}\right],B\right)$
$\displaystyle\leq~{}e^{AT}\times\max\left(\mathbb{E}\left[x_{0}^{8}\right],\mathbb{E}\left[v_{0}^{8}\right],B\right)\,.$

This establishes (i) with
$A_{1}=e^{AT}\max\left(\mathbb{E}\left[x_{0}^{8}\right],\mathbb{E}\left[v_{0}^{8}\right],B\right)$,
which is finite by Assumption 3, as the Langevin trajectory is initiated in
stationarity.

Assumption 3 implies that $\phi^{\prime}$ is $M$-Lipschitz and grows at most
linearly: $|\phi^{\prime}(x)|\leq M(|x|+D)$ for appropriate constants $M,D$,
so that $(\phi^{\prime}(x))^{8}\leq 8M^{8}(|x|^{8}+D^{8})$.

First, we deal with the bound on the eighth moment of the position. We use the
Multinomial theorem to expand the expression of $x_{i+1}^{8}$ in the OBABO
update (B.1). We then separate the zero-order term with respect to $h$ from
the rest and use the Hölder’s inequality to bound them. Note that the number
of terms is fixed and finite and recall that we are assuming $h<1$. Using the
fact that $\phi^{\prime}$ grows at most linearly, we get for appropriate
constants $A^{\prime},B^{\prime}>0$

$\displaystyle\begin{split}\mathbb{E}[x_{i+1}^{8}]&=\mathbb{E}\left[\left(x_{i}+he^{-\gamma
h/2}v_{i}+h\sqrt{1-e^{-\gamma
h}}\xi_{i}-\frac{h^{2}}{2}\phi^{\prime}(x_{i})\right)^{8}\right]\\\
&\leq\mathbb{E}[x_{i}^{8}]+\sum_{\begin{subarray}{c}j_{1}+j_{2}+j_{3}+j_{4}=8\\\
j_{1}<8\end{subarray}}h^{j_{2}+j_{3}+2j_{4}}\binom{8}{j_{1},j_{2},j_{3},j_{4}}\mathbb{E}\left[|x_{i}|^{j_{1}}|v_{i}|^{j_{2}}|\xi_{i}|^{j_{3}}|\phi^{\prime}(x_{i})|^{j_{4}}\right]\\\
&\leq\mathbb{E}[x_{i}^{8}]+(8!)h\sum_{\begin{subarray}{c}j_{1}+j_{2}+j_{3}+j_{4}=8\\\
j_{1}<8\end{subarray}}\left(\mathbb{E}\left[x_{i}^{8}\right]^{j_{1}}\mathbb{E}\left[v_{i}^{8}\right]^{j_{2}}\mathbb{E}\left[\xi_{i}^{8}\right]^{j_{3}}\mathbb{E}\left[\phi^{\prime}(x_{i})^{8}\right]^{j_{4}}\right)^{1/8}\\\
&\leq\mathbb{E}[x_{i}^{8}]+A^{\prime}h\max\left(\mathbb{E}[x_{i}^{8}],\mathbb{E}[v_{i}^{8}],B^{\prime}\right)\,.\end{split}$
(54)

We derive a similar identity for the update of the velocities where we use
Multinomial theorem to expand the expression of $v_{i+1}^{8}$ in the OBABO
update (B.1). Note that since $\xi_{i},\xi_{i+1}$ are standard Gaussian and
independent of each other and of $v_{i}$, there exists a standard Gaussian
$\xi_{i}^{\prime}$ that is independent of $v_{i}$ and satisfies $e^{-\gamma
h/2}\sqrt{1-e^{-\gamma h}}\xi_{i}+\sqrt{1-e^{-\gamma
h}}\xi_{i+1}=\sqrt{1-e^{-2\gamma h}}\xi^{\prime}_{i}$. Also note that the
inequality $1-e^{-t}\leq t$ implies that $\sqrt{1-e^{-2\gamma
h}}<\sqrt{2\gamma h}$. We separate the zero-order and half-order terms with
respect to $h$ and bound the remaining terms using the Lipschitz property of
$\phi^{\prime}$ and (54). The half-order term vanishes since
$\mathbb{E}[\xi^{\prime}_{i}]=0$. For appropriate constants $A,A^{\prime},B>0$
we get

$\displaystyle\begin{split}\mathbb{E}[v_{i+1}^{8}]~{}&=~{}\mathbb{E}\left[\left(e^{-\gamma
h}v_{i}-h\phi^{\prime}(x_{i})-\frac{h}{2}\left(\phi^{\prime}(x_{i+1})-\phi^{\prime}(x_{i})\right)+\sqrt{1-e^{-2\gamma
h}}\xi^{\prime}_{i}\right)^{8}\right]\\\
&\leq~{}\mathbb{E}[v_{i}^{8}]~{}+~{}8e^{-7\gamma h}\sqrt{1-e^{-2\gamma
h}}\mathbb{E}\left[v_{i}^{7}\right]\mathbb{E}\left[\xi^{\prime}_{i}\right]\\\
&+~{}A^{\prime}\sum_{\begin{subarray}{c}j_{1}+j_{2}+j_{3}+j_{4}=8\\\
j_{2}+j_{3}+\tfrac{1}{2}j_{4}\geq
1\end{subarray}}h^{j_{2}+j_{3}+\tfrac{1}{2}j_{4}}\mathbb{E}\left[|v_{i}|^{j_{1}}|\phi^{\prime}(x_{i})|^{j_{2}}\left|\phi^{\prime}(x_{i+1})-\phi^{\prime}(x_{i})\right|^{j_{3}}|\xi^{\prime}_{i}|^{j_{4}}\right]\\\
&\leq~{}\mathbb{E}[v_{i}^{8}]~{}+~{}Ah\max\left(\mathbb{E}\left[x_{i}^{8}\right],\mathbb{E}\left[v_{i}^{8}\right],B\right)\,.\end{split}$
(55)

Together, (54) and (55) imply (53) and establish (i).

We obtain (ii) from combining: (i), the Lipschitz property of $\phi^{\prime}$
and Cauchy’s inequality

$\displaystyle\mathbb{E}\left[\left(\phi^{\prime}(x_{L})^{2}-\phi^{\prime}(X_{hL})^{2}\right)^{2}\right]~{}$
$\displaystyle=~{}\mathbb{E}\left[\left(\phi^{\prime}(x_{L})-\phi^{\prime}(X_{hL})\right)^{2}\left(\phi^{\prime}(x_{L})+\phi^{\prime}(X_{hL})\right)^{2}\right]$
$\displaystyle\leq~{}\mathbb{E}\left[\left|\phi^{\prime}(x_{L})-\phi^{\prime}(X_{hL})\right|\left(|\phi^{\prime}(x_{L})|+|\phi^{\prime}(X_{hL})|\right)^{3}\right]$
$\displaystyle\leq~{}M^{4}\times\mathbb{E}\left[\left|x_{L}-X_{hL}\right|\left(|x_{L}|+|X_{hL}|+2D\right)^{3}\right]$
$\displaystyle\leq~{}M^{4}\times\mathbb{E}\left[\left|x_{L}-X_{hL}\right|^{2}\right]^{\frac{1}{2}}\times\mathbb{E}\left[\left(|x_{L}|+|X_{hL}|+2D\right)^{6}\right]^{\frac{1}{2}}$
$\displaystyle\leq~{}A_{2}h\,.$

The final bound follows for an appropriate constant by Proposition 3 and (i).

To show (iii) consider the following bound for $a,b\in\mathbb{R}$

$\displaystyle(a^{3}-b^{3})^{2}~{}$
$\displaystyle=~{}(a-b)^{2}(a^{2}+ab+b^{2})^{2}$
$\displaystyle=~{}(a-b)^{2}((a-b)^{2}+3(a-b)b+3b^{2})^{2}$
$\displaystyle\leq~{}3(a-b)^{6}+27(a-b)^{4}b^{2}+27(a-b)^{2}b^{4}\,.$

Fix $1\leq i\leq L$. Insert $a=h^{-1}(x_{i}-x_{i-1})$ and $b=v_{i-1}$ in the
above bound together with Hölder’s inequality gives us

$\mathbb{E}\left[\left(h^{-3}(x_{i}-x_{i-1})^{3}-v_{i-1}^{3}\right)^{2}\right]\,\leq\,27\,\mathbb{E}\left[\left(h^{-1}(x_{i}-x_{i-1})-v_{i-1}\right)^{6}\right]^{1/3}\mathbb{E}\left[v_{i-1}^{6}\right]^{2/3}\\\
+~{}27\,\mathbb{E}\left[\left(h^{-1}(x_{i}-x_{i-1})-v_{i-1}\right)^{6}\right]^{2/3}\mathbb{E}\left[v_{i-1}^{6}\right]^{1/3}~{}+~{}3\,\mathbb{E}\left[\left(h^{-1}(x_{i}-x_{i-1})-v_{i-1}\right)^{6}\right]\,.$

Hence, it is sufficient to show that for an appropriate constant
$A_{3}^{\prime}>0$

$\mathbb{E}\left[\left(h^{-1}(x_{i}-x_{i-1})-v_{i-1}\right)^{6}\right]~{}\leq~{}A_{3}^{\prime}h^{3}\,.$

Using the linear growth of $\phi^{\prime}$ and the inequality $1-e^{-t}\leq t$

$\displaystyle\mathbb{E}\left[\left(h^{-1}(x_{i}-x_{i-1})-v_{i-1}\right)^{6}\right]~{}$
$\displaystyle=~{}\mathbb{E}\Big{[}\Big{(}(e^{-\gamma
h/2}-1)v_{i-1}+\sqrt{1-e^{-\gamma
h}}\xi_{0}-\frac{h}{2}\phi^{\prime}(x_{i-1})\Big{)}^{6}\Big{]}$
$\displaystyle\leq~{}6(\gamma
h/2)^{6}\mathbb{E}\left[v_{i-1}^{6}\right]~{}+~{}6(\gamma
h)^{3}\mathbb{E}\left[\xi_{0}^{6}\right]+\frac{3h^{6}}{32}\left[\phi^{\prime}(x_{i-1})^{6}\right]$
$\displaystyle\leq~{}A^{\prime}_{3}h^{3}\,.$

To prove (iv) consider the following bound for $a,b\in\mathbb{R}$

$\displaystyle(a^{3}-b^{3})^{2}~{}$
$\displaystyle=~{}(a-b)^{2}(a^{2}+ab+b^{2})^{2}$
$\displaystyle\leq~{}(a-b)^{2}(|a|+|b|)^{4}~{}\leq~{}(a-b)^{2/3}(|a|+|b|)^{16/3}\,.$

Fix $1\leq i\leq L$. Inserting $a=V_{(i-1)h}$ and $b=v_{i-1}$ in the above
bound together with Hölder’s inequality gives us

$\mathbb{E}\left[\left(V_{(i-1)h}^{3}-v_{i-1}^{3}\right)^{2}\right]\leq\mathbb{E}\left[\left(V_{(i-1)h}-v_{i-1}\right)^{2}\right]^{1/3}\mathbb{E}\left[\left(|V_{(i-1)h}|+|v_{i-1}|\right)^{8}\right]^{2/3}\,.$

Proposition 3 implies that
$\mathbb{E}\left[\left(V_{(i-1)h}-v_{i-1}\right)^{2}\right]\leq
A^{\prime}_{4}h^{2}$ for an appropriate constant. Point (i) together with the
stationarity of process $V_{t}$ ensure the right factor is bounded.

###  C.4 Proofs of miscellaneous results

###### Proof of Proposition 6.

(i) This is a direct consequence of weak convergence in Theorem 2 applied to
the Lipschitz function $t\to 1\wedge e^{t}$. We use Proposition 2.4 in [64] to
evaluate $\mathbb{E}[1\wedge e^{Y}]$ for a normal random variable $Y$ on
$\mathbb{R}$.

(ii) By independence of the coordinates and the product structure (48)

$\displaystyle\mathbb{E}\Big{[}\big{(}f(\boldsymbol{X}^{n+1}(1))-$
$\displaystyle
f(\boldsymbol{X}^{n}(1))\big{)}^{2}\Big{]}=~{}\mathbb{E}\left[\left(f(\boldsymbol{x}_{L}(1))-f(\boldsymbol{x}_{0}(1))\right)^{2}1\wedge
e^{-\sum_{j=1}^{d}\Delta_{h,j}}\right]$
$\displaystyle=~{}\mathbb{E}\left[\left(f(\boldsymbol{x}_{L}(1))-f(\boldsymbol{x}_{0}(1))\right)^{2}\left(1\wedge
e^{-\sum_{j=1}^{d}\Delta_{h,j}}-1\wedge
e^{-\sum_{j=2}^{d}\Delta_{h,j}}\right)\right]$
$\displaystyle+~{}\mathbb{E}\left[\left(f(\boldsymbol{x}_{L}(1))-f(\boldsymbol{x}_{0}(1))\right)^{2}\right]\mathbb{E}\left[1\wedge
e^{-\sum_{j=2}^{d}\Delta_{h,j}}\right]$

We first show that the second term vanishes with increasing $d$. Using
Cauchy’s inequality and the fact that $t\mapsto 1\wedge e^{t}$ is
$1$-Lipschitz we argue

$\displaystyle\mathbb{E}\Big{[}\big{(}f(\boldsymbol{x}_{L}(1))-$
$\displaystyle f(\boldsymbol{x}_{0}(1))\big{)}^{2}\Big{(}1\wedge
e^{-\sum_{j=1}^{d}\Delta_{h,j}}-1\wedge
e^{-\sum_{j=2}^{d}\Delta_{h,j}}\Big{)}\Big{]}$
$\displaystyle\leq~{}\mathbb{E}\left[\left(f(\boldsymbol{x}_{L}(1))-f(\boldsymbol{x}_{0}(1))\right)^{4}\right]^{1/2}\mathbb{E}\left[\Delta^{2}_{h,j}\right]^{1/2}$
$\displaystyle\leq~{}C^{2}\mathbb{E}\left[\left(\boldsymbol{x}_{L}(1)-\boldsymbol{x}_{0}(1)\right)^{8}\right]^{1/4}\mathbb{E}\left[\left(|\boldsymbol{x}_{L}(1)|+|\boldsymbol{x}_{0}(1)|\right)^{8}\right]^{1/4}\mathbb{E}\left[\Delta^{2}_{h,j}\right]^{1/2}.$

The first two factors are finite by Lemma 4(i) and the third converges to zero
by Proposition 9.

Hence, as $\mathbb{E}\left[1\wedge e^{-\sum_{j=2}^{d}\Delta_{h,j}}\right]\to
a(\ell)$ by (i), we only need to show

$\mathbb{E}\left[\left(f(\boldsymbol{x}_{L}(1))-f(\boldsymbol{x}_{0}(1))\right)^{2}\right]\quad\xrightarrow{d\to\infty}\quad\mathbb{E}\left[\left(f(X_{T})-f(X_{0})\right)^{2}\right]\,.$

Note that for every dimension $d$ there exists a Langevin diffusion
$\boldsymbol{Z}_{t}$ linked to the numerical Langevin trajectory
$\boldsymbol{z}_{0:L}$ via (9). This ensures
$\boldsymbol{x}_{0}=\boldsymbol{X}_{0}$ and implies that

$\displaystyle\mathbb{E}\left[\left(f(\boldsymbol{x}_{L}(1))-f(\boldsymbol{x}_{0}(1))\right)^{2}\right]=$
$\displaystyle~{}\mathbb{E}\left[\left(f(\boldsymbol{x}_{L}(1))-f(\boldsymbol{X}_{T}(1))+f(\boldsymbol{X}_{T}(1))-f(\boldsymbol{X}_{0}(1))\right)^{2}\right]$
$\displaystyle=$
$\displaystyle~{}2\,\mathbb{E}\left[\left(f(X_{T})-f(X_{0})\right)\left(f(\boldsymbol{x}_{L}(1))-f(\boldsymbol{X}_{T}(1))\right)\right]$
$\displaystyle+$
$\displaystyle~{}\mathbb{E}\left[\left(f(X_{T})-f(X_{0})\right)^{2}\right]+\mathbb{E}\left[\left(f(\boldsymbol{x}_{L}(1))-f(\boldsymbol{X}_{T}(1))\right)^{2}\right]\,.$

We will show that the last term vanishes. By Cauchy’s inequality that means,
that the mixed term vanishes as well. By another Cauchy inequality and the
bound $|A-B|\leq|A|+|B|$

$\displaystyle\mathbb{E}\left[\left(f(\boldsymbol{x}_{L}(1))-f(\boldsymbol{X}_{T}(1))\right)^{2}\right]$
$\displaystyle\leq
C^{2}\mathbb{E}\left[\left(\boldsymbol{x}_{L}(1)-\boldsymbol{X}_{T}(1)\right)^{2}\left(\left|\boldsymbol{x}_{L}(1)\right|+\left|\boldsymbol{X}_{T}(1)\right|\right)^{2}\right]$
$\displaystyle\leq
C^{2}\mathbb{E}\left[\left|\boldsymbol{x}_{L}(1)-\boldsymbol{X}_{T}(1)\right|\left(\left|\boldsymbol{x}_{L}(1)\right|+\left|\boldsymbol{X}_{T}(1)\right|\right)^{3}\right]$
$\displaystyle\leq
C^{2}\Big{(}\mathbb{E}\big{[}\left(\boldsymbol{x}_{L}(1)-\boldsymbol{X}_{T}(1)\right)^{2}\big{]}\mathbb{E}\big{[}\left(\left|\boldsymbol{x}_{L}(1)\right|+\left|\boldsymbol{X}_{T}(1)\right|\right)^{6}\big{]}\Big{)}^{1/2}$

The first factor vanishes by Proposition 3 and the second is bounded by Lemma
4.

(iii) Again we synchronise the Langevin diffusion with the numerical Langevin
trajectory via (9), hence $\boldsymbol{X}^{n}(1)=\boldsymbol{x}_{0}(1)=X_{0}$.
We will prove that

$\mathbb{E}\left[f(\boldsymbol{X}^{n+1}(1))\right]\quad\xrightarrow{d\to\infty}\quad\mathbb{E}\left[f(X_{T}^{\prime})\right]$

for an arbitrary bounded Lipschitz function $f$. By Portmanteau theorem this
implies the result.

Note that since the function $t\mapsto 1\wedge e^{t}$ is $1$-Lipschitz,
Proposition 9 guarantees

$\mathbb{P}\left(U\text{ is between }1\wedge
e^{-\sum_{j=1}^{d}\Delta_{h,j}}\text{ and }1\wedge
e^{-\sum_{j=2}^{d}\Delta_{h,j}}\right)\\\ \leq\mathbb{E}\left[\left|1\wedge
e^{-\sum_{j=1}^{d}\Delta_{h,j}}-1\wedge
e^{-\sum_{j=2}^{d}\Delta_{h,j}}\right|\right]\leq\mathbb{E}\left[\left|\Delta_{h}\right|\right]\quad\xrightarrow{d\to\infty}\quad
0\,.$

Since $f$ is bounded this implies that

$\left|\mathbb{E}\left[f(\boldsymbol{x}_{L}(1))\mathds{1}_{U\leq 1\wedge
e^{-\sum_{j=1}^{d}\Delta_{h,j}}}\right]-\mathbb{E}\left[f(\boldsymbol{x}_{L}(1))\mathds{1}_{U\leq
1\wedge
e^{-\sum_{j=2}^{d}\Delta_{h,j}}}\right]\right|\quad\xrightarrow{d\to\infty}\quad
0$

and

$\left|\mathbb{E}\left[f(\boldsymbol{x}_{0}(1))\mathds{1}_{U>1\wedge
e^{-\sum_{j=1}^{d}\Delta_{h,j}}}\right]-\mathbb{E}\left[f(\boldsymbol{x}_{0}(1))\mathds{1}_{U>1\wedge
e^{-\sum_{j=2}^{d}\Delta_{h,j}}}\right]\right|\quad\xrightarrow{d\to\infty}\quad
0\,.$

Finally, this gives us

$\displaystyle\lim_{d\to\infty}\mathbb{E}\left[f(\boldsymbol{X}^{n+1}(1))\right]$
$\displaystyle=\lim_{d\to\infty}\left(\mathbb{E}\left[f(\boldsymbol{x}_{L}(1))\mathds{1}_{U\leq
1\wedge
e^{-\sum_{j=1}^{d}\Delta_{h,j}}}\right]+\mathbb{E}\left[f(\boldsymbol{x}_{0}(1))\mathds{1}_{U>1\wedge
e^{-\sum_{j=1}^{d}\Delta_{h,j}}}\right]\right)$
$\displaystyle=\lim_{d\to\infty}\left(\mathbb{E}\left[f(\boldsymbol{x}_{L}(1))\mathds{1}_{U\leq
1\wedge
e^{-\sum_{j=2}^{d}\Delta_{h,j}}}\right]+\mathbb{E}\left[f(\boldsymbol{x}_{0}(1))\mathds{1}_{U>1\wedge
e^{-\sum_{j=2}^{d}\Delta_{h,j}}}\right]\right)$
$\displaystyle=\lim_{d\to\infty}\left(\mathbb{E}\left[f(\boldsymbol{x}_{L}(1))\right]\mathbb{P}\left[U\leq
1\wedge
e^{-\sum_{j=2}^{d}\Delta_{h,j}}\right]+\mathbb{E}\left[f(\boldsymbol{x}_{0}(1))\right]\mathbb{P}\left[U>1\wedge
e^{-\sum_{j=2}^{d}\Delta_{h,j}}\right]\right)$
$\displaystyle=\mathbb{E}\left[f(X_{T})\right]a(\ell)+\mathbb{E}\left[f(X_{0})\right](1-a(\ell))~{}=~{}\mathbb{E}\left[f(X_{T}^{\prime})\right]\,.$

The first equality holds by construction of the accept-reject test, the third
by the product form assumption and the fourth by Proposition 3 and point (i).
∎

###### Proof of Proposition 7.

The proof relies on a similar argument as in the proof of [74, Theorem 3]. We
first establish that

$\displaystyle\begin{split}\mathbb{E}\left[\left(f(\boldsymbol{X}^{n+1}(1))-f(\boldsymbol{X}^{n}(1))\right)^{2}\right]~{}&=~{}\mathbb{E}\left[\left(f(\boldsymbol{x}_{L}(1))-f(\boldsymbol{x}_{0}(1))\right)^{2}1\wedge
e^{-\sum_{j=1}^{d}\Delta_{h,j}}\right]\\\
&\leq~{}2\mathbb{E}\left[\left(f(\boldsymbol{x}_{L}(1))-f(\boldsymbol{x}_{0}(1))\right)^{2}\right]\mathbb{E}\left[1\wedge
e^{-\sum_{j=2}^{d}\Delta_{h,j}}\right]\,.\end{split}$ (56)

The statement without the factor of two is apparent if $\Delta_{h,1}$ is
positive. If it is negative we use Proposition 10 and it is crucial that the
expression $(f(\boldsymbol{x}_{L}(1))-f(\boldsymbol{x}_{0}(1))^{2}$ is
symmetric with respect to the start $\boldsymbol{x}_{0}(1)$ and end
$\boldsymbol{x}_{L}(1)$ point of the trajectory:

$\displaystyle\mathbb{E}\left[\left(f(\boldsymbol{x}_{L}(1))-f(\boldsymbol{x}_{0}(1))\right)^{2}1\wedge
e^{-\sum_{j=1}^{d}\Delta_{h,j}}\mathds{1}_{(-\infty,0]}(\Delta_{h,1})\right]$
$\displaystyle\quad\leq~{}\mathbb{E}\left[\left(f(\boldsymbol{x}_{L}(1))-f(\boldsymbol{x}_{0}(1))\right)^{2}1\wedge
e^{-\sum_{j=2}^{d}\Delta_{h,j}}e^{-\Delta_{h,1}}\mathds{1}_{[0,\infty)}(-\Delta_{h,1})\right]$
$\displaystyle\quad=~{}\mathbb{E}\left[\left(f(\boldsymbol{x}_{0}(1))-f(\boldsymbol{x}_{L}(1))\right)^{2}1\wedge
e^{-\sum_{j=2}^{d}\Delta_{h,j}}\mathds{1}_{[0,\infty)}(\Delta_{h,1})\right]$
$\displaystyle\quad\leq~{}\mathbb{E}\left[\left(f(\boldsymbol{x}_{L}(1))-f(\boldsymbol{x}_{0}(1))\right)^{2}1\wedge
e^{-\sum_{j=2}^{d}\Delta_{h,j}}\right]\,.$

As $h\to 0$ we have
$\mathbb{E}[\left(f(\boldsymbol{x}_{L}(1))-f(\boldsymbol{x}_{0}(1))\right)^{2}]\to\Upsilon_{f}$;
see proof of Proposition 6(ii). Moreover, if $d^{1/4}h\to 0$, then
$d^{1/4}h\,\mathbb{E}[1\wedge\exp\\{-\sum_{j=2}^{d}\Delta_{h,j}\\}]\to 0$ as
$d$ increases, since the acceptance rate is bounded by one. Assume now
$d^{1/4}h\to\infty$. By [73, Theorem 7], Proposition 9 implies that

$\frac{\mathbb{E}\left[\Delta_{h}\right]}{\mathbb{E}\left[\Delta_{h}^{2}\right]}~{}\xrightarrow{h\to
0}~{}\frac{1}{2}\,.$

We split the probability space with respect to the event
$\mathcal{A}_{d}=\left\\{\sum_{j=2}^{d}\Delta_{h,j}\geq\tfrac{d-1}{2}\mathbb{E}[\Delta_{h}]\right\\}$.
Using the convergence above and Proposition 9 with $dh^{4}\to\infty$ we
obtain:

$\limsup_{d\to\infty}d^{1/4}h\mathbb{E}\left[1\wedge
e^{-\sum_{j=2}^{d}\Delta_{h,j}}\mathds{1}_{\mathcal{A}_{d}}\right]~{}\leq~{}\limsup_{d\to\infty}d^{1/4}he^{-\frac{d-1}{2}\mathbb{E}[\Delta_{h}]}\\\
\leq~{}\limsup_{d\to\infty}d^{1/4}he^{-\frac{d}{5}\mathbb{E}[\Delta^{2}_{h}]}=~{}\limsup_{d\to\infty}d^{1/4}he^{-\frac{dh^{4}}{5}h^{-4}\mathbb{E}[\Delta^{2}_{h}]}~{}\leq~{}\limsup_{d\to\infty}d^{1/4}he^{-dh^{4}\frac{\Sigma}{6}}~{}=~{}0\,.$

On the complement $\mathcal{A}_{d}^{c}$ we use Chebyshev’s inequality in
addition

$\displaystyle\limsup_{d\to\infty}$
$\displaystyle\,d^{1/4}h\,\mathbb{E}\left[1\wedge
e^{-\sum_{j=2}^{d}\Delta_{h,j}}\mathds{1}_{\mathcal{A}_{d}^{c}}\right]\leq\limsup_{d\to\infty}\,d^{1/4}h\,\mathbb{P}(\mathcal{A}_{d}^{c})$
$\displaystyle=\limsup_{d\to\infty}\,d^{1/4}h\,\mathbb{P}\bigg{(}\sum_{j=1}^{d}\left(\mathbb{E}[\Delta_{h,j}]-\Delta_{h,j}\right)\geq\frac{d-1}{2}\mathbb{E}[\Delta_{h}]\bigg{)}$
$\displaystyle\leq\limsup_{d\to\infty}\,d^{1/4}h\,\frac{4d\text{Var}[\Delta_{h}]}{d^{2}\mathbb{E}[\Delta_{h}]^{2}}\leq\limsup_{d\to\infty}\,d^{1/4}h\,\frac{17}{d\mathbb{E}[\Delta_{h}^{2}]}\leq\limsup_{d\to\infty}\frac{18}{\Sigma}\frac{1}{d^{3/4}h^{3}}=0\,.$

The last two conclusions together with (56) show the sub-optimality of scaling
of time-step $h$ different to $d^{-1/4}$. We conclude that $d^{-1/4}$ scaling
is optimal by noting that the claimed non-zero limit is achieved in the case
when $d^{1/4}h\to\ell$ by Proposition 6(ii).

Finally, to optimize over the choice of $\ell$ note that the function
$\ell\mapsto 2\ell\Psi\left(-\tfrac{1}{2}\ell^{2}\sqrt{\Sigma}\right)$ is
smooth and converges to zero both as $\ell\to 0$ and $\ell\to\infty$. Since it
is positive, its maximum must be attained at a stationary point. Substituting
$s=\ell^{2}\sqrt{\Sigma}/2$ we can find stationary points of $s\mapsto
2^{3/2}\Sigma^{-1/4}\sqrt{s}\Psi(-s)$ which is equivalent to finding solutions
of $\frac{1}{2}=s\frac{\psi(-s)}{\Psi(-s)}$, where $\psi$ is the density of
the standard Gaussian. There exists a unique solution $s^{*}$ since the
function $s\mapsto s\frac{\psi(-s)}{\Psi(-s)}$ is strictly increasing. Hence,
${\rm eff}(\ell)$ attains its greatest value at a specific value
$\ell^{*}=\sqrt{2s^{*}}\Sigma^{-1/4}$ that corresponds to an average
acceptance rate $2\Psi(-s^{*})$, which turns out to numerically equal $0.651$
to three decimal places. It also implies that

${\rm
eff}(\ell^{*})~{}=~{}2^{3/2}\sqrt{s^{*}}\Psi(-s^{*})\times\Sigma^{-1/4}~{}\approx~{}0.619219\times\Sigma^{-1/4}\,.$

∎

## References

  * [1] [author] Alamo, A.A. and Sanz-Serna, J. M.J. M. (2016). A technique for studying strong and weak local errors of splitting stochastic integrators. SIAM J. Numer. Anal. 54 3239–3257. 10.1137/16M1058765 3570281
  * [2] [author] Andrieu, ChristopheC., Lee, AnthonyA. and Livingstone, SamS. (2020). A general perspective on the Metropolis-Hastings kernel. Preprint. Available at arXiv:2012.14881.
  * [3] [author] Andrieu, ChristopheC. and Livingstone, SamuelS. (2021). Peskun-Tierney ordering for Markovian Monte Carlo: beyond the reversible scenario. Ann. Statist. 49 1958–1981. 10.1214/20-aos2008 4319237
  * [4] [author] Besag, JulianJ. (1994). Comments on “Representations of knowledge in complex systems” by U. Grenander and MI Miller. J. Roy. Statist. Soc. Ser. B 56 4\.
  * [5] [author] Beskos, AlexandrosA., Pillai, NateshN., Roberts, GarethG., Sanz-Serna, Jesus-MariaJ.-M. and Stuart, AndrewA. (2013). Optimal tuning of the hybrid Monte Carlo algorithm. Bernoulli 19 1501–1534. 10.3150/12-BEJ414 3129023
  * [6] [author] Betancourt, MichaelM. (2016). Identifying the optimal integration time in Hamiltonian Monte Carlo. Preprint. Available at arXiv:1601.00225.
  * [7] [author] Betancourt, MichaelM., Byrne, SimonS., Livingstone, SamS. and Girolami, MarkM. (2017). The geometric foundations of Hamiltonian Monte Carlo. Bernoulli 23 2257–2298. 10.3150/16-BEJ810 3648031
  * [8] [author] Bou-Rabee, NawafN. and Eberle, AndreasA. (2021). Mixing Time Guarantees for Unadjusted Hamiltonian Monte Carlo. Preprint. Available at arXiv:2105.00887.
  * [9] [author] Bou-Rabee, NawafN., Eberle, AndreasA. and Zimmer, RaphaelR. (2020). Coupling and convergence for Hamiltonian Monte Carlo. Ann. Appl. Probab. 30 1209–1250. 10.1214/19-AAP1528 4133372
  * [10] [author] Bou-Rabee, NawafN. and Owhadi, HoumanH. (2010). Long-run accuracy of variational integrators in the stochastic context. SIAM J. Numer. Anal. 48 278–297. 10.1137/090758842 2608370
  * [11] [author] Bou-Rabee, NawafN. and Sanz-Serna, Jesús MaríaJ. M. (2017). Randomized Hamiltonian Monte Carlo. Ann. Appl. Probab. 27 2159–2194. 10.1214/16-AAP1255 3693523
  * [12] [author] Bou-Rabee, NawafN. and Schuh, KatharinaK. (2020). Convergence of unadjusted Hamiltonian Monte Carlo for mean-field models. Preprint. Available at arXiv:2009.08735.
  * [13] [author] Bou-Rabee, NawafN. and Vanden-Eijnden, EricE. (2010). Pathwise accuracy and ergodicity of metropolized integrators for SDEs. Comm. Pure Appl. Math. 63 655–696. 10.1002/cpa.20306 2583309
  * [14] [author] Bussi, GiovanniG. and Parrinello, MicheleM. (2007). Accurate sampling using Langevin dynamics. Physical Review E 75 056707\.
  * [15] [author] Campos, Cédric M.C. M. and Sanz-Serna, J. M.J. M. (2015). Extra chance generalized hybrid Monte Carlo. J. Comput. Phys. 281 365–374. 10.1016/j.jcp.2014.09.037 3281978
  * [16] [author] Carpenter, BobB., Gelman, AndrewA., Hoffman, Matthew DM. D., Lee, DanielD., Goodrich, BenB., Betancourt, MichaelM., Brubaker, MarcusM., Guo, JiqiangJ., Li, PeterP. and Riddell, AllenA. (2017). Stan: A probabilistic programming language. Journal of statistical software 76 1–32.
  * [17] [author] Chen, YuansiY., Dwivedi, RaazR., Wainwright, Martin J.M. J. and Yu, BinB. (2020). Fast mixing of metropolized Hamiltonian Monte Carlo: benefits of multi-step gradients. J. Mach. Learn. Res. 21 Paper No. 92, 71. 4119160
  * [18] [author] Chen, ZongchenZ. and Vempala, Santosh S.S. S. (2019). Optimal convergence rate of Hamiltonian Monte Carlo for strongly logconcave distributions. 145 Art. No. 64, 12. 4012714
  * [19] [author] Cheng, XiangX., Chatterji, Niladri SN. S., Bartlett, Peter LP. L. and Jordan, Michael IM. I. (2018). Underdamped Langevin MCMC: A non-asymptotic analysis. In Conference on Learning Theory 300–323. PMLR.
  * [20] [author] Chewi, SinhoS., Lu, ChenC., Ahn, KwangjunK., Cheng, XiangX., Le Gouic, ThibautT. and Rigollet, PhilippeP. (2021). Optimal dimension dependence of the metropolis-adjusted langevin algorithm. In Conference on Learning Theory 1260–1300. PMLR.
  * [21] [author] Dalalyan, Arnak S.A. S. (2017). Theoretical guarantees for approximate sampling from smooth and log-concave densities. J. R. Stat. Soc. Ser. B. Stat. Methodol. 79 651–676. 10.1111/rssb.12183 3641401
  * [22] [author] Dalalyan, Arnak S.A. S. and Karagulyan, AvetikA. (2019). User-friendly guarantees for the Langevin Monte Carlo with inaccurate gradient. Stochastic Process. Appl. 129 5278–5311. 10.1016/j.spa.2019.02.016 4025705
  * [23] [author] Dalalyan, Arnak SA. S., Karagulyan, AvetikA. and Riou-Durand, LionelL. (2019). Bounding the error of discretized Langevin algorithms for non-strongly log-concave targets. Preprint. Available at arXiv:1906.08530.
  * [24] [author] Dalalyan, Arnak S.A. S. and Riou-Durand, LionelL. (2020). On sampling from a log-concave density using kinetic Langevin diffusions. Bernoulli 26 1956–1988. 10.3150/19-BEJ1178 4091098
  * [25] [author] Deligiannidis, GeorgeG., Paulin, DanielD., Bouchard-Côté, AlexandreA. and Doucet, ArnaudA. (2021). Randomized Hamiltonian Monte Carlo as scaling limit of the bouncy particle sampler and dimension-free convergence rates. Ann. Appl. Probab. 31 2612–2662. 10.1214/20-aap1659 4350970
  * [26] [author] Duane, SimonS., Kennedy, A. D.A. D., Pendleton, Brian J.B. J. and Roweth, DuncanD. (1987). Hybrid Monte Carlo. Phys. Lett. B 195 216–222. 10.1016/0370-2693(87)91197-x 3960671
  * [27] [author] Durmus, AlainA., Majewski, SzymonS. and Miasojedow, Bła˙zejB. a. z. (2019). Analysis of Langevin Monte Carlo via convex optimization. J. Mach. Learn. Res. 20 Paper No. 73, 46. 3960927
  * [28] [author] Durmus, AlainA. and Moulines, ÉricE. (2017). Nonasymptotic convergence analysis for the unadjusted Langevin algorithm. Ann. Appl. Probab. 27 1551–1587. 10.1214/16-AAP1238 3678479
  * [29] [author] Durmus, AlainA. and Moulines, ÉricE. (2019). High-dimensional Bayesian inference via the unadjusted Langevin algorithm. Bernoulli 25 2854–2882. 10.3150/18-BEJ1073 4003567
  * [30] [author] Durmus, AlainA., Moulines, ÉricE. and Saksman, EeroE. (2020). Irreducibility and geometric ergodicity of Hamiltonian Monte Carlo. Ann. Statist. 48 3545–3564. 10.1214/19-AOS1941 4185819
  * [31] [author] Dwivedi, RaazR., Chen, YuansiY., Wainwright, Martin J.M. J. and Yu, BinB. (2019). Log-concave sampling: Metropolis-Hastings algorithms are fast. J. Mach. Learn. Res. 20 Paper No. 183, 42. 4048994
  * [32] [author] Eberle, AndreasA., Guillin, ArnaudA. and Zimmer, RaphaelR. (2019). Couplings and quantitative contraction rates for Langevin dynamics. Ann. Probab. 47 1982–2010. 10.1214/18-AOP1299 3980913
  * [33] [author] Ethier, Stewart NS. N. and Kurtz, Thomas GT. G. (2009). Markov processes: characterization and convergence 282\. John Wiley & Sons.
  * [34] [author] Gardiner, C. W.C. W. (2004). Handbook of stochastic methods for physics, chemistry and the natural sciences, third ed. Springer Series in Synergetics 13\. Springer-Verlag, Berlin. 10.1007/978-3-662-05389-8 2053476
  * [35] [author] Green, Peter J.P. J. and Mira, AntoniettaA. (2001). Delayed rejection in reversible jump Metropolis-Hastings. Biometrika 88 1035–1053. 10.1093/biomet/88.4.1035 1872218
  * [36] [author] Hastings, W. K.W. K. (1970). Monte Carlo sampling methods using Markov chains and their applications. Biometrika 57 97–109. 10.1093/biomet/57.1.97 3363437
  * [37] [author] Hoffman, MatthewM., Radul, AlexeyA. and Sountsov, PavelP. (2021). An Adaptive-MCMC Scheme for Setting Trajectory Lengths in Hamiltonian Monte Carlo. In International Conference on Artificial Intelligence and Statistics 3907–3915. PMLR.
  * [38] [author] Hoffman, Matthew D.M. D. and Gelman, AndrewA. (2014). The no-U-turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. J. Mach. Learn. Res. 15 1593–1623. 3214779
  * [39] [author] Horowitz, Alan MA. M. (1991). A generalized guided Monte Carlo algorithm. Physics Letters B 268 247–252.
  * [40] [author] Karagulyan, AvetikA. and Dalalyan, ArnakA. (2020). Penalized Langevin dynamics with vanishing penalty for smooth and log-concave targets. Advances in Neural Information Processing Systems 33 17594–17604.
  * [41] [author] Kennedy, A. D.A. D. and Pendleton, BrianB. (2001). Cost of the generalised hybrid Monte Carlo algorithm for free field theory. Nuclear Phys. B 607 456–510. 10.1016/S0550-3213(01)00129-8 1850796
  * [42] [author] Lee, Yin TatY. T., Shen, RuoqiR. and Tian, KevinK. (2020). Logsmooth gradient concentration and tighter runtimes for metropolized hamiltonian monte carlo. In Conference on Learning Theory 2565–2597. PMLR.
  * [43] [author] Leimkuhler, BenedictB. and Matthews, CharlesC. (2013). Rational construction of stochastic numerical methods for molecular sampling. Appl. Math. Res. Express. AMRX 1 34–56. 10.1093/amrx/abs010 3040887
  * [44] [author] Livingstone, SamuelS., Betancourt, MichaelM., Byrne, SimonS. and Girolami, MarkM. (2019). On the geometric ergodicity of Hamiltonian Monte Carlo. Bernoulli 25 3109–3138. 10.3150/18-BEJ1083 4003576
  * [45] [author] Livingstone, SamuelS. and Zanella, GiacomoG. (2019). The Barker proposal: combining robustness and efficiency in gradient-based MCMC. Preprint. Available at arXiv:1908.11812.
  * [46] [author] Ma, Yi-AnY.-A., Chatterji, Niladri S.N. S., Cheng, XiangX., Flammarion, NicolasN., Bartlett, Peter L.P. L. and Jordan, Michael I.M. I. (2021). Is there an analog of Nesterov acceleration for gradient-based MCMC? Bernoulli 27 1942–1992. 10.3150/20-bej1297 4278799
  * [47] [author] Mangoubi, OrenO. and Smith, AaronA. (2019). Mixing of Hamiltonian Monte Carlo on strongly log-concave distributions 2: Numerical integrators. In The 22nd international conference on artificial intelligence and statistics 586–595. PMLR.
  * [48] [author] Mangoubi, OrenO. and Smith, AaronA. (2021). Mixing of Hamiltonian Monte Carlo on strongly log-concave distributions: continuous dynamics. Ann. Appl. Probab. 31 2019–2045. 10.1214/20-aap1640 4332690
  * [49] [author] Mangoubi, OrenO. and Vishnoi, Nisheeth KN. K. (2018). Dimensionally tight bounds for second-order Hamiltonian Monte Carlo. arXiv preprint arXiv:1802.08898. Preprint. Available at arXiv:1802.08898.
  * [50] [author] Mattingly, J. C.J. C., Stuart, A. M.A. M. and Higham, D. J.D. J. (2002). Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise. Stochastic Process. Appl. 101 185–232. 10.1016/S0304-4149(02)00150-3 1931266
  * [51] [author] Metropolis, NicholasN., Rosenbluth, Arianna WA. W., Rosenbluth, Marshall NM. N., Teller, Augusta HA. H. and Teller, EdwardE. (1953). Equation of state calculations by fast computing machines. The journal of chemical physics 21 1087–1092.
  * [52] [author] Meyn, Sean P.S. P. and Tweedie, R. L.R. L. (1993). Stability of Markovian processes III: Foster-Lyapunov criteria for continuous-time processes. Advances in Applied Probability 25 518–548. 10.2307/1427522 1234295
  * [53] [author] Milstein, G. N.G. N. and Tretyakov, M. V.M. V. (2004). Stochastic numerics for mathematical physics. Scientific Computation. Springer-Verlag, Berlin. 10.1007/978-3-662-10063-9 2069903
  * [54] [author] Mira, AntoniettaA. (2001). On Metropolis-Hastings algorithms with delayed rejection. Metron 59 231–241. 1889712
  * [55] [author] Monmarché, PierreP. (2021). High-dimensional MCMC with a standard splitting scheme for the underdamped Langevin diffusion. Electron. J. Stat. 15 4117–4166. 10.1214/21-ejs1888 4309974
  * [56] [author] Neal, Radford M.R. M. (2011). MCMC using Hamiltonian dynamics. 113–162. 2858447
  * [57] [author] Nelson, EdwardE. (1967). Dynamical theories of Brownian motion. Princeton University Press, Princeton, N.J. 0214150
  * [58] [author] Ollivier, YannY. (2009). Ricci curvature of Markov chains on metric spaces. J. Funct. Anal. 256 810–864. 10.1016/j.jfa.2008.11.001 2484937
  * [59] [author] Ottobre, MichelaM., Pillai, Natesh S.N. S., Pinski, Frank J.F. J. and Stuart, Andrew M.A. M. (2016). A function space HMC algorithm with second order Langevin diffusion limit. Bernoulli 22 60–106. 10.3150/14-BEJ621 3449777
  * [60] [author] Park, JoonhaJ. and Atchadé, YvesY. (2020). Markov chain Monte Carlo algorithms with sequential proposals. Stat. Comput. 30 1325–1345. 10.1007/s11222-020-09948-4 4137254
  * [61] [author] Pasarica, CristianC. and Gelman, AndrewA. (2010). Adaptively scaling the metropolis algorithm using expected squared jumped distance. Statist. Sinica 20 343–364. 2640698
  * [62] [author] Pillai, Natesh S.N. S., Stuart, Andrew M.A. M. and Thiéry, Alexandre H.A. H. (2012). Optimal scaling and diffusion limits for the Langevin algorithm in high dimensions. Ann. Appl. Probab. 22 2320–2356. 10.1214/11-AAP828 3024970
  * [63] [author] Ricci, AndreaA. and Ciccotti, GiovanniG. (2003). Algorithms for Brownian dynamics. Molecular Physics 101 1927–1931.
  * [64] [author] Roberts, G. O.G. O., Gelman, A.A. and Gilks, W. R.W. R. (1997). Weak convergence and optimal scaling of random walk Metropolis algorithms. Ann. Appl. Probab. 7 110–120. 10.1214/aoap/1034625254 1428751
  * [65] [author] Roberts, Gareth O.G. O. and Rosenthal, Jeffrey S.J. S. (1998). Optimal scaling of discrete approximations to Langevin diffusions. J. R. Stat. Soc. Ser. B Stat. Methodol. 60 255–268. 10.1111/1467-9868.00123 1625691
  * [66] [author] Roberts, Gareth O.G. O. and Rosenthal, Jeffrey S.J. S. (2004). General state space Markov chains and MCMC algorithms. Probab. Surv. 1 20–71. 10.1214/154957804100000024 2095565
  * [67] [author] Roberts, Gareth O.G. O. and Tweedie, Richard L.R. L. (1996). Exponential convergence of Langevin distributions and their discrete approximations. Bernoulli 2 341–363. 10.2307/3318418 1440273
  * [68] [author] Sanz-Serna, Jesus MaríaJ. M. and Zygalakis, Konstantinos C.K. C. (2021). Wasserstein distance estimates for the distributions of numerical approximations to ergodic stochastic differential equations. J. Mach. Learn. Res. 22 Paper No. 242, 37. 10.1080/14685248.2020.1855352 4329821
  * [69] [author] Scemama, AnthonyA., Lelièvre, TonyT., Stoltz, GabrielG., Cancès, EricE. and Caffarel, MichelM. (2006). An efficient sampling algorithm for variational Monte Carlo. The Journal of chemical physics 125 114105\.
  * [70] [author] Sherlock, ChrisC., Urbas, SzymonS. and Ludkin, MatthewM. (2021). Apogee to Apogee Path Sampler. Preprint. Available at arXiv:2112.08187.
  * [71] [author] Sohl-Dickstein, JaschaJ., Mudigonda, MayurM. and DeWeese, MichaelM. (2014). Hamiltonian Monte Carlo without detailed balance. In International Conference on Machine Learning 719–726. PMLR.
  * [72] [author] Tierney, LukeL. and Mira, AntoniettaA. (1999). Some adaptive Monte Carlo methods for Bayesian inference. Statistics in medicine 18 2507–2515.
  * [73] [author] Vogrinc, JureJ. and Kendall, Wilfrid S.W. S. (2021). Counterexamples for optimal scaling of Metropolis-Hastings chains with rough target densities. Ann. Appl. Probab. 31 972–1019. 10.1214/20-aap1612 4254501
  * [74] [author] Vogrinc, JureJ., Livingstone, SamuelS. and Zanella, GiacomoG. (2022). Optimal design of the Barker proposal and other locally-balanced Metropolis-Hastings algorithms. Preprint. Available at arXiv:2201.01123.
  * [75] [author] Wagoner, Jason AJ. A. and Pande, Vijay SV. S. (2012). Reducing the effect of Metropolization on mixing times in molecular dynamics simulations. The Journal of chemical physics 137 214105\.
  * [76] [author] Wang, ZiyuZ., Mohamed, ShakirS. and Freitas, NandoN. (2013). Adaptive hamiltonian and riemann manifold monte carlo. In International conference on machine learning 1462–1470. PMLR.


11institutetext: 1Center for Intelligent Machines, McGill University, Canada.
2MILA (Quebec AI institute), Canada.
3Montreal Neurological Institute, McGill University, Canada.
4Institute for Digital Communications, School of Engineering, University of
Edinburgh, UK.
5The Alan Turing Institute, UK
11email<EMAIL_ADDRESS>
# Debiasing Counterfactuals In the Presence of Spurious Correlations

Amar Kumar1,2 Nima Fathi1,2 Raghav Mehta1,2 Brennan Nichyporuk1,2 Jean-Pierre
R. Falet2,3 Sotirios Tsaftaris 4,5 Tal Arbel1,2

###### Abstract

Deep learning models can perform well in complex medical imaging
classification tasks, even when basing their conclusions on spurious
correlations (i.e. confounders), should they be prevalent in the training
dataset, rather than on the causal image markers of interest. This would
thereby limit their ability to generalize across the population.
Explainability based on counterfactual image generation can be used to expose
the confounders but does not provide a strategy to mitigate the bias. In this
work, we introduce the first end-to-end training framework that integrates
both (i) popular debiasing classifiers (e.g. distributionally robust
optimization (DRO)) to avoid latching onto the spurious correlations and (ii)
counterfactual image generation to unveil generalizable imaging markers of
relevance to the task. Additionally, we propose a novel metric, Spurious
Correlation Latching Score (SCLS), to quantify the extent of the classifier
reliance on the spurious correlation as exposed by the counterfactual images.
Through comprehensive experiments on two public datasets (with the simulated
and real visual artifacts), we demonstrate that the debiasing method: (i)
learns generalizable markers across the population, and (ii) successfully
ignores spurious correlations and focuses on the underlying disease pathology.

###### Keywords:

Biomarker Counterfactuals Debiasing Explainablity

##  1 Introduction

Deep learning models have shown tremendous success in disease
classificationbased on medical images, given their ability to learn complex
imaging markers across a wide population of subjects. These models can show
good performance and still be biased as they may focus on spurious
correlations in the image that are not causally related to the disease but
arise due to confounding factors - should they be common across the majority
of samples in the training dataset. As a result, the confounding predictive
image markers may not generalize across the population. For example, a deep
learning model was able to accurately detect COVID-19 from chest radiographs,
but rather than relying on pathological evidence, the model latched on to
spurious correlations such as medical devices or lettering in the image [3].
As a result, these image markers did not generalize across the population.

In order to safely deploy black-box deep learning models in real clinical
applications, explainability should be integrated into the framework so as to
expose the spurious correlations on which the classifier based its
conclusions. Popular post-hoc explainability strategies, such as Grad-CAM [16,
6, 19], SHAP [10], LIME [11] are not designed to expose the precise predictive
markers driving a classifier. Models that integrate counterfactual image
generation, along with black-box classifers [21, 2, 23], permit exposing the
predictive markers used by the classifier. However, should these methods
discover that the markers are indeed simply visual artifacts there are no
strategies to mitigate the resulting biases. Furthermore, although several
debiasing methods have been successfully implemented to account for
generalizability [1, 26, 17, 8, 27], they do not integrate explainability into
the framework in order to provide reasons for improved performance.

Therefore, the important question to be answered is - Can a model be trained
to disregard spurious correlations and identify generalizable predictive
disease markers?

Figure 1: Counterfactual (CF) image indicating that the classifier latched
onto spurious correlations (medical devices) when correctly predicting that
subject is sick (class: Pleural Effusion), due to their prevalence in the
training dataset for this class. (a) Chest radiograph of a sick subject with
several medical devices shown (cyan boxes), (b) Generated (CF) image, (c)
Difference heat map shows maximum change around the medical devices, rather
than indicating the correct markers for the disease.

In this paper, we propose the first end-to-end training framework for the
explainability of classifier and debiasing via counterfactual image
generation. We seek to discover imaging markers that reflect underlying
disease pathology and that generalize across subgroups. Extensive experiments
are performed on two different publicly available datasets - (i) RSNA
Pneumonia Detection Challenge and (ii) CheXpert [5]. To illustrate the goal,
Figure 1 shows an example from the contrived CheXpert dataset, where most of
the sick subjects have medical device(s) (e.g. a pacemaker) in their images
while most of the healthy subjects do not. As such, there exists a spurious
correlation between a confounding visual artifact (the medical devices) and
the disease. A classifier based on a standard optimization technique,
empirical risk minimization (ERM), incorrectly indicates the medical device as
a disease marker, as depicted by the counterfactual (CF). In this work, we
propose replacing ERM with a popular debiasing method, Group-DRO
(distributional robust optimization). This permits the classifier to focus on
the pathological image markers of the disease rather than on spurious
correlation(s). Additionally, we show that Group-DRO ignores the visual
artifact when making its decision, and generalizes across subgroups without
the spurious correlation. Since standard metrics to evaluate counterfactuals
do not indicate the region where the classifier focuses, we also propose a
novel metric, the Spurious Correlation Latching Score (SCLS), to measure the
degree to which the classifier latches onto spurious correlations. Our
experiments indicate an improvement (in terms of differences in classifier
outputs) of 0.68 and 0.54 in the SCLS using the Group-DRO classifier over the
ERM for each of the two datasets.

##  2 Methodology

We propose an end-to-end training strategy to explain the output of a
classifier. Here, we are considering a scenario where majority of the training
data encompasses a spurious correlation with the target label. However, there
is also a minority subgroup in the dataset that does not have any spurious
correlation with the target label i.e., if the classifier was to rely onto the
spurious correlation then the performance on these minority subgroups will be
poor. Also, the term ‘majority’ and ‘minority’ is based on the number of
samples in these groups. An overview of our approach is shown in Figure 2.

###  2.1 Classifier Explainability and Debiasing Via Counterfactual Image
Generation

Disease Classification Binary (e.g. "sick" or "healthy") classification of the
images is performed using either a standard classifier (ERM [24]), or a
classifier that mitigates biases across sub-groups (Group-DRO [18]). The ERM
classifier ($f_{ERM}$) is expected to be affected by the spurious correlation
present in the training dataset, as it minimizes the loss over the entire
training dataset and latching onto spurious correlation is a shortcut to
minimize the loss. Thus, it would not generalize across the minority subgroups
of the dataset [12, 20]. On the contrary, the DRO classifier ($f_{DRO}$) is
not expected to learn the spurious correlation as it considers the majority
and minority subgroups separately when optimizing the loss. Thus, it would
generalize well across all subgroups.

Generative model for synthesizing counterfactuals We develop an explainability
framework that integrates counterfactual image generation together with a
classifier during training. We adapted Cycle-GAN [25] as the generative model
for counterfactual image generation, chosen for its strong performance across
a variety of domains [13, 25]. A pre-trained, frozen binary classifier
($f_{ERM}$ or $f_{DRO}$) provides supervision to the generator. The proposed
architecture and optimization objectives (see Figure 2) are designed to
generate counterfactual images that adhere to the following common constraints
[14, 15, 7]: (i)Identity preservation: The counterfactual images resemble the
input images with minimal change; (ii) Classifier consistency: Counterfactual
images belong to the target class; (iii)Cycle consistency: When counterfactual
images are fed through the opposing generator, the output reverts to the
original image (see Figure 2).

During inference, based on the classifier’s decision (i.e., $f_{ERM}\text{ or
}f_{DRO}$) for the input image, we generate counterfactual images and analyze
the difference heatmap between the factual (input) and counterfactual
(synthesized) images. This interpretable heat map indicates the image markers
that contribute the most to changing the classifier’s decision.

Figure 2: Training procedure overview: The black-box classifier can be
$f_{ERM}$ or $f_{DRO}$ and provides supervision to maintain the correct target
class, $y_{t}$. Two U-Net generators, $G_{\pazocal{S}\pazocal{H}}$ and
$G_{\pazocal{H}\pazocal{S}}$, are employed to synthesize counterfactual
images, namely $x_{\pazocal{S}_{cf}}$ and $x_{\pazocal{H}_{cf}}$. The
discriminator $D_{\pazocal{H}}$ and $D_{\pazocal{S}}$ compares the
counterfactual images with the domain of healthy $\pazocal{H}$ and sick
$\pazocal{S}$ subjects respectively. Note, training a cycle-GAN requires
simultaneous use of two input images from the two distributions.

###  2.2 Metrics for Evaluating Counterfactuals: Accounting for Spurious
Correlations

Standard counterfactual evaluation metrics are structured so as to ensure that
the generated images (a) preserve the subject identity and thus penalize
generated counterfactual images that are significantly different from the
factual (original) images and (b) result in a maximal change in the class
label (e.g. from healthy to sick). Identity preservation is typically measured
by structural similarity index  (SSIM) [4] and Actionability [15, 14], defined
as $\mathbb{E}\left[\left\|x-x_{cf}\right\|_{L_{1}}\right]$ between factual
($x$) and counterfactual ($x_{cf}$) images. Here, a higher value for SSIM and
a lower value for Actionability would indicate better counterfactuals. The
counterfactual prediction gain (CPG) [15], defined as $|f(x)-f(x_{cf})|$,
indicates the degree of change in the classifier’s decision such that a higher
value of CPG indicates better counterfactuals.

While such metrics are required to measure the validity of the generated
counterfactuals, they do not assess whether the classifier latched onto
spurious correlations. For example, consider an image of a sick subject in the
presence of a spurious correlation. If the disease classifier, $f_{ERM}$,
latched onto the spurious correlation when identifying the subject as sick,
the corresponding counterfactual image (i.e., depicting a healthy subject)
would show changes in the area of the spurious correlation. In this case, all
three evaluation metrics mentioned above would determine that this is a valid
counterfactual image, based on high SSIM and low Actionability (shows minimal
changes made compared to the factual image) and high CPG (due to the
classifier decision changing from sick to healthy). However, the
counterfactual image shows changes in the area of the spurious correlation
rather than depicting the correct predictive image markers for the disease as
desired.

In order to indicate that the classifier is correct but for the wrong reasons,
we introduce a novel metric called Spurious Correlation Latching Score (SCLS)
defined as follows:

$\text{SCLS}=|d(x)-d(x_{cf})|.$ (1)

Here, $d(\cdot)$ is a separate classifier, trained to identify the presence of
spurious correlation in the image. In cases where the counterfactual image
makes changes in an area of spurious correlation, SCLS will be high, as the
$d(\cdot)$ will show a maximum change in its prediction between factual and
counterfactual images. On the other hand, if the counterfactual image does not
make changes in the area of the spurious correlation then SCLS will have a low
value. As such, this evaluation strategy will validate how well the
counterfactuals can help to determine that the classifier latched onto
spurious correlations.

##  3 Experiments and Results

###  3.1 Dataset and Implementation Details

We perform experiments on two publicly available datasets. The absence of
ground truth makes the validation of counterfactual images particularly
challenging. Therefore, to directly evaluate the quality of the generated
counterfactual images in the presence of spurious correlations, we modify a
publicly available dataset (RSNA Pneumonia Detection Challenge) by adding a
synthetic artifact to the majority of the sick images (90%). The majority of
the sick and few of the healthy subjects have an artifact in the image,
whereas the majority of the healthy and a few sick subjects do not have this
artifact. The spurious correlation (artifact) is a black dot of radius 9
pixels at the center of the image. Thus, there are a total of four subgroups
($majority_{S}$, $majority_{H}$, $minority_{S}$ and $minority_{H}$) in the
dataset with varying number of images: $majority_{S}$ and $majority_{H}$ are
majority subgroups (sick with artifact and healthy without artifact), while
$minority_{S}$ and $minority_{H}$ are minority subgroups (sick without
artifact and healthy with artifact). Henceforth, this dataset will be referred
to as Dataset 1.

| Disease | Image size | Classifier |  | # samples
---
[$majority_{S}$, $minority_{S}$,
$minority_{H}$, $majority_{H}$]
Dataset 1 | Pneumonia | 512 x 512 |  AlexNet [13] | 5413, 1526, 883, 7968
Dataset 2 |  | Pleural
---
Effusion
224 x 224 |  |  Resnet-50 [22]
---
(pre-trained)
2600, 260, 350, 3456
Table 1: Implementation details for the two datasets

We also show experiments on a subset of a publicly available dataset (CheXpert
[5]) with medical devices (visual artifacts), spuriously correlated with the
disease. Specifically, we extract the subset of images that have labels
“healthy” or “pleural effusion” (subjects with the presence of other diseases
are removed from the dataset). This dataset will be referred to as Dataset 2.
More details about both datasets are provided in Table 1. Note that both the
datasets are divided into training/validation/testing with 70/10/20 random
split. Example images for both datasets and all four subgroups are shown in
Figure 3.

Figure 3: Datasets 1 and 2 group division: The majority of the sick subjects
[$majority_{S}$] and the minority of healthy subjects [$minority_{H}$] have
visual artifacts (shown in cyan boxes). The majority of healthy subjects
[$majority_{H}$] and the minority of sick subjects [$minority_{S}$] do not
have visual artifacts. Top row: Simulated artifacts (black dots); Bottom row:
Real artifacts (medical devices).

###  3.2 Results

Classifier Evaluation For both datasets (Figure 4), the DRO-based classifier
($f_{DRO}$) performs better for the minority subgroups ($minority_{S}$ and
$minority_{H}$); indicating that it can better generalize to sub-populations
that do not have the same visual artifact as the majority subgroups. Both
classifiers perform similarly for the majority subgroups ($majority_{S}$ and
$majority_{H}$).

Figure 4: Performance of ERM ($f_{ERM}$) and DRO ($f_{DRO}$) based classifier
on a held out test set across all subgroups for both datasets. Notice that DRO
has improved performance on minority subgroup [$minority_{S}$ and
$minority_{H}$] showing improved generalizability across all subgroups.

Qualitative Counterfactual Evaluation Pneumonia in chest radiograph manifests
as increased brightness in some regions of the lungs. In dataset 1, when
examining the majority subgroup of sick subjects, the ERM-based classifier
latches onto the spurious correlation, as seen by the difference maps. On the
other hand, a DRO-based classifier focuses on the pathology of the disease,
indicated by darker intensity regions over the lungs, as shown in Figure 5.
The behavior of $f_{ERM}$ is also evident in the minority subgroup, where the
counterfactual for a healthy subject exhibits an enlarged artifact, wrongly
suggesting that the visual artifact serves as a disease marker. Pleural
effusion is characterized by the rounding of the costophrenic angle, augmented
lung opacity, and reduced clarity of the diaphragm and lung fissures [9]. For
the majority subgroup of sick subjects in Dataset 2, the counterfactual images
based on ERM remove the medical device rather than focusing on the disease. In
addition, for healthy subjects from the minority subgroup, maximum changes are
observed around the medical device. On the other hand, for the majority
subgroup, the DRO-based counterfactuals show changes around the expected areas
while preserving the medical device.

Figure 5: Qualitative comparison of counterfactual (CF) images generated with
ERM and DRO classifiers for both majority (top row) and minority (bottom row)
subgroups. The ERM-based CFs show significant changes in the areas of spurious
correlation (cyan boxes), whereas the DRO-based CFs show almost no changes in
the same areas. In contrast, significant changes can be seen in the expected
area of disease pathology (magenta boxes) in DRO-based CFs, while the ERM-
based CFs show little to no changes in these areas.

Quantitative Counterfactual Evaluation In Table 2, counterfactual images
generated by ERM and DRO show similar scores according to standard metrics:
SSIM, Actionability and CPG. As these metrics are not designed to quantify
whether the generated counterfactuals are affected by spurious correlations
(see Section 2.2), the quality of the counterfactuals is now examined based on
the proposed SCLS metric. The AUC of the classifier, $d$, trained to detect
the presence of artifacts is 1.0 and 0.82 for Dataset 1 and Dataset 2,
respectively. As indicated by the last row of Table 2, the ERM-based
classifier shows a high value (poor performance) for SCLS for both datasets.
On the other hand, the DRO-based classifier has a low value (good performance)
for SCLS for both datasets. These results corroborate the finding made by
visual comparison of the counterfactual images generated by the ERM and DRO
classifiers. Overall, both qualitative and quantitative evaluations indicate
that an ERM-optimized classifier latches on to the spurious correlation
prevalent in the dataset, while a DRO-optimized classifier can be trained to
successfully ignore the spurious correlation.

| Dataset 1 | Dataset 2
---|---|---
| ERM | DRO | ERM | DRO
Actionability $\downarrow$ | 7.68 $\pm$ 0.01 | 7.86 $\pm$ 0.01 | 4.93 $\pm$ 0.01 | 5.68 $\pm$ 0.04
SSIM $\uparrow$ | 98.03 $\pm$ 0.00 | 98.44 $\pm$ 0.01 | 98.21 $\pm$ 0.01 | 98.36 $\pm$ 0.01
CPG $\uparrow$ | 0.91 $\pm$ 0.04 | 0.96 $\pm$0.03 | 0.88 $\pm$ 0.07 | 0.89 $\pm$ 0.04
SCLS $\downarrow$ | 0.80 $\pm$ 0.08 | $\mathbf{0.12\pm 0.07}$ | 0.76 $\pm$ 0.09 | $\mathbf{0.22\pm 0.06}$
Table 2: Quantitative results to compare counterfactual images generated for
both datasets. A low SCLS value implies that the model ($f_{DRO}$ in this
case) did not latch onto the spurious correlation.

##  4 Conclusion

Safe deployment of black-box models requires explainability to disclose when
the classifier is basing its predictions on spurious correlations and is
therefore not generalizable. In this paper, we presented the first integrated
end-to-end training strategy for generating unbiased counterfactual images,
capitalizing on a DRO classifier to enhance generalization. Our experiments
based on two datasets demonstrate that, unlike standard ERM classifiers which
are susceptible to latching onto spurious correlations, the unbiased DRO
classifier performs significantly better for minority subgroups in terms of-
(a) the classifier performance and (b) the novel SCLS metric, which quantifies
the degree to which the classifier latches on to the spurious correlation as
depicted by the generated counterfactual images.

Current datasets typically do not provide the ground truth predictive markers
of interest. Future work will require localizing the predictive markers (e.g.
with bounding boxes) and determining the degree of overlap with the discovered
markers. Further, we intend to explore the power of alternative debiasing
techniques and their potential contribution to discovering generalizable image
markers.

####  4.0.1 Acknowledgements

The authors are grateful for funding provided by the Natural Sciences and
Engineering Research Council of Canada, the Canadian Institute for Advanced
Research (CIFAR) Artificial Intelligence Chairs program, the Mila - Quebec AI
Institute technology transfer program, Microsoft Research, Calcul Quebec, and
the Digital Research Alliance of Canada. S.A. Tsaftaris acknowledges the
support of Canon Medical and the Royal Academy of Engineering and the Research
Chairs and Senior Research Fellowships scheme (grant RCSRF1819 / 8 / 25), and
the UK’s Engineering and Physical Sciences Research Council (EPSRC) support
via grant EP/X017680/1.

## References

  * [1] Burlina, P., Joshi, N., Paul, W., Pacheco, K.D., Bressler, N.M.: Addressing artificial intelligence bias in retinal diagnostics. Translational Vision Science & Technology 10(2), 13–13 (2021)
  * [2] Cohen, J.P., Brooks, R., En, S., Zucker, E., Pareek, A., Lungren, M.P., Chaudhari, A.: Gifsplanation via latent shift: a simple autoencoder approach to counterfactual generation for chest X-rays. In: Medical Imaging with Deep Learning. pp. 74–104. PMLR (2021)
  * [3] DeGrave, A.J., Janizek, J.D., Lee, S.I.: AI for radiographic COVID-19 detection selects shortcuts over signal. Nature Machine Intelligence 3(7), 610–619 (2021)
  * [4] Hore, A., Ziou, D.: Image quality metrics: PSNR vs. SSIM. In: 2010 20th international conference on pattern recognition. pp. 2366–2369. IEEE (2010)
  * [5] Irvin, J., Rajpurkar, P., Ko, M., Yu, Y., Ciurea-Ilcus, S., Chute, C., Marklund, H., Haghgoo, B., Ball, R., Shpanskaya, K., et al.: CheXpert: A large chest radiograph dataset with uncertainty labels and expert comparison. In: Proceedings of the AAAI conference on artificial intelligence. vol. 33, pp. 590–597 (2019)
  * [6] Jiang, H., Xu, J., Shi, R., Yang, K., Zhang, D., Gao, M., Ma, H., Qian, W.: A multi-label deep learning model with interpretable grad-cam for diabetic retinopathy classification. In: 2020 42nd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC). pp. 1560–1563. IEEE (2020)
  * [7] Kumar, A., Hu, A., Nichyporuk, B., Falet, J.P.R., Arnold, D.L., Tsaftaris, S., Arbel, T.: Counterfactual image synthesis for discovery of personalized predictive image markers. In: Artificial Intelligence over Infrared Images for Medical Applications and Medical Image Assisted Biomarker Discovery: First MICCAI Workshop, AIIIMA 2022, and First MICCAI Workshop, MIABID 2022, Held in Conjunction with MICCAI 2022, Singapore, September 18 and 22, 2022, Proceedings. pp. 113–124. Springer (2022)
  * [8] Larrazabal, A.J., Nieto, N., Peterson, V., Milone, D.H., Ferrante, E.: Gender imbalance in medical imaging datasets produces biased classifiers for computer-aided diagnosis. Proceedings of the National Academy of Sciences 117(23), 12592–12594 (2020)
  * [9] Light, R.W.: Pleural effusion. New England Journal of Medicine 346(25), 1971–1977 (2002)
  * [10] Lundberg, S.M., Lee, S.I.: A unified approach to interpreting model predictions. Advances in neural information processing systems 30 (2017)
  * [11] Magesh, P.R., Myloth, R.D., Tom, R.J.: An explainable machine learning model for early detection of parkinson’s disease using lime on datscan imagery. Computers in Biology and Medicine 126, 104041 (2020)
  * [12] Mehta, R., Shui, C., Arbel, T.: Evaluating the fairness of deep learning uncertainty estimates in medical image analysis. In: Medical Imaging with Deep Learning (2023)
  * [13] Mertes, S., Huber, T., Weitz, K., Heimerl, A., André, E.: GANterfactual—counterfactual explanations for medical non-experts using generative adversarial learning. Frontiers in artificial intelligence 5, 825565 (2022)
  * [14] Mothilal, R.K., Sharma, A., Tan, C.: Explaining machine learning classifiers through diverse counterfactual explanations. In: Proceedings of the 2020 conference on fairness, accountability, and transparency. pp. 607–617 (2020)
  * [15] Nemirovsky, D., Thiebaut, N., Xu, Y., Gupta, A.: CounteRGAN: Generating realistic counterfactuals with residual generative adversarial nets. arXiv preprint arXiv:2009.05199 (2020)
  * [16] Panwar, H., Gupta, P., Siddiqui, M.K., Morales-Menendez, R., Bhardwaj, P., Singh, V.: A deep learning and grad-CAM based color visualization approach for fast detection of COVID-19 cases using chest X-ray and CT-scan images. Chaos, Solitons & Fractals 140, 110190 (2020)
  * [17] Ricci Lara, M.A., Echeveste, R., Ferrante, E.: Addressing fairness in artificial intelligence for medical imaging. Nature communications 13(1), 4581 (2022)
  * [18] Sagawa, S., Koh, P.W., Hashimoto, T.B., Liang, P.: Distributionally robust neural networks for group shifts: On the importance of regularization for worst-case generalization. In: International Conference on Learning Representations (2019)
  * [19] Selvaraju, R.R., Cogswell, M., Das, A., Vedantam, R., Parikh, D., Batra, D.: Grad-CAM: Visual explanations from deep networks via gradient-based localization. In: Proceedings of the IEEE international conference on computer vision. pp. 618–626 (2017)
  * [20] Shui, C., Szeto, J., Mehta, R., Arnold, D., Arbel, T.: Mitigating calibration bias without fixed attribute grouping for improved fairness in medical imaging analysis. arXiv preprint arXiv:2307.01738 (2023)
  * [21] Singla, S., Eslami, M., Pollack, B., Wallace, S., Batmanghelich, K.: Explaining the black-box smoothly—a counterfactual approach. Medical Image Analysis 84, 102721 (2023)
  * [22] Targ, S., Almeida, D., Lyman, K.: Resnet in resnet: Generalizing residual architectures. arXiv preprint arXiv:1603.08029 (2016)
  * [23] Thiagarajan, J.J., Thopalli, K., Rajan, D., Turaga, P.: Training calibration-based counterfactual explainers for deep learning models in medical image analysis. Scientific reports 12(1), 597 (2022)
  * [24] Vapnik, V.: Principles of risk minimization for learning theory. Advances in neural information processing systems 4 (1991)
  * [25] Zhu, J.Y., Park, T., Isola, P., Efros, A.A.: Unpaired image-to-image translation using cycle-consistent adversarial networks. In: Proceedings of the IEEE international conference on computer vision. pp. 2223–2232 (2017)
  * [26] Zong, Y., Yang, Y., Hospedales, T.: MEDFAIR: Benchmarking fairness for medical imaging. In: International Conference on Learning Representations (2023)
  * [27] Zou, J., Schiebinger, L.: AI can be sexist and racist—it’s time to make it fair (2018)