Je Remercie

Thomas Camille D'

Je suis très heureuse d'avoir ici l'occasion d'exprimer ma gratitude à ceux qui m'ont aidée et soutenue pendant ces quatre années de thèse.

Je suis profondément reconnaissante à mes directeurs Marc et Frédéric avec qui j'ai été honorée de travailler. Je remercie Marc pour m'avoir guidée sur les chemins de la géométrie riemannienne avec une patience et une gentillesse à nulles autres pareilles, et pour avoir toujours répondu par la négative à la question « est-ce que je te dérange ? ». Je remercie Frédéric pour son enthousiasme communicatif et ses encouragements, ses idées dont il n'est jamais avare, et pour m'avoir proposé ce beau sujet de thèse.

Je tiens à remercier Silvère Bonnabel, Eric Klassen et Xavier Pennec pour avoir accepté de rapporter cette thèse, ainsi que Stéphanie Allassonnière, Jérémie Bigot et Barbara Tumpach pour avoir bien voulu faire partie du jury. Je suis reconnaissante à Barbara pour son invitation à Lille et je la remercie d'avoir accepté de faire un long déplacement pour la soutenance.

Je remercie très chaleureusement Christian Léonard, qui m'a donné le premier le goût des probabilités, pour son soutien moral et mathématique et tout le temps qu'il a bien voulu me consacrer.

A few words in English to thank Martin Bauer and Martins Bruveris for allowing me to meet the "shape space" community in Vienna during the Infinite-dimensional Riemannian Geometry workshop, and for making it into a stimulating and friendly experience. I am particularly grateful to Martin for his warm welcome at this occasion, for his support during these past few years, and for sharing with me his wisdom in geometry as well as other subjects. I also thank Sergey Kushnarev for an exceptional culinary tour of Singapore during the Mathematics of Shapes and Applications workshop. I have a thought for Jakob who shared both experiences with me, as well as for Barbara, Pierre and Line. It wouldn't have been the same without them.

Je repasse en français pour remercier tous les doctorants avec qui j'ai eu la chance de partager cette expérience particulière qu'est la thèse. Je remercie mes voisins de bureau à Jussieu, anciens et nouveaux : Eric, Saad, i ii Alexandra, Omar, Yi, Candia, Nina et Pamela, ainsi que Olga qui n'était jamais très loin pour nous rappeler d'aller au GTT. Je remercie tout particulièrement Éric pour sa bonne humeur sans faille, les proverbes béninois, et pour avoir été mon coach de footing. Je remercie Pierre-Antoine, Nicolas et Thibault de s'être joints à nous pendant ces précieux moments de détente. Je remercie Sarah pour avoir organisé et accueilli avec brio le fameux Club Ciné des Thésards de P6, et pour avoir partagé avec moi des bons et des moins bons moments de la thèse depuis le tout début. Je remercie Marion, Hour, Yanis et Fabien pour les discussions à la machine à café de Limours, les bons moments passés à IRS à Cracovie, et les débats enflammés sur Paris, Bordeaux, Stalingrad, et j'en passe. À Bordeaux justement, j'ai eu la chance d'être là aussi très bien entourée. J'ai une pensée pour les anciens, JB, Marc, Alan, Marie (et leur nouvelle petite famille), Bruno, Delphine (et leurs poèmes), et Zoé, à qui je n'oublie pas que je dois toujours un déjeuner post-Apollo. La relève a été brillamment assurée par les nouveaux, et je remercie Nikola pour ses réceptions mondaines les plus prisées des Capucins, Thomas pour de nombreuses discussions cinéphiles (même si techniquement, nous sommes les cinéphiles, pas la discussion), mon frère de thèse Jonathan et son frère Nicolas, pour les bières bues à Bordeaux, Paris, Lille, et les pieds bleus. Je remercie Sami qui continue à me supporter depuis le stage par lequel tout ceci à commencé, et Elsa, la très vénérable, de m'avoir si bien préparé à ma future ville d'adoption (du moins provisoire). Je n'oublie pas Vasileios, Philippe, et tous les autres qui ont partagé ces bons moments à l'IMB avec moi. Enfin, je remercie ma chère Léo, qui m'a accueillie tous les mois et sans qui cette dernière année de thèse aurait été bien différente.

Je remercie mes amis parisiens pour toutes ces années à arpenter les rues de Paris (ou du moins de Belleville) ensemble, et en particulier Jeanne et Tristan pour la frontale paillette et (bientôt !) le mariage de l'ambiance, Martin pour les conseils beauté capillaire, Sylvain et Chacha pour un voyage post-thèse qui s'annonce épique.

Résumé général

Motivation Dans des domaines comme l'inférence statistique ou la théorie de l'information, il est parfois intéressant d'adopter un point de vue géométrique pour résoudre certains problèmes : c'est la stratégie de la géométrie de l'information. L'idée fondamentale est que chaque élément d'une famille paramétrique de densités de probabilités {p θ , θ ∈ Θ} peut être représenté par son paramètre θ dans la variété des paramètres Θ ⊂ R d , et que ces lois peuvent être comparées en munissant cette variété d'une structure métrique, celle induite par la métrique de Fisher. Dans le cas des gaussiennes univariées N (m, σ 2 ), la géométrie de Fisher n'est autre que la géométrie hyperbolique : l'espace des paramètres (m, σ) est en bijection avec le demi-plan de Poincaré. De façon similaire, il est possible de relier la géométrie de l'information des processus gaussiens centrés multivariés stationnaires à la géométrie hyperbolique. Ces processus sont représentés par leurs matrices de covariance, des matrices hermitiennes définies positives qui possèdent en outre une structure de Toeplitz (à cause de leur stationnarité). Celles-ci peuvent être paramétrées par un ensemble de coefficients dans le disque unité complexe, qui sont ceux d'un processus autorégressif à maximum d'entropie [START_REF] Parker | Maximum entropy spectral analysis[END_REF]. Autrement dit, ces coefficients sont estimés à partir d'une observation du processus de façon à rajouter le moins d'information possible par rapport à cette observation. Barbaresco a introduit une métrique entropique qui définit une géométrie à courbure négative dans l'espace de ces coefficients, celle du polydisque de Poincaré [5]. Si un signal gaussien stationnaire peut être représenté par un point dans le polydisque -représentation qui a déjà fait ses preuves en traitement du signal radar [START_REF] Arnaudon | Riemannian medians and means with applications to radar signal processing[END_REF] -alors nous pouvons représenter un signal gaussien localement stationnaire par une courbe. Comme pour bon nombre d'applications, il est intéressant de s'affranchir de la paramétrisation de ces courbes afin d'analyser le processus indépendamment de la vitesse de réalisation dans le temps. Ceci constitue notre motivation pour étudier des formes de courbes dans une variété, bien que le travail de cette thèse soit plus général. v vi Structure riemannienne sur les courbes Il y a naturellement de nombreuses façons de comparer des courbes et leurs formes dans une variété. Une possibilité est de voir l'espace de ces courbes M lui-même comme une variété (de dimension infinie), et de le munir d'une structure riemannienne. Deux courbes peuvent ainsi être comparées en utilisant la distance géodésique, mais surtout cela permet de définir la notion de déformation optimale d'une courbe en une autre, et de faire des statistiques sur un ensemble de courbes grâce à des méthodes géodésiques. Pour modéliser l'espace des « formes », une stratégie usuelle consiste à considérer l'espace quotient S induit par l'action du groupe des reparamétrisations sur l'espace des courbes. Dans S, deux courbes sont considérées comme identiques si elles ne diffèrent que par leur paramétrisation, c'est-à-dire si elles traversent les mêmes points de l'espace mais à des vitesses différentes. L'espace des courbes M, des formes S et le groupe des reparamétrisations forment un fibré principal, et si la métrique définie sur M est invariante par rapport à l'action du groupe des reparamétrisations, alors elle induit une structure riemannienne sur l'espace des formes. C'est donc cette propriété qui nous a guidés dans notre choix de métrique sur M.

