How to compute SVD of a huge sparse matrix? What is the best way to compute singular value decomposition (SVD) of a very large positive matrix (65M x 3.4M) where data is extremely sparse?
Less than 0.1% of the matrix is non zero. I need a way that:


*

*will fit into memory (I know that online methods exists)

*will be computed in a reasonable time: 3,4 days

*will be accurate enough however accuracy is not my main concern and I would like to be able to control how much resources I put into it.


It would be great to have a Haskell, Python, C# etc. library which implements it. I am not using mathlab or R but if necessary I can go with R.
 A: If you're willing to have a low-rank approximation (as you would with Lanczos-type algorithms and a limited number of singular vectors), an alternative is stochastic SVD.  You get similar accuracy and computational effort to things like irlba, but a much, much simpler implementation -- which is relevant if none of the available implementations handle precisely your situation.
If you start with an $N\times M$ matrix and you want a rank $k$ approximation, you need to be able to multiply your matrix by a dense matrix with, say, $M\times (k+10)$ entries and then do an ordinary SVD on the resulting $N\times (k+10)$ matrix. As with Lanczos-type algorithms, the computation only uses the matrix for multiplication (it's a "matrix-free" algorithm), so it will take advantage of sparsity and other structure.  (My use case was in genetics and the matrix was the product of a sparse matrix and a projection orthogonal to a low-rank subspace)
A: If it fits into memory, construct a sparse matrix in R using the Matrix package, and try irlba for the SVD.  You can specify how many singular vectors you want in the result, which is another way to limit the computation.
That's a pretty big matrix, but I've had very good results with this method in the past.  irlba is pretty state-of-the-art.  It uses the implicitly restarted Lanczos bi-diagonalization algorithm.
It can chew through the netflix prize dataset (480,189 rows by 17,770 columns, 100,480,507 non-zero entries) in milliseconds.  You dataset is ~ 200,000 times bigger than the Netflix dataset, so it take significantly longer than that.  It might be reasonable to expect that it could do the computation in a couple of days.
