Large-scale matrices emerging from stocks, genomes, web documents, web images and videos everyday bring new challenges in modern data analysis. Most efforts have been focused on manipulating, understanding and interpreting large-scale data matrices. In many cases, matrix factorization methods are employed for constructing parsimonious and informative representations to facilitate computation and interpretation. A principled approach is the truncated singular value decomposition (SVD) which finds the best low-rank approximation of a data matrix. Applications of SVD such as eigenfaces (Sirovich and Kirby, 1987; Turk and Pentland, 1991) and latent semantic analysis (Deerwester et al., 1990) have been illustrated to be very successful.

However, using SVD to find basis vectors and low-rank approximations has its limita- tions. As pointed out by Berry et al. (2005), it is often useful to find a low-rank matrix approximation which posses additional structures such as sparsity or nonnegativity. Since SVD or the standard QR decomposition for sparse matrices does not preserve sparsity in general, when the sparse matrix is large, computing or even storing such decompositions becomes challenging. Therefore it is useful to compute a low-rank matrix decomposition which preserves such structural properties of the original data matrix.

Another limitation of SVD is that the basis vectors resulting from SVD have little concrete meaning, which makes it very difficult for us to understand and interpret the data in question. An example of Drineas et al. (2008) and Mahoney and Drineas (2009) has well shown this viewpoint; that is, the vector [(12)height + (1/2)income], the sum of the significant uncorrelated features from a data set of people’s features, is not particularly informative. Kuruvilla et al. (2002) have also claimed: “it would be interesting to try to find basis vectors for all experiment vectors, using actual experiment vectors and not artificial bases that offer little insight.” Therefore, it is of great interest to represent a data matrix in terms of a small number of actual columns and/or actual rows of the matrix. Matrix column selection and the CUR matrix decomposition provide such techniques.

1.1 Matrix Column Selection

Column selection has been extensively studied in the theoretical computer science (TCS) and numerical linear algebra (NLA) communities. The work in TCS mainly focuses on choosing good columns by randomized algorithms with provable error bounds (Frieze et al., 2004; Deshpande et al., 2006; Drineas et al., 2008; Deshpande and Rademacher, 2010; Bout- sidis et al., 2011; Guruswami and Sinop, 2012). The focus in NLA is then on deterministic algorithms, especially the rank-revealing QR factorizations, that select columns by pivoting rules (Foster, 1986; Chan, 1987; Stewart, 1999; Bischof and Hansen, 1991; Hong and Pan, 1992; Chandrasekaran and Ipsen, 1994; Gu and Eisenstat, 1996; Berry et al., 2005). In this paper we focus on randomized algorithms for column selection.

Given a matrix , column selection algorithms aim to choose c columns of A to construct a matrix achieves the minimum. Here “” respectively represent the matrix spectral norm, the matrix Frobenius norm, and the matrix nuclear norm, and denotes the Moore-Penrose inverse of C. Since there are () possible choices of constructing C, selecting the best subset is a hard problem.

In recent years, many polynomial-time approximate algorithms have been proposed. Among them we are especially interested in those algorithms with multiplicative upper bounds; that is, there exists a polynomial function f(m, n, k, c) such that with columns selected from A the following inequality holds

with high probability (w.h.p.) or in expectation w.r.t. C. We call f the approximation factor. The bounds are strong when for an error parameter —they are known as relative-error bounds. Particularly, the bounds are called constant-factor bounds when f does not depend on m and n (Mahoney, 2011). The relative-error bounds and constant-factor bounds of the CUR matrix decomposition and the Nystr¨om approximation are similarly defined.

However, the column selection method, also known as the decomposition in some applications, has its limitations. For a large sparse matrix A, its submatrix C is sparse, but the coefficient matrix is not sparse in general. The CX decomposition suffices when is small in size. However, when m and n are near equal, computing and storing the dense matrix X in RAM becomes infeasible. In such an occasion the CUR matrix decomposition is a very useful alternative.

1.2 The CUR Matrix Decomposition

The CUR matrix decomposition problem has been widely discussed in the literature (Gor- einov et al., 1997a,b; Stewart, 1999; Tyrtyshnikov, 2000; Berry et al., 2005; Drineas and Mahoney, 2005; Mahoney et al., 2008; Bien et al., 2010), and it has been shown to be very useful in high dimensional data analysis. Particularly, a CUR decomposition algorithm seeks to find a subset of c columns of A to form a matrix , a subset of r rows to form a matrix , and an intersection matrix small. Accordingly, we use ˜A = CUR to approximate A.

Drineas et al. (2006) proposed a CUR algorithm with additive-error bound. Later on, Drineas et al. (2008) devised a randomized CUR algorithm which has relative-error bound w.h.p. if sufficiently many columns and rows are sampled. Mackey et al. (2011) established a divide-and-conquer method which solves the CUR problem in parallel. The CUR algorithms guaranteed by relative-error bounds are of great interest.

Unfortunately, the existing CUR algorithms usually require a large number of columns and rows to be chosen. For example, for an and a target rank min{m, n}, the subspace sampling algorithm (Drineas et al., 2008)—a classical CUR algorithm— requires ) columns and ) rows to achieve relative-error bound w.h.p. The subspace sampling algorithm selects columns/rows according to the statistical leverage scores, so the computational cost of this algorithm is at least equal to the cost of the truncated SVD of A, that is, O(mnk) in general. However, maintaining a large scale matrix in RAM is often impractical, not to mention performing SVD. Recently, Drineas et al. (2012) devised fast approximation to statistical leverage scores which can be used to speedup the subspace sampling algorithm heuristically—yet no theoretical results have been reported that the leverage scores approximation can give provably efficient subspace sampling algorithm.

The CUR matrix decomposition problem has a close connection with the column selection problem. Especially, most CUR algorithms such as those of Drineas and Kannan (2003); Drineas et al. (2006, 2008) work in a two-stage manner where the first stage is a standard column selection procedure. Despite their strong resemblance, CUR is a harder problem than column selection because “one can get good columns or rows separately” does not mean that one can get good columns and rows together. If the second stage is na¨ıvely solved by a column selection algorithm on , then the approximation factor will trivially beMahoney and Drineas, 2009). Thus, more sophisticated error analysis techniques for the second stage are indispensable in order to achieve relative-error bound.

1.3 The Nystr¨om Methods

The Nystr¨om approximation is closely related to CUR, and it can potentially benefit from the advances in CUR techniques. Different from CUR, the Nystr¨om methods are used for approximating symmetric positive semidefinite (SPSD) matrices. The methods approximate an SPSD matrix only using a subset of its columns, so they can alleviate computation and storage costs when the SPSD matrix in question is large in size. In fact, the Nystr¨om methods have been extensively used in the machine learning community. For example, they have been applied to Gaussian processes (Williams and Seeger, 2001), kernel SVMs (Zhang et al., 2008), spectral clustering (Fowlkes et al., 2004), kernel PCA (Talwalkar et al., 2008; Zhang et al., 2008; Zhang and Kwok, 2010), etc.

