domain: stat math

text:
met, Lennard, Mathias Jon Nicholas

Elodie, Jairo, Zhijie, Tom Luiś

Marcus

Libby

Thomas

Christian, Michael, Gefan Sofia Lily

Olivier Liwei

Thank You Alain

Pierre Klas

Keywords: geometric statistics, differential geometry, Riemannian geometry, sub-Riemannian geometry, mathematical statistics, machine learning

In this thesis, we propose new methods for dimension reduction based on differential geometry, that is, finding a representation of a set of observations in a space of lower dimension than the original data space. Methods for dimension reduction form a cornerstone of statistics, and thus have a very wide range of applications. For instance, a lower dimensional representation of a data set allows visualization and is often necessary for subsequent statistical analyses. In ordinary Euclidean statistics, the data belong to a vector space and the lower dimensional space might be a linear subspace or a non-linear submanifold approximating the observations. The study of such smooth manifolds, differential geometry, naturally plays an important role in this last case, or when the data space is itself a known manifold. Methods for analysing this type of data form the field of geometric statistics. In this setting, the approximating space found by dimension reduction is naturally a submanifold of the given manifold. The starting point of this thesis is geometric statistics for observations belonging to a known Riemannian manifold, but parts of our work form a contribution even in the case of data belonging to Euclidean space, R d .

An important example of manifold valued data is shapes, in our case discrete curves or surfaces. In evolutionary biology, researchers are interested in studying reasons for and implications of morphological differences between species. Shape is one way to formalize morphology. This application motivates the first main contribution of the thesis. We generalize a dimension reduction method used in evolutionary biology, phylogenetic principal component analysis (P-PCA), to work for data on a Riemannian manifold -so that it can be applied to shape data. P-PCA is a version of PCA for observations that are assumed to be leaf nodes of a phylogenetic tree. From a statistical point of view, the important property of such data is that the observations (leaf node values) are not necessarily independent. We define and estimate intrinsic weighted means and covariances on a manifold which takes the dependency of the observations into account. We then define phylogenetic PCA on a manifold to be the eigendecomposition of the weighted covariance in the tangent space of the weighted mean. We show that the mean estimator that is currently used in evolutionary biology for studying morphology corresponds to taking only a single step of our Riemannian gradient descent algorithm for the intrinsic mean, when the observations are represented in Kendall's shape space.

Our second main contribution is a non-parametric method for dimension reduction that can be used for approximating a set of observations based on a very flexible class of submanifolds. This method is novel even in the case of Euclidean data. The method works by constructing a subbundle of the tangent bundle on the data manifold M via local PCA. We call this subbundle the principal subbundle. We then observe that this subbundle induces a sub-Riemannian structure on M and we show that the resulting sub-Riemannian geodesics with respect to this structure stay close to the set of observations. Moreover, we show that sub-Riemannian geodesics starting from a given point locally generate a submanifold which is radially aligned with the estimated subbundle, even for non-integrable subbundles. Non-integrability is likely to occur when the subbundle is estimated from noisy data, and our method demonstrates that sub-Riemannian geometry is a natural framework for dealing which such problems. Numerical experiments illustrate the power of our framework by showing that we can achieve impressively large range reconstructions even in the presence of quite high levels of noise.

Méthodes riemanniennes et sous-riemanniennes pour la réduction de dimension

Nous proposons dans cette thèse de nouvelles méthodes de réduction de dimension fondées sur la géométrie différentielle. Il s'agit de trouver une représentation d'un ensemble d'observations dans un espace de dimension inférieure à l'espace d'origine des données. Les méthodes de réduction de dimension constituent la pierre angulaire des statistiques et ont donc un très large éventail d'applications. Dans les statistiques euclidiennes ordinaires, les données appartiennent à un espace vectoriel et l'espace de dimension inférieure peut être un sous-espace linéaire ou une sous-variété non linéaire approximant les observations. L'étude de telles variétés lisses, la géométrie différentielle, joue naturellement un rôle important dans ce dernier cas. Lorsque l'espace des données est lui-même une variété, l'espace approximant de dimension réduite est naturellement une sous-variété de la variété intiale. Les méthodes d'analyse de ce type de données relèvent du domaine des statistiques géométriques. Les statistiques géométriques pour des observations appartenant à une variété riemannienne sont le point de départ de cette thèse, mais une partie de notre travail apporte une contribution même dans le cas de données appartenant à l'espace euclidien, R d .

Les formes, dans notre cas des courbes ou des surfaces discrètes, sont un exemple important de données à valeurs dans les variétés. En biologie évolutive, les chercheurs s'intéressent aux raisons et aux implications des différences morphologiques entre les espèces. Cette application motive la première contribution principale de la thèse. Nous généralisons une méthode de réduction de dimension utilisée en biologie évolutive, l'analyse en composantes principales phylogénétiques (P-PCA), pour travailler sur des données à valeur dans une variété riemannienne -afin qu'elle puisse être appliquée à des données de forme. P-PCA est une version de PCA pour des observations qui sont les feuilles d'un arbre phylogénétique. D'un point de vue statistique, la propriété importante de ces données est que les observations ne sont pas indépendantes. Nous définissons et estimons des moyennes et des covariances intrinsèquement pondérées sur une variété qui prennent en compte cette dépendance des observations. Nous définissons ensuite l'ACP phylogénétique sur une variété comme la décomposition propre de la covariance pondérée dans l'espace tangent de la moyenne pondérée. Nous montrons que l'estimateur de moyenne actuellement utilisé en biologie évolutive pour étudier la morphologie correspond à ne prendre qu'une seule étape de notre algorithme de descente de gradient riemannien pour la moyenne intrinsèque, lorsque les observations sont représentées dans l'espace des formes de Kendall. Notre deuxième contribution principale est une méthode non paramétrique de réduction de dimension fondée sur une classe très flexible de sous-variétés qui est novatrice même dans le cas de données euclidiennes. Grâce à une PCA locale, nous construisons tout d'abord un sous-fibré du fibré tangent sur la variété des données que nous appelons le sous-fibré principal. Cette distribution (au sens géométrique) induit une structure sous riemannienne. Nous montrons que les géodésiques sousriemanniennes correspondantes restent proches de l'ensemble des observations et que l'ensemble des géodésiques partant d'un point donné génèrent localement une sousvariété qui est radialement alignée avec le sous-fibré principal, même lorsqu'il est non intégrables, ce qui apparait lorsque les données sont bruitées. Notre méthode démontre que la géométrie sous-riemannienne est le cadre naturel pour traiter de tels problèmes. Des expériences numériques illustrent la puissance de notre cadre en montrant que nous pouvons réaliser des reconstructions d'une extension importante, même en présence de niveaux de bruit assez élevés. Et vigtigt eksempel på data på en mangfoldighed er former, i vores tilfaelde diskrete kurver eller overflader. I evolutionsbiologi er forskere interesseret i at studere grunde til og implikationer af morfologiske forskelle mellem arter. Former er én måde at formalisere morfologi på. Denne anvendelse motiverer det første hovedbidrag i denne afhandling. We generaliserer en metode til dimensionsreduktion brugt i evolutionsbiologi, phylogenetisk principal component analysis (P-PCA), til at virke for data på en Riemannsk mangfoldighed -så den kan anvendes til observationer af former. P-PCA er en version af PCA for observationer som antages at vaere de yderste knuder i et phylogenetisk trae. Fra et statistisk synspunkt er den vigtige egenskab ved sådanne observationer at de ikke nødvendigvis er uafhaengige. We definerer og estimerer intrinsiske vaegtede middelvaerdier og kovarianser på en mangfoldighed, som tager højde for sådanne observationers afhaengighed. Vi definerer derefter phylogenetisk PCA på en mangfoldighed som egendekomposition af den vaegtede kovarians i tanget-rummet til den vaegtede middelvaerdi. Vi viser at estimatoren af middelvaerdien som pt. bruges i evolutionsbiologi til at studere morfologi svarer til at tage kun et enkelt skridt af vores Riemannske gradient descent algoritme for den intrinsiske middelvaerdi, når formerne repraesenteres i Kendall´s form-mangfoldighed.

Mots

Vores andet hovedbidrag er en ikke-parametrisk metode til dimensionsreduktion som kan bruges til at approksimere et data saet baseret på en meget flexibel klasse af submangfoldigheder. Denne metode er ny også i tilfaeldet med Euklidisk data. Metoden virker ved at konstruere et under-bundt af tangentbundet på datamangfoldigheden M via lokale PCA´er. Vi kalder dette underbundt principal underbundtet. Vi observerer at dette underbundt inducerer en sub-Riemannsk struktur på M og vi viser at sub-Riemannske geodaeter fra et givent punkt lokalt genererer en submangfoldighed som radialt flugter med det estimerede subbundt, selv for ikke-integrable subbundter. Ved støjfyldt data forekommer ikke-integrabilitet med stor sandsynlighed, og vores metode demonstrerer at sub-Riemannsk geometri er en naturlig tilgang til at håndtere dette. Numeriske eksperimenter illustrerer styrkerne ved metoden ved at vise at den opnår rekonstruktioner over store afstande, selv under høje niveauer af støj.

Nøgleord: geometrisk statistik, differentialgeometry, Riemannsk geometri, sub-Riemannsk geometri, matematisk statistik, maskinlaering

Introduction

Statistics on manifolds

Let x 1 , . . . , x N be a set of observations. The field of statistics traditionally assumes that such observations belong to a (finite dimensional) vector space.

That is, the methods rely on being able to add observations, x i + x j , and scale them,

x i • c, c ∈ R ≥0 .
However, this assumption of a vector space structure is prohibitive; important phenomena are naturally modelled as belonging to nonvector spaces. This thesis provides new methods for analyzing such non-linear data.

A simple example of non-linear data is observations on a surface in R 3 , such as the sphere. Such observations can be given to us in the form of vectors in

R 3 , i.e. x i = [(x i ) 1 , (x i ) 2 , (x i ) 3
], but prior knowledge might imply that each observation satisfies the equation

∥x i ∥ = (x i ) 2 1 + (x i ) 2 2 + (x i ) 2 3 = 1, i ∈ {1, . . . , N }, (1.1.1)
meaning that they are points on the unit sphere embedded in R 3 . Such data on the sphere, and its analogues in different dimensions, is studied in the field of directional statistics ([K. V. [START_REF] Mardia | Directional statistics[END_REF], [START_REF] Pewsey | Recent advances in directional statistics[END_REF]). The observations can be literal directions, e.g. wind directions, they can be observations of time points, or more generally any type of Euclidean observations that has been normalized, meaning that only directions matter, not magnitude (see [START_REF] Banerjee | Clustering on the Unit Hypersphere using von Mises-Fisher Distributions[END_REF] for applications to text classification, gene expression analysis and other domains). Applying traditional Euclidean statistical methods to such data, i.e. treating them as points in R 3 , disregards this prior knowledge and thus the output might not be sensible. For example, computing the Euclidean mean of the observations on the sphere yields a point that does not satisfy the constraint (1.1.1) -it is not a point on the sphere. There is a need for statistical methods which can incorporate prior knowledge about the space in which a set of observations might reside. The particular type of non-linear data spaces studied in this thesis are manifolds -of which the sphere is an example. In particular, we assume that the space containing the observations is a Riemannian manifold, implying 1. Introduction that e.g. distances between points can be computed. This type of statistics is called Riemannian geometric statistics. Given data on a Riemannian manifold, the field comprises methods for computing mean values, dimension reduction, regression, classification etc. (see [START_REF] Pennec | Riemannian geometric statistics in medical image analysis[END_REF]] and [S. F. Huckemann and Eltzner 2021] for an overview). Determining a set of constraints satisfied by the observations is a common way to arrive at a manifold data space in applications. Other examples of this are matrix manifolds such as the manifold of symmetric positive definite matrices,

