A Bootstrap Method for Spectral Statistics in High-Dimensional
Elliptical Models

At an algorithmic level, the proposed method is built on top of two estimators. The first is an estimator $\widehat{\varsigma}_{n}^{\,2}$ for the variance parameter $\varsigma_{n}^{2}=\operatorname{var}(\xi_{1}^{2})$. The second is an estimator $\tilde{\Sigma}_{n}$ for the diagonal matrix of population eigenvalues $\textup{diag}(\lambda_{1}(\Sigma_{n}),\dots,\lambda_{p}(\Sigma_{n}))$. Once these two estimators have been assembled, the method generates bootstrap data from an elliptical model parameterized in terms of $\widehat{\varsigma}_{n}^{\,2}$ and $\tilde{\Sigma}_{n}$. More specifically, the $i$th sample of each bootstrap dataset $\mathbf{x}_{1}^{*},\dots,\mathbf{x}_{n}^{*}$ is generated to be of the form

where $\xi_{i}^{*}$ is a non-negative random variable satisfying $\mathbb{E}((\xi_{i}^{*})^{2}|X)=p$ and $\operatorname{var}((\xi_{i}^{*})^{2}|X)=\widehat{\varsigma}_{n}^{\,2}$, and $\mathbf{u}_{i}^{*}\in\mathbb{R}^{p}$ is drawn uniformly from the unit sphere, independently of $\xi_{i}^{*}$. Then, a single bootstrap sample of a generic spectral statistic $T_{n}=\psi(\lambda_{1}(\widehat{\Sigma}_{n}),\dots,\lambda_{p}(\widehat{\Sigma}_{n}))$ is computed as $T_{n}^{*}=\psi(\lambda_{1}(\widehat{\Sigma}_{n}^{*}),\dots,\lambda_{p}(\widehat{\Sigma}_{n}^{*}))$, where $\widehat{\Sigma}_{n}^{*}$ is the sample covariance matrix of the bootstrap data.

To emphasize the modular role that $\widehat{\varsigma}_{n}^{\,2}$ and $\tilde{\Sigma}_{n}$ play in the bootstrap sampling process, we will provide the details for their construction later in Sections 2.2 and 2.3. With these points understood, the following algorithm shows that the method is very easy to implement.

The estimation of the parameter $\varsigma_{n}^{2}=\operatorname{var}(\xi_{1}^{2})$ is complicated by the fact that the random variables $\xi_{1},\dots,\xi_{n}$ are not directly observable in an elliptical model. It is possible to overcome this challenge with the following estimating equation, which can be derived from an explicit formula for $\operatorname{var}(\|\mathbf{x}_{1}\|_{2}^{2})$ that is given in Lemma D.1,

Based on this equation, our approach is to separately estimate each of the three moment parameters on the right hand side, denoted as

These parameters have the advantage that they can be estimated in a more direct manner, due to their simpler relations with the observations and the matrix $\Sigma_{n}$. Specifically, we use estimates defined according to

Substituting these estimates into (2.1) yields our proposed estimate for $\varsigma_{n}^{2}$

where $x_{+}=\max\{x,0\}$ denotes the non-negative part of any real number $x$. The consistency of this estimate will be established in Theorem 1, which shows $\frac{1}{p}(\widehat{\varsigma}_{n}^{\,2}-\varsigma_{n}^{2})\to 0$ in probability as $n\to\infty$.

The problem of estimating the eigenvalues of a population covariance matrix has attracted long-term interest in the high-dimensional statistics literature, and many different estimation methods have been proposed  [e.g. 12, 38, 7, 26, 25]. In order to estimate $\lambda_{1}(\Sigma_{n}),\dots,\lambda_{p}(\Sigma_{n})$ in our current setting, we modify the method of QuEST [26], which has become widely used in recent years. This choice has the benefit of making all aspects of our proposed method easy to implement, because QuEST is supported by a turnkey software package [27].

