<figure><img src="image_1.png" /></figure>

HAL
open science

# K-Means and Gaussian Mixture Models on Lie Groups: Application to Geometrical Clustering

Samy Labsir, Julien Lesouple, Jean-Yves Tourneret

► To cite this version:

Samy Labsir, Julien Lesouple, Jean-Yves Tourneret. K-Means and Gaussian Mixture Models on Lie Groups: Application to Geometrical Clustering. 2025. hal-05023486

HAL Id: hal-05023486

https://hal.science/hal-05023486v1

Preprint submitted on 7 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# K-Means and Gaussian Mixture Models on Lie Groups: Application to Geometrical Clustering

Samy Labsir, Julien Lesouple, Member, IEEE, and Jean-Yves Tourneret, Fellow, IEEE

***Abstract***—In this article, we derive and implement two new clustering algorithms dedicated to Lie groups, adapted from the well known K-Means algorithm and Gaussian mixture models. More precisely, the K-Means algorithm is reformalized by taking into account the fact that observations belong to a Lie group (LG) with an appropriate intrinsic metric and Gaussian mixture model are defined for data living in LGs. The consistency and the performance of the resulting are numerically validated for data lying on the LG *SE*(2), by comparison with state-of-the-art methods for synthetic data and pseudo-real data generated using an ultra-sound sensor model.

***Index Terms***—Clustering, estimation on Lie group, Gaussian distribution on Lie groups, ultra sound sensor.

## I. INTRODUCTION

Machine learning (ML) methods have become crucial and fundamental during the last decades in different fields of signal and image processing. Indeed, these methods have been intensively used to solve estimation, prediction and classification problems. For example, we can quote supervised methods such as neural networks and support vector machines that have shown interesting results for target radar classification [1] (for instance to classify different types of aircrafts) or pattern detection (for instance to detect anomalies in the trajectory of satellites) [2]. ML methods have also been used to solve unsupervised problems. As an example, it is interesting to note that the K-means algorithm [3] and Gaussian mixture models have good properties to segment radar images or classify targets observed by a teledetection system (RADAR or GNSS) [4] [5]. In some applications, the data to predict or classify is constrained by its geometrical properties, e.g., when performing estimation on the set of orthogonal matrices. For example, in applications such as navigation or computer vision, the sensor provides measurements depending on rotation matrices (associated with the object attitude) that are orthogonal matrices [6]. Another example is the set of symmetric positive definite matrices, which has been used in radar signal processing, since it contains the correlation matrix of the received signal in each detection cell [7], or in image processing to model the tensor structure of an image [8]. In order to be precise and consistent in the applications mentioned above, the dedicated prediction or classification algorithms have to take into account the constraints imposed by the data structure. In the state-of-the-art, ML methods using

Riemannian geometric structure are scarce but exist [9]. Deep neural networks based on the same structure have also been considered in [10].

In the context of this work, we consider Lie groups (LGs), more precisely Riemannian manifolds equipped with a group structure. Parameters lying in LGs are encountered in a plethora of applications in signal processing, robotics, or computer vision [11]. For instance, one can quote registration problems where the transformation between two unregistered images belongs to the *special Euclidean group*, denoted as *SE*(3), or the *similarity group* *SL*(3) [12]. Additionally, in several target tracking applications, it is often necessary to estimate the orientation of a mobile object constrained to the LG of three-dimensional rotation matrices, denoted as *SO*(3) [13]. *Contrary to Riemannian manifolds*, LGs have the advantage to leverage mathematical tools (especially statistics and optimization) making them suitable for statistical learning problems. More precisely, the group operation and smooth properties allows us to define LG integration and derivation operators in a simple way using Haar measure, which leads to the definition of probability density functions (pdfs) on Lie groups or derivation operators that can be used to derive optimization algorithms on LGs. Due to its simplicity of writing and sampling, the Gaussian distribution on Lie groups [14] is basically used for defining statistical models on LGs. ML methods on LGs exist mainly in a supervised framework for tasks such as estimation and filtering [15]–[17]. More precisely, optimization algorithms based on Newton or Gauss-Newton methods, filtering algorithms based on the Kalman or unscented Kalman filters or particle filters have been generalized to matrices belonging to LGs [13], [18]–[20]. On the other hand, clustering algorithms for matrices in LGs are crucial in a variety of applications [21]–[23]. One can mention the classification of rotation matrices lying in the LG *SO*(3), which is useful in radar signal processing when we want to discriminate two close targets characterized by their attitude [24]. Other applications can be found in computer vision and image processing, where several affine transformations linking a reference image to other images lying in the LG *SE*(3) have to be classified [25], [26].

## A. Related works

1) *Clustering on the Riemannian manifold of symmetric positive definite matrices:* some clustering methods have already been investigated for data lying in Riemannian manifolds. For instance, the K-Means and EM algorithms have been proposed in [27] and [28] to cluster covariance matrices. The

<footnote>S. Labsir is with IPSA/TéSA, Toulouse, France.</footnote>

<footnote>J. Lesouple is with Fédération ENAC ISAE-SUPAERO ONERA, Université de Toulouse, France.</footnote>

<footnote>J.-Y. Tourneret is with ENSEEIHT/IRIT/TéSA, Université de Toulouse, France.</footnote>
K-means algorithm investigated in [27] is based on a geodesic distance between symmetric positive definite (SPD) matrices $\{\mathbf{X}_i\}_{i=1}^N$ and a cluster centroid $\mathbf{M}$ (belonging to the manifold of SPD matrices) defined by:

$$d_{\text{geo}}(\mathbf{X}_i, \mathbf{M}) = \left\| \logm(\mathbf{X}_i^{-\frac{1}{2}} \mathbf{M} \mathbf{X}_i^{-\frac{1}{2}}) \right\|_F, \quad (1)$$

where $\logm(\cdot)$ is the matrix logarithm and $\|\cdot\|_F$ is the Frobenius norm. The cluster centroid $\mathbf{M}$ is determined by computing the Karcher mean [29]

$$\widehat{\mathbf{M}} = \underset{\mathbf{M}}{\operatorname{argmin}} \sum_{i=1}^{N} d_{\text{geo}}^{2}(\mathbf{X}_{i}, \mathbf{M}). \quad (2)$$

This optimization problem does not lead to a closed-from expression of $\widehat{\mathbf{M}}$ and has to be solved using a Riemannian gradient algorithm [30].

The EM algorithm studied in [28] enables a parametric clustering of SPD matrices. It is based on a Riemannian Gaussian mixture model, meaning that each observed matrix has the following distribution:

$$p(\mathbf{X}_i | \pi_{1:K}, \mathbf{M}_{1:K}, \sigma_{1:K}^2) = \sum_{k=1}^{K} \pi_k \mathcal{N}_{RG}(\mathbf{X}_i; \mathbf{M}_k, \sigma_k^2), \quad (3)$$

with $\pi_{1:K} = (\pi_1, ..., \pi_K)^T$ such that $\sum_{k=1}^K \pi_k = 1$, $\mathbf{M}_{1:K} = \{\mathbf{M}_1, ..., \mathbf{M}_K\}$, $\sigma_{1:K}^2 = (\sigma_1^2, ..., \sigma_K^2)^T$ and

$$\mathcal{N}_{RG}(\mathbf{X}_i; \mathbf{M}_k, \sigma_k^2) = \frac{1}{\eta(\sigma_k)} \exp \left( - \frac{d_{\text{geo}}^2(\mathbf{X}_i, \mathbf{M}_k)}{2\sigma_k^2} \right), \quad (4)$$

where $\eta(\sigma_k)$ is a normalization factor computed in [31]. As in the Euclidean case, the Riemannian EM algorithm on the manifold of SPD matrices consists in estimating $(\pi_{1:K}, \mathbf{M}_{1:K}, \sigma_{1:K}^2)$ using E and M-steps based on Riemannian tools. For more details, the reader is invited to consult [32]. Note that the algorithms studied in [27] and [28] do not use the LG structure of the data, which is an important difference with the methods proposed in this paper. Note also that the LG structure of the data allows LG integration and derivation operators to be defined in a simple way, which leads to the definition of pdfs on LGs or derivation operators that can be used to define optimization algorithms on LGs.

2) **Clustering on LGs**: classification and clustering on Lie groups is an emerging area of research. However, if supervised methods on Lie groups have received much attention [33], works on unsupervised algorithms are dedicated to specific applications [34]. Mean-shift algorithms have been developed for computer vision problems, particularly for image registration to cluster different kinds of motions [35]. These algorithms extend the standard Euclidean Mean-shift algorithm and are based on maximizing a non-parametric distribution. From a set data $\{\mathbf{X}_i\}_{i=1}^N$, $\mathbf{X}_i$ lying on a LG $G$, the cluster centroids used in mean-shift algorithms can be defined as

$$\widehat{\mathbf{M}}_{1:K} \in \underset{\mathbf{M} \in G}{\operatorname{argmax}} p(\{\mathbf{X}_i\}_{i=1}^N; \mathbf{M}). \quad (5)$$

More precisely, the nonparametric distribution of the data is defined by:

$$p(\{\mathbf{X}_i\}_{i=1}^N; \mathbf{M}) \propto \sum_{i=1}^{N} \kappa (\|\logm(\mathbf{M}^{-1}\mathbf{X}_i)\|^2), \quad (6)$$

where $\kappa: \mathbb{R} \to \mathbb{R}$ is a kernel function. The solutions of (5) can be determined by the following recursion:

$$\mathbf{M}^{(t+1)} = \mathbf{M}^{(t)} \operatorname{expm} (\mathbf{S}^{(t)}), \quad (7)$$

where $\mathbf{S}^{(t)}$ is a mean shift matrix depending on the quantities $\logm((\mathbf{M}^{(t)})^{-1} \mathbf{X}_i)$ for $i = 1, ..., N$.

Fuzzy clustering algorithms have recently been proposed for image segmentation, where a feature vector can be modeled by an affine transformation on $SE(3)$ [36]. Given a set of observations belonging to $SE(3)$, the cluster centroids $\mathbf{M}_1, ..., \mathbf{M}_K$ with $\mathbf{M}_k \in SE(3)$ corresponding to the image classes can be computed by defining weights $w_{i,k}$ for $k = 1, ..., K$ and $i = 1, ..., N$ indicating the probabilities that each pixel belongs to each cluster. The optimization problem solved in [36] is

$$\{\widehat{\mathbf{M}}_k\}_{k=1}^K = \underset{\{\mathbf{M}_k\}_{k=1}^K}{\operatorname{argmin}} \sum_{i=1}^N \sum_{k=1}^K w_{i,k} \|\logm(\mathbf{M}_k^{-1} \mathbf{X}_i)\|_F^2. \quad (8)$$

It is solved by a fuzzy C-Means algorithm on LGs, which leads to closed-form expressions of $\widehat{\mathbf{M}}_1, ..., \widehat{\mathbf{M}}_K$, assuming that

$$\logm(\mathbf{M}_k^{-1} \mathbf{X}_i) \approx \logm(\mathbf{X}_i) - \logm(\mathbf{M}_k). \quad (9)$$

