# Introduction to Riemannian Geometry and Geometric Statistics: from basic theory to implementation with Geomstats

Nicolas Guigui, Nina Miolane, Xavier Pennec

## ► To cite this version:

Nicolas Guigui, Nina Miolane, Xavier Pennec. Introduction to Riemannian Geometry and Geometric Statistics: from basic theory to implementation with Geomstats. Know Publishers, 16 (3), pp.182, 2023, Foundations and Trends® in Machine Learning, 978-1-63828-155-9. 10.1561/2200000098. hal-03766900v2

HAL Id: hal-03766900

https://inria.hal.science/hal-03766900v2

Submitted on 2 Feb 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Introduction to Riemannian Geometry and Geometric Statistics: From
Basic Theory to Implementation with Geomstats

Nicolas Guigui

Université Côte d'Azur, Inria

nicolas.guigui@inria.fr

Nina Miolane

University of California Santa Barbara

ninamiolane@ucsb.edu

Xavier Pennec

Université Côte d'Azur, Inria

xavier.pennec@inria.fr

February 2, 2023

Abstract

As data is a predominant resource in applications, Riemannian geometry is a
natural framework to model and unify complex nonlinear sources of data. However,
the development of computational tools from the basic theory of Riemannian geometry
is laborious. The work presented here forms one of the main contributions to the
open-source project geomstats, that consists of a Python package providing efficient
implementations of the concepts of Riemannian geometry and geometric statistics, both
for mathematicians and for applied scientists for whom most of the difficulties are
hidden under high-level functions. The goal of this monograph is two-fold. First, we
aim at giving a self-contained exposition of the basic concepts of Riemannian geometry,
providing illustrations and examples at each step and adopting a computational point
of view. The second goal is to demonstrate how these concepts are implemented in
Geomstats, explaining the choices that were made and the conventions chosen. The
general concepts are exposed and specific examples are detailed along the text. The
culmination of this implementation is to be able to perform statistics and machine
learning on manifolds, with as few lines of codes as in the wide-spread machine learning
tool scikit-learn. We exemplify this with an introduction to geometric statistics.

# 1 Introduction

Since the formal axiomization of Euclid in his famous *Elements* (dated around 300 BC),
geometry was considered as the properties of figures in the plane or in space. The abstract
notion of space as a mathematical object emerged in 1827 with C. F. Gauss' *Theorema
Egerium* proving that curvature is an intrinsic quantity of a surface, i.e. that can be
computed without reference to a "larger" embedding space. This notion was made precise by
the cornerstone work of Riemann (1868)<fnref n="1" /> built around the intuitive idea that a mathematical space results from varying a number of independent quantities, later identified as coordinates and formalized in the definition of a manifold by Whitney (1936). Riemannian Geometry (RG) is the study of such differentiable manifolds equipped with an inner product at each point that smoothly varies between points. This allows us to generalize the notions of angles, length and volumes, which can be integrated to global quantities highly coupled with the topology of the space.

Fruitful developments of these ideas allowed unifying previous examples of non-Euclidean geometries, that violate Euclid's parallel postulate (given a point and a straight-line, one and only one parallel straight-line can be drawn through the point). These ideas echoed with the developments of Lagrangian and Hamiltonian mechanics, and were instrumental in formalizing the modern theories of Physics, and especially Einstein's general relativity. They also made profound impact on many areas of mathematics such as group theory, representation theory, analysis, and algebraic and differential topology.

At the intersection of Physics and geometry, groups represent symmetries and transformations between states, and from the modern point of view of Klein's Erlangen program, the study of geometry boils down to studying the action of groups on a space, and their invariants. The work of Elie Cartan enabled significant progress in this direction.

Riemannian geometry has thus become a vast subject that is not usually taught before graduate education in mathematics or physics, and that requires familiarity with many concepts from differential geometry. Hence, although some books on the topic cover most of the pre-requisites and fundamental results of Riemannian geometry, the entry cost for applied mathematicians, computer scientists and engineers is high.

Nowadays, as data is a predominant resource in applications, Riemannian geometry is a natural framework to model and unify complex nonlinear sources of data. However, the development of computational tools from the basic theory of Riemannian geometry is laborous due to often high dimensional and non-exhaustive coordinate systems, nonlinear and intractable differential equations, etc. This monograph aims at providing the computational tools to perform statistics and machine learning on Riemannian manifolds to the wider data science community. The work presented here forms one of the main contributions to the open-source project `geomstats`, that consists in a Python package providing efficient implementations of the concepts of Riemannian geometry and geometric statistics, both for mathematicians and for applied scientists for whom most of the difficulties are hidden under high-level functions.

Other Python packages do exist and mainly fall under one of two following categories. On the one hand, there are the packages that focus on a single application, for instance on optimization: `Pymanopt` (Townsend et al., 2016), `Geoopt` (Becigneul and Ganea, 2019;

<footnote n="1">¹See Riemann (1873), the English translation by W.K. Clifford</footnote>
Kochurov et al., 2020), TensorFlow Riemann (Smirnov, 2021), and McTorch (Meghwanshi et al., 2018) or on deep learning such as PyG (Fey and Lenssen, 2019), where, in this case the geometry is often restricted to graph and mesh spaces. On the other hand, there are packages dedicated to a single manifold: PyRiemann on SPD matrices (Barachant, 2015), PyQuaternions on 3D rotations (Wynn, 2014), and PyGeometry on spheres, tori, 3D rotations and translations (Censi, 2012). Some other packages, like TheanoGeometry (Kühnel and Sommer, 2017) are not actively maintained anymore. There is therefore a need for a unified open-source implementation of differential geometry and associated learning algorithms for manifold-valued data.

The goal of this monograph is two-fold. First, we aim at giving a self-contained exposition of the basic concepts of Riemannian geometry, providing illustrations and examples at each step and adopting a computational point of view. We cover the basics of differentiable manifolds (Section 2), Riemannian manifolds (Section 3) and Lie groups (Section 4). Then we delve into more complex structures defined by invariance properties, in particular quotient metrics, and metrics on homogeneous and symmetric spaces (Section 5). Most proofs are omitted for brevity, but references to the proof of each statement are given. The interested reader may refer to the textbooks Lafontaine et al. (2004), Gallier and Quaintance (2020), Boumal (2023), and Lee (2003) for more details. Some mathematical definitions from the prerequisites can be found in the lexicon in Appendix A. The second goal is to demonstrate how these concepts are implemented in geomstats, explaining the choices that were made and the conventions chosen. The general concepts are exposed in Section 2.2, and detailed along the text and examples. The culmination of this implementation is to be able to perform statistics and machine learning on manifolds, with as few lines of codes as in the wide-spread machine learning tool scikit-learn. We exemplify this in Section 6 with a brief introduction to geometric statistics.

**Updates** This monograph was written for geomstats version 2.5.0 released in April 2022. There are two categories of code snippets in this monograph: those from the core of the package, that explain how it is implemented, and those for examples and figures. For updates on the former, we invite the interested reader to check the current master branch of the main geomstats repository<fnref n="2" />. In order to ensure compatibility of the latter with future releases of the package, we maintain a companion repository on Github<fnref n="3" />.

<footnote n="2">²https://github.com/geomstats/geomstats</footnote>

<footnote n="3">³https://github.com/geomstats/ftmal-paper</footnote>
# 2 Differentiable manifolds

## 2.1 Differentiable manifolds and tangent spaces

The differentiable manifold will be the structure underlying this entire monograph, yet its definition remains difficult for newcomers in the field. We start by that of an embedded manifold and generalize to the abstract case in a second part. The intuition behind the notion of smooth manifold is that around every point, it resembles the $d$-dimensional vector space $\mathbb{R}^d$ for some integer $d$, and the properties of $\mathbb{R}^d$ allow to define the notions of smooth functions, tangent vectors, etc. on the manifold.

### 2.1.1 Embedded manifolds

Let $N$ be a strictly positive integer. The fundamental examples of an embedded manifold are that of an open set* of $\mathbb{R}^N$, and a vector subspace $\mathbb{R}^d \times \{0\}_{N-d} \subset \mathbb{R}^N$ for $d \le N$, written $\mathbb{R}^d \times 0$ for short. These are ‘deformed’ via local diffeomorphisms* to obtain an embedded manifold, which can be thought of as a smooth surface in the ambient space. This was in fact one of the first motivations of the mathematical developments underlying the notion of manifold.

**Definition 2.1.1.** Let $d, N$ be integers with $1 \le d \le N$. Then a $d$-dimensional smooth embedded manifold in $\mathbb{R}^N$ is a non-empty subset $M$ of $\mathbb{R}^N$ such that for every point $p \in M$, there are two open subsets $U, V \subset \mathbb{R}^N$ with $p \in U$ and $0 \in V$, and a smooth diffeomorphism* $\varphi : U \to V$ such that $\varphi(U \cap M) = V \cap (\mathbb{R}^d \times \{0\}_{N-d})$.

In this definition, $\varphi$ may be called a local *chart*. Thankfully there are equivalent definitions that give greater insights into what makes $M$ a differentiable manifoldWe first need to define the notions of immersions and submersions.

**Definition 2.1.2.** Let $n \le p$ be two strictly positive integers, and $U \in \mathbb{R}^p, V \in \mathbb{R}^n$ two open sets.

* A differentiable* map $f: U \to V$ is called a *submersion* at $x \in U$ if $df_x$ is surjective*. We say that $f$ is a submersion if it is a submersion at every $x \in U$.

* A differentiable map $f: V \to U$ is called an *immersion* at $x \in V$ if $df_x$ is injective*. We say that $f$ is an immersion if it is an immersion at every $x \in V$.

The fundamental example of an immersion is the injection defined by $(x_1, \dots, x_n) \mapsto (x_1, \dots, x_n, 0, \dots, 0) \in \mathbb{R}^p$, while the projection $(x_1, \dots, x_p) \in \mathbb{R}^p \mapsto (x_1, \dots, x_n) \in \mathbb{R}^n$ is that of a submersion. In fact, one can show that up to a local change of variables (i.e. composition with a diffeomorphism), these maps are respectively the only immersions and

<footnote>*Defined in Appendix A.</footnote>
submersions. This results from the local inversion theorem, see Balzin (2021, Theorem 4.7) for a proof. We now have the following characterisation theorem (Gallier and Quaintance, 2020, Theorem 3.6):

**Theorem 2.1.1.** A nonempty subset $M \subseteq \mathbb{R}^N$ is a $d$-dimensional manifold if and only if any of the following conditions hold:

(1) (Local parametrization) For every $p \in M$, there are two open subsets $V \subseteq \mathbb{R}^d$ and $U \subseteq \mathbb{R}^N$ with $p \in U$ and $0 \in V$, and a smooth function $f : V \to \mathbb{R}^N$ such that $f(0) = p$, $f$ is a homeomorphism* between $V$ and $U \cap M$, and $f$ is an immersion at 0.

(2) (Local implicit function) For every $p \in M$, there exist an open set $U \in \mathbb{R}^N$ and a smooth map $f : U \to \mathbb{R}^{N-d}$ that is a submersion at $p$, such that $U \cap M = f^{-1}(\{0\})^2$.

(3) (Local graph) For every $x \in M$, there exist an open neighborhood $U \subseteq \mathbb{R}^N$ of $x$, a neighborhood $V \subseteq \mathbb{R}^d$ of 0 and a smooth map $f : V \to \mathbb{R}^{N-d}$ such that $U \cap M = \text{graph}(f)^3$.

The characterization (2) encodes the notion of constraint: a manifold can be understood as the set of points that verify a constraint defined by an implicit equation, given by the function $f$. Figure 2.1 (Left) gives an example where $f : x \in \mathbb{R}^2 \mapsto \|x\| \in \mathbb{R}$. This is one of the reasons manifolds are ubiquitous in applications, we will give many examples of this case. The other characterizations can also be understood as follows. The first (1) implies that at every point of the manifold, a coordinate system defined on $\mathbb{R}^d$ exists to parametrize the manifold around that point. This is depicted in Figure 2.1 (Right). Finally (3) is the most common to think of surfaces in $\mathbb{R}^3$ as sets of points $(x, y, f(x, y))$.

### Example 2.1: Hypersphere

The most simple manifold we will study is the hypersphere, or $d$-dimensional sphere. It is the set of unit-norm vectors of $\mathbb{R}^{d+1}$:

$$S^d = \{x \in \mathbb{R}^{d+1} \mid \|x\|_2^2 = 1\}.$$

Let $f : x \in \mathbb{R}^{d+1} \mapsto \|x\|^2 - 1 \in \mathbb{R}$. For $x \neq 0$, $df_x y = 2x^\top y$ is surjective, so Assertion (2) of Theorem 2.1.1 applies (with $U = \mathbb{R}^{d+1}$) and $S^d$ is a $d$-dimensional embedded manifold in $\mathbb{R}^{d+1}$. In dimension $d = 1$, this corresponds to the circle, and for $d = 2$ this is the usual sphere. Both cases are common to represent angles and directions in space, and as such appear in the field of directional statistics (Mardia and Jupp, 2009).

<footnote n="2">²Recall that the pre-image of a set $A$ by $f : E \to F$ is defined by $f^{-1}(A) = \{x \in E, f(x) \in A\}$.</footnote>

<footnote n="3">³The graph of $f$ is the set $\text{graph}(f) = \{(x, f(x)) | x \in V\}$.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 2.1:** Left: Representation of a manifold, definition with local $\varphi$ or by a local immersion $f$. Right: Representation of a manifold defined by a submersion, in this case the distance to the origin.</figcaption></figure>

### Example 2.2: Hyperbolic space

The fundamental counterpart to the hypersphere is the two-sheeted hyperboloid, defined by

$$H^d = \{x \in \mathbb{R}^{d+1} \mid -x_0^2 + \sum_{i=1}^{d} x_i^2 = -1\}.$$

It is one of the models of hyperbolic geometry, which is increasingly used to model hierarchical data, e.g. (Nickel and Kiela, 2017).

### Example 2.3: Special Orthogonal group

Matrix groups play an essential role in the theory of RG, and especially the special orthogonal group $SO(n)$, i.e. the set of unit determinant orthogonal matrices:

$$SO(n) = \{R \in M_n(\mathbb{R}) \mid R^\top R = I_n, \det(R) = 1\}.$$

Consider the map $f : \begin{cases} GL^+(n) & \longrightarrow & S(n) \\ A & \longmapsto & A^\top A - I_n \end{cases}$ where $GL^+(n) \in M_n(\mathbb{R}) \simeq \mathbb{R}^{n^2}$ is the open set of invertible squared matrices with positive determinant, $S(n)$ is the set of symmetric matrices of size $n$, a vector subspace of $M_n(\mathbb{R})$ of dimension $\frac{n(n+1)}{2}$. It is straightforward to show that the differential of $f$ at some $R$ is

$$df_R(H) = R^\top H + H^\top R,$$

and we can see that it is surjective for all $R \in SO(n)$, as for any $S \in S(n)$, $df_R(\frac{RS}{2}) = S$. As $SO(n) = f^{-1}(0)$, we conclude that $SO(n)$ is indeed an embedded manifold of dimension $\frac{n(n-1)}{2}$.
One can represent the motion of a rigid-body in the referential of its barycenter as a curve with values in SO(3), hence this group is widely used in e.g. robotics (Barczyk et al., 2015).

### Example 2.4: Stiefel manifold

A generalization of both hypersphere and special orthogonal group is the Stiefel manifold, defined as the set of orthonormal k-frames of $\mathbb{R}^n$. If we represent each vector $u_i$ of a k-frame $(u_1, \dots, u_k)$ as the $i^{th}$ column of a matrix $U$ (in the canonical basis of $\mathbb{R}^n$), then the Stiefel manifold can be seen as a subset in $M_{n,k}(\mathbb{R})$:

$$St(k, n) = \{U \in M_{n,k}(\mathbb{R}) | U^\top U = I_k\}.$$

As in the previous example, we can consider the map $f: U \mapsto U^\top U - I$ on an open subset of $M_{n,k}(\mathbb{R})$ and show that it is a submersion such that $S(k, n) = f^{-1}(0)$ to conclude that $S(k, n)$ is an embedded manifold of dimension $nk - \frac{k(k+1)}{2}$.

The Stiefel manifold naturally arises as the optimization domain in many problems related to matrix decompositions, in linear algebra, statistics, machine learning, computer vision, etc. see Absil et al. (2010, and references therein).

## 2.1.2 Manifolds

For generality, we now define a manifold in a more abstract way, i.e. as a topological space *that is not a priori embedded in some $\mathbb{R}^N$. The idea is still that a manifold is a space that can be covered by open sets that each look like (i.e. are diffeomorphic to an open set of) the usual space $\mathbb{R}^d$. Of course one can verify that embedded manifolds are indeed manifolds with this more general definition, and in fact, Whitney (1936) proved that any manifold can be smoothly embedded in a larger space, showing that the two concepts are equivalent. We first motivate the need for a more general definition of manifold by the example of the Kendall size-and-shape space.

### Example 2.5: Kendall size-and-shape space

The underlying idea is that a shape is what is left after removing the effects of translation and rotation. We first define the set of k landmarks of $\mathbb{R}^m$ as the space of $m \times k$ matrices $M_{m,k}(\mathbb{R})$. For $x \in M_{m,k}(\mathbb{R})$, let $x_i$ denote the columns of $x$, i.e. points of $\mathbb{R}^m$ and let $\bar{x}$ be their barycenter. We remove the effects of translation by

<footnote>*Defined in Appendix A.</footnote>
considering the matrix with columns $x_i - \bar{x}$ instead of $x$. Let $M_{m,k}^*(\mathbb{R})$ be the set of such centered matrices.

Now, in order to remove the effects of rotations, we would like to identify the landmark configurations that only differ by a rotation of all the landmarks. This defines an equivalence relation $\sim$:

$$ \forall x, y \in M_{m,k}^*(\mathbb{R}), \quad x \sim y \iff \exists R \in SO(m), y = Rx. $$

A shape thus corresponds to an equivalence class $[x]$ of landmark configurations, and we can define the size-and-shape space as the quotient of the landmark space by the equivalence relation $\sim$ (or equivalently by $SO(n)$):

$$ S\Sigma_m^k = \{[x] \mid x \in M_{m,k}^*(\mathbb{R})\}. $$

A quotient space of a manifold by another manifold may not even be a Hausdorff space, but we will give sufficient conditions in Section 4.4 to ensure that quotients resulting from a group action are indeed smooth manifolds.

In this case, $S\Sigma_m^k$ inherits a differentiable structure from the landmark space, that turns it into a smooth manifold, although we cannot see it explicitly as a subset of $\mathbb{R}^N$ for some $N$. Implementing tools to work with data on such spaces is a challenging task that we tackle in geomstats. This is the subject of Section 5. Kendall size-and-shape spaces are a ubiquitous framework for the statistical analysis of data arising from medical imaging, computer vision, biology, chemistry and many more domains (Dryden and Mardia, 2016).

**Definition 2.1.3 (Atlas).** Let $M$ be a topological space* and $k \ge 1$ an integer. A $C^k$-atlas $\mathcal{A}$ with values in $\mathbb{R}^d$ is a collection of pairs $(U, \varphi)$ called charts where $\varphi: U \to V$ is a homeomorphism between the open sets $U \subset M$ and $V \subset \mathbb{R}^d$, such that $M \subset \bigcup_{U \in \mathcal{A}} U$ and for any $(U, \varphi)$ and $(U', \varphi')$ in $\mathcal{A}$, the transition map

$$ \varphi' \circ \varphi^{-1} : \varphi(U \cap U') \to \varphi'(U \cap U') $$

is a $C^k$-diffeomorphism*.

Two atlases are $C^k$-compatible if their union is still a $C^k$-atlas. Compatibility defines an equivalence relation, and we will think of the equivalence class of an atlas whenever referring to one. There is a unique maximal atlas (for the inclusion) that contains a given atlas.

Note that the transition maps are defined between open sets of $\mathbb{R}^d$ (see Figure 2.2), the usual notions of differentiability are thus available, and will allow to define such notions on

<footnote>*Defined in Appendix A.</footnote>
manifolds. We will always consider the case $k = \infty$ and say that $C^\infty$ maps and atlases are smooth.

<figure><img src="image_3.png" /><figcaption>Figure 2.2: Transition maps.</figcaption></figure>

**Definition 2.1.4 (Differentiable manifold).** We call *differentiable manifold* of class $C^k$ and dimension $d$ any topological space $M$ that is Hausdorff* and second-countable* together with a maximal $C^k$-atlas $\mathcal{A}$ with values in $\mathbb{R}^d$.

We sometimes refer to the atlas of $M$ as its *differentiable structure*, and to this definition of manifold as the *intrinsic definition* as the charts are defined on $M$ rather than an extrinsic embedding space. A chart $(U, \varphi)$ defines a set of *local coordinates* on $U$ written $(x^1, \dots, x^d)$ for short, and defined by $x^i = \text{pr}_i(\varphi(x))$, where $\text{pr}_i$ is the projection on the $i^{th}$ coordinate of $\mathbb{R}^d$. Let us now use the intrinsic definition of a manifold to exemplify further how manifolds can be obtained from others.

### Example 2.6: Product manifold

Let $M, N$ be two manifolds with $(U_i, \varphi_i)_{i \in I}$, $(V_j, \psi_j)_{j \in J}$ their respective atlas. Define for $(i, j) \in (I \times J)$

$$ \phi_{ij}: \begin{cases} U_i \times V_j & \to \varphi_i(U_i) \times \psi_j(V_j) \\ (x, y) & \mapsto (\varphi_i(x), \psi_j(y)) \end{cases} $$

Then it is easy to check that $(U_i \times V_j, \varphi_{ij})_{(i,j) \in I \times J}$ is an atlas for the product space $M \times N$. This atlas does not depend on the choice of original atlases of $M, N$ (in the right equivalence class), and allows to endow $M \times N$ with the structure of manifold. We call it the product manifold of $M$ and $N$. Its dimension is $\dim(M \times N) = \dim M + \dim N$.

To give insights into the importance of the general notion of manifold, let us now consider two counter-examples.
### Example 2.7: Cusp and Node

For more details on these two examples see Gallier and Quaintance (2020, chapter 7). First, we consider the classic example of a space that is not a manifold: the nodal cubic, shown in Figure 2.3. It is the set of points

$$M_1 = \{(x, y) \in \mathbb{R}^2 | y^2 = x^2 - x^3\},$$

considered under the subset topology. The self-intersection at the origin does not preserve the topology of $\mathbb{R}$, so no homeomorphism can exist between $M$ and $\mathbb{R}$ around the origin. Thus, $M$ is not a manifold.

<figure><img src="image_4.png" /><figcaption>Figure 2.3: Nodal cubic $M_1$</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 2.4: Cuspidal curve $M_2$</figcaption></figure>

Secondly, we consider the cuspidal curve displayed on Figure 2.4 and defined as the set

$$M_2 = \{(x, y) \in \mathbb{R}^2 | y^2 = x^3\}.$$

We can define the maps $\varphi: (x, y) \in M_2 \mapsto y^{1/3} \in \mathbb{R}$ and $\psi: (x, y) \in M_2 \mapsto y \in \mathbb{R}$, that each define a smooth atlas on $M_2$ and endow it with a differentiable structure of manifold. However, the two atlases (constituted of a single chart) are not compatible, so they define different manifolds.

For the next notions that will be introduced, we will use the convenient setting of embedded manifolds, but all these notions can be generalized to the abstract case by using charts to recover functions defined between vector spaces.

#### 2.1.3 Tangent spaces and differentiable maps

We first define smooth curves on manifolds, using a local parametrization of $M$ and the notion of smooth function from $\mathbb{R}$ to $\mathbb{R}^d$. See Figure 2.5 for a representation.

The following definition does not depend on the choice of local parametrization. See Figure 2.6 for a representation.

**Definition 2.1.5 (Smooth curve).** Let $M$ be a $d$-dimensional manifold in $\mathbb{R}^N$. A smooth
<figure><img src="image_6.png" /><figcaption>Figure 2.5: Definition of a smooth curve on a manifold.</figcaption></figure>

curve $\gamma$ in $M$ is any function $\gamma: I \to M$ where $I$ is an open interval in $\mathbb{R}$, such that for any $t \in I$, $p = \gamma(t)$, there is a local chart $\varphi: U \to V$ of $M$ at $p$ and $\epsilon > 0$ such that $\varphi \circ \gamma: (t - \epsilon, t + \epsilon) \to \mathbb{R}^d$ is smooth.

This definition is extended to smooth curves defined on a closed interval $I = [a, b]$ of $\mathbb{R}$ by requiring that $\gamma$ be the restriction of some smooth curve defined on an open interval that contains $[a, b]$. As $\gamma: I \to M \subset \mathbb{R}^N$ is a curve in $\mathbb{R}^N$ and is differentiable, a tangent vector along $\gamma$ at some time $t \in I$ is obviously defined. This generalizes to tangent spaces to the manifold.

<figure><img src="image_7.png" /><figcaption>Figure 2.6: Left: Definition of a tangent vector as the derivative of a curve. Right: The collection of all tangent vectors forms a vector space.</figcaption></figure>

**Definition 2.1.6 (Tangent vector).** Let $d \le N \in \mathbb{N}$, $M$ be an embedded manifold in $\mathbb{R}^N$ of dimension $d$ and $p \in M$. A vector $v \in \mathbb{R}^N$ is tangent to $M$ at $p$ if there exists an open interval $I$ centered around 0, and a curve $\gamma: I \to M$ such that

$$ \gamma(0) = p \text{ and } \dot{\gamma}(0) = v. $$

We write $T_pM$ for the set of tangent vectors at $p$.

Recall that $\mathbb{R}^d \times \{0\}$ is a fundamental example of an embedded manifold, it is clear that the tangent space at any $p \in \mathbb{R}^d \times \{0\}$ is the whole $\mathbb{R}^d \times \{0\}$. From this case we
deduce the characterizations of tangent spaces equivalent to that of manifolds obtained in Theorem 2.1.1 (Paulin, 2007, Proposition 3.1).

**Theorem 2.1.2.** Let $M$ be a manifold in $\mathbb{R}^N$ of dimension $d$.

