Thursday, April 18, 2024

Riemann Curvature in Python

Target audience: Advanced
Estimated reading time: 7'

Beyond theoretical physics, the Riemann curvature tensor is essential for understanding motion planning in robotics, object recognition in computer vision, visualization and training of Physics-Inspired Neural Networks.
This article describes the curvature tensor, its derivative and its implementation in Python.


Table of contents
Introduction
Riemann curvature
       Curvature tensor
       Geodesic deviation
       Sectional curvature tensor
Implementation
       Hypersphere curvature
       Sectional curvature on hypersphere
References
Appendix
Follow me on LinkedIn

What you will learn: The intricacies of Riemannian metric curvature tensor and its implementation in Python using Geomstats library.

Notes: 

  • Environments: Python  3.10.10, Geomstats 2.7.0
  • This article assumes that the reader is somewhat familiar with differential and tensor calculus [ref 1, 2]. Please refer to the previous articles related to geometric learning [ref 4, 5, 6, 7].
  • Source code is available at  Github.com/patnicolas/Data_Exploration/manifolds
  • To enhance the readability of the algorithm implementations, we have omitted non-essential code elements like error checking, comments, exceptions, validation of class and method arguments, scoping qualifiers, and import statements.


Introduction

Geometric learning tackles the challenges posed by scarce data, high-dimensional spaces, and the requirement for independent representations in the creation of advanced machine learning models. The fundamental aim of studying Riemannian geometry is to comprehend and scrutinize the characteristics of curved spaces, which are not sufficiently explained by Euclidean geometry alone. 

In this article, we delve into the Riemann curvature tensors, which are based on parallel transport and the covariant derivative discussed in earlier articles.

This article is the eight part of our ongoing series focused on geometric learning. In this installment, we utilize the Geomstats Python library [ref. 3] and explore the hypersphere manifold, which was introduced in a previous piece, Geometric Learning in Python: Manifolds - Hypersphere  and is detailed in the Geomstats API [ref. 8]. 

I highly recommend watching the comprehensive series of 22 YouTube videos Tensor Calculus - Eigenchris  to familiarize yourself with fundamental concepts of differential geometry.  Summaries of my earlier articles on this topic can be found in the Appendix

Riemann curvature

Curvature tensor

Curvature measures the extent to which a geometric object like a curve strays from being straight or how a surface diverges from being flat. It can also be defined as a vector that incorporates both the direction and the magnitude of the curve's deviation. 

The Riemann curvature tensor is a widely used method for quantifying the curvature of Riemannian manifolds. It is calculated from the metric tensor and provides a tensor field at each point on the manifold.

Given a Riemann manifold M with a Levi-Civita connection and a tensor metric g, the Riemann curvature tensor related to vector fields X, Y and Z can be written:
\[\begin{matrix} \mathbf{\Omega (N) }\ \times  \ \mathbf{\Omega (N)} \ \times  \ \mathbf{\Omega (N)} \rightarrow \mathbf{\Omega (N)} \ \ \ \ \ \ \ \ \ \ \ \\ X, Y, Z \to R(X, Y)Z \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ R(X, Y)Z= \triangledown_{X} \triangledown_{Y}Z-\triangledown_{Y}\triangledown_{X}Z - \triangledown_{[X, Y]}Z \\ [X, Y] = X.\frac{\partial }{\partial Y} - Y.\frac{\partial }{\partial X} \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \end{matrix}\] [X, Y] is the lie bracket.
The Riemann curvature tensor can be written using the Christoffel symbols as: \[R_{jkl}^{i}=\frac{\partial \Gamma _{jl}^{i} } {\partial x^k} - \frac{\partial \Gamma _{jk}^{i} } {\partial x^l} + \Gamma_{rk}^{i} \Gamma_{jl}^{r} - \Gamma_{rl}^{i} \Gamma_{jk}^{r}\] This formula is complex and generate a lot of components for the curvature tensor. As the curvature has 4 indices, the number of curvature components is dim4
  • For a 2-dimension manifold (i.e., Plane, Sphere) 24=16 components
  • 3-dimension manifold  34= 81
  • For the space-time manifold 44= 256

