Notebook source code:
notebooks/01_foundations__manifolds.ipynb

Run it yourself on binder

\(\textbf{Lead Author: Adele Myers}\)

Inspired by: Guigui, Miolane, Pennec, 2022. Introduction to Riemannian Geometry and Geometric Statistics: from basic theory to implementation with Geomstats.

Note: Before you start to read this notebook, run the following code. This will import packages that will allow later code to run.

```
In [1]:
```

```
import geomstats.backend as gs
gs.random.seed(2020)
```

# 1. Introduction#

Welcome to the `Manifold`

class notebook! In this notebook, we will discuss manifolds and the `Manifold`

class in Geomstats. Geomstats is a software which analyzes data on manifolds, so naturally the `Manifold`

class is quite important.

You will learn:

What is a manifold?

Why do we care about manifolds?

What does the

`Manifold`

class do?How is the

`Manifold`

class structured?What is an open set?

What is a closed set?

What is a tangent space?

# 2. What is a Manifold?#

\(\textbf{Intuition:}\) When you are first learning, it can be a helpful starting point to think of a manifold as a surface. This surface can have any dimension and any shape as long as it is smooth (in the sense of being continuous and differentiable). (This is not a particularly precise definition, but it can be helpful starting point for building intuition.) For example, a hypersphere is a two dimensional manifold, and we will often use this manifold in examples.

## 2.1 Precise Mathematical Definition#

A manifold is a set of points that satisfy a specific set of constraints. More specifically, a nonempty subset M \(\subseteq \mathbb{R}^{N}\) is a d-dimensional manifold if and only if ANY of the following conditions hold:

(Local parametrization) For every \(p \in M\), there are two open subsets \(V \subseteq \mathbb{R}^{d}\) and \(U \subseteq \mathbb{R}^{N}\) with \(p \in U\) and \(0 \in V\). There is also a smooth function \(f: V \to \mathbb{R}^{N}\) such that \(f(0) = p\), where \(f\) is a homeomorphism between V and \(U \cap M\), and \(f\) is an immersion at 0.

(Local implicit function) For every \(p \in M\), there exists an open set \(U \in \mathbb{R}^{N}\) and a smooth map \(f: U \to \mathbb{R}^{N-d}\) that is a submersion at p, such that \(U \cap M = f^{-1}\)({0}).

(Local Graph) For every \(x \in M\), there exists an open neighborhood \(U \subseteq \mathbb{R}^{N}\) of \(x\), a neighborhood \(V \subseteq \mathbb{R}^{d}\) of 0 and a smooth map \(f: V \to \mathbb{R}^{N-d}\) such that \(U \cap M = graph(f)^{2}\)

## 2.2 Imprecise “Layman’s Terms” Definition:#

A \(\textbf{manifold}\) is a set of points that satisfy a specific set of constraints, and these constraints are given by the mathematical definition above. That means that if a set of points satisfies any one of the conditions, then it is a manifold. However, these conditions can be hard to digest. We can translate the above mathematical conditions into three different ways of describing manifolds:

(Local Parametrization) a manifold is a topological space that locally resembles Euclidean space near each point. For example, consider a two dimensional grid. We would not be able to deform this grid to have the shape of a sphere under any circumstance, but at each \(\textbf{local}\) point on the manifold, we can approximate the space around the point with a Euclidean grid.

(Local Implicit Function) a manifold can be understood as the set of points that verify a constraint defined by an implicit equation, given by the function \(f\). (see hypersphere example)

(Local Graph) a manifold can be understood as a d-dimensional surface described by d variables, i.e. by the “graph” of a smooth function \(f: (x_1,...,x_d) \to f(x_1,...,x_d)\). This local graph cannot be applied globally (to the whole manifold) because a function \(f\) must have only one output per set of inputs, and if we were to try to describe the whole manifold with such a graph, then many manifolds would have more than one output for one set of inputs.

The figure below provides a visualization of each of these definitions. The first figure is a visualization of “local parametrization”, the second is a visualization of “local implicit function” and the third is a visualization of “local graph”:

Each one of these definitions of a manifold is important because each one gives us an additional way to describe a manifold. This is useful because some manifolds are much easier described (or computationally implemented) by one of these definitions than the others. For example, any three of these definitions can be used to describe a hypersphere, but a hypersphere is most easily implemented using definition (2).

## 2.3 Hypersphere example:#

Here, we will consider how to prove that a hypersphere is a manifold using the second condition in the definition of a manifold above (2).

A \(\textbf{hypersphere}\) is any of a set of objects (n-dimensional spheres) resulting from the generalization of a one-dimensional circle and a two-dimensional sphere. The dimension of the manifold is n and is equal to the number of degrees of freedom of a point moving on In n-dimensional space, a hypersphere is the set of all points that are a given distance (called the radius) from a given point (called the center). For example, a 2-dimensional hypersphere in 3-dimensions describes all the points in 3D space that lie on the surface of a sphere. In other words, a hypersphere describes all of the points that lie in three dimensions that have two degrees of freedom \((\phi, \theta)\).

\(\textbf{How do we know that a hypersphere is a manifold?}\) Well, we know from the definition of a hypersphere that a hypersphere is the set of all points that are a given distance from the center of your coordinate system. For example in 3-dimensions, the surface \(S\) of a sphere with radius 1 can be described by the relation:

\(|x|^{2} = 1\)

