# Parsimonious Gaussian mixture models with piecewise-constant eigenvalue profiles

Tom Szwagier, Pierre-Alexandre Mattei, Charles Bouveyron, Xavier Pennec

► To cite this version:

Tom Szwagier, Pierre-Alexandre Mattei, Charles Bouveyron, Xavier Pennec. Parsimonious Gaussian mixture models with piecewise-constant eigenvalue profiles. Statistics and Computing, In press. hal-05140200v2

HAL Id: hal-05140200

https://inria.hal.science/hal-05140200v2

Submitted on 7 Nov 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.

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

Distributed under a Creative Commons Attribution 4.0 International License
Parsimonious Gaussian mixture models with
piecewise-constant eigenvalue profiles

Tom Szwagier¹*, Pierre-Alexandre Mattei², Charles Bouveyron², Xavier Pennec¹

¹*Université Côte d'Azur, INRIA, Epione Project Team, 2004 Rte des Lucioles, 06902, Sophia Antipolis, France.

²Université Côte d'Azur, INRIA, CNRS, LJAD, Maasai Project Team, 2004 Rte des Lucioles, 06902, Sophia Antipolis, France.

*Corresponding author(s). E-mail(s): tom.szwagier@inria.fr;

Contributing authors: pierre-alexandre.mattei@inria.fr; charles.bouveyron@inria.fr; xavier.pennec@inria.fr;

**Abstract**

Gaussian mixture models (GMMs) are ubiquitous in statistical learning, particularly for unsupervised problems. While full GMMs suffer from the overparameterization of their covariance matrices in high-dimensional spaces, spherical GMMs (with isotropic covariance matrices) certainly lack flexibility to fit certain anisotropic distributions. Connecting these two extremes, we introduce a new family of parsimonious GMMs with piecewise-constant covariance eigenvalue profiles. These extend several low-rank models like the celebrated mixtures of probabilistic principal component analyzers (MPPCA), by enabling any possible sequence of eigenvalue multiplicities. If the latter are prespecified, then we can naturally derive an expectation-maximization (EM) algorithm to learn the mixture parameters. Otherwise, to address the notoriously-challenging issue of jointly learning the mixture parameters and hyperparameters, we propose a componentwise penalized EM algorithm, whose monotonicity is proven. We show the superior likelihood-parsimony tradeoffs achieved by our models on a variety of unsupervised experiments: density fitting, clustering and single-image denoising.

**Keywords:** Gaussian Mixture Models, Parsimonious Models, Eigenvalues, EM Algorithm, Model Selection

# 1 Introduction

Gaussian mixture models (GMMs) (McLachlan and Peel 2000) are ubiquitous in statistical learning. Their probabilistic nature gives them many advantages against classical machine learning algorithms such as random forests (Breiman 2001), support vector machines (Boser et al. 1992), K-Means (Lloyd 1982), HDBSCAN (McInnes et al. 2017) and spectral clustering (Ng et al. 2001). They notably enable uncertainty estimates,

likelihood-based model selection, generative mod-
eling and they can handle missing data. Gaussian-
based models show good performances in low-to-
moderate dimensions and are therefore heavily
used in basic tasks such as classification or clus-
tering. However, the overparameterization of their
covariance matrices—*p*(*p* + 1)/2 parameters in
dimension *p*—is a clear issue in high dimensions,
where data often lack. Hence, there is a need for
simpler models.
Many works have attempted to overcome the high-dimensional issues of GMMs (see for instance the Chapter 8 of Bouveyron et al. (2019) for a review). A widely-used heuristic is to reduce the dimension of the data—e.g. via principal component analysis (Jolliffe 2002)—before fitting the mixture model. This practical method may however obviously discard some important information for the subsequent task; it may for instance mix some clusters that were well separated in the original space (Chang 1983). Other methods—referred to as parsimonious—consider GMMs in the full data space but with fewer parameters. The early works of Banfield and Raftery (1993); Celeux and Govaert (1995); Bensmail and Celeux (1996) consider GMMs with two simple parsimony-inducing principles: parameter sharing (e.g. constraining the mixture components to have the same covariance matrix) and covariance constraints (e.g. imposing diagonal or scalar covariance matrices). However, the former may not reduce enough the number of parameters in high dimensions, while the latter may reduce it too much. As a consequence, many intermediate parsimonious models later appeared, with low-rank (Ghahramani and Hinton 1996; McLachlan et al. 2003; McNicholas and Murphy 2008), sparse (Fop et al. 2019; Cappozzo et al. 2025), Markov random field (Finder et al. 2022), ultrametric (Cavicchia et al. 2024), or structured priors (Heaps and Jermyn 2024) for the covariance matrices, with an important impact (Fraley and Raftery 2002; Scrucca et al. 2016).

Among the parsimonious GMM family, the celebrated mixtures of probabilistic principal component analyzers (MPPCA) (Tipping and Bishop 1999a) assume that each cluster lies around a low-dimensional affine subspace, following a multivariate Gaussian distribution with a specific eigenvalue profile—the smallest eigenvalues being *constant* and equal to the noise level. This model however overlooks an important difficulty which is the one of *dimension selection*; the learning of the mixture parameters indeed requires to prespecify the intrinsic dimensions of the clusters, which are often unknown. The high-dimensional data clustering (HDDC) (Bouveyron et al. 2007a) and high-dimensional discriminant analysis (HDDA) (Bouveyron et al. 2007b) models address this issue by automatically adjusting the clusters dimensions over the iterations of the optimization algorithm,

based on Cattell's scree test (Cattell 1966). The theoretical guarantees of the learning algorithm are however not understood, and more practically, these models still require a hyperparameter as an input—the threshold for Cattell's scree test—which has an important impact on the choice of intrinsic dimensions.

Other methods for the automatic learning of intrinsic dimensions in parsimonious GMMs exist. A natural approach adopted in McNicholas and Murphy (2008) is the grid search combined with statistical model selection criteria such as the Akaike information criterion (AIC) (Akaike 1974), the Bayesian information criterion (BIC) (Schwarz 1978), or the integrated completed likelihood (ICL) (Biernacki et al. 2000). While such a grid search is computationally expensive, few alternatives have been proposed: Ghahramani and Beal (1999) used a sparsity-inducing prior together with variational inference. In contrast, the problem of estimating the intrinsic dimension of a single Gaussian has received a much wider attention, using for instance likelihood-based techniques (Bouveyron et al. 2011), Bayesian model selection (Bishop 1998; Bouveyron et al. 2020), or approximate cross-validation (Josse and Husson 2012).

Meanwhile, the building blocks of low-rank GMMs—the PPCA models of Tipping and Bishop (1999b)—have been extended into Gaussian models with piecewise-constant eigenvalue profiles under the name of *principal subspace analysis* (PSA) (Szwagier and Pennec 2025; Szwagier et al. 2025). The motivations are simple: if some sample eigenvalues are relatively close, then equalizing them slightly reduces the likelihood while *drastically* reducing the number of parameters. In that sense, the PSA models go beyond low-rank models like PPCA in terms of parsimony and flexibility—since any group of eigenvalues can be equalized and not just the smallest ones. The eigenvalues multiplicities are automatically learned by optimizing a penalized log-likelihood criterion—e.g. the BIC. While PSA raises promising perspectives for parsimonious modeling and interpretability, it is certainly not suited for *multi-modal* density modeling—which is the case of most real datasets—or classical learning tasks such as classification and clustering.

In this paper, we address all the limits raised in the previous paragraphs at once with a new
family of parsimonious GMMs termed mixtures of principal subspace analyzers (MPSA). These models consist in the multimodal extension of PSA models, or equivalently the parsimonious generalization of MPPCA models with arbitrary eigenvalues multiplicities. At fixed multiplicities, we propose an expectation-maximization (EM) algorithm (Dempster et al. 1977) for the unsupervised learning of mixture parameters. When the multiplicities are not prespecified, we propose an efficient and provably-monotonous variant of the EM algorithm which automatically adjusts the groups of eigenvalues to equalize within the iterations. The underlying idea is that increasing the penalized log-likelihood componentwise increases the mixture penalized log-likelihood. Importantly, this algorithm is hyperparameter-free, in the sense that the groups of eigenvalues to equalize do not depend on any hyperparameter to tune— compared to low-rank models like MPPCA and HDDC. As we shall discuss it, we believe that our EM variant can apply far beyond the scope of the current paper, providing new principles for integrated parameter estimation and model selection. We finally perform several experiments (density estimation, clustering, single-image denoising) which all show the superior likelihood-parsimony tradeoffs achieved by our models with respect to the classical GMMs.

## 2 Principal subspace analysis

PSA is a Gaussian model with a constraint on the covariance eigenvalues multiplicities. It has been originally introduced to identify and overcome a “curse of isotropy” that arises in PCA when eigenvalues are relatively close. In this section, we summarize these ideas through the lens of parsimonious Gaussian modeling. More details can be found in the original paper (Szwagier and Pennec 2025), including illustrations of the generative models, their pose structure, and the impact on the PCA methodology.

Let

$$C(p) := \left\{ \gamma \in \mathbb{N}_*^d : \sum_{k=1}^d \gamma_k = p, d \in \llbracket 1, p \rrbracket \right\} \quad (1)$$

be the set of compositions of $p \in \mathbb{N}_*$ (Bergeron et al. 1995). Let $\gamma \in C(p)$, $q_k := \sum_{k'=1}^k \gamma_{k'}$, and

let us write $q := q_{d-1}$. The PSA model of type $\gamma$ is defined as the parametric family of Gaussian densities $p(x | \mu, \Sigma) = N(x | \mu, \Sigma)$, where $\mu \in \mathbb{R}^p$ is the mean and $\Sigma \in \text{Sym}_+(p)$ is a covariance matrix with repeated eigenvalues of multiplicities $\gamma$. More precisely, one assumes that the covariance matrix takes the form $\Sigma = \sum_k \lambda_k \Pi_k$, where $\lambda_k \in \mathbb{R}_+$ is the k-th largest eigenvalue (of multiplicity $\gamma_k$) and $\Pi_k \in \text{Sym}(p)$ is the orthogonal projector (of rank $\gamma_k$) onto the associated eigenspace (of dimension $\gamma_k$). Such a set of covariance matrices with eigenvalue multiplicities $\gamma$ is written $\text{Sym}_+(\gamma)$.

By definition, PSA models generalize the celebrated PPCA models of Tipping and Bishop (1999b)—whose $q$ highest eigenvalues are (almost surely) distinct and $p-q$ remaining eigenvalues are equal to the noise level—with types of the form $\gamma = (1, \dots, 1, p-q)$.<fnref n="1" /> PSA models also generalize the isotropic PPCA models of Bouveyron et al. (2011)—whose eigenvalue profile forms two blocks, i.e. with types of the form $\gamma = (q, p-q)$.

PSA models are parameterized by the mean $\mu \in \mathbb{R}^p$ and the covariance $\Sigma \in \text{Sym}_+(\gamma)$. One could therefore think that the number of parameters is $p+p(p+1)/2$. However, the repeated eigenvalues make the actual number of parameters much smaller. Szwagier and Pennec (2025) show that the number of free parameters for a PSA model of type $\gamma$ is

$$\kappa(\gamma) = p + d + \frac{p^2 - \sum_{k=1}^d \gamma_k^2}{2}. \quad (2)$$

Therefore, equalizing $\gamma_k$ eigenvalues decreases the number of parameters quadratically ($O(\gamma_k^2)$).

Given a dataset $(x_i)_{i=1}^n \in \mathbb{R}^p$, the PSA model of type $\gamma$ that is the most likely to have generated the data (i.e. the maximum likelihood estimate, MLE) can be computed in closed-form (Szwagier and Pennec 2025). The optimal mean is the empirical mean $\hat{\mu} = \frac{1}{n} \sum_{i=1}^n x_i$ and the optimal covariance matrix $\hat{\Sigma}$ is a piecewise-constant version of the empirical covariance matrix, where the eigenvalues have been block-averaged according to $\gamma$. More precisely, if the empirical covariance $S := \frac{1}{n} \sum_{i=1}^n (x_i - \bar{x})(x_i - \bar{x})^\top$ eigendecomposes as $S = \sum_{j=1}^p \ell_j v_j v_j^\top$, where $\ell_1 \ge \dots \ge \ell_p > 0$

<footnote n="1">¹We will sometimes use the notation $1^a$, as an abbreviation for $1, \dots, 1$ (where 1 is repeated $a$ times).</footnote>
are the eigenvalues<fnref n="2" /> and $v_1 \perp \cdots \perp v_p$ are the eigenvectors, then

$$ \hat{\Sigma}(\gamma) = \sum_{k=1}^{d} \hat{\lambda}_k(\gamma) \hat{\Pi}_k(\gamma), \quad (3) $$

where

$$ \hat{\lambda}_k(\gamma) = \sum_{j=q_{k-1}+1}^{q_k} \frac{\ell_j}{\gamma_k}, \quad (4) $$

$$ \hat{\Pi}_k(\gamma) = \sum_{j=q_{k-1}+1}^{q_k} v_j v_j^\top \quad (5) $$

are respectively the block-averaged eigenvalues and their associated eigenspaces. The maximum log-likelihood then writes

$$ \ln \mathcal{L}(\hat{\mu}, \hat{\Sigma}(\gamma)) = -\frac{n}{2} \left( p \ln n + \sum_{k=1}^{d} \gamma_k \ln \hat{\lambda}_k(\gamma) + p \right). \quad (6) $$

The BIC is defined as

$$ \text{BIC}(\gamma) := \kappa(\gamma) \ln n - 2 \ln \mathcal{L}(\hat{\mu}, \hat{\Sigma}(\gamma)). \quad (7) $$

This model selection criterion results from an asymptotic approximation of the model evidence and it makes a tradeoff between model complexity and goodness-of-fit. Using the BIC, Szwagier and Pennec (2025) get the important result that two adjacent sample eigenvalues $\ell_j$ and $\ell_{j+1}$ should be equalized when their relative eigengap $\delta(\ell_j, \ell_{j+1}) := \frac{\ell_j - \ell_{j+1}}{\ell_j}$ is below a certain threshold $\delta(n)$ depending on the number of samples $n$ but not on the dimension $p$:

$$ \delta(n) := 2 \left( 1 - n^{2/n} + n^{1/n} \sqrt{n^{2/n} - 1} \right). \quad (8) $$

Through extensive empirical experiments, Szwagier and Pennec (2025) show that one rarely has enough samples to distinguish all the empirical eigenvalues, and that one should therefore transition from principal components to principal subspaces.

This important relative eigengap result suggests two simple methodologies to choose the groups of eigenvalues to equalize. The relative eigengap strategy consists in block-averaging the groups of eigenvalues whose relative eigengap is below the threshold $\delta(n)$. The hierarchical clustering strategy consists in a hierarchical clustering (Ward 1963) of the eigenvalues, yielding a family of models $M_h(\ell) \in \mathcal{C}(p)^p$ of decreasing complexity, from which one chooses the model with the lowest BIC. The detailed algorithm is given in Szwagier and Pennec (2025). While the relative eigengap strategy is fast ($O(p)$) compared to the hierarchical clustering strategy ($O(p^2)$), one can show that the former always selects a PSA model with greater or equal BIC, since the latter optimizes over a larger set.<fnref n="3" /> The more recent eigengap sparsity of Szwagier et al. (2025) relaxes the discrete BIC minimization problem into a continuous problem penalized by pairwise distances between sample eigenvalues (Tyler and Yi 2020).

While PSA has an important impact on the interpretability of principal components, its linear-Gaussian nature is too simplistic for real datasets. Moreover, it is not meant to be used for classical learning tasks such as classification and clustering. Since PSA is nothing but a parsimonious Gaussian model, it could naturally be extended into a parsimonious Gaussian mixture model. This is the purpose of next section. As we shall see, such a seemingly natural extension is actually not-so-straightforward from a technical point of view since practical difficulties emerge from the high-dimensional nature of the problem.

## 3 Mixtures of principal subspace analyzers

In this section, we introduce the MPSA models and propose an EM algorithm for learning the mixture parameters in the unsupervised context, with typical applications in density modeling and clustering. The proofs can be found in Appendix A.

<footnote n="2">²In the whole paper, for simplicity, we assume that the sample covariance matrices are always invertible. We make this assumption satisfied in practice by adding a covariance regularization term to the diagonal of the covariance matrices.</footnote>