Dans toute structure riemannienne, la distance entre deux points se mesure comme la longueur du plus court chemin (géodésique) menant de l'un à l'autre. Aussi surprenant que cela puisse paraître, nous savons depuis 2005 grâce à Michor et Mumford [START_REF] Peter | Vanishing geodesic distance on spaces of submanifolds and diffeomorphisms[END_REF] que la métrique L 2 induit une distance nulle sur l'espace des formes : on peut toujours relier deux formes par un chemin arbitrairement court pour cette métrique. Cela a motivé l'étude de métriques plus élaborées qui ne mesurent plus seulement les différences de positions mais aussi de vitesses, voire de dérivées supérieures. Dans cette thèse, nous proposons une métrique entre les vitesses qui généralise aux courbes dans une variété non-linéaire (comme l'espace hyperbolique) la structure riemannienne du "square root velocity (SRV) framework" introduite par Srivastava et al. pour les courbes planes [START_REF] Srivastava | Shape analysis of elastic curves in euclidean spaces[END_REF]. Celle-ci est très utilisée dans les applications car elle devient la métrique L 2 dans un système de coordonnées approprié. Nous introduisons une métrique élastique, qui prend une expression compacte dans les coordonnées SRV, où chaque courbe est représentée par son point de départ et son champ de vitesse correctement renormalisé. Contrairement à la généralisation proposée par Zhang et al. dans [START_REF] Zhang | Video-based action recognition using rate-invariant analysis of covariance trajectories[END_REF], où les calculs sont transportés dans un même espace tangent grâce au transport parallèle, notre structure riemannienne lit directement les informations dans la variété M et est donc plus directement dépendante de sa géométrie. Par rapport au cas linéaire, nous ajoutons un terme de position qui induit une structure de fibré vectoriel au-dessus de la variété M vue comme l'ensemble des points de vii départ des courbes. Si l'équation d'une géodésique n'est pas explicite dans le cas non-linéaire, l'équation différentielle vérifiée par les géodésiques prend néanmoins une forme facilement exploitable dans les coordonnées SRV et présente l'avantage d'être facile à résoudre numériquement, nous permettant de construire la carte exponentielle. Celle-ci donne la géodésique partant d'une courbe dans une direction donnée, autrement dit sa déformation optimale dans cette direction. Nous obtenons la géodésique déformant une courbe en une autre grâce à une méthode de tir géodésique, où la « direction de tir » par carte exponentielle est corrigée itérativement en réduisant l'écart à la courbe cible grâce à un champ de Jacobi.

Afin de comparer les formes de deux courbes, nous construisons les géodésiques horizontales du fibré principal π : M → S, qui sont les plus courts chemins entre les fibres. Pour ce faire, nous caractérisons la décomposition de l'espace tangent à M en sous-espaces horizontal et vertical. Les vecteurs verticaux sont ceux qui reparamétrisent la courbe sans changer sa forme, et les vecteurs horizontaux sont leurs orthogonaux par rapport à la métrique G. Comme précisé plus haut, l'invariance par reparamétrisation de G permet d'induire une structure riemannienne sur l'espace quotient S, telle que la projection π : M → S soit une submersion riemannienne. Les géodésiques pour cette métrique sont les projections des géodésiques horizontales de l'espace total M. Ainsi, pour trouver la géodésique entre les formes de deux courbes c 1 et c 2 , nous fixons la paramétrisation de la première et nous cherchons la paramétrisation de c 2 la plus proche, qui s'obtient comme l'extrémité de la géodésique horizontale reliant c 1 à la fibre de c 2 . Cette reparamétrisation fournit un « matching optimal » entre les deux courbes. Pour l'obtenir, nous décomposons chaque chemin de courbes en une partie horizontale composée avec un chemin de reparamétrisations, et proposons un algorithme qui, pour deux courbes c 1 et c 2 données, réduit itérativement la distance de c 1 à la fibre de c 2 en alternant entre tir géodésique et calcul de la partie horizontale de la géodésique obtenue. Cet algorithme s'appuie sur une décomposition canonique d'un chemin dans un fibré et permettrait plus généralement de construire les géodésiques de l'espace de base d'un fibré principal quelconque à partir de celles de l'espace total.

Discrétisation et simulations Dans un objectif d'implémentation, nous introduisons ensuite une discrétisation de la structure riemannienne sur M, qui est elle-même une structure riemannienne sur la variété de dimension finie M n des « courbes discrètes » données par n points. Ce modèle discret est élaboré pour des courbes à valeurs dans une variété de courbure sectionnelle viii K constante, et pourrait se généraliser aux courbes dans une variété symétrique. Nous montrons sa convergence vers le modèle continu quand n → ∞, plus précisément la convergence de l'énergie discrète d'un chemin de courbes discrètes vers l'énergie continue du chemin de courbes limite. Comme dans le cas continu, nous caractérisons les géodésiques pour cette structure métrique, décrivons la carte exponentielle et le tir géodésique. Nous proposons une version discrète de la structure quotient de l'espace des formes, et ainsi de l'algorithme de matching optimal. Ce dernier permet, étant données deux séries de points et leurs courbes sous-jacentes (par exemple, des interpolations optimales), de redistribuer les points sur la deuxième courbe de façon à minimiser la distance discrète à la série de points de la première courbe. Toute les équations sont données dans trois cas : courbure positive K = 1, courbure nulle K = 0 et courbure négative K = -1, c'est-à-dire, en dimension 2, pour des courbes dans la sphère, le plan et le demi-plan hyperbolique. Des simulations réalisées avec le logiciel Matlab permettent de remarquer que les géodésiques pour notre métrique ont tendance à « rétrécir » les courbes au milieu de leur déformation, quel que soit l'espace considéré. Les simulations réalisées avec l'algorithme de matching optimal montrent que les points de la courbe cible sont redistribués de façon « naturelle » par rapport à ceux de l'autre courbe, et que la géodésique horizontale obtenue est de même longueur quelles que soient les paramétrisations choisies. Cette longueur commune fournit la distance entre les formes des courbes, et les géodésiques horizontales obtenues pour divers jeux de paramétrisation se superposent pour former la géodésique entre ces formes.

Applications radar Enfin, nous donnons un exemple d'application de l'étude de formes dans une variété au traitement statistique du signal radar. Nous utilisons la structure riemannienne exposée précédemment pour comparer des courbes dans le polydisque de Poincaré représentant des signaux radar localement stationnaires. L'intérêt de cette approche est d'exploiter la structure riemannienne pour faire des statistiques entre des signaux radar, dans un but de détection ou de reconnaissance de cibles. En guise d'exemple, nous construisons grâce à un flot de Karcher la moyenne de Fréchet de plusieurs signatures d'hélicoptère simulées, correspondant à l'observation d'un même hélicoptère avec de petites variations de la vitesse de rotation du rotor, comme on peut l'observer dans des cas réels, par exemple lors d'une manoeuvre qui incline l'hélicoptère. Ceci permet de construire une signature de référence de cet objet, et de l'utiliser pour faire de la reconnaissance de cibles. ix x Introduction

Related work

Many applications require the measure of dissimilarity between two shapes. Whether one is concerned with the outline of an object, represented by a closed curve, or by an open curve drawing the evolution of a given time process in a certain space, say a manifold, it is often interesting to consider unparameterized versions of these curves. In the former case, this allows us to study the shape of the object, and in the latter, to deal with the time process regardless of speed or pace. There are different ways to define the shape of a curve and a great deal more possibilities to measure dissimilarity between two shapes. We will not attempt here to make an exhaustive list, only consider a few examples.

An intuitive way to compare two unparameterized curves is given by the Hausdorff distance, which controls the maximum value of the distance between a point on one curve and its closest neighbor on the other. If X and Y are two curves in a normed vector space (V, • ), it is given by In some cases however, as the one shown in Figure 1.1, the Hausdorff distance does not properly translate the notion of similarity as perceived by the human eye. This is because it does not take into account the continuity of the curves, which are simply seen as sets of points. This is remedied by the Fréchet distance [START_REF] Fréchet | Sur quelques points du calcul fonctionnel[END_REF], which considers all orientation-preserving matchings between two parameterized curves The infimum is taken over all continuous increasing functions ϕ, ψ that fix the extremities of the curves. A well-known interpretation of this involves a person walking his dog, one curve representing the trajectory of the dogowner and the other, the trajectory of the dog. Both can control their speed but never go backwards. The Fréchet distance is the size of the smallest leash necessary for a leash-minimizing walk. As a maximum distance, the Fréchet distance measures the discrepancy between two curves in two points only, instead of measuring dissimilarity along the whole curves. An integral version was defined in [START_REF] Buchin | On the computability of the fréchet distance between triangulated surfaces[END_REF], however it does not satisfy the triangular inequality.

Beyond computing distances between shapes, a desirable goal in many applications is to do statistics on a set of shapes, e.g. compute the mean, perform classification or even principal component analysis. For this purpose, considering shapes as elements of a shape manifold that we equip with a Riemannian structure provides a convenient framework. In this infinitedimensional shape manifold, points represent shapes and the distance between two shapes is given by the length of the shortest path linking them -the geodesic. This approach allows us to do more than simply compute distances : it enables us to define the notion of optimal deformation between two shapes -how does one shape optimally deform into another. Schematically, if the surface on the left-hand side of Figure 1.2 represents the manifold of shapes and we equip that manifold with a Riemannian metric, then the geodesic linking two shapes c and c 1 corresponds to the optimal (in the sense given by the metric) deformation from c to c 1 , as shown on the right-hand side of the figure. Another advantage of this approach is that it allows us to locally linearize the shape space at a given point, e.g. the barycenter of a set of shapes, by considering the tangent space at that point, and do statistics in that flat space using geodesic methods. Now, how can we rigorously define this shape manifold ? The idea of a shape space as a Riemannian manifold was first developed by Kendall [START_REF] David | Shape manifolds, procrustean metrics, and complex projective spaces[END_REF], who defines shapes as "what is left" of a curve after the effects of translation, rotation and changes of scale are filtered out. Indeed, the key idea that we will use here to mathematically define a shape is to quotient out the action of certain transformations, the choice of which differs according to the applications. The shapes considered by Kendall are represented by labelled points in the Euclidean space and the shape spaces are finite-dimensional.

More recent works on the other hand, deal with continuous curves with values in a Euclidean space or a nonlinear manifold, and infinite-dimensional shape spaces. There exist two distinct, complementary approaches to define the shape space [START_REF] Bauer | Overview of the geometries of shape spaces and diffeomorphism groups[END_REF]. One possibility is to consider the transitive left action of the group of spatial deformations G s -namely the set of diffeomorphisms of the ambient space -on the set of shapes (or unparameterized curves) S G s × S → S, (ϕ, c) → ϕ • c.

Equipping G s with a G s -invariant metric d Gs , we can measure the distance between two shapes c 1 and c 2 as the minimal cost of a deformation that transforms one shape into the other

d S (c 1 , c 2 ) = inf{d Gs (Id, ϕ) | ϕ ∈ G s , ϕ • c 1 = c 2 }.
This is the setting of the so-called outer metrics, because a deformation of a curve is realized by the deformation of the entire ambient space. The second approach consists in quotienting out the set G t of temporal deformations, which have an action of reparameterizing the curve. For this we consider the transitive right action of G t on the space of parameterized curves M,

M × G t → M, c → c • ϕ.
Given a G t -invariant Riemannian metric on M, we can measure the distance between the shapes of two curves as the distance between a fixed parameterization of one and an optimal reparametrization of the other where we use the same notation for a curve and its shape. We place ourselves in the setting of these so-called inner metrics, where the deformation of the curve does not affect the ambient space.

d S (c 1 , c 2 ) = inf{d M (c 1 , c 2 • ϕ) | ϕ ∈ G t }, (1.1)
To be more precise, we consider the infinite-dimensional manifold M of smooth parameterized curves

c : [0, 1] → M, t → c(t)
with values in a certain ambient space M , e.g. a Euclidean space or a manifold. Temporal deformations, or reparameterizations, of these curves are increasing diffeomorphisms of the time interval on which the curves are defined, i.e. elements of G t = Diff + ([0, 1]). We define a Riemannian metric on M, i.e. we equip each tangent space T c M of M at point c with a scalar product, denoted by

G c : T c M × T c M → R, (w, z) → G c (w, z).
The tangent vectors w, z are infinitesimal deformations of the curve c, and can be seen as infinitesimal vector fields along c, as shown schematically in the left-hand side of Figure 1.3. We require that the scalar product G be reparameterization invariant, i.e. that

G c•ϕ (w • ϕ, z • ϕ) = G c (w, z), ∀ϕ ∈ G t ,
so that the distance between two curves c 1 and c 2 does not change when we reparameterize them the same way

d M (c 1 • ϕ, c 2 • ϕ) = d M (c 1 , c 2 ), ∀ϕ ∈ G t .
This however does not guarantee that the obtained distance is completely independant of the parameterizations of the curves : if we reparameterize c 1 and c 2 in two distinct ways, then the distance does change. To illustrate this phenomenon, let us anticipate and show a simulation obtained using the reparameterization invariant metric defined in Chapter 2 and the algorithms presented in Chapter 3. In Figure 1.4, we show the optimal deformations (or geodesics) between two curves for three different pairs of parameterizations, represented by the way the points are distributed on the curves. These points indicate the position of the curve at regular times, so that on the left-hand side, the speed is constant on both curves, in the middle the speed decreases for both curves, and on the right-hand side the speed of the blue curve decreases while the speed of the black curve increases. We can see that the geodesic does not significantly change when the curves are reparameterized the same way, whereas it does when they are reparameterized differently. This is reflected in the lengths of these geodesics, and therefore in the distances between the curves, yielding

d M (c 1 , c 2 ) = d M (c 1 • ϕ, c 2 • ϕ) = d M (c 1 • ϕ, c 2 • ψ),
for ϕ = ψ. This explains the need to take the infimum over all reparameterizations in Equation (1.1), in order to compute the distance between the shapes of parameterized curves.

The question now is : how do we choose an appropriate reparameterization invariant metric G on M ? The literature for plane curves is quite abundant. As early as 1998, Younes considered in [START_REF] Younes | Computable elastic distances between shapes[END_REF] shapes of plane curves as elements of an infinite dimensional shape space and adopted a Riemannian point of view, leading to the computation of geodesic paths and optimal matchings between plane curves. Since then, convenient metrics to compare shapes have been actively studied both from a theoretical and an applied point of view. An important step was the discovery that the reparameterization invariant L 2 -metric,

G 0 c (w, z) = 1 0 w(t), z(t) |c (t)|dt, w, z ∈ T c M,
which can be written more compactly using the notation d = |c |dt to refer

CHAPTER 1. INTRODUCTION

to integration according to arc length,

G 0 c (w, z) = w, z d , w, z ∈ T c M,
induces a vanishing geodesic distance d S on the shape space [START_REF] Peter | Vanishing geodesic distance on spaces of submanifolds and diffeomorphisms[END_REF] : Michor and Mumford found that two shapes can be connected by an arbitrarily short path in that space. This was a motivation to explore more elaborate metrics, such as almost local metrics, Sobolev metrics, and elastic metrics, to name a few. Sobolev metrics are stronger versions of G 0 defined using linear combinations of higher order derivatives of the tangent vectors. Although definitions vary, a general form of a n-order Sobolev metric on the space of

R d -valued curves can be G n c (w, z) = n i=0 a i D i w, D i z d , w, z ∈ T c M,
where D denotes derivation with respect to arc length, i.e. D w = w /|c |, and the a i 's are constants. These metrics have been carefully studied in [START_REF] Acg Mennucci | Properties of sobolev-type metrics in the space of curves[END_REF] and [START_REF] Peter | An overview of the riemannian metrics on spaces of curves using the hamiltonian approach[END_REF], and continue to be [START_REF] Bruveris | Geodesic completeness for sobolev metrics on the space of immersed plane curves[END_REF], [START_REF] Nardi | Geodesics on shape spaces with bounded variation and sobolev metrics[END_REF], [START_REF] Bauer | Why use sobolev metrics on the space of curves[END_REF]. Many theoretical questions have been examined: completion of the space of curves for these metrics in the case n = 1, 2 were given in [START_REF] Acg Mennucci | Properties of sobolev-type metrics in the space of curves[END_REF] by Mennucci et al. -who also showed that the Fréchet distance is induced by the "Finsler L ∞ -metric", defined as the L ∞ -norm of the normal projection of the tangent vector -and the geodesic equations in the space of curves and the shape space were worked out in [START_REF] Peter | An overview of the riemannian metrics on spaces of curves using the hamiltonian approach[END_REF] by Michor and Mumford. The geodesic equation for first-order Sobolev metrics was found to be locally but not globally well-posed [START_REF] Peter | An overview of the riemannian metrics on spaces of curves using the hamiltonian approach[END_REF], while global existence of geodesics for Sobolev metrics G n with n ≥ 2 was proven in [START_REF] Bruveris | Geodesic completeness for sobolev metrics on the space of immersed plane curves[END_REF] by Bruveris et al. The first-order limit-case,

G 1,∞ c (w, z) = 1 L(c) D w, D z d , w, z ∈ T c M,
of which the metric of [START_REF] Younes | Computable elastic distances between shapes[END_REF] can be seen as a precursor, was examined in [START_REF] Younes | A metric on shape space with explicit geodesics[END_REF] by Younes et al. This elastic metric is mapped to an L 2 -metric using new shape coordinates, where a curve is represented by the complex square root of its speed. This gives explicit expressions for the geodesics and a closed form for the geodesic distance. In [START_REF] Mio | On shape of plane elastic curves[END_REF], Mio et al. consider a similar metric but put different weights on the tangential D w T , D z T and normal parts D w N , D z N of the arc-length derivatives of the tangent vectors, thereby defining a 2-parameter family of elastic metrics

G a,b c (w, z) = a 2 D w N , D z N + b 2 D w T , D z T d , w, z ∈ T c M.

RELATED WORK

As shown schematically on the right-hand side of Figure 1.3, the tangential part w T of an infinitesimal deformation w of a curve c, is simply its pointwise projection on the speed vector field of c, and the normal part is the remaining component w N = ww T . Parameters a and b respectively account for the degree of "stretching" and "bending" of the curve.

A particularly interesting choice of parameters is a = 1 and b = 1/2, since in that case the metric can also be mapped to an L 2 -metric. Indeed, Srivastava et al. introduced in [START_REF] Srivastava | Shape analysis of elastic curves in euclidean spaces[END_REF] another shape coordinate system, where a curve, provided its speed never vanishes, is represented by its speed renormalized by the square root of its (Euclidean) norm

c → q := c |c | .
This representation is referred to as the square root velocity representation. In these coordinates, the elastic distance between two curves is simply the L 2 -distance. In other words, the distance between two curves c 0 and c 2 is the pointwise distance between their renormalized speeds q 0 and q 1 ,

d 2 G 1, 1 2 (c 0 , c 2 ) = d 2 L 2 (q 0 , q 1 ) = 1 0 |q 1 (t) -q 0 (t)| 2 dt.
This property makes the G 1, 1 2 metric and the square root velocity framework particularly interesting for applications, as it is usually not easy to find explicit equations for the geodesics in shape analysis. An extension of this framework to any elastic metric with coefficients 4b 2 ≥ a 2 was proposed in [START_REF] Bauer | Constructing reparameterization invariant metrics on spaces of plane curves[END_REF], and several works focused on the existence of optimal reparameterizations in the SRV framework and how to compute them [START_REF] You | A riemannian approach for computing geodesies in elastic shape analysis[END_REF], [START_REF] Lahiri | Precise matching of pl curves in mathbbr n in the square root velocity framework[END_REF], [START_REF] Bruveris | Optimal reparametrizations in the square root velocity framework[END_REF].

The literature for similar metrics on spaces of manifold-valued curves is less profuse, although the topic in general is meeting growing interest [START_REF] Singh | A hierarchical geodesic model for diffeomorphic longitudinal shape analysis[END_REF], [START_REF] Schiratti | Mixed-effects model for the spatiotemporal analysis of longitudinal manifold-valued data[END_REF]. Focusing on inner metrics on the space of curves, Sobolev metrics for curves in a manifold have been studied in full generality by Bauer et al. in [START_REF] Bauer | Sobolev metrics on shape space of surfaces[END_REF]. It is shown that Sobolev metrics (of order one or higher) on manifoldvalued curves overcome the degeneracy of the L 2 -metric, and the geodesic equation is formulated in terms of the gradients of the metric with respect to itself. The square root velocity framework was extended to manifold-valued curves by Zhang et al. in [START_REF] Zhang | Video-based action recognition using rate-invariant analysis of covariance trajectories[END_REF] and applied to curves in the space of symmetric definite positive matrices for speech and movement recognition purposes. This same metric was further studied for the special case of spherical trajectories in [START_REF] Zhang | Phase-amplitude separation and modeling of spherical trajectories[END_REF]. Extensions to curves in a Lie group was proposed by Celledoni et al. in [START_REF] Celledoni | Shape analysis on lie groups with applications in computer animation[END_REF], and recently to homogeneous spaces in [START_REF] Celledoni | Shape analysis on homogeneous spaces[END_REF].

CHAPTER 1. INTRODUCTION

Here we are also concerned with the study of curves with values in nonlinear manifolds, in particular, but not restricted to, manifolds of constant sectional curvature such as the hyperbolic plane and the sphere. We use the square root velocity coordinates to define our metric structure on the space of manifold-valued curves, however this leads to a different metric than the one presented in [START_REF] Zhang | Video-based action recognition using rate-invariant analysis of covariance trajectories[END_REF]. We generalize the flat case, in the sense that in the particular case of plane curves, our metric gives the following distance between two curves

d 2 (c 0 , c 2 ) = |c 2 (0) -c 0 (0)| 2 + 1 0 |q 1 (t) -q 0 (t)| 2 dt.
Our motivation is also, to our knowledge, novel, as we use information geometry to represent locally stationary Gaussian radar signals as curves in the hyperbolic plane, and apply shape analysis to compare and average these curves and thereby, the radar signals they represent. A similar strategy is adopted in [START_REF] Ruiz | Radar detection for nonstationary doppler signal in one burst based on information geometry distance between paths[END_REF] using the Fréchet distance.

Motivation

If the shape analysis of plane or spherical curves can be easily motivated by applications involving shapes of 2D objects or trajectories on the Earth, the practical interest of considering curves in the hyperbolic plane is less obvious. And yet, as previously mentioned, the motivation of this work comes from the need to do statistics on curves in the 2-dimensional hyperbolic space, in the context of radar signal processing.

The hyperbolic geometry of Gaussians

Let us first say a few words about the link between hyperbolic geometry and Gaussian densities made by information geometry. Information geometry provides a geometrical approach to various fields such as statistical inference, information theory or signal processing [START_REF] Fréchet | Sur l'extension de certaines évaluations statistiques au cas de petits échantillons[END_REF], [START_REF] Burbea | Entropy differential metric, distance and divergence measures in probability spaces: A unified approach[END_REF], [START_REF] Amari | Differential geometry in statistical inference[END_REF]. The key idea is that elements of a parametric family of probability densities {p θ , θ ∈ Θ} can be seen as points in the manifold of parameters Θ, and compared through the definition of a Riemannian structure on that space. Intuitively, it is easy to see that the Euclidean metric is usually not a good choice. For example, if we identify univariate Gaussian distributions with their mean and standard deviation (m, σ) in the upper half-plane R × R * + , then we easily conceive that two densities N (m 1 , σ 1 ) and N (m 2 , σ 1 ) with different means but the same standard deviation get "closer" to each other as their common standard 1.2. MOTIVATION deviation increases, meaning that the distance between the points of coordinates (m 1 , σ 1 ) and (m 2 , σ 1 ) in the upper half-plane should be greater than the distance between the points (m 1 , σ 2 ) and (m 2 , σ 2 ) for σ 2 > σ 1 . A more pertinent Riemannian structure is given by the Fisher information metric. If the parameter θ ∈ R d is d-dimensional and E θ denotes the expected value with respect to the density p θ , the matrix form of the metric is given by the Fisher information matrix,

g ij (θ) = I(θ) ij = E θ ∂ ∂θ i ln p θ (X) ∂ ∂θ j ln p θ (X) , (1.2)
for any 1 ≤ i, j ≤ d and θ = (θ 1 , . . . , θ d ). This metric is chosen, among other reasons, because it has statistical meaning : in parameter estimation, the Fisher information measures the "amount of information" on the parameter θ contained in data sampled from the density p θ ; it also gives a fundamental limit to the precision at which one can estimate this θ, in the form of the Cramer-Rao bound. In the case of univariate Gaussian densities N (m, σ 2 ), Fisher geometry amounts to hyperbolic geometry. More precisely, the space of parameters (m, σ) equipped with the Fisher Information metric is in bijection with the hyperbolic half-plane via the change of variables (m, σ) → ( m √ 2 , σ). Indeed, with this rescaling of the mean, the Fisher Information matrix becomes

g(m, σ) = 2 1 σ 2 0 0 1 σ 2
, which, up to the factor 2, defines the Riemannian metric of the well-known hyperbolic half-plane. Note that this is coherent with the example given above, since in the hyperbolic half-plane the distance between the points of coordinates (m 1 , σ) and (m 2 , σ) decreases as σ increases for fixed values of m 1 , m 2 . As we will see in Chapter 4, hyperbolic geometry is closely linked to the Fisher geometry of multivariate Gaussians as well. The differential geometry of Gaussians has proved useful for applications, e.g. in image processing where in the image model, each pixel is represented by a univariate Gaussian distribution [START_REF] Angulo | Morphological processing of univariate gaussian distribution-valued images based on poincaré upperhalf plane representation[END_REF], and in radar signal processing [START_REF] Barbaresco | Interactions between symmetric cone and information geometries: Bruhat-tits and siegel spaces models for high resolution autoregressive doppler imagery[END_REF], [START_REF] Arnaudon | Riemannian medians and means with applications to radar signal processing[END_REF], [START_REF] Pilté | Tracking quality monitoring based on information geometry and geodesic shooting[END_REF] where the echo corresponding to an element of space is represented by a stationary Gaussian process.

Radar signal processing

The motivation behind the work of this thesis stems from radar signal processing. The basic functioning of a radar is well-known : a radar sends electromagnetic waves in a certain direction and deduces the presence of objects in CHAPTER 1. INTRODUCTION its environment through the echoes it receives. The distance to the object is deduced from the time interval between the emission of the radio wave and the reception of the echo, and the radial velocity of the object can be measured through the Doppler effect. Supposing the radar is looking in a fixed direction, it is useful to "discretize" the space in that direction by dividing it into several distance cells -as many as the range of the radio wave allows. In the laps of time between two emissions, we obtain an echo for each distance cell for which we can measure an information of phase r and of amplitude θ, obtaining a complex number z = re iθ . If the radar sends a burst of n pulses in a fixed direction, then we receive a vector of n complex data z = (z 1 , . . . , z n ) for each distance cell, that we will call a signature.

Target detection relies on the intuitive principle that an element of space contains a target if it is significantly different, in a way to be defined, from its environment. To know whether a certain distance cell contains a target or not, we compare its signature to those of its neighboring cells using a statistical test -eliminating the cells immediately adjacent which could also be affected by the target. This comparison allows us to remove the noise and clutter (echoes returned from the ground, the sea, or atmospheric turbulences) from consideration and isolate the echo due to the target. Besides the ability to detect and locate, the capacity to distinguish between different targets is also a crucial issue. This means being able to compare the signature of a cell suspected to contain a target to reference signatures which play the role of templates. These templates should correspond to the "average" signature observed for a certain class of targets -planes, helicopters, drones, or for a certain model of helicopter, plane, drone.

Classical CFAR (Constant False Alarm Rate) detection methods are based on a likelihood-ratio test taking into account the information of the environing cells and a threshold chosen such that the probability to have a false alarm, i.e. to decide on a target when there is none, stays constant. They usually take as input the output of FFT or Doppler filter banks applied to the measured signal. While these methods are satisfactory in many cases, they present low resolution limitations in the presence of dense, inhomogeneous clutter and when the number n of available observations is small. To remedy these low-resolution issues, Burg suggested as early as in the late 1960's a maximum entropy approach based on autoregressive processes [START_REF] Parker | Maximum entropy spectral analysis[END_REF], which was further developed and specifically applied to radar signal processing by Barbaresco [5]. In this model, a signature z = (z 1 , . . . , z n ) T ∈ C n measured for a given distance cell after a burst of n pulses, is assumed to be the realization of a centered, stationary circularly-symmetric [START_REF] Robert | Circularly-symmetric gaussian random vectors[END_REF] Gaussian vector Z = (Z 1 , . . . , Z n ) T which is therefore entirely described by its covariance matrix Σ = E(ZZ * ), where E is the expected value and Z T , Z * are respectively the transpose and transconjugate of Z. Since Z is considered stationary, its autocorrelation values r k = E(Z i Z i+k ), k = 0, . . . , n -1, depend only on the lag k, and Σ is a Toeplitz matrix

Σ =       r 0 r 1 • • • r n-1 r 1 . . . . . . . . . . . . . . . . . . r 1 r n-1 • • • r 1 r 0      
, where z is the conjugate of z. Burg showed that the maximum entropy process with respect to these autocorrelation constraints -i.e. the one that adds the fewest assumptions on the data -is an autoregressive process. In fact, there is a one-to-one correspondence between the Toeplitz covariance matrix and the coefficients of this autoregressive process [START_REF] Verblunsky | On positive harmonic functions[END_REF], [START_REF] William F Trench | An algorithm for the inversion of finite toeplitz matrices[END_REF]. The idea of Barbaresco is to represent each signature z = (z 1 , . . . , z n ) ∈ C n by the Toeplitz covariance matrix of the underlying signal, and to equip the space T + n of Toeplitz, hermitian, positive, definite matrices of size n with a Riemannian metric from information geometry. As we will see in Chapter 4, this strategy coupled with Burg's autoregressive approach is particularly interesting because the autoregressive coordinate space R + × D n-1 (D is the complex unit disk) of Toeplitz matrices becomes the Poincaré polydisk R + × D n-1 (D is the usual Poincaré disk) when equipped with a certain information geometry metric. Representing a stationary Gaussian radar signal by a point of the Poincaré polydisk allows us to compare, average and do statistics on stationary radar signals using the inherent Riemannian structure of the polydisk.

In this thesis, we are interested in configurations where the radar signal is non stationary, in the presence of very inhomogeneous clutter or when the target moves during the time interval of the burst. We make the assumption that the underlying process of a radar signature is locally stationary, and following the representation system introduced by Barbaresco, we identify each stationary portion with a point in the Poincaré polydisk. That way, statistical detection or recognition tests on locally stationary radar signals can be performed through the statistical study of the corresponding curves in the Poincaré polydisk.

CHAPTER 1. INTRODUCTION

Main results

In this section, we give a condensed overview of the main contributions and results of this thesis. In Chapter 2, we focus on the comparison of shapes of curves that take their values in a Riemannian manifold M . The curves c : [0, 1] → M that we consider are open and oriented, and we assume that their velocity c never vanishes. Their set is the space of smooth immersions

M = Imm([0, 1], M ) = {c ∈ C ∞ ([0, 1], M ), c (t) = 0 ∀t ∈ [0, 1]},
which is an open submanifold of the Fréchet manifold C ∞ ([0, 1], M ). We equip it with a Riemannian metric, i.e. a scalar product G c on each tangent space at point c ∈ M composed of infinitesimal vector fields along c,

G c : T c M × T c M → R, (w, z) → G c (w, z), T c M = {w ∈ C ∞ ([0, 1], T M ) : w(t) ∈ T c(t) M ∀t ∈ [0, 1]}.
This Riemannian structure is introduced in Section 2.2. We require that G be invariant with respect to the right action of the set of reparameterizations, i.e. increasing diffeomorphisms

ϕ ∈ Diff + ([0, 1]) of the time interval [0, 1], M × Diff + ([0, 1]) → M, (c, ϕ) → c • ϕ, (1.3)
so that it induces a Riemannian metric on the quotient space M/Diff + ([0, 1]) of unparameterized curves or shapes. Since the L 2 -metric is pathological [START_REF] Peter | Vanishing geodesic distance on spaces of submanifolds and diffeomorphisms[END_REF], we define a metric that takes into account first-order derivatives of the tangent vectors. We choose to extend to manifold-valued curves the metric of the square root velocity (SRV) framework [START_REF] Srivastava | Shape analysis of elastic curves in euclidean spaces[END_REF], where different coefficients are given to the normal and tangential parts of the tangent vector's derivative. It is defined using arc-length covariant derivation ∇ w = ∇ c w/|c |, its projection ∇ w T = ∇ w, v v on the unit velocity field v = c /|c | and the normal complement

∇ w N = ∇ w -∇ w T , G c (w, z) = w(0), z(0) + ∇ w N , ∇ z N + 1 4 ∇ w T , ∇ z T d . (1.4)
Arc-length derivation and integration d = |c |dt guarantee that G is reparameterization invariant, just like all elastic metrics.

G a,b c (w, z) = w(0), z(0) + a 2 ∇ w N , ∇ z N + b 2 ∇ w T , ∇ z T d . (1.5) 1.3. MAIN RESULTS

13

The choice of limiting the H 0 part to t = 0 is motivated by the extension of the flat case [START_REF] Mio | On shape of plane elastic curves[END_REF] -in which the metric contains only an H 1 term -as well as the computability, and induces a fiber bundle structure over the base manifold M seen as the set of origins of the curves. In the context of plane curves, the case a = 2b = 1 is interesting because the metric thereby defined is the L 2 -metric in the SRV coordinates. While we do not have such a nice property for curves in a nonlinear manifold, we obtain a simple form for G in that space of coordinates, which will prove convenient to derive the geodesic equation.

Notation. We denote by s the parameter of paths in the space of curves M and by t the parameter of a curve in M . For any path of curves s → c(s, •) the corresponding derivatives are denoted by c s = ∂c/∂s and c t = ∂c/∂t , and we also use the notations ∇ s = ∇ ∂c/∂s and ∇ t = ∇ ∂c/∂t . We represent each curve by its origin and its velocity field renormalized by the square root of its norm via the bijection M → T M, c → x := c(0), q := c t / |c t | .

(1.6)

Proposition 1.1. The metric (1.4) can be written in terms of the square root velocity coordinates (1.6) as

G c (w, z) = w(0), z(0) + 1 0 ∇ w(t) q, ∇ z(t) q dt, (1.7)
for any curve c ∈ M, and vector fields w, z ∈ T c M along c.

Remark 1.1. The term ∇ w q denotes the covariant derivative of the vector field q in the direction given by the vector field w. More precisely, if s → c(s, •) is a path of curves verifying c(0, t) = c(t) and c s (0, t) = w(t) for all t ∈ [0, 1], and q(s, t) = c t (s, t)/ |c t (s, t)| denotes its SRV representation, then

∇ w(t) q = ∇ s | s=0 q(s, t), t ∈ [0, 1],
where ∇ s denotes the covariant derivative along the curve s → c(s, t).

In Section 2.2.4, we give the geodesic distance induced by G.

Definition 1.2. The geodesic distance associated to the Riemannian metric G between two points c 0 , c 1 ∈ M is given by the length of the lengthminimizing path

[0, 1] s → c(s) ∈ M such that c(0) = c 0 and c(1) = c 1 , d(c 0 , c 1 ) = inf c(0)=c 0 ,c(1)=c 1 L(c) = inf c(0)=c 0 ,c(1)=c 1 1 0 G c s (s), c s (s) ds. (1.8) 14 CHAPTER 1. INTRODUCTION
We compare the geodesic distance induced by G to that of a different generalization of the SRV framework [START_REF] Zhang | Video-based action recognition using rate-invariant analysis of covariance trajectories[END_REF], where each curve is represented by its origin and its square root velocity field q parallel transported to the tangent plane to its origin. While they both coincide when the base manifold M is flat, we show that they differ by a curvature term when M is nonlinear. Proposition 1.2 (Geodesic distance). The geodesic distance induced by G between two curves c 0 , c 1 ∈ M is given by

d(c 0 , c 1 ) = inf c(0)=c 0 ,c(1)=c 1 1 0 |x s (s)| 2 + 1 0
|∇ s q(s, t)| 2 dt ds.

(1.9)

It can also be written as a function of the transported square root velocity function q(s)

∈ C ∞ ([0, 1], T c 0 (0) M ) of c(s), d(c 0 , c 1 ) = inf c(0)=c 0 ,c(1)=c 1 1 0 |x s (s)| 2 + 1 0 |q s (s, t) + Ω(s, t)| 2 dt ds, (1.10)
where Ω is a curvature term measuring the holonomy along a rectangle of infinitesimal width.

The compact form (1.7) of G in the SRV coordinates allows us to derive the geodesic equations in Section 2.3, by searching for the critical points of the energy functional E :

C ∞ ([0, 1], M) → R + , E(c) = 1 0 G c s (s), c s (s) ds.
(1.11)

Definition 1.3. The geodesics are (locally) the length-minimizing paths, i.e. those achieving the infimum in (1.8).

The Cauchy-Schwarz inequality gives L(c) 2 ≤ 2 E(c), and since the minimizers of the length L have velocity of constant G-norm, this is an equality for geodesics, which are therefore also the minimizers of the energy (1.11).

Proposition 1.3 (Geodesic equation). Let [0, 1]
s → c(s) ∈ M be a path of curves. It is a geodesic of M if and only if its SRV representation s → (x(s), q(s)) verifies

∇ s x s (s) + r(s, 0) = 0, ∀s (1.12a)
∇ 2 s q(s, t) + |q(s, t)| r(s, t) + r(s, t) T = 0, ∀t, s (1.12b)

MAIN RESULTS

where the vector field r depends on the curvature tensor R of the base manifold M r(s, t) = 1 t R(q, ∇ s q)c s (s, τ ) τ,t dτ, and w(s, τ ) τ,t denotes the parallel transport of the vector field w(s, •) along c(s, •) from c(s, τ ) to c(s, t).

Remark 1.2. In the flat case M = R d , the curvature term r vanishes and we obtain ∇ s x s (s) = 0, ∇ 2 s q(s, t) = 0 for all s and t : we recover the fact that the geodesic between two curves (x 0 , q 0 ) and (x 1 , q 1 ) in the SRV coordinates is composed of a straight line s → x(s) between the origins and an L 2 -geodesic s → q(s, •) between the renormalized speeds.

These equations are easily numerically solved to construct the exponential map, an algorithm that computes the geodesic path s → c(s) starting from c 0 ∈ M at speed w ∈ T c 0 M -i.e. that optimally "deforms" the curve c 0 in the direction given by w.

Algorithm 1.1 (Exponential map). Set c(0, t) = c 0 (t) and c s (0, t) = w(t) for all t ∈ [0, 1]. For k = 0, . . . , m -1, set s = kε with ε = 1/m and 1. update the position using the exponential map of the base manifold c(s + ε, t) = exp M c(s,t) (ε c s (s, t)), ∀t ∈ [0, 1], 2. compute acceleration ∇ s c s (s, 0) = ∇ s x s (s) at time t = 0 using Equation (1.12a) and propagate to ∇ s c s (s, t), t > 0, using Equation (1.12b),

3. update the speed c s (s + ε, t) = [ c s (s, t) + ε∇ s c s (s, t) ] s,s+ε , ∀t ∈ [0, 1].

Updating the position using the exponential map of the base manifold is a canonical choice to construct segments. Based on the exponential map, the geodesic between two curves c 0 , c 1 ∈ M is computed through geodesic shooting, which allows us to iteratively find the appropriate shooting direction w ∈ T c 0 M that will deform c 0 to c 1 through the exponential map defined by Algorithm 1.1 : exp M c 0 (w) = c 1 . At each step, the difference between the geodesic obtained by "shooting" in the current shooting direction w and the desired geodesic is measured through a Jacobi field J of initial value J(0) = 0 using the L 2 inverse exponential map log L 2 (most canonical choice). We then use the initial velocity ∇ s J(0) to update the shooting direction. Definition 1.4. A Jacobi field s → J(s) along a geodesic s → c(s) is a vector field measuring the variation between c and an infinitesimally close geodesic: there exists a family a → c(a, •) of geodesics such that c(0, s) = c(s) and J(s) = ∂ a c(0, s) for all s. Note that this algorithm could fail to converge if the curvature of the base manifold M is too high. This algorithm requires the characterization of the Jacobi fields for G on M, which is given in Section 2.3.3. In Section 2.4, we study the quotient structure of the space of shapes or unparameterized curves for general elastic metrics (1.5) and metric (1.4) in particular. Restricting to the set M f of immersions on which the diffeomorphism group acts freely, the right action (1.3) of the diffeomorphisms group defines a principal bundle π : M f → S f = M f /Diff + ([0, 1]), the fibers of which are the sets of all the curves that are identical modulo reparameterization. Definition 1.6. If G is a topological group, a principal G-bundle is a fiber bundle π : M → B with a continuous right action G × M → M such that G preserves the fibers and acts freely and transitively on them.

Any tangent vector w to M in c can then be decomposed as the sum of a vertical part w ver tangent to the fiber of c, which has an action of reparameterizing the curve without changing its shape, and a G-orthogonal horizontal part w hor . Whatever the metric G, the subspace of vertical vectors is given by

Ver c = ker T c π = mv = mc /|c | : m ∈ C ∞ ([0, 1], R), m(0) = m(1) = 0 .
The horizontal subspace on the other hand depends on the metric G.

MAIN RESULTS

Proposition 1.4 (Horizontal part of a vector). Let c ∈ M. Then h ∈ T c M is horizontal for the elastic metric G a,b if and only if

(a/b) 2 -1 ∇ t h, ∇ t v -∇ 2 t h, v + |c | -1 ∇ t c , v ∇ t h, v = 0.
In particular, for a = 2b = 1, horizontal vectors verify

3 ∇ t h, ∇ t v -∇ 2 t h, v + |c | -1 ∇ t c , v ∇ t h, v = 0.
The vertical and horizontal parts of a tangent vector w ∈ T c M are given by w ver = mv, w hor = wmv, where the real function m ∈ C ∞ ([0, 1], R) verifies m(0) = m(1) = 0 and the ordinary differential equation

m -∇ t c /|c |, v m -4|∇ t v| 2 m = ∇ 2 t w, v -3 ∇ t w, ∇ t v -∇ t c /|c |, v ∇ t w, v .
(1.13)

Characterizing the horizontal subspace associated to the quotient structure is key to compute geodesics between shapes. Definition 1.7. A horizontal path s → c(s) ∈ M in a principal bundle π : M → B is a path with horizontal speed c s (s) ∈ Hor c(s) at time s = 0, or equivalently at all time s ([36], §26.12).

Since G is reparameterization invariant, there exists a Riemannian metric G on the shape space S f such that π is a Riemannian submersion from (M f , G) to (S f , G). The geodesic s → c(s) for G between the shapes of two curves c 0 and c 1 is the projection c = π(c h ) of the horizontal geodesic s → c h (s) linking c 0 to an element c 1 • ϕ of the fiber of c 1 in M. To compute the optimal matching given by c h and ϕ, we decompose any path of curves s → c(s) in M into a horizontal path c hor composed with a path of reparameterizations ϕ,

c(s, t) = c hor (s, ϕ(s, t)) ∀s, t ∈ [0, 1].
(1.14) Definition 1.8. We call c hor the horizontal part of the path c with respect to the metric G.

Proposition 1.5. The horizontal part of a path of curves c is at most the same length as c L G (c hor ) ≤ L G (c).

The existence and uniqueness of the horizontal part of a path is given by the following result.

CHAPTER 1. INTRODUCTION

Proposition 1.6 (Horizontal part of a path). Let s → c(s) be a path in M. Then its horizontal part is given by c hor (s, t) = c(s, ϕ(s) -1 (t)), where the path of diffeomorphisms s → ϕ(s) is solution of the partial differential equation

ϕ s (s, t) = m(s, t)/|c t (s, t)| • ϕ t (s, t), (1.15)
with initial condition ϕ(0, •) = Id, and where m(s)

: [0, 1] → R, t → m(s, t) verifies m(s, 0) = m(s, 1) = 0 and is solution for all s of the ODE m tt -∇ t c t /|c t |, v m t -4|∇ t v| 2 m = ∇ 2 t c s , v -3 ∇ t c s , ∇ t v -∇ t c t /|c t |, v ∇ t c s , v .
Finally, we propose an optimal matching algorithm that reduces the distance between the initial curve c 0 and the fiber of the target curve c 1 at each step. The proof of convergence of this algorithm is work in progress. In Chapter 3, we assume that the base manifold M has constant sectional curvature K, and we give a detailed discretization of the Riemannian structure introduced in Chapter 2 that is itself a Riemannian structure on the finite-dimensional manifold M n+1 of discrete curves given by n + 1 points, n ∈ N * . This discrete model could more generally be extended to symmetric spaces. Its tangent space at a given point α = (x 0 , . . . , x n ) is given by

T α M n+1 = {w = (w 0 , . . . , w n ) : w k ∈ T x k M, ∀k}. We associate to each (α, w) ∈ T M n+1 a path of piecewise geodesic curves [0, 1] 2 (s, t) → c w (s, t) ∈ M such that c w (0, k n ) = x k and c w s (0, k n ) = w k for k = 0, . . . , n.
Then we define the scalar product between w and z in terms of the square root velocity representations q w,z (s) := c w,z t (s)/ |c w,z t (s)| of these paths of piecewise-geodesic curves

G n α (w, z) = w 0 , z 0 + 1 n n-1 k=0 ∇ s q w (0, k n ), ∇ s q z (0, k n ) . (1.16)
This is a discrete analog of (1.7) and it does not depend on the choices of c w and c z . Indeed, we can also obtain a discrete analog of (1.4).

MAIN RESULTS

19

Notation. For α = (x 0 , . . . , x n ) ∈ M n+1 we use the notations

τ k = log M x k x k+1 , q k = √ n τ k / |τ k |, v k = τ k /|τ k |, k = 0, . . . , n (1.17)
as well as

w k T = w k , v k v k and w N k = w k -w k T
to refer to the coordinates of the tangential and normal components of a tangent vector w ∈ T α M n+1 .

Proposition 1.7. The scalar product between two tangent vectors w, z ∈ T α M n+1 can also be written

G n α (w, z) = w 0 , z 0 + n-1 k=0 (D τ w) k N , (D τ z) k N + 1 4 (D τ w) k T , (D τ z) k T 1 |τ k | ,
where the map

D τ : T α M n+1 → T α M n+1 , w → D τ w = (D τ w) 0 , . . . , (D τ w) n is defined by (D τ w) k := 1 n ∇ t c w s (0, k n ) = (w k+1 -w k ) T + b -1 k (w k+1 -a k w k ) N ,
and the coefficients a k and b k take the following values depending on the sectional curvature K of the base manifold

M    a k = cosh |τ k |, b k = sinh |τ k |/|τ k |, if K = -1, a k = 1, b k = 1, if K = 0, a k = cos |τ k |, b k = sin |τ k |/|τ k |, if K = +1.
(1.18)

Remark 1.3. In the flat case our definition gives (D τ w) k = w k+1 -w k .

In the non-flat case (K = ±1), when the discretization gets "thinner", i.e. n → ∞ and

|τ k | → 0 while n|τ k | stays bounded for all 0 ≤ k ≤ n, we get (D τ w) k = n→∞ w k+1 -w k + o(1).
The main result of Chapter 3 establishes the convergence of the discrete model toward the continuous model, and is proven in Section 3.4. Definition 1.9. We say that α = (x 0 , . . . , Let s → c(s) be a C 1 -path of C 2 -curves with non vanishing derivative with respect to t. This path can be identified with an element (s,

x n ) ∈ M n+1 is the discretization of size n of c ∈ M when c( k n ) = x k for all k = 0, . . . , n. A path s → α(s) of discrete curves is the discretization of size n of a path of smooth curves s → c(s) when α(s) is the discretization of c(s) for all s.
t) → c(s, t) of C 1,2 ([0, 1] × [0, 1], M ) such that c t = 0. Consider the C 1 -path in M n+1 , s → α(s) = (x 0 (s), . . . , x n (s)), that is the discretization of size n of c. Then 20 CHAPTER 1. INTRODUCTION
there exists a constant λ > 0 that does not depend on c and such that for n large enough,

|E(c) -E n (α)| ≤ λ n (inf |c t |) -1 |c s | 2 2,∞ (1 + |c t | 1,∞ ) 3 ,
where E and E n are the energies with respect to metrics G and G n respectively and where

|c t | 1,∞ := |c t | ∞ + |∇ t c t | ∞ , |c s | 2,∞ := |c s | ∞ + |∇ t c s | ∞ + |∇ 2 t c s | ∞ ,
and

|w| ∞ := sup s,t∈[0,1]
|w(s, t)| denotes the supremum over both s and t of a vector field w along c.

We then give the geodesic equations for the discrete metric G n , in terms of the following notations.

Notation. For any discrete curve

α = (x 0 , . . . , x n ) ∈ M n+1 we define for all 0 ≤ k ≤ n the functions f k , g k : T x k M → T x k M , f k : w → w T + a k w N , g k : w → |q k |(2w T + b k w N ),
where the coefficients a k , b k are defined by (1.18), and denote by f

(-) k , g (-) k : T x k+1 M → T x k M
the maps obtained by post-composition with parallel transport along the geodesic linking x k to x k+1 .

Remark 1.4. When the discretization gets "thinner", we get for any fixed

w ∈ T x k+1 M , f k (w) = w + o(1/n) and g k (w) = |q k |(w + w T ) + o(1/n).
In the flat setting, these are always equalities.

Proposition 1.8 (Discrete geodesic equations).

A path in M n+1 , s → α(s) = (x 0 (s), . . . , x n (s)), is a geodesic for the metric G n if and only if its SRV representation s → x 0 (s), (q k (s)) k verifies the following differential equations

∇ s x 0 + 1 n R 0 + f (-) 0 (R 1 ) + . . . + f (-) 0 • • • • • f (-) n-2 (R n-1 ) = 0, ∇ 2 s q k + 1 n g (-) k R k+1 + f (-) k+1 (R k+2 ) + . . . + f (-) k+1 • • • • • f (-) n-2 (R n-1 ) = 0, (1.19) for all k = 0, . . . , n -1, where R k := R(q k , ∇ s q k )x k . Remark 1.5. If s → c(s) ∈ M is a C 1 path
of smooth curves and s → α(s) ∈ M n+1 the discretization of size n of c, the coefficients of the discrete geodesic equation (1.19) for α converge to the coefficients of the continuous geodesic equation (1.12) for c.

Just as in the continuous case, the discrete geodesic equations are exploited to build the exponential map and characterize the Jacobi fields on M n+1 in Section 3.2.2, allowing us to construct geodesics between "discrete curves" using geodesic shooting.

Finally, a discretization of the quotient shape space is given in Section 3.2.3, with the aim of finding the optimal matching between two discrete curves, i.e. to redistribute the n + 1 points on the target curve so as to minimize the discrete distance to the n + 1 points of the initial curve. We consider elements of the set M n of discrete curves of size n paired up with their underlying shapes -or equations, e.g. that of an optimal interpolation

M n := {(α, c) ∈ M n+1 × S : α ∈ Disc n (c)},
where Disc n (c) denotes the set of elements of M n+1 that are discretizations of smooth curves with shape c, in the sense of Definition 1.9, Disc n (c) := {α ∈ M n+1 : ∃c ∈ π -1 (c), α is the discretization of size n of c}.

To define the horizontal part of a path of discrete curves, we define a discrete analog of horizontality.

Definition 1.10. The discrete vertical and horizontal spaces in α are the following subsets of

T α M n+1 Ver n α := {mv : m = (m k ) k ∈ R n+1 , m 0 = m n = 0}, Hor n α := {h ∈ T α M n+1 : G n (h, mv) = 0 ∀m = (m k ) k ∈ R n+1 , m 0 = m n = 0}. Proposition 1.9 (Discrete horizontal space). Let α ∈ M n+1 and h ∈ T α M n+1 . Then h ∈ Hor n α if and only if it verifies (D τ h) k , v k -4 |τ k | |τ k-1 | (D τ h) k-1 , b -1 k-1 v k + ( 1 4 -b -1 k-1 )λ k-1 v k-1 = 0.
Any tangent vector w ∈ T α M n+1 can be uniquely decomposed into a sum w = w ver + w hor where w ver = mv ∈ Ver n α , w hor = wmv ∈ Hor n α and the

λ k m k+1 -1 + 4 |τ k | |τ k-1 | (b -2 k-1 + λ 2 k-1 ( 1 4 -b -2 k-1 )) m k + |τ k | |τ k-1 | λ k-1 m k-1 = (D τ w) k , v k -4 |τ k | |τ k-1 | b -1 k-1 (D τ w) k-1 , v k + ( 1 4 -b -1 k-1 )λ k-1 (D τ w) k-1 , v k-1 ,
with the notation

λ k := v k+1 , v k .
We then define a discrete analog of the reparameterization action, for a fixed integer p ∈ N * . Set N := np. Definition 1.11. Let (α = (x 0 , . . . , x n ), c) in M n . We call refinement of size N of (α, c) the discretization β = (y 0 , . . . , y n ) ∈ Disc N (c) of c such that y kp = x k for k = 0, . . . , n and the p -1 points {y , kp < < (k + 1)p} are distributed according to arc-length on c between y kp and y (k+1)p , for all k. Definition 1.12. A reparameterization of a discrete curve α = (x 0 , . . . , x n ) ∈ M n+1 is the result of the action of an increasing injection ϕ : {0, . . . , n} → {0, . . . , N } such that ϕ(0) = 0 and ϕ(n

) = N on (α, c) ∈ M n , ϕ ((x k ) k , c) := ((y ϕ(k) ) k , c),
where (y k ) k is the refinement of size N of ((x k ) k , c). The set of such increasing injections is denoted by Inj + (n, N ).

This definition of reparameterization in the discrete case simply boils down to redistributing the n + 1 points on c by choosing among the N + 1 points of the refinement of α, while preserving the order and keeping the extremities fixed.

Definition 1.13. The horizontal part (α hor , c) of a path s → (α(s), c(s)) ∈ M n is defined by (α(s), c(s)) := ϕ(s) (α hor (s), c(s)), ∀s ∈ [0, 1], where ϕ(s) ∈ Inj + (n, N ) verifies for all s ∈ [0, 1] ϕ s (s)(k) = m k (s) |nτ k (s)| ∆ϕ(s)(k), k = 0, . . . , n, (1.20) 1.3. MAIN RESULTS 23 with ∆ϕ(s)(k) = N/2(ϕ(s)(k + 1) -ϕ(s)(k -1)), 1 ≤ k ≤ N -1
and where m = (m k ) k is the norm of the vertical component of α (s) and verifies

λ k m k+1 -1 + 4 |τ k | |τ k-1 | (b -2 k-1 + λ 2 k-1 ( 1 4 -b -2 k-1 )) m k + |τ k | |τ k-1 | λ k-1 m k-1 = ∇ s τ k , v k -4 |τ k | |τ k-1 | × b -1 k-1 ∇ s τ k-1 , v k + ( 1 4 -b -1 k-1 )λ k-1 ∇ s τ k-1 , v k-1 .
(1.21)

Remark 1.6. Equation (1.20) defining the path of "reparameterizations" ϕ is merely a discretization of Equation (1.15). The recurrence relation (1.21) verified by the m k 's translates the fact that m(s)v(s) is the vertical component of α (s).

With these definitions, we are able to apply the optimal matching algorithm (Algorithm 1.3) to discrete curves. Results of simulations with curves in the hyperbolic half-plane M = H 2 , the plane M = R 2 and the sphere M = S 2 are given in Section 3.3.

Chapter 4 is dedicated to radar applications. We represent a locally stationary radar signal by the time series of the Toeplitz covariance matrices of its stationary portions, i.e. by a "discrete curve" in the space T + n of Toeplitz, hermitian, positive definite matrices. Definition 1.14. A Toeplitz matrix is a constant-diagonal matrix. A square Toeplitz matrix of size n ∈ N * is of the form

T =       t n t n+1 • • • t 2n-1 t n-1 . . . . . . . . . . . . . . . . . . t n+1 t 1 • • • t n-1 t n       .
To compare discrete curves in T + n , we equip it with a Riemannian metric. In Section 4.1.1, we first put aside the Toeplitz structure and consider the metric defined through the hessian of the entropy H on the space H + n of hermitian, positive, definite matrices Σ, seen as the set of covariance matrices of centered multivariate (circularly-symmetric) Gausian densities N (0, Σ)

ds 2 = -dΣ * Hess H(Σ)dΣ. (1.22)
We remind that this is the Legendre dual of the Fisher information metric, and therefore defines the same distance on

H + n . CHAPTER 1. INTRODUCTION Proposition 1.10. Let Σ ∈ H + n .
The dual potentials defined in the dual coordinate systems η = 1 2 Σ -1 and θ = Σ,

φ(η) = -n ln(π) -ln(det Σ), ψ(θ) = n ln(πe) + ln(det Σ) = H(Σ),
where H(Σ) is the entropy, define dual hessian metrics

g φ ij (η) = - ∂ 2 φ(η) ∂ ηi ∂η j , g ψ ij (θ) = - ∂ 2 ψ(θ) ∂ θi ∂θ j
which are the same tensor expressed in different coordinate systems and define the same local distance ds 2 φ = ds 2 ψ . They both coincide with the Fisher metric on H + n . In Section 4.1.2, we project metric (1.22) in the submanifold T + n of Toeplitz matrices in H + n , or more precisely in the equivalent coordinate space [START_REF] Verblunsky | On positive harmonic functions[END_REF][50] of the reflection coefficients {(P 0 , µ 0 , . . . , µ n-1 ) ∈ R * + × D n-1 }, and explain how we obtain the very convenient product metric of the Poincaré polydisk,

ds 2 = n dP 0 P 0 2 + n-1 k=1 (n -k) |dµ k | 2 (1 -|µ k | 2 ) 2 .
(1.23) Definition 1.15. We call Poincaré polydisk of size n ∈ N * the product manifold R * + × D n-1 , where D = {z ∈ C, |z| < 1} is the unit complex disk, equipped with metric (1.23). We denote it by R * + × D n-1 . This allows us to represent stationary centered Gaussian signals by elements of the Poincaré polydisk, and to compare locally stationary Gaussian signals by curves in the polydisk. In Section 4.2, we give an example of application of shape analysis to radar signal processing : we use the Riemannian framework detailed in Chapters 2 and 3 to compute the mean of locally stationary helicopter signatures.

Introduction

Computing distances between shapes of open or closed curves is of interest in many applications, from medical imaging to radar signal processing, as soon as one wants to compare, classify or statistically analyze trajectories or the outline of objects. While the shape of an organ or the trajectory of an object on a short distance can be modeled by a curve in the plane R 2 or in the ambient space R 3 , some applications provide curves in an intrinsically 25 CHAPTER 2. RIEMANNIAN FRAMEWORK nonlinear space. Simple examples in positive curvature include trajectories on the sphere where points represent positions on the earth, and a negativelycurved space of interest in signal processing is the hyperbolic half plane, which as we have seen in Section 1.2, coincides with the statistical manifold of Gaussian densities.

We consider open oriented curves in a Riemannian manifold M , that we represent by smooth immersions, i.e. curves with velocity that doesn't vanish. Their set is an open submanifold of the Fréchet manifold C ∞ ([0, 1], M ) (see [START_REF] Peter W Michor | Manifolds of differentiable mappings[END_REF], Theorem 10.4.)

M = Imm([0, 1], M ) = {c ∈ C ∞ ([0, 1], M ), c (t) = 0 ∀t ∈ [0, 1]}.
To compare or average elements of this space, we equip M with a Riemannian metric, i.e. we locally define a scalar product G on its tangent space T M. Elements w, z ∈ T c M of the tangent space in c ∈ M are infinitesimal deformations of c and can be seen as vector fields along the curve c in Mthis results from the so called "Exponential law" for smooth functions, see e.g. [START_REF] Kriegl | Aspects of the theory of infinite dimensional manifolds[END_REF], Theorem 5.6. -

T c M = {w ∈ C ∞ ([0, 1], T M ) : w(t) ∈ T c(t) M ∀t ∈ [0, 1]}.
An interesting property from the point of applications is reparametrization invariance, that is that the metric be the same at all points of M representing curves that are identical modulo reparametrization. Two curves are identical modulo reparametrization when they pass through the same points of M but at different speeds. Reparametrizations are represented by increasing diffeomorphisms ϕ : [0, 1] → [0, 1] (so that they preserve the end points of the curves), and their set is denoted by Diff

+ ([0, 1]). The metric G is reparametrization invariant if the action of Diff + ([0, 1]) is isometric for G G c•ϕ (w • ϕ, z • ϕ) = G c (w, z), (2.1)
for any c ∈ M, w, z ∈ T c M and ϕ ∈ Diff + ([0, 1]). This is often called the equivariance property, and it guarantees that the induced distance between two curves c 0 and c 1 does not change if we reparametrize them by the same

diffeomorphism ϕ d(c o • ϕ, c 1 • ϕ) = d(c 0 , c 1 ).
What's more, a reparametrization invariant metric on the space of curves induces a Riemannian metric G on the quotient space

S = Imm([0, 1], M )/Diff + ([0, 1]).

EXTENSION OF THE SRV FRAMEWORK

27

classically interpreted as the space of shapes or unparameterized curves, such that the natural projection π : M → S is a Riemannian submersion. With this definition, a shape is the equivalence class of all the curves that are identical modulo a change of parameterization. Under some restrictions on the space of curves, the projection π defines a principal bundle, and the horizontal geodesics of the total space -that is, those that are G-orthogonal to the fibers -project onto geodesics for G on the quotient space. The geodesic distances d on M and d on S are linked by

d (c 0 , c 1 ) = inf d (c 0 , c 1 • ϕ) | ϕ ∈ Diff + ([0, 1]) ,
where c 0 and c 1 denote the shapes of two given curves c 0 and c 1 , and d verifies the stronger property

d(c 0 • φ, c 1 • ψ) = d(c 0 , c 1 ),
for any reparametrizations φ, ψ ∈ Diff + ([0, 1]). This quotient structure motivates the choice of a reparametrization invariant metric on M.

In the following section, we introduce our metric as the pullback of a quite natural metric on the tangent bundle T M, and show that it induces a fiber bundle structure over the manifold M seen as the set of starting points of the curves. We give the induced geodesic distance and highlight the difference with respect to the generalization of the SRV framework introduced in [START_REF] Zhang | Video-based action recognition using rate-invariant analysis of covariance trajectories[END_REF]. In Section 2.3, we derive the geodesic equation and exploit them to solve the initial and boundary value problems. Finally, we address the quotient structure of the shape space in Section 2.4, and introduce an optimal matching algorithm.

Extension of the SRV framework to curves in a manifold

As mentioned in the short review of related work in Chapter 1, Riemannian metrics on the space of curves lying in the Euclidean space R n , and especially closed curves c : S 1 → R n (S 1 is the circle), have been widely studied. One first-order Sobolev metric where different weights are given to the tangential and normal parts of the derivative has proved particularly interesting for the applications ( [START_REF] Laga | Landmark-free statistical analysis of the shape of plant leaves[END_REF], [START_REF] Su | Statistical analysis of trajectories on riemannian manifolds: bird migration, hurricane tracking and video surveillance[END_REF])

G c (w, z) = D w N , D z N + 1 4 D w T , D z T d . (2.2) CHAPTER 2. RIEMANNIAN FRAMEWORK In that case c is a curve in R n , •, •
G a,b c (w, z) = a 2 D w N , D z N + b 2 D w T , D z T d ,
where parameters a and b respectively control the degree of bending and stretching of the curve. Srivastava et al. introduced in [START_REF] Srivastava | Shape analysis of elastic curves in euclidean spaces[END_REF] a convenient framework to study the case where a = 1 and b = 1/2 by showing that metric 2.2 could be obtained by pullback of the L 2 -metric via a simple transformation called the Square Root Velocity Function (SRVF), which associates to each curve its velocity renormalized by the square root of its norm. The SRV framework has known several extensions, to more general metrics -the general elastic metric G a,b with weights a and b satisfying 4b 2 ≥ a 2 [10]and to larger spaces of curves. Extension to curves in a Lie group can be achieved using translations [START_REF] Celledoni | Shape analysis on lie groups with applications in computer animation[END_REF], and to curves in a nonlinear manifold using parallel transport, so as to move the computations to the tangent space to the origins of the curves [START_REF] Su | Statistical analysis of trajectories on riemannian manifolds: bird migration, hurricane tracking and video surveillance[END_REF], [START_REF] Le Brigant | Reparameterization invariant metric on the space of curves[END_REF], [START_REF] Zhang | Video-based action recognition using rate-invariant analysis of covariance trajectories[END_REF]. In [START_REF] Le Brigant | Reparameterization invariant metric on the space of curves[END_REF] we considered the general elastic metric G a,b , but no Riemannian framework was given. In [START_REF] Zhang | Video-based action recognition using rate-invariant analysis of covariance trajectories[END_REF], a Riemannian framework is given for the case a = 1, b = 1/2, and the geodesic equations are derived. Here we also restrict to this particular choice of coefficients a and b for simplicity, but we propose another generalization of the SRV framework to manifold-valued curves. Instead of encoding the information of each curve within a tangent space at a single point as in [START_REF] Le Brigant | Reparameterization invariant metric on the space of curves[END_REF] and [START_REF] Zhang | Video-based action recognition using rate-invariant analysis of covariance trajectories[END_REF] using parallel transport, the distance is computed in the manifold itself which enables us to be more directly dependent on its geometry. Intuitively, the data of each curve is no longer concentrated at any one point, and so the energy of the deformation between two curves takes into account the curvature of the manifold along the entire "deformation surface", not just along the path traversed by the starting point of the curve.

Notations

Let (M, •, • ) be a Riemannian manifold. We first introduce a few notations. The norm associated to the Riemannian metric •, • is denoted by | • | and the Levi-Civita connection by ∇. We denote by exp M x : T x M → M the exponential map on M at point x ∈ M and by log M

x : M → T x M its inverse map. To avoid confusions, we will always denote by s the parameter of paths in the space of curves M and by t the parameter of a curve in M . If t → c(t) is a curve in M and t → w(t) ∈ T c(t) M a vector field along c, we denote by c t := dc/dt = c the derivative of c with respect to t and by ∇ t w := ∇ ct w, ∇ 2 t w := ∇ ct ∇ ct w the first and second order covariant derivatives of w along c. We use various notations depending on the context to denote parallel transport according to connection ∇. If u ∈ T c(t 1 ) M is a tangent vector to M in c(t 1 ), the parallel transport of u from c(t 1 ) to c(t 2 ) along c is denoted by P t 1 ,t 2 c (u), or when there is no ambiguity on the choice of the curve c, u t 1 ,t 2 , or even u to lighten notations in some cases. We associate to each curve c its renormalized speed vector field v := c /|c |, and to each vector field t → w(t) along c, its tangential and normal components w T := w, v v and w N := ww T . Finally, we identify a path of curves [0, 1] s → c(s) ∈ M with the function of two variables [0, 1] × [0, 1] (s, t) → c(s, t) ∈ M and denote by c s := ∂c/∂ s and c t := ∂c/∂ t its partial derivatives with respect to s and t, while ∇ s = ∇ ∂c/∂s and ∇ t = ∇ ∂c/∂t denote partial covariant derivatives.

Our metric on the space of curves

Let c : [0, 1] → M be a curve in M and w, z ∈ T c M two infinitesimal deformations. We consider the following first-order Sobolev metric on M

G c (w, z) = w(0), z(0) + ∇ w N , ∇ z N + 1 4 ∇ w T , ∇ z T d ,
where we integrate according to arc length d = |c (t)|dt, ∇ w = ∇ c w/|c | is the covariant derivative of w according to arc length, and ∇ w T = ∇ w, v v and ∇ w N = ∇ w -∇ w T are its tangential and normal components respectively. If M is a flat Euclidean space, we obtain the metric 2.2 studied in [START_REF] Srivastava | Shape analysis of elastic curves in euclidean spaces[END_REF], with an added term involving the origins. Without this extra term, the bilinear form G is not definite since it vanishes if w or z is covariantly constant along c. Here we show that G can be obtained as the pullback of a very natural metric Ĝ on the tangent bundle T M. We represent each curve c ∈ M by the pair formed by its starting point and its speed vector field renormalized by the square root of its norm, via the bijection

M → M × T M, c → x(c) := c(0), q(c) := c |c | .
When there is no ambiguity, we will use the lighter notations x := x(c) and q := q(c) to denote the SRV coordinates of a curve c. The inverse of this CHAPTER 2. RIEMANNIAN FRAMEWORK function is simply given by M × T M (x, q) → π M (q) ∈ M, if π M is the canonical projection T M → M that associates to each tangent vector its base point. In order to define Ĝ, we introduce the following projections from T T M to T M . Let ξ ∈ T (p,u) T M be an infinitesimal deformation of the pair formed by a point p ∈ M and a tangent vector u attached to p, and s → (x(s), U (s)) be a curve in T M that passes through (p, u) at time 0 at speed ξ, i.e. that deforms the pair (p, u) in the direction given by ξ. Then we define the vertical and horizontal projections linked to the fiber bundle structure π M :

T M → M , vp (p,u) : T (p,u) T M → T p M, ξ → ξ v := ∇ x (0) U, hp (p,u) : T (p,u) T M → T p M, ξ → ξ h := x (0).
The horizontal projection, which is simply the differential of the projection π M , corresponds to the way ξ "moves the base point x" and vertical projection to the way it "linearly moves u". They live in the tangent bundle T M and are not to be confused with the horizontal and vertical parts for this fiber bundle, which live in the double tangent bundle T T M . Furthermore, let us point out that according to these definitions, the very natural Sasaki metric ( [START_REF] Sasaki | On the differential geometry of tangent bundles of riemannian manifolds[END_REF], [START_REF] Sasaki | On the differential geometry of tangent bundles of riemannian manifolds ii[END_REF]) on the tangent bundle T M can be written g S (p,u) (ξ, η) = ξ h , η h + ξ v , η v . Now we can define the metric that we put on T M. Let us consider a curve c ∈ M, an infinitesimal deformation w ∈ T c M of c, and ξ, η ∈ T (c,w) T M infinitesimal deformations of the pair (c, w). Then for all time t ∈ [0, 1], ξ(t) and η(t) belong to T (c(t),w(t)) T M . We define

Ĝ(c,w) (ξ, η) = ξ(0) h , η(0) h + 1 0 ξ(t) v , η(t) v dt,
where ξ(t) h = hp(ξ(t)) ∈ T M and ξ(t) v = vp(ξ(t)) ∈ T M are the horizontal and vertical projections of ξ(t) for all t. Then we have the following result.

Proposition 2.1. The metric G on the space of curves M can be obtained as pullback of the metric Ĝ by the square root velocity function, that is

G c (w, z) = Ĝq(c) (T c q(w), T c q(z)) = w(0), z(0) + 1 0 ∇ w(t) q(c), ∇ z(t) q(c) dt,
for any curve c ∈ M and w, z ∈ T c M, where T c q is the differential of the function q : M → T M, c → q(c) at c. Remark 2.1. The term ∇ w q(c) denotes the covariant derivative of the vector field q(c) in the direction given by the vector field w. More precisely, if s → c(s, •) is a path of curves verifying c(0, t) = c(t) and c s (0, t) = w(t) for all t ∈ [0, 1], then

∇ w(t) q(c) = ∇ s | s=0 q c(s, t) , t ∈ [0, 1],
where ∇ s denotes the covariant derivative along the curve s → c(s, t).

Proof of Proposition 1. For any t ∈ [0, 1], the horizontal and vertical projections of T c q(w)(t) are given by T c q(w)(t) h = w(t) and T c q(w) v = ∇ w(t) q(c).

To prove this proposition, we just need to compute the latter. Let s → c(s, •) be a path in M such that c(0,

•) = c and c s (0, •) = w . Then ∇ w(t) q(c) = 1 |c | 1/2 ∇ w c + w |c | -1/2 c = 1 |c t | 1/2 ∇ s c t + ∂ s c t , c t -1/4 c t = 1 |c t | 1/2 ∇ t c s -1 2 c t , c t -5/4 ∇ s c t , c t c t = |c | 1/2 (∇ w) N + 1 2 (∇ w)
T , where we used twice the inversion ∇ s c t = ∇ t c s .

Fiber bundle over the starting points

The special role played by the starting point in the metric G induces a fiber bundle structure, where the base space is the manifold M , seen as the set of starting points of the curves, and the fibers are composed of the curves which have the same origin. The projection is then

π ( * ) : M → M, c → c(0).

It induces a decomposition of the tangent bundle in vertical and horizontal bundles

V ( * ) c M = ker T π ( * ) = { w ∈ T c M | w(0) = 0 } , H ( * ) c M = V ( * ) c M ⊥ G .
Proposition 2.2. We have the usual decomposition

T M = V ( * ) M ⊥ ⊕ H ( * ) M, the horizontal bundle H ( * )
c M consists of parallel vector fields along c, and π ( * ) is a Riemannian submersion for (M, G) and (M, •, • ).

CHAPTER 2. RIEMANNIAN FRAMEWORK

Proof. Let w ∈ T c M be a tangent vector. Consider w h the parallel vector field along c with initial value w h (0) = w(0). It is a horizontal vector for the previously described fiber bundle, since its vanishing covariant derivative along c assures that for any vertical vector l we have G c (w h , l) = 0. The difference w v = ww h between those two vectors has initial value 0 and so it is a vertical vector, which gives a decomposition of w into a horizontal part w h and a vertical part w v . The definition of H ( * ) M as the orthogonal complement of V ( * ) M guaranties that their sum is direct. Since T c π ( * ) (w) = w(0) and w(0) = w h (0) for all w ∈ T c M, the scalar product between the horizontal parts of two tangent vectors w, z is given by

G c (w h , z h ) = w h (0),z h (0) = w(0), z(0) = T c π ( * ) (w h ), T c π ( * ) (z h ) ,
which proves that T π ( * ) : H ( * ) M → T M is an isometry and so that π ( * ) is a Riemannian submersion.

Induced distance on the space of curves

Here we give an expression of the geodesic distance induced by the metric G. We show that it can be written similarly to the product distance given in [START_REF] Le Brigant | Reparameterization invariant metric on the space of curves[END_REF] and [START_REF] Zhang | Video-based action recognition using rate-invariant analysis of covariance trajectories[END_REF], with an added curvature term. Let us consider two curves c 0 , c 1 ∈ M, and a path of curves s → c(s,

•) linking them in M c(0, t) = c 0 (t), c(1, t) = c 1 (t), for all t ∈ [0, 1]. We denote by (x, q) : [0, 1] → M × T M, s → x(s) := x(c(s)) := c(s, 0), q(s) := q(c(s)) := c t (s)/|c t (s)| 1/2 ,
the image of this path of curves by the SRVF. Note that q is a vector field along the "surface" c in M . Let now q be the "raising" of q in the tangent space T c(0,0) M defined by the following parallel transport of q q(s, t)

= P s,0 c(•,0) • P t,0 c(s,•) (q(s, t)) . (2.3)
Notice that q is a "surface" in a vector space, as illustrated in Figure 2.1. Lastly, we introduce a vector field (a, τ ) → ω s,t (a, τ ) in M , which parallel translates q(s, t) along c(s, •) to its origin, then along c(•, 0) and back down again, as shown in Figure 2.1. More precisely

ω s,t (a, τ ) = P 0,τ c(a,•) • P s,a c(•,0) • P t,0 c(s,•) (q(s, t)) , (2.4)
for all b, s. That way the quantity ∇ a ω s,t (s, t) measures the holonomy along the rectangle of infinitesimal width shown in Figure 2.1.

Proposition 2.3. With the above notations, the geodesic distance induced by the Riemannian metric G between two curves c 0 and c 1 on the space M = Imm([0, 1], M ) of parameterized curves is given by

d(c 0 , c 1 ) = inf c(0)=c 0 ,c(1)=c 1 1 0 |x s (s)| 2 + 1 0 |∇ s q(s, t)| 2 dt ds. (2.5)
It can also be written as a function of the "raising" q of q in the tangent space T c 0 (0) M defined by 2.3,

d(c 0 , c 1 ) = inf c(0)=c 0 ,c(1)=c 1 1 0 |x s (s)| 2 + 1 0 |q s (s, t) + Ω(s, t)| 2 dt ds, (2.6)
where Ω is a curvature term measuring the holonomy along a rectangle of infinitesimal width

Ω(s, t) = P s,0 c(•,0) • P t,0 c(s,•) ∇ a ω s,t (s, t) = P s,0 c(•,0) t 0 P τ,0 c(s,•) R(c τ , c s )P t,τ c(s,•) q(s, t) dτ ,
and R is the curvature tensor of the manifold M and ω s,t is defined by 2.4. Remark 2.2. The second expression 2.6 highlights the difference with respect to the distance given in [START_REF] Le Brigant | Reparameterization invariant metric on the space of curves[END_REF] and [START_REF] Zhang | Video-based action recognition using rate-invariant analysis of covariance trajectories[END_REF]. In the first term under the square CHAPTER 2. RIEMANNIAN FRAMEWORK root we can see the velocity vector of the curve x(•) = c(•, 0) linking the two origins, and in the second the velocity vector of the curve q linking the TSRVF-images of the curves -Transported Square Root Velocity Function, as introduced by Su et al. in [START_REF] Su | Statistical analysis of trajectories on riemannian manifolds: bird migration, hurricane tracking and video surveillance[END_REF]. If instead we equip the tangent bundle TM with the metric

Ĝ w (ξ, ξ) = |ξ(0) h | 2 + 1 0 ξ(t) v - t 0 P τ,t c R(c , ξ h )P t,τ c q(t) dτ 2 dt,
for w ∈ T M and ξ, η ∈ T w T M, then the curvature term Ω vanishes and the geodesic distance on M becomes

d (c 0 , c 1 ) = inf c(0)=c 0 ,c(1)=c 1 1 0 |x s (s)| 2 + qs (s) 2 L 2 , ds, (2.7)
which corresponds exactly to the geodesic distance introduced by Zhang et al. in [START_REF] Zhang | Video-based action recognition using rate-invariant analysis of covariance trajectories[END_REF] on the space C = ∪ p∈M L 2 ([0, 1], T p M ). The difference between the two distances 2.5 and 2.7 resides in the curvature term Ω, which measures the holonomy along the rectangle of infinitesimal width shown in Figure 2.1, and arises from the fact that in the first one, we compute the distance in the manifold, whereas in the second, it is computed in the tangent space to one of the origins of the curves. Therefore, the first one takes more directly into account the geometry of the manifold between the two curves under comparison, since it reads the information directly in the manifold itself.

Remark 2.3. Let us briefly consider the flat case : if the manifold M is flat, e.g. M = R n , then the two distances 2.5 and 2.7 coincide. If two curves c 0 and c 1 in R n have the same starting point p, the first summand under the square root vanishes and the distance becomes the L 2 -distance between the two renormalized speed vector fields q 0 = q(c 0 ) and q 1 = q(c 1 ). If two R n -valued curves differ only by a translation, then the distance is simply the distance between their origins.

Remark 2.4. Note that this distance is only local in general, that is, only works for curves that are "close enough". Indeed, if we consider two curves c 1 , c 2 in M = R 2 with the same origin, the distance between them is the length of the L 2 geodesic between their SRV representations q 1 and q 2 in C ∞ ([0, 1], R 2 \{0}). If the minimizing geodesic between those two (in

C ∞ ([0, 1], R 2
)) passes through 0, then there is no minimizing geodesic between q 1 and q 2 in C ∞ ([0, 1], R 2 \{0}).

Proof of Proposition 3. Since G is defined by pullback of Ĝ by the SRVF, we know that the lengths of c in M and of (x = x(c), q = q(c)) in TM are equal and so that

d(c 0 , c 1 ) = inf c(0)=c 0 ,c(1)=c 1 1 0 Ĝ ∂ ∂s q(c(s)), ∂ ∂s q(c(s)) ds, with Ĝ ∂ ∂s q(c(s)), ∂ ∂s q(c(s)) = |x s (s)| 2 + 1 0 |∇ s q(s, t)| 2 dt.
To obtain the second expression of this distance we need to express ∇ s q as a function of the derivative qs . Let us fix t ∈ [0, 1], and consider the vector field ν along the surface (s, τ ) → c(s, τ ) that is parallel along all curves c(s, •) and takes value ν(s, t) = q(s, t) in τ = t for any s ∈ [0, 1], that is

ν(s, τ ) = P t,τ c(s,•) (q(s, t)) ,
for all s, τ ∈ [START_REF](s) of the horizontal part of α(s) is chosen to be the point y (s) whose value is closest to kp. Now we can go back to our initial problem[END_REF][START_REF] Amari | Differential geometry in statistical inference[END_REF]. With this definition we have ∇ s ν(s, t) = ∇ s q(s, t). Since ν(•, 0) : s → P t,0 c(s,•) (q(s, t)) is a vector field along c(•, 0), we can write ∇ s ν(s, 0) = ∇ s P t,0 c(s,•) q(s, t)

= P 0,s c(•,0) ∂ ∂s P s,0 c(•,0) • P t,0 c(s,•) (q(s, t)) = P 0,s c(•,0) qs (s, t).
Noticing that we additionally have ∇ τ ν(s, τ ) = 0 for all s, τ ∈ [0, 1], and using ∇ τ ∇ s ν = ∇ s ∇ τ ν + R(c τ , c s )ν, the covariant derivative in τ = t can be written

∇ s ν(s, t) = P 0,t c(s,•) (∇ s ν(s, 0)) + t 0 P τ,t c(s,•) (∇ τ ∇ s ν(s, τ )) dτ = P 0,t c(s,•) • P 0,s c(•,0) (q s (s, t)) + t 0 P τ,t c(s,•) R(c τ , c s )P t,τ c(s,•) q(s, t) dτ.
(2.8) Now let us fix s ∈ [0, 1] as well. Notice that the vector field ω s,t defined above as ω s,t (a, τ

) = P 0,τ c(a,•) • P s,a c(•,0) • P t,0 c(s,•) (q(s, t)) verifies ∇ τ ω s,t (s, τ ) = 0 ∀τ ∈ [0, 1], (2.9
)

∇ a ω s,t (a, 0) = 0 ∀a ∈ [0, 1].
(2.10)

CHAPTER 2. RIEMANNIAN FRAMEWORK

Note that unlike ν, we do not have ∇ a ω s,t (s, t) = ∇ s q(s, t) because ω s,t (a, t) = q(a, t) is only true for a = s. Using Equations 2.9 and 2.10 we get ∇ a ω s,t (s, t) = P 0,t c(s,•) ∇ a ω s,t (s, 0) +

t 0 P τ,t c(s,•) ∇ τ ∇ a ω s,t (s, τ ) dτ = t 0 P τ,t c(s,•) ∇ a ∇ τ ω s,t (s, τ ) + R(c τ , c s )ω s,t (s, τ ) dτ, = t 0 P τ,t c(s,•) R(c τ , c s )P t,τ c(s,•) q(s, t) dτ,
which is the same integral as the one in 2.8. Finally, since

|∇ s q(s, t)| = |∇ s ν(s, t)| = |P s,0 c(•,0) • P t,0 c(s,•) (∇ s ν(s, t)) |, we obtain |∇ s q(s, t)| = |q s (s, t) + P s,0 c(•,0) • P t,0 c(s,•) ∇ a ω s,t (s, t) | = qs (s, t) + P s,0 c(•,0) t 0 P τ,0 c(s,•) R(c τ , c s )P t,τ c(s,•) q(s, t) dτ
which gives Equation 2.6 and completes the proof.

Computing geodesics between curves

Geodesic equation

To be able to compute the distance given by 2.5 between two curves, we first need to compute the optimal deformation s → c(s, •) from one to the other.

That is, we need to characterize the geodesics of M for our metric. In order to do so, taking inspiration from [START_REF] Zhang | Video-based action recognition using rate-invariant analysis of covariance trajectories[END_REF], we use the variational principle, i.e. we search for the critical points of the energy functional E :

C ∞ ([0, 1], M) → R + , E(c) = 1 0 |x (s)| 2 + 1 0 |∇ s q(s, t)| 2 dt ds. (2.11)
The geodesics of M for G can be characterized in terms of their square root velocity representation as follows.

Proposition 2.4 (Geodesic equation). Let [0, 1] s → c(s) ∈ M be a path of curves. It is a geodesic of M if and only if its SRV representation s → (x(s), q(s)) verifies the following equations

∇ s x s (s) + r(s, 0) = 0, ∀s (2.12a) ∇ 2 s q(s, t) + |q(s, t)| r(s, t) + r(s, t) T = 0, ∀t, s, (2.12b)
where the vector field r is given by r(s, t) = 1 t R(q, ∇ s q)c s (s, τ ) τ,t dτ, and w(s, τ ) τ,t denotes the parallel transport of the vector field w(s, •) along c(s, •) from c(s, τ ) to c(s, t).

Remark 2.5. In the flat case M = R d , the curvature term r vanishes and we obtain ∇ s x s (s) = 0, ∇ 2 s q(s, t) = 0 for all s and t. We recover the fact that the geodesic between two curves (x 0 , q 0 ) and (x 1 , q 1 ) in the SRV representation space M × T M is composed of a straight line s → x(s) and an L 2 -geodesic s → q(s, •). This is illustrated in simulations of Section 3.3.

Proof. The path c is a geodesic if and only if it is a critical point of the energy functional (2.11). Let a → ĉ(a, •, •), a ∈ (-ε, ε), be a proper variation of the path s → c(s, •), meaning that it coincides with c in a = 0, and it preserves its end points

ĉ(0, s, t) = c(s, t) ∀s, t, ĉa (a, 0, t) = 0 ∀a, t, ĉa (a, 1, t) = 0 ∀a, t.
Then c is a geodesic of M if and only if d da a=0 E(ĉ(a, •, •)) = 0 for any proper variation ĉ. If we denote by E(a) = E(ĉ(a, •, •)), for a ∈ (-ε, ε), the energy of a proper variation ĉ, then we have

E(a) = 1 2 1 0 |ĉ s (a, s, 0)| 2 ds + 1 0 |∇ s q(a, s, t)| 2 dt ds,
where q = ĉt / |ĉ t | is the SRV representation of ĉ. Its derivative is given by

E (a) = 1 0
∇ a ĉs (a, s, 0), ĉs (a, s, 0) ds

+ 1 0 1 0 ∇ a ∇ s q(a, s, t), ∇ s q(a, s, t) dt ds.
Considering that the variation preserves the end points, integration by parts gives ∇ s ∇ a q, ∇ s q ds = -

1 0 ∇ s ∇ s q, ∇ a q ds, CHAPTER 2. RIEMANNIAN FRAMEWORK
and so the derivative E (a) can be written

- 1 0 ∇ s ĉs , ĉa | t=0 ds + 1 0 1 0 R(c a , c s )q + ∇ s ∇ a q, ∇ s q dt ds = - 1 0 ∇ s ĉs , ĉa | t=0 ds - 1 0 1 0
R(q, ∇ s q)ĉ s , ĉa + ∇ s ∇ s q, ∇ a q dt ds.

This quantity has to vanish in a = 0 for all proper variations ĉ

∇ s c s | t=0 , ĉa | a=0,t=0 ds + 1 0 1 0 R(q, ∇ s q)c s , ĉa | a=0 + ∇ s ∇ s q, ∇ a q| a=0 dt ds = 0.
We cannot yield any conclusions at this point, because the derivatives in a ĉa (0, s, t) and ∇ a q(0, s, t) cannot be chosen independently, since q is not any vector field along ĉ but its image via the square root velocity function.

Computing the covariant derivative of q = ĉt /|ĉ t | 1/2 according to a gives

∇ a q = |ĉ t | -1/2 (∇ a ĉt -1 2 ∇ a ĉt T ), and projecting both sides on v = c t /|c t | results in ∇ a qT = 1 2 |q| -1 ∇ a ĉt T .
We deduce

∇ a ĉt = |q| ∇ a q + ∇ a qT ,
and since ∇ t ĉa = ∇ a ĉt , we can express the variation ĉa as follows ĉa (0, s, t) = ĉa (0, s, 0) 0,t + t 0 |q(0, s, τ )| ∇ a q(0, s, τ ) + ∇ a q(0, s, τ ) T τ,t dτ.

Inserting this expression in the derivative of the energy we obtain the following, where we omit to write that the variations ĉ and q are always taken in a = 0 for the sake of readability,

1 0 ∇ s c s (s, 0), ĉa (s, 0) ds + 1 0 1 0 R(q, ∇ s q)c s (s, t), ĉa (s, 0) 0,t dtds + 1 0 1 0 R(q, ∇ s q)c s (s, t), t 0 |q(s, τ )| ∇ a q(s, τ ) + ∇ a q(s, τ ) T τ,t dτ dtds + 1 0 1 0 ∇ s ∇ s q(s, t), ∇ a q(s, t) dtds = 1 0 ∇ s c s (s, 0) + 1 0 R(q, ∇ s q)c s (s, τ ) τ,0 dτ , ĉa (s, 0) ds + 1 0 1 0 1 t R(q, ∇ s q)c s (s, τ ) τ,t , |q(s, t)|(∇ a q(s, t) + ∇ a q(s, t) T ) dτ dtds + 1 0 1 0 ∇ s ∇ s q(s, t), ∇ a q(s, t) dtds 2.3. COMPUTING GEODESICS BETWEEN CURVES 39 = 1 0 ∇ s c s (s, 0) + r(s, 0), ĉa (s, 0) ds + 1 0 1 0 ∇ s ∇ s q(s, t) + |q(s, t)|(r(s, t) + r(s, t) T ), ∇ a q(s, t) dtds = 0,
with the previously given definition of r. Since the variations ĉa (0, s, 0) and ∇ a q(0, s, t) can be chosen independently and take any value for all s and all t, we obtain the desired equations.

Exponential map

Now that we have the geodesic equation, we are able to describe an algorithm which allows us to compute the geodesic s → c(s, •) starting from a point c ∈ M at speed w ∈ T c M. This amounts to finding the optimal deformation of the curve c in the direction of the vector field w according to our metric. We initialize this path s → c(s, •) by setting c(0, •) = c and c s (0, •) = w, and we propagate it using iterations of fixed step ε > 0. The aim is, given c(s, •) and c s (s, •), to deduce c(s + ε, •) and c s (s + ε, •). The first is obtained by following the exponential map on the manifold M c(s + ε, t) = exp M c(s,t) (εc s (s, t)) , for all t ∈ [0, 1] and the second requires the computation of the acceleration

∇ s c s (s, •) c s (s + ε, t) = [c s (s, t) + ε∇ s c s (s, t)] s,s+ε ,
for all t ∈ [0, 1] where we use the notation w(s) s,s+ε = P s,s+ε c (w(s)) for the parallel transport of a vector field s → w(s) along a curve s → c(s) in M . If we assume that at time s we have c(s, •) and c s (s, •) at our disposal, then we can estimate c t (s, •) and ∇ t c s (s, •), and deduce q(s, •) = c t (s, •)/ |c t | as well as

∇ s q(s, •) = ∇ s c t |c t | (s, •) - 1 2
∇ s c t , c t |c t | 5/2 c t (s, •), (2.13)
using the fact that ∇ s c t = ∇ t c s . The variation ∇ s c s (s, •) can then be computed in the following way

∇ s c s (s, t) = ∇ s c s (s, 0) 0,t + t 0 ∇ 2 s c t (s, τ ) + R(c t , c s )c s (s, τ ) τ,t dτ (2.14)
for all t ∈ [0, 1], where ∇ s c s (s, 0) is given by equation 2.12a, the second order variation ∇ 2 s c t (s, •) is given by

∇ 2 s c t = |c t | 1/2 ∇ 2 s q + ∇ t c s , c t |c t | 2 ∇ t c s + ∇ 2 s q, c t |c t | 3/2 - 3 2 ∇ t c s , c t 2 |c t | 4 + |∇ t c s | 2 |c t | 2 c t , (2.15)
and ∇ 2 s q can be computed via equation 2.12b. Algorithm 2.1 (Exponential Map in M). Input : (c 0 , w) ∈ T M. Initialization : Set c(0, t) = c 0 (t) and c s (0, t) = w(t) for all t ∈ [0, 1]. Heredity : For k = 0, . . . , m -1, set s = kε with ε = 1/m and 1. compute for all t c t (s, t) = lim

δ→0 1 δ log M c(s,t) c(s, t + δ), ∇ t c s (s, t) = lim δ→0 1 δ c s (s, t + δ) t+δ,t -c s (s, t) ,
and compute q(s, t) = 1 √ |ct| c t (s, t) and ∇ s q(s, t) using Equation 2.13;

2. compute r(s, t) = 1 t R(q, ∇ s q)c s (s, τ ) τ,t dτ, and

∇ 2 s q(s, t) = -|q(s, t)| r(s, t) + r(s, t) T ,
and deduce ∇ 2 s c t (s, t) using equation 2.15 for all t ∈ [0, 1]; 3. initialize ∇ s c s (s, 0) = -r(s, 0) and compute ∇ s c s (s, •) using 2.14;

finally, for all

t ∈ [0, 1], set c(s + ε, t) = exp M c(s,t) (εc s (s, t)) , c s (s + ε, t) = [ c s (s, t) + ε∇ s c s (s, t) ] s,s+ε .
Output : c = exp M c 0 w. The last step needed to compute the optimal deformation between two curves c 0 and c 1 is to find the appropriate initial speed u, that is the one that will connect c 0 to c 1 . Since we do not have an explicit expression for this appropriate initial speed, we compute it iteratively using geodesic shooting.

Geodesic shooting and Jacobi fields

The aim of geodesic shooting is to compute the geodesic linking two points p 0 and p 1 of a manifold N , knowing the exponential map exp N . More precisely, the goal is to iteratively find the initial speed w p 0 ,p 1 such that exp N p 0 (w p 0 ,p 1 ) = p 1 .

An initial speed vector w ∈ T p 0 N is chosen, and is iteratively updated after evaluating the gap between the point p = exp N p 0 w obtained by taking the exponential map at point p 0 in w -that is, by "shooting" from p 0 in the direction w -and the target point p 1 . Assuming that the current point p is "not too far" from the target point p 1 , and that there exists a geodesic linking p 0 to p 1 , we can consider that the gap between p and p 1 is the extremity of a Jacobi field J : [0, 1] → N , s → J(s), in the sense that it measures the variation between the geodesics s → exp N p 0 (sw) and s → exp N p 0 (sw p 0 ,p 1 ). Since both geodesics start at p 0 , this Jacobi field has value J(0) = 0 in 0. Then, the current speed vector can be corrected by In the context of geodesic shooting between two curves c 0 and c 1 in M, the speed vector w can be initialized using the L 2 logarithm map, the inverse of the exponential map for the L 2 -metric (these maps are simply obtained by post-composition of mappings with the finite-dimensional maps exp M and log M ). That is, we set w = log L 2 c 0 (c 1 ). The L 2 logarithm map also allows us to approximate the gap between the current point and the target point. This amounts to minimizing the functional F (w) = dist L 2 (exp M c 0 (w), c 1 ). We summarize as follows. 2. estimate the gap j = log L 2 c (c 1 ); 3. if j L 2 > δ, set J(1) = j and w ← w + ∇ s J(0) where ∇ s J(0) = φ -1 (J( 1)) is computed using Algorithm 2.3, and go back to the first step; else, stop.

w ← w + ∇ s J(0),
Output : c approximation of the geodesic linking c 0 and c 1 .

The function φ : T c(0) M → T c(1) M is a linear bijection that maps the last value J(1) of a Jacobi field with initial value J(0) = 0 to the initial speed ∇ s J(0), and can be deduced from Algorithm 2.3. To find the inverse of this function, we consider the image of a basis of the tangent vector space T c(0) M. Now, let us characterize the Jacobi fields of M to obtain the function φ. A Jacobi field is a vector field that describes the way geodesics spread out on a manifold. Consider a → c(a, •, •), a ∈ (-ε, ε), a family of geodesics in M, that is for each a ∈ (-ε, ε), [0, 1] s → c(a, s, •) is a geodesic of M. Then for all a, c(a, •, •) verifies the geodesic equations ∇ s x s (a, s) + r(a, s, 0) = 0, ∀s (2.16a)

∇ 2 s q(a, s, t) + |q(a, s, t)| r(a, s, t) + r(a, s, t) T = 0, ∀t, s, (2.16b)
where (x = c| t=0 , q = c t / |c t |) is the SRV representation of c and r is given by r(a, s, t) = 1 t R(q, ∇ s q)c s (a, s, τ ) τ,t dτ.

Recall that we use the notation w T = w, v(a, s, t) v(a, s, t) with v = c t /|c t | for the tangential component of a tangent vector w ∈ T c(a,s,t) M . To characterize the way these geodesics spread out, we consider the Jacobi field CHAPTER 2. RIEMANNIAN FRAMEWORK

J : [0, 1] → T M, J(s, •) = ∂ ∂a a=0 c(a, s, •).
By decomposing ∇ 2 s J(s, 0) and ∇ t ∇ 2 s J(s, τ ) we can write the second order variation of J as

∇ 2 s J(s, t) = ∇ a ∇ s x s + R(x s , J)x s a=0,t=0 0,t + t 0 ∇ 2 s ∇ t J + R(c t , c s )∇ s J + ∇ s (R(c t , c s )J) a=0,t=τ τ,t
dτ.

(2.17) The term ∇ 2 s ∇ t J can be expressed as a function of ∇ 2 s ∇ a q by twice differentiating the equation ∇ a q = |c t | -1/2 (∇ a c t -1 2 ∇ a c t T ) according to s. This gives

∇ 2 s ∇ a q = ∇ 2 s (|c t | -1 2 )(∇ a c t -1 2 ∇ a c t T ) + 2∇ s (|c t | -1 2 ) ∇ s ∇ a c t -1 2 ∇ s ∇ a c t T -1 2 ∇ a c t , ∇ s v v -1 2 ∇ a c t , v ∇ s v + |c t | -1 2 ∇ 2 s ∇ a c t -1 2 ∇ 2 s ∇ a c t T -∇ s ∇ a c t , ∇ s v v -∇ s ∇ a c t , v ∇ s v -∇ a c t , ∇ s v ∇ s v -1 2 ∇ a c t , ∇ 2 s v -1 2 ∇ a c t , v ∇ 2 s v .
Since ∇ a c t = ∇ t c a = ∇ a J for a = 0, we know that the term we are looking for is ∇ 2 s ∇ t J = ∇ 2 s ∇ a c t . Since ∇ 2 s ∇ a q = ∇ a ∇ 2 s q + R(c s , J)∇ s q + ∇ s (R(c s , J)q), and noticing that W = Z-1 2 Z T is equivalent to Z = W +W T we get

∇ 2 s ∇ t J = W + W T , (2.18)
W = ∇ s ∇ t J, ∇ s v v + ∇ s ∇ t J, v ∇ s v + ∇ t J, ∇ s v ∇ s v + 1 2 ∇ t J, ∇ 2 s v v + 1 2 ∇ t J, v ∇ 2 s v + |c t | 1 2 ∇ a ∇ 2 s q -2∇ s |c t | -1 2 ∇ s ∇ t J -1 2 ∇ s ∇ t J T -1 2 ∇ t J, ∇ s v v -1 2 ∇ t J, v ∇ s v -∇ 2 s |c t | -1 2 ∇ t J -1 2 ∇ t J T + R(c s , J)∇ s q + ∇ s (R(c s , J)q) . (2.19)
The terms ∇ a ∇ s x s (0, s) and ∇ 2 a ∇ s q(0, s, τ ) for all τ ∈ [0, 1] can be obtained by differentiating the geodesic equations 2.16a and 2.16b

∇ a ∇ s x s + ∇ a r = 0, t = 0, s ∈ [0, 1], a = 0, ∇ a ∇ 2 s q + ∇ a |q| r + r T + |q| ∇ a r + ∇ a r T = 0, t, s ∈ [0, 1], a = 0.
The first one gives

∇ a ∇ s x s (0, s) = -∇ a r(0, s, 0), (2.20)
and for all s and t we get

∇ a ∇ 2 s q = -|c t | 1 2 ∇ a r + ∇ a r T -|c t | -1 2 r, ∇ t J v + r, v ∇ t J + 1 2 ∇ t J, v r -3r T .
(2.21)

The only term left to compute is the variation ∇ a r, which is by definition

∇ a r(0, s, t) = 1 t ∇ a V t (0, s, τ ) dτ,
if we define V t for any fixed t by

V t (a, s, τ ) = [R(q, ∇ s q)c s (a, s, τ )] τ,t , τ ∈ [t, 1].
Since the covariant derivative of V t in τ vanishes, we can write for any t ≤ τ ≤ 1

∇ a V t (0, s, τ ) = ∇ a V t (0, s, t) + τ t R(c t , J)| u=t (V t (0, s, u)) du.
Integrating this equation according to τ from t to 1 we obtain

∇ a r(0, s, t) = (1 -t)∇ a V t (0, s, t) + R(c t , J)| t 1 t (1 -τ )V t (0, s, τ ) dτ , (2.22
) where, since V t (0, s, t) = R(q, ∇ s q)c s (0, s, t), we get for τ = t ∇ a V t | t = ∇ J R(q, ∇ s q)c s +R(∇ a q, ∇ s q)c s +R(q, ∇ a ∇ s q)c s +R(q, ∇ s q)∇ s J,

(2.23) with finally

∇ a q = |c t | -1 2 ∇ t J -1 2 ∇ t J T , (2.24)
∇ a ∇ s q = |c t | -1 2 ∇ s ∇ t J -1 2 ∇ s ∇ t J T -1 2 ∇ t J, ∇ s v v (2.25) -1 2 ∇ t J, v ∇ s v + ∇ s |c t | -1 2 ∇ t J -1 2 ∇ t J T + R(J, c s )q.
We can notice that, however complicated, the numbered equations 2.17 to 2.25 when put together define a partial differential equation verified by the Figure 2.4: Optimal deformations between pairs of geodesics (in black) of the upper half-plane H, for our metric (in blue) and for the L 2 -metric (in green). The orientation of the right-hand curve is inverted in the second image compared to the first, and in the fourth compared to the third. Jacobi field J. They allow us to iteratively compute J(s + ε, •) and ∇ s J(s + ε, •), for a fixed step ε > 0, knowing J(s, •) and ∇ s J(s, •). Indeed, we can estimate ∇ t J(s, •) since J(s, t) is known for all t, as well as ∇ t ∇ s J(s, •) since ∇ s J(s, t) is known for all t, and finally ∇ s ∇ t J = ∇ t ∇ s J + R(c s , c t )J. Assuming that we are able to compute the covariant derivative ∇ J R of the curvature tensor, for example if we are in a symmetric space (then it is zero), we obtain an algorithm to compute the Jacobi fields in the space of curves. To summarize : Algorithm 2.3 (Jacobi fields in the space of curves in a symmetric space).

Input : c ∈ M, j 0 , w 0 ∈ T c 0 M. Initialization : Set J(0, t) = j 0 (t) and ∇ s J(0, t) = w 0 (t) for t ∈ [0, 1]. Heredity : For k = 0, . . . , m -1, set s = kε with ε = 1/m and 1. for all t, set

∇ t J(s, t) = lim δ→0 1 δ J(s, t + δ) t+δ,t -J(s, t) , ∇ t ∇ s J(s, t) = lim δ→0 1 δ ∇ s J(s, t + δ) t+δ,t -∇ s J(s, t) , ∇ s ∇ t J(s, t) = ∇ t ∇ s J(s, t) + R(c s , c t )J(s, t);
2. compute r(s, t) = 1 t R(q, ∇ s q)c s (s, τ ) τ,t dτ for all t ∈ [0, 1];

3. compute ∇ a q(0, s, t), ∇ s ∇ a q(0, s, t) and ∇ a V t (0, s, t) for all t using 2.24, 2.25 and 2.23, and deduce ∇ a r(0, s, t) using Equation 2.22;

4. compute ∇ a ∇ s x s (0, s) and ∇ a ∇ 2 s q(0, s, t) for all t using 2.20 and 2.21;

2.4. OPTIMAL MATCHING BETWEEN CURVES 47 5. compute W (s, t) using 2.19 and ∇ 2 s ∇ t J(s, t) = W (s, t) + W T (s, t), and deduce ∇ 2 s J(s, t) for all t using Equation 2.17;

6. finally, for all t ∈ [0, 1], set

J(s + ε, t) = [ J(s, t) + ε∇ s J(s, t) ] s,s+ε , ∇ s J(s + ε, t) = ∇ s J(s, t) + ε∇ 2 s J(s, t) s,s+ε .
Output : J(1).

The different steps of an iteration of the geodesic shooting algorithm, executed using Algorithms 2.1 and 2.3, are shown in Figure 2.3 for two curves of the hyperbolic half-plane H 2 , and examples of geodesic paths of curves in the same space are shown in Figure 2.4. The blue deformations are geodesics for metric G while the green ones are L 2 -geodesics, i.e. each pair of points are linked by a geodesic of the base manifold H 2 . We can already notice that our metric has a tendency to "shrink" the curves in the center of the deformation compared to the L 2 -metric. We can also see the influence of the orientation of the curves, on which the deformations depend. The necessary discretization effort to obtain these simulations is detailed in Chapter 3, and further simulations are given in Section 3.3, where we consider examples in zero and positive curvature as well.

Optimal matching between curves

Recall that our motivation to choose a reparameterization invariant metric was to induce a Riemannian structure on the quotient space S of curves modulo reparameterization, classically interpreted as the shape space. If the scalar product G c is the same at all points c ∈ M that project on the same "shape" -as it is the case for our metric -then it induces a Riemannian structure on the quotient space, which allows us to compare curves regardless of their parameterization. For closed curves, this amounts to considering only the outline of an object; for an open curve representing the evolution in time of a given process, this enables us to analyze it regardless of speed or pace.

Since the geodesics of the quotient space are the projected horizontal geodesics of the total space, solving the boundary value problem in the shape space can be achieved either through the construction of horizontal geodesics e.g. by minimizing the horizontal path energy, or by incorporating the optimal reparameterization of one of the boundary curves as a parameter in the CHAPTER 2. RIEMANNIAN FRAMEWORK optimization problem. The former has been applied to metrics where horizontality is easy to characterize, such as almost local metrics [START_REF] Bauer | Almost local metrics on shape space of hypersurfaces in n-space[END_REF], for which horizontality equals normality, while the latter has been used for Sobolev metrics [START_REF] Bauer | A numerical framework for sobolev metrics on the space of curves[END_REF] -for which finding the horizontal energy implies the much harder task of inverting a differential operator -as well as the elastic metric of the square root velocity framework [START_REF] Srivastava | Shape analysis of elastic curves in euclidean spaces[END_REF], [START_REF] Zhang | Video-based action recognition using rate-invariant analysis of covariance trajectories[END_REF]. For this particular metric, a quotient structure is carefully developed in [START_REF] Lahiri | Precise matching of pl curves in mathbbr n in the square root velocity framework[END_REF], where the authors prove that if at least one of two curves is piecewise-linear, then there exists a minimizing geodesic between the two, and give a precise algorithm to solve the matching problem ; in [START_REF] Bruveris | Optimal reparametrizations in the square root velocity framework[END_REF], it is proven that there always exists a minimizing geodesic between two C 1 plane curves. For elastic metrics in general, the authors in [START_REF] Barbara | Quotient elastic metrics on the manifold of arc-length parameterized plane loops[END_REF] restrict to arc-length parameterized curves and exploit the resulting simplifications to characterize the quotient structure.

Here we characterize the quotient structure for an elastic metric G a,b , and in particular for our metric, without restricting to arc-length parameterized curves. We introduce a simple algorithm allowing to compute geodesics of the quotient space using a canonical decomposition of a path in the associated principal bundle, thereby yielding optimal matchings between curves.

The quotient structure

We consider the quotient S = M/Diff + ([0, 1], M ) of the space of curves by the diffeomorphism group. This quotient is not a manifold, as it has singularities, i.e. points with non trivial isotropy group. If we get rid of these singularities and restrict ourselves to elements of M on which the diffeomorphism group acts freely, then the space of free immersions M f , the quotient shape space S f = M f /Diff + ([0, 1]) and the group of diffeomorphisms Diff + ([0, 1]) form a principal bundle, the fibers of which are the sets of all the curves that are identical modulo reparameterization, i.e. that project on the same "shape". We denote by π : M f → S f the projection of the fiber bundle and by c := π(c) ∈ S f the shape of a curve c ∈ M f . Any tangent vector w to M in c can then be decomposed as the sum

T c M f w = w ver + w hor ∈ Ver c ⊕ Hor c
of a vertical part w ver belonging to the vertical subspace, consisting of all vectors tangent to the fibers of M f over S f , i.e. those which have an action of reparameterizing the curve without changing its shape

Ver c = ker T c π = mv = mc /|c | : m ∈ C ∞ ([0, 1], R), m(0) = m(1) = 0 ,
and a horizontal part w hor belonging to the G-orthogonal of Ver c

Hor c = (Ver c ) ⊥ G = {h ∈ T c M f : G c (h, mv) = 0, ∀m ∈ C ∞ ([0, 1], R), m(0) = m(1) = 0} .
If G is constant along the fibers, i.e. verifies property (2.1), then there exists a Riemannian metric G on the shape space S f such that π is a Riemannian submersion from

(M f , G) to (S f , G), G c (w hor , z hor ) = G π(c) (T c π(w), T c π(z)) ,
where w hor and z hor are the horizontal parts of w and z, as well as the horizontal lifts of T c π(w) and T c π(z), respectively. This expression defines G in the sense that it does not depend on the choice of the representatives c, w and z ([36], §29.21). If a geodesic for G has a horizontal initial speed, then its speed vector stays horizontal at all times -we say it is a horizontal geodesic -and projects on a geodesic of the shape space for G ([36], §26.12). To compute the distance between two shapes c 0 and c 1 in the quotient space we choose a representative c 0 of c 0 and compute the distance (in M f ) to the closest representative of c 1

d (c 0 , c 1 ) = inf d (c 0 , c 1 • ϕ) | ϕ ∈ Diff + ([0, 1]) .
By definition, the distance in the quotient space allows us to compare curves regardless of parameterization

d c 0 • φ, c 1 • ψ = d (c 0 , c 1 ) , ∀φ, ψ ∈ Diff + ([0, 1]).
We now characterize the horizontal subspace for the general elastic metric G a,b , of which our metric G = G 1, 1 2 is a special case, which can be defined for manifold-valued curves by

G a,b c (w, z) = w(0), z(0) + 1 0 a 2 ∇ w N , ∇ z N + b 2 ∇ w T , ∇ z T d .
The following result gives the decomposition of a tangent vector associated to the quotient structure for this metric.

(a/b) 2 -1 ∇ t h, ∇ t v -∇ 2 t h, v + |c | -1 ∇ t c , v ∇ t h, v = 0.

CHAPTER 2. RIEMANNIAN FRAMEWORK

In particular, for a = 2b = 1 the horizontal subspace is given by

Hor c = {h ∈ T c M : ∀t ∈ [0, 1], 3 ∇ t h, ∇ t v -∇ 2 t h, v + |c | -1 ∇ t c , v ∇ t h, v = 0}.
The vertical and horizontal parts of a tangent vector w ∈ T c M are given by w ver = mv, w hor = wmv, where the real function m ∈ C ∞ ([0, 1], R) verifies m(0) = m(1) = 0 and the ordinary differential equation

m -∇ t c /|c |, v m -4|∇ t v| 2 m = ∇ 2 t w, v -3 ∇ t w, ∇ t v -∇ t c /|c |, v ∇ t w, v .
(2.26)

Proof. Let h ∈ T c M be a tangent vector. It is horizontal if and only if it is orthogonal to any vertical vector, that is any vector of the form mv with m ∈ C ∞ ([0, 1], R) such that m(0) = m(1) = 0. We have ∇ t (mv) = m v + m∇ t v and since ∇ t v, v = 0 we get ∇ t (mv) N = m∇ t v and ∇ t (mv) T = m v. The scalar product can then be written

G a,b c (h, mv) = h(0), m (0)v(0)
+ 1 0 (a 2 ∇ t h N , ∇ t (mv) N + b 2 ∇ t h T , ∇ t (mv) T ) dt |c | = 1 0 (a 2 m ∇ t h, ∇ t v + b 2 m ∇ t h, v ) dt |c | = 1 0 a 2 m ∇ t h, ∇ t v dt |c | - 1 0 b 2 m d dt ∇ t h, v |c | -1 dt = 1 0 m/|c | (a 2 -b 2 ) ∇ t h, ∇ t v -b 2 ∇ 2 t h, v + b 2 ∇ t c , v ∇ t h, v |c | -1 dt,
by integration by parts. The vector h is horizontal if and only if G c (h, mv) = 0 for all such m, and so multiplying by |c |/b 2 gives the desired equation. Now consider a tangent vector w and a real function m : [0, 1] → R such that m(0) = m(1) = 0. Then according to the above, wmv is horizontal if and only if it verifies

3 ∇ t (w -mv), ∇ t v -∇ 2 t (w -mv), v + |c | -1 ∇ t c , v ∇ t (w -mv), v = 0, i.e., since ∇ t v, v = 0, ∇ 2 t v, v = -|∇ t v| 2 and ∇ 2 t (mv) = m v + 2m ∇ t v + m∇ 2 t v, if 3 ∇ t w, ∇ t v -3|∇ t v| 2 m -∇ 2 t w, v + m -m|∇ t v| 2 + |c | -1 ∇ t c , v ( ∇ t w, v -m ) = 0,

OPTIMAL MATCHING BETWEEN CURVES

51

which is what we wanted.

Computing geodesics in the shape space

Recall that the geodesic path s → c(s) between the shapes of two curves c 0 and c 1 is the projection of the horizontal geodesic -if it existss → c h (s) linking c 0 to the fiber of c 1 in M, i.e. such that c h (0

) = c 0 , c h (1) ∈ π -1 (c 1 ) and ∂ s c h (s) ∈ Hor c h (s) for all s ∈ [0, 1], c = π(c h ).
The end point of c h then gives the optimal reparameterization c 1 • ϕ of the target curve c 1 with respect to the initial curve c 0 , i.e. such that

d(c 0 , c 1 ) = d(c 0 , c 1 • ϕ).
Here we propose a method to approach the horizontal geodesic c h . To that end we decompose any path of curves s → c(s) in M into a horizontal path composed with a path of reparameterizations, c(s

) = c hor (s) • ϕ(s), or equivalently c(s, t) = c hor (s, ϕ(s, t)) ∀s, t ∈ [0, 1], (2.27)
where the path [0, 1] s → c hor (s) ∈ M is such that c hor s (s) ∈ Hor c hor (s) for all s ∈ [0, 1], and [0, 1] s → ϕ(s) ∈ Diff + ([0, 1]) is a path of increasing diffeomorphisms. The horizontal and vertical parts of the speed vector of c can be expressed in terms of this decomposition. Indeed, by taking the derivative of (3.5) with respect to s and t we obtain

c s (s) = c hor s (s) • ϕ(s) + ϕ s (s) • c hor t (s) • ϕ(s), (2.28a)
c t (s) = ϕ t (s) • c hor t (s) • ϕ(s), (2.28b)
and so with v hor (s, t) := c hor t (s, t)/|c hor t (s, t)|, since ϕ t > 0, (2.28b) gives v(s) = v hor (s) • ϕ(s).

We can see that the first term on the right-hand side of Equation (2.28a) is horizontal. Indeed, for any path of real functions m :

[0, 1] → C ∞ ([0, 1], R), s → m(s, •) such that m(s, 0) = m(s, 1) = 0 for all s, since G is reparame- terization invariant we have G c hor s (s) • ϕ(s), m(s) • v(s) = G c hor s (s) • ϕ(s), m(s) • v hor (s) • ϕ(s) = G c hor s (s), m(s) • ϕ(s) -1 • v hor (s) = G c hor s (s), m(s) • v hor (s) with m(s) = m(s) • ϕ(s) -1
. Since m(s, 0) = m(s, 1) = 0 for all s, the vector m(s) • v hor (s) is vertical and its scalar product with the horizontal vector c hor s (s) vanishes. On the other hand, the second term on the right hand-side of Equation (2.28a) is vertical, since it can be written

ϕ s (s) • c hor t • ϕ(s) = m(s) • v(s), with m(s) = |c t (s)|ϕ s (s)/ϕ t (s) verifying m(s, 0) = m(s, 1) = 0 for all s.
Finally, the vertical and horizontal parts of the speed vector c s (s) are given by

c s (s) ver = m(s) • v(s) = |c t (s)|ϕ s (s)/ϕ t (s) • v(s),
(2.29a)

c s (s) hor = c s (s) -m(s) • v(s) = c hor s (s) • ϕ(s).
(2.29b) Definition 2.1. We call c hor the horizontal part of the path c with respect to the metric G.

Proposition 2.6. The horizontal part of a path of curves c is at most the same length as c L G (c hor ) ≤ L G (c).

Proof. Since the metric G is reparameterization invariant, the squared norm of the speed vector of the path c at time s ∈ [0, 1] is given by

c s (s, •) 2 G = c hor s (s, ϕ(s, •)) 2 G + |ϕ s (s, •)| 2 c hor t (s, ϕ(s, •) 2 G = c hor s (s, •) 2 G + |ϕ s (s, •)| 2 c hor t (s, •) 2 G , where • 2 G := G(•, •).
This gives c hor s (s) G ≤ c s (s) for all s and so L G (c hor ) ≤ L G (c). Now we will see how the horizontal part of a path of curves can be computed.

Proposition 2.7 (Horizontal part of a path). Let s → c(s) be a path in M. Then its horizontal part is given by c hor (s, t) = c(s, ϕ(s) -1 (t)), where the path of diffeomorphisms s → ϕ(s) is solution of the partial differential equation

ϕ s (s, t) = m(s, t)/|c t (s, t)| • ϕ t (s, t), (2.30)
with initial condition ϕ(0, •) = Id, and where m(s) : [0, 1] → R, t → m(s, t) verifies m(s, 0) = m(s, 1) = 0 and is solution for all s of the ordinary differential equation

m tt -∇ t c t /|c t |, v m t -4|∇ t v| 2 m = ∇ 2 t c s , v -3 ∇ t c s , ∇ t v -∇ t c t /|c t |, v ∇ t c s , v .
Proof. We have seen in Equation (2.29a) that the vertical part of c s (s) can be written as m(s) • v(s) where m(s) = |c t (s)|ϕ s (s)/ϕ t (s), and as the norm of the vertical part of c s (s), m(s) is solution of the ODE (2.26) for all s.

If we take the horizontal part of the geodesic linking two curves c 0 and c 1 , we will obtain a horizontal path linking c 0 to the fiber of c 1 which will no longer be a geodesic path. However this path reduces the distance between c 0 and the fiber of c 1 , and gives a "better" representative ĉ1 = c 1 • ϕ(1) of the target curve. By computing the geodesic between c 0 and this new representative ĉ1 , we are guaranteed to reduce once more the distance to the fiber. The algorithm that we propose simply iterates these two steps. Output : c hor , ĉ1 .

This algorithm yields an approximation of the horizontal geodesic between c 0 and the fiber of c 1 and of the optimal reparameterization ĉ1 of c 1 . Before we test Algorithms 2.2 and 2.4, we first introduce a formal discretization of the continuous model presented so far. ∇ s q w (0, t), ∇ s q z (0, t) dt.

(3.1)

This definition does not depend on the choice of c w and c z and we can reformulate this scalar product in terms of (covariant) derivatives of w and z. Indeed, the same computations as those developed in the proof of Proposition 2.1 give

∇ s q w (0, t) = |c | -1 (∇ t w N + 1 2 ∇ t w T
), and so we recover the original definition

G c (w, z) = w(0), z(0) + 1 0 ∇ w N , ∇ z N + 1 4 ∇ w T , ∇ z T d . (3.2)
With this formulation it is clear that G is reparameterization invariant, and that it belongs to the class of elastic metrics. Two curves c 0 , c 1 ∈ M can be compared using the geodesic distance induced by G, i.e. by computing the length of the shortest path of curves

[0, 1] s → c(s) ∈ M from c 0 to c 1 d G (c 0 , c 1 ) = inf {L(c) : c(0) = c 0 , c(1) = c 1 } . (3.3)
The length of a path c can be written in terms of its SRV representation as

L(c) = 1 0 |x s (s)| 2 + 1 0
|∇ s q(s, t)| 2 dt ds, and we have shown that the geodesic paths of M, i.e. those which achieve the infimum in (3.3), are the paths [0, 1] s → c(s) ∈ M for which the SRV representative s → (x(s), q(s)) verifies the equations

∇ s x s (s) + r(s, 0) =0, ∇ 2 s q(s, t) + |q(s, t)| r(s, t) + r(s, t) T =0, (3.4)
for all s, t ∈ [0, 1], where the vector field r is given by

r(s, t) = 1 t R(q, ∇ s q)c s (s, τ ) τ,t dτ, t ∈ [0, 1].
We exploited these geodesic equations to solve the initial value problem -in the form of the exponential map -and the boundary value problem using geodesic shooting, which additionally required the description of Jacobi fields.

SUMMARY OF THE CONTINUOUS MODEL 57

The space of unparameterized curves

Restricting to free immersions M f , we characterized the decomposition of the tangent bundle T M f associated to the principal bundle π :

M f → S f = M f /Diff + ([0, 1], M ) on the shape space T M f = Ver ⊕ Hor.
The vertical subspace is given by

Ver c = ker T c π = mv = mc /|c | : m ∈ C ∞ ([0, 1], R), m(0) = m(1) = 0 ,
and we showed that the G-orthogonal horizontal subspace Hor c = (Ver c ) ⊥ G for an elastic metric parameterized by a, b ∈ R + is

Hor c = {h ∈ T c M : ∀t ∈ [0, 1], (a/b) 2 -1 ∇ t h, ∇ t v -∇ 2 t h, v + |c | -1 ∇ t c , v ∇ t h, v = 0}
, and in particular for G,

Hor c = {h ∈ T c M : ∀t ∈ [0, 1], 3 ∇ t h, ∇ t v -∇ 2 t h, v + |c | -1 ∇ t c , v ∇ t h, v = 0}.
We introduced a simple algorithm that, given two curves c 0 and c 1 , computes the horizontal geodesic s → c h (s) -if it exists -linking c 0 to the fiber of c 1 in M. It is the horizontal lift at c 0 of the geodesic between the shapes of c 0 and c 1 , and its end point gives the optimal reparameterization c 1 • ϕ of c 1 with respect to c 0 . We decomposed any path of curves s → c(s) in M into a horizontal path c hor -that we refer to as the horizontal part of c -composed with a path of reparameterizations s → ϕ(s)

∈ Diff + ([0, 1]) c(s, t) = c hor (s, ϕ(s, t)) ∀s, t ∈ [0, 1], (3.5)
and showed that the path of diffeomorphisms s → ϕ(s) is solution of the partial differential equation

ϕ s (s, t) = m(s, t)/|c t (s, t)| • ϕ t (s, t), (3.6)
where m(s) : [0, 1] → R, t → m(s, t) is the vertical component of the velocity field c s (s), and as such verifies for all s the ODE

m tt -∇ t c t /|c t |, v m t -4|∇ t v| 2 m = ∇ 2 t c s , v -3 ∇ t c s , ∇ t v -∇ t c t /|c t |, v ∇ t c s , v .
Iteratively computing the geodesic c between c 0 and the current reparameterization c 1 • ϕ of c 1 , and its horizontal part c hor , we reduce the distance from c 0 to the fiber of c 1 at each step.

The discrete model

Applications usually give access to a finite number of observations of a continuous process and provide series of points instead of continuous curves. It is therefore important to discretize the model presented in Chapter 2. Here we develop a discrete model in the case where the base manifold M has constant sectional curvature K. This model could more generally be extended to symmetric spaces.

The Riemannian structure

We consider the product manifold M n+1 to model the space of "discrete curves" given by n + 1 points, for a fixed n ∈ N * . Its tangent space at a given point α = (x 0 , . . . , x n ) is given by

T α M n+1 = {w = (w 0 , . . . , w n ) : w k ∈ T x k M, ∀k}.
Assuming that there exists a connecting geodesic between x k and x k+1 for all k -which seems reasonable considering that the points x k should be "close" since they correspond to the discretization of a continuous curve -we use the following notations

τ k = log M x k x k+1 , q k = √ n τ k / |τ k |, v k = τ k /|τ k |, (3.7)
as well as w k T = w k , v k v k and w N k = w k -w k T to refer to the tangential and normal components of a tangent vector w k ∈ T x k M . Given a tangent vector w ∈ T α M n+1 , we consider a path of piecewise geodesic curves [0, 1] 2 (s, t) → c w (s, t) ∈ M such that c w (0, k n ) = x k for k = 0, . . . , n, c w (s, •) is a geodesic of M on the interval [ k n , k+1 n ] for all s ∈ [0, 1] and k -and in particular c w t (0, k n ) = nτ k -and such that c w s (0, k n ) = w k . Then we define the scalar product between w and z by

G n α (w, z) = w 0 , z 0 + 1 n n-1 k=0 ∇ s q w (0, k n ), ∇ s q z (0, k n ) . (3.8)
This definition is a discrete analog of (3.1), and just as in the continuous case, it does not depend on the choices of c w and c z . Indeed, we can also obtain a discrete analog of (3.2).

Proposition 3.1. The scalar product between two tangent vectors w, z ∈ T α M n+1 can also be written where the map D τ : T α M n+1 → T α M n+1 , w → D τ w = (D τ w) 0 , . . . , (D τ w) n is defined by

G n α (w, z) = w 0 , z 0 + n-1 k=0 (D τ w) k N , (D τ z) k N + 1 4 (D τ w) k T , (D τ z) k T 1 |τ k | ,
(D τ w) k := 1 n ∇ t c w s (0, k n ) = (w k+1 -w k ) T + b -1 k (w k+1 -a k w k ) N ,
and the coefficients a k and b k take the following values depending on the sectional curvature K of the base manifold Before we prove this proposition, let us recall a well-known result about Jacobi fields that will prove useful to derive the equations in the discrete case.

M    a k = cosh |τ k |, b k = sinh |τ k |/|τ k |, if K = -1, a k = 1, b k = 1, if K = 0, a k = cos |τ k |, b k = sin |τ k |/|τ k |, if K = +1.
Lemma 3.1. Let γ : [0, 1] → M be a geodesic of a manifold M of constant sectional curvature K, and J a Jacobi field along γ. Then the parallel transport of J(t) along γ from γ(t) to γ(0) is given by

J(t) t,0 = J T (0) + ãk (t)J N (0) + t ∇ t J T (0) + bk (t)∇ t J N (0), for all t ∈ [0, 1], where    ãk (t) = cosh (|γ (0)|t) , bk (t) = sinh(|γ (0)|t)/|γ (0)|, K = -1, ãk (t) = 1, bk (t) = t, K = 0, ãk (t) = cos (|γ (0)|t) , bk (t) = sin(|γ (0)|t)/|γ (0)|, K = +1.
Proof of Lemma 1. The proof is reminded in Appendix A.

Proof of Proposition 5. Let α ∈ M n+1 be a "discrete curve" and w, z ∈ T α M n+1 tangent vectors at α. Consider a path of piecewise geodesic curves s → c w (s) that verifies all the conditions given above to define G n (w, z), and set

(D τ w) k := 1 n ∇ t c w s (0, k n ).
Then by definition, the vector field

J k (u) = c w s (s, k+u n ), u ∈ [0, 1] is a Jacobi field along the geodesic linking x k to x k+1 , verifying J k (0) = w k , J k (1) = w k+1 and ∇ u J k (0) = (D τ w) k . Applying Lemma 3.1 gives w k+1 = w k T + a k w k N + (D τ w) k T + b k (D τ w) k N .
Taking the tangential part and then the normal parts on both sides gives

(w k+1 ) T = w k T + (D τ w) k T , (w k+1 ) N = a k w k N + b k (D τ w) k N .
and so

(D τ w) k = (D τ w) k T + (D τ w) k N = (w k+1 -w k ) T + b -1 k (w k+1 - a k w k ) N .
Finally, we observe that the covariant derivative involved in the definition of G n can be written

∇ s q w (0, k n ) = |c w t (0, k n )| -1 2 (∇ s c w t (0, k n ) -1 2 ∇ s c w t (0, k n ) T ) = |nτ k | -1 2 n(D τ w) k -1 2 n(D τ w) k T ),
i.e.

∇ s q w (0, k n ) = (n/|τ k |) 1/2 (D τ w) k N + 1 2 (D τ w) k T .
Injecting this into (3.8) gives the desired formula for the scalar product. Now we present our main result, that is, the convergence of the discrete model toward the continuous model. Definition 3.1. Let α = (x 0 , . . . , x n ) ∈ M n+1 be a discrete curve, and t → c(t) ∈ M a smooth curve. We say that α is the discretization of size n of c when c( k n ) = x k for all k = 0, . . . , n. A path s → α(s) = (x 0 (s), . . . , x n (s)) ∈ M n+1 of discrete curves is the discretization of a path of smooth curves s → c(s) ∈ M when α(s) is the discretization of c(s) for all s ∈ [0, 1], i.e. when x k (s) = c(s, k n ) for all s and k. We will still use this term if c is not smooth, and speak of the only path of piecewise-geodesic curves of which α is the discretization.

THE DISCRETE MODEL

61 Let [0, 1] s → α(s) = (x 0 (s), . . . , x n (s)) ∈ M n+1 be a path of discrete curves. Defining τ k (s) and q k (s) as in (3.7) for all s ∈ [0, 1], the path α can be represented by its SRV representation [0, 1] → M × T α M n+1 , s → x 0 (s), (q k (s)) 0≤k≤n-1 .

(3.10)

To compute the squared norm of its speed vector α (s), consider the path of piecewise geodesic curves

[0, 1] 2 (s, t) → c(s, t) ∈ M such that c(s, k n ) = x k (s) and c t (s, k n ) = nτ k (s)
for all s and k. Then, notice that we have

∇ s q(s, k n ) = ∇ s q k (s), (D τ α (s)) k = 1 n ∇ t c s (s, k n ) = ∇ s τ k (s), (3.11)
and so the squared norm of the speed vector of α can be expressed in terms of the SRV representation

G n (α (s), α (s)) = |x 0 (s)| 2 + 1 n n-1 k=0 |∇ s q k (s)| 2 .
In the following result, we show that if s → α(s) is the discretization of a path s → c(s) ∈ M of continuous curves, then its energy with respect to

G n , E n (α) = 1 2 1 0 |x 0 (s)| 2 + 1 n n-1 k=0 |∇ s q k (s)| 2 ds, (3.12)
gets closer to the energy (2.11) of c with respect to G as the size of the discretization grows.

Theorem 3.1 (Convergence of the discrete model to the continuous model).

Let s → c(s) be a C 1 -path of C 2 -curves with non vanishing derivative with respect to t. This path can be identified with an element (s, t) → c(s, t) of

C 1,2 ([0, 1] × [0, 1], M ) such that c t = 0. Consider the C 1 -path in M n+1 , s → α(s) = (x 0 (s), . . . , x n (s))
, that is the discretization of size n of c. Then there exists a constant λ > 0 that does not depend on c and such that for n large enough,

|E(c) -E n (α)| ≤ λ n (inf |c t |) -1 |c s | 2 2,∞ (1 + |c t | 1,∞ ) 3 ,
where E and E n are the energies with respect to metrics G and G n respectively and where |w(s, t)| denotes the supremum over both s and t of a vector field w along c. Remark 3.2. Note that since we assume that c is a C 1 -path of C 2 -curves, the following norms are bounded for i = 1, 2,

|c t | 1,∞ := |c t | ∞ + |∇ t c t | ∞ , |c s | 2,∞ := |c s | ∞ + |∇ t c s | ∞ + |∇ 2 t c s | ∞ , 62
|c t | ∞ , |c s | ∞ , |∇ i t c t | ∞ , |∇ i t c s | ∞ < ∞. Proof of Theorem 3.1.
The proof is put off to Section 3.4. Now that we have established a formal Riemannian setting to study discrete curves defined by a series of points, and that we have studied its link to the continuous model, we need to derive the equations of the corresponding geodesics and Jacobi fields to apply the methods described in Chapter 2. To simplify the formulas, we first introduce some notations.

Computing geodesics in the discrete setting

Notations The purpose of the notations that we introduce here is to lighten the equations derived in the rest of the chapter. For any discrete curve α = (x 0 , . . . , x n ) ∈ M n+1 we define for all 0 ≤ k ≤ n, using the coefficients a k and b k defined by (3.9) and (3.7), the functions f k , g k :

T x k M → T x k M , f k : w → w T + a k w N , g k : w → |q k |(2w T + b k w N ).
and for k = 0, . . . , n -1, the functions f

(-) k , g (-) k : T x k+1 M → T x k M by f (-) k = f k • P x k+1 ,x k γ k , g (-) k = g k • P x k+1 ,x k γ k
, where γ k denotes the geodesic between x k and x k+1 , which we previously assumed existed. Notice that when the discretization gets "thinner", that is n → ∞, |τ k | → 0 while n|τ k | stays bounded for all 0 ≤ k ≤ n, we get in the non flat setting, for any fixed

w ∈ T x k+1 M , f k (w) = w + o(1/n) and g k (w) = |q k |(w+w T )+o(1/n) -in
the flat setting, these are always equalities. Now if we consider a path s → α(s) = (x 0 (s), . . . , x n (s)) of discrete curves, we can define for each s the functions

f k (s), g k (s) : T x k (s) M → T x k (s) M, f k (s) (-) , g k (s) (-) : T x k+1 (s) M → T x k (s) M,
for 0 ≤ k ≤ n and 0 ≤ k ≤ n -1 respectively, corresponding to the discrete curve α(s). It is of interest for the rest of this chapter to compute the covariant derivatives of these maps with respect to s. Lemma 3.2. The first and second order covariant derivatives of f k and g k with respect to s are functions

T x k (s) M → T x k (s) M defined by ∇ s f k (w) = ∂ s a k w N + (1 -a k ) w, ∇ s v k v k + w, v k ∇ s v k , ∇ s g k (w) = ∂ s |q k |/|q k |g k (w) + |q k |∂ s b k w N + |q k |(2 -b k ) w, ∇ s v k v k + w, v k ∇ s v k , ∇ 2 s f k (w) = ∂ 2 s a k w N -2∂ s a k w, ∇ s v k v k + w, v k ∇ s v k + (1 -a k ) w, ∇ 2 s v k v k + 2 w, ∇ s v k ∇ s v k + w, v k ∇ 2 s v k , ∇ 2 s g k (w) = ∂ s ∂ s |q k |/|q k | g k (w) + ∂ s |q k |/|q k |∇ s g k (w) + |q k |(2 -b k ) w, ∇ 2 s v k v k + 2 w, ∇ s v k ∇ s v k + w, v k ∇ 2 s v k + ∂ s |q k |(2 -b k ) -2|q k |∂ s b k w, ∇ s v k v k + w, v k ∇ s v k + (∂ s |q k |∂ s b k + |q k |∂ 2 s b k )w N . Proof. For any vector field s → w(s) ∈ T x k (s) M along s → x k (s) we have by definition ∇ s f k (w) = ∇ s f k (w) + f k (∇ s w), ∇ s g k (w) = ∇ s g k (w) + g k (∇ s w), ∇ 2 s f k (w) = ∇ 2 s f k (w) + 2 ∇ s f k (∇ s w) + f k (∇ 2 s w), ∇ 2 s g k (w) = ∇ 2 s g k (w) + 2 ∇ s g k (∇ s w) + g k (∇ 2 s w). Noticing that ∇ s (w T ) = (∇ s w) T + w, ∇ s v k v k + w, v k ∇ s v k and ∇ s (w N ) =
∇ s w -∇ s (w T ), the formulas given in Lemma 3.2 result from simple calculation.

Using these functions, we can deduce the covariant derivatives of f

T x k+1 (s) M → T x k (s) M given by ∇ s f (-) k : w → (∇ s f k ) (-) (w) + f k R (Y k , τ k ) (w k+1 ) , ∇ s g (-) k : w → (∇ s g k ) (-) (w) + g k R (Y k , τ k ) (w k+1 ) , where (∇ s f k )(s) (-) = ∇ s f k (s) • P x k+1 (s),x k (s) γ k (s) , (∇ s g k )(s) (-) = ∇ s g k (s) • P x k+1 (s),x k (s) γ k (s) , Y k = (x k ) T + b k (x k ) N + 1 2 ∇ s τ k T + K 1 -a k |τ k | 2 ∇ s τ k N , (3.13)
and R is the curvature tensor of the base manifold.

Proof. The proof is given in Appendix A.

Geodesic equations and exponential map

With these notations, we can characterize the geodesics for metric G n . The geodesic equations can be derived in a similar way as in the continuous case, that is by searching for the critical points of the energy (3.12). We obtain the following characterization in terms of the SRV representation (3.10).

Proposition 3.2 (Discrete geodesic equations).

A path in M n+1 , s → α(s) = (x 0 (s), . . . , x n (s)), is a geodesic for metric G n if and only if its SRV representation s → x 0 (s), (q k (s)) k verifies the following differential equations

∇ s x 0 + 1 n R 0 + f (-) 0 (R 1 ) + . . . + f (-) 0 • • • • • f (-) n-2 (R n-1 ) = 0, ∇ 2 s q k + 1 n g (-) k R k+1 + f (-) k+1 (R k+2 ) + . . . + f (-) k+1 • • • • • f (-) n-2 (R n-1 ) = 0, (3.14)
for all k = 0, . . . , n -1, with the notations (3.7) and R k := R(q k , ∇ s q k )x k .

Proof. The proof is given in Appendix A.

Remark 3.3. Let [0, 1] s → c(s, •) ∈ M be a C 1 path of smooth curves and [0, 1] s → α(s) ∈ M n+1 the discretization of size n of c. We denote as usual by q := c t /|c t | 1/2 and (q k ) k their respective SRV representations. When n → ∞ and |τ k | → 0 while n|τ k | stays bounded for all 0 ≤ k ≤ n, the coefficients of the discrete geodesic equation (3.14) for α converge to the coefficients of the continuous geodesic equation (3.4) for c, i.e. 1), for all s ∈ [0, 1] and k = 0, . . . , n-1, where r n-1 = 0 and for k = 1, . . . , n-2,

