A Scalable Approach to Column-Based Low-Rank Matrix Approximation

1 Introduction

In multivariate data analysis, the data in question is usually defined as an m ⇥ n feature-instance matrix; that is, each column represents an instance of m features. Examples of feature-instance datasets include environmental measurements and genomes, words and documents, face images and people, etc. Typically, finding a parsimonious approximation of a data matrix assists in the efficient recognition of valuable information behind the instances. A classical approach is to use the singular value decomposition (SVD) to obtain a low rank approximation of the matrix, such as in the case of eigenface and eigengene.

Although the SVD method is widely used for low-rank approximation of a data matrix, it is hard to assign the resulting singular vectors meaning in terms of the field from which the data is generated. This makes the task of understanding and interpreting the data in question very difficult. Compared to SVD, column-based approximation methods, representing a data matrix in terms of a small number of actual columns and/or actual rows of the matrix, have been shown to have advantage of analyzing data [Goreinov et al., 1997a; 1997b; Tyrtyshnikov, 2000; Drineas, 2003; Drineas and Mahoney, 2005; Drineas et al., 2006].

⇤Z.Z. is the corresponding author, and he acknowledges support from the Natural Science Foundations of China (No. 61070239) and the Google visiting faculty program. This work was done when Y.P. and H.P. were research interns in Google Beijing. Y.P. was supported in part by the National Basic Research Program of China Grant 2011CBA00300, 2011CBA00301, the National Natural Science Foundation of China Grant 61033001, 61061130540.

Many algorithms, both randomized and deterministic, have been developed to find the proper columns of the columnbased approximation [Boutsidis et al., 2009; Kannan and Vempala, 2009; Boutsidis et al., 2011; Mahoney, 2011; C¸ ivril and Magdon-Ismail, 2012; Guruswami and Sinop, 2012]. Drineas et al. [2008] devised randomized algorithms with relative error by sampling sufficiently many columns. In particular, Drineas et al.'s algorithm has (1+✏) relative-error ratio with high probability (w.h.p.). Unfortunately, relativeerror bound algorithms for a general matrix often require the use of SVD as a subprocess when performing computation on the data matrix. This limits scalability of such algorithms.

To solve the scalability problem resulted from SVD, numerous algorithms have been proposed. While some approaches focus on speeding up underlying numerical iterations [Golub and Loan, 1996], most are numerically unstable and hard to implement. Randomized algorithms with fast running times and robust performance have enjoyed much interest in the research community, e.g., see Halko et al. [2011] for a review. A recent unpublished work claims to approximate low-rank SVD randomly in input sparsity time with relative error bound [Clarkson and Woodruff, 2012].

We propose a novel method to perform column-based approximation on large scale data matrices. Our work is motivated by the divide-and-combine methodology, which was also studied by Mackey et al. [2011] when improving scalability of both matrix completion and robust PCA. Informally, given some conditions, we divide the original matrix into smaller submatrices, factorize these submatrices in parallel, then combine the results from the factorization to obtain a final matrix consistent with the conditions (in this work, our conditions correspond to some specific number of columns).

Our framework is fully deterministic as long as the base algorithms used are deterministic. Our approach does not require the data matrix to satisfy any coherence conditions, so it applies to general matrices. We note that the size of matrices becomes m⇥ dpnke rather than m⇥ n, where k is a low rank parameter which is often regarded as a small constant in practice [Mahoney and Drineas, 2009]. Thus, our framework is scalable. In addition, any sparse input matrix can be handled by our algorithm flexibly. During the execution of the algorithm, only small-size dense matrices will be produced at any given step.

The rest of this paper is organized as follows. In Section 2 we present some preliminary backgrounds. In Section 3 we describe our approach. Section 4 establishes a relative error bound on general matrices. Section 5 gives an empirical analysis of our approach. Finally, we conclude our work in Section 6.

2 Preliminaries

