1 Introduction↩︎

Model order selection—the problem of estimating the number of statistically significant components in a noisy observation—is fundamental in signal processing, with applications in wireless communications [1], array processing [2], and many other areas [3], [4]. In the classical white-noise setting, the model order can be estimated from the eigenvalue structure of the sample covariance matrix via information-theoretic criteria such as the Akaike Information Criterion (AIC) [5], and related subspace methods have been applied to source localization [6], channel identification [7], and waveform estimation [8]. However, these methods become unreliable in the presence of large-dimensional and correlated data. Although some extensions have been proposed for correlated signals [9], [10], they cannot be generalized to arbitrary scenarios, and whitening requires knowledge of the noise covariance structure, which is typically unknown [11]. Moreover, in high-dimensional regimes, the sample covariance matrix no longer exhibits the same properties as in the classical setting, leading to poor estimation performance [12]–[14]. In the large-dimensional regime, where both the number of snapshots \(N\) and the signal dimension \(m\) grow to infinity with a constant ratio, Random Matrix Theory provides powerful tools for model order selection. In particular, methods based on the asymptotic behavior of the largest eigenvalues of the covariance matrix have been proposed [15]. RMT offers a rigorous probabilistic framework to analyze large random matrices; see [16] for a comprehensive review and, e.g., [13], [17]–[20] for applications to detection, MUSIC, radar, and hyperspectral imaging, respectively. When the noise is correlated, model order estimation remains possible, for instance by analyzing the spacing between eigenvalues [21], [22]. Moreover, many practical signal processing applications involve impulsive, heavy-tailed, or non-Gaussian perturbations that cannot be accurately modeled within the classical Gaussian framework. This motivates the use of more general distribution families together with robust covariance estimation techniques. When correlation and heavy-tailed behaviour occur jointly, both covariance estimation and eigenvalue-based model-order selection become significantly more challenging. To mitigate sensitivity to the noise distribution, robust model order selection methods have been proposed, notably in hyperspectral imaging [23], [24]. Nevertheless, these approaches often depend on unknown parameters [1] or are not well-suited for large-dimensional settings. Recent advances in RMT enable consistent estimation of covariance (or scatter) matrices for textured signals [25]. However, these results typically assume that the noise correlation structure is known and that the data can be whitened beforehand, which is rarely the case in practice. The present paper addresses both challenges simultaneously.

In this work, we consider noise following a Complex Elliptically Symmetric (CES) distribution [26], [27] (see also [28] for a comprehensive overview). CES distributions are widely used in signal processing due to their flexibility, as they can model a broad class of random signals beyond the Gaussian assumption. A CES random vector can be decomposed into two components: a texture and a speckle. Its second-order structure is characterized by a scatter matrix, which coincides with the covariance matrix up to a scaling factor. CES models have been successfully applied in various fields, including hyperspectral imaging [29], radar clutter modeling [30], and finance, particularly in portfolio optimization [31].

This paper proposes a robust two-stage framework that combines Toeplitz-structured scatter matrix estimation, whitening, and RMT-based eigenvalue thresholding to estimate the number of sources in correlated CES noise, extending [22], [25], [32] to left-sided correlation with an unknown scatter matrix under non-Gaussian noise.
Contributions:

(i) Almost sure consistency of Toeplitz-structured scatter matrix estimators under correlated CES noise, even in the presence of sources: SCM-based and Maronna-based (Theorems 1 and 3, Supplementary Material) and Tyler-based (Theorem 5 [33]).

(ii) Robust whitening via Toeplitz rectification of the SCM (Section 3-A) or Maronna’s \(M\)-estimator (Section 4-A), with asymptotic validity.

(iii) Explicit closed-form RMT eigenvalue thresholds: Maronna-based threshold \(t = \Phi_\infty(1+\sqrt{c})^2 / [\xi\,(1-c\,\Phi_\infty)]\) (Theorems 2 and 4) and Marčenko–Pastur upper edge \((1+\sqrt{c})^2\) for the distribution-free Tyler branch (Theorem 6).

Detailed proofs are provided in the Supplementary Material (Sections S-A and S-B).

Notations: Matrices are in bold and capital, vectors in bold. Let \(\mathbf{X}\) be a square matrix of size \(s\times s\), \(\left(\lambda_{i}(\mathbf{X})\right)_i\), \(i \in \{1,\ldots,s\}\), are the eigenvalues of \(\mathbf{X}\). \(\operatorname{tr}(\mathbf{X})\) is the trace of the matrix \(\mathbf{X}\). \(\left\Vert \mathbf{X} \right\Vert\) stands for the spectral norm. Let \(\mathbf{A}\) be a matrix, \(\mathbf{A}^T\) is the transpose of \(\mathbf{A}\) and \(\mathbf{A}^H\) the Hermitian transpose of \(\mathbf{A}\). \(\mathbf{I}_p\) is the \(p \times p\) identity matrix. For any \(m\)-vector \(\mathbf{x}\), \(\mathcal{L} :\, \mathbf{x} \mapsto \mathcal{L}(\mathbf{x})\) is the \(m \times m\) matrix defined as the Toeplitz operator: \(\left([\mathcal{L}(\mathbf{x})]_{i,j}\right)_{i\leq j} = x_{i-j}\) and \(\left([\mathcal{L}(\mathbf{x})]_{i,j}\right)_{i>j} = x_{i-j}^\ast\). For \(\mathbf{A} \in \mathbb{C}^{m \times m}\), the Toeplitz rectification \(\mathcal{T}(\mathbf{A})\) is defined as \(\mathcal{L}(\tilde{\mathbf{a}})\) with \(\tilde{a}_k = \displaystyle\frac{1}{m} \sum_{i-j=k} A_{i,j}\) for \(k \in \{0, \dots, m-1\}\). For any complex \(z\), \(z^{\star}\) is the conjugate of \(z\). \(\mathcal{R}e(\cdot)\) and \(\mathcal{I}m(\cdot)\) denote the real and imaginary parts. The notation \(\overset{a.s.}{\longrightarrow}\) means “converges almost surely”. \(C\) generally denotes a real positive constant.

2 Model and Assumptions↩︎

This section presents the signal model and the assumptions required for the theoretical analysis. We consider the following standard source-plus-noise model. Let \(\mathbf{y}_0, \ldots, \mathbf{y}_{N-1}\) be \(N\) observations of dimension \(m\), corresponding to \(p\) mixed sources embedded in additive CES noise: \[\mathbf{y}_i = \mathbf{M}\, \mathbf{s}_i + \sqrt{\tau_i} \,\mathbf{C}^{1/2} \, \mathbf{x}_i\,, \quad i \in \{0,\ldots,N-1\} \,, \label{modele}\tag{1}\] The first term represents \(p\) deterministic sources, while the second term models CES noise with texture \(\tau_i\) and scatter matrix \(\mathbf{C}\). We can rewrite this equation in a more compact form \(\mathbf{Y} = \mathbf{M}\, \mathbf{S} + \mathbf{C}^{1/2}\, \mathbf{X}\, \mathbf{T}^{1/2}\), where

• \(\mathbf{Y}=[\mathbf{y}_0,\ldots,\mathbf{y}_{N-1}] \in \mathbb{C}^{m\times N}\),

• \(\mathbf{M}\in \mathbb{C}^{m\times p}\) is the mixing matrix, whose columns \(\left\{\mathbf{m}_j\right\}_{j\in\{1,\ldots,p\}}\) represent the source signatures.

• The matrix \(\mathbf{S} = \left[\mathbf{s}_0, \ldots, \mathbf{s}_{N-1}\right] = \boldsymbol{\Gamma}^{1/2}\, \boldsymbol{\delta}^H \in \mathbb{C}^{p \times N}\) with \(\boldsymbol{\delta} \in\mathbb{C}^{N\times p}\) having i.i.d.\(\mathcal{CN}(0,1)\) entries, independent of \(\mathbf{X}\), contains the source amplitude where each entry \(s_{i,j}\) represents the contribution of source \(j\) in observation \(i\).

• \(\boldsymbol{\Gamma}\in\mathbb{C}^{p\times p}\) is Hermitian nonnegative definite,

• \(\mathbf{X}=[\mathbf{x}_0,\ldots,\mathbf{x}_{N-1}]\in\mathbb{C}^{m\times N}\) has i.i.d.entries with zero mean, unit variance, and sub-Gaussian tails,

• \(\mathbf{T}=\mathrm{diag}(\tau_0,\ldots,\tau_{N-1})\), with \(\{\tau_i\}\) i.i.d.positive random variables referred to as textures, with an unspecified distribution.

• \(\mathbf{C} = \mathcal{L}\bigl([c_0,\ldots,c_{m-1}]^T\bigr)\in \mathbb{C}^{m\times m}\) is a Hermitian nonnegative definite Toeplitz scatter matrix. This structural assumption facilitates estimation and is consistent with many practical applications.

In the sequel, we will consider the following assumptions:

The overall procedure is summarized in Fig. 1. The three whitening strategies correspond to the SCM-based, Maronna-based, and Tyler-based estimators, which respectively provide a baseline, a robust parametric, and a distribution-free approach.

3 SCM-Based Whitening with Robust Scatter Estimation↩︎

In this section, we exploit the consistency of the Sample Covariance Matrix (SCM) to perform a whitening step and subsequently estimate the model order using Maronna’s \(M\)-estimator. Model order estimation via Maronna’s \(M\)-estimator was studied in [25] for the spiked model with CES white noise. Here we address the more challenging scenario of correlated CES noise. For clarity, Table 1 summarizes the main matrices used throughout the paper.

3.1 Whitening Step↩︎

We set the following two matrices: \[\check{\mathbf{C}}_{SCM} = \mathcal{T}\!\left(\hat{\mathbf{C}}_{SCM}\right), \qquad \hat{\mathbf{C}}_{SCM} = \frac{1}{N}\,\mathbf{Y}\,\mathbf{Y}^H\, . \label{eq1}\tag{2}\]

The whitened observations are: \[\check{\mathbf{Y}}_{wSCM} = \check{\mathbf{C}}_{SCM}^{-1/2} \, \mathbf{M} \, \boldsymbol{\Gamma}^{1/2}\boldsymbol{\delta}^H + \check{\mathbf{C}}_{SCM}^{-1/2} \, \mathbf{C}^{1/2}\,\mathbf{X} \, \mathbf{T}^{1/2}\,. \label{scm95blanc}\tag{3}\]

3.2 Estimation of the Scatter Matrix↩︎

