class KendallShapeMetric(QuotientMetric):
    """Quotient metric on the shape space."""
    def __init__(self, k_landmarks, m_ambient):
        bundle = PreShapeSpace(k_landmarks, m_ambient)
        super().__init__(
            fiber bundle = bundle,
            dim = bundle.dim - int(m_ambient * (m_ambient - 1) / 2))
    def parallel_transport(
        self, tangent_vec_a, tangent_vec_b, base_point, n_steps=100,
        step='rk4'):
        """Compute the parallel transport of a tangent vec along a geodesic."""
        horizontal_a = self.fiber bundle horizontal projection(
            tangent_vec_a, base_point)
        horizontal_b = self.fiber bundle horizontal projection(
            tangent_vec_b, base_point)
    def force (state, time):
        gamma_t = self. ambient metric. exp (time * horizontal_b, base_point)
        speed = self. ambient metric. parallel_transport(
            horizontal_b, time * horizontal_b, base_point)
        coef = self. inner product (speed, state, gamma_t)
        normal = gs. einsum ('...ij->...ij', coef, gamma_t)
        align = gs. matmul (Matrices. transpose (speed), state)
        right = align - Matrices. transpose (align)
        left = gs. matmul (Matrices. transpose (gamma_t), gamma_t)
        skew = gs. linalg. solve sylvester (left, left, right)
        vertical = - gs. matmul (gamma_t, skew)
        return vertical - normal

    flow = integrate (force, horizontal_a, n_steps=n_steps, step=step)
    return flow[-1]

Example 5.8: Full rank correlation matrices

As a last example, we implement the metric on the set of Correlation matrices described in Thanwerdas and Pennec (2021) using our abstract FiberBundle and QuotientMetric classes.

The set of full-rank correlation matrices Corr(n) is a submanifold of SPD(n)
formed by matrices with unit diagonal. Define the action of positive diagonal matrices
on SPD(n) by congruence:

$$
\triangleright : (D, \Sigma) \in \mathrm{Diag}_{+}(n) \times \mathrm{SPD}(n) \mapsto D\Sigma D \in \mathrm{SPD}(n).
$$

This action is smooth, free, and proper (David and Gu, 2019). The quotient manifold
SPD(n)/Diag+(n) is thus well defined, and can be identified with Corr(n) by the
map [Σ] ↦ DΣ⁻¹/² ▷ Σ where DΣ is the diagonal matrix with coefficients Σᵢᵢ.

Let $\pi: \mathrm{SPD}(n) \to \mathrm{Corr}(n)$ be the canonical projection composed with this map,
i.e. $\Sigma \mapsto D_{\Sigma}^{-1/2} \triangleright \Sigma$, and consider the affine-invariant metric defined on $\mathrm{SPD}(n)$ by

$$
g_{\Sigma}(V, W) = \mathrm{tr}(\Sigma^{-1}V\Sigma^{-1}W).
$$
This metric will be detailed in Example 5.14 in Section 5.3. It is clear that it is invariant by the action of diagonal matrices.

Therefore we can define the quotient metric such that $\pi$ is a Riemannian submersion. The horizontal and vertical spaces are given Thanwerdas and Pennec (in 2021, Theorem 1). The submersion and its differential, the vertical projection and horizontal lifts are computed in closed form, so we are in the first scenario, and the quotient metric is simply

```python
class FullRankCorrelationAffineQuotientMetric(QuotientMetric):
    """Class for the quotient of the affine-invariant metric."""

    def __init__(self, n):
        super().__init__(
            fiber bundle=CorrelationMatricesBundle(n=n))
```

where the `FiberBundle` class is used to define `CorrelationMatricesBundle` as follows:

```python
class CorrelationMatricesBundle(SPDMatrices, FiberBundle):
    """Fiber bundle for the quotient metric on correlation matrices."""

    def __init__(self, n):
        super().__init__(
            n=n, base=CorrelationMatrices(n),
            ambient metric=SPDMetricAffine(n), group dim=n,
            group action=CorrelationMatrices.diag_action)
```
that invariant metrics on $G/H$ are in fact particular cases of the quotient metrics of the previous section. Indeed, a particular case of Theorem 5.1.2 (page 77) gives the following proposition.

**Proposition 5.2.1.** Let $G$ be a Lie group and $H \subseteq G$ a closed subgroup of $G$. Write $o = eH = H \in G/H$, and $\mathfrak{g}, \mathfrak{h}$ the Lie algebras of respectively $G, H$. Let $g$ be a left-invariant Riemannian metric on $G$, that is also right-invariant by $H$. Then there exists a unique Riemannian metric on $G/H$ that is invariant to the action of $G$ and such that $d\pi_e$ is an isometry between $\mathfrak{h}^{\perp} \subset \mathfrak{g}$ and $T_oG/H$. In fact, $\pi$ is a Riemannian submersion.

Indeed, Theorem 5.1.2 applies with $E = G$ and $G = H$ as the action of $H$ on $G$ is smooth, free and proper (by Theorem 4.4.2), and the metric on $G$ is invariant to this action. This defines a metric on $G/H$, and we can show that it is invariant to $G$. Indeed, for any $x, y \in G$, write $M_{yx^{-1}}: G/H \to G/H$ that maps some $[p]$ to $[yx^{-1}p]$, and $L_{yx^{-1}}: G \to G$ the left translation of $G$ by $yx^{-1}$. Then by definition of the action of $G$ on $M$, we have the following commutation

$$ \pi \circ L_{yx^{-1}} = M_{yx^{-1}} \circ \pi $$

Differentiating the above expression at $x$ and as $d\pi_x, d\pi_y$ and $d(L_{yx^{-1}})_x$ are isometries, we conclude that $d(M_{yx^{-1}})_{[x]}$ is also an isometry.

## 5.2.2 Existence

The above proposition thus states that we need a left-invariant metric on $G$ that is also right-invariant by $H$ to define an invariant metric on $G/H$. There are necessary and sufficient conditions on $G$ and $H$ for these to exist (see e.g. Gallier and Quaintance, 2020, Proposition 22.22).

**Proposition 5.2.2.** If $G$ acts faithfully on $G/H$ and if $\mathfrak{g}$ admits a decomposition $\mathfrak{g} = \mathfrak{h} \oplus \mathfrak{m}$ with $\mathrm{AD}_H(\mathfrak{m}) \subset \mathfrak{m}$, then $G$-invariant metrics on $G/H$ are in one-to-one correspondence with $\mathrm{AD}_H$-invariant scalar products on $\mathfrak{m}$. These exist if and only if the closure of the group $\mathrm{AD}_H(\mathfrak{m})$ is compact. Conversely, if $G/H$ admits a $G$-invariant metric, then $G$ admits a left-invariant metric which is right invariant under $H$, and the restriction of this metric to $H$ is bi-invariant. The decomposition of $\mathfrak{g}$ is given by $\mathfrak{m} = \mathfrak{h}^{\perp}$ in this case.

Note that if $H$ is connected, the condition $\mathrm{AD}_H(\mathfrak{m}) \subset \mathfrak{m}$ is equivalent to $[h, \mathfrak{m}] \subset \mathfrak{m}$. If $H$ is compact, then $\mathrm{AD}_H(\mathfrak{m})$ is compact, so a $G$-invariant metric exists. The above proposition suggests the following definitions due to Nomizu (1954).

**Definition 5.2.1 (Reductive Homogeneous space).** Let $G$ be a Lie group and $H$ a closed subgroup of $G$. Write $\mathfrak{g}, \mathfrak{h}$ their Lie algebra. We say that the homogeneous space $G/H$ is reductive if there exists some subspace $\mathfrak{m} \subset \mathfrak{g}$ such that

$$ \mathfrak{g} = \mathfrak{h} \oplus \mathfrak{m} \quad \text{and} \quad \mathrm{AD}_h(\mathfrak{m}) \subseteq \mathfrak{m} \quad \forall h \in H. $$
In this case, **m** is isomorphic to $T_oG/H$ via $d\pi_e$. In fact, **h** is the vertical subspace at $e$, and **m** is the horizontal space. The notion of reductive homogeneous space is important because it is a sufficient condition for $M = G/H$ to admit a $G$-invariant connection, whose geodesics are projections (by $\pi$) of the one-parameter subgroups of $G$. This is not sufficient however to obtain such property from the Levi-Civita connection of some metric, and an additional condition will be given in the next subsection. Accordingly with the decomposition $\mathfrak{g} = \mathfrak{h} \oplus \mathfrak{m}$, any element of the Lie algebra can be decomposed as the sum of an element of $\mathfrak{m}$ and an element of $\mathfrak{h}$. For any $x \in \mathfrak{g}$ we write this decomposition as $x = x_{\mathfrak{m}} + x_{\mathfrak{h}}$ where $x_{\mathfrak{m}} \in \mathfrak{m}$ and $x_{\mathfrak{h}} \in \mathfrak{h}$.

### Example 5.9: Invariant metric on Stiefel manifold

Recall from Example 4.17 (page 69) that with $G = SO(n)$ and

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

we obtain the Stiefel manifold $S(n, k) = G/H$, the set of orthonormal $k$-frames represented by $n \times k$ orthogonal matrices. The equivalence class of some $Q = (U, U_{\perp}) \in S(k, n)$ is $[Q] = QH = (U, U_{\perp}R)$ where $R$ is any matrix in $SO(n-k)$. The canonical projection is therefore given by $\pi : Q = (U, U_{\perp}) \mapsto U$ and the origin $o = eH$ of $G/H$ is $\pi(I_n) = I_k$. Thus we can write $\pi$ as the projection on the first $k$ columns. Moreover, $\mathfrak{g}$ is the set of skew-symmetric matrices, so that we have $\mathfrak{g} = \mathfrak{h} \oplus \mathfrak{m}$ with

$$\mathfrak{h} = \left\{ \begin{pmatrix} 0 & 0 \\ 0 & S \end{pmatrix}, S \in \mathfrak{so}(n-k) \right\}$$

$$\mathfrak{m} = \left\{ \begin{pmatrix} T & -A^{\top} \\ A & 0 \end{pmatrix}, T \in \mathfrak{so}(k) \mid A \in M_{n-k,k}(\mathbb{R}) \right\}$$

It is straightforward to check that $[\mathfrak{h}, \mathfrak{m}] \subseteq \mathfrak{m}$, and as $H$ is connected, $S(k, n)$ is a reductive homogeneous space.

As $H$ is compact, by Proposition 5.2.2 there exists a $G$-invariant metric on $S(k, n)$. Indeed, as the pullback of the Frobenius metric on $G$ is bi-invariant on $G$, hence on $H$, Proposition 5.2.1 applies and the metric on $S(k, n)$ is defined such that $\pi : Q = (U, U_{\perp}) \mapsto U$ is a Riemannian submersion. For $X = \begin{pmatrix} T & -A^{\top} \\ A & 0 \end{pmatrix} \in \mathfrak{m}$,
$d\pi_o X = \begin{pmatrix} T \\ A \end{pmatrix}$ and we obtain

$$ \begin{aligned} \langle \begin{pmatrix} T \\ A \end{pmatrix}, \begin{pmatrix} S \\ B \end{pmatrix} \rangle &\triangleq \frac{1}{2} \operatorname{tr} \left( \begin{pmatrix} T & -A^{\top} \\ A & 0 \end{pmatrix}^{\top} \begin{pmatrix} S & -B^{\top} \\ B & 0 \end{pmatrix} \right) \\ &= \frac{1}{2} \operatorname{tr}(T^{\top} S) + \operatorname{tr}(A^{\top} B) \end{aligned} \tag{5.11} $$

A more convenient way of representing a tangent vector at $U$ is by $X = US + (I - UU^{\top})A$. Then one can show that the $SO(n)$-invariant metric defined in Equation (5.11) can be written

$$ \langle X_1, X_2 \rangle_U = \operatorname{tr}(X_1(I_k - \frac{1}{2}UU^{\top})X_2). $$

For more on this metric, see Gallier and Quaintance (2020, Section 22.5).

### 5.2.3 Properties

Recall that for a reductive Lie algebra, we have the decomposition: $[X, Y] = [X, Y]_m + [X, Y]_h$. By using the formulas of the connection of an invariant metric on a Lie group (see Section 4.3) and of a Riemannian submersion (Section 5.1), one can show the following (see e.g. Gallier and Quaintance, 2020, Propositions 22.25 and 22.27).

**Theorem 5.2.1.** Let $G$ be a Lie group and $H$ a closed subgroup of $G$. If there exists an $AD_H$-invariant inner product on $\mathfrak{m}$, then the Levi-Civita connection of the induced metric on $G/H$ is given by $\nabla_X Y = -\frac{1}{2}[X, Y]_\mathfrak{m}$ and the geodesics are projections of one-parameter subgroups if and only if

$$ \langle x, [z, y]_\mathfrak{m} \rangle + \langle [z, x]_\mathfrak{m}, y \rangle = 0 \quad \forall x, y, z \in \mathfrak{m}. \tag{5.12} $$

In this case, we say that $G/H$ is naturally reductive.

This property allows to derive closed-form expressions for the geodesics in many cases, as one parameter subgroups are given by the matrix exponential in classical Lie groups. A similar behavior holds for curvature and parallel transport. Many formulas implemented in `geomstats` can be retrieved this way.
### Example 5.10: Stiefel Exponential map

Following Example 5.9 and applying Theorem 5.2.1, we can verify that the metric verifies Equation (5.12). Therefore, the exponential map is given by

