Benoît Combès

Marc Fournier

David N Kennedy

José Braga

Neil Roberts

Sylvain Prima

EM-ICP STRATEGIES FOR JOINT MEAN SHAPE AND CORRESPONDENCES ESTIMATION: APPLICATIONS TO STATISTICAL ANALYSIS OF SHAPE AND OF ASYMMETRY

Keywords: point set, mean shape, EM-ICP, shape analysis, symmetry, surface

In this paper, we propose a new approach to compute the mean shape of unstructured, unlabelled point sets with an arbitrary number of points. This approach can be seen as an extension of the EM-ICP algorithm, where the fuzzy correspondences between each point set and the mean shape, the optimal non-linear transformations superposing them, and the mean shape itself, are iteratively estimated. Once the mean shape is computed, one can study the variability around this mean shape (e.g. using PCA) or perform statistical analysis of local anatomical characteristics (e.g. cortical thickness, asymmetry, curvature). To illustrate our method, we perform statistical shape analysis on human osseous labyrinths and statistical analysis of global cortical asymmetry on control subjects and subjects with situs inversus.

INTRODUCTION

Algorithmic tools for groupwise anatomical studies are of great interest to characterise and compare anatomical structures over and between populations in an efficient and comprehensive way. Lots of efforts have been made to design methods allowing to build statistical shape models (i.e. defining and computing a mean shape M from the observations and characterising the main variations of the shapes around M ) when the shapes are represented by implicit surfaces, meshes or simple point sets. Recent works include methods based on i) the minimum description length (MDL) theory [START_REF] Heimann | 3d active shape models using gradient descent optimization of description length[END_REF]: one looks for the most compact model (i.e. involving as few parameters as possible) that best explains the population or on ii) the definition of the mean shape as the minimiser of a registration-like criterion [START_REF] Durrleman | Statistical models of sets of curves and surfaces based on currents[END_REF][START_REF] Chui | Unsupervised learning of an atlas from unlabeled point-sets[END_REF][START_REF] Hufnagel | Generation of a statistical shape model with probabilistic point correspondences and em-icp[END_REF][START_REF] Wang | Simultaneous nonrigid registration of multiple point sets and atlas construction[END_REF][START_REF] Liu | Models of normal variation and local contrasts in hippocampal anatomy[END_REF]; the subsequent model simply explains the variability of the population with respect to the distance (typically a mixture of data superimposition distance and of smoothness of the underlying transformations) assumed by the registration criterion.

Among this last class of methods, the Mean Shape EM-ICP (MS-EM-ICP) recently developed by Hufnagel and colleagues [START_REF] Hufnagel | Generation of a statistical shape model with probabilistic point correspondences and em-icp[END_REF] is based on the original EM-ICP developed by Granger and Pennec [START_REF] Granger | Multi-scale EM-ICP: A fast and robust approach for surface registration[END_REF] and more recently extended to handle non-linear deformations efficiently [START_REF] Combès | An efficient EM-ICP algorithm for symmetric consistent non-linear registration of point sets[END_REF]. This algorithm relies on a probabilistic modelling of the point-to-point correspondences that allows i) a pragmatic definition of the superimposition of two point sets and ii) to work on a relatively smooth cost function to minimise [START_REF] Granger | Multi-scale EM-ICP: A fast and robust approach for surface registration[END_REF]. Moreover, this approach is well-grounded (monotonic convergence), generic (no assumption on tessellation/topology/number of points), can deal with large point sets and can be easily improved in many ways (adding priors, ...). However, only a few works have been dedicated to the investigation of the ability to this approach to perform groupwise studies.

The purpose of this paper is to extend the previous work of Hufnagel and colleagues and show that the EM-ICP approach leads to particularly simple and powerful solutions for groupwise studies. In particular, its underlying probabilistic modelling allows to define the mean shape of a given population and the projection of individual pointwise information on the mean shape in a simple way. In Section 2, we first introduce the mean shape EM-ICP (MS-EM-ICP) algorithm as an extension of the classical Procrustean mean to unlabelled point sets. From a population of point sets X = {X 1 , . . . , X C }, the MS-EM-ICP allows to estimate both a mean shape M and a "virtual" correspondent of each point of M in each point set X c . Then we identify some flaws of the MS-EM-ICP and develop two alternative extensions of the MS-EM-ICP allowing to tackle them. Finally, we show that these algorithms allow to perform statistical shape analysis (e.g. PCA) and to project individual local characteristics (e.g. cortical thickness, asymmetry, curvature) on the common reference M . In Section 3, we show how the previous tools can be used to perform a statistical analysis of asymmetries over and between populations of anatomical structures. These tools entirely rely on extensions of the EM-ICP algorithm. Finally, we illustrate our methods on real applications for the shape analysis of osseous labyrinths (Section 4.2) and the comparison global of cortical asymmetry between two patients with situs inversus and a control population (Section 4.3).

MEAN SHAPE AND VIRTUAL CORRESPONDENCES TO THE MEAN SHAPE

In this section, we first present the mean shape EM-ICP (MS-EM-ICP) algorithm (Section 2.1.2) as an extension of the Procrustean mean to unlabelled point sets (Section 2.1.1). Then, we study its limitations (Section 2.1.3) and we propose two extensions aiming at tackling these limitations (Section 2.2.1 and Section 2.2.2). Finally, we notice that the MS-EM-ICP algorithm outputs some quantities allowing to compute virtual correspondences from each X c to M and to project efficiently pointwise information of each X c to M (Section 2.3).

Mean shape EM-ICP

Basics: labelled point sets

The Procrustean mean is a Fréchet mean for the space of labelled point sets. This space is defined as the orbit space of the noncoincident point set configurations under the action of the rigidbody (or similarity) transformation. In other words, the shape of an object is defined as the geometrical information that remains when one filters out translation, rotation (and optionally scale).

The distance from a shape

V 1 = {v 1 1 , . . . , v 1 N } to another shape V 2 = {v 2 1 , . . . , v 2
N } is considered as the sum of the squared Euclidean distances between each labelled point after an optimal rigid (or rigid plus scale) alignment of the shapes. In case of rigid-body alignment, this Procrustean distance writes:

d 2 (V 1 , V 2 ) = min T X i ||T (v 1 i ) -v 2 i || 2 , ( 1
)
where T is a rigid-body transformation. This distance d allows the definition of the Fréchet mean for d of a set of labelled point sets

V = {V 1 , . . . , V C }: M = arg min M C X c=1 d 2 (M, V c ). (2)
M can be computed using the classical Procrustean superimposition framework:

Init: choose an arbitrary reference point set V ref in V
Step 1: compute the optimal rigid body transformation T c for all c, T c = arg min

T c P i ||T c (v c i ) -v ref i || 2
Step 2: compute the optimal mean shape M for all i, mi = 1

C P c T c (v c i )
Then the variability of the population around the mean can be studied by analysing the residuals (V c -M ) (e.g. using PCA).

Basics : unlabelled point sets

However, in practice, defining and labelling adequate points (landmarks) characterising the surfaces of interest is often problematic. To cope with this limitation, one can study unlabelled point sets representing the surfaces under study. The classical Procrustean distance between two labelled point sets is what remains when the optimal superimposition has been applied, thus in order to extend this approach to unlabelled point sets, we have to design a measure of what remains when the optimal superimposition and the optimal labelling have been applied. The EM-ICP algorithm, designed as a solution to a ML problem for a rigid-body transform T superposing a point set T (X 1 ) to its noised version X 2 , is a pragmatic and efficient solution for this purpose. This ML estimation is performed by an EM algorithm that amounts to the alternated minimisation of the following criterion with respect to two unknowns T /A:

δ 2 (X 1 , X 2 ) = min A,T X x 1 i ∈X 1 X x 2 j ∈X 2 Ai,j||T (x 1 i ) -x 2 j || 2 (3) +2σ 2 X i,j
Ai,j log(Ai,j), with ∀i, j, Ai,j ≥ 0 and ∀j, X i Ai,j = 1 where:

• A = (Ai,j) is the match matrix encoding the a posteriori probabilities of correspondence between points of X 1 and of X 2 . This probabilistic interpretation of A is linked to the presence of the barrier function P i,j Ai,j log(Ai,j) [START_REF] Chui | Unsupervised learning of an atlas from unlabeled point-sets[END_REF]. In essence, the greater Ai,j, the more likely the point x 1 i ∈ X 1 to be the correspondent of the point x 2 j ∈ X 2 . σ 2 is the Gaussian noise variance on X. This fuzzy control on the correspondences allows to handle efficiently problems due to differences of sampling/number of points: we do not look for one-to-one correspondences between points of each shape but instead for "fuzzy" correspondences linking each point of X 2 to the surface represented by T (X 1 ).

• T is a rigid-body (or rigid+scale) transformation superposing

X 1 on X 2 .
Strictly speaking, this departure δ is not a distance as it satisfies neither the symmetry constraint nor the triangle inequality. Note that the symmetry constraint could be achieved by symmetrising the criterion [START_REF] Combès | An efficient EM-ICP algorithm for symmetric consistent non-linear registration of point sets[END_REF] but for a sake of clarity, we do not discuss this aspect here.

Then, if we consider a set of C unlabelled point sets (of different size) defining a population under study X = {X 1 , . . . , X C } and M = {m1, ..., mN } as the mean shape corresponding to the population X (we consider the number of points N fixed), M can be written as:

M = arg min M C X c=1 δ 2 (M, X c ) (4)
M = arg min M C X c=1 min A c ,T c h X x c i ∈X c X m j ∈M A c i,j ||T c (x c i ) -mj|| 2 (5) +2σ 2 X x c i ∈X c ,m j ∈M A c i,j log(A c i,j ) i ,
with ∀c, ∀j, ∀iA c i,j ≥ 0, and ∀c, ∀j,

X i A c i,j = 1,
where:

• for all c, A c = (A c i,j ) is the match matrix linking X c and M , • T c is either a rigid-body or a similarity transformation superposing X c on M .

This formulation can be seen as an extension of the well-known EM-ICP algorithm [START_REF] Granger | Multi-scale EM-ICP: A fast and robust approach for surface registration[END_REF] and can be viewed as a ML problem modelling M as a noised version of T c (X c ) (∀c) that can be solved using an EM algorithm. This EM algorithm writes:

Algo MS1: MS-EM-ICP Init: M is one of the point sets from X. ∀c, T c = Id Step 1: update Ãc for all c, for all i, j, Ãc i,j = exp(-|| T c (x c i )-mj || 2 /(2σ 2 )) P k exp(-|| T c (x c k )-mj || 2 /(2σ 2 ))
Step 2: update T c for all c, T c = arg min

T c P i P j Ãc i,j ||T c (x c i ) -mj|| 2
Step 3: update M for all j, mj = 1

C P c P i Ãc i,j T c (x c i )
As for the original EM-ICP algorithm, this algorithm is iterative and decreases the value of the criterion (and increases the corresponding likelihood function) until local convergence.

M is typically initialised using a point set from X that appears both of good quality (no hole, no segmentation artifacts) and a "typical" representative of X. Step 1 can be performed efficiently using a kd-tree. When T is a rigid/similarity transformation Step 2 has a closed-form solution. Step 3 is computationally tractable.

In this framework, one can interpret the mean shape as the one for which the sum of "fuzzy residuals" Ãc i,j || T c (x c i ) -mj|| 2 subject to optimal transformations T c and optimal match matrices Ãc is minimal. In this sense, it extends the classical Procrustean mean.

Particularly, if we note ṽc j =

P i Ãi,j T c (x c i ) then mj = 1 C P c ṽc j . We call ṽc j the "virtual correspondent" of mj in T c (X c ); the term "virtual" is coined because v c j is not guaranteed to be in T c (X c
). The algorithm MS1 can then be simply rewritten as:

Step 1: update correspondences (ṽ c j ) for all c, compute T c and Ãc for all c, for all j, ṽc j =

P i Ãi,j T c (x c i ) Step 2: update M : for each j, mj = 1 C P c ṽc j
For a better understanding of this process, one can regard it as:

M = arg min M,V X c d 2 (M, V c (X c , M )), (6)
where M and V c (., .) are labelled point sets with N points each and d(, .) is the Procrustean distance described in Eq. 1. Each labelled point set

V c (X c , M ) = (ṽ c j ) (c ∈ [1, . . . , C]
) is obtained from the rigid (or rigid+scale) registration of X c on M . Then the previous algorithm can be described as i) the estimation of V c knowing the population X c and M and ii) the estimation of M knowing V c .

Limitations

This framework appears seductively simple: it naturally extends the classical Procrustes mean shape while solving the correspondence problem between the point sets under study, using all the points of the surfaces, in a unified probabilistic framework, where a welldefined likelihood criterion can be maximised using a 3-step iterative algorithm.

However, it is deceptively attractive and inherently flawed. In this framework, the point-to-point correspondences are computed based on the implicit hypothesis that the point sets are identical up to a rigid-body (or rigid+scale) transformation and a Gaussian noise. As a consequence, the further the shape of the surfaces are different from each other, the worse the correspondences. This is illustrated on Fig. 1. In this case, the mean shape will then be computed based on erroneous matches, and the subsequent analysis of variability around this mean shape will be strongly biased. In other words, this algorithm turns out to be useless when there is a high variability within the set of point sets of interest.

To cope with this problem, following Hufnagel et al. [START_REF] Hufnagel | Generation of a statistical shape model with probabilistic point correspondences and em-icp[END_REF], we could replace the rigid-body/similarity transformations by more general affine transformations. This improves the estimation of correspondences, but with the deleterious effect of removing the potentially relevant differences between the surfaces under study (such as non-uniform scaling or stretching).

In the following, we propose to explore this idea further by improving the previous definition of the mean shape by considering the transformations T c (c ∈ [1, . . . , C]) as non-linear transformations subject to regularisations L(T c ). Interestingly, depending on the way we interpret the original MS-EM-ICP criterion, this improvement can lead to two different paradigms. These two paradigms are equivalent when T c is a simple linear transformation but quite different when dealing with non-linear deformations. In the two following paragraphs, we discuss these two paradigms.

Handling non-linear transformations

First solution: adding penalised non-linear deformations in δ

The natural solution to handle non-linear transformations consists in replacing the linear transformation in Eq 3 by a non-linear transformation T subject to regularisation L:

δ 2 (X 1 , X 2 ) = min A c ,T c X x i ∈X 1 X x j ∈X 2 Ai,j||T (x 1 i ) -x 2 j || 2 (7)
+ 2σ 2 X i,j
Ai,j log(Ai,j) + αL(T )

where α weighs the strength of the regularisation over the data attachment term (including the barrier function P i,j Ai,j log(Ai,j)). Roughly speaking, the subsequent mean shape M is then defined as the shape that minimises both the overall "mismatches" P

x c i ∈X c P m j ∈M A c i,j ||T c (x c i ) -mj|| 2
and the deformation energies αL(T c ) between each X c and M . The choice of the parameter α and of L is crucial as it has a strong impact on the definition of the mean shape and thus on the interpretation of the subsequent statistical model. Moreover, in order not to capture simple differences of pose, orientation (and optionally scale) of the point set X c in the model, the regulariser L must not penalise rigid-body (or similarity) transformations.

Similarly to the original MS-EM-ICP, this modified criterion can be minimised by the following algorithm:

Algo MS2: non linear MS-EM-ICP

Step 1: update Ãc for all c, for all i, j,

Ãc i,j = exp(-|| T c (x c i )-mj || 2 /(2σ 2 )) P k exp(-|| T c (x c k )-mj || 2 /(2σ 2 ))
Step 2: update T c : for all c, T c = arg min

T c P i,j Ãc i,j ||T c (x c i ) -mj|| 2 + αL(T c )
Step 3: update M : for all j, mj = 1

C P c P i Ãc i,j T c (x c i )
Step 1 can be performed efficiently using a kd-tree. Step 2 can be solved efficiently using the techniques we recently proposed in the registration context [START_REF] Combès | An efficient EM-ICP algorithm for symmetric consistent non-linear registration of point sets[END_REF]. This algorithm is quite close to the one proposed by Chui and colleagues [START_REF] Chui | Unsupervised learning of an atlas from unlabeled point-sets[END_REF]. It mainly differs by the fact that Chui and colleagues introduce a clustering task in Step 1 allowing to reduce the size of matrices A c and thus the computational burden during Step 2. We think that our solution is more natural and simpler than the one they propose.

Second solution: improving correspondence estimation using non linear deformations

Although the previous solution provides a first proposition, it importantly modifies the nature of the distance and thus the interpretation of the estimated mean with respect to the initial Procrustean point of view. Indeed, by modifying the deformation model as we did in Section 2.2.1, one both affects i) the estimation of correspondences and ii) the nature of the distance d considered as the sum of X 1 (in dotted line)

T 1 (X 1 ) and M T 1 (X 1 ) and M mj mj mj i A 1 i,j T 1 (xi) = i A 1 i,j T 1 inv (xi) T 1 inv (X 1 ) and M T 1 (X 1 ) and M mj v 1 j = i A 1 i,j T 1 (xi) v 1 j i A 1 i,j T 1 (xi)
and M (in solid line) with T 1 rigid with T 1 non linear Fig. 1. First row: Illustration of the classical MS-EM-ICP (MS1). From left to right: i) the point set X 1 and M , ii) T (X 1 ) and M (T 1 is a rigid-body transformation) and the resulting matching

P i A 1 i,j T 1 (xi) and iii) the resulting correspon- dences v 1 j = P i A 1 i,j T 1 (xi),
Second row: Illustration of the improved MS-EM-ICP (MS3). From left to right: i) T 1 (X 1 ) on M (T 1 is a non-linear smooth transformation) and the resulting matching

P i A 1 i,j T 1 (xi) and ii) the resulting correspondences v 1 j = P i A 1 i,j T 1 inv (xi).
"mismatches" P

x i ∈X c P m j ∈M A c i,j ||T c (x c i )
-mj|| 2 and of the deformation energies αL(T c ) between each X c and M .

An alternative solution would consist in improving the estimation of the correspondences v c j by introducing a non-linear deformation linking M and X c inside the criterion. However, to preserve the spirit of the original Procrustean mean, the nonlinear deformations must affect only the estimation of the matrices A c not to change the nature of the shape space (that is defined as an orbit space under the action of a simple rigid or similarity transformation). This view results in the following algorithm:

Algo MS3: Improved Procrustean mean on unlabelled point sets

Step 1: update correspondences (ṽ c j ): ∀c, compute T c and Ãc ∀c, compute T c inv as the rigid component of T c ∀c, ∀j, compute ṽc j =

P i Ãc i,j T c inv (x c i ) Step 2: update M : for each j, mj = 1 C P c ṽc j
Step 1 can be performed efficiently using the techniques proposed in [START_REF] Combès | An efficient EM-ICP algorithm for symmetric consistent non-linear registration of point sets[END_REF]. Step 2 is computationally tractable. This algorithm can be seen as a generalisation of MS1 but we did not find any proper criterion it minimises. However, we observe that the value P c d(M, V c (X, M )) 2 decreases throughout the iterations and that the obtained mean shape fits well the data X. An illustration of this strategy is given in the second row of Figure 1. By contrast with MS2, this algorithm provides a mean that does not depend directly on the regulariser L but on the quality of the correspondences obtained during the registration process. As a result, we think that this approach is less dependent on the design of L. On the other hand, by explicitly taking into account the regularisation into δ, MS2 is more suited to characterise large deformations linking a point set of the population X c to the mean.