<footnote n="3">³Indeed the (single-linkage) hierarchical clustering boils down to regrouping the eigenvalues sequentially according to their relative eigengap. Therefore, at some point, all the eigenvalue pairs below the relative eigengap threshold (8) are clustered, similarly to the relative eigengap strategy. In other words, the family $M_h(\ell) \in \mathcal{C}(p)^p$ contains the model output by the relative eigengap strategy.</footnote>
## 3.1 Mixture density and parameters

**Definition 1** (MPSA density) Let $C \in \mathbb{N}^*$ be a number of clusters. Let $\gamma := (\gamma_1, \dots, \gamma_C) \in \mathcal{C}(p)^C$ be a sequence of PSA types—i.e. for $c \in [1, C]$, $\gamma_c := (\gamma_{c1}, \dots, \gamma_{cd_c}) \in \mathbb{N}_*^{d_c}$ with $\sum_{k=1}^{d_c} \gamma_{ck} = p$ and $d_c \in [1, p]$. The MPSA model of types $\gamma$ is defined as the parametric family of probability densities

$$p(x | \theta) := \sum_{c=1}^{C} \pi_c \mathcal{N}(x | \mu_c, \Sigma_c), \quad (9)$$

where $\pi_c \in \mathbb{R}_+$, $\mu_c \in \mathbb{R}^p$ and $\Sigma_c \in \text{Sym}_+(\gamma_c)$. We denote the parameters $\theta := (\pi_c, \mu_c, \Sigma_c)_{c=1}^C \in \Theta(\gamma)$.

**Remark 1** (Generalization of parsimonious GMMs) As PSA extends and unifies PPCA with other parsimonious Gaussian models (like IPPCA, spherical and full-Gaussian models), MPSA extends and unifies MPPCA with other classical parsimonious GMMs belonging to Bouveyron et al. (2007a) and McNicholas and Murphy (2008).

Since no parameter is shared between the mixture components, the total number of parameters can be directly obtained from Equation (2).

**Proposition 1** (MPSA complexity) The number of parameters for the MPSA model of types $\gamma$ is

$$\kappa(\gamma) = C - 1 + \sum_{c=1}^{C} \kappa(\gamma_c). \quad (10)$$

## 3.2 Inference and EM algorithm

Let us consider the problem of mixture parameter estimation under the MPSA model of types $\gamma$. In the unsupervised setting, the class labels are unknown. The classical EM algorithm for mixture parameters estimation then consists in alternating between the computation of the expected labels given the mixture parameters (E-step) and the maximum likelihood estimation of the mixture parameters given the expected labels (M-step). The computational interest lies in the explicit solutions to the M-step, which can be performed componentwise. Moreover, the EM algorithm is guaranteed to increase the mixture log-likelihood $\ln \mathcal{L}(\theta) := \sum_{i=1}^n \ln p(x_i | \theta)$ at each step, notably by virtue of Gibbs' inequality. We summarize our EM algorithm for MPSA inference in Algorithm 1 and give the details of each step in the remainder of this subsection. The posterior probabilities

of class membership ($t_{ci} := \mathbb{P}(y_i = c | x_i)$), where $y_i \in [1, C]$ is a cluster indicator variable), which are explicit later, are gathered in a matrix $t \in \mathbb{R}^{C \times n}$ for conciseness in the pseudo-code.

**Algorithm 1** EM algorithm for MPSA

**Require:** $(x_i)_{i=1}^n \in \mathbb{R}^{p \times n}$, $C \in \mathbb{N}_*$, $\gamma \in \mathcal{C}(p)^C$

**Ensure:** $\hat{\theta} \in \Theta(\gamma) \triangleright$ mixture parameters
$y \leftarrow \text{KMeans}(x, C) \in [1, C]^n \triangleright$ K-Means init.
$t \leftarrow (\mathbf{1}_{y_i=c})_{c,i} \triangleright$ posterior proba.
**for** $s = 1, 2, \dots$ **do**
  $\hat{\theta} \leftarrow \text{M-step}(t, \gamma)$ \quad $\triangleright$ thm. 3
  $t \leftarrow \text{E-step}(\hat{\theta})$ \quad $\triangleright$ thm. 2
**end for**

**Remark 2** (Supervised setting) One can directly use Algorithm 1 for classification (supervised setting) by initializing the cluster indexes with the true labels rather than with K-Means, and running only one iteration.

### 3.2.1 E-step

Let us assume that a dataset $(x_i)_{i=1}^n$ has been sampled from an MPSA model of types $\gamma$ with known parameters $\theta$—i.e. according to the mixture density of Equation (9). Then the probability that a data point $x \in \mathbb{R}^p$ belongs to cluster $c \in [1, C]$ can be computed using Bayes' theorem:

$$\mathbb{P}(y = c | x, \theta) = \frac{\pi_c \mathcal{N}(x | \mu_c, \Sigma_c)}{\sum_{c'=1}^{C} \pi_{c'} \mathcal{N}(x | \mu_{c'}, \Sigma_{c'})} \quad (11)$$

We simplify this formula in the next theorem, referred to as E-step, based on the eigendecomposition of the covariance matrices.

**Theorem 2** (E-step) Consider an MPSA model of types $\gamma$ with parameters $\theta$. Let the covariance matrices decompose as $\Sigma_c = \sum_{k=1}^{d_c} \lambda_{ck} \Pi_{ck}$, where $\lambda_{ck}$ are the (ordered decreasing) eigenvalues and $\Pi_{ck}$ the associated eigenspaces, of dimension $\gamma_{ck}$. Then the probability that $x \in \mathbb{R}^p$ belongs to cluster $c \in [1, C]$ is

$$\mathbb{P}(y = c | x, \theta) = \frac{1}{\sum_{c'=1}^{C} \exp\left(\frac{K_c(x) - K_{c'}(x)}{2}\right)}, \quad (12)$$

where $K_c$ corresponds to the cost function
$$
K_c(x) = -2 \ln \pi_c + \sum_{k=1}^{d_c} \gamma_{ck} \ln \lambda_{ck} \\
+ \sum_{k=1}^{d_c} \frac{(x - \mu_c)^\top \Pi_{ck} (x - \mu_c)}{\lambda_{ck}}. \quad (13)
$$

The log-likelihood of a dataset $(x_i)_{i=1}^n$ then writes

$$
\ln \mathcal{L}(\theta) = \sum_{i=1}^{n} \ln \left( \sum_{c=1}^{C} \exp \left( -\frac{K_c(x_i) + p \ln 2\pi}{2} \right) \right). \quad (14)
$$

3.2.2 M-step

Let us now assume that a dataset $(x_i)_{i=1}^n$ has been
sampled from an MPSA model of types $\gamma$ with
unknown parameters. In the supervised context,
one has access to the class-labels $y_i \in [1, C]$, so
that the mixture parameters can be estimated by
maximizing the complete-data log-likelihood

$$
\ln \mathcal{L}_c(\hat{\theta}) := \sum_{c,i} \mathbf{1}_{y_i=c} \left( \ln \hat{\pi}_c + \ln \mathcal{N}(x_i | \hat{\mu}_c, \hat{\Sigma}_c) \right). \quad (15)
$$

In the unsupervised context, one does not
have access to the labels, but the previous
Theorem 2 (E-step) gives a probabilistic estimate
$t_{ci} := \mathbb{P}(y_i = c | x_i, \theta)$. Therefore, as classically
done, we estimate the mixture parameters by max-
imizing the *expected complete-data log-likelihood*

$$
\mathbb{E}[\ln \mathcal{L}_c(\hat{\theta})] = \sum_{c,i} t_{ci} \left( \ln \hat{\pi}_c + \ln \mathcal{N}(x_i | \hat{\mu}_c, \hat{\Sigma}_c) \right). \quad (16)
$$

**Theorem 3 (M-step)** Let $(x_i)_{i=1}^n$ be a dataset in $\mathbb{R}^p$, and $t_{ci} := \mathbb{P}(y_i = c | x_i)$ the expected class labels. The MPSA parameters that maximize the expected complete-data log-likelihood (16) are

$$
\hat{\pi}_c(t_c) = \frac{1}{n} \sum_{i=1}^{n} t_{ci}, \qquad (17)
$$

$$
\hat{\mu}_c(t_c) = \frac{1}{n\hat{\pi}_c} \sum_{i=1}^{n} t_{ci} x_i, \qquad (18)
$$

$$
\hat{\Sigma}_c(t_c, \gamma_c) = \sum_{k=1}^{d_c} \hat{\lambda}_{ck}(t_c, \gamma_c) \hat{\Pi}_{ck}(t_c, \gamma_c), \quad (19)
$$

where $\hat{\lambda}_{ck}(t_c, \gamma_c)$ and $\hat{\Pi}_{ck}(t_c, \gamma_c)$ are the block-averaged eigenvalues and eigenspaces of the covariance matrix

$$
S_c(t_c) := \frac{1}{n\hat{\pi}_c} \sum_{i=1}^{n} t_{ci} (x_i - \hat{\mu}_c) (x_i - \hat{\mu}_c)^{\top}. \quad (20)
$$

More precisely, if $\ell_{c1} \ge \cdots \ge \ell_{cp}$ are the eigenvalues of $S_c$ and $v_{c1} \perp \cdots \perp v_{cp}$ are the associated eigenvectors, then one has

$$
\hat{\lambda}_{ck}(t_c, \gamma_c) = \sum_{j=q_{c(k-1)}+1}^{q_{ck}} \frac{\ell_{cj}(t_c)}{\gamma_{ck}}, \quad (21)
$$

$$
\hat{\Pi}_{ck}(t_c, \gamma_c) = \sum_{j=q_{c(k-1)}+1}^{q_{ck}} v_{cj}(t_c) v_{cj}(t_c)^{\top}. \quad (22)
$$

Hence, the M-step consists in a classical Gaus-
sian mixture M-step (i.e. computing the local
responsibility-weighted covariance matrices S_c)
followed by a PSA inference step (i.e. block-
averaging the eigenvalues according the pre-
specified eigenvalues multiplicities γ_c).

3.3 EM complexity analysis

We now analyze the time complexity of the pro-
posed EM algorithm 1. For simplicity, we assume
that $n > p$. The case $n \le p$ can be tackled via
transpose tricks combined with covariance regu-
larization to avoid singular covariance matrices.

The K-Means clustering algorithm (Lloyd
1982) involved in the EM initialization has time
complexity O(CnpT), with T the number of
iterations, which is by default at most 300 in
the scikit-learn implementation (Pedregosa et al.
2011).

A classical E-step would have time complexity
O(Cp³) due to the covariance inversion. How-
ever, the formulas that we derive in Theorem 2
have several computational advantages. First,
they remove the necessity to numerically invert
the covariance matrices—which can be patho-
logical in high dimensions. Since the eigenvalues
λck ∈ ℝ and the eigenvectors Qck ∈ ℝp×γck
are computed in the M-step, the computa-
tional complexity of the cost Kc(x) is the one
of ∑k λck−¹(x − μc)ᵀΠck(x − μc), which can be
rewritten as ∑k λck−¹||Qck Ḥ (x − μc)||², that is
∑d_c=1 O(γckp) = O(p²). Therefore, the complexity
of the E-step is O(Cnp²). But this is not all. Since
(Πc1, ..., Πcdc) represents a flag of subspaces (Ye
et al. 2022)—i.e. mutually-orthogonal subspaces
spanning the whole space ℝp—then any orthog-
onal projection matrix Πck can be expressed as
Πck = I − ∑k'≠k Πck'. Notably, if one projec-
tor is of much greater rank than the other ones
(for instance the last one Πcdc, accounting for
the noise), then it can be replaced with lower dimensional terms in (13) and yield much more computationally-efficient formulas for the E-step. More precisely, if one has $q_c := p - \gamma_{cd_c} = O(1)$, then the computational complexity of the cost function $K_c$ drops to $O(q_c p) = O(p)$. Eventually, the complexity of the E-step in high dimensions will more likely be $O(Cnpq)$, with $q := \max_{c \in [1, C]} q_c$.

The M-step involves the eigenvalue decomposition (EVD) of the local responsibility-weighted covariance matrices $S_c$. As currently implemented for simplicity, computing both the covariance matrices and their EVDs results in a complexity of $O(Cnp^2) + O(Cp^3)$. Therefore, even though the eigenvalue equalization of Theorem 3 is extremely simple, it does not remove the need for the costly EVD. To make it more efficient in high dimensions, one can use classical numerical tricks such as the truncated singular value decomposition (SVD), which computes only a few leading eigenvectors. This perspective is investigated in Section B.4.

In Algorithm 1, the mixture model types $\gamma$ are assumed to be prespecified. Unless they are known, a classical way to choose these hyper-parameters is via cross-validation or statistical model selection techniques, implying running several EM algorithms (McNicholas and Murphy 2008). The important issue here is that the number of candidate MPSA models is $(2^{p-1})^C$. Indeed, each independent mixture component is parameterized by its type $\gamma_c$ which can take $2^{p-1}$ values, since each adjacent-eigenvalue pair can be either equalized or not. Therefore, an exhaustive model selection would imply running $O(2^C p)$ EM algorithms, which becomes computationally prohibitive for $p \gtrsim 10$—even though the EM itself may not be too much computationally costly after integrating the previously-exposed low-dimensional numerical tricks.

Drawing from the HDDC models of Bouveyron et al. (2007a), one could naturally be tempted to perform some kind of model selection at the component level (i.e. based on the PSA methodology of Szwagier and Pennec (2025) presented in Section 2) at each EM iteration. This would add a negligible time complexity to Algorithm 1 ($O(Cp^2)$ for the hierarchical clustering strategy), instead of a multiplicative factor. However, such a

variant of the EM algorithm would not enjoy any theoretical guarantee a priori. Consequently, we propose in the next section to reframe our intuitions as a penalized optimization problem, for which we show that a thoughtful variant of the EM algorithm yields a monotonous loss function.

# 4 A penalized EM algorithm

We propose to reframe the joint learning of the mixture parameters and hyperparameters (eigen-value multiplicities) as a penalized log-likelihood optimization problem:

$$ \arg\max_{\gamma \in C(p)^C \atop (\pi_c, \mu_c, \Sigma_c)_{c=1}^C \in \Theta(\gamma)} \sum_{i=1}^n \ln \left( \sum_{c=1}^C \pi_c N(x_i | \mu_c, \Sigma_c) \right) - \alpha \kappa(\gamma), \quad (23) $$

where $\alpha \in \mathbb{R}_+$ is a regularization hyperparameter (for instance $\alpha = \ln(n)/2$ for a BIC-like penalized log-likelihood).

## 4.1 Componentwise monotonicity implies mixture monotonicity

It is well known from the classical EM theory (Dempster et al. 1977) that increasing the expected complete-data log-likelihood $E_{y|x,\theta}[\ln \mathcal{L}_c(\cdot)]$ (16) increases the mixture log-likelihood $\ln \mathcal{L}(\cdot)$ (14). This is notably due to Gibbs’ inequality—i.e. the non-negativity of the relative entropy. In this work, we are not interested in increasing the log-likelihood but a penalized version of it (23), to enforce parsimonious modeling. As said earlier, we would like to select the eigenvalues multiplicities with componentwise model selection strategies drawn from the PSA methodology. However, such a variant of the EM algorithm would not enjoy a priori any guarantee on the evolution of the objective function (23) over the iterations. The following theorem states that a thoroughly-designed variant of the EM algorithm where we replace the M-step with a component-wise model selection guarantees the monotonicity of the penalized log-likelihood. We call this variant componentwise penalized EM (CPEM).

**Theorem 4 (CPEM)** Let $\gamma \in C(p)^C$, $\theta \in \Theta(\gamma)$, $t_{ci} := \mathbb{P}(y_i = c | x_i, \theta)$ (cf. Thm. 2) and $\Psi_c(\gamma_c) :=$
$\sum_{i=1}^{n} t_{ci} \ln \mathcal{N}(x_i | \hat{\mu}_c(t_c), \hat{\Sigma}_c(t_c, \gamma_c))$ (cf. Thm. 3). Let $\gamma' \in \mathcal{C}(p)^C$ such that $\forall c \in [1, C]$,

$$ \Psi_c(\gamma_c') - \alpha \kappa(\gamma_c') \geq \Psi_c(\gamma_c) - \alpha \kappa(\gamma_c), \quad (24) $$

and let $\theta(\gamma') := (\hat{\pi}_c(t_c), \hat{\mu}_c(t_c), \hat{\Sigma}_c(t_c, \gamma_c'))_{c \in [1, C]}$. Then one has:

