Showing posts with label Hypersphere. Show all posts
Showing posts with label Hypersphere. Show all posts

Saturday, June 29, 2024

Fisher Information Matrix

Target audience: Advanced
Estimated reading time: 7'
The Fisher Information Matrix plays a crucial role in various aspects of machine learning and statistics. Its primary significance lies in providing a measure of the amount of information that an observable random variable carries about an unknown parameter upon which the probability depends.


Table of contents
       Key elements
       Use cases
Follow me on LinkedIn

What you will learn: How to estimate and visualize the Fisher information matrix for Normal and Beta distributions on a hypersphere.

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]. Please refer to our previous articles related to geometric learning listed on Appendix.
  • Source code is available at  Github.com/patnicolas/Data_Exploration/Information Geometry
  • 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 10th installments of our ongoing series focused on geometric learning. It introduces some basic elements of information geometry as an extension of differential geometry. As with previous articles, we utilize the Geomstats Python library [ref. 2] 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. 

Here is a synapsis of this article
  1. Brief introduction to information geometry
  2. Overview and mathematical formulation of the Fisher information matrix
  3. Computation of the Fisher metric to Normal and Beta distributions
  4. Implementation in Python using the Geomstats library

Information geometry

Information geometry applies the principles and methods of differential geometry to problems in probability theory and statistics [ref 3]. It studies the manifold of probability distributions and provides a natural framework for understanding and analyzing statistical models.

Key elements

  • Statistical manifolds: Families of probability distributions are considered as a manifold, with each distribution representing a point on this manifold.
  • Riemannian metrics: The Fisher information metric is commonly used to define a Riemannian metric on the statistical manifold. This metric measures the amount of information that an observable random variable carries about an unknown parameter.
  • Divergence measures: Divergence measures like the Kullback-Leibler (KL) divergence, which quantify the difference between two probability distributions.
  • Connections and curvature: Differential geometry concepts such as affine connections and curvature are used to describe the geometric properties of statistical models (i.e. α-connection family).
  • Dualistic inference: Exponential and mixture connections provide a rich structure for statistical inference.

Use cases

Here is a non-exclusive list of application of information geometry
  • Statistical Inference: Parameter estimation, hypothesis testing, and model selection (i.e. Bayesian posterior distributions and in the development of efficient sampling algorithms like Hamiltonian Monte Carlo)
  • Optimization: Natural gradient descent method uses the Fisher information matrix to adjust the learning rate dynamically, leading to faster convergence compared to traditional gradient descent.
  • Finance: Modeling uncertainties and analyzing statistical properties of financial models.
  • Machine Learning: Optimization of learning algorithms (i.e. Understanding the EM algorithm used in statistical estimation for latent variable model)
  • Neuroscience: Neural coding and information processing in the brain by modeling neural responses as probability distributions.
  • Robotics: Development of probabilistic robotics, where uncertainty and sensor noise are modeled using probability distributions.
  • Information Theory: Concepts for encoding, compression, and transmission of information.

Fisher information matrix

The Fisher information matrix is a type of Riemannian metric that can be applied to a smooth statistical manifold [ref 4]. It serves to quantify the informational difference between measurements. The points on this manifold represent probability measures defined within a Euclidean probability space, such as the Normal distribution. Mathematically, it is represented by the Hessian of the Kullback-Leibler divergence.

Let's consider a statistical manifold with coordinates (or parameters) θ and its probability density functions over an interval X as follow:\[P= \left \{ p(x, \theta); \ x \in X \ \int_{R}^{} p(x, \theta) dx = 1\right \}\]The Fisher metric is a Riemann metric tensor defined as the expectation of the partial derivative of the negative log likelihood over two coordinates θ.\[g_{ij}(\theta) = -E\left [ \frac{\partial^2\ log\ p(x,\theta) }{\partial \theta_{i}\partial\theta_{j}} \right ] = - \int_{R}^{}{\frac{\partial^2\ log\ p(x,\theta)) }{\partial \theta_{i}\partial\theta_{j}}}p(x, \theta)dx\]

