Jérémie Boulanger

Nicolas Le Bihan

Stefan Catheline

Vincent Rossetto
email: vincent.rossetto@grenoble.cnrs.fr

Observation of a non-adiabatic geometric phase for elastic waves

We report the experimental observation of a geometric phase for elastic waves in a waveguide with helical shape. The setup reproduces the experiment by Tomita and Chiao (Phys. Rev. Lett. 57, 1986) that showed first evidence of a Berry phase, a geometric phase for adiabatic time evolution, in optics. Experimental evidence of non-adiabatic geometric has been reported in quantum mechanics. We have performed an experiment to observe the polarization transport of classical elastic waves. In a waveguide, these waves are polarized and dispersive. Whereas the wavelength is of the same order of magnitude as the helix' radius, no frequency dependent correction is necessary to account for the theoretical prediction. This shows that in this regime, the geometric phase results directly from geometry and not from a correction to an adiabatic phase.

Introduction

Polarization is a feature shared by several kinds of waves: light and elastic waves, for instance, have two transverse polarization modes [START_REF] Born | Principles of Optics[END_REF]. The polarization degrees of freedom are constrained to lie in the plane orthogonal to the direction of propagation. This constraint is responsible, in optics, for the existence of a geometric phase. Geometric phases of different kinds have been discovered after the Berry phase [START_REF] Berry | Quantal phase factors accompanying adiabatic changes[END_REF][START_REF] Hannay | Angle variable holonomy in adiabatic excursion of an integrable Hamiltonian[END_REF][START_REF] Aharonov | Phase change during a cyclic quantum evolution[END_REF][START_REF] Anandan | Geometric angles in quantum and classical physics[END_REF][START_REF] Zwanziger | Berry's phase[END_REF][START_REF] Shapere | Geometric phases in physics[END_REF][START_REF] Bhandari | Polarization of light and topological phases[END_REF].

The geometric phase of light was first experimentally observed by Tomita and Chiao [START_REF] Tomita | Observation of Berry's topological phase by use of an optical fiber[END_REF] using optical fibers. It was later suggested that a geometric phase should exist for any polarized waves [START_REF] Segert | Photon Berry's phase as a classical topological effect[END_REF]. This Letter discusses the case of elastic waves when the adiabatic conditions are not fulfilled, a situation which can not be reached in Optics. For polarized waves in a waveguide, the geometric phase differs from zero only if the shape of the waveguide is three-dimensional and is defined even if the time evolution is not cyclic [START_REF] Samuel | General setting for Berry's phase[END_REF].

The fundamental origin of geometric phases lies in the geometric description of the phase space. The existence of a geometric phase is related to the curvature, either local or global, of the phase space. In a first attempt to classify geometric phases, Zwanziger et al. distinguish adiabatic geometric phases [START_REF] Zwanziger | Berry's phase[END_REF], the main example of which is a spin in a magnetic field [START_REF] Berry | Quantal phase factors accompanying adiabatic changes[END_REF]. The geometric phase for the spin is defined if the system evolves adiabatically, such that transitions between spin states are negligible. The direction of the magnetic field must therefore evolve at a rate 1/T much smaller than the oscillation frequency between spin eigenstates.

Consider the case of waves propagating in a curved waveguide. The role of the magnetic field's direction is played by the direction of propagation and the phase between the spin eigenstates is the orientation of the linear polarization of the wave. The adiabatic approximation imposes that the evolution rate 1/T is much smaller than the wave frequency. In the first experimental evidence for the adiabatic geometric phase, performed in optics by Tomita and Chiao [START_REF] Tomita | Observation of Berry's topological phase by use of an optical fiber[END_REF], the frequency of light ν was indeed several orders of magnitude larger than 1/T . Photon spin flip is negligible, therefore the adiabaticity conditions are fulfilled. In Foucault's pendulum [START_REF] Bergmann | Foucault pendulum through basic geometry[END_REF], a renowned case of classical geometric phase, these conditions are fulfilled as well. Some geometric phases do not require adiabaticity, such as Pancharatnam phase [START_REF] Pancharatnam | Generalized theory of interference and its applications i: coherent pencils[END_REF] and the Aharonov-Anandan quantum phase [START_REF] Aharonov | Phase change during a cyclic quantum evolution[END_REF] or certain canonical classical angles [START_REF] Anandan | Geometric angles in quantum and classical physics[END_REF][START_REF] Berry | Classical non-adiabatic angles[END_REF]. Geometric phases have been observed in many fields of science and called with different names. In classical mechanics, the geometric phase for adiabatic invariants is often referred to as Hannay angle [START_REF] Hannay | Angle variable holonomy in adiabatic excursion of an integrable Hamiltonian[END_REF][START_REF] Berry | Geometric angle for rotated rotators, and the hannay angle of the world[END_REF], in knot theory and DNA physics, the name writhe is mostly used [START_REF] Agarwal | Computing the writhing number of a polygonal knot[END_REF][START_REF] Rossetto | Writhing geometry of open DNA[END_REF].

We consider from now on elastic waves, which polarization state can be represented as a combination of two linearly polarized states. These are classical waves, quantum transition between polarization eigenstates is not possible. The adiabaticity condition ν ≫ 1/T should therefore not be required to observe a geometric phase. In our experiment, indeed, the frequency and the evolution rate have the same order of magnitude. In this Letter, we briefly introduce the geometric phase in a purely geometric picture and compute its value along a helix. We present an experiment designed to measure the geometric phase of elastic waves in a helical waveguide with the condition ν ≃ 1/T . Without loss of generality, we only consider linear polarization.

Although the experiment we investigate has many common points with previous studies, some distinctions must be pointed out. Contrary to light polarization, the phase cannot be interpreted as a quantum phase difference between two eigenstates. There is no established classification of geometric phase but as the system we study is purely classical and neither cyclic nor adiabatic, the observed geometric phase cannot be rigorously identified as a Berry, Pancharatnam, Aharonov-Anandan or Wilczek-Zee phase, for instance. The classical geometric angles [START_REF] Anandan | Geometric angles in quantum and classical physics[END_REF][START_REF] Berry | Classical non-adiabatic angles[END_REF] follow from the action-angle representation of the system, which is valid for the motion of a material point of the waveguide, but does not rigorously apply to the elastic wave transport.

The geometric phase for polarized waves

A wave travelling along a straight path keeps a constant polarization along the trajectory; if the direction of propagation is not constant, as polarization is ascribed to remain in the orthogonal plane, it evolves along the path. The path-dependent transformation transporting polarization must be, for physical reasons, linear, reversible and continuous. There is only one transformation satisfying these requirements: parallel transport [START_REF] Segert | Photon Berry's phase as a classical topological effect[END_REF]. Along the path followed by the wave, the direction of propagation is represented as a point on the unit sphere, polarization is represented in the tangent plane to the sphere and transported in the sphere's tangent bundle, see Figure 1.

Let us consider a geodesic on the unit sphere, i.e. an arc of a great circle. If polarization is colinear or orthogonal to the great circle, it remains so along the geodesic to preserve symmetry. Any polarization is a linear combination of these two particular linear polarizations. The linearity of parallel transport implies, for linear polarizations, that if one rotates the initial polarization in the initial tangent plane, the parallel transported polarizations rotate by the same angle in all tangent planes of the tangent bundle. Parallel transport can be extended to any smooth and piecewise differentiable trajectory of the tangent vector on the unit sphere by discretizing the path into elementary arcs of great circles [START_REF] Rossetto | Writhing geometry of stiff polymers and scattered light[END_REF].

A vector parallel transported along a closed trajectory does not necessarily have the same orientation in the tangent plane at the starting point after one stride along the trajectory. The angle between the initial and final vectors is algebraically equal to the area enclosed by the trajectory on the sphere [START_REF] Nakahara | Geometry, topology and physics[END_REF]. In the Figure 1, this area is equal to π/2. It is also known as the anholonomy of the trajectory in the phase space.

Computation of the phase

All the distances, denoted by s, are given in arclength along the helix and the accelerometers position is taken as the origin. We call R the radius of the helix and P its pitch, and define L = (2πR) 2 + P 2 . The angle θ between the direction of propagation t(s) and the helix main axis is given by sin θ = P/L. The Frénet torsion τ = (2π cos θ)/L of the helix is constant. 2πs/L is the azimuth angle spanned by the tangent vector along a distance s on the helix. Fuller's theorem states that the geometric phase difference between two trajectories is equal to the area spanned by the tangent vector on the unit sphere during one smooth deformation of one trajectory into the other. We chose as reference the trajectory θ = π/2 (the limit case where the helix is a flat circle) for which it is known that the geometric phase vanishes because the trajectory is two-dimensional. After deformation of the reference into the helix, Fuller's theorem yields a phase proportional to 2πs/L and the proportionality coefficient is cos θ (given by the classic geometry formula for the spherical cap area). We obtain a geometric phase Φ(s) of:

Φ(s) = (2πs/L)(cos θ) = τ s. (1)
From an intrinsic point of view, parallel transport corresponds to transporting a polarization vector v while keeping it constant for an imaginary walker striding along the trajectory. In the language of differential geometry the vector's covariant derivative must equal zero along the trajectory [START_REF] Nakahara | Geometry, topology and physics[END_REF]. We obtain the equation of parallel transport:

Dv = dv -(dv • t) t = 0. ( 2
)
The components of v in the Frénet-Serret frame follow a differential equation describing a rotation at the rate τ , leading to a phase:

Φ(s) = τ s. (3)

Experiment

In order to reproduce the setup used by Tomita and Chiao [START_REF] Tomita | Observation of Berry's topological phase by use of an optical fiber[END_REF], we use a metallic spring as a waveguide for elastic waves, taken from a car's rear damper. It has a circular section of r = 13.5 mm making five coils of radius R = 75 ± 1 mm and with a pitch of P = 91.5 ± 1 mm. The direction of propagation makes a constant angle θ = 1.379 rad with the helix' main axis. As the crosssection of the helix is circular, the flexural modes are degenerated. The helix is suspended to two strings to isolate the system. We use two accelerometers (Brüel & Kjaer, type 4518-003) located at one end of the spring and record vibrations in two orthogonal directions (see Figure 3). The sampling frequency being set to 50 kHz, the information in the signal is available up to 25 kHz. Considering the accelerometers' power spectrum in their stability range 1kHz -25kHz (see Figure 4), frequencies above 5 kHz can be ignored. The signal is amplified (Brüel & Kjaer amplifier, type 2694) before signal processing.

