# WHITENING SPHERICAL GAUSSIAN MIXTURES IN THE LARGE-DIMENSIONAL REGIME

Mohammed Racim Moussa Boudjemaa, Alper Kalle, Xiaoyi Mai, José Henrique de Morais Goulart, Cédric Févotte

► To cite this version:

Mohammed Racim Moussa Boudjemaa, Alper Kalle, Xiaoyi Mai, José Henrique de Morais Goulart, Cédric Févotte. WHITENING SPHERICAL GAUSSIAN MIXTURES IN THE LARGE-DIMENSIONAL REGIME. 2025. hal-05349494

HAL Id: hal-05349494

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

Preprint submitted on 5 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
# WHITENING SPHERICAL GAUSSIAN MIXTURES IN THE LARGE-DIMENSIONAL REGIME

M. R. M. Boudjemaa,¹ A. Kalle,² X. Mai,³ J. H. de M. Goulart,¹ and C. Févotte¹

¹ IRIT, CNRS, Toulouse INP, Université de Toulouse

² List, CEA, Université Paris-Saclay

³ IMT, Université de Toulouse - Jean Jaurès.

## ABSTRACT

Whitening is a classical technique in unsupervised learning that can facilitate estimation tasks by standardizing data. An important application is the estimation of latent variable models via the decomposition of tensors built from high-order moments. In particular, whitening orthogonalizes the means of a spherical Gaussian mixture model (GMM), thereby making the corresponding moment tensor orthogonally decomposable, hence easier to decompose. However, in the large-dimensional regime (LDR) where data are high-dimensional and scarce, the standard whitening matrix built from the sample covariance becomes ineffective because the latter is spectrally distorted. Consequently, whitened means of a spherical GMM are no longer orthogonal. Using random matrix theory, we derive exact limits for their dot products, which are generally nonzero in the LDR. As our main contribution, we then construct a corrected whitening matrix that restores asymptotic orthogonality, allowing for performance gains in spherical GMM estimation.

**Index Terms—** Whitening, Gaussian mixture models, random matrix theory, spiked models, large-dimensional regime.

## 1 Introduction

Whitening is a classical technique in unsupervised learning and latent variable model estimation [1], [2]. Concretely, it exploits the latent geometry of the data by constructing a linear transformation that eliminates pairwise correlations and normalizes variances. This elementary yet powerful idea leads to a significant simplification of certain estimation tasks. For instance, whitening reduces independent component analysis (ICA) to the estimation of an orthogonal rotation of the whitened signals [2], [3]. Recent estimators of latent variable models based on decomposing tensors built from high-order moments also rely on whitening, as an intermediate step that makes tensor decomposition solvable by spectral methods [4], [5].

This strategy applies, in particular, to the estimation of spherical Gaussian mixture models (GMM) [5], [6], [7]. For concreteness, consider a spherical Gaussian mixture model (GMM) with K components, weights {$\omega_k$}$_{k=1}^K$ (positive, summing to one), means {$\mu_k$}$_{k=1}^K \subset \mathbb{R}^P$, and a common variance $\sigma^2 > 0$, so that an observation $\mathbf{x}_n$ is distributed as

$$ \mathbf{x}_n \sim \sum_{k=1}^{K} \omega_k \mathcal{N}(\mu_k, \sigma^2 \mathbf{I}_P). \quad (1) $$

As usual, one seeks to estimate the parameters $\sigma^2$ and {$\langle \omega_k, \mu_k \rangle$}$_{k=1}^K$

from N i.i.d. samples $\mathbf{x}_n$. The uncentered covariance matrix reads

$$ \Sigma := \mathbb{E}[\mathbf{x}_n \mathbf{x}_n^\top] = \sum_{k=1}^{K} \omega_k \mu_k \mu_k^\top + \sigma^2 \mathbf{I}_P = \mathbf{M}_2 + \sigma^2 \mathbf{I}_P, \quad (2) $$

where $\mathbf{M}_2 = \sum_{k=1}^K \omega_k \mu_k \mu_k^\top$ encodes the “inter-cluster” covariance. Since $\mathbf{M}_2$ carries information on the mixture structure, it makes sense to whiten the data with respect to this covariance, which amounts to constructing from $\Sigma$ a $K \times P$ whitening matrix $\mathbf{W}$ such that $\mathbf{W}\mathbf{M}_2\mathbf{W}^\top = \mathbf{I}_K$. By construction, such a linear map performs two things at once: (i) dimensionality reduction, since it maps data onto the K-dimensional “signal” subspace; (ii) orthogonalization of the whitened means $\mathbf{W}\mu_k$ (see Section 3 for more details).

Now, Hsu and Kakade showed in [5] that one can exploit third-order moments of $\mathbf{x}_n$ to build an order-3 tensor $\mathcal{M}_3$ having a structure analogous to $\mathbf{M}_2$, that is, $\mathcal{M}_3 = \sum_{k=1}^K \omega_k \mu_k^{\otimes 3}$, where $\mu_k^{\otimes 3} := \mu_k \otimes \mu_k \otimes \mu_k$. This fact alone is relevant for estimation, since this decomposition of $\mathcal{M}_3$ into rank-one tensors is notoriously unique under mild constraints, unlike that of $\mathbf{M}_2$. Yet, computing the above decomposition is generally challenging, and this is where whitening helps: after whitening the data, one obtains a corresponding (transformed) tensor which is orthogonally decomposable, namely

