# Apply PCA on very large sparse matrix

I am doing a text classification task with R, and I obtain a document-term matrix with size 22490 by 120,000 (only 4 million non-zero entries, less than 1% entries). Now I want to reduce the dimensionality by utilizing PCA (Principal Component Analysis). Unfortunately, R cannot handle this huge matrix, so I store this sparse matrix in a file in the "Matrix Market Format", hoping to use some other techniques to do PCA.

So could anyone give me some hints for useful libraries (whatever the programming language), which could do PCA with this large-scale matrix with ease, or do a longhand PCA by myself, in other words, calculate the covariance matrix at first, and then calculate the eigenvalues and eigenvectors for the covariance matrix.

What I want is to calculate all PCs (120,000), and choose only the top N PCs, who accounts for 90% variance. Obviously, in this case, I have to give a threshold a priori to set some very tiny variance values to 0 (in the covariance matrix), otherwise, the covariance matrix will not be sparse and its size would be 120,000 by 120,000, which is impossible to handle with one single machine. Also, the loadings (eigenvectors) will be extremely large, and should be stored in sparse format.

Thanks very much for any help !

Note: I am using a machine with 24GB RAM and 8 cpu cores.

• How sparse is the matrix? How do you use the resulting SVD? If you only need part of it you could probably approximate it far cheaper. – Arnold Neumaier May 23 '12 at 13:36
• @ArnoldNeumaier Excuse me, I forgot to add the sparse info. I've updated the post, together with my complete idea. – Ensom Hodder May 23 '12 at 15:34
• each of SLEPc, mahout and irlba suggested in the answers so far seem suitable for your problem. – Arnold Neumaier May 23 '12 at 18:32
• Why do you want to compute all 120k? It sounds like you just want those accounting for 90% of the variance, which should be much cheaper to compute. – Jed Brown May 24 '12 at 3:11
• @JedBrown Hey Jed, you're totally right! I am only interested on those who account for 90% variance, and also corresponding eigenvectors (for transforming the test dataset afterwards). Could you please let me know your cheaper methods ? – Ensom Hodder May 24 '12 at 8:39

I suggest the irlba package - it produces virtually the same results as svd, yet you can define a smaller number of singular values to solve for. An example, using sparse matrices to solve the Netflix prize, can be found here: http://bigcomputing.blogspot.de/2011/05/bryan-lewiss-vignette-on-irlba-for-svd.html

• Thanks for your comments. In fact, I'd watched that video and also tried irlba package yesterday, but it seemed that it could be used only for calculating a few singular values. However, as stated in the post, I want to calculate ALL singular values (120,000), so as to choose suitable number of PCs according to the variances they account for. In this case, I guess irlba is no longer suitable. – Ensom Hodder May 23 '12 at 19:46
• Can you use the results of SVD in a similar manner to PCA? Don't you need to center the data BEFORE doing the SVD, in order to perform PCA? – Zach Aug 28 '12 at 18:29
• @Zach - SVD is the main algorithm behind PCA (see prcomp - stat.ethz.ch/R-manual/R-patched/library/stats/html/prcomp.html). Centering of data is also standard procedure before subjecting to PCA, although there are a wide variety of options depending on your question (e.g. different types of scaling may also be applied). – Marc in the box Aug 29 '12 at 4:14
• How big of a deal is it if I do not center the data before SVD? I have a sparse matrix that fits into memory, but centering would make it dense and too big to fit into memory. – Zach Aug 29 '12 at 12:18
• @Zach - It really depends on how you want to relate your samples to each other. If you cant work with centered data due to memory limits, then I guess the decision has been made up for you. Generally, centering data has the PCA operate on a covariance matrix of the samples while centering and scaling of data has PCA operate on a correlation matrix. For more insight into these decisions, you may consider asking a question over at stats.stackexchange.com or search through the existing answers regarding PCA. – Marc in the box Aug 30 '12 at 5:41

I suggest using SLEPc to compute a partial SVD. See Chapter 4 of the User's Manual and the SVD man pages for details.

• Since he wants PCA he must center the data before computing the SVD. This will wreck the sparsity. Is there some way that SLEPc accommodates for this? – dranxo May 24 '12 at 2:53
• That is just sparse + low rank. SLEPc does not need matrix entries, just a linear operator, which can be applied as a sparse matrix plus a correction. – Jed Brown May 24 '12 at 3:09

I vote for mahout which is also good for other NLP/TA tasks and implements map/reduce.

• Yes, you are right, mahout is exactly in my road map. But I prefer to create a prototype with some "simple"(i suppose) techniques in advance. – Ensom Hodder May 23 '12 at 14:59

I'd suggest using an incremental singular value decomposition, of which there are many in the literature. For instance:

• the technical reports of Matthew Brand 1 and 2 are fairly easy to follow
• Chris Baker's master's thesis, his software IncPACK, and his later paper on incremental SVD method
• Bunch and Nielsen published the earliest known paper
• Papers by Hall on updating eigenvalue problems 1 and 2
• Sequential Karhunen-Loeve analysis by Levy, et al., which is basically the same thing

All of these approaches reduce to the following:

• calculate an SVD somehow (this step is trivial for a single column matrix)
• repeat until finished:
In your application, if you have an idea of where your singular value threshold for the top $N$ values will be, you can use that value to calculate a truncated SVD; if the threshold value is small enough, then the matrix you have to keep in memory will also be small (only the singular values above the threshold value are retained, along with their singular vectors; it's not even necessary to keep both left and right singular vectors, in Brand's algorithm).
Revolution R is a build of R that handles data sets that are larger than RAM. Use the function princomp.