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

Siyao<EMAIL_ADDRESS>[ Miles E.
<EMAIL_ADDRESS>[ Department of Statistics, University of
California, Davis

(2020)

###### Abstract

Although there is an extensive literature on the eigenvalues of high-
dimensional sample covariance matrices, much of it is specialized to Marčenko-
Pastur (MP) models—in which observations are represented as linear
transformations of random vectors with independent entries. By contrast, less
is known in the context of elliptical models, which violate the independence
structure of MP models and exhibit quite different statistical phenomena. In
particular, very little is known about the scope of bootstrap methods for
doing inference with spectral statistics in high-dimensional elliptical
models. To fill this gap, we show how a bootstrap approach developed
previously for MP models can be extended to handle the different properties of
elliptical models. Within this setting, our main theoretical result guarantees
that the proposed method consistently approximates the distributions of linear
spectral statistics, which play a fundamental role in multivariate analysis.
Lastly, we provide empirical results showing that the proposed method also
performs well for a variety of nonlinear spectral statistics.

62H99 ,

62F40,

62H25,

high-dimensional statistics,

bootstrap,

covariance matrices,

elliptical distributions,

###### keywords:

[class=MSC]

###### keywords:

††volume: 0††issue: 0

and

t1Supported in part by NSF Grant DMS 1915786

##  1 Introduction

The analysis of spectral statistics of sample covariance matrices is a major
research area within multivariate analysis, random matrix theory, high-
dimensional statistics, and related fields [5, 45, 43, 38]. If
$\mathbf{x}_{1},\dots,\mathbf{x}_{n}$ are centered i.i.d. observations in
$\mathbb{R}^{p}$ with a sample covariance matrix denoted by

$\widehat{\Sigma}_{n}=\frac{1}{n}\sum_{i=1}^{n}\mathbf{x}_{i}\mathbf{x}_{i}^{\top},$
(1.1)

then we say that a random variable $T_{n}$ is a spectral statistic if it has
the form
$T_{n}=\psi(\lambda_{1}(\widehat{\Sigma}_{n}),\dots,\lambda_{p}(\widehat{\Sigma}_{n}))$,
where
$\lambda_{1}(\widehat{\Sigma}_{n})\geq\cdots\geq\lambda_{p}(\widehat{\Sigma}_{n})$
are the sorted eigenvalues of $\widehat{\Sigma}_{n}$, and $\psi$ is a generic
real-valued function. Over the past two decades, there has been a tremendous
growth of interest in the behavior of spectral statistics in high-dimensional
settings where $p=p(n)$ grows so that $p/n$ converges to a positive constant
as $n\to\infty$. Likewise, statistical tools for approximating the
distributions of spectral statistics have been applied to high-dimensional
data in a broad range of domains, such as electrical engineering, finance, and
biology [37, 1, 9, 10].

Although the research on spectral statistics has dealt with many different
statistical models, two of the most influential ones have been _elliptical_
and _Marčenko-Pastur_ (MP) models. To be specific, we say that the random
vector $\mathbf{x}_{1}$ follows an elliptical model if it can be represented
in the form

$\mathbf{x}_{1}=\xi_{1}\Sigma_{n}^{1/2}\mathbf{u}_{1}$ (1.2)

where $\Sigma_{n}=\mathbb{E}(\mathbf{x}_{1}\mathbf{x}_{1}^{\top})$ is the
$p\times p$ population covariance, and
$(\mathbf{u}_{1},\xi_{1})\in\mathbb{R}^{p}\times[0,\infty)$ is a random vector
such that $\mathbf{u}_{1}$ and $\xi_{1}$ are independent, with
$\mathbf{u}_{1}$ being uniformly distributed on the unit sphere of
$\mathbb{R}^{p}$. Alternatively, we say that $\mathbf{x}_{1}$ follows an MP
model if

$\mathbf{x}_{1}=\Sigma_{n}^{1/2}\mathbf{z}_{1},$ (1.3)

where $\mathbf{z}_{1}\in\mathbb{R}^{p}$ is a random vector whose entries are
centered and independent.

At first sight, these two models may seem to be very similar, but this outward
appearance conceals some crucial differences in modelling capabilities. In
particular, it should be stressed that the entries of the random vector
$\xi_{1}\mathbf{u}_{1}$ in an elliptical model are correlated, which contrasts
with the independence of the entries of $\mathbf{z}_{1}$ in an MP model. Also,
since the scalar random variable $\xi_{1}$ is shared across all entries of
$\mathbf{x}_{1}$ in an elliptical model, this enhances the ability to capture
scenarios where the magnitudes of all entries of $\mathbf{x}_{1}$ move in the
same direction simultaneously. This is a key effect in some application
domains, such as in finance, where the entries of $\mathbf{x}_{1}$ correspond
to stock prices that can fall in tandem during a sharp market downturn.
Additional background on related merits of elliptical models can be found in
[11, 35, 15]. More generally, the multivariate analysis literature has placed
a longstanding emphasis on the benefits of elliptical models in fitting
various types of non-Gaussian data [13, 2, 16].

However, looking beyond the points just mentioned, MP models have played a
more dominant role than elliptical models in the literature on spectral
statistics in high dimensions. Consequently, the established body of high-
dimensional limit theory is much less complete for elliptical models. Indeed,
the challenge of extending results from MP models to elliptical ones has
become a prominent topic of ongoing research, which has led to important
advances in the limit theory for spectral statistics [e.g. 27, 18, 17, 44, 23,
28, 47]. As a matter of historical context, it also worth bearing in mind that
for some spectral statistics, it took many years for such extensions to be
established.

From the standpoint of statistical methodology, a corresponding set of gaps
exists between elliptical and MP models. These gaps are especially apparent in
the current state of bootstrap methods for high-dimensional data. In
particular, it is known from [32] that a form of parametric bootstrapping can
successfully approximate the distributions of spectral statistics in MP
models, whereas very little is known for elliptical models. Accordingly, our
primary goal in the current paper is to resolve this issue by developing a
parametric bootstrap method that is both theoretically and empirically
effective in high-dimensional elliptical models.

With regard to theory, we will focus on the class of _linear spectral
statistics_ , which have the form

$T_{n}(f)=\frac{1}{p}\sum_{j=1}^{p}f(\lambda_{j}(\widehat{\Sigma}_{n})),$
(1.4)

for a suitable real-valued function $f$. Beginning with the pathbreaking works
[22, 3] that established the earliest versions of the central limit theorem
for linear spectral statistics in high dimensions, these statistics have been
a perennial focus of research. Their importance is underscored by the fact
that they appear frequently throughout multivariate analysis, with some of the
most well-known examples being $\frac{1}{p}\textup{tr}(\widehat{\Sigma}_{n})$,
$\frac{1}{p}\textup{tr}(\widehat{\Sigma}_{n}^{2})$, and
$\frac{1}{p}\log\det(\widehat{\Sigma}_{n})$, among various other classical
statistics for testing hypotheses [45].

Motivated by these considerations, our main theoretical result (Theorem 2)
shows that the proposed bootstrap method consistently approximates the
distributions of linear spectral statistics when the underlying data are
elliptical and $p/n$ converges to a positive constant as $n\to\infty$. The
proof substantially leverages recent progress on the central limit theorem for
linear spectral statistics in elliptical models due to [17]. Also, an
intermediate step in the proof (Lemma A.4) shows that the well-known
eigenvalue estimation method QuEST is consistent in elliptical models—which
may be of independent interest, since it seems that QuEST’s consistency has
not previously been reported outside of MP models [25]. Moreover, Section 3.4
develops an application of Theorem 2 where we establish the asymptotic
validity of inference procedures related to the _stable rank_ parameter
$r_{n}=\textup{tr}(\Sigma_{n})^{2}/\textup{tr}(\Sigma_{n}^{2})$ of the
population covariance matrix
$\Sigma_{n}=\mathbb{E}(\mathbf{x}_{1}\mathbf{x}_{1}^{\top})$.

To address the empirical performance of the proposed method, Section 4
presents numerical results for a variety of model settings and statistics.
Most notably, these results show encouraging performance for _both linear and
nonlinear spectral statistics_. To put this point into perspective, it is
important to highlight the fact that asymptotic formulas for nonlinear
spectral statistics are typically developed on a case-by-case basis, and are
relatively scarce in comparison to those for linear spectral statistics. Even
when such formulas are available, they may be very different from those for
linear spectral statistics—as they may require the estimation of different
model parameters, or the implementation of different algorithms for numerical
evaluation. On the other hand, the bootstrap approach may be considered more
user-friendly, since it can be applied to different types of spectral
statistics in an automatic and unified way. Similarly, the bootstrap approach
can provide the user with the freedom to easily explore statistics that depend
on several linear spectral statistics in a complicated manner (which would
otherwise require intricate delta-method calculations), or statistics for
which formulas may not be available at all.

Notation and terminology. For a random object $Y$, the expression
$\mathcal{L}(Y)$ denotes its distribution, and $\mathcal{L}(Y|X)$ denotes its
conditional distribution given the observations
$\mathbf{x}_{1},\dots,\mathbf{x}_{n}$. Similarly, we use $|X$ when referring
to probabilities, expectations, and variances that are conditional on the
observations. The symbols $\xrightarrow{\mathbb{P}}$ and $\Rightarrow$
respectively denote convergence in probability and convergence in
distribution. For a set $A\subset\mathbb{R}^{k}$ and a number $\delta>0$, the
outer $\delta$-neighborhood of $A$ is defined as
$A^{\delta}=\\{x\in\mathbb{R}^{k}|\inf_{a\in A}\|x-a\|_{2}\leq\delta\\}$,
where $\|\cdot\|_{2}$ is the Euclidean norm. If $\mathbf{v}$ and $\mathbf{w}$
are random vectors in $\mathbb{R}^{k}$, then the Lévy-Prohorov metric
$d_{\textup{LP}}(\mathcal{L}(\mathbf{v}),\mathcal{L}(\mathbf{w}))$ between
their distributions is defined as the infimum over all numbers $\delta>0$ such
that the inequality $\mathbb{P}(\mathbf{v}\in A)\leq\mathbb{P}(\mathbf{w}\in
A^{\delta})+\delta$ holds for all Borel sets $A\subset\mathbb{R}^{k}$. If
$\nu,\nu_{1},\nu_{2},\dots$ is a sequence of random probability distributions
on $\mathbb{R}^{k}$, then the expression
$\nu_{n}\overset{\mathbb{P}}{\Rightarrow}\nu$ means that the sequence of
scalar random variables $d_{\textup{LP}}(\nu_{n},\nu)$ converges to 0 in
probability as $n\to\infty$. For two sequences of non-negative real numbers
$\\{a_{n}\\}$ and $\\{b_{n}\\}$, we write $a_{n}\lesssim b_{n}$ if there is a
constant $C>0$ not depending on $n$ such that $a_{n}\leq Cb_{n}$ holds for all
large $n$. When both of the relations $a_{n}\lesssim b_{n}$ and $b_{n}\lesssim
a_{n}$ hold, we write $a_{n}\asymp b_{n}$. The relation $a_{n}=o(b_{n})$ means
$a_{n}/b_{n}\to 0$ as $n\to\infty$, and the relation
$a_{n}=\mathcal{O}(b_{n})$ is equivalent to $a_{n}\lesssim b_{n}$. The
$k\times k$ identity matrix is denoted as $I_{k}$, and indicator function for
a condition $\cdots$ is denoted as $1\\{\cdots\\}$. Lastly, we use
$\mathbb{C}^{+}$ to refer to the set of complex numbers with positive
imaginary part.

##  2 Method

Conceptually, the proposed method is motivated by the fact that the standard
bootstrap (based on sampling with replacement) often performs poorly when it
is applied naively to high-dimensional data, except in special situations.
This general difficulty can be understood by noting that sampling with
replacement implicitly relies on the empirical distribution of the data as a
substitute for the true data-generating distribution. In other words, the
standard bootstrap attempts to approximate a $p$-dimensional distribution in a
fully non-parametric way, which can be challenging for even moderately large
values of $p$. For this reason, alternative bootstrap methods that sample from
parametric distributions have been advocated to improve upon the standard
bootstrap in high dimensions [e.g. 33, 32, 48], and this is the viewpoint that
we pursue here.

###  2.1 Bootstrap algorithm

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

$\mathbf{x}_{i}^{*}=\xi_{i}^{*}\tilde{\Sigma}^{1/2}_{n}\mathbf{u}_{i}^{*},$

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.

######  Algorithm 1 (Bootstrap for spectral statistics).

Input: The number of bootstrap replicates $B$, as well as $\tilde{\Sigma}_{n}$
and $\widehat{\varsigma}_{n}^{2}$.
For: $b=1,\dots,B$ do in parallel

  * 1.

Generate independent random variables $g_{1}^{*},\dots,g_{n}^{*}$ from a Gamma
distribution with mean $p$ and variance $\widehat{\varsigma}_{n}^{2}$, and
then put $\xi^{*}_{i}=\sqrt{g^{*}_{i}}$ for $i=1,\dots,n$.

  * 2.

Generate independent random vectors
$\mathbf{u}^{*}_{1},\ldots,\mathbf{u}^{*}_{n}\in\mathbb{R}^{p}$ from the
uniform distribution on the unit sphere.

  * 3.

Compute
$\mathbf{x}^{*}_{i}=\xi^{*}_{i}\tilde{\Sigma}_{n}^{1/2}\mathbf{u}^{*}_{i}$ for
$i=1,\dots,n$, and then form
$\widehat{\Sigma}_{n}^{*}=\frac{1}{n}\sum_{i=1}^{n}\mathbf{x}_{i}^{*}(\mathbf{x}_{i}^{*})^{\top}$.

  * 4.

Compute the spectral statistic
$T^{*}_{n,b}=\psi(\lambda_{1}(\widehat{\Sigma}_{n}^{*}),\ldots,\lambda_{p}(\widehat{\Sigma}_{n}^{*}))$.

end for

Return: The empirical distribution of $T^{*}_{n,1},\ldots,T^{*}_{n,B}.$

Remarks. One basic but valuable feature of the algorithm is that it can be
applied with equal ease to both linear and nonlinear spectral statistics.
Another notable property is that the algorithm allows for spectral statistics
to be uncentered. This contrasts with many other bootstrap methods that are
only designed to approximate the distributions of centered statistics, i.e.
$T_{n}-\mathbb{E}(T_{n})$, and in Section 4.3 our numerical results confirm
that the algorithm can perform well with uncentered statistics.

To comment on some more technical aspects of the algorithm, the Gamma
distribution is used in step 1 because it offers a convenient way to generate
non-negative random variables whose means and variances can be matched to any
pair of positive numbers. (If the event $\widehat{\varsigma}_{n}^{\,2}=0$
happens to occur, then the Gamma distribution in step 1 is interpreted as the
point mass at $p$, so that $\xi_{i}^{*}=\sqrt{p}$ for all $i=1,\dots,n$.)
Nevertheless, the choice of the Gamma distribution for generating
$g_{1}^{*},\dots,g_{n}^{*}$ is not required. Any other family of distributions
on $[0,\infty)$ parameterized by means $\mu$ and variances $\sigma^{2}$, say
$\\{G_{\mu,\sigma^{2}}\\}$, will be compatible with our bootstrap consistency
result in Theorem 2 if it satisfies the following two conditions. First, the
pair $(\mu,\sigma^{2})$ can be set to $(p,v)$ for any integer $p\geq 1$ and
real number $v>0$. Second, there exists some fixed $\varepsilon>0$ such that
the quantity $\int|\frac{t-p}{\sqrt{p}}|^{4+\varepsilon}dG_{p,\sigma^{2}}(t)$
remains bounded when $p$ diverges and $\sigma^{2}/p$ converges to a finite
limit.

###  2.2 Variance estimation

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 the 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,

$\varsigma_{n}^{2}\ =\
p(p+2)\frac{\operatorname{var}(\|\mathbf{x}_{1}\|_{2}^{2})-2\textup{tr}(\Sigma_{n}^{2})}{\textup{tr}(\Sigma_{n})^{2}+2\textup{tr}(\Sigma_{n}^{2})}+2p.$
(2.1)

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

$\begin{split}\alpha_{n}&=\textup{tr}(\Sigma_{n}^{2})\\\\[5.69046pt]
\beta_{n}&=\operatorname{var}(\|\mathbf{x}_{1}\|_{2}^{2})\\\\[5.69046pt]
\gamma_{n}&=\textup{tr}(\Sigma_{n})^{2}.\end{split}$

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

$\displaystyle\widehat{\alpha}_{n}$
$\displaystyle=\textup{tr}(\widehat{\Sigma}_{n}^{2})-\textstyle\frac{1}{n}\textup{tr}(\widehat{\Sigma}_{n})^{2},$
$\displaystyle\widehat{\beta}_{n}$
$\displaystyle=\frac{1}{n-1}\sum_{i=1}^{n}\Big{(}\|\mathbf{x}_{i}\|_{2}^{2}-\textstyle\frac{1}{n}\sum_{i^{\prime}=1}^{n}\|\mathbf{x}_{i^{\prime}}\|_{2}^{2}\Big{)}^{2},$
$\displaystyle\widehat{\gamma}_{n}$
$\displaystyle=\textup{tr}(\widehat{\Sigma}_{n})^{2}.$

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

$\displaystyle\widehat{\varsigma}_{n}^{\,2}=\bigg{(}p(p+2)\frac{\widehat{\beta}_{n}-2\widehat{\alpha}_{n}}{\widehat{\gamma}_{n}+2\widehat{\alpha}_{n}}+2p\bigg{)}_{\\!+}$
(2.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$.

###  2.3 Spectrum estimation

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,
36, 39, 7, 25, 24]. In order to estimate
$\lambda_{1}(\Sigma_{n}),\dots,\lambda_{p}(\Sigma_{n})$ in our current
setting, we modify the method of QuEST [25], 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 [26].

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

$\displaystyle\tilde{\lambda}_{j}=\min(\widehat{\lambda}_{\textup{Q},j},\widehat{b}_{n}),$
(2.3)

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

$\displaystyle\tilde{\Sigma}_{n}=\textup{diag}(\tilde{\lambda}_{1},\ldots,\tilde{\lambda}_{p}).$
(2.4)

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$.

##  3 Theoretical results

In this section, we present three theoretical guarantees for the proposed
method. Theorem 1 establishes appropriate notions of consistency for each of
the estimators $\widehat{\varsigma}_{n}^{\,2}$ and $\tilde{\Sigma}_{n}$.
Second, our main result in Theorem 2 shows that the bootstrap samples
generated in Algorithm 1 consistently approximate the distributions of linear
spectral statistics. Lastly, Theorem 3 demonstrates the asymptotic validity of
bootstrap-based inference procedures involving nonlinear spectral statistics.

###  3.1 Setup

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.

######  Assumption 1 (Data generating model).

As $n\to\infty$, the dimension $p$ grows so that the ratio $c_{n}=p/n$
satisfies $c_{n}\to c$ for some positive constant $c$ different from 1. For
each $i=1,\dots,n$, the observation $\mathbf{x}_{i}\in\mathbb{R}^{p}$ can be
represented as

$\mathbf{x}_{i}=\xi_{i}\Sigma_{n}^{1/2}\mathbf{u}_{i},$ (3.1)

where $\Sigma_{n}\in\mathbb{R}^{p\times p}$ is a deterministic non-zero
positive semidefinite matrix, and $(\mathbf{u}_{1},\xi_{1}),\dots,$
$(\mathbf{u}_{n},\xi_{n})$ are i.i.d. random vectors in
$\mathbb{R}^{p}\times[0,\infty)$ satisfying the following conditions: The
vector $\mathbf{u}_{1}$ is drawn from the uniform distribution on the unit
sphere of $\mathbb{R}^{p}$, and is independent of $\xi_{1}$. In addition, as
$n\to\infty$, the random variable $\xi_{1}^{2}$ satisfies
$\mathbb{E}(\xi_{1}^{2})=p$, as well as the conditions

$\displaystyle\operatorname{var}\Big{(}\textstyle\frac{\xi_{1}^{2}-p}{\sqrt{p}}\Big{)}=\tau+o(1)\
\ \ \ \ \ \text{ and}\ \ \ \ \ \
\mathbb{E}\,\Big{|}\textstyle\frac{\xi_{1}^{2}-p}{\sqrt{p}}\Big{|}^{4+\varepsilon}\,\lesssim\,1,$
(3.2)

for some fixed constants $\tau\geq 0$ and $\varepsilon>0$ that do not depend
on $n$.

Remarks. With regard to the limiting value $c$ for the ratio $p/n$, the single
case of $c=1$ is excluded so that we may employ certain facts about the QuEST
method that were established in [25]. From a practical standpoint, our
proposed method can still be used effectively when $p/n=1$, as shown in
Section 4. To address the moment conditions on $\xi_{1}^{2}$, a similar set of
conditions was used in [17] to establish a high-dimensional central limit
theorem for linear spectral statistics. However, our condition involving a
$4+\varepsilon$ moment on $(\xi_{1}^{2}-p)/\sqrt{p}$ replaces a corresponding
$2+\varepsilon$ moment condition in that work. The extra bit of integrability
is used here to show that the estimators $\widehat{\varsigma}_{n}^{2}$ and
$\tilde{\Sigma}_{n}$ have suitable asymptotic properties for ensuring
bootstrap consistency.

Examples. To illustrate that the conditions in (3.2) cover a substantial range
of situations, it is possible to provide quite a few explicit examples of
distributions for $\xi_{1}^{2}$ that are conforming:

  1. 1.

Chi-Squared distribution on $p$ degrees of freedom

  2. 2.

Poisson($p$)

  3. 3.

$(1-\tau)$Negative-Binomial$(p,1-\tau)$, for any $\tau\in(0,1)$

  4. 4.

Gamma$(p/\tau,1/\tau)$, for any $\tau>0$

  5. 5.

Beta-Prime$\left(\frac{p(1+p+\tau)}{\tau},\frac{1+p+2\tau}{\tau}\right)$, for
any $\tau>0$

  6. 6.

Log-
Normal$\left(\log(p)-\frac{1}{2}\log\big{(}1+\frac{\tau}{p}\big{)},\log\big{(}1+\frac{\tau}{p}\big{)}\right)$,
for any $\tau>0$

  7. 7.

$(p+2\beta)$ Beta$(p/2,\beta)$, for any $\beta>0$

It is also possible to give a more abstract class of examples that subsumes
some of the previous ones as special cases. In detail, the conditions in (3.2)
will hold for any $\tau\geq 0$ if $\xi_{1}^{2}=\sum_{j=1}^{p}z_{1j}^{2}$ for
some independent random variables $z_{11},\dots,z_{1p}$ satisfying

$\small\textstyle\frac{1}{p}\sum_{j=1}^{p}\mathbb{E}(z_{1j}^{2})=1,\quad\textstyle\frac{1}{p}\sum_{j=1}^{p}\operatorname{var}(z_{1j}^{2})\to\tau,\text{\quad
and \quad}\displaystyle\max_{1\leq j\leq
p}\mathbb{E}|z_{1j}|^{8+2\varepsilon}\,\lesssim\,p^{1+\frac{\varepsilon}{4}}.$
(3.3)

Further details for checking the validity of the previous examples, as well as
explicit parameterizations, are provided in Appendix E.

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

$H_{n}(t)\ =\
\frac{1}{p}\displaystyle\sum_{j=1}^{p}1\\{\lambda_{j}(\Sigma_{n})\leq t\\}.$
(3.4)

######  Assumption 2 (Spectral structure).

There is a limiting spectral distribution $H$ such that as $n\to\infty$,

$H_{n}\ \Rightarrow\ H,$ (3.5)

where the support of $H$ is a finite union of closed intervals, bounded away
from $0$ and $\infty$. Furthermore, there is a fixed compact interval in
$(0,\infty)$ containing the support of $H_{n}$ for all large $n$.

Remarks. Variations of these conditions on population eigenvalues are commonly
used throughout random matrix theory. This particular set of conditions was
used in [25] to establish theoretical guarantees for the QuEST estimation
method in the context of MP models.

###  3.2 Consistency of estimators

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

$\tilde{H}_{n}(t)=\frac{1}{p}\sum_{j=1}^{p}1\\{\tilde{\lambda}_{j}\leq t\\}.$
(3.6)

######  Theorem 1.

Under Assumptions 1 and 2, the following limits hold as $n\to\infty$,

$\displaystyle\textstyle\frac{1}{p}(\widehat{\varsigma}_{n}^{2}-\varsigma_{n}^{2})\
$ $\displaystyle\xrightarrow{\mathbb{P}}0,$ (3.7) $\displaystyle\tilde{H}_{n}$
$\displaystyle\xRightarrow{\mathbb{P}}H.$ (3.8)

Remarks. The limits (3.7) and (3.8) are proved in Appendices A.1 and A.2
respectively. Although these limits can be stated in a succinct form, quite a
few details are involved in their proofs. For instance, the analysis of
$\widehat{\varsigma}_{n}^{\,2}$ is based on extensive calculations with
polynomial functions of the quadratic forms $\|\mathbf{x}_{i}\|_{2}^{2}$ and
$\mathbf{x}_{i}^{\top}\mathbf{x}_{j}$ with $i\neq j$, as well as associated
mixed moments. The consistency of $\tilde{H}_{n}$ is also notable because it
requires showing the consistency of QuEST in elliptical models, and it seems
that the consistency of QuEST has not previously been reported outside of MP
models.

###  3.3 Consistency of bootstrap

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] 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 classical Marčenko-Pastur equation (3.9) over $\mathbb{C}^{+}$.
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

$m_{H^{\prime}\\!,c^{\prime}}=\int\frac{1}{t(1-c^{\prime}-c^{\prime}zm_{H^{\prime}\\!,c^{\prime}})-z}dH^{\prime}(t)$
(3.9)

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

$\vartheta_{n}(f)\ =\ \int f(t)dH_{n,c_{n}}(t).$ (3.10)

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 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.

######  Theorem 2.

Suppose Assumptions 1 and 2 hold, and let $\mathbf{f}=(f_{1},\ldots,f_{k})$ be
a fixed set of real-valued functions that are analytic on an open set
containing $[0,\infty)$. Under these conditions, if $T_{n,1}^{*}(\mathbf{f})$
is generated as in Algorithm 1 using the estimators
$\widehat{\varsigma}_{n}^{\,2}$ and $\tilde{\Sigma}_{n}$ defined by (2.2) and
(2.4), then the following limit holds as $n\to\infty$

$d_{\mathrm{LP}}\Big{(}\mathcal{L}\big{(}p\\{T_{n}(\mathbf{f})-\vartheta_{n}(\mathbf{f})\\}\big{)}\,,\,\mathcal{L}\big{(}p\\{T_{n,1}^{*}(\mathbf{f})-\tilde{\vartheta}_{n}(\mathbf{f})\\}\big{|}X\big{)}\Big{)}\
\xrightarrow{\mathbb{P}}\ 0.$ (3.11)

Remarks. The proof is given in Appendix B, and makes key use of a recently
developed central limit theorem for linear spectral statistics due to [17].
Regarding other aspects of the theorem, there are two points to discuss.
First, the assumption that the functions $f_{1},\dots,f_{k}$ are defined on an
open set containing $[0,\infty)$ has been made for technical simplicity. In
the setting where $p/n\to c>1$, this assumption is minor because
$f_{1},\dots,f_{k}$ must be defined at 0 due to the singularity of
$\widehat{\Sigma}_{n}$. Nevertheless, if $p/n\to c\in(0,1)$, then it is
possible to show that the theorem holds for analytic functions that are not
defined at 0, and our numerical results confirm that the proposed method can
successfully handle such cases. Second, the quantities
$\vartheta_{n}(\mathbf{f})$ and $\tilde{\vartheta}_{n}(\mathbf{f})$ are only
introduced so that the distributions appearing in (3.11) have non-trivial weak
limits, which facilitates the proof. Still, the proposed method can be applied
without requiring any particular type of centering.

###  3.4 Application to inference on stable rank with guarantees

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

$r_{n}=\frac{\ \textup{tr}(\Sigma_{n})^{2}}{\textup{tr}(\Sigma_{n}^{2})},$
(3.12)

which arises naturally in a plethora of situations [e.g. 4, 41, 42, 29, 34,
31]. Whenever $\Sigma_{n}$ is non-zero, this parameter always 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.

####  3.4.1 Confidence interval for stable rank

Our confidence interval for $r_{n}$ is constructed using the estimator

$\displaystyle\widehat{r}_{n}=\frac{\textup{tr}(\widehat{\Sigma}_{n})^{2}}{\textup{tr}(\widehat{\Sigma}_{n}^{2})-\frac{1}{n}\textup{tr}(\widehat{\Sigma}_{n})^{2}},$
(3.13)

which is taken to be equal to $n$ in the exceptional case that the denominator
is 0. (The denominator is always non-negative, and it can only vanish in the
asymptotically negligible situations where $\widehat{\Sigma}_{n}=0$, or where
$\widehat{\Sigma}_{n}$ has $n$ non-zero eigenvalues that are all equal.) It
should be noted too that $\widehat{r}_{n}$ differs from the naive plug-in rule
$\textup{tr}(\widehat{\Sigma}_{n})^{2}/\textup{tr}(\widehat{\Sigma}_{n}^{2})$,
since the extra term $\frac{1}{n}\textup{tr}(\widehat{\Sigma}_{n})^{2}$ in the
denominator serves to correct the bias of
$\textup{tr}(\widehat{\Sigma}_{n}^{2})$ in estimating
$\textup{tr}(\Sigma_{n}^{2})$.

To proceed, let $q_{1-\alpha}$ denote $(1-\alpha)$-quantile of the random
variable $(\widehat{r}_{n}-r_{n})/p$ for any fixed $\alpha\in(0,1)$, and
consider the interval
$\mathcal{I}_{n}=[\widehat{r}_{n}-pq_{1-\alpha/2},\widehat{r}_{n}-pq_{\alpha/2}]$.
Whenever the distribution of $(\widehat{r}_{n}-r_{n})/p$ is continuous, this
interval satisfies

$\mathbb{P}(r_{n}\in\mathcal{I}_{n})=1-\alpha.$ (3.14)

However, the quantiles $q_{\alpha/2}$ and $q_{1-\alpha/2}$ are unknown, and so
they must be estimated. This can be accomplished by generating bootstrap
samples of the following form in Algorithm 1,

$T^{*}_{n,1}=\frac{\textstyle\frac{1}{p}\textup{tr}(\widehat{\Sigma}^{*}_{n})^{2}}{\textup{tr}((\widehat{\Sigma}^{*}_{n})^{2})-\frac{1}{n}\textup{tr}(\widehat{\Sigma}_{n}^{*})^{2}}\,-\,\frac{\textstyle\frac{1}{p}\textup{tr}(\tilde{\Sigma}_{n})^{2}}{\textup{tr}(\tilde{\Sigma}_{n}^{2})},$

where $T_{n,1}^{*}$ is arbitrarily set to 0 if either denominator is 0.
Letting $\widehat{q}_{\alpha/2}$ and $\widehat{q}_{1-\alpha/2}$ denote the
respective $\alpha/2$ and $(1-\alpha/2)$-quantiles of
$\mathcal{L}(T_{n,1}^{*}|X)$, the proposed confidence interval is defined as

$\widehat{\mathcal{I}}_{n}\ =\
\Big{[}\widehat{r}_{n}-p\widehat{q}_{1-\alpha/2}\,,\,\widehat{r}_{n}-p\widehat{q}_{\alpha/2}\Big{]}.$
(3.15)

Below, Theorem 3 shows that as $n\to\infty$, the coverage probability
$\mathbb{P}(r_{n}\in\widehat{\mathcal{I}}_{n})$ converges to $1-\alpha$, as
desired.

####  3.4.2 Hypotheses related to stable rank

Screening data for PCA. Consider a scenario where a collection of different
datasets are to be screened for further investigation by principal components
analysis (PCA). In this situation, the datasets that should be discarded are
the ones that cannot be well summarized by a moderate number of principal
components. To put this in more quantitative terms, a dataset may be
considered unsuitable for PCA if the stable rank $r_{n}$ exceeds a certain
fraction of the full dimension $p$. That is, if $\frac{r_{n}}{p}>\epsilon_{0}$
for some fixed reference value $\epsilon_{0}\in(0,1)$. On the other hand, if
$\frac{r_{n}}{p}\leq\epsilon_{0}$, then the dataset may be retained. This
leads to considering the hypothesis testing problem

$\mathsf{H}_{0,n}:\textstyle\frac{r_{n}}{p}\leq\epsilon_{0}\text{ \ \ \ \ \ \
vs. \ \ \ \ \ \ }\mathsf{H}_{1,n}:\textstyle\frac{r_{n}}{p}>\epsilon_{0}.$
(3.16)

To develop a testing procedure, we may again consider the
$(1-\alpha)$-quantile $q_{1-\alpha}$ of the random variable
$(\widehat{r}_{n}-r_{n})/p$. This quantile serves as a conceptual basis for a
rejection criterion, because it satisfies the following inequality under the
null hypothesis

$\mathbb{P}\Big{(}\textstyle\frac{\widehat{r}_{n}}{p}-\epsilon_{0}>q_{1-\alpha}\Big{)}\
\leq\ \alpha.$ (3.17)

In other words, if $q_{1-\alpha}$ were known, then a level-$\alpha$ testing
procedure would result from using
$\frac{\widehat{r}_{n}}{p}-\epsilon_{0}>q_{1-\alpha}$ as a rejection
criterion. Accordingly, a bootstrap-based version of this procedure rejects
the null hypothesis when
$\frac{\widehat{r}_{n}}{p}-\epsilon_{0}>\widehat{q}_{1-\alpha}$, with the
quantile estimate $\widehat{q}_{1-\alpha}$ defined as before.

Testing for sphericity. One more example of a testing problem related to
$r_{n}$ is that of testing for sphericity,

$\mathsf{H}^{\prime}_{0,n}:\Sigma_{n}=\sigma^{2}I\text{ \ for some \
$\sigma^{2}>0$}\text{ \ \ \ vs. \ \ \
}\mathsf{H}^{\prime}_{1,n}:\Sigma_{n}\neq\sigma^{2}I\text{ \ for any \
$\sigma^{2}>0$}.$ (3.18)

The connection to the parameter $r_{n}$ arises from the fact that (3.18) is
equivalent to the problem $\mathsf{H}^{\prime}_{0,n}:\frac{r_{n}}{p}=1$ vs.
$\mathsf{H}^{\prime}_{1,n}:\frac{r_{n}}{p}<1$. This observation was used in
the influential paper [41] to develop a formula-based sphericity test. By
analogy with our discussion of the problem (3.16), this observation can also
be used to develop a bootstrap-based sphericity test. However, there is one
point of distinction in the current situation, which is that the eigenvalues
of $\Sigma_{n}$ no longer need to be estimated. The reason is that under
$\mathsf{H}_{0,n}^{\prime}$, the scale-invariance of the statistic
$(\widehat{r}_{n}-r_{n})/p$ causes it to behave as if $\Sigma_{n}=I_{p}$.
Consequently, to estimate the quantiles of the null distribution of
$(\widehat{r}_{n}-r_{n})/p$, Algorithm 1 can be run using
$\tilde{\Sigma}_{n}=I_{p}$. To summarize, if we let
$\widehat{q}_{\alpha}^{\prime}$ denote the resulting estimate for the
$\alpha$-quantile of the null distribution of $(\widehat{r}_{n}-r_{n})/p$,
then the rejection criterion is
$\frac{\widehat{r}_{n}}{p}-1\leq\widehat{q}_{\alpha}^{\prime}$.

The following result establishes the theoretical validity of the procedures
discussed in this subsection.

######  Theorem 3.

If Assumptions 1 and 2 hold, then as $n\to\infty$,

$\mathbb{P}\big{(}r_{n}\in\widehat{\mathcal{I}}_{n}\big{)}\ \to\ 1-\alpha.$
(3.19)

Furthermore, if $\mathsf{H}_{0,n}$ holds for all large $n$, then

$\mathbb{P}\Big{(}\textstyle\frac{\widehat{r}_{n}}{p}-\epsilon_{0}>\widehat{q}_{1-\alpha}\Big{)}\
\leq\ \alpha+o(1),$ (3.20)

or if $\mathsf{H}^{\prime}_{0,n}$ holds for all large $n$, then

$\mathbb{P}\Big{(}\textstyle\frac{\widehat{r}_{n}}{p}-1\leq\widehat{q}_{\alpha}^{\prime}\Big{)}\
\to\ \alpha.$ (3.21)

Remarks. This result provides a notable complement to Theorem 2, because it
demonstrates that bootstrap consistency can be established in tasks that are
based on a nonlinear spectral statistic, namely $\widehat{r}_{n}-r_{n}$. The
proof of this result is given in Appendix C, where it can be seen that the
limiting distribution of $\widehat{r}_{n}-r_{n}$ has a very complicated
dependence on the moments of $H$, due to the correlations between
$\textup{tr}(\widehat{\Sigma}_{n})$ and
$\textup{tr}(\widehat{\Sigma}_{n}^{2})$. In this way, the proof illustrates
the utility of the bootstrap, since the bootstrap enables the user to
completely bypass such complexity.

##  4 Numerical results

This section explores the empirical performance of the proposed bootstrap
method in three different ways. Sections 4.2 and 4.3 deal with bootstrap
approximations for linear and nonlinear spectral statistics, while Section 4.4
looks at procedures for doing inference on the stable rank parameter $r_{n}$.

###  4.1 Parameter settings

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

  1. (i).

Chi-Squared distribution on $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)$.

