# Eigengap Sparsity for Covariance Parsimony

Tom Szwagier, Guillaume Olikier, Xavier Pennec

► To cite this version:

Tom Szwagier, Guillaume Olikier, Xavier Pennec. Eigengap Sparsity for Covariance Parsimony. GSI'25 - 7th International Conference on Geometric Science of Information, Oct 2025, Saint-Malo (35400), France. hal-05047196v2

HAL Id: hal-05047196

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

Submitted on 11 Jul 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
# Eigengap Sparsity for Covariance Parsimony

Tom Szwagier(✉), Guillaume Olikier(✉), and Xavier Pennec(✉)

Université Côte d'Azur, and Inria, Epione Project Team, Sophia-Antipolis, France
{tom.szwagier,guillaume.olikier,xavier.pennec}@inria.fr

**Abstract.** Covariance estimation is a central problem in statistics. An important issue is that there are rarely enough samples $n$ to accurately estimate the $p(p+1)/2$ coefficients in dimension $p$. Parsimonious covariance models are therefore preferred, but the discrete nature of model selection makes inference computationally challenging. In this paper, we propose a relaxation of covariance parsimony termed “eigengap sparsity” and motivated by the good accuracy–parsimony tradeoffs of eigenvalue-equalization in covariance matrices. This penalty can be included in a penalized-likelihood framework that we propose to solve with a projected gradient descent on a monotone cone. The algorithm turns out to resemble an isotonic regression of mutually-attracted sample eigenvalues, drawing an interesting link between covariance parsimony and shrinkage.

**Keywords:** Covariance estimation · Parsimony · Eigengaps · Flag manifolds · Monotone cone · Isotonic regression.

## 1 Introduction

The *principle of parsimony*, also known as Occam's razor (“The simplest explanation is usually the best one.”), is central in statistics. It implies that for two competing statistical models with similar likelihoods, one should choose the model with fewer parameters. It can be quantitatively achieved via model selection criteria such as the Bayesian Information Criterion (BIC) [25], or the Akaike Information Criterion (AIC) [1]. Since the seminal LASSO paper [31], it has been common to investigate parsimonious model estimation as a regularized optimization problem with a *sparsity-inducing penalty* [2]. Sparsity involves $l^0$-norms and leads to hard combinatorial problems, which is why it is classically *relaxed* with $l^1$-norms, enjoying convergence and sparsity guarantees [2].

The principle of parsimony finds a particular interest in the fundamental problem of *covariance estimation* [23]. Indeed, in dimension $p$, covariance matrices have $p(p+1)/2$ independent parameters, which makes them often *over-parameterized* with respect to the number $n$ of available samples. The *sample covariance matrix*—which is the maximum likelihood estimate of the population covariance matrix—then suffers from many *sampling errors*. A typical high-dimensional phenomenon is the Marčenko–Pastur law [20], which describes an important dispersion of the sample eigenvalues around their true value. This *eigenvalue-spreading* phenomenon [14] is illustrated in the seminal paper of
Ledoit and Wolf [17, Fig. 4] and motivates the idea of shrinking eigenvalues towards a mean value. To that extent, many regularized covariance estimation methods have emerged, notably involving shrinkage [26,17,9,18], sparsity [13,10], Bayesian methods [36,21] and random matrix theory [15].

In this paper, we propose to tackle the problem of covariance estimation with a penalized-likelihood optimization on positive-definite matrices:

$$ \hat{\Sigma} \in \underset{\Sigma \in \mathrm{Sym}_+(p)}{\operatorname{argmin}} \left(-2 \ln \mathcal{L}(\Sigma) + \alpha \dim(\Sigma)\right), \qquad (1) $$

where $\mathcal{L}$: $\mathrm{Sym}_+(p) \to \mathbb{R}$ is a likelihood function (e.g. Gaussian [34,29], elliptical [8] or Gaussian mixture [33,7,27]), $\alpha \in \mathbb{R}_{\ge 0}$ is a constant (like $\alpha_n = \ln n$ for the BIC, $\alpha = 2$ for the AIC or any tuning hyperparameter) and $\dim: \mathrm{Sym}_+(p) \to \mathbb{Z}_{>0}$ is a penalty on the number of covariance parameters, i.e. the dimension of the submanifold of $\mathrm{Sym}_+(p)$ to which $\Sigma$ belongs, as properly defined later. There are many ways to partition the space of covariance matrices into submanifolds (e.g. via coefficient-wise sparsity [10] or rank constraints [2]). The one that we defend in this paper is the stratification of covariance matrices by the multiplicities of the eigenvalues [11]. This original notion of parsimony is strongly motivated by the recent principal subspace analysis methodology [29] (equalizing close eigenvalues yields a quadratic decrease in the number of parameters without much decreasing the likelihood) which is a natural extension of the seminal probabilistic principal component analysis (PPCA) of Tipping and Bishop [34].

