Year: 2020

Abstract:
We give a distributed protocol with nearly-optimal communication and number of rounds for Column Subset Selection with respect to the entrywise `1 norm (k-CSS1), and more generally, for the `p-norm with 1 ≤ p < 2. We study matrix factorization in `1-norm loss, rather than the more standard Frobenius norm loss, because the `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 ∈ Rd×n (n d) are distributed across the s servers. We give a protocol in this model with Õ(sdk) communication, 1 round, and polynomial running time, and which achieves a multiplicative k 1 p− 1 2 poly(log nd)-approximation to the best possible column subset. A key ingredient in our proof is the reduction to the `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 `1,2 norm. Finally, we implement our protocol and give significant practical advantages on real-world data analysis tasks.

text: 1 p− 1 2 poly(log nd)-approximation to the best possible column subset.
A key ingredient in our proof is the reduction to the `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 `1,2 norm. Finally, we implement our protocol and give significant practical advantages on real-world data analysis tasks.

heading: 1 INTRODUCTION
text: Column Subset Selection (k-CSS) is a widely studied approach for rank-k approximation and feature selection. In k-CSS, one seeks a small subset U ∈ Rd×k of k columns of a data matrix A ∈ Rd×n, typically n d, for which there is a right factor V such that |UV −A| is small under some norm | · |. k-CSS is a 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.
k-CSS has been extensively studied in the Frobenius norm (Guruswami & Sinop, 2012; Boutsidis et al., 2014; Boutsidis & Woodruff, 2017; Boutsidis et al., 2008) and operator norms (Halko et al., 2011; Woodruff, 2014). A number of recent works (Song et al., 2017; Chierichetti et al., 2017; Dan et al., 2019; Ban et al., 2019; Mahankali & Woodruff, 2020) studied this problem in the `p norm (k-CSSp) for 1 ≤ p < 2. The `1 norm is less sensitive to outliers, and better at handling missing data and non-Gaussian noise, than the Frobenius norm (Song et al., 2017). Specifically, the `1 norm leads to improved performance in many real-world applications, such as structure-from-motion (Ke & Kanade, 2005) and image denoising (Yu et al., 2012).
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 on multiple machines. Despite the flurry of recent work on k-CSSp, 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 a number of results in the distributed model are known, see, e.g., Altschuler et al. (2016); Balcan et al. (2015; 2016); Boutsidis et al. (2016). We consider a widely applicable model in the distributed setting, where s
servers communicate to a 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 a factor of 2 and an additive log s bits per message to identify the destination server.
We consider the column partition model, in which each column of A ∈ Rd×n is held by exactly one server. The column partition model is widely-studied and arises naturally in many real world scenarios such as federated learning (Farahat et al., 2013; Altschuler et al., 2016; Liang et al., 2014). In the column partition model, we typically have n d, i.e., A has many more columns than rows. Hence, we desire a protocol for distributed k-CSSp that has a communication cost that is only logarithmic in the large dimension n, as well as fast running time. In addition, it is important that our protocol only uses a small constant number of communication rounds (meaning back-and-forth exchanges between servers and the coordinator). Indeed, otherwise, the servers and coordinator would need to interact more, making the protocol sensitive to failures in the machines, e.g., if they go offline. Further, a 1-round 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 1.
In the following, we denote Ai∗ and A∗j as the i-th row and j-th column of A respectively, for i ∈ [d], j ∈ [n]. We denote AT as the subset of columns of A with indices in T ⊆ [n]. The entrywise `p-norm of A is |A|p = ( ∑d i=1 ∑n j=1 |Aij |p) 1 p . The `p,2 norm is defined as |A|p,2 =
( ∑d j=1 |A∗j | p 2) 1 p . We consider 1 ≤ p < 2. We denote the best rank-k approximation error for A in `p norm by OPT := minrank-kAk |A−Ak|p. Given an integer k > 0, we say U ∈ Rd×k, V ∈ Rk×n are the left and right factors of a rank-k factorization for A in the `p norm with approximation factor α if |UV −A|p ≤ α ·OPT.
Since general rank-k approximation in `1 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 k columns, such algorithms return a subset of Õ(k) columns of A, suppressing logarithmic factors in k or n. It is known that the best approximation factor to OPT that can be obtained through the span of a column subset of size Õ(k) is Ω(k1/2−γ) for p = 1 (Song et al., 2017) and Ω(k1/p−1/2−γ) for p ∈ (1, 2) (Mahankali & Woodruff, 2020), where γ is an arbitrarily small constant.
1We provide a detailed analysis of our streaming algorithm in Appendix E.
1.1 PREVIOUS APPROACHES TO k-CSSp IN THE DISTRIBUTED SETTING
If one only wants to obtain a good left factor U , and not necessarily a column subset of A, in the column partition model, one could simply sketch the columns of Ai by applying an oblivious sketching matrix S on each server. Each server sends Ai · S to the coordinator. The coordinator obtains U = AS as a column-wise concatenation of AiS’s. Song et al. (2017) showed that AS achieves an Õ( √ k) approximation to OPT, and this protocol only requires Õ(sdk) communication, O(1) rounds and polynomial running time. However, while AS is a 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 to the Lewis weights (Cohen & Peng, 2015) of AS to 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 Õ(k3/2) approximation to OPT (Song et 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-CSSp algorithms to the distributed setting under the column partition model. Existing works on polynomial time k-CSSp (Chierichetti et al., 2017; Song et al., 2019b; Dan et al., 2019; Mahankali & Woodruff, 2020) give bi-criteria algorithms, and are based on a recursive framework with multiple rounds, which is as follows: in each round, Õ(k) columns are selected uniformly at random, and with high probability, the selected columns can provide a good approximation to a constant fraction of all columns of A. Among the remaining columns that are not well approximated, Õ(k) columns are recursively selected until all columns of A are well approximated, resulting in a total of O(log n) rounds.
A naı̈ve extension of this bi-criteria k-CSSp framework to a distributed protocol requires O(log 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 O(log n) round selection procedure on the coordinator only. In order to do this, the coordinator needs to first collect all columns of A from the servers, but directly communicating all columns is prohibitive.
Alternatively, one could first apply k-CSSp on Ai to obtain factors Ui and Vi on each server, and then send the coordinator all of the Ui and Vi. The coordinator then column-wise stacks the UiVi to obtain U · V and selects Õ(k) columns from U · V . Even though this protocol applies to all p ≥ 1, it achieves a loose O(k2) approximation to OPT and requires a prohibitive O(n+ d) communication cost2. One could instead try to just communicate the matrices Ui to the coordinator, which results in much less communication, but this no longer gives a good approximation. Indeed, while each Ui serves as a good approximation locally, there may be columns that are locally not important, but become globally important when all of the matrices Ai are put together. What is really needed here is a small coreset Ci for each Ai so that if one concatenates all of the Ci to obtain C, any good column subset of the coreset C corresponds to a good column subset for A. Unfortunately, coresets for the entrywise `p-norm are not known to exist.

heading: 1.2 OUR CONTRIBUTIONS
text: Our Distributed Protocol We overcome these problems and propose the first efficient protocol for distributed k-CSSp (1 ≤ p < 2) in the column partition model that selects Õ(k) columns of A achieving an Õ(k1/p−1/2)-approximation to the best possible subset of columns and requires only Õ(sdk) communication cost, 1 round and polynomial time. Figure 1 gives an overview of the protocol. We note that our subset of columns does not necessarily achieve an Õ(k1/p−1/2)-approximation to OPT itself, although it does achieve such an approximation to the best possible subset of columns. Using the fact that there always exists a subset of columns providing an Õ(k1/p−1/2)-approximation to OPT (Song et al., 2017), we conclude that our subset of columns achieves an Õ(k2/p−1)-approximation to OPT. Recently, and independently of our work, Mahankali & Woodruff (2020) show how to obtain a subset of columns achieving an Õ(k1/p−1/2)-approximation to OPT itself; however, such a subset is found by uniformly sampling columns in O(log n) adaptive rounds using the recursive sampling
2We give this protocol and the analysis in Appendix A
framework above, and is inherently hard to implement in a distributed setting with fewer rounds. In contrast, our protocol achieves 1 round of communication, which is optimal.
We make use of a strong coreset, i.e., a sampled and reweighted subset of columns of each Ai that approximates the cost of all potential left factors of Ai, by first embedding all subspaces spanned by any subset of Õ(k) columns of A from `p-space to Euclidean space, to bypass the lack of strong coresets for the `p norm. We denote this new norm as `p,2 norm, which is the sum of the p-th powers of the `2 norms of the columns. To reduce the error incurred by switching to the `p,2-norm, we reduce the row dimension of A by left-multiplying by an oblivious sketching matrix S shared across servers, resulting in an overall approximation factor of only Õ(k1/p−1/2). Afterwards, each server sends its strong coreset to the coordinator. The coordinator, upon receiving the coresets from each server, runs an O(1)-approximate bi-criteria k-CSSp,2 algorithm to select the final column subset, giving an overall Õ(k1/p−1/2) approximation to the best column subset.
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 `p,2 norm (Sohler & Woodruff, 2018; Huang & Vishnoi, 2020) to distributed k-CSS. Furthermore, to show that our oblivious sketching step only increases the final approximation error by a logarithmic factor, we combine a 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 (2020), only consider a single subspace at a time.
Theoretical Guarantees and Empirical Benefits for Greedy k-CSS1,2 We also propose a 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 is at most (1 − ) minV |ASV − A|1,2 + |A|1,2. Similar error guarantees were known for the Frobenius norm (Altschuler et al., 2016), though nothing was known for the `1,2 norm. We also implement our protocol and experiment with distributed k-CSS1 on various real-world datasets. We compare the O(1)-approximate bi-criteria k-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.

heading: 2 PROBLEM SETTING

heading: 2.1 THE COLUMN PARTITION MODEL
text: We consider a model where there are s servers, the ith of which holds Ai ∈ Rd×ni , and a coordinator which initially does not hold any data. Each server talks only to the coordinator, via a 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 is O(log(snd)) bits. The overall data matrix A ∈ Rd×n is A = [A1, A2, . . . , As] (the column-wise concatenation of the Ai’s). Here, n is defined to be ∑s i=1 ni. Typically, in the column partition model, n d.
2.2 p-STABLE DISTRIBUTION AND p-STABLE RANDOM VARIABLES
Let Z,X1, . . . , Xd be random variables drawn i.i.d. from some distribution D. D is called p-stable if for an arbitrary vector v ∈ Rd, 〈v,X〉 = |v|pZ, where X = [X1, . . . , Xd]T . Random variables drawn from such a distribution are called p-stable random variables, which exist for p ∈ (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 θ ∈ [−π2 , π 2 ] and r ∈ [0, 1] are sampled uniformly at random, then, sin(pθ) cos1/p θ ( cos(θ(1−p)) ln( 1r ) ) 1−p p follows a p-stable distribution.

heading: 3 PRELIMINARIES FOR OUR PROTOCOL
text: We first note a standard relationship between the `p norm and the `p,2 norm.
Lemma 1. For a matrix A ∈ Rd×n and p ∈ [1, 2), |A|p,2 ≤ |A|p ≤ d 1 p− 1 2 |A|p,2.
3.1 `p-NORM OBLIVIOUS SKETCHING
We left-multiply A by an oblivious sketching matrix S with p-stable random variables so that we lose only an Õ(k 1 p− 1 2 ) approximation factor when we switch to the `p,2 norm.
The purpose of the next two lemmas is to show that we 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, |SATV − SA|p, 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 A ∈ Rd×n and k ∈ N. Let t = k · poly(log(nd)), and let S ∈ Rt×d be a matrix whose entries are i.i.d. standard p-stable random variables, rescaled by Θ(1/t 1 p ). Then, with probability 1− o(1), for all T ⊂ [n] with |T | = k · poly(log k) and for all V ∈ R|T |×n, |ATV −A|p ≤ |SATV − SA|p.
Next, we show an upper bound on the approximation error of k-CSSp on a sketched subset of columns, |SATV − SA|p, which holds for a fixed subset of columns AT and a fixed right factor V .
Lemma 3 (Sketched Error Upper Bound). Let A ∈ Rd×n and k ∈ N. Let t = k · poly(log(nd)), and let S ∈ Rt×d be a matrix whose entries are i.i.d. standard p-stable random variables, rescaled by Θ(1/t 1 p ). Then, for a fixed subset T ⊂ [n] of columns with |T | = k · poly(log k) and a fixed V ∈ R|T |×n, with probability 1− o(1), we have minV |SATV − SA|p ≤ minV O(log1/p(nd))|ATV − A|p.
3.2 STRONG CORESETS IN THE `p,2 NORM
To enable sub-linear communication cost in the number n of columns, the i-th server sends the coordinator a strong coreset of columns of SAi, which is a reweighted subset of the columns of SAi. Such strong coresets preserve the error incurred by any rank-k projection, up to a constant factor, in the `p,2 norm. The coreset of a matrix A ∈ Rd×n is usually denoted as AT , where T = D ·W ∈ Rn×t is a sampling and reweighting matrix and t is the number of columns to be included in the coreset. The sampling matrix D is a matrix with t columns where each column has only one 1 in the index of the column of A to be included in the coreset and 0 everywhere else. The reweighting matrix W is a diagonal t × t matrix with weights associated with each sample in the coreset.
Lemma 4 (Strong Coreset in `p,2 norm). Let A ∈ Rd×n, k ∈ N, p ∈ [1, 2), and , δ ∈ (0, 1). Then, in n · poly(k log n/ ) time, one can find a sampling and reweighting matrix T with O(d log d/ 2) · log(1/δ) columns such that, with probability 1− δ, for all rank-k matrices U ,
min rank-k V |UV −AT |p,2 = (1± ) min rank-k V |UV −A|p,2
AT is called a strong coreset of A.
3.3 POLYNOMIAL TIME, O(1)-APPROXIMATE BI-CRITERIA k-CSSp,2 After server i sends a strong coreset to the coordinator, the coordinator does k-CSS on a column-wise concatenation of these coresets, in the `p,2 norm rather than the `p norm. We give a polynomial time, O(1)-approximate bi-criteria k-CSSp,2 algorithm for p ∈ [1, 2).
Theorem 5 (Bicriteria O(1)-Approximation Algorithm for k-CSSp,2). Let A ∈ Rd×n and k ∈ N. There exists an algorithm that runs in (nnz(A)+d2) ·kpoly(log k) time and outputs a rescaled subset of columns U ∈ Rd×Õ(k) of A and a right factor V ∈ RÕ(k)×n for which V = minV |UV −A|p,2, such that with probability 1− o(1),
|UV −A|p,2 ≤ O(1) · min rank-k Ak |Ak −A|p,2
Our polynomial time bi-criteria k-CSSp,2 algorithm is based on that of Clarkson & Woodruff (2015). The main difference is that the algorithm of Clarkson & Woodruff (2015) outputs a subset with O(k2) columns due to the usage of `p leverage scores — we reduce the number of selected columns to Õ(k) by using `p Lewis weights (Cohen & Peng, 2015). Details are given in Appendix C.
Algorithm 1 An efficient protocol for bi-criteria k-CSSp in the column partition model Initial State: Server i holds matrix Ai ∈ Rd×ni , ∀i ∈ [s]. Coordinator: Generate a dense p-stable sketching matrix S ∈ Rk poly(log(nd))×d. Send S to all servers. Server i: Compute SAi. Let the number of samples in the coreset be t = O(kpoly(log(nd))). Construct a coreset of SAi under the `p,2 norm by applying a sampling matrix Di of size ni × t and a diagonal reweighting matrix Wi of size t× t. Let Ti = DiWi. Send SAiTi along with AiDi to the coordinator. Coordinator: Column-wise stack SAiTi to obtain SAT = [SA1T1, SA2T2, . . . , SAsTs]. Apply k-CSSp,2 on SAT to obtain the indices I of the subset of selected columns with size O(k · poly(log k)). Since Di’s are sampling matrices, the coordinator can recover the original columns of A by mapping indices I to AiDi’s. Denote the final selected subset of columns by AI . Send AI to all servers. Server i: Solve minVi |AIVi − Ai|p to obtain the right factor Vi. AI and V will be factors of a rankk · poly(log k) factorization of A, where V is the (implicit) column-wise concatenation of the Vi.
4 AN EFFICIENT PROTOCOL FOR DISTRIBUTED k-CSSp Theorem 6 (A Protocol for Distributed k-CSSp). In the column partition model, let A ∈ Rd×n be the data matrix whose columns are partitioned across s servers and suppose server i holds a subset of columns Ai ∈ Rd×ni , where n = ∑ i∈[s] ni. Then, given p ∈ [1, 2) and a desired rank k ∈ N, Algorithm 1 outputs a subset of columns AI ∈ Rd×kpoly(log(k)) in Õ(nnz(A)k + kd+ k3) time, such that with probability 1− o(1),
min V |AIV −A|p ≤ Õ(k1/p−1/2) min L⊂[n],|L|=k |ALV −A|p
Algorithm 1 uses 1 round of communication and Õ(sdk) words of communication.
Proof. Approximation Factor. In the following proof, let L ⊂ [n], |L| = k denote the best possible subset of k columns of A that gives the minimum k-CSSp cost, i.e., the cost minV |ALV − A|p achieves minimum. First, note that
min V |AIV −A|p ≤ |AIV ′ −A|p V ′ := arg min V |SAIV − SA|p,2 ≤ |SAIV ′ − SA|p By Lemma 2
= Õ(k 1 p− 1 2 )|SAIV ′ − SA|p,2 By Lemma 1, and S has k · poly(log(nd)) rows
SAI is the selected columns output from the bi-criteria O(1)-approximation k-CSSp,2 algorithm. Let (SAT )∗ denote the best rank k approximation to SAT . By Theorem 5,
Õ(k 1 p− 1 2 )|SAIV ′ − SA|p,2 ≤ Õ(k 1 p− 1 2 ) ·O(1)|(SAT )∗ − SAT |p,2
≤ Õ(k 1 p− 1 2 ) min
V |SALV − SAT |p,2
Note that SAT = [SA1T1, . . . , SAsTs] is a column-wise concatenation of all coresets of SAi, ∀i ∈ [s]. By Lemma 4,
(min V |SALV − SAT |pp,2)1/p = ( s∑ i=1 min Vi |SALVi − SAiTi|pp,2)1/p
= ( s∑ i=1 (1± )p min Vi |SALVi − SAi|pp,2)1/p = (1± )( s∑ i=1 min Vi |SALVi − SAi|pp,2)1/p
= (1± ) min V |SALV − SA|p,2
Hence,
Õ(k 1 p− 1 2 ) min
V |SALV − SAT |p,2 ≤ Õ(k
1 p− 1 2 ) min
V |SALV − SA|p,2
≤ Õ(k 1 p− 1 2 ) min
V |SALV − SA|p By Lemma 1
≤ Õ(k 1 p− 1 2 ) · log1/p(nd) min
V |ALV −A|p By Lemma 3
Therefore, we conclude: minV |AIV −A|p ≤ Õ(k 1 p− 1 2 ) minV |ALV −A|p.
Communication Cost. Sharing the dense p-stable sketching matrix S with all servers costsO(sdk · poly(log(nd))) communication (this can be removed with a shared random seed). Sending all coresets SAiTi (∀i ∈ [s]) and the corresponding columns AiDi to the coordinator costs Õ(sdk) communication, since each coreset contains only Õ(k) columns. Finally, the coordinator needs Õ(sdk) words of communication to send the Õ(k) selected columns to each server. Therefore, the overall communication cost is Õ(sdk), suppressing a 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 SAi’s takes O(nnz(A)k · poly(log(nd))) time. By Lemma 4, computing a single coreset SAiTi on server i takes time nipoly(k log ni) 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 ) + k2poly(log nd)) · kpoly(log k) ≤ k3poly(log(knd)) to find the set of selected columns. Since the number of selected columns is O(kpoly(log k)), it then takes the protocol O(kpoly(log k)) time to map the indices of the output columns from k-CSSp,2 to recover the original columns AI . Therefore, the overall running time for the protocol to find the subset of columns AI is Õ((nnz(A)k + kd + k3), suppressing a low degree polynomial dependency on log(knd). After the servers receive AI , it is possible to solve minVi |AIVi−Ai|p in Õ(nnz(AI)) + poly(d log n) time , ∀i ∈ [s] due to Wang & Woodruff (2019); Yang et al. (2018).
5 GREEDY k-CSS1,2 We propose a greedy algorithm, shown in Algorithm 2, for k-CSS1,2, which can be used in the place of the algorithm described in Theorem 5. 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-CSS2 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 et al. (2016), where in each iteration, rather than considering all columns of A as candidate additional columns of AT , we only sample a subset of the columns of A of size O(n log(1/δ)k ), and pick the column among those that improves the objective the most.
Theorem 7. Let A ∈ Rd×n be the data matrix and k ∈ N be the desired rank. Let AS be the best possible subset of k columns, i.e., AS = arg minAS minV |ASV − A|1,2. Let σ be the minimum non-zero singular value of the matrix B of normalized columns of AS , (the j-th column of B is B∗j = (AS)∗j/|(AS)∗j |2). Then, if T ⊂ [n] is the subset of columns selected by Greedy k-CSS1,2, the following holds with |T | = Ω( kσ2 2 ):
min V |ATV −A|1,2 ≤ (1− ) min S⊂[n],|S|=k,V ∈Rk×n |ASV −A|1,2 + |A|1,2
Similarly, if T ⊂ [n] is the subset of columns selected by Lazier-than-lazy Greedy k-CSS1,2, the following holds with |T | = Ω( kσ2 2 ) and δ = :
E[min V |ATV −A|1,2] ≤ (1− ) min S⊂[n],|S|=k,V ∈Rk×n |ASV −A|1,2 + |A|1,2
Algorithm 2 Lazier-than-lazy Greedy k-CSS1,2. This version of the greedy algorithm is based on Section 5.2 of Altschuler et al. (2016).
If we let |T | = Ω( kσ2 2 ), then the overall running time of Algorithm 2 is O( n log(1/ ) σ2 2 · F ), where F is the running time needed to evaluate minV |AT∪jV − A|1,2 for a fixed j ∈ T . We can get F = O( dk 2
σ4 4 + ndk σ2 2 ) by taking A † T∪j .
Since the error upper bound for greedy k-CSSp,2 depends on |A|1,2, it is not directly comparable to the error upper bound for the proposed k-CSSp,2 from Subsection 3.3, which achieves a multiplicative O(1)-approximation to the best rank-k approximation. We empirically compare the two versions of k-CSSp,2 for p = 1 in Section 6.

heading: 6 EXPERIMENTS
text: We implement our protocol for distributed k-CSSp in Algorithm 1, setting p = 1, which enables us to compare two subroutines on the coordinator: Regular k-CSS1,2 from Section 3.3 and Greedy k-CSS1,2 from Section 5. We compare our k-CSS1 protocol against a commonly applied baseline for `p low rank approximation (used by Song et al. (2019a); Chierichetti et al. (2017)): rank-k Singular Value Decomposition (SVD).
Datasets. We demonstrate the benefits of our k-CSS1 protocol on one synthetic data and one real-world application. We present a summary of the datasets, along with the number of servers s, the column distribution across servers and the rank k we consider for each dataset in Table 1. The synthetic dataset constructs a data matrix M ∈ R(k+n)×(k+n) such that the top left k × k submatrix is the identity matrix multiplied by n 3 2 , and the bottom right n × n submatrix has all 1’s. The optimal rank-k left factor consists of one of the last n columns along with k − 1 of the first k columns, incurring an error of n 3 2 in the `1 norm and an error n3 in the squared `2 norm. SVD, however, will not cover any of the last n columns, and thus will get an error of n2 in both the `1 and squared `2 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 a real-world application, term-document clustering, where k-CSS2 algorithm is previously applied (Mahoney & Drineas, 2009). TechTC3 contains 139 documents processed in bag-of-words representation with a dictionary of 18446 words. Such representation naturally results in a sparse matrix. k-CSSp is used to select the top k most representative words.
Hyperparameters. Experiment hyperparameters are summarized in Table 2. 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-CSS1,2. We denote the number of rows in the sparse embedding matrix of Õ(k) rows by sketch size, and the number of non-zero entries in each column of the sparse embedding matrix by sparsity.
3http://gabrilovich.com/resources/data/techtc/techtc300/techtc300.html
Setup. We conduct 15 trials for each experiment setting, and report the average `1 error, i.e. |Âk − A|1 for a Õ(k)-rank matrix Âk output by the algorithm that approximates the original data matrix A ∈ Rd×n, across all trials. However, we notice that computing the `1 error for our protocol involves computing an `1 regression on all columns of the data matrix, which takes prohibitively long time on the real-world dataset TechTC. Therefore, for TechTC, we proceed with an approximate computation as follows: for the output subset of columns AI , first i.i.d. sample 30% columns from A, denoted as AJ , with probability proportional to |A∗j |1, ∀j ∈ [n]; then compute a lower bound for the `1 error L = ∑ j∈J minv |AIv −A∗j |1; after that, sample 20|A|1 L columns from A, denoted as AK , with probability again proportional to |A∗j |1, ∀j ∈ [n]; finally, the `1 error is computed as L |A|1 · ∑ k∈K 1 20|A∗k|1 minv |AIv −A∗k|1. The approximated `1 error can be shown as an unbiased estimation of minV |AIV −A|1 with small variance by Hoeffding’s bound.
Results. We present our empirical results in Figure 2. The distributed protocol performs better using GREEDY k-CSS1,2 than REGULAR k-CSS1,2 on both datasets, and in other settings we include in the supplementary material.
We also conducted experiments on three different datasets, bcsstk13, isolet, and caltech-101, extensively comparing the approximation error as well as the running time for different algorithms. A comprehensive reporting of our results is given in Appendix G.

heading: 7 CONCLUSION
text: In this work, we give the first nearly-optimal communication and number of rounds protocol for distributed k-CSSp (1 ≤ p < 2) in the column partition model, which achieves Õ(k1/p−1/2)approximation to the best possible subset of columns, with Õ(sdk) communication cost, 1 rounds and polynomial time. To achieve a good approximation factor, we use dense p-stable sketching and work with the `p,2 norm, which enables us to use an efficient construction of strong coresets and an O(1)-approximation bi-criteria k-CSSp,2 algorithm. We further propose a greedy algorithm for kCSS1,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-CSS1,2 and regular k-CSS1,2. Our results empirically show that greedy k-CSS1,2 gives substantial improvements over regular k-CSS1,2 on real-world datasets. For future works, it is not known whether a Õ(k1/p−1/2)-approximation factor to the best possible subset of columns, is optimal for distributed k-CSSp (1 ≤ p < 2).

heading: B PRELIMINARIES FOR OUR UPPER BOUND PROOFS

heading: B.1 NORMS
text: Lemma 1. (Norm Relationships) For a matrix A ∈ Rd×n, |A|p,2 ≤ |A|p and |A|p ≤ d 1 p − 1
2 |A|p,2, where 1 ≤ p < 2.
Proof. Let x ∈ Rd. For 0 < p < r, |x|r ≤ |x|p ≤ d 1 p − 1 r |x|r
Let r = 2. Then we have |x|2 ≤ |x|p ≤ d 1 p − 1 2 |x|2
Note that |A|p,2 = (∑ j |A∗j | p 2 ) 1 p and |A|p = (∑ j |A∗j | p p ) 1 p .
Therefore,
|A|p,2 = (∑ j |A∗j |p2 ) 1 p ≤ (∑ j |A∗j |pp ) 1 p = |A|p
and
|A|p = (∑ j |A∗j |pp ) 1 p ≤ d 1 p − 1 2 (∑ j |A∗j |p2 ) 1 p = d 1 p − 1 2 |A|p,2
B.2 LOWER BOUND ON THE COST - NO CONTRACTION WHEN APPLYING A p-STABLE SKETCH
We show a lower bound on the approximation error for a sketched subset of columns, |SATV − SA|p, in terms of |ATV −A|p. The lower bound holds simultaneously for any arbitrary subset AT of chosen columns, and for any arbitrary right factor V .
We begin the proof by first showing that applying a dense p-stable sketch to a vector will not shrink its p-norm. This is done in Lemma 2.1. We further observe that although p-stable random variables are heavy-tailed, we can still bound their tail probabilities by applying Lemma 9 from Meng & Mahoney (2013). We note this in Lemma 2.2. Note that the Xi’s do not need to be independent in this lemma.
Equipped with Lemma 2.1, Lemma 2.2 and a net argument, we can now establish a lower bound on |SATV − SA|p. We first show in Lemma 2.3 that, with high probability, for any arbitrarily selected subset AT of columns and for an arbitrary column A∗j , the error incurred to fit SA∗j using the columns of SAT is no less than the error incurred to fit A∗j using the columns of AT . We then apply a union bound over all subsets T ⊂ [n] and columns j ∈ [n] to conclude our lower bound in Lemma 2.
Lemma 2.1. (No Contraction of p-stable Sketch) Given a matrix S ∈ Rt×m whose entries are i.i.d. p-stable random variables rescaled by Θ ( 1
t 1 p
) , where 1 ≤ p < 2, for any fixed y ∈ Rm, with probability 1− 1
et , the
following holds:
|Sy|p ≥ |y|p
Proof. By p-stability, we have |Sy|pp = ∑t i=1 ( |y|p |Zi|
t 1 p
)p , where |Zi| ≥ 0 are half p-stable random variables.
Since Pr[|Zi| = Ω(1)] > 12 , by applying a Chernoff bound (to the indicators 1|Zi|≥C for a sufficiently small constant C), we have ∑t i=1 |Zi| p = Ω(t) with probability 1− 1 et . Therefore, with probability 1− 1 et
, |Sy|p ≥ |y|p.
Lemma 2.2. (Upper Tail Inequality for p-stable Distributions) Let p ∈ (1, 2), and m > 3. For i ∈ [m], let Xi be a standard p-stable random variable, and let γi > 0 and γ = ∑m i=1 γi. Let X = ∑m i=1 γi|Xi|
p. Then, for any t ≥ 1, Pr[X ≥ tαpγ] ≤ 2 log(mt)t , where αp > 0 is a constant that is at most 2 p−1.
Proof. Lemma 9 from Meng & Mahoney (2013) for p ∈ (1, 2).
Lemma 2.3. (No Contraction for All Sketched Subsets and Columns) Let A ∈ Rd×n, and k ∈ N. Let t = k · poly(lognd), and let S ∈ Rt×d be a matrix whose entries are i.i.d. standard p-stable random variables, rescaled by Θ(1/t 1 p ). Finally, let m = k · poly(log k). Then, with probability 1− o(1), for all T ⊂ [n] with |T | = m, for all j ∈ [n], and for all y ∈ R|T |,
|AT y −A∗j |p ≤ |S(AT y −A∗j)|p
Proof. Step 1: We first extend Lemma 2.1 and argue that applying a p-stable sketching matrix S ∈ Rt×n will not shrink the norm |Sy|p ≥ |y|p (1 ≤ p < 2) simultaneously for all y in the column span of [AT , Aj ] =: AT,j , by a net argument.
In order to bound the p-norm of sketched vectors y in a net, we begin by showing that with high probability all entries in S are bounded. Let D = poly(mt). Consider the following two cases:
Case 1: p = 1. The entries of the 1-stable sketching matrix Sij are standard Cauchy random variables. Consider half-Cauchy random variables Xi,j = |Si,j |. The cumulative distribution function of half Cauchy random variables x is F (x) = ∫ x 0 2 π(t2+1) dt = 1 − Θ( 1 x
). Thus, for any i ∈ [t] and j ∈ [m], Pr[|Sij | ≤ D] = 1−Θ( 1
D ).
Case 2: p ∈ (1, 2). We apply the upper tail bound of p-stable random variables in Lemma 2.2. For any fixed i ∈ [t] and j ∈ [m], Pr[|Sij |p ≤ Dp] ≥ 1−Θ( log(t)Dp ), which implies Pr[|Sij | ≤ D] ≥ 1−Θ( 1 D ).
Therefore, for p ∈ [1, 2), by a union bound over all entries in S, if we define the event E1 to mean that for all i ∈ [t] and j ∈ [m], we simultaneously have |Sij | ≤ D, then Pr[E1] ≥ 1 − Θ(mtD ) by a union bound. The event E1 occurring implies that for any y ∈ Rd, since all entries in S are rescaled by O(1/t1/p),
|Sy|p ≤ |y|pt1/p|S|∞ ≤ D|y|p
Consider the unit `p ball B = {y ∈ Rd : |y|p = 1, ∃z ∈ Rm s.t. y = AT,jz} in the column span of AT,j . A subset N ⊂ B is a γ-net for B if for all y ∈ B there exists some u ∈ N such that |y − u|p ≤ γ, for some distance γ > 0. There exists such a net N for B of size |N | = ( 1
γ )O(m) by a standard greedy
construction, since the column span of AT,j has dimension at most m + 1. We choose γ = 1m2D , and thus |N | ≤ (m2D)O(m) = DO(m logm).
By applying Lemma 2.1, and a union bound over all vectors y ∈ N , we have that event E2: for all y ∈ N simultaneously, |Sy|p ≥ |y|p — E2 has probability at least 1− D O(m logm)
et .
Consider an arbitrary unit vector x ∈ B. There exists some y ∈ N such that |x−y|p ≤ γ = 1m2D . Conditioning on both events E1 and event E2, we have the following with probability 1−Θ(mtD )− DO(m logm) et :
|Sx|p ≥ |Sy|p − |S(x− y)|p Triangle Inequality ≥ |y|p − |S(x− y)|p By event E2 ≥ |y|p −D|(x− y)|p Implication of event E1 ≥ |y|p −Dγ By |x− y|p ≤ γ
= |y|p −O( 1
m2 )
= |x|p −O( 1
m2 ) |x|p = |y|p = 1
For a sufficiently large m, O( 1 m2 ) is at most 1 2 , and thus |x|p 2 = 1 2 ≥ 1 m2 . This implies |Sx|p ≥ |x|p − |x|p2 = |x|p 2 . We can rescale S by a factor of 2 so that |Sx|p ≥ |x|p.
We have shown that |Sy|p ≥ |y|p holds simultaneously for all unit vectors y in the column span of AT,j . By linearity, we conclude that |Sy|p ≥ |y|p (1 ≤ p < 2) holds simultaneously for all y in the column span of AT,j .
Step 2: Next, we apply a union bound over all possible subsets T ⊂ [n] of chosen columns from A and all possible single columnsA∗j for j ∈ [n], to argue that |S(AT y−A∗j)|p ≥ |AT y−A∗j |p holds for all y ∈ R|T | and all T ⊂ [n], j ⊂ [n] with high probability. Note that |T | = m = O(k · poly(log k)).
From Step 1, we have shown that event E2 fails with probability D O(m logm)
et . Thus event E2 fails over all
possible subsets T and all possible single columns A∗j with probability at most
DO(m logm)
et ·
( d
m− 1
) · d ≤ D O(m logm)
et · dO(m) · d
The above failure probability is o(1) as long as t = Θ ( log ( DO(m logm) · dO(m) · d) )) = Θ(O(m logm) log(D) +O(m) log d)
= Θ(O(kpoly(log k)) log(mt) +O(kpoly(log k) log(d)))
Thus it suffices to have t = kpoly(lognd) to have failure probability at most 1 eO(t) .
Now we condition on a single global event E1 and that event E2 holds for all possible T and A∗j . We conclude that with probability 1−Θ(mt
D )− 1 eO(t) , the following holds simultaneously for all T ⊂ [n] and for all A∗j for
which j ∈ [n]:
|AT y −A∗j |p ≤ |S(AT y −A∗j)|p.
Lemma 2. (Lower Bound for Sketched Error) Let A ∈ Rd×n and k ∈ N. Let t = k · poly(log(nd)), and let S ∈ Rt×d be a matrix whose entries are i.i.d. standard p-stable random variables, rescaled by Θ(1/t 1 p ). Then, with probability 1− o(1), for all T ⊂ [n] with |T | = k · poly(log k) and for all V ∈ R|T |×n,
|ATV −A|p ≤ |SATV − SA|p
Proof. Let yj denote the j-th column of V , where j ∈ [n]. By applying Lemma 2.3, and a union bound over all columns of V , for m = |T | = k · poly(log k), t = k · poly(log(nd)) and D = poly(mt) = poly(k log(nd)), the following holds with probability 1−Θ(mt
D )− n eO(t) = 1− o(1),
|ATV −A|p = ( n∑ j=1 |AT yj −Aj |pp) 1 p
≤ ( n∑ j=1 |S(AT yj −Aj)|pp) 1 p
= |SATV − SA|p

heading: B.3 UPPER BOUND ON THE COST
text: We show an upper bound on the approximation error of k-CSSp on a sketched subset of columns, |SATV − SAT |p, which holds for a fixed subset AT of columns and for the minimizing right factor V = arg minV |SATV − SA|p for that subset of columns.
We first adapt Lemma E.17 from Song et al. (2017) to establish an upper bound on the error |SATV − SA|p for any fixed V in Lemma 3.1. We then apply Lemma 3.1 to the minimizer V to conclude the upper bound in Lemma 3.
Lemma 3.1. (An Upper Bound on Norm of A Sketched Matrix) GivenA ∈ Rn×d and p ∈ [1, 2), andU ∈ Rn×k and V ∈ Rk×d, if S ∈ Rt×n is a dense p-stable matrix, whose entries are rescaled by Θ ( 1
t 1 p
) , then with
probability at least 1−O(1),
|SUV − SA|pp ≤ O(log(td))|UV −A|pp
Here, the failure probability O(1) can be arbitrarily small.
Proof. Lemma E.17 from Song et al. (2017).
Lemma 3. (Upper Bound on Sketched Error) Let A ∈ Rd×n and k ∈ N. Let t = k · poly(log(nd)), and let S ∈ Rt×d be a matrix whose entries are i.i.d. standard p-stable random variables, rescaled by Θ(1/t 1 p ). Then, for a fixed subset T ⊂ [n] of columns with |T | = k · poly(log k), with probability 1−O(1), we have
min V |SATV − SA|p ≤ min V O(log1/p(nd))|ATV −A|p
Here, the failure probability can be an arbitrarily small constant.
Proof. Let X∗1 = arg minX |SATX − SA|p and X∗2 = arg minX |ATX −A|p. By Lemma 3.1,
|SATX∗1 − SA|pp ≤ |SATX∗2 − SA|pp ≤ O(log(kpoly(logn)d))|ATX∗2 −A|pp ≤ O(log(nd))|ATX∗2 −A|pp
Therefore, min X |SATX − SA|p ≤ min X O(log1/p(nd))|ATX −A|p
.
B.4 STRONG CORESETS IN THE `p,2 NORM
Lemma 4 (Strong Coreset in `p,2 norm). Let A ∈ Rd×n, k ∈ N, p ∈ [1, 2), and , δ ∈ (0, 1). Then, in n · poly(k logn/ ) time, one can find a sampling and reweighting matrix T with O(d log d/ 2) · log(1/δ) columns such that, with probability 1− δ, for all rank-k matrices U ,
min rank-k V |UV −AT |p,2 = (1± ) min rank-k V |UV −A|p,2
AT is called a strong coreset of A.
Proof. We can obtain T with O(d(log d)/ 2) columns using the strong coreset construction from Lemma 16 in Sohler & Woodruff (2018). Note that the coreset construction for k-subspace approximation in Sohler & Woodruff (2018) aims at removing a dependence on d in the coreset size. The algorithm first finds a poly(k)dimensional subspace S by running a dimensionality reduction algorithm and constructs coresets in the lower dimensional subspace S, resulting in a coreset size of poly(k/ ). But in our case, we do not want our coreset size to have a polynomial dependency on k while a linear dependency on d suffices. Thus, instead of running their dimensionality reduction algorithm to find such a subspace S to project A to, we directly use the the column span of the input matrix A as the subspace S appended with a row of zeros to construct the B in Lemma 16 of Sohler & Woodruff (2018). Note that in the original algorithm, the appended column encodes the distances from the input matrix A to the subspace S, but in our case it is just all 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 Õ(k), suppressing a logarithmic dependence on k, 1
, 1 δ , using an additional O(nd) time, by combining Corollary 5.16 of Huang & Vishnoi (2020) and the importance sampling scheme from Stage 2 of Algorithm 1 in Huang & Vishnoi (2020). Furthermore, though not explicitly stated in Lemma 16 of Sohler & Woodruff (2018), the coreset size has a log( 1
δ ) dependence on the failure probability δ due to the importance sampling (Sohler & Woodruff, 2018;
Huang & Vishnoi, 2020).
Note that an accuracy of (1± ) is desired for our extension to the streaming model, described in Appendix E, and for both Algorithm 1 and Algorithm 7, we perform a union bound over strong coreset constructions, for which a log(1/δ) dependence on δ is sufficient.
C POLYNOMIAL TIME, O(1)-APPROXIMATE BI-CRITERIA k-CSSp,2 We give a detailed analysis on the polynomial time, O(1)-approximate k-CSSp,2 algorithm presented in Algorithm 4, which is based on ideas in Clarkson & Woodruff (2015). We first use a sparse embedding matrix S to obtain an O(1)-approximate left factor. We then use `p-Lewis weight sampling Cohen & Peng (2015) to select a subset of columns.
Algorithm 4 polynomial time, O(1)-approximation for k-CSSp,2 (1 ≤ p < 2) Input: The data matrix A ∈ Rd×n, rank k ∈ N Output: The left factor U ∈ Rd×Õ(k), the right factor V ∈ RÕ(k)×n such that |UV − A|p,2 ≤ O(1) minrank-kAk |Ak −A|p,2 S ← Õ(k)× d sparse embedding matrix, with sparsity s = poly(log k). S′ ← n× Õ(k) sampling matrix, each column of which is a standard basis vector chosen randomly according to the `p Lewis weights of columns of SA. return U ← AS′, V ← (AS′)†A {† denotes the Moore-Penrose pseudoinverse.}

heading: C.1 SPARSE EMBEDDING MATRICES
text: The sparse embedding matrix S ∈ RÕ(k)×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 ± 1√
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 σ be the hash function that determines the sign ± of each non-zero entry.
Applying the sparse embedding matrix S to A enables us to obtain a rank-k right factor that is at most a factor of O(1) worse than the best rank-k approximation error in the `p,2 norm. We adapt Theorem 32 from Clarkson & Woodruff (2015) to show this in Theorem 5.5. Notice that in Theorem 32 of Clarkson & Woodruff (2015), the number of rows required for S is O(k2), but this can be reduced to Õ(k) through a different choice of hyperparameters when constructing the sparse embedding matrix S.
We note two choices of hyperparameters, i.e., the number m of rows and sparsity s, 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 Theorem 5.1. We replace the hyperparameters from Theorem 5.2 and show in Lemma 5.3 that Õ(k) rows of S suffice to preserve certain desired properties. We then combine Lemma 5.3 and Lemma 5.4 adapted from Clarkson & Woodruff (2015), to conclude our result in Theorem 5.5, following the analysis from Clarkson & Woodruff (2015).
Theorem 5.1. (Theorem 3 from Nelson & Nguyen (2013)) For a sparse embedding matrix S ∈ Rm×n with sparsity s = 1 and a data matrix U ∈ Rn×d, let ∈ (0, 1). With probability at least 1− δ all singular values of SU are (1± ) as long as m ≥ δ−1(d2 + d)/(2 − 2)2. For the hash functions used to construct S, σ is 4-wise independent and h is pairwise independent.
Theorem 5.2. (Theorem 9 from Nelson & Nguyen (2013)) For a sparse embedding matrix S ∈ Rm×n with sparsity s = Θ(log3(d/δ)/ ) and a data matrix U ∈ Rn×d, let ∈ (0, 1). With probability at least 1− δ all singular values of SU are (1± ) as long as m = Ω(d log8(d/δ)/ 2). For the hash functions used to construct S, we have that σ, h are both Ω(log(d/δ))-wise independent.
Lemma 5.3. Let C be a constraint set and A ∈ Rn×d, B ∈ Rn×d ′ be two arbitrary matrices. For a sparse embedding matrix S ∈ Rm×n, there is m = O( d log8( d p+1 )
2(p+1) ), such that with constant probability, the following
holds for all X ∈ Rd×d ′
i) |S(AX −B)|p,2 ≥ (1− )|AX −B|p,2 ii) |S(AX∗ −B)|p,2 ≤ (1 + )|AX∗ −B|p,2, where X∗ = arg min
X∈C |AX −B|p,2
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.e., sparsity s 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 = 1, but requires m = O(d2) rows. We replace the construction by Theorem 5.2, where we pick δ = p+1. Now the sparsity s is larger but this construction reduces the number of rows required to m = Õ(d). Since both choices of hyperparameters to construct S result in bounded (1± ) singular values of SU for any data matrix U , the rest of the proof follows.
Lemma 5.4. Consider a data matrix A ∈ Rn×d. Let the best rank-k matrix in the `p,2 norm be Ak = arg minrank-kAk |Ak − A|p,2. For R ∈ R
d×m, if RT satisfies both of the following two conditions for all X ∈ Rn×n:
i) |RT (ATkX −AT )|p,2 ≥ (1− )|ATkX −AT |p,2 ii) |RT (ATkX∗ −AT )|p,2 ≤ (1 + )|ATkX∗ −AT |p,2, where X∗ = arg min
X |ATkX −AT |p,2
then
min rank-kX
|XRTATk −AT |pp,2 ≤ (1 + 3 )|A T k −AT |pp,2
Proof. Lemma 31 from Clarkson & Woodruff (2015).
Theorem 5.5. (`p,2-Low Rank Approximation) Let the data matrix be A ∈ Rd×n and k ∈ N be the desired rank. Let S ∈ Rm×d be a sparse embedding matrix with m = O(kpoly(log k)poly( 1 )) rows, and sparsity s = poly(log k). Then, the following holds with constant probability,
min rank-kX |XSA−A|p,2 ≤ (1 + 3 ) min rank-kAk |Ak −A|p,2
Proof. The proof is the same as the proof of Theorem 32 in Clarkson & Woodruff (2015), except that we adapt a different construction of the sparse embedding matrix S, which reduces the number of rows from O(k2) to Õ(k) with increased sparsity s.
Consider Ak = arg minrank-kAk |Ak − A|p,2. Let Vk be a basis for the column space of Ak. By applying Lemma 5.3 and Lemma 5.4 on the basis Vk, we conclude the above theorem by setting the number m of rows
to m = O ( k log8( k 2 )
4
) , and sparsity s = poly(log k) in the sparse embedding matrix S.

heading: C.2 LEWIS WEIGHTS SAMPLING
text: The `p Lewis weight is an inherent property of a matrix. By Cohen & Peng (2015), the unique set of `p Lewis weights w for a matrix is defined as follows: for the i-th column of a matrix M ∈ Rn×d, M∗j , is defined as w
2/p i = M T ∗j(MW 1−2/pMT )−1M∗j , where W is the diagonal matrix with Wii = wi, ∀i ∈ [d].
We first note that sampling by Lewis weights provides a subspace embedding in the `p-norm in Theorem 5.6. We further apply a randomized version of Dvoretzky’s Theorem in Theorem 5.7, which allows the embedding from the `p norm into a low-dimensional Euclidean subspace with very small distortion and thus enables us to switch between the `p norm and the `p,2 norm. Based on Theorem 5.6 and Theorem 5.7, we show in Theorem 5.8 that Lewis weights sampling provides a good subset of columns, on which the analysis of k-CSSp,2 is based.
Theorem 5.6. (`p-Lewis Weights Subspace Embedding) Let A ∈ Rn×d and t = Õ(d). For 1 ≤ p < 2, there exists a distribution (λ1, λ2, . . . , λn) on the rows of A such that if we generate a matrix S with t rows, each chosen independently as the ith standard basis vector times 1
(rλi) 1 p
with probability λi, then with probability
1− o(1), the following holds for all x ∈ Rd,
Ω(1)|Ax|p ≤ |SAx|p ≤ O(1)|Ax|p
Proof. Theorem 7.1 from Cohen & Peng (2015)
Theorem 5.7. (Randomized Dvoretzky’s Theorem) Let n ∈ N, and ∈ (0, 1). Let r = n 2 . Let G ∈ Rr×n be a random matrix whose entries are i.i.d. standard Gaussian random variables, rescaled by 1√
r . For r = n 2 , the
following holds with probability 1− e−Θ(n), for all y ∈ Rn,
|Gy|p = (1± )|y|2
Proof. This follows from Theorem 1.2 from Paouris et al. (2017).
Theorem 5.8. (Subset of Columns by Lewis Weights Sampling) Let A ∈ Rd×n. Let S ∈ Rm×d be a sparse embedding matrix, with m = O(k · poly(log k)poly( 1 )). Further, let S′ ∈ Rn×t be a sampling matrix whose columns are random standard basis vectors generated according to the `p Lewis weights of columns of SA, with t = k · poly(log k). Then, for X̂ = arg minrank-kX |XSAS′ −AS′|p,2, the following holds with constant probability,
|X̂SA−A|p,2 ≤ Θ(1) min rank-kAk |Ak −A|p,2
Proof. Let X∗ = arg minrank-kX∗ |X∗SA−A|p,2. By the triangle inequality,
|X̂SA−A|p,2 ≤ |X∗SA− X̂SA|p,2 + |X∗SA−A|p,2
Let us consider bounding |X∗SA− X̂SA|p,2.
By Lemma D.28 and Lemma D.29 from Song et al. (2017), for any column sampling matrix S and for any fixed matrix Y , it can be shown that E[|Y S|pp] = |Y S|pp. In our case, since S′ is a sampling matrix, we have E[|Y S′|pp] = |Y S′|pp for any fixed matrix Y .
Now let G ∈ RΘ(d)×d be a rescaled random matrix whose entries are i.i.d. standard Gaussian random variables as in Theorem 5.7. We apply Theorem 5.7 to transform between the `p space and the Euclidean space. Since transformation of both directions can be done with very small distortion, we obtain a Θ(1) approximation. With constant probability, we have
|X∗SA− X̂SA|p,2 = Θ(1)|G(X∗ − X̂)SA|p By Theorem 5.7
= Θ(1)|G(X∗ − X̂)SAS′|p By Theorem 5.6
= Θ(1)|(X∗ − X̂)SAS′|p,2 By Theorem 5.7 ≤ Θ(1) ( |X∗SAS′ −AS′|p,2 + |X̂SAS′ −AS′|p,2 ) Triangle Inequality
≤ Θ(1)|X∗SAS′ −AS′|p,2 Since X̂ = arg min rank-kX |XSAS′ −AS′|p,2 = Θ(1)|G(X∗SA−A)S′|p By Theorem 5.7 ≤ Θ(1)|G(X∗SA−A)|p By Markov Bound on E[|Y S′|pp] = |Y S′| = Θ(1)|X∗SA−A|p,2 By Theorem 5.7
Therefore,
|X̂SA−A|p,2 ≤ |X∗SA− X̂SA|p,2 + |X∗SA−A|p,2 ≤ Θ(1)|X∗SA−A|p,2 ≤ Θ(1) min
rank-kAk |A−Ak|p,2 By Theorem 5.5
C.3 ANALYSIS FOR k-CSSp,2 We are now ready to show an O(1) approximation factor and polynomial running time for the bi-criteria k-CSSp,2 presented in Algorithm 4, given in Theorem 5.
Theorem 5 (Bicriteria O(1)-Approximation Algorithm for k-CSSp,2). Let A ∈ Rd×n and k ∈ N. There exists an algorithm that runs in (nnz(A) + d2) · kpoly(log k) time and outputs a rescaled subset of columns U ∈ Rd×Õ(k) of A and a right factor V ∈ RÕ(k)×n for which V = minV |UV − A|p,2, such that with probability 1− o(1),
|UV −A|p,2 ≤ O(1) · min rank-k Ak |Ak −A|p,2
Proof. Approximation Factor First notice that the minimizer X̂ of |X̂SAS′−AS′|p,2 has to be in the column span of AS′. Thus we can write X̂ = (AS′)Y for some matrix Y . By Theorem 5.8,
|X̂SA−A|p,2 = |(AS′)Y SA−A|p,2 ≤ Θ(1) min rank-kAk |A−Ak|p,2
We denote Y SA = V . We take the left factor U = AS′ and solving for minV |UV −A|p,2 will give us a Θ(1) approximation to minrank-kAk |A−Ak|p,2. A good minimizer for the right factor V in the Euclidean space is V = (AS′)†A. This concludes our result. Notice that since S′ is a sampling matrix with Õ(k) columns, we get a rank-k left factor U as a subset of columns of A as desired.
Running time First notice that S is a sparse embedding matrix with poly(log k) non-zero entries. Thus computing SA takes time nnz(A) · poly(log k). By Cohen & Peng (2015), computing the Lewis weights of SA takes time nnz(SA) + poly(k) ≤ nnz(A)poly(log k) + poly(k), and computing the output left factor U = AS′ takes time nnz(A). Computing (AS′)† takes time d2 · kpoly(log k). Computing the right factor V = (AS′)†A takes nnz(A)kpoly(log k). Therefore, the overall running time is (nnz(A) + d2) · kpoly(log k).
D ANALYSIS FOR GREEDY k-CSS1,2 We propose a greedy algorithm for selecting columns in k-CSS1,2 presented in Algorithm ??. We give a detailed analysis on the first additive approximation compared to the best possible subset of columns for Greedy k-CSS1,2. Our analysis is inspired by the analysis of Greedy k-CSS2 for the Frobenius norm in Altschuler et al. (2016). We then describe how the running time of Algorithm ?? can be improved to Õ( n
σ · F ), where F is the running time
required to evaluate minV |AT∪jV −A|1,2 for a fixed j ∈ T (note that this running time is O(nk2 + knd) if we evaluate this by computing the pseudo-inverse of AT∪j). This improvement in the running time is obtained by randomly sampling candidate columns from T and adding the best of these randomly sampled columns to AT , rather than trying all columns of T — this method was previously used in Altschuler et al. (2016) for Frobenius norm column subset selection, where it is referred to as “Lazier-than-lazy Greedy,” and the general approach was first introduced in Mirzasoleiman et al. (2015). This version of the greedy algorithm is shown in Algorithm 2, and we show both versions of the greedy algorithm below for convenience.
Algorithm 5 Greedy k-CSS1,2. Here we denote the set of selected columns T from A by AT and the set of unselected columns by AT .
Input: The data matrix A ∈ Rd×n, the number of iterations r ≤ n. Output: A subset of columns AT from A, where |T | = r. AT ← ∅ for i = 1 to r do
Column j∗ ← arg minj∈AT (minV |AT∪jV −A|1,2) AT ← AT∪j∗
end for
Algorithm 6 Lazier-than-lazy Greedy k-CSS1,2. This version of the greedy algorithm is based on Section 5.2 of Altschuler et al. (2016).
Input: The data matrix A ∈ Rd×n, the number of iterations r ≤ n, a parameter δ ∈ (0, 1). Output: A subset of columns AT from A, where |T | = r. AT ← ∅ for i = 1 to r do T ← A subset of n log(1/δ)k columns of A, each selected uniformly at random (excluding the columns whose indices are in T ) Column j∗ ← arg minj∈AT (minV |AT∪jV −A|1,2) AT ← AT∪j∗ end for
First we analyze Algorithm 4, then show how this analysis can be extended to Algorithm 6. We first show in Lemma 7.2 an improvement of 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 Inequality, following the analysis in Altschuler et al. (2016). Finally, we conclude our analysis for Greedy k-CSS1,2 in Theorem 7.
Notation Consider the input matrix A ∈ Rd×n (n d). Let B be the matrix of normalized columns of A, where the j-th column of B is B∗j = A∗j/|A∗j |2. Let πT : Rd → Rd be the projection onto the column span of AT or equivalently BT . Let σmin(M) denote the minimum singular value of some matrix M .
To aid our analysis, we define a utility function as follows, inspired by Altschuler et al. (2016). For a subset T ⊂ [n] and a matrix M ∈ Rd×t (or a vector M ∈ Rd),
ΦM (T ) = |M |1,2 − |M − πTM |1,2 = t∑ i=1 ( |M∗i|2 − |M∗i − πTM∗i|1,2 ) = t∑ i=1 ΦM∗i(T )
Observe that as the number of columns selected and added to T increases, we get a more accurate estimation of M and thus the approximation error |M − πTM |1,2 decreases, which results in an increase in the utility function ΦM (T ).
Lemma 7.1. Let S, T ⊂ [n] be two sets of column indices, with |πSu|2 ≥ |πTu|2 for some vector u ∈ Rd. Then,
k∑ i=1 ( |πT ′iu| 2 2 − |πTu|22 ) ≥ σmin(BS)2 (|πSu|22 − |πTu|22)2 4|πSu|22
Proof. Lemma 2 from Altschuler et al. (2016), except that we replace the condition for S and T , i.e., Φu(S) ≥ Φu(T ) in Altschuler et al. (2016) with |πSu|2 ≥ |πTu|2. The two conditions are equivalent, since
Φu(S) ≥ Φu(T ) ⇔ |u− πSu|2 ≤ |u− πTu|2
⇔ |u|22 − |πSu|22 ≤ |u|22 − |πTu|22 ⇔ |πSu|2 ≥ |πTu|2
Lemma 7.2. (Utility Improvement by Projecting a Single Vector) Let S, T ⊂ [n] be two sets of column indices, with Φu(S) ≥ Φu(T ) for some vector u ∈ Rd. Let k = |S|. For i ∈ [k], let T ′i = T ∪ {i}. Then,
k∑ i=1 ( Φu(T ′ i )− Φu(T ) ) ≥ σmin(BS)2 (Φu(S)− Φu(T ))3 16Φu(S)2
Proof. We define a function for convenience in the analysis g : (−∞, |u|22] → R≥0 by g(x) = √ |u|22 − x. Note that |g′(x)| = 1 2 √ |u|22−x .
k∑ i=1 ( Φu(T ′ i )− Φu(T ) ) = k∑ i=1 ( |u− πTu|2 − |u− πT ′iu|2 ) By definition of Φ
= k∑ i=1 (√ |u|22 − |πTu|22 − √ |u|22 − |πT ′iu| 2 2 ) By Pythagorean Theorem
= k∑ i=1 ( g(|πTu|22)− g(|πT ′iu| 2 2) ) ≥ k∑ i=1 |g′(|πTu|22)| ( |πT ′iu| 2 2 − |πTu|22 ) By the Mean Value Theorem
= k∑ i=1
1 2 √ |u|22 − |πTu|22
( |πT ′iu| 2 2 − |πTu|22 ) = 1
2|u− πTu|2 k∑ i=1 ( |πT ′iu| 2 2 − |πTu|22 ) By Pythagorean Theorem
≥ 1 2|u− πTu|2 · σmin(BS)2 (|πSu|22 − |πTu|22)2 4|πSu|22 By Lemma 7.1 = σmin(BS) 2 (|u− πTu|22 − |u− πSu|22)2
8 · |πSu|22|u− πTu|2 By Pythagorean Theorem
= σmin(BS)
2
8 · |πSu|22|u− πTu|2 · (|u− πTu|2 − |u− πSu|2)2
· (|u− πTu|2 + |u− πSu|2)2
≥ σmin(BS)2 (Φu(S)− Φu(T ))2 · |u− πTu|2
8 · |πSu|22 Since |u− πSu|2 ≥ 0
We now lower bound |u−πT u|2|πSu|22 as follows.
|u− πTu|2 |πSu|22 = |u− πTu|2 |u|22 − |u− πSu|22 By Pythagorean Theorem
= |u− πTu|2 |u|2 − |u− πSu|2 · 1|u|2 + |u− πSu|2
= 1 Φu(S) · |u− πTu|2|u|2 + |u− πSu|2
By definition of Φ
≥ 1 2Φu(S) · |u− πTu|2|u|2 Since |u− πSu|2 ≤ |u|2
= 1 2Φu(S) · |u|2 − Φu(T )|u|2
By definition of Φ
= 1 2Φu(S) · ( 1− Φu(T )|u|2 ) ≥ 1
2Φu(S) · ( 1− Φu(T ) Φu(S) ) Since Φu(S) ≤ |u|2
= (Φu(S)− Φu(T ))
2Φu(S)2
Therefore,
k∑ i=1 ( Φu(T ′ i )− Φu(T ) ) ≥ σmin(BS)2 (Φu(S)− Φu(T ))2 · |u− πTu|2 8 · |πSu|22
≥ σmin(BS)2 (Φu(S)− Φu(T ))3
16Φu(S)2
Lemma 7.3. (Utility Improvement by Projecting a Matrix) Let A ∈ Rd×n, and T, S ⊂ [n] be two sets of column indices, with ΦA(S) ≥ ΦA(T ). Furthermore, let k = |S|. Then, there exists a column index i ∈ S such that
ΦA(T ∪ {i})− ΦA(T ) ≥ σmin(BS)2 (ΦA(S)− ΦA(T ))3
16kΦA(S)2
Proof. The proof mostly follows the proof of Lemma 1 in Altschuler et al. (2016). We combine Lemma 7.2 with Jensen’s inequality to conclude an improvement of the utility function with one additional column when projecting a matrix instead of a single column.
For j ∈ [n], we define δj = min(1, ΦA∗j (T )
ΦA∗j (S) ). Note that δj is 1 if the j-th column A∗j has a larger projection
onto BT than BS , and ΦA∗j (T ) ΦA∗j (S) otherwise. Let k = |S|. For i ∈ [k], let T ′i = T ∪ {i}.
1
σmin(BS)2 k∑ i=1 ( ΦA(T ′ i )− ΦA(T ) ) =
1
σmin(BS)2 n∑ j=1 k∑ i=1 ( ΦA∗j (T ′ i )− ΦA∗j (T ) ) By definition of Φ
≥ n∑ j=1 (1− δj)3 16 · ΦA∗j (S) By Lemma 7.2
= ΦA(S)
16 n∑ j=1 (1− δj)3 · ΦA∗j (S)∑n i=1 ΦA∗i(S)
Note ΦA(S) = n∑ i=1 ΦA∗i(S)
≥ ΦA(S) 16 ( n∑ j=1 (1− δj) · ΦA∗j (S)∑n i=1 ΦA∗i(S) )3 By Jensen’s Inequality
= 1
16ΦA(S)2 ( n∑ j=1 (1− δj) · ΦA∗j (S) )3
Since 1− δj ≥ 1− ΦA∗j (T )
ΦA∗j (S)
⇒ (1− δj) · ΦA∗j (S) ≥ ΦA∗j (S)− ΦA∗j (T )
≥ 1 16ΦA(S)2 ( n∑ j=1 (ΦA∗j (S)− ΦA∗j (T )) )3
= (ΦA(S)− ΦA(T ))3
16ΦA(S)2
Hence, k∑ i=1 ( ΦA(T ′ i )− ΦA(T ) ) ≥ σmin(BS)2 (ΦA(S)− ΦA(T ))3 16ΦA(S)2
This implies there is at least one column of BS , with index i ∈ S, such that when i is added to T , the utility function ΦA(T ) increases by at least 1k · σmin(BS) 2 (ΦA(S)−ΦA(T ))3 16ΦA(S) 2 .
Theorem 7. Let A ∈ Rd×n be the data matrix and k ∈ N be the desired rank. Let AS be the best possible subset of k columns, i.e., AS = arg minAS minV |ASV − A|1,2. Let σ be the minimum non-zero singular
value of the matrix B of normalized columns of AS , (the j-th column of B is B∗j = (AS)∗j/|(AS)∗j |2). Then, if T ⊂ [n] is the subset of columns selected by Greedy k-CSS1,2, the following holds with |T | = Ω( kσ2 2 ),
min V |ATV −A|1,2 ≤ (1− ) min S⊂[n],|S|=k,V ∈Rk×n |ASV −A|1,2 + |A|1,2
Proof. The proof follows the one of Theorem 1 in Altschuler et al. (2016).
Let Tt be the subset of columns of B selected by Greedy k-CSS1,2 after t iterations. Notice that T0 = ∅. In addition, define F = ΦA(S) = ΦA(S)− ΦA(T0), ∆0 = F , and ∆i = ∆02i for i ∈ N.
Let ∆i ≥ ΦA(S)− ΦA(Tt) ≥ ∆i+1 = ∆i2 . Our goal is to bound the number of iterations needed for the gap between ΦA(S) and ΦA(Tt) to become less than ∆i2 .
Consider some iteration s for which ΦA(S) − ΦA(Ts) ≥ ∆i2 . The improvement of the utility function after adding a column of B to Ts through greedy selection is at least the improvement of the utility function after adding the best column of BS to Ts. By Lemma 7.3, this is at least
σ2 (ΦA(S)− ΦA(T ))3 16kΦA(S)2 = σ2 · ∆ 3 i 16 · 8 · k · F 2 = σ2∆3i 128kF 2 .
If the gap ΦA(S)− ΦA(Tt) after t iterations is at most ∆i and at least ∆i2 , then after at most 64kF2
σ2∆2i iterations,
the gap becomes at most ∆i 2 .
We can use this to bound the number of iterations required for the gap to become at most F . Take N ∈ N such that ∆N+1 ≤ F ≤ ∆N . Then the number of iterations required for the gap to become at most ∆N+1 is at most
N∑ i=0 64kF 2 σ2∆2i = 64kF 2 σ2 N∑ i=0 1 ∆2i
= 64kF 2
σ2∆2N+1 N∑ i=0
1
4N+1−i Since ∆N+1 = ∆i 2N+1−i
≤ 256k 3σ2 2
Since ∆N+1 ≥ F
2 and N∑ i=0
1 4N+1−i ≤ 1 3
Therefore, after |T | = Ω( k δ2 2 ) iterations, we have
ΦA(S)− ΦA(T ) ≤ ΦA(S) ⇒ |A|1,2 − |A− πSA|1,2 − (|A|1,2 − |A− πTA|1,2) ≤ (|A|1,2 − |A− πSA|1,2) ⇒ |A− πTA|1,2 ≤ (1− )|A− πSA|1,2 + |A|1,2
Since S is the set of indices for the best possible subset of k columns, the above is equivalent to
min V |ATV −A|1,2 ≤ (1− ) min S⊂[n],|S|=k,V ∈Rk×n |ASV −A|1,2 + |A|1,2
We now analyze Algorithm 2, based on the analysis of the Lazier-than-lazy greedy heuristic in Altschuler et al. (2016). The first step is the following lemma based on Lemma 6 of Altschuler et al. (2016), 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 A ∈ Rd×n, and let T, S ⊂ [n] be two sets of column indices, with k := |S| and ΦA(S) ≥ ΦA(T ). Let T be a set of n log(1/δ)
k column indices of A, chosen uniformly at random from [n] \ T . Then,
E[max i∈T
ΦA(T ∪ {i})]− ΦA(T ) ≥ (1− δ) · σmin(BS)2 · (ΦA(S)− ΦA(T ))3
16kΦA(S)2
Proof. The proof is nearly identical to the proof of Lemma 6 of Altschuler et al. (2016) — we include the full proof for completeness. The first step in the proof is showing that T ∩ (S \T ) is nonempty with high probability. Then, by conditioning on T ∩ (S \ T ) being nonempty, we can show that the expected increase in utility is large. 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 T without replacement.
First, observe that Pr[T ∩ (S \ T ) = ∅] = O ( n log(1/δ) k )∏ t=1 ( 1− |S \ T | n− |T |
) = (
1− |S \ T | n− |T |
)O(n log(1/δ) k ) ≤ e− |S\T | n−|T | · n log(1/δ) k By 1− x ≤ e−x
≤ e− |S\T | log(1/δ)
k Because n− |T | < n meaning that
Pr[T ∩ (S \ T )] ≥ 1− e− |S\T | log(1/δ) k
= 1− δ |S\T | k
≥ (1− δ) |S \ T | k
Since |S \ T | ≤ k, and 1− δx ≥ (1− δ)x for x, δ ∈ [0, 1]
Therefore, E[max
i∈T ΦA(T ∪ {i})− ΦA(T )] ≥ Pr[T ∩ (S \ T ) 6= ∅] · E [
max i∈T
ΦA(T ∪ {i})− ΦA(T ) ∣∣∣T ∩ (S \ T ) 6= ∅]
≥ (1− δ) |S \ T | k
· E [
max i∈T
ΦA(T ∪ {i})− ΦA(T ) ∣∣∣T ∩ (S \ T ) 6= ∅]
≥ (1− δ) |S \ T | k
· E [
max i∈T
ΦA(T ∪ {i})− ΦA(T ) ∣∣∣|T ∩ (S \ T )| = 1]
(Since it is always better for T ∩ (S \ T ) to be larger)
= (1− δ) |S \ T | k
· ∑ i∈S\T (ΦA(T ∪ {i})− ΦA(T ))
|S \ T | = (1− δ) · ∑ i∈S(ΦA(T ∪ {i})− ΦA(T ))
|S| (Since ΦA(T ∪ {i}) = ΦA(T ) for i ∈ T )
≥ (1− δ) · 1|S| · σmin(BS) 2 (ΦA(S)− ΦA(T ))2 16ΦA(S)2
(See the proof of Lemma 7.3.) This proves the lemma in the case where the 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 of the probability mass from this instance of T to other sets T
′ ⊂ [n] \ T , which contain T but do not have repeated elements. This leads to the uniform distribution on subsets of [n] \ T with no repeated elements, i.e., the distribution that results from sampling without replacement.
Using this lemma, we analyze the convergence of Algorithm 2:
Theorem 8. Let A ∈ Rd×n be the data matrix and k ∈ N the desired rank. Let AS be the best subset of k columns, i.e., AS = arg minAS minV |ASV − A|1,2. Let σ be the minimum non-zero singular value of the matrix B of normalized columns of AS (meaning the j-th column of B is B∗j = (AS)∗j/|(AS)∗j |2). Then, if T ⊂ [n] is the subset of columns selected by Algorithm 2, the following holds if |T | = Ω( k
σ2 2 ):
E[min V |ATV −A|1,2] ≤ (1− ) min S⊂[n],|S|=k,V ∈Rk×n |ASV −A|1,2 + |A|1,2
Proof. The proof is uses the same strategy as that of Theorem 5 of Altschuler et al. (2016) (and Theorem 7 above), with minor modifications. Let Tt be the subset of columns of B selected by Algorithm 2 after t iterations (in particular, T0 = ∅). In addition, let F = ΦA(S) = ΦA(S)− ΦA(T0), ∆0 = F , and ∆i+1 = ∆i2 . Now, fix a time t such that for some i, ∆i ≥ ΦA(S) − ΦA(Tt) ≥ ∆i+1 = ∆i2 . Then, we bound the number of additional iterations t′ needed so that
E[ΦA(S)− ΦA(Tt+t′) | Tt] < ∆i+1
For convenience, for each k ≥ 0, define Ek := ΦA(Tt+k). Then, our goal is to find t′ such that ΦA(S)− Et′ < ∆i+1
However, observe that from Lemma 7.4 above, we obtain Ek+1 − Ek = E [ ΦA(Tt+k+1)− ΦA(Tt+k) ∣∣∣Tt] = E [ E [ ΦA(Tt+k+1)− ΦA(Tt+k)
∣∣Tt+k]∣∣∣Tt] By E[E[X|Y ]] = E[X] ≥ E [ (1− δ) · σmin(BS)2 · (ΦA(S)− ΦA(Tt+k))3
16kΦA(S)2 ∣∣∣Tt] By Lemma 7.4 = (1− δ) · σmin(BS)2
16kΦA(S)2 · E [ (ΦA(S)− ΦA(Tt+k))3|Tt ] ≥ (1− δ) · σmin(BS) 2
16kΦA(S)2 · ( E[ΦA(S)− ΦA(Tt+k)|Tt] )3 By Jensen’s Inequality
= (1− δ) · σmin(BS)2 · (E[ΦA(S)]− Ek)3
16kΦA(S)2
Now, suppose that ∆i ≥ ΦA(S) − Es ≥ ∆i+1, for s = 0, . . . , t′ − 1. Then, for all such s, Es+1 − Es ≥ (1−δ)σmin(BS)2∆3i+1
16kF2 . Summing these inequalities for s = 0, . . . , t′ − 1, we find that
Et′ − E0 ≥ (1− δ)σmin(BS)2
16kF 2 ·∆3i+1 · t′
and for the increase from E0 to Et′ to be greater than ∆i+1, it suffices to have
t′ ≥ 32kF 2
∆2i+1 · (1− δ)σmin(BS)2
In summary, if ΦA(S) − E[ΦA(Tt)] ≤ ∆i, then in at most 32kF 2
∆i+1·(1−δ)σmin(BS)2 iterations, ΦA(S) −
E[ΦA(Tt)] ≤ ∆i+1. Thus, if we let N ∈ N such that ∆N+1 ≤ √1−δF ≤ ∆N , then the number of iterations t needed to have ΦA(S)− E[ΦA(Tt)] < ∆N+1 is at most N∑ i=0 32kF 2 ∆2i+1 · (1− δ)σmin(BS)2 = 32kF 2 (1− δ)σmin(BS)2 N∑ i=0 1 ∆2i+1
= 32kF 2 (1− δ)σmin(BS)2 N∑ i=0 1 4N−i · 1 ∆2N+1
≤ 32kF 2 (1− δ)σmin(BS)2 · 4(1− δ) 2F 2 N∑ i=0 1 4N−i Since ∆N+1 ≥
F
2 √ 1− δ
= 128k
σmin(BS)2 2 N∑ i=0 1 4i
≤ 512k 3σmin(BS)2 2
Thus, after t = O( k σmin(BS) 2 2 ) iterations,
ΦA(S)− E[ΦA(Tt)] ≤ √
1− δ ΦA(S)
meaning
|A|1,2 − |A− πSA|1,2 − E[|A|1,2 − |A− πTA|1,2] ≤ √ 1− δ |A|1,2 − √ 1− δ |A− πSA|1,2
and rearranging gives E[|A− πTA|1,2] ≤ (
1− √ 1− δ
) |A− πSA|1,2 +
√ 1− δ |A|1,2
This completes the proof (note that we can select δ = , meaning 1√ 1−δ = O(1) for < 1 2 ).

heading: E EXTENSION TO THE STREAMING MODEL
text: In this section, we describe how our protocol in Algorithm 1 can be made into a 1-pass streaming algorithm for column subset selection in the `p norm. The algorithm is shown in Algorithm 7, and is analyzed in Theorem 9 below. The algorithm and its analysis follow the standard merge-and-reduce framework (see McGregor (2014)).
Theorem 9 (Analysis of Algorithm 7). Let A ∈ Rd×n and k ∈ N, and assume Algorithm 7 sees the columns of

heading: A one at a time in the stream S. Then, Algorithm 7 returns U ∈ Rd×Õ(k) such that
text: min V ∈RÕ(k)×n ‖UV −A‖p ≤ Õ(k1/p−1/2) min Ak rank k ‖Ak −A‖p
with probability 0.9. Moreover, the space complexity of Algorithm 7 is Õ(dk).
Proof. Let r = Õ(k) be the bi-criteria rank of Algorithm 4. Then, we will be repeatedly applying Lemma 4 with k being equal to r, i.e., 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, t) of L can be thought of as holding a coreset of the columns of A in some interval I in [n]. We prove the following intermediate lemma by induction on t:
Lemma 10. Let B be the concatenation of all the sketched columns in L (which have been reweighted by multiple applications of Lemma 4). Then, for all subspaces V ⊂ RÕ(k) of dimension at most r,
‖B − PVB‖p,2 = (
1± 1 logn
)t ‖SAI − PV SAI‖p,2
Proof. We proceed by induction on t. The lemma is clear when t = 0, since in that case, the columns of B are simply sketched columns of A which have not been re-weighted.
Now, suppose t > 0, and suppose the lemma holds for smaller values of t. Then, the sketched columns in Lmust have been obtained as follows: there previously existed two elements (L1, t− 1) and (L2, t− 1) of L, such that if B1 is the concatenation of the sketched columns in L1 and B2 is the concatenation of the sketched columns in L2, and T is a coreset for the concatenation B3 of B1 and B2, then B = B3T . The sketched columns in L1 and L2 form coresets for two intervals in [n], which we denote by I1 and I2 respectively. Applying Lemma 4, we find that
‖B − PVB‖pp,2 = (
1± 1 logn
)p ‖B3 − PVB3‖pp,2
= (
1± 1 logn
)p( ‖B1 − PVB1‖pp,2 + ‖B2 − PVB2‖ p p,2 ) = (
1± 1 logn
)p · (
1± 1 logn
)p(t−1)( ‖SAI1 − PV SAI1‖ p p,2 + ‖SAI2 − PV SAI2‖ p p,2 ) = (
1± 1 logn
)pt · ‖SAI − PV SAI‖pp,2
(1)
where the second and last equalities are because the pth power of the `p,2 norm of a matrix decomposes across the columns of the matrix, and the third equality is by the induction hypothesis. By taking pth roots, we find that
‖B − PVB‖pp,2 = (
1± 1 logn
)t ‖SAI − PV SAI‖pp,2
This proves the lemma.
Now, if L is the unique element of L remaining at the end of the for loop in Algorithm 7, let B be the concatenation of the sketched columns in L. Then, L is a coreset for all the columns of A, and by the above lemma, the distortion of L is at most (1± 1
logn )logn ∈ [ 1 e , e] (since t ≤ logn — note that, as in all applications
of the merge-and-reduce framework, the coresets contained in L over the course of the algorithm form a binary tree, with the leaf nodes being contiguous intervals of length Õ(k)). Hence, for all subspaces V of dimension at most r, ‖B − PVB‖p,2 = Θ(1)‖SA− PV SA‖p,2 and in particular, if M ∈ RÕ(k)×Õ(k) is the matrix formed by concatenating the columns of B selected by running Algorithm 4, then
‖SA− PMSA‖p,2 ≤ Θ(1) min T⊂[n],|T |≤k ‖SA− PTSA‖p,2
where PM on the left-hand side denotes the projection onto the column span of M , and PT on the right-hand side denotes the projection to the column span of SAT . Hence, by Lemma 1,
‖SA− PMSA‖p ≤ Õ(k1/p−1/2) min T⊂[n],|T |≤k ‖SA− PTSA‖p
meaning that if U ∈ Rd×Õ(k) is the matrix whose columns are the unsketched columns corresponding to M , then
min V ∈RÕ(k)×d ‖SA− SUV ‖p ≤ min T⊂[n],|T |≤k,V ∈Rk×d ‖SA− SATV ‖p
and by Lemmas 2 and 3, this means
min V ∈RÕ(k)×d ‖A− UV ‖p ≤ min T⊂[n],|T |≤k,V ∈Rk×d ‖A−ATV ‖p
Now, we analyze the space complexity of Algorithm 7. Note that at any iteration of the algorithm, L can hold at most logn lists of columns (since each element of L is a coreset corresponding to an interval of column indices in [n] of size 2k for some k ∈ N, and if two adjacent coresets are of the same size then they will have been merged, so all coresets in L are of different sizes). Each list of columns is of size at most Õ(k), and each column has d entries, meaning the amount of space used in any iteration is at most Õ(dk).

heading: F COMPARISON OF OUR PROTOCOL WITH THE DISTRIBUTED GREEDY PROTOCOL OF ALTSCHULER ET AL. (2016) FOR THE FROBENIUS NORM
text: In this section, we perform an empirical comparison of our protocol with the distributed greedy protocol for column subset selection in the Frobenius norm due to Altschuler et al. (2016).

heading: F.1 DISTRIBUTED GREEDY PROTOCOL OF ALTSCHULER ET AL. (2016)
text: We first recall the distributed protocol of Altschuler et al. (2016), in Algorithm 8. Note that it is a bi-criteria algorithm, i.e., more than k columns are selected, and in Altschuler et al. (2016) it is shown that this gives a good approximation relative to the optimal column subset for the Frobenius norm.
Naı̈vely, this would require a large communication cost, since the entire data matrix A would have to be communicated between the servers and the coordinator, in order for them to perform the calls GREEDY(A, Ti, 32kσ ) and GREEDY(A, T, 12k
σ ). Instead, Altschuler et al. (2016) uses the technique of Projection-Cost Preserving
sketches, developed in Cohen et al. (2015), in which a key result is as follows:
Theorem 11 (Projection-Cost Preserving Sketches - Theorem 4 of Cohen et al. (2015), as stated in Altschuler et al. (2016)). Let R be a random matrix with n rows and n′ = O( k+log( 1 δ )
2 ) columns, where each entry is independently and uniformly set to ± √
1 n′ . Then for any matrix A ∈ R d×n, with probability 1 − O(δ), the following holds: for some constant c ≥ 0, for any k-dimensional subspace U of Rd, if PU ∈ Rd×d is the corresponding projection matrix, then
(1− )|PUA|22 ≤ |PU (AR)|22 + c ≤ (1 + )|PUA|22
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 of Rd.
Hence, when each server performs the call GREEDY(A, Ti, 32kσ ), in place of A, it can use a projection-cost preserving sketchAR forA (and the coordinator can similarly useAR when it makes the call GREEDY(A, T, 12k
σ )).
This is because the calls to GREEDY repeatedly compute the cost of the projection of the matrix A onto various column subsets of Ti or T , and hence all that is needed is a way to efficiently compute the cost of the projection of A onto various subsets of size O( k
σ ).
AR can be computed without communicating the entire matrix A between the servers, as follows:
• First, the coordinator generates R ∈ Rn×n ′ (where n′ = O(k)) as specified in Theorem 11. This is sent to all the servers (which can be done with negligible cost if a random seed is sent, for example).
• Suppose the ith server holds Ai, which consists of columns with indices si through ti of A. Then, if Ri is the (ti−si+1)×n′ submatrix ofR consisting of rows si through ti, then the ith server can send AiRi to the coordinator. Since AiRi is a d×O(k) matrix, this step takes O(sdk) communication.
• Finally, by definition of matrix multiplication, AR = ∑s i=1 AiRi. The coordinator performs this
computation and sends AR to each server. This step also takes O(sdk) communication.
Algorithm 7 1-pass streaming algorithm for k-CSSp. L is a collection of coresets of columns. Over the course of the algorithm, each element L ∈ L will represent a contiguous subset of 2t columns, for some t ∈ [log n] — to determine when two coresets should be merged, we will also keep track of the size of each coreset in L. Hence, each element of L is of the form (L, t) where L is a list of sketched columns (and their corresponding unsketched columns) and t is the number of times this list has been involved in a merging operation.
Input: A stream S in which the columns of the data matrix A ∈ Rd×n arrive one at a time, and the target rank k ∈ N Output: The left factor U ∈ Rd×Õ(k) S ← An Õ(k)× d random matrix with i.i.d. standard p-stable entries ← 1logn δ ← O( 1n ) f ← Õ(k) L ← ∅ for Each column A∗j that arrives from S do
if L is empty. then L← {}, where {} is the empty list. else if The last element (J, t) of L is such that t = 0 (i.e., it has been merged 0 times). then L← J Remove (J, t) from the end of L. else L← {} end if L← L ∪ {(SA∗j , A∗j)} L ← L ∪ (L, 0)
/* Now, we merge coresets in L as much as possible. */ while True do
Exit this while loop if L has only 1 element. (L, t)← last element of L— remove this from L. (L′, t′)← last element of L— remove this from L. If t 6= t′, then add (L′, t′) back to L, and add (L, t) back to L. Then, exit this while loop. C ← A coreset of L ∪ L′, computed as specified in Lemma 4 — the k in the statement of Lemma 4 will be the bicriteria rank of Algorithm 4, which is Õ(k). The parameters δ 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-weighting them)). L ← (C, t + 1) — note that for each of the re-scaled columns that are included in C, we include their original indices in A as well.
end while end for
L← The unique element of L L′ ← The result of running Algorithm 4 on the sketched columns in L— for each of these sketched columns, store the corresponding unsketched column as well. U ← The d× Õ(k) matrix whose columns are the unsketched columns in L′ return U
Algorithm 8 Distributed Greedy Column Subset Selection for the Frobenius Norm (Algorithm 2 of Altschuler et al. (2016)). Throughout this protocol, GREEDY(A, T, r) denotes a single-machine procedure, which does the following: if A ∈ Rd×n is a data matrix, and T ∈ Rd×t is a set of columns (not necessarily of A) then a subset S of r columns of T is constructed iteratively over r steps, such that at each step, the new column of T to add to S is greedily chosen — that is, the chosen column increases |πSA|22 the most, or equivalently, decreases |A− πSA|22 the most. In other words, it is the same as our Algorithm ??, but for the Frobenius norm rather than the `1,2-norm.
Input: The data matrix A ∈ Rd×n, target rank k ∈ N, the number of servers s ∈ N. The columns of A are assumed to be randomly partitioned among servers T1, T2, . . . , Ts. Output: A subset of columns AT from A. |T | = O( kσ ), where if AOPT is the optimal subset of columns of A of size k, and the columns of AOPT are normalized to have unit `2 norm, then σ is the smallest singular value of AOPT . Si ← GREEDY(A, Ti, 32kσ ) for all i ∈ [s] (The i
th server performs this computation.) Each server sends its Si to the coordinator. T ← ∪si=1Si (This computation is done by the coordinator.) S ← GREEDY(A, T, 12kσ ) The coordinator returns arg maxS′∈{S,S1,S2,...,Ss} |πS′A| 2 2
In our implementation of Algorithm 8, 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 LAZIERTHAN-LAZY-GREEDY ALGORITHM. The difference between the GREEDY(A, T, r) algorithm and the LAZIER-THAN-LAZY-GREEDY(A, T, r) algorithm is as follows: while at each of the r iterations, GREEDY considers all columns of T and chooses the one that leads to the most improvement in the objective, LAZIERTHAN-LAZY-GREEDY samples |T | log( 1 δ )
r columns uniformly at random for some small δ (which we take in our
implementation to be 0.005), and out of those columns, chooses the one which leads to the greatest improvement (where |T | is the number of columns in T ). This leads to a significant speedup to GREEDY, and it was shown in Altschuler et al. (2016) that this does not significantly worsen the approximation guarantees that Altschuler et al. (2016) shows for GREEDY.

heading: F.2 SETUP
text: We compare our distributed protocol, in the case p = 1, to Algorithm 8. For both protocols, at the outset we fix the number of columns which are selected. Algorithm 8 is rewritten as Algorithm 9 to reflect this. In this section, we use k to denote the number of columns ultimately selected.
Algorithm 9 Distributed Greedy Column Subset Selection for the Frobenius Norm (Algorithm 2 of Altschuler et al. (2016)). For our empirical comparison, we fix the number of columns selected at the outset.
Input: The data matrix A ∈ Rd×n, desired number of columns k ∈ N, the number of servers s ∈ N. The columns of A are assumed to be randomly partitioned among servers T1, T2, . . . , Ts. Output: A subset of columns AT from A, with |T | = k. AR ← a PCP for A as discussed in the previous section. All servers and the coordinator have access to AR. Si ← GREEDY(AR, Ti, k) for all i ∈ [s] (The ith server performs this computation.) Each server sends its Si to the coordinator. T ← ∪si=1Si (This computation is done by the coordinator.) S ← GREEDY(AR, T, k) The coordinator returns arg maxS′∈{S,S1,S2,...,Ss} |πS′(AR)| 2 2
Note that in Algorithm 9, the coordinator now chooses as many columns as chosen by the servers, as opposed to Algorithm 8, where it chooses a 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 lesions, a 76×698 matrix dataset available at https://archive.ics.uci.edu/ ml/datasets/Gastrointestinal+Lesions+in+Regular+Colonoscopy.
• secom, a 591× 1567 matrix dataset. secom has missing entries, which we replace with 0s for the purposes of our experiments. secom is available at https://archive.ics.uci.edu/ml/ datasets/SECOM.

heading: F.2.1 PARAMETERS
text: We compare our distributed protocol for k-CSS1, using Greedy k-CSS1,2 as a subroutine, with Algorithm 9, for several choices of k — on gastro lesions, we let k ∈ {10, 20, 30}, while on secom, we let k ∈ {30, 60, 90, 120}. For our protocol, cauchy size is set to 2k and coreset size is set to 5k for both datasets (where cauchy size and coreset size have the same meanings as in the main body of this paper.) For Algorithm 9, the number of columns in the PCPs is set to 8k on gastro lesions, and 7k on secom.
The hyperparameters coreset size and cauchy size, 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 is allowed slightly more communication). To see this, let d be the number of rows in the data matrix A. If c is the number of columns in the PCPs, then the total communication used to transmit the PCPs between servers is 2scd, since the servers must first send their respective AiRi to the coordinator, and the coordinator then sends AR to all of the servers. By comparison, if r1 is the number of rows in the initial Cauchy matrix in our protocol, and r2 is 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 is (r1 + r2)sd. We choose r1, r2 and c so that 2scd is slightly higher than (r1 + r2)sd.
Note that if k is the number of columns ultimately selected, then our protocol uses 2dk additional bits of communication between the servers and the coordinator to recover the final k-subset of columns, while Algorithm 9 will use sdk communication to send the subsets of columns Si (of size k) from the servers to the coordinators. This is not included in our hyperparameter calculations.

heading: F.2.2 HOW TRIALS ARE CONDUCTED
text: With these choices of hyperparameters, we conduct 15 trials as follows. In each trial, if A ∈ Rd×n 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 9 to apply, the columns should be partitioned randomly). Using this partition across 2 servers, we run our protocol and Algorithm 9. For each protocol, once the subset of columns is computed, we perform multiple-response `1 regression to evaluate the error.
In the next section, we report the minimum error across the 15 trials for each dataset and for each value of k. We also report the mean error, along with the standard deviation. Finally, we compute the work and the span of our protocol and Algorithm 9 using Python’s time.process time() utility. This does not include the time taken to perform multiple-response `1 regression.
For secom, with k = 30 and k = 60, the trials were performed on a Late-2016 Macbook Pro with a 2.7 GHz Quad-Core Intel Core i7 processor and 16GM of memory. The rest were performed on a 2019 MacBook Pro with a 2.4 GHz Intel Core i5 processor and 8 GB 2133 MHz LPDDR3 memory.

heading: F.3 RESULTS
text: Our results are shown and discussed in Figures 3 (minimum and mean/standard deviation for `1 error) and 4 (average work and span across 15 trials).

