Riemannian Gaussian Mixtures

• Riemannian Gaussian mixtures are finite mixture models whose densities are defined intrinsically on manifolds, leveraging geodesic distances and affine-invariance.
• They employ an expectation–maximisation algorithm with weighted Riemannian barycentres, enabling efficient and robust clustering and parameter estimation.
• They demonstrate practical advantages in applications such as texture classification and denoising, unifying manifold geometry with statistical modeling.

Searching arXiv for papers on Riemannian Gaussian mixtures and closely related developments. Riemannian Gaussian mixtures are finite mixture models whose component densities are defined intrinsically on a Riemannian manifold rather than in a Euclidean vector space. In the foundational construction on the manifold $\mathcal{P}_m$ of $m\times m$ symmetric positive definite matrices, a Riemannian Gaussian distribution has density

$$p(Y\mid \mu,\sigma)=\frac{1}{\zeta(\sigma)}\,\exp\!\Big(-\frac{d^2(Y,\mu)}{2\,\sigma^2}\Big),$$

with respect to the Riemannian volume element, where $\mu\in\mathcal{P}_m$ is a location parameter, $\sigma>0$ is a dispersion parameter, and $d$ is the geodesic distance induced by the affine-invariant Rao–Fisher metric (Said et al., 2015). A Riemannian Gaussian mixture then takes the form

$$p(Y)=\sum_{k=1}^K \pi_k\,p(Y\mid \mu_k,\sigma_k),\qquad \pi_k>0,\ \sum_{k=1}^K\pi_k=1,$$

and is estimated by an expectation–maximisation procedure whose mean updates are weighted Riemannian barycentres (Said et al., 2015). More generally, related constructions have been developed on compact manifolds, symmetric spaces of structured covariance matrices, Lie groups such as $\mathrm{SE}(3)$, and hyperbolic space, which together position Riemannian Gaussian mixtures as a likelihood-based framework for clustering, density estimation, classification, denoising, and control on non-Euclidean domains (Jaffe et al., 9 Jun 2026).

1. Geometric setting on the manifold of SPD matrices

The classical reference setting is the space

$$\mathcal{P}_m=\{\,Y\in\mathbb{R}^{m\times m}: Y^\dagger=Y,\; x^\dagger Y x>0\ \text{for all}\ x\in\mathbb{R}^m\,\},$$

equipped with the affine-invariant, or Rao–Fisher, Riemannian metric

$$ds^2(Y)=\mathrm{tr}\,\big[\,Y^{-1}dY\,\big]^2, \qquad g_Y(U,V)=\mathrm{tr}\,\big[\,Y^{-1}U\,Y^{-1}V\,\big].$$

This metric is invariant under congruence transformations and inversion: $m\times m$0

$m\times m$1

hence the induced geodesic distance satisfies

$m\times m$2

The unique geodesic from $m\times m$3 to $m\times m$4 is

$m\times m$5

and the squared Rao distance is

$m\times m$6

The Riemannian volume element is

$m\times m$7

with an explicit polar-coordinate form in eigenvalue coordinates (Said et al., 2015).

This geometry is not an auxiliary device; it determines the model class. The density depends on geodesic distance and Riemannian volume rather than on a Euclidean proxy such as $m\times m$8 in a flat space. The data types motivating this choice are explicitly covariance-like or tensor-valued descriptors arising in medical imaging, computer vision, and radar signal processing (Said et al., 2015).

A recurrent misconception is that a “Gaussian on SPD matrices” is simply a Euclidean Gaussian applied to matrix logarithms. The comparison in the source material is explicit: log-Euclidean methods model $m\times m$9 in a Euclidean space, but they lose affine-invariance and the exact manifold geometry, whereas the Rao–Fisher metric is affine-invariant and respects the natural homogeneous structure of $$p(Y\mid \mu,\sigma)=\frac{1}{\zeta(\sigma)}\,\exp\!\Big(-\frac{d^2(Y,\mu)}{2\,\sigma^2}\Big),$$0 (Said et al., 2015).

2. Riemannian Gaussian distributions and the role of the normaliser

