Suppose that the conditions of Theorem 2 hold. Then,

$\max_{1\leq j\leq
k}\big{|}\mathbb{E}[\lambda_{j}(\widehat{\Sigma})]-\lambda_{j}(\Sigma)\big{|}\
\lesssim\ \frac{\log(n)\,\beta_{2q}\operatorname{tr}(\Sigma)}{n}$ (_i_)

and

$\max_{1\leq j\leq
k}\big{|}\operatorname{var}\big{(}\lambda_{j}(\widehat{\Sigma})\big{)}-\lambda_{j}(\Sigma)^{2}\,\Gamma_{jj}/n\big{|}\
\lesssim\
\frac{\log(n)\,\beta_{2q}^{2}\,\lambda_{1}(\Sigma)\operatorname{tr}(\Sigma)}{n^{3/2}}.$
(_ii_)

Also, there is a constant $c>0$ not depending on $n$ such that the events

$\max_{1\leq j\leq
k}\big{|}\mathbb{E}[\lambda_{j}(\widehat{\Sigma}^{\star})|X]-\lambda_{j}(\widehat{\Sigma})\big{|}\
\leq\ \frac{c\,\log(n)\,\beta_{2q}\operatorname{tr}(\Sigma)}{n}$ (_iii_)

and

$\max_{1\leq j\leq
k}\big{|}\operatorname{var}\big{(}\lambda_{j}(\widehat{\Sigma}^{\star})\big{|}X\big{)}-\lambda_{j}(\Sigma)^{2}\,\widehat{\Gamma}_{jj}/n\big{|}\
\leq\
\frac{c\,\log(n)\,\beta_{2q}^{2}\,\lambda_{1}(\Sigma)\operatorname{tr}(\Sigma)}{n^{3/2}}$
(_iv_)

each hold with probability at least $1-\frac{c}{n}$.

Proof. Part (_i_): Recalling the choice $q=5\log(kn)$ from the statement of
Theorem 2, it follows from Lemmas S2.1 and S2.2 that there is a constant $c$
not depending on $n$ such that the event

$\max_{1\leq j\leq
k}\big{|}\lambda_{j}(\widehat{\Sigma})-d_{j}(\tilde{\Sigma}[1,1])\big{|}\
\leq\ \frac{c\,\log(n)\,\beta_{2q}\operatorname{tr}(\Sigma)}{n}$ (S5.64)

holds with probability at least $1-\frac{c}{n^{4}}$. To simplify presentation,
let $\delta$ be a number of the form
$\delta=\frac{c}{n}\log(n)\beta_{2q}\operatorname{tr}(\Sigma)$ and define the
event

$\mathcal{E}_{j}=\Big{\\{}|\lambda_{j}(\widehat{\Sigma})-d_{j}(\tilde{\Sigma}[1,1])|\leq\delta\Big{\\}}$

for each $j=1,\dots,k$. Also, as a temporary short hand, let $U_{j}$ and
$V_{j}$ denote the random variables

$U_{j}=\lambda_{j}(\widehat{\Sigma})-\lambda_{j}(\Sigma)\text{ \ \ \ \ \ \ and
\ \ \ \ \ \ }V_{j}=d_{j}(\tilde{\Sigma}[1,1])-\lambda_{j}(\Sigma).$

Noting that $\mathbb{E}[d_{j}(\tilde{\Sigma}[1,1])]=\lambda_{j}(\Sigma)$, we
have

$\begin{split}\big{|}\mathbb{E}[\lambda_{j}(\widehat{\Sigma})]-\lambda_{j}(\Sigma)\big{|}&\
=\ |\mathbb{E}[U_{j}]-\mathbb{E}[V_{j}]|\\\\[5.69046pt] &\ \leq\
\|U_{j}-V_{j}\|_{2}\\\\[5.69046pt] &\ \leq\ \big{\|}|U_{j}-V_{j}|\cdot
1\\{\mathcal{E}_{j}\\}\big{\|}_{2}+\big{\|}|U_{j}-V_{j}|\cdot
1\\{\mathcal{E}_{j}^{c}\\}\big{\|}_{2}\\\\[5.69046pt] &\ \lesssim\
\delta+\|\lambda_{1}(\tilde{\Sigma})\|_{4}\,\big{(}\mathbb{P}(\mathcal{E}_{j}^{c})\big{)}^{1/4},\end{split}$
(S5.65)

where the fact
$d_{j}(\tilde{\Sigma}[1,1])=d_{j}(\tilde{\Sigma})\leq\lambda_{1}(\tilde{\Sigma})=\lambda_{1}(\widehat{\Sigma})$
has been used in the last step. Due to (S5.64), we have
$\mathbb{P}(\mathcal{E}_{j}^{c})\leq c/n^{4}$, and so it is adequate to derive
a simple upper bound on $\|\lambda_{1}(\tilde{\Sigma})\|_{4}$ using (S2.1) and
condition (7) as follows

$\displaystyle\|\lambda_{1}(\tilde{\Sigma})\|_{4}$ $\displaystyle\ \leq\
\lambda_{1}(\Sigma)+\big{\|}\|\tilde{\Sigma}-\Lambda\|_{\textup{op}}\big{\|}_{q},$
$\displaystyle\ \leq\ c\,\lambda_{1}(\Sigma).$ (S5.66)

Hence, using the last bound in (S5.65) completes the proof of Part (_i_).

Part (_ii_): First note that

$\begin{split}\big{|}\operatorname{var}\big{(}\lambda_{j}(\widehat{\Sigma})\big{)}-\lambda_{j}(\Sigma)^{2}\,\Gamma_{jj}/n\big{|}&\
=\ |\operatorname{var}(U_{j})-\operatorname{var}(V_{j})|\\\\[5.69046pt] &\
\leq\ 2\,(\|U_{j}\|_{2}+\|V_{j}\|_{2})\,\|U_{j}-V_{j}\|_{2}.\end{split}$
(S5.67)

Also note that the intermediate steps in Part (_i_) give

$\|U_{j}-V_{j}\|_{2}\ \lesssim\
\frac{\log(n)\,\beta_{2q}\operatorname{tr}(\Sigma)}{n}.$

From the proof of Lemma S2.3, we have the identity

$V_{j}\ =\ d_{j}(\tilde{\Sigma}[1,1])-\lambda_{j}(\Sigma)\ =\
\frac{1}{n}\sum_{i=1}^{n}(\Lambda_{k}W_{i})_{j},$

which leads to

$\displaystyle\|V_{j}\|_{2}^{2}$ $\displaystyle\ =\
\Big{\|}\frac{1}{n}\sum_{i=1}^{n}(\Lambda_{k}W_{i})_{j}\Big{\|}_{2}^{2}\
\lesssim\ \frac{\beta_{2}^{2}\,\lambda_{1}(\Sigma)^{2}}{n}.$

Also, we can bound $\|U_{j}\|_{2}^{2}$ as follows

$\displaystyle\|U_{j}\|_{2}^{2}$ $\displaystyle\ \leq\
2\|V_{j}\|_{2}^{2}+2\|U_{j}-V_{j}\|_{2}^{2}$ $\displaystyle\ \lesssim\
\frac{\beta_{2}^{2}\,\lambda_{1}(\Sigma)^{2}}{n}+\bigg{(}\frac{\log(n)\beta_{2q}\operatorname{tr}(\Sigma)}{\sqrt{n}}\bigg{)}^{2}\
\frac{1}{n}$ $\displaystyle\ \lesssim\
\frac{\beta_{2}^{2}\,\lambda_{1}(\Sigma)^{2}}{n},$ (S5.68)

where the third step uses the condition (7). Combining the last several bounds
completes the proof of Part (_ii_).

Part (_iii_): The proof is similar to that of Part (_i_). Lemmas S3.1 and S3.2
imply that the event

$\mathbb{P}\bigg{(}\Big{\|}\sqrt{n}\Big{(}{\bm{\lambda}}_{k}(\tilde{\Sigma}^{\star})-{\bm{d}}_{k}(\tilde{\Sigma}^{\star}[1,1])\Big{)}\Big{\|}_{\infty}\,\geq\,\frac{c\,\log(n)\beta_{2q}\operatorname{tr}(\Sigma)}{n^{1/2}}\,\bigg{|}\,X\bigg{)}\
\leq\ \frac{c}{n^{4}}$ (S5.69)

holds with probability at least $1-\frac{c}{n}$. Letting $\delta$ have the
same form as in Part (_i_), define the event

$\mathcal{E}_{j}^{\star}=\Big{\\{}|\lambda_{j}(\widehat{\Sigma}^{\star})-d_{j}(\tilde{\Sigma}^{\star}[1,1])|\leq\delta\Big{\\}}$

for each $j=1,\dots,k$. Also, define

$U_{j}^{\star}=\lambda_{j}(\widehat{\Sigma}^{\star})-d_{j}(\tilde{\Sigma}[1,1])\quad\text{and}\quad
V_{j}^{\star}=d_{j}(\tilde{\Sigma}^{\star}[1,1])-d_{j}(\tilde{\Sigma}[1,1])$

as the bootstrap counterparts of $U_{j}$ and $V_{j}$. To proceed, note that
$\mathbb{E}[d_{j}(\tilde{\Sigma}^{\star}[1,1])|X]=d_{j}(\tilde{\Sigma}[1,1])$,
and so

$\displaystyle\big{|}\mathbb{E}[\lambda_{j}(\widehat{\Sigma}^{\star})|X]-\lambda_{j}(\widehat{\Sigma})\big{|}$
$\displaystyle\ \leq\
|\mathbb{E}[U_{j}^{\star}|X]-\mathbb{E}[V_{j}^{\star}|X]|+|d_{j}(\tilde{\Sigma}[1,1])-\lambda_{j}(\widehat{\Sigma})|.$
(S5.70)

The first term on the right side can be handled similarly to (S5.65),

$\displaystyle|\mathbb{E}[U_{j}^{\star}|X]-\mathbb{E}[V_{j}^{\star}|X]|$
$\displaystyle\ \leq\
\mathbb{E}\big{[}(U_{j}^{\star}-V_{j}^{\star})^{2}1\\{\mathcal{E}_{j}^{\star}\\}\big{|}X\big{]}^{1/2}\
+\
\mathbb{E}\big{[}(U_{j}^{\star}-V_{j}^{\star})^{2}1\\{(\mathcal{E}_{j}^{\star})^{c}\\}\big{|}X\big{]}^{1/2}$
$\displaystyle\ \leq\
c\Big{(}\delta+\big{(}\mathbb{E}\big{[}\lambda_{1}^{4}(\tilde{\Sigma}^{\star})\,\big{|}\,X\big{]}\big{)}^{1/4}\,\big{(}\mathbb{P}(\cup_{j=1}^{k}(\mathcal{E}_{j}^{\star})^{c}|X)\big{)}^{1/4}\,\Big{)},$

where the fact
$d_{j}(\tilde{\Sigma}^{\star}[1,1])=d_{j}(\tilde{\Sigma}^{\star})\leq\lambda_{1}(\tilde{\Sigma}^{\star})$
has been used in the last step. Due to (S5.69), we have
$\mathbb{P}(\cup_{j=1}^{k}(\mathcal{E}_{j}^{\star})^{c}|X)\leq c/n^{4}$ with
probability at least $1-c/n$. Furthermore, it follows from Lemma S6.3 and
(S2.12) that the event

$\begin{split}\big{(}\mathbb{E}\big{[}\lambda_{1}(\tilde{\Sigma}^{\star})^{4}\,\big{|}\,X\big{]}\big{)}^{1/4}&\
\leq\
\lambda_{1}(\Sigma)+\|\tilde{\Sigma}-\Sigma\|_{\textup{op}}+\big{(}\mathbb{E}\big{[}\|\tilde{\Sigma}^{\star}-\tilde{\Sigma}\|_{\textup{op}}^{q}\,\big{|}\,X\big{]}\big{)}^{1/q}\\\\[5.69046pt]
&\ \leq\
\lambda_{1}(\Sigma)+c\,\lambda_{1}(\Sigma)\sqrt{\frac{\log(n)\beta_{q}{\tt{r}}(\Sigma)}{n}}\\\\[5.69046pt]
&\ \leq\ c\,\lambda_{1}(\Sigma)\end{split}$ (S5.71)

holds with probability at least $1-\frac{c}{n}$, where the last line has used
condition (7). Thus, the event

$\max_{1\leq j\leq
k}|\mathbb{E}[U_{j}^{\star}|X]-\mathbb{E}[V_{j}^{\star}|X]|\ \leq\
\frac{c\,\log(n)\,\beta_{2q}\operatorname{tr}(\Sigma)}{n}$

holds with probability at least $1-c/n$. For the second term on the right side
of (S5.70), note that (S5.64) implies that the event

$\max_{1\leq j\leq
k}|d_{j}(\tilde{\Sigma}[1,1])-\lambda_{j}(\widehat{\Sigma})|\ \leq\
\frac{c\,\log(n)\,\beta_{2q}\operatorname{tr}(\Sigma)}{n}$

holds with probability at least $1-c/n$. Applying the last several steps into
(S5.70) completes the proof of Part (_iii_).

Part (_iv_): The proof is essentially analogous to that of Part (_ii_), and so
the details are omitted.∎

######  Lemma S5.4.

Suppose that the conditions of Theorem 2 hold. Then,

$\max_{1\leq j\leq
k}\big{|}\operatorname{var}\big{(}h(\lambda_{j}(\widehat{\Sigma}))\big{)}-h^{\prime}(\lambda_{j}(\Sigma))^{2}\operatorname{var}\big{(}\lambda_{j}(\widehat{\Sigma})\big{)}\big{|}\
\lesssim\
\frac{\log(n)\,\beta_{2q}^{2}\,\lambda_{1}(\Sigma)\,\operatorname{tr}(\Sigma)}{n^{3/2}}.$
(_i_)

In addition, there is a constant $c>0$ not depending on $n$ such that the
events

$\max_{1\leq j\leq
k}\big{|}\operatorname{var}\big{(}h(\lambda_{j}(\widehat{\Sigma}^{\star}))\big{|}X\big{)}-h^{\prime}(\lambda_{j}(\widehat{\Sigma}))^{2}\operatorname{var}\big{(}\lambda_{j}(\widehat{\Sigma}^{\star})\big{|}X\big{)}\big{|}\
\leq\
\frac{c\,\log(n)\,\beta_{2q}^{2}\,\lambda_{1}(\Sigma)\,\operatorname{tr}(\Sigma)}{n^{3/2}}$
(_ii_)

and

$\max_{1\leq j\leq
k}\big{|}\operatorname{var}\big{(}h(\lambda_{j}(\widehat{\Sigma}^{\star}))\big{|}X\big{)}-\operatorname{var}\big{(}h(\lambda_{j}(\widehat{\Sigma}))\big{)}\big{|}\
\leq\
\frac{c\,\log(n)\,\beta_{2q}^{3}\,\lambda_{1}(\Sigma)\,\operatorname{tr}(\Sigma)}{n^{3/2}}$
(_iii_)

each hold with probability at least $1-\frac{c}{n}$.

Proof. Part (_i_): For each $j=1,\dots,k$ define the random variables

$S_{j}=h(\lambda_{j}(\widehat{\Sigma}))-h(\lambda_{j}(\Sigma))\quad\text{and}\quad
T_{j}=h^{\prime}(\lambda_{j}(\Sigma))\,\big{(}\lambda_{j}(\widehat{\Sigma})-\lambda_{j}(\Sigma)\big{)}.$
(S5.72)

In this notation, we have