$$ \ln \mathcal{L}(\theta(\gamma')) - \alpha \kappa(\gamma') \geq \ln \mathcal{L}(\theta) - \alpha \kappa(\gamma). \quad (25) $$

In other words, thoroughly updating the eigenvalue multiplicities componentwise via (24) guarantees the monotonicity of the objective function (23) over the iterations of the EM algorithm.

**Remark 3** (Simplification of componentwise inequalities) Using Theorem 2, one can show with a similar reasoning to Equation (6) that $\Psi_c(\gamma_c) = -\frac{n\pi_c}{2} \sum_{k=1}^{d_c} \gamma_{ck} \ln \hat{\lambda}_{ck}(t_c, \gamma_c) + \text{cst}$, which simplifies the implementation. For instance, with a BIC-like penalization ($\alpha = \ln(n)/2$), the problem of model selection (24) can simply be formulated as

$$ \arg\min_{\gamma_c \in \mathcal{C}(p)} \sum_{k=1}^{d_c} \gamma_{ck} \ln \hat{\lambda}_{ck}(t_c, \gamma_c) + \kappa(\gamma_c) \frac{\ln n}{n\pi_c}. \quad (26) $$

**Remark 4** (Impact on low-rank mixture models) Theorem 4 applies straightforwardly to other parsimonious GMMs such as the mixtures of PPCA of Tipping and Bishop (1999a), the HDDC models of Bouveyron et al. (2007a) and the parsimonious models of McNicholas and Murphy (2008) (which the MPSA models generalize), by choosing types of the form $\gamma_c = (1, \dots, 1, p - q_c)$. Therefore, our results implicitly address the important issue of automatic intrinsic dimension selection in low-rank mixture models that we raised in Section 1. At each step of the EM, for each mixture component $c$, we could select the intrinsic dimension $q_c$ that solves

$$ \arg\min_{q \in \{0, \dots, p-1\}} \sum_{j=1}^{q} \ln \ell_{cj} + (p-q) \ln \left( \frac{1}{p-q} \sum_{j=q+1}^{p} \ell_{cj} \right) + \left( pq - \frac{q(q-1)}{2} \right) \frac{\ln n}{n\pi_c}. \quad (27) $$

The point is that our MPSA models contain the MPPCA models, so that solving (26)—which is what we do in this paper—always yields more adapted models than what we could get from solving (27).

**Remark 5** (Novelty of CPEM) As far as we know, our componentwise penalized EM monotonicity theorem is new. We believe that it can have an important impact in the parsimonious mixture models community (Tipping and Bishop 1999a; McNicholas and Murphy

2008; Scrucca et al. 2016). Indeed, the problem of hyperparameter selection in mixture models is notoriously challenging and open, and it is often addressed heuristically without theoretical guarantees (Bouveyron et al. 2007a). Therefore, the monotonicity result that we obtain addresses this important issue. Let us note that the idea of integrating both parameter estimation and model selection in an EM algorithm is absolutely not new; it was notably proposed in the seminal work of Figueiredo and Jain (2002) for finite mixture models and theoretically investigated in Jiang et al. (2015) for missing data. Let us also note that our algorithm has little to do with the componentwise EM of Celeux et al. (2001), which updates the mixture parameters sequentially, one component per iteration.

## 4.2 Algorithm

Theorem 4 motivates us to design a variant of the EM algorithm where the componentwise quantities $\Psi_c(\cdot) - \alpha \kappa(\cdot)$ increase at each EM iteration. For that, at each step $s$, for each mixture component $c$, we consider a set of candidate models $\mathcal{M}(\gamma_c) \subset \mathcal{C}(p)$ and take the one that maximizes $\Psi_c(\cdot) - \alpha \kappa(\cdot)$. The general methodology is summarized in Algorithm 2.

**Algorithm 2** CPEM algorithm for MPSA

**Require:** $(x_i)_{i=1}^n \in \mathbb{R}^{p \times n}, C \in \mathbb{N}_*, \gamma_0 \in \mathcal{C}(p)^C$

**Ensure:** $\hat{\gamma} \in \mathcal{C}(p)^C, \hat{\theta} \in \Theta(\hat{\gamma})$

$y \leftarrow \text{KMeans}(x, C) \in [1, C]^n \triangleright \text{K-Means init.}$

$t \leftarrow (\mathbf{1}_{y_i=c})_{c,i} \triangleright \text{posterior proba.}$

$\hat{\gamma} \leftarrow \gamma_0 \triangleright \text{types init. (spherical or full)}$

**for** $s = 1, 2, \dots$ **do**

$\hat{\gamma} \leftarrow \arg\max_{\gamma_c \in \{\hat{\gamma}_c\} \cup \mathcal{M}(\hat{\gamma}_c)} \Psi_c(\gamma_c) - \alpha \kappa(\gamma_c) \triangleright \text{thm. 4}$

$\hat{\theta} \leftarrow \text{M-step}(t, \hat{\gamma}) \triangleright \text{thm. 3}$

$t \leftarrow \text{E-step}(\hat{\theta}) \triangleright \text{thm. 2}$

**end for**

**Remark 6** (Satisfying monotonicity conditions) In order to ensure that the theoretical monotonicity conditions from Theorem 4 are met at each iteration—regardless of the set of candidate models $\mathcal{M}(\hat{\gamma}_c)$—we include the current model $\hat{\gamma}_c$ in the set of candidate models. Indeed, not switching model at a CPEM iteration corresponds to a classical M-step from the EM algorithm and it naturally satisfies Equation (24) (with $\gamma'_c = \gamma_c$). Consequently, the selected model,
arg max$_{γ_c \in \{\hat{γ}_c\} \cup \mathcal{M}(\hat{γ}_c)} \Psi_c(γ_c) - α \kappa(γ_c)$ satisfies the theoretical monotonicity conditions too, a fortiori.

There are many possibilities for the choice of candidate models $\mathcal{M}(\gamma_c)$—the idea being of course to avoid considering all the possible models, i.e. $\mathcal{M}(\gamma_c) = \mathcal{C}(p)$, in which case we would have $2^{p-1}$ models to compare at each step $C$. Two ideas have been mentioned in Section 2: the relative eigen-gap strategy—for which $\mathcal{M}(\gamma_c)$ consists in the singleton containing the model obtained by clustering the eigenvalues of $S_c$ whose relative distance is below the threshold (8)—and the hierarchical clustering strategy—for which $\mathcal{M}(\gamma_c)$ consists in the union of models obtained by the hierarchical clustering of the eigenvalues $M_h(l_c) \in \mathcal{C}(p)^p$, the metric being the relative distance.

We now introduce two additional ideas, borrowing from classical paradigms of model selection and sometimes referred to as top-down and bottom-up approaches (see Bach et al. (2012) for uses in sparse optimization problems and Verbockhaven et al. (2024) for uses in neural architecture search). Let $\gamma := (\gamma_1, \dots, \gamma_d) \in \mathcal{C}(p)$. We say that $\gamma' \in \mathcal{C}(p)$ is an upper neighbor of $\gamma$ and write $\gamma' \in \mathcal{N}^+(\gamma)$ if it can be expressed as $\gamma' = (\gamma_1, \dots, \gamma_{k-1}, \gamma_{k1}, \gamma_{k2}, \gamma_{k+1}, \dots, \gamma_d)$ with $(\gamma_{k1}, \gamma_{k2}) \in \mathcal{C}(\gamma_k)$, for any $k \in \{1, \dots, d\}$. Similarly, we say that $\gamma' \in \mathcal{C}(p)$ is a lower neighbor of $\gamma$ and write $\gamma' \in \mathcal{N}^-(\gamma)$ if it can be expressed as $\gamma' = (\gamma_1, \dots, \gamma_{k-1}, \gamma_k + \gamma_{k+1}, \gamma_{k+2}, \dots, \gamma_d)$. One can easily show that $\#\mathcal{N}^+(\gamma) = p - d$ and $\#\mathcal{N}^-(\gamma) = d - 1$.<fnref n="4" /> The bottom-up strategy consists in starting from spherical GMMs $(\gamma_c = (p)^{\forall c})$ and increasing step by step the model complexity by finding the “split” in the eigenvalue profile that maximizes the penalized log-likelihood (i.e. taking $\mathcal{M}(\gamma_c) = \mathcal{N}^+(\gamma_c)$). The top-down strategy is similar but instead starts from the full GMMs $(\gamma_c = (1, \dots, 1)^{\forall c})$ and decreases the complexity at each step by equalizing the two adjacent groups of eigenvalues that maximize the penalized log-likelihood (i.e. taking $\mathcal{M}(\gamma_c) = \mathcal{N}^-(\gamma_c)$).

*Remark 7 (Up and down)* Early experiments showed that the above-defined bottom-up and top-down

strategies may get stuck in local optima due to their monotonous behavior in terms of model complexity. For instance, the bottom-up approach may increase too much the number of parameters during the first few iterations (similarly to Houdard et al. (2018), Figure 7) because one needs a lot of degrees of freedom to model the density when the mixture centers are initialized far from the true cluster centers. Once the true centers are found, one may realize that the clusters have a much simpler structure and one would like to reduce the number of parameters (as we can see on the bottom of Figure 2 for the red curve). To enable such late changes, we eventually take $\mathcal{M}(\gamma_c) = \mathcal{N}^+(\gamma_c) \cup \mathcal{N}^-(\gamma_c)$ for both bottom-up and top-down strategies.

*Remark 8 (Which heuristic to choose?)* Four strategies for updating the eigenvalue multiplicities during the CPEM algorithm have been proposed in this subsection, all of which satisfy the theoretical monotonicity conditions (cf. Remark 6). While theory shows that each strategy yields a non-decreasing sequence of penalized log-likelihoods, it does not show how much the final parameters increase the penalized log-likelihood—with respect to initialization, to more classical GMMs, or in between the different strategies. Gibbs’ inequality only implies that the penalized log-likelihood increases at least as much as its expected complete-data version. Therefore, one can expect complex strategies to outperform simple ones (for instance, the hierarchical clustering strategy should outperform the relative eigengap strategy, that it encompasses, as discussed in the end of Section 2) but one cannot really prove it. The empirical experiments of the next section aim to better understand the respective advantages of each strategy.

# 5 Experiments

In this section, we perform various density modeling and clustering experiments with our new MPSA models (MPSA-H for the hierarchical clustering strategy, MPSA-R for the relative eigengap strategy, MPSA-U for the bottom-up strategy and MPSA-D for the top-down strategy), and compare them to the more standard full Gaussian mixture model (GMM-F) and the spherical one (GMM-S). The code is publicly available on GitHub: https://github.com/tomszwagier/MPSA. Technical implementation details can be found in the appendix (sections B.4 and B.6).

We notably perform some experiments on a synthetic dataset following the MPSA density

<footnote n="4">⁴These definitions are connected to the notion of type refinement introduced in Szwagier and Pennec (2025), and they can be more-intuitively visualized via the Hasse diagram representation of PSA models from the paper.</footnote>
$p(x | \theta) = \sum_{c=1}^{C} \pi_c N(x | \mu_c, \Sigma_c)$ (9). While $C$ and $\pi_c$ are arbitrary and specified for each experiment, the choice of other hyperparameters follow certain rules that we describe in this paragraph. The mean vectors are sampled uniformly between some prespecified bounds, i.e. $\mu_c \sim U([-b, b]^p)$. The covariance matrices $\Sigma_c = Q_c \Lambda_c Q_c^\top$ are such that $Q_c$ is sampled uniformly (according to Haar's measure) in $O(p)$, and $\Lambda_c$ has a constant relative eigengap $\delta$ between groups of eigenvalues, i.e. $(\lambda_{cj} - \lambda_{c(j+1)}) / \lambda_{cj} = \delta$ if $j = q_{ck}$, else $\lambda_j = \lambda_{j+1}$. The latter implying that $\lambda_{cp}/\lambda_{c1} = (1 - \delta)^{d_c-1}$, one then only needs to specify the highest eigenvalue $\lambda_{c1}$, the signal-to-noise ratio $\lambda_{cp}/\lambda_{c1}$ and the type $\gamma_c$ to fully determine the law of the covariance matrix $S_c$. We can now dive into the experiments.

**Remark 9 (Choice of $\alpha$)** The penalty factor is set to $\alpha = \ln(n)/2$ in all the experiments, which corresponds to a BIC-like penalized log-likelihood. The choice of $\alpha$ certainly has an influence on the results— we can naturally expect the number of parameters of the selected model to decrease with $\alpha$. We conduct a sensitivity analysis in Section B.1, showing notably that one recovers the full GMM in the small $\alpha$ limit and the spherical GMM in the large $\alpha$ limit, and that the BIC tends to select (slightly) underparameterized models while the AIC tends to select (slightly) overparameterized models.

**Remark 10 (Choice of $C$)** The number of clusters is set to the true one (if known) in the following experiments. Indeed, this paper focuses on the parsimony induced by equalizing close covariance eigenvalues rather than reducing the number of mixture components. However, for completeness, we discuss in Section B.2 how $C$ can be inferred asymptotically via a model selection framework or alternative methodologies relying on Dirichlet-type priors on the mixture proportions.

## 5.1 CPEM in action for 2D data

This experiment aims at understanding how the mixture parameters as well as the types evolve during the learning (Algorithm 2) with a simple two-dimensional example. For that, we take $n = 1000$, $p = 2$, $C = 3$, $(\gamma_1, \gamma_2, \gamma_3) = ((1, 1), (2), (2))$ (i.e. one anisotropic and two isotropic Gaussians), $(\pi_1, \pi_2, \pi_3) = (0.4, 0.3, 0.3)$, $\mu_c \sim U([-8, 8]^p)$, $(\lambda_{11}, \lambda_{21}, \lambda_{31}) = (1, .5, .1)$, and $\lambda_{cp}/\lambda_{c1} = 0.01 \forall c$,

<figure><img src="image_2.png" /><figcaption>**Fig. 1** Two-dimensional density fitting with the MPSA-H. Top: evolution of the mixture parameters over the CPEM iterations. The dots represent the data points, the crosses represent the evolution of the centers $\mu_c$ and the ellipses (with increasing opacity) represent the evolution of the covariance matrices. Bottom: evolution of the penalized log-likelihood (left) and the number of parameters (right) over the CPEM iterations.</figcaption></figure>

and we fit an MPSA model with the hierarchical clustering strategy. We display in Figure 1 the scatter plot of the data and the learned mixture model, as well as the evolution of the penalized log-likelihood and the number of parameters.

We see that the clusters are well estimated, as well as their isotropic nature (for the gray and red ones) due to an appropriate MPSA type change occurring after a dozen iterations. Moreover, the penalized log-likelihood increases at each iteration as expected by theory (cf. Theorem 4). The number of parameters decreases over the iterations and ends up being the right one.<fnref n="5" />

## 5.2 Model selection strategies

The aim of this experiment is to study the learning algorithm 2 for different model selection strategies (MPSA-H/R/U/D) in higher dimension. To that extent, we plot in Figure 2 the penalized

<footnote n="5">⁵ $\kappa((1, 1), (2), (2)) = 2 + 5 + 3 + 3 = 13$ (cf. Proposition 1).</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2** Evolution of the penalized log-likelihood and the number of parameters over the iterations of the CPEM algorithm for the four model selection strategies (MPSA-H/R/U/D). Each transparent curve refers to one repetition while the thick opaque curve represents the average over the 10 independent repetitions.</figcaption></figure>

log-likelihood curves and the evolution of the number of parameters over the CPEM iterations, on a synthetic dataset with $n = 1000$, $p = 20$, $C = 3$, $(\gamma_1, \gamma_2, \gamma_3) = ((3, 5, 12), (8, 12), (20, ))$, $(\pi_1, \pi_2, \pi_3) = (0.4, 0.3, 0.3)$, $\mu_c = 0$, $(\lambda_{11}, \lambda_{21}, \lambda_{31}) = (1, .5, .1)$, and $\lambda_{cp}/\lambda_{c1} = 0.01 \forall c$, with 10 independent repetitions. We see that the penalized log-likelihoods are once again monotonous, as expected by theory (cf. Theorem 4).

Let us now comment the results based on our intuitions. All penalized log-likelihood curves seem to converge to the same optimum, after approximately 5 to 20 iterations. MPSA-D takes a much larger number of iterations to converge than the other strategies. However, the learning curves (penalized log-likelihood and number of

parameters) are quite regular and with small variance, compared to the other strategies. We guess that in higher dimensions, MPSA-D may badly behave at first iterations due to the initialization with full covariance matrices, and it can moreover require a large number of iterations to converge, since it can only decrease the model complexity one by one. Therefore, smarter initialization heuristics could be designed, but this is left to future research. MPSA-U also yields quite regular learning curves, although we can observe small oscillations of the number of parameters near the optimum (cf. Remark 7 in Section 4.2). We moreover guess that MPSA-U may get stuck in local optima due to the simplicity of the models at first iterations (Verbockhaven et al. 2024). MPSA-H and MPSA-R are more flexible in terms of model complexity—which may make them attain the optimum faster as we can see in Figure 2. However, we guess that their flexibility may also cause higher variance in performance and parsimony.

To deepen our understanding of the different strategies and compare them to the more standard GMMs, let us run some density modeling and clustering experiments at varying dimensions and distributions.

## 5.3 Penalized log-likelihood benchmark

The goal of this experiment is to compare our new parsimonious GMMs (MPSA-H/R/U/D) to the more-standard GMMs (full and spherical) in terms of penalized log-likelihood (23) $\sum_{i=1}^n \ln(\sum_{c=1}^C \hat{\pi}_c \mathcal{N}(x_i | \hat{\mu}_c, \hat{\Sigma}_c)) - \alpha \kappa(\hat{\gamma})$ (with $\alpha = \ln(n)/2$) under various settings. All datasets but the last one are synthetic, sampled according to the generative model presented in the beginning of this section. The hyperparameters are $C = 3$, $n = 1000$, $(\pi_1, \pi_2, \pi_3) = (0.4, 0.3, 0.3)$, $\mu_c \sim \mathcal{U}([-5, 5]^p)$, $(\lambda_{11}, \lambda_{21}, \lambda_{31}) = (3, 2, 1)$, $\lambda_{cp}/\lambda_{c1} = 0.01 \forall c$, and $p \in \{10, 100\}$. For the first two datasets (MPSA 10 and MPSA 100), one takes $\gamma_1, \gamma_2, \gamma_3 = (1, p-1), (1, 2, p-3), (1, 2, 4, p-7)$ for $p \in \{10, 100\}$. For the third and fourth dataset (Full 10 and Full 100), one takes $(\gamma_1, \gamma_2, \gamma_3) = ((1^p), (1^p), (1^p))$ for $p \in \{10, 100\}$, i.e. full Gaussian models with distinct population eigenvalues. The Skew 100 dataset has $(\gamma_1, \gamma_2, \gamma_3) = ((1^p), (1^p), (1^p))$ for $p = 100$ and a skew-normal distribution (instead of the Gaussian
one used in the previous experiments). The last dataset (Ionosphere) is real (Sigillito et al. 1989), with $C = 2$, $n = 351$ and $p = 34$. We report in Table 1 the average penalized log-likelihood reached at the end of the CPEM algorithm, over 10 independent $n$-samples of the synthetic datasets. The last dataset—which is real—is split using a stratified 10-fold cross-validation.

In the first setting (MPSA 10), all our MPSA models significantly outperform the classical GMMs. Since $p = 10$ and $n = 1000$, one probably has enough samples to estimate the full (anisotropic) covariance parameters, and the additional parsimony of our MPSA models yields a significantly higher penalized log-likelihood. In the second setting (MPSA 100), all our MPSA models again significantly outperform the classical GMMs. Since $p = 100$ and $n = 1000$, one probably does not have anymore enough samples to estimate the full covariance parameters, therefore the full GMM begins to severely drop in performance, while the spherical model is not flexible enough to model the data distribution due to its isotropic nature. In the third setting (Full 10), since the true model is full and we have enough samples, the full GMM is the best estimator, and our MPSA models are on par with it, while the spherical GMM is below. In the fourth setting (Full 100), the true model is still full but we don't have enough samples anymore, so that the full GMM gets worse and the spherical GMM becomes better. Except for MPSA-R which is on par with GMM-S, our other MPSA models are significantly better. The same kind of observation goes for the fifth setting (Skew 100). In the sixth setting (Ionosphere), our MPSA models significantly outperform the other GMMs.

To conclude, our models achieve (almost always significantly) the best likelihood-parsimony tradeoffs in all data regimes, are robust to non-Gaussian distributions and real data, and even significantly outperform the full GMM when the true model is full and we have a small number of samples. All our MPSA models are globally even in terms of final penalized log-likelihood, except for MPSA-R (the simplest and fastest model selection strategy) which seems to be slightly less adapted to high dimensional settings (Full 100 and Skew 100).

The previous experiments used the penalized log-likelihood (23) as a metric to evaluate the quality of the different GMMs. Beyond checking the monotonicity of the penalized log-likelihood over the CPEM iterations, these experiments showed that the MPSA models yield significantly higher penalized log-likelihoods than the classical GMMs—an empirical behavior that cannot be deduced from Theorem 4 alone. We believe that the penalized log-likelihood—as a tradeoff between goodness-of-fit and number of parameters—is a rather good indicator of the generalization properties of the learned models on unseen data, enjoying moreover statistical guarantees such as asymptotic consistency (Keribin 2000; Bai et al. 2018). That being said, it is clear that other interesting questions emanate from the study of the MPSA models, notably the one of density estimation: how well the MPSA models estimate the true underlying probability density. This is the purpose of the next subsection.

## 5.4 Density estimation benchmark

We perform the following density estimation experiment. We generate data from the MPSA generative model, with varying $n$, $p = 15$, $C = 3$, $(\gamma_1, \gamma_2, \gamma_3) = ((15), (5, 10), (5, 5, 5))$, $(\pi_1, \pi_2, \pi_3) = (0.4, 0.3, 0.3)$, $\mu_c = 10$, $(\lambda_{11}, \lambda_{21}, \lambda_{31}) = (3, 2, 1)$, and $\lambda_{cp}/\lambda_{c1} = 0.1 \forall c$. The goal is to estimate the mixture parameters, that are the mixture proportions, means and covariances $(\pi_c, \mu_c, \Sigma_c)_{c=1}^C \in \Theta(\gamma)$, as well as the MPSA types $\gamma \in C(p)^C$. For an increasing number of samples $n$, we display in Figure 3 the covariance errors $\sum_{c=1}^C \|\hat{\Sigma}_c - \Sigma_c\|_F^2$, the type errors $1 - \sum_{c=1}^C \mathbf{1}_{\{\hat{\gamma}_c\}}(\hat{\gamma}_c)/C$, the number of parameters $\kappa(\hat{\gamma})$ and the (averaged) log-likelihood on an independent test set $\ln p(X^{\text{test}} | \hat{\theta})/n^{\text{test}}$.<fnref n="6" /> The component indexes $c$ of the true and estimated mixture densities are matched by minimizing the Euclidean distance between the centers; we check a posteriori that each estimated component is matched to one unique true component—i.e., the map corresponds to a permutation. Each metric is averaged over 10 independent experiments.

<footnote n="6">⁶The Kullback-Leibler divergence—a classical density estimation metric—is not explicit for Gaussian mixture models (Hershey and Olsen 2007).</footnote>
<table><caption>**Table 1** Comparison of the penalized log-likelihoods $\sum_{i=1}^{n} \ln(\sum_{c=1}^{C} \hat{\pi}_c \mathcal{N}(x_i | \hat{\mu}_c, \hat{\Sigma}_c)) - \alpha \kappa(\hat{\gamma})$ reached by the MPSA-H/R/U/D models at the end of the CPEM algorithm, and the standard GMMs under various data distributions. The leading penalized log-likelihoods are in bold.</caption><thead><tr><th>Model</th><th>MPSA 10</th><th>MPSA 100</th><th>Full 10</th><th>Full 100</th><th>Skew 100</th><th>Ionosphere</th></tr></thead><tbody><tr><td>MPSA-H</td><td>-0.65 ± 0.06</td><td>46.00 ± 0.22</td><td>-07.59 ± 0.06</td><td>-093.19 ± 0.42</td><td>-093.98 ± 2.82</td><td>+06.17 ± 1.69</td></tr><tr><td>MPSA-R</td><td>-0.65 ± 0.06</td><td>46.00 ± 0.22</td><td>-07.59 ± 0.06</td><td>-100.02 ± 0.28</td><td>-099.09 ± 0.72</td><td>+04.76 ± 1.67</td></tr><tr><td>MPSA-U</td><td>-0.65 ± 0.06</td><td>46.03 ± 0.22</td><td>-07.59 ± 0.06</td><td>-092.25 ± 0.28</td><td>-092.73 ± 2.00</td><td>+04.59 ± 1.26</td></tr><tr><td>MPSA-D</td><td>-0.65 ± 0.06</td><td>46.03 ± 0.22</td><td>-07.59 ± 0.06</td><td>-093.11 ± 0.41</td><td>-093.98 ± 3.49</td><td>+04.90 ± 2.30</td></tr><tr><td>GMM-F</td><td>-0.92 ± 0.06</td><td>05.37 ± 0.23</td><td>-07.60 ± 0.06</td><td>-104.92 ± 0.23</td><td>-120.22 ± 9.18</td><td>-06.12 ± 1.74</td></tr><tr><td>GMM-S</td><td>-8.40 ± 0.17</td><td>15.11 ± 1.01</td><td>-11.67 ± 0.14</td><td>-100.02 ± 0.28</td><td>-099.67 ± 1.18</td><td>-17.15 ± 0.59</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Fig. 3** Density estimation with several GMMs. Several metrics are reported: the covariance errors, the type errors, the number of parameters and the (averaged) log-likelihood on an independent test set.</figcaption></figure>

To summarize the results, everything behaves as expected: the estimation errors decrease with the number of samples and tend to zero. Interestingly, the covariance estimation errors are significantly smaller for the MPSA models than for the spherical and full GMMs—even for $n = 1000$. This can be explained by looking more closely at the covariance eigenvalue profiles, which is what we do in Figure 4. We see that the sample eigenvalues (represented by the purple curve) are highly dispersed around the population eigenvalues (represented by the black curve). This eigenvalue spreading phenomenon is well known in high dimensions (see, for instance, Johnstone

and Paul (2018)). In contrast, the eigenvalues estimated by the spherical GMM have low variance, but high bias for the second and third mixture component. The MPSA models not only achieve a significantly better estimation of the population eigenvalues, but also well capture the block structure, i.e., the eigenvalue multiplicities, which is not trivial. The eigenvalue plots of Figure 4 therefore well explain the results from Figure 3 in terms of covariance and type estimation.

We can see moreover that the number of parameters increases with the number of samples and seems to converge to the accurate one, testifying of the BIC's asymptotic consistency.
<figure><img src="image_5.png" /><figcaption>**Fig. 4** Covariance eigenvalue estimation with several GMMs for *n* = 1000 (top) and *n* = 200 (bottom). The three mixture components have, respectively, eigenvalue multiplicities (15) (left), (5, 10) (middle) and (5, 5, 5) (right)—plotted in black. The four MPSA strategies cannot be visually distinguished on the first five plots, as they yield relatively similar results. On the sixth plot, MPSA-H/U/D can be distinguished, while MPSA-R cannot be distinguished from GMM-S.</figcaption></figure>

The spherical and full GMMs—respectively the most and least parsimonious GMMs belonging to the MPSA family—naturally bound the MPSA models on the plot of the number of parameters.

Finally, the test log-likelihood curves show that although the full GMM fails at estimating the covariance eigenvalues, it achieves similar density estimation results as the MPSA models, especially with respect to the spherical GMM.

Let us note that the estimation errors for the other mixture parameters—the proportions $\sum_{c=1}^{C}(\hat{\pi}_c - \pi_c)^2$ and the means $\sum_{c=1}^{C}\|\hat{\mu}_c - \mu_c\|_2^2$—were also evaluated: the errors also tend to 0 for all models. We did not display them since the differences between the GMMs were not significant. Indeed, since the clusters are well separated, the proportions and means are well estimated by all models, independently of their covariance structure.

The differences between the model selection strategies (H/R/U/D) are not stringent in terms of estimation, at least in the setting of the experiment and for *n* large enough, where the parameters are well estimated for all strategies. The MPSA-R strategy seems overall slightly worse

than the other strategies, as we can particularly see on the tasks of covariance, type and eigenvalue estimation.

## 5.5 Clustering benchmark

We now perform some clustering experiments, where the learning procedure is similar to the previous experiments, but the reported metric is the adjusted rand index (ARI) (Rand 1971; Hubert and Arabie 1985), which is the classical substitute for the accuracy in clustering problems. This scores varies between (approximately) 0 for a random clustering and (exactly) 1 for a perfect clustering. We add two competitors: HDDC (Bouveyron et al. 2007a) and spectral clustering (Ng et al. 2001). The first three datasets are synthetic. They all have in common that $C = 3$, $(\pi_1, \pi_2, \pi_3) = (0.4, 0.3, 0.3)$, $(\lambda_{11}, \lambda_{21}, \lambda_{31}) = (10, 1, .1)$, $\lambda_{cp}/\lambda_{c1} = 0.01 \forall c$. Now, more specifically, the first dataset (MPSA 10) has $n = 200$, $p = 10$, $\mu_c \sim U([-1, 1]^p)$, and $\gamma_1, \gamma_2, \gamma_3 = (1, p-1), (1, 2, p-3), (1, 2, 4, p-7)$. The second dataset (MPSA 50) has the same characteristics as the previous one but with
$\mu_c = 0$ (mixed clusters, which is the case of the next ones), $p = 50$ and $n = 1000$. The third one (Full 10) has $p = 10$, $n = 200$ and $(\gamma_1, \gamma_2, \gamma_3) = ((1^p), (1^p), (1^p))$.<fnref n="7" /> We then have three real datasets (Ionosphere (Sigillito et al. 1989), Wine (Aeberhard and Forina 1991), Breast (Wolberg et al. 1995)), respectively with $(n, p, C) = (351, 34, 2), (178, 13, 3), (569, 30, 2)$. The average clustering ARIs after 10 independent $n$-samples (or after a stratified 10-fold cross-validation for the real datasets) are reported in Table 2.

The results are obviously harder to comment than in the previous examples, since we evaluate the clustering ARI metric which is not the quantity that is optimized in the GMM learning algorithm (23). Nonetheless, we have several interesting observations. When the clusters centers are separated (MPSA 10), all the GMMs achieve reasonable clustering except for the spherical one which does not fit well the anisotropic clusters, and the spectral clustering which is below in average. Note that the high standard deviation is probably due to a few instances of the 10 repetitions where the randomly-sampled centers are close, or where the randomly-sampled covariance orientations make the clusters harder to separate. When the clusters are mixed (MPSA 50 and Full 10), the performances globally decrease since the setting is more complicated, and the spectral clustering even gets extremely bad (similar to random labeling). Our MPSA models are globally on par with the HDDC models on synthetic data (slightly better in average, but not significantly), although the HDDC models are run several times to tune the Cattell threshold hyperparameter. On real data, the MPSA models surprisingly yield better or even significantly better clustering performances.

The main conclusion of this experiment is that our MPSA models adapt well to several challenging synthetic clustering tasks, taking the most of the two extreme GMMs (full and spherical). A positive surprise is that our MPSA models compare more-than-well to the HDDC models (which they extend), although being fully unsupervised, while HDDC requires hyperparameter tuning. In

<figure><img src="image_6.png" /><figcaption>Fig. 5 Running time benchmark for increasing dimension.</figcaption></figure>

a sense, our models are both robust and remove the need to determine the intrinsic dimensions.

## 5.6 Running times

As seen in the previous experiments, the four strategies are rather comparable in terms of likelihood-parsimony tradeoffs, density estimation and clustering, with the MPSA-R strategy being overall slightly worse than the other strategies. This subsection aims to compare the different strategies from the running time perspective.

We generate data from the MPSA generative model, with $n = 1000$, $C = 3$, $(\gamma_1, \gamma_2, \gamma_3) = ((1, p-1), (1, 2, p-3), (1, 2, 4, p-7)), (\pi_1, \pi_2, \pi_3) = (0.4, 0.3, 0.3)$, $\mu_c = 5$, $(\lambda_{11}, \lambda_{21}, \lambda_{31}) = (3, 2, 1)$, and $\lambda_{cp}/\lambda_{c1} = 0.01 \forall c$. The experiments are repeated independently 10 times. We display the running times of the different GMMs as a function of the dimension $p$ in Figure 5.

The spherical GMM is significantly faster than all the other models—which makes sense since the estimation of the isotropic covariance matrices does not involve any matrix inversion or eigenvalue decomposition. The full GMM is faster than the MPSA models, in average, but it scales like the MPSA-R and MPSA-U strategies. The MPSA-H and MPSA-D strategies, for their part, are much slower, making them less adapted in high dimensions ($p \gtrsim 1000$ for MPSA-H and $p \gtrsim 100$ for MPSA-D). Let us now compare the four MPSA strategies.

<footnote n="7">⁷Note that we choose particularly challenging settings (notably with close—or even exactly similar—centers) to show how our models can adapt the covariance structures, rather than the centers.</footnote>
<table><caption>**Table 2** Comparison of the clustering ARIs (×100) reached by the MPSA-H/R/U/D models at the end of the CPEM algorithm, the standard GMMs, and the HDDC and spectral clustering models, for various challenging data distributions. “—” indicates an average ARI of 0 or less, meaning that the clustering is similar or worse than random. The leading clustering ARIs are in bold.</caption><thead><tr><th>Model</th><th>MPSA 10</th><th>MPSA 50</th><th>Full 10</th><th>Ionosphere</th><th>Wine</th><th>Breast</th></tr></thead><tbody><tr><td>MPSA-H</td><td><b>81 ± 23</b></td><td><b>49 ± 04</b></td><td><b>52 ± 05</b></td><td><b>56 ± 11</b></td><td>44 ± 04</td><td><b>80 ± 06</b></td></tr><tr><td>MPSA-R</td><td><b>81 ± 23</b></td><td><b>49 ± 04</b></td><td><b>53 ± 06</b></td><td>40 ± 09</td><td>42 ± 05</td><td><b>80 ± 06</b></td></tr><tr><td>MPSA-U</td><td><b>81 ± 23</b></td><td><b>48 ± 04</b></td><td><b>50 ± 05</b></td><td>36 ± 06</td><td><b>50 ± 11</b></td><td><b>83 ± 01</b></td></tr><tr><td>MPSA-D</td><td><b>81 ± 24</b></td><td><b>50 ± 04</b></td><td><b>49 ± 05</b></td><td><b>50 ± 12</b></td><td><b>51 ± 12</b></td><td><b>79 ± 06</b></td></tr><tr><td>GMM-F</td><td><b>80 ± 24</b></td><td><b>48 ± 02</b></td><td><b>48 ± 03</b></td><td>32 ± 06</td><td><b>52 ± 12</b></td><td>74 ± 10</td></tr><tr><td>GMM-S</td><td>37 ± 17</td><td>35 ± 10</td><td>43 ± 11</td><td>16 ± 02</td><td>39 ± 02</td><td>67 ± 01</td></tr><tr><td>HDDC</td><td><b>81 ± 24</b></td><td><b>48 ± 02</b></td><td><b>48 ± 04</b></td><td>29 ± 06</td><td>04 ± 04</td><td>29 ± 03</td></tr><tr><td>Spectral</td><td>60 ± 31</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr></tbody></table>

MPSA-R is significantly faster than the other strategies—at the cost of apparently lower modeling performances, as we could see in the previous experiments. MPSA-U is the second fastest strategy, scaling like MPSA-R, although possibly requiring more iterations than MPSA-H, as seen in Figure 2. MPSA-H is the third, as it involves a hierarchical clustering of the eigenvalues, which can become pretty costly in high dimensions (cf. Section 2). MPSA-D is the slowest—by far—as it starts from the full covariance model, and decreases the number of parameters gradually, which can take a lot of iterations to converge if the true model is parsimonious (cf. Figure 2). Hence, the empirical running times rather match the previously built intuitions.

The comparisons between the MPSA and GMMs are to nuance, as our code is certainly not as much optimized as the one of scikit-learn (Pedregosa et al. 2011), from which the classical GMMs are drawn. Moreover, the current version of the code does not include certain numerical tricks such as the previously evoked truncated SVD, which is investigated in Section B.4. We could expect well-engineered MPSA models to be significantly faster than full GMMs, as truncated SVDs are significantly faster than covariance inversion.

As a conclusion to the series of experiments performed in this section, we would probably recommend to adopt the MPSA-U strategy by default in general settings. Indeed, it is the second fastest strategy (scaling as the fastest one and being much faster than the third and fourth), and it has a good control on parsimony (as it is initialized with isotropic covariances and it increases

gradually the number of parameters throughout the CPEM iterations, by finding the most informative gaps in the eigenvalue profiles). In contrast, the other strategies can possibly select overparameterized models in the first iterations, which can negatively impact learning in high dimensions.

# 6 Application to single-image denoising

Beyond traditional learning tasks like clustering, GMMs provide a flexible method for density modeling and therefore come with many other applications. In this section, we briefly describe the application of MPSA models to single-image denoising. More details about the denoising framework can be found in Houdard et al. (2018).

## 6.1 Image denoising with GMMs

Many denoising methods rely on the self-similarity assumption: any small patch extracted from an image has “similar” patches in the same image (Efros and Leung 1999; Buades et al. 2005). This regularity and periodicity prior implies that patches extracted from natural images form well-defined clusters (see Lebrun et al. (2012) for more details). Such a redundancy enables to reduce the noise variance by averaging the similar patches. Under this self-similarity hypothesis, GMMs appear as natural models for the set of patches extracted from one image.

The celebrated work of Zoran and Weiss (2011) runs a GMM on a large database of natural image patches and use it for image denoising, deblurring and inpainting, outperforming all the
traditional methods. However, this approach is not suitable for the task of single-image denoising, where we just have one image as an input, that is moreover noisy; hence the need for parsimonious generative models. Consequently, Houdard et al. (2018) propose to model the set of overlapping patches extracted from one image according to a low-dimensional MPPCA density with a shared mixture noise. The resulting model, HDMI (high-dimensional mixture models for unsupervised image denoising), automatically updates the intrinsic mixture dimensions $q_c$ during the EM based on a fixed hyperparameter $\sigma^2$ that models the variance of the image noise. More precisely, the intrinsic dimensions are chosen such that

$$q_c = \operatorname{argmin}_{q \in [0, p]} \left| \sigma^2 - \frac{1}{p-q} \sum_{j=q+1}^{p} \ell_{cj} \right|, \quad (28)$$

where the $\ell_{cj}$ are the eigenvalues of the local responsibility-weighted covariance matrices $S_c$. The noise level $\sigma^2$ can either be known (supervised denoising) or unknown (unsupervised denoising). In the latter case, it is estimated a posteriori as a hyperparameter, by minimizing the BIC of the learned mixture, similarly to HDDC (Bouveyron et al. 2007a). The next paragraph describes the extension of HDMI with MPSA models.

Consider an image $X \in \mathbb{R}^{s_1 \times s_2}$ that has been corrupted by some isotropic Gaussian noise $E \sim N(0, \sigma^2 I_{s_1 s_2})$. Then the observed image can be written as $\tilde{X} = X + E$. We extract all the overlapping $s \times s$ patches from $\tilde{X}$, from left to right, top to bottom. This yields a dataset $(\tilde{x}_i)_{i=1}^n$ of $n = (s_1 - s + 1)(s_2 - s + 1)$ noisy patches embedded in $\mathbb{R}^p$ with $p = s^2$. A nice illustration of such a dataset is provided in the Figure 1 of Houdard et al. (2018). We assume that it follows an MPSA density (9), with the additional constraints $\lambda_{cd_c} = \sigma^2, \forall c$. We run Algorithm 2 to infer the mixture parameters $\hat{\theta}$. Proposition 5 then gives an explicit and efficient formula to denoise the patches.

**Proposition 5 (Denoising)** Let $\tilde{x} \in \mathbb{R}^p$ represent an $s \times s$ image patch following the noise model $\tilde{x} = x + \epsilon$, where $\epsilon \sim N(0, \sigma^2 I_p)$ is an isotropic Gaussian noise and $x \in \mathbb{R}^p$ is the original patch—sampled from an MPSA model of types $\hat{\gamma}$ and parameters $(\dots \hat{\pi}_c, \hat{\mu}_c, \hat{\Sigma}_c - \sigma^2 I_p \dots)$. Then the denoised patch

$\hat{x} := \mathbb{E}(x | \tilde{x})$ can be expressed as

$$\hat{x} = \sum_{c=1}^{C} w_c \left( \mu_c + \sum_{k=1}^{d_c-1} \left( 1 - \frac{\sigma^2}{\hat{\lambda}_{ck}} \right) \hat{\Pi}_{ck} (x - \mu_c) \right), \quad (29)$$

where $w_c = \mathbb{P}(y = c | \tilde{x})$ (cf. Theorem 2).

Similarly to Section 3.3, one can show that the use of the already-computed eigenvalues and eigenvectors and the substitution of the last eigenspace enables to switch from a computational complexity of $O(Cnp^3)$ (with classical covariance inversion) to $O(Cnp^2q)$ (the most expensive step being the computation of $\hat{\Pi}_{ck} = \hat{Q}_{ck}\hat{Q}_{ck}^\top$). The denoised patches are finally aggregated by performing a uniform reprojection—i.e., positioning back the denoised patches in the image from left to right, top to bottom and averaging them in the overlapping regions—yielding the denoised image $\tilde{X}$.

## 6.2 Denoising experiments

We illustrate the previously presented approach on a 512 $\times$ 512 image with (unknown) noise-deviation $\sigma_{\text{true}} = 30/255$ and patch size $s = 8$. We fit the MPSA-H/R/U models with $C \in \{3, 10, 20\}$ and compare them to the full, spherical, and (supervised) HDMI mixture models.<fnref n="8" /> We evaluate the quality of the denoising using the partial signal-to-noise ratio evaluation metric PSNR := $-10 \log_{10}(\text{MSE}(X, \hat{X}))$, where MSE stands for mean squared error, similarly to Houdard et al. (2018). We repeat the experiment 5 times independently and report the results in Table 3.

We see that the MPSA-H/R/U models always achieve significantly greater PSNR than the classical GMMs. Moreover, while the denoising performance of the MPSA and GMM-S models increase with the number of patches, the ones of GMM-F surprisingly decrease. This can be explained by the fact that increasing the number of clusters reduces the number of patches-per-cluster, so that the full GMM becomes overparameterized, resulting in poorer modeling. The supervised HDMI

<footnote n="8">⁸The HDMI models are learned in a supervised way (i.e. the true noise variance $\sigma_{\text{true}}^2$ is given as an input), which is supposed to give a considerable advantage to HDMI over its competitors. We indeed found out that the results of HDMI were very sensitive with respect to the hyperparameter $\sigma^2$, and knowing that, we did not want to put arbitrary values in the grid search, which would have badly affected its performances.</footnote>
<table><caption>**Table 3** Comparison of the denoising PSNR (↑) reached by MPSA-H/R/U models at the end of the CPEM algorithm, with the ones of standard GMMs and the HDMI model, on the Simpson image, for different numbers of clusters *C*. The leading PSNRs are in bold.</caption><thead><tr><th>Model</th><th><i>C</i> = 3</th><th><i>C</i> = 10</th><th><i>C</i> = 20</th></tr></thead><tbody><tr><td>MPSA-H</td><td><b>31.17</b> ± <b>0.02</b></td><td><b>31.74</b> ± <b>0.04</b></td><td><b>31.97</b> ± <b>0.04</b></td></tr><tr><td>MPSA-R</td><td><b>31.14</b> ± <b>0.02</b></td><td><b>31.72</b> ± <b>0.04</b></td><td><b>31.94</b> ± <b>0.04</b></td></tr><tr><td>MPSA-U</td><td><b>31.17</b> ± <b>0.02</b></td><td><b>31.76</b> ± <b>0.03</b></td><td><b>31.98</b> ± <b>0.03</b></td></tr><tr><td>HDMI(<i>σ</i><sub>true</sub><sup>2</sup>)</td><td><b>31.11</b> ± <b>0.01</b></td><td><b>31.72</b> ± <b>0.03</b></td><td><b>31.90</b> ± <b>0.02</b></td></tr><tr><td>GMM-F</td><td>30.73 ± 0.03</td><td>30.27 ± 0.05</td><td>29.33 ± 0.03</td></tr><tr><td>GMM-S</td><td>22.71 ± 0.01</td><td>25.76 ± 0.13</td><td>27.56 ± 0.04</td></tr></tbody></table>

model is on par or even significantly worse than the MPSA models, which is a good surprise since the MPSA models do not have access to the noise level.

Moreover, we display in Figure 6 the denoised images for the fifth independent repetition of the *C* = 20 experiment. We observe that the GMM-S denoised image is relatively blurry with respect to the other methods, and not much textured. This can be explained by the simplicity of the underlying model for the clusters, which is isotropic Gaussian, so that the denoised patches are essentially sparse convex combinations of a few prototypes.<fnref n="9" /> In contrast, the full GMM has too much degrees of freedom, so that it “over-fits” the noise.<fnref n="10" /> The MPSA models stand in the middle and achieve better denoising, both qualitatively (visually) and quantitatively (via the PSNR metric).

Finally, we propose in Figure 7 a heatmap-like visualization of the clusters directly on the image, for both (supervised) HDMI and MPSA-U models with *C* = 3. More precisely, we map each patch to the cluster it most likely belongs to ($\hat{c} = \arg\max_{c \in [1, C]} \mathbb{P}(y = c | \tilde{x})$), and color this patch according to the number of parameters of the underlying mixture component (which is a PPCA model for HDMI and a PSA model for MPSA-U). We see that the clusters well recover the different textures that are present in the image, the more parsimonious being related to uniform textures (like walls) and the more complex being

related to more contrasted textures (like edges and boundaries). Moreover, we see that the HDMI model requires much more parameters to model the patches (mostly for the textured zones but also for the uniform zones). More precisely, one has:

$$ \gamma_{\text{HDMI}} := ((1^9, 55), (1^{10}, 54), (1^{39}, 25)), $$

$$ \gamma_{\text{MPSA}} := ((1^5, 59), (1^6, 58), (1^9, 2, 1, 4, 2, 2, 44)), $$

which yields by virtue of Proposition 1:

$$ \kappa(\gamma_{\text{HDMI}}) = 2 + 605 + 660 + 1820 = 3087, $$

$$ \kappa(\gamma_{\text{MPSA}}) = 2 + 375 + 434 + 1140 = 1951. $$

This means that the uniform textures are attributed by HDMI an intrinsic dimension of respectively 9 and 10, while MPSA-U only requires respectively 5 and 6 dimensions. In contrast, the complex textures have an intrinsic dimension of 39 with the HDMI model, while they have an intrinsic dimension of 20(= 64 – 44) with the MPSA-U model.<fnref n="11" /> Finally, we find that the heatmap obtained with MPSA-U is cleaner visually than the one of HDMI, probably due to the added parsimony.

<footnote n="9">⁹The denoising formula (29) applied in the context of a spherical GMM indeed yields $\hat{x} = \sum_{c=1}^{C} w_c \mu_c$.</footnote>

<footnote n="10">¹⁰Examining the denoising formula (29), we see that if the estimated noise $\sigma^2$—consisting in the proportion-weighted average of the smallest eigenvalues of the $S_c$—is small with respect to the other eigenvalues of the $S_c$, then one roughly has $\hat{x} \approx \tilde{x}$.</footnote>

<footnote n="11">¹¹A legitimate question might be: does an HDMI model of similar complexity as the MPSA-U model yield smaller PSNRs? The only way to reduce the number of parameters of the HDMI model being to increase the estimated noise (cf. Eq. (28)), we tried several values of $\sigma$. For $\sigma = 40/255$, one gets the types (1, 63), (1, 63), (1, 63) and the PSNR 26.8; for $\sigma = 35/255$, one gets the types (1, 63), (1, 63), (1², 62) and the PSNR 27.8; for $\sigma = 32/255$, one gets the types (1, 63), (1, 63), (1⁷, 57) and the PSNR 30.7; for $\sigma = 30/255(=\sigma_{\text{true}})$, one gets the types (1⁹, 55), (1¹⁰, 54), (1³⁹, 25) and the PSNR 31.1. We notice an extreme sensitivity of the selected model complexity and PSNRs with respect to the choice of hyperparameter $\sigma$, but it seems that with similar numbers of parameters, HDMI and MPSA-U behave relatively similarly.</footnote>
<figure><img src="image_7.png" /><figcaption>Fig. 6 Denoised images and PSNR with different parsimonious GMMs with $C = 20$.</figcaption></figure>

*Remark 11* (Influence of *s* and *C* on denoising) In all the experiments of this section, the patch size is set to *s* = 8 and the number of components *C* varies heuristically according to the denoising task. We discuss further in depth the influence of *s* and *C* on the denoising performances and running times in Section B.3. Similar experiments can be found in Houdard et al. (2018).

# 7 Conclusion

We introduced a new family of parsimonious Gaussian mixture models with piecewise-constant covariance eigenvalue profiles. These were motivated by the important parsimony achieved by equalizing similar eigenvalues (Szwagier and Penec 2025), going beyond classical low-rank Gaussian mixture models (Tipping and Bishop 1999a; Bouveyron et al. 2007a). While the unsupervised learning of the mixture parameters was easily enabled via an EM algorithm, the automatic inference of the groups of eigenvalues to equalize was more technical. It was finally achieved via a *provably-monotonous componentwise penalized EM* algorithm. The numerical experiments

showed the superior likelihood-parsimony trade-offs achieved by our models when compared to the two extreme cases: the full GMM (with distinct covariance eigenvalues) and the spherical GMM (with equal covariance eigenvalues), and the effective application of our models to clustering, density modeling and denoising tasks.

Beyond our extensive experiments to analyze the learning algorithm, one could investigate more data regimes, variants and applications. We notably believe that our new CPEM algorithm has the potential to be applied to many classical low-rank mixture models, removing the need for post-training hyperparameter selection. Instead of an EM algorithm, one could also directly optimize the penalized log-likelihood with gradient descent techniques (Hosseini and Sra 2020), drawing from the recent relaxation of parsimony of Szwagier et al. (2025). The MPSA models could also be naturally applied to the numerous generative modeling tasks for which classical GMMs are often used. Their expressivity could be improved by drawing from the recent advances in GMMs and machine learning (Richardson and Weiss 2018; Loconte et al. 2023; Byrski et al. 2024). Robustness could be achieved by integrating our principles into
<figure><img src="image_8.png" /><figcaption>**Fig. 7** Heatmap-based representation of the clustering achieved by the HDMI and MPSA-U models on the patches of the noisy Simpson image, for *C* = 3. The patches are colored according to the number of parameters of the mixture component they belong to.</figcaption></figure>

the recent mixtures of high-dimensional elliptical distributions of Oudoumanessah et al. (2024).

Other potential applications of our parsimonious models are *manifold learning* and *surface reconstruction*. Indeed, it is often assumed that the observed data has been sampled from a low-dimensional manifold (up to isotropic Gaussian noise). The density could be fit with our MPSA

models, and the underlying shape could be recovered with classical surface reconstruction methods. The types of the mixture components could tell something about the local intrinsic dimensions of the underlying shape (Akhøj et al. 2023; Buet and Pennec 2023).

Finally, since the PSA models have been originally introduced to improve the interpretation of principal components with close eigenvalues,
one should definitely investigate the interpretation of the principal subspaces found by the MPSA models. We conclude this paper with a simple example on image patches. As discussed and experimented in Section 6, image patches are classically represented with Gaussian mixtures, and we believe that they carry some invariances that could be well modeled with piecewise-constant eigenvalue profiles. Indeed, many classical works have noticed that the principal components in natural image patch datasets have repeated eigenvalues, accounting for shared frequencies (see, for instance, Olshausen and Field (1996, FIG. 1), Mairal et al. (2010, Figure 5(a)) and Zoran and Weiss (2011, Figure 6)). Therefore, the MPSA models could reveal these hidden symmetries by detecting the multiple covariance eigenvalues in the mixture components, similarly to Szwagier and Pennec (2025). We illustrate this perspective by fitting an MPSA-U model with $C = 3$ on patches extracted from the clean Simpson image. The output model types are quite interesting, e.g. $\gamma = (2, 3, 4, 5, ...)$ for the first cluster. We display in Figure 8 the first few principal components associated with the first cluster, gathered according to their eigenvalue multiplicities, and then sample some points from the associated principal subspaces. This essentially boils down to generating convex combinations of the principal components with similar eigenvalues. We observe some interesting patterns within the principal subspaces, which can somewhat be interpreted as feature subspaces with (limited) frequency invariance.

**Acknowledgments.** This work was supported by the ERC grant #786854 G-Statistics from the European Research Council under the European Union's Horizon 2020 research and innovation program and by the French government through the 3IA Côte d'Azur Investments ANR-23-IACL-0001 managed by the National Research Agency.

# Appendix A Proofs

This section provides the proofs of the propositions and theorems from the main text.

## A.1 Proposition 1: number of parameters for the MPSA models

*Proof.* The parameters of the MPSA model of types $\gamma$ are $\pi_1, \mu_1, \Sigma_1, \dots, \pi_C, \mu_C, \Sigma_C$, with the constraints $\pi_c \in \mathbb{R}_{\ge 0}, \sum_{c=1}^C \pi_c = 1, \mu_c \in \mathbb{R}^p$ and $\Sigma_c \in \text{Sym}_+(\gamma_c)$. The number of free parameters for $\pi_1, \dots, \pi_C$ (belonging to the standard simplex) is $C-1$, since the last mixture proportion $\pi_C$ is automatically determined from the $C-1$ other ones, via $\pi_C = 1 - \sum_{c=1}^{C-1} \pi_c$. Then for each $c \in \{1, C\}$, the parameters $\mu_c$ and $\Sigma_c$ describe a PSA model of type $\gamma_c$. Therefore, the number of free parameters is $\kappa(\gamma_c)$, as defined in Equation (2). Since neither the means nor the covariance matrices are shared between the mixture components, the number of parameters for $\mu_1, \Sigma_1, \dots, \mu_C, \Sigma_C$ is simply the sum $\sum_{c=1}^C \kappa(\gamma_c)$. Eventually, since the mixture proportions are independent from the means and covariance matrices, the total number of parameters is $\kappa(\gamma) = C - 1 + \sum_{c=1}^C \kappa(\gamma_c)$. $\square$

## A.2 Theorem 2: E-step of the EM

*Proof.* The posterior log-probabilities can be rewritten using Bayes' theorem:

$$ \ln \mathbb{P}(y = c | x, \theta) = \ln (\pi_c \mathcal{N}(x | \mu_c, \Sigma_c)) - \ln \left( \sum_{c'=1}^{C} \pi_{c'} \mathcal{N}(x | \mu_{c'}, \Sigma_{c'}) \right). \quad (A1) $$

Let us now develop the operand using the eigen-decomposition of the sample covariance matrices, $\Sigma_c = \sum_{k=1}^{d_c} \lambda_{ck} \Pi_{ck}$, where $\lambda_{ck}$ are the (ordered decreasing) eigenvalues of $\Sigma_c$ and $\Pi_{ck}$ are the associated eigenspaces, of dimension $\gamma_{ck}$. One has $\ln |\Sigma_c| = \sum_{k=1}^{d_c} \gamma_{ck} \ln \lambda_{ck}$, and via classical trace tricks, one gets

$$ d_{\Sigma_c}^2(x, \mu_c) := (x - \mu_c)^{\top} \Sigma_c^{-1} (x - \mu_c), \quad (A2) $$

$$ = \text{tr} \left( \Sigma_c^{-1} (x - \mu_c) (x - \mu_c)^{\top} \right), \quad (A3) $$

$$ = \sum_{k=1}^{d_c} \text{tr} \left( \lambda_{ck}^{-1} \Pi_{ck} (x - \mu_c) (x - \mu_c)^{\top} \right), \quad (A4) $$

$$ = \sum_{k=1}^{d_c} \frac{(x - \mu_c)^{\top} \Pi_{ck} (x - \mu_c)}{\lambda_{ck}}. \quad (A5) $$

Consequently, one has

$$ \ln(\pi_c \mathcal{N}(x | \mu_c, \Sigma_c)) = \ln \pi_c - \frac{p}{2} \ln 2\pi - \frac{1}{2} \sum_{k=1}^{d_c} \left( \gamma_{ck} \ln \lambda_{ck} + \frac{(x - \mu_c)^{\top} \Pi_{ck} (x - \mu_c)}{\lambda_{ck}} \right). \quad (A6) $$
<figure><img src="image_9.png" /><figcaption>Samples from Principal Subspaces</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8 Illustration of the principal subspaces found by fitting an MPSA-U model to patches extracted from the clean Simpson image. The inferred eigenvalue multiplicities for the first component are $\gamma_1 = (2, 3, 4, 5, \dots)$, which means that the principal subspaces #1, #2, #3, #4 have respective dimensions 2, 3, 4 and 5. We display some samples from these principal subspaces (more precisely, convex combinations of the associated principal components). These can be interpreted as feature subspaces with (limited) frequency invariance, which are reminiscent of the findings of Szwagier and Pennec (2025) about the curse of isotropy.</figcaption></figure>

Let us define the cost function $K_c$ similarly to Bou-veyron et al. (2007a) via the following relationship: $\ln(\pi_c \mathcal{N}(x | \mu_c, \Sigma_c)) = -\frac{K_c(x)+p \ln 2\pi}{2}$. Re-injecting $K_c$ in the posterior probabilities yields

$$ \mathbb{P}(y = c | x, \theta) = \frac{\exp\left(-\frac{K_c(x)+p \ln 2\pi}{2}\right)}{\sum_{c'=1}^{C} \exp\left(-\frac{K_{c'}(x)+p \ln 2\pi}{2}\right)}, \quad (A7) $$

$$ = \frac{1}{\sum_{c'=1}^{C} \exp\left(\frac{K_c(x)-K_{c'}(x)}{2}\right)}. \quad (A8) $$