∇ s x 0 (s) = -r 0 (s) + o(1), ∇ 2 s q k (s) = -|q k (s)|(r k (s) + r k (s) T ) + o(
r k (s) := 1 n n-1 =k+1 P l n , k n c R(q, ∇ s q)c s (s, n ) → n→∞ r(s, k n ),
with the exception that the sum starts at = 0 for r 0 . More details on this can be found in Appendix A.

Remark 3.4. Just as in the continuous case, when the base manifold is a Euclidean space, the curvature terms R k 's vanish and we obtain

x 0 (s) = 0, q k (s) = 0, k = 0, . . . , n -1, ∀s ∈ [0, 1],
i.e. the geodesics are those for which the SRV representations are L 2 -geodesics.

Using equations (3.14) we can now build the exponential map, that is, an algorithm allowing us to approximate the geodesic of M n+1 starting from a point (x 0 0 , . . . ,

x 0 n ) ∈ M n+1 at speed (u 0 , . . . , u n ) with u k ∈ T x 0 k
M for all k = 0, . . . , n. In other words, we are looking for a path [0, 1] s → α(s) = (x 0 (s), . . . , x n (s)) such that x k (0) = x 0 k and x k (0) = u k for all k, and that verifies the geodesic equations (3.14). Assume that we know at time s ∈ [0, 1] the values of x k (s) and x k (s) for all k = 0, . . . , n. Then we propagate using

x k (s + ε) = log M x k (s) εx k (s), x k (s + ε) = x k (s) + ε∇ s x k (s) s,s+ε .
In the following proposition, we see how we can compute the acceleration

∇ s x k for each k. Proposition 3.3 (Exponential map in M n+1 ). Let [0, 1]
s → α(s) = (x 0 (s), . . . , x n (s)) be a geodesic path in M n+1 . For all s ∈ [0, 1], the coordinates of its acceleration ∇ s α (s) can be iteratively computed in the following way

∇ s x 0 = - 1 n R 0 + f (-) 0 (R 1 ) + . . . + f (-) 0 • • • • • f (-) n-2 (R n-1 ) , ∇ s x k+1 = ∇ s f k (x k ) + f k (∇ s x k ) + 1 n ∇ s g k (∇ s q k ) + 1 n g k (∇ 2 s q k ) + R(τ k , Y k )(x k+1 ),
for k = 0, . . . , n -1, where the R k 's are defined as in Proposition 3.2, the symbol • denotes the parallel transport from x k+1 (s) back to x k (s) along the geodesic linking them, the maps ∇ s f k and ∇ s g k are given by Lemma 3.2, Y k is given by Equation (3.13) and

∇ s τ k = (D τ α ) k , ∇ s v k = 1 |τ k | ∇ s τ k -∇ s τ k T , ∇ s q k = n |τ k | ∇ s τ k - 1 2 ∇ s τ k T , ∇ 2 s q k = - 1 n g (-) k R k+1 + f (-) k+1 (R k+2 ) + . . . + f (-) k+1 • • • • • f (-) n-2 (R n-1 ) .
Proof. The proof is given in Appendix A.

The equations of Proposition 3.3 allow us to iteratively construct a geodesic s → α(s) in M n+1 for metric G n from the knowledge of its initial conditions α(0) and α (0). The next step is to construct geodesics under boundary constraints, i.e. to find the shortest path between two elements α 0 and α 1 of M n+1 .

Jacobi fields and geodesic shooting As explained in Chapter 2 for the continuous model, we solve the boundary value problem using geodesic shooting. To do so, recall that we need to characterize the Jacobi fields for the metric G n , since these play a role in the correction of the shooting direction at each iteration of the algorithm. Jacobi fields are vector fields that describe the way that geodesics spread out in the Riemannian manifold: for any geodesic s → α(s) in M n+1 and Jacobi field s → J(s) along α, there exists a family of geodesics (-δ, δ) a → α(a, •) such that α(0, s) = α(s) for all s and

J(s) = ∂ ∂a a=0 α(a, s).
Proposition 3.4 (Discrete Jacobi fields). Let s → α(s) = (x 0 (s), . . . , x n (s)) be a geodesic path in M n+1 , s → J(s) = (J 0 (s), . . . , J n (s)) a Jacobi field along α, and (-δ, δ) a → α(a, •) a corresponding family of geodesics, in the sense just described. Then J verifies the second order linear ODE

∇ 2 s J 0 = R(x 0 , J 0 )x 0 - 1 n n-2 k=0 k =0 f (-) 0 • • • • • ∇ a f (-) • • • • • f (-) k (R k+1 ), - 1 n ∇ a R 0 + f (-) 0 (∇ a R 1 ) + . . . + f (-) 0 • • • • • f (-) n-2 (∇ a R n-1 ) ∇ 2 s J k+1 = f k (∇ 2 s J k ) + 2∇ s f k (∇ s J k ) + ∇ 2 s f k (J k ) + 1 n g k (∇ 2 s ∇ a q k ) + 2 n ∇ s g k (∇ s ∇ a q k ) + 1 n ∇ 2 s g k (∇ a q k ) + R(τ k , Y k ) R(Y k , τ k )(J k+1 ) + R(∇ s τ k , Y k )(J k+1 ) + R(τ k , ∇ s Y k )(J k+1 ) + 2R(τ k , Y k )(∇ s J k+1 ), for all 0 ≤ k ≤ n -1, where R k := R(q k , ∇ s q k )
x k and the various covariant derivatives according to a can be expressed as functions of J and ∇ s J,

∇ a R k = R ∇ a q k , ∇ s q k x k + R q k , ∇ s ∇ a q k + R(J, x k )q k x k + R q k , ∇ s q k )∇ s J k , ∇ a τ k = (D τ J) k , ∇ a v k = |τ k | -1 ∇ a τ k -∇ a τ k T , ∇ a q k = n/|τ k | ∇ a τ k -1 2 ∇ a τ k T , ∇ s ∇ a q k = n g k -1 (∇ s J k+1 ) + R(Y k , τ k )(J k+1 ) -∇ s f k (J k ) -f k (∇ s J k ) + n ∇ s g k -1 J k+1 -f k (J k ) , ∇ s ∇ s ∇ a q k = R(∇ s x k , J k )q k + R(x k , ∇ s J k )q k + 2R(x k , J k )∇ s q k - 1 n n-1 =k+1 g (-) k • f (-) k+1 • • • • • f (-) -1 (∇ a R ) - 1 n n-1 =k+1 -1 j=k g (-) k • • • • • ∇ a f (-) j • • • • • f (-) -1 (R ), ∇ s Y k = (∇ s x k ) T + b k (∇ s x k ) N + (1 -b k ) x k , ∇ s v k v k + x k , v k ∇ s v k + ∂ s b k (x k ) N + 1 2 (∇ s ∇ s τ k ) T + K 1 -a k |τ k | 2 (∇ s ∇ s τ k ) N + ∂ s K 1 -a k |τ k | 2 (∇ s τ k ) N + 1 2 -K 1 -a k |τ k | 2 ( ∇ s τ k , ∇ s v k v k + ∇ s τ k , v k ∇ s v k ),
with the notation conventions f