In the case of the Levi-Civita connection, the properties of Torsion free and Metric compatibility produce symmetries that zero out many of these values. These symmetries are known as 1-2 indices symmetry, Bianchi identity, 3-4 indices symmetry and Flip symmetry [ref 9].

Note: All the 16 components of the curvature tensor in 2-dimension for a flat space are 0, as expected.

Geodesic deviation

A key challenge is determining whether geodesics in a curved space converge or diverge, a concept known as geodesic deviation, as depicted in the accompanying diagram.

Fig 1 Illustration of divergence of geodesics

The convergence and divergence of the geodesics is computed using the Riemann curvature tensor established in the previous section, more specifically the sign of the product of the Riemann curvature R and the vector orthogonal to the tangent vector on the geodesics.\[[R(X, Y)Y].X \ \ \ > 0 \ divergence, \ \ \ <0 \ convergence\]The geodesic deviation depends on the selection of the tangent vectors. This dependency is eliminated through normalization known as the sectional curvature tensor.

Sectional curvature tensor

The curvature tensor is undoubtedly a complex concept accompanied by intimidating formulas. However, a simpler expression of the curvature tensor can be obtained through sectional curvature. Essentially, sectional curvature refers to the curvature of two-dimensional sections of our manifold.

As outlined in the previous section, it is necessary to normalize the geodesic deviation to eliminate any dependency on the basis vectors of the manifold.
Given any point p ∈ M, any tangent 2-plane K to M at base point p, and a curvature tensor R, the sectional curvature of K(X, Y) is the real number: \[K(X, Y)= \frac{\left \langle R(X,Y)Y, X \right \rangle}{\left \langle X, X \right \rangle \left \langle Y, Y \right \rangle - \left \langle X, Y \right \rangle^{2}}\]Note: The sectional curvature tensor is a prerequisite to the computation of the Ricci tensor which will be described in a future article.

Implementation

As in the previous article on Riemann geometry, we rely on the Geomstats library [ref 3] to evaluate the curvature tensor.

Hypersphere curvature

For a hypersphere (unit sphere) the Riemann curvature tensor can be expressed as: \[R_{ijkl}=g_{ik}g_{jl} - g_{il}g_{jk} \ \ \ \  \ R_{ikl}^{j} = \frac{R_{ijkl}}{g_{ij}}\]After taking into account the symmetry properties of Levi-Civita connection, only one of the 16 components of the curvature tensor is not null. For the intrinsic coordinates u, v \[R_{212}^{1}=sin(u)^{2}\] The computation of the curvature tensor implemented by the method curvature_tensor added to the class RiemannianConnection. It leverages the Geomstats curvature method of the class RiemannianMetric.
The computation requires 3 tangent unit vectors and the base point on the manifold.

def curvature_tensor(self, tgt_vectors: List[np.array], base_pt: np.array) -> np.array:
  
  if len(tgt_vectors) != 3 or any(vec is None for vec in tgt_vectors):
     raise GeometricException(f'Tangent vectors for the curvature {str(tgt_vectors)} is not properly defined')

  return self.riemannian_metric.curvature(
            tgt_vectors[0], 
            tgt_vectors[1], 
            tgt_vectors[2], 
            base_pt)

Let's apply this method to the computation of the Riemann  curvature for a sphere. 

hypersphere = Hypersphere(dim=2, equip=True, default_coords_type='extrinsic')
riemann_connection = RiemannianConnection(hypersphere, 'HyperSphere')
base_pt = np.array([1.5, 2.0, 1.6])

X = np.array([0.4, 0.1, 0.8])
Y = np.array([0.5, 0.1, -0.2])
Z = np.array([0.4, 0.9, 0.0])

curvature = riemann_connection.curvature_tensor([X, Y, Z], base_pt)
print(f'Curvature: {curvature}')

Output:
Curvature: [ 0.009 -0.004 -0.282]

Sectional curvature of hypersphere

The method sectional_curvature_tensor added to the class RiemannianConnection implements the computation of the sectional curvature tensor that requires two tangent unit vectors and a base point on the sphere. It invokes the sectional_curvature method of the Geomstats class RiemannianMetric.