$\displaystyle\max_{1\leq j\leq
k}\big{|}\operatorname{var}\big{(}h(\lambda_{j}(\widehat{\Sigma}))\big{)}-h^{\prime}(\lambda_{j}(\Sigma))^{2}\operatorname{var}\big{(}\lambda_{j}(\widehat{\Sigma})\big{)}\big{|}$
$\displaystyle\ =\ \max_{1\leq j\leq
k}\big{|}\operatorname{var}(S_{j})-\operatorname{var}(T_{j})\big{|}$
$\displaystyle\ \leq\ \max_{1\leq j\leq
k}2(\|S_{j}\|_{2}+\|T_{j}\|_{2})\|S_{j}-T_{j}\|_{2}.$ (S5.73)

As a way of handling $\|S_{j}-T_{j}\|_{2}$, first note that the argument used
to establish (S5.51) can also be used to show that the event

$\displaystyle\max_{1\leq j\leq k}|S_{j}-T_{j}|$ $\displaystyle\ \leq\
\frac{c\,\log(n)\,\beta_{2q}\operatorname{tr}(\Sigma)}{n}$ (S5.74)

holds with probability at least $1-\frac{c}{n^{4}}$. Using the bound (S5.74)
and an argument analogous to the one leading up to (S5.65), we obtain

$\|S_{j}-T_{j}\|_{2}\ \lesssim\
\frac{\log(n)\beta_{2q}\operatorname{tr}(\Sigma)}{n}+\frac{\|S_{j}\|_{4}+\|T_{j}\|_{4}}{n}.$
(S5.75)

With regard to $\|T_{j}\|_{4}$ we use (S5.66) to obtain the following
conservative but adequate bound,

$\displaystyle\|T_{j}\|_{4}$ $\displaystyle\ \lesssim\
\|\lambda_{j}(\widehat{\Sigma})-\lambda_{j}(\Sigma)\|_{4}$ $\displaystyle\
\leq\ \|\lambda_{j}(\widehat{\Sigma})\|_{4}+\lambda_{j}(\Sigma)$
$\displaystyle\ \lesssim\ \lambda_{1}(\Sigma).$ (S5.76)

To handle $\|S_{j}\|_{4}$, first note that the concavity of $h$ implies that
$S_{j}\leq T_{j}$ almost surely, and so

$0\ \leq\ h(\lambda_{j}(\widehat{\Sigma}))\ \leq\
T_{j}+h(\lambda_{j}(\Sigma)).$

In turn, this yields

$\displaystyle\|S_{j}\|_{4}$
$\displaystyle=\|h(\lambda_{j}(\widehat{\Sigma}))-h(\lambda_{j}(\Sigma))\|_{4}$
$\displaystyle\ \leq\ \|h(\lambda_{j}(\widehat{\Sigma}))\|_{4}\ +\
h(\lambda_{j}(\Sigma))$ $\displaystyle\ \leq\|T_{j}\|_{4}\ +\
2h(\lambda_{j}(\Sigma))$ $\displaystyle\ \lesssim\ \lambda_{1}(\Sigma).$
(S5.77)

So, substituting (S5.1) and (S5.1) into (S5.75) gives

$\|S_{j}-T_{j}\|_{2}\ \lesssim\
\frac{\log(n)\,\beta_{2q}\operatorname{tr}(\Sigma)}{n}.$ (S5.78)

Now we turn to bounding $\|T_{j}\|_{2}$ and $\|S_{j}\|_{2}$ in (S5.1). Using
Lemma S5.3 and the condition (7) we have

$\displaystyle\|T_{j}\|_{2}$ $\displaystyle\ \lesssim\
\sqrt{\operatorname{var}\big{(}\lambda_{j}(\widehat{\Sigma})\big{)}}+\big{|}\mathbb{E}[\lambda_{j}(\widehat{\Sigma})]-\lambda_{j}(\Sigma)\big{|}$
$\displaystyle\ \lesssim\
\frac{\lambda_{j}(\Sigma)\,\beta_{2}}{\sqrt{n}}+\Big{(}\frac{\log(n)\,\beta_{2q}\,{\tt{r}}(\Sigma)}{\sqrt{n}}\Big{)}\frac{\lambda_{1}(\Sigma)}{\sqrt{n}}$
$\displaystyle\ \lesssim\ \frac{\lambda_{1}(\Sigma)\,\beta_{2}}{\sqrt{n}}$

Likewise, we may bound $\|S_{j}\|_{2}$ as

$\displaystyle\|S_{j}\|_{2}$ $\displaystyle\ \leq\
\|T_{j}\|_{2}+\|S_{j}-T_{j}\|_{2}$ $\displaystyle\ \lesssim\
\frac{\lambda_{1}(\Sigma)\,\beta_{2}}{\sqrt{n}}+\Big{(}\frac{\log(n)\,\beta_{2q}\,{\tt{r}}(\Sigma)}{\sqrt{n}}\Big{)}\frac{\lambda_{1}(\Sigma)}{\sqrt{n}}$
$\displaystyle\ \lesssim\ \frac{\lambda_{1}(\Sigma)\,\beta_{2}}{\sqrt{n}},$

which completes the proof of Part (_i_).

Part (_ii_): Define random vectors $S^{\star},T^{\star}\in\mathbb{R}^{k}$ with
entries given by

$S_{j}^{\star}=h(\lambda_{j}(\widehat{\Sigma}^{\star}))-h(\lambda_{j}(\widehat{\Sigma}))\quad\text{and}\quad
T_{j}^{\star}=h^{\prime}(\lambda_{j}(\widehat{\Sigma}))\big{(}\lambda_{j}(\widehat{\Sigma}^{\star})-\lambda_{j}(\widehat{\Sigma})\big{)}.$

As an initial step, it can be shown that there is a constant $c>0$ not
depending on $n$ such that the event

$\mathbb{P}\bigg{(}\|S^{\star}-T^{\star}\|_{\infty}\ \geq\
\frac{c\log(n)\beta_{2q}\operatorname{tr}(\Sigma)}{n}\,\bigg{|}\,X\bigg{)}\
\leq\ \frac{c}{n^{4}}$ (S5.79)

holds with probability at least $1-\frac{c}{n}$. Verifying this is similar to
the proof of Lemma S3.1, and can be handled by showing that

$\|S^{\star}-T^{\star}\|_{\infty}\ \leq\
\frac{c\,\log(n)\,\beta_{2q}\,\operatorname{tr}(\Sigma)}{n}$

holds with probability at least $1-\frac{c}{n^{5}}$. In turn, this can be
shown using the condition (7) and the entrywise Taylor expansion

$\displaystyle
S_{j}^{\star}-T_{j}^{\star}=\big{(}\lambda_{j}(\widehat{\Sigma}^{\star})-\lambda_{j}(\widehat{\Sigma})\big{)}^{2}\,\int_{0}^{1}(1-t)\,h^{\prime\prime}\big{(}\lambda_{j}(\widehat{\Sigma})+t(\lambda_{j}(\widehat{\Sigma}^{\star})-\lambda_{j}(\widehat{\Sigma}))\big{)}\textup{d}t.$

along with (S3.34) and (S3.35).

To proceed with the rest of the proof, we can bound the left side of (_ii_) in
a manner that is similar to (S5.1),

$\small\max_{1\leq j\leq
k}\big{|}\operatorname{var}(S_{j}^{\star}|X)-\operatorname{var}(T_{j}^{\star}|X)\big{|}\leq\
\max_{1\leq j\leq
k}2\big{(}(\mathbb{E}[(S_{j}^{\star})^{2}|X])^{1/2}+(\mathbb{E}[(T_{j}^{\star})^{2}|X])^{1/2}\big{)}\big{(}\mathbb{E}\big{[}(S_{j}^{\star}-T_{j}^{\star})^{2}\,\big{|}\,X\big{]}\big{)}^{1/2}.$
(S5.80)

It follows from an argument analogous to (S5.65) and an application of (S5.79)
that the event

$\small\max_{1\leq j\leq
k}\big{(}\mathbb{E}\big{[}(S_{j}^{\star}-T_{j}^{\star})^{2}\,\big{|}\,X\big{]}\big{)}^{1/2}\
\leq\
\frac{c}{n}\,\bigg{(}\log(n)\,\beta_{2q}\operatorname{tr}(\Sigma)+\max_{1\leq
j\leq
k}\Big{(}(\mathbb{E}[(S_{j}^{\star})^{4}|X])^{1/4}+(\mathbb{E}[(T_{j}^{\star})^{4}|X])^{1/4}\Big{)}\bigg{)}$

holds with probability at least $1-c/n$. Bounds on the conditional fourth
moments can also be derived similarly to the way that the bounds (S5.1) and
(S5.1) were. Namely, by using (S2.12) and (S5.71), it can be shown that the
events

$\max_{1\leq j\leq k}(\mathbb{E}[(T_{j}^{\star})^{4}|X])^{1/4}\ \leq\
c\,\lambda_{1}(\Sigma)$

and

$\max_{1\leq j\leq k}(\mathbb{E}[(S_{j}^{\star})^{4}|X])^{1/4}\ \leq\
c\,\lambda_{1}(\Sigma)$

each hold with probability at least $1-\frac{c}{n}$. Hence, the event

$\max_{1\leq j\leq
k}\big{(}\mathbb{E}\big{[}(S_{j}^{\star}-T_{j}^{\star})^{2}\,\big{|}\,X\big{]}\big{)}^{1/2}\
\leq\ \frac{c\,\log(n)\,\beta_{2q}\operatorname{tr}(\Sigma)}{n}$

holds with probability at least $1-\frac{c}{n}$.

To address the conditional $L_{2}$ norms of $T_{j}^{\star}$ and
$S_{j}^{\star}$ in (S5.80), first note that

$\displaystyle(\mathbb{E}[(T_{j}^{\star})^{2}|X])^{1/2}$ $\displaystyle\ \leq\
h^{\prime}(\lambda_{j}(\Sigma))\big{(}\operatorname{var}(\lambda_{j}(\widehat{\Sigma}^{\star})|X)^{1/2}+\big{|}\mathbb{E}[\lambda_{j}(\widehat{\Sigma}^{\star})|X]-\lambda_{j}(\widehat{\Sigma})\big{|}\big{)},$

In turn, Lemmas S5.3 and S6.4 as well as the condition (7) imply that the
event

$\displaystyle\max_{1\leq j\leq k}(\mathbb{E}[(T_{j}^{\star})^{2}|X])^{1/2}$
$\displaystyle\ \leq\
\frac{c\,\lambda_{j}(\Sigma)\,\beta_{2}}{\sqrt{n}}+\Big{(}\frac{\log(n)\beta_{2q}{\tt{r}}(\Sigma)}{\sqrt{n}}\Big{)}\frac{c\,\lambda_{1}(\Sigma)}{\sqrt{n}}$
$\displaystyle\ \leq\ \frac{c\,\lambda_{1}(\Sigma)\,\beta_{2}}{\sqrt{n}}$

holds with probability at least $1-\frac{c}{n}$. Furthermore, this implies
that the event

$\displaystyle\max_{1\leq j\leq k}(\mathbb{E}[(S_{j}^{\star})^{2}|X])^{1/2}$
$\displaystyle\ \leq\
(\mathbb{E}[(T_{j}^{\star})^{2}|X])^{1/2}+\big{(}\mathbb{E}\big{[}(S_{j}^{\star}-T_{j}^{\star})^{2}\,\big{|}\,X\big{]}\big{)}^{1/2}$
$\displaystyle\ \leq\
\frac{c\,\lambda_{1}(\Sigma)\,\beta_{2}}{\sqrt{n}}+\Big{(}\frac{\log(n)\beta_{2q}{\tt{r}}(\Sigma)}{\sqrt{n}}\Big{)}\frac{c\,\lambda_{1}(\Sigma)}{\sqrt{n}}$
$\displaystyle\ \leq\ \frac{c\,\lambda_{1}(\Sigma)\,\beta_{2}}{\sqrt{n}}$

holds with probability at least $1-\frac{c}{n}$. The proof is completed by
combining the last several steps with (S5.80).

Part (_iii_): Using several applications of the triangle inequality, we have