On $$p(Y\mid \mu,\sigma)=\frac{1}{\zeta(\sigma)}\,\exp\!\Big(-\frac{d^2(Y,\mu)}{2\,\sigma^2}\Big),$$1, the Riemannian Gaussian distribution $$p(Y\mid \mu,\sigma)=\frac{1}{\zeta(\sigma)}\,\exp\!\Big(-\frac{d^2(Y,\mu)}{2\,\sigma^2}\Big),$$2 is defined by

$$p(Y\mid \mu,\sigma)=\frac{1}{\zeta(\sigma)}\,\exp\!\Big(-\frac{d^2(Y,\mu)}{2\,\sigma^2}\Big),$$3

A key property is that the normalising constant depends only on $$p(Y\mid \mu,\sigma)=\frac{1}{\zeta(\sigma)}\,\exp\!\Big(-\frac{d^2(Y,\mu)}{2\,\sigma^2}\Big),$$4 and $$p(Y\mid \mu,\sigma)=\frac{1}{\zeta(\sigma)}\,\exp\!\Big(-\frac{d^2(Y,\mu)}{2\,\sigma^2}\Big),$$5, not on $$p(Y\mid \mu,\sigma)=\frac{1}{\zeta(\sigma)}\,\exp\!\Big(-\frac{d^2(Y,\mu)}{2\,\sigma^2}\Big),$$6. Writing

$$p(Y\mid \mu,\sigma)=\frac{1}{\zeta(\sigma)}\,\exp\!\Big(-\frac{d^2(Y,\mu)}{2\,\sigma^2}\Big),$$7

the invariance of distance and volume yields

$$p(Y\mid \mu,\sigma)=\frac{1}{\zeta(\sigma)}\,\exp\!\Big(-\frac{d^2(Y,\mu)}{2\,\sigma^2}\Big),$$8

An exact polar-coordinate formula is available: $$p(Y\mid \mu,\sigma)=\frac{1}{\zeta(\sigma)}\,\exp\!\Big(-\frac{d^2(Y,\mu)}{2\,\sigma^2}\Big),$$9 where

$\mu\in\mathcal{P}_m$0

For $\mu\in\mathcal{P}_m$1, the normaliser has the closed form

$\mu\in\mathcal{P}_m$2

The paper introducing this formulation states that it gives an exact expression of the probability density function for the first time in the existing literature (Said et al., 2015).

The exact normaliser matters statistically. Because $\mu\in\mathcal{P}_m$3 is independent of $\mu\in\mathcal{P}_m$4, the maximum-likelihood estimator of the location parameter reduces to minimisation of a Fréchet functional. The source material explicitly identifies this independence as crucial for the MLE–Fréchet mean correspondence (Said et al., 2015). A plausible implication is that the model is not merely geometrically natural; it is also algorithmically structured in a way that cleanly separates location estimation from scale estimation.

The same pattern reappears in broader settings. On compact manifolds, a Riemannian Gaussian likelihood can be written as

$\mu\in\mathcal{P}_m$5

and on homogeneous manifolds the normaliser $\mu\in\mathcal{P}_m$6 does not depend on $\mu\in\mathcal{P}_m$7 (Jaffe et al., 9 Jun 2026). On hyperbolic space under the hyperboloid model, the isotropic density

$\mu\in\mathcal{P}_m$8

has a radial normaliser with an exact finite-sum representation involving the complementary error function (You, 27 Apr 2026). This suggests that the SPD case is part of a wider programme: defining Gaussian-like laws by geodesic distance and deriving the accompanying geometry-specific partition function exactly or numerically.

3. Maximum likelihood, Fréchet means, and the EM algorithm

For i.i.d. samples $\mu\in\mathcal{P}_m$9 on $\sigma>0$0, the log-likelihood is

$\sigma>0$1

Since $\sigma>0$2 does not depend on $\sigma>0$3, the maximum-likelihood estimator $\sigma>0$4 minimises

$\sigma>0$5

so

$\sigma>0$6

This is the empirical Riemannian centre of mass, or Fréchet mean. Existence and uniqueness follow from negative curvature of $\sigma>0$7 under the Rao–Fisher metric (Said et al., 2015).

The intrinsic gradient calculus is explicit. The Riemannian logarithm and exponential maps are