Conclusion

When modelling T c as linear transformations (without penalisation), the two developed algorithms MS2 and MS3 are equivalent. When using non-linear transformations, the two algorithms describe two different points of view. As of now, we did not investigate the respective merits of both frameworks for the estimation of statistical model and in the following, we choose to use the improved automatic Procrustean mean (MS3) to perform mean shape estimation and statistical shape analysis because of i) its simplicity of interpretation as opposed to MS2 (largely depending on the choice of L/α) and ii) the easiness of the underlying statistical analysis (that can be performed from the virtual correspondence similarly to what is done in the classical Procrustean framework). Finally, notice that the estimations given by the algorithms MS1, MS2 and MS3 are initialisation dependent and thus provide biased means (contrary to the classical Procrustean mean for labelled point sets). The sensitivity to initialisation is a key property to investigate in future works.

Virtual correspondence & projection

As previously mentioned, the match matrices allow to compute the "virtual" correspondent of each point of the mean shape M to each point of the X c s using the a posteriori probabilities of matching encoded in A c i,j (∀c):

v c j = P i A c i,j T c inv x c i .
In the following, we assume that the point sets X c are initially rigidly registered together and thus that Tinv = Id and v c j =

P i A c i,j x c i .
Another interesting use of this a posteriori probability concerns the projection of pointwise individual mapping on the mean shape M . Consider that each point x c i of X c contains a scalar information (e.g. cortical thickness, local asymmetry, curvature, ...) that we call s X c i . Then, the projection of S X c = (s X c i ) on M is given by s X c →M j = P i A c i,j s X c i . As ∀j, P i A c i,j = 1, the "interpolated" value at mj is a weighted mean of the values at different points of X c . This last property will be useful below for the projection of individual asymmetry maps on a reference template.

m j s X c →M j = i A c i,j s X c i x c i 1 x c i 2 x c i ... M X c
T c (X c ) (and subsequent A c ) M Fig. 2. Projection of scalar information on the mean shape: From left to right i) mean shape M and X c . ii) M and T c (X c ) (computed jointly with M ). iii) the subsequent correspondences (hot (resp. cold) colors indicate high (resp. low) values for A c i,j ) between each point mj and points of X c allow to compute the projected scalar information at point mj of M .

QUANTIFICATION OF ASYMMETRIES OVER/BETWEEN POPULATIONS

In this section, based on the previously described tools, we develop a pipeline for the quantification of asymmetries over or between populations. We consider that all the structures can be and are oriented in a common orthogonal frame consisting of three axes that we call anterior-posterior, left-right and head-foot.

Computation of individual asymmetry maps

Let X c ∈ X be a point set representing an anatomical structure under study. Its individual asymmetry maps are computed as follows:

1. computation of the approximate symmetry plane P of X c using a modified EM-ICP algorithm [START_REF] Combès | Automatic symmetry plane estimation of bilateral objects in point clouds[END_REF].

2. computation of the asymmetry field as the deformation field best superposing X c and SP (X c ) using another modified EM-ICP algorithm [START_REF] Combès | An efficient EM-ICP algorithm for symmetric consistent non-linear registration of point sets[END_REF].

3. computation of the 3 individual asymmetry maps by projecting each pointwise vector of the asymmetry field on the 3 coordinate axes. This allows to differentiate the anterior-posterior, left-right and head-foot components of the asymmetry field. This leads to 3 different scalar asymmetry maps for each point set X c . We call these asymmetry maps S X c ,AP , S X c ,LR and S X c ,HF .

Computation of a mean shape & projection of the asymmetry mappings

We consider C point sets X 1 , . . . , X C representing the C structures under study. Their individual asymmetry maps S X c ,AP S X c ,LR and S X c ,HF have to be put in a common geometry to be compared. For this purpose, we use the algorithm MS3 described above to iteratively compute:

1. the point set M representing the mean shape, 2. the match matrices {A 1 , . . . , A C } (describing the fuzzy point-to-point correspondences between the point sets X c and M ),

3. the transformations T c (best superposing the point sets X c and M ).

Once the mean point set M and the fuzzy match matrices A C are computed, we project each individual (scalar) asymmetry maps S X c ,AP , S X c ,LR and S X c ,HF on M which provides the normalised individual asymmetry maps S X i →AP , S X i →M,LR and S X i →M,HF . As pointed out in Section 2.3, the mapping at point mj ∈ M is defined as s X c →M,AP

j = P i A c ij s X c ,AP i
(and similarly for s X c →M,LR and s X c →M,HF ).

EXAMPLES OF APPLICATION

Implementation details

The transformation T is parametrised as a deformation field (i.e. T (x) = x + t(x)). The regularisation L(T ) = L(t) is chosen as a scalar Fourier based regulariser:

L(t = (t1, t2, t3) T ) = L(t1) + L(t2) + L(f3), with L(ti) = 1 (2π) 3 Z ∞ -∞ |t * i (ω)| 2 φ * (|ω|/b) dω,
where * is the Fourier transform operator, φ : IR → IR is an integrable function and b is a real positive rescaling factor. We choose φ as a Wu kernel. The resulting algorithmic solutions are given in a previous work [START_REF] Combès | An efficient EM-ICP algorithm for symmetric consistent non-linear registration of point sets[END_REF].

The parameters are initialised as: α = 400, σ 2 = 20 × S, b = 250 × S where S is the size (in metres) of the object; α, σ 2 , b are then decreased throughout the algorithm until they reach the respective values of 20, 5 × S and 150 × S.

Finally, notice that we add a second match matrix B that is column stochastic in addition to the first match matrix A that is row stochastic in the function δ (Eq. 3 and 7). This modification does not change the general minimisation strategy and tends to make the formulation of δ more symmetrical (particularly the matching process) and to improve the estimation of correspondences between points. We omitted these details in the previous sections for a sake of clarity.

Shape analysis on osseous labyrinths

The shape analysis of the inner ear, and especially of the vestibular system (composed of the three semicircular canals) could be relevant to understand diseases such as adolescent idiopathic scoliosis, which affects about 3% of children between 10 and 16 worldwide [START_REF] Zeng | Shape Analysis of Vestibular Systems in Adolescent Idiopathic Scoliosis Using Geodesic Spectra[END_REF]. More generally, the morphometric analysis of the bony/osseous labyrinth (that houses the membranous semicircular canals, the vestibule and the cochlea) has been shown to be key to the establishment of phylogenetic affinities between fossil hominids [START_REF] Spoor | The bony labyrinth of Neanderthals[END_REF].

We segmented 15 osseous labyrinths of modern Homo sapiens from CT images with amira.com, giving surfaces of about 30,000 points. We validated our non-linear registration method on osseous labyrinths in a previous paper [START_REF] Combès | An efficient EM-ICP algorithm for symmetric consistent non-linear registration of point sets[END_REF]. In Fig 3, we display the mean and two principal modes of variations around the mean obtained by analysing these 15 osseous labyrinths using MS3 and a principal component analysis on the residuals (M -V c ) (where V c are the virtual correspondents of M with respect to point set X c ). These results illustrate the ability of our algorithm to work on large and very convoluted structures with complicated topologies.

Global cerebral asymmetries of subjects with situs inversus

Subjects with situs inversus (SI) have the potential to provide unique clues into the developmental mechanisms underlying the brain anatomical asymmetry and its relationship with hemispheric dominance for language. Situs inversus is a very rare condition where all the visceral organs are on the opposite side of the body to where they would be expected and as if they were reflected in a mirror. Indeed, Kennedy [START_REF] Kennedy | Structural and functional brain asymmetries in human situs inversus totalis[END_REF] et al. report that the brain torque (reversed left-right differences in position and width of the frontal and occipital lobes) is reversed in subjects with situs inversus and that leftward anatomical asymmetry of language structures e.g. planum temporale is present in situs inversus. Moreover, Kennedy et al. suggest a left hemisphere functional dominance for language using fMRI [START_REF] Kennedy | Structural and functional brain asymmetries in human situs inversus totalis[END_REF]. Therefore, structurally there seems to be preserved local asymmetries of language-related structures in situs inversus but reversed torque and reversed functional dominance for language.

Material

Subjects: Data were available for two right-handed male subjects with situs inversus [START_REF] Kennedy | Structural and functional brain asymmetries in human situs inversus totalis[END_REF]. Both were aged 33 years. The control subjects are eleven right-handed male volunteers, ranging in age from 24 to 42 years (mean 28 ± 4.8). Preprocessing: For each subject's MRI, the following pipeline was applied. We first segmented grey matter and separated each hemisphere (surfer.nmr.mgh.harvard.edu). The mesh was smoothed to remove inter-subject variability arising from the gyri and sulci using a technique similar to the one proposed in [START_REF] Fournier | Surface-based method to evaluate global brain asymmetries:comparison between human and champanzee brains[END_REF]. By removing the gyri and the sulci, we concentrate our study on gross brain asymmetries.

Each generated mesh contains about 100K points.

Results