$\displaystyle\big{|}\operatorname{var}\big{(}h(\lambda_{j}(\widehat{\Sigma}^{\star}))\big{|}X\big{)}-\operatorname{var}\big{(}h(\lambda_{j}(\widehat{\Sigma})\big{)}\big{|}$
$\displaystyle\ \leq\
\big{|}\operatorname{var}\big{(}h(\lambda_{j}(\widehat{\Sigma}^{\star}))\big{|}X\big{)}-h^{\prime}(\lambda_{j}(\widehat{\Sigma}))^{2}\operatorname{var}\big{(}\lambda_{j}(\widehat{\Sigma}^{\star})\big{|}X\big{)}\big{|}$
$\displaystyle\ \ \ \ \ +\
\big{|}\operatorname{var}\big{(}h(\lambda_{j}(\widehat{\Sigma})\big{)}-h^{\prime}(\lambda_{j}(\Sigma))^{2}\operatorname{var}\big{(}\lambda_{j}(\widehat{\Sigma})\big{)}\big{|}$
$\displaystyle\ \ \ \ \ +\
\big{|}h^{\prime}(\lambda_{j}(\widehat{\Sigma}))^{2}\operatorname{var}\big{(}\lambda_{j}(\widehat{\Sigma}^{\star})|X\big{)}-h^{\prime}(\lambda_{j}(\Sigma))^{2}\operatorname{var}\big{(}\lambda_{j}(\widehat{\Sigma})\big{)}\big{|}$

The first and second terms on the right side have been handled by Parts (_ii_)
and (_i_) respectively. To handle the third term, we may use the triangle
inequality to obtain

$\displaystyle\big{|}h^{\prime}(\lambda_{j}(\widehat{\Sigma}))^{2}\operatorname{var}\big{(}\lambda_{j}(\widehat{\Sigma}^{\star})\big{|}X\big{)}-h^{\prime}(\lambda_{j}(\Sigma))^{2}\operatorname{var}\big{(}\lambda_{j}(\widehat{\Sigma})\big{)}\big{|}\
\leq\ \textup{M}_{j}+\textup{M}_{j}^{\prime},$

where the two terms on the right are defined as

$\displaystyle\textup{M}_{j}$ $\displaystyle\ =\
h^{\prime}(\lambda_{j}(\widehat{\Sigma}))^{2}\cdot\big{|}\operatorname{var}\big{(}\lambda_{j}(\widehat{\Sigma}^{\star})\big{|}X\big{)}-\operatorname{var}\big{(}\lambda_{j}(\widehat{\Sigma})\big{)}\big{|}$
$\displaystyle\textup{M}_{j}^{\prime}$ $\displaystyle\ =\
\operatorname{var}\big{(}\lambda_{j}(\widehat{\Sigma})\big{)}\cdot\big{|}h^{\prime}(\lambda_{j}(\widehat{\Sigma}))^{2}-h^{\prime}(\lambda_{j}(\Sigma))^{2}\big{|}$

Using (S2.12) and the mean value theorem, it can be shown that the event

$\max_{1\leq j\leq
k}\big{|}h^{\prime}(\lambda_{j}(\widehat{\Sigma}))^{2}-h^{\prime}(\lambda_{j}(\Sigma))^{2}\big{|}\
\leq\ c\,\sqrt{\frac{\log(n)\,\beta_{q}\,{\tt{r}}(\Sigma)}{n}}$

holds with probability at least $1-\frac{c}{n}$. Also, using Lemma S5.3 and
the condition (7), it follows that

$\displaystyle\max_{1\leq j\leq
k}\operatorname{var}\big{(}\lambda_{j}(\widehat{\Sigma})\big{)}$
$\displaystyle\ \lesssim\ \frac{\lambda_{1}(\Sigma)^{2}\,\beta_{2}^{2}}{n}.$

Hence, the event

$\max_{1\leq j\leq k}\textup{M}_{j}^{\prime}\ \leq\
\frac{c\,\lambda_{1}(\Sigma)^{2}\,\log(n)\,\beta_{2q}^{5/2}\,{\tt{r}}(\Sigma)}{n^{3/2}}$
(S5.81)

holds with probability at least $1-\frac{c}{n}$.

Regarding the quantity $\textup{M}_{j}$, observe that the bound

$\displaystyle\big{|}\operatorname{var}(\lambda_{j}(\widehat{\Sigma}^{\star})|X)-\operatorname{var}(\lambda_{j}(\widehat{\Sigma}))\big{|}$
$\displaystyle\ \leq\
\big{|}\operatorname{var}(\lambda_{j}(\widehat{\Sigma}^{\star})|X)-\lambda_{j}(\Sigma)^{2}\widehat{\Gamma}_{jj}/n\big{|}$
$\displaystyle\ \ \ \ \ +\ \ \
\big{|}\operatorname{var}(\lambda_{j}(\widehat{\Sigma}))-\lambda_{j}(\Sigma)^{2}\Gamma_{jj}/n\big{|}$
$\displaystyle\ \ \ \ \ +\ \ \
\frac{\lambda_{j}(\Sigma)^{2}\,\|\widehat{\Gamma}-\Gamma\|_{\textup{op}}}{n}$

holds almost surely. Consequently, it follows from Lemmas S5.3 and S6.4, that
the event

$\max_{1\leq j\leq k}\textup{M}_{j}\ \leq\
\frac{c\,\log(n)\,\beta_{2q}^{2}\,\lambda_{1}(\Sigma)\operatorname{tr}(\Sigma)}{n^{3/2}}$
(S5.82)

holds with probability at least $1-\frac{c}{n}$. Combining (S5.81) and (S5.82)
completes the proof of Part (_iii_).∎

######  Lemma S5.5.

Suppose that the conditions of Theorem 2 hold. Then,

$\min_{1\leq j\leq
k}\operatorname{var}\big{(}h(\lambda_{j}(\widehat{\Sigma}))\big{)}\ \gtrsim\
\textstyle\frac{\lambda_{1}(\Sigma)^{2}}{n},$

and there is a constant $c\geq 1$ not depending on $n$ such that the event

$\min_{1\leq j\leq
k}\operatorname{var}\big{(}h(\lambda_{j}(\widehat{\Sigma}^{\star}))\big{|}X\big{)}\
\geq\ \textstyle\frac{\lambda_{1}(\Sigma)^{2}}{cn}$

holds with probability at least $1-\frac{c}{n}$.

Proof. For any $j=1,\dots,k$, and any $t>0$, Chebyshev’s inequality and the
triangle inequality give

$\displaystyle\operatorname{var}\big{(}h(\lambda_{j}(\widehat{\Sigma}))\big{)}$
$\displaystyle\ \geq\
t^{2}\,\mathbb{P}\Big{(}\big{|}h(\lambda_{j}(\widehat{\Sigma}))-\mathbb{E}\big{[}h(\lambda_{j}(\widehat{\Sigma}))\big{]}\big{|}\geq
t\Big{)}$ $\displaystyle\ \geq
t^{2}\,\mathbb{P}\Big{(}\sqrt{n}\big{|}h(\lambda_{j}(\widehat{\Sigma}))-h(\lambda_{j}(\Sigma))\big{|}\geq\sqrt{n}\big{(}t+\big{|}\mathbb{E}\big{[}h(\lambda_{j}(\widehat{\Sigma}))\big{]}-h(\lambda_{j}(\Sigma))\big{|}\big{)}\Big{)}.$
(S5.83)

Next, let $S_{j}$ and $T_{j}$ be as defined in (S5.72). Applying Part (_i_) of
Lemma S5.3 along with (S5.78) shows that the bounds

$\displaystyle\big{|}\mathbb{E}\big{[}h(\lambda_{j}(\widehat{\Sigma}))\big{]}-h(\lambda_{j}(\Sigma))\big{|}$
$\displaystyle\ \leq\
|\mathbb{E}[T_{j}]|+|\mathbb{E}[S_{j}]-\mathbb{E}[T_{j}]|$ $\displaystyle\
\lesssim\
|h^{\prime}(\lambda_{j}(\Sigma))|\,\big{|}\mathbb{E}\big{[}\lambda_{j}(\widehat{\Sigma})-\lambda_{j}(\Sigma)\big{]}\big{|}+\|S_{j}-T_{j}\|_{2}$
$\displaystyle\
\lesssim\frac{\log(n)\,\beta_{2q}\operatorname{tr}(\Sigma)}{n},$

hold for every $j=1,\dots,k$. Hence, by taking
$t=\lambda_{1}(\Sigma)/\sqrt{n}$ in (S5.1), and using Lemma S5.1 with the
condition (7), there is a constant $c>0$ not depending on $n$ such that

$\displaystyle\operatorname{var}\big{(}h(\lambda_{j}(\widehat{\Sigma}))\big{)}$
$\displaystyle\ \geq\
t^{2}\,\mathbb{P}\Big{(}\sqrt{n}\big{|}h(\lambda_{j}(\widehat{\Sigma}))-h(\lambda_{j}(\Sigma))\big{|}\geq
c\Big{)}$ $\displaystyle\ \gtrsim\
\frac{\lambda_{1}(\Sigma)^{2}}{n}\bigg{\\{}\mathbb{P}\Big{(}h^{\prime}(\lambda_{j}(\Sigma))\,\zeta_{j}\geq
c\Big{)}-\frac{\log(n)\,\beta_{2q}^{3}{\tt{r}}(\Sigma)}{n^{1/2}}\bigg{\\}}$
$\displaystyle\ \gtrsim\ \frac{\lambda_{1}(\Sigma)^{2}}{n},$

for every $j=1,\dots,k$. Lastly, the proof for the corresponding lower bound
on $\operatorname{var}(h(\lambda_{j}(\widehat{\Sigma}^{\star}))|X)$ is
analogous and so the details are omitted. ∎

##  S6 Proof of Technical Lemmas

######  Lemma S6.1.

Suppose that Assumption 1 holds and let $q\geq 5\log(kn)$. Then,

$\big{\|}\|\tilde{\Sigma}[1,2]\|_{\emph{{op}}}^{2}\big{\|}_{q}\ \lesssim\
\frac{q\,\beta_{2q}\,\lambda_{1}(\Sigma)\operatorname{tr}(\Sigma)}{n^{1-3/(2q)}}.$

Proof. We will need two auxiliary matrices to extract $\tilde{\Sigma}[1,2]$
from $\tilde{\Sigma}$. Define matrices $\Pi_{1}\in\mathbb{R}^{k\times p}$ and
$\Pi_{2}\in\mathbb{R}^{p\times(p-k)}$ according to

$\displaystyle\Pi_{1}\ =\
\begin{bmatrix}I_{k}&\bm{0}\end{bmatrix}\quad\text{and}\quad\Pi_{2}\ =\
\begin{bmatrix}\bm{0}\\\ I_{p-k}\end{bmatrix},$

which allow us to write $\tilde{\Sigma}[1,2]=\Pi_{1}\tilde{\Sigma}\Pi_{2}$.
Next, let $\mathbb{B}^{k}$ and $\mathbb{B}^{p-k}$ denote the unit
$\ell_{2}$-balls in $\mathbb{R}^{k}$ and $\mathbb{R}^{p-k}$, and let
$\mathbb{T}=\mathbb{B}^{k}\times\mathbb{B}^{p-k}$. With this notation in hand,
it follows that

$\displaystyle\|\tilde{\Sigma}[1,2]\|_{\textup{op}}$ $\displaystyle\ =\
\sup_{(u,v)\in\mathbb{T}}u^{\top}\Pi_{1}\tilde{\Sigma}\Pi_{2}v$
$\displaystyle\ =\
\sup_{(u,v)\in\mathbb{T}}(U\Lambda^{1/2}\Pi_{1}^{\top}u)^{\top}\bigg{(}\frac{1}{n}\sum_{i=1}^{n}Z_{i}Z_{i}^{\top}\bigg{)}(U\Lambda^{1/2}\Pi_{2}v).$
(S6.84)

Observe that the vectors $U\Lambda^{1/2}\Pi_{1}^{\top}u$ and
$U\Lambda^{1/2}\Pi_{2}v$ both lie in the ellipsoid
$\mathcal{E}:=U\Lambda^{1/2}(\mathbb{B}^{p})$, and are orthogonal. Therefore,

$\|\tilde{\Sigma}[1,2]\|_{\textup{op}}\ \leq\
\sup_{\begin{subarray}{c}\small(\mathsf{u},\mathsf{v})\in\mathcal{E}\times\mathcal{E}\\\
\langle\mathsf{u},\mathsf{v}\rangle=0\end{subarray}}\frac{1}{n}\sum_{i=1}^{n}\langle\mathsf{u},Z_{i}\rangle\langle\mathsf{v},Z_{i}\rangle.$
(S6.85)

It will be convenient to write the summands in terms of the matrix
$\mathcal{Q}(\mathsf{u},\mathsf{v}):=\frac{1}{2}(\mathsf{u}\mathsf{v}^{\top}+\mathsf{v}\mathsf{u}^{\top})$,
namely

$\langle\mathsf{u},Z_{i}\rangle\langle\mathsf{v},Z_{i}\rangle=Z_{i}^{\top}\mathcal{Q}(\mathsf{u},\mathsf{v})Z_{i}.$

Under this definition, it can be checked that if $\mathsf{u}$ and $\mathsf{v}$
are any pair of orthogonal vectors, then
$\mathbb{E}[Z_{i}^{\top}\mathcal{Q}(\mathsf{u},\mathsf{v})Z_{i}]=0$. Hence, if
we subtract $\mathbb{E}[Z_{i}^{\top}\mathcal{Q}(\mathsf{u},\mathsf{v})Z_{i}]$
from the $i$th term in (S6.85) for each $i=1,\dots,n$, and subsequently drop
the constraint $\langle\mathsf{u},\mathsf{v}\rangle=0$ from the supremum, then
we obtain the bound

$\displaystyle\|\tilde{\Sigma}[1,2]\|_{\textup{op}}$ $\displaystyle\ \leq\
\sup_{\begin{subarray}{c}(\mathsf{u},\mathsf{v})\in\mathcal{E}\times\mathcal{E}\end{subarray}}\bigg{(}\frac{1}{n}\sum_{i=1}^{n}Z_{i}^{\top}\mathcal{Q}(\mathsf{u},\mathsf{v})Z_{i}-\mathbb{E}[Z_{i}^{\top}\mathcal{Q}(\mathsf{u},\mathsf{v})Z_{i}]\bigg{)}.$
(S6.86)

Next, for a given pair of vectors $\mathsf{u},\mathsf{v}\in\mathcal{E}$,
define two associated vectors

$\displaystyle\mathsf{w}$ $\displaystyle\ =\
\textstyle\mathsf{u}/2+\textstyle\mathsf{v}/2,$
$\displaystyle\tilde{\mathsf{w}}$ $\displaystyle\ =\
\textstyle\mathsf{u}/2-\textstyle\mathsf{v}/2,$

which satisfy the algebraic relation

$\mathcal{Q}(\mathsf{u},\mathsf{v})=\mathsf{w}\mathsf{w}^{\top}-\tilde{\mathsf{w}}\tilde{\mathsf{w}}^{\top}.$

To proceed, define the matrix $A\in\mathbb{R}^{p\times 2p}$ as the column
concatenation $A=\begin{pmatrix}U\Lambda^{1/2},U\Lambda^{1/2}\end{pmatrix}$,
and note that both $\mathsf{w}$ and $\tilde{\mathsf{w}}$ lie in the ellipsoid
$\mathcal{E}^{\prime}:=A(\mathbb{B}^{2p})$. As a result, if we write
$\xi_{i}=A^{\top}Z_{i}$, then we have

$\displaystyle\|\tilde{\Sigma}[1,2]\|_{\textup{op}}$ $\displaystyle\ \leq\
\sup_{(\mathsf{w},\tilde{\mathsf{w}})\in\mathcal{E}^{\prime}\times\mathcal{E}^{\prime}}\bigg{(}\bigg{|}\frac{1}{n}\sum_{i=1}^{n}\langle
Z_{i},\mathsf{w}\rangle^{2}-\mathbb{E}[\langle
Z_{i},\mathsf{w}\rangle^{2}]\bigg{|}+\bigg{|}\frac{1}{n}\sum_{i=1}^{n}\langle
Z_{i},\tilde{\mathsf{w}}\rangle^{2}-\mathbb{E}[\langle
Z_{i},\tilde{\mathsf{w}}\rangle^{2}]\bigg{|}\bigg{)}$ $\displaystyle\ \leq\
\sup_{w\in\mathbb{B}^{2p}}\bigg{|}\frac{2}{n}\sum_{i=1}^{n}\langle
A^{\top}Z_{i},w\rangle^{2}-\mathbb{E}[\langle
A^{\top}Z_{i},w\rangle^{2}]\bigg{|}$ $\displaystyle\ =\
\bigg{\|}\frac{2}{n}\sum_{i=1}^{n}\xi_{i}\xi_{i}^{\top}-\mathbb{E}[\xi_{1}\xi_{1}^{\top}]\bigg{\|}_{\textup{op}}.$
(S6.87)

Next, Lemma S7.3 implies that

$\footnotesize\bigg{\|}\Big{\|}\frac{1}{n}\sum_{i=1}^{n}\xi_{i}\xi_{i}^{\top}-\mathbb{E}[\xi_{i}\xi_{i}^{\top}]\Big{\|}_{\textup{op}}\bigg{\|}_{q}\
\leq\
c\,\bigg{(}\sqrt{\frac{q}{n^{1-3/q}}}\,\big{\|}\mathbb{E}[\xi_{1}\xi_{1}^{\top}]\big{\|}_{\textup{op}}^{1/2}\,\big{(}\mathbb{E}\|\xi_{1}\|_{2}^{2q}\big{)}^{\frac{1}{2q}}\bigg{)}\bigvee\bigg{(}\frac{q}{n^{1-3/q}}\,\big{(}\mathbb{E}\|\xi_{1}\|_{2}^{2q}\big{)}^{1/q}\bigg{)}.$
(S6.88)

We can further compute

$\displaystyle\big{\|}\mathbb{E}[\xi_{1}\xi_{1}^{\top}]\big{\|}_{\textup{op}}$
$\displaystyle\ =\
\big{\|}A^{\top}\mathbb{E}[Z_{1}Z_{1}^{\top}]A\big{\|}_{\textup{op}}$
$\displaystyle\ =\ \|A\|_{\textup{op}}^{2}$ $\displaystyle\ \lesssim\
\lambda_{1}(\Sigma),$ (S6.89)

and

$\displaystyle\big{(}\mathbb{E}\|\xi_{1}\|_{2}^{2q}\big{)}^{1/q}$
$\displaystyle\ =\ 2\big{\|}Z_{1}^{\top}U\Lambda U^{\top}Z_{1}\big{\|}_{q}$
$\displaystyle\ =\ 2\bigg{\|}\sum_{j=1}^{p}\lambda_{j}(\Sigma)\langle
u_{j},Z_{1}\rangle^{2}\bigg{\|}_{q}$ $\displaystyle\ \lesssim\
\operatorname{tr}(\Sigma)\beta_{q}.$ (S6.90)

To finish, we use the relation
$\big{\|}\|\tilde{\Sigma}_{12}\|_{\textup{op}}^{2}\big{\|}_{q}=\big{\|}\|\tilde{\Sigma}_{12}\|_{\textup{op}}\big{\|}_{2q}^{2}$.
Specifically, the previous two bounds can be substituted into (S6.87) and
(S6.88) while replacing $q$ with $2q$ and using the condition (7).∎

######  Lemma S6.2.

Suppose that Assumption 1 holds and let $q\geq 5\log(kn)$. Then, there is a
constant $c>0$ not depending on $n$ such that the event

$\Big{(}\mathbb{E}\big{[}\|\tilde{\Sigma}^{\star}[1,2]\|_{\textup{op}}^{2q}\,\big{|}\,X\big{]}\Big{)}^{1/q}\
\leq\
\frac{c\,q\,\beta_{2q}\,\lambda_{1}(\Sigma)\,\operatorname{tr}(\Sigma)}{n^{1-3/(2q)}}$

holds with probability at least $1-ce^{-q}$.

Proof. We will follow the same notation that was used in the proof of Lemma
S6.1. Repeating the argument from that proof up to (S6.86) and using $2q$ in
place of $q$ we have

$\displaystyle\Big{(}\mathbb{E}\big{[}\|\tilde{\Sigma}^{\star}[1,2]\|_{\textup{op}}^{2q}\big{|}X\big{]}\Big{)}^{\frac{1}{2q}}$
$\displaystyle\ \leq\
\Bigg{(}\mathbb{E}\bigg{[}\bigg{(}\sup_{\begin{subarray}{c}(\mathsf{u},\mathsf{v})\in\mathcal{E}\times\mathcal{E}\end{subarray}}\frac{1}{n}\sum_{i=1}^{n}(Z_{i}^{\star})^{\top}\mathcal{Q}(\mathsf{u},\mathsf{v})Z_{i}^{\star}-\mathbb{E}[Z_{i}^{\top}\mathcal{Q}(\mathsf{u},\mathsf{v})Z_{i}]\bigg{)}^{2q}\bigg{|}X\bigg{]}\Bigg{)}^{\frac{1}{2q}}$

Next, observe that

$\mathbb{E}\big{[}(Z_{1}^{\star})^{\top}\mathcal{Q}(\mathsf{u},\mathsf{v})Z_{1}^{\star}\,\big{|}\,X\big{]}\
=\
\frac{1}{n}\sum_{i=1}^{n}Z_{i}^{\top}\mathcal{Q}(\mathsf{u},\mathsf{v})Z_{i}$

and so the triangle inequality for the conditional $L_{2q}$ norm gives

$\displaystyle\Big{(}\mathbb{E}\big{[}\|\tilde{\Sigma}^{\star}[1,2]\|_{\textup{op}}^{2q}\big{|}X\big{]}\Big{)}^{\frac{1}{2q}}$
$\displaystyle\ \leq\
\Bigg{(}\mathbb{E}\bigg{[}\bigg{(}\sup_{\begin{subarray}{c}(\mathsf{u},\mathsf{v})\in\mathcal{E}\times\mathcal{E}\end{subarray}}\frac{1}{n}\sum_{i=1}^{n}(Z_{i}^{\star})^{\top}\mathcal{Q}(\mathsf{u},\mathsf{v})Z_{i}^{\star}-\mathbb{E}\big{[}(Z_{i}^{\star})^{\top}\mathcal{Q}(\mathsf{u},\mathsf{v})Z_{i}^{\star}\,\big{|}\,X\big{]}\bigg{)}^{2q}\bigg{|}X\bigg{]}\Bigg{)}^{\frac{1}{2q}}$
(S6.91) $\displaystyle\ \ \ \ \ \ \ \
+\sup_{\begin{subarray}{c}(\mathsf{u},\mathsf{v})\in\mathcal{E}\times\mathcal{E}\end{subarray}}\bigg{(}\frac{1}{n}\sum_{i=1}^{n}Z_{i}^{\top}\mathcal{Q}(\mathsf{u},\mathsf{v})Z_{i}-\mathbb{E}[Z_{i}^{\top}\mathcal{Q}(\mathsf{u},\mathsf{v})Z_{i}]\bigg{)}.$

With regard to the second term in the last bound, the proof of Lemma S6.1
shows (via Chebyshev’s inequality and condition (7)) that the event

$\sup_{\begin{subarray}{c}(\mathsf{u},\mathsf{v})\in\mathcal{E}\times\mathcal{E}\end{subarray}}\bigg{(}\frac{1}{n}\sum_{i=1}^{n}Z_{i}^{\top}\mathcal{Q}(\mathsf{u},\mathsf{v})Z_{i}-\mathbb{E}[Z_{1}^{\top}\mathcal{Q}(\mathsf{u},\mathsf{v})Z_{1}]\bigg{)}\
\leq\
c\sqrt{\frac{q\beta_{q}\lambda_{1}(\Sigma)\operatorname{tr}(\Sigma)}{n^{1-3/q}}}.$
(S6.92)

holds with probability at least $1-ce^{-q}$, for some constant $c>0$ not
depending on $n$. Therefore, the proof of the current lemma is complete once
we derive a similar bound for the first term on the right side of (S6.91).

Let $\xi_{i}^{\star}=A^{\top}Z_{i}^{\star}$, with $A$ as defined in the proof
of Lemma S6.1. Then, by following the argument leading up to (S6.87) and
applying Lemma S7.3, there is a constant $c>0$ not depending on $n$ such that

$\displaystyle\,\bigg{(}\mathbb{E}\bigg{[}\bigg{(}\sup_{\begin{subarray}{c}(\mathsf{u},\mathsf{v})\in\mathcal{E}\times\mathcal{E}\end{subarray}}$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}(Z_{i}^{\star})^{\top}\mathcal{Q}(\mathsf{u},\mathsf{v})Z_{i}^{\star}-\mathbb{E}\big{[}(Z_{1}^{\star})^{\top}\mathcal{Q}(\mathsf{u},\mathsf{v})Z_{1}^{\star}\,\big{|}\,X\big{]}\bigg{)}^{2q}\,\bigg{|}\,X\bigg{]}\bigg{)}^{\frac{1}{2q}}$
$\displaystyle\ \leq\
c\,\bigg{(}\mathbb{E}\bigg{[}\bigg{\|}\frac{1}{n}\sum_{i=1}^{n}\xi_{i}^{\star}(\xi_{i}^{\star})^{\top}-\mathbb{E}\big{[}\xi_{1}^{\star}(\xi_{1}^{\star})^{\top}\big{|}\,X\big{]}\bigg{\|}_{\textup{op}}^{2q}\,\bigg{|}\,X\bigg{]}\bigg{)}^{\frac{1}{2q}}$
$\displaystyle\ \leq\
c\,\bigg{(}\sqrt{\frac{2q}{n^{1-\frac{3}{2q}}}}\cdot\big{(}\big{\|}\mathbb{E}[\xi_{1}^{\star}(\xi_{1}^{\star})^{\top}\,|\,X]\big{\|}_{\textup{op}}\big{)}^{1/2}\cdot\big{(}\mathbb{E}\big{[}\|\xi_{1}^{\star}\|_{2}^{4q}\,\big{|}\,X\big{]}\big{)}^{\frac{1}{4q}}\bigg{)}\bigvee\bigg{(}\frac{q}{n^{1-\frac{3}{2q}}}\cdot\big{(}\mathbb{E}\big{[}\|\xi_{1}^{\star}\|_{2}^{4q}\,\big{|}\,X\big{]}\big{)}^{\frac{1}{2q}}\bigg{)}.$
(S6.93)

Next, the triangle inequality implies

$\displaystyle\big{\|}\mathbb{E}\big{[}\xi_{1}^{\star}(\xi_{1}^{\star})^{\top}\,\big{|}\,X\big{]}\big{\|}_{\textup{op}}$
$\displaystyle\ =\
\bigg{\|}\frac{1}{n}\sum_{i=1}^{n}(A^{\top}Z_{i})(A^{\top}Z_{i})^{\top}\bigg{\|}_{\textup{op}}$
$\displaystyle\ \leq\
\bigg{\|}\frac{1}{n}\sum_{i=1}^{n}(A^{\top}Z_{i})(A^{\top}Z_{i})^{\top}-A^{\top}A\bigg{\|}_{\textup{op}}+4\lambda_{1}(\Sigma),$

where we have used $\|A^{\top}A\|_{\textup{op}}\leq 4\lambda_{1}(\Sigma)$ and
$\mathbb{E}[(A^{\top}Z_{1})(A^{\top}Z_{1})^{\top}]=A^{\top}A$. By applying
Lemma S7.3 to the first term (along with the condition (7) and Chebyshev’s
inequality), it follows that the event

$\big{\|}\mathbb{E}\big{[}\xi_{1}^{\star}(\xi_{1}^{\star})^{\top}\,\big{|}\,X\big{]}\big{\|}_{\textup{op}}\leq
c\lambda_{1}(\Sigma)$ (S6.94)

holds with probability at least $1-e^{-q}$, for some constant $c>0$ not
depending on $n$.

It remains to develop an upper bound for
$\big{(}\mathbb{E}\big{[}\|\xi_{1}^{\star}\|_{2}^{4q}\,\big{|}\,X\big{]}\big{)}^{\frac{1}{2q}}$.
According to the definition of $\xi_{i}^{\star}$, we have

$\displaystyle\|\xi_{1}^{\star}\|_{2}^{2}$ $\displaystyle\ =\
(Z_{1}^{\star})^{\top}AA^{\top}Z_{1}^{\star}$ $\displaystyle\ =\
2(Z_{1}^{\star})^{\top}\Sigma Z_{1}^{\star},$

and so

$\displaystyle\big{(}\mathbb{E}\big{[}\|\xi_{1}^{\star}\|_{2}^{4q}\,\big{|}\,X\big{]}\big{)}^{\frac{1}{2q}}$
$\displaystyle\ =\ 2\bigg{(}\frac{1}{n}\sum_{i=1}^{n}(Z_{i}^{\top}\Sigma
Z_{i})^{2q}\bigg{)}^{\frac{1}{2q}}$ $\displaystyle\ =\ 2S,$

where the non-negative random variable $S$ is defined by the last line. For
any $t>0$, Chebyshev’s inequality implies

$\begin{split}\mathbb{P}(S\geq e\,t)&\ \leq\
\frac{e^{-2q}\,\mathbb{E}[S^{2q}]}{t^{2q}}\\\\[5.69046pt] &\ =\
\frac{e^{-2q}\,\|Z_{1}^{\top}\Sigma Z_{1}\|_{2q}^{2q}}{t^{2q}}\\\\[5.69046pt]
&\ =\
e^{-2q}\,\Bigg{(}\frac{1}{t}\bigg{\|}\displaystyle\sum_{j=1}^{p}\lambda_{j}(\Sigma)\langle
u_{j},Z_{1}\rangle^{2}\bigg{\|}_{2q}\Bigg{)}^{2q}\\\\[5.69046pt] &\ \leq\
e^{-2q}\,\Big{(}\frac{\beta_{2q}\operatorname{tr}(\Sigma)}{t}\Big{)}^{2q}.\end{split}$
(S6.95)

Combining the last few steps and using the choice
$t=\beta_{2q}\operatorname{tr}(\Sigma)$, it follows that there is a constant
$c>0$ not depending on $n$ such that the event

$\big{(}\mathbb{E}\big{[}\|\xi_{1}^{\star}\|_{2}^{4q}\,\big{|}\,X\big{]}\big{)}^{\frac{1}{2q}}\
\leq\ c\,\beta_{2q}\,\operatorname{tr}(\Sigma)$ (S6.96)

holds with probability at least $1-e^{-2q}$. Hence, we may substitute (S6.94)
and (S6.96) into (S6.93), and use the condition (7) to conclude that the event

$\small\bigg{(}\mathbb{E}\bigg{[}\bigg{(}\sup_{\begin{subarray}{c}(\mathsf{u},\mathsf{v})\in\mathcal{E}\times\mathcal{E}\end{subarray}}\frac{1}{n}\sum_{i=1}^{n}(Z_{i}^{\star})^{\top}\mathcal{Q}(\mathsf{u},\mathsf{v})Z_{i}^{\star}-\mathbb{E}\big{[}(Z_{1}^{\star})^{\top}\mathcal{Q}(\mathsf{u},\mathsf{v})Z_{1}^{\star}\,\big{|}\,X\big{]}\bigg{)}^{2q}\,\bigg{|}\,X\bigg{]}\bigg{)}^{\frac{1}{2q}}\
\leq\
c\,\sqrt{\frac{2q\,\beta_{2q}\,\lambda_{1}(\Sigma)\,\operatorname{tr}(\Sigma)}{n^{1-\frac{3}{2q}}}}$
(S6.97)

holds with probability at least $1-ce^{-q}$. Finally, the proof is completed
by combining (S6.97) and (S6.92) with (S6.91). ∎

######  Lemma S6.3.

Suppose that Assumption 1 holds and let $q\geq 5\log(kn)$. Then, there is a
constant $c>0$ not depending on $n$, such that the event

$\big{(}\mathbb{E}\big{[}\|\tilde{\Sigma}^{\star}-\tilde{\Sigma}\|_{\textup{op}}^{q}\,\big{|}\,X\big{]}\big{)}^{1/q}\
\leq\
c\,\sqrt{\frac{q\,\beta_{q}\,\lambda_{1}(\Sigma)\,\operatorname{tr}(\Sigma)}{n^{1-3/q}}}$

holds with probability at least $1-ce^{-q}$.

Proof. Letting $\xi_{1}^{\star}=\Lambda^{1/2}U^{\top}Z_{1}^{\star}$ and using
Lemma S7.3, we have

$\small\big{(}\mathbb{E}\big{[}\|\tilde{\Sigma}^{\star}-\tilde{\Sigma}\|_{\textup{op}}^{q}\,\big{|}\,X\big{]}\big{)}^{1/q}\
\leq\
c\cdot\bigg{(}\sqrt{\frac{q}{n^{1-3/q}}}\,\|\tilde{\Sigma}\|_{\textup{op}}^{1/2}\,\big{(}\mathbb{E}\big{[}\|\xi_{1}^{\star}\|_{2}^{2q}\,\big{|}\,X\big{]}\big{)}^{\frac{1}{2q}}\bigg{)}\bigvee\bigg{(}\frac{q}{n^{1-3/q}}\,\big{(}\mathbb{E}\big{[}\|\xi_{1}^{\star}\|_{2}^{2q}\,\big{|}\,X\big{]}\big{)}^{\frac{1}{q}}\bigg{)}.$
(S6.98)

Using (S2.1) in the proof of Lemma S2.1, as well as the condition (7), it
follows that the bound

$\displaystyle\|\tilde{\Sigma}\|_{\textup{op}}$ $\displaystyle\ \leq\
\|\tilde{\Sigma}-\Lambda\|_{\textup{op}}+\|\Lambda\|_{\textup{op}}$
$\displaystyle\ \leq\
c\,\lambda_{1}(\Sigma)\sqrt{\frac{q\,\beta_{q}\,{\tt{r}}(\Sigma)}{n^{1-3/q}}}+\lambda_{1}(\Sigma)$
$\displaystyle\ \leq\ c\,\lambda_{1}(\Sigma)$ (S6.99)

holds with probability at least $1-e^{-q}$. Also, recall that the argument
leading up to (S6.96) implies that the event

$\big{(}\mathbb{E}\big{[}\|\xi_{1}^{\star}\|_{2}^{2q}\,\big{|}\,X\big{]}\big{)}^{\frac{1}{2q}}\
\leq\ c\sqrt{\beta_{q}\operatorname{tr}(\Sigma)}$ (S6.100)

holds with probability at least $1-e^{-q}$, for some constant $c>0$ not
depending on $n$. Combining (S6) and (S6.100) with (S6.98) and the condition
(7) completes the proof. ∎

######  Lemma S6.4.

Suppose that Assumption 1 holds and let $q\geq 5\log(kn)$. Then, there is a
constant $c>0$ not depending on $n$ such that the bound

$\|\widehat{\Gamma}-\Gamma\|_{\textup{op}}\ \leq\
\frac{c\,q\,\beta_{2q}^{2}}{n^{1/2}}$ (S6.101)

holds with probability at least $1-e^{-q}$.

Proof. Let $W_{1},\dots,W_{n}$ and $\bar{W}$ be as defined at the beginning of
Appendix S3. Note that $\Gamma=\mathbb{E}[W_{1}W_{1}^{\top}]$, and that the
definition of $\widehat{\Gamma}$ in (S3.24) gives

$\widehat{\Gamma}=\frac{1}{n}\sum_{i=1}^{n}W_{i}W_{i}^{\top}-\bar{W}\bar{W}^{\top}.$

We may apply Lemma S7.3 to obtain

$\displaystyle\Big{(}\mathbb{E}\big{\|}\widehat{\Gamma}-\Gamma\big{\|}_{\textup{op}}^{q}\Big{)}^{1/q}$
$\displaystyle\ \leq\
\bigg{(}\mathbb{E}\bigg{\|}\frac{1}{n}\sum_{i=1}^{n}W_{i}W_{i}^{\top}-\mathbb{E}[W_{1}W_{1}^{\top}]\bigg{\|}_{\textup{op}}^{q}\bigg{)}^{1/q}\
+\ \Big{(}\mathbb{E}\|\bar{W}\bar{W}^{\top}\|_{\textup{op}}^{q}\Big{)}^{1/q}$
$\displaystyle\ \lesssim\
\bigg{(}\sqrt{\frac{q}{n^{1-3/q}}}\,\big{(}\|\mathbb{E}[W_{1}W_{1}^{\top}]\|_{\textup{op}}\big{)}^{1/2}\,\big{(}\mathbb{E}\|W_{1}\|_{2}^{2q}\big{)}^{1/(2q)}\bigg{)}\bigvee\bigg{(}\frac{q}{n^{1-3/q}}\,\big{(}\mathbb{E}\|W_{1}\|_{2}^{2q}\big{)}^{1/q}\bigg{)}\
+\ \Big{(}\mathbb{E}\|\bar{W}\|_{2}^{2q}\Big{)}^{1/q}.$ (S6.102)

It is straightforward to verify that

$\|\mathbb{E}[W_{1}W_{1}^{\top}]\|_{\textup{op}}\ \leq\
\operatorname{tr}(\Gamma)\ \lesssim\ \beta_{2}^{2}.$

Also, we have

$\displaystyle\big{(}\mathbb{E}\|W_{1}\|_{2}^{2q}\big{)}^{1/q}$
$\displaystyle\ =\ \bigg{\|}\sum_{j=1}^{k}\big{(}\langle
u_{j},Z_{1}\rangle^{2}-1\big{)}^{2}\bigg{\|}_{q}$ $\displaystyle\ \leq\
\sum_{j=1}^{k}\big{\|}\langle u_{j},Z_{1}\rangle^{2}-1\big{\|}_{2q}^{2}$
$\displaystyle\ \lesssim\,\beta_{2q}^{2}.$

To handle the term involving $\bar{W}$, observe that the previous bound and
Lemma S7.1 lead to

$\begin{split}\Big{(}\mathbb{E}\|\bar{W}\|_{2}^{2q}\Big{)}^{\frac{1}{2q}}&\
\lesssim\
q\Bigg{(}\big{(}\mathbb{E}\|\bar{W}\|_{2}^{2}\big{)}^{1/2}+\Big{(}\textstyle\sum_{i=1}^{n}\mathbb{E}\|\textstyle\frac{1}{n}W_{i}\|_{2}^{2q}\Big{)}^{\frac{1}{2q}}\Bigg{)}\\\\[5.69046pt]
&\ \lesssim\
q\Bigg{(}\textstyle\frac{\big{(}\mathbb{E}\|W_{1}\|_{2}^{2}\big{)}^{1/2}}{n^{1/2}}+\textstyle\frac{\big{(}\mathbb{E}\|W_{1}\|_{2}^{2q}\big{)}^{\frac{1}{2q}}}{n^{1-\frac{1}{2q}}}\Bigg{)}\\\\[5.69046pt]
&\ \lesssim\
q\Bigg{(}\textstyle\frac{\beta_{2}}{n^{1/2}}+\textstyle\frac{\beta_{2q}}{n^{1-\frac{1}{2q}}}\Bigg{)}\\\\[5.69046pt]
&\ \lesssim\ \frac{q\,\beta_{2q}}{n^{1/2}}\end{split}$ (S6.103)

Using condition (7), we may substitute the last few bounds into (S6) to obtain

$\Big{(}\mathbb{E}\big{\|}\widehat{\Gamma}-\Gamma\big{\|}_{\textup{op}}^{q}\Big{)}^{1/q}\
\lesssim\ \ \frac{q\,\beta_{2q}^{2}}{n^{1/2}}$ (S6.104)

Combining this result with Chebyshev’s inequality completes the proof.

∎

For the next lemma, let $Y_{1}^{\star}$ be as defined at the beginning of the
proof of Lemma S3.3.

######  Lemma S6.5.

Suppose that Assumption 1 holds and let $q\geq 5\log(kn)$. Then, there is a
constant $c>0$ not depending on $n$, such that the event

$\mathbb{E}\Big{[}\big{\|}(\Lambda_{k}\widehat{\Gamma}\Lambda_{k})^{-1/2}Y_{1}^{\star}\big{\|}_{2}^{3}\,\Big{|}\,X\Big{]}\
\leq\ c\,\beta_{3q}^{3}.$

holds with probability at least $1-ce^{-q}$.

Proof. Note that

$\displaystyle\big{\|}(\Lambda_{k}\widehat{\Gamma}\Lambda_{k})^{-1/2}Y_{1}^{\star}\big{\|}_{2}^{2}$
$\displaystyle\ =\
(W_{1}^{\star}-\bar{W})^{\top}\widehat{\Gamma}^{-1}(W_{1}^{\star}-\bar{W})$
$\displaystyle\ \leq\
\frac{(W_{1}^{\star}-\bar{W})^{\top}(W_{1}^{\star}-\bar{W})}{\lambda_{k}(\widehat{\Gamma})}$
$\displaystyle\ \leq\
\frac{2}{\lambda_{k}(\widehat{\Gamma})}\bigg{(}\sum_{j=1}^{k}\big{(}\langle
u_{j},Z_{1}^{\star}\rangle^{2}-1\big{)}^{2}+\bar{W}^{\top}\bar{W}\bigg{)}.$

This implies

$\begin{split}\mathbb{E}\Big{[}\big{\|}(\Lambda_{k}\widehat{\Gamma}\Lambda_{k})^{-1/2}Y_{1}^{\star}\big{\|}_{2}^{3}\Big{|}X\Big{]}&\
\leq\
\frac{c}{\lambda_{k}(\widehat{\Gamma})^{3/2}}\bigg{(}\frac{1}{n}\sum_{i=1}^{n}\Big{(}\sum_{j=1}^{k}\big{(}\langle
u_{j},Z_{i}\rangle^{2}-1\big{)}^{2}\Big{)}^{3/2}+(\bar{W}^{\top}\bar{W})^{3/2}\bigg{)}\\\\[5.69046pt]
&\ =\
\frac{c\,\big{(}T+(\bar{W}^{\top}\bar{W})^{3/2}\big{)}}{\lambda_{k}(\widehat{\Gamma})^{3/2}},\end{split}$
(S6.105)

where the first line has used the convexity of the function $x\mapsto
x^{3/2}$, and the non-negative random variable $T$ is defined in the second
line. By the triangle inequality for the $L^{q}$ norm, we have

$\begin{split}\|T\|_{q}&\ \leq\Bigg{(}\bigg{\|}\sum_{j=1}^{k}\big{(}\langle
u_{j},Z_{1}\rangle^{2}-1\big{)}^{2}\bigg{\|}_{3q/2}^{3q/2}\Bigg{)}^{1/q}\\\\[8.5359pt]
&\ \lesssim\ \max_{1\leq j\leq k}\Big{\|}\langle
u_{j},Z_{1}\rangle^{2}-1\Big{\|}_{3q}^{3}\\\\[5.69046pt] &\ \lesssim\
\beta_{3q}^{3}.\end{split}$ (S6.106)

So, by Chebyshev’s inequality, there is a constant $c>0$ not depending on $n$
such that the event

$T\,\leq\,c\,\beta_{3q}^{3}$

holds with probability at least $1-e^{-q}$. The bound (S6.103) in the proof of
Lemma S6.4 and the condition (7) also imply that there is a constant $c>0$ not
depending on $n$ such that the bound

$(\bar{W}^{\top}\bar{W})^{3/2}\ \leq\
\Big{(}\frac{c\,q\,\beta_{2q}}{n^{1/2}}\Big{)}^{3/2}\ \leq\ c$

holds with probability at least $1-e^{-q}$.

Now we turn to showing that $\lambda_{k}\big{(}\widehat{\Gamma}\big{)}$ is
greater than a positive constant with high probability. Due to Weyl’s
inequality we have $\lambda_{k}(\widehat{\Gamma})\ \geq\
\lambda_{k}(\Gamma)-\|\widehat{\Gamma}-\Gamma\|_{\textup{op}}$. Using Lemma
S6.4, Assumption 1.(c)., and the condition (7), it follows that there is a
constant $c>0$ not depending on $n$ such that the bound

$\lambda_{k}(\widehat{\Gamma})\geq c,$ (S6.107)

holds with probability at least $1-e^{-q}$. Combining the bounds (S6.106) and
(S6.107) with (S6.105) completes the proof. ∎

##  S7 Background Results

######  Lemma S7.1 (Theorem 1 in Talagrand [1989]).

Let $X_{1},\dots,X_{n}$ be independent centered random elements of a Banach
space with norm $\|\cdot\|$. Then, there is an absolute constant $c>0$ such
that the following inequality holds for any $q\geq 1$,

$\displaystyle\bigg{(}\mathbb{E}\Big{\|}\sum_{i=1}^{n}X_{i}\Big{\|}^{q}\bigg{)}^{1/q}\
\leq\
\frac{cq}{1+\log(q)}\bigg{(}\mathbb{E}\Big{\|}\sum_{i=1}^{n}X_{i}\Big{\|}+\Big{(}\mathbb{E}\displaystyle\max_{1\leq
i\leq n}\|X_{i}\|^{q}\Big{)}^{1/q}\bigg{)}.$

######  Lemma S7.2 (Weilandt’s inequality Eaton and Tyler [1991], Wielandt and
Meyer [1967]).