The Fisher information or Fisher-Rao metric quantifies the amount of information in the data regarding a parameter θ. The Fisher-Rao metric, an intrinsic measure, enables the analysis of a finite, n-dimensional statistical manifold M.\[ds=\sum_{i=1}^{p}{\sum_{j=1}^{p}}g_{ij}\theta^{i}\theta^{j}\]
The Fisher metric for the normal distribution θ = {μ, σ} is computed as:\[\mathfrak{I}(\mu, \sigma)=-\textit{E}_{x-p}\begin{bmatrix} \frac{\partial ^2\ log\ p(\theta)}{\partial \mu^2} & \frac{\partial ^2\ log\ p(\theta)}{\partial \mu \partial \sigma} \\ \frac{\partial ^2\ log\ p(\theta)}{\partial \sigma \partial \mu} & \frac{\partial ^2\ log\ p(\theta)}{\partial \sigma^2} \end{bmatrix} = \begin{bmatrix} \sigma^{-2} & 0\\ 0 & 2\sigma^{-2} \end{bmatrix}\]
The Fisher metric for the beta distribution  θ = {α, β} is computed as:\[\varphi (z)=\frac{d^2}{dz^2}\ log \ \Gamma (z)\]
\[\mathfrak{I(\alpha,\beta)}=-\textit{E}_{x-p}\begin{bmatrix} \frac{\partial ^2\ log\ p(\theta)}{\partial \alpha^2} & \frac{\partial ^2\ log\ p(\theta)}{\partial \alpha \partial \beta} \\ \frac{\partial ^2\ log\ p(\theta)}{\partial \beta \partial \alpha} & \frac{\partial ^2\ log\ p(\theta)}{\partial \beta^2} \end{bmatrix}\]
\[\mathfrak{I(\alpha,\beta)}=\begin{bmatrix} \varphi(\alpha)-\varphi(\alpha+\beta) & -\varphi(\alpha+\beta)\\ -\varphi(\alpha+\beta) & \varphi(\beta)-\varphi(\alpha+\beta) \end{bmatrix}\]

Implementation

We leverage the following classes defined in the previous articles:
Let's first define a base class for all distributions to be defined on a hypersphere [ref 5].

class GeometricDistribution(object):
    _ZERO_TGT_VECTOR = [0.0, 0.0, 0.0]

    def __init__(self) -> None:
        self.manifold = HypersphereSpace(True)


    def show_points(self, num_pts: int, tgt_vector: List[float] = _ZERO_TGT_VECTOR) -> NoReturn:
        # Random point generated on the hypersphere
        manifold_pts = self._random_manifold_points(num_pts, tgt_vector)
        
        # Exponential map used to project the tgt vector on the hypersphere
        exp_map = self.manifold.tangent_vectors(manifold_pts)

        for v, end_pt in exp_map:
            print(f'Tangent vector: {v} End point: {end_pt}')
        self.manifold.show_manifold(manifold_pts)

The purpose of the method show_points is to display the various data point with optional tangent vector on the hypersphere. The argument num_pts specifies the number of random points to be defined in the hypersphere. The tangent vector is displayed if the argument tgt_vector not defined as the origin (_ZERO_TGT_VECTOR).


Normal distribution

The class NormalHypersphere encapsulates the display of the normal distribution on the hypersphere. The constructor initialized the normal distribution implemented in the Geomstats library.
The method show_distribution display num_pdfs probability density function over a set of num_manifold_pts, manifold points on the hypersphere. This specific implementation uses only two points. The Fisher-Rao metric is computed using the metric.geodesic Geomstats method.
The metric is applied to 100 points along the geodesic between the two points A and B. Finally, the density functions, pdfs are computed by converting the metric values to the NormalDistribution.point_to_pdf Geomstats method.

class NormalHypersphere(GeometricDistribution):

    def __init__(self) -> None:
        from geomstats.information_geometry.normal import NormalDistributions

        super(NormalHypersphere, self).__init__()
        self.normal = NormalDistributions(sample_dim=1)


    def show_distribution(self, num_pdfs: int, num_manifold_pts: int) -> NoReturn:
        manifold_pts = self._random_manifold_points(num_manifold_pts)
        A = manifold_pts[0]
        B = manifold_pts[1]

        # Apply the Fisher metric for the two manifold points on a Hypersphere
        geodesic_ab_fisher = self.normal.metric.geodesic(A.location, B.location)
        t = gs.linspace(0, 1, 100)

        # Generate the various density functions associated to the Fisher metric between the
        # two points on the hypersphere
        pdfs = self.normal.point_to_pdf(geodesic_ab_fisher(t))
        x = gs.linspace(0.2, 0.7, num_pdfs)

        for i in range(num_pdfs):
            plt.plot(x, pdfs(x)[i, :]/20.0)   # Normalization factor
        plt.title(f'Normal distribution on Hypersphere')
        plt.show()