After summarizing these methods, we would like to mention some issues with these algorithms:

* The mean-shift algorithm directly assimilates the maximum of the distribution to the cluster centroid. However, on manifolds, probability density functions (pdfs) do not possess symmetry properties, and the modes differ from the geometric/intrinsic barycenter defined as:

$$\widehat{\mathbf{M}} = \underset{\mathbf{M} \in G}{\operatorname{argmin}} \sum_{i=1}^{N} d^2(\mathbf{M}, \mathbf{X}_i), \quad (10)$$

where $d(\cdot)$ is an LG metric.

* The fuzzy C-Means algorithm of [36] is interesting but does not assume a statistical distribution for the data taking into account their geometry. Moreover, it makes an approximation of (8), suggesting that the Lie group $G = SE(3)$ is commutative, thereby causing a loss of information about the group structure.

This work investigates methods that overcome these issues by designing clustering algorithms that allow for the computation of geometric centroids while incorporating information provided by a statistical model of the data.

## B. Contributions

The contributions of this paper are twofold: first, a K-means algorithm dedicated to LGs (denoted as LG-KM) is defined. This algorithm uses an original intrinsic metric allowing one to determine a similarity measure between a data lying in an LG to the intrinsic centroids of the data clusters. In order to solve the corresponding optimization problem, we propose to update the intrinsic centroid of each cluster by a Gauss-Newton algorithm in LGs. To our knowledge, it is the first time a K-means algorithm is proposed for data lying in LGs.
Second, we derive a new EM algorithm for Gaussian Mixture Models (GMMs) defined in LGs, denoted as LG-GMM-EM. We introduce a theoretical formalism allowing one to intrinsically estimate the LG-mean and LG-covariance of each LG Gaussian component characterizing the barycenter and the shape of non-elliptical clusters. More precisely, the expectation and maximization steps (E-step and M-step) of the EM algorithm are derived by taking into account the properties of distributions defined in LGs. Contrary to the conventional EM algorithm, the M-step of the proposed algorithm involves the computation of LG-means, which is conducted using a non-linear optimization algorithm for LG. In the Euclidean context, the EM algorithm with GMM models can be seen as a probabilistic fuzzy generalization of the soft K-means (fuzzy C-means) algorithm. In the LG context, in addition to better accounting for its geometry through modeling, the proposed LG-GMM-EM generalizes the fuzzy C-means on Lie groups proposed in [36], which was applied to image segmentation.

Compared to the algorithms studied in [27], [28], [36], the proposed algorithms are more generic because they can be applied to any data belonging to LGs such as the space of symmetric positive definite matrices (due to its LG structure [37], [38]), rotation matrices or affine transformations. Furthermore, the considered metric depending on an appropriate LG structure is well-suited to define uncertainties in LGs. The proposed K-means and EM algorithms for LGs are validated numerically in two simulation scenarios. The first scenario considers a synthetic GMM in SE(2). This model allows us to evaluate the algorithm performance in a controlled situation with an available ground truth. In a second step, a realistic application is considered. An ultrasound (US) sensor, lying in an indoor environment, provides several observations coming from different obstacles. One of the challenge in this application is to be able to associate every US observation to an obstacle by taking into account the specific dispersion of the measurements. This problem is classically encountered in the context of bundle adjustment [39] for US sensor calibration [40]. The US sensor is known to have a high uncertainty for its angle measurements. Moreover, spatial positions associated with US measurements are spread out on a curved shape [41]. Consequently, it is challenging to discriminate US observations because the particular geometry of the dispersion is difficult to be accounted with an Euclidean clustering framework. Since this kind of shape is characteristic of the uncertainty cloud of Gaussian distributions in the LG SE(2), we propose to embed US observations in this LG.

The paper is organized as follows: Section II reminds the theoretical background concerning LGs. Sections III and IV introduce the proposed clustering LG-KM and LG-GMM-EM algorithms. Section V is dedicated to experiments conducted to analyze the algorithm performance in two scenarios. In particular, the obtained results are compared with spectral clustering, a well-known state-of-the-art clustering method.

### C. Notations

In the following, $\|\mathbf{a}\|_{\mathbf{A}}$ denotes the Mahalanobis distance of the vector $\mathbf{a} \in \mathbb{R}^p$ with respect to the covariance matrix

$\mathbf{A} \in \mathbb{R}^{p \times p}$, $\mathbf{A}_{1:K}$ denotes the set of matrices $\{\mathbf{A}_k\}_{k=1}^K$ and $|\cdot|$ is the determinant of a matrix. A Lie group is denoted as $G$, $\text{Log}_{\mathcal{G}}^\nu(\cdot)$ is the group logarithm, $\text{Exp}_{\mathcal{G}}^\lambda(\cdot)$ is the group exponential, $\psi_G(\cdot)$ is the inverse left Jacobian of $G$, $l_G(\mathbf{X}, \mathbf{Y}) = \text{Log}_{\mathcal{G}}^\nu(\mathbf{X}^{-1}\mathbf{Y}) \ \forall (\mathbf{X}, \mathbf{Y}) \in G \times G$. Furthermore, $\mathbf{X}_i$ denotes the $i$th observation lying in $G$, $\mathbf{\Theta}$ is the unknown parameter vector of the LG-GMM-EM and $\mathcal{N}_{\mathcal{G}}^L(\cdot)$ denotes the left Gaussian distribution in $G$. Finally, $\mathbb{N}^*$ is the set of non-zero integer values.

## II. BACKGROUND ON LIE GROUPS

This section introduces the different tools necessary to derive the proposed clustering algorithms in LGs.

### A. Definitions

A matrix $LG \subset \mathbb{R}^{r \times r}$ (with $r \in \mathbb{N}^*$) is a set that is endowed with both the properties of a group<fnref n="1" /> and the structure of a smooth manifold [42]. Thus, for any $\mathbf{X} \in G$, a tangent vector space $T_XG$ at $\mathbf{X}$ can be defined. The tangent space at the identity $I_{r \times r}$ plays a specific role and is called the Lie algebra $\mathfrak{g}$. The dimension $m$ of the manifold $G$ is defined as the dimension of the Lie algebra $\mathfrak{g}$. The elements of $\mathfrak{g}$ and $G$ are linked by the exponential map, $\text{Exp}_G : \mathfrak{g} \to G$, and the logarithmic map, $\text{Log}_G : G \to \mathfrak{g}$. There exists an open neighborhood of $0_{r \times r}$ in $\mathfrak{g}$ and an open neighborhood of the identity in $G$ such that the restrictions of the exponential and logarithmic maps to these sets are diffeomorphisms that are inverse transformations (see [42], [43] for more details).

The Lie algebra $\mathfrak{g}$ and the Euclidean space $\mathbb{R}^m$ are linked by a linear isomorphism. If $\mathbf{a}$ is an element of $\mathbb{R}^m$, its image by this isomorphism $\mathbf{a} \in \mathfrak{g}$ is denoted as $\mathbf{a} = [\mathbf{a}]_G^\lambda$ and its inverse is $\mathbf{a} = [\mathbf{a}]_G^\nu$. The operator $[\cdot]^\nu$ is defined such that it transforms the basis elements of the Lie algebra (referred to as generators [44]) into the elements of the canonical basis of $\mathbb{R}^m$. The relationships between the LG, its Lie algebra and the related Euclidean space are shown in Fig. 1. For the sake of simplicity, we introduce the notations $\mathbf{a} = \text{Log}_{\mathcal{G}}^\nu(\mathbf{X}) = [\text{Log}_{\mathcal{G}}(\mathbf{X})]^\nu$ and $\mathbf{X} = \text{Exp}_{\mathcal{G}}^\lambda(\mathbf{a}) = \text{Exp}_{\mathcal{G}}([\mathbf{a}]_G^\lambda)$, with $\mathbf{X} \in G$ and $\mathbf{a} \in \mathbb{R}^m$.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Relationships between $G$, $\mathfrak{g}$ and $\mathbb{R}^m$.</figcaption></figure>

<footnote n="1">¹In this case, the group law reduces to the matrix product.</footnote>
## B. Estimation in Lie groups

Classically, an Euclidean estimator $\hat{\mathbf{x}}$ of the unknown parameter $\mathbf{x} \in \mathbb{R}^p$, constructed from the likelihood $p(\mathbf{z}|\mathbf{x})$ of the observed data vector $\mathbf{z}$, can be characterized by three pertinent statistical quantities<fnref n="2" />: its intrinsic mean $\mathbf{m}_{\hat{\mathbf{x}}}$ such as $\int_{\mathbf{z} \in \mathbb{R}^m} (\hat{\mathbf{x}} - \mathbf{m}_{\hat{\mathbf{x}}})p(\mathbf{z}|\mathbf{x})d\mathbf{z} = 0$, its bias $\int_{\mathbf{z} \in \mathbb{R}^m} (\mathbf{x} - \hat{\mathbf{x}})p(\mathbf{z}|\mathbf{x})d\mathbf{z}$, and its mean square error (MSE) $\int_{\mathbf{z} \in \mathbb{R}^m} (\mathbf{x} - \hat{\mathbf{x}})(\mathbf{x} - \hat{\mathbf{x}})^{\top} p(\mathbf{z}|\mathbf{x})d\mathbf{z}$. Consider a random observation $\mathbf{Z}$ that belongs to an LG denoted as $G'$ whose likelihood function $p(\mathbf{Z}|\mathbf{X})$ depends on an unknown parameter $\mathbf{X} \in G$. An estimator of $\mathbf{X}$ takes the form of a function $\hat{\mathbf{X}}(\mathbf{Z}) : G' \to G$. The error between $\mathbf{X}$ and $\hat{\mathbf{X}}$ can be gauged through the term $||\text{Log}_G^\nabla (\mathbf{X}^{-1}\hat{\mathbf{X}}) ||$ as illustrated in Fig. 2. It is important to emphasize here that this term delineates a geodesic distance for certain LGs, i.e., for semi-simple LGs. However, for the LGs SE(2) and SE(3) (that are considered in various applications), this error does not result from an LG metric. Despite this, it is used as a valuable indicator of the intrinsic trajectory between $\mathbf{X}$ to $\hat{\mathbf{X}}$ and is traditionally employed in the LG estimation literature [13], [15]. For the sake of simplicity, we introduce the notation $l_G(\mathbf{X}, \mathbf{Y}) = \text{Log}_G^\nabla (\mathbf{X}^{-1}\mathbf{Y}) \quad \forall (\mathbf{X}, \mathbf{Y}) \in G \times G$.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Illustration of the intrinsic error between $\mathbf{X}$ and $\hat{\mathbf{X}}$, which takes into account the curvature of an LG.</figcaption></figure>

## C. Optimization in Lie groups