Consider a real symmetric $p\times p$ matrix

$A=\begin{pmatrix}B&C\\\ C^{\top}&D\end{pmatrix},$

where $B$ is $k\times k$ and $D$ is $(p-k)\times(p-k)$. If
$\lambda_{k}(B)>\lambda_{1}(D)$, then

$0\leq\lambda_{j}(A)-\lambda_{j}(B)\leq\frac{\lambda_{1}(CC^{\top})}{\lambda_{j}(B)-\lambda_{1}(D)},\quad
j=1,\dots,k$

and

$0\leq\lambda_{p-k-i}(D)-\lambda_{p-i}(A)\leq\frac{\lambda_{1}(CC^{\top})}{\lambda_{k}(B)-\lambda_{p-k-i}(D)},\quad
i=0,\dots,p-k-1.$

######  Lemma S7.3 (Proposition 1 in Lopes et al. [2022+]).

Let $\xi,\dots,\xi_{n}\in\mathbb{R}^{p}$ be i.i.d. random vectors, let $q\geq
3$, and define the quantity

$r(q)=q\cdot\frac{\big{(}\mathbb{E}\|\xi_{1}\|_{2}^{2q}\big{)}^{1/q}}{\big{\|}\mathbb{E}[\xi_{1}\xi_{1}^{\top}]\big{\|}_{\textup{op}}}.$

