An `∞ Eigenvector Perturbation Bound and Its
Application to Robust Covariance Estimation
arXiv:1603.03516v2 [math.ST] 2 Jun 2017

Jianqing Fan∗, Weichen Wang and Yiqiao Zhong
Department of Operations Research and Financial Engineering, Princeton University

Abstract
In statistics and machine learning, we are interested in the eigenvectors (or singular
vectors) of certain matrices (e.g. covariance matrices, data matrices, etc). However,
those matrices are usually perturbed by noises or statistical errors, either from random
sampling or structural patterns. The Davis-Kahan sin θ theorem is often used to bound
the difference between the eigenvectors of a matrix A and those of a perturbed matrix
e = A + E, in terms of `2 norm. In this paper, we prove that when A is a low-rank and
A
incoherent matrix, the `∞ norm perturbation bound of singular vectors (or eigenvectors
√
√
in the symmetric case) is smaller by a factor of d1 or d2 for left and right vectors,
where d1 and d2 are the matrix dimensions. The power of this new perturbation result
is shown in robust covariance estimation, particularly when random variables have
heavy tails. There, we propose new robust covariance estimators and establish their
asymptotic properties using the newly developed perturbation bound. Our theoretical
results are verified through extensive numerical experiments.

Keywords: Matrix perturbation theory, Incoherence, Low-rank matrices, Sparsity, Approximate factor model.

∗

Address: Department of ORFE, Sherrerd Hall, Princeton University, Princeton, NJ 08544, USA, e-mail:
jqfan@princeton.edu, yiqiaoz@princeton.edu, weichenw@princeton.edu. The research was partially supported
by NSF grants DMS-1206464 and DMS-1406266 and NIH grants R01-GM072611-10.

1

1

Introduction

The perturbation of matrix eigenvectors (or singular vectors) has been well studied in
matrix perturbation theory (Wedin, 1972; Stewart, 1990). The best known result of eigenvector perturbation is the classic Davis-Kahan theorem (Davis and Kahan, 1970). It originally
emerged as a powerful tool in numerical analysis, but soon found its widespread use in other
fields, such as statistics and machine learning. Its popularity continues to surge in recent
years, which is largely attributed to the omnipresent data analysis, where it is a common
practice, for example, to employ PCA (Jolliffe, 2002) for dimension reduction, feature extraction, and data visualization.
The eigenvectors of matrices are closely related to the underlying structure in a variety
of problems. For instance, principal components often capture most information of data
and extract the latent factors that drive the correlation structure of the data (Bartholomew
et al., 2011); in classical multidimensional scaling (MDS), the centered squared distance
matrix encodes the coordinates of data points embedded in a low dimensional subspace
(Borg and Groenen, 2005); and in clustering and network analysis, spectral algorithms are
used to reveal clusters and community structure (Ng et al., 2002; Rohe et al., 2011). In
those problems, the low dimensional structure that we want to recover, is often ‘perturbed’
by observation uncertainty or statistical errors. Besides, there might be a sparse pattern
corrupting the low dimensional structure, as in approximate factor models (Chamberlain
and Rothschild, 1982; Stock and Watson, 2002) and robust PCA (De La Torre and Black,
2003; Candès et al., 2011).
A general way to study these problems is to consider
e = A + S + N,
A

(1)

where A is a low rank matrix, S is a sparse matrix, and N is a random matrix regarded
as random noise or estimation error, all of which have the same size d1 × d2 . Usually A is
regarded as the ‘signal’ matrix we are primarily interested in, S is some sparse contamination
whose effect we want to separate from A, and N is the noise (or estimation error in covariance
matrix estimation).
The decomposition (1) forms the core of a flourishing literature on robust PCA (Chandrasekaran et al., 2011; Candès et al., 2011), structured covariance estimation (Fan et al.,
2008, 2013), multivariate regression (Yuan et al., 2007) and so on. Among these works, a
standard condition on A is matrix incoherence (Candès et al., 2011). Let the singular value

2

decomposition be
A = U ΣV

T

=

r
X

σi ui viT ,

(2)

i=1

where r is the rank of A, the singular values are σ1 ≥ σ2 ≥ . . . ≥ σr > 0, and the matrices
U = [u1 , . . . , ur ] ∈ Rd1 ×r , V = [v1 , . . . , vr ] ∈ Rd2 ×r consist of the singular vectors. The
coherences µ(U ), µ(V ) are defined as
r

r

X
d1
U2 ,
max
µ(U ) =
r i j=1 ij

X
d2
µ(V ) =
V 2,
max
r i j=1 ij

(3)

where Uij and Vij are the (i, j) entry of U and V , respectively. It is usually expected that
µ0 := max{µ(U ), µ(V )} is not too large, which means the singular vectors ui and vi are
incoherent with the standard basis. This incoherence condition (3) is necessary for us to
separate the sparse component S from the low rank component A; otherwise A and S are
not identifiable. Note that we do not need any incoherence condition on U V T , which is
different from Candès et al. (2011) and is arguably unnecessary (Chen, 2015).
Now we denote the eigengap γ0 = min{σi − σi+1 : i = 1, . . . , r} where σr+1 := 0 for
notational convenience. Also we let E = S + N , and view it as a perturbation matrix
to the matrix A in (1). To quantify the perturbation, we define a rescaled measure as
p
p
τ0 := max{ d2 /d1 kEk1 , d1 /d2 kEk∞ }, where
kEk1 = max
j

d1
X

|Eij |,

kEk∞ = max
i

i=1

d2
X

|Eij |,

(4)

j=1

which are commonly used norms gauging sparsity (Bickel and Levina, 2008). They are also
p
operator norms in suitable spaces (see Section 2). The rescaled norms d2 /d1 kEk1 and
p
d1 /d2 kEk∞ are comparable to the spectral norm kEk2 := maxkuk2 =1 kEuk2 in many cases;
p
p
for example, when E is an all-one matrix, d2 /d1 kEk1 = d1 /d2 kEk∞ = kEk2 .
e also has the singular value decomposition:
Suppose the perturbed matrix A
e=
A

dX
1 ∧d2

σ
ei u
ei veiT ,

(5)

i=1

where σ
ei are nonnegative and in the decreasing order, and the notation ∧ means a ∧ b =
e = [e
min{a, b}. Denote U
u1 , . . . , u
er ], V = [e
v1 , . . . , ver ], which are counterparts of top r singular
vectors of A.
We will present an `∞ matrix perturbation result that bounds ke
ui − ui k∞ and ke
vi − vi k∞

3

up to sign.1 This result is different from `2 bounds, Frobenius-norm bounds, or the sin Θ
bounds, as the `∞ norm is not orthogonal invariant. The following theorem is a simplified
version of our main results in Section 2.
e = A+E and suppose the singular decomposition in (2) and (5). Denote
Theorem 1.1. Let A
γ0 = min{σi − σi+1 : i = 1, . . . , r} where σr+1 := 0. Then there exists C(r, µ0 ) = O(r4 µ20 )
such that, if γ0 > C(r, µ0 )τ0 , up to sign,
max ke
ui − ui k∞ ≤ C(r, µ0 )

1≤i≤r

τ
√0

γ0 d1

and

max ke
vi − vi k∞ ≤ C(r, µ0 )

1≤i≤r

where µ0 = max{µ(U ), µ(V )} is the coherence given after
p
p
max{ d2 /d1 kEk1 , d1 /d2 kEk∞ }.

τ
√0

γ0 d2

,

(3) and τ0

(6)
:=

When A is symmetric, the condition on the eigengap is simply γ0 > C(r, µ0 )kEk∞ . It
naturally holds for a variety of applications, where the low rank structure emerges as a
consequence of a few factors driving the data matrix. For example, in Fama-French factor
models, the excess returns in a stock market are driven by a few common factors (Fama and
French, 1993); in collaborative filtering, the ratings of users are mostly determined by a few
common preferences (Rennie and Srebro, 2005); in video surveillance, A is associated with
the stationary background across image frames (Oliver et al., 2000). We will have a detailed
discussion in Section 2.3.
The eigenvector perturbation was studied by Davis and Kahan (1970), where Hermitian
matrices were considered, and the results were extended by Wedin (1972) to general rectangular matrices. To compare our result with these classical results, assuming γ0 ≥ 2kEk2 , a
combination of Wedin’s theorem and Mirsky’s inequality (Mirsky, 1960) (the counterpart of
Weyl’s inequality for singular values) implies
√

2 2kEk2
.
max kvk − vek k2 ∨ kuk − u
ek k2 ≤
1≤k≤r
γ0

(7)

where a ∨ b := max{a, b}.
Yu et al. (2015) also proved a similar bound as in (7), and that result is more convenient to
use. If we are interested in the `∞ bound but naively use the trivial inequality kxk∞ ≤ kxk2 ,
we would have a suboptimal bound O(kEk2 /γ0 ) in many situations, especially in cases where
√
kEk2 is comparable to kEk∞ . Compared with (6), the bound is worse by a factor of d1
√
for uk and d2 for vk . In other words, converting the `2 bound from Davis-Kahan theorem
directly to the `∞ bound does not give a sharp result in general, in the presence of incoherent
1

‘Up to sign’ means we can appropriately choose an eigenvector or singular vector u to be either u or
−u in the bounds. This is becuase eigenvectors and singular vectors are not unique.

4

and low rank structure of A. Actually, assuming kEk2 is comparable with kEk∞ , for square
matrices, our `∞ bound (6) matches the `2 bound (7) in terms of dimensions d1 and d2 . This
√
√
√
is because kxk2 ≤ n kxk∞ for any x ∈ Rn , so we expect to gain a factor d1 or d2 in
those `∞ bounds. The intuition is that, when A has an incoherent and low-rank structure,
the perturbation of singular vectors is not concentrated on a few coordinates.
To understand how matrix incoherence helps, let us consider a simple example with no
matrix incoherence, in which (7) is tight up to a constant. Let A = d(1, 0, . . . , 0)T (1, 0, . . . , 0)
be a d-dimensional square matrix, and E = d(0, 1/2, 0, . . . , 0)T (1, 0, . . . , 0) of the same size. It
√
√
is apparent that γ0 = d, τ0 = d/2, and that v1 = (1, 0, . . . , 0)T , ve1 = (2/ 5, 1/ 5, 0, . . . , 0)T
up to sign. Clearly, the perturbation ke
v1 − v1 k∞ is not vanishing as d tends to infinity in
this example, and thus, there is no hope of a strong upper bound as in (6) without the
incoherence condition.
√
√
ek
The reason that the factor d1 or d2 comes into play in (7) is that, the error uk − u
(and similarly for vk ) spreads out evenly in d1 (or d2 ) coordinates, so that the `∞ error is
far smaller than the `2 error. This, of course, hinges on the incoherence condition, which in
essence precludes eigenvectors from aligning with any coordinate.
Our result is very different from the sparse PCA literature, in which it is usually assumed
that the leading eigenvectors are sparse. In Johnstone and Lu (2009), it is proved that there
is a threshold for p/n (the ratio between the dimension and the sample size), above which
PCA performs poorly, in the sense that he
v1 , v1 i is approximately 0. This means that the
principal component computed from the sample covariance matrix reveals nothing about the
true eigenvector. In order to mitigate this issue, in Johnstone and Lu (2009) and subsequent
papers (Vu and Lei, 2012; Ma, 2013; Berthet and Rigollet, 2013), sparse leading eigenvectors
are assumed. However, our result is different, in the sense that we require a stronger eigengap
condition γ0 > C(r, µ0 )kEk∞ (i.e. stronger signal), whereas in Johnstone and Lu (2009), the
eigengap of the leading eigenvectors is a constant times kEk2 . This explains why it is
plausible to have a strong uniform eigenvector perturbation bound in this paper.
We will illustrate the power of this perturbation result using robust covariance estimation
as one application. In the approximate factor model, the true covariance matrix admits a
decomposition into a low rank part A and a sparse part S. Such models have been widely
applied in finance, economics, genomics, and health to explore correlation structure.
However, in many studies, especially financial and genomics applications, it is well known
that the observations exhibit heavy tails (Gupta et al., 2013). This problem can be resolved
with the aid of recent results of concentration bounds in robust estimation (Catoni, 2012;
Hsu and Sabato, 2014; Fan et al., 2017), which produces the estimation error N in (1) with
an optimal entry-wise bound. It nicely fits our perturbation result, and we can tackle it
5

easily by following the ideas in Fan et al. (2013).
Here are a few notations in this paper. For a generic d1 by d2 matrix, the matrix maxnorm is denoted as kM kmax = maxi,j |Mij |. The matrix operator norm induced by vector `p
P1
norm is kM kp = supkxkp =1 kM xkp for 1 ≤ p ≤ ∞. In particular, kM k1 = maxj di=1
|Mij |;
Pd2
kM k∞ = maxi j=1 |Mij |; and k·k denotes the spectral norm, or the matrix 2-norm k·k2 for
simplicity. We use σj (M ) to denote the j th largest singular value. For a symmetric matrix
M , denote λj (M ) as its j th largest eigenvalue. If M is a positive definite matrix, then M 1/2
is the square root of M , and M −1/2 is the square root of M −1 .

2

The `∞ perturbation result

2.1

Symmetric matrices

