# Distributional Theory and Statistical Inference for Linear Functions of
Eigenvectors with Small Eigengaps

Joshua Agterberg Department of Electrical and Systems Engineering and
Department of Statistics and Data Science, University of Pennsylvania

###### Abstract

Spectral methods have myriad applications in high-dimensional statistics and
data science, and while previous works have primarily focused on $\ell_{2}$ or
$\ell_{2,\infty}$ eigenvector and singular vector perturbation theory, in many
settings these analyses fall short of providing the fine-grained guarantees
required for various inferential tasks. In this paper we study statistical
inference for linear functions of eigenvectors and principal components with a
particular emphasis on the setting where gaps between eigenvalues may be
extremely small relative to the corresponding spiked eigenvalue, a regime
which has been oft-neglected in the literature. It has been previously
established that linear functions of eigenvectors and principal components
incur a non-negligible bias, so in this work we provide Berry-Esseen bounds
for empirical linear forms and their debiased counterparts respectively in the
matrix denoising model and the spiked principal component analysis model, both
under Gaussian noise. Next, we propose data-driven estimators for the
appropriate bias and variance quantities resulting in approximately valid
confidence intervals, and we demonstrate our theoretical results through
numerical simulations. We further apply our results to obtain distributional
theory and confidence intervals for eigenvector entries, for which debiasing
is not necessary. Crucially, our proposed confidence intervals and bias-
correction procedures can all be computed directly from data _without_ sample-
splitting and are asymptotically valid under minimal assumptions on the
eigengap and signal strength. Furthermore, our Berry-Esseen bounds clearly
reflect the effects of both signal strength and eigenvalue closeness on the
estimation and inference tasks.

###### Contents

  1. 1 Introduction
    1. 1.1 Two Canonical Settings
    2. 1.2 Contributions
    3. 1.3 Organization
    4. 1.4 Notation
  2. 2 Statistical Models and Key Quantities
    1. 2.1 Matrix Denoising
    2. 2.2 Principal Component Analysis
  3. 3 Distributional Theory
    1. 3.1 Matrix Denoising
    2. 3.2 Principal Component Analysis
  4. 4 Statistical Inference
    1. 4.1 Matrix Denoising
    2. 4.2 Principal Component Analysis
  5. 5 Previous Work
  6. 6 Numerical Simulations
    1. 6.1 Matrix Denoising
    2. 6.2 PCA
  7. 7 Discussion
  8. 8 Technical Tools and Proofs of Main Results
    1. 8.1 Matrix Analytic Tools
    2. 8.2 Analysis for Matrix Denoising
      1. 8.2.1 Step 1: Characterizing the Bias
      2. 8.2.2 Step 2: Eigenvalue Concentration
      3. 8.2.3 Step 3: Finding the Leading-Order Expansion
      4. 8.2.4 Step 4: Completing the Proofs of Theorems 3.1 and 3.2
      5. 8.2.5 Step 5: Variance Approximation and Proof of Theorem 4.1
    3. 8.3 Analysis for PCA
      1. 8.3.1 Step 1: Characterizing the Bias
      2. 8.3.2 Step 2: Eigenvalue Concentration
      3. 8.3.3 Step 3: Finding the Leading-Order Expansion
      4. 8.3.4 Step 4: Completing the Proofs of Theorems 3.3 and 3.4
      5. 8.3.5 Step 6: Variance Approximation and Proof of Theorem 4.2
  9. A Proofs of Matrix Analysis Results
    1. A.1 Proof of Theorem 8.1
    2. A.2 Proof of Theorem 8.2
  10. B Proofs of Matrix Denoising Lemmas and Corollaries
    1. B.1 Proof of Lemma 8.1
    2. B.2 Proof of Lemma 8.2
    3. B.3 Proof of Lemma 8.3
    4. B.4 Proof of Lemma 8.4
    5. B.5 Proof of Lemma 8.5
    6. B.6 Proof of Lemma 8.6
    7. B.7 Proof of Corollary 3.1
    8. B.8 Proof of Corollary 4.1
  11. C Proofs of PCA Lemmas and Corollaries
    1. C.1 Proof of Lemma 8.8
    2. C.2 Proof of Lemma 8.9
    3. C.3 Proof of Lemma 8.10
    4. C.4 Proof of Lemma 8.11
    5. C.5 Proof of Lemma 8.12
    6. C.6 Proof of Lemma 8.13
    7. C.7 Proof of Lemma 8.14
    8. C.8 Proof of Corollary 3.2
    9. C.9 Proof of Corollary 4.2
  12. D Proofs of Additional Lemmas from Appendix B
    1. D.1 Proof of Lemma B.2
    2. D.2 Proof of Lemma B.3
    3. D.3 Proof of Lemma B.4
    4. D.4 Proof of Lemma B.5
    5. D.5 Proof of Lemma B.6
    6. D.6 Proof of Lemma B.7
    7. D.7 Proof of Lemma B.8
    8. D.8 Proof of Lemma B.9
  13. E Proofs of Additional Lemmas from Appendix C
    1. E.1 Proof of Lemma C.2
    2. E.2 Proof of Lemma C.3
    3. E.3 Proof of Lemma C.4
    4. E.4 Proof of Lemma C.5
    5. E.5 Proof of Lemma C.6
    6. E.6 Proof of Lemma C.7
    7. E.7 Proof of Lemma C.8
    8. E.8 Proof of Lemma C.9
    9. E.9 Proof of Lemma C.10

##  1 Introduction

Spectral methods, or algorithms and procedures that rely on eigenvectors,
singular vectors, and related quantities, are a useful collection of tools for
analyzing high-dimensional data, with applications to clustering (Löffler et
al., 2021), network analysis (Jin et al., 2023; Mao et al., 2020), ranking
(Chen et al., 2019, 2022), principal component analysis (Johnstone and Lu,
2009), nonconvex optimization (Chi et al., 2019), and tensor data analysis
(Zhang and Xia, 2018), to name a few. With the growing ubiquity of spectral
methods, there has been an increased interest in providing statistical
analyses of these techniques that take into account both the structural
properties of the problem at hand as well as the particularities of the noise,
in some cases resulting in provably optimal guarantees for various spectral-
based algorithms (Abbe et al., 2020, 2022; Ling, 2022; Nguyen and Zhang, 2023;
Zhang, 2022, 2023).

One common paradigm to analyze spectral methods is via classical deterministic
perturbation theory, which studies the $\ell_{2}$ errors between eigenvectors
of some underlying matrix $\bm{M}\in\mathbb{R}^{n\times n}$ and the
eigenvectors of a perturbed matrix $\bm{\widehat{M}}=\bm{M}+\bm{E}$ as a
function of the perturbation matrix $\bm{E}$ and the eigenvalues or singular
values of $\bm{M}$. By combining these techniques with tools from
nonsaymptotic random matrix theory one can derive high-probability or expected
$\ell_{2}$ error rates for eigenvectors and singular vectors (Cai and Zhang,
2018). Furthermore, there has been a growing interest in $\ell_{2,\infty}$
eigenspace and singular subspace perturbation theory (Cape et al., 2019b),
which quantifies the entrywise and row-wise fluctuations of subspaces through
careful probabilistic analysis of the noise. By carefully tailoring these
analyses to the distributional properties of the noise matrix $\bm{E}$,
estimation error rates can be obtained that go beyond what classical
perturbation theory offers. In some cases, these entrywise guarantees can
yield asymptotic normality for estimators of subspaces (Agterberg and Cape,
2023; Cape et al., 2019a), and these tools can be applied in various contexts,
resulting in statistically sound inferential procedures (Chen et al., 2021b;
Gao et al., 2023; Yan et al., 2021).

However, in many contexts entrywise perturbation and distributional theory is
still not sufficiently fine-grained to characterize the fluctuations of
individual eigenvectors. Therefore, in this paper we contribute to the
literature by studying statistical inference for a general linear function of
a single eigenvector. Explicitly, given a pre-specified unit vector
$\bm{a}\in\mathbb{R}^{n}$ and an eigenvector of interest
$\bm{u}_{j}\in\mathbb{R}^{n}$, we study estimating the linear form
$\bm{a}^{\top}\bm{u}_{j}$ with the empirical eigenvector
$\bm{\widehat{u}}_{j}$. There have been a few works in this area (Chen et al.,
2021a; Cheng et al., 2021; Koltchinskii and Xia, 2016; Koltchinskii and
Lounici, 2016; Koltchinskii et al., 2020), though we provide a more thorough
review of related work in Section 5. One key downside of existing theory is
that theoretical results often require highly pessimistic assumptions on the
_eigengap_ , the minimal separation between distinct eigenvalues. Therefore,
in this work we contribute to the literature by studying statistical inference
in regimes where the eigengap can be _significantly smaller_ than previously
assumed (possibly up to order $n^{-1/2}$ smaller, ignoring logarithmic terms).

Our analysis is motivated by the findings of Li et al. (2022) who provide
perturbation bounds for estimators of $\bm{a}^{\top}\bm{u}_{j}$ in the
presence of small eigengaps. They show that a form of “debiasing” is necessary
for the empirical estimate $\bm{a}^{\top}\bm{\widehat{u}}_{j}$ when
$|\bm{a}^{\top}\bm{u}_{j}|$ is large, and they provide a bias-corrected
estimator for $\bm{a}^{\top}\bm{u}_{j}$ that achieves the minimax rate (up to
some logarithmic factor), provided the noise variance is known _a priori_.
While their work demonstrated that it is possible to obtain minimax-optimal
estimators of linear functions of eigenvectors, in this work we go one step
further and analyze the distributional fluctuations of these estimators. We
provide data-driven confidence intervals for $\bm{a}^{\top}\bm{u}_{j}$, and
our results hold under nearly optimal signal-to-noise ratio conditions such
that consistent eigenspace estimation is possible. In addition, our estimators
can be computed directly from the data without the need for sample splitting.

In the following subsections we first describe the two models studied in this
work, and we then describe the primary contributions of this paper in greater
detail.

###  1.1 Two Canonical Settings

Throughout this work we focus on two fundamental statistical models: matrix
denoising and principal component analysis.

  * •

Matrix Denoising. Consider the setting where one is interested in the
eigenvectors of some symmetric rank $r$ matrix $\bm{S}\in\mathbb{R}^{n\times
n}$ corrupted with some symmetric Gaussian noise matrix $\bm{N}$, where
$\bm{N}_{ij}\sim\mathcal{N}(0,\sigma^{2})$ for $i<j$ and
$\bm{N}_{ii}\sim\mathcal{N}(0,2\sigma^{2})$ are independent (i.e., $\bm{N}$
belongs to the _Gaussian Orthogonal Ensemble_). Given a pre-specified
deterministic vector $\bm{a}$, in this paper we consider estimating
$\bm{a}^{\top}\bm{u}_{j}$, where $\bm{u}_{j}$ is the eigenvector associated to
the $j$’th largest in magnitude eigenvalue of $\bm{S}$.

  * •

Principal Component Analysis. Suppose one has $n$ i.i.d. observations
$\bm{X}_{i}\in\mathbb{R}^{p}$ with $\bm{X}_{i}\sim\mathcal{N}(0,\bm{\Sigma})$,
where $\bm{\Sigma}$ belongs to the class of _spiked covariance matrices_ ;
i.e., $\bm{\Sigma}=\bm{\Sigma}_{0}+\sigma^{2}\bm{I}_{p}$, where
$\bm{\Sigma}_{0}$ is a positive semidefinite rank $r$ matrix and $\bm{I}_{p}$
denotes the $p\times p$ identity. As in matrix denoising we consider
estimating $\bm{a}^{\top}\bm{u}_{j}$ for some pre-specified unit vector
$\bm{a}$, where $\bm{u}_{j}$ is the $j$’th eigenvector of $\bm{\Sigma}_{0}$
(and hence also $\bm{\Sigma}$).

These models have been studied explicitly in the literature previously. For
example, Bao et al. (2021); Knowles and Yin (2014); Xia (2021) study the
asymptotic properties of empirical eigenvectors and singular vectors, and Bao
et al. (2022); Koltchinskii and Lounici (2016); Koltchinskii et al. (2020)
study the asymptotic properties of principal components. However, unlike the
aforecited works, our emphasis is on providing finite-sample distributional
guarantees and data-driven uncertainty quantification in regimes where the
eigengaps (see Section 2 for the formal definition) can be significantly
smaller than considered in previous works.

Both of these models were studied in Li et al. (2022) who also consider
estimating $\bm{a}^{\top}\bm{u}_{j}$. However, unlike their work, in this work
we focus on distributional theory and statistical inference for estimators of
$\bm{a}^{\top}\bm{u}_{j}$, which requires novel analysis.

###  1.2 Contributions

Throughout this work we provide analogous results for both matrix denoising
and principal component analysis. Our main contributions are as follows:

  * •

We establish nonasymptotic distributional theory in the form of Berry-Esseen
theorems for linear functions of eigenvectors
$\bm{a}^{\top}\bm{\widehat{u}}_{j}$ in the presence of small eigengaps about
the quantity $\bm{a}^{\top}\bm{u}_{j}\bm{u}_{j}^{\top}\bm{\widehat{u}}_{j}$,
which is a biased centering term. We then provide similar results for a bias-
corrected estimator about the quantity $\bm{a}^{\top}\bm{u}_{j}$, where the
bias correction is computed directly from data assuming knowledge of the noise
variance $\sigma^{2}$. Our Berry-Esseen bounds clearly demonstrate the effect
of bias, signal strength, and small eigengaps on the rate of convergence. As a
corollary we provide entrywise distributional theory for individual
eigenvectors.

  * •

We then study the problem of statistical inference for
$\bm{a}^{\top}\bm{u}_{j}$, and we provide fully data-driven estimators for the
bias and variance without assuming knowledge of the noise variance
$\sigma^{2}$. We show that these estimators yield approximately valid
confidence intervals, and we apply these results to obtain approximately valid
uncertainty quantification for entries of individual eigenvectors.

Unlike previous work, our procedures are fully data-driven and do not require
any sample splitting, and to the best of our knowledge, these results are the
first to provide valid confidence intervals for linear functions of
eigenvectors in these settings under nearly optimal signal-to-noise ratio
conditions without any sample splitting. It is worth emphasizing that our
results continue to hold even in the “large eigengaps” regime, and our
distributional results reveal a novel phenomenon on the role of “close
eigenvalues” in the rate of convergence. Our proofs are based on higher-order
asymptotic expansions that allow us to identify leading-order terms and
isolate several residual terms that are closely related to certain “eigenvalue
bias” quantities. Together with several deterministic matrix-analytic results,
the form of the residual terms allow us to decouple statistical randomness
even when eigengaps may be small, enabling the sharp concentration required to
prove our main results.

###  1.3 Organization

The rest of this paper is organized as follows. In Section 2 we formally
define our models as well as several fundamental quantities that play a key
role in our theory. In Section 3 we state our distributional theory results
for both the biased and unbiased estimators, and Section 4 we describe our
inferential procedures and their associated theoretical guarantees. We cover
previous work in Section 5, we provide numerical simulations in Section 6, and
we discuss in Section 7. The proofs of our main distributional theory and
statistical inference results are given in Section 8. Proofs of additional
corollaries and lemmas are deferred to the appendices.

###  1.4 Notation

We let $[n]$ denote the set of integers $\\{1,\dots,n\\}$. For two numbers $a$
and $b$ we write $a\wedge b$ to mean the minimum of $a$ and $b$ and $a\vee b$
as the maximum. For two sequences $a_{n}$ and $b_{n}$, we write $a_{n}\lesssim
b_{n}$ or $a_{n}=O(b_{n})$ if there exists a universal constant $C>0$ such
that $a_{n}\leq Cb_{n}$, and we write $a_{n}\ll b_{n}$ if $a_{n}/b_{n}\to 0$
as $n\to\infty$. We write $a_{n}\asymp b_{n}$ and if $a_{n}\lesssim b_{n}$ and
$b_{n}\lesssim a_{n}$. We set $a_{n}=\Omega(b_{n})$ if $b_{n}=O(a_{n})$. We
write $a_{n}=\widetilde{O}(b_{n})$ if there exists some $c>0$ not depending on
$n$ such that $a_{n}\lesssim b_{n}\log^{c}(n)$, and similarly for
$a_{n}=\widetilde{\Omega}(b_{n})$. For a matrix $\bm{U}$ with orthonormal
columns we write $\bm{U}_{\perp}$ to denote the (not necessarily unique)
matrix satisfying $\bm{U}_{\perp}^{\top}\bm{U}=0$ such that
$[\bm{U},\bm{U}_{\perp}]$ is an orthogonal matrix. We write $\bm{I}_{n}$ for
the $n\times n$ identity matrix, $\|\cdot\|$ as the spectral norm for matrices
and Euclidean norm for vectors, $\|\cdot\|_{F}$ as the Frobenius norm on
matrices, and we let $\|\cdot\|_{\infty}$ denote the $\ell_{\infty}$ norm on
vectors. We let $\bm{e}_{i}$ denote the $i$’th standard basis vector whose
dimension will be clear from context. For a given vector $\bm{v}$, we let
$\bm{v}(i)$ denote its $i$’th entry. For a random variable $X$, we let
$\|X\|_{\psi_{2}}$ be its subgaussian Orlicz norm defined via
$\|X\|_{\psi_{2}}:=\inf\\{t>0:\mathbb{E}\exp(X^{2}/t^{2})\leq 2\\}$.

##  2 Statistical Models and Key Quantities

In this section we formulate our two statistical models and introduce several
important quantities that are particularly relevant in our main results.

###  2.1 Matrix Denoising

Suppose one is given a matrix $\bm{\widehat{S}}\in\mathbb{R}^{n\times n}$ of
the form

$\displaystyle\bm{\widehat{S}}$ $\displaystyle=\bm{S}+\bm{N}$ (1)

where $\bm{S}$ is a rank $r$ signal matrix with eigendecomposition given by

$\displaystyle\bm{S}$
$\displaystyle=\bm{U}\bm{\Lambda}\bm{U}^{\top}=\sum_{i=1}^{r}\lambda_{i}\bm{u}_{i}\bm{u}_{i}^{\top},$

and $\bm{N}$ is a symmetric noise matrix with independent entries given by

$\displaystyle\bm{N}_{ij}\sim\begin{cases}\mathcal{N}(0,\sigma^{2})&i<j,\\\
\mathcal{N}(0,2\sigma^{2})&i=j.\end{cases}$

The matrix $\bm{N}$ is said to belong to the _Gaussian Orthogonal Ensemble_
(GOE). A key property of the GOE is that it is rotationally invariant; i.e.,
given a matrix $\bm{R}$ satisfying
$\bm{RR}^{\top}=\bm{R}^{\top}\bm{R}=\bm{I}_{n}$, it holds that
$\bm{R}\bm{N}\bm{R}^{\top}$ has the same distribution as $\bm{N}$.

We define

$\displaystyle\lambda_{\min}=\min_{1\leq i\leq
r}|\lambda_{i}|;\qquad\lambda_{\max}=\max_{1\leq i\leq r}|\lambda_{i}|,$

representing the magnitude of the smallest and largest nonzero eigenvalues of
$\bm{S}$, and we set $\kappa$ as the _reduced condition number_
$\kappa=\lambda_{\max}/\lambda_{\min}$. Similarly, we let $\bm{\widehat{S}}$
have eigenvectors $\bm{\widehat{u}}_{k}$ and eigenvalues
$\widehat{\lambda}_{k}$. For a given index $j$ with $1\leq j\leq r$, we define
the eigenvalue gap (or eigengap) via

$\displaystyle\Delta_{j}$ $\displaystyle=\min_{k\neq j,1\leq k\leq
r}|\lambda_{j}-\lambda_{k}|,$

with the convention that $\Delta_{j}=\lambda_{\min}$ if $r=1$. We also set

$\displaystyle\Delta_{\min}=\min_{1\leq j\leq r}\Delta_{j}.$

We let $\bm{\widehat{U}}$ denote the $n\times r$ matrix of leading
eigenvectors of $\bm{\widehat{S}}$ organized according to the magnitude of the
eigenvalues of $\bm{\widehat{S}}$. We set $\bm{\widehat{u}}_{j}$ to be the
$j$’th column of $\bm{\widehat{U}}$.

We will use the following setting as a running example.

######  Example 1 (Gaussian Mixed-Membership Blockmodel).

Consider the case that one is given a matrix $\bm{S}$ of the form

$\displaystyle\bm{S}$ $\displaystyle=\frac{1}{n}\bm{Z}\bm{B}\bm{Z}^{\top},$

where $\bm{Z}\in[0,1]^{n\times 2}$ has rows that sum to one and
$\bm{B}=\bm{V}\mathrm{Diag}(\lambda+\varepsilon,\lambda)\bm{V}^{\top}$ for
some $2\times 2$ orthonormal matrix $\bm{V}$. Here the matrix
$\frac{1}{n}\bm{B}$ is interpreted as a matrix of “community means,” with
$\bm{Z}$ playing the role of the “membership matrix,” where each row of
$\bm{Z}$ denotes the membership of the $i$’th vertex. The assumption that each
row of $\bm{Z}$ sums to one corresponds to a setting wherein each vertex
belongs to a convex combination of communities, with “pure” membership
corresponding to the row being strictly $\\{0,1\\}$-valued. It is
straightforward to demonstrate that $\bm{S}$ has matrix of eigenvectors equal
to

$\displaystyle\bm{U}$
$\displaystyle=\bm{Z}(\bm{Z}^{\top}\bm{Z})^{-1/2}\bm{V}.$