The Nystr¨om methods approximate any SPSD matrix in terms of a subset of its columns. Specifically, given an SPSD matrix A, they require sampling c (< m) columns of A to construct an . Since there exists an permutation matrix that consists of the first c columns of , we always assume that C consists of the first c columns of A without loss of generality. We partition A and C as

where are of sizes , respectively. There are three models which are defined as follows.

• . The standard Nystr¨om approximation to A is

Here is called the intersection matrix. The matrix (is the best k-rank approximation to W, is also used as an intersection matrix for constructing approximations with even lower rank. But using results in a tighter approximation than using (

• Kumar et al., 2009). It selects a collection of t samples, each sample ), containing c columns of A. Then the ensemble method combines the samples to construct an approximation in the form of

where are the weights of the samples. Typically, the ensemble Nystr¨om method seeks to find out the weights by minimizing but effective strategy is to set the weights as

• (proposed in this paper). It is defined as

This model is not strictly the Nystr¨om method because it uses a quite different in- tersection matrix ) time to compute the Moore-Penrose inverse flops to compute matrix multiplications. The matrix multiplications can be executed very efficiently in multi-processor environment, so ideally computing the intersection matrix costs time only linear in m. This model is more accurate (which will be justified in Section 4.3 and 4.4) but more costly than the conventional ones, so there is a trade-off between time and accuracy when deciding which model to use.

Here and later, we call those which use intersection matrix , including the standard Nystr¨om and the ensemble Nystr¨om.

To generate effective approximations, much work has been built on the upper error bounds of the sampling techniques for the Nystr¨om method. Most of the work, for example, Drineas and Mahoney (2005), Li et al. (2010), Kumar et al. (2009), Jin et al. (2011), and Kumar et al. (2012), studied the additive-error bound. With assumptions on matrix coherence, better additive-error bounds were obtained by Talwalkar and Rostamizadeh (2010), Jin et al. (2011), and Mackey et al. (2011). However, as stated by Mahoney (2011), additive-error bounds are less compelling than relative-error bounds. In one recent work, Gittens and Mahoney (2013) provided a relative-error bound for the first time, where the bound is in nuclear norm.

However, the error bounds of the previous Nystr¨om methods are much weaker than those of the existing CUR algorithms, especially the relative-error bounds in which we are more interested (Mahoney, 2011). Actually, as will be proved in this paper, the lower error bounds of the standard Nystr¨om method and the ensemble Nystr¨om method are even much worse than the upper bounds of some existing CUR algorithms. This motivates us to improve the Nystr¨om method by borrowing the techniques in CUR matrix decomposition.

1.4 Contributions and Outline

The main technical contribution of this work is the adaptive sampling bound in Theorem 5, which is an extension of Theorem 2.1 of Deshpande et al. (2006). Theorem 2.1 of Deshpande et al. (2006) bounds the error incurred by projection onto column or row space, while our Theorem 5 bounds the error incurred by the projection simultaneously onto column space and row space. We also show that Theorem 2.1 of Deshpande et al. (2006) can be regarded as a special case of Theorem 5.

More importantly, our adaptive sampling bound provides an approach for improving CUR and the Nystr¨om approximation: no matter which relative-error column selection algorithm is employed, Theorem 5 ensures relative-error bounds for CUR and the Nystr¨om approximation. We present the results in Corollary 7.

Based on the adaptive sampling bound in Theorem 5 and its corollary 7, we provide a concrete CUR algorithm which beats the best existing algorithm—the subspace sampling algorithm—both theoretically and empirically. The CUR algorithm is described in Algorithm 2 and analyzed in Theorem 8. In Table 1 we present a comparison between our proposed CUR algorithm and the subspace sampling algorithm. As we see, our algorithm requires much fewer columns and rows to achieve relative-error bound. Our method is more scalable for it works on only a few columns or rows of the data matrix in question;

Table 1: Comparisons between our adaptive sampling based CUR algorithm and the best existing algorithm—the subspace sampling algorithm of Drineas et al. (2008).

Table 2: Lower bounds of the standard Nystr¨om method and the ensemble Nystr¨om method. The blanks indicate the lower bounds are unknown to us. Here m denotes the column/row number of the SPSD matrix, c denotes the number of selected columns, and k denotes the target rank.

in contrast, the subspace sampling algorithm maintains the whole data matrix in RAM to implement SVD.

Another important application of the adaptive sampling bound is to yield an algorithm for the modified Nystr¨om method. The algorithm has a strong relative-error upper bound: for a target rank k, by sampling columns it achieves relative-error bound in expectation. The results are shown in Theorem 10.

Finally, we establish a collection of lower error bounds of the standard Nystr¨om and the ensemble Nystr¨om that use as the intersection matrix. We show the lower bounds in Theorem 12 and Table 3; here Table 2 briefly summarizes the lower bounds in Table 3. From the table we can see that the upper error bound of our adaptive sampling algorithm for the modified Nystr¨om method is even better than the lower bounds of the conventional Nystr¨om methods.

The remainder of the paper is organized as follows. In Section 2 we give the notation that will be used in this paper. In Section 3 we survey the previous work on the randomized column selection, CUR matrix decomposition, and Nystr¨om approximation. In Section 4 we present our theoretical results and corresponding algorithms. In Section 5 we empirically evaluate our proposed CUR and Nystr¨om algorithms. Finally, we conclude our work in Section 6. All proofs are deferred to the appendices.

First of all, we present the notation and notion that are used here and later. We let denote the identity matrix, denote the 1 vector of ones, and 0 denote a zero vector or matrix with appropriate size. For a matrix -th column, and be a submatrix consisting of its i to j-th columns (

Let . The singular value decomposition (SVD) of A can be written as

where ) correspond to the top k singular values. We denote which is the best (or closest) rank-k approximation to A. We also use to denote the i-th largest singular value. When A is SPSD, the SVD is identical to the eigenvalue decomposition, in which case we have

We define the matrix norms as follows. Let (be the Frobenius norm, be the spectral norm, and be the nuclear norm. We always use represent

Based on SVD, the statistical leverage scores of the columns of A relative to the best rank-k approximation to A is defined as

We have that . The leverage scores of the rows of A are defined according to . The leverage scores play an important role in low-rank matrix approximation. Informally speaking, the columns (or rows) with high leverage scores have greater influence in rank-k approximation than those with low leverage scores.

Additionally, let be the Moore-Penrose inverse of A (Ben-Israel and Greville, 2003). When A is nonsingular, the Moore-Penrose inverse is identical to the matrix inverse. Given matrices is the projection of A onto the column space of is the projection of A onto the row space of Y.

Finally, we discuss the computational costs of the matrix operations mentioned above. For an general matrix ), it takes ) flops to compute the full SVD and O(mnk) flops to compute the truncated SVD of rank k (< n). The computation of ) flops. It is worth mentioning that, although multiplying an trix by an matrix runs in mnp flops, it can be easily performed in parallel (Halko et al., 2011). In contrast, implementing operations like SVD and QR decomposition in parallel is much more difficult. So we denote the time complexity of such a matrix multiplication by ), which can be tremendously smaller than O(mnp) in practice.

