An Efficient Protocol for Distributed Column Subset Selection in the Entrywise $\ell_p$ Norm
Authors: Shuli Jiang
Dongyu Li
Irene Mengze Li
Arvind V. Mahankali
David Woodruff
Venue:  (2021)
Keywords: Column Subset Selection
Distributed Learning

Abstract: We give a distributed protocol with nearly-optimal communication and number of rounds for Column Subset Selection with respect to the entrywise {$\ell_1$} norm ($k$-CSS$_1$), and more generally, for the $\ell_p$-norm with $1 \leq p < 2$. We study matrix factorization in $\ell_1$-norm loss, rather than the more standard Frobenius norm loss, because the $\ell_1$ norm is more robust to noise, which is observed to lead to improved performance in a wide range of computer vision and robotics problems.
In the distributed setting, we consider $s$ servers in the standard coordinator model of communication, where the columns of the input matrix $A \in \mathbb{R}^{d \times n}$ ($n \gg d$) are distributed across the $s$ servers. We give a protocol in this model with $\widetilde{O}(sdk)$ communication, $1$ round, and polynomial running time, and which achieves a multiplicative $k^{\frac{1}{p} - \frac{1}{2}}\poly(\log nd)$-approximation to the best possible column subset. A key ingredient in our proof is the reduction to the $\ell_{p,2}$-norm, which corresponds to the $p$-norm of the vector of Euclidean norms of each of the columns of $A$. This enables us to use strong coreset constructions for Euclidean norms, which previously had not been used in this context. This naturally also allows us to implement our algorithm in the popular streaming model of computation. We further propose a greedy algorithm for selecting columns, which can be used by the coordinator, and show the first provable guarantees for a greedy algorithm for the $\ell_{1,2}$ norm. Finally, we implement our protocol and give significant practical advantages on real-world data analysis tasks.

Under review as a conference paper at ICLR 2021
AN EFFICIENT PROTOCOL FOR DISTRIBUTED COLUMN SUBSET SELECTION IN THE ENTRYWISE lp NORM
Anonymous authors Paper under double-blind review
ABSTRACT
We give distributed protocol with nearly-optimal communication and number of rounds for Column Subset Selection with respect the entrywise & norm (k-CSS1), and more generally, for the no-norm with 1 S p < 2. We study matrix factorization f--morm loss, rather than the more standard Frobenius norm loss, because the l norm is more robust to noise, which is observed to lead to improved performance in a wide range of computer vision and roboties problems. In the distributed setting, we consider servers in the standard coordinator model of communication, where the columns the input matrix A E Rdxx (n > d) are distributed across the & servers. We give a protocol in this model with (ssdk) communication, round, and polynomial running time, and which achieves a multiplicative ; poly((og nd) aaprroxixinnton to the best possible column subset. A key ingredient in proof is the reduction to the -ormm which corresponds to the p-norm of the vector of Euclidean norms of each of the columns of A. This enables us use strong coreset constructions for Fuclidean norms, which previously had not been used in this context. This naturally also allows us to implement our algorithm in the popular streaming model of computation. We further propose greedy algorithm for selecting columns, which can be used by the coordinator, and show the first provable guarantees for greedy algorithm for the f1 norm. Finally, we implement our protocol and give significant practical advantages on real-world data analysis tasks.
INTRODUCTION
Column Subset Selection (k-CSS) is widely studied approach for rank-k approximation and feature selection. In k-CSS, one seeks a small subset U E iddk of k columns of a data matrix A E dx typically n >> d, for which there is right factor V such that /U/V A| small under some norm 1 k-CSS is special case of low rank approximation for which the left factor is an actual subset of columns. The main advantage of k-CSS over general low rank approximation is that the resulting factorization is more interpretable, as columns correspond to actual features while general low rank approximation takes linear combinations of such features. In addition, k-CSS preserves the sparsity of the data matrix A.
Boutsidis let al. 2014. Boutsidis & Woodrufi 2017 Boutsidis et al.; 2008) and operator norms (Halko et al.
k-CSS has been extensively studied in the Frobenius norm Gurruswmmi & Sinop 2012
Woodruff 2014). A number of recent works (Song et al. 2017] CCherichetti et al.) 2017 Dan 2019 Ban al. (219)] Mahankali & Woodruff 2020) studied this problem in the Ep norm (K-CSS,,) for I < 2. The & norm is less sensitive to outliers, and better at handling missing data and non-Gaussian noise, than the Frobenius norm (Song et al.) (0017). Specifically, the l norm leads to improved performance in many woal-world applications, such as structure- -fom--rooo-ooo------ motion -roo--moionn (Ke) Kanade 2005) and image denoising (Yu al., 2012).
2011
al.
Distributed low-rank approximation arises naturally when a dataset is too large to store on one machine, takes prohibitively long time for a single machine to compute a rank-k approximation, or is collected simultaneously multiple machines. Despite the flurry of recent work on k-CSS,, this problem remains largely unexplored in the distributed setting. This should be contrasted to Frobenius norm column subset selection and low rank approximation, for which number of results in the distributed model are known, see, e.g.. Altschuleeuler et al.) (2016); |Blcaan et al. (2015) 2016); Boutsidis et al. (2016). We consider widely applicable model in the distributed setting, where

---

