Showing posts with label scikit-learn. Show all posts
Showing posts with label scikit-learn. Show all posts

Wednesday, June 12, 2024

Principal Geodesic Analysis

Target audience: Advanced
Estimated reading time: 7'

Principal Component Analysis (PCA) is essential for dimensionality and noise reduction, feature extraction, and anomaly detection. However, its effectiveness is limited by the strong assumption of linearity. 
Principal Geodesic Analysis (PGA) addresses this limitation by extending PCA to handle non-linear data that lies on a lower-dimensional manifold.


Table of contents
       Tangent PCA
       Setup
       Euclidean PCA
Follow me on LinkedIn

What you will learn: How to implement principal geodesic analysis as extension of principal components analysis on a manifolds tangent space.

Notes

  • Environments: Python  3.10.10, Geomstats 2.7.0, Scikit-learn 1.4.2
  • This article assumes that the reader is somewhat familiar with differential and tensor calculus [ref 1]. Please refer to our previous articles related to geometric learning [ref 234].
  • 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

This article is the ninth installments of our ongoing series focused on geometric learning. As with previous articles, we utilize the Geomstats Python library [ref5] to implement concepts associated with geometric learning.

NoteSummaries of my earlier articles on this topic can be found in the Appendix

As a reminder, 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. 

This article revisits the widely used unsupervised learning technique, Principal Component Analysis (PCA), and its counterpart in non-Euclidean space, Principal Geodesic Analysis (PGA).
The content of this article is as follows:
  1. Brief recap of PCA
  2. Overview of key components of differential geometry
  3. Introduction to PCA on tangent space using the logarithmic map
  4. Implementation in Python using the Geomstats library

   Principal components

    Principal component analysis

Principal Component Analysis (PCA) is a technique for reducing the dimensionality of a large dataset, simplifying it into a smaller set of components while preserving important patterns and trends. The goal is to reduce the number of variables of a data set, while preserving as much information as possible.

For a n-dimensional data, PCA tries to put maximum possible information in the first component c0, then maximum remaining information in the second c1 and so on, until having something like shown in the scree plot below.
Finally, we select the p << n top first components such as:\[\left \{ c_{i} :\ \ \sum_{i=1}^{p} c_{i} < T \ and \ \ c_{1} \geqslant c_{2} .... \geqslant c_{p} \right \}\]
The principal components are actually the eigenvectors of the Covariance matrix.
The first thing to understand about eigenvectors and eigenvalues is that they always appear in pairs, with each eigenvector corresponding to an eigenvalue. Additionally, the number of these pairs matches the number of dimensions in the data.

For instance, in a 3-dimension space, the eigenvalues are extracted from the following 3x3 symmetric Covariance matrix:\[\begin{bmatrix} cov(x, x)) & cov(x,y) & cov(x, z)\\ cov(x, y) & cov(y,y) & cov(y,z) \\ cov(x, z) & cov(y, z) & cov(z, z) \end{bmatrix}\] as cov(a, b) = cov(b, a).
The eigenvectors of the Covariance matrix (principal components) are the directions of the axes where there is the most variance (most information).
Assuming n normalized data points xi, the first principal component (with most significant eigenvalue) is defined as: \[P_{1}= arg\max_{||p||=1} \sum_{i=1}^{n} \left ( p.x_{i} \right )^{2} \ \ \ (1)\]where '.' is the dot product.


Differential geometry

Extending principal components to differentiable manifolds requires basic knowledge of differential and Riemann geometry introduced in previous articles [ref  2, 3, 4].

Smooth manifold
A smooth (or differentiable) manifold is a topological manifold equipped with a globally defined differential structure. Locally, any topological manifold can be endowed with a differential structure by applying the homeomorphisms from its atlas and the standard differential structure of a vector space.

Tangent space
At every point P on a differentiable manifold, one can associate a tangent space, which is a real vector space that intuitively encompasses all the possible directions in which one can move tangentially through P. The elements within this tangent space at P are referred to as the tangent vectors, tgt_vector at P
This concept generalizes the idea of a vector originating from a specific point in Euclidean space. For a connected manifold, the dimension of the tangent space at any point is identical to the dimension of the manifold itself.

Geodesics
Geodesics are curves on a surface that only turn as necessary to remain on the surface, without deviating sideways. They generalize the concept of a "straight line" from a plane to a surface, representing the shortest path between two points on that surface.
Mathematically, a curve c(t) on a surface S is a geodesic if at each point c(t), the acceleration is zero or parallel to the normal vector:\[\frac{d^{2}}{dt^{2}}c(t) = 0 \ \ or \ \ \frac{d^{2}}{dt^{2}}c(t).\vec{n}=\frac{d^{2}}{dt^{2}}c(t)\]
Fig. 1 Illustration of a tangent vector and geodesic on a sphere

