Faster Subset Selection for Matrices and Applications

Introduction

Problem 1 occurs in numerous situations: column-based low-rank matrix approximation ; feature selection in $k$-means clustering ; optimal experiment design ; multipoint boundary value problems ; sparse solutions to least-squares regression ; sensor selection in a wireless network ; rank-deficient linear least squares , and rank-deficient non-linear least squares , to name just a few. We discuss some of these situations in Section 6.

However, our initial motivation for investigating Problem 1 was our observation that the combinatorial problem of finding a low-stretch spanning tree in an undirected graph corresponds to the Frobenius norm version of Problem 1. This connection is new and might be of independent interest.

We study three aspects of Problem 1: algorithms, lower bounds, and applications. We now summarize our contributions in each of these aspects.

In Section 3 we describe five different approximation algorithms for Problem 1. We suggest five different algorithms because no single algorithm has the lowest operation count; the choice of the most efficient algorithm depends on the actual values of $m$, $n$ and $k$. Our algorithms are considerably faster than the previously known algorithms, and they achieve the same or tighter approximation bounds. Table 1 summarizes the algorithms we propose, as well as previously known algorithms for Problem 1.

Notice, for example, that if $k=m-\alpha$, for some small integer $0<\alpha\leq 0.9(m-n+1)$, then the approximation bound is $1+10\alpha(m-n+1)^{-1}$.

Algorithm 2 in Corollary 8 is designed for the spectral norm case ($\xi=2$). It’s operation count is $O\left(mn^{2}+mn(m-k)\right)$ as well. It finds a subset ${\cal S}$ of cardinality $k$ such that

Similarly, if, for example, $k=n+1+\beta$, for some integer $\beta$ close to $m$ with $0<\beta<m-n+1$, then the approximation bound is $1+n+n(m-n-1)\beta^{-1}$.

This algorithm is inspired by recent results on approximate decompositions of the identity . Notice that, for example, if $k=\Theta(n)$, the approximation bound is $1+O(m/k)$.

Algorithm 5 finds a subset S of cardinality $k=n$ such that

for any $\eta>0$ chosen by the user. This bound is deterministic but the bound on the number of operations is probabilistic. Specifically, for any $0<\delta<1$, we show that with probability $1-\delta$, the operation count is $O\left(mn^{3}\log\delta^{-1}\log^{-1}{(1+\eta)}\right)$.

Our volume-sampling-based algorithm for the subset selection problem can be viewed as a complementary result to the volume-sampling-based algorithms designed before for low-rank matrix approximation . In low-rank matrix approximation, the subspace spanned by the columns that are selected by volume sampling contains a rank $k$ matrix that approximates the best rank $k$ matrix computed via the SVD; in our case, the objective is different but we show that volume sampling gives useful results as well.

1.2 Lower Bounds

1.3 Applications

In Section 5, we study the connection between low-stretch spanning trees and subset selection. Using a result by Spielman and Woo , we prove that the stretch of any tree in an undirected graph equals the Frobenius norm squared of the pseudo-inverse of the sampled matrix that arises by sampling columns from an orthonormal matrix which is a basis for the row space of the so-called node-by-edge incidence matrix of the graph. This incidence matrix contains as many columns as edges in the graph; so, sampling columns from this matrix corresponds to sampling edges from the graph. We then use this reduction to develop novel algorithms for constructing spanning trees with low stretch in undirected graphs. Unfortunately, our algorithms are worse than the available state-of-the-art . We believe, however, that the connection is interesting and might be useful to shed new light on the combinatorial problem of finding a low stretch spanning tree in an undirected graph.

In Section 6 we use the subset selection algorithms of this paper to design novel algorithms for three other problems involving sub-sampling: column-based low-rank matrix reconstruction, sparse solution of least-squares problems, and feature selection in $k$-means clustering.

2 Related Work

We now provide a comprehensive summary of known results regarding Problem 1 and we comment on two related subset selection problems studied in the literature.

and so on, until $m-k$ columns are removed.

A straightforward implementation of this idea requires $O(mn^{3}(m-k))$ operations. However, one can use the Sherman-Morrison formula for rank one updates to the inverse of a matrix and improve the operation count to $O(n^{3}+mn^{2}(m-k))$.

2.2 Rank Revealing Factorizations

