Every matrix in three acts

The spectral theorem is a gorgeous result with an annoying asterisk: it only speaks to symmetric matrices. Feed it a lopsided one — a data matrix that's tall and thin, a map that shears the plane, anything not equal to its own transpose — and it shrugs. Eigenvalues might go complex, eigenvectors might not be perpendicular, and a matrix that isn't even square has no eigen-decomposition to offer at all. For most of the nineteenth century that was where the story stopped. Then, in 1873, Eugenio Beltrami found a way around the asterisk, and Camille Jordan found the same thing independently the next year: every matrix, no exceptions, can be factored into a rotation, a diagonal stretch, and another rotation. The trick was to stop demanding that the matrix diagonalize itself and instead borrow the spectral theorem for a symmetric matrix built out of it. That factorization is the singular value decomposition, and it is arguably the single most useful thing in all of applied linear algebra.

Its promise is a sentence: whatever a matrix does to space, it does it in three acts. Rotate. Stretch each axis by its own factor. Rotate again. No shear left over, no mixing — once you're in the right coordinates going in and the right coordinates coming out, the middle is nothing but independent scaling. This essay is about why that's true, what the three acts look like, and why the numbers in the middle turn out to run half of machine learning.

The claim: rotate, stretch, rotate

Write it down and it's compact. Any real m \times n matrix A factors as

\boxed{\;A = U \Sigma V^\top\;}

where U (m \times m) and V (n \times n) are orthogonal — their columns are orthonormal, so they're pure rotations (or reflections), and their inverses are just their transposes — and \Sigma (m \times n) is diagonal with nonnegative entries \sigma_1 \ge \sigma_2 \ge \cdots \ge 0 down its main diagonal. Those diagonal entries are the singular values, and their descending order is a convention worth its weight: it means the first act of stretching is always the strongest.

Read the factorization right to left, the way it acts on a vector x. First V^\top rewrites x in a new orthonormal basis — a rotation into the map's input frame, the same move diagonalization made when it changed into a matrix's own coordinates. Then \Sigma scales each of those coordinates independently by its singular value, no cross-talk, because a diagonal matrix can't mix axes. Then U rotates the stretched result into the map's output frame. Three clean steps, and the awkward truth about the original A — that it sheared and skewed and turned all at once — dissolves into a rotation, a set of honest stretches, and a rotation. The columns of V are called the right singular vectors (the input axes), the columns of U the left singular vectors (the output axes).

Notice what this buys over the spectral theorem's A = Q\Lambda Q^\top. There, the same orthogonal Q appeared on both sides, which is exactly the demand a non-symmetric matrix can't meet. The SVD relaxes it: two different rotations, V going in and U coming out. That single degree of freedom is what lets the factorization exist for every matrix, including the rectangular ones that have no eigenvalues to speak of.

The circle becomes an ellipse

The cleanest way to see the three acts is to watch what A does to the unit circle — every vector of length one. A rotation (V^\top) sends the unit circle to itself; a circle has no favorite direction, so nothing visible happens yet, except that the two special input axes (the columns of V) get lined up with the coordinate axes. Then the diagonal stretch \Sigma pulls the circle out along those axes — by \sigma_1 one way, by \sigma_2 the other — and a circle stretched unequally along perpendicular axes is precisely an ellipse, sitting axis-aligned. Finally U rotates that ellipse into its final orientation.

So the whole content of a 2\times2 matrix, geometrically, is this: it maps the unit circle to an ellipse. The singular values \sigma_1, \sigma_2 are the lengths of the ellipse's semi-axes — how far the map can stretch a unit vector (\sigma_1, the longest reach) and how little it can (\sigma_2, the shortest). The columns of U point along those axes in the output. The columns of V are the two input directions that get mapped onto those axes: A v_i = \sigma_i u_i, the right singular vector going in, the left one coming out, scaled by the singular value between them. Everything the matrix can do to a length is bracketed between \sigma_1 and \sigma_2 — a fact that, generalized, is the whole theory of matrix norms and conditioning.

Where the singular values come from

Here's the sleight of hand Beltrami and Jordan pulled, and it leans entirely on the two essays before this one. You can't diagonalize a general A, but you can always form A^\top A, and that matrix is symmetric — so the spectral theorem applies to it without complaint. Better still, it's positive semidefinite: for any x, the number x^\top A^\top A\, x = \lVert Ax \rVert^2 \ge 0 is a squared length, which can't be negative, so all its eigenvalues are real and \ge 0. That's exactly the guaranteed-positive structure we need. Substitute the factorization and watch U cancel itself away:

A^\top A = (U\Sigma V^\top)^\top (U \Sigma V^\top) = V \Sigma^\top U^\top U \Sigma V^\top = V\, (\Sigma^\top\Sigma)\, V^\top,

using U^\top U = I. Read the right-hand side: it's the spectral decomposition of the symmetric matrix A^\top A. So V is its orthonormal eigenvectors, and \Sigma^\top\Sigma — a diagonal matrix of the \sigma_i^2 — holds its eigenvalues. In one line:

\sigma_i = \sqrt{\lambda_i(A^\top A)}.