Let's plot 2 randomly sampled data points associated with a tangent_vector on Hypersphere (1) then visualize 40 normalized normal probability density distributions (2).

normal_dist = NormalHypersphere()
num_points = 2
tangent_vector = [0.4, 0.7, 0.2]

         # 1. Display the 2 data points on the hypersphere
num_manifold_pts = normal_dist.show_points(num_points, tangent_vector)

        # 2. Visualize the 40 normal probabilities density functions
num_pdfs = 40
succeeded = normal_dist.show_distribution(num_pdfs, num_points)

Fig. 1 Two random data points on a Hypersphere with their tangent vectors 



Fig. 2 Visualization of Normal distribution between two random points on a hypersphere


Beta distribution

Let's wrap the evaluation of the Beta distribution on a hypersphere into the class BetaHypersphere that inherits GeometriDistribution. It leverages the BetaDistributions class in Geomstats. 

class BetaHypersphere(GeometricDistribution):

    def __init__(self) -> None:
        from geomstats.information_geometry.beta import BetaDistributions

        super(BetaHypersphere, self).__init__()
        self.beta = BetaDistributions()


    def show_distribution(self, num_manifold_pts: int, num_interpolations: int) -> NoReturn:

        # 1. Generate random points on Hypersphere using Von Mises algorithm
        manifold_pts = self._random_manifold_points(num_manifold_pts)
        t = gs.linspace(0, 1.1, num_interpolations)[1:]
        # 2. Define the beta pdfs associated with each
        beta_values_pdfs = [self.beta.point_to_pdf(manifold_pt.location)(t) for manifold_pt in manifold_pts]

        # 3. Generate, normalize and display each Beta distribution
        for beta_values in beta_values_pdfs:
            min_beta = min(beta_values)
            delta_beta = max(beta_values) - min_beta
            y = [(beta_value - min_beta)/delta_beta  for beta_value in beta_values]
            plt.plot(t, y)
        plt.title(f'Beta distribution on Hypersphere')
        plt.show()


The method show_distribution generates random points on the Hypersphere (1)  and compute the beta density function at these points using the Geomstats BetaDistributions.point_to_pdf (2).
The values generated by the pdfs are normalized then plotted (3)


Let's plot 10 randomly sampled data points on Hypersphere (1) then visualize 200 normalized beta probability density distributions (2).

beta_dist = BetaHypersphere()
        
num_interpolations = 200
num_manifold_pts = 10
    # 1. Display the 10 data points on the hypersphere
beta_dist.show_points(num_manifold_pts)

   # 2. Visualize the probabilities density functions with interpolation points
succeeded = beta_dist.show_distribution(num_manifold_pts, num_interpolations)
Fig. 3  10 random data points with on a Hypersphere


Fig. 4 Visualization of Beta distributions associated with 10 data points on hypersphere


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:

Wednesday, May 1, 2024

K-means on Riemann Manifolds

Target audience: Advanced
Estimated reading time: 7'

Traditional clustering models often fail complex datasets commonly found in advanced applications like medical imaging, 3D shape analysis, and natural language processing where data is highly interrelated.
K-means on manifolds respects the intrinsic geometry of the data, such as curvature and metric.


Table of contents

       Setup
       Euclidean space
       Hypersphere
Follow me on LinkedIn

What you will learn: How to apply k-means clustering on a Riemann manifold (Hypersphere) using Geomstats, contrasted with its implementation in Euclidean space, using scikit-learn library.

Notes

  • Environments: Python  3.10.10, Geomstats 2.7.0, Scikit-learn 1.4.2, Matplotlib 3.8.3
  • 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 2, 3, 4].
  • 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

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. 5] and explore the ubiquitous K-means clustering algorithm on the hypersphere manifold. The hypersphere which was introduced in a previous piece, Geometric Learning in Python: Manifolds - Hypersphere  and is detailed in the Geomstats API [ref. 6]. 

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

There are many benefits for clustering data on a manifold for complex data sets [ref 7]:
  • Grouping of dense, continuous non-linear data depends on the 'shape' of data
  • Projection to Euclidean space may introduce distortion
  • Loss, distances are better assessed and computed through geodesics than Euclidean metrics (i.e. sphere)

K-means


