Computing the Singular Value Decomposition for Sparse Matrices

The singular value decomposition is a powerful tool used in computational science and matrix computations. It is the most robust of all factorization methods, can be applied to "any" matrix, and reveals many key properties of a matrix. Applications include principal component analysis, signal processing, and statistics. The SVD can accurately reveal the rank, range, null space, and other key properties of a matrix. It can also be used to solve linear systems and least squares problems. That is, it can do anything that an LU, Cholesky, or QR factorization can do, and more. If the SVD can do so much, why use other methods at all? The reason is cost. It is much more expensive to compute than the other factorization methods, particularly for sparse matrices. In particular, prior to the work presented in this talk, there was no good way to compute all the singular values of a sparse matrix. Thus, in practice, "any" matrix means "any matrix you can store as a full or banded matrix."

Computing the SVD for a sparse matrix starts with a sparse QR factorization. Next, selective Givens rotations are applied to the skyline form of R, with intriguing combinatorial aspects. Let c(j) be the row index of the topmost nonzero entry in column j. Define a "corner" as a column j so that c(j-1) >= c(j) < c(j+1). A Givens rotation to annihilate a(c(j),j) causes a single fill-in entry in the skyline form, the entry R(j,j-1). Another Givens rotation between rows j-1 and j can then chase this fill-in to the right. If j is the rightmost corner, then the resulting chase causes at most one fill-in to the right. Thus, progress is ensured and the matrix is reduced to bidiagonal form. High performance is obtained by blocking together multiple sweeps of Givens rotations. All prior methods of reducing a general R to bidiagonal form destroy the sparsity of R. It is then straightforward to compute the singular values of the bidiagonal matrix. Naturally, the left and right singular vectors are normally completely nonzero, but our method can efficiently compute all the singular values of a sparse matrix, in no more space than what is required to hold the QR factor R in skyline form. Try s=svd(sparse(A)) in MATLAB (it fails); our algorithm addresses this limitation. In the time taken by s=svds(A) to iteratively compute a small fraction of the singular values, our algorithm can compute all of them.