# ๐ Get Started with `Manifolds.jl`

This is a short overview of `Manifolds.jl`

.

This tutorial is rendered from a pluto notebook, so you can also open the file ๐ tutorials/getstarted.jl in Pluto and work on this tutorial interactively.

As usual, if you want to install the package, just type

`] add Manifolds`

in Julia REPL or use

`using Pkg; Pkg.add("Manifolds");`

before the first use. Then load the package with

`using Manifolds, Random`

`Random.seed!(42);`

Since the package heavily depends on `ManifoldsBase.jl`

we will sometimes also link to the interface definition of functions in the interface and mark this with ๐. When referring to Wikipedia, the link is marked with ๐.

## Contents

Representations with and without charts (see notebook

`working-in-charts.jl`

, currently not rendered as a part of this documentation).

## Using the Library of Manifolds

`Manifolds.jl`

is first of all a library of manifolds, see the list in the menu here under โbasic manifoldsโ.

Let's look at three examples together with the first few functions on manifolds.

#### 1. The Euclidean space

The Euclidean space brings us (back) into linear case of vectors, so in terms of manifolds, this is a very simple one. It is often useful to compare to classical algorithms, or implementations.

`Mโ = Euclidean(3)`

Euclidean(3; field = โ)

Note that since a manifold is a type in Julia, we write it in CamelCase. Its parameters are first a dimension or size parameter of the manifold, sometimes optional is a field the manifold is defined over.

For example the above definition is the same as the real-valued case

`Mโ === Euclidean(3, field=โ)`

true

But we even introduced a short hand notation, since โ is also just a symbol/variable to use

`Mโ === โ^3`

true

And similarly here are two ways to create the manifold of vectors of length two with complex entries โ or mathematically the space $\mathbb C^2$

true

The easiest to check is the dimension of a manifold. Here we have three โdirections to walk intoโ at every point $p\in \mathbb R ^3$ so ๐ `manifold_dimension`

) is

`manifold_dimension(Mโ)`

3

#### 2. The hyperpolic space

The $d$-dimensional hyperbolic space is usually represented in $\mathbb R^{d+1}$ as the set of points $p\in\mathbb R^3$ fulfilling

$$p_1^2+p_2^2+\cdots+p_d^2-p_{d+1}^2 = -1.$$

`Mโ = Hyperbolic(2)`

Hyperbolic(2)

`manifold_dimension(Mโ)`

2

Here, a useful function is to check, whether some $pโ\mathbb R^3$ is a point on the manifold. We can check

`is_point(Mโ, [0, 0, 1])`

true

`is_point(Mโ, [1, 0, 1])`

false

Keyword arguments are passed on to any numerical checks, for example an absolute tolerance when checking the above equiality.

`is_point(Mโ, [0, 0, 1.001]; atol=1e-3)`

true

But in an interactive session an error message might be helpful. A positional (third) argument is present to activate this. Here we illustrate this with try-catch to keep the notebook as valid running code.

`@expect_error is_point(Mโ, [0, 0, 1.001], true) DomainError`

```
DomainError with -1.0020009999999997:
The point [0.0, 0.0, 1.001] does not lie on Hyperbolic(2) since its Minkowski inner product is not -1.
```

#### 3. The sphere

The sphere$\mathbb S^d$ is the $d$-dimensional sphere represented in its embedded form, that is unit vectors $p \in \mathbb R^{d+1}$ with unit norm $\lVert p \rVert_2 = 1$.

`Mโ = Sphere(2)`

Sphere(2, โ)

Here we can show a last nice check: `is_vector`

to check whether a tangent vector `X`

is a representation of a tangent vector $XโT_p\mathcal M$ to a point `p`

on the manifold.

This function has two positional asrguments, the first to again indicate whether to throw an error, the second to disable the check that `p`

is a valid point on the manifold. Usually this validity is essential for the tangent check, but if it was for example performed before, it can be turned off to spare time.

For example in our first example the point is not of unit norm

`is_vector(Mโ, [2, 0, 0], [0, 1, 1])`

false

But the orthogonality of `p`

and `X`

is still valid, so we get

`is_vector(Mโ, [2, 0, 0], [0, 1, 1], true, false)`

