arXiv:1607.06017v2 [math.OC] 26 Nov 2016

Doubly Accelerated Methods for
Faster CCA and Generalized Eigendecomposition
Zeyuan Allen-Zhu
zeyuan@csail.mit.edu
Institute for Advanced Study

Yuanzhi Li
yuanzhil@cs.princeton.edu
Princeton University

July 20, 2016∗

Abstract
We study k-GenEV, the problem of finding the top k generalized eigenvectors, and k-CCA,
the problem of finding the top k vectors in canonical-correlation analysis. We propose algorithms
LazyEV and LazyCCA to solve the two problems with running times linearly dependent on the
input size and on k.
Furthermore, our algorithms are doubly-accelerated : our running times depend only on the
square root of the matrix condition number, and on the square root of the eigengap. This is the
first such result for both k-GenEV or k-CCA.
We also provide the first gap-free results, which
√
provide running times that depend on 1/ ε rather than the eigengap.

1

Introduction

The Generalized Eigenvector (GenEV) problem and the Canonical Correlation Analysis (CCA)
are two fundamental problems in scientific computing, machine learning, operations research, and
statistics. Algorithms solving these problems are often used to extract features to compare largescale datasets, as well as used for problems in regression [17], clustering [10], classification [18],
word embeddings [11], and many others.
Given two symmetric matrices A, B ∈ Rd×d where B is positive definite. The GenEV problem
is to find generalized eigenvectors v1 , . . . , vd where each vi satisfies
 >
v Bv = 1
>
vi ∈ arg max v Av such that
v > Bvj = 0 ∀j ∈ [i − 1]
v∈Rd
def

The values λi = vi> Avi are known as the generalized eigenvalues, and it satisfies |λ1 | ≥ · · · |λd |.
Following the tradition of [13, 29], we assume without loss of generality that λi ∈ [−1, 1].

Given matrices X ∈ Rn×dx , Y ∈ Rn×dy and denoting by Sxx = n1 X > X, Sxy = n1 X > Y ,
Syy = n1 Y > Y , the CCA problem is to find canonical-correlation vectors {(φi , ψi )}ri=1 where r =
min{dx , dy } and each pair
 >
 >
φ Sxx φ = 1 ∧ φ> Sxx φj = 0 ∀j ∈ [i − 1]
(φi , ψi ) ∈ arg max
φ Sxy ψ
such that
ψ > Syy ψ = 1 ∧ ψ > Syy ψj = 0 ∀j ∈ [i − 1]
φ∈Rdx ,ψ∈Rdy
def

The values σi = φ>
i Sxy ψi ≥ 0 are known as the canonical-correlation coefficients, and it always
satisfies 1 ≥ σ1 ≥ · · · ≥ σr ≥ 0.
∗

First appeared on arXiv on this date. In this new version, we have stated more clearly why this paper has
outperformed relevant previous results, and included discussions for doubly-stochastic methods.

1

