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

usetruncated SVDs in various ways, andclaim(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**

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:

- Initialize random $\mathbf x$.
- Update $\mathbf x \gets \mathbf A\mathbf x$.
- Normalize $\mathbf x \gets \mathbf x / \|\mathbf x\|$.
- 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:

- Golub & Van Loan, Matrix Computations
- Trefethen & Bau, Numerical Linear Algebra
- Demmel, Applied Numerical Linear Algebra
- 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.

**Attribution***Source : Link , Question Author : John Doucette , Answer Author : Community*