Theory of Computing, Volume 2 (2006), pp. 225–247 http://theoryofcomputing.org

Matrix Approximation and Projective Clustering via Volume Sampling

Amit Deshpande^{∗} Luis Rademacher^{†} Santosh Vempala^{‡}

Grant Wang^{§}

Received: August 11, 2005; revised: September 28, 2006; published: October 12, 2006.

Abstract: Frieze, Kannan, and Vempala (JACM 2004) proved that a small sample of rows of a given matrix A spans the rows of a low-rank approximation D that minimizes ∥A−D∥_{F} within a small additive error, and the sampling can be done efficiently using just two passes over the matrix. In this paper, we generalize this result in two ways. First, we prove that the additive error drops exponentially by iterating the sampling in an adaptive manner (adaptive sampling). Using this result, we give a pass-efficient algorithm for computing a low-rank approximation with reduced additive error. Our second result is that there exist k rows of A whose span contains the rows of a multiplicative (k + 1)-approximation to the best rank-k matrix; moreover, this subset can be found by sampling k-subsets of rows from a natural distribution (volume sampling). Combining volume sampling with adaptive sampling yields the existence of a set of k + k(k + 1)/ϵ rows whose span contains the rows of a multiplicative (1 + ϵ)-approximation. This leads to a PTAS for the following NP-hard

∗Supported by NSF Award CCR-0312339. †Supported by NSF Award CCR-0312339. ‡Supported by NSF Award CCR-0312339 and a Guggenheim Foundation Fellowship. §Supported by NSF Award CCR-0312339.

ACM Classification: G.1.3, F.2.1

AMS Classification: 68W25, 65F15

Key words and phrases: algorithms, matrix approximation, projective clustering.

Authors retain copyright to their papers and grant “Theory of Computing” unlimited rights to publish the paper electronically and in hard copy. Use of the article is permitted as long as the author(s) and the journal are properly acknowledged. For the detailed

copyright statement, see http://theoryofcomputing.org/copyright.html.

\begin{matrix}c \\ ◯ \end{matrix} 2006 Amit Deshpande and Luis Rademacher and Santosh Vempala and Grant Wang

A. Deshpande, L. Rademacher, S. Vempala, G. Wang

projective clustering problem: Given a set P of points in ℝ^{d}, and integers j and k, find subspaces F_{1},...,F_{j}, each of dimension at most k, that minimize ∑_{p∈P} min_{i}d(p,F_{i})^{2}.

1 Introduction

1.1 Motivation

Let the rows of a matrix be points in a high-dimensional space. It is often of interest to find a lowdimensional representation. The subspace spanned by the top k right singular vectors of the matrix is a good choice for many applications. The problem of efficiently finding an approximation to this subspace has received much attention in the past decade [19, 15, 1, 14, 16]. In this paper, we give new algorithms for this problem and show existence of subspaces lying in the span of a small set of rows with better additive approximation as well as multiplicative approximation. At the heart of our analysis are generalizations of previous sampling schemes [19]. We apply these results to the general problem of finding j subspaces, each of dimension at most k, so as to minimize the sum of squared distances of each point to its nearest subspace, a measure of the “error” incurred by this representation; as a result, we obtain the first polynomial-time approximation scheme for this projective clustering problem [5, 32, 6, 3] when j and k are fixed.

The case of j = 1, i. e., finding a single k-dimensional subspace, is an important problem in itself and can be solved efficiently (for j ≥ 2, the problem is NP-hard [30], even for k = 1 [15]). The optimal projection is given by the rank-k matrix A_{k} = AY Y^{T} where the columns of Y are the top k right singular vectors of A and can be computed using the Singular Value Decomposition. Note that among all rank-k matrices D, A_{k}is the one that minimizes ∥A−D∥_{F}^{2} = ∑_{i,j}(A_{ij} −D_{ij})^{2}. The running time of this algorithm, dominated by the SVD computation of an m × n matrix, is O(min{mn^{2},nm^{2}}). Although polynomial, this is still too high for some applications.

For problems on data sets that are too large to store/process in their entirety, one can view the data as a stream of data items arriving in arbitrary order, and the goal is to process a subset chosen judiciously on the fly and then extrapolate from this subset. Motivated by the question of finding a faster algorithm, Frieze et al. [19] showed that any matrix A has k/ϵ rows whose span contains the rows of a rank-k approximation to A within additive error ϵ∥A∥_{F}^{2}. In fact, the subset of rows can be obtained as independent samples from a distribution that depends only on the norms of the rows. (In what follows, A^{(i)} denotes the ith row of A.)

Theorem 1.1 ([19]). Let S be a sample of s rows of an m × n matrix A, each chosen independently from the following distribution: row i is picked with probability

P_{i} = \frac{∥A^{(i)}∥^{2}}{∥A∥_{F}^{2}} .

If s ≥ k/ϵ, then span(S) contains the rows of a matrix Ã_{k}of rank at most k for which

E_{S}(∥A − Ã_{k}∥_{F}^{2}) ≤ ∥A − A_{k}∥_{F}^{2} + ϵ∥A∥_{F}^{2} .

Matrix Approximation and Projective Clustering via Volume Sampling

This can be turned into an efficient algorithm based on sampling [15].^{1} The algorithm makes one pass through A to figure out the sampling distribution and another pass to sample and compute the approximation. Its complexity is O(min{m,n}k^{2}/ϵ^{4}). These results lead us to the following questions:

(1) Can the error be reduced significantly by using multiple passes through the data? (2) Can we get multiplicative (1+ϵ)-approximations? (3) Do these sampling algorithms have any consequences for the general projective clustering problem?

1.2 Our results

We begin with the observation that the additive error term drops exponentially with the number of passes when the sampling is done adaptively. Thus, low-rank approximation is a natural problem for which multiple passes through the data are highly beneficial.

The idea behind the algorithm is quite simple. As an illustrative example, suppose the data consists of points along a 1-dimensional subspace of ℝ^{n} except for one point. The best rank-2 subspace has zero error. However, one round of sampling will most likely miss the point far from the line. So we use a two-round approach. In the first pass, we get a sample from the squared norm distribution. Then we sample again, but adaptively—we sample with probability proportional to the squared distance to the span of the first sample. We call this procedure adaptive sampling. If the lone far-off point is missed in the first pass, it will have a high probability of being chosen in the second pass. The span of the full sample now contains the rows of a good rank-2 approximation. In the theorem below, for a set S of rows of a matrix A, we denote by π_{S}(A) the matrix whose rows are the projection of the rows of A to the span of S.

Theorem 1.2. Let S = S_{1} ∪ ··· ∪ S_{t} be a random sample of rows of an m × n matrix A where, for j = 1,...,t, each set S_{j} is a sample of s rows of A chosen independently from the following distribution: row i is picked with probability

P_{i}^{(j)} = \frac{∥E_{j}^{(i)}∥^{2}}{∥E_{j}∥_{F}^{2}}

where E_{1} = A, E_{j} = A − π_{S_{1}∪··· ∪S_{j−1}} (A). Then for s ≥ k/ϵ, span(S) contains the rows of a matrix Ã_{k}of rank at most k such that

E_{S}(∥A − Ã_{k}∥_{F}^{2}) ≤ \frac{1}{1 − ϵ} ∥A − A_{k}∥_{F}^{2} + ϵ^{t}∥A∥_{F}^{2} .

The proof of Theorem 1.2 is given in Section 2.1. The resulting algorithm, described in Section 2.2 uses 2t passes through the data. Although the sampling distribution is modified t times, the matrix itself is not changed. This is especially significant when A is sparse as the sparsity of the matrix is maintained.

Theorem 1.2 raises the question of whether we can get a multiplicative approximation instead of an additive approximation. To answer this question, we generalize the sampling approach. For a set S of k points in ℝ^{n}, let Δ(S) denote the k-dimensional simplex formed by them along with the origin. We pick a random k-subset S of rows of A with probability proportional to vol(Δ(S))^{2}. This procedure, which we call volume sampling, is a generalization of the earlier sampling approach which picks single rows according to their squared norms.

1Frieze et al. [19] go further to show that there is an s × s submatrix for s = poly(k/ϵ) from which the low-rank approximation can be computed in poly(k,1/ϵ) time in an implicit form.

Theory of Computing, Volume 2 (2006), pp. 225–247 227

Theorem 1.3. Let S be a random k-subset of rows of a given matrix A chosen with probability

P_{S} = \frac{vol(Δ(S))^{2}}{∑_{T:|T |=k}vol(Δ(T))^{2}} .

Then Ã_{k}, the projection of A to the span of S, satisfies

E_{S}(∥A − Ã_{k}∥_{F}^{2}) ≤ (k + 1)∥A − A_{k}∥_{F}^{2} .