In Section 3.1 we present an adaptive sampling algorithm and its relative-error bound established by Deshpande et al. (2006). In Section 3.2 we highlight the near-optimal column selection algorithm of Boutsidis et al. (2011) which we will use in our CUR and Nystr¨om algorithms for column/row sampling. In Section 3.3 we introduce two important CUR algorithms. In Section 3.4 we introduce the only known relative-error algorithm for the standard Nystr¨om method.

3.1 The Adaptive Sampling Algorithm

Adaptive sampling is an effective and efficient column sampling algorithm for reducing the error incurred by the first round of sampling. After one has selected a small subset of columns (denoted ), an adaptive sampling method is used to further select a proportion of columns according to the residual of the first round, that is, . The approximation error is guaranteed to be decreasing by a factor after the adaptive sampling (Deshpande et al., 2006). We show the result of Deshpande et al. (2006) in the following lemma.

Lemma 1 (The Adaptive Sampling Algorithm) (Deshpande et al., 2006) Given a matrix columns of A, and define the residual . Additionally, for

We further sample columns i.i.d. from A, in each trial of which the i-th column is chosen with probability contain the sampled columns and let . Then, for any integer k > 0, the following inequality holds:

where the expectation is taken w.r.t.

We will establish in Theorem 5 a more general and more useful error bound for this adaptive sampling algorithm. It can be shown that Lemma 1 is a special case of Theorem 5.

3.2 The Near-Optimal Column Selection Algorithm

Boutsidis et al. (2011) proposed a relative-error column selection algorithm which requires only (1)) columns get selected. Boutsidis et al. (2011) also proved the lower bound of the column selection problem which shows that no column selection algorithm can achieve relative-error bound by selecting less than columns. Thus this algorithm is near optimal. Though an optimal algorithm recently proposed by Guruswami and Sinop (2012) attains the the lower bound, this algorithm is quite inefficient in comparison with the near-optimal algorithm. So we prefer to use the near-optimal algorithm in our CUR and Nystr¨om algorithms for column/row sampling.

The near-optimal algorithm consists of three steps: the approximate SVD via random projection (Boutsidis et al., 2011; Halko et al., 2011), the dual set sparsification algorithm

(Boutsidis et al., 2011), and the adaptive sampling algorithm (Deshpande et al., 2006). We describe the near-optimal algorithm in Algorithm 1 and present the theoretical analysis in Lemma 2.

Lemma 2 (The Near-Optimal Column Selection Algorithm) Given a matrix , a target rank . Algorithm 1 selects

columns of A to form a matrix , then the following inequality holds:

where the expectation is taken w.r.t. C. Furthermore, the matrix C can be obtained in

This algorithm has the merits of low time complexity and space complexity. None of the three steps—the randomized SVD, the dual set sparsification algorithm, and the adaptive sampling—requires loading the whole of A into RAM. All of the three steps can work on only a small subset of the columns of A. Though a relative-error algorithm recently proposed by Guruswami and Sinop (2012) requires even fewer columns, it is less efficient than the near-optimal algorithm.

3.3 Previous Work in CUR Matrix Decomposition

We introduce in this section two highly effective CUR algorithms: one is deterministic and the other is randomized.

3.3.1 The Sparse Column-Row Approximation (SCRA)

Stewart (1999) proposed a deterministic CUR algorithm and called it the sparse column-row approximation (SCRA). SCRA is based on the truncated pivoted QR decomposition via a quasi Gram-Schmidt algorithm. Given a matrix , the truncated pivoted QR decomposition procedure deterministically finds a set of columns umn pivoting, whose span approximates the column space of A, and computes an upper triangular matrix that orthogonalizes those columns. SCRA runs the same procedure again on to select a set of rows and computes the corresponding upper triangular matrix the resulting truncated pivoted QR decomposition. The intersection matrix is computed by . According to our experiments, this algorithm is quite effective but very time expensive, especially when c and r are large. Moreover, this algorithm does not have data-independent error bound.

3.3.2 The Subspace Sampling CUR Algorithm

Drineas et al. (2008) proposed a two-stage randomized CUR algorithm which has a relative-error bound with high probability (w.h.p.). In the first stage the algorithm samples c columns of A to construct C, and in the second stage it samples r rows from A and C simultaneously to construct . The sampling probabilities in the two stages are proportional to the leverage scores of A and C, respectively. That is, in the first stage the sampling probabilities are proportional to the squared -norm of the rows of ; in the second stage the sampling probabilities are proportional to the squared -norm of the rows of . That is why it is called the subspace sampling algorithm. Here we show the main results of the subspace sampling algorithm in the following lemma.

and a target rank , the subspace sampling algorithm selects columns and r = rows without replacement. Then

holds with probability at least 1 contains the rows of C with scaling. The running time is dominated by the truncated SVD of A, that is, O(mnk).

3.4 Previous Work in the Nystr¨om Approximation

In a very recent work, Gittens and Mahoney (2013) established a framework for analyzing errors incurred by the standard Nystr¨om method. Especially, the authors provided the first and the only known relative-error (in nuclear norm) algorithm for the standard Nystr¨om method. The algorithm is described as follows and, its bound is shown in Lemma 4.

Like the CUR algorithm in Section 3.3.2, the Nystr¨om algorithm also samples columns by the subspace sampling of Drineas et al. (2008). Each column is selected with probability with replacement, where are leverage scores defined in (3). After column sampling, C and W are obtained by scaling the selected columns, that is,

Here is a column selection matrix that -th column of A is the j-th column selected, and is a diagonal scaling matrix satisfying

SPSD matrix A and a target rank , the subspace sampling algorithm selects

columns without replacement and constructs C and W by scaling the selected columns. Then the inequality

We now present our main results. We establish a new error bound for the adaptive sampling algorithm in Section 4.1. We apply adaptive sampling to the CUR and modified Nystr¨om problems, obtaining effective and efficient CUR and Nystr¨om algorithms in Section 4.2 and Section 4.3 respectively. In Section 4.4 we study lower bounds of the conventional Nystr¨om methods to demonstrate the advantages of our approach. Finally, in Section 4.5 we show that our expected bounds can extend to with high probability (w.h.p.) bounds.

4.1 Adaptive Sampling

The relative-error adaptive sampling algorithm is originally established in Theorem 2.1 of Deshpande et al. (2006) (see also Lemma 1 in Section 3.1). The algorithm is based on the following idea: after selecting a proportion of columns from by an arbitrary algorithm, the algorithm randomly samples additional columns according to the residual . Here we prove a new and more general error bound for the same adaptive sampling algorithm.

Theorem 5 (The Adaptive Sampling Algorithm) Given a matrix matrix consist of , and define the residual Additionally, for

We further sample rows i.i.d. from A, in each trial of which the i-th row is chosen with probability contain the sampled rows and let . Then we have

where the expectation is taken w.r.t.

Remark 6 This theorem shows a more general bound for adaptive sampling than the original one in Theorem 2.1 of Deshpande et al. (2006). The original one bounds the error incurred by projection onto the column space of C, while Theorem 5 bounds the error incurred by projection onto the column space of C and row space of R simultaneously—such situation rises in problems such as CUR and the Nystr¨om approximation. It is worth pointing out that Theorem 2.1 of Deshpande et al. (2006) is a direct corollary of this theorem when