This paper is organized as follows. In section 2, we study the aforementioned stratification of the optimization space $\mathrm{Sym}_+(p)$, which enables to formally define the penalty dim in terms of eigenvalue multiplicities. We deduce that dim is constant on fixed-multiplicities submanifolds of $\mathrm{Sym}_+(p)$, which drastically simplifies the optimization problem. The bad news, however, is that solving the latter would require solving $2^{p-1}$ (unpenalized) optimization problems, which is computationally prohibitive, even in moderate dimensions. To that extent, we propose in section 3 a relaxation of the penalty. This $\ell^1$-relaxation is based on the intriguing result that the number of covariance parameters can be exactly rewritten as an $\ell^0$-norm of eigenvalue gaps. We will subsequently refer to our parsimony-inducing penalty as eigengap sparsity. We study the relaxed problem in the Gaussian setting in section 4 and propose a projected gradient descent algorithm on the monotone cone of eigenvalues to address it. The algorithm draws an interesting link between covariance parsimony and shrinkage methods [18], and it shows promising results in synthetic experiments in section 5.

## 2 Stratification of $\mathrm{Sym}_+(p)$ by Eigenvalue Multiplicities

Let $\Sigma \in \mathrm{Sym}_+(p)$ eigendecompose as $\Sigma = \sum_{j=1}^p \lambda_j v_j v_j^\top$ with $(\lambda_1, \dots, \lambda_p) \in \mathbb{R}^p$ the eigenvalues (ordered-decreasing and positive) and $(v_1, \dots, v_p) \in \mathbb{O}(p)$ some associated eigenvectors ($\mathbb{O}(p)$ denotes the orthogonal group). If all the eigenvalues are simple, then we need $p$ parameters to describe the eigenvalues and $\dim(\mathbb{O}(p)) = p(p-1)/2$ parameters to describe the eigenvectors, resulting in a
total of dim(Sym$_{+}$(p)) = p(p + 1)/2 covariance parameters to estimate. Otherwise, some eigenvalues are multiple, therefore we need fewer parameters not only for these eigenvalues, but also for the associated eigenspaces—since we quotient the rotational invariance of the eigenvectors within the eigenspace they span.

More precisely, let $C(p) := \{\gamma \in \mathbb{Z}_{>0}^d : \sum_{k=1}^d \gamma_k = p, d \in \llbracket 1, p \rrbracket\}$ be the set of compositions of $p$ [4]. Let $\gamma := (\gamma_1, \dots, \gamma_d) \in C(p)$, $q_k := \sum_{l=1}^k \gamma_l$, and

$$ \mathcal{K}_+(\gamma) := \{\lambda \in \mathbb{R}^p : \lambda_1 = \dots = \lambda_{q_1} > \lambda_{q_1+1} = \dots = \lambda_{q_2} > \dots > 0\} \quad (2) $$

be the set of ordered-decreasing piecewise-constant p-tuples of positive reals, whose pieces have respective sizes $\gamma_1, \gamma_2, \dots, \gamma_d$. By eigenvalue decomposition,

$$ \mathrm{Sym}_+(\gamma) := \{V \mathrm{diag}(\lambda) V^\top : \lambda \in \mathcal{K}_+(\gamma), V \in \mathrm{O}(p)\} \quad (3) $$

is the set of symmetric positive-definite matrices with eigenvalue multiplicities $\gamma_1, \dots, \gamma_d$. It is a submanifold of $\mathrm{Sym}_+(p)$ of dimension $d + (p^2 - \sum_{k=1}^d \gamma_k^2)/2$.

**Remark 1.** The theory behind this result is deep and has been worked out in [11]. Let us still give a concise justification. Let $\mathrm{Fl}(\gamma)$ be the flag manifold of type $\gamma$, i.e. the space of mutually-orthogonal linear subspaces $\mathcal{V}_1, \dots, \mathcal{V}_d$ of respective dimensions $\gamma_1, \dots, \gamma_d$ [37,28]. Let $\Pi_{\mathcal{V}_k}$ denote the orthogonal projection matrix onto $\mathcal{V}_k$. Then the eigendecomposition map $\lambda, \mathcal{V} \mapsto \sum_{k=1}^d \lambda_{q_k} \Pi_{\mathcal{V}_k}$ induces a diffeomorphism between $\mathcal{K}_+(\gamma) \times \mathrm{Fl}(\gamma)$ and $\mathrm{Sym}_+(\gamma)$. The submanifold and dimension properties then follow from those of $\mathcal{K}_+(\gamma)$ and $\mathrm{Fl}(\gamma)$.

