arXiv:1504.06729v4 [cs.DS] 12 Jul 2016

Optimal Principal Component Analysis in Distributed and
Streaming Models
Christos Boutsidis
New York, New York
christos.boutsidis@gmail.com

David P. Woodruff
IBM Research
Almaden, California
dpwoodru@us.ibm.com

Peilin Zhong
Institute for Interdisciplinary Information Sciences
Tsinghua University, Beijing, China
zpl12@mails.tsinghua.edu.cn
Abstract
We study the Principal Component Analysis (PCA) problem in the distributed and streaming
models of computation. Given a matrix A ∈ Rm×n , a rank parameter k < rank(A), and an
accuracy parameter 0 < ε < 1, we want to output an m × k orthonormal matrix U for which
kA − UUT Ak2F ≤ (1 + ε) · kA − Ak k2F ,

where Ak ∈ Rm×n is the best rank-k approximation to A. We show the following.

1. In the arbitrary partition
Psmodel of Kannan et al. (COLT 2014), each of s machines holds
a matrix Ai and A = i=1 Ai . Each machine should output U. Kannan et al. achieve
O(skm/ε) + poly(sk/ε) words (of O(log(nm)) bits) communication. We obtain the improved bound of O(skm) + poly(sk/ε) words, and show an optimal Ω(skm) lower bound
up to low order terms. This resolves an open question for high precision PCA. A poly(ε−1 )
dependence is known to be required, but we separate this dependence from m.
2. We bypass the above lower bound when A is φ-sparse in each column and each server
receives a subset of columns. Here we obtain an O(skφ/ε) + poly(sk/ε) word protocol.
Our communication is independent of the matrix dimensions, and achieves the guarantee
that each server, in addition to outputting U, outputs a subset of O(k/ε) columns of A
containing a U in its span (that is, we solve distributed column subset selection). We show
a matching Ω(skφ/ε) lower bound for distributed column subset selection. Achieving our
communication bound when A is sparse but not sparse in each column, is impossible.
3. In the streaming model in which the columns arrive one at a time, an algorithm of Liberty
(KDD, 2013) with an improved analysis by Ghashami and Phillips (SODA, 2014) shows
O(km/ε) “real numbers” of space is achievable in a single pass, which we first improve to
an O(km/ε)+poly(k/ε) word space upper bound. This almost matches a known Ω(km/ε)
bit lower bound of Woodruff (NIPS, 2014). We show with two passes one can achieve
O(km) + poly(k/ε) words of space and (up to the poly(k/ε) term and the distinction
between words versus bits) this is optimal for any constant number of passes.
4. In turnstile streams, in which we receive entries of A one at a time in an arbitrary order, we show how to obtain a factorization of a (1 + ε)-approximate rank-k matrix using
O((m+n)kε−1 ) words of space. This improves the O((m+nε−2 )kε−2 ) bound of Clarkson
and Woodruff (STOC 2009), and matches their Ω((m + n)kε−1 ) word lower bound.

Notably, our results do not depend on the condition number of A.

1 Introduction
In distributed-memory computing systems, such as, Hadoop [1] or Spark [3], Principal Component Analysis (PCA) and the related Singular Value Decomposition (SVD) of large matrices is
becoming very challenging. Machine learning libraries implemented on top of such systems, for
example mahout [2] or mllib [4], provide distributed PCA implementations since PCA is often
used as a building block for a learning algorithm. PCA is useful for dimension reduction, noise
removal, visualization, etc. In all of these implementations, the bottleneck is the communication;
hence, the focus has been on minimizing the communication cost of the related algorithms, and
not the computational cost, which is the bottleneck in more traditional batch systems.
The data matrix corresponding to the dataset in hand, e.g., a term-document matrix representing a text collection, or the Netflix matrix representing user’s ratings for different movies, could
be distributed in many different ways [46]. In this paper, we focus on the following so-called arbitrary partition and column partition models. In the arbitrary partition
model, a matrix A ∈ Rm×n
Ps
is arbitrarily distributed among s machines. Specifically, A = i=1 Ai where Ai ∈ Rm×n is held
by machine i. Unless otherwise stated, we always assume m ≤ n. In the column partition model,

A
A
.
.
.
A
a matrix A ∈ Rm×n is distributed column-wise among s < n machines:
A
=
;
1
2
s
P
here, for i = 1 : s, Ai is an m × wi column submatrix of A with i wi = n. Note that the column partition model is a special case of the arbitrary partition model. Thus, it is desirable to
prove upper bounds in the arbitrary partition model, and lower bounds in the column partition
model. Both models have been adapted by traditional numerical linear algebra software libraries
for distributed memory matrix computations [12].
A recent body of work [29, 10, 41, 31, 37, 13] (see Section 3 for a comparison) has focused
on designing algorithms which minimize the communication needed for each machine to output an m × k matrix U for which an m × k orthonormal matrix U for which kA − UUT Ak2F ≤
(1 + ε) · kA − Ak k2F , where Ak ∈ Rm×n is the best rank-k approximation to A. Each machine
should output the same matrix U which can then be used for downstream applications, such as
clustering, where one first projects the data to a lower dimensional subspace (see, e.g., [10]). The
protocol should succeed with large constant probability. The model is such that each of the s machines communicates with one machine, called the “coordinator” (which we sometimes refer to as
Server), but does not communicate with other machines. This is known as the coordinator model
and one can simulate arbitrary point-to-point communication with a constant factor overhead in
communication together with an additive O(log(mn)) bit overhead per message (see, e.g., [45]).

1.1 Our Results
The best upper bound in the arbitrary partition model is O((skmε−1 ) + s · poly(kε−1 )) words (of
O(log(mn)) bits) communication [37], while the only known lower bound for these problems is in
the arbitrary partition model and is Ω(skm) [37]. In high precision applications one may want to
set ε to be as small as possible and the leading order term of O(skmε−1 ) is undesirable. A natural
question is whether there is an O((skm) + s · poly(kε−1 )) word protocol. We note that there is a
fairly simple known Ω(ε−2 ) bit lower bound [54] (we present this in Lemma 75 for completeness),
so although one cannot achieve an O(log(1/ε)) dependence in the communication as one would
maybe expect given iterative algorithms for regression with running time O(log(1/ε)) (see section
2.6 of [55] for an overview), an O((skm) + s · poly(kε−1 )) would at least separate the dependence
of ε and m and allowing for much smaller ε with the same amount of communication.
There are many existing protocols using O(skmε−1 ) words of communication [29, 10, 41, 31,
37, 13], and the protocols work in very different ways: that of [29, 10] is based on coresets, while
1

that of [41, 31] is based on adapting a streaming algorithm to the communication setting, while
that of [37] is based on sketching, and that of [13] is based on alternating minimization and is
useful only under some assumptions on the condition number (see Section 3 for a more detailed
discussion of these protocols). It was thus natural to assume there should be a lower bound of
Ω(skmε−1 ).
Instead, we obtain a new upper bound in the arbitrary partition model and a matching lower
bound (up to lower order terms) in the column partition model.
Theorem 1 (Restatement of Theorem 43 and Theorem 82). Suppose matrix A ∈ Rm×n is partitioned
in the arbitrary-partition model (See Definition 14). For any 1 ≥ ε > 0, there is an algorithm which on
termination leaves the same orthonormal matrix U ∈ Rm×k on each machine such that kA − UUT Ak2F ≤
(1 + ε) · kA − Ak k2F holds with arbitrarily large constant probability. Further, the algorithm runs in
polynomial time with total communication complexity O(skm + s · poly(kε−1 )) words each containing
O(log(smnε−1 )) bits.
If matrix A ∈ Rm×n is distributed in the column-partition model (See Definition 15), then for any
positive constant C ≤ O(poly(skm)), the algorithm which can leave a orthonormal matrix U ∈ Rm×k on
each machine such that kA − UUT Ak2F ≤ C · kA − Ak k2F holds with constant probability has at least
Ω(skm log(skm)) bits of communication.
In some applications even an O(skm) word protocol may be too costly, as m could be very
large. One could hope to do for communication what recent work on input sparsity algorithms
[22, 43, 44, 16, 25] has done for computation, i.e., obtain protocols sensitive to the number of
non-zero entries of A. Many matrices are sparse, e.g., Netflix provides a training data set of
100, 480, 507 ratings that 480, 189 users give to 17, 770 movies. If users correspond to columns
in the matrix, then the average column sparsity is ≈ 200 non-zero elements (out of the 17, 770
possible coordinates).
Our second contribution is the first protocol depending on the number of non-zero entries of
A in the column partition model. Our communication is independent of the matrix dimensions.
Denote by φ the maximum number of non-zero elements of a column in A. When we say that A
is φ-sparse, we mean that every column of A has at most φ non-zero elements and nnz(A) ≤ φ · n.
Our protocol has the additional feature of leaving on each machine the same subset C of O(k/ε)
columns of A for which there exists an m × k orthonormal matrix U in the column span of C for
which kA − CC† Ak2F ≤ kA − UUT Ak2F ≤ (1 + ε) · kA − Ak k2F . This is known as the column subset
selection problem, which is useful since C may be sparse if A is, and also it can lead to better data
intepretability. To partially complement our protocol, we show our protocol is optimal for any
protocol solving the column subset selection problem. We summarize these results as follows.
Theorem 2 (Restatement of Theorem 63). Suppose a φ-sparse matrix A ∈ Rm×n is partitioned in
the column-partition model (See Definition 15). For any 1 ≥ ε > 0, there is an algorithm which on
termination leaves C ∈ Rm×c with c = O(k/ε) columns of A and an orthonormal matrix U on each
machine such that kA − CC† Ak2F ≤ kA − UUT Ak2F ≤ (1 + ε) · kA − Ak k2F holds with arbitrarily large
constant probability. Further,
the algorithm runs in polynomial time with total communication complexity

O skφε−1 + sk2 ε−4 words each containing of O(log(smnε−1 )) bits.
If a φ-sparse matrix A ∈ Rm×n is distributed in the column-partition model (See Definition 15),
then for any positive constant C ≤ O(poly(skφ)), the algorithm which can leave a orthonormal matrix
U ∈ Rm×k on each machine such that kA − UUT Ak2F ≤ C · kA − Ak k2F holds with constant probability
has at least Ω(skφ log(skφ)) bits of communication.
Remark 3. Our protocol has no dependence on condition number and works for arbitrary matrices. We
2

Upper bounds

Lower bounds

Definition 16 (arbitrary partition model)

O(skm + s · poly(k/ε))

(Theorem 43)

Ω(skm) (Theorem 1.2 in [37])

Definition 17 (column partition model)

O(skm + s · poly(k/ε))

(Theorem 43)

Ω(skm) (Theorem 82)

Definition 17 with sparsity φ = o(ε · m)

O(skφε−1 + s · poly(k/ε)) (Theorem 63)

Ω(skφ) (Corollary 83 )

Table 1: Communication upper/lower bounds for the Distributed PCA problems.

make a serious effort as shown in Section 6 to achieve this. Note that even if the entries of A are specified
using O(log(mn)) bits, the condition number can be 2−Ω(mn log(mn)) [7].
Remark 4. We did not discuss the running time of the algorithms, but we are interested in the fastest
possible algorithms with the minimal communication cost, and provide fast algorithms as well.
Remark 5. The hard instance in Section 10.2.2 implies an Ω(skm) lower bound even if the input matrix
is sparse overall but has skm non-zero positions in arbitrarily locations. Therefore, to obtain our bounds we
discuss the sparsity on each column instead of the overall sparsity of the matrix.
A model closely related to the distributed model of computation is the streaming model of
computation. The model we focused on is the so-called turnstile streaming model. In this model,
there is a stream of update operations and each operation indicates that the corresponding entry of
A should be incremented by a specific number. We present novel PCA algorithms in this model.
Our first one-pass algorithm improves upon the best existing streaming PCA algorithm [41, 31]
in two respects. First, the space of [41, 31] is described in “real numbers” while our space bound
(O(mk/ε + poly(k/ε)) - see Theorem 52) is in terms of words (we also bound the word size). This
matches an Ω(km/ε) bit lower bound for one-pass algorithms in [53], up to the distinction between
words versus bits and a low order term poly(k/ε). Second, our algorithm can be applied in the
turnstile streaming model which is stronger than column update streaming model in [41, 31].
Theorem 6 (Restatement of Theorem 52). Suppose A ∈ Rm×n is given by a stream of update operations
in the turnstile streaming model (See Definition 20). For any 1 ≥ ε > 0, there is an algorithm which
uses a single pass over the stream and on termination outputs an orthonormal matrix U ∈ Rm×k such
that kA − UUT Ak2F ≤ (1 + ε) · kA − Ak k2F holds with arbitrarily large constant
 probability. Further,
the algorithm runs in polynomial time with space of total O mk/ε + poly(kε−1 ) words each containing
O(log(smnε−1 )) bits.
A slight modification of the previous algorithm leads to a one-pass algorithm which can compute a factorization of a (1 + ε)-approximate rank-k matrix. The modified algorithm only needs
O((n + m)k/ε + poly(k/ε)) words of space which improves the O((m + nε−2 )kε−2 ) upper bound
in [21] and matches the Ω((n + m)k/ε) bit lower bound given by [21], up to the low order term
poly(k/ε).
Theorem 7 (Restatement of Theorem 53). Suppose A ∈ Rm×n is given by a stream of update operations
in the turnstile streaming model (See Definition 20). For any 1 ≥ ε > 0, there is an algorithm which uses
a single pass over the stream and on termination outputs a matrix A∗k ∈ Rm×n with rank(A∗k ) ≤ k such
that kA − A∗k k2F ≤ (1 + ε) · kA − Ak k2F holds with arbitrarily large constant probability. Further, the
algorithm runs in polynomial time with space of total O((m + n)k/ε + poly(kε−1 )) words each containing
O(log(smnε−1 )) bits.
3

Upper bounds

One-pass turnstile, Def. 21
One-pass turnstile, factorization, Def. 22

Two-pass turnstile, Def. 21

O(mkε−1 + poly(k, ε−1 ))

Lower bounds

(Theorem 52)

O((n + m)kε−1 + poly(k, ε−1 )) (Theorem 53)

O(mk + poly(k, ε−1 ))

(Theorem 55)

Ω(mkε−1 ) bits

[53]

Ω((n + m)kε−1 ) words [21]

Ω(mk) bits

[53]

Table 2: Space upper/lower bounds for the Streaming PCA problem.

We also show a two-pass algorithm which is an implementation of our distributed PCA protocol. It uses O((km)+poly(k/ε)) words of space, which up to the poly(k/ε) term and the distinction
between words versus bits, is optimal for any constant number of passes. A “next natural goal”
in [53] was to improve the lower bound of Ω(km) to a bound closer to the 1-pass Ω(km/ε) lower
bound established in that paper; our upper bound shows this is not possible.
Theorem 8 (Restatement of Theorem 55). Suppose A ∈ Rm×n is given by a stream of update operations
in the turnstile streaming model (See Definition 20). For any 1 ≥ ε > 0, there is an algorithm which uses
two passes over the stream and on termination outputs an orthonormal matrix U ∈ Rm×k such that kA −
UUT Ak2F ≤ (1 + ε) · kA − Ak k2F holds with arbitrarily large constant
probability. Further, the algorithm

runs in polynomial time with total space O mk + poly(kε−1 ) words each containing O(log(smnε−1 ))
bits.

1.2 Technical Overview
1.2.1

Upper Bounds

The Arbitrary Partition Model. Recall that the inputP
matrix A ∈ Rm×n is arbitrarily partitioned
m×n
into s matrices Ai ∈ R
, i.e., for i = 1, 2, . . . , s: A = si=1 Ai . By sketching on the left and right
2
2
by two small random sign matrices S ∈ RO(k/ε )×m and T ∈ Rn×O(k/ε ) , with O(sk 2 /ε4 ) words of
communication, the server can learn Ã = SAT which is a “sketch” of A. It suffices to compute
T
the best rank-k approximation Ãk to Ã. Suppose the SVD of Ãk = UÃk ΣÃk VÃ
. Then the server
k
can learn ATVÃk by the second round of communication which needs O(skm) words. We then
prove (see Lemma 28): kA − UUT Ak2F ≤ (1 + ε)kA − Ak k2F where U is an orthonormal basis of the
column space of ATVÃk . Notice that rank(U) ≤ k. Thus U is exactly what we want. Section 5.3
shows the details.
A technical obstacle in the analysis above is that it may require a large number of machine
words to specify the entries of VÃk , even if each entry in A is only a single word. In fact, one can
show (we omit the details) that the singular values of Ãk can be exponentially large in k/ε, which
means one would need to round the entries of VÃk to an additive exponentially small precision,
which would translate to an undesirable sm · poly(k/ε) bits of communication.
To counter this, we use the smoothed analysis of Tao and Vu [52] to argue that if we add small
random Bernoulli noise to A, then its minimum singular value becomes inverse polynomial in
n. Moreover, if the rank of A is at least 2k and under the assumption that the entries of A are
representable with a single machine word (so we can identify them with integers in magnitude
at most poly(nms/ε) by scaling), then we can show the additional noise preserves relative error
approximation. To our knowledge, this is the first application of this smoothing technique to
distributed linear algebra algorithms. Section 6 discusses the details of this approach.
4

On the other hand, if the rank of A is smaller than 2k, the smoothing technique need not preserve relative error. In this case though, we can learn a basis C ∈ Rm×O(k) for the column span
of A by multiplying by a pseudorandom matrix based on Vandermonde matrices. Since C has at
most 2k columns, we can efficiently communicate it to all servers using O(skm) communication.
At this point we set up the optimization problem minrank(X)≤k kCXCT A−AkF . The server cannot
learn CT A and A directly, as this would be Ω(nm) words, but we can choose additional “sketching matrices” Tlef t and Tright to instead solve minrank(X)≤k kTlef t (CXCT A − A)Tright kF , which
is a generalization of the subspace embedding technique to “affine embeddings” [53]. the server
learns Tlef t CXCT ATright and Tlef t ATright , which are small matrices, and can be sent to each machine. Each machine locally solves for the rank-k solution X∗ minimizing kTlef t CXCT ATright −
Tlef t ATright kF . Now each machine can find the same m × k orthonormal basis and output it
without further communication.
Sparse Matrices in the Column Partition Model. Our algorithm for sparse matrices is technically simpler than that in the arbitrary partition model. Section 8.2 presents a distributed PCA
algorithm that has O((φksε−1 ) + poly(sk/ε)) words of communication cost. Our idea in order to
take advantage of the sparsity in the input matrix A is to select and transfer to the coordinator
a small set of “good” columns from each sub-matrix Ai . Specifically, we design an algorithm
that first computes, in a distributed way, a matrix C̃ ∈ Rm×c with c = O(kε−1 ) columns of A, and
then finds U ∈ span(C̃) using a distributed, communication-efficient algorithm developed in [37].
The matrix C̃ is constructed using optimal algorithms for column sampling [17, 26, 24], extended
properly to the distributed case.
As mentioned, our algorithm actually solves the distributed column subset selection problem.
It builds on results for column subset selection in the batch setting [17, 26, 24], but gives a novel
analysis showing that in the distributed setting one can find columns as good as in the batch
setting. To the best of our knowledge, only heuristics for distributed column subset selection were
known [28]. We also optimize the time complexity of our algorithm, see Table 4.
Turnstile streaming model. Our one-pass streaming PCA algorithm starts by generating two
random sign matrices S ∈ RO(k/ε)×n and R ∈ Rm×O(k/ε) . As shown in Lemma 51,
min
kARXSA − Ak2F ≤ (1 + ε) · kA − Ak k2F
rank(X)≤k
T
The intuition of the above is as the following. Suppose the SVD of Ak is UAk ΣAk VA
where
k
Ak is the best rank-k approximation matrix to A. Since X can be chosen as ΣAk VAk , we have
minrank(X)≤k kUAk X − Ak2F = kAk − Ak2F . As shown in [21], S can provide a sketch for the
regression problem. Thus, kUAk X̃ − Ak2F ≤ (1 + ε) · minrank(X)≤k kUAk X − Ak2F where

X̃ = arg

min

rank(X)≤k

kSUAk X − SAk2F

Notice that X̃ is in the row space of span(SA). Then we focus on the regression problem:
min
kXSA − Ak2F .
rank(X)≤k
Similarly, we can use R to sketch the problem. Thus finally we only need to optimize
min
kARXSA − Ak2F .
rank(X)≤k
5

Upper bounds

Definition 18

Definition 19

O(skφε−1 )

Lower bounds

(Theorem 74)

Ω(skφε−1 ) (Theorem 90)



O skφε−1 + s · poly(k/ε) (Theorem 74)

Ω(skφε−1 ) (Corollary 91)

Table 3: Communication upper/lower bounds for the CSSP problems of Definitions 18 and 19.

However, it is too expensive to store the entire matrix A. Thus, we sketch on the left and right
by Tlef t and Tright . Since all of Tlef t AR, SATright and Tlef t ATright can be maintained in small
space, X∗ = arg minrank(X)≤k kTlef t (ARXSA−A)Tright k2F can be constructed after one pass over

the stream. Additionally, if AR is also maintained in the algorithm, we can get U ∈ Rm×k of which
columns are an orthonormal basis of span(ARX∗ ) satisfying kA − UUT Ak2F ≤ (1+ ε)·kA − Ak k2F .
T
,
Furthermore, if SA is also maintained, since it suffices to compute the SVD of X∗ = UX∗ ΣX∗ VX
∗
T
T = Tlef t ARUX∗ and K = VX
SAT
can
be
constructed
in
O((n
+
m)k)
words
of
space.
right
∗
Therefore, the algorithm can compute A∗k = TΣX∗ K = ARX∗ SA in O((n + m)k) words of
space. The algorithm then can output A∗k with total space O((n + m)k/ε + poly(k/ε)) words (see
Theorem 53).
Our two-pass streaming PCA algorithm is just an implemententation of our distributed PCA
algorithm in the streaming model.
1.2.2

Lower Bounds

Table 3 summarizes the matching communication lower bounds that we obtain for distributed
column-based matrix reconstruction. Theorem 90 proves a lower bound for the problem in Definition 18; then, a lower bound for the problem of Definition 19 follows immediately since this is a
harder problem (see Corollalry 91).
Distributed Column Subset Selection. Our most involved lower bound argument is in showing
the tightness for the distributed column subset selection problem, and is given in Theorem 90. To
illustrate the ideas, suppose k = 1. We start with the canonical hard matrix for column subset
selection, namely, an m × m matrix A whose first row is all ones, and remaining rows are a subset
of the identity matrix. Intuitively the best rank-1 approximation is very aligned with the first row
of A. However, given only o(1/ε) columns of the matrix, there is not a vector in the span which
puts a (1 − ε)-fraction of its mass on the first coordinate, so Ω(1/ε) columns are needed. Such a
matrix has φ = 2.
Obtaining a lower bound for general φ involves creating a large net of such instances. Suppose
we set m = φ and consider L̃A, where L̃ is formed by choosing a random φ × φ orthonormal
matrix, and then rounding each entry to the nearest integer multiple of 1/poly(n). We can show
that o(1/ε) columns of L̃A do not span a rank-1 approximation to L̃A. We also need a lemma
which shows that if we take two independently random orthonormal matrices, and round their
entries to integer multiples of 1/poly(n), then these matrices are extremely unlikely to share a
constant fraction of columns. The idea then is that if one server holds L̃A, and every other server
holds the zero matrix, then every server needs to output the same subset of columns of L̃A. By
construction of A, each column of L̃A is the sum of the first column of L̃ and an arbitrary column
of L̃, and since we can assume all servers know the first column of L̃ (which involves a negligible

6

O(sφ) amount of communication), this implies that each server must learn Ω(1/ε) columns of L̃.
The probability that random discretized orthonormal matrices share Ω(1/ε) columns is sufficiently
small that we can choose a large enough net for these columns to identify L̃ in that net, requiring
Ω(φ/ε) bits of communication per server, or Ω(sφ/ε) words in total. The argument for general k
can be viewed as a “direct sum theorem”, giving Ω(skφ/ε) words of communication in total.
Dense Matrices in the Column Partition Model. Our optimal lower bound for dense matrices
is technically simpler. It gives an Ω(skm) word communication lower bound in column partition
model for dense matrices. Theorem 82 argues that there exists an m × n matrix A such that for this
A, any k ≤ 0.99m, and any error parameter C with 1 < C < poly(skm), if there exists a protocol to
construct an m × k matrix U satisfying kA − UUT Ak2F ≤ C · kA − Ak k2F , with constant probability,
then this protocol has communication cost at least Ω(skm) words. This lower bound is proven for
a matrix A that has k fully dense columns. We should note that the lower bound holds even if
each machine only outputs the projection of Ai onto U, i.e., UUT Ai , rather than U itself. Note
that this problem is potentially easier, since given U, a machine could find UUT Ai .
The intuition of the proof is that machine 1 has a random matrix A1 ∈ Rm×k with orthonormal
columns (rounded to integer multiples of 1/poly(n)), while all other machines have a very small
multiple of the identity matrix. When concatenating the columns of these matrices, the best kdimensional subspace to project the columns onto must be very close to A1 in spectral norm.
Since machine i, i > 1, has tiny identity matrices, after projection it obtains a projection matrix
very close to the projection onto the column span of A1 . By choosing a net of m × k orthonormal
matrices, all of which pairwise have high distance in spectral norm, each machine can reconstruct
a random element in this net, which requires lots of information, and hence communication.
Our lower bound of Theorem 82 implies a lower bound of Ω(skφ) words for matrices in which
the column sparsity is φ (see Corollary 83) as a simple corollary.

1.3 Road Map
In the first ten pages, we have one additional section, Section 2, which gives an overview of our
algorithm in the arbitrary partition model.
In the Appendix, we then present more detail on prior results on distributed PCA algorithms
in Section 3. Section 4 introduces the notation and basic results from linear algebra that we use
in our algorithms. Section 5 presents our distributed PCA algorithm for arbitrary matrices for
Definition 16. The communication of the algorithm in this section can be stated only in terms of
“real numbers”. We resolve this issue in Section 6 where a modified algorithm has communication
cost bounded in terms of machine words. Section 7 discusses space-optimal PCA methods in
the streaming model of computation. Section 8 presents a distributed PCA algorithm for sparse
matrices in column partition model, while Section 9 extends this algorithm to a faster distributed
PCA algorithm for sparse matrices. Section 10 presents communication lower bounds.