Having established that \(\check{\mathbf{C}}_{SCM}\) consistently estimates \(\mathbf{C}\), we apply Maronna’s \(M\)-estimator [39] to the whitened observations \(\check{\mathbf{Y}}_{wSCM}\) in 3 to robustly estimate the scatter matrix of the decorrelated signal. The estimator \(\check{\boldsymbol{\Sigma}}_{SCM}\) is the unique solution of: \[\boldsymbol{\Sigma} = \frac{1}{N} \displaystyle\sum_{i=0}^{N-1} u\!\left(\frac{1}{m}\,\check{\mathbf{y}}_{wSi}^H\,\boldsymbol{\Sigma}^{-1}\, \check{\mathbf{y}}_{wSi}\right) \check{\mathbf{y}}_{wSi} \, \check{\mathbf{y}}_{wSi}^H \,, \label{eqsigma1}\tag{4}\] where \(\check{\mathbf{Y}}_{wSCM} = [\check{\mathbf{y}}_{wS0}, \ldots, \check{\mathbf{y}}_{wS\, N-1}]\) and \(u:[0,+\infty)\to (0,+\infty)\) satisfies:

(i) \(u\) is nonnegative, continuous, and non-increasing [40];

(ii) \(\phi(x) = x\,u(x)\) is increasing and bounded, with \(\displaystyle\lim_{x\to\infty}\phi(x) = \Phi_{\infty} > 1\);

(iii) \(c < \Phi_{\infty}^{-1}\) (cf.Assumption 2-(v)).

These conditions guarantee the existence and uniqueness of the solution [39], [40].

3.3 Model Order Selection↩︎

3.3.0.1 Threshold Derivation

To count sources, we study the eigenvalue distribution of \(\check{\boldsymbol{\Sigma}}_{SCM}\). A direct RMT analysis is intractable because the weight \(u\!\left(\displaystyle\frac{1}{m}\,\check{\mathbf{y}}_{wSi}^H\, \check{\boldsymbol{\Sigma}}_{SCM}^{-1}\,\check{\mathbf{y}}_{wSi}\right)\) depends on the very vector \(\check{\mathbf{y}}_{wSi}\) it weights, creating a circular dependency. To break this dependency, we introduce a surrogate matrix defined on the ideally whitened signal. We therefore work with an auxiliary matrix defined on the ideally whitened signal \(\mathbf{Y}_w = \mathbf{C}^{-1/2}\, \mathbf{Y}=[\mathbf{y}_{w0}, \ldots, \mathbf{y}_{w N-1}]\): \[\mathbf{Y}_w = \mathbf{C}^{-1/2}\,\mathbf{M}\,\boldsymbol{\Gamma}^{1/2}\, \boldsymbol{\delta}^H + \mathbf{X}\,\mathbf{T}^{1/2} \,. \label{scm95normal}\tag{5}\]

To overcome the dependence between the weights and the data, we introduce the surrogate matrix: \[\hat{\mathbf{S}} \triangleq \frac{1}{N}\sum_{i=0}^{N-1} v(\tau_i\,\xi)\,\mathbf{y}_{wi}\,\mathbf{y}_{wi}^H = \frac{1}{N}\sum_{i=0}^{N-1} v(\tau_i\,\xi)\,\tau_i\,\mathbf{x}_i\,\mathbf{x}_i^H. \label{matrixS}\tag{6}\] with \(g(x) = x/(1-c\,\phi(x))\), which is strictly increasing on \([0,+\infty)\) under Assumption 2-(v) (since \(c\,\phi(x) < c\,\Phi_\infty < 1\) for all \(x\)), and hence invertible on its range. Let \(\psi(x)=x\,v(x)\), and let \(\xi>0\) be the unique positive solution of \[\displaystyle\frac{1}{N}\sum_{i=0}^{N-1} \frac{\psi(\tau_i\,\xi)}{1+c\,\psi(\tau_i\,\xi)} = 1\, . \label{xi}\tag{7}\] Existence and uniqueness of \(\xi\) follow from the monotonicity of the left-hand side: since \(\psi(0)=0\) and \(\psi(\tau_i\,\xi)\to\Phi_\infty\) as \(\xi\to\infty\), the left-hand side of 7 increases continuously from \(0\) to \(\displaystyle\frac{\Phi_\infty}{1+c\,\Phi_\infty}\), which exceeds \(1\) under Assumption 2-(v) (since \(c < \Phi_\infty^{-1}\) implies \(\Phi_\infty(1-c)>1\)), guaranteeing existence and uniqueness of \(\xi\) (see [25] for details). The weights \(\{v(\tau_i\,\xi)\}\) are independent of \(\{\mathbf{x}_i\}\), which makes \(\hat{\mathbf{S}}\) more tractable. This decoupling is key to enabling a tractable RMT analysis.

It is proven in [25] that \(\left\|\hat{\boldsymbol{\Sigma}} - \hat{\mathbf{S}}\right\|\overset{a.s.}{\to}0\), where \(\hat{\boldsymbol{\Sigma}}\) is Maronna’s \(M\)-estimator applied to the ideally whitened data \(\mathbf{Y}_w\). The key insight is that \(\hat{\mathbf{S}}\) replaces the data-dependent weights \(u\!\left(\displaystyle\frac{1}{m}\,\mathbf{y}_{wi}^H\,\hat{\boldsymbol{\Sigma}}^{-1}\,\mathbf{y}_{wi}\right)\) by the oracle weights \(v(\tau_i\,\xi)\), which are independent of \(\mathbf{x}_i\), making \(\hat{\mathbf{S}}\) tractable for RMT analysis. The following theorem then extends this to the empirically whitened case.

3.3.0.2 Source-free upper bound

Since \(\psi(x)=x\,v(x)\) by definition, one has \(v(\tau_i\,\xi)\,\tau_i = v(\tau_i\,\xi)\,\tau_i\,\xi/\xi = \psi(\tau_i\,\xi)/\xi\). In the source-free case (\(\mathbf{M}=\mathbf{0}\)), Eq. 6 becomes \[\hat{\mathbf{S}} = \frac{1}{N}\sum_{i=0}^{N-1} \frac{\psi(\tau_i\,\xi)}{\xi}\,\mathbf{x}_i\mathbf{x}_i^H. \label{eq:Shat95sourcefree}\tag{8}\] Since \(\{\tau_i\}\) are independent of \(\{\mathbf{x}_i\}\) and \(0 < \psi(\tau_i\,\xi)/\xi \leq \Phi_\infty/\xi\) almost surely, a standard comparison argument gives \(\left\|\hat{\mathbf{S}}\right\| \leq \displaystyle\frac{\Phi_\infty}{\xi} \left\|\frac{1}{N}\,\mathbf{X}\mathbf{X}^H\,\right\|\) (see Supplementary Material S-A.2). More precisely, since \(v(\tau_i\,\xi)\leq\Phi_\infty/\xi\) and the weights are independent of \(\{\mathbf{x}_i\}\), the matrix inequality \(\hat{\mathbf{S}}\preceq \displaystyle\frac{\Phi_\infty}{\xi}\,\frac{1}{N}\,\mathbf{X}\,\mathbf{X}^H\) holds almost surely entry-wise in a quadratic form sense, giving \(\|\hat{\mathbf{S}}\|\leq \displaystyle\frac{\Phi_\infty}{\xi} \left\|\frac{1}{N}\,\mathbf{X}\,\mathbf{X}^H\right\|\). Applying the Bai–Silverstein almost sure upper bound [9], which guarantees that \(\lambda_{\max}\!\left(\displaystyle\frac{1}{N} \,\mathbf{X}\, \mathbf{X}^H\right)\overset{a.s.}{\longrightarrow}\left(1+\sqrt{c}\right)^2\), then yields: \[\left\Vert \hat{\mathbf{S}} \right\Vert \leq t, \qquad t = \frac{\Phi_{\infty}\,(1+\sqrt{c})^2}{\xi\,(1-c\,\Phi_{\infty})} > 0\,, \label{eqseuil}\tag{9}\] where \(c\,\Phi_\infty<1\) is guaranteed by Assumption 2-(v). The additional factor \((1-c\,\Phi_\infty)^{-1}\) in 9 arises from the definition \(v = u\,\circ \,g^{-1}\) with \(g(x) = x/(1-c\,\phi(x))\): one has \(v(x)\leq \Phi_\infty/\xi\) but more precisely \(\psi(x)=x\,v(x)\leq \Phi_\infty\) and \(\xi\geq \xi_{\min}>0\), so that \(v(\tau_i\,\xi)\,\tau_i = \psi(\tau_i\,\xi)/\xi \leq \Phi_\infty/\xi\), and the threshold follows from \(\|\hat{\mathbf{S}}\|\leq (\Phi_\infty/\xi)(1+\sqrt{c})^2\) combined with the normalization \(\xi\,(1-c\,\Phi_\infty)\geq \xi_{\min}\,(1-c\,\Phi_\infty)>0\) from the fixed-point equation 7 . The bound 9 holds almost surely as \(N,m\to\infty\); for fixed \((N,m)\), a finite-sample correction may be needed. In practice, \(\xi\) is unknown; a consistent leave-one-out estimator is \[\hat{\xi} = \displaystyle \frac{1}{N}\sum_{i=0}^{N-1}\frac{1}{m}\, \mathbf{y}_i^H\,\check{\boldsymbol{\Sigma}}_{(i)}^{-1}\,\mathbf{y}_i\, ,\] where \(\check{\boldsymbol{\Sigma}}_{(i)}\) is the leave-one-out approximation of the scatter estimator with observation \(i\) removed (see Section 6.1 for the explicit formula and efficient computation).

In practice, finite-sample exceedances above \(t\) may occur with empirical whitening; Theorem 1 guarantees asymptotic validity.

3.3.0.3 Model order estimator.