The population covariance matrix $\Sigma_{n}$ was selected from the three
options below:

  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/10)$ 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.2 Linear spectral statistics

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\\}$.

Design of experiments. For each parameter setting and choice of $p/n$, we
generated 5000 realizations of the dataset
$\mathbf{x}_{1},\dots,\mathbf{x}_{n}$, with $n=400$. These datasets allowed us
to compute 5000 realizations of the statistic $p(T_{n}(f)-\vartheta_{n}(f))$,
and we treated the empirical mean, standard deviation, and 95th percentile of
these 5000 realizations as ground truth for our parameters of interest. In
Tables 1 and 2, the ground truth values are reported in the first row of
numbers corresponding to each choice of $\Sigma_{n}$.

With regard to the bootstrap, we ran Algorithm 1 on the first 500 datasets
corresponding to each parameter setting and choice of $p/n$. 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 and choice of
$p/n$, we approximated $\vartheta_{n}(f)$ by averaging 30 realizations of
$\frac{1}{40p}\sum_{j=1}^{40p}f(\lambda_{j}(\frac{1}{40n}\sum_{i=1}^{40n}\xi_{i}^{2}\mathbf{u}_{i}\mathbf{u}_{i}^{\top}(I_{40}\otimes\Sigma_{n})))$,
where 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), (ii), and (iii), but with $40p$ replacing $p$. Also, we used one
realization of
$\frac{1}{40p}\sum_{j=1}^{40p}f(\lambda_{j}(\frac{1}{40n}\sum_{i=1}^{40n}(\xi_{i}^{*})^{2}\mathbf{u}_{i}^{*}(\mathbf{u}_{i}^{*})^{\top}(I_{40}\otimes\tilde{\Sigma}_{n})))$
to approximate $\tilde{\vartheta}_{n}(f)$ during every run of Algorithm 1,
where 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}$.