SP D(d) = {A ∈ R d×d | A T = A, x ̸ = 0 =⇒ x T Ax > 0}.
Such matrices constitutes the observations in e.g. the study of diffusion tensor images, used in medical imaging, where they form covariance matrices of Brownian motions describing the motion (diffusion) of water molecules in tissue (see e.g. [START_REF] Lazar | Statistical analysis of diffusion tensors in diffusion-weighted magnetic resonance imaging data: Commentary[END_REF]. Another example of a matrix manifold is the rotation group SO(3), which is abundant in engineering applications (see [START_REF] Chirikjian | Engineering applications of noncommutative harmonic analysis: with emphasis on rotation and motion groups[END_REF]). In [START_REF] Hauberg | Principal curves on Riemannian manifolds[END_REF]] it is used to model motion capture data.

Another approach to deriving a manifold representing a data set is by introducing equivalence relations, that is, to consider certain subsets of points in the original space as identical. The result is a quotient space, which is often a manifold. An example of this is the notion of a shape space, which is a central motivation for geometric statistics in general as well as for the work in this thesis in particular. Consider the Euclidean space R d•k ∼ = R d×k , interpreted as a set of matrix-valued points, where each x ∈ R d×k is a collection of k points in R d (the columns), called landmarks. These k points constituting x could for example be points along a discretized curve in R 2 , outlining some 'shape'. What do we mean by 'shape'? One answer is that a shape is that which is left after removing the effects of translation, scaling and rotation. This can be formalized as a set of equivalence relations, where two points x, x ′ ∈ R d×k are identified if there exists a translation, a scaling and a rotation such that applying these operations to each landmark in x makes them equal to the corresponding landmarks in x ′ . Identifying such points in R d×k leads to Kendall's shape space, which is a quotient manifold [START_REF] Kendall | Shape manifolds, procrustean metrics, and complex projective spaces[END_REF]]. Kendall's shape space is part of the theoretical foundation of the field of geometric morphometrics -the study of morphology (shape) in biology [START_REF] Klingenberg | Walking on Kendall's shape space: understanding shape spaces and their coordinate systems[END_REF]]. Another example of a shape manifold is the LDDMM framework [START_REF] Younes | Shapes and diffeomorphisms[END_REF]], that enables analysis of full shapes, e.g. continuous curves and surfaces. Such spaces of full shapes has been used in the context of computational anatomy and medical imaging for studying organs like the heart or brain (see [START_REF] Guigui | Cardiac motion modeling with parallel transport and shape splines[END_REF]] for an analysis of heart movement).

Lastly, we note that defining a manifold in terms of a set of constraints or a set of equivalence relations can be complementary views on the same object. For example, it turns out that Kendall's shape space for k = 3 landmarks in R 2 , interpreted as a space of triangle shapes, can be identified with the sphere. Thus, performing statistics on such triangle shapes is equivalent to performing statistics on the sphere.

In the next chapter, we introduce the subject of geometric statistics in more detail. We put a special emphasis on methods for dimension reduction, which is the main subject of the thesis. Our contributions to this problem are described in the following subsection.

Thesis organization and contributions

Chapter 2 In this chapter, we present background theory concerning the parts of Riemannian geometric statistics most closely related to our contributions in Chapters ?? and ??. We start out by elaborating on the notion of doing statistics for observations belonging to non-vector spaces. We then summarize the notions in Riemannian geometry that are fundamental to Riemannian geometric statistics, and subsequently explain broadly how the various operations are used, and computed, for this type of statistics. We devote some space to the notion of tangent space approximations, and statistics based on these, since they are used extensively both a tool and a benchmark in the methods we develop in chapters ?? and ??. The exposition of tangent space statistics in this Chapter, and in Appendix A.2 where more technical details are given, includes some points which are of both theoretical and computational importance but which has, to the best of our knowledge, not been explicitly presented in the literature on geometric statistics. The section is concluded with presentations of the two foundational methods for dimension reduction on Riemannian manifolds; tangent PCA and Principal geodesic analysis (PGA). The methods for dimension reduction that we present in chapters ?? and ?? are novel extensions and combinations of these two methods. In Appendix A.3 we elaborate on the complications with computing PGA (which is an open problem), and rewrite the objective function to something more tractable. In Appendix A.4 we derive a Taylor expansion of the PGA objective function which reduces to tangent PCA when only the 1st order term is considered and includes curvature in the higher order terms.

Introduction

Chapter 3 This chapter presents an extension of Phylogenetic principal component analysis (p-PCA) to manifold valued data. P-PCA is a version of PCA for observations that are assumed to be leaf nodes of a phylogenetic tree. Mathematically, a phylogenetic tree is a rooted, bifurcating tree graph, which in evolutionary biology represents evolutionary relationships between species. From a statistical point of view, the important property of such data is that the observations (leaf node values) are not necessarily independent. In particular, they are the endpoints of Brownian motions that are more or less coupled depending on how closely related the species are. P-PCA consists of eigendecomposition of the socalled phylogenetic covariance matrix, which is centered around the socalled phylogenetic root, a type of mean. These notions of mean and covariance takes the dependency of the observations into account. In evolutionary biology there is a need for applying this method to data on manifolds, in particular shapes. We therefore extend to a general Riemannian manifold the probabilistic model of Brownian motions structured according to a tree, and define intrinsic estimators of the phylogenetic root and covariance. We then define phylogenetic PCA on a Riemannian manifold to be eigendecomposition of this covariance matrix in the tangent space at the mean. We show simulation results on the sphere (Kendall's shape space for triangles), demonstrating well-behaved error distributions and fast convergence of estimators. Furthermore, we apply the method to a data set of mammal jaws, represented as points in the LDDMM landmark manifold. Lastly, we discuss the relationship between our approach and the approach taken to the same problem in geometric morphometrics -the field of biology analysing shapes (morphology). We show that the phylogenetic root estimator used in geometric morphometrics corresponds to taking a single step of our Riemannian gradient descent algorithm for the intrinsic root, when the observations are represented in Kendall's shape space.

Chapter 4 The method presented in Chapter ?? is an extension of tangent PCA. Tangent PCA is based on geodesics of the manifold M to which the observations belong, which is a generalization of straight lines to such a manifold. In this chapter, we describe a method for dimension reduction that is based on a more flexible class of curves. Furthermore, whereas the previous chapter presented an extension of a method that is already well defined in Euclidean space, the method presented in this chapter is novel also in the case of M = R d . The method provides a solution to a number of dimension reduction problems, such as; construction of a lower-dimensional submanifold approximating the observations, representation of the observations in a lower-dimensional Euclidean space, and metric learning, in the sense of estimating a distance metric on M which reflects a lower dimensional geometry. The method works by constructing a subbundle of the tangent bundle on M via local tangent PCA's, which we call the principal subbundle. We then observe that this subbundle induces a sub-Riemannian (SR) structure on M , and we show that sub-Riemannian geodesics with respect to this structure follow the point cloud and can be used to provide solutions to the aforementioned problems. For example, we show that it is possible to generate submanifolds consisting of such SR geodesics, even if the subbundle is non-integrable, which loosely means that it doesn't 'naturally' determine submanifolds of M . Non-integrability is likely to occur when the subbundle is estimated from noisy data, and our method demonstrates that sub-Riemannian geometry is a natural framework for dealing which such noise. Last, but not least, we implemented a class with methods for sub-Riemannian geometry in geomstats [START_REF] Miolane | Geomstats: A Python Package for Riemannian Geometry in Machine Learning[END_REF], making computations in sub-Riemannian geometry available for future applications in geometric and Euclidean statistics.

Chapter 5 We round off by giving a brief summary of the contributions of the thesis and describing current and future directions of work related to phylogenetic PCA and principal subbundles.

Publications and invited talks

The material presented in Chapter ?? is an extension of the paper, Tangent phylogenetic PCA [Akhøj, Pennec, and Sommer 2023], published in the proceedings of the Scandinavian Conference on Image Analysis, 2023.

The material presented in Chapter ?? is based on the preprint Principal subbundles for dimension reduction [START_REF] Akhøj | Principal subbundles for dimension reduction[END_REF]

Generalizing Euclidean statistical methods

Many Euclidean statistical methods can be reformulated in terms of operations that does not depend on the vector space structure, suggesting that they make sense on more general spaces. Often, there exist multiple formulations of the same method that are equivalent in Euclidean space, but which lead to different solutions in the more general space. One approach to generalizing a given method is to reformulate it in terms of distances, making it well-defined on a general

metric space. Example 2.1.1 (The empirical mean in R d ). Let {x i } i=1...N be observations in R and let d(x, y) = |x -y| be the Euclidean distance on R. Consider the optimization problem μ = argmin µ∈R N i=1 d(µ, x i ) 2 . (2.1.1)
Differentiating with respect to µ gives

d dµ N i=1 |µ -x i | 2 = 2µN -2 N i=1 x i and d 2 d 2 µ N i=1 |µ -x i | 2 = 2N > 0. Thus, the unique minimizer of (2.1.1) is μ = 1 N N i=1 x i , (2.1.2)
the arithmetic mean, or empirical first moment. Likewise, for x 1 , . . . , x N ∈ R d one can similarly show that the solution is given by μ = 1 N N i=1 x i . In Equation (2.1.2), vector space operations addition and scalar multiplication are applied to the observations, while the minimization problem (2.1.1) only applies a distance metric to the observations. The above formulation of the arithmetic mean as a distance minimization problem can be used as a definition of a mean on a general metric space (M, d), where M is some set and d : M × M → R ≥0 is a distance metric on M . The point here is that M is not necessarily a vector space. Definition 2.1.2 (Fréchet mean). [START_REF] Fréchet | Les éléments aléatoires de nature quelconque dans un espace distancié[END_REF]] Let (M, d) be a metric space, and let x 1 , . . . , x N be points on M . The Fréchet mean is defined as a minimizer of the following minimization problem,

μ = argmin µ∈M N i=1 d(µ, x i ) 2 .
(2.1.3)

Background: geometric statistics

The Fréchet mean optimization problem is an instance of the following more general problem, for observations x 1 , . . . , x N ∈ M ,

U ∈ argmin U ∈Q N i=1 d(x i , π U (x i )) 2 , (2.1.4)
where Q is some family of closed subsets of M and π is a projection of a point

x ∈ M to a subset U ∈ Q with respect to the distance metric d;

