The Singular Value Decomposition

Diagonalization asks whether a map is a pure scaling in one basis. For a rectangular matrix the question does not even make sense, input and output live in different spaces, and for a defective square matrix the answer is no.

Weaken the question and it always has an answer. Is the map a pure scaling from one orthonormal basis to another? Yes, for every real matrix without exception: rectangular, singular, defective, any shape at all.

What comes back is a ranked list. The singular values say how much the map stretches along each of a set of orthogonal directions, in decreasing order, so they measure how much of the matrix each direction accounts for. Keeping the largest few and discarding the rest gives the best approximation of that rank, which is why this decomposition underlies compression, denoising and principal component analysis.

Why the construction starts with $A^{\mathsf{T}}A$. Suppose the factorisation $A=U\mathrm{\Sigma }{V}^{\mathsf{T}}$ exists. Then

using $U^{\mathsf{T}}U=I$. That is a diagonalization of $A^{\mathsf{T}}A$ with eigenvalues ${\sigma }_i^2$ and eigenvectors the columns of $V$. So the right singular vectors and the squared singular values are forced. They are not a choice, they are read off a symmetric matrix that is guaranteed to be diagonalizable.

The guarantee comes from $A^{\mathsf{T}}A$ being symmetric, which the spectral theorem handles, and positive semidefinite: $x^{\mathsf{T}}{A}^{\mathsf{T}}Ax=\Vert Ax{\Vert }^2\ge 0$, so no eigenvalue is negative and every square root is real.

Why $u_i=A{v}_i/{\sigma }_i$ works. These are orthonormal when :

which is 0 for $i\ne j$ and 1 for $i=j$. The definition divides by ${\sigma }_i$ precisely so the result has norm 1, which is why it is unavailable when ${\sigma }_i=0$, those columns of $U$ must be filled by extending to an orthonormal basis instead.

What is unique and what is not. The singular values are unique, including their multiplicities, and they are always nonnegative and conventionally listed in decreasing order.

The singular vectors are not. Replacing $(u_i,v_i)$ by $(-u_i,-v_i)$ leaves ${\sigma }_i{u}_i{v}_i^{\mathsf{T}}$ unchanged, so sign flips are always available, but they must be made in matched pairs, since flipping only $v_i$ would change $A{v}_i$ without changing ${\sigma }_i{u}_i$. When a singular value repeats, any orthonormal basis of its subspace will serve, so there is a continuum of valid answers.

The practical consequence: an answer that differs from a reference by signs, or by a rotation within a repeated singular value's subspace, may be entirely correct. The verification $A{v}_i={\sigma }_i{u}_i$ decides the question, not agreement with a printed table.

Shapes. For $A$ of size $m\times n$: $U$ is $m\times m$, $\mathrm{\Sigma }$ is $m\times n$, $V$ is $n\times n$. The reduced form keeps only the first $r=\mathrm{rank}A$ columns of $U$ and $V$ and the $r\times r$ block of $\mathrm{\Sigma }$, which is what the rank-one expansion uses and what an implementation returns by default.

Theorem (existence). Every real $m\times n$ matrix has a singular value decomposition.

Proof. $A^{\mathsf{T}}A$ is symmetric, so by the spectral theorem it has an orthonormal basis of eigenvectors $v_1,\dots ,v_n$ with real eigenvalues ${\lambda }_i$. Each is nonnegative, since ${\lambda }_i={\lambda }_i\Vert v_i{\Vert }^2=v_i^{\mathsf{T}}{A}^{\mathsf{T}}A{v}_i=\Vert A{v}_i{\Vert }^2\ge 0$. Order them decreasingly and set ${\sigma }_i=\sqrt{{\lambda }_i}$, with $r$ of them nonzero.

For $i\le r$ define $u_i=A{v}_i/{\sigma }_i$; the calculation in the definition block shows these are orthonormal. Extend $\{u_1,\dots ,u_r\}$ to an orthonormal basis of ${\mathbb{R}}^m$, which is possible because any orthonormal set extends.

Now check $A=U\mathrm{\Sigma }{V}^{\mathsf{T}}$ by applying both sides to each $v_j$. For $j\le r$: the right side gives $U\mathrm{\Sigma }{e}_j={\sigma }_j{u}_j$, and the left gives $A{v}_j={\sigma }_j{u}_j$ by definition. For : $\Vert A{v}_j{\Vert }^2={\lambda }_j=0$, so $A{v}_j=0$, and the right side is also zero. Two linear maps agreeing on a basis are equal. $◼$

Nothing in the argument assumed $A$ square, invertible, or diagonalizable, which is why the decomposition exists for every real matrix.

Theorem (Eckart–Young). Let $A=\sum _{i=1}^r{\sigma }_i{u}_i{v}_i^{\mathsf{T}}$ and let $A_k$ be the sum truncated after $k$ terms. Then among all matrices $B$ with $\mathrm{rank}B\le k$,

is the minimum of $\Vert A-B{\Vert }_F$, and likewise in the spectral norm where the minimum is ${\sigma }_{k+1}$.

Why the error takes that form. The rank-one pieces are mutually orthogonal in the Frobenius inner product $⟨X,Y⟩=\mathrm{tr}(X^{\mathsf{T}}Y)$, because the $u_i$ are orthonormal and so are the $v_i$. The discarded terms therefore contribute their squared norms additively, and $\Vert {\sigma }_i{u}_i{v}_i^{\mathsf{T}}{\Vert }_F={\sigma }_i$. So truncation is a projection, and the error is the length of what was projected away. The same Pythagorean statement as in the orthogonality unit, one level up.

