Wednesday, September 18, 2024

Posts History

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.

Follow me on LinkedIn

2024

Monday, September 16, 2024

Operations on SO3 Lie Groups in Python

Target audience: Advanced
Estimated reading time: 5'
Lie groups play a crucial role in Geometric Deep Learning by modeling symmetries such as rotation, translation, and scaling. This enables non-linear models to generalize effectively for tasks like object detection and transformations in generative models.


Table of contents
       Lie manifolds
       Why Lie groups
       Implementation
       Inverse rotation
       Composition
Follow me on LinkedIn

What you will learn:  How to implement and evaluation key operations on Lie Special Orthogonal Group of dimension 3.

Notes

  • This post is a follow up on articles related to differential geometry and geometry learning [ref 1, 2, 3 and 4]
  • Environments: Python 3.11,  Matplotlib 3.9, Geomstats 2.8.0
  • Source code is available at  Github.com/patnicolas/Data_Exploration/Lie
  • 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 statement.

Disclaimer : A thorough tutorial and explanation of Lie groups, Lie algebras, and geometric priors for deep learning models is beyond the scope of this article. Instead, the following sections concentrate on experiments involving key elements and operations on Lie groups using the Geomstats Python library [ref 5].
 

Overview

Let's start with some basic definitions

Lie manifolds

A smooth manifold is a topological space that locally resembles Euclidean space and allows for smooth (infinitely differentiable) transitions between local coordinate systems. This structure allows for the use of calculus on the manifold. 

The tangent space at a point on a manifold is the set of tangent vectors at that point, like a line tangent to a circle or a plane tangent to a surface.
Tangent vectors can act as directional derivatives, where you can apply specific formulas to characterize these derivatives.

Fig. 1 Manifold with tangent space and exponential/logarithm maps

In differential geometry, a Lie group is a mathematical structure that combines the properties of both a group and a smooth manifold. It allows for the application of both algebraic and geometric techniques. As a group, it has an operation (like multiplication) that satisfies certain axioms (closure, associativity, identity, and invertibility) [ref 6].
A 'real' Lie group is a set G with two structures: G is a group and G is a (smooth, real) manifold. These structures agree in the following sense: multiplication (a.k.a. product or composition) and inversion are smooth maps.
A morphism of Lie groups is a smooth map which also preserves the group operation: f(gh) = f(g)f(h) and f(1) = 1.
Fig. 2 Manifold with tangent space and identity and group element 
(Courtesy A. Kirillov Jr Department of Mathematics SUNY at Stony Brook)

Why Lie groups

Lie groups have numerous practical applications in various fields:
  • Physics: They describe symmetries in classical mechanics, quantum mechanics, and relativity. 
  • Robotics: Lie groups model the motion of robots, particularly in the context of rotation and translation (using groups like SO(3) and SE(3)).
  • Control Theory: Lie groups are used in the analysis and design of control systems, especially in systems with rotational or symmetrical behavior.
  • Computer Vision: They help in image processing and 3D vision, especially in tasks involving rotations and transformations.
  • Differential Equations: Lie groups are instrumental in solving differential equations by leveraging symmetry properties.
In term of machine learning...
  • Geometric Deep Learning: Lie groups help capture rotational, translational, or scaling symmetries, making models more efficient and generalizable in tasks like image recognition and 3D object detection.
  • Generative Models: Lie groups allow for more structured latent spaces, enabling better control over transformations in generative models like GANs and VAEs.
  • Reinforcement Learning: They are used to model continuous action spaces, improving control over robotic systems.
  • Optimization: Lie groups help design efficient optimization techniques on curved spaces, like Riemannian manifolds.

Examples of Lie groups

Group of Invertible 2 x 2 matrices of real values
equation

Group of Invertible 2 x 2 matrices of complex values
equation

Group of Invertible 3 x 3 matrices of real values
equation

Special Unitary group 2 x 2 matrices with determinant 1
equation

Special Unitary group 3 x 3 matrices with determinant 1
equation

Special Orthogonal (3D rotation, 2x2 matrices) group
equation

Special Orthogonal (3D rotation, 3x3 matrices) group
equation

Special Euclidean group
equation


Special Orthogonal Group

The Special Orthogonal Group in 3 dimensions, SO(3) is the group of all rotation matrices in 3 spatial dimensions.
It can be defined by 3 rotation elements for each of the axis of rotation x, y, and z.

