implicit forward problem and rejecting if it fails; step [7.ii] corresponds to
solving the implicit backward problem and rejecting if it fails; and step
[7.iii] corresponds to the equality (18) and rejecting if it does not hold
numerically.

##  5 Numerical results

We illustrate in this section through two numerical examples the need to check
for $S$-reversibility to guarantee unbiased sampling. In particular, we show
that the Markov chains created via Algorithms 2 and 3 sample without bias the
target measure $\mu$ defined in (1) (and hence $\pi$ defined in (10) when
using the Hamiltonian function (9) as in RMHMC, see Section 3.1), while the
standard HMC and GHMC algorithms fail to do so for large time steps. We
consider sampling from a Boltzmann–Gibbs measure associated with a one-
dimensional double well potential in Section 5.1. We then sample from an
anisotropic two-dimensional potential in Section 5.2 and show the benefits of
using a position-dependent diffusion coefficient to faster reach stationarity.
The GSV numerical scheme (14) is used to integrate the Hamiltonian dynamics,
and consider the Hamiltonian function (9) and the momentum reversal map $S$
defined in (7). Note that, by Lemma 3, it suffices to check the
$S$-reversibility on the beginning/ending positions. Our implementation of the
GSV algorithm consists in using successively two Newton’s methods to first
fully solve the Euler A scheme, and then fully solve the Euler B scheme. This
corresponds to strategy (ii) in Remark 4. To check the invertibility of the
matrix $\nabla_{y}\Phi_{\Delta t}(x,y^{i})$ in step [6.i], we use the rank()
method from the Julia programming language, which amounts to counting the
number of singular values above a default threshold of
$m\varepsilon\sigma_{m}$ where $m$ is the size of the matrix, $\sigma_{m}$ its
largest singular value and $\varepsilon=2.22\times 10^{-16}$. The codes which
have been used to obtain the numerical results are available at the url
https://github.com/rsantet/RMHMC.

###  5.1 One dimensional confining double well potential

In this example, we show by computing histograms that not checking for
numerical reversibility induces a bias on the invariant probability measure
with respect to the target probability distribution, especially when time
steps become large; while using either Algorithms 2 or 3 with the numerical
flow with $S$-reversibility check (23) (implemented with Algorithm 7) allows
performing unbiased numerical sampling.

We consider a confining double well potential defined on $\mathbb{R}$ as

$V_{\sigma,h}(q)=q^{2}-1+\frac{h}{\sqrt{2\pi\sigma^{2}}}\exp\left(-\frac{q^{2}}{2\sigma^{2}}\right),$
(63)

and the target probability measure is the usual Gibbs measure $\pi$ defined by
(10). The parameters $\sigma$ and $h$ allows to vary the width and height of
the energy barrier at $q=0$. The potential (63) exhibits two wells as seen in
Figure 2, which enables us to evaluate the performance of our algorithm for
tackling multimodal target distributions. We choose $\sigma=0.2$, $h=1$, a
friction parameter $\gamma=1$ for the GHMC algorithm, and the diffusion
coefficient (which corresponds to the inverse mass in a Hamiltonian setting)

$D(q)=\left(\dfrac{1.5+\cos(\pi q)}{2}\right)^{2}\geqslant\dfrac{1}{16}>0.$

We run each simulation for $N_{\mathrm{iter}}=10^{7}$ time steps starting from
$q_{0}=-0.5$ and $p_{0}$ sampled from $\mathcal{N}(0,D(q_{0})^{-1})$. Newton’s
algorithm to solve the implicit steps has convergence tolerance thresholds of
$\eta_{\mathrm{Newton}}=\widetilde{\eta}_{\mathrm{Newton}}=10^{-12}$, and we
perform a maximum of $N_{\mathrm{Newton}}=100$ iterations for every Newton
algorithm loop. We use a tolerance threshold of $\eta_{\mathrm{rev}}=10^{-8}$
to check for $S$-reversibility. Observe that $\eta_{\mathrm{rev}}$ is somewhat
larger than $\eta_{\mathrm{Newton}},\widetilde{\eta}_{\mathrm{Newton}}$ to
take into account an error amplification phenomenon due to the iterations in
Newton’s method.

Figure 2: Confining double well potential with $\sigma=0.2$, $h=1$ and the
associated invariant measure proportional to $\mathrm{e}^{-V}$.

In our examples, not checking $S$-reversibility means that we use Algorithm 3
with $\psi_{\Delta t}^{\mathrm{rev}}$ replaced by

$\forall x\in\mathcal{X},\qquad\psi_{\Delta t}^{\mathrm{fwd}}(x)=\psi_{\Delta
t}(x)\mathbb{1}_{x\in\mathcal{A}_{\Delta
t}}+x\mathbb{1}_{x\notin\mathcal{A}_{\Delta t}}.$

In other words, only forward convergence is checked.

##### Observing the bias with histograms.

For $\Delta t$ small enough, checking for $S$-reversibility is not needed as
each implicit problem is well posed: the solution computed by the algorithm
satisfies $S$-reversibility, see Figure 3 (this is consistent with the results
from Section 4.1). When the time step increases, convergence problems for
Newton’s method appear and hence $S$-reversibility issues arise, see Figures 4
and 5. All these figures represent the stationary distributions of the Markov
chains: we have indeed checked that the obtained histograms do not vary for a
larger number of iterations. The fact that the histograms are even (as
expected from the symmetry of the potential) is also a good indication that
the stationary state has been reached since the initial condition is not
symmetric. These results illustrate the fact that $S$-reversibility checks are
needed to create a Markov chain which yields an unbiased sampling of the Gibbs
measure $\pi$.

Figure 3: Sampling with a time step $\Delta t=0.15$. The histogram on the left
is obtained thanks to Algorithm 3, while the one on the right is obtained
without checking $S$-reversibility. There is no significant difference between
the two choices. Figure 4: Same setting as Figure 3, with a time step $\Delta
t=0.69$. Significant bias can be observed when $S$-reversibility is not
checked. Figure 5: Same setting as Figure 3, with a time step $\Delta t=1.08$.
Significant bias can be observed when $S$-reversibility is not checked.

##### Rejection probabilities and weak order consistency.

We present in Figure 6 the rejection probabilities as a function of the time
step $\Delta t$. We also give the values for the time steps associated with
the histograms of Figures 3 to 5 in Table 1. For small time steps, there is no
rejection due to $S$-reversibility checks, and the global rejection
probability is only composed of the rejection due to the Metropolis–Hastings
acceptance/rejection step. The global rejection probability therefore scales
as $\mathrm{O}(\Delta t^{3})$ for $\Delta t$ small, and we observe numerically
that the rejection probability due to the $S$-reversibility is very small and
that $r_{3}\geqslant 3$ in any case. This is consistent with (42), and suggest
that the HMC algorithm (Algorithm 2) using the Hamiltonian function (9),
_i.e._ the RMHMC algorithm, yields a weakly consistent discretization of order
1 of the overdamped Langevin dynamics with a position-dependent diffusion
coefficient (33), following the result of Proposition 11.

For large time steps, the rejection due to $S$-reversibility checks become
prominent. In particular, forward convergence and the $S$-reversibility
condition (18) are difficult to satisfy. It is interesting to note that, when
there is convergence of the implicit forward problem, the implicit backward
problem very often converges as well, but not necessarily to the solution
compatible with $S$-reversibility. This highlights the need to implement
$S$-reversibility checks to perform unbiased sampling.

Time step | Forward | Backward |  $S$-reversibility | Metropolis–Hastings | Global
---|---|---|---|---|---
$0.15$ | $4.8\times 10^{-1}$ | $5.1\times 10^{-4}$ | $1.3\times 10^{-3}$ | $2.6$ | $3.1$
$0.69$ | $27$ | $0.5$ | $23.9$ | $13$ | $64$
$1.08$ | $34$ | $1.2$ | $44$ | $6.8$ | $86$
Table 1: Decomposition in percentages of the rejection probabilities for
various time steps $\Delta t$. “Forward” refers to non-convergent iterations
to solve the implicit forward problem; “Backward” refers to a convergence of
the implicit forward problem but a non-convergent implicit backward problem;
“S-reversibility” refers to convergent implicit forward and backward problems
for which (18) is not satisfied; upon acceptance in the three previous steps,
“Metropolis–Hastings” refers to a rejection in the acceptance/rejection step
[3.iii]. Finally, “Global” is the global rejection probability, namely the sum
of all the previous columns. Figure 6: Rejection probabilities as a function
of the time step $\Delta t\in[10^{-2},10]$. We refer to the caption of Table 1
for the meanings of the various labels.

###  5.2 Anisotropic two-dimensional potential

Figure 7: Invariant probability measure with density proportional to
$\mathrm{e}^{-V}$ for the anisotropic two-dimensional potential
$V(x,y)=100(x^{2}+y^{2}-1)^{2}$.

We demonstrate on a representative example that using a non-constant diffusion
coefficient helps to reach equilibrium faster. In particular, using the
Hamiltonian function defined in (9), we compare the use of a constant
diffusion (which corresponds to the usual HMC scheme) to a diffusion
coefficient which takes into account the anisotropy of the potential energy
landscape (using the RMHMC algorithm). In the first case, GSV is an explicit
scheme as the Hamiltonian function is separable, so that standard HMC results
apply and unbiased numerical sampling is performed. In the second case,
implicit schemes, hence $S$-reversibility checks, have to be used in order to
guarantee unbiasedness.

We consider a target probability measure $\pi$ with density proportional to
$\mathrm{e}^{-V(q)}$ with $q=(x,y)\in\mathbb{R}^{2}$ and
$V(x,y)=100(x^{2}+y^{2}-1)^{2}$, whose support is essentially an annulus in
$\mathbb{R}^{2}$, see Figure 7. The width of the annulus is very small. This
suggests that, when using an isotropic diffusion coefficient
$D(x,y)=\mathrm{I}_{2}$, the Metropolis–Hastings procedure will likely reject
the proposal as it typically ends up in high energy regions. Instead, one can
choose a diffusion coefficient which takes into account the geometry of the
potential energy landscape so that the proposal is more likely to be close to
the annulus. To that aim, we propose the following diffusion coefficient:

$D^{\mathrm{aniso}}(x,y)=\varepsilon\mathrm{I}_{2}+tt^{\mathsf{T}},$ (64)

where

$t=\frac{1}{\sqrt{x^{2}+y^{2}}}\begin{pmatrix}-y\\\ x\end{pmatrix},$

and $\varepsilon>0$ is a small parameter ensuring that the diffusion
coefficient is always positive definite. Looking at the overdamped Langevin
dynamics (33), the diffusion coefficient defined in (64) helps to project the
noise onto the tangent plane to the circle of center 0 and radius 1. In the
extended phase-space setting of (G)HMC, this is equivalent to drawing momenta
which are essentially tangent to the circle.

We compare the speeds of convergence towards equilibrium of two Markov chains:
one with the diffusion coefficient as defined in (64), and one with the
diffusion coefficient

$D^{\mathrm{iso}}(q)=(1+\varepsilon)\mathrm{I}_{2}.$ (65)

Both diffusion coefficients are normalized in order to have the same spectral
radius $1+\varepsilon$. To analyze the results, we consider the polar
coordinates $(x,y)=(r\cos\theta,r\sin\theta)$. As the potential exhibits a
spherical symmetry, the random variable $\theta$ should be uniformly
distributed on $[0,2\pi)$. In order to measure convergence, we thus compute
the total variation distance between the empirical distribution on $\theta$
and the uniform distribution on $[0,2\pi)$ after a fixed number of iterations
of the dynamics (in order to fix the computational cost), and average over a
given number of simulations, each time starting the dynamics randomly on the
circle of radius 1. A small total variation indicates that the dynamics
explores faster the annulus. In practice, the total variation is approximated
by discretizing the interval $[0,2\pi)$ into $N_{\mathrm{bins}}$ bins of equal
sizes, and comparing the probabilities of each bin. More precisely, for
$1\leqslant i\leqslant N_{\mathrm{bins}}$, we denote by
$B_{i}=[(i-1)/N_{\mathrm{bins}},i/N_{\mathrm{bins}})$ the $i$-th bin. The
total variation distance between two probability measures $\mu,\nu$, namely
$\mathrm{TV}(\mu,\nu)=\frac{1}{2}\int\mathrm{d}\left\lvert\mu-\nu\right\rvert$,
is numerically approximated as

$\mathrm{TV}_{N_{\mathrm{bins}}}(\mu,\nu)=\frac{1}{2}\frac{2\pi}{N_{\mathrm{bins}}}\sum_{i=1}^{N_{\mathrm{bins}}}\left\lvert\widehat{\mu}_{i}-\widehat{\nu}_{i}\right\rvert,$

where $\widehat{\mu}_{i},\widehat{\nu}_{i}\geqslant 0$ are proportional to the
fraction of points in the $i$-th bin, and normalized as

$\frac{2\pi}{N_{\mathrm{bins}}}\sum_{i=1}^{N_{\mathrm{bins}}}\widehat{\mu}_{i}=\frac{2\pi}{N_{\mathrm{bins}}}\sum_{i=1}^{N_{\mathrm{bins}}}\widehat{\nu}_{i}=1.$

In particular, for the uniform law on $[0,\pi)$, it holds
$\widehat{\mu}_{1}=\dots=\widehat{\mu}_{N_{\mathrm{bins}}}=1/(2\pi)$. Note
that the total variation distance is bounded by 1 with our convention.

We run the simulations for 32 values of $\Delta t$, evenly spread log-wise on
$[10^{-3},2]$, the number of iterations being fixed to
$N_{\mathrm{iter}}=10^{6}$ for each value of $\Delta t$. Each simulation is
performed $N_{\mathrm{real}}=10^{3}$ times. The approximation of the total
variation distances relies on $N_{\mathrm{bins}}=100$. The other parameters
are the same as in Section 5.1. The results are presented in Figure 8, where
we plot the mean total variation distance for each time step. The standard
deviation, quantified using the variance estimated over the independent
$N_{\mathrm{real}}$ realizations, is not shown in the plot as the error bars
are smaller than the symbol sizes.

Figure 8: Comparison between the isotropic and anisotropic diffusion
coefficients in log-log scale, using either Algorithms 2 (HMC) or 3 (GHMC).
“Aniso” corresponds to the diffusion coefficient defined by (64), “Iso” to the
one defined by (65). The red circles correspond to the minima of each curve,
_i.e._ the point (optimal time step, minimum average TV distance). Minimum
value/$\Delta t$ for Iso_HMC: $9.17\times 10^{-2}$, $\Delta t=6.46\times
10^{-2}$, Aniso_HMC: $3.58\times 10^{-2}$, $\Delta t=1.72\times 10^{-1}$, for
Iso_GHMC: $4.78\times 10^{-2}$, $\Delta t=3.96\times 10^{-2}$, for Aniso_GHMC:
$2.37\times 10^{-2}$, $\Delta t=1.05\times 10^{-1}$.

The results show that the optimal time step, _i.e._ the one which minimizes in
average the total variation distance, yields a lower error of the empirical
angle distribution with respect to the uniform distribution on $[0,2\pi)$ when
the anisotropic diffusion coefficient (64) is chosen. The error made is 2
times smaller in the HMC case and 2.5 times smaller in the GHMC case.
Moreover, for a given algorithm, the optimal time step using the anisotropic
diffusion coefficient is larger than the one with the isotropic diffusion
coefficient. This is in agreement with the fact that larger moves are being
accepted as the noise is projected back onto the circle.

Note also that the optimal minimum error corresponding to the GHMC method is
lower than the one corresponding to HMC. Hence, the additional computation
cost arising from integrating the Ornstein–Uhlenbeck process (namely steps
[3.i] and [3.v]) is compensated by a faster exploration of the configuration
space.

For small time steps, the total variation distance scales as $\Delta t^{-1}$
for the HMC algorithm and as $\Delta t^{-1/2}$ for the GHMC algorithm. This is
consistent with the facts that:

  1. (i)

at the continuous level, the error is of order $1/\sqrt{t_{\mathrm{sim}}}$
where $t_{\mathrm{sim}}$ is the simulation time (by a standard central limit
theorem for continuous time Markov processes);

  2. (ii)

GHMC is formally a weakly consistent discretization of the Langevin dynamics
(27) (see Section 3.2.2), so that in this case
$t_{\mathrm{sim}}=N_{\mathrm{iter}}\Delta t$;

  3. (iii)