$\sigma>0$8

and

$\sigma>0$9

An unweighted Karcher-type iteration therefore computes

$d$0

with backtracking line search and termination when $d$1 (Said et al., 2015).

The scale parameter is estimated from the unique solution of

$d$2

or equivalently

$d$3

where $d$4 is the inverse of the strictly increasing map $d$5 (Said et al., 2015).

For mixtures, the model is

$d$6

The expectation–maximisation algorithm consists of responsibilities

$d$7

followed by the M-step updates

$d$8

$d$9

$$p(Y)=\sum_{k=1}^K \pi_k\,p(Y\mid \mu_k,\sigma_k),\qquad \pi_k>0,\ \sum_{k=1}^K\pi_k=1,$$0

The mean update is thus a weighted Riemannian barycentre, implemented by the weighted Karcher flow

$$p(Y)=\sum_{k=1}^K \pi_k\,p(Y\mid \mu_k,\sigma_k),\qquad \pi_k>0,\ \sum_{k=1}^K\pi_k=1,$$1

with suitable backtracking (Said et al., 2015).

Standard EM monotonicity applies, and uniqueness of each weighted barycentre follows from negative curvature, so the mean subproblems are well-defined (Said et al., 2015). A plausible implication is that, unlike Euclidean mixture models with singular or flat covariance subproblems, the non-positive curvature of the underlying manifold contributes directly to the well-posedness of the location updates.

4. Statistical properties, computation, and implementation

The population mean parameter $$p(Y)=\sum_{k=1}^K \pi_k\,p(Y\mid \mu_k,\sigma_k),\qquad \pi_k>0,\ \sum_{k=1}^K\pi_k=1,$$2 of $$p(Y)=\sum_{k=1}^K \pi_k\,p(Y\mid \mu_k,\sigma_k),\qquad \pi_k>0,\ \sum_{k=1}^K\pi_k=1,$$3 is the Riemannian centre of mass of the distribution: it minimises

$$p(Y)=\sum_{k=1}^K \pi_k\,p(Y\mid \mu_k,\sigma_k),\qquad \pi_k>0,\ \sum_{k=1}^K\pi_k=1,$$4

The estimator $$p(Y)=\sum_{k=1}^K \pi_k\,p(Y\mid \mu_k,\sigma_k),\qquad \pi_k>0,\ \sum_{k=1}^K\pi_k=1,$$5 is consistent and asymptotically normal in the tangent space at $$p(Y)=\sum_{k=1}^K \pi_k\,p(Y\mid \mu_k,\sigma_k),\qquad \pi_k>0,\ \sum_{k=1}^K\pi_k=1,$$6, with covariance involving

$$p(Y)=\sum_{k=1}^K \pi_k\,p(Y\mid \mu_k,\sigma_k),\qquad \pi_k>0,\ \sum_{k=1}^K\pi_k=1,$$7

where $$p(Y)=\sum_{k=1}^K \pi_k\,p(Y\mid \mu_k,\sigma_k),\qquad \pi_k>0,\ \sum_{k=1}^K\pi_k=1,$$8 (Said et al., 2015).

Sampling from the identity-centred model is also explicit. If $$p(Y)=\sum_{k=1}^K \pi_k\,p(Y\mid \mu_k,\sigma_k),\qquad \pi_k>0,\ \sum_{k=1}^K\pi_k=1,$$9 is uniform on $\mathrm{SE}(3)$0 and $\mathrm{SE}(3)$1 has density proportional to

$\mathrm{SE}(3)$2

then

$\mathrm{SE}(3)$3

and $\mathrm{SE}(3)$4 (Said et al., 2015).

Computationally, the dominant cost is repeated matrix spectral calculus. For each sample–component pair, computing $\mathrm{SE}(3)$5 with

$\mathrm{SE}(3)$6