If sources are present with sufficiently high SNR, some eigenvalues of \(\check{\boldsymbol{\Sigma}}_{SCM}\) exceed \(t\). The estimated number of sources is defined as \(\hat{p} = \#\bigl\{k : \lambda_k(\check{\boldsymbol{\Sigma}}_{SCM}) > t\bigr\}\) where \(\lambda_1(\cdot) \geq \lambda_2(\cdot) \geq \cdots\) denotes the eigenvalues in decreasing order.

3.4 Results↩︎

The source-free configuration (\(\mathbf{M}=\mathbf{0}\)) is considered here in order to validate the theoretical noise eigenvalue upper bound 9 independently of signal subspace effects. The parameters are set to \(c=0.45\), \(m=900\), and \(N=2000\). Thus, \(\mathbf{Y} = \mathbf{C}^{1/2}\,\mathbf{X}\,\mathbf{T}^{1/2}\) with \(\mathbf{C} = \mathcal{L}\left((\rho^0,\rho^1,\ldots,\rho^{m-1})^T\right)\) where \(\rho=0.7\) and \(\mathbf{X}\) is a zero-mean complex Gaussian noise with identity scatter matrix. The texture matrix \(\mathbf{T}\) is a diagonal \(N\times N\)-matrix containing the \(\{\tau_i\}_{i\in\{0,\ldots,N-1\}}\) on its diagonal where \(\{\tau_i\}_i\) are i.i.d.inverse gamma distributed with mean equal to \(1\) and with a shape parameter equal to \(10\). The function \(u\) is here defined as \(u: x \mapsto \displaystyle\frac{1+\nu}{x+\nu}\) where \(\nu\) is a fixed parameter equal to \(0.1\) and corresponds to the Maximum Likelihood Estimation (MLE) for the chosen distribution of the \(\tau\) with \(\Phi_{\infty} = 1+\nu\).

Fig. 2 illustrates both the validity of Theorem 2 and the critical role of the whitening step. In the top panel, the data \(\mathbf{Y}\) have been whitened by \(\check{\mathbf{C}}_{SCM}^{-1/2}\): the log-eigenvalue distribution of \(\check{\boldsymbol{\Sigma}}_{SCM}\) (green) closely tracks that of the surrogate matrix \(\hat{\mathbf{S}}\) (blue), in agreement with the almost sure convergence ?? . All eigenvalues remain below the RMT threshold \(t\) (red vertical line), validating the source-free upper bound 9 with \(\Phi_\infty = 1+\nu\). The bottom panel shows the same data without whitening: the eigenvalues of \(\check{\mathbf{C}}_{SCM}^{1/2}\, \check{\boldsymbol{\Sigma}}_{SCM}\, \check{\mathbf{C}}_{SCM}^{1/2}\) (green) extend well beyond \(t\), confirming that the threshold 9 has no theoretical justification unless the correlation structure has been removed beforehand. Together, the two panels demonstrate that the robust \(M\)-estimation step effectively absorbs the impact of large texture values \(\tau_i\), while the whitening step is indispensable for calibrating the eigenvalue threshold.

Fig. 3 presents the log-eigenvalue distributions of \(\hat{\mathbf{S}}\) and \(\check{\boldsymbol{\Sigma}}_{SCM}\) for samples following a different CES distribution. The texture matrix \(\mathbf{T}\) is diagonal with i.i.d. entries \({\tau_i}_{i\in\{0,\ldots,N-1\}}\), where each \(\tau_i=t^2\) and \(t\) follows a Student-\(t\) distribution with \(100\) degrees of freedom and \(\nu=0.1\). Compared to Fig.2-Left, the eigenvalues are less concentrated around those of \(\hat{\mathbf{S}}\) and move closer to the threshold \(t\). As the distribution of \(\tau\) deviates from the one for which \(u\) corresponds to the MLE, the convergence in Theorem 2 becomes slower, making the results less reliable for fixed \(N\) and \(m\).

The key limitation of the SCM branch lies in the slow convergence rate of ?? when the texture distribution departs from the MLE design of \(u\), as illustrated in Fig. 3. This motivates Section 4, which replaces the SCM-based whitening by a Maronna \(M\)-estimator whitening, yielding improved finite-sample robustness at the cost of a mild additional computational overhead.

4 Robust Model Order Selection via \(M\)-Estimators↩︎

4.1 Whitening Step↩︎

Let \(\widetilde{\mathbf{C}}_{FP}\) be a scaled estimator of the scatter matrix \(\mathbf{C}\) defined by \(\widetilde{\mathbf{C}}_{FP} = \mathcal{T}\left(\hat{\mathbf{C}}_{FP}\right)\), where \(\hat{\mathbf{C}}_{FP}\) is the unique solution to the Maronna’s \(M\)-estimator [39]: \[\boldsymbol{\Sigma} = \frac{1}{N}\displaystyle\sum_{i=0}^{N-1} u\!\left(\frac{1}{m}\,\mathbf{y}_i^H\,\boldsymbol{\Sigma}^{-1}\,\mathbf{y}_i\right) \mathbf{y}_i\,\mathbf{y}_i^H \,.\] As in the previous section, \(u: [0,+\infty) \mapsto (0,+\infty)\) is nonnegative, continuous, and non-increasing. The following theorem stands for \(\widetilde{\mathbf{C}}_{FP}\):

As in the previous section, the samples \(\mathbf{Y}\) can then be whitened thanks to \(\check{\mathbf{C}}_{FP}^{-1/2}\). Let \(\check{\mathbf{Y}}_{wFP} = [\check{\mathbf{y}}_{wF0},\ldots, \check{\mathbf{y}}_{wF\,N-1}]\) be the whitened samples: \[\check{\mathbf{Y}}_{wFP} = \check{\mathbf{C}}_{FP}^{-1/2}\,\mathbf{M}\,\boldsymbol{\Gamma}^{1/2}\boldsymbol{\delta}^H + \check{\mathbf{C}}_{FP}^{-1/2}\,\mathbf{C}^{1/2}\,\mathbf{X}\,\mathbf{T}^{1/2} \,.\]

4.2 Robust Estimation of the Scatter Matrix↩︎

The robust estimation of the scatter matrix and the model order selection are performed as previously. The robust estimator of the scatter matrix of the whitened signal \(\check{\mathbf{Y}}_{wFP}\) is a Fixed-Point estimator denoted by \(\check{\boldsymbol{\Sigma}}_{FP}\) and defined as the unique solution of the equation: \[\boldsymbol{\Sigma} = \frac{1}{N}\displaystyle\sum_{i=0}^{N-1} u\!\left(\frac{1}{m}\,\check{\mathbf{y}}_{wFi}^H\,\boldsymbol{\Sigma}^{-1}\, \check{\mathbf{y}}_{wFi}\right) \check{\mathbf{y}}_{wFi}\,\check{\mathbf{y}}_{wFi}^H \,. \label{eqsigma2}\tag{10}\] Thus, \(\check{\boldsymbol{\Sigma}}_{FP}\) is a robust estimator of the scatter matrix of the whitened signal. Through Theorem 4 proposed below, Theorem 2 extends the consistency result of \(\check{\boldsymbol{\Sigma}}_{SCM}\) to the Maronna-whitened estimator \(\check{\boldsymbol{\Sigma}}_{FP}\).

4.3 Model Order Selection↩︎

The threshold \(t\) given in Equation 9 holds asymptotically, with the proof following the same perturbation argument as in Supplementary Material S-B, and can be used on the eigenvalues of \(\check{\boldsymbol{\Sigma}}_{FP}\) to estimate \(p\).

If sources with sufficiently high SNR are present, some eigenvalues of \(\check{\boldsymbol{\Sigma}}_{FP}\) exceed the threshold \(t\) defined in 9 . The estimated number of sources is \(\hat{p} = \#\bigl\{k : \lambda_k(\check{\boldsymbol{\Sigma}}_{FP}) > t\bigr\}\).

4.4 Results↩︎

Fig. 4 displays the log-eigenvalue distributions of \(\hat{\mathbf{S}}\) and \(\check{\boldsymbol{\Sigma}}_{FP}\) under the same simulation setup as Section 3 (\(c=0.45\), \(m=900\), \(N=2000\), \(\rho=0.7\), \(\nu=0.1\)), with the same weight function \(u(\cdot)\) and data whitened through \(\check{\mathbf{C}}_{FP}^{-1/2}\). In the left panel (\(\mathrm{InvGamma}\) textures, \(\Phi_\infty=1{+}\nu\)), the eigenvalue distribution of \(\check{\boldsymbol{\Sigma}}_{FP}\) closely tracks that of \(\hat{\mathbf{S}}\) and all eigenvalues remain below the threshold \(t\) 9 , validating Theorem 4. The convergence is noticeably tighter than in Fig. 2-Top (SCM branch), consistent with the faster convergence rate of ?? relative to ?? . The right panel (\(\mathrm{Student-t^2}\) textures, \(t\sim\mathrm{Student}\) with 100 degrees of freedom) shows that the Maronna-whitened estimator maintains reliable eigenvalue control even when the texture distribution departs from the MLE design of \(u\)—in contrast to the SCM branch (Fig. 3), where the same mismatch causes eigenvalues to approach \(t\). This robustness to texture distribution mismatch is a direct consequence of the down-weighting mechanism built into Maronna’s \(M\)-estimator during the whitening step itself.

Fig. 5 presents the same setup as in Fig. 4-left, but with three sources having signal-to-noise ratios of \([5, 10, 15]\,\)dB. It can be seen that three distinct groups of eigenvalues (10 Monte-Carlo runs) exceed the threshold, and that the noise eigenvalue distribution of \(\check{\boldsymbol{\Sigma}}_{FP}\) closely matches that of \(\hat{\mathbf{S}}\).

The results obtained with the robust method improve upon those of the previous section (e.g., Fig.3). The proposed approach is robust to the distribution of \(\tau\): even when the texture distribution differs from that assumed in the design of \(u\), the method remains reliable. This behavior stems from the robustness of the scatter matrix estimation. In contrast, without whitening, several eigenvalues exceed the threshold, preventing reliable conclusions or model order estimation. These results extend [22], [32] to the case of left-correlated noise. Moreover, the \(L_2\)-norm of the difference between the estimated scatter matrix and the SCM tends to zero as \(N,m\to\infty\) with constant ratio \(c\). Since many signal subspace rank estimation methods rely on the eigenvalues of the scatter matrix, the proposed estimator improves their consistency.

While Sections 3–4 require specifying a weight function \(u\) tuned to the texture distribution—a requirement that may be difficult to fulfill in practice—Section 5 eliminates this constraint by adopting Tyler’s \(M\)-estimator, which is distribution-free and achieves a universal Marčenko–Pastur threshold.

5 Tyler \(M\)-Estimator: A Distribution-Free Approach↩︎

A major practical challenge in robust signal processing is the lack of prior knowledge regarding the noise texture distribution. To circumvent this, we specialize our framework to Tyler’s \(M\)-estimator [42], which is distribution-free [43], and does not require specifying a weight function \(u\). Unlike Maronna’s framework of Sections III–IV, it was proved in [44] that the eigenvalue distribution of Tyler’s \(M\)-estimator converges, in the white CES case, to the Marčenko–Pastur law [45]; the detection threshold is therefore the upper edge \((1+\sqrt{c})^2\), which replaces 9 . The proposed method proceeds in two steps: (i) whitening via a Toeplitz-rectified Tyler estimator of \(\mathbf{C}\), and (ii) eigenvalue thresholding of the Tyler scatter matrix of the whitened signal at \((1+\sqrt{c})^2\) to infer the number of sources.

5.1 Signal Whitening↩︎

The first step of the procedure is to whiten the signal with a consistent estimation of \(\mathbf{C}\), the scatter matrix of the CES noise. The proposed estimator \(\hat{\mathbf{C}}_{TYL}\) is a Tyler \(M\)-estimator of the scatter matrix \(\mathbf{C}\) enforced to be Toeplitz-structured with the operator \(\mathcal{T} \left( \hat{\mathbf{C}}_{TYL} \right)\). The Tyler \(M\)-estimator \(\hat{\mathbf{C}}_{TYL}\) is defined as the unique solution of [42]: \[\boldsymbol{\Sigma} = \frac{m}{N} \displaystyle \sum_{i=0}^{N-1} \frac{\mathbf{y}_i \,\mathbf{y}_i^H}{\mathbf{y}_i^H \,\boldsymbol{\Sigma}^{-1} \,\mathbf{y}_i} \, . \label{eqTyler}\tag{11}\]

This equation requires an iterative approach to compute \(\hat{\mathbf{C}}_{TYL}\) as the final value of \(\boldsymbol{\Sigma}\). Note that the \(\mathbf{y}_i\) also contain the sources. The consistency of \(\mathcal{T} \left( \hat{\mathbf{C}}_{TYL} \right)\) is proven thanks to the following theorem:

Set \(\mathbf{\check{C}}_{TYL}=\mathcal{T} \left( \hat{\mathbf{C}}_{TYL}\right)\). The Tyler-branch whitened data are: \[\begin{align} \hat{\mathbf{Y}}_w &= \mathbf{\check{C}}_{TYL}^{-1/2}\,\mathbf{Y} \\ &= \mathbf{\check{C}}_{TYL}^{-1/2}\,\mathbf{M}\,\boldsymbol{\Gamma}^{1/2}\, \boldsymbol{\delta}^H + \mathbf{\check{C}}_{TYL}^{-1/2} \,\mathbf{C}^{1/2}\,\mathbf{X}\,\mathbf{T}^{1/2}. \end{align}\] with \(\hat{\mathbf{Y}}_w = [ \, \hat{\mathbf{y}}_{w0}, \, ... , \, \hat{\mathbf{y}}_{wN-1}]\).

5.2 Estimation of the Tyler Scatter Matrix↩︎

The signal being whitened, it is now possible to apply a Tyler \(M\)-estimator to threshold its eigenvalues. Let \(\check{\boldsymbol{\Sigma}}_{TYL}\) denote a Tyler estimation of the scatter matrix of the whitened data \(\hat{\mathbf{Y}}_w\), that is, \(\check{\boldsymbol{\Sigma}}_{TYL}\) is the unique solution if it exists of: \[\boldsymbol{\Sigma} = \frac{m}{N} \displaystyle \sum_{i=0}^{N-1} \frac{\hat{\mathbf{y}}_{wi} \,\hat{\mathbf{y}}_{wi} ^H}{\hat{\mathbf{y}}_{wi} ^H \,\boldsymbol{\Sigma}^{-1} \, \hat{\mathbf{y}}_{wi} } \, . \label{eqTylerbis}\tag{12}\] Existence and uniqueness of 12 require \(m < N\) (i.e., \(c_N < 1\)), which is guaranteed by Assumption 1 (specifically, \(c < 1\), which may be assumed without loss of generality in the large-dimensional regime; no constraint from Assumption 2-(v) is needed here), as well as the almost sure general position of the whitened observations \(\hat{\mathbf{y}}_{wi}\) in \(\mathbb{C}^m\); since the entries of \(\mathbf{X}\) have a density under Assumption 2-(iv), any \(m\) of them are linearly independent almost surely [42].

5.3 Eigenvalue Thresholding via the Marčenko–Pastur Law↩︎

Since \(\hat{\mathbf{S}}_w = \displaystyle\frac{1}{N}\,\mathbf{X}\,\mathbf{X}^H\), the Marčenko–Pastur theorem [45] guarantees that the Empirical Spectral Distribution (ESD) of \(\hat{\mathbf{S}}_w\) converges weakly to the Marčenko–Pastur law with ratio \(c\). Theorem 6 gives \(\left\|\check{\boldsymbol{\Sigma}}_{TYL} - \hat{\mathbf{S}}_w\right\| \overset{a.s.}{\to} 0\); since ESD convergence is continuous with respect to spectral-norm perturbations, the ESD of \(\check{\boldsymbol{\Sigma}}_{TYL}\) converges to the same Marčenko–Pastur law.

Eigenvalues of \(\check{\boldsymbol{\Sigma}}_{TYL}\) that exceed \(t=(1+\sqrt{c})^2\) are therefore attributed to signal components. Note that this threshold assumes the normalisation \(\operatorname{tr}\left(\check{\boldsymbol{\Sigma}}_{TYL}\right)=m\), consistent with 11 and 12 .

If sources are present with sufficiently high SNR, some eigenvalues of \(\check{\boldsymbol{\Sigma}}_{TYL}\) exceed \(t\). The estimated number of sources is \(\hat{p} = \#\bigl\{k : \lambda_k(\check{\boldsymbol{\Sigma}}_{TYL}) > t\bigr\}\).

5.4 Simulations↩︎

The efficiency of the whitening process is illustrated in Fig. 6. We simulate \(N=1000\) observations of dimension \(m=100\) from a correlated \(K\)-distributed process with shape parameter \(\nu=0.5\) and Toeplitz covariance matrix coefficient \(\rho=0.8\), whose \((i,j)\)-th entry is \(\rho^{|i-j|}\) for \(i,j=1,\ldots,m\). We then embed \(p=3\) informative sources with SNR \(=[5,10,15]\,\)dB into the non-Gaussian correlated noise and compare the resulting eigenvalue distributions with the Marčenko–Pastur upper bound. The eigenvalues are computed from: i)-Left the sample covariance matrix \(\hat{\mathbf{C}}_{SCM}\), ii)-Middle the Tyler \(M\)-estimator \(\hat{\mathbf{C}}_{TYL}\), and iii)-Right the Tyler \(M\)-estimator \(\check{\boldsymbol{\Sigma}}_{TYL}\) applied to whitened observations. The results show that the \(p=3\) factors (10 Monte Carlo trials) are clearly identifiable only after whitening.