2 Outline of our result in the arbitrary partition model
Here we give a brief outline of our result in the arbitrary patition model, deferring the full details
to the appendix.
We describe a fast distributed PCA algorithm with total communication O(msk) words plus
low order terms, which is optimal in the arbitrary partition model in the sense that an Ω̃(msk) bit

7

lower bound was given by [37]. The algorithm employs, in a novel way, the notion of projectioncost preserving sketches from [24]. In particular, whereas all previous [48, 23] dimension-reductionbased SVD methods reduce one dimension of the input matrix to compute some approximation to
the SVD, our method reduces both dimensions and computes an approximation to the SVD from a
small almost square matrix (we note that other work, such as work on estimating eigenvalues [8],
uses sketches in both dimensions, but it is not clear how to compute singular vectors using such
sketches). Unlike [37] which reduces only one dimension in the first communication round, we
do the reduction on both dimensions in the same round.
We first present the batch version of the algorithm which offers a new low-rank matrix approximation technique; a specific implementation of this algorithm offers a communication-optimal
distributed PCA algorithm (we also discuss in Section 7 a variant of this algorithm that offers
a two-pass space-optimal PCA method in the turnstile streaming model). Before presenting all
these new algorithms in detail, we present the relevant results from the previous literature that
we employ in the analysis.

2.1 Projection-cost preserving sketching matrices
In this section, we recap a notion of sketching matrices which we call “projection-cost preserving
sketching matrices”. A sketching matrix from this family is a linear matrix transformation and it
has the property that for all projections it preserves, up to some error, the difference between the
matrix in hand and its projection in Frobenius norm.
Definition 9 (Projection-cost preserving sketching matrices). We say that W ∈ Rn×ξ is an (ε, k)projection-cost preserving sketching matrix of A ∈ Rm×n , if for all rank-k orthogonal projection matrices
P ∈ Rm×m , it satisfies
(1 − ε)kA − PAk2F ≤ kAW − PAWk2F + c ≤ (1 + ε)kA − PAk2F
where c is a non-negative constant which only depends on A and W. We also call AW an (ε, k)-projectioncost preserving sketch of A.
Due to the following lemma, we know that a good rank-k approximation projection matrix of
(ε, k)-projection-cost preserving sketch AW also provides a good rank-k approximation to A.
Lemma 10 (PCA via Projection-Cost Preserving Sketches - Lemma 3 in [24]). Suppose W ∈ Rn×ξ is
∗
an (ε, k)-projection-cost preserving sketching matrix of A ∈ Rm×n . Let P̂ = arg minrank(P)≤k kAW −
∗
PAWk2F . For all P̂, ε′ satisfying rank(P̂) ≤ k, ε′ ≥ 0, if kAW− P̂AWk2F ≤ (1+ε′ )kAW− P̂ AWk2F ,
kA − P̂Ak2F ≤

1+ε
· (1 + ε′ )kA − Ak k2F
1−ε

[24] also provides several ways to construct projection-cost preserving sketching matrices. Because we mainly consider the communication, we just choose one which can reduce the dimension
as much as possible. Furthermore, it is also an oblivious projection-cost preserving sketching matrix.
Lemma 11 (Dense Johnson-Lindenstrauss matrix - part of Theorem 12 in [24]). √
For ε < √
1, suppose
n×ξ
each entry of W ∈ R
is chosen O(log(k))-wise independently and uniformly in {1/ ξ, −1/ ξ} where
ξ = O(kε−2 ) [21]. For any A ∈ Rm×n , with probability at least 0.99, W is an (ε, k)-projection-cost
preserving sketching matrix of A.
8

2.2 A batch algorithm for the fast low rank approximation of matrices
In this section, we describe a new method for quickly computing a low-rank approximation to
a given matrix. This method does not offer any specific advantages over previous such techniques [48, 23, 18]; however, this new algorithm can be implemented efficiently in the distributed
setting (see Section 5.3) and in the streaming model of computation (see Section 7); in fact we
are able to obtain communication-optimal and space-optimal results, respectively. For completeness as well as ease of presentation, we first present and analyze the simple batch version of the
algorithm. The algorithm uses the dense Johnson-Lindenstrauss matrix of Lemma 27 in order
to reduce both dimensions of A, before computing some sort of SVD to a poly(k/ε) × poly(k/ε)
matrix (see Step 2 in the algorithm below).
Consider the usual inputs: a matrix A ∈ Rm×n , a rank parameter k < rank(A), and an accuracy parameter 0 < ε < 1. The algorithm below returns an orthonormal matrix U ∈ Rm×k such
that
kA − UUT Ak2F ≤ (1 + ε) · kA − Ak k2F .
Algorithm
1. Construct two dense Johnson-Lindenstrauss matrices S ∈ Rξ1 ×m , T ∈ Rn×ξ2 with ξ1 = O(kε−2 ), ξ2 =
O(kε−2 ) (see Lemma 27).
2. Construct Ã = SAT.
T
3. Compute the SVD of Ãk = UÃk ΣÃk VÃ
(UÃk ∈ Rξ1 ×k , ΣÃk ∈ Rk×k , VÃk ∈ Rξ2 ×k ).
k

4. Construct X = ATVÃk
5. Compute an orthonormal basis U ∈ Rm×k for span(X) (notice that rank(X) ≤ k).

Theorem 29 later in this section analyzes the approximation error and the running time of the
previous algorithm. First, we prove the accuracy of the algorithm.
Lemma 12. The matrix U ∈ Rm×k with k orthonormal columns satisfies with probability at least 0.98:
kA − UUT Ak2F ≤ (1 + ε) · kA − Ak k2F .
Proof.
T
T
T
kÃ − Ãk k2F = kÃ − ÃVÃk VÃ
k2F = kSAT − SATVÃk VÃ
k2F = kTT AT ST − VÃk VÃ
TT AT ST k2F
k

k

k

T
. The second equality is by the
The first equality follows by the SVD of Ãk = UÃk ΣÃk VÃ
k

construction of Ã = SAT. The third equality is due to ∀M, kMk2F = kM T k2F .
Due to Lemma 27, with probability at least 0.99, ST is an (ε, k)-projection-cost preserving
sketch matrix of TT AT . According to Lemma 26,
T
T
TT AT k2F ≤
k2F = kTT AT − VÃk VÃ
kAT − ATVÃk VÃ
k

k

1+ε
· kAT − (AT)k k2F
1−ε

Observe that
T
T
kUUT AT − ATk2F = kXX† AT − ATk2F ≤ kXVÃ
− ATk2F = kATVÃk VÃ
− ATk2F
k

which is at most

1+ε
1−ε

· kAT − (AT)k k2F .
9

k

(1)

The first equality uses the fact that U is an orthonormal basis of span(X). The first inequality is
followed by ∀X, M, N, kXX† M − Mk2F ≤ kXN − Mk2F . The second equality uses the construction
that X = ATVÃk . The second inequality follows by Eqn (2).
Due to Lemma 27, with probability at least 0.99, T is an (ε, k)-projection-cost preserving sketch
matrix of A. Due to Lemma 26,
kUUT A − Ak2F ≤

(1 + ε)2
· kA − Ak k2F
(1 − ε)2

Due to union bound, the probability that ST is an (ε, k)-projection-cost preserving sketch matrix
of TT AT and T is an (ε, k)-projection-cost preserving sketch matrix of A is at least 0.98. Note that
(1+ε)2
is 1 + O(ε) when ε is small enough, so we can adjust ε here by a constant factor to show the
(1−ε)2
statement.
Next, we present the main theorem.
Theorem 13. The matrix U ∈ Rm×k with k orthonormal columns satisfies with probability at least 0.98:
kA − UUT Ak2F ≤ (1 + ε) · kA − Ak k2F .

The running time of the algorithm is

O nmkε−2 + mk2 ε−4 + poly(kε−1 ) .

Proof. The correctness is shown by Lemma 28. The running time is analyzed in Section 5.

2.3 The distributed PCA algorithm
Ps
Recall that the input matrix A ∈ Rm×n is partitioned arbitrarily as: A =
i Ai for i = 1 : s,
Ai ∈ Rm×n . The idea in the algorithm below is to implement the algorithm in Section 5.2 in the
distributed setting.
Input:
1. A ∈ Rm×n arbitrarily partitioned A =

2. rank parameter k < rank(A)

Ps
i

Ai for i = 1 : s, Ai ∈ Rm×n .

3. accuracy parameter ε > 0
Algorithm
1. Machines agree upon two dense Johnson-Lindenstrauss matrices S ∈ Rξ1 ×m , T ∈ Rn×ξ2 with ξ1 =
O(kε−2 ), ξ2 = O(kε−2 ) (see Lemma 27).
P
2. Each machine locally computes Ãi = SAi T and sends Ãi to the server. Server constructs Ã = i Ãi .
T
3. Server computes the SVD of Ãk = UÃk ΣÃk VÃ
(UÃk ∈ Rξ1 ×k , ΣÃk ∈ Rk×k , VÃk ∈ Rξ2 ×k ).
k

4. Server sends VÃk to all machines.
5. Each machine construct Xi = Ai TVÃk and sends Xi to the server. Server constructs X =
6. Server computes an orthonormal basis U ∈ Rm×k for span(X) (notice that rank(X) ≤ k).

P

i

Xi .

7. Server sends U to each machine.

Notice that in the first step, S and T can be described using a random seed that is O(log(k))wise independent due to Lemma 27.
The major remaining challenge here is to obtain small bit complexity, since this protocol could
have very large bit complexity. See Section 5 in the Appendix for the full details of our algorithm.
10

References
[1] http://hadoop.apache.org/.
[2] https://mahout.apache.org.
[3] https://spark.apache.org/.
[4] https://spark.apache.org/mllib/.
[5] D. Achlioptas. Database-friendly random projections: Johnson-lindenstrauss with binary
coins. J. Comput. Syst. Sci., 66(4):671–687, 2003.
[6] N. Ailon and E. Liberty. Fast dimension reduction using rademacher series on dual bch codes.
In ACM-SIAM Symposium on Discrete Algorithms (SODA), 2008.
[7] N. Alon and V. H. Vu. Anti-hadamard matrices, coin weighing, threshold gates, and indecomposable hypergraphs. J. Comb. Theory, Ser. A, 79(1):133–160, 1997.
[8] A. Andoni and H. L. Nguyen. Eigenvalues of a matrix in the streaming model. In Proceedings
of the Twenty-Fourth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA 2013, New
Orleans, Louisiana, USA, January 6-8, 2013, pages 1729–1737, 2013.
[9] Z.-J. Bai, R. H. Chan, and F. T. Luk. Principal component analysis for distributed data sets
with updating. In Advanced Parallel Processing Technologies, pages 471–483. Springer, 2005.
[10] M.-F. Balcan, V. Kanchanapally, Y. Liang, and D. Woodruff. Improved distributed principal
component analysis. arXiv preprint arXiv:1408.5823, to appear in NIPS, 2014.
[11] J. Batson, D. Spielman, and N. Srivastava. Twice-ramanujan sparsifiers. In Proceedings of the
41st annual ACM symposium on Theory of computing, pages 255–262. ACM, 2009.
[12] A. R. Benson, D. F. Gleich, and J. Demmel. Direct qr factorizations for tall-and-skinny matrices
in mapreduce architectures. In Big Data, 2013 IEEE International Conference on, pages 264–272.
IEEE, 2013.
[13] S. Bhojanapalli, P. Jain, and S. Sanghavi. Tighter low-rank approximation via sampling the
leveraged element. http://uts.cc.utexas.edu/˜bsrinadh/main.pdf, to appear in
SODA, 2015.
[14] L. S. Blackford, J. Choi, A. Cleary, E. D’Azevedo, J. Demmel, I. Dhillon, J. Dongarra, S. Hammarling, G. Henry, A. Petitet, et al. ScaLAPACK users’ guide, volume 4. siam, 1997.
[15] Y.-A. L. Borgne, S. Raybaud, and G. Bontempi. Distributed principal component analysis for
wireless sensor networks. Sensors, 2008.
[16] J. Bourgain, S. Dirksen, and J. Nelson. Toward a unified theory of sparse dimensionality
reduction in euclidean space. In Proceedings of the Forty-Seventh Annual ACM on Symposium
on Theory of Computing, STOC 2015, Portland, OR, USA, June 14-17, 2015, pages 499–508, 2015.
[17] C. Boutsidis, P. Drineas, and M. Magdon-Ismail. Near optimal column based matrix reconstruction. SIAM Journal on Computing (SICOMP), 2013.

11

[18] C. Boutsidis and A. Gittens. Improved matrix algorithms via the subsampled randomized
hadamard transform. SIAM Journal on Matrix Analysis and Applications, 34(3):1301–1340, 2013.
[19] C. Boutsidis and D. P. Woodruff. Optimal cur matrix decompositions. STOC, 2014.
[20] A. Chakrabarti and O. Regev. An optimal lower bound on the communication complexity of
gap-hamming-distance. SIAM Journal on Computing, 41(5):1299–1317, 2012.
[21] K. Clarkson and D. Woodruff. Numerical linear algebra in the streaming model. In Proceedings of the 41st annual ACM symposium on Theory of computing (STOC), 2009.
[22] K. L. Clarkson and D. P. Woodruff. Low rank approximation and regression in input sparsity
time. In In STOC, 2013.
[23] K. L. Clarkson and D. P. Woodruff. Low rank approximation and regression in input sparsity
time. ArxiV report: http://arxiv.org/pdf/1207.6365v4.pdf, 2013.
[24] M. Cohen, S. Elder, C. Musco, C. Musco, and M. Persu. Dimensionality reduction for k-means
clustering and low rank approximation. arXiv preprint arXiv:1410.6801, 2014.
[25] M. B. Cohen. Simpler and tighter analysis of sparse oblivious subspace embeddings. In
SODA, 2016.
[26] A. Deshpande, L. Rademacher, S. Vempala, and G. Wang. Matrix approximation and projective clustering via volume sampling. Theory of Computing, 2(12):225–247, 2006.
[27] A. Deshpande and S. Vempala. Adaptive sampling and fast low-rank matrix approximation.
In RANDOM - APPROX, 2006.
[28] A. K. Farahat, A. Elgohary, A. Ghodsi, and M. S. Kamel. Distributed column subset selection
on mapreduce. In Data Mining (ICDM), 2013 IEEE 13th International Conference on, pages 171–
180. IEEE, 2013.
[29] D. Feldman, M. Schmidt, and C. Sohler. Turning big data into tiny data: Constant-size coresets for k-means, pca and projective clustering. In SODA, pages 1434–1453. SIAM, 2013.
[30] S. Friedland and A. Torokhti. Generalized rank-constrained matrix approximations. SIAM
Journal on Matrix Analysis and Applications, 29(2):656–659, 2007.
[31] M. Ghashami and J. Phillips. Relative errors for deterministic low-rank matrix approximations. In SODA, 2013.
[32] G. H. Golub and C. F. Van Loan. Matrix computations, volume 3. JHU Press, 2012.
[33] M. Gu and S. C. Eisenstat. A divide-and-conquer algorithm for the bidiagonal svd. SIAM
Journal on Matrix Analysis and Applications, 16(1):79–92, 1995.
[34] R. Horn and C. Johnson. Matrix Analysis. Cambridge University Press, New York, 1985.
[35] E. Jessup and D. Sorensen. A parallel algorithm for computing the singular value decomposition of a matrix. Siam Journal on Matrix Analysis and Applications, 15(2):530–548, 1994.
[36] R. Kannan, S. Vempala, and D. P. Woodruff. Nimble algorithms for cloud computing. arXiv
preprint arXiv:1304.3162, v3, 2013.
12

[37] R. Kannan, S. S. Vempala, and D. P. Woodruff. Principal component analysis and higher
correlations for distributed data. In Proceedings of The 27th Conference on Learning Theory,
pages 1040–1057, 2014.
[38] M. Kapralov and K. Talwar. On differentially private low rank approximation. In Proceedings
of the Twenty-Fourth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA 2013, New
Orleans, Louisiana, USA, January 6-8, 2013, pages 1395–1414, 2013.
[39] D. Kempe and F. McSherry. A decentralized algorithm for spectral analysis. In Proceedings of
the thirty-sixth annual ACM symposium on Theory of computing, pages 561–568. ACM, 2004.
[40] Y. Li, H. L. Nguyên, and D. P. Woodruff. On sketching matrix norms and the top singular
vector. In Proceedings of the Twenty-Fifth Annual ACM-SIAM Symposium on Discrete Algorithms,
pages 1562–1581. SIAM, 2014.
[41] E. Liberty. Simple and deterministic matrix sketching. In KDD, pages 581–588. ACM, 2013.
[42] S. V. Macua, P. Belanovic, and S. Zazo. Consensus-based distributed principal component
analysis in wireless sensor networks. In Signal Processing Advances in Wireless Communications
(SPAWC), 2010 IEEE Eleventh International Workshop on, pages 1–5. IEEE, 2010.
[43] X. Meng and M. W. Mahoney. Low-distortion subspace embeddings in input-sparsity time
and applications to robust linear regression. In Proceedings of the forty-fifth annual ACM symposium on Theory of computing, pages 91–100. ACM, 2013.
[44] J. Nelson and H. L. Nguyên. Osnap: Faster numerical linear algebra algorithms via sparser
subspace embeddings. In Foundations of Computer Science (FOCS), 2013 IEEE 54th Annual
Symposium on, pages 117–126. IEEE, 2013.
[45] J. M. Phillips, E. Verbin, and Q. Zhang. Lower bounds for number-in-hand multiparty communication complexity, made easy. In Proceedings of the Twenty-Third Annual ACM-SIAM
Symposium on Discrete Algorithms, SODA 2012, Kyoto, Japan, January 17-19, 2012, pages 486–
501, 2012.
[46] J. Poulson, B. Marker, R. A. van de Geijn, J. R. Hammond, and N. A. Romero. Elemental: A
new framework for distributed memory dense matrix computations. ACM Transactions on
Mathematical Software (TOMS), 39(2):13, 2013.
[47] Y. Qu, G. Ostrouchov, N. Samatova, and A. Geist. Principal component analysis for dimension reduction in massive distributed data sets. In Proceedings of IEEE International Conference
on Data Mining (ICDM), 2002.
[48] T. Sarlos. Improved approximation algorithms for large matrices via random projections. In
IEEE Symposium on Foundations of Computer Science (FOCS), 2006.
[49] K. C. Sou and A. Ranzer. On generalized matrix approximation problem in the spectral norm.
Linear Algebra and its Applications, 436(7):2331–2341, 2012.
[50] F. Tisseur and J. Dongarra. A parallel divide and conquer algorithm for the symmetric eigenvalue problem on distributed memory architectures. SIAM Journal on Scientific Computing,
20(6):2223–2236, 1999.

13

[51] M. D. Vose. A linear algorithm for generating random numbers with a given distribution.
Software Engineering, IEEE Transactions on, 17(9):972–975, 1991.
[52] V. H. Vu and T. Tao. The condition number of a randomly perturbed matrix. In Proceedings
of the 39th Annual ACM Symposium on Theory of Computing, San Diego, California, USA, June
11-13, 2007, pages 248–255, 2007.
[53] D. Woodruff. Low rank approximation lower bounds in row-update streams. In Advances in
Neural Information Processing Systems, pages 1781–1789, 2014.
[54] D. Woodruff and P. Zhong. Distributed low rank approximation of implicit functions of a
matrix. In manuscript, 2015.
[55] D. P. Woodruff. Sketching as a tool for numerical linear algebra. Foundations and Trends in
Theoretical Computer Science, 10(1-2):1–157, 2014.

14

3 Related Work
Distributed PCA (or distributed SVD) algorithms have been investigated for a long time. One line
of work, developed primarily within the numerical linear algebra literature, studies such algorithms from the perspective of parallelizing existing standard SVD algorithms without sacrificing
accuracy. This approach aims at high accuracy implementations with the least possible communication cost. The distributed models of computation go typically beyond the column-partition
model and arbitrary-partition model that we study in this paper [46]. An extensive survey of this
line of work is out of the scope of our paper; we refer the reader to [35, 50, 33] and references
therein for more details as well as to popular software for distributed SVD such ScaLAPACK [14]
and Elemental [46].
Another line of work for distributed PCA algorithms has emerged within the machine learning
and datamining communities. Such algorithms have been motivated by the need to apply SVD or
PCA to extremely large matrices encoding enormous amounts of data. The algorithms in this line
of research are typically heuristic approaches that work well in practice but come with no rigorous
theoretical analysis. We refer the reader to [47, 42, 9, 15] for more details.
Finally, distributed PCA algorithms in the column-partition model have been recently studied within the theoretical computer science community. Perhaps the more intuitive algorithm
for PCA in this model appeared in [29, 10]: first, a number of left singular vectors and singular
values are computed in each machine; then, the server collects those singular vectors and concatenates them column-wise in a new matrix and then it computes the top k left singular vectors
of this “aggregate” matrix. It is shown in [29, 10] that if the number of singular vectors and singular values in the first step is O(kε−1 ), then, the approximation error in Frobenius norm is at
most (1 + ε) times the optimal Frobenius norm error; the communication cost is O(skmε−1 ) real
numbers because each of the s machines sends O(kε−1 ) singular vectors of dimension m; unfortunately, it is unclear how one can obtain a communication cost in terms of words/bits. A different
algorithm with the same communication cost (only in terms of real numbers since a word/bit
communication bound remained unexplored) is implicit in [41, 31] (see Theorem 3.1 in [31] and
the discussion in Section 2.2 in [41]). Kannan, Vempala and Woodruff proposed arbitrary-partition
model in [37]. They developed
a (1 + ε) Frobenius norm error algorithm with communication cost

−1
2
−4
O skmε + sk ε
words. Bhojanapalli, Jain, and Sanghavi in [13] developed an algorithm
that provides a bound with respect to the spectral norm. Their algorithm is based on sampling
elements from the input matrix, but the communication cost is prohibitive if n is large. The cost
also depends on the condition number of A. Moreover, to implement the algorithm, one needs to
know kA − Ak kF . Finally, [39] discussed a distributed implementation of the “orthogonal iteration” to compute eigenvectors of graphs. The model they consider is different, but perhaps their
protocol could be extended to our model.

4 Preliminaries
Definitions. We formally define the problems and the models here.
Definition 14 (Arbitrary-partition model).
P An m×n matrix A is arbitrarily partitioned into s matrices
Ai ∈ Rm×n , i.e., for i = 1, 2, . . . , s: A = si=1 Ai . There are s machines, and the i-th machine has Ai as
input. There is also another machine, to which we refer to as the “server”, which acts as the central coordinator. The model only allows communication between the machines and the server. The communication
cost of an algorithm in this model is defined as the total number of words transferred between the machines
and the server, where we assume each word is O(log(nms/ε)) bits.
15

Reference

kA − UUT Ak2 ≤

kA − UUT Ak2
F ≤

δ

Communication cost

Implicit in [29]

-

(1 + ε)kA − Ak k2
F

0

-

Theorem 2 in [10]

-

(1 + ε)kA − Ak k2
F

0

-

Theorem 6 in [10]

-

(1 + ε)kA − Ak k2
F

>0

-

Implicit in [41, 31]

-

(1 + ε)kA − Ak k2
F

0

O(mnkε−2 )

Thm 1.1 in [37]*

-

(1 + ε)kA − Ak k2
F

O(1)

Thm 5.1 in [13]

kA − Ak k2 + Γ

-

>0

Remark p. 11 [13]

kA − Ak k2 + Γ

-

>0

Theorem 43*

-

(1 + ε)kA − Ak k2
F

O(1)

- 

O skmε−1 + sk2 ε−4


5 −2
O sm + nk ε
∆


O sm + nk3 ε−2 ∆

Theorem 63

-

(1 + ε)kA − Ak k2
F

O(1)

O(skφε−1 + sk2 ε−4 )

Theorem 74

-

(1 + ε)kA − Ak k2
F

O(1)

O(skφε−1 + sk3 ε−5 )

O (mns
· poly(k, 1/ε))


O(skm + sk3 ε−5 )

Total
 number of arithmetic operations

O mn min{m, n} + mskε−1 min{m, skε−1 }


O mn min{m, n} + mskε−1 min{m, skε−1 }







2 2
3
log sk
log m
O nnz(A) + s m 4k + k m
6
ε
dε
ε

ε

O(poly(m, n, k, s, ε−1 ))


−2
2
O nnz(A) + δ −1 nk5 ε−2 σ1
(A)σk (A)


−2
2
O nnz(A) + δ −1 nk5 ε−2 σ1
(A)σk (A)


O mnkε−2 + mk2 ε−4 + poly(kε−1 ) .