(-) k+1 • . . . • f (-)
k-1 := Id, n-1 =n := 0 and with the maps

∇ a f (-) k (w) = (∇ a f k ) (-) (w) + f k R(Z k , τ k )(w k+1 ) , ∇ a g (-) k (w) = (∇ a g k ) (-) (w) + g k R(Z k , τ k )(w k+1 ) , ∇ s g k -1 (w) = ∂ s |q k | -1 |q k |g k -1 (w) +|q k | -1 ∂ s (b -1 k ) w N + |q k | -1 1/2 -b -1 k w, ∇ s v k v k + w, v k ∇ s v k , and Z k = J k T + b k J k N + 1 2 ∇ a τ k T + K 1 -a k |τ k | 2 ∇ a τ k N .
Proof. The proof is given in Appendix A.

The equations of Proposition 3.4 allow us to iteratively compute the Jacobi field J along a geodesic α -and in particular, its end value J(1) -from the knowledge of the initial conditions {J k (0), 0 ≤ k ≤ n} and {∇ s J k (0), 0 ≤ k ≤ n}. Indeed, if at time s ∈ [0, 1] we have J k (s) and ∇ s J k (s) for all k = 0, . . . , n, then we can propagate using

J k (s + ε) = J k (s) + ε∇ s J k (s) x k ,x k+1 , ∇ s J k (s + ε) = ∇ s J k (s) + ε∇ 2 s J k (s) x k ,x k+1 ,
where ∇ 2 s J k (s) is deduced from ∇ 2 s J k-1 (s) using Proposition 3.4. We can now apply Algorithm 2.2, where we replace the smooth geodesic equations (3.4) by the discrete geodesic equations (3.14) and we solve them using the exponential map described in Proposition 3.3. Notice that in M n+1 , the k th component of the L 2 -logarithm map between two elements α 0 = (x 0 0 , . . . , x 0 n ) and α