The previous result tells us that $\mathrm{Sym}_+(p)$ can be stratified into submanifolds $\mathrm{Sym}_+(\gamma)$ of dimension $d + (p^2 - \sum_{k=1}^d \gamma_k^2)/2$, as $\mathrm{Sym}_+(p) = \bigsqcup_{\gamma \in C(p)} \mathrm{Sym}_+(\gamma)$. Therefore, for any $\Sigma \in \mathrm{Sym}_+(p)$, one can define $\mathrm{dim}(\Sigma)$ as the dimension of the unique stratum $\mathrm{Sym}_+(\gamma)$ it belongs to. This implies that the penalized covariance estimation problem (1) can be replaced with a set of unpenalized optimization problems on smooth matrix manifolds [37].

In the Gaussian case, problem (1) even has an explicit solution, referred to as principal subspace analysis [29].<fnref n="1" /> On each stratum $\mathrm{Sym}_+(\gamma)$, the maximum likelihood estimate $\hat{\Sigma}(\gamma)$ consists in the eigenvalue decomposition of the sample covariance matrix followed by a $\gamma$-block-averaging of the eigenvalues. Therefore, problem (1) is equivalent to

$$ \underset{\gamma \in C(p)}{\operatorname{argmin}} \left( -2 \ln \mathcal{L}(\hat{\Sigma}(\gamma)) + \alpha \left( d + \frac{p^2 - \sum_{k=1}^d \gamma_k^2}{2} \right) \right), \quad (4) $$

which can be solved exactly as $C(p)$ is finite. Unfortunately, $\#C(p) = 2^{p-1}$, so that solving exactly the problem is computationally prohibitive, even in moderate dimensions. And this is not to mention non-Gaussian distributions, which do not even have an explicit maximum likelihood estimate of the covariance.

<footnote n="1">¹ Principal subspace analysis was originally introduced to address a *curse of isotropy* in PCA, in which close sample eigenvalues yield unstable principal components.</footnote>
### 3 Relaxation via Eigengap Sparsity

We want to relax problem (1) in order to make it computationally tractable. The first step consists in rewriting the parsimony-inducing penalty dim as $ℓ^0$-norms of spectral gaps.

**Theorem 1.** Let $\gamma := (\gamma_1, \dots, \gamma_d) \in C(p)$ and $\Sigma \in \text{Sym}_+(\gamma)$ with eigenvalues $\lambda := (\lambda_1, \dots, \lambda_p) \in \mathcal{K}_+(\gamma)$. Let $\delta: \mathbb{R}^2 \to \mathbb{R}_{\ge 0}$ s.t. $\delta(a, b) = 0$ iff $a = b$. One has

$$ \dim(\text{Sym}_+(\gamma)) = 1 + \|\delta(\lambda_s, \lambda_{s+1})_{1\le s<p}\|_0 + \|\delta(\lambda_s, \lambda_t)_{1\le s<t\le p}\|_0. \quad (5) $$

*Proof.* By definition of $\lambda \in \mathcal{K}_+(\gamma)$, the eigenvalues can be divided into $d$ blocks of equal eigenvalues, of respective sizes $\gamma_1, \dots, \gamma_d$. Let $q_k := \sum_{l=1}^k \gamma_l$. Then one has $\delta(\lambda_s, \lambda_t) = 0 \iff \exists k \in [1, d]: (s, t) \in [q_{k-1}, q_k]^2$. Consequently, one has

$$ \|\delta(\lambda_s, \lambda_{s+1})_{1\le s<p}\|_0 := \#\{s: \delta(\lambda_s, \lambda_{s+1}) \neq 0\} = \#\{k: \lambda_{q_k} \neq \lambda_{q_{k+1}}\} = d-1. $$

Similarly, one has

$$ \|\delta(\lambda_s, \lambda_t)_{1\le s<t\le p}\|_0 = \#\{(s, t): s < t, \lambda_s \neq \lambda_t\} = \sum_{k=1}^{d} \gamma_k \sum_{l=k+1}^{d} \gamma_l. $$

The last quantity can be rewritten by double summation inversion as

$$ \sum_{1 \le k < l \le d} \gamma_k \gamma_l = \sum_{k=1}^{d} \gamma_k \left( p - \gamma_k - \sum_{l=1}^{k-1} \gamma_l \right) = p^2 - \sum_{k=1}^{d} \gamma_k^2 - \sum_{1 \le l < k \le d} \gamma_k \gamma_l. $$