We record the waves generated by 32 equally spaced sources, which distance from the accelerometers ranges from s = 5 cm to s = 1.45 m. Making weak impacts on the metal spring creates bending waves that are linearly polarized. Impacts are made radially with respect to the helix. The source signal is generated manually, by gently hitting the waveguide with a hammer at the different source positions. Polarization depending only on the amplitude ratio measured by the accelerometers, it is not sensitive to the energy of the source.

The wave propagation modes in helical waveguide constitute a difficult problem. Solutions can only be found numerically [START_REF] Treyssède | Numerical investigation of elastic modes of propagation in helical waveguides[END_REF] and yield a complex mode structure. The fundamental mode is a flexural mode with linear polarization. Polarization losses appear rapidely after the first waveforms. In order to measure the geometric phase, we use only the first periods of oscillation of the signal, before depolarization is complete. In figure 2, the time windowed signals used for polarization orientation estimation are presented, together with the complete waveforms and associated polarization parametric plots. In Figure 4, we show an example of two parametric polarization plots for signals obtained with the source at different distances from the accelerometers. The relative angle between polarization orientations is the geometric phase difference.

Bending waves in the waveguide are dispersive. We therefore filtered the signal using non-overlapping bandpass filters centered at frequencies: 1.4 kHz, 2.5 kHz, 4.0 kHz and 5.8 kHz. These values correspond to maxima of the power spectrum displayed in figure 4 (left). The corresponding velocities are displayed in the Table 1.

The polarization orientation is obtained from the records using a principal component analysis [START_REF] Vidale | Complex polarization analysis of particle motion[END_REF]. This technique consists in obtaining the eigenvectors of the cross-correlation between the two orthogonal signals (recorded on the two accelerometers). The eigenvector with the largest eigenvalue gives the polarization direction. The geometric phase Φ(s) is thus the orientation of polarization measured when the source is located at distance s from the accelerometers (see Figure 3).

We perform several measurements for each position of the source, in order to reduce effects of external perturbations and variations of the source signal. In our results, we observe oscillations that are due to imperfections related to the ν(kHz) c(m.s -1 ) wavelength (m) τ (rad.m -1 ) η ǫ (rad) 1.4 accelerometers response and setup. Denoting by ǫ the shift from orthogonality between the accelerometers and 1 + η the ratio of their mechanical couplings, the principal components analysis yields:

Φ(s) = Φ(0) + τ s + sin Φ(s) η cos Φ(s) -ǫ sin Φ(s) . (4)
We obtain a linear coefficient τ ≃ 2.5 rad.m -1 for the signals filtered at four differents frequencies (bandwidth ±20 %). Numerical values are presented in Table 1.

Results and discussion

The theoretical Frénet torsion of the helix used in the experiment is τ = 2.49 ± 0.1 rad.m -1 . The fits performed from experimental data give an estimation of τ = 2.50 ± 0.1 rad.m -1 . We do not observe a significant dependency of τ with respect to the frequency or the velocity (see Table 1), which justifies the denomination "geometric phase" : the effect we observe is solely due to the geometry of the waveguide.

Apart from bending waves, compression waves and torsion waves can also propagate in the helix. In a straight rod, these propagation modes are not coupled and bending waves remain polarized at long times. In helices, the coupling between compression waves and torsion waves increases with curvature and torsion. Bending waves and torsional waves are therefore partially converted into each other, back and forth, during the propagation. This explains depolarization in our measurements and why we consider only the first oscillations of the signals.

We have filtered the polarized waves at four frequencies of the order of magnitude of the rate at which the propagation direction evolves 1/T = c/L. Such evolution rates imply that the regime is not adiabatic. Thanks to the dispersivity of the bending waves in the system, the adiabatic parameter c/L can be varied by changing the filtering central frequency. Therefore filtering at several frequencies is equivalent to explore the non-adiabatic regime. No significant variations of the results are observable in the frequency range, which is an experimental evidence of the non-adiabaticity of the geometric phase, complementary to the theoretical derivation. Converted into lengths, this means that the wavelength of the bending wave is of the same order of magnitude as the length of a helix coil.

One may argue that the conversions between modes observed in the experiment are the classical equivalent of the transition between spin eigenstates of the quantum problem. By nature, however, these transitions are different. In the quantum problem, the transitions are of the first order in time but for the elastic waves, depolarization processes are of the second order, because an intermediate, unpolarized state is necessary. According to the dispersion relation in a helical waveguide [START_REF] Treyssède | Numerical investigation of elastic modes of propagation in helical waveguides[END_REF], the intermediate state is mainly the torsional mode, and second the longitudinal compression mode.

Conclusion

We have demonstrated the existence of a geometric phase for elastic waves in a waveguide far from the adiabatic regime. We have investigated the results at several frequencies, or equivalently with several values of the adiabatic parameter, and observed no significant variations of the phase's value. This is an experimental evidence, in addition to the theoretical derivation, that the observed phase is then not an adiabatic geometric phase with non-adiabatic corrections, but a non-adiabatic geometric phase for classical systems. Because there is no need to invoke the adiabatic approximation to preserve the polarization state, the geometric phase should extend to all frequencies in smooth waveguides with circular sections, a interesting result for waves with very large wavelengths. In quantum mechanics, the Aharonov-Anandan phase [START_REF] Aharonov | Phase change during a cyclic quantum evolution[END_REF] shares the same non-adiabatic aspect as the non-quantum geometric phase studied in this letter.

In nature, polarized elastic waves, such as seismic S waves (shear waves) are observed under certain conditions, and the concept of geometric phase applies. The degree of polarization and the statistics of the polarization orientation in seismic signals created by a polarized source contain informations concerning the disorder of propagating medium that can be understood in terms of geometric phases.

Figure 1 :

 1 Figure 1: (Left) A trajectory represented on the sphere of tangent vector. Linear polarization is represented as a direction in the plane orthogonal to the direction of propagation and the trajectory is made of three arcs of great circles forming a right angle with each others. Both at the start and the end of the trajectory (upper pole), the tangent vector points upwards. One observes the anholonomy of parallel transport: The double arrow indicating the polarization direction, although it is parallel transported all along the trajectory back to its initial position, has different directions at the beginning and at the end of the displacement. (Right) Same trajectory in real space.

Figure 2 :

 2 Figure 2: Example of vertical (solid line) and horizontal (dashed line) displacement signals recorded during the experiment. (a): complete recorded waveforms associated to one source. (b): time windowed signals considered for polarization orientation estimation. (c): associated polarization parametric plot. The truncated (time windowed) signal is chosen in order to isolate direct linearly polarized wave from other modes and reflected waves. Only few first oscillations are considered.

Figure 3 :

 3 Figure 3: (Left) Setup of the experiment. Accelerometers are represented as cubes, they record acceleration in two orthogonal directions (small arrows). The large arrow symbolizes the source, a radial impact on the helix. The dots indicate some of the source positions. (Right) Geometry of the transformation from the reference trajectory to the helix trajectory. During the transformation the tangent vector for s fixed follows a meridian (vertical arrow on the figure) between colatitudes π/2 and θ. A displacement of s corresponds to an azimuth angle of 2πs/L (horizontal arrow). The area spanned by the tangent vector during the deformation is equal to the geometric phase difference computed by Fuller's formula.

Figure 4 :

 4 Figure 4: (Left) Power spectrum of the first oscillations of the signal corresponding to direct arrival of the flexural mode. Circles indicate the maxima used as central filtering frequencies in the signal analysis. (Right) Parametric plot of the first oscillations (direct arrival) for sources located at s = 0.36 m and s = 1.00 m. The relative angle between polarization orientations is the geometric phase difference.

Figure 5 :

 5 Figure 5: Measured geometric phase as a function of arc length along the helix for signals filtered at 1.4 kHz, 2.5 kHz, 4.0 kHz and 5.8 kHz. The dashed line represents the linear regression with correction for the coupling and orthogonality of the accelerometers. The origins of the curves are arbitrary.

Table 1 :

 1 Frequencies, velocities, wavelengths and fit parameters for the filtered signals. The values of τ have a confidence range of ±0.1 rad.m -1 .

		2 10 3	1.4	2.53	0.31	-0.03
	2.5	3 10 3	1.2	2.35	0.45	0.03
	4.0	4 10 3	1.0	2.45	0.39	0.05
	5.8	3.5 10 3	0.6	2.50	0.38	0.07

Acknowledgments

The authors would like to thank É. Larose for discussions and B. de Cacqueray for his experimental help. This work was partially funded by the ANR-JC08-313906 SISDIF and CNRS/PEPS-PTI grants.


Jérémie Bigot

Raúl Gouet

Alfredo Lopez

Geometric PCA of images

Keywords: Geometric variability, Principal component analysis, Mean pattern estimation, Fréchet mean, Image registration, Deformable models. AMS classifications: Primary 62H25, secondary 62H35. in Chile

We describe a method for analyzing the shape variability of images, called geometric PCA. Our approach is based on the use of deformation operators to model the geometric variability of images around a reference mean pattern. This leads to a new algorithm for estimating shape variability. Some numerical experiments on real images are proposed to highlight the benefits of this approach. The consistency of this procedure is also analyzed in statistical deformable models.

Introduction