Furthermore, the smallest eigenvalue of $\bm{B}$ is of order $\lambda$ as long
as $\lambda_{\min}(\bm{Z}^{\top}\bm{Z})\asymp n$, which can be interpreted as
stating that each community is of “approximately equal size.” For
concreteness, we will consider the following model. Suppose $\bm{Z}$ is the
matrix with the first $(1-\alpha)n/2$ rows equal to $(1,0)$, the second
$(1-\alpha)n/2$ rows equal to $(0,1)$, and the final $\alpha n$ rows equal to
$(.5,.5)$. Then

$\displaystyle\bm{Z}^{\top}\bm{Z}$
$\displaystyle=\begin{pmatrix}(1-\alpha)n/2&\frac{1}{4}\alpha n\\\
\frac{1}{4}\alpha n&(1-\alpha)n/2\end{pmatrix},$

so that
$\lambda_{\min}(\bm{Z}^{\top}\bm{Z})=n\big{(}\frac{1}{2}-\frac{3}{4}\alpha\big{)}$
and
$\lambda_{\max}(\bm{Z}^{\top}\bm{Z})=n\big{(}\frac{1}{2}-\frac{1}{4}\alpha\big{)}$.
Consequently, if $\alpha=o(1)$, then $\lambda_{\min}(\bm{S})\gtrsim\lambda$
and

$\displaystyle\Delta_{1}=\lambda_{\max}(\bm{S})-\lambda_{\min}(\bm{S})=\frac{1}{n}\lambda_{\max}\big{(}\bm{Z}^{\top}\bm{Z}\bm{B}\big{)}-\frac{1}{n}\lambda_{\min}\big{(}\bm{Z}^{\top}\bm{Z}\bm{B}\big{)}$
$\displaystyle\geq\frac{\varepsilon}{100},$

which holds as long as $\alpha\asymp c\frac{\varepsilon}{\lambda}$ and
$\varepsilon\ll\lambda$.

In this work we consider several related problems:

  * •

Distributional theory for $\bm{a}^{\top}\bm{\widehat{u}}_{j}$ Suppose one is
given a fixed (deterministic) vector $\bm{a}$ satisfying $\|\bm{a}\|=1$.
Consider estimating $\bm{a}^{\top}\bm{u}_{j}$ based on the observation
$\bm{\widehat{S}}$ for some $1\leq j\leq r$. We will study the approximate
Gaussianity of $\bm{a}^{\top}\bm{\widehat{u}}_{j}$ about a biased centering
term of the form
$\bm{a}^{\top}\bm{u}_{j}\bm{u}_{j}^{\top}\bm{\widehat{u}}_{j}$, which
necessitates some form of debiasing procedure.

  * •

Distributional theory for the debiased quantity: We then study the approximate
Gaussianity of the debiased quantity $\sqrt{1+b_{j}^{{\sf
MD}}}\bm{a}^{\top}\bm{\widehat{u}}_{j}$, where $b_{j}^{{\sf MD}}$ is a debias
parameter computed from the data assuming knowledge of the noise variance
$\sigma^{2}$.

  * •

Entrywise Distributional Theory: We also consider the particular case
$\bm{a}=\bm{e}_{i}$, where $\bm{e}_{i}$ is a standard basis vector, which
corresponds to the $i$’th entry of $\bm{u}_{j}$. In this case, under
_incoherence conditions_ (see Corollary 3.1 for the formal definition),
debiasing is not necessary.

  * •

Uncertainty quantification: In Section 4 we study the problem of obtaining
confidence intervals for $\bm{a}^{\top}\bm{u}_{j}$ as well as entrywise
confidence intervals without assuming knowledge of the noise variance
$\sigma^{2}$.

###  2.2 Principal Component Analysis

Consider the setting where one is given $n$ observations
$\bm{X}_{i}\in\mathbb{R}^{p}$ satisfying

$\displaystyle\bm{X}_{i}\sim\mathcal{N}(0,\bm{\Sigma}),$ (2)

where $\bm{\Sigma}$ has the form

$\displaystyle\bm{\Sigma}=\bm{\Sigma}_{0}+\sigma^{2}\bm{I}_{p},$

with

$\displaystyle\bm{\Sigma}_{0}=\bm{U}\bm{\Lambda}\bm{U}^{\top}=\sum_{i=1}^{r}\lambda_{i}\bm{u}_{i}\bm{u}_{i}^{\top},$

where each $\lambda_{i}$ satisfies
$0<\lambda_{r}\leq\lambda_{r-1}\leq\cdots\leq\lambda_{1}$. We denote

$\displaystyle\lambda_{\min}=\lambda_{r};\qquad\lambda_{\max}=\lambda_{1},$

and we define $\kappa$ as the _reduced condition number_
$\kappa:=\lambda_{\max}/\lambda_{\min}$. Similar to matrix denoising, for a
given index $j$ with $1\leq j\leq r$, we denote the eigengap associated to the
index $j$ via

$\displaystyle\Delta_{j}=\min_{k\neq j,1\leq k\leq
r}|\lambda_{j}-\lambda_{k}|,$

with $\Delta_{j}=\lambda_{\min}$ if $r=1$. We also set

$\displaystyle\Delta_{\min}=\min_{1\leq j\leq r}\Delta_{j}.$

We consider estimating $\bm{\Sigma}$ with the sample covariance
$\bm{\widehat{\Sigma}}$ defined via

$\displaystyle\bm{\widehat{\Sigma}}$
$\displaystyle\coloneqq\frac{1}{n}\sum_{i=1}^{n}\bm{X}_{i}\bm{X}_{i}^{\top}=\frac{\bm{XX}^{\top}}{n},$

where $\bm{X}\in\mathbb{R}^{p\times n}$ is the matrix whose columns are the
observations $\bm{X}_{i}$. We let $\bm{\widehat{\Sigma}}$ have eigenvectors
and eigenvalues $\\{\bm{\widehat{u}}_{k},\widehat{\lambda}_{k}\\}_{k=1}^{p}$.
We let $\bm{\widehat{U}}$ denote the $p\times r$ matrix of leading
eigenvectors of $\bm{\widehat{\Sigma}}$, and we set $\bm{\widehat{u}}_{j}$ to
be the $j$’th column of $\bm{\widehat{U}}$.

For concreteness is it useful to think of the following two settings.

######  Example 2 (“Diagonal“ versus “Dense” Covariance).

Suppose that $\bm{\Sigma}_{0}$ is the $p\times p$ matrix with first two
diagonal entries equal to $\lambda+\varepsilon$ and $\lambda$ respectively,
and all others zero. Then $\bm{u}_{i}=\bm{e}_{i}$, and
$\Delta_{1}=\Delta_{\min}=\varepsilon$. We refer to this as the _diagonal
covariance model_.

Now suppose that that $p$ is even and suppose further that
$\bm{\Sigma}_{0}=(\lambda+\varepsilon)\bm{u}_{1}\bm{u}_{1}^{\top}+\lambda\bm{u}_{2}\bm{u}_{2}^{\top}$,
where $\bm{u}_{1}$ is the vector with entries all equal to
$\frac{1}{\sqrt{p}}$ and $\bm{u}_{2}$ is the vector with entries
$\pm\frac{1}{\sqrt{p}}$ such that $\bm{u}_{1}^{\top}\bm{u}_{2}=0$. Then one
again has $\Delta_{1}=\Delta_{\min}=\varepsilon$. In this setting observe that
$\bm{\Sigma}_{0}$ is a dense block matrix, and we refer to this as the _dense
covariance model_.

As in the case for matrix denoising, we consider several related problems:

  * •

Distributional theory for $\bm{a}^{\top}\bm{\widehat{u}}_{j}$ and its debiased
counterpart: Suppose one is given a fixed (deterministic) vector $\bm{a}$
satisfying $\|\bm{a}\|=1$, and consider estimating $\bm{a}^{\top}\bm{u}_{j}$
based on the observation $\bm{\widehat{\Sigma}}$ for some $1\leq j\leq r$. We
will study the approximate Gaussianity of $\bm{a}^{\top}\bm{\widehat{u}}_{j}$
about the biased centering term of the form
$\bm{a}^{\top}\bm{u}_{j}\bm{u}_{j}^{\top}\bm{\widehat{u}}_{j}$, as well as the
approximate Gaussianity of the debiased quantity $\sqrt{1+b_{j}^{{\sf
PCA}}}\bm{a}^{\top}\bm{\widehat{u}}_{j}$, where $b_{j}^{{\sf PCA}}$ is a bias
parameter computed from the data assuming knowledge of the noise variance
$\sigma^{2}$.

  * •

Entrywise distributional theory and uncertainty quantification: We further
consider the entrywise distribution of $\bm{\widehat{u}}_{j}$, and we show
that debiasing is not strictly necessary in this setting. In addition, in
Section 4 we study how to obtain asymptotically valid confidence intervals for
$\bm{a}^{\top}\bm{u}_{j}$ as well as entrywise confidence intervals without
assuming knowledge of the noise variance $\sigma^{2}$.

##  3 Distributional Theory

With our models defined we are now prepared to discuss our first results,
namely, our Berry-Esseen bounds. In Section 3.1 we study matrix denoising and
in Section 3.2 we study PCA.

Throughout all of our results, we note that there is an eigenvector ambiguity
issue since the sign of an eigenvector can always be flipped. We therefore
assume without loss of generality that $\bm{\widehat{u}}_{j}$ is chosen such
that $\bm{\widehat{u}}_{j}^{\top}\bm{u}_{j}$ is nonnegative.

###  3.1 Matrix Denoising

The following result provides a distributional characterization for the plug-
in estimator $\bm{a}^{\top}\bm{\widehat{u}}_{j}$.

######  Theorem 3.1 (Approximate Gausssianity for Linear Functions of
Eigenvectors).

Consider the model (1), and suppose $\lambda_{j}$ is unique. Suppose that

$\displaystyle\lambda_{\min}/\sigma\geq
C_{0}\sqrt{rn}\qquad\Delta_{j}/\sigma\geq C_{0}r\log^{2}(n)$

where $C_{0}$ is some sufficiently large constant. Suppose further that $r\leq
c_{0}n/\log^{2}(n)$, where $c_{0}$ is a sufficiently small constant. Let
$\bm{a}$ be any determinstic unit vector such that $\bm{a}\neq\pm\bm{u}_{j}$.
Let $\bm{U}^{(k)}$ denote the matrix whose columns are the eigenvectors
corresponding to $\lambda_{k}$. Define

$\displaystyle(s_{\bm{a},j}^{{\sf MD}})^{2}\coloneqq 2\sum_{k\neq j,k\leq
r^{\prime}}\frac{\sigma^{2}\|\bm{a}^{\top}\bm{U}^{(k)}\|^{2}}{(\lambda_{j}-\lambda_{k})^{2}}+2\frac{\sigma^{2}\|\bm{U}_{\perp}^{\top}\bm{a}\|^{2}}{\lambda_{j}^{2}},$

where the summation is over all $r^{\prime}\leq r$ distinct eigenvalues. Let
$\Phi(z)$ denote the cumulative distribution function of a standard Gaussian
random variable. Then it holds that

$\displaystyle\sup_{z\in\mathbb{R}}\bigg{|}\mathbb{P}\bigg{\\{}\frac{1}{s_{\bm{a},j}^{{\sf
MD}}}\big{(}\bm{a}^{\top}\bm{\widehat{u}}_{j}-\bm{a}^{\top}\bm{u}_{j}\bm{u}_{j}^{\top}\bm{\widehat{u}}_{j}\big{)}\leq
z\bigg{\\}}-\Phi(z)\bigg{|}$ $\displaystyle\lesssim{\sf ErrMD}$

where

$\displaystyle{\sf
ErrMD}:=\underbrace{\frac{r\sqrt{n\log(n)}}{\lambda_{\min}/\sigma}}_{\text{Subspace
Estimation
Effect}}+\underbrace{\frac{r^{3/2}\log(n)}{\Delta_{j}/\sigma}}_{\text{Small
Eigengaps Effect}}+\ n^{-9}.$

The assumption that $\bm{a}\neq\pm\bm{u}_{j}$ is necessary, since otherwise
$\bm{a}^{\top}\bm{\widehat{u}}_{j}-\bm{a}^{\top}\bm{u}_{j}\bm{u}_{j}^{\top}\bm{\widehat{u}}_{j}\equiv
0$ and $s_{\bm{a},j}^{{\sf MD}}=0$, so the limiting distribution may no longer
be Gaussian. Therefore, this condition is not particularly stringent. In
addition, the factor $n^{-9}$ can be replaced with any negative power of $n$
at the cost of increasing the implicit constants, so in general it can be
ignored.

######  Remark 3.1 (Signal-Strength Condition).

Observe that the bound in Theorem 3.1 is only non-vacuous if
$\lambda_{\min}/\sigma\gg r\sqrt{n\log(n)}$. It has been demonstrated in Cai
and Zhang (2018) that $\lambda_{\min}/\sigma\gg\sqrt{rn}$ is both neccessary
and sufficient for consistent subspace estimation in $\|\cdot\|_{F}$ when
$r\ll n$; consequently the condition $\lambda_{\min}/\sigma\gg
r\sqrt{n\log(n)}$ is optimal up to the factor of $\sqrt{r\log(n)}$.

######  Remark 3.2 (Eigengap Conditions).

Again the bound is only non-vacuous if $\Delta_{j}/\sigma\gg r^{3/2}\log(n)$.
Together with our assumption on the eigengap $\Delta_{j}/\sigma\geq
C_{0}r\log^{2}(n)$, we see that a sufficient condition for normality is
$\Delta_{j}/\sigma\gg r^{3/2}\log^{2}(n)$. However, this eigengap condition is
significantly weaker than what is required in classical theory. For example,
the Davis-Kahan Theorem requires that

$\displaystyle\Delta_{j}\gtrsim\|\bm{N}\|\asymp\sigma\sqrt{n},$

which is more stringent by a factor of
$\widetilde{\Omega}\big{(}\sqrt{\frac{n}{r^{3}}}\big{)}$.

######  Remark 3.3 (Rate of Convergence).

Ignoring the additional factor of $n^{-9}$, when $r\asymp 1$, the rate of
convergence to asymptotic normality is of the form

$\displaystyle\max\bigg{\\{}\frac{\sqrt{n}}{\lambda_{\min}/\sigma},\frac{1}{\Delta_{j}/\sigma}\bigg{\\}},$

ignoring logarithmic terms, where we dub the first term as the “subspace
estimation effect” and the second term as the “small eigengaps effect.”
Consequently, Theorem 3.1 reveals a novel phenomenon – when the signal
strength is large but eigengaps are small, the primary difficulty in
estimation is separating the $j$’th eigenvector from the other leading $r-1$
eigenvectors, whereas when the eigengaps are large but the signal strength is
small, the primary difficulty lies in estimating the subspace $\bm{U}$, which
corresponds to separating the leading $r$ eigenvectors from the bottom $n-r$
eigenvectors.

######  Remark 3.4 (Biased Estimation).

Observe that the centering term in Theorem 3.1 is given by
$\bm{a}^{\top}\bm{u}_{j}\bm{u}_{j}^{\top}\bm{\widehat{u}}_{j}$. In general
$\bm{\widehat{u}}_{j}^{\top}\bm{u}_{j}$ need not equal one; in fact, in the
regime $\lambda_{\min}/\sigma\gg\sqrt{nr}$ it can be shown that
$\bm{\widehat{u}}_{j}^{\top}\bm{u}_{j}=1-o(1)$ and when
$\lambda_{\min}/\sigma\asymp\sqrt{nr}$ then
$\bm{\widehat{u}}_{j}^{\top}\bm{u}_{j}\approx c<1$ (e.g., Benaych-Georges and
Nadakuditi (2011)). Consequently, Theorem 3.1 demonstrates that
$\bm{a}^{\top}\bm{\widehat{u}}_{j}$ is a _biased_ estimate of
$\bm{a}^{\top}\bm{u}_{j}$ with bias given by
$\bm{a}^{\top}\bm{u}_{j}(1-\bm{\widehat{u}}_{j}^{\top}\bm{u}_{j})$, a
phenomenon that has also been demonstrated for asymmetric noise in Chen et al.
(2021a). When $|\bm{a}^{\top}\bm{u}_{j}|$ is sufficiently small relative to
$s_{\bm{a},j}^{{\sf MD}}$ this bias can be negligible (as we demonstrate in
the proof Corollary 3.1, which focuses on the case $\bm{a}=\bm{e}_{i}$), but
when $|\bm{a}^{\top}\bm{u}_{j}|\asymp 1$, this bias may be significant.

######  Remark 3.5 (Comparison to Previous Works).

Our work is most similar to the works Li et al. (2022) and Koltchinskii and
Xia (2016), who study perturbation bounds for the estimator
$\bm{a}^{\top}\bm{\widehat{u}}_{j}$ in similar contexts, and our assumptions
on the signal strength $\lambda_{\min}/\sigma$ match theirs up to logarithmic
terms and factors of $r$. However, only Li et al. (2022) consider the small
eigengaps regime. Our condition on $\Delta_{j}$ matches theirs up to
logarithmic terms and factors of $r$, while simultaneously providing a sharper
characterization of the estimator $\bm{a}^{\top}\bm{\widehat{u}}_{j}$ in terms
of nonasymptotic distributional theory.

In order to prove our results we must must identify the leading-order term,
compute its variance $(s_{\bm{a},j}^{{\sf MD}})^{2}$, and then demonstrate
that the residual terms are of order $o(s_{\bm{a},j}^{{\sf MD}})$. In
contrast, both of these previous works only needed to show that the difference
is of order $O(s_{\bm{a},j}^{{\sf MD}})$. Therefore, our analysis must extend
further (implicitly considering higher-order terms) to show that these
residuals are all sufficiently small even with small eigengaps. To accomplish
this goal we leverage heavily two deterministic matrix analysis results
(Theorems 8.1 and 8.2) that allow us to decouple several statistically
dependent quantities as well as isolate certain “eigenvalue bias” terms (see
Lemma 8.2). See Section 8.2.3 for further details on how these eigenvalue bias
terms arise in the analysis. Identifying how these eigenvalue bias terms arise
plays a key role in providing the sharp concentration guarantees required to
prove Theorem 3.1.

######  Remark 3.6 (Cosines).

Suppose $\bm{a}=\bm{u}_{l}$ for some $l\neq j$ with $l\leq r$, and suppose
$\lambda_{l}$ is unique. Then note that $(s_{\bm{u}_{l},j}^{{\sf
MD}})^{2}=\frac{\sigma^{2}}{(\lambda_{j}-\lambda_{l})^{2}}$. Theorem 3.1
demonstrates that as $n\to\infty$,

$\displaystyle\frac{|\lambda_{j}-\lambda_{l}|}{\sigma}\cos\theta(\bm{u}_{l},\bm{\widehat{u}}_{j})\to\mathcal{N}(0,1),$

as long as $\lambda_{\min}/\sigma\gg r\sqrt{n\log(n)}$ and
$\Delta_{j}/\sigma\gg r^{3/2}\log(n)$. To the best of our knowledge, our
result is the first establishing the convergence of the cosines between
leading eigenvectors when the order of the eigengaps can be significantly
smaller than the order of the spike.

######  Example 3 (Gaussian Mixed-Membership Blockmodel, Continued).

In the setting of Example 1 consider the case
$\varepsilon\asymp\sigma\mathrm{polylog}(n)$ and
$\lambda\asymp\sigma\sqrt{n}$. Then for any $\bm{a}$ such that
$\bm{a}\neq\pm\bm{u}_{j}$, it holds that
$\bm{a}^{\top}\bm{\widehat{u}}_{j}-\bm{a}^{\top}\bm{u}_{j}\bm{u}_{j}^{\top}\bm{\widehat{u}}_{j}$
is asymptotically Gaussian. In particular, consider testing whether $\alpha=0$
(i.e., testing whether all the vertices are “pure” or not). If $\alpha=0$, it
holds that $\bm{u}_{1}$ is the constant vector; and hence by taking $\bm{a}$
as the constant vector it holds that $\bm{a}^{\top}\bm{\widehat{u}}_{2}$ is
asymptotically Gaussian about zero if $\alpha=0$; otherwise
$\bm{a}^{\top}\bm{\widehat{u}}_{2}$ is asymptotically Gaussian about some
nonzero quantity.

Our next result provides distributional guarantees for a bias-corrected
estimate of $\bm{a}^{\top}\bm{u}_{j}$ assuming knowledge of the noise variance
$\sigma^{2}$.

######  Theorem 3.2 (Asymptotic Normality of Debiased Estimators – Matrix
Denoising).

Consider the setting of Theorem 3.1. Define

$\displaystyle b_{j}^{{\sf MD}}:$
$\displaystyle=\sum_{k>r}\frac{\sigma^{2}}{(\widehat{\lambda}_{j}-\widehat{\lambda}_{k})^{2}},$

where $\sigma^{2}$ is assumed known and $\widehat{\lambda}_{k}$ are the
ordered empirical eigenvalues of $\bm{\widehat{S}}$. Then it holds that

$\displaystyle\sup_{z\in\mathbb{R}}\bigg{|}\mathbb{P}\bigg{\\{}\frac{1}{s_{\bm{a},j}^{{\sf
MD}}}\big{(}\bm{a}^{\top}\bm{\widehat{u}}_{j}\sqrt{1+b_{j}^{{\sf
MD}}}-\bm{a}^{\top}\bm{u}_{j}\big{)}\leq z\bigg{\\}}-\Phi(z)\bigg{|}$
$\displaystyle\lesssim{\sf
ErrMD}+\frac{|\bm{a}^{\top}\bm{u}_{j}|}{s_{\bm{a},j}^{{\sf MD}}}{\sf
ErrBiasMD},$

where ${\sf ErrMD}$ and $s_{\bm{a},j}^{{\sf MD}}$ are as in Theorem 3.1 and

$\displaystyle{\sf
ErrBiasMD}:=\frac{\sigma^{2}r\log(n)}{\Delta_{j}^{2}}+\frac{\sigma^{2}\sqrt{n\log(n)}}{\lambda_{j}^{2}}.$