π U (x) ∈ argmin p∈U d(x, p), (2.1.5)
assumed to be unique for the considered points. The optimization problem (2.1.4) searches for a subset U which approximates the observations with respect to the distance d. This reduces to the Fréchet mean problem for Q = M , in which case π U (x i ) = U for any point U ∈ Q. The problem (2.1.4) can be called an unsupervised approximation problem -unsupervised because there are no independent variables on which the observations depend, as there is in e.g. a regression problem. We think of U as a subset that represents the data given a chosen set of constraints, expressed by Q. For example, elements of Q might be of lower dimension than the ambient space of observations, and it might be a linear subspace if the ambient space is a vector space. For a particular problem, the question is then which family of subsets Q ⊂ P(M ) yield solutions with desired properties, e.g. properties similar to the case of Euclidean space. Principal component analysis, the foundation of the methods presented in this thesis, is another important example of problem (2.1.4).

Example 2.1.3 (Principal component analysis).

Let {x i } i=1...N ⊂ R d be observations in R d and let µ be their arithmetic mean. The two core steps of principal component analysis (see e.g. [START_REF] Jolliffe | Principal component analysis: a review and recent developments[END_REF]) are the following, 1. Compute the empirical covariance matrix

K = 1 N -1 N i=1 (x i -µ)(x i -µ) T . (2.1.6)
2. Compute the eigendecomposition K = P ΛP -1 , where Λ is a diagonal matrix containing the eigenvalues λ 1 ≥ • • • ≥ λ d and P is a matrix containing the corresponding eigenvectors e 1 , . . . , e d as columns.

Assuming that λ k is strictly greater than λ k+1 , we call V k . . = span{e 1 , . . . , e k } the k'th eigenspace. This algebraically defined object can be characterized in a

Riemannian geometry different way, namely as the k-dimensional subspace that best approximates the observations, in the sense of minimizing the following problem

V k = argmin U ∈Gr(k,R d ) N i=1 d(x i , π U (x i )) 2 , i = 1 . . . d -1, (2.1.7)
where Gr(k, R d ) is the Grassmannian manifold of all k-dimensional subspaces of R d [START_REF] Bendokat | A Grassmann manifold handbook: Basic geometry and computational aspects[END_REF]. This problem is of the form

(2.1.4) for Q = Gr(k, R d ), in which case π U (x i ) is orthogonal projection of x i to linear subspace U ⊂ R d .
This formulation of eigenspaces suggests that we can generalize them to more general metric spaces (M, d). However, the parameter space Q, the Grassmannian manifold of linear subspaces of R d , only makes sense when the space of observations is a vector space. A natural idea is then to look for a family of subsets Q of M analogous to linear subspaces of a vector space. In the case of Riemannian manifolds, such subsets, in fact submanifolds, can be built up from collections of geodesics -a generalization of straight lines to manifolds. This leads to Principal geodesic analysis and its approximative version tangent PCA which will be described in Sections 2.4.1 and 2.4.2.

Minimal prerequisites for approaching a problem of the form (2.1.4) on some metric space (M, d) seems to be the ability to compute the distance metric d, the existence of a suitable family of subsets Q for the problem at hand, and last, but not least, the ability to solve optimization problems. In this thesis we study the case where M is a complete Riemannian manifold and where the distance metric d is the geodesic distance. In the next section, we give an overview of how the tools available on Riemannian manifolds enables us to solve problems in non-Euclidean statistics, such as (2.1.4).

Riemannian geometry

In this section, we give a brief overview of the central objects in Riemannian geometry used for Riemannian geometric statistics. For more details, we refer the reader to, e.g., [START_REF] Lee | Introduction to smooth manifolds[END_REF]] or [START_REF] Do Carmo | Riemannian geometry[END_REF].

A Riemannian metric g on M is a smoothly varying family of inner products on the tangent bundle T M . Being an inner product on each tangent space T x M, x ∈ M , g enables measurement of lengths and angles of tangent vectors. From this, it defines the length of a smooth path φ : [0, 1] → M as L(φ) . . = 1 0 g( γ(t), γ(t))dt, and a distance between two points as the infimum of 9 2. Background: geometric statistics lengths of paths joining them;

d(x, y) = inf L(φ) γ : [0, 1] → M is smooth γ(0) = x, γ(1) = y , (2.2.1)
for any x, y ∈ M .

A crucial fact about (finite dimensional) Riemannian geometry, as opposed to sub-Riemannian geometry, is that the metric determines a canonical isomorphism between T p M and its dual space T ⋆ p M , called the flat map,

Z : T p M → T ⋆ p M : v → v Z . The inverse of the flat map is called the sharp map, denoted \ . Given a real- valued function f : M → R we can compute its differential at a point p ∈ M , which is a covector d p f ∈ T ⋆ p M .
Given a Riemannian metric, the induced sharp map thus provides a canonical way to identify this cotangent vector with a tangent vector, called the Riemannian gradient of

f at p, ∇f | p = \ (d p f ) ∈ T p M .
Recall that a connection on M allows to define the covariant derivative of e.g. a vector field or a tensor field along another vector field. A vector field is parallel along a curve if its covariant derivative along the curve is zero. For any points

x, y ∈ M , any tangent vector v ∈ T x M and any smooth curve φ betweeen x and y, there exists a parallel vector field along the curve, called the parallel transport of v along φ. The parallel vector field solves the parallel transport equation, a 1st-order linear ODE. A Riemannian metric induces a canonical connection, the Levi-Civita connection, as the unique connection satisfying the two conditions of torsion being constantly zero, ∇ X Y -∇ Y X = [X, Y ], and compatibility with the metric. The latter condition is equivalent to the parallel transport operator Π y

x : T x M → T y M along any smooth curve between x, y ∈ M being an isometry ([Lee 2018], lemma 5.2).

Geodesics are curves whose tangent vectors are parallel with respect to a chosen connection. This implies that they satisfy the geodesic equation, a 2ndorder ODE derived from the parallel transport equation. When the connection is Levi-Civita, determined by a Riemannian metric, geodesics have the property of being locally length minimizing curves. In particular, the path distance d on M is a geodesic distance, meaning that a shortest path between x, y ∈ M is a geodesic w.r.t. the Levi-Civita connection. Furthermore, the Hopf-Rinow theorem states that (M, d) is a complete metric space if and only if M is geodesically complete, meaning that solutions to the geodesic equations exist for all time.

The exponential map is the data-to-solution map for the geodesic equation and carries an initial point p ∈ M and initial direction v ∈ T p M to the value of Basic tools for Riemannian geometric statistics the corresponding geodesic γ v p at time 1, exp p :

T p M → M : v → exp p (v) = γ v p (1). (2.2.2)
The derivative of the exponential map at v = 0 is the identity map, so the inverse function theorem implies that exp is locally a diffeomorphism. Let C p ⊂ T p M be the largest subset containing 0 on which exp p is a diffeomorphism, called the injectivity domain. Let C p = exp p (C p ) ⊂ M be its exponential image. The inverse of exp on C p is called the Riemannian logarithmic map,

log p : C p → T p M : x → log p (x) = exp -1 p (x). (2.2.3) Since T p M ∼ = R d , log p thus defines a chart on M , called a normal chart.
Lastly, a crucial property of the log is its relation to the geodesic distance metric,

∥ log x (y)∥ x = d(x, y), (2.2.4)
where ∥ • ∥ x is the norm induced by the Riemannian metric g x at x.

Basic tools for Riemannian geometric statistics

In this section, we give an overview of how the objects and operations described in the previous section are used in Riemannian geometric statistics.

First of all, we mention that there is not necessarily a canonical choice of metric for a certain application. In the case of e.g. the LDDMM shape space (Section ??), there is a natural parameterized family of metrics, but not a canonical way to choose among these. Other examples are the manifolds of correlation matrices and of covariance matrices, where there is a plethora of different families of metrics (see the PhD thesis [START_REF] Thanwerdas | Riemannian and stratified geometries on covariance and correlation matrices[END_REF]] for an overview). This degree of freedom can be seen as modelling flexibility or a nuissance, depending on the application. In other cases, such as Kendall's shape space, there is a natural choice of metric.

After a choice of metric, the geodesic distance induced by the Riemannian metric turns the manifold M into a metric space (M, d) with respect to which distance based statistics can be performed. In practice, the distance d(x, y),

x, y ∈ M , is computed via the logarithm, implicitly assuming that y is in the exponential image of the injectivity area of x.

The gradient induced by the Riemannian metric allows for formulating optimization algorithms, e.g. gradient descent, for minimizing objective functions such as (2.1.3) for the Fréchet mean. See e.g. [START_REF] Boumal | An introduction to optimization on smooth manifolds[END_REF]] or [Absil, Mahony, 2. Background: geometric statistics and Sepulchre 2008] and Section ?? where we discuss gradient descent for a weighted Fréchet mean.

