I delve into a diverse range of topics, spanning programming languages, machine learning, data engineering tools, and DevOps. Our articles are enriched with practical code examples, ensuring their applicability in real-world scenarios.
Showing posts with label Deep learning. Show all posts
Showing posts with label Deep learning. Show all posts
Thursday, May 2, 2024
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
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
[5] Geometric Learning in Python: Basics
[6] Geometric Learning in Python: Manifolds
[7] Geometric Learning in Python: Vector and Covector fields
[6] Geometric Learning in Python: Manifolds
[7] Geometric Learning in Python: Vector and Covector fields
[8] Geomstats API
-------------
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
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
Appendix
Here is the list of published articles related to geometric learning:
Geometric Learning in Python: Basics introduced differential geometry as an applied to machine learning and its basic components.
Geometric Learning in Python: Manifolds described manifold components such as tangent vectors, geodesics with implementation in Python for Hypersphere using the Geomstats library.
Geometric Learning in Python: Intrinsic Representation reviewed the various coordinates system using extrinsic and intrinsic representation.
Geometric Learning in Python: Vector and Covector fields described vector and covector fields with Python implementation in 2 and 3-dimension spaces.
Geometric Learning in Python: Vector Operators illustrated the differential operators, gradient, divergence, curl and laplacian using SymPy library.
Geometric Learning in Python: Functional Data Analysis described the key elements of non-linear functional data analysis to analysis curves, images, or functions in very high-dimensional spaces
Geometric Learning in Python: Manifolds described manifold components such as tangent vectors, geodesics with implementation in Python for Hypersphere using the Geomstats library.
Geometric Learning in Python: Intrinsic Representation reviewed the various coordinates system using extrinsic and intrinsic representation.
Geometric Learning in Python: Vector and Covector fields described vector and covector fields with Python implementation in 2 and 3-dimension spaces.
Geometric Learning in Python: Vector Operators illustrated the differential operators, gradient, divergence, curl and laplacian using SymPy library.
Geometric Learning in Python: Functional Data Analysis described the key elements of non-linear functional data analysis to analysis curves, images, or functions in very high-dimensional spaces
Geometric Learning in Python: Riemann Metric & Connection reviewed Riemannian metric tensor, Levi-Civita connection and parallel transport for hypersphere.
#geometriclearning #riemanngeometry #manifold #differential geometry #ai #python #geomstats
Wednesday, April 3, 2024
Geometric Learning in Python: Vector Operators
Target audience: Beginner
Estimated reading time: 5'
Physics-Informed Neural Networks (PINNs) are gaining popularity for integrating physical laws into deep learning models. Essential tools like vector operators, including the gradient, divergence, curl, and Laplacian, are crucial in applying these constraints.
Table of contents
What you will learn: How to implement the vector gradient, divergence, curl and laplacian operators in Python using SymPy library.
Notes:
- Environments: Python 3.10.10, SymPy 1.12, Matplotlib 3.8.2
- This article assumes that the reader is familiar with differential and tensor calculus.
- Source code is available at github.com/patnicolas/Data_Exploration/diffgeometry
- 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.
Note: This article is the 5th installment in our series on Geometric Learning in Python following
- Geometric Learning in Python: Basics introduces differential geometry as an applied to machine learning and its basic components.
- Geometric Learning in Python: Manifolds describes manifold components such as tangent vectors, geodesics with implementation in Python for Hypersphere using the Geomstats library.
- Geometric Learning in Python: Intrinsic Representation Reviews the various coordinates system using extrinsic and intrinsic representation.
- Geometric Learning in Python: Vector and Covector fields describes vector and covector fields with Python implementation in 2 and 3-dimension spaces.
The following description of vector differential operators leverages some of concept defined in the previous articles in the Geometric Learning in Python series and SymPy library.
SymPy is a Python library dedicated to symbolic mathematics. Its implementation is as simple as possible to be comprehensible and easily extensible with support for differential and integral calculus, Matrix operations, algebraic and polynomial equations, differential geometry, probability distributions, 3D plotting. The source code is available at github.com/sympy/sympy.git
In tensor calculus, a vector operator is a type of differential operator:
- The gradient transforms a scalar field into a vector field.
- The divergence changes a vector field into a scalar field.
- The curl converts a vector field into another vector field.
- The laplacian takes a scalar field and yields another scalar field.
Gradient
Consider a scalar field f
We create a class, VectorOperators, that wraps the following operators: gradient, divergence, curl and laplacian. The constructor initializes the function as an expression for which the various operators are applied to.
class VectorOperators(object):
def __init__(self, expr: Expr): # Expression for the input function
self.expr = expr
def gradient(self) -> VectorZero:
from sympy.vector import gradient
return gradient(self.expr, doit=True)Let's calculate the gradient vector for the function \[f(x,y,z)=x^{2}+y^{2}+z^{2} \] , as depicted in the plot below.
Fig. 1 3D visualization of function f(x,y,z) = x**2 + y**2 + z**2
The Python for the visualization of the function is described in Appendix (Visualization function)
The following code snippet computes the gradient of this function as \[\triangledown f(x,y,z) = 2x.\vec{i} + 2y.\vec{j} + 2y.\vec{k}\]
r = CoordSys3D('r')f = r.x*r.x+r.y*r.y+r.z*r.zvector_operator = VectorOperators(f) grad_f = vector_operator.gradient() # 2*r.x + 2*r.y+ 2*r.z
The function f is defined using the default Euclidean coordinates r. The gradient is depicted in the following plot implemented using Matplotlib module. The actual implementation is described in Appendix (Visualization gradient)
Fig. 2 3D visualization of grad_f(x,y,z) = 2x.i + 2y.j + 2z.k
Divergence
Divergence is a vector operator used to quantify the strength of a vector field's source or sink at a specific point, producing a signed scalar value. When applied to a vector F, comprising components X, Y, and Z, the divergence operator consistently yields a scalar result [ref 2]. \[div(F)=\triangledown .F=\frac{\partial X}{\partial x}+\frac{\partial Y}{\partial y}+\frac{\partial Z}{\partial z}\]
Let's implement the computation of the divergence as method of the class VectorOperators:
def divergence(self, base_vec: Expr) -> VectorZero:
from sympy.vector import divergence
div_vec = self.expr*base_vec
return divergence(div_vec, doit=True)Using the same instance vector_operators as with the gradient calculation:
divergence = vector_operators.divergence(r.i + r.j + r.k)
The execution of the code above produces the following:\[div(f(x,y,z)[\vec{i} + \vec{j}+ \vec{k}]) = 2x(y+z+xy)\]
Curl
In mathematics, the curl operator represents the minute rotational movement of a vector in three-dimensional space. This rotation's direction follows the right-hand rule (aligned with the axis of rotation), while its magnitude is defined by the extent of the rotation [ref 2]. Within a 3D Cartesian system, for a three-dimensional vector F, the curl operator is defined as follows: \[ \triangledown * \mathbf{F}=\left (\frac{\partial F_{z}}{\partial y}- \frac{\partial F_{y}}{\partial z} \right ).\vec{i} + \left (\frac{\partial F_{x}}{\partial z}- \frac{\partial F_{z}}{\partial x} \right ).\vec{j} + \left (\frac{\partial F_{y}}{\partial x}- \frac{\partial F_{x}}{\partial y} \right ).\vec{k} \]. Let's add the curl method to our class VectorOperator as follow:
def curl(self, base_vectors: Expr) -> VectorZero:
from sympy.vector import curl
curl_vec = self.expr*base_vectors
return curl(curl_vec, doit=True)
For the sake of simplicity let's compute the curl along the two base vectors j and k:
curl_f = vector_operator.curl(r.j+r.k) # j and k direction only
The execution of the curl method outputs: \[curl(f(x,y,z)[\vec{j} + \vec{k}]) = x^{2}\left ( z-y \right).\vec{i}-2xyz.\vec{j}+2xyz.\vec{k}\]
Laplacian
In mathematics, the Laplacian, or Laplace operator, is a differential operator derived from the divergence of the gradient of a scalar function in Euclidean space [ref 3]. The Laplacian is utilized in various fields, including calculating gravitational potentials, solving heat and wave equations, and processing images.
It is a second-order differential operator in n-dimensional Euclidean space and is defined as follows: \[\triangle f= \triangledown ^{2}f=\sum_{i=1}^{n}\frac{\partial^2 f}{\partial x_{i}^2}\]. The implementation of the method laplacian reflects the operator is the divergence (step 2) of the gradient (Step 1)
def laplacian(self) -> VectorZero:
from sympy.vector import divergence, gradient
# Step 1 Compute the Gradient vector
grad_f = self.gradient()
# Step 2 Apply the divergence to the gradient
return divergence(grad_f)
Once again we leverage the 3D coordinate system defined in SymPy to specify the two functions for which the laplacian has to be evaluated
r = CoordSys3D('r')
f = r.x*r.x + r.y*r.y + r.z*r.z
vector_operators = VectorOperators(f)
laplace_op = vector_operators.laplacian()
print(laplace_op) # 6
f = r.x*r.x*r.y*r.y*r.z*r.z
vector_operators = VectorOperators(f)
laplace_op = vector_operators.laplacian()
print(laplace_op) # 2*r.x**2*r.y**2 + 2*r.x**2*r.z**2 + 2*r.y**2*r.z**2
Outputs: \[\triangle (x^{2} + y^{2}+z^{2}) = 6\] \[\triangle (x^{2}y^{2}z^{2})=2(x^{2}y^2 +x^{2}z^{2}+y^{2}z^{2}) \]
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
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
Appendix
Visualization function
The implementation relies on
- Numpy meshgrid to setup the axis and values for the x, y and z axes.
- Matplotlib scatter method to display the values generated by the function f
The function to plot takes x, y and z grid values and to generate data f(x,y,z)
@staticmethod
def show_3D_function(f: Callable[[float, float, float], float], grid_values: np.array) -> NoReturn:
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
# Setup the grid with the appropriate units and boundary
x, y, z = np.meshgrid(grid_values, grid_values, grid_values)
# Apply the function f
data = f(x, y, z)
# Set up plot (labels, legend,..)
ax: Axes = self.__setup_3D_plot('3D Plot f(x,y,z) = x^2 + y^2 + z^2')
# Display the data along x, y and z using scatter plot
ax.scatter(x, y, z, c=data)
plt.show()
Visualization gradient
The 3 components of the gradient vector as passed as argument, grad_f.
The implementation relies on
- Numpy meshgrid to setup the axis and values for the x, y and z axes.
- Matplotlib quiver method to display the gradient vectors at each grid values
@staticmethod
def show_3D_gradient(grad_f: List[Callable[[float], float]], grid_values: np.array) -> NoReturn:
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
# Setup the grid with the appropriate units and boundary
x, y, z = np.meshgrid(grid_values, grid_values, grid_values)
ax = self.__setup_3D_plot('3D Plot Gradient 2x.i + 2y.j + 2z.k')
# Extract the gradient df/dx, df/dy and df/dz
X = grad_f[0](x)
Y = grad_f[1](y)
Z = grad_f[2](z)
# Display the gradient vectors as vector fields
ax.quiver(x, y, z, X, Y, Z, length=1.5, color='grey', normalize=True)
plt.show()
Subscribe to:
Posts (Atom)