Ergodic SDEs on submanifolds and related numerical sampling schemes
Wei Zhang





[1]Zuse Institute Berlin, Takustrasse 7, 14195 Berlin,
Germany.  Email: wei.zhang@fu-berlin.de

Ergodic SDEs on submanifolds and related numerical sampling schemes
    
  Wei Zhang 1 

    
===================================================================


  In many applications, it is often necessary to 
   sample the mean value of certain quantity with respect to a probability
   measure μ on the level set of a smooth function ξ:ℝ^d
   →ℝ^k, 1 ≤ k < d.
  A specially interesting case is the so-called conditional probability measure, which is
  useful in the study of free energy calculation and model reduction of diffusion processes. 
  By Birkhoff's ergodic theorem, one approach to estimate the mean value is to
  compute the time average along an infinitely long trajectory of an ergodic diffusion process on the level set whose invariant
  measure is μ. 
Motivated by the previous work of <cit.>, as well as
the work of <cit.>, in this paper we construct a family of ergodic diffusion processes on 
  the level set of ξ whose invariant measures coincide with the given one.
  For the conditional measure,  we propose a consistent numerical scheme which samples the conditional measure asymptotically.
  The numerical scheme doesn't require computing the second derivatives of ξ and the
  error estimates of its long time sampling efficiency are obtained.






§ INTRODUCTION


Many stochastic dynamical systems in real-world applications in
physics, chemistry, and biology often involve a large number of
degrees of freedom which evolve on vastly different time scales. Understanding
the behavior of these systems can be highly challenging due to the
high dimensionality and the existence of multiple time scales. 
To tackle these difficulties, the terminology reaction coordinate, or
collective variable, is often introduced
to help describe the essential dynamical behavior of complex systems <cit.>.

In various research topics, in particular those related to molecular dynamics, 
one often encounters the problem of computing the mean value of certain quantity on the level set 

    Σ=ξ^-1(0)={ x ∈ℝ^d  | ξ(x) = 0∈ℝ^k}

of a reaction coordinate function ξ : ℝ^d →ℝ^k, 1 ≤ k <d.  
Among different probability measures on Σ, the one defined by

    dμ_1 = 1/Z e^-β U[(∇ξ^T
    ∇ξ)]^-1/2 dν

is especially relevant in applications and is sometimes called the conditional probability measure
on Σ. In (<ref>), the parameter β > 0, U: ℝ^d →ℝ is a smooth function, Z is the normalization constant, 
∇ξ denotes the d × k Jacobian matrix of the map ξ, 
and ν is the surface measure on Σ induced from the Lebesgue measure on ℝ^d.  
The probability measure μ_1 has a probabilistic interpretation, and the numerical computation of the mean value 

    f = ∫_Σ f(x)  dμ_1(x)

for a function f on the level set is involved in various contexts, such as 
free energy calculations based on the thermodynamics integration
formula <cit.>.

Applying Birkhoff's ergodic
theorem, the mean value f can be approximated by the time average 
1/T∫_0^T f(X_s) ds
along a long trajectory of the process X_s which evolves on the level set
Σ and has the invariant measure μ_1. 
For this purpose, it is helpful to construct a diffusion process on the level set
with the correct invariant measure μ_1, i.e., to write down
the stochastic differential equation (SDE) of X_s in ℝ^d. 
While finding such a SDE is trivial in the linear reaction
coordinate case <cit.>, it is not obvious when the reaction
coordinate ξ is a nonlinear function of system's state. 

In the literature, the problem finding SDEs on the level set of the reaction coordinate
function with a given invariant measure has been considered in the 
study of free energy calculations <cit.>. Given a smooth function
U : ℝ^d →ℝ, 
the authors in <cit.> constructed a diffusion process
Y_s on Σ whose unique invariant measure is μ_2, given by 

    dμ_2 = 1/Z e^-β U dν .

It is also shown in <cit.> that this process Y_s can be obtained by projecting 
the dynamics 

    d Y_s     = -∇ U (Y_s)  ds +
      √(2β^-1) dW_s

from ℝ^d onto the level set Σ, where
 W_s = (W_s^1, ⋯, W_s^d)^T is a d-dimensional Brownian motion.
 The dynamics Y_s can be used to sample μ_2, and therefore to sample the
 conditional measure μ_1 in (<ref>) as well, by either modifying
 the potential U or reweighting the
function f according to the factor [(∇ξ^T ∇ξ)]^-1/2.
In a more recent work <cit.>, 
the authors studied the constrained Langevin dynamics, which evolves on the submanifold of the entire phase space 
including both position and momentum.
It is shown in <cit.> that 
 the position components of the constrained
Langevin dynamics has the marginal invariant measure which coincides with μ_2.
Therefore, it can also be used to compute the average f with
respect to the conditional measure μ_1 (by either modifying the potential
or reweighting f according to 
[(∇ξ^T ∇ξ)]^-1/2).
Detailed studies on the numerical schemes as well as applications of the constrained Langevin dynamics have been carried out in <cit.>. 

The same conditional probability measure μ_1 in (<ref>), as well as 
the average f in (<ref>), also plays an important role
in the study of the effective dynamics of diffusion processes
<cit.>.
As a generalization of the dynamics (<ref>), 
 the diffusion process 

    d Y^i_s     = -(a_ij∂ U/∂ x_j)
      (Y_s) ds +
      1/β∂ a_ij/∂ x_j(Y_s)  ds +
      √(2β^-1)σ_ij(Y_s)  dW^j_s ,    1 ≤ i ≤ d ,