true

But of course it is better to use a valid point in the first place

`is_vector(Mโ, [1, 0, 0], [0, 1, 1])`

true

`@expect_error is_vector(Mโ, [1, 0, 0], [0.1, 1, 1], true) DomainError`

```
DomainError with 0.1:
The vector [0.1, 1.0, 1.0] is not a tangent vector to [1, 0, 0] on Sphere(2, โ), since it is not orthogonal in the embedding.
```

To learn about how to define a manifold youself check out the How to define your own manifold tutorial of `ManifoldsBase.jl`

.

### Building more advanced manifolds

Based on these basic manifolds we can directly build more advanced manifolds.

The first one concerns vectors or matrices of data on a manifold, the PowerManifold.

`Mโ = Mโ^2`

PowerManifold(Hyperbolic(2), 2)

Then points are represented by arrays, where the power manifold dimension is added in the end. In other words โ for the hyperbolic manifold here, we have a matrix with 2 columns, where each column is a valid point on hyperbolic space.

`p = [0 0; 0 1; 1 sqrt(2)]`

3ร2 Matrix{Float64}: 0.0 0.0 0.0 1.0 1.0 1.41421

`[is_point(Mโ, p[:, 1]), is_point(Mโ, p[:, 2])]`

2-element Vector{Bool}: 1 1

But of course the method we used previously also works for power manifolds:

`is_point(Mโ, p)`

true

Note that nested power manifolds are combined into one as in

`Mโโ = Mโ^4`

PowerManifold(Hyperbolic(2), 2, 4)

which represents $2\times 4$ โ matrices of hyperbolic points represented in $3\times 2\times 4$ arrays.

We can โ alternatively โ use a power manifold with nested arrays

`Mโ
= PowerManifold(Mโ, NestedPowerRepresentation(), 2)`

PowerManifold(Sphere(2, โ), ManifoldsBase.NestedPowerRepresentation(), 2)

which emphasizes that we have vectors of length 2 that contain points, so we store them that way.

`pโ = [[0.0, 0.0, 1.0], [0.0, 1.0, 0.0]]`

2-element Vector{Vector{Float64}}: [0.0, 0.0, 1.0] [0.0, 1.0, 0.0]

Top unify both representations, elements of the power manifold can also be accessed in the classical indexing fashion, if we start with the corresponding manifold first. This way one can implement algorithms also independent of which representation is used.

`p[Mโ, 1]`

3-element Vector{Float64}: 0.0 0.0 1.0

`pโ[Mโ
, 2]`

3-element Vector{Float64}: 0.0 1.0 0.0

Another construtor is the ProductManifold to combine different manifolds. Here of course the order matters. First we construct these using $ร$

`Mโ = Mโ ร Mโ`

ProductManifold with 2 submanifolds: Hyperbolic(2) Sphere(2, โ)

Since now the representations might differ from element to element, we have to encapsulate these in their own type.

`pโ = Manifolds.ArrayPartition([0, 0, 1], [0, 1, 0])`

([0, 0, 1], [0, 1, 0])

Here `ArrayPartition`

taken from `RecursiveArrayTools.jl`

to store the point on the product manifold efficiently in one array, still allowing efficient access to the product elements.

`is_point(Mโ, pโ, true)`

true

But accessing single components still works the same.

`pโ[Mโ, 1]`

3-element Vector{Int64}: 0 0 1

Finally, also the TangentBundle, the manifold collecting all tangent spaces on a manifold is available as

TangentBundle(Sphere(2, โ))

## Implementing generic Functions

In this section we take a look how to implement generic functions on manifolds.

For our example here, we want to implement the so-called ๐ Bรฉzier curve using the so-called ๐ de-Casteljau algorithm. The linked algorithm can easily be generalised to manifolds by replacing lines with geodesics. This was for example used in [BergmannGousenbourger2018] and the following example is an extended version of an example from [AxenBaranBergmannRzecki2022].

The algorithm works recursively. For the case that we have a Bรฉzier curve with just two points, the algorithm just evaluates the geodesic connecting both at some time point $tโ[0,1]$. The function to evaluate a shortest geodesic (it might not be unique, but then a deterministic choice is taken) between two points `p`