O nnz(A) · log2 ( ns
) + (m + n) · s · poly( k
log( ns
) .
δ
ε
δ

Table 4:

Distributed PCA Algorithms in the column-partition model with s machines (see Definition 15). ∗ indicates that the
algorithm can also be applied in the arbitrary partition model (see Definition 14). A ∈ Rm×n has rank ρ, U ∈ Rm×k with k < ρ is
orthonormal, 0 < ε < 1, and each column of A contains at most φ ≤ m
 non-zero elements. δ is the failure probability. Finally, for
−1 .
notational convenience let Γ := εkA − Ak kF , ∆ := σ12 (A)σk−2 (A) log2 kAk2 kA − Ak k−1
F ε

Definition 15 (Column-partition model). An m × n matrix A is partitioned
arbitrarily
column-wise

P
m×w
i
into s blocks Ai ∈ R
, i.e., for i = 1, 2, . . . , s: A = A1 A2 . . . As . Here, wi = n. There are
s machines, and the i-th machine has Ai as input. There is also another machine, to which we refer to as the
“server”, which acts as the central coordinator. The model only allows communication between the machines
and the server. The communication cost of an algorithm in this model is defined as the total number of words
transferred between the machines and the server, where we assume each word is O(log(nms/ε)) bits.
Definition 16 (The Distributed Principal Component Analysis Problem in arbitrary partition
Psmodel).
m×n
Give an m×n matrix A arbitrarily partitioned into s matrices Ai ∈ R
(i : 1, 2, . . . , s): A = i=1 Ai ,
a rank parameter k < rank(A), and an accuracy parameter 0 < ε < 1, design an algorithm in the model
of Definition 14 which, upon termination, leaves on each machine a matrix U ∈ Rm×k with k orthonormal
columns such that kA − UUT Ak2F ≤ (1 + ε) · kA − Ak k2F , and the communication cost of the algorithm
is as small as possible.
Definition 17 (The Distributed Principal Component Analysis Problem in column partition model).
Given an m × n matrix A partitioned column-wise into s arbitrary blocks Ai ∈ Rm×wi (i : 1, 2, . . . , s):
A = A1 A2 . . . As , a rank parameter k < rank(A), and an accuracy parameter 0 < ε < 1, design an algorithm in the model of Definition 15 which, upon termination, leaves on each machine a matrix
U ∈ Rm×k with k orthonormal columns such that kA − UUT Ak2F ≤ (1 + ε) · kA − Ak k2F , and the
communication cost of the algorithm is as small as possible.
Definition 18 (The Distributed Column Subset Selection Problem). Given an m × n matrix A par
titioned column-wise into s arbitrary blocks Ai ∈ Rm×wi (i : 1, 2, . . . , s): A = A1 A2 . . . As , a
rank parameter k < rank(A), and an accuracy parameter 0 < ε < 1, design an algorithm in the model of
Definition 15 that, upon termination, leaves on each machine a matrix C ∈ Rm×c with c < n columns of
A such that kA − CC† Ak2F ≤ (1 + ε) · kA − Ak k2F , and
1. The number of selected columns c is as small as possible.
2. The communication cost of the algorithm is as small as possible.
Definition 19 (The Distributed Column Subset Selection Problem - rank k subspace version).
Given an m × n matrix Apartitioned column-wise into s arbitrary blocks Ai ∈ Rm×wi (i : 1, 2, . . . , s):
A = A1 A2 . . . As , a rank parameter k < rank(A), and an accuracy parameter 0 < ε < 1,
design an algorithm in the model of Definition 15 that, upon termination, leaves on each machine a matrix C ∈ Rm×c with c < n columns of A and a matrix U ∈ Rm×k with k orthonormal columns with
U ∈ span(C), such that kA − CC† Ak2F ≤ kA − UUT Ak2F ≤ (1 + ε) · kA − Ak k2F , and
16

1. The number of selected columns c is as small as possible.
2. The communication cost of the algorithm is as small as possible.
Definition 20 (Streaming model for Principal Component Analysis). Let all the entries in A ∈ Rm×n
initially be zeros. In the streaming model of computation, there is a stream of update operations that the
q th operation has form (iq , jq , xq ) which indicates that Aiq ,jq should be incremented by xq where iq ∈
{1, ..., m}, jq ∈ {1, ..., n}, xq ∈ R. An algorithm is allowed a single pass over the stream. At the end of
the stream the algorithm stores some information regarding A which we call a “sketch” of A. The space
complexity of an algorithm in this model is defined as the total number of words required to describe the
information the algorithm stores during the stream including the sketch. Each word is O(log(nms/ε))
bits.
Definition 21 (The Streaming Principal Component Analysis Problem). Given an m × n matrix A,
a rank parameter k < rank(A), and an accuracy parameter 0 < ε < 1, design an algorithm that, using
as little space as possible, first finds a sketch of A in the streaming model (see Definition 20) and then,
using only this sketch outputs U ∈ Rm×k with k orthonormal columns such that kA − UUT Ak2F ≤
(1 + ε) · kA − Ak k2F .
Definition 22 (The Streaming Principal Component Analysis Problem (factorization)). Given an
m × n matrix A, a rank parameter k < rank(A), and an accuracy parameter 0 < ε < 1, design an
algorithm that, using as little space as possible, first finds a sketch of A in the streaming model (see Definition 20) and then, using only this sketch outputs A∗k ∈ Rm×n with rank(A∗k ) ≤ k such that kA−A∗k k2F ≤
(1 + ε) · kA − Ak k2F .
Notation. A, B, . . . are matrices; a, b, . . . are column vectors. In is the n × n identity matrix;
0m×n is the m × n matrix of zeros; 1n is the n × 1 vector of ones; ei is the standard basis (whose
dimensionality will be clear from the context): the ith element of ei is one and the rest are zeros.
A(i) and A(j) denotes the ith column and jth row of A, respectively. Aij is the (i, j)th entry in A.
Sampling Matrices. Let A = [A(1) , . . . , A(n) ] ∈ Rm×n and let C = [A(i1 ) , . . . , A(ic ) ] ∈ Rm×c
consist of c < n columns of A. Note that we can write C = AΩ, where the sampling matrix is
Ω = [ei1 , . . . , eic ] ∈ Rn×c (here ei are the standard basis vectors in Rn ). If D ∈ Rc×c is a diagonal
matrix, then AΩD contains c columns of A rescaled with the corresponding elements in D. We
abbreviate S := ΩD, hence the matrix S ∈ Rn×c “samples” and “rescales” c columns from A.
P
Matrix norms. We use the Frobenius and the spectral matrix-norms: kAk2F = i,j A2ij ; kAk2 =
maxkxk2 =1 kAxk2 . kAkξ is used if a result holds for both norms ξ = 2 and ξ = F. The standard
submultiplicativity property for matrix norms implies that for any A and B: kABkξ ≤ kAkξ ·kBkξ .
The triangle inequality for matrix norms implies that kA + Bkξ ≤ kAkξ + kBkξ . A version of the
triangle inequality for the norms squared is: kA + Bk2ξ ≤ 2 · kAk2ξ + 2 · kBk2ξ . A version of the
matrix pythagorean theorem is: if AT B is the all-zeros matrix, then kA + Bk2F = kAk2F + kBk2F .
If V has orthonormal columns, then kAV T kξ = kAkξ , for any A. If P is a symmetric projection
matrix (i.e., P = PT and P2 = P) then, kPAkξ ≤ kAkξ , for any A.

17

Singular Value Decomposition (SVD) and Moore-Penrose Pseudo-inverse.
Rm×n with rank(A) = ρ is



 Σk
VkT
0
A = Uk Uρ−k
,
T
0 Σρ−k
Vρ−k
{z
}
|
|
{z
} | {z }
UA ∈Rm×ρ
ΣA ∈Rρ×ρ

The SVD of A ∈

T
VA
∈Rρ×n

with singular values σ1 ≥ . . . σk ≥ σk+1 ≥ . . . ≥ σρ > 0. Here, k < ρ. We will use σi (A) to denote
the i-th singular value of A when the matrix is not clear from the context. The matrices Uk ∈
Rm×k and Uρ−k ∈ Rm×(ρ−k) contain the left singular vectors of A, and, similarly, the matrices
Vk ∈ Rn×k and Vρ−k ∈ Rn×(ρ−k) contain the right singular vectors of A. It is well-known that
Ak = Uk Σk VkT = AVk VkT = Uk UkT A minimizes kA − Xkξ over P
all matrices X ∈ Rm×n of rank
2
2
2
at most k. Specifically, kA − Ak k2 = σk+1 (A) and kA − Ak kF = ρi=k+1 σi2 (A) (see [32]). A† =
T
n×m denotes the so-called Moore-Penrose pseudo-inverse of A ∈ Rm×n (here Σ−1
VA Σ−1
A
A UA ∈ R
is the inverse of ΣA ). By the SVD of A and A† , for all i = 1, . . . , ρ = rank(A) = rank(A† ):
σi (A† ) = 1/σρ−i+1 (A).

4.1 The best rank k matrix ΠξV,k (A) within a subspace V
Our analysis uses theory involving computing the best rank k approximation of a matrix A within
a given column subspace V. Let A ∈ Rm×n , let k < n be an integer, and let V ∈ Rm×c with
m×n is the best rank k approximation to A in the column span of V.
k < c < n. ΠF
V,k (A) ∈ R
Equivalently, we can write
ΠF
V,k (A) = VXopt ,
where
Xopt =
X∈R

argmin
kA − VXk2F .
:rank(X)≤k

c×n

In order to compute ΠF
V,k (A) given A, V, and k, one can use the following algorithm:
V = YΨ is a qr decomposition of V with Y ∈ Rm×c and Ψ ∈ Rc×c. This step requires O(mc2 )
arithmetic operations.
2: Ξ = Y T A ∈ Rc×n . This step requires O(mnc) arithmetic operations.
1:

T

Ξk = ∆Σ̃k Ṽk ∈ Rc×n is a rank k SVD of Ξ with ∆ ∈ Rc×k , Σ̃k ∈ Rk×k , and Ṽk ∈ Rn×k . This
step requires O(nc2 ) arithmetic operations.
4: Return Y∆∆T Y T A ∈ Rm×n of rank at most k.

3:

Notice that Y∆∆T Y T A ∈ Rm×n is a rank k matrix that lies in the column span of V. The next
lemma is a simple corollary of Lemma 4.3 in [21].
T

Lemma 23. Given A ∈ Rm×n , V ∈ Rm×c and an integer k, Y∆∆T Y T and QŨk Σ̃k Ṽk satisfy:
T

2
kA − Y∆∆T Y T Ak2F ≤ kA − Y∆Σ̃k Ṽk k2F = kA − ΠF
V,k (A)kF

The above algorithm requires O(mnc + nc2 ) arithmetic operations to construct Y, Ψ, and ∆. We will
denote the above procedure as [Y, Ψ, ∆] = BestColumnSubspaceSV D(A, V, k).
Proof. The equality was proven in Lemma 4.3 in [21]. To prove the inequality, notice that for any
matrix X : kA − Y∆(Y∆)† Ak2F ≤ kA − Y∆Xk2F . Also, (Y∆)† = ∆† Y† = ∆T Y T , because both
matrices are orthonormal.
18

Next, we state a version of the above lemma for the transpose of A, equivalently for the best
rank k approximation within the row space of a given subspace R. Let A ∈ Rm×n , let k < n be an
m×n is the best rank k approximation to A
integer, and let R ∈ Rc×n with k < c < m. ΠF
R,k (A) ∈ R
in the row span of R. Equivalently, we can write
ΠF
R,k (A) = Xopt R,
where
Xopt =

argmin
kA − XRk2F .
X∈Rm×c :rank(X)≤k

In order to compute ΠF
R,k (A) given A, R, and k, one can use the following algorithm:
RT = YZ is a qr decomposition of RT with Y ∈ Rn×c and Z ∈ Rc×c . This step requires O(nc2 )
arithmetic operations.
2: Ξ = AY ∈ Rm×c . This step requires O(mnc) arithmetic operations.
3: Ξk = Ũk Σ̃k ∆T ∈ Rm×c is a rank k SVD of Ξ with Ũk ∈ Rm×k , Σ̃k ∈ Rk×k , and ∆ ∈ Rc×k .
This step requires O(mc2 ) arithmetic operations.
4: Return AY∆∆T Y T ∈ Rm×n of rank at most k.
1:

Notice that AY∆∆T Y T ∈ Rm×n is a rank k matrix that lies in the row span of R. The next lemma
is a corollary of Lemma 23 when applied with A := AT and V := RT .
Lemma 24. Given A ∈ Rm×n , R ∈ Rc×n and an integer k, Y∆∆T Y T and Ũk Σ̃k ∆T satisfy:
2
kA − AY∆∆T Y T k2F ≤ kA − Ũk Σ̃k ∆T Rk2F = kA − ΠF
R,k (A)kF .

The above algorithm requires O(mnc + mc2 ) arithmetic operations to construct Y, Z, and ∆. We will
denote the above procedure as [Y, Z, ∆] = BestRowSubspaceSV D(A, R, k).

5 Distributed PCA in the arbitrary partition model
This section describes a fast distributed PCA algorithm with total communication O(msk) words
plus low order terms, which is optimal in the arbitrary partition model in the sense that an
Ω̃(msk) bits lower bound was given by [37]. The algorithm employs, in a novel way, the notion of
projection-cost preserving sketches from [24]. In particular, whereas all previous [48, 23] dimensionreduction-based SVD methods reduce one dimension of the input matrix to compute some approximation to the SVD, our method reduces both dimensions and computes an approximation to
the SVD from a small almost square matrix. Unlike [37] which reduces only one dimension in the
first communication round, we do the reduction on both dimensions in the same round.
We first present the batch version of the algorithm which offers a new low-rank matrix approximation technique; a specific implementation of this algorithm offers a communication-optimal
distributed PCA algorithm (we also discuss in Section 7 a variant of this algorithm that offers
a two-pass space-optimal PCA method in the turnstile streaming model). Before presenting all
these new algorithms in detail, we present the relevant results from the previous literature that
we employ in the analysis.

19

5.1 Projection-cost preserving sketching matrices
In this section, we recap a notion of sketching matrices which we call “projection-cost preserving
sketching matrices”. A sketching matrix from this family is a linear matrix transformation and it
has the property that for all projections it preserves, up to some error, the difference between the
matrix in hand and its projection in Frobenius norm.
Definition 25 (Projection-cost preserving sketching matrices). We say that W ∈ Rn×ξ is an (ε, k)projection-cost preserving sketching matrix of A ∈ Rm×n , if for all rank-k orthogonal projection matrices
P ∈ Rm×m , it satisfies
(1 − ε)kA − PAk2F ≤ kAW − PAWk2F + c ≤ (1 + ε)kA − PAk2F

where c is a non-negative constant which only depends on A and W. We also call AW an (ε, k)-projectioncost preserving sketch of A.
Due to the following lemma, we know that a good rank-k approximation projection matrix of
(ε, k)-projection-cost preserving sketch AW also provides a good rank-k approximation to A.
Lemma 26 (PCA via Projection-Cost Preserving Sketches - Lemma 3 in [24]). Suppose W ∈ Rn×ξ is
∗
an (ε, k)-projection-cost preserving sketching matrix of A ∈ Rm×n . Let P̂ = arg minrank(P)≤k kAW −
∗
PAWk2F . For all P̂, ε′ satisfying rank(P̂) ≤ k, ε′ ≥ 0, if kAW− P̂AWk2F ≤ (1+ε′ )kAW− P̂ AWk2F ,
kA − P̂Ak2F ≤

1+ε
· (1 + ε′ )kA − Ak k2F
1−ε

[24] also provides several ways to construct projection-cost preserving sketching matrices. Because we mainly consider the communication, we just choose one which can reduce the dimension
as much as possible. Furthermore, it is also an oblivious projection-cost preserving sketching matrix.
Lemma 27 (Dense Johnson-Lindenstrauss matrix - part of Theorem 12 in [24]). √
For ε < √
1, suppose
n×ξ
each entry of W ∈ R
is chosen O(log(k))-wise independently and uniformly in {1/ ξ, −1/ ξ} where
ξ = O(kε−2 ) [21]. For any A ∈ Rm×n , with probability at least 0.99, W is an (ε, k)-projection-cost
preserving sketching matrix of A.

5.2 A batch algorithm for the fast low rank approximation of matrices
In this section, we describe a new method for quickly computing a low-rank approximation to
a given matrix. This method does not offer any specific advantages over previous such techniques [48, 23, 18]; however, this new algorithm can be implemented efficiently in the distributed
setting (see Section 5.3) and in the streaming model of computation (see Section 7); in fact we
are able to obtain communication-optimal and space-optimal results, respectively. For completeness as well as ease of presentation, we first present and analyze the simple batch version of the
algorithm. The algorithm uses the dense Johnson-Lindenstrauss matrix of Lemma 27 in order
to reduce both dimensions of A, before computing some sort of SVD to a poly(k/ε) × poly(k/ε)
matrix (see Step 2 in the algorithm below).
Consider the usual inputs: a matrix A ∈ Rm×n , a rank parameter k < rank(A), and an accuracy parameter 0 < ε < 1. The algorithm below returns an orthonormal matrix U ∈ Rm×k such
that
kA − UUT Ak2F ≤ (1 + ε) · kA − Ak k2F .
Algorithm

20

1. Construct two dense Johnson-Lindenstrauss matrices S ∈ Rξ1 ×m , T ∈ Rn×ξ2 with ξ1 = O(kε−2 ), ξ2 =
O(kε−2 ) (see Lemma 27).
2. Construct Ã = SAT.
T
3. Compute the SVD of Ãk = UÃk ΣÃk VÃ
(UÃk ∈ Rξ1 ×k , ΣÃk ∈ Rk×k , VÃk ∈ Rξ2 ×k ).
k

4. Construct X = ATVÃk
5. Compute an orthonormal basis U ∈ Rm×k for span(X) (notice that rank(X) ≤ k).

Theorem 29 later in this section analyzes the approximation error and the running time of the
previous algorithm. First, we prove the accuracy of the algorithm.
Lemma 28. The matrix U ∈ Rm×k with k orthonormal columns satisfies with probability at least 0.98:
kA − UUT Ak2F ≤ (1 + ε) · kA − Ak k2F .
Proof.
T
T
T
k2F = kSAT − SATVÃk VÃ
k2F = kTT AT ST − VÃk VÃ
TT AT ST k2F
kÃ − Ãk k2F = kÃ − ÃVÃk VÃ
k

k

k

T
The first equality follows by the SVD of Ãk = UÃk ΣÃk VÃ
. The second equality is by the
k

construction of Ã = SAT. The third equality is due to ∀M, kMk2F = kM T k2F .
Due to Lemma 27, with probability at least 0.99, ST is an (ε, k)-projection-cost preserving
sketch matrix of TT AT . According to Lemma 26,
T
T
kAT − ATVÃk VÃ
k2F = kTT AT − VÃk VÃ
TT AT k2F ≤
k

k

1+ε
· kAT − (AT)k k2F
1−ε

(2)

Observe that
kUUT AT − ATk2F = kXX† AT − ATk2F
T
≤ kXVÃ
− ATk2F
k

T
= kATVÃk VÃ
− ATk2F
k

≤

1+ε
· kAT − (AT)k k2F
1−ε

The first equality uses the fact that U is an orthonormal basis of span(X). The first inequality is
followed by ∀X, M, N, kXX† M − Mk2F ≤ kXN − Mk2F . The second equality uses the construction
that X = ATVÃk . The second inequality follows by Eqn (2).
Due to Lemma 27, with probability at least 0.99, T is an (ε, k)-projection-cost preserving sketch
matrix of A. Due to Lemma 26,
kUUT A − Ak2F ≤

(1 + ε)2
· kA − Ak k2F
(1 − ε)2

Due to union bound, the probability that ST is an (ε, k)-projection-cost preserving sketch matrix
of TT AT and T is an (ε, k)-projection-cost preserving sketch matrix of A is at least 0.98. Note that
(1+ε)2
is 1 + O(ε) when ε is small enough, so we can adjust ε here by a constant factor to show the
(1−ε)2
statement.
Next, we present the main theorem.
21

Theorem 29. The matrix U ∈ Rm×k with k orthonormal columns satisfies with probability at least 0.98:
kA − UUT Ak2F ≤ (1 + ε) · kA − Ak k2F .

The running time of the algorithm is

O nmkε−2 + mk2 ε−4 + poly(kε−1 ) .

Proof. The correctness is shown by Lemma 28.
Running time.

Next, we analyze the running time of the algorithm:

1. There are a total of (ξ1 × m + ξ2 × n) entries of S and T. It is enough to generate them in
O((n + m)kε−2 ) operations.
2. We first compute AT with O(mnξ2 ) arithmetic operations. Then, we compute SAT with
O(ξ1 mξ2 ) arithmetic operations.
3. This step requires O(poly(k/ε)) operations since we compute the SVD of a O(poly(k/ε)) ×
O(poly(k/ε)) matrix Ã.
4. We already have AT from the second step. Hence, O(mξ2 k) additional arithmetic operations
suffice to compute ATVÃk
5. O(mk2 ) operations suffice to compute an orthonormal basis for X, e.g., with a QR factorization.

5.3 The distributed PCA algorithm
Ps
Recall that the input matrix A ∈ Rm×n is partitioned arbitrarily as: A =
i Ai for i = 1 : s,
m×n
Ai ∈ R
. The idea in the algorithm below is to implement the algorithm in Section 5.2 in the
distributed setting.
Input:
1. A ∈ Rm×n arbitrarily partitioned A =

2. rank parameter k < rank(A)

Ps
i

Ai for i = 1 : s, Ai ∈ Rm×n .

3. accuracy parameter ε > 0
Algorithm
1. Machines agree upon two dense Johnson-Lindenstrauss matrices S ∈ Rξ1 ×m , T ∈ Rn×ξ2 with ξ1 =
O(kε−2 ), ξ2 = O(kε−2 ) (see Lemma 27).
P
2. Each machine locally computes Ãi = SAi T and sends Ãi to the server. Server constructs Ã = i Ãi .
T
3. Server computes the SVD of Ãk = UÃk ΣÃk VÃ
(UÃk ∈ Rξ1 ×k , ΣÃk ∈ Rk×k , VÃk ∈ Rξ2 ×k ).
k

4. Server sends VÃk to all machines.
5. Each machine construct Xi = Ai TVÃk and sends Xi to the server. Server constructs X =
6. Server computes an orthonormal basis U ∈ R

m×k

7. Server sends U to each machine.

for span(X) (notice that rank(X) ≤ k).

P

i

Xi .

Notice that in the first step, S and T can be described using a random seed that is O(log(k))wise independent due to Lemma 27.
22

5.3.1

Main result

The theorem below analyzes the approximation error, the communication complexity, and the
running time of the previous algorithm. Notice that the communication cost of this algorithm
is only given in terms of “real numbers”. The only step where we can not bound the length of
a machine word is when the server communicates VÃk to all machines; and this is because the
entries of VÃk could be unbounded (see the discussion regarding the upper bounds in Section 1.2).
We resolve this issue in the following section.
Theorem 30. The matrix U ∈ Rm×k with k orthonormal columns satisfies w.p. 0.98:
kA − UUT Ak2F ≤ (1 + ε) · kA − Ak k2F .

(3)

The communication cost of the algorithm is
O(msk + s · poly(k/ε))
“real numbers” and the running time is of the order

O nmkε−2 + mk2 ε−4 + poly(kε−1 ) .

Proof. The matrix U - up to the randomness in the algorithm - is exactly the same matrix as in the
batch algorithm in Section 5.2, hence Theorem 29 proves Eqn. 3.
The algorithm communicates O(msk + s · poly(k/ε)) real numbers in total: O(s · poly(k/ε)) in
steps 2 and 4, and O(smk) in steps 5 and 7.
The operations in the algorithm are effectively the same operations as in the batch algorithm
in Section 5.2, hence the analysis of the running time in Theorem 29 shows the claim.

6 Obtaining bit complexity for the distributed PCA algorithm
For the algorithm in the previous section, we were only able to provide a communication upper
bound in terms of “real numbers”. In this section, we describe how to obtain a communication
upper bound in terms of words for the above protocol, where each word is O(log(mnsk/ε)) bits.
The basic idea is that we have a case analysis depending on the rank of the matrix A. If the
rank of A is less than or equal to 2k, we follow one distributed protocol and if the rank is at
least 2k we follow a different protocol. In Section 6.1, Section 6.2, and Section 6.3 we describe
the algorithm that tests the rank of a distributed matrix, and the two PCA protocols, respectively.
Then, in Section 6.4 we give the details of the overall algorithm and in Section 6.5 we give its
analysis.

6.1 Testing the rank of a distributed matrix
Lemma 31. Given A ∈ Rm×n and a rank parameter k < rank(A), there exists a distributed protocol
in the arbitrary partition model to test if the rank of A is less than or equal to 2k using O(sk2 ) words of
communication and succeeding with probability 1 − δ for an arbitrarily small constant δ > 0.
Proof. This is an immediate implementation of a streaming algorithm due to [21] for testing if an
n × n matrix A has rank at least 2k in the streaming model, using O(k2 ) words of space. In that
algorithm, there is a fixed 6nk/δ × n matrix H whose entries are integers of magnitude at most
poly(n), where δ > 0 is an arbitrarily small constant. The algorithm simply chooses 4k random
23

rows from H. Letting H′ be the 2k × n matrix of the first 2k random rows, and H′′ be the n × 2k
matrix whose columns are the next 2k randomly chosen rows, the algorithm just declares that A
has rank at least k iff the rank of H′ AH′′ is 2k.
The above streaming algorithm can be implemented in the distributed setting by having the
coordinator choose 4k random rows of the fixed, known matrix, and send the row identities to
each of the machines. This only takes O(sk) words of communication. Then machine i computes
H′ Ai H′′ , and returns this to the coordinator. The coordinator can then add these up to compute
H′ AH′′ and compute its rank. The total communication is O(sk 2 ) words and the protocol succeeds
with probability at least 1 − δ for an arbitrarily small constant δ > 0. Note that we can assume our
input matrix A, which is m × n, is a square n × n matrix by padding with all-zeros rows. Those
rows of course, will never get communicated in the implementation described above.

6.2 Distributed PCA protocol when rank(A) ≤ 2k
6.2.1

Subsampled Randomized Hadamard Transform and Affine Embeddings

Our algorithms use the following tool, known as “Subsampled Randomized Hadamard Transform” or SRHT for short, to implement efficiently fast dimension reduction in large matrices.
Definition 32 (Normalized Walsh–Hadamard Matrix). Fix an integer m = 2p , for p = 1, 2, 3, .... The
(non-normalized) m × m matrix of the Walsh–Hadamard transform is defined recursively as,




Hm/2 Hm/2
+1 +1
.
,
with
H2 =
Hn =
+1 −1
Hm/2 −Hm/2
1

The m × m normalized matrix of the Walsh–Hadamard transform is equal to H = m− 2 Hm ∈ Rm×m .
Definition 33 (Subsampled Randomized Hadamard Transform (SRHT) matrix). Fix integers ξ and
m = 2p with ξ < m and p = 1, 2, 3, .... An SRHT matrix is an ξ × m matrix of the form
r
m
T=
· RHD;
ξ
• D ∈ Rm×m is a random diagonal matrix whose entries are independent random signs, i.e. random
variables uniformly distributed on {±1}.
• H ∈ Rm×m is a normalized Walsh–Hadamard matrix (see Definition 32).
• R ∈ Rξ×m is a subset or r rows from the n × n identity matrix, where the rows are chosen uniformly
at random and without replacement.
The next lemma argues that an SRHT matrix is a so-called “affine embedding matrix”. The
SRHT is one of the possible choices of Lemma 32 in [23] that will satisfy the lemma.
Lemma 34 (Affine embeddings - Theorem 39 in [23]). Suppose G and H are matrices with m rows,
and G has rank at most r. Suppose T is a ξ × m SRHT matrix (see Definition 33) with ξ = O(r/ε2 ).
Then, with probability 0.99, for all X simultaneously:
(1 − ε) · kGX − Hk2F ≤ kT(GX − H)k2F ≤ (1 + ε) · kGX − Hk2F .
Finally, we note that matrix-vector multiplications with SRHT’s are fast.
Lemma 35 (Fast Matrix-Vector Multiplication, Theorem 2.1 in [6]). Given x ∈ Rm and ξ < m, one
can construct T ∈ Rξ×m and compute Tx in at most 2m log2 (ξ + 1)) operations.
24

