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

s 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 `exp`

onential 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 (`p`

and `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 `exp`

onential 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 `exp`

onential 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

- doCarmo1992
do Carmo, Manfredo

**Riemannian Geometry**, Birkhäuser Boston, 1992, ISBN: 0-8176-3490-8.