and its effective dynamics have been considered in <cit.>, where the matrix-valued coefficients σ, a : ℝ^d
→ℝ^d× d are related by a = σσ^T, such
that a is uniformly positive definite. Notice that, 
(<ref>) is written in component-wise form with Einstein's summation
convention (the same Einstein's summation convention will be used throughout this paper, whenever no ambiguity will arise), and it reduces to (<ref>) when σ=a =.
The infinitesimal generator of (<ref>) can be written as 

    ℒ = e^β U/β∂/∂ x_i(e^-β Ua_ij∂/∂ x_j) .

Under mild conditions on U, it is known that,
for any (smooth, uniformly positive definite) coefficient a, the dynamics
(<ref>) has the common unique invariant measure whose probability
density is 1/Ze^-β U with respect to the Lebesgue measure on
ℝ^d.

Motivated by these previous work, in this paper we try to answer the following
two questions.




(Q1)  Besides the process constructed in
  <cit.> that is closely related to (<ref>),
  can we obtain other diffusion processes  on Σ, which are probably
  related to (<ref>) involving the coefficients σ, a, and have the same invariant measure? 
In particular, can we construct SDEs on Σ whose invariant measure is
μ_1? 




(Q2)  Numerically, instead of sampling μ_2,
can we directly estimate the mean value in (<ref>) with respect to μ_1, preferably with a numerical algorithm that is easy to implement? 



The main contributions of the current work are related to the above questions
and are summarized below. First, concerning Question (Q1), 
in Theorem <ref> of Section <ref>, we will construct a family of diffusion processes on Σ which sample
either μ_1 or μ_2. In particular, we show that the diffusion process

    dX_s^i 
    =     - (Pa)_ij∂ U/∂ x_j ds + 1/β∂
    (Pa)_ij/∂ x_j ds + √(2β^-1)  P_j,i  dW_s^j ,    1 ≤ i ≤ d ,

evolves on Σ and the invariant measure is the conditional probability
measure μ_1 in (<ref>), where the projection map P and the invertible k × k symmetric matrix Ψ
are given by 

    P=- a∇ξΨ^-1∇ξ^T,  Ψ = ∇ξ^T a ∇ξ.

Correspondingly, the infinitesimal generator of (<ref>) is 

    ℒ = e^β U/β∂/∂ x_i(e^-β
      U(Pa)_ij∂/∂ x_j) ,

which should be compared to the infinitesimal generator in (<ref>).
Second, concerning Question (Q2), in Section <ref> we study a numerical algorithm which
estimates the mean value f in (<ref>). Specifically,
we propose to use the numerical scheme

    x^(l + 1/2)_i =    x^(l)_i + (-a_ij∂ U/∂
      x_j + 1/β∂ a_ij/∂ x_j)(x^(l))  h
      + √(2 β^-1h) σ_ij(x^(l)) η^(l)_j ,    1 ≤ i ≤
      d  ,
    
      x^(l+1) =   Θ(x^(l+1/2)) ,

with x^(0)∈Σ, and to approximate f by f_n = 1/n∑_l=0^n-1 f(x^(l)).
In (<ref>), h is the step-size, 
η^(l)=(η_1^(l), η_2^(l), ⋯,
η_d^(l))^T are independent d-dimensional standard Gaussian random
variables, and 
Θ(x) = lim_s→ +∞φ(x,s) is the limit of the flow map 

    dφ(x,s)/ds =    - (a∇ F) (φ(x,s)) ,  φ(x,0) = x,     ∀ x∈ℝ^d ,

with F(x)=1/2|ξ(x)|^2=1/2∑_α=1^kξ_α^2(x).
Following the approach developed in <cit.>,
in Theorem <ref>, we obtain the
estimates of the approximation error between f_n and
f. While different
constraint approaches have been proposed in the
literature <cit.>, to the best of the
author's knowledge, constraint using the flow map φ has not been studied yet. 

Let us comment on the two contributions mentioned above.
First, knowing the SDE (<ref>) and the expression
(<ref>) of its infinitesimal generator ℒ is helpful for analysis. In fact, in Section <ref>, 
the analysis of sampling error estimate of the scheme (<ref>) 
relies on Poisson equation on Σ related to ℒ in (<ref>). 
Furthermore, (<ref>) plays a role in the work <cit.> in analyzing the approximation quality of the effective dynamics, while SDE (<ref>)
has been used in <cit.> to study fluctuation relations and Jarzynski's equality for nonequilibrium systems.
Second, we emphasize that Θ(x) in the scheme (<ref>) can be evaluated by solving the ODE (<ref>) starting from x. 
Although Θ is defined as the limit when s → +∞,
in many cases the computational cost is not large, due to the exponential
convergence of the (gradient) flow (<ref>) to its limit,
particularly for the initial state x=x^(l+1/2) that is close to Σ. 
Furthermore, comparing to the direct (Euler-Maruyama) discretization of SDE
(<ref>) which may deviate from Σ and
requires second order derivatives of ξ, the scheme (<ref>)
satisfies x^(l)∈Σ for all l ≥ 0, and 
it doesn't require computing the second order derivatives of ξ. 
Therefore, we expect the numerical scheme
(<ref>)–(<ref>) is both stable and relatively easy to implement.
Readers are referred to Remark <ref>–<ref>
in Section <ref> and Example 1 in Section <ref> for further algorithmic discussions. 

In the following, we briefly explain the approach that we will use to study Question (Q1), as well as the idea behind the scheme (<ref>)–(<ref>).
Concerning Question (Q1), we take the manifold point of view by considering
ℝ^d as a Riemannian manifold ℳ=(ℝ^d, g) with the metric g=a^-1, defined by 

    g(u, v) = ⟨u, v⟩_g = u_i (a^-1)_ij
      v_j ,   ∀ u, v∈ℝ^d .

A useful observation is that, for ℒ in (<ref>), we have <cit.>

    ℒ f    = [-^ℳ(U + 1/2βln G) + 1/βΔ^ℳ] f ,
        ∀   f : ℝ^d →ℝ ,

where G=g, and ^ℳ, Δ^ℳ denote the gradient and the Laplacian-Beltrami operator on ℳ, respectively. 
Accordingly, (<ref>) can be written as a SDE on ℳ as

    dY_s = -^ℳ(U + 1/2βln G)  ds + √(2β^-1)
      dB_s ,

where B_s is the Brownian motion on ℳ <cit.>. 
Conversely, SDE (<ref>) can be seen as the equation of (<ref>)
under the (global) coordinate chart of ℳ. This equivalence allows us to study (<ref>) on ℝ^d by
the corresponding SDE (<ref>) on manifold ℳ. 
Comparing to (<ref>), one advantage of working with
the abstract equation (<ref>) is that the invariant measure
of (<ref>) can be recognized as
easily as in (<ref>), provided that we apply integration by parts
formula on the manifold ℳ. 

A family of ergodic SDEs on Σ (i.e., Question (Q1)) is obtained by taking the same manifold
point of view. Specifically, consider Σ as a submanifold of ℳ and
denote by ^Σ, Δ^Σ, B_s the gradient operator, the Laplacian and the Brownian motion 
(with generator 1/2Δ^Σ <cit.>) on Σ, respectively.
Since the infinitesimal generator of the SDE 

    dY_s = - ^Σ U ds + √(2 β^-1) dB_s

is ℒ = -^Σ U + 1/βΔ^Σ, 
under mild assumptions on U, it is straightforward to verify that dynamics (<ref>) 
evolves on Σ and has the unique invariant measure 1/Ze^-β Udν_g,
where ν_g is the surface measure on Σ induced from the metric g=a^-1 on ℳ=(ℝ^d, g).
 Therefore, answering Question (Q1) boils down to calculating the
 expression of (<ref>) under the coordinate
 chart of ℳ (not Σ). 
 This will be achieved by calculating the expressions
 of ^Σ, Δ^Σ under the coordinate chart of
 ℳ and then figuring out the
 relation between the two measures ν and ν_g.

Concerning the idea behind the numerical scheme (<ref>)–(<ref>),
we recall that one way to (approximately) sample μ_1 on Σ is to constrain
the dynamics (<ref>) in the neighborhood of Σ by adding an extra potential to it. 
This is often termed as softly constrained dynamics <cit.>
and has been widely used in applications.
In this context, one consider the dynamics

    d X^ϵ,i_s     = [- a_ij∂ U/∂
    x_j  -
    1/ϵa_ij∂/∂ x_j(1/2∑_α=1^k
      ξ_α^2) +
    1/β∂ a_ij/∂ x_j]  ds +
      √(2β^-1)σ_ij  dW^j_s ,

where ϵ > 0, 1 ≤ i ≤ d, based on the fact that the invariant measure
of (<ref>) converges weakly to μ_1, as ϵ→ 0.
The dynamics (<ref>) stays close to Σ most of the time, thanks to the existence of the extra constraint force.
Furthermore, only the first order
derivatives of ξ are involved. In spite of these nice properties, however, direct simulation of
(<ref>) is inefficient when ϵ is small, 
because the time step-size in numerical simulations becomes severely
limited due to the strong stiffness in the dynamics. 
Indeed, our numerical scheme is motivated in order to overcome the aforementioned drawback of the softly constrained dynamics (<ref>),  
and the scheme (<ref>)–(<ref>) can be viewed as a multiscale numerical method for
(<ref>), where the stiff and non-stiff terms in 
(<ref>) are handled separately <cit.>. 
In contrast to the previous work <cit.>, where the
convergence of (<ref>) was studied on a finite time interval, our result concerns the long time sampling efficiency of the discretized numerical scheme.

Before concluding this introduction, we compare the current work with several previous ones. Generally speaking, Monte Carlo samplers (based on
ergodicity) either on ℝ^d or on its submanifolds can be
classified into Metropolis-adjusted samplers and samplers without Metropolis step (unadjusted).
For Metropolis-adjusted methods, in particular, exploiting
Riemannian geometry structure to develop MCMC methods has been studied in <cit.>.
The authors there demonstrated that incorporating the geometry of the space into numerical methods can lead
to significant improvement of the sampling efficiency. 
In line with this development, in Section <ref> we will consider a concrete example
where a non-constant matrix a can help remove the stiffness in the sampling task.
On the other hand, despite of the common Riemannian manifold point of view in the current work and in <cit.>, 
the main difference is that the current work deals with sampling on the
submanifold Σ instead of the entire ℝ^d (or its domain).
The derivations in the current work are more involved mainly due to this difference.
Besides sampling on the entire space, Metropolis-adjusted samplers on submanifolds,
using either MCMC or Hybrid Monte Carlo, have been considered in several recent
work <cit.>. Reversible Metropolis random walk on submanifolds has been
constructed in <cit.>, which is then extended in <cit.> by allowing non-zero gradient forces in the proposal move. 
In contrast to these Metropolis-adjusted samplers, the numerical scheme
(<ref>)–(<ref>) in the current work is unadjusted (without Metropolis-step) and samples the conditional probability measure μ_1
when the step-size h→ 0. This means that in practice the step-size h should be chosen properly such that the discretization error is tolerable.
In this direction, we point out that unadjusted samplers on ℝ^d,
which naturally arise from discretizations of SDEs, have been well studied in
the
literature <cit.>.
The current work can be thought as a further step along this direction for sampling schemes on submanifolds,
by applying the machinery developed in <cit.>. 
Comparison between the scheme (<ref>) and the Metropolis-adjusted algorithm
in <cit.> can be found in
Remark <ref>, as well as in Example 2 in Section <ref>.
We also refer to <cit.> for related discussions.

The rest of the paper is organized as follows. In Section <ref>, we construct ergodic SDEs on Σ which
sample either μ_1 or μ_2. In Section <ref>, we study the numerical scheme
(<ref>)–(<ref>) and quantify its approximation error in estimating the mean value in
(<ref>). In Section <ref>, we demonstrate our results through concrete examples.
Conclusions and further discussions are made in Section <ref>.
Technical details related to the Riemannian manifold ℳ in Section <ref> are included in Appendix <ref>.
Proofs of the results in Section <ref> are collected in Appendix <ref>.

Finally, we conclude this introduction with the assumptions which will be
made (implicitly) throughout this paper.

  The matrix σ: ℝ^d →ℝ^d × d is
  both smooth and invertible at each x ∈ℝ^d. The matrix
  a=σσ^T is uniformly positive definite with uniformly bounded
  inverse a^-1. 
 


  The function ξ: ℝ^d →ℝ^k is C^2 smooth and the level set Σ is both connected and compact, such that 
  (∇ξ) = k at each x ∈Σ. 
  



§ SDES OF ERGODIC DIFFUSION PROCESSES ON Σ


In this section, we construct SDEs of ergodic processes on Σ that sample a given invariant measure. The main result of this section is Theorem <ref>, which shows that the invariant measure of the SDE (<ref>) in Introduction is the conditional probability measure
μ_1 in (<ref>). Readers who are mainly interested in numerical
algorithms can read Theorem <ref> and then directly jump to Section <ref>. 

First of all, let us point out that, the semigroup approach based on
functional inequalities on Riemannian manifolds is well developed 
to study the solution of Fokker-Planck equation towards equilibrium.  
One sufficient condition for the exponential convergence of the
Fokker-Planck equation (and therefore the ergodicity of the corresponding dynamics) is the famous Bakry-Emery criterion <cit.>. In particular, concrete conditions
are given in <cit.> which guarantee the exponential convergence to
the unique invariant measure. In the following, we will always assume that the
potential U ∈ C^∞(Σ) and the Bakry-Emery condition in <cit.> is satisfied.

Recall that ℳ=(ℝ^d, g), where g=a^-1 and ν_g is the surface measure on Σ induced from
ℳ. Matrices P, P_j,i are given in (<ref>), (<ref>) in Appendix <ref>, respectively.
For 1 ≤ i ≤ d, e_i denotes the vector whose ith component equals to 1 while all the other d-1 components equal to 0.
We refer the reader to Appendix <ref> for further details.
Let us first consider the 
probability measure μ on Σ given by dμ = 1/Z e^-β
U dν_g, where β > 0 and Z is a normalization constant. 
The following proposition is a direct application of
Proposition <ref> in Appendix <ref>.

  Consider the dynamics on ℝ^d which satisfies the Ito SDE 

    dY_s^i =    - (Pa)_ij∂[U - 1/2βln((
        a)^-1(∇ξ^T a ∇ξ))]/∂ x_j ds + 1/β∂
    (Pa)_ij/∂ x_j ds + √(2β^-1)  P_j,i  dW_s^j

for 1 ≤ i ≤ d, where W_s = (W_s^1, W_s^2, ⋯, W_s^d)^T is a d-dimensional Brownian motion.
Suppose Y_0 ∈Σ, then Y_s ∈Σ almost surely for s≥ 0. 
Furthermore, it has a unique invariant measure μ given by dμ =
1/Z e^-β U dν_g.



  Using (<ref>) and Proposition <ref> in Appendix <ref>, we know that the infinitesimal generator 
  of SDE (<ref>) is 
  
    ℒf = -⟨^Σ U, ^Σ f⟩_g + 1/βΔ^Σ f ,   ∀ f: Σ→ℝ ,

    where ^Σ, Δ^Σ are the gradient and Laplace-Beltrami operators on Σ of ℳ, respectively.
    Applying Ito's formula to ξ_α(Y_s), we have 
    
    dξ_α(Y_s) = ℒξ_α(Y_s) ds + √(2β^-1)∂ξ_α(Y_s)/∂ x_i P_j, i dW_s^j ,    1 ≤α≤ k .

    Using (<ref>), (<ref>) in Appendix <ref>, and the fact that ^Σξ_α =
    P^ℳξ_α=0, it is straightforward to verify that
    
    ℒξ_α = -⟨^Σ U, ^Σξ_α⟩_g + 1/β^Σ(^Σξ_α) = 0  ,
       ∂ξ_α/∂ x_i P_j, i = 0  ,
             1 ≤ j ≤ d ,

    on Σ, which implies dξ_α(Y_s) = 0, ∀ s ≥ 0.
    Since Y_0 ∈Σ, we conclude that ξ_α(Y_s) = ξ_α(Y_0) = 0 
    a.s. s ≥ 0, for 1 ≤α≤ k, and
    therefore Y_s ∈Σ for  s ≥ 0, almost surely. 

    Using the expression (<ref>) of ℒ and the integration
    by parts formula (<ref>) in Appendix <ref>, it is easy to see that μ is an invariant
    measure of the dynamics (<ref>).  The uniqueness is implied by the
    exponential convergence result established in <cit.>, since we assume Bakry-Emery condition is satisfied.


In the above, we have considered the level set Σ as a submanifold of ℳ=(ℝ^d, g). 
In applications, on the other hand, it is natural to view Σ as a submanifold of the
standard Euclidean space ℝ^d, with the surface measure ν on Σ that is induced
from the Euclidean metric on ℝ^d. 
As already mentioned in the Introduction, the following two probability measures 

    dμ_1 = 1/Z e^-β U[(∇ξ^T ∇ξ)]^-1/2 dν ,     dμ_2 = 1/Z e^-β U dν,

where Z denotes possibly different normalization constants,
are often interesting and arise in many
situations <cit.>.
In particular, μ_1 has a probabilistic interpretation and 
often appears in the study of free energy calculation and model reduction of
stochastic dynamics <cit.>. 
In order to construct processes which sample μ_1 or μ_2, we need to 
figure out the relations between the two surface measures ν_g and ν on
Σ. 

  Let ν_g, ν be the surface measures on Σ induced
  from the metric g=a^-1 and the Euclidean metric on ℝ^d,
  respectively. We have 
  
    dν_g =
        ( a)^-1/2[(∇ξ^T
        a∇ξ)/(∇ξ^T∇ξ)]^1/2
    dν  .




  Let x ∈Σ and v_1, v_2, ⋯, v_d-k be a basis
  of T_xΣ. 
  Assume that v_i = c_ije_j, where
  c = (c_ij) is a
  (d-k) × d matrix whose rank is d-k.
  Using the fact ⟨v_i, ^ℳξ_α⟩_g = 0 for 1≤ i ≤ d-k, 1 ≤α≤ k, we can deduce that
   c ∇ξ = 0. 
  Calculating the surface measures ν_g and ν under this basis, we
  obtain
  
    dν_g =
        [(ca^-1c^T)/(cc^T)]^1/2
        dν .

  To simplify the right hand side of (<ref>), we use the following equality 
  
    [       c; ∇ξ^T  a ][ c^T  ∇ξ ]
        = 
        [       cc^T          0; ∇ξ^T a c^T  ∇ξ^T a ∇ξ ]
        = 
        [ ca^-1;  ∇ξ^T ]
        a
        [ c^T  ∇ξ ] .

  After computing the determinants of the last two matrices above, we obtain
  
    (cc^T)  (∇ξ^T a ∇ξ) 
        =
        ( a) [
            [ ca^-1;  ∇ξ^T ][ c^T  ∇ξ ]]  
        = ( a)  [ ca^-1c^T  ca^-1∇ξ;        0   ∇ξ^T∇ξ ] .

  The conclusion follows after we substitute the above relation into
  (<ref>).


Applying Lemma <ref> and
Proposition <ref>, we can obtain ergodic processes whose invariant
measures are given in (<ref>). 

  Let μ_1, μ_2 be the two probability measures on Σ
  defined in (<ref>).
  Consider the dynamics X_s, Y_s on ℝ^d which satisfy the Ito SDEs 

    dX_s^i 
    =     - (Pa)_ij∂ U/∂ x_j ds + 1/β∂
    (Pa)_ij/∂ x_j ds + √(2β^-1)  P_j,i  dW_s^j ,

and

    dY_s^i =     - (Pa)_ij∂[U - 1/2βln(∇ξ^T ∇ξ)]/∂ x_j ds + 1/β∂
    (Pa)_ij/∂ x_j ds + √(2β^-1)  P_j,i  dW_s^j ,

for 1 ≤ i ≤ d, where β > 0 and W_s = (W_s^1, W_s^2,⋯, W_s^d)^T is a
d-dimensional Brownian motion. Suppose that X_0 , Y_0 ∈Σ, then X_s,
Y_s ∈Σ almost surely for s≥ 0. 
Furthermore, the unique invariant probability measures of the dynamics X_s and Y_s 
are μ_1 and μ_2, respectively.



  Applying Lemma <ref>, we can rewrite the probability measures
  μ_1, μ_2 as
  
    dμ_1 = 1/Z e^-β
        U[(∇ξ^T∇ξ)]^-1/2 dν
        = 1/Z[-β(U + 1/2βln(∇ξ^T a ∇ξ)/ a)]
        dν_g ,
       dμ_2 = 1/Z e^-β
        Ud ν= 1/Z[-β(U + 1/2βln(∇ξ^T a ∇ξ)/( a) (∇ξ^T∇ξ))]
        dν_g ,

where again Z denotes different normalization constants.
Applying Proposition <ref> to the two probability measures expressed in
(<ref>), we can conclude that both the dynamics X_s in
(<ref>) and Y_s in (<ref>) evolve on
the submanifold Σ, and their 
invariant probability measures are given by μ_1 and μ_2, respectively.


  Under Assumptions <ref>–<ref>, we can find a neighborhood 𝒪 of Σ, such
  that P can be extended to 𝒪. Furthermore, the relations in
   (<ref>) in Appendix <ref> are still satisfied in 𝒪. Due to this fact, 
   in (<ref>)–(<ref>) we can talk about the derivatives of P at states x∈Σ. 
   


  
    
  * Notice that, similar to (<ref>), the infinitesimal generator
      of X_s in (<ref>) can be written as 

    ℒ = e^β U/β∂/∂ x_i(e^-β
      U(Pa)_ij∂/∂ x_j) .

      Using (<ref>) and (<ref>) in Appendix <ref>, we can also verify the integration by parts formula
      
    ∫_Σ (ℒf) f'  dμ_1 = 
    	∫_Σ (ℒf') f  dμ_1 = 
    	- 1/β∫_Σ (Pa ∇ f)·∇ f'   dμ_1 ,

      for any two C^2 smooth functions f, f': Σ→ℝ.

  * Using Jacobi's formula <cit.> ∂ln (∇ξ^T∇ξ)/∂ x_j =
      (∇ξ^T∇ξ)^-1_αη∂
      (∇ξ_α^T∇ξ_η)/∂ x_j and (Pa)_ij∂_jξ_α = 0, 
      the equation (<ref>) can be simplified as 
      
    dY_s^i =     - (Pa)_ij∂ U/∂ x_j ds + 1/β Q_jl∂
    (Pa)_ij/∂ x_l ds + √(2β^-1)  P_j,i 
    	dW_s^j ,

      where the matrix Q =  - ∇ξ
      (∇ξ^T∇ξ)^-1∇ξ^T. 
  In the special case when g=a =, we have ν_g = ν
      and P_j, i = P_ji=Q_ji from (<ref>). 
      Accordingly, 
  we can write the dynamics (<ref>) as
  
    dY_s^i =    
    - P_ij∂ U/∂ x_j ds 
    + 1/β P_lj∂ P_li/∂ x_j ds +
    √(2β^-1)  P_ji  dW_s^j
    
    =    
    - P_ij∂ U/∂ x_j ds 
    -1/β (Ψ^-1)_αγ  P_lj (∂^2_ljξ_α) ∂_iξ_γ ds
    + √(2β^-1)  P_ji  dW_s^j ,
    
    =    
    - P_ij∂ U/∂ x_j ds 
    + 1/β H_i ds + √(2β^-1)  P_ji  dW_s^j ,

for 1 ≤ i ≤ d, where H=H_ie_i is the mean curvature vector of
Σ (see Proposition <ref> in Appendix <ref>).
In Stratonovich form, (<ref>) can be written as 
  
    dY_s^i =    -P_ij∂ U/∂ x_j ds  
        +
      √(2β^-1) P_ji∘ dW_s^j ,   1 ≤
    i ≤ d .

      In this case, our results are accordant with those in <cit.>.




The dynamics constructed in Proposition <ref> and Theorem <ref> are
reversible on Σ, in the sense that their infinitesimal generators are
self-adjoint with respect to their invariant measures. In fact, using the same
idea, we can construct non-reversible ergodic SDEs on Σ as well. 
We will only consider the conditional probability measure μ_1, since it is
more relevant in applications and the result is also simpler.

  Let μ_1 be the conditional probability measure on Σ defined in (<ref>).
  The vector field J=(J_1, J_2, ⋯, J_d)^T=J_ie_i, defined
  on x ∈Σ, satisfies 
  
    J(x) ∈ T_x Σ,     ∀ x ∈Σ ,
    
        P_ij∂ J_j/∂ x_i + J_j∂ P_ij/∂ x_i - β J_i ∂ U/∂ x_i = 0 .

  Consider the dynamics X_s on ℝ^d which satisfies the Ito SDE 

    dX_s^i 
    =     J_i ds - (Pa)_ij∂ U/∂ x_j ds + 1/β∂
    (Pa)_ij/∂ x_j ds + √(2β^-1)  P_j,i  dW_s^j ,

for 1 ≤ i ≤ d, where β > 0 and W_s = (W_s^1, W_s^2,⋯, W_s^d)^T is a
d-dimensional Brownian motion.
  Suppose that X_0 ∈Σ, then X_s ∈Σ almost surely for s≥ 0. 
Furthermore, the unique invariant probability measure of X_s is μ_1.


The proof can be found in Appendix <ref>.

  We make two remarks regarding the non-reversible vector J. 
  
    
  * Notice that, as tangent vectors acting on functions, we have Pe_j =
  P_ij∂/∂ x_i∈ T_x Σ. Therefore, the
  condition (<ref>) indeed only depends on the value of J on Σ.
  Supposing that J and U are defined in a neighborhood 𝒪
  of Σ (see Remark <ref>), the
  condition (<ref>) can be written equivalently as 
  
    J(x) ∈ T_x Σ,     ∀ x ∈Σ ,
       ∂/∂ x_i[(P_ijJ_j)e^-β U]=0
           ,    ∀ x  Σ .

    
  * Recall that, the non-reversible dynamics on ℝ^d 
  
    d Y^i_s     = J_i
      (Y_s) ds -(a_ij∂ U/∂ x_j)
      (Y_s) ds +
      1/β∂ a_ij/∂ x_j(Y_s)  ds +
      √(2β^-1)σ_ij(Y_s)  dW^j_s ,    1 ≤ i ≤ d ,

has the invariant probability density 1/Z e^-β U, if the
vector J=(J_1, J_2, ⋯,
J_d)^T satisfies 

    (J e^-β U) = ∂(J_i e^-β U)/∂ x_i =0  ,    ∀ x ∈ℝ^d .

Comparing (<ref>) with (<ref>), it is
clear that J=J|_Σ satisfies the condition
(<ref>) of Corollary <ref> and can be used to construct non-reversible
SDEs on Σ, provided that PJ = J in the neighborhood
𝒪. Roughly speaking, in this case the vector field
J is tangential to the level sets of ξ in 𝒪. In general cases, however, 
we can not simply take J = PJ to obtain non-reversible
processes on Σ which sample μ_1, since (<ref>)
may not be satisfied. 
We refer to Remark <ref> in Section <ref> for an
alternative idea to develop “non-reversible” numerical schemes.






§ NUMERICAL SCHEME SAMPLING THE CONDITIONAL MEASURE ON Σ


Given a smooth function f : Σ→ℝ on
the level set Σ, in this section we study the numerical scheme 
(<ref>)–(<ref>) in the Introduction, which allows us to 
numerically compute the average 

    f = ∫_Σ f(x) dμ_1(x)

with respect to the conditional probability measure μ_1 in (<ref>).

To motivate the numerical scheme, let us first introduce the softly constrained dynamics, which satisfies the SDE

    d X^ϵ,i_s     = [- a_ij∂ U/∂
    x_j  -
    1/ϵa_ij∂/∂ x_j(1/2∑_α=1^k
      ξ_α^2) +
    1/β∂ a_ij/∂ x_j](X^ϵ_s)  ds +
      √(2β^-1)σ_ij(X^ϵ_s)  dW^j_s ,

where ϵ > 0, 1 ≤ i ≤ d. It is straightforward to
verify that (<ref>) has a unique invariant measure 

    dμ^ϵ(x) = 1/Z^ϵexp[-β(U(x) +
        1/2ϵ∑_α=1^k ξ_α^2(x))] dx ,  ∀ x ∈ℝ^d ,

where Z^ϵ is the normalization constant. 
As ϵ→ 0, 
the authors in <cit.>
studied the convergence of the dynamics (<ref>) itself on a
finite time horizon in the case when a=σ = id and k=1.
Closely related problems have also been studied
in <cit.>. 
Since we are mainly interested in sampling the invariant measure, 
we record the following known convergence result of the measure
μ^ϵ to μ_1. We omit its proof since it is a standard
application of the co-area formula. 

  Let f : ℝ^d →ℝ be a bounded smooth function. 
  μ_ϵ is the probability measure in (<ref>) and μ_1 is the
  conditional probability measures on Σ defined in (<ref>).
  We have 
  
    lim_ϵ→ 0∫_ℝ^d f(x)  dμ^ϵ(x)
        = 1/Z∫_Σ f(x) e^-β U(x)[(∇ξ^T
        ∇ξ)(x)]^-1/2 dν(x) =
        ∫_Σ f(x) dμ_1(x) ,

  where Z is the normalization constant given by 
  
    Z = ∫_Σ e^-β U(x)[(∇ξ^T
        ∇ξ)(x)]^-1/2 dν(x) .

  


Lemma <ref> suggests that the softly constrained
dynamics (<ref>) with a small ϵ is a good candidate to sample μ_1 on Σ. However, direct simulation of
(<ref>) is probably inefficient when ϵ is small, 
because the time step-size in numerical simulations becomes 
limited due to the strong stiffness in the dynamics. 
The numerical scheme we will study below can be viewed as a multiscale numerical method for the
dynamics (<ref>). To explain the method, let us introduce the flow map 
φ : ℝ^d × [0, +∞) →ℝ^d, defined
  by

    dφ(x,s)/ds =    - (a∇ F) (φ(x,s)) ,  φ(x,0) = x,     ∀ x∈ℝ^d ,

where the function F is 

    F(x)=1/2|ξ(x)|^2=1/2∑_α=1^kξ_α^2(x) .

Under proper conditions <cit.>, one can define the limiting map of
φ as

    Θ(x) = lim_s→ +∞φ(x,s) ,     ∀  x ∈ℝ^d .

Since ∇ F|_Σ=0 and Σ is the set consisting of all global minima of F, it is clear
that Θ: ℝ^d →Σ and Θ(x) = x, for ∀ x ∈Σ.

With the map Θ, 
we propose to approximate the average f in (<ref>) by 

    f_n = 1/n∑_l=0^n-1 f(x^(l))  ,

where n is a large number and the states x^(l) are sampled from the numerical scheme

    x^(l + 1/2)_i =    x^(l)_i + (-a_ij∂ U/∂
    x_j   + 1/β∂ a_ij/∂ x_j)  h
      + √(2 β^-1h) σ_ijη^(l)_j ,    1 ≤ i ≤
      d  ,
    
      x^(l+1) =   Θ(x^(l+1/2)) ,

starting from x^(0)∈Σ. In (<ref>),  
h>0 is the time step-size, functions a, σ, U are evaluated at
x^(l), and η^(l)=(η_1^(l), η_2^(l), ⋯,
η_d^(l))^T are independent d-dimensional standard Gaussian random
variables, for 0 ≤ l < n-1. 

  We make two comments about the scheme (<ref>)–(<ref>).
  
    
  * Since the image of Θ is on Σ, the discrete dynamics
  x^(l) stays on Σ all the time. As in the case of the softly constrained
  dynamics (<ref>), the numerical scheme has the advantage that only the 1st order
  derivatives of ξ are needed.
    
  * When a=, the numerical scheme (<ref>) becomes 

    x^(l + 1/2) =    x^(l)-∇ U(x^(l)) h + √(2 β^-1h) η^(l) , 
    
      x^(l+1) =   Θ(x^(l+1/2)) .

    
  

  At each step l ≥ 0, one needs to
  compute Θ(x^(l+1/2)). This can be done by
  solving the ODE (<ref>) starting from x^(l+1/2), using 
  numerical integration methods such as Runge-Kutta methods.
  In the following remark, we discuss issues associated with the computation of the ODE flow map Θ.
  
    
    
  * Exploiting the gradient structure of the ODE (<ref>), we can in fact
 establish exponential convergence of the dynamics φ to its limit
 Θ, at least in the neighborhood 𝒪 of Σ. For instance, we refer to 
 <cit.> and <cit.>. 
 Here, for brevity, we point out that the exponential decay of F(φ(x,s)) can
 be easily obtained and is therefore a good candidate for the convergence criterion
 in numerical implementations. Actually, under Assumption <ref>–<ref>, we can suppose 
 
    z^TΨ(x)z ≥ c_0|z|^2,   ∀ z ∈ℝ^k ,  ∀ x∈𝒪,
 for some c_0 >0, where Ψ = ∇ξ^T a ∇ξ.
 Direct calculation gives 
 
    d F(φ(x,s))/ds = -(ξ_αΨ_αη ξ_η)(φ(x,s)) ≤
       -c_0|ξ(φ(x,s))|^2 = -2c_0 F(φ(x,s)) ,

 which implies that |ξ(φ(x,s))|^2= 2F(φ(x,s)) ≤ e^-2c_0 s |ξ(x)|^2. 
 In practice, suppose that we choose the condition |ξ(φ(x,s))| ≤ϵ_tol as the stop
	criterion of ODE solvers and set Θ(x) = φ(x,s_ode) when the
	condition is met at the time s_ode.
 Then the above analysis indicates that we need to integrate the ODE (<ref>)
	until the time s_ode=max{1/c_0(ln |ξ(x^(l+1/2))| +
 ln1/ϵ_tol),0},
 which grows logarithmically as ϵ_tol→ 0.
  Since x^(l+1/2) is likely to remain close to Σ when h is small, we
  can expect that Θ(x^(l+1/2)) can be computed up to
	sufficient accuracy with affordable numerical effort. 
      
  * As a complement of the discussion above, 
	we point out that adaptivity techniques (e.g., using adaptive step-sizes) can be used to
	accelerate the computation of the flow map Θ.
	For instance, instead of (<ref>), we can consider the ODE 

    dφ̅(x,s)/ds =    - (a∇ |ξ|^2-κ) (φ̅(x,s)) ,  φ̅(x,0) = x,

where 0 ≤κ < 1. In fact, from the identity 
	
    ∇ |ξ|^2-κ = (2-κ) |ξ|^-κ∇|ξ|^2/2 = (2-κ) ∑_α=1^k
    	  ξ_α/|ξ|^κ∇ξ_α,

	we know that ODE (<ref>) is related to ODE (<ref>)
	by a rescaling of the time s. Accordingly, for each x, the
	solution φ̅(x,·) coincides
	with φ(x,·) after a reparametrization and therefore can be
	used to compute the projection Θ(x) as well. Furthermore, 
	similar to (<ref>), in this case we have 
	
    d F(φ̅(x,s))/ds = -(2-κ)[|ξ|^-κ (ξ_αΨ_αη ξ_η)](φ̅(x,s)) ≤
    	  -c_0(2-κ) [2F(φ̅(x,s))]^1-κ/2 ,

	from which we obtain |ξ(φ̅(x,s))|^κ≤ |ξ(x)|^κ
	- 2c_0(2-κ)s, and therefore φ̅(x,s) reaches the state
	Θ(x) ∈Σ before the finite time s_ode =
	|ξ(x)|^κ/2c_0(2-κ).

In applications, Θ(x) can be computed by solving the ODE (<ref>)
	with a proper κ∈ [0, 1) (and decreasing step-sizes).
	From the above discussion, in particular the identity
	(<ref>), we know that 
	this is equivalent to solving the ODE
	(<ref>) using adaptive step-sizes. 
	We refer to Examples 1–2 in Section <ref>
	for numerical validation.
	
  


Our main result of this section concerns the approximation quality of the
mean value f by the running average f_n in (<ref>), in the case when h is small and n is large. 
For this purpose, it is necessary to study the properties of the limiting
flow map Θ, since it is involved in 
the numerical scheme (<ref>).
In fact, we have the following important result, which characterizes the
derivatives of Θ by the projection map P in (<ref>).
(We refer the reader to (<ref>)–(<ref>) in Appendix <ref> for properties of P.)
    
      Let Θ be the limiting flow map in (<ref>) and P be the projection map in (<ref>). At each x ∈Σ, we have 
    
    ∂Θ_i/∂ x_j =    P_ij , 
    
          a_lr∂^2 Θ_i/∂ x_l∂ x_r 
          =   ∂ (Pa)_il/∂ x_l - 
          P_il∂ a_lr/∂ x_r ,

for 1 ≤ i, j ≤ d. 
    
    
    The proof of Proposition <ref> can be found in Appendix <ref>.

Based on the above result, we are ready to quantify the approximation error between the
estimator f_n and the mean value f.

  Suppose that both the step-size h and the number of the total steps n are fixed.
  Assume that f : Σ→ℝ is a smooth function on
  Σ and f is its mean value defined in (<ref>) with respect to the measure μ_1.
  Consider the running average f_n in (<ref>),
  which is computed by simulating the numerical scheme (<ref>) with time step-size
  h>0. Let T=nh and C denote a generic positive constant that is independent of
  h, n. We have the following approximation results. 
  
  
  * |𝐄f_n - f| ≤ C(h+ 1/T).
  
  * 𝐄|f_n - f|^2 ≤ C(h^2 + 1/T).
  
  * For any 0 < ϵ < 1/2, there is an almost surely bounded positive random variable ζ(ω), such that
      |f_n - f| ≤ Ch+ ζ(ω)/T^1/2-ϵ , almost surely.
  
  

We present the proof of Theorem <ref> in
Appendix <ref>, since it is technical and the idea follows
the standard approach developed in <cit.>, where Poisson
equation played a crucial role. However, let us emphasize that, in contrast to
<cit.>, in the current setting we are working on the
submanifold Σ and, furthermore, the map Θ is involved in our
numerical scheme. In particular, in the proof we use the Poisson equation
related to the generator ℒ in (<ref>) of the process (<ref>),
based on the fact that μ_1 is the invariant measure of
(<ref>) (This is the place where
Theorem <ref> in Section <ref> is used in order
to establish Theorem <ref>).

  Theorem <ref> concerns the long time
  behavior of the scheme (<ref>) with a small
  step-size, i.e., large T and small h. This is often relevant in molecular dynamics simulations.
While the estimates of Theorem <ref> are
  stated in terms of the variables h
  and T, we should point out that the time T=nh and therefore it depends on both h and n.
  Alternatively (and more precisely), the estimates can be expressed using the independent
  variables h and n. For instance, for the mean square error estimate, we have 
  
    𝐄|f_n - f|^2 ≤ C(h^2 + 1/nh) .

  Therefore, for a fixed (large) total sample number n,
  we can conclude that the optimal upper bound in (<ref>)
  is 𝒪(n^-2/3) and is achieved when h=𝒪(n^-1/3).
  We refer to <cit.> for related discussions.
  

In applications, the conditional probability measure μ_1 often satisfies the
following Poincaré inequality <cit.>
       
    _μ_1(f) := ∫_Σ (f-f)^2 dμ_1 
    	 ≤ -1/K∫_Σ (ℒf) f dμ_1 =
    	 1/Kβ∫_Σ (Pa ∇ f)·∇ f dμ_1  ,

       for all f : Σ→ℝ such that the right hand
       side of the above inequality is finite, where K>0 is the  
Poincaré constant, ℒ is the infinitesimal generator 
(<ref>), and the identity (<ref>) in
Remark <ref> has been used. Under this condition, the mean square error
estimate in Theorem <ref> can be 
improved (i.e., the constant in front of the 𝒪(T^-1) term 
is small when K is large) and we have the following corollary (The proof is
in Appendix <ref>).

  Under the same assumptions in
  Theorem <ref> and further assuming that 
  μ_1 satisfies the Poincaré inequality
  (<ref>), we have 
  
    𝐄|f_n - f|^2 ≤2C_1 _μ_1(f)/KT + C_2(h^2 +
        h/T+1/T^2) ,

  where C_1 is any constant larger than 1, the constant C_2 
  depends on the choice of C_1 but is independent of both h and n.
  


  The idea of using the map Θ in the constraint step of the numerical
  scheme (<ref>) is motivated by the softly
  constrained (reversible) dynamics (<ref>). It is natural to
  consider whether certain “non-reversible” numerical scheme can
  be obtained using the same idea. In fact, let A∈ℝ^d×
  d be a constant skew-symmetric matrix such that A^T=-A. The softly constrained
  (non-reversible) dynamics 
  
    d X^ϵ,A, i_s     = [A_ij∂/∂
    x_j(U+1/2ϵ∑_α=1^kξ_α^2) - a_ij∂ U/∂
    x_j  -
    1/ϵa_ij∂/∂ x_j(1/2∑_α=1^k
      ξ_α^2) +
    1/β∂ a_ij/∂ x_j](X^ϵ,A_s)  ds
    
        + √(2β^-1)σ_ij(X^ϵ,A_s)  dW^j_s ,

indeed has the same invariant measure μ^ϵ in (<ref>).
Based on this fact, a reasonable guess of the “non-reversible” numerical
scheme that samples the conditional measure μ_1 is 
the multiscale method of (<ref>), i.e.,

    x^(l + 1/2)_i =    x^(l)_i + (A_ij∂ U/∂ x_j  
        -a_ij∂ U/∂ x_j   + 1/β∂
        a_ij/∂ x_j)(x^(l))  h
      + √(2 β^-1h) σ_ij(x^(l)) η^(l)_j ,    1 ≤ i ≤
      d  ,
    
      x^(l+1) =   Θ^A(x^(l+1/2)) ,

where Θ^A(x)=lim_s→ +∞φ^A(x,s) is the limit of the (non-gradient) flow map 

    dφ^A(x,s)/ds =    -((a-A)∇ F)(φ^A(x,s)) ,  φ^A(x,0) = x,     ∀ x∈ℝ^d ,

with the same function F in (<ref>).
We expect that the long time sampling error estimates of the numerical scheme
(<ref>) can be studied following the same approach of this
section as well.
For this purpose, however, it is necessary to handle the non-gradient term in the ODE (<ref>),
which brings difficulties when calculating the derivatives of the map
Θ^A (cf. Proposition <ref> as well as
its proof in Appendix <ref>). We will postpone the analysis in the future work
and readers are referred to Example 1 in Section <ref> for numerical
validation of the scheme (<ref>–(<ref>).
   


In the literature, reversible Metropolis random walk on submanifold <cit.> and Hybrid Monte Carlo algorithm <cit.> 
have been proposed to sample distributions on submanifolds. A comprehensive
comparison between our scheme (without Metropolis step) and these Metropolis-adjusted
approaches is a complicate task and goes beyond the scope of the current paper. 
We only discuss this issue briefly in the following remark. 

  We compare the following three aspects.  
  
    
  * Constraint step.  In the scheme
      (<ref>), the map Θ in (<ref>) is
      used to project the state x^(l+1/2) to the submanifold Σ.
      Under mild assumptions, the gradient structure of the ODE (<ref>)
allows to define Θ for all states at which the map ξ is C^2
      smooth (Assumption <ref>).
      Differently, in <cit.>, newly proposed states
      are projected back to the submanifold by solving a nonlinear system. 
      One usually uses Newton's method to find the solution of the
      system, with the hope that the convergence can be achieved within a few
      iteration steps (success), thanks to the quadratic convergence rate of Newton's method.
      In practice, however, it may happen that either the solution does not
      exist or, even if the solution exists, the Newton's method does not
      converge (due to local convergence). In these two cases, the constraint step
      ends without finding a new state on the submanifold (no success). Although the Markov chain samples the correct invariant distribution regardless whether the constraint step is
      successful or not <cit.>, the sampling efficiency is
      affected by the success rate of the constraint step. 
    
  * Computational complexity.
      Suppose the computational complexity of evaluating the d× k matrix
      ∇ξ is 𝒪(k· d) and n states are sampled in total.
      In each step of both the scheme (<ref>) and the Metropolis-adjusted
      methods in <cit.>, the
      major computational effort is devoted to the constraint step, i.e.,
      either computing the map Θ by integrating ODE or solving equations using Newton's method. 
      For the scheme (<ref>) with a=, the overall
      computational complexity of the constraint step is therefore 𝒪(n
      · k · d· s_ode/Δ s), where s_ode  and Δ s
      are the average final time (see
      Remark <ref>) and the average
      step-size in the ODE integration, respectively. 
      For the Metropolis-adjusted methods in <cit.>,
in each Newton iteration it is necessary to compute the matrix product ∇ξ^T(x)∇ξ(x') for two states x,x'. 
      Therefore, the overall computational complexity is 𝒪(n·
      k^2· d· N_iter), where N_iter is the average total
      iteration steps of Newton's method.
      In practice, one can implement the method in <cit.> in a way such that Newton's method 
      ends within a few Newton iterations, e.g., N_iter≤ 10.  On the
      other hand, the ODE integration in the scheme (<ref>) requires more
      iteration steps (for the examples in Section <ref>, 20-40 steps are needed), i.e., s_ode/Δ s ≥ N_iter.
      However, when comparing the computational cost of both methods, we
      should keep in mind that an ODE iteration step is generally
      cheaper than a Newton iteration step, since the latter involves 
      both matrix-matrix multiplication and solving linear systems. (The
      cost of solving k× k linear systems is not major and therefore for
      simplicity is not included in the estimation above.)
      While the computational cost of Newton's method is smaller when k is
      small, the ODE integration becomes faster for medium or large k.
      We refer to Example 2 in Section <ref> for
      numerical comparison on the computational cost of both methods.
    
  * Choice of step-size.
      To apply the scheme (<ref>), one usually
      chooses a suitably small step-size h and runs the scheme for sufficient many steps (large n).
      In concrete applications, one often needs to tune the step-size h,
      keeping in mind that a large h will lead to large bias, while a unnecessarily small
      h will result in large correlations. See Remark <ref>. 
      For the Metropolis-adjusted
      methods <cit.>, the
      choice of the step-size (in the proposal step) is in fact a more
      delicate issue. Although the Markov chain remains unbiased for large
      step-sizes, the sampling efficiency will be possibly limited due to a low
      acceptance rate in the Metropolis step. This issue has been discussed
      in <cit.>, where the performance 
      with different step-sizes has been numerically investigated.
      Besides the acceptance rate in the Metropolis step, the success rate of the 
      constraint step also depends on the step-size used in the proposal step.
      Taking the (d-1)-dimensional unit sphere 𝕊^d-1 as an
      example, it is not difficult to see that there won't be corresponding
      projected state on the sphere (i.e., the solution of the constraint
      equation does not exist), if
      the norm of the tangent vector generated in the proposal step
      (see <cit.>) is larger than one.
      This implies that the success rate of the constraint step will decrease
      when we increase the step-size in the proposal step.
      To summarize, it is important to choose the step-size in the
      Metropolis-adjusted method in <cit.> properly, such that both the acceptance rate in the
      Metropolis step and the success rate of the constraint step are not too small.
  
  


Before concluding, let us point out that the approach used in the above proof
allows us to study other numerical schemes on Σ as well. As an example,
we consider the projection from ℝ^d to Σ along geodesic
curves (instead of using the flow map (<ref>)–(<ref>)) defined by the metric g=a^-1 in (<ref>), i.e., the metric on ℳ=(ℝ^d, g).
Let d be the distance function on ℝ^d induced by the
metric g in (<ref>). We introduce the projection function

    Π(x) = {y | d(x,y) = d(x,Σ),  y ∈Σ} ,  ∀ x ∈ℝ^d .

Clearly, we have Π|_Σ=|_Σ. Given any x ∈Σ, there is a
neighborhood Ω⊂ℝ^d of x such that Π|_Ω is a single-valued map.
Furthermore, applying 
inverse function theorem, we can verify that Π is smooth on Ω.
Similar to Proposition <ref>, we need the following result which connects the derivatives of Π to the
projection map P in (<ref>). (Note that, comparing to the
derivatives of the map Θ in (<ref>), there is an extra term in the second equation of
(<ref>).)  Its proof is given in Appendix <ref>.

  Let Π = (Π_1, Π_2, ⋯, Π_d)^T: ℝ^d →Σ be the projection function in
  (<ref>), where Π_i : ℝ^d →ℝ are smooth functions, 1≤ i ≤ d. 
 For x ∈Σ∩Ω, we have 
 
    ∂Π_i/∂ x_j = P_ij , 
       a_lr∂^2 Π_i/∂ x_l∂ x_r  
         = 
          - P_il∂ a_lr/∂ x_r
          + ∂ (Pa)_il/∂ x_l
          +
          1/2 (Pa)_il∂lnΨ/∂
    	x_l ,

for 1 ≤ i,j ≤ d, where Ψ = ∇ξ^T a ∇ξ.
 

Now we are ready to study the numerical scheme 

    x^(l + 1/2)_i =    x^(l)_i + (-a_ij∂ U/∂
      x_j   + 1/β∂ a_ij/∂ x_j)(x^(l))  h
      + √(2 β^-1h) σ_ij(x^(l)) η^(l)_j ,    1 ≤ i ≤
      d  ,
    
      x^(l+1) =   Π(x^(l+1/2)) ,

where x^(0)∈Σ, and the map Π (instead of Θ) is used in each step to project the states x^(l+1/2) back to Σ.

  Assume that f : Σ→ℝ is a smooth function on
  Σ and f is its mean value 
  
    f = ∫_Σ f(x) dμ(x) ,

  with respect to the probability measure 
  
    dμ=1/Z e^-β U√( (∇ξ^T a
        ∇ξ)/ (∇ξ^T ∇ξ))  dν .

  Consider the running average f_n in (<ref>),
  which is computed by simulating the numerical scheme (<ref>) with time step-size
  h>0. Let T=nh and C denote a generic positive constant that is independent of
  h, T. We have the following approximation results. 
  
  
  * |𝐄f_n - f| ≤ C(h+ 1/T).
  
  * 𝐄|f_n - f|^2 ≤ C(h^2 + 1/T).
  
  * For any 0 < ϵ < 1/2, there is an almost surely bounded positive random variable ζ(ω), such that
    |f_n - f| ≤ Ch+ ζ(ω)/T^1/2-ϵ , almost surely.
  
  

We omit the proof since it resembles the proof of Theorem <ref>.

  For the projection map Π induced by a general metric g=a^-1
  or, equivalently, by a general (positive definite) matrix a,
  implementing the numerical scheme (<ref>) is not as easy as
  the numerical scheme (<ref>). We decide to omit
  the algorithmic discussions, due to the fact that the probability measure
  (<ref>) seems less relevant in applications.
  However, it is meaningful to point out that, when a=, the above result is relevant to the one in <cit.>. In this case, the probability
  measure in (<ref>) reduces to μ_2=1/Z e^-β Udν in (<ref>) and the numerical scheme
  (<ref>) can be formulated equivalently using
  Lagrange multiplier. We refer to <cit.> for comprehensive numerical details. 



§ NUMERICAL EXAMPLES


In this section, we study three concrete examples. In the first example, we investigate the different schemes in Section <ref>. In particular, the sampling performance of the constrained schemes using different maps Θ, Θ^A, and Π, as well as the performance of the unconstrained Euler-Maruyama discretization of the SDE (<ref>), will be compared.
In the second example, we compare the computational costs of the scheme (<ref>)
and the Metropolis-adjusted algorithm introduced in <cit.>.
In the last example, we show that in some cases it is helpful to consider non-constant matrices σ and a.
The C/C++ code used for producing the numerical results in the following examples is available at: <https://github.com/zwpku/sampling-on-levelset>.


  §.§.§ Example 1: Comparison of schemes using different projection maps

Let us define ξ: ℝ^2→ℝ by 

    ξ(x) = 1/2(x_1^2/c^2+x_2^2 -1) ,   ∀ x = (x_1, x_2)^T ∈ℝ^2 ,

with the constant c=3. The level set Σ={(x_1,
x_2)^T | x_1^2/c^2+x_2^2=1} is an ellipse in ℝ^2.
We have ∇ξ = (x_1/c^2,x_2)^T and therefore
(∇ξ^T∇ξ) = |∇ξ|^2 = x_1^2/c^4 + x_2^2.
For simplicity, we choose the potential U=0 and the matrices
a=σ=∈ℝ^2× 2.
The two probability measures in (<ref>) on Σ are 

    dμ_1 = 1/Z (x_1^2/c^4 + x_2^2)^-1/2 dν ,   
      dμ_2 = 1/Z  dν ,

where Z denotes two different normalization constants and ν is the
surface measure on Σ. Since Σ is a one-dimensional manifold, it
is helpful to consider the parametrization of Σ by

    x_1 = ccosθ ,   x_2 = sinθ ,

where the angle θ∈ [0, 2π].
 Applying the chain rule ∂/∂θ
= -csinθ∂/∂ x_1 + cosθ∂/∂
x_2, we can obtain the expressions of μ_1, μ_2  under this coordinate as 

    dμ_1 = 1/Z dθ ,   dμ_2 = 1/Z(c^2sin^2θ +
      cos^2θ)^1/2 dθ .

With these preparations, we proceed to study the following four numerical approaches.

  
  * Numerical scheme (<ref>) using Θ. Since U≡ 0 and a=, 
(<ref>) becomes  

    x^(l + 1/2) =    x^(l) + √(2 β^-1h) η^(l) , 
    
      x^(l+1) =   Θ(x^(l+1/2)) ,

where Θ(x) is the limit of the flow map φ, given by 

    ẏ_1(s) = -ξ(y(s)) y_1(s)/c^2 ,  ẏ_2(s) = -ξ(y(s)) y_2(s) ,      s ≥ 0 ,

starting from y(0) = x.
  
  * Numerical scheme (<ref>)–(<ref>) using Θ^A.
    Let us choose the skew-symmetric matrix  
    
    A=
        [    0  1/2; -1/2    0 ].

    Since U≡ 0 and a=, we have 

    x^(l + 1/2) =    x^(l) + √(2 β^-1h) η^(l) , 
    
      x^(l+1) =   Θ^A(x^(l+1/2)) ,

where Θ^A(x) is the limit of the flow map φ^A, given by

    ẏ_1(s) = -ξ(y(s)) (y_1(s)/c^2 - y_2(s)/2)  ,  ẏ_2(s) = -ξ(y(s)) (y_1(s)/2c^2 + y_2(s)) ,      s ≥ 0 ,

starting from y(0) = x.
  
  * Numerical scheme (<ref>) using Π.
Similarly, since U≡ 0 and a=, (<ref>) becomes  

    x^(l + 1/2) =    x^(l) + √(2 β^-1h) η^(l) , 
    
      x^(l+1) =   Π(x^(l+1/2)) ,

where Π is the projection map onto Σ, defined in (<ref>). 
  
  * Euler-Maruyama discretization of the SDE
    (<ref>). 
    Notice that we have Pa = P, and it is straightforward to compute 
    
    P_11 = c^4x_2^2/x_1^2+c^4x_2^2 ,   
          P_12 = P_21 = -c^2x_1x_2/x_1^2+c^4x_2^2 ,   
          P_22 = x_1^2/x_1^2+c^4x_2^2 .

  Therefore, discretizing (<ref>), we obtain 

    x^(l + 1)_1 =    x^(l)_1 + 1/βc^4(c^2-2) x_1x_2^2-c^2x_1^3/(x_1^2+c^4x_2^2)^2 h 
        + √(2 β^-1h) (c^4x_2^2/x_1^2+c^4x_2^2 η^(l)_1
        -c^2x_1x_2/x_1^2+c^4x_2^2 η^(l)_2) 
    
        x^(l + 1)_2 =    x^(l)_2 + 1/β(1-2c^2)x_1^2x_2-c^6x_2^3/(x_1^2+c^4x_2^2)^2 h 
        + √(2 β^-1h) (
        -c^2x_1x_2/x_1^2+c^4x_2^2 η^(l)_1
        +x_1^2/x_1^2+c^4x_2^2 η^(l)_2 ) .



Based on Theorems  <ref>, <ref>,  and
Remark <ref>, we
study the performance of the schemes (<ref>),  
(<ref>), and (<ref>) in sampling the conditional
measure μ_1, as well as the performance of the scheme (<ref>)
in sampling the measure μ_2. 

In the numerical experiment, we choose β=1.0 in each of the above schemes.
For the first scheme using Θ, we simulate (<ref>) 
for n=2× 10^7 steps with the step-size h=0.01. In each constraint step, 
Θ(x^(l+1/2)) is computed by solving the ODE
(<ref>) starting from y(0) = x^(l+1/2) until
the time when |ξ(y(s))| < 10^-8 is satisfied, using the 3rd order (Bogacki-Shampine) Runge-Kutta
(RK) method.  
The adaptivity technique in the second point of Remark <ref> is used with κ=0.5. 
The step-size for solving the ODE is set to Δ s = 0.1 initially and
is divided by 2.0 whenever we find that the value of |ξ| is not decreasing.
(The numerical error of Θ is 3.7 × 10^-7 on average, comparing to the
reference solution that is obtained by solving the ODE with κ=0 and
the fixed step-size Δ s= 0.001.)
On average, we observe that 25 iterations of the RK method are needed in
each constraint step in order to meet the criterion |ξ(y(s))| < 10^-8.

For the second scheme using Θ^A, we simulate (<ref>) 
for n=2× 10^7 steps with the step-size h=0.005. 
Notice that, a slightly smaller step-size h is used, because in this case the non-gradient ODE flow 
(<ref>) produces a drift force on the level set Σ.
In each constraint step, Θ^A(x^(l+1/2)) is computed by solving the ODE
(<ref>) in the same way (with the same parameters) as we did in the first scheme.
On average, we find that 23 iterations of the RK method are needed in
each constraint step in order to meet the criterion |ξ(y(s))| < 10^-8.

For the third scheme using Π, (<ref>) is simulated for n=2× 10^7 steps with the step-size h=0.01. Using the parametrization (<ref>), we have Π(x) = (ccosθ^*, sinθ^*)^T, where 

    θ^* = _θ∈ [0, 2π]((x_1-ccosθ)^2 +
      (x_2-sinθ)^2) ,    x=(x_1,x_2)^T .

Therefore, in each step, Π(x^(l+1/2)) is computed by solving (<ref>) 
using the simple gradient descent method. The step-size is fixed to Δ t
= 0.1 and the gradient descent iteration terminates when the derivative of the objective
function in (<ref>) has an absolute value that is less than 10^-8.
On average, it requires 32 gradient descent iterations in each step in order to meet the convergence criterion.

Let us make a comparison among the three schemes (<ref>), (<ref>) and (<ref>).
From Figure <ref> and Figure <ref>, 
we can see that the three maps Θ, Θ^A and Π indeed have different effects.
Roughly speaking, comparing to the projection map Π, both Θ
and Θ^A tend to map states towards one of the two vertices (± c,
0), where |∇ξ| are smaller, while Π^A introduces a further rotational force on Σ.
Based on the states generated from these three schemes,
in Figure <ref> we show the empirical probability densities of the parameter θ in (<ref>).  
From the agreement between the empirical densities and the densities
computed from the analytical expressions in (<ref>), we can
make the conclusion that 
the trajectories generated from the two schemes using Θ and Θ^A indeed sample the
probability measure μ_1, while the trajectory generated from the scheme
using Π samples μ_2. 

Lastly, concerning the fourth scheme, we simulate (<ref>) for n=10^7
steps using the step-size h = 0.0001.
In this case, we find that it is necessary to choose a small step-size h in order to keep the trajectory close to the level set Σ.
As can be seen from Figure <ref>, 
even with this smaller step-size h=0.0001,
the generated trajectory departs from the level set Σ. This
indicates the limited usefulness of the direct Euler-Maruyama
discretization of the SDE (<ref>) in long time simulations.






  §.§.§ Example 2: Numerical comparison with the
Metropolis-adjusted method on the special orthogonal group SO(11).

In this example, we compare the computational efficiency between our scheme 
(<ref>) using the flow map Θ and the Metropolis-adjusted method
introduced in <cit.>. We consider the special orthogonal
group SO(11), which consists of  orthogonal matrices of size 11 × 11 with determinant equals to 1. 
This example is taken from <cit.>. The authors there applied their method to the estimation of the mean value of the function f(x) = (x), 
i.e, the trace of the matrix x, where x follows the surface measure of SO(11).
The manifold SO(11) can be viewed as (one connected component of) the level set 
of the map ξ: ℝ^121→ℝ^66, which includes
all the row ortho-normality constraints. Readers are referred to the original work <cit.>
for a detailed introduction on the example.

In this numerical study, we implement both the scheme
(<ref>) and the (Metropolis-adjusted) algorithm
in <cit.> to estimate the mean value of (x).
Notice that, since (∇ξ^T∇ξ) is constant, the conditional
measure μ_1 in (<ref>) coincides with the surface measure of SO(11) when we choose the
potential U≡ 0.
In both cases, we generate n=10^6 samples on the same laptop (CPU: Intel Core i5,
2.60GHz, 4 cores; system: Ubuntu 18.04).
For the scheme (<ref>), we choose the step-size
h=0.022. The map Θ is computed
by integrating the ODE (<ref>) with a=, until
the condition |ξ(φ(x,s))|<10^-9 is satisfied.
To accelerate the ODE integration, we have applied the adaptivity technique in the second point of Remark <ref> with κ=0.5.
Starting from the initial step-size Δ s=0.2, the step-size used in the ODE
integration is divided by 2.0 whenever we find that the value of |ξ| is not decreasing.
(The numerical error of Θ is 2.0 × 10^-4 on average, comparing to the
reference solution that is obtained by solving the ODE with κ=0 and
the fixed step-size Δ s= 0.002.)
Furthermore, the new state will be discarded (and resampled) if its determinant equals to -1.
With these parameters, we observe that on average 37 Runge-Kutta iterations
are needed for each evaluation of the map Θ.
In total, it takes 2676.9 seconds to generate n=10^6
samples, while the estimated mean value of (x) is 3.8 × 10^-3
with a statistical error 3.9 × 10^-3.
For the algorithm in <cit.>, 
the maximal number of Newton steps is set to 10 and the proposal length
scale is chosen to be 0.257 [The roles of the proposal length scale
in <cit.> and the step-size h in the
scheme (<ref>) are different. 
The proposal length scale 0.257 used in this example corresponds to a step-size 0.033
(≈ 0.257^2/2) in the scheme (<ref>).].
In our experiment, we find that this proposal length scale (different from the
one used in <cit.>) leads to slightly smaller correlation time.
Within the entire computation, the success rate of the Newton's method (i.e.,
the rate that the Newton's method converges) is 67.2% and 
each time it takes 5-6 iterations 
on average for the Newton's method to reach convergence (the convergence
criteria is |ξ(x)| < 10^-9).
In total, it takes 7315.9 seconds to generate n=10^6 samples. The
estimated mean value is -3.8× 10^-3 and the statistical error is 3.8× 10^-3. 

The empirical density distributions of (x) using both the scheme
(<ref>) and the algorithm in <cit.> are shown in the left
plot of Figure <ref>, while the autocorrelation functions are
plotted in the right plot of Figure <ref>.
From these results, we can conclude that in this example both approaches provide similar
statistical estimations (The autocorrelation time using Metropolis-adjusted
method is slightly smaller with the above parameters.). At the same time, the total computational time using the scheme (<ref>)
is less than half of the computational time required by the algorithm in <cit.>.
For this example, although the average number of Newton steps in the latter algorithm is
smaller than the average number of ODE iterations, the computational cost of
each Newton step is indeed larger. We refer to Remark <ref> for the
comparison of computational complexity of both approaches.



  §.§.§ Example 3: Removing stiffness by choosing a non-constant matrix a

In this example, we choose the reaction coordinate function 

    ξ(x) = ξ(x_1, x_2,⋯, x_d) =
    1/2(x_1^2 + x_2^2 + ⋯ + x^2_d - 1) .

Correspondingly, the level set 
 
    Σ={(x_1, x_2,⋯, x_d) ∈ℝ^d |  x_1^2 + x_2^2 +⋯ + x_d^2 = 1}

is the (d-1)-dimensional unit sphere, and we have 

    ∇ξ = (x_1,x_2,⋯,x_d)^T,   ∇ξ^T∇ξ = ∑_i=1^d x_i^2 .

In the following, we give an example to show that in some applications it is
helpful to use a non-constant matrix a in the numerical scheme (<ref>). 
Briefly speaking, varying the matrix a properly allows to rescale the scheme along different directions. 
It has a preconditioning effect when different time scales (stiffness) exist. 

Consider d=3 and the potential U = θ^2/2ϵ, where
ϵ>0 is a small parameter, θ is the angle of the state x = (x_1, x_2, x_3) under the  
spherical coordinate system 

    x_1 = ρcosθcosφ ,    x_2 = ρcosθsinφ ,    x_3 = ρsinθ ,

where ρ≥ 0, θ∈ [-π/2,π/2], and φ∈
[0, 2π]. We can verify that 

    ∇θ = 1/ρ^2(-x_1x_3/(x_1^2+x_2^2)^1/2, 
      -x_2x_3/(x_1^2+x_2^2)^1/2, 
      (x_1^2+x_2^2)^1/2)^T .

Correspondingly, with the choice of σ=a=, the scheme (<ref>) is 

    x^(l + 1/2) =    x^(l)-1/ϵ (θ∇θ)(x^(l)) h + √(2 β^-1h) η^(l) , 
    
      x^(l+1) =   Θ(x^(l+1/2)) ,

where ∇θ is given in (<ref>). 
Notice that, the coefficients in (<ref>) are 𝒪(1/ϵ) when
ϵ is small.
In particular, it implies that sampling  the invariant measure using (<ref>)
will be inefficient when ϵ is small, since the step-size h
will be severely limited due to the large magnitude of the coefficients in (<ref>). 

On the other hand, based on the form of U and the expression
(<ref>), we consider the orthogonal vectors

    σ_1 =    (x_1 , x_2 , x_3)^T=∇ξ ,  σ_2 = (x_2 , -x_1 , 0)^T ,
    σ_3 =   (-√(ϵ) x_1x_3 /(x_1^2+x_2^2)^1/2, 
        -√(ϵ) x_2x_3/(x_1^2+x_2^2)^1/2, 
        √(ϵ) (x_1^2+x_2^2)^1/2)^T=√(ϵ)ρ^2
        ∇θ ,

and we define σ =
(σ_1,σ_2, σ_3) ∈ℝ^3 × 3.
Direct calculation shows that 

    a=σσ^T = 
      [ x_1^2 + x_2^2+ϵ x_1^2x_3^2/x_1^2+x_2^2             ϵ x_1x_2 x_3^2/x_1^2+x_2^2                           (1-ϵ) x_1x_3;             ϵ x_1x_2 x_3^2/x_1^2+x_2^2 x_1^2+x_2^2+ϵ
      x_2^2x_3^2/x_1^2+x_2^2                           (1-ϵ) x_2x_3;                           (1-ϵ) x_1x_3                           (1-ϵ) x_2x_3               x_3^2 + ϵ(x_1^2 + x_2^2) ] .

Correspondingly, using (<ref>), the scheme (<ref>) becomes 

    x^(l + 1/2)_i =    x^(l)_i+[-θ∂θ/∂ x_i+1/β∂ a_ij/∂ x_j](x^(l)) h +
        √(2 β^-1h) σ_ij(x^(l)) η^(l)_j ,    1
        ≤ i ≤ 3 ,
    
      x^(l+1) =   Θ(x^(l+1/2)) ,

where Θ(x) is the limit of the ODE flow 

    ẏ(s) = - ξ(y(s))(2ξ(y(s))+1) y(s) ,    y(0) = x .

Importantly, in contrast to (<ref>), the scheme
(<ref>)–(<ref>) is no longer stiff when ϵ is small.

Now we compare the numerical efficiency between the
schemes (<ref>) and (<ref>)–(<ref>).
First of all, since the surface measure on Σ satisfies dν =
cosθ dθ dφ, we know that the target measure is 

    dμ = 1/Z e^-βθ^2/2ϵ dν = 
      1/Z e^-βθ^2/2ϵcosθ  dθ dφ .

In the numerical study, we choose ϵ = 0.005 and generate n=10^7 states for both schemes.
For the scheme (<ref>) which corresponds to a=, 
we use both a small step-size h=0.0002 and a (relatively) larger step-size h=0.005, while
we choose a large step-size h=0.01 in the scheme (<ref>)–(<ref>) . The empirical probability
densities of the angles θ, φ for the two schemes are shown
in Figure <ref> and Figure <ref>, respectively.
From Figure <ref>, we see that the step-size h has to be small (h=0.0002) in (<ref>) in order to produce the
correct probability density of the angle θ (left plot). However, with such a small
h, the estimated empirical density of the angle φ (right plot) is
still noisy with n=10^7. On the other hand, for the scheme 
(<ref>)–(<ref>) which corresponds
to the matrix a in (<ref>),
Figure <ref> shows that the probability densities of both angles
θ, φ are well approximated using the large step-size h=0.01.
Therefore, we conclude that in this example choosing the non-constant matrix a in (<ref>) indeed 
helps improve the sampling efficiency.






§ CONCLUSIONS


Ergodic diffusion processes on a submanifold of ℝ^d and related numerical sampling schemes have been considered
in this work. A family of SDEs has been obtained whose invariant measures coincide with
the given probability measure on the submanifold. In particular, for the
conditional probability measure, we found that the corresponding SDEs
have a relatively simple form. We proposed and analyzed a consistent numerical scheme
which only requires 1st order derivatives of the reaction coordinate function.
Different sampling schemes on the submanifolds are numerically evaluated. 

The current work extends results in the literature and may further contribute to 
both the analysis and the development of numerical methods on
related problems, in particular problems in molecular dynamics such as free energy calculation and model
reduction of high-dimensional stochastic processes. 
Closely related to the current paper, the following topics could be considered. 
First, the “non-reversible” scheme (<ref>) is supported by a simple
numerical example but theoretical justification still needs to be investigated. This will be
considered in future following the approach described in Remark <ref>. 
Second, the constrained numerical schemes in the current work do
not involve system's momentum variables. In view of the
work <cit.>, it is interesting to study the
Langevin dynamics under different constraints (such as certain variants of the map Θ used in this work).
Third, there is a research interest in the literature to study the effective dynamics of
molecular systems along a given reaction coordinate ξ. The coefficients of the
effective dynamics are usually defined as averages on the level set of
ξ <cit.>. As an application of the numerical scheme proposed in this work, we will study
numerical algorithms to simulate the effective dynamics. This topic is related to the
heterogeneous multiscale methods <cit.> and the equation-free
approach <cit.> in the literature. 


§ ACKNOWLEDGEMENT

This work is funded by the Einstein Center of Mathematics (ECMath) through project CH21.
The author would like to thank Gabriel Stoltz for stimulating discussions on
constrained Langevin processes at the Institut Henri Poincaré - Centre
Émile Borel during the trimester “Stochastic Dynamics Out of Equilibrium”. 
The author appreciates the hospitality of this institution.
The author also thanks the anonymous referees for their valuable comments and
criticism which helped improve the manuscript substantially.


§ USEFUL FACTS ABOUT THE RIEMANNIAN MANIFOLD ℳ


In this section, we present technical details of Section <ref> related to the Riemannian 
manifold ℳ=(ℝ^d, g), where g=a^-1. 
The main result is Proposition <ref>, where we give the expression of
the Laplacian-Beltrami operator Δ^Σ on the level set Σ in (<ref>), viewed as
a submanifold of ℳ. Before that, we first introduce some
notations and quantities related to ℳ and Σ.
Readers are referred to <cit.>
for related discussions on general Riemannian manifolds. 

Under Assumption <ref>, given two vectors u = (u_1, u_2, ⋯, u_d)^T, v = (v_1, v_2, ⋯, v_d)^T,
we consider the space ℝ^d with the weighted inner product

    g(u, v) = ⟨u, v⟩_g = u_i (a^-1)_ij v_j .

The inner product in (<ref>) defines a Riemannian metric g on ℝ^d and 
we denote by ℳ = (ℝ^d, g)
the Riemannian manifold ℝ^d endowed with this metric.

Notice that ℳ as a manifold is quite special
(simple), in that 
it has a natural global coordinate chart which is given by the usual Euclidean coordinate.
Since we will always work with this coordinate, we will
not distinguish between tangent vectors (operators acting on functions) and their coordinate representations
(d-dimensional vectors). 
In particular,  e_i denotes the vector whose ith component equals
to 1 while all the other d-1 components equal to 0, where 1 ≤ i ≤ d.
At each point x ∈ℳ, vectors e_1, e_2, ⋯, e_d form a basis of the
tangent space T_xℳ and under this basis we have g=a^-1, as can be seen from (<ref>).

Denote by ^ℳ, ^ℳ the
gradient and the divergence operator on ℳ, respectively.
For any smooth function f: ℳ→ℝ, it is direct to verify that 

    ^ℳ f = g^ij∂ f/∂ x_je_i = (a∇ f)_i e_i  ,
 
where g^ij=(g^-1)_ij = a_ij, and ∇ f 
denotes the ordinary gradient operator for functions on the Euclidean space ℝ^d.
For simplicity, we will also write ∂_i f for the partial derivative with respect to x_i, and (a∇
f)_i to denote the ith component of the vector a∇ f, i.e., 
∂_i f = ∂ f/∂ x_i, and 
(a∇ f)_i = a_ij∂ f/∂ x_j = a_ij∂_jf.

The Laplace-Beltrami operator on ℳ is defined by
Δ^ℳ f = ^ℳ (^ℳ
f). Equivalently, we have Δ^ℳ f =
(^ℳf), where  ^ℳ is
the Hessian operator on ℳ. 
The integration by parts formula on ℳ has the form

    ∫_ℳ (Δ^ℳ f_1) f_2 dm = - ∫_ℳ⟨^ℳ f_1,
      ^ℳ f_2⟩_g  dm = ∫_ℳ
      (Δ^ℳ f_2)f_1 dm ,

for ∀ f_1,  f_2 ∈ C^∞_0(ℳ),
where dm= ( g)^1/2 dx = ( a)^-1/2 dx is the volume form, and C^∞_0(ℳ) consists of all smooth functions on ℳ with compact support.

Besides the vector basis e_1, e_2, ⋯, e_d, the vectors 

    σ_i =
    (σ_1i, σ_2i,⋯, σ_di)^T,    1 ≤ i ≤ d ,

will also be useful. Note that a=σσ^T=g^-1 implies ⟨σ_i,
σ_j⟩_g = (a^-1)_rlσ_riσ_lj = δ_ij.
In other words, σ_1, σ_2, ⋯, σ_d form
an orthonormal basis of T_xℳ at each x ∈ℳ.

Denote by ∇^ℳ the Levi-Civita connection on ℳ.
Given x ∈ℳ and a tangent vector v∈ T_xℳ, 
∇^ℳ_v
is the covariant derivative operator on ℳ along the vector v.
For two vectors u = (u_1, u_2, ⋯, u_d)^T, v = (v_1, v_2, ⋯, v_d)^T,
the Hessian of a smooth function f: ℳ→ℝ is defined as 

    ^ℳ f(u,v) = u(vf) -
      (∇^ℳ_uv) f =
      u_iv_j ^ℳ f(
      e_i,e_j) = u_iv_j (∂^2 f/∂
      x_i∂ x_j - Γ_ij^l ∂ f/∂ x_l) ,

  where  
  
    Γ_ij^l =   1/2 g^lr(∂ g_ir/∂ x_j
      + ∂ g_jr/∂ x_i - ∂ g_ij/∂
      x_r) 
      = 1/2 a_lr(∂ (a^-1)_ir/∂ x_j
      + ∂ (a^-1)_jr/∂ x_i - ∂ (a^-1)_ij/∂
      x_r) , 
         1 ≤ i,j,l ≤ d

  are the Christoffel's symbols defined by  
  ∇^ℳ_e_ie_j =
Γ^l_ije_l, for 1 ≤ i,j≤ d.

Now let us consider the level set 

    Σ=ξ^-1(0)={ x ∈ℳ=ℝ^d  | ξ(x)
      = 0∈ℝ^k}

of the C^2 function ξ : ℝ^d →ℝ^k with ξ = (ξ_1, ξ_2, ⋯, ξ_k)^T, 1 ≤ k < d. 
Applying regular value theorem <cit.>, we know that Σ is a (d-k)-dimensional submanifold of ℳ, under Assumption <ref>.

Given x ∈Σ and a vector v∈ T_xℳ, the orthogonal
projection operator (d × d matrix) P : ℝ^d →
T_xΣ is defined such that ⟨v-Pv, u⟩_g = 0,
for ∀ u∈ T_xΣ.
It is straightforward to verify that P=-
a∇ξΨ^-1∇ξ^T, or entry-wise, 

    P_ij = δ_ij -
      (Ψ^-1)_αγ (a ∇ξ_α)_i
      ∂_jξ_γ ,    1 ≤ i,j ≤ d ,

where Ψ is the invertible k × k symmetric matrix at each point x
∈Σ, given by

    Ψ_αγ = ⟨^ℳξ_α, ^ℳξ_γ⟩_g 
      =(∇ξ^T a ∇ξ)_αγ ,    1 ≤α, γ≤
      k .

In the above, ∇ξ denotes the d × k matrix with entries
(∇ξ)_iα =∂_iξ_α, for 1 ≤α≤ k, 1 ≤ i
≤ d. We can verify that 

    aP^T = Pa,    P^2 = P ,    P^T∇ξ_α = 0 ,   
     1 ≤α≤ k .

  Let us further assume that v∈ T_xΣ is a tangent vector of Σ at x. 
  Since {σ_i}_1 ≤ i ≤ d forms an orthonormal basis of the
  tangent space T_x ℳ,
we have v =
⟨v, σ_i⟩_g σ_i. 
Using the fact that Pv = v, we obtain
v
= ⟨v,
p_i⟩_g p_i,
where p_i = Pσ_i ∈ T_xΣ. If we denote p_i = P_i,je_j, then it follows from (<ref>) and (<ref>) that

    P_i,j = (Pσ)_ji = σ_ji -
        (Ψ^-1)_αγ (a ∇ξ_α)_j
      (σ^T∇ξ_γ)_i  
       P_l,iP_l,j = (Pa)_ij=(aP^T)_ij = a_ij - 
         (Ψ^-1)_αγ (a∇ξ_α)_i
      (a∇ξ_γ)_j ,

for 1≤ i,j ≤ d.

Let ^Σ, ^Σ, Δ^Σ, ^Σ denote the 
gradient operator, the divergence operator, the Laplace-Beltrami operator and the Hessian operator on Σ,
respectively. It is direct to check that
the Levi-Civita connection and the gradient operator on Σ 
are given by ∇^Σ=P∇^ℳ and
^Σ=P ^ℳ, respectively. 
 In particular, for f ∈ C^∞(Σ) and let f be its
 extension to ℳ such that f∈
 C^∞(ℳ) and f|_Σ = f, we have 

    ^Σ f =
      P^ℳf = Pa ∇f .

Let ν_g be the surface measure on Σ induced from
the metric g on ℳ.
We recall that the mean curvature vector H on Σ is defined such
that <cit.>

    ∫_Σ^Σv  dν_g = - ∫_Σ⟨ H, v⟩_g  
      dν_g  ,

for all vector fields v on ℳ. 

We have the following lemma, concerning the operators on Σ. 

  Let f ∈ C^∞(Σ) and f∈
  C^∞(ℳ) be its extension
  to ℳ. u∈Γ(Tℳ) is a tangent vector field
  on ℳ and we recall the vectors p_i = Pσ_i, 1
  ≤ i ≤ d.
We have  
  

  * ^Σu = ⟨∇^ℳ_p_iu,
  p_i⟩_g.

  * (^Σp_i)p_i+ P∇^ℳ_p_ip_i= 0.


  * Δ^Σ f = ∑_i=1^dp_i^2 f + (^Σp_i) p_i
    f
    = ∑_i=1^dp_i^2 f - (P
    ∇^ℳ_p_ip_i) f.

  * Δ^Σ f  
  = ^ℳf( p_i, p_i) + Hf,
where H is the mean curvature vector of the submanifold Σ.

  * In the special case when g=a=,
  we have (^Σp_i)p_i=
  P∇^ℳ_p_ip_i= 0, and Δ^Σ = ∑_i=1^d p_i^2.
  
  


  The first two assertions can be directly verified. Let us prove the last three assertions.
Let x ∈Σ and assume that v_i, 1 ≤ i ≤ d-k, is an orthonormal basis of T_xΣ. 
We have v_i = ⟨v_i, p_j⟩_g p_j.
For the third assertion, by definition, 

    Δ^Σ f =   ^Σ(^Σ f) =
      ^Σ(P ^ℳf ) = ^Σ(⟨^ℳf,
        p_i⟩_g p_i) 
    
        =   ^Σ((p_if )p_i) =
        ∑_i=1^dp_i^2 f + (^Σp_i) p_i
        f
        = ∑_i=1^dp_i^2 f - (P
        ∇^ℳ_p_ip_i) f ,

where the second assertion has been used in the last equality.

For the fourth assertion, 
starting from the third assertion, using the definition of
^ℳ in (<ref>),
and applying Proposition <ref> below, we obtain

    Δ^Σ f =    ∑_i=1^d p_i^2 f - P ∇^ℳ_p_ip_i f
    
       =   ^ℳf( p_i, p_i) + 
      [(I-P)∇^ℳ_p_ip_i] f
    = ^ℳf( p_i, p_i) + Hf .


For the last assertion, when g=,
we have Γ_ij^l ≡ 0, for ∀ 1≤ i,j,l ≤ d.  
  Also, it follows from (<ref>) that P_i,j=P_ij = P_ji and P_il
  P_lj = P_ij. 
  We obtain  
  
    (^Σp_i) p_i
    
        =   ⟨∇^ℳ_p_jp_i,
        p_j⟩_g p_i
    
        =    P_j, lP_j,j' P_i,i'⟨∇^ℳ_e_l(P_i, re_r),
        e_j'⟩_g e_i'
    
        =    P_lr∂ P_ir/∂ x_l  P_ii' e_i'
    
      =   [∂ P_il/∂ x_l  P_ii' -  ∂
    P_lr/∂ x_l P_ri'] e_i'
    
    =    0 ,

  and the other assertions follow accordingly.


  Let H be the mean curvature vector defined in (<ref>) on the submanifold Σ. 
  We have 
  
    H =    (I-P)∇^ℳ_p_ip_i
    
        =    
          -(Ψ^-1)_αγ[1/2
     (Pa)_ij (a∇ξ_α)_l∂ (a^-1)_ij/∂ x_l +
     P_il∂ (a∇ξ_α)_l/∂ x_i] 
           a∇ξ_γ .

  In the special case when g=a=, we have 
  
    H = P_jl∂ P_il/∂ x_je_i= 
        -[(Ψ^-1)_αγ  P_ij ∂^2_ijξ_α]
        ∇ξ_γ .

  


  Given a tangent vector field v on ℳ, 
  from the definition of P we have v = Pv + (Ψ^-1)_αγ⟨v,
  a∇ξ_γ⟩_g  a∇ξ_α.
  Since Pv is a tangent vector field on Σ, using (<ref>) and the
  divergence theorem on Σ, we know 
  
    ∫_Σ⟨ H, v⟩_g  dν_g 
        = -∫_Σ^Σ[(I - P)v]  dν_g = 
        -∫_Σ^Σ[(Ψ^-1)_αγ⟨v, a∇ξ_γ⟩_g  a∇ξ_α]  dν_g .

  For the first expression, we notice that ⟨ 
  (I - P)v, p_i⟩_g ≡ 0, 1 ≤ i
  ≤ d. Applying Lemma <ref>, we have 
  
    -∫_Σ^Σ[(I - P)v] dν_g 
    
        =    -∫_Σ⟨∇^ℳ_p_i[(I - P)v], p_i⟩_g  dν_g 
    
        =    -∫_Σp_i⟨ (I -
      P)v, p_i⟩_g  dν_g +
      ∫_Σ⟨ (I - P)v, ∇^ℳ_p_ip_i
      ⟩_g  dν_g 
        
    
        =   ∫_Σ⟨v, (I -
        P)∇^ℳ_p_ip_i⟩_g  dν_g  .

  Comparing the last equality above with (<ref>), we
  conclude that H=(I - P)∇^ℳ_p_ip_i.
  
    For the second expression, we notice that ⟨
    a∇ξ_α, p_i⟩_g=0, 
and also recall the expressions (<ref>), (<ref>) and (<ref>).
Applying Lemma <ref>, integrating by parts, and noticing the cancellation
of some terms, we can derive 
  
    ^Σ[(Ψ^-1)_αγ⟨v, a∇ξ_γ⟩_g  a∇ξ_α]
    
        =   ⟨∇^ℳ_p_i[(Ψ^-1)_αγ⟨v, a∇ξ_γ⟩_g  a∇ξ_α], p_i⟩_g 
    
    =    (Ψ^-1)_αγ⟨v, a∇ξ_γ⟩_g ⟨∇^ℳ_p_i (a∇ξ_α), p_i⟩_g 
    
    =    (Ψ^-1)_αγ⟨v, a∇ξ_γ⟩_g  
          P_i, jP_i,l⟨∇^ℳ_e_j((a∇ξ_α)_r e_r), e_l⟩_g 
    
      =    (Ψ^-1)_αγ⟨v, a∇ξ_γ⟩_g  
          (Pa)_jl[(a∇ξ_α)_r Γ^i_jr(a^-1)_il +
    ∂ (a∇ξ_α)_r/∂ x_j (a^-1)_lr] 
    
    =    (Ψ^-1)_αγ⟨v, a∇ξ_γ⟩_g  
          (Pa)_ij[1/2
      (a∇ξ_α)_l∂ (a^-1)_ij/∂ x_l +
    ∂ (a∇ξ_α)_l/∂ x_i (a^-1)_lj] .

  The second identity in (<ref>) is obtained after comparing the above expression with (<ref>) .

In the case g=a=, we have Γ_il^r ≡ 0, 
1 ≤ i,l,r ≤ d. 
It follows that

    ^Σ[(Ψ^-1)_αγ⟨v, a∇ξ_γ⟩_g  a∇ξ_α]
      = (Ψ^-1)_αγ⟨v, a∇ξ_γ⟩_g 
          P_ij ∂^2_ijξ_α

and we obtain that 
H=-[(Ψ^-1)_αγ  P_ij ∂^2_ijξ_α] ∇ξ_γ.
Using (<ref>) and (<ref>), we have 
    
    P_jl∂ P_il/∂ x_je_i
          = -P_jl∂((Ψ^-1)_αγ∂_lξ_α∂_iξ_γ)/∂ x_je_i
          =- [(Ψ^-1)_αγ  P_jl ∂^2_jlξ_α ∂_iξ_γ]
          e_i
          =H ,

and therefore the first expression in (<ref>) holds as well.

Next, we study the Laplace-Beltrami operator Δ^Σ on the submanifold
Σ. Clearly, Δ^Σ is self-adjoint and, similar to (<ref>), 
we have the integration by parts formula on Σ with respect to the
measure ν_g, as 

    ∫_Σ (Δ^Σ f_1) f_2 dν_g = - ∫_Σ⟨^Σ f_1, ^Σ f_2⟩_g  dν_g
      = ∫_Σ (Δ^Σ f_2)f_1 dν_g ,

for ∀ f_1,  f_2 ∈ C^∞(Σ).
The expression of Δ^Σ can be computed explicitly and this is the
content of the following proposition. 

  Let Σ be the submanifold of ℳ defined in
  (<ref>), P be the projection matrix in (<ref>),
  and Δ^Σ be the Laplace-Beltrami operator on Σ.   We have 
  
    Δ^Σ
        =    
        (Pa)_ij∂^2 /∂ x_i∂ x_j 
        + [∂ (Pa)_ij/∂ x_j + 1/2 (Pa)_ij∂ln(( a)^-1(∇ξ^T a ∇ξ))/∂
    x_j] ∂/∂ x_i .

  In the special case when g=a=, we have 
  
    Δ^Σ =   ∑_i=1^d p^2_i=  
        P_ij∂^2 /∂ x_i∂ x_j
    + P_lj∂ P_li/∂ x_j∂/∂ x_i
    
    =    P_ij∂^2 /∂ x_i∂ x_j
          + H_i ∂/∂ x_i ,

  where H = H_i e_i is the mean curvature vector of the submanifold
  Σ.  
  


  Let f ∈ C^∞(Σ) and f∈
  C^∞(ℳ) be its extension to ℳ.
  Using Lemma <ref> and Proposition <ref>, we have 
  
    Δ^Σ f  
        =   ^ℳf( p_r, p_r) + Hf
    
      =    P_r, jP_r,l^ℳf(
        e_j,
        e_l) + Hf
    
        =    (Pa)_jl(∂^2 f/∂
    x_j∂ x_l - Γ_jl^i ∂f/∂
    x_i)
       - 
    (Ψ^-1)_αγ[1/2
          (Pa)_jl (a∇ξ_γ)_r∂ (a^-1)_jl/∂ x_r +
          P_lr∂ (a∇ξ_γ)_r/∂ x_l] 
    (a∇ξ_α)_i ∂f/∂ x_i .

  Notice that we have already obtained the coefficients of the second order
  derivative terms. For the terms of the first order derivatives, let us denote 
  
    I_1=     -(Pa)_jlΓ_jl^i
    
        I_2 =    
    -1/2
        (Ψ^-1)_αγ
          (Pa)_jl
    (a∇ξ_α)_i (a∇ξ_γ)_r∂ (a^-1)_jl/∂ x_r
    
    I_3 =    - (Ψ^-1)_αγ
    P_lr (a∇ξ_α)_i ∂
          (a∇ξ_γ)_r/∂ x_l .

    Using the expression of Pa in (<ref>), the property Pa∇ξ_γ = 0, and integrating by parts,
    we easily obtain
    
    I_2=   1/2((Pa)_ir - a_ir) (Pa)_jl∂
      (a^-1)_jl/∂ x_r
    
      I_3 =   ∂ P_lr/∂ x_l(a_ir - (Pa)_ir) .

    For I_1, direct calculation using (<ref>) gives 
    
    I_1 =    - 1/2
        (Pa)_jl
        a_ir(∂ (a^-1)_lr/∂ x_j
        + ∂ (a^-1)_jr/∂ x_l - ∂
        (a^-1)_jl/∂ x_r)
    
      =    - (Pa)_jl  a_ir∂ (a^-1)_lr/∂ x_j
        + 1/2
        (Pa)_jl
      a_ir∂ (a^-1)_jl/∂ x_r
    
      =   ∂ (Pa)_ij/∂ x_j 
    - ∂ P_jr/∂ x_j a_ir
        + 1/2
        (Pa)_jl
      a_ir∂ (a^-1)_jl/∂ x_r .

    Therefore, 
    
    I_1 + I_2 + I_3
          =    ∂ (Pa)_ij/∂ x_j 
    -∂ P_lr/∂ x_l(Pa)_ir
        + 1/2
        (Pa)_jl
      a_ir∂ (a^-1)_jl/∂ x_r + 
    1/2((Pa)_ir - a_ir) (Pa)_jl∂
      (a^-1)_jl/∂ x_r
    
    =     ∂ (Pa)_ij/∂ x_j 
    -∂ P_lr/∂ x_l(Pa)_ir
    + 1/2(Pa)_ir (Pa)_jl∂ (a^-1)_jl/∂ x_r .

    Applying Lemma <ref> below to handle the last term above, we conclude
  
    I_1 + I_2 + I_3 =   ∂ (Pa)_ij/∂ x_j 
     -1/2(Pa)_ir∂ln a/∂ x_r
     +1/2(Pa)_ir∂ln Ψ/∂ x_r .


Finally, when g=a=, applying Lemma <ref>, we can obtain
  
    Δ^Σ f 
    =   ^ℳf( p_i, p_i) +
        Hf
    
      =    P_l, iP_l,j(∂^2 f/∂ x_i∂ x_j - Γ_ij^l ∂f/∂ x_l)
      + Hf
    
      =    P_ij∂^2 f/∂ x_i∂ x_j + H_i
      ∂f/∂ x_i .

The other equality in (<ref>) follows from Proposition <ref>.

We point out that the proof of Proposition <ref>
is indeed valid for a general Riemannian manifold ℳ and its level set Σ as well. In this case, 
(<ref>) holds true on a local coordinate of the manifold ℳ. 

The following identity has been used in the above proof, and will be
useful in Appendix <ref> as well.

  
    1/2(Pa)_ir  (Pa)_jl∂ (a^-1)_jl/∂ x_r 
         = -1/2(Pa)_ir∂ln a/∂ x_r
    + (Pa)_ir∂ P_lr/∂ x_l
    + 1/2(Pa)_ir∂lnΨ/∂ x_r .

    


  Using the expression of Pa in (<ref>), the relations 
  
    Pa∇ξ_γ = 0 ,  ∂ln a/∂ x_r = (a^-1)_jl∂ a_jl/∂ x_r ,   ∂lnΨ/∂ x_r = (Ψ^-1)_αγ∂Ψ_αγ/∂ x_r ,

and the integration by parts formula, we can compute
  
    1/2(Pa)_ir  (Pa)_jl∂ (a^-1)_jl/∂ x_r
    
      =    1/2(Pa)_ir(a_jl - 
      (Ψ^-1)_αγ (a∇ξ_α)_j (a∇ξ_γ)_l)
      ∂ (a^-1)_jl/∂ x_r
    
      =   
    -1/2(Pa)_ir∂ln a/∂ x_r
    - 1/2(Pa)_ir  
    (Ψ^-1)_αγ (a∇ξ_α)_l ∂^2_lrξ_γ 
    + 1/2(Pa)_ir  
    (Ψ^-1)_αγ∂_lξ_α ∂ (a∇ξ_γ)_l/∂ x_r
    
      =   
    -1/2(Pa)_ir∂ln a/∂ x_r
    - (Pa)_ir  
    (Ψ^-1)_αγ (a∇ξ_α)_l ∂^2_lrξ_γ 
    + 1/2(Pa)_ir∂lnΨ/∂ x_r
    
      =   
    -1/2(Pa)_ir∂ln a/∂ x_r
    + (Pa)_ir∂ P_lr/∂ x_l
    + 1/2(Pa)_ir∂lnΨ/∂ x_r .

  
  We conclude this section with the proof of Corollary <ref> in Section <ref>.

  Notice that the infinitesimal generator of (<ref>)
  can be written as 
  
    ℒ^J = J_i∂/∂ x_i + ℒ ,

  where ℒ is the infinitesimal generator of (<ref>).
  Using the fact J∈ T_x Σ, the same argument of Proposition <ref> implies that (<ref>) evolves on Σ as well. 
  Since μ_1 is invariant with respect to ℒ, to show
  the SDE (<ref>) has the same invariant measure, it is enough to verify that 
  
    ^Σ{Jexp[-β(U+1/2βln (∇ξ^T a ∇ξ)/ a)]} = 0 ,   ∀ x ∈Σ ,

  where we have used the expression of μ_1 in (<ref>). Applying
  the formula of ^Σ in Lemma <ref>, 
  we can compute the right hand side of (<ref>), as 
  
    ^Σ{Jexp[-β(U+1/2βln (∇ξ^T a ∇ξ)/ a)]}
    
        =   ⟨∇^ℳ_p_j{exp[-β(U+1/2βln (∇ξ^T a ∇ξ)/ a)]J_ie_i},
        p_j⟩_g 
    
          =    P_j,lP_j,r⟨∇^ℳ_e_l{exp[-β(U+1/2βln (∇ξ^T a ∇ξ)/ a)]J_ie_i},
          e_r⟩_g .

  which implies that (<ref>) is equivalent to 
  
    0 =    (Pa)_lr∂ J_i/∂ x_l (a^-1)_ir +
        (Pa)_lrJ_i Γ_li^r' (a^-1)_r'r - β (Pa)_lr J_i (a^-1)_ir∂ U/∂ x_l
        - 1/2 (Pa)_lr J_i (a^-1)_ir∂/∂
        x_l[ ln (∇ξ^T a ∇ξ)/ a]
    
    =    P_li∂ J_i/∂ x_l +
        P_lr'J_i Γ_li^r' - β P_li J_i ∂ U/∂ x_l 
         - 1/2 P_li J_i ∂/∂ x_l[ ln (∇ξ^T a ∇ξ)/ a] ,

  where Γ_li^r' are the Christoffel's symbols satisfying ∇^ℳ_e_le_i =
  Γ^r'_lie_r'. Using the expression (<ref>) of Γ_li^r', the fact J_i = P_ij J_j, and
  Lemma <ref>, we can
  further simplify the above equation 
  and obtain
  
    P_lr'J_i Γ_li^r' - 1/2 P_li J_i
        ∂/∂ x_l[ ln (∇ξ^T a ∇ξ)/ a]
    
        =   1/2 J_i(Pa)_lr∂ (a^-1)_lr/∂ x_i -
        1/2 P_li J_i ∂/∂ x_l[ ln
        (∇ξ^T a ∇ξ)/ a]
    
        =    J_j ∂ P_ij/∂ x_i .

  Therefore, we see that 
  (<ref>) is equivalent to the condition (<ref>).





§ PROOFS IN SECTION <REF>


In this section, we collect proofs of the various results in Section <ref>. 

First, we prove Proposition <ref>, which concerns the properties of the flow map Θ defined in 
(<ref>), (<ref>), and (<ref>).
While the approach of the proof is similar to the one in <cit.>,
here we consider the specific function F in (<ref>) 
and we will provide full details of the derivations.


      In this proof, we will always assume x ∈Σ. 
      For a function which only depends on the state and is evaluated at x,
      we will often omit its argument in order to keep the notations simple.
      Also notice that, repeated indices other than l and l'
      indicate that they are summed up, while for the indices l, l'
      we assume that they are fixed by default unless the summation operator is used explicitly.

      Since ∇
      F = 0 on Σ, from the equation (<ref>) we know that
      φ(x,s) ≡ x, ∀ s ≥ 0. 
      Let us Denote by ∇^2 F the Hessian matrix (on the standard Euclidean space)
      of the function F in (<ref>), i.e., 
     ∇^2 F =(∂^2_ij F)_1 ≤ i, j ≤ d.
      Since ξ(x) = 0∈ℝ^k, direct calculation gives 
    
    (a∇^2 F)_ij = a_ir∂^2 F/∂
          x_r∂ x_j = 
          (a∇ξ∇ξ^T)_ij ,      1 ≤ i, j ≤ d .

    Meanwhile, it is straightforward to verify that a∇^2 F satisfies 
    
    ⟨ a ∇^2 F u, v⟩_g = ⟨u, a
          ∇^2 F v⟩_g   ,  ∀ u, v∈ℝ^d , 
       ⟨ a ∇^2F u, u⟩_g = |∇ξ^Tu|^2 ≥
          0 ,   ∀ u∈ℝ^d ,
       (a∇^2 F) u= a∇ξ∇ξ^T u = 0 ,
             ∀ u∈ T_xΣ .

    Therefore, we can assume that a∇^2F has real (non-negative) eigenvalues 
	  
    λ_1=λ_2= ⋯=λ_d-k = 0 <
        λ_d-k+1≤⋯≤λ_d ,

     and the corresponding eigenvectors, denoted by  
    v_i = (v_i1, v_i2, ⋯, v_id)^T, 1 ≤ i ≤ d, 
    are orthonormal with respect to the inner product
    ⟨·,·⟩_g in (<ref>), such that
v_1, v_2, ⋯, v_d-k∈ T_xΣ.

The projection matrix P in (<ref>) can be expressed using the vectors 
v_i as
	  
    P_ij = ∑_l=1^d-k v_li (a^-1)_jr v_lr ,    1 ≤ i, j ≤ d ,

	  and we have 
	  
    ∑_l=1^d-k v_liv_lj = (Pa)_ij ,    a_ij - (Pa)_ij
    	    =∑_l=d-k+1^d v_liv_lj .

    It is also a simple fact that the eigenvalues of the k× k matrix
    Ψ=∇ξ^T a∇ξ are λ_d-k+1, λ_d-k+2,
    ⋯, λ_d, 
    with the corresponding eigenvectors given by ∇ξ^T v_d-k+1,
    ∇ξ^Tv_d-k+2, ⋯, ∇ξ^Tv_d. In
    particular, this implies
      
    ∏_i=d-k+1^d λ_i= (∇ξ^T a∇ξ) = Ψ .


      In the following, we study the ODE (<ref>) using the eigenvectors v_i.
    Differentiating the ODE (<ref>) twice, using the facts that 
    φ(x,s)≡ x, ∀ s ≥ 0, and ∇ F = 0 on Σ, we obtain     
	  
    d/ds∂φ_i/∂ x_j(x,s) =   
    	-(a_ir'∂^2 F/∂ x_r'∂ x_i')  ∂ φ_i'/∂ x_j (x,s)
    d/ds∂^2 φ_i/∂ x_j∂
          x_r(x,s) =   
    	-(2∂ a_ir'/∂ x_i'∂^2
    	  F/∂ x_r'∂ x_j' + a_ir'∂^3
    	  F/∂ x_r'∂ x_i'∂
    	    x_j') 
          ∂ φ_i'/∂ x_j (x,s)
          ∂ φ_j'/∂ x_r (x,s)
       -(a_ir'∂^2 F/∂ x_r'∂ x_i')
          ∂^2 φ_i'/∂ x_j ∂ x_r(x,s) ,

for s ≥ 0 and 1 ≤ i, j, r ≤ d.
      
	
  * The first equation of (<ref>) implies 
    
    d/ds(v_lj∂φ_i/∂
          x_j(x,s)) =   
    	-(a_ir'∂^2 F/∂ x_r'∂ x_i')
          ( v_lj∂ φ_i'/∂
          x_j(x,s)) ,   1 ≤ l ≤ d .

    Since φ(·, 0) is the identity map, we have 
    
    v_lj∂φ_i/∂
        x_j(x,0) = v_li,     s=0 .

  Because v_l is the eigenvector of a∇^2 F, we can 
  directly solve the solution of (<ref>)-(<ref>) and obtain
    
    v_lj∂ φ_i/∂ x_j(x,s) =
          e^-λ_l s v_li ⟺∂ φ_i/∂ x_j(x,s) =
    ∑_l=1^d e^-λ_l s v_li (a^-1)_jrv_lr ,     ∀ s ≥ 0 ,

    for 1 ≤ i, j ≤ d.
	  Sending s→ +∞, using (<ref>) and (<ref>), we obtain 
    
    ∂Θ_i/∂ x_j = lim_s→ +∞∂ φ_i/∂ x_j(x,s) = ∑_l=1^d-k v_li 
          (a^-1)_jrv_lr =P_ij .

	
  * We proceed to compute a_jr∂^2Θ_i/∂
	  x_j∂ x_r, 1 ≤ i≤ d. For this purpose,
    let us define
    
    A_l(x,s) =  (a^-1)_ij' v_lj' a_jr∂^2 φ_i/∂ x_j∂ x_r(x,s)  ,   1 ≤ l
          ≤ d ,
    ⟺     
           a_jr∂^2 φ_i/∂ x_j∂
           x_r(x,s) = ∑_l=1^d v_li A_l(x,s)  .

    Using the second equation of (<ref>), 
    the solution (<ref>), and the orthogonality
    of the eigenvectors, we can obtain
    
    d A_l/ds(x,s)
          =
          - ∑_l'=1^d[
    	2∂ a_ir/∂ x_i'∂^2 F/∂
    	x_r∂ x_j(a^-1)_ir' +       
          ∂^3 F/∂ x_r'∂ x_i'∂ x_j]
          v_l'i'v_l'j v_lr'  e^-2λ_l's 
          - λ_l A_l(x,s) ,

for 1 ≤ l ≤ d, from which we get 
    
    P_ii' 
          a_jr∂^2 φ_i'/∂ x_j∂ x_r(x,s) 
    
          =   ∑_l=1^d-k v_liA_l(x,s) 
    
          =    
          -∑_l=1^d-k∑_l'=1^d[
    	2∂ a_jr/∂ x_i'∂^2 F/∂
    	x_r∂ x_j'(a^-1)_jr' +       
          ∂^3 F/∂
        x_r'∂ x_i'∂ x_j']
        v_l'i'v_l'j'v_lr'v_li 
        e^-λ_l s∫_0^s e^(λ_l -2λ_l')u du 
    
          =    ∑_l=1^d-k∑_l'=1^d
    [2λ_l'∂ (a^-1)_i'r/∂ x_j
          - 
          ∂^3 F/∂
        x_r∂ x_i'∂ x_j]
          v_l'i'v_l'jv_lrv_li ∫_0^s e^-2λ_l' u du .

To further simplify the last expression above, 
 we differentiate the identity
    
    ∂^2 F/∂ x_i'∂ x_j
          v_l'i'v_l'j = λ_l' ,

where l' is fixed, 1 ≤ l' ≤ d, along the eigenvector v_l, which
gives
    
    ∂^3 F/∂ x_r∂ x_i'∂ x_j
          v_l'i'v_l'jv_lr
    
          =    -2 ∂^2 F/∂ x_i'∂ x_j∂ v_l'i'/∂ x_r v_l'j
          v_lr +
          ∂λ_l'/∂ x_r v_lr
    
          =   -2 λ_l'
           (a^-1)_i'r' v_l'r'∂ v_l'i'/∂ x_r 
          v_lr +
          ∂λ_l'/∂ x_r v_lr .

    Therefore, taking the limit s→ +∞, using 
    the relations (<ref>),
    (<ref>), and
    Lemma <ref> in Appendix <ref>, we can compute
    
    P_ii' a_jr∂^2 Θ_i'/∂ x_j∂ x_r
    
          =   lim_s→ +∞ P_ii' 
          a_jr∂^2 φ_i'/∂ x_j∂ x_r(x,s)
    
          =    lim_s→ +∞∑_l'=d-k+1^d∑_l=1^d-k[2λ_l'∂ (a^-1)_i'r/∂ x_j
    v_l'i'v_l'j
    +2 λ_l'
           (a^-1)_jr' v_l'r'∂ v_l'j/∂ x_r -
          ∂λ_l'/∂ x_r]
          v_lrv_li ∫_0^s e^-2λ_l' u du
    
          =   ∂ (a^-1)_i'r/∂ x_j(a_i'j -
          (Pa)_i'j) (Pa)_ir - 
          1/2∂ (a^-1)_jr'/∂ x_r(a_jr' -
          (Pa)_jr') (Pa)_ir
          - 1/2 (Pa)_ir∂lnΨ/∂ x_r
    
          =   
          -P_ir∂ a_jr/∂ x_j - ∂
    	P_jr/∂ x_j (Pa)_ir + P_ir∂ (Pa)_jr/∂ x_j 
    + ∂ P_jr/∂ x_j (Pa)_ir
    
    =    
          -P_ir∂ a_jr/∂ x_j + P_ir∂ (Pa)_jr/∂ x_j .


    On the other hand, differentiating the relation ξ(Θ(x)) ≡0 twice and using (<ref>), we get
    
    ∂ξ_γ/∂ x_i'∂^2
          Θ_i'/∂ x_j∂ x_r
          =-
          ∂^2
          ξ_γ/∂ x_i'∂ x_j'∂Θ_i'/∂ x_j∂Θ_j'/∂ x_r  = 
          - ∂^2
          ξ_γ/∂ x_i'∂ x_j' P_i'jP_j'r ,

       for 1 ≤γ≤ k.
    Therefore, using PaP^T=P^2a=Pa and Pa∇ξ_γ=0, we
    can compute
    
    (δ_ii'-P_ii')  
          a_jr∂^2 Θ_i'/∂ x_j∂ x_r
    
          =       (Ψ^-1)_αγ  (a∇ξ_α)_i
          ∂ξ_γ/∂ x_i' 
          a_jr∂^2 Θ_i'/∂ x_j∂ x_r
    
          =    
    - (Ψ^-1)_αγ
    (a∇ξ_α)_i
    ∂^2 ξ_γ/∂ x_i'∂ x_j' 
    P_i'jP_j'r a_jr
    
     =   
    - (Ψ^-1)_αγ
    (a∇ξ_α)_i
    (∂^2_i'j'ξ_γ)(Pa)_i'j'
    
          =    
    (Pa)_i'j'∂ P_ii'/∂ x_j' .

    Summing up (<ref>) and (<ref>), we conclude
    that
    
    a_jr∂^2 Θ_i/∂ x_j∂ x_r 
          = 
          ∂ (Pa)_ij/∂ x_j - 
          P_ir∂ a_rj/∂ x_j .

      
    
    Now, we prove Theorem <ref>.

  Since we follow the approach in <cit.>, we will only sketch the proof and will mainly focus on the differences.

  First of all, we introduce some notations. Let x^(l), l = 0, 1, ⋯, be the states generated from the
  numerical scheme (<ref>) and let 
  ψ be a function on Σ. We will adopt the abbreviations ψ^(l) = ψ(x^(l)), P^(l) = P(x^(l)), etc. 
For j≥ 1, D^jψ[u_1, u_2, ⋯, u_j] denotes the 
  jth order directional derivatives of ψ along the vectors u_1, u_2,
  ⋯, u_j, and |D^jψ|_∞ is the supremum norm of D^jψ on Σ.
  Similarly, D^jΘ[u_1,
  u_2, ⋯, u_j] denotes the d-dimensional vector whose ith component is 
D^jΘ_i[u_1, u_2, ⋯, u_j], for 1 ≤ i ≤ d.

Define the vector b^(l)=(b^(l)_1, b^(l)_2, ⋯, b^(l)_d)^T by 
  
    b_i^(l) = (-a_ij∂ U/∂ x_j   +
      1/β∂ a_ij/∂ x_j)(x^(l)) ,    1
      ≤ i ≤ d ,
 
  for l = 0, 1, ⋯, and set 
  
    δ^(l) = b^(l)h + √(2β^-1 h) σ^(l)η^(l) .

 We have 

    δ^(l) = x^(l+1/2)- x^(l) ,      
      x^(l+1) = Θ(x^(l+1/2)) =
    Θ(x^(l)+δ^(l)) .

Let ℒ be the infinitesimal generator of the SDE
  (<ref>) in Theorem <ref>, given by  

    ℒ =    - (Pa)_ij∂ U/∂
      x_j∂/∂ x_i + 1/β∂
      (Pa)_ij/∂ x_j∂/∂ x_i + 1/β
      (Pa)_ij∂^2/∂ x_i∂ x_j
    
      =   e^β U/β∂/∂ x_i(e^-β
      U(Pa)_ij∂/∂ x_j) ,

in Remark <ref>.
We consider the Poisson equation on Σ

    ℒψ = f - f .

The existence and the regularity of the solution ψ can be established under
  Assumption <ref>–<ref>, and the Bakry-Emery condition in
  Section <ref>.
Applying Taylor's theorem and using the fact that Θ(x^(l)) = x^(l) since x^(l)∈Σ, we have

    ψ^(l+1)
    
        =    (ψ∘Θ)(x^(l) + δ^(l))
    
        =   ψ^(l) + D(ψ∘Θ)^(l)[δ^(l)]
        + 1/2D^2(ψ∘Θ)^(l)[δ^(l), δ^(l)]
        +
        1/6D^3(ψ∘Θ)^(l)[δ^(l),δ^(l),δ^(l)]
        + R^(l) ,

where the reminder is given by 

    R^(l) = 
      1/6(∫_0^1 s^3 D^4(ψ∘Θ)(x^(l)+(1-s)δ^(l))
      ds)[δ^(l),δ^(l),δ^(l),δ^(l)] .


Now we apply Proposition <ref> to 
simplify the expression in (<ref>). Using the chain rule, the expressions (<ref>)–(<ref>), we can derive 

    ψ^(l+1)
    
      =    ψ^(l) + Dψ^(l)[P^(l)δ^(l)+1/2
    D^2Θ^(l)[δ^(l),δ^(l)]]
      + 1/2D^2ψ^(l)[P^(l)δ^(l), P^(l)δ^(l)]
       + 1/6D^3(ψ∘Θ)^(l)[δ^(l),δ^(l),δ^(l)] + R^(l)
    
      =    ψ^(l) + (ℒψ)^(l) h+ √(2β^-1h)  Dψ^(l)[(Pσ)^(l)η^(l)]
      + 
      h^2/2 Dψ^(l)[D^2Θ^(l)[
      b^(l), b^(l)]]
        + √(2β^-1) h^3/2 Dψ^(l)[D^2Θ^(l)[
    b^(l), σ^(l)η^(l)]] +
    h^2/2D^2ψ^(l)[P^(l)b^(l),P^(l)b^(l)]
       + √(2β^-1) h^3/2 D^2ψ^(l)[P^(l)b^(l), (Pσ)^(l)η^(l)] 
      + hDψ^(l)/β[D^2Θ^(l)[σ^(l)η^(l),σ^(l)η^(l)]-a^(l):∇^2 Θ^(l)]
         + h/β(D^2ψ^(l)[(Pσ)^(l)η^(l),
      (Pσ)^(l)η^(l)] -
      (Pa)^(l):D^2ψ^(l))
     + 1/6D^3(ψ∘Θ)^(l)[δ^(l),δ^(l),δ^(l)]
    + R^(l) ,

where in the last equation we added and subtracted some terms, and 
we used the identity 

    Dψ^(l)[P^(l)b^(l)+1/β
      a^(l):∇^2 Θ^(l)] + 1/β (Pa)^(l):∇^2ψ^(l) =
      (ℒψ)^(l) ,

which can be verified using Proposition <ref>,
(<ref>) and (<ref>).
In (<ref>), a:∇^2Θ is the vector whose ith component is given by
a_jr∂^2 Θ_i/∂ x_j∂ x_r, and
(Pa):∇^2 ψ is defined in a similar way. 

Summing up (<ref>) for l=0, 1, ⋯, n-1, dividing both
sides by T,
and using the Poisson equation (<ref>), gives

    f_n - f = 
    1/n∑_l=0^n-1 f(x^(l)) -
      f 
      =
      ψ^(n)-ψ^(0)/T + 1/T∑_i=1^5 M_i,n +
      1/T∑_i=1^4 S_i,n ,

where 

    M_1,n =    -√(2β^-1h)∑_l=0^n-1Dψ^(l)[(Pσ)^(l)η^(l)]  ,
    
        M_2,n =    -√(2β^-1)
      h^3/2∑_l=0^n-1 Dψ^(l)[D^2Θ^(l)[b^(l),
      σ^(l)η^(l)]] , 
    
      M_3,n =    -h/β∑_l=0^n-1
      Dψ^(l)[D^2Θ^(l)[σ^(l)η^(l),σ^(l)η^(l)]-
        a^(l):∇^2 Θ^(l)] ,
    
         M_4,n =    -
      √(2β^-1) h^3/2∑_l=0^n-1
      D^2ψ^(l)[P^(l)b^(l), (Pσ)^(l)η^(l)] , 
    
        M_5,n =   
      -h/β∑_l=0^n-1(D^2ψ^(l)[(Pσ)^(l)η^(l),
      (Pσ)^(l)η^(l)] - (Pa)^(l):D^2ψ^(l)) ,

and

    S_1,n =    -h^2/2∑_l=0^n-1
      Dψ^(l)[D^2Θ^(l)[b^(l),b^(l)]]  , 
    
        S_2,n =    -h^2/2∑_l=0^n-1
      D^2ψ^(l)[P^(l)b^(l), P^(l)b^(l)] ,
      
    
      S_3,n =    -∑_l=0^n-1 R^(l) ,     
      S_4,n = -1/6∑_l=0^n-1
      D^3(ψ∘Θ)^(l)[δ^(l),δ^(l),δ^(l)] .

Using (<ref>),  the last term S_4,n above 
can be further decomposed as

    S_4,n = M_0,n + S_0,n ,

where 

    M_0,n =    -√(2β^-1)/6 h^3/2∑_l=0^n-1(2/β
        D^3(ψ∘Θ)^(l)[σ^(l)η^(l),σ^(l)η^(l),σ^(l)η^(l)]+ 
        3 h D^3(ψ∘Θ)^(l)[b^(l),b^(l),
        σ^(l)η^(l)]) , 
    
        S_0,n =    - h^2/6∑_l=0^n-1(6/β
        D^3(ψ∘Θ)^(l)[b^(l),
        σ^(l)η^(l),σ^(l)η^(l)] + hD^3(ψ∘Θ)^(l)[
        b^(l),b^(l),b^(l)]) .

Notice that the terms M_i,n, 0 ≤ i ≤ 5, are all martingales and in
particular we have 𝐄 M_i,n = 0. Therefore, since the level set
Σ is compact (Assumption <ref>), the first conclusion follows from the estimates 

    |S_1,n| ≤ C|Dψ|_∞ h T ,      |S_2,n| ≤
       C|D^2ψ|_∞  h T ,
       𝐄|S_0,n| ≤ C|D^3ψ|_∞  h T  ,    𝐄|S_3,n| ≤ C|D^4ψ|_∞ h T  ,

while the second conclusion follows by squaring both sides of (<ref>) and using
the estimates 

    𝐄|S_0,n|^2 ≤ C h^2 T^2 |D^3ψ|^2_∞ ,     𝐄|S_3,n|^2 ≤ C h^2 T^2 |D^4ψ|^2_∞ ,    
       𝐄 |M_0,n|^2 ≤ Ch^2 T|D^3ψ|^2_∞ ,     𝐄 |M_1,n|^2 ≤ CT|Dψ|^2_∞ ,    𝐄 |M_2,n|^2 ≤ Ch^2 T|Dψ|^2_∞ ,
       𝐄 |M_3,n|^2 ≤ ChT|Dψ|^2_∞ ,      𝐄 |M_4,n|^2 ≤ Ch^2 T|D^2ψ|^2_∞ ,    𝐄 |M_5,n|^2 ≤ ChT|D^2ψ|^2_∞ .


As far as the third conclusion (pathwise estimate) is concerned, notice that (<ref>) implies

    | f_n - f|
      ≤   |ψ^(n)-ψ^(0)|/T + 1/T∑_i=0^5 |M_i,n| + 1/T∑_i=0^3 |S_i,n| 
    ≤   
      C(h + 1/T) + 1/T∑_i=0^5 |M_i,n| ,

where we have used the estimates (<ref>) for |S_1,n|, |S_2,n|, and the upper bounds

    |S_0,n|≤    Ch^2 ∑_l=0^n-1 |η^(l)|^2 + Ch^3 n ≤ ChT ,
           a.s. 
    
      |S_3,n|≤    Ch^2 ∑_l=0^n-1 |η^(l)|^4 + Ch^4 n ≤ ChT ,
           a.s.

which are implied by the strong law of large numbers for 1/n∑_l=0^n-1 |η^(l)|^4, when n → +∞.
Finally, we estimate the martingale terms M_i,n in (<ref>). Notice that, for any r ≥ 1, we can deduce the following upper bounds
(see <cit.>)

    1/T^2r𝐄|M_1,n|^2r≤C/T^r ,  1/T^2r𝐄|M_2,n|^2r≤Ch^2r/T^r , 
       1/T^2r𝐄|M_3,n|^2r≤Ch^r/T^r ,  1/T^2r𝐄|M_4,n|^2r≤Ch^2r/T^r , 
       1/T^2r𝐄|M_5,n|^2r≤Ch^r/T^r ,  1/T^2r𝐄|M_0,n|^2r≤Ch^2r/T^r ,

which give 

    𝐄(1/T∑_i=0^5 |M_i,n|)^2r≤C/T^2r∑_i=0^5 𝐄|M_i,n|^2r≤C/T^r .

Now, for any 0 < ϵ < 1/2, the Borel-Cantelli lemma 
implies that there is an almost surely bounded random variable ζ(ω), such that

    1/T∑_i=0^5 |M_i,n| ≤ζ(ω)/T^1/2 -
    ϵ .

Therefore, the third conclusion follows readily from 
(<ref>) and (<ref>).

Next, we prove Corollary <ref>.

  From the estimates in (<ref>), we know that it is only
  necessary to consider the term M_1,n in (<ref>).
  Recall that ψ solves the Poisson equation (<ref>) and we
  can assume ∫_Σψ dμ_1 = 0 without loss of generosity.
  Applying the Poisson equation, the Poincaré inequality, and the
  Cauchy-Schwarz inequality, we have the standard estimates
    
    ∫_Σψ^2 dμ_1 ≤    -1/K∫_Σ
        (ℒψ) ψ dμ_1 
    ≤   1/K[∫_Σ (ℒψ)^2
        dμ_1]^1/2( ∫_Σψ^2 dμ_1)^1/2
    
        =   1/K[∫_Σ (f-f)^2 dμ_1]^1/2( ∫_Σψ^2 dμ_1)^1/2 ,

  which implies
  
    (∫_Σψ^2 dμ_1)^1/2≤1/K[∫_Σ (f-f )^2 dμ_1]^1/2 ,
            
    - ∫_Σ (ℒψ) ψ dμ_1 
    ≤1/K∫_Σ (f-f )^2 dμ_1 .

  Since the term M_1,n in (<ref>) is a martingale, we
  have 
  
    1/T^2𝐄 |M_1,n|^2 = 2β^-1/T1/n∑_l=0^n-1𝐄[((Pa)^(l)∇ψ^(l))
        ·∇ψ^(l)] .

  Applying the first estimate in the conclusion of
  Theorem <ref>, using
  (<ref>) in Remark <ref>, as well as the
  estimate (<ref>), we obtain
  
    1/T^2𝐄 |M_1,n|^2 ≤   2β^-1/T∫_Σ (Pa∇ψ)·∇ψ  dμ_1 + C(h/T+1/T^2)
    ≤   2∫_Σ (f-f )^2 dμ_1/KT 
        + C(h/T+1/T^2) .

  The conclusion follows by squaring both sides of
  (<ref>), applying Young's inequality, and using the same argument of Theorem <ref>.


Finally, we prove Proposition <ref>, which concerns the properties of the projection map Π defined in (<ref>).

    
  For 1 ≤ l ≤ d, recall that
  p_l = (Pσ)_i'le_i'
  is the tangent vector field defined in Appendix <ref>
  such that 
  p_l∈ T_xΣ at each x ∈Σ.
  Since Π_i(x) = x_i for
  x ∈Σ, 1 ≤ i ≤ d, taking derivatives along p_l
  twice, we obtain
  
    ∂Π_i/∂ x_j (Pσ)_jl  = (Pσ)_il ,
        
       ∂^2 Π_i/∂ x_j∂ x_r 
        (Pσ)_jl(Pσ)_rl = (Pσ)_rl∂
          (Pσ)_il/∂ x_r - ∂Π_i/∂ x_j 
          (Pσ)_rl∂ (Pσ)_jl/∂ x_r .

Notice that, for a function which only depends on the state and is evaluated at x ∈Σ,
      we will often omit its argument in order to keep the notations simple.

 On the other hand, the vector
 σ_l - p_l = ((I-P)σ)_i'le_i'∈ (T_xΣ)^⊥ (the complement
 of the subspace T_xΣ in T_xℳ). Let
 ϕ(s) be the geodesic curve in ℳ such that ϕ(0) = x and ϕ'(0) = 
 σ_l - p_l.
 We have Π_i(ϕ(s)) = x_i, ∀ s ∈ [0, ϵ) for some
   ϵ > 0. Taking derivatives with respect to s twice, we obtain
   
    ∂Π_i/∂ x_j(ϕ(s))dϕ_j(s)/ds =
         0 , 
       ∂^2 Π_i/∂ x_j∂ x_r(ϕ(s))
         dϕ_j(s)/dsdϕ_r(s)/ds=-∂Π_i/∂ x_j(ϕ(s)) d^2ϕ_j(s)/ds^2 
         =∂Π_i/∂ x_j(ϕ(s)) Γ^j_rr'(ϕ(s)) dϕ_r(s)/dsdϕ_r'(s)/ds ,

   for 1 ≤ i ≤ d, where ϕ_j denotes the jth component of ϕ, and
   the geodesic equation of the curve ϕ has been used to obtain the last expression above. In
   particular, setting s=0, we obtain 
   
    ∂Π_i/∂ x_j(σ_jl - (Pσ)_jl) = 0 , 
       ∂^2 Π_i/∂ x_j∂ x_r(σ_jl - (Pσ)_jl)
    (σ_rl - (Pσ)_rl)
          =∂Π_i/∂ x_jΓ^j_rr'(σ_rl - (Pσ)_rl)
         (σ_r'l - (Pσ)_r'l) .

   Combining the first equations in both (<ref>) and
   (<ref>), we can conclude that ∂Π_i/∂ x_j
   = P_ij at x ∈Σ. 
   Since (<ref>) holds at any x ∈Σ, taking
   the derivative in the first equation of (<ref>) along the tangent vector 
   p_l∈ T_xΣ, we obtain 
   
    ∂^2 Π_i/∂ x_j∂ x_r(σ_jl - (Pσ)_jl) (Pσ)_rl = - 
         ∂Π_i/∂ x_j 
         (Pσ)_rl∂(σ_jl - (Pσ)_jl) /∂ x_r .

   Combining (<ref>), (<ref>) and (<ref>),
   using Lemma <ref> in Appendix <ref>,
   the expression in (<ref>),
   the relations
   
    (Pσ)_jl (Pσ)_rl = (PaP^T)_jr = (Pa)_jr , 
       (σ_rl - (Pσ)_rl) (σ_r'l -
    (Pσ)_r'l)= a_rr' - (Pa)_rr' ,

   and the integration by parts formula, we can compute
   
    ∂^2 Π_i/∂ x_j∂ x_r a_jr
       ∂^2 Π_i/∂ x_j∂ x_r(Pσ + (σ - Pσ))_jl(Pσ + (σ - Pσ))_rl
    
         =    ∂^2 Π_i/∂ x_j∂ x_r 
         (Pσ)_jl (Pσ)_rl
         + 
         2 ∂^2 Π_i/∂ x_j∂ x_r 
         (σ - Pσ)_jl (Pσ)_rl
         + ∂^2 Π_i/∂ x_j∂ x_r
    (σ - Pσ)_jl(σ - Pσ)_rl
    
         =    
    (Pσ)_rl∂ (Pσ)_il/∂ x_r - P_ij 
          (Pσ)_rl∂ (Pσ)_jl/∂ x_r 
          - 2P_ij (Pσ)_rl∂(σ_jl - (Pσ)_jl) /∂ x_r
       
         + P_ijΓ^j_rr'(σ_rl - (Pσ)_rl) (σ_r'l - (Pσ)_r'l) 
    
         =   [(Pσ)_rl∂ (Pσ)_il/∂ x_r + P_ij 
          (Pσ)_rl∂ (Pσ)_jl/∂ x_r 
          - 2P_ij (Pσ)_rl∂σ_jl/∂ x_r] 
          + P_ijΓ^j_rr'(a_rr' - (Pa)_rr')
    
         =   [2(Pσ)_rl∂ (Pσ)_il/∂ x_r -  
    (Pσ)_jl
          (Pσ)_rl∂ P_ij/∂ x_r 
          - 2 (Pσ)_rl∂ (Pσ)_il/∂ x_r 
          + 2 (Pσ)_rlσ_jl∂ P_ij/∂ x_r] 
       + (Pa)_ij∂ (a^-1)_lj/∂ x_r(a_lr-(Pa)_lr)
          -1/2 (Pa)_ij∂ (a^-1)_lr/∂ x_j(a_lr-(Pa)_lr)
    
         =    
      (Pa)_lj∂ P_ij/∂ x_l 
    +[ - P_ij∂ a_jl/∂ x_l
           - (Pa)_ij∂ P_lj/∂ x_l
          + P_il∂ (Pa)_lj/∂ x_j]
          +(Pa)_ij∂ P_lj/∂ x_l + 
    	1/2 (Pa)_ij∂lnΨ/∂
    	x_j
    
         =    
     - P_ij∂ a_jl/∂ x_l
          + ∂ (Pa)_ij/∂ x_j
          +
    	1/2 (Pa)_ij∂lnΨ/∂ x_j .

 

  
plainnat