equation

equation

equation

Implementation

Geomstats is a free, open-source Python library designed for conducting machine learning on data situated on nonlinear manifolds, an area known as Geometric Learning. This library offers object-oriented, thoroughly unit-tested features for fundamental manifolds, operations, and learning algorithms, compatible with various execution environments, including NumPyPyTorch, and TensorFlow (Overview Geomstats library).

The library is structured into two principal components:
  • geometry: This part provides an object-oriented framework for crucial concepts in differential geometry, such as exponential and logarithm maps, parallel transport, tangent vectors, geodesics, and Riemannian metrics.
  • learning: This section includes statistics and machine learning algorithms tailored for manifold data, building upon the scikit-learn framework.


First let's wrap the element of a point or matrix into a class, SO3Point with the following attributes:
  • group_element Point on the manifold
  • base_point on the manifold (Identity if undefined)
  • description which describes the rotation matrix
from dataclasses import dataclass

@dataclass
class SO3Point:
    group_element: np.array
    base_point: np.array                # Default identity matrix
    descriptor: AnyStr

Secondly, let's build a class ,LieSO3Group, that encapsulates the definition of the Special orthogonal group of dimension 3 and its related operations.
We specify two constructors:
  • __init__: Default constructor that create a new element in the SO3 manifold, group_element, using a tangent vector (3x3 rotation matrix) tgt_vector and identity, base_point.
  • build: Alternative constructor for which the tangent vector is a 9- element list and a tuple as base point.
import geomstats.backend as gs
from geomstats.geometry.special_orthogonal import SpecialOrthogonal