6.2.2

Generalized rank-constrained matrix approximations

Let M ∈ Rm×n , N ∈ Rm×c , L ∈ Rr×n , and k ≤ c, r be an integer. Consider the following optimization problem,
Xopt ∈
argmin
kM − NXLk2F .
X∈Rc×r ,rank(X)≤k
Then, the solution Xopt ∈ Rc×r with rank(Xopt ) ≤ k that has the minimum kXopt kF out of all
possible feasible solutions is given via the following formula,


T
T
L† .
MVL VL
Xopt = N† UN UN
k


T

T
T
T
MVL VL
∈
UN UN
MVL VL k ∈ Rm×n of rank at most k denotes the best rank k matrix to UN UN
m×n
R
. This result was proven in [30] (see also [49] for the spectral norm version of the problem).

6.2.3

The PCA protocol

Lemma 36. Suppose the rank ρ of A ∈ Rm×n satisfies ρ ≤ 2k, for some rank parameter k. Then, there is a
protocol for the Distributed Principal Component Analysis Problem in the arbitrary partition model using
O(smk + sk 2 /ε2 ) words of communication and succeeding with probability 1 − δ for an arbitrarily small
constant δ > 0.
Proof. The n×2k matrix H′′ chosen in the protocol of Lemma 31 satisfies that with probability 1−δ,
for an arbitrarily small constant δ > 0, the rank of AH′′ is equal to the rank of A if ρ < 2k. Indeed,
if this were not true, the algorithm could not be correct, as the rank of H′ AH′′ is at most the rank
of AH′′ (and the same algorithm can be used for any ρ < 2k). It follows that with probability 1 − δ,
the column span of AH′′ is equal to the column span of A. Hence, as in the protocol of Lemma
31, the coordinator learns the column space of A, which can be described with O(km) words. The
coordinator thus communicates this to all machines, using O(skm) total words of communication,
assuming the entries of A are integers of magnitude at most poly(mns/ε).
Let C = AH′′ , which is m × 2k. We can set up the optimization problem:
min
kCXCT A − AkF .
rank(X)≤k

(4)

Because the size of C is only m × 2k, every machine can know C by sending a total of O(smk)
words. Since the rank of C and A are small, we can sketch on the left and right using affine
embeddings Tlef t and Tright . Then machine i sends CT Ai Tright to the coordinator, together with
Tlef t Ai Tright . The coordinator computes CT ATright and Tlef t ATright by adding up the sketches,
and sends these back to all the machines. Each machine can then solve the optimization problem
min
kTlef t CXCT ATright − Tlef t ATright kF ,
rank(X)≤k

(5)

each obtaining the same X∗ which is the optimal solution to Eqn (5). Due to Lemma 34, X∗ is a
(1 + O(ε))-approximation to the best solution to Eqn (4). Finally, every machine outputs the same
orthonormal basis U ∈ Rm×k for CX∗ .
We can construct affine embedding matrices Tlef t ∈ Rξ1 ×m and Tright ∈ Rn×ξ2 with ξ1 =
O(k/ε2 ), ξ2 = O(k/ε2 ) (see Definition 33). The total communication of this protocol is O(skm +
sk 2 /ε2 ) words and the success probability can be made 1−δ for an arbitrarily small constant δ > 0.

25

6.3 Distributed PCA protocol when rank(A) > 2k
The idea here is more involved than in the previous subsection and in order to describe the algorithm we need several intermediate results.
6.3.1

Lower bounds on singular values of matrices with integers entries

The first lemma gives a lower bound on the singular values of a matrix with integer entries with
bounded magnitude.
Lemma 37. (Lemma 4.1 of [21], restated) If an m × n matrix A has integer entries bounded in magnitude
by γ, and has rank ρ = rank(A), then the k-th largest singular value A satisfies
σk ≥ (mnγ 2 )−k/(2(ρ−k)) .
Proof. In the proof of Lemma 4.1 of [21], equation (10), it is shown that if λk is the k-th largest
eigenvalue of AT A, then
λk ≥ (mnγ 2 )−k/(ρ−k) .
This implies the k-th singular value σk of A satisfies σk ≥ (mnγ 2 )−k/(2(ρ−k)) .

Next, we state two immediate corollaries of this lemma for future reference.
Corollary 38. If an m × n matrix A has integer entries bounded in magnitude by poly(mns/ε) and
kA − Ak kF > 0, then
kA − Ak kF > (mns/ε)−O(k) .
Proof. Since kA − Ak kF > 0, the rank ρ of A is at least k + 1, so by the preceding lemma,
kA − Ak kF ≥ σk+1 ≥ (poly(mns/ε))−k/2 ,
as desired.
Corollary 39. If an m × n matrix A has integer entries bounded in magnitude by poly(mns/ε) and
rank(A) ≥ 2k, then
kA − Ak kF > 1/poly(mns/ε).
Proof. This follows by plugging ρ = 2k into Lemma 37.
6.3.2

Lower bounds on singular values of integer-perturbed matrices

In this section, we describe a perturbation technique for matrices and provide lower bounds for
the smallest singular value of the perturbed matrix. We start with a theorem of Tao and Vu.
Theorem 40. (Theorem 2.5 of [52]) Let M be an n×n matrix with integer entries bounded in magnitude by
nC for a constant C > 0. Let Nn be a matrix with independent entries each chosen to be 1 with probability
1/2, and −1 with probability 1/2. Then, there is a constant B > 0 depending on C, for which
Pr[k(M + Nn )−1 k2 ≥ nB ] ≤ 1/n.
In words, the result indicates that the spectral norm of the inverse of the perturbed matrix is
bounded from below with high probability. We now describe a simple corollary of this result.

26

Corollary 41. Let M be an n × n matrix with integer entries bounded in magnitude by nC for a constant
C > 0. Let Nn be a matrix with independent entries each chosen to be 1/nD with probability 1/2, and
−1/nD with probability 1/2, where D > 0 is a constant. Then, there is a constant B > 0 depending on C
and D for which
Pr[k(M + Nn )−1 k2 ≥ nB ] ≤ 1/n.
Proof. This follows by Theorem 40 after replacing the constat C in that theorem with C + D, and
scaling by nD .
We need to generalize Corollary 41 to rectangular matrices since we will eventually apply this
perturbation technique to the matrix A to which we would like to compute a distributed PCA.
Lemma 42. Let M be an m × n matrix with integer entries bounded in magnitude by nC for a constant
C > 0, and suppose m ≤ n. Let Nm,n be a matrix with independent entries each chosen to be 1/nD with
probability 1/2 and −1/nD with probability 1/2, where D > 0 is a constant. Then, there is a constant
B > 0 depending on C and D for which
Pr[σm (M + Nm,n ) < 1/nB ] ≤ 1/n,
where σm (M + Nm,n ) denotes the smallest singular value of M + Nm,n .
Proof. Suppose we were to pad M with n − m zero rows, obtaining a square matrix M. Now
consider the n × n matrix Nn,n with independent entries each chosen to be 1/nD with probability
1/2 and −1/nD with probability 1/2. By Corollary 41, all singular values of M + Nn,n are at least
1/nB . Now consider a unit vector x ∈ Rn which is zero on all but its top m coordinates. Let
y ∈ Rm be the unit vector which agrees with x on its top m coordinates. Then,
1/nB ≤ kx(M + Nn,n )k2

= ky(M + Nm,n )k2 ,

where the inequality uses the lower bound on the singular values of M + Nn,n , which occurs with
probability at least 1 − 1/n, and the equality follows by definition of the matrices and vectors we
have defined. As y ∈ Rm can be chosen to be an arbitrary unit vector, it follows that σm (M +
Nm,n ) ≥ 1/nB , which completes the proof.

6.4 Description of algorithm
Using the above results, we are now ready to describe a distributed PCA algorithm whose communication cost can be bounded in terms of machine words and not just in terms of “real numbers”.
As in the algorithm in Section 5.3, we denote with Bi ∈ Rm×n the matrix that B1 arises from A1
after applying the Bernoulli perturbation technique
P discussed above in Lemma 42 and ∀i > 1, Bi
is equal to Ai . Using this notation, we have B := i Bi ∈ Rm×n . Notice that B exactly arises from
A after applying the such Bernoulli perturbation technique.
Input:
1. A ∈ Rm×n arbitrarily partitioned A =
2. rank parameter k < rank(A)

Ps
i

Ai for i = 1 : s, Ai ∈ Rm×n .

3. accuracy parameter ε > 0
Algorithm

27

1. Use the protocol of Lemma 31 with δ = 0.01 to test if the rank of A is less than 2k.
2. If rank(A) ≤ 2k, use the protocol of Lemma 36 to find some orthonormal U ∈ Rm×k .
3. If rank(A) > 2k,
(a) machine 1 locally and independently adds 1/nD with probability 1/2, and −1/nD with probability 1/2, to each of the entries of A1 , where D is the constant of Lemma 42. Note that this
effectively adds the matrix Nm,n of Lemma 42 to the entire matrix A. For notational convenience let B = A + Nm,n , B1 ∈ Rm×n be the local perturbed matrix of machine 1 and ∀i > 1,
Bi is just equal to Ai .
(b) Machines agree upon two dense Johnson-Lindenstrauss matrices S ∈ Rξ1 ×m , T ∈ Rn×ξ2 with
ξ1 = O(kε−2 ), ξ2 = O(kε−2 ) (see Lemma 27).
(c) Each machine locally computes B̃i = SBi T and sends B̃i to the server. Server constructs B̃ =
P
i B̃i .

T
(d) Server computes the SVD of B̃k = UB̃k ΣB̃k VB̃
(UB̃k ∈ Rξ1 ×k , ΣB̃k ∈ Rk×k , VB̃k ∈ Rξ2 ×k ).
k

(e) Now server rounds each of the entries in VB̃k to the nearest integer multiple of 1/nγ for a
sufficiently large constant γ > 0. Let the matrix after the rounding be V̂B̃k
(f) Server sends V̂B̃k to all machines.

(g) Each machine construct X̂i = Bi TV̂B̃k and sends X̂i to the server. Server constructs X̂ =
P
i X̂i .

(h) Server computes an orthonormal basis U ∈ Rm×k for span(X̂) (notice that rank(X̂) ≤ k), e.g.
with a QR factorization.
(i) Server sends U to each machine.

6.5 Main result
The theorem below analyzes the approximation error, the communication complexity, and the
running time of the previous algorithm. Notice that the communication cost of this algorithm is
given in terms of machine words.
Theorem 43. The matrix U ∈ Rm×k with k orthonormal columns satisfies with arbitrarily large constant
probability:
(6)
kA − UUT Ak2F ≤ (1 + ε) · kA − Ak k2F .

The communication cost of the algorithm
is O(msk + s · poly(k/ε)) words and the running time is

O nmkε−2 + mk2 ε−4 + poly(kε−1 ) .

6.6 Proof of Theorem 43

Proof of Eqn. 6 If rank(A) ≤ 2k, then Lemma 36 shows the claim. If rank(A) > 2k, then the
situation is more involved and we prove the approximation bound in the following subsection.
Communication Complexity Step 1 requires O(sk 2 ) words (Lemma 31). Step 2 requires O(smk+
sk 2 /ε2 ) words (Lemma 36). Step 3 requires O(msk + s · poly(k/ε)) words and this follows from the
analysis in Theorem 30. The only difference with Theorem 30 is that now all matrices communicated in the algorithm have real numbers which can be represented efficiently with one machine
word of at most O(log(mns/ε)) bits. To see this, note that each entry of V̂B̃k is bounded in magnitude by nO(B+C)/γ . Therefore, the entries of V̂B̃k and the entries of Xi = Bi TV̂B̃k can each be
described using O(log n) bits, i.e., a constant number of machine words.
28

Running time The operations in the third step in algorithm are effectively the same operations
as in the batch algorithm in Section 5.2 (plus some extra operations whose running time is not
asymptotically larger), hence the analysis of the running time in Theorem 29 shows the claim (the
number of operations in the first two steps of the algorithm is asymptotically less than the number
of operations in the third step).

6.7 Proof of Eqn. 6 if rank(A) > 2k
ξ×p be a matrix of which
Lemma 44 (Result in [21]). Suppose M ∈ Rp×q . For ε > 0, let√P ∈ R√
entries are O(q)-wise independently and uniformly chosen from {−1/ ξ, +1/ ξ} where ξ = O(q/ε2 ).
With arbitrarily large constant probability, P is an ε-subspace embedding matrix of column space of M.
Specifically, ∀x, it has
(1 − ε)kMxk2 ≤ kPMxk2 ≤ (1 + ε)kMxk2

Lemma 45 (Lemma C.2 in [40]). Suppose P is an ε-subspace embedding matrix of column space of
M ∈ Rp×q (q ≤ p). With arbitrarily large constant probability, ∀1 ≤ i ≤ q
(1 − ε)σi (PM) ≤ σi (M) ≤ (1 + ε)σi (PM)
where σi means the ith singular value.
It is suffice to show the following:
Lemma 46. With arbitrarily large constant probability, k(BTVB̃k )† k2 ≤ 8nB
Proof. Notice that BTVB̃k is an m × k matrix. With high probability, S is a ε-subspace embedding
matrix of column space of BTVB̃k due to Lemma 44. According to Lemma 45, σmin (BTVB̃k ) ≥
(1 − ε)σmin (SBTVB̃k ) where σmin means the minimum singular value. Without loss of generality, we assume ε < 1/2. Then σmin (BTVB̃k ) ≥ 12 σmin (SBTVB̃k ). Because VB̃k are the top
k right singular vectors of SBT, σmin (SBTVB̃k ) ≥ σmin (SBT). Applying Lemma 44 again,
TT is a constant error subspace embedding matrix of column space of BT ST . Combining with
Lemma 45, we can make σmin (SBT) ≥ 21 σmin (SB). Notice that Lemma 44 also implies that
q
ξ1 T
ξ1 which means ∀x ∈ Rξ1 , we can
m S is a constant error subspace embedding matrix of R
q
q
make kx ξm1 Sk2 > 12 kxk2 . Thus, ∀x ∈ Rξ1 , kx ξm1 SBk2 > 12 σmin (B)kxk2 . Due to Lemma 42,
q
q
σmin (B) > 1/nB which implies σmin (SB) > 12 · ξm1 σmin (B) > 12 ξm1 n1B > 12 n1B . To conclude
k(BTVB̃k )† k2 ≤ 1/σmin (BTVB̃k ) ≤ 8nB .
Let E = VB̃k − V̂B̃k , we have
kA − UUT AkF = kB − Nm,n − UUT (B − Nm,n )kF
≤ kB − UUT BkF + kNm,n kF
†

= kB − (BTV̂B̃k )(BTV̂B̃k ) BkF + kNm,n kF

≤ kB − (BT(VB̃k + E))(BTVB̃k )† BkF + kNm,n kF

≤ kB − (BTVB̃k )(BTVB̃k )† BkF + kBTE(BTVB̃k )† BkF + kNm,n kF

≤ (1 + O(ε))kB − Bk kF + kBTE(BTVB̃k )† BkF

≤ (1 + O(ε))kB − Bk kF

29

The first equality follows by using the relation B = A + Nm,n . The first inequality uses the triangle inequality for the Frobenius norm and the fact that I − UUT is a projector matrix and can
be dropped without increasing the Frobenius norm. The second equality uses the fact that both
†

UUT and (BTV̂B̃k )(BTV̂B̃k ) are the same projector matrices. The second inequality follows by

∀M, C, X, kM − CC† MkF ≤ kM − CXkF and the relation E = VB̃k − V̂B̃k . The third inequality
uses the triangle inequality. The fourth inequality uses that kB − Bk kF ≥ 1/poly(mns/ε) (follows
from Lemma 42) while kNm,n kF can be made 1/nt for an arbitrarily large integer t, and the fact
that kB − (BTVB̃k )(BTVB̃k )† BkF ≤ (1 + O(ε))kB − Bk kF (implied by replacing all A with B in

Lemma 28). The last inequality uses the fact kBTE(BTVB̃k )† BkF ≤ kBkF kTkF kEkF k(BTVB̃k )† kF kBkF

where kBkF , kTkF , k(BTVB̃k )† kF are poly(nms/ε) (Lemma 46 implies the bound of k(BTVB̃k )† kF )
and we can make |E|F be 1/np for an arbitrarily large p, and kB − Bk kF ≥ 1/poly(mns/ε).
Overall, after rescaling ε, we have
kA − UUT AkF ≤ (1 + ε)kB − Bk kF .

(7)

Finally, we need to relate kB − Bk kF to kA − Ak kF , which we do in the following lemma.
Lemma 47. Let A ∈ Rm×n be any matrix, Nm,n be the random matrix of Lemma 42, and B = A + Nm,n .
Let k < rank(A), ε > 0. Then, for arbitrarily large constant probability,
kB − Bk kF ≤ (1 + ε)kA − Ak kF .
Proof. Since B = A + Nm,n , the proof idea is to relate the singular values of B to the singular
values of A using Weyl’s inequality. Specifically, for i : 1 : m (recall we assume m ≤ n) we have
|σi (B) − σi (A)| ≤ kNm,n k2 ≤ 1/nD .
Rearranging terms in this inequality, and taking squares of the resulting relation we obtain (for
i = 1 : m),
2
σi2 (B) ≤ σi (A) + 1/nD .
(8)

We manipulate the term kB − Bk kF as follows:
v
u X
u m
σi2 (B)
kB − Bk kF = t
i=k+1

v
u X
u m
(σi (A) + 1/nD )2
≤ t
i=k+1

v
u X
u m
= t

σi2 (A) + 1/n2D + 2σi (A)/nD

i=k+1



v
u X
q
q
u m
2
2D
t
σi (A) + 1/n + 2σi (A)/nD
≤
i=k+1

≤ kA − Ak kF + ε · kA − Ak kF + ε · kA − Ak kF

≤ (1 + 2ε)kA − Ak kF .
30

The first inequality uses Eqn. (8). In the third inequality, we used
q
1/n2D ≤ εkA − Ak kF ,

which follows from Corollary 39 for a sufficiently large constant D. Also, we used
q
q
q
2σi (A)/nD ≤ 2σ1 (A)/nD ≤ poly(nms/ε)/nD ≤ εkA − Ak kF ,

where the second inequality follows because kAk2 ≤ poly(nms/ε) and the last inequality uses
again Corollary 39 for a sufficiently large constant D.
Completing the proof. Using Lemma 47 in Eqn (7), taking squares in the resulting inequality,
and rescaling ε concludes the proof.

7 Streaming Principal Component Analysis
In this section, we are interested in computing a PCA of a matrix in turnstile streaming model.
Specifically, there is a stream of update operations that the q th operation has form (iq , jq , xq ) which
indicates that the (iq , jq )th entry of A ∈ Rm×n should be incremented by xq where iq ∈ {1, ..., m}, jq ∈
{1, ..., n}, xq ∈ R. Initially, A is zero. In the streaming setting of computation, we are allowed only
one pass over the update operations, i.e., the algorithm “sees” each update operation one by one
and only once. Upon termination, the algorithm should return a matrix U with k orthonormal
columns which is a “good” basis for span(A) . In Section 7.1.1 below, we describe an algorithm
which gives a space-optimal streaming algorithm. Further more, we provide a variation of this
algorithm in Section 7.1.2 which can output A∗k satisfying kA − A∗k k2F ≤ (1 + ε) · kA − Ak k2F . It
meets the space lower bound shown in [21]. In Section 7.2, we relax the problem and we describe
a two-pass streaming algorithm which is a version of the algorithm of Section 5.2; unfortunately,
the space complexity of this algorithm can be bounded only in terms of “real numbers”. We fix
this in Section 7.3 where we describe an optimal two-pass streaming algorithm.
Inputs to the algorithms in Section 7.1, Section 7.2 and Section 7.3 are a stream of update
operations (i1 , j1 , x1 ), (i2 , j2 , x2 ), ..., (il , jl , xl ), a rank parameter k < rank(A), and an accuracy
parameter 0 < ε < 1.

7.1 One-pass streaming PCA
7.1.1

The algorithm which outputs U

In the following algorithm, the output should be U with orthogonal unit columns and satisfying
kA − UUT Ak2F ≤ (1 + ε) · kA − Ak k2F
Our algorithm uses the following property of random sign matrices:
Lemma 48 (Sketch for regression - Theorem 3.1 in [21]). Suppose both of A and B have m rows and
rank(A) ≤ k. Further more, if each entry of S ∈ Rξ×m is O(k)-wise independently chosen from {−1, +1}
where ξ = O(k/ε) and
X̃ = arg min kSAX − SBk2F
X

, with probability at least 0.99,
kAX̃ − Bk2F ≤ (1 + ε) · min kAX − Bk2F
X

31

Algorithm
1. Construct random sign sketching matrices S ∈ Rξ1 ×m with ξ1 = O(kε−1 ) and R ∈ Rn×ξ2 with
ξ2 = O(kε−1 ) (see Lemma 48)
2. Construct affine embedding matrices Tlef t ∈ Rξ3 ×m and Tright ∈ Rn×ξ4 with ξ3 = O(kε−3 ), ξ4 =
O(k/ε−3 ) (see Definition 33).
3. Initialize all-zeros matrices: M ∈ Rξ3 ×ξ4 , L ∈ Rξ1 ×ξ4 , N ∈ Rξ3 ×ξ2 , D ∈ Rm×ξ2 . We will maintain
M, L, N, D such that M = Tlef t ATright , L = SATright , N = Tlef t AR and D = AR.
4. For (iq , jq , xq ) := (i1 , j1 , x1 ), ..., (il , jl , xl ) (one pass over the stream of update operations)
(a) For all i = 1, ..., ξ3 , j = 1, ..., ξ4 , let Mi,j = Mi,j + (Tlef t )i,iq · xq · (Tright )jq ,j .

(b) For all i = 1, ..., ξ1 , j = 1, ..., ξ4 , let Li,j = Li,j + Si,iq · xq · (Tright )jq ,j .

(c) For all i = 1, ..., ξ3 , j = 1, ..., ξ2 , let Ni,j = Ni,j + (Tlef t )i,iq · xq · Rjq ,j .

(d) For all j = 1, ..., ξ2 , let Diq ,j = Miq ,j + xq · Rjq ,j
5. end

6. Construct X∗ = argmin kN·X·L−Mk2F . (Notice that X∗ = argmin kTlef t (ARXSA − A) Tright k2F .)
rank(X)≤k
rank(X)≤k
T
(UX∗ ∈ Rξ×k , ΣX∗ ∈ Rk×k VX∗ ∈ Rξ×k ); then, compute
7. Compute the SVD of X∗ = UX∗ ΣX∗ VX
∗

T = DUX∗ .
8. Compute an orthonormal basis U ∈ Rm×k for span(T).

Theorem 52 later in this section analyzes the approximation error, the space complexity, and
the running time of the previous algorithm. First, we prove a few intermediate results.
Lemma 49. For all matrices X ∈ Rξ2 ×ξ1 and with probability at least 0.98:
(1 − ε)2 kARXSA − Ak2F ≤ kTlef t (ARXSA − A) Tright k2F ≤ (1 + ε)2 kARXSA − Ak2F
Proof. Notice that rank(AR) ≤ ξ2 = O(kε−1 ). Then from Lemma 34 (G := AR and H := A),
with probability at least 0.99 for all Y ∈ Rξ2 ×n :
(1 − ε)kARY − Ak2F ≤ kTlef t (ARY − A)k2F ≤ (1 + ε)kARY − Ak2F
Replacing Y = XSA, for an arbitrary X ∈ Rξ2 ×ξ1 , we obtain that with probability at least 0.99 and
for all X ∈ Rξ2 ×ξ1 simultaneously:
(1 − ε)kARXSA − Ak2F ≤ kTlef t (ARXSA − A) k2F ≤ (1 + ε)kARXSA − Ak2F
Now notice that rank(SA) ≤ ξ1 = O(kε−1 ). Then from Lemma 34 (G := AT ST and H :=
T
m×ξ2 :
AT Tlef
t ), with probability at least 0.99 for all Z ∈ R
T
2
T
T T T
T T
2
T T T
T T
2
(1 − ε)kA T ST ZT − AT Tlef
t kF ≤ kTright (A S Z − A Tlef t )kF ≤ (1 + ε)kA S Z − A Tlef t kF
T
ξ2 ×ξ1 , we obtain that with probability at least
Replacing Z = XT RT AT Tlef
t for an arbitrary X ∈ R
0.99 and for all X simultaneously:

(1 − ε)2 k(ARXSA − A)T k2F ≤ k(Tlef t (ARXSA − A) Tright )T k2F ≤ (1 + ε)2 k(ARXSA − A)T k2F

32

Lemma 50. Let Xopt = arg minrank(X)≤k kARXSA−Ak2F with Xopt ∈ Rξ2 ×ξ1 . Then, with probability
at least 0.98
(1 + ε)2
min
kTlef t (ARXSA − A) Tright k2F ≤ kARX∗ SA − Ak2F ≤
kARXopt SA − Ak2F
2
(1
−
ε)
rank(X)≤k
Proof. From Lemma 49, we have that for all matrices X ∈ Rξ2 ×ξ1 and with probability at least 0.98:
(1 − ε)2 · kARXSA − Ak2F ≤ kTlef t (ARXSA − A) Tright k2F ≤ (1 + ε)2 · kARXSA − Ak2F
Applying this for X := Xopt , we obtain:
(e)

(1 − ε)2 · kARXopt SA − Ak2F ≤ kTlef t (ARXopt SA − A) Tright k2F ≤ (1 + ε)2 · kARXopt SA − Ak2F
Applying this for X := X∗ , we obtain:
(f )