First, we study `∞ perturbation for symmetric matrices (so d1 = d2 ). The approach we
study symmetric matrices will be useful to analyze asymmetric matrices, because we can
always augment a d1 × d2 rectangular matrix into a (d1 + d2 ) × (d1 + d2 ) symmetric matrix,
and transfer the study of singular vectors to the eigenvectors of the augmented matrix. This
augmentation is called Hermitian dilation. (Tropp, 2012; Paulsen, 2002)
Suppose that A ∈ Rd×d is an d-dimensional symmetric matrix. The perturbation matrix
e := A + E.
E ∈ Rd×d is also d-dimensional and symmetric. Let the perturbed matrix be A
Suppose the spectral decomposition of A is given by
A = [V, V⊥ ]

Λ1 0
0 Λ2

!
[V, V⊥ ]T =

r
X

λi vi viT +

i=1

X

λi vi viT ,

(8)

i>r

where Λ1 = diag{λ1 , . . . , λr }, Λ2 = diag{λr+1 , . . . , λn }, and where |λ1 | ≥ |λ2 | ≥ . . . ≥ |λn |.
P
Note the best rank-r approximation of A under the Frobenius norm is Ar := i≤r λi vi viT .2
e is
Analogously, the spectral decomposition of A
e=
A

r
X

ei vei veT +
λ
i

i=1

X

ei vei veT ,
λ
i

i>r

e1 | ≥ |λ
e 2 | ≥ . . . ≥ |λ
en |. Recall that kEk∞ given
and write Ve = [e
v1 , . . . , ver ] ∈ Rd×r , where |λ
by (4) is an operator norm in the `∞ space, in the sense that kEk∞ = supkuk∞ ≤1 kEuk∞ .
This norm is the natural counterpart of the spectral norm kEk2 := supkuk2 ≤1 kEuk2 .

2

This is a consequence of Wielandt-Hoffman theorem.

6

We will use notations O(·) and Ω(·) to hide absolute constants.3 The next theorem
bounds the perturbation of eigenspaces up to a rotation.
Theorem 2.1. Suppose |λr | − ε = Ω(r3 µ2 kEk∞ ), where ε = kA − Ar k∞ , which is the
approximation error measured under the matrix ∞-norm and µ = µ(V ) is the coherence of
V defined in (3). Then, there exists an orthogonal matrix R ∈ Rr×r such that

kVe R − V kmax = O

r5/2 µ2 kEk∞
√
(|λr | − ε) d


.

This result involves an unspecified rotation R, due to the possible presence of multiplicity
of eigenvalues. In the case where λ1 = · · · = λr > 0, the individual eigenvectors of V are
only identifiable up to rotation. However, assuming an eigengap (similar to Davis-Kahan
theorem), we are able to bound the perturbation of individual eigenvectors (up to sign).
Theorem 2.2. Assume the conditions in Theorem 2.1. In addition, suppose δ satisfies
δ > kEk2 , and for any i ∈ [r], the interval [λi − δ, λi + δ] does not contain any eigenvalues
of A other than λi . Then, up to sign,

max ke
vi − vi k∞
i∈[r]

= kVe − V kmax = O

r4 µ2 kEk∞
r3/2 µ1/2 kEk2
√ +
√
(|λr | − ε) d
δ d


.

To understand the above two theorems, let us consider the case where A has exactly
rank r (i.e., ε = 0), and r and µ are not large (say, bounded by a constant). Theorem
√
2.1 gives a uniform entrywise bound O(kEk∞ /|λr | d) on the eigenvector perturbation. As
a comparison, the Davis–Kahan sin Θ theorem (Davis and Kahan, 1970) gives a bound
√
O(kEk2 /|λr |) on kVe R − V k2 with suitably chosen rotation R.4 This is an order of d larger
than the bound given in Theorem 2.1 when kEk∞ is of the same order as kEk2 . Thus, in
scenarios where kEk2 is comparable to kEk∞ , this is a refinement of Davis-Kahan theorem,
because the max-norm bound in Theorem 2.1 provides an entry-wise control of perturbation.
Although kEk∞ ≥ kEk2 ,5 there are many settings where the two quantities are comparable;
for example, if E has a submatrix whose entries are identical and has zero entries otherwise,
then kEk∞ = kEk2 .
Theorem 2.2 provides the perturbation of individual eigenvectors, under a usual eigengap
3

We write a = O(b) if there is a constant C > 0 such that a < Cb; and a = Ω(b) if there is a constant
C > 0 such that a > C 0 b.
4
To see how the Davis-Kahan sin Θ theorem relates to this form, we can use the identity k sin Θ(Ve , V )k2 =
T
kVe Ve −V V T k2 (Stewart, 1990), and the (easily verifiable) inequality 2 minR kVe R−V k2 ≥ kVe Ve T −V V T k2 ≥
minR kVe R − V k2 where R is an orthogonal matrix.
5
Since kEk1 kEk∞ ≤ kEk22 (Stewart, 1990), the inequality follows from kEk1 = kEk∞ by symmetry.
0

7

√
assumption. When r and µ are not large, we incur an additional term O(kEk2 /δ d) in the
bound. This is understandable, since ke
vi − vi k2 is typically O(kEk2 /δ).
When the rank of A is not exactly r, we require that |λr | is larger than the approximation
error kA − Ar k∞ . It is important to state that this assumption is more restricted than the
eigengap assumption in the Davis-Kahan theorem, since kA − Ar k∞ ≥ kA − Ar k2 = |λr+1 |.
However, different from the matrix max-norm, the spectral norm k · k2 only depends on the
eigenvalues of a matrix, so it is natural to expect `2 perturbation bounds that only involve
λr and λr+1 . It is not clear whether we should expect an `∞ bound that involves λr+1 instead
of ε. More discussions can be found in Section 5.
We do not pursue the optimal bound in terms of r and µ(V ) in this paper, as the two
quantities are not large in many applications, and the current proof is already complicated.

2.2

Rectangular matrices

Now we establish `∞ perturbation bounds for general rectangular matrices. The results
here are more general than those in Section 1, and in particular, we allow the matrix A to be
e := A + E.
of approximate low rank. Suppose that both A and E are d1 × d2 matrices, and A
The rank of A is at most d1 ∧ d2 (where a ∧ b = min{a, b}). Suppose an integer r satisfies
r ≤ rank(A). Let the singular value decomposition of A be

A=

r
X

σi ui viT

+

i=1

dX
1 ∧d2

σi ui viT ,

i=r+1

where the singular values are ordered as σ1 ≥ σ2 ≥ . . . ≥ σd1 ∧d2 ≥ 0, and the unit vectors u1 , . . . , ud1 ∧d2 (or unit vectors v1 , . . . , vd1 ∧d2 ) are orthogonal to each other. We denote
U = [u1 , . . . , ur ] ∈ Rd1 ×r and V = [v1 , . . . , vr ] ∈ Rd2 ×r . Analogously, the singular value
e is
decomposition of A
dX
r
1 ∧d2
X
T
e
A=
σ
ei u
ei vei +
σ
ei u
ei veiT ,
i=1

i=r+1

e = [e
where σ
e1 ≥ . . . ≥ σ
ed1 ∧d2 . Similarly, columns of U
u1 , . . . , u
er ] ∈ Rd1 ×r and V =
[e
v1 , . . . , ver ] ∈ Rd2 ×r are orthonormal.
Define µ0 = max{µ(V ), µ(U )}, where µ(U ) (resp. µ(V )) is the coherence of U (resp. V ).
This µ0 will appear in the statement of our results, as it controls both the structure of left and
right singular spaces. When, specially, A is a symmetric matrix, the spectral decomposition
of A is also the singular value decomposition (up to sign), and thus µ0 coincides with µ
defined in Section 2.1.
Recall the definition of matrix ∞-norm and 1-norm of a rectangular matrix (4). Similar
8

to the matrix ∞-norm, k · k1 is an operator norm in the `1 space. An obvious relationship
between matrix ∞-norm and 1-norm is kEk∞ = kE T k1 . Note that the matrix ∞-norm
and 1-norm have different number of summands in their definitions, so we are motivated to
p
p
consider τ0 := max{ d1 /d2 kEk∞ , d2 /d1 kEk1 } to balance the dimensions d1 and d2 .
P
Let Ar = i≤r σi ui viT be the best rank-r approximation of A under the Frobenius norm,
p
p
and let ε0 = d1 /d2 kA − Ar k∞ ∨ d2 /d1 kA − Ar k1 , which also balances the two dimensions.
Note that in the special case where A is symmetric, this approximation error ε0 is identical
to ε defined in Section 2.1. The next theorem bounds the perturbation of singular spaces.
Theorem 2.3. Suppose that δ0 − ε0 = Ω(r3 µ20 τ0 ). Then, there exists orthogonal matrices
RU , RV ∈ Rr×r such that,
e RU − U kmax = O
kU



r5/2 µ20 τ0 
√
,
(σr − 0 ) d1

kVe RV − V kmax = O



r5/2 µ20 τ0 
√
.
(σr − 0 ) d2

Similar to Theorem 2.2, under an assumption of gaps between singular values, the next
theorem bounds the perturbation of individual singular vectors.
Theorem 2.4. Suppose the same assumption in Theorem 2.3. In addition, suppose δ0
satisfies δ0 > kEk2 , and for any i ∈ [r], the interval [σi − δ0 , σi + δ0 ] does not contain any
eigenvalues of A other than σi . Then, up to sign,

i∈[r]

max ke
vi − vi k∞
i∈[r]

1/2

r3/2 µ0 kEk2 
r4 µ20 τ0
√ +
√
,
(σr − ε0 ) d1
δ0 d 1
1/2

r4 µ20 τ0
r3/2 µ0 kEk2 
√ +
√
=O
.
(σr − ε0 ) d2
δ0 d 2

max ke
ui − ui k∞ = O



(9)
(10)

As mentioned in the beginning of this section, we will use dilation to augment all d1 × d2
matrices into symmetric ones with size d1 + d2 . In order to balance the possibly different
scales of d1 and d2 , we consider a weighted max-norm. This idea will be further illustrated
in Section 5.

2.3

Examples: which matrices have such structure?

In many problems, low-rank structure naturally arises due to the impact of pervasive latent factors that influence most observed data. Since observations are imperfect, the low-rank
structure is often ‘perturbed’ by an additional sparse structure, gross errors, measurement
noises, or the idiosyncratic components that can not be captured by the latent factors. We
give some motivating examples with such structure.
9

Panel data in stock markets. Consider the excess returns from a stock market over a
period of time. The driving factors in the market are reflected in the covariance matrix as
a low rank component A. The residual covariance of the idiosyncratic components is often
modeled by a sparse component S. Statistical analysis including PCA is usually conducted
e = Σ,
b which is perturbed from the true covariance
based on the estimated covariance matrix A
Σ = A + S by the estimation error N (Stock and Watson, 2002; Fan et al., 2013). In Section
3.1, we will develop a robust estimation method in the presence of heavy-tailed return data.
Video surveillance. In image processing and computer vision, it is often desired to separate moving objects from static background before further modeling and analysis (Oliver
et al., 2000; Hu et al., 2004). The static background corresponds to the low rank component
A in the data matrix, which is a collection of video frames, each consisting of many pixels
represented as a long vector in the data matrix. Moving objects and noise correspond to
the sparse matrix S and noise matrix N . Since the background is global information and
reflected by many pixels of a frame, it is natural for the incoherence condition to hold.
Wireless sensor network localization. In wireless sensor networks, we are usually interested in determining the location of sensor nodes with unknown position based on a few
(noisy) measurements between neighboring nodes (Doherty et al., 2001; Biswas and Ye,
2004). Let X be an r by n matrix such that each column xi gives the coordinates of each
node in a plane (r = 2) or a space (r = 3). Assume the center of the sensors has been
relocated at origin. Then the low rank matrix A = XT X, encoding the true distance information, has to satisfy distance constraints given by the measurements. The noisy distance
e after centering, equals to the sum of A and a matrix N consisting of measurement
matrix A
errors. Suppose that each node is a random point uniformly distributed in a rectangular
region. It is not difficult to see that with high probability, the top r eigenvalues of XT X
and their eigengap scales with the number of sensors n and the leading eigenvectors have a
bounded coherence.
In our theorems, we require that the coherence µ is not too large. This is a natural
structural condition associated with the low rank matrices. Consider the following very
simple example: if the eigenvectors v1 , . . . , vr of the low rank matrix A are uniform unit
√
vectors in a sphere, then with high probability, maxi kvi k∞ = O( log n), which implies µ =
O(log n). An intuitive way to understand the incoherence structure is that no coordinates
of v1 (or v2 , . . . vr ) are dominant. In other words, the eigenvectors are not concentrated on a
few coordinates.
In all our examples, the incoherence structure is natural. The factor model satisfies such
structure, which will be discussed in Section 3. In the video surveillance example, ideally,
when the images are static, A is a rank one matrix x1T . Since usually a majority of pixels
10

(coordinates of x) help to display an image, the vector x often has dense coordinates with
comparable magnitude, so A also has an incoherence structure in this example. Similarly,
in the sensor localization example, the coordinates of all sensor nodes are comparable in
magnitude, so the low rank matrix A formed by XT X also has the desired incoherence
structure.

2.4

Other perturbation results

Although the eigenvector perturbation theory is well studied in numerical analysis, there
is a renewed interest among statistics and machine learning communities recently, due to
the wide applicability of PCA and other eigenvector-based methods. In Cai and Zhang
(2016); Yu et al. (2015), they obtained variants or improvements of Davis-Kahan theorem
(or Wedin’s theorem), which are user-friendly in the statistical contexts. These results
assume the perturbation is deterministic, which is the same as Davis-Kahan theorem and
Wedin’s theorem. In general, these results are sharp, even when the perturbation is random,
as evidenced by the BBP transition (Baik et al., 2005).
However, these classical results can be suboptimal, when the perturbation is random and
the smallest eigenvalue gap λ1 − λ2 does not capture particular spectrum structure. For
example, Vu (2011); O’Rourke et al. (2013) showed that with high probability, there are
bounds sharper than the Wedin’s theorem, when the signal matrix is low-rank and satisfies
certain eigenvalue conditions.
In this paper, our perturbation results are deterministic, thus the bound can be suboptimal when the perturbation is random with certain structure (e.g. the difference between
sample covariance and population one for i.i.d. samples). However, the advantage of a deterministic result is that it is applicable to any random perturbation. This is especially useful
when we cannot make strong random assumptions on the perturbation (e.g., the perturbation
is an unknown sparse matrix). In Section 3, we will see examples of this type.

3

Application to robust covariance estimation

We will study the problem of robust estimation of covariance matrices and show the
strength of our perturbation result. Throughout this section, we assume both rank r and
the coherence µ(V ) are bounded by a constant, though this assumption can be relaxed. We
will use C to represent a generic constant, and its value may change from line to line.