In many applications observations are in the form of a set of n gray-level images y 1 , . . . , y n (e.g. in geophysics, biomedical imaging or in signal processing for neurosciences), which can be considered as square-integrable functions on a domain Ω, a convex subset of R d . Such data are generally two or three dimensional images. In many situations the observed images share the same structure. This may lead to the assumption that these observations are random elements, which vary around the same mean pattern. Estimating such a mean pattern and characterizing the modes of individual variations around this common shape, is of fundamental interest. Principal components analysis (PCA) is a widely used method to estimate the variations in intensity of images around the usual Euclidean mean ȳn = 1 n n i=1 y i . However, such data typically exhibit not only a classical source of photometric variability (a pixel intensity changes from one image to another) but also a (less standard) geometric source of variability in shape, which cannot be recovered by standard PCA. The goal of this paper is to describe a method for analyzing the shape variability of images, called geometric PCA. We present some numerical experiments to highlight the benefits of this approach in applications, and we study the consistency of this procedure in statistical deformable models.

PCA in a Hilbert space

First, let us introduce some tools and notations to perform a standard PCA in a Hilbert space that will be used throughout the paper. Let H a separable Hilbert space endowed with inner product •, • and associated norm • . Let Z be an H-valued random variable. If E Z < +∞, then Z has expectation EZ ∈ H which happens to be the unique element satisfying EZ, h = E Z, h , for all h ∈ H. If E Z 2 < +∞, then the (population) covariance operator K : H → H corresponding to Z is given by

Kh = E Z -EZ, h (Z -EZ) for h ∈ H.
Moreover, the operator K is self-adjoint, positive semidefinite and trace-class. Hence, K is compact, with nonnegative (population) eigenvalues (γ λ ) λ∈Λ and orthonormal (population) eigenvectors (φ λ ) λ∈Λ and such that Kh = λ∈Λ γ λ h, u λ u λ where Λ = {1, . . . , dim(H)}, if dim(H) < ∞ or Λ = N otherwise. If we assume that the eigenvalues are arranged in decreasing order γ 1 ≥ γ 2 ≥ . . . ≥ 0, the λ-th mode of variation of the random variable Z, is defined as the function EZ + ρ √ γ λ u λ , for λ ∈ Λ, where ρ is a real playing the role of a weight parameter (typically one takes ρ between -1 and 1). Thus, the PCA of Z is obtained by diagonalizing the covariance operator K. Now, let z 1 , . . . , z n ∈ H and Kn : H → H, their (empirical) covariance operator defined as

Kn h = 1 n n i=1 z i -zn , h (z i -zn ), for h ∈ H,
where, zn = 1 n n i=1 z i . As Kn is also self-adjoint, positive semidefinite and compact, it admits the decomposition Kn h = λ∈Λ γλ h, ûλ ûλ , where γ1 ≥ γ2 ≥ . . . ≥ 0 are the (empirical) eigenvalues, and (û λ ) λ∈Λ is the set of (empirical) orthonormal eigenvectors of Kn . Now, the λ-th (empirical) mode of variation of the data z 1 , . . . , z n is defined as zn + ρ √ γλ ûλ , for λ ∈ Λ. Thus, the empirical PCA of the data z 1 , . . . , z n is obtained by diagonalizing the operator Kn .

If H is finite-dimensional, diagonalizing Kn corresponds to the empirical PCA for vectors in a finite dimensional Euclidean space. If H = L 2 (Ω) = f : Ω → R, f 2 2 := Ω |f (x)| 2 dx < ∞ diagonalizing Kn is usually referred to as the method of functional PCA in nonparametric statistics (see e.g. [START_REF] Ramsay | Functional data analysis[END_REF] for an introduction to functional data analysis). Various authors (see e.g. [START_REF] Dauxois | Asymptotic theory for the principal component analysis of a vector random function: some applications to statistical inference[END_REF], [START_REF] Silverman | Smoothed functional principal components analysis by choice of norm[END_REF] and references therein), under standard probabilistic assumptions on the data, have studied the consistency of empirical PCA in Hilbert spaces and proposed sufficient conditions to ensure that the empirical eigenvalues and eigenvectors converge to the population eigenvalues and eigenvectors as n → +∞.

Empirical PCA, with H = L 2 (Ω), applied to a set of n images y 1 , . . . , y n , is a method to compute the principal directions of photometric variability of the y i 's around the usual Euclidean mean ȳn . However, in many situations, images also exhibit a large variability in shape, see the example of images of handwritten digits displayed in Figure 1 and Figure 2. In such cases, the standard Euclidean mean ȳn is not a satisfying estimator of the typical shape of each individual image, see Figure 1 In particular, the second empirical modes of variation are no longer a single digit but rather the superposition of two digits in different orientations.

It is well known that the usual Euclidean mean ȳn is the minimizer of the sum-of-squares Euclidean distances to each of the data points, namely ȳn = arg min

f ∈L 2 (Ω) 1 n n m=1 f -y m 2 2 .
The idea underlying PCA is that the Hilbert space L 2 (Ω), equipped with the standard inner product, is well suited to model natural images. However, the set of such objects (as those in Figures 12) typically cannot be considered as a linear sub-space of L 2 (Ω). Therefore, the Euclidean distance f 1 -f 2 2 is generally not well suited, since it is not adapted to the geometry of the set to which the images f 1 and f 2 truly belong. Actually one can see that the images in Figures 12have mainly a geometric variability in space, which is much more important than the photometric variability.

Organization of the paper

In Section 2, we define the geometric PCA of random images. We describe the numerical implementation of this method, and we illustrate its advantages over standard PCA for the analysis of handwritten digits. A review of PCA-like methods to analyze geometric variability of images is also given. In Section 3, we study the consistency of geometric PCA in statistical deformable models. We conclude the paper in Section 4 by a short discussion. All proofs are gathered in a technical Appendix.

Geometric PCA

For convenience, we prefer to present the ideas of geometric PCA under the assumption that the images are observed on a continuous domain Ω. In practice, such data are obviously observed on a discrete set of time points or pixels. However, assuming that the data are random elements of L 2 (Ω) is more convenient for dealing with statistical aspects of an inferential procedure, as it avoids the treatment of the bias introduced by any discretization of the domain Ω.

Grenander's pattern theory of deformable templates

Following the ideas of Grenander's pattern theory (see [START_REF] Grenander | Pattern Theory: From Representation to Inference[END_REF] for a recent overview), one may consider that the data y 1 , . . . , y n are obtained through the deformation of the same reference image. In this setting, images are treated as points in L 2 (Ω) and the geometric variations of the images are modeled by the action of Lie groups on the domain Ω. Recently, there has been a growing interest in Lie groups of transformations to model the geometric variability of images (see e.g. [BMTY05, TY05a, TY05b, You10] and references therein), and applications Grenander's pattern theory leads to the construction of non-Euclidean distances between images. In this paper, we propose to model shape variability through the use of deformation operators (acting on Ω) that are parameterized by a separable Hilbert space V, with inner product •, • . We also assume that Ω is equipped with a metric denoted by d Ω .

Definition 2.1. Let V be a Hilbert space. A deformation operator parameterized by V is a mapping ϕ : V × Ω → Ω such that, for any v ∈ V, the function x → ϕ(v, x) is a homeomorphism on Ω. Moreover, ϕ(0, •) is the identity on Ω and, for any v ∈ V, there exists

v * ∈ V such that ϕ -1 (v, •) = ϕ(v * , •).
In this paper, we will study as illustrative examples of deformation operators the cases of translations, rigid deformations and non-rigid deformations generated by stationary vector fields.

Translations: Let Ω = [0, 1) d , for some integer d ≥ 1 and V = R d . Let also ϕ(v, x) = (mod(x 1 + v 1 , 1), . . . , mod(x d + v d , 1)) , (2.1) with ϕ -1 (v, x) = (mod(x 1 -v 1 , 1), . . . , mod(x d -v d , 1)) ,
for all v = (v 1 , . . . , v d ) ∈ R d and x = (x 1 , . . . , x d ) ∈ Ω, where mod(a, 1) denotes the modulo operation between a real a and 1. Clearly, ϕ(0, •) is the identity in Ω and v * = -v. Moreover, it can be shown (see Section 3.3) that ϕ(v, •) is an homeomorphism.

Rigid deformations of 2D images:

Let Ω = R 2 and V = R × R 2 . Let also ϕ(v, x) = R α x + b, with ϕ -1 (v, x) = R -α (x -b), for all v = (α, b) ∈ R × R 2
, and x ∈ R 2 where R α is the rotation matrix of angle α and b ∈ R 2 defines a translation. Observe that, ϕ((0, 0), •) is the identity in Ω and v * = (-α, -R -α b).

Clearly, ϕ(v, •) is an homeomorphism.

Diffeomorphic deformations generated by stationary vector fields: Let Ω = [0, 1] d for some integer d ≥ 1, and V a separable Hilbert space of smooth vector fields such that V is continuously embedded the set of functions v : Ω → R d which are continuously differentiable and such v and its derivatives vanish at the boundary of Ω. For x ∈ Ω and v ∈ V, define ϕ(v, x) as the solution at time t = 1 of the following ordinary differential equation (O.D.E.)

dφ t dt = v(φ t ), (2.2)
with initial condition φ 0 = x ∈ Ω. It is well known (see e.g. [START_REF] Younes | Shapes and diffeomorphisms[END_REF]) that, for any v ∈ V, the function x → ϕ(v, x) is a C 1 diffeomorphism on Ω. The inverse of x → ϕ(v, x) is given by x → ϕ(-v, x), and thus v * = -v. Hence, ϕ(v, •) satisfies the conditions in Definition 2.1.

Registration

Registration is a widely used method in image processing that consist in geometric transforms of a set of images y 1 , . . . , y n ∈ L 2 (Ω), so that they be compared. This method can be described as an optimization problem which amounts to minimizing a dissimilarity functional between images.

Definition 2.2 (Dissimilarity functional). Let ϕ be a deformation operator, as described in Definition 2.1, v = (v 1 , . . . , v n ) ∈ V := V n and y = (y 1 , . . . , y n ), with y i ∈ L 2 (Ω), i = 1, . . . , n.

(a) The template dissimilarity functional corresponding to v, y and f ∈ L 2 (Ω) is defined as

M t (v, y, f ) := 1 n n i=1 Ω y i (ϕ(v i , x)) -f (x) 2 dx. (2.3) (b)
The groupwise dissimilarity functional corresponding to v and y is defined as