(1 − ε)2 · kARX∗ SA − Ak2F ≤ kTlef t (ARX∗ SA − A) Tright k2F ≤ (1 + ε)2 · kARXopt SA − Ak2F
Combining (f ), (e) along with the optimality of X∗ , formally using the relation:
kTlef t (ARX∗ SA − A) Tright k2F ≤ kTlef t (ARXopt SA − A) Tright k2F ,
shows the claim.
In words, the lemma indicates that in order to (1 + ε)-approximate the optimization problem
minrank(X)≤k kARXSA − Ak2F , it suffices to “sketch” the matrix ARXSA − A from left and right
with the matrices Tlef t and Tright . Recall that X∗ = argmin kTlef t (ARXSA − A) Tright k2F .
rank(X)≤k
Lemma 51. Let Xopt = arg minrank(X)≤k kARXSA − Ak2F . Then, with probability at least 0.98
kARXopt SA − Ak2F ≤ (1 + ε) · kA − Ak k2F ,
T
T
Proof. Suppose the SVD of A = UA ΣA VA
and Ak = UAk ΣAk VA
. Consider about the following
k
regression problem:
min kUAk X − Ak2F
rank(X)≤k

Since we can choose X to be

T
UA
A,
k

we have

min

rank(X)≤k

kUAk X − Ak2F ≤ kA − Ak k2F

Let X̃ = arg minrank(X)≤k kSUAk X − SAk2F . According to Lemma 48, with probability 0.99
kUAk X̃ − Ak2F ≤ (1 + ε)

min

rank(X)≤k

kUAk X − Ak2F ≤ (1 + ε)kA − Ak k2F

Since X̃ minimizes kSUAk X − SAk2F , the row space of X̃ is in the row space of SA (Otherwise,
we can project SUAk X into the row space of SA to get a better solution). We denote X̃ = WSA
where rank(W) ≤ k. Now, consider about another regression problem:
min

rank(X)≤k

k(WSA)T XT − AT k2F
33

If we choose X to be UAk , we have
min

rank(X)≤k

k(WSA)T XT − AT k2F ≤ kUAk X̃ − Ak2F ≤ (1 + ε)kA − Ak k2F

Let X̂ = arg minrank(X)≤k kRT (WSA)T XT − RT AT k2F . According to Lemma 48, with probability
0.99
T

kRT (WSA)T X̂ − RT AT k2F ≤ (1 + ε)

min

rank(X)≤k

k(WSA)T XT − AT k2F ≤ (1 + ε)2 kA − Ak k2F

Since X̂ minimizes kXWSAR − ARk2F , the column space of X̂ is in the column space of AR. We
denote X̂ = ARG where rank(G) ≤ k. Thus, we have
kARXopt SA − Ak2F ≤ kARGWSA − Ak2F ≤ (1 + ε)2 kA − Ak k2F
We scale the ε with a constant factor. The statement is shown by applying union bound.
Theorem 52. The matrix U ∈ Rm×k with k orthonormal columns satisfies w.p. 0.96:
kA − UUT Ak2F ≤ (1 + O (ε)) · kA − Ak k2F .

The space complexity of the algorithm is O mk/ε + poly(kε−1 ) words, and the running time for each
update operation is O(poly(kε−1 )) and the total running time is of the order O l · poly(kε−1 ) + mk2 ε−1
where l is total number of updates.
Proof. Because UUT projects the columns of A into the column space of Tlef t ARX∗ ,
kA − UUT Ak2F ≤ kTlef t (ARXSA − A) Tright k2F
According to Lemma 50 and Lemma 51,
kA − UUT Ak2F ≤ (1 + O (ε)) · kA − Ak k2F
To see the space complexity of the algorithm, observe that it suffices to maintain the matrices
in the fourth step of the algorithm. Furthermore, observe that by the end of the stream: M =
Tlef t ATright ∈ Rξ3 ×ξ4 , L = SATright ∈ Rξ1 ×ξ4 , N = Tlef t AR ∈ Rξ3 ×ξ2 and D = AR ∈ Rm×ξ2 .
Those matrices form the so called “sketch” of the algorithm.
Since all of ξ1 , ξ2 , ξ3 , ξ4 are O(kε−1 ), the running time for each update operation in the fourth
step is O(poly(kε−1 )). We can do the sixth step by using the result of Section 6.2.2. The computation takes poly(kε−1 ) running time. In the seventh step, computing SVD needs poly(kε−1 ),
and computing T needs O(mk2 ε−1 ). We can use O(mk2 ) to compute U in the final step, e.g. QR
decomposition.
7.1.2

A variation which outputs A∗k

We just slightly modify the previous algorithm in Section 7.1.1 to get the following one which can
output a matrix A∗k ∈ Rm×n satisfying
kA − A∗k k2F ≤ (1 + ε) · kA − Ak k2F
Algorithm

34

1. Construct random sign sketching matrices S ∈ Rξ1 ×m with ξ1 = O(kε−1 ) and R ∈ Rn×ξ2 with
ξ2 = O(kε−1 ) (see Lemma 48)
2. Construct affine embedding matrices Tlef t ∈ Rξ3 ×m and Tright ∈ Rn×ξ4 with ξ3 = O(kε−3 ), ξ4 =
O(k/ε−3 ) (see Definition 33).
3. Initialize all-zeros matrices: M ∈ Rξ3 ×ξ4 , L ∈ Rξ1 ×ξ4 , N ∈ Rξ3 ×ξ2 , D ∈ Rm×ξ2 and C ∈ Rξ1 ×n . We
will maintain M, L, N, D, C such that M = Tlef t ATright , L = SATright , N = Tlef t AR, D = AR
and C = SA.
4. For (iq , jq , xq ) := (i1 , j1 , x1 ), ..., (il , jl , xl ) (one pass over the stream of update operations)
(a) For all i = 1, ..., ξ3 , j = 1, ..., ξ4 , let Mi,j = Mi,j + (Tlef t )i,iq · xq · (Tright )jq ,j .

(b) For all i = 1, ..., ξ1 , j = 1, ..., ξ4 , let Li,j = Li,j + Si,iq · xq · (Tright )jq ,j .

(c) For all i = 1, ..., ξ3 , j = 1, ..., ξ2 , let Ni,j = Ni,j + (Tlef t )i,iq · xq · Rjq ,j .

(d) For all j = 1, ..., ξ2 , let Diq ,j = Miq ,j + xq · Rjq ,j
(e) For all i = 1, ..., ξ1 , let Ci,jq = Ci,jq + Si,iq · xq

5. end
6. Construct X∗ = argmin kN·X·L−Mk2F . (Notice that X∗ = argmin kTlef t (ARXSA − A) Tright k2F .)
rank(X)≤k
rank(X)≤k
T
(UX∗ ∈ Rξ×k , ΣX∗ ∈ Rk×k VX∗ ∈ Rξ×k ); then, compute
7. Compute the SVD of X∗ = UX∗ ΣX∗ VX
∗

T = DUX∗
T
C
K = VX
∗

8. Output A∗k = TΣX∗ K

Theorem 53. With probability at least 0.96:
kA − A∗k k2F ≤ (1 + O (ε)) · kA − Ak k2F .

The space complexity of the algorithm is O (m + n)k/ε + poly(kε−1 ) words, and the running time for
each update operation is O(poly(kε−1 )) and the total running time is of the order

O l · poly(kε−1 ) + (m + n)k2 ε−1 + mnk
where l is total number of updates.

Proof. Notice that A∗k = ARX∗ SA. According to Lemma 50 and Lemma 51, we have
kA − A∗k k2F ≤ (1 + O (ε)) · kA − Ak k2F .
The total space of “sketch” matrices M, L, N, D is the same as the algorithm in Section 7.1.1. The
maintenance of C needs additional O(nkε−1 ) space.
T
The running time is almost the same as before. Computation of K = VX
C needs additional
∗
O(nk2 ε−1 ), and computation of A∗k = TΣX∗ K needs additional O(mnk).

35

7.2 Two-pass streaming PCA with “real numbers” space complexity
A simple modification of the batch algorithm in Section 5.2 leads to a two-pass streaming algorithm with O(mk) + poly(k/ε) “real numbers” space complexity:
Algorithm
1. Construct two dense Johnson-Lindenstrauss matrices S ∈ Rξ1 ×m , T ∈ Rn×ξ2 with ξ1 = O(kε−2 ), ξ2 =
O(kε−2 ) (see Lemma 27).
2. Initialize all-zero matrices: Ã ∈ Rξ1 ×ξ2 , X ∈ Rm×k .
3. For (iq , jq , xq ) := (i1 , j1 , x1 ), ..., (il , jl , xl ) (first pass over the stream)
4. For all i = 1, ..., ξ1 , j = 1, ..., ξ2 , let Ãi,j = Ãi,j + Si,iq · xq · Tjq ,j .
5. end
T
6. Compute the SVD of Ãk = UÃk ΣÃk VÃ
(UÃk ∈ Rξ1 ×k , ΣÃk ∈ Rk×k , VÃk ∈ Rξ2 ×k ).
k

7. For (iq , jq , xq ) := (i1 , j1 , x1 ), ..., (il , jl , xl ) (second pass over the stream)
8. For all i = 1, ..., ξ2 , j = 1, ..., k, let Xiq ,j = Xiq ,j + xq · Tjq ,i · (VÃk )i,j .
9. end
10. Compute an orthonormal basis U ∈ Rm×k for span(X) (notice that rank(X) ≤ k).

The theorem below analyzes the approximation error, the space complexity, and the running
time of the previous algorithm. Notice that the space complexity of this algorithm is only given in
terms of “real numbers” - we resolve this issue in the following section.
Theorem 54. The matrix U ∈ Rm×k with k orthonormal columns satisfies w.p. 0.98:
kA − UUT Ak2F ≤ (1 + O (ε)) · kA − Ak k2F .
The space complexity of the algorithm is

O mk + poly(kε−1 )

“real numbers”, the running time of each update operation is O(poly(kε−1 )), and the total running time is
of the order

O l · poly(kε−1 ) + mk2
where l is the total number of update operations.

Proof. The matrix U - up to the randomness in the algorithms - is exactly the same matrix as in the
algorithm in Section 5.2, hence Theorem 29 shows the claim.
To see the space complexity of the algorithm, observe that it suffices to maintain the matrices in
the fourth and eighth steps of the algorithm. Furthermore, observe that by the end of the stream:
Ã = SAT,
and
X = ATVÃk
Since ξ1 , ξ2 are poly(kε−1 ), the running time for each update operation is poly(kε−1 ). Computing SVD in the sixth step needs poly(kε−1 ). Computing U in the final step needs O(mk 2 ).
36

7.3 Two-pass streaming PCA with bit space complexity
The previous algorithm could suffer from large space complexity in case we need a lot of machine
words to write down the entries of VÃk . To fix this issue we use the same idea as in the case of the
distributed PCA algorithm in Section 6. This leads to a two-pass streaming algorithm for PCA.
Again, as in the distributed case, we need to test if the rank of A is less than 2k, and then we use
one approach if rank(A) < 2k and another approach if rank(A) ≥ 2k. Both of these approaches
can be implemented with two passes. In the overall streaming PCA algorithm that we would like
to design, we can not wait for the algorithm that tests the rank to finish and then start running
one of the two PCA protocols, because this will lead to a three-pass algorithm (one pass to test
the rank and two passes for the actual PCA protocol). To keep the number of passes to two, we
just start running the PCA protocols in parallel with the protocol that tests the rank of the input
matrix. In the end of the first pass over the stream of update operations, we already know which
protocol to follow, and we just do this, disregarding the other protocol.
We already discussed the streaming version of the algorithm that tests the rank of the matrix
in Lemma 31. Below, we describe separately the two streaming PCA protocols.
7.3.1

Streaming PCA protocol when rank(A) ≤ 2k

The idea here is to implement a streaming version of the distributed PCA protocol in Lemma 36.
Algorithm
1. Construct an n × 2k matrix H′′ as in Lemma 36.

2. Construct affine embedding matrices Tlef t ∈ Rξ1 ×m and Tright ∈ Rn×ξ2 with ξ1 = O(k/ε2 ), ξ2 =
O(k/ε2 ) (see Definition 33).

3. Initialize all-zeros matrices: C ∈ Rm×2k , M ∈ Rξ1 ×ξ2 , L ∈ Rξ×ξ2 , N ∈ Rξ1 ×2k .
4. (iq , jq , xq ) := (i1 , j1 , x1 ), ..., (il , jl , xl ) (first pass over the stream)
(a) For all j = 1, ..., 2k, let Ciq ,j = Ciq ,j + xq · H′′jq ,j
5. end
6. For (iq , jq , xq ) := (i1 , j1 , x1 ), ..., (il , jl , xl ) (second pass over the stream)
(a) For all i = 1, ..., ξ1 , j = 1, ..., ξ2 , let Mi,j = Mi,j + (Tlef t )i,iq · xq · (Tright )jq ,j .

(b) For all i = 1, ..., 2k, j = 1, ..., ξ2 , let Li,j = Li,j + CT i,iq · xq · (Tright )jq ,j .
(c) For all i = 1, ..., ξ1 , j = 1, ..., 2k, let Ni,j = Ni,j + (Tlef t )i,iq · xq · H′′jq ,j .

7. end
8. Construct
X∗ = argmin kN · X · L − Mk2F := argmin kTlef t AH′′ X(H′′ )T AT ATright − Tlef t ATright k2F .
rank(X)≤k
rank(X)≤k
9. Compute an orthonormal basis U ∈ Rm×k for span(CX∗ ).

37

7.3.2

Streaming PCA protocol when rank(A) > 2k

The idea here is to implement a streaming version of step 3 of the algorithm in Section 6.4.
Algorithm
1. Construct two dense Johnson-Lindenstrauss matrices S ∈ Rξ1 ×m , T ∈ Rn×ξ2 with ξ1 = O(kε−2 ), ξ2 =
O(kε−2 ) (see Lemma 27).
2. Initialize all-zero matrices: B̃ ∈ Rξ1 ×ξ2 , X̂ ∈ Rm×k .

3. For p := 1, ..., ξ1 , q := 1, ..., ξ2 , i := 1, ..., m, j := 1, ..., n
4. Let ri,j = −1/nD with probability 1/2 and ri,j = 1/nD with probability 1/2.

5. B̃p,q = B̃p,q + Sp,i · ri,j · Tj,q
6. end

7. For (iq , jq , xq ) := (i1 , j1 , x1 ), ..., (il , jl , xl ) (first pass over the stream)
8. For all i = 1, ..., ξ1 , j = 1, ..., ξ2 , let B̃i,j = B̃i,j + Si,iq · xq · Tjq ,j .
9. end

T
10. Compute the SVD of B̃k = UB̃k ΣB̃k VB̃
(UB̃k ∈ Rξ1 ×k , ΣB̃k ∈ Rk×k , VB̃k ∈ Rξ2 ×k ). Then, we
k
round each of the entries of VB̃k to the nearest integer multiple of 1/nγ for a sufficient large constant
γ > 0. Let the matrix after the rounding be V̂B̃k .

11. For (iq , jq , xq ) := (i1 , j1 , x1 ), ..., (il , jl , xl ) (second pass over the stream)
12. For all i = 1, ..., ξ2 , j = 1, ..., k, let X̂iq ,j = X̂iq ,j + xq · Tjq ,i · (V̂B̃k )i,j .

13. end

14. Compute an orthonormal basis U ∈ Rm×k for span(X̂) (notice that rank(X̂) ≤ k).

7.3.3

Main result

The theorem below analyzes the approximation error, the space complexity, and the running time
of the previous algorithm. Notice that the space complexity of this algorithm is given in terms of
machine words.
Theorem 55. The matrix U ∈ Rm×k with k orthonormal columns satisfies with arbitrarily large constant
probability:
kA − UUT Ak2F ≤ (1 + ε) · kA − Ak k2F .

The space complexity of the algorithm is


O mk + poly(kε−1 )

machine words, the running time of each update operation is O(poly(kε−1 )), and the total running time is
of the order

O l · poly(kε−1 ) + mk2
where l is the total number of update operations.

Proof. The matrix U - up to the randomness in the algorithms - is exactly the same matrix as in the
algorithm in Theorem 43, hence the approximation bound in that theorem shows the claim.
To see the space complexity of the algorithm, observe that it suffices to maintain all the “sketch”
matrices in the both two protocols above.
Since ξ1 , ξ2 in both protocols are poly(kε−1 ), the running time for each update operation is
O(poly(kε−1 )). Additional O(mk 2 ) running time is caused by computing U in the final step.
38

8 Distributed PCA for sparse matrices in column-partition model
Recall that in the problem of Definition
15 we are given

P 1) an m × n matrix A partitioned columnwise as A = A1 A2 . . . As , with Ai ∈ Rm×wi ( wi = n); 2) a rank parameter k < rank(A);
3) an accuracy parameter ε > 0. We would like to design an algorithm that finds an m × k matrix
U with UT U = Ik and, upon termination, leaves this matrix U in each machine of the network.
The high level idea of our algorithm is to find, in a distributed way, a small set - O(kε−1 ) - of
columns from A and then find U in the span of those columns. To choose those O(kε−1 ) columns
of A in a distributed way we implement the following three-stage sampling procedure:
1. Local sampling: Each machine samples O(k) columns using the sampling algorithm from [24].
2. Global sampling: The server collects the columns from each machine and down-samples
them to O(k) columns using the deterministic algorithm from [17].
3. Adaptive sampling: the server sends back to each machine those O(k) columns; then, an
extra of O(kε−1 ) columns are selected randomly from the entire matrix A using [26].
We argue that if C̃ contains the columns selected with this three-stage approach, then, w.p. 0.99,
†

(A)k2F ≤ (1 + O(ε)) · kA − Ak k2F .
kA − C̃C̃ Ak2F ≤ kA − ΠF
C̃,k
(A) to be the rank k matrix that approximates A, we are not
Though we could have used ΠF
C̃,k
(A). To address this issue, using
familiar with any communication-efficient computation of ΠF
C̃,k
(A)k2F ;
an idea from [37], we compute U ∈ span(C̃) such that kA−UUT Ak2F ≤ (1+O(ε))kA−ΠF
C̃,k
this U can be calculated with small communication cost and it is sufficient for our purposes.
Before presenting the new algorithm in detail, we discuss results from previous literature that
we employ in the analysis.

8.1 Background material
8.1.1

Column sampling algorithms

Our distributed PCA algorithm in Section 8 samples columns from the input matrix in three stages:
1. Local sampling: O(k) columns are selected locally in each machine.
2. Global sampling: O(k) columns are selected in the server.
3. Adaptive sampling: an extra O(kε−1 ) columns are selected from the entire matrix A.
In the first stage, we use a sampling algorithm mentioned in [24]. Actually, the sampling
algorithm is the same as the Batson, Spielman, and Strivastava (BSS) spectral sparsification algorithm [11]. But Cohen et al. showed that a small set of columns sampled by BSS sampling
algorithm can provide a projection-cost preserving sketch. In the second stage, we use a deterministic algorithm developed in [17], which itself extends the Batson, Spielman, and Strivastava
(BSS) spectral sparsification algorithm [11]. For the actual algorithm in the first stage we defer
the reader to Theorem 15 in [24]. Lemma 56 states the result. And for the actual algorithm in the
second stage we defer the reader to Lemma 3.6 in [17]. Lemma 57 and Lemma 58 below present
the relevant results. In the third sampling stage, we use an adaptive sampling algorithm from [26].

39

Lemma 56 (BSS Sampling for Projection-Cost Preserving Sketch. Theorem 15 in [24]). Let A ∈
Rm×w be an arbitrary matrix. For any integer 0 < k < m, real numbers 0 < ε < 1, 0 < δ, there
is a randomized algorithm that runs in O(nnz(A) log(m/δ) + m · poly(k, ε, log(1/δ))) time, and constructs a “sampling/rescaling” w × O(k/ε2 ) matrix S with probability at least 1 − δ such that S is
an (ε, k)-projection-cost preserving sketching matrix of A. We denote this sampling procedure as S =
BssSamplingI(A, k, ε, δ).
In words, there exists a randomized algorithm that runs in near input-sparsity running time
can construct a projection-cost preserving sketch of a given matrix A by sampling and rescaling a
small subset of columns of A.
Lemma 57 (Dual Set Spectral-Frobenius Sparsification. Lemma 3.6 in [17]). Let V ∈ Rw×k be a
matrix with w > k and V T V = Ik . Let E ∈ Rm×w be an arbitrary matrix. Then, given an integer ℓ such
that k < ℓ ≤ w, there exists a deterministic algorithm that runs in O ℓwk2 + mw time, and constructs a
“sampling/rescaling” w × ℓ matrix S such that

 
p 2
σk2 V T S ≥ 1 − k/ℓ ,

kESk2F ≤ kEk2F .

Specifically, rank(V T S) = k. We denote this sampling procedure as S = BssSamplingII(V, E, ℓ).
In words, given V and E, there exists an algorithm to select (and rescale) a small number of
columns from E and rows from V such that: 1) the Frobenius norm squared of the sampled E is
less than the Frobenius norm squared of E; 2) the sampled V has rank equal to the rank of V; and

p 2
3) the smallest singular value squared of the sampled V is bounded from below by 1 − k/ℓ .

Lemma 58 (Constant factor column-based matrix reconstruction; Theorem 5 in [17]). Given matrix G ∈ Rm×α of rank ρ and a target rank k 1 , there exists a deterministic algorithm that runs in
O αm min{α, m} + αck2 time and selects c > k columns of G to form a matrix C ∈ Rm×c with
†

kG − CC

Gk2F

≤



 rank

X(G)
p −2
1 + 1 − k/c
·
σi2 (G).
i=k+1

The algorithm in this theorem finds C as C = GS, where S = BssSamplingII(V, G − GVV T , c)
and V ∈ Rα×k contains the top k right singular vectors of G. We denote this sampling procedure as
C = DeterministicCssF robenius(G, k, c).
In words, there exists a deterministic algorithm, running in polynomial time, to select any
number of c > k columns from the given matrix G, such that the residual error, in Frobenius
norm, from projecting G onto the span of the sampled columns is bounded from above with
respect to the residual error of the best rank k matrix from the SVD of G. Notice that
rank
X(G)
i=k+1

σi2 (G) = kG − Gk k2F ,

where Gk ∈ Rm×α has rank at most k and is computed via the SVD of G.
1

The original Theorem 5 in [17] has the assumption that k < ρ, but this assumption can be dropped having the result
unchanged. The only reason the assumption k < ρ exists is because otherwise column subset selection is trivial.

40

Before presenting the adaptive sampling theorem from [26] we introduce some useful notation
m×n is the
used in the lemma. Let A ∈ Rm×n , k < n, and V ∈ Rm×c with k < c < n. ΠF
V,k (A) ∈ R
best rank k approximation to A - wrt the Frobenius norm - in the column span of V. Equivalently,
ΠF
V,k (A) = VXopt , where
Xopt =

argmin
kA − VXk2F .
X∈Rc×n ,rank(X)≤k

Lemma 59 (Adaptive sampling; Theorem 2.1 of [26]). Given A ∈ Rm×n and V ∈ Rm×c1 (with
c1 ≤ n, m), define the residual matrix Ψ = A − VV† A ∈ Rm×n . For j = 1, . . . , n, let pj be a probability
distribution such that pj ≥ βkΨ(j) k22 /kΨk2F , for some 1 > β > 0, where Ψ(j) is the j-th column of the
matrix Ψ. Sample c2 columns from A in c2 i.i.d. trials, where in each trial the j-th column is chosen with
probability pj . Let C2 ∈ Rm×c2 contain the c2 sampled columns and let C = [V C2 ] ∈ Rm×(c1 +c2 )
contain the columns of V and C2 . Then, for any integer k > 0,
rank
i
h
X(A)


2
F
†
2
E kA − CC AkF ≤ E kA − ΠC,k (A)kF ≤
σi2 (A) +
i=k+1

k
· kA − VV† Ak2F .
β · c2

Given A and C, this method requires O(c1 mn) time to compute Ψ, another O(mn) time to compute
the probabilities pj ’s and another O(n + c2 ) time for the sampling step - using the method in [51]. In
total, the method requires O(c1 mn + c2 ) time to compute C2 . We denote this sampling procedure as
C2 = AdaptiveCols(A, V, c2 , β).
In words, given the matrix A and the subspace V, there exists a randomized algorithm to
sample an additional of c2 columns from A, based on probabilities from the residual error matrix
A − VV† A, such that residual error, in Frobenius norm, from projecting A onto the span of the
union of the columns of V and the sampled columns is bounded from above with respect to the
best rank k approximation to A, the residual A − VV† A, and the number of sampled columns c2 .
8.1.2

Distributed adaptive sampling

In our distributed PCA algorithm below, we also need to use a distributed version of the previous
adaptive sampling procedure. We provide some preliminary results for that task in this section.
Lemma 60. Suppose that the columns of an m × n matrix A are partitioned arbitrarily across the machines
into matrices A1 , . . . , As . Let C be an arbitrary m × r matrix. Consider the distribution p on n columns
in which
kaj − CC† aj k2F
,
pj =
kA − CC† Ak2F

where aj is the j-th column of A (j = 1 : n here).
For each i ∈ [s], let some value βi satisfies

kAi − CC† Ai k2F ≤ β i ≤ 2kAi − CC† Ai k2F .
For each j ∈ [n], if column aj is held on the i-th server (denoted aij ), then let
q j = Ps

βi

i′ =1 βi′

Then for each j ∈ [n],

·

kaij − CC† aij k22

kAi − CC† Ai k2F

pj /2 ≤ qj ≤ 2pj .
41

.

Proof. By definition of the β i , we have that

Hence,

kAi − CC† Ai k2F
2kAi − CC† Ai k2F
βi
P
≤
.
≤
s
2kA − CC† Ak2F
kA − CC† Ak2F
i′ =1 βi′
pj /2 ≤ qj ≤ 2pj .

Lemma 61. Suppose the coordinator has an m × r matrix C of columns of an m × n matrix A, where
r = O(k). Suppose the entries of A are integers bounded by poly(mns/ε) in magnitude, and let the
columns of A be partitioned arbitrarily across the servers into matrices A1 , . . . , As .
There is a protocol using O(skm) machine words of communication for the coordinator to learn values
i
β so that for all i ∈ [s],
kAi − CC† Ai k2F ≤ β i ≤ 2kAi − CC† Ai k2F .
Proof. The coordinator first sends C to all machines. The i-th server locally computes kAi −
CC† Ai k2F . If this number is 0, he/she sends 0 to the coordinator, and in this case β i satisfies
the requirement in the statement of the lemma.
Otherwise, consider the matrix Bi which is formed by concatenating the columns of Ai with
those of C. Then
kBi − CC† Bi k2F = kAi − CC† Ai k2F > 0,