We prove this theorem in Section 1.4. Moreover, the factor of k +1 is the best possible for a k-subset (Proposition 3.3). By combining Theorem 1.3 with the adaptive sampling idea from Theorem 1.2, we show that there exist O(k^{2}/ϵ) rows whose span contains the rows of a multiplicative (1 + ϵ)- approximation.

Theorem 1.4. For any m × n matrix A, there exist k + k(k + 1)/ϵ rows whose span contains the rows of a rank-k matrix Ã_{k}such that

∥A − Ã_{k}∥_{F}^{2} ≤ (1 + ϵ) ∥A − A_{k}∥_{F}^{2} .

The existence of a small number of rows containing a good multiplicative approximation is the key ingredient in our last result—a polynomial-time approximation scheme (PTAS) for the general projective clustering problem. This result makes a connection between matrix approximation and projective clustering. A key idea in the matrix approximation work of [14, 15, 19] is that, for any matrix, there is a small subset of its rows whose span contains a good approximation to the row space of the entire matrix. This is similar to the idea of core-sets [2], which have been studied in computing extent measures in computational geometry (and applied to a variant of the projective clustering problem [24]). Roughly speaking, a core-set is a subset of a point-set such that computing the extent measure on the core-set provides an approximation to the extent measure on the entire point set. An extent measure is just a statistic on the point set (for example, the diameter of a point set, or the radius of the minimum enclosing cylinder); typically, it measures the size of the point set or the minimum size of an object that encloses the point set [2].

We state the projective clustering problem using the notation from computational geometry: Let d(p,F) be the orthogonal distance of a point p to a subspace F. Given a set P of n points in ℝ^{d}, find j subspaces F_{1},...,F_{j}, each of dimension k, such that

𝒞(F_{1},...,F_{j}) = ∑_{p∈P} min_{i}d(p,F_{i})^{2}  (1.1)

is minimized. When subspaces are replaced by flats, the case k = 0 corresponds to the “j-means problem” in computational geometry.

Theorem 1.4 suggests an enumerative algorithm. The optimal set of k-dimensional subspaces induces a partition P_{1} ∪ ··· ∪ P_{j} of the given point set. Each set P_{i} contains a subset of size O(k^{2}/ϵ) in whose span lies a multiplicative (1+ϵ)-approximation to the optimal k-dimensional subspace for P_{i}. So we consider all possible combinations of j subsets each of size O(k^{2}/ϵ), and a δ-net of k-dimensional subspaces in the span of each subset. The δ-net depends on the points in each subset and is not just a grid, as is often the case. Each possible combination of subspaces induces a partition and we output the best of these. Since the subset size is bounded (and so is the size of the net), this gives a PTAS for the problem (see Section 4) when j and k are taken to be fixed constants.

Theorem 1.5. Given n points in ℝ^{d} and parameters B and ϵ, in time

d ( \frac{n}{ϵ})^{O(jk^{3}/ϵ)}

we can find a solution to the projective clustering problem which is of cost at most (1 + ϵ)B provided there is a solution of cost B.

1.3 Related work

The work of Frieze et al. [19] and Drineas et al. [15] introduced matrix sampling for fast low-rank approximation. Subsequently, an alternative sampling-based algorithm was given by Achlioptas and Mc-Sherry [1]. That algorithm achieves somewhat different bounds (see [1] for a detailed comparison) using only one pass. It does not seem amenable to the multipass improvements presented here. Bar-Yossef [9] has shown that the bounds of these algorithms for one or two passes are optimal up to polynomial factors in 1/ϵ.

These algorithms can also be viewed in the streaming model of computation [25]. In this model, we do not have random access to data; the data comes as a stream of data items in arbitrary order and we are allowed one or a few sequential passes over the data. Algorithms for the streaming model have been designed for computing frequency moments [7], histograms [22], etc. and have mainly focused on what can be done in one pass. There has been some recent work on what can be done in multiple passes [14, 18]. The “pass-efficient” model of computation was introduced in [25]. Our multipass algorithms fit this model and relate the quality of approximation to the number of passes. Feigenbaum et al. [18] show such a relationship for computing the maximum unweighted matching in bipartite graphs.

The Lanczos method is an iterative algorithm that is used in practice to compute the Singular Value Decomposition [20, 26]. An exponential decrease in an additive error term has also been proven for the Lanczos method under a different notion of additive error ([20, 26]). However, the exponential decrease in error depends on the gap between singular values. In particular, the following is known for the Lanczos method: after k iterations, each approximate singular value θ_{i}^{2} obeys:

θ_{i}^{2} ≥ σ_{i}^{2}− c^{k}C  (1.2)

where σ_{i} is the ith singular value. Both constants c and C depend on the gap between singular values; in particular, c is proportional to (σ_{2}^{2} − σ_{r}^{2})/(2σ_{1}^{2} − σ_{2}^{2} − σ_{r}^{2}) and C = σ_{1}^{2} − σ_{r}^{2}. The guarantee (1.2) can be transformed into an inequality:

∥A − Ã_{k}∥_{F}^{2} ≤ ∥A − A_{k}∥_{F}^{2} + c^{k}C  (1.3)

very similar to Theorem 1.2, but without the multiplicative error term for ∥A−A_{k}∥_{F}^{2}. However, note that when σ_{1}^{2} = σ_{2}^{2}, successive rounds of the Lanczos method fail to converge to σ_{i}^{2}. In the Lanczos method, each iteration can be implemented in one pass over A, whereas our algorithm requires two passes over A in each iteration. Kuczyński and Woźniakowski [27] prove that the Lanczos method, with a randomly chosen starting vector, outputs a vector v with Ã_{1}= π_{span(v)}(A) such that ∥Ã_{1}∥_{F}^{2} ≥ (1 − ϵ)∥A_{1}∥_{F}^{2} in log(n)/\sqrt{ϵ} iterations. However, this does not imply a multiplicative (1 + ϵ) error to ∥A − A_{1}∥_{F}^{2}.

While the idea of volume sampling appears to be new, in [21] it is proved that, using the rows and columns that correspond to the k × k submatrix of maximal volume, one can compute a rank-k approximation to the original matrix which differs in each entry by at most (k+1)σ_{k+1} (leading to much weaker \sqrt{mn}approximations for the Frobenius norm).

The results of our paper connect two previously separate fields: low-rank approximation and projective clustering. Algorithms and systems based on projective clustering have been applied to facial recognition, data-mining, and synthetic data [5, 32, 6], motivated by the observation that no single subspace performs as well as a few different subspaces. It should be noted that the advantage of a lowdimensional representation is not merely in the computational savings, but also the improved quality of retrieval. In [3], Agarwal and Mustafa consider the same problem as in this paper, and propose a variant of the j-means algorithm for it. Their paper has promising experimental results but does not provide any theoretical guarantees. There are theoretical results for special cases of projective clustering, especially the j-means problem (find j points). Drineas et al. [15] gave a 2-approximation to j-means using SVD. Subsequently, Ostrovsky and Rabani [31] gave the first randomized polynomial time approximation schemes for j-means (and also the j-median problem). Matoušek [29] and Effros and Schulman [17] both gave deterministic PTAS’s for j-means. Fernandez de la Vega et al. [11] describe a randomized algorithm with a running time of n(logn)^{O(1)}. Using the idea of core-sets, Har-Peled and Mazumdar [23] showed a multiplicative (1+ϵ)-approximation algorithm that runs in linear time for fixed j,ϵ. Kumar et al. [28] give a linear-time PTAS that uses random sampling. There is a PTAS for k = 1 (lines) as well [4]. Other objective functions have also been studied, e. g. sum of distances (j-median when k = 0, [31, 23]) and maximum distance (j-center when k = 0, [10]). For general k, Har-Peled and Varadarajan [24] give a multiplicative (1 + ϵ)-approximation algorithm for the maximum distance objective function. Their algorithm runs in time dn^{O(jk^{6}log(1/ϵ)/ϵ^{5})}and is based on core-sets (see [2] for a survey).

1.4 Previous versions of this paper

This paper is a journal version of the paper by the same set of authors presented at the 17th Annual Symposium on Discrete Algorithms (SODA 2006) [12]. This paper contains the same major results, but additionally provides more intuition and more complete proofs. A previous version [33] of this paper by a subset of three of the authors (L. Rademacher, S. Vempala, G. Wang) appeared as an MIT CSAIL technical report. That version contained a subset of the results in this paper. In particular, the notion of volume sampling and related results did not appear in the technical report.

1.5 Notation and preliminaries

Any m × n real matrix A has a singular value decomposition, that is, it can be written in the form

A = ∑_{i=1}^{r} σ_{i}u^{(i)}v^{(i)}^{T}