The subset selection problem that we study in this paper has deep connections, which we do not explain in detail, with the so-called Rank-Revealing QR (and also see for a summary of available RRQR algorithms) and Rank-Revealing LU factorizations.

By applying the inequality to $i=n$ we have the following bound,

By summing up the bounds on each singular value we get the following bound,

For $f>1$ and $k=n$, the operation count of this method is $O(mn^{2}\log_{f}m)$.

Rank revealing approaches can only be used to sample $k\leq n$ columns; extending these approaches to sample arbitrary $k\geq n$ columns, which is the focus of this paper, is not obvious.

2.3 Incoherent Subset Selection

2.4 Approximation via Convex Relaxation

2.5 Maximum-volume Subsets

Finally, note that the strong RRQR algorithm of finds a local maximum-volume subset. By local maximum-volume subset, we mean that the volume of the subset found is always bigger than the volume of any subset obtained by interchanging a single column.

2.6 Computational Complexity of Subset Selection

The computational complexity of finding a maximum volume subset was also investigated in the computational geometry literature. The problem is stated differently: finding a large simplex in a V-polytope. NP-hardness was established in , and exponential inapproximability was established in .

2.7 Variants of the Subset Selection Problem

Finally, we mention that all these algorithms have found many applications in numerous problems involving subsampling: least-squares regression ; column-based low-rank matrix approximation ; spectral graph sparsification ; and, dimensionality reduction in clustering .

2.8 Restricted Invertibility

Preliminaries

2 Sampling Columns

3 Singular Value Decomposition

4 Moore-Penrose Pseudo-inverse

5 Column Exchanges and Cramer’s rule

6 Volume Sampling

7 Other Known Results

In addition we use the following two known results.

Algorithms

This section describes an algorithm based on the same greedy removal strategy as in , but it is faster, since it exploits the SVD decomposition of the matrix and the ability to quickly update it. Additionally, our algorithm efficiently detects columns whose removal results in a rank deficient matrix, and avoids removing them (see the discussion in Section 1.2). We prove that our algorithm achieves the same approximation bounds as in . The proof of does not apply to our algorithm, since assumes implicitly that in all the iterations, removing a single column does not result in a rank deficient matrix.

Finding an index $j_{i}$ to remove. We then set ${\cal S}_{i}={\cal S}_{i-1}-\{j_{i}\}$.

Before proceeding to the analysis of Algorithm 1, we discuss two numerical stability issues that affect an actual implementation of Algorithm 1. Computing the terms in equation (2) might be problematic since the computation can potentially suffer from catastrophic cancellations when $\mbox{}\|{\mathbf{y}}^{(i-1)}_{r}\|_{2}\approx 1$. However, to find the minimizer we need to do only comparisons. That is, we need to be able to determine for two indices $g,h\in{{\cal S}_{i-1}}$ whether

or not. It is easy to verify that provided $\mbox{}\|{\mathbf{y}}^{(i-1)}_{g}\|_{2}<1$ and $\mbox{}\|{\mathbf{y}}^{(i-1)}_{h}\|_{2}<1$, the last equation holds if and only if

The last equation does not do any subtraction, so it does not suffer from catastrophic cancellations.

Another issue with equation (2) is that an index $h\in{{\cal S}_{i-1}}$ is a candidate minimizer only if $\mbox{}\|{\mathbf{y}}^{(i-1)}_{h}\|_{2}<1$. Under inexact arithmetic that will always be the case, even if removing the column results in a rank deficient system. This issue can be solved by replacing the test $\mbox{}\|{\mathbf{y}}^{(i-1)}_{h}\|_{2}<1$ with $\mbox{}\|{\mathbf{y}}^{(i-1)}_{h}\|_{2}<1-\tau$ for some small threshold $\tau$.

Before proceeding with the proof we state an auxiliary lemma. However, we defer the proof to Section 3.5 since this Lemma is a corollary of a Theorem that appears in that section.

So, we prove only the Frobenius norm bound.

First, we argue that for any $r\in{\cal S}_{i-1}$ the matrix

is singular. That can hold only if $\mbox{}\|{\mathbf{y}}^{(i-1)}_{r}\|_{2}=1$. Therefore, comparing the norm of ${\mathbf{y}}^{(i-1)}_{r}$ with $1$ is an efficient way (once we have an SVD) under exact arithmetic to detect if