$$ \mathrm{Exp}_U \left( \begin{pmatrix} S \\ A \end{pmatrix} \right) = (U, U_{\perp}) \exp \left( \begin{pmatrix} S & -A^{\top} \\ A & 0 \end{pmatrix} \right) \begin{pmatrix} I_k \\ 0 \end{pmatrix} $$

This expression can be simplified with a QR decomposition of A, and its inverse (the Logarithm) can be computed recursively Zimmermann (see 2017, for more details).

## 5.3 Symmetric spaces

To conclude this section, we briefly expose symmetric spaces, as they can be defined from a homogeneous space $G/H$ with an additional tool defined on $G$. This results in one of the most simple geometries, where the geodesics, parallel transport and curvature can be computed in closed form. We start with a more intrinsic definition and will connect with this description after. Similarly to that of homogeneous space, the structure of symmetric space does not necessarily require a Riemannian metric but only an affine connection. We first focus on the Riemannian case for simplicity, and give a few remarks on the more general case, referring to it as the affine symmetric case. This notion was introduced by Cartan (1926) who fully achieved a classification of symmetric spaces. The most complete reference is Helgason (1979), and a good exposition of the non Riemannian case is given in Kobayashi and Nomizu (1996b, Chapter XI) and Postnikov (2001, Chapters 4–10).

We first define the geodesic symmetries, the maps defined locally that revert geodesics.

**Definition 5.3.1 (Geodesic symmetry).** Let $(M, g)$ be a Riemannian manifold. Let $x \in M$ and $U \subset M$ be an open neighborhood of $x$ such that the exponential map is injective on $U$. The geodesic symmetry at $x$ is the map defined by

$$ s_x : \begin{cases} U \to M \\ y \mapsto \mathrm{Exp}_x(-\mathrm{Log}_x(y)) \end{cases} $$

It is clear that for any $x \in M$, $x$ is an isolated fixed point of $s_x$, and that $(ds_x)_x = -\mathrm{id}$, where $\mathrm{id}$ is the identity transformation of $T_x M$.

**Definition 5.3.2 (Locally symmetric space).** Let $(M, g)$ be a Riemannian manifold. $M$ is called locally symmetric if for any $x \in M$, the geodesic symmetry $s_x$ is an isometry.

**Remark 5.3.1.** This definition is valid in an affine space where the notion of isometry is replaced by the notion of affine map, i.e., that preserve the connection.
The property that defines locally symmetric spaces has direct consequences on its
curvature tensor (Kobayashi and Nomizu, 1996b, Chapter XI, Theorem 6.2).

**Theorem 5.3.1.** A Riemannian manifold $M$ is locally symmetric if and only if $\nabla R = 0$.

We shall say that the curvature of $M$ is covariantly constant, and this will simplify the
parallel transport on locally symmetric spaces.

**Definition 5.3.3 (Symmetric Space).** (M, g) is called (globally) symmetric if the geodesic symmetries are defined on the whole manifold M and are isometries.

The two notions are equivalent up to topological constrains, as stated in the following
theorem (Kobayashi and Nomizu, 1996b, Chapter XI, Theorems 6.3-6.4).

**Theorem 5.3.2.** A geodesically complete, simply connected, locally symmetric space is globally symmetric. Conversely, every globally symmetric space is geodesically complete.

We now come to the interesting theorem that relates symmetric spaces to homogeneous
spaces (Kobayashi and Nomizu, 1996b, Chapter XI, Theorem 6.5).

**Theorem 5.3.3.** The group of isometries of $M$ Isom($M$) is a Lie group that acts transitively on $M$. $M$ is thus a homogeneous space. Let $G$ be the largest connected group of isometries of $M$, and consider a reference point $o \in M$. Let $H$ its stabilizer by the action of $G$. Then $H$ is compact and $M \simeq G/H$.

Define $G$ as in the above theorem. We can define an additional structure on $G$ that will
be the essence that differentiate symmetric spaces from homogeneous spaces.

**Remark 5.3.2.**

* The composition of two geodesic symmetries belongs to $G$, and is called a transvection.

* Let $\gamma$ be a geodesic through $o$, then transvections $s_{\gamma(t)} \circ s_o$ are one-parameter subgroups of $G$.

**Theorem 5.3.4.** Let $M = G/H$ be a symmetric space, define on $G$ the map $\sigma : g \mapsto s_o \circ g \circ s_o$.

* $\sigma$ is involutive: $\sigma \circ \sigma = \text{Id}$,

* $\sigma$ is a group homeomorphism, i.e. $\sigma(g \circ h) = \sigma(g) \circ \sigma(h)$, thus an automorphism,

* The set $G^\sigma$ of fixed point of $\sigma$ is a closed subgroup, and with $G_0^\sigma$ its connected component, we have $G_0^\sigma \subset H \subset G^\sigma$ (implying $\dim(H) = \dim(G^\sigma)$).
Define $\mathfrak{t} = \{v \in \mathfrak{g}, d\sigma_e(v) = v\}$ and $\mathfrak{m} = \{v \in \mathfrak{g}, d\sigma_e(v) = -v\}$. Then $\mathfrak{t}$ coincides with the Lie algebra $\mathfrak{h}$ of $H$ and $G/H$ is naturally reductive with $\mathfrak{g} = \mathfrak{t} \oplus \mathfrak{m}$, i.e. $[\mathfrak{t}, \mathfrak{t}] \subset \mathfrak{t}$ and $[\mathfrak{t}, \mathfrak{m}] \subset \mathfrak{m}$, and we have the additional property:

$$[\mathfrak{m}, \mathfrak{m}] \subset \mathfrak{t}.$$

As a consequence (Theorem 5.2.1), the geodesics of a symmetric space are thus projections of one parameter subgroups. In fact symmetric spaces are the only homogeneous spaces with an involutive automorphism as described in Theorem 5.3.4.

**Definition 5.3.4 (Symmetric pair).** A symmetric pair is a triplet $(G, H, \sigma)$ where $G$ is a connected Lie group, $H$ a closed subgroup of $G$ and $\sigma$ is an involutive automorphism of $G$, such that its set of fixed points $G^\sigma$ satisfies $G_0^\sigma \subset H \subset G^\sigma$.

**Remark 5.3.3.**

* If there exists $a \in G$ s.t. $a^2 = e$, then $\sigma : g \mapsto a \circ g \circ a^{-1}$ is an involutive automorphism, and its set of fixed points $G^\sigma$ is a closed (normal) subgroup of $G$, so that $(G, G^\sigma, \sigma)$ forms a symmetric pair.

* The inversion $g \mapsto g^{-1}$ is an automorphism if and only if $G$ is commutative. This is very restrictive, and in general the involution is not the inversion map.

* In fact one can show (see e.g. Cheeger and Ebin, 1975, Proposition 3.37) that a simply connected group $G$ possesses an involutive automorphism $\sigma$ if and only if its Lie algebra admits a decomposition $\mathfrak{g} = \mathfrak{h} \oplus \mathfrak{m}$ with $[\mathfrak{h}, \mathfrak{h}] \subset \mathfrak{h}$, $[\mathfrak{h}, \mathfrak{m}] \subset \mathfrak{m}$ and $[\mathfrak{m}, \mathfrak{m}] \subset \mathfrak{h}$.

We now see how to recover a symmetric space from a symmetric pair. Recall that $\pi : g \in G \to gH \in G/H$ is the canonical projection of the quotient $G/H$, and that $G$ acts on $G/H$ by $g_1 \triangleright (g_2 H) = (g_1 g_2) H$. Let $o = eH = H$, and $s_o : gH \mapsto \sigma(g)H$, i.e. $s_o \circ \pi = \pi \circ \sigma$.

**Theorem 5.3.5.** Let $(G, H, \sigma)$ be a symmetric pair such that $H$ and $G_0^\sigma$ are compact. Define $\mathfrak{h} = \{v \in \mathfrak{g}, d\sigma_e(v) = v\}$ and $\mathfrak{m} = \{v \in \mathfrak{g}, d\sigma_e(v) = -v\}$. Then $M = G/H$ is naturally reductive with the decomposition $g = \mathfrak{h} \oplus \mathfrak{m}$. Furthermore, with the family of symmetries defined at any $x = gH \in M$ by

$$s_x = g \circ s_o \circ g^{-1},$$

$M$ is a globally symmetric space.

**Remark 5.3.4.** The assumption that $H$ and $G_0^\sigma$ are compact is sufficient to ensure the existence of a $G$-invariant metric on $G/H$ such that it is naturally reductive. It is not necessary however to show that $G/H$ is affine symmetric, but the connection may not be the Levi-Civita connection of any metric. If a $G$-invariant metric exists however on an affine symmetric space $M = G/H$, its Levi-Civita connection coincides with the connection of the affine symmetric structure (Kobayashi and Nomizu, 1996b, Chapter XI, Theorem 3.3).
The last theorem along with Remark 5.3.3 now make it easier to exemplify the notion of symmetric space.

### Example 5.11: Hypersphere as a symmetric space

Recall from Example 4.16 (page 68) that the hypersphere $S^d$ can be seen as the quotient of $G = SO(d+1)$ by $H \simeq SO(d)$ defined by

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

Define $J = \begin{pmatrix} -1 & 0 \\ 0 & I_d \end{pmatrix}$. Obviously $J^2 = I_{d+1}$ so that the map defined by $\sigma: P \in G \mapsto JPJ \in G$ is an involutive automorphism, and it is straightforward to check that $G^\sigma = H$. Then $(G, H, \sigma)$ is a symmetric pair and by Theorem 5.3.5, $S^d$ is a symmetric space. The Lie algebra is decomposed in $\mathfrak{so}(d+1) = \mathfrak{h} \oplus \mathfrak{m}$ with

$$\mathfrak{h} = \left\{ \begin{pmatrix} 0 & 0 \\ 0 & S \end{pmatrix}, S \in \mathfrak{so}(d) \right\}$$
$$\mathfrak{m} = \left\{ \begin{pmatrix} 0 & -u^\top \\ u & 0 \end{pmatrix}, u \in \mathbb{R}^d \right\}.$$

From the expressions of the Exp and Log map, we can compute the symmetry at $x \in S^d$. For any $y \in S^d$

$$s_x(y) = 2(\langle x, y \rangle)x - y.$$

Similarly, the upper hyperboloid and the Euclidean space $\mathbb{R}^d$ are also symmetric. Thus, all constant-curvature spaces are symmetric.

### Example 5.12: The Grassmannian as a symmetric space

The Grassmann manifold $Gr(k, n)$ (Example 4.18 page 69) is the set of $k$-dimensional subspaces of $\mathbb{R}^n$, and is identified with $O(n)/(O(k) \times O(n-k))$. Define $J = \begin{pmatrix} I_k & 0 \\ 0 & -I_{n-k} \end{pmatrix}$. Obviously $J^2 = I_n$ so that the map defined by $\sigma: P \in G \mapsto JPJ \in G$ is an involutive automorphism with fixed points

$$G^\sigma = \left\{ \begin{pmatrix} Q & 0 \\ 0 & R \end{pmatrix} \mid Q \in O(k), R \in O(n-k), \det(Q) \det(R) = 1 \right\},$$

and $G^\sigma = S(O(k) \times O(n-k))$ with $G_0^\sigma = SO(k) \times SO(n-k)$. Note that if we use $H = G_0^\sigma$, we restrict to oriented $k$-subspaces, and consider in this case the oriented
Grassmann manifold. In both cases, we have

$$
\mathfrak{h} = \left\{ \begin{pmatrix} S & 0 \\ 0 & T \end{pmatrix}, S \in \mathfrak{so}(k), T \in \mathfrak{so}(n-k) \right\}
$$
$$
\mathfrak{m} = \left\{ \begin{pmatrix} 0 & -A^\top \\ A & 0 \end{pmatrix}, A \in M_{n-k,k}(\mathbb{R}) \right\}.
$$

Then (G, H, σ) is a symmetric pair and by Theorem 5.3.5, Gr(k, n) is a symmet-
ric space. Given any P = QP_kQ^† ∈ Gr(k, n), the symmetry at P is s_P : ˜P ↦→
(QJQ^†)P̃(QJQ^†). The symmetric space structure allows to deduce many properties
of Gr(k, n), namely that it is geodesically complete, hence a complete metric space
(by the Hopf-Rinow theorem page 37), and its exponential map is surjective at all
points.

**Example 5.13: Stiefel manifold as a symmetric space**

We now give an example of homogeneous space that is not symmetric: the Stiefel
manifold St(k, n). Indeed recall that we have St(k, n) = G/H with G = SO(n) and
H = SO(n − k) and reductive decomposition so(n) = h ⊕ m with

$$
\mathfrak{h} = \left\{ \begin{pmatrix} 0 & 0 \\ 0 & T \end{pmatrix}, T \in \mathfrak{so}(n-k) \right\}
$$
$$
\mathfrak{m} = \left\{ \begin{pmatrix} S & -A^\top \\ A & 0 \end{pmatrix}, S \in \mathfrak{so}(k), A \in M_{n-k,k}(\mathbb{R}) \right\}
$$

We can check that [m, m] /∈ h, so that St(k.n) is not symmetric with this decomposi-
tion.

**Example 5.14: The Affine-Invariant metric**

Recall that SPD(n) is a homogenous space with (restricting to a connected
component) G = GL+(n) and H = SO(n) and canonical projection that coincides
with π : A ↦→ AA⊤. Define σ : A ∈ G ↦→ A−⊤. It is clear that σ is an involutive
automorphism with Gσ = H. Then (G, H, σ) is a symmetric pair and by Theorem 5.3.5,
SPD(n) is a symmetric space. From σ and π we deduce the symmetry for Σ, Λ ∈
SPD(n): sΣ(Λ) = ΣΛ−$^{1}$Σ.