The following section validates the three branches on simulated and real-world data and compares them against baseline methods.

6 Numerical Results and Comparisons↩︎

The theoretical guarantees of Sections 3–5 are validated through Monte Carlo simulations (Section 6.2), real hyperspectral images (Section 6.3), EEG source detection (Section 6.4), and portfolio optimization (Section 6.5).

6.1 Data-Driven Estimation of the Scaling Parameters↩︎

In practice, neither the textures \(\{\tau_i\}\) nor \(\xi\) are directly accessible. We describe here how both \(\xi\) and the scaling factor \(\mathbb{E}[v(\tau\,\xi)\,\tau]\)—required to form the consistent estimator \(\check{\mathbf{C}}_{FP} = \widetilde{\mathbf{C}}_{FP}/\mathbb{E}[v(\tau\,\xi)\, \tau]\) of the true scatter matrix \(\mathbf{C}\)—can be estimated from the data alone.

6.1.1 Empirical estimator of \(\xi\)↩︎

An empirical estimator of \(\xi\) based on the leave-one-out approximation [25] is used: \[\hat{\xi} = \frac{1}{N}\sum_{i=0}^{N-1} \frac{1}{m}\,\mathbf{y}_i^H\,\check{\boldsymbol{\Sigma}}^{-1}_{(i)}\, \mathbf{y}_i \,, \label{eq:xi95loo}\tag{13}\] where \(\check{\boldsymbol{\Sigma}}_{(i)}\) denotes the rank-one update of \(\check{\boldsymbol{\Sigma}}_{FP}\) with observation \(i\) removed: \[\check{\boldsymbol{\Sigma}}_{(i)} = \check{\boldsymbol{\Sigma}}_{FP} - \frac{1}{N}\,u\!\left(\frac{1}{m}\,\mathbf{y}_i^H\, \check{\boldsymbol{\Sigma}}_{FP}^{-1}\,\mathbf{y}_i\right) \mathbf{y}_i\,\mathbf{y}_i^H\, . \label{eq:Sigma95loo}\tag{14}\] It is shown in [25] that \(\xi - \hat{\xi} \overset{a.s.}{\longrightarrow} 0\).

Computing 13 naively would require inverting \(\check{\boldsymbol{\Sigma}}_{(i)}\) for each \(i\), at a cost of \(O(Nm^3)\). Applying the matrix inversion lemma to 14 instead yields the closed-form expression \[\frac{1}{m}\,\mathbf{y}_i^H\,\check{\boldsymbol{\Sigma}}_{(i)}^{-1}\, \mathbf{y}_i \;=\; t_i \;+\; \frac{u(t_i)\,m^2\,t_i^2}{N - u(t_i)\,m\,t_i}\, , \label{eq:SM95xi}\tag{15}\] where \(t_i = \displaystyle\frac{1}{m}\,\mathbf{y}_i^H\,\check{\boldsymbol{\Sigma}}_{FP}^{-1} \,\mathbf{y}_i\).

The estimator 13 therefore requires a single matrix inversion \(\check{\boldsymbol{\Sigma}}_{FP}^{-1}\), followed by \(N\) scalar corrections via 15 , at a total cost of \(O(m^3 + Nm^2)\).

6.1.2 Estimation of \(\mathbb{E}[v(\tau\,\xi)\,\tau]\)↩︎

A key observation underlying this estimator is that \(t_i\) converges almost surely to \(\tau_i\,\xi\). This can be established in four steps:

Step 1.

The surrogate matrix \(\hat{\mathbf{S}}\) defined in 6 admits the rank-one decomposition \[\hat{\mathbf{S}} = \hat{\mathbf{S}}_{(i)} + \frac{v(\tau_i\,\xi)}{N}\,\mathbf{y}_{wi}\,\mathbf{y}_{wi}^H,\] where \(\hat{\mathbf{S}}_{(i)}\) is the leave-one-out version of \(\hat{\mathbf{S}}\) with observation \(i\) removed. Applying the matrix inversion lemma to \(\hat{\mathbf{S}}^{-1}\) yields \[t_i = \frac{\dfrac{1}{m}\,\mathbf{y}_{wi}^H\, \hat{\mathbf{S}}_{(i)}^{-1}\,\mathbf{y}_{wi}}{1 + \dfrac{v(\tau_i\xi)}{N}\, \dfrac{1}{m}\,\mathbf{y}_{wi}^H\, \hat{\mathbf{S}}_{(i)}^{-1}\,\mathbf{y}_{wi}}. \label{eq:ti95SM}\tag{16}\]

Step 2.

Since observation \(i\) is absent from \(\hat{\mathbf{S}}_{(i)}\), the speckle vector \(\mathbf{x}_i\) (with \(\mathbf{y}_{wi}=\sqrt{\tau_i}\,\mathbf{x}_i\)) is independent of \(\hat{\mathbf{S}}_{(i)}\). The Hanson–Wright inequality (Lemma S-2 in the Supplementary Material; see also [37], [38]) gives, for any \(\varepsilon>0\), \[\begin{align} \mathbb{P}\Biggl( \biggl| \frac{1}{m}\,\mathbf{x}_i^{H}\,\hat{\mathbf{S}}_{(i)}^{-1}\,\mathbf{x}_i -\frac{1}{m}\operatorname{tr}\!\bigl(\hat{\mathbf{S}}_{(i)}^{-1}\bigr) \biggr| >\varepsilon \Biggr) \le 2\exp\!\Biggl( -\frac{C\,m\,\varepsilon^2}{\|\hat{\mathbf{S}}_{(i)}^{-1}\|^2} \Biggr). \end{align}\] Since \(\|\hat{\mathbf{S}}_{(i)}^{-1}\|\) is almost surely bounded, the right-hand side is summable in \(m\). The Borel–Cantelli lemma therefore implies \[\frac{1}{m}\,\mathbf{y}_{wi}^H\,\hat{\mathbf{S}}_{(i)}^{-1}\, \mathbf{y}_{wi} = \frac{\tau_i}{m}\,\mathbf{x}_i^H\,\hat{\mathbf{S}}_{(i)}^{-1} \,\mathbf{x}_i \overset{a.s.}{\longrightarrow} \frac{\tau_i}{m} \operatorname{tr}\!\left(\hat{\mathbf{S}}_{(i)}^{-1}\right). \label{eq:ti95LLN}\tag{17}\]

