Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
113,744 characters · 29 sections · 95 citation commands
How weak are weak factors? Uniform inference for signal strength in signal plus noise models
\address[Anna Bykhovskaya]{Duke University} \email{[email removed]}
\address[Vadim Gorin]{University of California at Berkeley} \email{[email removed]}
\address[Sasha Sodin]{The Hebrew University of Jerusalem and Queen Mary University of London} \email{[email removed]}
In the modern era researchers increasingly have access to high-dimensional data across a wide range of fields. These data are inevitably contaminated by various forms of error and noise, making the separation of meaningful structure from background noise a central challenge. To address this analysts commonly employ dimension-reduction techniques. The two dominant approaches are low-rank methods, which assume that the underlying signal lies in a lower-dimensional subspace, and sparsity-based methods, which assume that only a small subset of variables or parameters are truly relevant, i.e., nonzero. This paper adopts the low-rank perspective. For a discussion of settings where this assumption is appropriate, we refer to udell2019big, giannone2021economic, and thibeault2024low. In particular, giannone2021economic argue that numerous data sets in macroeconomics, microeconomics, and finance exhibit dense, rather than sparse, structures.
A prototypical example of a low rank setting is the factor model, where one observes an $N\times S$ data matrix $X$ and assumes that it can be decomposed as
where $F$ is an $S\times r$ matrix of factors, $L$ is an $N\times r$ matrix of factor loadings, and $LF^\mathsf T$ represents the low-rank signal of interest. The signal rank $r$ is small relative to the large dimensions $N$ and $S$. The remainder $\mathcal{E}$ is a noise matrix, often assumed to have i.i.d. mean-zero entries in the simplest setting.
The feasibility of consistently estimating the signal component $LF^\mathsf T$ from the observed data $X$ hinges on the strength of the signal, which can be quantified by the singular values of $LF^\mathsf T$. When these singular values are large, the signal is strong and estimation is reliable, as can be directly predicted from the form of (ref). As the signal weakens, the data $X$ becomes less informative, and below a certain critical threshold, accurate recovery becomes impossible. The relationships between the strength of the signal and feasibility of reconstruction of $LF^\mathsf T$ have been rigorously analyzed in a number of studies, see, e.g., stock2002forecasting,bai2002determining,bai2003inferential,paul2007asymptotics,onatski2012asymptotics,johnstone2018pca,bai2023approximate,fan2024can,barigozzi2024dynamic and references therein.
Given this behavior, applied work using factor models should begin by assessing the strength of the factors, since the validity of any inference on $L$ or $F$ critically depends on it. However, in practice, this step is often overlooked\footnote{This pattern is evident in the vast majority of approximately 120 papers that employ factor models or PCA-related techniques, published in the five leading economics journals between 2015 and 2025.}, and most studies tacitly assume that the factors are strong, without conducting any formal diagnostics.
This paper seeks to emphasize the importance of assessing signal strength in a broad class of “signal plus noise” models. To that end, we develop novel procedures for constructing confidence intervals for signal strength. Crucially, we do not assume that the signals are strong -- an assumption often unjustified in empirical applications. Instead, our analysis remains valid across the full range of regimes: strong, weak, and critical signals.
We analyze four classical high-dimensional statistical models: the factor model ((ref)), spiked sample covariance, the spiked Wigner model, and spiked canonical correlations. These correspond to three fundamental ensembles from random matrix theory: the Laguerre/Wishart ensemble for the first two, the Hermite/Gaussian/Wigner ensemble for the third, and the Jacobi ensemble for the fourth. Each model can be viewed as an instance of the signal plus noise framework, also known as spiked random matrices, a term originating with johnstone2001distribution, in which a low-rank signal matrix is embedded in a high-dimensional noisy environment. The goal is to detect and quantify the signal.
The signal in each model can be decomposed into a sum of rank-one components. Each component is characterized by a positive scalar (its strength) and one or two unit-norm vectors (its direction), depending on the setup. In this work we focus solely on the signal strength and do not consider inference on directions.
Our analysis is based on spectral methods, whereby signal strength is inferred from the eigenvalues of certain model-specific matrices. In all four setups a well-documented phase transition phenomenon arises: the signal strength can be consistently estimated (in the high-dimensional asymptotic regime with proportional growth of data dimensions) only when it exceeds a critical threshold, see jones1978eigenvalue, baik2006eigenvalues, onatski2012asymptotics, bao2019canonical and more references in Section (ref). When the signal strength falls below the threshold, only partial probabilistic information, such as asymptotics of the likelihood ratio test can be recovered, but reliable point estimation becomes impossible, see, e.g.\ onatski2013asymptotic,onatski2014signal,dobriban2017sharp,johnstone2020testing,el2020fundamental. The intermediate regime, where the signal strength is close to the threshold, is typically referred to as the “critical” regime. This regime is particularly challenging for inference.
In the super-critical case, where the strength is significantly above the threshold, the estimation procedure is quite straightforward: one takes the largest eigenvalue, applies to it a certain explicit function (see Section (ref) for the formulas) and gets the strength of the strongest signal. Repeating the same with the second, third, etc., eigenvalues one gets strengths of the further components of the signal and the only question is when to stop, i.e., after which step one should declare that the following signals are too weak and can not be recovered. There are many results in the literature proposing various algorithms to choose the stopping point. We further remark that for very strong signals the function one should apply to the eigenvalues is close to identity ($f(x)=x$), whereas for weaker signals the function exhibits stronger dependence on the model of interest.
Once point estimates of the signal strengths are obtained, the next natural question is how to quantify uncertainty -- specifically, how to construct confidence intervals for these estimates. The existing literature offers little guidance on this front -- particularly guidance that is consistent across models and signal strengths -- with most results focusing on strong signals. The technical challenge is rooted in the nonstandard asymptotic behavior of the eigenvalues near the phase transition threshold. While the fluctuations of the top eigenvalues are asymptotically Gaussian for well-separated (super-critical) signals, the limiting distribution becomes highly non-Gaussian and analytically intricate as the signal strength approaches the critical boundary (see baik2005phase, mo2012rank, and bloemendal2013limits for rigorous results in the sample covariance setting).
Our paper fills this gap by proposing a general procedure for constructing confidence intervals for signal strength. Remarkably, across all four models we study, the confidence intervals are characterized by a common limiting (stochastic) object, which we call the Airy--Green function and denote $\mathcal{G}(w)$. Our main contributions are: a rigorous construction of this function, a unified set of theorems linking it to the four canonical models, and tabulated confidence intervals based on $\mathcal{G}(w)$. The only model-specific components are a set of scaling constants, which we provide explicitly for each setting. In addition, our results imply a formula-free, bootstrap-type procedure for constructing confidence intervals.
In economics and finance it has long been observed that many data sets contain factors that are either non-informative or far from strong -- see e.g.\ giglio2023prediction and kim2024testing for overviews and extensive references. This concern is especially apparent in the vast “factor zoo” of potential variables proposed to explain stock returns. This empirical reality has motivated a line of theoretical research focused on inference for weaker factors. Broadly speaking, factors can be classified by their strength into three categories: strong (as in, e.g.\ bai2002determining,stock2002forecasting), semi-strong (as in, e.g.\ bai2023approximate,fan2024can), and weak\footnote{What we call semi-strong factors are sometimes referred to as weak, while weak factors may be termed weakly influential or extremely weak.} (as in, e.g., onatski2012asymptotics). The literature also includes statistical procedures for testing and distinguishing between these types of factors (see, in particular, kim2024testing). Over the past decades a growing body of research has focused specifically on factor strength, including contributions by chudik2011weak, bailey2016exponent, wang2017asymptotics, lettau2020estimating, cai2020limiting, bailey2021measurement, freyaldenhoven2022factor, uematsu2022estimation, and pesaran2025identifying.
In comparison to this literature, our main methodological contribution is a unified procedure for constructing confidence intervals for signal strength across all four models and all signal ranges, as presented in Section (ref). This approach does not rely on standard Gaussian quantiles, but instead uses a novel random transition process $\mathcal{T}(\Theta)$, whose quantiles are tabulated in Table (ref). Figure (ref) reveals that our procedure performs well for a wide range of different signals. The quality of approximations is excellent for critical and weak signals, and does not deteriorate for larger signals, hence, covering also the case of strong signals. In contrast, as shown in Figure (ref), the Gaussian approximation performs poorly near the critical threshold, making $\mathcal{T}(\Theta)$ essential for accurate inference in that regime. Our approach is reminiscent of the construction of uniform confidence intervals for autoregressive models in stock1991confidence,mikusheva2007uniform, where the non-standard asymptotics near the unit root are smoothly connected to the standard normal behavior in the stationary region.
A surprising finding is that the same transition process $\mathcal{T}(\Theta)$ governs all four models. In fact, the proofs in Section (ref) follow different paths depending on the model, and only in the final step does a structural identity emerge, revealing that all four asymptotic distributions coincide. Random matrix theory has many universality theorems, and based on our results, we predict that the same transition process $\mathcal{T}(\Theta)$ governs a much wider class of signal plus noise models, beyond the ones analyzed here.
Beyond quantifying uncertainty in signal strength, our framework also enables signal detection and the assessment of factor informativeness. Specifically, one can check whether the uniform confidence intervals include zero and the identification threshold, respectively. paul2007asymptotics,onatski2012asymptotics,benaych2012singular,BG_CCA show that, for weak signals, estimates of the signal direction are inconsistent. Asymptotically, the estimated direction is inclined at an angle $\phi$ relative to the true direction. There are two complementary cases: if $\phi= \frac{\pi}{2}$, then the estimated direction contains no information about the truth and can be discarded; if $\phi<\frac{\pi}{2}$, then information is present and can potentially be extracted. It turns out that $\phi$ depends on the signal strength, decreasing as the strength increases, and that the transition between these two cases occurs precisely at the identification threshold discussed above. This reinforces our results on estimating signal strength in the critical regime as a tool for distinguishing between these two cases in direction estimation.
From a mathematical perspective, we develop a new approach to analyzing critical spikes, grounded in perturbation theory equations that relate the eigenvalues of spiked and unspiked random matrices. This contrasts with earlier treatments of critical spikes in real symmetric matrices, which relied on Pfaffian point processes (as in mo2012rank) or on tridiagonal matrix models (as in bloemendal2013limits, bloemendal2016limits, lamarre2019edge). Our central technical contribution is to show that these perturbation equations admit a well-defined edge-scaling limit, which captures the asymptotic behavior of the largest eigenvalues. While our approach is novel in all four settings, we particularly emphasize the fourth -- canonical correlations -- where no prior results on critical spikes were available.
In Section (ref) and (ref) we establish this edge limit result under two key assumptions on the unspiked model: (i) the asymptotics of the largest eigenvalues converge to the Airy$_1$ point process, and (ii) a form of the local law holds for the Stieltjes transform near the spectral edge. These assumptions are known to hold for a wide range of random matrix ensembles, including the four models considered in this paper. A notable strength of our approach is its minimal reliance on model-specific structure: we require only the two inputs above.
We build on some of the ideas in aizenman2015ubiquity. In contrast, however, we focus on the limit at the spectral edge—rather than in the bulk—which requires subtracting diverging counterterms. Moreover, we establish convergence in a stronger topology, which allows us to work directly on the real axis; see Appendix (ref) for further details.
Section (ref) introduces the four main signal plus noise models. Section (ref) presents a unified procedure for constructing confidence intervals for signal strengths. Section (ref) lays out the theoretical foundations underlying this procedure. Section (ref) offers three empirical illustrations. Extensions are discussed in Section (ref). Section (ref) concludes. All proofs are in Appendices (ref) and (ref).
In this section we present the four models, beginning with the simplest case -- the spiked Wigner model -- then proceeding to sample covariance and factor models based on PCA, and concluding with canonical correlation analysis (CCA). Although PCA-based models are the most widely used in practice, we adopt this order because the formulas are simpler in the Wigner case, making the key ideas more transparent.
Suppose we observe an $N\times N$ matrix $\mathbf A$ of the form
where $r$ is fixed (not growing with $N$) and $\theta_1>\dots>\theta_r>0\in\mathbb R$ are the strengths of $r$ signals, with corresponding directions $\mathbf u^*_1,\dots,\mathbf u^*_r$, which are assumed to be orthonormal $N$--dimensional vectors. The noise matrix $\mathcal E$ is a (Wigner) matrix sampled from the Gaussian Orthogonal Ensemble, meaning that $\mathcal E=\frac{1}{\sqrt{2 N}}(\mathcal Z+\mathcal Z^\mathsf T)$, where $\mathcal Z$ is an $N\times N$ matrix of i.i.d. $\mathcal N(0,\sigma^2)$ entries (see Section (ref) for non-Gaussian setting). We assume that $\theta_i$ and $\mathbf u^*_i$ are unknown deterministic parameters; one could alternatively allow $\mathbf u^*_i$ to be random, provided they are independent of $\mathcal E$. Our goal is to estimate the signal strengths $\theta_1,\dots,\theta_r$.
We first assume that the variance of the underlying noise $\mathcal Z$, $\sigma^2$, is known and set it to $1$ by rescaling the model.\footnote{The prefactor $\frac{1}{\sqrt{2N}}$ in the definition of $\mathcal E$ ensures that its eigenvalues fill the interval $[-2,2]$ as $N\to\infty$.} In Section (ref) we discuss adjustments for the case of unknown $\sigma^2$.
One common application of the spiked Wigner framework is modeling symmetric interaction networks, such as economic or social activity among $N$ agents. Each rank‑one component $\theta_i \mathbf u^*_i (\mathbf u^*_i)^\top$ captures a latent structure in agent attributes $\mathbf u^*_i$, while the observed interactions are contaminated by noise $\mathcal E$. Low-rank approximations of this form underpin seminal network models including the stochastic block model of Holland1983stochastic, where communities are inferred from block‑structured adjacency matrices, and latent space models.
The following result establishes the threshold for the estimation of $\theta_i$ via spectral methods.
Informally, the proposition says that “good” recovery of $\theta_i$ from the largest eigenvalues is possible if and only if $\theta_i$ is larger than the critical value $\theta^c=1$. In this case, to estimate $\theta_i$, one should take $\lambda_i$ and apply the inverse of the mapping $\theta\mapsto\lambda(\theta)$, which is $\lambda\mapsto \frac{1}{2}\left(\lambda+\sqrt{\lambda^2-4}\right)$.
We assess the quality of estimating $\theta_i$ by constructing a confidence interval for it. Specifically, for each fixed $i$ and significance level $\alpha$ we aim to find endpoints $\theta^-_i(\lambda_i,N,\alpha),\, \theta^+_i(\lambda_i,N,\alpha)$ such that
where $\approx$ denotes an $N\to\infty$ approximation, which should be uniform over the model parameters $\theta_1,\dots,\theta_r$ and $\mathbf u_1^*,\dots,\mathbf u_r^*$ in (ref).
In principle, since we deal with multiple $\theta_i$ simultaneously, one could consider joint multi-dimensional confidence sets. However, due to the asymptotic independence of $\lambda_i$ in (ref), it is sufficient to construct separate intervals for each $\theta_i$, which is the approach we take.\footnote{In contrast, if $\theta_i$ coincide, then the limits in (ref) are neither Gaussian nor independent, cf.\ capitaine2012central.}
The asymptotics (ref) provides a way to construct confidence intervals by approximating $\theta_i$ in the argument of $V(\theta_i)$ with $\theta(\lambda_i)= \frac{1}{2}\left(\lambda_i+\sqrt{\lambda^2_i-4}\right)$ and then using Gaussian quantiles. This leads to the following formula for the confidence interval:
where $z_{\alpha/2}$ denotes the $\alpha/2$ quantile of $\mathcal N(0,1)$. E.g., to obtain a $95\%$ confidence interval for a single fixed $i$, we set $z_{\alpha/2}=1.96$.
The formula (ref) reveals a problem as $\theta\to 1$ (i.e., $\lambda\to 2$): the confidence intervals diverge due to the $\sqrt{\lambda_i^2-4}$ singularity in the denominator. However, Monte Carlo simulations in Figure (ref) indicate that no such explosion actually occurs. This suggests that the approximation error in the confidence interval (ref) becomes non-negligible when $\lambda_i$ is close to 2, making the formula unreliable in this regime. In contrast, our novel procedure, introduced in Section (ref), closely matches the simulations across all values of $\lambda_i$.
Second, we consider a deterministic $N\times N$ matrix
where $r$ is fixed, $\theta_1>\dots>\theta_r>\sigma^2$ are the signal strengths, and $\mathbf u^*_1,\dots,\mathbf u^*_r$ are orthonormal $N$-dimensional vectors representing $r$ signal directions. The eigenvalues of $\Omega$ are $\theta_1,\theta_2,\dots,\theta_r$, and $\sigma^2$ with multiplicity $(N-r)$. As before, we assume that $\sigma^2$ is known and set it to $1$ without loss of generality; adjustments for unknown $\sigma^2$ are discussed in Section (ref).
We observe an $N\times S$ data matrix $X$, whose columns are i.i.d. $\mathcal{N}(0,\Omega)$, and aim to estimate $\theta_1,\dots,\theta_r$ from the sample covariance matrix $\frac{1}{S} X X^\mathsf T$. This model has been central in statistics and random matrix theory since johnstone2001distribution; see johnstone2018pca for a comprehensive overview, historical context, and many practical examples. Typically, the dimension $S$ reflects multiple independent observations: across individuals, measurement points, time periods, etc. An exact analogue of (ref) holds in this setting as well.
As in the previous section, we can use this Gaussian approximation to construct confidence intervals for each $\theta_i$, yielding a modification of (ref). However, this approach faces the same issue: the intervals become unreliable as $\lambda_i$ approaches $\lambda_+$ and must be corrected.
For the third setup we consider a random $N\times S$ matrix $X$ defined by
where $r$ is a fixed small number, $\theta_1,\dots,\theta_r>0$ are the signal strengths, $\mathbf u^*_1,\dots,\mathbf u^*_r$ are $N$-dimensional orthonormal vectors of signal directions, called “loadings”\footnote{Sometimes $\{\sqrt{\theta_i}\mathbf u^*_i\}$ rather than $\{\mathbf u^*_i\}$ are referred to as loadings.}, and $\mathbf v^*_1,\dots,\mathbf v^*_r$ are $S$-dimensional orthonormal vectors called “factors”. The noise matrix $\mathcal E$ has independent $\mathcal N(0,\sigma^2)$ entries. For now and until Section (ref) we assume $\sigma^2$ to be known and set it to $1$.
Our goal is to estimate $\theta_1,\dots,\theta_r$ from the eigenvalues $\lambda_1\ge \lambda_2\ge \dots\ge \lambda_N$ of the sample covariance matrix $\frac{1}{S} X X^\mathsf T$. While the factor model has similarities to the spiked covariance model of the previous section, they are not equivalent, because we treat $\sqrt{\theta_i} \cdot \mathbf u^*_i (\mathbf v^*_i)^\mathsf T$ in (ref) as deterministic parameters (the models would have been equivalent up to shift $\theta_i\to \theta_i+\sigma^2$, if each $\sqrt{S}\mathbf v^*_i$ were a mean $0$ Gaussian vector with i.i.d.\ components). This distinction allows the factor model to capture complex structures along the $S$-dimension, which is essential in applications across finance, macroeconomics, natural sciences, and other fields. Once again, an analogue of (ref) holds.
As in Section (ref), the Gaussian approximation of $\lambda_i$ leads to two methods for constructing confidence intervals. An analogue of (ref) is $$ \theta_i\in \left[\theta(\lambda_i) - \frac{\sigma(\lambda_i)}{\sqrt{N}} z_{\alpha/2}, \theta(\lambda_i) + \frac{\sigma(\lambda_i)}{\sqrt{N}} z_{\alpha/2}\right], \qquad \text{ where } $$ $$ \theta(\lambda)=\frac{\lambda-1-\gamma^2+ \sqrt{(1+\gamma^2-\lambda)^2-4\gamma^2}}{2}, \quad \sigma(\lambda)= \frac{\sqrt{2 \gamma^2 (2 \theta(\lambda)+1+\gamma^2)(\theta(\lambda)^2-\gamma^2)}}{\sqrt{(1+\gamma^2-\lambda)^2-4\gamma^2}}. $$ There is also a direct analogue of the second Gaussian method described in Remark (ref). Figure (ref) compares these two Gaussian-based intervals with our new approach, presented in Section (ref). The comparison reveals the same key features as in the spiked Wigner model.
For the final setup we fix a small integer $r$ and parameters $1\ge \theta_1,\dots,\theta_r\ge 0$. We consider a deterministic symmetric positive-definite $(N+M)\times (N+M)$ matrix $\Omega$ that satisfies
where $A$ and $B$ are $N\times N$ and $M\times M$ matrices, respectively, $I_N$ and $I_M$ are identity matrices of $N\times N$ and $M\times M$ dimensions, respectively, and $\mathrm{diag}(\sqrt{\theta_1},\dots,\sqrt{\theta_r})$ is a rectangular matrix with $\sqrt{\theta_1},\dots,\sqrt{\theta_r}$ on the first $r$ elements of the main diagonal and $0$ everywhere else.
Let $\mathbf x$ be an $(N+M)$--dimensional Gaussian mean $0$ random vector with covariance $\Omega$, and let $\mathbf u$ and $\mathbf v$ denote its first $N$ and last $M$ coordinates, respectively. The parameters $\theta_1,\dots,\theta_r$ are the squared canonical correlations between $\mathbf u$ and $\mathbf v$; see BG_review, as well as classical statistics references such as thompson1984canonical,gittins1985canonical,anderson1958introduction,muirhead2009aspects for detailed introductions to canonical correlation analysis (CCA). Algorithmically, $\theta_i$ are the largest eigenvalues of the matrix $(\mathbb E \mathbf u \mathbf u^\mathsf T)^{-1} \mathbb E \mathbf u \mathbf v^\mathsf T (\mathbb E \mathbf v\mathbf v^\mathsf T)^{-1} \mathbb E \mathbf v \mathbf u^{\mathsf T}$.
Given $S$ independent samples of $\mathbf x$, we construct two matrices: the $N\times S$ matrix $\mathbf U$ has $S$ samples of $\mathbf u$ as columns and the $M\times S$ matrix $\mathbf V$ has $S$ samples of $\mathbf v$ as columns. The sample squared canonical correlations $\lambda_1\ge \lambda_2\ge \dots$ are the eigenvalues of the $N\times N$ matrix $(\mathbf U \mathbf U^\mathsf T)^{-1} \mathbf U \mathbf V^\mathsf T (\mathbf V \mathbf V^\mathsf T)^{-1} \mathbf V \mathbf U^\mathsf T$. Our goal is to estimate $\theta_1,\dots,\theta_r$ from these eigenvalues.
In typical applications CCA is used to explore dependencies between two data sets, for example, two sets of individual characteristics, brain measurements versus behavioral scores, or two groups of stocks. The parameter $\theta_i$ quantify the strength of these dependencies. Once again, an analogue of (ref) holds.
The same conclusion applies here: using (ref) and Gaussian quantiles we can construct confidence intervals for $\theta_i$ that perform well when $\lambda_i$ is bounded away from $\lambda_+$, but become inaccurate as $\lambda_i$ approaches $\lambda_+$ and, therefore, require correction.
In this section we present our algorithm for constructing confidence intervals and explain how they can be interpreted and used to distinguish between noise, non-informative signals, and meaningful signals. We begin by introducing the transition process $\mathcal T(\Theta)$ and its properties, and then show how to use it to construct confidence intervals. We present two approaches: The first one is based on pre-tabulated quantiles of $\mathcal T(\Theta)$. The second one relies on bootstrap-type methodology. The underlying theorems will be presented in Section (ref).
As highlighted in Figure (ref), the Gaussian limits in (ref), (ref), (ref), and (ref) ought to be replaced by a different limiting object, which we call the transition process $\mathcal T(\Theta)$. This is a random function of $\Theta \in \mathbb{R}$. Its formal definition is provided in Section (ref), while for the purposes of constructing confidence intervals, the key quantities of interest are the quantiles of its distribution, which may be computed as follows:
The transition process $\mathcal T(\Theta)$, with appropriate centering and scaling, can be used to approximate the fluctuations of the largest eigenvalues, leading to the following algorithm.
Theorem (ref) and Corollary (ref) establish that the approximation (ref) is valid both when $\theta$ is bounded away from the critical value $\theta^c$ and when $\theta$ is close to $\theta^c$. One can also show that the approximation remains valid as $\theta\to\infty$ (corresponding to strong signals)\footnote{In the case of CCA, $\theta$ represents a correlation and is therefore bounded by $1$.}; we omit the proof since the present work focuses on the weak-signal regime. As is evident from the simulations in Figure (ref), the quality of the approximation improves as $\theta \to \infty$; we provide one formal result in this direction at the end of Section (ref).
Our results further show that, in many cases, the approximations for different $\lambda_i$ are asymptotically independent. Consequently, we can use (ref) as a foundation for constructing confidence intervals.
We begin with the case where the noise variance $\sigma^2$ is known, as specified in Section (ref). A simple rescaling allows us to assume $\sigma^2 = 1$ without loss of generality. The algorithm then proceeds as follows:
The first step is to draw a histogram of all eigenvalues $\lambda_1,\lambda_2,\dots$. In the settings of (ref), (ref), (ref), or (ref), the histogram should resemble a known limiting shape; namely, the semicircle law, Marchenko-Pastur law, or Wachter law, depending on the model, as detailed in Table (ref), with parameters specified in Table (ref), see Appendix (ref) for more details. If the histogram is reminiscent of one of these shapes, we regard the modelling assumptions as valid and apply Procedure (ref) to construct confidence intervals. Section (ref) discusses possible extensions when the empirical histogram deviates from the expected limit shape.
For the second step, we choose a significance level $\alpha$ (or confidence level $1-\alpha$) and, using Section (ref), construct two deterministic functions $t_{\alpha/2,+}(\Theta)$ and $t_{\alpha/2,-}(\Theta)$ such that
Following (ref) and using the parameter choices from Table (ref), we rescale the functions $t_{\alpha/2,\pm}(\Theta)$ to obtain $\widehat{t}_{\pm}(\theta)$, defined as
For the third step, we fix an index $i$ and consider the $i$th largest eigenvalue $\lambda_i$, such that $\lambda_i > \lambda_+$. We then determine two numbers $\theta_-<\theta_+$ such that
The procedure amounts to plotting the functions $\theta\mapsto \widehat{t}_{\pm}(\theta)$ and finding their intersection with the horizontal line $y=\lambda_i$. The resulting $[\theta_-,\theta_+]$ serves as the confidence interval for the $i$th signal strength $\theta_i$. We have $\mathrm{Prob}(\theta_i\in [\theta_-,\theta_+])\to 1-\alpha$ as $N\to\infty$ by Corollary (ref).
There are two special cases to consider at this step. First, it may happen that no value $\theta_-$ satisfies $\widehat{t}_{+}(\theta-) = \lambda_i$. This occurs when the shifted and rescaled $\lambda_i$ falls below the $(1-\alpha/2)$ quantile of the Tracy-Widom distribution $F_1$. In this case the confidence interval becomes one-sided, and one should set $\theta_-=-\infty$ or, equivalently, to the lower bound of admissible values of $\theta_i$, that is $\theta_i\ge \sigma^2$ for the spiked covariance and $\theta_i\ge 0$ for the others. Second, it may happen that $\theta_-$ exists, but lies below the lower bound for admissible values of $\theta_i$. In this case $\theta_-$ should again be replaced by the appropriate lower bound. In terms of statistical consequences the two cases are equivalent.
An important and unexpected feature of our asymptotic results (see Theorem (ref) for details) is that for a fixed index $i$ the asymptotic approximation of $\lambda_i$ depends only on $\theta_i$, but not on other parameters of the models (ref), (ref), (ref), (ref), such as $\{\theta_j\}_{j\ne i}$ or vectors $\mathbf u^*_i$. Additionally, the rank $r$ does not enter into the formulas. Hence, instead of relying on formulas (ref) for the functions $ \widehat{t}_{\pm}(\theta)$, we can obtain them through a bootstrap procedure using only the $r=1$ case.
{\bf Alternative formula-free second step.} Consider one of the models (ref), (ref), (ref), or (ref) with the desired matrix sizes $N$, $S$, and $M$, and instead of the true rank set $r=1$. Fix the direction of the unique signal arbitrarily; for example, set $\mathbf u^*_1$ (and $\mathbf v^*_1$ for the factor model, or the canonical variables for CCA) to the first coordinate vector. The model then depends on a single remaining parameter, $\theta_1$. Discretize $\theta_1$ on a grid, and for each value compute the largest eigenvalue $\lambda_1$ of the corresponding model matrix ($\mathbf A$ for the spiked Wigner matrix, $\frac{1}{S} X X^\mathsf T$ for the spiked covariance and factor models, and $(\mathbf U \mathbf U^\mathsf T)^{-1} \mathbf U \mathbf V^\mathsf T (\mathbf V \mathbf V^\mathsf T)^{-1} \mathbf V \mathbf U^\mathsf T$ for CCA). By repeating sufficiently many Monte-Carlo simulations, compute the $\alpha/2$ and $(1-\alpha/2)$ quantiles of $\lambda_1$ for each $\theta_1$, which yield the desired $\widehat{t}_{-}(\theta)$ and $\widehat{t}_{+}(\theta)$. This procedure produces the thick gray curves in Figure (ref). Then proceed to the third step as in Section (ref).
An advantage of the algorithm in this section is that (ref) and Tables (ref) and (ref) are not required. However, this comes at the cost of running many Monte-Carlo simulations, making the implementation slower than that of the algorithm in Section (ref).
For the CCA setting in Section (ref) the asymptotics in (ref) does not depend on the noise covariance, i.e., the matrices $A$ and $B$ in (ref). In contrast, for the other three settings, Sections (ref), (ref), and (ref), the scaling depends on the noise variance, denoted by $\sigma^2$. The same holds for Theorem (ref), which underlies the algorithms for constructing confidence intervals in Sections (ref) and (ref). In particular, Procedure (ref) assumes $\sigma^2 = 1$. If $\sigma^2 \neq 1$ but is known, then the entries of the data matrix $\mathbf A$ or $X$ should be divided by $\sigma$ to reduce to the baseline case $\sigma^2 = 1$. If $\sigma^2$ is unknown, it must first be estimated.
We propose estimating the variance by discarding $25\%$ of the eigenvalues at both ends and matching sample moments to their theoretical values to solve for $\sigma^2$.
For the spiked Wigner model, let $\ell\approx 0.81$ denote the positive number such that
and set
Given eigenvalues $\lambda_1\ge \dots\ge \lambda_N$ of $\mathbf A$, we can form an estimate
The Wigner semicircle law for the GOE with explicit estimates for the remainders (see e.g., o2010gaussian), combined with the interlacing inequalities between the eigenvalues of $\mathbf A$ and $\mathbf B$ in (ref), as in Corollary (ref), can be used to show that
The scale of the random component in Theorem (ref) is much larger than the error term in (ref). Henc3, our confidence intervals are much wider than this error term and normalizing the data by $\widehat \sigma$ does not change the validity of the confidence intervals of Sections (ref), (ref).
For the spiked covariance and factor models, the procedure is analogous, but relies on the Marchenko-Pastur law (see Table (ref)) rather than the semicircle law. Fixing the parameter $\gamma^2=\frac{N}{S}\in [0,1)$, we define $\ell_-$ and $\ell_+$ as two positive numbers such that
and set
Given eigenvalues $\lambda_1\ge \dots\ge \lambda_N$ of $\frac{1}{S} X X^\mathsf T$, we can form an estimate
The Marchenko--Pastur law with explicit estimates for the reminders (see e.g., bourgade2022optimal), combined with the interlacing inequalities between the eigenvalues of spiked and unspiked models, as in Corollaries (ref) and (ref), can be used to show that
Once again, normalizing the data by dividing by $\widehat \sigma$ does not affect the validity of the confidence intervals constructed in the previous section.
For a discussion of alternative procedures for estimating $\sigma^2$ see, for example, kritchman2009non, shabalin2013reconstruction, gavish2014optimal, or ke2023estimation.
Confidence intervals play two key roles in the analyzing signal strength. First, they measure uncertainty: the narrower the interval, the more precisely the signal is estimated.
Second, they assess the informativeness of estimated signals. If the lower bound starts at $-\infty$ or at the minimal admissible value of $\theta$, (which is $\sigma^2$ for the spiked covariance and zero for other models), the signal may be spurious and reflect pure noise. Alternatively, if the interval is bounded away from the minimal admissible value, but contains the identification threshold $\theta^c$, then a signal exists, but we cannot reject that its strength falls below the cutoff. When $\theta \leq \theta^c$ the sample estimates of the directions $\mathbf u$ and $\mathbf v$ are asymptotically orthogonal to their true counterparts (see, e.g., paul2007asymptotics,onatski2012asymptotics,benaych2012singular,johnstone2018pca,BG_CCA), rendering the signal effectively non-informative.
The new asymptotics, which improves upon the Gaussian approximations (ref), (ref), (ref), and (ref), is based on a novel stochastic object we call the Airy--Green function. Its definition, along with the transition process $\mathcal T(\Theta)$ constructed from it, is presented in Section (ref). A discussion of its nature is provided in Section (ref). Theorem on the convergence of the eigenvalue distributions in four models towards this object is stated in Section (ref).
We begin by recalling the Airy$_1$ point process, a random sequence of points $\mathfrak a_1 \ge \mathfrak a_2 \ge \mathfrak a_3 \ge \dots$, which can be defined as the scaling limit of the largest eigenvalues of Wigner matrices.
Similar asymptotic results hold for the other models we consider. All existing formulae for the finite-dimensional distributions of $\{\mathfrak a_i\}_{i=1}^{\infty}$ are quite complicated and do not provide explicit distribution function, see, e.g., forrest. Nevertheless, the distribution can be sampled and tabulated, see bornemann2009numerical,vignette_largevars. In particular, Table (ref) lists quantiles of the Tracy–Widom distribution, which describes the law of $\mathfrak a_1$.
The following theorem defines the Airy--Green function $\mathcal G(w)$; see Section (ref) for the proof.
Note that any fixed $w \in \mathbb C$ is almost surely not in $\{\mathfrak a_j\}$, hence, for such $w$ the convergence holds almost surely. Eq.\ (ref) and Proposition (ref) imply that $\mathcal G(w)$ changes monotonically from $+\infty$ to $-\infty$ over the interval $[\mathfrak a_1,+\infty)$, allowing us to state the following key definition.
Figure (ref) and Table (ref) show the simulated quantiles for the random variables $\mathcal T(\Theta)$ as functions of $\Theta$, or equivalently, confidence intervals for $\Theta$ as a function of $\mathcal T$. These results are based on $MC=10^6$ Monte Carlo simulations of the $\sqrt{N}\times\sqrt{N}$ top-left corners of $N\times N$ tridiagonal matrices of dumitriu_edelman, with a perturbed $(1,1)$ matrix element and $N=10^8$; see edelman2005numerical and johnstone2021spin for justifications of this approach. The figure shows that the Gaussian approximation from Proposition (ref) performs well for large $\Theta$, but deteriorates near $\Theta = 0$.
The term “Airy” in the name $\mathcal G(w)$ refers to the Airy point process, whose points $\mathfrak a_i$ appear in its definition. The term “Green” stands from the tradition in random matrix theory to refer to matrix elements of the resolvent $(z I-D)^{-1}$ of a symmetric matrix $D$ as the Green's function. Via eigenvalue decomposition, the $(1,1)$ matrix element of $(z I -D)^{-1}$ is $\sum_i \frac{u_{1i}^2}{z-d_i}$, where $u_{1i}$ is the first coordinate of the $i$th normalized eigenvector of $D$ corresponding to the eigenvalue $d_i$, making it reminiscent of the sum in (ref).
The term “transition” in $\mathcal T(\Theta)$ refers to its role in capturing the transition between subcritical $\theta<\theta^c$ and supercritical $\theta>\theta^c$ behavior in (ref), (ref), (ref), (ref). This phenomenon is commonly known as the BBP phase transition, following baik2005phase.
There are two other approaches to the transition process $\mathcal T(\Theta)$ in the literature. One is based on the limit of tridiagonal matrix model: bloemendal2013limits,lamarre2019edge construct $\mathcal T(\Theta)$ as the largest eigenvalue of the Stochastic Airy Operator with $\Theta$--dependent boundary condition. Another approach, developed by mo2012rank using the framework of Pfaffian point processes, provides an integral representation for the one-dimensional marginal distribution of $\mathcal T(\Theta)$.
An advantage of our definition via the Airy--Green function is its robustness. Proving convergence to either of the two alternative definitions, requires finding delicate algebraic structures (tridiagonalization or Pfaffians) in the prelimit objects, which are not known in some cases (e.g., CCA). In contrast, our approach relies only on identifying the eigenvalues of a spiked model as solutions to an equation, that can be obtained in all spiked models via finite-rank perturbation theory.
The next theorem presents the asymptotics of the largest eigenvalues for all signal plus noise models of Section (ref).
Note that the two limit regimes (ref) and (ref) heuristically agree with each other: if one sets $\tilde \theta = \varepsilon N^{1/3}$ with a small $\varepsilon>0$, then using (ref) and Proposition (ref), we expect $$ \lambda_q \approx \lambda_++ N^{-2/3} \kappa_1 \mathcal T(\kappa_2 \tilde \theta)\approx \lambda_+ + \kappa_1 \kappa_2^2 \varepsilon^2 + 2N^{-1/2} \kappa_1 \sqrt{\varepsilon \kappa_2} \mathcal N(0,1). $$ If one sets $\theta_i=\theta^c+\varepsilon$, then Taylor expanding (ref) (noting $V(\theta^c)=\lambda'(\theta^c)=0$), we expect $$ \lambda_i\approx \lambda_+ + \frac{\varepsilon^2}{2} \lambda''(\theta^c) + N^{-1/2} \sqrt{\varepsilon V'(\theta^c)} \mathcal N(0,1). $$ Using (ref), we see that the last two asymptotic expansions are the same. In parallel, using Proposition (ref) we can combine two asymptotic regimes of Theorem (ref) into one (among several asymptotically equivalent formulas, we chose the one with the best finite sample performance):
The formula (ref) is a direct corollary of (ref), while (ref) combines (ref) and (ref) together. Indeed, when $\theta_i$ is bounded away from $\theta^c$, (ref) converts (ref) into $$ \lambda(\theta_i)-\frac{ \kappa_2^{3/2}}{2} \sqrt{V(\theta_i) (\theta_i-\theta^c)^3}+ \frac{\kappa_2^{-1/2}}{2} \sqrt{\frac{V(\theta_i)}{\theta_i-\theta^c}}\kappa_2^2 (\theta_i-\theta^c)^2 + \frac{\kappa_2^{-1/2}}{N^{1/2}} \sqrt{\frac{V(\theta_i)}{\theta_i-\theta^c}} \sqrt{\kappa_2(\theta_i-\theta^c)} \mathcal{N}(0,1), $$ which is readily seen to be equivalent to (ref). When $\theta_i$ is close to $\theta^c$, $\theta_i=\theta^c + N^{-1/3} \tilde \theta$, Taylor expanding $\lambda(\cdot)$ and $V(\cdot)$ near $\theta^c$, (ref) turns into an equivalent form of (ref): $$ \lambda_++\frac{\lambda''(\theta^c)}{2}(\theta_i-\theta^c)^2-\sqrt{V'(\theta_c)}\frac{ \kappa_2^{3/2}}{2} (\theta_i-\theta^c)^2+ \frac{\kappa_2^{-1/2}}{2 N^{2/3}} \sqrt{V'(\theta^c)} \mathcal T \Bigl(\kappa_2 N^{1/3} (\theta_i-\theta^c) \Bigr), $$
Since $\kappa_3/N=o(N^{-2/3})$, and therefore the choice of $\kappa_3$ does not affect the validity of the asymptotic formulas (ref) and (ref). These terms are introduced to improve the performance of the formulas for intermediate values of $N$, cf.\ Johnstone_Jacobi,ma2012accuracy,johnstone2012fast, which emphasize the importance of $1/N$ corrections for the practical applicability. The reasoning behind our choice of $\kappa_3$ is as follows. First, we require continuity at $\theta^c$; hence, (ref) and (ref) use the same $\kappa_3/N$. Second, we leverage additional information available at $q = r = 1$ and $\tilde \theta = -N^{1/3}\theta^c$ for the spiked Wigner, factor, and CCA models. (For the spiked covariance model, one instead takes $\tilde \theta = -N^{1/3}\gamma$ and adjusts the formula accordingly.) On one hand, combining (ref) with the asymptotic approximation (ref), we obtain
On the other hand, $\tilde \theta= - N^{1/3}\theta^c$ corresponds to $\theta=0$ (or $\theta = 1$ for the spiked covariance model with $\tilde \theta = -N^{1/3}\gamma$), meaning that in all four models, we are in the unspiked regime without a signal component. In this setting, the convergence of $\lambda_1$ to the Tracy–Widom distribution $\mathfrak a_1$ is well established, and $1/N$-order asymptotic corrections have been studied in Johnstone_Jacobi,ma2012accuracy,johnstone2012fast. From these works, one can extract
Equating (ref) with (ref) yields the formula for $\kappa_3$, as recorded in Table (ref). An interesting observation is that in each case the term $\frac{\kappa_1}{\kappa_2 \theta^c}$ in (ref) is twice the $\frac{1}{N}$ correction term in (ref), which leads to the $\frac{3}{2}$ coefficient appearing in $\kappa_3$ across all four models.
The proof of Theorem (ref) in Section (ref) begins with a rank-one perturbation equation, which expresses the eigenvalues in a signal plus noise model with $r$ spikes (the “target model”) as solutions to an algebraic equation involving a simpler model with $r - 1$ spikes (the “base model”). Analyzing the asymptotic behavior of this equation leads to the following conclusion, stated informally below:
On the technical level, the key novelty is in our ability to handle the most delicate case, when the spike is critical. If all spikes are subcritical or supercritical, the arguments are much simpler and follow ideas similar to those found in the references cited in Section (ref). Some special cases of Theorem (ref) can be handled by other methods, for example, the $r=1$ case for the spiked Wigner and spiked covariance models is addressed in mo2012rank,bloemendal2013limits; see also bloemendal2016limits,lamarre2019edge. However, we believe that the level of generality achieved here -- particularly our treatment of the factor model and CCA -- was not previously available in the literature and is beyond the reach of those alternative methods.
We present three examples that illustrate the application of the procedure described in Section (ref) to empirical data sets. In each case, the first step reveals a strong agreement between the histogram of eigenvalues and the corresponding theoretical curve—specifically, the Marchenko-Pastur law in the first two examples (factor models) and the Wachter law in the third (CCA). This suggests that the data aligns well with our modeling assumptions.
Industrial production (IP) accounts for more than $10\%$ of the United States' Gross Domestic Product (GDP), making it a significant component of total U.S. output. In this subsection we use data from andreou2019inference, which investigates whether IP constitutes a dominant factor in U.S. economic activity. The data set contains quarterly IP growth rates across 117 sectors, spanning the period from $1977:Q1$ to $2011:Q4$. We de-mean the data and standardize each sector to have unit sample variance, thereby working with the sample correlation matrix of IP.
Figure (ref) shows all eigenvalues of the standardized IP and highlights four of them ($3.75,\, 4.26,\, 6.12,\, 30.05$) that lie to the right of the theoretical Marchenko–Pastur upper edge, $\lambda_+ = 3.68$. The largest eigenvalue, $30.05$, stands out markedly and represents a strong “market” factor. The two smallest among the four, being close to the cutoff, may reflect spurious signals arising from noise. To assess their significance, we construct $95\%$ confidence intervals. These intervals are represented by vertical segments in Figure (ref) and are also summarized in Table (ref).
Notably, the interval for the fourth largest eigenvalue at $3.86$ differs substantially from what one would obtain using a Gaussian approximation (discussed in Section (ref) and at the end of Section (ref)), highlighting the importance of our new procedure. This interval intersects the critical identification threshold $\theta^c = 0.92$, indicating that we cannot reject the null hypothesis that it represents noise (or a non-informative signal). For the remaining eigenvalues, the null is rejected. The confidence interval for the largest eigenvalue is nearly identical under our method and the Gaussian approximation, whereas the differences grow as the eigenvalues and corresponding signal strengths decrease.
The discussion in the previous paragraph leads to the conclusion that the IP growth rate is driven by three factors: one strong factor and two weak factors. It is instructive to compare this result with the classical information criterion of bai2002determining, one of the most commonly used methods in applied research and one that is based on strong-signal asymptotics. Using $IC_{p_2}$ from bai2002determining, we identify only the largest market factor. If we instead use $IC_{p_1}$ from the same article, which imposes a smaller penalty on the number of factors, we select the two largest factors. As a result, informative weak factors are missed by these procedures, underscoring the importance of accounting for weak signals and motivating our proposed methodology.
Analyzing stock returns is essential for understanding market dynamics, evaluating investment performance, and guiding both individual and institutional investment strategies. A key statistical object in this context is the covariance matrix of stock returns, which plays a central role in portfolio optimization, such as in the Markowitz mean–variance framework. The vast “factor zoo” -- the large number of potential variables proposed to explain stock returns -- highlights the practical challenge of distinguishing meaningful factors from noise in high-dimensional settings (see cochrane2011presidential for an influential discussion). Here we demonstrate how our methodology can be applied to the sample covariance matrix of weekly S$\&$P$100$ stock returns. We use data from BG1, which covers 92 stocks over the period from January 1, 2010, to January 1, 2020.
Before turning to the empirical spectrum, it is important to address a potential concern regarding the applicability of PCA-based methods in financial data. Financial returns may exhibit heavy tails, complicating factor identification via PCA. For pure-noise covariance matrices, a sharp transition occurs at tail index $\alpha = 4$ yin1988limit,bai1988note: lighter tails keep the largest eigenvalue at the Marchenko–Pastur edge, while heavier tails can produce spurious spikes. Although short-horizon returns often have $\alpha \approx 3$, weekly and longer-horizon returns typically have $\alpha\geq 5$ gopikrishnan1999scaling,gabaix2009power,fan2017elements, as aggregation and the central limit theorem dampen extremes. This places weekly returns in a regime where spurious tail-driven spikes are less likely to arise.
Figure (ref) presents the full spectrum of the sample covariance matrix and identifies ten eigenvalues, $10^{-4}(365,\, 52.4,\, 41.7,\, 30.7,\, 26.1,\, 20.7,\, 17.6,\, 15.1,\, 13.9,\, 13.6)$, that lie to the right of the theoretical Marchenko–Pastur upper edge, $\lambda_+ = 13\times10^{-4}$. To better fit the empirical data, we adopt an effective parameter value $\gamma^2 = 0.4$, in contrast to the true value $N/S = 0.18$. This adjustment may reflect temporal dependence in the data, which effectively reduces the sample size and increases $\gamma^2 = 0.4$. We also set $\sigma = 0.02$ to align the overall variance of the eigenvalues in the data.
Figure (ref) and Table (ref) report the 95$\%$ confidence intervals for ten candidate signals. As in the previous example, the largest eigenvalue is much larger than the others and corresponds to the “market” factor. The two smallest eigenvalues among them yield intervals that intersect the identification threshold $\theta^c = 3.11\times 10^{-4}$, indicating that they cannot be statistically distinguished from being non-informative. This is also the region in which the two Gaussian approximations produce markedly different intervals, reflecting the fact that the variance term $V(\theta_i)$ is very close to zero. We therefore conclude that only eight of the ten observed spikes represent informative signals. In contrast, applying the information criterion $IC_{p_2}$ of bai2002determining selects only four of these eight factors, while $IC_{p_1}$ identifies five; both criteria miss many of the weaker factors.
Financial stocks are typically classified into cyclical and non-cyclical (defensive) categories, depending on whether their performance tracks economic business cycles. These groups are generally assumed to be uncorrelated, aside from exposure to a common “market” factor. BG_CCA identify three non-zero canonical correlations between these two groups, suggesting the presence of three common factors. Here we revisit their analysis using the same data set to assess whether these observed correlations reflect genuine signals or could instead be attributed to noise.
The data set comprises weekly returns for $80$ cyclical and $80$ defensive stocks, spanning the period from January 1, 2010, to January 1, 2020. Figure (ref) reproduces the canonical correlations reported by BG_CCA and augments it with $95\%$ confidence intervals. BG_CCA used the signal strengths to estimate the angles between true and estimated canonical variables. By incorporating our results, one can generate confidence intervals for these angles.
As shown in Figure (ref), the $95\%$ confidence intervals for three largest canonical correlations, $0.58$, $0.62$, and $0.89$, lie above the cutoff $\theta^c=0.18$, confirming them as true signals. In contrast, the confidence interval for the fourth largest value intersects the cutoff, indicating that it cannot be reliably distinguished from noise and is therefore classified as a non-informative component.
A comparison of Figures (ref), (ref), and (ref) reveals an interesting pattern: in the factor models of the first two figures, the confidence intervals widen as the signal strengths increase, whereas in the CCA setting, the intervals become narrower as the signals approach 1. Theoretically, this behavior in CCA can be attributed to the factor $(1-\theta)^2$ in $V(\theta)$, as shown in Table (ref).
In this section we discuss possible extensions of our results, focusing on non-Gaussian data and broader classes of models than those considered in Section (ref).
The four models in Section (ref) are based on Gaussian noise matrices. A natural generalization is to replace the Gaussian vectors with more general random vectors having the same mean and covariance. It is well known (see lee2014necessary,ding2018necessary,FanYang and the broader reviews deift2009random,tao2012random,erdos2017dynamical) that for pure noise models without signals, the distribution of the largest eigenvalues remains unchanged in many non-Gaussian settings. This raises the question of whether the robustness extends to the confidence intervals of Section (ref).
Figure (ref) shows the results of Monte Carlo simulations of confidence intervals for Wigner matrices with a single spike and non-Gaussian noise. These can be compared to the Gaussian case in Figure (ref). The thick gray line depicts sample confidence intervals based on $10^5$ simulations of $100 \times 100$ matrices. The noise matrix $\mathcal E$ in (ref) is still defined as $\mathcal E = \frac{1}{\sqrt{2N}}(\mathcal Z + \mathcal Z^\mathsf T)$, with $\mathcal Z$ having i.i.d. entries. We vary the distribution of these entries, each rescaled to have mean $0$ and variance $1$, to be: (i) uniform on $[0,1]$, (ii) Bernoulli with success probability $p=1/2$, or (iii) Binomial with parameters $n=2$, $p=0.3$. We also consider two signal vectors: a localized signal $\mathbf u^* = (1,0,\dots,0)^\mathsf T$ and a delocalized signal $\mathbf u^* = \tfrac{1}{\sqrt{N}}(1,1,\dots,1)^\mathsf T$. While in the Gaussian setting the choice of $\mathbf u^*$ is irrelevant due to rotational invariance, this is no longer true for general noise distributions.
Figure (ref) shows that, while the sample confidence intervals remain reasonably close to those constructed via the procedure in Section (ref), the agreement is notably better for delocalized signals. For the localized signal significant deviations appear at larger values of the observed eigenvalue $\lambda_1$. A heuristic explanation follows directly from the model $\mathbf A=\theta \cdot \mathbf u^* (\mathbf u^*)^\mathsf T + \mathcal E$ in (ref). When $\mathbf u^* = (1,0,\dots,0)^\mathsf T$, the parameter $\theta$ enters $\mathbf A$ only through its sum with the $(1,1)$ entry of $\mathcal E$, so the distribution of that single element directly affects any estimate of $\theta$. By contrast, when $\mathbf u^* = \tfrac{1}{\sqrt{N}}(1,1,\dots,1)^\mathsf T$, the projection of the noise onto the signal direction aggregates many independent entries. Thus, by the central limit theorem, the influence of the individual noise distribution diminishes as $N$ grows, leading to behavior indistinguishable from the Gaussian case. Similar robustness is expected for other delocalized signals.
For spiked Wigner matrices with supercritical signals ($\theta > \theta^c$), non-Gaussian cases have been rigorously analyzed in several papers. For localized signals, non-trivial dependence of the asymptotics of $\lambda_1$ on the distribution of the noise has been established in capitaine2009largest,capitaine2012central,pizzo2013finite,knowles2013isotropic,knowles2014. In contrast, universality of the limit for delocalized signals has been shown under various conditions (including different definitions of “delocalized”) in feral2007largest,benaych2011fluctuations,capitaine2012central,pizzo2013finite,renfrew2013finite,knowles2013isotropic,knowles2014. Also knowles2013isotropic explain that even for localized signals the dependence of the limit on the noise distribution is washed out as $\theta$ approaches $\theta^c$. Similar phenomenology is expected (and in some cases proven) for other signal plus noise models.
This discussion, together with the simulation results, leads us to conjecture that for all signal plus noise models of interest, the confidence intervals constructed via the procedure in Section (ref) remain good approximations even under non-Gaussian noise, particularly when the signal is delocalized. From a practical standpoint, this reinforces the robustness of our method, especially in light of giannone2021economic, who argue that many economic data sets are non-sparse and should therefore be modeled using delocalized signals.
In all four models from Section (ref) the empirical eigenvalue distributions take specific parametric forms, as summarized in Table (ref). Although our estimation procedure relies only on the largest eigenvalues, the formulas in Table (ref) were derived using empirical distributions in Table (ref).
In other signal plus noise models and some empirical data sets the eigenvalue histograms may differ substantially from the four cases we consider. For $\theta > \theta^c$, the fluctuations of the largest eigenvalues have been analyzed in alternative settings; see, e.g., benaych2011fluctuations,benaych2012singular,onatski2012asymptotics. These works develop analogues of (ref), (ref), and (ref), but the formulas for $\theta^c$, $\lambda_+$, $\lambda(\theta)$, and $V(\theta)$ become more complex. We conjecture that the confidence interval procedure from Section (ref) remains valid in such broader contexts, with updated model-dependent parameters from Table (ref).
One applied setting of particular interest is the approximate factor model, which resembles the setup in Section (ref) but allows the noise matrix $\mathcal E$ to have a more complex correlation structure rather than being i.i.d. To apply our confidence interval procedure here, proceed as follows: first, select a value for $\lambda_+$; then use all eigenvalues below $\lambda_+$ to estimate the empirical distribution of noise eigenvalues (replacing the parametric forms in Table (ref)); next, substitute this estimate into the formulas from onatski2012asymptotics for $\lambda(\theta)$ and $V(\theta)$; and finally, apply our method to construct confidence intervals for the eigenvalues exceeding $\lambda_+$. Choosing $\lambda_+$ optimally is delicate; one approach, suggested in onatski2010determining, is to select it based on the characteristic $\sqrt{x}$ behavior of the eigenvalue density near the edge -- a feature clearly visible in the four models of Table (ref) and present in many other cases.
The paper presents a unified framework for conducting inference on signal strength in high-dimensional signal plus noise models, with a particular focus on the critical regime where standard Gaussian approximations fail. We demonstrate that the limiting distribution of top eigenvalues is governed by a universal stochastic process, the transition process $\mathcal T(\Theta)$, whose quantiles can be tabulated and used to construct valid confidence intervals. This approach applies uniformly across four canonical models: spiked Wigner matrices, spiked sample covariance matrices, factor models, and canonical correlation analysis.
Our procedure is robust to both weak and critical signals, enabling practitioners to distinguish between informative and non-informative components without imposing assumptions on signal strength. Our methodology reveals a surprising universality: despite differences in the statistical structure of the models, the same transition process governs the fluctuations of their top eigenvalues. This suggests deeper underlying principles in high-dimensional inference and opens an avenue for future research in more general signal plus noise settings.