(1) If $U, V \subset \mathbb{R}^N$ are two open neighborhoods respectively around $p$ and $0$ in $\mathbb{R}^N$ and $f: U \to V$ is a diffeomorphism such that $f(p) = 0$ and $f(U \cap M) = V \cap (\mathbb{R}^d \times \{0\})$ then

$$T_p M = df_p^{-1}(\mathbb{R}^d \times \{0\}).$$

(2) (Local parametrization) If $U \subset \mathbb{R}^N$ is an open neighborhood of $p \in M$, $V \subset \mathbb{R}^d$ is an open neighborhood around $0$ and $f: V \to \mathbb{R}^N$ is a smooth function such that $f(0) = p$, $f$ is a homeomorphism between $V$ and $U \cap M$, and $f$ is an immersion at $0$, then<fnref n="2" />

$$T_p M = \mathrm{Im} \, df_0.$$

(3) (Local implicit function) If $U \subset \mathbb{R}^N$ is an open neighborhood around $p \in M$ and $f: U \to \mathbb{R}^{N-d}$ is a smooth map that is a submersion at $p$, such that $U \cap M = f^{-1}(\{0\})$, then<fnref n="3" />

$$T_p M = \ker df_p.$$

(4) (Local graph) If $U \subset \mathbb{R}^N$ is an open neighborhood of $p \in M$, $V \subset \mathbb{R}^d$ a neighborhood of $0$ and $f: V \to \mathbb{R}^{N-d}$ a smooth map such that $U \cap M = \mathrm{graph}(f)$ and $p = (0, f(0))$, then

$$T_p M = \mathrm{Im}\{v \mapsto (v, df_0(v))\}.$$

From (1) we can see that $T_p M$ is a linear subspace of $\mathbb{R}^N$ of dimension $d$. Tangent spaces thus provide local linearizations of the manifold, a property that will be useful as a first way to handle data on manifolds. The previous theorem allows computing the tangent spaces of the common manifolds seen in the previous section.

### Example 2.8: Tangent space of the hypersphere

Recall that the hypersphere is the embedded manifold defined by $S^d = f^{-1}(0)$ where $f: x \mapsto \|x\|^2 - 1$. This corresponds to (3) of Theorem 2.1.2, therefore for any $x \in S^d$

$$T_x S^d = \{v \in \mathbb{R}^{d+1} \mid \langle x, v \rangle = 0\}.$$

<footnote n="2">²Recall that the range of $f: E \to F$ is defined by $\mathrm{Im} \, f = \{f(x), x \in E\} \subset F$.</footnote>

<footnote n="3">³Similarly the kernel of $f: E \to F$ is defined by $\mathrm{ker} \, f = \{x \in E \mid f(x) = 0\} \subset E$.</footnote>
### Example 2.9: Tangent space of the hyperbolic space

Similarly, as the hyperbolic space is defined as $H^d = f^{-1}(0)$ where $f : x \mapsto -x_0^2 + \sum_{i=1}^d x_i^2 + 1$, we obtain for any $x \in H^d$

$$T_x H^d = \{v \in \mathbb{R}^{d+1} \mid -x_0 v_0 + \sum_{i=1}^d x_i v_i = 0\}.$$

### Example 2.10: Tangent space of Stiefel manifold

Similarly, for any $U \in St(k, n)$,

$$T_U St(k, n) = \{H \in M_{n,k}(\mathbb{R}) \mid U^\top H + H^\top U = 0\}.$$

### Example 2.11: Tangent space of SO(n)

Recall $SO(n) = f^{-1}(0)$ with $f: A \mapsto A^\top A - I_n$, and for any $R \in SO(n), H \in M_n(\mathbb{R})$ we have $df_R(H) = R^\top H + H^\top R$. Therefore for any $R \in SO(n)$

$$T_R SO(n) = \{H \in M_n(\mathbb{R}) \mid R^\top H + H^\top R = 0\}.$$

Note the special case $R = I_n$, then $T_{I_n} SO(n) = \text{Skew}(n)$, the set of skew-symmetric matrices of size $n$. The tangent space at the identity of a Lie group plays a particularly important role as will be exposed in Section 4.

We can now define the notion of smooth maps between manifolds and their differential.

**Definition 2.1.7 (Smooth map).** Let $M, Q$ be two manifolds of dimensions $d_1, d_2 \in \mathbb{N}$ embedded in $\mathbb{R}^N$. A function $f: M \to Q$ is smooth if for every $p \in M$, there are parametrizations $\varphi: V_1 \to U_1$ of $M$ at $p$ and $\psi: V_2 \to U_2$ of $Q$ at $f(p)$ such that $f(U_1) \subseteq U_2$ and

$$\psi^{-1} \circ f \circ \varphi: V_1 \to V_2 \quad \text{is smooth.}$$

Note that in the above definition $V_1, V_2$ are open sets respectively of $\mathbb{R}^{d_1}, \mathbb{R}^{d_2}$ so the notion of smooth map from $V_1$ to $V_2$ is already well known (see lexicon in Appendix A).

**Definition 2.1.8 (Differential).** Let $M, Q$ be two manifolds of dimensions $d_1, d_2 \in \mathbb{N}$ embedded in $\mathbb{R}^N$ and $f: M \to Q$ a smooth map. For any $p \in M$ and any $v \in T_p M$, let $\gamma$ be a smooth curve through $p$ such that $\dot{\gamma}(0) = v$ and define

$$df_p(v) = (f \circ \gamma)'(0).$$

This definition does not depend on the choice of curve $\gamma$ and the map $df_p: T_p M \to T_{f(p)} Q$ is called the differential or tangent map of $f$ at $p$. It is a linear map between tangent spaces.
It generalizes the differential map of a differentiable function defined from $\mathbb{R}^{N_1}$ to some $\mathbb{R}^{N_2}$ to functions defined on the manifold $M$ instead of the embedding space. It coincides with the original differential (Definition A.0.11 in the lexicon) when $M = \mathbb{R}^d$, hence the use of the same notation $df_p$. The set of real-valued smooth maps on $M$ will be particularly useful. For short, we denote it $C^\infty(M) \triangleq C^\infty(M, \mathbb{R})$. $C^\infty(M)$ is clearly an infinite dimensional $\mathbb{R}$-vector space, and with point-wise multiplication, an algebra.

Next, it is convenient to consider the set of all the tangent spaces at all points

$$TM = \bigsqcup_{x \in M} \{x\} \times T_x M = \{(x, v) | x \in M, v \in T_x M\}.$$

and its natural projection

$$\pi : \left\{ \begin{array}{l} TM \quad \to \quad M \\ (x, v) \quad \mapsto \quad x \end{array} \right. .$$

This space is called the *tangent bundle* of $M$, and one can show that if $M$ is a manifold of class $C^{k+1}$ and dimension $d$, then $TM$ is itself a manifold in $\mathbb{R}^N \times \mathbb{R}^N$, of class $C^k$ and dimension $2d$. The tangent bundle is the domain of definition of the differential of smooth functions defined on manifolds:

$$f : M \to Q, \quad df : TM \to TQ$$

It is also the space where *vector fields* are valued: a vector field $X$ is a smooth assignment of a tangent vector to each point of a manifold, i.e. $X : M \to TM$ such that $\forall p \in M, \pi \circ X(p) = p$. $X(p)$ will be written $X_p$ for convenience. Let $\Gamma(TM)$ denote the set of all vector fields (VF). It is clear that $\Gamma(TM)$ equipped with point-wise sum and multiplication by a scalar forms a vector space. Multiplication by a smooth function is also defined pointwise: for any $f \in C^\infty(M)$ and $X \in \Gamma(TM)$, $fX$ is the vector field such that

$$\forall p \in M, (fX)_p = f(p)X_p.$$

This turns the set of vector fields into a $C^\infty(M)$-module. For a smooth map $f : M \to \mathbb{R}$ and a vector field $X$, we write $X(f)$ the function defined at every $p$ by

$$X(f)(p) = df_p(X_p).$$

This leads to the following remark.

**Remark 2.1.1.** We defined vector fields as sections of the tangent bundle, i.e., maps $\sigma : M \to TM$ such that $\pi \circ \sigma = \text{Id}$. Alternatively, vector fields can be defined as *derivations* over the algebra $C^\infty(M)$ of smooth real valued functions. A derivation $X : C^\infty(M) \to C^\infty(M)$ is a linear map that satisfies the Leibniz rule:

$$\forall f, g \in C^\infty(M), \quad X(fg) = fX(g) + X(f)g. \qquad (2.1)$$
One can check that a vector field as defined above indeed defines a derivation. However, applying the “composition” of two vector fields to a function $f$ is not a derivation because of second-order derivatives of $f$. This leads to the definition of the Lie bracket of vector fields.

**Definition 2.1.9 (Lie bracket over $\Gamma(TM)$).** The Lie bracket of vector fields is defined as the map

$$ [\cdot, \cdot] : \left\{ \begin{array}{ccc} \Gamma(TM) \times \Gamma(TM) & \to & \Gamma(TM) \\ (X, Y) & \mapsto & f \mapsto X(Y(f)) - Y(X(f)) \end{array} \right. \qquad (2.2) $$

A useful tool to handle vector fields locally is to use a basis of $T_pM$ for $p$ in some open set $U$.

**Definition 2.1.10 (Frame).** Let $M$ be a $d$-dimensional manifold. For any open set $U \subseteq M$, a family of vector fields $(X_1, \dots, X_d)$ over $U$ is called a *frame* over $U$ if for every $p$ in $U$, $(X_1(p), \dots, X_d(p))$ is a basis of $T_pM$.

Any chart $(U, \varphi)$ defines a local frame that corresponds to its local coordinates: define the $i^{th}$ curve $\gamma_i : t \in \mathbb{R} \mapsto \varphi^{-1}(0, \dots, 0, t, 0, \dots, 0) \in U$ and $X_i(\gamma_i(t)) = \dot{\gamma}_i(t)$. The vector field $X_i$ defined on $U$ is often written $X_i = \frac{\partial}{\partial x^i}$ or simply $\partial_i$ and corresponds to $d\varphi^{-1}(e_i)$ where $(e_1, \dots, e_d)$ is the canonical basis of $\mathbb{R}^d$. Then the family $(X_1, \dots, X_d)$ is a local frame over $U$.

**Remark 2.1.2.** If a family $(X_1, \dots, X_d)$ is a frame over the whole manifold (i.e. $U = M$), we say that it is a global frame. Whether global frames exist depends on the topology of the manifold, and in that case the tangent bundle is called trivial, i.e. isomorphic to the direct product $M \times \mathbb{R}^N$. This is not the case of e.g. the sphere (of dimension 2), by the hairy ball theorem<fnref n="4" />.

Vector fields can be considered as infinitesimal generators of local maps, as we shall see in the following. These maps, called flows, usually supply strong information on global properties of the manifold. In this monograph, we will mainly focus on geodesic flows (Section 3.1), and flows of left-invariant vector fields on Lie groups (Section 4.2).

**Definition 2.1.11 (Integral curve).** Let $X \in \Gamma(M)$ and $p_0 \in M$. An *integral curve* for $X$ with initial condition $p_0$ is a curve $\gamma : I \to M$ such that

$$ \forall t \in I, \quad \dot{\gamma}(t) = X_{\gamma(t)} \quad \text{and} \quad \gamma(0) = p_0, $$

where $I = (a, b) \subseteq \mathbb{R}$ is an open interval containing 0.

An integral curve is thus a curve whose speed $\dot{\gamma}(t)$ coincide with $X$ at any point along the curve (see Figure 2.7). A collection of such integral curves is called a flow (Figure 2.8):

<footnote n="4">⁴https://en.wikipedia.org/wiki/Hairy\_ball\_theorem</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 2.7:** Example of an integral curve. The curve $\gamma$ is obtained by integrating the vector field $X$ (blue), meaning that $X_{\gamma(t)}$ is tangent to $\gamma$ at all times.</figcaption></figure>

**Definition 2.1.12 (Local flow).** Let $X \in \Gamma(M)$ and $p_0 \in M$. A local flow of $X$ at $p_0$ is a map $\phi : I \times U \to M$ where $I$ is an open interval containing 0 and $U$ is an open set in $M$ containing $p$, such that for every $p \in U$, the curve $t \mapsto \phi(t,p)$ is an integral curve of $X$ starting from $p$.

<figure><img src="image_9.png" /><figcaption>**Figure 2.8:** Example of a flow. The point $\phi(t,p)$ is reached at time $t$ by the integral curve of the vector field $X$.</figcaption></figure>

Thanks to the theory of ordinary differential equations (ODE), one can prove that for any vector field, there is a local flow defined around any point, and if two such flows are defined on overlapping domains, they coincide on the intersection. For $t \in I$, we write $\phi_t : x \mapsto \phi(t, x)$. It is clear that $\phi_0$ is the identity map, and for some $t \neq 0$, $\phi_t$ is a map defined locally on $M$. See Lafontaine et al. (2004, Proposition 1.55-56-58) for proofs.

**Proposition 2.1.1.** Let $X$ be a smooth vector field on $M$, $p_0 \in M$ and $\phi : I \times U \to M$ the local flow of $X$ at $p_0$. For any $s, t \in I$ and $x \in U$,

* if $\phi_s(x) \in U$ and $t+s \in I$, then $\phi_t \circ \phi_s(x) = \phi_{t+s}(x)$;

* $\phi_t$ is a local diffeomorphism;

* $\phi_t$ preserves $X$, in the sense that $\forall t \in I, \forall x \in U, d(\phi_t)_x(X_x) = X_{\phi_t(x)}$.

**Definition 2.1.13 (Complete).** We say that $X$ is *complete* if the domain of definition of its flow $(t, x) \mapsto \phi_t(x)$ is the entire $\mathbb{R} \times M$.

In that case $\phi_t$ is a diffeomorphism of $M$, and $(\phi_t)_{t \in \mathbb{R}}$ is a *one-parameter subgroup* of $\text{Diff}(M)$.

We now have the ingredients to introduce *geomstats*.
## 2.2 Implementation in geomstats

Now that the fundamental notion of manifold has been exposed, we can delve more into the architecture of the `geomstats` package, and summarize the choices that we made in its development. Firstly the package is organized in different modules that distinguish between the geometric and the statistical operations. There is thus a geometry module, that gathers all the implementations of manifolds, connections and Riemannian metrics, and a learning module where estimation algorithms are implemented in a generic fashion and take the geometric structure and the data as inputs. The goal is that all the learning algorithms can be run seamlessly on different manifolds, and with different metrics. Basic sampling schemes are available, and a more extended sampling module is currently being developed with the same spirit. A visualization module allows to plot data on the common manifolds in dimension two or three. In this section, we focus on the geometry module, and more specifically on the objects that represent manifolds. The Riemannian metric objects will be described in Section 3.1 along with the mathematical definitions. The statistical and learning tools will be described in Section 6.

The package is object-oriented in the sense that all the tools are implemented as classes that contain all the methods related to a tool. Object-Oriented programming (OOP) is a programming paradigm that consists in grouping properties and functions related to a common concept into an object, called a *class* in Python<fnref n="5" />

The geometry module of `geomstats` is organized by classes that each represent a geometric structure. To guarantee the consistency of all the classes, we implement an abstract parent `Manifold` class, and the actual implementations of the usual manifolds are subclasses of this parent class. The aim of abstract base classes is to provide the minimal skeleton of attributes and methods expected in its subclasses. The methods of the abstract class are thus declared but contain no implementation, and they are overridden by the subclasses.

**The manifold classes** A subclass of `Manifold` gathers the methods to work with data lying on the considered manifold. Note that these classes do not explicitly provide a representation of the manifold (as e.g. a triangulated surface in 2d), but the tools to handle points and tangent vectors. As we work with embedded manifolds in most cases, points and tangent vectors are themselves represented by multi-dimensional arrays. These are either `NumPy` arrays, or `TensorFlow` of `PyTorch` tensors according to the backend that is being used.

Mathematically, the first attribute of a manifold is of course its dimension, called `dim` for brevity throughout the package. Then we use an attribute to inform on the expected type of the point: whether vectors (for e.g. the hypersphere and hyperbolic space) or

<footnote n="5">⁵For more on OOP, we refer the reader to online tutorials such as the one written by the RealPython team.</footnote>
matrices (e.g. SPD matrices, the special orthogonal group) should be used. This is called `default_point_type`.

Furthermore, a `Manifold` in `geomstats` should always implement a method that evaluates whether a given element belongs to that manifold, and whether a given vector is a tangent vector at a given point. These are the `belongs` and `is_tangent` methods. For practical reasons we also add a `random_point` method, to generate random points that belong to the manifold (regardless of the distribution). This is useful in particular to test the methods and the learning algorithms. We obtain the following base class:

```python
class Manifold(abc.ABC):
    """Class for manifolds."""

    def __init__(
        self, dim, metric=None, default_point_type='vector', **kwargs):
        super().__init__(**kwargs)
        self.dim = dim
        self.default_point_type = default_point_type
        selfmetric = metric

    @abc.AbstractMethod
    def belongs(self, point, atol=gs.atol):
        """Evaluate if a point belongs to the manifold."""

    @abc.AbstractMethod
    def is_tangent(self, vector, base_point, atol=gs.atol):
        """Check whether the vector is tangent at base_point."""

    @abc.AbstractMethod
    def random_point(self, n Samples=1, bound=1.):
        """Sample random points on the manifold."""
```

Note that the `Manifold` class contains a metric attribute. This will be detailed in Section 3.1.

**Implementation trick 2.1.** The methods decorated (`@` symbol) with `abc.AbstractMethod` are declared as abstract methods of the class. A class that contains abstract methods cannot be instantiated. This constrains the developer to implement these functions explicitly when writing subclasses of `Manifold`.

**Two elementary classes of manifolds** Throughout the current section, we have met two elementary ways of defining a manifold, that correspond respectively to (1) and (2) of Theorem 2.1.1 (page 5):

1. As pre-image of a submersion $f: \mathbb{R}^N \to \mathbb{R}^{N-d}$. We refer to such space as a level-set. In this case, specifying $f$, it is straightforward to implement the `belongs` and `is_tangent` method by evaluating $f$ and its differential. It also makes sense to add projection and to_tangent methods from the embedding space to respectively the manifold and the tangent space at a point.
2. As open sets of a $d$-dimensional vector space, called ambient space. In this case, all the tangent spaces are identified with the ambient space. For consistency, we add a projection method that maps any $d$-dimensional vector to the manifold at a tolerance threshold $\epsilon > 0$ away from the boundary of the set (if there is one). This method is not uniquely defined and can be understood as a regularization method for inputs very close to the open set, which is helpful in learning algorithms. Therefore, the method `is_tangent` just checks if the input belongs to the ambient space, and we can add a method `to_tangent`, that calls the projection of the ambient space, to project to a tangent space, assuming the ambient space to be itself embedded in another space, or being a vector space, whose projection is the identity.

We thus implement two more abstract classes, for open sets and level sets respectively:

```python
class OpenSet(Manifold, abc.ABC):
    """Class for manifolds that are open sets of a vector space."""

    def __init__(self, dim, ambient_space, **kwargs):
        kwargssetdefault("shape", ambient_space.shape)
        super().__init__(dim=dim, **kwargs)
        self.ambient_space = ambient_space

    def is_tangent(self, vector, base_point, atol=gs.atol):
        """Check whether the vector is tangent at base_point."""
        return self.ambient_space.b belongs (vector, atol)

    def to_tangent(self, vector, base_point):
        """Project a vector to a tangent space of the manifold."""
        return self.ambient_space.projection (vector)

    def random_point(self, n Samples=1, bound=1.):
        """Sample random points on the manifold."""
        sample = self.ambient_space.random_point(n Samples, bound)
        return self.projection (sample)

@abc .abstractmethod
    def projection(self, point):
        """Project a point in ambient manifold on manifold."""
```

```python
class LevelSet(Manifold, abc.ABC):
    """Class for manifolds embedded in a vector space by a submersion."""

    def __init__(
        self, dim, embedding_space, submersion, value,
        tangent_submersion, **kwargs):
        super().__init__ (
            dim=dim, default_point_type=embedding_space.default_point_type,
            **kwargs)
        self embedding_space = embedding_space
        self embedding metric = embedding_space .metric
        self.submersion = submersion
        self.value = value
        self.tangent_submersion = tangent_submersion

    def belongs(self, point, atol=gs.atol):
        """Evaluate if a point belongs to the manifold."""
```
```
belongs = self Embedding_space.belongs(point, atol)
if not gs.any(belongs):
    return belongs
value = self.value
constraint = gs.isclose(self.submersion(point), value, atol=atol)
if value.ndim == 2:
    constraint = gs.all(constraint, axis=(-2, -1))
elif value.ndim == 1:
    constraint = gs.all(constraint, axis=-1)
return gs(logical_and(belongs, constraint))

def is_tangent(self, vector, base_point, atol=gs.atol):
    """Check whether the vector is tangent at base_point."""
    belongs = self embedding_space.belongs(vector, atol)
    tangent_sub_applied = self.tangent_submersion(vector, base_point)
    constraint = gs.isclose(tangent_sub_applied, 0., atol=atol)
    value = self.value
    if value.ndim == 2:
        constraint = gs.all(constraint, axis=(-2, -1))
    elif value.ndim == 1:
        constraint = gs.all(constraint, axis=-1)
    return gs(logical_and(belongs, constraint))

@abcimmethod
def projection(self, point):
    """Project a point in embedding space on the manifold."""
```

**Implementation trick 2.2.** All the methods of `geomstats` are vectorized<fnref n="6" />, in the sense that they can take as argument either one input, or a collection of inputs corresponding to multiple samples. It is very useful to use the `einsum` method for that purpose, with ellipses (`'...'`) that represent an optional additional dimension.

For example, the syntax `gs.einsum('...i->...i', coef, point)` performs scalar multiplication between a list of scalars (`coef`) and a list of points, but it also works for a single scalar and a single point.

To make sure that the attributes that represent the ambient or embedding space do implement the methods that are called in `OpenSet` and `LevelSet`, we also implemented an abstract `VectorSpace` class. Actual manifolds are then implemented as subclasses of the corresponding abstract manifold and must implement all the abstract methods. We give examples of each class below, using the Hypersphere and Stiefel manifold as level sets, and the Poincaré ball as open set.

### Example 2.12: Implementation of the hypersphere

The hypersphere is implemented as embedded in $\mathbb{R}^{d+1}$, so it is a subclass of `LevelSet`.

<footnote n="6">⁶Vectorization may also be referred to as *array programming* on e.g. Wikipedia.</footnote>
class Hypersphere (LevelSet):
    """Class for the n-dimensional hypersphere."""
    def __init__ (self, dim):
        super() .__init__ (
            dim=dim, embedding_space=Euclidean(dim + 1),
            metric=HypersphereMetric(dim),
            submersion= lambda x: gs .sum(x ** 2, axis=-1), value=1.,
            tangent_submersion= lambda v, x: 2 * gs .sum(x * v, axis=-1))
    def projection(self, point):
        """Project a point on the hypersphere."""
        norm = gs .lINALG .norm(point, axis=-1)
        if gs .any(norm < gs .ATOL ):
            logging.warning('0 cannot be projected to the hypersphere')
        return gs .einsum('... ,...i->...i', 1. / norm, point)
    def to_tangent(self, vector, base_point):
        """Project a vector to the tangent space."""
        sq_norm = gs .sum(base_point ** 2, axis=-1)
        innerProd = self .embeddingMetric .inner_product(base_point, vector)
        coef = innerProd / sq_norm
        return vector - gs .einsum('... ,...j->...j', coef, base_point)

Example 2.13: Implementation of the Stiefel manifold

The Stiefel manifold is implemented as embedded in the space of *n* × *p* matrices,
so it is a subclass of LevelSet. The derivation of the projection map can be found
in Absil and Malick (2012).
```python
class Stiefel(LevelSet):
    """Class for Stiefel manifolds St(n,p."""

    def __init__(self, n, p, **kwargs):
        if p > n:
            raise ValueError("p needs to be smaller than n.")

        dim = int(p * n - (p * (p + 1) / 2))
        matrices = Matrices(n, p)
        canonicalapa metric = StiefelCanonicalMetric(n, p)
        kwargsOGdefault("metric", canonical_apa metric)
        super(Stiefel, self).__init__(
            dim=dim,
            embedding_space=matrices,
            submersion=lambda x: matrices.mul matrices transpose(x), x),
            value=gs. eye(p),
            tangent_submersion=lambda v, x: 2
            * matrices.to_symmetric matrices.mul matrices transpose(x), v)),
            **kwargs
        )
        self.n = n
        self.p = p

    def to_tangent(self, vector, base_point):
        """Project a vector to a tangent space of the manifold."""
        aux = Matrices.mul Matrices transpose(base_point), vector
        sym_aux = Matrices.to_symmetric (aux)
        return vector - Matrices.mul (base_point, sym_aux)

    def projection(self, point):
        """Project a close enough matrix to the Stiefel manifold."""
        mat_u, ..., mat_v = gs.linalg.svd(point)
        return Matrices.mul (mat_u..., : , : self.p], mat_v)
```

## Example 2.14: Implementation of the Poincaré ball

The Poincare ball is one of the models of hyperbolic geometry, and is defined as the open unit disk of $\mathbb{R}^d$. It is then a subclass of `OpenSet`.

```python
class PoincareBall(OpenSet):
    """Class for the n-dimensional hyperbolic space."""

    def __init__(self, dim, scale=1):
        super().__init__(
            dim=dim, ambient_space=Euclidean(dim),
            metric=PoincareBallMetric(dim))

    def belongs(self, point, atol=gs.atol):
        """Test if a point belongs to the unit ball."""
        return gs.sum(point**2, axis=-1) < (1 - atol)
```

As stated above, the `projection` method is not well defined and we simply map any $d$-dimensional vector to the manifold at a tolerance threshold $\epsilon > 0$ away from the boundary of the set.
```python
def projection(self, point):
    """Project a point on the unit ball."""
    if point.shape[-1] != self.dim:
        raise NameError("Wrong dimension, expected ", self.dim)

    l2_norm = gs.linalg.norm(point, axis=-1)
    if gs.any(l2_norm >= 1 - gs.atol):
        projected_point = gs.einsum(
            '...j,...->...j', point * (1 - gs.atol), 1. / l2_norm)
        projected_point = -gs.maximum(-projected_point, -point)
    return projected_point

return point
```

Finally, we illustrate this with a diagram representing all the base classes for manifolds and all the manifolds on Figure 2.9. The abstract base classes are shown with black bounding boxes.

<figure><img src="image_10.png" /><figcaption>**Figure 2.9:** Architecture of the manifolds of `geomstats`. The abstract base classes are in black bounded boxes. Inheritance is shown by blue arrows. `Hyperbolic` is an exception as it is a common interface to the different representations of hyperbolic geometry.</figcaption></figure>

Inheritances occur in two cases:

* With an abstract base class as parent class;

* When both parent and child class represent the same manifold, as is the case in the `CorrelationBundle` of Example 5.8. A common interface `Hyperbolic` for the three representations of hyperbolic geometry follows this logic as any of the three representations can be chosen, but the instantiated object is either `Hyperboloid`, `PoincareBall` or `PoincareHalfSpace` as chosen.
Manifolds can then be composed to define other manifolds by products (Example 2.6) or quotients (see Section 5.1.4). For products, we create the class ProductManifold that takes existing manifolds to construct a new one, and computations on each manifold can be done in parallel:

class ProductManifold(Manifold):
    r"""Class for a product of manifolds M_1 \times ... \times M_n."""
    def __init__(self, manifolds, n_jobs=1):
        self.dims = [manifold.dim for manifold in manifolds]
        super().__init__(dim=sum(self.dims))

# 3 Riemannian manifolds

## 3.1 Riemannian metrics

We now introduce a new structure on a differentiable manifold: the Riemannian metric, that allows to define the length of a curve, a distance function, a volume form, etc. Note that this additional structure may not be canonical, raising the thorny question of choosing the metric for the applications.

Recall that TM is the tangent bundle of a smooth manifold M (defined page 14).

**Definition 3.1.1 (Riemannian metric).** Let $M$ be a smooth $d$-dimensional manifold. A Riemannian metric on $M$ (or $TM$) is a family $((\cdot, \cdot)_p)_{p \in M}$ of inner products* on each tangent space $T_pM$, such that $((\cdot, \cdot)_p)$ depends smoothly on $p$. More formally, for any chart $\varphi, U$, and frame $(X_1, \dots, X_n)$ on $U$, the maps

$$p \mapsto \langle X_i(p), X_j(p) \rangle_p \quad 1 \le i, j \le n,$$

are smooth. A pair $(M, (\cdot, \cdot))$ is called a Riemannian manifold.

A metric is often written $g = (g_p)_{p \in M}$, where $g_p$ is the symmetric, positive definite (SPD) matrix representing the inner-product in a chart, that is

$$g_{ij}(p) = \langle (\partial_i)_p, (\partial_j)_p \rangle_p.$$

Alternatively, we sometimes define a metric with the notation $g = f(dx_1, \dots, dx_n)$ where $f$ is the quadratic form associated with $g$, and $dx_i$ represent vector coordinates (as linear forms). For example the usual Euclidean metric is $g = \sum_{i=1}^n dx_i^2$. The following theorem ensures that a metric is indeed a general structure. A proof can be found e.g. in Lafontaine et al., 2004, Theorem 2.2.

<footnote>*Defined in Appendix A.</footnote>
**Theorem 3.1.1 (Existence).** Any smooth manifold admits a Riemannian metric.

A Riemannian metric also defines a norm on $TM$, as usually defined by an inner-product on each tangent space:

$$ \forall x \in M, \forall v \in T_x M, \quad \|v\|_x = \sqrt{g_x(v, v)} = \sqrt{\langle v, v \rangle_x}. $$

**Implementation in geomstats** In order to guarantee flexibility, we have decided to keep manifolds and metrics separated in different objects. Indeed, although there may exist a canonical Riemannian metric on a given manifold, the choice of metric is not always natural, and researchers have struggled to find criteria to choose the right metric for the application at hand. The aim of `geomstats` is to allow researchers to compare different metrics on their problem, and in the future, to allow to *learn*, or optimize the metric (Louis et al., 2019; Hauberg, 2019).

We create an abstract `RiemannianMetric` class in `geomstats` to gather the basic attributes and methods expected of a metric. The most general way of defining a metric is to provide the `metric_matrix` method, that is $x \mapsto (g_{ij}(x))_{1 \le i,j \le d}$. By default, we use the identity matrix for all points, resulting in the Euclidean metric.

```python
class RiemannianMetric(Connection):
    """Class for Riemannian and pseudo-Riemannian metrics."""
    def __init__(self, dim, signature=None):
        super().__init__(dim=dim)
        if signature is None:
            self signature = (dim, 0)
    def metric_matrix(self, base_point=None):
        """Inner product matrix at base point."""
        return gs.eye(self.dim)
    def inner_product(self, tangent_vec_a, tangent_vec_b, base_point=None):
        """Inner product between two tangent vectors at a base point."""
        inner produ = self叹 metric matrix (base_point)
        inner produ = gs.einsum('...j,...jk,...k->...', tangent_vec_a, innerprodu, tangent_vec_b)
        return innerprodu
    def squared_norm(self, vector, base_point=None):
        """Compute the square of the norm of a vector."""
        return self标的 product (vector, vector, base_point)
    def norm(self, vector, base_point=None):
        """Compute norm of a vector."""
        sq norm = self.squared_norm (vector, base_point)
        return gs.sqrt (sq norm)
```

**Remark 3.1.1.**

1. The above class inherits from `Connection`. We indeed chose a class for an affine connection as parent class to a Riemannian metric because it is a more general structure, as explained in Section 3.2.1.
2. The attribute `signature` refers to the signature of the inner product, in case it is only a non-degenerate bilinear form, not necessarily positive. In this case the metric is called a pseudo-Riemannian metric.

We also add a `metric` property to the `Manifold` class, meaning that it is an attribute of the class, that can be set externally, calling the `setter` that checks that the given argument is indeed an instance of a `RiemannianMetric` object. Of course, all manifolds studied in this monograph come with a default metric, but users can choose to use different metrics or implement new metrics with only a few minimal operations. When closed-form solutions are available, the generic methods are overridden. With a `metric` attribute, a `Manifold` actually becomes a Riemannian manifold, but we did not think relevant to implement another layer of abstract class for Riemannian manifolds (i.e. a `RiemannianManifold` class), as all necessary operations are either in the `Manifold` object if they don't depend on the metric, or in the `RiemannianMetric` object if they do.

```python
@property
def metric(self):
    """Riemannian Metric associated to the Manifold."""
    return self._metric

@metric setter
def metric(self, metric):
    if metric is not None:
        if not isinstance(metric, RiemannianMetric):
            raise ValueError(
                'The argument must be a RiemannianMetric object')
        if metric.dim != self.dim:
            metric.dim = self.dim
    self._metric = metric
```

### Example 3.1: Euclidean metric

Let $M = \mathbb{R}^d$ be the standard vector space of dimension $d$, that is trivially a smooth manifold, and consider its standard inner-product defined for all $x, y \in \mathbb{R}^d$ by

$$ \langle x, y \rangle_2 = \sum_{i=1}^{d} x_i y_i = x^{\top} y. $$

As $T_x \mathbb{R}^d = \mathbb{R}^d$, it defines a metric on $\mathbb{R}^d$, which is referred to as the Euclidean metric.

### Example 3.2: Product metric

Let $(M, g)$ and $(M', g')$ be two Riemannian manifolds, and recall from Example 2.6 (page 9) that the Cartesian product $M \times M'$ is a manifold. It is also a Riemannian manifold. Indeed, define the *product metric* $g \oplus g'$ as the map defined at any $(x, x') \in$
$M \times M'$ and $\forall(v, v')$, $(w, w') \in T_x M \times T_{x'} M'$ by

$$g \oplus g'_{(x,x')}(v, v'), (w, w') = g(v, w) + g'(v', w').$$