Hence, $\|\delta(\lambda_s, \lambda_t)_{1\le s<t\le p}\|_0 = (p^2 - \sum_{k=1}^d \gamma_k^2)/2$, which concludes the proof. □

Although relatively simple, Theorem 1 is quite remarkable. While one may have thought of penalizing the successive gaps between eigenvalues (with total-variation-like penalties [32,12]), our theorem conveys a greater message: to account for the parsimony of the *eigenspaces* too, one should also penalize the gaps between non-adjacent eigenvalues. In other words, covariance parsimony tends to make all the eigenvalues attracted to each other, which is reminiscent of the recent findings of Ledoit and Wolf in shrinkage estimation [18]. One can now naturally relax problem (1) by replacing the $ℓ^0$-norm with the $ℓ^1$-norm.

Let $\mathcal{K}_+(p) = \{(\lambda_1, \dots, \lambda_p) \in \mathbb{R}^p: \lambda_1 \ge \dots \ge \lambda_p > 0\}$ denote the positive monotone cone. Then the $ℓ^1$-relaxation of problem (1) is

$$ \underset{V \in \mathcal{O}(p)}{\operatorname{argmin}} \left( -2 \ln \mathcal{L} (V \operatorname{diag}(\lambda) V^\top) + \alpha \sum_{s=1}^{p-1} \left( \delta(\lambda_s, \lambda_{s+1}) + \sum_{t=s+1}^{p} \delta(\lambda_s, \lambda_t) \right) \right). \quad (6) $$

The relaxed problem (6) can be solved with classical constrained optimization algorithms, for instance alternating between the update of the eigenvectors in $V \in \mathcal{O}(p)$ and the eigenvalues in $\lambda \in \mathcal{K}_+(p)$. In the remainder of the paper, we focus on the fundamental case where $\mathcal{L}$ is a Gaussian likelihood.
## 4 Projected Gradient Descent for Gaussian Densities

Under the Gaussian assumption, the maximum likelihood estimate of the eigen-vectors is explicit. The relaxed optimization problem (6) can then be simplified.

**Proposition 1.** Let $X := (x_1, \dots, x_n) \in \mathbb{R}^{p \times n}$. Let $S := \frac{1}{n}XX^\top$ be the sample covariance matrix and $\ell_1 \ge \dots \ge \ell_p$ be its eigenvalues. Then the relaxed problem (6) is equivalent to

$$ \underset{\lambda \in \mathcal{K}_+(p)}{\operatorname{argmin}} \left( \sum_{j=1}^{p} \left( \ln \lambda_j + \frac{\ell_j}{\lambda_j} \right) + \frac{\alpha}{n} \sum_{s=1}^{p-1} \left( \delta(\lambda_s, \lambda_{s+1}) + \sum_{t=s+1}^{p} \delta(\lambda_s, \lambda_t) \right) \right). \quad (7) $$

*Proof.* Under the Gaussian assumption, the log-likelihood can be rewritten as $-2 \ln \mathcal{L}(V \operatorname{diag}(\lambda) V^\top) = n(p \ln(2\pi) + \sum_{j=1}^p \ln \lambda_j + \operatorname{tr}(\operatorname{diag}(\lambda)^{-1} V^\top S V))$ [29]. An optimal $\hat{V} \in O(p)$ is any sequence of ordered eigenvectors of $S$ (with decreasing eigenvalues) [29]. This implies that $\operatorname{tr}(\operatorname{diag}(\lambda)^{-1} \hat{V}^\top S \hat{V}) = \sum_{j=1}^p \lambda_j^{-1} \ell_j$. We conclude the proof by dividing by $n$ and removing the constant $p \ln(2\pi)$. $\square$

**Remark 2.** In the remainder of this paper, to ensure the compactness of the constraint set and the good conditioning of the covariance matrices [6], the cone $\mathcal{K}_+(p)$ is replaced with $\mathcal{K}_\varepsilon(p) := \{\lambda \in \mathbb{R}^p: 1/\varepsilon \ge \lambda_1 \ge \dots \ge \lambda_p \ge \varepsilon\}$ with $\varepsilon := 10^{-10}$. It turns out to have little effect in practice, although it sometimes does in high dimensions or when some population eigenvalues are small.