Pointwise results: We computed the collection of C (× 3) scalar normalised individual asymmetry maps S X 1 →M , . . . , S X C →M projected on a common mean mesh M . For each point of M and each of the 3 components of the asymmetry field, we computed the mean asymmetry and its standard deviation over the 11 subjects and we performed a pointwise t-test with the null hypothesis H0: "There is a perfect symmetry". We corrected the obtained p-values for multiple comparisons using permutation tests. Then, we compared the control and the SI populations by performing a pointwise permutation test (H0: both SI subjects belong to the control population) using the mean difference as a statistic (correcting p-values for multiple comparisons). In Figures 4 and6, we display components of the asymmetry field and views that appear the most significant. In Figure 5, we display the asymmetry maps for one of the two SI patients. Integrated results: For an easier interpretation of the results, we integrated the pointwise asymmetry values over the frontal and occipital lobes of each of the subjects. Both lobes were extracted following affine registration of the LONI LPBS40/SPM5 atlas [START_REF] Shattuck | Construction of a 3d probabilistic atlas of human cortical structures[END_REF] to each subject of each of the subjects. Then for each subject, we computed the average left-right ("width") and the posterior-anterior ("protrusion") components for both the frontal (F) and occipital (O) lobes. To be even more synthetic, for each subject we noted the lobe (L or R, i.e. left or right) that appears wider and that protrudes more than its counterpart. When the difference in width or protrusion between the hemispheres was low (< 10 -1 mm), we noted 0 instead of L or R. For both the width (L-R), protrusion (P-A) components of the asymmetry field over the frontal and occipital lobes), we performed a t-test with the null hypothesis H0: "There is a perfect symmetry". This integrated analysis was performed both on controls and SI subjects (Table 1 andTable 2).

Interpretation

Control population: At the population level, and at the chosen significance level (p=0.05, corrected for multiple comparisons) the left occipital lobe appears to be wider (or bending towards the other side) and protrudes more posteriously than the right occipital lobe (Fig. 4 and Table 1). This observation is in line with the literature. At last and interestingly, we also find an area in the parietal lobe with a small significant "vertical" (head-foot) asymmetry: the right side seems to be "higher" compared to its counterpart (Fig. 4). Table 1. Protrusion/width interhemispheric differences. We display these differences for the 11 subjects for occipital (O) and frontal (F) lobe. More details in the core of the text.

- SI population: The two SI subjects appear significantly different from the control subjects in the frontal lobes for the left-right component and in the occipital lobes for the antero-posterior component (Fig 6). Right occipital lobes protrude more anteriorly than the left in the two subjects, while it is the opposite in the frontal lobes. The left frontal lobes are also wider than the right counterpart and bend towards it. On the contrary, width measurements are different between the two subjects in the occipital lobes (Table 2). Interestingly, both SI and controls show significant, but opposite, protrusion asymmetries in the occipital lobes.

CONCLUSION

Our contributions are threefold. First, we analysed and extended previous works allowing to build a mean shape by using an EM-ICP variant. Second, we showed how the a posteriori probabilities of correspondences (encoded in the match matrices) estimated by the MS-EM-ICP algorithm allow to perform the statistical analysis of bilateral asymmetries. Third, we showed on some real applications that our improved automatic Procrustean mean (and subsequent statistical shape analysis) and our pipeline for asymmetry quantification allow to deal efficiently with large structures. Both methodologies (shape and asymmetry analysis) can deal with numerous applications. In particular, a more practical study comparing the global asymmetry of human and chimpanzee cortex is proposed in another ISBI 2011 paper [START_REF] Fournier | Surface-based method to evaluate global brain asymmetries:comparison between human and champanzee brains[END_REF].

An important point of our methodology is that there is no need to define a parametrisation on the surfaces to register, and no assumption on their topology. Moreover, there is no need to have one-toone correspondences between the surfaces to register, which allows to tackle the difference of sampling and the outliers in an efficient way. As previously mentioned, the process proposed to study analysis can easily be modified to analyse other local patterns such as cortical thickness, curvature, etc. Further works will mainly consist in investigating the respective merits of MS2 and MS3 for the estimation of statistical models and comparing them with existing methods [START_REF] Durrleman | Statistical models of sets of curves and surfaces based on currents[END_REF][START_REF] Heimann | 3d active shape models using gradient descent optimization of description length[END_REF][START_REF] Wang | Simultaneous nonrigid registration of multiple point sets and atlas construction[END_REF][START_REF] Liu | Models of normal variation and local contrasts in hippocampal anatomy[END_REF].

Fig. 3 .

 3 Fig. 3. Mean shape and first modes of variation on 15 osseous labyrinths. The magnitude of the deformations around the mean are mapped. The first mode of variation can be mainly interpreted as a change of size and of position of the lateral semicircular canal. The second mode of variation can be mainly interpreted as a change of size of the posterior semicircular canal.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Mean (first row) and p-value (second row) for each of the 3 components of the asymmetry field (11 subjects). From left to right: left-right, head-foot, posterior-anterior components of the asymmetry field. Only the most significant view for each of the 3 tests is displayed (From left to right; posterior, above and posterior)

Fig. 6 .

 6 Fig. 6. p-value of permutation test with H0: SI population and control population are identical From left to right: left-right and posterior-anterior components of the asymmetry field.

Table 2 .

 2 Integrated values for the two SI subjects and p-value (t-test) with H0: SI subjects belong to the control population.

		Subject 1 Subject 2	p
	O protrusion	R	R	p<0.05
	O width	L	R	p>0.05
	F protrusion	L	L	p<0.05
	F width	L	L	p<0.05


Arnaud Debussche
email: arnaud.debussche@bretagne.ens-cachan.fr

Erwan Faou
email: erwan.faou@inria.fr

Weak backward error analysis for SDEs

Keywords: backward error analysis, stochastic differential equations, exponential mixing, numerical scheme, Kolmogorov equation, weak error. MSC number: 65C30, 60H35, 37M25

We consider numerical approximations of stochastic differential equations by the Euler method. In the case where the SDE is elliptic or hypoelliptic, we show a weak backward error analysis result in the sense that the generator associated with the numerical solution coincides with the solution of a modified Kolmogorov equation up to high order terms with respect to the stepsize. This implies that every invariant measure of the numerical scheme is close to a modified invariant measure obtained by asymptotic expansion. Moreover, we prove that, up to negligible terms, the dynamic associated with the Euler scheme is exponentially mixing.

Résumé

Nous étudions la discrétisation d'une équation différentielle stochastique (EDS) par le schéma d'Euler. Dans le cas d'une EDS elliptique ou hypoelliptique nous montrons un résultat d'analyse d'erreur rétrograde : une fonctionnelle de la solution numérique est proche de la solution d'une équation de Kolmogorov modifiée à des ordres arbitrairement élevés par rapport au pas de discrètisation. On obtient ainsi que toute mesure invariante du schéma numérique est proche d'une mesure invariante modifiée obtenue par développement asymptotique. De plus, le schéma est exponentiellement mélangeant à des ordres arbitrairement élevés.

Introduction

In the last decades, backward error analysis has become one of the most powerful tool to analyze the long time behavior of numerical schemes applied to evolution equations. The main idea can be described as follows: Let us consider an ordinary differential equation of the form ẏ(t) = f (y(t)), where f : R n → R n is a smooth vector field, and denote by ϕ f t (y) the associated flow. By definition, a numerical method defines for a small time step τ an approximation Φ τ of the exact flow ϕ τ : We have for bounded y ∈ R n , Φ τ (y) = ϕ f τ (y) + O(τ r+1 ) where r is the order of the method. The idea of backward error analysis is to show that Φ τ can be interpreted as the exact flow ϕ fτ τ or a modified vector field defined as a series in powers of τ

f τ = f + τ r f r + τ r+1 f r+1 + • • • ,
where f ℓ , ℓ ≥ r are vector fields depending on the numerical method. In general, the series defining f τ does not converge, but it can be shown that for bounded y, we have for arbitrary N

Φ τ (y) = ϕ f N τ τ (y) + C N τ N ,
where f N τ is the truncated series:

f N τ = f + τ r f r + • • • + τ N f N .
Under some analyticity assumptions, the constant C N τ N can be optimized in N , so that the error term in the previous equation can be made exponentially small with respect to τ . Such a result is very important and has many applications in the case where f has some strong geometric properties, such as Hamiltonian or reversible structure. In this situation, and under some compatibility conditions on the numerical method Φ τ , the modified vector field f τ inherits the structure of f . For example if Φ τ is symplectic and f Hamiltonian, then f τ remains Hamiltonian. This has major consequences such as the preservation of a modified Hamiltonian over very long time (of order τ -N ) for the numerical solution, from which we can deduce long time stability results, existence of numerical invariant tori in the integrable case, etc... In the Hamiltonian case, this idea goes back to Moser [START_REF] Moser | Lectures on Hamiltonian systems[END_REF], but was applied later to symplectic integrator by Benettin & Giorgilli [START_REF] Benettin | On the Hamiltonian interpolation of near to the identity symplectic mappings with application to symplectic integration algorithms[END_REF], Hairer & Lubich [START_REF] Hairer | The life-span of backward error analysis for numerical integrators[END_REF] and Reich [START_REF] Reich | Backward error analysis for numerical integrators[END_REF]. Such results now form the core of the modern geometric numerical integration theory for which we refer to the classical textbooks [START_REF] Hairer | Geometric Numerical Integration[END_REF] and [START_REF] Leimkuhler | Simulating Hamiltonian dynamics[END_REF].

More recently, these ideas have been extended in some situations to Hamiltonian PDEs: First in the linear case [START_REF] Debussche | Modified energy for split-step methods applied to the linear Schrödinger equation[END_REF], and then in the semilinear case (nonlinear Schrödinger or wave equations), see [START_REF] Faou | Grébert Hamiltonian interpolation of splitting approximations for nonlinear PDEs[END_REF][START_REF] Faou | Geometric numerical integration of Hamiltonian PDEs and applications to computational quantum mechanics[END_REF].

As far as stochastic differential equations (SDEs) are concerned, this approach has not been developed very much so far. Let us recall that given a SDE in R d of the form dX = f (X)dt + σ(X)dW discretized by a numerical scheme -such as the Euler scheme for instancewith time step τ providing a discrete sequence (X p ) p∈ N , then the error can be measured in the strong or weak sense. Strong error means that X p is a pathwise approximation of X(pτ ), and it is well known that the Euler scheme has strong order 1/2. Under standard assumptions on f and σ, we have