If closed-form expressions for geodesics are not available for a particular manifold, we can approximate them by numerically integrating the geodesic equation. This, and a wide range of other computations on Riemannian manifolds, is available in Python libraries such as geomstats [START_REF] Miolane | Geomstats: A Python Package for Riemannian Geometry in Machine Learning[END_REF]] and gaxgeometry (http://bitbucket.org/stefansommer/jaxgeometry).

Computing the log is a crucial problem for practical geometric statistics, e.g. for computing the geodesic distance and normal charts. On a general Riemannian manifold it is not possible to characterize the injectivity domain neither theoretically nor numerically. When we write log p (x i ) it will thus be under a, usually implicit, assumption that x i ∈ C p . In the case of no additional theoretical knowledge, in geomstats and jaxgeometry the logarithm is computed via a 'shooting method', which amounts to searching for a geodesic whose endpoint is as close as possible to the desired endpoint, w.r.t. some easy-tocompute ambient distance. Such an optimal geodesic is found via a numerical optimization scheme on the initial tangent vector. This is a relatively expensive optimization problem. We discuss this further in Appendix A.3.

Tangent spaces are important for Riemannian geometric statistics, and normal charts are a particularly convenient representation of these. Since T p M is a vector space approximating M around p, it provides us with the possibility to do Euclidean statistics on the transformed dataset log p (x i ) ⊂ T p M in a way that approximates exact statistics on the manifold. This comes with the benefit of being more computationally efficient, as well as removing the need for reformulating a method in an intrinsic way. Equation (2.2.4) expresses that the distance between the base point x and any point y ∈ C x is preserved in a normal chart. However, if neither x nor y is the base point then the corresponding distance in the tangent space is distorted -see subsection 2.3.1 for details. Theoretically, the basis induced by a normal chart on T p M is orthonormal. This is convenient, since Euclidean statistical formulas usually assumes an orthonormal basis. However, if the log is computed by an optimization in the tangent space, as described above, the basis is arbitrary, which then has to be taken into account. This and other points are not fully adressed in the literature. We discuss it further in Section 2.4.1 and in Appendix A.2.

Lastly, we mention the role of the parallel transport map along a geodesic. It is often necessary to transport computations done in one tangent space to another, or collecting tangent vectors from different tangent spaces into a common one, in which a computation on the set can be performed (see an example of this in Appendix A.4 on Taylor approximation of PGA). We therefore need suitable maps between tangent spaces. If the manifold is a Lie group, the differential of the group multiplication map (left or right) can be used. On a general Riemannian manifold, the parallel transport map provides solution.

Distortion of distances in the tangent space

As mentioned, there is a trade-off between the accuracy of the exact computation of geodesic distances and the convenient properties of the tangent space distances, such as computational efficiency. This trade-off is a recurrent theme in this thesis, so we elaborate on it in this section. Equation (2.2.4) says that the distance between the origin in a tangent space T p M and each log p (x i ) ∈ T p M , as measured with respect to the norm ∥ • ∥ p , equals the geodesic distance between p and x i . Thus, the tangent space does not distort these radial distances. However, the distance ∥v 1 -v 2 ∥ p between arbitrary points v 1 , v 2 ∈ T p M does not equal the corresponding geodesic distance d(exp p (v 1 ), exp p (v 2 )) on the manifold. In particular, the distortion of this geodesic distance depends on the curvature around p in the following way (see e.g. [Do Carmo and Flaherty Francis 1992], Section 5 proposition 2.7),

d(exp p (t • v 1 ), exp p (t • v 2 )) 2 = t 2 ∥v 1 -v 2 ∥ 2 p (2.3.1) -t 4 1 3 ⟨R(v 1 , v 2 )v 1 , v 2 ⟩ + O(t 4 ),
where R is the curvature tensor and ⟨R(v 1 , v 2 )v 1 , v 2 ⟩ p therefore the sectional curvature with respect to the plane spanned by v 1 , v 2 in T p M . Thus, in the case of positive sectional curvature (such as for any point p on the sphere), the Euclidean distances in the tangent space gets underestimated, and vice versa in the case of negative curvature. The formula also implies that the smaller t is, i.e. the closer the points being compared are, the smaller is the distortion. For Riemannian geometric statistics this means that the more concentrated, in the sense of low variability, the data set is, the smaller is the error incurred by using a distance based Euclidean formula in the tangent space. What a sufficient level of concentration is, relative to an accepted level of error, depends on the curvature of the manifold in the vicinity of the observations.

2. Background: geometric statistics

Dimension reduction

In this section, we focus on the particular problem of dimension reduction on Riemannian manifolds. We give details on the two foundational methods, tangent PCA and principal geodesic analysis. The methods presented in Chapters ?? and ?? are variations and combinations of these two methods. In the last subsection, we give an overview of more recent methods.

Tangent PCA

Tangent PCA can be described roughly as performing PCA on the logtransformed observations log p (x 1 ), . . . , log p (x N ) ∈ T p M , where p might be the Fréchet mean. However, this description leaves out some important theoretical and practical details that are not described in the literature, so we do that below.

Our presentation is structured according to the two steps of PCA: step one, computing the empirical covariance matrix, and step two, the eigendecomposition.

Intrinsic definition of the empirical covariance matrix

In Euclidean statistics, there is a canonical choice of an orthonormal basis. Maybe for this reason, it is rarely mentioned that the empirical covariance matrix (2.1.6) depends on the choice of basis, in the following sense. Let v A and v B be representations of a vector v in R d with respect to bases A and B, and let Q be the change of basis matrix from A to B, so that v B = Qv A . Changing the basis of v changes an outer product in the following way,

v B v B T = (Qv A ) (Qv A ) T = Qv A v A T Q T .
In contrast, the coordinate representation of φ A = v A v A T , considered as a linear map, is changed by

φ B = Qφ A Q -1 = Qv A v A T Q -1 .
Thus, if and only if Q is an orthogonal change of basis matrix, i.e. Q T = Q -1 , does changing the basis of v A lead to the desired change of basis of v A v A T , considered as a linear map. Eigendecomposition of a linear map, represented by a matrix, is independent of the choice of basis, in the sense that changing the basis of the linear map leads to a corresponding change of basis of the eigenvectors (the eigenvalues are not affected by basis change). Recall that the empirical covariance matrix is a sum of outer products. Thus, PCA depends on the choice of basis, up to orthogonal basis change, since the empirical covariance matrix does.

On a Riemannian manifold there is no canonical basis for a tangent space T p M . However, from the Riemannian metric one can compute a, non-unique, orthonormal basis and represent the vectors log p (x i ) i=1..N w.r.t. to this basis. The empirical covariance matrix can then be computed w.r.t. this orthonormal basis, and the discussion above implies that the result of the subsequent eigendecomposition will be independent on the particular choice of orthonormal basis. However, we can do better than this; we can give a coordinatefree definition of empirical covariance, which has the additional advantage of allowing computations to be performed in any basis -removing the need for computing, and transforming the data to, an orthonormal basis.

Following [START_REF] Pennec | A Riemannian framework for tensor computing[END_REF], we define the empirical covariance matrix at p ∈ M as the tensor

Σ(p) = 1 N -1 N i=1 log p (x i ) ⊗ log p (x i ) ∈ T p M ⊗ T p M. (2.4.1)
After choosing a local chart around p, which induces a basis A for T p M , the coordinate representation of this tensor is

[Σ(p)] A = 1 N -1 N i=1 log p (x i ) A log p (x i ) A T [g p ] A , ( 2
(C p ∩ V k ) is a k-dimensional
embedded submanifold of M . This submanifold can be used as a k-dimensional approximation of the data. Since the image of a linear subspace under the exponential map is formed by a collection of geodesics, we interpret the resulting submanifold of M as an analogue of a linear subspace.

Tangent PCA approximates distances We now discuss the sense in which tangent PCA is approximative. As discussed in Example 2.1.3 eigenspaces of a linear map can be formulated as solutions to a distance minimization problem. In the setting of tangent PCA, eigenspace V k , k = 1 . . . d-1, solves the optimization problem

V k = argmin U ∈Gr(k,TpM ) N i=1 d(log p (x i ), π U (log p (x i ))) 2 , (2.4.3)
where d(x, y) = ∥x -y∥ 2 p is the distance between x, y ∈ T p M w.r.t. the norm induced by g p . Thus, V k is a linear subspace approximating the log-transformed observations in T p M . However, as discussed in the previous section, the fact that V k is close to the observations {log p (x i )} i=1...N in the tangent space does not imply that the submanifold exp p (V k ∩ C p ) is close to the observations {x i } i=1...k in M . More specifically, that d(log p (x i ), π V k (log p (x i ))) is small does not imply that d(x i , exp p (π V k (log p (x i ))) is small. To ensure this, one has to minimize the geodesic distances instead of the tangent space distances. This is done in principal geodesic analysis, described in the next section.

In Chapter ??, we derive a tensor representation of a particular type of covariance matrix used in evolutionary biology, and define 'tangent phylogenetic PCA' as eigendecomposition of the corresponding matrix. In Chapter ?? we define socalled local tangent PCA, which are combined in a certain way to construct approximating submanifolds that does not consist of (Riemannian) geodesics.

Principal geodesic analysis

In this section we assume that the base point is the Fréchet mean, µ, as it is presented in the original paper on PGA [T. P. [START_REF] Fletcher | Principal geodesic analysis for the study of nonlinear statistics of shape[END_REF]]. As we saw, tangent PCA approximates the observations with the submanifold exp µ (V k ∩ C µ ) where V k is found by minimizing certain distances in the tangent space T µ M . In principal geodesic analysis, V k ⊂ T p M is found by minimizing the corresponding geodesic distances on M . The V k 's are obtained in the following sequential way. Let V 1 . . = span{v 1 }, where v 1 minimizes the following objective function, (2.4.5) where d is the geodesic distance on (M, g), and π is projection with respect to this distance. The following subspaces

v 1 ∈ argmin v∈TµM 1 N N j=1 d(x j , π exp µ (span(v)∩Cµ) (x j )) 2 (2.4.4) = argmin v∈TµM 1 N N j=1 ∥ log xj (π exp µ (span(v)∩Cµ) (x j ))∥ 2 ,
V 2 . . = span(v 1 , v 2 ), . . . , V d-1 . . = span(v 1 , . . . , v d-1
) are found by optimizing for one basis vector v i at a time, (2.4.6) for i = 2 . . . d-1. For practical purposes, to reduce the dimension of the parameter space, it can be helpful to add constraints enforcing unit length and orthogonality to the optimization problem, i.e. ∥v∥ g = 1 and v i ⊥ span{V i-1 }. Regardless, the PGA optimization problem is high-dimensional and computationally intensive.

v i ∈ argmin v∈TµM 1 N N j=1 ∥ log xj (π exp µ (Vi-1∩Cµ) (x j ))∥ 2 ,
The exact dimension and complexity depends on how many operations are available closed-form, e.g. log or exp. In Appendix A.3 we elaborate on the complications with computing PGA, and rewrite the objective function to something more tractable. No publicly available implementation of PGA exists to this date (for work in this direction, see [START_REF] Sommer | Manifold valued statistics, exact principal geodesic analysis and the effect of linear approximations[END_REF]). On a general manifold, as opposed to Euclidean space, the described sequential procedure of optimizing for one basis vector at a time is not equivalent to directly solving for the optimal k-dimensional subspace of T p M . Likewise, in Euclidean space, the eigenspace can equivalently be defined as the solution to a variance maximization problem, whereas this is not equivalent to error (distance) minimization on a general manifold. In practice, the sequential approach leads to a simpler optimization problem, and the formulation in terms of distance minimization seems to be more stable than variance maximization [START_REF] Sommer | Manifold valued statistics, exact principal geodesic analysis and the effect of linear approximations[END_REF].

We end our presentation of PGA by deriving tangent PCA as an approximation of PGA. The tangent PCA objective function (2.4.3) can be derived from the PGA objective function by making the approximation log µ (π exp µ (span(v)∩Cµ) (x j )) ≈ v • ⟨v, log µ (x j )⟩ with circles of radius 1, while the arbitrary subspheres of PNS can have any radius. Thus, the class of approximating submanifolds of the sphere considered by PNS is larger than that of PGA and G-PCA.

Barycentric subspace analysis (BSA) [START_REF] Pennec | Barycentric subspace analysis on manifolds[END_REF]] is a method which works for data on an arbitrary Riemannian manifold. This method approximates the data by the generalization to Riemannian manifolds of barycentric subspaces. Like PNS, BSA does not depend on a base point. It depends on a chosen set of points, e.g. observations, generating the barycentric subspace. By adding or removing points, the result is a nested sequence of barycentric subspaces of the original manifold.

All the methods mentioned so far consider a certain parameterized class of curves or submanifolds, and choose among these by minimizing a sum-of-squarederrors objective function, based on geodesic distance (although for BSA, the SSE approach is just one option). This corresponds to our earlier general formulation of the dimension reduction problem (2.1.4) based on a family of subsets Q of the original space. We now mention some methods which are not based on a parameterized family of subsets or on minimizing SSE. A method which is completely non-parametric, in the sense of the family of approximating curves being infinite dimensional, is the method of principal flows [START_REF] Panaretos | Principal flows[END_REF]. The method fits local tangent PCA's and generates a curve whose tangent vectors are aligned with the first eigenvector -i.e. it moves in the direction that locally describes the most variation in the data. This method only works for approximating the observations by a curve, but was generalized to higher dimensional approximations in [START_REF] Yao | Principal Sub-manifolds[END_REF]. Our method of principal subbundles, described in Chapter ??, can also be considered a generalization of the principal flow to higher dimensions (in Section ?? we elaborate on the relationship between these methods). An approach which is not based on deterministic curves, let alone geodesics, is infinitesmal probabilistic principal component analysis [START_REF] Pennec | Riemannian geometric statistics in medical image analysis[END_REF]] and [START_REF] Pennec | Riemannian geometric statistics in medical image analysis[END_REF]], Chapter 10), which is based on anisotropic Brownian motions on the manifold. This method uses maximum likelihood to fit a covariance matrix of lower rank than the ambient manifold, say k, describing the motion of the Brownian motions. This enables, among other things, a representation of the observations in a k-dimensional Euclidean space.

Chapter 3

Conclusion and future work

3.2

Future work related to phylogenetic PCA . . . . . .

3.3

Future work related to principal subbundles . . . . .

Summary of contributions

In Chapter 2 and its associated appendix A, we gave background on geometric statistics on Riemannian manifolds, with a special focus on computational aspects and on the two foundational methods for dimension reduction; Tangent PCA and Principal geodesic analysis. We gave details, not contained in the literature, on how to do statistics in the tangent space in a way that does not depend on the chosen coordinates. A special attention was given to the empirical covariance tensor, central to tangent PCA. Furthermore, we discussed and rewrote the objective function of PGA from a computational point of view, and derived a Taylor expansion of it.

In Chapter ??, we formulated a version of Tangent PCA adapted to applications in evolutionary biology; Tangent phylogenetic PCA. This method assumes that the data are leaf-nodes of a phylogenetic tree, and takes the implied correlation structure of such data into account. We defined, and derived estimators for, the associated notions of phylogenetic mean and covariance. This work is a generalization to Riemannian manifolds of a method, phylogenetic PCA, for Euclidean data. In studies of morphology in evolutionary biology the method is being applied to Procrustes aligned point configurations representing landmark shapes. In our version, we give a Riemannian gradient descent algorithm for estimating the intrinsic phylogenetic mean, and show that the way the mean is currently estimated in morphometrics corresponds to taking only a single step of this algorithm, when the Riemannian manifold is Kendall's shape space.

In Chapter ??, we developed a method for dimension reduction on a Riemannian manifold M which approximates the data based on a very flexible class of curves, and submanifolds composed of such curves. The method works by constructing a subbundle of the tangent bundle on the manifold M via local PCA's. We call this subbundle the principal subbundle. This subbundle induces a sub-Riemannian (SR) structure on M , and we show that sub-Riemannian geodesics with respect to this structure stay close to the set of observations. We show that it is possible to approximate the data set by submanifolds consisting of such SR geodesics, even if the subbundle is non-integrable. In particular, we show that the image under the SR exponential map of the dual space to the subbundle at a base point µ is an embedding into M , yielding a submanifold of M whose dimension equals the rank of the subbundle. Non-integrability is likely to occur when the subbundle is estimated from noisy data, and our method demonstrates that sub-Riemannian geometry is a natural framework for dealing which such noise. We also contributed to the software library geomstats [START_REF] Miolane | Geomstats: A Python Package for Riemannian Geometry in Machine Learning[END_REF], with a class representing a sub-Riemannian metric and associated computations of geodesic etc. This makes computations in sub-Riemannian geometry available for future applications in geometric, and Euclidean, statistics.

Future work related to phylogenetic PCA

Using anisotopic normal distributions Sommer et al. defines an anisotropic normal distribution [START_REF] Sommer | Modelling anisotropic covariance using stochastic development and sub-Riemannian frame bundle geometry[END_REF], [START_REF] Pennec | Riemannian geometric statistics in medical image analysis[END_REF]) on a Manifold M , equipped with a connection or a Riemannian metric, as the time-t transition distribution of an anisotropic Brownian motion on M . This distribution has mean µ, called the diffusion mean, and covariance Σ if the initial point of the Brownian motion is µ and its covariance is Σ (more precisely Σ is a frame, so (µ, Σ) is a point in the frame bundle on M ). The mean and covariance can be estimated by maximum likelihood methods (see [START_REF] Sommer | Anisotropic distributions on manifolds: Template estimation and most probable paths[END_REF]] and [START_REF] Grong | Most probable paths for anisotropic Brownian motions on manifolds[END_REF]). The diffusion mean, defined through a stochastic process, is an alternative to the Fréchet mean, defined via the geodesic distance. The framework of anisotropic normal distributions (ANS) seems natural to use in relation to phylogenetic comparative methods (PCM) since the underlying model for the latter is exactly based on such Brownian motions. To adapt the ANS framework to the PCM setting, however, requires to modify the framework to deal with Brownian motions structured according to a tree-graph. This is ongoing work. An off-the-shelf solution is to exchange the weighted Fréchet mean with the weighted diffusion mean (see [START_REF] Sommer | Horizontal flows and manifold stochastics in geometric deep learning[END_REF] and [START_REF] Jensen | Mean Estimation on the Diagonal of Product Manifolds[END_REF]).