We denote the estimates produced by QuEST as $\widehat{\lambda}_{\textup{Q},1}\geq\dots\geq\widehat{\lambda}_{\textup{Q},p}$, and we will use modified versions of them defined by

for $j=1,\dots,p$, where we let $\widehat{b}_{n}=\lambda_{1}(\widehat{\Sigma}_{n})+1$. The modification is done for theoretical reasons, to ensure that $\tilde{\lambda}_{1},\dots,\tilde{\lambda}_{p}$ are asymptotically bounded, which follows from Lemma A.5. In addition, we define the diagonal $p\times p$ matrix associated with these estimates as

Later, in Theorem 1, we will show that the estimates $\tilde{\lambda}_{1},\dots,\tilde{\lambda}_{p}$ are consistent, in the sense that their empirical distribution converges weakly in probability to the correct limit as $n\to\infty$.

All of our theoretical analysis is framed in terms of a sequence of models indexed by $n$, so that all model parameters are allowed to vary with $n$, except when stated otherwise. The details of our model assumptions are given below in Assumptions 1 and 2.

In addition to Assumption 1, we need one more assumption dealing with the spectrum of the population covariance matrix $\Sigma_{n}$. To state this assumption, let $H_{n}$ denote the empirical distribution function associated with $\lambda_{1}(\Sigma_{n}),\dots,\lambda_{p}(\Sigma_{n})$, which is defined for any $t\in\mathbb{R}$ according to

Here, we establish the consistency of the estimators $\widehat{\varsigma}_{n}^{\,2}$ and $\tilde{\Sigma}_{n}$, defined in (2.2) and (2.4). The appropriate notion of consistency for $\tilde{\Sigma}_{n}$ is stated in terms of its empirical spectral distribution function, which is defined for any $t\in\mathbb{R}$ as

To develop our main result on the consistency of the proposed bootstrap method, it is necessary to recall some background facts and introduce several pieces of notation.

Under Assumptions 1 and 2, it is known from [6, Theorem 1.1] that an extended version of the classical Marčenko-Pastur Theorem holds for the empirical spectral distribution function $\widehat{H}_{n}(t)=\frac{1}{p}\sum_{j=1}^{p}1\{\lambda_{j}(\widehat{\Sigma}_{n})\leq t\}$. Namely, there is a probability distribution $\Psi(H,c)$ on $[0,\infty)$, depending only on $H$ and $c$, such that the weak limit $\widehat{H}_{n}\Rightarrow\Psi(H,c)$ occurs almost surely. In this statement, we may regard $\Psi(\cdot,\cdot)$ as a map that takes a distribution $H^{\prime}$ on $[0,\infty)$ and a number $c^{\prime}>0$ as input, and returns a new distribution $\Psi(H^{\prime}\!,c^{\prime})$ on $[0,\infty)$ whose Stieltjes transform $m_{H^{\prime}\!,c^{\prime}}(z)=\int\frac{1}{\lambda-z}d(\Psi(H^{\prime}\!,c^{\prime}))(\lambda)$ solves the Marčenko-Pastur equation (3.9) below. That is, for any $z\in\mathbb{C}^{+}$, the number $m_{H^{\prime}\!,c^{\prime}}=m_{H^{\prime}\!,c^{\prime}}(z)$ is the unique solution to the equation

within the set $\{m\in\mathbb{C}\,|\,(c^{\prime}-1)/z+c^{\prime}m\in\mathbb{C}^{+}\}$.

The map $\Psi(\cdot,\cdot)$ is relevant to our purposes here, because it determines a centering parameter that is commonly used in limit theorems for linear spectral statistics. In detail, if $H_{n,c_{n}}$ denotes a shorthand for the probability distribution $\Psi(H_{n},c_{n})$, and if $f$ is a real-valued function defined on the support of $H_{n,c_{n}}$, then the associated centering parameter is defined as