In `geomstats`, it is possible to define such product metrics from existing objects of the class `RiemannianMetric`, and operations for each metric can be run in parallel if necessary.

```python
class ProductRiemannianMetric(RiemannianMetric):
    """Class for product of Riemannian metrics."""

    def __init__(self, metrics, default_point_type='vector', n_jobs=1):
        self.n_METRICS = len metrics)
        self.dims = [metric dim for metric in metrics]
        signatures = [metric signature for metric in metrics]

        sig_pos = sum(sign[0] for sign in signatures)
        sig谓 neg = sum(sign[1] for sign in signatures)
        super().__init__(
            dim=sum(dims), signature=(sig_pos, sig谓 neg),
            default_point_type=default_point_type
        )
```

Let $(N, g)$ be a Riemannian manifold, $M$ a smooth manifold, and $f: M \to N$ a map. Define the *pull-back* metric $(f^*g)$ on $M$ by

$$(f^*g)_x : \left\{ \begin{array}{rcl} T_x M \times T_x M & \to & \mathbb{R} \\ (v, w) & \mapsto & g_{f(x)}(df_x v, df_x w) \end{array} \right. . \qquad (3.1)$$

Now suppose that $f$ is an immersion. If $(f^*g)_x$ is non-degenerate and of constant signature for all $x \in M$, then $(M, f^*g)$ is a Riemannian manifold.

**Definition 3.1.2 (Isometry).** Let $(M, g)$ and $(M', g')$ be two Riemannian manifolds, and $f: M \to M'$. Then $f$ is called an *isometry* if it is a bijection and $f^*g' = g$.

When $M = M'$ and $g = g'$, we write $\text{Isom}(M)$ the set of isometries of $M$, and Myers and Steenrod (1939) showed that it is a Lie group that acts smoothly on $M$.

Now, consider that $M \subseteq N = \mathbb{R}^d$ is a submanifold of $M$, $f = i$ is the inclusion map. Then $g$ is called the embedding metric and $i^*g$ is its restriction to $M$. This case appears in many examples in `geomstats`.

### Example 3.3: Metric on the hypersphere

The hypersphere $S^d \subset \mathbb{R}^{d+1}$ endowed with the restriction of the Euclidean metric to $S^d$ is a Riemannian manifold. We call this metric the standard spherical metric, and implement it by calling the embedding metric.
class HypersphereMetric(RiemannianMetric):
    """Class for the Metric on the Hypersphere."""
    def __init__(self, dim):
        super().__init__(
            dim=dim, signature=(dim, 0))
        self embedding \_metric = EuclideanMetric(dim + 1)
    def metric\_matrix(self, base\_point=None):
        """Inner-product matrix at a base point."""
        return gs. eye(self. dim + 1)
    def inner\_product(self, tangent\_vec\_a, tangent\_vec\_b, base\_point=None):
        """Inner-product of two tangent vectors at a base point."""
        return self. embedding\_metric. inner\_product(
            tangent\_vec\_a, tangent\_vec\_b, base\_point)

Let us now define the Lorentz bilinear form of $\mathbb{R}^{d+1}$. It is the canonical bilinear form with signature (1, d), i.e. for any $x, y \in \mathbb{R}^{d+1}$

$$ \langle x, y \rangle_{\mathcal{L}} = -x_0 y_0 + \sum_{i=1}^{d} x_i y_i \qquad (3.2) $$

And write $\| \cdot \|_{\mathcal{L}}$ for the associated quadratic form. This is the underlying embedding metric in the following example.

### Example 3.4: Hyperbolic metric

Using the Lorentz form, recall that $H^d \in \mathbb{R}^{d+1}$ is the set of points such that $\|x\|_{\mathcal{L}} = -1$. Consider now the open subset $H_+^d$ of $H^d$:

$$ H_+^d = \{ x \in \mathbb{R}^{d+1} \mid x_0 > 0, \|x\|_{\mathcal{L}} = -1 \}. $$

Now, consider any $x \in H_+^d$ and two tangent vectors $v, w \in T_x M$. By Example 2.9, this means that $\langle v, x \rangle_{\mathcal{L}} = \langle w, x \rangle_{\mathcal{L}} = 0$. As $\langle \cdot, \cdot \rangle_{\mathcal{L}}$ is negative definite on $\mathbb{R}x = \{\lambda x \mid \lambda \in \mathbb{R}\}$, and its signature is (1, d), it is positive definite on the orthogonal of $\mathbb{R}x$ for the Lorentz metric, i.e. on $T_x M$. We can thus conclude that $(H_+^d, \langle \cdot, \cdot \rangle_{\mathcal{L}})$ is a Riemannian manifold.

In the implementation, we use the abstract class and simply override the `inner_product` function.
class HyperbolicMetric(RiemannianMetric):
    """Class for the hyperbolic metric."""
    def __init__(self, dim):
        super().__init__(
            dim=dim, signature=(1, dim))
    def metric_matrix(self, base_point=None):
        """Inner product matrix at base point."""
        diagonal = gs.array([-1.] + [1.] * self.dim)
        return from_vector_to_diagonal_matrix(diagonal)
    def inner_product(self, tangent_vec_a, tangent_vec_b, base_point=None):
        """Inner product between two tangent vectors at a base point."""
        diagonal = gs.array([-1.] + [1.] * self.dim)
        return gs.sum(diagonal * tangent_vec_a, tangent_vec_b, axis=-1)

### Example 3.5: Frobenius metric

The analog of the Euclidean metric on matrix spaces is the *Frobenius inner product* defined by

$$ \forall A, B \in M_{m,n}(\mathbb{R}), \quad \langle A, B \rangle_F = \operatorname{tr}(A^\top B) = \sum_{i,j} A_{ij}B_{ij}. $$

where tr is the trace operator. We use the right-hand-side expression in *geomstats* to avoid computing a matrix product ($O(mn^2)$) operations against $O(mn)$).

Endowed with this metric (or its restriction), $M_n(\mathbb{R})$, $GL(n)$ (Example 4.1 page 45), $SO(n)$ (Example 2.3 page 6) and $SE(n)$ (Example 4.7 page 51) are Riemannian manifolds.

## 3.2 Affine connections and the Levi-Civita connection

A connection is an additional structure that can be defined independently of a Riemannian metric. It provides a way to compare tangent spaces from one point to another, by defining the notion of parallelism. For a detailed and historical account of the different approaches to defining connections, we refer to Marle (2005).

### 3.2.1 Connections

**Definition 3.2.1 (Connection).** Let $M$ be a smooth manifold. A *connection* on $M$ is an $\mathbb{R}$-bilinear map $\nabla : \Gamma(TM) \times \Gamma(TM) \to \Gamma(TM)$ that verifies for all $X, Y \in \Gamma(TM), \forall f \in C^\infty(M)$:

1. (Linearity of 1st argument) $\nabla_{fX}Y = f\nabla_XY,$
2. (Leibniz rule in 2nd argument) $\nabla_X(fY) = X(f)Y + f\nabla_X Y$.

The vector field $\nabla_X Y$ is called the *covariant derivative* of $Y$ w.r.t. $X$.

In fact, $(\nabla_X Y)_p$ only depends on the value of $X$ at $p$ and not in its neighborhood. In contrast, it does depend on $Y$ around $p$.

In local coordinates $(x^1, \dots, x^d)$, the *Christoffel symbols* are used to specify the connection. Recall that $(\partial_i)_i$ is a local frame, so we can decompose $\nabla$ in this basis, and define $(\Gamma^k_{ij})_{ijk}$ such that

$$ \nabla_{\partial_i}(\partial_j) = \Gamma^k_{ij} \partial_k, $$

where we used Einstein summation convention, meaning that a sum occurs along the indices that appear both in subscript and superscript, here $k$.

Two vector fields $X, Y \in \Gamma(TM)$ can be decomposed locally in coordinates, writing $X = X^i \partial_i$ and $Y = Y^i \partial_i$ where $X^i, Y^i$ are coordinate functions defined locally. Then using the properties of a connection,

$$ \begin{aligned} \nabla_X Y &= X^i \nabla_{\partial_i} (Y^j \partial_j) \\ \nabla_X Y &= X^i \left( \frac{\partial Y^j}{\partial x^i} \partial_j + Y^j \Gamma^k_{ij} \partial_j \right). \end{aligned} \quad (3.3) $$

As this formula shows, a connection provides a correction term when compared to the directional derivative of $Y$ with respect to $X$ in a chart.

### 3.2.2 Parallel transport and geodesics

We now focus on vector fields that are defined along a curve $\gamma : [a, b] \to M$, i.e. a smooth map $X : [a, b] \to TM$ such that at any $t \in [a, b]$, $X(t) \in T_{\gamma(t)}M$. Note that such vector field need not be defined on the whole manifold, but can be locally extended to an open set around every point. Thankfully, one can show that there exists a *covariant derivative* that coincides with $(\nabla_{\dot{\gamma}(t)}X)_{\gamma(t)}$ at all $t \in [a, b]$ and for any $X$ defined on a neighborhood of $\gamma(t)$ (Lafontaine et al., 2004, Theorem 2.68). We will skip these technicalities and admit that $\nabla_{\dot{\gamma}}X$ is well-defined for any vector field $X$ along $\gamma$.

We now define a central notion in geometric statistics: the parallel transport.

**Definition 3.2.2 (Parallel vector field).** Let $M$ be a smooth manifold and $\nabla$ a connection on $M$. For any curve $\gamma : [a, b] \to M$ in $M$, a vector field $X$ along $\gamma$ is *parallel* if

$$ \nabla_{\dot{\gamma}(t)} X(t) = 0. \quad (3.4) $$

In a local chart and using the Christoffel symbols and in particular using Equation (3.3), Equation (3.4) can be written $\forall t \in [a, b]$

$$ \dot{X}^k(t) + \Gamma^k_{ij} X^i(t) \dot{\gamma}^j(t) = 0. \quad (3.5) $$
From the properties of ODEs, one can prove the following existence and uniqueness property (Lafontaine et al., 2004, Proposition 2.72)

**Proposition 3.2.1.** Let $M$ be a smooth manifold and let $\nabla$ be a connection on $M$. For every $C^1$ curve $\gamma: [a, b] \to M$ in $M$, for every $t \in [a, b]$ and every $v \in T_{\gamma(t)}M$, there is a unique parallel vector field $X$ along $\gamma$ such that $X(t) = v$.

For such parallel vector field $X$ and $s \in [a, b]$, we thus call $X(s)$ the *parallel transport* of $v$ along $\gamma$ from $t$ to $s$, and write $X(s) = \Pi_{\gamma,t}^s v$. Another consequence of the properties of ODEs is that for any $s, t \in [a, b]$, $\Pi_{\gamma,t}^s$ is a linear isomorphism between the tangent spaces $T_{\gamma(t)}M$ and $T_{\gamma(s)}M$.

<figure><img src="image_11.png" /><figcaption>**Figure 3.1:** Representation of a parallel vector field $X$ (green) along a curve $\gamma$ (blue). The orientation of $X$ with respect to the speed $\dot{\gamma}$ of $\gamma$ stays the same.</figcaption></figure>

Intuitively, the parallel transport equation constrains the transported vector to keep a certain angle w.r.t the speed $\dot{\gamma}$ of the curve while moving along it, see Figure 3.1 for an illustration. Conversely, the connection $\nabla_X Y$ could be retrieved from the parallel transport $\Pi_{\gamma,t}^s Y$ when the curve is shrunk to a point in the direction of $X$, i.e. $s \to t$ and $\dot{\gamma}(t) = X$.

We now focus on a particular set of curves on $M$ for which the velocity is parallel transported along the curve.

**Definition 3.2.3.** Let $M$ be a smooth manifold endowed with a connection $\nabla$. A curve $\gamma: [a, b] \to M$ is said to be *autoparallel* and is called a *geodesic* of $(M, \nabla)$ if it satisfies for all $t \in [a, b]$

$$ \nabla_{\dot{\gamma}(t)} \dot{\gamma}(t) = 0. \qquad (3.6) $$

Note that $\nabla_{\dot{\gamma}} \dot{\gamma}$ can be interpreted as the covariant acceleration of $\gamma$, so that Equation (3.6) constrains geodesics to be zero acceleration curves. The flow of the geodesic equation is called geodesic flow, and is a fundamental example of dynamical system, that generalizes straight lines from Euclidean spaces.

The geodesic equation can be written in local coordinates like the parallel transport equation. A geodesic curve $\gamma$ satisfies for all times $t \in [a, b]$

$$ \ddot{\gamma}^k(t) + \Gamma^k_{ij} \dot{\gamma}^i(t) \dot{\gamma}^j(t) = 0. \qquad (3.7) $$
From the properties of second-order differential equations, for any $(x, v) \in TM$, there exists a maximal interval $I_{x,v} \subseteq \mathbb{R}$ such that $\gamma_{x,v}: I_{x,v} \to M$ is the unique geodesic that verifies $\gamma_{x,v}(0) = x$ and $\dot{\gamma}_{x,v}(0) = v$. Moreover, by homogeneity of the Equation (3.7), for any $s > 0$, $I_{x,sv} = \frac{1}{s}I_{x,v}$ and $\gamma_{x,sv}(t) = \gamma_{x,v}(st)$. We deduce that the set of vectors in $T_x M$ such that $1 \in I_{x,v}$ is non empty, open, and contains 0. This leads to the following definition

**Definition 3.2.4 (Exponential map).** Let $\nabla$ be a connection on a smooth manifold $M$. The map $(x, v) \mapsto \gamma_{x,v}(1)$ defined on the open set $\{(x, v) \in TM, 1 \in I_{x,v}\}$ and with values in $M$ is called the exponential map of $\nabla$. We say that $M$ is geodesically complete if the exponential map is defined on the whole of $TM$.

For any $x \in M$, we write the exponential map at $x$ $\text{Exp}_x : v \in T_x M \mapsto \gamma_{x,v}(1)$.

**Remark 3.2.1.** Note that although we introduced many notions in this paragraph, there were no examples. Indeed a connection or its Christoffel symbols are rarely explicit, except when the connection is compatible with a metric. We explain this notion of compatibility in the next paragraph, and will give several examples of geodesics and parallel transports.

### 3.2.3 The Levi-Civita connection

First, given two vector fields $Y, Z$, recall that $\langle Y, Z \rangle$ is a smooth function on $M$ so $X(\langle Y, Z \rangle)$ must be understood as $X(f) = df(X)$ for $f: p \mapsto \langle Y_p, Z_p \rangle_p$.

The following is considered the fundamental theorem of Riemannian geometry. It ensures that there exists a unique connection that is “compatible” with the metric. See Lafontaine et al. (2004, Theorem 2.51) for a proof.

**Theorem 3.2.1.** Let $(M, g)$ be a Riemannian manifold. There is a unique connection on $M$ that verifies for all $X, Y, Z \in \Gamma(TM)$

$$1. \text{ (Torsion-free)} \qquad \nabla_X Y - \nabla_Y X = [X, Y] \qquad (3.8)$$

$$2. \text{(Compatibility)} \qquad X(\langle X, Y \rangle) = \langle \nabla_X Y, Z \rangle + \langle Y, \nabla_X Z \rangle \qquad (3.9)$$

This connection is called the *Levi-Civita* connection and is determined by the Koszul formula:

$$2\langle\nabla_X Y, Z\rangle = X(\langle Y, Z\rangle) + Y(\langle X, Z\rangle) - Z(\langle X, Y\rangle) \\ - \langle Y, [X, Z]\rangle - \langle X, [Y, Z]\rangle - \langle Z, [Y, X]\rangle. \qquad (3.10)$$

The notion of compatibility is thus detailed in Equations (3.8) and (3.9). The former is quite general but ensures uniqueness of the Levi-Civita connection, it is called the zero-torsion condition. The latter can be understood as a Leibniz rule where $\langle Y, Z \rangle$ is seen as a product and the derivative is $\nabla_X$. More precisely, this condition means that the metric is
parallel with respect to the connection. Indeed, although we introduced a connection as a map on vector fields, it can be extended to tensors of any order. For a 0-order tensor, i.e. a smooth function $f$, we have $\nabla_X f = X(f)$, and for a (2, 0)-tensor such as $g$, we have (by generalising the Leibniz rule) for any $X, Y, Z \in \Gamma(TM)$

$$ (\nabla_X g)(Y, Z) = \nabla_X(g(Y, Z)) - g(\nabla_X Y, Z) - g(Y, \nabla_X Z). \\ = X(\langle Y, Z \rangle) - \langle \nabla_X Y, Z \rangle - \langle Y, \nabla_X Z \rangle \quad (3.11) $$

Therefore the combination of Equations (3.9) and (3.11) results in $\nabla g = 0$.

In the general case, the Levi-Civita connection is characterized locally by its Christoffel symbols. Let $(x^1, \dots, x^n)$ be a local coordinate chart, by definition

$$ \nabla_{\partial_i} \partial_j = \Gamma^k_{ij} \partial_k $$

The torsion-free condition (3.8) and Schwartz theorem ($[\partial_i, \partial_j] = 0$) thus imply the symmetry of the symbols: $\Gamma^k_{ij} = \Gamma^k_{ji}$ for all $1 \le i, j, k \le n$. Furthermore, together with (3.9), Koszul formula (3.10) translates into

$$ 2g(\nabla_{\partial_i} \partial_j, \partial_k) = \partial_i g_{jk} + \partial_j g_{ki} - \partial_k g_{ij}. $$

Writing $(g^{ij})_{ij} = (g_{ij})^{-1}$ for the inverse of the metric matrix, we obtain

$$ \Gamma^k_{ij} = \frac{1}{2} g^{lk} (\partial_i g_{jl} + \partial_j g_{li} - \partial_l g_{ij}) \quad (3.12) $$

Thus the Christoffel symbols can be computed from the metric $g$. This formula is rarely used by mathematicians for computations “by hand”, but we shall use it in geomstats to implement pull back metrics.

The following characterization allows to compute the Levi-Civita connection in the case of embedded manifolds equipped with the embedding metric, as is the case in many of our examples. See Lafontaine et al. (2004, Proposition 2.56) for a proof.

**Proposition 3.2.2.** Let $M$ be a Riemannian manifold and $N$ a submanifold of $M$ endowed with the induced metric (i.e., the restriction of the embedding metric). If $\nabla^M$ and $\nabla^N$ are the Levi-Civita connections on $M$ and $N$ respectively defined by the metric on $M$. Then for any vector field $X, Y \in \Gamma(TN)$, we have

$$ \nabla_X^N Y = (\nabla_X^M Y)^\|, \quad (3.13) $$

where $(v)^\|$ is the orthogonal projection of any $v \in T_pM$ onto $T_pN$ for any $p \in N$.