Step 3.

By Theorem 1 of [25], the normalized trace of \(\hat{\mathbf{S}}^{-1}_{(i)}\) converges almost surely: \[\frac{1}{m} \operatorname{tr}\!\left(\hat{\mathbf{S}}^{-1}_{(i)}\right) \overset{a.s.}{\longrightarrow} \xi, \label{eq:trace95xi}\tag{18}\] where \(\xi\) is the unique positive solution of 7 .

Step 4.

Substituting 17 and 18 into 16 , and using \(\psi(x)=x\,v(x)\), the denominator converges to \(1+c\,\psi(\tau_i\,\xi)> 0\) (since \(v\) is bounded and continuous). The two factors cancel, giving \(t_i \overset{a.s.}{\longrightarrow} \tau_i\,\xi\). Since \(t_i \overset{a.s.}{\longrightarrow} \tau_i\,\xi\), we have \(v(t_i) \overset{a.s.}{\longrightarrow} v(\tau_i\,\xi)\) and \(t_i\,v(t_i) \overset{a.s.}{\longrightarrow} \xi\,v(\tau_i\,\xi)\,\tau_i\), so that \[\widehat{\rho} \;\triangleq\; \frac{1}{\hat{\xi}}\,\frac{1}{N}\sum_{i=0}^{N-1} \psi(t_i) \;\overset{a.s.}{\longrightarrow}\; \mathbb{E}[v(\tau\,\xi)\,\tau] \triangleq \rho\,, \label{eq:Evtau95est}\tag{19}\] where \(t_i = \displaystyle\frac{1}{m}\,\mathbf{y}_i^H\,\check{\boldsymbol{\Sigma}}_{FP}^{-1} \,\mathbf{y}_i\).

6.2 Simulated Data↩︎

We evaluate the proposed approach through Monte Carlo simulations by comparing four methods: M1 (SCM-based whitening, \(\check{\boldsymbol{\Sigma}}_{SCM})\), M2 (Maronna-based whitening, \(\check{\boldsymbol{\Sigma}}_{FP}\)), M3 (direct application of Maronna’s M-estimator \(\hat{\mathbf{C}}_{FP}\) to the unwhitened observations \(\mathbf{Y}\)), and the classical AIC criterion.

Simulations use \(m=400\), \(N=2000\) (\(c=0.2\)), Toeplitz correlation \(\rho=0.7\), \(p=4\) sources with SNR \(\in[-15,20]\) dB, and CES noise with inverse-gamma or Student-\(t^2\) textures. Fig. 7 shows that M1 and M2 detect sources from \(\approx{-5}\) dB, with M2 achieving better accuracy at high SNR; both consistently outperform M3 and AIC across all SNR values.

6.3 Hyperspectral Images Application↩︎

We use Indian Pines and SalinasA (AVIRIS) and PaviaU (ROSIS) [46].

Results are in Table 2. M1 and M2 dramatically outperform M3 and AIC. For PaviaU and Cars, M1/M2 underestimate \(p\) due to (i) some sources having SNR below the threshold \(t\), and/or (ii) the empirical estimate of \(\mathbb{E}[v(\tau\,\xi)\,\tau]\) deviating from its true value due to the unknown texture distribution.

6.4 Estimation of the Model Order: EEG data↩︎

To evaluate the proposed framework in a realistic high-dimensional setting, we consider the publicly available MAMEM EEG SSVEP Dataset I [47], [48]. The dataset consists of high-density EEG recordings acquired from healthy subjects using \(m=256\) electrodes sampled at 250 Hz during steady-state visually evoked potential (SSVEP) experiments. Covariance matrices are estimated from non-overlapping 2-second windows (\(N=500\) samples), yielding a dimensional ratio \(c=m/N=0.512\), representative of modern high-dimensional EEG applications.

EEG signals exhibit spatial correlation, non-Gaussianity, and heavy-tailed fluctuations due to physiological and measurement artifacts. This motivates the use of robust covariance estimators and CES-type models, which have been shown to provide more accurate statistical descriptions of EEG observations in the presence of heavy-tailed fluctuations [49], [50].Although the Toeplitz structure only approximates EEG spatial covariance matrices, it provides an effective regularization for covariance estimation in short observation windows.

Each temporal segment is analyzed using the proposed SCM-, Maronna-, and Tyler-based estimators (Section 5), together with AIC, MDL [51], and the Kritchman–Nadler (KN) detector [52]. Table 3 reports the empirical mean and variance of the estimated model order \(\hat{p}\).

Figure 8 displays the eigenvalue distributions of the sample covariance matrix \(\hat{\mathbf{C}}_{\mathrm{SCM}}\) and the whitened covariance matrices \(\check{\mathbf{\Sigma}}_{\mathrm{SCM}}\), \(\check{\mathbf{\Sigma}}_{\mathrm{FP}}\), and \(\check{\mathbf{\Sigma}}_{\mathrm{TYL}}\). Despite the moderate sample size, the proposed asymptotic thresholds remain effective for \(c=0.512\).

For the \(K=235\) temporal windows, AIC and MDL systematically estimate a model order close to the observation dimension (\(\hat{p} \approx 255\)), while KN still yields unrealistically large estimates with high variability. These results reflect the mismatch between the assumptions underlying these methods and the correlated, non-Gaussian nature of EEG observations. In contrast, the proposed robust approaches provide stable estimates between 6.7 and 7.8 sources with variances below 2. As shown in Fig. 8, robust whitening compresses the bulk eigenvalue distribution while preserving a small number of dominant outliers, thereby improving signal/noise subspace separation. Overall, the results demonstrate the effectiveness of robust covariance whitening for model-order estimation in high-dimensional EEG data.

6.5 Estimation of the Model Order: Portfolio Optimization Application↩︎

The proposed model-order selection framework has been successfully applied to portfolio optimization problems. In particular, [53], [54] show that replacing standard covariance-based factor selection with the proposed robust approach improves both maximum diversification and global minimum variance portfolios, leading to lower out-of-sample risk and improved risk-adjusted performance.

6.5.1 Financial factor model↩︎

Let \(\mathbf{Y}=[\mathbf{y}_0,\ldots,\mathbf{y}_{N-1}]\) denote the \(m\times N\) matrix of asset returns. We assume that returns follow a \(p\)-factor model \[\mathbf{y}_t=\mathbf{M}\,\mathbf{s}_t+\sqrt{\tau_t}\,\mathbf{C}^{1/2}\,\mathbf{x}_t\, ,\] where \(\mathbf{M}\) is the factor loading matrix, \(\mathbf{s}_t\) contains the latent factors, \(\mathbf{x}_t\) is a zero-mean unitarily invariant random vector with \(\|\mathbf{x}_t\|^2=1\), \(\mathbf{C}\) is a Toeplitz-structured scatter matrix, and \(\tau_t\) models time-varying volatility. This formulation accounts for both cross-sectional dependence and heavy-tailed fluctuations commonly observed in financial returns.

Estimating the number of factors \(p\) is a central problem in portfolio construction and risk management, since an accurate separation of signal and noise subspaces directly impacts covariance estimation and portfolio allocation [55], [56].

6.5.2 Real financial data↩︎

We consider the publicly available 100 Portfolios Formed on Size and Book-to-Market dataset from the Kenneth R. French Data Library. The dataset consists of 100 value-weighted portfolios obtained from the intersection of 10 market-capitalization deciles and 10 Book-to-Market deciles. Monthly returns are available over several decades, providing a high-dimensional dataset with strong cross-sectional dependence.

The \(m=100\) portfolio return series are analyzed over rolling windows of \(N=150\) observations, corresponding to a dimensional ratio \(c=m/N=0.67\). Each window is processed using the proposed SCM-, Maronna-, and Tyler-based estimators (Section 5), together with AIC, MDL [51], and the Kritchman–Nadler (KN) detector [52]. Table 4 reports the empirical mean and variance of the estimated model order \(\hat{p}\).

Figure 9 displays the eigenvalue distributions of \(\hat{\mathbf{C}}_{\mathrm{SCM}}\), \(\check{\mathbf{\Sigma}}_{\mathrm{SCM}}\), \(\check{\mathbf{\Sigma}}_{\mathrm{FP}}\), and \(\check{\mathbf{\Sigma}}_{\mathrm{TYL}}\). The robust whitening procedures concentrate the bulk of the eigenvalues while preserving only a few dominant outliers associated with the factor structure.

In contrast with AIC, MDL, and KN, which tend to overestimate the number of factors, the proposed robust approaches consistently identify between three and five dominant factors with low variability across rolling windows. This behavior is consistent with the widely documented low-dimensional factor structure of equity returns and confirms the relevance of robust covariance whitening for model-order estimation in financial applications.

7 Conclusion↩︎

This paper has proposed and analysed a robust two-stage framework for model order selection under large-dimensional, correlated, non-Gaussian (CES) noise. Three estimation branches—SCM-based, Maronna \(M\)-estimator-based, and distribution-free Tyler-based—were developed within a unified Toeplitz whitening and RMT thresholding pipeline. Almost sure consistency of each branch was established, and an explicit closed-form eigenvalue threshold was derived, valid in the proportional-growth regime \(m/N \to c > 0\).

The proposed approach demonstrates strong performance across all experimental settings—synthetic data, real hyperspectral images, EEG recordings, and financial data—highlighting its robustness to correlation and heavy-tailed noise distributions. In contrast with existing methods, the framework requires no prior knowledge of the noise structure and remains effective in challenging high-dimensional regimes.

Beyond hyperspectral imaging, the proposed methodology is particularly relevant for high-dimensional financial data, where asset returns exhibit latent factor structures embedded in correlated and non-Gaussian noise. In this context, it provides a principled tool for robust factor selection and covariance estimation, with potential applications to EEG source detection, portfolio optimization, and risk management.

More broadly, this work opens several perspectives. On the theoretical side, extensions to non-Toeplitz structured noise (e.g., block-Toeplitz or low-rank-plus-Toeplitz covariance models) and to the spiked CES model with unknown \(p\) and unknown \(\mathbf{C}\) would be valuable. On the practical side, natural targets include radar clutter rank estimation, source localisation in antenna arrays, hyperspectral anomaly detection, and online (recursive) implementations for streaming financial data.

This Supplementary Material provides detailed proofs of Theorems 1–4 stated in the main paper. We follow the same notations as in the main paper throughout. Theorem and lemma numbers in the Supplementary Material are prefixed by “S-” to distinguish them from the main paper. Reference numbers in this document are local to the Supplementary Material.