This section presents some notational conventions that we will use in this paper. The norm function k · k⇠ can be either the Frobenius norm (denoted k · kF) or the spectral norm (denoted k · k2). Let In be the n ⇥ n identity matrix. For a matrix M 2 Rm⇥n , we denote its rank by rank(M) and its Moore-Penrose pseudoinverse by M+. Let σ1(M) ≥ σ2(M) ≥ ▷▷▷ ≥ σr(M) > 0 be r positive singular values of M. Let ui(M) 2 Rm and vi(M) 2 Rn be respectively the left singular vector and right singular vector of M associated with σi for all 1  i  r. Given a matrix C 2 Rm⇥r with r > k, ⇧⇠C↪k (M) minimizes the residual kM- cMk⇠ over all cM in the column space of C that have rank at most k. For column submatrices, we follow the Matlab notation: M(:↪J) represents the submatrix formed by columns in index set J of M. Note that we also use C and B to represent column submatrices if the base matrix is implicit according to the context. We often use k as the rank

parameter and r as the number of columns to be chosen.

We now introduce three matrix factorization algorithms that we will refer to as base algorithms for short.

2.1 Truncated Singular Value Decomposition

Singular Value Decomposition (SVD) is often used when we want to have a concise approximation of a data matrix. It is a well-known result that for M 2 Rm⇥n and k  min¶m↪ n♦,

The classical approach to such a best rank-k approximation of M is to truncate the result of the full SVD to obtain the first top k singular values and singular vectors. The running time for this approach is O(mnmin¶m↪ n♦). There are also other numerical methods such as the Krylov subspace method

[Golub and Loan, 1996] to calculate the best rank-k approximation in O(mnk) time.

2.2 Column-basedMatrix Factorization (CX)

In contrast to the k-rank approximation of the matrix in the SVD approach, the CX approach focuses on choosing a subset of r columns for some k  r  n, whose span should approximate M as close as possible to Mk. Suppose those columns we choose form a matrix C 2 Rm⇥r . We aim to reconstruct M via the column space of C. We note that a common way of reconstructing the matrix M is to compute CC+M, but the rank of this matrix may be larger than k. In order to ensure the reconstructed matrix also has rank k, Boutsidis et al. [2011] devised methods to find a matrix ⇧⇠C↪k (M). Unless mentioned explicitly, we assume the reconstructed matrix to be CC+M below. However, in our theoretical analysis on the error upper bound, we consider the case when the reconstructed matrix is ⇧⇠C↪k (M) because

Letting C be the matrix consisting of r columns of M, we want to relate either kM-CC+Mk⇠ or kM-⇧⇠C↪k (M)k⇠ to kM-Mkk⇠. Based on the relative-error approximation, we have kM-CC+Mk⇠  kM-⇧⇠C↪k (M)k⇠  ∆kMMkk⇠ for some ∆ which is a function of r, k and the size of the matrix. If an algorithm is randomized, the above bound is obtained in the expectation sense. Note that as r becomes larger, ∆ becomes smaller. This gives rise to a more accurate approximation. Please see Figure 1 of Guruswami and Sinop [2012] for comparison on the error bounds and the running times of various CX algorithms.

There are also algorithms that yield an additive error bound [Kannan and Vempala, 2009], i.e., kM-Mkk⇠ + "kMk⇠. These algorithms run in linear time by using some trivial sampling methods, but they result in a huge error when kMk⇠ is much larger than kM-Mkk⇠. In this paper we focus on algorithms that produce a relative-error bound.

2.3 Interpolative Decomposition (ID)

Interpolative Decomposition (ID) selects a collection of k columns from a rank-k matrix M 2 Rm⇥n that span the column space of M. Formally, ID finds a column index set J = ¶j1↪ ▷▷↪ jk♦ such that M = M(:↪J)X↪ where X 2 Rk⇥n satisfies X(:↪J) = Ik and has no entry with magnitude larger than two. When the Gu-Eisenstat deterministic algorithm is used to produce ID [Gu and Eisenstat, 1996; Cheng et al., 2005], the running time is close to O(mnk) for an m ⇥ n rank-k matrix, though it can be higher (i.e., O(mnmin¶m↪ n♦)) in rare cases [Halko et al., 2011].

3 Methodology

In this section, we first give the framework of our divide-andcombine column-based matrix approximation method, and then present its running time analysis and discuss the related work.

3.1 Framework

Given a matrix M 2 Rm⇥n , our purpose is to choose r columns to approximate Mk (k  rank(M), k  r  n). Our method consists of four stages which are described as follows.