We propose to solve problem (7) on $\mathcal{K}_\varepsilon(p)$ with a projected gradient descent (see e.g. [22] and references therein) summarized in Algorithm 1. A central observation is that the orthogonal projection onto the constraint set $\mathcal{K}_\varepsilon(p)$ corresponds to the well-known problem of isotonic regression (cf. [24] for a pedagogical explanation and illustration). This quadratic program can be exactly solved in linear time via the Pools Adjacent Violators Algorithm (PAVA) [5], which roughly boils down to sequentially investigating the eigenvalues and block-averaging the ones that violate the ordering constraints. It is followed by an eigenvalue thresholding to ensure that the eigenvalues fit the bounds in $\mathcal{K}_\varepsilon(p)$.

Consequently, at each iteration, Algorithm 1 somewhat makes a tradeoff between increasing the Gaussian likelihood and decreasing the eigengaps via

<caption>**Algorithm 1 Projected Gradient Descent for Gaussian Eigengap Sparsity**</caption>

**Input:** $\ell \in \mathcal{K}_0(p), Q \in O(p), \alpha \in \mathbb{R}_{\ge 0} \triangleright \text{sample eigenvalues/vectors, hyperparameter}$
$\lambda \leftarrow \ell \quad \triangleright \text{initialization with sample eigenvalues}$
**for** $i = 1, 2, \dots$ **do**
$g \leftarrow \frac{\lambda - \ell}{\lambda^2} + \frac{\alpha}{n} \partial_\lambda \sum_{s=1}^{p-1} (\delta(\lambda_s, \lambda_{s+1}) + \sum_{t=s+1}^{p} \delta(\lambda_s, \lambda_t)) \quad \triangleright \text{gradient (elementwise)}$
$\lambda \leftarrow \Pi_{\mathcal{K}_\varepsilon(p)}(\lambda - \beta g) \quad \triangleright \text{backtracking projected line search [22]}$
**end for**
**Output:** $\hat{\Sigma} := Q \operatorname{diag}(\lambda) Q^\top \quad \triangleright \text{optimal covariance matrix}$
a gradient step, and subsequently equalizes the eigenvalues that “intersect” (i.e. that cross the boundary of the monotone cone) via an isotonic regression. While isotonization has been historically employed by Stein as a post-shrinkage trick to avoid numerical issues caused by close eigenvalues (see e.g. [18]), it naturally arises in our Algorithm 1 as a projected gradient step. Therefore, while close eigenvalues are a *curse* in the classical shrinkage literature, they are a *blessing* in our framework since they get equalized and consequently bring parsimony.

The choice of eigengap function $\delta$ has a great impact on the results. While the most natural choice seems to be the *absolute eigengap* $\delta(a, b) := |a - b|$, we conjecture that it tends to make *all* the eigenvalues equal (in the Gaussian case), which is not always desired. Consequently, we choose to consider the *relative eigengap* $\delta(a, b) := (a - b)/a$ (for $a \ge b$), which is strongly justified by theory [29] and which also turns out to yield better experimental results.

# 5 Experiments

We run some synthetic experiments (with 10 repetitions) to compare our Eigen-gap Sparsity for Covariance Parsimony (ESCP, with $\alpha_n = \ln n$) to the original Principal Subspace Analysis (PSA) [29], which solves exactly problem (1) for a Gaussian density (cf. section 2). We also include for completeness the Sample Covariance Matrix (SCM) and the Ledoit-Wolf (LW) estimators [17], although they do not optimize the same objective function. Indeed, SCM optimizes the unpenalized-likelihood problem (1) (with $\alpha = 0$) and LW looks for a linear combination of the sample covariance and identity matrices that minimizes the expected Frobenius distance to the population covariance matrix.<fnref n="2" />

We sample $n$ points from a Gaussian distribution of covariance $\Sigma \in \text{Sym}_+(p)$ under three settings: (a) $n = 40, p = 20$, $\Sigma = I_{20}$, (b) $n = 200, p = 100$, $\Sigma = I_{100}$ and (c) $n = 400, p = 200$, $\Sigma = \text{diag}(10 I_{80}, 1 I_{80}, 0.1 I_{40})$, the last covariance matrix being denoted by $\Xi_{80,80,40}$ for short. The last setting is inspired from the baseline scenario of Ledoit and Wolf [18], which is claimed to be an interesting and difficult case. The first two settings have isotropic covariance matrices, therefore dimension 1, while the last setting has dimension $3+(200^2-80^2-80^2-40^2)/2 = 12803$. For each method, we report the penalized-likelihood $L_p(\hat{\Sigma}) := -2\ln \mathcal{L}(\hat{\Sigma}) + \alpha_n \text{dim}(\hat{\Sigma})$, the covariance estimation error $L_F(\hat{\Sigma}, \Sigma) := \|\hat{\Sigma} - \Sigma\|_F^2/p$ and the number of parameters $\text{dim}(\hat{\Sigma})$ (cf. Table 1). We also report in Fig. 1 the running times and the estimated scree plots.

