11footnotetext:  Department of Mathematics, University of Bergen, Norway.
<EMAIL_ADDRESS>
#  A uniformly accurate scheme for the numerical integration of penalized
Langevin dynamics

Adrien Laurent1

###### Abstract

In molecular dynamics, penalized overdamped Langevin dynamics are used to
model the motion of a set of particles that follow constraints up to a
parameter $\varepsilon$. The most used schemes for simulating these dynamics
are the Euler integrator in $\mathbb{R}^{d}$ and the constrained Euler
integrator. Both have weak order one of accuracy, but work properly only in
specific regimes depending on the size of the parameter $\varepsilon$. We
propose in this paper a new consistent method with an accuracy independent of
$\varepsilon$ for solving penalized dynamics on a manifold of any dimension.
Moreover, this method converges to the constrained Euler scheme when
$\varepsilon$ goes to zero. The numerical experiments confirm the theoretical
findings, in the context of weak convergence and for the invariant measure, on
a torus and on the orthogonal group in high dimension and high codimension.

Keywords:  constrained stochastic differential equations, penalized Langevin
dynamics, manifolds, uniform accuracy, invariant measure.

AMS subject classification (2020):  60H35, 70H45, 37M25.

##  1 Introduction

In molecular dynamics, the overdamped Langevin equation in $\mathbb{R}^{d}$ is
often used for modeling the behavior of a large set of particles in a high
friction regime. It is given by

$dX(t)=f(X(t))dt+\sigma dW(t),\quad X(0)=X_{0},$ (1.1)

where $f$ is a smooth Lipschitz function (typically of the form $f=-\nabla V$
for $V$ a smooth potential), $\sigma>0$ is a constant scalar, and $W$ is a
standard $d$-dimensional Brownian motion in $\mathbb{R}^{d}$ on a probability
space equipped with a filtration $(\mathcal{F}_{t})$ and fulfilling the usual
assumptions. If the particles are subject to smooth constraints
$\zeta\colon\mathbb{R}^{d}\rightarrow\mathbb{R}^{q}$, such as strong covalent
bonds between atoms or fixed angles in molecules, the dynamics follow the
constrained overdamped Langevin equation

$dX^{0}(t)=\Pi_{\mathcal{M}}(X^{0}(t))f(X^{0}(t))dt+\sigma\Pi_{\mathcal{M}}(X^{0}(t))\circ
dW(t),\quad X^{0}(0)=X_{0}\in\mathcal{M},$ (1.2)

where the solution lies on the manifold
$\mathcal{M}=\\{x\in\mathbb{R}^{d},\zeta(x)=0\\}$ with codimension $q$ thanks
to $\Pi_{\mathcal{M}}\colon\mathbb{R}^{d}\rightarrow\mathbb{R}^{d\times d}$,
the orthogonal projection on the tangent bundle of the manifold $\mathcal{M}$.

In physical applications, constrained systems are often used as a limit model
for stiff equations. For instance, in the dynamics of a diatomic molecule, the
distance between the two atoms oscillates around an average length, called the
bond length (see, for instance, [34, Sect. 1.2.1] on the interactions of
particles). One can work with a simpler constrained dynamics where the
distance between the atoms is fixed as a constraint, or with the original
(possibly stiff) dynamics in $\mathbb{R}^{d}$. We refer the reader to [41,
34], and references therein, for discussions on the use of constraints and
penalizations in molecular dynamics. For overdamped Langevin dynamics (1.1),
choosing a function $\widetilde{f}=-\nabla\widetilde{V}$ with the potential

$\widetilde{V}=V-\frac{\sigma^{2}}{4}\ln(\det(G))+\frac{1}{2\varepsilon}\left|\zeta\right|^{2}$

gives penalized Langevin dynamics [16] of the form

$dY^{\varepsilon}(t)=f(Y^{\varepsilon}(t))dt+\sigma
dW(t)+\frac{\sigma^{2}}{4}\nabla\ln(\det(G))(Y^{\varepsilon}(t))dt-\frac{1}{\varepsilon}(g\zeta)(Y^{\varepsilon}(t))dt,$
(1.3)

where we fix $Y^{\varepsilon}(0)=X_{0}$, the parameter $\varepsilon>0$ is
fixed with arbitrary size, $f=-\nabla V$,
$g=\nabla\zeta\colon\mathbb{R}^{d}\rightarrow\mathbb{R}^{d\times q}$, and
$G=g^{T}g\colon\mathbb{R}^{d}\rightarrow\mathbb{R}^{q\times q}$ is the Gram
matrix. It was shown in [16, Appx. C] that the solution $Y^{\varepsilon}$ of
(1.3) converges strongly to the solution $X^{0}$ of the constrained dynamics
(1.2) if $X_{0}\in\mathcal{M}$. The additional term
$\frac{\sigma^{2}}{4}\nabla\ln(\det(G))$ is a correction term (called the
Fixman correction) that is needed to obtain the convergence to the constrained
dynamics (1.2) (see [34, Sect. 3.2.3.4] and references therein). Thus, for
$\varepsilon$ small, the trajectory of the solution of (1.3) lies in the
vicinity of the manifold $\mathcal{M}$. This penalization can also appear
naturally when simulating Langevin dynamics with a stiff potential (see, for
instance, [45, Sect. 5.1]). One is then interested in numerical schemes that
are robust with respect to the parameter $\varepsilon$ and that lie on the
manifold $\mathcal{M}$ in the limit $\varepsilon\rightarrow 0$. In this paper,
we study the following similar penalized dynamics in $\mathbb{R}^{d}$ to
simulate trajectories in a vicinity of the manifold $\mathcal{M}$:

$\displaystyle dX^{\varepsilon}(t)$
$\displaystyle=f(X^{\varepsilon}(t))dt+\sigma
dW(t)+\frac{\sigma^{2}}{4}\nabla\ln(\det(G))(X^{\varepsilon}(t))dt-\frac{1}{\varepsilon}(gG^{-1}\zeta)(X^{\varepsilon}(t))dt,$
(1.4)

where $X^{\varepsilon}(0)=X_{0}$. It is a simpler version of (1.3) that also
evolves in a vicinity of the manifold $\mathcal{M}$ in the limit
$\varepsilon\rightarrow 0$. One result of this paper is the strong convergence
of the solution $X^{\varepsilon}$ of (1.4) to the solution $X^{0}$ of (1.2) if
$X_{0}\in\mathcal{M}$. We mention that in the deterministic setting, that is,
when $\sigma=0$, equation (1.4) is a singular perturbation problem, and it
converges to a differential algebraic equation (DAE) of index two in the limit
$\varepsilon\rightarrow 0$ (see [22, Chaps. VI-VII]). We propose in this
article a method that is robust with respect to the parameter $\varepsilon$
for solving equations of the form (1.4), and we leave the creation of robust
integrators for solving (1.3) for future work for the sake of clarity.

There are different ways to approximate the solution of the dynamics (1.4). A
strong approximation focuses on approximating the realization of a single
trajectory of (1.4) for a given realization of the Wiener process $W$. A weak
approximation approximates the average of functionals of the solution at a
fixed time $T$, that is, quantities of the form
$\mathbb{E}[\phi(X^{\varepsilon}(T))]$ for $\phi$ a smooth test function. In
addition, under growth and smoothness assumptions on the vector fields in
(1.4) (see, for instance, [23]), the dynamics (1.4) naturally satisfy an
ergodicity property; that is, there exists a unique invariant measure
$d\mu_{\infty}^{\varepsilon}$ in $\mathbb{R}^{d}$ that has a density
$\rho_{\infty}^{\varepsilon}$ with respect to the Lebesgue measure, such that
for all test functions $\phi$,

$\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}\phi(X(t))dt=\int_{\mathbb{R}^{d}}\phi(x)d\mu_{\infty}^{\varepsilon}(x)\quad\text{almost
surely}.$

An approximation for the invariant measure focuses on approximating the
average of a functional in the stationary state, that is, the quantity
$\int_{\mathbb{R}^{d}}\phi(x)d\mu_{\infty}^{\varepsilon}(x)$. This is a
computational challenge when the dimension $d$ is high, which is the case in
the context of molecular dynamics where the dimension is proportional to the
number of particles, as a standard quadrature formula becomes prohibitively
expensive in high dimension. We emphasize that the invariant measure
$\mu_{\infty}^{\varepsilon}$ becomes singular with respect to the Lebesgue
measure on $\mathbb{R}^{d}$ in the limit $\varepsilon\rightarrow 0$, and tends
weakly as $\varepsilon\rightarrow 0$ to $d\mu_{\infty}^{0}$, a measure that is
absolutely continuous to $d\sigma_{\mathcal{M}}$, the canonical measure on
$\mathcal{M}$ induced by the Euclidean metric of $\mathbb{R}^{d}$. In this
paper, we propose weak convergence results for a new uniformly accurate
integrator for solving (1.4), and numerical experiments in the weak context
and for the invariant measure, as we recall that a scheme of weak order $r$
automatically has order $p\geq r$ for the invariant measure (see, for
instance, [37]).

The most used discretization for solving (1.4) is the explicit Euler
integrator in $\mathbb{R}^{d}$ (see [15, 33, 34, 35], for instance),

$X_{n+1}=X_{n}+\sqrt{h}\sigma\xi+hf(X_{n})+h\frac{\sigma^{2}}{4}\nabla\ln(\det(G))(X_{n})-\frac{h}{\varepsilon}(gG^{-1}\zeta)(X_{n}).$
(1.5)

This integrator has weak order one of accuracy, but it faces some severe
stepsize restriction due to its instability, typically of the form
$h\ll\varepsilon$, in order to be accurate in the regime
$\varepsilon\rightarrow 0$. Since the solution $X^{\varepsilon}(t)$ of (1.4)
converges to the solution $X^{0}(t)$ of (1.2) when $\varepsilon\rightarrow 0$,
one can use integrators for the limit equation (1.2) and apply them to solve
the original problem (1.4) when $\varepsilon$ is close to zero. Indeed, one
can prove that the solution $X^{\varepsilon}(t)$ of (1.4) stays at distance
$\mathcal{O}(\sqrt{\varepsilon})$ of the constrained solution $X^{0}(t)$ of
(1.2) (see Theorem 2.3). Thus, if the timestep of the integrator is small
enough and satisfies $\varepsilon\ll h$, then this integrator is consistent
for solving (1.4). The alternative for the discretization of (1.2) on the
manifold of the explicit Euler scheme (1.5) is the constrained Euler scheme

$X^{0}_{n+1}=X^{0}_{n}+\sqrt{h}\sigma\xi+hf(X^{0}_{n})+g(X^{0}_{n})\lambda^{0}_{n+1},\quad\zeta(X^{0}_{n+1})=0,$
(1.6)

where $\lambda^{0}_{n+1}\in\mathbb{R}^{q}$ acts as a Lagrange multiplier, and
is entirely determined by the constraint $\zeta(X^{0}_{n+1})=0$. This
integrator has weak order one for solving (1.2) (see [34, Sect. 3.2.4]) and it
lies on the manifold $\mathcal{M}$. It is a consistent approximation of (1.4)
if $\varepsilon$ is close to zero and $\varepsilon\ll h$. This integrator is,
however, not appropriate for solving (1.4) if the size of $\varepsilon$ is of
the order of one, since the exact solution does not evolve in a neighborhood
of the manifold in this regime. We mention a few other techniques to integrate
numerically (1.4) or (1.2). In [28, 30], high order Runge-Kutta methods are
proposed for sampling the invariant measure in $\mathbb{R}^{d}$ and on
manifolds. The paper [36] presents a constrained integrator based on the
RATTLE scheme (see [42, 5, 22]) in the context of the underdamped Langevin
dynamics. Some of the previously cited discretizations can be combined with
Metropolis-Hastings rejection procedures [39, 24], as done, for instance, in
[20, 8, 35, 46, 36]. As for the Euler integrators (1.5)-(1.6), all of the
previously mentioned methods are consistent provided that the parameter
$\varepsilon$ and the timestep satisfy $h\ll\varepsilon$ in the context of
methods in $\mathbb{R}^{d}$, and satisfy $\varepsilon\ll h$ in the context of
methods on the manifold $\mathcal{M}$. When applied in the regime where
$\varepsilon$ and $h$ share the same order of magnitude, the accuracy of these
methods quickly deteriorates, and they may face stability issues.

In past decades, different solutions were proposed for treating the loss of
accuracy in the intermediate regime $h\sim\varepsilon$ in the context of
multiscale problems with the help of uniformly accurate (UA) methods. These
methods are capable of solving dynamics indexed by a possibly stiff parameter
$\varepsilon$ with an accuracy and a cost both independent of $\varepsilon$. A
uniformly accurate method is automatically asymptotic-preserving (AP), that
is, it converges in the two regimes $\varepsilon\rightarrow 0$ and
$\varepsilon\sim 1$, but the converse is not true in general. We refer the
reader to the review [25] and references therein for examples of AP
integrators for solving multiscale problems. We mention in particular the
paper [41], which proposes a penalized Hamiltonian dynamics, and AP
discretizations for solving it, and the paper [7] that gives an AP scheme for
the approximation of a class of multiscale SDEs. In [17, 44, 29],
trigonometric and multirevolution integrators are considered for solving
highly oscillatory SDEs (see also the deterministic works [38, 9, 10, 14]). We
mention the recent papers [11, 12, 13, 4] (see also the references therein)
that introduce uniformly accurate methods for solving a variety of multiscale
problems. There is a rich literature on AP and UA methods, but, to the best of
our knowledge, the problem we study here and the techniques we consider are
new. We propose in this paper a new consistent integrator with uniform
accuracy and uniform cost for solving penalized Langevin dynamics, that is, a
method for solving (1.4) whose accuracy and cost do not depend on the
parameter $\varepsilon$.

The article is organized as follows. Section 2 is devoted to the presentation
of the new integrator and to the main convergence results. In Section 3, we
build a weak asymptotic expansion of the solution of (1.4) that is uniform in
$\varepsilon$, and we use it for proving the uniform accuracy of our
integrator. We compare in Section 4 the new integrator with the explicit Euler
scheme in $\mathbb{R}^{d}$ (1.5) and the constrained Euler scheme (1.6) on
$\mathcal{M}$ in numerical experiments on a torus and on the orthogonal group
to confirm its order of convergence in the weak context and for sampling the
invariant measure. Finally, we present some possible extensions and future
work in Section 5.

##  2 Uniformly accurate integrator for penalized Langevin dynamics

In this section, we present the new uniformly accurate integrator and the main
convergence results of this paper. The proofs are postponed to Section 3. Let
us first lay down a few notations and assumptions. We assume in the rest of
the article that $\mathcal{M}=\\{x\in\mathbb{R}^{d},\zeta(x)=0\\}$ is a
compact and smooth manifold of codimension $q\geq 1$ embedded in
$\mathbb{R}^{d}$, where the constraints are given by the smooth map
$\zeta\colon\mathbb{R}^{d}\rightarrow\mathbb{R}^{q}$. We write
$g=\nabla\zeta\colon\mathbb{R}^{d}\rightarrow\mathbb{R}^{d\times q}$ and we
assume that the Gram matrix $G(x)=g^{T}(x)g(x)\in\mathbb{R}^{q\times q}$ is
invertible for all $x$ in $\mathcal{M}$. With these notations, the projection
$\Pi_{\mathcal{M}}$ on the tangent bundle is given by
$\Pi_{\mathcal{M}}(x)=I_{d}-G^{-1}(x)g(x)g^{T}(x)$. The test functions
typically belong to a subspace of
$\mathcal{C}^{\infty}_{P}(\mathbb{R}^{d},\mathbb{R})$, the vector space of
$\mathcal{C}^{\infty}$ functions $\phi(x)$ such that all partial derivatives
up to any order have a polynomial growth of the form