Optimization is the process of identifying the value of a quantity that minimizes or maximizes a specified criterion, typically defined in an Euclidean space. In cases where the optimization problem does not have a closed-form solution, iterative algorithms such as the well-known Newton, Gauss-Newton, or conjugate-gradient methods are classically considered. For the clustering problem considered in this work, we need to minimize cost functions depending on variables belonging to specific LGs. Consequently, it is important to explore optimization techniques that guarantee the sequence of iterates to remain within the admissible space of solutions.

Recent works have given prominence to algorithms that can be used for optimization on manifolds, offering a comprehensive overview in [45]. This interest has led to the development of algorithms specifically designed for LGs. For instance, an adapted version of the Newton algorithm for LGs was

proposed in [46] and [47]. Expanding on this foundation, a Gauss-Newton algorithm customized for LGs was introduced to optimize a quadratic criterion in [14], [48].

Let $\mathbf{X} \in G$ be a variable characterized by an intrinsic dimension $m$, and let $J(\mathbf{X}) : G \to \mathbb{R}$ denote a smooth function on $G$. In order to solve the following optimization problem:

$$ \hat{\mathbf{X}} = \arg \min_{\mathbf{X} \in G} J(\mathbf{X}), \quad (11) $$

the Newton algorithm generates a sequence of iterates $\mathbf{X}_l$ converging towards $\hat{\mathbf{X}}$ under specific conditions on the initial guess and the function $J$. At each iteration $l$, an optimal descent direction $\delta_{l+1} \in \mathbb{R}^m$ is determined by minimizing, with respect to $\delta$, a second-order approximation of $J(\mathbf{X}_l \text{Exp}_G^\wedge(\delta))$. The update step of the algorithm is expressed as follows:

$$ \mathbf{X}_{l+1} = \mathbf{X}_l \text{Exp}_G^\wedge(\delta_{l+1}) \quad \text{with} \quad \delta_{l+1} = -\mathbf{H}_{\hat{\mathbf{X}}_l}^{-1} \nabla J_{\hat{\mathbf{X}}_l}, \quad (12) $$

where

$$ \nabla J_{\mathbf{X}_l} = \frac{\partial(J(\mathbf{X}_l \text{Exp}_G^\wedge(\mathbf{a})))}{\partial \mathbf{a}} \Big|_{\mathbf{a}=0}, \quad (13) $$

$$ \mathbf{H}_{\mathbf{X}_l} = \frac{\partial^2(J(\mathbf{X}_l \text{Exp}_G^\wedge(\mathbf{a}_1) \text{Exp}_G^\wedge(\mathbf{a}_2)))}{\partial \mathbf{a}_1 \partial \mathbf{a}_2} \Big|_{\mathbf{a}_1=0, \mathbf{a}_2=0}. \quad (14) $$

Note that the symbol $\nabla J_{\mathbf{X}_l}$ denotes the gradient of the function $J$ evaluated at the point $\mathbf{X}_l$. This gradient is a column vector of size $m \times 1$, where $m$ represents the dimensionality of the space. Additionally, $\mathbf{H}_{\mathbf{X}_l}$ refers to the geometric Hessian of $J$ computed at $\mathbf{X}_l$, which is a square matrix of size $m \times m$. In the context of this paper, the definition of the gradient aligns with the Lie derivative of the function $J$, while the Hessian corresponds to the second covariant derivative. The computation of (14) involves an affine connection known as the symmetric Cartan-Schouten, as introduced in [47]. For an in-depth understanding of the underlying mathematical concepts, particularly those related to affine connections and Riemannian geometry, the reader can consult [49], [50].

Drawing parallels to the Euclidean scenario, assume that the criterion $J(\mathbf{X})$ can be expressed as the square norm of a vector-valued function:

$$ J(\mathbf{X}) = \frac{1}{2} \|\phi(\mathbf{X})\|_{\mathrm{C}}^{2} = \frac{1}{2} \phi(\mathbf{X})^{\top} \mathbf{C}^{-1} \phi(\mathbf{X}), \quad (15) $$

where $\phi$ is a mapping from $G$ to $\mathbb{R}^d$ with $d \ge m$, and $\mathbf{C}$ is a $d \times d$ positive definite matrix. In this case, it is possible to approximate this mapping using a quadratic form. By using a first-order expansion of $\phi$ around the current estimate at iteration $l$ (i.e., $\mathbf{X}_l$), $J(\mathbf{X})$ can be approximated as:

$$ \begin{aligned} J(\mathbf{X}) & \approx \frac{1}{2} \|\phi(\mathbf{X}_l) + \mathbf{J}_{\phi(\mathbf{X}_l)} (\mathbf{X} - \mathbf{X}_l)\|_{\mathrm{C}}^2 \\ & \approx \frac{1}{2} \phi(\mathbf{X}_l)^{\top} \mathbf{C}^{-1} \phi(\mathbf{X}_l) + \frac{1}{2} \phi(\mathbf{X}_l)^{\top} \mathbf{C}^{-1} \mathbf{J}_{\phi(\mathbf{X}_l)} (\mathbf{X} - \mathbf{X}_l) \\ & \quad + \frac{1}{2} (\mathbf{X} - \mathbf{X}_l)^{\top} \mathbf{J}_{\phi(\mathbf{X}_l)}^{\top} \mathbf{C}^{-1} \mathbf{J}_{\phi(\mathbf{X}_l)} (\mathbf{X} - \mathbf{X}_l). \end{aligned} \quad (16) $$

Note that $\mathbf{J}_{\phi(\mathbf{X}_l)}$ is the LG Jacobian matrix with dimension $d \times m$ of the mapping $\phi$ evaluated at $\mathbf{X}_l$, which can be derived using (13). After differentiating the approximation (16), the

<footnote n="2">²Note that other indicators could be considered (such as the median).</footnote>
Hessian matrix $H_{X_l}$ and the gradient $\nabla J_{X_l}$ of the objective function $J$ at $X_l$ can be approximated as follows:

$$H_{X_l} \approx J_{\phi(X_l)}^{\top} C^{-1} J_{\phi(X_l)}, \nabla J_{X_l} \approx J_{\phi(X_l)}^{\top} C^{-1} \phi(X_l). \quad (17)$$

By substituting these expressions into (12), a Gauss-Newton algorithm for LGs is obtained, denoted as LG-GN algorithm, as detailed in [48].

### D. Uncertainty in Lie groups

To conduct inference in LGs, it is necessary to define distributions on these manifolds. Among those introduced in the literature, our focus is on the concentrated Gaussian distribution (CGD), which can be seen as a generalization of the Euclidean Gaussian multivariate pdf. This distribution has relevant properties, e.g., allowing entropy minimization under certain constraints [51]. Samples from CGDs can be obtained thanks to a change of variable. Let $\epsilon$ be a zero-mean m-dimensional Gaussian vector with covariance matrix $\Sigma$, and $\mu$ a matrix defined on an LG $G$ with intrinsic dimension $m$. The left-CGD pdf of $\mathbf{X} = \mu \exp_G(\epsilon)$ is defined as:

$$p(\mathbf{X}) = \frac{A(\mathbf{X}, \mu)}{\sqrt{(2\pi)^m}} \exp \left( -\frac{1}{2} \|l_G(\mu, \mathbf{X})\|_{\Sigma}^2 \right), \quad (18)$$

with hyperparameters $\mu$ and $\Sigma$ and where $\|\cdot\|_{\Sigma}^2$ denotes the Mahalanobis distance. The multiplicative constant of this pdf is known to be [48]:

$$A(\mathbf{X}, \mu) = |\phi_G(-l_G(\mu, \mathbf{X})) \Sigma \phi_G(-l_G(\mu, \mathbf{X}))^{\top}|^{-\frac{1}{2}},$$

where $|\cdot|$ is the determinant of a matrix and $\phi_G(\cdot)$ is the left Jacobian matrix of $G$ [42] whose expression is given by:

$$\phi_G(\mathbf{a}) = \sum_{n=0}^{+\infty} \frac{ad_G(\mathbf{a})^n}{(n+1)!} \quad \forall \mathbf{a} \in \mathbb{R}^m, \quad (19)$$

with $ad_G(\mathbf{a}) = \mathbf{b} - \mathbf{a} - \mathbf{a}\mathbf{b}, \forall \mathbf{b} \in \mathbb{R}^m$ ($ad_G$ is the adjoint operator of the LG). When $\mathbf{X}$ is sufficiently close to $\mu$, $\phi_G(-l_G(\mu, \mathbf{X}))$ can be approximated by the identity matrix, which allows (18) to be approximated as follows:

$$p(\mathbf{X}) \approx \frac{1}{\sqrt{(2\pi)^m |\Sigma|}} \exp \left( -\frac{1}{2} \|l_G(\mu, \mathbf{X})\|_{\Sigma}^2 \right), \quad (20)$$

$$\triangleq \mathcal{N}_G^L(\mathbf{X}; \mu, \Sigma). \quad (21)$$

We can also define a sample from CGDs with the following formula $\mathbf{X} = \exp_G^\Lambda(\epsilon) \mu$. The resulting right-CGD distribution of $\mathbf{X}$ has pdf:

$$p(\mathbf{X}) \approx \frac{1}{\sqrt{(2\pi)^m |\Sigma|}} \exp \left( -\frac{1}{2} \|l_G(\mathbf{X}^{-1}, \mu^{-1})\|_{\Sigma}^2 \right), \quad (22)$$

$$\triangleq \mathcal{N}_G^R(\mathbf{X}; \mu, \Sigma). \quad (23)$$

### III. K-MEANS ALGORITHM FOR LIE GROUPS

This section studies a new K-Means algorithm for data belonging to an LG. As far as we know, K-means using the LG structure, and not the Riemannian ones, does not yet exist in literature. The main contribution is to define a new metric able to compare and classify two points lying in LGs. Contrary to the clustering algorithms mentioned in

section I-A, the computation of a cluster centroid in this LG is carried out by taking advantage of this metric through an LG optimisation algorithm, based on the second order information of the criterion to optimize. Furthermore, the resulting centroid will correspond to a geometric barycenter due to the inherent structure of the problem.

#### A. Euclidean K-means

The K-means algorithm in an Euclidean space consists of partitioning a set of Euclidean data $\{\mathbf{x}_i\}_{i=1}^N$ into $K$ clusters by computing the centroid associated with each cluster by:

$$\hat{c}_k = \underset{\mathbf{c}_k \in \mathbb{R}^p}{\operatorname{argmin}} \sum_{i \in \mathcal{P}_k} \|\mathbf{x}_i - \mathbf{c}_k\|_{\Sigma_k}^2, \quad (24)$$

where $\mathcal{P}_k$ is the set of index points forming the $k$th cluster of cardinal card($\mathcal{P}_k$), $\|\cdot\|_{\Sigma_k}$ is the Mahalanobis distance, and $\Sigma_k$ is a covariance matrix associated with the $k$th class.<fnref n="3" /> It is known that the centroids satisfying (24) can be found by minimizing the following criterion:

$$J(\mathbf{c}_{1:K}) = \sum_{i=1}^{N} \sum_{k=1}^{K} r_{i,k} \|\mathbf{x}_i - \mathbf{c}_k\|_{\Sigma_k}^2, \quad (25)$$