$$ \mathcal{M}_3(\mathbf{W}, \mathbf{W}, \mathbf{W}) := \sum_{k=1}^{K} \omega_k (\mathbf{W} \mu_k)^{\otimes 3} = \sum_{k=1}^{K} \omega_k^{-1/2} \mathbf{v}_k^{\otimes 3}, \quad (3) $$

with {$\mathbf{v}_k := \omega_k^{1/2} \mathbf{W} \mu_k$}$_{k=1}^K$ orthonormal (as proven in Section 3). Hence, the estimation problem reduces to an orthogonal tensor decomposition on $\mathbb{R}^K$, and thus is solvable by stable spectral algorithms; see [4], [5]. This whitening-plus-tensor-decomposition procedure can serve as a stand-alone estimator or as an initialization to a maximum likelihood estimation algorithm such as EM.

However, in practice the whitening matrix is built from the sample covariance matrix $\hat{\Sigma} := \frac{1}{N} \mathbf{x}_n \mathbf{x}_n^\top$, which is not always an accurate estimate of $\Sigma$ [8]. This happens in particular in many modern applications where data are high-dimensional and (relatively) scarce, a situation which is formalized by the so-called large-dimensional regime (LDR), where both P (data dimension) and N (sample size) tend to infinity, with fixed $P/N = c \in (0, \infty)$. In this (common) regime, random matrix theory (RMT) models $\hat{\Sigma}$ as a low-rank “signal” perturbed by isotropic noise (a spiked model) and predicts that it is spectrally distorted: its eigenvalues follow a Marchenko-Pastur law and its leading eigenvectors are rotated away from the true ones [9], [10], [11], [12], [13]. As a consequence, the whitening matrix $\mathbf{W}$ built from $\hat{\Sigma}$ fails to orthogonalize the GMM component means $\mu_k$, and the whitened tensor moment $\mathcal{M}_3$ loses its orthogonal structure, which in turn degrades tensor-based GMM estimators.

Our goal in this paper is twofold: (i) precisely quantifying the impact of these distortions on the whitening step employed in these

<footnote>This work is supported by the LabEx CIMI (ANR-11-LABX-0040) and the AI cluster ANITI (ANR-23-IACL-0002). Part of this work was conducted while A. Kalle was an MSc intern student at IRIT.</footnote>
estimators, when operating in the LDR; (ii) proposing a corrected whitening based on RMT results for mitigating this impact. Specifically, we show that standard whitening no longer orthogonalizes the GMM clusters' means $\mu_k$ in the LDR, by giving exact closed-form limits of the dot products between whitened means. We then construct an RMT-corrected whitening operator that reweighs empirical eigen-pairs to restore orthogonality of the transformed means. Finally, as an illustration, we plug the corrected whitening into the LEARNGMM algorithm of [5], one of the simplest of tensor-based spherical GMM estimators. Our numerical results clearly show the degradation of its performance in the LDR, and the improvement brought by our corrected whitening.

## 2 Model and assumptions

We will henceforth consider the spherical GMM specified in (1).

**Assumptions.** We work under the following regime:

(A1) **Large-dimensional regime (LDR).** The sample number $N$ and the dimension $P$ grow large at the same rate with fixed $P/N = c \in (0, \infty)$.

To avoid trivial cases where estimation becomes either too easy or impossible in this regime, we also impose the following assumptions:

(A2) **Non-degeneracy of signal subspace.** We assume that rank($M_2$) = $K$ for $M_2$ the low rank term in (2).

(A3) **Non-vanishing, non-degenerate components.** The variance $\sigma^2$, and the prior probabilities $\{\omega_k\}_{k=1}^K$ are all nonzero and constant with respect to $P$ and $N$, for fixed $K$.

(A4) **Constant signal strength.** The norms of the GMM means $\mu_k$ and their alignments $\langle \mu_i, \mu_j \rangle$ are assumed to be constant with respect to $N$ and $P$.

**Phase transition of sample covariance in the LDR.** Consider the following spectral decompositions for the covariance $\Sigma$ given in (2) and its empirical version $\hat{\Sigma} = \frac{1}{N}x_n x_n^\top$:

$$ \Sigma = \sum_{k=1}^{P} \lambda_k \mathbf{u}_k \mathbf{u}_k^\top, \quad \hat{\Sigma} = \sum_{k=1}^{P} \hat{\lambda}_k \hat{\mathbf{u}}_k \hat{\mathbf{u}}_k^\top, \quad (4) $$

where $\{ (\lambda_k, \mathbf{u}_k) \}_{k=1}^P$, $\{ (\hat{\lambda}_k, \hat{\mathbf{u}}_k) \}_{k=1}^P$ are respectively the eigenpairs of $\Sigma$, $\hat{\Sigma}$ with $\lambda_1 \ge \dots \ge \lambda_P$, $\hat{\lambda}_1 \ge \dots \ge \hat{\lambda}_P$. Using (2), we get