This proposition is particularly useful for embedded manifolds of $\mathbb{R}^N$ endowed with the ambient Euclidean metric, as it is straightforward to see (by uniqueness) that their Levi-Civita connection coincide with the directional derivative of functions from $\mathbb{R}^N$ to $\mathbb{R}^N$.
Therefore, one only needs to compute the projection to the tangent spaces of any point of an embedded manifold to compute its connection.

Finally, thanks to the compatibility of the metric with the connection, the parallel transport preserves the metric, as stated in the following proposition (Lafontaine et al., 2004, Proposition 2.74).

**Proposition 3.2.3.** Let $\gamma: I \to M$ be a smooth curve, and $s, t \in I$. Then the parallel transport map $\Pi_{\gamma,s}^t: T_{\gamma(s)}M \to T_{\gamma(t)}M$ along $\gamma$ for the Levi-Civita connection is an isometry, i.e.

$$ \forall v, w \in T_{\gamma(s)}M, \quad \langle \Pi_{\gamma,s}^t v, \Pi_{\gamma,s}^t w \rangle_{\gamma(t)} = \langle v, w \rangle_{\gamma(s)}. $$

This justifies the use of parallel transport in statistical procedures to move data from one reference point to another while preserving distances.

Moreover, the compatibility with the metric also ensures that isometries preserve the Levi-Civita connection and its geodesics. The following can be deduced from the Kozsul formula (3.10).

**Proposition 3.2.4.** Let $(M, g)$ be a Riemannian connection and $\nabla$ its Levi-Civita connection. Let $f \in \text{Isom}(M)$ be an isometry of $(M, g)$, then for any $X, Y \in \Gamma(TM)$ and $x \in M$, we have

$$ df \nabla_X Y = \nabla_{df_X} df Y \qquad (3.14) $$

$$ \mathrm{Exp}_{f(x)} \circ df_x = f \circ \mathrm{Exp}_x. \qquad (3.15) $$

We now come to the definition of distance and its link with geodesics.

## 3.3 Distance and Geodesics

### 3.3.1 Injectivity and parametrization

Let $(M, g)$ be a Riemannian manifold. As the Levi-Civita is uniquely defined, we call geodesic of $(M, g)$ any geodesic of its Levi-Civita connection. Similarly, we define:

**Definition 3.3.1 (Exp and Log maps).** We call *Riemannian exponential* the exponential map of the Levi-Civita connection, and define its *injectivity radius* $\text{inj}_M(x)$ at $x$ as the greatest $\epsilon > 0$ such that $\text{Exp}$ is a diffeomorphism on the open ball of radius $\epsilon$ of $T_x M$.

The following properties of the exponential map can be proved (Lafontaine et al., 2004, Proposition 2.88).

**Proposition 3.3.1.** Let $(M, g)$ be a Riemannian manifold and $x \in M$.

