Efficient volume sampling for row/column subset selection

Introduction

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 in the context of low-rank approximation of matrices. It is equivalent to sampling $k$-subsets of $\{1,\dotsc,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 . 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 highlights the limitations of PCA and gives experimental data on practical applications of low-rank approximation based on row/column-subset selection.

Row/column-subset selection and volume sampling

In other words, the rows of $A_{k}$ are projections of the rows of $A$ onto $\operatorname{span}\left(v_{i}\;:\;1\leq i\leq k\right)$. 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$.

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 . The first among these is the squared-length sampling of rows introduced by Frieze, Kannan and Vempala . Another sampling scheme due to Drineas, Mahoney and Muthukrishnan uses the singular values and singular vectors to decide the sampling probabilities. Later Deshpande, Rademacher, Vempala and Wang introduced volume sampling as a generalization of squared-length sampling.

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 . 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$.

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.

However, no efficient algorithm was known for volume sampling prior to this work. An algorithm mentioned in Deshpande and Vempala 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 ).

Boutsidis, Drineas and Mahoney 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 , and the previous best algorithmic result for row-subset selection in the spectral norm case was given by a result of Gu and Eisenstat on strong rank-revealing QR decompositions. The following theorem is a direct consequence of as pointed out in .

Moreover, this subset $S$ can be found in time $O\left((m+n\log_{f}n)n^{2}\right)$.

In the context of volume sampling, it is interesting to note that Pan has used an idea of picking submatrices of locally maximum volume (or determinants) for rank-revealing matrix decompositions. We refer the reader to for details.

The results of Goreinov, Tyrtyshnikov and Zamarashkin on pseudo-skeleton approximations of matrices look at submatrices of maximum determinants as good candidates for row/column-subset selection.

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.

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 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 , 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 $m$-by-$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 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:

This improves the $O(k\sqrt{\log k})$-approximation of Boutsidis, Drineas and Mahoney for the Frobenius norm case and matches the lower bound shown in Theorem 3 due to .

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+\epsilon)$-approximate volume sampling, using

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$.

Preliminaries and notation

Throughout the paper we assume $m\geq 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.

can be generalized into the following lemma.

where $\sigma_{1},\dotsc,\sigma_{m}$ are the singular values of $A$, i.e., eigenvalues of $AA^{T}$, and

is the characteristic polynomial of $AA^{T}$. Using $\operatorname{det}\left(xI-AA^{T}\right)=x^{m-n}\operatorname{det}\left(xI-A^{T}A\right)$, we can alternatively use $c_{m-k}(AA^{T})=c_{n-k}(A^{T}A)$ in the above formula, for $k\leq n$.

Let $\omega$ 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^{\omega}\log n)$ arithmetic operations [2, Section 16.6].

Here is another lemma that we will need about dividing determinants into products of two determinants.

To show this, we consider two cases. If $C_{S}C_{S}^{T}$ is singular, then both sides of the equality are zero. If $C_{S}C_{S}^{T}$ is invertible, then we can perform block Gaussian elimination and write

applied to $\begin{pmatrix}E&F\\ G&H\end{pmatrix}=C$. 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

so that $D_{T}=C_{T}-C_{T}C_{S}^{T}(C_{S}C_{S}^{T})^{-1}C_{S}$, and

Finally, a well-known lemma about how the determinant of a matrix changes under a rank-$1$ update.

Efficient volume sampling algorithms

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},\dotsc,X_{k})$ instead of $k$-subsets, where each of the $k!$ permutations of a $k$-subset is equally likely, i.e., for any $(i_{1},i_{2},\dotsc,i_{k})\in[m]^{k}$,

Then the marginal probabilities ${\sf Pr}\left(X_{t}=i\;|\;X_{1}=i_{1},\dotsc,X_{t-1}=i_{t-1}\right)$ have the following interpretation in terms of the coefficients of certain characteristic polynomials.

Let $(i_{1},\dotsc,i_{t-1})\in[m]^{t-1}$ such that ${\sf Pr}\left(X_{1}=i_{1},\dotsc,X_{t-1}=i_{t-1}\right)>0$, for a random $k$-tuple $(X_{1},X_{2},\dotsc,X_{k})$ from the extended volume sampling over $k$-tuples. Let $S=\{i_{1},\dotsc,i_{t-1}\}$, $B=A-\pi_{S}(A)$ and $C_{i}=B-\pi_{\{i\}}(B)=A-\pi_{S\cup\{i\}}(A)$. Then,

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.

Output: a subset $S$ of $k$ rows of $A$ picked with probability proportional to $\operatorname{det}\left(A_{S}A_{S}^{T}\right)$.

Initialize $S\leftarrow\emptyset$ and $B\leftarrow A$. For $t=1$ to $k$ do:

where $C_{i}=B-\pi_{\{i\}}(B)$ is a matrix obtained by projecting each row of $B$ orthogonal to $b_{i}$.

Pick $i$ with probability proportional to $p_{i}$. Let $S\leftarrow S\cup\{i\}$ and $B\leftarrow C_{i}$.

Now we show the correctness of the algorithm:

The probability that our volume sampling algorithm outlined above picks a $k$-subset $S$ is proportional to $\operatorname{det}\left(A_{S}A_{S}^{T}\right)$. This algorithm can be implemented with a cost of $O(km^{3}n+km^{\omega+1}\log m)$ arithmetic operations.

By Lemma 14, for any $i_{1},i_{2},\dotsc,i_{k}$ such that ${\sf Pr}\left(X_{1}=i_{1},\dotsc,X_{k}=i_{k}\right)$, the probability that our algorithm picks a sequence of rows indexed $i_{1},i_{2},\dotsc,i_{k}$ in that order is equal to