$$ \Sigma = \sum_{k=1}^{K} \sigma^2 \ell_k \mathbf{u}_k \mathbf{u}_k^\top + \sigma^2 \mathbf{I}_P, \quad \text{with } \ell_k = (\lambda_k - \sigma^2)/\sigma^2. \quad (5) $$

Under (A4), well-known results from RMT indicate a phase transition occurs at $\sqrt{c}$ for $\ell_k$: spikes with $\ell_k > \sqrt{c}$ generate outlier eigenvalues $\hat{\lambda}_k$ and nonzero asymptotic eigenvector alignments between $\mathbf{u}_k$ and $\hat{\mathbf{u}}_k$; spikes with $\ell_k \le \sqrt{c}$ are absorbed by the Marchenko-Pastur distribution bulk and the alignments $\langle \mathbf{u}_k, \hat{\mathbf{u}}_k \rangle$ vanish [14], [10].

**Theorem 1 (Spiked-covariance limits in the LDR; adapted from [14], [10]).** Consider the Gaussian mixture model (1), and the uncentered covariance $\Sigma$ defined in (2) with its eigenpairs $\{ (\lambda_k, \mathbf{u}_k) \}_{k=1}^P$ and those $\{ (\hat{\lambda}_k, \hat{\mathbf{u}}_k) \}_{k=1}^P$ of its empirical version $\hat{\Sigma}$ as specified in (4), and $\ell_k$ as introduced in (5). Then, under (A1)-(A4),

$$ \hat{\lambda}_k \xrightarrow{a.s.} \hat{\lambda}_k := 1_{\{\ell_k > \sqrt{c}\}} \beta_k \lambda_k + 1_{\{\ell_k \le \sqrt{c}\}} (1 + \sqrt{c})^2, \quad (6) $$

$$ \mathbf{a}^\top \hat{\mathbf{u}}_k \hat{\mathbf{u}}_k^\top \mathbf{b} \xrightarrow{a.s.} 1_{\{\ell_k > \sqrt{c}\}} \psi_k \mathbf{a}^\top \mathbf{u}_k \mathbf{u}_k^\top \mathbf{b}, \quad \forall \mathbf{a}, \mathbf{b} \in \mathbb{R}^P, \quad (7) $$

where $\beta_k = 1 + \frac{c}{\ell_k}$ and $\psi_k = 1 - \frac{\beta_k^{-1} \ell_k}{\beta_k}$.

**Remark 1.** Directly applying (7) to $\mathbf{a} = \mathbf{b} = \mathbf{u}_k$ gives the asymptotic (squared) alignment between $\hat{\mathbf{u}}_k$ and $\mathbf{u}_k$:

$$ \mathbf{u}_k^\top \hat{\mathbf{u}}_k \hat{\mathbf{u}}_k^\top \mathbf{u}_k = (\hat{\mathbf{u}}_k^\top \mathbf{u}_k)^2 \xrightarrow{a.s.} 1_{\{\ell_k > \sqrt{c}\}} \psi_k =: \zeta_k. $$

**Remark 2.** When $c > 0$ and $\ell_k > \sqrt{c}$, one still has $\psi_k < 1$, so empirical eigenvectors do not perfectly align with their population counterparts. In addition, sample eigenvalues are distorted: $\hat{\lambda}_k \xrightarrow{a.s.} \hat{\lambda}_k = \beta_k \lambda_k$ for recoverable spikes and $\hat{\lambda}_k \xrightarrow{a.s.} \sigma^2(1 + \sqrt{c})^2$ when $\ell_k \le \sqrt{c}$. Both these effects are precisely why whitening becomes ineffective in the LDR.

## 3 Whitening in large dimension

**Population whitening matrix.** According to (4), the low-rank term $M_2$ of $\Sigma$ introduced in (2) admits the spectral decomposition $M_2 = U_K \Gamma U_K^\top$ where $U_K := [\mathbf{u}_1, \dots, \mathbf{u}_K] \in \mathbb{R}^{P \times K}$ and $\Gamma = \text{Diag}(\lambda_1 - \sigma^2, \dots, \lambda_K - \sigma^2)$. Assumption (A1) implies that matrix $\Gamma \in \mathbb{R}^{K \times K}$ is non-singular. The population whitening matrix with respect to the inter-cluster covariance $M_2$ is properly defined as $\mathbf{W} := \Gamma^{-1/2} \mathbf{U}_K^\top \in \mathbb{R}^{K \times P}$. Indeed, with this definition,

$$ \mathbf{W} M_2 \mathbf{W}^\top = \Gamma^{-1/2} \mathbf{U}_K^\top (U_K \Gamma \mathbf{U}_K^\top) U_K \Gamma^{-1/2} = \mathbf{I}_K. $$

