Covariance test and universal bootstrap by operator norm

Building upon prior work, various approaches have been developed for analyzing $T$ and its variants in specific scenarios using random matrix theory. When $\bm{\Sigma}_{0}=\bm{I}_{p}$, the distribution of $T$ can be derived by examining the limiting behavior of the extreme eigenvalues of the sample covariance matrix $\bm{\hat{\Sigma}}$. The study of extreme spectral properties of high-dimensional covariance matrices has attracted substantial interest, resulting in a number of influential contributions, including Johnstone (2001); Bianchi et al. (2011); Onatski, Moreira and Hallin (2013); Johnstone and Paul (2018). Building on the progress of these studies, El Karoui (2007); Bao, Pan and Zhou (2015); Lee and Schnelli (2016); Knowles and Yin (2017) established that under conditions where the dimension $p$ and sample size $n$ grow proportionally and other certain regularity conditions, the limiting distribution of the largest eigenvalue of $\bm{\hat{\Sigma}}$ follows the Tracy–Widom law. In cases where $\bm{\Sigma}_{0}$ is invertible, Bao, Pan and Zhou (2015) proposed transforming the data $\bm{X}_{1},\cdots,\bm{X}_{n}$ into $\bm{\Sigma}_{0}^{-1/2}\bm{X}_{1},\cdots,\bm{\Sigma}_{0}^{-1/2}\bm{X}_{n}$, allowing for testing whether the covariance of $\bm{\Sigma}_{0}^{-1/2}\bm{X}_{i}$ is the identity matrix $\bm{I}_{p}$. This transformation results in the statistic $T$ taking the form of Roy’s largest root statistic

and the Tracy–Widom distribution results hold for $T^{\text{Roy}}$ when $p$ and $n$ grow proportionally, i.e., $\alpha=1$. However, similar distributional results for $T$ have not yet been obtained for general $\bm{\Sigma}_{0}$, even when $p/n$ is bounded. To address this gap, we leverage random matrix theory to characterize the extreme singular values for the matrix of the form $\hat{\bm{\Sigma}}+\bm{R}$ with a general matrix $\bm{R}$. When taking $\bm{R}=-\bm{\Sigma}_{0}$, this gives the limiting properties of $T$. Additionally, we introduce a novel technique, the universal bootstrap, designed to directly approximate the distribution of $T$ in the context of general $\bm{\Sigma}_{0}$.

The bootstrap method, a generic resampling technique to approximate the distribution of a statistic, was first introduced by Efron (1979). Originally designed for fixed-dimensional problems, the bootstrap method has recently been adapted for high-dimensional settings through substantial work. Chernozhukov, Chetverikov and Kato (2013, 2017); Lopes, Lin and Müller (2020); Lopes (2022b); Chernozhukov, Chetverikov and Koike (2023) developed Gaussian approximation methods and established the approximation rates for empirical bootstrap and multiplier bootstrap for the supremum norm of the sum of high-dimensional vectors. Moreover, Han, Xu and Zhou (2018); Lopes, Blandino and Aue (2019); Yao and Lopes (2021); Lopes (2022a); Lopes, Erichson and Mahoney (2023) explored bootstrap methods for spectral statistics of covariance matrices, while these works often assume a low intrinsic dimension or effective rank relative to the sample size $n$. Despite such advances in high-dimensional bootstrap methods, El Karoui and Purdom (2019) and Yu, Zhao and Zhou (2024) demonstrated that empirical and multiplier bootstraps can be inconsistent for $T$, even when $\bm{\Sigma}_{0}=\bm{I}_{p}$ if $p$ and $n$ grow proportionally and the eigenvalues of $\bm{\Sigma}$ and $\bm{\Sigma}_{0}$ do not decay.