The affine-invariant (AI) metric is defined as the quotient metric on G/H of the
left-invariant metric on G that coincides with the Frobenius metric at I. Its expression
is thus at any $\Sigma \in SPD(n)$, for all $V, W \in \text{Sym}(n)$

$$g_{\Sigma}(V, W) = \text{tr}(\Sigma^{-1}V\Sigma^{-1}W).$$

From the projection of one-parameter subgroups we deduce $\forall \Sigma, \Sigma_1, \Sigma_2 \in SPD(n), \forall W \in \text{Sym}(n)$

$$\text{Exp}_{\Sigma}(W) = \Sigma^{\frac{1}{2}} \exp(\Sigma^{-\frac{1}{2}} W \Sigma^{-\frac{1}{2}}) \Sigma^{\frac{1}{2}},$$

$$\text{Log}_{\Sigma_1}(\Sigma_2) = \Sigma_1^{\frac{1}{2}} \log(\Sigma_1^{-\frac{1}{2}} \Sigma_2 \Sigma_1^{-\frac{1}{2}}) \Sigma_1^{\frac{1}{2}},$$

where when not indexed, exp and log refer to the matrix operators. Finally, let

$$P_t = \Sigma^{\frac{1}{2}} \exp(\frac{t}{2} \Sigma^{-\frac{1}{2}} W \Sigma^{-\frac{1}{2}}) \Sigma^{-\frac{1}{2}}.$$

The parallel transport from $\Sigma$ along the geodesic with initial velocity $W \in \text{Sym}(n)$ of $V \in \text{Sym}(n)$ a time $t$ is (Yair et al., 2019)

$$\Pi_{0,W}^{t} V = P_t V P_t^{\top}.$$

We now focus on the case of Lie groups themselves, that can be seen as symmetric spaces. However, one must be cautious on the structure, either metric or affine that is used. Let $G$ be a connected Lie group. Consider the product group $\tilde{G} = G \times G$ with the involution $\sigma : (g, h) \mapsto (h, g)$. The subgroup of fixed points is the diagonal of $\tilde{G}$, $H = \{(g, g), g \in G\}$, and its Lie algebra is $\mathfrak{h} = \{(x, x), x \in \mathfrak{g}\} \simeq \mathfrak{g}$. We thus see from Proposition 5.2.2, that for $G = G \times G/G$ to be Riemannian symmetric hence reductive homogeneous, it must admit a bi-invariant metric. There are still three possible reductive decompositions $\mathfrak{m} \oplus \mathfrak{h}$ of the Lie algebra $\tilde{\mathfrak{g}} = \mathfrak{g} \times \mathfrak{g}$ of $\tilde{G}$:

$$\mathfrak{m} = \{(x, 0), x \in \mathfrak{g}\}$$

$$\mathfrak{m} = \{(0, x), x \in \mathfrak{g}\}$$

$$\mathfrak{m} = \{(x, -x), x \in \mathfrak{g}\}$$

With either of these, $G$ is homogeneous reductive, and each decomposition lead to a different connection, called respectively *left*, *right or mean* connection. These are known as the Cartan-Schouten connections (Lorenzi and Pennec, 2013). On can show that with the last connection, $G$ is an affine symmetric space with the symmetries $\forall g \in G$, $s_g : h \mapsto gh^{-1}g$ (Pennec et al., 2020, Theorem 5.8). If $G$ admits a bi-invariant metric, this connection is the Levi-Civita connection of the bi-invariant metric and $G$ is a Riemannian symmetric space. An example of this case is the group of rotation matrices $SO(n)$. On the contrary, $SE(n)$ does not admit any bi-invariant metric, so it does not admit a Riemannian symmetric structure that coincides with its Cartan-Schouten connection.
The properties of symmetric spaces allow to compute the geodesics and parallel transport in closed-form, by the projection of the one-parameter subgroups of the Lie group G. We rarely use explicitly this structure in geomstats, as all the results were derived and implemented case by case. The structure is however useful to compute the curvature. Indeed, using the identification of m with T<sub>o</sub>(M) induced by the restriction of dπ to m, we have for all X, Y, Z ∈ m

$$
R(X, Y)Z = -[[X, Y], Z]
$$

Moreover under topological conditions (M is simply connected and irreducible), O'Neill's formula (5.4) (page 77) simplifies and the sectional curvature is for u, v ∈ m orthogonal

$$
\kappa(u, v) = \frac{1}{2} \langle [[u, v], v], u \rangle - \frac{1}{2} \langle [[v, u], u], v \rangle.
$$

To conclude this section, symmetric spaces offer a useful framework for statistics on manifolds beyond the convenience of the closed-form solutions for the geodesics and parallel transport. Firstly, the normal distribution can be defined on all symmetric spaces, by generalizing the property that its maximum likelihood estimate coincides with the least-square problem of the Frechet mean (Said et al., 2018). The normalizing factor (or partition function) does not depend on the mean and can be computed in closed form thanks a decomposition of the Lie algebra.

Moreover, many computations such as interpolation or sampling can be performed in the subspace **m** of the Lie algebra and be projected back to the symmetric space as in Gawlik and Leok (2018), Munthe-Kaas et al. (2014), and Barp et al. (2019). Finally, there is a large theory of harmonic analysis on symmetric spaces (see Terras, 1988), and limit theorems for stochastic processes such as Brownian motion.

In this section, we focused on the common implementation of Riemannian metrics that arise from group actions. The first case is for invariant metrics on Lie groups themselves, for which the invariance allows to reformulate all the geometric operations as algebraic operations in the Lie algebra. In particular, we applied this reasoning to derive a parallel transport equation, that lead to a stable implementation in geomstats.

The second case is that of quotient metrics, that arise on the orbits of the group actions. We identified the key ingredients to a common implementation in geomstats, and focused on the Kendall shape spaces and Bures-Wasserstein metric to exemplify these. This formulation allowed to easily implement parallel transport on Kendall shape spaces, and is to the best of our knowledge the only open-source Python implementation of these spaces. The space of correlation matrices seen as the quotient of SPD manifolds by the action of diagonal matrices is also implemented in geomstats, and others spaces such as the space of positive semi-definite matrices will be implemented in the near future.

Two particular instances of quotient spaces were then introduced: reductive homogeneous spaces and symmetric spaces. They are in fact direct generalizations of Lie groups and
many operations from quotient spaces simplify. These concepts are however not used in
the present implementation, as their properties actually allows to derive more efficient
closed-form solutions.

# 6 Statistics and machine learning with Geomstats

To conclude this section, we demonstrate the use of geomstats to perform statistics on manifold-valued data. The strength of the package is that learning algorithms are defined as external estimators that take the geometric object as input. This happens thanks to the standardized interface of the classes, and ensures that all learning tools be available for all the manifolds and metrics. The obtained flexibility allows to compare the impact of the different metrics on the learning results. We give a brief introduction to geometric statistics, and the interested reader is referred to Pennec et al. (2020, Chapter 2) for the theoretical exposition. An introductory paper gathering more examples using geomstats was presented at the Scipy Conference 2020 in Austin, Texas (Miolane et al., 2020). More generally, statistical methods for objects living in stratified and infinite dimensional spaces that may not be manifolds are developed in the emerging domain of Object-Oriented Data Analysis (OODA) (Marron and Dryden, 2021).

## 6.1 Probability distributions and sampling

Given a probability measure, one can define random variables valued in a manifold $M$ as
follows.

**Definition 6.1.1.** Let $(\Omega, \mathcal{F}, \mathrm{Pr})$ be a probability space, where $\Omega$ is a sample space, $\mathcal{F}$ a $\sigma$-algebra* and $\mathrm{Pr}$ a probability measure*. Then a *random variable* in the Riemannian manifold $M$ is an $\mathcal{F}$-measurable function $X : \Omega \to M$.

Recall that a probability measure is an assignment of a size to each subset (e.g. length, area or volume in Euclidean spaces) such that the size of the entire space is $\mathrm{Pr}(\Omega) = 1$. A Riemannian metric offers a convenient framework to define probability distributions on manifolds, as it defines a volume measure, that in turn, allows to define the probability measure. This measure is used as reference-measure of probability densities.

Consider a Euclidean space $\mathbb{R}^n$, and an orthonormal basis $(e_1, \dots, e_n)$ of $\mathbb{R}^n$. Let $(a_1, \dots, a_m)$ be a set of vectors, and $A$ the matrix whose columns are formed by the coordinates of the $a_i$. Then the volume spanned by $(a_1, \dots, a_m)$ is given by $\det(A) = \sqrt{\det G}$ where $G = A^\top A$, i.e. $A_{ij} = \langle a_i, a_j \rangle$ (see Figure 6.1.)

<footnote>*Defined in Appendix A.</footnote>
<figure><img src="image_19.png" /><figcaption>Figure 6.1: Volume spanned by three vectors in $\mathbb{R}^3$.</figcaption></figure>

Now in a Riemannian manifold $(M, g)$, for any $x \in M$, one may choose an orthonormal basis $(a_1, \dots, a_d)$ of $T_x M$ with respect to $g$. $g$ is then represented by the SPD matrix $G_{ij} = \langle a_i, a_j \rangle$, so that the volume spanned by these vectors is $\sqrt{\det G}$.

**Definition 6.1.2 (Riemannian volume form).** An oriented Riemannian manifold $(M, g)$ has a natural volume form defined by $d\mathrm{Vol}(x) = \sqrt{\det(g(x))}dx$.

The volume form provides the way to define integrals on $M$, i.e. it defines a measure, so that some probability distributions can be expressed by their densities with respect to that measure.

**Definition 6.1.3.** The random variable $X$ has density $f$ if $\forall \mathcal{X} \in \mathcal{B}(M)$, $\mathrm{Pr}(X \in \mathcal{X}) = \int_{\mathcal{X}} f(y)dv(y)$ and $\mathrm{Pr}(M) = 1$.

### Example 6.1: Uniform distribution

Let $M$ be a compact Riemannian manifold, for example the hypersphere, the special orthogonal group of the Grassmann manifold. Then a uniform distribution on $M$ has density

$$f(x) = \frac{1}{\mathrm{Vol}(M)}.$$

### Example 6.2: Gaussian distribution

Let $M$ be a Riemannian Symmetric space. A Gaussian distribution with mean and concentration $(\mu, \Gamma)$ is defined by the density

$$f(x) = \alpha(\Gamma, \mu) \exp(-\log_{\mu}(x)^{\top}\Gamma\log_{\mu}(x)).$$

It is the entropy maximizing distribution (Pennec, 2006), and in the isotropic, the maximum likelihood estimator of $\mu$ coincides with the Frechet Mean (Pennec et al., 2020, Section 2.5.1). The normalizing constant is explicitly computed in negative
curvature spaces (Said et al., 2018).

**Sampling** A common task in statistics is to draw samples from a probability distribution. This might be motivated by inference tasks where the posterior distributions has constrained parameters, in testing goodness of fit for exponential families, or to generate samples of data to test our learning algorithms (Diaconis et al., 2013; Barp et al., 2019, and references therein). However, on a manifold this does not reduce to sampling from usual distributions even when a parametrization of the manifold is available, as the curvature of the space may deform or stretch the densities. This is illustrated in the following examples.

**Example 6.3: Uniform distribution on the sphere**

Consider the unit sphere $S^2 \subset \mathbb{R}^3$, with the spherical coordinates from the north pole $e_0 = (0,0,1): x = \sin(\phi)\cos(\theta), y = \sin(\phi)\sin(\theta), z = \cos(\phi)$. The volume element (here the area) is $d\text{Vol}(\theta, \phi) = \sin(\phi)d\theta d\phi$. A naive attempt at sampling from the uniform distribution on the sphere would be to sample $\theta$ uniformly in $[0, 2\pi)$ and $\phi$ uniformly in $[0, \pi)$.

However, near $\phi = 0$ and $\phi = \pi$, as $\sin(\phi)$ there is less density than with respect to the usual measure of the 2d-plane. We thus expect points to accumulate in a non-uniform way around the poles of the sphere.

In fact, if $Y$ be a Gaussian vector in $\mathbb{R}^d$ with mean 0 and covariance identity. Then $X = \frac{Y}{\|Y\|}$ is uniformly distributed in $S^d$. The two sampling strategies are tested with `geomstats` and shown on Figure 6.2.

```python
space = Hypersphere(2)
n Samples = 5000
uniform_param = gs.randomrand(
    n Samples, 2) * gs.pi * gs.array([1., 2.])[None, :]
naive_sample = space.spherical_to_extrinsic(uniform_param)
uniform_sample = space.randomuniform(n_sample)
```
<figure><img src="image_20.png" /><figcaption>**Figure 6.2:** Comparison of the naive sampling from a uniform distribution on the spherical coordinates (left) and a more appropriate scheme that respects the isotropy of the sphere (right). In the naive case, points accumulate near the poles and are sparser around the equator.</figcaption></figure>

### Example 6.4: Uniform distribution on $SO(n)$

Let $SO(n)$ be the group of unit determinant orthogonal matrices of size $n \times n$. Let $Y$ be a matrix of size $n \times n$ with independent standard normal distribution on each coefficient. The QR decomposition of $Y$ defines $Y = XR$ where $X \in SO(n)$ and $R$ upper triangular. Then $X$ is uniformly distributed on $SO(n)$ w.r.t. to the Haar measure, which coincides with the Riemannian measure (Eaton, 1983).

