# Sampling without replacement from a high-dimensional finite population

Jiang<EMAIL_ADDRESS>[ Shaochen
<EMAIL_ADDRESS>[ Yangchun
<EMAIL_ADDRESS>Corresponding author[ Wang
<EMAIL_ADDRESS>[[ url]  KLASMOE and School of
Mathematics $\&$ Statistics, Northeast Normal University, Changchun, 130024,
P.R. China.  School of Mathematics, South China University of Technology,
Guangzhou, 510640, P.R. China.  Department of Mathematics, Shanghai
University, Shanghai, 200444, P.R. China.  Department of Statistics and Data
Science, National University of Singapore, 117546, Singapore.  Some University
and Another University

###### Abstract

It is well known that most of the existing theoretical results in statistics
are based on the assumption that the sample is generated with replacement from
an infinite population. However, in practice, available samples are almost
always collected without replacement. If the population is a finite set of
real numbers, whether we can still safely use the results from samples drawn
without replacement becomes an important problem. In this paper, we focus on
the eigenvalues of high-dimensional sample covariance matrices generated
without replacement from finite populations. Specifically, we derive the
Tracy-Widom laws for their largest eigenvalues and apply these results to
parallel analysis. We provide new insight into the permutation methods
proposed by Buja and Eyuboglu in [Multivar Behav Res. 27(4) (1992) 509–540].
Simulation and real data studies are conducted to demonstrate our results.

Largest eigenvalue; Tracy-Widom law; sample covariance matrix; finite
population model; parallel analysis,

###### keywords:

and

##  1 Introduction

Let ${\mathbf{U}}_{p\times
N}=({\bm{u}}^{\top}_{1},\ldots,{\bm{u}}^{\top}_{p})^{\top}$ be a finite
population matrix of dimension $p$ and size $N$, where
${\bm{u}}_{i}=(u_{i1},\ldots,u_{iN})$, $1\leq i\leq p$, is the i-th row of
${\mathbf{U}}_{p\times N}$. Suppose that for each $i=1,\ldots,p$,
${{\bm{w}}}_{i}=(w_{i1},\ldots,w_{in})$ is generated by sampling without
replacement from the $N$ coordinates of ${\bm{u}}_{i}$ with equal probability.
Then, the sample covariance matrix from the finite population is given by

$\displaystyle{{\bm{S}}}_{n}=\frac{1}{n}{\bm{W}}_{n}{\bm{W}}_{n}^{\top},$
(1.1)

where
${\bm{W}}_{n}=({{\bm{w}}}_{1}^{\top},\ldots,{{\bm{w}}}_{p}^{\top})^{\top}$. In
this paper, our aim is to derive the asymptotic properties of the eigenvalues
of ${{\bm{S}}}_{n}$ in the large dimensional framework, that is, cases where
the dimension and sample size proportionally tend to infinity.

Sample covariance matrices have been extensively investigated since proposed
by Wishart (1928). For fixed sample dimension cases, one can see any
multivariate statistical analysis (MSA) textbook (e.g., (Anderson, 1984)).
However, when the dimension of the sample is large or even larger than the
sample size, these classical theories may no longer be applicable. Since
Marcenko and Pastur discovered the well-known global spectral distribution (MP
law) in their seminal work Marcenko and Pastur (1967), there has been much
literature devoted to the spectral property of the large-dimensional sample
covariance matrix and its variants. For the global properties of eigenvalues
of sample covariance matrices, see Bai and Silverstein (2010), where one can
also find many references related to this topic. For the individual
eigenvalue, the most interesting one is the largest eigenvalue, since it
accounts for as much of the variability in data as possible. The limiting
distribution of the largest eigenvalue of the large-dimensional sample
covariance matrix was first studied and shown to follow the Tracy-Widom (TW)
distribution by Johansson (2000) (complex Gaussian populations) and Johnstone
(2001) (real Gaussian populations) under the assumption that all components of
the population are independent and identically distributed (i.i.d.). For the
follow-up work, see El Karoui (2007); Onatski (2008); Bao et al. (2015); Lee
and Schnelli (2016).

The coordinates of the population are linear combinations of independent
random variables, which is a common condition in the above results. To relax
this condition, many researchers have considered more general cases. For
example, Bai and Zhou (2008) considered the global spectral distribution of
the sample covariance matrix without an independence structure in columns, Bao
et al. (2012) and Pillai and Yin (2012) considered the TW law of the largest
eigenvalue of the sample correlation matrix, Gao et al. (2017) considered the
linear spectral statistics of the sample correlation matrix, Bao (2019)
considered the TW law of the largest eigenvalue of the Kendall rank
correlation matrix, Hu et al. (2019) and Wen et al. (2020) considered the
linear spectral statistics and the largest eigenvalue of the sample covariance
matrix, respectively, in elliptical distributions, and Tian et al. (2021)
considered the largest eigenvalue of the sample covariance matrix generated by
VARMA.

It is worth noting that the sample covariance matrix from finite population
${{\bm{S}}}_{n}$, which is a classical statistic in MSA, is a kind of sample
covariance matrix without an independence structure. However, it does not
belong to any model presented in the above references. This is because (i) the
entries of each row of ${\bm{W}}_{n}$ are correlated and (ii) the populations
of different rows of ${\bm{W}}_{n}$ could be different. Here we introduce
Spearman’s rank correlation matrix and the permutation matrix in parallel
analysis (PA) as two examples of our matrix model.

  * (I):

Spearman’s rank correlation coefficient is a classical nonparametric measure
of correlation between two random variables. Suppose that we have $n$ i.i.d.
samples ${\bm{Y}}_{1},\ldots,{\bm{Y}}_{n}$, where ${\bm{Y}}_{k}$ is a
$p$-vector consisting of i.i.d. random variables
$Y_{1k},\ldots,Y_{pk},k=1,\ldots,n.$ Spearman’s rank correlation matrix is
defined by

$\displaystyle\mathfrak{R}_{p}=\left(\mathfrak{r}_{kl}\right)_{1\leq k,l\leq
p},$ (1.2)

where $\mathfrak{r}_{kl}$ is Spearman’s rank correlation coefficient between
the $k$-th and $l$-th rows of $\left[{\bm{Y}}_{1},\ldots,{\bm{Y}}_{n}\right]$.
Actually, $\mathfrak{R}_{p}$ can be expressed in the form of (1.1) with
$u_{ij}=\frac{\sqrt{12}}{\sqrt{n^{2}-1}}(j-\frac{n+1}{2})$, $i=1,\dots,p$,
$j=1,\dots,N$ and $n=N$. See Bai and Zhou (2008); Bao et al. (2015); Bao
(2017) for more details.

  * (II):

PA is a popular method that is used to select the number of factors in the
factor model. PA is an applicable method because the largest singular value of
the permuted data matrix decreases compared with that of the original value if
there is indeed at least one factor (see Dobriban and Owen (2019); Dobriban
(2020) for further explanation). Suppose the sample matrix
${\bm{B}}=(b_{ij})_{p\times n}$ satisfies the signal-plus-noise model, i.e.,
${\bm{B}}={\bm{S}}+\bm{\Psi}$, where ${\bm{S}}$ is a $p\times n$ signal matrix
of low rank and $\bm{\Psi}$ is a $p\times n$ noise matrix. In this case, the
permutation matrix in Dobriban and Owen (2019); Dobriban (2020) turns into
${\bm{W}}$ of (1.1) with $N=n$.

We highlight two main contributions of the present paper. First, we
theoretically investigate some spectral properties of high-dimensional
covariance matrices from finite populations, which have not attracted much
attention in modern statistics. Second, we practically characterize the
asymptotic distribution of the largest eigenvalue of the high-dimensional
standard permutation matrix in PA, which provides a direct procedure without
permutation.

The remainder of this paper is organized as follows. In Section 2, we
introduce the necessary notations and state the main results, which include
the global law, the edge universality and the fluctuation of the largest
eigenvalue of the sample covariance matrix ${{\bm{S}}}_{n}$. In Section 3, we
introduce some applications of our main results in high-dimensional
statistical inference to parallel analysis. Related simulations and real data
analyses are also conducted to demonstrate our results. Proofs of the main
theorems and some preparatory lemmas are given in Section 4.

##  2 Notations and statements of main results

###  2.1 Main results

To easily achieve our goal, we write ${\mathbf{U}}_{p\times N}$ as

$\displaystyle{\mathbf{U}}_{p\times
N}={\bm{T}}^{\frac{1}{2}}{\mathbf{V}}_{p\times N},$

where ${\bm{T}}=\text{{\rm diag}}\\{t_{1},\ldots,t_{p}\\}$ is a positive
definite $p$-dimensional diagonal matrix and ${\mathbf{V}}_{p\times
N}=(v_{jk})_{1\leq j\leq p,~{}1\leq k\leq N}$ is a $p\times N$ matrix.
${\bm{T}}$ is chosen such that the rows of ${\mathbf{V}}_{p\times N}$ are
standardized. The rows of ${\mathbf{U}}_{p\times N}$ and
${\mathbf{V}}_{p\times N}$ are denoted as ${\bm{u}}_{1},\ldots,{\bm{u}}_{p}$
and ${\bm{v}}_{1},\ldots,{\bm{v}}_{p}$, respectively. After scaling the
observations
${{\bm{y}}}_{i}={{\bm{w}}}_{i}/\sqrt{n},~{}{{\bm{x}}}_{i}={{\bm{y}}}_{i}/\sqrt{t_{i}},~{}i=1,\ldots,p$
and setting
${\bm{X}}_{n}=({{\bm{x}}}_{1}^{\top},\ldots,{{\bm{x}}}_{p}^{\top})^{\top}=(x_{jk})_{1\leq
j\leq p,~{}1\leq k\leq n}$,
${\bm{Y}}_{n}=({{\bm{y}}}_{1}^{\top},\ldots,{{\bm{y}}}_{p}^{\top})^{\top}=(y_{jk})_{1\leq
j\leq p,~{}1\leq k\leq n}$, we can rewrite (1.1) as

$\displaystyle{{\bm{S}}}_{n}={\bm{Y}}_{n}{\bm{Y}}_{n}^{\top}$
$\displaystyle={\bm{T}}^{\frac{1}{2}}{\bm{X}}_{n}{\bm{X}}_{n}^{\top}{\bm{T}}^{\frac{1}{2}}.$
(2.1)

Henceforth, we denote
$\lambda_{p}({\bm{A}})\leq\cdots\leq\lambda_{1}({\bm{A}})$ as the ordered
eigenvalues of a $p\times p$ Hermitian matrix ${\bm{A}}$. The empirical
spectral distribution (ESD) of ${\bm{T}}$ is

$\displaystyle
H_{n}(\lambda):=\frac{1}{p}\sum_{j=1}^{p}{\mathbb{I}}_{\\{\lambda_{j}({\bm{T}})\leq\lambda\\}},\quad\lambda\in\mathbb{R}$
(2.2)

and that of ${{\bm{S}}}_{n}$ is

$\displaystyle
F_{n}(\lambda):=\frac{1}{p}\sum_{j=1}^{p}{\mathbb{I}}_{\\{\lambda_{j}({{\bm{S}}}_{n})\leq\lambda\\}},\quad\lambda\in\mathbb{R}.$

Here and throughout the paper, ${\mathbb{I}}_{A}$ represents the indicator
function of event $A$. The Stieltjes transform of $F_{n}$ is given by

$\displaystyle m_{n}(z)=\int\frac{1}{x-z}{\rm d}F_{n}(x),$

where $z=E+i\eta\in\mathbb{C}^{+}$.

In addition, we need a crucial parameter $\xi_{+}$ that satisfies

$\displaystyle\int\left(\frac{\lambda\xi_{+}}{1-\lambda\xi_{+}}\right)^{2}{\rm
d}H_{n}(\lambda)=p/n.$ (2.3)

It is easy to check that the solution to (2.3) in
$[0,1/\lambda_{1}({\bm{T}}))$ is unique. See El Karoui (2007) for more
discussion of $\xi_{+}$. Now, we state our main assumptions as follows.

######  Assumption 2.1.

  * (i)

(On dimensionality) We assume

$\displaystyle{c_{n}:=p/n\rightarrow
c_{0}\in(0,\infty),~{}~{}y_{n}:=n/N\rightarrow y\in[0,1]}$ (2.4)

as $n\to\infty$. For simplicity, we assume that there are two positive
constants $c_{1}$, $C_{1}$ such that $c_{1}<c_{n}<C_{1}$ for all $n\geq 1$.

  * (ii)

(On ${\bm{T}}$) We assume that ${\bm{T}}$ is diagonal with
${\liminf}_{n}\lambda_{p}({\bm{T}})>0$,
${\limsup}_{n}\lambda_{1}({\bm{T}})<\infty$,

$\displaystyle\limsup_{n}\lambda_{1}({\bm{T}})\xi_{+}<1,$ (2.5)

and there exists a deterministic distribution $H$ such that $H_{n}\to H$ as
$n\to\infty$.

  * (iii)

(On ${\bm{V}}$) We assume that the elements in ${\mathbf{V}}_{p\times N}$
satisfy for each $j=1,\ldots,p$, $\sum_{k=1}^{N}v_{jk}=0$,
$N^{-1}\sum_{k=1}^{N}v_{jk}^{2}=1$. In addition, for any positive integer $m$,
there is a constant $C_{m}$ such that

$\displaystyle\frac{1}{N}\sum_{k=1}^{N}|v_{jk}|^{m}\leq C_{m},\ j=1,\ldots,p.$
(2.6)

######  Remark 2.1.

We remark that the aim of Assumption (2.6) is to make moments of
$\sqrt{n}x_{ij}$ bounded, i.e., $\mathbb{E}|\sqrt{n}x_{ij}|^{m}\leq C_{m}$,
which justifies the large deviation of the strong MP law and the moments’
comparison procedure. Details can be found in Section 4.

Our main results on the sample covariance matrix ${{\bm{S}}}_{n}$ from a
finite population are formulated as follows. First, we present the global law
of ${{\bm{S}}}_{n}$.

######  Theorem 2.2.

Suppose that ${{\bm{S}}}_{n}$ satisfies Assumption 2.1. Then, as $n\to\infty$,
$F_{n}$ almost surely converges to a probability distribution $F_{c_{0},H}$,
whose Stieltjes transform $m_{c_{0},H}=m_{c_{0},H}(z)$ is determined by

$\displaystyle m_{c_{0},H}=\int\frac{1}{t(1-c_{0}-c_{0}zm_{c_{0},H})-z}{\rm
d}H(t).$ (2.7)

We further denote by
${\bm{Q}}:=n^{-1}{\bm{T}}^{\frac{1}{2}}\mathcal{X}\mathcal{X}^{\top}{\bm{T}}^{\frac{1}{2}}$,
where $\mathcal{X}$ is a $p\times n$ data matrix with i.i.d. $N(0,1)$
variables. Our main theorem on the edge universality of ${{\bm{S}}}_{n}$ can
be formulated as follows, and its proof is deferred to Section 4. Let

$\displaystyle
E_{+}=\frac{1}{\xi_{+}}\left(1+c_{n}^{-1}\int\frac{\lambda\xi_{+}}{1-\lambda\xi_{+}}{\rm
d}H_{n}(\lambda)\right).$ (2.8)

######  Theorem 2.3.

Suppose that Assumption 2.1 holds, and $N=n$. Then, there exist positive
constants $\varepsilon$ and $\delta$ such that for any $s\in\mathbb{R}$,

$\displaystyle\mathbb{P}\left(n^{2/3}(\lambda_{1}({\bm{Q}})-E_{+})\leq
s-n^{-\varepsilon}\right)-n^{-\delta}$ $\displaystyle\leq$
$\displaystyle\mathbb{P}\left(n^{2/3}(\lambda_{1}({{\bm{S}}}_{n})-E_{+})\leq
s\right)$ $\displaystyle\leq$
$\displaystyle\mathbb{P}\left(n^{2/3}(\lambda_{1}({\bm{Q}})-E_{+})\leq
s+n^{-\varepsilon}\right)+n^{-\delta},$

holds for sufficiently large $n$.

One typical application of Theorems 2.2 and 2.3 lies in Spearman’s correlation
matrix (1.2). In this case, ${\bm{T}}$ is identity, and it is straightforward
to verify that Assumption 2.1 holds. Thus, we have the following corollary,
which was first derived in Bao (2017).

######  Corollary 2.4.

Suppose that $c_{n}:=p/n\rightarrow c_{0}\in(0,\infty)$. Then, the ESD of
Spearman’s rank correlation matrix $\mathfrak{R}_{p}$ almost surely converges
to the standard MP law, and

$\displaystyle
n^{\frac{2}{3}}\left(\frac{\lambda_{1}(\mathfrak{R}_{p})-\underline{E}_{+}}{\underline{\gamma}_{0}}\right)\Rightarrow
TW_{1},$

where $\underline{E}_{+}=(1+\sqrt{c_{n}})^{2}$ and
$\underline{\gamma}_{0}=c_{n}^{-\frac{1}{6}}\underline{E}_{+}^{\frac{2}{3}}$.
Here and in the sequel , $TW_{1}$ is the type-1 TW distribution,
$``\Rightarrow"$ means convergence in distribution.

However, one cannot apply Theorem 2.3 directly in PA because in PA, the data
matrix is random, and consequently, the condition $\sum_{k=1}^{N}v_{jk}=0$ is
not satisfied. To overcome this problem, we consider modified versions of
${{\bm{S}}}_{n}$, which fortunately work for all cases of $n\leq N$.

Now, we consider the general case of $n\leq N$. We first introduce the spatial
sign matrix. Let $\bar{{\bm{Y}}}={\rm
diag}\\{\frac{1}{n}\sum_{l=1}^{n}y_{1l},\dots,\frac{1}{n}\sum_{l=1}^{n}y_{pl}\\}{\bm{e}}_{p}{\bm{e}}_{p}^{\top}$,
${\bm{S}}_{P}=({\bm{Y}}_{n}-\bar{{\bm{Y}}})({\bm{Y}}_{n}-\bar{{\bm{Y}}})^{\top}$
and ${\bm{D}}_{2}^{(n)}={\rm
diag}\\{({\bm{S}}_{P})_{11},\dots,({\bm{S}}_{P})_{pp}\\}$. Here ${\bm{e}}_{p}$
represents the $p$-dimensional vector with elements all 1. Then, the spatial
sign matrix is defined as

$\displaystyle{\bm{P}}_{n}=({\bm{D}}_{2}^{(n)})^{-\frac{1}{2}}{\bm{S}}_{P}({\bm{D}}_{2}^{(n)})^{-\frac{1}{2}}.$
(2.9)

In addition, we provide an alternative assumption for (iii) of Assumption 2.1;
namely,

  * (iii’)

Assume that for any positive integer $m$, there is a constant $C_{m}$ such
that

$\displaystyle\frac{1}{N}\sum_{k=1}^{N}|v_{jk}|^{m}\leq
C_{m},~{}~{}~{}~{}j=1,\dots,p.$ (2.10)

Furthermore, we assume

$\displaystyle\frac{1}{N}\sum_{j=1}^{N}\Big{(}v_{ij}-\frac{1}{N}\sum_{k=1}^{N}v_{ik}\Big{)}^{2}\geq
c>0,$ (2.11)

for some constant $c$.

######  Remark 2.5.

Note that $\sum_{k=1}^{N}v_{jk}=0$ and $N^{-1}\sum_{k=1}^{N}v_{jk}^{2}=1$ in
(iii) are no longer needed because of the standardization in the spatial-sign
matrix. The aim of the new condition (2.11) is to bound the moments of
$\sqrt{n}\frac{y_{ij}-\frac{1}{n}\sum_{k=1}^{n}y_{ik}}{\sqrt{({\bm{S}}_{P})_{ii}}}$,
i.e.,

$\displaystyle\mathbb{E}\Big{|}\sqrt{n}\frac{y_{ij}-\frac{1}{n}\sum_{k=1}^{n}y_{ik}}{\sqrt{({\bm{S}}_{P})_{ii}}}\Big{|}^{m}\leq
C_{m},$

which will be used in the proof of Theorem 2.6.

######  Theorem 2.6.

In addition to (i) and (ii) of Assumption 2.1, we assume (iii’) holds. We
write ${\bm{Q}}_{1}:=\frac{1}{n}\mathcal{X}\mathcal{X}^{\top}$; then, there
exist positive constants $\varepsilon$ and $\delta$ such that for any
$s\in\mathbb{R}$,

$\displaystyle\mathbb{P}\left(n^{2/3}(\lambda_{1}({\bm{Q}}_{1})-\underline{E}_{+})\leq
s-n^{-\varepsilon}\right)-n^{-\delta}$ $\displaystyle\leq$
$\displaystyle\mathbb{P}\left(n^{2/3}(\lambda_{1}({\bm{P}}_{n})-\underline{E}_{+})\leq
s\right)$ $\displaystyle\leq$
$\displaystyle\mathbb{P}\left(n^{2/3}(\lambda_{1}({\bm{Q}}_{1})-\underline{E}_{+})\leq
s+n^{-\varepsilon}\right)+n^{-\delta}$

holds for sufficiently large $n$. Here, $\underline{E}_{+}$ is defined in
Corollary 2.4.

######  Remark 2.7.

The above results can be generalized to the joint distribution of the largest
several eigenvalues. More specifically, under the assumptions in Theorem 2.3,
there exist positive constants $\varepsilon$ and $\delta$ such that for any
$s_{1},\ldots,s_{k}\in\mathbb{R}$,

$\displaystyle\mathbb{P}\left(n^{2/3}(\lambda_{1}({\bm{Q}})-E_{+})\leq
s_{1}-n^{-\varepsilon},\ldots,n^{2/3}(\lambda_{k}({\bm{Q}})-E_{+})\leq
s_{k}-n^{-\varepsilon}\right)-n^{-\delta}$ $\displaystyle\leq$
$\displaystyle\mathbb{P}\left(n^{2/3}(\lambda_{1}({{\bm{S}}}_{n})-E_{+})\leq
s_{1},\ldots,n^{2/3}(\lambda_{k}({{\bm{S}}}_{n})-E_{+})\leq s_{k}\right)$
$\displaystyle\leq$
$\displaystyle\mathbb{P}\left(n^{2/3}(\lambda_{1}({\bm{Q}})-E_{+})\leq
s_{1}+n^{-\varepsilon},\ldots,n^{2/3}(\lambda_{1}({\bm{Q}})-E_{+})\leq
s_{k}+n^{-\varepsilon}\right)+n^{-\delta}$

holds for sufficiently large $n$. Correspondingly,
$\lambda_{1}({\bm{P}}_{n}),\ldots,\lambda_{k}({\bm{P}}_{n})$ have analogous
properties under the conditions of Theorem 2.6.

From Theorems 2.3 and 2.6, we have the following corollary on the largest
eigenvalues of ${{\bm{S}}}_{n},{\bm{P}}_{n}$.

######  Corollary 2.8.

Suppose that Assumption 2.1 holds. Denoting

$\displaystyle\gamma_{0}^{3}=\frac{1}{\xi_{+}^{3}}\left(1+c_{n}^{-1}\int\left(\frac{\lambda\xi_{+}}{1-\lambda\xi_{+}}\right)^{3}{\rm
d}H_{n}(\lambda)\right),$ (2.12)

we have

$\displaystyle
n^{\frac{2}{3}}\left(\frac{\lambda_{1}({{\bm{S}}}_{n})-E_{+}}{\gamma_{0}}\right)\Rightarrow
TW_{1}.$

Additionally, under the same assumptions in Theorem 2.6, we have

$\displaystyle
n^{\frac{2}{3}}\left(\frac{\lambda_{1}({\bm{P}}_{n})-\underline{E}_{+}}{\underline{\gamma}_{0}}\right)\Rightarrow
TW_{1}.$

The proof of Theorem 2.3 is deferred to Section 4. The proofs of Theorem 2.6
and Corollary 2.8 are deferred to Section A.3 of Supplementary Material (Hu et
al. (2022)).

###  2.2 Basic notions

What truly pertains to our discussion in the sequel is the nonasymptotic
version of $F_{c_{0},H}$ (i.e., $F_{c_{n},H_{n}}$), which can be obtained by
replacing $c_{0}$ and $H$ with $c_{n}$ and $H_{n}$ in $F_{c_{0},H}$,
respectively. More precisely, $F_{c_{n},H_{n}}$ is the corresponding
distribution function of the Stieltjes transform
$m_{c_{n},H_{n}}(z):=m_{c_{n},H_{n}}\in\mathbb{C}^{+}$ satisfying the
following self-consistent equation:

$\displaystyle
m_{c_{n},H_{n}}(z)=\int\frac{1}{t(1-c_{n}-c_{n}zm_{c_{n},H_{n}}(z))-z}{\rm
d}H_{n}(t).$ (2.13)

Define an $n\times n$ matrix

$\displaystyle\mathcal{S}_{n}={\bm{Y}}_{n}^{\top}{\bm{Y}}_{n},$ (2.14)

which shares the same nonzero eigenvalues with ${\bm{S}}_{n}$. Denoting the
ESD of $\mathcal{S}_{n}$ as $\underline{F}_{n}$, we have

$\underline{F}_{n}=c_{n}F_{n}+(1-c_{n})\mathbb{I}_{[0,\infty)},$

which converges almost surely to

$\displaystyle\underline{F}_{c_{0},H}=c_{0}F_{c_{0},H}+(1-c_{0})\mathbb{I}_{[0,\infty)}.$
(2.15)

Then, the n-dependent version of (2.15) is

$\displaystyle\underline{F}_{c_{n},H_{n}}=c_{n}F_{c_{n},H_{n}}+(1-c_{n})\mathbb{I}_{[0,\infty)}.$
(2.16)

We also obtain

$\displaystyle\underline{m}_{c_{n},H_{n}}=c_{n}m_{c_{n},H_{n}}-\frac{1-c_{n}}{z},$
(2.17)

which satisfies the following self-consistent equation:

$\displaystyle
z=-\frac{1}{\underline{m}_{c_{n},H_{n}}}+c_{n}\int\frac{t}{1+t\underline{m}_{c_{n},H_{n}}}\text{d}H_{n}(t).$
(2.18)

For simplicity, we briefly use the notation

$\displaystyle\tilde{m}(z):=m_{c_{n},H_{n}}(z),~{}~{}\tilde{F}:=F_{c_{n},H_{n}},$
$\displaystyle\tilde{\underline{m}}(z):=\underline{m}_{c_{n},H_{n}}(z),~{}~{}\underline{\tilde{F}}:=\underline{F}_{c_{n},H_{n}}.$

From the discussions in Silverstein and Choi (1995), we know that
$\underline{F}_{c_{n},H_{n}}$ has a continuous derivative
$\underline{\rho}_{0}$ on $\mathbb{R}\backslash\\{0\\}$. The inequality (2.5)
guarantees that the density $\underline{\rho}_{0}$ exhibits a square-root-type
behavior at the rightmost endpoint of its support; see Theorem 3.1 in Bao et
al. (2015). The rightmost boundary of the support of $\underline{\rho}_{0}$ is
$E_{+}$ defined in (2.8), that is,
$E_{+}=\inf\\{x\in\mathbb{R}:\underline{\tilde{F}}(x)=1\\}$. Moreover, the
parameter $\xi_{+}$ defined by (2.3) satisfies $\xi_{+}=-\lim_{\eta\to
0}\tilde{\underline{m}}(E_{+}+i\eta)=-\tilde{\underline{m}}(E_{+})$.

##  3 Applications and numerical studies

###  3.1 High-dimensional PA

As Spearman’s rank correlation matrix has been discussed in the unpublished
paper by Bao (2017), in this section, we only focus on the high-dimensional
PA. Recall the signal-plus-noise matrix ${\bm{B}}={\bm{S}}+\bm{\Psi}$, and
recall that $\bm{\Psi}$ has the following structure:

$\bm{\Psi}=\tilde{{\bm{T}}}^{\frac{1}{2}}\bm{\Phi},$

where $\tilde{{\bm{T}}}={\rm diag}\\{\tilde{t}_{1},\ldots,\tilde{t}_{p}\\}$ is
a diagonal matrix with $\|\tilde{{\bm{T}}}\|$ and $\|\tilde{{\bm{T}}}^{-1}\|$
bounded in $p$, $\bm{\Phi}$ is a $p\times n$ collection of i.i.d. noise
$\phi_{ij}$ satisfying

$\displaystyle{\mathbb{E}}\phi_{ij}=0,~{}~{}{\mathbb{E}}\phi_{ij}^{2}=1.$
(3.1)

Moreover, we assume that every $\phi_{ij}$ has a uniform subexponential decay,
that is, there exists a constant $\vartheta>0$ such that for $u>1$

$\displaystyle\mathbb{P}(|\phi_{ij}|>u)\leq\vartheta^{-1}\exp(-u^{\vartheta}).$
(3.2)

Denote $\bm{S}=(s_{ij})_{p\times n}$, and assume that for any positive integer
$m$, there exists a constant $C_{m}$ such that

$\displaystyle\mathbb{E}|s_{ij}|^{m}\leq C_{m},~{}~{}1\leq i\leq p,~{}1\leq
j\leq n.$ (3.3)

Recall the standard factor model
${\bm{B}}=\bm{H}\bm{\Lambda}^{\top}+\bm{\Psi}$, where $\bm{H}$ is the $p\times
r$ ($r$ is fixed) matrix containing the factor values $h_{ij}$ and
$\bm{\Lambda}$ is the $n\times r$ factor loading matrix with entries
$\lambda_{jk}$. The first term is the signal due to the factor component,
whose rank is at most $r$. Thus, the factor model falls into the signal-plus-
noise framework.

We start with a $p\times n$ data matrix ${\bm{B}}$ with measurements given by
$b_{ij}$, $i=1,\ldots,p,~{}j=1,\ldots,n$. We generate a matrix
${\bm{B}}_{\pi}$ by separately permuting the entries in each row of
${\bm{B}}$. Here $\pi=(\pi_{1},\pi_{2},\ldots,\pi_{p})$ is a permutation
array, which is a collection of permutations $\pi_{i}$ of $1,\ldots,n$. The
permutation $\pi_{i}$ permutes the $i$-th row of ${\bm{B}}$; hence,
${\bm{B}}_{\pi}$ has $({\bm{B}}_{\pi})_{ij}={\bm{B}}_{i\pi_{i}(j)}$ entries.

Buja and Eyuboglu (1992) suggested selecting the first factor if the top
singular value of ${\bm{B}}$ is larger than a fixed percentile of the top
singular values of the permuted matrices. Here, one can use the median or the
$95\%$ or $100\%$ percentiles. If the first factor is selected, then we select
the second factor when the second largest singular value of ${\bm{B}}$ is
larger than the second singular value of the permuted matrices. We repeat the
same procedure until a factor is not selected.

We try to calculate the quantiles theoretically. Based on the testing strategy
introduced by Dobriban (2020), we consider an analogous strategy to compare
the top singular value of ${\bm{B}}$ and a certain percentile of the top
singular values of the permuted matrices $(\hat{{\bm{B}}}_{1})_{\pi}$ or
$(\hat{{\bm{B}}}_{2})_{\pi}$, where

$\hat{{\bm{B}}_{1}}:={\bm{B}}-\bar{{\bm{B}}}~{}~{}\text{and}~{}~{}\hat{{\bm{B}}_{2}}:=({\bm{D}}_{2}^{(n)})^{-1/2}({\bm{B}}-\bar{{\bm{B}}}).$

$\bar{{\bm{B}}}$ and ${\bm{D}}_{2}^{(n)}$ are obtained in the same way as
similar matrices in the construction of the spatial sign matrix ${\bm{P}}_{n}$
in (2.9).

Let $\varrho_{n}(\hat{{\bm{B}}}_{\ell})$ be the largest eigenvalue of
$\frac{1}{n}(\hat{{\bm{B}}}_{\ell})_{\pi}(\hat{{\bm{B}}}_{\ell})_{\pi}^{\top}$
for a random permutation $\pi$, $\ell=1,2$. As the number of all the
permutations of $(\hat{{\bm{B}}}_{\ell})_{\pi}$ is $(n!)^{p}$, the probability
mass function of $\varrho_{n}(\hat{{\bm{B}}}_{\ell})$ given ${\bm{B}}$ is

$\mathbb{P}\left(\varrho_{n}(\hat{{\bm{B}}}_{\ell})=\\{\lambda_{1}(\frac{1}{n}(\hat{{\bm{B}}}_{\ell})_{\pi}(\hat{{\bm{B}}}_{\ell})_{\pi}^{\top}),\pi\in\Pi(\hat{{\bm{B}}}_{\ell})\\}_{k}|{\bm{B}}\right)=\frac{1}{(n!)^{p}},~{}k=1,\ldots,(n!)^{p},~{}\ell=1,2{\color[rgb]{1,0,0},}$

where $\Pi(\hat{{\bm{B}}}_{\ell})$ represents the set of all permutations of
$(\hat{{\bm{B}}}_{\ell})_{\pi}$ and $\\{A\\}_{k}$ denotes the $k$-th element
in a set $A$. For $\varrho_{n}(\hat{{\bm{B}}}_{\ell})$, we have the following
corollary, whose proof is deferred to the Supplementary Material.

######  Corollary 3.1.

Suppose that (i) and (ii) of Assumption 2.1 hold when ${\bm{T}}$ is replaced
by $({\bm{D}}_{2}^{(n)})^{-1}$. Moreover, assume that (3.1), (3.2) and (3.3)
hold. Then, we have that

$\displaystyle
n^{\frac{2}{3}}\left(\frac{\varrho_{n}(\hat{{\bm{B}}}_{2})-\underline{E}_{+}}{\underline{\gamma}_{0}}\right)\Rightarrow
TW_{1},$ (3.4)

and

$\displaystyle
n^{\frac{2}{3}}\left(\frac{\varrho_{n}(\hat{{\bm{B}}}_{1})-E_{+}}{\gamma_{0}}\right)\Rightarrow
TW_{1},$ (3.5)

where $\underline{E}_{+}$ and $\underline{\gamma}_{0}$ are defined in
Corollary 2.4, $E_{+}$ and ${\gamma_{0}}$ are obtained by (2.8) and (2.12),
respectively, when ${\bm{T}}$ is replaced by $({\bm{D}}_{2}^{(n)})^{-1}$.

######  Remark 3.2.

In simulation, $\underline{E}_{+}$ and $\underline{\gamma}_{0}$ are reassigned
to $((\sqrt{n-1}+\sqrt{p})/\sqrt{n})^{2}$ and
$(\sqrt{n-1}+\sqrt{p})(1/\sqrt{n^{2}(n-1)}+1/\sqrt{n^{2}p})^{1/3}$
respectively, which are obtained by Theorem 1.1 in Johnstone (2001). This is
due to:

  1. i

Simulations behave better, especially for small $n,p$.

  2. ii

The limiting distribution does not change.

From this corollary, we know that for any data matrix ${\bm{B}}$ shifted by
the mean matrix $\bar{{\bm{B}}}$, the normalized largest eigenvalue of its
permuted matrices tends to follow the TW law. Note that for a data matrix
${\bm{B}}$, the auxiliary matrix ${\bm{D}}_{2}^{(n)}$ is known, and thus the
parameters $E_{+}$ and $\gamma_{0}$ can be directly obtained. For
illustration, we generate the data matrix ${\bm{B}}=\bm{\Psi}$, where
$\tilde{{\bm{T}}}={\bm{I}}_{p}$, $\bm{\Phi}$ is a matrix consisting of i.i.d.
standard Gaussian random variables. The dimensional settings are
$p=400,800,1600,3200,4800$ and the sample size is $n=3p/2$. Each ${\bm{B}}$ is
permuted $5,000$ times. The empirical quantiles of
$\varrho_{n}(\hat{{\bm{B}}}_{2})$ of the permuted matrices are provided in
Table 1. From these results, we can find that although there is distortion for
the quantiles when $p$ and $n$ are small, this distortion fades away as $p$
and $n$ increase. This indicates that in practice, if the data at hand are
centralized, the critical values obtained by PA can be estimated by the
quantiles of the TW distribution when $p$ and $n$ are large. This will
eliminate the massive permutation computations, while also getting relatively
accurate quantiles. We will discuss the asymptotic distribution of the largest
eigenvalue of permuted sample covariance matrices without centralization, as
well as other parallel analysis methods in Section 5.

Table 1: The percentiles of the normalized largest eigenvalues for permuted matrices Percentile | TW law | $p=400$ | $p=800$ | $p=1600$ | $p=3200$ | $p=4800$
---|---|---|---|---|---|---
$5\%$ | -3.1880 | -3.4281 | -3.3621 | -3.3511 | -3.3664 | -3.3108
$50\%$ | -1.2680 | -1.6086 | -1.5405 | -1.4960 | -1.4475 | -1.3835
$95\%$ | 0.9765 | 0.5968 | 0.6401 | 0.7339 | 0.8246 | 0.8795

###  3.2 Real data example

We consider the HGDP dataset (e.g., Li et al. (2008) and Dobriban and Owen
(2019)). The purpose of collecting this dataset is to evaluate the diversity
in the patterns of genetic variation across the globe. We use the Centre
d’Etude du Polymorphisme Humain panel, in which SNP data were collected for
1043 samples representing 51 populations from Africa, Europe, Asia, Oceania
and America. The data can be obtained from
http://www.hagsc.org/hgdp/data/hgdp.zip.

We use the same processing pipeline as Dobriban and Owen (2019). The setting
can be found at github.com/dobriban/DPA. The dataset has $n=1043$
observations, and $p=9730$ SNPs on chromosome 22, which are standardized. Thus
we have a $p\times n$ data matrix ${\bm{B}}$. Figure 1 shows the numerical
behavior of the limiting distribution of $\varrho_{n}(\hat{{\bm{B}}}_{2})$
based on 5,000 independent permutations. We find that the numerical
distribution fits the theoretical distribution (i.e., the TW distribution)
comparatively well.

Figure 1: Numerical behavior of the permuted sample covariance matrix on the
HGDP data.

##  4 Technical proofs

###  4.1 Proof strategy

The proof of Theorem 2.3 is completed with the aid of a general strategy
presented in Knowles and Yin (2017) for the covariance type matrices, which
itself is an adaptation of the method originally described in Erdös et al.
(2012) for Wigner matrices. Generally, to prove the TW law for the largest
eigenvalue, first, one needs to prove a local law for the spectral
distribution, which controls the location of the eigenvalues on an optimal
local scale. Second, with the aid of the local law, one needs to perform Green
function comparison between the matrix of interest and a certain reference
matrix ensemble, whose edge spectral behavior is already known. In Erdös et
al. (2012), an extended criterion of the local law for the covariance type of
matrices with independent columns (or rows) was given, see Theorem 3.6 in
Pillai and Yin (2014). It allows one to relax the independence assumption on
the entries within columns (or rows) to a certain extent, as long as some
large deviation estimates hold for certain linear and quadratic forms of each
column (or row) of the data matrix; see Lemma 3.4 in Pillai and Yin (2014).

Recall the sample covariance matrix ${{\bm{S}}}_{n}$ defined in (1.1). Our
main task is to show a large deviation estimate for each row of
${\bm{X}}_{n}$. Once the large deviation estimate is established, we can use
the criterion in Proposition 5.1 in Wen et al. (2020) to conclude the weak
local law of ${{\bm{S}}}_{n}$ and Lemma 5.6 in Knowles and Yin (2017) for
strong local law. In fact, Green function comparison can be completed
similarly to the approach in Erdös et al. (2012) by choosing an appropriate
reference matrix ensemble. The reference matrix that is chosen in this paper
is the traditional unbiased estimator of the population covariance matrix.

###  4.2 Notations

For two positive quantities $A_{n}$ and $B_{n}$ that depend on $n$, we use the
notation $A_{n}\asymp B_{n}$ to mean that $C^{-1}A_{n}\leq B_{n}\leq CA_{n}$
for some positive constant $C>1$. We need the following probability comparison
definition from Erdös et al. (2013).

Let $\mathcal{X}\equiv\mathcal{X}^{(N)}$ and
$\mathcal{Y}\equiv\mathcal{Y}^{(N)}$ be two sequences of nonnegative random
variables. We say that $\mathcal{Y}$ stochastically dominates $\mathcal{X}$
if, for all (small) $\epsilon>0$ and (large) $D>0$,

$\displaystyle\mathbb{P}(\mathcal{X}^{(N)}>N^{\epsilon}\mathcal{Y}^{(N)})\leq
N^{-D},$ (4.1)

for sufficiently large $N>N_{0}(\epsilon,D)$, and we write
$\mathcal{X}\prec\mathcal{Y}$ or $\mathcal{X}=O_{\prec}(\mathcal{Y})$. When
$\mathcal{X}^{(N)}$ and $\mathcal{Y}^{(N)}$ depend on a parameter
$v\in\mathcal{V}$ (typically an index label or a spectral parameter), then
$\mathcal{X}\prec\mathcal{Y}$ uniformly in $v\in\mathcal{V}$, which means that
the threshold $N_{0}(\epsilon,D)$ can be chosen independently of $v$. We use
the symbols $O(\cdot)$ and $o(\cdot)$ for the standard big-O and little-o
notations, respectively. We use $c$ and $C$ to denote strictly positive
constants that do not depend on $N$. Their values may change from line to
line. For any matrix ${\bm{A}}$, we denote $\|{\bm{A}}\|$ as its operator
norm, while for any vector $\bm{v}$, we use $\|\bm{v}\|$ to denote its
$L_{2}$-norm. In addition, we use ${\bm{a}}_{i}$ to represent the $i$-th row
of matrix ${\bm{A}}$.

###  4.3 Global law of $S_{n}$

In this section, we prove Theorem 2.2 (Global law of $S^{(n)}$).

###### Proof of Theorem 2.2.

When ${\bm{T}}$ is identity, Bai and Zhou (2008) proved that $F_{n}$ is
asymptotically given by the standard MP law. For general ${\bm{T}}$, we have

$\displaystyle\mathbb{E}w_{ij}^{2}$
$\displaystyle=nt_{i}\mathbb{E}x_{ij}^{2}=t_{i},~{}~{}\mathbb{E}w_{i1}w_{i2}=nt_{i}\mathbb{E}x_{i1}x_{i2}=-\frac{t_{i}}{N-1},$
$\displaystyle\mathbb{E}w_{i1}^{2}w_{i2}w_{i3}$
$\displaystyle=\frac{1}{N(N-1)(N-2)}\sum_{1\leq s\neq t\neq l\leq
N}u_{is}^{2}u_{it}u_{il}$
$\displaystyle=\frac{2}{(N-1)(N-2)}(\sum_{s=1}^{N}u_{is}^{4}/N)-\frac{N}{(N-1)(N-2)}t_{i}^{2},$
$\displaystyle\mathbb{E}w_{i1}w_{i2}w_{i3}w_{i4}$
$\displaystyle=\frac{1}{N(N-1)(N-2)(N-3)}\sum_{1\leq s\neq t\neq l\neq m\leq
N}u_{is}u_{it}u_{il}u_{im}$
$\displaystyle=\frac{1}{(N-1)(N-2)(N-3)}(-6(\sum_{s=1}^{N}u_{is}^{4}/N)+3Nt_{i}^{2}),$
$\displaystyle\mathbb{E}(w_{i1}^{2}-\mathbb{E}w_{i1}^{2})(w_{i2}^{2}-\mathbb{E}w_{i2}^{2})$
$\displaystyle=\frac{1}{N(N-1)}\sum_{1\leq s\neq t\leq
N}u_{is}^{2}u_{it}^{2}-t_{i}^{2}$
$\displaystyle=\frac{1}{N-1}(t_{i}^{2}-\sum_{s=1}^{N}u_{is}^{4}/N).$

Therefore, (2.7) follows from Corollary 1.1 in Bai and Zhou (2008), which
completes the proof of Theorem 2.2. ∎

###  4.4 Local law of $S_{n}$

In this section, our goal is to prove a strong local law for matrix
${{\bm{S}}}_{n}$. To achieve this goal, we first establish some large
deviation estimates for certain linear and quadratic forms of ${\bm{x}}_{i}$,
which are rows of the matrix ${\bm{X}}_{n}$. We now introduce a lemma on the
stochastic domination of high-order moments.

######  Lemma 4.1.

(Lemma 7.1 in Benaych-Georges and Knowles (2016)) Suppose that the
deterministic control parameter $\Phi$ satisfies $\Phi\geq n^{-C}$ for some
constant $C>0$ and that for all $p$, there is a constant $C_{p}$ such that the
nonnegative random variable $X$ satisfies $\mathbb{E}X^{p}\leq n^{C_{p}}$.
Then, we have the following equivalence

$X\prec\Phi\Longleftrightarrow\mathbb{E}X^{n}\prec\Phi^{n}$

for any fixed $n\in\mathbb{N}.$

For simplicity, we set $\frac{1}{\sqrt{n}}{\bm{v}}_{i}=(v_{1},\ldots,v_{n})$
and ${\bm{x}}_{i}=(x_{1},\ldots,x_{n})$ as the $i$-th row of
$\frac{1}{\sqrt{n}}{\bm{V}}_{p\times n}$ and ${\bm{X}}_{n}$, respectively. We
emphasize that in this subsection, $N=n$.

####  4.4.1 Large deviation estimates

######  Lemma 4.2.

(Large deviation lemma). For any deterministic vector
${\bm{A}}={A_{i}}\in\mathbb{C}_{n}$ and deterministic matrix
${\bm{B}}={B_{ij}}\in\mathbb{C}_{n\times n}$, we have

$\displaystyle|\sum_{i=1}^{n}x_{i}A_{i}|\prec(\|{\bm{A}}\|^{2}/n)^{1/2},$
(4.2) $\displaystyle|\sum_{i=1}^{n}B_{ii}x_{i}^{2}-n^{-1}{\rm
Tr}{\bm{B}}|\prec n^{-1}(\sum_{i=1}^{n}|B_{ii}|^{2})^{1/2},$ (4.3)
$\displaystyle|\sum_{i\neq j}^{n}x_{i}B_{ij}x_{j}+\frac{1}{n(n-1)}\sum_{i\neq
j}B_{ij}|\prec n^{-1}(\sum_{i\neq j}^{n}|B_{ij}|^{2})^{1/2}.$ (4.4)

###### Proof of Lemma 4.2.

The proof of Lemma 4.2 is very similar to the proof of Proposition 2.1 in Bao
(2017).

By Chebyshev’s inequality

$\mathbb{P}(|x_{i}|\geq
n^{\epsilon}n^{-\frac{1}{2}})\leq\mathbb{E}\frac{|\sqrt{n}x_{i}|^{q}}{n^{q\epsilon}},$

and (2.10), we have $|x_{i}|=O_{\prec}(n^{-1/2})$.

For (4.2), we first construct a sequence of martingale differences. Define the
filtration

$\mathcal{F}_{0}=\varnothing,~{}~{}~{}\mathcal{F}_{l}:=\sigma\\{x_{1},\ldots,x_{l}\\},~{}~{}~{}l\in[1,n],$

and set

$\displaystyle\mathcal{M}_{l}=\sum_{j=1}^{n}A_{j}[\mathbb{E}(x_{j}|\mathcal{F}_{l})-\mathbb{E}(x_{j}|\mathcal{F}_{l-1})].$
(4.5)

Conditioned on $\mathcal{F}_{l}$, $x_{i}$ is uniformly distributed for all
$x_{i}\in\\{x_{l+1},\ldots,x_{n}\\}$ by sampling without replacement. Hence,
we derive

$\displaystyle\mathbb{E}(x_{j}|\mathcal{F}_{l})=(\sum_{i=1}^{n}v_{i}-\sum_{i=1}^{l}x_{i})/(n-l)=-\sum_{i=1}^{l}x_{i}/(n-l),~{}~{}~{}j\geq
l+1,$ (4.6)

where we used the fact that $\sum_{j=1}^{n}x_{j}=\sum_{j=1}^{n}v_{j}=0$.
Furthermore, because $\mathbb{E}(x_{j}|\mathcal{F}_{l})=x_{j},~{}j\in[1,l]$
and $\mathbb{E}(x_{j}|\mathcal{F}_{l-1})=x_{j},~{}j\in[1,l-1]$,
$\mathcal{M}_{l}$ can be rewritten as follows:

$\displaystyle\mathcal{M}_{l}=x_{l}A_{l}+\sum_{j=l+1}^{n}A_{j}\mathbb{E}(x_{j}|\mathcal{F}_{l})-\sum_{j=l}^{n}A_{j}\mathbb{E}(x_{j}|\mathcal{F}_{l-1})$
$\displaystyle=$ $\displaystyle
A_{l}x_{l}-x_{l}\frac{\sum_{j=l+1}^{n}A_{j}}{n-l}+A_{l}\frac{\sum_{k=1}^{l-1}x_{k}}{n-l+1}+\sum_{j=l+1}^{n}A_{j}\frac{-\sum_{k=1}^{l-1}x_{k}}{n-l}-\sum_{j=l+1}^{n}A_{j}\frac{-\sum_{k=1}^{l-1}x_{k}}{n-l+1}$
$\displaystyle=$
$\displaystyle(A_{l}-\frac{1}{n-l}\sum_{j=l+1}^{n}A_{j})(x_{l}+\frac{1}{n-l+1}\sum_{k=1}^{l-1}x_{k}).$
(4.7)

Using the bound $|x_{k}|=O_{\prec}(n^{-1/2})$ and the fact that

$\displaystyle|\sum_{k=1}^{l-1}x_{k}|\prec n^{-1/2}\min(l-1,n-l+1),$ (4.8)

we deduce that

$\displaystyle|\mathcal{M}_{l}|\prec
n^{-1/2}(|A_{l}|+\frac{1}{n-l}\sum_{j=l+1}^{n}|A_{j}|).$ (4.9)

By the Burkholder inequality, we derive that for any fixed positive integer
$q>1$,

$\displaystyle\mathbb{E}|\sum_{l=1}^{n}\mathcal{M}_{l}|^{q}\leq
K_{q}\mathbb{E}(\sum_{l=1}^{n}\mathcal{M}_{l}^{2})^{q/2},$ (4.10)

where $K_{q}$ is finite depending on $q$. Now, we estimate
$\sum_{l=1}^{n}\mathcal{M}_{l}^{2}$. According to (4.9), we have

$\displaystyle\sum_{l=1}^{n}\mathcal{M}_{l}^{2}=$
$\displaystyle\sum_{l=1}^{n-1}\mathcal{M}_{l}^{2}\prec
n^{-1}\sum_{l=1}^{n-1}[\frac{1}{(n-l)^{2}}(\sum_{j=l+1}^{n}|A_{j}|)^{2}+|A_{l}|^{2}]$
$\displaystyle\prec$ $\displaystyle
n^{-1}\sum_{l=1}^{n-1}[\frac{1}{(n-l)^{2}}(\sum_{j=l+1}^{n}|A_{j}|)^{2}+|A_{l}|^{2}]\prec
n^{-1}\sum_{l=1}^{n-1}(\frac{1}{n-l}\sum_{j=l+1}^{n}|A_{j}|^{2}+|A_{l}|^{2})$
$\displaystyle\prec$ $\displaystyle n^{-1}\log n\|A\|^{2}\prec
n^{-1}\|A\|^{2}.$ (4.11)

Plugging (4.4.1) into (4.10), using Lemma 4.1 and Markov’s inequality, we can
find (4.2).

For (4.3), we construct a similar sequence of martingale differences as

$\mathcal{N}_{l}=\sum_{j=1}^{n}B_{jj}[\mathbb{E}(x_{j}^{2}|F_{l})-\mathbb{E}(x_{j}^{2}|F_{l-1})].$
(4.12)

$\mathcal{N}_{n}=0$ since $\mathcal{F}_{n}=\mathcal{F}_{n-1}$. Conditioned on
$\mathcal{F}_{l}$, $x_{i}$ is uniformly distributed for all
$x_{i}\in\\{x_{l+1},\ldots,x_{n}\\}$. Hence,

$\displaystyle\mathbb{E}(x_{j}^{2}|\mathcal{F}_{l})=\frac{\sum_{i_{=}1}^{n}v_{i}^{2}-\sum_{i=1}^{l}x_{i}^{2}}{n-l}=\frac{1-\sum_{i=1}^{l}x_{i}^{2}}{n-l},~{}~{}~{}j\geq
l+1.$ (4.13)

Applying (4.13), we deduce that

$\displaystyle\mathcal{N}_{l}=$ $\displaystyle
B_{ll}x_{l}^{2}+\sum_{j=l+1}^{n}B_{jj}\mathbb{E}(x_{j}^{2}|\mathcal{F}_{l})-\sum_{j=l}^{n}B_{jj}\mathbb{E}(x_{j}^{2}|\mathcal{F}_{l-1})$
$\displaystyle=$
$\displaystyle(B_{ll}-\frac{1}{n-l}\sum_{j=l+1}^{n}B_{jj})(x_{l}^{2}+\frac{1}{n-l+1}(1-\sum_{k=1}^{l-1}x_{k}^{2})).$
(4.14)

Based on the fact that $|x_{k}|=O_{\prec}(n^{-1/2})$, we have
$|\frac{1}{n-l}(1-\sum_{k=1}^{l-1}x_{k}^{2})|\prec n^{-1}.$ Thus,

$\displaystyle|\mathcal{N}_{l}|\prec
n^{-1}(|B_{ll}|+\frac{1}{n-l}\sum_{j=l+1}^{n}|B_{jj}|).$ (4.15)

The remaining proof of (4.3) is the same as that for (4.2). Thus, we omit the
details.

Finally, we prove (4.4). Again, we first set

$\displaystyle\mathcal{R}_{l}=\sum_{i<j}B_{ij}[\mathbb{E}(x_{i}x_{j}|\mathcal{F}_{l})-\mathbb{E}(x_{i}x_{j}|\mathcal{F}_{l-1})].$
(4.16)

Considering the different cases $i\leq l-1,~{}j=l$, $i\leq l-1,~{}j\geq l+1$,
$i=l,~{}j\geq l+1$, $j>i\geq l+1$, and $i<j\leq l-1$, we rewrite
$\mathcal{R}_{l}=\mathcal{R}_{l1}+\mathcal{R}_{l2}+\mathcal{R}_{l3}+\mathcal{R}_{l4}+\mathcal{R}_{l5}$,
where

$\displaystyle\mathcal{R}_{l1}=\sum_{i=1}^{l-1}B_{il}x_{i}[x_{l}-\mathbb{E}(x_{l}|\mathcal{F}_{l-1})];~{}\mathcal{R}_{l2}=\sum_{i=1}^{l-1}\sum_{j=l+1}^{n}B_{ij}x_{i}[\mathbb{E}(x_{j}|\mathcal{F}_{l})-\mathbb{E}(x_{j}|\mathcal{F}_{l-1})];$
$\displaystyle\mathcal{R}_{l3}=\sum_{j=l+1}^{n}B_{lj}[\mathbb{E}(x_{j}|\mathcal{F}_{l})x_{l}-\mathbb{E}(x_{j}x_{l}|\mathcal{F}_{l-1})];$
$\displaystyle\mathcal{R}_{l4}=\sum_{i=l+1}^{n}\sum_{j=i+1}^{n}B_{ij}[\mathbb{E}(x_{i}x_{j}|\mathcal{F}_{l})-\mathbb{E}(x_{i}x_{j}|\mathcal{F}_{l-1})];$
$\displaystyle\mathcal{R}_{l5}=\sum_{i<j\leq
l-1}B_{ij}[\mathbb{E}(x_{i}x_{j}|\mathcal{F}_{l})-\mathbb{E}(x_{i}x_{j}|\mathcal{F}_{l-1})]=0.$

Applying (4.2) and setting
$A=(\underbrace{1,\ldots,1}_{l},\underbrace{0,\ldots,0}_{n-l})^{\top}$ and
$A=(\underbrace{0,\ldots,0}_{l},\underbrace{1,\ldots,1}_{n-l})^{\top}$, we
have

$\displaystyle|\sum_{j=1}^{l-1}x_{j}|=|\sum_{j=l}^{n}x_{j}|\prec
n^{-1/2}\min\\{\sqrt{l-1},\sqrt{n-l+1}\\}.$ (4.17)

Therefore, we have

$\displaystyle\mathbb{E}(x_{j}|\mathcal{F}_{l-1})\prec
1/\sqrt{n(n-l+1)},~{}~{}j\geq l.$ (4.18)

Moreover, for $i,j\geq l$ and $i\neq j$, similar to (4.6) and (4.13), we have

$\displaystyle\mathbb{E}(x_{i}x_{j}|\mathcal{F}_{l-1})$ $\displaystyle=$
$\displaystyle\frac{1}{(n-l+1)(n-l)}\\{\sum_{i\neq j}x_{i}x_{j}-\sum_{i,j\leq
l-1,i\neq j}x_{i}x_{j}-2\sum_{i\leq l-1,j\geq l}x_{i}x_{j}\\}$
$\displaystyle=$
$\displaystyle\frac{1}{(n-l+1)(n-l)}\\{(\sum_{i}x_{i})^{2}-\sum_{i}x_{i}^{2}-(\sum_{i\leq
l-1}x_{i})^{2}+\sum_{i\leq l-1}x_{i}^{2}$ $\displaystyle-2\sum_{i\leq
l-1}x_{i}(\sum_{j=1}^{n}x_{j}-\sum_{j\leq l-1}x_{j})\\}$ $\displaystyle=$
$\displaystyle\frac{1}{(n-l+1)(n-l)}\\{(\sum_{i\leq
l-1}x_{i})^{2}-(1-\sum_{i\leq l-1}x_{i}^{2})\\}.$ (4.19)

In addition, we find that

$\displaystyle|1-\sum_{i=1}^{l-1}x_{i}^{2}|\prec n^{-1}(n-l+1).$ (4.20)

Thus, by (4.17), (4.4.1) and (4.20), we derive the following for $l\in[1,n]$.

$\displaystyle\mathbb{E}(x_{i}x_{j}|\mathcal{F}_{l-1})\prec n^{-1}(n-l)^{-1}.$
(4.21)

By the Burkholder inequality, we obtain for any fixed positive integer $q>1$,

$\displaystyle\mathbb{E}|\sum_{l=1}^{n}\mathcal{R}_{l}|^{q}\leq
K_{q}\mathbb{E}(\sum_{l=1}^{n}\mathcal{R}_{l}^{2})^{q/2},$ (4.22)

and by the Minkowski inequality, we have

$\displaystyle(\mathbb{E}(\sum_{l=1}^{n}\mathcal{R}_{l}^{2})^{q/2})^{2/q}\leq\sum_{l=1}^{n}(\mathbb{E}|\mathcal{R}_{l}|^{q})^{2/q}.$
(4.23)

Hence, bounding $\mathbb{E}|\mathcal{R}_{l\alpha}|^{q}$ for $\alpha=1,2,3,4$
is sufficient.

For $\mathbb{E}|\mathcal{R}_{l1}|^{q}$,

$\displaystyle\mathbb{E}|\mathcal{R}_{l1}|^{q}=$
$\displaystyle\mathbb{E}(|\sum_{i=1}^{l-1}B_{il}x_{i}|^{q}|x_{l}-\mathbb{E}(x_{l}|\mathcal{F}_{l-1})|^{q})$
$\displaystyle\prec$ $\displaystyle
n^{-q/2}\mathbb{E}|\sum_{i=1}^{l-1}B_{il}x_{i}|^{q}\prec
n^{-q}(\sum_{i=1}^{l-1}|B_{il}|^{2})^{q/2},$ (4.24)

which is due to (4.2) and $x_{i}=O_{\prec}(n^{-1/2})$.

For $\mathbb{E}|\mathcal{R}_{l2}|^{q}$, from (4.6) and (4.8), when $l\leq n$,
we have

$\displaystyle|\mathbb{E}(x_{j}|\mathcal{F}_{l})-\mathbb{E}(x_{j}|\mathcal{F}_{l-1})|=\frac{1}{n-l}|x_{l}+\frac{\sum_{j=1}^{l-1}x_{j}}{n-l+1}|\prec\frac{1}{(n-l)\sqrt{n}}.$
(4.25)

Thus, we have

$\displaystyle|\mathcal{R}_{l2}|\prec$
$\displaystyle\frac{1}{(n-l)\sqrt{n}}\sum_{j=l+1}^{n}|\sum_{i=1}^{l-1}B_{ij}x_{i}|\prec\frac{1}{(n-l)n}\sum_{j=l+1}^{n}(\sum_{i=1}^{l-1}|B_{ij}|^{2})^{1/2}$
$\displaystyle\prec$
$\displaystyle\frac{\sqrt{n-l}}{(n-l)n}(\sum_{j=l+1}^{n}\sum_{i=1}^{l-1}|B_{ij}|^{2})^{1/2}\prec\frac{1}{n\sqrt{n-l}}(\sum_{j=l+1}^{n}\sum_{i=1}^{l-1}|B_{ij}|^{2})^{1/2}.$

Hence, we have

$\displaystyle\mathbb{E}|\mathcal{R}_{l2}|^{q}\prec
n^{-q}(n-l)^{-q/2}(\sum_{j=l+1}^{n}\sum_{i=1}^{l-1}|B_{ij}|^{2})^{q/2}.$
(4.26)

Next, we estimate $\mathbb{E}|\mathcal{R}_{l3}|^{q}$. Obviously, according to
(4.17) and (4.21),

$\displaystyle|\mathcal{R}_{l3}|\prec\sum_{j=l+1}^{n}|B_{lj}|/(n\sqrt{n-l})\prec(\sum_{j=l+1}^{n}|B_{lj}|^{2})^{1/2}/n,$

In other words,

$\displaystyle\mathbb{E}|\mathcal{R}_{l3}|^{q}\prec(\sum_{j=l+1}^{n}|B_{lj}|^{2})^{q/2}/n^{q}.$
(4.27)

Now, we estimate $\mathbb{E}|\mathcal{R}_{l4}|^{q}$. From (4.17), (4.4.1) and
(4.20), we derive that for $i,j>l$,

$\displaystyle|\mathbb{E}(x_{i}x_{j}|\mathcal{F}_{l})-\mathbb{E}(x_{i}x_{j}|\mathcal{F}_{l-1})|$
$\displaystyle=$ $\displaystyle|\frac{1}{(n-l-1)(n-l)(n-l+1)}(\\{(\sum_{i\leq
l-1}x_{i})^{2}-(1-\sum_{i\leq l-1}x_{i}^{2})\\})$
$\displaystyle+\frac{2}{(n-l)(n-l-1)}x_{l}\sum_{j=1}^{l}x_{j}|\prec\frac{1}{n(n-l)^{3/2}}.$

Hence, we obtain

$\displaystyle\mathcal{R}_{l4}\prec\frac{1}{n(n-l)^{3/2}}\sum_{i=l+1}^{n}\sum_{j=i+1}^{n}|B_{ij}|\prec\frac{1}{n(n-l)^{1/2}}(\sum_{i=l+1}^{n}\sum_{j=i+1}^{n}|B_{ij}|^{2})^{1/2}.$
(4.28)

Therefore,

$\displaystyle\mathbb{E}|\mathcal{R}_{l4}|^{q}\prec\frac{1}{n^{q}(n-l)^{q/2}}(\sum_{i=l+1}^{n}\sum_{j=i+1}^{n}|B_{ij}|^{2})^{q/2}.$
(4.29)

Finally, from (4.4.1), (4.26), (4.27) and (4.29), we conclude that

$\displaystyle(\mathbb{E}|\mathcal{R}_{l}|^{q})^{2/q}=$
$\displaystyle(\mathbb{E}(|\sum_{\alpha=1}^{4}\mathcal{R}_{l\alpha}|^{2})^{q/2})^{2/q}\leq
4(\mathbb{E}(\sum_{\alpha=1}^{4}\mathcal{R}_{l\alpha}^{2})^{q/2})^{2/q}\leq
4\sum_{\alpha=1}^{4}(\mathbb{E}|\mathcal{R}_{l\alpha}^{q}|)^{2/q}$
$\displaystyle\prec$ $\displaystyle
n^{-2}(\sum_{i=1}^{l-1}|B_{il}|^{2}+\sum_{j=l+1}^{n}|B_{lj}|^{2}$
$\displaystyle+\frac{1}{n-l}\sum_{j=l+1}^{n}\sum_{i=1}^{l-1}|B_{ij}|^{2}+\frac{1}{n-l}\sum_{i=l+1}^{n}\sum_{j=i+1}^{n}|B_{ij}|^{2})$
$\displaystyle\prec$ $\displaystyle
n^{-2}(\sum_{i=1}^{l-1}|B_{il}|^{2}+\sum_{j=l+1}^{n}|B_{lj}|^{2}+\frac{1}{n-l}\sum_{j=l+1}^{n}\sum_{i=1}^{j-1}|B_{ij}|^{2}).$

Since

$\displaystyle\sum_{l=1}^{n}\sum_{j=l+1}^{n}\sum_{i=1}^{j-1}\frac{1}{n-l}|B_{ij}|^{2}=\sum_{j=1}^{n}\sum_{i=1}^{j-1}(\sum_{l=1}^{j-1}\frac{1}{n-l})|B_{ij}|^{2}\leq
C\log
n(\sum_{j=1}^{n}\sum_{i=1}^{j-1})|B_{ij}|^{2}\prec\sum_{i<j}|B_{ij}|^{2},$

we conclude that

$\displaystyle\mathbb{E}|\sum_{l=1}^{n}\mathcal{R}_{l}|^{q}\prec
n^{-q}(\sum_{i<j}|B_{ij}|^{2})^{q/2}.$

Hence, by Markov’s inequality, we have

$\displaystyle|\sum_{i\neq j}x_{i}B_{ij}x_{j}-\mathbb{E}(\sum_{i\neq
j}x_{i}B_{ij}x_{j})|\prec n^{-1}(\sum_{i\neq j}|B_{ij}|^{2})^{1/2}.$ (4.30)

By the fact that $\mathbb{E}x_{i}x_{j}=-\frac{1}{n(n-1)},~{}i\neq j$, Lemma
4.2 is concluded. ∎

Furthermore, for each $i=1,\ldots,p$, we let
$\bm{\xi}_{i}=\\{\xi_{i1},\ldots,\xi_{in}\\}$ be a random vector with i.i.d.
components that are uniformly distributed on ${\bm{v}}_{i}$, i.e.,
$\mathbb{P}(\xi_{ik}=v_{ij})=n^{-1}$, $k=1,\ldots,n,~{}j=1,\ldots,n$. Note
that $x_{ij}$’s are also identically distributed as $\xi_{ij}$’s, but
$x_{ij}$’s are correlated. Set

$\displaystyle\tilde{{\bm{x}}}_{i}=\bm{\xi}_{i}\bm{\Sigma}^{\frac{1}{2}},$
(4.31)

where

$\displaystyle\bm{\Sigma}=\frac{n}{n-1}\bm{I}_{n}-\frac{1}{n-1}{\bm{e}}_{n}{\bm{e}}_{n}^{\top}.$
(4.32)

Let $\bm{\Xi}$ and $\tilde{{\bm{X}}}_{n}$ be the $p\times n$ matrices with
$\bm{\xi}_{i}$ and $\tilde{{\bm{x}}}_{i}$ as their $i$-th rows, respectively.
Introduce

$\displaystyle\tilde{{\bm{S}}}={\bm{T}}^{\frac{1}{2}}\tilde{{\bm{X}}}_{n}\tilde{{\bm{X}}}_{n}^{\top}{\bm{T}}^{\frac{1}{2}}={\bm{T}}^{\frac{1}{2}}\bm{\Xi}\bm{\Sigma}\bm{\Xi}^{\top}{\bm{T}}^{\frac{1}{2}}.$
(4.33)

Similarly, set $\tilde{{\bm{y}}}_{i}=\sqrt{t_{i}}\tilde{{\bm{x}}}_{i}$,
$\tilde{{\bm{Y}}}_{n}=(\tilde{{\bm{y}}}_{1}^{\top},\ldots,\tilde{{\bm{y}}}_{p}^{\top})^{\top}$
and rewrite (4.33) as

$\displaystyle\tilde{{\bm{S}}}=\tilde{{\bm{Y}}}_{n}\tilde{{\bm{Y}}}_{n}^{\top}.$
(4.34)

Observe that $\tilde{{\bm{S}}}$ is the classical sample covariance matrix in
statistics theory, which is centralized by sample mean and scaled by $n-1$,
although in random matrix theory the simplified model
${\bm{T}}^{\frac{1}{2}}\bm{\Xi}\bm{\Xi}^{\top}{\bm{T}}^{\frac{1}{2}}$ is
considered more often. Since
${\bm{T}}^{\frac{1}{2}}\bm{\Xi}\bm{\Sigma}\bm{\Xi}^{\top}{\bm{T}}^{\frac{1}{2}}$
is simply a rank one perturbation of
$\frac{n}{n-1}{\bm{T}}^{\frac{1}{2}}\bm{\Xi}\bm{\Xi}^{\top}{\bm{T}}^{\frac{1}{2}}$,
the empirical spectral distribution of $\tilde{{\bm{S}}}$ almost surely also
weakly converges to a probability distribution, whose Stieltjes transform is
given by (2.7).

Observe that the entries of $\bm{\Xi}$ are i.i.d.. Then, using the large
deviation estimates for linear and quadratic forms of the i.i.d. random
variables (e.g., Corollary B.3 in Erdös et al. (2012)), we see that (4.2),
(4.3) and (4.4) still hold if we replace ${\bm{x}}_{i}$ with
$\tilde{{\bm{x}}}_{i}$.

####  4.4.2 Strong local law for $S_{n}$

We denote $\gamma_{1}\geq\cdots\geq\gamma_{p\wedge n}$ as the ordered
$p$-quantiles of $F_{c_{n},H_{n}}$, i.e., $\gamma_{j}$ is the smallest real
number such that

$\displaystyle\int_{-\infty}^{\gamma_{j}}{\rm
d}F_{c_{n},H_{n}}(x)=\frac{p-j+1}{p},~{}j=1,\ldots p\wedge n.$ (4.35)

Recall matrix ${{\bm{S}}}_{n}$ defined in (2.1) and matrix $\tilde{{\bm{S}}}$
defined in (4.34). We introduce some intermediate matrices between
${{\bm{S}}}_{n}$ and $\tilde{{\bm{S}}}$. Starting from ${\bm{Y}}_{n}$, we
replace ${\bm{y}}_{i}$’s with $\tilde{{\bm{y}}}_{i}$’s one by one and obtain a
sequence of intermediate matrices

$\displaystyle{\bm{Y}}_{n}({\bm{Y}}_{0}),~{}~{}{\bm{Y}}_{1},\ldots,{\bm{Y}}_{l},{\bm{Y}}_{l+1},\ldots,{\bm{Y}}_{p-1},~{}~{}\tilde{{\bm{Y}}}_{n}({\bm{Y}}_{p}).$
(4.36)

Correspondingly, we set

$\displaystyle{\bm{S}}_{l}={\bm{Y}}_{l}{\bm{Y}}_{l}^{\top},~{}~{}$
$\displaystyle\bm{G}_{l}(z)=({\bm{S}}_{l}-z)^{-1},~{}~{}m_{l}(z)=\frac{1}{p}{\rm{\rm
Tr}}\bm{G}_{l}(z)$ (4.37)
$\displaystyle\mathcal{{\bm{S}}}_{l}={\bm{Y}}_{l}^{\top}{\bm{Y}}_{l},~{}~{}$
$\displaystyle\mathcal{G}_{l}(z)=(\mathcal{{\bm{S}}}_{l}-z)^{-1},~{}~{}\underline{m}_{l}(z)=\frac{1}{n}{\rm{\rm
Tr}}\mathcal{G}_{l}(z),$ (4.38)

When $l=0$, we have ${\bm{Y}}_{0}={\bm{Y}}_{n}$, $m_{0}=m_{n}$; then, we
simply write ${\bm{S}}_{0}$; $\mathcal{{\bm{S}}}_{0}$; $\bm{G}_{0}$;
$\mathcal{G}_{0}$; $\underline{m}_{0}$ as ${\bm{S}}$; $\mathcal{{\bm{S}}}$
$\bm{G}$; $\mathcal{G}$; and $\underline{m}_{n}$.

Furthermore, we introduce the following notation:

$\displaystyle\Lambda_{d}=\max_{k,l}|\bm{G}_{kl}-\delta_{kl}\tilde{m}_{k}|,~{}~{}\Lambda_{o}=\max_{k\neq
l}|\bm{G}_{kl}|,~{}~{}\Lambda=|m_{n}-\tilde{m}|=c_{n}^{-1}|\underline{m}_{n}-\tilde{\underline{m}}|,$
(4.39)

where $\delta_{kl}$ denotes the Kronecker delta, i.e. $\delta_{kl}=1$ if
$k=l$, and $\delta_{kl}=0$ if $k\neq l$, and
$\tilde{m}_{k}=\frac{1}{t_{k}(1-c_{n}-c_{n}z\tilde{m})-z}=-\frac{1}{z+z\tilde{\underline{m}}t_{k}}$.
The last equality is derived by (2.17).

For a small positive $\bm{c},\epsilon$ and sufficiently large $C_{+}$,
$(C_{+}>E_{+})$, we define the domain

$\displaystyle
D(\bm{c},\epsilon):=\\{z=E+i\eta\in\mathbb{C}^{+}:E_{+}-\bm{c}\leq E\leq
C_{+},~{}n^{-1+\epsilon}\leq\eta\leq 1\\}.$ (4.40)

The following lemma on $\tilde{\underline{m}}(z)$ is elementary.

######  Lemma 4.3.

(Theorem 3.1 in Bao et al. (2015)) Under Assumption 2.1, there is a $c>0$ such
that for any $z=E+i\eta\in D(c,0)$

$\displaystyle|\tilde{\underline{m}}(z)|\asymp 1,$ (4.41)
$\displaystyle{\rm{\rm
Im}}\tilde{\underline{m}}(z)\asymp\left\\{\begin{array}[]{ll}\sqrt{\kappa+\eta},&~{}~{}~{}~{}if~{}E\in[E_{+}-c,E_{+}+\eta),\\\
\frac{\eta}{\sqrt{\kappa+\eta}},&~{}~{}~{}~{}if~{}E\geq
E_{+}+\eta,\end{array}\right.$ (4.44)
$\displaystyle|1+t\tilde{\underline{m}}(z)|>0,~{}~{}~{}t\in[t_{p},t_{1}],$
(4.45)

where $\kappa\equiv\kappa(E)=|E-E_{+}|$.

We further define the control parameter
$\Omega\equiv\Omega(z)=\sqrt{\frac{{\rm{\rm
Im}}\tilde{m}(z)}{n\eta}}+\frac{1}{n\eta}.$ We claim that the following local
law holds.

######  Lemma 4.4.

Under Assumption 2.1, we have that for any sufficiently small $\epsilon>0$,

  1. 1.

(Entrywise local law):

$\Lambda_{d}(z)\prec\Omega(z)$

holds uniformly on $D(c,\epsilon)$.

  2. 2.

(Average local law):

$\Lambda(z)\prec\frac{1}{n\eta}$

holds uniformly on $D(c,\epsilon)$, where $c$ is the constant in Lemma 4.3.

  3. 3.

(Rigidity on the right edge): For $i\in[1,\delta p]$ with any sufficiently
small constant $\delta\in(0,1)$, we have

$|\lambda_{i}({\bm{S}}_{n})-\gamma_{i}|\prec
n^{-\frac{2}{3}}i^{-\frac{1}{3}}.$

All the above hold if we replace ${{\bm{S}}}_{n}$ with ${\bm{S}}_{l}$ for all
$l=1,\ldots,p$.

The proof of Lemma 4.4 is deferred to Section A.1 of the Supplementary
Material (Hu et al. (2022)).

###  4.5 Edge universality for $S_{n}$

In this section, we prove the edge universality of ${{\bm{S}}}_{n}$ by the
Green function comparison.

Recall the intermediate matrices defined in (4.36) and the notation introduced
in (4.37). We aim to show the following lemma.

######  Lemma 4.5.

Fix any $\gamma\in[1,\ldots,p]$, and choose $\epsilon$ as any sufficiently
small constant. Let $E,E_{1},E_{2}\in R$ satisfy $E_{1}<E_{2}$ and

$\displaystyle|E-E_{+}|,~{}|E_{1}-E_{+}|,~{}|E_{2}-E_{+}|\leq
n^{-\frac{2}{3}+\epsilon},$ (4.46)

and set $\eta_{0}=n^{-\frac{2}{3}-\epsilon}$. Define
$F:\mathbb{R}\rightarrow\mathbb{R}$ as a smooth function satisfying

$\displaystyle\max_{x\in\mathbb{R}}|F^{(l)}(x)|(|x|+1)^{-C}\leq
C,~{}~{}~{}l=1,2,3,4$

for some positive constant $C$. Then, there exists a constant $\delta>0$ such
that for sufficiently large $n$, we have

$\displaystyle|\mathbb{E}F(n\eta_{0}{\rm
Im}m_{\gamma}(z))-\mathbb{E}F(n\eta_{0}{\rm Im}m_{\gamma+1}(z))|\prec
n^{-1-\delta},~{}~{}~{}~{}~{}~{}~{}z=E+i\eta_{0},$ (4.47)

and

$\displaystyle|\mathbb{E}F(n\int_{E_{1}}^{E_{2}}{\rm
Im}m_{\gamma}(x+i\eta_{0}){\rm d}x)-\mathbb{E}F(n\int_{E_{1}}^{E_{2}}{\rm
Im}m_{\gamma+1}(x+i\eta_{0}){\rm d}x)|\prec n^{-1-\delta}.$ (4.48)

The proof of Lemma 4.5 can be found in Section A.2 of the Supplementary
Material (Hu et al. (2022)). Using Lemma 4.5, we can now prove Theorem 2.3.

###### Proof of Theorem 2.3.

Since the proof of Theorem 2.3 is similar to that of Theorem 3.3 in Wen et al.
(2020), we only outline the general strategy here. First, we attempt to
convert the problem to a comparison of the Stieltjes transform. As a result,
for any $\eta>0$, we define

$\vartheta_{\eta}(x)=\frac{\eta}{\pi(x^{2}+\eta^{2})}=\frac{1}{\pi}\operatorname{Im}\frac{1}{x-i\eta}.$

We notice that for any $a,b\in\mathbb{R}$ with $a\leq b$, the convolution of
$\mathbf{1}_{[a,b]}$ and $\vartheta_{\eta}$ applied to the eigenvalues
$\lambda_{i}$, $i=1,\dots,p$ yields that

$\sum_{i=1}^{p}\mathbf{1}_{[a,b]}*\vartheta_{\eta}(\lambda_{i})=\frac{p}{\pi}\int_{a}^{b}\operatorname{Im}m_{n}(x+i\eta){\rm
d}x.$

In terms of functional calculus notation, we have

$\displaystyle\sum_{i=1}^{p}\mathbf{1}_{[a,b]}(\lambda_{i})={\rm
Tr}\mathbf{1}_{[a,b]}({{\bm{S}}}_{n}),\quad\sum_{i=1}^{p}\mathbf{1}_{[a,b]}*\vartheta_{\eta}(\lambda_{i})={\rm
Tr}\mathbf{1}_{[a,b]}*\vartheta_{\eta}({{\bm{S}}}_{n}).$

For $a,b\in\mathbb{R}\cup\\{-\infty,\infty\\}$, define ${\cal
N}(a,b)=\\#(a\leq\lambda_{j}\leq b)$ as the number of eigenvalues of
${{\bm{S}}}_{n}$ in $[a,b]$. Then, we have the following lemma.

######  Lemma 4.6.

(Lemma 4.1 in Pillai and Yin (2014)) Let $\varepsilon>0$ be an arbitrarily
small number. Set $E_{\varepsilon}=E_{+}+n^{-2/3+\varepsilon}$,
$\ell_{1}=n^{-2/3-3\varepsilon}$ and $\eta_{1}=n^{-2/3-9\varepsilon}$. Then,
for any $E$ satisfying $|E-E_{+}|\leq\frac{3}{2}n^{-2/3+\varepsilon}$, it
holds with high probability that

$|{\rm Tr}\mathbf{1}_{[E,E_{\varepsilon}]}({{\bm{S}}}_{n})-{\rm
Tr}\mathbf{1}_{[E,E_{\varepsilon}]}*\vartheta_{\eta}({{\bm{S}}}_{n})|\leq
C(n^{-2\varepsilon}+{\cal N}(E-\ell_{1},E+\ell_{1})).$

Lemma 4.6 shows that ${\rm Tr}\mathbf{1}_{[a,b]}({{\bm{S}}}_{n})$ can be well
approximated by its smoothed version ${\rm
Tr}\mathbf{1}_{[a,b]}*\vartheta_{\eta}({{\bm{S}}}_{n})$ for $a,b$ around edge
$E_{+}$ so that the problem can be converted to a comparison of the Stieltjes
transform.

Then, from Lemma 4.5, Lemma 4.6, the square root behavior of
$\underline{\rho}_{0}$ (see Section 2.2), and the rigidity on the right edge
(see Lemma 4.4), one can see that

$\displaystyle\mathbb{P}\left(n^{2/3}(\lambda_{1}(\tilde{{\bm{S}}})-E_{+})\leq
s-n^{-\varepsilon}\right)-n^{-\delta}\leq\mathbb{P}\left(n^{2/3}(\lambda_{1}({{\bm{S}}}_{n})-E_{+})\leq
s\right)$ $\displaystyle\leq$
$\displaystyle\mathbb{P}\left(n^{2/3}(\lambda_{1}(\tilde{{\bm{S}}})-E_{+})\leq
s+n^{-\varepsilon}\right)+n^{-\delta},$ (4.49)

where $\tilde{{\bm{S}}}$ is defined in (4.33). It is known from the eigenvalue
sticking result (Theorem 2.7 of Bloemendal et al. (2016) and Section 11 in
Knowles and Yin (2017)) that the edge universality for $\tilde{{\bm{S}}}$
holds. Specifically,

$\displaystyle\mathbb{P}\left(n^{2/3}(\lambda_{1}({\bm{Q}})-E_{+})\leq
s-n^{-\varepsilon}\right)-n^{-\delta}\leq\mathbb{P}\left(n^{2/3}(\lambda_{1}(\tilde{{\bm{S}}})-E_{+})\leq
s\right)$ $\displaystyle\leq$
$\displaystyle\mathbb{P}\left(n^{2/3}(\lambda_{1}({\bm{Q}})-E_{+})\leq
s+n^{-\varepsilon}\right)+n^{-\delta},$ (4.50)

where ${\bm{Q}}$ is defined in Theorem 2.3. Combining (4.5) and (4.5), we can
conclude the proof of Theorem 2.3. ∎

##  5 Conclusion and discussion

In this paper, we focus on the eigenvalues of high-dimensional sample
covariance matrices generated without replacement from finite populations. The
largest eigenvalues are proven to follow the TW law. Based on the permutation
method introduced by Dobriban (2020), we consider an analogous strategy to
compare the top singular values of ${\bm{B}}$ and a certain percentile of the
top singular values of its permuted matrices after centralization, where the
latter can be attributed to our theoretical results.

Technically, for the PA, one might be more interested in the asymptotic
distribution of
$\lambda_{1}(\frac{1}{n}({\bm{B}})_{\pi}({\bm{B}})_{\pi}^{\top})$, that is the
largest eigenvalue of permuted sample covariance matrices without
centralization. We believe that it would still tend to the TW law as
$n,~{}p\to\infty$ by adding some more restrictions on ${\bm{S}}$. However, it
requires more technical work, thus we decide to leave it for future work.
Besides, Dobriban and Owen (2019) proposed derandomized PA (DPA), deflated DPA
(DDPA) and DDPA+ methods, which are faster and more reproducibly than the PA
method. Essentially, our theoretical results can explain why the DPA method is
a kind of parallel analysis from another perspective. For DDPA and DDPA+
methods, both counter the shadowing problem of DPA by removing some largest
factors one by one. We think that the idea of removing the poorly estimated
singular vectors can also improve the accuracy of the quantiles of the TW
distribution in practice. We will also leave the rigorous discussion of this
property for future work.

##  Acknowledgments

The authors would like to thank Zhigang Bao for allowing them to use the
method in his unpublished work Bao (2017). We also thank the editor, AE and
referee for their careful reading and invaluable comments, which have greatly
improved the quality of the paper.

## Funding

Jiang Hu’s research was supported by NSFC Nos. 12171078, 11971097, 12292980,
12292982. Shaochen Wang’s research was supported by NSFC No. 11801181,
Guangdong Basic and Applied Basic Research Foundation No. 2018A030310358.
Yangchun Zhang’s research was supported by Shanghai Sailing Program
21YF1413500. Wang Zhou’s research was partially supported by a grant
A-0004803-00-00 at the National University of Singapore.

## References

  * Anderson (1984) Anderson, T. W. (1984). An Introduction to Multivariate Statistical Analysis, 2nd ed. Wiley, New York.
  * Bai and Silverstein (2010) Bai, Z. D. and Silverstein, J. W. (2010) Spectral Analysis of Large Dimensional Random Matrices, 2nd edn. Springer.
  * Bai and Zhou (2008) Bai, Z. D. and Zhou, W. (2008) Large sample covariance matrices without independence structure in columns. Stat. Sin., 18, 425-442.
  * Bao (2017) Bao, Z. G. (2017). Tracy-Widom limit for Spearman’s rho. preprint.
  * Bao (2019) Bao, Z. G. (2019). Tracy-Widom limit for Kendall’s tau. Ann. Stat. 47(6) 3504–3532.
  * Bao et al. (2015) Bao, Z. G.,  Lin, L. C., Pan, G. M. and  Zhou, W. (2015). Spectral statistics of large dimensional Spearman’s rank correlation matrix and its application. Ann. Stat. 43(6) 2588–2623.
  * Bao et al. (2012) Bao, Z. G., Pan, G. M. and Zhow, W. (2012). Tracy-Widom law for the extreme eigenvalues of sample correlation matrices. Electron. J. Probab. 88(17) 1–32.
  * Bao et al. (2015) Bao, Z. G., Pan, G. M. and Zhow, W. (2015). Universality for the largest eigenvalue of sample covariance matrices with general population. Ann. Stat. 43(1) 382–421.
  * Bloemendal et al. (2016) Bloemendal, A., Knowles, A., Yau, H.-T and Yin, J. (2016). On the principal components of sample covariance matrices. Probab. Theory Relat. Fields 164(2) 459–552.
  * Buja and Eyuboglu (1992) Buja, A. and Eyuboglu, N. (1992). Remarks on parallel analysis. Multivariate behavioral research 27(4) 509–540.
  * Dobriban (2020) Dobriban, E. (2020). Permutation methods for factor analysis and PCA. Ann. Stat. 48(5) 2824–2847
  * Dobriban and Owen (2019) Dobriban, E. and Owen. A. B. (2019). Deterministic parallel analysis: an improved method for selecting factors and principal components. J. Roy. Statist. Soc. Ser. B 81 163–183.
  * El Karoui (2007) El Karoui, N. (2007). Tracy-Widom limit for the largest eigenvalue of a large class of complex sample covariance matrices. Ann. Probab. 35 663–714.
  * Erdös et al. (2013) Erdös, L., Knowles, A. and Yau, H.-T. (2013). Averaging fluctuations in resolvents of random band matrices. Ann. Henri Poincaáe 14 1837–1926.
  * Benaych-Georges and Knowles (2016) Benaych-Georges, F., Knowles, A. (2016). Lectures on the local semicircle law for Wigner matrices. https://doi.org/10.48550/arXiv.1601.04055
  * Erdös et al. (2012) Erdös, L., Yau, H.-T. and Yin, J. (2012). Rigidity of eigenvalues of generalized Wigner matrices. Adv. Math. 229(3) 1435–1515.
  * Gao et al. (2017) Gao, J., Han, X., Pan, G. M. and Yang, Y. R (2017). High dimensional correlation matrices: CLT and its applications. J. Roy. Statist. Soc. Ser. B 79 677-693.
  * Hu et al. (2019) Hu, J., Li, W. M., Liu, Z. and Zhou, W. (2019). High-dimensional covariance matrices in elliptical distributions with application to spherical test. Ann. Statist. 47 527–555.
  * Hu et al. (2022) Hu, J., Wang, S. C., Zhang, Y. C. and Zhou, W. (2022). Sampling without replacement from a high-dimensional finite population supplementary material.
  * Johansson (2000) Johansson, K. (2000). Shape fluctuations and random matrices. Comm. Math. Phys. 209 437–476.
  * Johnstone (2001) Johnstone, I. M. (2001). On the distribution of the largest eigenvalue in principal components analysis. Ann. Stat. 29 295–327.
  * Knowles and Yin (2017) Knowles, A. and Yin, J. (2017). Anisotropic local laws for random matrices. Probab. Theory Relat. Fields 169(1) 257–352.
  * Lee and Schnelli (2016) Lee, J. O. and Schnelli, K. (2016) Tracy-Widom distribution for the largest eigenvalue of real sample covariance matrices with general population. Ann. Appl. Probab., 26(6), 3786-3839.
  * Li et al. (2008) Li, J. Z., Absher, D. M., Tang, H., Southwick, A. M., Casto, A. M., Ramachandran, S., Cann, H. M., Barsh, G. S., Feldman, M., Cavalli-Sforza, L. L. and Myers, R. M. (2008) Worldwide human relationships inferred from genome-wide patterns of variation. Science 319 1100–1104.
  * Marcenko and Pastur (1967) Marcenko, V. A. and Pastur, L. A. (1967). Distribution of eigenvalues for some sets of random matrices. Mathematics of the USSR-Sbornik 1(4) 457\.
  * Onatski (2008) Onatski, A. (2008). The Tracy-Widom limit for the largest eigenvalues of singular complex Wishart matrices. Ann. Appl. Probab. 18 470–490.
  * Pillai and Yin (2012) Pillai, N. S. and Yin, J. (2012). Edge universality of correlation matrices. Ann. Stat. 40(3) 1737–1763.
  * Pillai and Yin (2014) Pillai, N. S. and Yin, J. (2014). Universality of covariance matrices. Ann. Appl. Probab. 24(3) 935–1001.
  * Silverstein and Choi (1995) Silverstein, J. and Choi, S. (1995). Analysis of the limiting spectral distribution of large dimensional random matrices. J. Multivariate Anal. 54(2) 295–309.
  * Tian et al. (2021) Tian, B. P., Zhang, Y. C. and Zhou, W.  (2021). Tracy-Widom law for the largest eigenvalue of sample covariance matrix generated by VARMA Random Matrices: Theory and Applications. 10 (02) 2150022
  * Wishart (1928) Wishart, J. (1928) The generalised product moment distribution in samples from a normal multivariate population. Biometrika 20A(1/2) 32–52.
  * Wen et al. (2020) Wen, J., Xie, J. H., Yu, L. and Zhou, W. (2020). Tracy-Widom limit for the largest eigenvalue of high-dimensional covariance matrices in elliptical distributions.  arXiv:1901.05166v2.