As discussed in Section 1.2, selecting good columns or rows separately does not ensure good columns and rows together for CUR and the Nystr¨om approximation. Theorem 5 is thereby important for it guarantees the combined effect column and row selection. Guaranteed by Theorem 5, any column selection algorithm with relative-error bound can be applied to CUR and the Nystr¨om approximation. We show the result in the following corollary.

a matrix , a target rank , and a column selection algorithm achieves relative-error upper bound by selecting columns. Then we have the following results for CUR and the Nystr¨om approximation.

(1) By selecting columns of A to construct rows to construct , both using algorithm , followed by selecting additional rows using the adaptive sampling algorithm to construct , the CUR matrix decomposition achieves relative-error upper bound in expectation:

(2) Suppose symmetric matrix. By selecting columns of A to construct and selecting columns of A to construct the adaptive sampling algorithm, the modified Nystr¨om method achieves relative-error upper bound in expectation:

Based on Corollary 7, we attempt to solve CUR and the Nystr¨om by adaptive sampling algorithms. We present concrete algorithms in Section 4.2 and 4.3.

4.2 Adaptive Sampling for CUR Matrix Decomposition

Guaranteed by the novel adaptive sampling bound in Theorem 5, we combine the near-optimal column selection algorithm of Boutsidis et al. (2011) and the adaptive sampling algorithm for solving the CUR problem, giving rise to an algorithm with a much tighter theoretical bound than existing algorithms. The algorithm is described in Algorithm 2 and its analysis is given in Theorem 8. Theorem 8 follows immediately from Lemma 2 and Corollary 7.

Theorem 8 (Adaptive Sampling for CUR) Given a matrix and a positive integer , the CUR algorithm described in Algorithm 2 randomly selects columns of A to construct , and then selects A to construct . Then we have

The algorithm costs time compute matrices C, U and R.

When the algorithm is executed in a single-core processor, the time complexity of the CUR algorithm is linear in mn; when executed in multi-processor environment where matrix multiplication is performed in parallel, ideally the algorithm costs time only linear in m+n. Another advantage of this algorithm is that it avoids loading the whole data matrix A into RAM. Neither the near-optimal column selection algorithm nor the adaptive sampling algorithm requires loading the whole of A into RAM. The most space-expensive operation throughout this algorithm is computation of the Moore-Penrose inverses of C and R, which requires maintaining an matrix or an matrix in RAM. To compute the intersection matrix , the algorithm needs to visit each entry of A, but it is not RAM expensive because the multiplication can be done by computing separately. The above analysis is also valid for the Nystr¨om algorithm in Theorem 10.

Remark 9 If we replace the near-optimal column selection algorithm in Theorem 8 by the optimal algorithm of Guruswami and Sinop (2012), it suffices to select columns and rows totally. But the optimal algorithm is less efficient than the near-optimal algorithm.

4.3 Adaptive Sampling for the Nystr¨om Approximation

Theorem 5 provides an approach for bounding the approximation errors incurred by projection simultaneously onto column space and row space. Thus this approach can be applied to solve the modified Nystr¨om method. The following theorem follows directly from Lemma 2 and Corollary 7.

Given a symmetric matrix and a target rank columns sampled by Algorithm 1 and columns sampled by the adaptive sampling algorithm, that is, with totally columns being sampled, the approximation error incurred by the modified Nystr¨om method is upper bounded by

The algorithm costs time in computing C and U.

Remark 11 The error bound in Theorem 10 is the only Frobenius norm relative-error bound for the Nystr¨om approximation at present, and it is also a constant-factor bound. If

Table 3: Lower bounds of the standard Nystr¨om method and the ensemble Nystr¨om method. The blanks indicate the lower bounds are unknown to us. Here m denotes the column/row number of the SPSD matrix, c denotes the number of selected columns, and k denotes the target rank.

one uses the optimal column selection algorithm of Guruswami and Sinop (2012), which is less efficient, the error bound is further improved: only columns are required. Furthermore, the theorem requires the matrix A to be symmetric, which is milder than the SPSD requirement made in the previous work.

This is yet the strongest result for the Nystr¨om approximation problem—much stronger than the best possible algorithms for the conventional Nystr¨om method. We will illustrate this point by revealing the lower error bounds of the conventional Nystr¨om methods.

4.4 Lower Error Bounds of the Conventional Nystr¨om Methods

We now demonstrate to what an extent our modified Nystr¨om method is superior over the conventional Nystr¨om methods (namely the standard Nystr¨om defined in (1) and the ensemble Nystr¨om in (2)) by showing the lower error bounds of the conventional Nystr¨om methods. The conventional Nystr¨om methods work no better than the lower error bounds unless additional assumptions are made on the original matrix A. We show in Theorem 12 the lower error bounds of the conventional Nystr¨om methods; the results are briefly summarized previously in Table 2.

To derive lower error bounds, we construct two adversarial cases for the Nystr¨om methods. To derive the spectral norm lower bounds, we use an SPSD matrix B whose diagonal entries equal to 1 and off-diagonal entries equal to 1). For the Frobenius norm and nuclear norm bounds, we construct an block diagonal matrix A which has k diagonal blocks, each of which is in size and constructed in the same way as B. For the lower bounds on is set to be constant; for the bounds on is set to be 1. The detailed proof of Theorem 12 is deferred to Appendix C.

Assume we are given an SPSD matrix and a target rank denote the best rank-k approximation to denote either the rank-c approximation to A constructed by the standard

Nystr¨om method in (1), or the approximation constructed by the ensemble Nystr¨om method in (2) with t non-overlapping samples, each of which contains c columns of A. Then there exists an SPSD matrix such that for any sampling strategy the approximation errors of the conventional Nystr¨om methods, that is, , are lower bounded by some factors which are shown in Table 3.

Remark 13 The lower bounds in Table 3 (or Table 2) show the conventional Nystr¨om methods can be sometimes very ineffective. The spectral norm and Frobenius norm bounds even depend on m, so such bounds are not constant-factor bounds. Notice that the lower error bounds do not meet if is replaced by , so our modified Nystr¨om method is not limited by such lower bounds.

4.5 Discussions of the Expected Relative-Error Bounds

The upper error bounds established in this paper all hold in expectation. Now we show that the expected error bounds immediately extend to w.h.p. bounds using Markov’s inequality. Let the random variable denote the error ratio, where

Then we have by the preceding theorems. By applying Markov’s inequality we have that

where s is an arbitrary constant greater than 1. Repeating the sampling procedure for t times and letting correspond to the error ratio of the i-th sample, we obtain an upper bound on the failure probability:

which decays exponentially with t. Therefore, by repeating the sampling procedure multiple times and choosing the best sample, our CUR and Nystr¨om algorithms are also guaranteed with w.h.p. relative-error bounds. It follows directly from (4) that, by repeating the sampling procedure for

holds with probability at least 1

For another instance, we let (1 + 1) times, the inequality

holds with probability at least 1

In Section 5.1 we empirical evaluate our CUR algorithms in comparison with the algorithms introduced in Section 3.3. In Section 5.2 we conduct empirical comparisons between the standard Nystr¨om and our modified Nystr¨om, and comparisons among three sampling algorithms. We report the approximation error incurred by each algorithm on each data set. The error ratio is defined by