and `q`

on a manifold `M`

๐ `shortest_geodesic(M, p, q, t)`

.

```
function de_Casteljau(M::AbstractManifold, t, pts::NTuple{2})
return shortest_geodesic(M, pts[1], pts[2], t)
end
```

de_Casteljau (generic function with 1 method)

```
function de_Casteljau(M::AbstractManifold, t, pts::NTuple)
p = de_Casteljau(M, t, pts[1:(end - 1)])
q = de_Casteljau(M, t, pts[2:end])
return shortest_geodesic(M, p, q, t)
end
```

de_Casteljau (generic function with 2 methods)

This works fine on the sphere, see this tutorial for an optimization task involving Bรฉzier curves.

Now on several manifolds the ๐ exponential map and its (locally defined) inverse, the logarithmic map might not be available in an implementation. So one way to generalise this, is the use of a retraction (see [AbsilMahonySepulchre2008], Def. 4.1.1 for details) and its (local) inverse.

The function itself is quite similar to the expponential map, just that ๐ `retract(M, p, X, m)`

has one further parameter, the type of retraction to take, so `m`

is a subtype of `AbstractRetractionMethod`

`m`

, the same for the ๐ `inverse_retract(M, p, q, n)`

with an `AbstractInverseRetractionMethod`

`n`

.

Now, thinking of a generic implementation, we would like to have a way to specify one, that is available. This can be done by using ๐ `default_retraction_method`

and ๐ `default_inverse_retraction_method`

, respectively. We implement

```
function generic_de_Casteljau(
M::AbstractManifold,
t,
pts::NTuple{2};
m::AbstractRetractionMethod=default_retraction_method(M),
n::AbstractInverseRetractionMethod=default_inverse_retraction_method(M),
)
X = inverse_retract(M, pts[1], pts[2], n)
return retract(M, pts[1], X, t, m)
end
```

generic_de_Casteljau (generic function with 1 method)

and for the recursion

```
function generic_de_Casteljau(
M::AbstractManifold,
t,
pts::NTuple;
m::AbstractRetractionMethod=default_retraction_method(M),
n::AbstractInverseRetractionMethod=default_inverse_retraction_method(M),
)
p = generic_de_Casteljau(M, t, pts[1:(end - 1)]; m=m, n=n)
q = generic_de_Casteljau(M, t, pts[2:end]; m=m, n=n)
X = inverse_retract(M, p, q, n)
return retract(M, p, X, t, m)
end
```

generic_de_Casteljau (generic function with 2 methods)

Note that on a manifold `M`

where the exponential map is implemented, the `default_retraction_method(M)`

returns ๐ `ExponentialRetraction`

, which yields that the `retract`

function falls back to calling `exp`

.

The same mechanism exists for ๐ `parallel_transport_to(M, p, X, q)`

and the more general ๐ `vector_transport_to(M, p, X, q, m)`

whose ๐ `AbstractVectorTransportMethod`

`m`

has a default defined by ๐ `default_vector_transport_method(M)`

.

## Allocating and in-place computations

Memory allocation is a critical performace issue when programming in Julia. To take this into account, `Manifolds.jl`

provides special functions to reduce the amount of allocations.

We again look at the ๐ exponential map. On a manifold `M`

the exponential map needs a point `p`

(to start from) and a tangent vector `X`

, which can be seen as direction to โwalk intoโ as well as the length to walk into this direction. In `Manifolds.jl`

the function can then be called with `q = exp(M, p, X)`

๐. This function returns the resulting point `q`

, which requires to allocate new memory.

To avoid this allocation, the function `exp!(M, q, p, X)`

๐ can be called. Here `q`

is allocated beforehand and is passed as the memory, where the result is returned in. It might be used even for interims computations, as long as it does not introduce side effects. Thas means that even with `exp!(M, p, p, X)`

the result is correct.

Let's look at an example.

We take another look at the `Sphere`

but now a high-dimensional one. We can also illustrate how to generate radnom points and tangent vectors.

`M = Sphere(10000)`

Sphere(10000, โ)

`pโ = rand(M);`

`X = rand(M; vector_at=pโ);`

Now looking at the allocations required we get