def sectional_curvature_tensor(self,
                              tgt_vec1: np.array,
                              tgt_vec2: np.array,
                              base_pt: Optional[np.array] = None) -> np.array:
 
   if len(tgt_vec1) != len(tgt_vec2):
      raise GeometricException(f'Dimension of tangent vectors for sectional curvature {len(tgt_vec1)} '
                                                f'and {len(tgt_vec2)} should be identical')

   return self.riemannian_metric.sectional_curvature(tgt_vec1, tgt_vec2, base_pt)

Our evaluation has 3 test cases for the tangent vectors:
  • V(x, y, z) & V(x, y, -z)
  • V(x, y, z) & -V(x, y, z)
  • V(x, y, z) & 2V(x, y, z)
hypersphere = Hypersphere(dim=2, equip=True, default_coords_type='intrinsic')
riemann_connection = RiemannianConnection(hypersphere, 'HyperSphere')
base_pt = np.array([1.5, 2.0, 1.6])

# Case 1 (X, Y=X with inverse on Z component)
X = np.array([0.4, 0.1, 0.8])
Y = np.array([0.4, 0.1, -0.8])
sec_curvature = riemann_connection.sectional_curvature_tensor(X, Y, base_pt)
print(f'Sectional curvature 1: {sec_curvature}')

# Case 2: (X, Y=-X)
X = np.array([0.4, 0.1, 0.8])
Y = np.array([-0.4, -0.1, -0.8])
sec_curvature = riemann_connection.sectional_curvature_tensor(X, Y, base_pt)
print(f'Sectional curvature 2: {sec_curvature}')

# Case 3:  (X, Y = 2.X)
X = np.array([0.4, 0.1, 0.8])
Y = np.array([0.8, 0.2, 1.6])
sec_curvature = riemann_connection.sectional_curvature_tensor(X, Y, base_pt)
print(f'Sectional curvature 3: {sec_curvature}')

Output
Sectional curvature: -1.0
Sectional curvature: 0.0
Sectional curvature: 0.0

References


-------------
Patrick Nicolas has over 25 years of experience in software and data engineering, architecture design and end-to-end deployment and support with extensive knowledge in machine learning. 
He has been director of data engineering at Aideo Technologies since 2017 and he is the author of "Scala for Machine Learning", Packt Publishing ISBN 978-1-78712-238-3 and Geometric Learning in Python Newsletter on LinkedIn.

Appendix

Here is the list of published articles related to geometric learning:

Keywords: geometriclearning, riemanngeometry, manifold, differential geometry, ai, python, geomstats




Posted by No comments:
Labels: , , , , , , , ,

Friday, April 12, 2024

Riemann Metric & Connection for Geometric Learning

Target audience: Advanced
Estimated reading time: 7'

Concepts from Riemannian geometry are fundamental in various fields such as Robotics, Computer Vision, Machine Learning, Quantitative Finance, and Medical Imaging. 
This article utilizes the Geomstats library to explore and experiment with metrics, connections, and parallel transport in these contexts.


Table of contents

       Vector fields
       Setup
       Riemann metric
Follow me on LinkedIn

What you will learn: How to implement Riemann metric and connection in Python using Geomstats library.


Notes: 
  • Environments: Python  3.10.10, Geomstats 2.7.0
  • This article assumes that the reader is somewhat familiar with differential and tensor calculus [ref 1, 2]. Please refer to the previous articles related to geometric learning [ref 5, 6, 7].
  • Source code is available at  Github.com/patnicolas/Data_Exploration/manifolds
  • To enhance the readability of the algorithm implementations, we have omitted non-essential code elements like error checking, comments, exceptions, validation of class and method arguments, scoping qualifiers, and import statements.

Introduction

Geometric learning addresses the difficulties of limited data, high-dimensional spaces, and the need for independent representations in the development of sophisticated machine learning models, including graph-based and physics-informed neural networks.

The primary goal of learning Riemannian geometry is to understand and analyze the properties of curved spaces that cannot be described adequately using Euclidean geometry alone. Riemannian geometry enables to describe the geometric structures of manifolds equipped with a metric, which defines the concept of distance and angle on these spaces.

This article is the seventh part of our ongoing series focused on geometric learning. In this installment, we utilize the Geomstats Python library [ref. 3] and explore the hypersphere manifold, which was introduced in a previous piece, Geometric Learning in Python: Manifolds - Hypersphere  and is detailed in the Geomstats API [ref. 4]. 