Comments on results. To explain the reading of the tables with an example, the
entries in the upper right corner of Table 1 show that in settings (i) and (1)
with $p/n=1.5$, the 95th percentile of the distribution
$\mathcal{L}(p(T_{n}(f)-\vartheta_{n}(f)))$ corresponding to $f(x)=x^{2}$ is
equal to 16.48, and the bootstrap estimate for the 95th percentile has a
mean(standard deviation) of 16.73(1.48). Overall, Tables 1 and 2 reveal that
the bootstrap estimates have low bias, and have standard deviations that tend
to be a small proportion of the target parameter.

Table 1: Results for $p(T_{n}(f)-\vartheta_{n}(f))$ with $f(x)=x^{2}$

|  | $p/n=0.5$ |  | $p/n=1$ |  | $p/n=1.5$
---|---|---|---|---|---|---
$\xi_{1}$ | $\Sigma_{n}$ | mean | sd | 95th |  | mean | sd | 95th |  | mean | sd | 95th
| 1 | 0.51 | 3.31 | 6.02 |  | 0.85 | 6.11 | 10.87 |  | 1.44 | 9.15 | 16.48
| 0.48(0.24) | 3.27(0.19) | 5.85(0.53) |  | 0.96(0.44) | 6.09(0.35) | 11.02(1) |  | 1.47(0.65) | 9.25(0.51) | 16.73(1.48)
(i) | 2 | 0 | 0.11 | 0.19 |  | 0.01 | 0.11 | 0.19 |  | 0 | 0.11 | 0.2
0(0.01) | 0.11(0.01) | 0.19(0.03) |  | 0(0.01) | 0.11(0.01) | 0.19(0.02) |  | 0(0.01) | 0.11(0.01) | 0.2(0.03)
| 3 | 2.16 | 12.88 | 23.12 |  | 4.02 | 24.12 | 44.03 |  | 5.8 | 36.13 | 65.72
| 1.9(0.88) | 12.85(0.73) | 23.2(2.07) |  | 3.94(1.75) | 24.27(1.35) | 43.95(3.82) |  | 5.8(2.77) | 37.08(2.2) | 67.09(6.1)
| 1 | 3.57 | 6.35 | 14.09 |  | 6.69 | 11.93 | 26.9 |  | 10.25 | 18.09 | 39.92
| 3.47(0.59) | 6.36(0.45) | 14.03(1.29) |  | 6.83(1) | 11.85(0.76) | 26.51(2.32) |  | 10.26(1.53) | 17.99(1.11) | 39.92(3.25)
(ii) | 2 | 0.01 | 0.12 | 0.21 |  | 0 | 0.11 | 0.2 |  | 0 | 0.11 | 0.19
0.01(0.01) | 0.12(0.01) | 0.2(0.03) |  | 0.01(0.01) | 0.11(0.01) | 0.2(0.03) |  | 0(0.01) | 0.11(0.01) | 0.2(0.03)
| 3 | 14.16 | 25.27 | 56 |  | 27.59 | 47.57 | 105.75 |  | 41.49 | 71.66 | 159.75
| 13.67(2.24) | 25.12(1.76) | 55.35(5.1) |  | 26.94(3.72) | 47.09(2.97) | 104.85(8.33) |  | 40.39(6.13) | 71.59(4.5) | 158.55(13.38)
| 1 | -0.46 | 1.13 | 1.44 |  | -0.95 | 2.13 | 2.61 |  | -1.47 | 3.12 | 3.71
| -0.47(0.08) | 1.14(0.07) | 1.4(0.17) |  | -0.96(0.14) | 2.15(0.11) | 2.58(0.33) |  | -1.43(0.21) | 3.17(0.15) | 3.81(0.45)
(iii) | 2 | 0.01 | 0.11 | 0.19 |  | 0 | 0.11 | 0.19 |  | 0 | 0.11 | 0.19
0(0.01) | 0.11(0.01) | 0.19(0.03) |  | 0(0.01) | 0.11(0.01) | 0.19(0.02) |  | 0(0.01) | 0.11(0.01) | 0.19(0.03)
| 3 | -1.93 | 4.49 | 5.41 |  | -3.89 | 8.66 | 10.62 |  | -5.62 | 13 | 15.33
| -1.85(0.33) | 4.56(0.25) | 5.68(0.69) |  | -3.78(0.57) | 8.69(0.48) | 10.57(1.26) |  | -5.7(0.9) | 12.87(0.68) | 15.52(1.88)

Table 2: Results for $p(T_{n}(f)-\vartheta_{n}(f))$ with $f(x)=x-\log(x)-1$

|  | $p/n=0.3$ |  | $p/n=0.5$ |  | $p/n=0.7$
---|---|---|---|---|---|---
$\xi_{1}$ | $\Sigma_{n}$ | mean | sd | 95th |  | mean | sd | 95th |  | mean | sd | 95th
| 1 | 0.18 | 0.35 | 0.77 |  | 0.34 | 0.64 | 1.4 |  | 0.6 | 1 | 2.24
| 0.17(0.03) | 0.35(0.02) | 0.74(0.05) |  | 0.34(0.06) | 0.63(0.03) | 1.37(0.09) |  | 0.58(0.09) | 1.01(0.05) | 2.25(0.14)
(i) | 2 | 0.18 | 0.83 | 1.55 |  | 0.36 | 1.16 | 2.25 |  | 0.59 | 1.53 | 3.1
0.19(0.22) | 0.79(0.26) | 1.49(0.65) |  | 0.36(0.45) | 1.11(0.37) | 2.19(1.05) |  | 0.62(0.7) | 1.51(0.44) | 3.1(1.42)
| 3 | 0.18 | 0.86 | 1.62 |  | 0.35 | 1.19 | 2.27 |  | 0.55 | 1.58 | 3.12
| 0.17(0.06) | 0.86(0.05) | 1.59(0.14) |  | 0.33(0.09) | 1.19(0.06) | 2.3(0.19) |  | 0.58(0.12) | 1.56(0.08) | 3.17(0.25)
| 1 | 1 | 0.37 | 1.62 |  | 1.74 | 0.65 | 2.79 |  | 2.56 | 1.02 | 4.25
| 1.04(0.11) | 0.37(0.02) | 1.66(0.12) |  | 1.79(0.17) | 0.66(0.03) | 2.87(0.19) |  | 2.64(0.22) | 1.04(0.05) | 4.35(0.25)
(ii) | 2 | 1.04 | 1.52 | 3.54 |  | 1.74 | 2.02 | 5.06 |  | 2.54 | 2.49 | 6.54
1.05(0.4) | 1.51(0.25) | 3.54(0.81) |  | 1.77(0.7) | 2.01(0.35) | 5.09(1.27) |  | 2.54(0.96) | 2.46(0.41) | 6.59(1.63)
| 3 | 1.02 | 1.63 | 3.78 |  | 1.78 | 2.15 | 5.34 |  | 2.57 | 2.63 | 6.94
| 1.04(0.16) | 1.6(0.12) | 3.68(0.36) |  | 1.8(0.22) | 2.11(0.14) | 5.31(0.45) |  | 2.61(0.28) | 2.59(0.16) | 6.88(0.53)
| 1 | -0.11 | 0.34 | 0.45 |  | -0.13 | 0.62 | 0.9 |  | -0.14 | 1 | 1.49
| -0.11(0.02) | 0.34(0.02) | 0.45(0.05) |  | -0.14(0.04) | 0.62(0.03) | 0.88(0.08) |  | -0.09(0.07) | 1.01(0.04) | 1.57(0.14)
(iii) | 2 | -0.12 | 0.36 | 0.49 |  | -0.14 | 0.64 | 0.93 |  | -0.08 | 1.01 | 1.57
-0.04(0.12) | 0.47(0.18) | 0.74(0.43) |  | 0.05(0.29) | 0.85(0.28) | 1.45(0.74) |  | 0.18(0.45) | 1.22(0.31) | 2.19(0.95)
| 3 | -0.11 | 0.38 | 0.52 |  | -0.14 | 0.65 | 0.93 |  | -0.07 | 1.03 | 1.66
| -0.11(0.03) | 0.38(0.02) | 0.51(0.05) |  | -0.14(0.05) | 0.65(0.03) | 0.93(0.09) |  | -0.09(0.07) | 1.03(0.05) | 1.6(0.14)

Table 3: Results for $\lambda_{1}(\widehat{\Sigma}_{n})$

|  | $p/n=0.5$ |  | $p/n=1$ |  | $p/n=1.5$
---|---|---|---|---|---|---
$\xi_{1}$ | $\Sigma_{n}$ | mean | sd | 95th |  | mean | sd | 95th |  | mean | sd | 95th
| 1 | 2.9 | 0.06 | 3 |  | 3.96 | 0.06 | 4.07 |  | 4.9 | 0.06 | 5.01
| 2.93(0.05) | 0.06(0.01) | 3.04(0.07) |  | 3.99(0.04) | 0.06(0.01) | 4.1(0.06) |  | 4.93(0.04) | 0.07(0.01) | 5.05(0.06)
(i) | 2 | 0.8 | 0.05 | 0.89 |  | 0.8 | 0.05 | 0.89 |  | 0.81 | 0.05 | 0.9
0.8(0.05) | 0.05(0.01) | 0.89(0.06) |  | 0.8(0.05) | 0.05(0.01) | 0.89(0.06) |  | 0.81(0.05) | 0.05(0.01) | 0.9(0.06)
| 3 | 5.76 | 0.11 | 5.95 |  | 7.91 | 0.12 | 8.12 |  | 9.81 | 0.13 | 10.03
| 5.81(0.08) | 0.12(0.02) | 6.02(0.11) |  | 7.98(0.09) | 0.13(0.02) | 8.21(0.13) |  | 9.87(0.1) | 0.14(0.03) | 10.12(0.14)
| 1 | 2.96 | 0.07 | 3.08 |  | 4.02 | 0.07 | 4.14 |  | 4.96 | 0.07 | 5.08
| 3.03(0.07) | 0.07(0.01) | 3.15(0.09) |  | 4.09(0.07) | 0.07(0.01) | 4.22(0.09) |  | 5.04(0.06) | 0.08(0.01) | 5.17(0.08)
(ii) | 2 | 0.8 | 0.06 | 0.89 |  | 0.8 | 0.05 | 0.89 |  | 0.8 | 0.05 | 0.9
0.8(0.05) | 0.05(0.01) | 0.89(0.06) |  | 0.8(0.06) | 0.05(0.01) | 0.89(0.06) |  | 0.81(0.06) | 0.05(0.01) | 0.9(0.06)
| 3 | 5.88 | 0.13 | 6.11 |  | 8.04 | 0.14 | 8.26 |  | 9.93 | 0.14 | 10.16
| 6.03(0.13) | 0.14(0.02) | 6.27(0.16) |  | 8.2(0.12) | 0.15(0.02) | 8.47(0.15) |  | 10.1(0.12) | 0.16(0.02) | 10.37(0.16)
| 1 | 2.88 | 0.05 | 2.97 |  | 3.94 | 0.06 | 4.04 |  | 4.88 | 0.06 | 4.99
| 2.89(0.03) | 0.06(0.01) | 2.98(0.05) |  | 3.95(0.03) | 0.06(0.01) | 4.06(0.05) |  | 4.89(0.03) | 0.06(0.01) | 5(0.05)
(iii) | 2 | 0.8 | 0.05 | 0.89 |  | 0.8 | 0.05 | 0.89 |  | 0.8 | 0.05 | 0.9
0.8(0.05) | 0.05(0.01) | 0.89(0.06) |  | 0.8(0.05) | 0.05(0.01) | 0.89(0.06) |  | 0.8(0.05) | 0.05(0.01) | 0.89(0.06)
| 3 | 5.72 | 0.1 | 5.9 |  | 7.88 | 0.12 | 8.08 |  | 9.77 | 0.12 | 9.99
| 5.74(0.07) | 0.11(0.02) | 5.94(0.11) |  | 7.9(0.08) | 0.12(0.02) | 8.12(0.13) |  | 9.79(0.08) | 0.13(0.02) | 10.01(0.11)

Table 4: Results for
$\lambda_{1}(\widehat{\Sigma}_{n})-\lambda_{2}(\widehat{\Sigma}_{n})$