* The differential of $\mathrm{Exp}_x$ at 0 is the identity.
• $(x, v) \mapsto (x, \mathrm{Exp}_x(v))$ is a smooth diffeomorphism from an open neighborhood of the null section (i.e. "$M \times \{0\}''$) in $TM$ to an open neighborhood of the diagonal of $M \times M$.

The open ball of radius $\mathrm{inj}_M(x)$ on which the Exp map is a diffeomorphism can be maximally extended into a larger star-shaped open set called the *injectivity domain* $\mathrm{Inj}(x) \subseteq T_x M$. The *Riemannian logarithm* is defined as the inverse of the Exp map on this *injectivity domain*. The injectivity radius of $M$, $\mathrm{inj}_M$, is then the smallest of all injectivity radii at $x$ for $x \in M$. It may not be finite, as in e.g. the hyperbolic space.

By Equation (3.9), the squared velocity $E(t) = \frac{1}{2}\|\dot{\gamma}(t)\|^2$ of any geodesic $\gamma$ is constant in time:

$$ \frac{d}{dt} \langle \dot{\gamma}(t), \dot{\gamma}(t) \rangle = 2 \langle \nabla_{\dot{\gamma}(t)} \dot{\gamma}(t), \dot{\gamma}(t) \rangle = 0 $$

We say that a curve $\gamma$ is parametrized with constant speed if the function $t \mapsto \|\dot{\gamma}(t)\|_{\gamma(t)}$ is constant. The above equation thus shows that geodesics are parametrized with constant speed.

Furthermore, if $N \subset M$ is an embedded manifold, (3.13) implies that the geodesics of $N$ are the curves in $M$ whose acceleration is normal to $N$.

### Example 3.6: Geodesics of the hypersphere

For any $x \in S^d \subset \mathbb{R}^{d+1}$, and $v \in T_x S^d$, define the curve on $S^d$ that parametrizes the great circle

$$ \gamma : t \mapsto \cos(t\|v\|)x + \sin(t\|v\|) \frac{v}{\|v\|} \qquad (3.16) $$

Then, it is clear that the acceleration $\ddot{\gamma}(t) = -\|v\|^2\gamma(t)$ is normal to $S^d$. By uniqueness this curve is the geodesic from $x$ with initial velocity $v$.

The Exp map at $x$, $v \mapsto \cos(t\|v\|)x + \sin(t\|v\|) \frac{v}{\|v\|}$ is thus well defined and is a diffeomorphism from the open ball of radius $\pi$ to $S^d \setminus \{-x\}$, i.e. onto the entire sphere except the antipodal point of $x$. Its inverse is defined for any $y \notin \{x, -x\} \in S^d$ by

$$ \mathrm{Log}_x(y) = \arccos(\langle y, x \rangle) \frac{y - \langle y, x \rangle x}{\|y - \langle y, x \rangle x\|} \qquad (3.17) $$

We implement the Riemannian Exp and Log maps in the *HypersphereMetric*.

**Implementation trick 3.1.** To improve numerical stability around 0, we use a Taylor approximation of the sinc and cosine functions for inputs smaller than $10^{-6}$. As both functions are even, the squared norm of $v$ is computed and taking the square-root is not necessary.
For a $C^1$-curve $\gamma : [a, b] \to M$, we now define its length $L$ and energy $E$ (also called action integral in physics) by:

$$L(\gamma) = \int_{a}^{b} \|\dot{\gamma}(t)\| dt, \quad E(\gamma) = \int_{a}^{b} \|\dot{\gamma}(t)\|^2 dt$$

Note that the length does not depend on the parametrization of the curve, while the energy does. Indeed, moving along a path from $a$ to $b$ does not require the same energy if the speed is increased, but the length remains the same. This finally leads to the definition of distance.

### 3.3.2 Distance and completeness

**Definition 3.3.2 (Riemannian distance).** Let $(M, g)$ be a Riemannian manifold, and $x, y \in M$. The Riemannian distance between $x$ and $y$ is the lower bound of the lengths of all piecewise smooth curves joining $x$ to $y$:

$$d(x, y) = \inf\{L(\gamma) \mid \gamma : I \to M \text{ piecewise } C^1, \gamma(0) = x, \gamma(1) = y\}.$$

We say that $\gamma$ is minimizing if $L(\gamma) = d(x, y)$.

If $M$ is connected*, this distance function is indeed a distance, and the induced topology coincides with that of $M$ (Paulin, 2014, Proposition 3.13). We say that a curve $\gamma : [a, b] \to M$ is parametrized by (or resp. proportional to) arc-length if $L(\gamma) = b - a$ (resp. $\exists \lambda > 0, L(\gamma) = \lambda(b - a)$). We now see that the geodesics of $(M, g)$ are the minimizing curves, and further minimize the total energy (Paulin, 2014, Proposition 3.14).

**Theorem 3.3.1.** Let $(M, g)$ be a Riemannian manifold, and $\gamma : [a, b] \to M$ a $C^1$ curve. The following assertions are equivalent

* $\gamma$ is a geodesic;

* $\gamma$ is parametrized with constant velocity and locally minimizing;

* $\gamma$ is locally energy minimizing.

In particular, as for any $(x, v) \in TM$, the curve $\gamma : t \in \text{Inj}(x) \mapsto \text{Exp}_x(tv)$ is a geodesic, it is locally length-minimizing, so that for any $y \in \text{Exp}_x(\text{Inj}(x))$, $d(x, y) = L(\gamma) = \|v\|$. By definition, $v = \text{Log}_x(y)$, so that

$$d(x, y) = \| \text{Log}_x(y) \|.$$

We add this as a method to the RiemannianMetric class of geomstats.

<footnote>*Defined in Appendix A.</footnote>
```python
def dist(self, point_a, point_b):
    """Riemannian distance between two points"""
    log = self.log(point_b, base_point=point_a)
    return self.norm(log, base_point=point_a)
```

Furthermore, the last point of Theorem 3.3.1 is a variational principle, i.e. the minimization of a function on the space of paths.

Recall from Definition 3.2.4 that a Riemannian manifold is called geodesically complete if the Exp map is defined on the whole of $TM$, meaning that geodesics are defined on all $\mathbb{R}$. On the other hand, the more usual notion of completeness of a metric space is defined as follows:

**Definition 3.3.3 (Complete metric space).** A metric space is called complete if every Cauchy sequence converges in that space.

**Example 3.7: Completeness of the half-plane**

Let the half-plane $P = \{(x, y) \in \mathbb{R}^2 \mid y > 0\}$ be equipped with the canonical Euclidean metric of $\mathbb{R}^2$. Then $P$ is obviously not geodesically complete, as geodesics are straight lines exiting $P$.

Now define the metric $g(x, y) = \frac{dx^2+dy^2}{y^2}$. ($P, g$) is now a geodesically complete Riemannian manifold.

The following theorem gives a sufficient condition for Riemannian manifolds to share the same nice properties as complete metric spaces.

**Theorem 3.3.2 (Hopf-Rinow).** Let $(M, g)$ be a Riemannian manifold. If $M$ is geodesically complete, then any two points of $M$ can be joined by a minimizing geodesic.

**Corollary 3.3.3.** Let $(M, g)$ be a connected Riemannian manifold, and $d$ the Riemannian distance. The following are equivalent

* $(M, g)$ is geodesically complete;

* every closed and bounded subset of $(M, d)$ is compact;

* the metric space $(M, d)$ is complete (as a metric space).

**Remark 3.3.1 (Infinite dimension).** The Hopf-Rinow theorem is not true in infinite dimension. Counter-examples have been constructed and exhibit manifolds that are complete as metric spaces but for which points exist that cannot be joined by a minimizing geodesic, or even by any geodesic at all. For more details, see Atkin (1975).

In the previous example of the hypersphere, geodesics are great circles and are known in closed-form. To finish this subsection, we detail our implementation of the Riemannian distance when no closed-form solution is available for the geodesics.
### 3.3.3 Numerical approximations in `geomstats`

In the most general case, recall that geodesics are defined by an ODE expressed either in a chart or in the ambient space (e.g. Equation (3.7)). Discrete integration methods can thus be used to approximate the exponential map, and optimization algorithms to compute the logarithm. Inspired by Kühnel et al. (2019), we implemented those methods in `geomstats`.

Consider the geodesic equation as a coupled system of first-order ODEs:

$$ egin{cases} v(t) = \dot{\gamma}(t) \\ \dot{v}(t) = f(v(t), \gamma(t), t) \end{cases} \qquad (3.18) $$

where $f$ is a smooth function. A method of the class `Connection` is used to compute the right-hand-side of (3.18), called `geodesic_equation` where the state variable is $(\gamma(t), \dot{\gamma}(t))$:

```python
def geodesic_equation(self, state, time):
    """Return the right-hide-side of the geodesic equation."""
    position, velocity = state
    gamma = self.christoffels(position)
    equation = -gs.einsum(...kij,...i,...j->...k', gamma, velocity, velocity)
    return gsstack([velocity, equation])
```

Given initial conditions, a first-order forward Euler scheme, or higher-order Runge-Kutta method can be used to integrate this system.

```python
STEP_FUNCTION5 = {'euler': 'euler_step', 'rk2': 'rk2_step'}

def euler_step(force, state, time, dt):
    """Compute one step of the Euler approximation."""
    derivatives = force(state, time)
    new_state = state + derivatives * dt
    return new_state

def rk2_step(force, state, time, dt):
    """Compute one step of the rk2 approximation."""
    k1 = force(state, time)
    k2 = force(state + dt / 2 * k1, time + dt / 2)
    new_state = state + dt * k2
    return new_state

def integrate(
    function, initial_state, end_time=1.0, n_steps=10, step='euler'):
    """Compute the flow of a vector field."""
    dt = end_time / n_steps
    states = [initial_state]
    step_function = globals()[STEP_FUNCTION5[step]]
    current_state = initial_state
    for i in range(n_steps):
        current_state = step_function(
            state = current_state, force = function, time = i * dt, dt = dt)
        states.append(current_state)
    return states
```

We can thus add the following method to the class `Connection` of `geomstats` to compute the exponential map
```python
def exp(self, tangent_vec, base_point, n_steps=10, step='euler', **kwargs):
    """Exponential map associated to the affine connection."""
    initial_state = gs.stack([base_point, tangent_vec])
    flow = integrate(
        self.geodesic_equation, initial_state, n_steps=n_steps, step=step)
    exp, velocity = flow[-1]
    return exp
```

On the other hand, the Log map corresponds to a boundary value problem (BVP), and an optimization procedure is required. We choose the geodesic shooting method, that solves the following problem (P) which corresponds to an energy minimization. Working in a convex neighborhood, it admits a unique minimizer $v^*$:

$$ \min d^2(\mathrm{Exp}_x(v), y) \qquad (\mathcal{P}) $$

As in a coordinate chart or in the embedding space, all norms are equivalent, we use the appropriate Euclidean metric (written $||.||_2$) for practical purpose.

$$ \| \mathrm{Exp}_x(v) - y \|_2^2. \qquad (\mathcal{P}') $$

and the problem is solved by gradient descent (GD) until a convergence tolerance $\epsilon$ is reached. We use the `minimize` function from the `scipy` function, and use automatic differentiation to compute the gradient of the exponential map.

```python
def log(self, point, base_point, n_steps=N_STEPS, step='euler',
        max_iter=25, verbose=False, tol=gs.atol):
    """Compute logarithm map associated to the affine connection."""
    maxShape = point.shape if point.ndim > base_point.ndim else \n
               base_point.shape

def objective(velocity):
    """Define the objective function."""
    velocity = gs.array(velocity, ndtype=base_point.dtype)
    velocity = gs Reshape(velocity, maxShape)
    delta = self .exp(velocity, base_point, n_steps, step) - point
    return gs .sum delta ** 2)

objective_withGrad = gs .autograd.value_andGrad(objective)
tangent_vec = gs .flatten (gs .random .rand (*maxShape))
res = minimize(
    objective_withGrad, tangent_vec, method='L-BFGS-B', jac=True,
    options={ 'disp': verbose, 'maxiter': max_iter}, tol=tol)

tangent_vec = gs .array(res.x, ndtype=base_point.dtype)
tangent_vec = gs .reshape(tangent_vec, maxShape)
return tangent_vec
```

## 3.4 Curvature

### 3.4.1 Definition and Properties

In a Euclidean space, a constant field is parallel along any curve. In a Riemannian manifold in general, parallel transport depends on the curve followed, and there may not exist fields that are parallel along all curves, not even locally. One can investigate the effect
of parallel transport along small closed curves. Consider a point $x \in M$ and the closed curves whose tangent velocities at that point span a subspace of dimension two. They introduce a deviation of parallel transport from the identity map of the tangent space at that point, and this deviation can be shown to depend only on a basis of this plane. This is due to the commuting properties of the covariant derivative, i.e. the difference between evaluating $\nabla_X$ after $\nabla_Y$ and vice-versa. This difference is not sufficient however to define a tensor. The following lemma gives a sufficient condition for a map to define a tensor on a manifold (Lafontaine et al., 2004, Proposition 1.114).

**Lemma 3.4.1 (Tensoriality).** Let $p \in \mathbb{N}$ and $A : \Gamma(TM)^p \to \Gamma(TM)$ be a $C^\infty(M)$-multilinear map, i.e. $\forall f_1, \dots, f_p \in C^\infty(M), \forall X_1, \dots, X_p \in \Gamma(TM), A(f_1 X_1, \dots, f_p X_p) = f_1 \dots f_p A(X_1, \dots, X_p)$. Then for any $x \in M$, the value of $A(X_1, \dots, X_p)$ at $x$ only depends on the value of the $X_i$'s at $x$.

Simple computations using the Leibniz rule from the definition of a connection (2) show that the following defines a $C^\infty$-multilinear map (see Lee (2018, Proposition 7.3) for the computations).

**Definition 3.4.1 (Curvature tensor).** Let $(M, \nabla)$ be a manifold equipped with a connection. The curvature tensor of $\nabla$ is defined as the map from $\Gamma(TM)^3$ to $\Gamma(TM)$ by

$$R(X, Y)Z = \nabla_X \nabla_Y Z - \nabla_Y \nabla_X Z - \nabla_{[X, Y]} Z \quad (3.19)$$

The curvature tensor of a Riemannian manifold $(M, g)$ is the curvature of its Levi-Civita connection.

In the above definition, the arguments $X, Y, Z$ are vector fields, but the tensoriality Lemma 3.4.1 allows to write at any $x \in M$, $R_x$ as a map defined on tangent vectors $u, v, w \in T_x M$.

### Example 3.8: Curvature of the Euclidean space

In a Euclidean space $\mathbb{R}^d$ with the canonical inner-product, the connection coincides with the directional derivative, which commutes. Therefore $R = 0$.

The following properties help cope with the complexity of this tensor. Proofs can be found in Lafontaine et al. (2004, Proposition 3.5), Gallier and Quaintance (2020, Proposition 16.1, 16.3).

**Proposition 3.4.1.** Let $(M, g)$ be a Riemannian manifold and $R$ its curvature tensor. The following properties hold
1. (Skew-symmetry) The map $(X, Y) \in \Gamma(TM)^2 \to R(X, Y)$ is skew-symmetric.

2. (Bianchi's identity) For any $X, Y, Z \in \Gamma(TM)$, we have

$$R(X, Y)Z + R(Y, Z)X + R(Z, X)Y = 0.$$

3. (Bianchi's second identity) For any $X, Y, Z \in \Gamma(TM)$, we have

$$\nabla_X R(Y, Z) + \nabla_Y R(Z, X) + \nabla_Z R(X, Y) = 0.$$

4. For any isometry $f \in \text{Isom}(M)$ and $X, Y, Z \in \Gamma(TM)$,

$$R(dfX, dfY)(dfZ) = df(R(X, Y)Z)$$

5. For all $X, Y \in \Gamma(TM)$, the field of linear maps $R(X, Y)$ is skew-symmetric with respect to the metric $g$, i.e. $\forall W, Z \in \Gamma(TM)$

$$g(R(X, Y)Z, W) = -g(Z, R(X, Y)W)$$

6. For all $X, Y, Z, W \in \Gamma(TM)$, we have

$$g(R(X, Y)Z, W) = g(R(Z, W)X, Y) = g(R(W, Z)Y, X)$$

Note that the compatibility of the metric (Equations (3.9), (3.8)) are required for the last two assertions.

### 3.4.2 Sectional curvature

The previous proposition allows to define the sectional curvature for any $x \in M$ and $u, v \in T_x M$ such that $u, v$ are not collinear

$$\kappa_x(u, v) = \frac{\langle R_x(u, v)v, u \rangle}{\|u\|^2 \|v\|^2 - \langle u, v \rangle^2} \quad (3.20)$$

From the above properties, the value of $\kappa_x$ in fact only depends on the plane spanned by $(u, v)$, i.e. for any non-vanishing linear combination $\alpha u + \beta v$, $\kappa(\alpha u + \beta v, v) = \kappa(u, v)$. In fact, this scalar quantity is enough to characterize the curvature of $M$. Indeed, one can prove from algebraic computations only (see e.g. Kobayashi and Nomizu (1996a, Chapter V, Proposition 1.2)) the following

**Lemma 3.4.2.** Let two quadrilinear mappings $A, B$ defined on a vector space $V$ that both verify the properties $\forall u, v, w, z \in V$

$$A(u, v, w, z) = -A(v, u, w, z) = -A(v, u, z, w)$$

$$A(u, v, w, z) + A(u, w, z, v) + A(u, z, v, w) = 0$$

and coincide for any two variables $u, v \in V$ in the sense: $A(u, v, u, v) = B(u, v, u, v)$. Then $A = B$.
This applies to the map $(u, v, w, z) \mapsto g(R(w, z)v, u)$ and one can thus show the following theorem (Lafontaine et al., 2004, Theorem 3.8).

**Theorem 3.4.3.** The sectional curvature determines the curvature tensor.

In particular, if for every $x \in M$, $\kappa$ has a constant value $\kappa_x$ on every planes of $T_x M$, then we can compute the curvature up to this constant. Indeed from the definition (3.20) we have for any $u, v \in T_x M$

$$ \langle R(u, v)v, u \rangle = \kappa_x (\|u\|^2 \|v\|^2 - \langle u, v \rangle^2) $$

Consider the maps $A : (u, v, w, z) \mapsto \langle R(u, v)w, z \rangle$ and $B$ defined by

$$ B(u, v, w, z) = (\langle u, w \rangle \langle v, z \rangle - \langle u, z \rangle \langle v, w \rangle) $$

It is clear that the hypotheses of Lemma 3.4.2 are verified for $A$ and $\kappa_x B$, so we conclude

$$ \langle R(u, v)w, z \rangle = (\langle u, w \rangle \langle v, z \rangle - \langle u, z \rangle \langle v, w \rangle) \kappa_x. \quad (3.21) $$

For example, Equation (3.21) is always valid in dimension $d = 2$, as there is only one plane in each tangent space.

We now state the surprising result of F. Schur that gives a sufficient condition for formula_3.21 to hold. A proof can be found in Kobayashi and Nomizu (1996a, Theorem 2.2).

**Theorem 3.4.4 (Schur).** Let $(M, g)$ be a connected Riemannian manifold of dimension $d$. If $d \ge 3$ and if the sectional curvature $\kappa$ does not depend on the plane but only on the point $x \in M$, then $\kappa$ is constant (i.e. does not even depend on $x$).

**Example 3.9: Sectional curvature of $S^d$ and $H_+^d$**

Recall that the rotation group $SO(d+1)$ acts on the hypersphere $S^d$, and suppose that $d \ge 2$. The stabilizer of this action on any $x \in S^d$ is isomorphic to $SO(d)$, whose action on $T_x S^d$ is transitive on 2-planes, i.e. any 2-plane can be mapped by a rotation to any other 2-plane.

By assertion 4. of Proposition 3.4.1, the sectional curvature is preserved from one plane to another by rotations, it is therefore constant on the whole $T_x M$, write this value $\kappa_x$. If $d \ge 3$, we can use Theorem 3.4.4 to conclude it is constant on $S^d$. Of course, it is also constant for $d = 2$, as we can use an isometry to map the curvature tensor from one point to any other.

The same argument, using this time the group $O(1, d)$ of isometries of the Lorentz product, that acts transitively on the hyperboloid $H_+^d$, shows that this space also has constant curvature.
<figure><img src="image_12.png" /><figcaption>**Figure 3.2:** Drawing of the distortion of geodesics compared to Euclidean tangent vectors. Arrows represent tangent vectors, thin black lines represent geodesics and we use the notation $x_v = \text{Exp}_x(v)$. Curvature modifies the geodesic distance $d$ compared with the distance between $rv$ and $rw$ in $T_x M$ (adapted from Paulin, 2014, Section 3.6.3).</figcaption></figure>

**Definition 3.4.2 (Constant curvature).** Let $(M, g)$ be a Riemannian manifold. $(M, g)$ is said to have *constant* (resp. *negative*, resp. *positive*) curvature if it has constant (resp. negative, resp. positive) sectional curvature.

In fact, the (complete simply connected) constant curvature spaces are all isometric to one of the Euclidean space (flat), the hypersphere (positive curvature) or the hyperbolic space (negative curvature) (see e.g. Lafontaine et al., 2004, Theorem 3.82). These three examples thus describe the entire class of constant curvature spaces (up to covering and isometry).

To finish this section, we may now explain how curvature modifies the correspondence between distances in tangent spaces and Riemannian distances on a manifold. Indeed, let $(M, g)$ be a Riemannian manifolds and consider two orthogonal tangent vectors $v, w$ at some $x \in M$. In the vector space $T_x M$, the distance (induced by the metric at $x$) between $rv$ and $rw$ for some $r > 0$ is of course $r\sqrt{2}$. Now map these vectors to the manifold using the exponential map, then the Riemannian distance between $\text{Exp}_x rv$ and $\text{Exp}_x rw$ may increase or decrease compared to the distance in the tangent space. Curvature is the fundamental tool that allows to quantify these variations, and the sign of the sectional curvature tells whether geodesics accumulate or grow apart (Figure 3.2).

**Theorem 3.4.5.** Let $(M, g)$ be a Riemannian manifold, $R$ its curvature tensor, $\kappa$ its sectional curvature, $x \in M$ and $\gamma_v, \gamma_w$ two geodesics starting from $x$ with initial velocities $v, w \in T_x M$. Then for $r \to 0$,

$$d(\gamma_v(r), \gamma_w(r))^2 = r^2 \|v - w\|^2 - \frac{r^4}{3} \langle R(v, w)w, v \rangle + O(r^5). \quad (3.22)$$

Moreover, if $v, w$ are orthonormal

$$d(\gamma_v(r), \gamma_w(r)) = \sqrt{2}r \left(1 - \frac{\kappa(v, w)}{12}r^2\right) + O(r^4). \quad (3.23)$$
We now understand the meaning of the sign of the sectional curvature: if $\kappa > 0$, geodesics get closer to one another, while if $\kappa < 0$, they grow apart. Using this theorem and formulas for geodesics on the hypersphere in Example 3.6, we can compute a Taylor expansion of $d(\gamma_v, \gamma_w)$ and identify the coefficients, to find that $\kappa = 1$. Similarly, for the hyperbolic space $\kappa = -1$.

Finally, we state two important theorems that give a hint about how the sign of the curvature determines the geometry. These theorems also show how curvature and topology may be intertwined. See Lafontaine *et al.* (2004, Theorem 3.87) for a proof.

**Theorem 3.4.6 (Cartan-Hadamard).** Let $(M, g)$ be a complete connected Riemannian manifold with non-positive sectional curvature. Then the exponential map is a Riemannian covering, i.e. a covering that is a local isometry. This means that if $M$ is simply connected, then $M$ is diffeomorphic to $\mathbb{R}^d$, and any two points are joined by a unique minimizing geodesic.

**Remark 3.4.1 (Infinite dimension).** This theorem carries over to infinite dimensional Riemannian Hilbert manifolds, i.e. whose atlases are valued in a Hilbert space (McAlpin, 1965).

On the contrary (Berger, 2003, Theorem 63):

**Theorem 3.4.7 (Particular case of Bonnet-Myers theorem).** Let $(M, g)$ be a complete connected Riemannian manifold with sectional curvature bounded below $\kappa > \frac{1}{r^2}$ for some $r > 0$. Then the diameter of $M$ (i.e. the largest distance between points in $M$) is bounded above by $\pi r$ and $M$ is compact.

This ends our exposition of Riemannian metrics. In the next section, we focus on a particular class of manifolds that play a fundamental role in geometry and in many applications.

# 4 Lie groups

We first define the notions of Lie groups, subgroups and algebras, then the exponential map of Lie groups, and finally introduce smooth actions of Lie groups and homogeneous spaces. The metrics that can be defined on Lie groups, and their implementation in `geomstats` are covered in Section 4.3. We restrict to finite-dimensional Lie groups and only give a few remarks about infinite dimension.
## 4.1 Lie groups, Lie algebras and Lie subgroups

**Definition 4.1.1 (Lie group).** A Lie group is a group* $(G, \cdot)$ such that $G$ is also a finite dimensional smooth manifold, and the group and differential structures are compatible, in the sense that the group law $\cdot$ and the inverse map $g \mapsto g^{-1}$ are smooth.

We generally omit $\cdot$ in the notation and simply write the group composition as a multiplication. Let $e$ denote the neutral element, or identity of $G$.

**Remark 4.1.1 (Infinite dimension).** Infinite dimensional groups with a smooth structure appear naturally, for example as the set of diffeomorphisms of a smooth manifolds. However, although some properties of finite dimensional Lie groups generalize to infinite dimension, many don't. See Milnor (1984) for a good treatment of the topic. Unless otherwise specified, we restrict our exposition to finite dimension.

For any $g \in G$, we define the left and right translation maps $L_g$ and $R_g$ by

$$L_g : h \in G \mapsto gh \in G \qquad R_g : h \in G \mapsto hg^{-1}.$$

By the definition of a Lie group, $L_g$ and $R_g$ are diffeomorphisms of $G$ and their differential are linear isomorphisms of tangent spaces. This means that all the tangent spaces of $G$ are identical to $T_eG$, and there is a canonical way of mapping vectors from $T_gG$ to $T_eG$, by $L_{g^{-1}}$ for any $g \in G$.

This fact, that could be rephrased as “the tangent bundle of a Lie group is trivial”, means that $TG$ is diffeomorphic to the product $G \times T_eG$, and is of fundamental importance to characterize the geometry of a Lie group, and is very handy for the implementation. Furthermore, the maps $g \mapsto L_g$ and $g \mapsto R_g$ are group homomorphisms* between $(G, \cdot)$ and $(\mathrm{Diff}(G), \circ)$, the group of diffeomorphisms of $G$ with the composition of maps as group law.

**Definition 4.1.2 (Invariant vector field).** A vector field $X \in \Gamma(TG)$ is left-invariant if $\forall g, h \in G$,

$$X_{gh} = dL_g X_h$$

Let $\mathcal{L}(G)$ denote the set of all left-invariant vector fields.

### Example 4.1: General Linear group

The real general linear group $GL(n)$ plays a fundamental role as all the groups implemented in `geomstats` and encountered in the applications are subgroups of $GL(n)$. It is defined as the set of invertible matrices of size $n$

$$GL(n) = \{A \in M_n(\mathbb{R}) | \det(A) \neq 0\},$$

<footnote>*Defined in Appendix A.</footnote>
and as an open set of the vector space $M_n(\mathbb{R}) \simeq \mathbb{R}^{n^2}$ it is a smooth manifold.

The group law is the matrix multiplication, that can be written as a polynomial of the matrices coefficients and is thus smooth. Thanks to Cramer's formula $A^{-1} = \det(A)^{-1}Co(A)^T$, where $Co$ denotes the matrix formed by all cofactors, the inversion map is also smooth. Therefore, $GL(n)$ is a Lie group.

**Definition 4.1.3 (Lie algebra).** A real Lie algebra is a real vector space $\mathfrak{g}$ equipped with a bilinear map $[\cdot, \cdot]$ that verifies

$$ \bullet \text{ (Skew-symmetry)} \quad \forall x, y \in \mathfrak{g}, \quad [x, y] = -[y, x], $$

$$ \bullet \text{ (Jacobi identity)} \quad \forall x, y \in \mathfrak{g}, \quad [x, [y, z]] + [y, [z, x]] + [z, [y, x]] = 0. $$

### Example 4.2: Vector fields

The space of vector fields $\Gamma(TM)$ of a manifold $M$ equipped with the Lie bracket of vector field defined in Equation (2.2) page 15 is an infinite dimensional Lie algebra.

### Example 4.3: Matrix algebra

The algebra of squared matrices $M_n(\mathbb{R})$ equipped with the commutator

$$ \forall A, B \in M_n(\mathbb{R}), \quad [A, B] = AB - BA $$

is a Lie algebra.

**Implementation in geomstats** As the Lie groups we are working with are all matrix groups, we implement an abstract class `MatrixLieGroup` to gather the properties of a Lie group: its identity, composition law, translation map, exponential and logarithm maps. Moreover, it is useful to have the Lie algebra implemented as a separate class and set as an attribute of the Lie group. This class is a subclass of `VectorSpace` and implements a `belongs` and `projection` method. This allows for example to write the methods `is_tangent` and `to_tangent` as follows:

```python
def tangentaramation(self, point, left_or_right='left', inverse=False):
    """Return the differential map of the right or left translation map."""
    point_ = selfpelverse(point) if inverse else point
    if left_or_right == 'left':
        return lambda tan_vec: selfcompose(point_, tan_vec)
    return lambda tan_vec: selfcompose(tan_vec, point_)

def is_tangent(self, vector, base_point=None, atol=gs.atol):
    """Check whether a vector is tangent at base point."""
```
```python
if base_point is None:
    base_point = self的身份

if gs.allclose(base_point, self的身份):
    tangent_vec_at_id = vector
else:
    tangent_vec_at_id = self.tangentrimulation
    base_point, inverse = True)(vector)
return self.lie_algebra.belongs(tangent_vec_at_id, atol)

def to_tangent(self, vector, base_point=None):
    """Project a vector onto the tangent space at a base point."""
    if base_point is None:
        return self.lie_algebra.projection(vector)
    tangent_vec_at_id = self.tangentrimulation
    base_point, inverse = True)(vector)
    projected = self.lie_algebra.projection(tangent_vec_at_id)
    return self.tangentstranslation(base_point)(projected)
```

### Example 4.4: Implementation of GL(n)

The general linear group (Example 4.1) is the prototype of matrix Lie group, and thus created as a subclass of MatrixLieGroup. It is defined as an open set of the matrix space, so it also inherits from OpenSet.

```python
class GeneralLinear (MatrixLieGroup, OpenSet):
    """Class for the general linear group GL(n)."""

    def __init__ (self, n, **kwargs):
        if 'dim' not in kwargs.keys():
            kwargs['dim'] = n ** 2
        super().__init__ (
            ambient_space = Matrices(n, n), n=n, **kwargs)

    def belongs(self, point, atol = gs.atol):
        """Check if a matrix is invertible and of size n."""
        has_right_size = self.ambient_space.belongs(point)
        if gs.all (has_right_size):
            det = gs.linalg. det (point)
            return gs.abs (det) > atol
        return has_right_size
```

Again this is an open set, so the projection performs a regularization by adding $\epsilon I_n$:

```python
def projection(self, point):
    """Project a matrix to the general linear group."""
    belongs = self.belongs(point)
    regularization = gs.einsum(
        '...ij->...ij', gs. where (~belongs, gs.atol, 0.), self身份)
    projected = point + regularization
    return projected
```

The set of left-invariant vector fields is fundamental when studying Lie groups it is closed under Lie brackets and boils down to the tangent space at the identity, as stated by the following theorem (Lee, 2003, Proposition 8.33 and Theorem 8.37).
**Theorem 4.1.1 (Left-invariant vector fields).** Let $G$ be a Lie group of finite dimension $d$.

(1) The space of left-invariant vector fields $\mathcal{L}(G)$ of a Lie group $G$ is a sub-algebra of $\Gamma(TG)$.
This means that the bracket of two left-invariant vector fields is also left-invariant.

(2) As a vector space, $\mathcal{L}(G)$ is isomorphic to $T_eG$, the tangent space at the identity of $G$.
This implies that $\mathcal{L}(G)$ is finite dimensional.

Indeed, the map $X \in \mathcal{L}(G) \mapsto X_e \in T_eG$ is linear and has inverse $x \mapsto \tilde{x}$ where we define $\tilde{x}$ to be the left-invariant vector field

$$ \tilde{x} : g \in G \mapsto dL_g x \in T_g G. $$

We can thus define the bracket on $T_eG: [x, y] \triangleq [\tilde{x}, \tilde{y}]_e$, which turns $T_eG$ into a Lie algebra that is isomorphic (as Lie algebras) to $\mathcal{L}(G)$. We call $T_eG$ the Lie algebra of $G$ and denote it $\mathfrak{g} = T_eG$.

Thankfully, a linear representation of $G$ allows to compute the Lie bracket without requiring the handling of vector fields. Indeed, define the conjugation map of $G$, for any $g \in G$ by $C_g : h \mapsto ghg^{-1}$, and its differential at the identity, called adjoint representation of $G$:

$$ \mathrm{AD}_g : \mathfrak{g} \mapsto \mathfrak{g}. $$

Then the differential of $g \mapsto \mathrm{AD}_g$, written ad is called the adjoint representation of $\mathfrak{g}$, and one can show that it coincides with the Lie bracket of $\mathfrak{g}$:

$$ \mathrm{ad}_x(y) = [x, y]. $$

**Example 4.5: General linear Lie algebra**

Recall from Example 4.1 that $GL(n)$ is a Lie group, whose differentiable structure comes from its embedding in $M_n(\mathbb{R})$. Thus, by (1) of Theorem 2.1.2 (page 12), its tangent space at the identity is the entire matrix space, and as the Lie algebra of $GL(n)$ it is written $\mathfrak{gl}(n)$. We now compute explicitely $\mathrm{AD}$ and $\mathrm{ad}$ and verify that the bracket of $\mathfrak{gl}(n)$ coincides with the commutator defined in Example 4.3.

The conjugation map is the restriction to $GL(n)$ of the linear map $h \mapsto ghg^{-1}$ of $M_n(\mathbb{R})$ for any $g \in GL(n)$. Its differential is thus for any $X \in M_n(\mathbb{R})$

$$ \mathrm{AD}_g(X) = gXg^{-1}. $$

Now consider a curve $c : (-\epsilon, \epsilon) \to GL(n)$ such that $c(0) = I_n$ and $c'(0) = X \in M_n(\mathbb{R})$.
For any $Y \in M_n(\mathbb{R})$ we have

$$ \mathrm{ad}_X Y = d(g \mapsto \mathrm{AD}_g(Y))_{I_n} = \left. \frac{d}{dt} \right|_{t=0} (c(t)Yc(t)^{-1}) = XY - YX $$
To summarize,

* The Lie algebra of a Lie group is its tangent space at identity $T_eG$ or equivalently, the algebra of left-invariant vector fields.

* The Lie bracket on $T_eG$ is defined via the adjoint representation and coincides with that of left-invariant vector fields.

* In the case of linear Lie algebras (i.e. subalgebras of $M_n(\mathbb{R})$), the Lie bracket coincides with the matrix commutator.

**Remark 4.1.2.** Note that right-invariant vector fields can be defined analogously as left-invariant vector fields, and their set forms a Lie algebra. However, the Lie bracket of right-invariant vector fields coincides with the opposite of the adjoint representation. For practical and historical reasons, right-invariant fields are used in infinite dimension while left-invariant fields are used in finite dimension.

As all the next examples will be subgroups of the general linear group, we state a result by von Neumann and Cartan known as the *closed-subgroup theorem* that gives a necessary and sufficient condition for a subgroup of a Lie group $G$ to be an embedded Lie subgroup, i.e. a Lie group with differential structure agreeing with that of $G$, and such that the inclusion map is smooth (Cartan, 1930, Section III, paragraph 26). See Lee (2003, Theorem 20.12) for a proof.

**Theorem 4.1.2 (Cartan-von-Neumann).** Let $G$ be a Lie group. Any closed subgroup $H \subset G$ of $G$ is a Lie subgroup of $G$. Conversely, any Lie subgroup of $G$ is closed.

One can also show that the Lie algebra of a subgroup $G$ is a Lie subalgebra of the Lie algebra $\mathfrak{g}$ of $G$. Thus, they share the same bracket. Therefore, all Lie algebras considered in *geomstats* use the matrix commutator as bracket.

### Example 4.6: Implementation of $SO(n)$

Recall that the special orthogonal group is the set of orthogonal matrices of positive determinant. It is clear that it is stable by multiplication, hence it is a subgroup of $GL(n)$.

Furthermore, $SO(n) = f^{-1}(I_n) \cap \det^{-1}(1)$, where $f: A \mapsto A^\top A$ and $\det$ are continuous maps. $SO(n)$ is thus closed in $M_n(\mathbb{R})$, and is hence a Lie subgroup of $GL(n)$.

Moreover, from Example 2.11 (page 13), its Lie algebra is

$$ \mathfrak{so}(n) \triangleq T_{I_n} \mathfrak{SO}(n) = \{A \in M_n(\mathbb{R}) | A^\top + A = 0\} = \text{Skew}(n). $$

In *geomstats*, $SO(n)$ is implemented as a subclass of the `MatrixLieGroup` class and as
it is embedded in $GL_+(n)$ (the subgroup of $GL(n)$ formed by matrices with positive determinant), it also inherits from the `LevelSet` class. This allows to inherit the composition method, but the inverse method is overridden by the transposition.

```python
class SpecialOrthogonalMatrices (MatrixLieGroup, LevelSet):
    """Class for special orthogonal group."""

    def __init__ (self, n):
        matrices = Matrices(n, n)
        gln = GeneralLinear(n, positive Det = True)
        super().__init__ (
            dim=int((n * (n - 1)) / 2), n=n, value=gs. eye(n),
            lie_algebra=SkewSymmetricMatrices(n=n), embedding_space=gln,
            submersion=lambda x: matrices.mul matr ices transpose(x), x),
            tangent_submersion=lambda v, x: 2 * matrices. to symmetric (
                matrices.mul (matrices. transpose(x), v))
        self.bi_invariant metric = Bi InvariantMetric(group = self)
        self .metric = self.bi_invariant metric

    @classmethod
    def inverse(cls, point):
        """Return the transpose matrix of point."""
        return cls. transpose(point)

    def projection(self, point):
        """Project a matrix on SO(n) by minimizing the Frobenius norm."""
        aux mat = self.submersion(point)
        inv sqrt mat = SymmetricMatrices. powerm (aux mat, - 1 / 2)
        rotation mat = Matrices.mul (point, inv sqrt mat)
        det = gs. linalg. det (rotation mat)
        return utils .flip determinant (rotation mat, det)
```

And the Lie algebra is implemented as follows.

```python
class SkewSymmetricMatrices (MatrixLieAlgebra):
    """Class for skew-symmetric matrices."""

    def __init__ (self, n):
        dim = int(n * (n - 1) / 2)
        super().__init__ (dim, n)
        self.ambient_space = Matrices(n, n)

    def belongs(self, mat, atol=gs.atol):
        """Evaluate if mat is a skew-symmetric matrix."""
        has_rightShape = self.ambient_space.belongs (mat)
        if has_rightShape:
            return Matrices. equal (mat, - Matrices. transpose (mat), atol=atol)
        return False

    @classmethod
    def projection(cls, mat):
        """Compute the skew-symmetric component of a matrix."""
        return 1 / 2 * (mat - Matrices. transpose (mat))
```
### Example 4.7: Implementation of SE(n)

We now introduce a group that is ubiquitous in applications. This group is defined as the set of direct isometries - or rigid-body transformations - of $\mathbb{R}^n$, i.e. the linear transformations of the affine space $\mathbb{R}^n$ that preserve its canonical inner-product. Such transformation $\rho$ can be decomposed in a rotation part and a translation part: $\rho(x) = Rx + u$, where $R \in SO(n)$ and $x, u \in \mathbb{R}^n$. Define

$$SE(n) = \{(R, u) | R \in SO(n), u \in \mathbb{R}^n\}.$$

Now, the composition of two isometries $\rho = (R, u), \rho'(R', u')$ remains an isometry:

$$\rho \circ \rho'(x) = RR'x + Ru' + u,$$

and can be written $\rho \circ \rho' = (RR', Ru' + u)$. This suggests the representation of $SE(n)$ in homogeneous coordinates

$$\rho = \begin{pmatrix} R & u \\ 0 & 1 \end{pmatrix} \in GL(n+1) \qquad (4.1)$$

The composition of isometries then corresponds to matrix multiplication, and $SE(n)$ is now a Lie subgroup of $GL(n+1)$. Its Lie algebra is

$$\mathfrak{se}(n) = \left\{ \begin{pmatrix} S & v \\ 0 & 0 \end{pmatrix} \mid S \in \text{Skew}(n), v \in \mathbb{R}^n \right\}.$$

Thus its dimension is $\frac{n(n+1)}{2}$. In geomstats, $SE(n)$ inherits from MatrixLieGroup and is embedded in $GL_+(n+1)$ with overridden inverse function. We used a utility function that builds a block matrix according to Equation (4.1) The submersion used is the map

$$\begin{pmatrix} R & u \\ v^\top & c \end{pmatrix} \in GL_+(n+1) \mapsto (R^\top R, v, c) \in M_n(\mathbb{R}) \times \mathbb{R}^n \times \mathbb{R},$$

so that $SE(n) = f^{-1}\{(I_n, 0, 1)\}$. It is not printed in the code below in the interest of space.
class SpecialEuclideanMatrices (MatrixLieGroup, LevelSet):
    """Class for special Euclidean group."""
    def __init__ (self, n):
        super().__init__ (
            n=n + 1, dim=int((n * (n + 1)) / 2),
            embedding_space=GeneralLinear(n + 1, positiveTrue),
            submersion =submersion, value=gs. eye(n + 1),
            tangent_submersion=tangent_submersion,
            lie_algebra=SpecialEuclideanMatrixLieAlgebra(n=n))
        self. rotations =SpecialOrthogonal(n=n)
        self. translations =Euclidean(dim=n)
        self. n =n
        self. left_canonical metric = \n
        SpecialEuclideanMatrixCannonicalLeftMetric(group =self)
        self. metric =self. left_canonical metric
    @property
    def identity(self):
        """Return the identity matrix."""
        return gs. eye(self. n + 1, self. n + 1)
    def inverse(self, point):
        """Return the inverse of a point."""
        n =self. n
        transposed_rot =self. transpose(point [...] ,:n, :n])
        translation =point[...] ,:n, -1]
        translation =gs. einsum(
            '...ij,...j', transposed_rot, translation)
        return homogeneous_representation(
            transposed_rot, -translation, point. shape)

**Remark 4.1.3.** Note that as a manifold, *SE*(n) defined in the above example corresponds to the product manifold *SO*(n) × ℝⁿ (Example 2.6 page 9), but as a Lie group, it is the semi-direct product *SO*(n) × ℝⁿ, because the rotation part acts on the translation part in the composition rule.

## 4.2 The exponential map

We now study the flow of left-invariant vector fields of a Lie group. For more details and proofs, we refer the reader to Gallier and Quaintance (2020, Chapter 18). It allows to canonically map the Lie algebra to its Lie group. First, the left-invariance translates into a commutation property of the flow with the left translation map:

**Proposition 4.2.1.** Let $G$ be a Lie group and $X \in \mathcal{L}(G)$. Then $X$ is complete, and if $(\phi_t)_t$ is its flow, then for any $t \in \mathbb{R}$ and $g, g' \in G$ we have

$$ \phi_t(gg') = g\phi_t(g'), \quad \text{i.e.} \quad \phi_t \circ L_g = L_g \circ \phi_t. $$

This allows the following:

**Definition 4.2.1 (Exponential map).** We call exponential map, and write $\exp : \mathfrak{g} \to G$ the map defined by $x \mapsto \phi_1(e)$ where $(\phi_t)_t$ is the flow of the left-invariant vector fields $\tilde{x}$.
**Remark 4.2.1.** This is the second definition of exponential map we encounter, where the first was the Riemannian exponential defined by the geodesic flow of a metric. We will refer to the one canonically defined on a Lie group as the *group exponential* and use a lowercase exp.

**Example 4.8: Exponential map of GL(n)**

The fundamental case is of course again that of the general linear group. Recall that its Lie algebra $\mathfrak{gl}(n) = M_n(\mathbb{R})$ is the set of square matrices, and the group law is the (linear) matrix multiplication. Therefore, a left-invariant vector field $\tilde{X}$ associated to $X \in \mathfrak{gl}(n)$ is defined by

$$g \mapsto gX$$

Let $\gamma$ be the integral curve from the identity $I_n$. This means that $\gamma$ is solution to the ODE defined on $\mathbb{R}$ with initial condition $\gamma(0) = I_n$

$$\gamma'(t) = \gamma(t)X$$

It is well known that the unique solution is the matrix exponential, defined by the series

$$e^{tX} = \sum_{k=0}^{\infty} \frac{t^k}{k!} X^k.$$

Thus $\exp(X) = e^X$.

We now state some of the fundamental properties of the exponential map (see e.g. Gallier and Quaintance, 2020, Proposition 18.6-18.7-18.13).

**Proposition 4.2.2.** For any Lie group $G$, the exponential map $\exp : \mathfrak{g} \to G$ is smooth and a local diffeomorphism at 0.

The inverse of the group exponential, defined locally, is called *logarithm map* and is valued in $\mathfrak{g} = T_eG$. It thus allows to map data defined on $G$ to its Lie algebra $\mathfrak{g}$, which is a vector space! This fact is at the basis of many algorithms that handle Lie group data.

**Remark 4.2.2 (Infinite dimension).** Proposition 4.2.2 is not true in infinite dimension. For example, if $G = \text{Diff}(M)$, the exp map is not even surjective in any neighborhood of the identity (Schmid, 2004).

**Proposition 4.2.3.** Let $G, H$ be a Lie group, and $f : G \to H$ a Lie group homomorphism. Then

$$f \circ \exp = \exp \circ df_e \qquad (4.2)$$

In particular, if $G = H$ and $f = C_g$, we have

$$\exp(t \mathrm{AD}_g(u)) = g \exp(tu) g^{-1} = C_g(\exp(tu)). \qquad (4.3)$$
See Lee (2003, Proposition 20.8) for a proof. The commutation property (4.2) can be depicted by saying that the following diagram commutes, meaning both paths leading from the top left space to the bottom right one are equivalent:

<figure><img src="image_13.png" /></figure>

Using $f = L_g$, or equivalently (thanks to Equation (4.3)) $f = R_g$, we define the exponential map at any point $g \in G$ so that it fulfills the above properties. For any $v \in T_g G$

$$ \exp_g(v) = L_g \circ \exp \circ (dL_g)^{-1} = R_{g^{-1}} \circ \exp \circ dR_g. $$

where differentiation is at $e$, and so that the group exp is now understood as $\exp_e$. In `geomstats` we use the implementation of the matrix exponential and logarithm from the backends and include it in the `MatrixLieGroup` class.

Moreover, we can recover the whole connected component of the identity from the Lie algebra:

**Theorem 4.2.1.** If $G$ is a Lie group and $G_0$ is the connected component of $e$, then $G_0$ is generated by $\exp(g)$.

In particular (see Gallier and Quaintance (2020, Theorem 1.6 and 1.12)), we have the following

**Proposition 4.2.4.** SO(n) and SE(n)

* The exponential map $\exp : \mathfrak{so}(n) \to SO(n)$ is surjective,

* The exponential map $\exp : \mathfrak{sc}(n) \to SE(n)$ is surjective.

Note that by applying the commutation property (4.2) to the inclusion map, the exponential maps of $SE(n)$ and $SO(n)$ coincide with the restriction of the exponential map of $GL(n)$, i.e. the matrix exponential.

Finally, recall that flows of complete vector fields are one-parameters subgroups of $Diff(M)$. Define now a one-parameter subgroup of $G$ as a Lie group homomorphism $t \in \mathbb{R} \mapsto \gamma_t \in G$, written $(\gamma_t)_{t \in \mathbb{R}}$. Then by a change of variable in the flow equation, it is clear that $t \mapsto \exp(tX)$ is a one-parameter subgroup for any $X \in g$. In fact, all one-parameter subgroups are of this form and the maps:

$$ \theta : (\gamma_t)_{t \in \mathbb{R}} \mapsto \frac{d}{dt}\Big|_{t=0} \gamma_t; \quad \theta' : X \mapsto (\exp(tX))_{t \in \mathbb{R}} $$

are inverse to each other, between $g$ and the set of one-parameter subgroups.

Before closing this section, we illustrate it with a representation of the one-parameter subgroups of $SE(2)$.
### Example 4.9: Curve in $SE(2)$

Recall that $SE(2)$ is the group of rotation-translation transformations, so that a smooth curve of $SE(2)$ can represent the motion of a rigid-body in a plane. An element of $SE(2)$ can be represented by the application of the rotation part to the canonical orthonormal frame of the 2d-plane, while the translation part is canonically mapped to a point in the 2d-plane. We use this representation in Figure 4.1 to show one-parameter subgroups. See Appendix B.3 for the code and different initial vectors.

<figure><img src="image_14.png" /><figcaption>**Figure 4.1:** One-parameter subgroup of $SE(2)$ for three different initial conditions (that share the same rotation part). The translation part corresponds to a circle, while the rotation occurs with constant angular velocity.</figcaption></figure>

## 4.3 Invariant metrics on Lie groups

### 4.3.1 Rationale

Different geometric structures are compatible with the group structure. For instance, the group exponential defined in Section 4 can be defined as the exponential map of a connection on $G$, implying that geodesics are one-parameter subgroups. This connection is known as the canonical Cartan connection, and is of practical interest as in linear groups geodesics can be computed in closed-form using the matrix exponential.

The canonical Cartan connection may however not be the Levi-Civita connection of any metric, for more details we refer the reader to Pennec *et al.* (2020, Chapter 5). In this section we focus on the study of invariant metrics on Lie groups, i.e. metrics for which the left (or right) translation map is an isometry. This case is fundamental in geometric
mechanics (Kolev, 2004) and has been studied in depth since the foundational papers of Arnold (1966) and Milnor (1976). Using left-invariant vector fields, one can compute explicitly the LC connection, allowing to rewrite the geodesic equation. This fundamental idea, known as Euler-Poincaré reduction, is that the geodesic equation can be expressed entirely in the Lie algebra thanks to the symmetry of left-invariance (Marsden and Ratiu, 2009), alleviating the burden of coordinate charts.

We derive here a similar reduction of the parallel transport equation, resulting in a stable and efficient implementation of parallel transport in geomstats. Finally, knowing the connection, one can also compute algebraically the curvature of the space. We exemplify these results with an anisotropic metric on the special Euclidean groups SE(2) and SE(3). Part of the material and results of this section were presented at the GSI 2021 conference in Guigui and Pennec (2021).

### 4.3.2 Definitions and computations of the connection

Let $G$ be a Lie group. Recall that the Lie algebra of left-invariant vector fields $\mathcal{L}(G)$ and the tangent space at the identity $\mathfrak{g}$ can be identified, and in this section we will write $\tilde{x}$ the left-invariant field generated by $x \in \mathfrak{g}$: $\forall g \in G, \tilde{x}_g = dL_g x$.

**Definition 4.3.1 (Invariant metric).** Let $G$ be a Lie group. A Riemannian metric $\langle \cdot, \cdot \rangle$ on $G$ is called *left-invariant* if the differential map of the left translation is an isometry between tangent spaces, that is

$$ \forall g, h \in G, \forall u, v \in T_g G, \quad \langle u, v \rangle_g = \langle dL_h u, dL_h v \rangle_{hg}. $$

It is thus uniquely determined by an inner product on the tangent space at the identity $T_e G = \mathfrak{g}$ of $G$. Similarly, a *right-invariant* metric is such that the differential of the right translation is an isometry. A metric is called *bi-invariant* if it is both left and right-invariant.

Moreover, let $(e_1, \dots, e_n)$ be an orthonormal basis of $\mathfrak{g}$, and the associated left-invariant vector fields $\tilde{e}_i$. As $dL_g$ is an isomorphism, $(\tilde{e}_{1,g}, \dots, \tilde{e}_{n,g})$ form a basis of $T_g G$ for any $g \in G$, so for any $X \in \Gamma(G)$ one can write $X_g = f^i(g)\tilde{e}_{i,g}$ where for $i = 1, \dots, n$ $g \mapsto f^i(g)$ is a smooth real-valued function on $G$. Any vector field on $G$ can thus be expressed as a linear combination of the $\tilde{e}_i$ with functional coefficients.

#### Example 4.10: Invariance metric on $SO(n)$

Consider the special orthogonal group endowed with the restriction of the Frobenius metric of Example 3.5. Let $P, Q \in SO(n)$ and $A, B \in T_P SO(n) = P \text{Skew}(n)$. Then it is clear that

$$ \langle dL_Q A, dL_Q B \rangle = \text{tr}((QA)^\top QB) = \text{tr}(A^\top B) = \langle A, B \rangle, $$

and similarly that $\langle dR_Q A, dR_Q B \rangle = \langle A, B \rangle$. Therefore, the Frobenius metric is
bi-invariant on $SO(n)$.

**Example 4.11: Invariance metric on $SE(n)$**

Now, consider the special Euclidean group endowed with the restriction of the Frobenius metric. As in the previous example, it is easy to show that it is left-invariant. However, it is not right-invariant. Indeed, let $g = (P, t) \in SE(n)$ and $x = (S, u), y = (T, v) \in \mathfrak{sc}(n) = \text{Skew}(n) \oplus \mathbb{R}^n$. Then

$$ \begin{aligned} \langle dR_g x, dR_g y \rangle &= \langle \begin{pmatrix} SP^\top & u - SP^\top t \\ 0 & 0 \end{pmatrix}, \begin{pmatrix} TP^\top & v - TP^\top t \\ 0 & 0 \end{pmatrix} \rangle \\ &= \text{tr}(S^\top T) + u^T v + t^\top P S^\top T P^\top t \\ &= \langle x, y \rangle + t^\top P S^\top T P^\top t \end{aligned} $$

So that one can find $t \neq 0$ and $S, T \in \text{Skew}(n)$ such that $\langle dR_g x, dR_g y \rangle \neq \langle x, y \rangle$. Therefore the metric is not bi-invariant. In fact one can show that there does not exist any bi-invariant metric on $SE(n)$ for $n \neq 1$, and there exists a bi-invariant pseudo-metric for $n = 3$ (Miolane and Pennec, 2015).

**Definition 4.3.2 (Dual adjoint).** Define the metric dual adjoint map on $\mathfrak{g}$ as the unique map that verifies

$$ \forall a, b, c \in \mathfrak{g}, \langle \mathrm{ad}_a^*(b), c \rangle = \langle b, \mathrm{ad}_a(c) \rangle = \langle [a, c], b \rangle. $$

As the bracket can be computed explicitly in the Lie algebra, so can $\mathrm{ad}^*$ thanks to the orthonormal basis of $\mathfrak{g}$. Now let $\nabla$ be the Levi-Civita connection associated to the metric. It is also left-invariant and can be characterized by a bi-linear form on $\mathfrak{g}$ that verifies $\forall x, y \in \mathfrak{g}$ (Pennec and Arsigny, 2013; Gallier and Quaintance, 2020):

$$ \alpha(x, y) \triangleq (\nabla_{\tilde{x}} \tilde{y})_e = \frac{1}{2} ([x, y] - \mathrm{ad}_x^*(y) - \mathrm{ad}_y^*(x)) \quad (4.4) $$

Indeed by the left-invariance, for two left-invariant vector fields $X = \tilde{x}, Y = \tilde{y} \in \mathcal{L}(G)$, the map $g \mapsto \langle X, Y \rangle_g$ is constant, so for any vector field $Z = \tilde{z}$ we have $Z(\langle X, Y \rangle) = 0$. Kozsul formula thus becomes

$$ \begin{aligned} 2\langle \nabla_X Y, Z \rangle &= \langle [X, Y], Z \rangle - \langle [Y, Z], X \rangle - \langle [X, Z], Y \rangle \\ 2\langle \nabla_X Y, Z \rangle_e &= \langle [x, y], z \rangle_e - \langle \mathrm{ad}_y(z), x \rangle_e - \langle \mathrm{ad}_x(z), y \rangle_e \\ 2\langle \alpha(x, y), z \rangle_e &= \langle [x, y], z \rangle_e - \langle \mathrm{ad}_y^*(x), z \rangle_e - \langle \mathrm{ad}_x^*(y), z \rangle_e. \end{aligned} \quad (4.5) $$

**Definition 4.3.3 (Structure constants).** Let $G$ be a Lie group and $g = (\cdot, \cdot)$ be a left-invariant metric on $G$. Let $(e_1, \dots, e_d)$ be an orthonormal basis of $\mathfrak{g}$ for $g$. Define the structure constants
as

$$C_{ij}^{k} = \langle [e_i, e_j], e_k \rangle.$$

### Example 4.12: Structure constants on SO(3)

An orthonormal basis of $\mathfrak{so}(3)$ endowed with the Frobenius metric is $\frac{1}{\sqrt{2}}(a_1, a_2, a_3)$ with

$$a_1 = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & -1 \\ 0 & 1 & 0 \end{pmatrix}$$

$$a_2 = \begin{pmatrix} 0 & 0 & 1 \\ 0 & 0 & 0 \\ -1 & 0 & 0 \end{pmatrix}$$

$$a_3 = \begin{pmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}.$$

And as the Lie bracket is the matrix commutator, we can compute the structure constants to obtain

$$C_{ij}^{k} = \frac{1}{\sqrt{2}} \text{ if } ijk \text{ is a direct cycle of } \{1, 2, 3\},$$

and 0 otherwise.

### Example 4.13: Structure constants on SE(3)

We define a left-invariant metric on the special Euclidean group by defining an inner product in its Lie algebra. Let the metric matrix at the identity be diagonal: $g = \text{diag}(1, 1, 1, \beta, 1, 1)$ for some $\beta > 0$, an anisotropy parameter. For $\beta = 1$, this metric coincides with the Frobenius metric. An orthonormal basis of the Lie algebra $\mathfrak{se}(3)$ for this metric is

$$e_1 = \begin{pmatrix} a_1 & 0 \\ 0 & 0 \end{pmatrix}$$ $$e_4 = \frac{1}{\sqrt{\beta}} \begin{pmatrix} 0 & \epsilon_1 \\ 0 & 0 \end{pmatrix}$$

$$e_2 = \begin{pmatrix} a_2 & 0 \\ 0 & 0 \end{pmatrix}$$ $$e_5 = \begin{pmatrix} 0 & \epsilon_2 \\ 0 & 0 \end{pmatrix}$$

$$e_3 = \begin{pmatrix} a_3 & 0 \\ 0 & 0 \end{pmatrix}$$ $$e_6 = \begin{pmatrix} 0 & \epsilon_3 \\ 0 & 0 \end{pmatrix}.$$

where $(a_1, a_2, a_3)$ is the basis of $\mathfrak{so}(3)$ defined above, and $(\epsilon_1, \epsilon_2, \epsilon_3)$ is the canonical basis of $\mathbb{R}^3$. As the Lie bracket is the usual matrix commutator, it is straightforward to compute

$$C_{ij}^{k} = \frac{1}{\sqrt{2}} \text{ if } ijk \text{ is a direct cycle of } \{1, 2, 3\}; \qquad (4.6)$$

$$C_{15}^{6} = -C_{16}^{5} = -\sqrt{\beta}C_{24}^{6} = \frac{1}{\sqrt{\beta}}C_{26}^{4} = \sqrt{\beta}C_{34}^{5} = -\frac{1}{\sqrt{\beta}}C_{35}^{4} = \frac{1}{\sqrt{2}}. \qquad (4.7)$$
and all others that cannot be deduced by skew-symmetry of the bracket are equal to 0.

The structure constants allow to compute ad* and formula (4.4) in practice:

$$ \alpha(e_i, e_j) = \nabla_{e_i} e_j = \frac{1}{2} \sum_k (C_{ij}^k - C_{jk}^i + C_{ki}^j) e_k. $$

Note however that Equation (4.4) gives the connection for left-invariant vector fields only. We will now generalize to any vector field defined along a smooth curve on G, using the left-invariant basis $(\tilde{e}_1, \dots, \tilde{e}_n)$.

Let $\gamma: [0, 1] \to G$ be a smooth curve, and $Y$ a vector field defined along $\gamma$. Write $Y = h^i \tilde{e}_i$, $\dot{\gamma} = f^i \tilde{e}_i$. Let us also define the left-angular velocities $\omega = dL_{\gamma}^{-1} \dot{\gamma} = f^i \circ \gamma e_i \in \mathfrak{g}$ and $\zeta = dL_{\gamma}^{-1} Y = h^j \circ \gamma e_j \in \mathfrak{g}$. Then the covariant derivative of $Y$ along $\gamma$ is

$$ \begin{aligned} \nabla_{\dot{\gamma}(t)} Y &= (f^i \circ \gamma)(t) \nabla_{\tilde{e}_i} (h^j \tilde{e}_j) \\ &= (f^i \circ \gamma)(t) \tilde{e}_i (h^j) \tilde{e}_j + (f^i \circ \gamma)(t) (h^j \circ \gamma)(t) (\nabla_{\tilde{e}_i} \tilde{e}_j)_{\gamma(t)} \\ dL_{\gamma(t)}^{-1} \nabla_{\dot{\gamma}(t)} Y &= (f^i \circ \gamma)(t) \tilde{e}_i (h^j) e_j \\ &\quad + (f^i \circ \gamma)(t) (h^j \circ \gamma)(t) dL_{\gamma(t)}^{-1} (\nabla_{\tilde{e}_i} \tilde{e}_j)_{\gamma(t)} \\ &= (f^i \circ \gamma)(t) \tilde{e}_i (h^j) e_j + (f^i \circ \gamma)(t) (h^j \circ \gamma)(t) \nabla_{e_i} e_j \end{aligned} $$

where the Leibniz formula and the invariance of the connection is used in $(\nabla_{\tilde{e}_i} \tilde{e}_j) = dL_{\gamma(t)} \nabla_{e_i} e_j$. Therefore, for $k = 1..n$

$$ \langle dL_{\gamma(t)}^{-1} \nabla_{\dot{\gamma}(t)} Y, e_k \rangle = (f^i \circ \gamma)(t) \tilde{e}_i (h^j) \langle e_j, e_k \rangle \\ + (f^i \circ \gamma)(t) (h^j \circ \gamma)(t) \langle \nabla_{e_i} e_j, e_k \rangle $$

but on the one hand

$$ \begin{aligned} \zeta(t) &= (h^j \circ \gamma)(t) e_j \\ \dot{\zeta}(t) &= (h^j \circ \gamma)'(t) e_j = dh_{\gamma(t)}^j \dot{\gamma}(t) e_j \\ &= dh_{\gamma(t)}^j ((f^i \circ \gamma)(t) \tilde{e}_{i,\gamma(t)}) e_j \\ &= (f^i \circ \gamma)(t) d_{\gamma(t)} h^j \tilde{e}_{i,\gamma(t)} e_j \\ &= (f^i \circ \gamma)(t) \tilde{e}_i (h^j)_{\gamma(t)} e_j \end{aligned} $$
and on the other hand, using Equation (4.4):

$$
\begin{align*}
& 2(f^i \circ \gamma)(h^j \circ \gamma)\langle \nabla_{e_i} e_j, e_k \rangle \\
&= (f^i \circ \gamma)(h^j \circ \gamma)\left(\langle [e_i, e_j], e_k \rangle - \langle \mathrm{ad}_{e_i}^* e_j, e_k \rangle - \langle \mathrm{ad}_{e_j}^* e_i, e_k \rangle\right) \\
&= (f^i \circ \gamma)(h^j \circ \gamma)\left(\langle [e_i, e_j], e_k \rangle - \langle [e_j, e_k], e_i \rangle - \langle [e_i, e_k], e_j \rangle\right) \\
&= \left(\left(\langle [f^i \circ \gamma)e_i, (h^j \circ \gamma)e_j \rangle, e_k \rangle\right.\right. \\
&\qquad \left. - \langle \left(\langle [h^j \circ \gamma)e_j, e_k \rangle, (f^i \circ \gamma)e_i \rangle\right.\right. \\
&\qquad \left. - \langle \left(\langle [f^i \circ \gamma)e_i, e_k \rangle, (h^j \circ \gamma)e_j \rangle\right)\right)\right) \\
&= \langle \left(\left(\omega, \zeta\right] - \mathrm{ad}_\omega^* \zeta - \mathrm{ad}_\zeta^* \omega\right), e_k \rangle = \langle 2\alpha(\omega, \zeta), e_k \rangle.
\end{align*}
$$

Thus, we obtain an algebraic expression for the covariant derivative of any vector field Y along a smooth curve γ. It is the main ingredient of this section.

**Lemma 4.3.1.** Let $G$ be a Lie group and $\nabla$ be the Levi-Civita connection of a left-invariant metric on $G$. Let $\gamma$ be a smooth curve on $G$, and $Y$ a vector field defined along $\gamma$. Consider the left-angular velocities $\omega = dL_{\gamma}^{-1}\dot{\gamma}$ and $\zeta = dL_{\gamma}^{-1}Y$. Then

$$dL_{\gamma(t)}^{-1} \nabla_{\dot{\gamma}(t)} Y(t) = \dot{\zeta}(t) + \alpha(\omega(t), \zeta(t)) \quad (4.8)$$

A similar expression can be found in Arnold (1966) and Gay-Balmaz et al. (2012). As all the variables of the right-hand side are defined in g, they can be computed with matrix operations and an orthonormal basis.

We now focus on two particular cases of (4.8) to derive the equations of geodesics and of parallel transport along a curve.

### 4.3.3 Geodesic equation

The first particular case is for $Y(t) = \dot{\gamma}(t)$. It is then straightforward to deduce from Equation (4.8) the Euler-Poincaré equation for a geodesic curve (e.g Kolev, 2004; Cendra et al., 1998). Indeed in this case, recall that $\omega = dL_{\gamma}^{-1}\dot{\gamma}$ is the left-angular velocity, $\zeta = \omega$ and $\alpha(\omega, \omega) = -\mathrm{ad}_{\omega}^{*}(\omega)$. Hence $\gamma$ is a geodesic if and only if $dL_{\gamma(t)}^{-1}\nabla_{\dot{\gamma}(t)}\dot{\gamma}(t) = 0$ i.e. setting the left-hand side of Equation (4.8) to 0. We obtain

$$\left\{
\begin{array}{l}
\dot{\gamma}(t) = dL_{\gamma(t)}\omega(t) \\
\dot{\omega}(t) = \mathrm{ad}_{\omega(t)}^{*}\omega(t).
\end{array}
\right.
\qquad (4.9)$$

#### Remark 4.3.1.

* A similar treatment of a right-invariant metric is straightforward. Indeed, in this case the Lie bracket of right-invariant vector fields is the opposite of the adjoint representation of the Lie algebra, therefore the expressions are all the same with $-\alpha$ instead of $\alpha$.
* One can show that the metric is bi-invariant if and only if the adjoint map is skew-symmetric (see Pennec and Arsigny (2013) or Gallier and Quaintance (2020, Prop. 20.7)). In this case $ad_{\omega}^{*}(\omega) = 0$ and Equation (4.9) coincides with the equation of one-parameter subgroups on $G$.

In the class *InvariantMetric* of *geomstats*, the *geodesic_equation* method is modified to implement Equation (4.9), and the integrators of Section 3.3 (page 38) are used to compute the exponential map.

**Example 4.14: Geodesics on SE(2)**

We use the same visualization as in Example 4.9 (page 55) to plot the geodesics of a left-invariant metric on $SE(2)$. We first compare the left and right-invariant metrics with same inner-product at the identity, and compare their geodesics to one-parameter subgroups. These curves are shown on Figure 4.2.

<figure><img src="image_15.png" /><figcaption>**Figure 4.2:** Visualization of geodesics and one-parameter subgroups on $SE(2)$.</figcaption></figure>

Next, we define the metric as in Example 4.13, by the bilinear form defined in $\mathfrak{se}(2)$ by the matrix $g = \text{diag}(1, \beta, 1)$ for some $\beta > 0$, an anisotropy parameter, and extended by left-invariance.
For $\beta = 1$, the metric coincides with the Frobenius metric, and there is no interaction between the rotation and translation part. This metric is thus isomorphic to the direct product metric of the bi-invariant metric on $SO(2)$ and the usual inner product of $\mathbb{R}^2$. The geodesics are thus straight lines on the translation part, and rotations with constant speed (as one-parameter subgroups) on the rotation part.

For $\beta \neq 1$, the rotation part remains the same, but there is a direction towards which less movement occurs (Figure 4.3).

<figure><img src="image_16.png" /><figcaption>**Figure 4.3:** Visualization of geodesics on $SE(2)$ with the same initial conditions and different values of $\beta$.</figcaption></figure>
### 4.3.4 Reduced parallel transport equation

The second case is for a vector field $Y$ that is parallel along the curve $\gamma$, that is, $\forall t$, $\nabla_{\dot{\gamma}(t)}Y(t) = 0$. Similarly to the geodesic equation, we deduce from Equation (4.8) the parallel transport equation expressed in the Lie algebra. To the best of our knowledge, this formulation of the parallel transport on Lie groups with an invariant metric is original.

**Theorem 4.3.2.** Let $\gamma$ be a smooth curve on $G$. The vector $Y$ is parallel along $\gamma$ if and only if it is solution to

$$ egin{cases} 
\omega(t) = dL_{\gamma(t)}^{-1} \dot{\gamma}(t) \\ 
Y(t) = dL_{\gamma(t)} \zeta(t) \\ 
\dot{\zeta}(t) = -\alpha(\omega(t), \zeta(t)) 
\end{cases} 
\quad (4.10) $$

Note that in order to parallel transport along a geodesic curve, Equations (4.9) and (4.10) are solved jointly. We add the corresponding method to the `InvariantMetric` class of `geomstats`.

```python
def geodesic_equation(self, state, _time):
    """Compute the right-hand side of the geodesic equation."""
    sign = 1. if self.left_or_right == 'left' else -1.
    basis = self нормalbasis(self.lie_algebra.basis)

    point, vector = state
    velocity = self.group.tangent map(
        point, left_or_right = self. left_or_right (vector)
    )
    coefficients = gs.array([self. structure. constant(
        vector, basis_vector, vector) for basis_vector in basis])
    acceleration = gs.einsum('i... , ijk->...jk', coefficients, basis)
    return gs. stack([velocity, sign * acceleration])

def parallel_transport(
    self, tangentvec_a, tangent. vec_b, base_point, n_steps=10,
    step='rk4', **kwargs):
    """Compute the parallel transport of a tangent vector along a geodesic."""
    group = self.group
    translation_map = group.tangent map(
        base_point, left_or_right = self. left_or_right, inverse = True)
    left angular vel_a = group. to_tangent (translation_map(tangent. vec_a))
    left angular vel_b = group. to_tangent (translation_map(tangent. vec_b))

    def acceleration (state, time):
        point, omega, zeta = state
        gamdot, omegaDot = self. geodesic equation (state[:2], time)
        zetadot = - self. connection_atidentity (omega, zeta)
        return gs. stack([gamdot, omegaDot, zetadot])

    if (base_point.ndim == 2 or base_point.shape[0] == 1) and \\
        (tangent. vec_a. ndim == 3 or tangent. vec_b. ndim == 3):
        n_sample = tangent. vec_a. shape[0] if tangent. vec_a. ndim == 3 else \\
            tangent. vec_b. shape[0]
        base_point = gs. stack([base_point] * n_sample)

    initial_state = gs. stack([base_point, left angular vel_b, left angular vel_a])
    flow = integrate (acceleration, initial_state, n_steps = n_steps, step = step, **kwargs)
    gamma, gammaDot, zeta_t = flow[-1]
    transported = group.tangent map(
        gamma, left_or_right = self. left_or_right, inverse = False) (zeta_t)
    return transported
```
### Example 4.15: Parallel transport in SE(3)

We exemplify the use of the reduced parallel transport equation on SE(3) endowed with an anisotropic left-invariant metric, as in Examples 4.13 and 4.14. The metric is defined by the matrix $g = \text{diag}(1, 1, 1, \beta, 1, 1)$ for some $\beta > 0$ at the identity and extended by left-invariance.

We randomly generate a point in $x \in SE(3)$ and two tangent vectors $v, w \in T_xSE(3)$, and transport the vector $v$ along the geodesic $t \mapsto \text{Exp}_x(tw)$.

The results are plotted in a loglog plot (Figure 4.4), with the error shown with respect to the number of steps used to integrate Equations (4.9) and (4.10), with an RK2 or RK4 scheme. As expected, the speed of convergence depends on the order of the scheme used: quadratic speed is reached for RK2 and speed of order 4 for RK4.

<figure><img src="image_17.png" /><figcaption>**Figure 4.4:** Norm of the absolute error represented with respect to the number of steps.</figcaption></figure>

When $\beta = 1$, it coincides with the product metric between the bi-invariant metric of SO(3) and of $\mathbb{R}^3$, so the geodesics and parallel transport are known in closed form. For $\beta > 1$, the curvature grows away from 0 (this is proved in Appendix B), and so does its covariant derivative, so that parallel transport becomes more complex to compute. We use $n = 1100$ steps in the discretization of $[0, 1]$ for integration scheme to compute a reference value that is then used to measure the error of the method for $10 \le n \le 1000$. We use different values of $\beta$ but keep the same initial vectors
regardless of the value of $\beta$.

### 4.3.5 Curvature

By definition of the left-invariance of the metric, for any $g \in G$, the left-translation map $L_g$ is an isometry. It is thus sufficient to compute the curvature tensor at the identity, that is, in the Lie algebra $\mathfrak{g}$, and to map it to any other point. Using an orthonormal basis and formula (4.4), this reduces to simple algebra. The covariant derivative of the curvature tensor can also be computed using Leibniz rule, for any $u, v, w, z \in \mathfrak{g}$:

$$ (\nabla_z R)(u, v)w = \nabla_z(R(u, v)w) - R(\nabla_z u, v)w - R(u, \nabla_z v)w - R(u, v)\nabla_z w \quad (4.11) $$

Of course this is also valid for right-invariant metrics. This is implemented in `geomstats` in the class `InvariantMetric` as follows.

```c
def curvature_at_identity(
    self, tangent_vec_a, tangent_vec_b, tangent_vec_c):
    """Compute the curvature at identity."""
    bracket = Matrices.bracket(tangent_vec_a, tangent_vec_b)
    bracket_term = self connection_at_identity(bracket, tangent_vec_c)
    left_term = self connection_at_identity(
        tangent_vec_a, self connection_at_identity(tangent_vec_b, tangent_vec_c))
    right_term = self connection_at_identity(
        tangent_vec_b, self connection_at_identity(tangent_vec_a, tangent_vec_c))
    return left_term - right_term - bracket_term

def sectional_curvature_at_identity(self, tangent_vec_a, tangent_vec_b):
    """Compute the sectional curvature at identity."""
    curvature = self.curvature_at_identity(tangent_vec_a, tangent_vec_b, tangent_vec_b)
    num = self inner product(curvature, tangent_vec_a)
    denom = (
        self.squared_norm(tangent_vec_a)
        * self.squared_norm(tangent_vec_a)
        - self inner product(tangent_vec_a, tangent_vec_b) ** 2
    )
    condition = gs.isclose(denom, 0.)
    denom = gs where(condition, 1., denom)
    return gs where(~condition, num / denom, 0.)

def curvature_derivative_at_identity(
    self, tangent_vec_a, tangent_vec_b, tangent_vec_c, tangent_vec_d):
    """Compute the covariant derivative of the curvature at identity."""
    first_term = self connection_at_identity(
        tangent_vec_a,
        self curvature_at_identity(
            tangent_vec_b, tangent_vec_c, tangent_vec_d))
    second_term = self curvature_at_identity(
        self connection_at_identity(tangent_vec_a, tangent_vec_b),
        tangent_vec_c,
        tangent_vec_d)
    third_term = self curvature_at_identity(
        tangent_vec_b,
        self connection_at_identity(tangent_vec_a, tangent_vec_c),
        tangent_vec_d)
```
fourth_term = self.curvature_at_identity(
    tangent_vec_b,
    tangent_vec_c,
    self connection_at_identity(tangent_vec_a, tangent_vec_d))

return first_term - second_term - third_term - fourth_term

If the metric is bi-invariant, it is well known that these formulas greatly simplify (La-
fontaine et al., 2004, Proposition 3 .17) or (Gallier and Quaintance, 2020, Proposition
20.19).

**Proposition 4.3.1.** If $G$ is a Lie group equipped with a bi-invariant metric, and if $X, Y, Z, T \in \mathcal{L}(G)$, then

$$
\langle R(X, Y)Z, T \rangle = \frac{1}{4} \langle [X, Y], [Z, T] \rangle \quad (4.12)
$$

In particular, $G$ has non-negative sectional curvature.

We now turn to the case of submersions defined by the action of a Lie group on a
manifold, and the properties of the induced quotient metric.

4.4 Group action and homogeneous spaces

To conclude this section, we define the notion of group action, that will be fundamental in
Section 5. It allows to model groups as generators of transformations of any manifold. In
Example 2.5 (page 7), the action of SO(m) on M<sub>m,k</sub>(R) is in fact the main factor to study
the geometry of the quotient space SΣ<sub>m</sub><sup>k</sup>.

**Definition 4.4.1 (Group action).** Given a set $M$ and a group $G$, a *left-action* of $G$ on $M$ is a function $\triangleright: G \times M \to M$, such that:

$$
\bullet \quad \forall g, h \in G, \forall x \in M, \quad g \triangleright (h \triangleright x) = (gh) \triangleright x,
$$

$$
\bullet \quad \forall x \in M, \quad e \triangleright x = x.
$$

If furthermore *M* is a smooth manifold and *G* is a Lie group, we say that this action is smooth if the map ▷ is smooth from the product manifold *G* × *M* to *M*. In this case, for any *g* ∈ *G*, *x* ↦ *g* ▷ *x* is a diffeomorphism of *M*, whose inverse is *x* ↦ *g*⁻¹ ▷ *x*. We call this map *left translation* by *g* and write it *L*<sub>*g*</sub> by analogy with the group law.

**Remark 4.4.1.** Note that the map $g \mapsto L_g$ is a group homomorphism between $G$ and $\mathrm{Diff}(M)$. In fact, if this mapping is injective, we say that the action is *faithful* and we can see $G$ as an immersed subgroup of $\mathrm{Diff}(M)$.

We say that an action is
* free if for all $g \in G$ and $x \in M$, if $g \triangleright x = x$, then $g = e$;

* faithful if for all $g \in G$, if $\forall x \in M, g \triangleright x = x$ then $g = e$;

* transitive if for all $x, y \in M$, there exists $g \in G$ such that $y = g \triangleright x$;

* proper if for all compact set $K \subset M$, the set $\{g \in G \mid K \cap gK \neq \emptyset\}$ is compact. This is always the case if $G$ is compact.

Further, The define for any $x \in M$ the *orbit* of $x$ as $[x] = G \triangleright x = \{g \triangleright x \mid g \in G\}$. As in Example 2.5 (page 7), the orbits are the equivalence classes of the following relation

$$x \sim y \iff \exists g \in G, y = g \triangleright x \iff y \in g \triangleright x.$$

Note that a free action means that $G_x = e$ for all $x \in M$. In contrast, if an action is faithful, then the map $x \mapsto g \triangleright x$ is different from the identity of $M$ for all $g \in G$. These notions play an essential role in determining the geometry of $M$. For example, we shall see in Kendall shape spaces that the action of $SO(n)$ is free if we remove certain points, that are otherwise considered as singularities.

**Definition 4.4.2 (Quotient space).** The set of orbits is denoted $M/G$ and is called *quotient* of $M$ by $G$. This set is obtained by identifying all the points in an orbit. Define the *canonical projection*

$$\pi : \begin{cases} M & \to M/G \\ x & \mapsto [x] \end{cases} \qquad (4.13)$$

We have the following sufficient conditions on the action to ensure that a quotient space is indeed a smooth manifold (Lee, 2003, Theorem 7.10).

**Theorem 4.4.1 (Quotient manifold).** Let $M$ be a smooth manifold, $G$ a Lie group with a smooth left-action on $M$ that is free and proper. Then there exist a unique differential structure on $M/G$ such that $M/G$ is a smooth manifold and the canonical projection (4.13) is a submersion.

Note that $M/G$ is equipped with the quotient topology, and the differential structure on $M/G$ can be derived by taking charts of $M$ that are *adapted* to the group action in the sense that the open sets intersect with only one orbit. See Lee (2003) for more details. Free and proper actions occur in many manifolds considered until here, e.g. the Stiefel manifold (Example 2.4 page 7).

We now focus on a particular case of the action of a subgroup of $G$ on $G$. At any $x \in M$, we can define the *isotropy* subgroup, or *stabilizer* of $x$ as

$$G_x = \{g \in G \mid g \triangleright x = x\}.$$
Then $G_x$ is a Lie subgroup of $G$. We consider the action of $G_x$ on $G$, or more generally of a Lie subgroup $H$ of $G$, obviously defined by the group law. We have (Gallier and Quaintance, 2020, Corollary 22.10):

**Theorem 4.4.2.** The action of a Lie subgroup $H$ of a Lie group $G$ on $G$ is free and proper.

Theorem 4.4.1 thus applies to the right action $(g, h) \in G \times H \mapsto gh \in G$ and $G/H$ is a smooth manifold such that $\pi: G \to G/H$ is a submersion. The orbits of this right action are the left cosets $\{gH | g \in G\}$ Note that $G$ acts on $G/H$ by $g_1 \triangleright g_2 H = g_1 g_2 H$, and it is clear that this action is transitive. In fact all transitive left actions yield to a quotient space of orbits from a right action.

**Definition 4.4.3 (Homogeneous space).** We say that $M$ is homogeneous if there exists a smooth transitive action of a Lie group $G$ on $M$.

**Remark 4.4.2.** Note that there might be several such Lie groups with a transitive action on a given homogeneous space $M$.

All homogeneous spaces correspond (up to a diffeomorphism) to quotient spaces of the form $G/H$, as stated in the following (Gallier and Quaintance, 2020, Theorem 22.13).

**Theorem 4.4.3.** Let $G$ be a connected Lie group acting smoothly and transitively on a smooth manifold $M$, so that $M$ is homogeneous. Then for any $x \in M$, writing $H = G_x$, the map

$$ \left\{ \begin{array}{lcl} G/H & \to & M \\ gH & \mapsto & g \triangleright x \end{array} \right. $$

is a diffeomorphism, so that $M \simeq G/H$.

Note that the choice of reference point $x$ above does not matter, as for two different choices $x$ and $y$, the action is transitive so there exists $g \in G$ such that $y = g \triangleright x$, and the stabilizers $G_x$ and $G_y = G_{g \triangleright x} = gG_x g^{-1}$ are conjugate, hence isomorphic.

We exemplify all these notions with the sphere, the Stiefel manifold and the manifold of symmetric positive definite (SPD) matrices, although other manifolds in geomstats such as hyperbolic spaces and the Grassmannian are also homogeneous.

**Example 4.16: Hypersphere as a homogeneous space**

The hypersphere embedded in $\mathbb{R}^{d+1}$ can be seen as a homogeneous space by considering the action of the special orthogonal group $SO(d+1)$ of the embedding space. For $R \in SO(d+1), x \in S^d$:

$$ R \triangleright x = Rx \in S^d \qquad (4.14) $$
Indeed, the sphere is stable by the action of rotation matrices, and this action is transitive. Consider a pole $x_0 = (1, 0, \dots, 0) \in S^d$. Its stabilizer is the set of rotations whose axis is $x_0$, i.e.

$$H = \left\{ \begin{pmatrix} 1 & 0 \\ 0 & R \end{pmatrix} \mid R \in SO(d) \right\} \simeq SO(d).$$

Then by Theorem 4.4.3, $S^d = SO(d+1)/SO(d)$. Intuitively this corresponds to identifying all the rotations that share the same axis where rotation axes are described by a unit vector, i.e., a point on the sphere.

### Example 4.17: Stiefel manifold as a homogenous space

Recall that the Stiefel manifold $St(k, n)$ is the set of orthonormal $k$-frames of $\mathbb{R}^n$. The special orthogonal group $SO(n)$ thus acts on $St(k, n)$ by

$$R \triangleright (u_1, \dots, u_k) = (Ru_1, \dots, Ru_k)$$

It is clear that this action is transitive. Furthermore, as a frame can be represented by an orthogonal matrix of size $n \times k$, the above action corresponds to matrix multiplication and is thus smooth.

For any $U \in St(k, n)$, the stabilizer of $U$ is

$$H = \left\{ \begin{pmatrix} I_k & 0 \\ 0 & R \end{pmatrix} \mid R \in SO(n-k) \right\} \simeq SO(n-k).$$

Then by Theorem 4.4.3, $St(k, n) = SO(n)/SO(n-k)$. Intuitively, this corresponds to the fact that a $k$-frame can be completed into an orthonormal basis of $\mathbb{R}^n$, i.e. a matrix in $SO(n)$, and matrices in $H$ map one such completion to all the other completions. Thus quotienting $SO(n)$ by $SO(n-k)$ amounts to identifying the orthonormal bases that agree on the first $k$ vectors.

### Example 4.18: Grassmannian as a homogenous space

Define the Grassmann manifold $Gr(n, k)$ as the set of subspaces of $\mathbb{R}^n$ of dimension $k$. On any $k$-dimensional subspace $U$ is defined a unique orthogonal projector, i.e. a linear map $p$ defined in $\mathbb{R}^n$ for which $p \circ p = p$ and $\mathrm{Im}(p) = U$ and $\mathrm{ker}(p) = U^\perp$. Any such projector is represented by a symmetric matrix $P$ of size $n$, rank $k$ and such
that $P^2 = P$. We thus adopt the representation of $Gr(n, k)$ as the set

$$Gr(n, k) = \{P \in \mathrm{Sym}(n) \mid P^2 = P \text{ and } \mathrm{rank}(P) = k\}.$$

Intuitively, any $k$ dimensional subspace can be represented by an equivalence class of orthonormal bases, so that the Grassmannian is a quotient space of the orthogonal group $O(n)$. To consider a connected group, let $G = SO(n)$, and let the action of $G$ on $Gr(k, n)$ correspond to a change of basis: $Q \triangleright P = QPQ^\top$. It is clear that any rank-$k$ projector can be represented by a matrix of the form

$$P_k = \begin{pmatrix} I_k & 0 \\ 0 & 0 \end{pmatrix}.$$

This means exactly that the action is transitive, and $Gr(k, n)$ is the orbit of $P_k$. Let $H$ be its stabilizer, i.e. the set

$$H = \left\{ \begin{pmatrix} Q & 0 \\ 0 & R \end{pmatrix} \mid Q \in SO(k), R \in SO(n-k) \right\} \\ \simeq SO(k) \times SO(n-k).$$

$Gr(k, n) = SO(n)/(SO(k) \times SO(n-k))$ is therefore a homogeneous space with canonical projection $\pi: Q \in SO(n) \mapsto QP_kQ^\top$ and dimension $k(n-k)$. Of course this manifold also be described as a quotient of $St(k, n)$ by $SO(p)$, where the action is by right multiplication and the projection is identified with $U \in St(k, n) \mapsto UU^\top \in Gr(k, n)$. This quotient can be more practical than that of $SO(n)$ when $n \gg k$. The Grassmann manifold is widely used in applications both in numerical problems such as low-rank matrix decomposition or optimization, and in higher-level applications in machine learning, computer vision and image processing. We refer to Bendokat et al. (2020) for a very complete exposition of this manifold.

### Example 4.19: SPD matrices as a homogenous space

Recall that a symmetric matrix $\Sigma$ is positive definite if it is invertible and $\forall x \in \mathbb{R}^n, x^\top\Sigma x \ge 0$. The set $SPD(n)$ of such matrices is an open set of the vector space $\mathrm{Sym}(n)$ of symmetric matrices, hence it is an embedded manifold and the canonical immersion id: $SPD(n) \to \mathrm{Sym}(n)$ defines a global chart. Furthermore, for any $\Sigma \in SPD(n)$ the tangent space $T_\Sigma SPD(n)$ is identified with $\mathrm{Sym}(n)$.

Now, define the action of $GL(n)$ on $SPD(n)$ by:

$$\triangleright: (A, \Sigma) \in GL(n) \times SPD(n) \mapsto A\Sigma A^\top.$$
This action is sometimes called action by *congruence*. It is smooth and transitive. Indeed, let $\Sigma = PDP^\top, \Sigma' = Q\Delta Q^\top \in SPD(n)$ be two SPD matrices with their eigenvalue decomposition given by the spectral theorem, then $A = Q\Delta^{-1/2}D^{-1/2}P^\top$ is such that $A \triangleright \Sigma = \Sigma'$. The manifold of SPD matrices is thus a homogeneous space.

Finally, the isotropy group of the identity matrix $I_n$ is the orthogonal group $O(n) \subset GL(n)$:

$$ 
\begin{aligned}
G_{I_n} &= \{A \in GL(n) \mid A \triangleright I_n = I_n\} \\
&= \{A \in GL(n) \mid AA^\top = I_n\} = O(n).
\end{aligned}
$$

It is indeed a closed subgroup of $GL(n)$, so its right action on $GL(n)$ is free and proper (Theorem 4.4.2), and by Theorem 4.4.3, $SPD(n)$ is isomorphic to the orbit space $GL(n)/O(n)$. The canonical projection $\pi : \begin{cases} GL(n) & \to SPD(n) \\ A & \mapsto AA^\top \end{cases}$ is a submersion.

In `geomstats`, we implement a class for symmetric matrices and a class for the manifold of SPD matrices that inherits from `OpenSet`.

```python
class SPDMatrices(OpenSet):
    """Class for the manifold of symmetric positive definite matrices."""

    def __init__(self, n):
        super().__init__(
            dim=int(n * (n + 1) / 2),
            ambient_space=SymmetricMatrices(n))
        self.n = n

    def belongs(self, mat, atol=gs.atol):
        """Check if a matrix is symmetric with positive eigenvalues."""
        is_symmetric = self.ambient_space.belongs mat, atol)
        eigvalues = gs.linalg.eigh mat, eigvals_only=True)
        is_positve = gs.all(eigvalues > 0, axis=-1)
        belongs = gs(logical_and(is_symmetric, is_positve))
        return belongs
```

As the covariance matrix of a multivariate random variable is a positive semi-definite matrix, this manifold is ubiquitous in applications such as signal processing, neuroscience, etc. (see Pennec et al., 2020, Chapter 3 and references therein), although some issues arise when the matrices are degenerate.

# 5 Metrics defined by invariance properties

As seen in the previous section, group actions play a key role in determining the geometric properties of a manifold. Designing operations that are either invariant or equivariant to a group action has been a fruitful source of improvement to deal with data lying on
manifolds (e.g. Pennec et al., 2006), and more recently to generalize convolutional neural
networks (Cohen et al., 2019, and references therein).

In this section, we focus on the implementation of Riemannian metrics that verify
invariance properties on quotient manifolds, defining quotient metrics.

In the second case (Section 5.1), we propose a generic framework to implement quotient
metrics, that are defined by a free and proper group action, and a metric that is invariant
to this group action. This framework is exemplified in the case of Kendall shape spaces,
and results in an efficient implementation that allows to compute parallel transport and
curvature. This is, up to our knowledge, the first open-source Python implementation of
such spaces. It was presented at the GSI 2021 conference in Guigui et al. (2021). Likewise,
we implement the quotient metric on correlation matrices studied by Thanwerdas and
Pennec (2021).

This section thus gathers two contributions to the formulation and implementation of
the parallel transport in the cases of invariant metrics and quotient metrics. Our original
implementations allow to apply these results to numerous spaces. For the completeness of
our exposure, we briefly present the particular cases of the invariant metrics on homogeneous
and symmetric spaces in Sections 5.2 and 5.3, and discuss the impact of these properties on
the implementation.

5.1 Submersions and quotient metrics

A particularly interesting setting is when there exists a submersion $\pi : E \to M$ between two Riemannian manifolds $E$ and $M$, such that $\pi$ is compatible with the two metrics. In that case, most of the geometry of $M$ can be deduced from the geometry of $E$ and the knowledge of $\pi$. Such submersions arise for homogeneous spaces, and for quotients by a Lie group, hence we shall revisit some examples from the previous sections with this viewpoint. We first introduce the general notions, then focus on quotient metrics. Metrics on homogeneous spaces will be treated in Section 5.2.

Although the mathematical notions exposed in this section are well known since the
works of O’Neill (1966) and Le and Kendall (1993), this work is original by the common
implementation of these structures in geomstats. This implementation allowed computing
parallel transport on Kendall shape spaces (Section 5.1.3), outperforming the state-of-the-art
implementations, and are available on other quotient spaces with limited additional software
development. The generic implementation is presented in Section 5.1.4, and constitutes one
of the main contributions of this section.

To exemplify this section, we extensively use the work done in collaboration with Elodie
Maignant on Kendall shape spaces, and presented at the GSI 2021 conference in Guigui
et al. (2021). We also use the Bures-Wasserstein metric, with results from Bhatia et al.
(2019) and following Thanwerdas and Pennec (2022).

### 5.1.1 Riemannian submersions

Throughout this section, let $(E, g)$ and $(M, h)$ be two Riemannian manifolds, and $\pi: E \to M$ be a submersion. Recall that this means that at every $p \in E$, $d\pi_p$ is surjective. We adopt the vocabulary of principal fiber bundles, referring to $E$ as the *total space*, and to $M$ as the *base manifold*. In particular, for every $x \in M$, $\pi^{-1}(x) \subset E$ is a submanifold in $E$ usually called *fiber* above $x$, and the tangent space at any $p \in \pi^{-1}(x)$ is $T_p\pi^{-1}(x) = \ker d\pi_p$.

**Definition 5.1.1 (Section).** A section $\sigma$ of $\pi: E \to M$ is a smooth map $\sigma: U \to E$ defined on an open set $U \subseteq M$ of $M$, such that $\forall x \in M$, $\pi(\sigma(x)) = x$. If $U = M$, we say that $\sigma$ is a *global section*.

A section allows to choose a point in a fiber above any $x$ in a small open set. In the sequel, we assume that $\pi$ always admits local sections. This is true for Riemannian submersions in particular by the inverse function theorem.

### 5.1.2 Quotient metric

**Definition 5.1.2 (Horizontal - Vertical subspaces).** Let $x \in M$ and $p \in \pi^{-1}(x)$. The *vertical subspace* of $T_pE$ is defined as $V_p = \ker d\pi_p$ and is the tangent space to the fiber through $p$. The *horizontal subspace* is its orthogonal complement in $T_pE$: $H_p = (V_p)^\perp$.

Any tangent vector $u \in T_pE$ can thus be decomposed into a vertical and a horizontal component, and we write ver and hor respectively the orthogonal projections on the vertical and horizontal subspace. Note that this decomposition depends on the metric on $E$ by the definition of the orthogonal complement for the horizontal subspace. Furthermore, we say that $u$ is horizontal if $u \in H_p$, and a curve $c: I \to E$ is called horizontal if at every time $t \in I$, the velocity of $c$ is horizontal: $c'(t) \in H_{c(t)}$.

**Example 5.1: Bures-Wasserstein metric**

Recall from Example 4.19 (page 70) that the map $A \in GL(n) \mapsto AA^\top \in SPD(n)$ is a smooth submersion. Its differential at any $A$ is for any $H \in M_n(\mathbb{R})$:

$$d\pi_A H = AH^\top + HA^\top,$$

and its kernel is $\ker d\pi_A = \{H | HA^\top \in \text{Skew}(n)\} = \text{Skew}(n)A^{-\top}$. This defines the vertical subspace $V_A$.

Now, consider $GL(n)$ endowed with the restriction of the Frobenius metric. The
orthogonal complement to $V_A$ is

$$
\begin{align*}
H_A &= \{H \in M_n(\mathbb{R}) \mid \forall B \in \text{Skew}(n), \langle H, BA^{-\top} \rangle = 0\} \\
&= \{H \in M_n(\mathbb{R}) \mid \forall B \in \text{Skew}(n), \text{tr}(BA^{-\top}H^{\top}) = 0\} \\
&= \{H \in M_n(\mathbb{R}) \mid HA^{-1} \in \text{Sym}(n)\} = \text{Sym}(n)A.
\end{align*}
$$

Note that in Section 5.3, we will consider the same manifolds and submersion, but instead endow GL(n) with a left-invariant metric. This will result in another geometry on the manifold of SPD matrices.

Recall that as $\pi$ is a submersion, $d\pi_p$ is a linear isomorphism between $H_p$ and $\text{Im}(\pi) = T_x M$. If this identification is additionally isometric, then geodesics and curvature on $M$ will be deduced from $E$.

**Definition 5.1.3 (Riemannian submersion).** Let $\pi : E \to M$ be a smooth submersion between two Riemannian manifolds. $\pi$ is called a *Riemannian submersion* if for every $p \in E$, the map $d\pi_p$ is an isometry between the horizontal subspace $H_p$ of $T_pE$ and $T_xM$ where $x = \pi(p)$.

If $\pi$ is also surjective, this allows to identify tangent vectors of $M$ to horizontal vectors of $E$.

**Definition 5.1.4 (Horizontal lift).** Let $\pi : E \to M$ be a Riemannian submersion that is surjective onto $M$. Let $X$ be a vector field on $M$, then the unique horizontal lift $\bar{X}$ in $E$ is defined such that for every $x \in M$ and every $p \in \pi^{-1}(x)$,

$$d\pi_p \bar{X}_p = X_x.$$

If we have a local section $\sigma$, we speak about the *corresponding* horizontal lift such that
$\forall (x, v) \in TM, \bar{v} \in H_{\sigma(x)} \subset T_{\sigma(x)}E$.

**Example 5.2: Horizontal lift of SPD matrices**

Following Example 5.1, it is known from linear algebra that the submersion $\pi : A \in GL(n) \mapsto AA^\top \in SPD(n)$ is surjective. To compute the horizontal lift of a tangent vector $X$ at $\Sigma = AA^\top$, we seek a symmetric matrix $S$ such that $d\pi_A(SA) = X$. This is true if and only if $AA^\top S + SAA^\top = X$, i.e. if $S \in \text{Sym}(n)$ solves the Sylvester equation

$$\Sigma S + S \Sigma = X. \qquad (5.1)$$

More precisely, using the eigenvalue decomposition of $\Sigma = PDP^\top$, $S$ solves Equa-
tion (5.1) iff

$$ \begin{align} \Sigma S + S\Sigma &= X \notag \\ \iff \quad &DP^\top SP + P^\top SPD = P^\top XP \notag \\ \iff \quad & (P^\top SP)_{ij} = \frac{(P^\top XP)_{ij}}{d_i + d_j}. \tag{5.2} \end{align} $$

Define the map $S_Σ : \text{Sym}(n) \to \text{Sym}(n)$ that solves Sylvester equation using Equation (5.2). It uniquely defines the horizontal lift $\bar{X}$ of $X$ at $A$ by $\bar{X} = S_{AA^\top}(X)A$.

**Proposition 5.1.1.** Similarly, if $c : [0, 1] \to M$ is a piecewise smooth curve in $M$, and $p \in \pi^{-1}(c(0))$, then there exists a unique curve $\bar{c} : [0, 1] \to E$ such that $\pi \circ \bar{c} = c$ and $\bar{c}' \in H$. The curve $\bar{c}$ is called horizontal lift of $c$.

**Geodesics** This leads to the main theorem of this section, that is usually attributed to O'Neill, and relates geodesics of $M$ with horizontal geodesics of $E$ (Gallier and Quaintance, 2020, Proposition 17.8).

**Theorem 5.1.1.** Let $\pi : E \to M$ be a Riemannian submersion between two Riemannian manifolds $E$ and $M$.

(1) If $\gamma$ is a geodesic in $E$ such that $c'(0)$ is horizontal, then $c$ is horizontal, and its projection $\gamma = \pi \circ c$ is a geodesic in $M$ of same length as $c$.

(2) For every $p \in E$, if $\gamma$ is a geodesic in $M$ such that $\gamma(0) = \pi(p)$, then there exists $\epsilon > 0$ such that there exists a unique horizontal lift $c$ of the restriction of $\gamma$ to $[-\epsilon, \epsilon]$ and $c$ is a geodesic of $E$ through $p$.

(3) If furthermore $\pi$ is surjective, for any vector fields $X, Y \in \Gamma(TM)$, $\langle \bar{X}, \bar{Y} \rangle = \langle X, Y \rangle \circ \pi$ (where $\bar{\cdot}$ is the horizontal lift).

The first property of this theorem can be written by the commutative diagram below:

$$ \begin{tikzcd} H_p \arrow[r, "\text{Exp}_p"] \arrow[d, "d\pi_p"'] & E \arrow[d, "\pi"'] \\ T_x M \arrow[r, "\text{Exp}_x"] & M \end{tikzcd} $$

where $p \in E$ and $x = \pi(p)$. Furthermore, as the restriction of $d\pi$ on horizontal spaces is an isomorphism, choosing a section $\sigma$ of $E$ and the corresponding horizontal lift $\bar{\cdot}$ allows to compute the exponential map of $M$ from the one of $E$:

$$ \forall x \in M, \forall v \in T_x M, \quad \text{Exp}_x(v) = \pi \circ \text{Exp}_{\sigma(x)}(\bar{v}) \tag{5.3} $$
**Remark 5.1.1.**

* One can also show that the connection $\nabla_{\bar{X}}\bar{Y}$ on $E$ verifies (see e.g. Lafontaine et al. (2004, Proposition 3.55))

$$ \nabla_{\bar{X}}\bar{Y} = \overline{\nabla_X Y} + \frac{1}{2} \mathrm{ver}[\bar{X}, \bar{Y}]. $$

Consequently, one cannot hope to obtain a similar commutation rule between the parallel transport in $E$ and the one in $M$. Indeed, if $\bar{Y}(t)$ is the horizontal lift of a parallel vector field $Y(t)$ along a curve $\gamma$ (whose horizontal lift is $\bar{\gamma}$), then $\nabla_{\bar{\gamma}'}\bar{Y}$ has a non-zero vertical component given above, which vanishes in the case of a geodesic. Computing parallel transport is thus not as straightforward as computing geodesics in this case.

* A Riemannian submersion shortens distances, i.e. for $p, q \in E$, and writing $d_E$ and $d_M$ the Riemannian distances on respectively $E$ and $M$, we have

$$ d_M(\pi(p), \pi(q)) \le d_E(p, q). $$

**Example 5.3: Bures-Wasserstein geodesics**

There is a unique metric that turns $\pi : A \in GL(n) \mapsto AA^\top \in SPD(n)$ into a Riemannian submersion between $GL(n)$ endowed with the Frobenius metric and $SPD(n)$. This metric is called the Bures-Wasserstein (BW) metric. Let us apply Theorem 5.1.1 to compute the geodesics of this metric. Recall that $(GL(n), \mathrm{Frobenius})$ is flat, so the geodesics are “straight” lines: $\mathrm{Exp}_A(tX) = A + tX$. Let $\Sigma = AA^\top \in SPD(n)$ and $X \in T_\Sigma SPD(n)$, then $\forall t > 0$

$$ \begin{aligned} \mathrm{Exp}_\Sigma(tX) &= \pi \circ \mathrm{Exp}_A(t\bar{X}) = \pi(A + tS_\Sigma(X)A) \\ &= (A + tS_\Sigma(X)A)(A + tS_\Sigma(X)A)^\top \\ &= AA^\top + t(S_\Sigma(X)\Sigma + \Sigma S_\Sigma(X)) + t^2S_\Sigma(X)\Sigma S_\Sigma(X) \\ &= AA^\top + tX + t^2S_\Sigma(X)\Sigma S_\Sigma(X) \end{aligned} $$

Note that a geodesic hits the boundaries of the manifold when $-1/t$ is in an eigenvalue of $S_\Sigma(X)$. $(SPD(n), \mathrm{BW})$ is therefore not complete. All the properties of this metric are derived using the equations of Riemannian submersions in Thanwerdas and Pennec (2022) and Thanwerdas (2022).

The relation between the connection of the total space and that of the base manifold also translates into a relation between the curvatures of each space.
**Curvature** O'Neill (1966) showed that the curvatures of the total space and of the base space were related, and that those of the base space could be computed using two fundamental tensors defined by the horizontal and vertical projections of the connection. We state here the main result, that shows that sectional curvature can only increase after a Riemannian submersion.

**Theorem 5.1.2 (O'Neill).** Let $\pi : E \to M$ be a Riemannian submersion, and $X, Y$ be orthonormal vector fields on $M$, with horizontal lifts $\bar{X}, \bar{Y}$. Then

$$ \kappa(X, Y) = \kappa(\bar{X}, \bar{Y}) + \frac{3}{4} \| \text{ver}[\bar{X}, \bar{Y}] \|^2. \qquad (5.4) $$

**Example 5.4: Curvature of the Bures-Wasserstein metric**

As $GL(n)$ with the Frobenius metric is flat, O'Neill's theorem implies that the space of SPD matrices endowed with the Bures-Wasserstein metric has non-negative curvature.

We now apply these results to the special case of submersions defined as the canonical projection to a quotient space.

**Metric** Recall from Theorem 4.4.1 (page 67) that if $G$ is a Lie group and $E$ a smooth manifold such that $G$ acts on $E$ and the action is smooth, free and proper, then the canonical projection $\pi : x \in E \mapsto [x] \in M = E/G$ is a submersion.

In this case, the action of $G$ allows to move in the fibers. More precisely, as the fibers are defined as orbits $G \triangleright p = \{g \triangleright p | g \in G\}$ for $p \in E$, the fibers are stable by the action of $G$, i.e. $\forall x \in M, \forall g \in G, \forall p \in \pi^{-1}(x), g \triangleright p \in \pi^{-1}(x)$, which can also be written $\pi(g \triangleright p) = \pi(p) = x$. Moreover, this action is transitive on fibers, i.e. for any $q \in E$ such that $\pi(q) = x$, there exists $g \in G$ such that $q = g \triangleright p$.

**Definition 5.1.5 (Invariant metric).** A Riemannian metric on $E$ is $G$-invariant if for any $g \in G$, $L_g : p \in E \mapsto g \triangleright p \in E$ is an isometry.

Now suppose that $E$ is equipped with a $G$-invariant metric. Then one can define a metric on the quotient manifold $M = E/G$ such that $\pi$ is a Riemannian submersion (Lafontaine et al., 2004, Proposition 2.28).

**Proposition 5.1.2 (Quotient metric).** Let $E$ be a smooth manifold, $G$ a Lie group acting smoothly, properly and freely on $G$, and $\langle \cdot, \cdot \rangle$ be a $G$-invariant metric on $E$. Let $\pi : E \to M$ be the canonical projection. Then there exists a unique Riemannian metric on $M$ such that $\pi$ is a Riemannian submersion. Let $x \in M$ and $p \in \pi^{-1}(x)$, for any $u, v \in T_x M$, it is given by

$$ \langle u, v \rangle_x \triangleq \langle \bar{u}, \bar{v} \rangle_p. \qquad (5.5) $$
In both cases, automatic differentiation is used to compute the gradient of the objective function.

```python
def align(self, point, base_point,
          max_iter=25, verbose=False, tol=gs.atol):
    """Align point to base_point by minimization."""
    ... # definitions of wrap for different cases
    objective_withGrad = gs.autograd.value_andGrad(
        lambda param: self.ambientMetric.squared_dist(
            wrap(param), base_point))
    init_param = gs flatten(gs.randomrand(*maxShape))
    res = minimize(
        objective_withGrad, init_param, method='L-BFGS-B', jac=True,
        options={ 'disp': verbose, 'maxiter': max_iter}, tol=tol)
    return wrap(res.x)
```

This map is particularly useful because it allows to compute the Log map of the quotient space from the one of the total space. Indeed, choosing a section $\sigma$ and the corresponding horizontal lift $\bar{\sigma}$, we can deduce from Theorem 5.1.1 the following relation:

$$ \forall x, y \in M, \mathrm{Log}_x(y) = d\pi_{\sigma(x)} \mathrm{Log}_x (\omega(x, y)). \qquad (5.7) $$

### Example 5.5: Bures-Wasserstein distance

Continuing Example 5.2 where the submersion $A \mapsto AA^\top$ from $GL(n)$ to $SPD(n)$, and explaining why the metric of Example 5.3 is correctly defined, we define the quotient metric of $GL(n)$ with the Frobenius metric by the orthogonal group. Indeed, consider the right-action of $O(n)$ on $GL(n)$ by matrix multiplication. As shown in Example 4.19 (page 70), the canonical projection of $GL(n)/O(n)$ coincides with $\pi: A \mapsto AA^\top$.

The metric hereby defined on the manifold of SPD matrices is called the Bures-Wasserstein metric (Bhatia et al., 2019). It arises naturally as the $L_2$-Wasserstein distance on multivariate centred normal random variables. It also corresponds to the Procrustes reflection size-and-shape metric (Dryden et al., 2009) taken on the square root SPD matrices. As $GL(n)$ is Euclidean under the Frobenius metric, the geodesics and curvature are known and can be projected to the quotient space. The minimization problem (5.6) can be solved in closed-form to obtain (Thanwerdas and Pennec, 2022; Thanwerdas, 2022).

$$ d(\Sigma, \Sigma')^2 = \mathrm{tr}(\Sigma) + \mathrm{tr}(\Sigma') - 2 \mathrm{tr}((\Sigma \Sigma')^{1/2}). $$

In the next subsection, we focus on Kendall shape spaces, as most of the geometry of these spaces can be computed using the results of this section. Our implementation in `geomstats` is, to the best of our knowledge, the first open-source Python implementation of Kendall shape spaces, and allows efficient computation of parallel transport.
### 5.1.3 Application to Kendall shape spaces

The following subsection is inspired by Nava-Yazdani et al. (2020) but constitutes an original contribution by putting all the code open source in our package geomstats, and by improving the implementation of parallel transport. It is a collaboration with Elodie Maignant and was presented at GSI 2021 in Guigui et al. (2021).

We revisit Example 2.5 (page 7) with the additional idea that the size of a set of landmarks may be filtered out to define a shape. The study of these spaces, including their mathematical structure, the properties of statistical distributions and estimation methods for shape data and their applications to many scientific fields goes back to the late 1970’s with the works of Mardia, Bookstein and Kendall among others. For historical notes on this research area and an introduction to the applications Dryden and Mardia (see 2016, Preface and Section 1.4). Recall that we consider the set of centred matrices of size $m \times k$ as the space of configurations (or landmarks). We assume that $m \ge 2$, and refer the reader to Le and Kendall (1993) for more details.

To further remove the effects of scaling, we restrict to non-zero $x$ (i.e. at least two landmarks are different), and divide $x$ by its Frobenius norm (written $\| \cdot \|$). This defines the pre-shape space

$$S_m^k = \{x \in M(m, k) \mid \sum_{i=1}^k x_i = 0, \|x\| = 1\},$$

which is identified with the hypersphere of dimension $m(k-1) - 1$. The pre-shape space is therefore a differential manifold whose tangent space at any $x \in S_m^k$ is given by

$$T_x S_m^k = \{w \in M(m, k) \mid \sum_{i=1}^k w_i = 0, \text{tr}(w^\top x) = 0\}.$$

Moreover, the rotation group $SO(m)$ acts on $S_m^k$ by matrix multiplication, and this action corresponds to applying a rotation to each landmark individually. As $SO(m)$ is compact, this action is proper. However, this action is not free everywhere if $m \ge 3$. This makes the orbit space

$$\Sigma_m^k = \{[x] \mid x \in S_m^k\} = S_m^k/\text{SO}(m)$$

a “differential manifold with singularities where the action is not free”, and these points correspond to matrices of rank $m-2$ or less (i.e. some landmarks are aligned).

By Theorem 4.4.1 (page 67), the canonical projection $\pi: x \mapsto [x]$ is a submersion. For any $x \in S_m^k$ and $A \in \text{Skew}(m)$, as for all $t \in \mathbb{R}$, $[\exp(tA)x] = x$, the curve $t \mapsto \exp(tA)x$ is a curve in the fiber through $x$, so the vertical space at $x$ is

$$V_x = \{Ax \mid A \in \text{Skew}(m)\} = \text{Skew}(m)x.$$
The pullback of Frobenius metric on the pre-shape space allows to define the horizontal spaces:

$$ 
\begin{aligned}
H_x &= \{w \in T_x S_m^k \mid \mathrm{Tr}(Axw^\top) = 0 \quad \forall A \in \mathrm{Skew}(m)\} \\
&= \{w \in T_x S_m^k \mid xw^\top \in \mathrm{Sym}(m)\},
\end{aligned}
$$

where $\mathrm{Sym}(m)$ is the space of symmetric matrices of size $m$. Lemma 1 from Nava-Yazdani et al. (2020) allows to compute the vertical component of any tangent vector.

**Lemma 5.1.3.** For any $x \in S_m^k$ and $w \in T_x S_m^k$, the vertical component of $w$ can be computed as $V_x(w) = Ax$ where $A$ solves the Sylvester equation:

$$ Axx^\top + xx^\top A = wx^\top - xw^\top. \qquad (5.8) $$

If $\mathrm{rank}(x) \ge m - 1$, the skew-symmetric solution $A$ of Equation (5.8) is unique.

In practice, the Sylvester equation can be solved by an eigenvalue decomposition of $xx^\top$ (see Example 5.2). This defines $\mathrm{ver}_x$, the orthogonal projection on $V_x$. As $T_x S_m^k = V_x \oplus H_x$, any tangent vector $w$ at $x \in S_m^k$ may be decomposed into a horizontal and a vertical component, by solving Equation (5.8) to compute $\mathrm{ver}_x(w)$, and then $\mathrm{hor}_x(w) = w - \mathrm{ver}_x(w)$.

As the Frobenius metric is invariant to the action of $SO(m)$, we can define the quotient metric on $\Sigma_m^k$, and this makes $\pi$ a Riemannian submersion. Furthermore, the Riemannian distances $d$ on $S_m^k$ and $d_\Sigma$ on $\Sigma_m^k$ are related by

$$ d_{\Sigma}(\pi(x), \pi(y)) = \inf_{R \in SO(m)} d(x, Ry). $$

The optimal rotation $R$ between any $x, y$ is unique in a subset $U$ of $S_m^k \times S_m^k$, which allows to define the align map $\omega: U \to S_m^k$ that maps $(x, y)$ to $Ry$. In this case, $d_\Sigma(\pi(x), \pi(y)) = d(x, \omega(x, y))$ and $x\omega(x, y)^\top \in \mathrm{Sym}(m)$. It is useful to notice that $\omega(x, y)$ can be directly computed by a pseudo-singular value decomposition of $xy^\top$ (Kendall and Le, 2009).

**Remark 5.1.2.** The alignment problem is similar to the canonical correlation analysis problem (CCA) between two data sets:

$$ \max_{U,V \in O(m)} \langle Ux, Vy \rangle. $$

Finally, in the case of the Kendall shape spaces, the quotient space cannot be seen explicitly as a submanifold of some $\mathbb{R}^N$. Moreover, the projection $\pi$ and its derivative $d\pi$ cannot be computed. However, the align map amounts to identifying the shape space with a local horizontal section of the pre-shape space, and thanks to the characteristics of Riemannian submersions mentioned in the previous subsections, all the computations can be done in the pre-shape space.
Let Exp, Log, and $d$ denote the operations of the pre-shape space $S_m^k$, that are given in Equation (3.16), (3.17) (page 35). We obtain from Theorem 5.1.1 for any $x, y \in S_m^k$ and $v \in T_x S_m^k$

$$ 
\begin{aligned}
\mathrm{Exp}_{\Sigma, [x]}(d\pi_x v) &= \pi(\mathrm{Exp}_x(\mathrm{hor}_x(v))), \\
\mathrm{Log}_{\Sigma, [x]}([y]) &= d\pi_x \mathrm{Log}_x(\omega(x, y)), \\
d\Sigma([x], [y]) &= d(x, \omega(x, y)).
\end{aligned}
$$

To end this section, we state Proposition 2 of Kim et al. (2021), that echoes with Remark 5.1.1 for the computation of parallel transport.

**Proposition 5.1.4.** Let $\gamma$ be a horizontal $C^1$-curve in $S_m^k$ and $v$ be a horizontal tangent vector at $\gamma(0)$. Assume that $\mathrm{rank}(\gamma(s)) \ge m-1$ except for finitely many $s$. Then the vector field $s \mapsto v(s)$ along $\gamma$ is horizontal and the projection of $v(s)$ to $T_{[\gamma(s)]}S_m^k$ is the parallel transport of $d\pi_x v$ along $[\gamma(s)]$ if and only if $s \mapsto v(s)$ is the solution of

$$ \dot{v}(s) = -\mathrm{tr}(\dot{\gamma}(s)v(s)^{\top})\gamma(s) + A(s)\gamma(s), \quad v(0) = v \quad (5.9) $$

where for every $s$, $A(s) \in \mathrm{Skew}(m)$ is the unique solution to

$$ A(s)\gamma(s)\gamma(s)^{\top} + \gamma(s)\gamma(s)^{\top}A(s) = \dot{\gamma}(s)v(s)^{\top} - v(s)\dot{\gamma}(s)^{\top}. \quad (5.10) $$

Equation (5.9) means that the covariant derivative of $s \mapsto v(s)$ along $\gamma$ must be a vertical vector at all times, defined by the matrix $A(s) \in \mathrm{Skew}(m)$. These equations can be used to compute parallel transport in the shape space. To compute the parallel transport of $d\pi_x w$ along $[\gamma]$, Kim et al. propose the following method: one first chooses a discretization time-step $\delta = \frac{1}{n}$, then repeat for every $s = \frac{i}{n}, i = 0 \dots n$

1. Compute $\gamma(s)$ and $\dot{\gamma}(s)$,

2. Solve the Sylvester Equation (5.10) to compute $A(s)$ and the r.h.s. of (5.9),

3. Take a discrete Euler step to obtain $\tilde{v}(s + \delta)$

4. Project $\tilde{v}(s + \delta)$ to $T_{\gamma(s)}S_m^k$ to obtain $\hat{v}(s + \delta)$,

5. Project to the horizontal subspace: $v(s + \delta) \leftarrow \mathrm{hor}(\hat{v}(s + \delta))$

6. $s \leftarrow s + \delta$

We notice that this method can be accelerated by a higher-order integration scheme, such as Runge-Kutta (RK) by directly integrating the system $\dot{v} = f(v, s)$ where $f$ is a smooth map given by Equations (5.9) and (5.10). In this case, steps 4. and 5. are not necessary. The precision and complexity of this method is then bound to that of the integration scheme used, as shown in Figure 5.1 for randomly generated orthogonal tangent vectors of unit norm.
<figure><img src="image_18.png" /><figcaption>**Figure 5.1:** Convergence speed of the integration of Equation (5.9) for two orthonormal initial vectors $v$ and $\dot{v}(0)$, with respect to the number of steps. The curve $\gamma$ is a horizontal geodesic here.</figcaption></figure>

### 5.1.4 Implementation

We now present an abstract class to construct a quotient metric. There are two scenarios. In the first one, a submersion $\pi: E \to M$ is given between the total space and the base manifold, together with a metric on $E$ whose restriction to horizontal spaces is preserved by $\pi$. This is the case of the Bures-Wasserstein metric (Examples 5.2, 5.3 and 5.5).
In the second scenario, we are only given

1. a total space $E$,

2. a group $G$ acting freely, properly and smoothly on $E$,

3. a metric on $E$ invariant to the action of $G$.

In this case, the base manifold cannot be represented explicitly, that is, as an embedded submanifold of $\mathbb{R}^N$. Thus, we cannot implement it as an `OpenSet` or `LevelSet` class on its own. However, we can use the properties of the canonical projection, that is a Riemannian submersion, to construct a metric on the quotient space. This metric is in fact defined on horizontal spaces of the total space, and the base manifold is locally identified with a section of the total space. This is the case of Kendall shape spaces (Section 5.1.3).

To model both cases, we use the structure of fiber bundle. The first scenario is a fiber bundle by definition. For the second scenario, we notice that $\pi: E \to M$ and $(G, \cdot)$ form a principal fiber bundle. Conversely, if we endow a principal fiber bundle with a $G$-invariant metric, we can define a quotient metric on $B$. We thus choose to implement an abstract `FiberBundle` class, as main ingredient of a `QuotientMetric` class.

As the goal is to compute the Riemannian Exp, Log and distance maps according to Equations (5.3), (5.7) and (5.6), we need to specify a section and the corresponding horizontal lift, and the align map. Then the quotient metric is implemented as follows

```python
class QuotientMetric(RiemannianMetric):
    """Quotient metric."""

    def __init__(self, fiberBundle: FiberBundle, dim: int = None):
        if dim is None:
            if fiberBundle.base is not None:
                dim = fiberBundle.base.dim
            elif fiberBundle.group is not None:
                dim = fiberBundle.dim - fiberBundle.group.dim
        else:
            raise ValueError('Either the base manifold, '
                             'its dimension, or the group acting on the '
                             'total space must be provided.')

    super().__init__(
        dim=dim, default_point_type=fiberBundle.default_point_type)
    self.fiberBundle = fiberBundle
    self.group = fiberBundle.group
    self.ambientMetric = fiberBundle.ambientMetric

    def inner_product(
        self, tangent_vec_a, tangent_vec_b, base_point=None,
        pointAbove=None):
        """Compute the inner-product of two tangent vectors at a base point."""
        if pointAbove is None:
            if base_point is not None:
                pointAbove = self.fiberBundle.lift(base_point)
        else:
            raise ValueError(
                'Either a point (of the total space) or a base point (of '
                'the quotient manifold) must be given.')
        hor_a = self.fiberBundle horizontal lift(tangent_vec_a, pointAbove)
```
```python
hor_b = self.fiberBundle horizontal_lift (tangent_vec_b, pointAbove)
return self.ambientMetric. inner_product (hor_a, hor_b, pointAbove)
```

```python
def exp(self, tangent_vec, base_point, **kwargs):
    """Compute the Riemannian exponential of a tangent vector."""
    lift = self.fiberBundle. lift (base_point)
    horizontal_vec = self.fiberBundle horizontal_lift (tangent_vec, lift)
    return self.fiberBundle.riemannian_submersion(
        self.ambientMetric. exp (horizontal_vec, lift))
```

```python
def log(self, point, base_point, **kwargs):
    """Compute the Riemannian logarithm of a point."""
    point_fiber = self.fiberBundle. lift (point)
    bp_fiber = self.fiberBundle. lift (base_point)
    aligned = self.fiberBundle. align (point_fiber, bp_fiber, **kwargs)
    return self.fiberBundle. tangent_riemannian_submersion(
        self.ambientMetric. log (aligned, bp_fiber), bp_fiber)
```

```python
def squared_dist(self, point_a, point_b, **kwargs):
    """Squared geodesic distance between two points."""
    lift_a = self.fiberBundle. lift (point_a)
    lift_b = self.fiberBundle. lift (point_b)
    aligned = self.fiberBundle. align (lift_a, lift_b, **kwargs)
    return self.ambientMetric. squared_dist (aligned, lift_b)
```

## Example 5.6: Implementation of the BW metric

In the setting of Examples 5.2, 5.3 and 5.5, that corresponds to the first scenario as $E = GL(n)$ and $M = SPD(n)$, we create a new class that corresponds to $GL(n)$ with the fiber bundle structure: Note that all computations can actually be carried out in closed form (see Bhatia et al., 2019), and this specific examples allows to test our general implementation. In this case we used the minimization procedure to compute the align map.
```python
class BuresWassersteinBundle(GeneralLinear, FiberBundle):
    def __init__(self, n):
        super().__init__(
            n=n, base=SPDMatrices(n), group=SpecialOrthogonal(n),
            ambientMetric=MatricesMetric(n, n))

    @staticmethod
    def riemannian_submersion(point):
        return Matrices.mul(point, Matrices transpose(point))

    @staticmethod
    def lift(point):
        return gs.lINALg.cholesky(point)

    def tangent_riemannian_submersion(self, tangent_vec, base_point):
        product = Matrices.mul(
            base_point, Matrices transpose(tangent_vec))
        return 2 * Matrices.to_symmetric(product)

    def horizontal_lift(self, tangent_vec, pointAbove=None, base_point=None):
        if base_point is None and pointAbove is not None:
            if pointAbove is not None:
                base_point = self.riemannian_submersion(pointAbove)
            else:
                raise ValueError(
                    'Either a point in the fiber or a base point in base manifold')
                    ' must be given.')
        sylvester = gs.lINALg.solve_sylvester(
            base_point, base_point, tangent_vec)
        return Matrices.mul(sylvester, pointAbove)
```

In the second scenario, we only compute in $E$ so both submersion and lift maps are the identity of $E$. This can be set by default and overridden as in Example 5.6. Furthermore, in that case, tangent vectors to $M$ are identified with horizontal vectors of $E$, so we need to compute the horizontal decomposition of a tangent vector to $E$. Thus, either the horizontal or vertical projection needs to be implemented, and the other one can be deduced.

In the first scenario, these are given by using $d\pi$ and the horizontal lift, as vertical vectors are in the kernel of $d\pi$ by definition.

This explains the following `FiberBundle` class (the align method is not duplicated from above in the interest of space):

```python
class FiberBundle(Manifold, abc.ABC):
    """Class for (principal) fiber bundles."""

    def __init__(
        self, dim: int, base: Manifold = None,
        group: LieGroup = None, ambientMetric: RiemannianMetric = None,
        group_action=None, **kwargs):
        super().__init__(dim=dim, **kwargs)
        self.base = base
        self.group = group
        self.ambientMetric = ambientMetric

        if group_action is None and group is not None:
            group_action = groupposed
        self.group_action = group_action

    @staticmethod
```
```python
def riemannian_submersion(point):
    """Project a point to base manifold."""
    return point

@staticmethod
def lift(point):
    """Lift a point to total space."""
    return point

def align(self, point, base_point,
          max_iter=25, verbose=False, tol=gs.atol):
    """Align point to base_point by optimization in the Lie algebra."""
    pass

def tangent_riemannian_submersion(self, tangent_vec, base_point):
    """Project a tangent vector to base manifold."""
    return self horizontal projection (tangent_vec, base_point)

def horizontallift(self, tangent_vec, pointAbove=None, base_point=None):
    """Lift a tangent vector to a horizontal vector in the total space."""
    if pointAbove is None:
        if base_point is not None:
            pointAbove = self .lift (base_point)
        else:
            raise ValueError(
                'Either a point (of the total space) or a base point (of'
                'the base manifold) must be given.')
    return self horizontal projection (tangent_vec, pointAbove)

def horizontal projection (self, tangent_vec, base_point):
    """Project to horizontal subspace."""
    try:
        return tangent_vec - self .vertical projection (
            tangent_vec, base_point)
    except (RecursionError, NotImplementedError):
        return self horizontal lift (
            self .tangent_riemannian submersion (
                tangent_vec, base_point), base_point)

def vertical projection (self, tangent_vec, base_point, **kwargs):
    """Project to vertical subspace."""
    try:
        return tangent_vec - self .horizontal projection (
            tangent_vec, base_point)
    except RecursionError:
        raise implementedError
```

Example 5.7: Kendall shape metric

With the above construction of the second scenario, Kendall shape metric is just
a subclass of the QuotientMetric class. We also add the parallel transport method
discussed above and used in Figure 5.1.