CODEWITHSUNDEEP
c++sparse-matrixsvd
physicsGuy
physicsGuy

Reputation: 3675

library for full SVD of sparse matrices

I want to do a singular value decomposition for large matrices containing a lot of zeros. In particular I need U and S, obtained from the diagonalization of a symmetric matrix A. This means that A = U * S * transpose(U^*), where S is a diagonal matrix and U contains all eigenvectors as columns.

I searched the web for c++ librarys that combine SVD and sparse matrices, but could only find libraries that find a few, but not all eigenvectors. Does anyone know if there is such a library?

Also after obtaining U and S I need to multiply them to some dense vector.

Upvotes: 0

Views: 1492

Answers (1)

BrunoLevy
BrunoLevy

Reputation: 2633

For this problem, I am using a combination of different techniques:

  • Arpack can compute a set of eigenvalues and associated eigenvectors, unfortunately it is fast only for high frequencies and slow for low frequencies

  • but since the eigenvectors of the inverse of a matrix are the same as the eigenvectors of a matrix, one can factor the matrix (using a sparse matrix factorization routine, such as SuperLU, or Choldmod if the matrix is symmetric). The "communication protocol" with Arpack only expects you to compute a matrix-vector product, so if you do a linear system solve using the factored matrix instead, then this makes Arpack fast for the low frequencies of the spectrum (do not forget then to replace the eigenvalue lambda by 1/lambda !)

  • This trick can be used to explore the entire spectrum, with a generalized transform (the transform in the previous point is refered as "invert" transform). There is also a "shift-invert" transform that allows one to explore an arbitrary portion of the spectrum and have fast convergence of Arpack. Then you compute (1/lambda + sigma) instead of lambda, when sigma is a "shift" (the transform is slightly more complicated than the "invert" transform, see the references below for a full explanation).

ARPACK: http://www.caam.rice.edu/software/ARPACK/

SUPERLU: http://crd-legacy.lbl.gov/~xiaoye/SuperLU/

The numerical algorithm is explained in my article that can be downloaded here: http://alice.loria.fr/index.php/publications.html?redirect=0&Paper=ManifoldHarmonics@2008

Sourcecode is available there: https://gforge.inria.fr/frs/download.php/file/27277/manifold_harmonics-a4-src.tar.gz

See also my answer to this question: https://scicomp.stackexchange.com/questions/20243/sparse-generalized-eigensolver-using-opencl/20255#20255

Upvotes: 1

Related Questions