|  | $p/n=0.5$ |  | $p/n=1$ |  | $p/n=1.5$
---|---|---|---|---|---|---
$\xi_{1}$ | $\Sigma_{n}$ | mean | sd | 95th |  | mean | sd | 95th |  | mean | sd | 95th
| 1 | 0.09 | 0.05 | 0.18 |  | 0.1 | 0.05 | 0.19 |  | 0.1 | 0.06 | 0.21
| 0.1(0.02) | 0.06(0.01) | 0.2(0.05) |  | 0.1(0.02) | 0.06(0.01) | 0.21(0.04) |  | 0.11(0.02) | 0.06(0.01) | 0.22(0.04)
(i) | 2 | 0.18 | 0.07 | 0.3 |  | 0.18 | 0.07 | 0.29 |  | 0.18 | 0.07 | 0.29
0.19(0.06) | 0.06(0.01) | 0.29(0.07) |  | 0.18(0.06) | 0.06(0.01) | 0.29(0.07) |  | 0.19(0.06) | 0.06(0.01) | 0.3(0.07)
| 3 | 0.18 | 0.1 | 0.35 |  | 0.19 | 0.1 | 0.38 |  | 0.2 | 0.11 | 0.41
| 0.2(0.04) | 0.11(0.02) | 0.4(0.08) |  | 0.22(0.05) | 0.12(0.03) | 0.44(0.1) |  | 0.23(0.05) | 0.13(0.03) | 0.47(0.11)
| 1 | 0.1 | 0.05 | 0.2 |  | 0.1 | 0.05 | 0.2 |  | 0.11 | 0.06 | 0.21
| 0.11(0.02) | 0.06(0.01) | 0.22(0.05) |  | 0.11(0.02) | 0.06(0.01) | 0.23(0.05) |  | 0.12(0.02) | 0.07(0.01) | 0.24(0.04)
(ii) | 2 | 0.18 | 0.07 | 0.29 |  | 0.18 | 0.07 | 0.29 |  | 0.18 | 0.07 | 0.29
0.18(0.06) | 0.06(0.01) | 0.29(0.07) |  | 0.19(0.06) | 0.06(0.01) | 0.29(0.07) |  | 0.19(0.06) | 0.06(0.01) | 0.3(0.07)
| 3 | 0.19 | 0.1 | 0.38 |  | 0.2 | 0.11 | 0.41 |  | 0.21 | 0.11 | 0.42
| 0.22(0.04) | 0.12(0.02) | 0.44(0.09) |  | 0.23(0.05) | 0.13(0.03) | 0.48(0.1) |  | 0.24(0.04) | 0.13(0.03) | 0.49(0.1)
| 1 | 0.09 | 0.05 | 0.17 |  | 0.09 | 0.05 | 0.19 |  | 0.1 | 0.06 | 0.21
| 0.09(0.02) | 0.05(0.01) | 0.19(0.04) |  | 0.1(0.02) | 0.05(0.01) | 0.2(0.03) |  | 0.11(0.01) | 0.06(0.01) | 0.21(0.03)
(iii) | 2 | 0.18 | 0.07 | 0.29 |  | 0.18 | 0.07 | 0.29 |  | 0.18 | 0.07 | 0.29
0.19(0.06) | 0.06(0.01) | 0.3(0.07) |  | 0.18(0.06) | 0.06(0.01) | 0.29(0.07) |  | 0.18(0.06) | 0.06(0.01) | 0.29(0.07)
| 3 | 0.17 | 0.09 | 0.34 |  | 0.19 | 0.1 | 0.38 |  | 0.2 | 0.11 | 0.41
| 0.19(0.04) | 0.1(0.02) | 0.38(0.08) |  | 0.21(0.05) | 0.11(0.03) | 0.42(0.1) |  | 0.18(0.06) | 0.06(0.01) | 0.29(0.07)

###  4.3 Nonlinear spectral statistics

This subsection presents numerical results for nonlinear spectral statistics
in a parallel manner to the previous subsection. The particular statistics
under consideration here are the largest sample eigenvalue
$\lambda_{1}(\widehat{\Sigma}_{n})$, and the leading eigengap
$\lambda_{1}(\widehat{\Sigma}_{n})-\lambda_{2}(\widehat{\Sigma}_{n})$. The
underlying experiments for these statistics were designed in the same manner
as before, and Tables 3 and 4 display the results in the same format. To a
large extent, the favorable patterns that were noted in the results for linear
spectral statistics are actually enhanced in the results for nonlinear
spectral statistics—in the sense that the bias and standard deviations of the
bootstrap estimates are generally smaller here than before.

###  4.4 Inference on stable rank

Confidence interval. Table 5 presents numerical results on the width and
coverage probability of the bootstrap confidence interval
$\widehat{\mathcal{I}}$ for the stable rank parameter $r_{n}$, as defined in
(3.15). The results were computed using experiments based on the same design
as in Section 4.2, with every interval having a nominal coverage probability
of $95\%$. Due to the fact that the width of the interval scales with the
value of $r_{n}$, we report the width as a percentage of $r_{n}$. To
illustrate a particular example, the upper right corner of Table 5 shows that
in settings (1) and (i) with $p/n=1.5$, the average width of the interval over
repeated experiments is 1.96% of $r_{n}$, with a standard deviation of 0.12%.
When $\Sigma_{n}$ is chosen as in cases 1 and 3, the width is typically quite
small as a percentage of $r_{n}$. In case 2, the percentage is larger and
ranges between 14% and 19%, which seems to occur because $r_{n}$ is much
smaller in case 2 compared to cases 1 and 3. Lastly, to consider coverage
probability, the table shows very good general agreement with the nominal
level, and only a handful of cases differ from the nominal level by more than
1%.

Table 5: Results for the stable rank confidence interval

|  | $p/n=0.5$ |  | $p/n=1$ |  | $p/n=1.5$
---|---|---|---|---|---|---
$\xi_{1}$ | $\Sigma_{n}$ |  width$/r_{n}$ $(\%)$ |  coverage $(\%)$ |  |  width$/r_{n}$ $(\%)$ |  coverage $(\%)$ |  |  width$/r_{n}$ $(\%)$ |  coverage $(\%)$
(i) | 1 | 1.98(0.12) | 94.40 |  | 1.96(0.08) | 94.30 |  | 1.96(0.12) | 96.20
2 | 14.44(1.07) | 93.60 |  | 17.04(1.29) | 95.20 |  | 18.68(1.41) | 92.40
3 | 1.98(0.12) | 92.20 |  | 1.98(0.11) | 95.20 |  | 1.97(0.12) | 95.60
(ii) | 1 | 2.07(0.13) | 95.80 |  | 2.00(0.12) | 93.20 |  | 1.98(0.12) | 95.40
2 | 14.60(1.10) | 95.00 |  | 17.08(1.24) | 95.60 |  | 18.60(1.36) | 93.60
3 | 2.06(0.13) | 94.80 |  | 2.01(0.12) | 94.60 |  | 1.98(0.12) | 95.80
(iii) | 1 | 1.97(0.13) | 95.00 |  | 1.96(0.08) | 94.70 |  | 1.96(0.12) | 95.40
2 | 14.43(1.09) | 95.20 |  | 16.98(1.25) | 94.80 |  | 18.61(1.48) | 95.40
3 | 1.98(0.12) | 95.20 |  | 1.97(0.12) | 94.00 |  | 1.97(0.12) | 93.40

Stable rank test. To conclude this section, we discuss numerical results for
an instance of the testing problem (3.16) given by

$\mathsf{H}_{0,n}:\textstyle\frac{r_{n}}{p}\leq 0.1\ \ \ \text{\ \ \ \ \ \ \
vs. \ \ \ \ \ \ \ }\ \ \ \mathsf{H}_{1,n}:\textstyle\frac{r_{n}}{p}>0.1.$
(4.1)

As a way to unify the study of level and power, we modified the experiments
from Section 4.2 as follows. We rescaled the leading 15 eigenvalues in setting
(1) to tune the ratio $r_{n}/p$ within the grid
$\\{0.098,0.985,\dots,0.105\\}$. More precisely, the eigenvalues of
$\Sigma_{n}$ were taken to be of the form $\lambda_{j}(\Sigma_{n})=4s/3$ for
$j\in\\{1,\dots,5\\}$, $\lambda_{j}(\Sigma_{n})=s$ for $j\in\\{6,\dots,15\\}$
and $\lambda_{j}(\Sigma_{n})=1$ for $j\in\\{16,\dots,p\\}$, with different
values of $s$ being chosen to produce values of $r_{n}/p$ matching the stated
gridpoints. Hence, gridpoints to the left of 0.1 correspond to
$\mathsf{H}_{n,0}$, and gridpoints to the right of $0.1$ correspond to
$\mathsf{H}_{1,n}$.

At each gridpoint, we performed experiments following the design of those in
Section 4.2, allowing for $p/n$ to take the values $0.5,1.0,1.5$, and allowing
the distribution of $\xi_{1}^{2}$ to be of the types (i), (ii), and (iii). For
each such parameter setting, we applied the relevant bootstrap test from
Section 3.4.2 at a 5% nominal level to 500 datasets, and then we recorded the
rejection rate over the 500 trials. Figures 1, 2, and 3 display the results by
plotting the rejection rate as a function of the ratio $r_{n}/p$. The separate
figures correspond to choices of the distribution $\mathcal{L}(\xi_{1})$, and
within each figure, three colored curves correspond to choices of $p/n$, as
indicated in the legend. In addition, all the plots include a dashed
horizontal line to signify the 5% nominal level.

An important feature that is shared by all the curves is that they stay below
the 5% level for essentially every value of $r_{n}/p<0.1$, which corroborates
our theoretical bound (3.20) in Theorem 3. Furthermore, when $r_{n}/p=0.1$,
the curves are mostly quite close to 5%, demonstrating that the testing
procedure is well calibrated. For values of $r_{n}/p>0.1$, the procedure
exhibits substantial power, with the curve corresponding to $p/n=0.5$
achieving approximately 100% power at $r_{n}/p=0.105$ for every choice of
$\mathcal{L}(\xi_{1})$. In the cases of $p/n\in\\{1.0,1.5\\}$, the procedure
still retains power, but with some diminution, as might be anticipated in
settings of higher dimension.

\begin{overpic}[width=281.85034pt]{stable_rank_distr4.eps} \put(32.0,75.0){
rejection rate vs. $r_{n}/p$ } \put(0.0,17.0){\rotatebox{90.0}{ {\small\ \ \
rejection rate (\%) \ \ }}} \put(45.0,-1.0){ $r_{n}/p\ (\%)$ } \end{overpic}
Figure 1:  Results for the hypothesis testing problem (4.1) when $\xi_{1}^{2}$
follows a Chi-Squared distribution on $p$ degrees of freedom.
\begin{overpic}[width=281.85034pt]{stable_rank_distr2.eps} \put(32.0,75.0){
rejection rate vs. $r_{n}/p$ } \put(0.0,17.0){\rotatebox{90.0}{ {\small\ \ \
rejection rate (\%) \ \ }}} \put(45.0,-1.0){ $r_{n}/p\ (\%)$ } \end{overpic}
Figure 2: Results for the hypothesis testing problem (4.1) when $\xi_{1}^{2}$
follows a Beta-Prime$\left(\frac{p(1+p+8)}{8},\frac{1+p+16}{8}\right)$
distribution. \begin{overpic}[width=281.85034pt]{stable_rank_distr3.eps}
\put(32.0,75.0){ rejection rate vs. $r_{n}/p$ }
\put(0.0,17.0){\rotatebox{90.0}{ {\small\ \ \ rejection rate (\%) \ \ }}}
\put(45.0,-1.0){ $r_{n}/p\ (\%)$ } \end{overpic} Figure 3: Results for the
hypothesis testing problem (4.1) when $\xi_{1}^{2}$ follows a
$(p+4)\textup{Beta}(p/2,2)$ distribution.

## References

  * Akemann et al. [2011] Akemann, G., Baik, J., and Di Francesco, P.  _The Oxford Handbook of Random Matrix Theory_.  Oxford, 2011.
  * Anderson [2003] Anderson, T. W.  _Introduction to Multivariate Statistical Analysis_.  Wiley, 2003.
  * Bai and Silverstein [2004] Bai, Z. and Silverstein, J. W.  CLT for linear spectral statistics of large-dimensional sample covariance matrices.  _Annals of Probability_ , 32(1A):553 – 605, 2004\.
  * Bai and Saranadasa [1996] Bai, Z. and Saranadasa, H.  Effect of high dimension: by an example of a two sample problem.  _Statistica Sinica_ , pages 311–329, 1996.
  * Bai and Silverstein [2010] Bai, Z. and Silverstein, J. W.  _Spectral Analysis of Large Dimensional Random Matrices_.  Springer, 2010.
  * Bai and Zhou [2008] Bai, Z. and Zhou, W.  Large sample covariance matrices without independence structures in columns.  _Statistica Sinica_ , pages 425–442, 2008.
  * Bai et al. [2010] Bai, Z., Chen, J., and Yao, J.  On estimation of the population spectral distribution from a high-dimensional sample covariance matrix.  _Australian & New Zealand Journal of Statistics_, 52(4):423–437, 2010.
  * Boucheron et al. [2013] Boucheron, S., Lugosi, G., and Massart, P.  _Concentration Inequalities: A Nonasymptotic Theory of Independence_.  Oxford, 2013.
  * Couillet and Debbah [2011] Couillet, R. and Debbah, M.  _Random Matrix Methods for Wireless Communications_.  Cambridge, 2011.
  * Edelman and Wang [2013] Edelman, A. and Wang, Y.  Random matrix theory and its innovative applications.  In _Advances in Applied Mathematics, Modeling, and Computational Science_ , pages 91–116. Springer, 2013.
  * El Karoui [2009] El Karoui, N.  Concentration of measure and spectra of random matrices: Applications to correlation matrices, elliptical distributions and beyond.  _Annals of Applied Probability_ , 19(6):2362–2405, 2009.
  * El Karoui [2008] El Karoui, N.  Spectrum estimation for large dimensional covariance matrices using random matrix theory.  _Annals of Statistics_ , 36(6):2757–2790, 2008\.
  * Fang et al. [1990] Fang, K.-T., Kotz, S., and Ng, K. W.  _Symmetric Multivariate and Related Distributions_.  Chapman and Hall, 1990.
  * Gao et al. [2017] Gao, J., Han, X., Pan, G., and Yang, Y.  High dimensional correlation matrices: The central limit theorem and its applications.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 79(3):677–693, 2017.
  * Gupta et al. [2013] Gupta, A. K., Varga, T., and Bodnar, T.  _Elliptically Contoured Models in Statistics and Portfolio Theory_.  Springer, 2013.
  * Härdle and Simar [2019] Härdle, W. K. and Simar, L.  _Applied Multivariate Statistical Analysis_.  Springer, 2019.
  * Hu et al. [2019a] Hu, J., Li, W., Liu, Z., and Zhou, W.  High-dimensional covariance matrices in elliptical distributions with application to spherical test.  _Annals of Statistics_ , 47(1):527–555, 2019a.
  * Hu et al. [2019b] Hu, J., Li, W., and Zhou, W.  Central limit theorem for mutual information of large mimo systems with elliptically correlated channels.  _IEEE Transactions on Information Theory_ , 65(11):7168–7180, 2019b.
  * Johnson et al. [1994] Johnson, N. L., Kotz, S., and Balakrishnan, N.  _Continuous Univariate Distributions, Volume 1, (2nd ed.)_.  Wiley, 1994.
  * Johnson et al. [1995] Johnson, N. L., Kotz, S., and Balakrishnan, N.  _Continuous Univariate Distributions, Volume 2, (2nd ed.)_.  Wiley, 1995.
  * Johnson et al. [1985] Johnson, W. B., Schechtman, G., and Zinn, J.  Best constants in moment inequalities for linear combinations of independent and exchangeable random variables.  _Annals of Probability_ , pages 234–253, 1985.
  * Jonsson [1982] Jonsson, D.  Some limit theorems for the eigenvalues of a sample covariance matrix.  _Journal of Multivariate Analysis_ , 12(1):1–38, 1982.
  * Jun et al. [2022] Jun, W., Jiahui, X., Long, Y., and Wang, Z.  Tracy-Widom limit for the largest eigenvalue of high-dimensional covariance matrices in elliptical distributions.  _Bernoulli_ , 28(4):2941–2967, 2022.
  * Kong and Valiant [2017] Kong, W. and Valiant, G.  Spectrum estimation from samples.  _Annals of Statistics_ , 45(5):2218–2247, 2017\.
  * Ledoit and Wolf [2015] Ledoit, O. and Wolf, M.  Spectrum estimation: A unified framework for covariance matrix estimation and PCA in large dimensions.  _Journal of Multivariate Analysis_ , 139:360–384, 2015\.
  * Ledoit and Wolf [2017] Ledoit, O. and Wolf, M.  Numerical implementation of the QuEST function.  _Computational Statistics & Data Analysis_, 115:199–223, 2017.
  * Li and Yao [2018] Li, W. and Yao, J.  On structure testing for component covariance matrices of a high dimensional mixture.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 80(2):293–318, 2018.
  * Li et al. [2022] Li, W., Wang, Q., Yao, J., and Zhou, W.  On eigenvalues of a high-dimensional spatial-sign covariance matrix.  _Bernoulli_ , 28(1):606–637, 2022.
  * Lopes [2016] Lopes, M. E.  Unknown sparsity in compressed sensing: Denoising and inference.  _IEEE Transactions on Information Theory_ , 62(9):5145–5166, 2016.
  * Lopes [2022+] Lopes, M. E.  Central limit theorem and bootstrap approximation in high dimensions: Near $1/\sqrt{n}$ rates via implicit smoothing.  _to appear: Annals of Statistics (arXiv:2009.06004)_ , 2022+.
  * Lopes and Yao [2022] Lopes, M. E. and Yao, J.  A sharp lower-tail bound for Gaussian maxima with application to bootstrap methods in high dimensions.  _Electronic Journal of Statistics_ , 16(1):58–83, 2022.
  * Lopes et al. [2019] Lopes, M. E., Blandino, A., and Aue, A.  Bootstrapping spectral statistics in high dimensions.  _Biometrika_ , 106(4):781–801, 2019.
  * Mammen [1993] Mammen, E.  Bootstrap and wild bootstrap for high dimensional linear models.  _Annals of statistics_ , 21(1):255–285, 1993\.
  * Marcus et al. [2022] Marcus, A. W., Spielman, D. A., and Srivastava, N.  Interlacing families III: Sharper restricted invertibility estimates.  _Israel Journal of Mathematics_ , 247(2):519–546, 2022.
  * McNeil et al. [2011] McNeil, A., Frey, R., and Embrechts, P.  _Quantitative Risk Management._ Princeton, 2011.
  * Mestre [2008] Mestre, X.  Improved estimation of eigenvalues and eigenvectors of covariance matrices using their sample estimates.  _IEEE Transactions on Information Theory_ , 54(11):5113–5129, 2008.
  * Patterson et al. [2006] Patterson, N., Price, A. L., and Reich, D.  Population structure and eigenanalysis.  _PLoS genetics_ , 2(12):e190, 2006.
  * Potters and Bouchaud [2020] Potters, M. and Bouchaud, J.-P.  _A First Course in Random Matrix Theory: For Physicists, Engineers and Data Scientists_.  Cambridge, 2020.
  * Rao et al. [2008] Rao, N. R., Mingo, J. A., Speicher, R., and Edelman, A.  Statistical eigen-inference from large Wishart matrices.  _Annals of Statistics_ , 36(6):2850–2885, 2008\.
  * Siegrist [2017] Siegrist, K.  _Probability, Mathematical Statistics, and Stochastic Processes_.  LibreTexts, 2017.
  * Srivastava [2005] Srivastava, M. S.  Some tests concerning the covariance matrix in high dimensional data.  _Journal of the Japan Statistical Society_ , 35(2):251–272, 2005.
  * Tang and Nehorai [2012] Tang, G. and Nehorai, A.  The stability of low-rank matrix reconstruction: a constrained singular value view.  _IEEE Transactions on Information Theory_ , 58(9):6079–6092, 2012.
  * Vershynin [2018] Vershynin, R.  _High-Dimensional Probability: An Introduction with Applications in Data Science_.  Cambridge, 2018.
  * Yang et al. [2021] Yang, X., Zheng, X., and Chen, J.  Testing high-dimensional covariance matrices under the elliptical distribution and beyond.  _Journal of Econometrics_ , 221(2):409–423, 2021\.
  * Yao et al. [2015] Yao, J., Zheng, S., and Bai, Z.  _Sample Covariance Matrices and High-Dimensional Data Analysis_.  Cambridge, 2015.
  * Yin et al. [1988] Yin, Y. Q., Bai, Z., and Krishnaiah, P. R.  On the limit of the largest eigenvalue of the large dimensional sample covariance matrix.  _Probability Theory and Related Fields_ , 78(4):509–521, 1988.
  * Zhang et al. [2022] Zhang, Y., Hu, J., and Li, W.  CLT for linear spectral statistics of high-dimensional sample covariance matrices in elliptical distributions.  _Journal of Multivariate Analysis_ , page 105007, 2022.
  * Zhao and Candès [2022] Zhao, Q. and Candès, E. J.  An adaptively resized parametric bootstrap for inference in high-dimensional generalized linear models.  _arXiv:2208.08944_ , 2022.