HMC is a weakly consistent discretization of order 1 of the overdamped
Langevin dynamics with a position-dependent diffusion coefficient (33) and
effective time step $h=\Delta t^{2}/2$ (see Proposition 11), so that in this
case $t_{\mathrm{sim}}=N_{\mathrm{iter}}h$.

In this small time step regime, we also observe that there is no difference
between using the isotropic or the anisotropic diffusion coefficient: the
simulation time is not long enough to see the benefits from using an
anisotropic diffusion coefficient. In fact, we see that the total variation is
very close to 1 for the HMC algorithm, hence a saturation effect emerges: the
average error simply converges to the total variation distance between the
initial condition and the target probability measure. This saturation is less
pronounced in GHMC thanks to the momentum variable, which allows for an
initial ballistic exploration.

For large time steps, we observe that there is no difference between HMC and
GHMC. Indeed, the proposal obtained by integrating (or attempting to
integrate) the Hamiltonian dynamics is almost always rejected. Eventually, the
algorithm only resamples the momenta, either by direct sampling in the HMC
case (see step [2.i]) or by partial resampling in the GHMC case (see steps
[3.i]–[3.v]).

##  6 Proof of various technical results

We provide in this section the proof of various technical results stated in
Sections 2.2, 2.4, 3.3, 4.1 and 4.2.

###  6.1 Proof of Proposition 2

##### The IMR numerical scheme.

Using the definition of the IMR numerical scheme (13), we readily compute

$\forall(x,x_{1})\in\mathcal{X}^{2},\qquad\nabla_{x}\Phi_{\Delta
t}^{\mathrm{IMR}}(x,x_{1})=\mathrm{I}_{d}-A(x,x_{1}),\qquad\nabla_{y}\Phi_{\Delta
t}^{\mathrm{IMR}}(x,x_{1})=\mathrm{I}_{d}+A(x,x_{1}),$

where $A(x,x_{1})=\frac{\Delta t}{2}J\nabla^{2}H(x,x_{1})$ with $J$ the matrix
defined in (2) and

$\nabla^{2}H(x,x_{1})=\begin{pmatrix}\nabla_{q}^{2}H\left(\frac{x+x_{1}}{2}\right)&\nabla_{p,q}^{2}H\left(\frac{x+x_{1}}{2}\right)\\\\[5.69046pt]
\nabla_{q,p}^{2}H\left(\frac{x+x_{1}}{2}\right)&\nabla_{p}^{2}H\left(\frac{x+x_{1}}{2}\right)\end{pmatrix}.$

Note that $\nabla^{2}H(x,x_{1})$ is symmetric since $H$ is $\mathcal{C}^{2}$.
Let $x\in\mathcal{A}_{\Delta t}$. Since $\nabla_{y}\Phi_{\Delta
t}^{\mathrm{IMR}}\left(x,\varphi_{\Delta t}^{\mathrm{IMR}}(x)\right)$ is
invertible, the matrix $\mathrm{I}_{d}+A\left(x,\varphi_{\Delta
t}^{\mathrm{IMR}}(x)\right)$ is invertible. In the remainder of the proof, we
drop the arguments for readability: the argument of the matrix $A$ is fixed to
$(x,\varphi_{\Delta t}^{\mathrm{IMR}}(x))$. Using the fact that
$J^{2}=-\mathrm{I}_{d}$, $J^{-1}=J^{\mathsf{T}}=-J$, $JA^{\mathsf{T}}=-AJ$ and
$A^{\mathsf{T}}J=-JA$, it holds
$J(\mathrm{I}_{d}-A)=(\mathrm{I}_{d}+A^{\mathsf{T}})J$. Thus, the matrix
$\mathrm{I}_{d}-A^{\mathsf{T}}$ is invertible and

$(\mathrm{I}_{d}-A)^{-\mathsf{T}}J=J(\mathrm{I}_{d}+A)^{-1}.$

By Lemma 1, it therefore holds

$\displaystyle\nabla\varphi_{\Delta}(x)^{\mathsf{T}}J\nabla\varphi_{\Delta
t}(x)$
$\displaystyle=(\mathrm{I}_{d}+A)^{\mathsf{T}}(\mathrm{I}_{d}-A)^{-\mathsf{T}}J(\mathrm{I}_{d}-A)^{-1}(\mathrm{I}_{d}+A)$
$\displaystyle=(\mathrm{I}_{d}+A^{\mathsf{T}})J(\mathrm{I}_{d}+A)^{-1}(\mathrm{I}_{d}-A)^{-1}(\mathrm{I}_{d}+A)$
$\displaystyle=J(\mathrm{I}_{d}-A)(\mathrm{I}_{d}+A)^{-1}(\mathrm{I}_{d}-A)^{-1}(\mathrm{I}_{d}+A)=J,$

where we used the fact that the matrices $(\mathrm{I}_{d}-A)^{-1}$ and
$(\mathrm{I}_{d}+A)^{-1}$ commute for the last equality. This shows that the
numerical flow $\varphi_{\Delta t}^{\mathrm{IMR}}=\chi_{\Delta
t}^{\mathrm{IMR}}$ is symplectic on its domain.

##### The GSV numerical scheme.

Since the composition of symplectic methods is symplectic, and the GSV
numerical scheme can be seen as the composition of the Euler B and Euler A
symplectic methods, it should not be surprising that it is also symplectic.
However, the proof still requires some care as we defined numerical flows on
subsets of the whole configuration space and constructed the GSV numerical
scheme by stacking the two Euler methods into a single vector (and not
composing them). Therefore, one has to check that a numerical solver for the
GSV numerical scheme allows us to define numerical flows for the Euler B and
Euler A schemes, and that we still recover a symplectic numerical flow with
our construction.

Using the definition of the GSV numerical scheme (15), we readily compute
(recall that differentiating with respect to $y$ corresponds to computing the
gradient with respect to the last $k$ coordinates):