where r is the rank of A and σ_{1} ≥ σ_{2} ≥ ···≥ σ_{r} ≥ 0 are called the singular values; {u^{(1)},...,u^{(r)}} ∈ ℝ^{m}, {v^{(1)},...,v^{(r)}} ∈ ℝ^{n} are sets of orthonormal vectors, called the left and right singular vectors, respectively. It follows that A^{T}u^{(i)} = σ_{i}v^{(i)} and Av^{(i)} = σ_{i}u^{(i)} for 1 ≤ i ≤ r.

The Frobenius norm of a matrix A ∈ ℝ^{m×n} having elements (a_{ij}) is denoted ∥A∥_{F} and is given by

∥A∥_{F}^{2} = ∑_{i=1}^{m} ∑_{j=1}^{n} a_{ij}^{2} .

It satisfies ∥A∥_{F}^{2} = ∑ _{i=1}^{r} σ_{i}^{2}.

For a subspace V ⊆ ℝ^{n}, let π_{V,k}(A) denote the best rank-k approximation (under the Frobenius norm) of A with its rows in V. Let

π_{k}(A) = π_{ℝ^{n},k}(A) = ∑_{i=1}^{k} σ_{i}u^{(i)}v^{(i)}^{T}

be the best rank-k approximation of A. Also π_{V}(A) = π_{V,n}(A) is the orthogonal projection of A onto V. When we say “a set (or sample) of rows of A” we mean a set of indices of rows, rather than the actual rows. For a set S of rows of A, let span(S) ⊆ ℝ^{n} be the subspace generated by those rows; we use the simplified notation π_{S}(A) for π_{span(S)}(A) and π_{S,k}(A) for π_{span(S),k}(A).

For subspaces V,W ⊆ ℝ^{n}, their sum is denoted V + W and is given by

V + W = {x + y ∈ ℝ^{n} : x ∈ V,y ∈ W } .

The following elementary properties of the operator π_{V} will be used:

• π_{V} is linear, that is, π_{V}(λA + B) = λ π_{V}(A) + π_{V}(B) for any λ ∈ ℝ and matrices A,B ∈ ℝ^{m×n}.

• If V,W ∈ ℝ^{n} are orthogonal linear subspaces, then π_{V +W}(A) = π_{V}(A) + π_{W}(A), for any matrix A ∈ ℝ^{m×n}.

For a random vector v, its expectation, denoted E(v), is the vector having as components the expected values of the components of v.

2 Improved approximation via adaptive sampling

2.1 Proof that adaptive sampling works

We will prove Theorem 1.2 in this section. It will be convenient to formulate an intermediate theorem as follows, whose proof is quite similar to one in [19].

Theorem 2.1. Let A ∈ ℝ^{m×n}, and V ⊆ ℝ^{n} be a vector subspace. Let E = A − π_{V}(A) and let S be a random sample of s rows of A from a distribution D such that row i is chosen with probability

P_{i} = \frac{∥E^{(i)}∥^{2}}{∥E∥_{F}^{2}} .  (2.1)

Then, for any nonnegative integer k,

E_{S}(∥A − π_{V +span(S),k}(A)∥_{F}^{2}) ≤ ∥A − π_{k}(A)∥_{F}^{2} + \frac{k}{s}∥E∥_{F}^{2} .

Proof. We define vectors w^{(1)},...,w^{(k)}∈ V +span(S) such that W = span{w^{(1)},...,w^{(k)}} and show that W is a good approximation to span{v^{(1)},...,v^{(k)}} in the sense that

E_{S}(∥A − π_{W}(A)∥_{F}^{2}) ≤ ∥A − π_{k}(A)∥_{F}^{2} + \frac{k}{s}∥E∥_{F}^{2} .  (2.2)

Recall that π_{k}(A) = π_{span{v^{(1)},...,v^{(k)}}}(A), i. e., span{v^{(1)},...,v^{(k)}} is the optimal subspace upon which to project. Proving (2.2) proves the theorem, since W ⊆ V + span(S).

To this end, define X_{l}^{(j)} to be a random variable such that for i = 1,...,m and l = 1,...,s,

X_{l}^{(j)} = \frac{u_{i}^{(j)}}{P_{i}}E^{(i)} = \frac{u_{i}^{(j)}}{P_{i}}(A^{(i)}− π_{V}(A^{(i)})) with probability P_{i} .

Note that X_{l}^{(i)} is a linear function of a row of A sampled from the distribution D. Let X^{(j)} = \frac{1}{s} ∑ _{l=1}^{s} X_{l}^{(j)}, and note that E_{S}(X^{(j)}) = E^{T}u^{(j)}.

For 1 ≤ j ≤ k, define:

w^{(j)} = π_{V}(A)^{T}u^{(j)} + X^{(j)} .  (2.3)

Then we have that E_{S}(w^{(j)}) = σ_{j}v^{(j)}. We seek to bound the second central moment of w^{(j)}, i. e., E_{S}(∥w^{(j)} − σ_{j}v^{(j)}∥^{2}). We have that

w^{(j)} − σ_{j}v^{(j)} = X^{(j)} − E^{T}u^{(j)} ,

which gives us

E_{S}(∥w^{(j)} − σ_{j}v^{(j)}∥^{2}) = E_{S}(∥X^{(j)} − E^{T}u^{(j)}∥^{2})
= E_{S}(∥X^{(j)}∥^{2}) − 2E_{S}(X^{(j)}) · E^{T}u^{(j)} + ∥E^{T}u^{(j)}∥^{2}
= E_{S}(∥X^{(j)}∥^{2}) − ∥E^{T}u^{(j)}∥^{2} .  (2.4)

We evaluate the first term in (2.4),

E_{S}(∥X^{(j)}∥^{2})= E_{S}( ∥ \frac{1}{s} ∑_{l=1}^{s} X_{l}^{(j)} ∥^{2})

= \frac{1}{s^{2}} ∑_{l=1}^{s} E_{S}(∥X_{l}^{(j)}∥^{2}) + \frac{2}{s^{2}} ∑_{1≤l_{1}<l_{2}≤s} E_{S}(X_{l_{1}}^{(j)} · X_{l_{2}}^{(j)})

= \frac{1}{s^{2}} ∑_{l=1}^{s} E_{S}(∥X_{l}^{(j)}∥^{2}) + \frac{s − 1}{s} ∥E^{T}u^{(j)}∥^{2} .  (2.5)

In (2.5) we used that X_{l_{1}}^{(j)} and X_{l_{2}}^{(j)} are independent. From (2.4) and (2.5) we have that

E_{S}(∥w^{(j)} − σ_{j}v^{(j)}∥^{2}) = \frac{1}{s^{2}} ∑_{l=1}^{s} E_{S}(∥X_{l}^{(j)}∥^{2}) − \frac{1}{s}∥E^{T}u^{(j)}∥^{2} .  (2.6)

The definition of P_{i} gives us

E_{S}(∥X_{l}^{(j)}∥^{2}) = ∑_{i=1}^{m} P_{i} \frac{∥u_{i}^{(j)}E^{(i)}∥}{P_{i}^{2}} ≤ ∥E∥_{F}^{2} .  (2.7)

Thus, we have obtained a bound on the second central moment of w^{(j)}:

E_{S}(∥w^{(j)} − σ_{j}v^{(j)}∥^{2}) ≤ \frac{1}{s}∥E∥_{F}^{2} .  (2.8)

With this bound in hand, we can complete the proof. Let y^{(j)} = w^{(j)}/σ_{j} for 1 ≤ j ≤ k, and consider the matrix F = A∑ _{i=1}^{k} v^{(i)}y^{(i)^{T}}. The row space of F is contained in W = span{w^{(1)},...,w^{(k)}}. Therefore, ∥A − π_{W}(A)∥_{F}^{2} ≤ ∥A − F ∥_{F}^{2}. We will use F to bound the error ∥A − π_{W}(A)∥_{F}^{2}.

By decomposing A − F along the left singular vectors u^{(1)},...,u^{(r)}, we can use the inequality (2.8) to bound ∥A − F ∥_{F}^{2}:

E_{S}(∥A − π_{W}(A)∥_{F}^{2}) ≤ E_{S}(∥A − F ∥_{F}^{2}) = ∑_{i=1}^{r} E_{S}(∥(A − F)^{T}u^{(i)}∥_{F}^{2})

= ∑_{i=1}^{k} E_{S}(∥σ_{i}v^{(i)}− w^{(i)}∥^{2}) + ∑_{i=k+1}^{r} σ_{i}^{2}

≤ \frac{k}{s}∥E∥_{F}^{2} + ∥A − π_{k}(A)∥_{F}^{2} .  (2.9)

We can now prove Theorem 1.2 inductively using Theorem 2.1.

Proof of Theorem 1.2. We will prove the inequality by induction on t. Theorem 1.1 gives us the base case t = 1.

For the inductive step, let E = A − π_{S_{1}∪··· ∪S_{t−1}} (A). By means of Theorem 2.1 with s ≥ k/ϵ we have that