Otherwise, the probability is zero because in the execution of the algorithm, $\left\|b_{i}\right\|=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^{\omega}\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_{i}^{T}$, $O(m^{2}n+m^{\omega}\log m)$ in total for $p_{i}$. Thus, volume sampling in $O(km^{3}n+km^{\omega+1}\log m)$. ∎

First subroutine for marginal probabilities

Compute the characteristic polynomial of $C_{i}^{T}C_{i}$ and output

$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}$,

So once we have $B^{T}B$, for each $i$, $C_{i}^{T}C_{i}$ can be computed in time $O(n^{2})$ and the characteristic polynomial of $C_{i}^{T}C_{i}$ can be computed in time $O(n^{\omega}\log n)$ [2, Section 16.6]. By Lemma 11, $c_{m-k+t}(C_{i}C_{i}^{T})=c_{n-k+t}(C_{i}^{T}C_{i})$ and hence, the above subroutine results into an $O(kmn^{\omega}\log n)$ time algorithm for volume sampling. ∎

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^{\omega}\log n)$ arithmetic operations, the total number of arithmetic operations is $O(kmn^{\omega}\log n)$. ∎

2 Efficient volume sampling using SVD

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.

Second subroutine for marginal probabilities.

In the real arithmetic model and given exact $U$ and $\Sigma$, 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})$.

Using the matrix determinant lemma (Lemma 13), the characteristic polynomial of $C_{i}C_{i}^{T}$ can be written as

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 $\sigma_{i}$’s. Thus, $f(x)$ and all the $g_{j}(x)$ for $1\leq j\leq m$ can be computed in time $O(mn^{2})$ and then using the above formula we get $c_{m-k+t}(C_{i}C_{i}^{T})$. ∎

3 Approximate volume sampling in nearly linear time

Magen and Zouzias 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 using $O(\epsilon/k)$ instead of $\epsilon$ in their original statement.

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+\epsilon)$-approximate volume sampling, using

Moreover, this can be implemented using only one pass over the matrix $A$ with extra space $m\log m\cdot k^{2}/\epsilon^{2}$. ∎

Derandomized row/column-subset selection

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},\dotsc,X_{k})$ where

where the expectation is over $(X_{1},\dotsc,X_{k})$.

Let us consider $i_{1},\dotsc,i_{t-1}$ for which ${\sf Pr}\left(X_{1}=i_{1},\dotsc,X_{t-1}=i_{t-1}\right)>0$. Let $S=\{i_{1},\dotsc,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.

Let $(i_{1},\dotsc,i_{t-1})\in[m]^{t-1}$ be such that ${\sf Pr}\left(X_{1}=i_{1},\dotsc,X_{t-1}=i_{t-1}\right)>0$ for a random $k$-tuple $(X_{1},X_{2},\dotsc,X_{k})$ from extended volume sampling. Let $S=\{i_{1},\dotsc,i_{t-1}\}$ and $B=A-\pi_{S}(A)$. Then

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.

Output: a subset $S$ of $k$ rows of $A$ with the guarantee

Initialize $S\leftarrow\emptyset$ and $B\leftarrow A$. For $t=1$ to $k$ do:

For $i=1$ to $m$ do: compute $c_{n-k+t-1}(C_{i}^{T}C_{i})$ and $c_{n-k+t}(C_{i}^{T}C_{i})$, where $C_{i}=B-\pi_{\{i\}}(B)$ is the matrix obtained by projecting each row of $B$ orthogonal to $b_{i}$.

Pick $i$ that minimizes $\left|c_{n-k+t-1}(C_{i}^{T}C_{i})\right|/\left|c_{n-k+t}(C_{i}^{T}C_{i})\right|$. Let $S\leftarrow S\cup\{i\}$ and $B\leftarrow C_{i}$.

By applying Lemma 21 to $S\cup\{i\}$ instead of $S$, as $C_{i}=B-\pi_{\{i\}}(B)=A-\pi_{S\cup\{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$,

The guarantee for spectral norm follows immediately from our guarantee in the Frobenius norm, just using properties of norms and the fact that $\operatorname{rank}(A-A_{k})\leq n-k$:

Moreover, this algorithm runs in time $O(kmn^{\omega}\log n)$ if we use the subroutine in Subsection 4.1 to compute the characteristic polynomial of $C_{i}C_{i}^{T}$ using that of $C_{i}^{T}C_{i}$. ∎

Lower bound for rank-111 spectral approximation using one row

Let $B$ be the best rank-$1$ approximation to $A$ whose rows lie in the span of $(1,\epsilon,0,\dotsc,0)$ (or for that matter, any fixed row of $A$). Then, we want to show that

$(1,1,\dotsc,1)$ is an eigenvector of $AA^{T}$ with eigenvalue $n+\epsilon^{2}$. Thus, $\sigma_{1}(A)=\sqrt{n+\epsilon^{2}}$. Observe that, by symmetry, all other singular values of $A$ must be equal, i.e., $\sigma_{2}(A)=\sigma_{3}(A)=\dotsc=\sigma_{n}(A)$. However,

Therefore, $\left\|A-A_{1}\right\|_{2}=\sigma_{2}(A)=\epsilon$.

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 $\operatorname{span}\left(a_{1}\right)$. We are interested in the singular values of $A-B$. For $i\geq 2$:

Again, $(0,1,1,\dotsc,1)$ is the top eigenvector of $(A-B)(A-B)^{T}$ and using this we get,

Discussion

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 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.

References