requires eigen-decomposition at cost $\mathrm{SE}(3)$7, so responsibilities and log maps typically cost $\mathrm{SE}(3)$8 per outer EM iteration, with additional inner iterations for barycentres (Said et al., 2015). Numerical recommendations include eigen-decomposition-based matrix log/exp, thresholding extremely small eigenvalues for stability, computing $\mathrm{SE}(3)$9 from the eigendecomposition $$\mathcal{P}_m=\{\,Y\in\mathbb{R}^{m\times m}: Y^\dagger=Y,\; x^\dagger Y x>0\ \text{for all}\ x\in\mathbb{R}^m\,\},$$0, and exploiting embarrassingly parallel structure across samples and components (Said et al., 2015).

The source material also emphasises offline computation of $$\mathcal{P}_m=\{\,Y\in\mathbb{R}^{m\times m}: Y^\dagger=Y,\; x^\dagger Y x>0\ \text{for all}\ x\in\mathbb{R}^m\,\},$$1 and its derivatives. Once pretabulated, the one-dimensional $$\mathcal{P}_m=\{\,Y\in\mathbb{R}^{m\times m}: Y^\dagger=Y,\; x^\dagger Y x>0\ \text{for all}\ x\in\mathbb{R}^m\,\},$$2 updates are negligible relative to the repeated matrix logarithms and exponentials (Said et al., 2015).

Later work on compact manifolds reveals a related computational tension. Empirical Bayes denoising with Riemannian Gaussian mixture models on compact manifolds relies on the identity

$$\mathcal{P}_m=\{\,Y\in\mathbb{R}^{m\times m}: Y^\dagger=Y,\; x^\dagger Y x>0\ \text{for all}\ x\in\mathbb{R}^m\,\},$$3

where $$\mathcal{P}_m=\{\,Y\in\mathbb{R}^{m\times m}: Y^\dagger=Y,\; x^\dagger Y x>0\ \text{for all}\ x\in\mathbb{R}^m\,\},$$4 is the mixture marginal density (Jaffe et al., 9 Jun 2026). However, the same work shows that cut locus singularities limit smoothness of the density and slow nonparametric rates relative to Euclidean Gaussian mixtures (Jaffe et al., 9 Jun 2026). This suggests that efficient computation and statistical regularity of Riemannian Gaussian mixtures depend not only on dimension, but also on global manifold geometry.

5. Applications and empirical performance

The most detailed empirical study in the source material concerns texture classification on the VisTex database. The setup uses $$\mathcal{P}_m=\{\,Y\in\mathbb{R}^{m\times m}: Y^\dagger=Y,\; x^\dagger Y x>0\ \text{for all}\ x\in\mathbb{R}^m\,\},$$5 images of size $$\mathcal{P}_m=\{\,Y\in\mathbb{R}^{m\times m}: Y^\dagger=Y,\; x^\dagger Y x>0\ \text{for all}\ x\in\mathbb{R}^m\,\},$$6, each partitioned into $$\mathcal{P}_m=\{\,Y\in\mathbb{R}^{m\times m}: Y^\dagger=Y,\; x^\dagger Y x>0\ \text{for all}\ x\in\mathbb{R}^m\,\},$$7 patches of size $$\mathcal{P}_m=\{\,Y\in\mathbb{R}^{m\times m}: Y^\dagger=Y,\; x^\dagger Y x>0\ \text{for all}\ x\in\mathbb{R}^m\,\},$$8 with overlap; $$\mathcal{P}_m=\{\,Y\in\mathbb{R}^{m\times m}: Y^\dagger=Y,\; x^\dagger Y x>0\ \text{for all}\ x\in\mathbb{R}^m\,\},$$9 patches per image are used for training and $$ds^2(Y)=\mathrm{tr}\,\big[\,Y^{-1}dY\,\big]^2, \qquad g_Y(U,V)=\mathrm{tr}\,\big[\,Y^{-1}U\,Y^{-1}V\,\big].$$0 for testing. The texture descriptors are SPD features in $$ds^2(Y)=\mathrm{tr}\,\big[\,Y^{-1}dY\,\big]^2, \qquad g_Y(U,V)=\mathrm{tr}\,\big[\,Y^{-1}U\,Y^{-1}V\,\big].$$1 with $$ds^2(Y)=\mathrm{tr}\,\big[\,Y^{-1}dY\,\big]^2, \qquad g_Y(U,V)=\mathrm{tr}\,\big[\,Y^{-1}U\,Y^{-1}V\,\big].$$2, and class-conditional mixtures are trained with $$ds^2(Y)=\mathrm{tr}\,\big[\,Y^{-1}dY\,\big]^2, \qquad g_Y(U,V)=\mathrm{tr}\,\big[\,Y^{-1}U\,Y^{-1}V\,\big].$$3 or $$ds^2(Y)=\mathrm{tr}\,\big[\,Y^{-1}dY\,\big]^2, \qquad g_Y(U,V)=\mathrm{tr}\,\big[\,Y^{-1}U\,Y^{-1}V\,\big].$$4 components by EM on $$ds^2(Y)=\mathrm{tr}\,\big[\,Y^{-1}dY\,\big]^2, \qquad g_Y(U,V)=\mathrm{tr}\,\big[\,Y^{-1}U\,Y^{-1}V\,\big].$$5 (Said et al., 2015).