######  Remark 3.7 (Bias Effect).

Observe that Theorem 3.2 exhibits an additional error reflecting the order of
$\bm{a}^{\top}\bm{u}_{j}$, whereas Theorem 3.1 does not depend on
$|\bm{a}^{\top}\bm{u}_{j}|$ except through the centering term. For simplicity
assume that $\bm{S}$ is positive semidefinite. Then it is straightforward to
demonstrate that

$\displaystyle\frac{|\bm{a}^{\top}\bm{u}_{j}|}{s_{\bm{a},j}^{{\sf MD}}}{\sf
ErrBiasMD}\lesssim\frac{\sigma\sqrt{r}\log(n)}{\Delta_{j}}\min_{k\neq j,k\leq
r}\frac{|\bm{a}^{\top}\bm{u}_{j}|}{|\bm{a}^{\top}\bm{u}_{k}|}\frac{|\lambda_{k}-\lambda_{j}|}{\Delta_{j}}+\frac{|\bm{a}^{\top}\bm{u}_{j}|}{\|\bm{U}_{\perp}^{\top}\bm{a}\|}\frac{\sigma\sqrt{n\log(n)}}{\lambda_{j}}.$

If it holds that $|\bm{a}^{\top}\bm{u}_{k}|\asymp|\bm{a}^{\top}\bm{u}_{j}|$
for all $k\neq j,k\leq r$ and $|\lambda_{k}-\lambda_{j}|\asymp\Delta_{j}$ for
all $k\neq j,k\leq r$, then Theorem 3.2 implies that

$\displaystyle\sup_{z\in\mathbb{R}}\bigg{|}\mathbb{P}\bigg{\\{}\frac{1}{s_{\bm{a},j}^{{\sf
MD}}}\big{(}\bm{a}^{\top}\bm{\widehat{u}}_{j}\sqrt{1+b_{j}^{{\sf
MD}}}-\bm{a}^{\top}\bm{u}_{j}\big{)}\leq z\bigg{\\}}-\Phi(z)\bigg{|}$
$\displaystyle\lesssim\frac{r\sqrt{n\log(n)}}{\lambda/\sigma}+\frac{r^{3/2}\log(n)}{\Delta_{j}/\sigma}$
$\displaystyle\quad+\frac{|\bm{a}^{\top}\bm{u}_{j}|}{\|\bm{U}_{\perp}^{\top}\bm{a}\|}\frac{\sigma\sqrt{n\log(n)}}{\lambda_{j}}+n^{-9}.$

The additional term incurs a factor of
$\frac{|\bm{a}^{\top}\bm{u}_{j}|}{\|\bm{U}_{\perp}^{\top}\bm{a}\|}$, which
quantifies the effect of $\bm{a}$ on the subspace corresponding to
$\bm{u}_{j}$ relative to the orthogonal complement of the leading subspace.
Consequently, this result can be understood as stating that
$\bm{a}^{\top}\bm{u}_{j}$ cannot be “too close” to $\pm 1$ to obtain
asymptotic normality with this plug-in estimator. In contrast, Theorem 3.1
allows $\bm{a}^{\top}\bm{u}_{j}$ to be arbitrarily close (but not equal) to
$\pm 1$ to obtain asymptotic normality.

We now turn our attention to the entrywise distribution of the individual
eigenvectors; i.e., when $\bm{a}=\bm{e}_{i}$. Under this setup, we have the
following corollary.

######  Corollary 3.1 (Entrywise Distributional Theory – Matrix Denoising).

In the setting of Theorem 3.2, let $\bm{a}=\bm{e}_{i}$; i.e., $\bm{a}$ is a
standard basis vector, and suppose that there exists a constant $\mu=O(1)$
such that

$\displaystyle\max_{k\leq r}|\bm{u}_{k}(i)|$
$\displaystyle\leq\mu\sqrt{\frac{1}{n}}.$

Suppose further that

$\max_{k\leq r,k\neq j}\frac{|\bm{u}_{j}(i)|}{|\bm{u}_{k}(i)|}=O(1).$

Let $s^{{\sf MD}}_{\bm{e}_{i},j}$ be as in Theorem 3.1. Then it holds that

$\displaystyle\sup_{z\in\mathbb{R}}\bigg{|}\mathbb{P}\bigg{\\{}\frac{1}{s^{{\sf
MD}}_{\bm{e}_{i},j}}\bigg{(}\bm{\widehat{u}}_{j}(i)-\bm{u}_{j}(i)\bigg{)}\leq
z\bigg{\\}}-\Phi(z)\bigg{|}$
$\displaystyle\lesssim\frac{r\sqrt{n\log(n)}}{\lambda_{\min}/\sigma}+\frac{r^{3/2}\log(n)}{\Delta_{j}/\sigma}+n^{-8}.$

The proof of this result can be found in Section B.7.

######  Remark 3.8 (Comparison to Previous Works).

The constant $\mu$ plays the role of the _incoherence constant_ that arises in
a number of problems related to entrywise estimation. The works Chen et al.
(2021b) and Agterberg et al. (2022) establish distributional theory for
individual rows of $\bm{\widehat{U}}\mathcal{O}-\bm{U}$, where $\mathcal{O}$
is an $r\times r$ orthogonal transformation that arises due to eigenvalue
multiplicity. Explicitly, a straightforward adaptation of their results yields
that

$\displaystyle\frac{|\bm{\Lambda}|}{\sigma}\bigg{(}\bm{\widehat{U}}\mathcal{O}-\bm{U}\bigg{)}_{i\cdot}\to\mathcal{N}(0,\bm{I}_{r}),$

when $r,\kappa,\mu_{0}=O(1)$, provided certain signal-strength conditions are
met. Therefore, these results imply that marginally

$\displaystyle\frac{|\lambda_{j}|}{\sigma}\bigg{(}\bm{\widehat{U}}\mathcal{O}-\bm{U}\bigg{)}_{ij}\to\mathcal{N}(0,1).$
(3)

Observe that Corollary 3.1 and (3) are of a similar form: both results
qualitatively provide the order of the fluctuations of entries of the $j$’th
eigenvector. Note that these previous works do not assume any conditions on
the eigengaps, in particular allowing eigenvalue multiplicity. However, in
these settings their results can fail to yield the precise distribution of
individual eigenvector entries when the eigengaps can be small, as (3) still
contains an additional orthogonal transformation that cannot in general be
eliminated (unless the eigengaps are sufficiently large). Therefore, our
result complements these previous results as Corollary 3.1 provides individual
eigenvector entry fluctuations without assuming stringent conditions on the
eigengaps. Our result can also be interpreted as a more “fine-grained” result
isolating the effect of “nearby” eigenvectors (as quantified by eigengap) on
individual eigenvectors.

######  Remark 3.9 (Extensions to Other “Negligible Bias” Regimes).

In Corollary 3.1 we focus only on the entrywise distribution, where
incoherence plays a crucial role providing the “unalignment” of
$\bm{a}=\bm{e}_{i}$ with $\bm{U}$. Corollary 3.1 can be generalized to
settings such that $\max_{k\leq r}|\bm{a}^{\top}\bm{u}_{k}|\lesssim n^{-1/2}$
with $|\bm{a}^{\top}\bm{u}_{j}|\lesssim|\bm{a}^{\top}\bm{u}_{k}|$. In such a
regime, bias correction is no longer necessary, as the bias is of order
$n^{-1/2}$, and hence is negligible. However, such a result requires the _a
priori_ knowledge that $\|\bm{a}^{\top}\bm{U}\|$ is sufficiently small, which
may not occur in practice without additional assumptions (such as
incoherence).

###  3.2 Principal Component Analysis

In this section we provide analogous results to those in Section 3.1 for
principal component analysis. The following result provides an analogue of
Theorem 3.1 for the principal component analysis model (2).

######  Theorem 3.3 (Approximate Gaussianity of Linear Functions of Principal
Components).

Consider the model in (2) and suppose $\lambda_{j}$ is unique. Suppose that
$r\leq c_{0}\frac{n}{\kappa^{2}\log^{4}(n\vee p)}$ for some sufficiently small
constant $c_{0}$, that $\log(p)\lesssim n$, and that

$\displaystyle\lambda_{\min}$ $\displaystyle\geq
C_{1}\max\bigg{\\{}\sigma^{2}\log^{2}(n\vee
p)\bigg{(}\frac{p}{n}+\sqrt{\frac{p}{n}}\bigg{)},\kappa\log^{4}(n\vee
p)\sigma\frac{p}{n}\Bigg{\\}};$ (4) $\displaystyle\Delta_{j}$
$\displaystyle\geq
C_{1}\big{(}\lambda_{\max}+\sigma^{2}\big{)}\sqrt{\frac{r}{n}}\log(n\vee p),$

where $C_{1}$ is a sufficiently large constant. Let $\bm{a}$ be any
determinstic unit vector such that $\bm{a}\neq\pm\bm{u}_{j}$. Let
$\bm{U}^{(k)}$ denote the eigenspace corresponding to $\lambda_{k}$. Define

$\displaystyle(s_{\bm{a},j}^{{\sf PCA}})^{2}$ $\displaystyle=\sum_{k\neq
j,k\leq
r^{\prime}}\frac{(\lambda_{j}+\sigma^{2})(\lambda_{k}+\sigma^{2})\|\bm{a}^{\top}\bm{U}^{(k)}\|^{2}}{n(\lambda_{j}-\lambda_{k})^{2}}+\frac{(\lambda_{j}+\sigma^{2})\sigma^{2}\|\bm{U}_{\perp}^{\top}\bm{a}\|^{2}}{n\lambda_{j}^{2}},$

where the summation is over all $r^{\prime}\leq r$ distinct eigenvalues. Then
it holds that

$\displaystyle\sup_{z\in\mathbb{R}}\bigg{|}\mathbb{P}\bigg{\\{}\frac{1}{s_{\bm{a},j}^{{\sf
PCA}}}\bigg{(}\bm{a}^{\top}\bm{\widehat{u}}_{j}-\bm{a}^{\top}\bm{u}_{j}\bm{u}_{j}^{\top}\bm{\widehat{u}}_{j}\bigg{)}\leq
z\bigg{\\}}-\Phi(z)\bigg{|}$ $\displaystyle\lesssim{\sf ErrPCA},$

where

$\displaystyle{\sf ErrPCA}$
$\displaystyle:=\underbrace{\kappa^{3/2}r\log^{3}(n\vee
p)\Bigg{(}\frac{\sigma^{2}}{\lambda_{j}}\bigg{(}\frac{p}{n}+\sqrt{\frac{p}{n}}+\sqrt{\frac{\log(n\vee
p)}{n}}\bigg{)}+\frac{\sigma}{\sqrt{\lambda_{j}}}\sqrt{\frac{p}{n}}\Bigg{)}}_{\text{Subspace
Estimation Effect}}$ $\displaystyle\quad+\underbrace{\kappa
r^{5/2}\log^{5/2}(n\vee
p)\frac{(\lambda_{\max}+\sigma^{2})}{\Delta_{j}\sqrt{n}}}_{\text{Small
Eigengaps Effect}}+\underbrace{\kappa^{2}r^{3/2}\log^{5/2}(n\vee
p)\frac{1}{\sqrt{n}}}_{\text{Parametric Effect}}.$

######  Remark 3.10 (Signal-Strength Condition).

Up to logarithmic terms, factors of $\kappa$, and factors of $r$, our
assumption on $\lambda_{\min}$ in (4) requires that

$\displaystyle\lambda_{\min}/\sigma^{2}\gtrsim\max\bigg{\\{}\frac{p}{n},\sqrt{\frac{p}{n}},\frac{1}{\sigma}\frac{p}{n}\bigg{\\}},$

which matches the assumption in Li et al. (2022). By Cai et al. (2013), the
minimax rate of estimation for $\bm{U}$ in Frobenius norm is of the form

$\displaystyle\frac{pr}{n}\frac{(\lambda_{\min}/\sigma^{2})+1}{(\lambda_{\min}/\sigma^{2})^{2}}\asymp\frac{pr}{n}\max\bigg{\\{}\frac{1}{\lambda_{\min}/\sigma^{2}},\frac{1}{(\lambda_{\min}/\sigma^{2})^{2}}\bigg{\\}}.$
(5)

Therefore, when $\sigma\gtrsim\max\\{\sqrt{\frac{n}{p}},1\\}$, our assumption
is optimal up to logarithmic terms. When
$\sigma\ll\max\\{\sqrt{\frac{n}{p}},1\\}$, our condition is not optimal, but
this regime falls under a setting where $\bm{\Sigma}$ is close to rank
degenerate, which is a simpler setting than when $\sigma$ is larger. This
regime may require a slightly different analysis to take into account the
small amount of noise, but the main results will likely not materially change.

######  Remark 3.11 (Eigengap Conditions).

In order for the bound to be non-vacuous, when $\kappa=O(1)$, our result
requires
$\Delta_{j}\gtrsim(\lambda_{\max}+\sigma^{2})\frac{r^{5/2}\log^{5/2}(n\vee
p)}{\sqrt{n}}$. In contrast, classical theory such as the Davis-Kahan Theorem
implicitly requires that

$\displaystyle\Delta_{j}\asymp(\lambda_{\max}+\sigma^{2}),$

when $r,\kappa=O(1)$. Consequently, when $r,\kappa=O(1)$ and $\log(p)\ll n$,
our assumption is significantly weaker by a factor of
$\widetilde{\Omega}(n^{-1/2})$, while still being amenable to the setting
$\Delta_{j}\asymp(\lambda_{\max}+\sigma^{2})$.

######  Remark 3.12 (Biased Estimation).

As in matrix denoising, the centering term in Theorem 3.3 is given by
$\bm{a}^{\top}\bm{u}_{j}\bm{u}_{j}^{\top}\bm{\widehat{u}}_{j}$. In previous
work (e.g., Paul (2007)) it has been established that
$\bm{\widehat{u}}_{1}^{\top}\bm{u}_{1}\to c\in(0,1)$ when
$\lambda_{\min}/\sigma^{2}\asymp 1$ and $p/n$ converges to a constant
(provided $\lambda_{\min}/\sigma^{2}$ is above a certain threshold), and
$\bm{\widehat{u}}_{1}^{\top}\bm{u}_{1}=1-o(1)$ when
$\lambda_{\min}/\sigma^{2}$ grows sufficiently quickly. Therefore, Theorem 3.3
demonstrates that $\bm{a}^{\top}\bm{\widehat{u}}_{j}$ is a biased estimate of
$\bm{a}^{\top}\bm{u}_{j}$, with bias given by
$\bm{a}^{\top}\bm{u}_{j}(1-\bm{\widehat{u}}_{j}^{\top}\bm{u}_{j})$. In Theorem
3.4 we will study the effect of bias estimation on the distributional theory.

######  Remark 3.13 (Rate of Convergence).

When $\kappa,r=O(1)$, it holds that

$\displaystyle\bigg{|}\mathbb{P}\bigg{\\{}\frac{1}{s_{\bm{a},j}^{{\sf PCA}}}$
$\displaystyle\bigg{(}\bm{a}^{\top}\bm{\widehat{u}}_{j}-\bm{a}^{\top}\bm{u}_{j}\bm{u}_{j}^{\top}\bm{\widehat{u}}_{j}\bigg{)}\leq
z\bigg{\\}}-\Phi(z)\bigg{|}$
$\displaystyle=\widetilde{O}\Bigg{\\{}\frac{(\lambda_{\max}+\sigma^{2})}{\Delta_{j}\sqrt{n}}+\frac{1}{\sqrt{n}}+\frac{\sigma^{2}}{\lambda_{j}}\bigg{(}\frac{p}{n}+\sqrt{\frac{p}{n}}\bigg{)}+\frac{\sigma}{\sqrt{\lambda_{j}}}\sqrt{\frac{p}{n}}\Bigg{\\}}.$

Therefore, the rate becomes

$\displaystyle\max\Bigg{\\{}\frac{(\lambda_{\max}+\sigma^{2})}{\Delta_{j}\sqrt{n}},\frac{1}{\sqrt{n}},\frac{\sigma^{2}}{\lambda_{j}}\sqrt{\frac{p}{n}},\frac{\sigma}{\sqrt{\lambda_{j}}}\sqrt{\frac{p}{n}}\Bigg{\\}},$

ignoring logarithmic terms. Similar to the matrix denoising setting, the
quantities $\frac{\sigma^{2}}{\lambda_{j}}\sqrt{\frac{p}{n}}$ and
$\frac{\sigma}{\lambda_{j}^{1/2}}\sqrt{\frac{p}{n}}$ quantify the “subspace
estimation effect,” or how easy it is to separate the leading $r$ principal
components from the remaining components, where the two different regimes
arise according to the relationship between $\lambda/\sigma^{2}$ and
$\frac{p}{n}$ (see the minimax rate given in (5)), and the quantity
$\frac{\lambda_{\max}+\sigma^{2}}{\Delta_{j}\sqrt{n}}$ quantifies the “small
eigengaps effect” via how well it is to separate the $j$’th eigenvector from
the other eigenvectors. Unlike matrix denoising, PCA exhibits an additional
factor of $\frac{1}{\sqrt{n}}$, representing the classical parametric rate
that cannot be overcome even when the noise is relatively small.

######  Remark 3.14 (Comparison to Previous Works).

Theorem 3.3 is most similar to the works Koltchinskii and Lounici (2016);
Koltchinskii et al. (2020), both of whom study estimating
$\bm{a}^{\top}\bm{u}_{j}$. In fact, in the particular setting that
$\sigma^{2}=o\bigg{(}\lambda_{\min}\frac{n}{p\times\mathrm{polylog}(n)}\bigg{)}$,
that $\lambda_{1}$ through $\lambda_{r}$ converge to some fixed distinct
values, that $r=O(1)$, and that $\Delta_{j}\gtrsim\lambda_{\max}+\sigma^{2}$,
our results resemble Theorem 7 of Koltchinskii and Lounici (2016). However,
our eigengap condition is significantly weaker as mentioned above, and our
results are nonasymptotic, thereby allowing $\kappa$ and $r$ to grow in $n$
and $p$. Furthermore, our asymptotic variance matches the variance obtained in
Koltchinskii et al. (2020), who demonstrate that this variance is optimal over
a sufficiently broad class of covariance matrices $\bm{\Sigma}$. Therefore,
Theorem 3.3 shows that the empirical estimator is _biased_ , but
_asymptotically efficient_.

Our proof techniques are also significantly different from these previous
works. In Koltchinskii and Lounici (2016) and Koltchinskii et al. (2020) the
authors use tools from holomorphic functional calculus to develop their main
results; however, when the eigengaps are permitted to be small these
techniques are no longer available. A key feature of the holomorphic
functional calculus is that it automatically yields the form of the leading-
order term that gives rise to the asymptotic normality. Therefore, in our
results we must identify the leading-order term directly and show that the
remaining terms are sufficiently small relative to the leading-order standard
deviation $s_{\bm{a},j}^{{\sf PCA}}$.

In addition, our assumption on $\Delta_{j}$ is comparable to that in Li et al.
(2022). However, while Li et al. (2022) need only show that the error is of
order $O(s_{\bm{a},j}^{{\sf PCA}})$, we must demonstrate that the residual
terms are of order $o(s_{\bm{a},j}^{{\sf PCA}})$. As in matrix denoising, we
find it useful to identify a certain form of “eigenvalue bias” (see Lemma 8.9)
that allows us to provide sharp bounds on residual terms. This
characterization together with our careful residual expansion allows us to
provide the stronger guarantees required to obtain Theorem 3.3.

######  Remark 3.15 (Cosines).

Assuming that $\kappa,r=O(1)$, take $\bm{a}=\bm{u}_{l}$ for $l\neq j$ with
$l\leq r$ in Theorem 3.3 to reveal that (assuming $\lambda_{l}$ is unique)

$\displaystyle\sqrt{n}\frac{|\lambda_{j}-\lambda_{l}|}{\sqrt{(\lambda_{j}+\sigma^{2})(\lambda_{l}+\sigma^{2})}}\cos\theta(\bm{u}_{l},\bm{\widehat{u}}_{j})\to\mathcal{N}(0,1),$

as long as $\lambda_{\min}$ satisfies

$\displaystyle\lambda_{\min}/\sigma^{2}\gg\widetilde{\Omega}\bigg{(}\max\bigg{\\{}\frac{p}{n},\sqrt{\frac{p}{n}},\frac{1}{\sigma}\frac{p}{n}\bigg{\\}}\bigg{)}$

and $\Delta_{j}$ satisfies

$\displaystyle\Delta_{j}\gg(\lambda_{\max}+\sigma^{2})\frac{1}{\sqrt{n}}.$

To the best of our knowledge this is the first result establishing the
asymptotic Gaussianity of the cosines of empirical principal components in the
small eigengaps regime.

######  Example 4 (“Diagonal” versus “Dense” Covariance, Continued).

Suppose that $\bm{\Sigma}$ is the diagonal covariance matrix in Example 2.
Then $(s_{\bm{a},1}^{{\sf PCA}})^{2}$ simplifies to

$\displaystyle(s_{\bm{a},1}^{{\sf PCA}})^{2}$
$\displaystyle=\frac{(\lambda+\varepsilon+\sigma^{2})(\lambda+\sigma^{2})|\bm{a}(2)|^{2}}{n\varepsilon^{2}}+\frac{(\lambda+\varepsilon+\sigma^{2})\sigma^{2}\|\bm{a}_{-1-2}\|^{2}}{n(\lambda+\varepsilon)^{2}},$