$\big{|}\phi^{(k)}(x)\big{|}\leq C(1+\left|x\right|^{K}),$

where the constants $C$ and $K$ are independent of $x\in\mathbb{R}^{d}$ (but
can depend on $k$), and where we denote by $\left|x\right|=(x^{T}x)^{1/2}$ the
Euclidean norm in $\mathbb{R}^{d}$. Similarly, we denote by
$\mathcal{C}^{p}_{P}(\mathbb{R}^{d},\mathbb{R})$ the space of
$\mathcal{C}^{p}$ functions whose partial derivatives up to order $p$ have
polynomial growth. Letting
$\varphi\colon\mathbb{R}^{d}\rightarrow\mathbb{R}^{d\times k}$, we use the
following notations for differentials: for all vectors $x$, $a^{1}$, …,
$a^{m}\in\mathbb{R}^{d}$, we denote

$\varphi^{(m)}(x)(a^{1},\dots,a^{m})=\sum_{i_{1},\dots,i_{m}=1}^{d}\frac{\partial^{m}\varphi}{\partial
x_{i_{1}}\dots\partial x_{i_{m}}}(x)\thinspace a^{1}_{i_{1}}\dots
a^{m}_{i_{m}}.$

For the sake of clarity, we will often drop the coefficient $x$, and if $m=1$,
we also use the notation $\varphi^{\prime}$ for the Jacobian matrix of
$\varphi$. Moreover, for $(e_{i})$ the canonical basis of $\mathbb{R}^{d}$, we
write

$\Delta\varphi(x)=\sum_{i=1}^{d}\varphi^{\prime\prime}(x)(e_{i},e_{i})\quad\text{and}\quad(\operatorname{div}\varphi(x))_{j}=\sum_{i=1}^{d}\varphi_{ij}^{\prime}(x)(e_{i}).$

In the rest of the paper, we make the following assumption in the spirit of
the regularity assumptions made in [16, Appx. C] and [40, Chap. 2].

######  Assumption 2.1.

The map $f$ is bounded, is Lipschitz and lies in $\mathcal{C}^{3}_{P}$. The
maps $g$ and $g^{\prime}$ are bounded in $\mathbb{R}^{d}$, and there exists
$c>0$ such that, for $x$, $y\in\mathbb{R}^{d}$,

$\left|G_{y}(x)\right|\geq c(1+\left|y\right|)^{-1},\quad\text{where}\quad
G_{y}(x)=\int_{0}^{1}g^{T}(x+\tau y)d\tau g(x).$ (2.1)

In addition, there exists a smooth change of coordinate

$\displaystyle\psi\colon$
$\displaystyle\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ $\displaystyle
x\mapsto\begin{pmatrix}\varphi(x)\\\ \zeta(x)\end{pmatrix}$

where $\varphi\colon\mathbb{R}^{d}\rightarrow\mathbb{R}^{d-q}$ satisfies
$\varphi^{\prime}g=0$. The map $\psi$ lies in $\mathcal{C}^{5}_{P}$ and is
invertible, $\psi^{\prime}$ and $\psi^{\prime\prime}$ are Lipschitz and there
exist two constants $c$, $C>0$ such that
$c\leq\left|\psi^{\prime}(x)\right|\leq C$.

######  Remark 2.2.

Assumption 2.1 is almost the same as the one given in [16, Appx. C] to prove
the strong convergence of the dynamics (1.3) to the constrained dynamics
(1.2). The difference lies in the additional estimate (2.1) that replaces the
weaker assumption that the Gram matrix $G(x)=G_{0}(x)$ is invertible on
$\mathcal{M}$. We use the estimate (2.1) for obtaining a uniform expansion of
the Lagrange multipliers in the new method and for proving that the new method
evolves in a neighborhood of the manifold (see Lemma 3.8). Note that the
existence of the change of coordinate $\psi$ is always valid in a neighborhood
of the smooth manifold $\mathcal{M}$. The same goes for the estimate (2.1) for
$x$ in a neighborhood of $\mathcal{M}$ and $y$ in a ball centered on zero.
Assumption 2.1 is valid in particular if $\mathcal{M}$ is a vector subspace of
$\mathbb{R}^{d}$, but it is quite restrictive. It would be interesting to
extend the results of this paper under simpler regularity assumptions made
only on the manifold, as numerical experiments hint that the results presented
in this paper still stand without global assumptions. This is matter for
future work.

Under Assumption 2.1, the problems (1.2) and (1.4) are well posed, and we
obtain the strong convergence of the penalized dynamics (1.4) to the
constrained dynamics (1.2).

######  Theorem 2.3.

Under Assumption 2.1, the solution $X^{\varepsilon}$ of the penalized dynamics
(1.4) converges strongly to $X^{0}$, the solution of the constrained dynamics
(1.2); that is, for all $t\leq T$, there exists a constant $C>0$ such that,
for all $\varepsilon>0$,

$\sup_{t\leq
T}\mathbb{E}\Big{[}\left|X^{\varepsilon}(t)-X^{0}(t)\right|^{2}\Big{]}\leq
C\varepsilon.$

Moreover, $\zeta(X^{\varepsilon}(t))$ satisfies

$\sup_{t\leq
T}\mathbb{E}\Big{[}\left|\zeta(X^{\varepsilon}(t))\right|^{2}\Big{]}\leq
C\varepsilon.$

This result was first introduced in [16, Appx. C] for slightly different
penalized dynamics. The proof is almost identical, but we present it in Appx.
A for the sake of completeness. In the deterministic setting with $\sigma=0$,
the convergence to the manifold is of order $1$ in $\varepsilon$ instead of
order $1/2$. We introduce the following additional assumption.

######  Assumption 2.4.

There exists $c>0$ such that, for $x$, $y\in\mathbb{R}^{d}$,

$\left|G_{y}(x)\right|\geq c,\quad\text{where}\quad
G_{y}(x)=\int_{0}^{1}g^{T}(x+\tau y)d\tau g(x).$

Assumption (2.4) is a stronger version of the inequality (2.1) and is in the
spirit of the concept of admissible Lagrange multipliers [36]. It is always
satisfied for $x$ in a neighborhood of the manifold $\mathcal{M}$ and $y$ in a
ball centered on zero if we assume that the Gram matrix $G(x)=G_{0}(x)$ is
invertible on $\mathcal{M}$. We emphasize that we do not need this assumption
for proving the uniform accuracy property of the new method, but we use it for
obtaining uniform estimates in the regime $\varepsilon\rightarrow 0$ and on
the numerical implementation of the uniformly accurate method.

We introduce the new integrator for approximating the penalized dynamics (1.4)
with cost and accuracy independent of the parameter $\varepsilon$, and a cost
comparable to that of the constrained Euler scheme (1.6) in terms of the
number of evaluations of the functions $f$, $\zeta$, $g$, and $g^{\prime}$.

New Method (Uniform discretization of penalized overdamped Langevin dynamics)

$X_{0}^{\varepsilon}=X_{0}$ $\in\mathcal{M}$

for $n\geq 0$ do

$\displaystyle X_{n+1}^{\varepsilon}$
$\displaystyle=X_{n}^{\varepsilon}+\sqrt{h}\sigma\xi_{n}+hf(X_{n}^{\varepsilon})+\frac{(1-e^{-h/\varepsilon})^{2}}{2}(g^{\prime}(gG^{-1}\zeta)G^{-1}\zeta)(X_{n}^{\varepsilon})$
$\displaystyle+\frac{\sigma^{2}\varepsilon}{8}(1-e^{-2h/\varepsilon})\nabla\ln(\det(G))(X_{n}^{\varepsilon})+g(X_{n}^{\varepsilon})\lambda^{\varepsilon}_{n+1},$
(2.2) $\displaystyle\zeta(X_{n+1}^{\varepsilon})$
$\displaystyle=e^{-h/\varepsilon}\zeta(X_{n}^{\varepsilon})+\sigma\sqrt{\frac{\varepsilon}{2}(1-e^{-2h/\varepsilon})}g^{T}(X_{n}^{\varepsilon})\xi_{n}$
$\displaystyle+\varepsilon(1-e^{-h/\varepsilon})\big{(}g^{T}f+\frac{\sigma^{2}}{4}g^{T}\nabla\ln(\det(G))+\frac{\sigma^{2}}{2}\operatorname{div}(g)\big{)}(X_{n}^{\varepsilon})$
$\displaystyle+\sigma^{2}\Big{(}\varepsilon(1-e^{-h/\varepsilon})-\sqrt{\frac{\varepsilon
h}{2}(1-e^{-2h/\varepsilon})}\Big{)}$
$\displaystyle\times\Big{(}\sum_{i=1}^{d}(g^{\prime}(gG^{-1}g^{T}e_{i}))^{T}gG^{-1}g^{T}e_{i}-\sum_{i=1}^{d}(g^{\prime}(e_{i}))^{T}gG^{-1}g^{T}e_{i}\Big{)}(X_{n}^{\varepsilon}).$

end for

The new method works in a way similar to the constrained Euler integrator
(1.6). Knowing the approximation $X_{n}^{\varepsilon}$ of
$X^{\varepsilon}(nh)$, we project a modified Euler step on a modified manifold
defined by the constraint given in (2.2) (in place of
$\zeta(X_{n+1}^{\varepsilon})=0$ for the constrained Euler integrator (1.6)).
We project the modified step in the direction $g(X_{n}^{\varepsilon})$ with
the help of a Lagrange multiplier $\lambda^{\varepsilon}_{n+1}$. For the
implementation of the method, one can use, for instance, a fixed point
iteration or a Newton method at each step to find the solution
$(X_{n+1}^{\varepsilon},\lambda^{\varepsilon}_{n+1})$ of the implicit system
of equations (2.2).

In order for the discretization (2.2) to be well-defined, we use bounded
random variables. The $\xi_{n}$ are independent and bounded discrete random
vectors that have the same moments as standard Gaussian random vectors up to
order four, in the spirit of [40, Chap. 2], that is, for instance, that their
components satisfy

$\mathbb{P}(\xi_{i}=0)=\frac{2}{3}\quad\text{and}\quad\mathbb{P}(\xi_{i}=\pm\sqrt{3})=\frac{1}{6},\quad
i=1,\dots,d.$ (2.3)

Note that using truncated Gaussian random variables would also work.

######  Remark 2.5.

The new integrator is related to the popular idea of backward error analysis
and modified equations for SDEs (see, for instance, [48, 1, 18, 26, 27]). The
idea is to define a projection method (see [21, Sect. IV.4]) with a modified
constraint in place of $\zeta(X_{n})=0$. Instead of evaluating the stiff term
$\frac{h}{\varepsilon}gG^{-1}\zeta$ as in the Euler scheme (1.5), we project a
modified step of the explicit Euler scheme in $\mathbb{R}^{d}$ on a manifold
that is close to $\mathcal{M}$ when $\varepsilon\ll 1$ and whose constraint is
given by a truncation of a uniform expansion of $\zeta(X^{\varepsilon})$. When
$\varepsilon\rightarrow\infty$, the expression of the constraint
$\zeta(X_{n+1}^{\varepsilon})$ in (2.2) tends to a truncated Taylor expansion
in $h$ around $X_{n}^{\varepsilon}$, while for $\varepsilon\rightarrow 0$,
$\zeta(X_{n+1}^{\varepsilon})$ tends to zero, which enforces that the
integrator lies on $\mathcal{M}$. These intuitions will be made rigorous in
Section 3.

######  Remark 2.6.

In the context of a manifold $\mathcal{M}$ of codimension $q=1$, the Gram
matrix $G(x)$ and $g(x)^{T}e_{i}=g_{i}(x)$ are real numbers, so that
$\nabla\ln(\det(G))=2G^{-1}g^{\prime}(g)$ and

$\sum_{i=1}^{d}(g^{\prime}(gG^{-1}g^{T}e_{i}))^{T}gG^{-1}g^{T}e_{i}=G^{-1}(g^{\prime}(g))^{T}g=\sum_{i=1}^{d}(g^{\prime}(e_{i}))^{T}gG^{-1}g^{T}e_{i}.$

The discretization (2.2) thus reduces to

$\displaystyle X_{n+1}^{\varepsilon}$
$\displaystyle=X_{n}^{\varepsilon}+\sqrt{h}\sigma\xi_{n}+hf(X_{n}^{\varepsilon})+\frac{(1-e^{-h/\varepsilon})^{2}}{2}(\zeta^{2}G^{-2}g^{\prime}(g))(X_{n}^{\varepsilon})$
$\displaystyle+\frac{\sigma^{2}\varepsilon}{4}(1-e^{-2h/\varepsilon})(G^{-1}g^{\prime}(g))(X_{n}^{\varepsilon})+g(X_{n}^{\varepsilon})\lambda^{\varepsilon}_{n+1},$
(2.4) $\displaystyle\zeta(X_{n+1}^{\varepsilon})$
$\displaystyle=e^{-h/\varepsilon}\zeta(X_{n}^{\varepsilon})+\sigma\sqrt{\frac{\varepsilon}{2}(1-e^{-2h/\varepsilon})}g^{T}(X_{n}^{\varepsilon})\xi_{n}$
$\displaystyle+\varepsilon(1-e^{-h/\varepsilon})(g^{T}f+\frac{\sigma^{2}}{2}G^{-1}g^{T}g^{\prime}(g)+\frac{\sigma^{2}}{2}\operatorname{div}(g))(X_{n}^{\varepsilon}).$

We present in the rest of the section the uniform accuracy property of the
discretization (2.2), and we show that the integrator converges to the
constrained Euler scheme (1.6) when $\varepsilon\rightarrow 0$. The different
convergence results are summarized by the following commutative diagram, where
$T=Nh$ is fixed. Note that, as we present a convergence result in $h$ that is
uniform in $\varepsilon$, the two arrows for the convergence in $h$ rely on
the same Theorem 2.7.

${\text{integrator }X^{\varepsilon}_{N}\text{ in
}\mathbb{R}^{d}}$${\text{integrator }X^{0}_{N}\text{ on
}\mathcal{M}}$${\text{solution }X^{\varepsilon}(T)\text{
of\leavevmode\nobreak\ \eqref{equation:Langevin_epsilon}}}$${\text{solution
}X^{0}(T)\text{ of\leavevmode\nobreak\
\eqref{equation:projected_Langevin}}}$$\scriptstyle{h\rightarrow
0}$$\scriptstyle{(Thm.\thinspace\leavevmode\nobreak\
\ref{theorem:uniform_consistency_UA_scheme})}$$\scriptstyle{\varepsilon\rightarrow
0}$$\scriptstyle{(Thm.\thinspace\leavevmode\nobreak\
\ref{theorem:convergence_integrator_epsilon})}$$\scriptstyle{h\rightarrow
0}$$\scriptstyle{(Thm.\thinspace\leavevmode\nobreak\
\ref{theorem:uniform_consistency_UA_scheme})}$$\scriptstyle{\varepsilon\rightarrow
0}$$\scriptstyle{(Thm.\thinspace\leavevmode\nobreak\
\ref{theorem:CV_penalized_constrained})}$

We now state the main result of this work, that is, the uniform accuracy of
the discretization (2.2) for approximating the solution of the penalized
Langevin dynamics (1.4).

######  Theorem 2.7.

Under Assumption 2.1, the integrator $(X_{n}^{\varepsilon})$ given by (2.2) is
a consistent uniformly accurate approximation of the solution
$X^{\varepsilon}(t)$ of the penalized Langevin dynamics (1.4); that is, for a
given test function $\phi\in\mathcal{C}^{5}_{P}$, there exist $h_{0}>0$, $C>0$
such that for all $\varepsilon>0$, $h\leq h_{0}$, the following estimate
holds:

$\left|\mathbb{E}[\phi(X_{n}^{\varepsilon})]-\mathbb{E}[\phi(X^{\varepsilon}(nh))]\right|\leq
C\sqrt{h},\quad n=0,1,\dots,N,\quad Nh=T.$ (2.5)

######  Remark 2.8.

Note that Theorem 2.7 states the uniform consistency, but not the uniform weak
order one, as one could expect. The discretization (2.2) has weak order one if
$\varepsilon=\varepsilon_{0}$ is fixed or in the limit $\varepsilon\rightarrow
0$, but an order reduction occurs in the intermediate regime
$0<\varepsilon<\varepsilon_{0}$, and the integrator only has weak order $1/2$
with respect to $h$ in general. For the sake of simplicity, we leave the
creation of uniformly accurate integrators of higher weak order for future
works.

We present the proof of Theorem 2.7 in Section 3. It relies on a weak
expansion in $h$ of the solution of (1.4) that is uniform in $\varepsilon$.
One could directly use this uniform expansion as an explicit numerical
integrator for solving (1.4). It would also yield a uniformly accurate scheme
and would not require one to solve a fixed point problem. However, in the
limit $\varepsilon\rightarrow 0$, this integrator would almost surely not stay
on the manifold. The crucial geometric property that the integrator lies on
the manifold when $\varepsilon\rightarrow 0$ is satisfied for the new method,
as stated in the following result.

######  Theorem 2.9.

Under Assumption 2.1, the integrator $(X_{n}^{\varepsilon})$ in (2.2)
converges to the Euler scheme on the manifold (1.6) when
$\varepsilon\rightarrow 0$; that is, for $h$ and $N$ fixed such that $T=Nh$,
there exists a constant $C_{h}>0$ that depends on $h$ but not on $\varepsilon$
such that

$\left|\zeta(X^{\varepsilon}_{n})\right|\leq C_{h}\sqrt{\varepsilon},\quad
n=0,1,\dots,N.$

In addition, if Assumption 2.4 is satisfied, then, for $h_{0}$ small enough,
for $h\leq h_{0}$ fixed, there exists a constant $C_{h}>0$ such that

$\left|X_{n}^{\varepsilon}-X_{n}^{0}\right|\leq C_{h}\sqrt{\varepsilon},\quad
n=0,1,\dots,N,\quad Nh=T.$ (2.6)

######  Remark 2.10.

In the deterministic context (i.e. when $\sigma=0$), the uniform accuracy of
the discretization (2.2) still holds, and the speed of convergence to the
manifold $\mathcal{M}$ of both the exact solution and the integrator are in
$\mathcal{O}(\varepsilon)$. To the best of our knowledge, the integrator given
by (2.2) is the first integrator with the uniform accuracy property for
solving the singular perturbation problem (1.4) with $\sigma=0$. However,
similar expansions that are uniform with respect to $\varepsilon$ are
presented in [22, Chaps. VI-VII] and references therein.

######  Remark 2.11.

Another widely used scheme on the manifold is the Euler scheme with implicit
projection direction,

$X^{0}_{n+1}=X^{0}_{n}+hf(X^{0}_{n})+\sqrt{h}\sigma\xi_{n}+g(X^{0}_{n+1})\lambda^{0}_{n+1},\quad\zeta(X^{0}_{n+1})=0,$
(2.7)

where the Lagrange multiplier $\lambda^{0}_{n+1}$ is determined by the
constraint $\zeta(X^{0}_{n+1})=0$. The uniformly accurate discretization given
in (2.2) can be modified so that it converges to the integrator (2.7) when
$\varepsilon\rightarrow 0$. It suffices to replace the first line of (2.2)

$\displaystyle X_{n+1}^{\varepsilon}$
$\displaystyle=X_{n}^{\varepsilon}+\sqrt{h}\sigma\xi_{n}+hf(X_{n}^{\varepsilon})-\frac{(1-e^{-h/\varepsilon})^{2}}{2}(g^{\prime}(gG^{-1}\zeta)G^{-1}\zeta)(X_{n}^{\varepsilon})$
$\displaystyle+\frac{\sigma^{2}}{4}\Big{(}\sqrt{2\varepsilon
h(1-e^{-2h/\varepsilon})}-\frac{\varepsilon}{2}(1-e^{-2h/\varepsilon})\Big{)}\nabla\ln(\det(G))(X_{n}^{\varepsilon})+g(X_{n+1}^{\varepsilon})\lambda^{\varepsilon}_{n+1},$

and to keep the same expansion for the constraint
$\zeta(X_{n+1}^{\varepsilon})$. The methodology for the uniform expansion of
the integrator that we present in Section 3.2 extends to this context, so that
the convergence results persist. Similarly, one could change the direction of
projection $g(X_{n}^{\varepsilon})$ into $g(Y_{n}^{\varepsilon})$, where
$Y_{n}^{\varepsilon}$ is any consistent one-step approximation of
$X_{n}^{\varepsilon}$, in the spirit of the class of projected Runge-Kutta
methods presented in [30]. Finding a class of uniformly accurate
discretizations that converge to a more general class of Runge-Kutta methods
on the manifold $\mathcal{M}$ is matter for future works.

The uniform discretization (2.2) is implicit and requires one to solve a fixed
point problem at each step with, for instance, a fixed point iteration or a
Newton method. The following result, in the spirit of [22, Chap. VII] for
deterministic DAEs and [30, Lemma 3.3] for the constrained dynamics (1.2),
confirms that the associated implicit system is not stiff, that is, that its
complexity does not depend on the stiff parameter $\varepsilon$.

######  Theorem 2.12.

Under Assumption 2.1, each step of the integrator $(X_{n}^{\varepsilon})$
given by (2.2) can be rewritten as a solution of a fixed point problem of the
form

$X_{n+1}^{\varepsilon}=F_{h}^{\varepsilon}(X_{n+1}^{\varepsilon}),$

where $F_{h}^{\varepsilon}\colon\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$
depends on $X_{n}^{\varepsilon}$, $\xi_{n}$, $h$, and $\varepsilon$. Moreover,
if Assumption 2.4 is satisfied, then there exists $h_{0}>0$ independent of
$\varepsilon$ such that for all $h\leq h_{0}$, $F_{h}^{\varepsilon}$ is a
uniform contraction, that is, there exists a positive constant $L<1$
independent of $h$ and $\varepsilon$ such that, for all $y_{1}$,
$y_{2}\in\mathbb{R}^{d}$,

$\left|F_{h}^{\varepsilon}(y_{2})-F_{h}^{\varepsilon}(y_{1})\right|\leq
L\left|y_{2}-y_{1}\right|.$

##  3 Weak convergence analysis

In this section, we present the uniform weak expansion and the stability
properties of the solution $X^{\varepsilon}(t)$ to the penalized dynamic (1.4)
and of the uniform integrator $(X_{n}^{\varepsilon})$ given in (2.2). We then
use these results to prove Theorem 2.7, Theorem 2.9 and Theorem 2.12. Let us
begin the analysis with a few technical lemmas and notations.

######  Lemma 3.1.

Let $(X_{n}^{\varepsilon})$ be given by (2.2), then, under Assumption 2.1,
there exists a constant $C_{0}>0$ independent of $X_{0}\in\mathcal{M}$,
$\varepsilon$, and $h$ such that, for all $n\geq 0$,

$(1-e^{-h/\varepsilon})\left|\zeta(X^{\varepsilon}_{n})\right|\leq
C_{0}\sqrt{h}.$ (3.1)

###### Proof.

As $f$, $\xi_{n}$, $g$, and $g^{\prime}$ are bounded, we obtain from the
definition of the integrator given by (2.2) that

$|\zeta(X^{\varepsilon}_{n+1})|\leq
e^{-h/\varepsilon}|\zeta(X^{\varepsilon}_{n})|+C\sqrt{\varepsilon(1-e^{-2h/\varepsilon})},$

where we used that the function $(1-e^{-x})/x$ is bounded for $x>0$. Thus, as
$\zeta(X_{0})=0$,

$\displaystyle|\zeta(X^{\varepsilon}_{n})|$ $\displaystyle\leq
e^{-nh/\varepsilon}|\zeta(X_{0})|+C\sqrt{\varepsilon(1-e^{-2h/\varepsilon})}\sum_{k=0}^{n-1}e^{-kh/\varepsilon}$
$\displaystyle\leq
C\frac{\sqrt{\varepsilon(1-e^{-2h/\varepsilon})}}{1-e^{-h/\varepsilon}}\leq
C_{0}\frac{\sqrt{h}}{1-e^{-h/\varepsilon}},$

where the constant $C_{0}$ does not depend on $X_{0}$, $\varepsilon$, $n$, and
$h$. This yields the estimate (3.1). ∎

The estimate (3.1) is a direct consequence of our choice of using bounded
random variables in (2.2). We shall use this estimate extensively in the rest
of this section. Thus, we denote by $\mathcal{M}^{\varepsilon}_{h}$ the set of
vectors $x\in\mathbb{R}^{d}$ that satisfy the estimate

$(1-e^{-h/\varepsilon})\left|\zeta(x)\right|\leq C_{0}\sqrt{h},$ (3.2)

where $C_{0}$ is the constant given in Lemma 3.1. The set
$\mathcal{M}^{\varepsilon}_{h}$ is a closed subset of $\mathbb{R}^{d}$ that
contains $\mathcal{M}$. The numerical scheme given by (2.2) takes values in
$\mathcal{M}^{\varepsilon}_{h}$. We mention that the convergence results are
still valid if the initial condition $X_{0}$ of (2.2) is chosen in
$\mathcal{M}^{\varepsilon}_{h}$ instead of $\mathcal{M}$.

As we aim at writing uniform expansions, we introduce the convenient notation
$R^{\varepsilon}_{h}(x)$ for any remainder that satisfies at least
$\left|\mathbb{E}[R^{\varepsilon}_{h}(x)]\right|\leq Ch^{3/2}$, where $C$ is
independent of $\varepsilon$, $h$, and $x$.

The following result serves as a technical tool for simplifying the
calculations in the uniform expansions of Subsection 3.1 and Subsection 3.2.
It is proved with elementary computations.

######  Lemma 3.2.

The Fixman correction can be rewritten in the following way:

$\frac{\sigma^{2}}{4}\nabla\ln(\det(G))=\frac{\sigma^{2}}{2}\sum_{i=1}^{d}g^{\prime}(e_{i})G^{-1}g^{T}e_{i}=\frac{\sigma^{2}}{2}\sum_{i=1}^{d}g^{\prime}(gG^{-1}g^{T}e_{i})G^{-1}g^{T}e_{i},$

where $(e_{i})$ is the canonical basis of $\mathbb{R}^{d}$.

###### Proof of Lemma 3.2.

Using that $G=g^{T}g$ is a symmetric matrix, that $g=\nabla\zeta$ is a
gradient (which implies $x^{T}g^{\prime}(y)=y^{T}g^{\prime}(x)$) and the
standard properties of the trace operator $\operatorname{Tr}$, we deduce that

$\displaystyle\partial_{j}(\ln(\det(G)))$
$\displaystyle=\operatorname{Tr}(G^{-1}\partial_{j}G)=2\operatorname{Tr}(G^{-1}(g^{\prime}(e_{j}))^{T}g)=2\operatorname{Tr}(g^{\prime}(e_{j})G^{-1}g^{T})$
$\displaystyle=2\sum_{i=1}^{d}e_{i}^{T}g^{\prime}(e_{j})G^{-1}g^{T}e_{i}=2\sum_{i=1}^{d}e_{j}^{T}g^{\prime}(e_{i})G^{-1}g^{T}e_{i},$

that is,
$\nabla\ln(\det(G))=2\sum_{i=1}^{d}g^{\prime}(e_{i})G^{-1}g^{T}e_{i}$. For the
second equality, we have

$\displaystyle\sum_{i=1}^{d}e_{j}^{T}g^{\prime}(gG^{-1}g^{T}e_{i})G^{-1}g^{T}e_{i}$
$\displaystyle=\sum_{i=1}^{d}e_{i}^{T}gG^{-1}g^{T}g^{\prime}(e_{j})G^{-1}g^{T}e_{i}$
$\displaystyle=\operatorname{Tr}(gG^{-1}g^{T}g^{\prime}(e_{j})G^{-1}g^{T})=\operatorname{Tr}(g^{\prime}(e_{j})G^{-1}g^{T})$
$\displaystyle=\sum_{i=1}^{d}e_{i}^{T}g^{\prime}(e_{j})G^{-1}g^{T}e_{i}=\sum_{i=1}^{d}e_{j}^{T}g^{\prime}(e_{i})G^{-1}g^{T}e_{i}.$

Hence we get the result. ∎

###  3.1 Uniform expansion of the exact solution

We consider the exact solution $X^{\varepsilon}(t)$ of the penalized Langevin
dynamics (1.4), with the initial condition $X_{0}=x$, that we assume is
deterministic for simplicity. Then, $X^{\varepsilon}$ satisfies the following
expansion in $h$ that is uniform with respect to $\varepsilon$.

######  Proposition 3.3.

Under Assumption 2.1, there exists $h_{0}>0$ such that for all $h\leq h_{0}$,
if $X^{\varepsilon}$ is the solution of the penalized Langevin dynamics (1.4)
starting at $x\in\mathcal{M}^{\varepsilon}_{h}$, then, for all
$\phi\in\mathcal{C}^{3}_{P}$, the following estimate holds:

$\left|\mathbb{E}[\phi(X^{\varepsilon}(h))]-\mathbb{E}[\phi(x+\sqrt{h}A^{\varepsilon}_{h}(x)+hB^{\varepsilon}_{h}(x))]\right|\leq
C(1+\left|x\right|^{K})h^{3/2},$ (3.3)

where $C$ is independent of $h$ and $\varepsilon$ and where the functions
$A^{\varepsilon}_{h}$ and $B^{\varepsilon}_{h}$ are given by

$\displaystyle A^{\varepsilon}_{h}$
$\displaystyle=\sigma\xi+\frac{e^{-h/\varepsilon}-1}{\sqrt{h}}gG^{-1}\zeta+\sigma\Big{(}\sqrt{\frac{\varepsilon}{2h}(1-e^{-2h/\varepsilon})}-1\Big{)}gG^{-1}g^{T}\xi,$
$\displaystyle B^{\varepsilon}_{h}$
$\displaystyle=f+\Big{(}\frac{\varepsilon}{h}(1-e^{-h/\varepsilon})-1\Big{)}gG^{-1}g^{T}f+\frac{\sigma^{2}}{2}\Big{(}\frac{\varepsilon}{h}(1-e^{-h/\varepsilon})-1\Big{)}gG^{-1}\operatorname{div}(g)$
$\displaystyle+\frac{\sigma^{2}\varepsilon}{8h}(1-e^{-2h/\varepsilon})\nabla\ln(\det(G))+\frac{\sigma^{2}\varepsilon}{8h}(1-e^{-h/\varepsilon})^{2}gG^{-1}g^{T}\nabla\ln(\det(G))$
$\displaystyle+\frac{1}{2h}(e^{-h/\varepsilon}-1)^{2}\Big{(}g^{\prime}(gG^{-1}\zeta)G^{-1}\zeta-
gG^{-1}g^{T}g^{\prime}(gG^{-1}\zeta)G^{-1}\zeta$ $\displaystyle-
gG^{-1}(g^{\prime}(gG^{-1}\zeta))^{T}gG^{-1}\zeta\Big{)}+\sigma^{2}\Big{(}1+\frac{\varepsilon}{h}(e^{-h/\varepsilon}-1)\Big{)}\sum_{i=1}^{d}gG^{-1}(g^{\prime}(e_{i}))^{T}gG^{-1}g^{T}e_{i}$
$\displaystyle+\frac{\sigma^{2}}{4}\Big{(}\frac{\varepsilon}{h}(e^{-2h/\varepsilon}-4e^{-h/\varepsilon}+3)-2\Big{)}\sum_{i=1}^{d}gG^{-1}(g^{\prime}(gG^{-1}g^{T}e_{i}))^{T}gG^{-1}g^{T}e_{i},$