The reported overall accuracies, averaged over $$ds^2(Y)=\mathrm{tr}\,\big[\,Y^{-1}dY\,\big]^2, \qquad g_Y(U,V)=\mathrm{tr}\,\big[\,Y^{-1}U\,Y^{-1}V\,\big].$$6 random splits, are as follows.

Method	$$ds^2(Y)=\mathrm{tr}\,\big[\,Y^{-1}dY\,\big]^2, \qquad g_Y(U,V)=\mathrm{tr}\,\big[\,Y^{-1}U\,Y^{-1}V\,\big].$$7	$$ds^2(Y)=\mathrm{tr}\,\big[\,Y^{-1}dY\,\big]^2, \qquad g_Y(U,V)=\mathrm{tr}\,\big[\,Y^{-1}U\,Y^{-1}V\,\big].$$8
Riemannian Gaussian mixture classifier	$$ds^2(Y)=\mathrm{tr}\,\big[\,Y^{-1}dY\,\big]^2, \qquad g_Y(U,V)=\mathrm{tr}\,\big[\,Y^{-1}U\,Y^{-1}V\,\big].$$9	$m\times m$00
Nearest neighbor with Riemannian centers	$m\times m$01	$m\times m$02
Wishart mixture classifier	$m\times m$03	$m\times m$04

The observations stated in the source are that the proposed Riemannian Gaussian mixture classifier outperforms both nearest-neighbor and Wishart-based classifiers, and that using mixtures with $m\times m$05 substantially improves accuracy over single-component models with $m\times m$06 (Said et al., 2015).

Related applications broaden the notion of “Riemannian Gaussian mixture” beyond SPD matrices. On compact manifolds, Riemannian Gaussian mixture models serve as the latent noise model in nonparametric empirical Bayes denoising on spheres and tori, including denoising gamma ray burst locations on $m\times m$07 and Ramachandran plots on $m\times m$08 (Jaffe et al., 9 Jun 2026). On $m\times m$09, mixtures of tangent-space Gaussians defined via Lie $m\times m$10 are used for offline reachability modelling and online reactive control in dynamic grasping, where feasible grasp poses are selected by gradient ascent on the mixture or log-mixture density (Choi et al., 2023). On hyperbolic space, finite mixtures of isotropic Riemannian Gaussians under the hyperboloid model are used as exploratory likelihood-based clustering tools for network embeddings, with exact EM and generalised EM algorithms (You, 27 Apr 2026).

These applications share a common pattern: discrete or Euclidean representations are replaced by a density whose modes, distances, and optimisation steps are intrinsic to the manifold carrying the data. This suggests that the main practical contribution of Riemannian Gaussian mixtures is not only better fit in curved spaces, but also the ability to turn geometry itself into a statistical prior.

6. Generalisations, comparisons, and related developments

The construction on $m\times m$11 generalises to other Riemannian symmetric spaces of non-positive curvature. A detailed extension has been developed for structured covariance matrices, including complex Hermitian positive-definite matrices, Toeplitz HPD matrices, and block-Toeplitz HPD matrices with Toeplitz blocks (Said et al., 2016). In that setting, the general density is

$m\times m$12