I highly recommend watching the comprehensive series of 22 YouTube videos Tensor Calculus - Eigenchris  to familiarize yourself with fundamental concepts of differential geometry. 
Summaries of my earlier articles on this topic can be found in the Appendix

Riemann geometry

First, let's review some basic concepts in Riemannian geometry that we described in previous articles [ref 5, 6, 7].

Vector fields

In tensor calculus a vector field is an assignment of a vector X to each point p on a Euclidean space. For instance, a vector field on a 3D sphere can be visualized as a collection of arrows with a given direction and length.

Fig. 1 Visualization of vector fields on a sphere in 3D Euclidean space

Formally, a vector field X(p) at point p is defined as: \[X: S \rightarrow \mathbb{R}^{n} \ \ \ \ \ X(p) = \sum_{i=1}^{n} X_{i}(p)\frac{\partial }{\partial x_{i}}\]
The gradient of a vector field on a differential manifold is defined as: \[\triangledown _{\vec{v}}f= \triangledown f. \vec{v}=\sum_{i=1}^{n}v^{i}\frac{\partial f}{\partial x^{i}} \ \ \ (1) \]

Covariant derivative

The gradient of a vector field is defined by its components along the basis vectors, which are orthogonal in Euclidean space. The covariant derivative, on the other hand, calculates a gradient along a specific vector, particularly the tangent vector of a manifold.

The covariant derivative measures the rate of change of vector fields, considering the variation of basis vectors. For instance, in flat space, the covariant derivative of a vector field corresponds to the ordinary derivative, calculating both the derivative of the vector components and the derivative of the basis vectors.
The covariant derivative can be expressed as a directional derivative in extrinsic coordinates as \[\triangledown_{\frac{\partial }{\partial u^{i}}}\vec{v}=\frac{\partial \vec{v}}{\partial u^{i}}-\vec{n} \ \ \ \ (2)\]
Fig. 2 Visualization of covariant derivative on a sphere from formula (2)


Parallel transport

The covariant derivative is often conceptualized as an affine connection because it links the tangent spaces (or planes) at two points on the manifold. Parallel transport involves moving a geometric vector or tensor along smooth curves within the manifold. This process is facilitated by the covariant derivative, which enables the transportation of manifold vectors along these curves, ensuring they remain parallel according to the connection.

Fig. 3 Illustration of parallel transport along smooth curves on a sphere

In a parallel transport the rate of change of the tangent vector as illustrated in the previous figure is 0. Therefore the covariant derivative along the tangent vector is 0.
\[\triangledown _{\frac{\partial }{\partial u^{i}}}\vec{v} = 0 \ \ \  (3)\] 
Fig. 4 Illustration of Christoffel symbols for the covariant derivative on a sphere

The covariant derivative can be expressed using the Christoffel symbols as\[\triangledown _{\frac{\partial }{\partial u^{i}}}\vec{v} =\left \lfloor \frac{\partial v^{k}}{\partial u^{i}} +v^{j}.\Gamma^{k}_{ij}\right \rfloor \vec{e_{k}} \ \ \ (4)\]

Metric & connection

A metric tensor is a structure (such as a tensor or matrix) on a manifold which facilitates the definition of distances, angles, and norms, in a manner comparable to the inner product in Euclidean space. Given a set of basis vector ei  the metric tensor gij and the inverse metric tensor are defined as \[g_{ij}=\vec{e_{i}}.\vec{e_{j}}=\left \langle \frac{\partial }{\partial x^{i}} \frac{\partial }{\partial x^{j}} \right \rangle \ \ (5) \ \ \ Inverse \ metric \ \mathfrak{g^{ij}} \ \ g_{ij}.\mathfrak{g^{ij}} =\delta _{ij}\]