hwere $\bm{a}_{-1-2}$ denotes the vector $\bm{a}$ with its first two
coordinates removed. The expression above reveals how the “spread” of the
entries of $\bm{a}$ affects the asymptotic variance, where less concentration
of the magnitude of $\bm{a}$ in the first two entries corresponds to a smaller
asymptotic variance. Furthermore, as $\varepsilon\to 0$ with $\bm{a}$ held
fixed, the asymptotic variance increases.

Now suppose that $\bm{\Sigma}$ is the “dense” covariance, and suppose that
$\bm{a}$ has its first $p/2$ entries nonnegative and its second $p/2$ entries
nonpositive. Then

$\displaystyle(s_{\bm{a},1}^{{\sf PCA}})^{2}$
$\displaystyle=\frac{(\lambda+\varepsilon+\sigma^{2})(\lambda+\sigma^{2})\|\bm{a}\|_{1}^{2}}{np\varepsilon^{2}}+\frac{(\lambda+\varepsilon+\sigma^{2})\sigma^{2}\|\bm{U}_{\perp}^{\top}\bm{a}\|^{2}}{n(\lambda+\varepsilon)^{2}}.$

If $\|\bm{a}\|_{1}\asymp 1$ (e.g., if $\bm{a}$ is sufficiently sparse), then
the second term can dominate the first term for $p$ sufficiently large.
However, if $\|\bm{a}\|_{1}\asymp\sqrt{p}$ (e.g., if $\bm{a}$ is sufficiently
dense with the magnitude of its entries approximately equally spread), then
the first term can dominate.

These two examples showcase how manifestly different phenomena result from
both the unit vector $\bm{a}$ and the underlying covariance matrix. Such a
result is reminiscent of the setting of statistical inference for linear
functions of high-dimensional coefficients in linear regression, where
different rates appear based on the sparsity of the linear function (Cai and
Guo, 2017).

We next consider the approximate Gaussianity of the debiased estimator
proposed in Li et al. (2022), assuming knowledge of the noise variance
$\sigma^{2}$.

######  Theorem 3.4 (Approximate Gaussianity of Debiased Estimator – PCA).

Suppose the conditions of Theorem 3.3 hold. Define the quantity $b_{j}^{{\sf
PCA}}$ via

$\displaystyle b_{j}^{{\sf PCA}}$
$\displaystyle:=\begin{cases}\dfrac{\widehat{\lambda}_{j}}{n+\sum_{r<i\leq
n}\frac{\widehat{\lambda}_{i}}{\widehat{\lambda}_{j}-\widehat{\lambda}_{i}}}{\sum}_{i:r<i\leq
n}\dfrac{\widehat{\lambda}_{i}}{(\widehat{\lambda}_{j}-\widehat{\lambda}_{i})^{2}}&n\geq
p;\\\
\dfrac{\frac{\sigma^{2}p}{n}}{\widehat{\lambda}_{j}-\frac{\sigma^{2}p}{n}}+\dfrac{\widehat{\lambda}_{j}}{\widehat{\lambda}_{j}-\frac{\sigma^{2}p}{n}}\dfrac{\widehat{\lambda}_{j}}{n+\sum_{i:r<i\leq
n}\frac{\widehat{\lambda}_{i}}{\widehat{\lambda}_{j}-\widehat{\lambda}_{i}}}{\sum}_{i:r<i\leq
n}\dfrac{\widehat{\lambda}_{i}-\frac{\sigma^{2}p}{n}}{(\widehat{\lambda}_{j}-\widehat{\lambda}_{i})^{2}}&n<p.\end{cases}$

Then it holds that

$\displaystyle\sup_{z\in\mathbb{R}}\bigg{|}\mathbb{P}\bigg{\\{}\frac{1}{s_{\bm{a},j}^{{\sf
PCA}}}\bigg{(}\bm{a}^{\top}\bm{\widehat{u}}_{j}\sqrt{1+b_{j}^{{\sf
PCA}}}-\bm{a}^{\top}\bm{u}_{j}\bigg{)}\leq z\bigg{\\}}-\Phi(z)\bigg{|}$
$\displaystyle\lesssim{\sf
ErrPCA}+\frac{|\bm{a}^{\top}\bm{u}_{j}|}{s_{\bm{a},j}^{{\sf PCA}}}{\sf
ErrBiasPCA},$

where ${\sf ErrPCA}$ and $s^{{\sf PCA}}_{\bm{a},j}$ are defined in Theorem
3.3, and ${\sf ErrBiasPCA}$ is defined via

$\displaystyle{\sf ErrBiasPCA}$
$\displaystyle:=\frac{(\lambda_{\max}+\sigma^{2})(\lambda_{j}+\sigma^{2})r\log(n\vee
p)}{\Delta_{j}^{2}n}$
$\displaystyle\quad+\frac{\sigma(\lambda_{j}+\sigma^{2})^{1/2}}{\lambda_{j}\sqrt{n}}\Bigg{(}\frac{\sigma\kappa\sqrt{r}\log^{2}(n\vee
p)}{\lambda_{\min}^{1/2}}\sqrt{\frac{p}{n}}+\frac{\sigma^{2}\sqrt{r\kappa}\log^{2}(n\vee
p)}{\lambda_{\min}}\sqrt{\frac{p}{n}}\Bigg{)}$

######  Remark 3.16 (Bias Effect).

As in matrix denoising, Theorem 3.4 incurs an additional “debiasing effect”
pertaining to the magnitude of $|\bm{a}^{\top}\bm{u}_{j}|$, whereas Theorem
3.3 does not depend on $|\bm{a}^{\top}\bm{u}_{j}|$ except through the
centering term. It is possible to demonstrate that

$\displaystyle\frac{|\bm{a}^{\top}\bm{u}_{j}|}{s^{{\sf PCA}}_{\bm{a},j}}{\sf
ErrBiasPCA}$
$\displaystyle\lesssim\frac{(\lambda_{\max}+\sigma^{2})\sqrt{r}\kappa\log(n\vee
p)}{\Delta_{j}\sqrt{n}}\min_{k\neq j,k\leq
r}\frac{|\bm{a}^{\top}\bm{u}_{j}|}{|\bm{a}^{\top}\bm{u}_{k}|}\frac{|\lambda_{k}-\lambda_{j}|}{\Delta_{j}}$
$\displaystyle\quad+\frac{|\bm{a}^{\top}\bm{u}_{j}|}{\|\bm{U}_{\perp}^{\top}\bm{a}\|}\Bigg{(}\frac{\sigma\kappa\sqrt{r}\log^{2}(n\vee
p)}{\lambda_{\min}^{1/2}}\sqrt{\frac{p}{n}}+\frac{\sigma^{2}\sqrt{r\kappa}\log^{2}(n\vee
p)}{\lambda_{\min}}\sqrt{\frac{p}{n}}\Bigg{)},$

where the term in parentheses is $o(1)$ under our noise assumption (4). This
additional quantity reflects the effect of $\bm{a}$ on the subspace
corresponding to $\bm{u}_{j}$ relative to the other subspaces.

Finally we consider the entrywise distribution of individual eigenvectors. The
following result provides an analogue of Corollary 3.1 for PCA. The proof can
be found in Section C.8.

######  Corollary 3.2 (Entrywise Distributional Theory).

In the setting of Theorem 3.2, let $\bm{a}=\bm{e}_{i}$; i.e. $\bm{a}$ is a
standard basis vector, and suppose there exists a constant $\mu=O(1)$ such
that

$\displaystyle\max_{k\leq r}|\bm{u}_{k}(i)|\leq\mu\sqrt{\frac{1}{p}}.$

Suppose further that

$\displaystyle\max_{k\neq j,k\leq
r}\frac{|\bm{u}_{j}(i)|}{|\bm{u}_{k}(i)|}=O(1).$

Let $s^{{\sf PCA}}_{\bm{e}_{i},j}$ be as in Theorem 3.3. Then it holds that

$\displaystyle\sup_{z\in\mathbb{R}}\bigg{|}\mathbb{P}\bigg{\\{}\frac{1}{s^{{\sf
PCA}}_{\bm{e}_{i},j}}\bigg{(}\bm{\widehat{u}}_{j}(i)-\bm{u}_{j}(i)\bigg{)}\leq
z\bigg{\\}}-\Phi(z)\bigg{|}\lesssim{\sf ErrPCA},$

where ${\sf ErrPCA}$ is defined in Theorem 3.3.

Just as in the case of matrix denoising, the form of
$\bm{s}_{\bm{e}_{i},j}^{{\sf PCA}}$ characterizes the effect of other
eigenvectors and eigenvalues on the asymptotic variance, whereas previous
works (e.g., Yan et al. (2021)) only hold up to some nonidentifiable
orthogonal transformation.

Note that Corollary 4.2 holds for the “dense” covariance model discussed in
Example 2, but it does not hold for the “diagonal” covariance model.

##  4 Statistical Inference

While the results in Section 3 yield approximate Gaussianity for the bias-
corrected estimator for $\bm{a}^{\top}\bm{u}_{j}$, they fall short of
providing the guarantees required to quantify the uncertainty in the estimate.
In this section we consider estimation and inference procedures that result in
approximately valid coverage.

###  4.1 Matrix Denoising

Recall that Theorem 3.2 requires knowledge of the noise variance $\sigma^{2}$
which is typically unknown in practice. We therefore define the following
estimator of $\widehat{\sigma}^{2}$, defined via

$\displaystyle\widehat{\sigma}^{2}:$
$\displaystyle=\binom{n}{2}^{-1}\|\mathcal{P}_{{\sf upper-
diag}}\big{(}\bm{\widehat{S}}-\bm{\widehat{S}}_{r}\big{)}\|_{F}^{2},$ (6)

where $\bm{\widehat{S}}_{r}$ is the best rank $r$ approximation of
$\bm{\widehat{S}}$, and where $\mathcal{P}_{{\sf upper-diag}}$ denotes the
projection onto the upper off-diagonal. Using this estimator of the noise
variance we can reliably estimate $b_{j}^{{\sf MD}}$ with
$\widehat{b_{j}^{{\sf MD}}}$, defined via

$\displaystyle\widehat{b_{j}^{{\sf
MD}}}:=\sum_{k>r}\frac{\widehat{\sigma}^{2}}{(\widehat{\lambda}_{j}-\widehat{\lambda}_{k})^{2}}.$
(7)

Next, when the eigenvalues are distinct, recall that $s^{{\sf MD}}_{\bm{a},j}$
is defined via

$\displaystyle\big{(}s^{{\sf MD}}_{\bm{a},j}\big{)}^{2}=2\sum_{k\neq j,k\leq
r}\frac{\sigma^{2}}{(\lambda_{k}-\lambda_{j})^{2}}(\bm{u}_{k}^{\top}\bm{a})^{2}+2\frac{\sigma^{2}}{\lambda_{j}^{2}}\|\bm{U}_{\perp}^{\top}\bm{a}\|^{2}.$

There are several quantities that require estimating; the first is
$\sigma^{2}$, for which we use $\widehat{\sigma}^{2}$ above. Another key
quantity that requires estimating is $(\bm{u}_{k}^{\top}\bm{a})^{2}$ for
$k\neq j$; however, by the theory in the previous section the plug-in
estimator $\bm{\widehat{u}}_{k}^{\top}\bm{a}$ will be biased for
$\bm{u}_{k}^{\top}\bm{a}$. Therefore, instead of using
$(\bm{\widehat{u}}_{k}^{\top}\bm{a})^{2}$ for our estimate of
$(\bm{u}_{k}^{\top}\bm{a})^{2}$, we propose to use
$(\bm{\widehat{u}}_{k}^{\top}\bm{a})^{2}(1+\widehat{b_{k}^{{\sf MD}}})$, where
$\widehat{b_{k}^{{\sf MD}}}$ is the estimate of the debiasing parameter for
the $k$’th eigenvector obtained via Eq. 7.

Finally, $s^{{\sf MD}}_{\bm{a},j}$ also requires estimating the leading $r$
eigenvalues. Given $k\leq r$, a natural estimate for $\lambda_{k}$ is the
empirical eigenvalue $\widehat{\lambda}_{k}$; however, our analysis (see Lemma
8.2) demonstrates that $\widehat{\lambda}_{k}$ is a biased estimate of
$\lambda_{k}$. Therefore, motivated by the form of this eigenvalue bias, we
propose to use a debiased estimate of $\lambda_{k}$, defined via

$\displaystyle\widecheck{\lambda}_{k}:=\widehat{\lambda}_{k}-\sum_{i>r}\frac{\widehat{\sigma}^{2}}{\widehat{\lambda}_{k}-\widehat{\lambda}_{i}}.$

With all these individual estimators in place, we then define

$\displaystyle\big{(}\widehat{s_{\bm{a},j}^{{\sf MD}}}\big{)}^{2}$
$\displaystyle=2\sum_{k\neq j,k\leq
r}\frac{\widehat{\sigma}^{2}}{(\widecheck{\lambda}_{j}-\widecheck{\lambda}_{k})^{2}}(\bm{\widehat{u}}_{k}^{\top}\bm{a})^{2}(1+\widehat{b_{k}^{{\sf
MD}}})+2\frac{\widehat{\sigma}^{2}}{\widecheck{\lambda}_{j}^{2}}\|\bm{\widehat{U}}_{\perp}^{\top}\bm{a}\|^{2}.$

This estimator is entirely data-driven and does not require any sample-
splitting. With this estimator we can define a $1-\alpha$ confidence interval
for $\bm{a}^{\top}\bm{u}_{j}$; the full procedure is summarized in Algorithm
1. The following result demonstrates the approximate validity of the resulting
confidence intervals.

0: Input matrix $\bm{\widehat{S}}$, rank $r$, index $j\in[r]$, unit vector
$\bm{a}$, desired confidence level $\alpha$.

1: Compute the eigendecomposition of
$\bm{\widehat{S}}=\big{[}\bm{\widehat{U}},\bm{\widehat{U}}_{\perp}\big{]}\bm{\widehat{\Lambda}}\big{[}\bm{\widehat{U}},\bm{\widehat{U}}_{\perp}\big{]}^{\top}$,
where $\bm{\widehat{U}}$ has columns consisting of the $r$ eigenvectors
corresponding to the largest in magnitude nonzero eigenvalues. Let
$\bm{\widehat{u}}_{k}$ denote the $k$’th column of $\bm{\widehat{U}}$ and
$\widehat{\lambda}_{k}$ denote the $k$’th largest in magnitude eigenvalue of
$\bm{\widehat{S}}$.

2: Set $\widehat{\sigma}^{2}$ via Eq. 6.

3: Set $\widehat{b_{k}^{{\sf MD}}}$ via Eq. 7 for $1\leq k\leq r$.

4: Compute, for $1\leq k\leq r$,
$\displaystyle\widecheck{\lambda}_{k}:=\widehat{\lambda}_{k}-\sum_{i>r}\frac{\widehat{\sigma}^{2}}{\widehat{\lambda}_{k}-\widehat{\lambda}_{i}}.$

5: Set  $\displaystyle\big{(}\widehat{s_{\bm{a},j}^{{\sf
MD}}}\big{)}^{2}:=2\sum_{k\leq r,k\neq
j}\frac{\widehat{\sigma}^{2}}{(\widecheck{\lambda}_{j}-\widecheck{\lambda}_{k})^{2}}(\bm{\widehat{u}}_{k}^{\top}\bm{a})^{2}(1+\widehat{b_{k}^{{\sf
MD}}})+2\frac{\widehat{\sigma}^{2}}{\widecheck{\lambda}_{j}^{2}}\|\bm{\widehat{U}}_{\perp}^{\top}\bm{a}\|^{2}.$

6: Compute ${\sf C.I.}^{(\alpha)}(\bm{\widehat{u}}_{j}^{\top}\bm{a})$ defined
via  $\displaystyle{\sf
C.I.}^{(\alpha)}(\bm{\widehat{u}}_{j}^{\top}\bm{a}):=\bigg{[}\bm{\widehat{u}}_{j}^{\top}\bm{a}\sqrt{1+\widehat{b_{j}^{{\sf
MD}}}}\pm\Phi^{-1}\big{(}1-\alpha/2\big{)}\widehat{s_{\bm{a},j}^{{\sf
MD}}}\bigg{]},$ where $\Phi^{-1}$ is the inverse of the cumulative density
function for a standard Gaussian random variable.

7: return Estimated confidence interval ${\sf
C.I.}^{(\alpha)}(\bm{\widehat{u}}_{j},\bm{a})$.

Algorithm 1 Confidence Interval for $\bm{a}^{\top}\bm{u}_{j}$ – Matrix
Denoising

######  Theorem 4.1 (Validity of Confidence Intervals for Matrix Denoising).

Consider the model (1), and suppose $\lambda_{1}$ through $\lambda_{r}$ are
unique. Suppose further that

$\displaystyle\lambda_{\min}/\sigma\geq
C_{0}\kappa^{2}r\sqrt{n}\log^{3/2}(n);\qquad\Delta_{\min}/\sigma\geq
C_{0}r^{3/2}\log^{5/2}(n)$

for some sufficiently large constant $C_{0}$. In addition, assume that $r\leq
c_{0}\frac{n}{\kappa^{2}\log^{2}(n)}$, where $c_{0}$ is a sufficiently small
constant. Suppose further that

$\displaystyle\max_{k\neq j,k\leq r}\frac{s_{\bm{a},k}^{{\sf
MD}}}{s_{\bm{a},j}^{{\sf MD}}}\frac{\sigma\sqrt{r}\log(n)}{\Delta_{j}}=o(1).$

Furthermore, assume that

$\displaystyle\frac{|\bm{a}^{\top}\bm{u}_{j}|}{s_{\bm{a},j}^{{\sf
MD}}}\bigg{(}\frac{\sigma^{2}r\log(n)}{\Delta_{j}^{2}}+\frac{\sigma^{2}\kappa\sqrt{rn\log(n)}}{\lambda_{\min}^{2}}\bigg{)}=o(1).$

Let ${\sf C.I.}^{(\alpha)}(\bm{\widehat{u}}_{j}^{\top}\bm{a})$ denote the
output of Algorithm 1. If either
$\|\bm{U}^{\top}\bm{a}\|>\|\bm{U}_{\perp}^{\top}\bm{a}\|$ with
$|\bm{a}^{\top}\bm{u}_{j}|\lesssim|\bm{a}^{\top}\bm{u}_{k}|$ for $k\neq
j,k\leq r$ or $\|\bm{U}_{\perp}^{\top}\bm{a}\|\geq\|\bm{U}^{\top}\bm{a}\|$, it
holds that

$\displaystyle\bigg{|}\mathbb{P}\bigg{\\{}\bm{u}_{j}^{\top}\bm{a}\in{\sf
C.I.}^{(\alpha)}(\bm{\widehat{u}}_{j}^{\top}\bm{a})\bigg{\\}}-(1-\alpha)\bigg{|}=o(1).$

To the best of our knowledge, Theorem 4.1 is the first result in the
literature providing asymptotically valid confidence intervals for linear
functions of eigenvectors in the matrix denoising context considered herein.
Moreover, our results allow $r$ and $\kappa$ to grow with $n$.

######  Remark 4.1 (Generality).

En route to the proof of Theorem 4.1 we actually prove a significantly more
general, albeit more complicated, result in the form of a Berry-Esseen bound
similar to Theorem 3.2. In particular, we show that

$\displaystyle\sup_{z\in\mathbb{R}}\bigg{|}\mathbb{P}\bigg{\\{}\frac{\bm{a}^{\top}\bm{\widehat{u}}_{j}\sqrt{1+\widehat{b_{j}^{{\sf
MD}}}}-\bm{a}^{\top}\bm{u}_{j}}{\widehat{s_{\bm{a},j}^{{\sf MD}}}}\leq
z\bigg{\\}}-\Phi(z)\bigg{|}\leq\zeta,$

where $\zeta$ is the quantity in Lemma 8.6. In Lemma 8.6 we demonstrate that
the conditions of Theorem 4.1 imply that $\zeta=o(1)$.

######  Remark 4.2 (Alignment of $\bm{a}$).

Observe that Theorem 4.1 requires different assumptions based on how well
$\bm{a}$ is aligned with the leading eigenspace: if
$\|\bm{U}_{\perp}^{\top}\bm{a}\|<\|\bm{U}^{\top}\bm{a}\|$, we further assume
that $|\bm{a}^{\top}\bm{u}_{j}|\lesssim|\bm{a}^{\top}\bm{u}_{k}|$, whereas if
$\|\bm{U}_{\perp}^{\top}\bm{a}\|\geq\|\bm{U}^{\top}\bm{a}\|$ we make no
further assumptions. To understand the intuition behind this condition, when
$\|\bm{U}_{\perp}^{\top}\bm{a}\|<\|\bm{U}^{\top}\bm{a}\|$, the vector $\bm{a}$
is more closely aligned with the leading eigenspace, and hence we must ensure
that $\bm{a}$ does not fall “too much” into the subspace spanned by
$\bm{u}_{j}$ (relative to $\bm{u}_{k}$ for $k\leq r$). When
$\|\bm{U}_{\perp}^{\top}\bm{a}\|\geq\|\bm{U}^{\top}\bm{a}\|$, the vector
$\bm{a}$ lies primarily in the subspace spanned by $\bm{U}_{\perp}$, and we
need not make any further assumptions (beyond the assumptions already made).

######  Remark 4.3 (Eigengap Conditions).

Note that Theorem 4.1 requires a condition on the minimum eigengap
$\Delta_{\min}$, wheras Theorems 3.1 and 3.2 only impose conditions on
$\Delta_{j}$. In essence, this extra assumption is imposed to ensure a
faithful estimate of $s_{\bm{a},j}^{{\sf MD}}$, which requires estimating
$(\bm{a}^{\top}\bm{u}_{k})^{2}$ for all $k\leq r$ with $k\neq j$, and hence
our estimation of these quantities will only be uniformly guaranteed when
$\Delta_{\min}$ is sufficiently large. Nonetheless, this condition on
$\Delta_{\min}$ is still significantly weaker than required by classical
theory. It may be possible to extend the result to settings with repeated
eigenvalues assuming _a priori_ knowledge of which eigenvalues are repeated.