Then, there is an absolute constant $c>0$ such that

$\bigg{(}\mathbb{E}\bigg{\|}\frac{1}{n}\sum_{i=1}^{n}\xi_{i}\xi_{i}^{\top}-\mathbb{E}[\xi_{1}\xi_{1}^{\top}]\bigg{\|}_{\textup{op}}^{q}\bigg{)}^{1/q}\
\leq\
c\cdot\big{\|}\mathbb{E}[\xi_{1}\xi_{1}^{\top}]\big{\|}_{\textup{op}}\cdot\bigg{(}\sqrt{\frac{r(q)}{n^{1-3/q}}}\bigvee\frac{r(q)}{n^{1-3/q}}\bigg{)}$

The following anti-concentration lemma originates from Nazarov [2003], and was
further elucidated in [Chernozhukov et al., 2017, Theorem 1].

######  Lemma S7.4 (Nazarov’s inequality).

Let $Y=(Y_{1},\dots,Y_{p})$ be a centered Gaussian random vector in
$\mathbb{R}^{p}$ and suppose that the parameter
$\underline{\sigma}^{2}=\min_{1\leq j\leq p}\mathbb{E}[Y_{j}^{2}]$ is
positive. Then for every $y\in\mathbb{R}^{p}$ and $\delta>0$,

$\mathbb{P}(Y\preceq y+\delta 1_{k})-\mathbb{P}(Y\preceq y)\ \leq\
\frac{\delta}{\underline{\sigma}}\big{(}\sqrt{2\log(p)}+2\big{)}.$

The following lemma is Bentkus’ multivariate Berry-Esseen theorem.

######  Lemma S7.5 (Theorem 1.1 in [Bentkus, 2003]).

Let $V_{1},\dots,V_{n}$ be i.i.d. random vectors $\mathbb{R}^{d}$, with zero
mean and identity covariance matrix. Furthermore, let $\zeta$ be a standard
Gaussian vector in $\mathbb{R}^{d}$, and let $\mathcal{A}$ denote the
collection of all Borel convex subsets of $\mathbb{R}^{d}$. Then, there is an
absolute constant $c>0$ such that

$\sup_{A\in\mathcal{A}}\Big{|}\mathbb{P}\Big{(}\frac{1}{\sqrt{n}}\sum_{i=1}^{n}V_{i}\in
A\Big{)}-\mathbb{P}(\zeta\in A)\Big{|}\ \leq\ \frac{c\cdot
d^{1/4}\cdot\mathbb{E}[\|V_{1}\|_{2}^{3}]}{n^{1/2}}.$

For the statement of Lemma S7.6 below, we need to introduce a bit of notation.
For any $r>0$ and set $A\subset\mathbb{R}^{p}$, define the outer
$r$-neighborhood as $A^{r}=\big{\\{}x\in\mathbb{R}^{p}\,|\,d(x,A)\leq r\\},$
where $d(x,A)=\inf\\{\|x-y\|\,|\,y\in A\\}$, and $\|\cdot\|$ is any norm on
$\mathbb{R}^{p}$. The corresponding inner $r$-neighborhood may be defined as
$A^{-r}=\big{\\{}x\in A\,|\,B(x,r)\subset A\\},$ where
$B(x,r)=\\{y\in\mathbb{R}^{p}|\|x-y\|\leq r\\}$.

######  Lemma S7.6 (Lemma 7.3 in [Lopes, 2020]).

Let $\|\cdot\|$ be any norm on $\mathbb{R}^{p}$, and let
$\zeta,\xi\in\mathbb{R}^{p}$ be any two random vectors. Then, the following
inequality holds for any Borel set $A\subset\mathbb{R}^{p}$, and any $r>0$,

$|\mathbb{P}(\zeta\in A)-\mathbb{P}(\xi\in A)|\ \leq\
\mathbb{P}\big{(}\xi\in(A^{r}\setminus A^{-r})\big{)}\ +\
\mathbb{P}\big{(}\|\zeta-\xi\|\geq r\big{)}.$