with $\xi$ a discrete bounded random vector that satisfies (2.3), and where
the functions $A^{\varepsilon}_{h}$ and $B^{\varepsilon}_{h}$ are bounded
uniformly in $\varepsilon$ and $h$ on $\mathcal{M}^{\varepsilon}_{h}$.

Note that for $q=1$, we obtain the first step of the uniformly accurate
discretization (2.4) by gathering all the terms of the form $gM$ with
$M\in\mathbb{R}^{q}$ of the weak approximation given in Proposition 3.3 in a
Lagrange multiplier $g\lambda^{\varepsilon}_{1}\in\mathbb{R}^{q}$ and by
adding the truncated expansion of $\zeta(X^{\varepsilon}(h))$ as a constraint.

The proof of Proposition 3.3 relies on the change of coordinate $\psi$ given
in Assumption 2.1. Instead of discretizing directly the penalized dynamics
(1.4), we first apply the change of coordinate $\psi$ and we derive an
expansion in time of $X^{\varepsilon}(t)$ that is uniform in the parameter
$\varepsilon$. The following result is used for proving Proposition 3.3.

######  Lemma 3.4.

With the same notations and assumptions as in Proposition 3.3, the following
estimates hold for all $h\leq h_{0}$ and all
$x\in\mathcal{M}^{\varepsilon}_{h}$:

$\displaystyle\mathbb{E}[|X^{\varepsilon}(h)-x|^{2}]^{1/2}$ $\displaystyle\leq
C\sqrt{h},$ (3.4)
$\displaystyle\mathbb{E}[|X^{\varepsilon}(h)-(x+\sqrt{h}\widehat{A}^{\varepsilon}_{h}(x))|^{2}]^{1/2}$
$\displaystyle\leq Ch,$ (3.5)
$\displaystyle\left|\mathbb{E}[\psi(X^{\varepsilon}(h))]-\mathbb{E}[\psi(x+\sqrt{h}\widehat{A}^{\varepsilon}_{h}(x)+hB^{\varepsilon}_{h}(x))]\right|$
$\displaystyle\leq Ch^{3/2},$ (3.6)

where $C$ is independent of $h$ and $\varepsilon$, $B^{\varepsilon}_{h}$ is
defined in Proposition 3.3, and $\widehat{A}^{\varepsilon}_{h}$ is given by

$\displaystyle\widehat{A}^{\varepsilon}_{h}$
$\displaystyle=\sigma\frac{W(h)}{\sqrt{h}}+\frac{e^{-h/\varepsilon}-1}{\sqrt{h}}gG^{-1}\zeta+\frac{\sigma}{\sqrt{h}}gG^{-1}g^{T}\int_{0}^{h}(e^{(s-h)/\varepsilon}-1)dW(s).$

###### Proof of Proposition 3.3.

The uniform bounds on $A^{\varepsilon}_{h}(X^{\varepsilon}(t))$ and
$B^{\varepsilon}_{h}(X^{\varepsilon}(t))$ are obtained straightforwardly by
using Assumption 2.1 and the fact that $x\in\mathcal{M}^{\varepsilon}_{h}$. We
prove the local weak order one of the approximation
$Y^{\varepsilon}(h)=x+\sqrt{h}\widehat{A}^{\varepsilon}_{h}(x)+hB^{\varepsilon}_{h}(x)$
given in Lemma 3.4. Let $\phi\in\mathcal{C}^{3}_{P}$ and
$\widetilde{\phi}=\phi\circ\psi^{-1}$; then a Taylor expansion around $x$
yields

$\displaystyle\Big{|}\mathbb{E}[\phi(X^{\varepsilon}(h))-\phi(Y^{\varepsilon}(h))]\Big{|}$
$\displaystyle=\Big{|}\mathbb{E}[\widetilde{\phi}\circ\psi(X^{\varepsilon}(h))]-\mathbb{E}[\widetilde{\phi}\circ\psi(Y^{\varepsilon}(h))]\Big{|}$
$\displaystyle\leq\Big{|}\mathbb{E}[\widetilde{\phi}^{\prime}(\psi(x))(\psi(X^{\varepsilon}(h))-\psi(Y^{\varepsilon}(h)))]\Big{|}$
$\displaystyle+\frac{1}{2}\Big{|}\mathbb{E}[\widetilde{\phi}^{\prime\prime}(\psi(x))(\psi(X^{\varepsilon}(h))-\psi(x),\psi(X^{\varepsilon}(h))-\psi(x))$
$\displaystyle-\widetilde{\phi}^{\prime\prime}(\psi(x))(\psi(Y^{\varepsilon}(h))-\psi(x),\psi(Y^{\varepsilon}(h))-\psi(x))]\Big{|}$
$\displaystyle+C(1+\left|x\right|^{K})h^{3/2}$
$\displaystyle\leq\Big{|}\widetilde{\phi}^{\prime}(\psi(x))\Big{|}\Big{|}\mathbb{E}[\psi(X^{\varepsilon}(h))-\psi(Y^{\varepsilon}(h))]\Big{|}$
$\displaystyle+\frac{1}{2}\Big{|}\widetilde{\phi}^{\prime\prime}(\psi(x))\Big{|}\mathbb{E}[|\psi(X^{\varepsilon}(h))+\psi(Y^{\varepsilon}(h))-2\psi(x)|^{2}]^{1/2}$
$\displaystyle\cdot\mathbb{E}[|\psi(X^{\varepsilon}(h))-\psi(Y^{\varepsilon}(h))|^{2}]^{1/2}+C(1+\left|x\right|^{K})h^{3/2}$

where we used (3.4), Assumption 2.1, and the bilinearity of
$\widetilde{\phi}^{\prime\prime}(\psi(x))$. With Lemma 3.4 and the regularity
properties of $\widetilde{\phi}$ and $\psi$, we get

$\left|\mathbb{E}[\phi(X^{\varepsilon}(h))]-\mathbb{E}[\phi(x+\sqrt{h}\widehat{A}^{\varepsilon}_{h}(x)+hB^{\varepsilon}_{h}(x))]\right|\leq
C(1+\left|x\right|^{K})h^{3/2}.$ (3.7)

In the spirit of [40, Chap. 2], we replace the random variable
$\widehat{A}^{\varepsilon}_{h}(x)$ by the random variable
$A^{\varepsilon}_{h}(x)$ that share the same expectation and covariance
matrix. Indeed, a calculation gives

$\displaystyle\operatorname{Cov}(\widehat{A}^{\varepsilon}_{h,i}(x),\widehat{A}^{\varepsilon}_{h,j}(x))$
$\displaystyle=\sigma^{2}\delta_{ij}+\frac{2\sigma^{2}}{h}(gG^{-1}g^{T})_{ij}\int_{0}^{h}(e^{(s-h)/\varepsilon}-1)ds$
$\displaystyle+\frac{\sigma^{2}}{h}\sum_{k=1}^{d}(gG^{-1}g^{T})_{ik}(gG^{-1}g^{T})_{jk}\int_{0}^{h}(e^{(s-h)/\varepsilon}-1)^{2}ds$
$\displaystyle=\sigma^{2}\delta_{ij}+\sigma^{2}(gG^{-1}g^{T})_{ij}\Big{(}\frac{\varepsilon}{2h}(1-e^{-2h/\varepsilon})-1\Big{)},$

where we used that $g^{T}g=G$ and the Itô isometry. On the other hand, a
similar calculation yields

$\displaystyle\operatorname{Cov}(A^{\varepsilon}_{h,i}(x),A^{\varepsilon}_{h,j}(x))$
$\displaystyle=\sigma^{2}\delta_{ij}+2\sigma^{2}\Big{(}\sqrt{\frac{\varepsilon}{2h}(1-e^{-2h/\varepsilon})}-1\Big{)}(gG^{-1}g^{T})_{ij}$
$\displaystyle+\sigma^{2}\Big{(}\sqrt{\frac{\varepsilon}{2h}(1-e^{-2h/\varepsilon})}-1\Big{)}^{2}\sum_{k=1}^{d}(gG^{-1}g^{T})_{ik}(gG^{-1}g^{T})_{jk}$
$\displaystyle=\sigma^{2}\delta_{ij}+\sigma^{2}(gG^{-1}g^{T})_{ij}\Big{(}\frac{\varepsilon}{2h}(1-e^{-2h/\varepsilon})-1\Big{)}.$

Replacing $\widehat{A}^{\varepsilon}_{h}(x)$ by $A^{\varepsilon}_{h}(x)$ in
the weak expansion (3.7) gives the estimate (3.3). ∎

The main ingredient of the proof of Lemma 3.4 is the decomposition of the
terms of the expansion in a part that stays on the tangent space and a part of
the form $gM$ with $M\in\mathbb{R}^{q}$ that is orthogonal to the tangent
space.

###### Proof of Lemma 3.4.

As $\psi^{-1}$ is Lipschitz, we have

$\displaystyle\mathbb{E}[|X^{\varepsilon}(h)-x|^{2}]^{1/2}$ $\displaystyle\leq
C\mathbb{E}[|\psi(X^{\varepsilon}(h))-\psi(x)|^{2}]^{1/2}$ $\displaystyle\leq
C\mathbb{E}[|\varphi(X^{\varepsilon}(h))-\varphi(x)|^{2}]^{1/2}+C\mathbb{E}[|\zeta(X^{\varepsilon}(h))-\zeta(x)|^{2}]^{1/2}.$

On the one hand, applying the Itô formula to $\varphi(X^{\varepsilon})$ yields

$\displaystyle\varphi(X^{\varepsilon}(h))$
$\displaystyle=\varphi(x)+\sigma\int_{0}^{h}\varphi^{\prime}(X^{\varepsilon}(s))dW(s)$
(3.8)
$\displaystyle+\int_{0}^{h}[\varphi^{\prime}f+\frac{\sigma^{2}}{4}\varphi^{\prime}\nabla\ln(\det(G))+\frac{\sigma^{2}}{2}\Delta\varphi](X^{\varepsilon}(s))ds,$

where the term in $\varepsilon$ vanishes as $\varphi^{\prime}g=0$ (see
Assumption 2.1). Assumption 2.1 allows us to write the uniform strong
expansion

$\mathbb{E}[|\varphi(X^{\varepsilon}(h))-\varphi(x)|^{2}]^{1/2}\leq
C\sqrt{h}.$

On the other hand, for $\zeta(X^{\varepsilon})$, we have

$d\zeta(X^{\varepsilon})=\sigma
g^{T}(X^{\varepsilon})dW+\Big{[}g^{T}f+\frac{\sigma^{2}}{4}g^{T}\nabla\ln(\det(G))+\frac{\sigma^{2}}{2}\operatorname{div}(g)-\frac{1}{\varepsilon}\zeta\Big{]}(X^{\varepsilon})dt.$

With the variation of constants formula, it rewrites into

$\displaystyle\zeta(X^{\varepsilon}(h))$
$\displaystyle=e^{-h/\varepsilon}\zeta(x)+\sigma\int_{0}^{h}e^{(s-h)/\varepsilon}g^{T}(X^{\varepsilon}(s))dW(s)$
(3.9)
$\displaystyle+\int_{0}^{h}e^{(s-h)/\varepsilon}\Big{[}g^{T}f+\frac{\sigma^{2}}{4}g^{T}\nabla\ln(\det(G))+\frac{\sigma^{2}}{2}\operatorname{div}(g)\Big{]}(X^{\varepsilon}(s))ds.$

As the integrands in (3.9) are bounded (using Assumption 2.1), we get

$\mathbb{E}[|\zeta(X^{\varepsilon}(h))-\zeta(x)|^{2}]^{1/2}\leq
C((e^{-h/\varepsilon}-1)^{2}\left|\zeta(x)\right|^{2}+h)^{1/2}\leq C\sqrt{h},$

where we used that $x\in\mathcal{M}^{\varepsilon}_{h}$. We thus get the
desired estimate (3.4). The estimate (3.5) is obtained with the same arguments
by keeping track of the terms of size $\mathcal{O}(\sqrt{h})$ in the
expansions.

We now prove the weak estimate (3.6). We denote for simplicity
$Y^{\varepsilon}(h)=x+\sqrt{h}\widehat{A}^{\varepsilon}_{h}(x)+hB^{\varepsilon}_{h}(x)$.
Let us first look at the approximation of $\varphi(X^{\varepsilon}(h))$. On
the one hand, applying the Itô formula to $\varphi(X^{\varepsilon}(t))$ gives

$\displaystyle\varphi(X^{\varepsilon}(h))$
$\displaystyle=\varphi(x)+h\varphi^{\prime}f(x)+h\frac{\sigma^{2}}{4}\varphi^{\prime}\nabla\ln(\det(G))(x)+h\frac{\sigma^{2}}{2}\Delta\varphi(x)+R^{\varepsilon}_{h}(x),$

where we used (3.4) and we put in $R^{\varepsilon}_{h}(x)$ all the terms that
are zero in average. On the other hand, an expansion in $h$ of
$\varphi(Y^{\varepsilon}(h))$ yields

$\displaystyle\varphi(Y^{\varepsilon}(h))$
$\displaystyle=\varphi+\sqrt{h}\varphi^{\prime}\widehat{A}^{\varepsilon}_{h}+h\Big{[}\varphi^{\prime}B^{\varepsilon}_{h}+\frac{1}{2}\varphi^{\prime\prime}(\widehat{A}^{\varepsilon}_{h},\widehat{A}^{\varepsilon}_{h})\Big{]}+R^{\varepsilon}_{h}$
$\displaystyle=\varphi+h\varphi^{\prime}f+\frac{\sigma^{2}\varepsilon}{8}(1-e^{-2h/\varepsilon})\varphi^{\prime}\nabla\ln(\det(G))+\frac{1}{2}(e^{-h/\varepsilon}-1)^{2}\varphi^{\prime}g^{\prime}(gG^{-1}\zeta)G^{-1}\zeta$
$\displaystyle+\frac{\sigma^{2}}{2}\varphi^{\prime\prime}(W(h),W(h))+\frac{1}{2}(e^{-h/\varepsilon}-1)^{2}\varphi^{\prime\prime}(gG^{-1}\zeta,gG^{-1}\zeta)$
$\displaystyle+\frac{\sigma^{2}}{2}\int_{0}^{h}\int_{0}^{h}(e^{(s-h)/\varepsilon}-1)(e^{(u-h)/\varepsilon}-1)\varphi^{\prime\prime}(gG^{-1}g^{T}dW(s),gG^{-1}g^{T}dW(u))$
$\displaystyle+\sigma^{2}\int_{0}^{h}(e^{(s-h)/\varepsilon}-1)\varphi^{\prime\prime}(W(h),gG^{-1}g^{T}dW(s))+R^{\varepsilon}_{h},$

where we use that $\varphi^{\prime}g=0$, we omit the dependence in $x$ for
conciseness, and we put in $R^{\varepsilon}_{h}(x)$ all the terms that are
zero in average. We now replace the random terms by their expectation,