Since $M_2 = \sum_{k=1}^K \omega_k \mu_k \mu_k^\top$, then $\sum_{k=1}^K \omega_k (\mathbf{W}\mu_k)(\mathbf{W}\mu_k)^\top = \mathbf{I}_K$. Using the rescaled and whitened means $\mathbf{v}_k = \omega_k^{1/2} \mathbf{W}\mu_k \in \mathbb{R}^K$ introduced in (3), the previous identity rewrites as $\sum_{k=1}^K \mathbf{v}_k \mathbf{v}_k^\top = \mathbf{I}_K$, which means that the family $\{\mathbf{v}_k\}_{k=1}^K$ is orthonormal. It follows that $\|\mathbf{W}\mu_k\|^2 = \omega_k^{-1}$ and $\langle \mathbf{W}\mu_i, \mathbf{W}\mu_j \rangle = 0$ for $i \neq j$.

**Whitening matrix estimation.** Let $\hat{M}_2 = \hat{\Sigma} - \hat{\sigma}^2 \mathbf{I}_P$ (for some chosen variance estimator $\hat{\sigma}^2$). Upon computing the spectral decomposition of $\hat{M}_2$ and keeping its $K$ dominant eigenpairs, $\hat{M}_2 \approx \hat{U}_K \hat{\Gamma} \hat{U}_K^\top$, one can estimate the whitening matrix as $\hat{\mathbf{W}} := \hat{\Gamma}^{-1/2} \hat{U}_K^\top \in \mathbb{R}^{K \times P}$, where the subscript $K$ indicates that only the $K$ leading eigenpairs are retained. In the classical regime ($P$ constant and $N \to \infty$), $\hat{\mathbf{W}} \xrightarrow{a.s.} \mathbf{W}$ and the orthogonality between whitened means is (asymptotically) verified. Yet, this property breaks down in the LDR, as seen next.

### 3.1 Residual alignment and phase transition

We now examine the loss of orthogonality between whitened means that arises in the LDR, as a consequence of Theorem 1. Given the above estimate of the whitening matrix $\hat{\mathbf{W}}$, the main quantity of interest is the empirical residual alignment between component means defined below:

$$ \rho_{i,j} := \frac{|\langle \hat{\mathbf{W}}\mu_i, \hat{\mathbf{W}}\mu_j \rangle|}{\|\hat{\mathbf{W}}\mu_i\| \|\hat{\mathbf{W}}\mu_j\|}, \quad (8) $$

for $i \neq j \in \{1, \dots, K\}$.

Since $\hat{\mathbf{W}}$ is built from a reduced spectral decomposition of $\hat{M}_2 = \hat{\Sigma} - \hat{\sigma}^2 \mathbf{I}_P$, it depends on the chosen variance estimator $\hat{\sigma}^2$. In the classical regime, simply taking $\hat{\sigma}^2 = \hat{\lambda}_p$ with any $p > K$ yields a consistent estimator of $\sigma^2$. Yet, this no longer holds in the LDR where $\hat{\lambda}_p$ with $p > K$ are distributed according to the Marchenko-Pastur law on the support $[\sigma^2(1 - \sqrt{c})^2, \sigma^2(1 + \sqrt{c})^2]$ [9]. Nevertheless, using Assumptions (A1) and (A4) one can show that the estimator $\hat{\sigma}^2 =$
$\frac{1}{NP} \sum_{n=1}^{N} \|\mathbf{x}_n\|^2$ is consistent, even when the data are not centered. Hence, in the following we consider that $\sigma^2$ is consistently estimated, focusing only on the effect of the spectral distortion described by Theorem 1.

**Lemma 1 (Limit of residual alignment).** Assume (A1)-(A4) and consider the same notation as in Theorem 1. Set $d_k := \mathbb{1}_{\ell_k \ge \sqrt{c}} \psi_k / (\lambda_k - \sigma^2) = \mathbb{1}_{\ell_k \ge \sqrt{c}} \psi_k / (\beta_k \lambda_k - \sigma^2)$, $\mathbf{D} := \text{Diag}(d_1, \dots, d_K)$ and $\mathbf{T} := \mathbf{D}^{1/2} \mathbf{U}_K^\top$. Then, we have

$$ \langle \hat{\mathbf{W}} \mu_i, \hat{\mathbf{W}} \mu_j \rangle \xrightarrow{a.s.} \langle \mathbf{T} \mu_i, \mathbf{T} \mu_j \rangle, \quad (9) $$

and, whenever both $\langle \mathbf{T} \mu_i, \mathbf{T} \mu_i \rangle > 0$ and $\langle \mathbf{T} \mu_j, \mathbf{T} \mu_j \rangle > 0$,

$$ \rho_{i,j} \xrightarrow{a.s.} \rho_{i,j}^\infty := \frac{|\langle \mathbf{T} \mu_i, \mathbf{T} \mu_j \rangle|}{\sqrt{\langle \mathbf{T} \mu_i, \mathbf{T} \mu_i \rangle} \sqrt{\langle \mathbf{T} \mu_j, \mathbf{T} \mu_j \rangle}}. \quad (10) $$

*Proof.* Notice first that