Similarly, let $\tilde{H}_{n,c_{n}}=\Psi(\tilde{H}_{n},c_{n})$ and $\tilde{\vartheta}_{n}(f)=\int f(t)d\tilde{H}_{n,c_{n}}(t)$. Also, in order to simplify notation for handling the joint distribution of several linear spectral statistics arising from a fixed set of functions $\mathbf{f}=(f_{1},\dots,f_{k})$, we write $T_{n}(\mathbf{f})=(T_{n}(f_{1}),\dots,T_{n}(f_{k}))$, and likewise for $T_{n,1}^{*}(\mathbf{f})$, $\vartheta_{n}(\mathbf{f})$, and $\tilde{\vartheta}_{n}(\mathbf{f})$.

As one more preparatory item, recall from page 1 that $d_{\textup{LP}}$ denotes the Lévy-Prohorov metric for comparing distributions on $\mathbb{R}^{k}$. This metric is a standard choice for formulating bootstrap consistency results, as it has the fundamental property of metrizing weak convergence.

When dealing with high-dimensional covariance matrices, it is often of interest to have a measure of the number of “dominant eigenvalues”. One such measure is the stable rank, defined as

which arises naturally in a plethora of situations [e.g. 4, 43, 44, 30, 36, 33]. Whenever $\Sigma_{n}$ is non-zero, this parameter satisfies $1\leq r_{n}\leq\textup{rank}(\Sigma_{n})$, and the equality $r_{n}=p$ holds if and only if $\Sigma_{n}$ is proportional to the identity matrix.

In this subsection, we illustrate how the proposed bootstrap method can be applied to solve some inference problems involving the parameter $r_{n}$. Our first example shows how to construct a confidence interval for $r_{n}$, and our subsequent examples deal with testing procedures related to $r_{n}$. Later, in Theorem 3, we establish the theoretical validity of the methods used in these examples—showing that the confidence interval has asymptotically exact coverage, and that the relevant testing procedures maintain asymptotic control of their levels.

In our simulations, we generated data from elliptical models that were parameterized as follows. The random variable $\xi_{1}^{2}$ was generated using four choices of distributions:

1. (i).

   Chi-Squared distribution with $p$ degrees of freedom
   

2. (ii).

   Beta-Prime$\left(\frac{p(1+p+8)}{8},\frac{1+p+16}{8}\right)$,
   

3. (iii).

   $(p+4)\textup{Beta}(p/2,2)$,
   

4. (iv).

   $\frac{9p}{10}\textup{F}(p,20)$,

where F$(d_{1},d_{2})$ denotes an F-distribution with $d_{1}$ and $d_{2}$ degrees of freedom. Note that cases (i), (iii), (iv), correspond respectively to a multivariate Gaussian distribution, a multivariate Pearson type II distribution, and a multivariate t-distribution with 20 degrees of freedom. Also, the numerical values appearing in (i)-(iv) were chosen to ensure the normalization condition $\mathbb{E}(\xi_{1}^{2})=p$.

The population covariance matrix $\Sigma_{n}$ was selected from five options:

1. 1.

   The eigenvalues of $\Sigma_{n}$ are $\lambda_{1}(\Sigma_{n})=\cdots=\lambda_{5}(\Sigma_{n})=\frac{4}{3}$ and $\lambda_{j}(\Sigma_{n})=1$ for $j\in\{6,\dots,p\}$. The $p\times p$ matrix of eigenvectors of $\Sigma_{n}$ is generated from the uniform distribution on orthogonal matrices.
   

2. 2.

   The eigenvalues of $\Sigma_{n}$ are $\lambda_{j}(\Sigma_{n})=\exp(-j/4)$ for $j\in\{1,\dots,20\}$, and $\lambda_{20}(\Sigma_{n})=\cdots=\lambda_{p}(\Sigma_{n})$. The eigenvectors are the same as in case 1.
   