where $r_{i,k} = 1$ if $i \in \mathcal{P}_k$ and $r_{i,k} = 0$ else. The K-means algorithm searches a minimum of $J$ by using two steps recursively, i.e., by building a set of values $\mathbf{c}_k^{(t)}$ such as, at each iteration $t$:

* The current cluster of each data $\mathbf{x}_i$ is determined

$$k_i^{(t)} = \underset{k=1,...,K}{\operatorname{argmin}} \|\mathbf{x}_i - \mathbf{c}_k^{(t)}\|_{\Sigma_k}^2 \quad \forall i = 1, ..., N, \quad (26)$$

where $\Sigma_k$ is an SPD matrix, allowing the dispersion between the set of observations $\{\mathbf{x}_i\}_{i=1}^N$ and $\mathbf{c}_k^{(t)}$ to be defined.

* The centroid $\mathbf{c}_k^{(t)}$ is updated for each cluster as:

$$\mathbf{c}_k^{(t+1)} = \frac{1}{\operatorname{card}(\mathcal{P}_k^{(t)})} \sum_{i \in \mathcal{P}_k^{(t)}} \mathbf{x}_i, \quad (27)$$

where $\mathcal{P}_k^{(t)} = \{k_i^{(t)}\}$ is the set of index points lying in the $k$th cluster at iteration $t$.

#### B. K-means algorithm in an LG

This section studies a K-means algorithm for data $\{\mathbf{X}_i\}_{i=1}^N$ lying in a matrix LG $G$, with dimension $m$ and scattered in $K$ clusters. By using the LG metric defined in Section II-B, the centroid of each cluster $K$ can be defined as the value of $M_k \in G$ verifying:

$$M_k = \underset{M \in G}{\operatorname{argmin}} \sum_{i \in \mathcal{P}_k} \|l_G(M, \mathbf{X}_i)\|_{\Sigma_k}^2, \quad (28)$$

where $\mathcal{P}_k$ is the set of points lying in the $k$th cluster of the dataset. Note that $M_k$ defines the intrinsic barycenter of the

<footnote n="3">³The K-means algorithm can be defined using different norms. The Mahalanobis distance is used in this work in order to take into account the geometric dispersion of the point clusters.</footnote>
kth cluster. Finding the centroid of each cluster reduces to minimize the following criterion

$$J(\mathbf{M}_{1:K}) = \sum_{i=1}^{N} \sum_{k=1}^{K} \|l_G(\mathbf{M}_k, \mathbf{X}_i)\|_{\Sigma_k}^2, \quad (29)$$

which can be achieved in two steps, for each iteration *t* of the algorithm:

1) First step: As for the Euclidean K-means, the first step determines for each point **X**<sub>*i*</sub>, the closest centroid in the geometric sense:

$$k_i^{(t)} = \underset{k=1,...,K}{\operatorname{argmin}} \|l_G(\mathbf{M}_k^{(t)}, \mathbf{X}_i)\|_{\Sigma_k}^2, \quad (30)$$

where $\Sigma_k$ is an SPD matrix quantifying the dispersion between $\mathbf{M}_k^{(t)}$ and the set of matrices $\{\mathbf{X}_i\}_{i=1}^N$.

2) Second step: The second step computes the values of the new centroids minimizing (28):

$$\mathbf{M}_k^{(t+1)} = \underset{\mathbf{M} \in \mathcal{P}_k^{(t)}}{\operatorname{argmin}} \sum_{i \in \mathcal{P}_k^{(t)}} \|l_G(\mathbf{M}, \mathbf{X}_i)\|_{\Sigma_k}^2 \quad (31)$$

where $\mathcal{P}_k^{(t)}$ is the current set of index points lying at iteration *t* in the *k*th cluster. Contrary to the Euclidean case, the computation of $\mathbf{M}_k^{(t+1)}$ is not trivial, because there is no closed-form expression for the solution of (31). To overcome this problem we can resort to iterative methods using numerical approaches, the minimization problem (31) can be rewritten as:

$$\mathbf{M}_k^{(t+1)} = \underset{\mathbf{M} \in \mathcal{M}}{\operatorname{argmin}} \|\phi_{\text{LG-KM}}(\mathbf{M})\|_{\Sigma_k}^2, \quad (32)$$

with

$$\phi_{\text{LG-KM}}(\mathbf{M}) = [l_G(\mathbf{M}, \mathbf{X}_{F(1)}), \dots, l_G(\mathbf{M}, \mathbf{X}_{F(\text{card}(\mathcal{P}_k^{(t)}))})]^{\top}, \quad (33)$$

and *F* is a function from $\{1, \dots, \text{card}(\mathcal{P}_k)\}$ to $\{1, \dots, N\}$ relating an element of the *k*th cluster to its label $i \in \{1, \dots, N\}$.

We propose to use the Gauss-Newton algorithm for LGs (denoted as LG-GN) to build a set of values $\{\mathbf{M}_k^{(s)}\}_{s=1}^S$ converging to $\mathbf{M}_k^{(t+1)}$, where *S* is fixed by the user such that $\|\phi_{\text{LG-KM}}(\mathbf{M}_k^{(T)})\|^2 - \|\phi_{\text{LG-KM}}(\mathbf{M}_k^{(S)})\|^2 \le \epsilon, \forall T \ge S, \forall \epsilon \ge 0$. Using the recursion (12) with the approximations (17), we propose to compute the LG Jacobian of $l_G(\mathbf{M}_k^{(s)}, \mathbf{X}_i)$ for all indices of the *k*th cluster, i.e., $\frac{\partial l_G(\mathbf{M}_k^{(s)} \text{Exp}_G^\Lambda(\delta), \mathbf{X}_i)}{\partial \delta} \Big|_{\delta=0}$, for $i \in \{F(1), \dots, F(\text{card}(\mathcal{P}_k^{(t)}))\}$. Using the Baker-Campbell-Hausdorff formula [52], the following result is obtained [53]:

$$\frac{\partial l_G(\mathbf{M}_k^{(s)} \text{Exp}_G^\Lambda(\delta), \mathbf{X}_i)}{\partial \delta} \Big|_{\delta=0} = -\psi_G(l_G(\mathbf{M}_k^{(s)}, \mathbf{X}_i)), \quad (34)$$

where $\psi_G(\cdot)$ is the inverse of the left Jacobian $\phi_G(\cdot)$ introduced in II-D. Consequently, the centroid $\mathbf{M}_k^{(s)}$ can be updated by the following recursion

$$\mathbf{M}_k^{(s+1)} = \mathbf{M}_k^{(s)} \text{Exp}_G^\Lambda \left( \underbrace{\delta^{(s+1)}}_{\in \mathbb{R}^m} \right), \forall s = 1, \dots, S, \quad (35)$$

with

$$\delta^{(s+1)} = - \left[ \left( \mathbf{J}^{(s)} \right)^{\top} \Sigma_k^{-1} \mathbf{J}^{(s)} \right]^{-1} \left( \mathbf{J}^{(s)} \right)^{\top} \Sigma_k^{-1} \phi_{\text{LG-KM}}(\mathbf{M}_k^{(s)}), \quad (36)$$

$$\mathbf{J}^{(s)} = \begin{bmatrix} -\psi_G(l_G(\mathbf{X}_{F(1)}, \mathbf{M}_k^{(s)}))^{\top} \\ \vdots \\ -\psi_G(l_G(\mathbf{X}_{F(\text{card}(\mathcal{P}_k^{(t)}))}, \mathbf{M}_k^{(s)}))^{\top} \end{bmatrix}. \quad (37)$$

Note that when the LG is commutative, $\psi_G(\cdot) = \mathbf{I}$ and

$$\mathbf{J}^{(s)} = \underbrace{\begin{bmatrix} -\mathbf{I} \\ \vdots \\ -\mathbf{I} \end{bmatrix}}_{\text{repeated } \text{card}(\mathcal{P}_k^{(t)}) \text{ times}}. \quad (38)$$

The pseudo-codes that can be used to implement the *K*-means algorithm for LGs, referred to as LG-KM, and the dedicated Gauss-Newton algorithm, referred to as LG-GN-KM, are provided in Algorithms 1 and 2.

## IV. PARAMETRIC CLUSTERING ON LGS USING GMMs

This section studies a new parametric clustering model based on GMMs for data belonging to an LG. The idea is to modify the Euclidean GMM to an LG-GMM by using the formalism of the concentrated Gaussian distribution defined in Section II. We propose a new EM algorithm, called LG-GMM-EM, allowing the maximum likelihood of the proposed GMM model to be approximated. It is interesting to highlight the differences between this algorithm and the EM algorithm developed on the Riemannian manifold of covariance matrices introduced in Section I-A. In addition to being less general as it applies to covariance matrices only, the Riemannian algorithm is derived using optimization tools and statistics that differ significantly from those used for Lie groups. The E-step of the proposed LG algorithm is very similar to the E-step of the conventional EM algorithm. However, the M-step allowing the parameters of the LG-GMM to be estimated is performed by estimating the intrinsic mean of each component of the LG-GMM, which exploits the fact that the data belong to an LG.

### A. Gaussian mixture models for LGs

Consider a set of observations $\{\mathbf{X}_i\}_{i=1}^N$ belonging to an LG *G*, that we want to classify in *K* clusters. The LG-GMM algorithm is based on the following pdf:

$$p(\mathbf{X}_i | \mathbf{\Theta}) = \sum_{k=1}^{K} \pi_k \mathcal{N}_G^L(\mathbf{X}_i; \mathbf{M}_k, \mathbf{\Sigma}_k), \quad (39)$$

where $\mathbf{\Theta} = \{\pi_{1:K}, \mathbf{M}_{1:K}, \mathbf{\Sigma}_{1:K}\}$ contains the unknown parameters of the GMM, i.e., $\mathbf{M}_k \in G$ and $\mathbf{\Sigma}_k \in \mathbb{R}^{m \times m}$ are the LG mean and covariance of the *k*th LG-Gaussian distribution, $\pi_k \in [0, 1]$ is the prior probability for the *k*th cluster such that $\sum_{k=1}^K \pi_k = 1$. Assuming independence
between the observations {$\mathbf{X}_i$}$_{i=1}^N$, the joint density of the GMM is obtained:

$$p(\mathbf{X}|\mathbf{\Theta}) \triangleq \prod_{i=1}^{N} p(\mathbf{X}_i|\mathbf{\Theta}), \quad (40)$$

leading to the following log-likelihood:

$$\log p(\mathbf{X}|\mathbf{\Theta}) = \sum_{i=1}^{N} \log \sum_{k=1}^{K} \pi_k \mathcal{N}_G^L(\mathbf{X}_i; \mathbf{M}_k, \mathbf{\Sigma}_k). \quad (41)$$

It is worth noting that this mixture model could be easily defined by using a right Gaussian distribution as defined in (23). However, the left formalism leads to simpler formulations.

## B. Maximum likelihood estimator: LG-GMM-EM algorithm