Under review as a conference paper at ICLR 2021
Sep4
66-T
A/k
Ar
Shit, ShoT
SA.7, A0,
40.
GAR
AA
Servers
Coordinator
An overview of the proposed protocol for distributed k-CSS, in the column partition model. Step 1: Server applies a dense p-stable sketching matrix S to reduce the row dimension of the data matrix A, S is shared between all servers. Step 2: Server i constructs a strong coresef for its sketched data matrix A, which a subsampled and reweighted set of columms of SA,, Server i then sends the coreset SA,,; as well as the corresponding unsketched, unweighted columns A;D; selected in the strong coreset S SA,T; to the coordinator. Step 3: The coordinator concatenates the SA,T, column-wist -wiee applies k-CSSpu 2 to the concatenated columns and computes the set of indices of the selected columns. Step 4: The coordinator recovers the set of selected columns Ay from the unsketched, unweighted columns D,'s through previously computed indices.
Figure 1:
servers communicate central coordinator via 2-way channels. This model can simulate arbitrary point-to-point communication by having the coordinator forward a message from one server to another; this increases the total communication by factor of 2 and additive log s bits per message to identify the destination server.
We consider the column partition model, in which each column of A E R is held by exactly one server. The column partition model widely-studied studied and arises naturally in many real world scenarios such as federated learning (Farahat et al. 2013 Altschuler et al. 2016 Liang et al., 40.4). In the column partition model, we typically have >> d, i.e., has many more columns than rows. Hence, we desire protocol for distributed k-CSS, that has a communication cost that is only logarithmic in the large dimension n, as well as fast running time. addition, it is important that protocol only uses small constant number of communication rounds (meaning buck-and-fortt exchanges between servers and the coordinator). Indeed, otherwise, the servers and coordinator would need to interact more, making the protocol sensitive failures in the machines, e.g., if they go offline. Further, a --ronnd protocol can naturally be adapted to an single pass streaming algorithm when we consider applications with limited memory and access to the data. In fact, our protocol can be easily extended to yield such a streaming algorithm
In the following, we denote A,a and A.j as the i-th row and j-th column of A respectively. for E [4].j E [n]- We denote Ar as the subset of columns of with indices in T S [n)]. The entrywise l, no-mo of A is |Alp = E,
The lp. norm is defined as |Alp.2 = (E ( We consider 1 S p < 2. We denote the best rank-k approximation error for A in 9p norm by OPT min A2 |A Axlp- Given an integer k > we say U d V E R are the left and right factors of a rank-k factorization for A in the lp norm with approximation factor JUV - Alp S a OPT.
Since general rank-k approximation in 4 norm is NP hard Gillis & Vavasis 2015), we follow previous work and consider bi-criteria k-CSS algorithms which obtain polynomial running time. Instead of outputting exactly l columns, such algorithms return a subset of Ō columns of A, suppressing logarithmic factors in le or n. It is known that the best approximation factor to OPT that can be obtained through the span of column subset of size Ö(k) is 1/21/227) for p 1 (Song et al. 2017) and 1(k21/p112211 for p E (1 (Mahankali & Woodruff 2020), where is an arbitrarily small constant.
We provide detailed analysis of our streaming algorithm in Appendix
2

---

Under review as a conference paper at ICLR 2021
1.1 PREVIOUS APPROACHES TO k- k- IN THE DISTRIBUTED SETTING
If one only wants to obtain good left factor U, and not necessarily column subset of A, in the column partition model, one could simply sketch the columns of A by applying an oblivious sketching matrix S on each server. Each server sends A; SS to the coordinator. The coordinator obtains U = -AS as a column-wise concatenation of A,5's. [Song et al. (2017) showed that AS achieves an ((vk) approximation to OPT, and this protocol only requires ((sdk) communication, O(1) rounds and polynomial running time. However, while AS is good left factor, it does not correspond to an actual subset of columns of A.
Obtaining a subset of columns that approximates A well with respect to the p-norm in a distributed setting is non-trivial. One approach due to Song et al. (2017) is to take the matrix AS described above, sample rows according the Lewis weights (Cohen & Peng, 2015 of AS get a right factor V. which is in the row span of A, and then use the Lewis weights of V to in turn sample columns of A. Unfortunately, this protocol only achieves a loose ((23/2) approximation to OPT (Song ctt al., 2017). Moreover, it is not known how to do Lewis weight sampling in a distributed setting. Alternatively, one could adapt existing single-machine k-CSS, algorithms to the distributed setting under the column partition model. Existing works on polynomial time k- k-CSSp (Chierichetti et al. 2017:(Song et al., [019bb] Dan al. | Mahankali & Woodruff] (220)) give bi-criteria algorithms, and are based on a recursive framework with multiple rounds, which is as follows: in each round, Ō(k) columns are selected uniformly at random, and with high probability, the selected columns can provide a good approximation to constant fraction of all columns of A. Among the remaining columns that are not well approximated, (kk) columns are recursively selected until all columns of A are well approximated. resulting in a total of Olog n) rounds.
A naive extension of this bi-criteria k-CSS, framework to distributed protocol requires Ollog n) rounds, as in each round, the servers and the coordinator need to communicate with each other in order to find the columns that are covered well and select from the remaining unselected columns. To reduce this to a single round, one might consider running the Olog n) round selection procedure on the coordinator only. In order to do this, the coordinator needs first collect all columns of A from the servers, but directly communicating all columns is prohibitive.
Alternatively, one could first apply A-CSS, on A; to obtain factors U2 and Vi each server, and then send the coordinator all of the U, and V. The coordinator then column-wise stacks the U.V; to obtain U- V and selects (kk columns from U V. Even though this protocol applies to all p 2 1, achieves a loose O(k²) approximation to OPT and requires prohibitive O(n + d) communication cos One could instead try to just communicate the matrices U, to the coordinator, which results in much less communication, but this no longer gives a good approximation. Indeed, while each U; serves as a good approximation locally, there may be columns that are locally not important, but become globally important when all of the matrices A, are put together. What really needed here a small coreset C; for each A; so that if one concatenates all of the C; to obtain C. any good column subset of the coreset C' corresponds to good column subset for A. Unfortunately, coresets for the entrywise -n-norm are not known to exist.
fl
1.2 OUR CONTRIBUTIONS
Our Distributed Protocol We overcome these problems and propose the first efficient protocol for distributed k-CSSp < in the column partition model that selects Ō(k) columns of A achiev- to the best possible subset of columns and requires only Ö(sdk)
ing an (1tt://=1/))aapprxximii) communication cost, 1 round and polynomial time. Figure gives overview of the protocol. We note that our subset of columns does not necessarily achieve an itself, although it does achieve such an approximation to the best possible subset of columns. Using to OPT
to OPT
the fact that there always exists subset of columns providing an http:/--t//)-appoxxxmie (Song et al. 2017), we conclude that our subset of columns achieves an ( - OPT. Recently, and independently of our work, Mahankali & Woodruff (2020) show how to obtain a subset of columns achieving an http://=/t"aapppppximtttto to OPT itself; however, such a subset found by uniformly sampling columns in Ollog n.) adaptive rounds using the recursive sampling
)-approximation to
We give this protocol and the analysis in Appendix
3

---

We consider a model where there are servers, the hth of which holds A, E Rdxm; and coordinator which initially does hold any data. Each server talks only the coordinator, via 2-way communication channel. The communication cost is the total number of words transferred between the servers and the coordinator over the course of the protocol. Each word (log(snd) bits. The Under review as a conference paper at ICLR 2021
framework above, and is inherently hard implement in a distributed setting with fewer rounds. In contrast, our protocol achieves round of communication, which is optimal.
We make of strong coreset, i.e., sampled and reweighted subset of columns of each A, that approximates the cost of all potential left factors of A, by first embedding all subspaces spanned by any subset of Ōkk) columns of A from Space to Euclidean space, to bypass the lack of strong coresets for the lp norm. We denote this new norm as l..2 norm, which is the sum of the p-th powers of the & norms of the columns. To reduce the error incurred by switching to the fp.2-norm, the row dimension of A by multiilli by an oblivious sketching matrix shared across servers, resulting in an overall approximation factor of only O(k2/p-1/2) Afterwards, each server sends its strong coreset the coordinator. The coordinator, upon receiving the coresets from each server, runs an (1)-approximmae bi-criteria k- k-CSSp2 p.2 algorithm to select the final column subset, giving an overall (1,1/P-11 12/)) approximation the best column subset.
reduce
We introduce several new technical ideas in the analysis of our protocol. Our work is the first to apply a combination of oblivious sketching in the p-norm via p stable random variables and strong coresets in the fp.2 norm (Sohler & Woodruff [0088] Huang & Vishnoi (2200) to distributed k-CSS. Furthermore, to show that our oblivious sketching step only increases the final approximation error by logarithmic factor, we combine net argument with a union bound over all possible subspaces spanned by column subsets of A of size O(k). Previous arguments involving sketching, such as those by Song et al. (2017); Ban et al. 2019); Mahankali & Woodruff (020),, only consider a single subspace at time.
Theoretical Guarantees and Empirical Benefits for Greedy k-CSS1,2 We also propose greedy algorithm to select columns in the k-CSS1.2 step of our protocol, and show the first additive error guarantee compared to the best possible subset As of columns, i.e., our cost at most (1 <) miny | Al1,2 + A|A/,, Similar error guarantees were known for the Frobenius norm (Altschuler et al.) 2016), though nothing was known for the 4,2 norm. We also implement our protocol and experiment with distributed k-CSSS, various real-world datasets. We compare the appppxppoe bi-criteria 1-CSS1,2 and the greedy k-CSS1,2 as different possible subroutines in our protocol, and show that greedy k-CSS1,2 yields an improvement in practice.
2
PROBLEM SETTING
2.1 THE COLUMN PARTITION MODEL
overall data matrix A Id dxdx is A = [Ai, A A.] (the column-wiss concatenation of the A('s). Here, n is defined be Ei ny, Typically, in the column partition model, n d.
2.2
P-STABLE DISTRIBUTION AND PTABEE RANDOM VARIABLES
Let Z, XI....... Xa be random variables drawn i..dd from some distribution D. D is called p-stable for an arbitrary vector u E Rd, (u, = |ul,Z, Z. where X = [X,,..,,X,, X.)". Random variables drawn from such distribution are called p-stable random variables, which exist for p E (0, 2]. Though there is no closed form expression for the p-stable distribution in general except for a few values of p. we can efficiently generate a single p-stable random variable in O(1) time using the following method due to Chambers et al. 1976): if 0 E 2 2 and r e [0, 1] are sampled uniformly follows a p-stable distribution.
at random, then, cost/Fg
sin(pe) (0s(
3 PRELIMINARIES FOR OUR PROTOCOL
We first note a standard relationship between the lp norm and the fp. norm.
Lemma 1. For a matrix A E Bddx and p (1, |Alm2 S IAlp S di |App22
4

---

Under review as a conference paper at ICLR 2021
3.1 f--NORM OBLIVIOUS SKETCHING
We left-multiply A by an oblivious sketching matrix S with p-stable random variables so that we lose only an Ö(k* 2) approximation factor when we switch to the f,.2 norm.
The purpose of the next two lemmas is show that can perform oblivious sketching while preserving the costs of all possible column subsets up to logarithmic factors. We first show a lower bound on the approximation error for a sketched subset of columns, AS - -SAlp which holds simultaneously for any arbitrary subset AT of chosen columns, and for any arbitrary right factor V.
Lemma 2 (Sketched Error Lower Bound). Let E R and & E N. Lett = k poly(log(dd)), d be matrix whose entries are i.i.d. standard p-stable random variables. rescaled by (1 (1/th) ). Then, with probability
and let S E B
(. for all T C [n] with [T] = k poly(log k) and for all SAArV SAlp-
V E R Tixm,
|Ar
S Alp
A
Next, we show an upper bound the approximation error of k-CSS, sketched subset of columns, |SATV- - SAlp. which holds for a fixed subset of columns Ar and a fixed right factor V.
Lemma 3 (Sketched Error Upper Bound). Let A R and is E N. Let k poly(log(nd)) and let S E Rxxd be a matrix whose entries are i...d standard p-stable random variables, rescaled by ((1/t²). Then, for a fixed subset T C [n]] of columns with [T] = k pooyolooyooo k) and fixed V e RITix,, with probability 1 - o(1), we have miny ASArV SAlp miny O(log' (Ilg)/M(nd)))) ArV Alp-
3.2 STRONG CORESETS THE lp.2 NORM
To enable sub-linear communication cost in the number of columns, the i-th server sends the coordinator a strong coreset of columns of SA, which is a reweighted subset of the columns of SA,- Such strong coresets preserve the error incurred by any rank-k projection, up to constant factor, in the fa.2 norm. The coreset of a matrix A E [ddx is usually denoted as AT', where D- W Rnxt is a sampling and reweighting matrix and is the number of columns to be included in the coreset. The sampling matrix D is matrix with t columns where each column has only one in the index of the column of A be included in the coreset and everywhere else. The reweighting matrix W is diagonal ( i matrix with weights associated with each sample in the coreset.
Lemma (Strong Coreset in f2.2 norm). Let A E Rdxn, k E E. e (1, and e, E (0,1). Then, in n poly(k log n/e) time, one can find sampling and reweighting matrix T with Od d/e2). log(1/0 /6) columns such that, with probability for all rank-k matrices U,
min IUV - -AT p.2 (1 1l nonko V
min ronk-k V
IUV - -Alp.2
AT is called strong coreset of A.
3.3 POLYNOMIAL TIME, O(L)-APPROXIMAT BI-CRITERIA k-CSSp.2
After server i sends strong coreset to the coordinator, the coordinator does k-CSS on column-wi woise concatenation of these coresets, the fp.2 norm rather than the lp norm. We give polynomial time, ((1))approximate bi-criteria k- k-CSSp.2 algorithm for p E (1.
Theorem 5 (Bicriteria 0(1) O(1)-Approximation Algorithm for k-CSSp.2). Let E Rxdx and i N. There exists algorithm that runs in (mnz(A) dP). kpoly(log k) time and outputs rescaled subset of columns U E gdx((k) of A and a factor V e (o(k))) for which V = miny IUV - AApp2 such that with probability 1 o(1)
IUV Alp2 S 0(1) min
|Ag
Alp.2
rand-k Ac
Our polynomial time bi-criteria k-CSSp.2 algorithm is based on that of Clarkson & Woodruff (2015).
The main difference that the algorithm of Clarkson Woodruff (2015) outputs a subset with (i²))
columns due to the usage of &, leverage scores we reduce the number of selected columns to Ö(k)
by using t Lewis weights (Cohen & Peng. 2015). Details are given in Appendix C
5

---

Under review as a conference paper at ICLR 2021
Algorithm I An efficient protocol for bi-criteria k-CSS, in the column partition model
Initial State: Server holds matrix A, Rdxxs Vi E [s]. Coordinator:
Generate dense p-stable sketching matrix S E R pgly Send 5 to all servers.
Server
Compute SA,,
Let the number of samples in the coreset be (Igpol)))o(nn)), Construct coreset of SA, under the fp.2 norm by applying a sampling matrix D; of size n, X t and diagonal reweighting matrix W; of size t. Let T; = W. Send SA,T; along with A,D; to the coordinator.
Coordinator:
Column-wise stack SA;7; to obtain SAT = [SA, ,
Apply k-CSSp, on SAT to obtain the indices I the subset of selected columns with size O(k poly(log k)).
Since Dj's are sampling matrices, the coordinator can recover the original columns of A by mapping indices I to A,D;'s
Denote the final selected subset of columns by As. Send A to all servers.
Server
Solve miny, |A,V, Ailp to obtain the right factor V₁ A₁ and will be factors of a rank- poly(log k) factorization of A, where is the (implicit) column-wist concatenation of the V,.
4 AN EFFICIENT PROTOCOL FOR DISTRIBUTED k CSC
Theorem (A Protocol for Distributed k-CSSp)) In the column partition model, let A E Rd be the data matrix whose columns are partitioned across servers and suppose server i holds subset of columns A₁ RRdxns where Eiela 0g- Then, given p E (1,2) and desired rank N, Algorithm ooupptt a subset of columns A₁ E (Rdxkpol((oo(()) in ((nz A)k + kd + k.3) time, such that with probability - o(1),
min LA,V Alp - 0(K//111/2)
min LCM
ALL Alp
Algorithm 1 Juses! round of communication and Ö(skk) words of communication.
Proof. Approximation Factor. In the following proof, let L C [m], = k denote the best possible subset of & columns A that gives the minimum k-CSS, cost, i.c. the cost miny |Ag V Alp achieves minimum. First, note that
min A,V - Alp V
< AA -App
V" = arg |SA, - -SA p.2
< VSA,V' - -SApp
By Lemma 2
= (IAA;-))SAA/' -SApp,
By Lemma I and has poly(log(nd)) rows
SAr the selected columns output from the bi-criteria approxiaiitnttoonnn k-CSS,.2 algorithm.
Let (SA denote the best rank l: approximation to SAT. By Theorem 5
Ö(k SAp,, s Ō(k+-T)- ((A) -SAT p.2
< )
-SAT p.2
Note that SAT [SAT Vi [s). By Lemma
SA.T. is a column-wise concatenation of all coresets of S.,,
(min SAAV
SAT 1/p p.2)
1/p
min SALV - SA,T)), 1/P V.

---

Under review as a conference paper at ICLR 2021
+ min [SA.V SAA.p.2) -1+ (1+
min |SALV SA..p.2)
Vi
+ min [SAVV - -SA
Hence,
0(K min SALV SAT 1 Õ(k -)) min |SALV -SAlp.2 V
V
< Ōkk -b) min |SA,V SAlp
By Lemma]
- Ō(kk- +) loog (nd) min AAL - Alp By Lemma
V
Therefore, conclude: miny |A, - Al, S Ō miny |Av AApp
Communication Cust. Sharing the dense p-stable sketching matrix S with all servers costs (sdk poly(log(nd)) communication (this can be removed with shared random seed). Sending all coresets S S,T,; (vi E [s]) and the corresponding columns A, to the coordinator costs (sdkk) communication, since each coreset contains only Ō(k) columns. Finally, the coordinator needs (sdk) words of communication to send the (kk) selected columns to each server. Therefore, the overall communication cost is (sdk)), suppressing logarithmic factor in n, d.
Running time. Since generating a single p-stable random variable takes O(1) time, generating the dense p-stable sketching matrix S takes O(dk poly(log(nd)) time. Computing all SA,'s takes (nnz(AA): poly(logy))) (nd)) time. By Lemma 4 computing a single coreset SA,T, on server i takes time n.poly(k log n.)) and thus computing all coresets across the servers takes at most time npoly(k log n.). By Theorem 5 the k-CSSp.2 algorithm takes time (nnz(SAT) + k² poly (log nd)) kpoly (log k) <. poloollllloooo (logy knd)) to find the set of selected columns. Since the number of selected columns O(pooly (log k)). then takes the protocol (kkoly (log k)) time to map the indices of the output columns from C-CSSp.2 to recover the original columns A,. Therefore, the overall running time for the protocol to find the subset of columns A, is (((nn((A)
kdd + k²)), suppressing low degree polynomial dependency on log log(knd). After the servers receive A,. it is possible to solve miny, A V/2 Aslp in ((nnz(Ar) + poly(d n) time, Vi E [8] due to Wang & Woodruff (2019) Yang et al. (2018).
5 GREEDY 4- 1-CSS1.2
We propose a greedy algorithm, shown in Algorithm 2 for k-CSS which can be used in the place of the algorithm described in Theorem The basic version of this algorithm, discussed in Appendix D performs k-CSS1,2 by simply selecting the additional column, among those of A, that reduces the approximation error the most at each iteration. Our analysis of that algorithm is inspired by the analysis of Greedy k-CSS, for the Frobenius norm in Altschuler et al. (2016). Here we provide the first additive error guarantee, compared to the best possible subset of columns, for the greedy k-CSS1, 2 algorithm. For a faster running time, we make use of the Lazier-than-lazy heuristic described in Section 5.2 of Altschuler al. (2016), where in each iteration, rather than considering all columns of A as candidate additional columns of Ar, we only sample subset of the columns of A of size O n log(1/0) ), and pick the column among those that improves the objective the most.
Theorem 7. Let A E Rdxın be the data matrix and ls E N be the desired rank. Let As be the best possible subset of l columns, i.e,, As = arg min Ag miny \AgV - Al1,2- Let a be the minimum non-zero singular value of the matrix of normalized columns of As, (the j-th column of B is B.j = (As)xş/ (As)))())))))))))((
(5>gl2) (2). Then, ifT C [n] is the subset of columns selected by Greedy k-CS51.2
the following holds with [T] ② g- }b=)):
min AT V A11,2 S c1
min ScIn
LASV -Al1.2 E|A/1.2
Similarly, T C [n]] is the subset of columns selected by Lazier-than-lazy Greedy k--SS1, the following holds with (T) 2( ( and =
Elmin Ar - A11,2] $ (1-c)
min AASV - -A/1,2 + c|Al1,2
7

---

Under review as a conference paper at ICLR 2021
Algorithm 2 Lazier-than-lazy Greedy k-CSS1,22 This version of the greedy algorithm is based on Section 5.2 of Altschulee al. (2016).
Input: The data matrix A E [2x= the number of iterations n, parameter ó (0, 1). Output: A subset of columns Ar from A, where [7] = r.
Ar 0
for i = 1 to do
T t- subset of (=log(1/6) columns of A, each selected uniformly random (excluding the
columns whose indices are in T)
Column j f arg min (miny |ATuv - 1
AT 4 ATUj"
end for
Datasct
Size
Coloma Distribation 1001. 1002
Rank [10, 20, 30)
synthetic TechTC
(2000 KJ (2000 K)) 139 18446
922 columns servers, 928 columns
{10, 30, 70, 100)
Table A summary of datasets used in the experiments.
If we let [2'] =
(}( then the overall running time of Algorithm 2 is O n log(1/e) a=e²
F), where
F is the running time needed to evaluate miny [Aruy V -Al1,2 for fixed j T. We can get
F = ( dk² ndk by taking Aluy-
ares +
Since the error upper bound for greedy k-CSSp.2 depends on |Al it not directly comparable to the error upper bound for the proposed k-CSSp.2 from Subsection 3.3 which achieves multiplicative (1) approximaationnnn to the best rank approximation. We empirically compare the two versions of k K-CSSp.2 2.2 2 for p in Section
EXPERIMENTS
We implement our protocol for distributed k-CSS,, in Algorithm D setting 1, which enables compare two subroutines the coordinator: Regular K-CSS1,2 from Section 3.3 and Greedy k-CSS1,2 from Section 5 We compare our k-CSS, protocol against a commonly applied baseline for Ep low rank approximation (used by Song et (2019)): Chierichetti et al. (2017)): rank-k Singular Value Decomposition (SVD).
Datasets. We demonstrate the benefits of our k-CSS) protocol on one synthetic data and one real-world application. We present summary of the datasets, along with the number of servers 8. the column distribution across servers and the rank l we consider for each dataset in Table []] The synthetic dataset constructs data matrix M E [(+++))(+++)) such that the top left k submatrix has all
submatrix the identity matrix multiplied by a and the bottom right I's. The optimal rank-k left factor consists of one of the last n columns along with k - 1 of the first l columns, incurring an eтror of nå in the & norm and an error n in the squared lz norm. SVD, however, will not cover any of the last n columns, and thus will get an error of n² in both the & and squared f2 norms. We set n 2000 and apply i.i.d. Gaussian noise to each entry with mean 0 and standard deviation 0.01. We consider real-world application, term-document clustering, where k-CSS, algorithm is previously applied (Mahoney & Drincas 2009). Te TechTd contains 139 documents processed in bag-of- words representation with a dictionary 18446 words. Such representation naturally results in sparse matrix. k-CSSp is used select the top k most representative words.
Hyperparameters. Experiment hyperparameters are summarized in Table We denote the number of rows in our 1-stable (Cauchy) sketching matrix by cauchy size, and the strong coreset size by coreset size. We have two additional hyperparameters for regular k--SSS, We denote the number of rows in the sparse embedding matrix of Ō(k) rows by sketch size, and the number of non-zero entries in each column the sparse embedding matrix by sparsity.
8

---

Under review as a conference paper at ICLR 2021
synthetle
TechTC
Greedy
Regular
Greedy
Regular
cauchy size
8k
40
corese: size
10
10k
min(4 x caucay size, 6k) min(4 X cauchy size, lik)
sketch size
1/3
k/3 mini5,
sparaity
min/5, k/3)
Table 2: A summarization of hyperparameters used for each dataset. Note that a too small coreset size compared to cauchy size will incur large l cost, and coreset size needs to be increased as rank k: increases.
Spsthete-Counder
TechfCC Seeuts
+000000
13:000
714044
140000
103000
20000
1000
00000
3000
43000
and SARCCAAAN
(a) synthetic
(b) TechCC
Results on synthetic and TechTC. The green line denotes Greedy k-CSS the orange lines denotes Regular k- k-CSS and the blue line denotes SVD.
Figure 2:
Setup. We conduct 15 trials for each experiment setting, and report the average & error, i.e. IÁg Al1 for a -(k)-rank matrix As output by the algorithm that approximates the original data matrix A E Rdx", across all trials. However, we notice that computing the & eтror for our protocol involves computing an & regression on all columns of the data matrix, which takes prohibitively long time the real-world dataset Tech TC. Therefore, for TechTC, we proceed with approximate computation as follows: for the output subset of columns A,, first i.i.d. sample 30% columns from A. denoted as A,, with probability proportional to |A, 1. Vj [n]: then compute lower bound for the & eгror L = ΣjE, min, |Are A.5l1: after that, sample 20 columns from A, denoted as Ag. with probability again proportional to /A.511. Vjj E [m]: finally, the ε error is computed as TAT Ekek K 20]A..]T min, |Are - A. The approximated & error can be shown as unbiased estimation of miny [A, - A/1 with small variance by Hoeffding's bound.
L.
The distributed protocol performs better using GREEDY k-CSS1,2 than REGULAR 1- both datasets, and in other settings include in the supplementary material.
Results. We present our empirical results in Figure 2
We also conducted experiments on three different datasets, bcsstk13, isolet. and caltech-101, extensively comparing the approximation error well as the running time for different algorithms. A comprehensive reporting of our results is given in Appendix G
CONCLUSION
In this work, we give the first nearly-optimal communication and number of rounds protocol for distributed k-CSS, (1 S < 2) in the column partition model, which achieves 0(k//-11/2) approximation to the best possible subset of columns, with ((sdk) communication cost, 1 rounds and polynomial time. To achieve a good approximation factor, we use dense p-stable sketching and work with the 4p.2 norm, which enables to an efficient construction of strong coresets and an ((1) apprxxxpooniiaaiion bi-criteria k-CSSp.2 p.2 algorithm. We further propose greedy algorithm for k- CSS1,2 and show the first additive error upper bound compared to the best possible subset of columns. We implement our distributed protocol using both greedy k- k-CSS1,2 and regular k-CSS1, Our results empirically show that greedy k-CSS 1 gives substantial improvements over regular J,CSS1,2 on real-world datasets. For future works, it is not known whether a Okk1/r 1/p////apprxximmatioo factor to the best possible subset of columns, is optimal for distributed k-CSS, (1 S p < 2).
9

---

Under review as a conference paper at ICLR 2021
REFERENCES
Jason Altschuler, Aditya Bhaskara, Gang Fu, Vahab Mirrokni, Afshin Rostamizadeh, and Morteza Zadimoghad- dam. Greedy column subset selection: New bounds and distributed algorithms. Proceedings of the 33rd International Conference on International Conference on Machine Learning Volume 48, ICML'16. pp. 2539-2548,, JMg..org, 2016.
FlriiaFlorina Balcan, Yingyu Liang, I.e Song, David P. Woodruff, and Bo Xie. Distrihuted kernel principal component analysis. CoRR, 103/1500306888, 2015.
Poria-Florina Balcan, Yingyu Liang. Le Song, David P. Woodruff, and Bo Xie. Communication efficient distributed kemee principal component analysis. In Balaji Krishnapuram, Mohak Shah, Alexander J. Smola, Charu C. Aggarwal, Dou Shen, and Rajeev Rastogi (eds.), Proceedings of the 22nd ACM SIGKDD Interna- tional Conference on Knowledge Discovery and Data Mining, San Francisco, CA, USA, August 13-17, 2016, pp. 725-734. ACM, 2016.
Frank Ban, Vijay Bhattiprolu, Karl Bringmann, Pavel Kolev, Euiwoong Lec, and David P. Woodruff. A PTAS for lowlow rank approximation. In Timothy M. Chan (cd.), Proceedings of the Thirtieth Anmual ACM- SIAM Symposium on Discrete Algorithuns, SODA 2019, San Diego, California, USA, January 6-9, 2019, pp. 747-766. SIAM, 2019. doi: 10.1137/1..7816111 1975482.47. htttttts:://do. org 1137/1
http://ooo.rrg/1/////
9788611975482.44
Christos Boutsidis and David Woodruff. Optimal cur matrix decompositions. SIAM Journal on Computing. 46 (2):543-589, 2017.
Christos Boutsidis, Michael Mahoney, and Petros Drineas. An improved approximation algorithm the column subset selection problem. CoRR, abs/0812.4293, 2008. BRLLhttp:/// httttp:://axxl..or////aaa/
4293
Christos Boutsidis, Petros Drineas, and Malik Magdon-Ismail. Near-optimal column-based matrix reconstruction. SIAM Journal on Computing, 43(2))687-717, 2014.
Christos Boutsidis, David P. Woodruff, and Peilin Zhong, Optimal principal component analysis in distributed and streaming models. In Daniel Wichs and Yishay Mansour (eds.). Proceedings of the 48th Annual ACM SIGACT Symposium on Theory of Computing, STOC 2016, Cambridge, MA, USA, June 18-21, 2016. pp. 236.249, ACM, 2016.
J. M. Chambers, C. L. Mallows, and B. W. Stuck. method for simulating stable random variables. Journal of the American Statistical Association, 71(354):3403344, 1976. doi: 10.1180001221459.11100 10480344. URL
https: 77 www. aannfoooonnn com/
Flavio Chierichetti, Screenivas Gollapudi, Ravi Kumar, Silvio Lattanzi, Rina Panigrahy, and David Woodruff. Algorithms for é, Low-Rank Approximation. 2017.
Kenneth L. Clarkson and David P. Woodruff. Input sparsity and hardness for robust subspace approximation. In Proceedings of the 2015 IEFE 56th Annual Symposium on Foundations of Computer Science (FOCS), FOCS '15. pp. 310-329, USA, 2015. IEFF Computer Society. ISBN 9781 467381918. doi: 10.1 109///0CS.0005..27
Michael B. Cohen and Richard Peng. Lp row sampling by lewis weights. In Proceedings of the Forty-Seventh Annual ACM Symposium Theory of Computing, STOC '15, pp. 183 192, New York, NY, USA, 2015. Association for Computing Machinery. ISBN 9781450335362. doi: 10.1 (41144776633992746556)
Michael B. Cohen, Sam Elder, Cameron Musco, Christopher Musco. and Madalina Persu. Dimensionality reduction for k-means clustering and low rank approximation. Rocco A. Servedio and Ronitt Rubinfeld (eds.), Proceedings of the Forty-Sevent Annual ACM on Symposium on Theory of Computing, STOC 2015, Portland, OR, USA, June 14-17, 2015, pp. 172. ACM, 2015. doi: 10. 10.114552746539.2746666 URL
https: : / /ccci. wrw./erg///// 1141 11w...111112225211111
2746569
Chen Dan, Hong Wang, Hongyang Zhang, Yuchen Zhou, and Pradeep K Ravikumar. Optimal analysis of subset- selection based Lp low-rank approximation. In H. Wallach. H. Larochelle. A. Beygelzimer. F. d'Alché-Buc. F. Fox, and R. Gareett (eds.), Advances Neural Information Processing Systems 32, pp. 2541-2552. Curran Associates, Inc., 2019.
Ahmed K. Farahat, Ahmed Elgohary, Ali Ghodsi, and Mohamed S. Kamel, Distributed column subset selection on mapreduce. In Hui Xiong. George Karypis, Bhavani M. Thuraisingham, Diane J. Cook, and Xindong Wu (eds.), 2013 IEEE 13th International Conference on Data Mining, Dallas, TX, USA, December 7- 10, 2013, pp. 171-180 IEEF Computer Society, 2013. doi: 10.1009/CDM222233.115 7/doi .rgg 10. 1109/ CDM, 2013. 155
_ httpp
10

---

Under review as a conference paper at ICLR 2021
Nicolas Gillis and Stephen A. Vavasis. On the complexity of robust PCA and é1-norm low-rank matrix approximation. CoRR, 1501.02996996 2015.
Venkatesan Guruswami and Ali Kemal Sinop. Optimal column- low-rank matrix reconstruction. In
Proceedings of the twenty-third annual ACM-SIAM symposium on Discrete Algorithms, pp. 1207-1214. SIAM, 2012.
Nathan Halko, Per-Gunnar Martinsson, and Joel A Tropp. Finding structure with randommess: Prohabilistic algorithms for constructing approximate matrix decompositions. SIAM review, 53(2):217-288 2011.
Lingxiao Huang and Nisheeth K. Vishnoi. Coresets for clustering in euclidean spaces: Importance sampling is nearly optimal. In Proceedings of the 52nd Annual ACM SIGACT Symposium Theory of Computing. STOC 2020, pp. 1416-1429, New York, NY. USA, 2020. Association for Computing Machinery. ISBN
9781450369794. doi: 10.1145/35771333384222 URL bttps://cole wtt.:oo....oogoooooo wttg:/..0000....31333333 B384296
Qifa Ke and Takeo Kanade. Robust 1 factorization in the presence of outliers and missing data by alternative convex programming. In 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR 2005), 20-26 June 2005, San Diego, CA, USA, pp. 739-743-97777 746. IEEE Computer Society, 2005. doi: 10. 10.1109CVPR.2000
URL URLnttps::// /do:://///o... org/ / 0. 1109/CVPR. htt.:///VR............ 309
Yingyu Liang Maria-Florina Balcan, Vandana Kanchanapally, and David P. Woodruff. Improved distributed principal component analysis. In Zoubin Ghahramani, Max Welling, Corinna Cortes, Neil D. Lawrence, and Kilian Q. Weinberger (eds.), Advances in Neural Information Processing Sys- tems 27: Annual Conference on Neural Information Processing Systems 2014, December 8-13 2014. Montreal, Quebee, Canada, pp. 3113 3121, 2014.
URL /ttt:::/papeee wwwppppppp.a.p.eee analysia
Arvind V. Mahankali and David P. Woodruff. Optimal & column subset selection and fast PTAS for low rank abs 2007 10307
approximation. CoRR, abs/2007 10307, 2020. URL httttp:://aaxxi..///
Michael Mahoney and Petros Drineas. Cur matrix decompositions for improved data analysis. Proceedings of the National Academy of Sciences of the United States of America, 106:697-702, 02 2009. doi: 10.1073/pnas. 0803205 106.
Andrew McGregor. Graph stream algorithms: A survey, SIGMOD Rec,, 43(1 900220 May 2014. ISSN 0163-5808. doi: 10. 00.1145/262769222676668 /tl.atpperr/b...orr/......e
7677644
Xiangrui Meng and Michael W. Mahoney. Low-distortion subspace embeddings in input-sparsity time and applications to robust linear regression. In Proceedings of the 45th Annual ACM Symposium on Theory of Compuring, pp. 91-100. ACM, 2013. doi: 10.1
Baharan Mirzasolciman, Ashwinkumar Badanidiyuru, Amin Karbasi, Jan Vondrák, and Andreas Krause. Lazier than lazy greedy. In Proceedings of the liventyNinthl AAAI Conference Artificial Intelligence, 15, pp. 1812-1818. AAAI Press, 2015. ISBN 0262511290
Jelani Nelson and Huy L. Nguyen. Osnap: Faster numerical lincar algebra algorithms via sparser subspace embeddings. In Proceedings of the 2013 IEEE 54th Annual Symposium on Foundations of Computer Science, FOCS '13, pp. 117-126, USA, 2013. IEEE Computer Society. ISBN 9780769551357. doi: 10.1109/FOCS. 2011..11.
Grigoris Paouris, Petros Valettas, and Joel Zinn. Random version of dvoretzky's 's theorem in C. Stochastic Processes and their Applications, 127(10):3187 3227, 2017. ISSN 0304-4149, doi: https://do..or//10. 1016/ j.spa.017..2..007
Christian Sohler and David P. Woodruff. Strong coresets for k-median and subspace approximation: Goodbye dimension. In Mikkel Thorup (cd.), 59th IEEE Annual Symposium on Foundations of Computer Science, FOCS 2018, Paris, France, October 7-9, 2018, pp. 802-813. IEEE Computer Society, 2018. doi: 10.1109) FOCS.200 8. 1.00081. URL/ttps: // crg. 1109// www.:0....1111.. 2000000002200011111111
Zhao Song. David Woodruff, and Peilin Zhong. Low rank approximation with entrywise 11-norm error. In Proceedings of the 49th Anmual ACM SIGACT Symposium on Theory of Computing, STOC 2017, pp- 688-701, New York, NY, USA, 2017. Association for Computing Machinery. ISBN 97819978100 450345286. doi: 10.1 011145/3055399.3055433
Zhao Song, David P. Woodruff, and Peilin Zhong. Average case column subset selection for entrywise & no-mr loss. In Advances in Neural Information Processing Systems 32; Annual Conference on Neural Information Processing Systems 2019, NeurlPS 2019, 3414 December 2019, Vancouver, BC, Canada, pp. 10111-10121, 2019a.
11

---

Under review as a conference paper at ICLR 2021
Zhao Song. David P. Woodruff, and Peilin Zhong. Towards zero-one law for column subset selection. In Advances in Neural Information Processing Systems 32: Anmual Conference on Neural Information Processing Systems 2019, NewrIPS 2019, 8-14 December 2019, Vancouver, BC, Canada, pp. 6120-6131, 6131, 2019b.
Ruosong Wang and David Woodruff. Tight bounds ip oblivious subspace embeddings. Proceedings of the Thirtieth Anmual ACM-SIAM Symposiun on Discrete Algorithms, SODA '19, pp. 1825-1843, USA, 2019. Society for Industrial and Applied Mathematics.
David P. Woodruff. Sketching as a tool for numerical linear algebra. Foundations and Trends in Theoretical Computer Science, 10(1 11-12))1-111 157, 2014. doi: 10.1561/440000060. URL htttps://ooi.....
[561/0400000060
Jiyan Yang, Yin-Lam Chow, Christopher Ré, and Michael W. Mahoney, Weighted sgd for lp regression with randomized preconditioning. Journal of Machine Learning Research, 18(211):1-43, 2018. URL
http: jmlr wwt.///g/pppp....e /vl w18/17--000 044. html
Linbin Yu, Miao Zhang, and Chris H. Q. Ding- An efficient algorithm for 11-norm principal component analysis. In 2012 IEEE International Conference on Acoustics, Speech and Signal Processing, ICASSP 20/2, Kyoto, Japan, March 25-30, 2012, pp. 1377-1380, IEEE, 2012, doi: 10. 109/CAASSP22222.28881444 URL
https: / /tw.///iii/// org/ / 10. 1109/ ICASS? 2012. 6288147.
12

---

Under review as a conference paper at ICLR 2021
A A HIGH COMMUNICATION COST PROTOCOL FOR k- k (p 1)
We start by describing protocol for distributed k-CSSp, which works for all 2 in the column partition model, and which achieves an O(k²) aapproximmaion the best rank- approximation, using O(1) rounds and polynomial time but requiring a communication cost that is linear in n + d. The inputs are a column- wise partitioned data matrix A E R g distributed across servers and a rank parameter & E N. Each server i holds part of the data matrix A, E R mdx=, Vi [s). and such that Ex 74
We use single machine, polynomial time bi-criteria &-CSS, algorithm as subroutine of the protocol, e.g.. Algorithm Chierichetti al. (2017). which selects a subset of Ō(k) columns Ar of the data matrix A Rdxa in polynomial time, for which minx |Ar. Alp I O(k) min Ag |A Aulp. Vp 2 1.
Algorithm 3 A protocol for k- C- (p 2 1)
Initial State: Server holds matrix A, Rdxx; Vi [s].
Server i:
Apply polynomial time bi-criteria k-CSS, on A, to obtain subset B, of columns as the left factor. Solve for the right factor V; arg miny, 1U;V; Ailp. Send U; and V; to the coordinator.
Coordinator:
Golumn-wiise concatenate the U, V. to obtain UV V U.Va]. Apply polynomial time bi-criteria k-CSS, algorithm on UV to obtain a subset C of columns. Send to each server.
Server i:
Solve minx, |CX; - Ailp to obtain the right factor.
Approximation Factor, Let UV denote the column-wise concatenation of the U.V. Let X arg miux [CX Alp. Then,
ICX" - PS S ICX" UV|2 + UV -Alp
By the triangle inequality By the O(p)aapproxiiatton of k-CSS,
< O(k) min
IUV -UVV) -VV Alp
(
S S(H)UV - Alp
= -0
U.V. - -..p))
< O(k) Σ O(k) min rank- A
IA Ailb)
By the apprrxxippoiiiooo of k-CSS,
0O(2 |4; - (A>)
A = arg min |A A'Ip makk-k 11-
= 0 A
Communication Cost. Since Us dxx(()) and I 0(0)xm sendingg U, and Vi costs Ö(sh)). Since C (dxx(()) sending from the coordinator to all servers costs Ö(skk). Thus overall communication cost is Ö(s(n + d))).
Running time. According to Chicrichetti et al. (017)) applying the k-CSS, algorithm and solving & regression can both be done in polynomial time. Thus the overall running time of the protocol is polynomial.
Problems with this protocol. Although this protocol works for all p 2 1, a communication cost that linearly depends the large dimension too high, and furthermore, the output C' is not subset of columns of A. because the protocol applies k- k-CSS,
concatenation of both the left factor U: and the right factor V. U; is subset of columns A, but V. not necessarily sampling matrix. One might wonder whether it is possible that each server only sends Ug and the coordinator then runs k-CSS, p on a concatenation of the U. This will not necessarily give a good approximation to min A, A Ax because the columns not selected in the U, locally on each server might become globally important. Finally, although it is possible to improve the approximation factor to Ō(k) by making use of an Õ Ō(VT) a approxxiiiiiio algorithm for fp-low rank approximation that also selects subset of columns (Mahankali & Woodruff] 2020), this protocol would still suffer from all of the aforementioned problems.
13