class LieSO3Group(object):
    dim = 3
    # Lie group as defined in Geomstats library
    lie_group =  SpecialOrthogonal(n=dim, point_type='vector', equip=False)
    identity = gs.eye(dim)     # Define identity


    def __init__(self, tgt_vector: np.array, base_point: np.array = identity) -> None
        self.tangent_vec = gs.array(tgt_vector)

        # Exp. a left-invariant vector field from a base point
        self.group_element = LieSO3Group.lie_group.exp(self.tangent_vec, base_point)
        self.base_point = base_point

    @classmethod
    def build(cls, tgt_vector: List[float], base_point: List[float] = None) -> Self:
        np_input = np.reshape(tgt_vector, (3, 3))
        np_point = np.reshape(base_point, (3, 3) if base_piont is not None 
                          else LieSo3Group.identity

        return cls(tgt_vector=np_input, base_point=np_point)


We use the constructors of the LieSO3Group class to generate two points on the SO3 manifold, so3_point1 and so3_point2. These points are represented as vectors in 3-dimensional Euclidean space for visualization purposes.
The first point use the identity as base point for the tangent space, while the second point use the coordinates:
[[-0.4, 0.5, -0.3], 
 [0.0, 1.0, -0.6],
 [0.2, 0.1, 0.0]]

so3_tangent_vec = [0.4, 0.3, 0.8, 0.2, 0.4, 0.1, 0.1, 0.2, 0.6]

# First vector set
so3_group = LieSO3Group.build(so3_tangent_vec)
so3_point1 = SO3Point(
     group_element=so3_group.group_element,
     base_point=LieSO3Group.identity,
     descriptor='3D Rotatio\n[0.4 0.3 0.8]\n[0.2 0.4 0.1]\n[0.1 0.2 0.6]')

# Second vectors set
base_point = [-0.4, 0.5, -0.3, 0.0, 1.0, -0.6, 0.2, 0.1, 0.0]
so3_group2 = LieSO3Group.build(so3_tangent_vec, base_point)
so3_point2 = SO3Point(
     group_element=so3_group2.group_element,
     base_point=base_point,
     descriptor='3D Rotation\n[-0.4, 0.5, -0.3]\n[0.0, 1.0, -0.6]\n[0.2, 0.1, 0.0]')

LieSO3Group.visualize_all([so3_point1, so3_point2])

The code for the visualization methods are available on Github at LieSO3Group

Fig. 3 Visualization of SO3 group element at identity vs. arbitrary base point

Inverse rotation

The first operation to illustrate is the inverse rotation using Geomstats library. It is implemented by the LieSO3Group method inverse.

def inverse(self) -> Self:
    inverse_group_element = LieSO3Group.lie_group.inverse(self.group_element)
    return LieSO3Group(inverse_group_element)

For simplicity sake, we create a SO3 group element at identity (base point 0, 0, 0) and visualize its inverse on the same base point.

 # Define the origin point on SO3 manifold given a 3x3-9 elements tangent space 
so3_tangent_vec = [0.4, 0.3, 0.8, 0.2, 0.4, 0.1, 0.1, 0.2, 0.6]
so3_group = LieSO3Group.build(so3_tangent_vec)
so3_point = SO3Point(
     group_element=so3_group.group_element,
     base_point=LieSO3Group.identity,
     descriptor='3D Rotation:\n[0.4 0.3 0.8]\n[0.2 0.4 0.1]\n[0.1 0.2 0.6]')
       
 # Inverse SO3 rotation matrix
so3_inv_group = so3_group.inverse()
inv_so3_point = SO3Point(
     group_element=so3_inv_group.group_element,
     base_point=LieSO3Group.identity,
     descriptor='3D Inverse Rotation:\n[-0.4 -0.3 -0.8]\n[-0.2 -0.4 -0.1]\n[-0.1 -0.2 -0.6]')

As expected, the inverse of an SO3 group element represents the inverse rotation in 3-dimensional space. This would not hold true if the base point were not the identity element.
Fig. 4 Visualization of the inverse of SO3 element at identity

Composition

As previously mentioned, a smooth manifold with the structure of a Lie group ensures that the product (or composition) of two elements also belongs to the manifold. The following LieSO3Group method, product, utilizes the `compose method of Geomstats library.

def product(self, lie_so3_group: Self) -> Self:
     composed_group_point = LieSO3Group.lie_group.compose(self.group_element lie_so3_group.group_element)
     return LieSO3Group(composed_group_element)


In this initial test, we use the inverse of the first element on the SO3 manifold as the second component for the composition.

# Second SO3 rotation matrix
so3_tangent_vec2 = [-x for x in so3_tangent_vec]
so3_group2 = LieSO3Group.build(tgt_vector=so3_tangent_vec2)
so3_point2 = SO3Point(
      group_element=so3_group2.group_element,
      base_point=LieSO3Group.identity
       descriptor='3D Rotation:\n[-0.4 -0.3 -0.8]\n[-0.2 -0.4 -0.1]\n[-0.1 -0.2 -0.6]')

so3_group_product = so3_group.product(so3_group2)
    

As expected, the composition of an element in the SO3 group with its inverse results in the null  rotation matrix.

Fig. 5 Visualization of the composition of a SO3 group element with its inverse at identity


The second evaluation consists of composing the initial SO3 group element with a 3x3 rotation matrix with identical elements.

Fig. 6  Visualization of the composition of two SO3 group elements at identity

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.









Tuesday, September 3, 2024

Fractal Dimension of Objects in Python

Target audience: Beginner
Estimated reading time: 5'

Are you finding it challenging to configure a convolutional neural network (CNN) for modeling 3D objects?

The complexity of a 3D object can pose a considerable challenge when tuning the parameters of a 3D convolutional neural network. Fractal analysis [ref 1] offers a way to measure the complexity of key features, volumes, and boundaries within the object, providing valuable insights that can help data scientists fine-tune their models for better performance.

Table of contents
Follow me on LinkedIn

What you will learn: How to evaluate the complexity of a 3-dimension object using fractal dimension index.

Notes

  • This article is a follow up on Fractal Dimension of Images in Python
  • Environments: Python  3.11,  Matplotlib 3.9.
  • Source code is available at  Github.com/patnicolas/Data_Exploration/fractal
  • 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

As described in a previous article, Fractal Dimension of Images in Python - Overview , a fractal dimension is a measure used to describe the complexity of fractal patterns or sets by quantifying the ratio of change in detail relative to the change in scale [ref 2].
Among the various approaches to estimate the fractal dimension, from variation, structure function methods, root mean square and R/S analysis, we selected the box counting method because of its simplicity and visualization capability.

Point cloud
To evaluate our method for calculating the fractal dimension index of an object, it's necessary to simulate or represent a 3D object. This is achieved by generating a cluster of random data points across the x, y, and z axes, commonly referred to as a point cloud.

Box counting method
The box counting method [ref 3] is described in a previous article, Box-counting method

If N is the number of measurements units for our counting boxes or cubesand eps the related scaling factor, the fractal dimension index is computed as\[ D=- \displaystyle \lim_{\epsilon  \to 0} \frac{log(N)}{log(\epsilon)} \simeq - \frac{log(N)}{log(eps)} \ \ with \ N=r^3  \]We will use the height of the square box r as our measurement unit for images.
The fractal dimension index varies from 2 for very simple object to 3 for objects with complex pattern.

Implementation

Let's define a class, FractalDimObject, that encapsulates the calculation of fractal dimension and the creation of wrapping boxes. 

This class provides two constructors:
  1. The default constructor, __init__, which accepts a 3D array xyz and a threshold value near 1.
  2. An alternative constructor, build, which generates a 3D array of shape (size, size, size) to simulate a 3-dimension object.
class FractalDimObject(object):
    def __init__(self, xyz: np.array, threshold: float) -> None:
        self.xyz = xyz
        self.threshold = threshold


    @classmethod
    def build(cls, size: int, threshold: float) -> Self:
        _xyz = np.zeros((size, size, size))

        # Create a 3D fractal-like structure such as cube
        for x in range(size):            # Width
            for y in range(size):        # Depth
                for z in range(size):    # Height
                    if (x // 2 + y // 2) % 2 == 0:    # Condition for non-zero values
                        _xyz[x, y, z] = random.gauss(size//2, size)

        return cls(_xyz, threshold)

The alternative constructor generates a test array for evaluation purposes. It starts by initializing the array with values of 0.0, and then assigns random Gaussian-distributed values to a specific subset of the array

The values used in define the 3D object is visualized below.

Fig. 1 Visualization of the point cloud representing 3D object


The __call__ method performs the fractal dimension calculation and tracks the relationship between box counts and sizes in three stages:
  1. It determines the box sizes for non-zero elements in the array.
  2. It counts the number of boxes for each size.
  3. It applies a linear regression to the logarithms of the box sizes and their corresponding counts.
def __call__(self) -> (np.array, List[int], List[int]):
     # Step 1 Extract the sizes of array
     sizes = self.__extract_sizes()
     sizes_list = list(sizes)
     sizes_list.reverse()

     # Step 2 Count the number of boxes of each size
     counts = [self.__count_boxes(int(size)) for size in sizes_list]

     # Step 3 Fit the points to a line log(counts) = a.log(sizes) + b
     coefficients = np.polyfit(np.log(sizes), np.log(counts), 1)
     return -coefficients[0], sizes, counts

The __extract_sizes method is responsible for generating the box sizes, as detailed in the Appendix. The implementation of __count_boxes, which counts the wrapping boxes for a given size, follows a similar approach to the method used in calculating the fractal dimension of images.

def __count_boxes(self, box_size: int) -> int:
     sx, sy, sz = self.xyz.shape
     count = 0
        
      for i in range(0, sz, box_size):
          for j in range(0, sy, box_size):
             for k in range(0, sz, box_size):
                  # Wraps the non-zero values (object) with boxes
                 data = self.xyz[i:i+box_size, j:j+box_size, k:k+box_size]

                 if np.any(data):    # For non-zero values (inside object)
                    count += 1
     return count


Evaluation

Let's compute the fractal dimension of the array representing a 3D object with an initial 3D sampling grid 1024 x 1024 x 1024 

import math

grid_size = 1024      # Grid size 
threshold = 0.92
        
fractal_dim_object = FractalDimObject.build(grid_size, threshold)
coefficient, counts, sizes = fractal_dim_object()
print(coefficient)
      

Output: 2.7456

Finally let's plot the profile of box sizes vs. box counts.


Fig. 2 Plot of box sizes vs box counts size = exp(dim*counts)

The plot reflects the linear regression of log (size) and log (counts).

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

Generate the array of box sizes used to compute the number of boxes to cover a given 3D object.

def __extract_sizes(self) -> np.array:
     # Remove values close to 1.0
     filtered = (self.xyz < self.threshold)
        
     # Minimal dimension of box size
     min_dim = min(filtered.shape)

     # Greatest power of 2 less than or equal to p
     n = 2 ** np.floor(np.log(min_dim) / np.log(2))
        
     # Extract the sizes
     size_x: int = int(np.log(n) / np.log(2))

     return np.arange(size_x, 1, -1) * 2