1 = (x 1 0 , . . . , x 1 n ) is given by log M x 0 k (x 1 k ). Algorithm 3.1 (Geodesic shooting in M n+1 ). Let (x 0 0 , . . . , x 0 n ), (x 1 0 , . . . , x 1 n ) ∈ M n+1 . Set u = log L 2 α 0 (α 1 )
and repeat until convergence : 1. compute the geodesic s → α(s) starting from α 0 at speed u using Proposition 3.3, 2. evaluate the difference j := log L 2 α(1) (α 1 ) between the target curve α 1 and the extremity α(1) of the obtained geodesic, 3. compute the initial derivative ∇ s J(0) of the Jacobi field s → J(s) along α verifying J(0) = 0 and J(1) = j, 4. correct the shooting direction u = u + ∇ s J(0).

Recall that the map ϕ :

T α(0) M n+1 → T α(1) M n+1 , ∇ s J(0) → J(1)
associating to the initial derivative ∇ s J(0) of a Jacobi field with initial value J(0) = 0 its end value J(1), is a linear bijection between two vector spaces which can be obtained using Proposition 3.4. Its inverse map can be computed by considering the image of a basis of T c(0) M n+1 .

A discrete analog of unparameterized curves

The final step in building our discrete model is to introduce a discretization of the quotient shape space. There seems to be no natural, intrinsic definition of the shape of a discrete curve, as by definition we are lacking information : we only have access to a finite number n+1 of points. Therefore to introduce our model, we will make the assumption that we know the equations of the underlying curves, that is, that for each discrete curve α, we have access to the shape c of the smooth curve c of which α is the discretization. In applications, if we don't have access to this information, we can set c to be the shape of an optimal interpolation. Recall that α = (x 0 , . . . ,