The following lemma is a consequence of Pinsker’s inequality and the proof of
Lemma A.7 in the paper Spokoiny and Zhilova [2015].

######  Lemma S7.7 (Spokoiny and Zhilova [2015]).

Let $\zeta$ and $\tilde{\zeta}$ be centered Gaussian vectors in
$\mathbb{R}^{k}$ with respective covariance matrices $C$ and $\tilde{C}$.
Also, suppose that $C$ is invertible, and let
$B=C^{-1/2}\tilde{C}C^{-1/2}-I_{k}$. Then, there is an absolute constant $c>0$
such that

$\sup_{t\in\mathbb{R}^{k}}\Big{|}\mathbb{P}(\zeta\preceq
t)-\mathbb{P}(\tilde{\zeta}\preceq t)\Big{|}\ \leq\
c\sqrt{k}\|B\|_{\textup{op}}.$

The last background lemma follows from the proof of [Lopes et al., 2020b,
Lemma D.3], Lemma S7.6, and Nazarov’s inequality (Lemma S7.4).

######  Lemma S7.8.

Let $U$, $V$, and $R$ be random vectors in $\mathbb{R}^{k}$ that satisfy
$U=V+R$. Also, let $W=(W_{1},\dots,W_{k})$ be a centered Gaussian random
vector in $\mathbb{R}^{k}$ and suppose that the parameter
$\underline{\sigma}^{2}=\min_{1\leq j\leq k}\mathbb{E}[W_{j}^{2}]$ is
positive. Then there is an absolute constant $c>0$, such that the following
bound holds for any $r>0$,

$\sup_{t\in\mathbb{R}^{k}}\big{|}\mathbb{P}(U{\ \preceq t})-\mathbb{P}(W{\
\preceq t})\big{|}\ \leq\ 3\cdot\sup_{t\in\mathbb{R}^{k}}\big{|}\mathbb{P}(V{\
\preceq t})-\mathbb{P}(W{\ \preceq
t})\big{|}+\frac{cr\sqrt{\log(k)}}{\underline{\sigma}}+\mathbb{P}(\|R\|_{\infty}\geq
r).$

##  S8 Additional Numerical Results

#### Nominal value of 95% in model (ii).

The following four figures are presented in the same manner as in the main
text for a 95% nominal value, except that they are based on simulation model
(ii).

\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_tau=00_log_lvl=095_plt_coverage_beta=07.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $\gamma=0.7$ \ \ \ \ }}
\put(-20.0,-5.0){\rotatebox{90.0}{ {\small\ \ \ log transformation \ \ }}}
\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_tau=00_log_lvl=095_plt_coverage_beta=10.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $\gamma=1.0$ \ \ \ \ }} \end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_tau=00_log_lvl=095_plt_coverage_beta=13.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $\gamma=1.3$ \ \ \ \ }} \end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_tau=10_power_10_lvl=095_plt_coverage_beta=07.eps}
\put(-20.0,-1.0){\rotatebox{90.0}{ {\small\ \ \ standardization \ \ \ }}}
\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_tau=10_power_10_lvl=095_plt_coverage_beta=10.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_tau=10_power_10_lvl=095_plt_coverage_beta=13.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_select_power_05_lvl=095_plt_coverage_beta=07.eps}
\put(-21.0,1.0){\rotatebox{90.0}{\ $\sqrt{\ \ }$}}
\put(-20.0,-3.0){\rotatebox{90.0}{ { \ \ \ \ \ \ \small transformation \ \ }
}}
\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_select_power_05_lvl=095_plt_coverage_beta=10.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_select_power_05_lvl=095_plt_coverage_beta=13.eps}
\end{overpic}

Figure 5: (Simultaneous coverage probability versus $n$ in simulation model
(ii) with a polynomial decay profile). The plotting scheme is the same as
described in the caption of Figure 1 in the main text.

\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_tau=00_log_lvl=095_plt_coverage_beta=07.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $\delta=0.7$ \ \ \ \ }}
\put(-20.0,-5.0){\rotatebox{90.0}{ {\small\ \ \ log transformation \ \ }}}
\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_tau=00_log_lvl=095_plt_coverage_beta=08.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $\delta=0.8$ \ \ \ \ }} \end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_tau=00_log_lvl=095_plt_coverage_beta=09.eps}
\par\put(25.0,80.0){ \textul{\ \ \ \ $\delta=0.9$ \ \ \ \ }} \par\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_tau=10_power_10_lvl=095_plt_coverage_beta=07.eps}
\put(-20.0,-1.0){\rotatebox{90.0}{ {\small\ \ \ standardization \ \ \ }}}
\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_tau=10_power_10_lvl=095_plt_coverage_beta=08.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_tau=10_power_10_lvl=095_plt_coverage_beta=09.eps}
\par\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_select_power_05_lvl=095_plt_coverage_beta=07.eps}
\put(-21.0,1.0){\rotatebox{90.0}{\ $\sqrt{\ \ }$}}
\put(-20.0,-3.0){\rotatebox{90.0}{ { \ \ \ \ \ \ \small transformation \ \ }
}}
\par\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_select_power_05_lvl=095_plt_coverage_beta=08.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_select_power_05_lvl=095_plt_coverage_beta=09.eps}
\par\end{overpic}

Figure 6: (Simultaneous coverage probability versus $n$ in simulation model
(ii) with an exponential decay profile). The plotting scheme is the same as
described in the caption of Figure 2 in the main text.

\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_tau=00_log_lvl=095_plt_width_beta=07.eps}
\put(25.0,80.0){ \textul{ \ \ \ \ $\gamma=0.7$ \ \ \ \ }}
\put(-20.0,-5.0){\rotatebox{90.0}{ {\small\ \ \ log transformation \ \ }}}
\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_tau=00_log_lvl=095_plt_width_beta=10.eps}
\put(25.0,80.0){\textul{\ \ \ \ $\gamma=1.0$ \ \ \ \ }} \end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_tau=00_log_lvl=095_plt_width_beta=13.eps}
\put(25.0,80.0){\textul{\ \ \ \ $\gamma=1.3$ \ \ \ \ }} \end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_tau=10_power_10_lvl=095_plt_width_beta=07.eps}
\put(-20.0,-1.0){\rotatebox{90.0}{ {\small\ \ \ standardization \ \ \ }}}
\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_tau=10_power_10_lvl=095_plt_width_beta=10.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_tau=10_power_10_lvl=095_plt_width_beta=13.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_select_power_05_lvl=095_plt_width_beta=07.eps}
\put(-21.0,1.0){\rotatebox{90.0}{\ $\sqrt{\ \ }$}}
\put(-20.0,-3.0){\rotatebox{90.0}{ { \ \ \ \ \ \ \small transformation \ \ }
}}
\par\par\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_select_power_05_lvl=095_plt_width_beta=10.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_select_power_05_lvl=095_plt_width_beta=13.eps}
\end{overpic}

Figure 7: (Average width versus $n$ in simulation model (ii) with a polynomial
decay profile). The plotting scheme is the same as described in the caption of
Figure 3 in the main text.

\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_tau=00_log_lvl=095_plt_width_beta=07.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $\delta=0.7$ \ \ \ \ }}
\put(-20.0,-5.0){\rotatebox{90.0}{ {\small\ \ \ log transformation \ \ }}}
\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_tau=00_log_lvl=095_plt_width_beta=08.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $\delta=0.8$ \ \ \ \ }} \end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_tau=00_log_lvl=095_plt_width_beta=09.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $\delta=0.9$ \ \ \ \ }} \end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_tau=10_power_10_lvl=095_plt_width_beta=07.eps}
\put(-20.0,-1.0){\rotatebox{90.0}{ {\small\ \ \ standardization \ \ \ }}}
\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_tau=10_power_10_lvl=095_plt_width_beta=08.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_tau=10_power_10_lvl=095_plt_width_beta=09.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_select_power_05_lvl=095_plt_width_beta=07.eps}
\put(-21.0,1.0){\rotatebox{90.0}{\ $\sqrt{\ \ }$}}
\put(-20.0,-3.0){\rotatebox{90.0}{ { \ \ \ \ \ \ \small transformation \ \ }
}}
\par\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_select_power_05_lvl=095_plt_width_beta=08.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_select_power_05_lvl=095_plt_width_beta=09.eps}
\par\end{overpic}

Figure 8: (Average width versus $n$ in simulation model (ii) with an
exponential decay profile). The plotting scheme is the same as described in
the caption of Figure 4 in the main text.

#### Nominal value of 90% in models (i) and (ii).

The following eight figures are presented in the same manner as in the main
text, except that they use a nominal value of 90%, and are based on both
models (i) and (ii).

\begin{overpic}[width=125.74689pt]{new_plots/power_elliptical_sqrtexp/power_elliptical_sqrtexp_sim_tau=00_log_lvl=09_plt_coverage_beta=07.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $\gamma=0.7$ \ \ \ \ }}
\put(-20.0,-5.0){\rotatebox{90.0}{ {\small\ \ \ log transformation \ \ }}}
\end{overpic}\begin{overpic}[width=121.41306pt]{new_plots/power_elliptical_sqrtexp/power_elliptical_sqrtexp_sim_tau=00_log_lvl=09_plt_coverage_beta=10.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $\gamma=1.0$ \ \ \ \ }} \end{overpic}

\begin{overpic}[width=121.41306pt]{new_plots/power_elliptical_sqrtexp/power_elliptical_sqrtexp_sim_tau=00_log_lvl=09_plt_coverage_beta=13.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $\gamma=1.3$ \ \ \ \ }} \end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/power_elliptical_sqrtexp/power_elliptical_sqrtexp_sim_tau=10_power_10_lvl=09_plt_coverage_beta=07.eps}
\put(-20.0,-1.0){\rotatebox{90.0}{ {\small\ \ \ standardization \ \ \ }}}
\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/power_elliptical_sqrtexp/power_elliptical_sqrtexp_sim_tau=10_power_10_lvl=09_plt_coverage_beta=10.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/power_elliptical_sqrtexp/power_elliptical_sqrtexp_sim_tau=10_power_10_lvl=09_plt_coverage_beta=13.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/power_elliptical_sqrtexp/power_elliptical_sqrtexp_sim_select_power_05_lvl=09_plt_coverage_beta=07.eps}
\put(-21.0,1.0){\rotatebox{90.0}{\ $\sqrt{\ \ }$}}
\put(-20.0,-3.0){\rotatebox{90.0}{ { \ \ \ \ \ \ \small transformation \ \ }
}}
\par\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/power_elliptical_sqrtexp/power_elliptical_sqrtexp_sim_select_power_05_lvl=09_plt_coverage_beta=10.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/power_elliptical_sqrtexp/power_elliptical_sqrtexp_sim_select_power_05_lvl=09_plt_coverage_beta=13.eps}
\end{overpic}

Figure 9: (Simultaneous coverage probability versus $n$ in simulation model
(i) with a polynomial decay profile). In each panel, the $y$-axis measures
$\mathbb{P}(\cap_{j=1}^{5}\\{\lambda_{j}(\Sigma)\in\widehat{\mathcal{I}}_{j}\\})$
based on a nominal value of 90%, and the $x$-axis measures $n$. The colored
curves correspond to the different values of $p$, indicated in the legend. The
three rows and three columns correspond to labeled choices of transformations
and values of the eigenvalue decay parameter $\gamma$.

\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_tau=00_log_lvl=09_plt_coverage_beta=07.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $\gamma=0.7$ \ \ \ \ }}
\put(-20.0,-5.0){\rotatebox{90.0}{ {\small\ \ \ log transformation \ \ }}}
\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_tau=00_log_lvl=09_plt_coverage_beta=10.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $\gamma=1.0$ \ \ \ \ }} \end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_tau=00_log_lvl=09_plt_coverage_beta=13.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $\gamma=1.3$ \ \ \ \ }} \end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_tau=10_power_10_lvl=09_plt_coverage_beta=07.eps}
\put(-20.0,-1.0){\rotatebox{90.0}{ {\normalsize\ \ \ standardization \ \ \ }}}
\end{overpic}\begin{overpic}[width=121.41306pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_tau=10_power_10_lvl=09_plt_coverage_beta=10.eps}
\end{overpic}

\begin{overpic}[width=121.41306pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_tau=10_power_10_lvl=09_plt_coverage_beta=13.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_select_power_05_lvl=09_plt_coverage_beta=07.eps}
\put(-21.0,1.0){\rotatebox{90.0}{\ $\sqrt{\ \ }$}}
\put(-20.0,-3.0){\rotatebox{90.0}{ { \ \ \ \ \ \ \small transformation \ \ }
}}
\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_select_power_05_lvl=09_plt_coverage_beta=10.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_select_power_05_lvl=09_plt_coverage_beta=13.eps}
\end{overpic}

Figure 10: (Simultaneous coverage probability versus $n$ in simulation model
(ii) with a polynomial decay profile). The plotting scheme is the same as
described in the caption of Figure 9 above.

\begin{overpic}[width=125.74689pt]{new_plots/exp_elliptical_sqrtexp/exp_elliptical_sqrtexp_sim_tau=00_log_lvl=09_plt_coverage_beta=07.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $\delta=0.7$ \ \ \ \ }}
\put(-20.0,-5.0){\rotatebox{90.0}{ {\small\ \ \ log transformation \ \ }}}
\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/exp_elliptical_sqrtexp/exp_elliptical_sqrtexp_sim_tau=00_log_lvl=09_plt_coverage_beta=08.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $\delta=0.8$ \ \ \ \ }} \end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/exp_elliptical_sqrtexp/exp_elliptical_sqrtexp_sim_tau=00_log_lvl=09_plt_coverage_beta=09.eps}
\par\put(25.0,80.0){ \textul{\ \ \ \ $\delta=0.9$ \ \ \ \ }} \par\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/exp_elliptical_sqrtexp/exp_elliptical_sqrtexp_sim_tau=10_power_10_lvl=09_plt_coverage_beta=07.eps}
\put(-20.0,-1.0){\rotatebox{90.0}{ {\small\ \ \ standardization \ \ \ }}}
\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/exp_elliptical_sqrtexp/exp_elliptical_sqrtexp_sim_tau=10_power_10_lvl=09_plt_coverage_beta=08.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/exp_elliptical_sqrtexp/exp_elliptical_sqrtexp_sim_tau=10_power_10_lvl=09_plt_coverage_beta=09.eps}
\par\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/exp_elliptical_sqrtexp/exp_elliptical_sqrtexp_sim_select_power_05_lvl=09_plt_coverage_beta=07.eps}
\put(-21.0,1.0){\rotatebox{90.0}{\ $\sqrt{\ \ }$}}
\put(-20.0,-3.0){\rotatebox{90.0}{ { \ \ \ \ \ \ \small transformation \ \ }
}}
\par\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/exp_elliptical_sqrtexp/exp_elliptical_sqrtexp_sim_select_power_05_lvl=09_plt_coverage_beta=08.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/exp_elliptical_sqrtexp/exp_elliptical_sqrtexp_sim_select_power_05_lvl=09_plt_coverage_beta=09.eps}
\par\end{overpic}

Figure 11: (Simultaneous coverage probability versus $n$ in simulation model
(i) with an exponential decay profile). The plotting scheme is the same as
described in the caption of Figure 9 above, except that the three columns
correspond to values of the eigenvalue decay parameter $\delta$.