This simply says that all the points must be a distance of 1 away from the center. (see above figure (b) and take a = 1).

We can now define a function

\(f(x) = |x|^{2} - 1\)

We know that this function \(f\) will always equal zero for all points that lie on the surface \(S\) because in order for a point to fall on the surface \(S\) the condition \(|x|^{2} = 1\) must be true. In other words,

\(x \in S \iff f(x) = 0\)

which tells us that

\(x \in S \iff x \in f^{-1}\)({0})

This last line matches the definition of a manifold: \(M = f^{-1}\)({0}) where \(M\) is the set of points \(x\) that satisfy the condition \(|x|^{2} = 1\). Therefore, the set of points that satisfy the condition \(|x|^{2} = 1\) form a manifold.

## 2.4 More examples of manifolds:#

A hypersphere is one type of manifold, but there are many other types of manifolds that are commonly used and seen in nature. We will show a few examples here to help you build intuition about what a manifold is and what a manifold can look like.

For example, a torus (the surface of a donut) is a manifold.

A klein bottle (the surface of the shape shown below) is also a manifold.

Perhaps a more familiar example for those in physics: Minkowski space (or Minkowski spacetime) is a combination of three-dimensional Euclidean space and time into a four-dimensional manifold (where x, y, z, and time are the dimensions of the space). This manifold is theorized to take one of three forms. A 4-dimensional manifold of: (a) flat curvature, (b) positive curvature, or (c) negative curvature.

# 3. Why do we care about manifolds?#

\(\textbf{Manifolds are important because data in nature "naturally falls on manifolds"}\), and as we discussed in the introduction section, knowing the manifold that a data set belongs to may give you more predictive power and a better understanding of the data’s evolution.

\(\textbf{What does it mean for data to "naturally fall on a manifold"?}\) In nature, data are often subject to constraints, and these constraints force the data to lie on manifolds. For example, consider the position of cities on the earth.

The cities are subject to the following constraints: 1) they cannot fly above the surface of the earth because gravity holds them down and 2) they cannot sink down into the earth because the surface of the earth holds them up. Therefore, they are constrained to move (or not move) on the surface of a 3-dimensional sphere. This space that they are confined to exist in should sound familiar– it is the space of a 2-dimensional hypersphere!

Other data falls on manifolds in similar ways: the data is subject to certain constraints, which forces it to fall on a manifold.

# 4. What does the `Manifold`

class do?#

The `Manifold`

class describes different types of manifolds. The `Manifold`

class and its subclasses implement methods that establish the properties of different types of manifolds, and the manifolds that exist under these subclasses inherit the properties implemented in their respective subclasses.

Note: The words “class”, “subclass”, “methods” refer to object oriented programming.

The manifold class also provides ways of checking whether a vector is \(\textit{tangent}\) to a manifold and ways of projecting a vector to the \(\textit{tangent space}\) of a manifold. In order to better describe what this means, we will now define \(\textit{tangent vectors}\) and \(\textit{tangent spaces}\).

## 4.1 Tangent Vectors#

Here we will provide a brief reminder of tangent vectors as a way to lay the groundwork for discussing tangent spaces.

A \(\textbf{tangent vector}\) is a vector that is “tangent” to a curve or surface at a given point. When a vector is “tangent” to a curve, this means that the vector has the same slope as the curve does at that point.

## 4.2 Tangent Spaces#

Now we will introduce the concept of a “tangent space”. The tangent space at a certain point on a manifold is comprised of all of the possible tangent vectors that exist at that point. For example, if you are considering the tangent space on a curve, then the only possible tangent vectors are tangent vectors that point forward and backward along the line (a). However, if you are considering the tangent space at a point on a surface, then the tangent vectors can point forward, backward, left, right, and everywhere in between, and the set of all tangent vectors forms a plane (b).

Thus, the tangent space of a 1-dimensional manifold (curve) is also one dimensional, and the tangent space of a 2-dimensional manifold (a 2-dimensional surface) is also 2-dimensional.

Similarly, for every n-dimensional manifold, there exists an n-dimensional tangent space at each point on the manifold, and the tangent space is comprised of all possible tangent vectors on that manifold.

# 5. How is the `Manifold`

class structured?#

The hierarchical structure of the classes inheriting from the `Manifold`

parent class is as follows (this Figure is a courtesy of Nicolas Guigui):

As discussed in the previous section, one of the primary purposes of the `Manifold`

class is to hold information about various types of manifolds. Rules that are universally true for all manifolds are implemented in methods in the parent class `Manifold`

. Rules that are true for some types of manifolds are implemented in the subclasses of `Manifold`

: `LevelSet`

, `VectorSpaceOpenSet`

, `FiberBundle`

, `ProductManifold`

, `VectorSpace`

, `MatrixLieAlgebra`

, and `MatrixLieGroup`

.
Specific types of manifolds are described in methods within these subclasses.

In this notebook, we will focus on describing the subclasses pertinent to the geometry module of geomstats: `LevelSet`

, `VectorSpaceOpenSet`

, `ProductManifold`

and `VectorSpace`

.

In the following subsections, we will discuss the methods and ideas implemented in the parent class and its subclasses.

## 5.1 The Parent Class: `Manifold`

#

The Manifold parent class is an abstract base class which provides the minimal skeleton of attributes and methods expected in its subclasses. Note that the methods of the abstract parent class are declared, but they contain no implementation, and they are overridden by the subclasses. The properties that are declared in the `Manifold`