where ˜A = CUR for the CUR matrix decomposition, ˜for the standard Nystr¨om method, and ˜for the modified Nystr¨om method.

We conduct experiments on a workstation with two Intel Xeon 2.40GHz CPUs, 24GB RAM, and 64bit Windows Server 2008 system. We implement the algorithms in MATLAB R2011b, and use the MATLAB function ‘svds’ for truncated SVD. To compare the running time, all the computations are carried out in a single thread by setting ‘maxNumCompThreads(1)’ in MATLAB.

5.1 Comparison among the CUR Algorithms

In this section we empirically compare our adaptive sampling based CUR algorithm (Algorithm 2) with the subspace sampling algorithm of Drineas et al. (2008) and the deterministic sparse column-row approximation (SCRA) algorithm of Stewart (1999). For SCRA, we use the MATLAB code released by Stewart (1999). As for the subspace sampling algorithm, we compute the leverages scores exactly via the truncated SVD. Although the fast approximation to leverage scores (Drineas et al., 2012) can significantly speedup subspace sampling, we do not use it because the approximation has no theoretical guarantee when applied to subspace sampling.

Table 4: A summary of the data sets for CUR matrix decomposition.

We conduct experiments on four UCI data sets (Frank and Asuncion, 2010) which are summarized in Table 4. Each data set is represented as a data matrix, upon which we apply the CUR algorithms. According to our analysis, the target rank k should be far less than m

Figure 1: Results of the CUR algorithms on the Enron data set.

Figure 2: Results of the CUR algorithms on the Dexter data set.

Figure 3: Results of the CUR algorithms on the Farm Ads data set.

Figure 4: Results of the CUR algorithms on the Gisette data set.

and n, and the column number c and row number r should be strictly greater than k. For each data set and each algorithm, we set k = 10 or 50, and c = ak, r = ac, where a ranges in each set of experiments. We repeat each of the two randomized algorithms 10 times, and report the minimum error ratio and the total elapsed time of the 10 rounds. We depict the error ratios and the elapsed time of the three CUR matrix decomposition algorithms in Figures 1, 2, 3, and 4.

We can see from Figures 1, 2, 3, and 4 that our adaptive sampling based CUR algorithm has much lower approximation error than the subspace sampling algorithm in all cases. Our adaptive sampling based algorithm is better than the deterministic SCRA on the Farm Ads data set and the Gisette data set, worse than SCRA on the Enron data set, and comparable to SCRA on the Dexter data set. In addition, the experimental results match our theoretical analysis in Section 4 very well. The empirical results all obey the theoretical relative-error upper bound

As for the running time, the subspace sampling algorithm and our adaptive sampling based algorithm are much more efficient than SCRA, especially when c and r are large. Our adaptive sampling based algorithm is comparable to the subspace sampling algorithm when c and r are small; however, our algorithm becomes less efficient when c and r are large. This is due to the following reasons. First, the computational cost of the subspace sampling algorithm is dominated by the truncated SVD of A, which is determined by the target rank k and the size and sparsity of the data matrix. However, the cost of our algorithm grows with c and r. Thus, our algorithm becomes less efficient when c and r are large. Second, the truncated SVD operation in MATLAB, that is, the ‘svds’ function, gains from sparsity, but our algorithm does not. The four data sets are all very sparse, so the subspace sampling algorithm has advantages. Third, the truncated SVD functions are very well implemented by MATLAB (not in MATLAB language but in Fortran/C). In contrast, our algorithm is implemented in MATLAB language, which is usually less efficient than Fortran/C.

5.2 Comparison among the Nystr¨om Algorithms

In this section we empirically compare our adaptive sampling algorithm (in Theorem 10) with some other sampling algorithms including the subspace sampling of Drineas et al. (2008) and the uniform sampling, both without replacement. We also conduct comparison between the standard Nystr¨om and our modified Nystr¨om, both use the three sampling algorithms to select columns.

We test the algorithms on three data sets which are summarized in Table 5. The experiment setting follows Gittens and Mahoney (2013). For each data set we generate a radial basis function (RBF) kernel matrix A which is defined by

where are data instances and is a scale parameter. Notice that the RBF kernel is dense in general. We set 2 or 1 in our experiments. For each data set with different

Table 5: A summary of the data sets for the Nystr¨om approximation. In the second tabular

settings of , we fix a target rank k = 10, 20 or 50 and vary c in a very large range. We will discuss the choice of in the following two paragraphs. We run each algorithm for 10 times, and report the the minimum error ratio as well as the total elapsed time of the 10 repeats. The results are shown in Figures 5, 6, and 7.

Table 5 provides useful implications on choosing the target rank denotes ratio that is not captured by the best rank-k approximation to the RBF kernel, and the parameter has an influence on the ratio the RBF kernel can be well approximated by a low-rank matrix, which implies that (i) a small k suffices when is large, and (ii) k should be set large when is small. So the settings (us take the RBF kernel in the Abalone data set as an example. When approximation well captures the kernel, so k can be safely set as small as 10; when the target rank k should be set large, say larger than 50, otherwise the approximation is rough.

The standard deviation of the leverage scores reflects whether the advanced importance sampling techniques such as the subspace sampling and adaptive sampling are useful. Figures 5, 6, and 7 show that the advantage of the subspace sampling and adaptive sampling over the uniform sampling is significant whenever the standard deviation of the leverage scores is large (see Table 5), and vise versa. Actually, as reflected in Table 5, the parameter influences the homogeneity/heterogeneity of the leverage scores. Usually, when small, the leverage scores become heterogeneous, and the effect of choosing “good” columns is significant.

The experimental results also show that the subspace sampling and adaptive sampling algorithms significantly outperform the uniform sampling when c is reasonably small, say c < 10k. This indicates that the subspace sampling and adaptive sampling algorithms are good at choosing “good” columns as basis vectors. The effect is especially evident on the

Figure 5: Results of the Nystr¨om algorithms on the RBF kernel in the Abalone data set.

Figure 6: Results of the Nystr¨om algorithms on the RBF kernel in the Wine Quality data set.

Figure 7: Results of the Nystr¨om algorithms on the RBF kernel in the Letters data set.

RBF kernel with the scale parameter 2, where the leverage scores are heterogeneous. In most cases our adaptive sampling algorithm achieves the lowest approximation error among the three algorithms. The error ratios of our adaptive sampling for the modified Nystr¨om are in accordance with the theoretical bound in Theorem 10; that is,

As for the running time, our adaptive sampling algorithm is more efficient than the subspace sampling algorithm. This is partly because the RBF kernel matrix is dense, and hence the subspace sampling algorithm costs ) time to compute the truncated SVD.

Furthermore, the experimental results show that using as the intersection matrix (denoted by “modified” in the figures) always leads to much lower error than using (denoted by “standard”). However, our modified Nystr¨om method costs more time to compute the intersection matrix than the standard Nystr¨om method costs. Recall that the standard Nystr¨om costs ) time to compute and that the mod-ified Nystr¨om costs ) time to compute . So the users should make a trade-off between time and accuracy and decide whether it is worthwhile to sacrifice extra computational overhead for the improvement in accuracy by using the modified Nystr¨om method.

