# 🚀 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 📖.

## 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
This is how the Error we expect here looks like
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
This is how the Error we expect here looks like
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.

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, pts, 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 AbstractRetractionMethodm, the same for the 🔗 inverse_retract(M, p, q, n) with an AbstractInverseRetractionMethodn.

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, pts, n)
return retract(M, pts, 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.

## 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

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.

AbsilMahonySepulchre2008

Absil, P.-A., Mahony, R. and Sepulchre R., Optimization Algorithms on Matrix Manifolds Princeton 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