class are properties that all types of manifold must possess. For example, the following methods and attributes are implemented in `Manifold`

:

`dim`

: \(\textit{attribute}\). the dimension of the manifold. “How many coordinates are necessary to fully describe the manifold?”`belongs()`

: \(\textit{method}\). evaluates whether a given element belongs to that manifold`is_tangent()`

: \(\textit{method}\). evaluates whether a given vector is a tangent vector at a given point`random_point()`

: \(\textit{method}\). generates a random point that lies on the manifold

While the abstract methods in Manifold do not contain any implementation, the methods of the subclasses of `Manifold`

, such as the `Hypersphere`

subclass, \(\textit{are}\) implemented and can be run. We will now exemplify this in the following section

### 5.1.1 Examples of Using `Manifold`

’s Attributes and Methods in the Subclass: `Hypersphere`

#

\(\textbf{Attributes:}\)

`dim`

: If we build a hypersphere of dimension 2 with the following code, we can check that `sphere.dim`

gives back 2. Run the following code to verify this:

```
In [2]:
```

```
from geomstats.geometry.hypersphere import Hypersphere
sphere = Hypersphere(dim=2)
print(f"The dimension of the sphere is {sphere.dim}")
```

```
The dimension of the sphere is 2
```

\(\textbf{Methods:}\)

`belongs`

: We can re-use the sphere we just built (called “sphere”) and verify that the point (0, 0, 1) belongs to that sphere (it is the north pole). Run the following code to verify this using Geomstats.

```
In [3]:
```

```
import geomstats.backend as gs
sphere.belongs(gs.array([0, 0, 1]))
```

```
Out [3]:
```

```
True
```

`is_tangent`

: The vector (1, 1, 0 ) is tangent to the sphere at the north pole, since it does not have a vertical component (last component is equal to 0). Run the following code to verify this using Geomstats.

```
In [4]:
```

```
sphere.is_tangent(vector=gs.array([1, 1, 0]), base_point=gs.array([0, 0, 1]))
```

```
Out [4]:
```

```
True
```

`random_point`

: Now, we will use `random_point`

to generate a random point, and then we will use `belongs`

to prove that this random point belongs to the sphere.

```
In [5]:
```

```
from geomstats.geometry.hypersphere import Hypersphere
rp = Hypersphere.random_point(sphere)
sphere.belongs(rp)
```

```
Out [5]:
```

```
True
```

### 5.1.2 The Full `Manifold`

Class Code#

You can see all of the methods in the `Manifold`

parent class by running the following code. Observe the abstract methods denoted with “@abc.abstractmethod” that do not contain any implementation, but serve as template for the subclasses to overwrite. The code of `Manifold`

can also be found here.

```
In [6]:
```

```
import inspect
from geomstats.geometry.manifold import Manifold
for line in inspect.getsourcelines(Manifold)[0]:
line = line.replace("\n", "")
print(line)
```

