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