Re-injecting $K_c$ in the log-likelihood yields

$$ \text{LL}(\theta) = \sum_{i=1}^{n} \ln \left( \sum_{c=1}^{C} \pi_c \mathcal{N}(x_i | \mu_c, \Sigma_c) \right), \quad (A9) $$

$$ = \sum_{i=1}^{n} \ln \left( \sum_{c=1}^{C} \exp \left( -\frac{K_c(x_i) + p \ln 2\pi}{2} \right) \right). \quad (A10) $$

Remark 12 The second term in (A1) is independent of $c \in [1, C]$. Therefore, only the first term intervenes in the cluster assignment problem, which can be rewritten as $y = \arg \max_{c \in [1, C]} \ln(\pi_c \mathcal{N}(x | \mu_c, \Sigma_c))$. Then

one has $y = \arg \min_{c \in [1, C]} K_c(x)$, which makes the implementation more efficient.

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

### A.3 Theorem 3: M-step of the EM

*Proof* The optimal mixture proportions are

$$ \hat{\pi}_1, \dots, \hat{\pi}_C = \arg \max_{\substack{\pi_c \ge 0 \\ \sum_{c=1}^C \pi_c = 1}} \sum_{c=1}^C \left( \ln \pi_c \sum_{i=1}^n t_{ci} \right), \quad (A11) $$