The singular values are the square roots of the eigenvalues of A^\top A. That's the whole secret. The nonnegativity we needed is free, because you can't take a real square root of a negative number and the eigenvalues of A^\top A are never negative — the positive-definite essay's guarantee, cashed in exactly when we need it. By the mirror-image calculation, A A^\top = U(\Sigma\Sigma^\top)U^\top, so the columns of U are the eigenvectors of A A^\top, sharing the same nonzero eigenvalues. And once you have V and the nonzero \sigma_i, the left singular vectors come for free from the defining relation itself: u_i = A v_i / \sigma_i. There's no sign ambiguity to agonize over — you define u_i that way, and A = U\Sigma V^\top then holds on the nose.

Rank, counted honestly

How many of the singular values are actually nonzero? Call that number r. Since \Sigma does the only stretching in the sandwich and rotations never collapse anything, a direction survives the map with positive length exactly when its singular value is positive. So

\operatorname{rank}(A) = \#\{\, i : \sigma_i > 0 \,\},

the count of nonzero singular values. A singular value of zero is a direction the map annihilates — an input axis that \Sigma multiplies by nothing, flattening it into the null space. This is a far sturdier way to measure rank than counting pivots, because it comes with a dial: in floating-point arithmetic a "zero" singular value is really a tiny one, and how tiny you're willing to tolerate before calling a direction dead is the numerical rank. The SVD doesn't just tell you the rank; it ranks the directions by how much they matter, from \sigma_1 down.

Eckart–Young: the best low-rank approximation

That ranking is where the SVD stops being elegant and starts being indispensable. Break the factorization into a sum, one term per singular value:

A = \sum_{i=1}^{r} \sigma_i\, u_i v_i^\top.

Each term \sigma_i u_i v_i^\top is a rank-1 matrix — an outer product, the crudest possible building block — weighted by its singular value. The sum rebuilds A exactly, but the terms are sorted by importance, biggest \sigma first. So what if you keep only the top k and throw the rest away?

A_k = \sum_{i=1}^{k} \sigma_i\, u_i v_i^\top.

In 1936 Carl Eckart and Gale Young proved the thing that makes this matter: A_k is the best possible rank-k approximation of A. Not a good one — the optimal one, closer to A (in the least-squares/Frobenius sense, and the operator-norm sense too) than any other matrix of rank k could ever be. And the error you eat by truncating is exactly the singular values you dropped: \lVert A - A_k \rVert^2 = \sigma_{k+1}^2 + \cdots + \sigma_r^2. If the singular values decay fast — and for real data they almost always do — then a handful of terms reconstructs almost everything, and you've compressed a full matrix into a few rank-1 pieces. That one theorem is the engine under image compression, noise filtering, latent-semantic indexing, recommender systems, and the whole idea that a big matrix of data usually wants to be low-rank. Keep the loud directions, drop the quiet ones, and you've thrown away nothing that mattered.

The pseudo-inverse: a fair inverse for everything

The SVD also hands you an inverse for matrices that have no business having one. If A were square and invertible, A^{-1} = V\Sigma^{-1}U^\top — just invert each piece and reverse the order, and inverting the diagonal \Sigma means reciprocating each \sigma_i. But if A is rectangular, or square but rank-deficient, some \sigma_i are zero and you can't divide by them. The Moore–Penrose pseudo-inverse does the only sensible thing: invert the singular values you can, and leave the zero ones zero.

A^+ = V \Sigma^+ U^\top, \qquad \Sigma^+ = \operatorname{diag}\!\left(\tfrac{1}{\sigma_1}, \ldots, \tfrac{1}{\sigma_r}, 0, \ldots, 0\right).

Reciprocate every nonzero singular value; send every zero one to zero (not infinity — that's the whole discipline of it). The result A^+ is the closest thing to an inverse a general matrix can have, and it's exactly what solves the least-squares problem: the minimum-norm best-fit solution to Ax = b is x = A^+ b. The normal equations were one route to that answer; the pseudo-inverse is the same answer, spelled out in singular values, and it stays numerically sane precisely where the normal equations fall apart — when the columns are nearly dependent and some \sigma_i drifts toward zero.

Where this goes

The singular value decomposition is the spectral theorem freed of its asterisk: every matrix, symmetric or not, square or not, is a rotation, a set of independent stretches, and a rotation. Hold onto the three-act picture and a startling amount of applied mathematics turns out to be footnotes to it — best low-rank approximation, matrix norms and conditioning, the pseudo-inverse, the numerical rank. The most direct sequel is principal component analysis: point the SVD at a centered data matrix and its right singular vectors are the principal directions, its singular values the standard deviations along them, and Eckart–Young becomes the promise that the top few components capture the most variance possible. PCA is the SVD wearing a statistician's hat.

But the fastest way to believe that a matrix is three rotations-and-a-stretch is to watch one happen. The SVD playground lets you set the four entries of a 2\times2 matrix and plays out the three acts on the unit circle: V^\top turning the marked axes, \Sigma pulling the circle open into an axis-aligned ellipse, U swinging it into its final pose — and it reads off \sigma_1, \sigma_2 as the semi-axes it just drew. Flip the rank-1 switch and watch \sigma_2 clamp to zero: the ellipse collapses to a bare line segment, the best rank-1 shadow of your matrix, Eckart–Young made visible. It's one thing to read that every matrix factors into a circle and an ellipse. It's another to drag the circle open and watch the ellipse fall out.