######  Example 5 (Gaussian Mixed-Membership Blockmodel, Continued).

In the setting of Example 1 again consider the case
$\varepsilon\asymp\sigma\mathrm{polylog}(n)$ and
$\lambda\asymp\sigma\sqrt{n}\mathrm{polylog}(n)$. Then to test whether
$\alpha=0$ one may form a confidence interval using Algorithm 1 with $\bm{a}$
equal to the constant vector and rejecting if $0$ is not in the interval.
Under the null hypothesis $\alpha=0$, this test is consistent, and under local
alternatives such that $\frac{|\bm{a}^{\top}\bm{u}_{j}|}{s_{\bm{a},j}^{{\sf
MD}}}\gg 1$, this test will have power converging to one. Further analysis is
needed for when such local alternatives hold, but we leave this to future
work.

Finally we return to the entrywise setting considered in Corollary 3.1, where
we recall that no debiasing parameter is needed. The following result
demonstrates the validity of entrywise confidence intervals obtained without
the extra debiasing step under a slightly different condition than given in
Theorem 4.1. The proof of this result can be found in Section B.8.

######  Corollary 4.1 (Entrywise Confidence Intervals – Matrix Denoising).

Consider the model (1) and suppose $\lambda_{1}$ through $\lambda_{r}$ are
unique. Suppose further that

$\displaystyle\lambda_{\min}\geq C_{0}\sigma
r\sqrt{n}\kappa^{2}\log^{3/2}(n);\qquad\Delta_{\min}\geq C_{0}\sigma
r^{3/2}\kappa\log^{5/2}(n)$

for some sufficiently large constant $C_{0}$. In addition, assume that $r\leq
c_{0}\frac{n}{\kappa^{2}\log^{2}(n)}$, where $c_{0}$ is a sufficiently small
constant. Assume that there exists a constant $\mu=O(1)$ such that

$\displaystyle\max_{k\leq r}|\bm{u}_{k}(i)|\leq\mu\sqrt{\frac{1}{n}}.$

Suppose further that

$\max_{k,l\leq r}\frac{|\bm{u}_{l}(i)|}{|\bm{u}_{k}(i)|}=O(1).$

Let ${\sf C.I.}^{(\alpha)}(\bm{\widehat{u}}_{j}(i))$ denote the output of
Algorithm 2. Then it holds that

$\displaystyle\bigg{|}\mathbb{P}\bigg{\\{}\bm{u}_{j}(i)\in{\sf
C.I.}^{(\alpha)}(\bm{\widehat{u}}_{j}(i))\bigg{\\}}-(1-\alpha)\bigg{|}=o(1).$

In essence, Corollary 4.1 delivers valid entrywise confidence intervals for
empirical eigenvectors when the magnitude of the entries of each of the
leading eigenvectors are of similar order. Existing work has primarily
delivered _row-wise_ confidence intervals for eigenspaces, and hence our
result can be understood as a more “fine-grained” result providing confidence
intervals for entries of individual eigenvectors while still allowing for the
setting where eigengaps can be small. While we assume that $\mu=O(1)$ it is
possible that the analysis can be extended to the setting $\mu=o(\sqrt{n})$,
but we leave this for future work.

0: Input matrix $\bm{\widehat{S}}$, rank $r$, index $j\in[r]$, index
$i\in[n]$, desired confidence level $\alpha$.

1: Compute the eigendecomposition of
$\bm{\widehat{S}}=\big{[}\bm{\widehat{U}},\bm{\widehat{U}}_{\perp}\big{]}\bm{\widehat{\Lambda}}\big{[}\bm{\widehat{U}},\bm{\widehat{U}}_{\perp}\big{]}^{\top}$,
where $\bm{\widehat{U}}$ has columns consisting of the $r$ eigenvectors
corresponding to the largest in magnitude nonzero eigenvalues. Let
$\bm{\widehat{u}}_{k}$ denote the $k$’th column of $\bm{\widehat{U}}$ and
$\widehat{\lambda}_{k}$ denote the $k$’th largest in magnitude eigenvalue of
$\bm{\widehat{S}}$.

2: Set $\widehat{\sigma}^{2}$ via Eq. 6.

3: Set $\widehat{b_{k}^{{\sf MD}}}$ via Eq. 7 for $1\leq k\leq r$, $k\neq j$.

4: Compute, for $1\leq k\leq r$,
$\displaystyle\widecheck{\lambda}_{k}:=\widehat{\lambda}_{k}-\sum_{i>r}\frac{\widehat{\sigma}^{2}}{\widehat{\lambda}_{k}-\widehat{\lambda}_{i}}.$

5: Set  $\displaystyle\big{(}\widehat{s_{\bm{e}_{i},j}^{{\sf
MD}}}\big{)}^{2}:=2\sum_{k\leq r,k\neq
j}\frac{\widehat{\sigma}^{2}}{(\widecheck{\lambda}_{j}-\widecheck{\lambda}_{k})^{2}}(\bm{\widehat{u}}_{k}(i))^{2}(1+\widehat{b_{k}^{{\sf
MD}}})+2\frac{\widehat{\sigma}^{2}}{\widecheck{\lambda}_{j}^{2}}\|\bm{\widehat{U}}_{\perp}^{\top}\bm{e}_{i}\|^{2}.$

6: Compute ${\sf C.I.}^{(\alpha)}(\bm{\widehat{u}}_{j}(i))$ defined via
$\displaystyle{\sf
C.I.}^{(\alpha)}(\bm{\widehat{u}}_{j}(i)):=\bigg{[}\bm{\widehat{u}}_{j}(i)\pm\Phi^{-1}\big{(}1-\alpha/2\big{)}\widehat{s_{\bm{e}_{i},j}^{{\sf
MD}}}\bigg{]},$ where $\Phi^{-1}$ is the inverse of the cumulative density
function for a standard Gaussian random variable.

7: return Estimated confidence interval ${\sf
C.I.}^{(\alpha)}(\bm{\widehat{u}}_{j}(i))$.

Algorithm 2 Entrywise Confidence Intervals – Matrix Denoising

###  4.2 Principal Component Analysis

Similar to the matrix denoising setting, in Theorem 3.4 we assume that the
noise variance $\sigma^{2}$ is assumed known, which is not true in practice.
We therefore consider two different estimators based on the aspect ratio
$\frac{p}{n}$. Define

$\displaystyle\widehat{\sigma}^{2}:=\begin{cases}\dfrac{{\sf
Tr}\big{(}\bm{\widehat{\Sigma}}-\bm{\widehat{\Sigma}}_{r}\big{)}}{p-r}&p\geq\dfrac{n}{\log^{4}(n\vee
p)};\\\ \widehat{\lambda}_{r+1}&p<\dfrac{n}{\log^{4}(n\vee p)}.\end{cases}$
(8)

The intuition of this estimator can be understood as follows: when $p\ll n$,
the eigenvalues of $\bm{\widehat{\Sigma}}$ are already consistent for those of
$\bm{\Sigma}$, and hence $\widehat{\lambda}_{r+1}$ will approximate
$\sigma^{2}$. However, when $p$ is relatively large, the approximation fails,
and a more careful approximation must be done by leveraging the independence
of the noise, which becomes more effective when $p$ is large. The choice of
the cutoff $n/\log^{4}(n\vee p)$ is primarily theoretical; it is possible that
a similar argument will work with any choice for $p$ strictly smaller than
$n$.

With the estimator $\widehat{\sigma}^{2}$ defined, we can also define the
estimated debiasing parameter $\widehat{b_{j}^{{\sf PCA}}}$,defined via

$\displaystyle\widehat{b_{j}^{{\sf
PCA}}}:=\begin{cases}\dfrac{\widehat{\lambda}_{j}}{n+\sum_{r<i\leq
n}\frac{\widehat{\lambda}_{i}}{\widehat{\lambda}_{j}-\widehat{\lambda}_{i}}}{\sum}_{i:r<i\leq
n}\dfrac{\widehat{\lambda}_{i}}{(\widehat{\lambda}_{j}-\widehat{\lambda}_{i})^{2}}&n\geq
p;\\\
\dfrac{\frac{\widehat{\sigma}^{2}p}{n}}{\widehat{\lambda}_{j}-\frac{\widehat{\sigma}^{2}p}{n}}+\dfrac{\widehat{\lambda}_{j}}{\widehat{\lambda}_{j}-\frac{\widehat{\sigma}^{2}p}{n}}\dfrac{\widehat{\lambda}_{j}}{n+\sum_{i:r<i\leq
n}\frac{\widehat{\lambda}_{i}}{\widehat{\lambda}_{j}-\widehat{\lambda}_{i}}}{\sum}_{i:r<i\leq
n}\dfrac{\widehat{\lambda}_{i}-\frac{\widehat{\sigma}^{2}p}{n}}{(\widehat{\lambda}_{j}-\widehat{\lambda}_{i})^{2}}&n<p.\end{cases}$
(9)

Next, recall that

$\displaystyle(s_{\bm{a},j}^{{\sf PCA}})^{2}$ $\displaystyle=\sum_{k\neq
j,k\leq
r}\frac{(\lambda_{j}+\sigma^{2})(\lambda_{k}+\sigma^{2})(\bm{u}_{k}^{\top}\bm{a})^{2}}{n(\lambda_{j}-\lambda_{k})^{2}}+\frac{(\lambda_{j}+\sigma^{2})\sigma^{2}\|\bm{U}_{\perp}^{\top}\bm{a}\|^{2}}{n\lambda_{j}^{2}}.$

As in the matrix denoising setting, estimating the asymptotic variance
requires estimating $(\bm{a}^{\top}\bm{u}_{k})$ for $k\leq r$ and
$\lambda_{k}$ for $k\leq r$. As will be made clear in our analysis (and also
observed in Li et al. (2022)), the empirical eigenvalues
$\widehat{\lambda}_{j}$ exhibit a multiplicative bias (see Lemma 8.9).
Therefore, we propose to use a debiased estimator for
$\lambda_{k}+\sigma^{2}$, defined via:

$\displaystyle\widecheck{\lambda}_{k}:$
$\displaystyle=\frac{\widehat{\lambda}_{k}}{1+\widehat{\gamma}^{{\sf
PCA}}(\widehat{\lambda}_{k})},$

where

$\displaystyle\widehat{\gamma}^{{\sf PCA}}(\widehat{\lambda}_{k})$
$\displaystyle=\frac{1}{n}\sum_{i=r+1}^{\min(p-r,n)}\frac{\widehat{\lambda}_{i}}{\widehat{\lambda}_{k}-\widehat{\lambda}_{i}}.$

Under our assumptions we can demonstrate that $\widecheck{\lambda}_{k}$ is a
strong estimator of $\lambda_{k}+\sigma^{2}$. To estimate $\lambda_{j}^{2}$,
we simply use the estimator $\widecheck{\lambda}_{j}-\widehat{\sigma}^{2}$. If
$\widehat{\sigma}^{2}\approx\sigma^{2}$, then it holds that
$\widecheck{\lambda}_{j}-\widehat{\sigma}^{2}\approx\lambda_{j}$.

Finally, in the asymptotic variance we also need to estimate
$(\bm{a}^{\top}\bm{u}_{k})^{2}$ for $k\leq r$ and $k\neq j$. As in the matrix
denoising case, the plug-in estimators for $\bm{a}^{\top}\bm{u}_{k}$ may not
be sufficiently close to merit a strong estimate. We therefore propose to use
the debiasing parameter $\widehat{b_{k}^{{\sf PCA}}}$ to form the estimator
$(\bm{a}^{\top}\bm{\widehat{u}}_{k})^{2}(1+\widehat{b_{k}^{{\sf PCA}}})$,
which, from the previous analysis, is a strong estimator of
$(\bm{a}^{\top}\bm{u}_{k})^{2}$ under conditions on $\Delta_{\min}$.

Combining these estimators, we set

$\displaystyle\big{(}\widehat{s^{{\sf PCA}}_{\bm{a},j}}\big{)}^{2}:$
$\displaystyle=\sum_{k\neq
j}\frac{\widecheck{\lambda}_{k}\widecheck{\lambda}_{j}(\bm{a}^{\top}\bm{\widehat{u}}_{k})^{2}\big{(}1+\widehat{b_{k}^{{\sf
PCA}}}\big{)}}{(\widecheck{\lambda}_{j}-\widecheck{\lambda}_{k})^{2}}+\frac{\widehat{\sigma}^{2}(\widecheck{\lambda}_{j})\|\bm{\widehat{U}}_{\perp}^{\top}\bm{a}\|^{2}}{(\widecheck{\lambda}_{j}-\widehat{\sigma}^{2})^{2}},$

and we use this estimator to construct a $1-\alpha$ confidence interval for
$\bm{a}^{\top}\bm{u}_{j}$. The full procedure is summarized in Algorithm 3.
The following result demonstrates the approximate validity of the resulting
confidence intervals.

0: Sample data $\bm{X}\in\mathbb{R}^{p\times n}$, rank $r$, index $j\in[r]$,
unit vector $\bm{a}\in\mathbb{R}^{p}$, desired confidence level $\alpha$.

1: Compute the eigendecomposition of
$\bm{\widehat{\Sigma}}=\frac{\bm{XX}^{\top}}{n}=\big{[}\bm{\widehat{U}},\bm{\widehat{U}}_{\perp}\big{]}\bm{\widehat{\Lambda}}\big{[}\bm{\widehat{U}},\bm{\widehat{U}}_{\perp}\big{]}^{\top}$,
where $\bm{\widehat{U}}$ has columns consisting of the $r$ eigenvectors
corresponding to the largest nonzero eigenvalues. Let $\bm{\widehat{u}}_{k}$
denote the $k$’th column of $\bm{\widehat{U}}$ and $\widehat{\lambda}_{k}$
denote the $k$’th largest eigenvalue of $\bm{\widehat{\Sigma}}$.

2: Set $\widehat{\sigma}^{2}$ via Eq. 8.

3: Set $\widehat{b_{k}^{{\sf PCA}}}$ via Eq. 9 for $1\leq k\leq r$.

4: Compute, for $1\leq k\leq r$,  $\displaystyle\widehat{\gamma}^{{\sf
PCA}}(\widehat{\lambda}_{k}):=\frac{1}{n}\sum_{i=r+1}^{\min(p-r,n)}\frac{\widehat{\lambda}_{i}}{\widehat{\lambda}_{k}-\widehat{\lambda}_{i}}$

5: Compute, for $1\leq k\leq r$,
$\displaystyle\widecheck{\lambda}_{k}:=\frac{\widehat{\lambda}_{k}}{1+\widehat{\gamma}^{{\sf
PCA}}(\widehat{\lambda}_{k})}.$

6: Set  $\displaystyle\big{(}\widehat{s_{\bm{a},j}^{{\sf
PCA}}}\big{)}^{2}:=\sum_{k\leq r,k\neq
j}\frac{\widecheck{\lambda}_{k}\widecheck{\lambda}_{j}(\bm{a}^{\top}\bm{\widehat{u}}_{k})(1+\widehat{b_{k}^{{\sf
PCA}}})}{n(\widecheck{\lambda}_{j}-\widecheck{\lambda}_{k})^{2}}+\frac{\widehat{\sigma}^{2}\widecheck{\lambda}_{j}\|\bm{\widehat{U}}_{\perp}^{\top}\bm{a}\|^{2}}{(\widecheck{\lambda}_{j}-\widehat{\sigma}^{2})}.$

7: Compute ${\sf C.I.}^{(\alpha)}(\bm{\widehat{u}}_{j}^{\top}\bm{a})$ defined
via  $\displaystyle{\sf
C.I.}^{(\alpha)}(\bm{\widehat{u}}_{j}^{\top}\bm{a}):=\bigg{[}\bm{\widehat{u}}_{j}^{\top}\bm{a}\sqrt{1+\widehat{b_{j}^{{\sf
PCA}}}}\pm\Phi^{-1}\big{(}1-\alpha/2\big{)}\widehat{s^{{\sf
PCA}}_{\bm{a},j}}\bigg{]}$

8: return Estimated confidence interval ${\sf
C.I.}^{(\alpha)}(\bm{\widehat{u}}_{j}^{\top}\bm{a})$.

Algorithm 3 Confidence Interval for $\bm{a}^{\top}\bm{u}_{j}$ – PCA

######  Theorem 4.2 (Validity of Confidence Intervals – PCA).

Consider the model in (2) and suppose $\lambda_{1}$ through $\lambda_{r}$ are
unique. Suppose further that

$\displaystyle\lambda_{\min}$ $\displaystyle\geq
C_{1}\max\bigg{\\{}\sigma^{2}\kappa^{3}r^{2}\log^{7}(n\vee
p)\bigg{(}\frac{p}{n}+\sqrt{\frac{p}{n}}\bigg{)},\sigma\kappa\log^{4}(n\vee
p)\frac{p}{n}\bigg{\\}};$ $\displaystyle\Delta_{\min}$ $\displaystyle\geq
C_{1}\big{(}\lambda_{\max}+\sigma^{2})\frac{r^{5/2}}{\sqrt{n}}\log^{4}(n\vee
p),$

where $C_{1}$ is some sufficiently large constant. In addition, assume that
$r^{3}\leq c_{1}\frac{n}{\kappa^{6}\log^{12}(n\vee p)}$, where $c_{1}$ is a
sufficiently small constant, and that $\log(p)\lesssim n^{1/2}$. Suppose that

$\displaystyle\max_{k\neq j,k\leq r}\frac{s_{\bm{a},k}^{{\sf
PCA}}}{s_{\bm{a},j}^{{\sf
PCA}}}\frac{(\lambda_{\max}+\sigma^{2})\sqrt{r}\log^{2}(n\vee
p)}{\Delta_{j}\sqrt{n}}=o(1).$

Finally, assume that

$\displaystyle\frac{|\bm{a}^{\top}\bm{u}_{j}|}{s_{\bm{a},j}^{{\sf PCA}}}$
$\displaystyle\Bigg{[}\frac{(\lambda_{\max}+\sigma^{2})(\lambda_{j}+\sigma^{2})r\log(n\vee
p)}{\Delta_{j}^{2}n}$
$\displaystyle\quad+\frac{\sigma(\lambda_{j}+\sigma^{2})^{1/2}}{\lambda_{j}\sqrt{n}}\bigg{(}\frac{\sigma\kappa\sqrt{r}\log^{2}(n\vee
p)}{\lambda_{\min}^{1/2}}\sqrt{\frac{p}{n}}+\frac{\sigma^{2}\sqrt{r\kappa}\log^{2}(n\vee
p)}{\lambda_{\min}}\sqrt{\frac{p}{n}}\bigg{)}$
$\displaystyle\qquad+\frac{\sigma^{2}p}{\lambda_{\min}n}\frac{\sqrt{r}\kappa\sqrt{\log(n\vee
p)}}{n}\Bigg{]}=o(1).$

Then if either $\|\bm{U}^{\top}\bm{a}\|>\|\bm{U}_{\perp}^{\top}\bm{a}\|$ with
$|\bm{a}^{\top}\bm{u}_{j}|\lesssim|\bm{a}^{\top}\bm{u}_{k}|$ for $k\neq
j,k\leq r$ or $\|\bm{U}^{\top}\bm{a}\|\leq\|\bm{U}_{\perp}^{\top}\bm{a}\|$,
then it holds that

$\displaystyle\bigg{|}\mathbb{P}\bigg{\\{}\bm{a}^{\top}\bm{u}_{j}\in{\sf
C.I.}^{(\alpha)}(\bm{a}^{\top}\bm{\widehat{u}}_{j})\bigg{\\}}-(1-\alpha)\bigg{|}=o(1).$

The assumptions in Theorem 4.2 are similar to those of Theorem 4.1. As in the
setting therein, different alignment of the vector $\bm{a}$ results in
qualitatively different phenomena.

######  Remark 4.4 (Generality).

Similar to the matrix denoising setting, we first prove a more general result
in the form of a Berry-Esseen bound. See Lemma 8.10 for details.

######  Remark 4.5 (Optimality and Statistical Efficiency).

As proven in Koltchinskii et al. (2020), the limiting variance
$(s_{\bm{a},j}^{{\sf PCA}})^{2}$ is optimal for a broad class of covariance
matrices with “large” eigengaps. As a result, our confidence intervals are
effectively the shortest possible asymptotically valid confidence intervals
for a subclass of the class of covariance matrices studied therein. Moreover,
unlike Koltchinskii et al. (2020), our confidence intervals do not require
_any_ sample splitting, whereas Koltchinskii et al. (2020) require separating
some small (in $n$) sample to estimate the debiasing parameter. While
theoretically both results demonstrate their respective asymptotic optimality,
in practice any amount of sample splitting still results in a small loss of
finite-sample efficiency.

Finally, we return to the entrywise context discussed in Corollary 3.2. The
following result demonstrates the validity of entrywise confidence intervals
under incoherence conditions. The proof can be found in Section C.9.

######  Corollary 4.2 (Entrywise Confidence Intervals – PCA).

Consider the model (2) and suppose $\lambda_{1}$ through $\lambda_{r}$ are
unique. Suppose further that

$\displaystyle\lambda_{\min}$ $\displaystyle\geq
C_{1}\max\bigg{\\{}\sigma^{2}\kappa^{3}r^{2}\log^{7}(n\vee
p)\bigg{(}\frac{p}{n}+\sqrt{\frac{p}{n}}\bigg{)},\sigma\kappa\log^{4}(n\vee
p)\frac{p}{n}\bigg{\\}};$ $\displaystyle\Delta_{\min}$ $\displaystyle\geq
C_{1}\big{(}\lambda_{\max}+\sigma^{2})\kappa^{3/2}\frac{r^{5/2}}{\sqrt{n}}\log^{4}(n\vee
p),$