Among the array of unsupervised learning algorithms, K-means stands out as one of the most well-known. This algorithm has a straightforward goal: to divide the data space so that data points within the same cluster are as similar as possible (intra-cluster similarity), and data points in different clusters are as dissimilar as possible (inter-cluster similarity). K-means aims to identify a predetermined number of clusters in an unlabeled dataset. It employs an iterative approach to finalize the clustering, which depends on the number of clusters specified by the user (denoted by the variable K).

Given K, Ck clusters each with a centroid mk, the input data xi is distributed across the cluster so to minimize the reconstruction error:\[Rerr(K)=\max_{C_k}\sum_{k=1}^{K}\sum_{x_i\in C_k}^{}\left \| x_i-m_k\right \|^2\]

Clustering data on a manifold

To assess and contrast the k-means model in both Euclidean space and on a hypersphere, it's necessary to create clustered data. This involves a two-step process:
  1. Generate a template cluster by employing a random generator on the manifold.
  2. Generate 4 clusters from the template using a special orthogonal Lie group in 3-dimensional space, SO(3).

Randomly generated manifold data

Let's evaluate and compare the following random generators for data points on the hypersphere we introduced in a previous article [ref 8]
  • Uniform distribution
  • Uniform distribution with constraints
  • von Mises-Fisher distribution

Uniform distribution
We start with the basic random uniform generator over interval [0, 1].\[r=rand_{[0,1]}(x)\]The data points for the 4 clusters are visualized in the following plot.
4-Cluster random generation using uniform distribution


Constrained uniform random generator.
In this scenario, we constrain random values r on each of the 3 dimension (or axis) within a sub-interval [ai, bi].\[r=rand_{[0,1]}(x)\ \ \ a_i < r_i < b_i\]

4-Cluster random generation using constrained uniform distribution

von mises-Fisher random generator
This approach relies on a generative mixture-model approach to clustering directional data based on the von Mises-Fisher distribution [ref 9].
Given a d-dimensional unit random vector x on a hypersphere of dimension d-1, the d-variate von Moses-Fisher distribution is defined by the following probability density distribution:\[f(x|\mu , \kappa )=C_d(\kappa).e^{\kappa \mu^Tx} \ \ \  \ C_d(\kappa)=\frac{\kappa^{\frac{d}{2}-1}}{2\pi^{\frac{d}{2}}I_{\frac{d}{2}-1}(\kappa)}\]Id is the Bessel function and Cd is the normalization factor.

I
4-Cluster random generation using Von Mises-Fisher distribution

As anticipated, using a pure uniform random generator distributes data evenly across the hypersphere, rendering it ineffective for evaluating KMeans.
Instead, we will employ the von Mises-Fisher distribution and a constrained uniform random generator to more effectively analyze the performance of KMeans on a Riemann manifold.

Synthetic clusters using SO(3)

We leverage the SO(3) Lie group to replicate the randomly generated cluster.

Although the discussion of Lie groups and special orthogonal group in 3-dimensional space is beyond the scope of this article, here is a short summary:
In differential geometry, Lie groups play a crucial role by connecting the concepts of algebra and geometry. A Lie group is a mathematical structure that is both a group and a differentiable manifold. This means that the group operations of multiplication and taking inverses are smooth (differentiable), and it allows the application of calculus within the group structure.

The Special Orthogonal Lie group in 3-dimension space SO(3) is simply a group of 3 x 3 orthogonal matrices with determinant = 1. These represent rotation in
n-dimensional space and form a compact Lie group.

Implementation

Setup 

Let's wraps the random generators and k-means training methods in a class, KMeansOnManifold.
The von Mises-Fisher generator for the data in the initial cluster is initialized with a mean, _mu and a kappa arbitrary value. The constrained uniform random generator, accepts random values in each dimension x: [-1, -0.35], y: [0.3, 1] and z: [-1. 0.4].
The SO(3) Lie group is initialized without metric (equip=False) to generate the 4 synthetic clusters.

from geomstats.geometry.hypersphere import Hypersphere
from geomstats.geometry.special_orthogonal import SpecialOrthogonal