E_{S_{t}}(∥A − π_{S_{1}∪··· ∪S_{t},k}(A)∥_{F}^{2}) ≤ ∥A − π_{k}(A)∥_{F}^{2} + ϵ∥E∥_{F}^{2} .

Combining this inequality with the fact that ∥E∥_{F}^{2} ≤ ∥A − π_{S_{1}∪··· ∪S_{t−1},k}(A)∥_{F}^{2} we get

E_{S_{t}}(∥A − π_{S_{1}∪··· ∪S_{t},k}(A)∥_{F}^{2}) ≤ ∥A − π_{k}(A)∥_{F}^{2} + ϵ∥A − π_{S_{1}∪··· ∪S_{t−1},k}(A)∥_{F}^{2} .  (2.10)

Taking the expectation over S_{1},...,S_{t−1}, and using the induction hypothesis for t − 1 gives the result:

E_{S}(∥A − π_{S_{1}∪··· ∪S_{t},k}(A)∥_{F}^{2}) ≤ ∥A − π_{k}(A)∥_{F}^{2} + ϵ E_{S_{1},...,S_{t−1}}(∥A − π_{S_{1}∪··· ∪S_{t−1},k}(A)∥_{F}^{2})  (2.11)
≤ ∥A − π_{k}(A)∥_{F}^{2} + ϵ ( \frac{1}{1 − ϵ} ∥A − π_{k}(A)∥_{F}^{2} + ϵ^{t−1}∥A∥_{F}^{2})  (2.12)
= \frac{1}{1 − ϵ} ∥A − π_{k}(A)∥_{F}^{2} + ϵ^{t}∥A∥_{F}^{2} .  (2.13)

2.2 Algorithm

In this section, we present the multipass algorithm for low-rank approximation. We first describe it at a conceptual level and then give the details of the implementation.

Informally, the algorithm will find an approximation to the best rank-k subspace (the span of v^{(1)},...,v^{(k)}) by first choosing a sample T of s random rows with probabilities proportional to the squared norms of the rows (as in Theorem 1.1). Then we focus on the space orthogonal to the span of the chosen rows, that is, we consider the matrix E = A − π_{T}(A), which represents the error of our current approximation, and we sample s additional rows with probabilities proportional to the squared norms of the rows of E. We consider the union of this sample with our previous sample, and we continue adding samples in this way, up to the number of passes that we have chosen. Theorem 1.2 gives a bound on the error of this procedure.

Fast Approximate SVD

Input: A ∈ ℝ^{m×n}, integers k ≤ m, t, error parameter ϵ > 0. Output: A set of k vectors in ℝ^{n}.

1. Let S = ∅, s = k/ϵ.

2. Repeat t times:

(a) Compute the probabilities P_{i} = ∥E^{(i)}∥^{2}/∥E∥_{F}^{2} for i = 1,...,m. (One pass.)

(b) Let T be a sample of s rows of A according to the distribution that assigns probability P_{i} to row i. (Another pass.)

(c) Let S = S ∪ T.

3. Let h_{1},...,h_{k} be the top k right singular vectors of π_{S}(A).

In what follows, let M be the number of non-zeros of A.

Theorem 2.2. In 2t passes over the data, where in each pass the entries of the matrix arrive in arbitrary order, the algorithm finds vectors h_{1},...,h_{k} ∈ ℝ^{n} such that with probability at least 3/4 their span V satisfies

∥A − π_{V}(A)∥_{F}^{2} ≤ ( 1 + \frac{4ϵ}{1 − ϵ} ) ∥A − π_{k}(A)∥_{F}^{2} + 4ϵ^{t}∥A∥_{F}^{2} .  (2.14)

The running time is O( M \frac{kt}{ϵ}+ (m + n)\frac{k^{2}t^{2}}{ϵ^{2}}) .

Proof. For the correctness, observe that π_{V}(A) is a random variable with the same distribution as π_{S,k}(A) as defined in Theorem 1.2. Also, ∥A − π_{S,k}(A)∥_{F}^{2} − ∥A − π_{k}(A)∥_{F}^{2} is a nonnegative random variable and Theorem 1.2 gives a bound on its expectation:

E_{S}(∥A − π_{S,k}(A)∥_{F}^{2} − ∥A − π_{k}(A)∥_{F}^{2}) ≤ \frac{ϵ}{1 − ϵ} ∥A − π_{k}(A)∥_{F}^{2} + ϵ^{t}∥A∥_{F}^{2} .  (2.15)

Markov’s inequality applied to this variable gives that with probability at least 3/4

∥A − π_{V}(A)∥_{F}^{2} − ∥A − π_{k}(A)∥_{F}^{2} ≤ \frac{4ϵ}{1 − ϵ} ∥A − π_{k}(A)∥_{F}^{2} + 4ϵ^{t}∥A∥_{F}^{2} ,  (2.16)

which implies inequality (2.14).

We will now bound the running time. We maintain a basis of the rows indexed by S. In each iteration, we extend this basis orthogonally with a new set Y of vectors, so that it spans the new sample T. The residual squared norm of each row, ∥E^{(i)}∥^{2}, as well as the total, ∥E∥_{F}^{2}, are computed by subtracting the contribution of π_{T}(A) from the values that they had during the previous iteration. In each iteration, the projection onto Y needed for computing this contribution takes time O(Ms). In iteration i, the computation of the orthonormal basis Y takes time O(ns^{2}i) (Gram-Schmidt orthonormalization of s vectors in ℝ^{n} with reference to an orthonormal basis of size at most s(i + 1)). Thus, the total time in iteration i is O(Ms + ns^{2}i); with t iterations, this is O(Mst + ns^{2}t^{2}). At the end of Step 2 we have π_{S}(A) in terms of our basis (an m × st matrix). Finding the top k singular vectors in Step 3 takes time O(ms^{2}t^{2}). Bringing them back to the original basis takes time O(nkst). Thus, the total running time is O(Mst + ns^{2}t^{2}+ ms^{2}t^{2}+ nkst) or, in other words, O(Mkt/ϵ + (m + n)k^{2}t^{2}/ϵ^{2}).

3 Volume sampling and multiplicative approximation

We begin with a proof of Theorem 1.3, namely that volume sampling leads to a multiplicative (k + 1)- approximation (in expectation).

Proof of Theorem 1.3. For every S ⊆ {1,...,m}, let Δ_{S} be the simplex formed by {A^{(i)} : i ∈ S} and the origin. We wish to bound E_{S}(∥A − Ã_{k}∥_{F}^{2}) which can be written as follows:

E_{S}(∥A − Ã_{k}∥_{F}^{2}) = \frac{1}{∑_{S,|S|=k} vol_{k}(Δ_{S})^{2}} ∑_{S,|S|=k}vol_{k}(Δ_{S})^{2}∥A − π_{S}(A)∥_{F}^{2} .  (3.1)

For any (k +1)-subset S = {i_{1},...,i_{k+1}} of rows of A, we can express vol_{k+1}(Δ_{S})^{2} in terms of vol_{k}(Δ_{T})^{2} for T = {i_{1},...,i_{k}}, along with the squared distance from A^{(i_{k+1} )} to H_{T}, where H_{T} is the linear subspace spanned by {A^{(i)} : i ∈ T }:

vol_{k+1}(Δ_{S})^{2} = \frac{1}{(k + 1)^{2}} vol_{k}(Δ_{T})^{2}d(A^{(i_{k+1} )},H_{T})^{2} .  (3.2)

Summing over all subsets S of size k + 1:

∑_{S,|S|=k+1} vol_{k+1}(Δ_{S})^{2} = \frac{1}{k + 1} ∑_{T,|T |=k}∑_{j=1}^{m} \frac{1}{(k + 1)^{2}} vol_{k}(Δ_{T})^{2} d(A^{(j)},H_{T})^{2}

= \frac{1}{(k + 1)^{3}} ∑_{T,|T |=k}vol_{k}(Δ_{T})^{2} ∑_{j=1}^{m} d(A^{(j)},H_{T})^{2}

= \frac{1}{(k + 1)^{3}} ∑_{T,|T |=k}vol_{k}(Δ_{T})^{2}∥A − π_{T}(A)∥_{F}^{2} ,  (3.3)

where in the last step we noted that ∑ _{j=1}^{m} d(A^{(j)},H_{T})^{2} = ∥A − π_{T}(A)∥_{F}^{2}. By substituting this equality in (3.1) and applying Lemma 3.1 (proved next) twice, we have:

E_{S}(∥A − Ã_{k}∥_{F}^{2}) = \frac{1}{∑_{T,|T |=k}vol_{k}(Δ_{T})^{2}} ( (k + 1)^{3} ∑_{S,|S|=k+1} vol_{k+1}(Δ_{S})^{2})