To address the inconsistency issues associated with high-dimensional bootstraps based on Gaussian approximation, motivated by the concept of universality, a topic of recent interest in random matrix theory (Hu and Lu, 2022; Montanari and Saeed, 2022), we develop a novel resampling method informed by universality principles. Specifically, traditional bootstrap methods (Chernozhukov, Chetverikov and Kato, 2013, 2017; Chernozhukov, Chetverikov and Koike, 2023), depending on the high-dimensional central limit theorem, treat the large sample covariance matrix $\bm{\hat{\Sigma}}=\sum_{i=1}^{n}\bm{X}_{i}\bm{X}_{i}^{T}/n$ as a sum of random matrices and approximate its distribution with a Gaussian matrix sharing the same covariance as $\bm{\hat{\Sigma}}$. Nonetheless, as the dimension $p$ increases, the accuracy of this approximation diminishes, eventually leading to inconsistency, as shown by El Karoui and Purdom (2019). This inconsistency indicates that in high-dimensional contexts, the large matrix structures dominate the central limit theorem’s applicability, causing the sample covariance matrix $\bm{\hat{\Sigma}}$ to diverge from Gaussian behavior. To circumvent this limitation, the universal property leverages high-dimensional structures. Broadly, the universality suggests that the asymptotic distribution of $T$ only depends on the distribution of $\bm{X}_{1},\cdots,\bm{X}_{n}$ through their first two moments. This allows us to construct the universal bootstrap statistic by substituting independent Gaussian samples $\bm{Y}_{1},\dots,\bm{Y}_{n}\sim\mathcal{N}(\bm{0},\bm{\Sigma})$ in place of $\bm{X}_{1},\dots,\bm{X}_{n}$ in the definition of (1),

where $\bm{\hat{\Sigma}}^{\text{ub}}=\sum_{i=1}^{n}\bm{Y}_{i}\bm{Y}_{i}^{T}/n$. The key insight here is that while the universal bootstrap matrix $\bm{\hat{\Sigma}}^{\text{ub}}$ need not approximate a Gaussian matrix, it effectively replicates the structure of $\bm{\hat{\Sigma}}$, which is disrupted by the empirical and multiplier bootstraps. Although this method relies on sophisticated random matrix theory, the implementation remains straightforward.

We summarize our contribution as three-fold. First, we establish the anisotropic local law for matrices of the form $\hat{\bm{\Sigma}}+\bm{R}$ within both high- and ultra-high-dimensional settings (2), where $\bm{R}$ may have positive, negative, or zero eigenvalues. Erdős, Yau and Yin (2012); Knowles and Yin (2017) demonstrated that the anisotropic local law is essential for proving eigenvalue rigidity and universality. The anisotropic local law for $\hat{\bm{\Sigma}}$ has been previously established by Knowles and Yin (2017) for $\alpha=1$ and by Ding and Wang (2023); Ding, Hu and Wang (2024) for $\alpha>1$. However, the existence of $\bm{R}$ alters the structure of the sample covariance matrix $\hat{\bm{\Sigma}}$, rendering existing techniques invalid. Additionally, due to the potential ultra-high dimension, each block of the Green function matrix (see (21) for definition) may converge at different rates, making the bound derived in Ding and Wang (2023) suboptimal. To address these challenges, we introduce a new ”double Schur inversion” technique to represent the Green function matrix using local blocks. We also define an auxiliary parameter-dependent covariance matrix and establish the corresponding parameter-dependent Marčenko-Pastur law. This parameter-dependent structure enables us to prove the entry-wise local law, which represents a weaker version of the anisotropic local law. Furthermore, we present an improved representation of the anisotropic local law, enhancing the results of Ding and Wang (2023) and Ding, Hu and Wang (2024). This unified structure provides a valuable tool for establishing the anisotropic local law.

Second, we establish the universality results for the statistic $T$, based on which the universal bootstrap procedure is introduced. Leveraging this universality, we demonstrate that the empirical distribution generated by the universal bootstrap in (4) effectively approximates the distribution of $T$. Additionally, the asymptotic distribution of $T$ is shown to be independent of the third and fourth moments of $\bm{X}_{1},\cdots,\bm{X}_{n}$, allowing us to focus primarily on the covariance. In contrast, we find that two commonly used norms, the Frobenius norm and supremum norm, lack this universality, as their asymptotic distributions explicitly depend on all fourth moments of the data. Therefore, standardized statistics based on these norms requires fourth-moment estimation, which is generally more complex than estimating the covariance. This difference highlights a key advantage of employing the operator norm for covariance testing.