8 Proofs of Theorem 1 and Theorem 3↩︎

The proofs follow [32]. Throughout, \(c\) denotes the limit of \(c_N=m/N\to c>0\). We use the Toeplitz spectral-norm bound from Lemma 4.1 of [57]: for a Hermitian Toeplitz matrix \(\mathbf{A}=\mathcal{L}((a_0,\ldots,a_{m-1})^T)\) with absolutely summable coefficients, \[\|\mathbf{A}\|\leq M_f =\sup_{\lambda\in[0,2\pi)}\!\left|\sum_{k=1-m}^{m-1}a_k e^{i\lambda k}\right|. \label{eq:toeplitz95bound}\tag{20}\]

S-A.1Proof of Theorem 1 (Consistency of \(\check{\mathbf{C}}_{SCM}\))↩︎

8.0.0.1 Proof outline

We bound \(\left\|\mathcal{T}\left(\displaystyle\frac{1}{N}\mathbf{Y}\mathbf{Y}^H\right)) - \mathbf{C}\right\|\) into four spectral terms—see 22 —corresponding to: (1) stochastic fluctuations of the noise, (2) bias of the noise term, (3) cross noise–signal interference, and (4) pure signal contribution. Each term is shown to vanish almost surely via, respectively: (1) Hanson–Wright concentration and Bernstein’s inequality (Lemmas 4–6), (2) direct computation, (3) bilinearity and conditional centering, and (4) absolute summability of the columns of \(\mathbf{M}\) (Assumption 3-(i) of the main paper).

8.0.0.2 Proof

Since \(\mathbb{E}[\tau]=1\) (Assumption 2), the target identity is \(\mathcal{T}\!\left(\displaystyle\frac{1}{N}\mathbf{Y}\mathbf{Y}^H\right) \overset{a.s.}{\to}\mathbf{C}\). We decompose the error into noise, cross, and signal terms and show each vanishes a.s.: the noise term via Hanson–Wright and Bernstein concentration, the cross term via bilinearity in \(\mathbf{X}\) and \(\boldsymbol{\delta}\), and the signal term via the absolute summability of \(\mathbf{M}\) (Assumption 3). As in the main body, let \(\mathbf{Y}=\mathbf{M}\,\boldsymbol{\Gamma}^{1/2}\, \boldsymbol{\delta}^H +\mathbf{C}^{1/2}\,\mathbf{X}\,\mathbf{T}^{1/2}\) and let \(\mathcal{T}\) be the Toeplitz operator defined in the introduction. Under Assumptions 1, 2 and 3, since \(\mathcal{T}\!\left(\displaystyle \frac{1}{N}\mathbf{Y}\, \mathbf{Y}^H\right)\) and \(\mathbb{E}[\tau]\,\mathbf{C}\) are both Toeplitz matrices, one can write thanks to 20 : \[\left\|\mathcal{T}\!\left(\displaystyle \frac{1}{N}\mathbf{Y}\, \mathbf{Y}^H\right) -\mathbb{E}[\tau]\,\mathbf{C}\right\| \leq \!\sup_{\lambda\in[0,2\pi)} \left|\hat{\gamma}_m(\lambda) -\mathbb{E}[\tau]\,\gamma_m(\lambda)\right|\, , \label{eq:thm195main}\tag{21}\] where \(\gamma_m(\lambda)=\displaystyle\sum_{k=1-m}^{m-1}c_{k}\,e^{ik\lambda}\), \(\hat{\gamma}_m(\lambda)=\displaystyle\sum_{k=1-m}^{m-1}\check{c}_{k}\,e^{ik\lambda}\) with \(c_{-k}=c_k^\star\) and \(\check{c}_{-k}=\check{c}_k^\star\).
The following lemma is essential for the development of the proof.

Decomposing \(\hat{\gamma}_m(\lambda)\) according to the three contributions of \(\mathbf{Y}=\mathbf{M}\,\mathbf{S}+\mathbf{C}^{1/2}\,\mathbf{X}\,\mathbf{T}^{1/2}\), one defines \(\hat{\gamma}_m(\lambda) = \hat{\gamma}^{\mathrm{sign}}_m(\lambda) + \hat{\gamma}^{\mathrm{cross}}_m(\lambda) + \hat{\gamma}^{\mathrm{noise}}_m(\lambda)\) with \[\begin{align} \hat{\gamma}^{\mathrm{sign}}_m(\lambda) &= \frac{1}{N}\,\mathbf{d}_m^H(\lambda)\,\mathbf{M}\, \boldsymbol{\Gamma}^{1/2}\,\boldsymbol{\delta}^H\,\boldsymbol{\delta}\, \boldsymbol{\Gamma}^{1/2}\,\mathbf{M}^H\,\mathbf{d}_m(\lambda)\,,\\ \hat{\gamma}^{\mathrm{cross}}_m(\lambda) &= \frac{2}{N}\,\mathcal{R}e\!\left(\mathbf{d}_m^H(\lambda) \,\mathbf{C}^{1/2}\,\mathbf{X}\,\mathbf{T}^{1/2}\,\boldsymbol{\delta}\, \boldsymbol{\Gamma}^{1/2}\,\mathbf{M}^H\,\mathbf{d}_m(\lambda)\right),\\ \hat{\gamma}^{\mathrm{noise}}_m(\lambda) &= \frac{1}{N}\,\mathbf{d}_m^H(\lambda)\, \mathbf{C}^{1/2}\,\mathbf{X}\,\mathbf{T}\,\mathbf{X}^H\,\mathbf{C}^{1/2} \,\mathbf{d}_m(\lambda)\,. \end{align}\]

Equation 21 then leads to: \[\begin{align} \left\|\mathcal{T}\!\left(\frac{1}{N}\mathbf{Y}\mathbf{Y}^H\right) -\mathbb{E}[\tau]\,\mathbf{C}\right\| &\le \underbrace{ \sup_\lambda\left| \hat{\gamma}^{\mathrm{noise}}_m(\lambda) -\mathbb{E}\!\left[ \hat{\gamma}^{\mathrm{noise}}_m(\lambda) \right] \right| }_{\text{Term 1}} \nonumber\\ & +\underbrace{ \sup_\lambda\left| \mathbb{E}\!\left[ \hat{\gamma}^{\mathrm{noise}}_m(\lambda) \right] -\mathbb{E}[\tau]\,\gamma_m(\lambda) \right| }_{\text{Term 2}} \nonumber\\ & +\underbrace{ \sup_\lambda\left| \hat{\gamma}^{\mathrm{cross}}_m(\lambda) \right| }_{\text{Term 3}} +\underbrace{ \sup_\lambda\left| \hat{\gamma}^{\mathrm{sign}}_m(\lambda) \right| }_{\text{Term 4}}. \label{eq:thm195four} \end{align}\tag{22}\] We now analyse each term of 22 and show each term tends to zero almost surely.

Term 1: \(\displaystyle\sup_\lambda\left|\hat{\gamma}^{\mathrm{noise}}_m(\lambda) -\mathbb{E}[\hat{\gamma}^{\mathrm{noise}}_m(\lambda)]\right|\). We introduce auxiliary lemmas.

Let \(\lfloor\cdot\rfloor\) be the floor function. Choosing \(\beta>3/2\), \(I=\{0,\ldots,\lfloor N^\beta\rfloor-1\}\), and \(\lambda_i=2\pi i/\lfloor N^\beta\rfloor\) for \(i\in I\), one obtains: \[\sup_{\lambda\in[0,2\pi)} \left| \hat{\gamma}^{\mathrm{noise}}_m(\lambda) -\mathbb{E}\!\left[\hat{\gamma}^{\mathrm{noise}}_m(\lambda)\right] \right| \leq \chi_1+\chi_2+\chi_3 \,, \label{eq:chi95decomp}\tag{24}\] where \[\begin{align} \chi_1 &\triangleq \max_{i\in I}\sup_{\lambda\in[\lambda_i,\lambda_{i+1}]} \left| \hat{\gamma}^{\mathrm{noise}}_m(\lambda) -\hat{\gamma}^{\mathrm{noise}}_m(\lambda_i) \right|, \\[0.5em] \chi_2 &\triangleq \max_{i\in I} \left| \hat{\gamma}^{\mathrm{noise}}_m(\lambda_i) -\mathbb{E}\!\left[ \hat{\gamma}^{\mathrm{noise}}_m(\lambda_i) \right] \right|, \\[0.5em] \chi_3 &\triangleq \max_{i\in I}\sup_{\lambda\in[\lambda_i,\lambda_{i+1}]} \left| \mathbb{E}\!\left[ \hat{\gamma}^{\mathrm{noise}}_m(\lambda_i) \right] -\mathbb{E}\!\left[ \hat{\gamma}^{\mathrm{noise}}_m(\lambda) \right] \right|. \end{align}\]

The convergence of each term is established through the following lemmas.

These three inequalities prove that, Term 1, \(\forall\, x\in\mathbb{R}^{+\star}\) tends to zero almost surely: \[P\Bigl( \sup_{\lambda\in[0,2\pi)} \bigl| \hat{\gamma}^{\mathrm{noise}}_m(\lambda) -\mathbb{E}[\hat{\gamma}^{\mathrm{noise}}_m(\lambda)] \bigr| >x \Bigr) \xrightarrow[N\to\infty]{} 0 .\]

While Term 1 describes the purely stochastic fluctuations of the noise, the second term, Term 2, captures the interaction between the signal signatures and the estimated noise subspace.

Term 2. By definition, \(\hat{\gamma}^{\mathrm{noise}}_m(\lambda)\) involves only the noise component \(\mathbf{C}^{1/2}\mathbf{X}\mathbf{T}^{1/2}\) of \(\mathbf{Y}\). Since \(\{\tau_j\}\) and \(\{\mathbf{x}_j\}\) are independent and \(\mathbb{E}[\mathbf{x}_j\mathbf{x}_j^H]=\mathbf{I}_m\), a direct computation gives: \[\begin{align} \mathbb{E}\!\left[\hat{\gamma}^{\mathrm{noise}}_m(\lambda)\right] &\!= \!\frac{1}{N}\sum_{j=0}^{N-1}\! \mathbb{E}[\tau_j] \mathbf{d}_m^H(\lambda)\mathbf{C}^{1/2} \mathbb{E}\!\left[\mathbf{x}_j\mathbf{x}_j^H\right] \mathbf{C}^{1/2}\mathbf{d}_m(\lambda), \\ &= \mathbb{E}[\tau]\, \mathbf{d}_m^H(\lambda)\,\mathbf{C}\,\mathbf{d}_m(\lambda) = \mathbb{E}[\tau]\,\gamma_m(\lambda)\,, \end{align}\] where we used \(\mathbb{E}[\tau_j]=\mathbb{E}[\tau]=1\) (Assumption 2-(iii)) and \(\mathbb{E}[\mathbf{x}_j\mathbf{x}_j^H]=\mathbf{I}_m\) (Assumption 2-(iv)). Hence Term 2 vanishes identically.