```
class Manifold(abc.ABC):
r"""Class for manifolds.
Parameters
----------
dim : int
Dimension of the manifold.
shape : tuple of int
Shape of one element of the manifold.
Optional, default : None.
intrinsic : bool
Coordinate type.
equip : bool
If True, equip space with default metric.
Attributes
----------
point_ndim : int
Dimension of point array.
"""
def __init__(
self,
dim,
shape,
intrinsic=True,
equip=True,
):
geomstats.errors.check_integer(dim, "dim")
if not isinstance(shape, tuple):
raise ValueError("Expected a tuple for the shape argument.")
self.dim = dim
self.shape = shape
self.intrinsic = intrinsic
self.point_ndim = len(self.shape)
if equip:
self.equip_with_metric()
def equip_with_metric(self, Metric=None, **metric_kwargs):
"""Equip manifold with a Riemannian metric.
Parameters
----------
Metric : RiemannianMetric object or instance or ScalarProductMetric instance
If None, default metric will be used.
"""
if Metric is None:
out = self.default_metric()
if isinstance(out, tuple):
Metric, kwargs = out
kwargs.update(metric_kwargs)
metric_kwargs = kwargs
else:
Metric = out
if inspect.isclass(Metric):
self.metric = Metric(self, **metric_kwargs)
else:
if self.metric._space is not self:
raise ValueError(
"Cannot equip space with metric instantiated with another space."
)
self.metric = Metric
return self
def equip_with_group_action(self, group_action):
"""Equip manifold with group action.
Parameters
----------
group_action : str
Group action.
"""
self.group_action = group_action
return self
def equip_with_quotient(self):
"""Equip manifold with quotient structure.
Creates attributes `quotient` and `fiber_bundle` or `aligner` (
`aligner` is used in quotient contexts where the notion
of fiber bundle is not defined.).
Returns
-------
quotient : Manifold or None
Quotient space equipped with a quotient metric.
"""
if not _QuotientStructureRegistry.has_quotient(self):
raise ValueError("No quotient structure defined for this manifold.")
FiberBundle_, QuotientMetric_ = (
_QuotientStructureRegistry.get_fiber_bundle_and_quotient_metric(
self,
)
)
fiber_bundle = FiberBundle_(total_space=self)
if hasattr(fiber_bundle, "riemannian_submersion"):
self.fiber_bundle = fiber_bundle
else:
self.aligner = fiber_bundle
if QuotientMetric_ is None:
return
self.quotient = self.new(equip=False)
self.quotient.equip_with_metric(QuotientMetric_, total_space=self)
return self.quotient
@abc.abstractmethod
def belongs(self, point, atol=gs.atol):
"""Evaluate if a point belongs to the manifold.
Parameters
----------
point : array-like, shape=[..., *point_shape]
Point to evaluate.
atol : float
Absolute tolerance.
Optional, default: backend atol.
Returns
-------
belongs : array-like, shape=[...,]
Boolean evaluating if point belongs to the manifold.
"""
@abc.abstractmethod
def is_tangent(self, vector, base_point=None, atol=gs.atol):
"""Check whether the vector is tangent at base_point.
Parameters
----------
vector : array-like, shape=[..., *point_shape]
Vector.
base_point : array-like, shape=[..., *point_shape]
Point on the manifold.
atol : float
Absolute tolerance.
Optional, default: backend atol.
Returns
-------
is_tangent : bool
Boolean denoting if vector is a tangent vector at the base point.
"""
@abc.abstractmethod
def to_tangent(self, vector, base_point=None):
"""Project a vector to a tangent space of the manifold.
Parameters
----------
vector : array-like, shape=[..., *point_shape]
Vector.
base_point : array-like, shape=[..., *point_shape]
Point on the manifold.
Returns
-------
tangent_vec : array-like, shape=[..., *point_shape]
Tangent vector at base point.
"""
@abc.abstractmethod
def random_point(self, n_samples=1, bound=1.0):
"""Sample random points on the manifold according to some distribution.
If the manifold is compact, preferably a uniform distribution will be used.
Parameters
----------
n_samples : int
Number of samples.
Optional, default: 1.
bound : float
Bound of the interval in which to sample for non compact manifolds.
Optional, default: 1.
Returns
-------
samples : array-like, shape=[..., *point_shape]
Points sampled on the manifold.
"""
def regularize(self, point):
"""Regularize a point to the canonical representation for the manifold.
Parameters
----------
point : array-like, shape=[..., dim]
Point.
Returns
-------
regularized_point : array-like, shape=[..., *point_shape]
Regularized point.
"""
return gs.copy(point)
def random_tangent_vec(self, base_point=None, n_samples=1):
"""Generate random tangent vec.
This method is not recommended for statistical purposes,
as the tangent vectors generated are not drawn from a
distribution related to the Riemannian metric.
Parameters
----------
n_samples : int
Number of samples.
Optional, default: 1.
base_point : array-like, shape={[n_samples, *point_shape], [*point_shape,]}
Point.
Returns
-------
tangent_vec : array-like, shape=[..., *point_shape]
Tangent vec at base point.
"""
if (
n_samples > 1
and base_point is not None
and base_point.ndim > len(self.shape)
and n_samples != len(base_point)
):
raise ValueError(
"The number of base points must be the same as the "
"number of samples, when the number of base points is different from 1."
)
batch_size = () if n_samples == 1 else (n_samples,)
return self.to_tangent(
gs.random.normal(size=batch_size + self.shape), base_point
)
def projection(self, point):
"""Project a point to the manifold.
Parameters
----------
point: array-like, shape[..., *point_shape]
Point.
Returns
-------
point: array-like, shape[..., *point_shape]
Point.
"""
if self.intrinsic:
return gs.copy(point)
raise NotImplementedError("`projection` is not implemented yet")
```

## 5.2 `VectorSpaceOpenSet`

#

Earlier in the notebook, we were able to say that a set of points is a manifold if it satisfied one of three constraints. We also said that every manifold can be described by any three of these definitions, and the choice of definition is merely a question of which definition is most convenient. `VectorSpaceOpenSet`

provides a way of describing manifolds with local parametrization, which was labeled (1) on our definition list.

One such way to describe a manifold is with the concept of an \(\textbf{Open Set}\): a manifold is the open sets of a d-dimensional vector space, called \(\textbf{ambient space}\).

### 5.2.1 What is an Open Set?#

Intuitively, an open set is a group of numbers that does not include points on the boundary of whatever they are describing. For example, if you were to take the set of all points between a and b but \(\textbf{not}\) include the values a and b, this would be an open set (shown in figure a). If you were to take the set of all points between a and b and include the values a and b, this would be a closed set (shown in figure b)

The above examples showed open and closed sets in one dimension. Similarly, in two dimensions, open sets can be defined as sets which do not contain their boundaries. For example, the inside of the sphere, i.e. the ball without its boundary, is a manifold that is an open set. The image below shows an example of an open set (a) and a closed set (b) in two dimensions.

### 5.2.2 What Methods are Implemented in `VectorSpaceOpenSet`

?#

If we know that a manifold is conveniently described as an open set, then some of the manifold’s abstract methods can be rewritten in a specific form. For example, `VectorSpaceOpenSet`

implements the methods:

`projection()`

: a method to project any d-dimensional vector to the manifold.`is_tangent()`

: checks whether the input vector is tangent at the input point.`to_tangent()`

: projects a vector to a tangent space of the manifold.

Note that we do not (yet) specify which manifold we are talking about (whether it is a sphere or another surface), we are just saying that we are looking at some manifold that can be described as an open set

Run the code below to see the contents of the `VectorSpaceOpenSet`

class.

```
In [7]:
```

```
import inspect
from geomstats.geometry.base import VectorSpaceOpenSet
for line in inspect.getsourcelines(VectorSpaceOpenSet)[0]:
line = line.replace("\n", "")
print(line)
```