Third, we perform size and power analyses for the statistics $T$ and $T^{\text{Roy}}$ and propose a combined approach to enhance testing power. We also develop a generalized universality theorem to show the consistency of the universal bootstrap for statistics based on extreme eigenvalues, including this combined statistic. For size analysis, we demonstrate that the universality results apply to $T^{\text{Roy}}$. As a byproduct, we extend the Tracy–Widom law for the largest eigenvalue of sample covariance with general entries under the ultra-high dimension regime, addressing a long-standing gap since it was only proven for Gaussian entries by Karoui (2003). For power performance, we analyze both $T$ and $T^{\text{Roy}}$ within the generalized spiked model framework from Bai and Yao (2012); Jiang and Bai (2021). We show that $T^{\text{Roy}}$ performs better in worst-case scenarios while $T$ excels on average performance. To further enhance the power, we propose a new combined statistic, $T^{\text{Com}}$, supported by a generalized universality theorem that confirms the validity of the universal bootstrap for $T^{\text{Com}}$. This also serves as a theoretical guarantee for the universal bootstrap applicable to a broader range of other statistics based on extreme eigenvalues. Extensive simulations validate these findings, highlighting the superior performance of our combined statistics across a range of scenarios.

Throughout the paper, we reserve boldfaced symbols for vectors and matrices. For a complex number $z\in\mathbb{C}$, we use $\Re z$ and $\Im z$ for its real part and imaginary part, respectively. For a vector $\bm{u}\in\mathbb{R}^{p}$, we use $\|\bm{u}\|=\sqrt{\sum_{i=1}^{p}u_{i}^{2}}$ for its Euclidean norm. For a matrix $\bm{A}=(a_{ij})_{M\times N}$, we use $\|\bm{A}\|_{\text{op}}$, $\|\bm{A}\|_{\text{F}}=:\sqrt{\sum_{i=1}^{M}\sum_{j=1}^{N}a_{ij}^{2}}$, and $\|\bm{A}\|_{\text{sup}}=:\sup_{1\leq i\leq M,1\leq j\leq N}|a_{ij}|$ to represent its operator norm, Frobenius norm and the supremum norm, respectively. Denote the singular values of $\bm{A}$ by $\sigma_{1}(\bm{A})\geq\sigma_{2}(\bm{A})\geq\cdots\geq\sigma_{r}(\bm{A})$, where $r=\min\{M,N\}$. When $M=N$, we use $\lambda_{1}(\bm{A}),\lambda_{2}(\bm{A}),\cdots\lambda_{M}(\bm{A})$ for the eigenvalues of $\bm{A}$. The trace of a square matrix $\bm{A}$ is denoted by $\text{tr}(\bm{A})=\sum_{i=1}^{p}a_{ii}$. For two sequences $\{a_{n}\}_{n=1}^{\infty}$, $\{b_{n}\}_{n=1}^{\infty}$, we use $a_{n}\lesssim b_{n}$ or $a_{n}=O(b_{n})$ to show there exists a constant $C$ not depending on $n$ such that $|a_{n}|\leq Cb_{n}$ for all $n\in\mathbb{N}$. Denote $a_{n}=o(b_{n})$ if $a_{n}/b_{n}\to 0$, and $a_{n}\asymp b_{n}$ if both $a_{n}\lesssim b_{n}$ and $b_{n}\lesssim a_{n}$.