with $m\times m$13 independent of $m\times m$14 by transitive isometry-group invariance, and the MLE of the location parameter again coincides with the Riemannian barycentre (Said et al., 2016). The same paper gives mixture-model EM algorithms, density-estimation procedures, and classification rules for these structured covariance spaces (Said et al., 2016).

Hyperbolic space provides a different generalisation. Under the hyperboloid model,

$m\times m$15

the isotropic density

$m\times m$16

admits exact EM and generalised EM algorithms. The location update is again a weighted Fréchet mean, while the inverse-scale update is a one-dimensional strictly convex profile problem (You, 27 Apr 2026). The unrestricted likelihood is singular because a component can collapse onto a single data point with vanishing scale, but a constrained mixture estimator exists on a compact parameter set (You, 27 Apr 2026). This clarifies that likelihood degeneracy, familiar from Euclidean finite mixtures, is not removed by moving to negative curvature.

The denoising work on compact manifolds shifts attention from parameter estimation to posterior inference. There, the Riemannian Gaussian mixture density

$m\times m$17

enters a manifold Tweedie–Eddington identity,

$m\times m$18

and motivates the tangential Bayes denoiser

$m\times m$19

which is nearly Bayes-optimal in a low-noise regime (Jaffe et al., 9 Jun 2026). This suggests that Riemannian Gaussian mixtures are not only endpoint models for clustering, but also generative objects from which one can derive score fields and shrinkage rules.

A further distinction emerges when comparing intrinsic mixture models with optimisation on parameter manifolds. Several papers study Euclidean-data Gaussian mixtures whose covariance matrices are updated by Riemannian optimisation on the SPD cone, using affine-invariant or Bures–Wasserstein metrics (Hosseini et al., 2015, Sembach et al., 2021). Those works are about the geometry of parameter space, whereas the intrinsic Riemannian Gaussian mixture of (Said et al., 2015) models data that themselves lie on a manifold. The two directions intersect technically, but they solve different problems.

A common misconception is therefore to conflate any Gaussian mixture using SPD-valued covariances with an intrinsic Riemannian Gaussian mixture. The source material distinguishes these cases: intrinsic models define densities by geodesic distance and Riemannian volume on the data manifold, whereas parameter-manifold approaches use Riemannian geometry to optimise ordinary Euclidean GMMs (Said et al., 2015, Hosseini et al., 2015).

7. Advantages, limitations, and conceptual significance

Several advantages are explicit in the source material. On $m\times m$20, Riemannian Gaussian mixtures are invariant under congruence, treat the manifold geometry exactly, and have unique barycentres because of negative curvature (Said et al., 2015). Their EM updates have a closed-form structure up to the one-dimensional scale update through the exact $m\times m$21 (Said et al., 2015). On structured covariance spaces, product decompositions of the geometry can make distances, barycentres, and normalisers tractable even for large matrices (Said et al., 2016). On hyperbolic space, strong geodesic convexity of the weighted Fréchet functional guarantees existence and uniqueness of single-component weighted estimators (You, 27 Apr 2026).

The principal limitations are also clear. Computation is dominated by repeated logarithm and exponential evaluations on manifolds such as $m\times m$22 (Said et al., 2015). Exact or numerical evaluation of the normalising constant is required for scale estimation (Said et al., 2015). On compact manifolds, cut locus singularities degrade smoothness and lead to genuinely nonparametric estimation rates that are slower than in Euclidean Gaussian mixture models (Jaffe et al., 9 Jun 2026). On hyperbolic space, unrestricted finite-mixture likelihoods are singular and therefore require scale constraints or related regularisation (You, 27 Apr 2026).

The broader significance of the framework lies in its unification of geometry and statistics. Riemannian Gaussian mixtures formalise the widely used notion of a Riemannian centre of mass as a maximum-likelihood estimator, provide finite-mixture analogues of Gaussian modelling on curved domains, and support downstream procedures including Bayes classification, barycentric clustering, empirical Bayes denoising, reachability analysis, and reactive control (Said et al., 2015, Jaffe et al., 9 Jun 2026, Choi et al., 2023). A plausible implication is that their long-term role is analogous to that of Euclidean Gaussian mixtures: not a single algorithm, but a reusable statistical primitive whose exact implementation is adapted to the geometry of the manifold at hand.