`@allocated exp(M, pโ, X)`

80112

While if we have the memory

`qโ = zero(pโ);`

There are no new memory allocations necessary if we use the in-place function.

`@allocated exp!(M, qโ, pโ, X)`

0

This methodology is used for all functions that compute a new point or tangent vector. By default all allocating functions allocate memory and call the in-place function. This also means that if you implement a new manifold, you just have to implement the in-place version.

## Decorating a manifold

As you saw until now, an ๐ `AbstractManifold`

describes a Riemannian manifold. For completeness, this also includes the chosen Riemannian metric tensor or inner product on the tangent spaces.

In `Manifolds.jl`

these are assumed to be a โreasonable defaultโ. For example on the `Sphere(n)`

we used above, the default metric is the one inherited from restricting the inner product from the embedding space onto each tangent space.

Consider a manifold like

`Mโ = SymmetricPositiveDefinite(3)`

SymmetricPositiveDefinite(3)

which is the manifold of $3ร3$ matrices that are symmetric and positive definite. which has a default as well, the affine invariant `LinearAffineMetric`

, but also has several different metrics.

To switch the metric, we use the idea of a ๐ decorator pattern-like approach. Defining

`Mโโ = MetricManifold(Mโ, BuresWassersteinMetric())`

MetricManifold(SymmetricPositiveDefinite(3), Manifolds.BuresWassersteinMetric())

changes the manifold to use the `BuresWassersteinMetric`

.

This changes all functions that depend on the metric, most prominently the Riemannian matric, but also the exponential and logarithmic map and hence also geodesics.

All functions that are not dependent on a metric โ for example the manifold dimension, the tests of points and vectors we already looked at, but also all retractions โ stay unchanged. This means that for example

`[manifold_dimension(Mโโ), manifold_dimension(Mโ)]`

2-element Vector{Int64}: 6 6

both calls the same underlying function. On the other hand with

`pโ
, Xโ
= one(zeros(3, 3)), [1.0 0.0 1.0; 0.0 1.0 0.0; 1.0 0.0 1.0]`

([1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0], [1.0 0.0 1.0; 0.0 1.0 0.0; 1.0 0.0 1.0])

but for example the exponential map and the norm yield different results

`[exp(Mโ, pโ
, Xโ
), exp(Mโโ, pโ
, Xโ
)]`

2-element Vector{Matrix{Float64}}: [4.194528049465325 0.0 3.194528049465325; 0.0 2.718281828459045 0.0; 3.194528049465325 0.0 4.194528049465328] [2.5 0.0 1.5; 0.0 2.25 0.0; 1.5 0.0 2.5]

`[norm(Mโ, pโ
, Xโ
), norm(Mโโ, pโ
, Xโ
)]`

2-element Vector{Float64}: 2.23606797749979 1.118033988749895

Technically this done using Traits โ the trait here is the `IsMetricManifold`

trait. Our trait system allows to combine traits but also to inherit properties in a hierarchical way, see ๐ here for the technical details.

The same approach is used for

specifying a different connection

specifying a manifold as a certain quotient manifold

specifying a certain embedding

specify a certain group action

Again, for all of these, the concrete types only have to be used if you want to do a second, different from the details, property, for example a second way to embed a manfiold. If a manifold is (in its usual representation) an embedded manifold, this works with the default manifold type already, since then it is again set as the reasonable default.

## References

AbsilMahonySepulchre2008

Absil, P.-A., Mahony, R. and Sepulchre R.,

Optimization Algorithms on Matrix ManifoldsPrinceton University Press, 2008, doi: 10.1515/9781400830244open access

AxenBaranBergmannRzecki2022

Axen, S. D., Baran, M., Bergmann, R. and Rzecki, K:

Manifolds.jl: An Extensible Julia Framework for Data Analysis on Manifolds, arXiv preprint, 2022, 2106.08777

BergmannGousenbourger2018

Bergmann, R. and Gousenbourger, P.-Y.:

A variational model for data fitting on manifolds by minimizing the acceleration of a Bรฉzier curve. Frontiers in Applied Mathematics and Statistics, 2018. doi: 10.3389/fams.2018.00059, arXiv: 1807.10090