An affine connection that incorporates a metric is described as a bilinear map. Given a manifold M, a tangent spaceTM and vector fields space (TM). \[\begin{matrix} \Gamma(T_{M})\ X \ \Gamma(T_{M})\rightarrow \Gamma(T_{M}) \\ \triangledown : \ \ \ (X, Y) \rightarrow \triangledown_{X}Y \ \ \ \ \ \ \ \ \ \end{matrix}\]For any infinitely differentiable and continuous function f, \[\triangledown_{fX}Y=f\triangledown_{X}Y \ \ \ \ \triangledown_{X}(fY)=\frac{\partial f}{\partial X}Y+\triangledown_{fX}Y\]An affine connection that preserves the Riemann metric under parallel transport and is symmetric (commutative covariant derivative) is known as the Levi-Civita connection.

Levi-Civita coefficients

The Levi-Civita connection is a type of metric connection that is compatible with the metric tensor of a Riemannian manifold. Its covariant derivative of the metric tensor with respect to the Levi-Civita connection is zero. The property of parallel transport ensures that the connection preserves the inner product structure defined by the metric tensor.
It is related to the fundamental Theorem of riemannian geometry which states "There is a unique connection (or covariant derivative) that is Torsion free and has metric compatibility"
The Levi-Civita connection is essential is the computation of geodesics, curvature tensor and parallel transport on Riemannian manifolds.

The coefficients of the Levi-Civita used to compute the covariant derivative is defined by the Christoffel symbols. Using the intrinsic coordinates, the Christoffel symbols are computed from the metric tensors g as \[\Gamma _{jk}^{m}=\frac{1}{2}\mathfrak{g}^{im}\left ( \frac{\partial }{\partial u^{k}} .g_{ij} +\frac{\partial }{\partial u^{j}} .g_{ki} - \frac{\partial }{\partial u^{i}} .g_{jk} \right ) \ \ \ (6) \]


Use case: hypersphere

We will illustrate the implementation of Riemann metric and Levi-Civita connection on the hypersphere space we introduced in a previous article Geometric Learning in Python: Manifolds

Setup

As a reminder, we reuse the structure of data points on a manifold defined by the class ManifoldPoint described in a previous post ManifoldPoint definition. 

@dataclass
class ManifoldPoint:
    id: AnyStr
    location: np.array
    tgt_vector: Optional[List[float]] = None
    geodesic: Optional[bool] = False

As in previous articles, we utilize the Geomstats Python library to explore, implement, and assess the Riemannian manifold and its connection. For ease of use, we encapsulated the Riemann metric within a class named RiemannianConnection as follows:

import geomstats.backend as gs
from geomstats.geometry.riemannian_metric import RiemannianMetric
from geomstats.geometry.base import LevelSet
from typing import AnyStr, Optional
from manifoldpoint import ManifoldPoint
from geometricexception import GeometricException


class RiemannianConnection(object):

    def __init__(self, space: LevelSet, manifold_type: AnyStr):
        self.riemannian_metric = RiemannianConnection.__get_metric(space)

        self.manifold_descriptor = f'{manifold_type}\nDimension: {space.dim}\nShape: {space.shape}' \
                                   f'\nCoordinates type: {space.default_coords_type}'
        # Add 1 dimension for extrinsic coordinates
        self.ndim = space.dim+1 if space.default_coords_type == "extrinsic" else space.dim

The Riemann metric is inferred from the type of manifold (or level set), space. The variable ndim captures the dimension of tangent vectors in either intrinsic or extrinsic coordinates.


Riemann metric

Let's implement the inner product of two tangent vectors, tgt_vec1 and tgt_vec2 at a base point base_pt defined in formula (4) for the RiemannianConnection class. 

def inner_product(self, tgt_vec1: np.array, tgt_vec2: np.array, base_pt: np.array) -> np.array:
  if len(tgt_vec1) != len(tgt_vec2):
      raise GeometricException(f'Inner product of vector size {len(tgt_vec1)} and vector size {len(tgt_vec1)}')

  return 0 if len(tgt_vec1) == 0 or len(tgt_vec2) == 0 \
            else self.riemannian_metric.inner_product(tgt_vec1, tgt_vec2, base_pt)

The inner product on the hypersphere is computed for a vector and its negative value.

dim = 2
coordinates_type = 'extrinsic'

hypersphere = Hypersphere(dim=dim, equip=True, default_coords_type=coordinates_type)
riemann_connection = RiemannianConnection(hypersphere, 'hypersphere')

vector = np.array([0.4, 0.1, 0.8])
base_point = np.array([0.4, 0.1, 0.0])