= \frac{1}{∑_{T,|T |=k}vol_{k}(Δ_{T})^{2}} ( \frac{(k + 1)}{(k!)^{2}} ∑_{1≤t_{1}<··· <t_{k+1} ≤n} σ_{t_{1}}^{2}··· σ_{t_{k+1}}^{2})

≤ \frac{1}{∑_{T,|T |=k}vol_{k}(Δ_{T})^{2}} ( \frac{(k + 1)}{(k!)^{2}} ∑_{1≤t_{1}<··· <t_{k}≤r} σ_{t_{1}}^{2}··· σ_{t_{k}}^{2} ∑_{j=k+1}^{r} σ_{j}^{2})  (3.4)

= \frac{(k + 1)}{∑_{T,|T |=k}vol_{k}(Δ_{T})^{2}} ( ∑_{T:|T |=k}vol_{k}(Δ_{T})^{2} ∑_{j=k+1}^{r} σ_{j}^{2})

= (k + 1)∥A − A_{k}∥_{F}^{2} .

Lemma 3.1.

∑_{S,|S|=k}vol_{k}(Δ_{S})^{2} = \frac{1}{(k!)^{2}} ∑_{1≤t_{1}<t_{2}<··· <t_{k}≤r} σ_{t_{1}}^{2}σ_{t_{2}}^{2}··· σ_{t_{k}}^{2}

where σ_{1},σ_{2},···,σ_{r} > 0 = σ_{r+1} = ···= σ_{n} are the singular values of A.

Proof. Let A_{S} be the sub-matrix of A formed by the rows {A^{(i)} : i ∈ S}. Then we know that the volume of the k-simplex formed by these rows is given by vol_{k}(Δ_{S}) = \frac{1}{k!} \sqrt{det(A_{S}A_{S}^{T})}. Therefore,

∑_{S,|S|=k}vol_{k}(Δ_{S})^{2} = \frac{1}{(k!)^{2}} ∑_{S,|S|=k}det(A_{S}A_{S}^{T}) = \frac{1}{(k!)^{2}} ∑_{B: principal k-minor of AA^{T}}det(B) .  (3.5)

Let det(AA^{T}− λI) = (−1)^{m}λ^{m}+ c_{m−1}λ^{m−1} + ··· + c_{0} be the characteristic polynomial of AA^{T}. From basic linear algebra we know that the roots of this polynomial are precisely the eigenvalues of AA^{T}, i. e., σ_{1}^{2},σ_{2}^{2},...,σ_{r}^{2} and 0 with multiplicity (m − r). Moreover the coefficient c_{m−k} can be expressed in terms of these roots as:

c_{m−k} = (−1)^{m−k} ∑_{1≤t_{1}<t_{2}<··· <t_{k}≤r} σ_{t_{1}}^{2}σ_{t_{2}}^{2}··· σ_{t_{k}}^{2} .  (3.6)

But we also know that c_{m−k} is the coefficient of λ^{m−k} in det(AA^{T}− λI), which means

c_{m−k} = (−1)^{m−k} ∑_{B: principal k-minor of AA^{T}}det(B)  (3.7)

(see e. g., [8]; we prove it next as Proposition 3.2). This gives us our desired result.

Proposition 3.2. Let the characteristic polynomial of M ∈ ℝ^{m×m} be

det(M − λI_{m}) = λ^{m}+ c_{m−1}λ^{m−1} + ··· + c_{0} .

Then

c_{m−k} = (−1)^{m−k} ∑_{B: principal k-minor of AA^{T}}det(B)  for 1 ≤ k ≤ m .

Proof. First, it is clear that c_{0} = det(M). Next, let M′= M − λI, and S_{m} be the set of permutations of {1,2,...,m}. The sign of a permutation sgn(τ), for τ ∈ Perm([m]), is equal to 1 if it is a product of an even number of transpositions and −1 otherwise. For a subset S of rows, we denote the submatrix of entries (M_{i,j})_{i,j∈S} by M_{S}.

det(M − λI_{m}) = det(M′) = ∑_{τ∈Perm([m])}sgn(τ)M′_{1,τ(1)}M′_{2,τ(2)}··· M′_{m,τ(m)} .  (3.8)

The term c_{m−k}λ^{m−k} is the sum over τ which fix some set S ⊆ [m] of size (m − k), and the elements ∏ _{i∈S} M′_{i,i} contribute (−1)^{m−k}λ^{m−k} and the coefficient comes from the constant term in

∑_{τ∈Perm([m]−S)}sgn(τ)∏_{i∉S}M′_{i,τ(i)} .

Each term in this sum is the c_{0} term of a principal minor of M and so the sum is equal to

∑_{S,|S|=m−k} det(M_{[m]−S}) .

Hence

c_{m−k} = (−1)^{m−k} ∑_{S,|S|=m−k} det(M_{[m]−S}) = (−1)^{m−k} ∑_{B: principal k-minor of AA^{T}}det(B) .  (3.9)

The bound proved in Theorem 1.3 is in fact asymptotically tight:

Proposition 3.3. Given any ϵ > 0, there exists a (k + 1) × (k + 1) matrix A such that for any k-subset S of rows of A,

∥A − π_{S,k}(A)∥_{F}^{2} ≥ (1 − ϵ) (k + 1) ∥A − A_{k}∥_{F}^{2} .

Proof. The tight example consists of a matrix with k + 1 rows which are the vertices of a regular kdimensional simplex lying on the affine hyperplane {X_{k+1} = α} in ℝ^{k+1}. Let A^{(1)}, A^{(2)}, …, A^{(k+1)} be the vertices with the point p = (0,0,...,0,α) as their centroid. For α small enough, the best k dimensional subspace for these points is given by {X_{k+1} = 0} and

∥A − A_{k}∥_{F}^{2} = (k + 1)α^{2} .  (3.10)

Consider any k-subset of rows from these, say S = {A^{(1)},A^{(2)},...,A^{(k)}}, and let H_{S} be the linear subspace spanning them. Then,

∥A − π_{S,k}(A)∥_{F}^{2} = d(A^{(k+1)},H_{S})^{2} .  (3.11)

We claim that for any ϵ > 0, α can be chosen small enough so that

d(A^{(k+1)},H_{S}) ≥ \sqrt{(1 − ϵ)} (k + 1)α .  (3.12)

Choose α small enough so that d(p,H_{S}) ≥ \sqrt{(1 − ϵ)} α. Now

\frac{d(A^{(k+1)},H_{S})}{d(p,H_{S})}= \frac{d(A^{(k+1)},conv(A^{(1)},...,A^{(k)}))}{d(p,conv(A^{(1)},...,A^{(k)}))}= k + 1  (3.13)

since the points form a simplex and p is their centroid. The claim follows. Hence,

∥A − π_{S,k}(A)∥_{F}^{2} = d(A^{(k+1)},H_{S})^{2} ≥ (1 − ϵ) (k + 1)^{2} α^{2} = (1 − ϵ) (k + 1) ∥A − A_{k}∥_{F}^{2} .  (3.14)

Next, we prove Theorem 1.4. This theorem follows by interpreting Theorem 1.3 as an existence theorem and applying Theorem 2.1.

Proof of Theorem 1.4. By Theorem 1.3, there exists a k-subset S_{1} of rows of A such that

∥A − π_{S_{1}}(A)∥_{F}^{2} ≤ (k + 1)∥A − A_{k}∥_{F}^{2} .  (3.15)

Now, applying Theorem 2.1 with V = span(S_{1}) and s = k(k + 1)/ϵ we get that a random sample S_{2} of the rows of A (according to the specified distribution) satisfies

E_{S_{2}}(∥A − π_{V +span(S_{2}),k}(A)∥_{F}^{2}) ≤ (1 + ϵ)∥A − A_{k}∥_{F}^{2}  (3.16)

so there exists a subset of the rows achieving the expectation. Since V + span(S_{2}) = span(S_{1}∪ S_{2}), and |S_{1}∪ S_{2}| = k + k(k + 1)/ϵ, we have the desired result.

4 Application: Projective clustering

In this section, we give a polynomial-time approximation scheme for the projective clustering problem described in Section 1.2. We note that a simple multiplicative (k + 1)-approximation follows from Theorem 1.3 with running time O(dn^{jk}). Let V_{1},...,V_{j} be the optimal subspaces partitioning the point set into P_{1} ∪ ···∪ P_{j}, where P_{i} is the subset of points closest to V_{i}. Theorem 1.3 tells us that each P_{i} contains a subset S_{i} of k points whose span V′_{i}is a (k + 1)-approximation, i. e.,

∑_{p∈P_{i}}d(p,V′_{i})^{2} ≤ (k + 1) ∑_{p∈P_{i}}d(p,V_{i})^{2} .

We can find the S_{i}’s by simply enumerating all possible subsets of k points, considering j of them at a time, and taking the best of these. This leads to the complexity of dn^{jk}.

