The matrix could be as large as $2500\times 2500$, what is the best algorithm to do that, is there some algorithm that is easy to write a program, is there any convenient packages for that?
There is a survey of decomposition algorithms in the first table of this NIPS paper.
It lists modern algorithms (with links to known implementations), including the stochastic decomposition of Halko et al., arguably the state-of-the-art method today.
You ask for convenient programming packages but don't state your platform or language of choice. Assuming it's:
- redsvd (in-core) very clean and elegant, efficient implementation
- pca.m by Mark Tygert, one of the co-authors of the stochastic method.
Your problem isn't particularly big though, so I guess any existing package (using the iterative Lanczos algorithm) will do fine, eigen decompositions have been around for a while.
Take a look at A Survey of Software for Sparse Eigenvalue Problems by Hernández et al.
I don't know much about eigenvalues or what they are applicable to, but R seems to have a built in function for this purpose named
eigen(). Calculating the eigenvalues & eigenvectors for a 2500 * 2500 matrix took ~ 1 minute on my machine.
> sampData <- runif(6250000, 0, 2) > x <- matrix(sampData, ncol = 2500, byrow = TRUE) > system.time(eigen(x)) user system elapsed 79.74 2.90 65.69
This question has also come up on Stack Overflow.
2500x2500 is not such a large problem. Even without leveraging the sparsity the SVD implementation of scipy.linalg is able to decompose it in less than a minute. See my answer to a related questions for more details.
For larger problems you will need to use the sparsity explicitly. The gensim project my help you for middle size problems that fit on a single computer but not in RAM and the mahout implementation is able to deal with sparse matrices that don't even fit on a single hard-drive.