x n ) is the discretization of size n of t → c(t) if c(k/n) = x k for all 0 ≤ k ≤ n.
For an element c of the shape space S = M/Diff + ([0, 1]), we denote by Disc n (c) := {α ∈ M n+1 : ∃c ∈ π -1 (c), α is the discretization of size n of c},

THE DISCRETE MODEL

the set of its discretizations, i.e. the set of elements of M n+1 that are discretizations of smooth curves with shape c. Recall that π is the natural projection M → S. We denote by M n the set of discrete curves of size n paired up with their underlying shapes

M n := {(α, c) ∈ M n+1 × S : α ∈ Disc n (c)}.
The goal, for two elements (α 0 , c 0 ) and (α 1 , c 1 ) of M n , is to redistribute the n + 1 points on c 1 to minimize the discrete distance to the n + 1 points α 0 on c 0 , i.e. to find the optimal discretization α opt

1 of c 1 α opt 1 = argmin{d n (α 0 , α), α ∈ Disc n (c 1 )}, (3.15)
where d n is the geodesic distance associated to the discrete metric G n . We want to approach this reparameterization α opt 1 using Algorithm 2.4, i.e. by iteratively computing the "horizontal part" of the geodesic linking α 0 to an iteratively improved discretization α1 of c 1 . To define the horizontal part of a path of discrete curves, we need two things : a notion of horizontality, and a notion of reparameterization. Let us start with the former. We define the discrete vertical and horizontal spaces in α as the following subsets of

T α M n+1 Ver n α := {mv : m = (m k ) k ∈ R n+1 , m 0 = m n = 0}, Hor n α := {h ∈ T α M n+1 : G n (h, mv) = 0 ∀m = (m k ) k ∈ R n+1 , m 0 = m n = 0}
, where v = (v k ) k is still defined by (3.7). Similarly to the continuous case, we can show the following result.

(D τ h) k , v k -4 |τ k | |τ k-1 | (D τ h) k-1 , b -1 k-1 v k + ( 1 4 -b -1 k-1 )λ k-1 v k-1 = 0.
Any tangent vector w ∈ T α M n+1 can be uniquely decomposed into a sum w = w ver + w hor where w ver = mv ∈ Ver n α , w hor = wmv ∈ Hor n α and the components (m k ) k verify m 0 = m 1 = 0 and the following recurrence relation

λ k m k+1 -1 + 4 |τ k | |τ k-1 | (b -2 k-1 + λ 2 k-1 ( 1 4 -b -2 k-1 )) m k + |τ k | |τ k-1 | λ k-1 m k-1 = (D τ w) k , v k -4 |τ k | |τ k-1 | b -1 k-1 (D τ w) k-1 , v k + ( 1 4 -b -1 k-1 )λ k-1 (D τ w) k-1 , v k-1 ,
with the notation

λ k := v k+1 , v k .
{y , kp < < (k + 1)p} are distributed according to arc-length on c between y kp and y (k+1)p for all k. In other words, β is the discretization of size N of the only parameterized curve c ∈ π -1 (c) of which α is the discretization of size n and which is parametrized by arc length on the segments c| [k/n,(k+1)/n] , 0 ≤ k ≤ n.

Definition 3.2. We call the discrete curve β ∈ Disc N (c) the refinement of size N of (α, c).

The discrete analogs of the increasing diffeomorphims of the continuous case are defined as increasing injections ϕ : {0, . . . , n} → {0, . . . , N } such that ϕ(0) = 0 and ϕ(n) = N . Their set is denoted by Inj + (n, N ). We then define the discrete analog of reparameterizing as the action :

Inj + (n, N ) × M n → M n , ϕ ((x k ) k , c) := ((y ϕ(k) ) k , c),
where (y k ) k is the refinement of size N of ((x k ) k , c). Note that the action of ϕ is non transitive. This definition of reparameterization in the discrete case simply boils down to redistributing the n + 1 points on c by choosing among the N + 1 points of the refinement of α, while preserving the order and keeping the extremities fixed. We can now define the horizontal part of a path of discrete curves.

ϕ s (s)(k) = m k (s) |nτ k (s)| ∆ϕ(s)(k), k = 0, . . . , n, (3.16
)

with ∆ϕ(s)(k) = N/2(ϕ(s)(k + 1) -ϕ(s)(k -1)), 1 ≤ k ≤ N -1
and where m = (m k ) k is the norm of the vertical component of α (s) and verifies

λ k m k+1 -1 + 4 |τ k | |τ k-1 | (b -2 k-1 + λ 2 k-1 ( 1 4 -b -2 k-1 )) m k + |τ k | |τ k-1 | λ k-1 m k-1 = ∇ s τ k , v k -4 |τ k | |τ k-1 | × b -1 k-1 ∇ s τ k-1 , v k + ( 1 4 -b -1 k-1 )λ k-1 ∇ s τ k-1 , v k-1 .
(3.17) to the different geometries. The flat case allows us to validate our geodesic shooting algorithm. Indeed, we have a closed form for the geodesics in that case as they are simply the projections in M n+1 of the L 2 -geodesics between the square root velocity representations of the curves, as stated in Remark 3.4. We can see in Figure 3.2 that these exact geodesics, shown in red, are very close to the geodesics obtain by geodesic shooting, shown in blue. In the three settings (negative, flat and positive curvature), we can observe that our metric has a tendency to "shrink" the curve as it optimally deforms from one state to another, compared to the L 2 -metric.

We then tested Algorithm 3.2 to obtain optimal matchings between discrete curves. Results obtained on curves in the hyperbolic half-plane are shown in Figure 3.4 and on curves in the sphere in Figure 3.5. We fix a pair of curves that are identical modulo translation and parameterization, and consider 8 different combinations of parameterizations. We always fix the parameterization of the curve on the left-hand side, while searching for the optimal reparameterization of the curve on the right-hand side. The points are either "evenly distributed" along the latter, or along the former. In each case, the geodesic between the initial parameterized curves is shown in blue, and the horizontal geodesic obtained as output of the optimal matching algorithm is shown in red. We can see that for both the hyperbolic plane and the sphere, the red horizontal geodesics obtained as outputs of Algorithm 2.4 redistribute the points along the right-hand side curve in the way that seems natural : similarly to the distribution of the points on the left curve. The length of the regular and horizontal geodesics shown in Figure 3.4 are given in Table 3.1, in the same order as the corresponding images. For both examples, the horizontal geodesics are always shorter than the initial geodesics, as expected, and have always approximatively the same length. This common length is the distance between the underlying shapes. We also find that the underlying shapes of the horizontal geodesics are very similar. Testing other combinations of parameterizations for the same pair of shapes than the one showed in the bottom half of Figure 3.4 and overlapping the obtained horizontal geodesics gives an idea of the geodesic in the shape space, as shown in Figure 3.6.

Proof of Theorem 3.1

We conclude this chapter with the proof of Theorem 3.1. Let us first remind the result.

Theorem 3.1 (Convergence of the discrete model to the continuous model).

Let s → c(s) be a C 1 -path of C 2 -curves with non vanishing derivative with respect to t. This path can be identified with an element (s, t) → c(s, t) of

C 1,2 ([0, 1] × [0, 1], M ) such that c t = 0. Consider the C 1 -path in M n+1 , s → α(s) = (x 0 (s), . . . , x n (s))
, that is the discretization of size n of c. Then there exists a constant λ > 0 that does not depend on c and such that for n large enough,

|E(c) -E n (α)| ≤ λ n (inf |c t |) -1 |c s | 2 2,∞ (1 + |c t | 1,∞ ) 3 ,
where E and E n are the energies with respect to metrics G and G n respec-CHAPTER 3. DISCRETIZATION tively and where

|c t | 1,∞ := |c t | ∞ + |∇ t c t | ∞ , |c s | 2,∞ := |c s | ∞ + |∇ t c s | ∞ + |∇ 2 t c s | ∞ , and |w| ∞ := sup s,t∈[0,1]
|w(s, t)| denotes the supremum over both s and t of a vector field w along c.

Proof of Theorem 1. To prove this result, we introduce the unique path ĉ of piecewise geodesic curves of which α is the n-discretization. It is obtained by linking the points x 0 (s), p 1 (s), . . . , x n (s) of α by pieces of geodesics for all times s ∈ [0, 1]

ĉ(s, k n ) = c(s, k n ) = x k (s), ĉ(s, •)| [ k n , k+1
n ] is a geodesic, for k = 0, . . . , n. Then the difference between the energy of the path of curves E(c) and the discrete energy of the path of discrete curves E n (α) can be controlled in two steps :

|E(c) -E n (α)| ≤ |E(c) -E(ĉ)| + |E(ĉ) -E n (α)|.
Step 1. We first consider the difference between the continuous energies of the smooth and piecewise geodesic curves

|E(c) -E(ĉ)| = 1 0 1 0 |∇ s q(s, t)| 2 -|∇ s q(s, t)| 2 dt ds ≤ 1 0 1 0 |∇ s q(s, t)| 2 -|∇ s q(s, t)| 2 dt ds ≤ 1 0 1 0 (|∇ s q(s, t)| + |∇ s q(s, t)|) • ∇ s q(s, t) t, k n -∇ s q(s, t) t, k n |dt ds.
Note that the parallel transports ∇ s q(s, t) t, k n and ∇ s q(s, t) t, k n are performed along different curvesc(s, •) and ĉ(s, •) respectively. Let us fix s ∈ [0, 1], 0 ≤ k ≤ n and t ∈ k n , k+1 n . From now on we will omit "s" in the notation w(s, t) to lighten notations. Using the notation w (t) := w(t) t, k n to denote the parallel transport of a vector field w from t to k n along its baseline curve, the difference we need to control is

|∇ s q -∇ s q | = |c t | -1 2 (∇ s c t -1 2 ∇ s c t T ) -|ĉ t | -1 2 (∇ s ĉt -1 2 ∇ s ĉt T ) = (∇ s c t -1 2 ∇ s c t T ) (|c t | -1 2 -|ĉ t | -1 2 ) + |ĉ t | -1 2 (∇ s c t -∇ s ĉt ) -1 2 (∇ s c t -∇ s ĉt ) T .
Since |w -1 2 w T | ≤ |w| for any vector w, we can write