whose solutions are given by $\hat{\pi}_c = \frac{1}{n} \sum_{i=1}^n t_{ci}$. For each $c \in [1, C]$, let us write $S_c := \frac{1}{n\hat{\pi}_c} \sum_{i=1}^n t_{ci} (x_i - \mu_c)(x_i - \mu_c)^\top \in \text{Sym}(p)$. Then one can rewrite the objective function as

$$ \mathbb{E}_{y|x, \theta} \left[ \ln \mathcal{L}_c \left( (\hat{\pi}_c, \mu_c, \Sigma_c) \right) \right] = \text{const} \\ - \frac{1}{2} \sum_{c=1}^{C} n \hat{\pi}_c \left( \ln |\Sigma_c| + \text{tr} \left( \Sigma_c^{-1} S_c \right) \right). \quad (A12) $$

Consequently, the M-step on the mean and covariance parameters can be decomposed into the following
independent minimization problems for $c \in [1, C]$:

$$ \hat{\mu}_c, \hat{\Sigma}_c = \arg \min_{\mu_c \in \mathbb{R}^p} \ln |\Sigma_c| + \text{tr}(\Sigma_c^{-1} S_c). \quad (A13) \\ \Sigma_c \in \text{Sym}_+(\gamma_c) $$

The solutions to such problems are given in the maximum likelihood estimation theorem of Szwagier and Pennec (2025). One has $\hat{\mu}_c = \frac{1}{n\pi_c} \sum_{i=1}^n t_{ci} x_i$ and $\hat{\Sigma}_c = \sum_{k=1}^{d_c} \hat{\lambda}_{ck}(t_c, \gamma_c) \hat{\Pi}_{ck}(t_c, \gamma_c)$, where $\hat{\lambda}_{ck}$ and $\hat{\Pi}_{ck}$ are the block-averaged eigenvalues and eigenspaces of $S_c$. $\square$