There are many observations to make. First, SCM performs poorly in terms of covariance estimation ($L_F$), even for $p = 20$. This can be explained by the eigenvalue-spreading phenomenon [17], which we indeed observe in Fig. 1 (Middle)—where the sample eigenvalues are highly spread around their true value, which is 1. LW is much better in terms of covariance estimation ($L_F$, which is what it is designed for), although it may have a higher penalized-likelihood

<footnote n="2">² The code is written in Python 3.9 and run on an Intel® Core™ i7-10850H CPU with 16 GB of RAM: https://github.com/tomszwagier/eigengap-sparsity.</footnote>
<table><caption>Table 1: Evaluation of the penalized-likelihood $L_p(\hat{\Sigma})$, the covariance estimation error $L_F(\hat{\Sigma}, \Sigma)$ and the number of parameters $\dim(\hat{\Sigma})$ for different covariance estimation methods: Sample Covariance Matrix (SCM), Ledoit-Wolf (LW) [17], Principal Subspace Analysis (PSA) [29] and Eigengap Sparsity for Covariance Parsimony (ESCP) on multivariate Gaussian datasets with varying $(n, p, \Sigma)$.</caption><thead><tr><th colspan="4">(a) (40, 20, I<sub>20</sub>)</th><th colspan="4">(b) (200, 100, I<sub>100</sub>)</th><th colspan="4">(c) (400, 200, &Xi;<sub>80,80,40</sub>)</th></tr><tr><th>Method</th><th>L<sub>p</sub></th><th>L<sub>F</sub></th><th>dim</th><th>Method</th><th>L<sub>p</sub></th><th>L<sub>F</sub></th><th>dim</th><th>Method</th><th>L<sub>p</sub></th><th>L<sub>F</sub></th><th>dim</th></tr></thead><tbody><tr><td>SCM</td><td>33.2</td><td>0.544</td><td>210</td><td>SCM</td><td>202</td><td>0.50807</td><td>5050</td><td>SCM</td><td>532</td><td>9.9</td><td>20100</td></tr><tr><td>LW</td><td>31.3</td><td>0.008</td><td>126</td><td>LW</td><td>193</td><td>0.00030</td><td>3535</td><td>LW</td><td>661</td><td>6.7</td><td>20100</td></tr><tr><td>PSA</td><td>20.4</td><td>0.008</td><td>3</td><td>PSA</td><td>-</td><td>-</td><td>-</td><td>PSA</td><td>-</td><td>-</td><td>-</td></tr><tr><td>ESCP</td><td>20.5</td><td>0.004</td><td>3</td><td>ESCP</td><td>100</td><td>0.00004</td><td>1</td><td>ESCP</td><td>455</td><td>2.3</td><td>13508</td></tr></tbody></table>

$(L_p)$ since it shrinks the eigenvalues (so it reduces the likelihood with respect to SCM which is the maximum likelihood estimator) while not necessarily reducing the number of parameters. As expected, PSA yields the best penalized-likelihood $(L_p)$ since it finds its exact minimum. However, its running time is extremely large (cf. Fig. 1, Left), which makes it computationally prohibitive for $p \gtrsim 20$. In contrast, ESCP yields similar performance as PSA in drastically lower time, which makes it possible to run in higher dimensions. In dimensions 100 and 200 (Table 1 (b, c) and Fig. 1 (Middle, Right)) ESCP has a lower penalized-likelihood $(L_p)$ than both SCM and LW, which is not much surprising since it optimizes a relaxation of $L_p$. The good surprise is on the covariance estimation results: ESCP achieves a better estimation in average than LW although it is not designed to minimize $L_F$. This can be explained by looking more closely at Fig. 1. One can indeed see that ESCP almost perfectly recovers the true covariance matrix $I_{100}$ (b) and recovers well the piecewise-constant-eigenvalue matrix $\Xi_{80,80,40}$ (c), which we believe to be a challenging problem that cannot be prop-

<figure><img src="image_2.png" /><figcaption>Fig. 1: Evaluation of the covariance estimation methods for varying $(n, p, \Sigma)$. (Left) Running times in seconds, for increasing dimension $p$ in logarithmic scale. (Middle, Right) Scree plots of estimated eigenvalues in the settings (b) and (c).</figcaption></figure>
erly addressed by shrinkage methods yet.<fnref n="3" /> Putting this into perspective with the
curse of isotropy phenomenon [29], we expect the underlying principal subspaces
to be much more stable and interpretable than the sample eigenvectors output
by shrinkage methods. Hence, investigating covariance estimation through the
prism of covariance parsimony (induced by the equalization of close eigenvalues)
seems like a promising perspective.