$\forall(x,x_{1},x_{2})\in\mathcal{X}^{3},\qquad\left\\{\begin{aligned}
&\nabla_{x}\Phi_{\Delta
t}^{\mathrm{GSV}}(x,x_{1},x_{2})=\begin{pmatrix}\nabla_{x}\Phi_{\Delta
t/2}^{\mathrm{Euler\,B}}(x,x_{1})\\\ 0_{d}\end{pmatrix}\in\mathbb{R}^{2d\times
d},\\\ &\nabla_{y}\Phi_{\Delta
t}^{\mathrm{GSV}}(x,x_{1},x_{2})=\begin{pmatrix}\nabla_{y}\Phi_{\Delta
t/2}^{\mathrm{Euler\,B}}(x,x_{1})&0_{d}\\\ \nabla_{x}\Phi_{\Delta
t/2}^{\mathrm{Euler\,A}}(x_{1},x_{2})&\nabla_{y}\Phi_{\Delta
t/2}^{\mathrm{Euler\,A}}(x_{1},x_{2})\end{pmatrix}\in\mathbb{R}^{2d\times
2d}.\end{aligned}\right.$ (66)

Let $\chi_{\Delta t}^{\mathrm{GSV}}\colon\mathcal{A}_{\Delta
t}^{\mathrm{GSV}}\to\mathcal{X}^{2}$ be a numerical solver for the GSV
numerical scheme. For $x\in\mathcal{A}_{\Delta t}^{\mathrm{GSV}}$, if
$\chi_{\Delta t}(x)=(x_{1},x_{2})$, we define $\varphi_{\Delta
t/2}^{\mathrm{Euler\,B}}(x)=x_{1}$ and $\varphi_{\Delta
t/2}^{\mathrm{Euler\,A}}(x_{1})=x_{2}$. The map $\varphi_{\Delta
t/2}^{\mathrm{Euler\,B}}\colon\mathcal{A}_{\Delta
t}^{\mathrm{GSV}}\to\mathcal{X}$ (respectively the map $\varphi_{\Delta
t/2}^{\mathrm{Euler\,A}}\colon\varphi_{\Delta
t/2}^{\mathrm{Euler\,B}}\left(\mathcal{A}_{\Delta
t}^{\mathrm{GSV}}\right)\to\mathcal{X}$) defines a numerical flow for the
numerical scheme $\Phi_{\Delta t/2}^{\mathrm{Euler\,B}}$ (respectively for the
numerical scheme $\Phi_{\Delta t/2}^{\mathrm{Euler\,A}}$). Indeed, if
$x\in\mathcal{A}_{\Delta t}^{\mathrm{GSV}}$, then it holds

$\Phi_{\Delta t}^{\mathrm{GSV}}\left(x,\varphi_{\Delta
t/2}^{\mathrm{Euler\,B}}(x),\varphi_{\Delta
t/2}^{\mathrm{Euler\,A}}\circ\varphi_{\Delta
t/2}^{\mathrm{Euler\,B}}(x)\right)=0,$

which implies that $\Phi_{\Delta
t/2}^{\mathrm{Euler\,B}}\left(x,\varphi_{\Delta
t/2}^{\mathrm{Euler\,B}}(x)\right)=0$ and $\Phi_{\Delta
t/2}^{\mathrm{Euler\,A}}\left(\varphi_{\Delta
t/2}^{\mathrm{Euler\,B}}(x),\varphi_{\Delta
t/2}^{\mathrm{Euler\,A}}\circ\varphi_{\Delta
t/2}^{\mathrm{Euler\,B}}(x)\right)=0$. Moreover, the second equality in (66)
shows that the matrix $\nabla_{y}\Phi_{\Delta
t}^{\mathrm{GSV}}\left(x,\chi_{\Delta t}^{\mathrm{GSV}}(x)\right)$ is
invertible if and only if the matrices $\nabla_{y}\Phi_{\Delta
t/2}^{\mathrm{Euler\,B}}\left(x,\varphi_{\Delta
t/2}^{\mathrm{Euler\,B}}(x)\right)$ and $\nabla_{y}\Phi_{\Delta
t/2}^{\mathrm{Euler\,A}}\left(\varphi_{\Delta
t/2}^{\mathrm{Euler\,B}}(x),\varphi_{\Delta
t/2}^{\mathrm{Euler\,A}}\circ\varphi_{\Delta
t/2}^{\mathrm{Euler\,B}}(x)\right)$ are invertible. Finally, notice that
$\varphi_{\Delta t/2}^{\mathrm{Euler\,B}}\left(\mathcal{A}_{\Delta
t}^{\mathrm{GSV}}\right)$ is open and nonempty as a simple consequence of the
implicit function theorem.

We now compute the gradient of $\chi_{\Delta t}^{\mathrm{GSV}}$ using Lemma 1.
We drop the use of arguments of functions for readability: we fix
$x\in\mathcal{A}_{\Delta t}^{\mathrm{GSV}}$, so that $\left(x,\varphi_{\Delta
t/2}^{\mathrm{Euler\,B}}(x),\varphi_{\Delta
t/2}^{\mathrm{Euler\,A}}(x)\right)$ is the argument of $\Phi_{\Delta
t}^{\mathrm{GSV}}$, $\left(x,\varphi_{\Delta
t/2}^{\mathrm{Euler\,B}}(x)\right)$ of $\Phi_{\Delta
t/2}^{\mathrm{Euler\,B}}$, $\left(\varphi_{\Delta
t/2}^{\mathrm{Euler\,B}}(x),\varphi_{\Delta
t/2}^{\mathrm{Euler\,A}}\circ\varphi_{\Delta
t/2}^{\mathrm{Euler\,B}}(x)\right)$ of $\Phi_{\Delta
t/2}^{\mathrm{Euler\,A}}$, and $x$ of $\chi_{\Delta t}^{\mathrm{GSV}}$,
$\varphi_{\Delta t/2}^{\mathrm{Euler\,B}}$ and $\varphi_{\Delta
t/2}^{\mathrm{Euler\,A}}$.

First, notice that with straightforward computations, it holds

$\left(\nabla_{y}\Phi_{\Delta
t}^{\mathrm{GSV}}\right)^{-1}=\begin{pmatrix}\left(\nabla_{y}\Phi_{\Delta
t/2}^{\mathrm{Euler\,B}}\right)^{-1}&0_{d}\\\
0_{d}&\left(\nabla_{y}\Phi_{\Delta
t/2}^{\mathrm{Euler\,A}}\right)^{-1}\end{pmatrix}\begin{pmatrix}\mathrm{I}_{d}&0_{d}\\\
-\nabla_{x}\Phi_{\Delta t/2}^{\mathrm{Euler\,A}}\left(\nabla_{y}\Phi_{\Delta
t/2}^{\mathrm{Euler\,B}}\right)^{-1}&\mathrm{I}_{d}\end{pmatrix}.$ (67)

Then multiplying (67) with the first equation in (66), one obtains

$\displaystyle\nabla\chi_{\Delta t}^{\mathrm{GSV}}$
$\displaystyle=-\begin{pmatrix}\left(\nabla_{y}\Phi_{\Delta
t/2}^{\mathrm{Euler\,B}}\right)^{-1}&0_{d}\\\
0_{d}&\left(\nabla_{y}\Phi_{\Delta
t/2}^{\mathrm{Euler\,A}}\right)^{-1}\end{pmatrix}\begin{pmatrix}\nabla_{x}\Phi_{\Delta
t/2}^{\mathrm{Euler\,B}}\\\ -\nabla_{x}\Phi_{\Delta
t/2}^{\mathrm{Euler\,A}}\left(\nabla_{y}\Phi_{\Delta
t/2}^{\mathrm{Euler\,B}}\right)^{-1}\nabla_{x}\Phi_{\Delta
t/2}^{\mathrm{Euler\,B}}\end{pmatrix}$
$\displaystyle=-\begin{pmatrix}\left(\nabla_{y}\Phi_{\Delta
t/2}^{\mathrm{Euler\,B}}\right)^{-1}\nabla_{x}\Phi_{\Delta
t/2}^{\mathrm{Euler\,B}}\\\ -\left(\nabla_{y}\Phi_{\Delta
t/2}^{\mathrm{Euler\,A}}\right)^{-1}\nabla_{x}\Phi_{\Delta
t/2}^{\mathrm{Euler\,A}}\left(\nabla_{y}\Phi_{\Delta
t/2}^{\mathrm{Euler\,B}}\right)^{-1}\nabla_{x}\Phi_{\Delta
t/2}^{\mathrm{Euler\,B}}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\nabla\varphi_{\Delta
t/2}^{\mathrm{Euler\,B}}\\\ \nabla\varphi_{\Delta
t/2}^{\mathrm{Euler\,A}}\nabla\varphi_{\Delta
t/2}^{\mathrm{Euler\,B}}\end{pmatrix},$

where we used Lemma 1 for the last equality. In particular,

$\nabla\varphi_{\Delta t}^{\mathrm{GSV}}=\nabla\varphi_{\Delta
t/2}^{\mathrm{Euler\,A}}\nabla\varphi_{\Delta t/2}^{\mathrm{Euler\,B}}.$

Since the two numerical flows $\varphi_{\Delta t/2}^{\mathrm{Euler\,B}}$ and
$\varphi_{\Delta t/2}^{\mathrm{Euler\,A}}$ are symplectic on their domains
(using similar arguments as above and [16, Theorem VI.3.3]), it holds

$\displaystyle\left(\nabla\varphi_{\Delta
t}^{\mathrm{GSV}}\right)^{\mathsf{T}}J\nabla\varphi_{\Delta t}$
$\displaystyle=\left(\nabla\varphi_{\Delta
t/2}^{\mathrm{Euler\,B}}\right)^{\mathsf{T}}\left(\nabla\varphi_{\Delta
t/2}^{\mathrm{Euler\,A}}\right)^{\mathsf{T}}J\nabla\varphi_{\Delta
t/2}^{\mathrm{Euler\,A}}\nabla\varphi_{\Delta t/2}^{\mathrm{Euler\,B}}$
$\displaystyle=\left(\nabla\varphi_{\Delta
t/2}^{\mathrm{Euler\,B}}\right)^{\mathsf{T}}J\nabla\varphi_{\Delta
t/2}^{\mathrm{Euler\,B}}=J.$

This shows that the numerical flow $\varphi_{\Delta t}^{\mathrm{GSV}}$ is
symplectic on its domain.

###  6.2 Proof of Lemma 6

This proof is an adaptation of the proof of [26, Lemma 2.3]. As both
$\varphi_{\Delta t}$ and $S$ are continuous and $\mathcal{A}_{\Delta t}$ is
open, $\mathcal{A}_{\Delta t}\bigcap\psi_{\Delta t}^{-1}(\mathcal{A}_{\Delta
t})$ is an open set of $\mathcal{X}$. Let us prove (25) by contradiction.
Consider a non-empty path connected component $C$ of $\mathcal{A}_{\Delta
t}\bigcap\psi_{\Delta t}^{-1}(\mathcal{A}_{\Delta t})$, an element $X_{0}\in
C$ such that

$\chi_{\Delta t}\circ\psi_{\Delta t}(X_{0})=\Bigl{(}S\circ\chi_{\Delta
t,k-1}(X_{0}),\dots,S\circ\chi_{\Delta t,1}(X_{0}),S(X_{0})\Bigr{)}.$

and assume that there exists $X_{1}\in C$ such that

$\chi_{\Delta t}\circ\psi_{\Delta t}(X_{1})\neq\Bigl{(}S\circ\chi_{\Delta
t,k-1}(X_{1}),\dots,S\circ\chi_{\Delta t,1}(X_{1}),S(X_{1})\Bigr{)}.$ (68)

Let us connect the elements $X_{0}$ and $X_{1}$ in $C$ by a continuous path
$[0,1]\ni\theta\mapsto X_{\theta}$. Define

$\theta_{\star}=\sup\left\\{\theta\geqslant 0\,\middle|\,\begin{aligned}
&\forall s\in[0,\theta],\\\ &\chi_{\Delta t}\circ\psi_{\Delta
t}(X_{s})=\Bigl{(}S\circ\chi_{\Delta t,k-1}(X_{s}),\dots,S\circ\chi_{\Delta
t,1}(X_{s}),S(X_{s})\Bigr{)}\end{aligned}\right\\}.$

By continuity of $\chi_{\Delta t}$, $S$ and $\theta\mapsto X_{\theta}$, the
value $\theta_{\star}$ is actually a maximum, and it holds $\theta_{\star}<1$
by (68). By definition of the supremum, there exists a sequence
$(\theta_{n})_{n\geqslant 0}$ with values in $(\theta_{\star},1)$ such that
$\theta_{n}\to\theta_{\star}$ as $n\to+\infty$ and

$\forall n\geqslant 0,\qquad\chi_{\Delta t}\circ\psi_{\Delta
t}(X_{\theta_{n}})\neq\Bigl{(}S\circ\chi_{\Delta
t,k-1}(X_{\theta_{n}}),\dots,S\circ\chi_{\Delta
t,1}(X_{\theta_{n}}),S(X_{\theta_{n}})\Bigr{)}.$

We can therefore define the $dk$-dimensional unit vectors

$\tau_{n}=\frac{\chi_{\Delta t}\circ\psi_{\Delta
t}(X_{\theta_{n}})-\Bigl{(}S\circ\chi_{\Delta
t,k-1}(X_{\theta_{n}}),\dots,S\circ\chi_{\Delta
t,1}(X_{\theta_{n}}),S(X_{\theta_{n}})\Bigr{)}}{\left\lVert\chi_{\Delta
t}\circ\psi_{\Delta t}(X_{\theta_{n}})-\Bigl{(}S\circ\chi_{\Delta
t,k-1}(X_{\theta_{n}}),\dots,S\circ\chi_{\Delta
t,1}(X_{\theta_{n}}),S(X_{\theta_{n}})\Bigr{)}\right\rVert},$ (69)

since the denominator is positive. By compactness of the $dk$-dimensional unit
ball and upon extraction (keeping the same notation for the extracted
subsequence), there exists a subsequence $(\theta_{n})_{n\geqslant 0}$ and a
unit vector $\tau_{\star}$ such that $\tau_{n}\to\tau_{\star}$ as
$n\to+\infty$. In particular, it holds
$\left\lVert\tau_{\star}\right\rVert=1$. We can then rewrite (69) as

$\chi_{\Delta t}\circ\psi_{\Delta
t}(X_{\theta_{n}})=\Bigl{(}S\circ\chi_{\Delta
t,k-1}(X_{\theta_{n}}),\dots,S\circ\chi_{\Delta
t,1}(X_{\theta_{n}}),S(X_{\theta_{n}})\Bigr{)}+\varepsilon_{n}\tau_{n},$

with $\left\lVert\tau_{n}\right\rVert=1$ and

$\varepsilon_{n}=\left\lVert\chi_{\Delta t}\circ\psi_{\Delta
t}(X_{\theta_{n}})-\Bigl{(}S\circ\chi_{\Delta
t,k-1}(X_{\theta_{n}}),\dots,S\circ\chi_{\Delta
t,1}(X_{\theta_{n}}),S(X_{\theta_{n}})\Bigr{)}\right\rVert\xrightarrow[n\to+\infty]{}0$

by continuity of the mapping $\theta\mapsto X_{\theta}$ and the fact that the
latter quantity vanishes at $\theta_{\star}$. Using the fact that
$\psi_{\Delta t}(X_{\theta_{n}})\in\mathcal{A}_{\Delta t}$ and (16), and then
performing a first order Taylor expansion on $\Phi_{\Delta t}$,

$\displaystyle 0$ $\displaystyle=\Phi_{\Delta t}\left(\psi_{\Delta
t}(X_{\theta_{n}}),\chi_{\Delta t}\circ\psi_{\Delta
t}(X_{\theta_{n}})\right),$ (70) $\displaystyle=\Phi_{\Delta
t}\left(\psi_{\Delta t}(X_{\theta_{n}}),\Bigl{(}S\circ\chi_{\Delta
t,k-1}(X_{\theta_{n}}),\dots,S\circ\chi_{\Delta
t,1}(X_{\theta_{n}}),S(X_{\theta_{n}})\Bigr{)}\right)$
$\displaystyle\quad+\varepsilon_{n}\nabla_{y}\Phi_{\Delta t}\left(\psi_{\Delta
t}(X_{\theta_{n}}),\Bigl{(}S\circ\chi_{\Delta
t,k-1}(X_{\theta_{n}}),\dots,S\circ\chi_{\Delta
t,1}(X_{\theta_{n}}),S(X_{\theta_{n}})\Bigr{)}\right)\tau_{n}$
$\displaystyle\quad+\mathrm{O}\left(\varepsilon_{n}\right).$

The first term on the right-hand side of (70) vanishes. Indeed, from the
$S$-reversibility of $\Phi_{\Delta t}$ and (16),

$\displaystyle\Phi_{\Delta t}\left(\psi_{\Delta
t}(X_{\theta_{n}}),\Bigl{(}S\circ\chi_{\Delta
t,k-1}(X_{\theta_{n}}),\dots,S\circ\chi_{\Delta
t,1}(X_{\theta_{n}}),S(X_{\theta_{n}})\Bigr{)}\right)$
$\displaystyle=\Phi_{\Delta t}\Bigl{(}X_{\theta_{n}},S\circ S\circ\chi_{\Delta
t,1}(X_{\theta_{n}}),\dots,S\circ S\circ\chi_{\Delta
t,k-1}(X_{\theta_{n}}),S\circ S\circ\varphi_{\Delta
t}(X_{\theta_{n}})\Bigr{)}$ $\displaystyle=\Phi_{\Delta
t}\Bigl{(}X_{\theta_{n}},\chi_{\Delta t}(X_{\theta_{n}})\Bigr{)}=0.$

Dividing (70) by $\varepsilon_{n}\neq 0$, and taking the limit as
$n\to+\infty$, we obtain

$\nabla_{y}\Phi_{\Delta t}\Bigl{(}\psi_{\Delta
t}(X_{\theta_{\star}}),S\circ\chi_{\Delta
t,k-1}(X_{\theta_{\star}}),\dots,S\circ\chi_{\Delta
t,1}(X_{\theta_{\star}}),S(X_{\theta_{\star}})\Bigr{)}\tau_{\star}=0.$ (71)

Now, the matrix $\nabla_{y}\Phi_{\Delta t}(\psi_{\Delta
t}(X_{\theta_{\star}}),S\circ\chi_{\Delta
t,k-1}(X_{\theta_{\star}}),\dots,S\circ\chi_{\Delta
t,1}(X_{\theta_{\star}}),S(X_{\theta_{\star}}))$ is invertible as, from the
definition of $\theta_{\star}$,

$\displaystyle\nabla_{y}\Phi_{\Delta t}(\psi_{\Delta
t}(X_{\theta_{\star}}),S\circ\chi_{\Delta
t,k-1}(X_{\theta_{\star}}),\dots,S\circ\chi_{\Delta
t,1}(X_{\theta_{\star}}),S(X_{\theta_{\star}}))\phantom{\tau_{\star}=0.}$
$\displaystyle\qquad=\left(\nabla_{y}\Phi_{\Delta t}\right)\left(\psi_{\Delta
t}(X_{\theta_{\star}}),\chi_{\Delta t}\circ\psi_{\Delta
t}(X_{\theta_{\star}})\right),$

the latter matrix being invertible by the fact that $\psi_{\Delta
t}(X_{\theta_{\star}})\in\mathcal{A}_{\Delta t}$ and the definition of
$\mathcal{A}_{\Delta t}$. As $\tau_{\star}$ is a nonzero vector
($\left\lVert\tau_{\star}\right\rVert=1$), Equation (71) leads to a
contradiction. This shows that there is no element $X_{1}\in C$ such that

$\chi_{\Delta t}\circ\psi_{\Delta t}(X_{1})\neq\Bigl{(}S\circ\chi_{\Delta
t,k-1}(X_{1}),\dots,S\circ\chi_{\Delta t,1}(X_{1}),S(X_{1})\Bigr{)},$

which concludes the proof of (25).

###  6.3 Proof of Lemma 10

We use the following notation: for $q\in\mathcal{O}$ and for any
$\xi_{1},\xi_{2}\in\mathbb{R}^{m}$,

$\displaystyle\nabla D(q)=\left(\frac{\partial D}{\partial
q_{i}}(q)\right)_{1\leqslant i\leqslant m}\in\mathbb{R}^{m^{2}\times m},$
$\displaystyle\xi_{1}^{\mathsf{T}}\nabla
D(q)\xi_{2}=\left(\xi_{1}^{\mathsf{T}}\frac{\partial D}{\partial
q_{i}}(q)\xi_{2}\right)_{1\leqslant i\leqslant m}\in\mathbb{R}^{m},$
$\displaystyle\Tr\left(D(q)^{-1}\nabla
D(q)\right)=\left(\Tr\left(D(q)^{-1}\frac{\partial D}{\partial
q_{i}}(q)\right)\right)_{1\leqslant i\leqslant m}\in\mathbb{R}^{m},$
$\displaystyle F(q)=-\nabla V(q)+\frac{1}{2}\Tr\left(D(q)^{-1}\nabla
D(q)\right)\in\mathbb{R}^{m}.$

When we write a remainder term as $c(q,p)=\mathrm{O}(\Delta t^{\delta})$
(where $c(q,p)$ is either a vector or a matrix), it means that there exist
$\Delta t_{\star}>0$ and $M_{\alpha}\in\mathbb{R}_{+}$ (which depend on
$(q,p)$) such that

$\forall\Delta t\in(0,\Delta t_{\star}],\qquad\left\lVert
c(q,p)\right\rVert\leqslant M_{\alpha}\Delta t^{\alpha}.$ (72)

In all this section, we fix $(q^{n},p^{n})\in\mathcal{A}_{\Delta t}$.

####  6.3.1 GSV numerical scheme

As $\nabla\left(\det\,D\right)=\left(\det\,D\right)\Tr\left(D^{-1}\nabla
D\right)$, the GSV numerical scheme for the Hamiltonian function (9) writes,
since $(q^{n},p^{n})\in\mathcal{A}_{\Delta t}$,

$\left\\{\begin{aligned} \widetilde{p}^{n+1/2}&=p^{n}-\dfrac{\Delta
t}{2}\left(-F(q^{n})+\dfrac{1}{2}\left(\widetilde{p}^{n+1/2}\right)^{\mathsf{T}}\nabla
D(q^{n})\widetilde{p}^{n+1/2}\right),\\\
\widetilde{q}^{n+1}&=q^{n}+\dfrac{\Delta
t}{2}\left(D\left(\widetilde{q}^{n+1}\right)+D(q^{n})\right)\widetilde{p}^{n+1/2},\\\
\widetilde{p}^{n+1}&=\widetilde{p}^{n+1/2}-\dfrac{\Delta
t}{2}\left(-F\left(\widetilde{q}^{n+1}\right)+\dfrac{1}{2}\left(\widetilde{p}^{n+1/2}\right)^{\mathsf{T}}\nabla
D\left(\widetilde{q}^{n+1}\right)\widetilde{p}^{n+1/2}\right).\end{aligned}\right.$

Therefore,

$\displaystyle\widetilde{q}^{n+1}$ $\displaystyle=q^{n}+\dfrac{\Delta
t}{2}\left(D(q^{n})+D\left(\widetilde{q}^{n+1}\right)\right)\widetilde{p}^{n+1/2}$
$\displaystyle=q^{n}+\dfrac{\Delta
t}{2}\left(D(q^{n})+D\left(\widetilde{q}^{n+1}\right)\right)\left[p^{n}-\dfrac{\Delta
t}{2}\left(-F(q^{n})+\dfrac{\left(\widetilde{p}^{n+1/2}\right)^{\mathsf{T}}\nabla
D(q^{n})\widetilde{p}^{n+1/2}}{2}\right)\right].$ (73)

We expand the latter equality in powers of $\Delta t$ up to order 4 by
consecutively looking at the terms of the right-hand side of (73). We write

$D\left(\widetilde{q}^{n+1}\right)=D(q^{n})+\nabla
D(q^{n})\odot\left(\widetilde{q}^{n+1}-q^{n}\right)+\dfrac{1}{2}\nabla^{2}D(q^{n})\odot\left(\widetilde{q}^{n+1}-q^{n}\right)^{\otimes
2}+\mathrm{O}(\Delta t^{3}),$ (74)

where we use the following notation: for any $u\in\mathbb{R}^{m}$,

$\displaystyle\nabla D(q)\odot u=\displaystyle{\sum_{i=1}^{m}\dfrac{\partial
D}{\partial q_{i}}(q)u_{i}}\in\mathbb{R}^{m\times m},$
$\displaystyle\nabla^{2}D(q)\odot u^{\otimes
2}=\displaystyle{\sum_{i,j=1}^{m}}\dfrac{\partial^{2}D}{\partial q_{i}\partial
q_{j}}(q)u_{i}u_{j}\in\mathbb{R}^{m\times m}.$

Note that the remainder in (74) satisfies (72): this follows from Taylor’s
formula with integral remainder together with the implicit function theorem to
obtain estimates on $\left\lVert\widetilde{q}^{n+1}-q^{n}\right\rVert$ as done
in [16]. The latter estimate relies on the fact that
$(\widetilde{q}^{n+1},\widetilde{p}^{n+1})$ is close to $(q^{n},p^{n})$ (see
(36) and Assumption 3 since $D$ and its derivatives grow at most polynomially.
This reasoning is used everytime we write $\mathrm{O}(\Delta t^{\alpha})$
below.

One thus obtains

$\displaystyle\dfrac{\Delta
t}{2}\left(D(q^{n})+D\left(\widetilde{q}^{n+1}\right)\right)p^{n}$
$\displaystyle=\Delta t~{}D(q^{n})p^{n}$ $\displaystyle\quad+\dfrac{\Delta
t^{2}}{4}\left[\nabla
D(q^{n})\odot\left(\left(D(q^{n})+D(\widetilde{q}^{n+1})\right)\widetilde{p}^{n+1/2}\right)\right]p^{n}$
(75) $\displaystyle\quad+\dfrac{\Delta
t^{3}}{4}\left[\nabla^{2}D(q^{n})\odot\left(D(q^{n})p^{n}\right)^{\otimes
2}\right]p^{n}$ $\displaystyle\quad+\mathrm{O}(\Delta t^{4}).$

Further expanding the second term on the right-hand side of (75) we obtain

$\displaystyle\dfrac{\Delta t^{2}}{4}\left\\{\nabla
D(q^{n})\odot\left[\left(D(q^{n})+D\left(\widetilde{q}^{n+1}\right)\right)\widetilde{p}^{n+1/2}\right]\right\\}p^{n}$
$\displaystyle=\dfrac{\Delta t^{2}}{4}\Bigg{\\{}\nabla
D(q^{n})\odot\left[2D(q^{n})+\nabla
D(q^{n})\odot\left(\widetilde{q}^{n+1}-q^{n}\right)+\mathrm{O}(\Delta
t^{2})\right]$ $\displaystyle\quad\times\left[p^{n}-\dfrac{\Delta
t}{2}\left(-F(q^{n})+\dfrac{\left(\widetilde{p}^{n+1/2}\right)^{\mathsf{T}}\nabla
D(q^{n})\widetilde{p}^{n+1/2}}{2}\right)\right]\Bigg{\\}}p^{n},$
$\displaystyle\begin{aligned} =~{}&\dfrac{\Delta t^{2}}{2}\left(\nabla
D(q^{n})\odot\left[D(q^{n})p^{n}\right]\right)p^{n}\\\ &-\dfrac{\Delta
t^{3}}{4}\left\\{\nabla
D(q^{n})\odot\left[D(q^{n})\left(-F(q^{n})+\dfrac{\left(p^{n}\right)^{\mathsf{T}}\nabla
D(q^{n})p^{n}}{2}\right)\right]\right\\}p^{n}\\\ &+\dfrac{\Delta
t^{3}}{4}\Bigg{\\{}\nabla D(q^{n})\odot\Big{[}\big{(}\nabla
D(q^{n})\odot\left[D(q^{n})p^{n}\right]\big{)}p^{n}\Big{]}\Bigg{\\}}p^{n}\\\
&+\mathrm{O}(\Delta t^{4}).\end{aligned}$

Note that $\nabla D=\nabla D^{\mathsf{T}}$ (since $D(q)$ is a symmetric
matrix) so that

$\displaystyle\left(\widetilde{p}^{n+1/2}\right)^{\mathsf{T}}\nabla
D(q^{n})\widetilde{p}^{n+1/2}$
$\displaystyle=\left(p^{n}\right)^{\mathsf{T}}\nabla D(q^{n})p^{n}$
$\displaystyle\quad-\Delta t\left(p^{n}\right)^{\mathsf{T}}\nabla
D(q^{n})\left[-F(q^{n})+\dfrac{\left(p^{n}\right)^{\mathsf{T}}\nabla
D(q^{n})p^{n}}{2}\right]+\mathrm{O}(\Delta t^{2}).$

The last term in (73) can therefore be expanded as

$\displaystyle\dfrac{\Delta t^{2}}{4}$
$\displaystyle\left(D(q^{n})+D\left(\widetilde{q}^{n+1}\right)\right)\left(-F(q^{n})+\dfrac{\left(\widetilde{p}^{n+1/2}\right)^{\mathsf{T}}\nabla
D(q^{n})\widetilde{p}^{n+1/2}}{2}\right)$ $\displaystyle=\dfrac{\Delta
t^{2}}{2}D(q^{n})\left(-F(q^{n})+\dfrac{\left(p^{n}\right)^{\mathsf{T}}\nabla
D(q^{n})p^{n}}{2}\right)$ $\displaystyle\quad-\dfrac{\Delta
t^{3}}{4}D(q^{n})\left[\left(p^{n}\right)^{\mathsf{T}}\nabla
D(q^{n})\left(-F(q^{n})+\dfrac{\left(p^{n}\right)^{\mathsf{T}}\nabla
D(q^{n})p^{n}}{2}\right)\right]$ $\displaystyle\quad+\dfrac{\Delta
t^{3}}{4}\left[\nabla
D(q^{n})\odot\left(D(q^{n})p^{n}\right)\right]\left[-F(q^{n})+\dfrac{\left(p^{n}\right)^{\mathsf{T}}\nabla
D(q^{n})p^{n}}{2}\right]$ (76) $\displaystyle\quad+\mathrm{O}(\Delta t^{4}).$

By gathering the previous equalities (75) and (76) in (73),

$\displaystyle\widetilde{q}^{n+1}$ $\displaystyle=q^{n}+\Delta
t\,D(q^{n})p^{n}$ $\displaystyle\quad+\dfrac{\Delta
t^{2}}{2}\left\\{\left[\nabla
D(q^{n})\odot\left[D(q^{n})p^{n}\right]\right]p^{n}-D(q^{n})\left[-F(q^{n})+\dfrac{\left(p^{n}\right)^{\mathsf{T}}\nabla
D(q^{n})p^{n}}{2}\right]\right\\}$ $\displaystyle\quad+\Delta
t^{3}\left[\begin{array}[]{l}-\left\\{\dfrac{\nabla
D(q^{n})}{4}\odot\left[D(q^{n})\left(-F(q^{n})+\dfrac{\left(p^{n}\right)^{\mathsf{T}}\nabla
D(q^{n})p^{n}}{2}\right)\right]\right\\}p^{n}\\\\[8.5359pt]
+\left\\{\dfrac{\nabla D(q^{n})}{4}\odot\left[\left(\nabla
D(q^{n})\odot\left[D(q^{n})p^{n}\right]\right)p^{n}\right]\right\\}p^{n}\\\\[8.5359pt]
+\dfrac{D(q^{n})}{4}\left\\{\left(p^{n}\right)^{\mathsf{T}}\nabla
D(q^{n})\left[-F(q^{n})+\dfrac{\left(p^{n}\right)^{\mathsf{T}}\nabla
D(q^{n})p^{n}}{2}\right]\right\\}\\\\[8.5359pt]
+\dfrac{1}{4}\left\\{\nabla^{2}D(q^{n})\odot\left(D(q^{n})p^{n}\right)^{\otimes
2}\right\\}p^{n}\\\\[8.5359pt] -\dfrac{1}{4}\left[\nabla
D(q^{n})\odot\left[D(q^{n})p^{n}\right]\right]\left[-F(q^{n})+\dfrac{\left(p^{n}\right)^{\mathsf{T}}\nabla
D(q^{n})p^{n}}{2}\right]\end{array}\right]$ (82)
$\displaystyle\quad+\mathrm{O}(\Delta t^{4}),$ $\displaystyle=q^{n}+\Delta
t\,\mathcal{F}_{1}(q^{n},p^{n})+\dfrac{\Delta
t^{2}}{2}\mathcal{F}_{2}(q^{n},p^{n})+\Delta
t^{3}\mathcal{F}_{3}(q^{n},p^{n})+\mathrm{O}(\Delta t^{4}),$ (83)

where the functions $\mathcal{F}_{1},\mathcal{F}_{2},\mathcal{F}_{3}$ are
defined by comparing the power expansion in $\Delta t$ between the first and
second equality of (83). The remainder term in (83) satisfies (38), which
shows that (37) holds for the GSV numerical scheme.

Since $p^{n}\sim\mathcal{N}(0,D(q^{n})^{-1})$ in step [2.i] of Algorithm 2, we
replace the momenta by $p^{n}=D(q^{n})^{-1/2}\mathrm{G}^{n}$ where
$\mathrm{G}^{n}\sim\mathcal{N}(0,\mathrm{I}_{m})$. Then,

$\left\\{\begin{aligned}
\mathcal{F}_{1}(q^{n},D(q^{n})^{-1/2}\mathrm{G}^{n})&=D(q^{n})^{1/2}\mathrm{G}^{n},\\\
\mathcal{F}_{2}(q^{n},D(q^{n})^{-1/2}\mathrm{G}^{n})&=-D(q^{n})\nabla
V(q^{n})+\dfrac{D(q^{n})}{2}\Tr\left(D(q^{n})^{-1}\nabla D(q^{n})\right)\\\
&\phantom{=~{}}-\dfrac{D(q^{n})}{2}\left[\left[D(q^{n})^{-1/2}\mathrm{G}^{n}\right]^{\mathsf{T}}\nabla
D(q^{n})\left[D(q^{n})^{-1/2}\mathrm{G}^{n}\right]\right]\\\
&\phantom{=~{}}+\left(\nabla
D(q^{n})\odot\left[D(q^{n})^{1/2}\mathrm{G}^{n}\right]\right)D(q^{n})^{-1/2}\mathrm{G}^{n}.\end{aligned}\right.$

The functions $\mathcal{F}_{1}$ and $\mathcal{F}_{3}$ being odd in $p^{n}$
hence in $\mathrm{G}^{n}$, and the function $\mathcal{F}_{2}$ being even in
$p^{n}$ hence in $\mathrm{G}^{n}$, it holds

$\left\\{\begin{aligned}
\mathbb{E}_{\mathrm{G}^{n}}\left[\mathcal{F}_{1}(q^{n},D(q^{n})^{-1/2}\mathrm{G}^{n})\right]&=0,\\\
\mathbb{E}_{\mathrm{G}^{n}}\left[\mathcal{F}_{3}(q^{n},D(q^{n})^{-1/2}\mathrm{G}^{n})\right]&=0,\\\
\mathbb{E}_{\mathrm{G}^{n}}\left[\mathcal{F}_{1}(q^{n},D(q^{n})^{-1/2}\mathrm{G}^{n})\otimes\mathcal{F}_{2}(q^{n},D(q^{n})^{-1/2}\mathrm{G}^{n})\right]&=0,\end{aligned}\right.$
(84)

where here $\otimes$ simply denotes the outer product. Finally, to show that
(39) holds, let us compute the various terms in
$\mathbb{E}_{\mathrm{G}^{n}}\left[\mathcal{F}_{2}(q^{n},D(q^{n})^{-1/2}\mathrm{G}^{n}))\right]$.
For $1\leqslant i\leqslant m$, it holds

$\displaystyle\mathbb{E}_{\mathrm{G}^{n}}\left(\left[\left\\{\nabla
D(q^{n})\odot\left(D(q^{n})^{1/2}\mathrm{G}^{n}\right)\right\\}D(q^{n})^{-1/2}\mathrm{G}^{n}\right]_{i}\right)$
$\displaystyle\quad=\mathbb{E}_{\mathrm{G}^{n}}\left(\sum_{k=1}^{m}\left[\nabla
D(q^{n})\odot\left(D(q^{n})^{1/2}\mathrm{G}^{n}\right)\right]_{i,k}\left[D(q^{n})^{-1/2}\mathrm{G}^{n}\right]_{k}\right)$
$\displaystyle\quad=\mathbb{E}_{\mathrm{G}^{n}}\left(\sum_{k=1}^{m}\sum_{j=1}^{m}\partial_{q_{j}}D_{i,k}(q^{n})\left[D(q^{n})^{1/2}\mathrm{G}^{n}\right]_{j}\left[D(q^{n})^{-1/2}\mathrm{G}^{n}\right]_{k}\right)$
$\displaystyle\quad=\sum_{k=1}^{m}\sum_{j=1}^{m}\partial_{q_{j}}D_{i,k}(q^{n})\mathbb{E}_{\mathrm{G}^{n}}\left(\left[D(q^{n})^{1/2}\mathrm{G}^{n}\right]_{j}\left[D(q^{n})^{-1/2}\mathrm{G}^{n}\right]_{k}\right)$
$\displaystyle\quad=\sum_{k=1}^{m}\sum_{j=1}^{m}\partial_{q_{j}}D_{i,k}(q^{n})\underbrace{\sum_{\ell=1}^{m}\left[D(q^{n})^{1/2}\right]_{j,\ell}\left[D(q^{n})^{-1/2}\right]_{k,\ell}}_{=~{}\left[D(q^{n})^{1/2}D(q^{n})^{-\mathsf{T}/2}\right]_{j,k}=~{}\delta_{j,k}}$
$\displaystyle\quad=\sum_{k=1}^{m}\partial_{q_{k}}D_{i,k}(q^{n})=\left[\operatorname{div}D(q^{n})\right]_{i}.$

Additionally,

$\displaystyle\mathbb{E}_{\mathrm{G}^{n}}\left(\left[D(q^{n})\left\\{\left(D(q^{n})^{-1/2}\mathrm{G}^{n}\right)^{\mathsf{T}}\nabla
D(q^{n})\left(D(q^{n})^{-1/2}\mathrm{G}^{n}\right)\right\\}\right]_{i}\right)$
$\displaystyle\quad=\mathbb{E}_{\mathrm{G}^{n}}\left(\sum_{k=1}^{m}D_{i,k}(q^{n})\left[\left(D(q^{n})^{-1/2}\mathrm{G}^{n}\right)^{\mathsf{T}}\nabla
D(q^{n})\left(D(q^{n})^{-1/2}\mathrm{G}^{n}\right)\right]_{k}\right)$
$\displaystyle\quad=\mathbb{E}_{\mathrm{G}^{n}}\left(\sum_{k=1}^{m}D_{i,k}(q^{n})\sum_{j,\ell=1}^{m}\partial_{q_{k}}D_{j,\ell}(q^{n})\left[D(q^{n})^{-1/2}\mathrm{G}^{n}\right]_{j}\left[D(q^{n})^{-1/2}\mathrm{G}^{n}\right]_{\ell}\right)$
$\displaystyle\quad=\sum_{k=1}^{m}D_{i,k}(q^{n})\sum_{j,\ell=1}^{m}\partial_{q_{k}}D_{j,\ell}(q^{n})\underbrace{\mathbb{E}_{\mathrm{G}^{n}}\left[\left[D(q^{n})^{-1/2}\mathrm{G}^{n}\right]_{j}\left[D(q^{n})^{-\mathsf{T}/2}\mathrm{G}^{n}\right]_{\ell}\right]}_{=~{}\sum_{\alpha=1}^{m}\left[D(q^{n})^{-1/2}\right]_{j,\alpha}\left[D(q^{n})^{-\mathsf{T}/2}\right]_{\alpha,\ell}=~{}\left[D(q^{n})^{-1}\right]_{j,\ell}}$
$\displaystyle\quad=\sum_{k=1}^{m}D_{i,k}(q^{n})\sum_{j,\ell=1}^{m}\partial_{q_{k}}\left[D(q^{n})\right]_{j,\ell}\left[D(q^{n})^{-1}\right]_{j,\ell}.$

By noticing that

$\displaystyle\left[D(q^{n})\Tr\left(D(q^{n})^{-1}\nabla
D(q^{n})\right)\right]_{i}$
$\displaystyle=\sum_{k=1}^{m}D_{i,k}(q^{n})\left[\Tr\left(D(q^{n})^{-1}\nabla
D(q^{n})\right)\right]_{k}$
$\displaystyle=\sum_{k=1}^{m}D_{i,k}(q^{n})\Tr\left(D(q^{n})^{-1}\partial_{q_{k}}D(q^{n})\right)$
$\displaystyle=\sum_{k=1}^{m}D_{i,k}(q^{n})\sum_{j,\ell=1}^{m}\partial_{q_{k}}\left[D(q^{n})\right]_{j,\ell}\left[D(q^{n})^{-1}\right]_{j,\ell},$

we thus conclude that

$\mathbb{E}_{\mathrm{G}^{n}}\left[\mathcal{F}_{2}(q^{n},D(q^{n})^{-1/2}\mathrm{G}^{n})\right]=-D(q^{n})\nabla
V(q^{n})+\operatorname{div}D(q^{n}),$

which is exactly (39). This concludes the proof of Lemma 10 for the GSV
numerical scheme.

####  6.3.2 IMR numerical scheme

The IMR numerical scheme for the Hamiltonian function (9) writes

$\left\\{\begin{aligned} \widetilde{q}^{n+1}&=q^{n}+\Delta
tD\left(\frac{q^{n}+\widetilde{q}^{n+1}}{2}\right)\left(\frac{p^{n}+\widetilde{p}^{n+1}}{2}\right),\\\
\widetilde{p}^{n+1}&=p^{n}-\Delta
t\left(-F\left(\frac{q^{n}+\widetilde{q}^{n+1}}{2}\right)+\frac{1}{2}\left(\frac{p^{n}+\widetilde{p}^{n+1}}{2}\right)^{\mathsf{T}}\nabla
D\left(\frac{q^{n}+\widetilde{q}^{n+1}}{2}\right)\left(\frac{p^{n}+\widetilde{p}^{n+1}}{2}\right)\right).\end{aligned}\right.$

Introduce $y^{n}=(q^{n},p^{n})$. We start by rewriting the update of the
position as

$\displaystyle\widetilde{q}^{n+1}$ $\displaystyle=q^{n}+\Delta
t~{}D\left(\dfrac{q^{n}+\widetilde{q}^{n+1}}{2}\right)\left(\dfrac{p^{n}+\widetilde{p}^{n+1}}{2}\right),$
$\displaystyle=q^{n}+A(\Delta t,y^{n},\widetilde{y}^{n+1})-B(\Delta
t,y^{n},\widetilde{y}^{n+1}),$

where the term $A(\Delta t,y^{n},\widetilde{y}^{n+1})$ is

$\displaystyle A(\Delta t,y^{n},\widetilde{y}^{n+1})$ $\displaystyle=\Delta
t~{}D\left(\dfrac{q^{n}+\widetilde{q}^{n+1}}{2}\right)p^{n}$
$\displaystyle=\Delta t~{}D(q^{n})p^{n}$ $\displaystyle\quad+\dfrac{\Delta
t^{2}}{2}\left\\{\nabla
D(q^{n})\odot\left[D\left(\dfrac{q^{n}+\widetilde{q}^{n+1}}{2}\right)\left(\dfrac{p^{n}+\widetilde{p}^{n+1}}{2}\right)\right]\right\\}p^{n}$
$\displaystyle\quad+\dfrac{\Delta
t^{3}}{8}\left[\nabla^{2}D(q^{n})\odot\left(D(q^{n})p^{n}\right)^{\otimes
2}\right]p^{n}$ $\displaystyle\quad+\mathrm{O}(\Delta t^{4}),$

and the term $B(\Delta t,y^{n},\widetilde{y}^{n+1})$ is

$\dfrac{\Delta
t^{2}}{2}D\left(\frac{q^{n}+\widetilde{q}^{n+1}}{2}\right)\left[-F\left(\dfrac{q^{n}+\widetilde{q}^{n+1}}{2}\right)+\dfrac{1}{2}\left(\dfrac{p^{n}+\widetilde{p}^{n+1}}{2}\right)^{\mathsf{T}}\nabla
D\left(\dfrac{q^{n}+\widetilde{q}^{n+1}}{2}\right)\left(\dfrac{p^{n}+\widetilde{p}^{n+1}}{2}\right)\right].$

An expansion in powers of $\Delta t$ of the second term in $A$ gives

$\displaystyle\dfrac{\Delta t^{2}}{2}\left\\{\nabla
D(q^{n})\odot\left[D\left(\dfrac{q^{n}+\widetilde{q}^{n+1}}{2}\right)\left(\dfrac{p^{n}+\widetilde{p}^{n+1}}{2}\right)\right]\right\\}p^{n}$
$\displaystyle=\dfrac{\Delta t^{2}}{2}\left\\{\nabla
D(q^{n})\odot\left[D(q^{n})p^{n}\right]\right\\}p^{n}$
$\displaystyle\quad+\dfrac{\Delta t^{3}}{4}\left\\{\nabla
D(q^{n})\odot\left[\left(\nabla
D(q^{n})\odot\left[D(q^{n})p^{n}\right]\right)p^{n}\right]\right\\}p^{n}$
$\displaystyle\quad-\dfrac{\Delta t^{3}}{4}\left\\{\nabla
D(q^{n})\odot\left[D(q^{n})\left(-F(q^{n})+\left(p^{n}\right)^{\mathsf{T}}\dfrac{\nabla
D(q^{n})}{2}p^{n}\right)\right]\right\\}p^{n}$
$\displaystyle\quad+\mathrm{O}(\Delta t^{4}).$

To expand the terms present in $B(\Delta t,y^{n},\widetilde{y}^{n+1})$, we
need to expand $F\left(\dfrac{q^{n}+\widetilde{q}^{n+1}}{2}\right)$ in powers
of $\Delta t$. As $F(q)=-\nabla V(q)+\Tr(D(q)^{-1}\nabla D(q))/2$, we first
write

$\nabla V\left(\frac{q^{n}+\widetilde{q}^{n+1}}{2}\right)=\nabla
V(q^{n})+\frac{\Delta t}{2}\nabla^{2}V(q^{n})D(q^{n})p^{n}+\mathrm{O}(\Delta
t^{2}).$

As for the second term
$\Tr\left[D\left(\dfrac{q^{n}+\widetilde{q}^{n+1}}{2}\right)^{-1}\nabla
D\left(\dfrac{q^{n}+\widetilde{q}^{n+1}}{2}\right)\right]$, it holds

$\displaystyle D\left(\dfrac{q^{n}+\widetilde{q}^{n+1}}{2}\right)^{-1}$
$\displaystyle=D(q^{n})^{-1}-D(q^{n})^{-1}\left[\nabla
D(q^{n})\odot\left[\dfrac{\widetilde{q}^{n+1}-q^{n}}{2}\right]\right]D(q^{n})^{-1}+\mathrm{O}(\Delta
t^{2})$ $\displaystyle=D(q^{n})^{-1}-\dfrac{\Delta
t}{2}D(q^{n})^{-1}\left[\nabla
D(q^{n})\odot\left(D(q^{n})p^{n}\right)\right]D(q^{n})^{-1}+\mathrm{O}(\Delta
t^{2}),$

and

$\nabla D\left(\dfrac{q^{n}+\widetilde{q}^{n+1}}{2}\right)=\nabla
D(q^{n})+\dfrac{\Delta
t}{2}\nabla^{2}D(q^{n})\odot\left[D(q^{n})p^{n}\right]+\mathrm{O}(\Delta
t^{2}).$

Therefore,

$\displaystyle\dfrac{\Delta
t^{2}}{2}D\left(\dfrac{q^{n}+\widetilde{q}^{n+1}}{2}\right)\left[-F\left(\dfrac{q^{n}+\widetilde{q}^{n+1}}{2}\right)+\dfrac{1}{2}\left(\dfrac{p^{n}+\widetilde{p}^{n+1}}{2}\right)^{\mathsf{T}}\nabla
D\left(\dfrac{q^{n}+\widetilde{q}^{n+1}}{2}\right)\left(\dfrac{p^{n}+\widetilde{p}^{n+1}}{2}\right)\right]$
$\displaystyle=\dfrac{\Delta
t^{2}}{2}D(q^{n})\left(-F(q^{n})+\dfrac{1}{2}\left(p^{n}\right)^{\mathsf{T}}\nabla
D(q^{n})p^{n}\right)$ $\displaystyle\quad+\dfrac{\Delta
t^{3}}{4}\left[\Bigl{[}\nabla
D(q^{n})\odot\left[D(q^{n})p^{n}\right]\Bigr{]}\left(-F(q^{n})+\dfrac{1}{2}\left(p^{n}\right)^{\mathsf{T}}\nabla
D(q^{n})p^{n}\right)\right]$ $\displaystyle\quad+\dfrac{\Delta
t^{3}}{4}D(q^{n})\nabla^{2}V(q^{n})D(q^{n})p^{n}$
$\displaystyle\quad-\dfrac{\Delta
t^{3}}{8}D(q^{n})\Tr\left[D(q^{n})^{-1}\Bigl{(}\nabla^{2}D(q^{n})\odot\left[D(q^{n})p^{n}\right]\Bigr{)}\right]$
$\displaystyle\quad+\dfrac{\Delta
t^{3}}{8}D(q^{n})\Tr\left[\Bigl{(}D(q^{n})^{-1}\left[\nabla
D(q^{n})\odot\left(D(q^{n})p^{n}\right)\right]D(q^{n})^{-1}\Bigr{)}\nabla
D(q^{n})\right]$ $\displaystyle\quad-\dfrac{\Delta
t^{3}}{4}D(q^{n})\left[\left(p^{n}\right)^{\mathsf{T}}\nabla
D(q^{n})\left(-F(q^{n})+\dfrac{1}{2}\left(p^{n}\right)^{\mathsf{T}}\nabla
D(q^{n})p^{n}\right)\right]$ $\displaystyle\quad+\dfrac{\Delta
t^{3}}{8}D(q^{n})\left[\left(p^{n}\right)^{\mathsf{T}}\Bigl{\\{}\nabla^{2}D(q^{n})\odot\left[D(q^{n})p^{n}\right]\Bigr{\\}}p^{n}\right]$
$\displaystyle\quad+\mathrm{O}(\Delta t^{4}).$

By gathering the previous equalities, one finally obtains

$\displaystyle\widetilde{q}^{n+1}$ $\displaystyle=q^{n}+\Delta
t\,D(q^{n})p^{n}$ $\displaystyle\quad+\dfrac{\Delta
t^{2}}{2}\left\\{\left[\nabla
D(q^{n})\odot\left[D(q^{n})p^{n}\right]\right]p^{n}-D(q^{n})\left[-F(q^{n})+\dfrac{\left(p^{n}\right)^{\mathsf{T}}\nabla
D(q^{n})p^{n}}{2}\right]\right\\}$ $\displaystyle\quad+\Delta
t^{3}\left[\begin{array}[]{l}+\dfrac{1}{4}\left\\{\nabla
D(q^{n})\odot\left[\left(\nabla
D(q^{n})\odot\left[D(q^{n})p^{n}\right]\right)p^{n}\right]\right\\}p^{n}\\\\[5.69046pt]
-\dfrac{1}{4}\left\\{\nabla
D(q^{n})\odot\left[D(q^{n})\left(-F(q^{n})+\left(p^{n}\right)^{\mathsf{T}}\dfrac{\nabla
D(q^{n})}{2}p^{n}\right)\right]\right\\}p^{n}\\\\[5.69046pt]
+\dfrac{1}{4}\left[\nabla
D(q^{n})\odot\left[D(q^{n})p^{n}\right]\right]\left(-F(q^{n})+\left(p^{n}\right)^{\mathsf{T}}\dfrac{\nabla
D(q^{n})}{2}p^{n}\right)\\\\[5.69046pt]
+\dfrac{1}{4}D(q^{n})\nabla^{2}V(q^{n})D(q^{n})p^{n}\\\\[5.69046pt]
-\dfrac{1}{8}D(q^{n})\Tr\left[D(q^{n})^{-1}\left(\nabla^{2}D(q^{n})\odot\left[D(q^{n})p^{n}\right]\right)\right]\\\\[5.69046pt]
+\dfrac{1}{8}D(q^{n})\Tr\left[\left(D(q^{n})^{-1}\left[\nabla
D(q^{n})\odot\left(D(q^{n})p^{n}\right)\right]D(q^{n})^{-1}\right)\nabla
D(q^{n})\right]\\\\[5.69046pt]
-\dfrac{1}{4}D(q^{n})\left[\left(p^{n}\right)^{\mathsf{T}}\nabla
D(q^{n})\left(-F(q^{n})+\left(p^{n}\right)^{\mathsf{T}}\dfrac{\nabla
D(q^{n})}{2}p^{n}\right)\right]\\\\[5.69046pt]
+\dfrac{1}{8}D(q^{n})\left[\left(p^{n}\right)^{\mathsf{T}}\left\\{\nabla^{2}D(q^{n})\odot\left[D(q^{n})p^{n}\right]\right\\}p^{n}\right]\\\\[5.69046pt]
+\dfrac{1}{8}\left[\nabla^{2}D(q^{n})\odot\left(D(q^{n})p^{n}\right)^{\otimes
2}\right]p^{n}\end{array}\right]$ $\displaystyle\quad+\mathrm{O}(\Delta
t^{4})$ $\displaystyle=q^{n}+\Delta
t\,\mathcal{F}_{1}(q^{n},p^{n})+\dfrac{\Delta
t^{2}}{2}\mathcal{F}_{2}(q^{n},p^{n})+\Delta
t^{3}\mathcal{F}_{3}(q^{n},p^{n})+\mathrm{O}(\Delta t^{4}),$

where $\mathcal{F}_{1}$ and $\mathcal{F}_{2}$ are the same functions as in
(83), whereas $\mathcal{F}_{3}$ is different but still odd in
$\mathrm{G}^{n}$. In particular,
$\mathbb{E}_{\mathrm{G}^{n}}\left[\mathcal{F}_{3}(q^{n},D(q^{n})^{-1/2}\mathrm{G}^{n})\right]=0$
and
$\mathbb{E}_{\mathrm{G}^{n}}\left[\mathcal{F}_{2}(q^{n},D(q^{n})^{-1/2}\mathrm{G}^{n})\right]=-D(q^{n})\nabla
V(q^{n})+\operatorname{div}D(q^{n})$ following the same computations as in the
previous section. This concludes the proof for the IMR numerical scheme.

###  6.4 Proof of Proposition 11

First note that, if $f:\mathcal{X}\to\mathcal{X}$ grows at most polynomially,
with derivatives growing at most polynomially, then (35) implies that there
exist $\alpha\in\mathbb{N}$ and $K\in\mathbb{R}_{+}$ such that

$\forall q\in\mathcal{O},\qquad\left\lvert\Delta
t^{\delta}f\left(q,D(q)^{-1/2}\mathrm{G}\right)\right\rvert\leqslant K\Delta
t^{\delta}\left(1+\left\lVert\mathrm{G}\right\rVert^{\alpha}+\left\lVert
q\right\rVert^{\alpha}\right).$

We write such remainder terms as $\Delta
t^{\delta}f\left(q,D(q)^{-1/2}\mathrm{G}\right)=\mathrm{O}(\Delta
t^{\delta})$. In particular,

$\mathbb{E}_{\mathrm{G}}\left[\left\lvert\Delta
t^{\delta}f\left(q,D(q)^{-1/2}\mathrm{G}\right)\right\rvert^{r}\right]\leqslant
K_{r}\Delta t^{\delta r}\left(1+\left\lVert q\right\rVert^{\alpha r}\right),$

with $K_{r}\in\mathbb{R}_{+}$ for any $r\geqslant 0$, which we may also write
$\mathbb{E}_{\mathrm{G}}\left[\left\lvert\Delta
t^{\delta}f\left(q,D(q)^{-1/2}\mathrm{G}\right)\right\rvert^{r}\right]=\mathrm{O}\left(\Delta
t^{\delta r}\right)$. We now prove that, for any $\Delta t\in(0,\Delta
t_{\star}]$, it holds

$\mathbb{E}_{\mathrm{G}}\left[f\left(q,D(q)^{-1/2}\mathrm{G}\right)\right]=\mathbb{E}_{\mathrm{G}}\left[f\left(q,D(q)^{-1/2}\mathrm{G}\right)\mathbf{1}_{\left(q,D(q)^{-1/2}\mathrm{G}\right)\in\mathcal{B}_{\Delta
t}}\right]+\mathrm{O}\left(\Delta t^{2(r_{3}+1)}\right).$ (85)

Indeed, consider $K\in\mathbb{R}_{+}$ and $\alpha\in\mathbb{N}$ such that, for
all $q\in\mathcal{O}$, it holds $\left\lvert
f\left(q,D(q)^{-1/2}\mathrm{G}\right)\right\rvert\leqslant K(1+\left\lVert
q\right\rVert^{\alpha}+\left\lVert\mathrm{G}\right\rVert^{\alpha})$. Then,
using (42) with $a=0$ and $a=\alpha$,

$\displaystyle\left\lvert\mathbb{E}_{\mathrm{G}}\left[f\left(q,D(q)^{-1/2}\mathrm{G}\right)\mathbf{1}_{\left(q,D(q)^{-1/2}\mathrm{G}\right)\in\mathcal{X}\setminus\mathcal{B}_{\Delta
t}}\right]\right\rvert$ $\displaystyle\quad\leqslant K(1+\left\lVert
q\right\rVert^{\alpha})\mathbb{E}_{\mathrm{G}}\left[\mathbf{1}_{\left(q,D(q)^{-1/2}\mathrm{G}\right)\in\mathcal{X}\setminus\mathcal{B}_{\Delta
t}}\right]+K\mathbb{E}_{\mathrm{G}}\left[\left\lVert\mathrm{G}\right\rVert^{\alpha}\mathbf{1}_{\left(q,D(q)^{-1/2}\mathrm{G}\right)\in\mathcal{X}\setminus\mathcal{B}_{\Delta
t}}\right]$ $\displaystyle\quad\leqslant K(1+\left\lVert
q\right\rVert^{\alpha})C_{0}\Delta
t^{2(r_{3}+1)}+KC_{\alpha}\left(1+\left\lVert
q\right\rVert^{b_{\alpha}}\right)\Delta t^{2(r_{3}+1)}=\mathrm{O}\left(\Delta
t^{2(r_{3}+1)}\right).$

We now fix $f\in\mathscr{F}_{\mathrm{pol}}$ and $q\in\mathcal{O}$. The
discrete position evolution operator $P_{\Delta t}$ using one step of
Algorithm 2 with either the IMR or the GSV numerical scheme acts on $f$ as

$\left(P_{\Delta t}f\right)(q)=\mathbb{E}^{q}\left[f(q^{1})\right].$

Denote by $A_{\Delta t}(q,\mathrm{G})$ the acceptance probability in the
Metropolis–Hastings step [2.iii], which is defined for
$\left(q,D(q)^{-1/2}\mathrm{G}\right)\in\mathcal{A}_{\Delta t}$ by

$A_{\Delta t}(q,\mathrm{G})=\min\left(1,\mathrm{e}^{-\Delta
H\left(q,\mathrm{G}\right)}\right),$

where

$\Delta H\left(q,\mathrm{G}\right)=H\left(\varphi_{\Delta
t}\left(q,D(q)^{-1/2}\mathrm{G}\right)\right)-H\left(q,D(q)^{-1/2}\mathrm{G}\right).$

If $\mathrm{G}$ is such that
$\left(q,D(q)^{-1/2}\mathrm{G}\right)\in\mathcal{X}\setminus\mathcal{B}_{\Delta
t}$ or if the proposal is rejected in the Metropolis–Hastings procedure, then
$q^{1}=q$. Therefore,

$\displaystyle(P_{\Delta t}f-f)(q)$
$\displaystyle=\mathbb{E}^{q}\left[\left(f(q^{1})-f(q)\right)\mathbb{1}_{\left(q,D(q)^{-1/2}\mathrm{G}\right)\in\mathcal{B}_{\Delta
t}}A_{\Delta t}(q,\mathrm{G})\right]$ $\displaystyle=\begin{aligned}
&\mathbb{E}^{q}\left[\left(f(q^{1})-f(q)\right)\mathbb{1}_{\left(q,D(q)^{-1/2}\mathrm{G}\right)\in\mathcal{B}_{\Delta
t}}\right]\\\ &+\mathbb{E}^{q}\left[\left(A_{\Delta
t}\left(q,\mathrm{G}\right)-1\right)\left(f(q^{1})-f(q)\right)\mathbb{1}_{\left(q,D(q)^{-1/2}\mathrm{G}\right)\in\mathcal{B}_{\Delta
t}}\right].\end{aligned}$ (86)

We first focus on the first term of (86). We expand

$\displaystyle f(q^{1})-f(q)$ $\displaystyle=\nabla
f(q)\cdot\left(q^{1}-q\right)+\frac{1}{2}\nabla^{2}f(q):\left(q^{1}-q\right)^{\otimes
2}$
$\displaystyle\quad+\frac{1}{6}\nabla^{3}f(q):\left(q^{1}-q\right)^{\otimes
3}+\frac{1}{6}\int_{0}^{1}(1-\theta)^{3}\nabla^{4}f((1-\theta)q^{1}+\theta
q):(q^{1}-q)^{\otimes 4}\mathrm{d}\theta,$

where we denoted, for any positions $(q,\tilde{q})$ in $\mathcal{O}$,

$\nabla^{3}f(q):\tilde{q}^{\otimes
3}=\sum_{\alpha,\beta,\gamma=1}^{m}\frac{\partial^{3}f}{\partial_{q_{\alpha}}\partial_{q_{\beta}}\partial_{q_{\gamma}}}(q)\,\tilde{q}_{\alpha}\tilde{q}_{\beta}\tilde{q}_{\gamma},$

and similarly for $\nabla^{4}f$.

Since $\mathcal{B}_{\Delta t}\subset\mathcal{A}_{\Delta t}$, we may plug the
expansion (37) into the first term of (86). In order to write the remainder
$\Delta t^{4}\mathcal{R}_{\Delta t}(q,p)$ of the expansion (37) as
$\mathrm{O}(\Delta t^{4})$ when integrating with respect to $\mathrm{G}$, we
need to show that $\mathcal{R}_{\Delta t}(q,p)$ is bounded by a polynomial
function of $\left\lVert q\right\rVert$ and $\left\lVert p\right\rVert$. We
show more generally that $\varphi_{\Delta t}(q,p)$ (which is the proposal of
$(q,p)$ and not simply the proposal of the position) is bounded by a
polynomial function of $\left\lVert q\right\rVert$ and $\left\lVert
p\right\rVert$. Indeed, let $(q,p)\in\mathcal{X}$ and
$g\colon\mathbb{R}_{+}\to\mathcal{X}$ defined by $g(t)=\phi_{t}(q,p)$ where
$\phi_{t}$ is the exact flow of (2) with initial condition
$(q_{0},p_{0})=(q,p)$. Note first that $g(t)$ is well defined for all times
$t\in\mathbb{R}_{+}$. Indeed, the solution to (2) is well defined for small
times by the Cauchy–Lipschitz theorem, and globally defined for all times
using the Hamiltonian function as a Lyapunov function. To make the latter
point precise, we show in fact the stronger result that the norm of
$g(s)=(q_{s},p_{s})$ is controlled by a polynomial function of $\left\lVert
q\right\rVert$ and $\left\lVert p\right\rVert$. Since the Hamiltonian is
exactly preserved by $\phi_{s}$, is holds $H(q_{s},p_{s})=H(q,p)$, which can
be rewritten as

$V(q_{s})-\frac{1}{2}\ln\det
D(q_{s})+\frac{1}{2}p_{s}^{\mathsf{T}}D(q_{s})p_{s}=V(q)-\frac{1}{2}\ln\det
D(q)+\frac{1}{2}p^{\mathsf{T}}D(q)p.$

By Assumption 3, the right-hand side is bounded by a polynomial function of
$\left\lVert q\right\rVert$ and $\left\lVert p\right\rVert$. Using that $V$ is
bounded from below and that $D$ is uniformly bounded (from above and from
below), it follows that $\left\lVert p_{s}\right\rVert$ is bounded by a
polynomial function of $\left\lVert q\right\rVert$ and $\left\lVert
p\right\rVert$. The first equation of (2), namely

$\frac{\mathrm{d}}{\mathrm{d}t}q_{t}=D(q_{t})p_{t},$

then implies that $\left\lVert q_{s}\right\|$ is bounded by a polynomial
function of $\left\lVert q\right\rVert$ and $\left\lVert p\right\rVert$
(uniformly on bounded intervals). Combined with (41), this shows that
$\varphi_{\Delta t}(q,p)$ is also bounded by a polynomial function of
$\left\lVert q\right\rVert$ and $\left\lVert p\right\rVert$. Similar
reasonings underly remainder terms appearing below.

Using (85), and in view of (84) and (39),

$\displaystyle\mathbb{E}_{\mathrm{G}}\left[\left(f(q^{1})-f(q)\right)\mathbb{1}_{\left(q,D(q)^{-1/2}\mathrm{G}\right)\in\mathcal{B}_{\Delta
t}}\right]$ $\displaystyle\quad=\frac{\Delta t^{2}}{2}\nabla
f(q)\cdot\left(-D(q)\nabla V(q)+\operatorname{div}D(q)\right)$
$\displaystyle\qquad+\frac{\Delta
t^{2}}{2}\mathbb{E}_{\mathrm{G}}\left[\nabla^{2}f(q):\left(\mathcal{F}_{1}\left(q,D(q)^{-1/2}\mathrm{G}\right)\right)^{\otimes
2}\right]+\mathrm{O}\left(\Delta t^{2(\min(1,r_{3})+1)}\right).$

The expectation in the latter term can be computed as

$\displaystyle\mathbb{E}_{\mathrm{G}}\left(\nabla^{2}f(q):\left(\mathcal{F}_{1}\left(q,D(q)^{-1/2}\mathrm{G}\right)\right)^{\otimes
2}\right)$
$\displaystyle=\mathbb{E}_{\mathrm{G}}\left(\sum_{i,j=1}^{m}\frac{\partial^{2}f}{\partial
q_{i}\partial
q_{j}}(q)\left[D(q)^{1/2}\mathrm{G}\right]_{i}\left[D(q)^{1/2}\mathrm{G}\right]_{j}\right)$
$\displaystyle=\mathbb{E}_{\mathrm{G}}\left(\sum_{i,j=1}^{m}\frac{\partial^{2}f}{\partial
q_{i}\partial
q_{j}}(q)\sum_{k,\ell=1}^{m}\left[D(q)^{1/2}\right]_{i,k}\mathrm{G}_{k}\left[D(q)^{1/2}\right]_{j,\ell}\mathrm{G}_{\ell}\right)$
$\displaystyle=\sum_{i,j=1}^{m}\frac{\partial^{2}f}{\partial q_{i}\partial
q_{j}}(q)\sum_{k=1}^{m}\left[D(q)^{1/2}\right]_{i,k}\left[D(q)^{1/2}\right]_{j,k}$
$\displaystyle=\sum_{i,j=1}^{m}\frac{\partial^{2}f}{\partial q_{i}\partial
q_{j}}(q)\left[D(q)\right]_{i,j}=D(q):\nabla^{2}f(q).$

This shows that

$\mathbb{E}_{\mathrm{G}}\left[\left(f(q^{1})-f(q)\right)\mathbb{1}_{\left(q,D(q)^{-1/2}\mathrm{G}\right)\in\mathcal{B}_{\Delta
t}}\right]=\frac{\Delta t^{2}}{2}\mathcal{L}f(q)+\mathrm{O}\left(\Delta
t^{2(\min(1,r_{3})+1)}\right),$

where $\mathcal{L}$ is the generator of the dynamics (33) defined by (34).

We now focus on the second term of (86). In view of (41), it holds
$\varphi_{\Delta t}(q,p)=\phi_{\Delta t}(q,p)+\Delta
t^{3}\widetilde{R}_{\Delta t}(q,p)$, where $\widetilde{R}_{\Delta t}$ is a
function growing at most polynomially. Since the Hamiltonian is exactly
preserved by $\phi_{t}$, there exists, for any $\Delta t\in(0,\Delta
t_{\star}]$ and $(q,p)\in\mathcal{A}_{\Delta t}$, a real number
$\theta_{\Delta t}(q,p)\in[0,1]$ such that

$\displaystyle H(\varphi_{\Delta t}(q,p))=H(q,p)+\Delta t^{3}\nabla
H\left(\phi_{\Delta t}(q,p)+\Delta t^{3}\theta_{\Delta
t}(q,p)\widetilde{R}_{\Delta t}(q,p)\right)\cdot\widetilde{R}_{\Delta
t}(q,p).$

Using $0\leqslant
1-\min(1,\mathrm{e}^{-x})\leqslant\max(0,x)\leqslant\left\lvert
x\right\rvert$, it holds for $\left(q,p\right)\in\mathcal{A}_{\Delta t}$ with
$p=D(q)^{-1/2}\mathrm{G}$,

$\displaystyle 0\leqslant 1-A_{\Delta t}(q,\mathrm{G})$
$\displaystyle\leqslant\Delta t^{3}\left\lvert\nabla H\left(\phi_{\Delta
t}\left(q,p\right)+\Delta t^{3}\theta_{\Delta t}(q,p)\widetilde{R}_{\Delta
t}\left(q,p\right)\right)\cdot\widetilde{R}_{\Delta
t}\left(q,p\right)\right\rvert,$ (87)

the right-hand side being bounded by a polynomial function of $\left\lVert
q\right\rVert$ and $\left\lVert p\right\rVert$. Similarly, for
$\left(q,D(q)^{-1/2}\mathrm{G}\right)\in\mathcal{A}_{\Delta t}$,

$f(q^{1})-f(q)=\Delta t\nabla f(q)\cdot D(q)^{1/2}\mathrm{G}+\mathrm{O}(\Delta
t^{2}).$

Plugging first the previous estimate and (87) into the second term of (86),
and then using (85) to remove the indicator function, one obtains

$\displaystyle\left\lvert\mathbb{E}_{\mathrm{G}}\left[\left(A_{\Delta
t}\left(q,\mathrm{G}\right)-1\right)\left(f(q^{1})-f(q)\right)\mathbb{1}_{\left(q,D(q)^{-1/2}\mathrm{G}\right)\in\mathcal{B}_{\Delta
t}}\right]\right\rvert$
$\displaystyle\qquad\qquad\leqslant\mathbb{E}_{\mathrm{G}}\left[\left\lvert
A_{\Delta t}\left(q,\mathrm{G}\right)-1\right\rvert\left\lvert
f(q^{1})-f(q)\right\rvert\mathbb{1}_{\left(q,D(q)^{-1/2}\mathrm{G}\right)\in\mathcal{B}_{\Delta
t}}\right]=\mathrm{O}(\Delta t^{4}).$

We therefore conclude that

$(P_{\Delta t}f-f)(q)=\frac{\Delta
t^{2}}{2}\mathcal{L}f(q)+\mathrm{O}\left(\Delta
t^{2(\min(1,r_{3})+1)}\right),$

so that $P_{\Delta t}f-\mathrm{e}^{\Delta
t^{2}\mathcal{L}/2}f=\mathrm{O}\left(\Delta t^{2(\min(1,r_{3})+1)}\right)$ in
view of (40), which leads to (43).

###  6.5 Proof of Proposition 16

Let $x=(q,p)\in\mathcal{X}$. The configurations
$x_{1}=(q_{1},p_{1})\in\mathcal{X}$ and $x_{2}=(q_{2},p_{2})\in\mathcal{X}$
such that $\Phi_{\Delta t}^{\mathrm{GSV}}(x,x_{1},x_{2})=0$ satisfy

$\left\\{\begin{aligned} x_{1}&=x+\frac{\Delta t}{2}J\nabla H(q,p_{1}),\\\
x_{2}&=x_{1}+\frac{\Delta t}{2}J\nabla H(q_{2},p_{1}),\end{aligned}\right.$

where $J$ is defined in (2). Following the same reasoning and using the same
notation as in the proof of Proposition 15, fix $R>0$. For $\Delta
t\left\lVert\nabla
H\right\rVert_{\mathcal{C}^{0}\left(B_{x}^{2R}\right)}\leqslant 2R$, define

$\mathcal{F}_{\Delta
t}^{x}\colon\left\\{\begin{array}[]{ccl}B_{x}^{R}&\to&B_{x}^{R}\\\
x_{1}&\mapsto&x+\dfrac{\Delta t}{2}J\nabla
H\left(q,p_{1}\right),\end{array}\right.$

and for $x_{1}\in B_{x}^{R}$,

$\mathscr{F}_{\Delta
t}^{x_{1}}\colon\left\\{\begin{array}[]{ccl}B_{x_{1}}^{R}&\to&B_{x_{1}}^{R}\\\
x_{2}&\mapsto&x_{1}+\dfrac{\Delta t}{2}J\nabla
H\left(q_{2},p_{1}\right).\end{array}\right.$

The map $\mathcal{F}_{\Delta t}^{x}$ is well-defined since $\frac{\Delta
t}{2}\left\lVert\nabla
H\right\rVert_{\mathcal{C}^{0}(B_{x}^{R})}\leqslant\frac{\Delta
t}{2}\left\lVert\nabla H\right\rVert_{\mathcal{C}^{0}(B_{x}^{2R})}\leqslant
R$. The map $\mathscr{F}_{\Delta t}^{x_{1}}$ is also well-defined since, if
$x_{1}\in B_{x}^{R}$ and $x_{2}\in B_{x_{1}}^{R}$, then $x_{2}\in B_{x}^{2R}$
so that $\frac{\Delta t}{2}\left\lVert\nabla
H\right\rVert_{\mathcal{C}^{0}(B_{x_{1}}^{R})}\leqslant\frac{\Delta
t}{2}\left\lVert\nabla H\right\rVert_{\mathcal{C}^{0}(B_{x}^{2R})}\leqslant
R$. It is then easy to show that, under the condition

$\Delta
t\frac{\left\lVert\nabla^{2}H\right\rVert_{\mathcal{C}^{0}(B_{x}^{2R})}}{2}<1,$

the two maps $\mathcal{F}_{\Delta t}^{x}$ and $\mathscr{F}_{\Delta t}^{x_{1}}$
are contractions. Define $F_{1}\colon\mathcal{X}\to(0,+\infty]$ as

$F_{1}(x)=\left\\{\begin{array}[]{ll}\min\left(\dfrac{2R}{\left\lVert\nabla
H\right\rVert_{\mathcal{C}^{0}\left(B_{x}^{2R}\right)}},\dfrac{2}{\left\lVert\nabla^{2}H\right\rVert_{\mathcal{C}^{0}\left(B_{x}^{2R}\right)}}\right)&\text{if
}\left\lVert\nabla
H\right\rVert_{\mathcal{C}^{0}\left(B_{x}^{2R}\right)}\left\lVert\nabla^{2}H\right\rVert_{\mathcal{C}^{0}\left(B_{x}^{2R}\right)}\neq
0,\\\ +\infty&\text{otherwise},\end{array}\right.$

and $F_{2}\colon\mathcal{X}\to(0,+\infty]$ as

$F_{2}(x)=\min\left(F_{1}(x),\min\limits_{x_{1}\in
B_{x}^{R}}\left\\{F_{1}(x_{1})\right\\}\right).$

Note that $F_{1}$ (respectively $F_{2}$) is continuous on
$\operatorname{dom}F_{1}=\left\\{x\in\mathcal{X},F_{1}(x)<+\infty\right\\}$
(respectively on
$\operatorname{dom}F_{2}=\left\\{x\in\mathcal{X},f_{2}(x)<+\infty\right\\}$).
Note also that $F_{2}(x)$ is positive (and can be infinite). If $\Delta
t<F_{2}(x)$, there exists by the Banach fixed point theorem a unique $x_{1}\in
B_{x}^{R}$ such that $x_{1}-\mathcal{F}_{\Delta t}^{x}(x_{1})=\Phi_{\Delta
t/2}^{\mathrm{Euler\,B}}(x,x_{1})=0$ and $x_{2}\in B_{x_{1}}^{R}\subset
B_{x}^{2R}$ such that $x_{2}-\mathscr{F}_{\Delta
t}^{x_{1}}(x_{2})=\Phi_{\Delta t/2}^{\mathrm{Euler\,A}}(x_{1},x_{2})=0$, so
that $\Phi_{\Delta t}^{\mathrm{GSV}}(x,x_{1},x_{2})=0$. We then denote by
$\chi_{\Delta t}^{\mathrm{GSV},\mathrm{FP}}(x):=(x_{1},x_{2})$. This naturally
defines the numerical flow $\chi_{\Delta t}^{\mathrm{GSV},\mathrm{FP}}$ on the
open set

$\mathcal{A}_{\Delta t}^{\mathrm{GSV},\mathrm{FP}}:=F_{2}^{-1}(\Delta
t,+\infty)\cup\operatorname{Int}F_{2}^{-1}\left\\{+\infty\right\\}.$

For $\Delta t>0$ small enough, this set is nonempty. Indeed, if for all
$x\in\mathcal{X}$, $F_{2}(x)=+\infty$, then $\mathcal{A}_{\Delta
t}^{\mathrm{GSV},\mathrm{FP}}=\mathcal{X}$. Otherwise, there exists
$x_{0}\in\mathcal{X}$ such that $F(x_{0})$ is finite. In this case, we choose
$\Delta t_{\star}=F(x_{0})/2>0$ (which will later be reduced but will still
remain positive), and $\mathcal{A}_{\Delta
t_{\star}}^{\mathrm{GSV},\mathrm{FP}}$ is nonempty. Moreover, for any $\Delta
t\in(0,\Delta t_{\star}]$, it holds $\mathcal{A}_{\Delta
t}^{\mathrm{GSV},\mathrm{FP}}\subset\mathcal{A}_{\Delta
t_{\star}}^{\mathrm{GSV},\mathrm{FP}}$. We next define the map
$G\colon\mathcal{X}\to(0,+\infty]$ by

$G(x)=\min\left(F_{2}(x),\min\limits_{y\in
B_{x}^{2R}}\left\\{F_{2}(S(y))\right\\}\right).$

Note that $G(x)$ is positive (and can be infinite). Assume that $\Delta
t<G(x)$ for some $x\in\mathcal{X}$. Then $x_{2}:=\varphi_{\Delta
t}^{\mathrm{GSV},\mathrm{FP}}(x)$ and $(x_{1},x_{2}):=\chi_{\Delta
t}^{\mathrm{GSV},\mathrm{FP}}(x)$ are well-defined. Moreover, since in this
case $\Delta t<F_{1}(S(x_{2}))$, there exists a unique $x^{\prime}_{1}\in
B_{S(x_{2})}^{R}$ such that $\Phi_{\Delta
t/2}^{\mathrm{Euler\,B}}(S(x_{2}),x^{\prime}_{1})=0$. Since $\Phi_{\Delta
t}^{\mathrm{GSV}}$ is $S$-reversible (see Proposition 4), it holds
$\Phi_{\Delta t}^{\mathrm{Euler\,B}}(S(x_{2}),S(x_{1}))=0$. Using the fact
that $S$ is an isometry, $\left\lVert
S(x_{2})-S(x_{1})\right\rVert=\left\lVert x_{2}-x_{1}\right\rVert\leqslant R$.
Thus, by the uniqueness of the zero of $y\mapsto\Phi_{\Delta
t}^{\mathrm{Euler\,B}}(S(x_{2}),y)$ in $B_{S(x_{2})}^{R}$, one can conclude
that $x^{\prime}_{1}=S(x_{1})$. The same reasoning shows that
$x^{\prime}_{2}=S(x)$. Therefore, it holds

$\chi_{\Delta t}^{\mathrm{GSV},\mathrm{FP}}\circ S\circ\varphi_{\Delta
t}^{\mathrm{GSV},\mathrm{FP}}(x)=\left(S\circ\chi_{\Delta
t,1}^{\mathrm{GSV},\mathrm{FP}}(x),S(x)\right).$

Following the same reasoning as for the construction of $\mathcal{A}_{\Delta
t}^{\mathrm{GSV},\mathrm{FP}}$, there exists $\Delta t_{\star}>0$ such that
for any $\Delta t\in(0,\Delta t_{\star}]$, the map $\chi_{\Delta
t}^{\mathrm{GSV},\mathrm{FP}}$ is $S$-reversible. Then the set
$\mathcal{B}_{\Delta t}^{\mathrm{GSV},\mathrm{FP}}$ defined by (54) is open by
Lemma 6 and nonempty as it contains $G^{-1}(\Delta t,+\infty]$.

To conclude, it remains to prove that $\chi_{\Delta
t}^{\mathrm{GSV},\mathrm{FP}}\colon\mathcal{A}_{\Delta
t}^{\mathrm{GSV},\mathrm{FP}}\to\mathcal{X}^{2}$ is $\mathcal{C}^{1}$ for
$\Delta t\in(0,\Delta t_{\star}]$ and that the second statement of (53) holds.
Note that, for $x\in\mathcal{A}_{\Delta t}^{\mathrm{GSV},\mathrm{FP}}$, one
has $\Phi_{\Delta t}^{\mathrm{GSV}}\left(x,\chi_{\Delta
t}^{\mathrm{GSV},\mathrm{FP}}(x)\right)=0$, and the matrices

$\left\\{\begin{aligned} \nabla_{y}\Phi_{\Delta
t/2}^{\mathrm{Euler\,B}}&=\mathrm{I}_{d}-\frac{\Delta
t}{2}\begin{pmatrix}0_{m}&\nabla_{p}^{2}H\\\
0_{m}&-\nabla_{p,q}^{2}H\end{pmatrix},\\\\[5.69046pt] \nabla_{y}\Phi_{\Delta
t/2}^{\mathrm{Euler\,A}}&=\mathrm{I}_{d}-\frac{\Delta
t}{2}\begin{pmatrix}\nabla_{q,p}^{2}H&0_{m}\\\
-\nabla_{p}^{2}H&0_{m}\end{pmatrix},\end{aligned}\right.$

evaluated respectively at $\left(x,\chi_{\Delta
t,1}^{\mathrm{GSV},\mathrm{FP}}(x)\right)$ and at $\left(\chi_{\Delta
t,1}^{\mathrm{GSV},\mathrm{FP}}(x),\chi_{\Delta
t,2}^{\mathrm{GSV},\mathrm{FP}}(x)\right)$, are invertible since $\Delta
t<2/\left\lVert\nabla^{2}H\right\rVert_{\mathcal{C}^{0}(B_{x}^{2R})}$ so that
the matrix

$\nabla_{y}\Phi_{\Delta
t}^{\mathrm{GSV}}=\begin{pmatrix}\nabla_{y}\Phi_{\Delta
t/2}^{\mathrm{Euler\,B}}&0_{d}\\\\[5.69046pt] \star&\nabla_{y}\Phi_{\Delta
t/2}^{\mathrm{Euler\,A}}\end{pmatrix}$

evaluated at $\left(x,\chi_{\Delta t}^{\mathrm{GSV},\mathrm{FP}}(x)\right)$ is
invertible. Therefore, the map $\chi_{\Delta t}^{\mathrm{GSV},\mathrm{FP}}$ is
$\mathcal{C}^{1}$ by the implicit function theorem, and (53) holds. This
concludes the proof.

###  6.6 Proof of Proposition 18

We follow the strategy of the proof of [26, Proposition 2.5]. We start by
stating and proving two useful lemmas, and introduce to this end some
notation. Since $\mathcal{X}$ is an open subset of $\mathbb{R}^{d}$, there
exists a sequence $(K_{n})_{n\geqslant 1}$ of compact subsets of
$\mathbb{R}^{d}$ such that $\mathcal{X}=\bigcup_{n\geqslant 1}K_{n}$. For
instance, when $\mathcal{X}=\mathcal{O}\times\mathbb{R}^{m}$ with
$\mathcal{O}\subsetneq\mathbb{R}^{m}$, one can choose

$K_{n}=\left\\{x\in\mathcal{X},\,\left\lVert x\right\rVert\leqslant n\text{
and
}\operatorname{dist}(x,\mathbb{R}^{d}\setminus\mathcal{X})\geqslant\frac{1}{n}\right\\}.$

######  Lemma 21.

For any $1\leqslant n,m<+\infty$, the set

$\mathcal{G}_{\Delta t}^{n,m}=\left\\{(x,y)\in K_{n}\times
K_{n}^{k},\,\Phi_{\Delta t}(x,y)=0,\,\nabla_{y}\Phi_{\Delta t}(x,y)\text{ is
invertible and }\left\lVert\nabla_{y}\Phi_{\Delta
t}(x,y)^{-1}\right\rVert\leqslant m\right\\}.$

is compact.

###### Proof.

Since $\mathcal{G}_{\Delta t}^{n,m}\subset K_{n}\times K_{n}^{k}$ and $K_{n}$
is compact, we only have to show that the set is closed. Let
$((x_{j},y_{j}))_{j\geqslant 1}$ be a sequence of elements of
$\mathcal{G}_{\Delta t}^{n,m}$ converging to
$(x,y)\in\mathcal{X}\times\mathcal{X}^{k}$. Since $K_{n}$ is closed, $(x,y)\in
K_{n}\times K_{n}^{k}$. Moreover, $\Phi_{\Delta t}(x,y)=0$ follows from the
continuity of $\Phi_{\Delta t}$. Now, upon extraction, there exists a matrix
$B\in\mathbb{R}^{dk}\times\mathbb{R}^{dk}$ such that

$\left(\nabla_{y}\Phi_{\Delta
t}(x_{j},y_{j})\right)^{-1}\xrightarrow[j\to+\infty]{}B,$

with $\left\lVert B\right\rVert\leqslant m$. For all $j\geqslant 1$, one has

$\nabla_{y}\Phi_{\Delta t}(x_{j},y_{j})\left(\nabla_{y}\Phi_{\Delta
t}(x_{j},y_{j})\right)^{-1}=\mathrm{I}_{dk},$

so that from the $\mathcal{C}^{1}$ regularity of $\Phi_{\Delta t}$, taking the
limit $j\to+\infty$ implies that

$\nabla_{y}\Phi_{\Delta t}(x,y)B=\mathrm{I}_{dk}.$

Thus, the matrix $\nabla_{y}\Phi_{\Delta t}(x,y)$ is invertible, its inverse
being $B$, and it holds $\left\lVert\left(\nabla_{y}\Phi_{\Delta
t}(x,y)\right)^{-1}\right\rVert\leqslant m$. It follows that
$(x,y)\in\mathcal{G}_{\Delta t}^{n,m}$ and $\mathcal{G}_{\Delta t}^{n,m}$ is a
compact set. ∎

Note that from Assumption 5, the set $\mathcal{G}_{\Delta t}^{n,m}$ is
nonempty for $n,m\geqslant 1$ sufficiently large.

######  Lemma 22.

There exists an open subset $\mathcal{D}_{\Delta t}^{n,m}$ such that
$\mathcal{G}_{\Delta t}^{n,m}\subset\mathcal{D}_{\Delta t}^{n,m}$, and a
$\mathcal{C}^{1}$ function $\Lambda_{\Delta t}^{n,m}\colon\mathcal{D}_{\Delta
t}^{n,m}\to\mathcal{X}^{k}$ such that

$\forall(x,y)\in\mathcal{D}_{\Delta t}^{n,m},\qquad\Phi_{\Delta
t}(x,y+\Lambda_{\Delta t}^{n,m}(x,y))=0$

and

$\forall(x,y)\in\mathcal{D}_{\Delta t}^{n,m},\qquad\nabla_{y}\Phi_{\Delta
t}(x,y+\Lambda_{\Delta t}^{n,m}(x,y))\text{ is invertible}.$

Moreover, the set $\mathcal{D}_{\Delta t}^{n,m}$ can be chosen such that there
exists $\widetilde{\alpha}(n,m)>0$ for which,

$\displaystyle\forall(x,y)\in\mathcal{D}_{\Delta t}^{n,m},$
$\displaystyle\,\left\lVert\Lambda_{\Delta
t}^{n,m}(x,y)\right\rVert<\widetilde{\alpha}(n,m)\text{ and}$
$\displaystyle\forall\lambda\in\mathcal{X}^{k}\setminus\left\\{\Lambda_{\Delta
t}^{n,m}(x,y)\right\\},\,\Phi_{\Delta
t}(x,y+\lambda)=0\Longrightarrow\left\lVert\lambda\right\rVert\geqslant\widetilde{\alpha}(n,m).$
(88)

###### Proof.

Fix $n,m\geqslant 1$ and a couple $(x_{0},y_{0})\in\mathcal{G}_{\Delta
t}^{n,m}$. Let
$\mathcal{U}_{0}\subset\mathcal{X}\times\mathcal{X}^{k}\times\mathcal{X}^{k}$
be an open neighborhood of $(x_{0},y_{0},0)$ such that

$\forall(x,y,\lambda)\in\mathcal{U}_{0},\qquad y+\lambda\in\mathcal{X}^{k}.$

This neighborhood exists since $\mathcal{X}$ is an open set. Let us consider
the function

$F\colon\left\\{\begin{aligned} \mathcal{U}_{0}&\to\mathcal{X}^{k}\\\
(x,y,\lambda)&\mapsto\Phi_{\Delta t}(x,y+\lambda).\end{aligned}\right.$

This is a $\mathcal{C}^{1}$ function such that $F(x_{0},y_{0},0)=0$ and the
matrix

$\nabla_{\lambda}F(x_{0},y_{0},0)=\nabla_{y}\Phi_{\Delta t}(x_{0},y_{0})$

is invertible. Therefore, by the implicit function theorem, there exist
$\eta(x_{0},y_{0})>0$, $\alpha(x_{0},y_{0})>0$ and a $\mathcal{C}^{1}$
function $\Lambda_{\Delta
t}^{n,m,0}\colon\mathcal{B}(x_{0},\eta(x_{0},y_{0}))\times\mathcal{B}(y_{0},\eta(x_{0},y_{0}))\to\mathcal{B}(0,\alpha(x_{0},y_{0}))$
such that
$\mathcal{B}(x_{0},\eta(x_{0},y_{0}))\times\mathcal{B}(y_{0},\eta(x_{0},y_{0}))\times\mathcal{B}(0,\alpha(x_{0},y_{0}))\subset\mathcal{U}_{0}$
and

$\forall(x,y,\lambda)\in\mathcal{B}(x_{0},\eta(x_{0},y_{0}))\times\mathcal{B}(y_{0},\eta(x_{0},y_{0}))\times\mathcal{B}(0,\alpha(x_{0},y_{0})),\\\
\Phi_{\Delta t}(x,y+\lambda)=0\iff\lambda=\Lambda_{\Delta t}^{n,m,0}(x,y).$
(89)

By a continuity argument, upon reducing $\eta(x_{0},y_{0})$, one can assume
that the matrix $\nabla_{y}\Phi_{\Delta t}(x,y+\Lambda_{\Delta
t}^{n,m,0}(x,y))$ is invertible for all
$(x,y)\in\mathcal{B}(x_{0},\eta(x_{0},y_{0}))\times\mathcal{B}(y_{0},\eta(x_{0},y_{0}))$,
since this is true for $(x,y)=(x_{0},y_{0})$. Note that, by the definition of
$F$, for any
$(y,\lambda)\in\mathcal{B}(y_{0},\eta(x_{0},y_{0}))\times\mathcal{B}(0,\alpha(x_{0},y_{0}))$,
it holds $y+\lambda\in\mathcal{X}^{k}$ so that $\Phi_{\Delta
t}(x,y+\Lambda_{\Delta t}^{n,m,0}(x,y))$ is well-defined whenever
$(x,y)\in\mathcal{B}(x_{0},\eta(x_{0},y_{0}))\times\mathcal{B}(y_{0},\eta(x_{0},y_{0}))$.

Remark now that $\bigcup_{(x_{0},y_{0})\in\mathcal{G}_{\Delta
t}^{n,m}}\mathcal{B}(x_{0},\eta(x_{0},y_{0}))\times\mathcal{B}(y_{0},\eta(x_{0},y_{0}))$
is an open over of the compact set $\mathcal{G}_{\Delta t}^{n,m}$, from which
it is possible to extract a finite cover:

$\mathcal{G}_{\Delta
t}^{n,m}\subset\bigcup_{i=1}^{\ell^{n,m}}\mathcal{B}(x_{i},\eta_{0}(x_{i},y_{i}))\times\mathcal{B}(y_{i},\eta_{0}(x_{i},y_{i}))=:\mathcal{D}_{\Delta
t}^{n,m},$

where $(x_{i},y_{i})_{1\leqslant i\leqslant\ell^{n,m}}$ are elements of
$\mathcal{G}_{\Delta t}^{n,m}$. Note that the set $\mathcal{D}_{\Delta
t}^{n,m}$ is an open set contained in $\mathcal{X}\times\mathcal{X}^{k}$. Let
us define the function $\Lambda_{\Delta t}^{n,m}\colon\mathcal{D}_{\Delta
t}^{n,m}\to\mathcal{X}^{k}$ as follows: for $(x,y)\in\mathcal{D}_{\Delta
t}^{n,m}$ and $1\leqslant i\leqslant\ell^{n,m}$ for which
$(x,y)\in\mathcal{B}(x_{i},\eta(x_{i},y_{i}))\times\mathcal{B}(y_{i},\eta(x_{i},y_{i}))$,

$\Lambda_{\Delta t}^{n,m}(x,y)=\Lambda_{\Delta t}^{n,m,i}(x,y).$

Let us show that this function is well-defined. If
$(x,y)\in\mathcal{B}(x_{i},\eta(x_{i},y_{i}))\times\mathcal{B}(y_{i},\eta(x_{i},y_{i}))\cap\mathcal{B}(x_{j},\eta(x_{j},y_{j}))\times\mathcal{B}(y_{j},\eta(x_{j},y_{j}))$
with $i\neq j$, we can assume without loss of generality that
$\alpha(x_{i},y_{i})\leqslant\alpha(x_{j},y_{j})$. Then

$\Lambda_{\Delta
t}^{n,m,i}(x,y)\in\mathcal{B}(0,\alpha(x_{i},y_{i}))\subset\mathcal{B}(0,\alpha(x_{j},y_{j}))$

with $\Phi_{\Delta t}\left(x,y+\Lambda_{\Delta t}^{n,m,i}(x,y)\right)=0$. By
(89), this shows that $\Lambda_{\Delta t}^{n,m,i}(x,y)=\Lambda_{\Delta
t}^{n,m,j}(x,y)$.

Let us now define

$\widetilde{\alpha}(n,m)=\min\limits_{1\leqslant
i\leqslant\ell^{n,m}}\alpha(x_{i},y_{i})>0.$ (90)

Upon reducing $\mathcal{D}_{\Delta t}^{n,m}$ to a smaller open set still
containing $\mathcal{G}_{\Delta t}^{n,m}$, one can assume that for all
$(x,y)\in\mathcal{G}_{\Delta t}^{n,m}$, it holds $\Lambda_{\Delta
t}^{n,m}(x,y)\in\mathcal{B}(0,\widetilde{\alpha}(n,m))$ since $\Lambda_{\Delta
t}^{n,m}$ is $\mathcal{C}^{1}$ and $\Lambda_{\Delta t}^{n,m}=0$ on
$\mathcal{G}_{\Delta t}^{n,m}$. The property (88) is then a consequence of
(89). ∎

We are now in position to prove Proposition 18.

###### Proof of Proposition 18.

Let us define the set

$\mathcal{D}_{\Delta t}^{\mathrm{imp}}=\bigcup_{n,m\geqslant
1}\mathcal{D}_{\Delta t}^{n,m}.$

The set $\mathcal{D}_{\Delta t}^{\mathrm{imp}}$ is open and contained in
$\mathcal{X}\times\mathcal{X}^{k}$ as it is the union of open sets contained
in $\mathcal{X}\times\mathcal{X}^{k}$. By construction, the set
$\mathcal{D}_{\Delta t}^{\mathrm{imp}}$ contains $\mathcal{G}_{\Delta t}$
since $\mathcal{G}_{\Delta t}=\bigcup_{n,m\geqslant 1}\mathcal{G}_{\Delta
t}^{n,m}$.

Let us define the function $\Lambda_{\Delta
t}^{\mathrm{imp}}\colon\mathcal{D}_{\Delta
t}^{\mathrm{imp}}\to\mathcal{X}^{k}$ as follows: for
$(x_{0},y_{0})\in\mathcal{D}_{\Delta t}^{\mathrm{imp}}$, there exist
$n_{0},m_{0}\geqslant 1$ such that $(x_{0},y_{0})\in\mathcal{D}_{\Delta
t}^{n_{0},m_{0}}$; then,

$\Lambda_{\Delta t}^{\mathrm{imp}}(x,y)=\Lambda_{\Delta
t}^{n_{0},m_{0}}(x,y).$

The function $\Lambda_{\Delta t}^{\mathrm{imp}}$ is well-defined. Indeed, if
$n,m,n^{\prime},m^{\prime}\geqslant 1$ are such that
$(x,y)\in\mathcal{D}_{\Delta t}^{n,m}\cap\mathcal{D}_{\Delta
t}^{n^{\prime},m^{\prime}}$, we can assume without loss of generality that
$\widetilde{\alpha}(n,m)\leqslant\widetilde{\alpha}(n^{\prime},m^{\prime})$
where $\widetilde{\alpha}(n,m)$ is defined by (90). Then,

$\Lambda_{\Delta
t}^{n,m}(x,y)\in\mathcal{B}(0,\widetilde{\alpha}(n,m))\subset\mathcal{B}(0,\widetilde{\alpha}(n^{\prime},m^{\prime})),$

with $\Phi_{\Delta t}(x,y+\Lambda_{\Delta t}^{n,m}(x,y))=0$, so that (89)
implies that $\Lambda_{\Delta t}^{n,m}(x,y)=\Lambda_{\Delta
t}^{n^{\prime},m^{\prime}}(x,y)$. The last property of Proposition 18 follows
from (88) by considering $\mathcal{V}_{0}=\mathcal{D}_{\Delta
t}^{n_{0},m_{0}}$ and $\alpha_{0}=\widetilde{\alpha}(n_{0},m_{0})$ when
$(x,y)\in\mathcal{D}_{\Delta t}^{n_{0},m_{0}}$. This concludes the proof. ∎

##### Acknowledgements.

The works of T.L., R.S. and G.S. benefit from fundings from the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation program (project EMC2, grant agreement No 810367), and from the
Agence Nationale de la Recherche through the grants ANR-19-CE40-0010-01
(QuAMProcs) and ANR-21-CE40-0006 (SINEQ). This project was initiated as T.L.
was a visiting professor at Imperial College of London (ICL), with a visiting
professorship grant from the Leverhulme Trust. The Department of Mathematics
at ICL and the Leverhulme Trust are warmly thanked for their support.

## References

  * [1] A. Abdulle, G. A. Pavliotis, and G. Vilmart.  Accelerated convergence to equilibrium and reduced asymptotic variance for Langevin dynamics using Stratonovich perturbations.  C. R. Math. Acad. Sci. Paris, 357(4):349–354, 2019.
  * [2] H. C. Andersen.  Rattle: A “velocity” version of the shake algorithm for molecular dynamics calculations.  Journal of Computational Physics, 52(1):24–34, 1983.
  * [3] R. Balian.  From Microphysics to Macrophysics: Methods and Applications of Statistical Physics. Vol. II.  Springer, Berlin, 2007.
  * [4] N. Bou-Rabee and E. Vanden-Eijnden.  Pathwise accuracy and ergodicity of metropolized integrators for SDEs.  Communications on Pure and Applied Mathematics, 63(5):655–696, 2010\.
  * [5] J. A. Brofos and R. R. Lederman.  On numerical considerations for Riemannian Manifold Hamiltonian Monte Carlo.  arXiv preprint, 2111.09995, 2021.
  * [6] E. Cancès, F. Legoll, and G. Stoltz.  Theoretical and numerical comparison of some sampling methods for molecular dynamics.  M2AN Math. Model. Numer. Anal., 41(2):351–389, 2007.
  * [7] R. Douc, E. Moulines, P. Priouret, and P. Soulier.  Markov Chains.  Springer Series in Operations Research and Financial Engineering. Springer, Cham, 2018.
  * [8] S. Duane, A. D. Kennedy, B. J. Pendleton, and D. Roweth.  Hybrid Monte Carlo.  Physics Letters B, 195(2):216–222, 1987.
  * [9] A. B. Duncan, T. Lelièvre, and G. A. Pavliotis.  Variance reduction using nonreversible langevin samplers.  Journal of Statistical Physics, 163(3):457–491, 2016.
  * [10] A. Durmus, E. Moulines, and E Saksman.  Irreducibility and geometric ergodicity of Hamiltonian Monte Carlo.  Ann. Statist., 48(6):3545–3564, 2020.
  * [11] D. J. Evans and G. Morriss.  Statistical Mechanics of Nonequilibrium Liquids.  Cambridge University Press, 2 edition, 2008.
  * [12] Y. Fang, J. M. Sanz-Serna, and R. D. Skeel.  Compressible Generalized Hybrid Monte Carlo.  The Journal of Chemical Physics, 140(17):174108, 2014.
  * [13] M. Fathi and G. Stoltz.  Improving dynamical properties of metropolized discretizations of overdamped Langevin dynamics.  Numerische Mathematik, 136(2):545–602, 2017.
  * [14] D. Frenkel and B. Smit.  Understanding Molecular Simulation.  Academic Press, Inc., USA, 2nd edition, 2001.
  * [15] M. Girolami and B. Calderhead.  Riemann manifold Langevin and Hamiltonian Monte Carlo methods.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73(2):123–214, 2011.
  * [16] E. Hairer, C. Lubich, and G. Wanner.  Geometric Numerical Integration, volume 31 of Springer Series in Computational Mathematics.  Springer-Verlag, Berlin, second edition, 2006.
  * [17] W. K. Hastings.  Monte Carlo sampling methods using Markov chains and their applications.  Biometrika, 57(1):97–109, 1970.
  * [18] A. M. Horowitz.  A generalized guided Monte Carlo algorithm.  Physics Letters B, 268(2):247–252, 1991.
  * [19] J. A. Izaguirre and S. S. Hampton.  Shadow hybrid Monte Carlo: an efficient propagator in phase space of macromolecules.  Journal of Computational Physics, 200(2):581–604, 2004.
  * [20] M. Jardat, O. Bernard, P. Turq, and G. R. Kneller.  Transport coefficients of electrolyte solutions from Smart Brownian dynamics simulations.  The Journal of Chemical Physics, 110(16):7993–7999, 1999.
  * [21] M. Kopec.  Weak backward error analysis for Langevin process.  BIT, 55(4):1057–1103, 2015.
  * [22] B. Leimkuhler, C. Matthews, and G. Stoltz.  The computation of averages from equilibrium and nonequilibrium Langevin molecular dynamics.  IMA J. Numer. Anal., 36(1):13–79, 2016.
  * [23] B. Leimkuhler and S. Reich.  Simulating Hamiltonian Dynamics.  Cambridge Monographs on Applied and Computational Mathematics. Cambridge University Press, Cambridge, 2005.
  * [24] T. Lelièvre, G. Pavliotis, G. Robin, R. Santet, and G. Stoltz, 2022.  in preparation.
  * [25] T. Lelièvre, M. Rousset, and G. Stoltz.  Free Energy Computations: A Mathematical Perspective.  Imperial College Press, 2010.
  * [26] T. Lelièvre, M. Rousset, and G. Stoltz.  Hybrid Monte Carlo methods for sampling probability measures on submanifolds.  Numerische Mathematik, 143(2):379–421, 2019.
  * [27] T. Lelièvre and G. Stoltz.  Partial differential equations and stochastic methods in molecular dynamics.  Acta Numer., 25:681–880, 2016.
  * [28] T. Lelièvre, G. Stoltz, and W. Zhang.  Multiple projection Markov chain Monte Carlo algorithms on submanifolds.  IMA Journal of Numerical Analysis, 2022.
  * [29] S. Livingstone, M. Betancourt, S. Byrne, and M. Girolami.  On the geometric ergodicity of Hamiltonian Monte Carlo.  Bernoulli, 25(4A):3109 – 3138, 2019.
  * [30] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller.  Equation of state calculations by fast computing machines.  The Journal of Chemical Physics, 21(6):1087–1092, 1953.
  * [31] S. P. Meyn and R. L. Tweedie.  Markov Chains and Stochastic Stability.  Cambridge University Press, Cambridge, second edition, 2009.
  * [32] G. N. Milstein and M. V. Tretyakov.  Stochastic Numerics for Mathematical Physics.  Scientific Computation. Springer-Verlag, Berlin, 2004.
  * [33] P. Monmarché.  An entropic approach for Hamiltonian Monte Carlo: the idealized case.  arXiv preprint, 2209.13405, 2022.
  * [34] J. R. Munkres.  Topology.  Prentice Hall, Inc., Upper Saddle River, NJ, 2nd edition, 2000.
  * [35] R. M. Neal.  Probabilistic Inference Using Markov Chain Monte Carlo Methods.  Technical Report CRG-TR-93-1. Department of Computer Science, University of Toronto, 1993.
  * [36] R. Poncet.  Generalized and hybrid Metropolis–Hastings overdamped Langevin algorithms.  arXiv preprint, 1701.05833, 2017.
  * [37] L. Rey-Bellet and K. Spiliopoulos.  Irreversible Langevin samplers and variance reduction: a large deviations approach.  Nonlinearity, 28(7):2081, 2015.
  * [38] L. Rey-Bellet and K. Spiliopoulos.  Variance reduction for irreversible Langevin samplers and diffusion on graphs.  Electronic Communications in Probability, 20:1–16, 2015.
  * [39] L. Rey-Bellet and K. Spiliopoulos.  Improving the convergence of reversible samplers.  J. Stat. Phys., 164(3):472–494, 2016.
  * [40] C. P. Robert.  The Bayesian Choice.  Springer Texts in Statistics. Springer, New York, second edition, 2007\.
  * [41] C. P. Robert and G. Casella.  Markov chains.  In Christian P. Robert and George Casella, editors, Monte Carlo Statistical Methods, Springer Texts in Statistics, pages 205–265. Springer, New York, NY, 2004.
  * [42] G. O. Roberts and J. S. Rosenthal.  Optimal scaling of discrete approximations to Langevin diffusions.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 60(1):255–268, 1998.
  * [43] G. O. Roberts and O. Stramer.  Langevin diffusions and Metropolis-Hastings algorithms.  Methodol. Comput. Appl. Probab., 4:337–357 (2003), 2002.
  * [44] P. J. Rossky, J. D. Doll, and H. L. Friedman.  Brownian dynamics as smart Monte Carlo simulation.  The Journal of Chemical Physics, 69(10):4628–4633, 1978.
  * [45] J-P. Ryckaert, G. Ciccotti, and H. J. C. Berendsen.  Numerical integration of the Cartesian equations of motion of a system with constraints: Molecular dynamics of n-alkanes.  Journal of Computational Physics, 23(3):327–341, 1977.
  * [46] C. Schütte.  Habilitation Thesis.  Freie Universität Berlin, 1999.
  * [47] G. Stoltz and Z. Trstanova.  Langevin dynamics with general kinetic energies.  Multiscale Model. Simul., 16(2):777–806, 2018.
  * [48] D. Talay and L. Tubaro.  Expansion of the global error for numerical schemes solving stochastic differential equations.  Stochastic Anal. Appl., 8(4):483–509 (1991), 1990.
  * [49] M. E. Tuckerman.  Statistical Mechanics: Theory and Molecular Simulation.  Oxford University Press, Oxford, 2010.
  * [50] E. Zappa, M. Holmes-Cerfon, and J. Goodman.  Monte Carlo on manifolds: sampling densities and integrating functions.  Communications on Pure and Applied Mathematics, 71(12):2609–2647, 2018.
  * [51] E. Zeidler.  Applied Functional Analysis, volume 108 of Applied Mathematical Sciences.  Springer-Verlag, New York, 1995.
  * [52] B. J. Zhang, Y. M. Marzouk, and K. Spiliopoulos.  Geometry-informed irreversible perturbations for accelerated convergence of Langevin dynamics.  Statistics and Computing, 32(5):78, September 2022.