In this paper we have built a novel and more general relative-error bound for the adaptive sampling algorithm. Accordingly, we have devised novel CUR matrix decomposition and Nystr¨om approximation algorithms which demonstrate significant improvement over the classical counterparts. Our relative-error CUR algorithm requires only columns and ) rows selected from the original matrix. To achieve relative-error bound, the best previous algorithm—the subspace sampling algorithm—requires c = ) columns and ) rows. Our modified Nystr¨om method is differ-ent from the conventional Nystr¨om methods in that it uses a different intersection matrix. We have shown that our adaptive sampling algorithm for the modified Nystr¨om achieves relative-error upper bound by sampling only (1)) columns, which even beats the lower error bounds of the standard Nystr¨om and the ensemble Nystr¨om. Our proposed CUR and Nystr¨om algorithms are scalable because they need only to maintain a small fraction of columns or rows in RAM, and their time complexities are low provided that matrix multiplication can be highly efficiently executed. Finally, the empirical comparison has also demonstrated the effectiveness and efficiency of our algorithms.

This work has been supported in part by the Natural Science Foundations of China (No. 61070239) and the Scholarship Award for Excellent Doctoral Student granted by Chinese Ministry of Education.

For the sake of self-contained, we attach the dual set sparsification algorithm and describe some implementation details. The deterministic dual set sparsification algorithm is established by Boutsidis et al. (2011) and severs as an important step in the near-optimal column selection algorithm (described in Lemma 2 and Algorithm 1 in this paper). We show the dual set sparsification algorithm algorithm in Algorithm 3 and its bounds in Lemma 14, and we also analyze the time complexity using our defined notation.

contain the columns of an arbitrary matrix (k < n) be a decompositions of the identity, that is, . Given an integer r with k < r < n, Algorithm 3 deterministically computes a set of weights most r of which are non-zero, such that

The weights can be computed deterministically in

Here we mention some implementation issues of Algorithm 3 which were not described in detail by Boutsidis et al. (2011). In each iteration the algorithm performs once eigenvalue decomposition: is guaranteed to be SPSD in each iteration. Since

(can be efficiently computed based on the eigenvalue decomposition of . With the eigenvalues at hand, ) can also be computed directly.

The algorithm runs in r iterations. In each iteration, the eigenvalue decomposition of requires comparisons in Line 6 each requires ). Moreover, computing ). Overall, the running time of Algorithm 3 is at most

The near-optimal column selection algorithm described in Lemma 2 has three steps: randomized SVD via random projection which costs the dual set sparsification algorithm which costs time, and the adaptive sampling algorithm which costs time. Therefore, the near-optimal column selection algorithm costs totally

We present the proofs of Theorem 5, Corollary 7, Theorem 8, and Theorem 10 in Appendices B.1, B.2, B.3, and B.4, respectively.

B.1 The Proof of Theorem 5

Theorem 5 can be equivalently expressed in Theorem 15. In order to stick to the column space convention throughout this paper, we prove Theorem 15 instead of Theorem 5.

Theorem 15 (The Adaptive Sampling Algorithm) Given a matrix matrix consist of columns of A, and define the residual

where -th column of the matrix B. Sample further columns from trials, where in each trial the i-th column is chosen with probability contain the sampled columns and contain the columns of both , all of which are columns of A. Then the following inequality holds:

where the expectation is taken w.r.t.

Proof With a little abuse of symbols, we use bold uppercase letters to denote random matrices and bold lowercase to denote random vectors, without distinguishing between random matrices/vectors and non-random matrices/vectors.

We denote the j-th column of , and the (i, j)-th entry of . Define random vectors such that for

Notice that is a linear function of a column of A sampled from the above defined distribution. We have that

Then we let

According to the construction of , we define the columns of . Note that all the random vectors lie in the subspace span(We define random vectors

where the second equality follows from Lemma 16; that is, the top right singular vectors of . Then we have that any set of random vectors lies in span(] be a random matrix, we have that span(). The expectation of

The expectation of

To complete the proof, we denote

where -th largest singular value of is the corresponding left singular vector of . The column space of F is contained in span()), and thus

We use F to bound the error

