What fast algorithms exist for computing truncated SVD?

Question

Possibly off topic here, but there exist several (one, two) related questions already.

Poking around in the literature (or a google search for Truncated SVD Algorithms) turns up a lot of papers that use truncated SVDs in various ways, and claim (frustratingly, often without citation) that there are fast algorithms for computing it, but no one seems to be pointing at what those algorithms are.

The only thing I can find is a single randomized algorithm, used in the redSVD library.

What I'd like to see is a set of exact and inexact algorithms, suitable for understanding how the systems work (but not necessarily for actually implementing them of course!).

Does anyone have a good reference for this sort of thing?

Answer 1

Very broadly speaking, there are two approaches to compute eigenvalue or singular value decompositions. One approach is to diagonalize the matrix and this essentially yields the whole eigenvalue / singular value decomposition (the whole eigenvalue spectrum) at the same time, see some overview here: What are efficient algorithms to compute singular value decomposition (SVD)? The alternative is to use an iterative algorithm that yields one (or several) eigenvectors at a time. Iterations can be stopped after the desired number of eigenvectors has been computed.

I don't think there are iterative algorithms specifically for SVD. This is because one can compute SVD of a $n\times m$ matrix $\mathbf B$ by doing an eigendecomposition of a square symmetric $(n+m)\times(n+m)$ matrix $$\mathbf A=\left(\begin{array}{cc}0 & \mathbf B\\\mathbf B^\top & 0\end{array}\right).$$ Therefore instead of asking what algorithms compute truncated SVD, you should be asking what iterative algorithms compute eigendecomposition: $$\text{algorithm for truncated SVD} \approx \text{iterative algorithm for eigendecomposition}.$$

The simplest iterative algorithm is called power iteration and is indeed very simple:

1. Initialize random $\mathbf x$.
2. Update $\mathbf x \gets \mathbf A\mathbf x$.
3. Normalize $\mathbf x \gets \mathbf x / \|\mathbf x\|$.
4. Goto step #2 unless converged.

All the more complex algorithms are ultimately based on the power iteration idea, but do get quite sophisticated. Necessary math is given by Krylov subspaces. The algorithms are Arnoldi iteration (for square nonsymmetric matrices), Lanczos iteration (for square symmetric matrices), and variations thereof such as e.g. "implicitly restarted Lanczos method" and whatnot.

You can find this described in e.g. the following textbooks:

1. Golub & Van Loan, Matrix Computations
2. Trefethen & Bau, Numerical Linear Algebra
3. Demmel, Applied Numerical Linear Algebra
4. Saad, Numerical Methods for Large Eigenvalue Problems

All reasonable programming languages and statistic packages (Matlab, R, Python numpy, you name it) use the same Fortran libraries to perform eigen/singular-value decompositions. These are LAPACK and ARPACK. ARPACK stands for ARnoldi PACKage, and it's all about Arnoldi/Lanczos iterations. E.g. in Matlab there are two functions for SVD: svd performs full decomposition via LAPACK, and svds computes a given number of singular vectors via ARPACK and it is actually just a wrapper for an eigs call on the "square-ized" matrix.

Update

Turns out there are variants of Lanczos algorithm that are specifically tailored to perform SVD of a rectangular matrix $\mathbf B$ without explicitly constructing a square matrix $\mathbf A$ first. The central term here is Lanczos bidiagonalization; as far as I understand, it is essentially a trick to perform all the steps of Lanczos iterations on $\mathbf A$ directly on $\mathbf B$ without ever constructing $\mathbf A$ and thus saving space and time.

There is a Fortran library for these methods too, it's called PROPACK:

The software package PROPACK contains a set of functions for computing the singular value decomposition of large and sparse or structured matrices. The SVD routines are based on the Lanczos bidiagonalization algorithm with partial reorthogonalization (BPRO).

However, PROPACK seems to be much less standard than ARPACK and is not natively supported in standard programming languages. It is written by Rasmus Larsen who has a large 90-page long 1998 paper Lanczos bidiagonalization with partial reorthogonalization with what seems a good overview. Thanks to @MichaelGrant via this Computational Science SE thread.

Among the more recent papers, the most popular seems to be Baglama & Reichel, 2005, Augmented implicitly restarted Lanczos bidiagonalization methods, which is probably around the state of the art. Thanks to @Dougal for giving this link in the comments.

Update 2

There is indeed an entirely different approach described in detail in the overview paper that you cited yourself: Halko et al. 2009, Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. I don't know enough about it to comment.

Answer 2

I just stumbled on the thread via googling fast SVDs, so I'm trying to figure out things myself, but maybe you should look into adaptive cross approximation (ACA).

I don't really know what problem is like or what you need, but if your matrix $M$ is calculated from smooth functions, and you just need an approximate separated representation $M=\sum_{i=0}^k U_i\cdot V^T_i$ and not really a "proper" SVD, ACA algorithm has (almost) linear computational complexity ($N\times N$ matrix then it is almost $O(N)$). So it's really fast; unfortunately many people use the word "fast" lightly.

Again, it depends on your problem if that works. In many cases I personally encounter, the ACA is a very useful numerical tool.

Note: I wanted to write this as a comment, but because I just created this account I don't have enough reputation for comments... But posting works.

Answer 3

Here's a technique I have used successfully in the past for computing a truncated SVD (on the Netflix dataset). It is taken from this paper. In a collaborative filtering setting, I should note that most of the values are missing and the point is to predict them, so to use truncated SVD to solve such a problem, you have to use a technique that works under that condition. A short description:

1. Before you do anything, fit a simple model (e.g., global mean + column and row constant values), and only once you have done that should you move on to using truncated SVD to fit the residuals.
2. Initialize a random vector of length k (where that's the rank you're truncating to) to each row and column (to each movie and user in the Netflix case).
3. Hold the row vectors fixed and update the column vectors to minimize error w.r.t. the known entries in the matrix. The procedure is given in matlab code in the paper.
4. Hold the column vectors fixed and update the row vectors in an analogous way.
5. Repeat 3 & 4 until you converge or are getting good enough results.