|∇ s q -∇ s q | ≤ |∇ s c t | • |c t | -1/2 -|ĉ t | -1/2 + |ĉ t | -1/2 |∇ s c t -∇ s ĉt |. ( 3
|c t (t) t, k n -ĉt (t) t, k n | ≤ |c t (t) t, k n -c t ( k n )| + |c t ( k n ) -ĉt ( k n )|.
The first term is smaller than

1/n • |∇ t c t | ∞ .
To bound the second term, we place ourselves in a local chart

(φ, U ) centered in c( k n ) = c(s, k n ), such that c([0, 1] × [0, 1]) ⊂ U . After identification
with an open set of R d -where d is the dimension of the manifold M -using this chart, we get

|c t ( k n ) -ĉt ( k n )| ≤ c t ( k n ) -n c( k+1 n ) -c( k n ) + ĉt ( k n ) -n c( k+1 n ) -c( k n ) .
Since a geodesic locally looks like a straight line (see e.g. [START_REF] Émery | Fonctions convexes et semimartingales dans une variété[END_REF]) there exists a constant λ 1 such that

ĉt ( k n ) -n(c( k+1 n ) -c( k n )) ≤ λ 1 c( k+1 n ) -c( k n ) 2 ,
and so

|c t ( k n ) -ĉt ( k n )| ≤ 1 2n |c tt | ∞ + λ 1 n |c t | 2 ∞ .
The second derivative in t of the coordinates of c in the chart (U, φ) can be written c tt = ∇ t c t -Γ ij c t i c t j for = 1, . . . , d, and so there exists a constant

λ 2 such that |c tt | ≤ λ 2 |∇ t c t | ∞ + |c t | 2 ∞ , and
c t -ĉ t ≤ λ 3 n |c t | 1,∞ + |c t | 2 1,∞ . (3.19)
This means that for n large enough, we can write e.g.

1 2 inf |c t | ≤ |ĉ t | ≤ 3 2 |c t | ∞ . (3.20)
From (3.19) we can also deduce that

|c t | -1 2 -|ĉ t | -1 2 = |c t | -|ĉ t | |c t | 1 2 + |ĉ t | 1 2 ≤ |c t -ĉ t | |c t | 1 2 + |ĉ t | 1 2 ≤ λ 3 n (inf |c t |) -1 2 |c t | 1,∞ + |c t | 2 1,∞ . (3.21) When n → ∞, a k (1) → 1, b k (1) → 1, and since 0 ≤ nt -k ≤ 1, a k (nt - k) → 1, b k (nt -k) →
|∇ 2 t ĉs | ∞ ≤ 9λ 4 2 |c t | 2 ∞ |c s | ∞ ,
and so the difference (3.22) can be bounded by

|∇ s c t -∇ s ĉt | ≤ 3 2n |∇ 2 t c s | ∞ + 9λ 4 2 |c t | 2 ∞ |c s | ∞ λ 5 n |c s | 2,∞ 1 + |c t | 2 1,∞ , (3.26
|∇ s q -∇ s q | ≤ λ 3 n (inf |c t |) -1 2 |c s | 2,∞ |c t | 1,∞ + |c t | 2 1,∞ + λ 5 √ 2 n (inf |c t |) -1 2 |c s | 2,∞ 1 + |c t | 2 1,∞ , |∇ s q -∇ s q | ≤ λ 6 n (inf |c t |) -1 2 |c s | 2,∞ (1 + |c t | 1,∞ ) 2 , (3.27)
for some constant λ 6 . To conclude this first step, let us bound the sum

|∇ s q| + |∇ s q| = |c t | -1 2 |∇ s c t -1 2 ∇ s c t T | + |ĉ t | -1 2 |∇ s ĉt -1 2 ∇ s ĉt T | ≤ (inf |c t |) -1 2 |∇ t c s | ∞ + √ 2(inf |c t |) -1 2 |∇ t ĉs | ∞ . (3.28)
Taking the derivative according to t on both sides of (3.25), we get since

n|τ k (s)| = |ĉ t (s, k n )|, ∇ t ĉs (t) t, k n = |ĉ t ( k n )|e k (nt -k)c s ( k n ) N + n c s ( k+1 n ) -c s ( k n ) T + n a k (nt-k) b k (1) c s ( k+1 n ) -a k (1)c s ( k n ) N , since derivation of the coefficients give b k (u) = a k (u) and a k (u) = |τ k |e k (u) = 1 n |ĉ t ( k n )|e k (u)
, where

e k (u) =    sinh (|τ k |u) , if K = -1, 0 if K = 0, -sin (|τ k |u) , if K = +1.
Since the coefficients e k (ntk), a k (ntk)/b k (1) and a k (1) are bounded for n large enough, and since |ĉ t | ≤ 3 2 |c t | ∞ , we can write for some constant λ 7 ,

|∇ t ĉs | ∞ ≤ λ 7 (|ĉ t | ∞ |c s | ∞ + |∇ t c s | ∞ ) ≤ 3λ 7 2 |c s | 2,∞ (1 + |c t | 1,∞ ) . (3.29)
Inserting this into (3.28) gives

|∇ s q| + |∇ s q| ≤ (inf |c t |) -1/2 |∇ t c s | ∞ + 3λ 7 √ 2 |c s | 2,∞ (1 + |c t | 1,∞ ) ≤ λ 8 (inf |c t |) -1/2 |c s | 2,∞ (1 + |c t | 1,∞ ).
(3.30)

Finally, we are able to bound the difference between the energies of the smooth and piecewise-geodesic paths by combining Equations (3.27) and (3.30)

|E(c) -E(ĉ)| ≤ λ 6 λ 8 n (inf |c t |) -1 |c s | 2 2,∞ (1 + |c t | 1,∞ ) 3 .
Step 2. Let us now consider the difference of energy between the path of piecewise geodesic curves and the path of discrete curves. Since ∇ s q k (s) = ∇ s q(s, k n ) for all s ∈ [0, 1] and 0 ≤ k ≤ n, we can write

|E(ĉ) -E n (α)| = 1 0 1 0 |∇ s q(s, t)| 2 dt - 1 n n-1 k=0 |∇ s q k (s)| 2 ds ≤ n-1 k=0 1 0 k+1 n k n |∇ s q(s, t)| 2 -|∇ s q(s, k n )| 2 dt ds ≤ n-1 k=0 1 0 k+1 n k n |∇ s q(s, t)| + |∇ s q(s, k n )| •
∇ s q(s, t) t, k n -∇ s q(s, k n ) dt ds. We fix once again s ∈ [0, 1], 0 ≤ k ≤ n and t ∈ k n , k+1 n . As in step 1, we will omit "s" in most notations. Since |ĉ t (t)| = |ĉ t ( k n )|, we get

∇ s q(t) t, k n -∇ s q( k n ) ≤ |ĉ t ( k n )| -1 2 ∇ s ĉt (t) t, k n -∇ s ĉt ( k n ) -1 2 ∇ s ĉt (t) t, k n -∇ s ĉt ( k n ) T ≤ |ĉ t ( k n )| -1 2 ∇ s ĉt (t) t, k n -∇ s ĉt ( k n ) .
Considering once again the Jacobi field

J(u) := ĉs ( k+u n ), u ∈ [0, 1], along the geodesic γ(u) = ĉ( k+u n ), Equation (3.24) gives ĉs (t) t, k n = c s ( k n ) T + a k (nt -k)c s ( k n ) N + (t -k n )∇ t ĉs ( k n ) T + b k (nt -k) 1 n ∇ t ĉs ( k n ) N .
Recall that b k (u) = a k (u) and a k (u) = |τ k |e k (u), and so taking the derivative with respect to t and decomposing

∇ t ĉs ( k n ) T = ∇ t ĉs ( k n ) -∇ t ĉs ( k n ) N , we obtain ∇ t ĉs (t) t, k n -∇ t ĉs ( k n ) = |ĉ t ( k n )|e k (nt -k) • c s ( k n ) N + a k (nt -k) -1 ∇ t ĉs ( k n ) N .
Noticing that e k (nt-k) (nt-k)|τ k | → 1 and a k (nt-k)-1 (nt-k)|τ k | → 0 when n → ∞, we can deduce that for n large enough,

|e k (nt -k)| ≤ 2(nt -k)|τ k | ≤ 2|τ k | = 2 n |c t | ≤ 2 n |c t | ∞ , |a k (nt -k) -1| ≤ (nt -k)|τ k | ≤ |τ k | = 1 n |c t | ≤ 1 n |c t | ∞ . This gives k+1 n k n ∇ t ĉs (s, t) t, k n -∇ t ĉs (s, k n 2 )| dt ≤ 2 n 2 |c t | 2 ∞ |c s | ∞ + |c t | ∞ |∇ t ĉs | ∞ .
Recall from (3.29) and (3.30) that

|∇ t ĉs | ∞ ≤ 3λ 7 2 |c s | 2,∞ (1 + |c t | 1,∞ ) , |∇ s q| ∞ ≤ 3λ 7 √ 2 (inf |c t |) -1 2 |c s | 2,∞ (1 + |c t | 1,∞ ) ,
and so

k+1 n k n |∇ s q(t)| + |∇ s q( k n )| • |∇ s q(t) t, k n -∇ s q( k n )| dt ≤ 2|∇ s q| ∞ √ 2(inf |c t |) -1 2 k+1 n k n ∇ t ĉs (t) t, k n -∇ t ĉs ( k n )| dt ≤ 6λ 7 (inf |c t |) -1 |c s | 2,∞ (1 + |c t | 1,∞ ) 2 n 2 |c t | 2 ∞ |c s | ∞ + |c t | ∞ 3λ 7 2 |c s | 2,∞ (1 + |c t | 1,∞ ) ≤ λ 9 n 2 (inf |c t |) -1 |c s | 2 2,∞ (1 + |c t | 1,∞ ) 3 .
Finally, we obtain

|E(ĉ) -E n (α)| ≤ λ 9 n (inf |c t |) -1 |c s | 2 2,∞ (1 + |c t | 1,∞ ) 3 ,
which completes the proof.

Chapter 4

Radar applications

Abstract

This final chapter is dedicated to radar applications. We represent locally stationary radar signals by time series of Toeplitz matrices, corresponding to the covariant matrices of their stationary portions. We equip the space of hermitian positive definite matrices with a metric defined as the hessian of the entropy which is in fact the Fisher information metric, and project it on the submanifold of matrices with Toeplitz structure. In the autoregressive coordinates of Toeplitz matrices -used by Burg in the context of maximum entropy spectral analysis -this metric becomes a product metric on the Poincaré polydisk. Statistical analysis of locally stationary radar signals can then be achieved by exploiting the Riemannian structure of the polydisk. As an example, we compute the mean of simulated helicopter signatures.

Geometric approach to radar signal processing

As mentioned in the introduction, Barbaresco proposed in [5] a new representation tool for radar signals based on information geometry and Burg's maximum entropy spectral analysis [START_REF] Parker | Maximum entropy spectral analysis[END_REF], with the aim of improving the resolution with respect to the classical CFAR (Constant False Alarm Rate) methods. In this model, a signature z = (z 1 , . . . , z n ) T ∈ C n measured for a given distance cell after a burst of n pulses, is assumed to be the realization of a centered, stationary circularly-symmetric [START_REF] Robert | Circularly-symmetric gaussian random vectors[END_REF] Gaussian vector Z = (Z 1 , . . . , Z n ) T which is therefore entirely described by its covariance matrix Σ = E(ZZ * ), where E is the expected value and Z T , Z * are the transpose and transconjugate of Z. Recall that the stationarity hypothesis implies 86 CHAPTER 4. RADAR APPLICATIONS the Toeplitz structure of Σ, since the autocorrelation values r k = E(Z i Z i+k ), k = 0, . . . , n -1, depend only on the lag k

Σ =       r 0 r 1 • • • r n-1 r 1 . . . . . . . . . . . . . . . . . . r 1 r n-1 • • • r 1 r 0       .
Representing each signature z = (z 1 , . . . , z n ) ∈ C n by the Toeplitz covariance matrix of the underlying signal, we can compare, average and do statistics on a set of radar signatures using a metric structure on the set T + n of Toeplitz, hermitian, positive, definite matrices of size n. The choice of a Riemannian structure on the set of parameters (Toeplitz covariance matrices) of a family of probability distributions (centered stationary multivariate Gaussians) is at the heart of information geometry.

Information geometry of multivariate Gaussians

Let us put aside the stationarity and the associated Toeplitz structure for now. The metric that we consider on the space H + n of hermitian positive definite matrices is the hessian of minus the entropy. Barbaresco has shown in [START_REF] Barbaresco | Koszul information geometry and souriau geometric temperature/capacity of lie group thermodynamics[END_REF] that it is also the Legendre dual of the Fisher Information metric on the set of dual coordinates -which as we will see are the inverse covariance matrices in the case of multivariate Gaussians -and as such it defines the same metric structure [START_REF] Amari | Information geometry of divergence functions[END_REF]. For the sake of completeness we remind the proof in this section. The dual differential geometry of multivariate Gaussians with zero mean was investigated in [START_REF] Ohara | Dualistic differential geometry of positive definite matrices and its applications to related problems[END_REF], and for non-zero mean it can be found in [START_REF] Barbaresco | Geometric theory of heat from souriau lie groups thermodynamics and koszul hessian geometry: Applications in information geometry for exponential families[END_REF]. Consider the centered multivariate circularly-symmetric Gaussian density N (0, Σ) of covariance matrix Σ ∈ H + n defined for z ∈ C n by

p Σ (z) = 1 π n det Σ exp -z * Σ -1 z . (4.1)
The simpler form of the complex multivariate Gaussian with respect to the real multivariate Gaussian is due to Gauss' integral taking a more concise form over the complex plane than the real line Its entropy is given by H(Σ) := -E Σ ln p Σ (Z) = n ln(πe) + ln(det Σ).

GEOMETRIC APPROACH

87

We consider the following metric on H + n , where each element Σ is identified with its vectorization Vec(Σ) ∈ C n 2 , the column vector obtained by stacking the columns of Σ on top of one another,

ds 2 = -dΣ * Hess H(Σ)dΣ. (4.2)
This metric can be written as the Legendre dual of the Fisher metric on a different space of coordinates. To see this, let us write the density (4.1) under the canonical form

p Σ (z) = exp (-η * T (z)) C n exp (-η * T (z)) dz , z ∈ C n , (4.3)
where the parameters η and T (z) are the following vectorizations in

C n 2 η = 1 2 Vec(Σ -1 ), T (z) = Vec(zz * ).
This writing is the exponential form of a centered multivariate gaussian except for the minus that we choose to leave out of the variables. We define the following potential, the logarithm of the so-called Koszul-Vinberg Characteristic Function (KVCF) [START_REF] Barbaresco | Koszul information geometry and souriau geometric temperature/capacity of lie group thermodynamics[END_REF],

φ(η) := -ln C n exp(-η * T (z)) dz = -n ln(π) -ln(det Σ).
Then, noticing that the gradient of φ with respect to η is given by

θ := Grad η φ = E Σ T (Z) = Σ,
and since the logarithm of the density (4.3) gives for all z

∈ C n ln p Σ (z) = -η * T (z) + φ(η), (4.4)
we see that the entropy is simply the Legendre dual of φ

H(Σ) = η * θ -φ(η) =: ψ(θ),
and so the metric (4.2) that we put on H + n can be written as minus the hessian of the Legendre dual ψ of the KVCF φ

ds 2 = - ∂ 2 ψ(θ) ∂ θi ∂θ j d θi dθ j .
On the other hand, we know that if η → p η (z) is twice differentiable for all z -which is the case here -then the Fisher information can also be written

I(η) i,j = -E η ∂ 2 ∂ ηi ∂η j ln p Σ (Z) ,
where the η i 's are the coordinates of η ∈ C n 2 , and so we deduce from (4.4) that the Fisher information matrix is equal to minus the hessian of the KVCF

I(η) = -Hess φ(η).
We obtain the following result Proposition 4.1. Let Σ ∈ H + n . The dual potentials defined in the dual coordinate systems η = 1 2 Σ -1 and θ = Σ,

φ(η) = -n ln(π) -ln(det Σ), ψ(θ) = n ln(πe) + ln(det Σ) = H(Σ),
where H(Σ) is the entropy, define dual hessian metrics

g φ ij (η) = - ∂ 2 φ(η) ∂ ηi ∂η j , g ψ ij (θ) = - ∂ 2 ψ(θ) ∂ θi ∂θ j
which are the same tensor expressed in different coordinate systems and define the same local distance ds 2 φ = ds 2 ψ .

They both coincide with the Fisher metric on H + n .

Proof. We have already seen that the Fisher information matrix coincides with minus the hessian of the KVCF, and so it is straightforward that g φ defines the Fisher metric

ds 2 φ = i,j g φ ij (η)dη i dη j = i,j I(η) ij dη i dη j ,
Let us remind why metric g ψ , defined from the dual potential ψ, represents the same metric. From θ = Grad η φ we obtain

dθ i = n j=1 ∂ 2 φ(η) ∂ ηi ∂η j dη j = - n j=1 g φ ij (η)dη j ,

GEOMETRIC APPROACH

89

and an analogous expression for dη i using η = Gradθ ψ. In matrix form, dθ = g φ dη, dη = g ψ dθ,

and so we get that the matrices g φ and g ψ are inverse of one another and since they are hermitian,

ds 2 ψ = d θT g ψ dθ = dθ T dη = (g φ dη) T dη = dη T g φ dη = ds 2 φ .
Now let us recall the stationary hypothesis of the Gaussian signals. The covariance matrices that we consider have an additional Toeplitz structure and live in the submanifold T + n of H + n . As we will see in the following section, a Toeplitz matrix can be parameterized by n coefficients of the product space R + ×D n-1 , where D is the unit complex disk. In order to respect the Toeplitz structure, we consider the metric induced by the Fisher metric (4.2) in the submanifold T + n .

Autoregressive parameterization of Toeplitz matrices

Consider an element of T + n , i.e. a Toeplitz, hermitian, positive definite matrix of size n

Σ n =       r 0 r 1 • • • r n-1 r 1 . . . . . . . . . . . . . . . . . . r 1 r n-1 • • • r 1 r 0      
, and denote the extracted Toeplitz matrices of size k = 1, . . . , n -1 by

Σ k = (r j-i ) 1≤i,j≤k , with r -i := r i .
In his thesis [START_REF] Parker | Maximum entropy spectral analysis[END_REF], Burg showed that the stationary process Z , ∈ Z, of maximum entropy -i.e. that adds the least amount of information -under the k + 1 autocorrelation constraints E(Z Z +i ) = r i , i = 0, . . . , k, is an autoregressive process of order k, i.e. such that any component can be expressed as a linear combination of the k previous ones

Z = - k i=1 a (k) i Z -i + ε (k) , (4.5)
where the ε (k) 's are independent identically distributed centered Gaussians of variance P k ∈ R * + . The last coefficient of the autoregressive model of order k is called the k th reflection coefficient and is denoted by

µ k := a (k) k , k = 1, . . . , n -1.
From (4.5) we can easily deduce the matrix form of the Yule-Walker equations, which gives a relation between R k+1 and the coefficients of the autoregressive process (a

(k) 1 , . . . , a (k) k , P k )       r 0 r 1 • • • r k r 1 . . . . . . . . . . . . . . . . . . r 1 r k • • • r 1 r 0       •       1 a (k) 1 . . . a (k) k       =       P k 0 . . . 0       . (4.6)
Computing the determinant of Σ k+1 using the elements of its last column and the associated minors, as well as the k + 1 relations given by (4.6), gives

det Σ k+1 = (-1) 2(k+1) r 0 + a (k) 1 r 1 + . . . + a (k) k det Σ k , (4.7)
and so we obtain

P k = det Σ k+1 det Σ k .
The autoregressive coefficients (a

(k) 1 , . . . , a (k)
k ), or equivalently, the reflexion coefficients (µ 1 , . . . , µ k ), can be obtained by the classical Levinson recursion: the coefficients of order k + 1 are deduced from those of order k by noticing that

      r 0 r 1 • • • r k+1 r 1 . . . . . . . . . . . . . . . . . . r 1 r k+1 • • • r 1 r 0       •         1 0 a (k) 1 ā(k) k . . . . . . a (k) k ā(k) 1 0 1         =        P k e k 0 0 . . . . . . 0 0 e k P k        , where e k = r k+1 + a (k) k r k + . . . + a (k) 1 r 1 .
The trick is to apply the leftmultiplication on both sides by a 2 × 2 matrix M such that the right-hand side becomes

       P k e k 0 0 . . . . . . 0 0 e k P k        M =        P k+1 0 0 0 . . . . . . 0 0 0 P k+1        .
Removing the zeros in the center on both sides, we see that M is given by

M = P k+1 P 2 k -|e k | 2 P k -e k -e k P k .
We obtain the equality

        1 0 a (k) 1 ā(k) k . . . . . . a (k) k ā(k) 1 0 1         M =         1 µ k+1 a (k+1) 1 ā(k+1) k . . . . . . a (k+1) k ā(k+1) 1 µ k+1 1        
, which gives the equations of the Levinson recursion, that is P 0 = r 0 and for k = 0, . . . , n -2,

µ k+1 = - r k+1 + a (k) k r k + . . . + a (k) 1 r 1 P k , a (k+1) i = a (k)
i + µ k+1 • ā(k) k+1-i , i = 1, . . . , k, P k+1 = P k (1 -|µ k+1 | 2 ).
Φ : T + n → R * + × D n-1 , Σ n → (P 0 , µ 1 , . . . , µ n-1 ).
From (4.8) it is easy to see that this map is a bijection, and in fact it is a diffeomorphism between the two submanifolds of R 2n-1 [START_REF] Verblunsky | On positive harmonic functions[END_REF], [START_REF] William F Trench | An algorithm for the inversion of finite toeplitz matrices[END_REF]. This coordinate system for Toeplitz matrices is very convenient as it gives a product form to the projection of the Fisher metric (4.2) in T + n . Recall that this metric is given by minus the hessian of the entropy H(Σ) = n ln(πe) + ln(det Σ).

In the reflection coefficient coordinates (β 1 , β 2 , . . . , β n ) = (P 0 , µ 1 , . . . , µ n-1 ) this is written

ds 2 = - ∂ 2 ln(det Σ) ∂ βi ∂β j d βi dβ j , 1 ≤ i, j ≤ n.
observations, or performing classification. This can be very useful in target detection as well as target recognition. As an example, we compute the mean of a set of simulated helicopter signatures. We choose to focus on the second reflection coefficient and compute the Fréchet mean μ2 of m curves tracing the evolutions t → µ k 2 (t), 1 ≤ k ≤ m of the second reflection coefficient of m signals. The Fréchet mean, also called intrinsic mean, of a set of points in a metric space is defined as the point that minimizes the sum of the squared distances to each point of the set. Here, the metric space is the infinite-dimensional manifold M = Imm([0, 1], D) of open immersions in the Poincaré disk equipped with metric (2.2), and the mean is

μ2 = argmin µ∈M m k=1 d µ, µ k 2 2 ,
if d the geodesic distance on M. Since it is defined as a minimizer of a functional, this intrinsic mean can be found by a gradient descent type procedure called a Karcher flow, summarized as follows.

Algorithm 4.1 (Mean of a set of curves). Input : (µ k 2 (t), 1 ≤ k ≤ m). Initialize μ2 and repeat until convergence :

1. for k = 1, . . . , m, compute the geodesic γ k (s) linking μ2 to µ k 2 using geodesic shooting (Algorithm 2.2) and its initial tangent vector γ (j) s (0), 2. update the mean μ2 in the direction of the sum of the initial speed vectors using the exponential map (Algorithm 2.1)

μ2 ← exp M μ2 1 m m k=1 γ k s (0) .
Output : μ2 (t).

We give results for 4 sets of m = 11 curves. These signals are generated using the helicopter signature simulator, and correspond to the observation at 4 different times of 11 helicopters which differ only in their rotor rotation speeds. We consider small variations (less than 1%) around the mean value of 390 RPM (rotations per minute), and show the obtained curves in Figure 4.3. Theses curves are shown in the hyperbolic half-plane representation, which is equivalent to the Poincaré disk in terms of geometry. In each case, the red extremity of the colormap corresponds to the helicopter with the highest rotation speed, the blue extremity to the lowest rotation speed, and the mean curve is shown in black. This can be used to construct a "reference signature" for a given type of helicopter, for target recognition purposes. by w k+1 the vector field along the curve x k obtained by parallel transporting back the vector w k+1 (s) along the geodesic γ k (s, •) for all s ∈ [0, 1], i.e. w k+1 (s) = P 1,0 γ k (s,•) (w k+1 (s)). We have

∇ s f (-) k (w k+1 ) = ∇ s f k (w k+1 ) = ∇ s f k (w k+1 ) + f k ∇ s (w k+1 ) , (A.1)
and so we need to compute ∇ s (w k+1 ). Let V (s, t)

:= P 1,t γ k (s,•) (w k+1 ) so that ∇ s V (s, 1) = ∇ s w k+1 and ∇ s V (s, 0) = ∇ s (w k+1 ), then ∇ s V (s, 1) 1,0 = ∇ s V (s, 0) + 1 0 ∇ t ∇ s V (s, t) t,0 dt, = ∇ s V (s, 0) + 1 0 R(∂ t γ k t,0 , ∂ s γ k t,0 )V (s, t) t,0 dt, since ∇ t V = 0, and where ∂ t γ k (s, t) t,0 = τ k (s). We get, since ∇R = 0, (∇ s w k+1 ) = ∇ s (w k+1 ) + R τ k , 1 0 ∂ s γ k t,0 dt (w k+1 ). (A.2)
To find an expression for ∂ s γ k t,0 , we consider the Jacobi field J(t) := ∂ s γ k (s, t) along the geodesic t → γ k (s, t). The vector field J verifies J(0) = x k (s), J(1) = x k+1 (s), ∇ t J(0) = ∇ s τ k (s), where the last equality results from the inversion ∇ t ∂ s γ k (s, 0) = ∇ s ∂ t γ k (s, 0) and ∂ t γ k (s, 0) = τ k (s). Applying Lemma 3.1 gives, for all k = 0, . . . , n -1, ∂ s γ k (s, t) t,0 = x k (s) T + a k (s, t) x k (s) N + t ∇ s τ k (s) T + b k (s, t)∇ s τ k (s) N . where Y k is defined by

Y k = (x k ) T + b k (x k ) N + 1 2 ∇ s τ k T + K 1 -a k |τ k | 2 ∇ s τ k N ,
where we integrate by parts to obtain the third line from the second. The goal is to express ∂ a x k in terms of ∂ a x 0 and ∇ a q , = 0, • • • , k. That way, the only elements that depend on a once we take a = 0 are (∂ a x 0 , ∇ a q 0 , • • • , ∇ a q n-1 ) which can be chosen independently to be whatever we want. Let us fix 0 ≤ k ≤ n -1 and s ∈ [0, 1] and consider the path of geodesics a → γ k (a, •) such that γ k (a, 0) = x k (a, s), γ k (a, 1) = x k+1 (a, s) and ∂ t γ k (a, 0) = τ k (a, s) = log x k (a,s) (x k+1 (a, s)). Then by definition, for each a ∈ [0, 1], t → J(a, t) := ∂ a γ k (a, t) is a Jacobi field along the geodesic t → γ k (a, t) of M , and so Lemma 3.1 gives

∂ a x k+1 = ∂ a x k T + a k ∂ a x k N + ∇ a τ k T + b k ∇ s τ k N , (A.4)
where ∂ a x k+1 denotes the parallel transport of ∂ a x k+1 from x k+1 (s) to x k (s) along the geodesic. Differentiation of q k = √ n τ k /|τ k | gives

∇ s q k = √ n |τ k | -1/2 ∇ s τ k -1 2 ∇ s τ k T ,
and taking the tangential part on both sides yields ∇ s q k T = √ n |τ k | -1/2 1 2 ∇ s τ k T , and so finally

∇ s τ k = |τ k | 1/2 / √ n ∇ s q k + ∇ s q k T = |q k |/n ∇ s q k + ∇ s q k T .
Injecting this in (A.4) and noticing that f k (w), z = w, f k (z) and g k (w), z = w, g k (z) for any pair of vectors w, z gives ∂ a x k+1 = f k (∂ a x k ) + 1 n g k (∇ a q k ), (A.5)

w k+1 , ∂ a x k+1 = f (-) k (w k+1 ), ∂ a x k + 1 n g (-)
k (w k+1 ), ∇ a q k , (A.6)

for any tangent vector w k+1 ∈ T x k+1 M . From equation (A.6) we can deduce, for k = 1, . . . , n,

w k , ∂ a x k = f (-) 0 • • • • • f (-) k-1 (w k ), ∂ a x 0 + 1 n k-1 =0 g (-) • f (-) +1 • • • • • f (-) k-1 (w k ), ∇ a q .
With the notation R k := R(q k , ∇ s q k )x k we get

R k , ∂ a x k = f (-) 0 • • • • • f (-) k-1 (R k ), ∂ a x 0 + 1 n k-1 =0 g (-) • f (-) +1 • • • • • f (-) k-1 (R k ), ∇ a q ,
and we can then write the derivative of the energy (E n ) (0) for a = 0 in the following way

- 1 0 ∇ s x 0 + 1 n n-1 k=0 f (-) 0 • • • • • f (-) k-1 (R k ), ∂ a x 0 + 1 n 2 n-1 k=1 k-1 =0 g (-) • f (-) +1 • • • • • f (-) k+1 (R k ), ∇ a q + 1 n n-1 k=0 ∇ 2 s q k , ∇ a q k ds,
where in the first sum we use the notation convention

f 0 • • • • • f -1 := Id.
Noticing that the double sum can be rewritten

n-2 =0 n-1 k= +1 g (-) • f (-) +1 • • • • • f (-) k-1 (R k ), ∇ a q ,
we obtain for (E n ) (0)

- 1 0 ∇ s x 0 + 1 n n-1 k=0 f (-) 0 • • • • • f (-) k-1 (R k ), ∂ a x 0 + 1 n n-1 k=0 ∇ 2 s q k + 1 n n-1 =k+1 g (-) k • f (-) k+1 • • • • • f (-)
-1 (R ), ∇ a q k ds,

(A.7)
where in the last sum we use the convention n-1 =n = 0. Since this quantity has to vanish for any choice of (∂ a x 0 (0, •), ∇ a q 0 (0, •), . . . , ∇ a q n-1 (0, •)), the geodesic equations for the discrete metric are

∇ s x 0 + 1 n n-1 k=0 f (-) 0 • • • • • f (-) k-1 (R k ) = 0, ∇ 2 s q k + 1 n n-1 =k+1 g (-) k • f (-) k+1 • • • • • f (-) -1 (R ) = 0,
for all k = 0, . . . , n -1, with the conventions n-1 =n = 0 and s → α(s) = (x 0 (s), . . . , x n (s)) be a geodesic path in M n+1 . For all s ∈ [0, 1], the coordinates of its acceleration ∇ s α (s) can be iteratively computed in the following way

f 0 • • • • • f -1 := Id.
∇ s x 0 = - 1 n R 0 + f (-) 0 (R 1 ) + . . . + f (-) 0 • • • • • f (-) n-2 (R n-1 ) , ∇ s x k+1 = ∇ s f k (x k ) + f k (∇ s x k ) + 1 n ∇ s g k (∇ s q k ) + 1 n g k (∇ 2 s q k ) + R(τ k , Y k )(x k+1 ),
for k = 0, . . . , n -1, where the R k 's are defined as in Proposition 3.2, the symbol • denotes the parallel transport from x k+1 (s) back to x k (s) along the geodesic linking them, the maps ∇ s f k and ∇ s g k are given by Lemma 3.2, Y k is given by Equation (3.13) and

∇ s τ k = (D τ α ) k , ∇ s v k = 1 |τ k | ∇ s τ k -∇ s τ k T , ∇ s q k = n |τ k | ∇ s τ k - 1 2 ∇ s τ k T , ∇ 2 s q k = - 1 n g (-) k R k+1 + f (-) k+1 (R k+2 ) + . . . + f (-) k+1 • • • • • f (-)
n-2 (R n-1 ) .

Proof. For all s ∈ [0, 1], we initialize ∇ s x k (s) for k = 0 using the first geodesic equation in (3.14); the difficulty lies in deducing ∇ s x k+1 (s) from ∇ s x k (s). Just as we have previously obtained (A.5), we can obtain by replacing the derivatives with respect to a by derivatives with respect to s x k+1 = x k T + a k x k N + ∇ s q k T + b k ∇ s τ k N , (A.8)

x k+1 = f k (x k ) + 1 n g k (∇ s q k ),
and by differentiating with respect to s

∇ s x k+1 = ∇ s f k (x k ) + f k (∇ s x k ) + 1 n ∇ s g k (∇ s q k ) + 1 n g k (∇ 2 s q k ).
(A.9) We have already computed (A.3) the covariant derivative of a vector field s → w k+1 (s) ∈ T x k (s) M and so we can write

∇ s x k+1 = ∇ s x k+1 + R(Y k , τ k ) x k+1 ,
where Y k is defined by Equation (3.13). Together with Equation (A.9), this gives the desired equation for ∇ s x k+1 . Finally, ∇ s τ k = (D τ α ) k results directly from (A.8), ∇ 2 s q k is deduced from the second geodesic equation and the remaining equations follow from simple computation. Proposition 3.4 (Discrete Jacobi fields) Let s → α(s) = (x 0 (s), . . . , x n (s)) be a geodesic path in M n+1 , s → J(s) = (J 0 (s), . . . , J n (s)) a Jacobi field along α, and (-δ, δ) a → α(a, •) a corresponding family of geodesics, in the sense just described. Then J verifies the second order linear ODE ∇ 2 s J 0 = R(x 0 , J 0 )x 0 -

1 n n-2 k=0 k =0 f (-) 0 • • • • • ∇ a f (-) • • • • • f (-) k (R k+1 ), - 1 n ∇ a R 0 + f (-) 0 (∇ a R 1 ) + . . . + f (-) 0 • • • • • f (-) n-2 (∇ a R n-1 ) ∇ 2 s J k+1 = f k (∇ 2 s J k ) + 2∇ s f k (∇ s J k ) + ∇ 2 s f k (J k ) + 1 n g k (∇ 2 s ∇ a q k ) + 2 n ∇ s g k (∇ s ∇ a q k ) + 1 n ∇ 2 s g k (∇ a q k ) + R(τ k , Y k ) R(Y k , τ k )(J k+1 ) + R(∇ s τ k , Y k )(J k+1 ) + R(τ k , ∇ s Y k )(J k+1 ) + 2R(τ k , Y k )(∇ s J k+1 ),
for all 0 ≤ k ≤ n -1, where R k := R(q k , ∇ s q k )x k and the various covariant derivatives according to a can be expressed as functions of J and ∇ s J,

∇ a R k = R ∇ a q k , ∇ s q k x k + R q k , ∇ s ∇ a q k + R(J, x k )q k x k + R q k , ∇ s q k )∇ s J k , ∇ a τ k = (D τ J) k , ∇ a v k = |τ k | -1 ∇ a τ k -∇ a τ k T , ∇ a q k = n/|τ k | ∇ a τ k -1 2 ∇ a τ k T , ∇ s ∇ a q k = n g k -1 (∇ s J k+1 ) + R(Y k , τ k )(J k+1 ) -∇ s f k (J k ) -f k (∇ s J k ) + n ∇ s g k -1 J k+1 -f k (J k ) , ∇ s ∇ s ∇ a q k = R(∇ s x k , J k )q k + R(x k , ∇ s J k )q k + 2R(x k , J k )∇ s q k - 1 n n-1 =k+1 g (-) k • f (-) k+1 • • • • • f (-) -1 (∇ a R ) - 1 n n-1 =k+1 -1 j=k g (-) k • • • • • ∇ a f (-) j • • • • • f (-) -1 (R ), ∇ s Y k = (∇ s x k ) T + b k (∇ s x k ) N + (1 -b k ) x k , ∇ s v k v k + x k , v k ∇ s v k + ∂ s b k (x k ) N + 1 2 (∇ s ∇ s τ k ) T + K 1 -a k |τ k | 2 (∇ s ∇ s τ k ) N + ∂ s K 1 -a k |τ k | 2 (∇ s τ k ) N + 1 2 -K 1 -a k |τ k | 2 ( ∇ s τ k , ∇ s v k v k + ∇ s τ k , v k ∇ s v k ),
with the notation conventions f

∇ a f (-) k (w) = (∇ a f k ) (-) (w) + f k R(Z k , τ k )(w k+1 ) , ∇ a g (-) k (w) = (∇ a g k ) (-) (w) + g k R(Z k , τ k )(w k+1 ) , ∇ s g k -1 (w) = ∂ s |q k | -1 |q k |g k -1 (w) +|q k | -1 ∂ s (b -1 k ) w N + |q k | -1 1/2 -b -1 k w, ∇ s v k v k + w, v k ∇ s v k , and
Z k = J k T + b k J k N + 1 2 ∇ a τ k T + K 1 -a k |τ k | 2 ∇ a τ k N .
Proof. For all a ∈ (-δ, δ), α(a, •) verifies the geodesic equations (3.14). Tak-ing the covariant derivative of these equations according to a we obtain

∇ a ∇ s ∂ s x 0 + 1 n n-1 k=0 ∇ a f (-) 0 • • • • • f (-)
k-1 R(q k , ∇ s q k )∂ s x k = 0, (A.10)

∇ a ∇ 2 s q k + 1 n n-1 =k+1 ∇ a g (-) k • f (-) k+1 • • • • • f (-)
-1 R(q , ∇ s q )∂ s x = 0.

(A.11)

Since for a = 0, ∇ a ∇ s ∂ s x 0 = ∇ 2 s J 0 + R(J 0 , ∂ s x 0 )∂ s x 0 , we get

∇ 2 s J 0 = R(∂ s x 0 , J 0 )∂ s x 0 - 1 n n-1 k=0 ∇ a f (-) 0 • • • • • f (-) k-1 R(q k , ∇ s q k )∂ s x k ,
and the differentiation

∇ a f (-) 0 • • • • • f (-) k-1 (R k ) = f (-) 0 • • • • • f (-) k-1 (∇ a R k ) + k-1 =0 f (-) 0 • • • • • ∇ a f (-) • • • • • f (-) k-1 (R k )
gives the desired equation for ∇ 2 s J 0 . Now we will try to deduce ∇ 2 s J k+1 from (A.11). If J k+1 (s) denotes the parallel transport of the vector J k+1 (s) from x k+1 (s) back to x k (s) along the geodesic that links them, we know from (A.5) that J k+1 = f k (J k ) + 1 n g k (∇ a q k ). (A.12)

We also know from (A.3) that (∇ s J k+1 ) = ∇ s (J k+1 ) + R(τ k , Y k )(J k+1 ), (A. [START_REF] Bauer | Sobolev metrics on shape space of surfaces[END_REF] and by iterating

(∇ 2 s J k+1 ) = ∇ s (∇ s J k+1 ) + R τ k , Y k (∇ s J k+1 ) = ∇ 2 s (J k+1 ) + ∇ s R(τ k , Y k )(J k+1 ) + R τ k , Y k (∇ s J k+1
) Developping and injecting Equation (A.12) in the latter gives

(∇ 2 s J k+1 ) = ∇ 2 s f k (J k ) + 1 n ∇ 2 s g k (∇ a q k ) + R(∇ s τ k , Y k )(J k+1 ) + R(τ k , ∇ s Y k )(J k+1 ) + R(τ k , Y k )(R(Y k , τ k )(J k+1 )) + 2R τ k , Y k (∇ s J k+1 ) .
Developping the covariant derivatives ∇ 2 s f k (J k ) and ∇ 2 s g k (∇ a q k ) gives the desired formula. Now let us explicit the different terms involved in these differential equations. Since ∇R = 0 and ∇ a ∂ s x k = ∇ s ∂ a x k , we have ∇ a R k = R(∇ a q k , ∇ s q k )∂ s x k + R(q k , ∇ a ∇ s q k )∂ s x k + R(q k , ∇ s q k )∇ s J k = R ∇ a q k , ∇ s q k x k + R q k , ∇ s ∇ a q k + R(J, x k )q k x k + R q k , ∇ s q k )∇ s J k .

By taking the inverse of (A.12) we get

∇ a q k = ng -1 k J k+1 -f k (J k ) ,
and taking the derivative according to s on both sides and injecting Equation (A.13) gives

∇ s ∇ a q k = n g k -1 (∇ s J k+1 ) + R(Y k , τ k )(J k+1 ) -∇ s f k (J k ) -f k (∇ s J k ) + n ∇ s g k -1 J k+1 -f k (J k ) .
To obtain ∇ 2 s ∇ a q k , notice that ∇ 2 s ∇ a q k = ∇ s ∇ a ∇ s q k + ∇ s R(∂ s x k , J k )q k , = ∇ a ∇ 2 s q k + R(∂ s x k , J k )∇ s q k + ∇ s R(∂ s x k , J k )q k , and injecting Equation (A.11) with

∇ a g (-) k • f (-) k+1 • • • • • f (-) -1 (R ) = g (-) k • f (-) k+1 • • • • • f (-) -1 (∇ a R ) + -1 j=k g (-) k • • • • • ∇ a f (-) j • • • • • f (-) -1 (R ),
gives us the desired formula. ∇ s Y k results from simple differentiation, and differentiating the maps f (-) k and g (-) k with respect to a is completely analogous to the the computations of Lemma 3.3. Finally, the inverse of g k is given by g k -1 :

T x k M → T x k M , g k -1 : w → |q k | -1 b -1 k w + 1 2 -b -1 k w T ,
and since

∇ s (w T ) = (∇ s w) T + w, ∇ s v k v k + w, v k ∇ s v k ,
it is straightforward to verify that ∇ s g k -1 (w) = ∇ s g k -1 (w) -g k -1 (∇ s w)

Résumé : Probabilités sur les espaces de chemins et dans les espaces métriques associés via la géométrie de l'information ; applications radar Nous nous intéressons à la comparaison de formes de courbes lisses prenant leurs valeurs dans une variété riemannienne M . Dans ce but, nous introduisons une métrique riemannienne invariante par reparamétrisations sur la variété de dimension infinie des immersions lisses dans M . L'équation géodésique est donnée et les géodésiques entre deux courbes sont construites par tir géodésique. La structure quotient induite par l'action du groupe des reparamétrisations sur l'espace des courbes est étudiée. À l'aide d'une décomposition canonique d'un chemin dans un fibré principal, nous proposons un algorithme qui construit la géodésique horizontale entre deux courbes et qui fournit un matching optimal. Dans un deuxième temps, nous introduisons une discrétisation de notre modèle qui est elle-même une structure riemannienne sur la variété de dimension finie M n+1 des "courbes discrètes" définies par n + 1 points, où M est de courbure sectionnelle constante. Nous montrons la convergence du modèle discret vers le modèle continu, et nous étudions la géométrie induite. Des résultats de simulations dans la sphère, le plan et le demi-plan hyperbolique sont donnés. Enfin, nous donnons le contexte mathématique nécessaire à l'application de l'étude de formes dans une variété au traitement statistique du signal radar, où des signaux radars localement stationnaires sont représentés par des courbes dans le polydisque de Poincaré via la géométrie de l'information.

Mots-clés : Étude de formes, variété riemannienne, matching entre courbes, géométrie de l'information.

Abstract : Probability on the spaces of curves and the associated metric spaces via information geometry; radar applications

We are concerned with the comparison of the shapes of open smooth curves that take their values in a Riemannian manifold M . To this end, we introduce a reparameterization invariant Riemannian metric on the infinite-dimensional manifold of these curves, modeled by smooth immersions in M . We derive the geodesic equation and solve the boundary value problem using geodesic shooting. The quotient structure induced by the action of the reparametrization group on the space of curves is studied. Using a canonical decomposition of a path in a principal bundle, we propose an algorithm that computes the horizontal geodesic between two curves and yields an optimal matching. In a second step, restricting to base manifolds of constant sectional curvature, we introduce a detailed discretization of the Riemannian structure on the space of smooth curves, which is itself a Riemannian metric on the finite-dimensional manifold M n+1 of "discrete curves" given by n + 1 points. We show the convergence of the discrete model to the continuous model, and study the induced geometry. We show results of simulations in the sphere, the plane, and the hyperbolic halfplane. Finally, we give the necessary framework to apply shape analysis of manifold-valued curves to radar signal processing, where locally stationary radar signals are represented by curves in the Poincaré polydisk using information geometry.

Keywords : Shape analysis, Riemannian manifold, optimal matching between curves, information geometry.

c 1 : 2 CHAPTERFigure 1 . 1 :

 1211 Figure 1.1: Two polygonal curves with small Hausdorff distance.

c 1

 1 (ϕ(t)) -c 2 (ψ(t)) .

Figure 1 . 2 :

 12 Figure 1.2: Optimal deformation between two shapes.

Figure 1 . 3 :

 13 Figure 1.3: Infinitesimal deformation of a curve and its decomposition into tangential and normal parts.

Figure 1 . 4 :

 14 Figure 1.4: Optimal deformations for a reparameterization invariant metric.

Definition 1 . 1 .

 11 We call elastic the two parameter-family of Riemannian metrics defined on M for any a, b ∈ R and c ∈ M, w, z ∈ T c M by

Algorithm 1 . 2 (

 12 Geodesic shooting). Let c 0 , c 1 ∈ M. Set w = log L2 c 0 (c 1 ) and repeat until convergence : 1. compute the geodesic s → c(s) = exp M c 0 (s w) using Algorithm 1.1, 2. evaluate the difference j := log L 2 c(1) (c 1 ) between the target curve c 1 and the extremity c(1) of the obtained geodesic, 3. compute the initial derivative ∇ s J(0) of the Jacobi field s → J(s) along c verifying J(0) = 0 and J(1) = j, 4. correct the shooting direction w = w + ∇ s J(0).

Definition 1 . 5 .

 15 The shape of a curve c 0 ∈ M is the equivalence class c 0 := {c ∈ M, ∃ϕ ∈ Diff + ([0, 1]), c = c 0 • ϕ} of all curves that can be obtained by reparameterizing c 0 . The set of shapes is the quotient space S = M/Diff + ([0, 1]).

Algorithm 1 . 3 (

 13 Optimal matching). Let c 0 , c 1 ∈ M. Set ĉ1 = c 1 and repeat until convergence: 1. Construct the geodesic s → c(s) linking c 0 to ĉ1 using Algorithm 1.2. 2. Compute the horizontal part s → c hor (s) of c and set ĉ1 = c hor (1).

Theorem 1 . 1 (

 11 Convergence of the discrete model to the continuous model).

CHAPTER 1 .

 1 INTRODUCTIONcomponents (m k ) k verify m 0 = m 1 = 0 and the following recurrence relation

  denotes the Euclidean metric on R n , D w = w /|c | is the derivation of w according to arc length, D w T = D w, v v is the projection of D w on the unit speed vector field v = c /|c |, and D w N = D w -D w T . Recall that this metric belongs to the class of elastic metrics, defined for any weights a, b ∈ R + as

Figure 2 . 1 :

 21 Figure 2.1: Illustration of the distance between two curves c 0 and c 1 in the space of curves M.

1 0∇ 0 ∇ s ĉs , ĉa ds 1 0

 100 a ĉs , ĉs ds = -1

as shown schematically in Figure 3 . 3 .

 33 Let us briefly explain why. If c(a, •), a ∈ (-ε, ε), is a family of geodesics starting from the same point p 0 at different speeds w(a) ∈ T p 0 N , i.e. c(a, s) = exp N p 0 (sw(a)), and J(s) = c a (0, s), s ∈ [0, 1] measures the way that these geodesics spread out, then we have ∇ s J(0) = ∇ s c a (0, 0) = ∇ a c s (0, 0) = ∇ a | a=0 ∂s| s=0 exp N p 0 (sw(a)) = ∇ a w(0).

Algorithm 2 . 2 (Figure 2 . 2 :

 2222 Figure 2.2: Geodesic shooting in the space of curves M.

Figure 2 . 3 :

 23 Figure 2.3: Steps of the first iteration of the geodesic shooting algorithm for two geodesic curves of the upper half-plane.

Proposition 2 . 5 (

 25 Horizontal part of a vector). Let c ∈ M be a smooth immersion. Then h ∈ T c M is horizontal for the elastic metric G a,b if and only if

Algorithm 2 . 4 (

 24 Constructing horizontal geodesics). Input : c 0 , c 1 ∈ M. Set ĉ1 = c 1 and repeat until convergence: 1. construct the geodesic s → c(s) between c 0 and ĉ1 (e.g. using geodesic shooting), 2. compute the horizontal part s → c hor (s) of c and set ĉ1 = c hor (1).

56 CHAPTER 3

 563 . DISCRETIZATION With these notations, the Riemannian metric G of Chapter 2 can be written at point c ∈ M for two tangent vectors w, z ∈ T c M as G c (w, z) = w(0), z(0) + 1 0

3. 2 .Figure 3 . 1 :

 231 Figure 3.1: A piecewise geodesic curve.

(3. 9 ) 3 . 1 .

 931 Remark Notice that in the flat case our definition gives (D τ w) k = w k+1 -w k . In the non-flat case (K = ±1), when the discretization gets "thinner", i.e. n → ∞ and |τ k | → 0 while n|τ k | stays bounded for all 0 ≤ k ≤ n, we get (D τ w) k = n→∞ w k+1 -w k + o(1).

CHAPTER 3 .

 3 DISCRETIZATIONand |w| ∞ := sup s,t∈[START_REF](s) of the horizontal part of α(s) is chosen to be the point y (s) whose value is closest to kp. Now we can go back to our initial problem[END_REF][START_REF] Amari | Differential geometry in statistical inference[END_REF]

Lemma 3 . 3 .

 33 k . Denoting by γ k (s) the geodesic of M linking x k (s) to x k+1 (s) for all s ∈ [0, 1] and 0 ≤ k ≤ n -1, we have the following result. The covariant derivatives of the functions f

Proposition 3 . 5 (

 35 Discrete horizontal space). Let α ∈ M n+1 and h ∈ T α M n+1 . Then h ∈ Hor n α if and only if it verifies

Definition 3 . 3 .

 33 The horizontal part (α hor , c) of a path s → (α(s), c(s)) ∈ M n is defined by(α(s), c(s)) := ϕ(s) (α hor (s), c(s)), ∀s ∈ [0, 1],where ϕ(s) ∈ Inj + (n, N ) verifies for all s ∈[START_REF](s) of the horizontal part of α(s) is chosen to be the point y (s) whose value is closest to kp. Now we can go back to our initial problem[END_REF][START_REF] Amari | Differential geometry in statistical inference[END_REF]

Figure 3 . 2 :

 32 Figure 3.2: Geodesics between parameterized curves in R 2 obtained by geodesic shooting (in blue) compared to the exact geodesics (in red).

CHAPTER 3 .

 3 DISCRETIZATION (a) The hyperbolic half-plane (b) The plane (c) The sphere

Figure 3 . 3 :

 33 Figure 3.3: Geodesics between parameterized curves obtained by geodesic shooting (in blue) compared to the L 2 -geodesic (in green) in (a) negative (b) zero and (c) positive curvature. We can see that our metric has a tendency to "shrink" the curve as it deforms from one state to the other, in comparison to the L 2 -metric.

Figure 3 . 4 :Figure 3 . 5 : 1 :

 34351 Figure 3.4: Geodesics between parameterized curves in H 2 (blue) and corresponding horizontal geodesics (red). The parameterization of the curve on the left-hand side is fixed, and the horizontal geodesic redistributes the points of the curve on the right-hand side to yield an optimal matching.

Figure 3 . 6 :

 36 Figure 3.6: Geodesics between the shapes of two curves in H 2 (in black) obtained by osuperimposing the horizontal geodesics between several sets of parameterizations of the curves. These horizontal geodesics are parameterizations of the same shape geodesic.

  ) for some constant λ 5 . Injecting (3.20), (3.21) and (3.26) in Equation (3.18) we obtain

C e -z z dz = R 2 e 0 2π 0 e

 200 -(x 2 +y 2 ) dxdy = +∞ -r 2 rdrdθ = π.

( 4 . 8 )

 48 The last equation tells us that the reflection coefficients live in the complex unit disk D|µ k | < 1, k = 1, . . . , n -1,and we obtain a coordinate map of T + n in the product space R * + × D n-1

CHAPTER

Figure 4 . 3 :

 43 Figure 4.3: Computation of the mean curve (in black) for 2 sets of 11 curves representing the evolution of the second reflection coefficient of 11 simulated helicopter radar signatures in the hyperbolic half-plane. Each signature corresponds to a different rotor speed. The mean is computed using Algorithm 4.1.

  s, t) = cosh (|τ k (s)|t) , b k (s, t) = sinh (|τ k (s)|t) /|τ k (s)| if K = -1, a k (s, t) = 1, b k (s, t) = 1 if K = 0, a k (s, t) = cos (|τ k (s)|t) , b k (s, t) = sin (|τ k (s))|t) /|τ k (s)| if K = +1.Integrating this and injecting it in (A.2) gives∇ s (w k+1 ) = (∇ s w k+1 ) + R (Y k , τ k ) (w k+1 ), (A.3)