M g (v, y) := 1 n n i=1 Ω y i (ϕ(v i , x)) - 1 n n j=1 y j (ϕ(v j , x)) 2 dx.
(2.4)

Template registration of the images y 1 , . . . , y n onto some known template f ∈ L 2 (Ω), is defined as the problem of minimizing the criterion given by the dissimilarity functional (2.3), with respect to v in

V µ := {v = (v 1 , . . . , v n ), v i ∈ V µ } ,
where V µ := {v ∈ V : v ≤ µ}, for some regularization parameter µ ≥ 0 whose choice has to be discussed. Note that imposing the constraint v i ≤ µ allows to explicitly control the norm of the vector v i which is generally proportional to the distance between the deformation ϕ(v i , •) and the identity.

On the other hand, groupwise registration of y 1 , . . . , y n is defined as the problem of minimizing the functional (2.4) with respect to v in U ⊆ V µ . Two possible choices for U , defined in terms of linear constraints on v, are

U 0 := v ∈ V µ , n m=1 v m = 0 and U 1 := v ∈ V µ , v 1 = 0 . (2.5)
Choosing U = U 0 amounts to imposing that the deformation parameters (v 1 , . . . , v n ) used to align the data have an empirical mean equal to zero, while taking U = U 1 corresponds to choosing y 1 as a reference template onto which y 2 , . . . , y n will be aligned.

Geometric PCA applied to a set of images y = (y 1 , . . . , y n ) is the following two step procedure. In the first step, one applies either a template or a groupwise registration, which leads to the computation of

v ∈ arg min v∈Vµ M t (v, y, f ) or v ∈ arg min v∈U M g (v, y).
(2.6)

In the second step, a standard PCA is carried out on v = (v 1 , . . . , vn ), based on the following covariance operator

Kn v = 1 n n i=1 vi -v n , v (v i -v n ), for v ∈ V, (2.7) with v n = 1 n n i=1 vi . This operator admits the decomposition Kn v = λ∈Λ κλ v, φλ φλ , (2.8)
where κ1 ≥ κ2 ≥ . . . ≥ 0 and ( φλ ) λ∈Λ are the eigenvalues and orthonormal eigenvectors of Kn . We now state the definition of geometric PCA of a set of images.

Definition 2.3 (Geometric PCA). Let ϕ be a deformation operator parametrized by V, as described in Definition 2.1. Let (κ λ , φλ ) λ∈Λ be the eigenvalues and orthonormal eigenvectors of the operator Kn in (2.7). For λ ∈ Λ, the λ-th empirical mode of geometric variation of the data y 1 , . . . , y n is the homeomorphism ψλ : Ω → Ω defined by

ψλ (x) = ϕ -1 (v n + κλ φλ , x), x ∈ Ω.
(2.9)

We also denote ψλ,ρ (x) = ϕ -1 (v n + ρ √ κλ φλ , x), where ρ ∈ R is a weighting value.

After the registration step, we obtain a set of deformed images

y 1 • ϕ(v 1 , •), . . . , y n • ϕ(v n ,
•), each of them aligned either with respect to f in the case of template registration, or with respect to f := 1 n n j=1 y j (ϕ(v j , x)), in the case of groupwise registration. Hence, in the case of template registration, f • ψλ can be used to visualize the λ-th mode of geometric variation of the data. Similarly, in the case of groupwise registration, one uses f • ψλ . Note that f can be interpreted as a mean pattern image. Moreover, the computation of f is closely related to the notion of Fréchet mean of images, recently studied in [START_REF] Bigot | On the consistency of frchet means in deformable models for curve and image analysis[END_REF], from a statistical point of view.

Numerical implementation and application of geometric PCA to handwritten digits data

In this section we explain in detail the implementation of geometric PCA, in the case of groupwise registration, using the class of diffeomorphic deformation described in Section 2.1. The method is applied to a set of n = 30 images, defined on the domain Ω = [0, 1] 2 , taken from the Mnist data base of handwritten digits [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF].

Specification of the Hilbert space of parameter V

We choose V as the vector space of functions from Ω to R 2 , generated by a B-Splines basis of functions, because they have good properties for approximating continuous functions and implementing efficient computations [START_REF] Unser | B-Spline signal processing: Part I-Theory[END_REF][START_REF] Unser | B-Spline signal processing: Part II-Efficient design and applications[END_REF]. Let {b k : Ω → R, k = 1, . . . , p} denote a set of bi-dimensional tensor product B-Splines, with knots defined on a regular grid of Ω, and p some integer whose choice has to be discussed. We define V as the space of vector fields of the form v = p k=1 ṽk b k , where ṽk = (ṽ

(1) k , ṽ(2) k ) ∈ R 2 , k = 1, . . . , p. We denote by v (1) , v (2) : Ω → R the coordinates of v ∈ V, i.e., v(x) = (v (1) (x), v (2) (x)) for x ∈ Ω.
Note that the dimension of V is 2p and that a basis is given by {(b 1 , 0), . . . , (b p , 0), (0, b 1 ), . . . , (0, b p )}.

(2.10)

We endow V with the inner product

u, v := u (1) , v (1) L + u (2) , v (2) L , u, v ∈ V, where u (1) , v (1) L := Ω Lu (1) (x)Lv (1) (x)dx, u (2) , v (2) L := Ω Lu (2) (x)Lv (2) (x)
dx and L is a differential operator. As suggested in [START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF] we take L = γI + α∆, where I is the identity operator, ∆ is the laplacian operator and γ, α are positive scalars. By using the basic properties of differentiation and integration of B-Splines [UAE93a], we derive an explicit formula for computing the inner product in V, that can be implemented using convolution filters. By an adequate design of the B-Spline grid, we ensure that the values of v and its derivatives are zero at the boundary of Ω.

Minimization of the dissimilarity functional D g

In the case of groupwise registration, we minimize the dissimilarity functional (2.4) over the set U 0 defined in (2.2). Thanks to the above choice for V, the minimization of the criterion (2.4) has to be performed over a subset of R 2p . In order to take into account the constraint v i ≤ µ, 1 ≤ i ≤ n, we use a logarithmic barrier approach to obtain an approximate solution.

Then, for the minimization, we use a gradient descent algorithm, with an adaptive step. Such an algorithm requires the computation of the deformation operator ϕ : V × Ω → Ω and its gradient, with respect to the coefficients ṽk = (ṽ

(1)
k , ṽ(2) k ) that parameterize the vector field v. In the case of diffeomorphic deformation operators, ϕ(v, x) corresponds to the solution at time t = 1 of the ODE (2.2). We solve the O.D.E. using a forward Euler integration scheme. For a comparison of different methods for solving such O.D.E., we refer to [START_REF] Bossa | Algorithms for computing the group exponential of diffeomorphisms: Performance evaluation[END_REF]. It can be shown (see Lemma 2.1 in [START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF]), that the gradient of ϕ with respect to the ṽk = (ṽ

(1)
k , ṽ(2) k )'s has a closed-form expression.

Spectral decomposition of the empirical covariance operator

Let v1 , . . . , vn be the vector fields in V obtained after the registration step described above. Recall that the empirical covariance operator of the vi 's is defined as Kn

v = 1 n n i=1 vi -v n , v (v i - v n ), v ∈ V.
In what follows, we describe how to perform the spectral decomposition of Kn .