```
class VectorSpaceOpenSet(OpenSet, abc.ABC):
"""Class for manifolds that are open sets of a vector space.
In this case, tangent vectors are identified with vectors of the embedding
space.
Parameters
----------
embedding_space: VectorSpace
Embedding space that contains the manifold.
"""
def is_tangent(self, vector, base_point=None, atol=gs.atol):
"""Check whether the vector is tangent at base_point.
Parameters
----------
vector : array-like, shape=[..., *point_shape]
Vector.
base_point : array-like, shape=[..., *point_shape]
Point on the manifold.
atol : float
Absolute tolerance.
Optional, default: backend atol.
Returns
-------
is_tangent : bool
Boolean denoting if vector is a tangent vector at the base point.
"""
is_tangent = self.embedding_space.belongs(vector, atol)
if base_point is not None and base_point.ndim > vector.ndim:
return gs.broadcast_to(is_tangent, base_point.shape[: -self.point_ndim])
return is_tangent
def to_tangent(self, vector, base_point=None):
"""Project a vector to a tangent space of the manifold.
Parameters
----------
vector : array-like, shape=[..., *point_shape]
Vector.
base_point : array-like, shape=[..., *point_shape]
Point on the manifold.
Returns
-------
tangent_vec : array-like, shape=[..., *point_shape]
Tangent vector at base point.
"""
tangent_vec = self.embedding_space.projection(vector)
if base_point is not None and base_point.ndim > vector.ndim:
return gs.broadcast_to(tangent_vec, base_point.shape)
return tangent_vec
```

## 5.3 `LevelSet`

#

### 5.3.1 What is a Level Set?#

Another elementary class of manifolds are \(\textbf{Level Sets}\). A level set is the set of values \(x\) for which a function f(x) is equal to a given constant. In other words, a level set is a set of curves for which the function describing a manifold is constant along that curve.

In the same way that `VectorSpaceOpenSet`

is an implementation of the first definition of a manifold (Local Parametrization), `LevelSet`

is an implementation of the second definition of a manifold (Local Implicit Function). a level set is a set of points for which the function \(f\) takes the exact same value. This value is called the “level”, and does not need to be a scalar, it could also be a vector.

For example, consider a hypersphere in three dimensional space. Each of the concentric spheres is a 2-dimensional manifold, each corresponding to a different level \((r1, r2,r3)\).

You can see here that the subclass `Hypersphere`

is indeed implemented as a `LevelSet`

.

### 5.3.2 `LevelSet`

in Geomstats#

You can run the code below to see the contents of the `LevelSet`

class. The methods of LevelSet are similar to the methods of VectorSpaceOpenSet , but their implementation is different. Recall from the (Local Implicit Function) definition, which `LevelSet`

is an implementation of, that a manifold of this type must satisfy:

“For every \(p \in M\), there exists an open set \(U \in \mathbb{R}^{N}\) and a smooth map \(f: U \to \mathbb{R}^{N-d}\) that is a submersion at p, such that \(U \cap M = f^{-1}\)({0}).”

`VectorSpaceOpenSet`

is an implementation of the first manifold definition (Local Parametrization), and therefore need not follow the (Local Implicit Function) rule above. This is the reason that, `VectorSpaceOpenSet`

methods and `LevelSet`

methods are implemented differently. As an example of these different implementations: observe the implementation of the belongs methods in `LevelSet`

. For a general level set, in order to verify if a point belongs to the level set, we should verify
that the (Local Implicit Function) definition constraint is met, which is done with the line

`constraint = gs.isclose(self.submersion(point), value, atol=atol)`

```
In [8]:
```

```
import inspect
from geomstats.geometry.base import LevelSet
for line in inspect.getsourcelines(LevelSet)[0]:
line = line.replace("\n", "")
print(line)
```