Logarithmic map
Given two points P and Q on a manifold, the vector on the tangent space v from P to Q is defined as: \[\left \| log_{P}(Q) \right \| = \left \| v \right \|\]

Tangent PCA

On a manifold, tangent spaces (or plane) are local euclidean space for which PCA can be computed. The purpose of Principal Geodesic Analysis is to project the principal components on the geodesic using the logarithmic map (inverse exponential map).
Given a mean m of n data points x[i] on a manifold with a set of geodesics at m, geod, the first principal component on geodesics is defined as:\[\begin{matrix} P_{1}=arg\max_{\left \| v \right \| = 1}\sum_{i=1}^{n}\left \langle v.log_{m}\pi _{geod}(x_{i}) \right \rangle^{2} \\ \pi_{geod}(x_{i}) = arg\max_{g \in geod} \left \| log_{m}(x_{i})-log_{m}(g) \right \|^2 \end{matrix}\]
Fig. 2 Illustration of a tangent vector and geodesic on a sphere

The mean point on the data point x[I] on the manifold can be defined as either a barycenter or a Frechet Mean [ref 6]. Our implementation in the next section relies on the Frechet mean.

  Implementation

For the sake of simplicity, we illustrate the concept of applying PCA on manifold geodesics using a simple manifold, Hypersphere we introduced in a previous article, Differentiable Manifolds for Geometric Learning: Hypersphere

    Setup

Let's encapsulate the evaluation of principal geodesics analysis in a class HyperspherePCA. We leverage the class HyperphereSpace [ref 7] its implementation of random generation of data points on the sphere.

ff from manifolds.hyperspherespace import HypersphereSpace
I   import numpy as np

  
    class HyperspherePCA(object):
  def __init__(self):
     self.hypersphere_space = HypersphereSpace(equip=True)

  def sample(self, num_samples: int) -> np.array:
     return self.hypersphere_space.sample(num_samples)


    Euclidean PCA