Getting a PTAS will be a bit more complicated. Theorem 1.4 implies that there exists a set P̂_{i}⊆ P_{i} of size k + k(k + 1)/ϵ in whose span lies an approximately optimal k-dimensional subspace W_{i}. We can enumerate over all combinations of j subsets, each of size k + k(k + 1)/ϵ to find the P̂_{i}, but we cannot enumerate the infinitely many k-dimensional subspaces lying in the span of P̂_{i}. One natural approach to solve this problem would be to put a finite grid down in a unit ball in the span of P̂_{i}. The hope would be that there are k grid points whose span G is “close” to W_{i}, since each basis vector for W_{i} is close to a grid point. However, this will not work; consider a point p very far from the origin. Although the distance between a basis vector and a grid point might be small, the error induced by projecting p onto a grid point is proportional to its distance to the origin, which could be too large.

The problem described above suggests that a grid construction must be dependent on the point set P_{i}. Our grid construction considers grid points in the span of P̂_{i}, but instead of a uniform grid in a unit ball, we consider grid points at bounded distance from each p ∈ π_{span(P̂_{i})}(P_{i}), i. e., the points in P_{i} projected to the span of P̂_{i}. This avoids the problem of points far from the origin, since there are grid points around each point. Note that we only put grid points around projected points. This is because we seek a subspace “close” to W_{i}, which itself lies in the span of P̂_{i}; W_{i} and any subspace lying in the span of P̂_{i} incur the same error for the component of a point orthogonal to the span of P̂_{i}. In Lemma 4.1, we show that there exists a subspace spanned by k points in our grid that is not much worse than W_{i}. The lemma is stated for a general point set, but we apply it to the projected points in Theorem 1.5.

The algorithm is given below.

Algorithm Cluster

Input: P ⊆ ℝ^{d}, error parameter 0 < ϵ < 1, and upper bound B on the optimal cost. Output: A set {F_{1},...,F_{j}} of j k-dimensional subspaces.

1. Set δ = \frac{ϵ\sqrt{B}}{16jk\sqrt{(1+ \frac{ϵ}{2})n}} , R = \sqrt{(1 + \frac{ϵ}{2})B} + 2δk.

2. For each subset T of P of size j(k + 2k(k + 1)/ϵ):

(a) For each equipartition T = T_{1}∪ ··· ∪ T_{j}:

i. For each i, construct a δ-net D_{i} with radius R for the projection of P to the span

of T_{i}. ii. For each way of choosing j subspaces F_{1},...,F_{j}, where F_{i} is the span of k points

from D_{i}, compute the cost 𝒞(F_{1},...,F_{j}).

3. Report the subspaces F_{1},...,F_{j} of minimum cost 𝒞(F_{1},...,F_{j}).

In Step 2(a)i, we construct a δ-net D_{i}. A δ-net D with radius R for S is a set such that for any point q for which d(q,p) ≤ R, for some p ∈ S, there exists a g ∈ D such that d(q,g) ≤ δ. The size of a δ-net is exponential in the dimension of S. This is why it is crucial that we construct the δ-net for P projected to the span of T_{i}. By doing so, we reduce the dimension from d to O(k^{2}/ϵ). The correctness of the algorithm relies crucially on the next lemma.

Lemma 4.1. Let δ > 0. Let P be a point set, with |P| = n and W be a subspace of dimension k. Let D be a δ-net with radius R for P, satisfying

R ≥ \sqrt{∑_{p∈P}d(p,W)^{2}}+ 2δk .

Then there exists a subspace F spanned by k points in D such that:

∑_{p∈P}d(p,F)^{2} ≤ ∑_{p∈P}d(p,W)^{2}+ 4k^{2}nδ^{2}+ 4kδ ∑_{p∈P}d(p,W) .  (4.1)

Proof. We construct the subspace F in k steps. Let F_{0} = W. Inductively, in step i, we choose a point p_{i} and rotate F_{i−1} so that it includes a grid point g_{i} around p_{i}. The subspace resulting from the last rotation, F_{k}, is the subspace F with the bound promised by the lemma. To prove that (4.1) holds, we prove the following inequality for any point p ∈ P going from F_{i−1} to F_{i}

d(p,F_{i}) ≤ d(p,F_{i−1} ) + 2δ .  (4.2)

Summing over the k steps, squaring, and summing over n points, we have the desired result.

Let G_{1} = {0}. G_{i} will be the span of the grid points {g_{1},g_{2},...,g_{i−1} }. We describe how to construct the rotation R_{i}. Let p_{i} ∈ P maximize ∥π_{G_{i}^{⊥}}(π_{F_{i−1}}(p_{i}))∥ and let g_{i} ∈ D minimize d(π_{F_{i−1}}(p_{i}),g_{i}). The point p_{i} is chosen as the furthest point from the origin in the subspace of F_{i−1} orthogonal to G_{i}. The grid point g_{i} is the point closest to p_{i} in F_{i−1} . Consider the plane Z defined by

π_{G_{i}^{⊥}}(g_{i}), π_{G_{i}^{⊥}}(π_{F_{i−1}}(p_{i})), and 0 .

Let θ be the angle between π_{G_{i}^{⊥}}(g_{i}) and π_{G_{i}^{⊥}}(π_{F_{i−1}}(p_{i})). Let R_{i} be the rotation in the plane Z by the angle θ, and define F_{i} = R_{i}F_{i−1} . Set G_{i+1} = G_{i}+ span{g_{i}}. By choosing the rotation dependent on p_{i}, we ensure that no point moves more than p_{i}. This allows us to prove (4.2) for all points p.

Now we prove inequality (4.2). We do so by proving the following inequality by induction on i for any point p:

d(π_{F_{i−1}}(p),R_{i}π_{F_{i−1}}(p)) ≤ 2δ .  (4.3)

Note that this proves (4.2) by applying the triangle inequality, since:

d(p,F_{i})≤ d(p,F_{i−1} ) + d(π_{F_{i−1}}(p),π_{F_{i}}(p))  (4.4)
≤ d(p,F_{i−1} ) + d(π_{F_{i−1}}(p),R_{i}π_{F_{i−1}}(p)) .  (4.5)

The base case of the inequality, i = 1, is trivial. Consider the inductive case; here, we are bounding the distance between π_{F_{i−1}}(p) and R_{i}π_{F_{i−1}}(p). It suffices to bound the distance between these two points in the subspace orthogonal to G_{i}, since the rotation R_{i} is chosen orthogonal to G_{i}. That is,

d(π_{F_{i−1}}(p),R_{i}π_{F_{i−1}}(p)) ≤ d(π_{G_{i}^{⊥}}(π_{F_{i−1}}(p)),R_{i}π_{G_{i}^{⊥}}(π_{F_{i−1}}(p))) .  (4.6)

Now, consider the distance between a point π_{G_{i}^{⊥}}(π_{F_{i−1}}(p)) and its rotation, R_{i}π_{G_{i}^{⊥}}(π_{F_{i−1}}(p)). This distance is maximized when ∥π_{G_{i}^{⊥}}(π_{F_{i−1}}(p))∥ is maximized, so we have, by construction, that the maximum value is achieved by p_{i}:

d(π_{G_{i}^{⊥}}(π_{F_{i−1}}(p)),R_{i}π_{G_{i}^{⊥}}(π_{F_{i−1}}(p))) ≤ d(π_{G_{i}^{⊥}}(π_{F_{i−1}}(p_{i})),R_{i}π_{G_{i}^{⊥}}(π_{F_{i−1}}(p_{i}))) .  (4.7)

By the triangle inequality we have:

d(π_{G_{i}^{⊥}}(π_{F_{i−1}}(p_{i})),R_{i}π_{G_{i}^{⊥}}(π_{F_{i−1}}(p_{i}))) ≤ d(π_{G_{i}^{⊥}}(π_{F_{i−1}}(p_{i})),π_{G_{i}^{⊥}}(g_{i})) + d(π_{G_{i}^{⊥}}(g_{i}),R_{i}π_{G_{i}^{⊥}}(π_{F_{i−1}}(p_{i}))) .

To bound the first term, d(π_{G_{i}^{⊥}}(π_{F_{i−1}}(p_{i})),π_{G_{i}^{⊥}}(g_{i})), note that

d(π_{G_{i}^{⊥}}(π_{F_{i−1}}(p_{i})),π_{G_{i}^{⊥}}(g_{i})) ≤ d(π_{F_{i−1}}(p_{i}),g_{i}) .

We show that π_{F_{i−1}}(p_{i}) is within a ball of radius R around p_{i}; this implies

d(π_{F_{i−1}}(p_{i}),g_{i}) ≤ δ  (4.8)

by construction of the δ-net around p_{i}. We have:

d(p_{i},π_{F_{i−1}}(p_{i}))≤ d(p_{i},F_{0}) + ∑_{j=1}^{i−2} d(π_{F_{j}}(p_{i}),π_{F_{j+1}} (p_{i}))

≤ \sqrt{∑_{p∈P}d(p,W)^{2}}+ ∑_{j=1}^{i−2} d(π_{F_{j}}(p_{i}),R_{j+1} π_{F_{j}}(p_{i}))

≤ \sqrt{∑_{p∈P}d(p,W)^{2}}+ 2δk ≤ R .

The third line uses the induction hypothesis.

Now we bound the second term, d(π_{G_{i}^{⊥}}(g_{i}),R_{i}π_{G_{i}^{⊥}}(π_{F_{i−1}}(p_{i}))). Note that R_{i}π_{G_{i}^{⊥}}(π_{F_{i−1}}(p_{i})) is just a rescaling of π_{G_{i}^{⊥}}(g_{i}) and that ∥π_{G_{i}^{⊥}}(π_{F_{i−1}}(p^{∗}))∥ = ∥R_{i}π_{G_{i}^{⊥}}(π_{F_{i−1}}(p^{∗}))∥, since rotation preserves norms. The bound on the first term implies that ∥π_{G_{i}^{⊥}}(g_{i})∥ ≥ ∥π_{G_{i}^{⊥}}(π_{F_{i−1}}(p^{∗}))∥ − δ, so

d(π_{G_{i}^{⊥}}(g_{i}),R_{i}π_{G_{i}^{⊥}}(π_{F_{i−1}}(p_{i}))) ≤ δ .  (4.9)

Combining (4.8) and (4.9), we have proved (4.3).

Now, we are ready to prove Theorem 1.5, which includes the correctness of the algorithm.

Proof of Theorem 1.5. Assume that the optimal solution is of value at most B. Let V_{1},...,V_{j} be the optimal subspaces, and let P_{1},...,P_{j} be the partition of P such that P_{i} is the subset of points closest to V_{i}. Theorem 1.4 implies that there exists S_{i} ⊆ P_{i} of size at most k + 2k(k + 1)/ϵ such that there is a k-dimensional subspace W_{i} in the span of S_{i} with

∑_{p∈P_{i}}d(p,W_{i})^{2} ≤ ( 1 + \frac{ϵ}{2}) ∑_{p∈P_{i}}d(p,V_{i})^{2} ≤ ( 1 + \frac{ϵ}{2}) B .  (4.10)

Consider π_{span(S_{i})}(P_{i}), the projection of P_{i} to span(S_{i}). We want to apply Lemma 4.1 to π_{span(S_{i})}(P_{i}) and W_{i} with radius R and δ as in the algorithm. Note that the optimal solution is of value at most B, so we have that:

\sqrt{∑_{p∈π_{span(S_{i})}(P_{i})}d(p,W_{i})^{2}}+ 2δk ≤ \sqrt{∑_{p∈P_{i}}d(p,W_{i})^{2}}+ 2δk ≤ \sqrt{( 1 + \frac{ϵ}{2}) B}+ 2δk = R .  (4.11)

Theory of Computing, Volume 2 (2006), pp. 225–247 241

Let F_{i} be the subspace spanned by k points from the δ-net D_{i} for π_{span(S_{i})}(P_{i}) promised by Lemma 4.1. For every i, we have that:

∑_{p∈π_{span(S_{i})}(P_{i})}d(p,F_{i})^{2}≤ ∑_{p∈π_{span(S_{i})}(P_{i})}d(p,W_{i})^{2}+ 4k^{2}nδ^{2}+ 4kδ ∑_{p∈π_{span(S_{i})}(P_{i})}d(p,W_{i})

≤ ∑_{p∈π_{span(S_{i})}(P_{i})}d(p,W_{i})^{2}+ \frac{ϵ}{4j} B + 4kδ( n ∑_{p∈π_{span(S_{i})}(P_{i})}d(p,W_{i})^{2}) ^{1/2}

≤ ∑_{p∈π_{span(S_{i})}(P_{i})}d(p,W_{i})^{2}+ \frac{ϵ}{2j} B .  (4.12)

Now, for any point p ∈ P_{i}, we can decompose the squared distance from p to F_{i} as follows:

d(p,F_{i})^{2} = d(π_{span(S_{i})}(p),F_{i})^{2}+ d(π_{span(S_{i})^{⊥}} (p),F_{i})^{2} .

The same decomposition can be done for d(p,W_{i})^{2}. Now, since F_{i} and W_{i} both lie in the span of S_{i}, we have the following for any point p ∈ P_{i}: d(π_{span(S_{i})^{⊥}} (p),F_{i})^{2} = d(π_{span(S_{i})^{⊥}} (p),W_{i})^{2}. Applying this to (4.12) and (4.10), we have:

∑_{p∈P_{i}}d(p,F_{i})^{2} ≤ ( 1 + \frac{ϵ}{2}) ∑_{p∈P_{i}}d(p,V_{i})^{2}+ \frac{ϵ}{2j} B .

Let S = ⋃_{i}S_{i}. The algorithm will enumerate S in Step 2a, and it will enumerate the partition S_{1}∪··· ∪S_{j} in Step 2a. In Step 2(a)i, the algorithm will, for each i, construct a δ-net D_{i} for π_{span(S_{i})}(P). Lastly, in Step 2(a)ii it will consider the subspaces F_{1},...,F_{j} whose existence is proven above. The cost associated with this solution is:

𝒞(F_{1},...,F_{j}) ≤ ∑_{i=1}^{j} ∑_{p∈P_{i}}d(p,F_{i})^{2} ≤ ( 1 + \frac{ϵ}{2}) ∑_{i=1}^{j} ∑_{p∈P_{i}}d(p,V_{i})^{2}+ \frac{ϵ}{2}B = (1 + ϵ)B .

The number of subsets of size k + 2k(k + 1)/ϵ enumerated by the algorithm is at most \binom{n}{2(k+1)^{2}/ϵ}^{j}. A δ-net D with radius R for a point set Q of dimension d is implemented by putting a box with side length 2R of grid width δ/\sqrt{d} around each point in Q. Let X be the set of grid points in the box around a point p. The number of subspaces in each δ-net D_{i} is therefore at most the number of j subspaces that one can choose for a partition T_{1}∪ ··· ∪ T_{j} is (n|X|)^{jk}. The computation for projecting points, finding a basis, and determining the cost of a candidate family of subspaces takes time O(ndjk). The cardinality of X is (for ϵ < 1):

|X| = ( \frac{2R}{δ/\sqrt{2(k + 1)^{2}/ϵ}} ) ^{2(k+1)^{2}/ϵ} ≤ ( O( jk\sqrt{\frac{n}{ϵ}}) ) ^{2(k+1)^{2}/ϵ} .  (4.13)

Therefore, the running time of the algorithm is at most O(ndjk) n^{2j(k+1)^{2}/ϵ} (n|X|)^{jk} = d ( \frac{n}{ϵ})^{O(jk^{3}/ϵ)}.

Theory of Computing, Volume 2 (2006), pp. 225–247 242

5 Conclusions and subsequent work

Theorem 1.4 was further improved in [13] to show that for any real matrix, there exist O(k/ϵ + klogk) rows whose span contains the rows of a multiplicative (1 + ϵ)-approximation to the best rank-k matrix. Using this subset of O(k/ϵ + klogk) rows, the exponent in the running time of the projective clustering algorithm decreases from O(jk^{3}/ϵ) to O(jk^{2}/ϵ). It would be interesting to know if we can compute this set of O(k/ϵ + klogk) rows efficiently in a small number of passes. It would also be nice to see other applications of volume sampling.

In a recent result, Sarlos [34] proved that, using random linear combinations of rows instead of a subset of rows, we can compute a multiplicative (1 + ϵ)-approximation using only 2 passes.

References

[1] * D. Achlioptas and F. McSherry: Fast computation of low rank matrix approximations. In Proc. 33rd STOC., pp. 611–618. ACM Press, 2001. [STOC:380752.380858]. 1.1, 1.3

[2] * P. Agarwal, S. Har-Peled, and K. Varadarajan: Geometric approximations via coresets. Combinatorial and Computational Geometry - MSRI Publications, 52:1–30, 2005. 1.2, 1.3

[3] * P. Agarwal and N. Mustafa: k-means projective clustering. In Proc. Principles of Database Systems (PODS’04), pp. 155–165. ACM Press, 2004. [PODS:1055558.1055581]. 1.1, 1.3

[4] * P. Agarwal, C. Procopiuc, and K. Varadarajan: Approximation algorithms for a k-line center. Algorithmica, 42:221–230, 2005. [Algorithmica:p321725768012505]. 1.3

[5] * R. Agarwal, J. Gehrke, D. Gunopulos, and P. Raghavan: Automatic subspace clustering of high dimensional data for data mining applications. Data Mining and Knowledge Discovery, 11:5–33, 2005. [Springer:m18077t9134v55n2]. 1.1, 1.3