The analysis of Term 3 is more involved as it represents the residual interference after whitening. Its convergence is established by leveraging the concentration of quadratic forms and the properties of the Toeplitz rectification operator.

Term 3: \(\displaystyle\sup_\lambda|\hat{\gamma}^{\mathrm{cross}}_m(\lambda)|\).

Recall: \[\hat{\gamma}^{\mathrm{cross}}_m(\lambda) = \frac{2}{N}\,\mathcal{R}e\!\left(\mathbf{d}_m^H(\lambda)\, \mathbf{C}^{1/2}\,\mathbf{X}\,\mathbf{T}^{1/2}\,\boldsymbol{\delta}\, \boldsymbol{\Gamma}^{1/2}\,\mathbf{M}^H\,\mathbf{d}_m(\lambda)\right)\,.\] Define \(\mathbf{H}(\lambda) \;=\; \mathbf{M}^H\, \mathbf{d}_m(\lambda)\,\mathbf{d}_m^H(\lambda)\, \mathbf{C}^{1/2} \in \mathbb{C}^{p \times m}\) and \(\mathbf{K} \;=\; \mathbf{T}^{1/2}\,\boldsymbol{\Gamma}^{1/2} \in \mathbb{C}^{N \times p}\). Then: \[\begin{align} \hat{\gamma}^{\mathrm{cross}}_m(\lambda) &= \frac{2}{N}\, \mathrm{Re}\!\left[ \operatorname{tr}\!\left( \mathbf{H}(\lambda)\, \mathbf{X}\, \mathbf{K}\, \boldsymbol{\delta}^H \right) \right]\, , \nonumber\\ &= \frac{2}{N}\, \mathrm{Re}\!\left[ \mathrm{vec}^{T}\!\left( \mathbf{H}(\lambda)\,\mathbf{X} \right)\, \mathrm{vec}\!\left( \mathbf{K}\,\boldsymbol{\delta} \right) \right]. \label{eq:cross95compact} \end{align}\tag{29}\] The key point is that this expression is bilinear in the two independent random matrices \(\mathbf{X}\) and \(\boldsymbol{\delta}\); combined with \(\mathbb{E}[\mathbf{X}]=\mathbf{0}\), this drives convergence to zero.

We decompose the proof into two parts: a Lipschitz reduction to a finite grid (control of \(\chi_1\)), and a pointwise then uniform almost sure convergence over the grid (control of \(\chi_2\)).

Reduction to a finite grid.
Let \(\beta > 3/2\), \(I = \{0,\ldots,\lfloor N^\beta\rfloor-1\}\), and \(\lambda_i = 2\pi i/\lfloor N^\beta\rfloor\). Write: \[\sup_\lambda \left|\hat{\gamma}^{\mathrm{cross}}_m(\lambda)\right| \le \rho_1 + \rho_2\, ,\] where \[\begin{align} &\rho_1 = \max_{i\in I} \sup_{\lambda\in[\lambda_i,\lambda_{i+1}]} \left| \hat{\gamma}^{\mathrm{cross}}_m(\lambda) -\hat{\gamma}^{\mathrm{cross}}_m(\lambda_i) \right|\, , \nonumber\\ &\rho_2 = \max_{i\in I} \left| \hat{\gamma}^{\mathrm{cross}}_m(\lambda_i) \right|\, . \label{eq:cross95grid} \end{align}\tag{30}\]

Combining Lemmas 7 and 8 in 30 : \[\sup_{\lambda\in[0,2\pi)}\left|\hat{\gamma}^{\mathrm{cross}}_m(\lambda)\right| \leq \rho_1+\rho_2\xrightarrow{a.s.}0\,,\] completing the proof that Term 3 tends to zero almost surely.

Term 4: \(\displaystyle\sup_\lambda|\hat{\gamma}^{\mathrm{sign}}_m(\lambda)|\). We recall that: \[\hat{\gamma}^{\mathrm{sign}}_m(\lambda) = \frac{1}{N}\,\mathbf{d}_m^H(\lambda)\,\mathbf{M}\, \boldsymbol{\Gamma}^{1/2}\,\boldsymbol{\delta}^H\,\boldsymbol{\delta}\, \boldsymbol{\Gamma}^{1/2}\,\mathbf{M}^H\,\mathbf{d}_m(\lambda)\, .\] where \(\mathbf{M}\in\mathbb{C}^{m\times p}\), \(\boldsymbol{\delta}\in\mathbb{C}^{N\times p}\), \(\boldsymbol{\Gamma}\in\mathbb{C}^{p\times p}\) is Hermitian nonnegative definite, and \(\mathbf{d}_m(\lambda)\in\mathbb{C}^m\) with \(\|\mathbf{d}_m(\lambda)\|=1\).

The key structural difference with the cross term is that \(\hat{\gamma}^{\mathrm{sign}}_m(\lambda)\) is quadratic in \(\boldsymbol{\delta}\) (and not bilinear in \(\mathbf{X}\) and \(\boldsymbol{\delta}\)). There is therefore only a single source of randomness (\(\boldsymbol{\delta}\)), and the reasoning is more straightforward but different.

Step 1: reduction to a sum of i.i.d.complex products. Define the deterministic vectors: \(\mathbf{h}(\lambda) \triangleq \mathbf{M}^H\mathbf{d}_m(\lambda) \in\mathbb{C}^p\) and \(\mathbf{w}(\lambda) \triangleq \boldsymbol{\Gamma}\,\mathbf{h}(\lambda) = \boldsymbol{\Gamma}\,\mathbf{M}^H\,\mathbf{d}_m(\lambda)\in\mathbb{C}^p\), satisfying \(\|\mathbf{h}(\lambda)\|\leq\|\mathbf{M}\|<\infty\) and \(\|\mathbf{w}(\lambda)\|\leq\|\boldsymbol{\Gamma}\|\,\|\mathbf{M}\| \triangleq K < \infty\) uniformly in \(\lambda\) (Assumption 3).

Denote by \(\boldsymbol{\delta}_j\in\mathbb{C}^p\) the \(j\)-th row of \(\boldsymbol{\delta}\), viewed as a column vector. The \(\{\boldsymbol{\delta}_j\}_{j=0}^{N-1}\) are i.i.d. \(\mathcal{CN}(\mathbf{0},\mathbf{I}_p)\). Define the scalar random variables \(w_j(\lambda) \triangleq \boldsymbol{\delta}_j^T\,\mathbf{w}(\lambda) = \displaystyle\sum_{i=1}^p \delta_{ji}\,w_i(\lambda)\) and \(z_j(\lambda) \triangleq \boldsymbol{\delta}_j^T\,\mathbf{h}(\lambda) = \displaystyle\sum_{i=1}^p \delta_{ji}\,h_i(\lambda)\). These are two distinct linear forms of the same vector \(\boldsymbol{\delta}_j\), hence correlated. Decomposing \(\boldsymbol{\delta}^H\,\boldsymbol{\Gamma}\, \boldsymbol{\delta} = \displaystyle\sum_{j=0}^{N-1}\boldsymbol{\delta}_j^*\,\boldsymbol{\delta}_j^T\, \boldsymbol{\Gamma}\) and applying to \(\mathbf{h}(\lambda)\): \[\hat{\gamma}^{\mathrm{sign}}_m(\lambda) = \frac{1}{N}\,\mathbf{h}^H(\lambda)\,\boldsymbol{\delta}^H\, \boldsymbol{\Gamma}\,\boldsymbol{\delta}\,\mathbf{h}(\lambda) = \frac{1}{N}\sum_{j=0}^{N-1} z_j^\ast(\lambda)\,w_j(\lambda)\,. \label{eq:sign95product}\tag{31}\] Note that \(\hat{\gamma}^{\mathrm{sign}}_m(\lambda)\in\mathbb{R}\) since it is a diagonal entry of the Hermitian matrix \(\displaystyle\frac{1}{N}\,\mathbf{d}_m^H(\lambda)\, \mathbf{Y}\,\mathbf{Y}^H\, \mathbf{d}_m(\lambda)\).

Having expressed \(\hat{\gamma}^{\mathrm{sign}}_m(\lambda)\) as an average of i.i.d. complex products 31 , we next compute their joint distribution, which will reveal the mean and provide the sub-exponential bound needed for concentration.

Step 2: joint distribution of \((z_j, w_j)^T\). Since \(z_j\) and \(w_j\) are both linear forms of \(\boldsymbol{\delta}_j\sim\mathcal{CN}(\mathbf{0},\mathbf{I}_p)\): \[\begin{pmatrix}z_j(\lambda)\\ w_j(\lambda)\end{pmatrix} = \begin{pmatrix}\mathbf{h}^T(\lambda)\\ \mathbf{w}^T(\lambda)\end{pmatrix} \boldsymbol{\delta}_j \sim \mathcal{CN}\!\left(\mathbf{0},\,\boldsymbol{\Sigma}_{zw}\right),\] where the joint covariance matrix is: \[\boldsymbol{\Sigma}_{zw} = \begin{pmatrix} \|\mathbf{h}(\lambda)\|^2 & \mathbf{h}^H(\lambda)\,\boldsymbol{\Gamma}\,\mathbf{h}(\lambda) \\ \mathbf{h}^H(\lambda)\,\boldsymbol{\Gamma}^H\,\mathbf{h}(\lambda) & \|\mathbf{w}(\lambda)\|^2 \end{pmatrix}\,.\] In particular, the cross-covariance satisfies: \[\begin{align} \mathbb{E}\!\left[z_j^\ast(\lambda)\,w_j(\lambda)\right] &= \mathbf{h}^H(\lambda)\, \boldsymbol{\Gamma}\, \mathbf{h}(\lambda) \nonumber\\ &= \mathbf{d}_m^H(\lambda)\, \mathbf{M}\, \boldsymbol{\Gamma}\, \mathbf{M}^H\, \mathbf{d}_m(\lambda). \label{eq:sign95mean} \end{align}\tag{32}\] We now show this expectation vanishes as \(m, N\to\infty\).