\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_tau=00_log_lvl=09_plt_coverage_beta=07.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $\delta=0.7$ \ \ \ \ }}
\put(-20.0,-5.0){\rotatebox{90.0}{ {\small\ \ \ log transformation \ \ }}}
\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_tau=00_log_lvl=09_plt_coverage_beta=08.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $\delta=0.8$ \ \ \ \ }} \end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_tau=00_log_lvl=09_plt_coverage_beta=09.eps}
\par\put(25.0,80.0){ \textul{\ \ \ \ $\delta=0.9$ \ \ \ \ }} \par\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_tau=10_power_10_lvl=09_plt_coverage_beta=07.eps}
\put(-20.0,-1.0){\rotatebox{90.0}{ {\small\ \ \ standardization \ \ \ }}}
\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_tau=10_power_10_lvl=09_plt_coverage_beta=08.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_tau=10_power_10_lvl=09_plt_coverage_beta=09.eps}
\par\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_select_power_05_lvl=09_plt_coverage_beta=07.eps}
\put(-21.0,1.0){\rotatebox{90.0}{\ $\sqrt{\ \ }$}}
\put(-20.0,-3.0){\rotatebox{90.0}{ { \ \ \ \ \ \ \small transformation \ \ }
}}
\par\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_select_power_05_lvl=09_plt_coverage_beta=08.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_select_power_05_lvl=09_plt_coverage_beta=09.eps}
\par\end{overpic}

Figure 12: (Simultaneous coverage probability versus $n$ in simulation model
(ii) with an exponential decay profile). The plotting scheme is the same as
described in the caption of Figure 9 above.

\begin{overpic}[width=125.74689pt]{new_plots/power_elliptical_sqrtexp/power_elliptical_sqrtexp_sim_tau=00_log_lvl=09_plt_width_beta=07.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $\gamma=0.7$ \ \ \ \ }}
\put(-20.0,-5.0){\rotatebox{90.0}{ {\small\ \ \ log transformation \ \ }}}
\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/power_elliptical_sqrtexp/power_elliptical_sqrtexp_sim_tau=00_log_lvl=09_plt_width_beta=10.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $\gamma=1.0$ \ \ \ \ }} \end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/power_elliptical_sqrtexp/power_elliptical_sqrtexp_sim_tau=00_log_lvl=09_plt_width_beta=13.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $\gamma=1.3$ \ \ \ \ }} \end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/power_elliptical_sqrtexp/power_elliptical_sqrtexp_sim_tau=10_power_10_lvl=09_plt_width_beta=07.eps}
\put(-20.0,-1.0){\rotatebox{90.0}{ {\small\ \ \ standardization \ \ \ }}}
\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/power_elliptical_sqrtexp/power_elliptical_sqrtexp_sim_tau=10_power_10_lvl=09_plt_width_beta=10.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/power_elliptical_sqrtexp/power_elliptical_sqrtexp_sim_tau=10_power_10_lvl=09_plt_width_beta=13.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/power_elliptical_sqrtexp/power_elliptical_sqrtexp_sim_select_power_05_lvl=09_plt_width_beta=07.eps}
\put(-21.0,1.0){\rotatebox{90.0}{\ $\sqrt{\ \ }$}}
\put(-20.0,-3.0){\rotatebox{90.0}{ { \ \ \ \ \ \ \small transformation \ \ }
}}
\par\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/power_elliptical_sqrtexp/power_elliptical_sqrtexp_sim_select_power_05_lvl=09_plt_width_beta=10.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/power_elliptical_sqrtexp/power_elliptical_sqrtexp_sim_select_power_05_lvl=09_plt_width_beta=13.eps}
\end{overpic}

Figure 13: (Average width versus $n$ in simulation model (i) with a polynomial
decay profile). In each of the nine panels, the $y$-axis measures the average
width
$\mathbb{E}[|\widehat{\mathcal{I}}_{1}|+\cdots+|\widehat{\mathcal{I}}_{5}|]/5$,
and the $x$-axis measures $n$. The colored curves correspond to the different
values of $p=10$, $50$, $100$, $200$, indicated in the legend. The three rows
and three columns correspond to labeled choices of transformations and values
of the eigenvalue decay parameter $\gamma$.

\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_tau=00_log_lvl=09_plt_width_beta=07.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $\gamma=0.7$ \ \ \ \ }}
\put(-20.0,-5.0){\rotatebox{90.0}{ {\small\ \ \ log transformation \ \ }}}
\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_tau=00_log_lvl=09_plt_width_beta=10.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $\gamma=1.0$ \ \ \ \ }} \end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_tau=00_log_lvl=09_plt_width_beta=13.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $\gamma=1.3$ \ \ \ \ }} \end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_tau=10_power_10_lvl=09_plt_width_beta=07.eps}
\put(-20.0,-1.0){\rotatebox{90.0}{ {\small\ \ \ standardization \ \ \ }}}
\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_tau=10_power_10_lvl=09_plt_width_beta=10.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_tau=10_power_10_lvl=09_plt_width_beta=13.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_select_power_05_lvl=09_plt_width_beta=07.eps}
\put(-21.0,1.0){\rotatebox{90.0}{\ $\sqrt{\ \ }$}}
\put(-20.0,-3.0){\rotatebox{90.0}{ { \ \ \ \ \ \ \small transformation \ \ }
}}
\par\par\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_select_power_05_lvl=09_plt_width_beta=10.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/power_iid_gaussian/power_iid_gaussian_sim_select_power_05_lvl=09_plt_width_beta=13.eps}
\end{overpic}

Figure 14: (Average width versus $n$ in simulation model (ii) with a
polynomial decay profile). The plotting scheme is the same as described in the
caption of Figure 13 above.

\begin{overpic}[width=125.74689pt]{new_plots/exp_elliptical_sqrtexp/exp_elliptical_sqrtexp_sim_tau=00_log_lvl=09_plt_width_beta=07.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $\delta=0.7$ \ \ \ \ }}
\put(-20.0,-5.0){\rotatebox{90.0}{ {\small\ \ \ log transformation \ \ }}}
\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/exp_elliptical_sqrtexp/exp_elliptical_sqrtexp_sim_tau=00_log_lvl=09_plt_width_beta=08.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $\delta=0.8$ \ \ \ \ }} \end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/exp_elliptical_sqrtexp/exp_elliptical_sqrtexp_sim_tau=00_log_lvl=09_plt_width_beta=09.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $\delta=0.9$ \ \ \ \ }} \end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/exp_elliptical_sqrtexp/exp_elliptical_sqrtexp_sim_tau=10_power_10_lvl=09_plt_width_beta=07.eps}
\put(-20.0,-1.0){\rotatebox{90.0}{ {\small\ \ \ standardization \ \ \ }}}
\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/exp_elliptical_sqrtexp/exp_elliptical_sqrtexp_sim_tau=10_power_10_lvl=09_plt_width_beta=08.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/exp_elliptical_sqrtexp/exp_elliptical_sqrtexp_sim_tau=10_power_10_lvl=09_plt_width_beta=09.eps}
\par\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/exp_elliptical_sqrtexp/exp_elliptical_sqrtexp_sim_select_power_05_lvl=09_plt_width_beta=07.eps}
\put(-21.0,1.0){\rotatebox{90.0}{\ $\sqrt{\ \ }$}}
\put(-20.0,-3.0){\rotatebox{90.0}{ { \ \ \ \ \ \ \small transformation \ \ }
}}
\par\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/exp_elliptical_sqrtexp/exp_elliptical_sqrtexp_sim_select_power_05_lvl=09_plt_width_beta=08.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/exp_elliptical_sqrtexp/exp_elliptical_sqrtexp_sim_select_power_05_lvl=09_plt_width_beta=09.eps}
\par\end{overpic}

Figure 15: (Average width versus $n$ in simulation model (i) with an
exponential decay profile). The plotting scheme is the same as described in
the caption of Figure 13 above, except that the three columns correspond to
values of the eigenvalue decay parameter $\delta$.

\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_tau=00_log_lvl=09_plt_width_beta=07.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $\delta=0.7$ \ \ \ \ }}
\put(-20.0,-5.0){\rotatebox{90.0}{ {\small\ \ \ log transformation \ \ }}}
\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_tau=00_log_lvl=09_plt_width_beta=08.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $\delta=0.8$ \ \ \ \ }} \end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_tau=00_log_lvl=09_plt_width_beta=09.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $\delta=0.9$ \ \ \ \ }} \end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_tau=10_power_10_lvl=09_plt_width_beta=07.eps}
\put(-20.0,-1.0){\rotatebox{90.0}{ {\small\ \ \ standardization \ \ \ }}}
\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_tau=10_power_10_lvl=09_plt_width_beta=08.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_tau=10_power_10_lvl=09_plt_width_beta=09.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_select_power_05_lvl=09_plt_width_beta=07.eps}
\put(-21.0,1.0){\rotatebox{90.0}{\ $\sqrt{\ \ }$}}
\put(-20.0,-3.0){\rotatebox{90.0}{ { \ \ \ \ \ \ \small transformation \ \ }
}}
\par\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_select_power_05_lvl=09_plt_width_beta=08.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/exp_iid_gaussian/exp_iid_gaussian_sim_select_power_05_lvl=09_plt_width_beta=09.eps}
\par\end{overpic}

Figure 16: (Average width versus $n$ in simulation model (ii) with an
exponential decay profile). The plotting scheme is the same as described in
the caption of Figure 15 above.

##  S9 Illustration with stock market data

Within the context of finance, PCA is often applied to stock market return
data for the purposes of risk analysis and portfolio selection [Fabozzi et
al., 2007, Ruppert and Matteson, 2015]. Here, we look at several high-
dimensional datasets of stock market returns to illustrate how the bootstrap
can be applied to do inference on parameters of interest in PCA.

Starting from one dataset of S&P 500 returns during the period February 2013
to December 2017 [Nugent, 2017], we isolated four distinct datasets in the
following way. First, we ranked the 500 stocks based on their average monthly
trading volume over the stated time period. Second, we selected four subsets
of the 500 stocks, corresponding to the top 50, 150, 200, and 300 members of
the ranked list. Third, for each stock, we extracted its biweekly log returns
over the time period, resulting in $118$ log return values per stock. (The use
of log returns rather than ordinary returns is a standard practice in finance
[Ruppert and Matteson, 2015].) Altogether, this produced four data matrices of
size $n\times p$ with the same number of rows $n=118$, but differing numbers
of columns $p=50,150,200,300$. In addition to being high-dimensional, these
datasets also conform with our interest in settings that are well suited to
PCA, since the empirical effective rank satisfies
${\tt{r}}(\widehat{\Sigma})\leq 4$ for every dataset.

###  S9.1 Inference on population eigenvalues

When PCA is used to analyze stock market returns, the leading eigenvectors and
eigenvalues of the population covariance matrix $\Sigma$ have special
interpretations. Namely, the eigenvector corresponding to
$\lambda_{1}(\Sigma)$ is often viewed as representing an overall “market
portfolio”, while subsequent eigenvectors represent “principal portfolios”,
which produce returns that are uncorrelated with the overall market return
[Laloux et al., 2000]. Also, the eigenvalues can be interpreted as the
variances (or volatilities) of the returns associated with the principal
portfolios. For this reason, the population eigenvalues are important for risk
assessment, and so it is of interest to quantify the uncertainty in these
unknown parameters.

For each of the four datasets described above, we applied the bootstrap method
with square-root transformation from Section 3.2 to construct simultaneous
confidence intervals for the leading ten eigenvalues
$\lambda_{1}(\Sigma),\dots,\lambda_{10}(\Sigma)$. The bootstrap intervals are
plotted in Figure 17, based on a simultaneous coverage probability of 95%,
with a black dot representing the sample eigenvalue
$\lambda_{j}(\widehat{\Sigma})$ in the $j$th interval for $j=1,\dots 10$. Upon
close inspection, it can be seen that $\lambda_{j}(\widehat{\Sigma})$ tends to
sit slightly above the midpoint of the $j$th interval. This is encouraging,
because it means that the bootstrap intervals are able to counteract the well-
known phenomenon that the leading sample eigenvalues tend to be biased upwards
in high-dimensional settings [Yao et al., 2015, Ch.11]. In addition, as a way
to gain extra empirical support for the bootstrap intervals, we carried out
the following exercise with estimates of
$\lambda_{1}(\Sigma),\dots,\lambda_{10}(\Sigma)$ computed via the method of
QuEST [Ledoit and Wolf, 2015], which is designed for use in high-dimensional
settings, and has been adopted frequently in the literature. Specifically, we
verified that the QuEST estimate of $\lambda_{j}(\Sigma)$ was contained in the
$j$th bootstrap interval for every $j=1,\dots,10$ and $p=50,150,200,300$.
Hence, this makes it more plausible that the bootstrap intervals also contain
the population eigenvalues.

To comment further on the numerical results in Figure 17, first note that in
every panel, the interval for $\lambda_{1}(\Sigma)$ is well separated from the
intervals for $\lambda_{2}(\Sigma),\dots,\lambda_{10}(\Sigma)$, while there is
substantial overlap among the latter intervals. This type of situation occurs
frequently when PCA is applied to stock market return data, and this is
generally interpreted to mean that the overall behavior of the market has a
much more dominant effect on returns than other types of economic factors
Laloux et al. [2000], Ruppert and Matteson [2015]. A second observation is
that the bootstrap intervals can provide some additional insight into the
relationship between $\lambda_{2}(\Sigma)$ and $\lambda_{3}(\Sigma)$. On one
hand, a user who only looks at the sample eigenvalues
$\lambda_{2}(\widehat{\Sigma})$ and $\lambda_{3}(\widehat{\Sigma})$ might be
tempted to conclude that there is a clear difference between the population
eigenvalues $\lambda_{2}(\Sigma)$ and $\lambda_{3}(\Sigma)$. On the other
hand, a user who looks at the overlap of the second and third intervals would
have more information to see that the difference between $\lambda_{2}(\Sigma)$
and $\lambda_{3}(\Sigma)$ might actually be negligible. Lastly, one more
aspect of Figure 17 to mention is that the relative positions of the ten
intervals stay approximately the same for each of the four dimensions
$p=50,150,200,300$. Given that the empirical effective rank satisfies
${\tt{r}}(\widehat{\Sigma})\leq 4\ll p$ for every dataset, this makes sense
from the standpoint of our theoretical results, which indicate that the
bootstrap should be relatively insensitive to the ambient dimension compared
to the effective rank.

\begin{overpic}[width=173.44534pt]{new_plots/data_analysis/p=50/power_05_95_eig.eps}
\par\put(31.0,80.0){ \textul{\ \ \ \ $p=50$ \ \ \ \ }}
\end{overpic}\begin{overpic}[width=173.44534pt]{new_plots/data_analysis/p=150/power_05_95_eig.eps}
\put(31.0,80.0){ \textul{\ \ \ \ $p=150$ \ \ \ \ }} \end{overpic}

\begin{overpic}[width=173.44534pt]{new_plots/data_analysis/p=200/power_05_95_eig.eps}
\put(31.0,80.0){ \textul{\ \ \ \ $p=200$ \ \ \ \ }}
\end{overpic}\begin{overpic}[width=173.44534pt]{new_plots/data_analysis/p=300/power_05_95_eig.eps}
\put(31.0,80.0){ \textul{\ \ \ \ $p=300$ \ \ \ \ }} \end{overpic}

Figure 17: (Simultaneous confidence intervals for
$\lambda_{1}(\Sigma),\dots,\lambda_{10}(\Sigma)$.) In each panel, the $y$-axis
corresponds to the magnitude of eigenvalues, and the $x$-axis corresponds to
the index $j=1,\dots,10$. The intervals are based on a simultaneous coverage
probability of 95%, and the black dots represent the sample eigenvalues
$\lambda_{j}(\widehat{\Sigma})$ for each $j$.

###  S9.2 Inference on proportions of explained variance