First let's implement the traditional PCA algorithm for 3 dimension (features) data set using scikit-learn library with two methods:
  • euclidean_pca_components to extract the 3 eigenvectors along the 3 eigenvalues
  • euclidean_pca_transform to project the evaluation data onto 3 directions defined by the eigenvectors.
    from sklearn.decomposition import PCA

   
    @staticmethod
    def euclidean_pca_components(data: np.array) -> (np.array, np.array):
   num_components = 3
   pca = PCA(num_components)
   pca.fit(data)
   return (pca.singular_values_, pca.components_)

    @staticmethod
    def euclidean_pca_transform(data: np.array) -> np.array:
   num_components = 3
   pca = PCA(num_components)
 
        return pca.fit_transform(data

We compute the principal components on 256 random 3D data points on the hypersphere.

nu num_samples = 256
pca_hypersphere = HyperspherePCA()

# Random data on Hypersphere
data = pca_hypersphere.sample(num_samples)
eigenvalues, components = HyperspherePCA.euclidean_pca_components(data)
transformed = pca_hypersphere.euclidean_pca_transform(data)        

print(f'\nPrincipal components:\n{components}\nEigen values: {eigenvalues')
print(f'\nTransformed data:\n{transformed}')

Output
Principal components:
[[ 0.63124842 -0.7752775  -0.02168467]
 [ 0.68247094  0.54196625  0.49041411]
 [ 0.36845467  0.32437229 -0.8712197 ]]
Eigen values: [9.84728751 9.0183337  8.65835924]

Important note: The 3 eigenvalues are similar because the input data is random.

    Tangent PCA on hypersphere 

The private method __tangent_pca computes principal components on the tangent plane using the logarithmic map. As detailed in the previous section, this implementation employs the Frechet mean as the base point on the manifold, which is the argument of the Geomstats method fit.

The method tangent_pca_components extracts the principal components computed using the logarithmic map. Finally the method tangent_pca_transform projects input data along the 3 principal components, similarly to its Euclidean counterpart, euclidean_pca_transform .

   from geomstats.learning.pca import TangentPCA
from geomstats.learning.frechet_mean import FrechetMean


def tangent_pca_components(self, data: np.array) -> np.array:
    tgt_pca = self.__tangent_pca(data)
    return tgt_pca.components_

def tangent_pca_transform(self, data: np.array) -> np.array:
    tgt_pca = self.__tangent_pca(data)
    tgt_pca.transform(data)

    
def __tangent_pca(self, data: np.array) -> TangentPCA:
     sphere = self.hypersphere_space.space
     tgt_pca = TangentPCA(sphere)
       
          # We use the Frechet mean as center of data points on the hypersphere
     mean = FrechetMean(sphere, method="default")
     mean.fit(data)
          # Frechet mean estimate
     estimate = mean.estimate_
        
          # Invoke Geomstats fitting method for PCA on tangent plane
          tgt_pca.fit(data, base_point=estimate)
     return tgt_pca

    We evaluate the principal geodesic components on hypersphere using the same 256 points randomly generated on the hypersphere and compared with the components on the Euclidean space.

nu  num_samples = 256
pca_hypersphere = HyperspherePCA()
        
# Random data point generated on Hypersphere
data = pca_hypersphere.sample(num_samples)
        
_, components = HyperspherePCA.euclidean_pca_components(data)
tangent_components = pca_hypersphere.tangent_pca_components(data)
print(f'\nEuclidean PCA components:\n{components}\nTangent Space PCA components:\n{tangent_components}')

Output:
Euclidean PCA components:
[[ 0.63124842 -0.7752775  -0.02168467]
 [ 0.68247094  0.54196625  0.49041411]
 [ 0.36845467  0.32437229 -0.8712197 ]]
Tangent Space PCA components:
[[ 0.8202982  -0.45814582  0.34236423]
 [ 0.45806882  0.16783651 -0.87292832]
 [-0.34246724 -0.87288792 -0.34753831]]

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:

Monday, May 13, 2024

Logistic Regression on Riemann Manifolds

Target audience: Advanced
Estimated reading time: 6'

Traditional linear models in machine learning, such as logistic regression, struggle to grasp the complex characteristics of data in very high dimensions. 
Symmetric Positive Definite manifolds improve the output quality of logistic regression by enhancing feature representation in a lower-dimensional space.


Table of content
        Logistic regression
        SPD manifolds
        Logarithmic map
        Setup
        Data generation
        Manifold creation
        Euclidean space
Follow me on LinkedIn

What you will learn: How to implement and validate a binary logistic regression classifier on SPD manifolds using affine invariant and log Euclidean metrics.

Notes

  • Environments: Python  3.10.10, Geomstats 2.7.0, Scikit-learn 1.4.2
  • This article assumes that the reader is somewhat familiar with differential and tensor calculus [ref 1]. Please refer to our previous articles related to geometric learning [ref 234].
  • 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

This article is the eight installments of our ongoing series focused on geometric learning. In this installment, we utilize the Geomstats Python library [ref5].
NoteSummaries of my earlier articles on this topic can be found in the Appendix

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. 

Using logistic regression for classification on low-dimensional data manifolds offers several benefits:
  • Simplicity and interpretability: The model provides clear insights into the relationship between the input features and the probability of belonging to a certain class.
  • Efficiency: On low-dimensional manifolds, logistic regression is computationally efficient.
  • Good performance in linearly separable cases: The logistic regression performs exceptionally well if the data in the low-dimensional manifold is linearly separable.
  • Robustness to overfitting: In lower-dimensional spaces, the risk of a simpler model such as the logistic regression to overfit is generally reduced.
  • Support for non-linear boundaries: Although linear, the logistic regression can handle non-linear boundaries in low-dimensional space than Euclidean space.

This article relies on the Symmetric Positive Definite (SPD) Lie group of matrices as our manifold for evaluation. We will introduce, review or describe:
  1. Logistic regression as a binary classifier
  2. SPD  matrices
  3. Logarithms and exponential maps on manifolds introduced in (Geometric Learning in Python: Manifolds)
  4. Riemannian metrics associated to SPD
  5. Implementation of binary logistic regression using Scikit-learn and Geomstats Python libraries
  6. Verification using randomly generated SPDs and cross-validation.

Logistic regression on manifolds

Logistic regression

Let's review the ubiquitous binary logistic regression. For a set of two classes C = {0, 1} the probability of predicting the correct class given a features set x and a model weights w is defined by sigmoid, sigm transform:\[p(C|\mathbf{x},\mathbf{w})=p^k(1-p)^{1-k}\ \ \ \ p =sigm(w_{0}+\mathbf{w}^{T}\mathbf{x})= \frac{1}{1-e^{-w_{0}-\mathbf{w}^{T}\mathbf{x}}}\]The binary classifier is then defined as C := 1 <=> p(C=1|x, w) >= 0.5 and C := 0 <=> p(C=1|x, w) < 0.5.
For an introduction to basic logistic regression and its implementation for beginners, check out this detailed guide: Logistic Regression Explained and Implemented in Python (ref 6).

SPD manifolds

Let's introduce our manifold defined as the group of symmetric positive definite (SPD) matrices. 
SPD matrices are used in a wide range of applications:
  • Diffusion tensor imaging (analysis of diffusion of molecules and proteins)
  • Brain connectivity
  • Dimension reduction through kernels
  • Robotics and dynamic systems
  • Multivariate principal component analysis
  • Spectral analysis and signal reconstruction
  • Partial differential equations numerical methods
  • Financial risk management.

A square matrix A is symmetric if it is identical to its transpose, meaning that if aaij are the entries of A, then aaij=ajiaij. This implies that A can be fully described by its upper triangular elements.
aij

A square matrix A is positive definite if, for every non-zero vector b, the product bTAb >= 0.

If a matrix A is both symmetric and positive definite, it is referred to as a symmetric positive definite (SPD) matrix. This type of matrix is extremely useful and appears in various real-world applications. A prominent example in statistics is the covariance matrix, where each entry represents the covariance between two variables (with diagonal entries indicating the variances of individual variables). Covariance matrices are always positive semi-definite (meaning bTAb0), and they are positive definite if the covariance matrix has full rank, which occurs when each row is linearly independent from the others.
The collection of all SPD matrices of size n×n forms a manifold.

Logarithmic map 

As discussed in Geometric Learning in Python: Manifolds, the exponential & logarithmic maps and parallel transportation are crucial for Riemannian approaches in machine learning. On any manifold, distance, dist(x, y) are defined as geodesics that correspond to straight lines in Euclidean space.

Let consider a vector x to y on a tangent space at point y. Operations on the point on the manifold  rely on the exponential map (projection) onto the manifold. 
The table below show the equivalent operations between Euclidean space and manifolds.

                          Operations                    Euclidean                        Manifold
In the case of the binary logistic regression, the prediction on the manifold is defined by the exponential map expx: \[p(y|\mathbf{x}, \mathbf{w})=exp_{x}(y, sigm(w_{0}+\mathbf{w}^T.\mathbf{x}))\]Let select two Riemannian metrics for the SPD manifolds [ref 7].


Affine invariant Riemannian metric

For any two symmetric positive definite (SPD) matrices 𝐴 and 𝐵, the Affine Invariant Riemannian Metric (AIRM) between them is defined as: \[d(A, B)=\left \| log\left ( A^{-1/2}BA^{-1/2} \right ) \right \|_F\]

Log-euclidean Riemannian metric

Given a symmetric positive definite matrix SPD at point S and a tangent space TsSPD, the logarithmic and exponential maps can be expressed as: \[\begin{matrix} log_{s}(f(s))=D_{log(s)}\ exp.\left (log(f(s) -log(s) \right )\\ exp_{s}\left ( T_{f(s)} \right )=exp\left ( log(s)+D_{s}log.T_{f(s))} \right ) \end{matrix}\]
Fig. 1 Illustration of the log-euclidean metric for SPD


Implementation

Setup

First, let's create a data class, SPDTestData that encapsulates the training features X and label y. This class will be used to validate our implementation of logistic regression on SPD manifolds using various metrics, as well as in Euclidean space.

@dataclass
class SPDTestData:
  X: np.array      # Features
  y: np.array       # Labels

  def flatten(self) -> NoReturn:
      shape = self.X.shape
      self.X = self.X.reshape(shape[0], shape[1]*shape[2])

The flatten method vectorizes each 2-dimension SPD matrix entry in the training set to be processed by the scikit-learn cross validation function.

We wrap the generation of random data, SPD manifolds and the evaluation of various metrics in the class BinaryLRManifold.

class BinaryLRManifold(object):
   def __init__(self, n_features: int, n_samples: int):
      self.n_features = n_features
      self.n_samples = n_samples


   def generate_random_data(self) -> SPDTestData:
      y = np.stack([np.random.randint(0, 2) for _ in range(self.n_samples)])
      X = np.stack([self.__generate_spd_data() for _ in range(self.n_samples)])
  
    return SPDTestData(X, y)

The generation of the labeled training set uses the numpy random values generation method.

Data generation

The method __generate_spd_data create symmetric positive definite n_features x n_features matrices by computing eigenvalues for the diagonal component, reducing the upper triangle values and replicating them to the lower triangle.

def __generate_spd_data(self) -> np.array:
   epsilon = 1e-6
        
   mat = np.random.rand(self.n_features, self.n_features)
   mat = (mat + mat.T)/2
        
   eigenvalues = np.linalg.eigvals(mat)
   min_eigen = np.min(eigenvalues)
   if min_eigen <= 0:
      mat += (np.eye(self.n_features)*(-min_eigen + epsilon))

   return mat

Manifold generation

Creating an SPD matrix is straightforward:
  1. Instantiate the Geomstats SPDMatrices class
  2. Equip it with a Riemannian metric.
def create_spd(self, riemannian_metric: RiemannianMetric) -> SPDMatrices:
   spd = SPDMatrices(self.n_features, equip=False)
   spd.equip_with_metric(riemannian_metric)

   return spd

The following code snippet creates two SPD manifolds with affine metric and log Euclidean Riemannian metrics.

from geomstats.geometry.spd_matrices import SPDAffineMetric, SPDLogEuclideanMetric

n_samples = 10000
n_features = 16

binary_lr_on_spd = BinaryLRManifold(n_features, n_samples)

# Create a SPD matrix and assigned a Affine metric
spd = binary_lr_on_spd.create_spd(SPDAffineMetric)

# Create a SPD matrix and assigned a Log Euclidean metric
spd = binary_lr_on_spd.create_spd(SPDLogEuclideanMetric)


Validation

The initial phase involves verifying our implementation of the metrics related to SPD manifolds. This is achieved by calculating the cross-validation score for SPD matrices containing random values between [0, 1] and ensuring that the average score approximates 0.5.

Euclidean space

We utilize the logistic regression class and the cross_validate method from Scikit-learn once the contents of the matrix have been converted into a vector form.

from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_validate

@staticmethod
def evaluate_euclidean(spd_data: SPDTestData) -> Dict[AnyStr, np.array]:
   model = LogisticRegression()
   
   # Reduce the matrix into a vector for sklearn cross-validation
   spd_data.flatten()

   return cross_validate(model, spd_data.X, spd_data.y)


The test code used a training set of 6000 samples of 16 x 16 (256) SPD matrices. The binary logistic regression in the Euclidean space as a mean cross validation score of 0.487 instead of 0.5.

n_samples = 6000
n_features = 16

binary_lr_on_spd = BinaryLRManifold(n_features, n_samples)
train_data = binary_lr_on_spd.generate_random_data()
print(f'Training data shape: {train_data.shape()}')

result_dict = binary_lr_on_spd.evaluate_euclidean(train_data)
mean_test_score = np.mean(result_dict["test_score"])
print(f'Cross validation: {result_dict["test_score"]} with mean: {mean_test_score}')

Output
Cross validation: [0.478 0.513 0.497 0.474 0.471] with mean: 0.487

Classification on SPD manifold

To utilize scikit-learn's cross-validation features, the SPD matrix must first be differentiated on its tangent space before applying logistic regression. These two steps are executed using a scikit-learn pipeline.

@staticmethod
def evaluate_spd(spd_data: SPDTestData, spd_matrices: SPDMatrices) -> Dict[AnyStr, np.array]:
   from geomstats.learning.preprocessing import ToTangentSpace
   from sklearn.pipeline import Pipeline

   pipeline = Pipeline(
       steps=[ ('features', ToTangentSpace(space = spd_matrices)),
                    ('classifier', LogisticRegression())
       ]
   )
   return cross_validate(pipeline, spd_data.X, spd_data,y)

We employ the same training setup as used in the evaluation of logistic regression in Euclidean space, but we apply the log Euclidean (SPDLogEuclideanMetric) and affine invariant (SPDAffineInvariant) metrics. The mean values of the cross-validation scores are 0.492 and 0.5, respectively, which significantly improve upon the results from the Euclidean scenario.

n_samples = 6000
n_features = 16

binary_lr_on_spd = BinaryLRManifold(n_features, n_samples)
train_data = binary_lr_on_spd.generate_random_data()

spd = binary_lr_on_spd.create_spd(SPDLogEuclideanMetric)
result_dict = binary_lr_on_spd.evaluate_spd(train_data, spd)
mean_test_score = np.mean(result_dict["test_score"])

print(f'Cross validation: {result_dict["test_score"]} with mean: {mean_test_score}')


Output for Log Euclidean metric
Cross validation: [0.495  0.504 0.498 0.491 0.470] with mean: 0.492

Output for affine invariant metric
Cross validation: [0.514  0.490 0.490 0.490 0.504] with mean: 0.500

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:

#geometriclearning #riemanngeometry #manifold #ai #python #geomstats #Liegroups #kmeans