(9)

since the columns of Bi in C contribute a cost of 0. However, since Bi contains the columns of C
and its cost is non-zero, it implies the rank of Bi is at least r + 1. By Corollary 38, this implies
kBi − CC† Bi k2F > (mns/ε)−O(k) ,
which by (9) gives the same lower bound on kAi − CC† Ai k2F . Note also that
kAi − CC† Ai k2F ≤ poly(mns/ε),
since kAk2F ≤ poly(mns/ε). This implies if the i-th machine sends β i to the coordinator, where β i
is the nearest power of 2 to kAi − CC† Ai k2F , there will only be O(k log(mns/ε)) possible values
of β i , and hence each can be specified using O(log k + log log(mns/ε)) bits, i.e., a single machine
word. This completes the proof.
8.1.3

Low-rank matrix approximations within a subspace

The final stage of our distributed PCA algorithm below finds U ∈ span(C̃) such that the error
of the residual matrix A − UUT A is “small”. To implement this step, we employ an algorithm
developed in [37].
Lemma 62. Let A ∈ Rm×n be the input matrix and V ∈ Rm×c be the input subspace. We further assume
that for some rank parameter k < c and accuracy parameter 0 < ε < 1 :
2
2
kA − ΠF
V,k (A)kF ≤ (1 + O(ǫ))kA − Ak kF .

Let V = YΨ be a qr decomposition of V with Y ∈ Rm×c and Ψ ∈ Rc×c . Let Ξ = Y T AW T ∈ Rc×ξ ,
√
√
where W T ∈ Rn×ξ with ξ = O(c/ε2 ), each element of which is chosen i.i.d. to be {+1/ n, −1/ n} with
42

probability 1/2. Let ∆ ∈ Rc×k contain the top k left singular vectors of Ξ. Then, with probability at least
1 − e−c ,
kA − Y∆∆T Y T Ak2F ≤ (1 + ε)kA − Ak k2F .
Y, and ∆ can be computed in O(mnξ) time. We denote this procedure as
[Y, ∆] = ApproxSubspaceSV D(A, V, k, ε).
Proof. This result was essentially proven inside the proof of Theorem 1.1 in [36]. Specifically, the
error bound proven in [36] is for the transpose of A (also Y, ∆ are denoted with U, V in [36]).
As for the running time, given A, V, and k, one can compute (i) Y in O(mc2 ) time; (ii) Ξ in
O(mnξ + mcξ) time; and (iii) ∆ in O(c2 ξ) time.
In words, the lemma indicates that given the matrix A and the subspace V such that
2
2
kA − ΠF
V,k (A)kF ≤ (1 + O(ǫ))kA − Ak kF ,

for some small ε > 0, there exists a randomized algorithm to compute matrices Y ∈ span(V)
and ∆ such that the residual error, in Frobenius norm, from projecting A onto the span of Y∆ is
also bounded by (1 + O(ǫ))kA − Ak k2F . Notice that ∆ contains the top k left singular vectors of
a “sketched” version of Y T A, a property which makes the computation particularly effective in
terms of running time.

8.2 Detailed description of the algorithm
Input:
1. A ∈ Rm×n partitioned column-wise A = A1

A2

2. rank parameter k < rank(A)


P
. . . As ; for i = 1 : s, Ai ∈ Rm×wi ; i wi = n.

3. accuracy parameter ε > 0
Algorithm
1. Local Column Sampling
(a) For each sub-matrix Ai ∈ Rm×wi , compute Ci ∈ Rm×ℓ containing ℓ = O(k) columns from
Ai as follows: Ci = Ai Si . Here, Si has dimensions wi × ℓ and is constructed as follows:
Si = BssSamplingI(Ai , 4k, 1/2, 0.01/s) (see Lemma 56).
(b) Machine i sends Ci to the server.
2. Global Column Sampling
(a) Server constructs m × (s · ℓ) matrix G containing (s · ℓ) actual columns from A as follows:
G = C1 C2 . . . Cs . Then, server constructs C ∈ Rm×c1 via choosing c1 = 4k columns
from G as follows: C = DeterministicCssF robenius(G, k, c1 ) (see Lemma 58).

(b) Server sends C to all the machines.
3. Adaptive Column Sampling

(a) Machine i computes Ψi = Ai − CC† Ai ∈ Rm×wi and then computes βi as it was described in
Lemma 61. Machine i sends βi to server.
(b) Server computes probability distribution gi = Pβiβi . Server samples i.i.d. with replacement
i
⌈ 50k/ε ⌉ samples (machines) from gi . Then, server determines numbers ti (i = 1, 2, . . . , s),
where ti is the number of times the ith machine was sampled. It sends the ti ’s to the machines.

43

(c) Machine i computes probabilities qji = kxk22 /kΨi k2F (j = 1 : wi ), where x is the jth column of
Ψi . And now machine i samples ti samples fromP
it’s local probability distribution and sends
the corresponding columns to the server. Let c2 = i ti = ⌈ 50k/ε ⌉.

(d) Server collects the columns and assigns
them to Ĉ ∈ Rm×c2 . Server constructs C̃ to be the

m × (c1 + c2 ) matrix: C̃ = C; Ĉ . Let c = c1 + c2 = 4k + ⌈ 50k/ε ⌉.

4. Rank-k matrix in the span of C̃

(a) Server sends C̃ to all the machines and each machine computes (the same) qr factorization of
C̃: C̃ = YR where Y ∈ Rm×c has orthonormal columns and R ∈ Rc×c is upper triangular.

1/2. Implicitly
(b) Machine i generates W̃i ∈ Rξ×wi with ξ = O(c/ε2 ) to be i.i.d.
 {+1, −1} w.p.ξ×n
, with W̃ ∈ R
. Machine i
all machines together generate W̃ = W̃1 W̃2 . . . W̃s
T

computes Hi = C̃ (Ai WiT ) ∈ Rc×ξ . Machine i sends Hi to the server.
P
(c) Server computes Ξ = si=1 Hi ∈ Rc×ξ and sends this back to all the machines. Now machines

√
compute Ξ := R−1 · Ξ · 1/ n(:= Y T AWT ), where
∈ Rξ×n is
W √= W1 W2 . . . Ws
√
random matrix each element of which is {+1/ n, −1/ n} w.p. 1/2, and then they compute
∆ ∈ Rc×k to be the top k left singular vectors of Ξ. Each machine computes U = Y · ∆ ∈ Rm×k .

Discussion. A few remarks are necessary for the last two stages of the algorithm. The third stage
(adaptive column sampling), implements the adaptive sampling method of Lemma 59. To see this,
note that each column in A (specifically the jth column in Ai ) is sampled with probability
gi · qji ≥

1 ||xij ||22
,
·
2 kΨk2F

where Ψ = A − CC† A. This follows from Lemma 60. Overall, this stage effectively constructs Ĉ
such that (see Lemma 59)
Ĉ = AdaptiveCols(A, C, c2 , 1/2).
The last stage in the algorithm implements the algorithm in Lemma 62. To see this, note that
W satisfies the properties in the lemma. Hence,
[Y, ∆] = ApproxSubspaceSV D(A, C̃, k, ε).

8.3 Main Result
The theorem below analyzes the approximation error, the communication complexity, and the
running time of the previous algorithm.
Theorem 63. The matrix C̃ ∈ Rm×c with c = O(k/ε) columns of A satisfies w.p. at least 0.98,
†

(A)k2F ≤ (1 + O(ε)) · kA − Uk UkT Ak2F .
kA − C̃C̃ Ak2F ≤ kA − ΠF
C̃,k

(10)

The matrix U ∈ Rm×k with k orthonormal columns satisfies with probability at least 0.97,
kA − UUT Ak2F ≤ (1 + O(ε)) · kA − Uk UkT Ak2F .

(11)

Let each column of Ahas at most φ non-zero elements. Then, the communication cost of the algorithm is
O skφε−1 + sk2 ε−4 and the running time is O (mns · poly(k, 1/ε)) .
44

8.4 Proof of Theorem 63
8.4.1

Proof of Eqn. 10

From Lemma 59:
h

†

E kA − C̃C̃

Ak2F

i

h

≤ E kA −

(A)k2F
ΠF
C̃,k

i

≤

rank
X(A)

σi2 (A) +

i=k+1

2ε
· kA − CC† Ak2F .
50

(12)

Therefore, we need to bound kA − CC† Ak2F now.
Lemma 64. With probability at least 0.99, matrix G is a (4k, 1/2)-projection-cost preserving sketch of A.
Proof. For each i, Ci = Ai Si is a (4k, 1/2)-projection-cost preserving sketch of Ai with probability
at least 1 − 0.01/s. Due to union bound, with probability at least 1 − 0.01, ∀i, Ci is a (4k, 1/2)projection-cost preserving sketch of Ai . Since Ci is a (4k, 1/2)-projection-cost preserving sketch of
Ai , we can assume ci is a constant which is independent from projection matrix P which has rank
at most 4k such that
3
1
kAi − PAi k2F ≤ kCi − PCi k2F + ci ≤ kAi − PAi k2F
2
2
Then, for any projection matrix P which has rank at most 4k we have:
s

1
1X
kAi − PAi k2F
kA − PAk2F =
2
2
i=1

≤

s
X
i=1

kCi − PCi k2F + ci

s

≤

3X
kAi − PAi k2F
2



i=1

3
≤ kA − PAk2F
2

Consider the right hand side of the first inequality, it is equal to
kG −
Let constant c =

Ps

i=1 ci ,

PGk2F

+

s
X

ci

i=1

we have:

1
3
kA − PAk2F ≤ kG − PGk2F + c ≤ kA − PAk2F
2
2

Therefore, G is a (4k, 1/2)-projection-cost preserving sketch of A.
For convenience, we use PG to denote a rank-k projection matrix which provides the best
rank-k approximation Gk = PG G to matrix G. We also use PA to denote a rank-k projection

45

matrix which provides the best rank-k approximation Ak = PA A to matrix A. Thus,
1
kA − CC† Ak2F ≤ kG − CC† Gk2F + c
2
≤ 5kG − PG Gk2F + c
≤ 5kG − PA Gk2F + c


≤ 5 kG − PA Gk2F + c
15
kA − PA Ak2F
≤
2
The first inequality is due to Lemma 64. G is a (4k, 1/2)-projection cost preserving sketch of A
and C has rank at most 4k. The second inequality is based on Lemma 58. The third inequality is
because PG provides the best rank-k approximation to G. The fourth inequality used the fact c is
non-negative. The fifth inequality held because G is a (4k, 1/2)-projection cost preserving sketch
of A. Thus,
kA − CC† Ak2F ≤ 15kA − Ak k2F
Combining with equation 12, we have:
h

E kA − C̃C̃

†

Ak2F

i

h

≤ E kA −

(A)k2F
ΠF
C̃,k

i

30
≤ (1+ ε)
50

rank
X(A)
i=k+1

σi2 (A) = (1+O(ε))kA−Ak k2F (13)

Concluding the proof. Consider about
rank
i rank
i
h
h
X(A)
X(A)
†
2
2
≤
(A)
+
O(ε)
σ
σi2 (A).
(A)k
E kA − C̃C̃ Ak2F ≤ E kA − ΠF
i
F
C̃,k
i=k+1

(14)

i=k+1

Prank(A)
The expectation is taken w.r.t. the randomness in constructing C̃; hence, the term i=k+1 σi2 (A)
is a constant w.r.t. the expectation operation, which means that Eqn. 14 implies the bound:


rank
rank
X(A)
X(A)
F
2
2


E kA − ΠC̃,k (A)kF −
σi (A) ≤ O(ε) ·
σi2 (A).
i=k+1

(A)k2F −
Let Y := kA − ΠF
C̃,k

Prank(A)
i=k+1

i=k+1

σi2 (A). Y is a random variable with Y ≥ 0, because

(A)) ≤ k. Markov’s inequality on Y implies that with arbitrary constant probability,
rank(ΠF
C̃,k
kA −

ΠF
(A)k2F
C̃,k

kA −

(A)k2F
ΠF
C̃,k

−

rank
X(A)

σi2 (A)

≤

rank
X(A)

σi2 (A) +

i=k+1

≤ O(ε) ·

rank
X(A)

σi2 (A),

rank
X(A)

σi2 (A).

i=k+1

equivalently,

†

i=k+1

O(ε) ·

(A)k2F concludes the proof.
Using kA − C̃C̃ Ak2F ≤ kA − ΠF
C̃,k
46

i=k+1

8.4.2

Proof of Eqn. 11

We would like to apply Lemma 62 for the matrix U ∈ Rm×k in the algorithm. Note that we
have already proved that the matrix C̃ in the algorithm satisfies with probability at least 0.99:
(A)k2F ≤ (1 + O(ε))kA − Ak k2F . Lemma 62 and a simple union bound conclude the
kA − ΠF
C̃,k
proof.
8.4.3

Running time

Next, we give an analysis of the arithmetic operations required in various steps in the algorithm.
1. Local Column Sampling: O(nnz(A) log(ms) + m · poly(k, ε, log(s))) arithmetic operations in
total.
(a) O(nnz(Ai ) log(ms) + m · poly(k, ε, log(s))) for each Si

(b) -

2. Global Column Sampling: O(sk 3 + skm) arithmetic operations in total.
(a) O(sk 3 + skm) to construct C.
(b) 3. Adaptive Sampling O(sk 2 m + mnk + k/ε) arithmetic operations in total.
(a) O(k2 m) for each C† locally; O(wi mk) for each Ψi and O(mwi ) for each ψi .
(b) O(s) in total.
(c) O(k/ε) in total using the method of [51] locally.
(d) O(k/ε).
4. Rank-k matrix in span(C̃): O(mnk/ε3 + smk 2 /ε4 + sk 3 /ε5 ) arithmetic operations in total.
(a) O(smk 2 /ε2 ) in total
T

(b) O(mwi k/ε3 ) for each Ai WiT and another O(mk2 ε−4 ) for each C̃ (Ai WiT )
(c) O(sk 2 ε−3 ) to find Ξ in the server; then another O(sk 3 ε−3 + sk 3 ε−5 )) to update Ξ locally in each machine and another O(sk 3 ε−5 ) to find ∆ locally in each machine; and
O(smk 2 /ε) to find U locally in each machine.
8.4.4

Communication Complexity

Assume that we can represent each element in the input matrix A with at most b bits. Also, we
assume that one word has length b bits and b = O(log(mns/ε)).
1. Local Column Sampling: O(skφ) words.
(a) (b) O(skφ) elements of A.
2. Global Column Sampling:O(skφ) words.
(a) 47

(b) O(skφ) elements of A.
3. Adaptive Sampling: O(s + φk/ε) words.
(a) s integers each of which is representable with O(log k+log log(mns/ε)) (from Lemma 61).
(b) s integers (the ti ’s) each with magnitude at most n, hence representable with b bits.
(c) O(φk/ε) elements of A.
(d) 4. Best rank-k matrix in the span of C̃: O(sφk/ε + sk 2 ε−4 ) words.
(a) O(skφε−1 ) elements of A.
(b) O(sk 2 ε−4 ) numbers each of which can be represented with b bits.
(c) O(sk 2 ε−4 ) numbers each of which can be represented with b bits.
In total the communication complexity is O(skφ/ε + sk 2 /ε4 ) words.

9 Faster Distributed PCA for sparse matrices
We now explain how to modify certain steps of the distributed PCA algorithm in Section 8 in order
to improve the total running time spent to compute the matrix U; this, at the cost of increasing the
communication cost slightly. Specifically, we replace the parts 2-(a), 3-(a), and 4-(b) with similar
procedures, which are almost as accurate as the original procedures but run in time proportional
to the number of the non-zero elements of the underlying matrices. Before presenting the new
algorithm in detail, we discuss results from previous literature that we employ in the analysis. We
remark that in the analysis below we have not attempted to optimize the constants.

9.1 Background material
9.1.1

Constant probability sparse subspace embeddings and sparse SVD

First, we present the so-called sparse subspace embedding matrices of Clarskon and Woodruff [23].
Definition 65. [Sparse Subspace Embedding [23]] We call W ∈ Rξ×n a sparse subspace embedding of
dimension ξ < n if it is constructed as follows, W = Ψ · Y, where
• h : [n] → [ξ] is a random map so that for each i ∈ [n], h(i) = ξ ′ , for ξ ′ ∈ [ξ] w.p. 1/ξ.

• Ψ ∈ Rξ×n is a binary matrix with Ψh(i),i = 1, and all remaining entries 0.
• Y ∈ Rn×n is a random diagonal matrix, with each diagonal entry independently chosen to be +1 or
−1, with equal probability.
Such sparse subspace embedding matrices have favorable properties which we summarize
below:
Lemma 66 ([23]). Let AT ∈ Rn×m have rank ρ and let W ∈ Rξ×n be a randomly chosen sparse subspace
embedding with dimension ξ = Ω(ρ2 ε−2 ), for some 0 < ε < 1. Then, 1) computing AW T requires
O(nnz(A)) time; 2) nnz(AW T ) ≤ nnz(A); and 3) with probability at least 0.99, and for all vectors
y ∈ Rm simultaneously,
(1 − ε)kA T yk22 ≤ kWA T yk22 ≤ (1 + ε)kA T yk22 .
48

Notice that the third property in the lemma fails with constant probability. Based on the above
sparse subspace embeddings, Clarkson and Woodruff [23] described a low-rank matrix approximation algorithm for A that runs in time nnz(A) plus low-order terms.
Lemma 67 (Theorem 47 in [23]; for the exact statement see Lemma 3.4 in [19]). Given A ∈ Rm×n
of rank ρ, a target rank 1 ≤ k < ρ, and 0 < ǫ ≤ 1, there exists a randomized algorithm that computes
Z ∈ Rn×k with ZT Z = Ik and with probability at least 0.99,
kA − AZZT k2F ≤ (1 + ǫ) kA − Ak k2F .

The proposed algorithm requires O nnz(A) + (m + n) · poly(k, ε−1 ) arithmetic operations. We denote
this procedure as Z = SparseSV D(A, k, ε).
Notice that the algorithm in the lemma fails with constant probability.
9.1.2

High probability sparse SVD

Lemma 67 above presents a randomized algorithm to compute quickly an SVD of a sparse matrix;
however, this lemma fails with constant probability. Since in our fast distributed PCA algorithm
we will use the lemma multiple times (specifically, s times, since we need such a fast SVD locally
in each machine), the overall failure probability would be of the order O(s) (by a union bound).
To boost this failure probability, we develop a simple scheme to obtain a high probability result.
Before we proceed, we quickly review the Johnson Lindestrauss transform, which we need in
the analysis.
Lemma 68. [Theorem 1 in [5] for fixed ε = 1/2] Let B ∈ Rm×n . Given β > 0, let r = (1/2)4+2β
2 −(1/2)3 log n.
√
r×m
Construct a matrix S ∈ R
, each element of which is a random variable which takes values ±1/ r with
equal probability. Let B̃ = SB. Then, if bi and b̃i denote the ith column of B and B̃, respectively, with
probability at least 1− n−β , and for all i = 1, ...n, (1− 21 )kbi k22 ≤ kb̃i k22 ≤ (1+ 12 )kbi k22 . Given B, it takes
O(nnz(B) log n) arithmetic operations to construct B̃. We will denote this procedure as B̃ = JLT (B, β).
Now, we are fully equipped to design the “high-probability” analog of Lemma 67:
Lemma 69. Given A ∈ Rm×n of rank ρ, a target rank 1 ≤ k < ρ, and 0 < δ, ǫ ≤ 1, there exists a
randomized algorithm that computes Z ∈ Rn×k with ZT Z = Ik and with probability at least 1 − δ− 1/n,
kA − AZZT k2F ≤ 3 (1 + ǫ) kA − Ak k2F . The proposed algorithm requires O nnz(A) log2 ( nδ ) + n ·
poly(k, ε, log( nδ )) arithmetic operations. We denote this procedure as Z = SparseSV DBoosting(A, k, ε, δ).
Proof. The idea is to use the algorithm in Lemma 67 and generate i.i.d. r = O(log( 1δ )) matrices
Z1 , Z2 , . . . , Zr , where Zi = SparseSV D(A, k, ε), and kA − AZi ZiT k2F ≤ (1 + ε)kA − Ak k2F , with
probability at least 0.99 for a single Zi .
Now, we could have chosen Z := Zi which minimizes kA − AZi ZiT k2F , and this would have
implied that kA − AZZT k2F ≤ (1 + ε)kA − Ak k2F , with probability 1 − δ (via a simple chernoff
argument), but evaluating the errors kA − AZi ZiT k2F is computationally intensive for the overall
running time that we would like to obtain.
Selection of Z. Instead, we use the Johnson Lindestrauss transform to speedup this computation
(evaluating the errors). Hence, for every Zi , we (implicitly) set Bi := A − AZi ZiT and then we
compute B̃i = JLT (Bi , 1) (see Lemma 68). Now, we choose the Zi that minimizes kB̃i k2F .
49

Correctness. For a fixed Zi , from Lemma 68 it follows that with probability 1 − 1/n it is: 21 kA −
AZi ZiT k2F ≤ kB̃i k2F ≤ 23 kA−AZi ZiT k2F , where the failure probability is over the randomness of S in
Lemma 68. Now, for some Zi with probability at least 0.99 it is kA− AZi ZiT k2F ≤ (1+ ε)kA − Ak k2F ,
where the failure probability is over the randomness in constructing Zi . Overall, given that n > 1,
with probability at least 0.49 (over the randomness of both S and Zi ; the failure probability follows
via a union bound) for each single Zi it is: kB̃i k2F ≤ 23 (1 + ε)kA − Ak k2F . Hence, since r = O(log 1δ )
and we have picked Z as the Zi that minimizes kB̃i k2F , it follows that with probability at least 1− δ,
kB̃i k2F ≤ 23 (1+ε)kA −Ak k2F . Combining with 21 kA−AZi ZiT k2F ≤ kB̃i k2F , and using a union bound,
it follows that with probability 1 − 1/n − δ: kA − AZi ZiT k2F ≤ 62 (1 + ε)kA − Ak k2F .

Running time. The following costs occur O(log( 1δ )) times: (i) O nnz(A) + (m + n) · poly(k, ε−1 )
to find Zi (via Lemma 67); (ii) O(nnz(A) log n + nk log(n)) to evaluate each cost, i.e., find kB̃i k2F .
9.1.3

Fast column sampling techniques

Section 8.1.1 presented the column sampling algorithms that we used in our distributed PCA
algorithm in Section 8. Next, we present the “fast” analogs of those algorithms. To design such
fast analogs we employ the sparse subspace embedding matrices in the previous section. All the
results in this section developed in [19]. There is a small difference in Lemma 70 below, hence we
present a short proof for completeness. Specifically, the original result in [19] has a constant failure
probability; here, via using standard arguments, we extended this result to a high probability
bound.
Lemma 70 (Input-Sparsity-Time Dual-Set Spectral-Frobenius Sparsification [19].). Let V ∈ Rw×k
be a matrix with w > k and V T V = Ik . Let E ∈ Rm×w be an arbitrary matrix. Let δ be a failure
probability parameter.
For i = 1, 2, . . . , r = O(log 1δ ), let Bi = EWiT ∈ Rm×ξ , where Wi ∈ Rξ×w is a randomly chosen
sparse subspace embedding with ξ = O(k2 /ε2 ) < w, for some 0 < ε < 1 (see Definition 65 and Lemma 66),
and run the algorithm of Lemma 57 with V, Bi , and some ℓ with k < ℓ ≤ w, to construct Si ∈ Rw×ℓ .
Choose some S from the Si ’s - see the proof for the details.
2


p 2

1+ε
· kEk2F .
Then, with probability at least 1 − δ, σk2 V T S ≥ 1 − k/ℓ and kESk2F ≤ 1−ε


The matrix S can be computed in O log 1δ · nnz(A) + ℓwk2 + mk 2 /ε2 + mℓ time. We denote this
procedure as S = BssSamplingSparse(V, A, r, ε, δ).
Proof. The algorithm constructs the Si ’s as follows, Si = BssSampling(V, Bi , r) - see Lemma 57.

p 2

From that lemma, we also have that σk2 V T Si ≥ 1 − k/ℓ and kBiT Sk2F ≤ kBiT k2F , i.e.,

kWi AT Si k2F ≤ kWi AT k2F . Since Wi is a subspace embedding, from Lemma 66 we have that with
probability at least 0.99 and for all vectors y ∈ Rn simultaneously, (1 − ε) kAT yk22 ≤ kWi AT yk22 .
Apply this r times for y ∈ Rn being columns from Si ∈ Rn×r and take a sum on the resulting
inequalities: (1 − ε) kAT Si k2F ≤ kWAT Si k2F ; also, apply this n times for the basis vectors in Rn
and take a sum on the resulting inequalities: kWA T k2F ≤ (1 + ε) kA T k2F . Combining all these
· kA T k2 .
inequalities together, we conclude that with probability at least 0.99, kAT Si k2F ≤ 1+ε
1−ε p F2

To summarize, for each Si and with probability at least 0.99, σk2 V T Si ≥ 1 − k/ℓ and
2

· kEk2F .
kESi k2F ≤ 1+ε
1−ε
50

Selection of S. Now that we have r matrices Si each of which satisfying the above bounds with
probability 0.99 we select one of them as follows: 1) we compute all σk2 V T Si and sort them from
the largest to the smallest value; 2) we compute all kAT Si k2F and sort them from the smallest to
the largest value; 3) we choose an Si which occurs in the top 2/3 fraction of each of the lists.

Correctness. Let Xi be an indicator random variable: Xi = 1 if the ith matrix Si satisfies σk2 V T Si ≥

p 2
P
1 − k/ℓ , and Xi = 0, otherwise. Then, Pr[Xi ] ≥ 0.99. Let X = ri Xi . The expected value of
X is µX = E [X] = 0.99r. A standard Chernoff bound applied on X gives: Pr[X ≤ (1 − α)µX ] ≤
µX α
e− 2 . Let (1 − α)0.99r = 0.75r, i.e., α = 0.14/0.99. Applying this to the Chernoff bound, we get:
Pr[X ≤ 0.75r] ≤ e−O(r) . Replacing r = O(log(1/δ)): Pr[X ≤ 0.75r] ≤ δ/2. This means precisely
p 2
 