Geometric morphometrics in relation to recent developments in computational geometric statistics

We see potential in revisiting the geometric morphometrics (GM) literature in view of the recent developments in computational differential geometry and geometric statistics. The basis of the methods in GM was developed before software libraries such as geomstats and jaxgeometry made available intrinsic computations on manifolds, in particular for geometric statistics. As discussed in Chapter ??, the field of geometric morphometrics builds on a wide range of approximations of intrinsic computations, e.g. different ways to compute local linear approximations to Kendall's shape space. It is possible that, in light of the mentioned recent computational developments, some of these are superfluous, or that some of them can be exchanged for more suitable, intrinsic methods. We hope that our work on phylogenetic PCA is one step in this direction.

A broader picture In continuation of the previous subsection on GM, we state the following long term goal: to formulate on Riemannian manifolds those phylogenetic comparative methods that are relevant for morphometrics and where intrinsic differential geometry can make a contribution. Such methods could be other basic statistical operations than PCA, e.g. phylogenetic generalized least squares regression [START_REF] Symonds | A primer on phylogenetic generalised least squares[END_REF]. Or it could be in the direction of allowing more flexible probabilistic models, e.g. based on processes more general than Brownian motions, and where the parameters are not necessarily the same throughout the tree (see [START_REF] Harmon | Phylogenetic comparative methods[END_REF]] and [START_REF] Mitov | Automatic generation of evolutionary hypotheses using mixed Gaussian phylogenetic models[END_REF]). Or it could be in the direction of uncertainty quantification, e.g. how does uncertainty in tree estimation propagate to uncertainty in the estimated parameters.

Choosing landmark positions that are comparable (homologuous) between species can be difficult, and the annotation of such datasets involves costly manual labour. Therefore it seems desirable to also develop phylogenetic comparative methods for 'full shapes', i.e. shapes defined by continuous curves and surfaces instead of landmarks, e.g. using the LDDMM framework for continuous curves or surfaces (see [START_REF] Younes | Shapes and diffeomorphisms[END_REF]]) or other shape representations (see [START_REF] Salili-James | Classifying organisms and artefacts by their outline shapes[END_REF] for a comparison of different representations in the context of shape classification). See [START_REF] Mitteroecker | Thirty years of geometric morphometrics: Achievements, challenges, and the ongoing quest for biological meaning-fulness[END_REF] for a discussion of potential pros and cons of 'landmark-free' shape representations in the context of biology and geometric morphometrics. However, geometric statistics is less developed for such infinite-dimensional manifolds compared to the finite dimensional setting treated in this thesis.

Future work related to principal subbundles

Evaluation We illustrated the framework of principal subbundles on a number of experiments in Chapter 4, but it was not fully benchmarked with respect to competing methods for manifold reconstruction, dimension reduction and metric learning. This is our most immediate direction of future work.

Improving computations

As mentioned in Appendix B.2.1, we have so far experimented with two ways of integrating the Hamiltonian equations; standard Euler integration and symplectic Euler. We plan to test other, e.g. higherorder, symplectic integrators, with the aim of decreasing computational cost (by allowing for increased step size) while keeping the error tolerable. One candidate could be the second-order leapfrog scheme (see e.g. [START_REF] Hairer | Geometric integration of ordinary differential equations on manifolds[END_REF]]), which is

Conclusion and future work

the most commonly used integrator in the field of Hamiltonian Monte-Carlo according to [START_REF] Betancourt | A conceptual introduction to Hamiltonian Monte Carlo[END_REF].

A second strategy for speeding up computations is to limit the number of points in which the local PCA is computed, in the following way. We suggest to precompute local PCA's only at observations, and interpolate between these by local averages. This strategy is used in surface reconstruction methods based on normals (in [START_REF] Kazhdan | Poisson surface reconstruction[END_REF] Poisson surface reconstruction is computed in this way). In that setting, the interpolation is simply a local average of normal vectors. In the higher dimensional case we need to compute an average of k-dimensional subspaces (see [START_REF] Marrinan | Finding the subspace mean or median to fit your need[END_REF] for an overview of different such averages).

A third strategy for speeding up computations is to lower the number of observations that are considered in each local PCA when integrating the Hamiltonian equations. Let p be a point on the ambient manifold N and choose a radius r > 0. For a moment, we will consider only observations that are inside the geodesic ball B r (p) ⊂ N . Choose a weight threshold ϵ > 0 such that the distance corresponding to this weight, r 0 = K -1 α (ϵ), is below r, where K α is the kernel function in of the weighted second moment, Definition ??. Then it holds that B r0 (q) ⊆ B r (p) for any q such that d(p, q) ≤ r -r 0 . This holds in particular for any point along a geodesic q t = exp p (tv) up to time t ≤ r -r 0 if ∥v∥ = 1. This means that, relative to the threshold ϵ, the observations contained in B r (p) are sufficient for computing the principal subbundle (local PCA) at points along this geodesic until time t ′ = r -r 0 . At q t ′ a new ball B r (q t ′ ) can be computed, and the process can be repeated. We think of this process as changing between 'charts' B r (•) containing a relevant subset of observations for a given time interval. We have presented the strategy for computing a single geodesic, but we expect a similar procedure to be useful when computing neighbouring geodesics generating a principal submanifold (Algorithm ??).

Extensions In this section, we mention two alternative ways to use a principal subbundle. The common theme is to consider other classes of curves than sub-Riemannian geodesics. This is the subject of ongoing work and is therefore only described cursorily.

One extension is to consider horizontal stochastic processes instead of sub-Riemannian geodesics. E.g. a sub-Riemannian equivalent of Brownian motions. The time-t transition distribution of such processes can act as an analogue to anisotropic normal distributions; they are probability distributions concentrated around the data. To perform maximum likelihood estimation for such distributions we need the sub-Riemannian equivalent of bridge-processes (a stochastic process conditioned on hitting a specified point), which is work-inprogress. Such bridge processes are also expected to be of independent interest, e.g. for generarative models, as a process that moves between two specified observations while following the point cloud.

