How is the SVD of a matrix computed in practice

Question

How does MATLAB, for instance, calculate the SVD of a given matrix? I assume the answer probably involves computing the eigenvectors and eigenvalues of A*A'. If that is the case, I would also like to know how does it compute those?

Answer 1

This is typically done using the Golub-Reinsch algorithm, and no, it doesn't involve computing eigenvalues and eigenvectors of $AA^{T}$.

See

G. H. Golub and C. Reinsch. Singular Value Decomposition and Least Squares Solutions. Numerische Mathematik 14:403-420, 1970.

This material is discussed in many textbooks on numerical linear algebra.

Answer 2

Apart from the (now classical) Golub-Reinsch paper Brian notes in his answer (I have linked to the Handbook version of the paper), as well as the (also now classical) predecessor paper of Golub-Kahan, there have been a number of important developments in computing the SVD since then. First, I have to summarize how the usual method works.

The idea in computing the SVD of a matrix is qualitatively similar to the method used for computing the eigendecomposition of a symmetric matrix (and, as noted in the OP, there is an intimate relationship between them). In particular, one proceeds in two stages: the transformation to a bidiagonal matrix, and then finding the SVD of a bidiagonal matrix. This is completely analogous to the procedure of first reducing a symmetric matrix to tridiagonal form, and then computing the eigendecomposition of the resulting tridiagonal.

For computing the SVD of a bidiagonal matrix, one particularly interesting breakthrough was the paper by Jim Demmel and Velvel Kahan, which demonstrated that one can compute even the tiny singular values of a bidiagonal matrix with good accuracy, by suitably modifying the method initially proposed in Golub-Reinsch. This was then followed by the (re?)discovery of the dqd algorithm, which is a descendant of the old quotient-difference algorithm of Rutishauser. (Beresford Parlett gives a nice discussion here.) If memory serves, this is now the preferred method used internally by LAPACK. Apart from this, it has always been possible to derive SVD versions of developments in the solution of symmetric eigenproblems; for instance, there is an SVD version of divide-and-conquer, as well as an SVD version of the old Jacobi algorithm (which may be more accurate in some circumstances).

As for bidiagonalization, one improved method was outlined in Barlow's paper, which requires a bit more work than the original procedure of Golub and Reincsh, but yields more accurate bidiagonal matrices.