that 3/4 of the fraction of the Si ’s satisfy with probability 1 − δ/2 : σk2 V T Si ≥ 1 − k/ℓ .
Similarly, let Yi be an indicator random variable: Yi = 1 if the ith matrix Si satisfies kESi k2F ≤
2
P
1+ε
· kEk2F , and Yi = 0, otherwise. Then, Pr[Yi ] ≥ 0.99. Let Y = ri Xi . The expected value of
1−ε
Y is µY = E [Y ] = 0.99r. A standard Chernoff bound applied on Y gives: Pr[Y ≤ (1 − α)µY ] ≤
µY α
e− 2 . Let (1 − α)0.99r = 0.75r, i.e., α = 0.14/0.99. Applying this to the Chernoff bound, we get:
Pr[Y ≤ 0.75r] ≤ e−O(r) . Replacing r = O(log(1/δ)): Pr[Y ≤ 0.75r] ≤ δ/2. This means precisely
2

1+ε
· kEk2F .
that 3/4 of the fraction of the Si ’s satisfy w.p. 1 − δ/2 : kESi k2F ≤ 1−ε
Now a simple union bound indicates that with probability at least 1 − δ, the following two
events happen simultaneously: X > 0.75r and Y > 0.75r. Since we select an element Si on top of
the two sorted lists, it follows that this element should have Xi = Yi = 1.


Running time. The following costs occur O(log( 1δ )) times: (i) O(nnz(A)) to find Bi . (i) O(ℓwk2 +

mξ) to find Si (via Lemma 57); and (iii) O(k 2 ℓ + mℓ) to evaluate each cost and, i.e., find σk2 V T Si
and kESi k2F . Additionally, O(log( 1δ ) log log( 1δ )) time in total is spent in sorting.
The following Lemma is the “fast” analog of Lemma 58. The failure probability in the lemma is
constant, but this is sufficient for our purposes since we apply this lemma only once in the global
column sampling step of the distributed PCA algorithm in Section 9.2.
Lemma 71 (Input-Sparsity-Time constant factor column-based matrix reconstruction; Lemma 6.3
2
in [19]). Given matrix G ∈ Rm×α of rank ρ and a target
 rank k , there exists a randomized algorithm
−1
that runs in O nnz(A) · log α + m · poly(log α, k, ε ) time and selects c = 4k columns of G to form a
matrix C ∈ Rm×c such that with probability at least 0.69:
kG − CC† Gk2F ≤ 4820 ·

rank
X(G)

σi2 (G).

i=k+1

We denote this procedure as C = DeterministicCssF robeniusSparse(G, k, c).
Finally, the lemma below presents the “fast” analog of Lemma 59. The failure probability in
the lemma is, again, constant; we have not attempted to obtain a high probability bound since we
employ this lemma only once in the algorithm in Section 9.2.
2

The original Lemma 6.3 in [19] has the assumption that k < ρ, but this assumption can be dropped having the
result unchanged. The only reason the assumption k < ρ exists is because otherwise column subset selection is trivial.

51

Lemma 72 (Input-sparsity-time Adaptive Sampling [19]). Given A ∈ Rm×n and V ∈ Rm×c1 (with
c1 ≤ n, m), define the residual matrix Ψ = A − VV† A ∈ Rm×n . Let Ψ̃ = JLT (Ψ, 1). For j = 1, . . . , n,

let pj be a probability distribution such that, for some positive constant β ≤ 1, pj ≥ βkΨ̃
(j)

(j) 2
k2 /kΨ̃k2F ,

where Ψ̃ is the j-th column of the matrix Ψ̃. Sample c2 columns from A in c2 i.i.d. trials, where in each
trial the j-th column is chosen with probability pi . Let C2 ∈ Rm×c2 contain the c2 sampled columns and
let C = [V C2 ] ∈ Rm×(c1 +c2 ) contain the columns of V and C2 . Then, for any integer k > 0, and with
probability 0.9 − n1
2
2
kA − ΠF
C,k (A)kF ≤ kA − Ak kF +

30 · k
kA − VV† Ak2F .
β · c2

Given A and V, the algorithm takes O(nnz(A) log n + mc1 log n + mc21 ) time to find C2 . We denote this
sampling procedure as C2 = AdaptiveColsSparse(A, V, c2 , β).
9.1.4

Fast Low-rank matrix approximations within a subspace

Finally, we present the fast analog of the result in Section 8.1.3. The failure probability in the
lemma is, again, constant; we have not attempted to obtain a high probability bound since we
employ this lemma only once in the algorithm in Section 9.2.
Lemma 73. Let A ∈ Rm×n be the input matrix and V ∈ Rm×c be the input subspace. We further assume
that for some rank parameter k < c and accuracy parameter 0 < ε < 1 :
2
2
kA − ΠF
V,k (A)kF ≤ (1 + O(ε))kA − Ak kF .

Let V = YΨ be a qr decomposition of V with Y ∈ Rm×c and Ψ ∈ Rc×c . Let Ξ = Y T AW T ∈ Rc×ξ ,
where W T ∈ Rn×ξ with ξ = O(c/ε2 ), is a sparse subspace embedding matrix (see Definition 65 and
Lemma 66). Let ∆ ∈ Rc×k contain the top k left singular vectors of Ξ. Then, with probability at least 0.99,
kA − Y∆∆T Y T Ak2F ≤ (1 + ε)kA − Ak k2F .
Y and ∆ can be computed in O(nnz(A) + mcξ) time. We denote this procedure as
[Y, ∆] = ApproxSubspaceSV DSparse(A, V, k, ε).
Proof. This result was proven inside the proof of Theorem 1.5 in [36]. Specifically, the error bound
proven in [36] is for the transpose of A (also Y, ∆ are denoted with U, V in [36]). The only requirement for the embedding matrix W (denoted with P in the proof of Theorem 1.5 in [36]) is to be a
subspace embedding for Y T A, in the sense that W is a subspace embedding for A in Lemma 66.
Since our choice of W with ξ = O(c2 /ε2 ) satisfies this requirement we omit the details of the proof.
The running time is O(nnz(A)+mcξ): one can compute (i) Y in O(mc2 ); (ii) Ξ in O(nnz(A)+mcξ);
and (iii) ∆ in O(cξ min{c, ξ}). The failure probability 0.01 from Lemma 66.

9.2 Detailed description of the algorithm
This algorithm is very similar to the algorithm in Section 8.2; we only replace the parts 2-(a), 3-(a),
and 4-(b) with faster procedures.
Input:
1. A ∈ Rm×n partitioned column-wise A = A1

A2

2. rank parameter k < rank(A)

52


P
. . . As ; for i = 1 : s, Ai ∈ Rm×wi ; i wi = n.

3. accuracy parameter ε > 0
4. failure probability δ
Algorithm
1. Local Column Sampling
(a) For each sub-matrix Ai ∈ Rm×wi , compute Ci ∈ Rm×ℓ containing ℓ = O(k) columns from
Ai as follows: Ci = Ai Si . Here, Si has dimensions wi × ℓ and is constructed as follows:
Si = BssSamplingI(Ai , 4k, 1/2, 0.01/s) (see Lemma 56).
(b) Machine i sends Ci to the server.
2. Global Column Sampling
(a) Server constructs m × (s · ℓ) matrix G containing (s · ℓ) actual columns from A as follows:
G = C1 C2 . . . Cs . Then, server constructs C ∈ Rm×c1 via choosing c1 = 4k columns
from G as follows: C = DeterministicCssF robeniusSparse(G, k, c1 ) (see Lemma 71).

(b) Server sends C to all the machines.
3. Adaptive Column Sampling

(a) Server initializes the random seed and communicates it to all the machines. Each machine constructs the same S ∈ Rr×m , with r = (1/2)4+2
2 −(1/2)3 log n, each element of which is a random vari√
able which takes values ±1/ r with equal probability. Machine i finds Ψ̃i = SAi − SCC† Ai ∈
Rm×wi and then computes βi as it was described in Lemma 61. Machine i sends βi to server.
(b) Server computes probability distribution gi = Pβiβi . Server samples i.i.d. with replacement
i
⌈ 50k/ε ⌉ samples (machines) from gi . Then, server determines numbers ti (i = 1, 2, . . . , s),
where ti is the number of times the ith machine was sampled. It sends the ti ’s to the machines.
(c) Machine i computes probabilities qji = kxk22 /kΨ̃i k2F (j = 1 : wi ), where x is the jth column of
Ψ̃i . And now machine i samples ti samples fromP
it’s local probability distribution and sends
the corresponding columns to the server. Let c2 = i ti = ⌈ 50k/ε ⌉.

(d) Server collects the columns and assigns
them to Ĉ ∈ Rm×c2 . Server constructs C̃ to be the

m × (c1 + c2 ) matrix: C̃ = C; Ĉ . Let c = c1 + c2 = 4k + ⌈ 50k/ε ⌉.

4. Rank-k matrix in the span of C̃

(a) Server sends C̃ to all the machines and each machine computes (the same) qr factorization of
C̃: C̃ = YR where Y ∈ Rm×c has orthonormal columns and R ∈ Rc×c is upper triangular.

(b) Server initializes the random seed and sends this to each machine, such that each machine
generates the same matrix Ψ ∈ Rξ×n , for ξ = O(c2 /ε2 ), (see Definition 65). Machine i generates
Di ∈ Rn×wi such that implicitly D = D1 D2 . . . Ds is an n × n matrix each element of
ξ×wi
.
which is ±1, with probability 1/2. Each machine constructs implicitly
i ∈ R
 Wi = Ψ · Dξ×n
and W
Implicitly all machines together generate W = W1 W2 . . . Ws , with W ∈ R
T

is a subspace embedding as in Definition 65. Machine i computes Hi = C̃ (Ai WiT ) ∈ Rc×ξ .
Machine i sends Hi to the server.
Ps
(c) Server computes Ξ = i=1 Hi ∈ Rc×ξ and sends this back to all the machines. Now machines
compute Ξ := R−1 · Ξ(:= Y T Ai WiT ), and then they compute ∆ ∈ Rc×k to be the top k left
singular vectors of Ξ. Each machine computes U = Y · ∆ ∈ Rm×k .

53

Discussion. A few remarks are necessary for the last two stages of the algorithm. The third stage
(adaptive column sampling), implements the adaptive sampling method of Lemma 72. To see this
note that each column in A is sampled with probability
qji · gi ≥

1
· ||x||22 /kΨ̃k2F ,
2

where x is this column in Ψ̃. This follows from Lemma 60. Overall, this stage constructs Ĉ such
that Ĉ = AdaptiveColsSparse(A, C, c2 , 1/2).
The last stage in the algorithm implements the algorithm in Lemma 73. To see this, note that
W satisfies the properties in the lemma. Hence,
[Y, ∆] = ApproxSubspaceSV DSparse(A, C̃, k, ε).

9.3 Main result
The theorem below analyzes the approximation error, the communication complexity, and the
running time of the previous algorithm.
Theorem 74. The matrix C̃ ∈ Rm×c with c = O(k/ε) columns of A satisfies w.p. 0.59 −

s+1
n



rank
X(A)
(A)k2F ≤ (1 + O(ε)) · 
kA − C̃C̃ Ak2F ≤ kA − ΠF
σi2 (A) .
C̃,k
†

− 2δ:
(15)

i=k+1

The matrix U ∈ Rm×k with k orthonormal columns satisfies w.p. 0.58 −


kA − UUT Ak2F ≤ (1 + O(ε)) · 

rank
X(A)
i=k+1

s+1
n

− 2δ:



σi2 (A) .

(16)

Let each column of Ahas at most φ non-zero elements. Then, the communication cost of the algorithm is
O skφε−1 + sk3 ε−5 and the running time is
 ns 

 ns 
O nnz(A) · log2
+ (m + n)s · poly(k, ε−1 , log
.
δ
δ

9.4 Proof of Theorem 74
9.4.1

Proof of Eqn. 15

From Lemma 72, we have that with probability at least 0.9 − n1 :
†

(A)k2F ≤
kA − C̃C̃ Ak2F ≤ kA − ΠF
C̃,k

rank
X(A)
i=k+1

σi2 (A) +

60ε
· kA − CC† Ak2F .
50

(17)

According to Lemma 64, G is a (4k, 1/2)-projection-cost preserving sketch of A. It means that there
exists a constant c which is independent from projection matrix P and we have the following:
3
1
kA − PAk2F ≤ kG − PGk2F + c ≤ kA − PAk2F
2
2
54

For convenience, we use PG to denote a rank-k projection matrix which provides the best rank-k
approximation Gk = PG G to matrix G. We also use PA to denote a rank-k projection matrix
which provides the best rank-k approximation Ak = PA A to matrix A. With probability at least
0.68, the following will be held,
1
kA − CC† Ak2F ≤ kG − CC† Gk2F + c
2
≤ 4820kG − PG Gk2F + c
≤ 4820kG − PA Gk2F + c


≤ 4820 kG − PA Gk2F + c
4820
≤
kA − PA Ak2F
2
The first inequality is due to Lemma 64. G is a (4k, 1/2)-projection cost preserving sketch of A
and C has rank at most 4k. The second inequality is based on Lemma 71. The third inequality is
because PG provides the best rank-k approximation to G. The fourth inequality used the fact c is
non-negative. The fifth inequality held because G is a (4k, 1/2)-projection cost preserving sketch
of A. Thus,
kA − CC† Ak2F ≤ 4820kA − Ak k2F
Concluding the proof. Combining with Equation 17, we have:
†

kA − C̃C̃ Ak2F ≤ kA − ΠF
(A)k2F ≤
C̃,k
9.4.2

rank
X(A)
i=k+1

σi2 (A) + O(ε) ·

rank
X(A)

σi2 (A).

(18)

i=k+1

Proof of Eqn. 16

We would like to apply Lemma 73 for the matrix U ∈ Rm×k in the algorithm. Note that we already
proved that the matrix C̃ in the algorithm satisfies with probability 0.59 − s+1
n − 2δ:
kA − ΠF
(A)k2F ≤ (1 + O(ε))kA − Ak k2F .
C̃,k
Lemma 73 and a simple union bound conclude the proof.
9.4.3

Running time

1. Local Column Sampling: O(nnz(A) log(ms) + m · poly(k, ε, log(s))) arithmetic operations in
total.
(a) O(nnz(Ai ) log(ms) + m · poly(k, ε, log(s))) for each Si

(b) -


2. Global Column Sampling: O nnz(A) · log k + m · poly(k, ε−1 ) .

(a) O nnz(A) · log k + m · poly(k, ε−1 ) to construct C - from Lemma 71.
(b) -

55

3. Adaptive Sampling: O(nnz(A) · log n + m · s · poly(k, ε−1 , log n)).
(a) First, we analyze the cost to compute the matrix Ψ̃i . O(k2 m) for each C† locally; then,
we compute (SAi ) − (((SCi )C†i )Ai ). The costs are: SAi := D takes O(nnz(Ai )) time,
SCi := G takes O(nnz(Ai ) log(n)) time, GC†i := H takes O(m log(n)k/ε2 ) time, HA :=
L takes O(nnz(Ai ) log(n)) time, and D − L takes O(n log n) time. So, the total time to
compute one Ψ̃i is O(nnz(Ai ) log n + m · poly(k, ε−1 , log n)).
Also, there is a cost of computing ψi , which is O(log n · wi ) arithmetic operations.

(b) O(s) in total.

(c) O(k/ε) in total using the method of [51] locally.
(d) O(k/ε).
4. Rank-k matrix in span(C̃): O(nnz(A) + m · s · poly(k, ε−1 )) arithmetic operations in total.
(a) O(smk 2 /ε2 ) in total.
(b) O(sn) in total to generate s times the same Ψ. Then another O(n) in total to generate
the Di ’s. For all Ai WiT we need O(nnz(A)) operations and then another O(mk 3 ε−7 ) for
each Y T · (Ai WiT )
(c) O(sk 3 ε−5 ) to find Ξ; then another O(sk3 ε−3 + sk4 ε−6 )) to update Ξ locally in each
machine and another O(sk 4 ε−6 ) to find ∆; and O(smk 2 /ε) to find U (s times).

9.4.4

Communication Complexity

Assume that we can represent each element in the input matrix A with at most b bits. Also, we
assume that one word has length b bits and b = O(log(mns/ε)).
1. Local Column Sampling: O(skφ) words.
(a) (b) O(skφ) elements of A.
2. Global Column Sampling:O(skφ) words.
(a) (b) O(skφ) elements of A.
3. Adaptive Sampling: O(s + φk/ε) words.
(a) s integers each of which is representable with O(log k+log log(mns/ε)) (from Lemma 61).
(b) s integers (the ti ’s) each with magnitude at most n, hence representable with b bits.
(c) O(φk/ε) elements of A.
(d) 4. Best rank-k matrix in the span of C̃: O(sφk/ε + sk 3 ε−5 ) words.
(a) O(skφε−1 ) elements of A.
(b) O(sk 3 ε−5 ) numbers each of which can be represented with b bits.
(c) O(sk 3 ε−5 ) numbers each of which can be represented with b bits.
In total the communication complexity is O(skφ/ε + sk 3 /ε5 ) words.
56

10 Lower bounds
Only in Section 10.1, we discuss the result in arbitrary partition model. In Section 10.2, Section 10.3
and Section 10.4, we discuss results in column partition model which provides stronger lower
bounds than in arbitrary partition model.

10.1 Dependence on Ω(ε−2 ) bits for distributed PCA in arbitrary partition model
Lemma 75 (Theorem 8 in [54]). There are two machines, and each of them hold a matrix Ai ∈ Rm×n .
Let A = A1 + A2 . Suppose ε2 ≤ 2n, the communication of computing a rank-k matrix U with constant
probability such that
kA − UUT Ak2F ≤ (1 + ε) · kA − Ak k2F
needs Ω(1/ε2 ) bits.

Proof. (from [54]). The reduction is from GHD problem.
Lemma 76 (Theorem 1.1 in [20]). Each of two machines has an n-bit vector. There is a promise: either
√
√
the inputs are at Hamming distance less than n/2 − c n or greater than n/2 + c n. If they want to
distinguish these two cases with constant probability, the communication required Ω(n) bits.
Without loss of generality, we assume 1/ε2 is an integer, and machine 1 and machine 2 have
2
x, y ∈ {−1, 1}1/ε respectively. There is a promise that either xT y < −2/ε or xT y > 2/ε holds.
Consider about the following protocol:
1. Machine 1 constructs A1 ∈ R(k+1)×(1/ε


xT ε √0
 0
2


A1 = 
0
 0
 ...
...

0
0

0
√ 0

2 +k)

, and machine 2 constructs A2 ∈ R(k+1)×(1/ε

...
...

0
0

...

...
...

0

...

0
...
√

2(1+ε)
ε



2(1+ε)
ε













A2 = 



y T ε 0 0 ... 0
0
0 0 ... 0
0
0 0 0 0
... ... ... ... ...
0
0 0 ... 0

2 +k)

:








2. The server computes U such that
kA − UUT Ak2F ≤ (1 + ε) · kA − Ak k2F
and sends U to both of two machines, where A = A1 + A2 .
3. Let P = Ik+1 − UUT . Each machine construct v = P(1) /kP(1) k2
4. Each machine checks whether v12 < 21 (1 + ε). If yes, return the case xT y > 2/ε. Otherwise
return the case xT y < 2/ε
Suppose the server successfully computes U. Since the rank of A is at most k + 1:
kA − UUT Ak2F = kPAk2F
= kvAk22

= vAAT v T
57

An observation is that





AAT = 



kx + yk22 ε2 0
0
2
0
0
...
...
0
0

0
0
2(1+ε)
ε2

...
0

...
...
...
...
...

0
0
0
...
2(1+ε)
ε2

We have kA − Ak k2F ≤ 2. Then,

k+1

2(1 + ε) X 2
vi =
ε2
i=3









2(1 + ε)
(1 − v12 − v22 )
ε2

< 2(1 + ε)

We have v12 + v22 > 1 − ε2 . When xT y > 2/ε, kx + yk22 ε2 ≥ 2 + 4ε
(1 + ε)kA − Ak k2F = 2(1 + ε)

> v12 (2 + 4ε) + 2v22

> 2 − 2ε2 + 4εv12
So, v12 < 12 (1 + ε). When xT y < −2/ε, kx + yk22 ε2 ≤ 2 − 4ε
(1 + ε)kA − Ak k2F = kx + yk22 ε2 (1 + ε)

> v12 kx + yk22 ε2 + 2v22

= 2(v12 + v22 ) − (2 − kx + yk22 ε2 )v12
> 2(1 − ε2 ) − (2 − kx + yk22 ε2 )v12

So,
v12 >

1
2(1 − ε2 ) − kx + yk22 ε2 (1 + ε)
≥ (1 + ε)
2
2
2
2 − kx + yk2 ε

Therefore, machines can distinguish these two cases. The only communication is computing U.
This cost shoulde be the same as Ω(1/ε2 ) bits lower bound of the gap hamming distance problem.

10.2 Lower bounds for distributed PCA on dense matrices
This section provides a communication cost lower bound for the Distributed PCA problem of
Definition 17. Specifically, we describe the construction of an m × n matrix A (hard instance), and
formally argue that for this A, any k ≤ 0.99m, and any error parameter C with 1 < C < poly(skm),
if there exists some algorithm to construct an m × k matrix U such that, with constant probability,
kA − UUT Ak2F ≤ C · kA − Ak k2F , then this algorithm has communication cost Ω(skm) words.
10.2.1

Preliminaries

We use the notation Gk,m to denote the set of k-dimensional subspaces of Rm , which we identify
with corresponding projector matrices Q ∈ Rm×m , with rank(Q) = k ≤ m, onto the subspaces,
i.e.,
Gk,m = {Q ∈ Rm×m : Q2 = Q, rank(Q) = k ≤ m.}
58

We also need a description of a subset from Gk,m:
Cδk (P) = {Q ∈ Gk,m : kP − Qk2 ≤ δ}.
This is the set of those projectors Q ∈ Gk,m which differ in operator norm by at most δ from
another fixed projector P ∈ Gk,m . In our analysis below we also need a result from [38]:

Theorem 77. (Net Bound - Corollary 5.1 of [38]) For any m, k, δ > 0, there is a family N = {P1 , . . . , PN },
N = 2Ω(k(m−k) log(1/δ)) , where Pi ∈ Gk,m and Cδk (Pi ) ∩ Cδk (Pj ) = ∅, for all i 6= j.

Theorem 77 proves the existence of a large, but finite, set N of projection matrices, such that if
one considers the ball of matrices of operator norm at most δ centered at each matrix in N , then
these balls are disjoint. Our hard instance matrix A is constructed from some member in N .
10.2.2

Hard instance construction

Recall that in the column-partition model of Definition 15 there are s servers holding matrices
A1 , . . . , As , respectively,
P where Ai has m rows and some subset of wi 3columns of some m × n
matrix A. Notice that wi = n. First of all, for arbitrary m, s we assume that n ≥ sm. Below, we
describe a specific construction for a matrix A.
First of all, we set δ = 1/(skm), the parameter to be used in Theorem 77. Next, fix a family N
of subspaces of Gk,m with the guarantees of Theorem 77. Let PR = RRT be a uniformly random
member of N , where R ∈ Rm×k has orthonormal columns (any projector can be expressed as
such a matrix, where the columns of R span the subspace that P projects onto). The entries of
A1 are equal to the entries of R, each rounded to the nearest integer multiple of 1/B, where
B = poly(skm) is a large enough parameter specified later. We denote this rounded matrix with
k
R̃ ∈ Rm×k . So, if αij is the (i, j)th entry in A1 , then αij = kmin
B , with kmin = argmin |Rij − B |. Each
k∈Z

Ai for i = 2, 3, . . . , s − 1 is 1/B times the m × m identity matrix. Finally, As is the m × t matrix of
all zeros with t = n − (s − 1)m − k. I.e.,

A = R̃ B1 Im B1 Im . . . B1 Im 0m×t .

10.2.3

Intermediate results

First, we give a bound regarding the best rank k approximation for the above matrix A.
Lemma 78. For the matrix A ∈ Rm×n in Section 10.2.2, and any k ≤ 0.99m : kA − Ak k2F <

sm
.
B2

Proof. We have,
kA − Ak k2F ≤ kA − A1 A†1 Ak2F
X
=
k(I − A1 A†1 )Ai k2F
i>1

≤

X
i>1

kAi k2F

m
= (s − 2) · 2
B
sm
<
,
B2
As, otherwise we can choose a value s′ < s so that n − m ≤ s′ m ≤ n and apply the argument of Theorem 82 with
s replaced with s′ . Then, the lower bound in Theorem 82 will then be Ω(s′ mk), which assuming n ≥ 2d, is an Ω(n)
communication lower bound
3

59

where the first inequality uses the fact that A1 A†1 is a matrix of rank at most k, the first equality
uses the fact that (I − A1 A†1 )Ai is the all-zeros matrix, the second inequality uses that a projector
cannot increase a unitarily invariant norm, and the third inequality follows by construction.
Next, we bound, in the operator norm, the difference of A1 from the matrix R. The lemma
follows from the fact that A1 ∈ Rm×k is obtained by R ∈ Rm×k by rounding each entry to the
nearest integer multiple of 1/B, and then summing the squared differences across all km entries.
Lemma 79. (Precision Lemma) kA1 − Rk2 ≤
Proof. kA1 − Rk22 ≤ kA1 − Rk2F ≤

√
km
B .

km
.
B2

Next, we prove a pure linear algebraic result. The following lemma captures the fact that if
some P ∈ Gk,m is close to QP (in Frobenius norm) for Q ∈ Gk,m , then P is also close to Q (in
Frobenius norm). We will use this lemma later for P = PR .
Lemma 80. (Error Measure) Let P, Q ∈ Gk,m with kP − QPk2F ≤ ∆. Then kP − Qk2F ≤ 2∆.
Proof. By the matrix Pythagorean theorem,
k = kQk2F = kQPk2F + kQ(I − P)k2F .

(19)

Also by the matrix Pythagorean theorem,
k = kPk2F = kQPk2F + k(I − Q)Pk2F ,
and so using the premise of the lemma, kQPk2F ≥ k − ∆. Combining with (19),
kQ(I − P)k2F ≤ ∆.

(20)

Hence,
kP − Qk2F = k(P − Q)Pk2F + k(P − Q)(I − P)k2F
= kP − QPk2F + kQ − QPk2F
≤ ∆ + kQ(I − P)k2F

≤ 2∆,

