How to implement your own manifold

Introduction

This tutorial explains how to implement a manifold using the ManifoldsBase.jl interface. We assume that you are familiar with the basic terminology on Riemannian manifolds, especially the dimension of a manifold, the exponential map, and the inner product on tangent spaces. To read more about this you can for example check ^{[doCarmo1992]}, Chapter 3, first.

Furthermore, we will look at a manifold that is isometrically embedded into a Euclidean space.

In general you need just a data type (struct) that inherits from AbstractManifold to define a manifold. No function is per se required to be implemented. However, it is a good idea to provide functions that might be useful to others, for example check_point and check_vector, as we do in this tutorial.

In this tutorial we will

• model the manifold
• implement two tests, so that points and tangent vectors can be checked for validity, for example also within ValidationManifold,
• implement two functions, the exponential map and the manifold dimension.
• decorate the manifold with an embedding to gain further features.

The next two sections are just a technical detail and the necessary imports to extend the functions defined in this interface.

Technical preliminaries

There are only two small technical things we need to explain at this point before we get started. First of all our AbstractManifold{𝔽} has a parameter 𝔽. This parameter indicates the number_system the manifold is based on, for example ℝ for real manifolds, which is short for RealNumbers() or ℂ for complex manifolds, a shorthand for ComplexNumbers().

Second, this interface usually provides both an allocating and a mutating variant of each function, for example for the exponential map implemented below this interface provides exp(M,p,X) to compute the exponential map and exp!(M, q, p, X) to compute the exponential map in the memory provided by q, mutating that input. the convention is, that the manifold is the first argument – in both function variants – the mutating variant then has the input to be mutated in second place, and the remaining parameters are again the same (pand X here). We usually refer to these two variants of the same function as the allocating (exp) function and the mutating (exp!) one.

The convention for this interface is to document the allocation function, which by default allocates the necessary memory and calls the mutating function. So the convention is to just implement the mutating function, unless there is a good reason to provide an implementation for both. For more details see the design section on mutating and non-mutating functions

Startup

As a start, let's load ManifoldsBase.jl and import the functions we consider throughout this tutorial.

using ManifoldsBase, LinearAlgebra, Test
import ManifoldsBase: check_point, check_vector, manifold_dimension, exp!, inner, representation_size, get_embedding
import Base: show

We load LinearAlgebra for some computations. Test is only loaded for illustrations in the examples.

We import the mutating variant of the exponential map, as just discussed above.

The manifold

The manifold we want to implement here a sphere, with radius $r$. Since this radius is a property inherent to the manifold, it will become a field of the manifolds struct. The second information, we want to store is the dimension of the sphere, for example whether it's the 1-sphere, i.e. the circle, represented by vectors $p\in\mathbb R^2$ of norm $r$ or the 2-sphere in $\mathbb R^3$ of radius $r$. Since the latter might be something we want to dispatch on, we model it as a parameter of the type. In general the struct of a manifold should provide information about the manifold, which are inherent to the manifold or has to be available without a specific point or tangent vector present. This is – most prominently – all information required to determine the manifold dimension.

Note that this a slightly more general manifold than the Sphere in Manifolds.jl

For our example we define the following struct. While a first implementation might also just take AbstractManifold{ℝ} as supertype, we directly take AbstractDecoratorManifold{ℝ}, which will be useful later on. For now it does not make a difference.

"""
    ScaledSphere{N} <: AbstractDecoratorManifold{ℝ}

Define an `N`-sphere of radius `r`. Construct by `ScaledSphere(radius,n)`.
"""
struct ScaledSphere{N} <: AbstractDecoratorManifold{ManifoldsBase.ℝ} where {N}
    radius::Float64
end
ScaledSphere(radius, n) = ScaledSphere{n}(radius)
Base.show(io::IO, M::ScaledSphere{n}) where {n} = print(io, "ScaledSphere($(M.radius),$n)")

Here, the last line just provides a nicer print of a variable of that type. Now we can already initialize our manifold that we will use later, the $2$-sphere of radius $1.5$.

S = ScaledSphere(1.5, 2)

ScaledSphere(1.5,2)

Checking points and tangents

Points on a manifold are usually represented as vector, matrices or more generally arrays. Since we consider vectors of a certain norm (and space dimension), our points are vectors. For an arbitrary vector we would first like to check, that it is a valid point on the manifold. For this one can use the function is_point. This is a function on layer 1 which handles special cases as well cases, so it should not be implemented directly by a user of this interface. The functions that have to be implemented can be found on layer 3. Generically, for both is_point and is_vector, this layer contains a function to check correct size of an array, called check_size For the test of points the function to implement is check_point which we actually will implement, analogously there exists also check_vector. These functions return nothing if the point (vector, size) is a correct/valid and returns an error (but not throw it) otherwise. This is usually a DomainError.

We have to check two things: that a point p is a vector with N+1 entries and its norm is the desired radius.

A first thing we have specify is how points and tangent vectors are represented, that is we have to specify their representation_size

representation_size(::ScaledSphere{N}) where {N} = N+1