### Example 6.5: Uniform distribution on Stiefel manifold

Recall the Stiefel manifold is the set of orthonormal $k$-frames in $\mathbb{R}^n$. It can be represented by the set of matrices $\{U \in \mathbb{R}^{n \times k} | U^\top U = I_k\}$. Let $Z$ be a matrix with i.i.d standard normal distribution on the entries. Then $X = Z(Z^\top Z)^{-1/2}$ is uniformly distributed (Chikuse, 2003, Theorem 2.2.1). Note that $X$ is a factor of the polar decomposition of $Z$, representing the *orientation* of $Z$.

### Example 6.6: Uniform distribution on Grassmannians

The Grassmann manifold is the set of $k$-dim subspaces of $\mathbb{R}^n$. It can be represented by the set of projection matrices of size $n \times n$, i.e. symmetric rank-$k$ $P$ such that $P^2 = P$. Let $Z$ be a matrix with i.i.d standard normal distribution on the entries. Then $X = Z(Z^\top Z)^{-1/2}Z$ is uniformly distributed (Chikuse, 2003, Theorem 2.2.2).
Note that $X = YY^\top$ for $Y$ uniformly distributed on the Stiefel manifold.

These examples are particular cases where simple recipes were available to sample from the uniform distribution. They are implemented in `geomstats`. Sampling from non-uniform distribution is usually intractable and we resort to simulation methods from the literature. In particular, in `geomstats` we focused on rejection sampling, that consists in sampling from a proposal distribution whose density is greater than the target density, and accepting the samples with probability corresponding to the ratio of the two densities<fnref n="2" />. We do not detail this procedure in this monograph but refer to Wood (1994) and Hauberg (2018) for examples of rejection sampling algorithms for non-uniform distributions on the hypersphere.

## 6.2 Distance-based algorithms

A large class of machine learning algorithms only require computing a matrix where each entry is the distance between a pair of samples. This is the case for examples of the nearest neighbors methods or hierarchical clustering. These can readily be used by computing the distance matrix with a `geomstats` metric. We give below a toy example on the sphere, with data generated with two von Mises distributions (Example 6.7).

**Example 6.7: Hierarchical clustering on the sphere**

The coaches of the French skiing Olympic team studies the solid angles of inclination of their skier during a slalom. To first filter left and right turns, they use a hierarchical clustering algorithm on the sphere of solid angles.

```python
space = Hypersphere(2)
metric = space.metric
```

To test the method, they first generate a toy dataset with von Mises-Fisher distributions

```python
n Clusters = 2
n Samples = 50

left = sphere.random_von_mises_fisher(kappa=10, n Samples=n samples)
right = -sphere.random_von_mises_fisher(kappa=10, n samples=n samples)
dataset = gs.concatenate((left, right), axis=0)
```

A `geomstats` model with the Riemannian distance can then be used to compute the distances. We use `n clusters=None`, `distance_threshold=0` in order to compute the full tree and plot the dendrogram shown on the right panel of Figure 6.3 (as in the

<footnote n="2">²See the Wikipedia page for more details.</footnote>
scikit-learn example).

```r
model = AgglomerativeHierarchicalClustering(
  n clusters=None, distance=metric.dist, distance_threshold=0)
model.fit(dataset)
plot_dendrogram(model)
```

The cluster assignments can also be computed by setting `n_cluster=2`. These are shown on Figure 6.3 (Left).

```r
model = AgglomerativeHierarchicalClustering(n clusters, distance=metric.dist)
model.fit labels)
clustering labels = model.labels_
```

<figure><img src="image_21.png" /><figcaption>**Figure 6.3:** Toy experiment of a hierarchical clustering experiment. The data is sampled from two von Mises-Fisher distributions (left). A dendrogram is computed from the distance matrix (right).</figcaption></figure>

## 6.3 The Fréchet mean

Let $(M, g)$ be a Riemannian manifold and let $x_1, \dots, x_n \in M$ be independent identically distributed (i.i.d.) sample data points. The Euclidean sample mean $\bar{x} = \frac{1}{n} \sum_{i=1}^n x_i$ is not defined unless we consider $M \subset \mathbb{R}^N$, but in this case $\bar{x}$ may not lie on $M$.

**Definition, Existence, Uniqueness** In fact, the mean is characterized by the property that it minimises the sum of squared distances to the data points. This property was used by Fréchet (1948) to generalize the notion of mean to metric spaces, and later by Hermann Karcher in Riemannian manifolds. For a historical note and corresponding references, we refer to Karcher (2014).

**Definition 6.3.1 (Fréchet mean).** Let $(M, g)$ be a Riemannian manifold with Riemannian
distance function $d$ and let $x_1, \dots, x_n \in M$ be an i.i.d. data set. The sample Fréchet mean is defined as the set of minimizers of the sum-of-squared distances, i.e.,

$$ \bar{x} = \arg \min_{x \in M} \sum_{i=1}^{n} d(x_i, x)^2. $$

With some nuance in the requirements on the minimum (local or global), the Fréchet mean is also referred to as Karcher mean or Riemannian barycenter, or even Wasserstein barycenter in optimal transport. Note that as it is defined by an optimization, the Fréchet mean may not exist, or not be unique, depending on the properties of the distance function. For instance, the completeness of the metric is a sufficient condition to guarantee the existence of the Fréchet mean of a finite set of points (Pennec et al., 2020, Chapter 2).

**Theorem 6.3.1.** Let $M$ be a complete metric space. Then the Fréchet mean of any finite set of points $x_1, \dots, x_n$ exists.

On the other hand, uniqueness depends on the convexity of the distance function and can be related to the sign of the curvature (Karcher, 1977; Kendall, 1990; Afsari, 2011).

**Theorem 6.3.2.** Let $M$ be a complete Riemannian manifold with sectional curvature bounded above by some $\delta$, and let $\mathrm{inj}(M)$ be its injectivity radius. If samples $x_1, \dots, x_n \in M$ are contained in a geodesic ball of radius

$$ r = \frac{1}{2} \min \left( \mathrm{inj}(M), \frac{\pi}{\sqrt{\delta}} \right) $$

then the Fréchet mean $\bar{x}$ is unique.

Note that in the above, $\frac{\pi}{\sqrt{\delta}}$ is interpreted as $\infty$ if $\delta \le 0$. This means that in complete manifolds of non-positive curvature, the Fréchet mean is always uniquely defined, as in the Euclidean case. On the other hand, this is not the case in spaces of positive curvature. For example, on the sphere, Theorem 6.3.2 ensures that the mean is unique if all the data points lie on the same open hemisphere, but the simple data set constituted by two antipodal points ($\{x, -x\}$ for any $x \in S^d$) exemplifies non-uniqueness of the mean: a great circle in this case.

The asymptotic properties of the Fréchet mean were studied in Bhattacharya and Patrangenaru (2005), where a law of large numbers and a central limit theorem are established, and prove the relevance of this notion of mean.

**Characterization and Estimation** Suppose that the data points $x_i$, $i = 1, \dots, n$ lie in a ball $B$ of radius $r < \frac{\mathrm{inj}(M)}{2}$ as in Theorem 6.3.2. Then the squared distance is convex and is obtained by measuring the length of a minimizing geodesic and can be written $d^2(x, y) = \| \mathrm{Log}_x(y) \|_x^2$ for any $x, y \in B$. One can show that the gradient of $y \mapsto d^2(x, y)$
is $-2 \log_y(x) \in T_y M$ (see the first variation formula in Lafontaine et al., 2004, Theorem 3.31).

**Remark 6.3.1.** Recall that the gradient of a function $f: \mathbb{R}^d \to \mathbb{R}$ is defined as the adjoint to its differential, i.e.

$$ \forall (x, v) \in TM, \quad \langle \mathrm{grad}_x(f), v \rangle \triangleq df_x v. $$

This definition extends to functions defined on Riemannian manifolds, but the gradient thus depends on the metric, and is sometimes referred to as the Riemannian gradient. It is also named natural gradient in information geometry, in the case of manifolds of parameters of families of distributions with the Fisher-Rao metric. If the manifold is embedded in some $\mathbb{R}^N$ and equipped with the pullback metric, it can be obtained by projecting the usual gradient to the tangent space at the point where the function is being differentiated, similarly to the Levi-Civita connection.

By definition, the Fréchet mean must be a critical point of the sum-of-squared distances function, so that it verifies

$$ \sum_{i=1}^{n} \mathrm{Log}_{x}(y) = 0. \qquad (6.1) $$

A simple strategy to estimate the sample Fréchet mean is thus to use a fixed-point iterative algorithm until Equation (6.1) is verified. An update at iteration $t$ is performed along a geodesic:

$$ x^{t+1} = \mathrm{Exp}_{x^t}(\gamma \sum_{i=1}^{n} \mathrm{Log}_{x^t}(y)), $$

where $\gamma$ is a step size. Higher-order methods use adaptive step size, approximates of the Hessian function or its link with curvature to improve the performance of the estimation. Some of these are implemented in the `FrechetMean` module of `geomstats`.

### Example 6.8: Fréchet mean on the sphere

The following is an illustrative example with simulated data. A scientific expedition of 15 boats is collecting samples from the Pacific ocean to measure the quantity of micro-plastics in the water. Every month, they gather all the samples in one boat that goes back to their home harbour for the analyses. As they have been wandering around the Pacific for some time, we sample their positions from a spherical normal distribution around their last meeting point.

```python
space = Hypersphere(2)
n Samples = 15
last_meeting_point = gs.array([0., 1., 0.])
samples = space.random_riemannian_normal(
    mean=last_meeting_point, precision=10, n Samples=15)
```
The geometers onboard compute the Fréchet mean of their positions to choose the meeting spot that minimises the sum-of-squared distances to this gathering (supposing no wind or currents).

estimator = FrechetMean(space metric)
estimator.fit samples)
new_meeting_point = estimator.estimate_

Their position (black) and the next meeting point (red) are displayed in Figure 6.4. The trajectories they have to follow are drawn in green, and are geodesics.

<figure><img src="image_22.png" /><figcaption>Figure 6.4: Fréchet mean on the sphere</figcaption></figure>

The Fréchet mean $\bar{x}$ is key to define many other statistical and learning tools, e.g. the sample covariance matrix:

$$ \Sigma = \frac{1}{n} \sum_{i=1}^{n} \mathrm{Log}_{\bar{x}} x_i \otimes \mathrm{Log}_{\bar{x}} x_i. \qquad (6.2) $$

One can then parametrize densities from their mean and covariance, as in Chevallier and Guigui (2020). The proposed model allows to define a simple density estimation scheme.

One can also generalize the K-means clustering algorithm and the minimum distance to mean (MDM) classification algorithm with Fréchet means. These are implemented in the classes RiemannianKMeans and RiemannianMinimumDistanceToMeanClassifier in geomstats. Moreover, one can use the Fréchet mean to linearize the data by lifting it to the tangent space at the Fréchet mean, i.e. to consider the transformed data set ($\tilde{x}_1, \dots, \tilde{x}_n$) where

$$ \tilde{x}_i = \mathrm{Log}_{\bar{x}}(x_i). $$
This is equivalent to using a Taylor expansion at the first order of the Riemannian exponential map around 0: $x_i = \text{Exp}_{\bar{x}}(\tilde{x}_i) = \bar{x} + \tilde{x}_i + O(t^2)$. As $T_{\bar{x}}M$ is a vector space, all the usual statistical and machine learning tools can be used off-the-shelf on the transformed data set. This is implemented as a Transformer from the scikit-learn package in the ToTangent class of geomstats, and can be used in scikit-learn's pipeline.

## 6.4 Generalizations of PCA

The other fundamental tools to analyse data in vector spaces are the sample covariance matrix, and Principal Component Analysis (PCA). The aim of PCA is to find the sequence of subspaces such that data projected on these subspaces have maximum variance, or equivalently minimum reconstruction error (i.e. sum of squared distance to the original points). This equivalence is no longer true in Riemannian manifolds, as Pythagorean theorem is not true.

We describe the principal geodesic analysis (PGA), that aims at minimizing the reconstruction error when projecting data on a geodesic submanifold (Fletcher et al., 2004). In the simplest forward fashion, the mean is first computed, and the first component is a geodesic from the mean. This is close to a slightly different procedure called Geodesic PCA (Huckemann et al., 2010). The projection of a data point $x$ on a geodesic $\gamma$ with initial velocity $v$ at the mean $\bar{x}$ is defined by

$$ \pi_{\bar{x},v}(x) = \arg \min_{t \in \mathbb{R}} d^2(x, \text{Exp}_{\bar{x}}(tv)). $$

There is no guarantee that this projection exists, but if the data is not too spread out, one can hope that there exists a convex neighborhood where the exponential map is injective that contains $x$ and a portion of the geodesic. The projection is solved by a gradient descent, as in the case of the Fréchet mean, where the gradient of the exponential is computed by automatic differentiation.

The next step is to minimize the overall reconstruction error, given a dataset $x_1, \dots, x_n$, we look for the initial velocity of a geodesic as:

$$ v^* = \arg \min_{v \in T_{\bar{x}}M} \sum_{i=1}^{n} d^2(x_i, \pi_{\bar{x},v}(x_i)). \qquad (6.3) $$

This time computing the gradient of the objective function seems more complicated, as $n$ minimization problems are already solved by a gradient descent to evaluate the function. However Ablin et al. (2020) show that the gradient of a function defined by a minimum could be computed efficiently by automatic differentiation of the gradient descent approximation. We exemplify this in the example of the sphere. Of course this example is particularly simple as both Exp and Log maps are computed with closed-form solutions whose derivatives are
known, and some work remains to use automatic differentiation when the Log map itself is
hard to compute.