is singular. This justifies the restriction $\mbox{}\|{\mathbf{y}}^{(i-1)}_{r}\|_{2}<1$ in equation (2).

We proceed with some calculations. Fix an index $r\in{\cal S}_{i-1}$. If

These calculations, alongside the observation that

We now use this fact to establish the Frobenius norm approximation bound (recall that the spectral norm bound follows immediatly from the Frobenius norm bound). We will show using induction that

Since our algorithm returns ${\cal S}_{m-k}$ the claim follows from (3).

Equation (3) trivially holds for $i=0$. Assume it holds for $i-1$. We now show it holds for $i$. Note that the cardinality of ${\cal S}_{i-1}$ is $m-i+1$. Lemma 7 ensures that there exists a subset ${\cal T}_{i}\subset{\cal S}_{i-1}$ of cardinality $m-i$ such that

2 Deterministic Greedy Removal (spectral norm)

Let ${\cal S}\subseteq[m]$ be the set found by the algorithm, and let $\bar{\cal S}=[m]-{\cal S}$. Corollary 2 of asserts that

3 Deterministic Greedy Selection

The algorithm of this section builds the set ${\cal S}$ by iteratively adding columns to it, after starting with the empty set. It uses a deterministic algorithm presented in , which is, in turn, a generalization of an algorithm from . In particular, we use Lemma 10 from .

We refer the reader to for the full description of the algorithm. Lemma 9 implies that one can sample from two different set of vectors ${\cal V}=\{{\mathbf{v}}_{1},\ldots,{\mathbf{v}}_{m}\}$ and ${\cal U}=\{{\mathbf{u}}_{1},\ldots,{\mathbf{u}}_{m}\}$, and control simultaneously the smallest singular value of the matrix formed from the sampled vectors from the first set, and the largest singular value of the matrix formed from the sampled vectors from the second set.

We now present the analysis of Algorithm 3.

We first prove the approximation bound, and then bound the number of operations.

However, $\sum_{i=1}^{m}s_{i}{\mathbf{e}}_{i}{\mathbf{e}}_{i}^{\textsc{T}}=\hbox{\rm diag}(s_{1},\dots,s_{m})\in\R^{m\times m}$, a diagonal matrix containing the weights $s_{i}$’s in its main diagonal; so, $\max_{i}s_{i}\leq\left(1+\sqrt{{m\over k}}\right)^{2}$. Lemma 9 also guarantees that

4 Randomized Selection

for $i=1,\dots,m$. The set ${\cal S}$ is formed by non-uniformly, and independently, sampling $k$ numbers from ${1,\dots,m}$ with replacement. In each trial, $i$ is sampled with probability

We now present the analysis of Algorithm 4.

We first prove the approximation bound, and then bound the number of operations.

$(a)$ follows by replacing the values taken by the vector ${\mathbf{x}}$. $(b)$ follows by replacing the value for the probabilities $p_{i}$’s. $(c)$ follows by the fact that $\tau_{j}\geq\mbox{}\|{\mathbf{y}}_{j}\|_{2}^{2}$, for all $j=1,...,m$. $(d)$ follows by replacing the value for the parameters $\tau_{i}$’s. $(e)$ follows by simple algebra.

We are now ready to apply Lemma 11 for the random vector ${\mathbf{y}}$ described above. An immediate application of this Lemma ($M=\sqrt{2n}$, $\epsilon=1/2$) and our bound on $k$ give that with probability at least $1-\delta$,

which we proved in Eqn. $(f)$ in the previous calculations. $(f)$ follows by simple algebra.

To conclude the analysis of the approximation bound, notice that at the end of Theorem 10, we implicitly proved that

5 Volume based bounds and algorithms

We start with Lemma 13, which establishes the connection between determinants and subset selection.

We first prove the equality in the Lemma ($\cal S$ has cardinality k unless otherwise stated),

5.2 Random Subsets Chosen via Volume Sampling

We now conclude the first part of the proof as follows,

The above bound is obtained as follows. Recall the two inequalities proved above for any ${\cal T}$:

We can now prove the following corollary, which was previously stated as Lemma 7.

If there exists a subset ${\cal S}$ of columns of cardinality $k$ such that these columns are linearly dependent then (4) might be a strict inequality. For example, let