Step 3: the mean tends to zero. By 32 , we have: \(\mathbb{E}\!\left[\hat{\gamma}^{\mathrm{sign}}_m(\lambda)\right] = \mathbf{d}_m^H(\lambda)\,\mathbf{M}\,\boldsymbol{\Gamma}\, \mathbf{M}^H\,\mathbf{d}_m(\lambda)\). Since the columns of \(\mathbf{M}\) are absolutely summable (Assumption 3), for each \(j \in \{1,\ldots,p\}\), we can rewrite \(\left|[\mathbf{M}^H\,\mathbf{d}_m(\lambda)]_j\right|^2\) as: \[\frac{1}{m}\left|\sum_{k=0}^{m-1}[\mathbf{M}]_{k,j}\,e^{-ik\lambda}\right|^2\, \leq \frac{1}{m}\left(\sum_{k=0}^{m-1}|[\mathbf{M}]_{k,j}|\right)^2 \xrightarrow{m\to\infty} 0\,,\] since \(\displaystyle\sum_k |[\mathbf{M}]_{k,j}| < \infty\) by Assumption 3-(i), so \(\bigl(\displaystyle\sum_k |[\mathbf{M}]_{k,j}|\bigr)^2\) is a finite constant independent of \(m\). Therefore, \(\left|[\mathbf{M}^H\,\mathbf{d}_m(\lambda)]_j\right|^2 \to 0\) uniformly in \(\lambda\), and: \[\left|\mathbb{E}\!\left[\hat{\gamma}^{\mathrm{sign}}_m(\lambda)\right]\right| \leq \|\boldsymbol{\Gamma}\|\,\|\mathbf{M}^H\, \mathbf{d}_m(\lambda)\|^2 \xrightarrow{N\to\infty} 0\,,\] uniformly in \(\lambda\). The mean, therefore, tends to zero, but we must also control the fluctuations \(\displaystyle\frac{1}{N}\sum_j \xi_j(\lambda)\) around this vanishing mean. The key insight is that each \(\xi_j(\lambda)\) is sub-exponential (as a centered product of two sub-Gaussian variables), enabling a Bernstein bound.

Step 4: Sub-exponential concentration. Recall \(\|X\|_{\psi_r}=\inf\{t>0:\mathbb{E}[e^{|X|^r/t^r}]\leq 2\}\) for \(r\in\{1,2\}\); \(X\) is sub-Gaussian iff \(\|X\|_{\psi_2}<\infty\).

Let \(\xi_j(\lambda) \triangleq z_j^\ast(\lambda)\,w_j(\lambda) - \mathbb{E}[z_j(\lambda)^\ast \, w_j(\lambda)]\). Since \(z_j\), \(w_j\) are Gaussian linear forms (Step 2), they are sub-Gaussian with \(\|z_j(\lambda)\|_{\psi_2}\leq C\|\mathbf{M}\|\) and \(\|w_j(\lambda)\|_{\psi_2}\leq C\,\|\boldsymbol{\Gamma}\|\|\mathbf{M}\|\) uniformly in \(\lambda\) [37]. Their product \(\xi_j(\lambda)\) is sub-exponential [37]: \[\begin{align} \|\xi_j(\lambda)\|_{\psi_1} &\le b \triangleq C\,\|\mathbf{M}\|^2\, \|\boldsymbol{\Gamma}\|, \\ \mathbb{E}\!\left[|\xi_j(\lambda)|^2\right] &\le \nu^2 \triangleq C\,\|\mathbf{M}\|^4\, \|\boldsymbol{\Gamma}\|^2. \end{align}\] uniformly in \(j\) and \(\lambda\). Bernstein’s inequality [37] then gives for any \(t>0\): \[\label{eq:sign95bernstein95final} P\!\left(\left|\frac{1}{N}\sum_{j=0}^{N-1}\xi_j(\lambda)\right|>t\right) \leq 2\exp\!\left(-CN\min\!\left(\frac{t^2}{\nu^2},\frac{t}{b}\right)\right)\,,\tag{33}\] which is summable in \(N\) for any fixed \(t>0\). Recalling that \(\hat{\gamma}^{\mathrm{sign}}_m(\lambda) = \displaystyle \frac{1}{N}\sum_{j=0}^{N-1} \xi_j(\lambda) + \mathbb{E}\!\left[\hat{\gamma}^{\mathrm{sign}}_m(\lambda)\right]\), and using Step 3, which ensures that \(\mathbb{E}\!\left[\hat{\gamma}^{\mathrm{sign}}_m(\lambda)\right] \longrightarrow 0\), we obtain that, for any fixed \(\lambda\), \(\hat{\gamma}^{\mathrm{sign}}_m(\lambda) \xrightarrow{a.s.} 0\) by the Borel–Cantelli lemma.

Pointwise almost sure convergence is established. We extend it uniformly to \(\lambda\in[0,2\pi)\) via a discretisation argument analogous to those used above.

Step 5: uniform control over \([0,2\pi)\). For \(|\lambda-\lambda'|\leq 2\pi/\lfloor N^\beta\rfloor\) and by setting \(A_j = |z_j^\ast(\lambda)| \,|w_j(\lambda)-w_j(\lambda')|\) and \(B =|z_j^\ast(\lambda)-z_j^\ast(\lambda')| \,|w_j(\lambda')|\), we can rewrite \(\left| \hat{\gamma}^{\mathrm{sign}}_m(\lambda) - \hat{\gamma}^{\mathrm{sign}}_m(\lambda') \right|\) as \[\Bigl| \frac{1}{N} \sum_{j=0}^{N-1} \bigl( z_j^\ast(\lambda)w_j(\lambda) - z_j^\ast(\lambda')w_j(\lambda') \bigr) \Bigr| \le \frac{1}{N} \sum_{j=0}^{N-1}(A_j+B_j)\, .\] Since \(z_j(\lambda)-z_j(\lambda') = \boldsymbol{\delta}_j^T\,\mathbf{M}^H \left(\mathbf{d}_m(\lambda) -\mathbf{d}_m(\lambda')\right)\) and similarly for \(w_j\), and since \(\lambda\mapsto\mathbf{d}_m(\lambda)\) is Lipschitz with constant \((m-1)/\sqrt{m}\): \[\begin{align} &\max\!\left( \left|z_j(\lambda)-z_j(\lambda')|, |w_j(\lambda)-w_j(\lambda')\right| \right) \\ &\qquad\le K\, \|\boldsymbol{\delta}_j\| \|\mathbf{M}\| \frac{m-1}{\sqrt m} \,|\lambda-\lambda'|\,. \end{align}\] Bounding \(|z_j(\lambda)|\leq\|\boldsymbol{\delta}_j\|\|\mathbf{M}\|\) and \(|w_j(\lambda)|\leq\|\boldsymbol{\delta}_j\|K\), one obtains: \[\left|\hat{\gamma}^{\mathrm{sign}}_m(\lambda) -\hat{\gamma}^{\mathrm{sign}}_m(\lambda')\right| \leq \frac{C\,\|\mathbf{M}\|^2 K}{N} \,\frac{(m-1)}{\sqrt{m}\,\lfloor N^\beta\rfloor} \,\frac{1}{N}\sum_{j=0}^{N-1}\|\boldsymbol{\delta}_j\|^2\,.\] Since \(\sum_j\|\boldsymbol{\delta}_j\|^2/N\overset{a.s.}{\to}p\) (Law of Large Numbers, each \(\|\boldsymbol{\delta}_j\|^2\) has mean \(p\)) and \((m-1)/(\sqrt{m}\lfloor N^\beta\rfloor)=O(N^{1/2-\beta})\to 0\) for \(\beta>1/2\), \(\displaystyle \max_{i\in I}\sup_{\lambda\in[\lambda_i,\lambda_{i+1})} \left|\hat{\gamma}^{\mathrm{sign}}_m(\lambda) -\hat{\gamma}^{\mathrm{sign}}_m(\lambda_i)\right| \overset{a.s.}{\longrightarrow} 0\). Combining with the pointwise almost sure convergence of Step 4 and a union bound over the \(\lfloor N^\beta\rfloor\) grid points (valid since the bound in 33 is summable in \(N\)), we obtain \(\displaystyle \sup_{\lambda\in[0,2\pi)}\hat{\gamma}^{\mathrm{sign}}_m(\lambda) \overset{a.s.}{\longrightarrow} 0\).

S-A.2Proof of Theorem 3 (Consistency of \(\widetilde{\mathbf{C}}_{FP}\) and \(\check{\mathbf{C}}_{FP}\))↩︎

The proof follows the same idea. Recall that \(\widetilde{\mathbf{C}}_{FP}=\mathcal{T} \left(\hat{\mathbf{C}}_{FP}\right)\) and \(\check{\mathbf{C}}_{FP}=\widetilde{\mathbf{C}}_{FP}/ \mathbb{E}[v(\tau\,\xi)\,\tau]\). The equation to prove is \(\left\|\widetilde{\mathbf{C}}_{FP} -\mathbb{E}[v(\tau\,\xi)\,\tau]\,\mathbf{C}\right\| \xrightarrow{a.s.}0\), from which \(\|\check{\mathbf{C}}_{FP} - \mathbf{C}\| \xrightarrow{a.s.}0\) follows immediately by dividing by the positive constant \(\mathbb{E}[v(\tau\,\xi)\,\tau]\). This splits as: \[\begin{align} \left\|\mathcal{T}\left(\hat{\mathbf{C}}_{FP}\right) -\mathbb{E}[v(\tau\,\xi)\,\tau]\,\mathbf{C}\right\| & \leq \left\|\mathcal{T}\left(\hat{\mathbf{C}}_{FP}-\hat{\mathbf{S}}\right)\right\| \nonumber \\ & +\left\|\mathcal{T}\left(\hat{\mathbf{S}}\right) -\mathbb{E}[v(\tau\,\xi)\,\tau]\,\mathbf{C}\right\| \, , \nonumber \\ &\triangleq \, \eta_1+\eta_2\, , \label{eq:thm395split} \end{align}\tag{34}\] where \(\hat{\mathbf{S}}=\displaystyle \frac{1}{N}\sum_{i=0}^{N-1} v(\tau_i\, \xi)\,\mathbf{y}_{wi}\, \mathbf{y}_{wi}^H\).

Part 1: convergence of \(\eta_1=\left\|\mathcal{T}\left(\hat{\mathbf{C}}_{FP} -\hat{\mathbf{S}}\right)\right\|\).

It is proven in [25] that \(\left\|\hat{\mathbf{C}}_{FP}-\hat{\mathbf{S}}\right\|\xrightarrow{a.s.}0\). Since the Toeplitz rectification satisfies \(\|\mathcal{T}(\mathbf{A})\| \leq \displaystyle\sup_\lambda |\mathbf{d}_m^H(\lambda)\mathbf{A} \mathbf{d}_m(\lambda)| \leq \|\mathbf{A}\|\) for any Hermitian matrix \(\mathbf{A}\), we have: \[\begin{align} \left\|\mathcal{T}\left(\hat{\mathbf{C}}_{FP}-\hat{\mathbf{S}}\right)\right\| &\leq\sup_{\lambda\in[0,2\pi)}\left| \mathbf{d}_m^H(\lambda)\,\left(\hat{\mathbf{C}}_{FP}-\hat{\mathbf{S}}\right)\, \mathbf{d}_m(\lambda)\right| \, , \\ & \leq\left\|\hat{\mathbf{C}}_{FP}-\hat{\mathbf{S}}\right\| \xrightarrow{a.s.}0\, . \end{align}\]

Part 2: convergence of \(\eta_2 = \left\|\mathcal{T}\!\left(\hat{\mathbf{S}}\right) - \mathbb{E}[v(\tau\,\xi)\,\tau]\,\mathbf{C}\right\|\).

Since \(\eta_1,\eta_2\xrightarrow{a.s.}0\), Theorem 3 follows.0◻

9 Proofs of Theorems 2 and 4↩︎