Jérémie Bigot
email: jeremie.bigot@isae.fr

Benjamin Charlier
email: benjamin.charlier@math.univ-toulouse.fr

Consistent estimation of a mean planar curve modulo similarities

Keywords: Shape spaces, Perturbation models, Consistency, Procrustes Means, Nonparametric inference, High-dimensional data, Linear smoothing. AMS classifications: Primary 62G08, secondary 62H11

We consider the problem of estimating a mean planar curve from a set of J random planar curves observed on a k-points deterministic design. We study the consistency of a smoothed Procrustean mean curve when the observations obey a deformable model including some nuisance parameters such as random translations, rotations and scaling. The main contribution of the paper is to analyze the influence of the dimension k of the data and of the number J of observed configurations on the convergence of the smoothed Procrustean estimator to the mean curve of the model. Some numerical experiments illustrate these results.

Introduction

In this paper, we are interested in the statistical treatment of random planar curves observed through rigid deformations of the plane. In many fields of interest, such data appear as, for instance, contours extracted from digital images or level sets of a real function defined on the plane. In handwriting recognition problems one typically compares curves that describe letters, digits or signatures and the acquisition process often create some ambiguity of location, size and orientation. Our aim is then to study an estimation procedure for a mean curve from a sample of J noisy and discretized planar curves observed through translation, rotations and scaling. The group generated by this set of transformations is usually called similarity group of the plane and in the sequel, an element of this group will be called a deformation.

A deformable model coming from statistical shape analysis

Deformable model In many practical cases of interest, data are collected through a computer device such as a digital camera or an image scanner. In order to take this into account in our model, we assume that observations at hand are discretized versions of continuous planar curves. Each observation Y is then given by a set of k points of the plane called a configuration. It can be written as a k × 2 real matrix

Y =    Y 1 . . . Y k    ∈ R k×2 ,
where Y ℓ = Y

(1) ℓ , Y

(2) ℓ ∈ R 2 with ℓ = 1, . . . , k. A degenerated configuration is a configuration composed of k times the same point b = (b (1) , b (2) ) ∈ R 2 . In tensorial notation, such a configuration is written 1 k ⊗ b ∈ R k×2 where 1 k denotes the (column) vector of R k with all entries equal to one and ⊗ denotes the tensor product. From now on, we assume that the observations satisfy the following functional regression model for j = 1, . . . , J,

Y j = e a * j (f + ζ j )R α * j + 1 k ⊗ b * j , with R α * j = cos(α * j ) -sin(α * j ) sin(α * j ) cos(α * j ) , (1.1)
where the unknown mean pattern f ∈ R k×2 has been obtained by sampling on an equi-spaced design a planar curve f = (f (1) , f (2) ) : [0, 1] -→ R 2 satisfying f (0) = f (1). It means that we have