inner_product = riemann_connection.inner_product(vector, -vector, base_point   #. -0.81

Output:
Inner product: -0.81


Riemann connection

Parallel transport
Let's implement the formula (5) by adding the method parallel_transport to the class RiemannianConnection. Its implementation invoked the Geomstats method, RiemannianMetric.parallel_transport. The parallel transport is initialized from the manifold base point, manifold_base_pt, to an optional end point, manifold_end_pt. The direction for the parallel transport can be optionally specified with the tangent vector at the base point as default value.

def parallel_transport(
   self,
   manifold_base_pt: ManifoldPoint,
   manifold_end_pt: Optional[ManifoldPoint] = None,
   direction: Optional[np.array] = None) -> np.array:
        
   if manifold_base_pt.ndim() != self.ndim:
      raise GeometricException(f'Base pt dimension {manifold_base_pt.ndim()} should be {self.ndim}')

   return self.riemannian_metric.parallel_transport(
       manifold_base_pt.tgt_vector,
       manifold_base_pt.location,
       direction,
       manifold_end_pt.location)

The evaluation consists of one manifold base point with a tangent vector and another manifold end point with a translation +0.4

vector = np.array([0.5, 0.1, 2.1])
base_pt = np.array([1.5, 2.0, 1.6])

manifold_base_pt = ManifoldPoint(id='base', location=base_pt, tgt_vector=vector
manifold_end_pt = ManifoldPoint(id='end',location=base_pt + 0.4)

parallel_trans = riemannian_connection.parallel_transport(manifold_base_pt, manifold_end_pt)
print(f'Parallel transport: {parallel_trans}')

Output:
Parallel transport: [[-1.20097215 -2.09879718  0.29946284]]


Christoffel symbols 
Lastly, let's compute the Christoffel symbols for the hypersphere as a given point of intrinsic coordinates (u, v).\[\begin{bmatrix} \Gamma _{11}^{1} & \Gamma _{12}^{1}\\ \Gamma _{11}^{1} & \Gamma _{22}^{1} \end{bmatrix} = \begin{bmatrix} 0 & 0\\ 0 & - \frac{1}{2}sin\left ( 2u \right ) \end{bmatrix} \begin{bmatrix} \Gamma _{11}^{2} & \Gamma _{12}^{2}\\ \Gamma _{11}^{2} & \Gamma _{22}^{2} \end{bmatrix} = \begin{bmatrix} cot(u) & 0\\ 0 & cot(u) \end{bmatrix} \ (7) \]We add the method levi_civita_coefficients to the class RiemannianConnection. Its implementation of this formula leverages the Geomstats, RiemannianMetric.christoffels method.

def levi_civita_coefficients(self, base_pt: np.array) -> np.array:
   if len(base_pt) != self.ndim:
       raise GeometricException(f'Base pt dimension {len(base_pt)} should be {self.ndim}')
   
   return self.riemannian_metric.christoffels(base_pt)

Let's validate the formula (6) for the Christoffel symbols on the Hypersphere with u = 1.5 and v = 2.0

base_pt = np.array([1.5, 2.0, 1.6])
levi_civita_coefs = riemann_connection.levi_civita_coefficients(base_pt)
print(f'Levi-Civita coefficients:\n{str(levi_civita_coefs)}')

Output:
Levi-Civita coefficients:
[[[ 0.          0.        ]
  [ 0.         -0.07056   ]]

 [[ 0.          0.07091484]
  [ 0.07091484  0.        ]]]

We validate that 0.5*sin(2u) = 0.0756 and cot(u) = 0.0709148

References




-------------
Patrick Nicolas has over 25 years of experience in software and data engineering, architecture design and end-to-end deployment and support with extensive knowledge in machine learning. 
He has been director of data engineering at Aideo Technologies since 2017 and he is the author of "Scala for Machine Learning", Packt Publishing ISBN 978-1-78712-238-3 and Geometric Learning in Python Newsletter on LinkedIn.


Appendix


Keywords: geometriclearning, riemanngeometry, manifold, differential geometry, ai, python, geomstats
Posted by No comments:
Labels: , , , , , ,
Subscribe to: Posts (Atom)