The proportions of explained variance, denoted
$\pi_{j}(\Sigma)=\sum_{i=1}^{j}\lambda_{i}(\Sigma)/\operatorname{tr}(\Sigma)$
for $j=1,\dots,p$, often play a decisive role in applications of PCA, since
they form of the basis of standard decision rules for selecting an appropriate
number of components. To complement our previous example dealing with the
population eigenvalues $\lambda_{1}(\Sigma),\dots,\lambda_{10}(\Sigma)$, this
subsection looks instead at inference with simultaneous confidence intervals
for the parameters $\pi_{1}(\Sigma),\dots,\pi_{10}(\Sigma)$. As before, the
bootstrap method with square-root transformation from Section 3.2 was used to
construct the intervals based on a simultaneous coverage probability of 95%.
The results are given in Figure 18, with black dots showing the locations of
the empirical proportions $\pi_{j}(\widehat{\Sigma})$ within the $j$th
interval.

Due to the fact that the proportions $\pi_{1}(\Sigma),\dots,\pi_{p}(\Sigma)$
are unknown, one of the most widely used rules for selecting the number of
components is to choose the smallest number $k$ for which
$\pi_{k}(\widehat{\Sigma})$ exceeds a given threshold. Although this rule may
be appropriate when dealing with low-dimensional data, it is known in the
literature that this rule can be unreliable in high-dimensional settings,
because it tends to select too few components [Ledoit and Wolf, 2015]. One way
of avoiding this pitfall is to consider the following simple modification,
based on simultaneous bootstrap confidence intervals for the proportions: If
$\widehat{l}_{j}$ denotes the lower endpoint of the $j$th interval, then the
bootstrap-based rule selects the smallest number $k$ for which
$\widehat{l}_{k}$ exceeds the threshold. Consequently, if the intervals
perform properly with a simultaneous coverage probability of 95%, then the
bootstrap-based rule will select a sufficient number of components with at
least 95% probability.

\begin{overpic}[width=173.44534pt]{new_plots/data_analysis/p=50/power_05_95_top_eig.eps}
\par\put(31.0,80.0){ \textul{\ \ \ \ $p=50$ \ \ \ \ }}
\end{overpic}\begin{overpic}[width=173.44534pt]{new_plots/data_analysis/p=150/power_05_95_top_eig.eps}
\put(31.0,80.0){ \textul{\ \ \ \ $p=150$ \ \ \ \ }} \end{overpic}

\begin{overpic}[width=173.44534pt]{new_plots/data_analysis/p=200/power_05_95_top_eig.eps}
\put(31.0,80.0){ \textul{\ \ \ \ $p=200$ \ \ \ \ }}
\end{overpic}\begin{overpic}[width=173.44534pt]{new_plots/data_analysis/p=300/power_05_95_top_eig.eps}
\put(31.0,80.0){ \textul{\ \ \ \ $p=300$ \ \ \ \ }} \end{overpic}

Figure 18: (Simultaneous confidence intervals for the proportions of explained
variance $\pi_{1}(\Sigma),\dots,\pi_{10}(\Sigma)$.) In each panel, the
$y$-axis corresponds to the magnitude of proportions, and the $x$-axis
corresponds to the index $j=1,\dots,10$. The intervals are based on a
simultaneous coverage probability of 95%, and the black dots represent the
empirical proportions $\pi_{j}(\widehat{\Sigma})$ for each $j$.

To illustrate the difference between the two rules, a red line corresponding
to a particular threshold of 0.4 has been drawn in each panel of Figure 18.
(The value of 0.4 has no special importance, and is used only for ease of
presentation.) In each of the four cases $p=50,150,200,300$, the original rule
selects the respective values $k=2,3,3,4$. By contrast, the bootstrap-based
rule selects the respective values $k=4,6,7,8$, and hence, it clearly
counteracts the problem of selecting too few components. Moreover, from a
financial standpoint, it makes sense that extra components are needed as $p$
increases, because as a wider variety of stocks are included in the data,
there is greater opportunity for the returns to be influenced by economic
factors that are not captured by the previously leading components.

##  S10 Computational cost

The cost to compute a single bootstrap sample of
$\bm{\lambda}_{k}(\widehat{\Sigma}^{\star})-\bm{\lambda}_{k}(\widehat{\Sigma})$
can be broken into two steps. The first step is to sample $n$ points with
replacement from the original set of $n$ observations, which has a cost of
$O(n\log(n))$. The second step consists of computing the largest $k$
eigenvalues of $\widehat{\Sigma}^{\star}$. This can be done by computing the
largest $k$ singular values of the $n\times p$ matrix of resampled
observations, which has a cost of $\mathcal{O}(npk)$ [Halko et al., 2011].
(Note that the largest $k$ eigenvalues of $\widehat{\Sigma}$ only need to be
computed once, before any resampling is done.) So, the cost to compute $B$
bootstrap samples of
$\bm{\lambda}_{k}(\widehat{\Sigma}^{\star})-\bm{\lambda}_{k}(\widehat{\Sigma})$
on a single processor is $\mathcal{O}(B(n\log(n)+npk))$. However, it is common
to compute bootstrap samples in a parallel manner, across say $m$ processors,
and in this case the cost per processor becomes
$\mathcal{O}(\frac{B}{m}(n\log(n)+npk)$. In order to simplify this expression,
it is natural to consider a scenario where $B/m=\mathcal{O}(1)$ and
$\log(n)=\mathcal{O}(p)$, which leads to a cost per processor that is
$\mathcal{O}(npk)$.

###  S10.1 Empirical computational cost

Table 1 displays the time, in seconds, to compute a single bootstrap sample of
$\bm{\lambda}_{k}(\widehat{\Sigma}^{\star})-\bm{\lambda}_{k}(\widehat{\Sigma})$
for different choices of $n$ and $p$. Each entry reflects an average over
$1000$ trials with data generated under simulation model (i). The computations
were done on a single Intel Xeon E5-2699v3 processor. Notably, even when
$(n,p)=(500,200)$, the computing time for each bootstrap sample is on the
order of just $10^{-3}$ seconds. Hence, it is possible to generate hundreds
bootstrap samples within about 1 second.

| $p$
---|---
$n$ | 10 | 50 | 100 | 200
10 | 1.9e-04 | 4.2e-04 | 8.7e-04 | 2.6e-03
50 | 2.0e-04 | 4.4e-04 | 9.0e-04 | 2.8e-03
100 | 2.1e-04 | 4.6e-04 | 9.6e-04 | 2.9e-03
200 | 2.2e-04 | 5.1e-04 | 1.1e-03 | 3.0e-03
500 | 2.4e-04 | 6.5e-04 | 1.5e-03 | 4.0e-03
Table 1: (Average time, in seconds, to compute one bootstrap sample from
simulated data.) The rows correspond to the values of $n=10,50,100,200,500$,
and the columns correspond to the values of $p=10,50,100,200$.

Table 2 displays analogous computing times for
$\bm{\lambda}_{k}(\widehat{\Sigma}^{\star})-\bm{\lambda}_{k}(\widehat{\Sigma})$,
in seconds, based on the stock market data. Here, the computations were done
on a single 3.5 GHz Dual-Core Intel Core i7 processor. The results show that
there is little difference in computing time compared to the setting of
synthetic data in Table 1.

| $p$
---|---
$n$ | 50 | 150 | 200 | 300
118 | 3.8e-04 | 1.4e-03 | 2.0e-03 | 4.0e-03
Table 2: (Average time, in seconds, to compute one bootstrap sample from the
stock market data.) The single row corresponds to the value of $n=118$, and
the columns correspond to the values of $p=50,150,200,300$.

##  S11 Additional discussion on Assumption 1(b)

Recall that Assumption 1(b) requires the condition $\min_{1\leq j\leq
k}(\lambda_{j}(\Sigma)-\lambda_{j+1}(\Sigma))\gtrsim\lambda_{1}(\Sigma)$,
which ensures that there are gaps between the leading eigenvalues
$\lambda_{1}(\Sigma),\dots,\lambda_{k+1}(\Sigma)$. To inspect whether or not
this type of condition holds in practice, the paper [Hall et al., 2009]
proposes a diagnostic method that constructs a preliminary set of
_conservative_ simultaneous confidence intervals for
$\lambda_{1}(\Sigma),\dots,\lambda_{k+1}(\Sigma)$. In exchange for their
conservatism, these intervals have the property that they are not sensitive to
the existence of gaps. (See also [Lopes et al., 2022+] for theoretical
analysis related to this technique.) If none of the intervals overlap, then
the user may conclude that the eigenvalues are adequately separated. However,
if some of the intervals do overlap, then the paper [Hall et al., 2009]
recommends that an adjusted form of bootstrapping be used to approximate the
distribution of the statistic
$\bm{\lambda}_{k}(\widehat{\Sigma})-\bm{\lambda}_{k}(\Sigma)$.

###  S11.1 Sensitivity analysis

To study the sensitivity of the bootstrap to Assumption 1(b), we now discuss
some numerical experiments with varying gaps between the leading population
eigenvalues. The first three eigenvalues were specified as
$(\lambda_{1}(\Sigma),\lambda_{2}(\Sigma),\lambda_{3}(\Sigma))=(1+g,1,1-g)$
for a gap parameter $g\in\\{0,0.1,0.2\\}$, and the remaining eigenvalues were
chosen to follow the polynomial decay profile $\lambda_{j}(\Sigma)=j^{-1}$ for
$j\geq 4$. Next, we applied the bootstrap (as in Section 3) to construct
simultaneous confidence intervals for
$(\lambda_{1}(\Sigma),\dots,\lambda_{5}(\Sigma))$ based on a nominal level of
95%. Figure 19 contains a grid of plots in which the columns correspond to
increasing values of the gap parameter $g$, and the rows correspond to the
three transformation rules considered in Section 3. These plots are based on
data generated from simulation model (i), and an analogous set of plots based
on simulation model (ii) are given in Figure 20.

As expected, Figures 19 and 20 show that the coverage accuracy of the
bootstrap confidence intervals improves as the gap parameter increases. In the
case when $g=0$, the coverage accuracy is poor for all three transformation
rules, even at large sample sizes. Next, when $g=0.1$, the coverage is mostly
accurate for sample sizes $n\geq 400$, but at smaller sample sizes the
coverage generally falls below the desired level. Lastly, when $g=0.2$, the
coverage becomes more accurate when $n<400$, especially when the square-root
transformation is used.

\begin{overpic}[width=125.74689pt]{new_plots/extra_sim/elliptical_sqrtexp/power_tied_elliptical_sqrtexp_sim_tau=00_log_lvl=095_plt_coverage_beta=10.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $g=0$ \ \ \ \ }}
\put(-20.0,-5.0){\rotatebox{90.0}{ {\small\ \ \ log transformation \ \ }}}
\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/extra_sim/elliptical_sqrtexp/power_close=01_elliptical_sqrtexp_sim_tau=00_log_lvl=095_plt_coverage_beta=10.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $g=0.1$ \ \ \ \ }} \end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/extra_sim/elliptical_sqrtexp/power_close=02_elliptical_sqrtexp_sim_tau=00_log_lvl=095_plt_coverage_beta=10.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $g=0.2$ \ \ \ \ }} \end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/extra_sim/elliptical_sqrtexp/power_tied_elliptical_sqrtexp_sim_tau=10_power_10_lvl=095_plt_coverage_beta=10.eps}
\put(-20.0,-1.0){\rotatebox{90.0}{ {\small\ \ \ standardization \ \ \ } }}
\par\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/extra_sim/elliptical_sqrtexp/power_close=01_elliptical_sqrtexp_sim_tau=10_power_10_lvl=095_plt_coverage_beta=10.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/extra_sim/elliptical_sqrtexp/power_close=02_elliptical_sqrtexp_sim_tau=10_power_10_lvl=095_plt_coverage_beta=10.eps}
\par\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/extra_sim/elliptical_sqrtexp/power_tied_elliptical_sqrtexp_sim_select_power_05_lvl=095_plt_coverage_beta=10.eps}
\put(-21.0,1.0){\rotatebox{90.0}{\ $\sqrt{\ \ }$}}
\put(-20.0,-3.0){\rotatebox{90.0}{ { \ \ \ \ \ \ \small transformation \ \ }}}
\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/extra_sim/elliptical_sqrtexp/power_close=01_elliptical_sqrtexp_sim_select_power_05_lvl=095_plt_coverage_beta=10.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/extra_sim/elliptical_sqrtexp/power_close=02_elliptical_sqrtexp_sim_select_power_05_lvl=095_plt_coverage_beta=10.eps}
\end{overpic}

Figure 19: (Simultaneous coverage probability versus $n$ in simulation model
(i) with the decay profile:
$(\lambda_{1})(\Sigma),\lambda_{2}(\Sigma),\lambda_{3}(\Sigma))=(1+g,1,1-g)$
and $\lambda_{j}(\Sigma)=j^{-1}$ for $j\geq 4$). In each panel, the $y$-axis
measures
$\mathbb{P}(\cap_{j=1}^{5}\\{\lambda_{j}(\Sigma)\in\widehat{\mathcal{I}}_{j}\\})$
based on a nominal level of 95%, and the $x$-axis measures $n$. The colored
curves correspond to the different values of $p$, indicated in the legend.

\begin{overpic}[width=125.74689pt]{new_plots/extra_sim/iid_gaussian/power_tied_iid_gaussian_sim_tau=00_log_lvl=095_plt_coverage_beta=10.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $g=0$ \ \ \ \ }}
\put(-20.0,-5.0){\rotatebox{90.0}{ {\small\ \ \ log transformation \ \ }}}
\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/extra_sim/iid_gaussian/power_close=01_iid_gaussian_sim_tau=00_log_lvl=095_plt_coverage_beta=10.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $g=0.1$ \ \ \ \ }} \end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/extra_sim/iid_gaussian/power_close=02_iid_gaussian_sim_tau=00_log_lvl=095_plt_coverage_beta=10.eps}
\put(25.0,80.0){ \textul{\ \ \ \ $g=0.2$ \ \ \ \ }} \end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/extra_sim/iid_gaussian/power_tied_iid_gaussian_sim_tau=10_power_10_lvl=095_plt_coverage_beta=10.eps}
\put(-20.0,-1.0){\rotatebox{90.0}{ {\small\ \ \ standardization \ \ \ } }}
\par\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/extra_sim/iid_gaussian/power_close=01_iid_gaussian_sim_tau=10_power_10_lvl=095_plt_coverage_beta=10.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/extra_sim/iid_gaussian/power_close=02_iid_gaussian_sim_tau=10_power_10_lvl=095_plt_coverage_beta=10.eps}
\par\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/extra_sim/iid_gaussian/power_tied_iid_gaussian_sim_select_power_05_lvl=095_plt_coverage_beta=10.eps}
\put(-21.0,1.0){\rotatebox{90.0}{\ $\sqrt{\ \ }$}}
\put(-20.0,-3.0){\rotatebox{90.0}{ { \ \ \ \ \ \ \small transformation \ \ }}}
\end{overpic}\begin{overpic}[width=125.74689pt]{new_plots/extra_sim/iid_gaussian/power_close=01_iid_gaussian_sim_select_power_05_lvl=095_plt_coverage_beta=10.eps}
\end{overpic}

\begin{overpic}[width=125.74689pt]{new_plots/extra_sim/iid_gaussian/power_close=02_iid_gaussian_sim_select_power_05_lvl=095_plt_coverage_beta=10.eps}
\end{overpic}

Figure 20: (Simultaneous coverage probability versus $n$ in simulation model
(ii) with the decay profile:
$(\lambda_{1})(\Sigma),\lambda_{2}(\Sigma),\lambda_{3}(\Sigma))=(1+g,1,1-g)$
and $\lambda_{j}(\Sigma)=j^{-1}$ for $j\geq 4$). The plotting scheme is the
same as described in the caption of Figure 19.