Remark 3 . 3 |R 0 - 0 •

 3300 Let [0, 1] s → c(s, •) ∈ M be a C 1 path of smooth curves and [0, 1] s → α(s) ∈ M n+1 the discretization of size n of c. We denote as usual by q := c t /|c t | 1/2 and (q k ) k their respective SRV representations. When n → ∞ and|τ k | → 0 while n|τ k | stays bounded for all 0 ≤ k ≤ n, since ∇ s τ k = (D τ α ) k = (x k+1 -x k ) T + b -1 k (x k+1 -x k ) N and b -1 k → 1. Finally, we can see that |r 0 (s) -r 0 (s)| ≤ 1 n . . . • f (-) (R +1 ) -P +1 n ,0 c (R +1 )goes to 0 when n → ∞. We can show in a similar way that |g(-) k (r k ) -|q k |(r k + r k T )| → 0 when n → ∞.Proposition 3.3 (Discrete exponential map) Let [0, 1]

1 :

 1 = Id, n-1 =n := 0 and with the maps

  .18) Let us first consider the difference |c t -ĉ t |. Since ĉt (t) t, k n = ĉt ( k n ), we can write

  1 also. Therefore, for n large enough we can see that |ĉ s | ≤ λ 4 |c s | ∞ for some constant λ 4 . Injecting this into (3.23) gives

  4. RADAR APPLICATIONS

	4.2. HIGH RESOLUTION SPECTRAL ESTIMATION				97
	0.952	0.958
					mu389.5
					mu389.6
	0.9515	0.957		mu389.7
					mu389.8
					mu389.9
		0.956		mu390.0
	0.951				mu390.1
					mu390.2
		0.955		mu390.3 mu390.4
	0.9505				mu390.5
		0.954		mu m
	0.95
		0.953
	0.9495
		0.952 mu389.5
		mu389.6
	0.949	mu389.7
		0.951 mu389.8
		mu389.9
	0.9485	0.95 mu390.0 mu390.1
		mu390.2
		mu390.3
	0.948	mu390.4 0.949 mu390.5
		mu m
	0.9475	0.948
	-6	-8 -4		-6 -2		-4	0	-2	2	0	4		2	6	4	6
														×10	-3	×10 -3
	0.956	0.96
																mu389.5
																mu389.6
		0.958														mu389.7 mu389.8
	0.954															mu389.9
																mu390.0
																mu390.1
		0.956														mu390.2
																mu390.3
	0.952															mu390.4 mu390.5
		0.954														mu m
	0.95	0.952
		0.95											mu389.5 mu389.6
	0.948												mu389.7
													mu389.8
													mu389.9
		0.948											mu390.0
													mu390.1
	0.946												mu390.2 mu390.3
		0.946											mu390.4
													mu390.5
													mu m
	0.944	0.944
	-4	-3	-6	-2	-4 -1	0	-2	1	0	2	2	3	4	4	5	6	8
														×10	-3	×10	-3
	Figure 4.2: Computation of the mean curve (in black) for 2 sets of 11 curves
	representing the evolution of the second reflection coefficient of 11 simulated heli-
	copter radar signatures in the hyperbolic half-plane. Each signature corresponds to
	a different rotor speed. The mean is computed using Algorithm 4.1.

Remerciements APPENDIX B. GEOMETRY OF THE HYPERBOLIC PLANE

Bibliography

Chapter 3

Discretization

Abstract

In this chapter we assume that the base manifold M has constant sectional curvature K. We introduce a detailed discretization of the Riemannian structure presented in Chapter 2, which is itself a Riemannian structure on the product manifold M n+1 of "discrete curves" given by n + 1 points. We show that the discrete energy of a discretization of size n of a path of smooth curves converges to the continuous energy as n → ∞. We derive the geodesic equation, describe the discrete exponential map, the discrete Jacobi fields, and consider a discrete analog of the quotient structure. This enables us to show simulations in the hyperbolic half-plane (K = -1), the plane (K = 0) and the sphere (K = 1), solving in each case the boundary-value problem in the space of parameterized curves as well as the space of unparameterized curves. This chapter corresponds to most of the work of [START_REF] Le | A discrete framework to find the optimal matching between manifold-valued curves[END_REF].

Summary of the continuous model

The space of smooth parameterized curves

In Chapter 2, we considered smooth parameterized immersions c : [0, 1] → M , which we each represented by the pair formed by their starting point x = c(0) and renormalized speed vector field q = c / |c |. The renormalization in q allowed us to define a reparameterization invariant metric by pullback. Here we give an equivalent definition, for which the formulation of a discrete analog is more straightforward. For any tangent vector w ∈ T c M, consider a path of curves s → c w (s) ∈ M such that c w (0) = c and c w s (0) = ∂c w /∂s(0) = w. We denote by q w = c w t /|c w t | 1/2 the square root velocity representation of c w .

CHAPTER 3. DISCRETIZATION

Proof. Let h ∈ T α M be a tangent vector. It is horizontal if and only if it is orthogonal to any vertical vector, that is any vector of the form mv with m = (m k ) k ∈ R n+1 such that m 0 = m n = 0. Recall that by definition (D τ w) k := (w k+1 -w k ) T + b -1 k (w k+1 -a k w k ) N , and so with the notation λ k := v k+1 , v k , we get

The scalar product between h and mv is then

Changing the indices in the first sum and taking into account that m 0 = m n = 0, we obtain

Since this is true for all such m the summand is equal to zero for all k and we get the desired equation. The decomposition of a tangent vector w into a vertical part mv and a horizontal part w -

and so

We can decompose

To bound |ĉ s | we apply Lemma 3.1 to the Jacobi field J : [0, 1] u → ĉs (s, k+u n ) along the geodesic γ(u) = ĉ(s, k+u n ), that is

where, since γ (0) = 1 n ĉt (s, k n ) = τ k (s), the coefficients are defined by

This gives J(1)

Injecting this into (3.24), we obtain since u = ntk and ĉs (s

and so the metric on R * + × D n-1 is given by

that is

(4.9)

We obtain a product metric, in which we recognize a metric on R + and, up to a constant, the metric of the Poincaré disk D

Induced with metric (4.9), the coordinate space R * + × D n-1 becomes the product manifold R + × D n-1 = (R + , ds 2 0 ) × (D, ds 2 1 ) × . . . × (D, ds 2 n-1 ), that we call the Poincaré polydisk. It is important to note that while metric (4.2) is the Fisher metric on the space H + n of hermitian positive definite matrices of size n, the product metric (4.9) on the space of reflection coefficients that is used in practice for radar applications [START_REF] Barbaresco | Interactions between symmetric cone and information geometries: Bruhat-tits and siegel spaces models for high resolution autoregressive doppler imagery[END_REF], [START_REF] Arnaudon | Riemannian medians and means with applications to radar signal processing[END_REF], is not the Fisher information metric on the submanifold of stationary, centered multivariate Gaussian densities, as pointed out by Yang in [START_REF] Yang | Medians of probability measures in Riemannian manifolds and applications to radar target detection[END_REF]. It can rather be seen as the restriction of the Fisher metric to that submanifold.

High resolution spectral estimation of locally stationary radar signals

Now let us give some perspectives on the utility of shape analysis of open curves in a manifold to radar signal processing. We have seen that a stationary centered Gaussian process can be represented by a point in the Poincaré polydisk R + × D n-1 . As mentioned in the introduction, there are situations where the stationarity hypothesis is not appropriate, such as in the presence of very inhomogeneous clutter or when the target moves during the time interval of the burst. In our model, we assume that the underlying process of

HIGH RESOLUTION SPECTRAL ESTIMATION

93 a radar signature is locally stationary, and we are interested in its temporal modulations. If we represent each stationary portion by an element of the polydisk R + × D n-1 , then the locally stationary process is represented by a series of points, or "discrete curve", in the polydisk. Another way to see it, is that each point in the Poincaré polydisk represents the high-resolution autoregressive spectrum of a stationary portion of the radar signal. Estimating time-series of high-resolution spectra, or spectrograms, is at the heart of Micro-Doppler analysis, and studying their evolution in time gives valuable information on the temporal modulations of the signal. Here we apply the Riemannian framework detailed in the previous chapters to the spectral analysis of locally stationary helicopter signatures. The data we use for this example is synthetic data generated by a sim-94 CHAPTER 4. RADAR APPLICATIONS ulator of helicopter signatures. Using this simulator, we obtain a signature z = (z 1 , . . . , z N ) ∈ C N that simulates the reflected signal received by a fixed radar antenna after sending a burst of N pulses in the direction of a fixed helicopter. We assume that the underlying process is centered, locally stationary and Gaussian, and for each position of a gliding window of size n < N (to be adjusted), we estimate a high resolution autoregressive spectrum for the centered stationary Gaussian process Z (i) corresponding to the selected portion of observation z (i) = (z i , . . . , z i+n-1 ),

S

(n-1)

In Figure 4.1, we show the difference of resolution between the FFT spectrum and the Capon spectrum, which is obtained from the autoregressive spectra as follows

S

(n-1)

The coefficients (a

n-1 ) of each autoregressive spectrum (corresponding to a stationary portion of the signal) can be computed from the data z (i) using the regularized Burg algorithm (see e.g. [START_REF] Arnaudon | Riemannian medians and means with applications to radar signal processing[END_REF]), which simultaneously gives the reflection coefficients (P 0 , µ 1 , . . . , µ n-1 ) -as shown by the Levinson recursion, the computation of the former is done through the calculation of the latter. This gives us a time series of Nn + 1 spectra, as well as the corresponding time series of reflection coefficients in the polydisk

, which correspond to observations of the "real" evolution (P 0 (t), µ 1 (t), . . . , µ n-1 (t)) of the locally stationary signal. With this choice of representation, comparing two vectors of radar observations -or the associated spectrograms -can be carried out by computing the distance between the two corresponding curves in the product manifold R * + × D, which, thanks to the product metric, is the same as comparing their components separately -that is, pairwise comparing the evolutions of each reflection coefficient µ k (t) in the Poincaré disk. The evolution in the Poincaré disk of the first four reflection coefficients of a signature of size N = 40 with stationary portions of size n = 8 is shown as an example in Figure 4.1d.

Representing a radar signature as an element of the infinite dimensional Riemannian manifold of curves in R * + × D and applying the Riemannian framework presented in the previous chapters allows us to do statistics on these objects, such as defining the mean, median and variance of a set of Appendix A

Proofs of Chapter 3

Lemma 3.1 Let γ : [0, 1] → M be a geodesic of a manifold M of constant sectional curvature K, and J a Jacobi field along γ. Then the parallel transport of J(t) along γ from γ(t) to γ(0) is given by Proof. As a Jacobi field along γ, J satisfies the well-known equation

If M is flat, we get ∇ 2 t J(t) = 0 and so J(t) = J(0) + t∇ t J(0). If not, we can decompose J in the sum J = J T + J N of two vector fields that parallel translate along γ, by projecting it in the basis

t J(t), γ (t) = 0 and γ is parallel along γ, we get by integrating twice that

the normal component J N is also a Jacobi field, that is it verifies

And since M has constant sectional curvature K, for any vector field w along γ we have

and the differential equation verified by J N can be rewritten

Since the speed of the geodesic γ has constant norm, the solution to that differential equation is of the form

where the functions ãk (t)'s and bk (t)'s depend on the value of K as defined in the lemma.

if K is the sectional curvature of the base manifold.

Proof. Fix 0 ≤ k ≤ n and let w k+1 : s → w k+1 (s) be a vector field along the curve x k+1 : s → x k+1 (s). By definition,

Consider the path of gedesics s → γ k (s, •) such that for all s ∈ [0, 1], γ k (s, 0) = x k (s), γ k (s, 1) = x k+1 (s) and t → γ k (s, t) is a geodesic. We denote APPENDIX A. PROOFS OF CHAPTER 3 and injecting this in (A.1) finally gives,

which is what we wanted. The covariant derivative ∇ s g (-)

k (w k+1 ) can be computed in a similar way. Proposition 3.2 (Discrete geodesic equations) A path in M n+1 , s → α(s) = (x 0 (s), . . . , x n (s)), is a geodesic for metric G n if and only if its SRV representation s → x 0 (s), (q k (s)) k verifies the following differential equations

for all k = 0, . . . , n -1, with the notations (3.7) and R k := R(q k , ∇ s q k )x k .

Proof. We consider a variation (-δ, δ) a → α(a, •) = (x 0 (a, •), . . . , x n (a, •)) of this curve which coincides with α for a = 0, i.e. α(0, s) = α(s) for all s ∈ [0, 1], and which preserves the end points of α, i.e. α(a, 0) = α(0) and α(a, 1) = α(1) for all a ∈ (-δ, δ). The energy of this variation with respect to metric G n can be seen as a real function of the variable a and is given by

and its derivative (E n ) (a) with respect to a is

the coefficients of the discrete geodesic equation (3.14) for α converge to the coefficients of the continuous geodesic equation (3.4) for c, i.e.

∇ s x 0 (s) = -r 0 (s) + o(1),

for all s ∈ [0, 1] and k = 0, . . . , n-1, where r n-1 = 0 and for k = 1, . . . , n-2,

with the exception that the sum starts at = 0 for r 0 .

Proof. This is due to three arguments : (1) at the limit,

parallel transport along a piecewise geodesic curve uniformly converges to the parallel transport along the limit curve, and (3) the discrete curvature term R k (s) converges to the continuous curvature term R(q, ∇ s q)c s (s, k n ) for all k. Indeed, let ĉ be the unique piecewise geodesic curve of which α is the discretization, i.e. c k n = ĉ k n = x k for all k = 0, . . . , n and ĉ is a geodesic on each segment k n , k+1 n . Defining

the geodesic equations can be written in terms of the vectors rk

We can show that for any 0 ≤ k ≤ ≤ n -2 and any vector

Since |a j -1|/|τ k | 2 → 0 when n → ∞ and n|τ k | stays bounded, we have for all 0 ≤ j ≤ n and n large enough |a j -1| ≤ 1 n 2 , and using the fact that parallel transport along a piecewise geodesic curve uniformly converges to the parallel transport along the limit curve, we get

when n → ∞. Now, denoting by R(s, t) := R(q, ∇ s q)c s (s, t)

the curvature term involved in the continuous geodesic equations, we have since

and since the portion of c(s, •) on the segment [ k n , k+1 n ] is close to a geodesic at the limit, |nτ k -c t ( k n )| → 0 when n → ∞, and so does |q k (s) -q( k n )|. Similarly,

The last term can be bounded, for n large enough, by

gives

Appendix B

Geometry of the hyperbolic half-plane

Here we give a few tools which are necessary to work in the hyperbolic halfplane representation. Along with the Poincaré disk, the Klein model and others, the hyperbolic half-plane H = {z = x + iy ∈ C, y > 0} is one of the representations of two-dimensional hyperbolic geometry. The Riemannian metric is given by

This means that the scalar product between two tangent vectors u = u 1 +iu 2 and v = v 1 + iv 2 at a point z = x + iy is

Using the usual formula (see e.g. [START_REF] Perdigao | Riemannian geometry[END_REF]) to compute the Christoffel symbols, we can easily compute the covariant derivative of a vector field

where

Let us now remind a well-known expression [START_REF] Perdigao | Riemannian geometry[END_REF] for the Riemann curvature tensor in a manifold of constant sectional curvature. Recall that H has constant sectional curvature K = -1. can be written R(X, Y )Z = K ( Y, Z X -X, Z Y ) .

For the algorithms described above, we need to be able to compute the geodesic starting from a point p ∈ H at speed u 0 ∈ T p H -in other words, the exponential map u → exp H p (u) -as well as the geodesic linking two points p and q, with the associated initial vector speed -the inverse q → log H p (q) of the exponential map. The geodesics of the hyperbolic half-plane are vertical segments and half-circles whose origins are on the x-axis, as shown in Figure B.1, and they can be obtained as images of the vertical geodesic following the y-axis by a Moebius transformation z → az+b cz+d , with adbc = 1. Proposition B.2 (Geodesics of H and logarithm map). Let z 0 = x 0 + iy 0 and z 1 = x 1 + iy 1 be two elements of H.

• If x 0 = x 1 , then the geodesic going from z 0 to z 1 is the segment γ(t) = iy(t) with y(t) = y 0 e t ln y 1 y 0 , and the logarithm map is given by log H z 0 (z 1 ) = iy 0 ln y 1 y 0 .

• If x 0 = x 1 , the geodesic is given by γ(t) = x(t) + iy(t) with

where the coefficients of the Moebius transformation can be deduced from the center x Ω and the radius R of the semi-circle going through The logarithm map is in turn given by

Proof. The geodesic γ(t) = x(t) + iy(t) linking two points vertically aligned z 0 = x 0 +iy 0 and z 1 = x 0 +iy 1 is a vertical segment γ(t) = iy(t). It verifies the geodesic equation ∇ γ(t) γ(t) = 0. Using the expression B.1 of the covariant derivative of a vector field in H, this gives the equation ÿy = ẏ2 , which can be rewritten as ÿ ẏ = ẏ y . Integrating twice, we find that y(t) = y 0 e t ln y 1 y 0 . Now if z 1 = x 1 + iy 1 with x 1 = x 0 , the geodesic γ is the image by a Moebius transformation z → az+b cz+d (with adbc = 1) of a vertical line γ(t) = iȳ(t), which gives

We know that γ describes a semi-circle Ω whose origin x Ω is on the x-axis, and that one end of the vertical line γ is sent to the point a/c and the other to the point b/d. This implies that the center of the semi-circle is half-way between the two x Ω = ad+bc 2cd , and that the radius is R = 1 2cd . These two equations as well as the condition adbc = 1 gives a system of equations for the coefficients a, b, c and d, which, if we choose to set d = 1, yields the desired expressions. If the extremity z 0 is sent by the inverse of the obtained Moebius transformation on iȳ 0 , and z 1 on iȳ 1 , then the segment corresponding to the portion of γ linking z 0 to z 1 is γ(t) = ȳ0 e t ln ȳ1 ȳ0 . Taking the derivative of B.2 in t = 0 gives the logarithm map.

We now give the exponential map in H.

Proposition B.3 (Exponential map in H)

. Let z 0 = x 0 + iy 0 be an element of H and u 0 = ẋ0 + i ẏ0 a tangent vector. Then the exponential map is given by exp H z 0 (u 0 ) = γ(1), where Proof. The exponential map uses the same equations as the logarithm map with the difference that u 0 is known instead of z 1 . The proof is very similar to the the proof of Proposition B.2 and is not detailed here.

Finally, we give the expression of parallel transport along a geodesic in the hyperbolic plane. Proposition B.4 (Parallel transport in H). Let t → γ(t) be a curve in H with coordinates x(t), y(t), and u 0 ∈ T γ(t 0 ) H a tangent vector. The parallel transport of u 0 along γ from t 0 to t is given by u(t) = y(t) y(t 0 ) cos θ(t 0 , t) sin θ(t 0 , t) -sin θ(t 0 , t) cos θ(t 0 , t) u 0 , where θ(t i , t f ) =

y(τ ) dτ . If γ is a vertical segment then θ(t i , t f ) = 0, and if it is a portion of a circle, we get θ(t i , t f ) = 2 (arg(d + icȳ(t f )) -arg(d + icȳ(t i ))) , where the coefficients c and d of the Moebius transformation can be computed as explained previously, and γ = iȳ is the pre-image of γ by that transformation. Now that we have these explicit formulas at our disposal, we are able to test the algorithms described above in the simple case where the base manifold M has constant sectional curvature K = -1. Note that computations are further simplified by the existence of a global chart.

Proof. Parallel transporting a vector u 0 ∈ T γ(t 0 ) H along a curve γ from t 0 to t is gives a vector field u satisfying ∇ γ(t) u = 0 and u(t 0 ) = u 0 . Using Equation B.1, this can be rewritten u = Au where

A is of the form aI + bK where I is the identity matrix and K = 0 1 -1 0 and so the solution is u(t) = u(t 0 ) exp The matrix B(t) is diagonalizable and therefore its exponential can be written exp B(t) = e a(t 0 ,t) cos θ(t 0 , t) sin θ(t 0 , t)

sin θ(t 0 , t) cos θ(t 0 , t) ,

where a(t 0 , t) = ln y(t) y(t 0 ) and θ(t 0 , t) = t t 0 ẋ(τ ) y(τ ) dτ . This gives us the desired formula.