---

Under review as a conference paper at ICLR 2021
B PRELIMINARIES FOR OUR UPPER BOUND PROOFS
B.1 NORMS
Lemma 1. (Norm Relationships) For a matrix A E Rdd" |Alp.2 S |Alp and |Alp - - |A where
1 1sp<2.
Proof. Let R. For <p<r,
[x) F2lp S d $1/11
Let 2. Then we have
|x/2 |0lp S dd |0//2
[, |Anr1B) $ and || Eg lmm)13)
Note that Alp.2 =
Therefore,
Σ 1Aus1
S Σ
|Anil 1b
-1A1p
and
1
Alp
Σ
|Asl1
S
dd
B.2 LOWER BOUND ON THE COST No CONTRACTION WHEN APPLYING P-STABLE SKETCH
We show a lower bound on the approximation error for sketched subset columns, |SArV SAlp. in terms |ArV - Alp. The lower bound holds simultaneously for any arbitrary subset AT chosen columns, and for any arbitrary right factor V.
We begin the proof first showing that applying a dense p-stable sketch vector will not shrink its p-norm. This done Lemma 2.1 We further observe can still bound their tail probabilities by applying Lemma 9 from Meng & Mahoney (2013). We note this in Lemma 2.2 Note that the X,'s do not need to be independent in this lemma.
although p-stable random variables are heavy-tailed, we
Equipped with Lemma Lemma 2.2 and a net argument, we can now establish lower bound on (SArV SAlp. We first show in Lemma 2.3 that, with high probability, for any arbitrarily selected subset Ar of columns and for an arbitrary column A.j. the error incurred to fit SA., using the columns of SAr is no less than the error incurred fit A.j using the columns of Ar. We then apply union bound over all subsets T C [n] and columns [n] conclude lower bound in Lemma 2
Lemma 2.1. (No Contraction of p-stable Sketch) Given a matrix S E xR whose entries are i.i.d. p-stable
random variables rescaled by θ following holds:
where S p < 2, for any fixed R." with probubility 1 2/, the 87
1
1Sylp 2 lulz
Proof. By p-stability, we have |Sylß = Ei= (lb,, (=1) where [4, 2 0 are half' p-stable random variables.
Since Pr[[Z]] 2(1) > 2. by applying Chernoff bound (to the indicators 22,,2 for sufficiently small constant C), we have Σim I2.1" = 12(t) with probability 1 - 2. Therefore, with probability 0 |Syl, 2 lylp-
Lemma 2.2. (Upper Tail Inequality for p-stable Distributions) Letp (1,2), and m > 3. E [m], let X; be standard p-stable random variable, and let 7/s > O and 7 = Σm Let X = E X,|x, Then, for
any t 2 1, Pr(X 2 top7) < 2log(mt), where > 0 is a constant that is at most 2
Proof. Lemma from Meng & Mahoney 2013) for p (1,
14

---

Under review as a conference paper at ICLR 2021
Lemma 2.3. (No Contraction for All Sketched Subsets and Columns) Let A Id and k M. Let t poly(log nd), and let S E Rtx be a matrix whose entries are i.i.d. standard p-stable random variables, rescaled by @(1/t²). ). Finally, let m = k. poly(log k). Then, with probability 1 - for all T C [r]] with [T]] m. for all j [n], and for all y E!
Ary - A.sle S (S(ATV- - A.5)lp
We first extend Lemma 2.1 and argue that applying a p-stable sketching matrix S E Rtx= will not shrink the norm |Syle 2 vla (1 < < 2) simultancously for all in the column span of [Ar, A,l =: Ar. by net argument.
Proof. Step I:
In order bound the p-norm of sketched vectors y in net, we begin by showing that with high probability all entries in S are bounded. Let poly(mt). Consider the following two cases:
Case 1: p 1. The entries the 1-stable sketching matrix Siy are standard Cauchy random variables. Consider half-Cauchy random variables Xi, = |S The cumulative distribution function of half Cauchy random variables is F(x) = - नान di 1 ( Thus, for any i E [t] and [mm], Pr[[S;] s D]
1 1-(())-
Case p E (1, We apply upper tail bound of p-stable random variables Lemma 2.2 For any fixed i E and j E [m]. Pr|[Styl" I D*] 2 1- O loui))) which implies |r|] DJ 2 (
Therefore, for p E [1, 2). by a union bound over all entries in S. we define the event E to mean that for all i E [t] and j E [rn], we simultaneously have |Styl I D, then Pr[Ei] 2 1 O by union bound. The event E occurring implies that for any y E _2 since all entries in S are rescaled by ((1/t//))
| S 1y/ptt/p|s S) S Dlyle
3z E TT' s.t. y Ar.z2) in the column span of Ar.j. subset N B is net for for all y E there exists some N such that by - u/p S 7. for some distance > 0. There exists such net A for B of size [V (5) by standard greedy construction, since the column span of AT, has dimension most 1. We choose m- and thus
Consider the unit L, ball B (y R 1ulp
[W] < (m² D))(n) pm log ln)
By applying Lemma 2.1 and union bound over all vectors y N, we have that event & y e N simultancously, |Sylp 2 lulp Eg has probability at least 1 p00m log m)
Consider arbitrary unit vector E B. There exists some y E N such that 12- -12 S m Conditioning on both events & and event Eg, we have the following with probability 1 - O ml) - Du
|Sal p
|Sylp |S(x v)lp
Triangle Incquality By event &g Implication of event &
2 !ylp |5(x y)
2 1ulp D (xr 1))
2 1ælp DY
By jæ flp S
= Iulp O
= [ O m
|x
1l =
For sufficiently large 0( is at most 1, and thus
|ele 2
= m.. This implics |Sa| 2 |2/p
lalp
[xle We can rescale S by a factor of that [Sx], 2 |x|p-
2
We have shown that [Sylp 2 [ylp holds simultancously for all unit vectors y the column span Ar.. By linearity, we conclude that [ 2 lvlp (1 S p < 2) holds simultaneously for all in the column span of Ar.--
Step 2: Next, we apply union bound over all possible subsets T' C [n] of chosen columns from A and all possible single columns A., for j E [n], to argue that [S(Ary- - A=;)lp p 2 |Arg-A,, holds for all E R!7 and all T' < [m], C in] with high probability. Note that [T] = m = O(k: poly(log k)).
From Step 1, we have shown that event &g fails with probability pe log m)
Thus event &g fails over all
possible subsets T and all possible single columns A.j with prohahility most
Om log m)
d 1
Dmm log m)
d
((m)
et
15

---

Under review as a conference paper at ICLR 2021
The above failure probability is (11 long
= θ log Do lo do(n) d))
(OO (m log m) log + ( log d)
(O(polyyy)og k)) log(mnt) + Olkooly(log k) log(d))
Thus suffices to have = Apoly(log nd) have failure prohahility at most दoं...
Now condition single global event & and that event & holds for all possible T and A.g. We conclude that with probability 1 1- O p e0(s) the following holds simultancously for all T C [n) and for A., for which j [m]:
Ary - Auple & S(Ary Anj)lr.
Lemma 2. (Lower Bound for Sketched Error) Let A Bdx= and & E N. Let k. poly(log(nd)) and let S E R be matrix whose entries are i.i.d standard p-stable random variables, rescaled by ((1/t=) Then, with probability 1 for all C [n] with T = k poly(log k) and all V e B
|ATV Alp |SArV SALp
Proof. Let yy denote the j-th column of V, where j E [n]. By applying Lemma 2.3 and union bound over all columns of V. for = [T] = k poly(log k), = k poly(log(nn)) (nd)) and D poly(mf) = poly((: (ogd))))) the following holds with probability 1 - H( ml .oं9
o(1).
| Ale =
ArY - AイE)
1=1
|S(ATU) A))
1=1
= SATV SAlP
B.3 UPPER BOUND ON THE COST
We show upper bound the approximation
of k-CSSp sketched subset of columns, |SArV SAT which holds for fixed subset Ar of columns and for the minimizing right factor V arg miny [SArV- - SAppp for that subset of columns.
We first adapt Lemma E.17 from Song al. 2017) to establish an upper bound on the error |SArV - SAI for any fixed V in Lemma 3.1 We then apply Lemma 3.1 the minimizer V conclude the upper bound in Lemma
Lemma 3.1. (An Upper Bound on Norm of Sketched Matrix) Given A E Bn and E [1, 2), and E Raxk and V E paxd, if S E Rexa is a dense p-stable matrix, whose entries are rescaled by θ then with probability at least 1 O(1),
SSUV SAA < (llg(td))) AAE
Herr, the failure probability O(1) be arbitrarily small.
Proof. Lemma E.17 from Song et al. 2017
Lemma 3. (Upper Bound Sketched Error) Let E Rdx= and k N. Let = k poly(log(nd)), and let
gptxd a matrix whose entries are i.i.d. standard p-stable random variables, rescaled by @(1/t²). Then, for a fixed subset T' C [n) of columns with [T' = k poly(log k). with probability - O(1), we have
min SAAV SAlp < min
Alp
Here, the failure probubility can be un arbitrarily small constant.
16