# 6 Conclusion

We formulated the problem of parsimonious covariance estimation as a penalized-
likelihood optimization with an eigengap-sparsity-inducing penalty. The latter
was derived from the l¹-relaxation of intrinsic dimensionality in the stratifica-
tion of symmetric positive-definite matrices by the multiplicities of the eigenval-
ues. We proposed a projected gradient descent algorithm for the Gaussian case,
which boils down to iteratively updating the ordered-eigenvalues via a penalized-
likelihood gradient-step and equalizing the intersecting eigenvalues via an iso-
tonic regression which automatically induces parsimony. We illustrated the good
performance of our method in terms of speed (with respect to classical model
selection [29]), parsimony and covariance estimation (with respect to the sample
covariance matrix and shrinkage methods [17]).

The eigengap sparsity penalty is only at an early stage of research and many
perspectives arise: evaluating the performance of our method on more diverse
settings (eigenvalue profiles, non-Gaussian distributions) and choices of hyperpa-
rameters (regularization strength α and eigengap function δ); making the learn-
ing algorithm more robust and efficient in high dimensions using flag manifold
optimization [37,30]; studying the theoretical guarantees of the algorithm; con-
sidering Riemannian metrics on Sym_+(p) or K_+(p) beyond the embedded one
etc. One could also extend our penalized-likelihood estimation framework be-
yond covariance matrices, considering for instance graph Laplacians and Hes-
sians, whose parsimony is very important [19,16].

Our proposed approach draws some bridges with active areas of research.
First, the eigengap sparsity can be of interest to the sparse optimization com-
munity [2]. Our algorithm moreover surprisingly induces sparsity via an isotonic
regression instead of classical thresholding operators like in LASSO [31]. Second,
the eigengap penalty results from a relaxation of covariance parsimony, which
addresses the eigenvalue-spreading phenomenon studied in pioneering shrinkage
methods [17]. Interestingly, while most shrinkage methods focus on the eigen-
values, our framework for covariance estimation includes both the eigenvalues
and the flag of eigenspaces. Moreover, while low-rank models like PPCA [34]
only group the smallest eigenvalues (considered as noise), our models may also
group some of the largest eigenvalues, therefore better taking into account the
eigenvalue-spreading phenomenon. Third, our work has interesting links with the
elasso method from David E. Tyler and Mengxi Yi [35] and follow-up works [3].

<footnote n="3">³ As we can see, the average number of parameters (13508) is higher than the true one (12803), meaning that a few sample eigenvalues were not properly equalized.</footnote>
Although the motivations, the penalties and the optimization algorithms are quite different from ours, these methods also result in an automatic grouping of the eigenvalues. We believe that they could be seamlessly integrated in our framework by modifying the eigengap function $\delta$ and taking advantage of automatic differentiation.

**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.

**Disclosure of Interests.** The authors have no competing interests to declare that are relevant to the content of this article.

## References

1. Akaike, H.: A new look at the statistical model identification. IEEE Transactions on Automatic Control **19**(6), 716-723 (1974)

2. Bach, F., Jenatton, R., Mairal, J., Obozinski, G.: Optimization with sparsity-inducing penalties. Foundations and Trends® in Machine Learning **4**(1), 1-106 (2012)

3. Basiri, S., Ollila, E., Drašković, G., Pascal, F.: Fusing Eigenvalues. In: ICASSP 2019. pp. 4968-4972 (2019)

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

5. Best, M.J., Chakravarti, N.: Active set algorithms for isotonic regression; A unifying framework. Mathematical Programming **47**(1), 425-439 (1990)

6. Bickel, P.J., Levina, E.: Regularized estimation of large covariance matrices. The Annals of Statistics **36**(1), 199-227 (2008)

7. Bouveyron, C., Girard, S., Schmid, C.: High-dimensional data clustering. Computational Statistics & Data Analysis **52**(1), 502-519 (2007)

8. Chen, Y., Wiesel, A., Hero, A.O.: Robust Shrinkage Estimation of High-Dimensional Covariance Matrices. IEEE TSP **59**(9), 4097-4107 (2011)

9. Donoho, D., Gavish, M., Johnstone, I.: Optimal shrinkage of eigenvalues in the spiked covariance model. The Annals of Statistics **46**(4), 1742-1778 (2018)