11

3.1

PCA in spiked covariance model

To initiate our discussions, we first consider sub-Gaussian random variables. Let X =
(X1 , . . . , Xd ) be a random d-dimensional vector with mean zero and covariance matrix
Σ=

r
X

λi vi viT + σ 2 Id := Σ1 + Σ2 ,

(λ1 ≥ . . . ≥ λr > 0),

(11)

i=1

and X be an n by d matrix, whose rows are independently sampled from the same distribution. This is the spiked covariance model that has received intensive study in recent years.
b = XT X/n. Viewing the empirical covariance matrix
Let the empirical covariance matrix be Σ
as its population version plus an estimation error, we have the decomposition

1
T
b
X X−Σ ,
Σ = Σ1 + Σ2 +
n
which is a special case of the general decomposition in (1). Here, Σ2 is the sparse component,
and the estimation error XT X/n−Σ is the noise component. Note that v1 , . . . , vr are just the
top r leading eigenvectors of Σ and we write V = [v1 , . . . , vr ]. Assume the top r eigenvectors
b are denoted by vb1 , . . . , vbr . We want to find an `∞ bound on the estimation error vbi − vi
of Σ
for all i ∈ [r].
When the dimension d is comparable to or larger than n, it has been shown by Johnstone
and Lu (2009) that the leading empirical eigenvector vb1 is not a consistent estimate of the true
eigenvector v1 , unless we assume larger eigenvalues. Indeed, we will impose more stringent
conditions on λi ’s in order to obtain good `∞ bounds.
Assuming the coherence µ(V ) is bounded, we can easily see Var(Xj ) ≤ σ 2 + Cλ1 /d for
some constant C. It follows from the standard concentration result (e.g., Vershynin (2010))
that if rows of X contains i.i.d sub-Gaussian vectors and log d = O(n), then with probability
greater than 1 − d−1 ,
r
λ1  log d
1 T
2
.
(12)
k X X − Σkmax ≤ C σ +
n
d
n
b − Σ1 as E. If the conditions in Theorem 2.2
To apply Theorem 2.2, we treat Σ1 as A and Σ
are satisfied, we will obtain
√
√
max kb
vk − vk k∞ = O(kEk∞ /(λr d) + kEk2 /(δ d)).

1≤k≤r

12

(13)

Note there are simple bounds on kEk∞ and kEk2 :
kEk2 ≤ kEk∞

r
n

log d o
1
≤ σ 2 + d k XT X − Σkmax ≤ C 1 + dσ 2 + λ1
.
n
n

By assuming a strong uniform eigengap, the conditions in Theorem 2.2 are satisfied, and the
bound in (13) can be simplified. Define the uniform eigengap as
γ = min{λi − λi+1 : 1 ≤ i ≤ r},
Note that γ ≤ min{λr , δ}, so if γ > C(1 + dσ 2 + λ1
max kb
vk − vk k∞

1≤k≤r

λr+1 := 0.

p
log d/n), we have

 1 + dσ 2 + λ plog d/n 
 kEk 
1
∞
√
√
= OP
,
= OP
γ d
γ d

p
In particular, when λ1  γ and γ  max{1, σ 2 d log d/n}, we have
 1 
max kb
vk − vk k∞ = oP √ .
1≤k≤r
d
The above analysis pertains to the structure of sample covariance matrix. In the following
subsections, we will estimate the covariance matrix using more complicated robust procedure.
Our perturbation theorems in Section 2 provide a fast and clean approach to obtain new
results.

3.2

PCA for robust covariance estimation

The usefulness of Theorem 2.2 is more pronounced when the random variables are heavytailed. Consider again the covariance matrix Σ with structure (11). Instead of assuming subGaussian distribution, we assume there exists a constant C > 0 such that maxj≤d EXj4 < C,
i.e. the fourth moments of the random variables are uniformly bounded.
Unlike sub-Gaussian variables, there is no concentration bound similar to (12) for the
empirical covariance matrix. Fortunately, thanks to recent advances in robust statistics
(e.g., Catoni (2012)), robust estimate of Σ with guaranteed concentration property becomes
possible. We shall use the method proposed in Fan et al. (2017). Motivated by the classical
M -estimator of Huber (1964), Fan et al. (2017) proposed a robust estimator for each element
b by solving a Huber loss based minimization problem
of Σ,
b ij = argminµ
Σ

n
X
t=1

13

lα (Xti Xtj − µ),

(14)

where lα is the Huber loss defined as

2α|x| − α2 , |x| ≥ α,
lα (x) =
x2 ,
|x| ≤ α.
p
The parameter α is suggested to be α = nv 2 / log(−1 ) for  ∈ (0, 1), where v is assumed
p
to satisfy v ≥ maxij Var(Xi Xj ). If log(−1 ) ≤ n/8, Fan et al. (2017) showed


r

b ij − Σij | ≤ 4v
P |Σ

log(−1 ) 
≥ 1 − 2.
n

From this result, the next proposition is immediate by taking  = d−3 .
Proposition 3.1. Suppose that there is a constant C with maxj≤d EXj4 < C. Then with
probability greater than 1 − d−1 (1 + d−1 ), the robust estimate of covariance matrix with
p
α = 3nv 2 log(d) satisfies
r
b − Σkmax ≤ 4v 3 log d ,
kΣ
n
p
where v is a pre-determined parameter assumed to be no less than maxij Var(Xi Xj ).
This result relaxes the sub-Gaussianity assumption by robustifying the covariance estimate. It is apparent that the `∞ bound in the previous section is still valid in this case. To
be more specific, suppose µ(V ) is bounded by a constant. Then, (13) holds for the PCA
p
based on the robust covariance estimation. When λ1  γ and γ  max{1, σ 2 d log d/n},
we again have
max kb
vk − vk k∞

1≤k≤r

 1 + dσ 2 + λ plog d/n 
 1 
1
√
= OP
= oP √ .
γ d
d

√
Note that an entrywise estimation error op (1/ d) necessarily implies consistency of the
estimated eigenvectors, since we can easily convert an `∞ result into an `2 result. The
minimum signal strength (or magnitude of leading eigenvalues) for such consistency is shown
to be σ 2 d/n under the sub-Gaussian assumption (Wang and Fan, 2017+).
If the goal is simply to prove consistency of vbk , the strategy of using our `∞ perturbation
bounds is not optimal. However, there are also merits: our result is nonasymptotic; it holds
for more general distributions (beyond sub-Gaussian distributions); and its entrywise bound
gives stronger guarantee. Moreover, the `∞ perturbation bounds provide greater flexibility
for analysis, since it is straightforward to adapt analysis to problems with more complicated
structure. For example, the above discussion can be easily extended to a general Σ2 with
14

bounded kΣ2 k∞ rather than a diagonal matrix.

3.3

Robust covariance estimation via factor models

In this subsection, we will apply Theorem 2.2 to robust large covariance matrix estimation
for approximate factor models in econometrics. With this theorem, we are able to extend the
data distribution in factor analysis beyond exponentially decayed distributions considered
by Fan et al. (2013), to include heavy-tailed distributions.
Suppose the observation yit , say, the excess return at day t for stock i, admits a decomposition
i ≤ d, t ≤ n,
(15)
yit = bTi ft + uit ,
where bi ∈ Rr is the unknown but fixed loading vector, ft ∈ Rr denotes the unobserved
factor vector at time t, and uit ’s represent the idiosyncratic noises. Let yt = (y1t , . . . , ydt )T
and ut = (u1t , . . . , udt )T so that yt = Bft + ut , where B = (b1 , . . . , bd )T ∈ Rd×r . Suppose that
ft and ut are uncorrelated and centered random vectors, with bounded fourth moments, i.e.,
the fourth moments of all entries of ft and ut are bounded by some constant. We assume
{ft , ut } are independent for t, although it is possible to allow for weak temporal dependence
as in Fan et al. (2013). From (15), we can decompose Σ = Cov(yt ) into a low rank component
and a residual component:
Σ = BB T + Σu ,
(16)
where Σu := Cov(ut ). To circumvent the identifiability issue common in latent variable
models, here we also assume, without loss of generality, Cov(ft ) = Ir and that B T B is a
diagonal matrix, since rotating B will not affect the above decomposition (16).
We will need two major assumptions for our analysis: (1) the factors are pervasive in
the sense of Definition 3.1, and (2) there is a constant C > 0 such that kΣ−1
u k2 , kΣu k2 ≤ C,
which are standard assumptions in the factor model literature. The pervasive assumption is
reasonable in financial applications, since the factors have impacts on a large fraction of the
outcomes (Chamberlain and Rothschild, 1982; Bai, 2003). If the factor loadings {bi }di=1 are
regarded as random realizations from a bounded random vector, the assumption holds (Fan
et al., 2013).
Definition 3.1. In the factor model (15), the factors are called pervasive if there is a constant
C > 0 such that kBkmax ≤ C and the eigenvalues of the r by r matrix B T B/d are distinct
and bounded away from zero and infinity.
Let {λi , vi }ri=1 be the top r eigenvalues and eigenvectors of Σ, and similarly, {λi , v i }ri=1 for
BB T . In the following proposition, we show that pervasiveness is naturally connected to the
15

incoherence structure. This connects well between the econometrics and machine learning
literatures and provide a good interpretation on the concept of the incoherence. Its proof
can be found in the appendix.
Proposition 3.2. Suppose there exists a constant C > 0 such that kΣu k ≤ C. The factors
ft are pervasive if and only if the coherence µ(V ) for V = (v1 , . . . , vr ) ∈ Rd×r is bounded by
some constant, and λi = λi (Σ)  d for i ≤ r so that min1≤i6=j≤r |λi − λj |/λj > 0.
Our goal is to obtain a good covariance matrix estimator by exploiting the structure (16).
Our strategy is to use a generalization of the principal orthogonal complement thresholding
(POET) method proposed in Fan et al. (2013). The generic POET procedure encompasses
three steps:
b1 , . . . , λ
br ), Vb = (b
b Λ
b = diag(λ
(1) Given three pilot estimators Σ,
v1 , . . . , vbr ) respectively for
true covariance Σ, leading eigenvalues Λ = diag(λ1 , . . . , λr ) and leading eigenvectors
b u:
V = (v1 , . . . , vr ), compute the principal orthogonal complement Σ
bu = Σ
b − Vb Λ
b Vb T .
Σ

(17)

b u to obtain thresholded estimate Σ
b>
(2) Apply the correlation thresholding to Σ
u defined as
follows:
(
b u,ij ,
Σ
i=j
b>
Σ
,
(18)
u,ij =
b u,ij )I(|Σ
b u,ij | ≥ τij ), i 6= j
sij (Σ
where sij (·) is the generalized shrinkage function (Antoniadis and Fan, 2001; Rothman et al., 2009) and τij = τ (σ̂u,ii σ̂u,jj )1/2 is an entry-dependent threshold. τ will be
determined later in Theorem 3.1. This step exploits the sparsity of Σu .
b > = Vb Λ
b Vb T + Σ
b>
(3) Construct the final estimator Σ
u.
The key feature in the above procedure lies in the flexibility of choosing the pilot estimab according to data generating distribution. Typically
tors in the first step. We will choose Σ
b However, Λ̂ and V̂ in general
we can use λ̂i , v̂i for i ≤ r as the eigenvalues/vectors of Σ.
b and can be obtained separately via
do not have to come from the spectral information of Σ
different methods.
To guide the selection of proper pilot estimators, Fan et al. (2017+) provided a high
level sufficient condition for this simple procedure to be effective, and its performance is
P
gauged, in part, through the sparsity level of Σu , defined as md := maxi≤d j≤d |Σu,ij |q .
When q = 0, md corresponds to the maximum number of nonzero elements in each row of
Σu . For completeness, we present the theorem given by Fan et al. (2017+) in the following.
16

p
√
Theorem 3.1. Let wn =
log d/n + 1/ d. Suppose there exists C > 0 such that
b b b
kΣ−1
u k, kΣu k ≤ C and we have pilot estimators Σ, Λ, V satisfying
p
b − Σkmax = O( log d/n),
kΣ
p
bi /λi − 1| = O( log d/n),
|λ
√
kb
vi − vi k∞ = O(wn / d).

(19)
(20)
(21)

Under the pervasiveness condition of the factor model (15), with τ  wn , if md wn1−q = o(1),
the following rates of convergence hold with the generic POET procedure:

and



b > − Σu k2 = O md w1−q = k(Σ
b > )−1 − Σu −1 k2 ,
kΣ
u
n
u

(22)

 
b > − Σkmax = O wn ,
kΣ

 √d log d
>
b
+ md wn1−q ,
kΣ − ΣkΣ = O
n

>
−1
−1
1−q
b
k(Σ ) − Σ k2 = O md wn
,

(23)

where kAkΣ = d−1/2 kΣ−1/2 AΣ−1/2 kF is the relative Frobenius norm.
√
We remark that the additional term 1/ d in wn , is due to the estimation of unobservable
factors and is negligible when the dimensional d is high. The optimality of the above rates
of convergence is discussed in details in Fan et al. (2017+). Theorem 3.1 reveals a profound
deterministic connection between the estimation error bound of the pilot estimators with the
rate of convergences of the POET output estimators. Notice that the eigenvector estimation
error is under the `∞ norm, for which our `∞ perturbation bounds will prove to be useful.
b to be the
In this subsection, since we assume only bounded fourth moments, we choose Σ
robust estimate of covariance matrix Σ defined in (14). We now invoke our `∞ bounds to
show that the spectrum properties (eigenvalues and eigenvectors) are stable to perturbation.
b into a form such that Theorem 2.2 can be invoked:
Let us decompose Σ
b=
Σ

r
X