where (6) is due to that ) is orthogonal to (F both lie on the space spanned by the right singular vectors of decompose , obtaining that

where (7) follows from Lemma 16 and (8) follows from (5).

contain the top right singular vectors of R. Then the projection of A onto the row space of . Let the thin SVD of . Then the compact SVD of

According to the definition, -th column of (the column space of is orthogonal to . Finally, since , we have that is orthogonal to , that is, (which directly proves the lemma.

B.2 The Proof of Corollary 7

Since C is constructed by columns of A and the column space of C is contained in the column space of A, we have rank(. Consequently, the assumptions of Theorem 5 are satisfied. The assumptions in turn imply

and . It then follows from Theorem 5 that

which yields the error bound for CUR matrix decomposition.

When the matrix A is symmetric, the matrix consists of the rows A, and thus we can use Theorem 15 (which is identical to Theorem 5) to prove the error bound for the Nystr¨om approximation. By replacing R in Theorem 15 by , we have that

given by Lemma 17, we have that

Hence

, the following inequality holds:

denote the projection of A onto the column space of C, and ¯denote the projector onto the space orthogonal to the column space of C. It has been shown by Halko et al. (2011) that, for any matrix then the following inequalities hold:

Accordingly, is the projection of A onto the row space of further have that

and

where the last equalities follow from Since span(, which proves the lemma.

B.3 The Proof of Theorem 8

The error bound follows directly from Lemma 2 and Corollary 7. The near-optimal column selection algorithm costs time to construct C and time to construct . Then the adaptive sampling algorithm costs time to construct . Computing the Moore-Penrose inverses of time. The multiplication of ) time. So the total time complexity is

B.4 The Proof of Theorem 10

The error bound follows immediately from Lemma 2 and Corollary 7. The near-optimal column selection algorithm costs select ) columns of . Then the adaptive sampling algorithm costs time to select ) columns construct Finally it costs time to construct the intersection matrix So the total time complexity is

In Appendix C.1 we construct two adversarial cases which will be used throughout this appendix. In Appendix C.2 we prove the lower bounds of the standard Nystr¨om method. In Appendix C.3 we prove the lower bounds of the ensemble Nystr¨om method. Theorems 20, 21, 22, 24, and 25 are used for proving Theorem 12.

C.1 Construction of the Adversarial Cases

We now consider the construction of adversarial cases for the spectral norm bounds and the Frobenius norm and nuclear norm bounds, respectively.

We construct an positive definite matrix B as follows:

where 1). It is easy to verify 0 for any nonzero . We show some properties of B in Lemma 18.

be the best rank-k approximation to the matrix B defined in (9). Then we have that

where 1

Proof The squared Frobenius norm of B is

Then we study the singular values of B. Since B is SPSD, here we do not distinguish between its singular values and eigenvalues. The spectral norm, that is, the largest singular value, of B is

where the maximum is attained when is the top singular vector of B. Then the projection of B onto the subspace orthogonal to

Then for all j > 1, the j-th top eigenvalue and eigenvector , that is, the singular value and singular vector, of B satisfy

where the last equality follows from

for all 1 . Finally we have that

which complete our proofs.

C.1.2 The Adversarial Case for The Frobenius Norm and Nuclear Norm Bounds

Then we construct another adversarial case for proving the Frobenius norm and nuclear norm bounds. Let matrix with diagonal entries equal to one and off-diagonal entries equal to and we construct an block diagonal matrix A as follows:

be the best rank-k approximation to the matrix A defined in (10). Then we have that

Lemma 19 can be easily proved using Lemma 18.

C.2 Lower Bounds of the Standard Nystr¨om Method

with diagonal entries equal to one and off-diagonal entries equal to , the approximation error incurred by the standard Nystr¨om method is lower bounded by

Furthermore, the matrix (

Proof The matrix B is partitioned as in (9). The residual of the Nystr¨om approximation is

where is nonsingular when . We apply the Sherman-Morrison-Woodbury formula

to compute

According to the construction, matrix with all entries equal to follows that ) matrix with all entries equal to

Then we obtain that

which proves the Frobenius norm of the residual.

Now we compute the spectral norm of the residual. Based on the results above we have that

Similar to the proof of Lemma 18, it is easily obtained that is the top singular vector of the SPSD matrix (1 , so the top singular value is

which proves the spectral norm bound because It is also easy to show the rest singular values obey

Thus we have, for

The nuclear norm of the residual

The theorem follows from equalities (14), (15), and (16).

Now we use the matrix A constructed in (10) to show the Frobenius norm and nuclear norm lower bound. The bound is stronger than the one in Theorem 20 by a factor of k.

SPSD matrix A defined in (10), the approximation error incurred by the standard Nystr¨om method is lower bounded by

where k < m is an arbitrary positive integer.

Proof Let C consist of c column sampled from consist of columns sampled from the i-th block diagonal matrix in A. Without loss of generality, we assume ˆof the first columns of B, and accordingly ˆconsists of the top left Thus

Then it follows from Theorem 20 that

where ˆminimized when ˆ. Finally we have that

by which the Frobenius norm bound follows.

Since the matrices are all SPSD by Theorem 20, so the matrix (is also SPSD. We have that

where the former inequality follows from Theorem 20, and the latter inequality follows by minimizing w.r.t. subjecting to

Theorem 22 There exists an SPSD matrix A such that the approximation error incurred by the standard Nystr¨om method is lower bounded by

where k < m is an arbitrary positive integer.

Proof For the spectral norm bound we use the matrix A constructed in (9) and set then it follows directly from Lemma 18 and Theorem 20. For the Frobenius norm and nuclear norm bounds, we use the matrix A constructed in (10) and set 1, then it follows directly from Lemma 19 and Theorem 21.

C.3 Lower Bounds of the Ensemble Nystr¨om Method

The ensemble Nystr¨om method (Kumar et al., 2009) is previously defined in (2). To derive lower bounds of the ensemble Nystr¨om method, we assume that the t samples are non-overlapping. According to the construction of the matrix B in (9), each of the t non-overlapping samples are equally “important”, so without loss of generality we set the t samples with equal weights:

Lemma 23 Assume that the ensemble Nystr¨om method selects a collection of t samples, each sample columns of B without overlapping. For an matrix B with all diagonal entries equal to one and off-diagonal entries equal to the approximation error incurred by the ensemble Nystr¨om method is lower bounded by

where ˜. Furthermore, the matrix (

Proof We use the matrix B constructed in (9). It is easy to check that so we use the notation W instead. We assume that the samples contain the firs tc columns of B and each sample contains neighboring columns, that is,

If a sample C contains the first c columns of B, then

where is a permutation matrix. As was shown in Equation (12), () matrix with all entries equal to

Based on the properties of the matrix , we study the values of the entries of . We can express it as

and then a discreet examination reveals that can be partitioned into four kinds of regions as illustrated in Figure 8. We annotate the regions in the figure and summarize the values of entries in each region in the table below. (Region 1 and 4 are further partitioned into diagonal entries and off-diagonal entries.)

Figure 8: An illustration of the matrix for the ensemble Nystr¨om method where B is defined in (9). Here we set ensemble Nystr¨om method without overlapping, the matrix be partitioned into four regions as annotated.

Now we do summation over the entries of to compute its squared Frobenius norm:

Furthermore, since the matrices are all SPSD by Theorem 20, so their sum is also SPSD. Then the SPSD property of () follows from (18). Therefore, the nuclear norm of () equals to the matrix trace, that is,

which proves the nuclear norm bound in the lemma.

Theorem 24 Assume that the ensemble Nystr¨om method selects a collection of t samples, each sample columns of A without overlapping. For a the matrix A defined in (10), the approximation error incurred by the ensemble Nystr¨om method is lower bounded by

Proof According to the construction of is also block diagonal. We denote it by . Akin to (17), we have

Thus the approximation error of the ensemble Nystr¨om method is

where the inequality follows from Lemma 23, and the last equality follows from and kp = m. The summation in the last equality equals to

Here the inequality holds because the function is minimized when Finally we have that

which proves the nuclear norm bound in the theorem.

Theorem 25 Assume that the ensemble Nystr¨om method selects a collection of t samples, each sample columns of A without overlapping. Then there exists an SPSD matrix A such that the relative-error ratio of the ensemble Nystr¨om method is lower bounded by

Proof The theorem follows directly from Theorem 24 and Lemma 19 by setting

A. Ben-Israel and T. N. E. Greville. Generalized Inverses: Theory and Applications. Second Edition. Springer, 2003. 7

M. W. Berry, S. A. Pulatova, and G. W. Stewart. Algorithm 844: computing sparse reduced- rank approximations to sparse matrices. ACM Transactions on Mathematical Software, 31(2):252–269, 2005. 2, 3

J. Bien, Y. Xu, and M. W. Mahoney. CUR from a sparse optimization viewpoint. In Advances in Neural Information Processing Systems (NIPS). 2010. 3

C. H. Bischof and P. C. Hansen. Structure-preserving and rank-revealing QR-factorizations. SIAM Journal on Scientific and Statistical Computing, 12(6):1332–1350, 1991. 2

C. Boutsidis, P. Drineas, and M. Magdon-Ismail. Near-optimal column-based matrix recon- struction. CoRR, abs/1103.0995, 2011. 2, 8, 9, 12, 25

T. F. Chan. Rank revealing QR factorizations. Linear Algebra and Its Applications, 88: 67–82, 1987. 2

S. Chandrasekaran and I. C. F. Ipsen. On rank-revealing factorisations. SIAM Journal on Matrix Analysis and Applications, 15(2):592–622, 1994. 2

P. Cortez, A. Cerdeira, F. Almeida, T. Matos, and J. Reis. Modeling wine preferences by data mining from physicochemical properties. Decision Support Systems, 47(4):547–553, 2009. 20

S. Deerwester, S. T. Dumais, G. W. Furnas, T. K. Landauer, and R. Harshman. Indexing by latent semantic analysis. Journal of The American Society for Information Science, 41(6):391–407, 1990. 1

A. Deshpande and L. Rademacher. Efficient volume sampling for row/column subset selec- tion. In Proceedings of the 51st IEEE Annual Symposium on Foundations of Computer Science (FOCS), pages 329–338, 2010. 2

A. Deshpande, L. Rademacher, S. Vempala, and G. Wang. Matrix approximation and projective clustering via volume sampling. Theory of Computing, 2(2006):225–247, 2006. 2, 5, 8, 9, 11

P. Drineas and R. Kannan. Pass-efficient algorithms for approximating large matrices. In Proceeding of the 14th Annual ACM-SIAM Symposium on Dicrete Algorithms (SODA), pages 223–232, 2003. 3

P. Drineas and M. W. Mahoney. On the Nystr¨om method for approximating a gram matrix for improved kernel-based learning. Journal of Machine Learning Research, 6:2153–2175, 2005. 3, 5

P. Drineas, R. Kannan, and M. W. Mahoney. Fast Monte Carlo algorithms for matrices III: computing a compressed approximate matrix decomposition. SIAM Journal on Computing, 36(1):184–206, 2006. 3

P. Drineas, M. W. Mahoney, and S. Muthukrishnan. Relative-error CUR matrix decompo- sitions. SIAM Journal on Matrix Analysis and Applications, 30(2):844–881, September 2008. 2, 3, 6, 10, 16, 19

P. Drineas, M. Magdon-Ismail, M. W. Mahoney, and D. P. Woodruff. Fast approximation of matrix coherence and statistical leverage. In International Conference on Machine Learning (ICML), 2012. 3, 16

L. V. Foster. Rank and null space calculations using matrix decomposition without column interchanges. Linear Algebra and its Applications, 74:47–71, 1986. 2

C. Fowlkes, S. Belongie, F. Chung, and J. Malik. Spectral grouping using the Nystr¨om method. IEEE Transactions on Pattern Analysis and Machine Intelligence, 26(2):214– 225, 2004. 4

A. Frank and A. Asuncion. UCI machine learning repository, 2010. URL http://archive. ics.uci.edu/ml. 16, 20

A. Frieze, R. Kannan, and S. Vempala. Fast Monte Carlo algorithms for finding low-rank approximations. Journal of the ACM, 51(6):1025–1041, November 2004. ISSN 0004-5411. 2

A. Gittens and M. W. Mahoney. Revisiting the Nystr¨om method for improved large-scale machine learning. arXiv preprint arXiv:1303.1849, 2013. 5, 10, 19

S. A. Goreinov, E. E. Tyrtyshnikov, and N. L. Zamarashkin. A theory of pseudoskeleton approximations. Linear Algebra and Its Applications, 261:1–21, 1997a. 3

S. A. Goreinov, N. L. Zamarashkin, and E. E. Tyrtyshnikov. Pseudo-skeleton approxima- tions by matrices of maximal volume. Mathematical Notes, 62(4):619–623, 1997b. 3

M. Gu and S. C. Eisenstat. Efficient algorithms for computing a strong rank-revealing QR factorization. SIAM Journal on Scientific Computing, 17(4):848–869, 1996. 2

V. Guruswami and A. K. Sinop. Optimal column-based low-rank matrix reconstruction. In Proceedings of the 23rd Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), 2012. 2, 8, 9, 13, 14

I. Guyon, S. Gunn, A. Ben-Hur, and G. Dror. Result analysis of the NIPS 2003 feature selection challenge. Advances in Neural Information Processing Systems (NIPS), 2004. 16

N. Halko, P.-G. Martinsson, and J. A. Tropp. Finding structure with randomness: proba- bilistic algorithms for constructing approximate matrix decompositions. SIAM Review, 53(2):217–288, 2011. 7, 8, 30

Y. P. Hong and C. T. Pan. Rank-revealing QR factorizations and the singular value de- composition. Mathematics of Computation, 58(197):213–232, 1992. 2

R. Jin, T. Yang, and M. Mahdavi. Improved bound for the Nystr¨om method and its application to kernel classification. CoRR, abs/1111.2262, 2011. 5

S. Kumar, M. Mohri, and A. Talwalkar. Ensemble Nystr¨om method. In Advances in Neural Information Processing Systems (NIPS), 2009. 4, 5, 36

S. Kumar, M. Mohri, and A. Talwalkar. Sampling methods for the Nystr¨om method. Journal of Machine Learning Research, 13:981–1006, 2012. 5

F. G. Kuruvilla, P. J. Park, and S. L. Schreiber. Vector algebra in the analysis of genome- wide expression data. Genome Biology, 3:research0011–research0011.1, 2002. 2

M. Li, J. T. Kwok, and B.-L. Lu. Making large-scale Nystr¨om approximation possible. In International Conference on Machine Learning (ICML), 2010. 5

L. Mackey, A. Talwalkar, and M. I. Jordan. Divide-and-conquer matrix factorization. In Advances in Neural Information Processing Systems (NIPS). 2011. 3, 5

M. W. Mahoney. Randomized algorithms for matrices and data. Foundations and Trends in Machine Learning, 3(2):123–224, 2011. 2, 5

M. W. Mahoney and P. Drineas. CUR matrix decompositions for improved data analysis. Proceedings of the National Academy of Sciences, 106(3):697–702, 2009. 2, 3

M. W. Mahoney, M. Maggioni, and P. Drineas. Tensor-CUR decompositions for tensor- based data. SIAM Journal on Matrix Analysis and Applications, 30(3):957–987, 2008. 3

C. Mesterharm and M. J. Pazzani. Active learning using on-line algorithms. In ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD), 2011. 16

D. Michie, D. J. Spiegelhalter, and C. C. Taylor. Machine learning, neural and statistical classification. 1994. 20

L. Sirovich and M. Kirby. Low-dimensional procedure for the characterization of human faces. Journal of the Optical Society of America A, 4(3):519–524, Mar 1987. 1

G. W. Stewart. Four algorithms for the the efficient computation of truncated pivoted QR approximations to a sparse matrix. Numerische Mathematik, 83(2):313–323, 1999. 2, 3, 9, 16

A. Talwalkar and A. Rostamizadeh. Matrix coherence and the Nystr¨om method. arXiv preprint arXiv:1004.2008, 2010. 5

A. Talwalkar, S. Kumar, and H. Rowley. Large-scale manifold learning. In IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2008. 4

M. Turk and A. Pentland. Eigenfaces for recognition. Journal of Cognitive Neuroscience, 3 (1):71–86, 1991. 1

E. E. Tyrtyshnikov. Incomplete cross approximation in the mosaic-skeleton method. Computing, 64:367–380, 2000. 3

C. Williams and M. Seeger. Using the Nystr¨om method to speed up kernel machines. In Advances in Neural Information Processing Systems (NIPS), 2001. 4

K. Zhang and J. T. Kwok. Clustered Nystr¨om method for large scale manifold learning and dimension reduction. IEEE Transactions on Neural Networks, 21(10):1576–1587, 2010. 4

K. Zhang, I. W. Tsang, and J. T. Kwok. Improved Nystr¨om low-rank approximation and error analysis. In International Conference on Machine Learning (ICML), 2008. 4