The paper is organized as follows. Section 2 provides an overview of the proposed method and presents our main results on the consistency of the proposed universal bootstrap procedure. Section 3 describes the key tools from random matrix theory and applies the universality theorem to covariance testing problems. Section 4 presents the power analysis of the universal bootstrap procedure. Section 5 provides simulation results and a real data example, demonstrating the numerical performance of our operator norm-based statistics. Detailed technical proofs are provided in the supplementary material. The data and codes are publicly available in a GitHub repository (https://github.com/zhang-guoyu/universal_bootstrap).

Our results rely on the universality from the random matrix theory. To proceed, we first introduce some preliminary results relevant to our analysis. Denote $\bm{X}_{i}=\bm{\Sigma}^{1/2}\bm{Z}_{i}$, where $\bm{Z}_{i}$ has zero mean and identity covariance matrix for $i=1,\cdots,n$. We accordingly define $\bm{Z}=(\bm{Z}_{1},\cdots,\bm{Z}_{n})^{T}\in\mathbb{R}^{n\times p}$. The primary matrix of interest is $\bm{\hat{\Sigma}}-\bm{\Sigma}_{0}$. This inspires us to consider the matrix of a general form

where $\bm{M}^{\prime}_{n}=(np)^{-1/2}\ \bm{\Sigma}^{1/2}\bm{Z}^{T}\bm{Z}\bm{\Sigma}^{1/2}+\bm{R}^{\prime}$ and $\bm{R}^{\prime}=\phi^{-1/2}\ \bm{R}$. Here, $\bm{R}$ is a symmetric matrix with eigenvalues that may be positive, negative, or zero, and $\bm{M}^{\prime}$ is normalized to address cases where $\alpha>1$, as per Ding and Wang (2023). We impose some assumptions to determine the limit of the Stieltjes transform $m_{\bm{M}^{\prime}_{n}}(z)=\text{tr}(\bm{M}_{n}^{\prime}-z\bm{I}_{p})^{-1}/p$ of $\bm{M}^{\prime}_{n}$.

Assumptions 1 and 2 are standard and often appear in the random matrix literature. Notably, Assumption 1 relies only on independence, without assuming identical distribution as in Qiu, Li and Yao (2023). While all moments exist under Assumption 1, this condition could be relaxed as discussed in Ding and Yang (2018). We do not pursue this here. We also permit $\bm{\Sigma}$ to be singular, a less restrictive condition than the invertibility assumed in Ding and Wang (2023); Ding, Hu and Wang (2024). This generalization enables testing singular $\bm{\Sigma}_{0}$, where $T^{\text{Roy}}$ is invalid but $T$ remains applicable. Lastly, Assumption 3, imposed for technical requirements, necessitates that $\bm{\Sigma}$ and $\bm{R}$ share the same eigenvectors. The same Assumption is required in the signal-plus-noise model as in Zhou, Bai and Hu (2023); Zhou et al. (2024).

We introduce the following deterministic equivalence matrix

where $\phi=p/n$ and $e(z)$ is the fixed point of the equation $e(z)=\text{tr}\big{(}\bm{R}^{\prime}+\phi^{-1/2}(1+\phi^{1/2}e(z))^{-1}\bm{\Sigma}-z\bm{I}_{p}\big{)}^{-1}\bm{\Sigma}/p$. Using $\bar{\bm{Q}}_{p}(z)$, we define the associated Stieltjes transform $\bar{m}_{p}(z)$ as $\bar{m}_{p}(z)=\text{tr}(\bar{\bm{Q}}_{p}(z))/p$. According to Couillet, Debbah and Silverstein (2011), $e(z)$ and $\bar{\bm{Q}}_{p}(z)$ are well-defined for every $z\in\mathbb{C}^{+}$, and $\bar{m}_{p}(z)$ serves as the Stieltjes transform of a measure $\rho$ on $\mathbb{R}$. We further denote $E_{+}$ and $E_{-}$ as the endpoints of $\rho$. Formally, if we define $\text{supp}(\rho)=\left\{x\in\mathbb{R}\ :\ \rho([x-\epsilon,x+\epsilon])>0,\ \forall\epsilon>0\right\}$, we have $E_{+}=\sup\text{supp}(\rho)$ and $E_{-}=\inf\text{supp}(\rho)$. Combining results of Knowles and Yin (2017) for the $\alpha=1$ case and Ding and Wang (2023) for the $\alpha>1$ case, it follows that $E_{+}\asymp\phi^{1/2}$ and $|E_{-}|\lesssim 1$. Intuitively, the largest eigenvalue $\lambda_{1}(\bm{M}^{\prime}_{n})$ approaches $E_{+}$, while $\lambda_{p}(\bm{M}^{\prime}_{n})$ approaches $E_{-}$. The next subsection characterizes the fluctuations of $\lambda_{1}(\bm{M}^{\prime}_{n})-E_{+}$ and $\lambda_{p}(\bm{M}^{\prime}_{n})-E_{-}$.

In this subsection, we establish the anisotropic local law, which forms the foundation for our universality results. Specifically, we aim to describe the fluctuations of $\lambda_{1}(\bm{M}^{\prime}_{n})-E_{+}$ and $\lambda_{p}(\bm{M}^{\prime}_{n})-E_{-}$. This requires us to analyze the convergence of the Stieltjes transform near the endpoints $E_{+}$ and $E_{-}$. We define the local domains $\mathcal{D}_{+}$ and $\mathcal{D}_{-}$ as

for a fixed parameter $0<\tau<1$. For $\alpha=1$, let $\mathcal{D}=\mathcal{D}_{+}\cup\mathcal{D}_{-}$, and for $\alpha>1$, define $\mathcal{D}=\mathcal{D}_{+}$. We will focus on the behavior of the Green function on $\mathcal{D}$. This definition ensures that for $\alpha=1$, both the largest and smallest eigenvalues of $\bm{M}_{n}^{\prime}$ are controlled, so that $\sigma_{1}(\bm{M}_{n}^{\prime})=\max\{|\lambda_{1}(\bm{M}_{n}^{\prime})|,|\lambda_{p}(\bm{M}_{n}^{\prime})|\}$ can also be controlled. While for the case $\alpha>1$, we always have $|E_{+}|\gg|E_{-}|$, thus only the largest eigenvalue requires attention.

To facilitate the presentation of the anisotropic local law, we assume in this subsection that

Under (17), we can rewrite $\bm{M}^{\prime}_{n}$ as $\bm{M}^{\prime}_{n}=\bm{\Sigma}^{1/2}\left((np)^{-1/2}\ \bm{Z}^{T}\bm{Z}+\bm{R}^{{}^{\prime\prime}}\right)\bm{\Sigma}^{1/2}$ where $\bm{R}^{{}^{\prime\prime}}=\bm{\Sigma}^{-1/2}\bm{R}^{{}^{\prime}}\bm{\Sigma}^{-1/2}$. Since $\bm{M}^{\prime}_{n}$ is quadratic in $\bm{X}$, we follow Knowles and Yin (2017) to define the linearization matrix $\bm{H}(z)$ and the corresponding linearized Green function $\bm{G}(z)$

Notice that $\bm{G}(z)$ is a large $(n+2p)\times(n+2p)$ matrix. We also define the deterministic equivalent Green function as

where $\tilde{m}(z)=-z^{-1}(1+\phi^{1/2}e(z))^{-1}$. The anisotropic local law aims to control the difference $\bm{G}(z)-\bar{\bm{G}}(z)$ for $z\in\mathcal{D}$. A crucial observation is that defining the parameter $z$-dependent covariance $\bm{\Sigma}(z)=z\bm{\Sigma}(z\bm{I}_{p}-\bm{R}^{\prime})^{-1}$ with the corresponding measure $\pi^{z}=\sum_{i=1}^{p}\delta_{\lambda_{i}(\bm{\Sigma}(z))}/p$, where $\delta_{z}$ is the Dirac point measure at $z$, we have the following $z$-dependent deformed Marčenko-Pastur law,

This result mirrors the form of the deformed Marčenko-Pastur law in Ding and Wang (2023), with $\bm{\Sigma}(z)$ as the covariance matrix. This demonstrates that the effect of $\bm{R}$ can be represented by turning $\bm{\Sigma}$ into the $z$-dependent covariance $\bm{\Sigma}(z)$. This insight simplifies our proof and presentation. For the denominator in (26), we impose the following technical assumption.

We provide several remarks regarding Assumption 4. Informally, condition (27) ensures that the extreme eigenvalues of $\bm{\Sigma}$ do not spread near the endpoints $E_{\pm}$, thereby preventing spikes outside the support of $\rho$. Similar assumptions have been made in the literature for the universality of $\hat{\bm{\Sigma}}$, using the Stieltjes transform in place of $\tilde{m}$ and $\lambda_{i}(\bm{\Sigma})$ in place of $\lambda_{i}(\bm{\Sigma}(E_{\pm}))$, as in Bao, Pan and Zhou (2015); Knowles and Yin (2017). This aligns with the intuition that the effect of $\bm{R}$ can be expressed through the transformation from $\bm{\Sigma}$ to $\bm{\Sigma}(z)$. Moreover, (27) is automatically satisfied when $\phi=p/n\to\infty$, i.e. $\alpha>1$. We thus impose Assumption 4 only in the case $\alpha=1$.

To state our main results, we define the concept of stochastic domination, introduced by Erdős, Knowles and Yau (2013) and widely applied in random matrix theory (Knowles and Yin, 2017). For two families of non-negative random variables $A=\left(A^{(n)}(u)\ :\ n\in\mathbb{N},u\in U^{(n)}\right)$ and $B=\left(B^{(n)}(u)\ :\ n\in\mathbb{N},u\in U^{(n)}\right)$, where $U^{(n)}$ is $n$-dependent parameter set, we say $A$ is stochastically dominated by $B$ uniformly in $u$ if for any $\epsilon>0$ and $D>0$, we have $\sup_{u\in U^{(n)}}\mathbb{P}\left(A^{(n)}(u)>n^{\epsilon}B^{(n)}(u)\right)\leq n^{-D}$ for large enough $n\geq n_{0}(\epsilon,D)$. We use the notation $A\prec B$ to represent this relationship. When $A$ is a family of negative or complex random variables, we also write $A\prec B$ or $A=O_{\prec}(B)$ to indicate $|A|\prec B$. With these definitions, we state the anisotropic local law as follows.

The anisotropic local law provides a delicate characterization of the discrepancy between the random Green function and its deterministic counterpart, yielding a precise bound to analyze the behavior of the extreme eigenvalues of $\bm{M}_{n}^{\prime}$. When $\bm{R}=\bm{0}$, similar results have been established for $\alpha=1$ in Knowles and Yin (2017) and for $\alpha\geq 1$ in Ding and Wang (2023); Ding, Hu and Wang (2024). As demonstrated in Ding and Wang (2023), the results for general $\alpha\geq 1$ hold without requiring the dimension $p$ and sample size $n$ to grow proportionally, resulting in different convergence rates for the blocks of $\bm{G}(z)$. Ding and Wang (2023) separately provides convergence rates for each block and offers a coarse bound for the anisotropic local law

In this study, we enhance these findings by introducing the weight matrix $\bm{\Phi}$, allowing us to express the anisotropic local law (28) in a more compact form. This reformulation refines the convergence rate and simplifies our proof.

Building on the anisotropic local law theorem, we proceed to establish our universality result. Our first step involves deriving key implications from Theorem 3.1. Recall that $\rho$ is the measure on $\mathbb{R}$ whose Stieltjes transform is the deterministic equivalence $\bar{m}_{p}(z)$. We define the quantile sequence of $\rho$ as $w_{1}\geq\cdots\geq w_{p}$ such that $\int_{w_{i}}^{+\infty}\rho(x)\mathrm{d}x=(i-1/2)/p$ for $i=1,\cdots,p$. We aim to demonstrate that the eigenvalues $\lambda_{i}(\bm{M}_{n}^{\prime})$ of $\bm{M}_{n}^{\prime}$ remain close to $w_{i}$ for $i=1,\cdots,p$.

Theorem 3.2 establishes the rigidity of the eigenvalues of $\bm{M}_{n}^{\prime}$, with two main implications. First, when $\alpha=1$, equation (30) shows that the largest and smallest $k$-th eigenvalues lie within an $n^{-2/3+\epsilon}$-neighborhood of $w_{p}$, respectively, for any $\epsilon>0$. Furthermore, we have $|w_{1}-E_{+}|,|w_{p}-E_{-}|\prec n^{-2/3}$, leading to $|\lambda_{i}(\bm{M}_{n}^{\prime})-E_{+}|,|\lambda_{p-i}(\bm{M}_{n}^{\prime})-E_{-}|\prec n^{-2/3}$ uniformly for $1\leq i\leq k$. Second, with Theorem 3.2, we demonstrate that $w_{i}\asymp\phi^{1/2}\to+\infty$ for $1\leq i\leq n$, revealing a fast $n^{-2/3+\epsilon}$ approximation rate for the diverging quantities $\lambda_{i}(\bm{M}_{n}^{\prime})$ and $w_{i}$.

To establish the universality result, we provide bounds for the discrepancy between the distribution of $\lambda_{1}(\bm{M}_{n}^{\prime})$ and its Gaussian counterpart. We denote by $\mathbb{P}^{\text{Gau}}$ and $\mathbb{E}^{\text{Gau}}$ the probability and expectation under the additional assumption that $\{Z_{ij}\}_{1\leq i\leq n,1\leq j\leq p}$ are independent standard Gaussian variables. With this notation, we present the following result.