## Appendices

In Appendices A, B, and C, we give the proofs of Theorems 1, 2, and 3
respectively. Background results are given in Appendix D, and additional
details related to the examples in Section 3.1 are given in Appendix E.

##  Appendix A Proof of Theorem 1

Our consistency guarantees for $\widehat{\varsigma}_{n}^{2}$ and
$\tilde{H}_{n}$ are proven separately in the next two subsections.

###  A.1 Consistency of $\widehat{\varsigma}_{n}^{\,2}$: Proof of (3.7) in
Theorem 1

Define the parameter

$\displaystyle\tau_{n}=\frac{(p+2)(\beta_{n}-2\alpha_{n})}{\gamma_{n}+2\alpha_{n}}+2,$
(A.1)

and the estimate

$\displaystyle\widehat{\tau}_{n}=\frac{(p+2)(\widehat{\beta}_{n}-2\widehat{\alpha}_{n})}{\widehat{\gamma}_{n}+2\widehat{\alpha}_{n}}+2.$
(A.2)

Based on the definitions of $\varsigma_{n}^{2}$ and
$\widehat{\varsigma}_{n}^{2}$, note that

$\textstyle\frac{\varsigma_{n}^{2}}{p}=\tau_{n}\ \ \ \text{ and }\ \ \ \
\textstyle\frac{\widehat{\varsigma}_{n}^{2}}{p}=\max\\{\widehat{\tau}_{n},0\\}.$

Since the function $x\mapsto\max\\{x,0\\}$ is 1-Lipschitz and $\tau_{n}\geq
0$, it suffices to show
$\widehat{\tau}_{n}-\tau_{n}\xrightarrow{\mathbb{P}}0$.

In Lemmas A.1, A.2, and A.3 given later in this subsection, the following
three limits will be established,

$\displaystyle\frac{\widehat{\alpha}_{n}}{\alpha_{n}}\
\xrightarrow{\mathbb{P}}\ 1$ (A.3)
$\displaystyle\frac{\widehat{\gamma}_{n}}{\gamma_{n}}\
\xrightarrow{\mathbb{P}}\ 1$ (A.4)
$\displaystyle(p+2)\frac{\widehat{\beta}_{n}-\beta_{n}}{\gamma_{n}+2\alpha_{n}}\
\xrightarrow{\mathbb{P}}\ 0.$ (A.5)

Due to the ratio-consistency of $\widehat{\alpha}_{n}$ and
$\widehat{\gamma}_{n}$, as well as the fact that
$\frac{(p+2)\alpha_{n}}{\gamma_{n}+2\alpha_{n}}\asymp 1$ holds under
Assumption 2, it is straightforward to check that

$\frac{\gamma_{n}+2\alpha_{n}}{\widehat{\gamma}_{n}+2\widehat{\alpha}_{n}}\xrightarrow{\mathbb{P}}1\
\ \ \ \text{ and }\ \ \ \
\frac{(p+2)\alpha_{n}}{\gamma_{n}+2\alpha_{n}}-\frac{(p+2)\widehat{\alpha}_{n}}{\widehat{\gamma}_{n}+2\widehat{\alpha}_{n}}\xrightarrow{\mathbb{P}}0.$
(A.6)

Therefore,

$\begin{split}\widehat{\tau}_{n}-\tau_{n}&=\frac{(p+2)\widehat{\beta}_{n}}{\gamma_{n}+2\alpha_{n}}\left(\frac{\gamma_{n}+2\alpha_{n}}{\widehat{\gamma}_{n}+2\widehat{\alpha}_{n}}-1\right)+(p+2)\frac{\widehat{\beta}_{n}-\beta_{n}}{\gamma_{n}+2\alpha_{n}}+\bigg{(}2\frac{(p+2)\alpha_{n}}{\gamma_{n}+2\alpha_{n}}-2\frac{(p+2)\widehat{\alpha}_{n}}{\widehat{\gamma}_{n}+2\widehat{\alpha}_{n}}\bigg{)}\\\
&=\frac{(p+2)\beta_{n}}{\gamma_{n}+2\alpha_{n}}o_{\mathbb{P}}(1)+o_{\mathbb{P}}(1),\end{split}$
(A.7)

where we have applied (A.5) twice in the second step. Under Assumption 1, we
have $\mathbb{E}(\xi_{1}^{4})=p^{2}+\tau p+o(p)$ and so Lemma D.1 gives

$\displaystyle\beta_{n}$ $\displaystyle\ =\
\textstyle\frac{\mathbb{E}(\xi_{1}^{4})}{p(p+2)}\Big{(}\textup{tr}(\Sigma_{n})^{2}+2\textup{tr}(\Sigma_{n}^{2})\Big{)}\
-\ \textup{tr}(\Sigma_{n})^{2}$ $\displaystyle\ =\ \textstyle\frac{p^{2}+\tau
p+o(p)}{p(p+2)}\Big{(}\textup{tr}(\Sigma_{n})^{2}+2\textup{tr}(\Sigma_{n}^{2})\Big{)}\
-\ \textup{tr}(\Sigma_{n})^{2}$ $\displaystyle\ =\
\Big{(}\textstyle\frac{(\tau-2)p+o(p)}{p(p+2)}\Big{)}\textup{tr}(\Sigma_{n})^{2}+2(1+o(1))\textup{tr}(\Sigma_{n}^{2})$
$\displaystyle\ \lesssim\ p.$

Consequently, we have $\frac{(p+2)\beta_{n}}{\gamma_{n}+2\alpha_{n}}\lesssim
1$, and applying this to (A.7) completes the proof of (3.7) in Theorem 1.∎

######  Lemma A.1.

If Assumptions 1 and 2 hold, then as $n\rightarrow\infty$

$\frac{\widehat{\alpha}_{n}}{\alpha_{n}}\ \xrightarrow{\mathbb{P}}\ 1.$

###### Proof.

Recall that $\alpha_{n}=\textup{tr}(\Sigma_{n}^{2})$ and that
$\widehat{\alpha}_{n}=\textup{tr}(\widehat{\Sigma}_{n}^{2})-\frac{1}{n}\textup{tr}(\widehat{\Sigma}_{n})^{2}$.
The two terms in the estimate can be expanded as

$\begin{split}\textup{tr}(\widehat{\Sigma}_{n}^{2})&\ =\
\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}(\mathbf{x}_{i}^{\top}\mathbf{x}_{j})^{2}\\\\[5.69046pt]
\frac{1}{n}\textup{tr}(\widehat{\Sigma}_{n})^{2}&\ =\
\frac{1}{n^{3}}\sum_{i=1}^{n}\sum_{j=1}^{n}(\mathbf{x}_{i}^{\top}\mathbf{x}_{i})(\mathbf{x}_{j}^{\top}\mathbf{x}_{j}),\end{split}$
(A.8)

which leads to the algebraic relation

$\begin{split}\frac{\widehat{\alpha}_{n}-\alpha_{n}}{\alpha_{n}}&\ =\
\bigg{(}\frac{1}{\alpha_{n}n^{2}}\sum_{i\neq
j}(\mathbf{x}_{i}^{\top}\mathbf{x}_{j})^{2}\ -\ 1\bigg{)}\ +\
\frac{1}{\alpha_{n}n^{3}}\sum_{i\neq
j}\Big{(}(\mathbf{x}_{i}^{\top}\mathbf{x}_{i})^{2}-(\mathbf{x}_{i}^{\top}\mathbf{x}_{i})(\mathbf{x}_{j}^{\top}\mathbf{x}_{j})\Big{)}\\\\[5.69046pt]
&\ =:\ A_{n}+B_{n}.\end{split}$ (A.9)

In the remainder of the proof, we will show that $A_{n}$ and $B_{n}$ are both
$o_{\mathbb{P}}(1)$. We begin with the analysis of $B_{n}$, since it is
simpler. Note that $B_{n}$ is always non-negative, since it can be rewritten
as

$B_{n}\ =\
\textstyle\frac{1}{\alpha_{n}n^{3}}\sum_{i>j}(\mathbf{x}_{i}^{\top}\mathbf{x}_{i}-\mathbf{x}_{j}^{\top}\mathbf{x}_{j})^{2},$
(A.10)

and so it suffices to show that $\mathbb{E}(B_{n})=o(1)$. Furthermore, the
expectation of $B_{n}$ can be computed directly as

$\begin{split}\mathbb{E}(B_{n})&\ =\
\textstyle\frac{1}{\alpha_{n}n^{3}}\sum_{i>j}2\,\operatorname{var}(\mathbf{x}_{i}^{\top}\mathbf{x}_{i})\\\\[5.69046pt]
&\ =\
\textstyle\frac{n(n-1)}{\alpha_{n}n^{3}}\operatorname{var}(\mathbf{x}_{1}^{\top}\mathbf{x}_{1})\\\\[5.69046pt]
&\ =\
\frac{n(n-1)}{\alpha_{n}n^{3}}\bigg{[}\Big{(}\textstyle\frac{\mathbb{E}(\xi_{1}^{4})}{p(p+2)}-1\Big{)}\textup{tr}(\Sigma_{n})^{2}\
+\ 2\frac{\mathbb{E}(\xi_{1}^{4})}{p(p+2)}\textup{tr}(\Sigma_{n}^{2})\bigg{]}\
\ \ \ \ \text{(Lemma \ref{lem:quadform})}\\\\[5.69046pt] &\ =\
\frac{n(n-1)}{\alpha_{n}n^{3}}\bigg{[}\mathcal{O}\big{(}\textstyle\frac{1}{p}\big{)}\textup{tr}(\Sigma_{n})^{2}\
+\
2\Big{(}1+\mathcal{O}\big{(}\textstyle\frac{1}{p}\big{)}\Big{)}\textup{tr}(\Sigma_{n}^{2})\bigg{]}\
\ \ \ \ \ \ \text{(Assumption \ref{Data generating model})}\\\\[5.69046pt] &\
\lesssim\ \textstyle\frac{1}{n},\end{split}$ (A.11)

where the last step uses the fact that $\alpha_{n}\asymp p$ holds under
Assumption 2. Thus, $\mathbb{E}(B_{n})=o(1)$, as needed.

Now we handle the term $A_{n}$ by showing that $\mathbb{E}(A_{n}^{2})=o(1)$.
It is helpful to start by noting that if $i\neq j$, then

$\mathbb{E}((\mathbf{x}_{i}^{\top}\mathbf{x}_{j})^{2})=\textup{tr}(\Sigma_{n}^{2}),$
(A.12)

which can be checked by a direct calculation. (See (D.5) for additional
details.) Consequently, if we expand out the square $A_{n}^{2}$ and then take
the expectation, it follows that

$\mathbb{E}(A_{n}^{2})\ =\
\frac{1}{\alpha_{n}^{2}n^{4}}\mathbb{E}\bigg{[}\Big{(}\sum_{i\neq
j}(\mathbf{x}_{i}^{\top}\mathbf{x}_{j})^{2}\Big{)}^{2}\bigg{]}\ -\
\frac{2n(n-1)}{n^{2}}\ +\ 1.$ (A.13)

Next, we compute the second moment on the right as

$\mathbb{E}\bigg{[}\Big{(}\sum_{i\neq
j}(\mathbf{x}_{i}^{\top}\mathbf{x}_{j})^{2}\Big{)}^{2}\bigg{]}\ =\
\displaystyle\sum_{i\neq j}\sum_{l\neq
k}\mathbb{E}\Big{[}(\mathbf{x}_{i}^{\top}\mathbf{x}_{j})^{2}(\mathbf{x}_{l}^{\top}\mathbf{x}_{k})^{2}\Big{]}.$
(A.14)

In Lemma D.2 , it is shown that the terms in the double sum can be controlled
as follows

$\displaystyle\mathbb{E}((\mathbf{x}_{i}^{\top}\mathbf{x}_{j})^{4})$
$\displaystyle\ \asymp\ \alpha_{n}^{2}\ \ \ \ \ \ \ \ \ \ \text{$(i,j)$
distinct}$ (A.15)
$\displaystyle\mathbb{E}((\mathbf{x}_{i}^{\top}\mathbf{x}_{j})^{2}(\mathbf{x}_{i}^{\top}\mathbf{x}_{k})^{2})$
$\displaystyle\ \asymp\ \alpha_{n}^{2}\ \ \ \ \ \ \ \ \ \ \text{$(i,j,k)$
distinct}$ (A.16)
$\displaystyle\mathbb{E}((\mathbf{x}_{i}^{\top}\mathbf{x}_{j})^{2}(\mathbf{x}_{l}^{\top}\mathbf{x}_{k})^{2})$
$\displaystyle\ =\ \alpha_{n}^{2}\ \ \ \ \ \ \ \ \ \ \text{$(i,j,k,l)$
distinct}.$ (A.17)

Note that in the double sum (A.14), the number of terms involving 4 distinct
indices is equal to $n(n-1)(n-2)(n-3)$, whereas the number of terms involving
either 2 or 3 distinct indices is $\mathcal{O}(n^{3})$. Applying these
observations to (A.13) we have

$\begin{split}\mathbb{E}(A_{n}^{2})&\ =\
\frac{1}{\alpha_{n}^{2}n^{4}}n(n-1)(n-2)(n-3)\alpha_{n}^{2}+\mathcal{O}\Big{(}\frac{n^{3}\alpha_{n}^{2}}{\alpha_{n}^{2}n^{4}}\Big{)}\
-\ \frac{2n(n-1)}{n^{2}}\ +\ 1\\\\[5.69046pt] &\ =\
\big{(}1+\mathcal{O}(\textstyle\frac{1}{n})\big{)}\ +\
\mathcal{O}(\textstyle\frac{1}{n})\ -\
2(1+\mathcal{O}(\textstyle\frac{1}{n}))\ +\ 1\\\\[5.69046pt] &\ \lesssim\
\textstyle\frac{1}{n},\end{split}$ (A.18)

which completes the proof. ∎

######  Lemma A.2.

If Assumptions 1 and 2 hold, then as $n\to\infty$

$\frac{\widehat{\gamma}_{n}}{\gamma_{n}}\ \xrightarrow{\mathbb{P}}\ 1.$

###### Proof.

Recall $\gamma_{n}=\textup{tr}(\Sigma_{n})^{2}$ and
$\widehat{\gamma}_{n}=\textup{tr}(\widehat{\Sigma}_{n})^{2}$, and note that
the algebraic identity $a^{2}-b^{2}=2b(a-b)+(a-b)^{2}$ gives

$\displaystyle\mathbb{E}|\widehat{\gamma}_{n}-\gamma_{n}|$ $\displaystyle\ =\
\mathbb{E}\Big{|}2\textup{tr}(\Sigma_{n})\big{(}\textup{tr}(\widehat{\Sigma}_{n})-\textup{tr}(\Sigma_{n})\big{)}+\big{(}\textup{tr}(\widehat{\Sigma}_{n})-\textup{tr}(\Sigma_{n})\big{)}^{2}\Big{|}$
$\displaystyle\ \leq\
2\textup{tr}(\Sigma_{n})\big{(}\operatorname{var}(\textup{tr}(\widehat{\Sigma}_{n}))\big{)}^{1/2}+\operatorname{var}(\textup{tr}(\widehat{\Sigma}_{n})).$