Another extension aims at generating curves that are attracted to the point cloud in a different way than the sub-Riemannian geodesics described in Chapter ??. Sub-Riemannian geodesics determined by the principal subbundle are constrained to move horizontally w.r.t. the subbundle. I.e. they are subjected to constraints on their velocity vectors (as well as the constraint of being locally length-minimizing). When these curves are initialized within the point cloud they remain within it -at least in principle, and up to numerical error. Apart from this, they are not attracted to the data. It might be useful to introduce, in some way, an 'attraction term' that forces a curve initialized outside of the data to move towards it. This could also act as a correction to the deviation from the point cloud that is sometimes incurred by the numerical integration of the geodesic equations. We are currently considering multiple ways of adding such a data attraction term. One strategy is to add a term to the metric matrix which goes to infinity when evaluated outside of the point cloud, similar to what is done in [Hauberg 2018]. This changes the geometry of the space, in the sense of changing the metric. Another approach is to consider more general Hamiltonians than the sub-Riemannian Hamiltonian (Section ??)). This does not change the geometry, but only the solutions to the Hamiltonian equations, which are then no longer geodesics.

Statistics in the tangent space need to be computed for the same base point p ∈ M , such as in the gradient descent for a Fréchet mean (Algorithm ??, Chapter ??).

A.2 Statistics in the tangent space

A general strategy for performing statistics on data belonging to a Riemannan manifold, x 1 , . . . , x N ∈ M , is to map the data to a tangent space T µ M at a well chosen base point µ ∈ M , and apply Euclidean statistical methods in this vector space. To minimize distortion of distances due to tangent space linearization (see Section 2.3.1), it is preferable that the base point is relatively close to the observations. A Fréchet mean is a good candidate, since it is a minimizer of the sum of squared distances to the observations -note however that due to curvature it can be well outside the support of the data (this is the case for data distributed along a great arc on the sphere, for example). Having chosen µ, the next step is to map the observations to T µ M via the Riemannian logarithm, i.e. the data is represented as

log µ (x 1 ), . . . , log µ (x N ) ∈ T µ M.
The Euclidean statistical method can then be applied to the transformed dataset. For some methods, it is a reasonable final step to map the result back to M via the Riemannian exponential. In this way, e.g. a regression line in T µ M is mapped to a curve on M , or a subspace of T µ M found by e.g. PCA is mapped to a submanifold of M .

When following this strategy, one needs to take into account the fact that there is no canonical basis for a tangent space, and that Euclidean formulas usually assume the data to be represented in an orthonormal coordinate system. In various presentations of tangent space statistics (e.g. the tangent PCA approximation of PGA in [T. P. [START_REF] Fletcher | Principal geodesic analysis for the study of nonlinear statistics of shape[END_REF]) is it implicitly assumed that the observations are represented in a basis induced by a normal chart. This basis is orthonormal -the local representation of the Riemannian metric is the identity matrix. However, if the chart does not induce an orthonormal basis on the tangent space, one needs to change to one, or take the non-orthonormal basis into account. As described in Appendix A.1, when the log is computed numerically by the optimization problem (A.1.1), the basis chosen for the tangent space is arbitrary, not necessarily orthonormal.

To find an orthonormal basis for T µ M , one can do the following. Let p ∈ M and let g be the Riemannian metric on M , represented in local coordinates around p. Let LL T = g -1 (µ) be the cholesky decomposition of the cometric. Then the columns of L form an orthonormal basis for T p M .

In the next section, we focus on the particular case of tangent PCA, i.e. Euclidean PCA performed in a tangent space.

A.2.1 Expressing the second moment in coordinates

The geometric statistics literature is unclear regarding the meaning of the covariance matrix defined as the tensor product (2.4.1), and what its coordinate expression is. Therefore, in this section, we show how to derive the coordinate expression (2.4.2) from (2.4.1).

For some v, u ∈ T p N , the expression v ⊗ u can be identified with an endomorphism on T p N . Its coordinate representation is thus a d × d matrix. There seems to be some confusion about this in the geometric statistics literature, so we give details below. For the remainder of this section we denote a Riemannian metric by h, to align the notation with Chapter ??, where multiple metrics are at play. Lemma A.2.1. Let (N , h) be a Riemannian manifold, and u, v ∈ T p N . Given a choice of basis for T p N , the tensor v ⊗ u ∈ T p N ⊗ T p N can be expressed in coordinates as

vu T h p ∈ R d×d , (A.2.1)
where u, v ∈ R d×1 are the vectors and h p ∈ R d×d is the Riemannian metric represented w.r.t. the chosen basis.

Proof. The tensor v ⊗ u is an element of the tensor product space T p N ⊗ T p N . After choosing a Riemannian metric, there is a canonical isomorphism between T p N and its dual space, T ⋆ p N , given by the Riemannian flat map,

Z : T p N → T ⋆ p N : u → h p (u, •) := u Z . Thus T p N ⊗ T p N ∼ = T p N ⊗ T ⋆ p N
, where elements of the latter space are denoted (1, 1) tensors. Furthermore, there is a canonical isomorphism, independent of a Riemannian metric,

T p N ⊗ T ⋆ p N ∼ = End(T p N ),
Computing principal geodesic analysis

≈ argmin α∈R N v∈Tµ M W ∈R d×N 1 N N j=1 ∥W :j ∥ 2 g + λd 0 (exp exp µ (αj v) (W :j ), x j ) 2 , (A.3.8)
where W is a matrix whose columns are the log-vectors, one for each observation.

λ is the same trade-off parameter as in the log-optimization problem (A.1.1), controlling the emphasis on minimizing approximation error in the log relative to minimizing the length of the corresponding geodesic. d 0 is a distance that can be computed easily -henceforth we will assume that d 0 (x, y) = ∥x -y∥ Eucl is the Euclidean distance and that the manifold is embedded in a Euclidean space.

The interpretation of the optimization problem (A.3.8) can be problematic, in the following sense. An optimization scheme will produce a sequence of parameter candidates in the parameter space R (d+1)N +d , but for each update of parameters (α j ) 0 ∈ R, j = 1 . . . N , and v 0 ∈ T µ M the vector W :j ∈ R d is used as a tangent vector belonging to a different tangent space, T exp µ ((αj )0v0) M . On a general Riemannian manifold (that is not e.g. a Lie group), there is no canonical way to identify tangent spaces. One way to map between different tangent spaces is the parallel transport map, Π y x : T x M → T y M . This map is an isometry, meaning that it preserves angles and lengths of vectors. Parallel transport depends on a chosen curve between points x, y, which we will assume to be the length-minimizing geodesic (assumed unique) between x, y. Therefore, we can use parallel transport to ensure that the vectors W :j , j = 1, . . . , N , represent tangent vectors in the same tangent space, namely T µ M . Let α j v . . = exp µ (α j v) ∈ M . Then this amounts to the following version of the optimization problem,

argmin α∈R N v∈Tµ M W ∈R d×N 1 N N j=1 ∥W :j ∥ 2 g + λ∥exp αj v (Π αj v µ W :j ) -x j ∥ 2 Eucl , (A.3.9)
The parameter space of the final optimization problem (A.3.9) for the first principal direction v 1 is of dimension (d + 1)N + d. For the i'th principal component, solving (2.4.6), the projection problem causes the projection parameter α to grow in dimension; α ∈ R i•N , since a point in an i dimensional subspace is parameterized by i coordinates. Thus the parameter space dimension for the i'th component is (d + i)N + d. Evaluating the objective function involves solving the parallel transport equation once and computing two exponential maps, i.e. integrating the geodesic equations twice. The authors have implemented this scheme with promising results for lower dimensional parameter spaces, but more work needs to be done in testing it for higher dimensions.

The map (B.1.1) takes a point p ∈ N to a vector field at time 1 satisfying equations (B.1.4)-(B.1.9). For each p, the initial conditions are

x = f (p) Q = log h f (p) (p) W = O(p)
all of which depend smoothly on p, if p / ∈ Cut(f (p)). Since the solution to the ODE depends smoothly on the initial conditions, and since the initial conditions depends smoothly on p, the vector field (B.1.1) is smooth. * Proof. As in the Euclidean case, we want to prove the existence of a smooth frame around every point p ∈ S ′ α,k spanning the subbundle locally around p. We will make use of the corresponding result for N = R d , in a chart. In order to do this, we need to make sure that all of the involved maps are smooth as a function of p. speed. We compared ordinary Euler integration to semi-implicit Euler (see e.g. [START_REF] Hairer | Geometric numerical integration[END_REF]), a first-order symplectic integrator, and found the Hamiltonian to be better preserved using ordinary Euler integration in our experiments.

B.3 Choosing the kernel range α and bundle rank k

Firstly, note that these parameters can be considered to be a modelling choice, expressing the scale at which we want to analyze the data -what scale of variation to take into account. However, one can aim to find the 'lowest level of variation that is not due to random noise'. Secondly, note that the 'optimal' value of one hyperparameter depends on the value of the other. Since the rank k takes a finite number of values k ∈ {1, . . . , d -1}, we suggest to start by estimating this. See [START_REF] Bac | Scikit-dimension: a python package for intrinsic dimension estimation[END_REF]] for a survey and benchmarking of different methods. Given an estimated k, we suggest to select a range parameter for which the separation between eigenvalues λ k and λ k+1 is the most clear on average. The optimal kernel range depends on the level of noise and the rate of change of the affine subspace E p as a function of p, which, in the case of the manifold hypothesis, is an expression of the curvature of the underlying manifold. A fast varying E calls for a smaller α, while high levels of noise as well as a lower number of observations calls for a larger α.

B.4 Algorithm for combining principal submanifolds for 2D surface reconstruction

In this section, we present an algorithm for combining principal submanifolds {M k µj (r j )} j=1..l based at different base points µ j , j = 1 . . . l. In this case, k = 2 and we'll write M µj instead of M 2 µj . Given a point x ∈ R 3 , the algorithm first projects x to a set of nearest principal submanifolds, and then represents x as a weighted average of these projections, weighted by the SR distance between a projection and its corresponding base point. The point x can e.g. be an observation, x ∈ {x i } i=1..N , or a point in a principal submanifold, x ∈ M µj . The algorithm can then be run for each point x in {x i } i=1..N or in M µj , j = 1..l.

The point sets representing principal submanifolds M µj (r j ), j = 1 . . . l, are generated by Algorithm 1. For each point p ∈ M µj (r j ), we assume that the corresponding initial cotangent η(p) ∈ E ⋆ µj has been stored. Apart from the hyperparameters of the principal subbundle and submanifolds, the algorithm needs a 'threshold parameter' ϵ > 0. x will not be projected to .3 shows a single computed geodesic, found by solving the log problem log p (q), for p, q and observations as described in section??. The distance d E (p, q) is estimated as the length of the computed geodesic. The blue points are observations on the 4-dimensional sphere embedded in R 5 0 projected to R 3 .

.3: Illustration of a single computed geodesic found by solving the log problem log p (q), for p, q and observations as described in section ??. The blue points are observations on the 4-dimensional sphere embedded in R 50 projected to R 3 . The orange point is the initial point p, the red points are points on the geodesic, the green point is the target point q.