$\displaystyle\varphi(Y^{\varepsilon}(h))$
$\displaystyle=\varphi+h\varphi^{\prime}f+\frac{\sigma^{2}\varepsilon}{8}(1-e^{-2h/\varepsilon})\varphi^{\prime}\nabla\ln(\det(G))+h\frac{\sigma^{2}}{2}\Delta\varphi$
$\displaystyle+\frac{1}{2}(e^{-h/\varepsilon}-1)^{2}(\varphi^{\prime}(g^{\prime}(gG^{-1}\zeta)G^{-1}\zeta)+\varphi^{\prime\prime}(gG^{-1}\zeta,gG^{-1}\zeta))$
$\displaystyle+\frac{\sigma^{2}}{2}\int_{0}^{h}(e^{(s-h)/\varepsilon}-1)^{2}ds\sum_{i=1}^{d}\varphi^{\prime\prime}(gG^{-1}g^{T}e_{i},gG^{-1}g^{T}e_{i})$
$\displaystyle+\sigma^{2}\int_{0}^{h}(e^{(s-h)/\varepsilon}-1)ds\sum_{i=1}^{d}\varphi^{\prime\prime}(e_{i},gG^{-1}g^{T}e_{i})+R^{\varepsilon}_{h}.$

Letting $M\in\mathbb{R}^{q}$, we differentiate the equality
$\varphi^{\prime}(gM)=0$. We obtain that for any $M\in\mathbb{R}^{q}$ and
$v\in\mathbb{R}^{d}$ we have
$\varphi^{\prime}(g^{\prime}(v)M)+\varphi^{\prime\prime}(gM,v)=0$. We deduce
that

$\varphi^{\prime}(g^{\prime}(gG^{-1}\zeta)G^{-1}\zeta)+\varphi^{\prime\prime}(gG^{-1}\zeta,gG^{-1}\zeta)=0.$

Applying Lemma 3.2, we get by a direct calculation that

$\varphi(Y^{\varepsilon}(h))=\varphi(x)+h\varphi^{\prime}f(x)+h\frac{\sigma^{2}}{4}\varphi^{\prime}\nabla\ln(\det(G))(x)+h\frac{\sigma^{2}}{2}\Delta\varphi(x)+R^{\varepsilon}_{h}(x),$

which gives the desired estimate

$\left|\mathbb{E}[\varphi(X^{\varepsilon}(h))]-\mathbb{E}[\varphi(x+\sqrt{h}\widehat{A}^{\varepsilon}_{h}(x)+hB^{\varepsilon}_{h}(x))]\right|\leq
Ch^{3/2}.$ (3.10)

For the one-step approximation of $\zeta(X^{\varepsilon}(h))$, the Itô formula
and the variation of constants formula yield

$\zeta(X^{\varepsilon}(h))=e^{-h/\varepsilon}\zeta(x)+\varepsilon(1-e^{-h/\varepsilon})(g^{T}f+\frac{\sigma^{2}}{4}g^{T}\nabla\ln(\det(G))+\frac{\sigma^{2}}{2}\operatorname{div}(g))(x)+R^{\varepsilon}_{h}(x).$

For $\zeta(Y^{\varepsilon}(h))$, using $g^{T}g=G$, we get with the same
arguments as for $\varphi(Y^{\varepsilon}(h))$ that

$\displaystyle\zeta(Y^{\varepsilon}(h))$
$\displaystyle=\zeta+\sqrt{h}g^{T}\widehat{A}^{\varepsilon}_{h}+h\Big{[}g^{T}B^{\varepsilon}_{h}+\frac{1}{2}(g^{\prime}(\widehat{A}^{\varepsilon}_{h}))^{T}\widehat{A}^{\varepsilon}_{h}\Big{]}+R^{\varepsilon}_{h}$
$\displaystyle=e^{-h/\varepsilon}\zeta+\varepsilon(1-e^{-h/\varepsilon})g^{T}f+\frac{\sigma^{2}}{2}\Big{(}\varepsilon(1-e^{-h/\varepsilon})-h\Big{)}\operatorname{div}(g)+\frac{\sigma^{2}}{2}(g^{\prime}(W(h)))^{T}W(h)$
$\displaystyle+\frac{\sigma^{2}\varepsilon}{4}(1-e^{-h/\varepsilon})^{2}g^{T}\nabla\ln(\det(G))+\sigma^{2}\Big{(}h+\varepsilon(e^{-h/\varepsilon}-1)\Big{)}\sum_{i=1}^{d}(g^{\prime}(e_{i}))^{T}gG^{-1}g^{T}e_{i}$
$\displaystyle+\frac{\sigma^{2}}{4}\Big{(}\varepsilon(e^{-2h/\varepsilon}-4e^{-h/\varepsilon}+3)-2h\Big{)}\sum_{i=1}^{d}(g^{\prime}(gG^{-1}g^{T}e_{i}))^{T}gG^{-1}g^{T}e_{i}$
$\displaystyle+\frac{\sigma^{2}}{2}\int_{0}^{h}\int_{0}^{h}(e^{(s-h)/\varepsilon}-1)(e^{(u-h)/\varepsilon}-1)(g^{\prime}(gG^{-1}g^{T}dW(u)))^{T}gG^{-1}g^{T}dW(s)$
$\displaystyle+\sigma^{2}\int_{0}^{h}\int_{0}^{h}(e^{(s-h)/\varepsilon}-1)(g^{\prime}(W(h)))^{T}gG^{-1}g^{T}dW(s)+R^{\varepsilon}_{h}.$

We now replace the stochastic integrals by their expectations (putting the
remainders in $R^{\varepsilon}_{h}$), and we use Lemma 3.2 to simplify the
expansion. It yields

$\displaystyle\zeta(Y^{\varepsilon}(h))$
$\displaystyle=e^{-h/\varepsilon}\zeta+\varepsilon(1-e^{-h/\varepsilon})(g^{T}f+\frac{\sigma^{2}}{4}g^{T}\nabla\ln(\det(G))+\frac{\sigma^{2}}{2}\operatorname{div}(g))+R^{\varepsilon}_{h},$

which implies

$\left|\mathbb{E}[\zeta(X^{\varepsilon}(h))]-\mathbb{E}[\zeta(x+\sqrt{h}\widehat{A}^{\varepsilon}_{h}(x)+hB^{\varepsilon}_{h}(x))]\right|\leq
Ch^{3/2}.$ (3.11)

Combining the inequalities (3.10) and (3.11) gives the desired weak estimate
(3.6). ∎

To end this subsection, we recall the growth properties of $X^{\varepsilon}$
that will be of use in Subsection 3.3. For this particular result, we add the
dependency in the initial condition $x$ of the exact solution of (1.4) with
the notation $X^{\varepsilon}(t,x)$.

######  Lemma 3.5.

Under Assumption 2.1, for $\phi\in\mathcal{C}^{5}_{P}$ and $t\leq T$ fixed,
the function
$\widetilde{\phi}^{\varepsilon}(x)=\mathbb{E}[\phi(X^{\varepsilon}(t,x))]$
lies in $\mathcal{C}^{3}_{P}$ with constants independent of $t$ and
$\varepsilon$.

###### Proof.

The standard theory (see, for instance, the textbook [40]) gives
$\widetilde{\phi}^{\varepsilon}\in\mathcal{C}^{3}$. With the regularity
assumptions on $\psi$ and its derivatives, it is sufficient to prove that
$\mathbb{E}[\phi(Y^{\varepsilon}(t,y))]$ is in $\mathcal{C}^{3}_{P}$, where
$Y^{\varepsilon}(t,y)=\psi(X^{\varepsilon}(t,\psi^{-1}(x)))$ (replacing $\phi$
by $\phi\circ\psi^{-1}$). We recall from the proof of Lemma 3.4 that
$\varphi(X^{\varepsilon}(t,x))$ satisfies the integral formulation (3.8) and
that $\zeta(X^{\varepsilon}(t,x))$ satisfies (3.9). Putting together these
equations, we deduce that $Y^{\varepsilon}(t,y)$ satisfies an equation of the
form

$Y^{\varepsilon}(t,y)=A^{\varepsilon}(t)y+\int_{0}^{t}A^{\varepsilon}(t-s)F(Y^{\varepsilon}(s,y))ds+\int_{0}^{t}A^{\varepsilon}(t-s)G(Y^{\varepsilon}(s,y))dW(s),$

where $F$ and $G$ are in $\mathcal{C}^{3}_{P}$ and do not depend on
$\varepsilon$, and

$A^{\varepsilon}(t)=\begin{pmatrix}I_{d-q}&0\\\
0&e^{-t/\varepsilon}I_{q}\end{pmatrix}.$

As $\left|A^{\varepsilon}(t)\right|\leq C$, the process $Y^{\varepsilon}(t,y)$
satisfies $\mathbb{E}[\left|Y^{\varepsilon}(t,y)\right|^{2p}]\leq
C(1+\left|y\right|^{K})$, with $C$ and $K$ independent of $\varepsilon$. Thus
we have

$\mathbb{E}[\phi(Y^{\varepsilon}(t,y))]\leq
C(1+\mathbb{E}[\left|Y^{\varepsilon}(t,y)\right|^{K}])\leq
C(1+\mathbb{E}[\left|y\right|^{K}]).$

For the derivatives, we recall from [19] that
$Z^{\varepsilon}(t,y)=\partial_{y}Y^{\varepsilon}(t,y)$ satisfies the equation

$\displaystyle Z^{\varepsilon}(t,y)$
$\displaystyle=A^{\varepsilon}(t)I_{d}+\int_{0}^{t}A^{\varepsilon}(t-s)F^{\prime}(Y^{\varepsilon}(s,y))Z^{\varepsilon}(t,y)ds$
$\displaystyle+\int_{0}^{t}A^{\varepsilon}(t-s)G^{\prime}(Y^{\varepsilon}(s,y))Z^{\varepsilon}(t,y)dW(s).$

Applying the same arguments as for $Y^{\varepsilon}(t,y)$ yields that
$Z^{\varepsilon}(t,y)$ has bounded moments of all order and that
$\mathbb{E}[\phi(Z^{\varepsilon}(t,y))]\leq
C(1+\mathbb{E}[\left|y\right|^{K}])$. The same methodology extends to
$\partial_{y}^{2}Y^{\varepsilon}(t,y)$ and
$\partial_{y}^{3}Y^{\varepsilon}(t,y)$. ∎

###  3.2 Uniform expansion and bounded moments of the numerical solution

In this subsection, we show that the integrator given by (2.2) has bounded
moments of all order, that it lies on the manifold $\mathcal{M}$ in the limit
$\varepsilon\rightarrow 0$, and that it satisfies the same local weak uniform
expansion as the exact solution of (1.4) (see Proposition 3.3).

First, the integrator given by (2.2) satisfies the following bounded moments
property.

######  Proposition 3.6.

Under Assumption 2.1, $(X_{n}^{\varepsilon})$ has bounded moments of any order
along time; i.e., for all timestep $h\leq h_{0}$ small enough such that $Nh=T$
is fixed, for all integer $m\geq 0$,

$\sup_{n\leq N}\mathbb{E}[\left|X_{n}^{\varepsilon}\right|^{2m}]\leq C_{m},$

where the constant $C>0$ is independent of $\varepsilon$ and $h$.

The integrator $(X_{n}^{\varepsilon})$ given in (2.2) also satisfies a uniform
local expansion that is similar to its continuous counterpart presented in
Proposition 3.3.

######  Proposition 3.7.

Under Assumption 2.1, there exists $h_{0}>0$ such that for all $h\leq h_{0}$,
if $(X^{\varepsilon}_{n})$ is the numerical discretization given by (2.2)
beginning at $x\in\mathcal{M}^{\varepsilon}_{h}$ (assumed deterministic for
simplicity), then, for all test functions $\phi\in\mathcal{C}^{3}_{P}$, the
following estimate holds:

$\left|\mathbb{E}[\phi(X^{\varepsilon}_{1})]-\mathbb{E}[\phi(x+\sqrt{h}A^{\varepsilon}_{h}(x)+hB^{\varepsilon}_{h}(x))]\right|\leq
C(1+\left|x\right|^{K})h^{3/2},$ (3.12)

where $C$ is independent of $h$ and $\varepsilon$, and $A^{\varepsilon}_{h}$
and $B^{\varepsilon}_{h}$ are the functions given in Proposition 3.3.

To prove Proposition 3.6 and Proposition 3.7, we rely on the following lemma,
whose proof is postponed to the end of this subsection. We emphasize that an
inequality of the form $\left|\zeta(x)\right|\leq C$ does not imply in general
that $x$ stays close to $\mathcal{M}$. That is why we rely in Lemma 3.8 on an
estimate of the Lagrange multipliers (using the inequality (2.1)). This
estimate ensures that the method evolves in a neighborhood of the manifold.

######  Lemma 3.8.

Under Assumption 2.1 and if $x\in\mathcal{M}^{\varepsilon}_{h}$, there exists
$h_{0}>0$ such that, for all timestep $h\leq h_{0}$, the one-step
approximation $X_{1}^{\varepsilon}$ and the Lagrange multiplier
$\lambda^{\varepsilon}_{1}$ in the discretization (2.2) satisfy

$\left|X_{1}^{\varepsilon}-x\right|\leq
C\sqrt{h},\qquad\lambda^{\varepsilon}_{1}=\sqrt{h}G^{-1}(x)\lambda_{1,(1/2)}^{\varepsilon}+hG^{-1}(x)\lambda_{1,(1)}^{\varepsilon}+R^{\varepsilon}_{h}(x),$

where $|\lambda_{1,(1/2)}^{\varepsilon}|\leq C$,
$|\lambda_{1,(1)}^{\varepsilon}|\leq C$, and $|R^{\varepsilon}_{h}(x)|\leq
Ch^{3/2}$ with $C$ independent of $\varepsilon$, $h$ and $x$.

For proving Proposition 3.6, we apply the change of variable $\psi$ and we
adapt the standard methodology presented in [40, Lemma 1.1.6 & Lemma 2.2.2].

###### Proof of Proposition 3.6.

We derive from (2.2) that

$\displaystyle\left|\mathbb{E}[\zeta(X_{n+1}^{\varepsilon})-e^{-h/\varepsilon}\zeta(X_{n}^{\varepsilon})|X_{n}^{\varepsilon}]\right|$
$\displaystyle\leq Ch,$
$\displaystyle\left|\zeta(X_{n+1}^{\varepsilon})-e^{-h/\varepsilon}\zeta(X_{n}^{\varepsilon})\right|$
$\displaystyle\leq C\sqrt{h}.$

We prove that $\zeta(X_{n})$ has bounded moments by induction on $n$. The
binomial formula yields

$\displaystyle\mathbb{E}[\left|\zeta(X_{n+1}^{\varepsilon})\right|^{2m}]$
$\displaystyle=\mathbb{E}[\left|e^{-h/\varepsilon}\zeta(X^{\varepsilon}_{n})+\zeta(X^{\varepsilon}_{n+1})-e^{-h/\varepsilon}\zeta(X^{\varepsilon}_{n})\right|^{2m}]$
$\displaystyle\leq
e^{-2mh/\varepsilon}\mathbb{E}[\left|\zeta(X_{n}^{\varepsilon})\right|^{2m}]$
$\displaystyle+C\mathbb{E}[\left|\zeta(X_{n}^{\varepsilon})\right|^{2m-1}\left|\mathbb{E}[\zeta(X_{n+1}^{\varepsilon})-e^{-h/\varepsilon}\zeta(X_{n}^{\varepsilon})|X_{n}^{\varepsilon}]\right|]$
$\displaystyle+C\sum_{k=2}^{2m}\mathbb{E}[\left|\zeta(X_{n}^{\varepsilon})\right|^{2m-k}\left|\zeta(X_{n+1}^{\varepsilon})-e^{-h/\varepsilon}\zeta(X_{n}^{\varepsilon})\right|^{k}]$
$\displaystyle\leq\mathbb{E}[\left|\zeta(X_{n}^{\varepsilon})\right|^{2m}]+C(1+\mathbb{E}[\left|\zeta(X_{n}^{\varepsilon})\right|^{2m}])h.$

Following [40, Lemma 1.1.6], as $X_{0}\in\mathcal{M}$ is bounded, it implies
that $\mathbb{E}[\left|\zeta(X_{n})\right|^{2m}]$ is bounded uniformly in
$n=0,\dots,N$ and $\varepsilon$.