E sup p=0,...,[T /τ ] X p -X(t p ) k ≤ c 1 (k, T )τ 1/2 , k ≥ 1, T > 0.
In this work, we consider another error which is often more important. We investigate the weak error which concerns the law of the solution. The Euler scheme has weak order 1. Under suitable smoothness assumptions on f , σ and ϕ : R d → R (see for instance [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF][START_REF] Milstein | Stochastic Numerics for Mathematical Physics[END_REF][START_REF] Talay | Probabilistic numerical methods for partial differential equations: elements of analysis[END_REF]):

|E (ϕ(X p )) -E (ϕ(X(t p )))| ≤ c 2 (ϕ, T )τ, p = 0, . . . , [T /τ ], T > 0.
An attempt has been made by Shardlow [START_REF] Shardlow | Modified equations for stochastic differential equations[END_REF] to extend the backward error analysis to this context. He has shown that the construction of a modified SDE associated with the Euler scheme can be performed, but only at the first step, ie for N = 2, and only for additive noise, ie when σ(X) does not depend on X. In this case, he is able to write down a modified SDE:

d X = f ( X)dt + σ( X)dW such that E (ϕ(X p )) -E ϕ( X(t p )) ≤ c 3 (ϕ, T )τ 2 , p = 0, . . . , [T /τ ], T > 0.
He explains that for multiplicative noise or higher order, there are too many conditions to be satisfied by the coefficients of the modified equations.

In this paper we take another approach, and build a modified equation not at the level of the SDE, but at the level of the generator associated with the process solution of the SDE. It is well known that given ϕ : R d → R and denoting by X(t, x) the solution of the SDE satisfying X(0) = x, the function u(t, x) = E(ϕ(X(t, x))) satisfies the Kolmogorov equation

∂ t u(t, x) = L(x, ∂ x )u(t, x),
where L is the order 2 Kolmogorov operator associated with the SDE.

In the case of the Euler method applied to a SDE, we show that with the numerical solution, we can associate modified Kolmogorov operator of the form

L(τ, x, ∂ x ) = L(x, ∂ x ) + τ L 1 (x, ∂ x ) + τ L 2 (x, ∂ x ) + • • •
where L ℓ , ℓ ≥ 1 are some modified operator of order 2ℓ + 2. Again, the series does not converge but truncated series:

L (N ) (τ, x, ∂ x ) = L(x, ∂ x ) + τ L 1 (x, ∂ x ) + • • • + τ L N (x, ∂ x ) are considered.
Note that in contrast with the classical case, we do not have a modified SDE and cannot straightforwardly define a solution to the modified equation

∂ t v N (t, x) = L (N ) (τ, x, ∂ x )v N (t, x).
However, in the case where the SDE is elliptic or hypoelliptic, we can build an approximated solution v (N ) such that

E (ϕ(X p )) -v (N ) (pτ, x) ≤ c 4 (ϕ, T, N )τ N , p = 0, . . . , [T /τ ], T > 0.
Furthermore, using the exponential convergence to equilibrium, we prove that in fact the constant c 4 does not depend on T so that we have an approximation result valid on very long times. We also show that there exist a modified invariant measure for L (N ) (τ, x, ∂ x ).

We can then use this weak backward error analysis to prove that the numerical solution X p , p ≥ 1 obtained by the Euler scheme is exponentially mixing up to some very small error, and for all times. This is typically a geometric numerical integration result in the sense that we prove the persistence of a qualitative property of the exact flow (exponential mixing) to the numerical approximation, over long times.