class KMeansOnManifold(object):

   def __init__(self, num_samples: int, num_clusters: int, random_gen: AnyStr):
     # Step 1: Initialize the manifold
     self.hypersphere = Hypersphere(dim=2, equip=True)

     # Step 2: Generate a single cluster with random data points on hypersphere
     match random_gen:
        case 'random_von_mises_fisher':
             # Select a pivot or mean value
             _mu = self.hypersphere.random_uniform(n_samples=1)
             # Generate the cluster
             cluster = self.hypersphere.random_von_mises_fisher(
                      mu=_mu[0], 
                      kappa=60, 
                      n_samples=num_samples, 
                      max_iter=200)
        
        case 'random_riemann_normal':
             cluster = self.hypersphere.random_riemannian_normal(n_samples=num_samples, max_iter=300)
        
        case 'random_uniform':
            cluster = self.hypersphere.random_uniform(n_samples=num_samples)
        
        case 'constrained_random_uniform'
            # Generate random values with constrains on each dimension.
            y = [x for x in self.hypersphere.random_uniform(n_samples=100000)
            if x[0] <= -0.35 and x[1] >= 0.3 and x[2] <= -0.40]
            cluster = np.array(y)[0:num_samples]

        case _:
            raise ValueError(f'{random_gen} generator is not supported')
        
     # Step 3: Generate other clusters using SO(3) manifolds
     so3_lie_group = SpecialOrthogonal(3, equip=False)
     # Generate the clusters
     self.clusters = [cluster @ so3_lie_group.random_uniform() for _ in range(num_clusters)]



Data in Euclidean space

The data class, KMeansCluster encapsulates the output (centroid and label) of the training of the k-means algorithm on the synthetic clustered data.

@dataclass
class KMeansCluster:
    center: np.array
    label: np.array

We rely k-means implementation in scikit-learn library [ref 10] (class KMeans) on the  to identify the clusters in the Euclidean space, selecting the elkan algorithm, and k-means++ initialization.

def euclidean_clustering(self) -> List[KMeansCluster]:
   from sklearn.cluster import KMeans

   kmeans = KMeans(
       n_clusters=len(self.clusters), 
       init='k-means++', 
       algorithm='elkan',
       max_iter=140)
  
   # Create a data set from points in clusters
   data = np.concatenate(self.clusters, axis=0)
   kmeans.fit(data)

   # Extract centroids and labels
   centers = kmeans.cluster_centers_
   labels = kmeans.labels_

   return [KMeansCluster(center, label) for center, label in zip(centers, labels)]

Output:
Cluster Center: [ 0.56035023 -0.4030522   0.70054776], Label: 0
Cluster Center: [-0.1997325  -0.38496744  0.8826764 ], Label: 2
Cluster Center: [0.04443849 0.86749237 0.46118632], Label: 3
Cluster Center: [-0.83876485 -0.45621187  0.23570083], Label: 1

Clearly, this implementation of k-means was not able to identify the proper clusters

Data on hypersphere

The methods to train k-means on the hypersphere uses the same semantic as its sklearn counterpart. It leverage the Geomstats, RiemannianKMeans class method.

def riemannian_clustering (self) -> List[KMeansCluster]:
    from geomstats.learning.kmeans import RiemannianKMeans

    # Invoke the Geomstats Riemann Means
    kmeans = RiemannianKMeans(space=self.hypersphere, n_clusters=len(self.clusters))
    
    # Build the data set from the clustered data points
    data = gs.concatenate(self.clusters, axis =0)
        
    kmeans.fit(data)

    # Extract predictions, centroids and labels
centers = kmeans.centroids_ labels = kmeans.labels_ return [KMeansCluster(center, label) for center, label in zip(centers, labels)]


Similar to k-means in Euclidean space, we identify the centroids for 4 clusters using 500 randomly generated samples.

num_samples = 500
num_clusters = 4
kmeans = KMeansOnManifold(num_samples, num_clusters,  'random_von_mises_fisher')
kmeans_cluster = kmeans.riemannian_clustering()


Output:
500 random samples on 4 clusters with von-mises-Fisher distribution
Cluster Center: [ 0.17772496 -0.36363422  0.91443097], Label: 2
Cluster Center: [ 0.44403679  0.06735507 -0.89347335], Label: 0
Cluster Center: [ 0.85407911 -0.50905801  0.10681211], Label: 3
Cluster Center: [ 0.90899637  0.02635062 -0.41597025], Label: 1

500 random samples on 4 clusters with constrained uniform distribution
Cluster Center: [-0.05344069 -0.91613807  0.3972847 ], Label: 1
Cluster Center: [ 0.6796575   0.39400079 -0.61873181], Label: 2
Cluster Center: [ 0.51799972 -0.67116261 -0.530299 ], Label: 0
Cluster Center: [ 0.49290501 -0.45790221 -0.73984473], Label: 3

Note: The labels are arbitrary indices assigned to each cluster for the purpose of visualization and validation against true labels.

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, ai, python, geomstats, Liegroups, kmeans