$$ \langle \hat{\mathbf{W}} \mu_i, \hat{\mathbf{W}} \mu_j \rangle = \mu_i^\top \hat{\mathbf{U}}_K \hat{\Gamma}^{-1} \hat{\mathbf{U}}_K^\top \mu_j = \sum_{k=1}^K \frac{(\mu_i^\top \hat{\mathbf{u}}_k)(\hat{\mathbf{u}}_k^\top \mu_j)}{\hat{\lambda}_k - \hat{\sigma}^2}. $$

The asymptotic limit of $(\mu_i^\top \hat{\mathbf{u}}_k)(\hat{\mathbf{u}}_k^\top \mu_j)$ is given by letting $a = \mu_i$, $b = \mu_j$ in (7). Retrieving also the limiting value of $\hat{\lambda}_k$ from (6) and recalling that $\hat{\sigma}^2 \xrightarrow{a.s.} \sigma^2$, we get

$$ \langle \hat{\mathbf{W}} \mu_i, \hat{\mathbf{W}} \mu_j \rangle \xrightarrow{a.s.} \sum_{k=1}^{K} \frac{\mathbb{1}_{\ell_k \ge \sqrt{c}} \psi_k (\mu_i^\top \mathbf{u}_k) (\mathbf{u}_k^\top \mu_j)}{\hat{\lambda}_k^2 - \sigma^2} \\ = \sum_{k=1}^{K} d_k (\mu_i^\top \mathbf{u}_k) (\mathbf{u}_k^\top \mu_j) = \langle \mathbf{T} \mu_i, \mathbf{T} \mu_j \rangle, $$

thus recovering (9). □

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

**Remark 3 (Classical regime).** Letting $c \to 0$, we get $\psi_k \to 1$ and $\hat{\lambda}_k \to \lambda_k$. Hence $\mathbf{T} \to \mathbf{W}$ in the classical regime, implying $\rho_{i,j} \to 0$ for $i \neq j$.

**Example 1.** Consider the case $K = 2$ with $\omega_1 = \omega_2 = 1/2$, and suppose that $\mu_1, \mu_2$ have unit norm and alignment $\langle \mu_1, \mu_2 \rangle = \cos \theta \in [0, 1)$. With some algebraic calculation, we get $\lambda_1 = \cos^2(\theta/2) + \sigma^2$, $\lambda_2 = \sin^2(\theta/2) + \sigma^2$, from which it follows that $\ell_1 = \cos^2(\theta/2)/\sigma^2$, $\ell_2 = \sin^2(\theta/2)/\sigma^2$. When $\ell_1, \ell_2 > \sqrt{c}$, applying Lemma 1 leads to

$$ \rho_{1,2} \xrightarrow{a.s.} \frac{q_1 - q_2}{q_1 + q_2}, \quad \text{where} \quad q_k = \frac{\ell_k^4 - c\ell_k^2}{\ell_k(\ell_k + c)(\ell_k^2 + c\ell_k + c)}. $$

It can be seen from the above equation that the alignment $\rho_{1,2}$ between the whitened means $\hat{\mathbf{W}}\mu_1, \hat{\mathbf{W}}\mu_1$ is asymptotically non-zero in the LDR with $c > 0$, unless $\ell_1 = \ell_2$, in which case the original means $\mu_1, \mu_2$ are already orthogonal, as $\ell_1 = \ell_2$ implies $\theta/2 = \pi/4$, meaning $\theta = \pi/2$.

A numerical example is shown in Figure 1. We generate $N = 5000$ observations $\mathbf{x}_n$ following (1) with $K = 2$ and $P = 500$. Here, $\omega_1 = 0.567$ and $\omega_2 = 0.433$ and the means $\mu_1, \mu_2$ are unit-norm vectors satisfying $\mu_1^\top \mu_2 = 1/4$. We compute the estimated squared alignments $\hat{\zeta}_k = (\hat{\mathbf{u}}_k^\top \mathbf{u}_k)^2$ and compare them with the theoretical limiting quantities $\zeta_k$. The estimated residual alignment $\rho_{1,2}$ between $\hat{\mathbf{W}}\mu_1, \hat{\mathbf{W}}\mu_2$ is computed by (8) and compared to its theoretical counterpart $\rho_{1,2}^\infty$ given in (10). This is repeated for different levels of signal-to-noise ratio (SNR), here defined as $\text{SNR} := 1/\sigma^2$.

<figure><img src="image_3.png" /><figcaption>Fig. 1. Residual and eigenvector alignments in the LDR ($K = 2$).</figcaption></figure>

The estimated quantities are averaged over of 25 independent Monte Carlo runs (of $N = 5000$ samples each). Gray dashed vertical lines mark the two critical SNRs $S_k^* := (\sqrt{c} + 1)/\lambda_k$ dictated by the spiked model. For $1/\sigma^2 < S_1^*$, both spikes are subcritical and no signal can be found ($\zeta_1 = \zeta_2 = 0$). For $S_1^* < 1/\sigma^2 < S_2^*$, only the first eigenvector $\mathbf{u}_1$ is partially recovered ($\zeta_1 > 0 = \zeta_2$) and the whitened means become collinear. For $1/\sigma^2 > S_2^*$, both $\mathbf{u}_1, \mathbf{u}_2$ are partially recovered ($\zeta_1, \zeta_2 > 0$) but the whitened means still have nonzero residual alignment as predicted by Lemma 1. The observations closely track the theoretical curves, apart from some oscillation around the second phase transition, as usual in finite dimension.