Contents

  Summary of contributions . . . . . . . . . . . . . . .

  The tangent mean map m : N → N and the tensor field p → Σ α (p) ∈ T p N ⊗ T p N is smooth if each logarithm log h p (x i ), i = 1 . . . N , is smooth as a function of the base point p ∈ N . This is ensured by the cut locus conditions inS ′ α,k .Assuming smoothness of Σ α , we now consider charts (U, φ) on N and (O, ϕ)on T N ⊗ T N , U ⊂ R d , φ : U → φ(U ) ⊂ N , respectively O ⊂ R d×d , ϕ : O → φ(O) ⊂ T N ⊗T N (identifying each T p N ⊗T p Nwith the space of endomorphisms on T p N , cf. Section A.2.1), around a point p ∈ N and φ(p) ∈ T N ⊗ T N . In this chart,f := ϕ -1 • Σ α,k • m • φ is a smooth map from R d to R d×d . Eigendecomposition of the matrix f (p ′ ), p ′ ∈U , is independent of the basis and thus of the choice of charts. As shown in the proof of Proposition ??, there exists a smooth framep ′ → {X 1 (p ′ ), . . . , X k (p ′ )}, X i (p ′ ) ∈ R d , defined on some open subset U ⊂ R d around φ -1 (p) s.t. span{X 1 (p ′ ), . . . , X k (p ′ )} = V k (f (p ′ )), ∀p ′ ∈ U ,where the right hand side is the eigenspace of f (p ′ ) corresponding to the largest k eigenvalues. We have thus shown the existence of a smooth frame on φ(U ) ⊂ N spanning the corresponding eigenspaces of Σ α • m at every point of φ(U ).

Figure B. 1 :

 1 Figure B.1: Supplementary figure to Figure ??, zooming in on the region where the two submanifolds overlap. The left, beige submanifold in this figure is the purple one in Figure ??, recolored to be able to distinguish more clearly the two submanifolds.

  denne afhandling praesenteres nye metoder til dimensionsreduktion, baseret på differential geometri. Det vil sige metoder til at finde en repraesentation af et datasaet i et rum af lavere dimension end det opringelige rum. Sådanne metoder spiller en helt central rolle i statistik, og har et meget bredt anvendelsesområde. En laveredimensional repraesentation af et datasaet tillader visualisering og er ofte nødvendigt for efterfølgende statistisk analyse. I traditionel, Euklidisk statistik ligger observationerne i et vektor rum, og det lavere-dimensionale rum kan vaere et lineaert underrum eller en ikke-lineaer undermangfoldighed som approksimerer observationerne. Studiet af sådanne glatte mangfoldigheder, differential geometri, spiller en vigtig rolle i sidstnaevnte tilfaelde, eller hvis rummet hvori observationerne ligger i sig selv er en mangfoldighed. Metoder til at analysere observationer på en mangfoldighed udgør feltet geometrisk statistik. I denne kontekst er det approksimerende rum, fundet via dimensionsreduktion, naturligt en submangfoldighed af den givne mangfoldighed. Udgangspunktet for denne afhandling er geometrisk statistik for observationer på en a priori kendt Riemannsk mangfoldighed, men dele af vores arbejde udgør et bidrag selv i tilfaeldet med observationer i Euklidisk rum, R d .

-clés: statistiques géométriques, géométrie différentielle, géométrie Riemannienne, géométrie sous-Riemannienne, statistique mathématique, apprentissage automatique

Riemannske og sub-Riemannske metoder til dimensionsreduktion I

  which has been submitted to the journal Foundations of computational mathematics in July 2023. In April 2023, this work was presented by the author at the workshop Statistics of shapes and Geometry of Shape Spaces, at the Max-Planck-Institute for Mathematics in the Sciences, Leipzig. In May 2023 it was presented by the author at the Analysis and PDE Seminar at the university of Bergen, Department of mathematics.

	Chapter 2
	Background: geometric statistics

Contents 2.1 Generalizing Euclidean statistical methods . . . . . 2.2 Riemannian geometry . . . . . . . . . . . . . . . . . 2.3 Basic tools for Riemannian geometric statistics . . . 2.3.1 Distortion of distances in the tangent space 2.4 Dimension reduction . . . . . . . . . . . . . . . . . . 2.4.1 Tangent PCA . . . . . . . . . . . . . . . . 2.4.2 Principal geodesic analysis . . . . . . . . . 2.4.3 Other methods for dimension reduction on Riemannian manifolds . . . . . . . . . . .

  Background: geometric statistics in the tangent space. The log transformed observations can be projected to eigenspaces, π V k log p (x i ) , and mapped back to M via the exponential map, exp p π V k log p (x i ) . Since the exponential map is a diffeomorphism on the tangent cut locus C p it holds that exp p

	is the Fréchet mean then the covariance (2.4.2) equals the ordinary Euclidean
	empirical covariance formula (2.1.6).
	Eigendecomposition in the tangent space Given this intrinsic definition
	of empirical covariance on a Riemannian manifold, we can move to step 2 of

.4.2)

where

[•]

A denotes representation of an object with respect to basis A. See Appendix A.2.1 for a derivation of this coordinate expression. Thus, the coordinate representation of the tensor Σ(p) is a d × d matrix as we would expect. If the basis A is orthonormal with respect to metric g, then [g p ] A is the identitity matrix. Thus, the common description of tangent PCA as eigendecomposition of the matrix 1

N -1 N i=1 log p (x i ) A log p (x i ) A

T is only correct if the basis A is orthonormal. In theory this is the case when the basis is induced by a normal chart. However, if the log-map is computed numerically via an optimization scheme it is not necessarily the case (see Appendix A.1). Lastly, note that if (M, g) is Euclidean space then log p (x) = x -p, so that if p PCA, the eigendecomposition. Tangent PCA consists simply of performing eigendecomposition of Σ(p). The eigendecomposition yields eigenvectors e 1 , . . . , e d ∈ T p M and thus eigenspaces V k . . = span (e 1 , . . . , e k ) ⊂ T p M 2.

Acknowledgments

This PhD was carried out in France and in Denmark, with the time split equally between each place. The work was supported by the European Research Council (ERC) under the EU Horizon 2020 research and innovation program (grant agreement G-Statistics No. 786854). I'm grateful to have met and worked with Xavier, Stefan, James and Erlend. Thank you for being generous, inspiring and knowledgeable.

which leads to d(x j , π exp µ (span(v)∩Cµ) (x j )) ≈ ∥ log µ (x j ) -log µ (π exp µ (span(v)) (x j ))∥ p (2.4.7) ≈ ∥ log µ (x j ) -v • ⟨v, log µ (x j )⟩∥ p (2.4.8)

≈ ∥ log µ (x j ) -π span(v) (log µ (x j )))⟩∥ p (2.4.9) where the first approximation is based on Equation (2.3.1), and π span(v) is orthogonal projection in T p M w.r.t. the norm induced by g µ . The latter expression is exactly what tangent PCA minimizes, cf. Equation (2.4.3). In Appendix A.4 we derive a Taylor expansion of the objective function of PGA whose first order term corresponds to the tangent PCA objective function -this is another way of deriving tangent PCA as an approximation to PGA. The higher order terms in the Taylor expansion takes the curvature into account.

Since PGA is based on computing exact distances in M , as opposed to approximated distances in T µ M , it is more sensitive to the curvature of M than PGA. However, it is not sensitive to potential non-geodesic variation in the dataset. This is in line with it being a manifold generalization of PCA, which assumes that the data is well approximated by linear subspaces. The method we present in Chapter ?? picks up on non-geodesic variation in the data, by combining PGA and local tangent PCA's.

Other methods for dimension reduction on Riemannian manifolds

We round off by a short survey of other, more advanced, methods for dimension reduction on Riemannian manifolds. For another, short review, see Section 5 of [S. F. [START_REF] Huckemann | Data analysis on nonstandard spaces[END_REF].

The most direct extension of PGA is Geodesic principal component analysis (G-PCA) [S. [START_REF] Huckemann | Intrinsic shape analysis: Geodesic PCA for Riemannian manifolds modulo isometric Lie group actions[END_REF], which also approximates the observations by geodesics, but which furthermore optimizes for the base point µ. The method of Principal nested spheres (PNS) [START_REF] Jung | Generalized PCA via the backward stepwise approach in image analysis[END_REF] works for data belonging to spheres in any dimension (see also principal nested symmetric spaces [START_REF] Curry | Principal symmetric space analysis[END_REF] for a generalization). The observations are approximated by subspheres of the original sphere, minimizing a sum-ofsquared-errors (SSE) criterion where the error is measured by geodesic distance. By iteratively fitting a subsphere of one dimension less, the output is a nested sequence of subspheres (submanifolds). In comparison, the result of PGA is a nested sequence of submanifolds that are constructed by adding a dimension in each iteration. Geodesics on the sphere are great arcs, which can be identified

Appendices

Appendix for Chapter 2 A.1 Computational aspects of geometric statistics

For doing statistics on a Riemannian manifold we need to be able to compute operations such as those described in Section 2.3. For a given manifold, closed form expressions may or may not be available -most often numerical approximations are needed at some stage of a computation. Fortunately, a large effort has been put into implementing numerical differential geometry and geometric statistics in recent years, with software libraries such as geomstats [START_REF] Miolane | Geomstats: A Python Package for Riemannian Geometry in Machine Learning[END_REF] and jaxgeometry (http://bitbucket.org/stefansommer/jaxgeometry).

Two main computational bottlenecks for geometric statistics is computation of the exponential and logarithmic maps. The former is computed by numerical integration of the geodesic equation in a chart, either as a 2nd order ODE in the tangent bundle or as two coupled first order ODE's in the cotangent bundle (i.e. the Hamiltonian geodesic equations, see ??). The logarithm can be computed by a shooting method, i.e. by solving an optimization problem in the tangent space that searches for the initial tangent vector minimizing some distance between the endpoint of the candidate geodesic and the desired endpoint. This distance should be efficient to compute, e.g. an ambient Euclidean distance if we have available an embedding of the manifold in a Euclidean space. In more detail, the optimization problem is the following,

The term λ∥v∥ g facilitates finding a length-minimizing geodesic, with the hyperparameter λ > 0 determining how much emphasis is placed on this length-minimizing property compared to the 'matching term' d(exp p (v), y). In Geomstats and in Jax geometry, this optimization problem is solved using the BFGS algorithm [START_REF] Wright | Numerical optimization[END_REF]] with automatic differentiation. Thus, the objective function in the optimization problem for the log map involves computing the exponential map, i.e. solving the geodesic equations. This is computed at least once at each step of the optimization algorithm. This has led us to develop a discrete approximation of the log map, presented in the following section. We use this approximation in Chapter ??, e.g. when computing projections to a principal submanifold in Section ??, and we expect it to be of wider interest for computations on Riemann manifolds.

Lastly, we note that when the log is computed by solving (A.1.1), the basis for the tangent space is arbitrary, in particular it is not necessarily orthonormal as one might expect if considering the computed log to be a point in a normal chart.

A.1.1 Discretization of the log map

We propose a discrete approximation of the Riemannian logarithmic map. The idea is to pick a finite set of relatively evenly distributed tangent vectors on the unit sphere in T p M . For each such tangent vector we integrate the geodesic equations up to a chosen time point r > 0. For each point along a geodesic initialized by one of these tangent vectors, we store its time index and its position in M . The collection of all such positions forms a discretized geodesic ball of radius r in M . For a given point x ∈ M , log p (x) is then approximated by the tangent vector corresponding to the nearest point on the generated grid. 'Nearest' is here defined w.r.t. a distance metric d 0 that should be efficient to compute, e.g. an ambient Euclidean metric, as for the objective function (A.1.1). The procedure is described in 3 steps below.

1. Generate an even grid S d ⊂ S d ⊂ T p M on the unit sphere (w.r.t. the Riemannian meric g) in the tangent space at p, and a grid of timepoints [0, r] ⊂ [0, r], for some final time r > 0.

2. Generate a discrete ball B p of radius r > 0 around p on M by

For each point y v0 t0 . . = exp p (t 0 v 0 ) ∈ B p , store the parameters of the corresponding initial tangent vector, i.e. t 0 , v 0 .

3. Define the discrete log, log p (y), by log p (y) ∈ argmin

where λ is a trade-off parameter, as described in relation to Equation (A.1.1).

The optimization problem (A.1.2) is discrete. It can be vectorized, such that the only computationally heavy part is the initialization steps 1 and 2. We expect the discrete log to be the most useful in situations where 'many' logs, log p (•),

where End(T p N ) is the space of endomorphisms on T p N . This isomorphism is given by the map Φ which takes an endomorphism A to the (1, 1) tensor Φ(A) that acts on w ∈ T p N and η ∈ T ⋆ p N by Φ(A)(w, η) = η(Aw). The linear map corresponding to a (1,1) tensor of the form

After choosing a basis for T p N , the tangent vectors v, w can be represented as column vectors v, w ∈ R d×1 . The flat map can be represented by the matrix h p , which is the matrix representation of the Riemannian metric at p. After identifying covectors with row vectors (i.e. coordinate representations of linear maps from T p N to R), u Z can be represented as the row vector u Z = (h p u) T ∈ R 1×d . This acts on w by u Z (w) = (h p u) T w. Thus, w.r.t. some chosen basis, the matrix representation of our desired endomorphism is given by

A.2.1.1 Verifying independence of the coordinate system

Let Q be the change-of-basis matrix from basis a of T p N to basis b. Then

, where these bases are dual to a, b. Thus, the change of basis of tangent

The flat map Z is a linear map from T p N to T ⋆ p N , so if (h p ) a is its representation w.r.t. bases a and a ⋆ , then its representation w.r.t. bases b and b ⋆ is computed as

We verify that the change of basis of the individual elements u, v, h p matches the change of basis of the matrix (as a linear map) (A.2.1):

As opposed to this, the expression v b u T b does not transform properly under basis change:

if the basis change matrix is orthogonal, meaning that it only rotates the basis.

A.3 Computing principal geodesic analysis

In this section, we elaborate on the optimization problem (A.3.1) that defines principal geodesic analysis (PGA). We will see that, if no closed form expressions for the Riemannian log and projection w.r.t. to geodesic distance are available, then it is a triply nested optimization problem. As of this date, to the best of our knowledge, there exists no publicly available implementation of PGA. The first principal direction solves

The projection problem π can be rephrased as

The log is given by the optimization problem (A.1.1). The outer optimization for v ∈ T µ M thus makes this a triply nested optimization problem. Automatic differentiation [START_REF] Margossian | A review of automatic differentiation and its efficient implementation[END_REF]] provides an efficient way to compute gradients for solving optimization problems, and packages like geomstats and jaxgeometry relies heavily on it. At the time of writing, automatic differentiation of such 'implicitly' defined functions, where the output is the solution to an optimization problem, is still at an early stage of development. An alternative approach is to rewrite the problem into a large global optimization instead of a nested one.

We first move out the inner optimization for the projection,

We then move out the inner optimization for the log,

A.4 Taylor-approximation of PGA

In this section we derive a Taylor approximation of the PGA objective function.

The expansion contains 3 terms, the first of which corresponds to the objective function of tangent PCA. The subsequent terms takes curvature into account.

The expansion is derived by expressing the PGA objective function in terms of the socalled neighbouring log map, which is a map with a known Taylor expansion. The neighbouring log and its Taylor expansion was first introduced by Pennec in [START_REF] Pennec | Riemannian geometric statistics in medical image analysis[END_REF]] building on work by [START_REF] Gavrilov | Algebraic properties of covariant derivative and composition of exponential maps[END_REF]], [Gavrilov 2007]).