Divide Given an integer t > 0 such that bn◁tc ≥ k, we split M into t l-column submatrices in any order 1 , i.e., M = [M(1) ↪ ▷▷▷↪M(t) ]. Unlike the method of Mackey et al. [2011], our approach does not need to randomly permute the order of columns. Thus, determinism of our approach only depends on the base algorithms.

• Approximation We perform t rank-k approximation algorithms (e.g., the best rank-k approximation through SVD) on ¶M(1) ↪ ▷▷▷↪M(t) ♦ in parallel. It yields t rank-k
• corresponds to columns in cM(i) . Ensemble We first form an ensemble matrix B = M(:↪J1[▷▷▷[Jt) 2 Rm⇥s where s = kt. 2 Next, we run the column based base approximation algorithm on B. 3 We accordingly have B ⇡ CC+B where C contains r columns of M. C is the final output, i.e., M ⇡ CC+M.
• matrices ¶cM(1) ↪ ▷▷▷↪ cM(t) ♦. Selection In parallel, we run the ID algorithm on each cM(i) (a total of t parallel instances). This results in, each i, cM(i) = bC(i) bX(i) where bX(i) is a k ⇥ l matrix such that no entry of bX(i) has magnitude larger than two. We let Ji denote the set of indices of the columns in bC(i) that
• Adaptive Sampling (Optional) The purpose of this step is to sample more columns to C to obtain a tighter theoretical error bound. We can omit this step in practice. This step was introduced in [Deshpande et al., 2006]. At first, we computeM0 asM0 =M-CC+M. Since C is a small matrix, the pseudoinverse of C and the matrix multiplications can be computed very fast. Let pi = km0i k22 ◁kM0 k2F where m0i is the i-th column of M0 . Then, sample further s0 columns fromMin s0 i.i.d. trials, where in each trial the i-th column is chosen with probability pi. Finally, add these s0 columns to C to obtain the new output. The choice of the parameter s0 is discussed in Section 4.

Note that the rank-k approximation of M(i) is generally returned in factorized form, i.e., cM(i) = bR(i) bT(i) where bR(i) 2 Rm⇥k and bT(i) 2 Rk⇥l . Thus, when the ID algorithm is performed on bT(i) , we find an index set Ji such that bT(i) = bT(i) (:↪Ji) bX(i) . From this, we know cM(i) = bR(i) bT(i) (:↪Ji) bX(i) = cM(i) (:↪Ji) bX(i) = bC(i) bX(i) . The trick was devised by Halko et al. [2011] to perform the ID on the smaller matrix. We also employ this trick in the selection stage.

2 s may be less than r if t is small. If so, we just choose whole B, yielding the same error bound in the theoretical analysis.

1 For simplicity, assume that n mod t = 0, and hence, l = n◁t. Note that for an arbitrary n and t, M can always be partitioned into t submatrices, each with either bn◁tc or dn◁te columns.

3 Ideally, we should use bB = [ bC1↪ ▷▷▷↪ bCt]. However, the bCi's are not real columns of M. In the theoretical analysis, we will show that even our approach works on the B, we obtain a good approximation.

Additionally, note that while B in the ensemble stage can be used as prediction of the column space of M, it may contain too many columns and become large when M is split into more matrices. Thus, in the ensemble stage B should be resampled to obtain r columns in order to maintain consistency with the base CX algorithm.

It is worth pointing out that our framework can also handle the case that input matrix M is sparse. Particularly, in the approximation stage, we can employ a suitable algorithm for the sparse matrix. As for the selection stage, there is no problem because we only need to do ID on a small dense matrix bT(i) as mentioned earlier. Finally, in the ensemble stage, the matrix B is sparse because it is a submatrix of M.

3.2 Running Time Analysis

Let the running time be TC(m↪ n) for the base CX algorithm, TI(m↪ n) for the ID, and TL(m↪ n) for a base rank-k matrix approximation method. Then the total running time for choosing r columns to approximate M 2 Rm⇥n of rank k is

Note that by introducing parallelism the factor t in the term t[TL(m↪ l)+TI(k↪ l)] becomes negligible. In practice, our approach extends the CX algorithms to also work on large matrices that do not fit into the main memory. Our algorithm performs extremely fast with a relative-error bound when k ⌧ min¶m↪ n♦ (which is common in practice). For example, when t = dpn◁ke, we see that all matrices we need to handle are of size at most m⇥dpkne, much smaller than the original m⇥ n matrix. In addition, when kt is still large we can further perform the divide-and-combine strategy on B to obtain an approximation which is also relative-error bounded in terms of Theorems 1 and 2 given in Section 4.

3.3 RelatedWork

It is worth noting that Halko et al. [2011] proposed a compressed factorization method of large-scale matrices. The method uses a randomized technique to get the sketch of the column space of the matrix. Then it projects the matrix on the sketched space and performs factorization. Our approach keeps the main advantages of this method, while introducing several differences. First, we can use any underlying rank-k approximation (including the method in [Halko et al., 2011] and [Clarkson and Woodruff, 2012]), ID and CX algorithms. Second, our approach does not depend on randomness. Finally, our approach is easily parallelized and resolves the memory problem of a large-scale matrix.

There is also an approximation method so-called CUR similar to CX [Drineas et al., 2008]. Our approach adapts to the calculation of CUR decomposition by doing CX on M and M>, respectively, to obtain a column submatrix C and a row submatrix R. We approximate M with CC+MR+R = CUR. A similar error bound to the CUR decomposition under our approach is easily deduced by the equation (6) in [Mahoney and Drineas, 2009] when a bound of the CX algorithm is given. We note that the bound on CUR may be a factor of 2 larger than that of the CX algorithm.

4 Theoretical Analysis

In this section, we present theoretical results that establish a relative error bound of our approach. We first give a brief discussion on the type of the bound. Next, we introduce a main theorem about the error bound. Finally, we analyze our framework. The analysis shows that our approach enjoys a relative-error bound as long as base algorithms used also give relative-error bounds.

4.1 Main Results

In this subsection, we present two theorems, one based on the Frobenius norm and the other based on the spectral norm. The theorems show that the ensemble stage can obtain an accurate result on a block low rank matrix when the base CX algorithm yields a bounded approximation.

Theorem 1. Let t = n◁l and l ≥ k. Suppose L = [L(1) ↪ ▷▷▷↪L(t) ] 2 Rm⇥n where L(i) 2 Rm⇥l and rank(L(i) )  k. We apply an ID algorithm on L(i) to obtain L(i) = bC(i) bX(i) for i = 1↪ ▷ ▷ ▷ ↪ t. Denote bB = [ bC(1) ↪ ▷▷▷↪ bC(t) ] 2 Rm⇥s . For a fixed matrix B and a random matrix C 2 Rm⇥r , k  r  s, and EC[kbB-⇧FC↪k (B)kF]  ∆F , then EC[kL - ⇧FC↪k (L)kF]  2pkl∆F .

maxi k bX(i) kF erty of the spectral norm. (b) follows from the fact that the magnitude of entries in X(i) is bounded by 2.

Theorem 2. We use the same notation as in Theorem 1. If C holds the property that EC[kbC-⇧2C↪k (bC)k2]  ∆S, then EC[kL - ⇧2C↪k (L)k2]  2pkl∆S .

Once we replace the Frobenius norm with the spectral norm, the proof is immediately obtained from that of Theorem 1.

4.2 Theoretical Analysis for CX

We also give two theorems respectively based on the Frobenius norm and the spectral norm, which ensure our approach yields a relative-error bound. Note that ∆⇤ stated in the following theorems is the error parameter given by the base algorithm which depends on r, k and the size of the data matrix. Please refer to Figure 1 in [Guruswami and Sinop, 2012] and [Halko et al., 2011] for examples.

Theorem 3. Let t = n◁l and l ≥ k. Based on the same notation as in Section 3.1, the base rank-k approximation algorithm returns a matrix cM(i) of rank at most k such that EcM(i) [kM(i) - cM(i) k2F ]  ∆2LF kM(i) -M(i) k k2F for some ∆LF ≥ 1. The base CX algorithm yields EC[kB - ⇧FC↪k (B)kF♣B]  ∆CFkB-BkkF for ∆CF ≥ 1 by choosing r columns C from B. Then EC↪cM[kM-⇧FC↪k (M)kF]  [2pkl(∆LF + ∆CF ) + ∆LF ]kM - MkkF.

Proof. Split Mk = [N(1) ↪ ▷▷▷↪N(t) ] as the same way as M= [M(1) ↪ ▷▷▷↪M(t) ]. and note that rank(N(i) )  k for each i. Then, EcM[kM - cMk2F ] = EcM hPti=1 kM(i) - cM(i) k2F i

When cM is fixed, we know bB and B are determined because the ID algorithm is deterministic (see Section 2.3). Consider the error introduced by replacing bB with B, and perform the CX approximation on B. That is,

 Pti=1 ∆2LF kM(i) - M(i) k k2F (a)  ∆2LF Pti=1 kM(i) - N(i) k2F = ∆2LF kM-Mkk2F . (a) follows because M(i) k is the best rank-k approximation of M(i) , and rank(N(i) )  k. It follows that EcM[kM- cMkF]  qEcM[kM-cMk2F ]  ∆LF kM - MkkF.

 kM-cMkF +∆CFkM-MkkF. (a) follows because B and bB are the submatrix of M and cM respectively. (b) follows because B is always the submatrix of M. If we apply Theorem 1 on block low-rank matrix cM = [cM(1) ↪ ▷▷▷↪cM(t) ] and use the bound on EC[kbB⇧FC↪k (B)kF♣cM], we have EC[kcM - ⇧FC↪k (M)♣cMkF]  2pkl(kM-cMkF +∆CFkM-MkkF). Finally,

(a) follows from rank(⇧FC↪k (cM))  k. Theorem 4. Use the same settings as in Theorem 3, the base rank-k approximation algorithm returns a matrix cM(i) of rank at most k such that EcM(i) [kM(i) - cM(i) k22 ]  ∆2LS kM(i) -M(i) (k) k22 for some ∆LS ≥ 1. If the base CX algorithm yields EC[kB - ⇧2C↪k (B)k2♣B]  ∆CSkB - Bkk2 for ∆CS ≥ 1 by choosing r columns C from B, then EcM↪C[kM - CC+Mk2]  [2pkl(pt∆LS + ∆CS) + pt∆LS ]kM - Mkk2.

The proof is almost identical to that of Theorem 3. We omit the details for the sake of brevity.

4.3 Discussion

We highlight four points on the theoretical results. (1) We give an example to interpret the above theorems intuitively. If we use the best rank-k approximation on submatrices, Theorem 3 shows that we obtain a [2pkl(2 + ") + 1]kM-MkkF upper error bound on the Frobenius norm, where " is the error parameter determined by the base CX algorithm, and l is the number of columns in each submatrix after the divide step. We actually get a relative-error type bound the same as the base CX algorithm. It means if kM-MkkF is small (i.e., the numerical rank of M is k), then our approach approximate the matrix almost exactly. As noted before, the adaptive sampling can be used to obtain a better theoretical bound. By Lemma 16 of Boutsidis et al. [2011], if the parameter s0 is at least k[2pkl(2+")+1]◁", we can make sure that our framework has the same theoretical error bound as the base CX algorithm. For example, let us consider the case that base CX algorithm is Theorem 5 of Boutsidis et al. [2011]. In order to get a (1+")kM-MkkF error bound on the Frobenius norm, we need to choose ⇥(kpkl◁") columns. (2) In general, in order to obtain a (1 + ")kM - MkkF error bound on the Frobenius norm, we should let the parameter s0 in the adaptive sampling be at least k[2pkl(∆LF +∆CF )+∆LF ]◁". (3) There is a trade-off between the approximation error bound and the speed-up. If l is small (i.e., the number of submatrices is large), the error bound is small. But the running time of the Ensemble stage is longer as the matrix B becomes larger. (4) For flexibility, we analyze the approach with randomized algorithms. If all the base algorithms are deterministic, our consequence automatically becomes deterministic (just throw the expectation away). Note that our analysis only gives the error bound in expectation. It can be adapted to the highprobability low-error form by the tail inequality such as the Markov inequality.

5 Empirical Evaluation

The main purpose of this section is to evaluate the accuracy and speed of our approach on four different types of data. In particular, we want to evaluate the speed of our approach and compare the error to the conventional approach that runs the CX algorithm directly on the original data matrices.

5.1 Set-up

For the rank-k approximation algorithm, we use the truncated SVD (SVD command in Matlab with truncation). For the ID algorithm, we choose Algorithm 5 in [Gu and Eisenstat, 1996]. For the CX algorithm, we choose two: a deterministic one (Theorem 2 of Boutsidis et al. [2011], denoted Det), and a randomized one (Algorithm 1 of Drineas et al. [2008], denoted Rand), using the EXACTLY(C) sampling scheme [Drineas et al., 2008] and boosting 20 times. We split the m ⇥ n input matrix into t = dpn◁ke pieces for parallelization. We write the code in Matlab with the PARFOR feature to implement parallelization. We measure the error via the ratio of the Frobenius-norm error over the best rank-k approximation, i.e., reconstruction error ratio ⇥ , kM- CC+MkF◁kM-MkkF. The time measured is wall-clock time. We do not perform the adaptive sampling in all experiments. We evaluate the first three types of data by using Matlab 7.10.

5.2 Simulated Data

We evaluate our framework on two types of random matrices mentioned in [C¸ ivril and Magdon-Ismail, 2012]: Log and Scaled Random. We set m = n = 5000 and the rank parameter k = 15 for the reconstruction error comparison. Note that all reported results are averages over five random examples of the same size. The reconstruction error and the running time are presented in Figure 1 and Table 1 respectively. These figures show that the error is not amplified too much by our framework, and the speed-up is super-linear. Surprisingly, the error of our framework is smaller than the base CX algorithm on Log Data. We think it is due to the property of this type of data.

5.3 Collaborative Filtering Data

A collaborative filtering dataset naturally forms a rating matrix. We use the CX method to approximate the matrix. We evaluate our approach on two public datasets: the Jester dataset [Goldberg et al., 2001] and the MovieLens 10M dataset 4 . To generate suitable matrices, we perform some preprocessing for the two datasets. For the Jester dataset, we only choose the people who rate all 100 jokes to form a 100 ⇥ 14116 matrix. For the MovieLens 10M dataset, we set the unrated entry as a random number uniformly distributed in [0↪ 5] to form a 10K ⇥ 10K matrix. As shown in Figure 2, the error resulted from the framework is almost the same as the direct performance. In addition, the running time is reduced from 130s to 5s for the two base CX algorithms.

4 http://www.grouplens.org

5.4 Human Face Data

In image processing, a matrix can be constructed from a set of images with each image reshaped as a column vector [Peng et al., 2011]. Moreover, such a matrix is usually found to be of low rank. We now use the CX algorithm to approximate the matrix. In particular, we study our approach on the cropped images of the Extended Yale Face Database B [Georghiades et al., 2001; Lee et al., 2005] which forms a 32256 ⇥ 2424 matrix. The performance result is shown in Figure 3. The running time is reduced to 6s from 146s after employing our framework.

5.5 PicasaWeb Dataset

We use images from the PicasaWeb Dataset [Wang et al., 2012]. After having extracted the HMAX features (1000 dimension) [Serre et al., 2007] and deleted some images, we form a 1.44 million by 1000 dense data-feature matrix. We use the truncated SVD as the rank-k approximation algorithm and Rand as the base CX algorithm. In the other word, our algorithm is the same as Framework-Rand. Instead of using Matlab, we implement our framework with the distributed computing framework FlumeJava [Chambers et al., 2010] and the matrix algorithms library Eigen [Guennebaud et al., 2010]. Our method takes less than 10 minutes by 200 12-cores machines to choose 100 columns with k = 15.

6 Conclusion

We have proposed a divide-and-combine approach for handling the scalability problem of column-based approximation methods on large-scale matrices. Our intrinsically parallel approach leads to a significant speed-up on large-scale matrices, and holds a nice error bound theoretically and empirically. Our analysis applies to a general matrix; that is, no additional conditions such as matrix coherence are required. Theoretical analysis and empirical results demonstrate the effectiveness, efficiency and practicality of our approach.