How to Implement a Manifold
This tutorial illustrates, how to implement your very first manifold. We start from the very beginning and cover the basic ideas of the interface provided by ManifoldsBase.jl
interface.
Preliminaries
We will use a simple example in this tutorial, since the main focus here is to illustrate how to define a manifold. We will use the sphere of radius $r$ embedded in $\mathbb R^{d+1}$, i.e. all vectors of length $r$. Formally we define
\[\mathbb S_r^d := \bigl\{ p \in \mathbb R^{d+1} \big| \lVert p \rVert = r \bigr\}\]
When defining a Riemannian manifold mathematically, there is several things to keep in mind, for example the metric imposed on the tangent spaces. For this interface we assume these things to be given implicitly for a first implementation, but they can be made more precise when necessary.
The only thing we have to be aware of for now is the number_system
, i.e. whether our manifold is a real-valued or a complex-valued manifold. The abstract type all manifolds inherit from, the AbstractManifold
{𝔽}
has this number system as a parameter. The usual parameter we will use are the RealNumbers
()
, which have a short hand in ManifoldsBase.jl
, namely ℝ
. The second one are the ComplexNumbers
()
, or ℂ
for short.
using LinearAlgebra, ManifoldsBase
using ManifoldsBase: ℝ
Defining a manifold
A manifold itself is a struct
that is a subtype of AbstractManifold
and should contain. We usually recommend to also document your new manifold.
Since the Sphere
is already a name used within Manifolds.jl
, let’s use a slightly more specific name. We define
"""
ScaledSphere <: AbstractManifold{ℝ}
Define a sphere of fixed radius
# Fields
* `dimension` dimension of the sphere
* `radius` the radius of the sphere
# Constructor
ScaledSphere(dimension,radius=1.0)
Initialize the manifold to a certain `dimension` and `radius`,
which by default is set to `1.0`
"""
struct ScaledSphere <: AbstractManifold{ℝ}
dimension::Int
radius::Float64
end
And we can directly use this manifold and set
M = ScaledSphere(2,1.5)
ScaledSphere(2, 1.5)
Functions I: Manifold properties
While the interface provides a lot of possible functions to define for your manifold, you only need to define those that are necessary for your implementation. If you are using other packages depending on ManifoldsBase.jl
you will often just get a “Method not defined” and sometimes an ambiguity error indicating that a function is missing that is required for a certain task.
We can first start with a technical function which internally ist often used. Any of our points or tangent vectors is represented as a $(d+1)$-dimensional vector. This is internally often used when allocating memory, see representation_size
. It returns a tuple representing the size of arrays for valid points. We define
import ManifoldsBase: representation_size
representation_size(M::ScaledSphere) = (M.dimension+1,)
Similarly, we can implement the function returning the dimension of the manifold, cf. manifold_dimension
as
import ManifoldsBase: manifold_dimension
manifold_dimension(M::ScaledSphere) = M.dimension
and we can now easily use them to access the dimension of the manifold
manifold_dimension(M)
2
Functions II: Verifying Points and tangent vectors
The first two functions we want to define are those to check points and tangent vectors for our manifold. Let’s first clarify what the tangent space looks like. The directions “we can walk into” from a point $p\in \mathbb S_r^d$ are all $X$ that are orthogonal to $p$, which is the plane/vector space tangent to the sphere. Formally
\[T_p\mathbb S_r^d := \bigl\{ X \in \mathbb R^{d+1} \big| \langle p, X \rangle = 0 \bigr\}, \qquad p \in \mathbb S_r^d\]
to verify either p
or X
one uses is_point
(M,p)
and is_vector
(M, p, X)
respectively. Since both involve some automatic options and possibililities, for example whether to throw an error or just return false, both mention that the actual functions to implement are check_point
and check_vector
, which both do not throw but return an error if something is wrong.
In principle we would have to check two properties, namely that the size of p
is of correct size M.dimension+1
and that its norm is M.radius
. Luckily, by defining representation_size
the first check is automatically done already – actually for both points and vectors. We define
import ManifoldsBase: check_point
function check_point(M::ScaledSphere, p; kwargs...)
if !isapprox(norm(p), M.radius; kwargs...)
return DomainError(norm(p), "The norm of $p is not $(M.radius).")
end
return nothing
end
And we can directly test the function. To see all 3 failing ones, we switch from errors to warnings in the check
is_point(M, [1.5, 0.0], error=:warn) # wrong size
is_point(M, [1.0, 0.0, 0.0], error=:warn) # wrong norm
is_point(M, [1.5, 0.0, 0.0], error=:warn) # on the manifold, returns true
┌ Warning: DomainError with (2,):
│ The point [1.5, 0.0] can not belong to the manifold ScaledSphere(2, 1.5), since its size (2,) is not equal to the manifolds representation size ((3,)).
└ @ ManifoldsBase ~/work/ManifoldsBase.jl/ManifoldsBase.jl/src/ManifoldsBase.jl:799
┌ Warning: DomainError with 1.0:
│ The norm of [1.0, 0.0, 0.0] is not 1.5.
└ @ ManifoldsBase ~/work/ManifoldsBase.jl/ManifoldsBase.jl/src/ManifoldsBase.jl:812
true
similarly for vectors, we just have to implement the orthogonality check.
import ManifoldsBase: check_vector
function check_vector(M::ScaledSphere, p, X; kwargs...)
if !isapprox(dot(p,X), 0.0; kwargs...)
return DomainError(
dot(p,X),
"The tangent vector $X is not orthogonal to $p."
)
end
return nothing
end
and again, the high level interface can be used to display warning for vectors not fulfilling the criterion, but we can also activate a check for the point using the last positional argument
is_vector(M, [1.5, 0.0, 0.0], [0.0, 1.0]; error=:warn) # wrong size
is_vector(M, [1.5, 0.0, 0.0], [1.0, 1.0, 0.0]; error=:warn) # not orthogonal norm
is_vector(M, [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], true; error=:warn) # point not valid
is_vector(M, [1.5, 0.0, 0.0], [0.0, 1.0, 0.0], true; error=:warn) # returns true
┌ Warning: DomainError with (2,):
│ The tangent vector [0.0, 1.0] can not belong to the manifold ScaledSphere(2, 1.5), since its size (2,) is not equal to the manifodls representation size ((3,)).
└ @ ManifoldsBase ~/work/ManifoldsBase.jl/ManifoldsBase.jl/src/ManifoldsBase.jl:883
┌ Warning: DomainError with 1.5:
│ The tangent vector [1.0, 1.0, 0.0] is not orthogonal to [1.5, 0.0, 0.0].
└ @ ManifoldsBase ~/work/ManifoldsBase.jl/ManifoldsBase.jl/src/ManifoldsBase.jl:896
┌ Warning: DomainError with 1.0:
│ The norm of [1.0, 0.0, 0.0] is not 1.5.
└ @ ManifoldsBase ~/work/ManifoldsBase.jl/ManifoldsBase.jl/src/ManifoldsBase.jl:812
true
Functions on Manifolds III: The exponential map and a retraction.
For the final group of functions, we want to implement the exp
onential map and a retract
ion. Both are ways to “move around” on the manifold, that is, given a point $p$ and a tangent vector indicating a “walking direction”, the two functions we define will return a new point $q$ on the manifold.
For functions that compute a new point or tangent vector, ManifoldsBase.jl
always provides two varinats: One that allocates new memory and one, that allows to provide the memory, the result should be returned in – to spare memory allocations.
Let’s first take a look at what the exponential map is defined like. We follow the shortest curves, that is great arcs, on the sphere. Formally we have
\[\exp_p X = \cos\Bigl(\frac{1}{r}\lVert X \rVert\Bigr)p + \sin\Bigl(\frac{1}{r}\lVert X \rVert\Bigr)\frac{X}{\lVert X \rVert}.\]
In fact, from the two functions above, exp
(M, p, X)
and exp!
(M, q, p, X)
, that works in place of q
, we only have to implement the second. The first one, exp
by default falls back to allocating memory and calling the secnod. Sp exp
should only be defined, if there is a special reason for. Furthermore, we usually do not verify/check inputs to spare time. If a user feels insecure, they could for example use the ValidationManifold
wrapper which adds explicitly checks of inputs and outputs.
We define
import ManifoldsBase: exp!
function exp!(M::ScaledSphere, q, p, X)
nX = norm(X)
if nX == 0
q .= p
else
q .= cos(nX/M.radius)*p + M.radius*sin(nX/M.radius) .* (1/nX) .* X
end
return q
end
and we can directly test our function starting in the north pole and “waling down” to the equator
q = exp(M, [0.0, 0.0, 1.5], [0.75π, 0.0, 0.0])
3-element Vector{Float64}:
1.5
0.0
9.184850993605148e-17
but we also get the other variants for free, for example
q2 = zero(q)
exp!(M, q2, [0.0, 0.0, 1.5], [0.75π, 0.0, 0.0])
q2
3-element Vector{Float64}:
1.5
0.0
9.184850993605148e-17
or the one that shortens or enlongates the path by a factor, for example, if we walk twice the distance, we reach the opposite point
exp!(M, q2, [0.0, 0.0, 1.5], [0.75π, 0.0, 0.0], 2.0)
q2
3-element Vector{Float64}:
1.8369701987210297e-16
0.0
-1.5
Of course we can easliy check that both points we computed still lie on the sphere
all([is_point(M, q), is_point(M, q2)])
true
Since the exponential map might in general be expensive, we can do a similar implementation with the ProjectionRetraction
. Here, we really have to take into account, that the interface is designed with 3 levels in mind: While the actual function we would call in the end is retract(M, p, X, ProjectionRetraction())
(or its !
variant), we actually have to implement retract_project!(M, q, p, X, t)
for technical details, that are a bit beyond this introductionary tutorial. In short this split avoids ambiguity errors for decorators of the manifolds. We define
import ManifoldsBase: retract_project!
function retract_project!(M::ScaledSphere, q, p, X, t)
q .= (p + t*X) .* (M.radius/norm(p + t*X))
return q
end
retract_project! (generic function with 2 methods)
And to test also this function, and avoiding allocations at the same time, we call
retract!(M, q, [0.0, 0.0, 1.5], [0.75π, 0.0, 0.0], ProjectionRetraction())
3-element Vector{Float64}:
1.2653454121031529
0.0
0.8055439082194726
Finally, there is default_retraction_method
to specify which is the default retraction to use. By default this is
default_retraction_method(M)
ExponentialRetraction()
But we can easily specify this for our manifold as well, for example defining
import ManifoldsBase: default_retraction_method
default_retraction_method(::ScaledSphere) = ProjectionRetraction()
default_retraction_method (generic function with 6 methods)
Then
default_retraction_method(M)
ProjectionRetraction()
and retract without a method specified would always fall back to using the projection retraction instead of the exponential map. Note that for compatibilty there is the AbstractRetractionMethod
called ExponentialRetraction
which makes retract
fall back to calling exp
.
Technical Details
This notebook was rendered with the following environment
Pkg.status()
Status `~/work/ManifoldsBase.jl/ManifoldsBase.jl/tutorials/Project.toml`
[7073ff75] IJulia v1.25.0
[3362f125] ManifoldsBase v0.15.17 `~/work/ManifoldsBase.jl/ManifoldsBase.jl`
[91a5bcdd] Plots v1.40.8