**Example 6.9: Principal geodesic analysis on earth**

The following is an illustrative example with simulated data. The boats of the scientific exploration of Example 6.8 are in fact too busy collecting data and cannot make route to the Fréchet mean. Thankfully, the geometers onboard are expert users of `geomstats` and compute a PGA of their positions, so that one boat can visit the rest of the fleet by sailing on the principal geodesic. This is illustrated on Figure 6.5.

<figure><img src="image_23.png" /><figcaption>**Figure 6.5:** Example of PGA on the sphere. The target data are the black dots, the geodesics that realize the projection of each sample point to the principal geodesic are shown in green, the first geodesic subspace or principal geodesic is in gray, and the mean in red.</figcaption></figure>

A general simplification if the data is not too far from the mean is to approximate the projection operator by the Log map: $\operatorname{Log}_{\bar{x}}(\pi_{\bar{x}}, v(x_i)) \approx \operatorname{Log}_{\bar{x}}(x_i)$. With this approximation, PGA resumes to computing the logarithms of all the data points and computing a usual PCA of the obtained tangent vectors. This procedure is more generally called *tangent PCA* and is implemented in the `TangentPCA` class of `geomstats`.

**6.5 Geodesic Regression**

Similarly to the mean, the linear regression can be generalized to *geodesic regression* on Riemannian manifolds by solving a least-square fitting problem. Given target points $y_1, \dots, y_n \in M$, and data $t_1, \dots, t_n \in \mathbb{R}$, we seek the geodesic that best approximates the
data:

$$ \min_{(p,v) \in TM} \sum_{i=1}^{n} d^2(\text{Exp}_p(t_i v), y_i), \qquad (6.4) $$

where d and Exp are the Riemannian distance and exponential maps. When the metric is Euclidean, this coincides with the usual linear regression problem. However, there is no closed-form solution in general and the problem must be solved by optimization.

To differentiate the objective function of Equation (6.4), one needs to compose the gradients of the squared distance with that of the exponential map. The gradient of the squared distance is proportional to the Riemannian logarithm, as noted above. On the other hand, the gradient of the Exp map is usually computed via Jacobi fields. However, to avoid implementing those in geomstats we chose to leverage automatic differentiation tools to compute the extrinsic gradient, and to project it to the right tangent space, as explained in Remark 6.3.1.

To solve the optimization problem, either a Riemannian gradient descent, or an extrinsic one with scipy's solver is used.

### Example 6.10: Geodesic Regression above earth

The following is an illustrative example with simulated data. The European Space Agency sends a maintenance mission on their geo-stationary satellites. As these are all at the same altitude, they lie on a sphere, and the geometers seek an optimal trajectory.

The maintenance should be done as close as possible to the end of the fuel tanks, which are known from previous missions. One vessel will leave earth, reach the first point of the trajectory called point $\gamma$, and split in two parts that will each go in opposite directions given by $\beta \in T_{\gamma}S^2$.

We use data from the previous mission to generate random satellite positions:

```python
space = Hypersphere(2)
n Samples = 50
data = gs.randomrand(n Samples)
data -= gs.mean(data)

previous gamma = space.randomuniform()
beta = space.to_tangent(5. * gs.randomrand(3), previous gamma)
target = space.metrics.exp(data[:, None] * beta, previous gamma)
```

And add some noise because satellites did not stay on the previous mission's trajectory

```python
normal noise = gs.random.normal(size=(n samples, 3))
noise = space.to_tangent(normal noise, target) / gs.pi / 2
target = space.metrics.exp(noise, target)
```

The optimal trajectory is computed by fitting a geodesic regression model with the
Riemannian gradient descent:

```
gr = GeodesicRegression(space, algorithm='riemannian')
gr.fit(data, target, compute_training_score=True)
gamma_new, beta_new = gr.intercept_, gr.coef_
```

We can measure the mean squared error (MSE) with respect to the previous $\gamma$ and $\beta$, and the determination coefficient from the noise level. We obtain MSEs around $10^{-3}$ and $R^2 = 0.96$. The results are shown on Figure 6.6.

<figure><img src="image_24.png" /><figcaption>**Figure 6.6:** Example of geodesic regression on the sphere. The target data are the blue dots, the fitted values are in green, the regression geodesic in black, and the intercept ($\gamma$) in red.</figcaption></figure>

# 7 Conclusion

In this monograph, we introduced notions of differential geometry that form the essential building blocks of geometric statistics. Our developments target an audience with a general background in mathematics. We focused on embedded manifolds of $\mathbb{R}^N$, and highlighted the different ways of defining a manifold and a Riemannian metric. We detailed the notions of curvature, geodesic, distance and parallel transport. We also presented the concepts of Lie groups, group actions and quotient space, with an emphasis on Riemannian metrics that are invariant to a group action.

These differential geometric notions are key to the architecture of the `geomstats` library. They drove our recent contributions to the package, that aimed at making it more faithful to mathematical theory, more robust and more modular. We exemplified the geometric concepts with the most common manifolds encountered in mathematical text books as well as in applications. Our examples include code snippets using `geomstats` to demonstrate how
to leverage differential geometry in practical use cases.

Lastly, we gave an introduction to geometric statistic tools such as the Fréchet mean, principal geodesic analysis and geodesic regression. We illustrated these on toy examples using synthetic datasets on manifolds. These geometric statistical tools are implemented in the `geomstats` package with a common high-level interface, following the Scikit-Learn syntax. Consequently, geometric statistics become available to any data scientist on a wide variety of practical problems.

We hope that the concepts presented here will drive scientists to use, and contribute to, geometric statistics with the `geomstats` library. Future developments of the `geomstats` package will integrate additional statistical methodologies published by the geometric statistics community, such as wrapped Gaussian processes, multivariate and polynomial geodesic regression. We will also extend the scope of `geomstats` and include computational methods for information geometry, a field at the intersection of geometry and statistics closely related to geometric statistics. Another module on stratified spaces such as graph and tree spaces is currently being developed, and raises fundamental methodological questions as these spaces are not smooth manifolds but unions of smooth manifolds equipped with a distance function. Together, these advances aim to provide mathematically-grounded foundations for computational geometric statistics.

# 8 Acknowledgment

This work was funded by the ERC grant Nr. 786854 G-Statistics from the European Research Council under the European Union's Horizon 2020 research and innovation program. It was also supported by the French government through the 3IA Côte d'Azur Investments ANR-19-P3IA-0002 managed by the National Research Agency. The authors are grateful to Adele Myers for her corrections and suggestions to improve this monograph.
List of Examples

<table><tr><td>2.1</td><td>Hypersphere</td><td>5</td></tr><tr><td>2.2</td><td>Hyperbolic space</td><td>6</td></tr><tr><td>2.3</td><td>Special Orthogonal group</td><td>6</td></tr><tr><td>2.4</td><td>Stiefel manifold</td><td>7</td></tr><tr><td>2.5</td><td>Kendall size-and-shape space</td><td>7</td></tr><tr><td>2.6</td><td>Product manifold</td><td>9</td></tr><tr><td>2.7</td><td>Cusp and Node</td><td>10</td></tr><tr><td>2.8</td><td>Tangent space of the hypersphere</td><td>12</td></tr><tr><td>2.9</td><td>Tangent space of the hyperbolic space</td><td>13</td></tr><tr><td>2.10</td><td>Tangent space of Stiefel manifold</td><td>13</td></tr><tr><td>2.11</td><td>Tangent space of SO(n)</td><td>13</td></tr><tr><td>2.12</td><td>Implementation of the hypersphere</td><td>20</td></tr><tr><td>2.13</td><td>Implementation of the Stiefel manifold</td><td>21</td></tr><tr><td>2.14</td><td>Implementation of the Poincaré ball</td><td>22</td></tr><tr><td>3.1</td><td>Euclidean metric</td><td>26</td></tr><tr><td>3.2</td><td>Product metric</td><td>26</td></tr><tr><td>3.3</td><td>Metric on the hypersphere</td><td>27</td></tr><tr><td>3.4</td><td>Hyperbolic metric</td><td>28</td></tr><tr><td>3.5</td><td>Frobenius metric</td><td>29</td></tr><tr><td>3.6</td><td>Geodesics of the hypersphere</td><td>35</td></tr><tr><td>3.7</td><td>Completeness of the half-plane</td><td>37</td></tr><tr><td>3.8</td><td>Curvature of the Euclidean space</td><td>40</td></tr><tr><td>3.9</td><td>Sectional curvature of S<sup>d</sup> and H<sub>+</sub><sup>d</sup></td><td>42</td></tr><tr><td>4.1</td><td>General Linear group</td><td>45</td></tr><tr><td>4.2</td><td>Vector fields</td><td>46</td></tr><tr><td>4.3</td><td>Matrix algebra</td><td>46</td></tr><tr><td>4.4</td><td>Implementation of GL(n)</td><td>47</td></tr><tr><td>4.5</td><td>General linear Lie algebra</td><td>48</td></tr><tr><td>4.6</td><td>Implementation of SO(n)</td><td>49</td></tr><tr><td>4.7</td><td>Implementation of SE(n)</td><td>51</td></tr><tr><td>4.8</td><td>Exponential map of GL(n)</td><td>53</td></tr><tr><td>4.9</td><td>Curve in SE(2)</td><td>55</td></tr><tr><td>4.10</td><td>Invariance metric on SO(n)</td><td>56</td></tr><tr><td>4.11</td><td>Invariance metric on SE(n)</td><td>57</td></tr><tr><td>4.12</td><td>Structure constants on SO(3)</td><td>58</td></tr><tr><td>4.13</td><td>Structure constants on SE(3)</td><td>58</td></tr><tr><td>4.14</td><td>Geodesics on SE(2)</td><td>61</td></tr></table>
<table>
  <tbody>
    <tr>
      <td>4.15 Parallel transport in SE(3)</td>
      <td>64</td>
    </tr>
    <tr>
      <td>4.16 Hypersphere as a homogenous space</td>
      <td>68</td>
    </tr>
    <tr>
      <td>4.17 Stiefel manifold as a homogenous space</td>
      <td>69</td>
    </tr>
    <tr>
      <td>4.18 Grassmannian as a homogenous space</td>
      <td>69</td>
    </tr>
    <tr>
      <td>4.19 SPD matrices as a homogenous space</td>
      <td>70</td>
    </tr>
    <tr>
      <td>5.1 Bures-Wasserstein metric</td>
      <td>73</td>
    </tr>
    <tr>
      <td>5.2 Horizontal lift of SPD matrices</td>
      <td>74</td>
    </tr>
    <tr>
      <td>5.3 Bures-Wasserstein geodesics</td>
      <td>76</td>
    </tr>
    <tr>
      <td>5.4 Curvature of the Bures-Wasserstein metric</td>
      <td>77</td>
    </tr>
    <tr>
      <td>5.5 Bures-Wasserstein distance</td>
      <td>79</td>
    </tr>
    <tr>
      <td>5.6 Implementation of the BW metric</td>
      <td>85</td>
    </tr>
    <tr>
      <td>5.7 Kendall shape metric</td>
      <td>87</td>
    </tr>
    <tr>
      <td>5.8 Full rank correlation matrices</td>
      <td>88</td>
    </tr>
    <tr>
      <td>5.9 Invariant metric on Stiefel manifold</td>
      <td>92</td>
    </tr>
    <tr>
      <td>5.10 Stiefel Exponential map</td>
      <td>94</td>
    </tr>
    <tr>
      <td>5.11 Hypersphere as a symmetric space</td>
      <td>97</td>
    </tr>
    <tr>
      <td>5.12 The Grassmannian as a symmetric space</td>
      <td>97</td>
    </tr>
    <tr>
      <td>5.13 Stiefel manifold as a symmetric space</td>
      <td>98</td>
    </tr>
    <tr>
      <td>5.14 The Affine-Invariant metric</td>
      <td>98</td>
    </tr>
    <tr>
      <td>6.1 Uniform distribution</td>
      <td>102</td>
    </tr>
    <tr>
      <td>6.2 Gaussian distribution</td>
      <td>102</td>
    </tr>
    <tr>
      <td>6.3 Uniform distribution on the sphere</td>
      <td>103</td>
    </tr>
    <tr>
      <td>6.4 Uniform distribution on SO(n)</td>
      <td>104</td>
    </tr>
    <tr>
      <td>6.5 Uniform distribution on Stiefel manifold</td>
      <td>104</td>
    </tr>
    <tr>
      <td>6.6 Uniform distribution on Grassmannians</td>
      <td>104</td>
    </tr>
    <tr>
      <td>6.7 Hierarchical clustering on the sphere</td>
      <td>105</td>
    </tr>
    <tr>
      <td>6.8 Fréchet mean on the sphere</td>
      <td>108</td>
    </tr>
    <tr>
      <td>6.9 Principal geodesic analysis on earth</td>
      <td>111</td>
    </tr>
    <tr>
      <td>6.10 Geodesic Regression above earth</td>
      <td>112</td>
    </tr>
  </tbody>
</table>
List of Figures