It is a folklore that solving CCA exactly can be reduced to solving GenEV exactly, if one defines
> , 0]] ∈ Rd×d for d def
B = diag{Sxx , Syy } ∈ Rd×d and A = [[0, Sxy ]; [Sxy
= dx + dy , see Lemma 2.3.

Despite the fundamental importance and the frequent necessity in applications, there are few
results on obtaining provably efficient algorithms for GenEV and CCA until very recently. In the
breakthrough result of Ma, Lu and Foster [20], they proposed to study algorithms to find top k
generalized eigenvectors (k-GenEV) or top k canonical-correlation vectors (k-CCA). They designed
an alternating minimization algorithm whose running time is only linear in terms of the number of
non-zero elements of the matrix (that we denote by nnz(A) for a matrix A in this paper), and also
nearly-linear in k. Such algorithms are very appealing because in real-life applications, it is often
only relevant to obtain top correlation vectors, as opposed to the less meaningful vectors in the
directions where the datasets do not correlate. Unfortunately, the method of Ma, Lu and Foster
has a running time that linearly scales with κ and 1/gap, where
• κ ≥ 1 is the condition number of matrix B in GenEV, or of matrices X > X, Y > Y in CCA; and
• gap ∈ [0, 1) is the (relative) eigengap

λk −λk+1
λk

in GenEV, or

σk −σk+1
σk

in CCA.

These parameters are usually not constants and scale with the problem size.
Challenge 1: Acceleration
For many easier scientific computing problems, we are able to design algorithms that have accelerated dependencies on κ and 1/gap. As two concrete examples, k-PCA can be solved with a running
√
time linearly in 1/ gap as opposed to 1/gap [16]; computing B −1 w for a vector w can be solved in
√
time linearly in κ as opposed to κ, where κ is the condition number of matrix B [9, 23, 27].
Therefore, can we obtain doubly-accelerated methods for k-GenEV and k-CCA, meaning
√
√
that the running times linearly scale with both κ and 1/ gap? Before this paper, for the general
case k > 1, the method of Ge et al. [15] made acceleration possible for parameter κ, but not for
parameter 1/gap (see Table 2).
Challenge 2: Gap-Freeness
Since gap can be even zero in the extreme case, can we design algorithms that do not scale with
1/gap? Recall that this is possible for the easier task of k-PCA. The block Krylov method [22] runs
√
√
in time linear in 1/ ε as opposed to 1/ gap, where ε is the approximation ratio.
There is no gap-free result previously known for k-GenEV or k-CCA even for k = 1.
Challenge 3: Stochasticity
For matrix-related problems, one can often obtain a stochastic running time which requires some
notations to describe. Consider the simple task of computing B −1 w for some vector w, where recall
√
accelerated methods solve this problem with a running time linearly in κ for κ being the condition
number of B. If B = n1 X > X is given as the covariance matrix of X ∈ Rn×d , then (accelerated)
p
stochastic gradient methods can be used to compute B −1 w in a time linearly in (1+ κ0 /n) instead


√
maxi∈[n] {kXi k2 }
of κ, where κ0 =
∈
κ,
nκ
and Xi is the i-th row of X. (See Lemma 2.6.) Since
λ
(B)
min
p
√
1 + κ0 /n ≤ O( κ), stochastic methods are no slower than non-stochastic ones.
Therefore, can we obtain a similar (but doubly-accelerated!) stochastic method for k-CCA?1
Note that, if the doubly-accelerated requirement is dropped, this task is easier and indeed possible,
see Ge et al. [15]. However, since their stochastic method is not doubly-accelerated, in certain
parameter regimes, it can run slower than non-stochastic ones (even for k = 1, see Table 1).
1

Note that a similar problem can be also asked for k-GenEV when A and B are both given in their covariance
matrix forms. We refrain from doing it in this paper for notational simplicity.

2

Remark. In general, if designed properly,
• Accelerated results are better because they are never slower than non-accelerated ones.
• Gap-free results are better because they imply gap-dependent ones.2

• Stochastic results are better because they are never slower than non-stochastic ones.

1.1

Our Main Results

We provide algorithms LazyEV and LazyCCA that are doubly-accelerated, gap-free, and stochastic.3
For the general k-GenEV problem, our LazyEV can be implemented to run in time
 knnz(B)√κ knnz(A) + k 2 d 
 knnz(B)√κ knnz(A) + k 2 d 
e
e
√
√
+
or O
+
O
√
√
gap
gap
ε
ε
in the gap-dependent and gap-free cases respectively. Since our running time only linearly depends
√
√
√
on κ and gap (resp. ε), our algorithm LazyEV is doubly-accelerated.
For the general k-CCA problem, our LazyCCA can be implemented to run in time
 knnz(X, Y ) · 1 + pκ0 /n + k 2 d 
 knnz(X, Y ) · 1 + pκ0 /n + k 2 d 
e
e
√
O
or O
√
gap
ε

in the gap-dependent and gap-free cases respectively. Here, nnz(X, Y ) = nnz(X) + nnz(Y ) and
maxi {kXi k2 ,kYi k2 }
κ0 = 2λmin
(diag{Sxx ,Syy }) where Xi or Yi is the i-th row vector of X or Y . Therefore, our algorithm
LazyCCA is doubly-accelerated and stochastic.
We fully compare our running time with prior work in Table 1 and Table 2, and summarize our
main contributions below.
• For k > 1, we have outperformed all relevant prior works (see Table 2). Moreover, no known
method was doubly-accelerated even in the non-stochastic setting.
• For k ≥ 1, we have obtained the first gap-free running times.

• Even for k = 1, we have outperformed most of the state-of-the-arts (see Table 1).
Other Contributions. Besides the aforementioned running time improvements, we summarize
some other virtues of our algorithms as follows:
• For GenEV, our LazyEV distinguishes positive generalized eigenvalues from negative ones. For
instance, if A has two generalized eigenvectors v1 , v2 with respect to B, one with eigenvalue
λ and the other with −λ. Then, previous result GenELin only finds the subspace spanned by
v1 , v2 but cannot distinguish v1 from v2 .
• For CCA with k > 1, previous result CCALin only outputs the subspace spanned by the top k
correlation vectors but not identify which vector gives the highest correlation and so on. Our
LazyCCA provides per-vector guarantees on all the top k correlation vectors, see Theorem 6.3.
• Our LazyEV and LazyCCA reduce the underlying non-convex problem to multiple calls of
quadratic minimization. Since quadratic minimization is a well-studied convex optimization
problem, many efficient and robust algorithms can be found. In contrast, previous results for
the k > 1 case rely on more sophisticated nonconvex optimization; and the previous work
of [29] —although uses convex optimization to solve 1-CCA— requires one to work with a
sum-of-non-convex function which is usually less efficient to minimize.
2
3

If a method depends on 1/ε then one can choose ε = gap and this translates to a gap-dependent running time.
Recalling Footnote 1, for notational simplicity, we only state our k-GenEV result in non-stochastic running time.

3

Problem

Paper
GenELin [15]

1-GenEV

LazyEV Theorem 4.3
LazyEV Theorem 4.4

Problem

Paper
AppGrad [20]
CCALin [15]
ALS [29]
SI [29]
CCALin [15]
ALS [29]

1-CCA

LazyCCA Theorem 6.2
LazyCCA Theorem 6.3
SI [29]

Running time (× for beaten) gap-free? negative EV?
√

e nnz(B) κB + nnz(A)
O
no
no
×
gap
gap
√

κB
e nnz(B)
√
√
O
+ nnz(A)
no
yes
gap
gap
√

κB
e nnz(B)
√
√
+ nnz(A)
yes
yes
O
ε
ε
Running time (× for beaten) gap-free?

e κ
nnz(X, Y ) · O
×
no
gap
√ 
e κ
nnz(X, Y ) · O
no
×
gap
√ 
e κ2
nnz(X, Y ) · O
×
no
gap
√

κ
e √
nnz(X, Y ) · O
no
×
gap·σ1
√ 0 
1+
κ
/n
e
no
nnz(X, Y ) · O
×
gap
√
0 /n 
1+
κ
e
no
nnz(X, Y ) · O
×
gap2
√ 0 
e 1+√ κ /n
no
nnz(X, Y ) · O
√gap0 
e 1+ √κ /n
yes
nnz(X, Y ) · O
ε


√
1/4
e 1 + √κ0 /n
nnz(X, Y ) · O
gap·σ1
no
(see Remark 3)

LazyCCA Theorem 6.2
LazyCCA Theorem 6.3

e 1+
nnz(X, Y ) · O
e 1+
nnz(X, Y ) · O

Table 1: Performance comparison on 1-GenEV and 1-CCA.
(B)
2
In GenEV, gap = λ1λ−λ
∈ [0, 1] and κB = λλmax
> 1.
1
min (B)
In CCA, gap =

σ1 −σ2
σ1

∈ [0, 1], κ =

λmax (diag{Sxx ,Syy })
λmin (diag{Sxx ,Syy })

> 1, κ0 =

√
κ0 /n1/4 
√
gap·σ1
√
κ0 /n1/4 
√
ε·σ1

2 maxi {kXi k2 ,kYi k2 }
λmin (diag{Sxx ,Syy })

stochastic?
no
no
no
no
yes
yes
yes
yes
doubly

no

doubly

yes

doubly

∈ [κ, 2nκ], and σ1 ∈ [0, 1].

Remark 1. Stochastic methods depend on modified condition number κ0 ; the reason κ0 ∈ [κ, 2nκ] is in Def. 2.4.
p
Remark 2. All non-stochastic CCA methods in this table have been outperformed because 1 + κ0 /n ≤ O(κ).

Remark 3. Doubly-stochastic methods are not necessarily interesting. We discuss them in Section 1.2.

Remark 4. Some CCA methods have a running time dependency on σ1 ∈ [0, 1], and this is intrinsic and cannot be
removed. In particular, if we scale the data matrix X and Y , the value σ1 stays the same.
Remark 5. The only (non-doubly-stochastic) doubly-accelerated
p
p method before our work is SI [29] (for 1-CCA
only). Our LazyEV is faster than theirs by a factor Ω( nκ/κ0 × 1/σ1 ). Here, nκ/κ0 ≥ 1/2 and 1/σ1 ≥ 1 are two
scaling-invariant quantities usually much greater than 1.

4

Problem

Paper
GenELin [15]

k-GenEV
Problem

e
O

e
LazyEV Theorem 4.3 O
e
LazyEV Theorem 4.4 O
Paper

AppGrad [20]
CCALin [15]
CCALin [15]
k-CCA

Running time
√
knnz(B) κB
gap
√
knnz(B) κB
√
gap
√
knnz(B) κB
√
ε

(× for beaten) gap-free? negative EV?
d
+ knnz(A)+k
no
no
×
gap
2 
d
√
+ knnz(A)+k
no
yes
gap
2 
d
√
+ knnz(A)+k
yes
yes
ε
2

Running time
e
O
e
O
e
O

e
LazyCCA Theorem 6.2 O
e
LazyCCA Theorem 6.3 O

2

(× for beaten) gap-free?

d

knnz(X,Y )·κ+k
(local
gap
√
2 
knnz(X,Y )· κ+k d
gap

√
knnz(X,Y )· 1+ κ0 /n +k2 d 
gap

√
knnz(X,Y )· 1+ κ0 /n +k2 d 
√
gap


√

knnz(X,Y )· 1+
√
ε

κ0 /n

+k

2

d

conv.)

×

no

no

×

no

no

×

no

yes

no

yes

yes

yes

no

doubly

yes

doubly


κ0
e knnz(X, Y )· 1+ √
LazyCCA Theorem 6.2 O
gap·σk ·(nnz(X,Y )/kd)1/4
√

κ0
e knnz(X, Y ) · 1 + √
LazyCCA Theorem 6.3 O
ε·σk ·(nnz(X,Y )/kd)1/4

Table 2: Performance comparison on k-GenEV and k-CCA.
λ −λ
(B)
In GenEV, gap = k λkk+1 ∈ [0, 1] and κB = λλmax
> 1.
min (B)
In CCA, gap =

σk −σk+1
σk

∈ [0, 1], κ =

λmax (diag{Sxx ,Syy })
λmin (diag{Sxx ,Syy })

> 1, κ0 =

√

2 maxi {kXi k2 ,kYi k2 }
λmin (diag{Sxx ,Syy })

stochastic?

∈ [κ, 2nκ], and σk ∈ [0, 1].

Remark 1. Stochastic methods depend on a modified condition number κ0 . The reason κ0 ∈ [κ, 2nκ] is in Def. 2.4.
p
Remark 2. All non-stochastic CCA methods in this table have been outperformed because 1 + κ0 /n ≤ O(κ).
Remark 3. Doubly-stochastic methods are not necessarily interesting. We discuss them in Section 1.2.

1.2

Our Side Results on Doubly-Stocahstic Methods

Recall that when considering acceleration, there are two parameters κ and 1/gap. One can also
design stochastic methods with respect to both parameters κ and 1/gap, meaning that
√ 0 c
κ /n
with a running time proportional to 1 + √gap
(constant c is usually 1/2)
√ 0
√
κ /n
κ
1
instead of √gap
+ √gap (stochastic) or √gap
(non-stochastic). We call such methods doublystochastic.
Unfortunately, doubly-stochastic methods are usually slower than stochastic ones. Take 1CCA as an example.
best stochastic running time (obtained exclusively by us) for 1-CCA is
√ 0 The

1+
κ
/n
e √
nnz(X, Y ) · O
. In contrast, if one uses a doubly-stochastic method —either [29] or our
gap
√

e 1 + √κ0 /n1/4 . Therefore, for 1-CCA,
LazyCCA— the running time becomes nnz(X, Y ) · O
gap·σ1
doubly-stochastic methods are faster than stochastic ones only when

κ0
σ1

≤ o(n1/2 ) .

The above condition is usually not satisfied. For instance, (1) κ0 is usually around n for most
interesting data-sets, cf. the experiments of [25]; (2) κ0 is between n1/2 and 100n in all the CCA
experiments of [29]; and (3) by Def. 2.4 it satisfies κ0 ≥ d so κ0 cannot be smaller than o(n1/2 )
unless d  n1/2 .4 Even worse, parameter σ1 ∈ [0, 1] is usually much smaller than 1. Note that σ1
is scaling invariant: even if one scales X and Y up by the same factor, σ1 remains unchanged.
4

Note that item (3) κ0 ≥ d may not hold in the more general setting of CCA, see Remark 2.5.

5

Nevertheless, in order to compare our LazyCCA framework with all relevant prior works, we
also obtain doubly-stochastic running times for k-CCA. Our running time matches that of [29]
when k = 1, and no doubly-stochastic running time for k > 1 was known before our work.

1.3

Other Related Works

For the easier task of PCA and SVD, the first gap-free result was obtained by Musco and Musco [22]
(or in the online setting by [3]), the first stochastic result was obtained by Shamir [26], and the
first accelerated stochastic result was obtained by Garber et al. [13, 14]. The shift-and-invert
preconditioning technique of Garber et al. is also used in this paper.
For another related problem PCR (principle component regression), we recently obtained an
accelerated method [4] as opposed the previously non-accelerated one [12]; however, the acceleration
techniques in [4] are not relevant to this paper.
For GenEV and CCA, many scalable algorithms have been designed recently [19–21, 28, 30].
However, as summarized by the authors of CCALin, these cited methods are more or less heuristics
and do not have provable guarantees. Furthermore, for k > 1, the AppGrad result of [20] only provides local convergence guarantees and thus requires a warm-start whose computational complexity
is not discussed in their paper.
Finally, our algorithms on GenEV and CCA are based on finding vectors one-by-one, which is
advantageous in practice because one does not need k to be known and can stop the algorithm
whenever the eigenvalues (or correlation values) are too small. Known approaches for k > 1 cases
(such as GenELin, CCALin, AppGrad) find all k vectors at once, therefore requiring k to be known
beforehand. As a separate note, these known approaches do not need the user to know the desired
accuracy a priori but our LazyEV and LazyCCA algorithms do.

2

Preliminaries

For a vector x we denote by kxk or kxk2 the Euclidean norm of x. Given a matrix A we denote
by kAk2 and kAkF respectively the spectral and Frobenius norms of A. For q ≥ 1, we denote by
kAkSq the Schatten q-norm of A. We write A  B if A, B are symmetric and A − B is positive
semi-definite (PSD), and write A  B if A, B are symmetric but A − B is positive definite (PD).
We denote by λmax (M ) and λmin (M ) the largest and smallest eigenvalue of a symmetric matrix
M , and by κM the condition number λmax (M )/λmin (M ) of a PSD matrix M .
def

Throughout this paper, for a matrix A ∈ Rn×d , we define nnz(A) = max{n, d, N } where N is the
number of non-zero entries of A. For two matrices X, Y , we denote by nnz(X, Y ) = nnz(X)+nnz(Y ),
and by Xi or Yi the i-th row vector of X or Y . We also use poly(x1 , x2 , . . . , xt ) to represent
a quantity that is asymptotically at most polynomial in terms of variables x1 , . . . , xt . Given a
column orthonormal matrix U ∈ Rn×k , we denote by U ⊥ ∈ Rn×(n−k) the column orthonormal
matrix consisting of an arbitrary basis in the space orthogonal to the span of U ’s columns.
Given a PSD matrix B and a vector v, the value v > Bv is the B-inner product. Two vectors
v, w satisfying v > Bw = 0 are known as B-orthogonal. Given a PSD matrix B, we denote by B −1
the Moore-Penrose pseudoinverse of B which is also PSD, and denote by B 1/2 the matrix square
root of B (satisfying B 1/2  0). All occurrences of B −1 , B 1/2 and B −1/2 are for analysis purpose
only. When implementing our algorithms, it only requires one to multiply B to a vector.

6

Definition 2.1 (GenEV). Given symmetric matrices A, B ∈ Rd×d where B is positive definite.
The generalized eigenvectors of A with respect to B are v1 , . . . , vd , where each vi is
)
(
n v > Bv = 1
vi ∈ arg max v > Av such that
v > Bvj = 0 ∀j ∈ [i − 1]
v∈Rd
The corresponding generalized eigenvalues λ1 , . . . , λd satisfy λi = vi> Avi which is possibly negative.

Following the tradition of [13, 29], we assume without loss of generality that λi ∈ [−1, 1].

Definition 2.2 (CCA). Given X ∈ Rn×dx , Y ∈ Rn×dy , letting Sxx = n1 X > X, Sxy = n1 X > Y ,
Syy = n1 Y > Y , the canonical-correlation vectors are {(φi , ψi )}ri=1 where r = min{dx , dy } and ∀i:
)
(
n φ> S φ = 1 ∧ φ> S φ = 0 ∀j ∈ [i − 1]
xx
xx j
>
(φi , ψi ) ∈ arg max
φ Sxy ψ such that
ψ > Syy ψ = 1 ∧ ψ > Syy ψj = 0 ∀j ∈ [i − 1]
φ∈Rdx ,ψ∈Rdy

The corresponding canonical-correlation coefficients σ1 , . . . , σr satisfy σi = φ>
i Sxy ψi ∈ [0, 1].

We emphasize here that σi always lies in [0, 1] and is scaling-invariant. When dealing with a CCA
problem, we also denote by d = dx + dy .
Lemma 2.3 (CCA to GenEV). Given a CCA problem with matrices X ∈ Rn×dx , Y ∈ Rn×dy , and
r
suppose the canonical-correlation
 coefficients are {(φi , ψi , σi )}i=1 where r = min{dx , dy }.


vectors and
0 Sxy
>
Sxy
0

Sxx 0
0 Syy

. Then, the GenEV problem of A with respect to B has
on
n  
−φi
φi
,
. The re2r eigenvalues {±σi }ri=1 and corresponding generalized eigenvectors
ψi
ψi
i=1
maining dx + dy − 2r eigenvalues are zeros.
Define A =

and B =

Definition 2.4. In CCA, let A and B be as defined in Lemma 2.3. We define condition numbers
2
2
def λmax (B)
def 2 maxi {kXi k , kYi k }
κ=
and κ0 =
(it must satisfy κ0 ∈ [κ, 2nκ] and κ0 ≥ d.)
λmin (B)
λmin (B)
P
Proof. We have λmax (B) ≤ Tr(B) = n1 i kXi k2 + kYi k2 ≤ 2 maxi {kXi k2 , kYi k2 } and thus κ0 ≥ κ.
Suppose without loss of generality that kX1 k2 = maxi {kXi k2 , kYi k2 }; then, we have λmax (B) ≥
1
1
2
2
2
0
2
2
n kX1 k = n maxi {kXi k , kYi k }. This implies κ ≤ 2nκ. Finally, 2 maxi {kXi k , kYi k } ≥ Tr(B) ≥
0
dλmin (B) and therefore κ ≥ d.

Remark 2.5. We have followed the very original definition of CCA [15, 20] by letting Sxx = n1 X > X
and Syy = n1 Y > Y . In contrast, one prior work [29] considered the slightly more general version
Sxx = γx I + n1 X > X and Syy = γy I + n1 Y > Y for some γx , γy ≥ 0. All of the results in this paper
continue to hold in this more general setting, but we refrain from doing so for notational simplicity.
(The only difference is that the parameter κ0 will no longer satisfy κ0 ≥ d in Def. 2.4.)
Lemma 2.6. Given matrices X ∈ Rn×dx , Y ∈ Rn×dy , let A and B be as defined in Lemma 2.3.
For every w ∈ Rd , Katyusha method [1] finds a vector w0 ∈ Rd satisfying kw0 − B −1 Awk ≤ ε

p

κkwk2 
.
in time
O nnz(X, Y ) · 1 + κ0 /n · log
ε

7

Algorithm 1 AppxPCA± (A, M, δ× , ε, p)

Input: A, an approximate matrix inversion method; M ∈ Rd×d , a symmetric matrix satisfying
−I  M  I; δ× ∈ (0, 0.5], a multiplicative error; ε ∈ (0, 1), a numerical accuracy parameter;
and p ∈ (0, 1), the confidence parameter.
1: w
b0 ← RanInit(d);s ← 0; λ(0)
+ δ
 wb0 is a random unit vector, see Def. 3.2
×;
 ← 136dθ
288dθ
2: m1 ← 4 log
, m2 ← log p2 ε ;
 θ is the parameter of RanInit, see Def. 3.2
p2
δ× m1
δ × m 2
ε
1
and εe2 ← 8m2 48
3: ε
e1 ← 64m1 48
4: repeat
 it satisfies m1 = T PM (8, 1/32, p) and m2 = T PM (2, ε/4, p), see Lemma B.1
5:
s ← s + 1;
6:
for t = 1 to m1 do
bt − (λ(s−1) I − M )−1 w
bt−1 ≤ εe1 ;
7:
Apply A to find w
bt satisfying w
8:

9:

10:
11:
12:
13:
14:
15:
16:
17:
18:
19:

bm1 k;
wa ← w
bm1 /kw
Apply A to find va satisfying va − (λ(s−1) I − M )−1 wa ≤ εe1 ;
for t = 1 to m1 do
bt − (λ(s−1) I + M )−1 w
bt−1 ≤ εe1 ;
Apply A to find w
bt satisfying w
wb ← w
bm1 /kw
bm1 k;
Apply A to find vb satisfying vb − (λ(s−1) I + M )−1 wb ≤ εe1 ;
∆(s) ←

1
2

1
max{wa> va ,wb> vb }−e
ε1
δ× λ(s)
12

·

and λ(s) ← λ(s−1) −

∆(s)
2 ;

until ∆(s) ≤
f ← s;
if the last wa> va ≥ wb> vb then
for t = 1 to m2 do
bt − (λ(f ) I − M )−1 w
bt−1 ≤ εe2 ;
Apply A to find w
bt satisfying w
def

return (+, w) where w = w
bm2 /kw
bm2 k.
21: else
22:
for t = 1 to m2 do
bt − (λ(f ) I + M )−1 w
bt−1 ≤ εe2 ;
23:
Apply A to find w
bt satisfying w
20:

def

return (−, w) where w = w
bm2 /kw
bm2 k.
25: end if
24:

3

Leading Eigenvector via Two-Sided Shift-and-Invert

In this section we define AppxPCA± , the multiplicative approximation algorithm for computing
the two-sided leading eigenvector of a symmetric matrix using the shift-and-invert preconditioning
framework [13, 14]. Our pseudo-code Algorithm 1 is a modification of Algorithm 5 in [13].
The main differences between AppxPCA± and Algorithm 5 of [13] are two-fold. First, given a
symmetric matrix M , AppxPCA± simultaneously considers an upper-bounding shift together with
a lower-bounding shift, and try to invert both λI − M and λI + M . This allows us to determine
approximately how close λ is to the largest and the smallest eigenvalues of M , and decrease λ
accordingly; in the end, it outputs an approximate eigenvector of M that corresponds to a negative
eigenvalue if needed. Second, we provide a multiplicative-error guarantee rather than additive as
originally appeared in [13]. Without this multiplicative-error guarantee, our final running time will
1 5
1
depend on gap·λmax
(M ) rather than gap . Of course, we believe the bulk of the credit for conceiving
5

This is why the SI method [29] depends on

1
gap·σ1

in Table 1.

8

AppxPCA± belongs to the original authors of [13, 14].
Theorem 3.1 (AppxPCA± ). Let M ∈ Rd×d be a symmetric matrix with eigenvalues 1 ≥ λ1 ≥
· · · ≥ λd ≥ −1 and corresponding eigenvectors u1 , . . . , ud . Let λ∗ = kM k2 = max{λ1 , −λd }. With
probability at least 1 − p, AppxPCA± produces a pair (sgn, w) satisfying
X
(w> ui )2 ≤ ε and w> M w ≥ (1 − δ× /2)(1 − 3ε)λ∗ , and
if sgn = +, then
i∈[d],λi ≤(1−δ× /2)λ∗

if sgn = −, then

X

i∈[d],λi ≥−(1−δ× /2)λ∗

(w> ui )2 ≤ ε

and

w> M w ≤ −(1 − δ× /2)(1 − 3ε)λ∗ .

Furthermore, the total number of oracle calls to A is O(log(1/δ× )m1 + m2 ), and each time we call
A it satisfies that

λmax (λ(s) I−M ) λmax (λ(s) I+M )
,
λmin (λ(s) I−M ) λmin (λ(s) I+M )

∈ [1, δ96
] and
×

1
1
,
λmin (λ(s) I−M ) λmin (λ(s) I+M )

≤

48
δ × λ∗ .

We remark here that, unlike the original shift-and-invert method which chooses a random
(Gaussian) unit vector in Line 1 of AppxPCA± , we allow this initial vector to be generated from an
arbitrary θ-conditioned random vector generator, defined as follows:
Definition 3.2. An algorithm RanInit(d) is a θ-conditioned random vector generator if w =
RanInit(d) is a d-dimensional unit vector and, for every p ∈ (0, 1), every unit vector u ∈ Rd , with
2
probability at least 1 − p, it satisfies (u> w)2 ≤ p9dθ .
This modification is needed in order to obtain our efficient implementations of GenEV and CCA
algorithms. One can construct θ-conditioned random vector generator as follows:
def

1/2

B
v
Proposition 3.3. Given PSD matrix B ∈ Rd×d , if we set RanInit(d) = (v>
where v is a
Bv)0.5
random Gaussian vector, then RanInit(d) is a θ-conditioned random vector generator for θ = κB .

Proof of Proposition 3.3. We have
(u> w)2 =

Tr(uu> vBv > ) ¬ Tr(uu> vBv > ) ­ λmin (B) · Tr(uu> vv > )
≥
≥
= θ(u> v)2 .
λmax (B)
λmax (B)
v > Bv

Above, ¬ isbecause v > Bv ≤ λmax (B) · kvk22 = λmax (B), and ­ follows from the fact that vBv > 
v λmin (B)I v > = λmin (B)vv > . Finally, using for instance [8, Lemma 5], it holds with probability

at least 1 − p that (u> v)2 ≥

4



p2
9d .

LazyEV: Our Algorithm for Generalized Eigendecomposition

In this section, we propose LazyEV (see Algorithm 2) to compute approximately the k “leading”
eigenvectors corresponding to the k largest absolute eigenvalues of some symmetric matrix M ∈
Rd×d . Later, we shall solve the k-GenEV problem by setting M = B −1/2 AB −1/2 and using LazyEV
to find the top k leading eigenvectors of M , which correspond to the top k leading generalized
eigenvectors of A with respect to B.
Our algorithm LazyEV is formally stated in Algorithm 2. It applies k times AppxPCA± , each
time with a multiplicative error δ× /2, and projects the matrix M into the orthogonal space with
respect to the obtained leading eigenvector. We state our main approximation theorem below.

9

Algorithm 2 LazyEV(A, M, k, δ× , εpca , p)

Input: A, an approximate matrix inversion method; M ∈ Rd×d , a matrix satisfying −I  M  I;
k ∈ [d], the desired rank; δ× ∈ (0, 1), a multiplicative error; εpca ∈ (0, 1), a numerical accuracy
parameter; and p ∈ (0, 1), a confidence parameter.
1: M0 ← M ; V0 = [];
2: for s = 1 to k do
3:
vs0 ← AppxPCA± (A, Ms−1 , δ× /2, εpca , p/k);
>
0
> )v 0 / (I − V
⊥
 project vs0 to Vs−1
4:
vs ← (I − Vs−1 Vs−1
s−1 Vs−1 )vs ;
s
5:
Vs ← [Vs−1 , vs ];
6:
Ms ← (I − vs vs> )Ms−1 (I − vs vs> )
 we also have Ms = (I − Vs Vs> )M (I − Vs Vs> )
7: end for
8: return Vk .
Theorem 4.1 (approximation of LazyEV). Let M ∈ Rd×d be a symmetric matrix with eigenvalues
λ1 , . . . , λd ∈ [−1, 1] and corresponding eigenvectors u1 , . . . , ud , and assume |λ1 | ≥ · · · ≥ |λd |.

1
1|
6
For every k ∈ [d], δ× , p ∈ (0, 1), there exists some εpca ≤ O poly(δ× , |λ|λk+1
| , d ) such that

LazyEV outputs a (column) orthonormal matrix Vk = (v1 , . . . , vk ) ∈ Rd×k which, with probability at
least 1 − p, satisfies all of the following properties. (Denote by Ms = (I − Vs Vs> )M (I − Vs Vs> ).)
(a) Correlation guarantee: kVk> U k2 ≤ ε,
where U = (uj , . . . , ud ) and j is the smallest index satisfying |λj | ≤ (1 − δ× )kMk−1 k2 .
|λk+1 |
(b) Spectral norm guarantee: |λk+1 | ≤ kMk k2 ≤ 1−δ
.
×
1
>
(c) Rayleigh quotient guarantee: (1 − δ× )|λk | ≤ |vk M vk | ≤ 1−δ
|λk |.
×
(d) Schatten-q norm guarantee: for every q ≥ 1,
d

(1 + δ× )2  X q 1/q
(1 + δ× )2
min
k(I −V V > )M (I −V V > )kSq
λ
kMk kSq ≤
=
i
2
2
d×k
>
(1 − δ× )
(1 − δ× ) V ∈R ,V V =I

.

i=k+1

The next theorem states that, if M = B −1/2 AB −1/2 , then LazyEV can be implemented without
ever needing to compute B 1/2 or B −1/2 .
Theorem 4.2 (running time of LazyEV). Let A, B ∈ Rd×d be two symmetric matrices satisfying B  0 and −B  A  B. Suppose M = B −1/2 AB −1/2 and RanInit(d) is the random vector generator defined in Proposition 3.3 with respect to B. Then, the computation of
← B −1/2 LazyEV(A, M, k, δ× , εpca , p) can be implemented to run in time


2
e knnz(B)+k
√ d+kΥ where Υ is the time to multiply B −1 A to a vector,7 or
• O
δ×
 √

2d
e k κB nnz(B)+knnz(A)+k
√
• O
if we use Conjugate gradient to multiply B −1 A to a vector.

V

δ×

e notation hides polylogarithmic factors with respect to 1/εpca , 1/δ× , 1/p, κB , d.
Above, the O

Our main theorems immediately imply the following corollaries (proved in Appendix E.2):

6
The complete specifications of εpca is included in Appendix E. Since our final running time only depends on
log(1/εpca ), we have not attempted to improve the constants in this polynomial dependency.
7
e
More precisely, to compute (B −1 A)w for some vector w with error ε where log(1/ε) = O(1).

10

Theorem 4.3 (gap-dependent k-GenEV). Let A, B ∈ Rd×d be two symmetric matrices satisfying
B  0 and −B  A  B. Suppose the generalized eigenvalue and eigenvector pairs of A with
|λ |−|λ
|
respect to B are {(λi , ui )}di=1 , and it satisfies 1 ≥ |λ1 | ≥ · · · ≥ |λd |. Let gap = k |λk |k+1 ∈ [0, 1]
be the relative eigengap. For fixed ε, p > 0, consider the output


ε4 ·gap 
−1/2 LazyEV A, B −1/2 AB −1/2 , k, gap, O
←
B
,
p
∈ Rd×k .
k
k3 (σ1 /σ )4

V

Then, defining

k

W = (uk+1, . . . , ud), we have with probability at least 1 − p:
Vk>B Vk = I and kVk>BW
k2 ≤ ε .


e
Moreover, our running time is O

√
k κB nnz(B)+knnz(A)+k2 d
√
gap

.

Theorem 4.4 (gap-free k-GenEV). In the same setting as Theorem 4.3, for ε, p > 0, consider


 
5
def
( 1 , . . . , k ) = k ← B −1/2 LazyEV A, B −1/2 AB −1/2 , k, ε, O k3 d4 (σ1ε/σ )12 , p .

v

v

V

V V
>
kB

Then, with probability at least 1 − p: it satisfies
k = I and
h
i
|λs |
∀s ∈ [k] : s> A s ∈ (1 − ε)|λs |,
and
max
1−ε
w∈Rd ∧w> BVk =0

 √
2d
e k κB nnz(B)+knnz(A)+k
√
.
Moreover, our running time is O

v v

k+1

w > Aw
w> B w

≤

|λk+1 |
.
1−ε

ε

5

Ideas Behind Theorems 4.1 and 4.2

Our LazyEV algorithm reduces the problem of finding generalized eigenvectors to finding regular
eigenvectors of M = B −1/2 AB −1/2 . In Section 5.1 we discuss how to ensure accuracy: that is, why
does LazyEV guarantee to find approximately the top absolute eigenvectors of M ; and in Section 5.2
we discuss how to implement LazyEV without ever needing to compute B 1/2 or B −1/2 .

5.1

Ideas Behind Theorem 4.1: Approximation Guarantee of GenEV

Our approximation guarantee in Theorem 4.1 is a natural generalization of the recent work on
fast iterative methods to find the top k eigenvectors of a PSD matrix M [2]. That method is
called LazySVD. At a high level, LazySVD finds the top k eigenvectors of M one by one but only
approximately. Starting with M0 = M , in the s-th iteration where s ∈ [k], LazySVD computes
approximately the leading eigenvector of matrix Ms−1 (using shift-and-invert [13]) and call it vs .
Then, LazySVD projects Ms ← (I − vs vs> )Ms−1 (I − vs vs> ) and proceeds to the next iteration.
While the algorithmic idea of LazySVD is simple, the analysis requires some careful linear algebraic lemmas. Most notably, if vs is an approximate leading eigenvector of Ms−1 , then one needs to
prove that the small eigenvectors of Ms−1 somehow still “embed” into that of Ms after projection.
This is achieved by a gap-free variant of the Wedin theorem plus a few other technical lemmas, and
we recommend interested readers to see the high-level overview section of [2].
In this paper, to relax the assumption that M is PSD, and to find leading eigenvectors whose
absolute eigenvalues are large, we have to make some non-trivial changes in the algorithm and
the analysis. On the algorithm side, LazyEV replaces the use of the shift-and-invert protocol in
LazySVD with our two-sided variant developed in Section 3. On the analysis side, we have to make
sure all lemmas properly deal with negative eigenvalues: for instance, if we perform a projection
M 0 ← (I−vv > )M (I−vv > ) where v correlates by at most ε with all eigenvectors of M whose absolute
eigenvalues are smaller than a threshold µ, then, after the projection, we need to prove that these
eigenvectors can be approximately “embedded” into the eigenspace spanned by all eigenvectors of
11

M 0 whose absolute eigenvalues are smaller than µ+τ . The approximation of this embedding should
depend on ε, µ and τ . See Lemma C.4 in the appendix.
The proof of Theorem 4.1 is in Appendix E, and the matrix algebraic lemmas are in Appendix C.

5.2

Proof of Theorem 4.2: Fast Implementation of GenEV

We can implement LazyEV efficiently without the necessity of computing B 1/2 or B −1/2 . In each
iteration of LazyEV, we call AppxPCA± and compute a vector vs0 . We do not explicitly store vs0 ,
but rather write it as vs0 = B 1/2 s0 and store only s0 ∈ Rd . We shall later ensure that AppxPCA±
outputs s0 directly. Similarly, we also write vs = B 1/2 s and only store s . All together, we do
not explicitly compute Vs , but instead write Vs = B 1/2 s and only keep track of s ∈ Rd×s .

v

v

v

v
V

v

V

V

Now, the computation of vs becomes the B-projection into the s−1 space:


>
>
0
>
0 >
B
k(I − Vs−1 Vs−1
)vs0 k = kB 1/2 s0 − B 1/2 s−1 s−1
B s0 k =
s − s−1 s−1 B s


>
v
vs0 − Vs−1Vs−1
B vs0
>
>
>
>
k(I − Vs−1 Vs−1
)vs0 k · vs = B −1/2 (I − Vs−1 Vs−1
)vs0 = B −1/2 (I − Vs−1 Vs−1
)B 1/2 vs0 = vs0 − Vs−1 Vs−1
B vs0

v

V V

v V V

v

and this can be implemented to run in O(kd + nnz(B)) time. Finally, we write
Ms = (I − Vs Vs> )B −1/2 AB −1/2 (I − Vs Vs> ) = B −1/2 (I − B

VsVs>)A(I − VsVs>B)B−1/2

and only pass it implicity to AppxPCA± (without directly computing this matrix).

w

w

To implement AppxPCA± , we again write all vectors w
bt = B 1/2 t and only store t . Thus,
> B
1/2
the normalization wa ← w
bm1 /kw
bm1 k2 becomes the B-normalization a ← m1 /( m
m1 )
1
±
which runs in O(nnz(B)) time. Recall that AppxPCA makes a polylogarithmic number of calls to
the matrix inversion subroutine A, each time requesting to approximately invert either λI − Ms or
λI + Ms . Let us only focus on inverting λI − Ms and the other case is similar. We write

def
= B −1/2 λI − Ms B 1/2 = λI − (I − s s> B)B −1 A(I − s s> B) .

N

VV

Now, the accuracy requirement in

AppxPCA±

w

w

w

w

VV

becomes

wt satisfying kB1/2wt − (λI − Ms)−1B1/2wt−1k ≤ εe
find wt satisfying kB 1/2 wt − B 1/2 N−1 wt−1 k ≤ εe
p
find wt satisfying kwt − N−1 wt−1 k ≤ εe/ λmax (B)
find

⇐⇒
⇐=

Using accelerated gradient descent (see Theorem D.1), we can reduce this approximate inversion
−1
0
t ←
t−1 to T times of approximate matrix-vector
p multiplication (i.e., w ← w) for T =
√
8
e κ(λI−M ) ). We can further derive that T = O(1/
e
O(
δ× ) owing to Theorem 3.1. Notice that
s
0
each time we compute w ← w it suffices to compute it to an additive accuracy kw0 − wk ≤ ε
e
where the error satisfies log(1/ε) = O(1).

w

N w

N

N

N

N

VV

VV

Finally, the matrix-vector multiplication w = λw−(I − s s> B)B −1 A(I − s s> B)w consists
of two rank-s B-projections which run in time O(nnz(B) + kd), plus the time needed to multiply
B −1 A to a vector. This proves that LazyEV can be implemented so that
p
e
δ× ) times,
• It computes matrix-vector multiplication of the form w0 ← B −1 Aw a total of O(k/
e
each time to an accuracy ε where log(1/ε) = O(1);
p 
e (knnz(B) + k 2 d)/ δ× time.
• The rest of the computation costs a total of O
This finishes the proof of the first item of the theorem.
8

This reduction would be obvious if we required the matrix-vector multiplication to be exact, and for instance
Chebyshev method serves for exactly this purpose. However, in order to relax the multiplication to be approximate,
and without incurring an error that blows up exponentially with T , we build our own inexact variant of the accelerate
gradient descent method AGDinexact in Appendix D that could be of independent interest.

12

1/2

Algorithm 3 LazyCCA(A, M, k, δ× , εpca , p)

Input: A, an approximate matrix inversion method;
M ∈ Rd×d , a matrix satisfying −I  M  I;
k ∈ [d], the desired rank;
δ× ∈ (0, 1), a multiplicative error;
εpca ∈ (0, 1), a numerical accuracy parameter; and
p ∈ (0, 1), a confidence parameter.
1: M0 ← M ;
2: V0 = [];
3: for s = 1 to k do
4:
vs0 ← AppxPCA± (A, Ms−1 , δ× /2, εpca , p/k);
>
0
> )v 0 / (I − V
⊥
 project vs0 to Vs−1
5:
vs ← (I − Vs−1 Vs−1
s−1 Vs−1 )vs ;
s
6:
write vs =√(ξs0 , ζs0 ) where ξs0 ∈ Rdx √and ζs0 ∈ Rdy ;
7:
ξs ← ξs0 /( 2kξs0 k2 ) and ζs ← ζs0 /( 2kζs0 k2 );
h
 ξ −ξ i
s
s
;
8:
Vs ← Vs−1 ,
ζs ζs


9:
Ms ← I − 2diag(ξs ξs> , ζs ζs> ) Ms−1 I − 2diag(ξs ξs> , ζs ζs> )
 or equivalently, Ms = (I − Vs Vs> )M (I − Vs Vs> )
10: end for
11: return Vk .
As for the second item, we simply notice that whenever we want to compute w0 ← B −1 Aw, we
−1
can first compute Aw in time O(nnz(A)), and then use Conjugate gradient [27] to compute
 B
√
e κB · nnz(B) where
applied to this vector. The running time of Conjugate gradient is at most O
e
the O factor hides a logarithmic factor on the accuracy.


6

LazyCCA: Our Algorithm for Canonical Correlation Analysis

We propose LazyCCA (see Algorithm 3), a variant of LazyEV that is specially designed for matrices
M of the form M = B −1/2 AB −1/2 , where A and B come from a CCA instance following Lemma 2.3.
More specifically, recall from Lemma 2.3 that the eigenvectors of matrices M arising from CCA
instances are symmetric: if (ξ, ζ) is a normalized eigenvector of M with eigenvalue σ where ξ ∈
Rdx , ζ ∈ Rdy , then (−ξ, ζ) is also a normalized eigenvector but with √
eigenvalue −σ. Furthermore,
since (ξ, ζ) is orthogonal to (−ξ, ζ), we must have kξk = kζk = 1/ 2. Our LazyCCA method is
designed to ensure such symmetry and orthogonality as well. When an approximate eigenvector
√
vs = (ξs0 , ζs0 ) is obtained, we re-scale the pair to ξs and ζs where both of them have norm exactly 1/ 2
(see Line 7 of LazyCCA). Then, we simultaneously add two (orthogonal) approximate eigenvectors
(ξs , ζs ) and (−ξs , ζs ) to the column orthonormal matrix Vs .
Remark 6.1. This re-scaling step, together with the fact that we find vector pairs one by one,
allows us to provide per-vector guarantee on the obtained approximate correlation vectors (see
Theorem 6.3). This is in contrast to CCALin which is based on subspace power method so can
only find the subspace spanned by the top k correlation vectors but not distinguish them.
We prove three main theorems in the appendix:
• A main theorem for the approximation guarantee (see Theorem F.1).

The main “delta” between the proofs of Theorem F.1 and Theorem 4.1 is to show that, after
13

re-scaling, the vector (ξs , ζs ) is also an approximate leading eigenvector of Ms−1 . In particular,
its Rayleigh quotient can only become better after scaling (see (F.3) in the appendix).
• A main theorem for the stochastic running time (see Theorem F.2).

The main “delta” between the proofs of Theorem F.2 and Theorem 4.2 is to replace the use of
conjugate gradient with the Katyusha method to multiply B −1 A to a vector (see Lemma 2.6).

• A main theorem for the doubly-stochastic running time guarantee (see Theorem F.3).

The main “delta” between the proofs of Theorem F.3 and Theorem F.2 is to use accelerated
SVRG as opposed to Katyusha to compute matrix inverse.

We state below the final statements on LazyCCA.
Theorem 6.2 (gap-dependent k-CCA). Let X ∈ Rn×dx , Y ∈ Rn×dy be two matrices with
canonical-correlation coefficients 1 ≥ σ1 ≥ · · · σr ≥ 0 and the corresponding correlation vecσ −σ
> , 0]]
tors {(φi , ψi )}ri=1 . Let gap = k σkk+1 ∈ [0, 1] be the relative gap, and define A = [[0, Sxy ]; [Sxy
and B = diag(Sxx , Syy ) following Def. 2.2. For every ε, p > 0, consider the output
 ±φ0 . . . ±φ0 


√ −1/2
def
ε4 ·gap 
−1/2 AB −1/2 , ε, gap, O
1
k
=
←
2B
LazyCCA
A,
B
,
p
.
k
k3 (σ1 /σk )4
ψ10 . . . ψk0

V

Then, letting

Vφ = (φ01, . . . , φ0k ), Vψ = (ψ10 , . . . , ψk0 ), Wφ = (φk+1, φk+2, . . . ),

and

Wψ = (ψk+1, ψk+2, . . . ) ,

we have with probability at least 1 − p:

Vφ ∈ Rd ×k satisfies Vφ>SxxVφ = I and kVφ>SxxWφk2 ≤ ε ,
Vψ ∈ Rd ×k satisfies Vψ>Syy Vψ = I and kVψ>Syy Wψ k2 ≤ ε .
x

y

√


knnz(X,Y )·(1+ κ0 /n)+k2 d
e
√
The (stochastic) running time is O
. The doubly-stochastic running time
gap
√



√


1/4
0
κ0
e nnz(X, Y ) · 1 + √κ /n
e knnz(X, Y ) · 1 + √
for
is O
for
1-CCA,
and
O
gap·σ1
gap·σ ·(nnz(X,Y )/kd)1/4
k

k-CCA as long as nnz(X, Y ) ≥ kd.

Theorem 6.3 (gap-free k-CCA). In the same setting as Theorem 6.2, for ε, p > 0, consider
 ±φ0 . . . ±φ0 


√
ε4 ·gap 
1
k
= k ← 2B −1/2 LazyCCA A, B −1/2 AB −1/2 , ε, gap, O k3 (σ
,
p
.
4
0
0
1 /σk )
ψ1 . . . ψk

V

Vφ = (φ01, . . . , φ0k ) ∈ Rd ×k and Vψ = (ψ10 , . . . , ψk0 ) ∈ Rd ×k , with probability at least 1 − p,
• Vφ> Sxx Vφ = I, Vψ> Syy Vψ = I;

Letting

y

x

• (1 − ε)σi ≤ |φ0i Sxy
n ψi | ≤ (1 + ε)σi for every i ∈ [k]; and

o
• maxφ∈Rdx ,ψ∈Rdy φ> Sxy ψ φ> Sxx φ = 0 ∧ ψ > Syy ψ = 0 ≤ (1 + ε)σk+1 .
√


knnz(X,Y )·(1+ κ0 /n)+k2 d
e
√
The (stochastic) running time is O
. The doubly-stochastic running time
ε
√



√


1/4
0
κ /n
κ0
e nnz(X, Y ) · 1 + √
e knnz(X, Y ) · 1 + √
is O
for
1-CCA,
and
O
for
ε·σ1
ε·σ ·(nnz(X,Y )/kd)1/4

V

V

k

k-CCA as long as nnz(X, Y ) ≥ kd.

Remark 6.4. The doubly-stochastic running times for the general k > 1 setting can be slightly
faster (but notationally more involved) than the ones stated above. For instance, we are aware of
a proof that gives running time

√

1/4 1/4 

max k min{σ1 , σ1 κ/d}, 1
d
κ0
e
O knnz(X, Y ) · 1 + √
·
1/4
gap · σk
(nnz(X, Y ))
14

for the gap-dependent case of k-CCA, or the same formula but replacing gap with ε for the gapfree case. We refrain from proving it in full because as discussed in Section 1.2, doubly-stochastic
running times are not necessarily interesting.

Appendix
A

Proof of Lemma 2.6
def

Proof of Lemma 2.6. First of all, computing B −1 Aw is equivalent to minimizing f (x) = 21 x> Bx −
x> Aw. Suppose we write x = (x1 , x2 ) and w = (w1 , w2 ) where x1 , w1 ∈ Rdx , x2 , w2 ∈ Rdy , then
1
one can rewrite f as f (x) = 2n
kXx1 − Y w2 k22 + kXw1 − Y x2 k22 + C where C is a fixed constant.
Therefore, one can also write
n
1 X
f (x) =
(hXi , x1 i − hYi , w2 i)2 + (hYi , x2 i − hXi , w1 i)2
2n
i=1

where each Xi is a row vector of X and Yi is a row vector of Y . In such a case, we observe that

• f (x) is an average of 2n smooth functions, where function (hXi , x1 i − hYi , w2 i)2 is smooth
with parameter 2kXi k2 (meaning Hessian bounded by kXi k2 in spectral norm) and (hYi , x2 i −
hXi , w1 i)2 is smooth with parameter 2kYi k2 . In other words, each function has a smoothness
parameter at most 2 maxi {kXi k2 , kYi k2 }.

• f (x) is at least λmin (B) = min{λmin (Sxx ), λmin (Syy )} strongly convex, meaning the Hessian
∇2 f (x) has no eigenvalue less than this quantity.

For such reason, one can apply the convergence theorem of Katyusha [1] to find some x such that
p

0
(x∗ ) 
f (x) − miny f (y) ≤ εe in time O nnz(X, Y ) · 1 + κ0 /n log f (x )−f
where x0 is an arbitrary
εe
starting vector fed into Katyusha and x∗ is the exact minimizer. If we choose x0 to be the zero
vector, it is easy to verify that f (x0 ) − f (x∗ ) ≤ O(λmax (B) · kwk2 ).

It is not hard to see that an additive εe minimizer of f (x) implies an ε-approximate solution for
the inverse kw0 − B −1 Awk ≤ ε where ε2 = 2e
ε/λmin (B). This finishes the proof of the running time
in Lemma 2.6.


B
B.1

Proof for Section 3: Two-Sided Shift-and-Invert
Inexact Power Method

In this subsection we review some classical convergence lemmas regarding power method and its
inexact variant. These lemmas almost directly follow from previous results such as [13, 15], and
are more similar to [2]. We skip the proofs in this paper.
Consider power method that starts with a random unit vector w0 ← RanInit(d) and apply
wt ← M wt−1 /kM wt−1 k iteratively.
Lemma B.1 (Exact Power Method). Let M be a PSD matrix with eigenvalues λ1 ≥ · · · ≥ λd
and the correpsonding eigenvectors u1 , . . . , ud . Fix an error tolerance ε > 0, parameter κ ≥ 1, and
failure probability p > 0, define
lκ
 9dθ m
T PM (κ, ε, p) =
log 2
2
p ε

15

Then, with probability at least 1 − p it holds that ∀t ≥ T PM (κ, ε, p):
X
(wt> ui )2 ≤ ε and wt> M wt ≥ (1 − 1/κ − ε)λ1 .
i∈[d],λi ≤(1−1/κ)λ1

Lemma B.2 (Lemma 4.1 of [13]). Let M be a PSD matrix with eigenvalues λ1 ≥ · · · λd . Fix an
accuracy parameter εe > 0, and consider two update sequences
w
b0∗ = w0 ,

w
b0 = w0 ,

Then, defining wt = w
bt /kw
bt k and
where

2
Γ(M, t) = t
λd
def



∗
∀t ≥ 1 : w
bt∗ ← M w
bt−1

∀t ≥ 1 : w
bt satisfies kw
bt − M w
bt−1 k ≤ εe,

wt∗

=w
bt∗ /kw
bt∗ k, it satisfies

kwt − wt∗ k ≤ εe · Γ(M, t),

t,
if λ1 = 1;
(λt1 − 1)/(λ1 − 1), if λ1 6= 1.

and we have Γ(M, t) ≤ 2t ·

max{1, λt1 }
λtd

Theorem B.3 (Inexact Power Method). Let M be a PSD matrix with eigenvalues λ1 ≥ · · · ≥ λd
and the corresponding eigenvectors u1 , . . . , ud . With probability at least 1 − p it holds that, for every
ε ∈ (0, 1) and every t ≥ T PM (κ, ε/4, p), if wt is generated by the power method with per-iteration
ε
error εe = 4Γ(M,t)
, then
X
(wt> ui )2 ≤ ε and wt> M wt ≥ (1 − 1/κ − ε)λ1 .
i∈[d],λi ≤(1−1/κ)λ1

B.2

Proof of Theorem 3.1

We prove Theorem 3.1 by first showing the following lemma. Most of these properties are analogous
to their original variants in [13, 14], but here we take extra care also on negative eigenvalues and
thus allowing M to be non-PSD.
Lemma B.4 (useful properties of AppxPCA± ).
letting λ∗ = kM k2 ):

With probability at least 1 − p, it holds that (by

1
1
and εe1 ≤ 32Γ((λ(s−1) I+M
for each iteration s ≥ 1;
32Γ((λ(s−1) I−M )−1 ,m1 )
)−1 ,m1 )
ε
ε
and εe2 ≤ 4Γ((λ(f ) I+M )−1 ,m ) when the repeat-until loop is over;
4Γ((λ(f ) I−M )−1 ,m2 )
2
3 (s−1)
∗ ) ≤ ∆(s) ≤ λ(s−1) − λ∗ and 1 (λ(s−1) − λ∗ ) ≤ λ(s) − λ∗ for each iteration
(λ
−
λ
4
2

(a) εe1 ≤
(b) εe2 ≤

(c) 0 ≤
and
× (f ) δ× ∗
(d) λ(f ) − λ∗ ∈ [ δ48
λ , 13 λ ] when the repeat-until loop is over.
(e) when the repeat-until loop is over,
10
if wa> va ≥ wb> vb then λ(f ) − λmax (M ) ≤ (λ(f ) − λ∗ ); or
3
10
if wa> va ≤ wb> vb then λ(f ) + λmin (M ) ≤ (λ(f ) − λ∗ ) .
3
def

s ≥ 1;

def

Proof. We denote by C (s) = (λ(s) I − M )−1 and by D(s) = (λ(s) I + M )−1 for notational simplicity.
Below we prove all the items by induction for a specific iteration s ≥ 2 assuming that the items of
the previous s − 1 iterations are true. The base case of s = 1 can be verified similar to the general
arguments after some notational changes. We omitted the proofs of the base case s = 1.
(a) Recall that
Γ(C (s−1) , t) ≤ 2t ·

max{1, λmax (C (s−1) )t }
λmin (C (s−1) )t

and
16

Γ(D(s−1) , t) ≤ 2t ·

max{1, λmax (D(s−1) )t }
λmin (D(s−1) )t

On one hand, we have λmax (C (s−1) ) =

1
λ(s−1) −λ∗

≤

2
λ(s−2) −λ∗

≤

2
∆(s−1)

using Lemma B.4.c of

× (s−1)
the previous iteration. Combining this with the termination criterion ∆(s−1) ≥ δ12
λ
, we
1
24
(s−1)
(s−1)
≥
) = λ(s−1) −λ
have λmax (C
) ≤ δ λ(s−1) . On the other hand, we have λmin (C
(M )

1

×

1

min

(s−1) , t) ≤ 2t(48/δ )t . It
×
2λ(s−1) . Combining the two bounds we conclude that Γ(C
δ×  m 1
1
1
.
is now obvious that εe1 ≤ 32Γ(C (s−1) ,m ) is satisfied because εe1 = 64m1 48
λ(s−1) +λ∗

≥

1

Similarly, on one hand, we have λmax (D(s−1) ) =
2

∆(s−1)

1
λ(s−1) +λmin (M )

≤

1
λ(s−1) −λ∗

≤

2
λ(s−2) −λ∗

≤

using Lemma B.4.c of the previous iteration. Combining this with the termination

δ× (s−1)
, we have λmax (D(s−1) ) ≤ δ λ24
(s−1) . On the other hand, we have
12 λ
×
1
≥ 2λ 1 . Combining the two bounds we conclude
λmin (D(s−1) ) = λ(s−1) +λ1 (M ) ≥ λ(s−1)
+λ∗
(s−1)
max
1
is satisfied.
that Γ(D(s−1) , t) ≤ 2t(48/δ× )t . It is now obvious that εe1 ≤ 32Γ(D(s−1)
,m )

criterion ∆(s−1) ≥

1

Γ(C (f ) , t)

(b) The same analysis as in the proof of Lemma B.4.a suggests that
≤ 2t(48/δ× )t and
ε
(f
)
t
Γ(D , t) ≤ 2t(48/δ× ) . These immediately imply εe2 ≤ 4Γ(C (f ) ,m ) and εe2 ≤ 4Γ(D(fε ) ,m )
2
2
δ× m2
ε
because εe2 = 8m
2 48
(c) Because Lemma B.4.a holds for the current iteration s we can apply Theorem B.3 (with ε =
1/16 and κ = 16) and get
7
7
wa> C (s−1) wa ≥ λmax (C (s−1) ) and wb> D(s−1) wb ≥ λmax (D(s−1) ) .
8
8
±
By the definition of v in AppxPCA and the Cauchy-Schwartz inequality it holds that
 

wa> va = wa> C (s−1) wa + wa> va − C (s−1) wa ∈ wa> C (s−1) wa − εe1 , wa> C (s−1) wa + εe1 , and
 

wb> vb = wb> D(s−1) wb + wb> vb − D(s−1) wb ∈ wb> D(s−1) wb − εe1 , wb> D(s−1) wb + εe1 .

Combining the above equations we have
i
h7
ε1 , λmax (C (s−1) )
wa> va − εe1 ∈ λmax (C (s−1) ) − 2e
8
h3
i 3 
1
⊆ λmax (C (s−1) ), λmax (C (s−1) ) = , 1 · (s−1)
, and
4
4
λ
− λmax (M )
h3
i 3 
1
.
(B.1)
wb> vb − εe1 ⊆ λmax (D(s−1) ), λmax (D(s−1) ) = , 1 · (s−1)
4
4
λ
+ λmin (M )


def
In other words, ∆(s) = 34 · max{w> v 1,w> w }−eε ∈ 43 (λ(s−1) − λ∗ ), λ(s−1) − λ∗ because λ∗ =
a

max{λmax (M ), −λmin (M )}.

a

b

b

1

At the same time, our update rule λ(s) = λ(s−1) − ∆(s) /2 ensures that λ(s) − λ∗ = λ(s−1) −
(s−1)
∗
λ∗ − ∆(s) /2 ≥ λ(s−1) − λ∗ − λ 2 −λ = 12 (λ(s−1) − λ∗ ).

(f )
(f )
(d) The upper bound holds because λ(f ) − λ∗ = λ(f −1) − ∆2 − λ∗ ≤ 43 − 12 ∆(f ) ≤ 5δ×72λ where
the first inequality follows from Lemma B.4.c of this last iteration, and the second inequality
λ(f )
. Simply rewriting this inequality we have
follows from our termination criterion ∆(f ) ≤ δ×12
5δ× /72
δ× ∗
(f
)
∗
∗
λ − λ ≤ 1−5δ× /72 λ < 13 λ .
The lower bound is because using Lemma B.4.c (of this and the previous iteration) we have

(f −1) ¬ δ× λ(f −1)
λ(f )
λ(f ) − λ∗ ≥ 41 λ(f −2) − λ∗ ≥ ∆ 4
≥
≥ δ×48
. Here, inequality ¬ is because
48
∆(f −1) >

δ× λ(f −1)
12

due to the termination criterion.

17

(e) We only prove the case when wa> va ≥ wb> vb and the other case is similar. We compute that

∆(f )
∆(f ) ¬ 4 (f −1)
≤ (λ
+ λmin (M )) −
2
3
2
(f
)
­
4
∆
4
δ× λ(f )
= (λ(f ) + λmin (M )) +
≤ (λ(f ) + λmin (M )) +
3
2
3
24
® 4
¯ 10
(f )
(f )
∗
(f )
∗
≤ (λ + λmin (M )) + 2(λ − λ ) ≤ (λ − λ ) .
3
3
>
>
Above, ¬ is from (B.1) together with the fact that wa va ≥ wb vb ; ­ is using the termination
λ(f ) − λmax (M ) = λ(f −1) − λmax (M ) −

criterion ∆(f ) ≤

δ× λ(f )
12 ;

® is from Lemma B.4.d

Finally, since the success of Theorem B.3 only depends on the randomness of w
b0 , we have that with
probability at least 1 − p all the above items are satisfied.

We are now ready to prove Theorem 3.1.

Proof of Theorem 3.1. We only focus on the case when sgn = + and the other case is similar. It
follows from Theorem B.3 (with κ = 2) that, letting µi = 1/(λ(f ) −λi ) be the i-th largest eigenvalue
of the matrix (λ(f ) I − M )−1 , then
X
(w> ui )2 ≤ ε .
i∈[d],µi ≤µ1 /2

(f ) − λ∗ )
Note that if an index i ∈ [d] satisfies λ∗ − λi ≥ δ2× λ∗ , then we must have λ∗ − λi ≥ 13
2 (λ
δ
× ∗
(f ) − λ∗ ).
owing to λ(f ) − λ∗ ≤ 13
λ from Lemma B.4.d. This further implies that λ(f ) − λi ≥ 15
2 (λ
15
3
(f
)
(f
)
(f
)
Plugging in Lemma B.4.e we further have λ − λi ≥ 2 · 10 (λ − λ1 ) > 2(λ − λ1 ). Using the
definition of µi , we must have µ1 /2 > µi . In sum, we also have
X
(w> ui )2 ≤ ε .
i∈[d],λi ≤(1−δ× /2)λ∗

On the other hand,
>

w Mw =

d
X
i=1

λi (w> ui )2 ≥ −ελ∗ +

X

λi (w> ui )2

i∈[d],λi >(1−δ× /2)λ∗

≥ −ελ∗ + (1 − δ× /2)λ∗ ·
∗

X

(w> ui )2

i∈[d],λi >(1−δ× /2)λ∗

≥ −ελ + (1 − δ× /2)(1 − ε)λ∗ ≥ (1 − δ× /2)(1 − 3ε)λ∗ .
The number of oracle calls to A is determined by the number of iterations in the repeat-until
loop. It is easy to verify that there are at most O(log(1/δ× )) such iteartions, so the total number
of oracle calls to A is only O(log(1/δ× )m1 + m2 ).
As for the condition number, each time we call A we have

λ(s) − λd
λ(s) + λ∗
2λ(s)
λmax (λ(s) I − M )
≤
≤
≤
λmin (λ(s) I − M )
λ(s) − λ1
λ(s) − λ∗
λ(s) − λ∗
If s = 0 then we have
λ(s)
δ× λ(s+1) /12

≤

12
δ×

λ(0)
λ(0) −λ∗

≤

1+δ×
δ×

and

λmax (λ(s) I + M )
2λ(s)
≤
λmin (λ(s) I + M )
λ(s) − λ∗

because λ∗ ≤ 1. If s ≤ f − 2 then we have

λ(s)
λ(s) −λ∗

≤

λ(s)
∆(s+1)

≤

where the first inequality follows from Lemma B.4.c, the second inequality follows

from the stopping criterion, and the third inequality follows from the monotonicity of λ(s) . If
(s)
(s)
λ(s)
2λ(s)
s = f − 1 then we have λ(s)
≤ λ(s−1)
≤ 2λ
≤ δ 2λ
= δ24
where the first two inequalities
−λ∗
−λ∗
∆(s)
λ(s) /12
×
×

18

follow from Lemma B.4.c and the third inequality follows from our stopping criterion. If s = f
(λ(s) I−M )
λ(s)
96
≤ δ48
owing to Lemma B.4.d. In all cases we have λλmax(λ
then we have λ(s)
(s) I−M ) ≤ δ× and
−λ∗
×
min

λmax (λ(s) I+M )
λmin (λ(s) I+M )

≤

96
δ× .

(s)
1
= λ (λλ(s) I−M )
λmin (λ(s) I−M )
min
from λ(s) ≥ λ∗ . Similarly, we also

Finally, we have
inequality follows

C

·

1
λ(s)

have

λ(s)
· λ1∗ ≤ δ×48λ∗
λ(s) −λ∗
1
≤ δ×48λ∗ .
λmin (λ(s) I+M )

≤

where the first



Main Matrix-Algebra Lemmas

In this section we provide some necessary lemmas on matrix algebra that shall become essential
for our proof of Theorem 4.1. Many of these lemmas are analogous to those ones used in the SVD
algorithm by the same authors of this paper [2], however, we need some extra care in this paper
because the underlying matrix M is no longer PSD.
Proposition C.1. Let A, B be two (column) orthonormal matrix such that for η ≥ 0,
A> BB > A  (1 − η)I

Then we have: there exists a matrix Q, kQk2 ≤ 1 such that
√
kA − BQk2 ≤ η

Proof. Since A> A = I and A> BB > A  (1 − η)I, we know that A> B ⊥ (B ⊥ )> A  ηI. By the fact
that
A = (BB > + B ⊥ (B ⊥ )> )A = BB > A + B ⊥ (B ⊥ )> A
we can let Q = B > A and obtain
kA − BQk2 ≤ kB ⊥ (B ⊥ )> Ak2 ≤

C.1

√



η .

Approximate Projection Lemma

The next lemma states that, projecting a symmetric matrix M into the orthogonal space of Vs ∈
Rd×s is almost equivalent to projecting it into the orthogonal space of Qs ∈ Rd×s , if Qs is the
projection of Vs into the orthogonal space of U but kVs> U k is small. This lemma is obvious if
“small” means zero correlation: if Vs were completely orthogonal to U then Qs would equal to Vs ,
so projecting M into the orthogonal space of Vs would be equivalent to that of Qs . However, even
in the inexact scenario, this argument is true.
Lemma C.2. Let M be a symmetric matrix with (not necessarily sorted) eigenvalues λ1 , . . . , λd
and the corresponding (normalized) eigenvectors u1 , . . . , ud ∈ Rd . For every k ≥ 1, define U ⊥ =
(u1 , . . . , uk ) ∈ Rd×k and U = (uk+1 , . . . , ud ) ∈ Rd×(d−k) . For every ε ∈ (0, 21 ), let Vs ∈ Rd×s be a
column orthogonal matrix such that kVs> U k2 ≤ ε, define Qs ∈ Rd×s to be an arbitrary orthogonal
basis of the column span of U ⊥ (U ⊥ )> Vs , then we have:

 
 
 

>
>
>
≤ 13εkM k2 .
I − Qs Q>
M
I
−
Q
Q
−
I
−
V
V
M
I
−
V
V
s s
s s
s s
s
2

Proof of Lemma C.2. Since Qs is an orthogonal basis of the column span of U ⊥ (U ⊥ )> Vs , there is
a matrix R ∈ Rs×s such that
Qs = U ⊥ (U ⊥ )> Vs R
Using the fact that Q>
s Qs = I, we have:
(U ⊥ (U ⊥ )> Vs R)> (U ⊥ (U ⊥ )> Vs R) = I =⇒ R> Vs> U ⊥ (U ⊥ )> Vs R = I .
19

By the fact that Vs> Vs = I and U ⊥ (U ⊥ )> + U U > = I, we can rewrite the above equality as:


R> I − Vs> U U > Vs R = I
(C.1)

From our lemma assumption, we have: kVs> U k2 ≤ ε, which implies 0  Vs> U U > Vs  ε2 I.
Putting this into (C.1), we obtain:
1
4 
I  R> R 
I  1 + ε2 I
2
1−ε
3

>
The above inequality directly implies that I  RR  1 + 43 ε2 I. Therefore,
>
Qs Q>
s − Vs Vs

=
=
≤
≤

2

⊥

⊥ >

>

U (U ) Vs RR Vs> U ⊥ (U ⊥ )> − Vs Vs>
⊥

⊥ >

>

⊥

⊥ >

U (U ) Vs RR

Vs> U ⊥ (U ⊥ )>

>

U (U ) Vs (RR −
RR> − I

2

− (U (U ⊥ )> + U U > )Vs Vs> (U ⊥ (U ⊥ )> + U U > )

I)Vs> U ⊥ (U ⊥ )>

+ U > Vs Vs> U

2

2

⊥

2

+ UU

+ 2 Vs> U U > Vs

>

Vs Vs> U U >

2

⊥

2

⊥ >

+ 2 U (U ) Vs Vs> U U >

2

4 2
19
≤
ε + ε2 + 2ε < ε .
3
6
Finally, we have

 
 
 

>
>
>
I − Qs Q>
M
I
−
Q
Q
−
I
−
V
V
M
I
−
V
V
s s
s s
s s
s




2

>
>
>
>
>
>
≤ 2 Qs Qs − Vs Vs M + Qs Qs − Vs Vs M Qs Qs
+ Qs Q>
−
V
V
M Vs Vs>
s s
s
2

2

19 × 4
εkM k2 < 13εkM k2 .
≤
6

C.2

2

1/2

2



Gap-Free Wedin Theorem

Lemma C.3 (two-sided gap-free Wedin theorem). For ε ≥ 0, let A, B be two symmetric matrices
such that kA − Bk2 ≤ ε. For every µ ≥ 0, τ > 0, let U be column orthonormal matrix consisting
of eigenvectors of A with absolute eigenvalues ≤ µ, let V be column orthonormal matrix consisting
of eigenvectors of B with absolute eigenvalues ≥ µ + τ , then we have:
ε
kU > V k ≤
.
τ
Proof of Lemma C.3. We write A and B in terms of eigenvalue decomposition:
e0 V 0> ,
e > + V 0Σ
A = U ΣU > + U 0 Σ0 U 0> and B = V ΣV

where U 0 is orthogonal to U and V 0 is orthogonal to V . Letting R = A − B, we obtain:
ΣU > = U > A = U > (B + R)

e + U > RV
=⇒ ΣU > V = U > BV + U > RV = U > V Σ
e −1 = U > V + U > RV Σ
e −1 .
=⇒ ΣU > V Σ

Taking spectral norm on both sides, we obtain:
e −1 k2 ≥ kΣU > V Σ
e −1 k2 ≥ kU > V k2 − kU > RV Σ
e −1 k2 .
kΣk2 kU > V k2 kΣ
This can be simplified to

ε
µ
kU > V k2 ≥ kU > V k2 −
,
µ+τ
µ+τ
20

and therefore we have kU > V k2 ≤

C.3

ε
τ



as desired.

Eigenvector Projection Lemma

Our next technical lemma studies the projection of a matrix M into the orthogonal direction of
a vector v, where v has little correlation with M ’s leading eigenvectors below some threshold µ
(denoted by U ). The conclusion of the lemma says that, after the projection, if we study the leading
eigenvectors of M 0 = (I − vv > )M (I − vv > ) below some threshold µ + τ and denote it by V1 , then
U approximately embeds into V1 , meaning that although V1 could be of a larger dimension of U ,
however, there exists a matrix Q with spectral norm no more than 1 such that kU − V1 Qk2 is small.
Lemma C.4. Let M ∈ Rd×d be a symmetric matrix with eigenvalues λ1 , . . . , λd and corresponding
eigenvectors u1 , . . . , ud . Suppose |λ1 | ≥ · · · ≥ |λd |. Define U = (uj+1 , . . . , ud ) ∈ Rd×(d−j) to be the
matrix consisting of all eigenvectors with absolute eigenvalues ≤ µ. Let v ∈ Rd be a unit vector
such that kv > U k2 ≤ ε ≤ 1/2, and define

 

0
>
>
M = I − vv M I − vv
.

Then, denoting by [V2 , V1 , v] ∈ Rd×d the unitary matrix consisting of (column) eigenvectors of
where V1 consists of eigenvectors with absolute eigenvalue ≤ µ + τ , then there exists a matrix
Q with spectral norm kQk2 ≤ 1 such that
r
169ε2 kM k22
kU − V1 Qk2 ≤
+ ε2 .
τ2

M 0,

⊥

⊥ >

(U ) v
Proof of Lemma C.4. Using Lemma C.2, let q = kUU⊥ (U
be the projection of v to U ⊥ , we
⊥ )> vk
2
know that

 
 

 
≤ 13εkM k2 .
I − qq > M I − qq > − I − vv > M I − vv >


>

Denote I − qq M I − qq
with eigenvalue λj+1 , . . . , λd .


>

2

as M 00 . We know that uj+1 , . . . , ud are still eigenvectors of M 00

Apply Lemma C.3 on A = M 00 , U and B = M 0 , V = V2 , we obtain:
13εkM k2
kU > V2 k2 ≤
.
τ
This implies that


169ε2 kM k22
>
>
>
>
> >
2
U V1 V1 U = I − U V2 V2 U − U vv U  1 −
−ε I ,
τ2
where the inequality uses the assumption kv > U k2 ≤ ε.

Apply Proposition C.1 to A = U and B = V1 , we conclude that there exists a matrix Q,
kQk2 ≤ 1 such that
r
169ε2 kM k22
kU − V1 Qk2 ≤
+ ε2 .

τ2

D

Matrix Inversion via Approx Accelerated Gradient Descent

N

Given a positive definite matrix , it is well-known that one can reduce the (approximate) matrix
inversion problem −1 χ to multiple computations of the matrix-vector multiplication (i.e., of the
form w0 ← w). In particular, Chebyshev method [9] uses the so-called Chebyshev polynomial for

N

N

21

Algorithm 4 AGDinexact (f, x0 , T )
Input: f an L-smooth and σ-strongly convex function;
x0 some initial point; and
T the number of iterations.
Output: yT .
√
√2
1: τ ←
, η ← τ1L .
 τ = O( √Lσ ) and η = O( √1σL )
1+ 8L/σ+1
2: y0 ← x0 , z0 ← x0 .
3: for k ← 0 to T − 1 do
4:
xk+1 ← τ zk + (1 − τ )yk .
e (xk+1 ) satisfying k∇f
e (xk+1 ) − ∇f (xk+1 )k2 ≤ εe.
5:
Compute approximate gradient ∇f
1 e
6:
yk+1 ← xk+1 − L ∇f (xk+1 )

1
e (xk+1 )
7:
zk+1 ← 1+ησ
zk + ησxk+1 − η ∇f
8: end for
9: return yT .
this purpose, and the number of matrix-vector multiplications is determined by the degree of that
polynomial.
In this section, we revisit this problem by allowing matrix-vector multiplications to be computed
only approximately. We emphasize that this is not a simple task in general. If matrix inversion
is reduced to T matrix-vector multiplications, then a standard analysis implies that each of these
multiplications must be computed up to a very small error 2−Ω(T ) . If the actual matrix-vector
multiplication subroutine has a logarithmic dependency on the error in its running time, then we
will have a total running time at least quadratically dependent on T .9
To avoid such an exponentially accuracy loss, we abandon known results (such as Chebyshev
method) and design our own method. We prove the following theorem in this section:

N

Theorem D.1. Given a matrix , we reduce the problem of computing ξ ←
computations w0 ← w as follows.

N

N

N

N−1χ to multiple

If
satisfies
= B −1/2 N B 1/2 where N and B are both d × d positive definite matrices, for
every εe > 0 and χ ∈ Rd , in order to obtain ξ satisfying kξ − −1 χk ≤ εe,

e √κN times, and
• it suffices to compute w0 ← w only O

N

N
• each time of accuracy kw0 − Nwk ≤ O(1/poly(κB , εe, λmin (N ))).

Our reduction is based on an inexact variant of the accelerated gradient descent (AGD) method
originally put forward by Nesterov [23], which relies on some convex optimization techniques and
can be proved using the linear-coupling framework [5]. We prove this inexact AGD result in
Appendix D.1. Our final proof of Theorem D.1 is included in Appendix D.2.

D.1

Inexact Accelerated Gradient Descent

We study an inexact version of the classical accelerated gradient descent (AGD) method, and our
pseudocode is presented in Algorithm 4. The difference between our method and known AGD
e (xk+1 ) in each
methods is that we only require the algorithm to know an approximate gradient ∇f
e (xk+1 ) − ∇f (xk+1 )k2
iteration k, as opposed to the exact full gradient ∇f (xk+1 ). We require k∇f
to be upper bounded by some parameter εe in each iteration. Our next convergence theorem states
9
Indeed, for instance in the ALS algorithm of [29] for solving CCA, the authors obtained a running time proportional to 1/gap2 although there are only 1/gap iterations.

22

that this inexact AGD method only incurs an additive loss proportional to O(e
ε2 ).
Theorem D.2 (inexact AGD). If f (x) is L-smooth and σ-strongly convex, then AGDinexact (f, x0 , T )
produces an output yT satisfying
εe2 
f (yT ) − f (x∗ ) ≤ O(1) · (1 − τ )T (f (x0 ) − f (x∗ )) + O
,
σ
p
√
where τp
= Ω( σ/L). In other words, if the approximate gradient oracle satisfies εe ≤ O( εσ) and
T = O( L/σ · log(1/ε)), then we have f (yT ) − f (x∗ ) ≤ ε.

Theorem D.2 can be proved using the linear-coupling framework of [5]. In this framework,
accelerated methods are analyzed by a gradient descent lemma (Lemma D.3 below), a mirror
descent lemma (Lemma D.4 below), and a coupling step (Lemma D.5 and D.6 below).
Lemma D.3 (gradient descent). f (yk+1 ) ≤ f (xk+1 ) −

1
2
2L k∇f (xk+1 )k2

+

εe2
2L .

Proof. Abbreviating xk+1 by x and yk+1 by y, the smoothness property of function f (·) tells us
L
f (y) − f (x) ≤ h∇f (x), y − xi + ky − xk2 .
2
∇f (x)+χ
where kχk2 ≤ εe, we have
Now, since y − x = − L
L
−1
1
h∇f (x), y − xi + ky − xk22 =
h∇f (x), ∇f (x) + χi +
h∇f (x) + χ, ∇f (x) + χi
2
L
2L
1
εe2
≤ − k∇f (x)k22 +
.

2L
2L

Since our update on z can be written in the following minimization form, known as mirrordescent form in optimization literatures:
n1
o
(i)
e (xk+1 ), zi + ησ kz − xk+1 k22 .
zk+1 = min
kz − zk k22 + ηh∇f
(D.1)
z
2
2
It implies the following classical lemma (see for instance [6, Lemma 5.4]):
Lemma D.4 (mirror descent). For every u ∈ Rn ,
e (xk+1 ), zk+1 − ui − ησ kxk+1 − uk22 ≤ − 1 kzk − zk+1 k22 + 1 kzk − uk22 − 1 + ησ kzk+1 − uk22 .
ηh∇f
2
2
2
2
The following inequality is a nature linear combination of the two lemmas above:

Lemma D.5 (coupling 1). For every u ∈ Rn ,
ησ
ηh∇f (xk+1 ), zk − ui −
ku − xk+1 k22
2
 1
1 + ησ/2
η
η2
≤ η 2 L f (xk+1 ) − f (yk+1 ) + kzk − uk22 −
kzk+1 − uk22 + εe2 ( + ) .
2
2
σ
2

Proof. Combining Lemma D.3 and Lemma D.4 we deduce that for each i ∈ [n],
ησ
hη∇f (xk+1 ), zk − ui −
kxk+1 − uk22
2
e (xk+1 ), zk+1 − ui + εeηkzk+1 − uk2 − ησ kxk+1 − uk2
≤ hη∇f (xk+1 ), zk − zk+1 i + hη ∇f
2
2
¬
1
1
1 + ησ
≤ hη∇f (xk+1 ), zk − zk+1 i − kzk − zk+1 k22 + kzk − uk22 −
kzk+1 − uk22 + εeηkzk+1 − uk2
2
2
2
­ η2
1
1 + ησ/2
εe2 η
≤ k∇f (xk+1 )k22 + kzk − uk22 −
kzk+1 − uk22 +
2
2
2
σ
23

®
 1
1 + ησ/2
η
η2
≤ η 2 L f (xk+1 ) − f (yk+1 ) + kzk − uk22 −
kzk+1 − uk22 + εe2 ( + ) .
2
2
σ
2
Above, ¬ uses Lemma D.4, ­ uses the Young’s inequality which states 2ha, bi ≤ kak2 + kbk2 , ®
uses Lemma D.3.


Taking into account xk+1 = τ zk + (1 − τ )yk and the convexity of f (·), we can rewrite some
terms of Lemma D.5 and obtain
Lemma D.6 (coupling 2).
0≤

(1 − τ )η
η
1
1 + ησ/2
η η2
(f (yk )−f (x∗ ))− (f (yk+1 )−f (x∗ ))+ kzk −x∗ k22 −
kzk+1 −x∗ k22 + εe2 ( + )
τ
τ
2
2
σ 2

Proof.
¬

η(f (xk+1 ) − f (x∗ ))

ησ ∗
kx − xk+1 k22
2
ησ ∗
= ηh∇f (xk+1 ), xk+1 − zk i + ηh∇f (xk+1 ), zk − x∗ i −
kx − xk+1 k22
2
ησ ∗
­ (1 − τ )η
h∇f (xk+1 ), yk − xk+1 i + ηh∇f (xk+1 ), zk − x∗ i −
kx − xk+1 k22
=
τ
2
® (1 − τ )η
 1
1 + ησ/2
η
η2 
(f (yk ) − f (xk+1 )) + η 2 L f (xk+1 ) − f (yk+1 ) + kzk − uk22 −
kzk+1 − uk22 + εe2
+
.
≤
τ
2
2
σ
2
Above, ¬ is owing to the strong convexity of f (·), ­ uses the fact that xk+1 = τ zk + (1 − τ )yk , and
® uses the convexity of f (·) as well as Lemma D.5 with the choice of u = x∗ . Recall η = τ1L , we
arrive at the desired inequality.

≤ ηh∇f (xk+1 ), xk+1 − x∗ i −

We are now ready to prove Theorem D.2.
Proof of Theorem D.2. We choose τ =

√2

1+

8L/σ+1

1
1−τ .

∈ [0, 1), and this choice ensures that 1 + ησ/2 =

Under these parameter choices, Lemma D.6 becomes



τ
1 η
τ
kzk+1 −x∗ k22 ≤ (1−τ ) (f (yk )−f (x∗ ))+
kzk −x∗ k22 +e
ε2 τ +
f (yk+1 )−f (x∗ ) +
2η(1 − τ )
2η(1 − τ )
σ 2
Telescoping it for all iterations k = 0, 1, . . . , T − 1, we conclude that


τ
εe2 
1 η
f (yT )−f (x∗ ) ≤ (1−τ )T f (y0 )−f (x∗ )+ kz0 −x∗ k22 +e
.
ε2 ( + ) ≤ O(1)·(1−τ )T (f (x0 )−f (x∗ ))+O
2η
σ 2
σ
where the last inequality is because (i) x0 = y0 = z0 , (ii) τ /η = O(σ) and (iii) the strong convexity
of f (·) which implies f (x0 ) − f (x∗ ) ≥ σ2 kx0 − x∗ k22 .


D.2

Proof of Theorem D.1

Proof of Theorem D.1. We first verify accuracy. Since
kξ −

p

N−1χk ≤ εe ⇐= kB1/2ξ − B1/2N−1χk ≤ εe · pλmin(B)
⇐= kB 1/2 ξ − N −1 B 1/2 χk ≤ εe ·

λmin (B) ,

N

(D.2)

−1 χk

it suffices to find ξ to satisfy (D.2) in order to satisfy the accuracy requirement kξ −
≤ εe.
>
def 1 >
1/2
∗
∗
−1
Define f (x) = 2 x N x − B χ x and let x be its minimizer. Then it satisfies x = N B 1/2 χ
24

and f (x) − f (x∗ ) = 12 (x − x∗ )> N (x − x∗ ). For this reason, it suffices to find an approximate
minimizer of f (x) satisfying
1
εe2
(x − x∗ )> N (x − x∗ ) = f (x) − f (x∗ ) ≤ λmin (B)λmin (N ) =: εe0
(D.3)
2
2
2
because if we let ξ = B −1/2 x then the above inequality implies 21 kx − x∗ k2 ≤ εe2 · λmin (B) which is
the same as (D.2). In sum, we can call AGDinexact to find an approximate minimizer x with additive
error no more than εe0 , and then defining ξ = B −1/2 x gives a solution of ξ satisfying kξ − −1 χk ≤ εe.

N

AGDinexact .

We now focus on the actual implementation of
If we choose x0 = 0 as the initial
1/2
1/2
1/2
vector, we can write xk , yk , zk implicitly as xk = B
k , yk = B
k , zk = B
k (thus only keep
track of k , k , k ) throughout the algorithm. Under these notations, we claim that it suffices to
perform matrix vector multiplication p
on
(i.e., of the form w0 ←  w) for at most O(T ) times on
those implicit vectors where T = O λmax (N )/λmin (N ) log(1/e
ε0 ) is the number of iterations of
inexact
AGD
according to Theorem D.2.

This is so because ∇f (xk ) = N xk − B 1/2 χ = B 1/2
k − χ and therefore for instance yk+1 ←
e (xk+1 ) can be implemented as k+1 ← k+1 − 1 ( k+1 − χ) so only matrix-vector
xk+1 − L1 ∇f
L
0 ←
multiplication on
is needed. In
addition,
as
long
as
each
w
w is computed to an additive
p

0
error kw − wk ≤ O εe· λmin (N ) λmin (B)/λmax (B) , we can use B 1/2 (w0 − χ) as the approximate
p
gradient which is different from the true gradient ∇f (xk ) by an additive amount O( λmin (N )e
ε0 ).
This satisfies the approximation require of Theorem D.2, and thus the accuracy guarantee provided
by Theorem D.2 is satisfied.


x

x y z

N

E
E.1

y
N

N

Nx
x

y

N

y

Nx
N

Proof for Section 4: GenEV Theorems
Proof of Theorem 4.1

In this section we prove Theorem 4.1 formally.
Theorem 4.1 (restated). Let M ∈ Rd×d be a symmetric matrix with eigenvalues λ1 , . . . , λd ∈
[−1, 1] and corresponding eigenvectors u1 , . . . , ud . Suppose without loss of generality that |λ1 | ≥
· · · ≥ |λd |.

Suppose k ∈ [d], δ× , p ∈ (0, 1). Then, LazyEV outputs a (column) orthonormal matrix Vk =
(v1 , . . . , vk ) ∈ Rd×k which, with probability at least 1 − p, satisfies all of the following properties.
(Denote by Mk = (I − Vk Vk> )M (I − Vk Vk> ).)
4

ε δ×
>
(a) Core lemma: if εpca ≤ 212 k3 (|λ
2 , then kVk U k2 ≤ ε, where U = (uj , . . . , ud ) is the
1 |/|λk |)
(column) orthonormal matrix and j is the smallest index satisfying |λj | ≤ (1 − δ× )kMk−1 k2 .
5
δ×
, then |λk+1 |
228 k3 (|λ1 |/|λk+1 |)6
5
δ×
εpca ≤ 228 k3 (|λ1 |/|λk+1 |)6 , then

(b) Spectral norm guarantee: if εpca ≤

(c) Rayleigh quotient guarantee: if
1
1−δ× |λk |.

(d) Schatten-q norm guarantee: for every q ≥ 1, if εpca ≤
kMk kSq ≤

d
(1 + δ× )2  X q 1/q
(1 + δ× )2
λ
=
i
(1 − δ× )2
(1 − δ× )2
i=k+1

≤ kMk k2 ≤

(1 − δ× )|λk | ≤ |vk> M vk | ≤

5
δ×
,
228 k3 d4/q (|λ1 |/|λk+1 |)6

min

V ∈Rd×k ,V > V =I



then

k(I−V V > )M (I−V V > )kSq

Proof of Theorem 4.1. Let Vs = (v1 , . . . , vs ), so we can write
Ms = (I − Vs Vs> )M (I − Vs Vs> ) = (I − vs vs> )Ms−1 (I − vs vs> ) .
25

|λk+1 |
1−δ× .

.

We first claim that kMs−1 k2 ≥ |λs | for every s = 1, . . . , k. This can be proved by the Cauchy
2
> )M 2 (I − V
>
2
interlacing theorem. Indeed, Ms−1
= (I − Vs−1 Vs−1
s−1 Vs−1 ) is a projection of M into
2
a d − s + 1 dimensional space, and therefore its largest eigenvalue kMs−1 k2 should be at least as
large as |λs |2 , the s-th largest eigenvalue of M 2 . In other words, we have shown kMs−1 k2 ≥ |λs |.
b = kMk−1 k2 ≥ |λk |.
(a) Define λ

Note that all column vectors in Vs are automatically eigenvectors of Ms with eigenvalues zero.
For analysis purpose only, let Ws be the column matrix of eigenvectors in Vs⊥ of Ms that have
s
b where τs def
= 2k
absolute eigenvalues in the range [0, (1 − δ× + τs )λ],
δ× . We now show that for
every s = 0, . . . , k, there exists a matrix Qs such that kU − Ws Qs k2 is small and kQs k2 ≤ 1.
We will do this by induction.
In the base case: since τ0 = 0, we have W0 = U by the definition of U . We can therefore define
Q0 to be the identity matrix.
For every s = 0, 1, . . . , k − 1, suppose there exists a matrix Qs with kQs k2 ≤ 1 that satisfies
kU − Ws Qs k2 ≤ ηs for some ηs > 0, we construct Qs+1 as follows.

0
0> W k2 ≤ ε
First we observe that AppxPCA± outputs a vector vs+1
satisfying kvs+1
s 2
pca and
0>
2
kvs+1 Vs k2 ≤ εpca with probability at least 1 − p/k. This follows from Theorem 3.1 (usb ⊆ [0, (1 − δ× /2)λ],
b together with the fact that
ing M = Ms ) because [0, (1 − δ× + τs )λ]
0
b
into Vs⊥ , we have
kMs k2 ≥ kMk−1 k2 ≥ λ. Now, since vs+1 is the projection of vs+1
>
kvs+1
Ws k22 ≤

0> W k2
0> W k2
kvs+1
kvs+1
εpca
s 2
s 2
=
≤
< 1.5εpca .
0
0
1 − εpca
k(I − Vs Vs> )vs+1
k22
1 − kVs> vs+1
k22

(E.1)

Next we apply Lemma C.4 with M = Ms , M 0 = Ms+1 , U = Ws , V = Ws+1 , v = vs+1 ,
b and τ = (τs+1 − τs )λ.
b We obtain a matrix Q
e s , kQ
e s k2 ≤ 1 such that10
µ = (1 − δ× + τs )λ,
s
√
b 2 · 1.5εpca
32λ1 k εpca
169(λ1 /λ)
f
+ εpca <
,
kWs − Ws+1 Qs k2 ≤
λk δ×
(τs+1 − τs )2
and this implies that
√
32λ1 k εpca
f
f
kWs+1 Qs Qs − U k2 ≤ kWs+1 Qs Qs − Ws Qs k2 + kWs Qs − U k2 ≤ ηs +
.
λk δ×
fs Qs we know that kQs+1 k2 ≤ 1 and
Let Qs+1 = Q
kWs+1 Qs+1 − U k2 ≤ ηs+1

√
32λ1 k εpca
.
= ηs +
λk δ×

def

Therefore, after k-iterations of LazyEV, we obtain:
kWk Qk − U k2 ≤ ηk =

√
32λ1 k 2 εpca
λk δ×

Multiply U > from the left, we obtain kU > Wk Qk − Ik2 ≤ ηk . Since kQk k2 ≤ 1, we must have
σmin (U > Wk ) ≥ 1 − ηk (here σmin denotes the smallest singular value). Therefore,
U > Wk Wk> U  (1 − ηk )2 I .

Since Vk and Wk are orthogonal of each other, we have
U > Vk Vk> U  U > (I − Wk Wk> )U  I − (1 − ηk )2 I  2ηk I
10

Technically speaking, to apply Lemma C.4 we need U = Ws to consist of all eigenvectors of Ms with absolute eigenvalues ≤ µ. However, we only defined Ws to be such eigenvectors that are also orthogonal to Vs . It is
straightforward to verify that the same result of Lemma C.4 remains true because vs+1 is orthogonal to Vs .

26

Therefore,
1/4

kVk> U k2 ≤ 8

(|λ1 |/|λk |)1/2 kεpca
1/2

δ×

≤ε .

(b) The statement is obvious when k = 0. For every k ≥ 1, the lower bound is obvious. We prove
|λk+1 |
the upper bound by contradiction. Suppose that kMk k2 > 1−δ
. Then, since kMk−1 k2 ≥
×
kMk k2 and therefore |λk+1 |, . . . , |λd | < (1 − δ× )kMk−1 k2 , we can apply Theorem 4.1.a of
def
the current k to deduce that kVk> U>k k2 ≤ ε where U>k = (uk+1 , . . . , ud ). We now apply
Lemma C.2 with Vs = Vk and U = U>k , we obtain a matrix Qk ∈ Rd×k whose columns are
spanned by u1 , . . . , uk and satisfy

 
 
 

>
>
< 16|λ1 |ε .
I − Qk Q>
M I − Vk Vk>
k M I − Qk Qk − I − Vk Vk
2


> has spectral
However, our assumption says that the second matrix I − Vk Vk> M I − V
k Vk

>
norm at least |λk+1 |/(1−δ× ), but we know that I − Qk Q>
k M I − Qk Qk has spectral norm
|λk+1 |
1−δ×

− |λk+1 | ≤ 16|λ1 |ε

δ×
16|λ1 |/|λk+1 |

(which is satisfied

exactly |λk+1 | due to the definition of Qk . Therefore, we must have
due to triangle inequality.
In other words, by selecting ε in Theorem 4.1.a to satisfy ε ≤

by our assumption on εpca ), we get a contradiction so can conclude that kMk k2 ≤

|λk+1 |
1−δ× .

(c) We compute that
¬

|vk> M vk | = |vk> Mk−1 vk | ≥

­ |v 0> Mk−1 v 0 |
|vk0> Mk−1 vk0 |
k
k
≥
√
2
> )v 0 k2
(1
−
ε
)
k(I − Vk−1 Vk−1
pca
k 2
®

≥

1 − εpca
(1 − δ× /2)kMk−1 k2 ≥ (1 − δ× )kMk−1 k2 . (E.2)
√
(1 − εpca )2

⊥ , ­ is because kV > v 0 k2 ≤ ε
Above, ¬ is because vk is the projection of vk0 into Vk−1
pca following
k−1 k 2
the same reason as (E.1), and ® is owing to Theorem 3.1. Next, since kMk−1 k2 ≥ |λk |, we
automatically have |vk> M vk | ≥ (1 − δ× )|λk |. On the other hand, |vk> M vk | = |vk> Mk−1 vk | ≤
|λk |
kMk−1 k2 ≤ 1−δ
where the last inequality is owing to Theorem 4.1.b.
×
def

(d) Since kVk> U k2 ≤ εc = 8

1/4

(|λ1 |/|λk |)1/2 kεpca
1/2
δ×

from the analysis of Theorem 4.1.a, we can apply

Lemma C.2 to obtain a (column) orthogonal matrix Qk ∈ Rd×k such that
kMk0 − Mk k2 ≤ 16|λ1 |εc ,

>
where Mk0 = (I − Qk Q>
k )M (I − Qk Qk )
def

(E.3)

Suppose U = (ud−p+1 , . . . , ud ) is of dimension d × p, that is, there are exactly p eigenvalues of
M whose absolute value is ≤ (1 − δ× )kMk−1 k2 . Then, the definition of Qk in Lemma C.2 tells
us U > Qk = 0 so Mk0 agrees with M on all the eigenvalues and eigenvectors {(λj , uj )}dj=d−p+1
because an index j satisfies |λj | ≤ (1−δ× )kMk−1 k2 if and only if j ∈ {d−p+1, d−p+2, . . . , d}.
Denote by µ1 , . . . , µd−k the eigenvalues of Mk0 excluding the k zero eigenvalues in subspace Qk ,
and assume without loss of generality that {µ1 , . . . , µp } = {λd−p+1 , . . . , λd }. Then,
kMk0 kqSq =
¬

≤

d−k
X
i=1

|µi |q =

d
X

i=d−p+1

p
X
i=1

|µi |q +

d−k
X

i=p+1

|µi |q =
­

i=d−p+1

|λi |q + (d − k − p)kMk0 kq2 ≤

27

d
X

d
X

|λi |q +

i=d−p+1

d−k
X

i=p+1

|µi |q

|λi |q + (d − k − p)(kMk k2 + 16|λ1 |εc )q

®

≤

d
X

i=d−p+1

q
 |λ |
k+1
+ 16|λ1 |εc
|λi |q + (d − k − p)
(1 − δ× )

Above, ¬ is because each |µi | is no greater than kMk0 k2 , and ­ is owing to (E.3), and ® is
|λ
|δ×
(and this is indeed
because of Theorem 4.1.b. Suppose we choose εc so that εc ≤ k+1
16λ1
satisfied by our assumption on εpca ), then we can continue and write
kMk0 kqSq ≤

d
X

i=d−p+1

|λi |q + (d − k − p)
d
X

¯

≤

i=d−p+1

(1 + δ× )q
|λk+1 |q
(1 − δ× )q

|λi |q +

d−p
d
(1 + δ× )q X
(1 + δ× )q X
q
|λ
|
≤
|λi |q .
i
(1 − δ× )2q
(1 − δ× )2q
i=k+1

i=k+1

Above, ¯ is because for each eigenvalue λi where i ∈ {k + 1, k + 2, . . . , d − p}, we have
|λi | > (1 − δ× )kMk−1 k2 ≥ (1 − δ× )|λk | ≥ (1 − δ× )|λk+1 |. Finally, using (E.3) again we have
kMk kSq ≤ kMk0 kSq + kMk − Mk0 kSq ≤ kMk0 kSq + d1/p kMk − Mk0 k2
d
1/q
1 + δ×  X
q
≤
+ 16d1/p |λ1 |εc
|λ
|
i
(1 − δ× )2
i=k+1

As long as εc ≤

δ× |λk+1 |
,
16d1/p λ1

we have
kMk kSq ≤

d
1/q
(1 + δ× )2  X
q
|λ
|
i
(1 − δ× )2
i=k+1

δ× λk+1
as desired. Finally, we note that εc ≤ 16d
is
1/p λ
1

note that minV ∈Rd×k ,V > V =I k(I − V V > )M (I − V

easily from Cauchy interlacing theorem.

E.2

satisfied with our assumption on εpca , and

Pd
q 1/q which follows
V > )kSq =
i=k+1 |λi |



Proofs of Theorems 4.3 and 4.4

Vk = LazyEV(· · · ) to be the direct output of LazyEV. The
V Vk = I.

Proof of Theorem 4.3. Define Vk = B 1/2
column orthogonality of Vk implies k> B

It is clear from the definition of generalized eigenvectors that B 1/2 u1 , . . . , B 1/2 ud are eigenvectors
def
of M = B −1/2 AB −1/2 with eigenvalues λ1 , . . . , λd . Applying Theorem 4.1.a, we have: kVk> U k2 ≤ ε
where U = (B 1/2 uj , . . . , B 1/2 ud ) is a (column) orthonormal matrix and j is the smallest index
satisfying |λj | ≤ (1 − δ× )kMk−1 k2 . Since it satisfies kMk−1 k2 ≥ |λk |, we have
|λk+1 | = |λk |(1 − gap) = |λk |(1 − δ× ) ≤ (1 − δ× )kMk−1 k2 .

Therefore, j must be equal to k +1 according to its definition, so we have U = B 1/2
k k> B k2 = kVk> U k2 ≤ ε.

V W

W. This implies

The running time statement comes directly from Theorem 4.2 by putting in the parameters.



Vk = LazyEV(· · · ) to be the direct output of LazyEV. The
V Vk = I.

Proof of Theorem 4.4. Define Vk = B 1/2
column orthogonality of Vk implies k> B

It is clear from the definition of generalized eigenvectors that B 1/2 u1 , . . . , B 1/2 ud are eigenvectors
def
of M = B −1/2 AB −1/2 with eigenvalues λ1 , . . . , λd . Applying Theorem 4.1.b, we have: (I −
Vk Vk> )B −1/2 AB −1/2 (I − Vk Vk> )

2

≤

|λk+1 |
1−ε .

Next, for every vector
28

w ∈ Rd that is B-orthogonal to

Vk , that is, w>BVk = 0, we can define w = B1/2w and we know w is orthogonal to Vk . We can
def

apply the above spectral upper bound and get

w>Aw = w>B−1/2AB−1/2w = w>(I − Vk Vk>)B−1/2AB−1/2(I − Vk Vk>)w
as desired. At the same time, denoting by vs = B 1/2
∀s ∈ [k] :

vs>Avs

= vs> M vs

|
w w · |λ1 k+1
−ε

|λk+1 |
= >B
1−ε
s , Theorem 4.1.c implies that
h
|λs | i
∈ (1 − ε)|λs |,
.
1−ε
≤ kwk22 ·

v

The running time statement comes directly from Theorem 4.2 by putting in the parameters.

F
F.1



Proof for Section 6: CCA Theorems
The Main Convergence Theorem

Since LazyCCA only admits minor changes on top of LazyEV, the next theorem is an almost identical
copy of Theorem 4.1. To make this paper concise, instead of reproving Theorem F.1 line by line,
we here only sketch the main changes needed in the new proof.
Theorem F.1. Let M = B −1/2 AB −1/2 ∈ Rd×d be a symmetric matrix where A and B are
matrices coming from a CCA instance using Lemma 2.3. Suppose M has eigenvalues λ1 , . . . , λd ∈
[−1, 1] and corresponding eigenvectors u1 , . . . , ud . Suppose without loss of generality that |λ1 | ≥
· · · ≥ |λd |.

1|
For every k ∈ [d], δ× , p ∈ (0, 1), there exists some εpca ≤ O poly(δ× , |λ|λk+1
, d1 ) such that
LazyCCA outputs a (column) orthonormal matrix Vk = (v1 , . . . , v2k ) ∈ Rd×2k which, with probability at least 1−p, satisfies all of the following properties. (Denote by Ms = (I−Vs Vs> )M (I−Vs Vs> ).)
(a) Correlation guarantee: kVk> U k2 ≤ ε,
where U = (uj , . . . , ud ) and j is the smallest index satisfying |λj | ≤ (1 − δ× )kMk−1 k2 .
|λ2k+1 |
(b) Spectral norm guarantee: |λ2k+1 | ≤ kMk k2 ≤ 1−δ
.
×
1
> Mv | ≤
>
(c) Rayleigh quotient guarantee: (1 − δ× )|λ2k | ≤ |v2k−1 M v2k−1 | = |v2k
2k
1−δ× |λ2k |.
(d) Schatten-q norm guarantee: for every q ≥ 1,
d

(1 + δ× )2  X q 1/q
(1 + δ× )2
kMk kSq ≤
min
k(I−V V > )M (I−V V > )kSq
λ
=
i
(1 − δ× )2
(1 − δ× )2 V ∈Rd×2k ,V > V =I
i=2k+1

Proof sketch of Theorem F.1. Recall that when a vector vs ∈ Rd is obtained in iteration s of LazyEV,
the proof of Theorem 4.1 suggest that the following two properties hold
vs> Ms−1 vs ≥ (1 − δ× )kMs−1 k2

and

vs> Ws−1

2
2

≤ 1.5εpca .

(F.1)

(The first property is shown in (E.2), and the second property is shown in (E.1). Recall that Ws−1
is the column orthonormal matrix containing all eigenvectors of Ms−1 whose absolute eigenvalues
are below some threshold.) Then, the proof of Theorem 4.1 proceeds by heavily relying on (F.1).
In our LazyCCA, after
this same√vector vs , we write it as vs = (ξs0 , ζv0 ) and perform
√ obtaining
0
0
block-scaling ξs = ξs /( 2kξs k) and ζs = ζs0 /( 2kζs0 k), see Line 7 of LazyCCA. Therefore, in order
for the same proof of Theorem 4.1 to hold, we need to show that this new vector (ξs , ζs ) satisfies

29

.

the same properties up to constants:
 ξ 
 ξ >
 ξ >
2
s
s
s
(F.2)
Ws−1 ≤ 12.5εpca .
Ms−1
≥ (1 − δ× )kMs−1 k2 and
ζs
ζs
ζs
2

 ±ξ · · · ±ξ
1
s−1
. Since vs is orthogonal to all vectors in the column span of
Suppose Vs−1 =
ζ1 · · · ζs−1
Vs−1 according to Line 5 of LazyCCA, we automatically have ξs> ξi = 0 and ζs> ζi for all i ∈ [s − 1].
We also have kξs k2 + kζs k2 = 1/2 + 1/2 = 1 so the new vector (ξs , ζs ) has Euclidean norm 1.
As for the first property in (F.2), we observe that the new vector (ξs , ζs ) enjoys an (absolute)
Rayleigh quotient value that is no worse than the original vs = (ξs0 , ζs0 ). This is so because (without
loss of generality we consider vs> Ms−1 vs > 0):
 ξ   ξ >  ξ 
 ξ >

­
¬
s
s
s
s
−1/2
−1/2
= 2ξs> Sxx
Sxy Syy
ζs
M
=
Ms−1
ζs
ζs
ζs
ζs
±

1
1
1
®
°
−1/2
−1/2 0 ¯
= 0
· ξs0> Sxx
Sxy Syy
ζs =
· vs> M vs =
· vs> Ms−1 vs ≥ vs> Ms−1 vs .
0
0
0
0
0
kξs kkζs k
2kξs kkζs k
2kξs kkζs k
(F.3)
Above, ¬ is because (ξs , ζs ) is orthogonal to Vs−1 ; ­ is by the definition of M = B −1/2 AB −1/2 as
well as the definition of A and B; ® is by the definitions of ξs and ζs ; ¯ is by vs = (ξs0 , ζs0 ) and
again by the definition of M ; ° follows from the fact that vs is orthogonal to Vs−1 ; and ± follows
from AM-GM together with the fact that kξs0 k2 + kζs0 k2 = kvs k2 = 1. This finishes proving the first
property in (F.2) because the original vector vs satisfies |vs> Ms−1 vs | ≥ (1 − δ× )kMs−1 k2 according
to (F.1).
We make an additional observation here: kξs0 k2 and kζs0 k2 must be in the range [0.06, 0.94] before
scaling. Indeed, suppose for instance kξs0 k2 = c for some c ∈ [0, 1]. Then, it satisfies 2kξs0 kkζs0 k =
 ξ >
 ξ 
p
s
s
2 c(1 − c) and therefore (F.3) becomes
Ms−1
≥ √1−δ× kMs−1 k2 , meaning that
2 c(1−c)
ζs
ζs
√
√
√1−δ× ≤ 1. If δ× ≤ 1/2, this implies c − 1/2 ∈ [− 3/4, 3, 4] and thus c ∈ [0.06, 0.94].
2

c(1−c)

As for the second property in (F.2), for every matrix Ws−1 that is in the proof of Theorem 4.1.a,
its columns are all eigenvectors of Ms−1 whose absolute eigenvalues are below some threshold, so
 ±a . . . ±a 
1
t
must consist of only symmetric vectors in this CCA setting: that is, Ws−1 =
.11
b1 . . . bt
According to (F.1) we already know kvs> Ws−1 k2 ≤ 1.5εpca , which implies
kvs> Ws−1 k2 =

t
X
i=1

(ξs0> ai +ζs0> bi )2 +(ξs0> ai −ζs0> bi )2 = 2kξs0> (a1 , . . . , at )k2 +2kζs0> (b1 , . . . , bt )k2 ≤ 1.5εpca .

Now we can compute
 ξ >
2
s
Ws−1 = 2kξs> (a1 , . . . , at )k2 + 2kζs> (b1 , . . . , bt )k2
ζs
2

0.5
≤
· 2kξs0> (a1 , . . . , at )k2 + 2kζs0> (b1 , . . . , bt )k2 ≤ 12.5εpca ,
0.06

This is because matrix Ms is always of the form D−1/2 CD1/2 where D = diag{D1 , D2 } is block diagonal and
positive definite, while C = [[0, C1 ]; [C1> , 0]] has only zero on its two block diagonal locations. The same proof of
Lemma 2.3 shows that the eigenvectors of D−1/2 CD1/2 must be symmetric. In fact, to be precise, Ws may also
contain some eigenvectors corresponding to zero eigenvalues. However, adding them will make our notations heavier,
so we refrain from doing that in this sketched proof.
11

30

where the first inequality is because kξs0 k2 , kζs0 k2 ≥ 0.06.

This finishes proving the two properties in (F.2), so Theorem F.1 holds after plugging the rest
of the proof of Theorem 4.1 in but changing constants slightly.


F.2

Fast Implementation of LazyCCA: Stochastic

Theorem F.2 (stochastic running time of LazyCCA). Let X ∈ Rn×dx , Y ∈ Rn×dy be two matrices,
and define A and B according to Lemma 2.3. Suppose M = B −1/2 AB −1/2 , and RanInit(d) is
the random vector generator defined in Proposition 3.3, and we want to compute matrix
←
−1/2
B
LazyCCA(A, M, k, δ× , εpca , p). Then, this procedure can be implemented to run in time


2
e knnz(B)+k
√ d+kΥ where Υ is the time needed to compute B −1 Aw for a vector w to an
• O

V

δ×

e
accuracy ε where log(1/ε) = O(1),
or

 knnz(X,Y ) 1+√κ0 /n +k2 d 
e
√
• O
if we simply use Katyusha to compute B −1 Aw.
δ×

Above, κ =

λmax (B)
λmin (B)

=

max{λmax (Sxx ),λmax (Syy )}
min{λmin (Sxx ),λmin (Syy )} ,

and κ0 =

2 maxi∈[n] {kXi k2 ,kYi k2 }
λmin (B)

∈ [κ, 2nκ].

Proof. The proof of the first item is almost identical to the proof of Theorem 4.2 so ignored here.
As for the second item, it follows from the first item together with the running time of Katyusha
in Lemma 2.6.


F.3

Fast Implementation of LazyCCA: Doubly Stochastic

Theorem F.3 (doubly-stochastic running time of LazyCCA). Let X ∈ Rn×dx , Y ∈ Rn×dy be
two matrices, and define A and B according to Lemma 2.3. Suppose M = B −1/2 AB −1/2 , and
RanInit(d) is the random vector generator defined in Proposition 3.3, and we want to compute
← B −1/2 LazyCCA(A, M, k, δ× , εpca , p). Then, this procedure can be implemented to run in time
√


κ0 /n1/4 
e nnz(X, Y ) · 1 + √
if k = 1 and we use accelerated SVRG as the method A;
• O
δ× σ1

√

κ0
e
√
• O knnz(X, Y ) · 1 +
if nnz(X, Y ) ≥ kd and we use accelerated SVRG
1/4

V

as the method A.

Above, κ =

λmax (B)
λmin (B)

=

δ× σk (nnz(X,Y )/kd)

max{λmax (Sxx ),λmax (Syy )}
min{λmin (Sxx ),λmin (Syy )} ,

and κ0 =

2 maxi∈[n] {kXi k2 ,kYi k2 }
λmin (B)

∈ [κ, 2nκ].

To prove Theorem F.3, let us recall from the proof of Theorem 4.2 (see Section 5.2) that it
suffices to bound the time needed to compute −1 w where

def
= B −1/2 λI − Ms B 1/2 = λI − (I − s s> B)B −1 A(I − s s> B) ,

N

N

VV

VV

and we only need to compute it poly-logarithmic number of times for each s = 1, . . . , k.

N

N

Observe now that −1 w = (B )−1 Bw, so it suffices to bound the time needed to compute
(B )−1 applied to a vector, and we have (dropping the subscript of
for cleanness)

N

B
F.3.1

N = λB − (I − BVV

>

)A(I −

VV

>

V

B) .

Special Case of 1-CCA

V

When k = 1, the matrix
is empty and we want to compute (λB − A)−1 w for some vector w.
def
This is equivalent to minimizing f (z) = 21 z > (λB − A)z − w> z. We analyze its running time as
follows, and a similar analysis has also appeared in [29, Section 3.2.3].
31

P
Using the definition of A and B, we can write f (z) = n1 ni=1 fi (z) where
def 1
fi (z) = λ(hXi , z1 i2 + hYi , z2 i2 ) − 2hXi , z1 ihYi , z2 i + w> z ,
2
where we have denoted by z = (z1 , z2 ) for z1 ∈ Rdx and z2 ∈ Rdy , and by Xi , Yi the i-th row vector
of matrix X and Y respectively. The smoothness of fi (z) is given by (using the fact that λ ≤ 2
which comes from the definition of AppxPCA± )
k∇2 fi (z)k2 ≤ O(1) · max{kXi k2 , kYi k2 } .

The strong convexity of f (z) is given by

δ× σ1
σ1 δ× λmin (B)
I)B 1/2 
·I ,
48
48
∗
where the first inequality is because λmin (λI − B −1/2 AB −1/2 ) ≥ δ×48λ from Theorem 3.1, as well as
λ∗ ≥ σ1 = λmax (B −1/2 AB −1/2 ) from the definition of the algorithm.
∇2 f (z) = B 1/2 (λI − B −1/2 AB −1/2 )B 1/2  B 1/2

Therefore, by applying the accelerated SVRG method on minimizing the sum-of-non-convex
objective f (z) [7, 24], we can find an εe minimizer of f (z) in time
√



 pmax {kX k2 , kY k2 } 
0

kwk
kwk 
i
i
i
2
e nnz(X, Y )· 1+ p κ
e nnz(X, Y )· 1+ p
·log
=O
.
·log2
O
εe
εe
σ1 δ× n1/4
σ1 δ× λmin (B)n1/4

Finally, similar to the proof of Lemma 2.6, an εe minimizer of f (z) implies (λB − A)−1
w − z ≤ ε
e nnz(X, Y ) ·
when ε2 = 96e
ε/(σ1 δ× λmin (B)). Putting this back we obtain a final running time of O

√

0
1 + √ κ 1/4 , ignoring log factors.
σ 1 δ× n

F.3.2

Genearl Case of k-CCA

V is not empty, one can carefully check that to compute
(λB − (I − B VV> )A(I − VV> B))−1 w = (λB − A + (B VV> A + AVV> B − B VV> AVV> B)−1 w
For the general case when

P
it suffices to minimize f 0 (z) = n1 ni=1 fi0 (z) for
1
def
def
>
>
>
A+A
B−B
A
fi0 (z) = fi (z) + z > Qz where Q = B
2
and fi (z) is the same as defined in the k = 1 case of the previous subsection.

VV

Claim F.4. k∇2 fi0 (z)k2 ≤ O(1) · max{kXi k2 , kYi k2 }
Proof. Since

VV

and

∇2 f 0 (z) 

VV VV>B

σk δ× λmin (B)
48

,

·I .

1/2
V
1/2
kAVk2 ≤ kB 1/2 k2 kB −1/2 AB −1/2 k2 kB 1/2 Vk2 ≤ kBk2 σ1
kV> AV> k2 ≤ kV> B 1/2 k2 kB −1/2 AB −1/2 k2 kB 1/2 Vk2 ≤ σ1

kB k2 ≤ kB 1/2 k2 kB −1/2 V k2 ≤ kBk2

we can easily compute that kQk2 ≤ O(1) · kBk2 ≤ O(1) · Tr(B) ≤ O(1) · maxi {kXi k2 , kYi k2 }.
Therefore, we have the same smoothness property as in the previous subsection:
k∇2 fi0 (z)k2 ≤ O(1) · max{kXi k2 , kYi k2 } .

The strong convexity of f 0 (z) is given by

∇2 f 0 (z) = B 1/2 (λI − Ms )B 1/2  B 1/2

32

δ× σk
σk δ× λmin (B)
I)B 1/2 
·I ,
48
48

∗

where the first inequality is because (1) λmin (λI −Ms ) ≥ δ×48λ from Theorem 3.1, (2) λ∗ ≥ λmax (Ms )
from the definition of the algorithm, as well as (3) λmax (Ms ) ≥ σk from the fact that Ms is exactly
B −1/2 AB −1/2 but projecting out at most k − 1 dimensions.

Claim F.5. With a preprocessing time O(nnz(X, Y )), we can compute (∇fi0 (z))w for an arbitrary
vector w and an arbitrary index i ∈ [n] in time O(kd), and (∇f 0 (z))w in time O(nnz(X, Y ) + kd).
Proof. Note that (∇fi (z))w only costs running time O(d) and thus it suffices to focus on the
computation of Qw. We first note that one can pre-compute A , B in time O(nnz(X, Y )) with
the help from the previous iteration. This is so because, suppose
= [ 0 , v] where v is the new
column vector added in the current iteration. Then, A = [A 0 , Av] and this additional column
Av can be computed in time O(nnz(X, Y )). Next, we write
 > 
 > 
 > 

>
Q= B
A + A
B − B
A
B

V V
V V
V

V

V V

V V

V V

VV

and every matrix between parenthesis of the above formulation has either only O(k) rows or O(k)
columns. Therefore, computing Qw costs a total running time at most O(kd) and so is (∇fi0 (z))w.

Finally, we have (∇f 0 (z))w = (λB − A)w + Qw and therefore computing (∇f 0 (z))w costs time
O(nnz(X, Y ) + kd).

Since computing (∇fi0 (z))w requires time O(kd) and computing (∇f 0 (z))w requires time O(nnz(X, Y )+
kd), we can apply SVRG to minimize this sum-of-non-convex function f 0 (z) [7, 24], with a running
time




2
2 2
0 2

e nnz(X, Y ) + kd + (max{kXi k , kYi k }) kd = O
e nnz(X, Y ) + (κ ) kd .
O
2 σ2λ
2 σ2
2
δ×
δ×
k min (B)
k

Now, as long as nnz(X, Y ) ≥ kd, we can apply acceleration on top of SVRG [7, 24] to given an
accelerated running time
√



κ0
e nnz(X, Y ) · 1 +
O
.
1/2 1/2
δ× σk (nnz(X, Y )/(kd))1/4
Taking into account all the iterations s = 1, . . . , k concludes the proof of Theorem F.3.

F.4

Proofs of Theorems 6.2 and 6.3

The two corollaries follow from Theorem F.1 for the similar reason as Theorem 4.3 and Theorem 4.4
following from Theorem 4.1.

V

W

V

W

Proof Theorem 6.2. The approximation guarantees k φ> Sxx φ k2 ≤ ε and k ψ> Syy ψ k2 ≤ ε follow from Theorem F.1.a, and the running time follows from Theorem F.2 and Theorem F.3. 
n
Proof of Theorem 6.3. The approximation guarantees maxφ∈Rdx ,ψ∈Rdy φ> Sxy ψ φ> Sxx φ = 0 ∧
o
ψ > Syy ψ = 0 ≤ (1 + ε)σk+1 follow from Theorem F.1.b and the definition of M , and the approx-

V

V

imation guarantee (1 − ε)σi ≤ |φ0i Sxy ψi | ≤ (1 + ε)σi follows from Theorem F.1.c and the fact that
|λ2i−1 | = |λ2i | = σi . The running time follows from Theorem F.2 and Theorem F.3.


References
[1] Zeyuan Allen-Zhu. Katyusha: The First Direct Acceleration of Stochastic Gradient Methods.
ArXiv e-prints, abs/1603.05953, March 2016.
33

[2] Zeyuan Allen-Zhu and Yuanzhi Li. Even Faster SVD Decomposition Yet Without Agonizing
Pain. In NIPS, 2016.
[3] Zeyuan Allen-Zhu and Yuanzhi Li. Fast Global Convergence of Online PCA. ArXiv e-prints,
abs/1607.07837, July 2016.
[4] Zeyuan Allen-Zhu and Yuanzhi Li. Faster Principal Component Regression via Optimal Polynomial Approximation to sgn(x). ArXiv e-prints, abs/1608.04773, August 2016.
[5] Zeyuan Allen-Zhu and Lorenzo Orecchia. Linear Coupling: An Ultimate Unification of Gradient and Mirror Descent. In Proceedings of the 8th Innovations in Theoretical Computer
Science, ITCS ’17, 2017. Full version available at http://arxiv.org/abs/1407.1537.
[6] Zeyuan Allen-Zhu, Peter Richtárik, Zheng Qu, and Yang Yuan. Even faster accelerated coordinate descent using non-uniform sampling. In ICML, 2016.
[7] Zeyuan Allen-Zhu and Yang Yuan. Improved SVRG for Non-Strongly-Convex or Sum-of-NonConvex Objectives. In ICML, 2016.
[8] Sanjeev Arora, Satish Rao, and Umesh V. Vazirani. Expander flows, geometric embeddings
and graph partitioning. Journal of the ACM, 56(2), 2009.
[9] Owe Axelsson. A survey of preconditioned iterative methods for linear systems of algebraic
equations. BIT Numerical Mathematics, 25(1):165–187, 1985.
[10] Kamalika Chaudhuri, Sham M Kakade, Karen Livescu, and Karthik Sridharan. Multi-view
clustering via canonical correlation analysis. In ICML, pages 129–136, 2009.
[11] Paramveer Dhillon, Dean P Foster, and Lyle H Ungar. Multi-view learning of word embeddings
via cca. In NIPS, pages 199–207, 2011.
[12] Roy Frostig, Cameron Musco, Christopher Musco, and Aaron Sidford. Principal Component
Projection Without Principal Component Analysis. In ICML, 2016.
[13] Dan Garber and Elad Hazan. Fast and simple PCA via convex optimization. ArXiv e-prints,
September 2015.
[14] Dan Garber, Elad Hazan, Chi Jin, Sham M. Kakade, Cameron Musco, Praneeth Netrapalli,
and Aaron Sidford. Robust shift-and-invert preconditioning: Faster and more sample efficient
algorithms for eigenvector computation. In ICML, 2016.
[15] Rong Ge, Chi Jin, Sham M. Kakade, Praneeth Netrapalli, and Aaron Sidford. Efficient Algorithms for Large-scale Generalized Eigenvector Computation and Canonical Correlation Analysis. In ICML, 2016.
[16] Gene H. Golub and Charles F. Van Loan. Matrix Computations. The JHU Press, 4th edition,
2012.
[17] Sham M Kakade and Dean P Foster. Multi-view regression via canonical correlation analysis.
In Learning theory, pages 82–96. Springer, 2007.
[18] Nikos Karampatziakis and Paul Mineiro. Discriminative features via generalized eigenvectors.
In ICML, pages 494–502, 2014.
34

[19] Yichao Lu and Dean P Foster. Large scale canonical correlation analysis with iterative least
squares. In NIPS, pages 91–99, 2014.
[20] Zhuang Ma, Yichao Lu, and Dean Foster. Finding linear structure in large datasets with
scalable canonical correlation analysis. In ICML, pages 169–178, 2015.
[21] Tomer Michaeli, Weiran Wang, and Karen Livescu. Nonparametric canonical correlation analysis. arXiv preprint, abs/1511.04839, 2015.
[22] Cameron Musco and Christopher Musco. Randomized block krylov methods for stronger and
faster approximate singular value decomposition. In NIPS, pages 1396–1404, 2015.
[23] Yurii Nesterov. A method of solving a convex programming problem with convergence rate
O(1/k 2 ). In Doklady AN SSSR (translated as Soviet Mathematics Doklady), volume 269, pages
543–547, 1983.
[24] Shai Shalev-Shwartz. SDCA without Duality, Regularization, and Individual Convexity. In
ICML, 2016.
[25] Shai Shalev-Shwartz and Tong Zhang. Accelerated Proximal Stochastic Dual Coordinate Ascent for Regularized Loss Minimization. In Proceedings of the 31st International Conference
on Machine Learning, ICML 2014, pages 64–72, 2014.
[26] Ohad Shamir. A Stochastic PCA and SVD Algorithm with an Exponential Convergence Rate.
In ICML, pages 144—-153, 2015.
[27] Jonathan Richard Shewchuk. An introduction to the conjugate gradient method without the
agonizing pain, 1994.
[28] Weiran Wang and Karen Livescu. Large-scale approximate kernel canonical correlation analysis. arXiv preprint, abs/1511.04773, 2015.
[29] Weiran Wang, Jialei Wang, Dan Garber, and Nathan Srebro. Efficient Globally Convergent
Stochastic Optimization for Canonical Correlation Analysis. In NIPS, 2016.
[30] Daniela M Witten, Robert Tibshirani, and Trevor Hastie. A penalized matrix decomposition,
with applications to sparse principal components and canonical correlation analysis. Biostatistics, page kxp008, 2009.

35