---

Under review as a conference paper at ICLR 2021
Prouf. Let Xi arg |SArX- SAIp and X3 arg ming |Ar Alp. By Lemma 3.1
|SArXi -SA - A SAIE
(llog(nd)) AIE
Therefore,
min ISAr SAlp < min O(log"' (nd))) Alp
B.4 STRONG CORESETS IN THE p.22 NORM
Lemma 4 (Strong Coreset in fp.2 norm). Let A E xdx", ...nkk E (1,2), and c, E (0. Then, in poly(k log n/c) time, can find sampling and reweighting matrix T with O(d log dfe)) log(1/6) columns such that, with probability 6, for all rank-k matrices U,
min IUV AT (1 ronk V
min IUV Alp.2
rank
AT is called strong coreset of A.
Proof. We can obtain with (d((ogd columns using the strong coreset construction from Lemma 16 in SShler & Woodruff (0188). Note that the coreset construction for k-subspace approximation in Sohler & Woodruff (2018) aims removing dependence on d in the coreset size. The algorithm first finds poly(k)- dimensional subspace S by running dimensionality reduction algorithm and constructs coresets in the lower dimensional subspace S. resulting in coreset size poly(k/c), But our case, we not want our coreset size to have a polynomial dependency on & while a linear dependency on d suffices. Thus, instead of running their dimensionality reduction algorithm to find such a subspace $ to project A to, we directly use the the column span of the input matrix A as the subspace S appended with row of zeros to construct the B in Lemma 16 of (Shhler& & Woodruff (2018). Note that in the original algorithm, the appended column encodes the distances from the input matrix the subspace S. but our case it just 0's. Then the guarantees and running time claimed above immediately follow from Lemma 16 of Sohler & Woodruff (2018).
We note that the size of our coreset for k-subspace approximation can be further reduced to Ō(k). suppressing a logarithmic dependence on k, c: using additional ((nd) time, by combining Corollary 5.16 of Huang & Vishnoi (2020) and the importance sampling scheme from Stage of Algorithm I in Huang & Vishnoi (202 Furthermore, though not explicitly stated in Lemma 16 of Sohler & Woodruff (2018). the coreset size log(E) dependence on the failure probability ó due to the importance sampling (Sohler & Woodruff 2018 Huang & Vishnoi 2020).
Note that accuracy of (1 + <) is desired for extension the streaming model, described in Appendix B for both Algorithm and Algorithm| 7 we perform union bound over strong coreset constructions, for which a log(1/8) dependence on 8 is sufficient.
C POLYNOMIAL TIME, 0(1) APPPAAPPRO BI-CRITERIA k- &-CSSpp2
We give a detailed analysis on the polynomial time, appprpppoxeaa k-CSSp.2 algorithm presented in Algorithm 4 which is based idcas in Clarkson & Woodruff (015). We first sparse embedding matrix
to obtain O(1) aaprroxmtte left factor. We then 2,- Lewis weight sampling Cohen & Peng 2015 select subset of columns.
Algorithm 4 polynomial time, 0(1) appaaapproinnttoon for k-CSSp.2 (1 S < 2)
Input: The data matrix A E d rank & E N
Output: The left factor U (dx((()) the right factor V E I3 x such that IUV Alp.2 S ((1) miD As |Ag- Alp.2
S Ō(k) sparse embedding matrix, with sparsity s poly k).
S' t n x (kk) sampling matrix, each column of which is standard basis vector chosen randomly according to the lp Lewis weights of columns of SA.
return U + AS', V + (AS')'A { denotes the Moore- Penrose pseudoinverse.)
17

---

Under review as a conference paper at ICLR 2021
C.1 SPARSE EMBEDDING MATRICES
The sparse embedding matrix S xx())))d of Nelson Nguyen (2013). and used by Clarkson Woodruff (2015), is constructed as follows: each column of S has exactly s non-zero entries chosen in uniformly random locations, Each non-zero entry is a random value + <s with equal probability. s is also called the sparsity of S. Let h be the hash function that picks the location of the non-zero entries in each column of S and o be the hash function that determines the sign = of each con-zero entry.
Applying the sparse embedding matrix S to A enables to obtain rank-k right factor that is most factor of O(1) worse than the best rank-k approximation error in the 2,22 norm. We adapt Theorem 32 from Clarkson & Woodruff (2015) to show this in Theorem 5.5 Notice that in Theorem of Clarkson & Woodruff (2015).
the number of rows required for S is O(k*), but this can be reduced to O(k) through different choice of hyperparameters when constructing the sparse embedding matrix S.
We note two choices hyperparameters, i.e., the number m of rows and sparsity 8, of S in Theorem 5.1 and Theorem 5.2 both of which give the same result. The proof of Theorem 32 from Clarkson & Woodruff (2015
uses the hyperparameters from Thcorem 5.1 We replace the hyperparameters from Theorem 5.2 and show in Lemma 5.3 that (A) rows of S suffice to preserve certain desired properties. We then combine Lemma 5.3
and Lemma 5.4 adapted from Clarkson & Woodruff 2013). to conclude our result in Theorem 5.5
following
the analysis from Clarkson & Woodruft (2015).
Theorem 5.1. (Theorem from Nelson Nguyen 2013
For sparse embedding matrix S E Rmxx with
sparsity s = and a data matrix E R xn let E (0,
With probability at least 1 8 all singular values
of SU are (1+ e) as long as m. 2 8 (d dd) 4-wise independent and h is pairwise independent.
For the hash functions used to construct S, a is
Theorem 5.2. (Theorem 9 from Nelson & Nguyen 2013
For sparse embedding matrix S E m= with sparsity s = (log (d/8)/c) /c) and a data matrix U E Rudd e E (0, With probability at least ô all singular values of SU are (1 +) as long as m = (d (d/6)/e²). For the hash functions used to construct S. we have that o, h are both ((llg(d/0))-ww
independent.
Lemma 5.3. Let C be a construint set and A E Rand faadddp
B E Ruxdd be two arbitrary matrices. For a spurse
dlog 1 pert
embedding matrix S RmXx" there is m = 0 21++1)
), such that with constant probability, the following
holds for all X E gaxdd
i) SSAX B) 2 (1 (1-1))) c)1xx Blp.2
it SAAX - B)nn² & (1 + +) - Bl, where X" arg min |AX Blp.
X
Proof. The proof is the same as the proof of Lemma 29 from Clarkson & Woodruff 2015), except that we use a different choice of hyperparameters in constructing S, i.c., sparsity and the number m of rows. In the proof of Lemma 29 from Clarkson & Woodruff (2015), the construction of S follows Theorem 5.1 where the sparsity s I, but requires m O(d²) rows. replace the construction by Theorem 5.2 where we pick 8 eo+1,
Now the sparsity s larger but this construction reduces the number of rows required to m Ödd). Since both choices of hyperparameters to construct result in bounded (1 + e) singular values SU for any data matrix U. the rest of the proof follows.
Lemma 5.4. Consider a data matrix A E Raxd_ Let the best rank-k matrix in the tg2 norm be As arg min A2 LAk - Alp.2. For R Rdxm, if RD satisfies both of the following two conditions for all X E Roxe;
(AE N22 (1 -c)1 A² A²Ib.2
i) R² (AfX - bb2 S (1 + +e] x A lp.2, AFlp.2, where X* -argmin
min |A[ x A
then
min A ATE $ (1 3c) - -AT 1. Kunk X
Proof. Lemma 31 Clomkson & Woodruff 2015).
Theorem 5.5. (ow Rank Approximation) Let the data matrix be Rdx= and &s E N be the desired rank. Let B be sparse embedding matrix with m O(kpoly(log k)poly( ()) rows, and sparsity = poly(log k). Then, the following holds with constant probability,
min IXSA - Alp.2 I (1+ 3e)
min AAu p.2 Aq
18

---

Under review as a conference paper at ICLR 2021
Prouf. The proof is the same the proof of Theorem 32 in Clarkson & Woodruff (015), except that adapt a different construction of the sparse embedding matrix S, which reduces the number of rows from O(k*)
(kk) with increased sparsity
Consider Ag arg min mink Ag. |Ag - Alp.2. V2 he a basis for column space of Ag. By applying Lemma a.3 Lemma 5.4 the basis Vs. we conclude the above theorem by setting the number m of rows
to 172 tog (2)
and sparsity = poly(log k) the sparse embedding matrix
C.2 LEWIS WEIGHTS SAMPLING
The tp Lewis weight is an inherent property of matrix. By Cohen & Peng (2015). the unique set of fp Lewis weights for matrix is defined follows: for the i-th column of matrix M gaxd M.j. defined as u (., wtwpl///ppg y Mag. where W is the diagonal matrix with W.. w,
[d].
We first note that sampling by Lewis weights provides subspace embedding in the é, -norm in Theorem 5.6 We further apply randomized version of Dvoretzky's Theorem in Theorem 5.7 which allows the embedding from the fp norm into low-dimensional Euclidean subspace with very small distortion and thus enables us to switch between the bp norm and the 222 norm. Based on Theorem [56 and Theorem [5.7] show in Theorem 5.88 that Lewis weights sampling provides good subset columns, on which the analysis of k-CSSp.2 based.
Theorem 5.6. Lee-L Weights Subspace Embedding) Let A E Raxd and Ö(d). For < < < there exists distribution (A1, Au..... >.. on the rows of A such that if we generate a matrix with rows, each chosen independenly as the th standard basis vector times (r)a)²
with probability Xa then with probability
1 - o(1), the following holds for all CR.
(21) I |SAx|, S ((1)
Proof. Theorem 7.1 from Cohen & Peng 2015
Theorem 5.7. (Randomized Dvoretzky's Theorem) Let n. E N. and E (0, 1). Letr = 5. Let x be random matrix whose entries are i.i.d. standard Gaussian random variables, rescaled by + For 2, the following holds with probability 1 - 1 (((n), for all y E R!..
Gyy (1 = e)ly12
Proof. This follows from Theorem 1.2 from Paouris et al. 2017
Theorem 5.8. (Subset of Columns by Lewis Weights Sampling) Let A E Rdxa, Let S Raxd be sparse embedding matrix, with m O(k poly(log k k))pl((?) .)). Further, let S' Ruxt be a sampling matrix whose columns are random standard basis vectors generated according to the Ep Lewis weights of columns of S A. with k poly(log k). Then, for X- arg miiu b XXSAS AAS p.2, the following holds with constant probability,
İSA Alp.2 S O(1)
min At
AAe Alp
Proof. Let X = arg min X X SA Alp.2. 2. By the triangle inequality,
I&SA - A Alp.2 <
X SA &SA + [X SA Alp.2
Let consider bounding SA X.S.p...
By Lemma D.28 and Lemma D.29 from Song al. 2011 for any column sampling matrix S and for any our case, since S' is sampling matrix, we have
fixed matrix Y. can be shown that E [ = E|JyS s']g]] = SY for any fixed matrix Y.
rescaled random matrix whose entries are i.i.d. standard Gaussian random variables as in Theorem 5.7] We apply Theorem 5.7 transform between the ép space and the Euclidean space. Since transformation of both directions can be done with very small distortion, we obtain @(1) approximation. With constant probability, have
Now let d(d)xd
SX*SA- XSA p.2 (ONG(X"-XSAAA
By Theorem 5.7
19

---

Under review as a conference paper at ICLR 2021
(X)) X)SAAs
By Theorem 5.6
=(I))(XX X]SA5"1p. 1p.2
By Theorem 5.7
s (I SAS' - AS" 18 + IXSAS - -A52)1.))
Triangle Inequality
s e(1]IXT SAS" - AS 12.2
Since 8 arg min [XSAS* - AAS"| pp.2
mak-1
#(L)GG(X" SAA-AAA) 1
By Theorem 5.7
- ((1)GG( -CA >)>-
By Markov Bound on E S |Y S
+(1)]X SA- Alp.2
By Theorem 5.7
Therefore,
IXSA - Alp.2 1 IX'SA-XXXAA + +X*SA SA AAppa
$ (()]X""S A/,.2
< (1) min IA - Au
By Theorem 5.5
rank- Ax
C.3 ANALYSIS FOR k-CSS
We are now ready to show O(1) approximation factor and polynomial running time for the bi-criteria C-CSSp.2 presented in Algorithm given in Theorem
Theorem 5 (Bicriteria O(1)-Approximatio Algorithm for k-CSSp.2). Let A E g and k E N. There exists algorithm that nuss in (nnz( A) + df) -kpolylogg k) time and outputs rescaled subset of columns U E (xx((() of A and right factor V RO(k) for which V minv IUV Alp.2. such that with probability (1))
IUV Alp.2 S 0(1)
min |Aa Alp.2
sand-k Ag.
Proof. Approximation Factor First notice that the minimizer X of ISSAS AS p2 has to be in the column span of AS". Thus we can write X = Y for some matrix Y. By Theorem 5.8
P&SA - Alp.2 = ((ASS) Ale.2 S @(1) min IA Aslp.2 rank Ag
We denote YSA= V. We take the left factor U = AS" and solving for miny IUV Alp.2 will give ((1) approximation min A, |4 Aslp. A good minimizer for the right factor V in the Euclidean space V (AS')*A A. This concludes our result. Notice that since S' is sampling matrix with Ö(k) columns, get a rank-k left factor U as a subset columns of A as desired.
Running time First notice that S a sparse embedding matrix with poly(log k) non-zero entries. Thus computing SA takes time (n)(A) poly(log ) By Cohen & (2015), computing the Lewis weights of SA takes time anz( (SA) + poly(k) < nnz(A))poly(log takes time nnz( A). Computing ( takes time d² kpooyylog k). Computing the right factor (AS")' A takes (n))(())pol((lM) (log k). Therefore, the overall running time is (nnz( A) + d²) kpoly(log k).
poly(k), and computing the output left factor U = AS
D ANALYSIS FOR GREEDY k:- J-CSSS1.
We propose a greedy algorithm for selecting columns k-CSS presented in Algorithm ??. We give a detailed analysis on the first additive approximation compared to the best possible subset of columms for Greedy k-CSS 1.2. Our analysis is inspired by the analysis of Greedy k-CSSa for the Frobenius norm in Altschuler et al. (0016). We then describe how the running time of Algorithm ?? can be improved to Ō( F)), where F' is the running time required to evaluate miny [ATuy V A)1,2 for a fixed j E T (note that this running time is Onx + knd) we evaluate this by computing the pseudo-inverse of ATuj). This improvement in the running time is obtained by randomly sampling candidate columns from T and adding the best these randomly sampled columns to Ar, rather than trying all columns of 1' this method was previously used in Altschuule et al. (2016) for Frobenius norm column subset selection, where it is referred to as "Lazier-than than-lazy Greedy," and the general approach was first introduced in Mirzasoleiman al. (2015). This version of the greedy algorithm is shown in Algorithm 2 and show both versions of the greedy algorithm helow for convenience.
20

---

Under review as a conference paper at ICLR 2021
Algorithm 5 Greedy k-CSS),2 1,2. Here we denote the set of selected columns 7 from A by Ar and the set of unselected columns by Ar.
Input: The data matrix A E R the number of iterations r S 7. Output: A subset of columns Ar from A, where [7]
Ar + 0
for i to do Column j* + arg min Jeay (miny |ATuı - 1122))
Ar + Arus
end for
Algorithm Lazier-than-lazy Greedy k-CSS1,2 This version the greedy algorithm is based on
Section 5.2 of Altschuler et al. (2016).
Input: The data matrix A E Rddx, the number of iterations r S n, a parameter E (0, 1).
Output: A subset of columns Ar from A, where [T] r.
Ar + 0
for i 1 to do
T + A subset of (lg((1 columns of A, each selected uniformly at random (excluding the columns whose indices are T)
Column j" + arg minjea, (miny |Aru, - -1122)
AT t- Aruj
end for
extended Algorithm 6 We first show in Lemma 7.2 an improvement the utility function with one additional column when projecting a single vector, based on Lemma 7.1 from Altschuler et al. (2016). We then show an improvement of the utility function when projecting a matrix in Lemma 7.3 by applying Lemma 7.2 and Jensen's 's Inequality, following the analysis in Altschuler al. (016)) Finally, conclude analysis for Greedy 1-CSS1,2 Theorem
First analyze Algorithm then show how this analysis
Notation Consider the input matrix A E Rx (n d). Let B be the matrix of normalized columns of A. projection onto the column span
where the j-th column
B.j = A.5/ [A,5l2. Let mr R -> Rd of Ar or equivalently Br. Let omi (M) denote the minimum singular value of some matrix M.
To aid our analysis, we define utility function as follows, inspired by Altschuler et al. (2016). For subset [n] and matrix E Rext (or vector E IRd),
PM (T) = |M IM - MTM Σ (IMail |M.s - vr.li Σ (T)
Observe that as the number of columns selected and added to I' increases, we get more accurate estimation of M and thus the approximation error |M कт M|1,2 decreases, which results in an increase in the utility function 4 x5 (T).
Lemma 7.1. Let S,T C [rx] be two sets of column indices, with I 2 | for some vector u Rd. Then,
( ( 2 mmin (Bs) [msu u3))
4|
Proof. Lemma from Altschuler et al.
except that we replace the condition for S and T', i.e. P 2 T 2 Irru/2 The two conditions equivalent, since
2016
(T') Altschuuel et al. 2016) with
(s,(S) 2 (T)
ju msul2 Iu Tru/2
|u/3
|T 3 < |u|2 - 1|
/su 2 Trul2
21

---

Under review as a conference paper at ICLR 2021
Lemma 7.2. (Unilit Improvement by Projecting Single Vector) Let S,T [Tx] be two sets of column indices, with (S) 2 (T) for some vector u E R.. Let lo = [S]. For i E [k], let T. TU {i}. Then,
(T.(T:) - (T)) 2 mmin (Bs)² (2) D (T))
(S)²
Proef. We define function for convenience in the analysis (-o0, |u[3] - R20 by g(x) = vluB
Note that gg [ar)] 2V/u13=2
(T,()) (T)(
Σ [I TTuU22 wu To.
By definition of +
(Vlul3 krul3 Vlul3 क.र 3
By Pythagorean Theorem
(o(ITru()) - ((
1o'((TT ((trw(/) ( |Tru(E)
By the Mcan Value Theorem
Σ 2 2//ulp |rrul
| (
By Pythagorean Theorem
2/u Truu 1=1
u Tmen (Bs) (TT m
2 2/u - TT u/2
By Lemma 7.1
mmin (Bs) (/u
|uu s3)
By Pythagorean Theorem
8 TTS
TT 4/2 f/u
Omin (Bs)²
8 | TTU/2
(Iu - iT u/2 - |u1 T(/22)²
(lu - uTul2 + Iu T
mmn (Bs) (B0) (S)
- {T))))² |u =
Since |u Tgul2 2
|
We now lower bound as follows.
Ju-xyull Ingulg
Trul2
Tru l2
By Pythagorean Theorem
u
a Tsu
2
mrul2 la
1122
|u T
[u2 +
1u
TsU
1
By definition of
(S)
|1/2
u
Iu
Trul2
Since Ju rsulp S |u/a
(8)
|M/2
|u/2
(T)
By definition of
24.(()
1u/2
(T)
(S)
|u/2
(1
T. (T)
Since (.. S Ju/2
(S)
o. (S))
22

---

Under review as a conference paper at ICLR 2021
(S) (a(1)) (5)2
Therefore,
k
(a((T)) - pu(TT) 2 min (Bs)² (2,(S) (f.(S (T)) ( -/u--lu mr
8 |
min ( (Bg) (2n(
( (T)))
16)
Lemma 7.3. (Utility Improvement by Projecting Matrix) Let A E Rd*x", and T, C [m]] be two sets of column indices, with 1 (S) 2 4 * (T). Furthermore, let = |S|. Then, there exists a column index E such that
A (i}) xA(T) 2 Omi (Bs)² (#A(S) ( (T()))
A())
Proof. The proof mostly follows the proof of Lemma in Altschuuer al. (2016). We combine Lemma 7.2 with Jensen's 's inequality conclude an improvement the utility function with one additional column when projecting matrix instead a single column.
For [n], define 6, min(1, 2.., (T) TA., (S) Note that 6, 1 if the j-th column A., has larger projection
onto than Bs. and $A.; (2) B.. (3) otherwise. Let (S). For e [k], let T7 (U).
(Qx(2)) - -Qx
1
P (11) - 2 2 (1))
By definition of +
Omin (Bs)
Omin (Bs)2
1=1
6,)²
A (5)
By Lemma 7.2
16
① (S)
4 A (S)
(1 - 8,)3
Note 4 =
AA. (S)
16
Σn (A,,
j=1
(
$ (S)
(1 61)
By Jensen's Inequality
16
Σ 4A (S)
1=
(1 6z) 2A. (5) 3
4 P (T) Since 1 1-dj 11
(S)
A.A. (S)
1=1
+1 - 8,) 4 A,, (S)
2 DA., (S) D (T)
160 AS)2 Σ (A,, )A,
(DA(S) (T))
16)2
Hence,
Σ - -2a(T) 2 mmn (Bs)² ($) (T))
-MA(T)))((TTT)
164 (S))
This implics there least column of Bs, with index E S, such that when is added T' the utility (8)((8)4441))) 164. x(8)
function 4 (T) incrcases by at least 1 Tmii (Bs)²
Theorem 7. Let A E Rdxn be the data matrix and & E N. be the desired rank. Let Ag be the best possible subset of &s columns, Le., As arg min Ag miny A Al1.2- Let be the minimum non-zero singular
23

---

Under review as a conference paper at ICLR 2021
value the matrix B of normalized columns of As. (the j-th column B B., (ds)=;/((As)) (1s))) Then, T C [r]] is the subset of columns selected by Greedy k-CSS1,2. the following holds with [2]] 92( (2(==)))
min LAV Al1.2 s (1 c)
min -cl.).18-kk.Vee VeRkkk |AsV - Al1, + +
Proof. The proof follows the one of Theorem in Altschuler et al. (016)).
Let be the subset of columns of selected by Greedy k-CSS after t iterations. Notice that To = Ø. In addition, define F' ( @A((S) @a(Tb), De F, and A Au for i N.
A. 2 (S) (T)) 2 A.+1 2 Our goal bound the number iterations needed the gap between (S) and 4 (T,) to become less than A.
Consider some iteration adding a column of to T, through greedy selection is at least the improvement of the utility function after adding the best column of Bş to T.. By Lemma 7.3
which 4 ( 1 (T( 2 A. The improvement the utility function after
this is least
a ($) (x(1))
o² A
164 (S))
K.
128
the gap (S) (T,) after iterations at most A, and at least A then after at most 64kp2 eAA?
iterations,
the gap becomes most A.
We can use this to bound the number of iterations required for the gap to become at most eF. Take N E N such that A++1 oF < AN. Then the number of iterations required for the gap become most A AN+1 is at most
64kF GZA?
64kF
N
1
A
E=
F4
A
Since AN+1
o² A N41
4N+1-1
2N+1-1
256k
Since AN+1 2 cF
and Σ
4N+1-1 S
1
30
2
Therefore, after [T] = 9 x iterations, we have
a (T) v x$)
|A ne Al1.1 (lAh122 IA Ar c(A41.2 - IA As.)1.
A A11 (1 -)] गड A/1.1 <|Al1.2
Since S is the set of indices for the best possible subset of & columns, the above is equivalent
min A/1,2 I (1 -)
min Vee
|AsV Al1,2 + <|A/1,
We now analyze Algorithm based on the analysis of the Lazier- than- lazy greedy heuristic in Altschuler et al. (2016). The first step is the tollowing lemma based on Lemma 6 Altschule al. (0116), which shows that in expectation, the utility improves by a large amount on each iteration.
Lemma 7.4 (Expected Increase in Utility - Based on Lemma 6 of Altschuler et al. 2016)). ). Let E Rddn. and ler T.S C [n] be two sets of column indices, with & == |S and (S) 2 (T))). Ler T be set of 1og(1/8 column indices of A, chosen uniformly at random from [r]] T. Then, k
E Elmax de A (t)) (2 (1 5) mmnn (Bs) (DA(S) (T)) (S)2
24

---

Under review as a conference paper at ICLR 2021
Prouf. The proof nearly identical to the proof of Lemma 6 of Altschuler et al. proof for completeness. The first step the proof showing that ( \ I') is nonempty with high probability, show that the expected increase in utility is large.
2016
include the full
Then, by conditioning 7 n (S 7) being nonempty, For the purpose of this analysis, we assume that the columns of T are sampled independently with replacement. At the end of the proof, we discuss sampling the columns of 'I without replacement.
First, observe that
(
S T T
P = =
S T
21 (o)(11
(1/
S
By 1 - 28 S =
18 1/001/4)
I e
Because n - 1 n
meaning that
18 1oc(1/6)
(s)
ISV
2 (1 5) SS T
Since IS \ TI < k. and 1 8 2 (1 5) for x. 5 E [0,
Therefore,
Elmax Φ (TU (1)) (T)
EE
2 P 7) + Ø| E max max muN.:./rutt00--.r////11111
[S\
(1
5)
F
max IET
- (T
k
S \ (1 5)
max SET
C - (T
k
(Since is always better for T n (5' to be larger)
(1 - 8) S Eves a4 (TU (i)) (T))
k
| [
EIEs(ON(TU(R)) - (T))
(1 - 5)
[S]
(Since ( (T) for E T)
1 2 (1-
Jmn (Bs)² "mmn(Bs)
(S)
(T))²
[S]
(S)²
(See the proof of Lemma 7.3
This proves the lemma the case where columns are sampled with replacement. Now, we discuss what happens when sampling without replacement. Note that the expected increase in utility can only be higher if the columns of T are sampled without replacement, Intuitively, this is because if T has some repeated columns. then it is always better to replace those repeated columns with other columns of A. Thus, for each instance of T'' where some columns are sampled multiple times, we can "move" all the probability mass from this instance of to other sets T' [n] T. which contain T' but do not have repeated elements. This leads to the uniform distribution on subsets of [m) \ with no repeated elements, i.c., the distribution that results from sampling without replacement.
Using this lemma, we analyze the convergence of Algorithm 2
Theorem 8. Ler A Rxd be the data matrix and E N the desired rank. Let Ag be the best subset of l
columns, Le., As = arg min miny |AsV - A|12- Let a be the minimum non-zero singular value of the matrix B of normalized columns of Ag (meaning the j-th column of B is B.; = (As)-,/ / (As).sl2) Then, if
T [n]] is the subset of columns selected by Algorithm
2 the following holds if [T] 02 ((,,))))
E' G'min A - -11,2] < (1= -)
min =t.
|AsV - A/1.2 + <|A/1,2
25

---

Under review as a conference paper at ICLR 2021
Prouf. The proof is uses the same strategy that of Theorem 5 of Altschule et al. (2016) (and Theorem above), with minor modifications. Let T. be the subset of columns selected by Algorithm 2 after iterations (in particular, To 2). In addition, let F = 4 ( (S) (TT)), Do F. and A.+1 = a/2 Now,
fix a time t such that for some i, A: 2 ( (1i) 2 A.41 2.2 Then, we bound the number of additional iterations t' needed so that
(S) 4 (T(+e) Ti] A
For convenience, for each & 2 0, define Eb :- 2 (Ti+2). Then, our goal is to find 1' such that
4 Ev < Ait1
However, observe that from Lemma 7.4 above, we obtain
Ex+1 Es- E 4 (1+) (TT(T++)) T:
E E 21 (++1)) 2x (Ti++)) T(tt Ti
By B|E[X B[X]
(BA(S) (a(Te+)))
2 E (1 - 8) mmnn (Bs)²
Ti
By Lemma 7.4
AS))
(1 5) amen (Bs)
Fl (((()(()) 4b (T(++)))
16k (S))
8) Jmin (Bs)²
[F($A(S) (A(T+++)
By Jensen's 's Inequality
a
= (1 (1-0) Jmin (B²)² (E( - Ex)
16k A (5)2
Now, suppose that A, 2 (S) - E, 2 Act1. for s = 0,... I' - 1. Then, for all such s, E.s1 - Eg 2
(1-8)mmin (Bs) 16kF2
Summing these inequalities for 0, .,... 1' 1, we find that
Fy F. 2 (1 m (Bs)² A +'' 16k F2
and for the increase from Eg to Ep to be greater than A.n. it suffices to have
3
t' 2 A (1 (1- 5)0 (Bs)²
In summary, if (S) (Ti) S A,, then in at most 11p (1:-(1110mma 3)0 (Ba) iterations, (S
32
AT))(T))] S Ai+1. Thus, if we let N E N such that AN+1 S V F S AN. then the number of iterations needed to have (S) - El (T:)) < AN+1 is at most
F2kF2
3
111 (1 mm BB)²
(1 - 8) 8) Bs))
A
32kF2
(1 8) 6) )3
NN
A
1=
F2k
1(1
6)
N
1
eF 2
Σ
Since ANT1
8) (Bs)²
F2
-N
1=0
N
128k
Σ
omno (Bs)
4
i=D
512k
3amin 22
Thus, after = 0
Fmire (pg) iterations.
6
*x - ((Tx) S
-1 xS(S)
c
meaning
|A/1. - |A - A/12 - F/1| - 1A - - S c
A|/1. - v1 V1-51A
14 =s Al1,2
v1 0 6
and rearranging gives
B कz A11.2] < (1 - 1 A VI-8 A/1,
e
This completes the proof (note that can select e, meaning فعا = 0(1) for e < 5
26

---

Under review as a conference paper at ICLR 2021
E EXTENSION TO THE STREAMING MODEL
In this section, describe how our protocol in Algorithm T can be made into 1-pass streaming algorithm for column subset selection in the l, norm. The algorithm shown Algorithm 7 and is analyzed in Theorem
below. The algorithm and its analysis follow the standard merge-and-reduce framework (see McGregor (014))
Theorem 9 (Analysis of Algorithm 7.
Ler A E Baxn and k E N, and assume Algorithm 7
sees the columns of
A one at a time in the stream S. Then, Algorithm 7 returns U E Rdxx(() such that
min VeRO(e)xn
|UVV - Alb S 0(k//11/2)
min || Alla
Au rank
with probability 0.9. Moreover, the space complexity of Algorithm is Ö(dk).
Proof. Let (()) be the bi-criteria rank of Algorithm 4 Then, will he repeatedly applying Lemma with & being equal to r, i.c., we will create coresets which preserve the errors when projecting onto all subspaces of dimension r.
Now, at every iteration of Algorithm 7 each element (L, of L can be thought of as holding coreset of the columns of A in some interval I in [rd. We prove the following intermediate lemma by induction t:
Lemma 10. Let B be the concatenation of all the sketched columns L (which have been reweighted by multiple applications of Lemma 4). Then, for all subspaces V C COOK) of dimension most
BB - P. = (1 (1 log ml |SSA, Py.AAllb.2
Proof. We proceed by induction on f. The lemma is clear when = 0, since that case, columns of B are simply sketched columns of A which have not been re-weighted.
Now, suppose > 0, and suppose the lomma holds for smaller values of f. Then, the sketched columns in L must have been obtained as follows: there previously existed two elements (L1, f - 1) and (L2, 1) of L, such that if B3 is the concatenation of the sketched columns in L1 and B3 the concatenation of the sketched columns in coreset for concatenation B3 B3 and B2, then and L2 form coresets for two intervals in [n). which we denote by I and 1z respectively. Applying Lemma
L2, and
B3T. The sketched columns L
we find that
P
B - Pv Blp, =
| - Pv Ball,.2
log
+
"(1B: -1 PvB: 11 Pv B2l
log
+
(1 log (ISSAr PvSSAr 11..2 + +|SAA1 S.rr (1.2)
p(e-1)
(l..2)
log
=
+ log
||SA, PvSArll
(I)
where the second and last equalities are because the p" power of the fp.2 norm of a matrix decomposes across the columns of the matrix, and the third equality is by the induction hypothesis. By taking 8t roots, find that
BB - PvBlE = (1 log m
) [|SAA Sv.
This proves the lemma.
Now, if is the unique element L remaining at the end concatenation of the sketched columns in L. Then, is a coreset for all the columns of A, and by the above lemma, the distortion of / at most (1 + log I jlog" E [t, e] (since t S log nogn-me nn---o note that, as in all applications over the course of algorithm form a binary tree, with the leaf nodes being contiguous intervals of length (k:)). Hence, for all subspaces V of dimension most r,
for loop Algorithm 7 let be
of the merge-and-reduce framework, the coresets contained
B Bl S(1] - Pv
and in particular, if M E R ((())x(()) is the matrix formed by concatenating the columns of B selected by running Algorithm 4 then
ISA PM @(1)
min Tis
|SA P,SSAl|p.2
27

---

Under review as a conference paper at ICLR 2021
where Par the left-hand side denotes the projection the column span of M, and Pr the right-hand side denotes the projection to the column span of S Ar. Hence, by Lemma
SSA SMSA I 0(k 12-1/22 min TCM. St
|SA - PrSAl1>
meaning that if U E (dx((() is the matrix whose columns are the unsketched columns corresponding to M, then
min Vd
SA - SUVII -
min Ve
I[AA 15r VlIp
and by Lemmas
and
this means
min VER
||A-UVII S
min I AzVIlp
Now, we analyze the space complexity of Algorithm| 7 Note that any iteration of the algorithm, C can hold at most log n lists of columns (since each element of L 1s coresct corresponding
interval of column indices [n] of size 2 for some E N, and if two adjacent coresets are of the same size then they will have been merged, all coresets in L are of different sizes). Each list of columns is of size most (A)) and each column has d entries, meaning the amount of space used in any iteration is at most (dk)).
F COMPARISON OF OUR PROTOCOL WITH THE DISTRIBUTED GREEDY PROTOCOL OF ALTSCHULER ET AL. (2016) FOR THE FROBENIUS NORM
this section, we perform empirical comparison our protocol with the distributed greedy protocol for column subset selection in the Frobenius norm due to Alischuler et al. (2016).
F.1 DISTRIBUTED GREEDY PROTOCOL OF ALTSCHULEER ET AL.. (2016)
We first recall the distributed protocol of Altschuler et al. (2016), in Algorithm 8 Note that it is bi-criteria algorithm, i.e., more than & columns are selected, and in Altschuler et al. (016)) it is shown that this gives good approximation relative to the optimal column subset for the Frobenius norm.
Naively, this would require large communication cost, since the entire data matrix would have be com- municated between the servers and the coordinator, in order for them to perform the calls GREEDY (4, T; 32k and GREEDY A, T, 12A ). Instead, Altschuler et al. (2016) uses the technique of Projection- Preserving
sketches, developed in Cohen ct al. 2015), in which key result is as follows:
Theorem 11 (Projection-Costt Preserving Sketches Theorem 4 of Cohen et al. 2015), stated in Altschuler
al. 2016 Let R random matrix with n rows and n'' O (+
columns, where each entry
independently and uniformly set to 1 v # Then for any matrix E Rdx= with probability O6)) the
following holds: for some constant 2 for any k-dimensional subspace U of of Pu E gadd is the corresponding projection matrix, then
(1 - P S (
that is, AR is a Projection Cost Preserving sketch for A. In other words, R preserves the cost of all projections of A onto rank-k subspaces R...
Hence, when each server performs the call GREEDY T,, 32%), ), in place of A, it serving sketch AR for A (and the coordinator similarly use AR when it makes the call GREEDY(A, T, 12k )). This is because the calls to GREEDY repeatedly compute the cost of the projection of matrix onto various column subsets of T; or T. and hence that needed is way to efficiently compute the cost the projection of A onto various subsets of size (
projection-cost pre-
AR can be computed without communicating the entire matrix A between the servers, as follows:
First, the coordinator generates R E =xan (where n' = O(k)) specified in Theorem 11 This is sent all the servers (which can be done with negligible cost if a random seed is sent, for example).
Suppose the th server holds A.. which consists columns with indices s; through
A. Then, R; is the (f, f + x suhmatrix R consisting rows 8; through f;, then the 7th server can send A R; to the coordinator, Since A, R2 is a d x ((k) matrix, this step takes ((sck) communication.
Finally, by definition of matrix multiplication, AR Σi A,R. The coordinator performs this computation and sends AR each server. This step also takes (sdkk) communication.
28

---

Under review as a conference paper at ICLR 2021
Algorithm 7 1-pass streaming algorithm for k-CSSp. is collection of coresets of columns. Over the course of the algorithm, each element for some E [log ni] to determine when two coresets should be merged, we will also keep track of the size of each coreset L. Hence, each element of is the form (L, 1) where
will represent a contiguous subset of 2 columns,
list of sketched columns (and their corresponding unsketched columns) and t is the number of times this list has been involved a merging operation.
which the columns of the data matrix Rd arrive one at time, and the target rank & E N
Input: stream
Output: The left factor U Rdx((())
S An O(k) d random matrix with i.i.d. standard p-stable entries
log - )
f + O(k) C+ Ø
for Each column Auy that arrives from S do
if' L empty. then
L + (}. where O is the empty list.
else if The last element (J,t) t) of L is such that = 0 (i.e., it has been merged O times). then
L J
Remove (J,t) from the end of C.
else
L 0
end if
L LU ((SA.j: A.y))
C + (U,
1* Now, we merge coresets in L as much as possible. */ while True do
Exit this while loop if L has only element.
tL,t) + last element of L remove this from
(L',L') <- last element of L - remove this from C.
If 3f /', then add (L', 1') back to L, and add (I,, back to Then, exit this while loop.
C A coreset of L, U T., computed specified in Lemma [ the & in the statement of Lemma 4 will be the bicriteria rank of Algorithm 4 which is (k)). The parameters 5 and will be as specified at the beginning of this algorithm. (Only compute the coreset of the sketched columns for those which are included in the coreset, include the corresponding unsketched columns as well (but without re- re-weighting them)).
+ (C, + 1) note that for each of the re-scaled columns that are included in C, include their original indices in A as well. end while end for
L + The unique element of C L' t The result of rumning Algorithm 4 on the sketched columns for each of these sketched columns, store the corresponding unsketched column well.
U The d Ö(k) matrix whose columns are the unsketched columns in L' return U
29

---

Under review as a conference paper at ICLR 2021
Algorithm 8 Distributed Greedy Column Subset Selection for the Frobenius Norm (Algorithm 2 of Altschuler et al. (0166). Throughout this protocol, GREEDY A, T.r) denotes a single-machine procedure, which does the following: if A E Rd is data matrix, and T R is set of columns (not necessarily of A) then a subset S of r columns of T is constructed iteratively over r steps, such that is, the chosen column
that at each step, the new column of I to add to S is greedily chosen increases g A13 the most, or equivalently, decreases |A - TT A the most. In other words, it is the same as Algorithm ??, but for the Frobenius norm rather than the 1122mmomm.
The data matrix A E _ target rank k E N, the number of servers & E N. The columns of A are assumed to be randomly partitioned among servers T1, T2...., ,T...
Input:
Output: A subset of columns Ar from A. [T] O(*) where if Aopr is the optimal subset of columns of A of size k, and the columns of Aopr are normalized to have unit &g norm, then a is the smallest singular value of AOPT.
S: < GREEDY A, Tú 32k for all i E [s] (The jth server performs this computation.)
Each server sends its St to the coordinator.
T - S; (This computation is done by the coordinator.)
S + (A, T, 12k
The coordinator returns arg
A
we make use of Projection Cost Preserving sketches (PCPs). Another optimization described in Altschuler et al. (2016), which we use in our implementation, is the LAZIER- THAN-LAZY-GREEDY ALGORITHM. The difference between the GREEDY (A, T,r) algorithm and the SEATTR-THAN-LAZV-RRRERE (A, r) algorithm as follows: while at each the iterations, GREEDY considers all columns of T' and chooses the one that leads to the most improvement in the objective, LAZIER-
In our implementation of Algorithm 8,
THAN-LAZY- GREEDY samples (og()) columns uniformly random for some small 8 (which we take in implementation to be 0.005), and out those columns, chooses one which leads to greatest improvement (where [T' the number of columns in 2). This leads to a significant speedup to GREEDY, and was shown in Altschuler et al. (2016) that this does not significantly worsen the approximation guarantees that Altschuler et al. (2016) shows for GREEDY.
F.2 SETUP
We compare our distributed protocol, in the case 1. to Algorithm 8 For both protocols, at the outset we fix the number of columns which are selected. Algorithm is rewritten Algorithm to reflect this. In this section, we use & denote the number of columns ultimately selected.
Algorithm 9 Distributed Greedy Column Subset Selection for the Frobenius Norm (Algorithm 2 of Altschuler et al. (0066)). For our empirical comparison, we fix the number of columns selected at the outset.
Input: The data matrix A Rdxa, desired number of columns & E N, the number of servers E N. The columns of A are assumed to be randomly partitioned among servers Ti, T2..... T..
Output: A subset of columns AT from A, with [T]
AR < PCP for A discussed in the previous section. All servers and the coordinator have access to AR.
S; + GREEDY (AAR, T,, kt) for all E [s] (The 1th server performs this computation.)
Each server sends its S; to the coordinator.
T S (This computation is done by the coordinator.)
S + GREEDY(AL) T)
The coordinator returns arg "'S. TT (AR) 2
Note that in Algorithm 9. the coordinator now chooses as many columns as chosen by the servers, as opposed to Algorithm g where it chooses somewhat smaller number of columns note that this cannot harm the performance of the algorithm.
We compare our protocol to Algorithm 9 using the following datasets:
gastro_les 76 698 matrix dataset available awaishhealiiccpss//arhhw... ml/datasets//aatt inal htt.r///eee..........aooaaooooen
.ddu/
30

---

Under review as a conference paper at ICLR 2021
secom, 591 1567 matrix dataset. secom has missing entries, which we replace with Os for the purposes of our experiments. secom available at htttt:://acciive. uucs.ucii... wwt.:uuc..eeuuuueeeeee edu/m1/ datasets/5BCOm
F.2.1 PARAMETERS
We compare our distributed protocol for K-CSS1, using Greedy K-CSS1,2 2 as a subroutine, with Algorithm
9 for several choices of & on gastro.lesions, we let & E {10, 20, 30), while on secom, we let For our protocol, cauchy ,si size is set to 2k and coreset_si ze is set to 5k for both datasets (where cauchy size and coreset size ze have the same meanings as in the main body of this
K E (30, 60, 90, 120].
paper.) For Algorithm 9 secom.
the number of columns in the PCPs set 8k on gastro lesions ins, and 7h on
The hyperparameters coreset_siee and cauchy-si -szee and the number of columns in the PCPs, are set this way so that the amount of communication that each algorithm is allowed is roughly equal (Algorithm 9 allowed slightly more communication). To see this, let d he the numher of rows in the data matrix A. If e is the number of columns in the PCPs, then the total communication used to transmit the PCPs hetween servers is 2sod, since the servers must first send their respective A, Rg to the coordinator, and the coordinator then sends A to all of the servers. By comparison, if ri is the number of rows in the initial Cauchy matrix in our protocol, and v2 the number of columns in each coreset sent by the servers to the coordinator. then the total communication required to transmit the Cauchy matrix and the coresets (r1 + ra)sd. We choose rı, r2 and co that 2scd slightly higher than (r1 + r2)ad.
Note that if & is the number of columns ultimately selected, then protocol uses 2dk additional bits of communication between the servers and the coordinator to recover the final ku-ubset of columns, while Algorithm will use adk communication to send the subsets of columms S: (of size k) from the servers to the coordinators. This is not included in our hyperparameter calculations.
F.2.2 How TRIALS ARE CONDUCTED
With these choices of hyperparameters, we conduct 15 trials as follows. In each trial, if A E RR is our data matrix, then we shuffle the columns and then divide them equally between 2 servers (since for the theoretical guarantees of Algorithms 8 to apply, the columns should be partitioned randomly). Using this partition across 2 servers, perform multiple-respons. & regression to evaluate the error.
our protocol and Algorithm 9 For each protocol, once the subset of columns computed,
In the next section, we report the minimum error across the 15 trials for each datasct and for cach value of k. We also report the mean error, along with the standard deviation. Finally, compute the work and the span of our protocol and Algorithm 9 taken to perform multiple-response ε regression.
using Python's time. prwcesssttnnee () utility. This does include the time
For secom, with k = 30 and k = 60, the trials were performed on a L Macbook Pro with a 2.7 GHz Quad-Cor Intel Core 17 processor and 16GM of memory. The rest were performed on 2019 MacBook Pro with a 2.4 GHz Intel Core i5 processor and GB 2133 MHz LPDDR3 memory.
F.3 RESULTT
Our results are shown and discussed in Figures (minimum and mean/standard deviation for & error) and 4 (average work and span across 15 trials).
G FULL. EXPERIMENTAL RESULTS
For our protocol, we considered various additional hyperparameter settings on real-world datasets (ocsst k isolet and images in the Caltech- 101 dataset) these settings are shown in Tables 3 4 and As hefore, cauchy size is the number of rows in our initial Cauchy matrix, sent to the servers at the beginning the protocol, and coreset a iize is the size of the strong coreset sent by each server to the coordinator. For regular k-CSS1,2. sketch size is the number of rows in the sparse embedding matrix, and sparsity is the number of nonzero entries in each column of the sparse embedding matrix.
For each these settings, we display the minimum error, as well as the mean error and the standard deviation for each setting and each rank, as shown in Figure 5 We also display the average work/span of each setting, for each rank, as shown in Figure 6 We observe that not only does greedy k-CSS1,2 2 perform better than all settings of regular k-CSS1,2 in minimum and mean errors across multiple trials, it also has smaller variance in performance.
G.1 DETAILS OF WORK/SPAN COMPUTATION
In F work and span were recorded as the time taken (using Python's t time. process.t... st ime () utility) to compute the column subset in our distributed protocol. In particular, it does not include the time spent performing regression obtain the right factor and consequently the entry-wise f errors.
31

---

Under review as a conference paper at ICLR 2021
Nin Enor on qasaro Jesions
Erao on gastro estons
103300
200000
#0000
1500033
00000
100003
50003
10
20
(a) rastro. min error
(b) gastro. inn: mean std error
Min Error
Mean
1900000
Fratocel
3500000
Pnotacoi
0000000
6200000
4000030
1003320
4000000
3000000
3000300
2000000
3000000
Park
Hant
(c) secom: min
secom: mean & std error
Figure 3: Plots (a) and (c) show minimum l error across 15 trials on gastro_les ons and secom respectively, while plots (b) and (d) show the mean & errors and the standard deviation. For ranks 90 and 120, our protocol affords a 3% improvement in minimum error secom, while for rank 30 on gast ro..rrolleese ionss our protocol gives over 40% improvement. Note that the mean error for both protocols is noticeably higher gastro lesleiiinn in the case of our protocol, the error in the case k- = 30 is distorted by two trials with E errors 205827 and 134945 respectively.
Sotting Number
sketch size
5k
k/2
8/3
R75 (//(2
X/3 min(5. A:33
R75
sparsity
min 36)
min(20,5k))
mini 3k)
min 40, 5k)
min k/2)
min(2,
min/5. A/2)
min(5, m/3)
Table 3: All hyperparameters used on bcsstk13, when regular k-CSS is used. In this table, k denotes the number of columns ultimately selected. In all settings (including when greedy k-CSS not included here, is used), cauchy si ze is either 5k or 8k, and coreset size is 5k. The setting numbers are shown on top in the following error plots, Setting is used to refer to our protocol when greedy k-CSS1, is used. Note that sparsity can be at most sketch size, since sketch size is the number of rows of the sparse embedding matrix, while sparsity is the number of nonzero entries in any column. (There is a slight typo in Table 2 in Section 8 of our main paper, where coreset si ze is given as 10k.)
Setting Namber
10
sketch size
K75 min(2, A/S)
K72
K/3
K75
sparsity
min(20, 36)
min(20),6 [k] mini (k) min (4/44) 5k)
min 1/22)
min 2,n322
mini6,
min A/3)
min(G.
Table 4: All hyperparameters used on isolet, when regular k-CSS1,2 is used. all settings (including when greedy k: k-CSS1,2 2, not included here, is used), cauchy size is 4k; and coreset size is 4k. The setting numbers are shown on top - in the following error plots, Setting 0 is used to refer to our protocol when greedy k-CSS1, 1 is used.
32

---

Under review as a conference paper at ICLR 2021
Average work on gestro Josiens Using time. timel
Average span gastro lesions Using time ame.timal
Delatte Busiiru
Frebenss Frotocal Our mrsocol
Fraberi Prthooi
Our rrotcce
20
18
Hower
(a)
Average Work
(b) gastro.lesions: Average Span
Aveezgg work on seccm Likig time I mmme)
Average span on speam Living time time()
3000
1400
1200
2000
1000
800
200
Trelean Protooo our Pradcct
Tacerasx Freloxo our mosacol
J83
(c) secom: Average Work
(d) secom: Average Span error
Figure 4: Plots (a) and (b) show the average work and span respectively for gastro ro les ons, while plots (c) and (d) show the average work and span for secom. As expected, our protocol using GREEDY k-CSS takes more time - to our knowledge, there is not yet an optimization similar to LAZIER-THAN-LAZ GREEDY for the l fnorm... Nevertheless, running time is less important than communication in the distributed setting.
33

---

Under review as a conference paper at ICLR 2021
1014
Beoste13
2224
Sceetic15
1.75
2.8
Bank
bosstkl3 min error
(b) besstkls mcan std error
Isolet
kolle:
eropeeoo emor M00 330
0.0
0.
Rane
kank
(c) 1sol eet min error
(d) 8.so1 et: mean & std error
Cahteh101
Pattectd
1.0
1.2
tr0yyloo DIT0 MI
8+
127
(c) caltech-100 101: min error
(f) cawtechh-10 101: mean & std error
Figure 5: Results on besstk13, 3, isolet, and caltech-101 from top to bottom. The left plots show minimum deviation in error. In all plots, the first bar denotes SVD. the second bar denotes greedy k-CSS1.2 and the rest of the bars denote all settings of regular k-CSS 1,2, at all ranks on the axis 10 through 60.
across all 15 trails; the right plots show the corresponding mean and standard
34

---

Under review as a conference paper at ICLR 2021
105
Bcostk13
3costi13
Averaae Span
Rine
(a) bosstk13 3: work
(b) bosstk13 span
Isokt
Iolle:
122
PNtegg Noh
0.-
0.
24
0.
0.
20
30
40
00
Punk
Pank
(c) isolet: work
(d) isolet: span
Caltec1011
Caltech101
30
255
2.99
2.0
10
Park
(e) caltech-1011 work
(f) caltech-101:
Figure 6: Results on besstk13, isolet, and caltech-100 101 from top to bottom. The left plots show average work in seconds across all 15 trails; the right plots show the corresponding average span. In all plots, the first bar denotes greedy k-CSS1,, and the rest of the bars denote all settings of regular i-- 1,, at all ranks the axis 10 through
35

---

Under review as a conference paper at ICLR 2021
Setting Number
k73
k75 miin A/G)
476 minn k/5)
spereity
min
min A/3)
mintG A/31
Table 5: All hyperparameters used on isolet, when regular k-CSS is used. In all settings (including when greedy k-CSSS,2 2, not included here, is used), cauchy size 4k, and coreset size is 4k. The setting numbers are shown on top in the following error plots, Setting 0 is used refer to our protocol when greedy k-CSS 1,2 is used.
Figure shows that encounter tradeoff between accuracy and running time when choosing hetween these suhroutines, Both lead the overall communication cost, and since accuracy is of more interest than nunning time in distributed setting, (empirically) preferable to use Greedy k-CSS1,2 within the protocol.
G.2 ADDITIONAL DETAILS
All experiments on the Calt Malttech-100 dataset were run on a L016 Macbook Pro with 2.7 GHz Quad- Core Intel Core 17 processor and 16GM of memory. All experiments on bosstk13 3 and isolet were run on AWS z xlarxe instance with Deep Learning AMI (Amazon I. Linux 2) Version 29.0.
36


Reviews:
[{"review_id": "oV8uZ3CDIL6", "reviewer": "AnonReviewer4", "summary": "", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "SUMMARY:\n\nThe paper considers the column subset selection (CSS) problem, which has received considerable attention in numerical linear algebra. It considers a distributed variant of CSS in the $\\ell_p$ norm, where $p \\in [1,2)$. Despite the attention this problem has received previously, it seems like this is a new setting that has not been considered before. The paper primarily provides theoretical results, but also does some experiments. \n\nI think this is an interesting paper and problem. However, there are a few things I'd like to see addressed before I can give a higher score. There are a few details in the proof of Theorem 6 that are unclear. Moreover, the readability of the paper could be improved. I provide further details below.\n\n\nADVANTAGES:\n\n- CSS is an important problem in numerical linear algebra.\n- The setting under consideration seems to be new.\n- The paper provides both theoretical and experimental results.\n\n\nCONCERNS/QUESTIONS:\n\n- Section 1 would be easier to follow if the definitions in Section 2.1 were given when the concepts are first used in Section 1 instead.\n\n- The p-stable distribution is discussed throughout the paper, but doesn't seem to be defined anywhere, not even in the supplement. It would be helpful if the p-stable distribution was defined when it is first mentioned. Also, in the proof of Lemma 3.1, it would be helpful if you could point to a reference that contain the relevant facts about p-stable distributions.\n\n- In Section 1.2, the first paragraph, you make a distinction between an approximation to OPT and an approximation to the best possible subset of columns. Can you clarify the distinction between these two types of approximation?\n\n- The numbering of the lemmas is confusing. The numberings used in the main paper are reused in the supplement, but for different results. For example, Lemma 1 in the main paper is called Lemma 2 in the supplement, and Lemma 2 in the main paper is called Lemma 3 in the supplement. This makes the proof of Theorem 6 confusing.\n\n- On the 4th line in Section 1.1., you say that the coordinator computes $\\sum_i A_i S = A S$. But $A$ is of size $d \\times n$, and each $A_i$ is of size $d \\times n_i$, so this equality isn't right. Do you mean to say that $\\sum_i \\tilde{A}_i S = A S$, where each $\\tilde{A}_i$ is of size $d \\times n$ and with all entries zero except those columns corresponding to the columns in $A_i$ (or something to that effect)?\n\n- In the proof of Theorem 6, when you apply Lemma 5 the first time: Don't you need both $V'$ and $M$ to be defined in terms of the $p,2$-norm rather than the $p$ norm for Lemma 5 to apply? In other words, should they be defined as \n\n$V' := \\arg \\min_V |S A_I V - S A|_{p,2}$ \n\nand \n\n$M := \\arg \\min_M | S A_I M - S A T |_{p,2}$?\n\n- In the proof of Theorem 6, when you apply Alg. 4: It is not clear how Alg. 4 from the supplement is applied here. Alg. 4 outputs factors U, V such that UV approximates A. What is $(SAT)^*$? Is $(SAT)^* = UV$ here?\n\n- In the proof of Theorem 6, when you apply Lemma 5 the second time: It looks like this uses the first equality in Lemma 5 in the supplement, which requires $P_2^*$ to be a projection matrix. How can we know that it is? Moreover, even if the sampling matrices $T_i$ are chosen to satisfy Lemma 5, how do we know that the matrix $T$ obtained from the $T_i$'s in Alg. 1 also satisfies Lemma 5?\n\n- What is the run time of Algorithm 2? Does finding each $j^*$ require evaluating the objective for every choice of $j \\in \\overline{T}$, or is there a better way to do this?\n\n- In Section 6, in the Setup paragraph, you say that you report the minimum error over 15 trials. Is there a reason for reporting the minimum rather than the mean or median? Is the mean/median also competitive?\n\n- The paper ends abruptly with no conclusion.\n\n\nMINOR CONCERNS/QUESTIONS:\n\n- In the sentence right after Lemma 2 in the main paper, should $|S A_T V - S A_T|_p$ be $|S A_T V - S A|_p$?\n\n- In Algorithm 2, should the argmin be over $j \\in \\overline{T}$ rather than $j \\in A_{\\bar{T}}$?\n\n- In Tables 3, 4 and 5 of the supplement, should k-CSS_{!,2} be k-CSS_{1,2} (i.e., a '1' instead of an exclamation point)?\n\n\n###################################\n\nUpdate:\n\nThe authors have addressed the concerns I had in my initial review. I have raised the score from 6 to 7.", "extra_scores": {}, "extra_text": {}}, {"review_id": "gdFJWIC6U8L", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The authors present a distributed protocol for the column subset selection problem under the $\\ell_p$ norms for $1\\leq p <2$. The model of computation is a standard coordinator model of communication where the input matrix A is column partitioned to $s$ servers. The main contribution consists of a protocol that requires $O(sdk)$ communication, single round and polynomial running time and returns a $k^{1/p - 1/2}$ multiplicative approximation to the best possible column subset.\n\n#### Reasons for score:\nOverall, I vote for a borderline / slightly below acceptable but I am open to discussion about my decision. My main concerns are (a) the polynomial time complexity of the method makes the proposed protocol hard to be useful in real applications, and (b) the experimental evaluation is done on very small matrices (i.e., 2000 x 2000 shape) as well as small number of servers. On the positive side, the paper is a solid theoretical work on the distributed CSSP problem.\n\n#### Strong points:\n* Well written paper with a clear contribution statement\n* Related work is up-to-date (to the best of my knowledge)\n* Concise algorithm description and corresponding theoretical guarantees.\n\n\n#### Concerns: \n* The polynomial running time guarantee makes the practicality of the proposed algorithm marginal. Do you (implicitly) assume that you have limited space per server as well? Please make such assumption explicit.\n* Experimental section need sto be extended with larger datasets and more servers. How does the distributed protocol compares with a centralized one? The size of the datasets here suggest the following question.\n* Why is the distributed setting is studied here? If the resulting algorithm is polynomial, why not solving the problem in a single server?\n\nMinor comments:\n* Explicitly state what is hidden under the notation $\\tilde{O}$.\n*Lemma 3: the probability constant $0.999$ seems arbitrary here, please replace it with $0<\\delta<1$ failure probability parameter.\n* Lemma $4$: What is a \"sampling and reweighting matrix T\"? Is this defined in the main text?\n* After Theorem 5: give a reference to \"Lewis weights\".\n* Before Section 5 (Runtime): Why generating a p-state random variable takes constant time?\n* Is all communication bounds in bits?\n* Figure 2 is not easy to read, please make it larger and increase the font on both axis.", "extra_scores": {}, "extra_text": {}}, {"review_id": "Q0ViRejXKA6", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper focused $\\ell_p$-norm error analysis of distributed column subset selection problems. The proposed distributed algorithm guarantee an $\\tilde O(k^{1/p-1/2} )$-approximation to the best subset of columns with $\\tilde O(sdk)$ communication cost per round and polynomial time.  The lower bound analysis show that such communication complexity is near optimal. The design and analysis of the algorithm is based on f oblivious sketching and strong coresets, which leads to the overall approximation factor is only $\\tilde O(k^{1/p-1/2} )$. The authors also extended the main idea to address $\\ell_{1,2}$ norm error. \n\nI have not checked all details of the proof, but the results and the main ideas look reasonable. I think the theoretical contribution of this paper is valuable, however, the contents is somewhat limited to the audience of ICLR. There are some questions/comments:\n\n1. I strongly encourage the authors apply the proposed algorithm to address the real-world problems such as computer vision, image processing which is mentioned in introduction.\n\n2. Since the column select problem focus on the matrix $A\\in\\mathbb{R}^{d\\times n}$ with $n\\gg d$, the empirical studied should also contain such case. It is observed that proposed method performs better than SVD on both synthetic dataset and some real world dataset. I think one of the reason may be the $l_p$ norm for $1\\leq p<2$ is not a unitary invariant norm and SVD do not give an optimal low-rank approximation. Hence, SVD may be not a good deterministic algorithm as baseline. \n\n3. The presentation could be improved: a) The operations on Server and coordinates should be described in multiple lines.  b) The curves in figure 2 should be displayed with different types of lines.", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: Clarity: Well written paper with a clear contribution statement; related work is up-to-date; concise algorithm description and corresponding theoretical guarantees. However, the presentation could be still improved.

Significance: The polynomial running time guarantee makes the practicality of the proposed algorithm marginal. Experimental results do not back up strongly the significance of the algorithm.

Main pros:
- Solid theoretical work on the distributed CSSP problem.

Main cons:
- The reviewers point out the significance of the experimental results, beyond the theoretical contribution. For the ICLR audience, real, large-scale experiments, that really dictate that the proposed solutions is (if not the only) one of the few solutions to follow, are necessary. The reviewers highlight that the theoretical results need to be applied to really large-scale scenarios (e.g., the problems considered in the paper can definitely be handled by a single computer, and no distributed implementation is required). 
- Τhe polynomial time complexity of the method makes the proposed protocol hard to be useful in real applications.
- How does the distributed protocol compares with a centralized one? This is not fully addressed in the rebuttal. 
Decision: Reject