<table><tr><td>2.1</td><td>Representation of a manifold.</td><td>6</td></tr><tr><td>2.2</td><td>Transition maps.</td><td>9</td></tr><tr><td>2.3</td><td>Nodal cubic M<sub>1</sub></td><td>10</td></tr><tr><td>2.4</td><td>Cuspidal curve M<sub>2</sub></td><td>10</td></tr><tr><td>2.5</td><td>Definition of a smooth curve on a manifold.</td><td>11</td></tr><tr><td>2.6</td><td>Tangent vector and tangent space.</td><td>11</td></tr><tr><td>2.7</td><td>Example of an integral curve</td><td>16</td></tr><tr><td>2.8</td><td>Example of a flow</td><td>16</td></tr><tr><td>2.9</td><td>Architecture of the manifolds of geomstats.</td><td>23</td></tr><tr><td>3.1</td><td>Parallel vector field</td><td>31</td></tr><tr><td>3.2</td><td>Impact of the curvature on the distance function</td><td>43</td></tr><tr><td>4.1</td><td>One-parameter subgroups of SE(2).</td><td>55</td></tr><tr><td>4.2</td><td>Geodesics and one-parameter subgroups of SE(2).</td><td>61</td></tr><tr><td>4.3</td><td>Visualization of geodesics on SE(2).</td><td>62</td></tr><tr><td>4.4</td><td>Error curves for the parallel transport on SE(3)</td><td>64</td></tr><tr><td>5.1</td><td>Error curves for parallel transport on Kendall spaces</td><td>83</td></tr><tr><td>6.1</td><td>Volume spanned by three vectors in R<sup>3</sup>.</td><td>102</td></tr><tr><td>6.2</td><td>Samples from the uniform distribution on the sphere</td><td>104</td></tr><tr><td>6.3</td><td>Hierarchical clustering on the sphere</td><td>106</td></tr><tr><td>6.4</td><td>Fréchet mean on the sphere</td><td>109</td></tr><tr><td>6.5</td><td>Principal Geodesic Analysis on the sphere</td><td>111</td></tr><tr><td>6.6</td><td>Geodesic regression on the sphere.</td><td>113</td></tr></table>
A Lexicon

**Definition A.0.1 (Topology).** Let $M$ be a set and $\mathcal{P}(M)$ the set of subsets of $M$. Then a topology on $M$ is a set $\Theta \in \mathcal{P}(M)$ such that

• $\emptyset \in \Theta$ and $M \in \Theta$

• $\forall U, V \in \Theta, U \cap V \in \Theta$

• $C \subseteq \Theta \implies \bigcup_{U \in C} U \in \Theta$

The sets in $\Theta$ are called *open sets*, and a set $S$ is said to be *closed* if and only if $M \setminus S \in \Theta$. Such a pair $(M, \Theta)$ is called a *topological space*.

**Definition A.0.2 (Hausdorff).** A topological space $(M, \Theta)$ is said to be *Hausdorff* if, for any two distinct points $p, q \in M$, there exist open neighborhoods of $p$ and $q$ with empty intersection.

**Definition A.0.3.** A topological space $(M, \Theta)$ is second-countable if there exists some countable collection $\mathcal{A} = \{U_i\}_{i \in \mathbb{N}}$ of open sets of $M$ such that any open set can be written as a union of elements of $\mathcal{A}$

**Definition A.0.4 (Connected).** A topological space $(M, \Theta)$ is said to be *connected* unless there exist two non empty open sets $A, B \in \Theta$ such that $A \cap B \neq \emptyset$ and $M = A \cup B$.

Equivalently, $M$ is connected if and only if the only subsets that are both open and closed are $M$ itself and the empty set $\emptyset$.

**Definition A.0.5 (Group).** A group is a couple $(G, \cdot)$ where $G$ is a nonempty set, and $\cdot: G \times G \to G$ is a map such that

• $\exists e \in G, \forall g \in G, e \cdot g = g \cdot e = g,$

• $\forall g \in G, \exists h, g \cdot h = h \cdot g = e$. We define the inverse of $g$ for $\cdot$ as $g^{-1} = h$ in this case.

**Definition A.0.6 (Homomorphism).** Les $(G, \cdot)$ and $(H, \bullet)$ be two groups, and $f: G \to H$ be a map. Then $f$ is a *homomorphism* if for any $x, y \in G$, $f(x \cdot y) = f(x) \cdot f(y)$.

**Definition A.0.7 (Algebra).** An *algebra* over a field $K$ is a vector space $(A, +, \cdot)$ over $K$ equipped with a bilinear multiplicative law $\otimes: A \times A \to A$ such that

• (distributivity) $\forall x, y, z \in A, (x+y) \otimes z = x \otimes z + y \otimes z$ and $z \otimes (x+y) = z \otimes x + z \otimes y$

• (compatibility with scalars) $\forall a, b \in K, \forall x, y \in A, (ax) \otimes (by) = (ab)x \otimes y$.

**Definition A.0.8 (Injective-Surjective-Bijective map).** Let $E, F$ be two sets and $f: E \to F$ a map between $E$ and $F$. Then we say that