The maximum likelihood estimator (MLE) of the GMM for LGs introduced in the previous section is the value of $\mathbf{\Theta}$ that maximizes (41), i.e.,

$$\hat{\mathbf{\Theta}}_{\text{MLE}} = \underset{\mathbf{\Theta}}{\operatorname{argmax}} \sum_{i=1}^{N} \log \sum_{k=1}^{K} \pi_k \mathcal{N}_G^L(\mathbf{X}_i; \mathbf{M}_k, \mathbf{\Sigma}_k). \quad (42)$$

This maximization being difficult to be done analytically, we propose to investigate an EM algorithm dedicated to GMMs on LGs referred to as LG-GMM-EM.

1) **Estimation problem:** Similarly to the Euclidean framework, we introduce a discrete latent variable $z_i \in \{1, \dots, K\}$ [54] indicating to which cluster the ith observation belongs with $\pi_k \triangleq \mathbb{P}(z_i = k | \mathbf{\Theta})$. The EM algorithm for the GMM (42) maximizes the function $Q(\mathbf{\Theta}, \mathbf{\Theta}^{(t)})$ defined by:

$$Q(\mathbf{\Theta}, \mathbf{\Theta}^{(t)}) = \mathbb{E}_{z_i \mid \mathbf{\Theta}^{(t)}, \mathbf{X}}[\log p(\mathbf{z}, \mathbf{X} | \mathbf{\Theta})], \quad (43)$$

$$\triangleq \sum_{k=1}^{K} \sum_{i=1}^{N} \gamma_{i,k}(\mathbf{\Theta}^{(t)}) \log \mathbb{P}(z_i = k, \mathbf{X}_i | \mathbf{\Theta}), \quad (44)$$

where $\mathbf{\Theta}^{(t)} = \{\pi_{1:K}^{(t)}, \mathbf{M}_{1:K}^{(t)}, \mathbf{\Sigma}_{1:K}^{(t)}\}$ is the current value of the GMM parameters at iteration t, and $\gamma_{i,k}(\mathbf{\Theta}^{(t)}) = \mathbb{P}(z_i = k | \mathbf{\Theta}^{(t)}, \mathbf{X}_i)$. The function $Q(\mathbf{\Theta}, \mathbf{\Theta}^{(t)})$ is maximized with respect to $\mathbf{\Theta}$ to obtain the current parameter at iteration t + 1:

$$\mathbf{\Theta}^{(t+1)} = \underset{\mathbf{\Theta}}{\operatorname{argmax}} Q(\mathbf{\Theta}, \mathbf{\Theta}^{(t)}). \quad (45)$$

The E-step and M-step of the proposed EM algorithm are detailed below.

2) **E-step:** The E-step consists in computing $\gamma_{i,k}(\mathbf{\Theta}^{(t)}) = \mathbb{P}(z_i = k | \mathbf{\Theta}^{(t)}, \mathbf{X}_i)$. The Bayes rule leads to

$$\gamma_{i,k}(\mathbf{\Theta}^{(t)}) = \frac{p(\mathbf{X}_i|z_i = k, \mathbf{\Theta}^{(t)})\mathbb{P}(z_i = k|\mathbf{\Theta}^{(t)})}{\sum_{k=1}^{K} p(\mathbf{X}_i|z_i = k, \mathbf{\Theta}^{(t)})\mathbb{P}(z_i = k|\mathbf{\Theta}^{(t)})}, \quad (46)$$

The density $p(\mathbf{X}_i|z_i = k, \mathbf{\Theta}^{(t)})$ is the distribution of $\mathbf{X}_i$ assuming that it lies in the cluster k:

$$p(\mathbf{X}_i|z_i = k, \mathbf{\Theta}^{(t)}) = p(\mathbf{X}_i|\mathbf{M}_k^{(t)}, \mathbf{\Sigma}_k^{(t)}) = \mathcal{N}_G^L(\mathbf{X}_i; \mathbf{M}_k^{(t)}, \mathbf{\Sigma}_k^{(t)}). \quad (47)$$

Since

$$\mathbb{P}(z_i = k | \mathbf{\Theta}^{(t)}) = \pi_k^{(t)}, \quad (48)$$

we obtain

$$\gamma_{i,k}(\mathbf{\Theta}^{(t)}) = \frac{\mathcal{N}_G^L(\mathbf{X}_i; \mathbf{M}_k^{(t)}, \mathbf{\Sigma}_k^{(t)}) \pi_k^{(t)}}{\sum_{k=1}^{K} \mathcal{N}_G^L(\mathbf{X}_i; \mathbf{M}_k^{(t)}, \mathbf{\Sigma}_k^{(t)}) \pi_k^{(t)}}. \quad (49)$$

3) **M-step:** The M-step consists of maximizing (44) with respect to $\mathbf{\Theta}$, i.e.,

$$\mathbf{\Theta}^{(t+1)} = \underset{\mathbf{\Theta}}{\operatorname{argmax}} Q(\mathbf{\Theta}, \mathbf{\Theta}^{(t)}), \quad (50)$$

which needs to compute $\log \mathbb{P}(z_i = k, \mathbf{X}_i | \mathbf{\Theta})$. Using the definition of conditional distributions, one obtains:

$$\mathbb{P}(z_i = k, \mathbf{X}_i | \mathbf{\Theta}) = p(\mathbf{X}_i | z_i = k, \mathbf{\Theta}) \underbrace{\mathbb{P}(z_i = k | \mathbf{\Theta})}_{\pi_k}. \quad (51)$$

Using the expression of the left CGD pdf (20), we obtain:

$$\log \mathbb{P}(z_i = k, \mathbf{X}_i | \mathbf{\Theta}) = -\frac{1}{2} \|l_G(\mathbf{M}_k, \mathbf{X}_i)\|_{\Sigma_k}^2 \\ -\frac{1}{2} \log (|\Sigma_k|) + \log \pi_k \quad (52)$$

Substituting (52) in (44) yields the following result:

$$Q(\mathbf{\Theta}, \mathbf{\Theta}^{(t)}) = \sum_{k=1}^{K} \sum_{i=1}^{N} \gamma_{i,k}(\mathbf{\Theta}^{(t)}) \log(\pi_k) \\ - \frac{1}{2} \sum_{k=1}^{K} \sum_{i=1}^{N} \{\log (|\Sigma_k|) + \|l_G(\mathbf{M}_k, \mathbf{X}_i)\|_{\Sigma_k}^2\}. \quad (53)$$

Note that the computation of $Q$ is similar to the Euclidean algorithm, since the classical rules of probability can be applied to LGs. Particularly, we observe that the determination of $Q$ reduces to compute the “responsabilities” $\gamma_{i,k}(\mathbf{\Theta}^{(t)})$. However, contrary to the conventional EM algorithm, it is not possible to obtain a closed-form expression of the LG-mean estimators $\{\tilde{\mathbf{M}}_k\}_{k=1}^K$ from (53). An outstanding result is that, conditionally to $\mathbf{M}_k$, it is possible to derive analytical expressions of $\{\mathbf{\Sigma}_k\}_{k=1}^K$ and $\{\pi_k\}_{k=1}^K$ maximizing $Q(\mathbf{\Theta}, \mathbf{\Theta}^{(t)})$, which allows the optimization of $Q(\mathbf{\Theta}, \mathbf{\Theta}^{(t)})$ to be simplified significantly. Based on the previous comments, in order to solve (50), we propose to initialize the matrices $\mathbf{M}_k$, $k = 1, \dots, K$ by $\mathbf{M}_k^{(t)}$ and to repeat the following two steps until the differences between consecutive estimates is less than $10^{-1}$:

1) For fixed matrices $\mathbf{M}_k$ resulting from the previous iteration, compute the values of $\mathbf{\Sigma}_k$ and $\pi_k$ obtained by cancelling the Euclidean gradients of $Q(\mathbf{\Theta}, \mathbf{\Theta}')$ under the constraint $\sum_{k=1}^K \pi_k = 1$, which leads to:

$$\Sigma_k^{(t+1)} = \frac{\sum_{i=1}^{N} \gamma_{i,k}(\mathbf{\Theta}^{(t)}) l_G(\mathbf{M}_k^{(t)}, \mathbf{X}_i) l_G(\mathbf{M}_k^{(t)}, \mathbf{X}_i)^\top}{\sum_{i=1}^{N} \gamma_{i,k}(\mathbf{\Theta}^{(t)})}, \\ \pi_k^{(t+1)} = \frac{\sum_{i=1}^{N} \gamma_{i,k}(\mathbf{\Theta}^{(t)})}{\sum_{i=1}^{N} \sum_{k=1}^{K} \gamma_{i,k}(\mathbf{\Theta}^{(t)})}. \quad (54)$$
2) Replace the values of $\Sigma_k$ in $Q(\Theta, \Theta^{(t)})$ by the values of $\Sigma_k^{(t+1)}$ obtained in step 1) and determine the values of $M_{1:K}$ maximizing the resulting criterion, i.e., solve the following LG-optimization problem:

$$ M_{1:K}^{(t+1)} = \underset{M_{1:K}}{\operatorname{argmin}} \sum_{k=1}^{K} \sum_{i=1}^{N} \|l_G(M_k, X_i)\|_{\Sigma_k^{(t+1)} \gamma_{i,k}(\Theta^{(t)})}^2 \quad (55) $$

This minimization is complicated due to the non linearity of $\mathrm{Log}_G^\nabla(\cdot)$. In order to use the Gauss-Newton algorithm on LGs, we propose to rewrite the problem with the following compact form:

$$ M_{1:K}^{(t+1)} = \underset{M_{1:K}}{\operatorname{argmin}} \|\phi_{\mathrm{LG-GMM-EM}}(M_{1:K})\|_{\Sigma_N^{(t+1)}}^2, \quad (56) $$

with:

$$ \phi_{\mathrm{LG-GMM-EM}}(\mathbf{M}_{1:K}) = \begin{bmatrix} l_G(\mathbf{M}_1, \mathbf{X}_1) \sqrt{\gamma_{1,1}(\mathbf{\Theta}^{(t)})} \\ \vdots \\ l_G(\mathbf{M}_K, \mathbf{X}_N) \sqrt{\gamma_{N,K}(\mathbf{\Theta}^{(t)})} \end{bmatrix}, \quad (57) $$

$$ \Sigma_N^{(t+1)} = \begin{bmatrix} \Sigma^{(t+1)} & \cdots & 0 \\ \vdots & \ddots & 0 \\ 0 & \cdots & \Sigma^{(t+1)} \end{bmatrix}, \quad (58) $$

$$ \Sigma^{(t+1)} = \begin{bmatrix} \Sigma_1^{(t+1)} & \cdots & 0 \\ \vdots & \ddots & 0 \\ 0 & \cdots & \Sigma_K^{(t+1)} \end{bmatrix}, \quad (59) $$

where the notation $\Sigma_N^{(t+1)}$ has been used to highlight the fact that this matrix has $N$ diagonal blocks equal to the same matrix $\Sigma^{(t+1)}$. Thus, the optimization of the LG cluster means can be performed with an LG-GN iteration in the LG $G' = G \times \dots \times G = G^K$ in order to estimate the parameter