Let ṽi = (ṽ

(1) i , ṽ(2) i ) with ṽ(1) i = (ṽ (1)
i,1 , . . . , ṽ(1) i,p ) and ṽ(2

) i = (ṽ (2)
i,1 , . . . , ṽ(2) i,p ) being the coefficients of vi with respect to the base (2.10), i.e. vi = p k=1 (ṽ

(1) i,k , ṽ(2) i,k )b k .
We identify the Hilbert space V with R 2p endowed by the inner product ũ, ṽ := ũΣũ t , ũ, ṽ ∈ R 2p , (

where Σ is a 2p × 2p matrix with entries Σ j,k = Σ j+p,k+p := b j , b k L for j, k = 1, . . . , p and Σ j,k := 0 in the other cases. Hence, the operator Kn can be identified with a 2p × 2p matrix Kn , given by Kn := 1 n ṽṽ t Σ, where ṽ is the 2p × n matrix with i-th column equals to ṽt i -1 n n j=1 ṽt j . The matrix Σ is symmetric, hence admits a diagonalization Σ = P ΛP with Λ diagonal matrix and P t P = P P t = I. The idea now is to reduce the problem to a standard diagonalization of the symmetric matrix M := 1 n Λ 1 2 P t ṽṽ t P Λ 1 2 , namely we find the decomposition M = W DW t with D diagonal matrix and W t W = W W t = I. We obtain the following spectral decomposition of Kn respect to the inner product (2.11) Kn = U DU t Σ,

where U := P t W . Remark that the columns of U are orthonormal vectors in R 2p with respect to the inner product (2.11). Indeed, it holds that U t ΣU = I. Finally, we define κλ as the λ-th elements of the diagonal matrix D and we let φλ :

= p k=1 (U k,λ , U k+p,λ )b k , for λ = {1, . . . , 2p}. It can be checked that ( φλ ) 2p
λ=1 are orthonormal vectors of V, and we thus obtain that Kn

v = 2p λ=1 κλ v, φλ L φλ . If we assume that κ1 ≥, . . . , ≥ κ2p , then ψλ,ρ = ϕ -1 (v n + ρ κλ φλ , •)
is the λ-th empirical mode of geometric variation according to Definition 2.3.

Application of geometric PCA to handwritten digits data and comparison with standard PCA

We now describe the application of geometric PCA to handwritten digits, taken from the Mnist data base [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF], based on the numerical framework we had described so far. Also, we illustrate the benefits of geometric PCA over standard PCA.

For the B-Spline base of V, we choose a B-Spline degree equals to 3, as it provides a good trade-off among smoothness and the size of the support. The number of B-Spline knots is p = 81 arranged in a 9 × 9 regular grid. Such value of p provides a fine B-Spline grid with respect to an image size of 28 × 28. For defining the differential operator L, we take γ = 100 and α = 1. Note that γ >> α in order to compensate for scaling factor associated to the inter knot spacing. Finally, we choose a regularization parameter µ = 10. This value was determined experimentally, by trying to obtain a good compromise among regularity of the vector fields and matching of the images, at the registration step.

For each available digit (from 0 to 9), we took n = 30 images and we carried out a geometric PCA on each of these image sets. In this database, for each digit, one observes a large source of geometrical variability that can be modeled by diffeomorphic deformations. We proceed by groupwise registration, as there are no reference images available. To illustrate the advantages of our procedure, we have also carried out a standard PCA of each digit, which amounts to analyzing the photometric variability of the data. Thus we compute ȳn + ρ γλ ûλ , the λ-th standard empirical mode of photometric variation of the data as described in Section 1.1. Figures 3, 4 and5 show the geometric modes of variations by displaying the images

f • ψλ,ρ where f (x) = 1 n n j=1 y j (ϕ(v j , x)),
with λ = 1, 2 and ρ = 2, -2. Results using the standard PCA are also displayed in the same Figures. We observe that geometric PCA better reflects the main modes of variability of the digits. In contrast, standard PCA fails in several cases in representing the geometric variability of some digits, and it results sometimes in a blurring of the images. Also, it can be seen that f is a much better mean pattern of the data than the Euclidean mean ȳn .

PCA-like methods for analyzing geometric variability

The idea of applying a PCA analysis to the resulting transformation parameters after registration of a set of images is at the core of several methods to estimate the geometrical variability of images. In [START_REF] Rueckert | Automatic construction of 3d statistical deformation models using non-rigid registration[END_REF], a linear and finite dimensional space of non rigid transformations is considered as the admissible space of deformations on which PCA is carried out. This is the so-called statistical deformation model that succeed Coote's active shape models. An important limitation of this approach is the lack of invertibility of the deformations. In several cases, the invertibility is a desirable property from the point of view of physical (for instance when analyzing shape variability of a determined organ) and mathematical modelling. Within the context of linear space of deformations, the non-invertibility issue had been addressed by enforcing the positivness of the Jacobian determinant [START_REF] Haber | Image registration with guaranteed displacement regularity[END_REF][START_REF] Sdika | A fast nonrigid image registration with constraints on the jacobian using large scale constrained optimization[END_REF]. However, such methods are unsuited for further statistical analysis, as statistical procedures on resulting transformation parameters (such the empirical Euclidean mean), does not lead to invertible transforms. Moreover, the inverse transforms do not belong to the initial space of transformations.

More recently, diffeomorphisms have been used to model deformations between images in the context of Grenander's pattern theory [START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF][START_REF] Hernandez | Registration of anatomical images using paths of diffeomorphisms parameterized with stationary vector field flows[END_REF][START_REF] Ashburner | A fast diffeomorphic image registration algorithm[END_REF]. In this framework, the set of admissible diffeomorphic deformations is considered as a Riemmaniann manifold, and thus first and second order statistics analysis on manifolds [START_REF] Fletcher | Principal geodesic analysis for the study of nonlinear statistics of shape[END_REF][START_REF] Pennec | Intrinsic statistics on Riemannian manifolds: basic tools for geometric measurements[END_REF] can be applied to perform statistical analysis of diffeomorphic deformations [WBR + 07, ACPA06, BHO07, HBO07]. Such diffeomorphisms are constructed as solution of a differential equation governed by a time dependent vector field [START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF]. The optimal diffeomorphism after the registration of two images can be fully characterized by the initial point in time (or equivalently by the initial momentum) of the associated time dependent vector field. This important momentum conservation property allows to carry out PCA on the Hilbert space of initial momentums [WBR + 07].

A particular sub-class of diffeomorphic deformations is of the set of diffeomorphisms generated by an O.D.E. governed by stationary vector fields. In this way, diffeomorphisms are directly characterized by vector fields in a Hilbert space and statistical analysis like PCA can be carried on these vector fields. Compared to the case of diffeomorphisms generated by nonstationary vector fields, the computational cost of the registration is considerably smaller, while keeping comparable accuracy according to experimental results reported in [START_REF] Hernandez | Registration of anatomical images using paths of diffeomorphisms parameterized with stationary vector field flows[END_REF][START_REF] Ashburner | A fast diffeomorphic image registration algorithm[END_REF]. Moreover, diffeomorphisms generated by stationary vector fields have some desirable properties that allow simple and fast computations of exponential and logarithms, that is, computing the diffeomorphism associated to a vector field and vice versa. Hence, PCA methods for manifolds can be implemented to analyze shape variability of diffeomorphic deformations generated by stationary vector fields, see [ACPA06, BHO07, HBO07].

Consistency of geometric PCA in statistical deformable models

In the last decades, in the framework of Grenander's pattern theory, there has been a growing interest in the use of first order statistics for the computation of mean pattern from a set of images [AAT07, AKT10, BC11, BGV09, BGL09] and in the construction of consistent procedures. However, there is not so much work in the statistical literature on the consistency of second order statistics for the analysis of geometric variability of images. In particular, the convergence of such procedures in simple statistical models has generally not been established. We study the consistency of geometric PCA, in the context of the following statistical deformable model

Y i (x) = f * (ϕ -1 (V i , x)) + W i (x), x ∈ Ω, i = 1, . . . , n, (3.1)
where

• f * is an unknown mean pattern belonging to L 2 (Ω),

• ϕ is a deformation operator associated to a Hilbert space V, (in the sense of Definition 2.1), equipped with inner product •, • and induced norm • ,

• V 1 , . . . , V n are independents copies of V , a zero-mean, square-integrable V-valued random variable (i.e. EV = 0 and E V 2 < ∞),

• There exists µ > 0 (regularization parameter) such that P(V ∈ V µ ) = 1,

• > 0 is a noise level parameter,

• W 1 , . . . , W n are independents copies of a zero mean Gaussian process W ∈ L 2 (Ω), such that E W 2 2 = 1,

• (V 1 , . . . , V n ) and (W 1 , . . . , W n ) are mutually independent.
Additionally, we assume that the eigenvalues κ 1 ≥ κ 2 ≥ . . . ≥ 0 of the population covariance operator K, defined as

Kv = E V, v V , v ∈ V, (3.2)
have algebraic multiplicity 1, i.e. κ 1 > κ 2 > . . . ≥ 0. This implies that the λ-th eigen-gap, defined as δ λ := min λ ∈Λ\{λ} |κ λ -κ λ |, is strictly positive, for any λ ∈ Λ.

Observe that the function f * in (3.1) models the common shape of the Y i 's. The W i 's represent the individual variations in intensity of the data around the mean pattern f * , and thus correspond to a classical source of variability that could be analyzed by standard PCA. To the contrary, the random elements ϕ -1 (V i , •) model deformations of the domain Ω, and thus correspond to a source of shape variability in the data. Definition 3.1 (Population geometric modes of variations). Let K be the population covariance operator, defined in (3.2), with (population) eigenvalues κ 1 > κ 2 > . . . ≥ 0 and (population) orthonormal eigenvectors φ 1 , φ 2 , . . .. For λ ∈ Λ, the λ-th population mode of geometric variation of the random variable V is the homeomorphism ψ λ : Ω → Ω defined by

ψ λ (x) = ϕ -1 ( √ κ λ φ λ , x), x ∈ Ω.
In this paper we say that geometric PCA is a consistent procedure if, for data Y = (Y 1 , . . . , Y n ) following model (3.1), and for all λ ∈ Λ, the λ-th empirical mode of geometric variation ψλ (see equation (2.9)) tends to the λ-th population mode of geometric variation ψ λ , as n → +∞ and → 0, in a sense to be made precise later on. In this context, the empirical modes of geometric variation are obtained from the eigenvalues κλ and the eigenvectors φλ of the empirical covariance operator

Kn v = 1 n n i=1 Vi -V n , v ( Vi -V n ), for v ∈ V, (3.3)
where ( V1 , . . . , Vn ) belongs to arg min v∈Vµ M t (v, Y , f * ) or arg min v∈U M g (v, Y ). Consequently, from now on, empirical eigenvalues, eigenvectors and modes of geometric variations will be considered as random elements.

Remark 3.1. Observe that, for template registration, where v ∈ arg min v∈Vµ M t (v, y, f * ), each coordinate vi depends only on v i . This fact implies that V1 , . . . , Vn are i.i.d. However, this is not the case for groupwise registration, where vi may depend on all the v i 's.

The asymptotic n → +∞ is rather natural and correspond to the setting of a growing number of images. On the other hand, the setting → 0 corresponds to the analysis of the influence of the additive term W i in model (3.1). The main result of this section is that geometric PCA is consistent only in the double asymptotic setting n → +∞ and → 0. This result illustrates the fact that the photometric perturbations W i , i = 1, . . . , n in model (3.1) have to be sufficiently small in order to recover the geometric modes of variation. Definition 3.2. A deformation operator ϕ (see Definition 2.1) is said to be µ-regular if there exists µ > 0 such that

(i) Ω f 2 (ϕ -1 (v, x))dx ≤ A µ Ω f 2 (x)dx, (3.4)
for all f ∈ L 2 (Ω), v ∈ V µ and some constant A µ > 0; and

(ii) the mapping v → ϕ(v, •) from V µ to C(Ω, Ω) is continuous,
where C(Ω, Ω) is the space of continuous function from Ω to Ω, endowed with the metric d C (ψ, φ) := sup x∈Ω d Ω (ψ(x), φ(x)).

Note that if ϕ(v, •) is sufficiently smooth such that the determinant of its Jacobian matrix is bounded, that is |det(J(ϕ(v, x))| ≤ A µ for all v ∈ V µ and x ∈ Ω, then (3.4) follows from a change of variable.

Finally, before stating our consistency results, we define convergence in probability in the double asymptotic setting n → ∞, → 0. Let X n, , X n , X , X, n = 1, 2, . . . , > 0 random variables with values on a metric space (S, d). The notation plim X n, = X n stands for

d(X n, , X n ) → 0 in probability as → 0; plim n X n, = X denotes d(X n, , X ) → 0 in probabil- ity as n → ∞. Finally, plim n, X n, = X is equivalent to plim n plim X n, = plim plim n X n, = X.
In this paper, all equalities and inequalities involving random variables are understood in the almost sure sense. We require the following definition, in our main results: for u, v ∈ V,

sin(u, v) := 1 -u/ u , v/ v 2 .

Case of template registration

* : V µ → L 2 (Ω), defined by ϕ * (v) := f * • ϕ -1 (v, •) for v ∈ V µ , is injective and its inverse ϕ * -1 : ϕ * (V µ ) → V µ is continuous. Then plim n, κλ = κ λ and plim n,
sin 2 ( φλ , φ λ ) = 0, for all λ ∈ Λ.

Remark 3.2. Observe that, under the hypothesis of Theorem 3.1, the λ-th empirical mode of geometric variation ψλ converges the population mode of geometric variation ψ λ in probability, when n → +∞ and → 0, as elements of (C(Ω, Ω), d C ). Indeed, this result follows from the continuity of the mapping v → ϕ(v, •), which is guaranteed by the µ-regularity of ϕ.

We show below how a stronger regularity assumption on ϕ, allows one to obtain rates of convergence for κλ and φλ , via a concentration inequality that depends explicitly on n and .

Theorem 3.2. Under the hypothesis of Theorem 3.1 and if ϕ * -1 is uniformly Lipschitz (in the sense that u -v 2 ≤ L(f * , µ) ϕ * (u) -ϕ * (v) 2 2 , for every u, v ∈ V µ and some constant L(f * , µ) > 0 depending only on f * and µ), then

P |κ λ -κ λ | 2 > C(f * , µ) max(h(u, n, ) + h(u, n, ); g(u, n)) ≤ exp(-u),
for any u > 0, where C(f * , µ) > 0 is a constant depending only on f * and µ; h(u, n, ) =

2 1 + 2 u n + 2 u n and g(u, n) = u n + u 2 n 2 + u n 2 .
Take now u * > 0 such that

C(f * , µ) max(h(u * , n, ) + h(u * , n, ); g(u * , n)) < (δ λ /2) 2 , then for any 0 < u ≤ u * P sin 2 ( φλ , φ λ ) > (2/δ λ ) 2 C(f * , µ) max(h(u, n, ) + h(u, n, ); g(u, n)) ≤ 2 exp(-u).

Case of groupwise registration

In order to prove consistency for groupwise registration, we require model (3.1) to satisfy the set of identifiability assumptions, shown below. For u, v ∈ V = V n , let

D g (u, v) := M g (u, (f * 1 , . . . , f * n )) = 1 n n i=1 Ω   f * i (ϕ(u i , x)) - 1 n n j=1 f * j (ϕ(u j , x))   2 dx, (3.5)
where

f * i (x) := f * (ϕ -1 (v i , x)), x ∈ Ω, i = 1, . . . , n. Observe that, D g (u, V ) = M g (u, Y ) when Y = (Y 1 , . . . , Y n ) follows model (3.1) with = 0.
Definition 3.3 (g-identifiability). Model (3.1) is said to be g-identifiable if (i) there exists a measurable function u * : V µ → U such that for every η > 0 there exists a constant C > 0 not depending on n, with D g (u, v) -D g (u * , v) > C, for every u ∈ U satisfying d2 (u * , u) > η, and

(ii) plim n d2 (u * (V ), V ) = 0, where d2 (u, v) := 1 n n i=1 u i -v i 2 , for u, v ∈ V.
Observe that condition (i) above implies that, for every v ∈ V µ , D g (u, v) has a unique measurable minimizer u * (v) on U .

Theorem 3.3. Let Y = (Y 1 , . . . , Y n ) i.i.d. observations of model (3.1)
, with deformation operator ϕ and regularization parameter µ. Let κλ and φλ , λ ∈ Λ, be the empirical eigenvalues and eigenvectors corresponding to groupwise registration of Y . Suppose that ϕ is µ-regular and that (3.1) is g-identifiable. Then plim n, κλ = κ λ and plim n, sin 2 ( φλ , φ λ ) = 0, for all λ ∈ Λ.

Remark 3.3. Observe that, as in the case of template registration, it can be shown that under the hypothesis of Theorem 3.3, ψλ converges to ψ λ in probability, as n → +∞ and → 0.

Translation operators

We study the applicability of Theorems 3.1, 3.2 and 3.3 to translation operator ϕ given by (2.1). In this case, Ω = [0, 1) d , for some integer d ≥ 1, is equipped with the distance d Ω (x, y) :=

d k=1 min{|x k -y k |, 1 -|x k -y k |}, for x = (x 1 , . . . , x d ) ∈ Ω. Let also V = R d ,
be equipped with the usual Euclidean inner product. Now, let us show that ϕ is a deformation operator in the sense of Definition 2.1. It holds that ϕ(0, •) is the identity in Ω and ϕ -1 (v, •) = ϕ(-v, •). Last, from (i) in Lemma C.1, it follows that the mapping ϕ(v, •) is continuous for all v ∈ V. Moreover, we prove that ϕ is µ-regular, for all µ > 0: for (i) in Definition 3.2, take f ∈ L 2 (Ω) and consider its periodic extension

f per to R d . Then (3.4) is a consequence of f (ϕ(v, x)) = f per (x + v) which holds for all v ∈ R d , x ∈ Ω.
Finally, condition (ii) in Definition 3.2, follows from (ii) in Lemma C.1.

We impose further conditions on model (3.1) implying that ϕ * -1 , defined in Theorem 3.1 is Lipschitz. Let θ k , k = 1, . . . , d, be the low frequency Fourier coefficients of the template f * , that is

θ k := Ω f * (x)e -i2πx k dx = 0 for all 1 ≤ k ≤ d. (3.6)
Lemma 3.1. Suppose that f * is such that θ k = 0 for all 1 ≤ k ≤ d and let µ < 1/2, then ϕ * is injective and ϕ * -1 is uniformly Lipschitz.

Hence, if θ k = 0 for all 1 ≤ k ≤ d and µ < 1/2, the hypotheses of Theorems 3.1 and 3.2 are verified. Thus the geometric PCA is consistent, in the case of template registration, with translation operator. Observe that hypotheses of Lemma 3.1 imply that translation invariant templates f * are excluded.

We now turn our attention to groupwise registration. We have to impose further conditions on model (3.1) ensuring g-identifiability, so that Theorem 3.3 applies. The set of deformations parameters U ⊂ V µ over which M g (v, y) will be minimized, is U = U 0 , given in (2.2). We have the following. Proposition 3.1. Suppose θ k = 0 for all 1 ≤ k ≤ d, and that

P(V ∈ [-ρ, ρ] d ) = 1 with ρ = min( µ 2 , µ √ d ) and 0 < µ < 1 12 . Then D g (u, v) -D g (u * (v), v) ≥ C(f * , µ) d2 (u, u * (v)), for all u ∈ U , (3.7)
where u * (v

) := v 1 -1 n n i=1 v i , . . . , v n -1 n n i=1 v n and C(f * , µ) > 0 is a constant depending only on f * and µ.
Remark that in Proposition 3.1, D g (u * (v), v) = 0. This shows that D g (u, v) is bounded below by a quadratic functional.

We are now ready to prove g-identifiability under the hypotheses of the previous proposition. Observe that (i) in Definition 3.3 follows at once from (3.7). For (ii) note that d2 (u * (V ), V ) = 1 n n i=1 V i . Hence, given that EV = 0, from Bernstein's inequality for bounded random variables in a Hilbert space (see e.g. [Bos00], Theorem 2.6) we conclude that, for any η > 0,

P (d(u * (V ), V ) > η) ≤ 2 exp - nη 2 2E V 2 + µ 3 η .
Therefore d2 (u * (V ), V ) converges in probability to 0 as n → +∞. Finally, having checked the g-identifiability of the model, we conclude that the geometric PCA is consistent, in the case of groupwise registration, with translation operator.

Conclusion and discussion

The contribution of this paper is two folds. First, the use of deformation operators (as introduced in this paper) provides a general framework for modelling and analyzing the shape variability of a set of images. As a particular case, it allows the use of diffeomorphic deformations parametrized by stationary vector fields. In the case of diffeomorphims computed with nonstationary vector fields as in [START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF], the link with our framework is not straightforward. Indeed, in this setting, there are two possibilities for defining the deformation operators. One can parameterizing them either by the Hilbert space of time-dependent vector fields, or by the Hilbert space of initial velocities. Both cases are rather complex from the analytical and the computational point of views and treating them is beyond the scope of this paper. In contrast, due to its analytical and numerical tractability, we have preferred to focus on diffeomorphic deformation operators parametrized by stationary vector fields.

The second contribution of this paper is the study of the consistency of geometric PCA in statistical deformable models which, to the best of our knowledge, has not been investigated before. We also hope that the methods presented in this paper will stimulate further investigation into the development of consistent statistical procedures for the analysis of shape variability.

A Preliminary technical results

In this sub-section, we give a deviation in probability of sup v∈Vµ |D g (v) -M g (v, Y )| under appropriate assumptions on the deformation operators and the additive noise in model (3.1), where M g (v, Y ) and D g (v) are defined in (2.4) and (3.5) respectively.

Lemma A.1. Consider model (3.1), with µ-regular deformation operator ϕ. Let

Q(v) = 2 n n i=1 Ω W 2 i (ϕ(v i , x))dx, v ∈ V µ . (A.1)
Then, for any s > 0,

P sup v∈Vµ Q(v) ≥ A µ h(s, n, ) ≤ exp(-s),
where A µ is given in Definition 3.2 (i) and h(s, n, ) = 2 1 + 2 s n + 2 s n .

Proof. From Definition 3.2 (i), we have, for v ∈ V µ ,

Q(v) ≤ A µ 2 n n i=1 Ω W 2 i (x)dx. (A.2)
Let g ∈ L 2 (Ω) and K W g(x) = Ω k(x, y)g(y)dy be the covariance operator of the random process W , where k(x, y) = EW (x)W (y) for x, y ∈ Ω. Then, there exist orthonormal eigenfunctions (φ λ ) λ∈Λ in L 2 (Ω) and positive eigenvalues (w λ ) λ∈Λ , such that K W φ λ = w λ φ λ , with w 1 ≥ w 2 ≥ . . . ≥ 0 and Λ = {1, 2, . . .}. For any 1 ≤ i ≤ n, the Gaussian process W i can thus be decomposed as

W i = λ∈Λ w 1/2 λ ξ i,λ φ λ , 20
where ξ i,λ = w -1/2 λ W i , φ λ 2 is a Gaussian variable with zero mean and variance 1, such that Eξ i,λ ξ i,λ = 0 for λ = λ . Therefore, W i 2 2 = +∞ λ=1 w λ ξ 2 i,λ where ξ i,k , i = 1 . . . , n, k ≥ 1 are i.i.d. standard Gaussian random variables. We have, from the assumptions on W , E W i 2 2 = +∞ λ=1 w λ = 1 < +∞, and one can thus consider the following centered random variable

Z = n i=1 +∞ λ=1 w λ (ξ 2 i,λ -1).
Let 0 < t < (2w 1 ) -1 . Since the generating function of a χ 2 random variable, with one degree of freedom, is E e sξ 2 i,λ

= (1 -2s) -1/2 , for s > 0, it follows that log E e tZ = -n

+∞ λ=1 tw λ + 1 2 log(1 -2tw λ ) . (A.3)
Then, using the inequality -s -1 2 log(1 -2s) ≤ s 2 1-2s , which holds for all 0 < s < 1 2 , from (A.3) we obtain

log E e tZ ≤ n +∞ λ=1 t 2 w 2 λ 1 -2tw λ ≤ t 2 n 1 -2tw 1 +∞ λ=1 w λ 2 = t 2 n 1 -2tw 1 < ∞.
Arguing e.g. as in [START_REF] Birgé | Minimum contrast estimators on sieves: exponential bounds and rates of convergence[END_REF], the above inequality implies that, for any s > 0,

P Z > 2w 1 s + 2 √ ns ≤ exp(-s). (A.4) By (A.2), it follows that sup v∈Vµ Q(v) ≤ A µ 2 n n i=1 Ω W 2 i (x)dx = A µ 2 n n i=1 +∞ λ=1 w λ ξ 2 i,λ = A µ 2 n (Z + n).
Hence, it follows from (A.4) that

P sup v∈Vµ Q(v) > A µ 2 n n + 2w 1 s + 2 √ ns ≤ exp(-s),
for any s > 0. The conclusion follows noting that w 1 ≤ E W 2 2 = 1.

Lemma A.2. Consider model (3.1), with µ-regular deformation operator ϕ. Let

D t (v) := 1 n n i=1 Ω (f * i (ϕ(v i , x)) -f * (x)) 2 dx, v ∈ V, (A.5) with f * i (x) := f * (ϕ -1 (V i , x)), x ∈ Ω, i = 1, . . . , n. Then P sup v∈Vµ |D t (v) -M t (v, Y , f * )| > C h(s, n, ) + h(s, n, ) ≤ exp(-s), s > 0,
where M t is defined in (2.3), C > 0 is a constant depending only on f * and µ, and h(s, n, ) is defined in Lemma A.1.

Proof.

For v ∈ V = V n , let R(v) = 2 1 n n i=1 Ω (f * i (ϕ(v i , x)) -f * (x))W * i (ϕ(v i , x))dx.
For any v ∈ V µ , we have the decomposition

M t (v, Y , f * ) = D t (v) + Q(v) + R(v), (A.6)
where Q is defined in (A.1). By applying the Cauchy-Schwartz inequality in L 2 (Ω) and in

R n we obtain R(v) ≤ 2 D t (v) Q(v). Also, the µ-regularity of ϕ implies D t (v) ≤ 4A 2 µ f * 2 2 , and therefore R(v) ≤ 4A µ f * 2 Q(v). Now, using the decomposition (A.6), one obtains sup v∈Vµ |D t (v) -M t (v, Y , f * )| ≤ max (1, 4A µ f * 2 ) sup v∈Vµ Q(v) + sup v∈Vµ Q(v)
and the result follows from Lemma A.1.

Lemma A.3. Consider model (3.1), with µ-regular deformation operator ϕ. Then

P sup v∈Vµ |D g (v, V ) -M g (v, Y )| > C h(s, n, ) + h(s, n, ) ≤ exp(-s), s > 0,
where M g and D g are defined in (2.4) and (3.5) respectively; C > 0 is a constant, depending only on f * and µ, and h(s, n, ) is defined in Lemma A.1.

Proof. Let

Q g (v) = 2 n n i=1 Ω   W i (ϕ(v i , x)) - 1 n n j=1 W j (ϕ(v j , x))   2 dx, v ∈ V and R(v) = 2 1 n n i=1 Ω   1 n n j=1 f * j (ϕ(v j , x)) -f * i (ϕ(v i , x))   ×   1 n n j=1 W * j (ϕ(v j , x)) -W * i (ϕ(v i , x))   dx, v ∈ V.
Then, for any v ∈ V µ , we have the decomposition

M g (v, Y ) = D g (v) + Q g (v) + R(v). (A.7)
From the Cauchy-Schwartz inequality in L 2 (Ω) and in R n , we have R(v) ≤ 2 D g (v) Q g (v). Also, from the µ-regularity of ϕ, we obtain

D g (v) ≤ 1 n n i=1 Ω f * i (ϕ -1 (v i , x)) 2 dx ≤ A 2 µ f * 2 2 . So R(v) ≤ 2A µ f * 2 Q g (v)
. Now, using the decomposition (A.7), one obtains

sup v∈Vµ |D g (v) -M g (v, Y )| ≤ max (1, 2A µ f * 2 ) sup v∈Vµ Q g (v) + sup v∈Vµ Q g (v)
and the result follows from the fact that Q g ≤ Q (see (A.1)) and Lemma A.1.

Remark A.1. Observe that Lemma A.2 and Lemma A.3 imply

plim sup v∈Vµ |D t (v) -M t (v, Y , f * )| = plim sup v∈Vµ |D g (v, V ) -M g (v, Y )| = 0, (A.8) and plim n, sup v∈Vµ |D t (v) -M t (v, Y , f * )| = plim n, sup v∈Vµ |D g (v, V ) -M g (v, Y )| = 0. (A.9)
The proofs of Theorems 3.1, 3.2 and 3.3 rely on the following two propositions that establish the consistency of the registration procedures described in Section 2.2.

P d2 ( V , V ) > C h(s, n, ) + h(s, n, ) ≤ exp(-s), s > 0,
where C > 0 is a constant, depending only on f * , µ and h(s, n, ) is defined in Lemma A.1.

Proof. Observe that D t (V ) = 0 and so

D t (v) = D t (v) -D t (V ) ≤ 2 sup v∈Vµ |M t (v, Y , f * ) -D t (v)|, v ∈ V µ , (A.10)
where D t and M t are defined in (A.5) and (2.3) respectively. On the other hand, from Definition 3.2

(i), 1 n n i=1 ϕ * (V i ) -ϕ * (v i ) 2 2 ≤ A µ D t (v), v ∈ V µ .
Hence,

1 n n i=1 ϕ * (V i ) -ϕ * ( Vi ) 2 2 ≤ 2A µ sup v∈Vµ |M t (v, Y , f * ) -D t (v)|. (A.11)
We proceed now to prove part (i). From (A.11) and (A.8) we have plim ϕ * (V i ) -ϕ * ( Vi ) 2 2 = 0, that is, plim ϕ * ( Vi ) = ϕ * (V i ) for i = 1, . . . , n. From the continuity of ϕ * -1 we have plim Vi = V i , for i = 1, . . . , n, therefore

plim n plim 1 n n i=1 V i -Vi 2 = 0. Now, the fact that V 1 -V1
2 is bounded by 2µ and tends to 0 in probability as → 0, implies that E V 1 -V1 2 → 0 as → 0. Noting that (V i -Vi ) i≥1 are i. Proof. Let u * (V ) be the unique minimizer of D g (u, V ) on U , which exists because the model is g-identifiable. Since (by definition) V ∈ arg min v∈U M g (v, Y ), one obtains that

D g ( V , V ) -D g (u * , V ) ≤ 2 sup u∈U |M g (u, Y ) -D g (u, V )| ≤ 2 sup u∈Vµ |M g (u, Y ) -D g (u, V )|.
Therefore, from (A.9) and the g-identifiability of the model, we have plim n, d2 ( V , u * ) = 0. Also, the g-identifiability implies that plim n, d2 (u * , V ) = 0. Finally, the conclusion follows from the inequality d2 ( V , V ) ≤ 2 d2 ( V , u * ) + 2 d2 (u * , V ).

From the above inequality, combined with B.4 and B.5, we obtain plim

B.3 Proof of Theorem 3.3

Proof. We proceed similarly as in the proof of Theorem 3.1. In the case of groupwise registration, inequalities (B.1), (B.2) and (B.3) are still valid, and can be combined with Proposition A.2 to obtain plim n, Kn = K. The rest of proof is identical to that of Theorem 3.1. .

C Technical results for translation operators

Let 0 ≤ α < 1/4. Using a second order Taylor expansion and the mean value theorem, one has that cos(2πu) ≤ 1 -C(α)|u| 2 for any real u such that |u| ≤ α, with C(α) = 2π 2 cos(2πα). Under the hypothesis of the proposition, one has that u

(k) m -V (k) m -u (k) m -V (k) m
≤ 2(µ+ρ) < 1/4. Therefore, for α = 2(µ + ρ), it follows that

1 n n m=1 e i2π u (k) m -V (k) m 2 ≤ 1 n + 2 n 2 n-1 m=1 n m =m+1 1 -C(α) u (k) m -V (k) m -u (k) m -V (k) m 2 ≤ 1 - 2 n 2 n-1 m=1 n m =m+1 C(α) u (k) m -V (k) m -u (k) m -V (k) m 2 .

  (a) and Figure 2 (a), and standard PCA does not meaningfully reflect the modes of variability of the data, see Figure 1 (b),(c) and Figure 2 (b),(c).

Figure 1 :

 1 Figure 1: Digit 1. Two first rows: a sample of 8 images out of n = 100 taken from the Mnist data base [LBBH98] -Last row: standard Euclidean mean, first and second empirical modes of variation.

Figure 2 :

 2 Figure 2: Digit 5. Two first rows: a sample of 8 images out of n = 100 taken from the Mnist data base [LBBH98] -Last row: standard Euclidean mean, first and second empirical modes of variation.

Figure 3 :

 3 Figure 3: Visualization of the geometric PCA and standard PCA results for digits 1 to 3. For each digit, the superior row from left to right corresponds to ȳ, ȳn -2 √ γ1 û1 , ȳn + 2 √ γ1 û1 , ȳn -2 √ γ2 û2 and ȳn + 2 √ γ2 û2 respectively. The inferior row from left to right corresponds to f , f • ψ1,-2 , f • ψ1,2 , f • ψ2,-2 and f • ψ2,2 . We observe, in several cases, that standard PCA results (first row of each digit) do not recover well shapes of the digits and produce a blurring of the images. In contrast, geometric PCA results (second row of each digit) recover the geometric features of the digits.

Figure 4 :

 4 Figure4: Visualization of the geometric PCA and standard PCA results for digits 4 to 6. For each digit, the superior row from left to right corresponds to ȳ, ȳn -2 √ γ1 û1 , ȳn + 2 √ γ1 û1 , ȳn -2 √ γ2 û2 and ȳn + 2 √ γ2 û2 respectively. The inferior row from left to right corresponds to f , f • ψ1,-2 , f • ψ1,2 , f • ψ2,-2 and f • ψ2,2 . We observe, in several cases, that standard PCA results (first row of each digit) do not recover well shapes of the digits and produce a blurring of the images. In contrast, geometric PCA results (second row of each digit) recover the geometric features of the digits.

Figure 5 :

 5 Figure5: Visualization of the geometric PCA and standard PCA results for digits 7 to 9. For each digit, the superior row from left to right corresponds to ȳ, ȳn -2 √ γ1 û1 , ȳn + 2 √ γ1 û1 , ȳn -2 √ γ2 û2 and ȳn + 2 √ γ2 û2 respectively. The inferior row from left to right corresponds to f , f • ψ1,-2 , f • ψ1,2 , f • ψ2,-2 and f • ψ2,2 . We observe, in several cases, that standard PCA results (first row of each digit) do not recover well shapes of the digits and produce a blurring of the images. In contrast, geometric PCA results (second row of each digit) recover the geometric features of the digits.

Theorem 3. 1 .

 1 Let Y = (Y 1 , . . . , Y n ) i.i.d. observations of model (3.1), with deformation operator ϕ and regularization parameter µ. Let κλ and φλ , λ ∈ Λ, be the empirical eigenvalues and eigenvectors corresponding to template registration of Y , with f = f * . Suppose that ϕ is µ-regular and the mapping ϕ

Proposition A. 1 .

 1 Let V ∈ arg min v∈Vµ M t (v, Y , f * ) the parameters obtained from template registration of Y on f * . Then, (i) under the hypotheses of Theorem 3.1, plim n, d2 ( V , V ) = 0, and (ii) under the hypotheses of Theorem 3.2,

  i.d. (see Remark 3.1), we conclude from the weak law of large number that plim plimn 1 n n i=1 V i -Vi 2 = 0 thus proving part (i).For (ii), note that inequality A.11 and that ϕ * -1 is uniformly Lipschitz, with constantL(f * , µ) > 0, imply d2 (u, v) ≤ 2A µ L(f * , µ) sup v∈V |M t (v, Y , f * ) -D t (v)|and the result follows from Lemma A.2.Proposition A.2. Let V ∈ arg min v∈U M g (v, Y ) the parameters obtained from groupwise registration of Y . Then, under the hypotheses of Theorem 3.3, plim n, d2 ( V , V ) = 0.

Lemma C. 1 .

 1 Let ϕ be defined by (2.1), then (i) d Ω (ϕ(v, x), ϕ(v, y)) = d Ω (x, y) for all x, y ∈ Ω and v ∈ V.Then, by the mean value theorem, we have |e-i2πu k -e -i2πv k | 2 = | cos(2πu k ) -cos(2πv k )| 2 + | sin(2πu k ) -sin(2πv k )| 2 ≥ (2π) 2 cos 2 (ρ)|u k -v k | 2 , for any 0 ≤ u k , v k ≤ ρ. Hence, f * (ϕ -1 (u, x)) -f * (ϕ -1 (v, x)) 2 2 ≥ (2π) 2 cos 2 (ρ) min 0≤k≤d |θ (k) | 2 d k=1 |u k -v k | 2 .C.2 Proof of Proposition 3.1Proof. Remark that D g , defined in (3.5), has the following expression in the Fourier domain:D g (u, i2π ,V j -u j -θ e -i2π ,Vm-um 2   , u ∈ V. (C.3)For u ∈ U 0 we haveD g (u, V )

  sin 2 ( φλ , φ λ ) = 0.B.2 Proof of Theorem 3.2Proof. Combining (B.1), (B.2) and (B.3) with Proposition A.1 (ii), we obtainP Kn -K 2 HS > C max(h(s, n, ) + h(s, n, ); g(s, n)) ≤ 2 exp(-s), s > 0,where C > 0 is a constant depending only on f * and µ and g(s, n) = s n + s 2 n 2 + s P |κ λ -κ λ | 2 > C max(h(s, n, ) + h(s, n, ); g(s, n)) ≤ 2 exp(-s), s > 0.Take now s * > 0 such that C max(h(s * , n, ) + h(s * , n, ); g(s * , n)) < (δ λ /2) 2 . Then, thanks to (B.5) and (A.13) we obtain, for any 0 < s ≤ s

	2
	n	. Hence,
	from (A.12) we obtain

n, * , 1 -2 exp(-s) ≤ P Kn -K 2 HS < C max(h(s, n, ) + h(s, n, ); g(s, n)) = P Kn -K 2 HS < C max(h(s, n, ) + h(s, n, ); g(s, n)), δλ > δ λ /2 ≤ P (1/ δλ ) 2 Kn -K 2 HS < (2/δ λ ) 2 C max(h(s, n, ) +

h(s, n, ); g(s, n))

≤ P sin 2 ( φλ , φ λ ) < (2/δ λ ) 2 C max(h(s, n, ) + h(s, n, ); g(s, n)) .

Acknowledgments

The authors acknowledge the support of the French Agence Nationale de la Recherche (ANR) under reference ANR-JCJC-SIMI1 DEMOS, and the financial support from the French embassy

In what follows, HS denotes the Hilbert-Schmidt norm of operators on a Hilbert space H. Recall that, given an orthonormal basis {e j } j≥1 of H, the Hilbert-Schmidt norm of an operator K is defined as K 2 HS = j,k K(e j ), e k 2 .

Lemma A.4. Let H be a separable Hilbert space with inner product •, • and induced norm • . Let {u i } n i=1 , {v i } n i=1 in B r = {h ∈ H : h ≤ r} for some r > 0. Define the covariance operators

Proof.

where L * denotes the adjoint the operator of

Then, after some simple calculations, we get

Hence, since u i ∈ B r , i = 1, . . . , n, one has

Similarly, since i ≤ u i + v i ≤ 2r,

, which completes the proof.

The following theorem follows from the theory developed in [START_REF] Bhatia | Perturbation of spectral subspaces and solution of linear operator equations[END_REF][START_REF] Davis | The rotation of eigenvectors by a perturbation[END_REF].

Theorem A.1. Let H be a separable Hilbert space endowed with inner product •, • H . Let A, Â : H → H be self-adjoint Hilbert-Schmidt operators on H, with eigenvalues/eigenvectors pairs (κ λ , φ λ ) λ≥1 and (κ λ , φλ ) λ≥1 respectively. Then,

B Proofs of main results

B.1 Proof of Theorem 3.1

Proof. Let Kn be the sample covariance operator

The first term in the right hand side of inequality (B.1) can be controlled by using Lemma A.4 and noting that

Let us now bound the second term in the right hand side of (B.1). To do so, remark that V , V i ≤ µ, i = 1, . . . , n, and, thanks to a Bernstein's inequality for Hilbert-Schmidt operators (see e.g. [START_REF] Bosq | Nonparametric statistics for stochastic processes[END_REF], Chapter 3), it follows that From (A.12), it holds that δ λ ≤ δλ + max λ |κ λ -κλ | ≤ δλ + Kn -K HS . By (B.4), it follows that lim n,

Recalling that, from the specification of model (3.1), we have δ λ > 0, hence inequality (A.13) implies

Remark that, for any a ∈ R, there exists a unique k(a) ∈ Z such that mod(a, 1) = a+k(a). Then mod(a, 1) -mod(b, 1) = a -b + k(a) -k(b).

Take a, b ∈ R such that |a -b| < 1 and assume that a ≥ b. Since a -b ∈ [0, 1) and mod(a, 1)mod(b, 1) ∈ [-1, 1], we obtain that

We conclude that for a ≥ b min{|mod(a, 1) -mod(b, 1)|, 1 -|mod(a, 1) -mod(b, 1)|} = min{|b -a|, 1 -|b -a|}. (C.1) Because of the symmetry in the expression above, we conclude that (C.1) is valid for any a, b ∈ R, such that |a -b| < 1.

For the sake of simplicity, let us prove the lemma in the one-dimensional case (i.e. d = 1), where d Ω (x, y) := min{|x -y|, 1 -|x -y|}. Take x, y ∈ Ω and u, v ∈ V. Part (i) is directly implied by (C.1), taking a := x + v and b := y + v. For part (ii), note that C.1 Proof of Lemma 3.1

Hence, using the lower bound (C.4), it follows that, for u ∈ U 0 ,

The following identity is obtained from elementary algebraic manipulations and the fact that