When eigenvalues are not enough
By now you can chase eigenvalues with the power method, sharpen them with inverse iteration, and grind a whole spectrum out with the QR algorithm. But eigenvalues have two stubborn limitations. First, they only exist for square matrices — ask for the eigenvalues of a 1000-by-3 data matrix and the question is meaningless. Second, even for a square matrix the eigenvectors need not be orthogonal, and a defective matrix may not have a full set of them at all. The geometry can be slanted and tangled, exactly the kind of thing that wrecks numerical confidence.
The singular value decomposition (SVD) sidesteps all of this. It says that every matrix A, of any shape m-by-n, factors as A = U S V^T, where U and V are orthogonal (their columns are perpendicular unit vectors) and S is rectangular-diagonal with non-negative entries down the diagonal. No defectiveness, no complex numbers, no slanted axes — just two honest rotations and a pure stretch sandwiched between them.
Where the singular values come from
Here is the bridge back to eigenvalues. Form the square symmetric matrix A^T A. It is symmetric and positive semidefinite, so it has real, non-negative eigenvalues and a full orthonormal set of eigenvectors — that is exactly the well-behaved symmetric eigenproblem. Its eigenvectors are the columns of V, and its eigenvalues are the squares of the singular values: if A^T A v = lambda v then sigma = sqrt(lambda). The singular values are never negative because they are square roots of non-negative numbers.
A v_i = sigma_i u_i (input axis -> stretched output axis)
sigma_1 >= sigma_2 >= ... >= sigma_r > 0 = sigma_{r+1} = ...
rank(A) = r = number of nonzero singular values
A = U S V^T = sum_i sigma_i (u_i v_i^T)Notice the last line: A is the sum of rank-1 layers sigma_i times (u_i v_i^T), ordered from the strongest stretch to the weakest. That ordering is the whole reason the SVD is so useful in practice — it ranks the directions of A by importance, which is precisely what the next guide on low-rank approximation will exploit. The largest sigma is the matrix's 2-norm; the ratio sigma_1 / sigma_n is its condition number for inversion. The SVD literally hands you the conditioning of the problem.
How it is actually computed
The standard algorithm is a clever two-phase echo of how the QR algorithm tames a general matrix. Phase one squeezes A into a simpler shape exactly, with a fixed number of orthogonal moves; phase two iterates on that simple shape until the singular values fall out. Because every move is an orthogonal transformation, lengths and angles are preserved, which is the secret to the SVD's famous numerical reliability.
- Bidiagonalize. Using Householder reflections from the left and the right, reduce A to an upper-bidiagonal matrix B (nonzeros only on the diagonal and the one above it). This is Golub-Kahan bidiagonalization, and it finishes in a fixed O(m n^2) flops — no iteration, just clean orthogonal sweeps that do not change the singular values.
- Iterate to diagonal. Run an implicit, shifted QR-style sweep on the bidiagonal B, working with B implicitly so A^T A is never formed. Each sweep shaves the off-diagonal toward zero; with good shifts it converges cubically, so a handful of sweeps per singular value usually suffices.
- Read off and clean up. The converged diagonal entries are the singular values; take their absolute values and sort them in decreasing order. Accumulate the Householder and rotation moves into U and V so you recover the full A = U S V^T, then the smallest singular values that are below the round-off floor are honestly treated as numerically zero.
This is what the routine behind a one-line `svd(A)` call really does, and it is backward stable: the computed factors are the exact SVD of a matrix very close to A, differing by about machine epsilon times the norm of A. That is the realistic gold standard — remember from the conditioning rung that accuracy = conditioning times stability, so even this beautifully stable algorithm cannot recover singular values that an ill-conditioned A has already buried below the noise.
Why everyone reaches for it
The SVD is the Swiss-army knife of numerical linear algebra because the ordered singular values answer so many questions at once. The numerical rank is just the count of singular values above a sensible threshold — far more trustworthy than chasing exact zeros, since on a computer nothing is ever exactly zero. The condition number is sigma_1 / sigma_r, computed honestly without inverting anything. And when a least-squares problem is rank-deficient, the SVD builds the pseudoinverse by simply inverting the nonzero singular values and leaving the tiny, untrustworthy ones alone.
A tiny worked picture: take the 2-by-2 matrix that stretches the x-axis by 3 and the y-axis by 1, then rotates the result by 45 degrees. Its eigenvalues look complicated, but its singular values are exactly 3 and 1 — the SVD sees through the rotation straight to the true stretch factors. That is the everyday miracle: the SVD reports the genuine geometry of a map, undistorted by whatever spin happens to be wrapped around it.
There is a cost to be honest about. The full SVD takes O(m n^2) flops and stores dense U and V, which is fine for a 1000-by-50 matrix but hopeless for a sparse matrix with a million rows. For those, you do not want all the singular values, only the largest few — and that is a job for the Krylov eigensolvers from the previous rung, run on A^T A or its bidiagonal cousin, never forming the giant matrix at all. The dense SVD and the sparse iterative SVD are different tools for different sizes; reaching for the wrong one is the most common SVD mistake.
What to carry forward
Hold on to one image and one warning. The image: any matrix is rotate-stretch-rotate, and the singular values are the stretch factors, ordered from most to least important. The warning: compute the SVD through bidiagonalization plus implicit iteration, never by forming A^T A, because squaring throws away half your digits. With that ordered list of stretch factors in hand, you are perfectly set up for the finale — where keeping only the top few of them gives the best possible low-rank approximation, the mathematics behind PCA, image compression, and recommender systems.
One last honesty check, the same one that has shadowed every guide on this ladder: the SVD is computed in floating point, so its singular values and vectors are approximate, accurate to about machine epsilon times sigma_1. A backward-stable routine guarantees you the exact SVD of a nearby matrix, but it cannot resolve two singular values that sit closer together than the round-off floor, and it cannot un-bury directions that an ill-conditioned A has already drowned. The SVD is the most trustworthy tool in the toolbox — it is just not magic, and knowing exactly where its honesty ends is what separates using it from understanding it.