There is only one set of cardinality $n$ that has positive volume (i.e., the set of columns is full rank): ${\cal T}=[n]$. Since this is the only set with positive probability we have

5.3 Volume Sampling Subset Selection

We now present the analysis of Algorithm 5.

For every $0<\delta<1$, the algorithm will terminate after $O\left(mn^{3}\log{(1/\delta)}/\log{(1+\eta)}\right)$ operations with probability of at least $1-\delta$.

We first prove the approximation bound, and then bound the number of operations.

Each iteration (line 3) takes $O(n^{3}m)$. Combining this with the analysis in the previous paragraph reveals that for any $0<\delta<1$, Algorithm 5 will finish after $O\left(mn^{3}\log{(1/\delta)}/\log{(1+\eta)}\right)$ operations with probability of at least $1-\delta$.

Lower Bounds

We first state two known results that will be used in our proof.

As $\alpha\rightarrow 0$, the bound in the above theorem is $m/k-1$. If $k=(1+\Omega(1))n$ then the upper bound of the deterministic algorithm of Theorem 10 asymptotically matches this lower bound. The upper bounds of the algorithms in the Theorems 6, 12, and 16 and Corollary 8 are - asymptotically - slightly worse. However, if $k=(1+o(1))n$ there is a gap between the lower bound and the best upper bound.

2 Lower bound for the Frobenius norm version of the subset selection problem

To prove the bound we use Theorem 19 and Lemma 7 from .

As $\alpha\rightarrow 0$ and $k=O\left(n\right)$, this bound is $m/k-O(1)$. If $k=(1+\Omega(1))n$ the Frobenius norm bounds in Theorems 6 and 10 asymptotically match this lower bound. However, if $k=(1+o(1))n$ there is a gap between the lower bound and the best upper bound. There is also a gap when $k=\omega(n)$. We believe that the gap for $k=\omega(n)$ is the result of looseness in the lower bound, but we were unable to prove a tighter bound than Theorem 21.

Low-stretch Spanning Trees and Subset Selection

The problem of finding a low-stretch spanning tree is the problem of finding a spanning tree $T$ of $G$ such that $\text{\rm St}_{T}(G)$ is minimized, among all possible spanning trees of $G$. Let $\text{\rm St}(n)=\max_{G\in G_{n}}\min_{T}\text{\rm St}_{T}(G)$, where $G_{n}$ is the family of graphs with $n$ vertices. The following bounds are known: $\text{\rm St}(n)=\Omega(m\log n)$ ; $\text{\rm St}(n)=O(m\log n\cdot\log\log n\cdot(\log\log\log n)^{3})$ . In this section we show that finding a low-stretch spanning tree is in fact an instance of the Frobenius norm version of Problem 1.

Finding a low stretch spanning tree has quite a few uses. One important application is the solution of symmetric diagonally dominant (SDD) linear systems of equations. Boman and Hendrickson were the first to suggest the use of low-stretch spanning trees to build preconditioners for SDD matrices. Spielman and Teng later showed how to use low stetch spanning trees to solve SDD systems using a nearly linear amount of operations. The currently most efficient algorithm for solving SDD systems uses a low stretch spanning tree as well. One of the many obstacles in generalizing these algorithms for wider classes of matrices (e.g., finite-element matrices) is the lack of an equivalent concept, like the stretch, for such matrices. By studying the purely linear-algebraic nature of the low-stretch spanning tree problem (i.e. the Frobenius norm version of Problem 1), our hope is to glean new insights on how to generalize the concept of low-stretch trees, or to substitute it with something else.

Other applications of low-stretch spanning trees include: Alon-Karp-Peleg-West game, MCT approximation and message-passing model. See for details.

Next, we show that finding a low-stretch spanning tree is an instance of subset selection. We first relate graphs to matrices.

We conducted some simple experiments with our greedy removal algorithm. In the first experiment, we used greedy removal to generate a spanning tree $T_{n}$ of the complete graph $K_{n}$ on $n$ with equal weights vertices, for $n=10,11,\dots,50$. We then computed the stretch of $T_{n}$. We found that $\text{\rm St}_{T_{n}}(K_{n})\approx 0.6m\log^{2}n$. We then repeated this experiment with random weights on the edges of $K_{n}$. We found that in almost all runs, $\text{\rm St}_{T_{n}}(K_{n})\approx 0.3m\log^{2}n$. These values are much better than our theoretical bounds, and are closer to what it is possible to find using state-of-the-art algorithms for low-stretch trees. These experiments, although far from exhaustive, suggest that our theoretical worst-case upper bounds for greedy removal are rather pessimistic for the matrices relevant to finding a low-stretch spanning tree.