• $f$ is *injective* if for every $x, x' \in E, x \neq x' \implies f(x) \neq f(x')$,

• $f$ is *surjective* if for every $y \in F$, there exists $x \in E$ such that $y = f(x)$,

• $f$ is *bijjective* if it is both injective and surjective.
**Definition A.0.9 (Continuous map).** Let $E, F$ be two topological spaces. $f: E \to F$ is continuous of $C^0$ if for every open set $U \subset F$, its preimage $f^{-1}(U)$ by $f$ is an open set of $E$.

**Definition A.0.10 (Homeomorphism).** Let $f: E \to F$ be a map between two topological spaces. $f$ is called a homeomorphism if it has the following properties:

* $f$ is a bijection,

* $f$ is continuous,

* the inverse $f^{-1}$ of $f$ is continuous.

**Definition A.0.11 (Differential map).** Let $p, n \in \mathbb{N}$, $f: U \subset \mathbb{R}^n \to \mathbb{R}^p$ a map defined on an open set $U$, and $x_0 \in U$. We say that $f$ is differentiable at $x_0$ if there exists a linear map $L$ defined in $\mathbb{R}^n$ such that

$$ \forall h, \quad f(x_0 + h) = f(x_0) + L(h) + o(\|h\|). $$

In that case, $L$ is unique and is called the differential of $f$ at $x_0$, and written $df_{x_0}$.

**Definition A.0.12 (Class $C^k$).** Let $p, n \in \mathbb{N}$, $f: U \subset \mathbb{R}^n \to \mathbb{R}^p$ a map defined on an open set $U$. $f$ is $C^1$ if it is differentiable on $U$ and the map $x \mapsto df_x$ is continuous on $U$. Similarly we say that $f$ is $C^k$ or of class $C^k$ for $k \in \mathbb{N} \cup \{\infty\}$ if $f$ is $k$-times differentiable.

**Definition A.0.13 (C$^k$-diffeomorphism).** Let $k \in (\mathbb{N} \setminus \{0\}) \cup \{\infty\}$ and let $f: U \to V$ be a map between two open sets of $\mathbb{R}^n$. $f$ is called a *diffeomorphism of class C$^k$* if it has the following properties:

* $f$ is a bijection,

* $f$ is of class $C^k$,

* the inverse $f^{-1}$ of $f$ is $C^k$.

**Definition A.0.14 (Inner product).** Let $E$ be a real vector space. An inner product is a symmetric positive-definite bilinear map $\langle \cdot, \cdot \rangle: E \times E \to \mathbb{R}$, i.e. $\forall x, y \in E, \langle x, y \rangle = \langle y, x \rangle$ and $\langle x, x \rangle > 0 \iff x \neq 0$.

**Definition A.0.15 (σ-algebra).** Let $M$ be a set and $\mathcal{P}(M)$ the set of subsets of $M$. A subset $\Sigma \subseteq \mathcal{P}(M)$ is called σ-algebra if it has the three following properties:

* It is closed under complementation: for any set $S \in \Sigma$, $M \setminus S \in \Sigma$;

* $M$ is an element of $\Sigma$: $M \in \Sigma$;

* $\Sigma$ is closed under countable unions: $\forall (S_i)_{i \in \mathbb{N}}, \bigcup_i S_i \in \Sigma$.

**Definition A.0.16 (Borel σ-algebra).** Let $M$ be a topological set. The Borel σ-algebra $\mathcal{B}(M)$ on $M$ is the smallest σ-algebra that contains all the open sets of $M$.

**Definition A.0.17 (Probability measure).** Let $\mathcal{F}$ be a σ-algebra over a set $\Omega$. A probability measure $\text{Pr}$ is a function $\text{Pr}: \mathcal{F} \to [0, 1]$ such that:

* $\text{Pr}$ is σ-additive: $\forall (S_i)_{i \in \mathbb{N}}, \text{Pr}(\bigcup_i S_i) \in \Sigma_i \text{Pr}(S_i)$.

* $\text{Pr}$ has unit mass: $\text{Pr}(\Omega) = 1$
# B SE(n) with an anisotropic metric

## B.1 Geodesics

Below is the code to plot geodesics of SE(n) with the anisotropic metric described in Example 4.14. The code will be updated on the dedicated Github repository<fnref n="1" />.

```python
import matplotlib.pyplot as plt

import geomstatsisa backend as gs
import geomstatsvisualization as visualization
from geomstats.algebra_utils import from_vector_to_diagonal_matrix
from geomstats_geometry invariant metric import InvariantMetric
from geomstats_geometry.special_euclidean import SpecialEuclidean

SE2_GROUP = SpecialEuclidean(n=2, point_type='matrix')
N_STEPS = 15

def main():
    """Plot geodesics on SE(2) with different structures."""
    theta = gs.pi / 4
    initial_tangent_vec = gs.array([  
        [0., - theta, 1.],  
        [theta, 0., 1.],  
        [0., 0., 0.]])
    t = gs.linspace(0, 1., N_STEPS + 1)
    tangent_vec = gs.einsum('t,ij->tij', t, initial_tangent_vec)

    fig = plt.figure(figsize=(10, 10))
    maxs_x = []
    mins_y = []
    maxs = []
    for i, beta in enumerate([1., 2., 3., 5.]):
        ax = plt.subplots(2, 2, i + 1)
        metric mat = from_vector_to_diagonal_matrixgs.array([1, beta, 1.]))
        metric = InvariantMetric(SE2_GROUP, metric mat, point_type='matrix')
        points = metric .exp(tangent_vec, base point=SE2_GROUP .identity)
        ax = visualization .plot(  
            points, ax=ax, space='SE2_GROUP', color='black',  
            label='r' + '$\beta={}$'.format[beta))
        mins_y.append(min points [:, 1, 2]))
        maxs .append(max points [:, 1, 2]))
        maxs_x.append(max points [:, 0, 2]))
        plt.legend(loc='best')

    for ax in fig .axes:
        x_lim_inf, _ = plt .xlim()
        x_lim = [x_lim_inf, 1.1 * max(maxs_x)]
        y_lim = [min(mins_y) - .1, max(maxs) + .1]
        ax .set_ylim(y_lim)
        ax .set_xlim(x_lim)
        ax .setAspect('equal')
        plt .savefig('../figures/geo-se2.eps', bbox inches='tight', pad inches=0)
        plt .show()

if __name__ == '__main__':
    main()
```

<footnote n="1">¹https://github.com/geomstats/ftmal-paper</footnote>
## B.2 Curvature

This metric is used already in Example 4.13. From the structure constants and Equation (4.4), we can compute the associated Christoffel symbols at identity for the frame $(\tilde{e}_1, \dots, \tilde{e}_6)$. Let $\tau = (\sqrt{\beta} + \frac{1}{\sqrt{\beta}})$. We obtain

$$ \Gamma_{ij}^{k} = \frac{1}{2\sqrt{2}} \text{ if } ijk \text{ is a cycle of } [1,2,3], \qquad (B.1) $$

$$ \Gamma_{15}^{6} = -\Gamma_{16}^{5} = -\frac{2}{\tau}\Gamma_{24}^{6} = \frac{2}{\tau}\Gamma_{26}^{4} = \frac{2}{\tau}\Gamma_{34}^{5} = -\frac{2}{\tau}\Gamma_{35}^{4} = \frac{1}{\sqrt{2}}, \qquad (B.2) $$

and all the others are null.

**Lemma B.2.1.** $(SE(3), g)$ is locally symmetric, i.e. $\nabla R = 0$, if and only if $\beta = 1$.

We can now prove Lemma B.2.1, formulated as: $(SE(3), g)$ is locally symmetric, i.e. $\nabla R = 0$, if and only if $\beta = 1$. This is valid for any dimension $d \ge 2$ provided that the metric matrix $G$ is diagonal, of size $d(d+1)/2$, with ones everywhere except one coefficient of the translation part.

*Proof.* For $\beta = 1$, $(SE(d), g)$ is isometric to $(SO(d) \times \mathbb{R}^d, g_{rot} \oplus g_{trans})$. As the product of two symmetric spaces is again symmetric, $(SE(d), g)$ is symmetric.

We prove the contraposition of the necessary condition. Let $\beta \neq 1$. We give $i, j, k, l$ s.t. $(\nabla_{e_i} R)(e_j, e_k)e_l \neq 0$:

$$ \begin{aligned} (\nabla_{e_3} R)(e_3, e_2)e_4 &= \nabla_{e_3}(R(e_3, e_2)e_4) - R(e_3, \nabla_{e_3} e_2)e_4 - R(e_3, e_2)\nabla_{e_3} e_4 \\ &= \nabla_{e_3}(R(e_3, e_2)e_4) + \frac{1}{\sqrt{2}}R(e_3, e_1)e_4 - \frac{\tau}{2\sqrt{2}}R(e_3, e_2)e_5. \end{aligned} $$
And from the above

$$
\begin{align*}
R(e_3, e_2)e_4 &= \nabla_{e_3} \nabla_{e_2} e_4 - \nabla_{e_2} \nabla_{e_3} e_4 - \nabla_{[e_3, e_2]} e_4 \\
&= -\frac{\tau}{2\sqrt{2}} \nabla_{e_3} e_6 - \nabla_{e_2} e_5 + \frac{1}{\sqrt{2}} \nabla_{e_1} e_4 \\
&= 0.
\end{align*}
$$

$$
\begin{align*}
R(e_3, e_1)e_4 &= \nabla_{e_3} \nabla_{e_1} e_4 - \nabla_{e_1} \nabla_{e_3} e_4 - \nabla_{[e_3, e_1]} e_4 \\
&= -\frac{\tau}{2\sqrt{2}} \nabla_{e_1} e_5 - \frac{1}{\sqrt{2}} \nabla_{e_2} e_4 \\
&= -\frac{\tau}{4} e_6 + \frac{\tau}{4} e_6 = 0.
\end{align*}
$$

$$
\begin{align*}
R(e_3, e_2)e_5 &= \nabla_{e_3} \nabla_{e_2} e_5 - \nabla_{e_2} \nabla_{e_3} e_5 - \nabla_{[e_3, e_2]} e_5 \\
&= \frac{\tau}{2\sqrt{2}} \nabla_{e_2} e_4 + \frac{1}{\sqrt{2}} \nabla_{e_1} e_5 \\
&= -\frac{\tau^2}{8} e_6 + \frac{1}{2} e_6 = \frac{1}{2} \left(1 - \frac{\tau^2}{4}\right) e_6.
\end{align*}
$$

And therefore

$$
\begin{align*}
\beta \neq 1 \implies \tau &= \left(\sqrt{\beta} + \frac{1}{\sqrt{\beta}}\right) \neq 2 \\
&\implies (\nabla_{e_3} R)(e_3, e_1)e_4 = -\frac{\tau}{4\sqrt{2}}\left(1 - \frac{\tau^2}{4}\right)e_6 \neq 0,
\end{align*}
$$

which proves Lemma B.2.1.

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

B.3 One parameter subgroups

Note that these don’t depend on the choice of the metric.

import matplotlib.pyplot as plt
import geomstatsisa as gs
import geomstats visualization as visualization
from geomstats-geometry-special-euclidean import SpecialEuclidean

SE2_GROUP = SpecialEuclidean(n=2, point_type='matrix')
N_STEPS = 30
end_time = 2.7

theta = gs.pi / 3
initial_tangent_vecs = gs.array([ [[0., -theta, 2.], [theta, 0., 2.], [0., 0., 0.]], [[0., -theta, 1.2], [theta, 0., 1.2], [0., 0., 0.]], [[0., -theta, 1.6], [theta, 0., 1.6], [0., 0., 0.]]])
t = gs.linspace(-end_time, end_time, N_STEPS + 1)
```python
fig = plt.figure(figsize=(6, 6))
for tv, col in zip(initial_tangentamo, ['black', 'y', 'g']):
    tangentamo = gs.einsum('t, ij->tij', t, tv)
    groupGeo_points = SE2_GROUP .exp(tangentamo)
    ax = visualization .plot(
        groupGeo_points, space='SE2_GROUP', color=col)
    ax = visualization .plot(
        gs. eye(3) [None, : , :], space='SE2_GROUP', color='slategray')
    ax.setAspect('equal')
    ax .axis("off")
    plt .savefig( '../figures/exponential_se2.eps')
    plt .show()
```
References

Ablin, P., G. Peyré, and T. Moreau (2020). "Super-Efficiency of Automatic Differentiation for Functions Defined as a Minimum". In: *Proceedings of the 37th International Conference on Machine Learning*. Ed. by D. Hal III and A. Singh. Vol. 119. *Proceedings of Machine Learning Research*. PMLR. 32–41. URL: https://proceedings.mlr.press/v119/ablin20a.html.

Absil, P.-A., R. Mahony, and R. Sepulchre (2010). "Optimization On Manifolds: Methods and Applications". In: *Recent Advances in Optimization and Its Applications in Engineering*. Ed. by M. Diehl, F. Glineur, E. Jarlebring, and W. Michiels. Berlin, Heidelberg: Springer. 125–144. DOI: 10.1007/978-3-642-12598-0_12.

Absil, P.-A. and J. Malick (2012). "Projection-like Retractions on Matrix Manifolds". *SIAM Journal on Optimization*. 22(1): 135–158. DOI: 10.1137/100802529.

Afsari, B. (2011). "Riemannian Lp Center of Mass: Existence, Uniqueness, and Convexity". *Proceedings of the American Mathematical Society*. 139(2): 655–673. DOI: 10.1090/S0002-9939-2010-10541-5.

Arnold, V. (1966). "Sur la géométrie différentielle des groupes de Lie de dimension infinie et ses applications à l'hydrodynamique des fluides parfaits". *Annales de l'institut Fourier*. 16(1): 319–361. DOI: 10.5802/aif.233.

Atkin, C. J. (1975). "The Hopf-Rinow Theorem Is False in Infinite Dimensions". *Bulletin of the London Mathematical Society*. 7(3): 261–266. DOI: 10.1112/blms/7.3.261.

Balzin, E. (2021). "Lectures for MAA 306 Course "Topics in Differential Geometry"". June: 102. URL: http://eduard-balzin.perso.math.cnrs.fr/Differential_geometry_lecture_notes_070621.pdf.

Barachant, A. (2015). "pyRiemann v0.2.2". Zenodo. DOI: 10.5281/zenodo.593816.

Barczyk, M., S. Bonnabel, J.-E. Deschaud, and F. Goulette (2015). "Invariant EKF Design for Scan Matching-Aided Localization". *IEEE Transactions on Control Systems Technology*. 23(6): 2440–2448. DOI: 10.1109/TCST.2015.2413933.

Barp, A., A. Kennedy, and M. Girolami (2019). "Hamiltonian Monte Carlo on Symmetric and Homogeneous Spaces via Symplectic Reduction". DOI: 10.48550/arXiv.1903.02699.

Becigneul, G. and O. Ganea (2019). "Riemannian Adaptive Optimization Methods". In: *7th International Conference on Learning Representations (ICLR)*. URL: https://openreview.net/pdf?id=r1eiqi09K7.

Bendokat, T., R. Zimmermann, and P.-A. Absil (2020). "A Grassmann Manifold Handbook: Basic Geometry and Computational Aspects". *Technical Report* No. UCL-INMA-2020.07. arXiv: 2011.13699.

Berger, M. (2003). *A Panoramic View of Riemannian Geometry*. Berlin, Heidelberg: Springer. DOI: 10.1007/978-3-642-18245-7.
Bhatia, R., T. Jain, and Y. Lim (2019). "On the Bures-Wasserstein Distance between Positive Definite Matrices". *Expositiones Mathematicae*. 37(2): 165-191. doi: 10.1016/j.exmath.2018.01.002.

Bhattacharya, R. and V. Patrangenaru (2005). "Large Sample Theory of Intrinsic and Extrinsic Sample Means on Manifolds—II". *The Annals of Statistics*. 33(3): 1225-1259. doi: 10.1214/009053605000000093.

Boumal, N. (2023). *An Introduction to Optimization on Smooth Manifolds*. Cambridge: Cambridge University Press. URL: https://www.cambridge.org/core/books/an-introduction-to-optimization-on-smooth-manifolds/EAF2B35457B7034AC747188DC2FFC058 (accessed on 01/06/2023).

Cartan, E. (1926). "Sur une classe remarquable d'espaces de Riemann". *Bulletin de la Société Mathématique de France*. 54: 214-264. doi: 10.24033/bsmf.1105.

Cartan, E. (1930). *La théorie des groupes finis et continus et l'analysis situs*. Vol. 42. *Mémorial des sciences mathématiques*. Gauthier-Villars et cie. URL: https://books.google.fr/books?id=nv8SAQAAMAAJ.

Cendra, H., D. D. Holm, J. E. Marsden, and T. S. Ratiu (1998). "Lagrangian Reduction, the Euler-Poincaré Equations, and Semidirect Products". In: *Geometry of Differential Equations*. Vol. 186. *American Mathematical Society Translations: Series 2*. American Mathematical Society. 1-25. doi: 10.1090/trans2/186.

Censi, A. (2012). "PyGeometry: Library for Handling Various Differentiable Manifolds." URL: https://github.com/AndreaCensi/geometry.

"Chapter 3 Homogeneous Spaces" (1975). In: *Comparison Theorems in Riemannian Geometry*. Ed. by J. Cheeger and D. G. Ebin. Vol. 9. *North-Holland Mathematical Library*. Elsevier. 55-79. doi: 10.1016/S0924-6509(08)70208-3.

Chevallier, E. and N. Guigui (2020). "A Bi-Invariant Statistical Model Parametrized by Mean and Covariance on Rigid Motions". *Entropy*. 22(4): 432. doi: 10.3390/e22040432.

Chikuse, Y. (2003). *Statistics on Special Manifolds. Lecture Notes in Statistics*. New York: Springer-Verlag. doi: 10.1007/978-0-387-21540-2.

Cohen, T. S., M. Geiger, and M. Weiler (2019). "A General Theory of Equivariant CNNs on Homogeneous Spaces". In: *Advances in Neural Information Processing Systems*. Vol. 32. Curran Associates, Inc. URL: https://proceedings.neurips.cc/paper/2019/hash/b9cfe8b6042cf759dc4c0cccb27a6737-Abstract.html (accessed on 01/06/2023).

David, P. and W. Gu (2019). "A Riemannian Structure for Correlation Matrices". *Operators and Matrices*. 13(3): 607-627. doi: 10.7153/oam-2019-13-46.

Diaconis, P., S. Holmes, and M. Shahshahani (2013). "Sampling from a Manifold". *Advances in Modern Statistical Theory and Applications: A Festschrift in honor of Morris L. Eaton*. Jan.: 102-125. doi: 10.1214/12-IMSCOLL1006.
Dryden, I. L., A. Koloydenko, and D. Zhou (2009). "Non-Euclidean Statistics for Covariance Matrices, with Applications to Diffusion Tensor Imaging". *Annals of Applied Statistics*. 3(3): 1102-1123. DOI: 10.1214/09-AOAS249.

Dryden, I. L. and K. V. Mardia (2016). *Statistical Shape Analysis: With Applications in R*. Wiley Series in Probability and Statistics. John Wiley & Sons. DOI: 10.1002/9781119072492.

Eaton, M. L. (1983). *Multivariate Statistics: A Vector Space Approach*. Probability and Statistics Series. Wiley. URL: https://books.google.fr/books?id=1CvvAAAAMAAJ.

Fey, M. and J. E. Lenssen (2019). "Fast Graph Representation Learning with PyTorch Geometric". Apr. arXiv: 1903.02428.

Fletcher, P., C. Lu, S. Pizer, and S. Joshi (2004). "Principal Geodesic Analysis for the Study of Nonlinear Statistics of Shape". *IEEE Transactions on Medical Imaging*. 23(8): 995-1005. DOI: 10.1109/TMI.2004.831793.

Fréchet, M. (1948). "Les éléments aléatoires de nature quelconque dans un espace distancié". *Annales de l'institut Henri Poincaré*. 10(4): 215-310. URL: http://www.numdam.org/item/AIHP_1948__10__4__215__0/ (accessed on 05/20/2021).

Gallier, J. and J. Quaintance (2020). *Differential Geometry and Lie Groups: A Computational Perspective*. Geometry and Computing. Springer International Publishing. DOI: 10.1007/978-3-030-46040-2.

Gawlik, E. S. and M. Leok (2018). "Interpolation on Symmetric Spaces Via the Generalized Polar Decomposition". *Foundations of Computational Mathematics*. 18(3): 757-788. DOI: 10.1007/s10208-017-9353-0.

Gay-Balmaz, F., D. D. Holm, D. M. Meier, T. S. Ratiu, and F.-X. Vialard (2012). "Invariant Higher-Order Variational Problems II". *Journal of Nonlinear Science*. 22(4): 553-597. DOI: 10.1007/s00332-012-9137-2.

Guigui, N., E. Maignant, A. Trouvé, and X. Pennec (2021). "Parallel Transport on Kendall Shape Spaces". In: *GSI 2021 - 5th Conference on Geometric Science of Information*. Vol. 12829. Lecture Notes in Computer Science. Springer, Cham. 103-110. DOI: 10.1007/978-3-030-80209-7_12.

Guigui, N. and X. Pennec (2021). "A Reduced Parallel Transport Equation on Lie Groups with a Left-Invariant Metric". In: *GSI 2021 - 5th Conference on Geometric Science of Information*. Vol. 12829. Lecture Notes in Computer Science. Springer, Cham. 119-126. DOI: 10.1007/978-3-030-80209-7_14.

Hauberg, S. (2018). "Directional Statistics with the Spherical Normal Distribution". In: *2018 21st International Conference on Information Fusion (FUSION)*. 704-711. DOI: 10.23919/ICIF.2018.8455242.

Hauberg, S. (2019). "Only Bayes Should Learn a Manifold (on the Estimation of Differential Geometric Structure from Data)". Sept. arXiv: 1806.04994.
Helgason, S. (1979). *Differential Geometry, Lie Groups, and Symmetric Spaces*. Vol. 80. *Pure and Applied Mathematics*. Academic Press. URL: https://www.sciencedirect.com/bookseries/pure-and-applied-mathematics/vol/80/suppl/C.

Huckemann, S., T. Hotz, and A. Munk (2010). “Intrinsic Shape Analysis: Geodesic Pca for Riemannian Manifolds Modulo Isometric Lie Group Actions”. *Statistica Sinica*. 20(1): 1–58. URL: https://www.jstor.org/stable/24308976 (accessed on 01/09/2023).

Karcher, H. (1977). “Riemannian Center of Mass and Mollifier Smoothing”. *Communications on Pure and Applied Mathematics*. 30(5): 509–541. DOI: 10.1002/cpa.3160300502.

Karcher, H. (2014). “Riemannian Center of Mass and so Called Karcher Mean”. July. arXiv: 1407.2087.

Kendall, W. S. (1990). “Probability, Convexity, and Harmonic Maps with Small Image I: Uniqueness and Fine Existence”. *Proceedings of the London Mathematical Society*. s3-61(2): 371–406. DOI: 10.1112/plms/s3-61.2.371.

Kendall, W. S. and H. Le (2009). “10 Statistical Shape Theory”. In: *New Perspectives in Stochastic Geometry*. Ed. by W. S. Kendall and I. Molchanov. Oxford University Press. 348–373. DOI: 10.1093/acprof:oso/9780199232574.003.0010.

Kim, K.-R., I. L. Dryden, H. Le, and K. E. Severn (2021). “Smoothing Splines on Riemannian Manifolds, with Applications to 3D Shape Space”. *Journal of the Royal Statistical Society: Series B (Statistical Methodology)*. 83(1): 108–132. DOI: 10.1111/rssb.12402.

Kobayashi, S. and K. Nomizu (1996a). *Foundations of Differential Geometry, Volume 1*. *Wiley Classics Library*. Wiley.

Kobayashi, S. and K. Nomizu (1996b). *Foundations of Differential Geometry, Volume 2*. Wiley.

Kochurov, M., R. Karimov, and S. Kozlukov (2020). “Geoopt: Riemannian Optimization in PyTorch”. DOI: 10.48550/arXiv.2005.02819. arXiv: 2005.02819 [cs].

Kolev, B. (2004). “Lie Groups and Mechanics: An Introduction”. *Journal of Nonlinear Mathematical Physics*. 11(4): 480–498. DOI: 10.2991/jnmp.2004.11.4.5.

Kühnel, L. and S. Sommer (2017). “Computational Anatomy in Theano”. *Graphs in Biomedical Image Analysis, Computational Anatomy and Imaging Genetics*. Sept.: 164–176. DOI: 10.1007/978-3-319-67675-3_15.

Kühnel, L., S. Sommer, and A. Arnaudon (2019). “Differential Geometry and Stochastic Dynamics with Deep Learning Numerics”. *Applied Mathematics and Computation*. 356(Sept.): 411–437. DOI: 10.1016/j.amc.2019.03.044.

Lafontaine, J., S. Gallot, and D. Hulin (2004). *Riemannian Geometry*. Universitext. Springer Verlag. DOI: 10.1007/978-3-642-18855-8.

Le, H. and D. G. Kendall (1993). “The Riemannian Structure of Euclidean Shape Spaces: A Novel Environment for Statistics”. *The Annals of Statistics*. 21(3): 1225–1271. URL: https://www.jstor.org/stable/2242196 (accessed on 11/09/2020).
Lee, J. M. (2003). *Introduction to Smooth Manifolds*. Vol. 218. *Graduate Texts in Mathematics*. New York, NY: Springer Science & Business Media. doi: 10.1007/978-1-4419-9982-5.

Lee, J. M. (2018). *Introduction to Riemannian Manifolds*. Vol. 176. *Graduate Texts in Mathematics*. Cham: Springer International Publishing. doi: 10.1007/978-3-319-91755-9.

Lorenzi, M. and X. Pennec (2013). "Geodesics, Parallel Transport & One-Parameter Subgroups for Diffeomorphic Image Registration". *International Journal of Computer Vision*. 105(2): 111–127. doi: 10.1007/s11263-012-0598-4.

Louis, M., R. Couronné, I. Koval, B. Charlier, and S. Durrleman (2019). "Riemannian Geometry Learning for Disease Progression Modelling". In: *Information Processing in Medical Imaging*. Springer, Cham. 542–553. doi: 10.1007/978-3-030-20351-1_42.

Mardia, K. V. and P. E. Jupp (2009). *Directional Statistics*. *Wiley Series in Probability and Statistics*. John Wiley & Sons. doi: 10.1002/9780470316979.

Marle, C.-M. (2005). "The Works of Charles Ehresmann on Connections: From Cartan Connections to Connections on Fibre Bundles". In: *Geometry and Topology of Manifolds*. Vol. 76. Bedlewo, Poland. 65–86. URL: https://hal.archives-ouvertes.fr/hal-00940427 (accessed on 05/08/2021).

Marron, J. and I. Dryden (2021). *Object Oriented Data Analysis*. Chapman & Hall/CRC Monographs on Statistics & Applied Probability. CRC Press. URL: https://books.google.fr/books?id=X2yjzgEACAAJ.

Marsden, J. E. and T. S. Ratiu (2009). "Mechanical Systems: Symmetries and Reduction". In: *Encyclopedia of Complexity and Systems Science*. Ed. by R. A. Meyers. New York, NY: Springer. 5482–5510. doi: 10.1007/978-0-387-30440-3_326.

McAlpin, J. H. (1965). "Infinite Dimensional Manifolds and Morse Theory". *PhD thesis*. Columbia University.

Meghwanshi, M., P. Jawanpuria, A. Kunchukuttan, H. Kasai, and B. Mishra (2018). "McTorch, a Manifold Optimization Library for Deep Learning". doi: 10.48550/arXiv. 1810.01811.

Milnor, J. (1976). "Curvatures of Left Invariant Metrics on Lie Groups". *Advances in Mathematics*. 21(3): 293–329. doi: 10.1016/S0001-8708(76)80002-3.

Milnor, J. (1984). "Remarks on Infinite-Dimensional Lie Groups". *Relativity, groups and topology*. 2. URL: http://inis.iaea.org/Search/search.aspx?orig_q=RN:16043191 (accessed on 07/06/2021).

Miolane, N. and X. Pennec (2015). "Computing Bi-Invariant Pseudo-Metrics on Lie Groups for Consistent Statistics". *Entropy*. 17(4). doi: 10.3390/e17041850.

Miolane, N., N. Guigui, H. Zaatiti, C. Shewmake, H. Hajri, D. Brooks, A. Le Brigant, J. Mathe, B. Hou, Y. Thanwerdas, S. Heyder, O. Peltre, N. Koep, Y. Cabanes, T. Gerald, P. Chauchat, B. Kainz, C. Donnat, S. Holmes, and X. Pennec (2020). "Introduction to Geometric Learning in Python with Geomstats". In: *SciPy 2020 - 19th Python in Science Conference*. 48. doi: 10.25080/Majora-342d178e-007.
Munthe-Kaas, H. Z., G. R. W. Quispel, and A. Zanna (2014). “Symmetric Spaces and Lie Triple Systems in Numerical Analysis of Differential Equations”. *BIT Numerical Mathematics*. 54(1): 257–282. DOI: 10.1007/s10543-014-0473-5.

Myers, S. B. and N. E. Steenrod (1939). “The Group of Isometries of a Riemannian Manifold”. *Annals of Mathematics*. 40(2): 400–416. DOI: 10.2307/1968928.

Nava-Yazdani, E., H.-C. Hege, T. J. Sulllivan, and C. von Tycowicz (2020). “Geodesic Analysis in Kendall’s Shape Space with Epidemiological Applications”. *Journal of Mathematical Imaging and Vision*. 62(4): 549–559. DOI: 10.1007/s10851-020-00945-w.

Nickel, M. and D. Kiela (2017). “Poincaré Embeddings for Learning Hierarchical Representations”. In: *Advances in Neural Information Processing Systems 30*. Ed. by I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett. Curran Associates, Inc. 6338–6347. URL: http://papers.nips.cc/paper/7213-poincare-embeddings-for-learning-hierarchical-representations.pdf.

Nomizu, K. (1954). “Invariant Affine Connections on Homogeneous Spaces”. *American Journal of Mathematics*. 76(1): 33. DOI: 10.2307/2372398.

O'Neill, B. (1966). “The Fundamental Equations of a Submersion”. *Michigan Mathematical Journal*. 13(4): 459–469. DOI: 10.1307/mmj/1028999604.

Paulin, F. (2007). “Géométrie différentielle élémentaire”. URL: https://www.imo.universite-paris-saclay.fr/~paulin/notescours/cours_geodiff.pdf.

Paulin, F. (2014). “Géométrie Riemannienne, cours de second année de mastère”. URL: https://www.imo.universite-paris-saclay.fr/~paulin/notescours/cours_georiem.pdf.

Pennec, X. (2006). “Intrinsic Statistics on Riemannian Manifolds: Basic Tools for Geometric Measurements”. *Journal of Mathematical Imaging and Vision*. 25(1): 127–154. DOI: 10.1007/s10851-006-6228-4.

Pennec, X. and V. Arsigny (2013). “Exponential Barycenters of the Canonical Cartan Connection and Invariant Means on Lie Groups”. In: *Matrix Information Geometry*. Ed. by F. Nielsen and R. Bhatia. Berlin, Heidelberg: Springer Berlin Heidelberg. 123–166. DOI: 10.1007/978-3-642-30232-9_7.

Pennec, X., P. Fillard, and N. Ayache (2006). “A Riemannian Framework for Tensor Computing”. *International Journal of Computer Vision*. 66(1): 41–66. DOI: 10.1007/s11263-005-3222-z.

Pennec, X., S. Sommer, and T. Fletcher, eds. (2020). *Riemannian Geometric Statistics in Medical Image Analysis*. Vol. 3. The Elsevier and MICCAI Society Book Series. Elsevier. DOI: 10.1016/C2017-0-01561-6.

Postnikov, M. M. (2001). *Geometry VI: Riemannian Geometry*. Encyclopaedia of Mathematical Sciences, Geometry. Berlin Heidelberg: Springer-Verlag. URL: //www.springer.com/us/book/9783540411086 (accessed on 12/21/2018).
Riemann, B. (1868). "Über Die Hypothesen, Welche Der Geometrie Zu Grunde Liegen. (Mitgetheilt Durch r. Dedekind)". *Abhandlungen der Königlichen Gesellschaft der Wissenschaften in Göttingen*. 13: 133–152. URL: http://eudml.org/doc/135760.

Riemann, B. (1873). "On the Hypotheses Which Lie at the Bases of Geometry". *Nature*. 8(184): 36–37. DOI: 10.1038/008036a0.

Said, S., H. Hajri, L. Bombrun, and B. C. Vemuri (2018). "Gaussian Distributions on Riemannian Symmetric Spaces: Statistical Learning With Structured Covariance Matrices". *IEEE Transactions on Information Theory*. 64(2): 752–772. DOI: 10.1109/TIT.2017.2713829.

Schmid, R. (2004). "Infinite Dimensional Lie Groups with Applications to Mathematical Physics". *Journal of Geometry and Symmetry in Physics*. 1(Jan.): 54–120. DOI: 10.7546/jgsp-1-2004-54-120.

Smirnov, O. (2021). "TensorFlow RiemOpt: A Library for Optimization on Riemannian Manifolds". July. arXiv: 2105.13921.

Terras, A. (1988). *Harmonic Analysis on Symmetric Spaces and Applications, Vol. II*. New York: Springer-Verlag. DOI: 10.1007/978-1-4612-3820-1.

Thanwerdas, Y. (2022). "Riemannian and Stratified Geometries on Covariance and Correlation Matrices". *Thesis*. Université Côte d'Azur. URL: https://hal.archives-ouvertes.fr/tel-03698752.

Thanwerdas, Y. and X. Pennec (2021). "Geodesics and Curvature of the Quotient-Affine Metrics on Full-Rank Correlation Matrices". In: *GSI 2021 - 5th Conference on Geometric Science of Information. Proceedings of Geometric Science of Information*. Paris, France. DOI: 10.1007/978-3-030-80209-7_11.

Thanwerdas, Y. and X. Pennec (2022). "O(n)-Invariant Riemannian Metrics on SPD Matrices". *Linear Algebra and its Applications*. Dec. DOI: 10.1016/j.laa.2022.12.009.

Townsend, J., N. Koep, and S. Weichwald (2016). "Pymanopt: A Python Toolbox for Optimization on Manifolds Using Automatic Differentiation". *Journal of Machine Learning Research*. 17(137): 1–5. DOI: 10.5555/2946645.3007090.

Whitney, H. (1936). "Differentiable Manifolds". *Annals of Mathematics*. 37(3): 645–680. DOI: 10.2307/1968482.

Wood, A. T. A. (1994). "Simulation of the von Mises Fisher Distribution". *Communications in Statistics - Simulation and Computation*. 23(1): 157–164. DOI: 10.1080/03610919408813161.

Wynn, K. (2014). "PyQuaternions: A Fully Featured, Pythonic Library for Representing and Using Quaternions". URL: https://github.com/KieranWynn/pyquaternion.

Yair, O., M. Ben-Chen, and R. Talmon (2019). "Parallel Transport on the Cone Manifold of SPD Matrices for Domain Adaptation". In: *IEEE Transactions on Signal Processing*. Vol. 67. 1797–1811. DOI: 10.1109/TSP.2019.2894801.
Zimmermann, R. (2017). "A Matrix-Algebraic Algorithm for the Riemannian Logarithm on the Stiefel Manifold under the Canonical Metric". *SIAM Journal on Matrix Analysis and Applications*. 38(2): 322-342. doi: 10.1137/16M1074485.