f = f ( ℓ k ) k ℓ=1 := f (1) ( ℓ k ), f (2) ( ℓ k ) k ℓ=1 ∈ R k×2 .
The error terms ζ j ∈ R k×2 , j = 1, . . . , J are independent copies of a random perturbation ζ in R k×2 with zero expectation. For j = 1, . . . , J, the scaling, rotation and translation parameters (a * j , α * j , b * j ) ∈ R × [-π, π[×R 2 are independent and identically distributed (i.i.d) random variables independent of the random perturbations ζ j .

Model (1.1) is a deformable model in the sens of [Bigot & Charlier, 2011] and our aim is to estimate the underlying mean curve f : [0, 1] -→ R 2 from Y 1 , . . . , Y J ∈ R k×2 . More precisely, we study the influence on the estimation procedure of the number J of configurations at hand and the number k of discretization points composing the configurations Y j 's. Note that the model (1.1) has been introduced by Goodall in [Goodall, 1991] but with a fixed number k of points in each configuration and nonrandom nuisance parameters (a * j , α * j , b * j )'s. This framework has been highly popular in the statistical shape community, see [START_REF] Kent | [END_REF], Le, 1998, Lele, 1993]. Note that the configuration f was called a population mean in [Goodall, 1991] or a perturbation mean in [Huckemann, 2011].

Shape of a configuration. Since the seminal work of Kendall [Kendall, 1984], one considers that the shape of a configuration Y ∈ R k×2 is "what remains when translations, rotations and scaling are filtered out". More precisely, two configurations Y 1 , Y 2 ∈ R k×2 are said to have the same shape if there exists a vector (a, α, b)

∈ R × [-π, π[×R 2 such that Y 2 = e a Y 1 R α + 1 k ⊗ b, with R α = cos(α) -sin(α) sin(α) cos(α) , (1.2)
The Kendall's Shape Space is the quotient space modulo this equivalent relation. It is usually denoted by Σ k 2 and defined as the set of normalized (i.e. centered and scaled to size one) configurations quotiented by the rotation of the plane, see part 2.2 for further details.

The above definition can be trivially extended to the case of planar curves by replacing the configuration f ∈ R k×2 by the planar curve f : [0, 1] -→ R 2 . Thence, we call shape of f the set of planar curves that can be written e a f R α + 1 ⊗ b for some (a, α, b)

∈ R × [-π, π[×R 2 .
1.2 Estimation of the mean curve and of the shape of the mean curve Consistency in the Shape Space. In shape analysis, an important issue is the computation of a sample mean shape from a set of J random planar configurations Y 1 , . . . , Y J ∈ R k×2 satisfying model (1.1) and the study of its consistency as the number of samples J goes to infinity (k remaining fixed). According to Goodall [Goodall, 1991] a sample mean pattern f ∈ R k×2 computed from Y 1 , . . . , Y J is said to be consistent if, as J → ∞, it has asymptotically the same shape than the mean pattern f . In this framework, the deformations parameters are considered as nuisance parameters that contain no informations. That is why the data are first normalized (i.e. centered and scale to unit size) without changes in the statistical analysis. The study of consistent procedures to estimate the shape of the mean pattern f using this approach has been considered by various authors [START_REF] Kent | [END_REF], Le, 1998, Lele, 1993, Huckemann, 2011]. In this setting, sample mean patterns obtained by a Procrustes procedure have received a special attention. In particular, it is shown in [START_REF] Kent | [END_REF], Le, 1998] that, in the very specific case of isotropic perturbations ζ, the so-called full and partial Procrustes sample means are consistent estimators of the shape of f . Nevertheless, these estimators can be inconsistent for non-isotropic perturbations. Therefore, it is generally the belief that consistent statistical inference based on Procrustes analysis is restricted to very limited assumptions on the distribution of the data, see also [START_REF] Dryden | Statistical shape analysis[END_REF], Kendall et al., 1999, Huckemann, 2011] for further discussions.

Estimation of the mean curve. Our approach differ from the one developed in Shape analysis in two main aspects. First, we assume that the unknown mean pattern f ∈ R k×2 in model (1.1) has been obtained by sampling a planar curve f : [0, 1] -→ R 2 on an equi-spaced design. This allows us to study the influence of the number k of points composing each configuration on the estimation of f . Note that if f is sufficently regular and k is large, to estimate f ∈ R k×2 (i.e. the values of f on the design) roughly amounts to estimate f . The second difference concerns the randomness of the deformation parameters

(a * j , α * j , b * j ) ∈ R × [-π, π[×R 2 , j = 1, .
. . , J, that we want to estimate rather than to consider as nuisance parameters. Indeed, the values of the deformations parameters can be informative in some cases : assume that the size of the data are very similar except for one. This difference in size may be due to a (possibly) relevant factor and if the data are normalized as in Shape analysis this information is lost.

Under a suitable smoothness assumption on f , we are able to estimate consistently, with an asymptotic in k only, the curve f j := e a * j f R α * j +1⊗b * j from Y j for each j = 1, . . . , J. Note that, by definition, the f j 's have the same shape as f . When J is fixed, this smoothing step allows us to estimate consistently the deformations parameters (a * j , α * j , b * j ) J j=1 with a Procrustes matching step. Now, if we assume that the deformations parameters (a * j , α * j , b * j ) are centered random variables, we show that it is possible to recover the true mean shape f when min{J, k} → ∞. These results are rigorously stated in the next section, see Theorem 1.1 below.

Main contribution

Our estimating procedure is composed of two steps. First, we perform a dimension reduction step by projecting the data into a low-dimensional space of R k×2 to eliminate the influence of the random perturbations ζ j . Then, in a second step, we apply Procrustes analysis in this low-dimensional space to obtain a consistent estimator of f ∈ R k×2 .

The reduction dimension step is based on an appropriate smoothness assumptions on f . Let L > 0, s > 0 and define the Sobolev ball of radius L and degree s as

H s (L) = f = (f (1) , f (2) ) ∈ L 2 ([0, 1], R 2 ), m∈Z (1 + |m| 2 ) s (|c m (f (1) )| 2 + |c m (f (2) )| 2 ) < L ) (1.3) where c m (f ) = (c m (f (1) ), c m (f (2) )) = 1 0 f (1) (t)e -i2πmt dt, 1 0 f (2) (t)e -i2πmt dt ∈ C 2 is the m-th Fourier coefficient of f = (f (1) , f (2) ) ∈ L 2 ([0, 1], R 2 ), for m ∈ Z. Assumption 1. The curve f : [0, 1] -→ R 2 is closed, i.e. f (0) = f (1) ∈ R 2 ,
and belongs to H s (L) for some L > 0 and s > 3 2 . Moreover, the k

× 2 matrix f = f ( ℓ k )
k ℓ=1 is of rank two, i.e. the configuration f is not degenerated.

Assumption 1 implies that f is not reduced to a point, continuously differentiable and is equal to its Fourier series. We introduce the following k × k matrix

A λ = 1 k 0≤|m|≤λ e i2πm ℓ-ℓ ′ k k ℓ,ℓ ′ =1
.

(1.4)

The matrix A λ is the smoothing matrix corresponding to a discrete Fourier low pass filter with frequency cutoff λ ∈ N. It is a projection matrix in a sub-space V λ of R k of dimension 2λ + 1.

Then, we project the data on V λ × V λ ⊂ R k×2 , and we estimate the scaling, rotation and translation parameters in model (1.1) using M-estimation as follows: denote the scaling parameters by a = (a 1 , . . . , a J ) ∈ R J , the rotation parameters by α = (α 1 , . . . , α J ) ∈ R J and the translation parameters by b = (b 1 , . . . , b J ) ∈ R 2J , and introduce the functional,

M λ (a, α, b) = 1 Jk J j=1 e -a j A λ (Y j -1 k ⊗ b j )R -α j - 1 J J j ′ =1 e -a j ′ A λ (Y j ′ -1 k ⊗ b j ′ )R -α j ′ 2 R k×2 , (1.5)
where

• R k×2 is the standard Euclidean norm in R k×2 . An M-estimator of a * , α * , b * = a * 1 , . . . , a * J , α * 1 , . . . , α * J , b * 1 , . . . , b * J ∈ R J × [-π, π[ J ×R 2J is given by (â λ , αλ , bλ ) ∈ argmin (a,α,b)∈Θ 0 M λ (a, α, b), (1.6)
where âλ , αλ , bλ = âλ 1 , . . . , âλ

J , αλ 1 , . . . , αλ J , bλ 1 , . . . , bλ J ∈ R J × [-π, π[ J ×R 2J and Θ 0 = (a, α, b) ∈ [-A, A] J × [-A, A] J × R 2J : J j=1 a j = 0, J j=1 α j = 0 and J j=1 b j = 0 , (1.7)
with A, A > 0 being parameters whose values will be discussed below.

Finally, the mean pattern f is estimated by the following smoothed Procrustes mean

f λ = 1 J J j=1 e -â λ j A λ Y j -1 k ⊗ bλ j R -αλ j . (1.8)
To analyze the convergence of the estimator f λ to the mean pattern f , let us introduce some regularity conditions on the covariance structure of the random variable

ζ in R k×2 . Let ζ = (ζ (1) , ζ (2) ) = (ζ (1) 1 , . . . , ζ (1) k , ζ (2) 1 , . . . , ζ (2) k ) ∈ R 2k be the vectorized version of ζ = (ζ (1) ℓ , ζ (2) ℓ ) k ℓ=1 ∈ R k×2 . Assumption 2. The random variable ζ is a centered Gaussian vector in R 2k with covariance matrix Σ k ∈ R 2k×2k . Let γ max (k) be the largest eigenvalue of Σ k . Then, lim k→∞ γ max (k)k -2s 2s+1 = 0,
where s is the smoothness parameter defined in Assumption 1.

For example, the isotropic Gaussian error model corresponds to Σ k = Id k and γ max (k) = 1 satisfies Assumption 2. If there exists correlations terms (i.e. non-zero off-diagonal entries in Σ k ), then the level of the perturbation ζ has to be sufficiently small. A simple model is the case where [Gray, 2006]) and thus Assumption 2 is again satisfied. The following theorem is the main result of the paper.

Σ k = [S(|ℓ -ℓ ′ |)] k ℓ,ℓ ′ =1 for some function S : R -→ R satisfying R |S(t)| dt < +∞ implying that γ max (k) ≤ ℓ∈Z |S(ℓ)| (see Lemma 4.11 in
Theorem 1.1. Consider model (1.1) and suppose that Assumptions 1 and 2 hold. Suppose also that the random variables (a * , α * , b * ) are bounded and belong to

[-A 2 , A 2 ] J × [-A 2 , A 2 ] J × [-B, B] 2J for some 0 < A, B and 0 < A < π 4 and that λ(k) = k 1 2s+1 .
• For any J ≥ 2 there exists

f Θ 0 = e a 0 f R α 0 + 1 k ⊗ b 0 for some (a 0 , α 0 , b 0 ) ∈ R × [-π, π[×R 2 and a function V 1 (k, x
) such that for any x > 0,

P 1 k f λ -f Θ 0 2 R k×2 ≥ V 1 (k, x) ≤ e -x ,
(1.9)

with V 1 (k, x) → 0 when k → ∞ and x remains fixed.

• Suppose, in addition, that the random variables

(a * , α * , b * ) have zero expectation in [-A 2 , A 2 ] J × [-A 2 , A 2 ] J × [-B, B] 2J with A, A < 0.1.
Then, there exists a function V 2 (J, k, x) such that for any x > 0, (1.10) where V 2 (J, k, x) → 0 when min{J, k} → ∞ and x remains fixed.

P 1 k f λ -f 2 R k×2 ≥ V 2 (J, k, x) ≤ e -x ,
Statement (1.9) means that, under mild assumptions on the covariance structure of the error terms ζ j , it is possible to consistently estimate the shape of the mean curve f when the number of observations J is fixed and the number k of discretization points increases. Note that (a 0 , α 0 , b 0 ) depends on J and is given by formula (3.3) in Section 3. The function V 1 (k, x) is explicitly given in Section 5.2. To obtain statement (1.10), we assume the condition A, A < 0.1 which means that the random scaling and rotations in model (1.1) are not too large. Also, it is assumed that random scaling, rotations and translations have zero expectation, meaning that the deformations parameters in model (1.1) are centered around the identity. Then, under such assumptions, statement (1.10) shows that one can consistently estimate the true mean curve f when both the sample size J and the number of landmarks k go to infinity. Again, the function V 2 (J, k, x) is explicitly given in Section 5.2. These results are consistent with those obtained in [Bigot & Charlier, 2011], where we have studied the consistency of Fréchet means in deformable models for signal and image processing.

Organization of the paper

In Section 2, we recall some properties on the similarity group of the plane, and we describe its action on the mean pattern f . Then, we discuss General Procrustes Analysis (GPA) and we compare it to our approach. In Section 3 we discuss some identifiability issues in model (1.1). The estimating procedure is described in detail in Section 4. Consistency results are given in Section 5. Some experiments in Section 6 illustrate the numerical performances of our approach. All the proofs are gathered in a technical appendix.

2 Group structure and Generalized Procrustes Analysis

The similarity group

Group action First let us introduce some notations and definitions that will be useful throughout the paper. The similarity group of the plane is the group (G, .) generated by isotropic scaling, rotations and translations. The identity element in G is denoted e and the inverse of g ∈ G is denoted by g -1 . We parametrize the group G by a scaling parameter a ∈ R, an angle α ∈ [-π, π[ and a translation b ∈ R 2 , and we make no difference between g ∈ G and its parametrization

(a, α, b) ∈ R × [-π, π[×R 2 . For all g 1 = (a 1 , α 1 , b 1 ), g 2 = (a 2 , α 2 , b 2 ) ∈ R × [-π, π[×R 2 we have g 1 .g 2 = (a 1 , α 1 , b 1 ).(a 2 , α 2 , b 2 ) = (a 1 + a 2 , α 1 + α 2 , e a 1 b 2 R α 1 + b 1 ), g -1 1 = (a 1 , α 1 , b 1 ) -1 = (-a 1 , -α 1 , -e -a 1 b 1 R -α 1 ),
(2.1) e = (0, 0, 0).

The action of G onto R k×2 is given by the mapping (g, x) -→ g.x := e a xR α +1 k ⊗b, for g = (a, α, b) ∈ G and x ∈ R k×2 . Note that we use the same symbol "." for the composition law of G and its action on R k×2 . This action can also be defined on L 2 ([0, 1], R 2 ) by replacing f by f ∈ L 2 ([0, 1], R 2 ). Coming back to the final dimensional case, let

1 k×2 = 1 k ⊗ R 2
be the two dimensional linear subspace of R k×2 consisting of degenerated configurations, i.e. configurations composed of k times the same landmarks. The orthogonal subspace 1 ⊥ k×2 is the set of centered configurations. We have the orthogonal decomposition R k×2 = 1 ⊥ k×2 ⊕ 1 k×2 , and for any configuration x ∈ R k×2 we write x = x 0 + x ∈ 1 ⊥ k×2 ⊕ 1 k×2 . We call x 0 the centered configuration of x and

x = 1 k ⊗ 1 k k ℓ=1 x (1) ℓ , 1 k k ℓ=1 x
(2) ℓ the degenerated configuration associated to x, see Figure 1 for an illustration.

Orbit, stabilizer and section Given a configuration x in R k×2 , the orbit of x is defined as the set

G.x := {g.x, g ∈ G} ⊂ R k×2 .
This set is also called the shape of x. The orbit of any degenerated configuration is the entire subspace 1 k×2 . Note also that the linear subspace 1 k×2 is stable by the action of G, and that the action of G on 1 k×2 is not free, meaning that for any x ∈ 1 k×2 the equality g 1 .x = g 2 .x does not imply that

g 1 = g 2 . Now, if x ∈ R k×2 \ 1 k×2 is a non-degenerated configuration of k landmarks, its orbit G.x is a sub-manifold of R k×2 \ 1 k×2 of dimension dim(G) = 4.
Given a configuration x ∈ R k×2 , the stabilizer I(x) is the closed subgroup of G which leaves x invariant, namely 1) , x(2) ) ∈ 1 k×2 is a degenerated configuration, its stabilizer is non trivial and is equal to

I(x) = {g ∈ G : g.x = x}. If x = 1 k ⊗ (x (
I(x) = {(a, α, (x (1) , x(2) ) -e a (x (1) , x(2) )R α ), a ∈ R, α ∈ [-π, π[}. If x ∈ R k×2 \ 1 k×2
is a nondegenerated configuration, its stabilizer I(x) is trivial, i.e. is reduced to the identity {e}. The action of G is said free if the stabilizer of any point is reduced to the identity. Hence, the action of G is free on the set of non-degenerated configurations of k-ads in R 2 . A section of the orbits of G is a subset of R k×2 containing a unique element of each orbit. A wellknown example of section for the similarity group acting on R k×2 \ 1 k×2 is the so-called Bookstein's coordinates system (see e.g. [START_REF] Dryden | Statistical shape analysis[END_REF]] p. 27).

Kendall's shape space and Generalized Procrustes analysis

Shape space Let x ∈ R k×2 \ 1 k×2 be a non-degenerated configuration. Let H = Id k -1 k 1 k 1 ′ k
be a centering matrix. The effect of translation can be eliminated by centering the configuration x using the matrix H (see [START_REF] Dryden | Statistical shape analysis[END_REF] for other centering methods), while the effect of isotropic scaling is removed by projecting the centered configuration on a unit sphere, which yields to the so-called pre-shape x 0 of x defined as

x 0 = Hx Hx R k×2 ∈ R k×2 .
Consider now the pre-shape sphere defined by S k 2 := x 0 , x ∈ R k×2 \ 1 k×2 and see Figure 1 for an illustration. Note that this normalization of the planar configurations amounts to choose a section for

1 ⊥ k×2 G.x 1 k×2 x x x0 x 0
Figure 1: Three orbits of the action of the similarity group G are represented in blue. The space of centered configurations 1 ⊥ k×2 is the green plane. The preshape sphere S k 2 is the red circle. For a particular x ∈ R k×2 , the centered version is x 0 and the centered and normalized version is x 0 . The degenerated configuration associated to x is x.

the action of the group generated by the translation and scaling in the plane. The Kendall's shape space is then defined as the quotient of S k 2 by the group SO(2) of rotations of the plane, namely

Σ k 2 := S k 2 /SO(2) = [x 0 ] : x 0 ∈ S k 2 with [x 0 ] = x 0 R α , α ∈ [-π, π[ .
The space Σ k 2 can be endowed with a Riemannian structure and we refer to [START_REF] Kendall | Shape and shape theory[END_REF] for a detail discussion on its geometric properties.

Let us briefly recall the definition of the so-called partial and full Procrustes distances on S k 2 . The partial Procrustes distance is defined on the pre-shape sphere S k 2 as

d 2 P (x 0 , y 0 ) = inf α∈[-π,π[ x 0 -y 0 R α 2 R k×2 , x, y ∈ S k 2 .
Hence, it is the (Euclidean) distance between the orbits [x 0 ] = SO(2).x 0 and [y 0 ] = SO(2).y 0 with x 0 , y 0 ∈ S k 2 . Let now H be the group of transformations of the plane generated by scaling and rotations. The action of h ∈ H on the centered configuration x 0 is defined as h.x 0 := e a x 0 R α where h = (a, α) ∈ R × [-π, π[. The Full Procrustes distance is then defined as

d 2 F (x 0 , y 0 ) = inf h∈H x 0 -h.y 0 2 R k×2 , x, y ∈ S k 2 .
Generalized Procrustes analysis The full Procrustes sample mean ŶF of Y 1 , . . . , Y J (see e.g. [Goodall, 1991[START_REF] Dryden | Statistical shape analysis[END_REF]) is defined by

ŶF = argmin x 0 ∈S k 2 J j=1 d 2 F (Y 0 j , x 0 ). (2.2)
The partial Procrustes mean ŶP is defined in the same way by replacing d F by d P in (2.2). Thence, this two Procrustes means are Fréchet mean either on (S k 2 , d F ) or (S k 2 , d P ) endowed with the empirical measure µ J = J j=1 δ Y j . In practice, there are several way to compute the full Procrustes mean. In [Kent, 1992] the author used complex coordinates and expressed the full Procrustes mean as the biggest eigenvalue of a symmetrical positive definite complex matrix, see [START_REF] Dryden | Statistical shape analysis[END_REF]] result 3.2 and [Bhattacharya & Patrangenaru,2 The full Procrustes mean ŶF can also be approximated by General Procrustes procedure which amounts to use the following identity

ŶF = 1 J J j=1 ĥj • Y 0 j (2.3)
where ĥ1 , . . . , ĥJ are the argmins of the functional

M (h 1 , . . . , h J ) = 1 J J j=1 h j .Y 0 j -1 J J j ′ =1 h j ′ .Y 0 j ′ 2 R k×2 subject to the constraint 1 J J j=1 h j .Y 0 j 2 R k×2 = 1.
The configurations ĥ.Y j 's are known as Full Procrustes fits and the ĥ can be explicitly computed by using a singular value decomposition. In practice, one can use the iterative General Procrustes algorithm to compute Ŷ, see [START_REF] Dryden | Statistical shape analysis[END_REF]] pages 90-91.

Discussion on the double asymptotic setting and comparison with GPA

Asymptotic settings When random planar curves (such as digits or letters for instance) are observed, a natural framework for statistical inference is an asymptotic setting in the number J of curves. This setting means that increasing the number of curves at hand should help to compute a more accurate empirical mean curve. Unfortunately, consistency results of Procrustes type procedures are reduced to the very specific case of an isotropic perturbation. In this paper, we show that increasing the number k of discretization points will ensure a consistent estimation of a mean shape in more general cases.

Consider model (1.1) where k > 2 is fixed and the random perturbation ζ is isotropic, (see Proposition 1 in [Le, 1998] for a precise definition of isotropy for random variables belonging to R k×2 ). In this framework, it has been proved in [START_REF] Kent | [END_REF] that the functional

M J (x 0 ) = 1 J J j=1 d 2 F (Y j , x 0 ) defined on R k×2 converge uniformly in probability to the functional M (x 0 ) = E(d 2 F (Y 1 , x 0
)) which admits a unique minimum at [f 0 ]. These two facts imply that [ ŶF ] converges almost surely to [f 0 ] as J → +∞. In [Le, 1998] the author used a Fréchet mean approach to show the consistency of ŶF with a slightly more general kind of noise. Finally, note that the estimator ŶP defined above is also studied in [START_REF] Kent | [END_REF], Le, 1998] with similar consistency results.

When the random perturbation ζ in model (1.1) is non-isotropic, it has been argued in [START_REF] Kent | [END_REF]] that the Procrustes estimator ŶF can be arbitrarily inconsistent when the signal-to-noise ratio decrease. The heuristic presented by the authors suggests that the main phenomenon that prevent the Procrustes estimator to be consistent is the fact that the functional M do not attain its minimum at [f 0 ]. In section 4.2 of [Huckemann, 2011] the author makes this remark clear as he gives an explicit example : given a mean pattern f 0 , the idea is to increase the level of noise ζ until the argmin of the functional M which was initially equals to [f 0 ] jumps abruptly to another point. This phenomenon seems to be linked to the geometry of the sphere and properties of the Fréchet mean.

Comparison with GPA Hence, it is commonly the belief that Procrustes sample means can be inconsistent when considering convergence in Σ k 2 and the asymptotic setting J → +∞. Nevertheless, the above discussion suggests that a sufficient condition to ensure the consistency of Procrutres type estimators is to control the level of non-isotropic noise. That is why we introduced a pre-smoothing step that takes advantage of increasing the number k of points composing each configuration in order to ensure more general consistency results. Therefore, our approach and GPA share some similarities. They are both based on the estimation of scaling, rotation and translation parameters by a Procrustean procedure which leads to the M-estimators (1.6) which is related to (2.2). To compute a sample mean shape, this M-estimation step is then followed by a standard empirical mean in R k×2 of the aligned data using these estimated deformation parameters, see equations (1.8) and (2.3).

However, one of the main differences between the approach developed in this paper and GPA is the choice of the normalization of the data. In GPA, the deformation parameters ĥ1 , . . . , ĥJ are computed so that the full Procrustes sample mean ŶF belongs to the pre-shape sphere S k 2 , see the constraint appearing in (2.2). Therefore, the computation of ĥ1 , . . . , ĥJ is somewhat independent of any assumption on the true parameters (a * , α * , b * ) in model (1.1). In this paper, to ensure the wellposedness of the problem (1.6), we chose to compute the estimator (â λ , αλ , bλ ) by minimizing the matching criterion (1.5) on the constrained set Θ 0 . The choice of the constraints in Θ 0 is motivated by the hypothesis that the true deformation parameters (a * , α * , b * ) in (1.1) are centered around the identity.

Identifiability conditions

Recall that in model (1.1), the random deformations acting on the the mean pattern f are parametrized by a vector (a

* , α * , b * ) = (a * 1 , . . . , a * J , α * 1 , . . . , α * J , b * 1 , . . . , b * J ) in R J × [-π, π[ J ×R 2J . Assumption 3. Let 0 < A, B and 0 < A < π be three real numbers. The deformation parameters (a * j , α * j , b * j ), are i.i.

d random variables with zero expectation and and taking their values in

Θ * = - A 2 , A 2 × - A 2 , A 2 × [-B, B] 2 . Let Θ * = -A 2 , A 2 J × -A 2 , A 2 J × [-B, B] 2J . Under Assumption 3, we have (a * , α * , b * ) ∈ Θ * .
Note that the compactness of Θ * (and thus of Θ * ) is an essential condition to ensure the consistency of our procedure. Indeed, the estimation of the deformation parameters (a * , α * , b * ) and the mean pattern f is based on the minimization of the criterion (1.5). If there were no restriction on the amplitude of the scaling parameter, the degenerate solution a j = -∞ for all j = 1, . . . , J is always a minimizer of (1.5). Therefore, the minimization has to be performed under additional compact constraints.

The deterministic criterion

D Let (a, α, b) ∈ R J × [-π, π[ J ×R 2J
and consider the following criterion,

D(a, α, b) = 1 Jk J j=1 (g -1 j .g * j ).f - 1 J J j ′ =1 (g -1 j ′ .g * j ′ ).f 2 R k×2 . (3.1)
where g j = (a j , α j , b j ) and g * j = (a * j , α * j , b * j ) for all j = 1, . . . , J. The criterion D is a version without noise of the criterion M λ defined at (1.5). The estimation procedure described in Section 1.1 is based on the convergence of the argmins of M λ toward the argmin of D when k goes to infinity. As a consequence, choosing identifiability conditions amounts to fix a subset Θ 0 of R J × [-π, π[ J ×R 2J on which D has a unique argmin. In the rest of this section, we determine the zeros of D, and then we fix a convenient constraint set Θ 0 that contains a unique point at which D vanishes.

The criterion D clearly vanishes at

(a * , α * , b * ) ∈ R J × [-π, π[ J ×R 2J
. This minimum is not unique since easy algebra implies that

D(a, α, b) = 0 ⇐⇒ (g -1 j .g * j ).f = (g -1 j ′ .g * j ′ ).
f , for all j, j ′ = 1, . . . , J.

Suppose now that f / ∈ 1 2×k is a non-degenerated planar configuration. In Section 2.1, we have seen that the action of G on f is free, that is, the stabilizer I(f ) is reduced to the identity. Thus, we obtain,

D(a, α, b) = 0 ⇐⇒ g -1 j .g * j = g -1 j ′ .g * j ′ , for all j, j ′ = 1, . . . , J ⇐⇒      a * j -a j = a * j ′ -a j ′ , α * j -α j = α * j ′ -α j ′ , e -a j (b * j -b j )R -α j = e -a j ′ (b * j ′ -b j ′ )R -α j ′ for all j, j ′ = 1, . . . , J
We have proved the folowing result,

Lemma 3.1. Let f ∈ R k×2 be a non-degenerated configuration of k-ads in the plane, i.e. f / ∈ 1 2×k . Then, D(a, α, b) = 0 if and only if (a, α, b) belongs to the set (a * , α * , b * ) * G = (a * , α * , b * ) * (a 0 , α 0 , b 0 ), (a 0 , α 0 , b 0 ) ∈ R × [-π, π[×R 2 , where (a * , α * , b * ) * (a 0 , α 0 , b 0 ) = (a * 1 +a 0 , . . . , a * J +a 0 , α * 1 +α 0 , . . . , α * J +α 0 , e a * 1 b 0 R α * 1 +b * 1 , . . . , e a * J b 0 R α * J + b * J ) ∈ R J × [-π, π[ J ×R 2J .
Remark 1. Lemma 3.1 is simpler than it appears. By reordering the entries of the vector (a * , α * , b * ) there is an obvious correspondence between (a * , α * , b * ) ∈ Θ * and (g * 1 , . . . , g * J ) ∈ G J via the parametrization of the similarity group defined in Section 2.1. Hence, Lemma 3.1 tells us that the criterion D vanishes for all the vectors (a, α, b) ∈ R J × [-π, π[ J ×R 2J corresponding to the subset of the group G J given by

(g * 1 , . . . , g * J ) * G = {(g * 1 .g 0 , . . . , g * J .g 0 ), g 0 ∈ G} ⊂ G J .
The " * " notation is nothing else than the right composition by a same g 0 ∈ G of all the entries of a (g 1 , . . . , g J ) ∈ G J . Hence the subset (g * 1 , . . . , g * J ) * G can be interpreted as the orbit of (g * 1 , . . . , g * J ) ∈ G J under the (right) action of G. Indeed, G acts naturally by (right) composition on the all the coordinates of an element of G J .

The constraint set Θ 0

By Lemma 3.1, the set Θ 0 must intersect at a unique point, say

(a * Θ 0 , α * Θ 0 , b * Θ 0 ), each set (a * , α * , b * ) * G. It is convenient to choose Θ 0 to be of the form Θ 0 = R J × [-π, π[ J ×R 2J ∩ L 0 where L 0 is a linear space of R 4J .
The linear space L 0 must be chosen so that for any

(a * , α * , b * ) in Θ * , there exists a unique point (a * Θ 0 , α * Θ 0 , b * Θ 0 ) in Θ 0 that can be written as (a * Θ 0 , α * Θ 0 , b * Θ 0 ) = (a * , α * , b * ) * (a 0 , α 0 , b 0 ) for some (a 0 , α 0 , b 0 ) ∈ R × [-π, π[×R 2 .
Remark 2. As we have seen in Remark 1, the set (a * , α * , b * ) * G can be interpreted as an orbit of the action of G on G J . In this terminology, the set Θ 0 can be viewed as a section of the orbits. Indeed, the section is the set of representatives (a * Θ 0 , α * Θ 0 , b * Θ 0 ) of each orbit. See Figure 2 for an illustration. Let us consider a choice of Θ 0 motivated by the fact that, under Assumption 3, the random deformation parameters have zero expectation. In this setting, it is natural to impose that the estimated deformation parameters sum up to zero by choosing L 0 = 1 ⊥ 4J , which is the orthogonal of the linear space 1 4J = 1 4J .R ⊂ R 4J . Such a choice leads to the set Θ 0 defined equation (1.7) i.e.

Θ 0 = {(a, α, b) ∈ Θ J , (a 1 + . . . + a J , α 1 + . . . + α J , b 1 + . . . + b J ) = 0}. Now, let us show that for any (a * , α * , b * ) ∈ Θ * there exists a unique (a * Θ 0 , α * Θ 0 , b * Θ 0 ) = (a * , α * , b * ) * (a 0 , α 0 , b 0 ) ∈ Θ 0 for some (a 0 , α 0 , b 0 ) ∈ R × [-π, π[×R 2
. This amounts to solve the following equations

   a * 1 + a 0 + . . . + a * J + a 0 = 0, α * 1 + α 0 + . . . + α * J + α 0 = 0, e a * 1 b 0 R α * 1 + b * 1 + . . . + e a * J b 0 R α * J + b * J = 0. (3.2)
After some computations, we obtain that equations (3.2) are satisfied if and only if

(a 0 , α 0 , b 0 ) = (-ā * , -ᾱ * , - b * (e a * R α * ) -1 ), (3.3) where ā * = 1 J J j=1 a * j ∈ R, ᾱ * = 1 J J j=1 α * j ∈ R, b * = 1 J J j=1 b * j ∈ R 2 and e a * R α * = 1 J J j=1 e a * j R α * j is a 2 × 2 invertible matrix. Therefore, (a * Θ 0 , α * Θ 0 , b * Θ 0 ) is uniquely given by ([a * Θ 0 ] j , [α * Θ 0 ] j , [b * Θ 0 ] j ) = (a * j -ā * , α * j -ᾱ * , b * j -e a * j ( b * (e a * R α * ) -1 )R α * j ), for j = 1, . . . , J. Θ1 Θ * (a * , α * , b * ) (a * Θ 1 , α * Θ 1 , b * Θ 1 ) (a * Θ 0 , α * Θ 0 , b * Θ 0 ) Θ0 (a * , α * , b * ) * G
Figure 2: Choice of identifiability conditions when J = 2.

Remark 3. Another possible approach is to fix, say the first observation as a reference, meaning that the criterion D could be optimized on the following subspace of

R J × [-π, π[ J ×R 2J Θ 1 = {(a, α, b) ∈ [-A, A] J × [-A, A] 2J × R 2 , (a 1 , α 1 , b 1 ) = 0}.
With such a choice, for any

(a * , α * , b * ) ∈ Θ J , the j-th coordinate of (a * Θ 1 , α * Θ 1 , b * Θ 1 ) = (a * , α * , b * ) * (a 0 , α 0 , b 0 ) is given by ([a * Θ 1 ] j , [α * Θ 1 ] j , [b * Θ 1 ] j ) = (a * j -a * 1 , α * j -α * 1 , b * j -e a * j -a * 1 b * 1 R α * j -α * 1 ),
where

(a 0 , α 0 , b 0 ) = (a * 1 , α * 1 , b * 1 ) -1 = (-a * 1 , -α * 1 , -e -a * 1 b * 1 R -α * 1 )
. A graphical illustration of the choice of identifiability conditions for J = 2 is given in Figure 2.

4 The estimating procedure

A dimension reduction step

We use Fourier filtering to project the data into a low-dimensional space as follows. Assume for convenience that k is odd. For

x = (x ′ 1 , . . . , x ′ k ) ′ ∈ R k×2 and m = -k-1 2 , . . . , k-1 2 , let c m (x) = k ℓ=1 x ℓ e -i2πm ℓ k = k ℓ=1 x (1) ℓ e -i2πm ℓ k , k ℓ=1 x (2) ℓ e -i2πm ℓ k ∈ C 2 with x ℓ = (x (1) ℓ , x (2)
ℓ ), be the m-th (discrete) Fourier coefficient of x. Let λ ∈ {1, . . . , k-1 2 } be a smoothing parameter, and define for each Y j the smoothed shapes

f λ j = 1 k 0≤|m|≤λ c m (Y j )e i2πm ℓ k k ℓ=1 = A λ Y j ∈ R k×2 .
In Section 2.1, we have shown that the similarity group is not free on the subset 1 k×2 of degenerated configurations composed of k identical landmarks, see Section 2.1. That is why we are going to treat separately the subspace 1 k×2 and 1 ⊥ k×2 by considering the matrices

Ā = 1 k 1 k 1 ′ k and A λ 0 = 1 k 0<|m|≤λ e i2πm ℓ-ℓ ′ k k ℓ,ℓ ′ =1 . (4.1)
Remark that Ā is a projection matrix on the one dimensional sub-space

V := 1 k .R = {c1 k : c ∈ R} of R k . The matrix A λ 0 is a projection matrix in a (trigonometric) sub-space V λ 0 of dimension 2λ. Note that it is included in the linear space V⊥ = {x ∈ R k : x ′ 1 k = 0}. Hence, V λ 0 × V λ 0 is a linear subspace of 1 ⊥
k×2 which is the space of the centered configurations. We have,

A λ = A λ 0 + Ā and V λ = V λ 0 ⊕ V. (4.2)
Thus, we can write the smoothed shape f λ j as

f λ j = A λ Y j = A λ 0 Y j + ĀY j ∈ ( V × V) ⊕ (V λ 0 × V λ 0 )
where

V × V = 1 k×2 and V λ 0 ×V λ 0 ⊂ 1 ⊥ k×2 .
In other words, A λ 0 Y j is the smoothed centered configuration associated to Y j and ĀY j is the degenerated configuration given by the Euclidean mean of the k landmarks composing Y j . Finally, remark that the low pass filter and the action of similarity group commute, that is, we have for all g ∈ G and

f ∈ R k×2 g.(A λ f ) = e a A λ f R α + 1 k ⊗ b = A λ (e a f R α + 1 k ⊗ b) = A λ (g.f ).

Estimation of the deformation parameters

Recall that the estimator (â λ , αλ , bλ ) of (a * , α * , b * ) is defined by the optimization problem (1.6). Nevertheless, as it is suggested by the discussion of Sections 2.1 and 4.1, one can carry out the estimation process in two steps. First, we estimate the rotation and scaling parameters on the space

V λ 0 ×V λ 0 ⊂ 1 ⊥ k×2
of the centered configurations. We then use these estimators to estimate the translation parameters which act on V × V = 1 k×2 . Note that this procedure is equivalent to the optimization problem (1.6) as shown by Lemma 4.1 below.

Estimation of rotations and scaling. Define

M λ 0 (a, α) = 1 Jk J j=1 e -a j A λ 0 Y j R -α j - 1 J J j ′ =1 e -a j ′ A λ 0 Y j ′ R -α j ′ 2 R k×2
.

Let 1 ⊥ J = {(a 1 , . . . , a J ) ∈ R J , a 1 + . . . , a J = 0} and consider the space Θ a,α

0 = ([-A, A] J ∩ 1 ⊥ J ) × ([-A, A] J ∩ 1 ⊥ J ).
Then, estimators of the rotation and scaling parameters are given by

(â λ , αλ ) ∈ argmin (a,α)∈Θ a,α 0 M λ 0 (a, α). (4.3)
Estimation of translations. Now that we have computed estimators of the rotation and scaling parameters, let us define the criterion,

M (a, α, b) = 1 Jk J j=1 e -a j Ā(Y j -1 k ⊗ b j )R -α j - 1 J J j ′ =1 e -a j ′ Ā(Y j ′ -1 k ⊗ b j ′ )R -α j ′ 2 R k×2
and the space

Θ b 0 = {(b (1) 1 , . . . , b (1)
J , b

(2) 1 , . . . , b

(2) (4.4)

J ) ∈ R 2J , b (
We emphasis that the estimators of the translation parameters depend on the estimated rotation and scaling parameters. It is shown in the proof of Lemma 4.1 below that we have an explicit expression of bλ given by bλ = ( Ȳ1 -e a 1 Ȳ(e a R α )

-1 R α 1 , . . . , ȲJ -e a J Ȳ(e a R α ) -1 R α J ) (4.5)
where

e a R α = 1 J J j=1 e a j R α j ∈ R 2×2 , Ȳj = 1 k k ℓ=1 [Y j ] ℓ ∈ R 2 , that is ĀY j = 1 k ⊗ Ȳj ∈ R k×2 is a degenerated configuration, and Ȳ = 1 J J j=1 Ȳj ∈ R 2 .
This two steps procedure is equivalent to the optimization problem (1.6) as we have the following decompositions M λ (a, α, b) = M λ 0 (a, α) + M (a, α, b) and Θ 0 = Θ a,α 0 × Θ b 0 , implying following result (see the Appendix for a detailed proof)

Lemma 4.1. (â λ , αλ , bλ ) ∈ argmin (a,α,b)∈Θ 0 M λ (a, α, b) ⇐⇒        (â λ , αλ ) ∈ argmin (a,α)∈Θ a,α 0 M λ 0 (a, α) bλ = argmin b∈Θ b 0 M (â λ , αλ , b).

Consistency results

In what follows, C, C 0 , C 1 denote positive constants whose value may change from line to line. The notation C(•) specifies the dependency of C on some quantities.

Consistent estimation of the deformation parameters

Rotation and scaling. Recall that the rotation and scaling parameters are estimated separately on the smoothed and centered observations. We have the following result, Theorem 5.1. Consider model (1.1) and suppose that Assumptions 1 and 2 hold and that Assumption 3 is verified with max{A, A} < 0.1. Consider (â λ , αλ ) the estimators defined in equation

(4.3). If λ = k 1 2s+1
then for all x > 0 we have

P 1 J (â λ , αλ ) -(a * , α * ) 2 R 2J ≥ C 1 (L, s, A, A, f 0 )A 1 (k, J, x) + C 2 (A, A)A 2 (J, x) ≤ 4e -x , (5.1)
where

C 1 (L, s, A, A, f 0 ), C 2 (A, A) > 0 are positive constants independent of k and J, A 1 (k, J, x) = F k -2s 2s+1 +F V 1 (k, J, x) with V 1 (k, J, x) = 3γ max (k)k -2s 2s+1 1 + 2 x Jk + 2 x Jk and A 2 (J, x) = x 3J + 2x J 2
, where F (u) = u + √ u, for all u ≥ 0.

Remark that a direct consequence of Theorem 5.1 is the consistency of (â λ , αλ ) to (a * , α * ) when min{k, J} → ∞. Indeed, we have lim u→0 F (u) = 0 and under Assumption 2 and for any fixed J ≥ 1 and x > 0, the term V 1 (k, J, x) tends to zero as k goes to infinity. Hence for any x > 0 and J ≥ 1, we have lim

k→∞ A 1 (k, J, x) = 0.
Under the same hypothesis as in Theorem 5.1 but without the bounds on A and A, Proposition A.1 then ensures the convergence of (â λ , αλ ) to (a * Θ 0 , α * Θ 0 ) as J remains fixed and k → ∞. Now, for all x > 0, we have lim

J→∞ A 2 (J, x) = 0.
Thus a double asymptotic min{k, J} → ∞ ensures that 1 J (â λ , αλ ) -(a * , α * ) 2 R 2J tends to 0 in probability.

Translation parameters. We have the following result, Theorem 5.2. Consider the hypothesis and notations of Theorem 5.1 and the estimator bλ given by formula (4.4). Then we have for all x > 0,

P 1 J bλ -b * 2 R 2J ≥ C 3 (L, s, A, A, B, f )A 3 (k, J, x) + C 4 (A, A, B)A 2 (J, x) ≤ 9e -x , (5.2) where C 3 (L, s, A, A, B, f ), C 4 (A, A, B) > 0 are positive constants independent of k and J, A 3 (k, J, x) = F (k -2s 2s+1 ) + F (V 1 (k, J, x)) + V 2 (k, J, x) with V 2 (k, J, x) = 4 γmax(k) k 1 + 2 x J + x J .
Similar comments to those made after Theorem 5.1 are still valid here. For any J ∈ N and x > 0, we have lim k→∞ A 3 (k, J, x) = 0, since V 2 (k, J, x) tends to 0 as k goes to infinity by Assumption 2. Under the same hypothesis as in Theorem 5.2 but without the bounds on A and A, Proposition A.3 ensures the convergence in probability of b λ to b * Θ 0 with only an asymptotic in k. In the double asymptotic setting min{k, J} → ∞, Theorem 5.2 ensure the consistency of bλ to the true value b * of the translation parameters.

Consistent estimation of the mean shape

Theorem 5.3. Consider model (1.1) and suppose that Assumptions 1 and 2 hold and that Assumption 3 is verified with max{A, A} < 0.1. Consider the estimator f λ defined in (1.8) and let λ = k 1 2s+1 . Then, we have for all x > 0,

P 1 k f λ -f 2 R k×2 ≥ C(L, s, A, A, B, f ) A 1 (k, J, x) + A 3 (k, J, x) + A 2 (J, x) ≤ 14e -x , where C(L, s, A, A, B, f ) > 0 is a constant independent of k and J, A 1 (k, J, x) = F k -2s 2s+1 + F V 1 (k, J, x) with V 1 (k, J, x) = 3γ max (k)k -2s 2s+1 1 + 2 x Jk + 2 x Jk , A 3 (k, J, x) = V 2 (k, J, x) + F (k -2s 2s+1 ) + F (V 1 (k, J, x)) with V 2 (k, J, x) = 4 γmax(k) k 1 + 2 x J + x J and A 2 (J, x) = 2x J + x 3J 2 .
The terms A 1 (k, J, x), A 3 (k, J, x) and A 2 (J, x) that appear in the statement of Theorem 5.3 are the same to those appearing in Theorems 5.1 and 5.2. As a consequence, we have 1 k f λf 2 R k×2 → 0 in probability when min{k, J} → ∞ and Theorem 5.3 gives rates of convergences of f λ to f thanks to a concentration inequality.

Theorem 5.3 is similar to Theorems 5.1 and 5.2 and there is no need to assume an extra bound on A and A to ensure the convergence in probability of f λ when J is fixed and k → ∞, see the discussion following Theorems 5.1 and 5.2. Let

f Θ 0 = e ā * (f + 1 k ⊗ b * (e a * R α * ) -1 )R ᾱ * (5.3) where ā * = 1 J J j=1 a * j ∈ R, ᾱ * = 1 J J j=1 α * j ∈ R, b * = 1 J J j=1 b * j ∈ R 2 and e a * R α * = 1 J J j=1 e a * j R α * j
is an invertible 2 × 2 matrix, see also formula (3.3). A slight modification of the proof of Theorem 5.3 gives the following inequality,

P 1 k f λ -f Θ 0 2 ≥ C(L, s, A, A, B, f ) A 1 (k, J, x) + A 3 (k, J, x) ≤ 14e -x .
This yields to statements (1.9) in Theorem 1.1.

6 Numerical experiments

Description of the data

We make here some numerical simulations to show the effect of the dimension k and the number J of observations on the estimation of the deformation parameters and the mean pattern with data generated by model (1.1). Different types of noise are considered. For all t ∈ [0, 1], let f (t) = (10 sin 2 (πt) + cos(10π) + 20, 2 sin(6πt) -11 sin 2 (πt) + 12 exp(-25(t -0.4) 2 ) + 1).

Θ = [-1 4 , 1 4 ] × [-1 2 , 1 2 ] × [-1, 1] 2 .
The law of the deformation parameters is supposed to be unknown a priori and the minimization is performed on the constraint set

Θ 0 = (a, α, b) ∈ [-1, 1] J × [-1, 1] J × [-5, 5] 2J , J j=1 a j = 0,
(k) 1 , ζ (k) 2 ) ′ ∈ R 2k is a centered Gaussian vector of variance Σ 1 = 4Id 2k .
We have γ max (k) = 4 and this correspond to an isotropic Gaussian noise as in [START_REF] Kent | [END_REF], Le, 1998], see Figure 4.

Weakly correlated noise : the random variable ζ is a centered Gaussian vector of variance Σ 2 with

Σ 1 2 2 = Id 2 ⊗ 1 2 exp -|ℓ-ℓ ′ | 100 2 k ℓ,ℓ ′ =1
Hence, Σ 2 is a Toeplitz matrix and it follows from classical matrix theory, see e.g. [START_REF] Horn | [END_REF], that γ max (k) is bounded (here γ max ≤ 80). See Figure 5.

Highly correlated noise : the random variable ζ is a centered Gaussian vector of variance

Σ 3 = Id 2 ⊗ P Diag ℓ 2 2k , ℓ = 1, . . . , k P ′
where P is an arbitrary matrix in SO(k). Hence, in this case γ max (k) = k 4 and the level of noise increase with k, see Figure 6.

Description of the procedure

The estimation procedure follows the guidelines described in Section 4. We are testing the effect of the number J of observations and the number k of landmarks on the estimation of the parameters of interest of model (1.1). All the simulations are performed with J = 10, 100, 500 and k = 20, 50, 100, 1000, 3000 Moreover, estimations are done without and with the pre-smoothing step. In the former case we have λ(k) = k 2 , that is, there is no reduction of the dimension. In the latter case, the smoothing parameter λ(k) is fixed manually to ensure a proper reconstruction of the mean pattern f . Note that we need λ > 7 to get correct results and we took λ 20 = λ 50 = λ 100 = 7, λ 1000 = 11 and λ 3000 = 25.

We use a quasi-Newton trust-region based algorithm to solve the optimization problems (4.3) and (4.4). The formula for the gradient is given in (B.1). All the computations are performed with Matlab.

Results: estimation of the mean pattern

For of the 30 of model (1.1) with the possible values of k and J, we compute the quantities

1 k f λ -f * Θ 0 2
R k×2 where f λ corresponds to the smoothed Procrustes mean of the observations defined in (1.8) and,

1 k f -f * Θ 0 2
R k×2 where f is the (non smoothed) Procrustes mean of the data. Recall that f Θ 0 is defined by formula (5.3).

Boxplots of the results are given in Figures 7,8 and 9 for the different kinds of error terms described in Section 6.1. In the figures, the abscissa represents the different values of the number k of landmarks and boxplots in red correspond to = 10 observations, in green to J = 100 observations and in blue to J = 500 observations.

The estimation of the mean pattern with the white noise error term is given by Figure 7. In Figure 7a, for a fixed k, the non-smoothed version

1 k f -f * 2 R k×2 decreases when J increases. Moreover, the values of 1 k f -f * 2 R k×2
remain stable when J remains fixed and k increases. Recall that this framework corresponds to the isotropic Gaussian noise described in [START_REF] Kent | [END_REF]. The simulations seem to confirm their conclusions and show that in this framework the dimension k is not preponderant. In Figure 7b, the smoothed version

1 k f λ -f 2
R k×2 decreases when J and k increase. The main difference with the non-smoothed estimation is the convergence to 0 of 1 k f λf 2 R k×2 when J remains fixed and k increases.

In Figure 8, the results of the estimation of the mean pattern are plotted for the weakly correlated noise term. Figure 8a shows us a similar behavior of the non-smoothed Procrustes mean but with nondecreasing values of

1 k f λ -f * 2
R k×2 when k increases and J remains fixed. In Figure 8b, the smoothed Procrustes mean converges as k goes to infinity and the bigger J is the faster the convergence is.

The results of the estimations of the mean pattern with the highly correlated noise are presented Figure 9. The results that appear in Figure 9a are quite different compared to those presented Figures 7a and8a. The estimation seems to be worst when k increases and J remains fixed. The reason is that the level of noise, measured by γ max (k) is increasing with k. The smoothing step is efficient and the estimations presented Figure 9b have a similar behavior to those given in Figures 7b and8b.

A Proofs

A.1 Proof of Lemma 4.1

Using the decomposition (4.2), we obtain the following identity

M λ (a, α, b) = M λ 0 (a, α) + M (a, α, b). (A.1)
Note that we have used the fact that the subspaces V λ 0 and V are orthogonal and that Thence, thanks to the decomposition (A.1) and the fact that M (a, α, b(a, α)) = 0, we have argmin

Ā1 k ⊗ b j = 1 k ⊗ b j , for all b j ∈ R 2 . Let us also introduce the notation Ȳj = 1 k k ℓ=1 [Y j ] ℓ ∈ R 2 , that is ĀY j = 1 k ⊗ Ȳj ∈ R k×2 is a degenerated configuration, and Ȳ = 1 J J j=1 Ȳj ∈ R 2 . For a fixed (a, α) ∈ R J × [-π, π[ J , the functional b -→ M (a, α, b) vanishes if and only if there exists a b 0 ∈ R 2 such that e -a j ( ĀY j -1 k ⊗ b j )R -α j = 1 k ⊗ b 0 for all j = 1, . . . , J. Therefore, for this fixed (a, α), there is a unique point b = b(a, α) := ( Ȳ1 -e a 1 Ȳ(e a R α ) -1 R α 1 , . . . , ȲJ -e a J Ȳ(e a R α ) -1 R α J ) with
(a,α,b)∈Θ 0 M λ (a, α, b) = argmin (a,α)∈Θ a,α 0 M λ 0 (a, α) , b argmin (a,α)∈Θ a,α 0 M λ 0 (a, α) ,
and the claim is proved.

A.2 Proof of Theorem 5.1

For all (â λ , αλ ) ∈ Θ a,α 0 and (a * , α * ) ∈ [-A, A] × [-A, A],
we have the following inequality

1 J (â λ , αλ ) -(a * , α * ) 2 R 2J ≤ 2 J (â λ , αλ ) -(a * Θ 0 , α * Θ 0 ) 2 R 2J + 2 J (a * Θ 0 , α * Θ 0 ) -(a * , α * ) 2 R 2J .
The proof of Theorem 5.1 is a direct consequence of Proposition A.1 and Lemma A.2 below which control the convergence in probability of the two terms in the right hand side in the preceding inequality.

Proposition A.1. Consider model (1.1) and suppose that Assumptions 1 and 2 hold and that Assumption 3 is verified with A,

A < 0.1. If λ = λ(k) = k 1 2s+1 then there exists a constant C(L, s, A, A, f 0 ) such that for all x > 0 P 1 J (â λ , αλ ) -(a * Θ 0 , α * Θ 0 ) 2 R 2J ≥ C(L, s, A, A, f 0 ) F k -2s 2s+1 + F V 1 (k, J, x) ≤ 2e -x , where V 1 (k, J, x) = 3γ max (k)k -2s 2s+1 1 + 2 x Jk + 2 x Jk and F : R + -→ R, with F (u) = u + √ u.
The proof of Proposition A.1 is postponed to Section A.5. The following lemma is a direct consequence of Bernstein's inequality for bounded random variables, see e.g. Proposition 2.9 in [Massart, 2007].

Lemma A.2. Suppose that Assumption 3 holds and that the random variables (a * j , α * j ), j = 1, . . . , J have zero expectation in [-A, A] × [-A, A]. Then, for any x > 0, we have

P 1 J (a * Θ 0 , α * Θ 0 ) -(a * , α * ) 2 R 2J ≥ C(A, A) 2x J + x 3J 2 ≤ 4e -x , where C(A, A) = 4 max{A 2 , A 2 }.
A.3 Proof of Theorem 5.2

The proof of Theorem 5.2 follows the same guideline as the proof of Theorem 5.1. Consider the inequality

1 J bλ -b * 2 R 2J ≤ 2 J bλ -b * Θ 0 2 R 2J + 2 J b * Θ 0 -b * 2 R 2J .
Theorem 5.2 is now a direct consequence of Proposition A.3 and Lemma A.4.

Proposition A.3. Under the hypothesis of Proposition A.1, there exists a constant C(L, s, A, A, B, f ) > 0 such that for all x > 0,

P 2 J bλ -b * Θ 0 2 R 2J ≥ C(L, s, A, A, B, f ) F (k --2s 2s+1 ) + F (V 1 (k, J, x)) + V 2 (k, J, x) ≤ 5e -x , where V 1 (k, J, x) = 3γ max (k)k -2s 2s+1 1 + 2 x Jk + 2 x Jk , V 2 (k, J, x) = 4 k γ max (k) 1 + 2 x J + x J and F (u) = u + √ u, u ≥ 0.

The proof of Proposition A.3 is postponed to

A < π 4 and that the random variables (a * j , α * j , b * j ) j = 1, . . . , J have zero expectation in Θ ⊂ R 4J . For any x > 0 , we have

P 2 J b * Θ 0 -b * 2 R 2J ≥ C(A, A, B) 2x J + x 3J 2 ≤ 4e -x , where C(A, A, B) = 8B 2 e 4A (cos A -sin A) -2 .
Proof. This result is a consequence of the Bernstein's inequality for bounded random variable. We have

b * Θ 0 -b * R 2J = J j=1 e a * j b * (e a R α ) -1 R α * j 2 R 2 ≤ e A √ J b * (e a * R α * ) -1 R 2 ,
where e a * R α * = 1 J J j=1 e a * j R α * j is an invertible 2 × 2 matrix whose smallest eigenvalue is greater than e -A (cos A -sin A) > 0 as A < π 4 . To see this, remark that the eigenvalues of e a * R α * are 1 J J j=1 e a * j (cos α * j ± i sin α * j ) and we have 1

J J j=1 e a * j (cos α * j ± i sin α * j ) ≥ e -A (cos A -sin A) > 0. We now have 1 √ J b * Θ 0 -b * R 2J ≤ C(A, A) b * R 2 ,
where C(A, A) = e 2A (cos A -sin A) -1 . Finally, for all u > 0 we have P( 1

√ J b * Θ 0 -b * R 2J ≥ u) ≤ P(C(A, A) b * R 2 ≥ u)
and a Bernstein type inequality (see e.g. Proposition 2.9 in [Massart, 2007])

gives us P b * R 2 ≥ 2B 2x J + x 3J ≤ 4e -x which yields P 1 J b * Θ 0 -b * 2 R 2J ≥ C(A, A, B) 2x J + x 3J 2 ≤ 4e -x ,
where C(A, A, B) = 4B 2 e 4A (cos A -sin A) -2 .

A.4 Proof of Theorem 5.3

Recall the notations introduced Section 2.1: for ĝλ j = (â λ j , αλ j , bλ j ) we have ĝλ j .f = e âλ j f R αλ j + bλ j . Then, we have

1 k f λ -f 2 R k×2 = 1 k 1 J J j=1 (ĝ λ j ) -1 .(A λ Y j ) -f 2 R k×2 ≤ 2 kJ J j=1 A λ (ĝ λ j ) -1 .Y j -f 2 R k×2 + 2 k A λ f -f 2 R k×2 ≤ 4 kJ J j=1 A λ (g λ j ) -1 .g * j .f -f 2 R k×2 + e a * j -â λ j A λ ζ j R α * j -αλ j 2 R k×2 V + 2 k A λ f -f 2 R k×2 B .
The rest of the proof is devoted to control the terms B and V. The term B is controlled by the bias of the low pass filter. According to Lemma B.1, if Assumption 1 holds and by choosing the optimal frequency cutoff λ = λ(k) = k 1 2s+1 , there exists a constant C(L, s) > 0 such that

B = 2 k A λ f -f 2 R k×2 ≤ C(L, s)k -2s 2s+1 . (A.2)
The term V contains two expressions. To bound Lemma B.4. More precisely, we all j = 1, . . . , J

1 k A λ (g λ j ) -1 .g * j .f -f 2 R k×2 ≤ 1 k (g λ j ) -1 .g * j .f -f 2 R k×2 ≤ C(A, f ) (a * j -âλ j , α * j -αλ j , e âλ j (b * j -bλ j )R -αλ j ) 2 R 4 ≤ C(A, f ) (a * j -âλ j , α * j -αλ j ) 2 R 2 + C(A, f ) b * j -bλ j 2 R 2
We can now use Theorem 5.1 and 5.2 to derive the following concentration inequality,

P 4 kJ J j=1 A λ (g λ j ) -1 .g * j .f -f 2 R k×2 ≥ C(L, s, A, A, B, f ) A 1 (k, J, x) + A 3 (k, J, x) + A 2 (J, x) ≤ 13e -x (A.3)
where C(L, s, A, A, B, f ) > 0 is a constant independent of k and J and A 1 , A 2 , A 3 are defined in the statement of Theorem 5.1 and 5.2. The second term contained in V is treated by equation (A.10) below. Hence, formulas (A.3) and (A.10) yield

P V ≥ C(L, s, A, A, B, f ) A 1 (k, J, x) + A 3 (k, J, x) + A 2 (J, x) ≤ 14e -x , (A.4)
for some constant C(L, s, A, A, B, f ) > 0. Putting together equations (A.2) and (A.4) gives

P 1 k f λ -f 2 R k×2 ≥ C(L, s, A, A, B, f ) A 1 (k, J, x) + A 3 (k, J, x) + A 2 (J, x) ≤ 14e -x .
for some constant C(L, s, A, A, B, f ) > 0. The proof of Theorem 5.3 is completed.

A.5 Proof of Proposition A.1

The mean pattern f can be decomposed as

f = f + f 0 ∈ 1 k×2 ⊕ 1 ⊥ k×2 .
Then, f 0 is the centered version of f and we can consider the criterion,

D 0 (a, α) = 1 Jk J j=1 e a * j -a j f 0 R α * j -α j - 1 J J j ′ =1 e a * j ′ -a j ′ f 0 R α * j ′ -α j ′ 2 R k×2
.

We now have,

(â λ , αλ ) = argmin (a,α)∈Θ a,α 0 M λ 0 (a, α) and (a * Θ 0 , α * Θ 0 ) = argmin (a,α)∈Θ a,α 0 D 0 (a, α).
Then, the convergence of (â λ , αλ ) to (a

* Θ 0 , α * Θ 0 ) is guaranteed if (a * Θ 0 , α * Θ 0
) is uniquely defined and if there is a uniform convergence in probability of M λ 0 to D 0 , see e.g. [van der Vaart, 1998]. This is the aim of Lemmas A.5 and A.6 below.

Lemma A.5. Let f be a non-degenerated configuration in R k×2 , i.e. f / ∈ 1 k×2 . Then, the argmin of D 0 on Θ a,α 0 is unique and denoted by (a

* Θ 0 , α * Θ 0 ) = (a j -ā * , α j -ᾱ * ), where ā * = 1 J J j=1 a * j , ᾱ * = 1 J J j=1 α * j .
Proof. As f is a non-degenerated configuration, we have f 0 = 0. Thus, the stabilizer of f is reduced to the identity, see Section 2.1. Then D 0 (a, α) = 0 if and only if there exists

(a 0 , α 0 ) ∈ R 2 such that (a, α) = (a * , α * ) * (a 0 , α 0 ) = (a * 1 +a 0 , α * 1 +α 0 , . . . , a J +a 0 , α * J +α 0 ). By choosing (a 0 , α 0 ) = (-ā * , -ᾱ * ) we have J j=1 (a * j , α * j ) * (a 0 , α 0 ) = 0. That is (a * Θ 0 , α * Θ 0 ) = (a * , α * ) * (a 0 , α 0 ) ∈ Θ a,α 0 .
We now show the uniform convergence in probability, and let F : R -→ R, with F (u) = u + √ u. For any x > 0 we have

P sup (a,α)∈Θ a,α 0 M λ 0 (a, α) -D 0 (a, α) ≥ C(L, s, A, f 0 ) F k -2s 2s+1 + F V (k, J, x) ≤ 2e -x
where C(L, s, A,

f 0 ) = e 2A max 2 √ k f 0 R k×2 , √ 2e A √ k f 0 R k×2 , 2 and V (k, J, x) = 3γ max (k)k -2s 2s+1 1 + F 2x Jk .
Proof. Let us write the following decomposition,

M λ 0 (a, α) = 1 Jk J j=1 e a * j -a j f 0 R α * j -α j - 1 J J j ′ =1 e a * j ′ -a j ′ f 0 R α * j ′ -α j ′ 2 R k×2 (A.5) + 1 Jk J j=1 e a * j -a j (A λ 0 f -f 0 )R α * j -α j - 1 J J j ′ =1 e a * j ′ -a j ′ (A λ 0 f -f 0 )R α * j ′ -α j ′ 2 R k×2 (A.6) + 2 Jk J j=1 e a * j -a j f 0 R α * j -α j - 1 J J j ′ =1 e a * j ′ -a j ′ f 0 R α * j ′ -α j ′ , e a * j -a j (A λ 0 f -f 0 )R α * j -α j - 1 J J j ′ =1 e a * j ′ -a j ′ (A λ 0 f -f 0 )R α * j ′ -α j ′ R k×2 (A.7) + 1 Jk J j=1 e a * j -a j A λ 0 ζ j R α * j -α j - 1 J J j ′ =1 e a * j ′ -a j ′ A λ 0 ζ j ′ R α * j ′ -α j ′ 2 R k×2 (A.8) + 2 Jk J j=1 e a * j -a j A λ 0 f R α * j -α j - 1 J J j ′ =1 e a * j ′ -a j ′ A λ 0 f R α * j ′ -α j ′ , e a * j -a j A λ 0 ζ j R α * j -α j - 1 J J j ′ =1 e a * j ′ -a j ′ A λ 0 ζ j ′ R α * j ′ -α j ′ R k×2 (A.9)
Then, criterion M λ 0 is viewed as a perturbed version of the criterion D 0 (a, α) = (A.5),

M λ 0 (a, α) = D 0 (a, α) + B λ 0 (a, α) + V λ 0 (a, α)
where the bias term is B λ 0 (a, α) = (A.6) + (A.7) and the variance term is V λ 0 (a, α) = (A.8) + (A.9).

The bias term. We have for any (a, α) ∈ Θ a,α 0 ,

(A.6) = 1 Jk J j=1 e a * j -a j (A λ 0 f -f 0 )R α * j -α j 2 R k×2 - 1 J J j ′ =1 e a * j ′ -a j ′ (A λ 0 f -f 0 )R α * j ′ -α j ′ 2 R k×2 ≤ e 2A 1 k A λ 0 f -f 0 2 R k×2
.

A double application of Cauchy-Schwarz inequality implies that, e * j -a j α * j -α j -

1 J J j ′ =1 e a * j ′ -a j ′ f 0 R α * j ′ -α j ′ R k×2 e a * j -a j (A λ 0 f -f 0 )R α * j -α j - 1 J J j ′ =1 e a * j ′ -a j ′ (A λ 0 f -f 0 )R α * j ′ -α j ′ R k×2 ≤ 2 k 1 J J j=1 e 2A f 0 R α * j -α j 2 R k×2 1 2 1 J J j=1 e 2A (A λ 0 f -f 0 )R α * j -α j 2 R k×2 1 2 = 2e 2A 1 k f 0 R k×2 A λ 0 f -f 0 R k×2 .
Finally, by using Lemma B.1 we have,

sup (a,α)∈Θ a,α |B λ 0 (a, α)| ≤ C 1 (L, s, A, f ) k -2s 2s+1 + k -s 2s+1 .
where

C 1 (L, s, A, f ) = C(L, s)e 2A max 2 1 √ k f 0 R k×2 , 1 .
The variance term. First, the term (A.9) is by the Cauchy-Schwarz inequality controlled by 2e 2A 1 .8). The term (A.8) is bounded from above by 1

√ k f 0 R k×2 (A
Jk J j=1 e 2A A λ 0 ζ j 2 R k×2 .
To derive an upper bound in probability, note that we have the following equality in law,

J j=1 A λ 0 ζ j 2 R k×2 = ξ ′ B ξ, with B = Id J ⊗ Σ 1 2 Id 2J ⊗ (A λ 0 ) ′ A λ 0 Id J ⊗ Σ 1 2
∈ R 2Jk×2Jk and ξ = (ξ 1 , . . . , ξ 2Jk ) ′ is a centered Gaussian vector of variance Id 2Jk . We have tr(Id J ⊗ Σ) ≤ 2Jkγ max (k) and tr((A λ 0 ) ′ A λ 0 ) = 2λ+1 k . Using a classical concentration inequality for quadratic form of multivariate Gaussian random variables, see e.g. [Laurent & Massart, 2000] Lemma 1, we have for all x > 0, P ξ ′ Bξ ≥ 2Jkγ max (k)

2λ+1 k + 2γ max (k) 2λ+1 k √ x2Jk + 4xγ max (k) 2λ+1
k ≤ e -x , which yields together with formula (A.2),

P 1 Jk J j=1 e 2A A λ 0 ζ j 2 R k×2 ≥ 6e 2A γ max (k)k -2s 2s+1 1 + 2 x Jk + 2 x Jk ≤ e -x .
(A.10)

Hence we have,

P sup (a,α)∈Θ a,α |V λ 0 (a, α) + B λ 0 (a, α)| ≥ C(L, s, A, f 0 ) k -2s 2s+1 + k -s 2s+1 + V 1 (k, J, x) + V 1 (k, J, x) ≤ 2e -x , where C(L, s, A, f 0 ) = e 2A max 2 √ k f 0 R k×2 , 1 √ k f 0 R k×2 √ 2e A , 2 and V 1 (k, J, x) = 3 γmax(k) k 2s 2s+1 1 + 2 x Jk + 2 x Jk .
For a fixed J ∈ N, the convergence of the M-estimator (â λ , αλ ) to (a * Θ 0 , α * Θ 0 ) when k → ∞ is guaranteed by Lemma A.5 and A.6, see e.g. [van der Vaart, 1998]. Nevertheless, we are able to give a rate of convergence and non-asymptotic bounds in k and J by using the classical inequality,

|D 0 (â λ , αλ ) -D 0 (a * Θ 0 , α * Θ 0 )| ≤ 2 sup (a,α)∈Θ a,α 0 |D 0 (a, α) -M λ 0 (a, α)|.
This, together with Lemma A.7 below will prove Proposition A.1.

Lemma A.7. Assume that A, A < 0.1. There such that for all (a, α) ∈ Θ a,α we have

D 0 (a, α) -D 0 (a * Θ 0 , α * Θ 0 ) ≥ C(A, A, f 0 ) 1 J (a -a * Θ 0 , α -α * Θ 0 ) 2 R 2J , where C(A, A, f 0 ) = C(A, A) 1 k f 0 2 R k×2 . Proof. By definition, given a (a * , α * ) ∈ [-A, A] J × [-A, A] J , the point (a * Θ 0 , α * Θ 0 ) is the unique minimum of D 0 on Θ a,α 0 = [-A, A] J ∩ 1 ⊥ J × [-A, A] J ∩ 1 ⊥ J .
Then, for all (a, α) ∈ Θ a,α 0 , there exists a c = c(a, α) ∈ Θ a,α 0 such that the Taylor expansion of D 0 at (a * Θ 0 , α * Θ 0 ) can be written,

D 0 (a, α) -D 0 (a * Θ 0 , α * Θ 0 ) = 1 2 (a -a * Θ 0 , α -α * Θ 0 ) ′ [∇ 2 D 0 (a * Θ 0 , α * Θ 0 )](a -a * Θ 0 , α -α * Θ 0 ) + 1 6 [∇ 3 D 0 (c)](a -a * Θ 0 , α -α * Θ 0 ).
Let δ = max{A, A}. This, together with Lemma B.2 and B.3 imply that

D 0 (a, α) -D 0 (a * Θ 0 , α * Θ 0 ) ≥ 1 2 (a -a * Θ 0 , α -α * Θ 0 ) ′ [∇ 2 D 0 (a * Θ 0 , α * Θ 0 )](a -a * Θ 0 , α -α * Θ 0 ) -δ 40 6 e 2A f 0 2 R k×2 1 Jk (a -a * Θ 0 , α -α * Θ 0 ) 2 R 2J ≥ f 0 2 R k×2 1 Jk (a -a * Θ 0 , α -α * Θ 0 ) 2 R 2J e -2A -δ 40 6 e 2A .
Hence, one can choose δ > 0 sufficiently small such that e -2A -δ 40 6 e 2A is strictly positive for all J and k. For example, we have e -2δ -δ 40 6 e 2δ > 0, if δ < 0.1. Then, using such a δ it follows that for all (a, α) ∈ Θ a,α 0 ,

D 0 (a, α) -D 0 (a * Θ 0 , α * Θ 0 ) ≥ C(A, A) 1 k f 0 2 R k×2 1 J (a -a * Θ 0 , α -α * Θ 0 ) 2 R 2J .
The proof of Proposition A.1 is almost done. Remark that Lemma A.7 ensures that for all u ≥ 0 we have P

1 J (â λ , αλ ) -(a * Θ 0 , α * Θ 0 ) 2 R 2J ≥ u ≤ P 2 C(A,A,f 0 ) sup (a,α)∈Θ a,α 0
M λ 0 (a, α) -D 0 (a, α) ≥ u . Lemma A.6 ensures that there is a constant C(L, s, A, A, f 0 ) such that for all x > 0,

P 1 J (â λ , αλ ) -(a * Θ 0 , α * Θ 0 ) 2 R 2J ≥ C(L, s, A, A, f 0 )(F (k -2s 2s+1 ) + F (V 1 (k, J, x)) ≤ 2e -x .
where

V 1 (k, J, x) = γ max (k)k -2s 2s+1 1 + 2 x Jk + 2 x Jk and F : R -→ R, with F (u) = u + √ u.

A.6 Proof of Proposition A.3

First of all remark that, thanks to formulas (3.3) and (4.

* Θ 0 = (b * 1 -e a * 1 b * (e a * R α * ) -1 R α * 1 , . . . , b * J -e a * J b * (e a * R α * ) -1 R α * J )
,

where e a * R α * = 1 J J j=1 e a * j R α * j ∈ R 2×2 , b * = 1 J J j=1 b * j ∈ R 2 and bλ = ( Ȳ1 -e âλ 1 Ȳ(e âλ R αλ ) -1 R αλ 1 , . . . , ȲJ -e âλ J Ȳ(e âλ R αλ ) -1 R αλ J )
where

e âλ R αλ = 1 J J j=1 e âλ j R αλ j ∈ R 2×2 , Ȳj = 1 k k ℓ=1 [Y j ] ℓ = e a * j f R α * j + b * j + e a * j ζj R α * j ∈ R 2 where f = 1 k k ℓ=1 f ℓ ∈ R 2 , ζj = 1 k k ℓ=1 [ζ j ] ℓ ∈ R 2 and Ȳ = 1 J J j=1 Ȳj = f (e a * R α * ) + b * + 1 J J j=1 e a * j R . Thence, we have * Θ 0 -bλ 2 R 2J ≤ 1 J J j=1 f (e a * j R α * j -e âλ j (e a * R α * )(e âλ R αλ ) -1 R αλ j ) 2 R 2 (A.11) + b * (e a * j (e a * R α * ) -1 R α * j -e âλ j (e âλ R αλ ) -1 R αλ j ) 2 R 2 (A.12) + e a * j ζj R α * j 2 R 2 + e âλ j (e a * ζR α * )(e âλ R αλ ) -1 R αλ j 2 R 2 (A.13)
where (e a * ζR α

* ) = 1 J J j=1 e a * j ζj R α * j ∈ R 2 .
The rest of the proof is devoted to the control of the terms (A.11), (A.12) and (A.13).

In this section, we denotes by • op the operator norm of a 2 × 2 matrix, i.e. A op = |γ max (A)| where γ max (A) denotes the largest eigenvalue of a matrix A ∈ R 2×2 . Note by the way that the eigenvalues of the matrix 1 J J j=1 e a j R α j ∈ R 2×2 are 1 J J j=1 e a j (cos(α j ) ± i sin(α j )) for any J and (a 1 , . . . , a J , α 1 , . . . , α J ) ∈ [-A, A] J × [-A, A] J . It yields e -a j R -α j op ≤ e A and (e a R α ) -1 op ≤ e A (cos(A) -sin(A)) -1 which is a positive real number since by hypothesis we have A < π 4 . We are now able to derive an upper bound for (A.11),

(A.11) ≤ f 2 R 2 1 J J j=1 e a * j R α * j -e âλ j (e a * R α * )(e âλ R αλ ) -1 R αλ 2 op ≤ e 2A f 2 R 2 1 J J j=1 e a * j R α * j (e -â λ j R -αλ j ) -(e a * R α * )(e âλ R αλ ) -1 2 op ≤ 2e 2A f 2 R 2 1 J J j=1 e -â λ j R -αλ j 2 op e a * j R α * j -e âλ j R αλ j 2 op + 2e 2A f 2 R 2 (e âλ R αλ ) -1 2 op (e a * R α * ) -(e âλ R αλ ) 2 op
It is now easy to show that there exists a constant C(A, A, f ) > 0 independent of k and J such that

(A.11) ≤ C(A, A, f ) 1 J J j=1 (a * j -âλ j , α * j -αλ j ) 2 R 2 . (A.14)
The term (A.12) is very similar to the term (A.11) and we have

(A.12) ≤ C(A, A, B) 1 J J j=1 (a * j -âλ j , α * j -αλ j ) 2 R 2 , (A.15)
for some constant C(A, A, B) > 0 independent of k and J. By using formula (A.14), (A.15) and Proposition A.1 together, there exists a constant C(L, s, A, A, B, f ) independent of k and J such that for all x > 0,

P (A.11) + (A.12) ≥ C(L, s, A, A, B, f ) F (k -2s 2s+1 ) + F (V 1 (k, J, x)) ≤ 4e -x , (A.16) where V 1 (k, J, x) is defined in statement of proposition A.1 and F (u) = u + √ u for u ≥ 0.
The term (A.13) can be bounded in probability as follows. First of all, remark that there exists a constant C(A, A) > 0 such that

(A.13) ≤ C(A, A) 1 J J j=1 ζj 2 R 2 .
Then, for all j = 1, .

. ζj = k ℓ=1 ([ζ (1) ℓ , [ζ (2) j ] ℓ ) ∈ R 2 can be written ( ζ1 j , ζ2 j ) ′ = 1 k 1 ′ k 0 ′ k 0 ′ k 1 ′ k ([ζ (1) j ] 1 , . . . , [ζ (1) j ] k , [ζ (2) j ] 1 , . . . , [ζ (2) j ] k ) ′
where 0 k is a column vector of k zeros. The random vector ( ζ1 j , ζ2 j ) ′ is a two dimensional centered Gaussian vector of variance

V = 1 k 2 1 ′ k 0 ′ k 0 ′ k 1 ′ k Σ 1 k 0 k 0 k 1 k . The 2 × 2 matrix V is of trace less or equal to 2 k γ max (k). Hence, the random variable 1 J J j=1 ζj 2 R 2 = 1 J J j=1 ( ζ(1) j , ζ(2) j )( ζ(1) j , ζ (
2) j ) ′ has the same probability distribution as 1 J ξ ′ [Id J ⊗ V]ξ where ξ is a centered Gaussian vector of variance Id 2J . A standard concentration inequality of χ 2 distribution (see e.g. [Laurent & Massart, 2000]

Lemma 1) is P ξ ′ [Id J ⊗ V]ξ ≥ J 4 k γ max (k) + 4 k γ max (k) √ 2Jx + x 4 k γ max (k) ≤ e -x
for any x > 0. It yields that there exists a constant C(A, A) > 0 such that for all x > 0

P (A.13) ≥ C(A, A) γ max (k) k 1 + 2 x J + x J ≤ e -x .
(A.17)

To end the proof remark that for all u ≥ 0 we have P( .11) + (A.12) + (A.13) ≥ u). This together with (A.16) and (A.17) yield that there exists a constant C(L, s, A, A, B, f ) > 0 independent of k and J such that for all x > 0 we have

1 J b * Θ 0 -bλ 2 R 2J ≥ u) ≤ P((A
P 1 J b * Θ 0 -bλ 2 R 2J ≥ C(L, s, A, A, B, f ) F (k --2s 2s+1 ) + F (V 1 (k, J, x)) + V 2 (k, J, x) ≤ 5e -x , where V 2 (k, J, x) = γmax(k) k 1 + 2 x J + x
J and the proof of Proposition A.3 is completed.

B Technical Lemma

Lemma B.1. Assume that Assumption 1 holds, i.e. f ∈ H s (L) with s > 3 2 (see

(1.3)) and f = (f ( ℓ k )) k ℓ=1 ∈ R k×2 . If λ(k) = k 1 2s+1
then there exists a constant C(L, s) such that for all f ∈ H s (L) we have

1 k A λ f -f 2 R k×2 ≤ C(L, s)k -2s 2s+1 ,
where A λ is the projection matrix defined in (4.2).

Proof. Recall the notations introduced Sections 1.1 and 4.1 :

c m (f ) = (c m (f (1) ), c m (f (2) )) ∈ C 2 is the m-th Fourier coefficient of f ∈ L 2 ([0, 1], R 2 ) and c m (f ) = (c (1) m (f ), c (2)
m (f )) ∈ C 2 is the m-th discrete Fourier coefficient of f ∈ R k×2
. Thus, we have by Parseval's equality

1 k A λ f -f 2 R k×2 = 1 k 1 k |m|>λ c m (f )e i2πm ℓ k k ℓ=1 2 R k×2 = 1 k 2 |m|>λ c m (f ) 2 C 2 .
where for all c = (c (1) , c

(2) ) ∈ C 2 , c 2 C 2 = |c (1) | 2 + |c (2) | 2 . It yields 1 k A λ f -f 2 R k×2 ≤ 2 |m|>λ 1 k c m (f ) -c m (f ) 2 C 2 + 2 |m|>λ c m (f ) 2 C 2 .
Firstly, Lemma 1.10 in [Tsybakov, 2009]

ensures that if s > 1 2 then 1 k c m (f (i) ) -c m (f (i) ) ≤ C(L, s)k -s+ 1 2
for any m ∈ N and i = 1, 2. Secondly, equation (1.87) in [Tsybakov, 2009]

ensure that if λ(k) = k 1 2s+1 then |m|>λ |c m (f (i) )| 2 ≤ C(L, s)k -2s 2s+1 . Thence, there is a constant C(L, s) independent of k such that 1 k A λ f -f 2 R k×2 ≤ C(L, s)(k 2-2s + k -2s 2s+1 ). If s ≥ 3 2 then we have 1 2 k×2 ≤ C(L, s)k .
We now give a general expression of the gradient and of the Hessian of the criterion D. Recall that we have

D(a, α, b) = 1 J J j=1 e -aj (e a * j f R α * j + 1 k ⊗ (b * j -b j ))R -αj + 1 J J j ′ =1 e -a j ′ (e a * j ′ f R α * j ′ + 1 k ⊗ (b * j ′ -b j ′ ))R -α j ′ 2 R k×2 where (a, α, b) = (a 1 , . . . , a J , α 1 , . . . , α J , b 1 , . . . , b J ) ∈ R J × [-π, π[ J ×R 2J .
To shorten the formulas below we note g j = (g

(1) j , g (2)
j , g

j , g

j ) = (a j , α j , b j ), that is g

(1) j = a j , g

(2) j = α j , g

(3) j = b

(1) j and g (4)

j = b (2) j . Let f g j = g -1 .g * j .f = e -a j (e a * j f R α * j + 1 k ⊗ (b * j -b j ))R -α j
, and for all j 1 = 1, . . . , J and p 1 = 1, . . . , 4,

∂ g (p 1 ) j 1 D(a, α, b) = 2 Jk ∂ g (p 1 ) j 1 f g j 1 , f g j 1 - 1 J J j ′ =1 f g j ′ R k×2 . (B.1)
The second order derivatives are .g 0 , ∂ g (p 2 )

∂ g (p 2 ) j 2 ∂ g (p 1 ) j 1 D(a, α, b) = - 2 J 2 k ∂ g (p 1 ) j 1 f g j 1 , ∂ g (p 2 ) j 2 f g j 2 R k×2 , if j 1 = j 2 , (B.2) ∂ g (p 2 ) j 1 ∂ g (p 1 ) j 1 D(a, α, b) = 2 Jk ∂ g (p 2 ) j 1 ∂ g (p 1 ) j 1 f g j 1 , f g j 1 - 1 J J j ′ =1 f g j ′ R k×2 + 2 Jk - 2 J 2 k ∂ g (p 1 ) j 1 f g j 1 , ∂ g (p 2 ) j 1 f g j 1 R k×2 . (B.
j 2 f g * j 2 .g 0 R k×2 , if j 1 = j 2 , 2 Jk -2 J 2 k ∂ g (p 1 ) j 1 f g * j 1
.g 0 , ∂ g (p 2 )

j 2 f g * j 1
.g 0 R k×2

.

(B.4)

Lemma B.2. The smallest eigenvalue of ∇ 2 D 0 (a * Θ 0 , α * Θ 0 ) restricted to the subset Θ 0 is greater than e -2A 2 Jk f 0 2 R k×2 . Proof. In this proof, h = (a, α) ∈ R J ×[-π, π[ J and f 0 h j = e a * j -a j f 0 R α * j -α j . We have ∂ a j 1 f 0

h j 1 = -f 0 h j 1 and ∂ α j 1 f 0 h j 1 = f 0 h j 1 R -π
2 for all j 1 = 1, . . . , J. By using Formulas (B.4), the Hessian of D 0 at the point (a * Θ 0 , α Θ 0 ) = (a * -ā * , α * -ᾱ * ) ∈ Θ a,α 0 is given by,

∂ α j 1 ∂ α j 2 D 0 (a * -ā * , α * -ᾱ * ) = ∂ a j 1 ∂ a j 2 D 0 (a * -ā * , α * -ᾱ * ) = -2 J 2 k e ā * f 0 2 R k×2 , if j 1 = j 2 , 2 Jk -2 J 2 k e ā * f 0 2 R k×2
and the second order cross derivatives are 0. The Hessian of D 0 at (a * -ā * , α * -ᾱ * ) can be written,

∇ 2 D 0 (a * -ā * , α * -ᾱ * ) = 2 J 2 k e ā * f 0 2 R k×2 JId J -1 J×J 0 0 JId J -1 J×J , (B.5)
where Id J is the identity in R J and 1 J is the J × J matrix of ones. The eigenvalues of JId J -1 J×J are 0 with eigenvector 1 J and J with eigenspace 1 ⊥ J . It yields that on Θ a,α 0 the smallest eigenvalue of ∇ 2 D 0 (a * -ā * , α * -ᾱ * ) is e 2ā * 2 Jk f 0 2 R k×2 . To finish the proof, remark that ā * ≤ A.

Lemma B.3. Let δ = max{A, A}. For . . a J , α . . . , α J ) ∈ Θ a,α 0 and (a, α) ∈ R 2J we have

∇ 3 D 0 (c) (a, α) ≤ 40δe 2A f 0 2 R k×2 1 Jk (a, α) 2 R 2J .
Proof. In this proof, h = (a, α) ∈ R J × [-π, π[ J , that is h

(1) j = a j is the parameter of scaling and h

(2) j = α j is the parameter of rotation. We have f 0 h j = e a * j -a j f 0 R α * j -α j . Then, from equations (B.2) and (B.3), it follows that for all j 1 , j 2 , j 3 = 1, . . . , J and p 1 , p 2 , p 3 = 1, . . . , 2,

∂ h (p 3 ) j 3 ∂ h (p 2 ) j 2 ∂ h (p 1 ) j 1 D 0 (a, α) = 0,
if j 1 = j 2 and j 2 = j 3 and j 1 = j 3 , ∂ h (p 3 )

j 2 ∂ h (p 2 ) j 1 ∂ h (p 1 ) j 1 D 0 (a, α) = - 2 J 2 k ∂ h (p 2 ) j 1 ∂ h (p 1 ) j 1 f 0 h j 1 , ∂ h (p 3 ) j 2 f 0 h j 2 R k×2 , if j 1 = j 2 , ∂ h (p 3 ) j 1 ∂ h (p 2 ) j 1 ∂ h (p 1 ) j 1 D 0 (a, α) = 2 Jk ∂ h (p 3 ) j 1 ∂ h (p 2 ) j 1 ∂ h (p 1 ) j 1 f 0 h j 1 , f 0 h j 1 - 1 J J j ′ =1 f 0 h j ′ R k×2 + 2 Jk - 2 J 2 k ∂ h (p 2 ) j 1 ∂ h (p 1 ) j 1 f 0 h j 1 , ∂ h (p 3 ) j 1 f 0 h j 1 R k×2 + ∂ h (p 3 ) j 1 ∂ h (p 1 ) j 1 f 0 h j 1 , ∂ h (p 2 ) j 1 f 0 h j 1 R k×2 + ∂ h (p 3 ) j 1 ∂ h (p 2 ) j 1 f 0 h j 1 , ∂ h (p 1 ) j 1 f 0 h j 1 R k×2
By Cauchy-Schwarz inequality we have,

∂ h (p 2 ) j 1 ∂ h (p 1 ) j 1 f 0 h j 1 , ∂ h (p 3 ) j 2 f 0 h j 2 R k×2 ≤ ∂ h (p 2 ) j 1 ∂ h (p 1 ) j 1 f 0 h j 1 R k×2 ∂ h (p 3 ) j 2 f 0 h j 2 R k×2 ≤ e 2A f 0 2 R k×2 . (B.6)
and

∂ h (p 3 ) j 1 ∂ h (p 2 ) j 1 ∂ h (p 1 ) j 1 f 0 h j 1 , f 0 h j 1 - 1 J J j ′ =1 f 0 h j ′ R k×2 ≤ ∂ h (p 3 ) j 1 ∂ h (p 2 ) j 1 ∂ h (p 1 ) j 1 f 0 h j 1 R k×2 f 0 h j 1 - 1 J J j ′ =1 f 0 h j ′ R k×2 ≤ 2e 2A f 0 2 R k×2 (B.7)
For κ = (κ 1 , . . . , κ 2J ) ∈ N 2J , denote by |κ| = κ 1 + . . . + κ 2J and

(∂ h ) κ = (∂ a 1 ) κ 1 (∂ α 1 ) κ 2 . . . (∂ a J ) κ 2J -1 (∂ α J ) κ 2J .
Then, the differential of order 3 of D 0 at c ∈ Θ a,α 0 applied at (a, α) ∈ R 2J writes as

∇ 3 D 0 (c) (a, α) = |κ|=3 (∂ h ) κ D 0 (c)h κ
where h κ = a κ 1 1 α a 2J -1 J α

. This formula together with equations (B.6) and (B.7) give,

∇ 3 D 0 (c) (a, α) = 2 p 1 ,p 2 ,p 3 =1 J j 1 =1 ∂ h (p 3 ) j 1 ∂ h (p 2 ) j 1 ∂ h (p 1 ) j 1 D 0 (c)h (p 1 ) j 1 h (p 2 ) j 1 h (p 3 ) j 1 + 3 J j 1 =j 2 =1 ∂ h (p 3 ) j 2 ∂ h (p 2 ) j 1 ∂ h (p 1 ) j 1 D 0 (c)h (p 1 ) j 1 h (p 2 ) j 1 h (p 3 ) j 2 ≤ 2e 2A 1 k f 0 2 R k×2 2 p 1 ,p 2 ,p 3 =1 4 J J j 1 =1 |h (p 1 ) j 1 h (p 2 ) j 1 h (p 3 ) j 1 | + 6 J 2 J j 1 =j 2 =1 |h (p 1 ) j 1 h (p 2 ) j 1 h (p 3 ) j 2 | ≤ 2δe 2A 1 k f 0 2 R k×2 2 p 1 ,p 2 =1 4 J J j 1 =1 |h (p 1 ) j 1 h (p 2 ) j 1 | + 6(J -1) J 2 J j 1 =1 |h (p 1 ) j 1 h (p 2 ) j 1 | ≤ 20δe 2A 1 k f 0 2 R k×2 1 J J j=1 2 p 1 ,p 2 =1 |h (p 1 ) j h (p 2 ) j | = 20δe 2A 1 k f 0 2 R k×2 1 J J j=1 2 p 1 =1 |h (p 1 ) j | 2 ≤ 40δe 2A 1 k f 0 2 R k×2 1 J J j=1 2 p 1 =1 |h (p 1 ) j | 2 .
And the claim is now proved.

Lemma B.4. For all f ∈ R k×2 and (a 1 , α 1 , b 1 ), (a

2 , α 2 , b 2 ) ∈ [-A, A] × [-A, A] × [-B, B] 2 , let g i .f = e a i f R α i + 1 k ⊗ b i , i = 1, 2
. Then, we have

1 k g 1 .f -g 2 .f 2 R k×2 ≤ C(A, f ) (a 1 , α 1 , b 1 ) -(a 2 , α 2 , b 2 ) 2 R 4 ,
where C(A, f ) = 2 max{4e 4A 1 k f 2 R k×2 , 1}. Proof. We have

1 k g 1 .f -g 2 .f 2 R k×2 ≤ 2e 2A 1 k e a 1 -a 2 f R α 1 -α 2 -f 2 R k×2 + 2 1 k 1 k ⊗ (b 1 -b 2 ) 2 R k×2 (B.8) Let now F (a, α) = 1 √ k e a f R α -f R k×2 . We have |∂ a F (a, α)| = 1 √ k e a f Rα-f R k×2 | e a f R α , e a f R α - f R k×2 | ≤ e A √ k f R k×2 and |∂ α F (a, α)| = 1 √ k e a f Rα-f R k×2 | e a f R α+ π 2 , e a f R α -f R k×2 | ≤ e A √ k f R k×2 . The Euclidean norm of the gradient of F satisfies ∇F (a, α) R 2 = |∂ a F (a, α)| 2 + |∂ α F (a, α)| 2 ≤ √ 2e A 1 √ k f R k×2 . Since we have |F (a, α)| = |F (a, α) -F (0, 0)| ≤ √ 2e A 1 √ k f R k×2 (a, α) R 2 , equation (B.8) yields 1 k g 1 .f -g 2 .f 2 R k×2 ≤ 2 max 4e 4A 1 k f 2 R k×2 , 1 |a 1 -a 2 | 2 + |α 1 -α 2 | 2 + |b (1) 1 -b (1) 2 | 2 + |b (2) 1 -b (2) 2 | 2 ,
which concludes the proof.

  λ , αλ , b).

Figure 3 :

 3 Figure 3: (a) Plot of the mean pattern f = (f (1) , f (2) ) used in the simulations with k = 1024. (b)The first coordinates f(1) . (c) The second coordinates f(2) .

  Recall our notations: the error term is denoted by ζ = (ζ (1) , ζ (2) ) ∈ R k×2 and the vectorized version of ζ is denoted by ζ ∈ R 2k . The simulations were run with three different kinds of noise. White noise : the random variable ζ = (ζ

Figure 4 :Figure 5 :

 45 Figure 4: Example of data generated by model (1.1) with white noise.

Figure 6 :

 6 Figure 6: Example of observations generated by model (1.1) with the highly correlated noise.

  f * 2 R k×2 (estimation with smoothing). The missing boxplots for k = 20 belong to the range [0.7, 1.9].

Figure 7 :

 7 Figure7: Estimation of the mean pattern with the white noise. Boxplot in red correspond to J = 10, in green to J = 100 and in blue to J = 500.

Figure 8 :

 8 Figure8: Estimation of the mean pattern with the weakly correlated noise. Boxplot in red correspond to J = 10, in green to J = 100 and in blue to J = 500.

  of 1 k f -f * 2R k×2 (estimation without smoothing). The missing red boxplot for k = 20 and J = 10 belongs to the range[70, 110].

  f * 2 R k×2 (estimation with smoothing). The missing boxplots for k = 20 belong to the range [0.7, 1.9].

Figure 9 :

 9 Figure9: Estimation of the mean pattern with the highly correlated noise. Boxplot in red correspond to J = 10, in green to J = 100 and in blue to J = 500.

  5), we have explicit expressions of b * Θ 0 = argmin b∈Θ b 0 D(a * , α * , b) and of bλ = argmin b∈Θ b 0 M (â λ , αλ , b). Indeed, we have, b

  3)The expressions of the gradient and the Hessian of D simplify on the set (a * , α * , b * ) * G, see Lemma 3.1. For anyg 0 = (a 0 , α 0 , b 0 ) ∈ R × [-π, π[×R 2 , we have f g * f g * j ′ .g 0 = e -a 0 (f -1 k ⊗ b 0 )R -α 0 -1 J J j ′ =1 e -a 0 (f -1 k ⊗ b 0 )R -α 0 = 0. It yields that for all (a 0 , α 0 , b 0 ) ∈ R × [-π, π[×R2 we have ∇D((a * , α * , b * ) * (a 0 , α 0 , b 0 )) = 0, and ∂ g (p 2 ) a * , α * , b * ) * (a 0 , α 0 , b 0 ))

Acknowledgements

The authors acknowledge the support of the French Agence Nationale de la Recherche (ANR) under reference ANR-JCJC-SIMI1 DEMOS.