Note that v (N ) is in fact constructed as a truncated series: v (N ) = N n=0 τ n v n and that v 0 = u is the solution of the Kolmogorov equation. Therefore, our result provides an expansion of the error as in [START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic differential equations[END_REF] (see also [START_REF] Bally | The law of the Euler scheme for stochastic differential equations (I): convergence rate of the distribution function[END_REF][START_REF] Bally | The law of the Euler scheme for stochastic differential equations (II): convergence rate of the density[END_REF]). However, the expansion is different here.

Error estimates on long times for elliptic and hypoelliptic SDEs have already been proved. In [START_REF] Mattingly | Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise[END_REF][START_REF] Talay | Second order discretization schemes of stochastic differential systems for the computation of the invariant law[END_REF][START_REF] Talay | Probabilistic numerical methods for partial differential equations: elements of analysis[END_REF][START_REF] Talay | Stochastic Hamiltonian dissipative systems: exponential convergence to the invariant measure, and discretization by the implicit Euler scheme[END_REF], it is shown that for a sufficiently small time step the Euler scheme defines an ergodic process and that the invariant measure of the Euler scheme is close to the invariant measure of the SDE. In [START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic differential equations[END_REF], the first term of an expansion of the invariant measure of the Euler scheme with respect to τ is also given. In our work, we provide the expansion at any order.

We emphasize that in our result, there is no particular smallness assumption on the stepsize τ used to define the numerical solution. In particular, the discrete process is not supposed to have a unique invariant measure, as in [START_REF] Mattingly | Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise[END_REF] or [START_REF] Talay | Second order discretization schemes of stochastic differential systems for the computation of the invariant law[END_REF][START_REF] Talay | Probabilistic numerical methods for partial differential equations: elements of analysis[END_REF][START_REF] Talay | Stochastic Hamiltonian dissipative systems: exponential convergence to the invariant measure, and discretization by the implicit Euler scheme[END_REF]. This is also the case in the recent work [START_REF] Mattingly | Convergence of numerical time-averaging and stationary measures via Poisson equations[END_REF]. There, it is shown that given an elliptic or hypoelliptic SDE, the ergodic averages provided by the Euler scheme are asymptotically close to the average of the invariant measure of the SDE. Higher order schemes also considered. The main tool in [START_REF] Mattingly | Convergence of numerical time-averaging and stationary measures via Poisson equations[END_REF] is the ellipticity or hypoellipticity of the Poisson equation, ie the equation L(x, ∂ x )u = g.

As in [START_REF] Mattingly | Convergence of numerical time-averaging and stationary measures via Poisson equations[END_REF], we consider the case where the SDE is set on the torus T d . This simplifies the presentation and the main ideas are not hidden by technical difficulties. In the same spirit, we only study the Euler scheme. In a forthcoming article, we will present more realistic applications of our method for SDEs set on R n with polynomial growth coefficients under suitable assumptions, and for more general schemes. As an example, we will treat the Langevin equation as in [START_REF] Talay | Stochastic Hamiltonian dissipative systems: exponential convergence to the invariant measure, and discretization by the implicit Euler scheme[END_REF].

Preliminaries

We consider the stochastic differential equation

dX(t) = f (X(t))dt + σ(X(t))dW,
where the unknown X = (X i ) i=1...,d lives in the d-dimensional torus T d . Also, f = (f i ) d i=1 and σ = (σ i ℓ (x)) i=1,...,d,ℓ=1,...,m are smooth vector fields periodic in x ∈ T d . The process (W 1 (t), . . . , W m (t)) is a m-dimensional standard Wiener process over a probability space (Ω, F, P) endowed with a filtration (F t ) t≥0 . Using these notations, we rewrite the equation as:

dX i (t) = f i (X(t))dt + m ℓ=1 σ i ℓ (X(t))dW ℓ (t), i = 1, . . . , d.
In all the paper, smooth functions means C ∞ functions. Given a smooth function ψ defined on T d , we denote by

ψ C k its norm in C k (T d , R). We also denote by ψ ∞ = sup x∈T d |ψ(x)|. For a multiindex k = (k 1 , . . . , k d ) ∈ N d , we set |k| = k 1 + • • • + k d and ∂ k ψ(x) = ∂ |k| ψ(x) ∂x k 1 1 • • • ∂x k d d , x ∈ T d . Therefore ψ C k := sup j=(j 1 ,...,j d ) |j|≤k |∂ k ψ(x)|.
We also define the semi-norm

|ψ| C k := sup j=(j 1 ,...,j d ) 1≤|j|≤k |∂ k ψ(x)|.
In the following, we assume that f and σ are smooth and since we are working with a stochastic differential equation on the torus, standard theorems give existence and uniqueness of a solution for any initial data X(0) = x ∈ T d . We denote this solution by X(t, x), t ≥ 0. Also, since we chose to work on the torus, we do not have any problem of possible unbounded moments and this solution has clearly all moments finite.

We denote by L(x, ∂ x ) the Kolmogorov generator associated with the stochastic equation:

L(x, ∂ x )v(x) = f i (x)∂ i v(x) + a ij (x)∂ ij v(x),
where we use the summation convention for repeated indices and

∂ i = ∂ x i and a ij (x) := 1 2 m ℓ=1 σ i ℓ (x)σ j ℓ (x).
It is well known that the Kolmogorov equation:

du dt = L(x, ∂ x )u, x ∈ T d , t > 0, u(0, x) = ϕ(x), x ∈ T d , (2.1)
with periodic boundary conditions has a unique solution for a smooth function ϕ and that for all x ∈ T d :

u(t, x) = E(ϕ(X(t, x))).
Moreover, this solution is smooth. In the following, we write: u(t) = P t ϕ so that (P t ) t≥0 is the transition semigroup associated with the Markov process (X(t, x)) t≥0, x∈T d . Note that we use the standard identification u(t) = u(t, •).

We wish to investigate the approximation properties of the Euler scheme for long times. We need assumptions on the long time behavior of the law of the solutions of (2.1), ie of the law of the Markov process. We assume the following mixing properties:

[H1] There exists a C ∞ (T d , R) function ρ ≥ 0 such that L(x, ∂ x ) * ρ(x) = 0 and T d ρ(x)dx = 1. (2.2)
In other words, the measure ρ(x)dx is invariant by X(t, x).

[H2] Let g ∈ C ∞ (T d , R), and assume that T d g(x)dx = 0. Then there exists a unique function µ

(x) ∈ C ∞ (T d , R) such that L(x, ∂ x ) * µ(x) = g(x)
, and

T d µ(x)ρ(x)dx = 0. (2.
3)

[H3] Let u(t, x) be the solution of (2.1). Assume that T d ϕ(x)ρ(x)dx = 0, then there exists a constant λ and, for each

k ∈ N, a polynomial p k (t), such that if ϕ(x) ∈ C ∞ (T d , R) we have the estimates ∀ t ≥ 0, u(t, •) C k ≤ p k (t)e -λt ϕ C k . (2.4)
These hypothesis are usually satisfied under elliptic or hypoelliptic assumptions on the operator L(x, ∂ x ). The reader may find in [START_REF] Nualart | Malliavin Calculus and Related Topics[END_REF] conditions to ensure [H1]. We also refer to [START_REF] Bally | The law of the Euler scheme for stochastic differential equations (I): convergence rate of the distribution function[END_REF][START_REF] Mattingly | Convergence of numerical time-averaging and stationary measures via Poisson equations[END_REF] for a general definition of hypoelliticity and applications to numerical schemes. Combining kernel estimates for hypoelliptic diffusion ( [START_REF] Bally | The law of the Euler scheme for stochastic differential equations (I): convergence rate of the distribution function[END_REF][START_REF] Kusuoka | Applications of the Malliavin Calculus, part II[END_REF]) and exponential convergence to equilibrium, [H3] can be proved. Note that similar estimates are used in [START_REF] Mattingly | Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise[END_REF][START_REF] Talay | Second order discretization schemes of stochastic differential systems for the computation of the invariant law[END_REF][START_REF] Talay | Probabilistic numerical methods for partial differential equations: elements of analysis[END_REF][START_REF] Talay | Stochastic Hamiltonian dissipative systems: exponential convergence to the invariant measure, and discretization by the implicit Euler scheme[END_REF], where specific examples are considered. Finally, we mention that these hypothesis can be proved to be fulfilled using partial differential equations techniques (see [START_REF] Hörmander | The analysis of linear partial differential operators[END_REF]).

For a smooth function ψ, we set

ψ = T d ψ(x)ρ(x)dx.
Note that by [H3], we have for any solution of (2.1) and k ∈ N

∀ t ≥ 0, u(t, •) -ϕ C k ≤ p k (t)e -λt ϕ -ϕ C k . (2.5)
Now for a small time step τ > 0 and x ∈ T d , we consider the Euler method defined, for i = 1, . . . , d, by X 0 = x and the formula

X i n+1 = X i n + τ f i (X n ) + σ i ℓ (X n )(W ℓ ((n + 1)τ ) -W ℓ (nτ )), (2.6)
for n ≥ 0. Our main result can be stated as follows:

Theorem 2.1 Let N and τ 0 > 0 be fixed. Then there exists a modified smooth density

µ N (x) = ρ(x) + τ µ 1 (x) + • • • + τ N µ N (x)
such that T d µ N (x)dx = 1, a constant C N and a polynomial P N (t) such that the following holds: For all smooth function function ϕ(x) on T d , we have

∀ p ∈ N, Eϕ(X p ) - T d ϕdµ N ∞ ≤ P N (t p )e -λtp + C N τ N ϕ C 8N+2 , (2.7)
where for all p, t p = pτ and dµ N (x) = µ N (x)dx.

This result can be viewed as a discrete version of (2.5). Note that it implies that all the invariant measure of the numerical process X p are close to dµ N up to a very small error term C N τ N .

Using this result, we can also recover the weak convergence result

∀ p ≥ 1, Eϕ(X p ) - T d ϕdρ ∞ ≤ C(ϕ)τ
for some constant C(ϕ) depending on ϕ, and where we set dρ(x) = ρ(x)dx. This can be compared with [START_REF] Mattingly | Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise[END_REF][START_REF] Mattingly | Convergence of numerical time-averaging and stationary measures via Poisson equations[END_REF][START_REF] Talay | Second order discretization schemes of stochastic differential systems for the computation of the invariant law[END_REF][START_REF] Talay | Probabilistic numerical methods for partial differential equations: elements of analysis[END_REF]. As in [START_REF] Mattingly | Convergence of numerical time-averaging and stationary measures via Poisson equations[END_REF], the only assumption made on τ is that τ ≤ τ 0 where τ 0 is any fixed number. The influence of τ 0 is only reflected in the constants in the right-hand side -we can for example take τ 0 = 1. In particular, we do not assume that X p has a unique invariant measure -something that would be guaranteed only if τ is small enough. We also recover an expansion of the invariant measure as in [START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic differential equations[END_REF].

In the next sections, the constants appearing in the estimate depend in general on bounds on derivatives of f and g defining the SDE. They will also depend in general on τ 0 and N , but not on ϕ.

Asymptotic expansion of the weak error

We have the formal expansion for small t:

u(t, x) = ϕ(x) + tL(x, ∂ x )ϕ(x) + t 2 2 L(x, ∂ x ) 2 ϕ(x) + • • • + t n n! L(x, ∂ x ) n ϕ(x) + • • •
This is just obtained by Taylor expansion in time since by (2.1)

d n dt n u(t, x) = L(x, ∂ x ) n u(t, x) and in particular d n dt n u(0, x) = L(x, ∂ x ) n ϕ(x).
Since the solution u of the Kolmogorov equation is smooth and has its derivatives bounded in terms of the initial data ϕ, the above formal expansion can be justified and we have the following proposition whose proof is easy and left to the reader:

Proposition 3.1 Assume that ϕ ∈ C ∞ (T d , R
), and let τ 0 > 0. Then for all N , there exist a constant C N such that for all τ < τ 0 ,

u(τ, x) - N n=0 τ n n! L(x, ∂ x ) n ϕ(x) ∞ ≤ C N τ N +1 |ϕ| C 2N+2 . (3.1)
With the Euler scheme defined in (2.6), we associate the continuous process

Xi x (t) = X i n + (t -nτ )f i (X n ) + σ i ℓ (X n )(W ℓ (t) -W ℓ (nτ )), t ∈ [nτ, (n + 1)τ ], (3.2)
and X(nτ ) = X n . We thus have X n+1 = Xx (t n+1 ). The process (3.2) satisfies the equation d Xi

x (t) = f i (X n )dt + σ i ℓ (X n )dW ℓ (t), t ∈ [nτ, (n + 1)τ ]. (3.3)
Clearly, (X n ) defines a discrete in time homogeneous Markov process but Xx is not Markov.

In this work, we are only interested in the distributions of the solutions and of their approximation. We now examine in detail the first time step and its approximation properties in terms of the law. By Markov property, it is sufficient to then obtain information at all steps. Next result gives an expansion similar to Proposition 3.1 for the Euler process. Theorem 3.2 Then for all n ≥ 1, there exist operators A n (x, ∂ x ) of order 2n, such that for all N ≥ 1, there exist a constant

C N satisfying E ϕ( Xx (τ )) - N n=0 τ n A n (x, ∂ x )ϕ(x) ∞ ≤ C N τ N +1 |ϕ| C 2N+2 , (3.4
)

for all ϕ ∈ C ∞ (T d , R) and x ∈ T d .
Proof. Using (3.3) and the Itô formula, we get for t ≤ τ , dϕ( Xx (τ )) = L(x, ∂ x )ϕ( Xx (t)) + σ i ℓ (x)∂ i ϕ( Xx (t))dW ℓ (x), or equivalently,

ϕ( Xx (t)) = ϕ(x) + t 0 L(x, ∂ x )ϕ( Xx (s))ds + t 0 σ i ℓ (x)∂ i ϕ( Xx (s))dW ℓ (s). (3.5)
Note that the last term is a martingale. We define the operator

R 0,ℓ (x, ∂ x ) = σ i ℓ (x)∂ i .
We have for all s and all x,

L(x, ∂ x )ϕ( Xx (s)) = (f i (x)∂ i ϕ)( Xx (s)) + (a ij (x)∂ ij ϕ)( Xx (s)).
Hence applying (3.5) to ∂ i ϕ and ∂ ij ϕ, we obtain

L(x, ∂ x )ϕ( Xx (s)) = L(x, ∂ x )ϕ(x) + s 0 f i (x)f j (x)∂ ij ϕ( Xx (σ)) + f i (x)a nℓ (x)∂ nℓi ϕ( Xx (σ))dσ + s 0 a ij (x)f n (x)∂ ijn ϕ( Xx (σ)) + a ij (x)a nℓ (x)∂ ijnℓ ϕ( Xx (σ))dσ + s 0 f i (x)σ j ℓ (x)∂ ij ϕ( Xx (σ))dW ℓ (σ) + s 0 a in (x)σ j ℓ (x)∂ inj ϕ( Xx (σ))dW ℓ (σ).
We set A 0 = I, A 1 = L and plugg this in (3.5) to obtain

ϕ( Xx (t)) = ϕ(x) + tA 1 (x, ∂ x )ϕ(x) + t 0 s 0 A 2 (x, ∂ x )ϕ( Xx (σ))dσds + t 0 s 0 R 1,ℓ (x, ∂ x )ϕ( Xx (σ))dW ℓ (σ)ds + t 0 R 0,ℓ (x, ∂ x )ϕ( Xx (σ))dW ℓ (σ).
where

A 2 (x, ∂ x ) = f i (x)f j (x)∂ ij + f i (x)a nℓ (x)∂ nℓi + a ij (x)f n (x)∂ ijn + a ij (x)a nℓ (x)∂ ijnℓ
is an operator of order 4, and

R 1,ℓ (x, ∂ x ) = f i (x)σ j ℓ (x)∂ ij + a in (x)σ j ℓ (x)
∂ inj are operators of order 3. Taking the expectation so that the last two term disappear, we easily deduce the result for N = 1. Let us now prove recursively that there exist operators A n (x, ∂ x ) of order 2n and R n,ℓ (x, ∂ x ) of order 2n + 1, such that

ϕ( Xx (t)) = ϕ(x) + tL(x, ∂ x )ϕ(x) + N n=2 t n A n (x, ∂ x )ϕ(x) + t 0 • • • s N 0 A N +1 (x, ∂ x )ϕ( Xx (s N +1 ))ds 1 • • • ds N +1 + N n=1 t 0 • • • sn 0 R n,ℓ (x, ∂ x )ϕ( Xx (s n+1 ))ds 1 • • • ds n dW ℓ (s n+1 ).
(3.6) Note the expectation of the last term vanishes so that (3.6) easily implies (3.4). To prove (3.6), assume that A N +1 and R N,ℓ are known, and let us decompose A N +1 as A N +1 (x, ∂ x ) = A j N +1 (x)∂ j , where j = (j 1 , . . . , j m ) are multiindices (with the summation convention) and A j N +1 smooth functions of x. Such decompistion is easy to write for N = 1 or 2. We apply (3.5) to ∂ j ϕ( Xx (s N +1 )) for a given muti-index j, and obtain

A j N +1 (x)∂ j ϕ( Xx (s N +1 )) = A j N +1 (x)∂ j ϕ(x) + s N+1 0 A j N +1 (x)f n (x)∂ n ∂ j ϕ( Xx (s N +2 ))ds N +2 + s N+1 0 A j N +1 (x)a nℓ (x)∂ nℓ ∂ j ϕ( Xx (s N +2 ))ds N +2 + s N+1 0 A j N +1 (x)σ n ℓ (x)∂ n ∂ j ϕ( Xx (s N +2 ))dW ℓ (s N +2 ).
We thus choose

A N +2 (x) = A j N +1 (x)f n (x)∂ n ∂ j + A j N +1 (x)a ij (x)∂ ij ∂ j . and R N +1,ℓ (x) = A j N +1 (x)σ n ℓ (x)∂ n ∂ j ,

Modified generator 4.1 Formal series analysis

Let us now consider τ as fixed. We want to construct a formal series

L(τ ; x, ∂ x ) = L(x, ∂ x ) + τ L 1 (x, ∂ x ) + • • • τ n L n (x, ∂ x ) + • • • (4.1)
with operator coefficients L n (x, ∂ x ) smooth on T d , and such that formally the solution v(t, x) at time t = τ of the equation

∂ t v(t, x) = L(τ ; x, ∂ x )v(t, x), v(0, x) = ϕ(x)
coincides in the sense of asymptotic expansion with the approximation of the transition semigroup E(ϕ( Xx (τ ))) studied in the previous section. In other words, we want to have the equality in the sense of asymptotic expansion in powers of τ

exp(τ L(τ ; x, ∂ x ))ϕ(x) = ϕ(x) + n≥1 τ n A n (x, ∂ x )ϕ(x),
where the operators A n (x, ∂ x ) are defined in Theorem 3.2.

Formally, this equation can be written

exp(τ L(τ ; x, ∂ x )) -Id = τ A(τ ) (4.2)
where

A(τ ) = n≥1 τ n-1 A n . We have exp(τ L(τ ; x, ∂ x )) -Id = τ L(τ ; x, ∂ x ) n≥0 τ n (n + 1)! L(τ ; x, ∂ x ) n .
Note that the (formal) inverse of the series is given by

n≥0 τ n (n + 1)! L(τ ; x, ∂ x ) n -1 = n≥0 B n n! τ n L(τ ; x, ∂ x ) n .
where the B n are the Bernoulli numbers: see for instance [START_REF] Hairer | Geometric Numerical Integration[END_REF][START_REF] Faou | Geometric numerical integration of Hamiltonian PDEs and applications to computational quantum mechanics[END_REF] and [START_REF] Debussche | Modified energy for split-step methods applied to the linear Schrödinger equation[END_REF] for a similar analysis involving operators. Hence equations (4.1), (4.2) are equivalent in the sense of formal series to

L(τ ; x, ∂ x ) = ℓ≥0 B ℓ ℓ! τ ℓ L(τ ; x, ∂ x ) ℓ A(τ ) = n≥0 τ n   A n+1 + n ℓ=1 B ℓ ℓ! n 1 +•••+n ℓ +n ℓ+1 =n-ℓ L n 1 • • • L n ℓ A n ℓ+1 +1   .
(4.3) Identifying the right hand sides of (4.1) and (4.3), we get the following recursion formula

L n = A n+1 + n ℓ=1 B ℓ ℓ! n 1 +•••+n ℓ +n ℓ+1 =n-ℓ L n 1 • • • L n ℓ A n ℓ+1 +1 . (4.4)
Each of the terms of the above sum is an operator of order 2n + 2 with smooth coefficients and therefore L n is also an operator of order 2n + 2 with smooth coefficients.

Note that (4.2) gives immediately the inverse relation of this formal series equation:

A n = n ℓ=1 1 ℓ! n 1 +•+n ℓ =n-ℓ L n 1 • • • L n ℓ . (4.5)
Moreover, we have clearly

L n (x, ∂ x )1 = 0,
where 1 denote the constant function equal to 1.

Approximate solution of the modified flow

For a given N , we have constructed in the previous section an operator

L (N ) (τ ; x, ∂ x ) = L(x, ∂ x ) + N n=1 τ n L n (x, ∂ x ). (4.6)
In order to perform weak backward error analysis and estimate recursively the modified invariant law of the numerical process, we should be able to define a solution v N (t, x) of the modified flow

∂ t v N (t, x) = L (N ) (τ ; x, ∂ x )v N (t, x), v N (0, x) = ϕ(x). (4.7)
However, in our situation we do not know whether this equation has a solution. This is in contrast with standard backward error analysis where the modified flow can always be defined. The goal of the next proposition is to give a proper definition of the modified flow (4.7).

Theorem 4.1 Let ϕ be a smooth functions on T d . For all n ∈ N, there exist smooth functions v ℓ (t, x), defined for all times t ≥ 0, and such that for all t ≥ 0 and n ∈ N,

∂ t v n (t, x) -Lv n (t, x) = n ℓ=1 L ℓ v n-ℓ (t, x), (4.8)
with initial conditions v 0 (0, x) = ϕ(x) and v n (0, x) = 0 for n ≥ 1. For all N ≥ 0, setting

v (N ) (t, x) = N k=0 τ N v n (t, x), (4.9)
then the following holds: (i) There exists a constant C N such that for all time t ≥ 0, and all τ ≥ 0,

Ev (N ) (t, Xx (τ )) -v (N ) (t + τ, x) ∞ ≤ C N τ N +1 sup s∈(0,τ ) n=0,...,N |v n (t + s)| C 4N+2 . (4.10) (ii) For fixed τ 0 > 0, there exists a constant C N such that for all τ ≤ τ 0 , Eϕ( Xx (τ )) -v (N ) (τ, x) ∞ ≤ C N τ N +1 ϕ C 4N+2 .
Proof. For n = 0, the equation (4.8) implies v 0 (t, x) = u(t, x), the solution of (2.1). Let n ∈ N and assume that v j (t, x) are constructed for j = 1 . . . , n -1. Let

F n (t, x) := n ℓ=1 L ℓ v n-ℓ (t, x) (4.11)
be the right-hand side in (4.8). Then v n is uniquely defined and given by the formula

v n (t) = t 0 P t-s F n (s)ds, t ≥ 0. (4.12)
By [H3], it is not difficult to check that v n , n ∈ N are smooth and that for all

t ≥ 0, v n (t) C k ≤ C(t) ϕ C k+4n , k ∈ N, (4.13)
where the constant C(t) depends on t, k, n and on the coefficient of the equation.

Clearly, this type of estimate can be improved in the elliptic case. This proves the first part of the Theorem.

To prove (i), we consider a fixed time t, and define the functions w n (x, s) := v n (t + s) for s ≥ 0 and n ∈ N. By definition, these functions satisfy the relation

∂ s w n (s, x) -Lw n (s, x) = n ℓ=1 L ℓ w n-ℓ (s, x), w n (0, x) = v n (t, x).
Let us consider the successive time derivatives of the functions w n (s, x). We have using the definition of w n , for all s ≥ 0,

∂ 2 s w n (s, x) = n ℓ=0 L ℓ ∂ s w n-ℓ (s, x) = n k=0 ℓ 1 +ℓ 2 =k L ℓ 1 L ℓ 2 w n-k (s, x),
and we see by induction that for all m ≥ 1 and s ≥ 0

∂ m s w n (s, x) = ℓ 1 +•••+ℓ m+1 =n L ℓ 1 • • • L ℓm w ℓ m+1 (s, x). (4.14)
Using the fact that the operators L ℓ are of order ℓ + 2 with no terms of order zero, we see that there exists a constant C depending on n and m, such that

∂ m s w n (s, x) ∞ ≤ C sup k=0,...,n |w k (s)| C 2k+2m .
Now let us consider the Taylor expansion of w n (τ ), for τ ≤ τ 0 and n = 0, . . . , N ,

w n (τ, x) = N -n m=0 τ m m! ∂ m t (0, x) + τ 0 s N -m (N -m)! ∂ N -m+1 t w n (s, x)ds = N -n m=0 τ m m! ℓ 1 +•••+ℓ m+1 =n L ℓ 1 • • • L ℓm w ℓ m+1 (0, x) + R N,n (τ, x).
Using the bounds on the time derivatives of w n (s, x), we obtain that for all τ ≥ 0 and all n = 0, . . . , N ,

R N,n (τ, x) ∞ ≤ Cτ N -m+1 sup s∈(0,τ ) n=0,...,N |w n (s, x)| C 4N+2 ,
for some constant depending on N , m. After summation in n, and using the expression (4.5) of the operators A n and the definition of w n , we get

v (N ) (t + τ, x) = N n=0 τ n n m=0 A m v n-m (t, x) + R N (t, τ, x) where R N (t, τ, x) ∞ ≤ C N τ N +1 sup s∈(0,τ ) n=0,...,N |v n (t + s, x)| C 4N+2 .
To conclude, we use (3.4) applied to ϕ = v (N ) (t, x), and we easily verify that Ev (N ) (τ, Xx (τ )) satisfy the same asymptotic expansion. The second estimate (ii) is then a consequence of (i) with t = 0 and (4.13).

Note that in the previous theorem, we have constructed a function v (N ) (t, x) which is an approximate solution of (4.7). More precisely, we can easily show that we have for all time t ≥ 0,

∂ t v (N ) (t, x) = L (N ) (τ ; x, ∂ x )v (N ) (t, x) + R (N ) (t, x), v (N ) (0, x) = ϕ(x), where R (N ) (t, x) = - ℓ 1 ,ℓ 2 =0,...,N ℓ 1 +ℓ 2 >N τ ℓ 1 +ℓ 2 L ℓ 1 v ℓ 2 (t, x)
is of order O(τ N +1 ).

5 Asymptotic expansion of the invariant measure and long time behavior

We now analyze the long time behavior of the solution of the modified equation (4.7). In the following, for a given operator B(x, ∂ x ), we denote by B(x, ∂ x ) * its adjoint with respect to the L 2 product. We start by an asymptotic expansion of a formal invariant measure for the numerical scheme.

Proposition 5.1 Let (L n ) n≥0 be the collection of operators defined recursively by (4.4). There exists a collection of functions (µ n (x)) n≥0 such that µ 0 (x) = ρ(x),

T d µ n (x)dx = 0 for n ≥ 1, and for all n ≥ 0,

L * 0 µ n = - n ℓ=1 (L ℓ ) * µ n-ℓ . (5.1)
Let N ≥ 0 be fixed and L (N ) (τ ; x, ∂ x ) the operator defined by (4.6). Then the function

µ (N ) (τ ; x) = ρ(x) + N n=1 τ n µ n (x) ∈ C ∞ (T d , R) satisfies T d µ (N ) (τ ; x)dx = 1, and
L (N ) (τ ; x, ∂ x ) * µ (N ) (τ ; x) = G (N ) (τ ; x),
with for all k and all τ ∈ [0, τ 0 ],

G (N ) (τ ; x) C k ≤ C N,k τ N +1 and T d G (N ) (τ ; x)dx = 0,
where C N,k depends on N and k.

Proof. Assume that µ 0 = ρ and µ j (x) are known, for j = 0, . . . , n -1 with n ≥ 1. Consider the equation (5.1) given by

L * 0 µ n = - n ℓ=1 (L ℓ ) * µ n-ℓ = G n .
Note that the right-hand side G n (x) is a smooth function satisfying

T d G n (x)dx = - n ℓ=1 T d (L ℓ ) * µ n-ℓ dx = - n ℓ=1 T d µ n-ℓ L ℓ 1dx = 0,
where 1 denotes the constant function equal to 1, which as already seen is in the Kernel of all the L ℓ .

Using Hypothesis [H2], we easily obtain the existence of a C ∞ function µ n satisfying (5.1) and T d µ n (x)dx = 0. This shows the first part of the Proposition. We then write

L (N ) (τ ; x, ∂ x ) * µ (N ) (τ ; x) = 2N n=0 τ n ℓ 1 +ℓ 2 =n ℓ i ≤N (L ℓ 1 ) * µ ℓ 2 = 2N n=N +1 τ n ℓ 1 +ℓ 2 =n ℓ i ≤N (L ℓ 1 ) * µ ℓ 2 =: G (N ) (τ ; x),
and we easily verify that G (N ) satisfies the hypothesis of the Proposition, owing to the fact that L ℓ 1 = 0 for all ℓ ≥ 0.

Proposition 5.2 For all n and k there exists a polynomial P k,n (t) such that for all t ≥ 0,

v n (t, x) - T d ϕdµ n C k ≤ P k,n (t)e -λt ϕ -ϕ C k+4n . (5.2)
Proof. Using the fact that µ 0 = ρ and v 0 = u, wee see that estimate (5.2) is satisfied for n = 0 (see Equation (2.5)). Let n ≥ 1 and assume that v j , j = 0, . . . , n -1 satisfy for k ∈ N, t ≥ 0:

v j (t, x) - T d ϕ(x)µ j (x)dx C k ≤ P k,j (t)e -tλ ϕ(x) -ϕ C k ,
for some polynomial P k,j . Let us set

c n (t) = n m=0 T d v n-m (t, x)µ m (x)dx.
We claim that c n (t) does not depend on time. Indeed:

n m=0 ∂ t T d v n-m (t, x)µ m (x)dx = n m=0 ∂ t T d v m (t, x)µ n-m (x)dx = n m=0 m ℓ=0 T d L m-ℓ v ℓ (t, x)µ n-m (x)dx = n ℓ=0 n m=ℓ T d v ℓ (t, x)L * m-ℓ µ n-m (x)dx = n ℓ=0 T d v ℓ (t, x) n-ℓ m=0 L * m µ n-ℓ-m (x)dx = 0,
by definition of the coefficients µ n , see (5.1). Note that, thanks to the smoothness properties of all the functions, the computation above is easily justified. We deduce:

T d ∂ t v n (t, x)ρ(x)dx = - n m=1 T d ∂ t v n-m (t, x)µ m (x)dx. ( 5.3)
Next, we compute the average of F n . By (4.8), (4.11) and (5.3), we have

F n (t) = T d F n (t, x)ρ(x)dx = T d ∂ t v n (t, x)ρ(x)dx - T d Lv n (t, x)ρ(x)dx = T d ∂ t v n (t, x)ρ(x)dx = - n m=1 T d ∂ t v n-m (t, x)µ m (x)dx.
We rewrite (4.12) as follows

v n (t, x) = t 0 F n (s) ds + t 0 P t-s (F n (s, x) -F n (s) )ds.
Using the previous expression obtained for F n (s) and recalling the initial data for v n , we deduce that Note that, since L ℓ , ℓ ∈ N is a differential operator of order 2ℓ + 2 with smooth coefficients and containing no zero order terms, we have The conclusion follows.

v n (t, x) = - n m=1 T d
We give now our main result concerning the long time behavior of the numerical solution:

Theorem 5.3 Let τ 0 and N be fixed. Then there exists C N and a polynomial P N (t) such that the following holds: Let X p be the discrete process defined by (2.6), then we have for p ≥ 0, τ ≥ 0 and smooth function ϕ(x)

∀ p ∈ N, E ϕ(X p ) -v (N ) (t p , x) ∞ ≤ C N τ N ϕ C 8N+2 , (5.4)
where for all p, t p = pτ . Moreover, we have

∀ p ∈ N, Eϕ(X p ) - T d ϕdµ (N )
∞ ≤ P N (t p )e -λtp + C N τ N , ϕ C 8N+2 , (5.5) where dµ (N ) (x) = µ (N ) (x)dx.

Proof. For all p, with t j = jτ , we have

E ϕ(X p ) -v (N ) (t p , x) = E v (N ) (0, X p ) -v (N ) (t p , x) = E p-1 j=0
E X p-j-1 v (N ) (t j , X p-j ) -v (N ) (t j+1 , X p-j-1 ) .

Here we have used the notation E X p-j-1 for the conditional expectation with respect to the filtration generated by X p-j-1 . By the Markov property of the Euler process at times t j : E X p-j-1 v (N ) (t j , X p-j ) -v (N ) (t j+1 , X p-j-1 ) = E X p-j-1 v (N ) (t j , XX p-j-1 (τ )) -v (N ) (t j+1 , X p-j-1 ) .

Using (4.10) with t = t j , and Proposition 5.2, we deduce that Q N (t j )e -λt j for some constant C N and polynomial Q N (t). We have used: |v n (t j +s, x)| 4N +2 = |v n (t j + s, x) -T d ϕdµ n | 4N +2 . We conclude by using the fact that for a fixed constant γ > 0, we have where the constant C depends on γ and τ 0 . This shows (5.4). The second estimate is a consequence of Proposition 5.2.

E ϕ(X p ) -v (N ) (t p , x) ∞ ≤ C N τ N +1

  and obtain(3.6) with N + 1 replaced by N + 2.

v 0 P 0 P

 00 n-m (t, x)µ m (x)dx + T d ϕ(x)µ n (x)dx + t t-s (F n (s, x) -F n (s) )ds. Then, using T d µ m (x)dx = 0, m ∈ N, we get v n (t, x)-T d ϕ(x)µ n (x)dx = n m=1 T d (v n-m (t, x)-T d ϕ(x)µ n-m (x)dx)µ m (x)dx + t t-s (F n (s, x) -F n (s) )ds.

F 0 p

 0 n (s) -F n (s) C k ≤ c k,ℓ n-1 ℓ=0 v ℓ (s) -T d v ℓ (s, x)µ ℓ (x)dx C k+2(n-ℓ)+2 .Then by [H3]v n (t, x) -T d ϕ(x)µ n (x)dx C k ≤ n m=1 v n-m (t, x) -T d ϕ(x)µ n-m (x)dx ∞ T d |µ m (x)|dx + t k (t -s)e -λ(t-s) F n (s, x) -F n (s) C k ds,and using the recursion assumptionv n (t, x) -T d ϕ(x)µ n (x)dx C k ≤ n m=1c n,m P 0,n-m (t)e -tλ ϕ(x) -ϕ ∞ + t -s)P k+2(n-ℓ)+2,ℓ (s)e -λt ds ϕ(x) -ϕ k+4n .

  τ ) n=0,...,N |v n (t j + s, x)| 4N +2 ≤ C N τ N +1 ϕ C 8N+2 p-1 j=0

p- 1 j=0e

 1 -γjτ ≤ 1 1 -e -γτ ≤ C τ