Improving CUR Matrix Decomposition and the Nyström Approximation via Adaptive Sampling

Introduction

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 limitations. 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 $[(1/2)\textrm{age}-(1/\sqrt{2})\textrm{height}+(1/2)\textrm{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.

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

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 $c$ $(\geq k)$ columns selected from ${\bf A}$ the following inequality holds

with high probability (w.h.p.) or in expectation w.r.t. ${\bf C}$. We call $f$ the approximation factor. The bounds are strong when $f=1+\epsilon$ for an error parameter $\epsilon$—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öm approximation are similarly defined.

2 The CUR Matrix Decomposition

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 $m{\times}n$ matrix ${\bf A}$ and a target rank $k\ll\min\{m,n\}$, the subspace sampling algorithm (Drineas et al., 2008)—a classical CUR algorithm—requires ${\mathcal{O}}(k\epsilon^{-2}\log k)$ columns and ${\mathcal{O}}(k\epsilon^{-4}\log^{2}k)$ 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 ${\bf A}$, that is, ${\mathcal{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 ${\bf A}^{T}$, then the approximation factor will trivially be $\sqrt{2}f$It is because $\|{\bf A}-{\bf C}{\bf U}{\bf R}\|_{F}^{2}=\|{\bf A}-{\bf C}{\bf C}^{\dagger}{\bf A}+{\bf C}{\bf C}^{\dagger}{\bf A}-{\bf C}{\bf C}^{\dagger}{\bf A}{\bf R}^{\dagger}{\bf R}\|_{F}^{2}=\|({\bf I}-{\bf C}{\bf C}^{\dagger}){\bf A}\|_{F}^{2}+\|{\bf C}{\bf C}^{\dagger}({\bf A}-{\bf A}{\bf R}^{\dagger}{\bf R})\|_{F}^{2}\leq\|{\bf A}-{\bf C}{\bf C}^{\dagger}{\bf A}\|_{F}^{2}+\|{\bf A}-{\bf A}{\bf R}^{\dagger}{\bf R}\|_{F}^{2}\leq 2f^{2}\|{\bf A}-{\bf A}_{k}\|_{F}^{2}$, where the second equality follows from $({\bf I}-{\bf C}{\bf C}^{\dagger})^{T}{\bf C}{\bf C}^{\dagger}=0$. (Mahoney and Drineas, 2009). Thus, more sophisticated error analysis techniques for the second stage are indispensable in order to achieve relative-error bound.

3 The Nyström Methods

The Nyström approximation is closely related to CUR, and it can potentially benefit from the advances in CUR techniques. Different from CUR, the Nyström 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öm 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öm methods approximate any SPSD matrix in terms of a subset of its columns. Specifically, given an $m{\times}m$ SPSD matrix ${{\bf A}}$, they require sampling $c$ ($<m$) columns of ${{\bf A}}$ to construct an $m\times c$ matrix ${{\bf C}}$. Since there exists an $m{\times}m$ permutation matrix $\Pi$ such that \mbox{\boldmath\Pi\unboldmath}{{\bf C}} consists of the first $c$ columns of \mbox{\boldmath\Pi\unboldmath}{{\bf A}}\mbox{\boldmath\Pi\unboldmath}^{T}, we always assume that ${{\bf C}}$ consists of the first $c$ columns of ${\bf A}$ without loss of generality. We partition ${\bf A}$ and ${\bf C}$ as

where ${\bf W}$ and ${\bf A}_{21}$ are of sizes $c\times c$ and $(m{-}c)\times c$, respectively. There are three models which are defined as follows.

The Standard Nyström Method. The standard Nyström approximation to ${\bf A}$ is

Here ${\bf W}^{\dagger}$ is called the intersection matrix. The matrix $({\bf W}_{k})^{\dagger}$, where $k\leq c$ and ${\bf W}_{k}$ is the best $k$-rank approximation to ${\bf W}$, is also used as an intersection matrix for constructing approximations with even lower rank. But using ${\bf W}^{\dagger}$ results in a tighter approximation than using $({\bf W}_{k})^{\dagger}$ usually.

The Ensemble Nyström Method (Kumar et al., 2009). It selects a collection of $t$ samples, each sample ${{\bf C}^{(i)}}$, ($i=1,\cdots,t$), containing $c$ columns of ${\bf A}$. Then the ensemble method combines the samples to construct an approximation in the form of

The Modified Nyström Method (proposed in this paper). It is defined as

This model is not strictly the Nyström method because it uses a quite different intersection matrix ${\bf C}^{\dagger}{\bf A}({\bf C}^{\dagger})^{T}$. It costs ${\mathcal{O}}(mc^{2})$ time to compute the Moore-Penrose inverse ${\bf C}^{\dagger}$ and $m^{2}c$ 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 ${\bf W}^{\dagger}$ or $({\bf W}_{k})^{\dagger}$ the conventional Nyström methods, including the standard Nyström and the ensemble Nyström.

To generate effective approximations, much work has been built on the upper error bounds of the sampling techniques for the Nyström 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öm 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öm method and the ensemble Nyström method are even much worse than the upper bounds of some existing CUR algorithms. This motivates us to improve the Nyström method by borrowing the techniques in CUR matrix decomposition.

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öm approximation: no matter which relative-error column selection algorithm is employed, Theorem 5 ensures relative-error bounds for CUR and the Nyström 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; 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öm method. The algorithm has a strong relative-error upper bound: for a target rank $k$, by sampling \frac{2k}{\epsilon^{2}}\big{(}1+o(1)\big{)} 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öm and the ensemble Nyström that use ${\bf W}^{\dagger}$ 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öm method is even better than the lower bounds of the conventional Nyström methods.This can be valid because the lower bounds in Table 2 do not hold when the intersection matrix is not ${\bf W}^{\dagger}$.

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öm approximation. In Section 4 we present our theoretical results and corresponding algorithms. In Section 5 we empirically evaluate our proposed CUR and Nyström algorithms. Finally, we conclude our work in Section 6. All proofs are deferred to the appendices.

Notation

where ${{{\bf U}}_{{\bf A},k}}$ ($m{\times}k$), {{\mbox{\boldmath\Sigma\unboldmath}}_{{\bf A},k}} ($k{\times}k$), and ${{{\bf V}}_{{\bf A},k}}$ ($n{\times}k$) correspond to the top $k$ singular values. We denote {\bf A}_{k}={{{\bf U}}_{{\bf A},k}}{{\mbox{\boldmath\Sigma\unboldmath}}_{{\bf A},k}}{\bf V}_{{\bf A},k}^{T} which is the best (or closest) rank-$k$ approximation to ${\bf A}$. We also use $\sigma_{i}({\bf A})=\sigma_{{\bf A},i}$ to denote the $i$-th largest singular value. When ${\bf A}$ is SPSD, the SVD is identical to the eigenvalue decomposition, in which case we have ${{{\bf U}}_{{\bf A}}}={{{\bf V}}_{{\bf A}}}$.

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

Previous Work

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öm 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öm method.

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 ${\bf C}_{1}$), an adaptive sampling method is used to further select a proportion of columns according to the residual of the first round, that is, ${\bf A}-{\bf C}_{1}{\bf C}_{1}^{\dagger}{\bf A}$. 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.

where the expectation is taken w.r.t. ${\bf C}_{2}$.

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.

2 The Near-Optimal Column Selection Algorithm

Boutsidis et al. (2011) proposed a relative-error column selection algorithm which requires only $c={2k\epsilon^{-1}}(1{+}o(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 $c=k\epsilon^{-1}$ 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öm 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.

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 ${\bf A}$ into RAM. All of the three steps can work on only a small subset of the columns of ${\bf 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 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.2 The Subspace Sampling CUR Algorithm

Given an ${m\times n}$ matrix ${\bf A}$ and a target rank $k\ll\min\{m,n\}$, the subspace sampling algorithm selects $c={\mathcal{O}}(k\epsilon^{-2}\log k\log(1/\delta))$ columns and $r=$ {\mathcal{O}}\big{(}c\epsilon^{-2}\log c\log(1/\delta)\big{)} rows without replacement. Then

holds with probability at least $1-\delta$, where ${\bf W}$ contains the rows of ${\bf C}$ with scaling. The running time is dominated by the truncated SVD of ${\bf A}$, that is, ${\mathcal{O}}(mnk)$.

4 Previous Work in the Nyström Approximation

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

Given an $m\times m$ SPSD matrix ${\bf A}$ and a target rank $k\ll m$, the subspace sampling algorithm selects

columns without replacement and constructs ${\bf C}$ and ${\bf W}$ by scaling the selected columns. Then the inequality

holds with probability at least $0.6-\delta$.

Main Results

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öm problems, obtaining effective and efficient CUR and Nyström algorithms in Section 4.2 and Section 4.3 respectively. In Section 4.4 we study lower bounds of the conventional Nyström 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.

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 ${\bf A}$ to form ${\bf C}_{1}$ by an arbitrary algorithm, the algorithm randomly samples additional $c_{2}$ columns according to the residual ${\bf A}-{\bf C}_{1}{\bf C}_{1}^{\dagger}{\bf A}$. Here we prove a new and more general error bound for the same adaptive sampling algorithm.

where the expectation is taken w.r.t. ${\bf R}_{2}$.

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 ${\bf C}$, while Theorem 5 bounds the error incurred by projection onto the column space of ${\bf C}$ and row space of ${\bf R}$ simultaneously—such situation rises in problems such as CUR and the Nyström approximation. It is worth pointing out that Theorem 2.1 of Deshpande et al. (2006) is a direct corollary of this theorem when ${\bf C}={\bf A}_{k}$ (i.e., $c=n$, $\rho=k$, and ${\bf C}{\bf C}^{\dagger}{\bf A}={\bf A}_{k}$).

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öm 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öm approximation. We show the result in the following corollary.

By selecting $c\geq C(k,\epsilon)$ columns of ${\bf A}$ to construct ${\bf C}$ and $r_{1}=c$ rows to construct ${\bf R}_{1}$, both using algorithm ${{\mathcal{A}}_{\textrm{col}}}$, followed by selecting additional $r_{2}=c/\epsilon$ rows using the adaptive sampling algorithm to construct ${\bf R}_{2}$, the CUR matrix decomposition achieves relative-error upper bound in expectation:

where {\bf R}=\big{[}{\bf R}_{1}^{T},{\bf R}_{2}^{T}\big{]}^{T} and ${\bf U}={\bf C}^{\dagger}{\bf A}{\bf R}^{\dagger}$.

Suppose ${\bf A}$ is an $m\times m$ symmetric matrix. By selecting $c_{1}\geq C(k,\epsilon)$ columns of ${\bf A}$ to construct ${\bf C}_{1}$ using ${{\mathcal{A}}_{\textrm{col}}}$ and selecting $c_{2}=c_{1}/\epsilon$ columns of ${\bf A}$ to construct ${\bf C}_{2}$ using the adaptive sampling algorithm, the modified Nyström method achieves relative-error upper bound in expectation:

where {\bf C}=\big{[}{\bf C}_{1},{\bf C}_{2}\big{]} and {\bf U}={\bf C}^{\dagger}{\bf A}\big{(}{\bf C}^{\dagger}\big{)}^{T}.

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

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.

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 $m{\times}n$ data matrix ${\bf A}$ into RAM. Neither the near-optimal column selection algorithm nor the adaptive sampling algorithm requires loading the whole of ${\bf A}$ into RAM. The most space-expensive operation throughout this algorithm is computation of the Moore-Penrose inverses of ${\bf C}$ and ${\bf R}$, which requires maintaining an $m{\times}c$ matrix or an $r{\times}n$ matrix in RAM. To compute the intersection matrix ${\bf C}^{\dagger}{\bf A}{\bf R}^{\dagger}$, the algorithm needs to visit each entry of ${\bf A}$, but it is not RAM expensive because the multiplication can be done by computing ${\bf C}^{\dagger}{\bf a}_{j}$ for $j=1,\cdots,n$ separately. The above analysis is also valid for the Nyström algorithm in Theorem 10.

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 $c=k\epsilon^{-1}(1+o(1))$ columns and $r=c\epsilon^{-1}(1+\epsilon)$ rows totally. But the optimal algorithm is less efficient than the near-optimal algorithm.

3 Adaptive Sampling for the Nyström 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öm method. The following theorem follows directly from Lemma 2 and Corollary 7.

The error bound in Theorem 10 is the only Frobenius norm relative-error bound for the Nyström approximation at present, and it is also a constant-factor bound. If one uses the optimal column selection algorithm of Guruswami and Sinop (2012), which is less efficient, the error bound is further improved: only $c=\frac{k}{\epsilon^{2}}(1+o(1))$ columns are required. Furthermore, the theorem requires the matrix ${\bf 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öm approximation problem—much stronger than the best possible algorithms for the conventional Nyström method. We will illustrate this point by revealing the lower error bounds of the conventional Nyström methods.

4 Lower Error Bounds of the Conventional Nyström Methods

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

The lower bounds in Table 3 (or Table 2) show the conventional Nyström 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 ${\bf W}^{\dagger}$ is replaced by ${\bf C}^{\dagger}{\bf A}({\bf C}^{\dagger})^{T}$, so our modified Nyström method is not limited by such lower bounds.

5 Discussions of the Expected Relative-Error Bounds

where $s$ is an arbitrary constant greater than $1$. Repeating the sampling procedure for $t$ times and letting $X_{(i)}$ 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öm 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-\delta$.

For instance, we let $s=1+\log({1}/{\delta})$, then by repeating the sampling procedure for $t\geq 1+1/{\epsilon}$ times, the inequality

holds with probability at least $1-\delta$.

For another instance, we let $s=2$, then by repeating the sampling procedure for $t\geq(1+1/{\epsilon})\log(1/\delta)$ times, the inequality

holds with probability at least $1-\delta$.

Empirical Analysis

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öm and our modified Nyström, 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

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.

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

2 Comparison among the Nyström 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öm and our modified Nyström, 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 ${\bf A}$ which is defined by

where ${\bf x}_{i}$ and ${\bf x}_{j}$ are data instances and $\sigma$ is a scale parameter. Notice that the RBF kernel is dense in general. We set $\sigma=0.2$ or $1$ in our experiments. For each data set with different settings of $\sigma$, we fix a target rank $k=10$, $20$ or $50$ and vary $c$ in a very large range. We will discuss the choice of $\sigma$ and $k$ 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 $k$. In Table 5, $\frac{\|{\bf A}-{\bf A}_{k}\|_{F}}{\|{\bf A}\|_{F}}$ denotes ratio that is not captured by the best rank-$k$ approximation to the RBF kernel, and the parameter $\sigma$ has an influence on the ratio ${\|{\bf A}-{\bf A}_{k}\|_{F}}/{\|{\bf A}\|_{F}}$. When $\sigma$ is large, the RBF kernel can be well approximated by a low-rank matrix, which implies that (i) a small $k$ suffices when $\sigma$ is large, and (ii) $k$ should be set large when $\sigma$ is small. So the settings ($\sigma=1$, $k=10$) and ($\sigma=0.2$, $k=50$) are more reasonable than the rest. Let us take the RBF kernel in the Abalone data set as an example. When $\sigma=1$, the rank-$10$ approximation well captures the kernel, so $k$ can be safely set as small as $10$; when $\sigma=0.2$, 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 $\sigma$ influences the homogeneity/heterogeneity of the leverage scores. Usually, when $\sigma$ is 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 RBF kernel with the scale parameter $\sigma=0.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öm 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 ${\mathcal{O}}(m^{2}k)$ time to compute the truncated SVD.

Furthermore, the experimental results show that using ${\bf U}={\bf C}^{\dagger}{\bf A}({\bf C}^{\dagger})^{T}$ as the intersection matrix (denoted by “modified” in the figures) always leads to much lower error than using ${\bf U}={\bf W}^{\dagger}$ (denoted by “standard”). However, our modified Nyström method costs more time to compute the intersection matrix than the standard Nyström method costs. Recall that the standard Nyström costs ${\mathcal{O}}(c^{3})$ time to compute ${\bf U}={\bf W}^{\dagger}$ and that the modified Nyström costs ${\mathcal{O}}(mc^{2})+T_{\textrm{Multiply}}(m^{2}c)$ time to compute ${\bf U}={\bf C}^{\dagger}{\bf A}({\bf C}^{\dagger})^{T}$. 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öm method.

Conclusion

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öm approximation algorithms which demonstrate significant improvement over the classical counterparts. Our relative-error CUR algorithm requires only $c={2k}{\epsilon^{-1}}(1+o(1))$ columns and $r={c}{\epsilon^{-1}}(1{+}\epsilon)$ rows selected from the original matrix. To achieve relative-error bound, the best previous algorithm—the subspace sampling algorithm—requires $c={\mathcal{O}}(k\epsilon^{-2}\log k)$ columns and $r={\mathcal{O}}(c\epsilon^{-2}\log c)$ rows. Our modified Nyström method is different from the conventional Nyström methods in that it uses a different intersection matrix. We have shown that our adaptive sampling algorithm for the modified Nyström achieves relative-error upper bound by sampling only $c={2k}{\epsilon^{-2}}(1{+}o(1))$ columns, which even beats the lower error bounds of the standard Nyström and the ensemble Nyström. Our proposed CUR and Nyström 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.

A The Dual Set Sparsification Algorithm

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.

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: {\bf A}_{\tau}={\bf W}\mbox{\boldmath\Lambda\unboldmath}{\bf W}^{T}. Here ${\bf A}_{\tau}$ is guaranteed to be SPSD in each iteration. Since

$({\bf A}_{\tau}-(L_{\tau}+1){\bf I}_{k})^{q}$ can be efficiently computed based on the eigenvalue decomposition of ${\bf A}_{\tau}$. With the eigenvalues at hand, $\phi(L,{\bf A}_{\tau})$ can also be computed directly.

B Proofs of the Adaptive Sampling Bounds

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.

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.

where the expectation is taken w.r.t. ${\bf C}_{2}$.

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.

Notice that ${\bf x}_{j,(l)}$ is a linear function of a column of ${\bf A}$ sampled from the above defined distribution. We have that

Then we let ${\bf x}_{j}=\frac{1}{c_{2}}\sum_{l=1}^{c_{2}}{\bf x}_{j,(l)}$, we have

The expectation of $\|{\bf w}_{j}-{\bf A}{\bf v}_{j}\|_{2}^{2}$ is

We use ${\bf F}$ to bound the error $\|{\bf A}{\bf R}^{\dagger}{\bf R}-{\bf C}{\bf C}^{\dagger}{\bf A}{\bf R}^{\dagger}{\bf R}\|_{F}^{2}$. That is,

where (B.1) is due to that ${\bf A}({\bf I}-{\bf R}^{\dagger}{\bf R})$ is orthogonal to $({\bf I}-{\bf C}{\bf C}^{\dagger}){\bf A}{\bf R}^{\dagger}{\bf R}$. Since ${\bf A}{\bf R}^{\dagger}{\bf R}$ and ${\bf F}$ both lie on the space spanned by the right singular vectors of ${{\bf A}{\bf R}^{\dagger}{\bf R}}$ (i.e., $\{{\bf v}_{j}\}_{j=1}^{\rho}$), we decompose ${\bf A}{\bf R}^{\dagger}{\bf R}-{\bf F}$ along $\{{\bf v}_{j}\}_{j=1}^{\rho}$, obtaining that

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

B.2 The Proof of Corollary 7

and $c/r_{2}=\epsilon$. It then follows from Theorem 5 that

which yields the error bound for CUR matrix decomposition.

When the matrix ${\bf A}$ is symmetric, the matrix ${\bf C}_{1}^{T}$ consists of the rows ${\bf A}$, and thus we can use Theorem 15 (which is identical to Theorem 5) to prove the error bound for the Nyström approximation. By replacing ${\bf R}$ in Theorem 15 by ${\bf C}_{1}^{T}$, we have that

where the expectation is taken w.r.t. ${\bf C}_{2}$. Together with the inequality

Given an $m{\times}m$ matrix ${\bf A}$ and an $m{\times}c$ matrix ${\bf C}=[{\bf C}_{1},{\bf C}_{2}]$, the following inequality holds:

B.3 The Proof of Theorem 8

B.4 The Proof of Theorem 10

C Proofs of the Lower Error Bounds

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öm method. In Appendix C.3 we prove the lower bounds of the ensemble Nyström method. Theorems 20, 21, 22, 24, and 25 are used for proving Theorem 12.

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 $m{\times}m$ positive definite matrix ${\bf B}$ as follows:

Let ${\bf B}_{k}$ be the best rank-$k$ approximation to the matrix ${\bf B}$ defined in (9). Then we have that

Proof The squared Frobenius norm of ${\bf B}$ is

Then we study the singular values of ${\bf B}$. Since ${\bf B}$ is SPSD, here we do not distinguish between its singular values and eigenvalues.

The spectral norm, that is, the largest singular value, of ${\bf B}$ is

where the maximum is attained when ${\bf x}=\frac{1}{\sqrt{m}}{\bf 1}_{m}$. Thus ${\bf u}_{1}=\frac{1}{\sqrt{m}}{\bf 1}_{m}$ is the top singular vector of ${\bf B}$. Then the projection of ${\bf B}$ onto the subspace orthogonal to ${\bf u}_{1}$ is

Then for all $j>1$, the $j$-th top eigenvalue $\sigma_{j}$ and eigenvector ${\bf u}_{j}$, that is, the singular value and singular vector, of ${\bf B}$ satisfy

where the last equality follows from ${\bf u}_{j}\perp{\bf u}_{1}$, that is, ${\bf 1}_{m}^{T}{\bf u}_{j}=0$. Thus $\sigma_{j}=1-\alpha$, and

for all $1\leq k<m$. Finally we have that

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 ${\bf B}$ be a $p\times p$ matrix with diagonal entries equal to one and off-diagonal entries equal to $\alpha$. Let $m=kp$ and we construct an $m\times m$ block diagonal matrix ${\bf A}$ as follows:

Let ${\bf A}_{k}$ be the best rank-$k$ approximation to the matrix ${\bf A}$ defined in (14). Then we have that

Lemma 19 can be easily proved using Lemma 18.

C.2 Lower Bounds of the Standard Nyström Method

For an $m\times m$ matrix ${\bf B}$ with diagonal entries equal to one and off-diagonal entries equal to $\alpha\in[0,1)$, the approximation error incurred by the standard Nyström method is lower bounded by

Proof The matrix ${\bf B}$ is partitioned as in (9). The residual of the Nyström approximation is

where $\xi=2$, $F$, or $*$. Since ${\bf W}=(1-\alpha){\bf I}_{c}+\alpha{\bf 1}_{c}{\bf 1}_{c}^{T}$ is nonsingular when $\alpha\in[0,1)$, so ${\bf W}^{\dagger}={\bf W}^{-1}$. We apply the Sherman-Morrison-Woodbury formula

to compute ${\bf W}^{{\dagger}}$, yielding

According to the construction, ${\bf B}_{21}$ is an $(m{-}c)\times c$ matrix with all entries equal to $\alpha$, it follows that ${\bf B}_{21}{\bf W}^{\dagger}{\bf B}_{21}^{T}$ is an $(m{-}c){\times}(m{-}c)$ matrix with all entries equal to

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 $\frac{1}{\sqrt{m-c}}{\bf 1}_{m-c}$ is the top singular vector of the SPSD matrix $(1-\alpha){\bf I}_{m{-}c}+(\alpha-\eta){\bf 1}_{m{-}c}{\bf 1}_{m{-}c}^{T}$, so the top singular value is

It is also easy to show the rest singular values obey

The theorem follows from equalities (18), (19), and (20).

Now we use the matrix ${\bf A}$ constructed in (14) 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$.

For the $m\times m$ SPSD matrix ${\bf A}$ defined in (14), the approximation error incurred by the standard Nyström method is lower bounded by

where $k<m$ is an arbitrary positive integer.

Proof Let ${\bf C}$ consist of $c$ column sampled from ${\bf A}$ and $\hat{{\bf C}}_{i}$ consist of $c_{i}$ columns sampled from the $i$-th block diagonal matrix in ${\bf A}$. Without loss of generality, we assume $\hat{{\bf C}}_{i}$ consists of the first $c_{i}$ columns of ${\bf B}$, and accordingly $\hat{{\bf W}}_{i}$ consists of the top left $c_{i}\times c_{i}$ block of ${\bf B}$. Thus {\bf C}=\mathsf{BlkDiag}\big{(}\hat{{\bf C}}_{1},\cdots,\hat{{\bf C}}_{k}\big{)} and {\bf W}=\mathsf{BlkDiag}\big{(}\hat{{\bf W}}_{1},\cdots,\hat{{\bf W}}_{k}\big{)}.

where $\hat{p}=p+\frac{1-\alpha}{\alpha}$ and $\hat{c_{i}}=c_{i}+\frac{1-\alpha}{\alpha}$. Since $\sum_{i=1}^{k}\hat{c}_{i}=c+\frac{1-\alpha}{\alpha}k\triangleq\hat{c}$, the term $\sum_{i=1}^{k}\hat{c}_{i}^{-2}$ is minimized when $\hat{c}_{1}=\cdots=\hat{c}_{k}$. Thus $\sum_{i=1}^{k}\hat{c}_{i}^{-2}=k\frac{k^{2}}{\hat{c}^{2}}=k^{3}\hat{c}^{-2}$. Finally we have that

by which the Frobenius norm bound follows.

where the former inequality follows from Theorem 20, and the latter inequality follows by minimizing w.r.t. $c_{1},\cdots,c_{k}$ subjecting to $c_{1}+\cdots+c_{k}=c$.

There exists an $m{\times}m$ SPSD matrix ${\bf A}$ such that the approximation error incurred by the standard Nyström method is lower bounded by

where $k<m$ is an arbitrary positive integer.

Proof For the spectral norm bound we use the matrix ${\bf A}$ constructed in (9) and set $\alpha\rightarrow 1$, then it follows directly from Lemma 18 and Theorem 20. For the Frobenius norm and nuclear norm bounds, we use the matrix ${\bf A}$ constructed in (14) and set $\alpha\rightarrow 1$, then it follows directly from Lemma 19 and Theorem 21.

C.3 Lower Bounds of the Ensemble Nyström Method

The ensemble Nyström method (Kumar et al., 2009) is previously defined in (2). To derive lower bounds of the ensemble Nyström method, we assume that the $t$ samples are non-overlapping. According to the construction of the matrix ${\bf 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: $\mu^{(1)}=\cdots=\mu^{(t)}=\frac{1}{t}$.

Assume that the ensemble Nyström method selects a collection of $t$ samples, each sample ${{\bf C}^{(i)}}$ ($i=1,\cdots,t$) contains $c$ columns of ${\bf B}$ without overlapping. For an $m\times m$ matrix ${\bf B}$ with all diagonal entries equal to one and off-diagonal entries equal to $\alpha\in[0,1)$, the approximation error incurred by the ensemble Nyström method is lower bounded by

Proof We use the matrix ${\bf B}$ constructed in (9). It is easy to check that ${{\bf W}^{(1)}}=\cdots={{\bf W}^{(t)}}$, so we use the notation ${{\bf W}}$ instead. We assume that the samples contain the firs $tc$ columns of ${\bf B}$ and each sample contains neighboring columns, that is,

If a sample ${{\bf C}}$ contains the first $c$ columns of ${\bf B}$, then

otherwise, after permuting the rows and columns of ${\bf B}-{\bf C}{\bf W}^{\dagger}{{\bf C}}^{T}$, we get the same result:

where $\Pi$ is a permutation matrix. As was shown in Equation (16), ${\bf B}_{21}{\bf W}^{\dagger}{\bf B}_{21}^{T}$ is an $(m{-}c){\times}(m{-}c)$ matrix with all entries equal to

where the last inequality follows from \frac{c(m-c)}{t}=\frac{c}{t}\Big{(}(m-2c+\frac{c}{t})+(c-\frac{c}{t})\Big{)}\geq\frac{c}{t}\Big{(}m-2c+\frac{c}{t}\Big{)}.

which proves the nuclear norm bound in the lemma.

Assume that the ensemble Nyström method selects a collection of $t$ samples, each sample ${{\bf C}^{(i)}}$ ($i=1,\cdots,t$) contains $c$ columns of ${\bf A}$ without overlapping. For a the matrix ${\bf A}$ defined in (14), the approximation error incurred by the ensemble Nyström method is lower bounded by

Proof According to the construction of ${\bf A}$ in (14), the $i$-th sample ${\bf C}^{(i)}$ is also block diagonal. We denote it by {\bf C}^{(i)}=\mathsf{BlkDiag}\big{(}\hat{{\bf C}}^{(i)}_{1},\cdots,\hat{{\bf C}}^{(i)}_{k}\big{)}. Akin to (44), we have

Thus the approximation error of the ensemble Nyström method is

where the inequality follows from Lemma 23, and the last equality follows from $\sum_{j=1}^{k}c_{j}=c$ and $kp=m$. The summation in the last equality equals to

Here the inequality holds because the function is minimized when $c_{1}=\cdots=c_{k}=c/k$. Finally we have that

which proves the Frobenius norm bound in the theorem.

which proves the nuclear norm bound in the theorem.

Assume that the ensemble Nyström method selects a collection of $t$ samples, each sample ${{\bf C}^{(i)}}$ ($i=1,\cdots,t$) contains $c$ columns of ${\bf A}$ without overlapping. Then there exists an $m{\times}m$ SPSD matrix ${\bf A}$ such that the relative-error ratio of the ensemble Nyström method is lower bounded by

Proof The theorem follows directly from Theorem 24 and Lemma 19 by setting $\alpha\to 1$.

References