2 Maximum weight spanning trees and maximum volume subsets

The subset of columns ${\cal S}$ that maximizes the volume also maximizes $\prod_{e\in H({\cal S})}w(e)$. This trivially implies that the corresponding tree is a maximum weight spanning tree. So, for edge-incidence matrices one can use an efficient maximum weight spanning tree algorithm to find the maximum volume subset of columns efficiently. The bound we obtain is $\text{\rm St}_{T}(G)<(n-1)(m-n+2)$. We are unaware of any other analysis of the stretch of a maximum weight spanning tree, but this bound can be easily proven using much simpler arguments.

3 Low-stretch spanning trees via volume sampling

Let $G$ be a weighted undirected connected graph, and let ${\cal T}$ be a random spanning tree, where tree $T$ is sampled with relative probability $\prod_{e\in T}w(e)$. Then,

One can use VolumeSample from to generate such a spanning tree in $O(n^{3}m)$ operations. However, the problem of generating a sample from $\Gamma(G)$ is a well studied problem, and there exists algorithms that can generate a random spanning tree faster than $O(n^{3}m)$. See for a short review.

4 Towards better bounds for low-stretch spanning trees

State-of-the-art algorithms for finding low stretch spanning trees attain theoretical worst-case bounds that are better than the ones we obtain for a general matrix. We now provide a preliminary explanation for this gap.

Other Applications

2 Sparse Solutions to Least-squares Regression Problems

However, sometimes a different regularization is sought: requiring the solution vector to be sparse. That is, we are interested in constructing a vector ${\mathbf{x}}_{k}\in\R^{m}$ that has at most $k$ non-zeros, for some $k$. Since truncated SVD is arguably the most natural regularizer, it makes sense to compare ${\mathbf{x}}_{k}$ to ${\mathbf{x}}_{svd(r)}$ for some $r\leq k$. More specifically, we are interested in bounds of the form,

The idea of obtaining sparse solutions with approximation bounds of the above type can be traced to . Currently, the best deterministic method is in ($k>r$) with

We omit the proof since it follows immediately by combining Lemma 3 from with Corollary 8 in our paper. The algorithm is deterministic and the operation count is $O(dm\min\{d,m\}+mr(m-k))$.

3 Feature Selection in k𝑘k-Means Clustering

The deterministic algorithm of Corollary 8 can also be used for deterministic feature selection in $k$-means clustering. We refer the reader to for an introduction to this problem. Theorem 4 of gives such a polynomial-time deterministic unsupervised feature selection algorithm, which selects features from the data and then rescales them. Using Corollary 8, one can design a deterministic unsupervised feature selection algorithm without rescaling. We omit the details, since the algorithm is similar to the one described for sparse least squares, and the analysis is a combination of Lemma 10 from with Corollary 8. The approximation bound that is obtained is comparable to the bound in .

Open Problems and Future Directions

Several interesting questions remain unanswered and we leave them for future investigation. First, is the Frobenius-norm version of Problem 1 NP-hard? Second, is it possible to close the existing gaps between lower and upper bounds for Problem 1? Third, is it possible to extend the Strong Rank Revealing QR method of to sample arbitrary $k\geq n$ columns? Fourth, is it possible to extend the polynomial implementations of volume sampling in to sample arbitrary number of columns from short-fat matrices? Finally, is it possible to derandomize the algorithm of Theorem 16?

Acknowledgements

We would like to thank the two anonymous referees and the editor for their numerous comments and suggestions; Ioannis Koutis for bringing to our attention; Petros Drineas, Frank De Hoog, Ilse Ipsen, Sivan Toledo, and Mark Tygert for many useful discussions and suggestions on a preliminary draft of this work; and Anastasios Zouzias for pointing out the connection between the restricted invertibility line of research and ours.

The authors acknowledge the support from XDATA program of the Defense Advanced Research Projects Agency (DARPA), administered through Air Force Research Laboratory contract FA8750-12-C-0323.

References