Next, observe that
$\operatorname{var}(\textup{tr}(\widehat{\Sigma}_{n}))=\frac{1}{n}\operatorname{var}(\mathbf{x}_{1}^{\top}\mathbf{x}_{1})$,
and that the calculation in (A.11) shows
$\operatorname{var}(\mathbf{x}_{1}^{\top}\mathbf{x}_{1})\lesssim p$. Combining
this with the fact that $\gamma_{n}\asymp p^{2}$ under Assumption 2, we have

$\frac{\mathbb{E}|\widehat{\gamma}_{n}-\gamma_{n}|}{\gamma_{n}}\ \lesssim\
\frac{\textup{tr}(\Sigma_{n})+1}{p^{2}}$ (A.19)

which leads to the stated result.∎

######  Lemma A.3.

If Assumptions 1 and 2 hold, then as $n\to\infty$,

$(p+2)\frac{\widehat{\beta}_{n}-\beta_{n}}{\gamma_{n}+2\alpha_{n}}\
\xrightarrow{\mathbb{P}}\ 0.$ (A.20)

###### Proof.

Recall that $\beta_{n}=\operatorname{var}(\|\mathbf{x}_{1}\|_{2}^{2})$ and
$\widehat{\beta}_{n}=\frac{1}{n-1}\sum_{i=1}^{n}\Big{(}\|\mathbf{x}_{i}\|_{2}^{2}-\textstyle\frac{1}{n}\sum_{i^{\prime}=1}^{n}\|\mathbf{x}_{i^{\prime}}\|_{2}^{2}\Big{)}^{2}$.
It is clear that $\mathbb{E}(\widehat{\beta}_{n})=\beta_{n}$, and so it
suffices to show

$\textstyle\frac{(p+2)^{2}}{(\gamma_{n}+2\alpha_{n})^{2}}\operatorname{var}(\widehat{\beta}_{n})=o(1).$

Since Assumption 2 implies $\gamma_{n}+2\alpha_{n}\asymp p^{2}$, it remains to
show $\operatorname{var}(\widehat{\beta}_{n})=o(p^{2})$.

A standard bound on the variance of $\widehat{\beta}_{n}$ is

$\operatorname{var}(\widehat{\beta}_{n})\ \lesssim\
\textstyle\frac{1}{n}\,\mathbb{E}\Big{|}\left\|\mathbf{x}_{1}\right\|_{2}^{2}-\textup{tr}(\Sigma_{n})\Big{|}^{4}.$
(A.21)

To bound the right side of (A.21), observe that

$\displaystyle\mathbb{E}\Big{|}\\!\left\|\mathbf{x}_{1}\right\|_{2}^{2}-\textup{tr}\left(\Sigma_{n}\right)\\!\Big{|}^{4}$
$\displaystyle\ =\
\mathbb{E}\left|\xi_{1}^{2}\left(\mathbf{u}_{1}^{\top}\Sigma_{n}\mathbf{u}_{1}-\textstyle\frac{1}{p}\textup{tr}(\Sigma_{n})\right)+\left(\textstyle\frac{1}{p}\xi_{1}^{2}-1\right)\textup{tr}(\Sigma_{n})\right|^{4}$
$\displaystyle\ \lesssim\
\mathbb{E}(\xi_{1}^{8})\cdot\mathbb{E}\left|\mathbf{u}_{1}^{\top}\Sigma_{n}\mathbf{u}_{1}-\textstyle\frac{1}{p}\textup{tr}(\Sigma_{n})\right|^{4}+\textup{tr}(\Sigma_{n})^{4}\cdot\mathbb{E}\left|\textstyle\frac{1}{p}\xi_{1}^{2}-1\right|^{4}.$

Due to Lemma 5 in [14] we have

$\mathbb{E}\left|\mathbf{u}_{1}^{\top}\Sigma_{n}\mathbf{u}_{1}-\textstyle\frac{1}{p}\textup{tr}(\Sigma_{n})\right|^{4}\
\lesssim\ \textstyle\frac{1}{p^{2}}.$

Also, Assumption 1 implies

$\mathbb{E}\left|\textstyle\frac{1}{p}\xi_{1}^{2}-1\right|^{4}\ \lesssim\
\textstyle\frac{1}{p^{2}}\ \ \ \ \ \text{ and }\ \ \ \ \
\mathbb{E}(\xi_{1}^{8})\,\lesssim\,p^{4}.$ (A.22)

Applying the last several observations to (A.21) implies
$\operatorname{var}(\widehat{\beta}_{n})\lesssim p$, which yields the stated
result. ∎

###  A.2 Consistency of $\tilde{H}_{n}$: Proof of (3.8) in Theorem 1

For each $t\in\mathbb{R}$, denote the empirical distribution function of the
QuEST eigenvalue estimates as

$\widehat{H}_{\textup{Q},n}(t)=\frac{1}{p}\sum_{j=1}^{p}1\\{\widehat{\lambda}_{\textup{Q},j}\leq
t\\}.$

It follows from Lemma A.4 below that the limit

$\widehat{H}_{\textup{Q},n}\Rightarrow H$ (A.23)

holds in probability as $n\to\infty$. So, to prove (3.8), it is sufficient to
show
$\sup_{t\in\mathbb{R}}|\widehat{H}_{\textup{Q},n}(t)-\tilde{H}_{n}(t)|\xrightarrow{\mathbb{P}}0$.

Since $\widehat{\lambda}_{\textup{Q},j}$ and $\tilde{\lambda}_{j}$ can only
disagree when $\widehat{\lambda}_{\textup{Q},j}>\widehat{b}_{n}$, we have

$\begin{split}\sup_{t\in\mathbb{R}}|\widehat{H}_{\textup{Q},n}(t)-\tilde{H}_{n}(t)|&\
\leq\
\frac{1}{p}\sum_{j=1}^{p}1\\{\widehat{\lambda}_{\textup{Q},j}\neq\tilde{\lambda}_{j}\\}\\\\[5.69046pt]
&\ =\
\frac{1}{p}\sum_{j=1}^{p}1\\{\widehat{\lambda}_{\textup{Q},j}>\widehat{b}_{n}\\}\\\\[5.69046pt]
&\ =\ 1-\widehat{H}_{\textup{Q},n}(\widehat{b}_{n}).\end{split}$ (A.24)

Let $u$ denote the upper endpoint of the support of the distribution
associated with $H$, and fix any $\epsilon\in(0,\frac{1}{4})$. Since
$\widehat{H}_{\textup{Q},n}$ is a non-decreasing function, we have

$\begin{split}\sup_{t\in\mathbb{R}}|\widehat{H}_{n}(t)-\tilde{H}_{n}(t)|&\
\leq\ 1-\widehat{H}_{\textup{Q},n}(u+\epsilon)1\\{\widehat{b}_{n}\geq
u+\epsilon\\}.\end{split}$ (A.25)

By the definition of $u$, the value $u+\epsilon$ is a continuity point of $H$
with $H(u+\epsilon)=1$, and so the limit (A.23) implies
$\widehat{H}_{\textup{Q},n}(u+\epsilon)=1+o_{\mathbb{P}}(1)$. Hence, in order
to show that the right side of (A.25) converges to 0 in probability, it is
enough to show that $\mathbb{P}(\widehat{b}_{n}\geq u+\epsilon)\to 1$.

We will handle this remaining task by showing instead that
$\mathbb{P}(\widehat{b}_{n}<u+\epsilon)\to 0$. Recall that
$\widehat{b}_{n}=\lambda_{1}(\widehat{\Sigma}_{n})+1$ and note that the limit
$H_{n}\Rightarrow H$ implies that $\lambda_{1}(\Sigma_{n})\geq u-\epsilon$
must hold for all large $n$. Therefore, we have

$\begin{split}\mathbb{P}(\widehat{b}_{n}<u+\epsilon)&\ \leq\
\mathbb{P}\Big{(}\lambda_{1}(\widehat{\Sigma}_{n})\leq\lambda_{1}(\Sigma_{n})+2\epsilon-1\Big{)}\end{split}$
(A.26)

for all large $n$. To derive an upper bound on the last probability, we will
replace $\lambda_{1}(\widehat{\Sigma}_{n})$ with a smaller random variable,
and then rearrange the event. Let $\mathbf{v}_{1}$ denote an eigenvector of
$\Sigma_{n}$ corresponding to $\lambda_{1}(\Sigma_{n})$ with
$\|\mathbf{v}_{1}\|_{2}=1$. Defining the random matrix
$W_{n}=\frac{1}{n}\sum_{i=1}^{n}\xi_{i}^{2}\mathbf{u}_{i}\mathbf{u}_{i}^{\top}$,
the variational representation of $\lambda_{1}(\widehat{\Sigma}_{n})$ gives
$\lambda_{1}(\widehat{\Sigma}_{n})\geq\lambda_{1}(\Sigma_{n})\mathbf{v}_{1}^{\prime}W_{n}\mathbf{v}_{1}$.
Also note that our choice of $\epsilon$ ensures $2\epsilon-1\leq-1/2$. This
yields the following bounds for all large $n$,

$\begin{split}\mathbb{P}(\widehat{b}_{n}<u+\epsilon)&\ \leq\
\mathbb{P}\Big{(}\lambda_{1}(\Sigma_{n})|\mathbf{v}_{1}^{\top}W_{n}\mathbf{v}_{1}-1|\geq
1/2\Big{)}\\\\[5.69046pt] &\ \leq\
4\lambda_{1}(\Sigma_{n})^{2}\operatorname{var}\big{(}\mathbf{v}_{1}^{\top}W_{n}\mathbf{v}_{1}\big{)}\\\\[5.69046pt]
&\ =\
\textstyle\frac{4\lambda_{1}(\Sigma_{n})^{2}}{n}\,\operatorname{var}\big{(}\xi_{1}^{2}\mathbf{u}_{1}^{\top}(\mathbf{v}_{1}\mathbf{v}_{1}^{\top})\mathbf{u}_{1}\big{)}\\\\[5.69046pt]
&\ =\
\textstyle\frac{4\lambda_{1}(\Sigma_{n})^{2}}{n}\,\left(3\frac{\mathbb{E}(\xi_{1}^{4})}{p(p+2)}-1\right)\
\ \ \ \ \ \ \ \ \ \ \text{(Lemma~{}\ref{lem:quadform})}\\\\[0.0pt] &\
\lesssim\ \textstyle\frac{1}{n},\end{split}$ (A.27)

where the last step uses Assumptions 1 and 2.∎

######  Lemma A.4.

Suppose that Assumptions 1 and 2 hold. Then, the following limit holds almost
surely as $n\to\infty$

$\widehat{H}_{\textup{Q},n}\Rightarrow H.$ (A.28)

###### Proof.

In the paper [25, pp.381-382], the limit (A.28) is established in the context
of an MP model where $p/n\to c\in(0,\infty)\setminus\\{1\\}$, and the
eigenvalues of the population covariance matrix $\Sigma_{n}$ satisfy
Assumption 2. To adapt the proof from [25] to our current setting, it is
sufficient to show that the following two facts hold under Assumptions 1 and
2. First, there is a constant $C>0$ such that the bound
$\limsup_{n\to\infty}\lambda_{1}(\widehat{\Sigma}_{n})\leq C$ holds almost
surely, which we prove later using a truncation argument in Lemmas A.5 and
A.6. Second, the random distribution function
$\widehat{H}_{n}(t)=\frac{1}{p}\sum_{j=1}^{p}1\\{\lambda_{j}(\widehat{\Sigma}_{n})\leq
t\\}$ satisfies $\widehat{H}_{n}\Rightarrow\Psi(H,c)$ almost surely, where we
recall that the distribution $\Psi(H,c)$ is defined near equation (3.9). The
validity of this second fact under our current assumptions is a consequence of
Theorem 1.1 in [6].∎

###  A.3 Boundedness of sample eigenvalues

The following lemma shows that the variable
$\widehat{b}_{n}=\lambda_{1}(\widehat{\Sigma}_{n})+1$ used in the definition
of $\tilde{\lambda}_{1},\dots,\tilde{\lambda}_{p}$ is asymptotically bounded.
This is needed to prove our bootstrap consistency result in Theorem 2, as well
as Lemma A.4 on the consistency of the QuEST method in the context of
elliptical models.

######  Lemma A.5.

Under Assumptions 1 and 2, there is a constant $C>0$ such that the bounds

$\displaystyle\limsup_{n\to\infty}\lambda_{1}({\widehat{\Sigma}_{n}})$
$\displaystyle\leq C\text{ \ \ \ \ \ \ and \ \ \ \ \ \
}\limsup_{n\to\infty}\lambda_{1}({\tilde{\Sigma}_{n}})\leq C$ (A.29)

hold almost surely.

###### Proof.

Since the estimates $\tilde{\lambda}_{j}=\lambda_{j}(\tilde{\Sigma}_{n})$ are
bounded above by $\widehat{b}_{n}=\lambda_{1}(\widehat{\Sigma}_{n})+1$ for all
$j\in\\{1,\dots,p\\}$, it is enough to focus on the first inequality in
(A.29). Define a sequence of truncated random variables
$\breve{\xi}_{i}=\xi_{i}\,1\\{|\xi_{i}^{2}-p|<\sqrt{pn}\\}$ for $i=1,\dots,n$,
as well as the following truncated version of $\widehat{\Sigma}_{n}$,

$\breve{\Sigma}_{n}=\frac{1}{n}\sum_{i=1}^{n}\breve{\xi}_{i}^{2}\,\Sigma_{n}^{1/2}\mathbf{u}_{i}\mathbf{u}_{i}^{\top}\Sigma_{n}^{1/2}.$
(A.30)

Lemma A.6 below shows that
$\mathbb{P}(\widehat{\Sigma}_{n}\neq\breve{\Sigma}_{n}\ \text{i.o.})=0$.
Consequently, it suffices to show there is a constant $C>0$ such that
$\limsup_{n\to\infty}\lambda_{1}(\breve{\Sigma}_{n})\leq C$ holds almost
surely.

Since the vectors $\mathbf{u}_{1},\dots,\mathbf{u}_{n}$ are uniformly
distributed on the unit sphere of $\mathbb{R}^{p}$, we may express them as
$\mathbf{u}_{i}=\mathbf{z}_{i}/\|\mathbf{z}_{i}\|_{2}$ for a sequence of
i.i.d. standard Gaussian vectors $\mathbf{z}_{1},\dots,\mathbf{z}_{n}$ in
$\mathbb{R}^{p}$. This yields

$\breve{\Sigma}_{n}=\frac{1}{n}\sum_{i=1}^{n}\frac{\breve{\xi}_{i}^{2}/p}{\|\mathbf{z}_{i}\|_{2}^{2}/p}\Sigma_{n}^{1/2}\mathbf{z}_{i}\mathbf{z}_{i}^{\top}\Sigma_{n}^{1/2}.$
(A.31)

By construction, we have $\max_{1\leq i\leq n}\breve{\xi}_{i}^{2}/p\leq
1+1/\sqrt{c_{n}}$ for all $n\geq 1$. Also, using standard tail bounds for the
$\chi_{p}^{2}$ distribution and the Borel-Cantelli lemma, it is
straightforward to show that $\liminf_{n\to\infty}\min_{1\leq i\leq
n}\|\mathbf{z}_{i}\|_{2}^{2}/p$ is at least $1/2$ almost surely. Taken
together, these observations imply there is a constant $C>0$ such that that
the bound

$\lambda_{1}(\breve{\Sigma}_{n})\ \leq\
C\,\lambda_{1}(\Sigma_{n})\,\lambda_{1}\Big{(}\textstyle\frac{1}{n}\sum_{i=1}^{n}\mathbf{z}_{i}\mathbf{z}_{i}^{\top}\Big{)}$
(A.32)

holds almost surely for all large $n$. In addition, note that $\sup_{n\geq
1}\lambda_{1}(\Sigma_{n})\lesssim 1$ holds by Assumption 2. Lastly, it is
known from [46, Theorem 3.1] that the limit

$\lim_{n\to\infty}\lambda_{1}\Big{(}\textstyle\frac{1}{n}\sum_{i=1}^{n}\mathbf{z}_{i}\mathbf{z}_{i}^{\top}\Big{)}\
=\ (1+\sqrt{c})^{2}$

holds almost surely, which completes the proof.∎

######  Lemma A.6.

Suppose that Assumption 1 holds, and let $\breve{\Sigma}_{n}$ be as defined in
(A.30). Then,

$\mathbb{P}(\widehat{\Sigma}_{n}\neq\breve{\Sigma}_{n}\ \textup{i.o.})=0.$
(A.33)

###### Proof.

We adapt a classical argument from [46]. For a fixed number $n$, the matrices
$\widehat{\Sigma}_{n}$ and $\breve{\Sigma}_{n}$ can only disagree if
$|\xi_{i}^{2}-p|\geq\sqrt{pn}$ for at least one $i=1,\dots,n$, and so

$\mathbb{P}(\widehat{\Sigma}_{n}\neq\breve{\Sigma}_{n}\ \textup{i.o.})\ \leq\
\lim_{j\to\infty}\mathbb{P}\bigg{(}\bigcup_{n=j}^{\infty}\bigcup_{i=1}^{n}\\{|\xi_{i}^{2}-p|\geq\sqrt{pn}\\}\bigg{)}.$
(A.34)

Next, we partition the values of $n$ into the intervals
$[2^{m-1},2^{m}),[2^{m},2^{m+1}),\dots$ and take a union bound across the
intervals, yielding