heading: G FULL EXPERIMENTAL RESULTS
text: For our protocol, we considered various additional hyperparameter settings on real-world datasets (bcsstk13, isolet and 5 images in the Caltech-101 dataset) — these settings are shown in Tables 3, 4, and 5. As before, cauchy size is the number of rows in our initial Cauchy matrix, sent to the servers at the beginning of the protocol, and coreset size 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 of 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 perform better than all settings of regular k-CSS1,2 in minimum and mean errors across multiple trials, it also has a smaller variance in performance.

heading: G.1 DETAILS OF WORK/SPAN COMPUTATION
text: In Figure 6, work and span were recorded as the time taken (using Python’s time.process time() utility) to compute the column subset in our distributed protocol. In particular, it does not include the time spent performing `1 regression to obtain the right factor and consequently the entry-wise `1 errors.
Figure 5 shows that we encounter a tradeoff between accuracy and running time when choosing between these two subroutines. Both lead to the same overall communication cost, and since accuracy is of more interest than running time in the distributed setting, it is (empirically) preferable to use Greedy k-CSS1,2 within the protocol.

heading: G.2 ADDITIONAL DETAILS
text: All experiments on the Caltech-101 dataset were run on a Late-2016 Macbook Pro with a 2.7 GHz QuadCore Intel Core i7 processor and 16GM of memory. All experiments on bcsstk13 and isolet were run on an AWS z1d.xlarge instance with Deep Learning AMI (Amazon Linux 2) Version 29.0.

Review summaries:
The 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. 

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. 

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.