where $C_{1}$ is some sufficiently large constant. In addition, assume that
$r^{3}\leq c_{1}\frac{n}{\kappa^{6}\log^{12}(n\vee p)}$, that $r\leq
c_{1}p^{1/2}$, where $c_{1}$ is a sufficiently small constant, and that
$\log(p)\lesssim n^{1/2}$. Assume there exists a constant $\mu=O(1)$ such that

$\displaystyle\max_{k\leq r}|\bm{u}_{k}(i)|\leq\mu\sqrt{\frac{1}{p}}.$

Suppose further that

$\displaystyle\max_{k,l\leq r}\frac{|\bm{u}_{l}(i)|}{|\bm{u}_{k}(i)|}=O(1).$

Let ${\sf C.I.}^{(\alpha)}(\bm{\widehat{u}}_{j}(i))$ denote the output of
Algorithm 4. Then it holds that

$\displaystyle\bigg{|}\mathbb{P}\bigg{\\{}\bm{u}_{j}(i)\in{\sf
C.I.}^{(\alpha)}(\bm{\widehat{u}}_{j}(i))\bigg{\\}}-(1-\alpha)\bigg{|}=o(1).$

Concretely, Corollary 4.2 enables one to obtain asymptotically valid entrywise
confidence intervals when the corresponding entries of the leading
eigenvectors are of similar order to each other.

0: Sample data $\bm{X}\in\mathbb{R}^{p\times n}$, rank $r$, index $j\in[r]$,
unit vector $\bm{a}\in\mathbb{R}^{p}$, desired confidence level $\alpha$.

1: Compute the eigendecomposition of
$\bm{\widehat{\Sigma}}=\frac{\bm{XX}^{\top}}{n}=\big{[}\bm{\widehat{U}},\bm{\widehat{U}}_{\perp}\big{]}\bm{\widehat{\Lambda}}\big{[}\bm{\widehat{U}},\bm{\widehat{U}}_{\perp}\big{]}^{\top}$,
where $\bm{\widehat{U}}$ has columns consisting of the $r$ eigenvectors
corresponding to the largest nonzero eigenvalues. Let $\bm{\widehat{u}}_{k}$
denote the $k$’th column of $\bm{\widehat{U}}$ and $\widehat{\lambda}_{k}$
denote the $k$’th largest eigenvalue of $\bm{\widehat{\Sigma}}$.

2: Set $\widehat{\sigma}^{2}$ via Eq. 8.

3: Set $\widehat{b_{k}^{{\sf PCA}}}$ via Eq. 9 for $1\leq k\leq r$, $k\neq j$.

4: Compute, for $1\leq k\leq r$,  $\displaystyle\widehat{\gamma}^{{\sf
PCA}}(\widehat{\lambda}_{k}):=\frac{1}{n}\sum_{i=r+1}^{\min(p-r,n)}\frac{\widehat{\lambda}_{i}}{\widehat{\lambda}_{k}-\widehat{\lambda}_{i}}$

5: Compute, for $1\leq k\leq r$,
$\displaystyle\widecheck{\lambda}_{k}:=\frac{\widehat{\lambda}_{k}}{1+\widehat{\gamma}^{{\sf
PCA}}(\widehat{\lambda}_{k})}.$

6: Set  $\displaystyle\big{(}\widehat{s_{\bm{e}_{i},j}^{{\sf
PCA}}}\big{)}^{2}:=\sum_{k\leq r,k\neq
j}\frac{\widecheck{\lambda}_{k}\widecheck{\lambda}_{j}(\bm{\widehat{u}}_{k}(i))(1+\widehat{b_{k}^{{\sf
PCA}}})}{n(\widecheck{\lambda}_{j}-\widecheck{\lambda}_{k})^{2}}+\frac{\widehat{\sigma}^{2}\widecheck{\lambda}_{j}\|\bm{\widehat{U}}_{\perp}^{\top}\bm{e}_{i}\|^{2}}{(\widecheck{\lambda}_{j}-\widehat{\sigma}^{2})}.$

7: Compute ${\sf C.I.}^{(\alpha)}(\bm{\widehat{u}}_{j}(i))$ defined via
$\displaystyle{\sf
C.I.}^{(\alpha)}(\bm{\widehat{u}}_{j}(i)):=\bigg{[}\bm{\widehat{u}}_{j}(i)\pm\Phi^{-1}\big{(}1-\alpha/2\big{)}\widehat{s^{{\sf
PCA}}_{\bm{e}_{i},j}}\bigg{]}$

8: return Estimated confidence interval ${\sf
C.I.}^{(\alpha)}(\bm{\widehat{u}}_{j}(i))$.

Algorithm 4 Entrywise Confidence Intervals for $\bm{u}_{j}$ – PCA

##  5 Previous Work

The analysis of spectral methods has a rich history, going back to $\ell_{2}$
perturbation theory such as the Davis-Kahan Theorem (Davis and Kahan, 1970) or
Wedin’s Theorem (Wedin, 1972), and much of this theory now forms the basis of
book-level treatments (e.g. Stewart and Sun (1990); Kato (1995)). Classical
deterministic matrix perturbation theory is remarkably effective in a number
of statistical settings, and improvements on the deterministic theory have
been obtained (Cai and Zhang, 2018; Luo et al., 2021; Yu et al., 2015; Zhang,
2022), as well as improvements tailored explicitly for various types of random
noise (O’Rourke et al., 2018; Koltchinskii and Xia, 2016; Nadler, 2008).
Notably, a number of works have studied the $\ell_{2}$ perturbation in
principal component analysis (Cai et al., 2021b, 2013; Koltchinskii and
Lounici, 2017b). Many of these works focus on providing bounds for estimated
subspaces (or for eigenvectors with large eigengaps), though we single out
O’Rourke et al. (2023), which provides optimal $\ell_{2}$ perturbation bounds
for individual eigenvectors in a matrix denoising context under similar
eigengap conditions as in this work.

Beyond $\ell_{2}$ perturbation theory, a number of works have also considered
$\ell_{2,\infty}$ perturbation theory, which quantifies the entrywise (or row-
wise) fluctuations of eigenvectors (or subspaces). The works Cape et al.
(2019b); Fan et al. (2018); Eldridge et al. (2018); Damle and Sun (2020);
Bhardwaj and Vu (2023) provide deterministic $\ell_{2,\infty}$ bounds for the
leading subspaces of low-rank matrices and apply their results to problems
arising in high-dimensional statistics. Furthermore, by taking into account
both the structural assumptions and the probabilistic nature of the noise, a
number of works have also derived $\ell_{2,\infty}$ perturbation bounds
tailored to different probabilistic settings. For example, Lei (2019)
considers symmetric matrices with independent noise, Cai et al. (2021a) and
Zhou and Chen (2023) study highly unbalanced matrices, Agterberg and Sulam
(2022) focus on matrices with additional sparse structure, Agterberg and Zhang
(2022) examine problems arising from tensor data analysis, and Chen et al.
(2019, 2022) study ranking problems, though this list is still incomplete. A
general survey on applications of $\ell_{2,\infty}$ perturbation theory in
various statistical contexts can be found in Chen et al. (2021b).

The $\ell_{2,\infty}$ subspace perturbation bounds considered in previous
works are a special type of the more general linear functions studied in this
work. The works Koltchinskii and Xia (2016) and Xia and Zhou (2019) study
perturbation bounds for linear functions of eigenvectors, with the latter
focusing on the “unbalanced” setting, wherein the column dimension is
significantly larger than the row dimension. Similarly, Koltchinskii and
Lounici (2016) provides perturbation bounds for linear functions of principal
components under Gaussian noise in a general Hilbert space. However, these
previous results all operate in the “large eigengaps” regime, implicitly
requiring that the eigengap $\Delta_{j}$ is of similar order to the eigenvalue
spike $\lambda_{r}$. The works Cheng et al. (2021); Chen et al. (2021a); Li et
al. (2022) all provide perturbation bounds for linear functions of
eigenvectors in the “small eigengaps” regime, though both Cheng et al. (2021)
and Chen et al. (2021a) require _asymmetric_ noise, something that we do not
consider herein.

This work is also closely related to a number of works studying the asymptotic
theory for low-rank matrix models. For example, Xia (2021, 2019); Bao et al.
(2021) consider asymptotic theory for the $\sin\Theta$ distances in matrix
denoising, and Bao et al. (2022); Koltchinskii and Lounici (2017b) study the
asymptotic distribution of $\sin\Theta$ distances for the spiked principal
component analysis model. These works also pertain to the so-called
“supercritical” regime in random matrix theory, for which there are a number
of works studying the asymptotics for eigenvalues and eigenvectors in matrix
denoising (Ding, 2020; Capitaine, 2018; Benaych-Georges and Nadakuditi, 2011;
Benaych-Georges and Nadakuditi, 2012; Knowles and Yin, 2014; Capitaine et al.,
2012, 2009) and PCA (Ding and Yang, 2021; Paul, 2007; Bloemendal et al., 2016;
Baik et al., 2005; Ding, 2021; Onatski, 2012; Johnstone and Lu, 2009), though
these lists are incomplete.

Turning to asymptotic distributional theory, the works Koltchinskii and
Lounici (2017a); Ding and Yang (2021); Bao et al. (2022) study the asymptotics
for bilinear forms of spectral projectors in spiked covariance models, and
Koltchinskii et al. (2020) study efficient estimation in these settings. In a
matrix denoising setting, Fan et al. (2020) study linear functions of
eigenvectors of symmetric matrices when eigenvalues diverge. However, in the
aforementioned works, the emphasis has primarily been on the “large eigengaps”
regime, and the results typically do not hold for finite samples.

Finally, statistical inference in low-rank matrix models is a relatively
nascent field, with only a few results in specific contexts (Agterberg et al.,
2022; Yan et al., 2021; Carpentier et al., 2019; Silin and Fan, 2020;
Carpentier and Nickl, 2015; Carpentier et al., 2018; Chernozhukov et al.,
2023; Cai et al., 2023; Xia and Yuan, 2021). Considering estimating linear
functionals, Koltchinskii et al. (2020) show that it is possible to find
asymptotically efficient confidence intervals, though their estimators still
rely on sample splitting. Notably, the work Cheng et al. (2021) provides
inferential procedures for linear functions of eigenvectors under asymmetric,
heteroskedastic noise.

##  6 Numerical Simulations

In this section we demonstrate our theory through simulations. Each result is
generated through 200 Monte Carlo iterations.

###  6.1 Matrix Denoising

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 2: Empirical (dotted) and theoretical (solid) ellipses for the quantity
$\bm{a}^{\top}\bm{\widehat{u}}_{1}-\bm{a}^{\top}\bm{u}_{1}\bm{u}_{1}^{\top}\bm{\widehat{u}}_{1}$
with varying signal-strength and eigengaps under the matrix denoising model
(1). From top to bottom the eigengap increases, and from left to right the
signal strength increases.

In this section we demonstrate our theoretical results for matrix denoising.

Setup: In all simulations we set $n=200$ and $\sigma^{2}=1$, and $r=3$. To
generate the matrix $\bm{S}$, we first define a vector $\bm{v}$ with entries
consisting of $n$ independent $U(-1,1)$ random variables, and we define
$\bm{u}_{1}=\frac{\bm{v}}{\|\bm{v}\|}$. We then draw $\bm{u}_{2}$ by drawing a
new vector $\bm{v}^{\prime}$ with $U(-1,1)$ entries and setting
$\bm{u}_{2}=\big{(}\bm{I}-\bm{u}_{1}\bm{u}_{1}^{\top}\big{)}\frac{\bm{v}^{\prime}}{\|\bm{v}^{\prime}\|}$
to ensure orthogonality. Finally we set $\bm{u}_{3}$ by drawing another vector
$\bm{v}^{\prime\prime}$ with $U(-1,1)$ entries and setting
$\bm{u}_{3}=\big{(}\bm{I}-\bm{u}_{1}\bm{u}_{1}^{\top}-\bm{u}_{2}\bm{u}_{2}^{\top}\big{)}\frac{\bm{v}}{\|\bm{v}\|}$.
Finally, for the varying levels of signal-to-noise ratio
$\lambda_{\min}/\sigma\equiv\lambda_{\min}$ and eigengap $\Delta_{1}$, we
define $\lambda_{1}=5\lambda_{\min}$, $\lambda_{2}=\lambda_{1}-\Delta_{1}$,
$\lambda_{3}=\lambda_{\min}$. We first draw the matrix $\bm{S}$ before running
simulations, and we keep it fixed throughout all simulations. The matrix
$\bm{N}$ is re-drawn at each Monte Carlo iteration.

Distributional Theory: Plotted in Fig. 2 is the empirical histogram of the
random variable
$\bm{a}^{\top}\bm{\widehat{u}}_{1}-\bm{a}^{\top}\bm{u}_{1}\bm{u}_{1}^{\top}\bm{\widehat{u}}_{1}$,
with solid line corresponding to the theoretical distribution (mean zero and
variance $s^{{\sf MD}}_{\bm{a},1}$) and dotted line corresponding to the
empirical distribution, with $\bm{a}$ equal to the constant vector with
entries $\frac{1}{\sqrt{n}}$. By Theorem 3.1 the approximate Gaussianity is
governed by the eigengap $\Delta_{1}/\sigma$ and the signal strength
$\lambda_{\min}/\sigma$. From left to right we consider increasing signal
strength and from top to bottom we consider increasing eigengap, so that top
left has the smallest eigengap and signal strength, and bottom right has the
largest eigengap and signal strength. Furthermore, by Theorem 3.1, the
asymptotic variance decreases with larger eigengap and signal strength, and
observe that the theoretical and empirical histograms become narrower as these
quantities increase. We also ran simulations for the difference
$\bm{a}^{\top}\bm{\widehat{u}}_{1}\sqrt{1+b_{j}^{{\sf
MD}}}-\bm{a}^{\top}\bm{u}_{1}$, but the figures are not materially different.

Approximate Coverage Rates: We also consider approximate coverage rates for
$\bm{a}^{\top}\bm{u}_{1}$ using Algorithm 1 as well as the individual entry
$\bm{u}_{1}(1)$ using Algorithm 2. Similar to the previous setting we examine
increasing eigengaps and signal-strengths with $\alpha=.05$. The column
denoted “Mean” corresponds to the number of times the true value of either
$\bm{a}^{\top}\bm{u}_{1}$ or $\bm{u}_{1}(1)$ was in the outputted confidence
interval, and the column denoted “Std” denotes the standard deviation of this
value. Due to the sign ambiguity, we multiplied $\bm{a}^{\top}\bm{u}_{1}$ and
$\bm{u}_{1}(1)$ by ${\sf sign}(\bm{u}_{1}^{\top}\bm{\widehat{u}}_{1})$, as we
assume this quantity is positive throughout our theory. In both cases the
coverage remains close to $.95$.

Confidence Intervals for Matrix Denoising
---
Coverage Rates for $\bm{a}^{\top}\bm{u}_{1}$
$\Delta_{1}$ | $\lambda_{\min}$ | Mean | Std
$\log(n)$ | $2\sqrt{n}$ | 0.955 | 0.015
$\log(n)$ | $4\sqrt{n}$ | 0.980 | 0.010
$\log(n)$ | $6\sqrt{n}$ | 0.990 | 0.007
$2\log(n)$ | $2\sqrt{n}$ | 0.935 | 0.017
$2\log(n)$ | $4\sqrt{n}$ | 0.965 | 0.013
$2\log(n)$ | $6\sqrt{n}$ | 0.975 | 0.011
$3\log(n)$ | $2\sqrt{n}$ | 0.900 | 0.021
$3\log(n)$ | $4\sqrt{n}$ | 0.965 | 0.013
$3\log(n)$ | $6\sqrt{n}$ | 0.990 | 0.007

Confidence Intervals for Matrix Denoising
---
Coverage Rates for $\bm{u}_{1}(1)$
$\Delta_{1}$ | $\lambda_{\min}$ | Mean | Std
$\log(n)$ | $2\sqrt{n}$ | 0.960 | 0.014
$\log(n)$ | $4\sqrt{n}$ | 0.98 | 0.010
$\log(n)$ | $6\sqrt{n}$ | 0.995 | 0.005
$2\log(n)$ | $2\sqrt{n}$ | 0.925 | 0.019
$2\log(n)$ | $4\sqrt{n}$ | 0.98 | 0.010
$2\log(n)$ | $6\sqrt{n}$ | 0.970 | 0.012
$3\log(n)$ | $2\sqrt{n}$ | 0.885 | 0.023
$3\log(n)$ | $4\sqrt{n}$ | 0.94 | 0.017
$3\log(n)$ | $6\sqrt{n}$ | 0.970 | 0.012

Table 1: Matrix Denoising empirical coverage rates for both
$\bm{a}^{\top}\bm{u}_{1}$ (left) and $\bm{u}_{1}(1)$ (right) using Algorithms
1 and 2 respectively for varying signal $\lambda_{\min}$ and eigengap
$\Delta_{1}$. Here $\bm{a}\equiv\frac{1}{\sqrt{n}}$. The column “Mean”
represents the empirical probability of coverage averaged over $200$ Monte
Carlo iterations, and the column “Std” denotes the standard deviation of this
coverage rate.

###  6.2 PCA

We next consider the PCA model.

Setup:  In all simulations we set $n=300$, $p=200$ and $r=3$, and we draw
$\bm{u}_{1},\bm{u}_{2},\bm{u}_{3}$ the same way as in the matrix denoising
simulations. We keep $\sigma^{2}=1$, and for a given eigengap $\Delta_{1}$ and
signal strength $\lambda_{\min}/\sigma^{2}\equiv\lambda_{\min}$, we set
$\lambda_{1}=5\lambda_{\min}$, $\lambda_{2}=\lambda_{1}-\Delta_{1}$, and
$\lambda_{3}=\lambda_{\min}$.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 4: Empirical (dotted) and theoretical (solid) ellipses for the quantity
$\bm{a}^{\top}\bm{\widehat{u}}_{1}-\bm{a}^{\top}\bm{u}_{1}\bm{u}_{1}^{\top}\bm{\widehat{u}}_{1}$
with varying signal-strength and eigengaps under the matrix denoising model
(2) From top to bottom the eigengap increases, and from left to right the
signal strength increases.

Distributional Theory: Plotted in Fig. 4 is the empirical histogram of the
random variable
$\bm{a}^{\top}\bm{\widehat{u}}_{1}-\bm{a}^{\top}\bm{u}_{1}\bm{u}_{1}^{\top}\bm{\widehat{u}}_{1}$,
with solid line corresponding to the theoretical distribution (mean zero and
variance $s_{\bm{a},1}^{{\sf PCA}}$), and dotted line corresponding to the
empirical distribution with $\bm{a}$ equal to the constant vector with entries
$\frac{1}{\sqrt{p}}$. From left to right we consider increasing signal
strength and from top to bottom we consider increasing eigengap. Similar to
matrix denoising, the Gaussian approximation improves and the theoretical
variance decreases as the eigengap and signal strength increase, which we see
demonstrated in the figure.

Approximate Coverage Rates: We also consider approximate coverage rates for
$\bm{a}^{\top}\bm{u}_{1}$ using Algorithm 3 as well as the individual entry
$\bm{u}_{1}(1)$ using Algorithm 4. Similar to the matrix denoising we examine
increasing eigengaps and signal-strengths with $\alpha=.05$. The column
denoted “Mean” corresponds to the number of times the true value of either
$\bm{a}^{\top}\bm{u}_{1}$ or $\bm{u}_{1}(1)$ was in the outputted confidence
interval, and the column denoted “Std” denotes the standard deviation of this
value. Due to the sign ambiguity, we multiplied $\bm{a}^{\top}\bm{u}_{1}$ and
$\bm{u}_{1}(1)$ by ${\sf sign}(\bm{u}_{1}^{\top}\bm{\widehat{u}}_{1})$, as we
assume this quantity is positive throughout our theory. Observe that
empirically the coverage improves as the signal strength and eigengap
increases.

PCA
---
Coverage Rates for $\bm{a}^{\top}\bm{u}_{1}$
$\frac{\Delta_{1}}{\lambda_{\min}+\sigma^{2}}$ | $\lambda_{\min}/\sigma^{2}$ | Mean | Std
$\frac{\log(n)}{\sqrt{n}}$ | $2\log(n)$ | 0.935 | 0.017
$\frac{\log(n)}{\sqrt{n}}$ | $4\log(n)$ | 0.92 | 0.019
$\frac{\log(n)}{\sqrt{n}}$ | $6\log(n)$ | 0.93 | 0.018
$2\frac{\log(n)}{\sqrt{n}}$ | $2\log(n)$ | 0.860 | 0.025
$2\frac{\log(n)}{\sqrt{n}}$ | $4\log(n)$ | 0.90 | 0.021
$2\frac{\log(n)}{\sqrt{n}}$ | $6\log(n)$ | 0.92 | 0.019
$3\frac{\log(n)}{\sqrt{n}}$ | $2\log(n)$ | 0.885 | 0.023
$3\frac{\log(n)}{\sqrt{n}}$ | $4\log(n)$ | 0.92 | 0.019
$3\frac{\log(n)}{\sqrt{n}}$ | $6\log(n)$ | 0.92 | 0.019