$$ M = \begin{bmatrix} \mathbf{M}_1 & \cdots & 0 \\ \vdots & \ddots & 0 \\ 0 & \cdots & \mathbf{M}_K \end{bmatrix}. \quad (60) $$

The resulting recursions of the Gauss-Newton algorithm at iteration $s$ are

$$ \mathbf{P}^{(s)} = \mathbf{J}^{(s)} (\Sigma_N^{(t+1)})^{-1} \mathbf{J}^{(s)}, \quad (61) $$

$$ \delta^{(s)} = (\mathbf{P}^{(s)})^{-1} (\mathbf{J}^{(s)})^{\top} (\Sigma_N^{(t+1)})^{-1} \phi_{\mathrm{LG-GMM-EM}}(\mathbf{M}^{(s)}), \quad (62) $$

$$ \mathbf{M}^{(s+1)} = \mathbf{M}^{(s)} \mathrm{Exp}_{G'}^{\Lambda} (-\delta^{(s)}), \quad (63) $$

with $\delta^{(s)} \in \mathbb{R}^{mK}$,

$$ \mathbf{J}^{(s)} = [(\mathbf{J}^{(1)})^{\top}, \dots, (\mathbf{J}^{(N)})^{\top}]^{\top}, \quad (64) $$

and $\mathbf{J}^{(i)}$ is the following $K \times K$ block diagonal matrix:

$$ \mathbf{J}^{(i)} = \begin{bmatrix} \mathbf{A}_{1,i} & \cdots & 0 \\ \vdots & \ddots & 0 \\ 0 & \cdots & \mathbf{A}_{K,i} \end{bmatrix}, \quad (65) $$

whose elements are $\mathbf{A}_{k,i} = \psi_G(l_G(\mathbf{M}_k^{(s)}, \mathbf{X}_i)) \sqrt{\gamma_{i,k}(\mathbf{\Theta}^{(t)})}$, for $k = 1, \dots, K$, leading to Algorithms 3 and 4.

## V. SIMULATIONS

This section evaluates the performance of the proposed K-means and EM algorithms using two kinds of experiments. The first experiments are conducted with observations belonging to a mixture of LGs that we want to cluster with an unsupervised method. The numerical consistency and precision of the proposed algorithms LG-KM and LG-GMM-EM is assessed in terms of clustering performance. The second scenario considers a tracking application where observations are acquired by an ultrasound sensor. This sensor measures its distance and angle between two static objects at different instants. The problem then consists of separating the measurements into two clusters associated with the two objects.

### A. First scenario: synthetic LG mixture model

The first experiments consider a set of observations lying in the LG $SE(2)$ corresponding to the set of 2D rigid-body transformations. More precisely, each measurement can be characterized by the following parameters

$$ \mathbf{X} = \begin{bmatrix} \mathbf{R} & \mathbf{p} \\ \mathbf{0} & 1 \end{bmatrix}, \mathbf{R} \in SO(2), \mathbf{p} \in \mathbb{R}^2, \quad (66) $$

where $\mathbf{R}$ is a rotation matrix and $\mathbf{p}$ is a translation vector. A set of $N$ observations simulated according to a mixture of Gaussian distributions in $SE(2)$ is then considered:

$$ \mathbf{X}_i | \mathbf{M}_k, \mathbf{\Sigma}_k \sim \sum_{k=1}^{K} \pi_k \mathcal{N}_{SE(2)}^{L}(\mathbf{M}_k, \mathbf{\Sigma}_k), \quad (67) $$

with $i = 1, \dots, N$. The parameters of the $K = 3$ different clusters are summarized in Table I. In this scenario, we want to generate data having a curved shape in each cluster, in order to observe the benefits of the LG modelling. Since each cluster is controlled by the coefficients of the matrix $\Sigma_k$, we propose the following process to build it. First, we generate Gaussian vectors $\mathbf{E}_i \sim N(0, \mathrm{diag}(\boldsymbol{\sigma}_k))$ with $\boldsymbol{\sigma}_k$ the cluster variance vector modelling the curvature of the shape. A vector $\boldsymbol{\mu}$ modelling the size of the cluster is then defined. The covariance matrices are then constructed as follows:

$$ \mathbf{\Sigma}_k \triangleq \mathbf{\Sigma} = \frac{1}{N} \sum_{i=1}^{N} l_{SE(2)}(\boldsymbol{\mu} \cdot \mathbf{E}_i) l_{SE(2)}(\boldsymbol{\mu}, \mathbf{E}_i)^{\top}, \quad (68) $$

where $l_{SE(2)}$ corresponds to the operator $l_G$ when $G = SE(2)$.

<table><caption>TABLE I</caption><thead><tr><th>Cluster means</th><th></th></tr></thead><tbody><tr><td>M<sub>1</sub> = Exp<sup>Λ</sup><sub>SE(2)</sub> ([π/4, 2, 1])</td><td></td></tr><tr><td>M<sub>2</sub> = Exp<sup>Λ</sup><sub>SE(2)</sub> ([0.1, 0, 1])</td><td></td></tr><tr><td>M<sub>3</sub> = Exp<sup>Λ</sup><sub>SE(2)</sub> ([0.2, 0, 0])</td><td></td></tr><tr><td>Number of clusters</td><td>K = 3</td></tr><tr><td>Cluster variance</td><td>σ = [0.6, 1, 1]<sup>T</sup></td></tr><tr><td>Cluster size</td><td>μ = Exp<sup>Λ</sup><sub>SE(2)</sub> ([0, 1, 1]<sup>T</sup>)</td></tr><tr><td>Cluster weights</td><td>π<sub>1</sub>, π<sub>2</sub>, π<sub>3</sub> = 1/3</td></tr><tr><td>Number of points</td><td>N = 300 (100 points per cluster)</td></tr></tbody></table>

SIMULATION PARAMETERS FOR THE FIRST EXPERIMENT.
In an applicative context, the three point clouds can correspond to the positions of three objects moving in a turn (e.g., a car or an autonomous robot). Each point cloud has a banana-shape whose degree of curvature is characterized by the covariance matrix $\Sigma_k$, as illustrated in Fig. 3.

<figure><img src="image_4.png" /><figcaption>Fig. 3. Example of observations generated according to (67) with the simulation parameters of Table I.</figcaption></figure>

The first simulations analyze the performance of the LG-KM algorithm to cluster the three point clouds. In order to initialize this algorithm, we assume that the cluster centroids $\{M_k^{(0)}\}_{k=1}^K$ are randomly localized in the region of interest. Fig. 4 shows the initialization and the evolution of the clusters along the iterations. The estimated cluster centroids are very close to the true centroids after $n \sim 14$ iterations only. Furthermore, the algorithm is able to classify every point by taking into account the cluster geometry. Fig. 5 displays the evolution of the criterion (29) versus iterations, confirming a fast convergence of the algorithm.

<figure><img src="image_5.png" /><figcaption>Fig. 5. Evolution of the K-means criterion (29) versus iterations.</figcaption></figure>

The next experiments study the performance of the LG-GMM-EM algorithm for a mixture of Gaussian LG data. The simulation parameters were chosen as in Table I with random initializations of $M_k^{(0)}$ and initial covariance matrices equal to the identity matrix. The initial values of the weights are $\pi_k^{(0)} = 1/3$. To visualize the clustering performance, we display the evolution of the clusters at different iterations by computing a confidence shape in a geometrical way.

<figure><img src="image_6.png" /><figcaption>(a) Initialization.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(c) Iteration 4.</figcaption></figure>

<caption>(d) Iteration 14.</caption>

<figure><img src="image_8.png" /><figcaption>(a) Initial iteration of LG-GMM-EM.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(b) Iteration 4 of LG-GMM-EM.</figcaption></figure>

<caption>(c) Iteration 14 of LG-GMM-EM.</caption>

<caption>Fig. 6. Evolution of the LG-GMM-EM clustering along the iterations.</caption>
<figure><img src="image_10.png" /><figcaption>Fig. 7. Evolution of the log-likelihood versus the number of iterations.</figcaption></figure>

As for the conventional EM algorithm, a confidence region can be built using the Mahalanobis distance of the data. More precisely, we define a geometrical Mahalanobis distance allowing us to assign the observations $\mathbf{X}_i$ to one of the clusters and to determine a confidence shape specific to the SE(2) geometry. At each iteration $l$, the distance between the point $\mathbf{X}_i$ and the centroid $\mathbf{M}_k^{(t)}$ is defined by:

$$d(\mathbf{M}_k^{(t)}, \mathbf{X}_i; \mathbf{\Sigma}_k^{(t)}) = l_{SE(2)}(\mathbf{M}_k^{(t)}, \mathbf{X}_i)^{\top} (\mathbf{\Sigma}_k^{(t)})^{-1} l_{SE(2)}(\mathbf{M}_k^{(t)}, \mathbf{X}_i). \quad (69)$$

The confidence shape for the $k$th cluster containing all the points $\mathbf{X}_i$ of the $k$th class is such that:

$$d(\mathbf{M}_k^{(t)}, \tilde{\mathbf{X}}_i; \mathbf{\Sigma}_k^{(t)}) \le \gamma, \quad (70)$$

where $\gamma$ is a pre-defined threshold which is set empirically.

Fig. 6 shows the confidence shapes of each cluster superimposed to their estimated centroids, for different iterations of the LG-GMM-EM algorithm. The true point clouds can be recovered after 10 iterations for this example. The estimated centroids are also approaching the true centroids along the iterations. Finally, one can observe that every point cloud associated with a given cluster is located inside the corresponding confidence shape. This visual behaviour can be confirmed by plotting the log-likelihood (41) versus the number of iterations. Figure 7 shows an almost linear variation of the log-likelihood during the first iterations and a stabilization after 10 iterations. Finally, it is worth noting that the LG-GMM-EM algorithm converges slightly faster than the LG-KM algorithm (10 iterations against 14 iterations). This result is consistent because the observation model is theoretically adapted to the formalism of the LG-GMM-EM through the mixture model (39).

The LG-GMM-EM algorithm has a better precision than LG-KM, at the price of a higher computational complexity, as shown in Table II. This is due to the computation of $\delta$ in the Gauss-Newton algorithm requiring the inversion of the $3NK \times 3NK$ matrix $\tilde{\Sigma}_N$ used in (58). Note that the matrix $\Sigma'_N$ is not used in the LG-KM algorithm.

<table><caption>TABLE II COMPUTATIONAL COSTS FOR LG-KM AND LG-GMM-EM ALGORITHMS.</caption><thead><tr><th>Algorithm</th><th>Average running time for one iteration (s)</th></tr></thead><tbody><tr><td>LG-KM</td><td>1.3</td></tr><tr><td>LG-GMM-EM</td><td>5.2</td></tr></tbody></table>

## B. Second scenario: application to ultrasound data

1) Context: An US sensor located in an indoor environment sends an US wave that is backscattered to the sensor after reaching one of several objects (supposed to be static) located in the region of interest. We assume that several measurements of the same objects are available, which means that each US wave is reflected several times $t_i$, $i = 1, \dots, N$ before reaching the sensor. As illustrated in Fig. 8, the US sensor recovers a reflected signal from which we can extract several quantities:

* the distance $d_i^{(k)}$ between the sensor with position $\mathbf{p}_r$ and the $k$th object with position $\mathbf{p}^{(k)}$:

$$d_i^{(k)} = \|\mathbf{p}^{(k)} - \mathbf{p}_r\| + \epsilon_{d,i} \quad \epsilon_{d,i} \sim \mathcal{N}(0, \sigma_d^2), \quad (71)$$

* the incidence angle $\theta_i^{(k)}$ between the sensor and the $k$th object:

$$\theta_i^{(k)} = \operatorname{atan}\left(\frac{(\mathbf{p}^{(k)} - \mathbf{p}_r)_1}{\|\mathbf{p}^{(k)} - \mathbf{p}_r\|}\right) + \epsilon_{\theta,i} \quad \epsilon_{\theta,i} \sim \mathcal{N}(0, \sigma_{\theta}^2). \quad (72)$$

Note that the 2D Cartesian positions of the objects ($x_i^{(k)}, y_i^{(k)}$) can be recovered from these distances and angles since

$$x_i^{(k)} = d_i^{(k)} \cos(\theta_i^{(k)}), \quad y_i^{(k)} = d_i^{(k)} \sin(\theta_i^{(k)}). \quad (73)$$

<figure><img src="image_11.png" /><figcaption>Fig. 8. Ultrasound sensor measuring distances and angles of objects yielding object positions.</figcaption></figure>

Due to the sensor properties, the uncertainties on the measured angles $\{\theta_i\}_{i=1}^N$ are large, which induces measured positions that are spread according to a curved shape as illustrated in Fig. 9 for $K = 2$ objects with the following simulation parameters:

<table><caption>TABLE III SIMULATION PARAMETERS FOR THE SECOND SCENARIO.</caption><tr><td>K</td><td>2</td></tr><tr><td>σ<sub>θ</sub></td><td>0.6 rad</td></tr><tr><td>σ<sub>d</sub></td><td>1 m</td></tr><tr><td>P<sub>1</sub></td><td>7; 7 m</td></tr><tr><td>P<sub>2</sub></td><td>9; 9 m</td></tr><tr><td>P<sub>r</sub></td><td>0; 0 m</td></tr></table>
<figure><img src="image_12.png" /><figcaption>Fig. 9. Point cloud of $N = 30$ US measurements for $K = 2$ detected objects.</figcaption></figure>

When several temporal measurements are detected by the US sensor, it is challenging to determine the origin of each measurement because the point clouds are curved due to the geometry of this problem. Assigning the different measurements to the objects is fundamental for bundle adjustment, e.g., in the context of Simultaneous Localization And Mapping (SLAM) [55]. Indeed, in this context, data are collected temporally from different unknown landmarks. Thus, it is crucial to associate the measurements to a landmark to estimate the robot position.

To take into account the specific shape of the observed point clouds {$\theta_i^{(k)}, x_i^{(k)}, y_i^{(k)}$}$_{i=1}^N$, we propose to embed them in the LG SE(2) in order to create the observation $\mathbf{X}_{i,k}$:

$$ \mathbf{X}_{i,k} = \begin{bmatrix} \cos(\theta_i^{(k)}) & -\sin(\theta_i^{(k)}) & x_i^{(k)} \\ \sin(\theta_i^{(k)}) & \cos(\theta_i^{(k)}) & y_i^{(k)} \\ 0 & 0 & 1 \end{bmatrix} \quad (74) $$

where $\mathbf{X}_{i,k}$ belongs to the $k$th cluster.

In order to cluster the points of Fig. 9, we perform the LG-KM but also the LG-GMM-EM algorithm with a right Gaussian distribution formalism. Indeed, the right distribution is better suited for fitting curved shapes because it naturally provides curved shape through the right product between the mean and the noise as follows

$$ \mathbf{X}_{i,k} = \operatorname{Exp}^{\wedge}_{SE(2)} (\epsilon_i) \mathbf{M}_k, \quad \epsilon_i \sim \mathcal{N}(\mathbf{0}, \Sigma_k), \quad (75) $$

where {$\mathbf{M}_k, \Sigma_k$} $\in$ SE(2) $\times$ $\mathbb{R}^{p \times p}$ are the unknown parameters.

Both algorithms are initialized with $K = 2$ and:

$$ \mathbf{M}_1^{(0)} = \begin{bmatrix} \operatorname{Exp}^{\wedge}_{SO(2)} (\alpha_1) & \mathbf{0}_{2 \times 1} \\ \mathbf{0}_{1 \times 2} & 1 \end{bmatrix}, \quad (76) $$

$$ \mathbf{M}_2^{(0)} = \begin{bmatrix} \operatorname{Exp}^{\wedge}_{SO(2)} (\alpha_2) & \mathbf{0}_{2 \times 1} \\ \mathbf{0}_{1 \times 2} & 1 \end{bmatrix}, \quad (77) $$

with $\alpha_1 = 0$, $\alpha_2 = 0.4$. The LG-GMM-EM algorithm was run with initial covariance matrices $\Sigma_k^{(0)} = \mathbf{I}$ and initial weights $\pi_k^{(0)} = 1/2$. Figs. 10 and 11 illustrate the good clustering performance of the LG-KM algorithm. Indeed, the centroids of the different clusters converge towards the unknown natural intrinsic barycenters of the two curved shapes. This convergence implies an excellent clustering after $l = 10$ iterations only, with only two errors in the whole point cloud. This good performance is confirmed by the evolution of the log-likelihood versus the number of iterations in Fig. 12.

<figure><img src="image_13.png" /><figcaption>Fig. 10. Initial iteration of LG-KM.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 11. Evolution of the LG-KM clustering along the iterations.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 12. Evolution of the log-likelihood (41) versus the number of iterations.</figcaption></figure>

Figs. 13 and 14 display the evolution of the LG-GMM-EM clusters versus the number of iterations by drawing confidence shapes. The cluster shapes are well-reconstructed as well as the intrinsic centroids after $l = 9$ iterations only. The convergence speed is slightly faster for LG-GMM-EM than for LG-KM, as previously. The variation of the log-likelihood of the data versus the number of iterations displayed in Fig. 15 shows a fast convergence to a stable value.
<figure><img src="image_16.png" /><figcaption>Fig. 13. Initial iteration of LG-GMM-EM</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 14. Evolution of the LG-GMM-EM clustering along the iterations</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 15. Evolution of the log-likelihood (41) versus the number of iterations.</figcaption></figure>

2) *Comparison with spectral clustering:* In order to validate the proposed LG-GMM-EM algorithm for a mixture of distributions in an LG, this section compares its clustering performance with spectral clustering [56]. Spectral clustering (SC) is one of the benchmarks for unsupervised clustering, which was developed for non-elliptic and non-linear point clouds. The SC algorithm can be summarized as follows:

1) build a similarity graph $\mathbf{W} = \{w_{i,j}\}_{(i,j)\in[1:N]^2}$ and the associated vertex matrix $\mathbf{D} = \text{diag}(d_1, ..., d_N)$, with $d_i = \sum_{j=1}^N w_{i,j}$, $i = 1, ..., N$ (related to the distances between the observed vectors),

2) compute the Laplacian matrix $\mathbf{L} = \mathbf{W} - \mathbf{D}$ and its eigenvalue decomposition,

3) cluster the $N$ rows of the previous eigenvector matrix thanks to a conventional K-means algorithm.

Note that the curvature of the dataset modelled by (71) and (72) is directly related to the standard deviation $\sigma_{\theta}$. The higher $\sigma_{\theta}$, the higher this curvature, the relation between these two quantities being highly non-linear. Thus, the LG-modelling seems to be relevant to cluster this dataset displayed in Fig. 9. Fig. 16 shows the classification error rate of this dataset as a function of $\sigma_{\theta}$. This rate is globally constant for the LG-GMM-EM method that is able to learn the particular shape of both clusters. Conversely, the performance of the SC method decreases when the curvature of the point cloud increases. Note that for high value of $\sigma_{\theta}$, the algorithm classifies one point out of two in the wrong cluster. To conclude, SC is not always able to separate highly curved clusters contrary to the proposed LG-GMM-EM method.

<figure><img src="image_19.png" /><figcaption>Fig. 16. Evolution of the classification error rate versus $\sigma_{\theta}$.</figcaption></figure>

## VI. CONCLUSIONS AND PERSPECTIVES

This paper introduced and studied two new clustering algorithms for data lying in Lie groups (LGs) extending the classical K-Means and EM algorithms. These new clustering algorithms exploited an intrinsic metric related to the LG structure and were based on mixture of Gaussian distributions in LGs. The optimization algorithm proposed to estimate the unknown model parameters for K-means and EM was based on a Gauss-Newton algorithm adapted to data belonging to LGs. Both algorithms were validated by numerical simulations showing better clustering performance than existing conventional clustering algorithms. The perspectives of these works are numerous. A direct application would be to extend the problem to a 3D scenario, i.e., when $G = SE(3)$. The proposed clustering algorithm could be also adapted to statistical models with missing variables belonging to LGs, which requires to define a modified LG-GMM-EM algorithm with latent variables lying in LGs. Also, it could be relevant to adapt the proposed algorithms to cluster covariance matrices benefiting from their Lie group structure and the log-Euclidean metric.

## REFERENCES

[1] S. M. A. Bhuiyan and J. F. Khan, “A Simple Deep Learning Network for Target Classification,” in *Proc. SoutheastCon 2019*, Huntsville, AL, USA, April 2019, pp. 1–5.

[2] W. Jiang, Y. Ren, Y. Liu, and J. Leng, “Artificial Neural Networks and Deep Learning Techniques Applied to Radar Target Detection: A Review,” *Electronics*, vol. 11, no. 1, p. 156, 2022.

[3] J. B. MacQueen, “Some Methods for Classification and Analysis of Multivariate Observations,” in *Proc. 5th Berkeley Symp. Math. Stat. Probability*, vol. 1, Berkeley, CA, USA, Jun.-Jul. 1965, pp. 281–297.
[4] E. Niharika, H. Adeeba, A. S. R. Krishna, and P. Yugander, "K-means based noisy sar image segmentation using median filtering and otsu method," in *2017 International Conference on IoT and Application (ICIOT)*, 2017, pp. 1-4.

[5] S. Moss and E. R. Hancock, "Registering incomplete radar images using the em algorithm," *Image and Vision Computing*, vol. 15, no. 8, pp. 637-648, 1997, British Machine Vision Conference.

[6] A. Raskaliyev, S. H. Patel, T. M. Sobh, and A. Ibrayev, "GNSS-Based Attitude Determination Techniques—A Comprehensive Literature Survey," *IEEE Access*, vol. 8, pp. 24873-24886, 2020.