$\begin{split}\mathbb{P}(\widehat{\Sigma}_{n}\neq\breve{\Sigma}_{n}\
\textup{i.o.})&\ \leq\
\lim_{j\to\infty}\sum_{m=j}^{\infty}\mathbb{P}\Bigg{(}\max_{\underset{2^{m-1}\leq
n<2^{m}}{1\leq i\leq
n}}|\textstyle\frac{\xi_{i}^{2}-p}{\sqrt{p}}|\,\geq\,\sqrt{2^{m-1}}\Bigg{)}.\end{split}$
(A.35)

For a generic sequence of random variables $Y_{1},\dots,Y_{N}$ and number
$q\geq 1$, recall the standard maximal inequality

$\mathbb{E}\Big{|}\\!\max_{1\leq i\leq N}Y_{i}\Big{|}^{q}\ \leq\ N\max_{1\leq
i\leq N}\mathbb{E}|Y_{i}|^{q}.$ (A.36)

Since the number of pairs $(i,n)$ in the maximum in (A.35) is at most
$2^{2m}$, if we apply Chebyshev’s inequality to (A.35) and use the condition
$\mathbb{E}|\textstyle\frac{\xi_{1}^{2}-p}{\sqrt{p}}|^{4+\varepsilon}\lesssim
1$ from Assumption 1, then for each $m\geq 1$ we have

$\begin{split}\mathbb{P}\Bigg{(}\max_{\underset{2^{m-1}\leq n\leq 2^{m}}{1\leq
i\leq
n}}|\textstyle\frac{\xi_{i}^{2}-p}{\sqrt{p}}|\,\geq\,\sqrt{2^{m-1}}\Bigg{)}&\
\leq\
\frac{2^{2m}}{\big{(}\sqrt{2^{m-1}}\big{)}^{4+\varepsilon}}\\\\[5.69046pt] &\
\lesssim\ 2^{-m\varepsilon/2}.\end{split}$ (A.37)

Hence, we may insert this bound into (A.35) to conclude that

$\mathbb{P}(\widehat{\Sigma}_{n}\neq\breve{\Sigma}_{n}\ \textup{i.o.})\
\lesssim\ \lim_{j\to\infty}\sum_{m=j}^{\infty}2^{-m\varepsilon/2}\ =\ 0,$

which completes the proof.∎

##  Appendix B Proof of Theorem 2

Let $\mathbf{z}\in\mathbb{R}^{k}$ denote a Gaussian random vector to be
described in a moment, and consider the triangle inequality

$d_{\mathrm{LP}}\Big{(}\mathcal{L}(p\big{\\{}T_{n}(\mathbf{f})-\vartheta_{n}(\mathbf{f})\big{\\}}\big{)}\,,\,\mathcal{L}\big{(}p\big{\\{}T_{n,1}^{*}(\mathbf{f})-\tilde{\vartheta}_{n}(\mathbf{f})\big{\\}}\big{|}X\big{)}\Big{)}\
\leq\ \textup{I}_{n}\ +\ \textup{II}_{n},$ (B.1)

where we define

$\displaystyle\textup{I}_{n}$ $\displaystyle\ =\
d_{\mathrm{LP}}\Big{(}\mathcal{L}\big{(}p\big{\\{}T_{n}(\mathbf{f})-\vartheta_{n}(\mathbf{f})\big{\\}}\big{)},\mathcal{L}(\mathbf{z})\Big{)}$
(B.2) $\displaystyle\textup{II}_{n}$ $\displaystyle\ =\
d_{\mathrm{LP}}\Big{(}\mathcal{L}(\mathbf{z})\,,\,\mathcal{L}\big{(}p\big{\\{}T_{n,1}^{*}(\mathbf{f})-\tilde{\vartheta}_{n}(\mathbf{f})\\}\big{|}X\big{)}\Big{)}.$
(B.3)

To handle the terms $\textup{I}_{n}$ and $\textup{II}_{n}$, we will apply a
central limit theorem for linear spectral statistics established in [17],
which relies on the following two conditions when $c_{n}\to c>0$ as
$n\to\infty$:

  * (a).

The elliptical model in Assumption 1 holds with the conditions in (3.2) being
replaced by

$\mathbb{E}(\xi_{1}^{2})=p,\ \ \ \
\operatorname{var}\Big{(}\textstyle\frac{\xi_{1}^{2}-p}{\sqrt{p}}\Big{)}=\tau+o(1)\
\ \ \ \ \text{ and }\ \ \ \
\mathbb{E}\Big{|}\textstyle\frac{\xi_{1}^{2}-p}{\sqrt{p}}\Big{|}^{2+\varepsilon}\lesssim
1,$ (B.4)

for some constants constants $\tau\geq 0$ and $\varepsilon>0$ that do not
depend on $n$.

  * (b).

There is a distribution $H$ such that $H_{n}\Rightarrow H$ as $n\to\infty$,
and $\lambda_{1}(\Sigma_{n})\lesssim 1$.

Under these conditions, Theorem 2.2 in [17] ensures there exists a Gaussian
distribution $\mathcal{L}(\mathbf{z})$ depending only on
$(\mathbf{f},H,c,\tau)$ such that $\textup{I}_{n}\to 0$ as $n\to\infty$.

To finish the proof, we must show $\textup{II}_{n}\xrightarrow{\mathbb{P}}0$.
This can be done using the mentioned central limit theorem, but some extra
considerations are involved, due to the fact that $\textup{II}_{n}$ is random.
It is sufficient to show that for any subsequence $J\subset\\{1,2,\dots\\}$,
the limit $\textup{II}_{n}\to 0$ holds almost surely along a further
subsequence of $J$. Since the bootstrap data are generated from an elliptical
model that is parameterized in terms of $\mathcal{L}(\xi_{1}^{*}|X)$ and
$\tilde{\Sigma}_{n}$, this amounts to verifying that
$\mathcal{L}(\xi_{1}^{*}|X)$ and $\tilde{H}_{n}$ satisfy analogues of (a) and
(b) almost surely along a subsequence of $J$.

To proceed, recall that Algorithm 1 is designed so that
$\mathbb{E}((\xi_{1}^{*})^{2}|X)=p$ and
$\operatorname{var}((\xi_{1}^{*})^{2}|X)=\widehat{\varsigma}_{n}^{2}$. Also,
note that under Assumption 1, the parameter
$\varsigma_{n}^{2}=\operatorname{var}(\xi_{1}^{2})$ satisfies
$\varsigma_{n}^{2}/p\to\tau$. Consequently, Theorem 1 implies
$\widehat{\varsigma}_{n}^{2}/p=\tau+o_{\mathbb{P}}(1)$, and so there is a
subsequence $J^{\prime}\subset J$ along which the limit
$\widehat{\varsigma}_{n}^{2}/p\to\tau$ holds almost surely. In other words,
the limit

$\operatorname{var}\\!\Big{(}\textstyle\frac{(\xi_{1}^{*})^{2}-p}{\sqrt{p}}\Big{|}X\Big{)}\
\to\ \tau$ (B.5)

holds almost surely along $J^{\prime}$. Moreover, since
$\mathcal{L}((\xi_{1}^{*})^{2}|X)$ is a Gamma distribution with mean $p$ and
variance $\widehat{\varsigma}_{n}^{2}$, Lemma E.4 implies that the bound

$\Big{(}\mathbb{E}\Big{(}\big{|}\textstyle\frac{(\xi_{1}^{*})^{2}-p}{\sqrt{p}}\big{|}^{6}\Big{|}X\Big{)}\Big{)}^{1/6}\leq
C\big{(}\textstyle\frac{\widehat{\varsigma}_{n}}{\sqrt{p}}+\textstyle\frac{\widehat{\varsigma}_{n}^{2}}{p^{3/2}}\big{)}$
(B.6)

holds almost surely and so the left side of (B.6) is bounded almost surely
along $J^{\prime}$.

With regard to the empirical spectral distribution $\tilde{H}_{n}$ associated
with the matrix $\tilde{\Sigma}_{n}$ in Algorithm 1, it satisfies the limit
$\tilde{H}_{n}\Rightarrow H$ almost surely along a further subsequence
$J^{\prime\prime}\subset J^{\prime}$, due to Theorem 1. In addition, Lemma A.5
ensures there is a constant $C>0$ such that
$\limsup_{n\to\infty}\lambda_{1}(\tilde{\Sigma}_{n})\leq C$ almost surely.
Altogether, it follows that $\mathcal{L}(\xi_{1}^{*}|X)$ and $\tilde{H}_{n}$
simultaneously satisfy analogues of (a) and (b) almost surely along
$J^{\prime\prime}$, which implies $\textup{II}_{n}\xrightarrow{\mathbb{P}}0$.∎

##  Appendix C Proof of Theorem 3

The proof is based on viewing $\widehat{r}_{n}$ as a nonlinear function of
$T_{n}(\mathbf{f})$ when the components of $\mathbf{f}=(f_{1},f_{2})$ are
taken to be $f_{1}(x)=x$ and $f_{2}(x)=x^{2}$. In this case, the centering
parameter $\vartheta_{n}(\mathbf{f})$ defined by (3.10) reduces to

$\vartheta_{n}(\mathbf{f})\ =\
\left(\textstyle\frac{1}{p}\textup{tr}(\Sigma_{n})\,,\,\textstyle\frac{1}{p}\textup{tr}(\Sigma_{n}^{2})+\frac{1}{np}\textup{tr}(\Sigma_{n})^{2}\right),$
(C.1)

as established in Lemma 2.16 of [45]. Consequently, by considering the
function

$g_{n}(x_{1},x_{2})=\textstyle\frac{x_{1}^{2}}{x_{2}-c_{n}x_{1}^{2}},$ (C.2)

we have the key relation

$\widehat{r}_{n}-r_{n}\ =\
p\big{\\{}g_{n}(T_{n}(\mathbf{f}))-g_{n}(\vartheta_{n}(\mathbf{f}))\big{\\}}.$
(C.3)

We can also develop a corresponding relation for the bootstrap statistic
$\widehat{r}_{n}^{*}-\tilde{r}$, defined by

$\widehat{r}_{n}^{*}\ =\
\textstyle\frac{\textup{tr}(\widehat{\Sigma}^{*}_{n})^{2}}{\textup{tr}((\widehat{\Sigma}^{*}_{n})^{2})-\frac{1}{n}\textup{tr}(\widehat{\Sigma}^{*}_{n})^{2}}\quad\text{
\ \ and \ \ }\quad\tilde{r}_{n}\ =\
\textstyle\frac{\textup{tr}(\tilde{\Sigma}_{n})^{2}}{\textup{tr}(\tilde{\Sigma}_{n}^{2})},$
(C.4)

where both quantities are arbitrarily taken to be equal to $n$ when their
denominators are 0. Due to Lemma 2.16 in [45] and our choice of $\mathbf{f}$,
the definition of $\tilde{\vartheta}_{n}(\mathbf{f})$ below (3.10) implies

$\tilde{\vartheta}_{n}(\mathbf{f})\ =\
\left(\textstyle\frac{1}{p}\textup{tr}(\tilde{\Sigma}_{n})\,,\,\textstyle\frac{1}{p}\textup{tr}(\tilde{\Sigma}_{n}^{2})+\frac{1}{np}\textup{tr}(\tilde{\Sigma}_{n})^{2}\right).$
(C.5)

Likewise, the bootstrap version of (C.3) is

$\widehat{r}_{n}^{*}-\tilde{r}_{n}\ =\
p\big{\\{}g_{n}(T_{n}^{*}(\mathbf{f}))-g_{n}(\tilde{\vartheta}_{n}(\mathbf{f}))\big{\\}}.$
(C.6)

Using the relations developed above, the next phase of the proof shows that
the distributions $\mathcal{L}(\widehat{r}_{n}-r_{n})$ and
$\mathcal{L}(\widehat{r}_{n}^{*}-\tilde{r}_{n}|X)$ approach the same limit, as
recorded below.

######  Proposition 1.

Under Assumptions 1 and 2, there is a non-degenerate Gaussian random variable
$\zeta$ such that the following limits hold as $n\to\infty$

$\displaystyle\mathcal{L}(\widehat{r}_{n}-r_{n})$ $\displaystyle\ \Rightarrow\
\mathcal{L}(\zeta)$ (C.7)
$\displaystyle\mathcal{L}(\widehat{r}_{n}^{*}-\tilde{r}_{n}|X)$
$\displaystyle\ \xRightarrow{\mathbb{P}}\ \mathcal{L}(\zeta).$ (C.8)

This proposition will be proved at the end of the section. Since the limiting
distribution $\mathcal{L}(\zeta)$ is continuous, Pólya’s Theorem implies

$\displaystyle\sup_{t\in\mathbb{R}}\Big{|}\mathbb{P}(\widehat{r}_{n}-r_{n}\leq
t)-\mathbb{P}(\widehat{r}_{n}^{*}-\tilde{r}_{n}\leq t|X)\Big{|}\
\xrightarrow{\mathbb{P}}\ 0.$ (C.9)

Due to this uniform limit and the continuity of the distribution
$\mathcal{L}(\zeta)$, standard arguments can be used to show that the
quantiles of $\mathcal{L}(\widehat{r}_{n}^{*}-\tilde{r}_{n}|X)$ are
asymptotically equivalent to those of $\mathcal{L}(\widehat{r}_{n}-r_{n})$.
(For example, see the proof of Lemma 10.4 in [30].) More precisely, if we note
that the quantile estimate $\widehat{q}_{1-\alpha}$ defined in Section 3.4 is
the same as the $(1-\alpha)$-quantile of
$\mathcal{L}\big{(}\frac{1}{p}(\widehat{r}_{n}^{*}-\tilde{r}_{n})|X\big{)}$,
then as $n\to\infty$ we have the limit

$\mathbb{P}\big{(}\textstyle\frac{1}{p}(\widehat{r}_{n}-r_{n})>\widehat{q}_{1-\alpha}\big{)}\to\alpha.$
(C.10)

This limit directly implies the first two statements (3.19) and (3.20) in
Theorem 3. Regarding the third statement (3.21), if
$\mathsf{H}_{0,n}^{\prime}$ holds for all large $n$, then the use of
$\tilde{\Sigma}_{n}=I_{p}$ instead of
$\tilde{\Sigma}_{n}=\textup{diag}(\tilde{\lambda}_{1},\dots,\tilde{\lambda}_{p})$
does not change argument leading up (C.10). Consequently, if
$\mathsf{H}_{0,n}^{\prime}$ holds for all large $n$, then we have
$\mathbb{P}\big{(}\frac{1}{p}(\widehat{r}_{n}-1)\leq\widehat{q}_{\alpha}^{\prime}\big{)}\to\alpha$
as $n\to\infty$, which completes the proof.∎

Proof of Proposition 1. As in the preceding paragraphs, we continue to take
$\mathbf{f}=(f_{1},f_{2})$ with $f_{1}(x)=x$, and $f_{2}(x)=x^{2}$. To begin
the argument, note that the proof of Theorem 2 shows there exists a Gaussian
random vector $\mathbf{z}\in\mathbb{R}^{2}$ such that as $n\to\infty$,

$\displaystyle\mathcal{L}\big{(}p\big{\\{}T_{n}(\mathbf{f})-\vartheta_{n}(\mathbf{f})\big{\\}}\big{)}$
$\displaystyle\ \Rightarrow\ \mathcal{L}(\mathbf{z})$ (C.11)
$\displaystyle\mathcal{L}\big{(}p\big{\\{}T_{n,1}^{*}(\mathbf{f})-\tilde{\vartheta}_{n}(\mathbf{f})\\}\big{|}X\big{)}$
$\displaystyle\ \xRightarrow{\mathbb{P}}\ \mathcal{L}(\mathbf{z}).$ (C.12)

Next, to introduce some notation, we refer to the $j$th moment of the
distribution $H$ as

$\phi_{j}=\textstyle\int t^{j}dH(t),$ (C.13)

and we use these moments to define the parameter

$\vartheta(\mathbf{f})=(\phi_{1},\phi_{2}+c\phi_{1}^{2}).$ (C.14)

This parameter arises in the following limits as $n\to\infty$,

$\vartheta_{n}(\mathbf{f})\to\vartheta(\mathbf{f})\ \ \ \ \ \ \text{ and }\ \
\ \ \ \
\tilde{\vartheta}_{n}(\mathbf{f})\xrightarrow{\mathbb{P}}\vartheta(\mathbf{f}).$
(C.15)

To see why these limits hold, first note that by Assumption 2, there is a
compact interval containing support of $H_{n}$ for all large $n$, and by
Theorem 1 as well as Lemma A.5, the same statement holds almost surely for
$\tilde{H}_{n}$. Moreover, Assumption 2 and Theorem 1 ensure the limits
$H_{n}\Rightarrow H$ and $\tilde{H}_{n}\xRightarrow{\mathbb{P}}H$, and so it
follows that the moments of $H_{n}$ and $\tilde{H}_{n}$ converge to those of
$H$. Then, combining these facts with the formulas (C.1) and (C.5) yields the
limits (C.15).

In light of the relations (C.3) and (C.6), and the limits (C.15), we will
expand the function $g_{n}$ around $\vartheta(\mathbf{f})$ and apply the delta
method. This is justified because the gradient $\nabla g_{n}$ has the
following continuity property. Namely, if we let
$g(x_{1},x_{2})=\frac{x_{1}^{2}}{x_{2}-cx_{1}^{2}}$, let $\mathcal{U}$ be a
sufficiently small open neighborhood of $\vartheta(\mathbf{f})$, and let
$\\{(x_{1,n},x_{2,n})\\}$ be any sequence of points within $\mathcal{U}$ that
converges to $\vartheta(\mathbf{f})$, then we have the limit