Definition A.4.1 (The neighbouring log map). Let p be a point on M , and let v, w be tangent vectors in T p M . Let Π y x (v) be the parallel transport of v from T x M to T y M along the (assumed unique) length-minimizing geodesic from

x ∈ M to y ∈ M . Let p u . . = exp p (u). Then the neighbouring log map is given by

The Taylor expansion of the neighbouring log is given by

where R is the curvature tensor and ∇ • R its covariant derivative. See [START_REF] Pennec | Riemannian geometric statistics in medical image analysis[END_REF]] for details.

We now formulate the PGA objective function, for the first principal component, in terms of the neighbouring log map. The objective function for the subsequent components can be found similarly. Below, we will use the notation x j . . = log µ (x j ), and ∥ • ∥ p for the norm on T p M induced by that Riemannian metric, and ⟨•, •⟩ p for the Riemannian metric at p ∈ M .

Taylor-approximation of PGA = argmin

Equation (A.4.4) uses the fact that

and Equation (A.4.6) uses the fact that parallel transport is an isometry. Thus, the expansion consists of 3 terms. Including only the first order term yields the objective function of tangent PCA. The higher order terms take curvature into account.

The complexity of the approximated optimization problem is comparable to 'tangent PCA', thus relatively inexpensive. In particular, the log-map needs to be solved only once for every observation. In contrast to this, the exact PGA problem is a nested optimization problem, as we discussed in Section A.3.

Appendix for Chapter ??

B.1 Proofs

B.1.1 Smoothness of the principal subbundle

We show smoothness first on R d and then on a Riemannian manifold (N , h). The proof of the latter utilizes the former result in a chart, as well as smoothness results for the involved maps, which are only non-trivial in the manifold case. *

Proof. Let p ∈ R d \ S α,k be arbitrary. We will show that there exists a local frame of smooth vector fields spanning the subspace E α,k p ′ at every point p ′ on an open set U around p. By Lemma 10.32 in [START_REF] Lee | Introduction to smooth manifolds[END_REF]], this is equivalent to the subbundle being smooth on R d \ S α,k .

The eigenvalues of Σ α (p) at p are

where only λ k and λ k+1 are assumed to be different. Since Σ α : R d → R d×d is a smooth map, Theorem 3.1 of [START_REF] Sun | Multiple eigenvalue sensitivity analysis[END_REF]] implies that there exists an open set

Since each λi is continuous, there exists an open subset U ⊂ B on which the ordering λ1 (p ′ ) ≥ • • • ≥ λd (p ′ ) holds for all p ′ ∈ U , and where λi (p ′ ) = λj (p ′ ) is only possible for i, j s.t. λi (p) = λj (p). In particular λi (p ′ ) < λk+1 (p ′ ) for all i < k + 1 and p ′ ∈ U .

Theorem 3.2 of [START_REF] Sun | Multiple eigenvalue sensitivity analysis[END_REF]] now says that there exists a frame of analytic vector fields p → {X 1 (p), . . . , X k (p)} such that, for all p ′ ∈ U ,

denotes the eigenspace of Σ α (p ′ ) corresponding to eigenvalues λ1 (p ′ ), . . . , λk (p ′ ), which is exactly the principal subbundle subspace

To show that the principal subbundle on a Riemannian manifold is smooth, we need a result on smoothness of a certain map involving parallel transport.

Proofs

Lemma B.1.1. Let the map f : N → N and the vector field O on N be smooth. Let Π y x : T x N → T y N denote parallel transport along the (assumed unique) length-minimizing geodesic from x to y. Then the vector field

is smooth for every p / ∈ Cut(f (p)).

Proof. For x, y ∈ N , the parallel transported vector Π y x W ∈ T y N of W ∈ T x N along a curve γ : (0, 1) → N is the value at time 1 of a vector field V along γ satisfying the linear initial value problem (an ODE)

where Γ k ij , i, j, k ∈ {1, . . . , d} are the Christoffel symbols determined by the metric h. See [START_REF] Lee | Introduction to smooth manifolds[END_REF]], Section 4, for details.

If γ is a geodesic with initial velocity Q ∈ T x N then it is a solution to the geodesic equations (equations (B.1.5) and (B.1.6), below). In this case, we can write the parallel transport equation and the geodesic equations as a single, coupled, ODE:

Note that the equation for V is coupled with the equations for γ and U , but not vice versa, so that, in practice, the whole path γ can be computed first, and then subsequently V . This is again a linear initial value problem, and the fundamental theorem for ODE's states that solutions exist, and depend smoothly on the initial conditions Q, W, x. This shows smoothness of the parallel transport operator in the case where γ((0, 1)) is contained in a single chart. For the more general case, we refer to the technique used in the proof of Proposition 4.32 in [START_REF] Lee | Introduction to smooth manifolds[END_REF] for showing that solutions found on individual charts overlap smoothly.

Notes on implementation

The last thing we need to take account of is the parallel transport map. Since parallel transport is an isometry, it holds that

, where {F } i=1..k is any other frame spanning the same subspace as {X} i=1..k at p ′ . Thus, the parallel transported frame X spans the same subspace as the parallel transported eigenvectors {e i } i=1...k at p ′ (the X i 's are not necessarily eigenvectors, as explained in [START_REF] Sun | Multiple eigenvalue sensitivity analysis[END_REF]. By Lemma B.1.1, the map p → Π p m(p) V (p) is smooth, for a smooth vector field V . We have thus shown that the principal subbundle at p is spanned by a smooth frame around p.

B.2 Notes on implementation

At each step of the integration of a geodesic, eigenvectors needs to be computed at the current position p. This involves evaluating the kernel K α (|x i -p|) for all i = 1..N . For large datasets, we suggest to do this using libraries specialized at such kernel-operations, such as KEOPS, as well as automatically filtering out points far away from p whose weight will be close 0 anyway. We have not had the need to implement these optimizations in order to run the examples of Section ??.

The integration of the L geodesics in the algorithm for the principal submanifold can be parallelized; the computation of each one is independent from the rest.

To speed up computations further, we suggest to compute the sub-Riemannian metric at p ∈ R d as a weighted mean of the metric computed at a finite number of points around p, possibly at every observation. This is similar to the approach used to compute Poisson surface reconstruction [START_REF] Kazhdan | Poisson surface reconstruction[END_REF] (where the surface normal at p ∈ R 3 is evaluated as the mean of normals computed nearby p) and the learned Riemannian metrics in [START_REF] Hauberg | A geometric take on metric learning[END_REF]. In this way, the derivatives of the metric, and therefore the Hamiltonian, can be computed closed form, removing the need for automatic differentiation.

B.2.1 Choice of integration scheme

The integration of Hamilton's equations can be done using a symplectic integration scheme which aims at keeping the Hamiltonian constant. A constant hamiltonian is equivalent to constant speed, cf. eq. (??). This is desired because the computation of curve length and distance via eq. (??) assumes constant Supplementary figures principal submanifold M µj if the distance between x and its projection xj to M µj is greater than ϵ. Thus, the size of ϵ should be comparable to an estimate of the noise-level in the point cloud.

The algorithm is the following.

1. Project to each submanifold: project x to each M µj (r j ), j = 1..l, wrt.

Euclidean distance, i.e. find the closest point in M µj (r j ) w.r.