10. Friedman, J., Hastie, T., Tibshirani, R.: Sparse inverse covariance estimation with the graphical lasso. Biostatistics **9**(3), 432-441 (2008)

11. Groisser, D., Jung, S., Schwartzman, A.: Geometric foundations for scaling-rotation statistics on symmetric positive definite matrices: Minimal smooth scaling-rotation curves in low dimensions. Electronic Journal of Statistics **11**(1), 1092-1159 (2017)

12. Harchaoui, Z., Lévy-Leduc, C.: Multiple Change-Point Estimation With a Total Variation Penalty. Journal of the American Statistical Association **105**(492), 1480-1493 (2010)

13. Huang, J.Z., Liu, N., Pourahmadi, M., Liu, L.: Covariance matrix selection and estimation via penalised normal likelihood. Biometrika **93**(1), 85-98 (2006)

14. Johnstone, I.M., Paul, D.: PCA in High Dimensions: An Orientation. Proceedings of the IEEE **106**(8), 1277-1292 (2018)

15. Karoui, N.E.: Spectrum estimation for large dimensional covariance matrices using random matrix theory. The Annals of Statistics **36**(6), 2757-2790 (2008)
16. Lambert, M., Bonnabel, S., Bach, F.: The limited-memory recursive variational Gaussian approximation (L-RVGA). *Statistics and Computing* **33**(70) (2023)

17. Ledoit, O., Wolf, M.: A well-conditioned estimator for large-dimensional covariance matrices. *Journal of Multivariate Analysis* **88**(2), 365–411 (2004)

18. Ledoit, O., Wolf, M.: Quadratic shrinkage for large covariance matrices. *Bernoulli* **28**(3), 1519–1547 (2022)

19. Lombaert, H., Grady, L., Polimeni, J.R., Cheriet, F.: FOCUSR: feature oriented correspondence using spectral regularization–a method for precise surface matching. *IEEE TPAMI* **35**(9), 2143–2160 (2013)

20. Marčenko, V.A., Pastur, L.A.: Distribution of eigenvalues for some sets of random matrices. *Mathematics of the USSR-Sbornik* **1**(4), 457–483 (1967)

21. Minka, T.: Automatic Choice of Dimensionality for PCA. In: NeurIPS. vol. 13. MIT Press (2000)

22. Olikier, G., Waldspurger, I.: Projected gradient descent accumulates at Bouligand stationary points. *SIAM Journal on Optimization* **35**(2), 1004–1029 (2025)

23. Pourahmadi, M.: Covariance Estimation: The GLM and Regularization Perspectives. *Statistical Science* **26**(3), 369–387 (2011)

24. Salmon, J.: Isotonic regression, https://josephsalmon.eu/blog/isotonic/

25. Schwarz, G.: Estimating the Dimension of a Model. The Annals of Statistics **6**(2), 461–464 (1978)

26. Stein, C.M.: Estimation of a covariance matrix (1975), 39th Annual Meeting IMS

27. Szwagier, T., Mattei, P.A., Bouveyron, C., Pennec, X.: Parsimonious Gaussian mixture models with piecewise-constant eigenvalue profiles (2025). https://doi.org/10.48550/arXiv.2507.01542

28. Szwagier, T., Pennec, X.: Rethinking the Riemannian Logarithm on Flag Manifolds as an Orthogonal Alignment Problem. In: Geometric Science of Information. pp. 375–383. Springer (2023)

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

30. Szwagier, T., Pennec, X.: Nested subspace learning with flags (2025). https://doi.org/10.48550/arXiv.2502.06022

31. Tibshirani, R.: Regression Shrinkage and Selection via the Lasso. Journal of the Royal Statistical Society. Series B (Methodological) **58**(1), 267–288 (1996)

32. Tibshirani, R., Saunders, M., Rosset, S., Zhu, J., Knight, K.: Sparsity and Smoothness Via the Fused Lasso. Journal of the Royal Statistical Society Series B: Statistical Methodology **67**(1), 91–108 (2005)

33. Tipping, M.E., Bishop, C.M.: Mixtures of Probabilistic Principal Component Analyzers. *Neural Computation* **11**(2), 443–482 (1999)

34. 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)

35. Tyler, D.E., Yi, M.: Lassoing eigenvalues. *Biometrika* **107**(2), 397–414 (2020)

36. Yang, R., Berger, J.O.: Estimation of a Covariance Matrix Using the Reference Prior. *The Annals of Statistics* **22**(3), 1195–1211 (1994)

37. Ye, K., Wong, K.S.W., Lim, L.H.: Optimization on flag manifolds. *Mathematical Programming* **194**(1), 621–660 (2022)