## 3.2 Corrected whitening

Lemma 1 shows that the empirical whitening map $\hat{\mathbf{W}}$ fails to orthogonalize the means in the LDR because both eigenvalues and eigenvector alignments are distorted according to (6)–(7), implying that the coordinates of the whitened means are biased by some scalar factors. However, these factors can be predicted from (6)–(7). Hence, one can perform the following “scalar” correction: replace each empirical eigenvalue $\hat{\lambda}_k$ by a consistent estimator of $\lambda_k$ and compensate for eigenvector alignment $\hat{\mathbf{u}}_k^\top \mathbf{u}_k$ shrinkage by dividing it by a consistent estimator of $\psi_k$. Similarly corrected estimators based on RMT results have been derived in other contexts (see, e. g., [15], [16]).

Concretely, we proceed as follows. We may express $\lambda_k$ and $\beta_k$ as functions of $\ell_k$ in (6), and then express $\ell_k$ as a function of $\hat{\lambda}_k$ (when $\ell_k > \sqrt{c}$), which is tantamount to choosing the positive root of $\ell_k^2 + (1 + c - \hat{\lambda}_k/\sigma^2)\ell_k + c = 0$. By replacing $\sigma$ by $\hat{\sigma}$ and $\hat{\lambda}_k$ by $\hat{\lambda}_k$, we obtain the expression of the corrected spike:

$$ \hat{\ell}_k^c = \frac{1}{2} \left( \hat{\lambda}_k / \hat{\sigma}^2 - (1+c) + \sqrt{\left( \hat{\lambda}_k / \hat{\sigma}^2 - (1+c) \right)^2 - 4c} \right). $$

Then the corrected estimator of $\lambda_k$ is $\hat{\lambda}_k^c = \hat{\sigma}^2 (1 + \hat{\ell}_k^c)$. Next, define the finite-$(P, N)$ corrected versions of the overlap factors:

$$ \hat{\beta}_k^c = 1 + \frac{c}{\hat{\ell}_k^c}, \quad \hat{\psi}_k^c = 1 - \frac{\hat{\beta}_k^c - 1}{\hat{\beta}_k^c} \frac{1 + \hat{\ell}_k^c}{\hat{\ell}_k^c}. $$

Using these quantities, we can correct both eigenvalue inflation and eigenvector shrinkage via the diagonal scaling

$$ \Phi := \text{Diag} \left( \frac{1_{\{\hat{\ell}_1^c > \sqrt{c}\}}}{\sqrt{(\hat{\lambda}_1^c - \hat{\sigma}^2) \hat{\psi}_1^c}}, \dots, \frac{1_{\{\hat{\ell}_K^c > \sqrt{c}\}}}{\sqrt{(\hat{\lambda}_K^c - \hat{\sigma}^2) \hat{\psi}_K^c}} \right), $$
by simply defining the corrected whitening matrix as $\hat{\mathbf{W}}^c := \Phi \hat{\mathbf{U}}_K^\top \in \mathbb{R}^{K \times P}$. These definitions directly lead to a proof of the following theorem, which is our main result.

**Theorem 2** (Asymptotic dot products under corrected whitening). Assume (A1)-(A4) and consider the same notation as in Theorem 1. Let $\mathbf{H} := \text{Diag}(\mathbb{1}_{\{\ell_1 > \sqrt{c}\}}, \dots, \mathbb{1}_{\{\ell_K > \sqrt{c}\}})$. Then, for any $i, j \in \{1, \dots, K\}$,

$$ \langle \hat{\mathbf{W}}^c \mu_i, \hat{\mathbf{W}}^c \mu_j \rangle \xrightarrow{\text{a.s.}} \mu_i^\top \mathbf{U}_K \mathbf{H} \Gamma^{-1} \mathbf{U}_K^\top \mu_j. $$

In particular, if $\ell_k > \sqrt{c}$ for all $k = 1, \dots, K$, then

$$ \langle \hat{\mathbf{W}}^c \mu_i, \hat{\mathbf{W}}^c \mu_j \rangle \xrightarrow{\text{a.s.}} \langle \mathbf{W} \mu_i, \mathbf{W} \mu_j \rangle = \omega_i^{-1} \delta_{ij}. $$

# 4 Application to spherical GMM

We now illustrate the impact of our results, and in particular of our corrected whitening matrix, on the performance of a tensor-based algorithm for spherical GMM estimation. For concreteness, we will take the LEARNGMM algorithm of [5], which is arguably one of the simplest algorithms of this kind. In LEARNGMM, after the estimated whitening matrix $\hat{\mathbf{W}}$ is obtained from $\{\mathbf{x}_n\}_{n=1}^N$, a new set $\{\mathbf{x}_n\}_{n=N+1}^{2N}$ of $N$ independent samples is required to estimate the whitened tensor moment $\mathcal{M}_3(\mathbf{W}, \mathbf{W}, \mathbf{W})$ defined in (3). Letting $\xi_n = \hat{\mathbf{W}}\mathbf{x}_n$ for $n \in \{N+1, \dots, 2N\}$, the estimator $\hat{\mathcal{M}}_3(\hat{\mathbf{W}}, \hat{\mathbf{W}}, \hat{\mathbf{W}})$ of $\mathcal{M}_3(\mathbf{W}, \mathbf{W}, \mathbf{W})$ in LEARNGMM is given by