Using Lemma 3.8 and the equality $\varphi^{\prime}g=0$, a direct calculation
gives

$\displaystyle\left|\mathbb{E}[\varphi(X_{n+1}^{\varepsilon})-\varphi(X_{n}^{\varepsilon})|X_{n}^{\varepsilon}]\right|$
$\displaystyle\leq Ch,$ (3.13)
$\displaystyle\left|\varphi(X_{n+1}^{\varepsilon})-\varphi(X_{n}^{\varepsilon})\right|$
$\displaystyle\leq C\sqrt{h}.$ (3.14)

Following the same methodology as for
$\mathbb{E}[\left|\zeta(X_{n})\right|^{2m}]$, the estimates (3.13)-(3.14)
imply that the quantity $\mathbb{E}[\left|\varphi(X_{n})\right|^{2m}]$ is
bounded uniformly in $n=0,\dots,N$ and $\varepsilon$. Then, as $\psi^{-1}$ is
Lipschitz, we have

$\mathbb{E}[\left|X_{n}^{\varepsilon}\right|^{2m}]\leq
C(1+\mathbb{E}[\left|\psi(X_{n}^{\varepsilon})\right|^{2m}])\leq
C(1+\mathbb{E}[\left|\varphi(X_{n}^{\varepsilon})\right|^{2m}]+\mathbb{E}[\left|\zeta(X_{n}^{\varepsilon})\right|^{2m}])\leq
C.$

Hence we get the result. ∎

We obtain the uniform expansion of the numerical solution by writing
explicitly a uniform expansion of the Lagrange multiplier
$\lambda^{\varepsilon}_{1}$, in the spirit of [34, Lemma 3.25].

###### Proof of Proposition 3.7.

Using Lemma 3.8 and Assumption 2.1, we obtain

$X^{\varepsilon}_{1}=x+\sqrt{h}\Big{[}\sigma\xi+(gG^{-1})(x)\lambda_{1,(1/2)}^{\varepsilon}\Big{]}+R_{h}^{\varepsilon}(x),$

where the remainder satisfies $\left|R_{h}^{\varepsilon}(x)\right|\leq Ch$.
The constraint is then given by

$\displaystyle\zeta(X^{\varepsilon}_{1})$
$\displaystyle=\zeta(x)+\sqrt{h}\Big{[}\sigma
g^{T}(x)\xi+\lambda_{1,(1/2)}^{\varepsilon}\Big{]}+R_{h}^{\varepsilon}(x).$

On the other hand, we get from the definition of the integrator (2.2) that

$\displaystyle\zeta(X^{\varepsilon}_{1})$
$\displaystyle=\zeta(x)+\sqrt{h}\Big{[}\frac{e^{-h/\varepsilon}-1}{\sqrt{h}}\zeta(x)+\sigma\sqrt{\frac{\varepsilon}{2h}(1-e^{-2h/\varepsilon})}g^{T}(x)\xi\Big{]}+R_{h}^{\varepsilon}(x).$

By identifying the two terms in $\sqrt{h}$ in the expansions of
$\zeta(X^{\varepsilon}_{1})$, we deduce the value of
$\lambda_{1,(1/2)}^{\varepsilon}$, that is,

$\lambda_{1,(1/2)}^{\varepsilon}=\frac{e^{-h/\varepsilon}-1}{\sqrt{h}}\zeta+\sigma\Big{(}\sqrt{\frac{\varepsilon}{2h}(1-e^{-2h/\varepsilon})}-1\Big{)}g^{T}\xi.$

The expression of $\lambda_{1,(1/2)}^{\varepsilon}$ and Lemma 3.8 give

$\displaystyle X^{\varepsilon}_{1}$
$\displaystyle=x+\sqrt{h}A^{\varepsilon}_{h}(x)+h\Big{[}f(x)+\frac{(1-e^{-h/\varepsilon})^{2}}{2h}(g^{\prime}(gG^{-1}\zeta)G^{-1}\zeta)(x)$
$\displaystyle+\frac{\sigma^{2}\varepsilon}{8h}(1-e^{-2h/\varepsilon})\nabla\ln(\det(G))(x)+(gG^{-1})(x)\lambda_{1,(1)}^{\varepsilon}\Big{]}+R_{h}^{\varepsilon}(x),$
(3.15)

where $\left|R_{h}^{\varepsilon}(x)\right|\leq Ch^{3/2}$. We then compute the
expansion of $\zeta(X^{\varepsilon}_{1})$, and we compare it with the
definition of the integrator (2.2) to obtain the expression of
$\lambda_{1,(1)}^{\varepsilon}$. Inserting this expression in (3.15) gives

$X^{\varepsilon}_{1}=x+\sqrt{h}A^{\varepsilon}_{h}(x)+hB^{\varepsilon}_{h}(x)+R_{h}^{\varepsilon}(x),$

where the remainder satisfies
$\left|\mathbb{E}[R_{h}^{\varepsilon}(x)]\right|\leq Ch^{3/2}$ and
$\left|R_{h}^{\varepsilon}(x)\right|\leq Ch$. A Taylor expansion of
$\phi(X^{\varepsilon}_{1})$ around
$x+\sqrt{h}A^{\varepsilon}_{h}(x)+hB^{\varepsilon}_{h}(x)$ yields the estimate
(3.12). ∎

The proof of Lemma 3.8 mainly relies on the estimates (2.1) and (3.2). We
refer the reader to [30, Lemma 3.3] and [22, Chap. VII] for similar proofs
where explicit expressions of Lagrange multipliers are derived.

###### Proof of Lemma 3.8.

Let $x\in\mathcal{M}^{\varepsilon}_{h}$. For brevity, we rewrite the
discretization (2.2) as

$\displaystyle X_{1}^{\varepsilon}$
$\displaystyle=x+\sqrt{h}Y_{(1/2)}^{\varepsilon}(x)+hY_{(1)}^{\varepsilon}(x)+g(x)\lambda^{\varepsilon}_{1},$
$\displaystyle\zeta(X_{1}^{\varepsilon})$
$\displaystyle=\zeta(x)+\sqrt{h}\zeta_{(1/2)}^{\varepsilon}(x)+h\zeta_{(1)}^{\varepsilon}(x),$

where the functions $Y_{(1/2)}^{\varepsilon}$, $Y_{(1)}^{\varepsilon}$,
$\zeta_{(1/2)}^{\varepsilon}$, and $\zeta_{(1)}^{\varepsilon}$ are given by

$\displaystyle Y_{(1/2)}^{\varepsilon}$ $\displaystyle=\sigma\xi,$
$\displaystyle Y_{(1)}^{\varepsilon}$
$\displaystyle=f+\frac{(1-e^{-h/\varepsilon})^{2}}{2h}(g^{\prime}(gG^{-1}\zeta)G^{-1}\zeta)+\frac{\sigma^{2}\varepsilon}{8h}(1-e^{-2h/\varepsilon})\nabla\ln(\det(G)),$
$\displaystyle\zeta_{(1/2)}^{\varepsilon}$
$\displaystyle=\frac{e^{-h/\varepsilon}-1}{\sqrt{h}}\zeta+\sigma\sqrt{\frac{\varepsilon}{2h}(1-e^{-2h/\varepsilon})}g^{T}\xi,$
$\displaystyle\zeta_{(1)}^{\varepsilon}$
$\displaystyle=\frac{\varepsilon}{h}(1-e^{-h/\varepsilon})\big{(}g^{T}f+\frac{\sigma^{2}}{4}g^{T}\nabla\ln(\det(G))+\frac{\sigma^{2}}{2}\operatorname{div}(g)\big{)}$
$\displaystyle+\sigma^{2}\Big{(}\frac{\varepsilon}{h}(1-e^{-h/\varepsilon})-\sqrt{\frac{\varepsilon}{2h}(1-e^{-2h/\varepsilon})}\Big{)}\sum_{i=1}^{d}((g^{\prime}(gG^{-1}g^{T}e_{i}))^{T}gG^{-1}g^{T}e_{i})$
$\displaystyle-\sigma^{2}\Big{(}\frac{\varepsilon}{h}(1-e^{-h/\varepsilon})-\sqrt{\frac{\varepsilon}{2h}(1-e^{-2h/\varepsilon})}\Big{)}\sum_{i=1}^{d}((g^{\prime}(e_{i}))^{T}gG^{-1}g^{T}e_{i}).$

Using Assumption 2.1 and the estimate (3.2), the following uniform estimates
hold:

$\left|Y_{(i)}^{\varepsilon}(x)\right|\leq
C,\quad\left|\zeta_{(i)}^{\varepsilon}(x)\right|\leq C,\quad i\in\\{1/2,1\\}.$

The fundamental theorem of calculus yields

$\zeta(X_{1}^{\varepsilon})-\zeta(x)=\int_{0}^{1}g^{T}(x+\tau(X_{1}^{\varepsilon}-x))d\tau(X_{1}^{\varepsilon}-x)=\sqrt{h}\zeta_{(1/2)}^{\varepsilon}(x)+h\zeta_{(1)}^{\varepsilon}(x).$

Substituting $X_{1}^{\varepsilon}-x$ then gives

$\int_{0}^{1}g^{T}(x+\tau(X_{1}^{\varepsilon}-x))d\tau(\sqrt{h}Y_{(1/2)}^{\varepsilon}(x)+hY_{(1)}^{\varepsilon}(x)+g(x)\lambda^{\varepsilon}_{1})=\sqrt{h}\zeta_{(1/2)}^{\varepsilon}(x)+h\zeta_{(1)}^{\varepsilon}(x).$

Using Assumption 2.1, we get the following explicit expression of
$\lambda^{\varepsilon}_{1}$:

$\displaystyle\lambda^{\varepsilon}_{1}$
$\displaystyle=\sqrt{h}G_{X_{1}^{\varepsilon}-x}^{-1}(x)\Big{(}\zeta_{(1/2)}^{\varepsilon}(x)-\int_{0}^{1}g^{T}(x+\tau(X_{1}^{\varepsilon}-x))d\tau
Y_{(1/2)}^{\varepsilon}(x)\Big{)}$ (3.16)
$\displaystyle+hG_{X_{1}^{\varepsilon}-x}^{-1}(x)\Big{(}\zeta_{(1)}^{\varepsilon}(x)-\int_{0}^{1}g^{T}(x+\tau(X_{1}^{\varepsilon}-x))d\tau
Y_{(1)}^{\varepsilon}(x)\Big{)}.$

Then, the growth assumption (2.1) on $G_{y}^{-1}(x)$ allows us to write

$\left|\lambda^{\varepsilon}_{1}\right|\leq
C\sqrt{h}(1+\left|X_{1}^{\varepsilon}-x\right|)\quad\text{and}\quad\left|X_{1}^{\varepsilon}-x\right|\leq
C\sqrt{h}(1+\left|X_{1}^{\varepsilon}-x\right|).$

Hence, for $h\leq h_{0}$ small enough, we deduce that
$\left|X_{1}^{\varepsilon}-x\right|\leq C\sqrt{h}$ and

$\lambda^{\varepsilon}_{1}=\sqrt{h}G^{-1}(x)(\zeta_{(1/2)}^{\varepsilon}-g^{T}Y_{(1/2)}^{\varepsilon})(x)+R^{\varepsilon}_{h},$
(3.17)

where $\left|R^{\varepsilon}_{h}\right|\leq Ch$. For the term of size
$\mathcal{O}(h)$, we first deduce from (3.17) that

$\left|X_{1}^{\varepsilon}-x-\sqrt{h}(Y_{(1/2)}^{\varepsilon}+gG^{-1}\zeta_{(1/2)}^{\varepsilon}-gG^{-1}g^{T}Y_{(1/2)}^{\varepsilon})(x)\right|\leq
Ch.$

By using this estimate in (3.16), a Taylor expansion yields the desired
expansion of $\lambda^{\varepsilon}_{1}$. ∎

###  3.3 Proofs of the convergence theorems

Now that we have the local uniform expansion of the exact solution and the
numerical scheme, as well as the stability property of Proposition 3.6, we are
able to prove the main convergence theorems.

###### Proof of Theorem 2.7.

We derive the global weak consistency (2.5) with techniques similar to the
ones presented in [40, Chap. 2]. We denote by $X^{\varepsilon}(t,x)$ the
solution of the penalized dynamics with initial condition $x$ and
$X_{n}^{\varepsilon}(x)$ the numerical solution with initial condition $x$.
For $x\in\mathcal{M}^{\varepsilon}_{h}$, Proposition 3.3 and Proposition 3.7
yield

$\left|\mathbb{E}[\phi(X^{\varepsilon}(h,x))-\phi(X_{1}^{\varepsilon}(x))|x]\right|\leq
C(1+\left|x\right|^{K})h^{3/2},$

where $\phi\in\mathcal{C}^{3}_{P}$. Lemma 3.5 gives that
$\phi_{n}(x)=\mathbb{E}[\phi(X^{\varepsilon}((n-1)h,x))|x]$ is in
$\mathcal{C}^{3}_{P}$. We rewrite the global error, given by
$E^{\varepsilon}_{h}=\left|\mathbb{E}[\phi(X^{\varepsilon}(T,X_{0}))-\phi(X_{N}^{\varepsilon}(X_{0}))]\right|$,
with a telescopic sum,

$\displaystyle E^{\varepsilon}_{h}$
$\displaystyle\leq\sum_{n=1}^{N}\left|\mathbb{E}[\phi(X^{\varepsilon}(nh,X^{\varepsilon}_{N-n}(X_{0})))-\phi(X^{\varepsilon}((n-1)h,X^{\varepsilon}_{N-n+1}(X_{0})))]\right|$
$\displaystyle\leq\sum_{n=1}^{N}\left|\mathbb{E}[\phi_{n}(X^{\varepsilon}(h,X^{\varepsilon}_{N-n}(X_{0})))-\phi_{n}(X^{\varepsilon}_{1}(X^{\varepsilon}_{N-n}(X_{0})))]\right|$
$\displaystyle\leq\sum_{n=1}^{N}\mathbb{E}[\left|\mathbb{E}[\phi_{n}(X^{\varepsilon}(h,x))-\phi_{n}(X^{\varepsilon}_{1}(x))|x=X^{\varepsilon}_{N-n}(X_{0})]\right|]$
$\displaystyle\leq\sum_{n=1}^{N}C(1+\mathbb{E}[\left|X^{\varepsilon}_{N-n}(X_{0})\right|^{K}])h^{3/2}\leq
Ch^{1/2},$

where we used the bounded moments property of Proposition 3.6 and
$X_{n}^{\varepsilon}\in\mathcal{M}^{\varepsilon}_{h}$ (Lemma 3.1). ∎

###### Proof of Theorem 2.9.

We obtain straightforwardly from the expression of the integrator (2.2) that
$\left|\zeta(X_{n}^{\varepsilon})\right|\leq C_{h}\sqrt{\varepsilon}$. Using
this estimate and the notation introduced in the proof of Lemma 3.8, we
observe that

$\left|Y_{(i)}^{\varepsilon}(X_{n}^{\varepsilon})-Y_{(i)}^{0}(X_{n}^{\varepsilon})\right|\leq
C_{h}\sqrt{\varepsilon},\quad\left|\zeta_{(i)}^{\varepsilon}(X_{n}^{\varepsilon})\right|\leq
C_{h}\sqrt{\varepsilon},\quad i\in\\{1/2,1\\},$

where $Y_{(1/2)}^{0}(x)=\sigma\xi_{n}$ and $Y_{(1)}^{0}(x)=f(x)$. The Lagrange
multiplier given by (3.16) therefore satisfies
$\left|\lambda^{\varepsilon}_{n+1}-\widetilde{\lambda^{\varepsilon}_{n+1}}\right|\leq
C_{h}\sqrt{\varepsilon}$, where