b − Σ),
λi v i v Ti + Σu + (Σ

i=1

b is viewed as A,
e the low-rank part Pr λi v i v Ti , which is also BB T , is viewed as A,
where Σ
i=1
and the remaining terms are treated as E. The following results follow immediately.
Proposition 3.3. Assume that there is a constant C > 0 such that kΣu k ≤ C. If the factors
17

bi , vbi
are pervasive, then with probability greater than 1 − d−1 , we have (19) – (21) hold with λ
b for i ≤ r. In addition, (22) and (23) hold.
as the leading eigenvalues/vectors of Σ
The inequality (19) follows directly from Proposition 3.1 under the assumption of
bounded fourth moments. It is also easily verifiable that (20), (21) follow from (19) by
√
Weyl’s inequality and Theorem 2.2 (noting that kΣu k∞ ≤ dkΣu k). See Section 3.2 for
more details.
Note that in the case of sub-Gaussian variables, sample covariance matrix and its leading
eigenvalues/vectors will also serve the same purpose due to (12) and Theorem 2.2 as discussed
in Section 3.1.
We have seen that the `∞ perturbation bounds are useful in robust covariance estimation,
and particularly, they resolve a theoretical difficulty in the generic POET procedure for factor
model based covariance matrix estimation.

4
4.1

Simulations
Simulation: the perturbation result

In this subsection, we implement numerical simulations to verify the perturbation bound
in Theorem 2.2. We will show that the error behaves in the same way as indicated by our
theoretical bound.
In the experiments, we let the matrix size d run from 200 to 2000 by an increment
of 200. We fix the rank of A to be 3 (r = 3). To generate an incoherence low rank
matrix, we sample a d × d random matrix with iid standard normal variables, perform
singular value decomposition, and extract the first r right singular vectors v1 , v2 , . . . , vr . Let
V = (v1 , . . . , vr ) and D = diag(rγ, (r − 1)γ, . . . , γ) where γ as before represents the eigengap.
Then, we set A = V DV T . By orthogonal invariance, vi is uniformly distributed on the unit
sphere Sd−1 . It is not hard to see that with probability 1 − O(d−1 ), the coherence of V
√
µ(V ) = O( log d).
We consider two types of sparse perturbation matrices E: (a) construct a d × d matrix
E0 by randomly selecting s entries for each row, and sampling a uniform number in [0, L] for
each entry, and then symmetrize the perturbation matrix by setting E = (E0 + E0T )/2; (b)
pick ρ ∈ (0, 1), L0 > 0, and let Eij = L0 ρ|i−j| . Note that in (b) we have kEk∞ ≤ 2L0 /(1 − ρ),
and thus we can choose suitable L0 and ρ to control the `∞ norm of E. This covariance
structure is common in cases where correlations between random variables depend on their
“distance” |i − j|, which usually arises from autoregressive models.

18

0.4

0.35
−1.1

0.3
−2

log (er r )

er r

0.25

0.2

−2.9

−3.8

0.15
−4.7
0.1
−5.6
0.05
−6.5
0
200

600

1000

1400

1800

5.5

6

6.5

7

7.5

log d

d

Figure 1: The left plot shows the perturbation error of eigenvectors against matrix size d
ranging from 200 to 2000, with different eigengap γ. The right plot shows log(err) against
log(d). The slope is around −0.5. Blue lines represent γ = 10; red lines γ = 50; green lines
γ = 100; and black lines γ = 500. We report the largest error over 100 runs.
The perturbation of eigenvectors is measured by the element-wise error:
err := max min kηi vei − vi k∞ ,
1≤i≤r ηi ∈{±1}

e = A + E in the descending order.
where {e
vi }ri=1 are the eigenvectors of A
To investigate how the error depends on γ and d, we generate E according to mechanism
(a) with s = 10, L = 3, and run simulations in different parameter configurations: (1) let
the matrix size d range from 200 to 2000, and choose the eigengap γ in {10, 50, 100, 500}
√
(Figure 1); (2) fix the product γ d to be one of {2000, 3000, 4000, 5000}, and let the matrix
size d run from 200 to 2000 (Figure 2).
To find how the errors behave for E generated from different methods, we run simulations
as in (1) but generate E differently. We construct E through mechanism (a) with L =
10, s = 3 and L = 0.6, s = 50, and also through mechanism (b) with L0 = 1.5, ρ = 0.9 and
L0 = 7.5, ρ = 0.5 (Figure 3). The parameters are chosen such that kEk∞ is about 30.
In Figure 1 – 3, we report the largest error based on 100 runs. Figure 1 shows that the
error decreases as d increases (the left plot); and moreover, the logarithm of the error is linear
√
in log(d), with a slope −0.5, that is, err ∝ 1/ d (the right plot). We can take the eigengap
19

0.03

75

65
0.025

55

√
er r × γ d

er r

0.02

0.015

0.01

45

35

25

15
0.005

5
0
200

600

1000

1400

1800

200

600

d

1000

1400

1800

d

Figure 2: The left plot shows the perturbation
error of eigenvectors against matrix size d
√
ranging from 200 to 2000, when
γ
d
is
kept
fixed,
with different values.
The right plot
√
√
shows
by γ d against d. Blue lines√represent γ d = 2000; red lines
√ the error multiplied √
γ d = 3000; green lines γ d = 4000; and black lines γ d = 5000. We report the largest
error over 100 runs.
γ into consideration and characterize the relationship in a more refined way. In Figure 2,
it is clear that err almost falls on the same horizontal line for different configurations of
√
√
d and γ, with γ d fixed. The right panel clearly indicates that err × γ d is a constant,
√
and therefore err ∝ 1/(γ d). In Figure 3, we find that the errors behave almost the same
regardless of how E is generated. These simulation results provide stark evidence supporting
the `∞ perturbation bound in Theorem 2.2.

4.2

Simulation: robust covariance esitmation

We consider the performance of the generic POET procedure in robust covariance estimation in this subsection. Note that the procedure is flexible in employing any pilot estimators
b Λ,
b Vb satisfying the conditions (19) – (21) respectively.
Σ,
We implemented the robust procedure with four different initial trios: (1) the sample
b S with its leading r eigenvalues and eigenvectors as Λ
b S and Vb S ; (2) the Huber’s
covariance Σ
b R given in (14) and its top r eigen-structure estimators Λ
b R and Vb R ; (3)
robust estimator Σ
b K with its corresponding Λ
b K and Vb K ; (4) lastly, we
the marginal Kendall’s tau estimator Σ
20

−1.1
−1

−2.2

−1.8

−3.2

−2.6

−4.2

−3.4

−5.2

−4.2
−6.3
−5

log (er r )

−7.3
−5.8
5.5

6

6.5

7

7.5

−1.2

5.5

6

6.5

7

7.5

5.5

6

6.5

7

7.5

−1.4

−2.2
−2.3
−3.1
−3.1
−4.1
−4
−5.1
−4.9
−6
−5.7
−7
−6.6
5.5

6

6.5

7

7.5

log d

Figure 3: These plots show log(err) aginst log(d), with matrix size d ranging from 200 to
2000 and different eigengap γ. The perturbation E is generated from different ways. Top
left: L = 10, s = 3; top right: L = 0.6, s = 50; bottom left: L0 = 1.5, ρ = 0.9; bottom right:
L0 = 7.5, ρ = 0.5. The slopes are around −0.5. Blue lines represent γ = 10; red lines γ = 50;
green lines γ = 100; and black lines γ = 500. We report the largest error over 100 runs.
use the spatial Kendall’s tau estimator to estimate the leading eigenvectors instead of the
marginal Kendall’ tau, so Vb K in (3) is replaced with Ve K . We need to briefly review the two
b K and Ve K .
types of Kendall’s tau estimators here, and specifically give the formula for Σ
Kendall’s tau correlation coefficient, for estimating pairwise comovement correlation, is
defined as
X
2
τ̂jk :=
sgn((ytj − yt0 j )(ytk − yt0 k )) .
(24)
n(n − 1) t<t0
Its population
 expectation is related to the Pearson correlation via the transform rjk =
sin π2 E[τ̂jk ] for elliptical distributions (which are far too restrictive for high-dimensional


π
applications). Then r̂jk = sin 2 τ̂jk is a valid estimation for the Pearson correlation rjk .
q
q
R
b
b
b
b R ) containing the robustly estimated stanLetting R = (r̂jk ) and D = diag( Σ11 , . . . , Σ
dd

21

dard deviations, we define the marginal Kendall’s tau estimator as
bK = D
bR
bD
b.
Σ

(25)

b we still use the robust variance estimates from Σ
b R.
In the above construction of D,
The spatial Kendall’s tau estimator is a second-order U-statstic, defined as
e K :=
Σ

X (yt − yt0 )(yt − yt0 )T
2
.
n(n − 1) t<t0
kyt − yt0 k22

(26)

e K . It has been shown by Fan et al.
Then Ve S is constructed by the top r eigenvectors of Σ
b K and its top r eigenvalues Λ
b K satisfy (19)
(2017+) that under elliptical distribution, Σ
and (20) while Ve S suffices to conclude (21). Hence Method (4) indeed provides good initial
b K attains (19) for elliptical
estimators if data are from elliptical distribution. However, since Σ
distribution, by similar argument for deriving Proposition 3.3 based on our `∞ pertubation
b K is also valid for the generic POET
bound, Vb K consisting of the leading eigenvectors of Σ
procedure. For more details about the two types of Kendall’s tau, we refer the readers to
Fang et al. (1990); Choi and Marden (1998); Han and Liu (2014); Fan et al. (2017+) and
references therein.
In summary, Method (1) is designed for the case of sub-Gaussian data; Method (3)
and (4) work under the situation of elliptical distribution; while Method (2) is proposed in
this paper for the general heavy-tailed case with bounded fourth moments without further
distributional shape constraints.
We simulated n samples of (ftT , uTt )T from two settings: (a) a multivariate t-distribution
with covariance matrix diag{Ir , 5Id } and various degrees of freedom (ν = 3 for very heavy
tail, ν = 5 for medium heavy tail and ν = ∞ for Gaussian tail), which is one example
of the elliptical distribution (Fang et al., 1990); (b) an element-wise iid one-dimensional t
distribution with the same covariance matrix and degrees of freedom ν = 3, 5 and ∞, which
is a non-elliptical heavy-tailed distribution.
Each row of coefficient matrix B is independently sampled from a standard normal distribution, so that with high probability, the pervasiveness condition holds with kBkmax =
√
O( log d). The data is then generated by yt = Bft + ut and the true population covariance
matrix is Σ = BB T + 5Id .
For d running from 200 to 900 and n = d/2, we calculated errors of the four robust
estimators in different norms. The tuning for α in minimization (14) is discussed more
p
throughly in Fan et al. (2017). For the thresholding parameter, we used τ = 2 log d/n.
b>
b > −1 − Σ−1 k and
The estimation errors are gauged in the following norms: kΣ
u − Σu k, k(Σ )
22

2.5
error ratio
1.0 1.5 2.0
0.5
0.0

0.0

0.5

error ratio
1.0 1.5 2.0

2.5

3.0

Spectral norm error of Sigma_u (IQR)

3.0

Spectral norm error of Sigma_u (Median)

200

400

600

800

200

400

600

800

Spectral norm error of inverse Sigma (Median)

Spectral norm error of inverse Sigma (IQR)
2.5
error ratio
1.0 1.5 2.0
0.5
0.0

0.0

0.5

error ratio
1.0 1.5 2.0

2.5

3.0

p

3.0

p

200

400

600

800

200

400

600

800

Relative Frobenius norm error of Sigma (Median)

Relative Frobenius norm error of Sigma (IQR)

error ratio
0.8
0.4
0.0

0.0

0.4

error ratio
0.8

1.2

p

1.2

p

200

400

600

800

200

p

400

600

800

p

Figure 4: Error ratios of robust estimates against varying dimension. Blue lines represent
errors of Method (2) over Method (1) under different norms; black lines errors of Method
(3) over Method (1); red lines errors of Method (4) over Method (1). (ftT , uTt ) is generated
by multivariate t-distribution with df = 3 (solid), 5 (dashed) and ∞ (dotted). The median
errors and their IQR’s (interquartile range) over 100 simulations are reported.
b > − ΣkΣ as shown in Theorem 3.1. The two different settings are separately plotted in
kΣ
b S in Method
Figures 4 and 5. The estimation errors of applying sample covariance matrix Σ
(1) are used as the baseline for comparison. For example, if relative Frobenius norm is used to
b > )(k) − ΣkΣ /k(Σ
b > )(1) − ΣkΣ will be depicted for k = 2, 3, 4, where
measure performance, k(Σ
b > )(k) are generic POET estimators based on Method (k). Therefore if the ratio curve
(Σ
moves below 1, the method is better than naive sample estimator (Fan et al., 2013) and
23

vice versa. The more it gets below 1, the more robust the procedure is against heavy-tailed
randomness.
The first setting (Figure 4) represents a heavy-tailed elliptical distribution, where we
expect Methods (2), (3), (4) all outperform the POET estimator based on the sample covariance, i.e. Method (1), especially in the presence of extremely heavy tails (solid lines for
ν = 3). As expected, all three curves under various measures show error ratios visibly smaller
than 1. On the other hand, if data are indeed Gaussian (dotted line for ν = ∞), Method (1)
has better behavior under most measures (error ratios are greater than 1). Nevertheless, our
robust Method (2) still performs comparably well with Method (1), whereas the median error
ratios for the two Kendall’s tau methods are much worse. In addition, the IQR (interquartile
range) plots reveal that Method (2) is indeed more stable than two Kendall’s tau Methods
(3) and (4). It is also noteworthy that Method (4), which leverages the advantage of spatial
Kendall’s tau, performs more robustly than Method (3), which solely base its estimation of
the eigen-structure on marginal Kendall’s tau.
The second setting (Figure 5) provides an example of non-elliptical distributed data. We
can see that the performance of the general robust Method (2) dominates the other three
methods, which verifies the benefit of robust estimation for a general heavy-tailed distribution. Note that Kendall’s tau methods do not apply to distributions outside the elliptical
family, excluding even the element-wise iid t distribution in this setting. Nonetheless, even
in the first setting where the data are indeed elliptical, with proper tuning, the proposed
robust method can still outperform Kendall’s tau by a clear margin.

5
5.1

Proof Organization of Main Theorems
Symmetric Case

√
For shorthand, we write τ = kEk∞ , and κ = d kEV kmax . An obvious bound for κ is
√
κ ≤ rµ τ (by Cauchy-Schwarz inequality). We will use these notations throughout this
subsection.
Recall the spectral decomposition of A in (8). Expressing E in terms of the column
vectors of V and V⊥ , which form an orthogonal basis in Rn , we write
[V, V⊥ ]T E[V, V⊥ ] =:

E11 E12
E21 E22

!
.

(27)

T
Note that E12 = E21
since E is symmetric. Conceptually, the perturbation results in a

24

2.5
error ratio
1.0 1.5 2.0
0.5
0.0

0.0

0.5

error ratio
1.0 1.5 2.0

2.5

3.0

Spectral norm error of Sigma_u (IQR)

3.0

Spectral norm error of Sigma_u (Median)

200

400

600

800

200

400

600

800

Spectral norm error of inverse Sigma (Median)

Spectral norm error of inverse Sigma (IQR)
2.5
error ratio
1.0 1.5 2.0
0.5
0.0

0.0

0.5

error ratio
1.0 1.5 2.0

2.5

3.0

p

3.0

p

200

400

600

800

200

400

600

800

Relative Frobenius norm error of Sigma (Median)

Relative Frobenius norm error of Sigma (IQR)

error ratio
0.8
0.4
0.0

0.0

0.4

error ratio
0.8

1.2

p

1.2

p

200

400

600

800

200

p

400

600

800

p

Figure 5: Error ratios of robust estimates against varying dimension. Blue lines represent
errors of Method (2) over Method (1) under different norms; black lines errors of Method (3)
over Method (1); red lines errors of Method (4) over Method (1). (ftT , uTt ) is generated by
element-wise iid t-distribution with df = 3 (solid), 5 (dashed) and ∞ (dotted). The median
errors and their IQR’s (interquartile range) over 100 simulations are reported.
rotation of [V, V⊥ ], and we write a candidate orthogonal basis as follows:
V := (V + V⊥ Q)(Ir + QT Q)−1/2 ,

V ⊥ := (V⊥ − V QT )(Id−r + QQT )−1/2 ,

(28)

where Q ∈ R(d−r)×r is to be determined. It is straightforward to check that [V , V ⊥ ] is
e , V ⊥ ) is a block
an orthogonal matrix. We will choose Q in a way such that (V , V ⊥ )T A(V
25

T
e = 0. Substituting (28) and simplifying the equation, we obtain
diagonal matrix, i.e., V ⊥ AV

Q(Λ1 + E11 ) − (Λ2 + E22 )Q = E21 − QE12 Q.

(29)

The approach of studying perturbation through a quadratic equation is known (see Stewart
(1990) for example). Yet, to the best of our knowledge, existing results study perturbation
under orthogonal-invariant norms (or unitary-invariant norms in the complex case), which
includes a family of matrix operator norms and Frobenius norm, but excludes the matrix
max-norm. The advantages of orthogonal-invariant norms are pronounced: such norms of a
symmetric matrix only depend on its eigenvalues regardless of eigenvectors; moreover, with
suitable normalization they are consistent in the sense kABk ≤ kAk · kBk. See Stewart
(1990) for a clear exposition.
The max-norm, however, does not possess these important properties. An imminent issue
is that it is not clear how to relate Q to V⊥ Q, which will appear in (29) after expanding E
according to (27), and which we want to control. Our approach here is to study Q := V⊥ Q
directly through a transformed quadratic equation, obtained by left multiplying V⊥ to (29).
Denote H = V⊥ E21 , Q = V⊥ Q, L1 = Λ1 + E11 , L2 = V⊥ (Λ2 + E22 )V⊥T . If we can find an
appropriate matrix Q with Q = V⊥ Q, and it satisfies the quadratic equation
Q L1 − L2 Q = H − QH T Q,

(30)

then Q also satisfies the quadratic equation (29). This is because multiplying both sides of
(30) by V⊥T yields (29), and thus any solution Q to (30) with the form Q = V⊥ Q must result
in a solution Q to (29).
e , V ⊥ ) is a block diagonal
Once we have such Q (or equivalently Q), then (V , V ⊥ )T A(V
matrix, and the span of column vectors of V is a candidate space of the span of first r
eigenvectors, namely span{e
v1 , . . . , ver }. We will verify the two spaces are identical in Lemma
5.3. Before stating that lemma, we first provide bounds on kQkmax and kV − V kmax .
Lemma 5.1. Suppose |λr | − ε > 4rµ(τ + 2rκ). Then, there exists a matrix Q ∈ R(d−r)×r
such that Q = V⊥ Q ∈ Rd×r is a solution to the quadratic equation (30), and Q satisfies
√
kQkmax ≤ ω/ d. Moreover, if rω < 1/2, the matrix V defined in (28) satisfies
√
√
kV − V kmax ≤ 2 µ ωr/ d .

(31)

Here, ω is defined as ω = 8(1 + rµ)κ/(|λr | − ε).
The second claim of the lemma (i.e., the bound (31)) is relatively easy to prove once
the first claim (i.e., the bound on kQkmax ) is proved. To understand this, note that we can
26

T

T

rewrite V as V = (V + Q)(Ir + Q Q)−1/2 , and kQ Qkmax can be controlled by a trivial
T
inequality kQ Qkmax ≤ dkQk2max ≤ w2 . To prove the first claim, we construct a sequence of
matrices through recursion that converges to the fixed point Q, which is a solution to the
quadratic equation (30). For all iterates of matrices, we prove a uniform max-norm bound,
0
which leads to a max-bound on kQkmax by continuity. To be specific, we initialize Q = 0,
t
and given Q , we solve a linear equation:
t

t

Q L1 − L2 Q = H − Q H T Q ,
t+1

(32)
t

and the solution is defined as Q . Under some conditions, the iterate Q converges to a
limit Q, which is a solution to (30). The next general lemma captures this idea. It follows
from Stewart (1990) with minor adaptations.
Lemma 5.2. Let T be a bounded linear operator on a Banach space B equipped with a norm
k · k. Assume that T has a bounded inverse, and define β = kT −1 k−1 . Let ϕ : B → B be a
map that satisfies
kϕ(x)k ≤ ηkxk2 ,

and

kϕ(x) − ϕ(y)k ≤ 2η max{kxk, kyk}kx − yk

(33)

for some η ≥ 0. Suppose that B0 is a closed subspace of B such that T −1 (B0 ) ⊆ B0 and
ϕ(B0 ) ⊆ B0 . Suppose y ∈ B0 that satisfies 4ηkyk < β 2 . Then, the sequence initialized with
x0 = 0 and iterated through
xk+1 = T −1 (y + ϕ(xk )),

k≥0

(34)

converges to a solution x? to T x = y+ϕ(x). Moreover, we have x? ⊆ B0 , and kx? k ≤ 2kyk/β.
To apply this lemma to the equation (30), we view B as the space of matrices Rd×r
with the max-norm k · kmax , and B0 as the subspace of matrices of the form V⊥ Q where
Q ∈ R(d−r)×r . The linear operator T is set to be the T (Q) = Q L1 − L2 Q, and the map ϕ is
set to be the quadratic function ϕ(Q) = −QH T Q. Roughly speaking, under the assumption
of Lemma 5.2, the nonlinear effect caused by ϕ is weak compared with the linear operator
T . Therefore, it is crucial to show T is invertible, i.e. to give a good lower bound on
kT −1 k−1
max = inf kQkmax =1 kT (Q)kmax . Since the norm is not orthogonal-invariant, a subtle
issue arises when A is not of exact low rank, which will be discussed at the end of the
subsection.
t
If there is no perturbation (i.e., E = 0), all the iterates Q are simply 0, so V is identical
to V . If the perturbation is not too large, the next lemma shows that the column vectors of
V span the same space as span{e
v1 , . . . , ver }.
27

In other words, with a suitable orthogonal matrix R, the columns of V R are ve1 , . . . , ver .
Lemma 5.3. Suppose |λr | − ε > max{3τ, 64(1 + rµ)r3/2 µ1/2 κ}. Then, there exists an orthogonal matrix R ∈ Rr×r such that the column vectors of V R are ve1 , . . . , ver .
Proof of Theorem 2.1. It is easy to check that under the assumption of Theorem 2.1,
the conditions required in Lemma 5.1 and Lemma 5.3 are satisfied. Hence, the two lemmas
imply Theorem 2.1.
To study the perturbation of individual eigenvectors, we assume, in addition to the condition on |λr |, that λ1 , . . . , λr satisfy a uniform gap, (namely δ > kEk2 ). This additional
assumption is necessary, because otherwise, the perturbation may lead to a change of relative order of eigenvalues, and we may be unable to match eigenvectors from the order of
e
eigenvalues. Suppose R ∈ Rr×r is an orthogonal matrix such that V R are eigenvectors of A.
Now, under the assumption of of Theorem 2.1, the column vectors of Ve and V R are identical
up to sign, so we can rewrite the difference Ve − V as
Ve − V = V (R − Ir ) + (V − V ).

(35)

We already provided a bound on kV − V kmax in Lemma 5.1. By the triangular inequality,
we can derive a bound on kV kmax . If we can prove a bound on kR − Ir kmax , it will finally
leads to a bound on kVe − V kmax . In order to do so, we use the Davis-Kahan theorem to
obtain an bound on he
vi , vi i for all i ∈ [r]. This will lead to a max-norm bound on R − Ir
(with the price of potentially increasing the bound by a factor of r). The details about the
proof of Theorem 2.2 are in the appendix.
We remark that, we assume conditions on |λr | −  in Theorem 2.1 and Theorem 2.2,
which are only useful in cases where |λr | > kA − Ar k∞ . Ideally, we would like to have results
with assumptions only involving λr and λr+1 , since Davis-Kahan theorem only requires a
gap in neighboring eigenvalues. Unfortunately, unlike orthogonal-invariant norms that only
depend on the eigenvalues of a matrix, the max-norm k · kmax is not orthogonal-invariant,
and thus it also depends on the eigenvectors of a matrix. For this reason, it is not clear
whether we could obtain a lower bound on kT −1 k−1
max using only the eigenvalues λr and λr+1
so that we could apply Lemma 5.2. The analysis appears to be difficult if we do not have a
bound on kT −1 k−1
max , considering that even in the analysis of linear equations, we also need
invertibility, condition numbers, etc.

28

5.2

Asymmetric Case

Let Ad , E d be d1 + d2 square matrices defined as
d

A :=

0 A
AT 0

!

0 E
ET 0

d

,

E :=

!
.

ed := Ad + E d . This augmentation of an asymmetric matrix into a symmetric
Also denote A
one is called Hermitian dilation. Here the superscript d means the Hermitian dilation. We
ed .
also use this notation to denote quantities corresponding to Ad and A
An important observation is that
0 A
AT 0

!

ui
± vi

!

ui
± vi

= ± σi

!
.

From this identity, we know that Ad have nonzero eigenvalues ± σi where 1 ≤ i ≤ rank(A),
and its corresponding eigenvectors are (uTi , ± viT )T . For a given r, we stack these (normalized)
eigenvectors with indices i ∈ [r] into a matrix V d ∈ R(d1 +d2 )×2r :
1
V d := √
2

U U
V −V

!
.

Through the augmented matrices, we can transfer eigenvector results for symmetric matrices
to singular vectors of asymmetric matrices. However, we cannot directly invoke the results
proved for symmetric matrices, due to an issue about the coherence of V d : when d1 and d2
are not comparable, the coherence µ(V d ) can be very large even when µ(V ) and µ(U ) are
bounded. To understand this, consider the case where r = 1, d1  d2 , and all entries of U
√
√
are O(1/ d1 ), and all entries of V are O(1/ d2 ). Then, the coherences µ(U ) and µ(V ) are
O(1), but µ(V d ) = O((d1 + d2 )/d2 )  1.
This unpleasant issue about the coherence, nevertheless, can be tackled if we consider a
different matrix norm. In order to deal with the different scales of d1 and d2 , we define the
weighted max-norm for any matrix M with d1 + d2 rows as follows:
√
kM kw :=

d1 Id1
0
√
0
d2 Id2

!
M

.

(36)

max

√
In other words, we rescale the top d1 rows of M by a factor of d1 , and rescale the bottom
√
d2 rows by d2 . This weighted norm serves to balance the potential different scales of d1
and d2 .
29

The proofs of theorems in Section 2.2 will be almost the same with those in the symmetric
case, with the major difference being the new matrix norm. Because the derivation is slightly
repetitive, we will provide concise proofs in the appendix . Similar to the decomposition in
(2.1),
!
!
0
A
−
A
0
A
r
r
=: Adr + (Ad − Adr ),
+
Ad =
T
T
T
0
A − Ar
Ar 0
where Adr is has rank 2r. Equivalently,
Adr

=

r
X

σi (uTi , viT )T (uTi , viT )

−

r
X

σi (uTi , −viT )T (uTi , −viT ).

i=1

i=1

Analogously, we define notations in (28)–(30) and use d in the superscript to signify
that they are augmented through Hermitian dilation. It is worthwhile to note that
Λd1 = diag{σ1 , . . . , σr , −σr , . . . , −σ1 }, and that min{| ± σi | : i ∈ [r]} = σr (a similar
p
p
quantity as |λr |). Recall µ0 = µ(U ) ∨ µ(V ), τ0 =
d1 /d2 kEk∞ ∨ d2 /d1 kEk1 and
p
p
d1 /d2 kA − Ar k∞ ∨ d2 /d1 kA − Ar k1 . In the proof, we will also use κ0 =
ε0 =
√
√
max{ d1 kEV kmax , d2 kE T U kmax }, which is a quantity similar to κ.
The next key lemma, which is parallel to Lemma 5.1, provides a bound on the solution
d
Q to the quadratic equation
d

d

d

d

d

d

Q L1 − L2 Q = H d − Q (H d )T Q .

(37)

Lemma 5.4. Suppose σr − ε0 > 16rµ0 (τ0 + rκ0 ). Then, there exists a matrix Qd ∈
d
R(d1 +d2 −2r)×2r such that Q = V⊥d Qd ∈ R(d1 +d2 )×2r is a solution to the quadratic equation
d
d
d
(37), and Q satisfies kQ kw ≤ ω0 . Moreover, if rω0 < 1/2, the matrix V defined in (28)
satisfies
√
d
kV − V d kw ≤ 6 µ0 rω0 .
(38)
Here, ω0 is defined as ω0 = 8(1 + rµ0 )κ0 /3(σr − ε0 ).
In this lemma, the bound (38) bears a similar form to (31): if we consider the max-norm,
d
the first d1 rows of V − V d correspond to the left singular vectors ui ’s, and they scale with
√
1/ d1 ; and the last d2 rows correspond to the right singular vectors vi ’s, which scale with
√
1/ d2 . Clearly, the weighted max-norm k · kw indeed helps to balance the two dimensions.

30

A

Proofs for Section 2.1

Denote the column span of a matrix M by span(M ). Suppose two matrices M1 , M2 ∈
R
(m ≤ n) have orthonormal column vectors. It is known that (Stewart, 1990)
n×m

d(M1 , M2 ) := kM1 M1T − M2 M2T k2 = k sin Θ(M1 , M2 )k2 .

(39)

where Θ(M1 , M2 ) are the canonical angles between span(M1 ) and span(M2 ). Recall the
√
notations defined in (27), and also recall κ = dkEV kmax , Λ1 = diag{λ1 , . . . , λr }, Λ2 =
diag{λr+1 , . . . , λn }, L1 = Λ1 + E11 , L2 = V⊥ (Λ2 + E22 )V⊥T and H = V⊥ E21 . The first lemma
bounds kHkmax .
Lemma A.1. We have the following bound on kHkmax :
√
kHkmax ≤ (1 + rµ)κ/ d.
Proof of Lemma A.1. Using the definition E21 = V⊥T EV in (27), we can write H =
V⊥ V⊥T EV . Since the columns of V and V⊥ form an orthogonal basis in Rd , clearly
V V T + V⊥ V⊥T = Id .

(40)

By Cauchy-Schwarz inequality and the definition of µ, for any i, j ∈ [d],
T

|(V V )ij | =

r
X
k=1

|Vik Vjk | ≤

r
X

1/2
Vik2

k=1

·

r
X
k=1

Vjk2

1/2

≤

rµ
.
d

Using the identity (40) and the above inequality, we derive
kHkmax ≤ kEV kmax + kV V T EV kmax
≤ (1 + dkV V T kmax ) kEV kmax ≤ (1 + rµ)kEV kmax ,
which completes the proof.
√
Lemma A.2. If |λr | > κr µ, then L1 is an invertible matrix. Furthermore,
inf

kQ0 kmax =1

kQ0 L1 − L2 Q0 kmax ≥ |λr | − 3rµ(τ + rκ) − ε ,

where Q0 is an d × r matrix.

31

(41)

Proof of Lemma A.2. Let Q0 be any d × r matrix with kQ0 kmax = 1. Note
Q0 L1 − L2 Q0 = Q0 Λ1 + Q0 E11 − L2 Q0 .
We will derive upper bounds on Q0 E11 and L2 Q0 , and a lower bound on Q0 Λ1 . Since
E11 = V T EV by definition, we expand Q0 E11 and use a trivial inequality to derive
kQ0 E11 kmax ≤ d kQ0 V T kmax kEV kmax .

(42)

By Cauchy-Schwarz inequality and the definition of µ in (3), for i, j ∈ [d],
T

|(Q0 V )ij | ≤

r
X
k=1

|(Q0 )ik Vjk | ≤

r
X

1/2
(Q0 )2ik

k=1

r
X

1/2
Vjk2

≤

√

r
r·

k=1

rµ
,
d

√
Substituting kEV kmax = κ/ d into (42), we obtain an upper bound:
√
kQ0 E11 kmax ≤ κr µ .

(43)

To bound L2 Q0 = (V⊥ E22 V⊥T + (A − Ar ))Q0 , we use the identity (40) and write
V⊥ E22 V⊥T Q0 = V⊥ V⊥T EV⊥ V⊥T Q0 = (Id − V V T )E(Id − V V T )Q0 .
Using two trivial inequalities kEQ0 kmax ≤ kEk∞ kQ0 kmax = kEk∞ and kV T Q0 kmax ≤
√
kV T k∞ kQ0 kmax ≤ d, we have
kE(Id − V V T )Q0 kmax ≤ kEQ0 kmax + rkEV kmax kV T Q0 kmax
√
≤ kEk∞ + r d kEV kmax = τ + rκ .
In the proof of Lemma A.1, we showed kV V T kmax ≤ rµ/d. Thus,
kV⊥ E22 V⊥T Q0 kmax ≤ (1 + d kV V T kmax ) · kE(Id − V V T )Q0 kmax ≤ (1 + rµ)(τ + rκ) .
Moreover, k(A − Ar )Q0 kmax ≤ kA − Ar k∞ kQ0 kmax = ε. Combining the two bounds,
kL2 Q0 kmax ≤ (1 + rµ)(τ + rκ) + ε.

(44)

It is straightforward to obtain a lower bound on kQ0 Λ1 kmax : since there is an entry of Q0 ,

32

say (Q0 )ij , that has an absolute value of 1, we have
kQ0 Λ1 kmax ≥ |(Q0 )ij λj | ≥ |λr |.

(45)

To show L1 is invertible, we use (42) and (45) to obtain
√
kQ0 L1 kmax ≥ kQ0 Λ1 kmax − kQ0 E11 kmax ≥ |λr | − κr µ .
√
When |λr | − κr µ > 0, L1 must have full rank, because otherwise we can choose an approT
priate Q0 in the null space of L1 so that Q0 L1 = 0, which is a contradiction. To prove the
second claim of the lemma, we combine the lower bound (45) with upper bounds (43) and
(44) to derive
kQ0 L1 − L2 Q0 kmax ≥ kQ0 L1 kmax − kQ0 E11 kmax − kL2 Q0 kmax
√
≥ |λr | − κr µ − (1 + rµ)(τ + rκ) − ε
≥ |λr | − 3rµ(τ + rκ) − ε ,
which is exactly the desired inequality.
Next we prove Lemma 5.2. This lemma follows from Stewart (1990), with minor changes
that involves B0 . We provide a proof for the sake of completeness.
Proof of Lemma 5.2. Let us write α = kyk for shorthand and recall β = kT −1 k−1 . As
the first step, we show that the sequence {xk }∞
k=0 is bounded. By construction in (34), we
bound kxk+1 k using kxk k:
kxk+1 k ≤ kT −1 k(kyk + kϕ(xk )k) ≤

α η
+ kxk k2 .
β β

We use this inequality to derive an upper bound on {xk } for all k. We define ξ0 = 0 and
ξk+1 =

α η 2
+ ξ ,
β β k

k ≥ 0,

then clearly kxk k ≤ ξk (which can be shown by induction). It is easy to check (by induction)
2
that the sequence {ξk }∞
k=1 is increasing. Moreover, since 4αη < β , the quadratic function
φ(ξ) =

α η 2
+ ξ ,
β β

33

has two fixed points (namely solutions to φ(ξ) = ξ), and the smaller one satisfies
ξ? =

2α
2α
p
.
<
β
β + β 2 − 4ηα

If ξk < ξ? , then ξk+1 = φ(ξk ) ≤ φ(ξ? ) = ξ? . Thus, by induction, all ξk are bounded by ξ? .
This implies kxk k ≤ ξ? < 2α/β. The next step is to show that the sequence {xk } converges.
Using the recursive definition (34) again, we derive
kxk+1 − xk k ≤ kT −1 kkϕ(xk ) − ϕ(xk−1 )k
≤ 2β −1 η max{kxk k, kxk−1 k}kxk − xk−1 k
4αη
≤ 2 kxk − xk−1 k.
β
Since 4αη/β 2 < 1, the sequence {xk }∞
k=0 is a Cauchy sequence, and convergence is secured.
Let x? ∈ B be the limit. It is clear by assumption that xk ∈ B0 implies xk+1 ∈ B0 , so x? ∈ B0
and kx? k ≤ 2α/β by continuity.
The final step is to show x? is a solution to T x = y+φ(x). Because {xk }∞
k=0 is bounded and
∞
?
φ satisfies (33), the sequence {φ(xk )}k=0 converges to φ(x ) by continuity and compactness.
The linear operator T is also continuous, so we can take limits on both sides of T xk+1 =
y + φ(xk ), we conclude that x? is a solution to T x = y + φ(x).
With all the preparations, we are now ready to present the key lemma. As discussed in
Section 5, we set
B0 := {Q ∈ Rd×r : Q = V⊥ Q for some Q ∈ R(d−r)×r }.
which is a subspace of B = Rd×r . Consider the matrix max-norm k · kmax in B.
Lemma A.3. Suppose |λr | − ε > 4rµ(τ + 2rκ). Then there exists a solution Q ∈ B0 to the
equation (30) with
8(1 + rµ)κ
kQkmax ≤
√ .
|λr | − ε d
Proof of Lemma A.3. We will invoke Lemma 5.2 and apply it to the quadratic equation (30).
To do so, we first check the conditions required in Lemma 5.2.
Let the linear operator T be T Q = QL1 − L2 Q. By Lemma A.2, T has a bounded
inverse, and β := kT −1 k−1
max is bounded from below:
β ≥ |λr | − 3rµ(τ + rκ) − ε .
34

(46)

Let us define ϕ by ϕ(Q) = QH T Q. To check the inequalities in (33), observe that
√
kϕ(Q)kmax ≤ rdkQkmax kHkmax kQkmax ≤ (1 + rµ)κr d kQk2max
where we used Lemma A.1. We also observe
kϕ(Q1 ) − ϕ(Q2 )kmax = kQ1 H T (Q1 − Q2 ) + (Q1 − Q2 )H T Q2 kmax
≤ rdkQ1 kmax kHkmax kQ1 − Q2 kmax + rdkQ1 − Q2 kmax kHkmax kQ2 kmax
√
≤ 2(1 + rµ)κr d max{kQ1 kmax , kQ2 kmax }kQ1 − Q2 kmax .
√
Thus, if we set η = (1 + rµ)κr d, then inequalities required in (33) are satisfied. For any
Q with Q = V⊥ Q ∈ B0 , obviously ϕ(Q) = V⊥ QH T Q ∈ B0 . To show T −1 (Q) ∈ B0 , let
Q0 = T −1 (Q) and observe that
Q0 L1 − L2 Q0 = Q ∈ B0 .
By definition, we know L2 Q0 = V⊥ (E22 + Λ2 )V⊥T Q0 ∈ B0 , so we deduce Q0 L1 ∈ B0 . Our
√
assumption implies |λr | > κr µ, so by Lemma A.2, the matrix L1 is invertible, and thus
Q0 ∈ B0 . The last condition we check is 4ηkHkmax < β 2 . By Lemma A.1 and (46), this is
true if

2
4(1 + rµ)2 κ2 r < |λr | − 3rµ(τ + rκ) − ε .
The above inequality holds when |λr | > 4rµ(V )(τ + 2rκ) + ε. Under this condition, we have,
by Lemma 5.2,
kQkmax ≤

2(1 + rµ)κ
8(1 + rµ)κ
√ ,
√ ≤
|λr | − 3rµ(τ + rκ) − ε d
(|λr | − ε) d

where, the second inequality is due to 3rµ(τ + rκ) ≤ 3(|λr | − ε)/4.
The next lemma is a consequence of Lemma A.3. We define, as in Lemma 5.1, that
ω = 8(1 + rµ)κ/(|λr | − ε).
Lemma A.4. If rω 2 < 1/2, then
T

k(Ir + Q Q)−1/2 − Ir kmax ≤ rω 2 ,

3
T
k(Ir + Q Q)−1/2 kmax ≤ .
2

Proof of Lemma A.4. By the triangular inequality, the second inequality is immediate from
T
the first one. To prove the first inequality, suppose the spectral decomposition of Q Q is
T
T
Q Q = U Σ U , where Σ = diag{λ1 , . . . , λr } where λ1 ≥ . . . ≥ λr , and U = [u1 , u2 , . . . , ur ]
35

T

where u1 , . . . , ur are orthonormal vectors in Rr . Since Q Q has nonnegative eigenvalues, we
have λr ≥ 0. Using these notations, we can rewrite the matrix as
T

−1/2

(Ir + Q Q)

− Ir =

r
X


(1 + λi )−1/2 − 1 ui uTi .

i=1
T

Note that λ1 ≤ kQ Qk ≤ rdkQk2max ≤ rω 2 , which implies λ1 < 1/2. It is easy to check that
1+|x| ≥ (1+x)−1/2 ≥ 1−|x| whenever |x| < 1/2. From this fact, we know |(1+λi )−1/2 −1| ≤
λi ≤ rω 2 . Using Cauchy-Schwarz inequality, we deduce that for any j, k ∈ [d],
T

[(Ir + Q Q)

−1/2

− Ir ]jk ≤

r
X

(1 + λi )−1/2 − 1 · U ji U ki

i=1
2

≤ rω ·

r
X
i=1

2 1/2
U ji

r
X

2

U ki

1/2

i=1

2

≤ rω .
This leads to the desired max-norm bound.
Proof of Lemma 5.1. The first claim of the lemma (the existence of Q and its max-norm
bound) follows directly from Lemma A.3. To prove the second claim, we split V − V into
two parts:

V − V = V (Ir + QT Q)−1/2 − Ir + V⊥ Q(Ir + QT Q)−1/2


T
T
−1/2
− Ir + Q(Ir + Q Q)−1/2 ,
= V (Ir + Q Q)

(47)

where we used identity V⊥T V⊥ = Id−r . Note that rω < 1/2 implies rω 2 = (rω)2 /r < 1/(4r) <
1/2. Thus, we can use Lemma A.4 and derive



T

−1/2

V (Ir + Q Q)

− Ir

r



≤
max

r2 µ
T
k(Ir + Q Q)−1/2 − Ir kmax ≤
d

r

µ 2 2
r ω .
d

where we used Cauchy-Schwarz inequality. Using the above inequality and the bound on
kQkmax (namely, the first claim in the lemma),
r

µ 2 2
T
r ω + rkQkmax k(Ir + Q Q)−1/2 kmax
d
√
√ 2 2
≤ ( µ ω r + 3ωr/2)/ d.

kV − V kmax ≤

36

Simplifying the bound using rω ≤ 1/2 and a trivial bound µ ≥ 1, we obtain (31).
Proof of Lemma 5.3. Using the identity in (39), it follows from Davis-Kahan sin Θ theorem (Davis and Kahan, 1970) and Weyl’s inequality that
d(Ve , V ) ≤

kEk2
,
δr − kEk2

when δr > kEk2 , where δr = |λr | − |λr+1 |. Since λr+1 ≤ kA − Ar k2 ≤ ε and kEk2 ≤ τ , the
condition |λr | − ε > 3τ implies δr > 3kEk2 . Hence, we have d(Ve , V ) < 1/2. Moreover,
T

d(V , V ) = kV V − V V T k2 ≤ kV (V − V )T k2 + k(V − V )V T k2
√
≤ 2kV − V k2 ≤ 2 rd kV − V kmax
√
≤ 4r3/2 µ ω,
√
where we used a trivial inequality kM k2 ≤ kM kF ≤ rd kM kmax for any M ∈ Rd×r . Under
√
the condition |λr | −  > 64(1 + rµ)r3/2 µ1/2 κ, it is easy to check that 4r3/2 µ ω ≤ 1/2. Thus,
we obtain d(V , V ) < 1/2. By the triangular inequality,
d(Ve , V ) ≤ d(Ve , V ) + d(V , V ) < 1.
e , V ⊥ ) is a block diagonal matrix, span(V ) is the same as the subspace
Since (V , V ⊥ )T A(V
e We claim that span(V ) = span(e
spanned by r eigenvectors of A.
v1 , . . . , ver ). Otherwise, there
e1 , . . . , λ
er
exists an eigenvector u ∈ span(V ) whose associated eigenvalue is distinct from λ
(since δr > 3kEk2 ), and thus u is orthogonal to ve1 , . . . , ver . Therefore,
T

T

k(Ve Ve T − V V )uk2 = kV V uk2 = kuk2 .
This implies d(Ve , V ) ≥ 1, which is a contradiction.
Proof of Theorem 2.2. We split Ve − V into two parts—see (35). In the following, we first
obtain a bound on kR − Ir kmax , which then results in a bound on kVe − V kmax .
Under the assumption of the theorem, rω < 1/2, so
p
√
√
√
√
kV kmax ≤ kV − V kmax + kV kmax ≤ (2 µ rω)/ d + rµ/ d ≤ 2 rµ/d .

(48)

T
T
To bound kR − Ir kmax , we rewrite R as R = V V R = V Ve . Expand V according to (28),
T

R = (Ir + Q Q)−1/2 (V + Q)T Ve .
37

√
T
Let us make a few observations: (a) kQ Ve kmax ≤ dkQkmax ≤ ω by Cauchy-Schwarz
T
inequality; (b) kVe T V kmax ≤ 1 by Cauchy-Schwarz inequality again; and (c) k(Ir +Q Q)−1/2 −
Ir kmax ≤ rω 2 by Lemma A.4. Using these inequalities, we have
T

kR − (V + Q)T Ve kmax ≤ rk(Ir + Q Q)−1/2 − Ir kmax k(V + Q)T Ve kmax
≤ r2 ω 2 (1 + ω) .

(49)

Furthermore, by Davis-Kahn sin Θ theorem (Davis and Kahan, 1970) and Weyl’s inequality,
for any i ∈ [r],
p
kEk2
.
(50)
sin θ(vi , vei ) = 1 − hvi , vei i2 ≤
δ − kEk2
when δ > kEk2 (δ is defined in Theorem 2.2). This leads to the bound sin θ(vi , vei ) ≤ 2kEk2 /δ
(which is a simplified bound). This is because when δ ≥ 2kEk2 , the bound is implied by
(50); when δ < 2kEk2 , the bound trivially follows from sin θ(vi , vei ) ≤ 1. We obtain, up to
sign, for i ≤ r,
p
p
2kEk2
1 − hvi , vei i ≤ 1 − hvi , vei i2 ≤
.
(51)
δ
In other words, each diagonal entry of Ir −V T Ve , namely 1−hvi , vei i, is bounded by 4kEk22 /δ 2 .
P
Since {e
vi }ri=1 are orthonormal vectors, we have 1 − hvi , vei i2 ≥ i0 6=i hvi , vei0 i2 ≥ hvi , vej i2 for
any i 6= j, which leads to bounds on off-diagonal entries of V T Ve − Ir . We will combine the
two bounds. Note that when δ ≥ 2kEk2 ,
kV T Ve − Ir kmax ≤ max

 4kEk22 2kEk2
,
δ2
δ

=

2kEk2
;
δ

and when δ < 2kEk2 , kV T Ve − Ir kmax is trivially bounded by 1 (up to sign), which is trivially
bounded by 2kEk2 /δ. In either case, we deduce
kV T Ve − Ir kmax ≤

2kEk2
.
δ

T
Using the bounds in (49) and (52) and kQ Ve kmax ≤ ω, we obtain
T
kR − Ir kmax ≤ kR − (V + Q)T Ve kmax + kV T Ve − Ir kmax + kQ Ve kmax

≤ r2 ω 2 (1 + ω) +

2kEk2
+ω.
δ

38

(52)

We use the inequality rω < 1/2 to simplify the above bound:
1 1
kR − Ir kmax ≤ r2 ω 2 (1 + ω) + ω + 2kEk2 /δ ≤ ( + + 1)rω + 2kEk2 /δ
2 4
≤ 2rω + 2kEk2 /δ .

(53)

We are now ready to bound kVe − V kmax . In (35), we use the bounds (48), (53), (31) to
obtain
kVe − V kmax = kV (R − Ir ) + (V − V )kmax ≤ rkV kmax kR − Ir kmax + kV − V kmax
p
√
√
≤ 2r rµ/d (2rω + 2kEk2 /δ) + 2r µ ω/ d
(4r5/2 µ1/2 + 2rµ1/2 )ω 4r3/2 µ1/2 kEk2
√
√
+
d
δ d
48(1 + rµ)r5/2 µ1/2 κ 4r3/2 µ1/2 kEk2
√
√
≤
+
.
(|λr | − ε) d
δ d
≤

Using a trivial inequality κ ≤

√

rµ τ , the above bound leads to

kVe − V kmax = O

B



r3/2 µ1/2 kEk2 
r 4 µ2 τ
√ +
√
.
(|λr | − ε) d
δ d

Proofs for Section 2.2

Recall the definitions of µ0 , τ0 , κ0 and ε0 in Section 5.2. Similar to the symmetric case,
we will use the following easily verifiable inequalities.
κ0 ≤

√

rµ0 τ0 ,

p
1/2
p
kEk2 ≤
d1 /d2 kEk∞ · d2 /d1 kEk1
≤ τ0 .

(54)

Lemma B.1. Parallel to Lemma A.1, we have
kH d kw ≤ (1 + rµ0 )κ0 ,
where κ0 =

√
√
d1 kEV kmax ∨ d2 kE T U kmax as defined.

Proof of Lemma B.1. Recall H d = V⊥d (V⊥d )T E d V d = E d V d − V d (V d )T E d V d .

39

Note

V d (V d )T = diag(U U T , V V T ) and kU U T kmax ≤ rµ(U )/d1 , kV V T kmax ≤ rµ(V )/d2 . Thus,
kH d kw ≤ kE d V d kw + kV d (V d )T E d V d kw
≤ (1 + d1 kU U T kmax ∨ d2 kV V T kmax ) kE d V d kw ≤ (1 + rµ0 )κ0 .

√
d
Lemma B.2. Parallel to Lemma A.2, if σr > 2κ0 r µ0 , then L1 is a non-degenerate matrix.
Furthermore, we have the following bound
inf

kQd0 kw =1

d

d

kQd0 L1 − L2 Qd0 kw ≥ σr − 4rµ0 (τ0 + rκ0 ) − ε0 ,

(55)

where Qd0 ∈ R(d1 +d2 )×2r .
d
kw ≤
Proof of Lemma B.2. Following similar derivations with Lemma A.2, we have kQd0 E11
√
2κ0 r µ0 , and for any matrix Qd0 ∈ R(d1 +d2 )×2r with kQd0 kw = 1,
d

kL2 Qd0 kw = k(Ad − Adr )Qd0 + V⊥d (V⊥d )T E d V⊥d (V⊥d )T Qd0 kw ≤ ε0 + (1 + rµ0 )(τ0 + rκ0 ).
This can be checked by expressing Qd0 as a block matrix and expand the matrix multiplication.
In particular, one can verify that (i) k(Ad − Adr )Qd0 kw ≤ ε0 ; (ii) For any matrix M with
d1 + d2 rows, kV d (V d )T M kw ≤ rµ0 kM kw ; (iii) kE d Qd0 kw ≤ τ0 ; (iv) kE d V d (V d )T Qd0 kw ≤ rκ0 .
Moreover, kQd0 Λd1 kw ≥ σr kQd0 kw ≥ σr . Thus,
inf

kQd0 kw =1

d

d

kQd0 L1 − L2 Qd0 kw ≥ σr − 4rµ0 (τ0 + rκ0 ) − ε0 ,
d

which is the desired inequality in the lemma. In addition, L1 is non-degenerate if σr >
√
2κ0 r µ0 > 0.
d

Lemma B.3. Parallel to Lemma A.3, there is a solution Q ∈ B0 to the system (37) such
that if σr − ε0 > 16rµ0 (τ0 + rκ0 ), then
d

kQ kw ≤

8(1 + rµ0 )κ0
.
3(σr − ε0 )

Proof of Lemma B.3. We again invoke Lemma 5.2. Let B be the space R(d1 +d2 )×2r equipped
with the weighted max-norm k · kw . We also define B0 as a subspace of B consisting of
matrices of the form V⊥d Qd where Qd has size (d1 + d2 − 2r) × 2r. Let the linear operator
d
d d
d d
T d be T d Q := Q L1 − L2 Q . First notice from Lemma B.2, T d is a linear operator with
40

bounded inverse, i.e., β := k(T d )−1 k−1
w is bounded from below by
β ≥ σr − 4rµ0 (τ0 + rκ0 ) − ε0 .
d

d

d

Let ϕ be a map given by ϕ(Q ) = Q (H d )T Q . Note that H d ∈ B. Using the (easily
verifiable) inequality
kM1 M2T M3 kw ≤ 2rkM1 kw kM2T M3 kmax ≤ 4rkM1 kw kM2 kw kM3 kw

∀ M1 , M2 , M3 ∈ B, (56)

we derive, by the bound on kH d kw (Lemma B.1), that
d

d

d

kϕ(Q )kw ≤ 4rkH d kw kQ k2w ≤ 4r(1 + rµ0 )κ0 kQ k2w .
Moreover, using the inequality (56) and the bound on kH d kw (Lemma B.1),
d

d

d

d

d

d

d

d

kϕ(Q1 ) − ϕ(Q2 )kw ≤ 4rkQ1 kw kH d kw kQ1 − Q2 kw + 4rkQ1 − Q2 kw kH d kw kQ2 kw
d

d

d

d

≤ 8r(1 + rµ0 )κ0 max{kQ1 kw , kQ2 kw }kQ1 − Q2 kw .
Thus, we can choose η = 4r(1 + rµ0 )κ0 , and the condition (33) in Lemma 5.2 is satisfied.
To ensure 4ηkH d kw < β 2 , it suffices to require (again by Lemma B.1),

2
16r(1 + rµ0 )2 κ20 < σr − 4rµ0 (τ0 + rκ0 ) − ε0 .
It is easily checkable that the above inequality holds when σr − ε0 > 16rµ0 (τ0 + rκ0 ). Under
this condition, by Lemma 5.2,
d

kQ kw ≤

2(1 + rµ0 )κ0
2(1 + rµ0 )κ0
8(1 + rµ0 )κ0
2kH d kw
≤
≤
≤
,
β
σr − 4rµ(τ0 + rκ0 ) − ε0
σr − ε0 − (σr − ε0 )/4
3(σr − ε0 )

which completes the proof.
d

Proof of Lemma 5.4. The first claim of the lemma (existence of Q and its max-norm
d
bound) follows from Lemma B.3. To prove the second claim, we split V − V d into two
parts:


d
d
d
d
d
d
V − V d = V d (I2r + (Q )T Q )−1/2 − I2r + Q (I2r + (Q )T Q )−1/2 ,

(57)

√
Note κ0 ≤ τ0 rµ0 (see (54)). It can be checked that the condition σr − ε0 > 16rµ0 (τ0 + rκ0 )
d
d
d
d
d
d
implies rω0 < 1/3. Since k(Q )T Q kmax ≤ 2ω02 and k(Q )T Q k2 ≤ 2rk(Q )T Q kmax ≤

41

4rω02 < 1/2, similar to Lemma A.4, we have
d

d

d T

d −1/2

kmax ≤ 3/2.

d T

d −1/2



k(I2r + (Q )T Q )−1/2 − I2r kmax ≤ 4rω02 ,
k(I2r + (Q ) Q )

(58)

This yields
d



d

d

d

kV − V kw = kV (I2r + (Q ) Q )
− I2r kw + kQ (I2r + (Q )T Q )−1/2 kw
p
√
d
(59)
≤ 2r2 µ0 4rω02 + 2r · 3/2 · kQ kw ≤ 8 µ0 ω02 r2 + 3 ω0 r
√
√
√
≤ 8 µ0 ω0 r/3 + 3 µ0 ω0 r ≤ 6 µ0 rω0 .
(60)
d

d

1/2

Lemma B.4. Suppose σr − ε0 > max{16rµ0 (τ0 + rκ0 ), 64r3/2 µ0 (1 + rµ0 )κ0 }. Then, there
exists an orthogonal matrix RV ∈ Rr×r (or RU ) such that the column vectors of V RV (and
e
U RU ) are the top r right (and left) singular vectors of A.
Proof of Lemma B.4. Similar to the proof of Lemma 5.3, we will prove d(V , V ) < 1/2 and
d(Ve , V ) ≤ 1/2, which would then imply that V and Ve are the same only up to an orthogonal
e , and we will leave out its proof.
transformation. The same is true for U and U
By Weyl’s inequality for singular values (also known as Mirsky’s theorem (Mirsky, 1960)),
for any i, |e
σi − σi | ≤ kEk. By Wedin’s perturbation bounds for singular vectors (Wedin,
1972),
kEk2
.
d(Ve , V ) ≤
σr − kEk2
Note that kEk2 ≤ τ0 (see (54)) Under the assumption in the lemma, clearly σr − ε0 > 3τ0 ,
√
d
and we have d(Ve , V ) ≤ 1/2. Moreover, by Lemma 5.4, we have kV − V d kw ≤ 6rω0 µ0 .
d
Note that each column vector of V d and V are (d1 + d2 )-dimensional. Looking at the last
p
d2 dimensions, we have kV − V kmax ≤ 6rω0 µ0 /d1 .
p
1/2
d(V , V ) ≤ 2kV − V k ≤ 2 rd1 kV − V kmax ≤ 12r3/2 µ0 ω0 .
Under the assumption of the lemma, d(V , V ) ≤ 1/2. Therefore, we deduce d(Ve , V ) = 0, and
conclude that there exists an orthogonal matrix RV ∈ Rr×r such that Ve = V RV .

Proof of Theorem 2.3. Lemma 5.4, together with Lemma B.4, implies Theorem 2.3.

42

Proof of Theorem 2.4. Similar to the proof of Theorem 2.2, we first split the difference
Ve d − V d :
d
d
(61)
Ve d − V d = V (Rd − I2r ) + (V − V d ).
To bound the first term, note that under our assumption, rω0 < 1/3 (derived in the proof
√
d
of Lemma 5.4), it is easy to check kV kw ≤ 3 rµ0 . We rewrite the matrix Rd as
d

d

d

Rd = (I2r + (Q )T Q )−1/2 (V d + Q )T Ve d .
d

Notice that k(Q )T Ve d kmax ≤

√
√
d
2kQ kw ≤ 2 ω0 , k(Ve d )T V d kmax ≤ 1 and

d

d

d

d

kRd − (V d + Q )T Ve d kmax ≤ 2rk(I2r + (Q )T Q )−1/2 − I2r kmax k(V d + Q )T Ve d kmax
√
≤ 8r2 ω02 (1 + 2 ω0 ) .
where we used (58). Following the same derivations as in the proof of Theorem 2.2, and using
the (easily verifiable) fact kE d k2 = kEk2 , we can bound k(V d )T Ve d − I2r kmax by 2kEk2 /δ0 .
Thus, using rω0 ≤ 1/3, under δ0 > 2kEk2 , we have
d
d
kRd − I2r kmax ≤ kRd − (V d + Q )T Ve d kmax + k(V d )T Ve d − I2r kmax + k(Q )T Ve d kmax
√
√
2kEk2 √
2kEk2
+ 2 ω0 < 4 2 rω0 +
.
(62)
≤ 8r2 ω02 (1 + 2 ω0 ) +
δ0
δ0

Finally, in order to bound kVe d − V d kw , we use (60), (61) and (62), and derive
d

d

d

d

kVe d − V d kw = kV (Rd − I2r ) + (V − V d )kw ≤ 2rkV kw kRd − I2r kmax + kV − V d kw

2kEk2 
kEk2  √
√  √
√
≤ 6r rµ0 4 2 rω0 +
+ 6r µ0 ω0 ≤ 40r5/2 ω0 + 12r3/2
· µ0
δ0
δ0
1/2
1/2
107r5/2 µ0 (1 + rµ0 )κ0 12r3/2 µ0 kEk2
≤
+
σr − ε0
δ0
 r 4 µ2 τ

3/2 1/2
r µ0 kEk2
0 0
+
.
=O
σ r − ε0
δ0
This completes the proof.

C

Proofs for Section 3

Proof of Proposition 3.2. Note first by Weyl’s inequality, |λi − λi | ≤ kΣu k ≤ C. So this
implies that λi = λi (B T B)  d if and only if λi = λi (Σ)  d for i ≤ r. And furthermore the
eigenvalues of B T B/d are distinct if and only if min1≤i6=j≤r |λi (Σ) − λj (Σ)|/λj (Σ) > 0.
43

To prove the equivalency of bounded kBkmax and bounded coherence. We first prove
the necessary condition. Again from Weyl’s inequality, λi (Σ) ≤ C for i ≥ r + 1. If µ(V ) is
P
bounded, Σii must also be bounded, since Σii ≤ rj=1 vij2 λj (Σ) + λr+1 (Σ) ≤ C(µ(V ) + 1).
Therefore kbi k2 ≤ kbi k2 + (Σu )ii = Σii implies kBkmax is bounded. Namely, the factors are
pervasive.
On the contrary, if pervasiveness holds, we need to prove that µ(V ) is bounded. Let
B = (eb1 , . . . , ebr ). Obviously λi = kebi k2  d and v i = ebi /kebi k. Without loss of generality,
√
assume λ̄i ’s are decreasing. So kv i k∞ ≤ kBkmax /kebi k ≤ C/ d and µ(V ) ≤ C where V =
(v 1 , . . . , v r ). By Theorem 2.2,
kv i − vi k∞ ≤ C

kΣu k∞
√ ,
γ d

where γ = min{λi − λi+1 : 1 ≤ i ≤ r}  d with the convention λr+1 = 0. Hence, we have
√
kvi k∞ ≤ C/ d, which implies bounded coherence µ(V ).

References
Antoniadis, A. and Fan, J. (2001). Regularization of wavelet approximations. Journal
of the American Statistical Association 96.
Bai, J. (2003). Inferential theory for factor models of large dimensions. Econometrica 71
135–171.
Baik, J., Ben Arous, G. and Péché, S. (2005). Phase transition of the largest eigenvalue
for nonnull complex sample covariance matrices. Annals of Probability 1643–1697.
Bartholomew, D. J., Knott, M. and Moustaki, I. (2011). Latent variable models and
factor analysis: A unified approach, vol. 904. John Wiley & Sons.
Berthet, Q. and Rigollet, P. (2013). Optimal detection of sparse principal components
in high dimension. The Annals of Statistics 41 1780–1815.
Bickel, P. J. and Levina, E. (2008). Covariance regularization by thresholding. The
Annals of Statistics 2577–2604.
Biswas, P. and Ye, Y. (2004). Semidefinite programming for ad hoc wireless sensor network
localization. In Proceedings of the 3rd international symposium on Information processing
in sensor networks. ACM.

44

Borg, I. and Groenen, P. J. (2005). Modern multidimensional scaling: Theory and
applications. Springer Science & Business Media.
Cai, T. T. and Zhang, A. (2016). Rate-optimal perturbation bounds for singular subspaces
with applications to high-dimensional statistics. arXiv preprint arXiv:1605.00353 .
Candès, E. J., Li, X., Ma, Y. and Wright, J. (2011). Robust principal component
analysis? Journal of the ACM (JACM) 58 11.
Catoni, O. (2012). Challenging the empirical mean and empirical variance: a deviation
study. In Annales de l’Institut Henri Poincaré, Probabilités et Statistiques, vol. 48. Institut
Henri Poincaré.
Chamberlain, G. and Rothschild, M. (1982). Arbitrage, factor structure, and meanvariance analysis on large asset markets.
Chandrasekaran, V., Sanghavi, S., Parrilo, P. A. and Willsky, A. S. (2011).
Rank-sparsity incoherence for matrix decomposition. SIAM Journal on Optimization 21
572–596.
Chen, Y. (2015). Incoherence-optimal matrix completion. IEEE Transactions on Information Theory 61 2909–2923.
Choi, K. and Marden, J. (1998). A multivariate version of kendall’s τ . Journal of
Nonparametric Statistics 9 261–293.
Davis, C. and Kahan, W. M. (1970). The rotation of eigenvectors by a perturbation. iii.
SIAM Journal on Numerical Analysis 7 1–46.
De La Torre, F. and Black, M. J. (2003). A framework for robust subspace learning.
International Journal of Computer Vision 54 117–142.
Doherty, L., Pister, K. S. and El Ghaoui, L. (2001). Convex position estimation in
wireless sensor networks. In INFOCOM 2001. Twentieth Annual Joint Conference of the
IEEE Computer and Communications Societies. Proceedings. IEEE, vol. 3. IEEE.
Fama, E. F. and French, K. R. (1993). Common risk factors in the returns on stocks
and bonds. Journal of Financial Economics 33 3–56.
Fan, J., Fan, Y. and Lv, J. (2008). High dimensional covariance matrix estimation using
a factor model. Journal of Econometrics 147 186–197.

45

Fan, J., Li, Q. and Wang, Y. (2017). Robust estimation of high-dimensional mean
regression. Journal of Royal Statistical Society, Series B 79 247–265.
Fan, J., Liao, Y. and Mincheva, M. (2013). Large covariance estimation by thresholding
principal orthogonal complements. Journal of the Royal Statistical Society: Series B 75
1–44.
Fan, J., Liu, H. and Wang, W. (2017+). Large covariance estimation through elliptical
factor models. Annals of Statitistics to appear.
Fan, J., Wang, W. and Zhong, Y. (2017). Robust covariance estimation for approximate
factor models. Journal of Econometrics .
Fang, K.-T., Kotz, S. and Ng, K. W. (1990). Symmetric multivariate and related
distributions. Chapman and Hall.
Gupta, A. K., Varga, T. and Bodnar, T. (2013). Elliptically contoured models in
statistics and portfolio theory. Springer.
Han, F. and Liu, H. (2014). Scale-invariant sparse PCA on high-dimensional meta-elliptical
data. Journal of the American Statistical Association 109 275–287.
Hsu, D. and Sabato, S. (2014). Heavy-tailed regression with a generalized median-ofmeans. In Proceedings of the 31st International Conference on Machine Learning (ICML14).
Hu, W., Tan, T., Wang, L. and Maybank, S. (2004). A survey on visual surveillance of
object motion and behaviors. Systems, Man, and Cybernetics, Part C: IEEE Transactions
on Applications and Reviews 34 334–352.
Huber, P. J. (1964). Robust estimation of a location parameter. The Annals of Mathematical Statistics 35 73–101.
Johnstone, I. M. and Lu, A. Y. (2009). On consistency and sparsity for principal
components analysis in high dimensions. Journal of the American Statistical Association
104 682–693.
Jolliffe, I. (2002). Principal component analysis. Wiley Online Library.
Ma, Z. (2013). Sparse principal component analysis and iterative thresholding. The Annals
of Statistics 41 772–801.

46

Mirsky, L. (1960). Symmetric gauge functions and unitarily invariant norms. The quarterly
journal of mathematics 11 50–59.
Ng, A. Y., Jordan, M. I., Weiss, Y. et al. (2002). On spectral clustering: Analysis
and an algorithm. Advances in neural information processing systems 2 849–856.
Oliver, N. M., Rosario, B. and Pentland, A. P. (2000). A bayesian computer vision
system for modeling human interactions. IEEE Transactions on Pattern Analysis and
Machine Intelligence 22 831–843.
O’Rourke, S., Vu, V. and Wang, K. (2013). Random perturbation of low rank matrices:
Improving classical bounds. arXiv preprint arXiv:1311.2657 .
Paulsen, V. (2002). Completely bounded maps and operator algebras, vol. 78. Cambridge
University Press.
Rennie, J. D. and Srebro, N. (2005). Fast maximum margin matrix factorization for
collaborative prediction. In Proceedings of the 22nd international conference on Machine
learning. ACM.
Rohe, K., Chatterjee, S. and Yu, B. (2011). Spectral clustering and the highdimensional stochastic blockmodel. The Annals of Statistics 1878–1915.
Rothman, A. J., Levina, E. and Zhu, J. (2009). Generalized thresholding of large
covariance matrices. Journal of the American Statistical Association 104 177–186.
Stewart, G. W. (1990). Matrix perturbation theory .
Stock, J. and Watson, M. (2002). Forecasting using principal components from a large
number of predictors. Journal of the American Statistical Association 97 1167–1179.
Tropp, J. A. (2012). User-friendly tail bounds for sums of random matrices. Foundations
of Computational Mathematics 12 389–434.
Vershynin, R. (2010). Introduction to the non-asymptotic analysis of random matrices.
arXiv preprint arXiv:1011.3027 .
Vu, V. (2011). Singular vectors under random perturbation. Random Structures & Algorithms 39 526–538.
Vu, V. Q. and Lei, J. (2012). Minimax rates of estimation for sparse pca in high dimensions.
arXiv preprint arXiv:1202.0786 .
47

Wang, W. and Fan, J. (2017+). Asymptotics of empirical eigen-structure for ultra-high
dimensional spiked covariance model. Annals of Statitistics to appear.
Wedin, P.-Å. (1972). Perturbation bounds in connection with singular value decomposition. BIT Numerical Mathematics 12 99–111.
Yu, Y., Wang, T. and Samworth, R. J. (2015). A useful variant of the davis–kahan
theorem for statisticians. Biometrika 102 315–323.
Yuan, M., Ekici, A., Lu, Z. and Monteiro, R. (2007). Dimension reduction and
coefficient estimation in multivariate linear regression. Journal of the Royal Statistical
Society: Series B (Statistical Methodology) 69 329–346.

48