```
class LevelSet(Manifold, abc.ABC):
"""Class for manifolds embedded in a vector space by a submersion.
Parameters
----------
intrinsic : bool
Coordinates type.
"""
def __init__(self, intrinsic=False, shape=None, **kwargs):
self.embedding_space = self._define_embedding_space()
if shape is None:
shape = self.embedding_space.shape
super().__init__(intrinsic=intrinsic, shape=shape, **kwargs)
@abc.abstractmethod
def _define_embedding_space(self):
"""Define embedding space of the manifold.
Returns
-------
embedding_space : Manifold
Instance of Manifold.
"""
@abc.abstractmethod
def submersion(self, point):
r"""Submersion that defines the manifold.
:math:`\mathrm{submersion}(x)=0` defines the manifold.
Parameters
----------
point : array-like, shape=[..., *point_shape]
Returns
-------
submersed_point : array-like
"""
@abc.abstractmethod
def tangent_submersion(self, vector, point):
"""Tangent submersion.
Parameters
----------
vector : array-like, shape=[..., *point_shape]
point : array-like, shape=[..., *point_shape]
Returns
-------
submersed_vector : array-like
"""
def belongs(self, point, atol=gs.atol):
"""Evaluate if a point belongs to the manifold.
Parameters
----------
point : array-like, shape=[..., *point_shape]
Point to evaluate.
atol : float
Absolute tolerance.
Optional, default: backend atol.
Returns
-------
belongs : array-like, shape=[...,]
Boolean evaluating if point belongs to the manifold.
"""
belongs = self.embedding_space.belongs(point, atol)
if not gs.any(belongs):
return belongs
submersed_point = self.submersion(point)
n_batch = gs.ndim(point) - len(self.shape)
axis = tuple(range(-len(submersed_point.shape) + n_batch, 0))
if gs.is_complex(submersed_point):
constraint = gs.isclose(submersed_point, 0.0 + 0.0j, atol=atol)
else:
constraint = gs.isclose(submersed_point, 0.0, atol=atol)
if axis:
constraint = gs.all(constraint, axis=axis)
return gs.logical_and(belongs, constraint)
def is_tangent(self, vector, base_point, atol=gs.atol):
"""Check whether the vector is tangent at base_point.
Parameters
----------
vector : array-like, shape=[..., *point_shape]
Vector.
base_point : array-like, shape=[..., *point_shape]
Point on the manifold.
atol : float
Absolute tolerance.
Optional, default: backend atol.
Returns
-------
is_tangent : bool
Boolean denoting if vector is a tangent vector at the base point.
"""
belongs = self.embedding_space.is_tangent(vector, base_point, atol)
if not gs.any(belongs):
return belongs
submersed_vector = self.tangent_submersion(vector, base_point)
n_batch = len(get_batch_shape(self.point_ndim, base_point, vector))
axis = tuple(range(-len(submersed_vector.shape) + n_batch, 0))
constraint = gs.isclose(submersed_vector, 0.0, atol=atol)
if axis:
constraint = gs.all(constraint, axis=axis)
return gs.logical_and(belongs, constraint)
def intrinsic_to_extrinsic_coords(self, point_intrinsic):
"""Convert from intrinsic to extrinsic coordinates.
Parameters
----------
point_intrinsic : array-like, shape=[..., *point_shape]
Point in the embedded manifold in intrinsic coordinates.
Returns
-------
point_extrinsic : array-like, shape=[..., *embedding_space.point_shape]
Point in the embedded manifold in extrinsic coordinates.
"""
raise NotImplementedError("intrinsic_to_extrinsic_coords is not implemented.")
def extrinsic_to_intrinsic_coords(self, point_extrinsic):
"""Convert from extrinsic to intrinsic coordinates.
Parameters
----------
point_extrinsic : array-like, shape=[..., *embedding_space.point_shape]
Point in the embedded manifold in extrinsic coordinates,
i. e. in the coordinates of the embedding manifold.
Returns
-------
point_intrinsic : array-lie, shape=[..., *point_shape]
Point in the embedded manifold in intrinsic coordinates.
"""
raise NotImplementedError("extrinsic_to_intrinsic_coords is not implemented.")
```

## 5.4 `VectorSpace`

#

A \(\textit{vector space}\) is a special case of manifold where a vector exists at each point.

Many operations defined on manifolds become simpler when the manifold is actually a vector space. Additionally, the ambient spaces used to define manifolds as open or closed sets are often vector spaces. Thus, we provide an implementation `VectorSpace`

.

This class does not provide another way of describing a manifold. This class is an abstract class that makes sure that the ambient/embedding vector space is compatible with the methods that are called in `VectorSpaceOpenSet`

and `LevelSet`

.

Actual manifolds are implemented as subclasses of this abstract method and must implement all the abstract methods defined in this class.

You can run the code below to see the contents of the `VectorSpace`

class. Observe that, in a `VectorSpace`

manifold, the manifold is comprised of vectors – and more specifically, tangent vectors. Thus, in the implementation, we can see that in the `is_tangent()`

method: in order to check that a vector is a tangent to the manifold, we only need to check that the vector belongs to the vector space.

```
In [9]:
```

```
import inspect
from geomstats.geometry.base import VectorSpace
for line in inspect.getsourcelines(VectorSpace)[0]:
line = line.replace("\n", "")
print(line)
```

