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:808
┌ 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:821

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:892
┌ 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:905
┌ 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:821

true

Functions on Manifolds III: The exponential map and a retraction.

For the final group of functions, we want to implement the exponential map and a retraction. 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.26.0
  [3362f125] ManifoldsBase v0.15.21 `~/work/ManifoldsBase.jl/ManifoldsBase.jl`
  [91a5bcdd] Plots v1.40.9