This already finishes the size check which check_size performs by default (based on the representation size).

If something has to only hold up to precision, we can pass that down, too using the kwargs..., so all three check_ functions should usually have these in their signature. For our manifold we have to check that the norm of a point p is approximately the specified radius.

function check_point(M::ScaledSphere{N}, p; kwargs...) where {N}
    if !isapprox(norm(p), M.radius; kwargs...)
        return DomainError(norm(p), "The norm of $p is not $(M.radius).")
    end
    return nothing
end

Similarly, we can verify, whether a tangent vector X is valid. Its size is again already checked using check_size, so the only remaining property to verify is, that X is orthogonal to p. We can again use the kwargs.

function check_vector(M::ScaledSphere, p, X; kwargs...)
    if !isapprox(dot(p,X), 0.0; kwargs...)
        return DomainError(dot(p,X), "The tangent $X is not orthogonal to $p.")
    end
    return nothing
end

Note that the function is_vector even can check that the base point of X (the p the tangent space belongs to), can be checked for validity, see its keyword argument check_base_point, so within check_vector this can be (implicitly) assumed to hold.

to test points we can now use

is_point(S, [1.0,0.0,0.0]) # norm 1, so not on S, returns false
@test_throws DomainError is_point(S, [1.5,0.0], true) # only on R^2, throws an error.
p = [1.5,0.0,0.0]
X = [0.0,1.0,0.0]
# The following two tests return true
[ is_point(S, p); is_vector(S,p,X) ]

2-element Vector{Bool}:
 0
 0

Functions on the manifold

For the manifold_dimension we have to just return the N parameter

manifold_dimension(::ScaledSphere{N}) where {N} = N
manifold_dimension(S)

2

Note that we can even omit the variable name in the first line since we do not have to access any field or use the variable otherwise.

To implement the exponential map, we have to implement the formula for great arcs, given a start point p and a direction X on the $n$-sphere of radius $r$ the formula reads

\[\exp_p X = \cos\Bigl(\frac{1}{r}\lVert X \rVert\Bigr)p + \sin\Bigl(\frac{1}{r}\lVert X \rVert\Bigr)\frac{r}{\lVert X \rVert}X.\]

Note that with this choice we for example implicitly assume that the manifold is equipped with that certain metric. This is the default within this interface. To distinguish different metrics, see MetricManifold in Manifolds.jl for more details. Since we here only consider one metric, we do not have to specify that.

An implementation of the mutation version, see the technical note for the naming and reasoning, reads

function exp!(M::ScaledSphere{N}, q, p, X) where {N}
    nX = norm(X)
    if nX == 0
        q .= p
    else
        q .= cos(nX/M.radius)*p + M.radius*sin(nX/M.radius) .* (X./nX)
    end
    return q
end

A first easy check can be done taking p from above and any vector X of length 1.5π from its tangent space. The resulting point is opposite of p, i.e. -p and it is of course a valid point on S.

q = exp(S, p, [0.0,1.5π,0.0])
[isapprox(p, -q); is_point(S, q)]

2-element Vector{Bool}:
 1
 0

Adding an isometric embedding

Since the sphere is isometrically embedded, we do not have to implement the inner(M,p,X,Y) for tangent vectors X, Y in the tangent space at p , but we can “delegate” it to the embedding. The embedding is the Euclidean. The same manifold with a little smaller feature set is available in ManifoldsBase.jl as DefaultManifold for testing purposes.

using ManifoldsBase: DefaultManifold, IsIsometricEmbeddedManifold
import ManifoldsBase: active_traits, get_embedding
using ManifoldsBase: merge_traits

Now we can activate a decorator by specifying that the sphere has the IsIsometricEmbeddedManifold trait for the functions f on our scaled sphere manifold by writing

active_traits(f, ::ScaledSphere, args...) = merge_traits(IsIsometricEmbeddedManifold())

and then specifying that said embedding is the DefaultManifold.

get_embedding(::ScaledSphere{N}) where {N} = DefaultManifold(N+1)

Now metric related functions are passed to this embedding, for example the inner product. It now works by using the inner product from the embedding, so we can compute the inner product by calling inner

X = [0.0, 0.1, 3.0]
Y = [0.0, 4.0, 0.2]
# returns 1.0 by calling the inner product in DefaultManifold(3)
inner(S, p, X, Y)

1.0

Conclusion

You can now just continue implementing further functions from ManifoldsBase.jl but with just exp! you for example already have

• geodesic the (not necessarily shortest) geodesic emanating from p in direction X.
• the ExponentialRetraction, that the retract function uses by default.

For the shortest_geodesic the implementation of a logarithm log, or just log! is sufficient.

Sometimes a default implementation is provided; for example if you implemented inner, the norm is defined. You should overwrite it, if you can provide a more efficient version. For a start the default should suffice. With log! and inner you get the distance, and so.

In summary with just these few functions you can already explore the first things on your own manifold. Whenever a function from Manifolds.jl requires another function to be specifically implemented, you get a reasonable error message.

Literature