$\widetilde{\lambda^{\varepsilon}_{n+1}}=-G_{X_{n+1}^{\varepsilon}-X_{n}^{\varepsilon}}^{-1}(X_{n}^{\varepsilon})\int_{0}^{1}g^{T}(X_{n}^{\varepsilon}+\tau(X_{n+1}^{\varepsilon}-X_{n}^{\varepsilon}))d\tau(\sigma\sqrt{h}\xi_{n}+hf(X_{n}^{\varepsilon})).$

Similarly to (3.16), the Lagrange multiplier $\lambda^{0}_{n+1}$ of the
constrained Euler integrator (1.6) is given by

$\lambda^{0}_{n+1}=-G_{X_{n+1}^{0}-X_{n}^{0}}^{-1}(X_{n}^{0})\int_{0}^{1}g^{T}(X_{n}^{0}+\tau(X_{n+1}^{0}-X_{n}^{0}))d\tau(\sigma\sqrt{h}\xi_{n}+hf(X_{n}^{0})).$

Using that $G_{y}^{-1}(x)$ is bounded, a straightforward calculation shows
that $G_{y}^{-1}(x)$ is Lipschitz in $x$, $y\in\mathbb{R}^{d}$, that is, there
exists a constant $L>0$ such that

$\left|G_{y_{1}}^{-1}(x_{1})-G_{y_{2}}^{-1}(x_{2})\right|\leq
L(\left|x_{1}-x_{2}\right|+\left|y_{1}-y_{2}\right|),\quad
x_{1},x_{2},y_{1},y_{2}\in\mathbb{R}^{d}.$

Thus, we get the estimate

$\left|\lambda^{\varepsilon}_{n+1}-\lambda^{0}_{n+1}\right|\leq
C_{h}(\sqrt{\varepsilon}+\left|X_{n}^{\varepsilon}-X_{n}^{0}\right|)+C\sqrt{h}\left|X_{n+1}^{\varepsilon}-X_{n+1}^{0}\right|,$

and, as $\lambda^{\varepsilon}_{n+1}$ and $\lambda^{0}_{n+1}$ are bounded
uniformly in $\varepsilon$, we have

$\left|g(X_{n}^{\varepsilon})\lambda^{\varepsilon}_{n+1}-g(X_{n}^{0})\lambda^{0}_{n+1}\right|\leq
C_{h}(\sqrt{\varepsilon}+\left|X_{n}^{\varepsilon}-X_{n}^{0}\right|)+C\sqrt{h}\left|X_{n+1}^{\varepsilon}-X_{n+1}^{0}\right|.$

From the definitions of $X_{n+1}^{\varepsilon}$ in (2.2) and $X_{n+1}^{0}$ in
(1.6), we deduce that

$\left|X_{n+1}^{\varepsilon}-X_{n+1}^{0}\right|\leq
C_{h}(\sqrt{\varepsilon}+\left|X_{n}^{\varepsilon}-X_{n}^{0}\right|)+C\sqrt{h}\left|X_{n+1}^{\varepsilon}-X_{n+1}^{0}\right|.$

If $h\leq h_{0}$ is small enough, we obtain

$\left|X_{n+1}^{\varepsilon}-X_{n+1}^{0}\right|\leq
C_{h}(\sqrt{\varepsilon}+\left|X_{n}^{\varepsilon}-X_{n}^{0}\right|).$

We deduce the estimate (2.6) by iterating this inequality for $n\leq N=T/h$. ∎

###### Proof of Theorem 2.12.

We take over the notations and the expression (3.16) of the Lagrange
multiplier $\lambda^{\varepsilon}_{1}$ that we used in the proof of Lemma 3.8.
Without loss of generality, we concentrate on the first step of the algorithm
with $X_{0}=x$. Replacing $\lambda^{\varepsilon}_{1}$ by the explicit formula
(3.16) in (2.2) yields that $X_{1}^{\varepsilon}$ is a fixed point of the
following map:

$\displaystyle F_{h}^{\varepsilon}(y)$
$\displaystyle=x+\sqrt{h}Y_{(1/2)}^{\varepsilon}(x)+hY_{(1)}^{\varepsilon}(x)+g(x)G_{y-x}^{-1}(x)\Big{[}\sqrt{h}\zeta_{(1/2)}^{\varepsilon}(x)+h\zeta_{(1)}^{\varepsilon}(x)$
$\displaystyle-\int_{0}^{1}g^{T}(x+\tau(y-x))d\tau(\sqrt{h}Y_{(1/2)}^{\varepsilon}(x)+hY_{(1)}^{\varepsilon}(x))\Big{]}.$

For $y_{1}$, $y_{2}\in\mathbb{R}^{d}$, Assumption 2.1 gives

$\displaystyle\left|F_{h}^{\varepsilon}(y_{2})-F_{h}^{\varepsilon}(y_{1})\right|$
$\displaystyle\leq C\sqrt{h}\left|y_{2}-y_{1}\right|,$

where we used that $G_{y}^{-1}(x)$ is Lipschitz in $x$, $y\in\mathbb{R}^{d}$.
We deduce that $F_{h}^{\varepsilon}$ is a uniform contraction for $h\leq
h_{0}$ small enough. ∎

##  4 Numerical experiments

In this section, we perform numerical experiments to confirm the theoretical
findings, on a torus in $\mathbb{R}^{3}$ and on the orthogonal group in high
dimension and codimension, in the spirit of the experiments in [46, 47, 30].

###  4.1 Uniform approximation for the invariant measure on a torus

We consider the example in codimension one of a torus in $\mathbb{R}^{3}$. We
apply the new method given by the discretization (2.4) for sampling the
invariant measure of (1.4) for different steps $h$ and parameters
$\varepsilon$, and we compare it with the Euler integrators (1.5) in
$\mathbb{R}^{d}$ and (1.6) on the manifold. We recall that for dynamics of the
form (1.4), the weak convergence implies the convergence for the invariant
measure (see [43]). The numerical experiments of this subsection hint that the
uniform accuracy property also extends to the approximation of the invariant
measure, but we leave the mathematical analysis for the invariant measure for
future work. We consider the constraint
$\zeta(x)=(x_{1}^{2}+x_{2}^{2}+x_{3}^{2}+R^{2}-r^{2})^{2}-4R^{2}(x_{1}^{2}+x_{2}^{2})$,
with $R=3$ and $r=1$, and we choose the map $f(x)=-25(x_{1}-R+r,x_{2},x_{3})$,
with $\sigma=\sqrt{2}$, the test function $\phi(x)=\left|x\right|^{2}$,
$M=10^{7}$ trajectories, the final time $T=10$, and the initial condition
$X_{0}=(R-r,0,0)$. Increasing the value of $T$ does not modify the computed
averages, which hints that we reached the equilibrium. The factor $25$ in $f$
confines the solution in a reasonably small neighborhood of the torus, which
allows a faster convergence to equilibrium and to take fewer trajectories. We
compute the Monte-Carlo estimator
$\widebar{J}=\frac{1}{M}\sum_{m=1}^{M}\phi(X_{N}^{(m)})\simeq\mathbb{E}[\phi(X_{N})]$,
where $X_{n}^{(m)}$ is the $m$-th realization of the integrator at time
$t_{n}=nh$, and $N$ is an integer satisfying $Nh=T$. We compare this
approximation with a reference value of $\int_{\mathbb{R}^{d}}\phi
d\mu_{\infty}^{\varepsilon}$ computed with the uniformly accurate integrator,
by using a timestep $h_{\text{ref}}=2^{-12}$. In the case of the constrained
Euler scheme (1.6), it amounts to comparing an approximation of
$\int_{\mathcal{M}}\phi d\mu_{\infty}^{0}$ with the reference value of
$\int_{\mathbb{R}^{d}}\phi d\mu_{\infty}^{\varepsilon}$. We observe in Figure
1 that the accuracy of the constrained integrator (1.6) for solving the
unconstrained problem (1.4) deteriorates when $\varepsilon$ grows larger, as
$\mu_{\infty}^{\varepsilon}$ deviates from $\mu_{\infty}^{0}$. The explicit
Euler scheme (1.5) faces stability issues when $\varepsilon\rightarrow 0$. The
accuracy of the new method for solving (1.4) does not deteriorate depending on
$\varepsilon$, and it shares a behavior similar to the constrained Euler
scheme (1.6) when $\varepsilon\rightarrow 0$, which is in agreement with
Theorem 2.9. The right graph of Figure 1 shows that the behavior of the error
in $\varepsilon$ is the same for any fixed value of $h$. This is a numerical
confirmation of the uniform accuracy property of the discretization (2.2) (in
the spirit of the numerical experiments in [11]), as stated in Theorem 2.7.
For any fixed $\varepsilon$, the right graph of Figure 1 also shows that the
error decreases when $h\rightarrow 0$. A plot of the error against $h$ for a
fixed $\varepsilon$ (not included for conciseness) shows a slope of order one
(see Remark 2.8).

Figure 1: Error for sampling the invariant measure of penalized Langevin
dynamics on a torus in $\mathbb{R}^{3}$ of the uniform discretization (2.2)
and the Euler integrators (1.5) and (1.6) for different values of
$\varepsilon$ with $h=2^{-9}T$ (left), and error curves versus $\varepsilon$
of the uniformly accurate method for different timesteps $h=2^{-i}T$ and
$i=6,\dots,10$ (right), with the final time $T=10$, the maps
$f(x)=-25(x_{1}-R+r,x_{2},x_{3})$, $\phi(x)=\left|x\right|^{2}$, and
$M=10^{7}$ trajectories.

###  4.2 Weak approximation on the orthogonal group

We apply the uniformly accurate method on a compact Lie group (in the spirit
of the numerical experiments in [46, 47]) to see how it performs in high
dimension and codimension. We choose the orthogonal group
$\operatorname{O}(m)=\\{M\in\mathbb{R}^{m\times m},M^{T}M=I_{m}\\}$, seen as a
submanifold of $\mathbb{R}^{m^{2}}$ of codimension $q=m(m+1)/2$. We compare
the explicit Euler scheme (1.5), the constrained Euler scheme (1.6), and the
new method on $\mathcal{M}=\operatorname{O}(m)$ for $m=2,\dots,5$ with the
parameters $\varepsilon=0.005$, $T=1$, and $h=2^{-7}$. Note that, as $h$ and
$\varepsilon$ share the same order of magnitude, the explicit Euler scheme
(1.5) can face stability issues, and the solution does not lie on the manifold
$\mathcal{M}$. Thus, we are in the regime where the convergence results for
both Euler schemes (1.5)-(1.6) do not apply. We choose $f=-\nabla V$, where
$V$ is given by

$V(x)=50\operatorname{Tr}((x-I_{m^{2}})^{T}(x-I_{m^{2}}))$ (4.1)

with the parameters $\sigma=\sqrt{2}$, $X_{0}=I_{d}$,
$\phi(x)=\operatorname{Tr}(x)$, and $M=10^{6}$ trajectories. The reference
solution for $J(m)=\mathbb{E}[\phi(X^{\varepsilon}(T))]$ is computed with the
uniformly accurate integrator with $h_{\text{ref}}=2^{-9}$. With the factor
$50$ in the potential (4.1), the trajectories stay close to $I_{m^{2}}$, and
$J(m)$ is close to $\phi(I_{m^{2}})=m$. This choice of factor permits one to
explore a reasonably small area of $\operatorname{O}(m)$, in order to avoid
zones close to $\mathcal{M}$ where the Gram matrix $G$ has a bad condition
number or is singular, and to reduce the number of trajectories needed. We
observe numerically that replacing the factor $50$ by $1$ in (4.1) induces a
severe timestep restriction. We present the results of the experiment in Table
1. We omit the results for the explicit Euler scheme (1.5) as the method is
inaccurate in this regime (error of size $1$). We observe that, in the regime
where $h$ and $\varepsilon$ share the same order of magnitude, the uniformly
accurate integrator performs significantly better than the Euler schemes (1.6)
and (1.5) for solving the problem (1.4). In this regime, the Euler method
(1.5) faces stability issues, and it is inappropriate to use the constrained
Euler scheme (1.6) as the solution $X^{\varepsilon}(t)$ of (1.4) is not close
to the solution $X^{0}(t)$ of (1.2). Moreover, the cost in time of the new
method stays the same in average for any value of $\varepsilon$ (results not
included for conciseness). This confirms numerically the uniform cost of
solving the fixed point problem (2.2), as stated in Theorem 2.12.

$m$ |  $\dim(\mathcal{M})$ |  $q$ | $J(m)$ |  $\widebar{J}_{\text{UA}}$ | error of $\widebar{J}_{\text{UA}}$ |  $\widebar{J}_{\text{EC}}$ | error of $\widebar{J}_{\text{EC}}$
---|---|---|---|---|---|---|---
$2$ | 1 | 3 |  $2.00934$ |  $2.00619$ |  $3.1\cdot 10^{-3}$ |  $1.99165$ |  $1.8\cdot 10^{-2}$
$3$ | 3 | 6 |  $3.01458$ |  $3.00821$ |  $6.4\cdot 10^{-3}$ |  $2.97460$ |  $4.0\cdot 10^{-2}$
$4$ | 6 | 10 |  $4.02050$ |  $4.00972$ |  $1.1\cdot 10^{-2}$ |  $3.94846$ |  $7.2\cdot 10^{-2}$
$5$ | 10 | 15 |  $5.02669$ |  $5.00842$ |  $1.8\cdot 10^{-2}$ |  $4.91298$ |  $1.1\cdot 10^{-1}$
Table 1: Numerical approximation of
$J(m)=\mathbb{E}[\phi(X^{\varepsilon}(T))]$ for $2\leq m\leq 5$ with the
estimator $\widebar{J}=M^{-1}\sum_{k=1}^{M}\phi(X_{N}^{(k)})$, where $(X_{n})$
is given by the uniform discretization (2.2) for $\widebar{J}_{\text{UA}}$ and
by the constrained Euler scheme (1.6) for $\widebar{J}_{\text{EC}}$ with their
respective errors. The average is taken over $M=10^{6}$ trajectories with the
potential (4.1), $\phi(x)=\operatorname{Tr}(x)$, the final time $T=1$, the
stiff parameter $\varepsilon=0.005$, and the timestep $h=2^{-7}$.

##  5 Conclusion and future work

In this work, we presented a new method for the weak numerical integration of
penalized Langevin dynamics evolving in the vicinity of manifolds of any
dimension and codimension. On the contrary of the other existing
discretizations, the accuracy of the proposed integrator is independent of the
size of the stiff parameter $\varepsilon$. Moreover, its cost does not depend
on $\varepsilon$, and it converges to the Euler scheme on the manifold when
$\varepsilon\rightarrow 0$. Throughout the analysis, we gave an expansion in
time of the solution to the penalized Langevin dynamics that is uniform in
$\varepsilon$, as well as new tools for the study of stochastic projection
methods for solving stiff SDEs.

Multiple questions arise from the work presented in this paper, with many of
great interest for physical applications. First, it would be interesting to
get convergence results with weaker assumptions or to develop uniformly
accurate integrators for different penalized dynamics with the same limit when
$\varepsilon\rightarrow 0$ such as the original penalized dynamics (1.3) (see
[16]). One could build integrators for penalized dynamics of the form

$dX^{\varepsilon}=f(X^{\varepsilon})dt+\sigma
dW+\frac{\sigma^{2}}{4}\nabla\ln(\det(G))-\frac{1}{\varepsilon}(gG^{-1}\zeta_{1})(X^{\varepsilon})dt-\frac{1}{\nu}(gG^{-1}\zeta_{2})(X^{\varepsilon})dt,$