[6] * C. Aggarwal, C. Procopiuc, J. Wolf, P. Yu, and J. Park: Fast algorithms for projected clustering. In Proc. Conf. on Management of Data (SIGMOD’99), pp. 61–72. ACM Press, 1999. [SIGMOD:304182.304188]. 1.1, 1.3

[7] * N. Alon, Y. Matias, and M. Szegedy: The space complexity of approximating the frequency moments. J. Computer and System Sciences, 58:137–147, 1999. [JCSS:10.1006/jcss.1997.1545]. 1.3

[8] * M. Artin: Algebra. Prentice-Hall, 1991. 3

[9] * Z. Bar-Yossef: Sampling lower bounds via information theory. In Proc. 35th STOC., pp. 335–344. ACM Press, 2003. [STOC:780542.780593]. 1.3

[10] * M. Bădoiu, S. Har-Peled, and P. Indyk: Approximate clustering via core-sets. In Proc. 34th STOC., pp. 250–257. ACM Press, 2002. [STOC:509907.509947]. 1.3 [11] * W. Fernandez de la Vega, M. Karpinski, C. Kenyon, and Y. Rabani: Approxi-

mation schemes for clustering problems. In Proc. 35th STOC., pp. 50–58. ACM Press, 2003. [STOC:780542.780550]. 1.3

[12] * A. Deshpande, L. Rademacher, S. Vempala, and G. Wang: Matrix approximation and projective clustering via volume sampling. In Proc. 17th ACM-SIAM Symp. on Discrete Algorithms (SODA’06), pp. 1117–1126. SIAM, 2006. [SODA:1109557.1109681]. 1.4

[13] * A. Deshpande and S. Vempala: Adaptive sampling and fast low-rank matrix approximation. In Proc. 10th Internat. Workshop on Randomization and Computation (RANDOM’06), pp. 292– 303. Springer, 2006. 5

[14] * P. Drineas and R. Kannan: Pass efficient algorithms for approximating large matrices. In Proc. 14th ACM-SIAM Symp. on Discrete Algorithms (SODA’03), pp. 223–232. SIAM, 2003. [SODA:644108.644147]. 1.1, 1.2, 1.3

[15] * P. Drineas, R. Kannan, A. Frieze, S. Vempala, and V. Vinay: Clustering large graphs via the singular value decomposition. Machine Learning, 56:9–33, 2004. [ML:u424k6nn6k622788]. 1.1, 1.1, 1.2, 1.3, 1.3

[16] * P. Drineas, R. Kannan, and M. Mahoney: Fast Monte Carlo algorithms for matrices II: Computing a low-rank approximation to a matrix. SIAM J. on Computing, 36:158–183, 2006. [SICOMP:10.1137/S0097539704442696]. 1.1

[17] * M. Effros and L. J. Schulman: Rapid near-optimal VQ design with a deterministic data net. In Proc. Int. Symp. on Information Theory, p. 298. IEEE, 2004. [ISIT:10.1109/ISIT.2004.1365336]. 1.3

[18] * J. Feigenbaum, S. Kannan, A. McGregor, S. Suri, and J. Zhang.: On graph problems in a semi-streaming model. Theoretical Computer Science, 348:207–216, 2005. [TCS:10.1016/j.tcs.2005.09.013]. 1.3

[19] * A. Frieze, R. Kannan, and S. Vempala: Fast Monte Carlo algorithms for finding low-rank approximations. Journal of the ACM, 51:1025–1041, 2004. [JACM:1039488.1039494]. 1.1, 1.1, 1, 1.2, 1.3, 2.1

[20] * G. H. Golub and C. F. van Loan: Matrix Computations. Johns Hopkins, third edition, 1996. 1.3

[21] * S. A. Goreinov and E. E. Tyrtyshnikov: The maximal-volume concept in approximation by low-rank matrices. Contemporary Mathematics, 280:47–51, 2001. 1.3

[22] * S. Guha, N. Koudas, and K. Shim: Approximation and streaming algorithms for histogram construction problems. ACM Transactions on Database Systems, 31:396–438, 2006. [TODS:1132863.1132873]. 1.3 [23] * S. Har-Peled and S. Mazumdar: On coresets for k-means and k-median clustering. In

Proc. 36th STOC., pp. 291–300. ACM Press, 2004. [STOC:1007352.1007400]. 1.3

[24] * S. Har-Peled and K. Varadarajan: Projective clustering in high dimensions using core-sets. In Proc. 18th Symp. on Comp. Geometry (SOCG), pp. 312–318. ACM Press, 2002. [SOCG:513400.513440]. 1.2, 1.3

[25] * M. Henzinger, P. Raghavan, and S. Rajagopalan: Computing on data streams. External Memory Algorithms – DIMACS Series in Discrete Mathematics and Computer Science, 50:107– 118, 1999. 1.3

[26] * D. Kempe and F. McSherry: A decentralized algorithm for spectral analysis. In Proc. 36th STOC., pp. 561–568. ACM Press, 2004. [STOC:1007352.1007438]. 1.3

[27] * J. Kuczyński and H. Woźniakowski: Estimating the largest eigenvalues by the power and Lanczos algorithms with a random start. SIAM J. Matrix Anal. Appl., 13:1094–1122, 1992. [SIMAX:10.1137/0613066]. 1.3

[28] * A. Kumar, Y. Sabharwal, and S. Sen.: A simple linear time (1 + ϵ)-approximation algorithm for k-means clustering in any dimensions. In Proc. 45th FOCS, pp. 454–462. IEEE Computer Society Press, 2004. [FOCS:10.1109/FOCS.2004.7]. 1.3

[29] * J. Matoušek: On approximate geometric k-clustering. Discrete and Computational Geometry, 24:61–84, 2000. [Springer:xawxbau59gtjtvv6]. 1.3

[30] * N. Megiddo and A. Tamir: On the complexity of locating linear facilities in the plane. Operations Research Letters, 13:194–197, 1982. [Elsevier:10.1016/0167-6377(82)90039-6]. 1.1

[31] * R. Ostrovsky and Y. Rabani: Polynomial time approximation schemes for geometric minsum median clustering. Journal of the ACM, 49:139–156, 2002. [JACM:506147.506149]. 1.3

[32] * C. Procopiuc, P. Agarwal, T. Murali, and M. Jones: A Monte Carlo algorithm for fast projective clustering. In Proc. Conf. on Management of Data (SIGMOD’02), pp. 418–427. ACM Press, 2002. [SIGMOD:564691.564739]. 1.1, 1.3

[33] * L. Rademacher, S. Vempala, and G. Wang: Matrix approximation and projective clustering. Technical Report 2005–018, MIT CSAIL, 2005. 1.4

[34] * T. Sarlós: Improved approximation algorithms for large matrices via random projection. In Proc. 47th FOCS, pp. 143–152. IEEE Computer Society Press, 2006. 5 AUTHORS

Amit Deshpande 2-342 77 Mass. Ave., Cambridge, MA 02139, USA. amitd mit edu

http://www.mit.edu/~amitd/

Luis Rademacher 2-331 77 Mass. Ave. Cambridge, MA 02139, USA. lrademac math mit edu

http://www-math.mit.edu/~lrademac/

Santosh Vempala College of Computing Georgia Institute of Technology Atlanta, GA 30332, USA. vempala cc gatech edu

http://www-math.mit.edu/~vempala/

Grant Wang Yahoo! 701 First Avenue Sunnyvale, CA 94089, USA. gjw alum mit edu

http://theory.csail.mit.edu/~gjw/

ABOUT THE AUTHORS

Amit Deshpande is a graduate student in Applied Math at MIT. He did his undergraduate studies at Chennai Mathematical Institute (CMI) in India. Apart from math and theory, he enjoys north Indian classical music.

Luis Rademacher is a Ph. D. candidate in the Department of Mathematics at MIT, supervised by Santosh Vempala. His research interests include game theory, matrix approximation, computational lower bounds, and the intersection between geometry and algorithms. He grew up in Chile and enjoys music as a hobby.

Santosh Vempala is a professor in the College of Computing and director of the newly formed Algorithms and Randomness Center at Georgia Tech. His research interests, ironically, are in algorithms, randomness, and geometry. He graduated from CMU in 1997 following the advice of Avrim Blum and was at MIT till 2006 except for a year as a Miller fellow at UC Berkeley. He gets unreasonably excited when a phenomenon that appears complex from one perspective turns out to be simple from another.

Grant Wang graduated from MIT in August 2006 with a Ph. D. in computer science. His advisor was Santosh Vempala. He graduated from Cornell University with a B. S. in Computer Science in 2001. His research interests are in algorithms, machine learning, and data mining. As of September 2006, he is working at Yahoo!