Corollary (rank). $\mathrm{rank}A$ is the number of nonzero singular values, since $A_r=A$ and no smaller truncation reproduces it.

Corollary (invariants). $\Vert A{\Vert }_F^2=\sum _i{\sigma }_i^2$, and for square $A$, $|detA|=\prod _i{\sigma }_i$. The second because $det(U\mathrm{\Sigma }{V}^{\mathsf{T}})=detUdet\mathrm{\Sigma }det{V}^{\mathsf{T}}$ with the orthogonal determinants equal to $\pm 1$. Both serve as cheap checks on a computed decomposition.

1. Form $A^{\mathsf{T}}A$. It is $n\times n$, symmetric, and positive semidefinite.
2. Find its eigenvalues ${\lambda }_1\ge \cdots \ge {\lambda }_n\ge 0$ and orthonormal eigenvectors $v_1,\dots ,v_n$. These are the columns of $V$.
3. Take ${\sigma }_i=\sqrt{{\lambda }_i}$, in the same order. Let $r$ be the count of nonzero ones.
4. For $i\le r$, set $u_i=A{v}_i/{\sigma }_i$. These are automatically orthonormal.
5. If , extend $\{u_1,\dots ,u_r\}$ to an orthonormal basis of ${\mathbb{R}}^m$ for the remaining columns of $U$.
6. Assemble $U$, $\mathrm{\Sigma }$ and $V$ with matching index order.
7. Verify $A{v}_i={\sigma }_i{u}_i$ for each $i\le r$, which is cheaper than multiplying out $U\mathrm{\Sigma }{V}^{\mathsf{T}}$ and catches the same errors.

Choosing $A^{\mathsf{T}}A$ or $A{A}^{\mathsf{T}}$. Both give the same nonzero eigenvalues. Use whichever is smaller: for a $100\times 3$ matrix, $A^{\mathsf{T}}A$ is $3\times 3$ and $A{A}^{\mathsf{T}}$ is $100\times 100$. Starting from $A{A}^{\mathsf{T}}$ gives the $u_i$ directly, and then $v_i=A^{\mathsf{T}}{u}_i/{\sigma }_i$.

Checks worth running.

The first two catch most arithmetic slips before any singular vector is computed, and they cost far less than discovering the error at step 7.

The rank-one expansion. Once assembled,

where each $u_i{v}_i^{\mathsf{T}}$ is an outer product. A full matrix of rank 1, not a scalar. Truncating after $k$ terms gives the best rank-$k$ approximation, with Frobenius error $\sqrt{{\sigma }_{k+1}^2+\cdots +{\sigma }_r^2}$.

Why this route is for hand computation only. Forming $A^{\mathsf{T}}A$ squares the condition number, so small singular values are computed with badly degraded accuracy. Numerical libraries bidiagonalise $A$ directly and never form the product. The mathematics is identical; the arithmetic is not.

Step 1: form $A^{\mathsf{T}}A$.

Symmetric, as it must be.

Step 2: its eigenvalues. Trace 50, determinant $625-400=225$, so $p(\lambda )={\lambda }^2-50\lambda +225$. The discriminant is $2500-900=1600={40}^2$, giving

Step 3: the singular values.

Check before going further. ${\sigma }_1^2+{\sigma }_2^2=45+5=50=\mathrm{tr}(A^{\mathsf{T}}A)=9+0+16+25$. And ${\sigma }_1{\sigma }_2=\sqrt{225}=15=|detA|=|15-0|$. Both hold, so the singular values are right whatever happens next.

Step 4: right singular vectors. For $\lambda =45$: , giving $v_1\propto (1,1)$. For $\lambda =5$: , giving $v_2\propto (1,-1)$. Normalising,

orthonormal since $⟨v_1,v_2⟩=\frac{1}{2}(1-1)=0$.

Step 5: left singular vectors.

Orthonormality. $⟨u_1,u_2⟩=\frac{1}{10}(3-3)=0$, and each has norm $\frac{1}{\sqrt{10}}\sqrt{1+9}=1$.

Step 6: assemble.

Step 7: verify. $A{v}_1=\frac{1}{\sqrt{2}}(3,9{)}^{\mathsf{T}}$ and ${\sigma }_1{u}_1=3\sqrt{5}\cdot \frac{1}{\sqrt{10}}(1,3{)}^{\mathsf{T}}=\frac{3}{\sqrt{2}}(1,3{)}^{\mathsf{T}}=\frac{1}{\sqrt{2}}(3,9{)}^{\mathsf{T}}$. Similarly $A{v}_2=\frac{1}{\sqrt{2}}(3,-1{)}^{\mathsf{T}}={\sigma }_2{u}_2$. Multiplying out $U\mathrm{\Sigma }{V}^{\mathsf{T}}$ returns exactly.

The rank-one expansion.

That is $A_1$, the best rank-1 approximation of $A$. The error is

exactly as Eckart–Young predicts: discarding one term costs the singular value discarded, and no other rank-1 matrix does better.

A note on signs. Choosing $v_1\propto -(1,1)$ instead would give $u_1$ with both signs flipped, and the products ${\sigma }_1{u}_1{v}_1^{\mathsf{T}}$ would be identical. An answer differing from this one by matched sign flips is correct; what is not correct is flipping one of a pair without the other.