$$ \begin{aligned} \hat{\mathcal{M}}_3(\hat{\mathbf{W}}, \hat{\mathbf{W}}, \hat{\mathbf{W}}) &:= \frac{1}{N} \sum_{n=N+1}^{2N} \xi_n^{\otimes 3} \\ &\quad - \hat{\sigma}^2 \sum_{k=1}^K (\tilde{\xi} \otimes \mathbf{e}_k \otimes \mathbf{e}_k + \mathbf{e}_k \otimes \tilde{\xi} \otimes \mathbf{e}_k + \mathbf{e}_k \otimes \mathbf{e}_k \otimes \tilde{\xi}), \end{aligned} \quad (11) $$

where $\tilde{\xi} = \frac{1}{N} \sum_{n=N+1}^{2N} \xi_n$ and $\mathbf{e}_k$ denotes the $k$th canonical basis vector of $\mathbb{R}^K$. The authors of [5] then argued that $\hat{\mathcal{M}}_3(\hat{\mathbf{W}}, \hat{\mathbf{W}}, \hat{\mathbf{W}}) \to \mathcal{M}_3(\mathbf{W}, \mathbf{W}, \mathbf{W})$ as $N \to \infty$ with fixed $P$. One can then apply a spectral decomposition algorithm to $\mathcal{M}_3(\mathbf{W}, \mathbf{W}, \mathbf{W})$, thereby recovering the vectors $\mathbf{v}_k$, which can be transformed back into estimates of the component means $\mu_k$. In particular, the decomposition procedure used in LEARNGMM consists in computing the spectral decomposition of a symmetric matrix obtained by contracting $\mathcal{M}_3(\mathbf{W}, \mathbf{W}, \mathbf{W})$ with a random vector.

## 4.1 Corrected whitened third moment

In the LDR, the estimator $\hat{\mathcal{M}}_3(\hat{\mathbf{W}}, \hat{\mathbf{W}}, \hat{\mathbf{W}})$ of (11) no longer converges to the whitened moment tensor $\mathcal{M}_3(\mathbf{W}, \mathbf{W}, \mathbf{W})$, as a consequence of Lemma 1. Nonetheless, we claim in the following that this property can be recovered by simply replacing the whitening matrix by our corrected version $\hat{\mathbf{W}}^c$. Specifically, define the corrected estimator $\hat{\mathcal{M}}_3^c(\hat{\mathbf{W}}^c, \hat{\mathbf{W}}^c, \hat{\mathbf{W}}^c)$ by modifying the formula (11) as follows:

$$ \begin{aligned} \hat{\mathcal{M}}_3^c(\hat{\mathbf{W}}^c, \hat{\mathbf{W}}^c, \hat{\mathbf{W}}^c) &:= \frac{1}{N} \sum_{n=N+1}^{2N} (\xi_n^c)^{\otimes 3} \\ &\quad - \sum_{k=1}^K \hat{\sigma}^2 \left( \tilde{\xi}^c \otimes \mathbf{e}_k \otimes \mathbf{e}_k + \mathbf{e}_k \otimes \tilde{\xi}^c \otimes \mathbf{e}_k + \mathbf{e}_k \otimes \mathbf{e}_k \otimes \tilde{\xi}^c \right), \end{aligned} \quad (12) $$

where $\xi_n^c = \hat{\mathbf{W}}^c \mathbf{x}_n$ and $\tilde{\xi}^c = \frac{1}{N} \sum_{n=N+1}^{2N} \xi_n^c$.

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** GMM means estimation with and without corrected whitening.</figcaption></figure>

**Theorem 3** (Consistency of the corrected whitened third moment). Assume (A1)-(A4) and consider the same notation as in Theorem 1. If $\ell_k > \sqrt{c}$ for all $k \in \{1, \dots, K\}$, then

$$ \hat{\mathcal{M}}_3^c(\hat{\mathbf{W}}^c, \hat{\mathbf{W}}^c, \hat{\mathbf{W}}^c) \xrightarrow{\text{a.s.}} \mathcal{M}_3(\mathbf{W}, \mathbf{W}, \mathbf{W}), $$

where $\mathcal{M}_3(\mathbf{W}, \mathbf{W}, \mathbf{W})$ is as defined by (3).

**Remark 4.** While consistency of the low-dimensional object $\mathcal{M}_3(\mathbf{W}, \mathbf{W}, \mathbf{W})$ is ensured by using $\hat{\mathbf{W}}^c$, the component means $\mu_k$ are not consistently estimated because the signal subspace span$\{\mathbf{u}_1, \dots, \mathbf{u}_K\}$ cannot be perfectly recovered in the LDR, as per Theorem 1. Still, it can be shown that our correction allows recovering their projection onto the estimated subspace $\hat{\mathbf{U}}_K \hat{\mathbf{U}}_K^\top \mu_k$, which is the best that can be achieved when whitening is used.