where the first equality follows by the matrix Pythagorean theorem, the second equality uses
P2 = P (since P is a projector), the third equality uses the bound in the premise of the lemma, and
the fourth inequality uses (20).
10.2.4

Main Argument

Before presenting the main theorem, we give an intermediate technical lemma.
Lemma 81. (Implication of Correctness) Let A ∈ Rm×n be the hard instance matrix in Section 10.2.2 and
let A be distributed in s machines in the way it was described in Section 10.2.2:

A = R̃ B1 Im B1 Im . . . B1 Im 0m×t .

Suppose n = Ω(sm), k < rank(A), and 1 < C < poly(skm). A is revealed to the machines by just
describing the entries in A and without specifying any other detail regarding the construction. Given
this A and assuming further that each machine knows some projector matrix Q ∈ Rm×m with rank at
most k such that kA − QAkF ≤ CkA − Ak kF , there is a deterministic algorithm (protocol) which, upon
termination, leaves on each machine the matrix R ∈ Rm×k which was used to construct A1 .
60

1
Proof. First, we prove that kRRT − Qk2 < 2skm
. Then, using this bound, we describe a protocol
that deterministically reveals the matrix R which was used to construct A1 .
, which in
Using Lemma 78 and the premise in the lemma, it follows that: kA − QAk2F ≤ C sm
B2
particular implies that

kA1 − QA1 k2F ≤ C

sm
.
B2

(21)

We further manipulate the term kA1 − QA1 kF as follows:
kA1 − QA1 kF = k(I − Q)A1 kF

= k(I − Q)RRT + (I − Q)(A1 − RRT )kF

≥ k(I − Q)RRT kF − k(I − Q)(A1 − RRT )kF

≥ k(I − Q)RRT kF − kI − QkF k(A1 − RRT )k2
√
√
km √
T
≥ k(I − Q)RR kF −
( m + k),
B

where the first inequality is the triangle inequality, the second inequality uses sub-multiplicativity,
and the third inequality uses the triangle inequality, i.e., that kI − QkF ≤
kIkF + kQk
, and Lemma
√ F
√
T
Csm
2m k
79. Combining this bound with (21), it follows that k(I − Q)RR kF ≤ B + B . At this point
we would like to apply Lemma 80 to conclude that kRRT − QkF is small. To do so, we need RRT
and Q to be of rank k. While R has rank k by construction, Q may not. However, increasing the
rank of Q cannot increase the error kA−QAkF . Hence, if rank(Q) < k, given Q, the protocol in the
premise of the lemma allows each server to locally add standard basis vectors to the column space
of Q until the column space of Q becomes k. This involves no communication and each server
ends up with the same new setting of Q, which for simplicity we still denote with Q. Therefore,
rank(Q) = k. Applying Lemma 80, it now follows that
√ !
√
√
Csm 2m k
T
T
+
.
kRR − Qk2 ≤ kRR − QkF ≤ 2
B
B
1
.
By setting B to be a sufficiently large poly(skm), we have kRRT − Qk2 < 2skm
i
Therefore, given Q, by Theorem 77 there is a unique matrix P in N for which Q ∈ Cδk (Pi ),
1
1
with δ = skm
. Indeed, since kR − Qk2 ≤ skm
< δ/2, there cannot be two matrices Pi and Pj for
i 6= j for which Q ∈ Cδk (Pi ) and Q ∈ Cδk (Pj ) since then by the triangle inequality kPi − Pj k2 < δ,
contradicting the construction of N . Hence, it follows that each machine can deterministically
identify this Pi , by enumeration over N , and therefore compute R, via, for example, an SVD.

We are now fully equipped to present the main argument about the communication cost lower
bound for distributed PCA.
Theorem 82. (Main) Let A ∈ Rm×n be the hard instance matrix in Section 10.2.2 and let A be distributed
in s machines in the way it was described in Section 10.2.2:

A = R̃ B1 Im B1 Im . . . B1 Im 0m×t .

Suppose n = Ω(sm), k ≤ .99m, and 1 < C < poly(skm). Assume that there is an algorithm (protocol)
which succeeds with probability at least 2/3 in having the i-th server output QAi , for all i, where Q ∈
Rm×m is a projector matrix with rank at most k such that kA − QAkF ≤ CkA − Ak kF . Then, this
algorithm requires Ω(skm log(skm)) bits of communication.
61

Further, the bound holds even if the input matrices Ai to the servers have all entries which are integer
multiples of 1/B, for a value B = poly(skm), and bounded in magnitude by 1, and therefore all entries can
be specified with O(log(skm)) bits.
Note that assuming a word size of Θ(log(skm)) bits, we obtain an Ω(skm) word lower bound.
Proof. If the protocol succeeds, then, since each server outputs QAi , and Ai is 1/B times the
identity matrix for all i = 2, 3, ...s − 1, each of those servers can compute Q. One of those servers
can now send this Q to the first and the last server. It follows by Lemma 81, that each server can
identify PR and R.
Letting E be the event that the protocol succeeds, and letting Πi be the ordered sequence of all
incoming and outgoing messages at the i-th server, it follows that
I(Πi ; Q | E) = Ω(log |N |) = Ω(km log(smk)),
where I(X; Y | F) = H(X | F)−H(X | Y, F) is the mutual information between random variables
X and Y conditioned on F. To see the first equality in the above, observe that I(Πi ; Q | E) = H(Q |
E) − H(Q | Πi , E), and conditioned on E, there is a uniquely determined matrix PR each server
identifies, which is uniform over a set of size |N |, and so H(Q | E) ≥ H(PR ) = Ω(log2P
|N |).
Here, H(X) is the Shannon entropyP
of a random variable X, given by H(X) = x Pr[X =
x] log(1/ Pr[X = x]), and H(X | Y ) = y Pr[Y = y]H(X | Y = y) is the conditional Shannon
entropy.
Hence, if Z is an indicator random variable which is 1 if and only if E occurs, then using that
I(X; Y ) ≥ I(X; Y | W ) − H(W ) for any random variables X, Y, and W , and that I(X; Y | W ) =
P
w Pr[W = w]I(X; Y | W = w), we have
I(Πi ; Q) ≥ I(Πi ; Q | Z) − 1

≥ I(Πi ; Q | Z = 1) Pr[Z = 1] − 1
= Ω(km log(smk)),

mapping X to Πi , Y to Q, and W to Z in the mutual information bound above. It follows that
Ω(km log(smk)) ≤ I(Πi ; Q) ≤ H(Πi ) ≤ EPR ,rand [|Πi |],
where |Πi | denotes the length, in bits, of the sequence Πi , and rand is the concatenation of the
private randomness of all s players Note that the entropy of Πi is a lower bound on the expected encoding length of |Πi |, by the Shannon coding theorem. By linearity of expectation,
P
kΠi | = Ω(skm log(smk)), which implies there exists an R and setting of rand for which
Pi ER,rand
i
i kΠ | = Ω(skm log(smk)). It follows that the total communication is Ω(smk log(smk)) bits.

10.3 Lower bounds for distributed PCA on sparse matrices

Applying the previous theorem with m = φ gives a communication lower bound for sparse matrices.
Corollary 83. Let A ∈ Rφ×n be the hard instance matrix in Section 10.2.2 (with m replaced with φ in
noting the number of rows in A):

A = R̃ B1 Iφ B1 Iφ . . . B1 Iφ 0φ×t .

Suppose n = Ω(sφ), k ≤ .99φ, and 1 < C < poly(skφ). Now, for arbitrary m consider the matrix
Â ∈ Rm×n which has A in the top part and the all-zeros (m − φ) × n matrix in the bottom part. Â admits
the same column partition as A after padding with zeros:
62

Â =



R̃

1
B Iφ

1
B Iφ

0(m−φ)×k 0(m−φ)×φ 0(m−φ)×φ


1
...
0φ×t
B Iφ
.
. . . 0(m−φ)×φ 0(m−φ)×t

We denote the new sub-matrices with Âi ∈ Rm×wi .
Assume that there is an algorithm (protocol) which succeeds with probability at least 2/3 in having the
i-th server output Q̂Âi , for all i, where Q̂ ∈ Rm×m is a projector matrix with rank at most k such that
kÂ − Q̂ÂkF ≤ CkÂ − Âk kF . Then, this algorithm requires Ω(skφ log(skφ)) bits of communication.
Note that assuming a word size of Θ(log(skφ)) bits, we obtain an Ω(skφ) word lower bound.

10.4 Lower bounds for distributed column-based matrix reconstruction
In this section we develop a communication lower bound for the Distributed Column Subset Selection Problem of Definition 19.
10.4.1

A Net of discretized Deshpande-Vempala hard instances

We start by describing the construction of a special sparse matrix from the work of Desphande and
Vempala [27], which was further studied by [17]. This matrix gives a lower bound (on the number
of columns need to be selected in order to achieve a certain accuracy) for the standard column
subset selection problem (see the following Fact). Suppose φ ≥ 2k/ǫ. Consider a (φ + 1) × φ matrix
B for which the i-th column is e1 + ei+1 , where ei is the i-th standard basis vector in Rφ+1 . Now
define an ((φ + 1)k) × (φk) matrix A with k blocks along the diagonal, where each block is equal
to the matrix B. Let n = (φ + 1)k.
Fact 84. (Proposition 4 in [27], extended in Theorem 37 of [17]) Let A be as above. If P is an n×n projection
operator onto any subset of at most k/(2ǫ) columns of A, then kA − PAk2F > (1 + ǫ)kA − Ak k2F .
To prove a communication lower bound for distributed column subset selection, we will need
to transform the above matrix A by multiplying by a rounded orthonormal matrix, as follows. We
also need to prove a similar lower bound as in the previous fact for this transformed matrix.
Lemma 85. Let A, k, φ, 1/ε, and n be as above. For any n × n orthonormal matrix L, let L̃ denote
the matrix formed by rounding each of the entries of L to the nearest integer multiple of M , where M ≤
1/((φ + 1)k)c , for a sufficiently large constant c > 0. If P is an n × n projection operator onto any subset
of at most k/(6ǫ) columns of L̃A, then kL̃A − PL̃Ak2F > (1 + ǫ)kL̃A − [L̃A]k k2F .
Proof. We will show that, assuming 1/ε < n/3, if P is an n × n projection operator onto any subset
of at most k/(2ǫ) columns of L̃A, then kL̃A − PL̃Ak2F > (1 + ǫ/3)kL̃A − [L̃A]k k2F . The lemma will
then follow by replacing ε with 3ε.
Suppose S is a subset of k/(2ǫ) columns of L̃A for which
kL̃A − PL̃Ak2F ≤ (1 + ǫ/3)kL̃A − [L̃A]k k2F ,

(22)

where P is the projection onto the columns in S and [L̃A]k is the best rank-k approximation to L̃A
The proof strategy is to relate the left hand side of (22) to kA − QAk2F where Q is the projection
onto the columns of A indexed by the same index set S, and simultaneously to relate the right
hand side of (22) to kA − Ak k2F .
We start by looking at kL̃A − [L̃A]k kF . Notice that L̃A = LA + EA, where E is a matrix whose
entries are all bounded in magnitude by M . Therefore, L̃A = LA + F, where F = EA and
1/2
1/2
kFk2 ≤ kEkF kAkF ≤ n2 M 2
· (2(n − 1))1/2 ≤ 2n3 M 2
,
63

where we have used kAkF ≤ (2(n − 1))1/2 . By Weyl’s inequality (Corollary 7.3.8 of [34]), for all
i ≤ φk,
1/2
|σi (LA) − σi (LA + F)| ≤ kFk2 ≤ 2n3 M 2
,

where σi (X) denotes the i-th singular value of a matrix X. Since L is orthonormal, σi (LA) = σi (A)
for all i. It follows that
X
kL̃A − [L̃A]k k2F =
σi2 (L̃A)
i>k

X
1/2 2
=
σi (LA) ± 2n3 M 2
i>k

X
1/2 2
=
σi (A) ± 2n3 M 2
i>k

= kA − Ak k2F ± O(nkAk2 n3 M 2

= kA − Ak k2F ± O(1/n2 ),

1/2

+ n4 M 2 )
(23)

where the final line follows for a sufficiently small choice of M ≤ 1/nc , using that kAk2 ≤ kAkF =
√
O( n).
We now look at kL̃A − PL̃AkF . By the triangle inequality,
kL̃A − PL̃AkF ≥ kLA − PLAkF − kEAkF − kPEAkF

≥ kLA − PLAkF − kEkF kAkF − kPk2 kEkF kAkF

= kA − LT PLAkF − O(1/n3 ),

(24)

where the equality uses that multiplying by LT preserves norms, as well as that kEkF ≤ 1/poly(n)
for an arbitrarily small poly(n) by making M ≤ 1/nc sufficiently small, whereas kPk2 ≤ 1 and
√
kAkF = O( n). We claim LT PL = Q, where Q is the projection onto the columns in A indexed
by the set S. To see this, let U be an orthonormal basis for the columns in A spanned by the index
set S, so that P = LUUT LT . Then LT PL = UUT = Q. Hence, using that kA − QAkF ≤ kAkF =
√
O( n) and squaring (24), we have,
kL̃A − PL̃Ak2F ≥ kA − QAk2F − O(1/n2 ).

(25)

Combining (22), (23), and (25),
kA − QAk2F ≤ (1 + ǫ/3)kA − Ak k2F + O(1/n2 ).

(26)

Note that kA − Ak kF ≥ 1 for any k < n, since if we remove the top row of A, creating a matrix B,
then kB − Bk kF ≤ kA − Ak kF , yet B is the identity matrix, and so kB − Bk kF ≥ 1. It follows that
the additive O(1/n2 ) error in (26) translates into a relative error, provided k < n and for 1/ε < n/3.
Hence, kA − QAk2F ≤ (1 + ǫ)kA − Ak k2F , and therefore by Fact 84, we have |S| > k/(2ǫ).
10.4.2

Column sharing of discretized orthonormal matrices

We need the following technical lemma about sampling random orthonormal matrices, concerning the number of columns they have in common after rounding.

64

Lemma 86. Let 1 ≤ r ≤ n, and let M ≤ 1/nc for a sufficiently large absolute constant c > 0. Suppose we
choose two independently random orthogonal matrices A ∈ Rn×n and B ∈ Rn×n from the Haar measure
(so they each have orthonormal columns and orthonormal rows). We then round each of the entries in A to
the nearest integer multiple of M , obtaining a matrix Ã, and similarly obtain B̃. Let E denote the event that
there exist r distinct indices i1 , . . . , ir ∈ [n] and r distinct indices j1 , . . . , jr ∈ [n] for which Ã∗,iℓ = B̃∗,jℓ
for all ℓ = 1, 2, . . . , r, where Ã∗,iℓ denotes the iℓ -th column of Ã. Then, Pr[E] ≤ e−Θ(nr log M ) .
Proof. We fix a subset S of r distinct indices i1 , . . . , ir ∈ [n] and a subset T of r distinct indices
j1 , . . . , jr ∈ [n] and show that the probability Ã∗,iℓ = B̃∗,jℓ simultaneously for all ℓ ∈ [r] is at most
e−Θ(nr log M ) . It then follows by a union bound that
   
 ne 2r
n
n
−Θ(nr log M )
e−Θ(nr log M ) ≤ e−Θ(nr log M ) ,
Pr[E] ≤
·
·e
≤
r
r
r
as desired.
Since A and B are independent and their distribution is permutation-invariant, we can assume
iℓ = jℓ = ℓ for ℓ ∈ [r]. Let F be the event that Ã∗,iℓ = B̃∗,jℓ occurs for all ℓ ∈ [r]. If F occurs, then
√
we need kA∗,iℓ − B∗,jℓ k∞ ≤ M for all ℓ ∈ [r], which implies kA∗,iℓ − B∗,jℓ k2 ≤ nM for all ℓ ∈ [r].
Let Ar be the leftmost r columns of A, and Br the leftmost r columns of B. Then, if F occurs
r

kA −

Br k22

r

≤ kA −

Br k2F

=

r
X
ℓ=1

kA∗,ℓ − B∗,ℓ k22 ≤ rnM 2 .

Note that Ar and Br have orthonormal columns, and so PAr = (Ar )(Ar )T and PBr = (Br )(Br )T
are projection matrices. Then,
kAr − Br k2 ≥ k(Ar )(Ar )T − Br (Ar )T k2 /2 + k(Ar )(Br )T − (Br )(Br )T k2 /2

= k(Ar )(Ar )T − (Ar )(Br )T k2 /2 + k(Ar )(Br )T − (Br )(Br )T k2 /2

≥ k(Ar )(Ar )T − (Br )(Br )T k2 /2,

(27)

where the first inequality follows since multiplying by a matrix on the right with orthonormal
rows cannot increase the spectral norm, the equality follows by taking transposes, and the second
inequality follows by the triangle inequality.
We need another result from [38].
Theorem 87. (Probability Bound - Claim 5.2 of [38]) Suppose we choose two independent random subspaces Y and S of Rn each of dimension r from the Haar measure, where we denote the projection operators
onto the subspaces by YY T and SST , respectively, where Y and S have r orthonormal columns. Then for
any δ ∈ (0, 1), Pr[kYY T − SST k2 ≤ δ] ≤ e−Θ(r(n−r) log(1/δ)) .
Combining (27) with Theorem 87, we have that for r ≤ n/2, Pr[F] ≤ e−Θ(rn log(M )) , since by
√
taking M ≤ 1/nc sufficiently small, we can ensure kAr − Br k2 ≤ rnM ≤ M 1/2 , and so 2M 1/2
upper bounds k(Ar )(Ar )T − (Br )(Br )T k2 .
Note that if r > n/2, then in particular we still have that Ã∗,iℓ = B̃∗,jℓ for all ℓ ∈ [r/2], and
therefore in this case we also have Pr[F] ≤ e−Θ(rn log M ) . Hence, Pr[E] ≤ e−Θ(rn log M ) , which
completes the proof.
We now extend Lemma 86 to a k-fold version.

65

Lemma 88. Let M ≤ 1/nc for a sufficiently large absolute constant c > 0. Suppose we independently
choose k pairs of matrices Az , Bz ∈ Rn×n , z ∈ [k], each with orthonormal columns (and hence also
orthonormal rows). We then round each of the entries in each Az and Bz to the nearest integer multiple
of M , obtaining matrices Ãz , and B̃z . Let r = (r1 , . . . , rk ) ∈ (Z≥0 )k and let E be the event that for each
z ∈ [k], there exist distinct indices i1 , . . . , irz ∈ [n] and j1 , . . . , jrz ∈ [n] for which Ãz ∗,iℓ = B̃z ∗,jℓ for all
Pk
ℓ = 1, 2, . . . , rz . Then, Pr[E] ≤ e−Θ(n z=1 rz log M ) .
Proof. This follows by Lemma 86, and the fact that the matrices A1 , . . . , Ak , B1 , . . . , Bk are jointly
independent. Here we also use a union bound.
Corollary 89. Suppose 1 ≤ k < c0 n for an absolute constant c0 > 0. Let r = (r1 , . . . , rk ) ∈ (Z≥0 )k . Let
M ≤ 1/nc for a sufficiently large constant c > 0. There exists a set N of eΘ(nt log M ) k-tuples (Ã1 , . . . , Ãk )
of matrices each of dimensions n × n with entries that are integer multiples of M and bounded in magnitude
P
by 1, such that for every vector r = (r1 , . . . , rk ) ∈ (Z≥0 )k with kz=1 rz ≤ t, it holds that for all distinct ktuples (Ã1 , . . . , Ãk ), (B̃1 , . . . , B̃k ), there exists a z ∈ [k] for which there are no sets S = {i1 , . . . , irz }, T =
{j1 , . . . , jz } ⊂ [n] for which Ã∗,iℓ = B̃∗,jℓ for all ℓ = 1, 2, . . . , rz .
Proof. For a fixed vector r and a random choice of (A˜1 , . . . , A˜k ), (B˜1 , . . . , B˜k ), by the guarantee of
Lemma 88, the complement of the event in the statement of the corollary happens with probability
P
at most e−Θ(nt log M ) . The number of vectors r with kz=1 rz ≤ t is at most tk , and tk e−Θ(nt log M ) ≤
e−Θ(nt log M ) provided k < c0 n for an absolute constant c0 > 0, so we can apply a union bound
over all pairs in a random choice of eΘ(nt log M ) k-tuples.
10.4.3

Hard instance construction

We have s servers holding matrices A1 , . . . , A
Ps , respectively, where Ai has m rows and a subset
of wi columns of an m × n matrix A, where wi = n. For our lower bound we assume 2/ε ≤ φ
and kφ ≤ m. Note that for ε less than a sufficiently small constant, this implies that k < c0 kφ, and
so the assumption in Corollary 89 is valid (since the n of that corollary is equal to (k + 1)φ We set
M = 1/(mn)c for a sufficiently large constant c > 0.
We will have A1 being an m × k(φ − 1) matrix for which all but the first kφ rows are zero. We
assume 2/ε ≤ φ. On the first kφ rows, A1 is block diagonal containing k blocks, each block being
a φ × (φ − 1) matrix.
To specify A1 , let the parameters t and n of Corollary 89 equal k/(2ε) and φ, respectively. That
corollary gives us a net N of eΘ(kφ(log(mn))/ε) k-tuples of matrices each of dimensions φ × φ with
entries that are integer multiple of M and bounded in magnitude by 1.
For a random k-tuple (B1 , . . . , Bk ) in N , we create an m × m block diagonal matrix B, which is
0 on all but its first kφ rows and first kφ columns. On its first kφ rows, it is a block diagonal matrix
B whose blocks along the diagonal are B1 , . . . , Bk , respectively. Our matrix A1 is then the matrix
product of B and an matrix D, where D is m × k(φ − 1) is formed by taking the kφ × k(φ − 1)
matrix of Fact 84 and padding it with m − kφ rows which are all zeros. By Lemma 85, no k/(2ǫ)
columns of A1 contain a k-dimensional subspace in their span which is a (1 + ǫ)-approximation to
the best rank-k approximationt to A1 .
Each Ai for i = 2, 3, . . . , s − 1 is the m × m matrix of all zeros. Finally, As is the m × t matrix
 of
all zeros with t = n − (s − 1)m − k(φ − 1), i.e., A = A1 0m×m 0m×m . . . 0m×m 0m×t . By
construction, each column of A has at most φ non-zero entries, since the only non-zero columns
are those in A1 , and each column in A1 has the form Bz∗,1 + Bz∗,i for a z ∈ [k] and an i ∈ [φ], so it
has at most φ non-zero entries.

66

10.4.4

Main theorem

In the Distributed Column Subset Selection Problem in Definition 18, a correct protocol should
have each of the s machines simultaneously output a matrix C ∈ Rm×c with c < n columns of A
such that
kA − CC† Ak2F ≤ (1 + ε)kA − Ak k2F .
Theorem 90. Assume 2/ε ≤ φ and kφ ≤ min(m, n). Then, any, possibly randomized, protocol Π
which succeeds with probability at least 2/3 in solving the Distributed Column Subset Selection Problem
of Definition 18, under the promise that each column of A has at most φ non-zero entries which are integer
multiples of 1/(mn)c and bounded in magnitude by 1, for a constant c > 0, requires Ω(sφk(log(mn))/ε)
bits of communication. If the word size is Θ(log(mn)), this implies the communication cost of the porticol
is Ω(sφk/ε) words.
Proof. By construction of A1 , whenever Π succeeds, each of the s machines outputs the same set
C of at least t ≥ k/(2ǫ) columns of A1 .
Let Πi denote the transcript (sequence of all messages) between the coordinator and the i-th
player. We look at the information Πi reveals about A1 , given the event E that Π succeeds, that is,
I(Πi ; A1 | E).
A critical observation is that any column C∗,i of C, is equal to Bz∗,1 + Bz∗,i for some z ∈ [k],
where Bz is as defined above (i.e., A1 = B · D). Hence, given Bz∗,1 and C∗,i , one can reconstruct
Bz∗,i . Now, for any random variable W ,
I(Πi ; A1 | E) ≥ I(Πi ; A1 | E, W ) − H(W ).

(28)

We let W = (B1∗,1 , . . . , Bk∗,1 ), and observe that H(W ) ≤ log2 ((2M )kφ ) since there are 2M choices
for each coordinate of each Bz∗,1 . Hence, H(W ) = O(kφ log(mn)). Note that I(Πi ; A1 | E, W ) =
H(A1 ) = log2 (|N |), since conditioned on the protocol succeeding, and given W , by Corollary 89
we can reconstruct A1 . By (28), it follows that
(Πi ; A1 | E) ≥ Θ(φk(log(mn))/ǫ) − O(kφ log(mn)) = Ω(φk(log(mn))/ǫ).
Note that if Z is an indicator random variable that indicates that E occurs, we have H(Z) ≤ 1, and
so
I(Πi ; A1 ) ≥ I(Πi ; A1 | E) Pr[E] − 1 = Ω(φk(log(mn))/ǫ),

since Pr[E] ≥ 2/3.
It follows that

Ω(φk(log(mn))/ǫ) ≤ I(Πi ; A1 ) ≤ H(Πi ) ≤ EA1 ,rand [|Πi |],
where |Πi | denotes the length, in bits, of the sequence Πi of messages, and rand is the concatenation of the private randomness of all s machines. Note that the entropy of Πi is a lower bound on
the expected encoding length of |Πi |, by the Shannon coding theorem. By linearity of expectation,
X
EA1 ,rand |Πi | = Ω(φk(log(mn))/ǫ),
i

which implies there exists an A1 and setting of rand for which
X
|Πi | = Ω(sφk(log(mn))/ǫ).
i

It follows that the total communication is Ω(sφk(log(mn))/ǫ) bits.
67

Corollary 91. Assume 2/ε ≤ φ and kφ ≤ min(m, n). Then, any, possibly randomized, protocol Π
which succeeds with probability at least 2/3 in solving the Distributed Column Subset Selection Problem
- rank k subspace version (see Definition 19), promise that each column of A has at most φ non-zero
entries which are integer multiples of 1/(mn)c and bounded in magnitude by 1, for a constant c > 0,
requires Ω(sφk(log(mn))/ε) bits of communication. If the word size is Θ(log(mn)), this implies the total
communication across all s machines is Ω(sφk/ε) words.
Proof. As any such protocol is also a protocol for the Distributed Column Subset Selection Problem of Definition 18 (the problem in Definition 19 is more general/difficult than the problem of
Definition 18), the proof follows immediately from Theorem 90.

68