PCA
---
Coverage Rates for $\bm{u}_{1}(1)$
$\frac{\Delta_{1}}{\lambda_{\min}+\sigma^{2}}$ | $\lambda_{\min}/\sigma^{2}$ | Mean | Std
$\frac{\log(n)}{\sqrt{n}}$ | $2\log(n)$ | 0.895 | 0.022
$\frac{\log(n)}{\sqrt{n}}$ | $4\log(n)$ | 0.810 | 0.028
$\frac{\log(n)}{\sqrt{n}}$ | $6\log(n)$ | 0.855 | 0.025
$2\frac{\log(n)}{\sqrt{n}}$ | $2\log(n)$ | 0.935 | 0.017
$2\frac{\log(n)}{\sqrt{n}}$ | $4\log(n)$ | 0.945 | 0.016
$2\frac{\log(n)}{\sqrt{n}}$ | $6\log(n)$ | 0.925 | 0.019
$3\frac{\log(n)}{\sqrt{n}}$ | $2\log(n)$ | 0.930 | 0.018
$3\frac{\log(n)}{\sqrt{n}}$ | $4\log(n)$ | 0.905 | 0.021
$3\frac{\log(n)}{\sqrt{n}}$ | $6\log(n)$ | 0.900 | 0.021

Table 2: PCA empirical coverage rates for both $\bm{a}^{\top}\bm{u}_{1}$
(left) and $\bm{u}_{1}(1)$ (right) using Algorithm 1 and Algorithm 2
respectively for varying signal $\lambda_{\min}/\sigma^{2}$ and eigengap
$\Delta_{1}$. Here $\bm{a}\equiv\frac{1}{\sqrt{p}}$. The column “Mean”
represents the empirical probability of coverage averaged over $200$ Monte
Carlo iterations, and the column “Std” denotes the standard deviation of this
coverage rate.

##  7 Discussion

In this work we have considered estimating linear forms of the form
$\bm{a}^{\top}\bm{u}_{j}$ for a pre-specified unit vector $\bm{a}$. We have
seen how small eigengaps and signal strength affect the approximate
Gaussianity of both the plugin and debiased estimators, and we have provided
asymptotically valid and statistically efficient confidence intervals for
$\bm{a}^{\top}\bm{u}_{j}$. All of our results hold under nearly minimal
signal-strength conditions, and our proposed procedures are fully data-driven
and do not require any sample splitting.

There are a number of potential future works. First, our analysis relies
heavily on the assumption of Gaussian noise, which may not hold in general. It
is of interest to develop similar statistical theory under general noise
mechanisms such as subgaussian or heteroskedastic noise. Furthermore, our
results are likely suboptimal with respect to logarithmic terms, factors of
$r$, and factors of the condition number $\kappa$, and it would be interesting
to study the optimal dependence on these parameters. For example, in matrix
denoising, our analysis requires that
$\Delta_{j}/\sigma\gg\mathrm{polylog}(n)$ assuming $r,\kappa=O(1)$. Does
asymptotic normality hold even when $\Delta_{j}\asymp\sigma$ provided certain
necessary signal-strength conditions are met?

Beyond these immediate future directions, it would also be of interest to
study similar problems in other settings. For example, the work Xia et al.
(2022) studies statistical inference for linear functions of tensor singular
vectors; it would be interesting to see if similar guarantees can be developed
when the tensor singular value gaps are of similarly small order as considered
herein. Finally, the work Xia (2019) considers confidence regions for singular
subspaces in low-rank trace regression; it would be of interest to develop
similar guarantees for individual singular subspaces in the presence of small
singular value gaps.

##  8 Technical Tools and Proofs of Main Results

In this section we prove our key main results; namely, Theorems 3.1, 3.2, 3.3,
3.4, 4.1 and 4.2. Our main results rely on two entirely deterministic
expansions, which we discuss in the next section. We then prove the results
for matrix denoising and for PCA in Sections 8.2 and 8.3 respectively. Proofs
of additional lemmas and corollaries can be found in the appendix.

###  8.1 Matrix Analytic Tools

The bulk of our analysis is built off the following deterministic error
decomposition that yields the precise form of the leading-order term.

######  Theorem 8.1 (First-Order Expansion).

Let $\\{\bm{\widehat{u}}_{j},\widehat{\lambda}_{j}\\}$ denote the eigenvectors
and eigenvalues of the symmetric matrix $\bm{\widehat{M}}:=\bm{M}+\bm{E}$,
where the symmetric matrix $\bm{M}$ has eigenvectors and eigenvalues
$\bm{u}_{j}$ and $\lambda_{j}$ respectively. Suppose $\lambda_{j}$ is unique.
Then the following expansion is valid always:

$\displaystyle\bm{a}^{\top}\bm{\widehat{u}}_{j}-\bm{a}^{\top}\bm{u}_{j}\bm{u}_{j}^{\top}\bm{\widehat{u}}_{j}-\sum_{k\neq
j}\frac{\bm{u}_{j}^{\top}\bm{E}\bm{u}_{k}}{\lambda_{j}-\lambda_{k}}\bm{a}^{\top}\bm{u}_{k}$
$\displaystyle=\sum_{k\neq
j}\frac{(\bm{\widehat{u}}_{j}-\bm{u}_{j})^{\top}\bm{E}\bm{u}_{k}}{\lambda_{j}-\lambda_{k}}\bm{a}^{\top}\bm{u}_{k}$
$\displaystyle\quad+\sum_{k\neq
j}\frac{\lambda_{j}-\widehat{\lambda}_{j}}{\lambda_{j}-\lambda_{k}}\bm{\widehat{u}}_{j}^{\top}\bm{u}_{k}\bm{u}_{k}^{\top}\bm{a}.$

###### Proof.

See Section A.1. ∎

The particular form of the first-order error decomposition was motivated by
the proofs of Koltchinskii and Lounici (2016) who derive their results using
holomorphic functional calculus. In the setting considered herein the complex-
analytic techniques of Koltchinskii and Lounici (2016) do not apply as they
implicitly require that $\Delta_{j}\gtrsim\|\bm{E}\|$. In addition, a similar
first-order decomposition was used as an intermediate step in Jirak and Wahl
(2023), though their focus is on providing perturbation bounds and limit
theorems in non-Gaussian covariance models.

We will also use extensively the following matrix analysis identity that
allows us to quantify the influence of the null space of a rank $r$ matrix
$\bm{M}$ on the $j$’th estimated eigenvector $\bm{\widehat{u}}_{j}$.

######  Theorem 8.2.

Let $\bm{\widehat{M}}=\bm{M}+\bm{E}$, and let
$\\{\bm{\widehat{u}}_{j},\widehat{\lambda}_{j}\\}$ denote the eigenvectors and
eigenvalues of $\bm{\widehat{M}}$ and $\bm{u}_{j},\lambda_{j}$ denote the
eigenvectors and eigenvalues of $\bm{M}$. Suppose $\bm{M}$ is rank $r$ with
orthogonal complement $\bm{U}_{\perp}$, such that
$\bm{M}=\bm{U}\bm{\Lambda}\bm{U}^{\top}$. Then it holds that

$\displaystyle\bm{\widehat{u}}_{j}^{\top}\bm{U}_{\perp}$
$\displaystyle=\bm{\widehat{u}}_{j}^{\top}\bm{U}\bm{U}^{\top}\bm{E}\bm{U}_{\perp}^{\top}\big{(}\widehat{\lambda}_{j}\bm{I}_{n-r}-\bm{U}_{\perp}^{\top}\bm{E}\bm{U}_{\perp}\big{)}^{-1}$
$\displaystyle\equiv\bm{\widehat{u}}_{j}^{\top}\bm{u}_{k}\bm{u}_{k}^{\top}\bm{E}\bm{U}_{\perp}\big{(}\widehat{\lambda}_{j}\bm{I}_{n-r}-\bm{U}_{\perp}^{\top}\bm{E}\bm{U}_{\perp}\big{)}^{-1}$
$\displaystyle\quad+\bm{\widehat{u}}_{j}\bm{U}^{-k}(\bm{U}^{-k})^{\top}\bm{E}\bm{U}_{\perp}\big{(}\widehat{\lambda}_{j}\bm{I}_{n-r}-\bm{U}_{\perp}^{\top}\bm{E}\bm{U}_{\perp}\big{)}^{-1},$

provided the inverse is defined. Here $\bm{U}^{-k}$ corresponds to the matrix
$\bm{U}$ with its $k$’th column removed.

###### Proof.

See Section A.2. ∎

Theorem 8.2 allows for a decoupling of randomness as well as aids in
identifying the appropriate eigenvalue bias term. As an immediate corollary to
this result, we see that whenever $|\lambda_{j}|>3\|\bm{E}\|$ it holds that
$\widehat{\lambda}_{j}\bm{I}_{n-r}-\bm{U}_{\perp}^{\top}\bm{E}\bm{U}_{\perp}$
is invertible because

$\displaystyle|\widehat{\lambda}_{j}|\geq|\lambda_{j}|-\|\bm{E}\|>\frac{2|\lambda_{j}|}{3}>2\|\bm{E}\|\geq
2\|\bm{U}_{\perp}^{\top}\bm{E}\bm{U}_{\perp}\|,$

which also implies that
$|\lambda_{\min}\big{(}\widehat{\lambda}_{j}\bm{I}_{n-r}-\bm{U}_{\perp}^{\top}\bm{E}\bm{U}_{\perp}\big{)}|\geq\frac{|\lambda_{j}|}{3}$.
As a consequence,

$\displaystyle\|\bm{\widehat{u}}_{j}^{\top}\bm{U}_{\perp}\|$
$\displaystyle\leq\|\bm{U}^{\top}\bm{E}\bm{U}_{\perp}\|\|\big{(}\widehat{\lambda}_{j}-\bm{U}_{\perp}^{\top}\bm{E}\bm{U}_{\perp}\big{)}^{-1}\|$
$\displaystyle\leq\frac{3\|\bm{E}\|}{|\lambda_{j}|}.$

This result can be interpreted as a variant of the Davis-Kahan Theorem
demonstrating the error along the $j$’th subspace scales with
$|\lambda_{j}|^{-1}$.

###  8.2 Analysis for Matrix Denoising

Before embarking on the proof we first state several preliminary facts that
are standard from nonasymptotic random matrix theory. We will use these facts
repeatedly in our proofs. Without loss of generality we assume that
$\lambda_{j}$ is positive; otherwise repeat the analysis with
$-\bm{\widehat{S}}$ and $-\bm{S}$.

######  Fact 1.

It holds that $\|\bm{N}\|\lesssim\sigma\sqrt{n}$ with probability at least
$1-\exp(-cn)$. See, for example Vershynin (2018).

######  Fact 2.

It holds that
$\|\bm{U}^{\top}\bm{N}\bm{U}\|\lesssim\sigma\big{(}\sqrt{r}+\sqrt{\log(n)}\big{)}$
with probability at least $1-O(n^{-10})$. This can be proven via standard
$\varepsilon$-net arguments, or one can appeal directly to the fact that
$\bm{U}^{\top}\bm{N}\bm{U}$ is a symmetric $r\times r$ Gaussian Wigner matrix.

######  Fact 3.

By virtue of Fact 1 and Weyl’s inequality, it holds that
$4|\lambda_{j}|/3+C\sigma\sqrt{n}\geq|\widehat{\lambda}_{j}|\geq|\lambda_{j}|-C\sigma\sqrt{n}\geq
2|\lambda_{j}|/3$ for all $j\leq r$ and
$|\widehat{\lambda}_{j}|\lesssim\sigma\sqrt{n}$ for all $j\geq r+1$.
Therefore, on the event $\|\bm{N}\|\lesssim\sigma\sqrt{n}$, it holds that
$|\widehat{\lambda}_{j}|\in\big{[}2|\lambda_{j}|/3,4|\lambda_{j}|/3\big{]}$.

######  Fact 4.

By virtue of Fact 3, it holds that $|\widehat{\lambda}_{j}|\geq
2|\lambda_{\min}|/3\geq\|\bm{N}\|\geq\|\bm{U}_{\perp}^{\top}\bm{N}\bm{U}_{\perp}\|$,
and hence $\lambda\bm{I}_{n-r}-\bm{U}_{\perp}^{\top}\bm{N}\bm{U}_{\perp}$ is
invertible for all $|\lambda|\geq 2\lambda_{\min}/3$. Specifically,
$\widehat{\lambda}_{j}\bm{I}_{n-r}-\bm{U}_{\perp}^{\top}\bm{N}\bm{U}_{\perp}$
is invertible on the event $\|\bm{N}\|\lesssim\sigma\sqrt{n}$. In addition,
$\|\big{(}\widehat{\lambda}_{j}\bm{I}_{n-r}-\bm{U}_{\perp}^{\top}\bm{N}\bm{U}_{\perp}\big{)}^{-1}\|\lesssim\lambda_{j}^{-1}$.

With these preliminary results out of the way we are now prepared to start the
proof of Theorems 3.1 and 3.2 and Theorem 4.1. The first part of our proof
will be to characterize both the bias
$\bm{\widehat{u}}_{j}^{\top}\bm{u}_{j}-1$ as well as the leading eigenvalues,
both of which will be used in the subsequent analysis. These steps rely
heavily on the preliminary analysis in Li et al. (2022).

Throughout our proofs of our distributional theory results we assume that
$\lambda_{1}$ through $\lambda_{r}$ are unique. If not, the extensions are
only more notationally cumbersome and not significantly different.

####  8.2.1 Step 1: Characterizing the Bias

Theorem 8.1 yields an expansion for $\bm{a}^{\top}\bm{\widehat{u}}_{j}$ about
$\bm{a}^{\top}\bm{u}_{j}\bm{u}_{j}^{\top}\bm{\widehat{u}}_{j}$, which
contributes to the bias $\bm{u}_{j}^{\top}\bm{\widehat{u}}_{j}-1$ (recall that
we assume that $\bm{u}_{j}^{\top}\bm{\widehat{u}}_{j}\geq 0$). Therefore, the
following lemma describes the bias.

######  Lemma 8.1 (Bias Characterization - Matrix Denoising).

Instate the conditions in Theorem 3.1. Define

$\displaystyle b_{j}^{{\sf MD}}:$
$\displaystyle=\sum_{k>r}\frac{\sigma^{2}}{(\widehat{\lambda}_{j}-\widehat{\lambda}_{k})^{2}}.$

Then with probability at least $1-O(n^{-10})$ it holds that

$\displaystyle|1-(\bm{\widehat{u}}_{j}^{\top}\bm{u}_{j})^{2}|$
$\displaystyle\lesssim\frac{\sigma^{2}n}{\lambda_{j}^{2}}+\frac{\sigma^{2}r\log(n)}{\Delta_{j}^{2}};$
$\displaystyle\big{|}1-\sqrt{1+b_{j}^{{\sf
MD}}}\bm{u}_{j}^{\top}\bm{\widehat{u}}_{j}\big{|}$
$\displaystyle\lesssim\frac{\sigma^{2}r\log(n)}{\Delta_{j}^{2}}+\frac{\sigma^{2}\sqrt{n\log(n)}}{\lambda_{j}^{2}}.$

###### Proof.

See Section B.1. ∎

Lemma 8.1 also implies that:

  * •

The term
$\bm{u}_{j}^{\top}\bm{\widehat{u}}_{j}\approx\frac{1}{\sqrt{1+b_{j}^{{\sf
MD}}}}$ up to some small error term.

  * •

The bias satisfies $b_{j}^{{\sf
MD}}\lesssim\frac{\sigma^{2}n}{\lambda_{j}^{2}}$ (by Facts 1 and 4).

####  8.2.2 Step 2: Eigenvalue Concentration

While Theorem 8.1 is useful for determining the first-order error, it is often
not quite enough to prove that the remaining two terms are sufficiently small
order. In order to further analyze these residuals, we will introduce an
additional “centering” term associated to the eigenvalue bias. Therefore, the
following lemma characterizes this eigenvalue bias as well as its empirical
approximation.

######  Lemma 8.2 (Eigenvalue Concentration - Matrix Denoising).

Instate the conditions of Theorem 3.1. Define

$\displaystyle\gamma^{{\sf
MD}}(\widehat{\lambda}_{j}):=\sigma^{2}{\sf{Tr}}\bigg{(}\big{(}\widehat{\lambda}_{j}\bm{I}_{n-r}-\bm{U}_{\perp}^{\top}\bm{N}\bm{U}_{\perp}\big{)}^{-1}\bigg{)};\qquad\widehat{\gamma}^{{\sf
MD}}(\widehat{\lambda}_{j}):=\sum_{k>r}\frac{\sigma^{2}}{\widehat{\lambda}_{j}-\widehat{\lambda}_{k}}.$

Then with probability at least $1-O(n^{-10})$ it holds that

$\displaystyle|\widehat{\lambda}_{j}-\lambda_{j}-\gamma^{{\sf
MD}}(\widehat{\lambda}_{j})|$
$\displaystyle\lesssim\sigma\big{(}\sqrt{r}+\sqrt{\log(n)}\big{)}:=\delta^{{\sf
MD}}.$

Furthermore, with this same probability it holds that

$\displaystyle|\gamma^{{\sf MD}}(\widehat{\lambda}_{j})-\widehat{\gamma}^{{\sf
MD}}(\widehat{\lambda}_{j})|$
$\displaystyle\lesssim\frac{\sigma^{2}r}{\lambda_{j}}.$

###### Proof.

See Section B.2. ∎

In words, Lemma 8.2 demonstrates that:

  * •

The empirical eigenvalues $\widehat{\lambda}_{j}$ concentrate around
$\lambda_{j}+\gamma^{{\sf MD}}(\widehat{\lambda}_{j})$ with error of order
$\sigma(\sqrt{r}+\sqrt{\log(n)})$.

  * •

Moreover, when $\sigma^{2}$ is known, $\gamma^{{\sf
MD}}(\widehat{\lambda}_{j})$ can be estimated directly from data via
$\widehat{\gamma}^{{\sf MD}}(\widehat{\lambda}_{j})$ with error at most
$O\bigg{(}\frac{\sigma^{2}r}{\lambda_{j}}\bigg{)}$.

####  8.2.3 Step 3: Finding the Leading-Order Expansion

Now that we have characterized the bias
$1-\bm{u}_{j}^{\top}\bm{\widehat{u}}_{j}$ and the empirical eigenvalues, we
are prepared to obtain the leading-order expansion for the plugin estimator
$\bm{a}^{\top}\bm{\widehat{u}}_{j}$ about
$\bm{a}^{\top}\bm{u}_{j}\bm{u}_{j}^{\top}\bm{\widehat{u}}_{j}$. Invoke Theorem
8.1 (with $\bm{M}=\bm{S}$ and $\bm{E}=\bm{N}$) to obtain that

$\displaystyle\bm{a}^{\top}\bm{\widehat{u}}_{j}-$
$\displaystyle\bm{a}^{\top}\bm{u}_{j}\bm{u}_{j}^{\top}\bm{\widehat{u}}_{j}-\sum_{k\neq
j}\frac{\bm{u}_{j}^{\top}\bm{N}\bm{u}_{k}}{\lambda_{j}-\lambda_{k}}\bm{u}_{k}^{\top}\bm{a}$
$\displaystyle=\sum_{k\neq
j}\frac{(\bm{\widehat{u}}_{j}-\bm{u}_{j})^{\top}\bm{N}\bm{u}_{k}}{\lambda_{j}-\lambda_{k}}\bm{a}^{\top}\bm{u}_{k}+\sum_{k\neq
j}\frac{\lambda_{j}-\widehat{\lambda}_{j}}{\lambda_{j}-\lambda_{k}}\bm{\widehat{u}}_{j}^{\top}\bm{u}_{k}\bm{u}_{k}^{\top}\bm{a}$
$\displaystyle=\sum_{k\neq j,k\leq
r}\frac{(\bm{\widehat{u}}_{j}-\bm{u}_{j})^{\top}\bm{N}\bm{u}_{k}-\gamma^{{\sf
MD}}(\widehat{\lambda}_{j})\bm{\widehat{u}}_{j}^{\top}\bm{u}_{k}}{\lambda_{j}-\lambda_{k}}\bm{a}^{\top}\bm{u}_{k}+\frac{(\bm{\widehat{u}}_{j}-\bm{u}_{j})^{\top}\bm{N}\bm{U}_{\perp}\bm{U}_{\perp}^{\top}\bm{a}}{\lambda_{j}}$
$\displaystyle\quad+\sum_{k\neq j,k\leq
r}\frac{\lambda_{j}-\widehat{\lambda}_{j}+\gamma^{{\sf
MD}}(\widehat{\lambda}_{j})}{\lambda_{j}-\lambda_{k}}\bm{\widehat{u}}_{j}^{\top}\bm{u}_{k}\bm{u}_{k}^{\top}\bm{a}+\frac{\lambda_{j}-\widehat{\lambda}_{j}}{\lambda_{j}}\bm{\widehat{u}}_{j}^{\top}\bm{U}_{\perp}\bm{U}_{\perp}^{\top}\bm{a}$
$\displaystyle:=\mathcal{R}^{{\sf MD}}_{1}+\mathcal{R}^{{\sf
MD}}_{2}+\mathcal{R}^{{\sf MD}}_{3}+\mathcal{R}^{{\sf MD}}_{4},$

where we have subtracted $\gamma^{{\sf
MD}}(\widehat{\lambda}_{j})\bm{\widehat{u}}_{j}^{\top}\bm{u}_{k}$ from the
term containing $(\bm{\widehat{u}}_{j}-\bm{u}_{j})^{\top}\bm{N}\bm{u}_{k}$ and
added it to the term containing $\lambda_{j}-\widehat{\lambda}_{j}$ for $k\leq
r$. To understand the intuition behind adding and subtracting this quantity,
observe that by Theorem 8.2,