$\nabla g_{n}(x_{1,n},x_{2,n})\ \to\ \nabla g(\vartheta(\mathbf{f}))\ =\
\left(\textstyle\frac{2\phi_{1}}{\phi_{2}}+\frac{2c\phi_{1}^{3}}{\phi_{2}^{2}},-\frac{\phi_{1}^{2}}{\phi_{2}^{2}}\right).$
(C.16)

So, applying the delta method to the relations (C.3) and (C.6) and the weak
limits (C.11) and (C.12) gives

$\displaystyle\mathcal{L}(\widehat{r}_{n}-r_{n})$ $\displaystyle\ \Rightarrow\
\mathcal{L}(\nabla g(\vartheta(\mathbf{f}))^{\top}\mathbf{z})$ (C.17)
$\displaystyle\mathcal{L}(\widehat{r}_{n}^{*}-\tilde{r}_{n}\big{|}X\big{)}$
$\displaystyle\ \xRightarrow{\mathbb{P}}\ \mathcal{L}(\nabla
g(\vartheta(\mathbf{f}))^{\top}\mathbf{z}).$ (C.18)

Now, it remains to show that the variance of the Gaussian random variable
$\nabla g(\vartheta(\mathbf{f}))^{\top}\mathbf{z}$ is positive. Letting $K$
denote the $2\times 2$ covariance matrix of $\mathbf{z}$, it is shown below
equation 2.10 in [17] that the entries of $K$ are given by

$\displaystyle K_{11}$ $\displaystyle=2c\phi_{2}+c(\tau-2)\phi_{1}^{2},$
$\displaystyle K_{12}$
$\displaystyle=4c\phi_{3}+4c^{2}\phi_{1}\phi_{2}+2c(\tau-2)\phi_{1}(c\phi_{1}^{2}+\phi_{2})$
$\displaystyle K_{22}$
$\displaystyle=8c\phi_{4}+4c^{2}\phi_{2}^{2}+16c^{2}\phi_{1}\phi_{3}+8c^{3}\phi_{1}^{2}\phi_{2}+4c(\tau-2)(c\phi_{1}^{2}+\phi_{2})^{2}.$

Combining this with the formula for $\nabla g(\vartheta(\mathbf{f}))$ in
(C.16), a direct but lengthy computation of the quadratic form $\nabla
g(\vartheta(\mathbf{f}))^{\top}K\,\nabla g(\vartheta(\mathbf{f}))$ yields

$\begin{split}\operatorname{var}\\!\Big{(}\nabla
g(\vartheta(\mathbf{f}))^{\top}\mathbf{z}\Big{)}&\ =\
\textstyle\frac{4c\phi_{1}^{2}}{\phi_{2}^{4}}\Big{(}c\phi_{1}^{2}\phi_{2}^{2}+2\phi_{1}^{2}\phi_{4}+2\phi_{2}^{3}-4\phi_{1}\phi_{2}\phi_{3}\Big{)}.\end{split}$
(C.19)

To see that the variance is positive, it suffices to check that
$2\phi_{1}^{2}\phi_{4}+2\phi_{2}^{3}-4\phi_{1}\phi_{2}\phi_{3}$ is non-
negative. Since this quantity can be rewritten as
$2(\phi_{1}\phi_{4}^{1/2}-\phi_{2}^{3/2})^{2}+4\phi_{1}\phi_{2}(\phi_{4}^{1/2}\phi_{2}^{1/2}-\phi_{3})$,
its non-negativity follows from the observation that
$\phi_{3}\leq\phi_{4}^{1/2}\phi_{2}^{1/2}$ is an instance of the Cauchy-
Schwarz inequality. ∎

##  Appendix D Background results

######  Lemma D.1 ([17], Lemma A.1).

Let $\xi_{1}\in\mathbb{R}$ and $\mathbf{u}_{1}\in\mathbb{R}^{p}$ satisfy the
conditions in Assumption 1, and fix any $C\in\mathbb{R}^{p\times p}$. Then,

$\displaystyle\operatorname{var}(\xi_{1}^{2}\mathbf{u}_{1}^{\top}C\mathbf{u}_{1})=\frac{\mathbb{E}(\xi_{1}^{4})}{p(p+2)}(\textup{tr}(C)^{2}+2\textup{tr}(C^{2}))-\textup{tr}(C)^{2}.$
(D.1)

######  Lemma D.2.

Let $\mathbf{x}_{1},\ldots,\mathbf{x}_{n}$ be the observations defined in
Assumption 1, and let $i,j,k,l$ be four distinct indices in $\\{1,\dots,n\\}$.
Then,

$\displaystyle\mathbb{E}((\mathbf{x}_{i}^{\top}\mathbf{x}_{j})^{4})\ $
$\displaystyle\asymp\ \textup{tr}(\Sigma^{2}_{n})^{2}$ (D.2)
$\displaystyle\mathbb{E}((\mathbf{x}_{i}^{\top}\mathbf{x}_{j})^{2}(\mathbf{x}_{i}^{\top}\mathbf{x}_{k})^{2})\
$ $\displaystyle\asymp\ \textup{tr}(\Sigma^{2}_{n})^{2}$ (D.3)
$\displaystyle\mathbb{E}((\mathbf{x}_{i}^{\top}\mathbf{x}_{j})^{2}(\mathbf{x}_{l}^{\top}\mathbf{x}_{k})^{2})\
$ $\displaystyle=\ \textup{tr}(\Sigma^{2}_{n})^{2}.$ (D.4)

###### Proof.

Since the observations are i.i.d. and centered, we have

$\begin{split}\mathbb{E}((\mathbf{x}_{i}^{\top}\mathbf{x}_{j})^{2}(\mathbf{x}_{l}^{\top}\mathbf{x}_{k})^{2})&=(\operatorname{var}(\mathbf{x}_{i}^{\top}\mathbf{x}_{j}))^{2}\\\\[5.69046pt]
&=(\text{cov}(\mathbf{x}_{i}^{\top}\mathbf{x}_{j}\,,\,\mathbf{x}_{i}^{\top}\mathbf{x}_{j}))^{2}\\\
&=\Big{(}\sum_{r=1}^{p}\sum_{s=1}^{p}\mathbb{E}(\mathbf{x}_{ir}\mathbf{x}_{jr}\mathbf{x}_{is}\mathbf{x}_{js})\Big{)}^{2}\\\\[5.69046pt]
&=\Big{(}\sum_{r=1}^{p}\sum_{s=1}^{p}((\Sigma_{n})_{rs})^{2}\Big{)}^{2},\end{split}$
(D.5)

which establishes (D.4). The two other assertions can be shown using
conditional expectation. For the statement (D.2), we use Lemma D.1 to obtain

$\displaystyle\small\mathbb{E}((\mathbf{x}_{i}^{\top}\mathbf{x}_{j})^{4})$
$\displaystyle\ =\
\mathbb{E}(\mathbb{E}((\mathbf{x}_{i}^{\top}\mathbf{x}_{j}\mathbf{x}_{j}^{\top}\mathbf{x}_{i})^{2}|\mathbf{x}_{j}))$
$\displaystyle\ =\
\frac{\mathbb{E}(\xi_{1}^{4})}{p(p+2)}\mathbb{E}\bigg{(}\textstyle\big{(}\textup{tr}(\Sigma_{n}^{1/2}\mathbf{x}_{j}\mathbf{x}_{j}^{\top}\Sigma_{n}^{1/2})\big{)}^{2}+2\textup{tr}\Big{(}\big{(}\Sigma_{n}^{1/2}\mathbf{x}_{j}\mathbf{x}_{j}^{\top}\Sigma_{n}^{1/2}\big{)}^{2}\Big{)}\bigg{)}$
$\displaystyle\ =\
3\textstyle\frac{\mathbb{E}(\xi_{1}^{4})}{p(p+2)}\mathbb{E}((\mathbf{x}_{j}^{\top}\Sigma_{n}\mathbf{x}_{j})^{2})$
$\displaystyle\ =\
3\left(\frac{\mathbb{E}(\xi_{1}^{4})}{p(p+2)}\right)^{2}(\textup{tr}(\Sigma_{n}^{2})^{2}+2\textup{tr}(\Sigma_{n}^{4}))\
\ \ \ \ \ \ \ \ \ \ \ \ \ $ $\displaystyle\ \asymp\
\textup{tr}(\Sigma_{n}^{2})^{2}.$

The argument for (D.3) is similar. ∎

##  Appendix E Discussion of examples in Section 3.1

This section provides detailed information related to examples of the random
variable $\xi_{1}^{2}$ stated in Section 3.1. We give explicit
parameterizations, and we check that the distributions satisfy the conditions
in Assumption 1. The only three examples we do not individually cover are the
Chi-Squared, Poisson, and Negative-Binomial distributions, because they can be
represented in the form $\xi_{1}^{2}=\sum_{j=1}^{p}z_{1j}^{2}$ for independent
random variables $z_{11},\dots,z_{1p}$, and such examples are collectively
addressed by the following lemma.

######  Lemma E.1.

Suppose $\xi_{1}^{2}=\sum_{j=1}^{p}z_{1j}^{2}$ for some independent random
variables $z_{11},\dots,z_{1p}$ satisfying

$\textstyle\frac{1}{p}\sum_{j=1}^{p}\mathbb{E}(z_{1j}^{2})=1,\quad\textstyle\frac{1}{p}\sum_{j=1}^{p}\operatorname{var}(z_{1j}^{2})=\tau+o(1)\text{\quad
and \quad}\displaystyle\max_{1\leq j\leq
p}\mathbb{E}|z_{1j}|^{8+2\varepsilon}\,\lesssim\,p^{1+\frac{\varepsilon}{4}}$
(E.1)

as $n\to\infty$, for some fixed constants $\tau\geq 0$ and $\varepsilon>0$ not
depending on $n$. Then, $\xi_{1}$ satisfies the conditions in Assumption 1.

###### Proof.

It is only necessary to show
$\mathbb{E}\big{|}(\xi_{1}^{2}-\mathbb{E}(\xi_{1}^{2}))/\sqrt{p}\big{|}^{4+\varepsilon}\lesssim
1$. Using Rosenthal’s inequality [21] to bound the $L^{4+\varepsilon}$ norm of
a sum of independent centered random variables, we have

$\begin{split}\Big{(}\mathbb{E}\big{|}\xi_{1}^{2}-\mathbb{E}(\xi_{1}^{2})\big{|}^{4+\varepsilon}\Big{)}^{1/(4+\varepsilon)}&\
=\
\Big{(}\mathbb{E}\big{|}\sum_{j=1}^{p}z_{1j}^{2}-\mathbb{E}(z_{1j}^{2})\big{|}^{4+\varepsilon}\Big{)}^{1/(4+\varepsilon)}\\\\[5.69046pt]
&\ \lesssim\
\sqrt{\textstyle\sum_{j=1}^{p}\operatorname{var}(z_{1j}^{2})}+p^{1/(4+\varepsilon)}\max_{1\leq
j\leq
p}(\mathbb{E}|z_{1j}^{2}-\mathbb{E}(z_{1j}^{2})|^{4+\varepsilon})^{1/(4+\varepsilon)}\\\\[5.69046pt]
&\ \lesssim\sqrt{p}+p^{1/4}\max_{1\leq j\leq
p}(\mathbb{E}|z_{1j}|^{8+2\varepsilon})^{1/(4+\varepsilon)}\\\\[5.69046pt] &\
\lesssim\ \sqrt{p},\end{split}$ (E.2)

which completes the proof. ∎

Beta distribution The Beta$(a,b)$ distribution with parameters $a,b>0$ has a
density function that is proportional to $x^{a-1}(1-x)^{b-1}$ for $x\in(0,1)$.

######  Lemma E.2.

If $\beta>0$ is fixed with respect to $n$ and
$\xi_{1}^{2}\sim(p+2\beta)$Beta$(p/2,\beta)$, then $\xi_{1}$ satisfies the
conditions in Assumption 1.

###### Proof.

The mean and variance are given by $\mathbb{E}(\xi_{1}^{2})=p$ and
$\operatorname{var}(\xi_{1}^{2})=\frac{4\beta p}{p+2\beta+2}$, and in
particular we have $\operatorname{var}(\xi_{1}^{2})/p\to 0$ as $p\to\infty$.
Based on Equation (25.14) in [20], the higher moments of $\xi_{1}^{2}$ are

$\displaystyle\mathbb{E}(\xi_{1}^{2k})=\prod_{j=0}^{k-1}\Big{(}p+\textstyle\frac{4j\beta}{p+2\beta+2j}\Big{)}.$

Using the formula relating central moments to ordinary moments

$\mathbb{E}\big{|}\textstyle\frac{\xi_{1}^{2}-p}{\sqrt{p}}\big{|}^{6}\ =\
\frac{1}{p^{3}}\displaystyle\sum_{j=0}^{6}\textstyle\binom{6}{j}(-1)^{6-j}\mathbb{E}(\xi_{1}^{2})^{6-j}\mathbb{E}(\xi_{1}^{2j}),$
(E.3)

it can be checked that that $\mathbb{E}|\frac{\xi_{1}^{2}-p}{\sqrt{p}}|^{6}$
is a rational function of $\beta$ that converges pointwise to 0 as
$p\to\infty$. This implies the $(4+\varepsilon)$-moment condition in (3.2).∎

Beta-Prime distribution. A random variable $W$ is said to follow a Beta-
Prime$(a,b)$ distribution with parameters $a,b>0$ if it can be expressed as
$W=\frac{U}{1-U}$ with $U\sim$ Beta$(a,b)$.

######  Lemma E.3.

If $\tau>0$ and $\xi_{1}^{2}\sim$ Beta-
Prime$\left(\frac{p(1+p+\tau)}{\tau},\frac{1+p+2\tau}{\tau}\right)$, then
$\xi_{1}$ satisfies the conditions in Assumption 1.

###### Proof.

For any positive integer $k<(1+p+2\tau)/\tau$, the random variable
$\xi_{1}^{2}$ satisfies the following moment formula [40, §5.18],

$\displaystyle\mathbb{E}(\xi_{1}^{2k})=\prod_{j=1}^{k}\frac{p(1+p+\tau)+(j-1)\tau}{1+p+(2-j)\tau}.$
(E.4)

Due to $2<(1+p+2\tau)/\tau$, this gives $\mathbb{E}(\xi_{1}^{2})=p$ and
$\mathbb{E}(\xi_{1}^{4})=p(p+\tau)$, which yield
$\operatorname{var}\Big{(}\frac{\xi_{1}^{2}-p}{\sqrt{p}}\Big{)}=\tau$. Also,
using the formula (E.3), it can be checked that
$\mathbb{E}|\frac{\xi_{1}^{2}-p}{\sqrt{p}}|^{6}$ converges pointwise to a
polynomial function of $\tau$ as $p\to\infty$. This implies the
$(4+\varepsilon)$-moment condition in (3.2).∎

Gamma distribution. For $\alpha,\beta>0$, we parameterize the
Gamma$(\alpha,\beta)$ distribution so that its density function is
proportional to $x^{\alpha-1}e^{-\beta x}$ when $x\in(0,\infty)$.

######  Lemma E.4.

If $\tau>0$ and $\xi_{1}^{2}\sim\textup{Gamma}(p/\tau,1/\tau)$, then $\xi_{1}$
satisfies the conditions in Assumption 1.

###### Proof.

The conditions $\mathbb{E}(\xi_{1}^{2})=p$ and
$\operatorname{var}\Big{(}\frac{\xi_{1}^{2}-p}{\sqrt{p}}\Big{)}=\tau$ follow
directly from the stated parameterization. Also, Theorem 2.3 in [8] gives

$\displaystyle\big{(}\mathbb{E}|\xi_{1}^{2}-p|^{6}\big{)}^{1/6}\leq
C(\sqrt{\tau p}+\tau)$ (E.5)

for an absolute constant $C>0$, which implies the $(4+\varepsilon)$-moment
condition in (3.2).∎

Log-Normal distribution. A positive random variable $Y$ is said to follow a
Log-Normal$(\mu,\sigma^{2})$ distribution if $\log(Y)\sim N(\mu,\sigma^{2})$.

######  Lemma E.5.

If $\tau>0$ and $\xi_{1}^{2}\sim$ Log-
Normal$\left(\\!\log(p)-\frac{1}{2}\log\big{(}1+\frac{\tau}{p}\big{)},\log\big{(}1+\frac{\tau}{p}\big{)}\\!\right)$,
then $\xi_{1}$ satisfies the conditions in Assumption 1.

###### Proof.

Equations (14.8a)-(14.8b) in [19] show that if $Y\sim\textup{Log-
Normal}(\mu,\sigma^{2})$, then $\mathbb{E}(Y)=e^{\mu+\sigma^{2}/2}$ and
$\operatorname{var}(Y)=e^{2\mu+\sigma^{2}}\big{(}e^{\sigma^{2}}-1\big{)}$.
Consequently, the stated choice of $\xi_{1}^{2}$ satisfies
$\mathbb{E}(\xi_{1}^{2})=p$ and $\operatorname{var}(\xi_{1}^{2})/p=\tau$.
Equation (14.8c) in [19] shows that the central 6th moment of $Y$ is

$\displaystyle\mathbb{E}|Y-\mathbb{E}(Y)|^{6}$
$\displaystyle=e^{6\mu+3\sigma^{2}}\sum_{j=0}^{6}(-1)^{j}\tbinom{6}{j}e^{\sigma^{2}(6-j)(5-j)/2},$

and a direct calculation reveals that
$\mathbb{E}|\frac{\xi_{1}^{2}-p}{\sqrt{p}}|^{6}$ is a polynomial function of
$\tau$ whose coefficients are bounded as $p\to\infty$. This ensures that
$\xi_{1}$ satisfies the $(4+\varepsilon)$-moment condition in (3.2).