3. 3.

   The matrix $\Sigma_{n}$ has entries of the form $(\Sigma_{n})_{ij}=(\frac{1}{10})^{|i-j|}+1\{i=j\}$.
   

4. 4.

   The matrix $\Sigma_{n}$ has entries of the form $(\Sigma_{n})_{ij}=(\frac{1}{10})1\{i\neq j\}+1\{i=j\}$.
   

5. 5.

   The eigenvalues of $\Sigma_{n}$ are $\lambda_{1}(\Sigma_{n})=5$ and $\lambda_{j}(\Sigma_{n})=1$ for $j\in\{2,\dots,p\}$. The eigenvectors are the same as in case 1.

Our experiments for linear spectral statistics were based on the task of using the proposed bootstrap to estimate three parameters of $\mathcal{L}(p(T_{n}(f)-\vartheta_{n}(f)))$: the mean, standard deviation, and 95th percentile. We considered two choices for the function $f$, namely $f(x)=x^{2}$ and $f(x)=x-\log(x)-1$. In the first case, we selected the ratio $p/n$ so that $p/n\in\{0.5,1,1.5\}$, and in the second case we used $p/n\in\{0.3,0.5,0.7\}$.

With regard to the bootstrap, we ran Algorithm 1 on the first 500 datasets corresponding to each parameter setting. Also, we generated $B=250$ bootstrap samples of the form $p(T_{n,1}^{*}(f)-\tilde{\vartheta}_{n}(f))$ during every run. As a result of these runs, we obtained 500 different bootstrap estimates of the mean, standard deviation, and 95th percentile of $\mathcal{L}(p(T_{n}(f)-\vartheta_{n}(f)))$. In Tables 1 and 2, we report the empirical mean and standard deviation (in parenthesis) of these 500 estimates in the second row of numbers corresponding to each choice of $\Sigma_{n}$.

One more detail to mention is related to the computation of $\vartheta_{n}(f)$ and $\tilde{\vartheta}_{n}(f)$. For each parameter setting, we approximated $\vartheta_{n}(f)$ as follows. We averaged 30 realizations of $\frac{1}{40p}\sum_{j=1}^{40p}f(\lambda_{j}(\frac{1}{40n}\sum_{i=1}^{40n}\xi_{i}^{2}\mathsf{\Sigma}_{n}^{1/2}\mathbf{u}_{i}\mathbf{u}_{i}^{\top}\mathsf{\Sigma}_{n}^{1/2}))$, where $\mathsf{\Sigma}_{n}=I_{40}\otimes\Sigma_{n}$ is of size $40p\times 40p$, each $\mathbf{u}_{i}$ was drawn from the uniform distribution on the unit sphere of $\mathbb{R}^{40p}$, and each $\xi_{i}^{2}$ was generated as in (i)-(iv), but with $40p$ replacing $p$. For the bootstrap samples, we computed one realization of the statistic $\frac{1}{40p}\sum_{j=1}^{40p}f(\lambda_{j}(\frac{1}{40n}\sum_{i=1}^{40n}(\xi_{i}^{*})^{2}\tilde{\mathsf{\Sigma}}_{n}^{1/2}\mathbf{u}_{i}^{*}(\mathbf{u}_{i}^{*})^{\top}\tilde{\mathsf{\Sigma}}_{n}^{1/2}))$ to approximate $\tilde{\vartheta}_{n}(f)$ during every run of Algorithm 1, where $\tilde{\mathsf{\Sigma}}_{n}=I_{40}\otimes\tilde{\Sigma}_{n}$ is of size $40p\times 40p$, each $\mathbf{u}_{i}^{*}$ was drawn from the uniform distribution on the unit sphere of $\mathbb{R}^{40p}$, and each $(\xi_{i}^{*})^{2}$ was drawn from a Gamma distribution with mean $40p$ and variance $40\widehat{\varsigma}_{n}^{2}$.