```
class VectorSpace(Manifold, abc.ABC):
"""Abstract class for vector spaces.
Parameters
----------
shape : tuple
Shape of the elements of the vector space. The dimension is the
product of these values by default.
"""
def __init__(self, shape, dim=None, **kwargs):
if dim is None:
dim = math.prod(shape)
super().__init__(dim=dim, shape=shape, **kwargs)
self._basis = None
def belongs(self, point, atol=gs.atol):
"""Evaluate if the point belongs to the vector space.
This method checks the shape of the input point.
Parameters
----------
point : array-like, shape=[.., *point_shape]
Point to test.
atol : float
Unused here.
Returns
-------
belongs : array-like, shape=[...,]
Boolean evaluating if point belongs to the space.
"""
belongs = self.shape == point.shape[-self.point_ndim :]
shape = point.shape[: -self.point_ndim]
if belongs:
return gs.ones(shape, dtype=bool)
return gs.zeros(shape, dtype=bool)
def is_tangent(self, vector, base_point=None, atol=gs.atol):
"""Check whether the vector is tangent at base_point.
Tangent vectors are identified with points of the vector space so
this checks the shape of the input vector.
Parameters
----------
vector : array-like, shape=[..., *point_shape]
Vector.
base_point : array-like, shape=[..., *point_shape]
Point in the vector space.
atol : float
Absolute tolerance.
Optional, default: backend atol.
Returns
-------
is_tangent : array-like, shape=[...,]
Boolean denoting if vector is a tangent vector at the base point.
"""
belongs = self.belongs(vector, atol)
if base_point is not None and base_point.ndim > vector.ndim:
return gs.broadcast_to(belongs, base_point.shape[: -self.point_ndim])
return belongs
def to_tangent(self, vector, base_point=None):
"""Project a vector to a tangent space of the vector space.
This method is for compatibility and returns vector.
Parameters
----------
vector : array-like, shape=[..., *point_shape]
Vector.
base_point : array-like, shape=[..., *point_shape]
Point in the vector space
Returns
-------
tangent_vec : array-like, shape=[..., *point_shape]
Tangent vector at base point.
"""
tangent_vec = self.projection(vector)
if base_point is not None and base_point.ndim > vector.ndim:
return gs.broadcast_to(tangent_vec, base_point.shape)
return tangent_vec
def random_point(self, n_samples=1, bound=1.0):
"""Sample in the vector space with a uniform distribution in a box.
Parameters
----------
n_samples : int
Number of samples.
Optional, default: 1.
bound : float
Side of hypercube support of the uniform distribution.
Optional, default: 1.0
Returns
-------
point : array-like, shape=[..., dim]
Sample.
"""
size = (self.dim,)
if n_samples != 1:
size = (n_samples,) + size
return bound * (gs.random.rand(*size) - 0.5) * 2
def random_tangent_vec(self, base_point=None, n_samples=1):
"""Generate random tangent vec.
This method is not recommended for statistical purposes, as the
tangent vectors generated are not drawn from a distribution related
to the Riemannian metric.
Parameters
----------
n_samples : int
Number of samples.
Optional, default: 1.
base_point : array-like, shape={[n_samples, *point_shape], [*point_shape,]}
Point.
Returns
-------
tangent_vec : array-like, shape=[..., *point_shape]
Tangent vec at base point.
"""
if base_point is None:
return self.random_point(n_samples)
if (
n_samples > 1
and base_point.ndim > self.point_ndim
and n_samples != base_point.shape[0]
):
raise ValueError(
"The number of base points must be the same as the "
"number of samples, when the number of base points is different from 1."
)
if n_samples == 1 and base_point.ndim > self.point_ndim:
n_samples = base_point.shape[0]
return self.random_point(n_samples)
@property
def basis(self):
"""Basis of the vector space."""
if self._basis is None:
self._basis = self._create_basis()
return self._basis
@abc.abstractmethod
def _create_basis(self):
"""Create a canonical basis."""
```

## 5.5 `ProductManifold`

#

New manifolds can also be created by composing existing manifolds. A product manifold combines two manifolds \(M_1\) and \(M_2\) to create a new manifold \(M_3\), and more generally \(N\) manifolds \(M1, ..., M_N\) can be combined to create a new manifold \(M\).

### 5.5.1 Example#

For example, we can create a product of two spheres as:

```
In [10]:
```

```
from geomstats.geometry.product_manifold import ProductManifold
sphere1 = Hypersphere(dim=2)
sphere2 = Hypersphere(dim=2)
product_of_two_spheres = ProductManifold([sphere1, sphere2])
```

If we generate a random point on this product manifold, we will get 6 coordinates: the first 3 represent a point on the first sphere, and the last 3 represent a point on the second sphere:

```
In [11]:
```

```
product_of_two_spheres.random_point()
```

```
Out [11]:
```

```
array([[-0.3861891 , -0.52946814, -0.75532871],
[-0.14573203, 0.15373756, 0.97730596]])
```

You can run the code below to see the contents of the `ProductManifold`

class.

```
In [12]:
```

```
import inspect
from geomstats.geometry.product_manifold import ProductManifold
for line in inspect.getsourcelines(ProductManifold)[0]:
line = line.replace("\n", "")
print(line)
```