[7] M. Arnaudon, F. Barbaresco, and L. Yang, "Riemannian Medians and Means With Applications to Radar Signal Processing," *IEEE J. Sel. Topics Signal Process.*, vol. 7, no. 4, pp. 595-604, 2013.

[8] R. Rosu, M. Donias, L. Bombrun, S. Said, O. Regniers, and J.-P. Da Costa, "Structure Tensor Riemannian Statistical Models for CBIR and Classification of Remote Sensing Images," *IEEE Trans. Geosci. Remote Sens.*, vol. 55, no. 1, pp. 248-260, 2017.

[9] A. Barachant, S. Bonnet, M. Congedo, and C. Jutten, "Classification of Covariance Matrices Using a Riemannian-Based Kernel for BCI Applications," *Neurocomput.*, vol. 112, p. 172-178, July 2013.

[10] D. Brooks, O. Schwander, F. Barbaresco, J.-Y. Schneider, and M. Cord, "Riemannian Batch Normalization for SPD Neural Networks," in *Proc. NeurIPS 2019*, Vancouver, Canada, Dec. 2019, pp. 15489-15500.

[11] J. Solà, J. Deray, and D. Atchuthan, "A Micro Lie Theory for State Estimation in Robotics," 2021.

[12] L. Rodrigues de Almeida, G. A. Giraldi, and M. Bernardes Vieira, "Rigid Registration of Point Clouds Based on Indirect Lie Group Approach," in *Proc. SVR 2019*, Rio de Janeiro, Brazil, Oct. 2019, pp. 130-139.

[13] A. Barrau and S. Bonnabel, "Intrinsic Filtering on Lie groups With Applications to Attitude Estimation," *IEEE Trans. Autom. Control*, vol. 60, no. 2, pp. 436-449, 2015.

[14] T. D. Barfoot and P. T. Furgale, "Associating Uncertainty With Three-Dimensional Poses for Use in Estimation Problems," *IEEE Trans. Robot.*, vol. 30, pp. 679-693, 2014.

[15] S. Bonnabel and A. Barrau, "An intrinsic Cramér-Rao bound on Lie groups," in *Proc. GSI 2015*, Palaiseau, France, Oct. 2015, pp. 664-672.

[16] J. Ćesić, I. Marković, M. Bukal, and I. Petrović, "Extended information filter on matrix Lie groups," *Automatica*, vol. 82, pp. 226-234, 2017.

[17] G. Bourmaud, R. Mégret, A. Giremus, and Y. Berthoumieu, "From Intrinsic Optimization to Iterated Extended Kalman Filtering on Lie Groups," *J. Math. Imag. Vis.*, vol. 55, no. 3, pp. 284-303, Jan. 2016.

[18] S. Labsir, A. Giremus, B. Yver, and T. Benoudiba-Campanini, "Joint Shape and Centroid Position Tracking of a Cluster of Space Debris by Filtering on Lie Groups," *Signal Processing*, vol. 183, p. 108027, 2021.

[19] M. Brossard, S. Bonnabel, and J.-P. Condomines, "Unscented Kalman filtering on Lie groups," in *Proc. IROS 2017*, Vancouver, Canada, Sept. 2017, pp. 2485-2491.

[20] C. Chahbazian, N. Merlinge, K. Dahia, B. Winter-Bonnet, A. Blanc, and C. Musso, "Generalized Laplace Particle Filter on Lie Groups Applied to Ambiguous Doppler Navigation," in *Proc. IROS 2022*, Kyoto, Japan, Oct. 2022, pp. 2387-2394.

[21] P. Galdi, A. Serra, and R. Tagliaferri, "Rotation clustering: A consensus clustering approach to cluster gene expression data," in *Proc. WILF 2016*, Naples, Italy, Dec. 2016, pp. 229-238.

[22] R. Arora, "On Learning Rotations," in *Proc. NIPS 2009*, Vancouver, Canada, Dec. 2009, pp. 55-63.

[23] R. Arora and W. A. Sethares, "An Efficient and Stable Algorithm for Learning Rotations," in *Proc. ICPR 2010*, Istanbul, Turkey, Aug. 2010, pp. 2993-2996.

[24] W. Jiang, Y. Wang, Y. Li, Y. Lin, and W. Shen, "Radar Target Characterization and Deep Learning in Radar Automatic Target Recognition: A Review," *Remote Sensing*, vol. 15, no. 15, 2023.

[25] A. Sibiryakov and M. Bober, "Image Registration using RST-Clustering and its Application in Remote Sensing," in *Proc. SPIE 2006*, Stockholm, Sweden, Sept. 2006.

[26] J. Tang, "Image Registration Using Clustering Algorithm," *Adv. Mater. Res.*, vol. 108-111, pp. 63-68, May 2010.

[27] Y. Cabanes, F. Barbaresco, M. Arnaudon, and J. Bigot, "Toeplitz Hermitian Positive Definite Matrix Machine Learning Based on Fisher Metric," in *Proc. GSI 2019*, Toulouse, France, Aug. 2019, pp. 261-270.

[28] I. Ilea, L. Bombrun, S. Said, and Y. Berthoumieu, "Fisher Vector Coding for Covariance Matrix Descriptors Based on the Log-Euclidean and Affine Invariant Riemannian Metrics," *J. Imag.*, vol. 4, no. 7, p. 85, 2018.

[29] D. A. Bini and B. Iannazzo, "Computing the Karcher mean of symmetric positive definite matrices," *Linear Algebra and its Applications*, vol. 438,

no. 4, pp. 1700-1710, 2013, 16th ILAS Conference Proceedings, Pisa 2010.

[30] J. Hu, R. Ao, A. M.-C. So, M. Yang, and Z. Wen, "Riemannian natural gradient methods," *SIAM Journal on Scientific Computing*, vol. 46, no. 1, pp. A204-A231, 2024.

[31] S. Said, L. Bombrun, Y. Berthoumieu, and J. H. Manton, "Riemannian gaussian distributions on the space of symmetric positive definite matrices," *IEEE Transactions on Information Theory*, vol. 63, no. 4, pp. 2153-2170, 2017.

[32] I. Ilea, L. Bombrun, S. Said, and Y. Berthoumieu, "Covariance Matrices Encoding Based on the Log-Euclidean and Affine Invariant Riemannian Metrics," in *2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (CVPRW)*, 2018, pp. 506-50609.

[33] M. Lu and F. Li, "Survey on Lie group machine learning," *Big Data Mining and Analytics*, vol. 3, no. 4, pp. 235-258, 2020.

[34] M. Zheng, "Clustering of image covariance matrices on lie group manifold," *Autom. Control Comput. Sci.*, vol. 58, no. 2, p. 209-215, April 2024.

[35] R. Subbarao and P. Meer, "Nonlinear mean shift for clustering over analytic manifolds," *Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition*, vol. 1, pp. 1168-1175, 01 2006.

[36] H. Sun, L. Liu, and F. Li, "A Lie group semi-supervised FCM clustering method for image segmentation," *Pattern Recognition*, vol. 155, p. 110681, 2024.

[37] V. Arsigny, P. Fillard, X. Pennec, and N. Ayache, "Geometric Means in a Novel Vector Space Structure on Symmetric Positive-Definite Matrices," *SIAM J. Matrix Anal. Appl.*, vol. 29, no. 1, pp. 328-347, 2007.

[38] V. Arsigny, O. Commowick, X. Pennec, and N. Ayache, "A Log-Euclidean Framework for Statistics on Diffeomorphisms," in *Proc. MICCAI 2006*, Copenhagen, Denmark, Oct. 2006, pp. 924-931.

[39] A. P. Bustos, T.-J. Chin, A. Eriksson, and I. Reid, "Visual SLAM: Why Bundle Adjust?" in *Proc. ICRA 2019*, Montreal, Canada, May 2019, pp. 2385-2391.

[40] D. Gualda, J. Ureña, J. C. García, E. García, and J. Alcalá, "Simultaneous Calibration and Navigation (SCAN) of Multiple Ultrasonic Local Positioning Systems," *Inf. Fusion*, vol. 45, pp. 53-65, 2019.

[41] G. Kovács and S. Nagy, "Ultrasonic Sensor Fusion Inverse Algorithm for Visually Impaired Aiding Applications," *Sensors*, vol. 20, no. 13, p. 3682, 2020.

[42] J. Selig, *Geometrical Methods in Robotics*. Springer, 1996.

[43] B. C. Hall, *Lie Groups, Lie Algebras, and Representations: an Elementary Introduction*. McGraw-Hill Higher Education, 2004.

[44] G. S. Chirikjian, in *Stochastic Models, Information Theory and Lie Groups*, vol. 2. Springer, 2010, ch. 12, p. 84.

[45] P. Absil, R. Mahony, and R. Sepulchre, *Optimization Algorithm on Matrix Manifolds*. Princeton, 2008.

[46] B. Owren and B. Welfert, "The Newton iteration on Lie Groups," *BIT Numerical Math.*, vol. 40, no. 1, pp. 121-145, 2000.

[47] J. Manton and R. Mahony, "The geometry of the Newton method on non-compact Lie group," *J. Global Optimization*, vol. 23, pp. 309-327, 2002.

[48] G. Bourmaud, R. Mégret, A. Giremus, and Y. Berthoumieu, "From Intrinsic Optimization to Iterated Extended Kalman Filtering on Lie Groups," *J. Math. Imag. Vis.*, vol. 55, no. 3, pp. 284-303, Jan. 2016.

[49] X. Pennec, "Intrinsic statistics on Riemannian manifolds: Basic tools for geometric measurements," *J. Math. Imag. Vis.*, vol. 25, pp. 127-155, 2006.

[50] X. Pennec, "Bi-invariant means on Lie groups with Cartan-Schouten connections," in *Proc. GSI 2013*, Paris, France, 2013, pp. 59-67.

[51] G. S. Chirikjian, "Information Theory on Lie Groups and Mobile Robotics Applications," in *Proc. ICRA 2010*, Anchorage, AK, USA, May 2010, pp. 2751-2757.

[52] J. Faraut, "Analysis on Lie Groups: an Introduction." Cambridge, 2008, ch. 3, p. 46.

[53] S. Labsir, A. Renaux, J. Vilà-Valls, and E. Chaumette, "Cramér-Rao Bound on Lie Groups with Observations on Lie Groups: Application to SE(2)," in *Proc. ICASSP 2023*, Rhodes Island, Greece, June 2023, pp. 1-5.

[54] C. M. Bishop, *Pattern Recognition and Machine Learning (Information Science and Statistics)*. Berlin, Heidelberg: Springer-Verlag, 2006.

[55] P. Venkatesan, N. Pavitra, and R. Mohan, "Performing SLAM using Low-Cost Sensors for Autonomous Navigation in Household Environments," in *Proc. 12CT 2019*, Bombay, India, March 2019, pp. 1-5.

[56] U. V. Luxburg, "A Tutorial on Spectral Clustering," *Stat. Comput.*, vol. 17, no. 4, pp. 395-416, 2007.