## A.4 Theorem 4: monotonicity of the CPEM

*Proof* Let $\mathbb{E}_{y|x,\theta}[\ln \mathcal{L}_c(\cdot)]$ be the expected complete-data log-likelihood (16). Then one has

$$ \mathbb{E}_{y|x,\theta}[\ln \mathcal{L}_c(\theta(\gamma'))] = \sum_{c=1}^{C} \Psi_c(\gamma_c') + \text{const}, \quad (A14) $$

and

$$ \kappa(\gamma') = \sum_{c=1}^{C} \kappa(\gamma_c') + \text{const.} \quad (A15) $$

Therefore, summing the inequality (24) from the theorem statement over $c \in [1, C]$ and adding the constants yields

$$ \begin{aligned} \mathbb{E}_{y|x,\theta}[\ln \mathcal{L}_c(\theta(\gamma'))] & - \alpha \kappa(\gamma') \\ & \geq \mathbb{E}_{y|x,\theta}[\ln \mathcal{L}_c(\theta(\gamma))] - \alpha \kappa(\gamma), \end{aligned} \quad (A16) $$

where $\theta(\gamma) := (\hat{\pi}_c(t_c), \hat{\mu}_c(t_c), \hat{\Sigma}_c(t_c, \gamma_c))_c \in \Theta(\gamma)$. Additionally, by definition of $\theta(\gamma)$, that is $\theta(\gamma) := \arg \max_{\theta' \in \Theta(\gamma)} \mathbb{E}_{y|x,\theta}[\ln \mathcal{L}_c(\theta')]$, one has

$$ \begin{aligned} \mathbb{E}_{y|x,\theta}[\ln \mathcal{L}_c(\theta(\gamma))] & - \alpha \kappa(\gamma) \\ & \geq \mathbb{E}_{y|x,\theta}[\ln \mathcal{L}_c(\theta)] - \alpha \kappa(\gamma). \end{aligned} \quad (A17) $$

The two previous equations (A16, A17) simply mean that a classical EM iteration for the MPSA model of types $\gamma$ increases the expected complete-data penalized log-likelihood $\mathbb{E}_{y|x,\theta}[\ln \mathcal{L}_c(\cdot)] - \alpha \kappa(\gamma)$, and that an adapted change of types from $\gamma$ to $\gamma'$ increases this quantity even more. Using Gibbs' inequality as classically done in the proof of the EM monotonicity (Dempster et al. 1977), one gets

$$ \begin{aligned} \ln \mathcal{L}(\theta(\gamma')) - \ln \mathcal{L}(\theta) \\ \geq \mathbb{E}_{y|x,\theta}[\ln \mathcal{L}_c(\theta(\gamma'))] - \mathbb{E}_{y|x,\theta}[\ln \mathcal{L}_c(\theta)]. \end{aligned} \quad (A18) $$

We conclude the proof by subtracting $\alpha(\kappa(\gamma') - \kappa(\gamma))$ on both sides and using (A16, A17):

$$ \begin{aligned} (\ln \mathcal{L}(\theta(\gamma')) - \alpha \kappa(\gamma')) - (\ln \mathcal{L}(\theta) - \alpha \kappa(\gamma)) \geq \\ (\mathbb{E}_{y|x,\theta}[\ln \mathcal{L}_c(\theta(\gamma'))] - \alpha \kappa(\gamma')) \\ - (\mathbb{E}_{y|x,\theta}[\ln \mathcal{L}_c(\theta)] - \alpha \kappa(\gamma)) \geq 0. \end{aligned} \quad (A19) \quad \square $$

## A.5 Proposition 5: GMM-based single-image denoising

*Proof* The Proposition 1 of Houdard et al. (2018) proves that $\mathbb{E}(x | \tilde{x}) = \sum_{c=1}^{C} w_c \psi_c(\tilde{x})$ with $\psi_c(\tilde{x}) := \mu_c + (I_p - \sigma^2 \Sigma_c^{-1})(\tilde{x} - \mu_c)$. Then, using the eigendecomposition of the sample covariance matrices, $\Sigma_c = \sum_{k=1}^{d_c} \lambda_{ck} \Pi_{ck}$, where $\lambda_{ck}$ are the (ordered decreasing) eigenvalues of $\Sigma_c$ and $\Pi_{ck}$ are the associated eigenspaces, one can rewrite:

$$ \begin{aligned} \psi_c(\tilde{x}) &= \mu_c + \left( I_p - \sigma^2 \sum_{k=1}^{d_c} \lambda_{ck}^{-1} \Pi_{ck} \right) (\tilde{x} - \mu_c) \quad (A20) \\ &= \mu_c + \left( \sum_{k=1}^{d_c} \Pi_{ck} - \sigma^2 \sum_{k=1}^{d_c} \lambda_{ck}^{-1} \Pi_{ck} \right) (\tilde{x} - \mu_c), \quad (A21) \end{aligned} $$

$$ = \mu_c + \sum_{k=1}^{d_c} \left(1 - \frac{\sigma^2}{\lambda_{ck}}\right) \Pi_{ck}(\tilde{x} - \mu_c), \quad (A22) $$

$$ = \mu_c + \sum_{k=1}^{d_c-1} \left(1 - \frac{\sigma^2}{\lambda_{ck}}\right) \Pi_{ck}(\tilde{x} - \mu_c). \quad (A23) $$

Let us insist on the last line, which is due to the constraints $\lambda_{cd_c} = \sigma^2$ for all $c \in [1, C]$ and which shows that the patches are actually projected to the low-dimensional subspaces $V_c := \text{Span}(\sum_{k=1}^{d_c-1} \Pi_{ck})$. This represents an important computational gain in high dimensions (from $O(p^3)$ to $O(p^2)$) for the same reasons as outlined in Section 3.3. $\square$

# Appendix B Experiments

## B.1 Penalty factor $\alpha$

In all the experiments, we set the penalty factor to $\alpha = \ln(n)/2$. This corresponds to a BIC-like penalized log-likelihood (23), which therefore enjoys theoretical foundations and asymptotic consistency guarantees (Schwarz 1978; Keribin 2000; Bai et al. 2018). Notably, as we can see it in the density estimation experiments of Section 5.4, the estimated mixture parameters converge to the true ones when the number of samples goes to infinity. We also chose this criterion because this is a very classical one in model-based clustering (see, for instance, Scrucca et al. (2016); Bouveyron et al. (2019)).

Alternatives criteria could be the Akaike information criterion (AIC) (Akaike 1974) or the integrated completed likelihood (ICL) (Biernacki et al. 2000). These criteria might not enjoy the
<figure><img src="image_12.png" /><figcaption>**Fig. B1** Influence of the regularization hyperparameter $\alpha$ on the number of parameters of the inferred mixture models. The true number of parameters is $\kappa(\gamma) = 421$.</figcaption></figure>

same theoretical guarantees and result in different empirical results. Notably, the AIC ($\alpha = 1$) is known to favor overparameterized models with respect to the BIC ($\alpha = \ln(n)/2$). While AIC and BIC directly express under the form of Equation (23)—i.e., as a log-likelihood penalized by the number of samples—ICL has a slightly different expression, which makes it non-integratable in the current framework, especially if we want monotonicity guarantees such as the ones of Theorem 4.

To account for the influence of the penalty factor, we report in Figure B1 the evolution of the number of parameters of the learned mixture model (via the CPEM algorithm) as a function of $\alpha$.<fnref n="12" /> We see that the number of parameters indeed decreases with the penalty factor $\alpha$, going from the highest complexity (corresponding to the full GMM) to the lowest one (corresponding to the spherical GMM). The true number of parameters lies in between the ones selected by the BIC (underparameterized) and the AIC (overparameterized).

While the AIC and the BIC do not depend on the dimension $p$, alternative criteria such as the AICc (Hurvich and Tsai 1989)—a small-sample

version of the AIC—do. Some additional empirical results about the impact of the model selection criterion on the number of parameters of the PSA models are given in the appendix of Szwagier and Pennec (2025). The AICc tends to select even more parsimonious models than the BIC for small-to-moderate sample sizes.

## B.2 Number of components C

In contrast to many works from the GMM literature that reduce the number of mixture parameters by reducing the number of components $C$, our work controls parsimony through the covariance structure of the mixture components. While $C$ has a *linear* influence on the number of parameters, the eigenvalue multiplicities $\gamma_{ck}$ have a *quadratic* influence on the number of parameters (see Equation (2)). Therefore, in high dimensions, we believe that it is much more interesting to reduce the number of parameters via the eigenvalue multiplicities than via the number of components (unless $C$ is very large).

Except for the image denoising experiments in Section 6, this paper voluntarily involves experiments where the number of components $C$ is well specified—i.e., equal to the true number of clusters. This choice is motivated by the willingness to highlight the impact of the covariance structure on the learned mixtures—that is the main message of the paper—rather than the impact of the number of clusters.

In the current version of the code, the number of components $C$ can be updated during learning via several simple methodologies, for instance removing empty clusters, or the ones with less than a prespecified number of points. $C$ can also be estimated *a posteriori*, as commonly done, via model selection criteria such as the BIC (see, for instance, Bouveyron et al. (2019, Section 2.6)). As an illustration, we compare in Figure B2 the penalized log-likelihood of an MPSA-U model with $C = 2$, $C = 3$ and $C = 4$, under the setting of the previous experiment—with three clusters. As we can see, the underfitted model ($C = 2$) always has a significantly smaller penalized log-likelihood than the other models. Surprisingly, the overfitted model ($C = 4$) achieves a significantly larger penalized log-likelihood than the true one ($C = 3$) for $n \lesssim 100$. Looking more closely at the estimated parameters, it seems that under

<footnote n="12">¹²The experimental details are the following ones: the data is generated from an MPSA density with $n = 1000$, $p = 20$, $C = 3$, $(\gamma_1, \gamma_2, \gamma_3) = ((5, 15), (5, 5, 10), (5, 5, 5, 5))$, $(\pi_1, \pi_2, \pi_3) = (0.4, 0.3, 0.3)$, $\mu_c = 10$, $(\lambda_{11}, \lambda_{21}, \lambda_{31}) = (3, 2, 1)$, and $\lambda_{cp}/\lambda_{c1} = 0.1 \forall c$. The number of parameters is therefore $\kappa(\gamma) = 421$. The experiments are repeated independently 10 times.</footnote>
<figure><img src="image_13.png" /><figcaption>**Fig. B2** Influence of the number of components *C* on the penalized log-likelihood, under the MPSA-U model. The true number of clusters is *C* = 3.</figcaption></figure>

the overfitted model, two clusters are well esti-
mated, while the third one is split in two. This
overparameterization naturally yields a larger log-
likelihood, which is apparently not compensated
by the penalty on the number of parameters—
notably since the two overfitting components have
lower number of parameters than the true one.
Therefore, penalized log-likelihood criteria such as
the BIC may have some limitations in selecting
the true number of components, when the num-
ber of samples is small. For $n \gtrsim 100$, the true
model ($C = 3$) gets a significantly larger penal-
ized log-likelihood and is therefore selected over
the overfitted one ($C = 4$), testifying of the BIC's
asymptotic consistency.

Beyond the two simple ideas proposed above
(removing unpopulated clusters and using model
selection criteria), we believe that the choice of the
number of clusters could be thoroughly addressed
via Dirichlet-type priors on the mixture propor-
tions. Such priors are commonly used in the
context of overfitted mixture models (see, notably,
the classical works of Figueiredo and Jain (2002);
Rousseau and Mengersen (2011); Malsiner-Walli
et al. (2016)). They enable to automatically
remove underpopulated clusters within the M-
step, without losing the theoretical monotonicity
guarantees of the EM algorithm.

B.3 Patch size *s* and number of
components *C* in denoising

Except from Table 3, the main text’s denoising
experiments do not assess the impact of hyper-
parameters that are unrelated to the covariance
structure on the results. This is the goal of this
subsection.

The size of the patches naturally influences the
denoising results. Intuitively, if the patches are too
small, then they all look similar, while if they are
too large, then they all look different—especially
in textured zones. This has notable repercussions
on the denoising performances, since the investi-
gated methods precisely rely on a faithful model-
ing of the groups of similar patches extracted from
one image (see Section 6.1). Figure B3 reports the
influence of the patch size on the quality of denois-
ing (PSNR). We can see that it first increases
with the patch size, up until approximately 8 × 8
(which is the size chosen in the experiments of
Section 6) before decreasing. We can also see that
the computing time naturally increases with the
patch size, since the dimension of the dataset
(s²) increases quadratically with the patch size s.
Therefore, beyond the denoising PSNR, the run-
ning time is another metric to take into account
for choosing the patch size.

Similarly, the number of components *C* does
influence the denoising results. Intuitively, the
goodness-of-fit increases with the number of com-
ponents *C*. If the latter is too small, then we are
in a situation of underfitting; the noise is over-
estimated and the denoising leads to a blurry
image. In contrast, if the number of components
is too large, then we are in a situation of overfit-
ting; the noise is underestimated and the denoising
does not modify much the input patch, result-
ing in a sharp yet noisy image; the limit case
being for *C* = *n*, where each point is in its own
cluster, the estimated noise is zero and the denois-
ing leaves the noisy patches unchanged. These
intuitions are (partially) confirmed empirically in
Figure B3, which reports the influence of the num-
ber of clusters on the PSNR and running times.<fnref n="13" />
Note that the running time increases linearly with
*C*, therefore it might impact the practical choice

<footnote n="13">¹³The decrease in denoising PSNR for large *C* has been
observed in Houdard et al. (2018, Section 4.1), but it appears
to require an extremely large number of components to be
significant, therefore we chose not to report it in Figure B3.</footnote>
<figure><img src="image_14.png" /><figcaption>**Fig. B3** Influence of the patch size *s* (top) and the number of components *C* (bottom) on the denoising quality (PSNR) and running time, under the MPSA-U model.</figcaption></figure>

of *C* beyond the denoising performance. Notably, MPSA-U already reaches satisfying PSNRs for *C* = 5.

## B.4 Truncated SVD

The code for MPSA inference has been partially optimized in order to run in practical time. However, it could largely enjoy more engineering, especially for very large dimensions. One potential improvement is the truncated SVD (TSVD) evoked in Section 3.3.

The current M-step for MPSA inference implies a full EVD of the responsibility-weighted covariance matrices $S_c$. Indeed, all the *eigenvalues* are needed to select the appropriate eigenvalue multiplicities. However, not all the *eigenvectors* are needed, but just the leading ones (as explained

in Section 3.3). Hence, an idea can be to first compute all the eigenvalues (without the associated eigenvectors), then choose the appropriate eigenvalue multiplicities $\hat{\gamma}_c = (\hat{\gamma}_{c1}, \dots, \hat{\gamma}_{cd_c})$ and finally compute the $\hat{q}_c = p - \hat{\gamma}_{cd_c}$ leading eigenvectors (associated to the $\hat{q}_c$ largest eigenvalues) via truncated SVD. While such a trick should not change the outcome in theory—with respect to a full EVD—it actually does in practice, since the truncated SVD relies on different solvers from the classical EVD. More details can be found in the documentation of the scipy package (Virtanen et al. 2020): see `scipy.linalg.eigh` for the eigenvalue decomposition and `scipy.sparse.linalg.svds` for the truncated SVD. The randomized SVD of Halko et al. (2011) is another popular solver.

Another simple idea is to impose an upper bound on the intrinsic dimensions of the mixture components (i.e. $\hat{q}_c \le q_{\max} \forall c$). In that case, one can directly run a truncated SVD with $q_{\max}$ components without the need to compute all the eigenvalues a priori. This approach can be expected to be more efficient than the previous one, as it does not require to compute all the eigenvalues first. However, imposing such an upper bound on the intrinsic dimensions restricts the set of candidate models and may therefore lead to suboptimal results.

To compare these different approaches, we sample $n = 1000$ points from a multivariate Gaussian distribution with zero mean and covariance diag (10, 8, 6, 4, 2, 0.1$^{p-5}$), for varying $p$. We display the running times of the different approaches in Figure B4. We see that the full SVD on the data matrix is overall slower than the other methodologies, including the full EVD on the covariance matrix. For very large dimensions ($p \gtrsim 1000$), the methods involving first an eigenvalue decomposition without the eigenvectors and second a truncated SVD become significantly faster than the full EVD. The methods directly imposing an upper bound on the intrinsic dimensions are, for their part, significantly faster than the previous ones. Hence, we recommend integrating TSVD-based approaches in very large dimensions and directly imposing an upper bound on the intrinsic dimensions if speed and parsimony matter more than precise inference.
<figure><img src="image_15.png" /><figcaption>**Fig. B4** Comparison of the running times for different eigenvalue decomposition methodologies. EVD: full eigenvalue decomposition of the covariance matrix. SVD: full singular value decomposition of the data matrix. TSVD: truncated SVD of the data matrix, with five components (random: randomized SVD from Halko et al. (2011); arpack: scipy solver which uses an Arnoldi iteration). EVD + TSVD: eigenvalue decomposition of the covariance matrix (without the eigenvectors) followed by TSVD with five components.</figcaption></figure>

## B.5 Benchmark

In our benchmark experiments, we mostly focus on comparing our MPSA models to the classical full and spherical GMMs—and to the less-classical HDDC (Bouveyron et al. 2007a) and HDMI (Houdard et al. 2018). This seemingly restrictive choice is actually well motivated: we want to assess the impact of equalizing close eigenvalues in covariance matrices on mixture modeling. We do not aim to achieve state-of-the-art results on density estimation, clustering and denoising tasks. Therefore, we believe that the best way to convey the main message of the paper (that is the interest of equalizing close eigenvalues in high-dimensional, low-sample settings) is to compare our models to baseline models that are particular cases of ours, and whose distinctive feature is the covariance structure *alone*. Indeed, if we compared our models to other baseline models (for density estimation, clustering or denoising), then the difference in performance could be explained by other features than the covariance structure, which would drown the core message of the paper.

## B.6 Additional details

This subsection provides more details about the experiments. The code is written in Python 3.9 and run on an Intel® Core™ i7-10850H CPU with 16 GB of RAM. The code for the EM algorithm is a modification of Mathieu Fauvel's code (license: GPL-3.0), itself inspired from the R package of Bergé et al. (2012) (license: GPL-2). Some parts of the code related to PSA inference and model selection were inspired from the code of Szwagier and Pennec (2025) (license: MIT). The code for the denoising pipeline is a modification of Antoine Houdard's code (license: MIT).

The used packages are the following: scikit-learn (Pedregosa et al. 2011) (license: BSD 3), numpy (Harris et al. 2020) (license: BSD), scipy (Virtanen et al. 2020) (license: BSD 3), matplotlib (Hunter 2007) (license: custom), ucimlrepo (license: MIT), scikit-image (Walt et al. 2014) (license: BSD 3). The code for sampling from a skew-normal distribution comes from Gregory Gundersen's blog post (license: not specified).

Some datasets come from the UCI Machine Learning repository: wine (Aeberhard and Forina 1991) (license: CC BY 4.0), breast (Wolberg et al. 1995) (license: CC BY 4.0) and ionosphere (Sillito et al. 1989) (license: CC BY 4.0). The Simpson image, which is notably used in (Houdard et al. 2018) comes from the website page of Antoine Houdard (license: not specified).

## References

* Akhøj, M., Benn, J., Grong, E., Sommer, S., Pennec, X.: Principal subbundles for dimension reduction. arXiv (2023). https://doi.org/10.48550/arXiv.2307.03128

* Aeberhard, S., Forina, M.: Wine. UCI Machine Learning Repository (1991). https://doi.org/10.24432/C5PC7J

* Akaike, H.: A new look at the statistical model identification. IEEE Transactions on Automatic Control 19(6), 716-723 (1974) https://doi.org/10.1109/TAC.1974.1100705

* Bergé, L., Bouveyron, C., Girard, S.: HDclassif: An R Package for Model-Based Clustering and Discriminant Analysis of High-Dimensional
Data. Journal of Statistical Software **46**(6)
(2012) https://doi.org/10.18637/jss.v046.i06

Bergeron, F., Bousquet-Melou, M., Dulucq, S.: Standard Paths in the Composition Poset. Annales des sciences mathematiques du Quebec **19**(2), 139–151 (1995)

Bensmail, H., Celeux, G.: Regularized Gaussian Discriminant Analysis through Eigenvalue Decomposition. Journal of the American Statistical Association **91**(436), 1743–1748 (1996) https://doi.org/10.1080/01621459.1996.10476746

Bai, Z., Choi, K.P., Fujikoshi, Y.: Consistency of AIC and BIC in estimating the number of significant components in high-dimensional principal component analysis. The Annals of Statistics **46**(3) (2018) https://doi.org/10.1214/17-AOS1577

Biernacki, C., Celeux, G., Govaert, G.: Assessing a mixture model for clustering with the integrated completed likelihood. IEEE Transactions on Pattern Analysis and Machine Intelligence **22**(7), 719–725 (2000) https://doi.org/10.1109/34.865189

Bouveyron, C., Celeux, G., Girard, S.: Intrinsic dimension estimation by maximum likelihood in isotropic probabilistic PCA. Pattern Recognition Letters **32**(14), 1706–1713 (2011) https://doi.org/10.1016/j.patrec.2011.07.017

Buades, A., Coll, B., Morel, J.M.: A non-local algorithm for image denoising. In: 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05), vol. 2, pp. 60–65 (2005). https://doi.org/10.1109/CVPR.2005.38

Bouveyron, C., Celeux, G., Murphy, T.B., Raftery, A.E.: Model-Based Clustering and Classification for Data Science: With Applications in R, (2019). https://doi.org/10.1017/9781108644181

Bouveyron, C., Girard, S., Schmid, C.: High-dimensional data clustering. Computational Statistics & Data Analysis **52**(1), 502–519 (2007) https://doi.org/10.1016/j.csda.2007.02.009

Bouveyron, C., Girard, S., Schmid, C.: High-Dimensional Discriminant Analysis. Communications in Statistics - Theory and Methods **36**(14), 2607–2623 (2007) https://doi.org/10.1080/03610920701271095

Boser, B.E., Guyon, I.M., Vapnik, V.N.: A training algorithm for optimal margin classifiers. In: Proceedings of the Fifth Annual Workshop on Computational Learning Theory. COLT ’92, pp. 144–152 (1992). https://doi.org/10.1145/130385.130401

Bishop, C.: Bayesian PCA. In: NeurIPS, vol. 11 (1998)

Bach, F., Jenatton, R., Mairal, J., Obozinski, G.: Optimization with Sparsity-Inducing Penalties. Foundations and Trends® in Machine Learning **4**(1), 1–106 (2012) https://doi.org/10.1561/2200000015

Bouveyron, C., Latouche, P., Mattei, P.-A.: Exact dimensionality selection for Bayesian PCA. Scandinavian Journal of Statistics **47**(1), 196–211 (2020) https://doi.org/10.1111/sjos.12424

Buet, B., Pennec, X.: Flagfolds. arXiv (2023). https://doi.org/10.48550/arXiv.2305.10583

Banfield, J.D., Raftery, A.E.: Model-Based Gaussian and Non-Gaussian Clustering. Biometrics **49**(3), 803–821 (1993) https://doi.org/10.2307/2532201

Breiman, L.: Random Forests. Machine Learning **45**(1), 5–32 (2001) https://doi.org/10.1023/A:1010933404324

Byrski, K., Tabor, J., Spurek, P.: Gaussian model for closed curves. Expert Systems with Applications **249**, 123615 (2024) https://doi.org/10.1016/j.eswa.2024.123615

Cattell, R.B.: The Scree Test For The Number Of Factors. Multivariate Behavioral Research **1**(2), 245–276 (1966) https://doi.org/10.1207/s15327906mbr0102.10
Cappozzo, A., Casa, A., Fop, M.: Sparse Model-Based Clustering of Three-Way Data via Lasso-Type Penalties. Journal of Computational and Graphical Statistics, 1–21 (2025) https://doi.org/10.1080/10618600.2024.2429705

Celeux, G., Chretien, S., Forbes, F., Mkhadri, A.: A Component-Wise EM Algorithm for Mixtures. Journal of Computational and Graphical Statistics 10(4), 697–712 (2001) https://doi.org/10.1198/106186001317243403

Celeux, G., Govaert, G.: Gaussian parsimonious clustering models. Pattern Recognition 28(5), 781–793 (1995) https://doi.org/10.1016/0031-3203(94)00125-6

Chang, W.-C.: On Using Principal Components Before Separating a Mixture of Two Multivariate Normal Distributions. Journal of the Royal Statistical Society. Series C (Applied Statistics) 32(3), 267–275 (1983) https://doi.org/10.2307/2347949

Cavicchia, C., Vichi, M., Zaccaria, G.: Parsimonious ultrametric Gaussian mixture models. Statistics and Computing 34(3), 108 (2024) https://doi.org/10.1007/s11222-024-10405-9

Dempster, A.P., Laird, N.M., Rubin, D.B.: Maximum Likelihood from Incomplete Data via the EM Algorithm. Journal of the Royal Statistical Society. Series B (Methodological) 39(1), 1–38 (1977) https://doi.org/10.1111/j.2517-6161.1977.tb01600.x

Efros, A.A., Leung, T.K.: Texture synthesis by non-parametric sampling. In: Proceedings of the Seventh IEEE International Conference on Computer Vision, vol. 2, pp. 1033–1038 (1999). https://doi.org/10.1109/ICCV.1999.790383

Figueiredo, M.A.T., Jain, A.K.: Unsupervised learning of finite mixture models. IEEE Transactions on Pattern Analysis and Machine Intelligence 24(3), 381–396 (2002) https://doi.org/10.1109/34.990138

Fop, M., Murphy, T.B., Scrucca, L.: Model-based clustering with sparse covariance matrices. Statistics and Computing 29(4), 791–819 (2019) https://doi.org/10.1007/s11222-018-9838-y

Fraley, C., Raftery, A.E.: Model-Based Clustering, Discriminant Analysis, and Density Estimation. Journal of the American Statistical Association 97(458), 611–631 (2002) https://doi.org/10.1198/016214502760047131

Finder, S.E., Treister, E., Freifeld, O.: Effective Learning of a GMRF Mixture Model. IEEE Access 10, 7289–7299 (2022) https://doi.org/10.1109/ACCESS.2022.3141911

Ghahramani, Z., Beal, M.J.: Variational inference for Bayesian mixtures of factor analysers. In: NeurIPS, pp. 449–455 (1999)

Ghahramani, Z., Hinton, G.E.: The EM algorithm for mixtures of factor analyzers. Technical report, Technical Report CRG-TR-96-1, University of Toronto (1996)

Hubert, L., Arabie, P.: Comparing partitions. Journal of Classification 2(1), 193–218 (1985) https://doi.org/10.1007/BF01908075

Houdard, A., Bouveyron, C., Delon, J.: High-Dimensional Mixture Models for Unsupervised Image Denoising (HDMI). SIAM Journal on Imaging Sciences 11(4), 2815–2846 (2018) https://doi.org/10.1137/17M1135694

Heaps, S.E., Jermyn, I.H.: Structured prior distributions for the covariance matrix in latent factor models. Statistics and Computing 34(4), 143 (2024) https://doi.org/10.1007/s11222-024-10454-0

Halko, N., Martinsson, P.G., Tropp, J.A.: Finding Structure with Randomness: Probabilistic Algorithms for Constructing Approximate Matrix Decompositions. SIAM Review 53(2), 217–288 (2011) https://doi.org/10.1137/090771806

Harris, C.R., Millman, K.J., Walt, S.J., Gommers, R., Virtanen, P., Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N.J., Kern, R., Picus, M., Hoyer, S., Kerkwijk, M.H., Brett, M., Haldane, A., Río, J.F., Wiebe, M., Peterson, P., Gérard-Marchant, P., Sheppard, K., Reddy, T., Weckesser, W., Abbasi, H., Gohlke, C., Oliphant, T.E.: Array programming with NumPy. Nature 585(7825), 357–362 (2020) https://doi.org/10.1038/s41586-020-2649-2
Hershey, J.R., Olsen, P.A.: Approximating the Kullback Leibler Divergence Between Gaussian Mixture Models. In: 2007 IEEE International Conference on Acoustics, Speech and Signal Processing - ICASSP '07, vol. 4, pp. 317–320 (2007). https://doi.org/10.1109/ICASSP.2007.366913

Hosseini, R., Sra, S.: An alternative to EM for Gaussian mixture models: batch and stochastic Riemannian optimization. Mathematical Programming **181**(1), 187–223 (2020) https://doi.org/10.1007/s10107-019-01381-4

Hurvich, C.M., Tsai, C.-L.: Regression and time series model selection in small samples. Biometrika **76**(2), 297–307 (1989) https://doi.org/10.1093/biomet/76.2.297

Hunter, J.D.: Matplotlib: A 2D Graphics Environment. Computing in Science & Engineering **9**(3), 90–95 (2007) https://doi.org/10.1109/MCSE.2007.55

Josse, J., Husson, F.: Selecting the number of components in principal component analysis using cross-validation approximations. Computational Statistics & Data Analysis **56**(6), 1869–1879 (2012) https://doi.org/10.1016/j.csda.2011.11.012

Jiang, J., Nguyen, T., Rao, J.S.: The E-MS Algorithm: Model Selection With Incomplete Data. Journal of the American Statistical Association **110**(511), 1136–1147 (2015) https://doi.org/10.1080/01621459.2014.948545

Jolliffe, I.T.: Principal Component Analysis. Springer Series in Statistics, (2002). https://doi.org/10.1007/b98835

Johnstone, I.M., Paul, D.: PCA in High Dimensions: An Orientation. Proceedings of the IEEE **106**(8), 1277–1292 (2018) https://doi.org/10.1109/JPROC.2018.2846730

Keribin, C.: Consistent Estimation of the Order of Mixture Models. Sankhyā: The Indian Journal of Statistics, Series A (1961-2002) **62**(1), 49–66 (2000)

Lebrun, M., Colom, M., Buades, A., Morel,

J.M.: Secrets of image denoising cuisine. Acta Numerica **21**, 475–576 (2012) https://doi.org/10.1017/S0962492912000062

Lloyd, S.: Least squares quantization in PCM. IEEE Transactions on Information Theory **28**(2), 129–137 (1982) https://doi.org/10.1109/TIT.1982.1056489

Loconte, L., Sladek, A.M., Mengel, S., Trapp, M., Solin, A., Gillis, N., Vergari, A.: Subtractive Mixture Models via Squaring: Representation and Learning. In: The Twelfth International Conference on Learning Representations (2023)

Mairal, J., Bach, F., Ponce, J., Sapiro, G.: Online Learning for Matrix Factorization and Sparse Coding. Journal of Machine Learning Research **11**(2), 19–60 (2010)

McInnes, L., Healy, J., Astels, S.: hdbscan: Hierarchical density based clustering. Journal of Open Source Software **2**(11), 205 (2017) https://doi.org/10.21105/joss.00205

McNicholas, P.D., Murphy, T.B.: Parsimonious Gaussian mixture models. Statistics and Computing **18**(3), 285–296 (2008) https://doi.org/10.1007/s11222-008-9056-0

McLachlan, G.J., Peel, D.: Finite Mixture Models. Wiley Series in Probability and Statistics, (2000). https://doi.org/10.1002/0471721182

McLachlan, G.J., Peel, D., Bean, R.W.: Modelling high-dimensional data by mixtures of factor analyzers. Computational Statistics & Data Analysis **41**(3), 379–388 (2003) https://doi.org/10.1016/S0167-9473(02)00183-4

Malsiner-Walli, G., Frühwirth-Schnatter, S., Grün, B.: Model-based clustering based on sparse finite Gaussian mixtures. Statistics and Computing **26**(1), 303–324 (2016) https://doi.org/10.1007/s11222-014-9500-2

Ng, A., Jordan, M., Weiss, Y.: On Spectral Clustering: Analysis and an algorithm. In: NeurIPS, vol. 14 (2001)

Oudoumanessah, G., Coudert, T., Lartizien, C., Dojat, M., Christen, T., Forbes, F.: Scalable
magnetic resonance fingerprinting: Incremental
inference of high dimensional elliptical mixtures
from large data volumes (2024). https://hal-
science/hal-04836864

Olshausen, B.A., Field, D.J.: Emergence of
simple-cell receptive field properties by learn-
ing a sparse code for natural images. Nature
**381**(6583), 607–609 (1996) https://doi.org/10.
1038/381607a0

Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., Duchesnay, E.: Scikit-learn: Machine Learning in Python. Journal of Machine Learning Research **12**(85), 2825–2830 (2011)

Rand, W.M.: Objective Criteria for the Evaluation
of Clustering Methods. Journal of the Amer-
ican Statistical Association **66**(336), 846–850
(1971) https://doi.org/10.1080/01621459.1971.
10482356

Rousseau, J., Mengersen, K.: Asymptotic
Behaviour of the Posterior Distribution in
Overfitted Mixture Models. Journal of the
Royal Statistical Society Series B: Statistical
Methodology **73**(5), 689–710 (2011) https:
//doi.org/10.1111/j.1467-9868.2011.00781.x

Richardson, E., Weiss, Y.: On GANs and GMMs.
In: NeurIPS, vol. 31 (2018)

Schwarz, G.: Estimating the Dimension of
a Model. The Annals of Statistics **6**(2),
461–464 (1978) https://doi.org/10.1214/aos/
1176344136

Scrucca, L., Fop, M., Murphy, T.B., Raftery, A.E.: mclust 5: Clustering, Classification and Density Estimation Using Gaussian Finite Mixture Models. The R journal **8**(1), 289–317 (2016)

Szwagier, T., Olikier, G., Pennec, X.: Eigengap
Sparsity for Covariance Parsimony. In: Nielsen,
F., Barbaresco, F. (eds.) Geometric Science of
Information, pp. 50–59. Springer, Cham (2025).
https://doi.org/10.1007/978-3-032-03921-7_6

Szwagier, T., Pennec, X.: The curse of isotropy: from principal components to principal subspaces. arXiv (2025). https://doi.org/10.48550/arXiv.2307.15348

Sigillito, V.G., Wing, S.P., Hutton, L.V., Baker, K.B.: Ionosphere. UCI Machine Learning Repository (1989). https://doi.org/10.24432/C5W01B

Tipping, M.E., Bishop, C.M.: Mixtures of Probabilistic Principal Component Analyzers. Neural Computation **11**(2), 443–482 (1999) https://doi.org/10.1162/089976699300016728

Tipping, M.E., Bishop, C.M.: Probabilistic Principal Component Analysis. Journal of the Royal Statistical Society Series B: Statistical Methodology **61**(3), 611–622 (1999) https://doi.org/10.1111/1467-9868.00196

Tyler, D.E., Yi, M.: Lassoing eigenvalues.
Biometrika **107**(2), 397–414 (2020) https://doi.
org/10.1093/biomet/asz076

Virtanen, P., Gommers, R., Oliphant, T.E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., Walt, S.J., Brett, M., Wilson, J., Millman, K.J., Mayorov, N., Nelson, A.R.J., Jones, E., Kern, R., Larson, E., Carey, C.J., Polat, I., Feng, Y., Moore, E.W., VanderPlas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero, E.A., Harris, C.R., Archibald, A.M., Ribeiro, A.H., Pedregosa, F., Mulbregt, P.: SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature Methods **17**(3), 261–272 (2020) https://doi.org/10.1038/s41592-019-0686-2

Verbockhaven, M., Rudkiewicz, T., Chevallier, S., Charpiat, G.: Growing Tiny Networks: Spotting Expressivity Bottlenecks and Fixing Them Optimally. Transactions on Machine Learning Research (2024)

Ward, J.H.: Hierarchical Grouping to Optimize
an Objective Function. Journal of the Amer-
ican Statistical Association **58**(301), 236–244
(1963) https://doi.org/10.1080/01621459.1963.
10500845
Wolberg, W.H., Mangasarian, O.L., Street, W.N.: Breast Cancer Wisconsin (Diagnostic). UCI Machine Learning Repository (1995). https://doi.org/10.24432/C5DW2B

Walt, S.v.d., Schönberger, J.L., Nunez-Iglesias, J., Boulogne, F., Warner, J.D., Yager, N., Gouil-lart, E., Yu, T.: scikit-image: image processing in Python. PeerJ 2, 453 (2014) https://doi.org/10.7717/peerj.453

Ye, K., Wong, K.S.-W., Lim, L.-H.: Optimization on flag manifolds. Mathematical Programming 194(1), 621–660 (2022) https://doi.org/10.1007/s10107-021-01640-3

Zoran, D., Weiss, Y.: From learning models of natural image patches to whole image restoration. In: 2011 International Conference on Computer Vision, pp. 479–486 (2011). https://doi.org/10.1109/ICCV.2011.6126278