$\displaystyle(\bm{\widehat{u}}_{j}-$
$\displaystyle\bm{u}_{j})^{\top}\bm{N}\bm{u}_{k}$
$\displaystyle=(\bm{\widehat{u}}_{j}-\bm{u}_{j})\bm{U}\bm{U}^{\top}\bm{N}\bm{u}_{k}+\bm{\widehat{u}}_{j}^{\top}\bm{U}_{\perp}\bm{U}_{\perp}^{\top}\bm{N}\bm{u}_{k}$
$\displaystyle=(\bm{\widehat{u}}_{j}-\bm{u}_{j})\bm{U}\bm{U}^{\top}\bm{N}\bm{u}_{k}+\bm{\widehat{u}}_{j}^{\top}\bm{U}\bm{U}^{\top}\bm{N}\bm{U}_{\perp}^{\top}\big{(}\widehat{\lambda}_{j}\bm{I}_{n-r}-\bm{U}_{\perp}^{\top}\bm{N}\bm{U}_{\perp}\big{)}^{-1}\bm{U}_{\perp}^{\top}\bm{N}\bm{u}_{k}$
$\displaystyle=\underbrace{(\bm{\widehat{u}}_{j}-\bm{u}_{j})\bm{U}\bm{U}^{\top}\bm{N}\bm{u}_{k}}_{\text{approximately
mean-
zero}}+\underbrace{\bm{\widehat{u}}_{j}^{\top}\bm{U}^{-k}(\bm{U}^{-k})^{\top}\bm{N}\bm{U}_{\perp}^{\top}\big{(}\widehat{\lambda}_{j}\bm{I}_{n-r}-\bm{U}_{\perp}^{\top}\bm{N}\bm{U}_{\perp}\big{)}^{-1}\bm{U}_{\perp}^{\top}\bm{N}\bm{u}_{k}}_{\text{approximately
mean-zero}}$
$\displaystyle\quad+\underbrace{\bm{\widehat{u}}_{j}^{\top}\bm{u}_{k}\bm{u}_{k}^{\top}\bm{N}\bm{U}_{\perp}^{\top}\big{(}\widehat{\lambda}_{j}\bm{I}_{n-r}-\bm{U}_{\perp}^{\top}\bm{N}\bm{U}_{\perp}\big{)}^{-1}\bm{U}_{\perp}^{\top}\bm{N}\bm{u}_{k}}_{\text{approximately
mean $\gamma^{{\sf MD}}(\widehat{\lambda}_{j})$}}.$

The first quantity can be shown to be sufficiently small using the fact that
$\bm{U}^{\top}\bm{N}\bm{u}_{k}$ is an $r$-dimensional random vector.
Similarly, the second term consists of products containing
$(\bm{U}^{-k})^{\top}\bm{N}\bm{U}_{\perp}$,
$\bm{U}_{\perp}^{\top}\bm{N}\bm{U}_{\perp}$, and
$\bm{U}_{\perp}^{\top}\bm{N}\bm{u}_{k}$, all of which are mutually independent
by virtue of the rotational invariance of the Gaussian distribution. However,
the remaining term
$\bm{u}_{k}^{\top}\bm{N}\bm{U}_{\perp}^{\top}\big{(}\widehat{\lambda}_{j}-\bm{U}_{\perp}^{\top}\bm{N}\bm{U}_{\perp}\big{)}^{-1}\bm{U}_{\perp}^{\top}\bm{N}\bm{u}_{k}$
is _not_ approximately mean-zero, and in fact can be shown to concentrate
about $\gamma^{{\sf MD}}(\widehat{\lambda}_{j})$, which is what motivates
adding and subtracting this quantity.

The remainder of the argument requires demonstrating that $\mathcal{R}^{{\sf
MD}}_{1}+\mathcal{R}^{{\sf MD}}_{2}+\mathcal{R}^{{\sf
MD}}_{3}+\mathcal{R}^{{\sf MD}}_{4}=o(s_{\bm{a},j}^{{\sf MD}})$. The following
lemma accomplishes this. The proof of this lemma uses Theorem 8.2 and Lemmas
8.2 and 8.1 extensively.

######  Lemma 8.3.

Under the conditions of Theorem 3.1, with probability at least $1-O(n^{-9})$
it holds that

$\displaystyle|\mathcal{R}^{{\sf MD}}_{1}+\mathcal{R}^{{\sf
MD}}_{2}+\mathcal{R}^{{\sf MD}}_{3}+\mathcal{R}^{{\sf MD}}_{4}|$
$\displaystyle\lesssim\bigg{(}\frac{\sigma
r\sqrt{n\log(n)}}{\lambda_{\min}}+\frac{\sigma
r^{3/2}\log(n)}{\Delta_{j}}\bigg{)}s_{\bm{a},j}^{{\sf MD}}$
$\displaystyle={\sf ErrMD}\times s_{\bm{a},j}^{{\sf MD}}.$

###### Proof.

See Section B.3. ∎

Here we slightly abuse notation by letting ${\sf ErrMD}$ absorb the implicit
constants.

####  8.2.4 Step 4: Completing the Proofs of Theorems 3.1 and 3.2

We have shown thus far that with probability at least $1-O(n^{-9}),$

$\displaystyle\frac{1}{s_{\bm{a},j}^{{\sf
MD}}}\bigg{(}\bm{a}^{\top}\bm{\widehat{u}}_{j}-\bm{a}^{\top}\bm{u}_{j}\bm{u}_{j}^{\top}\bm{\widehat{u}}_{j}\bigg{)}=\frac{1}{s_{\bm{a},j}^{{\sf
MD}}}\sum_{k\neq
j}\frac{\bm{u}_{j}^{\top}\bm{N}\bm{u}_{k}}{\lambda_{j}-\lambda_{k}}\bm{u}_{k}^{\top}\bm{a}+{\sf
ErrMD}$

In addition, it is straightforward to observe that the leading-order term
$\sum_{k\neq
j}\frac{\bm{u}_{j}^{\top}\bm{N}\bm{u}_{k}}{\lambda_{j}-\lambda_{k}}\bm{u}_{k}^{\top}\bm{a}$
is Gaussian with variance $(s_{\bm{a},j}^{{\sf MD}})^{2}$, as
$\bm{u}_{j}^{\top}\bm{N}\bm{u}_{k}$ is independent from
$\bm{u}_{j}^{\top}\bm{N}\bm{u}_{l}$ for $l\neq k$ by rotational invariance. As
a consequence, letting $\Phi(\cdot)$ denote the CDF of the Gaussian
distribution, it holds that

$\displaystyle\bigg{|}\mathbb{P}\bigg{\\{}\frac{1}{s_{\bm{a},j}^{{\sf MD}}}$
$\displaystyle\big{(}\bm{a}^{\top}\bm{\widehat{u}}_{j}-\bm{a}^{\top}\bm{u}_{j}\bm{u}_{j}^{\top}\bm{\widehat{u}}_{j}\big{)}\leq
z\bigg{\\}}-\Phi(z)\bigg{|}$
$\displaystyle=\bigg{|}\mathbb{P}\bigg{\\{}\frac{1}{s_{\bm{a},j}^{{\sf
MD}}}\bigg{(}\sum_{k\neq
j}\frac{\bm{u}_{j}^{\top}\bm{N}\bm{u}_{k}}{\lambda_{j}-\lambda_{k}}\bm{u}_{k}^{\top}\bm{a}+\mathcal{R}^{{\sf
MD}}_{1}+\mathcal{R}^{{\sf MD}}_{2}+\mathcal{R}^{{\sf
MD}}_{3}+\mathcal{R}^{{\sf MD}}_{4}\bigg{)}\leq z\bigg{\\}}-\Phi(z)\bigg{|}$
$\displaystyle\leq\bigg{|}\mathbb{P}\bigg{\\{}\frac{1}{s_{\bm{a},j}^{{\sf
MD}}}\bigg{(}\sum_{k\neq
j}\frac{\bm{u}_{j}^{\top}\bm{N}\bm{u}_{k}}{\lambda_{j}-\lambda_{k}}\bm{u}_{k}^{\top}\bm{a}\bigg{)}\leq
z\pm{\sf ErrMD}\bigg{\\}}-\Phi(z)\bigg{|}+O(n^{-9})$
$\displaystyle\leq\bigg{|}\Phi(z\pm{\sf ErrMD})-\Phi(z)\bigg{|}+O(n^{-9})$
$\displaystyle\lesssim{\sf ErrMD}+n^{-9},$

where the final line follows from the Lipschitz property of $\Phi$, where
$\Phi(z\pm x)$ is interpreted as applying the result to both $\Phi(z+x)$ and
$\Phi(z-x)$ separately. This completes the proof of Theorem 3.1.

We now prove Theorem 3.2. By Lemma 8.1 with probability at least
$1-O(n^{-10})$ it holds that

$\displaystyle|\sqrt{1+b_{j}^{{\sf
MD}}}\bm{\widehat{u}}_{j}^{\top}\bm{u}_{j}-1|$
$\displaystyle\lesssim\frac{\sigma^{2}r\log(n)}{\Delta_{j}^{2}}+\frac{\sigma^{2}\sqrt{n\log(n)}}{\lambda_{j}^{2}}={\sf
ErrBiasMD},$

where again we have absorbed implicit constants. As a result,

$\displaystyle\bigg{|}\mathbb{P}\bigg{\\{}$
$\displaystyle\frac{1}{s_{\bm{a},j}^{{\sf
MD}}}\bigg{(}\bm{a}^{\top}\bm{\widehat{u}}_{j}\sqrt{1+b_{j}^{{\sf
MD}}}-\bm{a}^{\top}\bm{u}_{j}\bigg{)}\leq z\bigg{\\}}-\Phi(z)\bigg{|}$
$\displaystyle\leq\bigg{|}\mathbb{P}\bigg{\\{}\frac{1}{s_{\bm{a},j}^{{\sf
MD}}}\bigg{(}\bm{a}^{\top}\bm{\widehat{u}}_{j}\sqrt{1+b_{j}^{{\sf
MD}}}-\bm{a}^{\top}\bm{u}_{j}\bm{u}_{j}^{\top}\bm{\widehat{u}}_{j}\sqrt{1+b_{j}^{{\sf
MD}}}$
$\displaystyle\qquad+\bm{a}^{\top}\bm{u}_{j}\big{(}\bm{u}_{j}^{\top}\bm{\widehat{u}}_{j}\sqrt{1+b_{j}^{{\sf
MD}}}-1\big{)}\bigg{)}\leq z\bigg{\\}}-\Phi(z)\bigg{|}$
$\displaystyle\leq\bigg{|}\mathbb{P}\bigg{\\{}\frac{1}{s_{\bm{a},j}^{{\sf
MD}}}\bigg{(}\bm{a}^{\top}\bm{\widehat{u}}_{j}\sqrt{1+b_{j}^{{\sf
MD}}}-\bm{a}^{\top}\bm{u}_{j}\bm{u}_{j}^{\top}\bm{\widehat{u}}_{j}\sqrt{1+b_{j}^{{\sf
MD}}}\bigg{)}\leq z\pm\bm{a}^{\top}\bm{u}_{j}{\sf
ErrBiasMD}\bigg{\\}}-\Phi(z)\bigg{|}$ $\displaystyle\quad+O(n^{-9})$
$\displaystyle\leq\bigg{|}\mathbb{P}\bigg{\\{}\frac{1}{s_{\bm{a},j}^{{\sf
MD}}}\bigg{(}\bm{a}^{\top}\bm{\widehat{u}}_{j}-\bm{a}^{\top}\bm{u}_{j}\bm{u}_{j}^{\top}\bm{\widehat{u}}_{j}\bigg{)}\leq\frac{z}{\sqrt{1+b_{j}^{{\sf
MD}}}}\pm\frac{\bm{a}^{\top}\bm{u}_{j}{\sf ErrBiasMD}}{\sqrt{1+b_{j}^{{\sf
MD}}}}\bigg{\\}}-\Phi(z)\bigg{|}+O(n^{-9})$
$\displaystyle\leq\bigg{|}\mathbb{P}\bigg{\\{}\frac{1}{s_{\bm{a},j}^{{\sf
MD}}}\bigg{(}\bm{a}^{\top}\bm{\widehat{u}}_{j}-\bm{a}^{\top}\bm{u}_{j}\bm{u}_{j}^{\top}\bm{\widehat{u}}_{j}\bigg{)}\leq\frac{z}{\sqrt{1+b_{j}^{{\sf
MD}}}}\pm\frac{\bm{a}^{\top}\bm{u}_{j}{\sf ErrBiasMD}}{s_{\bm{a},j}^{{\sf
MD}}\sqrt{1+b_{j}^{{\sf MD}}}}\bigg{\\}}$
$\displaystyle\quad-\Phi\bigg{(}\frac{z}{\sqrt{1+b_{j}^{{\sf
MD}}}}\pm\frac{\bm{a}^{\top}\bm{u}_{j}{\sf ErrBiasMD}}{s_{\bm{a},j}^{{\sf
MD}}\sqrt{1+b_{j}^{{\sf MD}}}}\bigg{)}\bigg{|}$
$\displaystyle\quad+\bigg{|}\Phi(\frac{z}{\sqrt{1+b_{j}^{{\sf
MD}}}})-\Phi\bigg{(}\frac{z}{\sqrt{1+b_{j}^{{\sf
MD}}}}\pm\frac{\bm{a}^{\top}\bm{u}_{j}{\sf ErrBiasMD}}{s_{\bm{a},j}^{{\sf
MD}}\sqrt{1+b_{j}^{{\sf MD}}}}\bigg{)}\bigg{|}+O(n^{-9})$
$\displaystyle\lesssim{\sf ErrMD}+\frac{|\bm{a}^{\top}\bm{u}_{j}|{\sf
ErrBiasMD}}{s_{\bm{a},j}^{{\sf MD}}\sqrt{1+b_{j}^{{\sf MD}}}}+n^{-9}$
$\displaystyle\lesssim{\sf ErrMD}+\frac{|\bm{a}^{\top}\bm{u}_{j}|{\sf
ErrBiasMD}}{s_{\bm{a},j}^{{\sf MD}}}+n^{-9},$

where we have applied Theorem 3.1 in the penultimate line. This completes the
proof of Theorem 3.2.

####  8.2.5 Step 5: Variance Approximation and Proof of Theorem 4.1

In order to prove Theorem 4.1 we need to demonstrate that the estimated
variance

$\displaystyle(\widehat{s}_{\bm{a},j}^{{\sf
MD}})^{2}:=2\widehat{\sigma}^{2}\sum_{k\leq r,k\neq
j}\frac{(\bm{\widehat{u}}_{k}^{\top}\bm{a})^{2}(1+\widehat{b_{k}^{{\sf
MD}}})}{(\widecheck{\lambda}_{k}-\widecheck{\lambda}_{j})^{2}}+2\frac{\widehat{\sigma}^{2}\|\bm{\widehat{U}}_{\perp}^{\top}\bm{a}\|^{2}}{\widecheck{\lambda}_{j}^{2}}$

yields a strong estimate of $(s_{\bm{a},j}^{{\sf MD}})^{2}$. This analysis
requires two steps. First we study the approximated bias $\widehat{b_{k}^{{\sf
MD}}}$.

######  Lemma 8.4.

Let $k\in[r]$ be fixed. Under the conditions of Theorem 4.1, with probability
at least $1-O(n^{-10})$ it holds that

$\displaystyle\bigg{|}\sqrt{1+\widehat{b_{k}^{{\sf MD}}}}-\sqrt{1+b_{k}^{{\sf
MD}}}\bigg{|}\lesssim\kappa\frac{\sigma^{2}\sqrt{rn}}{\lambda_{\min}^{2}}.$

###### Proof.

See Section B.4. ∎

As a consequence,

$\displaystyle\bigg{|}\frac{\bm{a}^{\top}\bm{\widehat{u}}_{j}\sqrt{1+\widehat{b_{j}^{{\sf
MD}}}}-\bm{a}^{\top}\bm{u}_{j}}{s^{{\sf
MD}}_{\bm{a},j}}-\frac{\bm{a}^{\top}\bm{\widehat{u}}_{j}\sqrt{1+b_{j}^{{\sf
MD}}}-\bm{a}^{\top}\bm{u}_{j}}{s^{{\sf MD}}_{\bm{a},j}}\bigg{|}$
$\displaystyle\lesssim\frac{|\bm{a}^{\top}\bm{\widehat{u}}_{j}|}{s^{{\sf
MD}}_{\bm{a},j}}\kappa\frac{\sigma^{2}\sqrt{rn}}{\lambda_{\min}^{2}}.$

Furthermore, we note that the analysis leading up to the proof of Theorem 3.2
implies that

$\displaystyle\big{|}\bm{a}^{\top}\bm{\widehat{u}}_{j}\big{|}\lesssim\big{|}\bm{a}^{\top}\bm{u}_{j}\big{|}+\sqrt{\log(n)}s^{{\sf
MD}}_{\bm{a},j}+{\sf ErrMD}\times s^{{\sf
MD}}_{\bm{a},j}\lesssim|\bm{a}^{\top}\bm{u}_{j}|+\sqrt{\log(n)}s^{{\sf
MD}}_{\bm{a},j}$

with probability at least $1-O(n^{-8})$, provided that ${\sf
ErrMD}\lesssim\sqrt{\log(n)}$, which holds by the assumptions in Theorem 4.1,
though we will verify this directly in Lemma 8.6. Therefore, with this same
probability,

$\displaystyle\bigg{|}\frac{\bm{a}^{\top}\bm{\widehat{u}}_{j}\sqrt{1+\widehat{b_{j}^{{\sf
MD}}}}-\bm{a}^{\top}\bm{u}_{j}}{s^{{\sf MD}}_{\bm{a},j}}$
$\displaystyle-\frac{\bm{a}^{\top}\bm{\widehat{u}}_{j}\sqrt{1+b_{j}^{{\sf
MD}}}-\bm{a}^{\top}\bm{u}_{j}}{s^{{\sf MD}}_{\bm{a},j}}\bigg{|}$
$\displaystyle\lesssim\frac{|\bm{a}^{\top}\bm{u}_{j}|}{s^{{\sf
MD}}_{\bm{a},j}}\kappa\frac{\sigma^{2}\sqrt{rn}}{\lambda_{\min}^{2}}+\sqrt{\log(n)}\kappa\frac{\sigma^{2}\sqrt{rn}}{\lambda_{\min}^{2}}.$
(10)

Next we study the approximated variance $(\widehat{s^{{\sf
MD}}_{\bm{a},j}})^{2}$, which is accomplished through the following lemma.

######  Lemma 8.5.

Instate the conditions in Theorem 4.1. If
$\|\bm{U}_{\perp}^{\top}\bm{a}\|<\|\bm{U}^{\top}\bm{a}\|$, then with
probability at least $1-O(n^{-8})$ it holds that

$\displaystyle\bigg{|}\frac{\widehat{s_{\bm{a},j}^{{\sf
MD}}}}{s_{\bm{a},j}^{{\sf MD}}}-1\bigg{|}$
$\displaystyle\lesssim\kappa\sqrt{\frac{r}{n}}+\max_{k\neq j}\frac{s^{{\sf
MD}}_{\bm{a},k}}{s^{{\sf
MD}}_{\bm{a},j}}\frac{\sigma\sqrt{r\log(n)}}{\Delta_{j}}+\frac{\sigma^{2}r\log(n)}{\Delta_{\min}^{2}}+\frac{\sigma^{2}\kappa\sqrt{rn\log(n)}}{\lambda_{\min}^{2}}$
$\displaystyle\quad+\frac{\sigma\sqrt{n}}{\lambda_{\min}}\kappa^{2}\max_{k\neq
j,k\leq
r}\frac{(\bm{a}^{\top}\bm{u}_{j})^{2}}{(\bm{a}^{\top}\bm{u}_{k})^{2}}.$

If instead $\|\bm{U}^{\top}\bm{a}\|\leq\|\bm{U}_{\perp}^{\top}\bm{a}\|$, then
with this same probability,

$\displaystyle\bigg{|}$ $\displaystyle\frac{\widehat{s_{\bm{a},j}^{{\sf
MD}}}}{s_{\bm{a},j}^{{\sf
MD}}}-1\bigg{|}\lesssim\kappa\sqrt{\frac{r}{n}}+\max_{k\neq j}\frac{s^{{\sf
MD}}_{\bm{a},k}}{s^{{\sf
MD}}_{\bm{a},j}}\frac{\sigma\sqrt{r\log(n)}}{\Delta_{j}}+\frac{\sigma^{2}r\log(n)}{\Delta_{\min}^{2}}+\frac{\sigma^{2}\kappa\sqrt{rn\log(n)}}{\lambda_{\min}^{2}}+\frac{\sigma\sqrt{n}}{\lambda_{\min}}.$

###### Proof.

See Section B.5. ∎

Let the error in Lemma 8.5 be denoted as ${\sf ErrCIMD}$. As a consequence of
this lemma, with probability at least $1-O(n^{-8})$,

$\displaystyle\bigg{|}$
$\displaystyle\frac{\bm{a}^{\top}\bm{\widehat{u}}_{j}\sqrt{1+\widehat{b_{j}^{{\sf
MD}}}}-\bm{a}^{\top}\bm{u}_{j}}{\widehat{s^{{\sf
MD}}_{\bm{a},j}}}-\frac{\bm{a}^{\top}\bm{\widehat{u}}_{j}\sqrt{1+b_{j}^{{\sf
MD}}}-\bm{a}^{\top}\bm{u}_{j}}{s^{{\sf MD}}_{\bm{a},j}}\bigg{|}$
$\displaystyle\lesssim\bigg{|}\frac{\bm{a}^{\top}\bm{\widehat{u}}_{j}\sqrt{1+\widehat{b_{j}^{{\sf
MD}}}}-\bm{a}^{\top}\bm{u}_{j}}{s^{{\sf
MD}}_{\bm{a},j}}-\frac{\bm{a}^{\top}\bm{\widehat{u}}_{j}\sqrt{1+b_{j}^{{\sf
MD}}}-\bm{a}^{\top}\bm{u}_{j}}{s^{{\sf MD}}_{\bm{a},j}}\bigg{|}$