where $\varepsilon$ and $\nu$ do not share the same order of magnitude, or for
constrained dynamics with a penalized term. One could also build a uniformly
accurate numerical scheme with high order in the weak context, or just in the
context of the invariant measure (in the spirit of the works [6, 31, 2, 3, 32,
28, 30] where numerical schemes of high order for the invariant measure and
weak order one were introduced). Postprocessors [45] proved to be an efficient
tool for reaching high order for the invariant measure without increasing the
cost of the method and could be used in this context. Moreover, the order
conditions presented in [28] for Runge-Kutta methods for solving Langevin
dynamics in $\mathbb{R}^{d}$ both in the weak sense and for the invariant
measure do not match with the order conditions for solving Langevin dynamics
constrained on the manifold $\mathcal{M}$, as presented in [30]. It would be
interesting to create a unified class of high order Runge-Kutta methods with
the same order conditions in $\mathbb{R}^{d}$, on the manifold $\mathcal{M}$
and in the vicinity of $\mathcal{M}$. The discretizations presented in this
paper could also be combined with Metropolis-Hastings rejection procedures
[39, 24], in the spirit of the works [20, 8, 35, 46, 36], in order to get an
exact approximation for the invariant measure with a rejection rate that does
not deteriorate in the regime $\varepsilon\rightarrow 0$.

Acknowledgements. The author would like to thank Gilles Vilmart for helpful
discussions. This work was partially supported by the Swiss National Science
Foundation, grants No. 200020_184614, No. 200021_162404 and No. 200020_178752.
The computations were performed at the University of Geneva on the Baobab
cluster using the Julia programming language.

## References

  * [1] A. Abdulle, D. Cohen, G. Vilmart, and K. C. Zygalakis.  High weak order methods for stochastic differential equations based on modified equations.  SIAM J. Sci. Comput., 34(3):A1800–A1823, 2012.
  * [2] A. Abdulle, G. Vilmart, and K. C. Zygalakis.  High order numerical approximation of the invariant measure of ergodic SDEs.  SIAM J. Numer. Anal., 52(4):1600–1622, 2014.
  * [3] A. Abdulle, G. Vilmart, and K. C. Zygalakis.  Long time accuracy of Lie-Trotter splitting methods for Langevin dynamics.  SIAM J. Numer. Anal., 53(1):1–16, 2015.
  * [4] I. Almuslimani, P. Chartier, M. Lemou, and F. Méhats.  Uniformly accurate schemes for oscillatory stochastic differential equations.  Submitted, 2021.
  * [5] H. C. Andersen.  Rattle: A “velocity” version of the Shake algorithm for molecular dynamics calculations.  J. Comput. Phys., 52(1):24–34, 1983.
  * [6] N. Bou-Rabee and H. Owhadi.  Long-run accuracy of variational integrators in the stochastic context.  SIAM J. Numer. Anal., 48(1):278–297, 2010.
  * [7] C.-E. Bréhier and S. Rakotonirina-Ricquebourg.  On asymptotic preserving schemes for a class of stochastic differential equations in averaging and diffusion approximation regimes.  Multiscale Model. Simul., 20(1):118–163, 2022.
  * [8] M. Brubaker, M. Salzmann, and R. Urtasun.  A family of MCMC methods on implicitly defined manifolds.  In Artificial intelligence and statistics, pages 161–172, 2012\.
  * [9] M. Calvo, L. O. Jay, J. I. Montijano, and L. Rández.  Approximate compositions of a near identity map by multi-revolution Runge-Kutta methods.  Numer. Math., 97(4):635–666, 2004.
  * [10] M. Calvo, J. I. Montijano, and L. Rández.  On explicit multi-revolution Runge-Kutta schemes.  Adv. Comput. Math., 26(1-3):105–120, 2007.
  * [11] P. Chartier, N. Crouseilles, M. Lemou, and F. Méhats.  Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations.  Numer. Math., 129(2):211–250, 2015.
  * [12] P. Chartier, M. Lemou, F. Méhats, and G. Vilmart.  A new class of uniformly accurate numerical schemes for highly oscillatory evolution equations.  Found. Comput. Math., 20(1):1–33, 2020.
  * [13] P. Chartier, M. Lemou, and L. Trémant.  Uniformly accurate numerical schemes for a class of dissipative systems.  Submitted, 2020.
  * [14] P. Chartier, J. Makazaga, A. Murua, and G. Vilmart.  Multi-revolution composition methods for highly oscillatory differential equations.  Numer. Math., 128(1):167–192, 2014.
  * [15] G. Ciccotti, R. Kapral, and E. Vanden-Eijnden.  Blue moon sampling, vectorial reaction coordinates, and unbiased constrained dynamics.  ChemPhysChem, 6(9):1809–1814, 2005.
  * [16] G. Ciccotti, T. Lelievre, and E. Vanden-Eijnden.  Projection of diffusions on submanifolds: application to mean force computation.  Comm. Pure Appl. Math., 61(3):371–408, 2008.
  * [17] D. Cohen and M. Sigg.  Convergence analysis of trigonometric methods for stiff second-order stochastic differential equations.  Numer. Math., 121(1):1–29, 2012.
  * [18] A. Debussche and E. Faou.  Weak backward error analysis for SDEs.  SIAM J. Numer. Anal., 50(3):1735–1752, 2012.
  * [19] I. I. Gikhman and A. V. Skorokhod.  Stochastic differential equations.  In The Theory of Stochastic Processes III, pages 113–219. Springer, 2007.
  * [20] M. Girolami and B. Calderhead.  Riemann manifold Langevin and Hamiltonian Monte Carlo methods.  J. R. Stat. Soc. Ser. B Stat. Methodol., 73(2):123–214, 2011.  With discussion and a reply by the authors.
  * [21] E. Hairer, C. Lubich, and G. Wanner.  Geometric numerical integration, volume 31 of Springer Series in Computational Mathematics.  Springer-Verlag, Berlin, second edition, 2006.  Structure-preserving algorithms for ordinary differential equations.
  * [22] E. Hairer and G. Wanner.  Solving ordinary differential equations. II, volume 14 of Springer Series in Computational Mathematics.  Springer-Verlag, Berlin, 2010.  Stiff and differential-algebraic problems, Second revised edition, paperback.
  * [23] R. Z. Hasminskii.  Stochastic stability of differential equations, volume 7 of Monographs and Textbooks on Mechanics of Solids and Fluids: Mechanics and Analysis.  Sijthoff & Noordhoff, Alphen aan den Rijn—Germantown, Md., 1980.  Translated from the Russian by D. Louvish.
  * [24] W. K. Hastings.  Monte Carlo sampling methods using Markov chains and their applications.  Biometrika, 57(1):97–109, 1970.
  * [25] S. Jin.  Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: a review.  Riv. Math. Univ. Parma (N.S.), 3(2):177–216, 2012.
  * [26] M. Kopec.  Weak backward error analysis for Langevin process.  BIT Numer. Math., 55(4):1057–1103, 2015.
  * [27] M. Kopec.  Weak backward error analysis for overdamped Langevin processes.  IMA J. Numer. Anal., 35(2):583–614, 2015.
  * [28] A. Laurent and G. Vilmart.  Exotic aromatic B-series for the study of long time integrators for a class of ergodic SDEs.  Math. Comp., 89(321):169–202, 2020.
  * [29] A. Laurent and G. Vilmart.  Multirevolution integrators for differential equations with fast stochastic oscillations.  SIAM J. Sci. Comput., 42(1):A115–A139, 2020.
  * [30] A. Laurent and G. Vilmart.  Order conditions for sampling the invariant measure of ergodic stochastic differential equations on manifolds.  Found. Comput. Math., 22(3):649–695, 2022.
  * [31] B. Leimkuhler and C. Matthews.  Rational construction of stochastic numerical methods for molecular sampling.  Appl. Math. Res. Express. AMRX, 2013(1):34–56, 2013.
  * [32] 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.
  * [33] T. Lelièvre, C. Le Bris, and E. Vanden-Eijnden.  Analyse de certains schémas de discrétisation pour des équations différentielles stochastiques contraintes.  C. R. Math. Acad. Sci. Paris, 346(7-8):471–476, 2008.
  * [34] T. Lelièvre, M. Rousset, and G. Stoltz.  Free energy computations.  Imperial College Press, London, 2010.  A mathematical perspective.
  * [35] T. Lelièvre, M. Rousset, and G. Stoltz.  Langevin dynamics with constraints and computation of free energy differences.  Math. Comp., 81(280):2071–2125, 2012.
  * [36] T. Lelièvre, M. Rousset, and G. Stoltz.  Hybrid Monte Carlo methods for sampling probability measures on submanifolds.  Numer. Math., 143(2):379–421, 2019.
  * [37] J. C. Mattingly, A. M. Stuart, and D. J. Higham.  Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise.  Stochastic Process. Appl., 101(2):185–232, 2002.
  * [38] B. Melendo and M. Palacios.  A new approach to the construction of multirevolution methods and their implementation.  Appl. Numer. Math., 23(2):259–274, 1997.
  * [39] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller.  Equation of state calculations by fast computing machines.  J. Chem. Phys., 21(6):1087–1092, 1953.
  * [40] G. N. Milstein and M. V. Tretyakov.  Stochastic numerics for mathematical physics.  Scientific Computation. Springer-Verlag, Berlin, 2004.
  * [41] P. Plecháč and M. Rousset.  Implicit mass-matrix penalization of Hamiltonian dynamics with application to exact sampling of stiff systems.  Multiscale Model. Simul., 8(2):498–539, 2009/10.
  * [42] J.-P. Ryckaert, G. Ciccotti, and H. J. Berendsen.  Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes.  J. Comput. Phys., 23(3):327–341, 1977.
  * [43] 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.
  * [44] G. Vilmart.  Weak second order multirevolution composition methods for highly oscillatory stochastic differential equations with additive or multiplicative noise.  SIAM J. Sci. Comput., 36(4):A1770–A1796, 2014.
  * [45] G. Vilmart.  Postprocessed integrators for the high order integration of ergodic SDEs.  SIAM J. Sci. Comput., 37(1):A201–A220, 2015.
  * [46] E. Zappa, M. Holmes-Cerfon, and J. Goodman.  Monte Carlo on manifolds: sampling densities and integrating functions.  Comm. Pure Appl. Math., 71(12):2609–2647, 2018.
  * [47] W. Zhang.  Ergodic SDEs on submanifolds and related numerical sampling schemes.  ESAIM Math. Model. Numer. Anal., 2019.
  * [48] K. C. Zygalakis.  On the existence and the applications of modified equations for stochastic differential equations.  SIAM J. Sci. Comput., 33(1):102–130, 2011.

##  Appendix A Proof of Theorem 2.3

In this section, we prove the convergence of the penalized Langevin dynamics
(1.4) to the constrained dynamics (1.2) when $\varepsilon\rightarrow 0$, as
stated in Theorem 2.3. The proof uses techniques and arguments similar to
those in [16, Appx. C]. However, since we rescaled the stiff term in (1.4),
there is no need for a change of time to prove the convergence to the
constrained dynamics.

###### Proof of Theorem 2.3.

In the orthogonal coordinates system given by Assumption 2.1, equation (1.4)
becomes

$\displaystyle d\varphi(X^{\varepsilon}(t))$
$\displaystyle=(\varphi^{\prime}f)(X^{\varepsilon}(t))dt+\sigma\varphi^{\prime}(X^{\varepsilon}(t))dW(t)$
$\displaystyle+\frac{\sigma^{2}}{4}(\varphi^{\prime}\nabla\ln(\det(G)))(X^{\varepsilon}(t))dt+\frac{\sigma^{2}}{2}\Delta\varphi(X^{\varepsilon}(t))dt,$
$\displaystyle d\zeta(X^{\varepsilon}(t))$
$\displaystyle=(g^{T}f)(X^{\varepsilon}(t))dt+\sigma
g^{T}(X^{\varepsilon}(t))dW(t)+\frac{\sigma^{2}}{4}(g^{T}\nabla\ln(\det(G)))(X^{\varepsilon}(t))dt$
$\displaystyle+\frac{\sigma^{2}}{2}\operatorname{div}(g)(X^{\varepsilon}(t))dt-\frac{1}{\varepsilon}\zeta(X^{\varepsilon}(t))dt,$

where we used that $\varphi^{\prime}g=0$ and that $\zeta^{\prime}=g^{T}$.
Therefore, $\zeta(X^{\varepsilon}(t))$ satisfies

$\displaystyle\zeta(X^{\varepsilon}(t))$
$\displaystyle=\sigma\int_{0}^{t}e^{(s-t)/\varepsilon}g^{T}(X^{\varepsilon}(s))dW$
$\displaystyle+\int_{0}^{t}e^{(s-t)/\varepsilon}\Big{[}g^{T}f+\frac{\sigma^{2}}{4}g^{T}\nabla\ln(\det(G))+\frac{\sigma^{2}}{2}\operatorname{div}(g)\Big{]}(X^{\varepsilon}(s))dt,$

and a bound on $\zeta(X^{\varepsilon}(t))$ is given by

$\displaystyle\mathbb{E}[\left|\zeta(X^{\varepsilon}(t))\right|^{2}]$
$\displaystyle\leq
C\int_{0}^{t}e^{2(s-t)/\varepsilon}dt+C\Big{(}\int_{0}^{t}e^{(s-t)/\varepsilon}dt\Big{)}^{2}\leq
C\varepsilon.$ (A.1)

On the other hand, we rewrite (1.2) in the orthogonal coordinates as

$\displaystyle d\varphi(X^{0}(t))$
$\displaystyle=(\varphi^{\prime}f)(X^{0}(t))dt+\sigma\varphi^{\prime}(X^{0}(t))dW$
$\displaystyle+\frac{\sigma^{2}}{4}(\varphi^{\prime}\nabla\ln(\det(G)))(X^{0}(t))dt+\frac{\sigma^{2}}{2}\Delta\varphi(X^{0}(t))dt,$
$\displaystyle d\zeta(X^{0}(t))$ $\displaystyle=0,$

where we used that

$\frac{1}{2}\varphi^{\prime}\nabla\ln(\det(G))+\Delta\varphi=\sum_{i=1}^{d}\varphi^{\prime}(\Pi_{\mathcal{M}}^{\prime}(\Pi_{\mathcal{M}}e_{i})e_{i})+\sum_{i=1}^{d}\varphi^{\prime\prime}(\Pi_{\mathcal{M}}e_{i},\Pi_{\mathcal{M}}e_{i}).$

With Assumption 2.1, we obtain

$\mathbb{E}[\left|\varphi(X^{\varepsilon}(t))-\varphi(X^{0}(t))\right|^{2}]\leq
C\int_{0}^{t}\mathbb{E}[\left|X^{\varepsilon}(s)-X^{0}(s)\right|^{2}]ds,$

and from (A.1) we deduce that

$\mathbb{E}[\left|\zeta(X^{\varepsilon}(t))-\zeta(X^{0}(t))\right|^{2}]\leq
C\varepsilon.$

Thus, $\psi(X^{\varepsilon}(t))$ satisfies

$\mathbb{E}[\left|\psi(X^{\varepsilon}(t))-\psi(X^{0}(t))\right|^{2}]\leq
C\int_{0}^{t}\mathbb{E}[\left|X^{\varepsilon}(s)-X^{0}(s)\right|^{2}]ds+C\varepsilon,$

and, as $\psi^{-1}$ is Lipschitz, we have

$\mathbb{E}[\left|X^{\varepsilon}(t)-X^{0}(t)\right|^{2}]\leq
C\int_{0}^{t}\mathbb{E}[\left|X^{\varepsilon}(s)-X^{0}(s)\right|^{2}]ds+C\varepsilon,$

which, with the use of the Gronwall lemma, gives the desired estimate. ∎