```
class ProductManifold(_IterateOverFactorsMixins, Manifold):
"""Class for a product of manifolds M_1 x ... x M_n.
In contrast to the classes NFoldManifold, Landmarks, or DiscretizedCurves,
the manifolds M_1, ..., M_n need not be the same, nor of
same dimension, but the list of manifolds needs to be provided.
Parameters
----------
factors : tuple
Collection of manifolds in the product.
point_ndim : int or None
If None, defaults to 1, unless all factors have the same shape.
"""
def __init__(self, factors, point_ndim=None, equip=True):
factors = tuple(factors)
factor_dims = [factor.dim for factor in factors]
dim = sum(factor_dims)
shape = _find_product_shape(factors, point_ndim)
intrinsic = all(factor.intrinsic for factor in factors)
if not intrinsic:
factor_embedding_spaces = [
(
manifold.embedding_space
if hasattr(manifold, "embedding_space")
else manifold
)
for manifold in factors
]
# TODO: need to revisit due to removal of scales
self.embedding_space = ProductManifold(
factor_embedding_spaces, point_ndim, equip=False
)
cum_index = (
gs.cumsum(factor_dims)[:-1]
if intrinsic
else self.embedding_space._cum_index
)
super().__init__(
factors=factors,
cum_index=cum_index,
pool_outputs=True,
has_mixed_fields=_has_mixed_fields(factors),
dim=dim,
shape=shape,
intrinsic=intrinsic,
equip=equip,
)
@staticmethod
def default_metric():
"""Metric to equip the space with if equip is True."""
return ProductRiemannianMetric
def _pool_outputs_from_function(self, outputs):
"""Collect outputs for each product to be returned.
If each element of the output is a boolean array of the same shape, test along
the list whether all elements are True and return a boolean array of the same
shape.
Otherwise, if each element of the output has a shape compatible with points of
the corresponding factor, an attempt is made to map the list of points to a
point in the product by embed_to_product.
Parameters
----------
outputs : list
A list of outputs which must be pooled
Returns
-------
pooled_output : array-like, shape {(...,), (..., self.shape)}
"""
# TODO: simplify after cleaning gs.squeeze
all_arrays = gs.all([gs.is_array(factor_output) for factor_output in outputs])
if (
all_arrays
and _all_equal([factor_output.shape for factor_output in outputs])
and gs.all([gs.is_bool(factor_output) for factor_output in outputs])
or (not all_arrays)
):
outputs = gs.stack([gs.array(factor_output) for factor_output in outputs])
outputs = gs.all(outputs, axis=0)
return outputs
try:
return self.embed_to_product(outputs)
except geomstats.errors.ShapeError:
raise RuntimeError(
"Could not combine outputs - they are not points of the individual"
" factors."
)
except ValueError:
raise RuntimeError(
"Could not combine outputs, probably because they could"
" not be concatenated or stacked."
)
def belongs(self, point, atol=gs.atol):
"""Test if a point belongs to the manifold.
Parameters
----------
point : array-like, shape=[..., {dim, embedding_space.dim, \
[n_manifolds, dim_each]}]
Point.
atol : float,
Tolerance.
Returns
-------
belongs : array-like, shape=[...,]
Boolean evaluating if the point belongs to the manifold.
"""
belongs = self._iterate_over_factors("belongs", {"point": point, "atol": atol})
return belongs
def regularize(self, point):
"""Regularize the point into the manifold's canonical representation.
Parameters
----------
point : array-like, shape=[..., {dim, embedding_space.dim, \
[n_manifolds, dim_each]}]
Point to be regularized.
Returns
-------
regularized_point : array-like, shape=[..., {dim, embedding_space.dim, \
[n_manifolds, dim_each]}]
Point in the manifold's canonical representation.
"""
regularized_point = self._iterate_over_factors("regularize", {"point": point})
return regularized_point
def random_point(self, n_samples=1, bound=1.0):
"""Sample in the product space from the product distribution.
The distribution used is the product of the distributions used by the
random_point methods of each individual factor manifold.
Parameters
----------
n_samples : int, optional
Number of samples.
bound : float
Bound of the interval in which to sample for non compact manifolds.
Optional, default: 1.
Returns
-------
samples : array-like, shape=[..., {dim, embedding_space.dim, \
[n_manifolds, dim_each]}]
Points sampled from the manifold.
"""
samples = self._iterate_over_factors(
"random_point", {"n_samples": n_samples, "bound": bound}
)
return samples
def random_tangent_vec(self, base_point, n_samples=1):
"""Sample on the tangent space from the product distribution.
The distribution used is the product of the distributions used by the
random_tangent_vec methods of each individual factor manifold.
Parameters
----------
base_point : array-like, shape=[..., n, n]
Base point of the tangent space.
Optional, default: None.
n_samples : int
Number of samples.
Optional, default: 1.
Returns
-------
samples : array-like, shape=[..., {dim, embedding_space.dim, \
[n_manifolds, dim_each]}]
Points sampled in the tangent space of the product manifold at base_point.
"""
samples = self._iterate_over_factors(
"random_tangent_vec", {"base_point": base_point, "n_samples": n_samples}
)
return samples
def projection(self, point):
"""Project a point onto product manifold.
Parameters
----------
point : array-like, shape=[..., {dim, embedding_space.dim, \
[n_manifolds, dim_each]}]
Point in product manifold.
Returns
-------
projected : array-like, shape=[..., {dim, embedding_space.dim, \
[n_manifolds, dim_each]}]
Projected point.
"""
projected_point = self._iterate_over_factors("projection", {"point": point})
return projected_point
def to_tangent(self, vector, base_point):
"""Project a vector to a tangent space of the manifold.
Parameters
----------
vector : array-like, shape=[..., dim]
Vector.
base_point : array-like, shape=[..., dim]
Point on the manifold.
Returns
-------
tangent_vec : array-like, shape=[..., dim]
Tangent vector at base point.
Notes
-----
The tangent space of the product manifold is the direct sum of
tangent spaces.
"""
tangent_vec = self._iterate_over_factors(
"to_tangent", {"base_point": base_point, "vector": vector}
)
return tangent_vec
def is_tangent(self, vector, base_point=None, atol=gs.atol):
"""Check whether the vector is tangent at base_point.
The tangent space of the product manifold is the direct sum of
tangent spaces.
Parameters
----------
vector : array-like, shape=[..., dim]
Vector.
base_point : array-like, shape=[..., dim]
Point on the manifold.
Optional, default: None
atol : float
Absolute tolerance.
Optional, default: backend atol.
Returns
-------
is_tangent : bool
Boolean denoting if vector is a tangent vector at the base point.
"""
is_tangent = self._iterate_over_factors(
"is_tangent", {"base_point": base_point, "vector": vector, "atol": atol}
)
return is_tangent
```

# Conclusion#

Key takeaways from this notebook:

There are three ways to define a manifold, and each of these three definitions provide a different option for manifold implementation.

Learning to work with manifolds is useful because many datasets naturally lie on manifolds.

The

`Manifold`

class stores information about many types of manifolds. Each of the subclasses implements a different \(\textit{type}\) of manifold (for example,`LevelSet`

), and specific manifolds (for example,`Hypersphere`

) are implemented within these subclasses.