Efficient volume sampling for row/column subset selection

Volume sampling, i.e., picking k-subsets of the rows of any given matrix with probabilities proportional to the squared volumes of the simplicies defined by them, was introduced in [5] in the context of low-rank approximation of matrices. It is equivalent to sampling k-subsets of {1, . . . , m} with probabilities proportional to the corresponding k by k principal minors of any given m by m positive semidefinite matrix. In the context of low-rank approximation, volume sampling is related to a problem called row/column-subset selection [1]. Most large data sets that arise in search, microarray experiments, computer vision, data mining etc. can be thought of as matrices where rows and columns are indexed by objects and features, respectively (or vice versa), and we need to pick a small subset of features that are dominant. For example, while studying gene expression data biologists want a small subset of genes that are responsible for a particular disease. Usual dimension reduction techniques such as principal component analysis (PCA) or random projection fail to do this as they typically output singular vectors or random vectors which are linear combinations of a large number of feature vectors. A recent article by Mahoney and Drineas [18] highlights the limitations of PCA and gives experimental data on practical applications of low-rank approximation based on row/column-subset selection. While dealing with large matrices in practice, we seek smaller or low-dimensional representations of them which are close to them but can be computed and stored efficiently. A popular notion for low-dimensional representation of matrices is low-rank matrices, and the most popular metrics used to measure the closeness of two matrices are the Frobenius or Hilbert-Schmidt norm (i.e., the square root of the sum of squares of entries of their difference) and the spectral norm (i.e., the largest singular value of their difference). The singular value decomposition (SVD) tells us that any matrix A ∈ R m×n can be written as where σ 1 ≥ . . . ≥ σ m ≥ 0, u i ∈ R m are orthonormal and v i ∈ R n are orthonormal. Moreover, the nearest rank-k matrix to A, let us call it A k , under both the Frobenius and the spectral norm, is given by In other words, the rows of A k are projections of the rows of A onto span (v i : 1 ≤ i ≤ k). Because of this, most dimension reduction techniques based on the singular value decomposition, e.g., principal component analysis (PCA), are interpreted as giving v i 's as the dominant vectors, which happen to be linear combinations of a large number of the rows or feature vectors of A. The row-subset selection problem we consider in this paper is: Can we pick a k-subset of the rows (or feature vectors) of A ∈ R m×n so that projecting onto their span is almost as good as projecting onto span (v i : 1 ≤ i ≤ k)? Several row-sampling techniques have been considered in the past as an approximate but faster alternative to the singular value decomposition, in the context of streaming algorithms and large data sets that cannot be stored in random access memory [8,7,5]. The first among these is the squared-length sampling of rows introduced by Frieze, Kannan and Vempala [8]. Another sampling scheme due to Drineas, Mahoney and Muthukrishnan [7] uses the singular values and singular vectors to decide the sampling probabilities. Later Deshpande, Rademacher, Vempala and Wang [5] introduced volume sampling as a generalization of squared-length sampling. where a i denotes the i-th row of A, A S ∈ R k×n denotes the row-submatrix of A given by rows with indices i ∈ S, and conv (•) denotes the convex hull. The application of volume sampling to low-rank approximation and, more importantly, to the row-subset selection problem, is given by the following theorem shown in [5]. It says that picking a subset of k rows according to volume sampling and projecting all the rows of A onto their span gives a (k + 1)-approximation to the nearest rank-k matrix to A. when S is picked according to volume sampling, π S (A) ∈ R m×n denotes the matrix obtained by projecting all the rows of A onto span (a i : i ∈ S), and A k is the matrix of rank k closest to A under the Frobenius norm. As we will see later, this easily implies Theorem 2 gives only an existence result for row-subset selection and we also know a matching lower bound that says this is the best we can possibly do. Theorem 3. [5] For any ǫ > 0, there exists a matrix A ∈ R (k+1)×k such that picking any k-subset S of its rows gives However, no efficient algorithm was known for volume sampling prior to this work. An algorithm mentioned in Deshpande and Vempala [6] does k!-approximate volume sampling in time O(kmn), which means that plugging it in Theorem 2 can only guarantee (k + 1)!-approximation instead of (k + 1). Finding an efficient algorithm for volume sampling is mentioned as an open problem in the recent monograph on spectral algorithms by Kannan and Vempala (see Section 7.4 of [15]). Boutsidis, Drineas and Mahoney [1] gave an alternative approach to row-subset selection (without going through volume sampling) and here is a re-statement of the main theorem from their paper which uses columns instead of rows. Row/column-subset selection problem is related to rank-revealing decompositions considered in linear algebra [12,19], and the previous best algorithmic result for row-subset selection in the spectral norm case was given by a result of Gu and Eisenstat [12] on strong rank-revealing QR decompositions. The following theorem is a direct consequence of [12] as pointed out in [1]. Theorem 5. Given A ∈ R m×n , an integer k ≤ n and f ≥ 1, there exists a k-subset S of the columns of A such that Moreover, this subset S can be found in time In the context of volume sampling, it is interesting to note that Pan [19] has used an idea of picking submatrices of locally maximum volume (or determinants) for rank-revealing matrix decompositions. We refer the reader to [19] for details. The results of Goreinov, Tyrtyshnikov and Zamarashkin [10,11] on pseudo-skeleton approximations of matrices look at submatrices of maximum determinants as good candidates for row/columnsubset selection. Theorem 6. [10] If A ∈ R m×n can be written as where A 11 ∈ R k×k is the k by k submatrix of A of maximum determinant. Then, Because of this relation between row/column-subset selection and the related ideas about picking submatrices of maximum volume, C ¸ivril and Magdon-Ismail [3,4] looked at the problem of picking a k-subset S of rows of a given matrix A ∈ R m×n such that det A S A T S is maximized. They show that this problem is NP-hard [3] and moreover, it is NP-hard to even approximate it within a factor of 2 ck , for some constant c > 0 [4]. This is interesting in the light of our results because we show that even though finding the row-submatrix of maximum volume is NP-hard, we can still sample them with probabilities proportional to their volumes in polynomial time. Our main result is a polynomial time algorithm for exact volume sampling. In Section 4, we give an outline of our Algorithm 1, followed by two possible subroutines given by Algorithms 2 and 3 that could be plugged into it. Theorem 7 (polynomial-time volume sampling). The randomized algorithm given by the combination of the algorithm outlined in Algorithm 1 with Algorithm 2 as its subroutine, when given a matrix A ∈ R m×n and an integer 1 ≤ k ≤ rank(A), outputs a random k-subset of the rows of A according to volume sampling, using O(kmn ω log n) arithmetic operations. The basic idea of the algorithm is as follows: instead of picking a k-subset, pick an ordered ktuple of rows according to volume sampling (i.e., volume sampling suitably extended to all k-tuples such that for any fixed k-subset, all its k! permutations are all equally likely). We observe that the marginal distribution of the first coordinate of such a random tuple can be expressed in terms of coefficients of the characteristic polynomials of AA T and B i B T i , where B i ∈ R m×n is the matrix obtained by projecting each row of A orthogonal to the i-th row a i . Using this interpretation, it is easy to sample the first index of the k-tuple with the right marginal probability. Now we project the rows of A orthogonal to the chosen row and repeat to pick the next row, until we have picked k of them. The algorithm just described informally, if implemented as stated, would have a polynomial dependence in m, n and k, for some low-degree polynomial. We can do better and get a linear dependence in m by working with A T A in place of AA T and computing the projected matrices using rank-1 updates (Theorem 7), while still having a polynomial time guarantee and sampling exactly. It would be even faster to perform rank-1 updates to the characteristic polynomial itself, but that requires the computation of the inverse of a polynomial matrix (Proposition 18), and it is not clear to us at this time that there is a fast enough exact algorithm that works for arbitrary matrices. Jeannerod and Villard [14] give an algorithm to invert a generic n-by-n matrix with entries of degree d, with n a power of two, in time O(n 3 d). This would lead to the computation of all marginal probabilities for one row in time O(n 3 + mn 2 ) (a variation of Algorithm 3 and its analysis). Instead, if we are willing to be more practical, while sacrificing our guarantees, then we can perform rank-1 updates to the characteristic polynomial by using the singular value decomposition (SVD). In [9], an algorithm with cost O(min{mn 2 , m 2 n}) arithmetic operations is given for the singular value decomposition but the SVD cannot be computed exactly and we do not know how its error propagates in our algorithm which uses many such computations. If the SVD of an mby-n matrix can be computed in time O(T svd ), this leads to a nearly-exact algorithm for volume sampling in time O(kT svd + kmn 2 ). See Proposition 18 for details. Volume sampling was originally defined in [5] to prove Theorem 2, in particular, to show that any matrix A contains k rows in whose span lie the rows of a rank-k approximation to A that is no worse than the best in the Frobenius norm. Efficient volume sampling leads to an efficient selection of k rows that satisfy this guarantee, in expectation. In Section 5, we use the method of conditional expectations to derandomize this selection. This gives an efficient deterministic algorithm (Algorithm 4) for row-subset selection with the following guarantee in the Frobenius norm. This guarantee immediately implies a guarantee in the spectral norm, as follows: Theorem 8 (deterministic row subset selection). Deterministic Algorithm 4, when given a matrix A ∈ R m×n and an integer 1 ≤ k ≤ rank(A), outputs a k-subset S of the rows of A, using O(kmn ω log n) arithmetic operations, such that This improves the O(k √ log k)-approximation of Boutsidis, Drineas and Mahoney [1] for the Frobenius norm case and matches the lower bound shown in Theorem 3 due to [5]. The superlinear dependence on n might be too slow for some applications, while it might be acceptable to perform volume sampling or row/column-subset selection approximately. Our volume sampling algorithm (Algorithm 1) can be made faster, while losing on the exactness, by using the idea of random projection that preserves volumes of subsets. Magen and Zouzias [17] have the following generalization of the Johnson-Lindenstrauss lemma: for m points in R n there exists a random projection of them into R d , where d = O k 2 log m/ǫ 2 , that preserves the volumes of simplices formed by subsets of k or fewer points within 1 ± ǫ. Therefore, we get a (1 ± ǫ)approximate volume sampling algorithm that requires O(mnd)-time to do the random projection (by matrix multiplication) and then O(mnd ω log d) time for volume sampling on the new m-by-d matrix (according to Theorem 7). Theorem 9 (fast volume sampling). Using random projection for dimensionality reduction, the polynomial time algorithm for volume sampling mentioned in Theorem 7 (i.e., Algorithm 1 with Algorithm 2 as its subroutine), gives (1 + ǫ)-approximate volume sampling, using arithmetic operations. Finally, we show a lower bound for row/column-subset selection in the spectral norm that almost matches our upper bound in terms of the dependence on n. where π {i} (A) ∈ R n×(n+1) is the matrix obtained by projecting each row of A onto the span of its i-th row a i . For m ∈ N, let [m] denote the set {1, . . . , m}. For any matrix A ∈ R m×n , we denote its rows by a 1 , a 2 , . . . , a m ∈ R n . For S ⊆ [m], let A S be the row-submatrix of A given by the rows with indices in S. By span (S) we denote the linear span of {a i : i ∈ S} and let π S (A) ∈ R m×n be the matrix obtained by projecting each row of A onto span (S). Hence, Aπ S (A) ∈ R m×n is the matrix obtained by projecting each row of A orthogonal to span (S). Throughout the paper we assume m ≥ n. This assumption is not needed most of the time, but justifies sometimes working with A T A instead of AA T and, more generally, some choices in the design of our algorithms. It is also partially justified by our use of a random projection as a preprocessing step that makes n small. The singular values of A ∈ R m×n are defined as the positive square-roots of the eigenvalues of AA T ∈ R m×m (or A T A ∈ R n×n , up to some extra singular values equal to zero), and we denote them by where σ 1 , . . . , σ m are the singular values of A, i.e., eigenvalues of AA T , and Let ω be the exponent of the arithmetic complexity of matrix multiplication. We use that there is an algorithm for computing the characteristic polynomial of an n-by-n matrix using O(n ω log n) arithmetic operations [2,Section 16.6]. Here is another lemma that we will need about dividing determinants into products of two determinants. Proof. Without loss of generality, we can reduce ourselves to the case where S ∪ T is all rows of the given matrix: Then what we want to prove can be rewritten as: To show this, we consider two cases. If C S C T S is singular, then both sides of the equality are zero. If C S C T S is invertible, then we can perform block Gaussian elimination and write Writing the determinants of the block-triangular matrices gives Now, the projection of the rows of a matrix K onto the row-space of a matrix L can be written as This completes the proof. Finally, a well-known lemma about how the determinant of a matrix changes under a rank-1 update. We first outline our volume sampling algorithm to convince the reader that volume sampling can be done in polynomial time. In the subsequent subsections, we give improved subroutines to get faster implementations of the same idea. The main idea behind our algorithm is based on Lemma 14 about the marginal probabilities encountered in volume sampling. To explain this, it is more convenient to look at volume sampling defined as a distribution on k-tuples (X 1 , X 2 , . . . , X k ) instead of k-subsets, where each of the k! permutations of a k-subset is equally likely, i.e., for any (i Then the marginal probabilities Pr ) have the following interpretation in terms of the coefficients of certain characteristic polynomials. Proof. since the extra terms in the sum are all zero With this lemma in hand, let us consider the following outline of our algorithm. We will later give two more efficient implementations of this outline, depending on how the p i 's are computed. Input: a matrix A ∈ R m×n and 1 ≤ k ≤ rank(A). Output: a subset S of k rows of A picked with probability proportional to det A S A T S . 1. Initialize S ← ∅ and B ← A. For t = 1 to k do: (a) For i = 1 to m compute: where C i = Bπ {i} (B) is a matrix obtained by projecting each row of B orthogonal to b i . (b) Pick i with probability proportional to p i . Let S ← S ∪ {i} and B ← C i . 2. Output S. Proposition 15. The probability that our volume sampling algorithm outlined above picks a ksubset S is proportional to det A S A T S . This algorithm can be implemented with a cost of O(km 3 n+ km ω+1 log m) arithmetic operations. Proof. By Lemma 14, for any i 1 , i 2 , . . . , i k such that Pr (X 1 = i 1 , . . . , X k = i k ), the probability that our algorithm picks a sequence of rows indexed i 1 , i 2 , . . . , i k in that order is equal to Otherwise, the probability is zero because in the execution of the algorithm, b i = 0 for some step t. This proves the correctness of our algorithm. Given that one can compute the characteristic polynomial of an m-by-m matrix in O(m ω log m) (see Section 3), our outline can be implemented with the following count of arithmetic operations: for every t and i, O(m 2 n) to compute C i C T i , O(m 2 n + m ω log m) in total for p i . Thus, volume sampling in O(km 3 n + km ω+1 log m). Here we present the first (faster) subroutine for computing the marginal probabilities p i 's within the volume sampling algorithm outlined in Section 4. The two main ideas behind this subroutine are: (1) We can work with Assuming m ≥ n, this saves on running time. ( 2) Each C i is a rank-1 update of B and therefore, once we have B T B, it can be used to compute all C T i C i efficiently. Input: B ∈ R m×n . Output: p 1 , p 2 , . . . , p m . For i = 1 to m do: 1. Compute the matrix C T i C i ∈ R n×n by the following formula 2. Compute the characteristic polynomial of C T i C i and output Proposition 16. For any given B ∈ R m×n , the Algorithm 2 above computes p 1 , . . . , p m in O(mn ω log n) arithmetic operations. Proof. B T B can be computed in time O(mn 2 ). Observe that since C i is obtained by projecting each row of B orthogonal to b i , and therefore, So once we have B T B, for each i, C T i C i can be computed in time O(n 2 ) and the characteristic polynomial of and hence, the above subroutine results into an O(kmn ω log n) time algorithm for volume sampling. Theorem 17 (same as Theorem 7). The randomized algorithm given by the combination of the algorithm outlined in Algorithm 1 with Algorithm 2 as its subroutine, when given a matrix A ∈ R m×n and an integer 1 ≤ k ≤ rank(A), outputs a random k-subset of the rows of A according to volume sampling, using O(kmn ω log n) arithmetic operations. Proof. The proof follows by combining Proposition 15 and Proposition 16, and since we compute all the p i 's simultaneously in each round in O(mn ω log n) arithmetic operations, the total number of arithmetic operations is O(kmn ω log n). Taking further the idea that each C i is a rank-1 update of B, we can give a faster algorithm based on the singular value decomposition of B. Given the singular value decomposition of a matrix and using the matrix determinant lemma (Lemma 13), one can give a precise formula for how the characteristic polynomial changes under a rank-1 update. Using this subroutine in the volume sampling algorithm outlined in Section 4 we get an algorithm for nearly-exact volume sampling (depending on the precision of the computed SVD) in time O(kT svd + kmn 2 ), where T svd is the running time of SVD on an m-by-n matrix. 2. Compute the polynomial products ), and 3. For i = 1 to m output: Proposition 18. In the real arithmetic model and given exact U and Σ, using the Algorithm 3 as a subroutine inside Algorithm 1 outlined for volume sampling, we get an algorithm for volume sampling. If T svd is the running time for computing the singular value decomposition of m-by-n matrices, the algorithm runs in time O(kT svd + kmn 2 ). Proof. Using the matrix determinant lemma (Lemma 13), the characteristic polynomial of C i C T i can be written as Thus, Once we have the singular value decomposition of B, f (x) and g j (x) can all be computed in time O(n 2 ) using polynomial products. This is because there are at most n non-zero σ i 's. Thus, f (x) and all the g j (x) for 1 ≤ j ≤ m can be computed in time O(mn 2 ) and then using the above formula we get c m-k+t (C i C T i ). Magen and Zouzias [17] showed that the random projection lemma of Johnson and Lindenstrauss can be generalized to preserve volumes of subsets after embedding. Here is a restatement of Theorem 1 of [17] using O(ǫ/k) instead of ǫ in their original statement. and there is a mapping f : R n → R d such that for all S ⊆ [m] such that |S| ≤ k, where à ∈ R m×d has its i-th row as f (a i ). Moreover, f is a linear mapping given by multiplication with a random n by d matrix with i.i.d. Gaussian entries, so computing à takes time O(mnd). Theorem 20 (same as Theorem 9). Using random projection for dimensionality reduction, the polynomial time algorithm for volume sampling mentioned in Theorem 7 (i.e., Algorithm 1 with Algorithm 2 as its subroutine), gives (1 + ǫ)-approximate volume sampling, using arithmetic operations. Proof. Using Theorem 19 and doing volume sampling of k-subsets of rows from à gives (1 + ǫ)approximation to the volume sampling of k-subsets of rows from A. This can be done in two steps: first, we compute à using matrix multiplication in time O(mnd) and second, we do volume sampling on à using the algorithm from Subsection 4.1. Overall, it takes time O(mnd + kmd ω log d), which is equal to Moreover, this can be implemented using only one pass over the matrix A with extra space m log m• k 2 /ǫ 2 . Our derandomized row-subset selection algorithm is based on a derandomization of the volume sampling algorithm in Section 4, using the method of conditional expectations. Again, it may be easier to consider volume sampling extended to random k-tuples (X 1 , . . . , X k ) where From Theorem 2 we know that where the expectation is over (X 1 , . . . , X k ). Let us consider i 1 , . . . , i t-1 for which Pr (X 1 = i 1 , . . . , X t-1 = i t-1 ) > 0. Let S = {i 1 , . . . , i t-1 } and look at the conditional expectation. The following lemma shows that these conditional expectations have an easy interpretation in terms of the coefficients of certain characteristic polynomials, and hence can be computed efficiently. Lemma 21. Let (i 1 , . . . , i t-1 ) ∈ [m] t-1 be such that Pr (X 1 = i 1 , . . . , X t-1 = i t-1 ) > 0 for a random k-tuple (X 1 , X 2 , . . . , X k ) from extended volume sampling. Let S = {i 1 , . . . , i t-1 } and B = A-π S (A). Then Proof. Knowing the above lemma, it is easy to derandomize our algorithm outlined for volume sampling. In each step, we just compute the new conditional expectations for each additional i, and finally pick the i that minimizes the conditional expectation. Input: a matrix A ∈ R m×n and 1 ≤ k ≤ rank(A). Output: a subset S of k rows of A with the guarantee 1. Initialize S ← ∅ and B ← A. For t = 1 to k do: Theorem 22 (same as Theorem 8). Deterministic Algorithm 4, when given a matrix A ∈ R m×n and an integer 1 ≤ k ≤ rank(A), outputs a k-subset S of the rows of A, using O(kmn ω log n) arithmetic operations, such that Proof. By applying Lemma 21 to S ∪ {i} instead of S, as C i = Bπ {i} (B) = Aπ S∪{i} (A), we see that the step t of our algorithm picks i that minimizes . The correctness of our algorithm follows immediately from observing that in each step t, and that we started with F . The guarantee for spectral norm follows immediately from our guarantee in the Frobenius norm, just using properties of norms and the fact that rank( 2 . Moreover, this algorithm runs in time O(kmn ω log n) if we use the subroutine in Subsection 4.1 to compute the characteristic polynomial of C i C T i using that of C T i C i . 6 Lower bound for rank-1 spectral approximation using one row Here we show a lower bound for row/column-subset selection. We prove that there is a matrix A ∈ R n×(n+1) such that using the span of any single row of it, we can get only Ω( √ n)-approximation in the spectral norm for the nearest rank-1 matrix to A. This can be generalized to a similar Ω( √ n) lower bound for general k by using a matrix with k block-diagonal copies of A. Theorem 23 (same as Theorem 10). There exists a matrix A ∈ R n×(n+1) such that Let B be the best rank-1 approximation to A whose rows lie in the span of (1, ǫ, 0, . . . , 0) (or for that matter, any fixed row of A). Then, we want to show that We first compute the singular values of A, i.e., the positive square roots of the eigenvalue of AA T ∈ R n×n . . (1, 1, . . . , 1) is an eigenvector of AA T with eigenvalue n + ǫ 2 . Thus, σ 1 (A) = √ n + ǫ 2 . Observe that, by symmetry, all other singular values of A must be equal, i.e., σ 2 (A) = σ 3 (A) = . . . = σ n (A). However, Therefore, A -A 1 2 = σ 2 (A) = ǫ. Now denote the i-th row of A by a i . By definition, the i-th row of B is the projection of a i onto span (a 1 ). We are interested in the singular values of A -B. For i ≥ 2: -ǫ 1 + ǫ 2 , 0, ǫ (i + 1)-th coord. , 0, . . . , 0    . Thus, (A -B)(A -B) T ∈ R n×n can be written as Again, (0, 1, 1, . . . , 1) is the top eigenvector of (A -B)(A -B) T and using this we get, Therefore, We analyzed efficient algorithms for volume sampling that can be used for row/column subset selection. Here are some ideas for future investigation suggested by this work: • It would be interesting to explore how these algorithmic ideas are related to determinantal sampling [16,13] and, in particular, the generation of random spanning trees. • Find practical counterparts of the algorithms discussed here. In particular, we do not analyze the numerical stability of our algorithms. • Is there an efficient algorithm for volume sampling based on random walks? This question is inspired by MCMC as well as random walk algorithms for the generation of random spanning trees.