## 4.2 Numerical illustration

We consider the estimation of the means in a spherical GMM with $K=2$, $N=2500$, equal weights $\omega_1 = \omega_2 = 0.5$ and random unit-norm means satisfying $\mu_1^\top \mu_2 = 0.5$. We compare the performance of the estimators described in Section 4.1, with and without corrected whitening, by computing $\|\hat{\mu}_k - \mu_k\|^2$ for $k=1, 2$. Fig. 2 displays the averaged results over 60 Monte Carlo realizations of the dataset, for $P=2500$ (left) and $P=500$ (right). The dashed vertical lines indicate the critical SNR levels defined in Section 3.1. As predicted, both methods perform poorly before the second phase transition $S_2^*$, since the eigenvectors $\mathbf{u}_1, \mathbf{u}_2$ cannot be both weakly recovered in this regime. Above $S_2^*$, both spikes are recoverable, yet the original algorithm using standard whitening still incurs in substantial errors due to a high residual alignment between whitened means. By contrast, the corrected whitening matrix restores (asymptotic) orthogonality, thereby significantly reducing estimation errors. At very high SNR, the residual alignment decays and both methods perform nearly identically. As seen in Fig. 2 (right), significant improvement with the corrected whitening is obtained on a smaller interval for $P=500$, since $c$ then gets closer to the classical regime $c \to 0$.

# 5 Conclusion

By exploiting results pertaining to spiked covariance matrix models coming from random matrix theory, we precisely characterized the loss of orthogonality of the whitened GMM means that takes place in the large-dimensional regime. This result allowed us to propose a simple scalar correction that restores their asymptotic orthogonality, leading to a significant performance improvement, especially under low signal-to-noise ratio, as illustrated by our numerical experiments.
# 6 References

[1] C. M. Bishop, *Pattern Recognition and Machine Learning*, Springer, 2006.

[2] P. Comon, "Independent component analysis, a new concept?," *Signal Processing*, vol. 36, no. 3, pp. 287-314, 1994.

[3] A. Hyvärinen and E. Oja, "Independent component analysis: Algorithms and applications," *Neural Networks*, vol. 13, no. 4-5, pp. 411-430, 2000.

[4] A. Anandkumar, R. Ge, D. J. Hsu, S. M. Kakade, and M. Telgarsky, "Tensor decompositions for learning latent variable models," *Journal of Machine Learning Research*, vol. 15, no. 1, pp. 2773-2832, 2014.

[5] D. J. Hsu and S. M. Kakade, "Learning mixtures of spherical Gaussians: Moment methods and spectral decompositions," in *Proc. Conference on Innovations in Theoretical Computer Science*, 2013, pp. 11-20.

[6] A. Moitra and G. Valiant, "Settling the polynomial learnability of mixtures of Gaussians," in *Proc. Annual Symposium on Foundations of Computer Science*. IEEE, 2010, pp. 93-102.

[7] A. T. Kalai, A. Moitra, and G. Valiant, "Efficiently learning mixtures of two Gaussians," in *Proc. ACM Symposium on Theory of Computing*, 2010, pp. 553-562.

[8] O. Ledoit and M. Wolf, "A well-conditioned estimator for large-dimensional covariance matrices," *Journal of Multivariate Analysis*, vol. 88, no. 2, pp. 365-411, 2004.

[9] Z. Bai and J. W. Silverstein, *Spectral Analysis of Large Dimensional Random Matrices*, vol. 20, Springer, 2010.

[10] R. Couillet and Z. Liao, *Random Matrix Methods for Machine Learning*, Cambridge University Press, 2022.

[11] I. M. Johnstone, "On the distribution of the largest eigenvalue in principal components analysis," *Annals of Statistics*, vol. 29, no. 2, pp. 295-327, 2001.

[12] F. Benaych-Georges and R. R. Nadakuditi, "The eigenvalues and eigenvectors of finite, low rank perturbations of large random matrices," *Advances in Mathematics*, vol. 227, no. 1, pp. 494-521, 2011.

[13] D. Paul, "Asymptotics of sample eigenstructure for a large dimensional spiked covariance model," *Statistica Sinica*, vol. 17, no. 4, pp. 1617-1642, 2007.

[14] J. Baik and J. W. Silverstein, "Eigenvalues of large sample covariance matrices of spiked population models," *Journal of Multivariate Analysis*, vol. 97, no. 6, pp. 1382-1408, 2006.

[15] R. R. Nadakuditi, "Optshrink: An algorithm for improved low-rank signal matrix denoising by optimal, data-driven singular value shrinkage," *IEEE Transactions on Information Theory*, vol. 60, no. 5, pp. 3002-3018, 2014.

[16] M. Gavish and D. L. Donoho, "The optimal hard threshold for singular values is $4/\sqrt{3}$," *IEEE Transactions on Information Theory*, vol. 60, no. 8, pp. 5040-5053, 2014.