Singular value decomposition (finite-dim)

Intuition [Beginner]

A linear map sends the unit disk in its input space to an ellipse — possibly squashed, possibly stretched, possibly rotated. The singular value decomposition is the recipe for that ellipse. It says every linear map, no matter how complicated its matrix looks, is built out of three simple ingredients in order: rotate the input, stretch along the perpendicular coordinate axes by independent factors, rotate the output.

The stretch factors are the singular values. They are non-negative numbers, ordered from largest to smallest. The largest singular value is the longest semi-axis of the output ellipse — the biggest stretch the map performs anywhere on the unit disk. The smallest non-zero singular value is the shortest semi-axis — the smallest stretch that does not collapse to a point. Zero singular values record directions the map collapses to a point.

Every matrix, square or rectangular, real or complex, admits this factorisation. The eigenvalue picture from 01.01.08 requires the matrix to be square, and even then the eigenvectors may fail to be orthogonal or fail to span the space. The singular value picture has none of those defects. It applies to every matrix, the singular vectors are always orthogonal, and they always span the input and output spaces.

Visual [Beginner]

The picture shows the unit disk in the plane, then four stages of the map $A=U\Sigma V^{\ast }$. Stage one is the original unit disk with its two coordinate axes drawn. Stage two is the rotation by $V^{\ast }$: the disk is unchanged but the axes have been turned. Stage three is the stretch by $\Sigma$: the disk has become an ellipse with semi-axes ${\sigma }_1$ and ${\sigma }_2$, the singular values. Stage four is the final rotation by $U$: the ellipse has been turned to its final orientation. The two singular values are the semi-axis lengths of the final ellipse.

The geometry is the whole story. Any linear map factors as rotate, stretch along perpendicular axes, rotate. The number of non-zero stretch factors is the rank of the map. The largest stretch factor is the biggest number by which the map can magnify a vector — the operator norm — and the ratio of the largest to the smallest non-zero stretch is the condition number, the amplification of relative error in numerical computation.

Worked example [Beginner]

Take the two-by-two matrix

$$A=\left(\begin{array}{cc}3& 0\\ 0& 2\end{array}\right).$$

This matrix already stretches the horizontal axis by $3$ and the vertical axis by $2$, with no rotation. Its singular value decomposition is the simplest possible: $A=I\cdot \Sigma \cdot I^{\ast }$ with $\Sigma =A$ itself, $I$ the identity. The singular values are ${\sigma }_1=3$ and ${\sigma }_2=2$. The right singular vectors are the columns of the input identity, $(1,0)$ and $(0,1)$; the left singular vectors are the columns of the output identity, the same two vectors.

A less obvious example. Take

$$A=\left(\begin{array}{cc}1& 2\\ 0& 1\end{array}\right).$$

To find the singular values, compute $A^{T}A$:

$$A^{T}A=\left(\begin{array}{cc}1& 0\\ 2& 1\end{array}\right)\left(\begin{array}{cc}1& 2\\ 0& 1\end{array}\right)=\left(\begin{array}{cc}1& 2\\ 2& 5\end{array}\right).$$

The eigenvalues of $A^{T}A$ are the roots of $t^2-6t+1=0$, namely $t=3\pm 2\sqrt{2}$. The singular values are the non-negative square roots: ${\sigma }_1=\sqrt{3+2\sqrt{2}}\approx 2.41$ and ${\sigma }_2=\sqrt{3-2\sqrt{2}}\approx 0.41$. The product ${\sigma }_1{\sigma }_2=\sqrt{(3+2\sqrt{2})(3-2\sqrt{2})}=\sqrt{9-8}=1$ equals the absolute determinant of $A$, as it must.

What this tells us. The matrix $A$ has eigenvalue $1$ (with algebraic multiplicity $2$) but its action on the unit disk is highly non-uniform: it stretches one direction by a factor of $2.41$ and compresses another by a factor of $0.41$. The eigenvalues hide this asymmetry. The singular values expose it directly.

Check your understanding [Beginner]

Formal definition [Intermediate+]

Let $\mathbb{F}$ denote $\mathbb{R}$ or $\mathbb{C}$, and write $A^{\ast }$ for the conjugate transpose of a matrix $A$ over $\mathbb{F}$ (so $A^{\ast }=A^T$ over $\mathbb{R}$). A singular value decomposition of an $m\times n$ matrix $A$ over $\mathbb{F}$ is a factorisation

$$A=U\Sigma V^{\ast }$$

with:

• $U\in M_m(\mathbb{F})$ unitary, i.e., $U^{\ast }U=U{U}^{\ast }=I_m$;
• $V\in M_n(\mathbb{F})$ unitary, i.e., $V^{\ast }V=V{V}^{\ast }=I_n$;
• $\Sigma \in M_{m\times n}(\mathbb{R})$ a rectangular "diagonal" matrix, meaning ${\Sigma }_{ij}=0$ for $i\ne j$, with $({\Sigma }_{ii}{)}_{i=1}^{\min (m,n)}=({\sigma }_1,{\sigma }_2,\dots ,{\sigma }_{\min (m,n)})$ a decreasing sequence of non-negative real numbers ${\sigma }_1\ge {\sigma }_2\ge \cdots \ge {\sigma }_r>0={\sigma }_{r+1}=\cdots ={\sigma }_{\min (m,n)}$, where $r=\mathrm{rank}A$.

The numbers ${\sigma }_1,\dots ,{\sigma }_{\min (m,n)}$ are the singular values of $A$; the non-zero ones are ${\sigma }_1,\dots ,{\sigma }_r$. The columns of $V$ are the right singular vectors of $A$ — an orthonormal basis of the input space ${\mathbb{F}}^n$ — and the columns of $U$ are the left singular vectors — an orthonormal basis of the output space ${\mathbb{F}}^m$. The pairing $A{v}_i={\sigma }_i{u}_i$ for $i\le r$ links the two bases ^{[Horn, R. A. & Johnson, C. R. — Matrix Analysis (2nd ed.)]}.

Relation to eigenvalues of $A^{\ast }A$ and $A{A}^{\ast }$. The Gram matrices $A^{\ast }A\in M_n(\mathbb{F})$ and $A{A}^{\ast }\in M_m(\mathbb{F})$ are Hermitian and positive semidefinite, hence their eigenvalues are non-negative real numbers. The singular values of $A$ are the non-negative square roots of the eigenvalues of $A^{\ast }A$ (equivalently of $A{A}^{\ast }$), counted with multiplicity. Concretely, if $V$ is an orthonormal eigenbasis for $A^{\ast }A$ with eigenvalues ${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_n\ge 0$, then ${\sigma }_i=\sqrt{{\lambda }_i}$, the columns of $V$ are the right singular vectors, and the left singular vectors are obtained by setting $u_i=A{v}_i/{\sigma }_i$ when ${\sigma }_i>0$.

Reduced versus full SVD. The factorisation written above is the full SVD with $U$ of size $m\times m$ and $V$ of size $n\times n$. The reduced (or thin) SVD discards the columns of $U$ and $V$ that pair with zero singular values: writing $r=\mathrm{rank}A$ and letting $U_r,V_r$ be the matrices whose columns are the first $r$ left and right singular vectors, ${\Sigma }_r=\mathrm{diag}({\sigma }_1,\dots ,{\sigma }_r)\in M_r(\mathbb{R})$ the diagonal of non-zero singular values, the reduced SVD is

$$A=U_r{\Sigma }_r{V}_r^{\ast }=\sum _{i=1}^r{\sigma }_i{u}_i{v}_i^{\ast }.$$

The dyadic sum on the right exhibits $A$ as a sum of $r$ rank-$1$ matrices, ordered by decreasing singular value. This is the form used by Eckart-Young: truncating the sum at $k<r$ produces the best rank-$k$ approximation to $A$ (master tier).

Uniqueness. The singular values ${\sigma }_1\ge {\sigma }_2\ge \cdots \ge {\sigma }_{\min (m,n)}$ are uniquely determined by $A$. The unitary matrices $U$ and $V$ are unique up to (i) replacing $u_i,v_i$ by $e^{i{\theta }_i}{u}_i,e^{i{\theta }_i}{v}_i$ for any phases ${\theta }_i\in \mathbb{R}$ on each pair with ${\sigma }_i$ a simple non-zero singular value (over $\mathbb{R}$ the phase is $\pm 1$), and (ii) replacing the columns $\{u_i,v_i{\}}_{i\in I}$ corresponding to a repeated singular value ${\sigma }_i={\sigma }_j$ for $i,j\in I$ by a unitary recombination within that block, and a corresponding unitary recombination on the kernel and cokernel (for $\sigma =0$).

Operator norm and Frobenius norm. The operator norm (or spectral norm) of $A$ is $\Vert A{\Vert }_{\mathrm{op}}={\sup }_{\Vert x\Vert =1}\Vert Ax\Vert ={\sigma }_1$, the largest singular value. The Frobenius norm is $\Vert A{\Vert }_F=({\sum }_{ij}|A_{ij}{|}^2{)}^{1/2}=({\sum }_i{\sigma }_i^2{)}^{1/2}$, the root-sum-square of the singular values. Both norms are unitarily invariant: $\Vert UA{V}^{\ast }\Vert =\Vert A\Vert$ for any unitaries $U,V$, which makes them functions of the singular values alone.

Counterexamples to common slips

• Singular values are not eigenvalues. For the Jordan block $A=\left(\begin{array}{cc}1& 2\\ 0& 1\end{array}\right)$, the only eigenvalue is $1$ with algebraic multiplicity $2$, but the singular values are $\sqrt{3+2\sqrt{2}}\approx 2.41$ and $\sqrt{3-2\sqrt{2}}\approx 0.41$. The two notions coincide for Hermitian positive semidefinite matrices and otherwise differ.

• The convention ${\sigma }_1\ge {\sigma }_2\ge \cdots$ is universal in the modern literature, but the matrix $\Sigma$ in the factorisation $A=U\Sigma V^{\ast }$ is not uniquely determined as a matrix on the nose — only its diagonal entries are (up to the ordering convention). Permutations of singular values force corresponding column permutations on $U$ and $V$.

• The full SVD has $U$ of size $m\times m$ and $V$ of size $n\times n$, both square unitary. For a tall matrix ($m>n$), the matrix $\Sigma$ is $m\times n$ with the diagonal entries in its top $n\times n$ block and zero rows below; for a wide matrix ($m<n$), the diagonal entries fill an $m\times m$ block on the left and zero columns on the right. The reduced SVD keeps only the non-zero structure.

• The eigenvalue decomposition $A=PD{P}^{-1}$ of a square diagonalisable matrix is, in general, not an SVD. The matrix $P$ need not be unitary, and $D$ may contain negative or complex entries. The SVD coincides with the eigendecomposition iff $A$ is Hermitian positive semidefinite, in which case $P$ can be chosen unitary and $D$ has non-negative real diagonal.

Key theorem with proof [Intermediate+]

Theorem (existence of the singular value decomposition). Let $A$ be an $m\times n$ matrix over $\mathbb{F}=\mathbb{R}$ or $\mathbb{C}$, and let $r=\mathrm{rank}A$. There exist unitary matrices $U\in M_m(\mathbb{F})$ and $V\in M_n(\mathbb{F})$ and an $m\times n$ rectangular diagonal matrix $\Sigma$ with real diagonal entries ${\sigma }_1\ge {\sigma }_2\ge \cdots \ge {\sigma }_r>0={\sigma }_{r+1}=\cdots ={\sigma }_{\min (m,n)}$ such that $A=U\Sigma V^{\ast }$. ^{[Horn, R. A. & Johnson, C. R. — Matrix Analysis (2nd ed.)]}

Proof. The argument has three steps: build the right singular vectors from the spectral theorem on $A^{\ast }A$, define the left singular vectors as the images $A{v}_i/{\sigma }_i$ rescaled, and complete to a full unitary by orthogonal extension.

Step 1: Right singular vectors and singular values. The matrix $A^{\ast }A\in M_n(\mathbb{F})$ is Hermitian, since $(A^{\ast }A{)}^{\ast }=A^{\ast }(A^{\ast }{)}^{\ast }=A^{\ast }A$, and positive semidefinite: for any $v\in {\mathbb{F}}^n$, $v^{\ast }(A^{\ast }A)v=(Av{)}^{\ast }(Av)=\Vert Av{\Vert }^2\ge 0$, with equality iff $Av=0$, iff $v\in \ker A$. By the finite-dimensional spectral theorem for Hermitian operators [01.01.08, master tier], $A^{\ast }A$ has an orthonormal eigenbasis with real non-negative eigenvalues. Let ${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_n\ge 0$ be these eigenvalues in decreasing order, and let $v_1,v_2,\dots ,v_n\in {\mathbb{F}}^n$ be a corresponding orthonormal eigenbasis: $A^{\ast }A{v}_i={\lambda }_i{v}_i$ and $v_i^{\ast }{v}_j={\delta }_{ij}$.

Define ${\sigma }_i=\sqrt{{\lambda }_i}$ for $i=1,\dots ,n$. The number ${\sigma }_i$ is the $i$-th singular value of $A$. Let $r$ be the largest index with ${\sigma }_r>0$, so ${\sigma }_1\ge \cdots \ge {\sigma }_r>0={\sigma }_{r+1}=\cdots ={\sigma }_n$. The rank of $A$ equals $r$ because $\ker (A^{\ast }A)=\ker A$ (the equality $v^{\ast }{A}^{\ast }Av=\Vert Av{\Vert }^2=0$ iff $Av=0$) and the rank of $A$ is $n-\dim \ker A=n-\dim \ker A^{\ast }A=$ number of non-zero eigenvalues of $A^{\ast }A=r$.

Step 2: Left singular vectors via the image map. For each $i$ with ${\sigma }_i>0$, define the $i$-th left singular vector by

$$u_i=\frac{1}{{\sigma }_i}A{v}_i\in {\mathbb{F}}^m.$$

The vectors $u_1,\dots ,u_r$ are orthonormal. To verify: for any $i,j$ with $1\le i,j\le r$,

$$u_j^{\ast }{u}_i=\frac{1}{{\sigma }_i{\sigma }_j}(A{v}_j{)}^{\ast }(A{v}_i)=\frac{1}{{\sigma }_i{\sigma }_j}{v}_j^{\ast }(A^{\ast }A)v_i=\frac{{\lambda }_i}{{\sigma }_i{\sigma }_j}{v}_j^{\ast }{v}_i=\frac{{\sigma }_i^2}{{\sigma }_i{\sigma }_j}{\delta }_{ij}={\delta }_{ij},$$

using $A^{\ast }A{v}_i={\lambda }_i{v}_i={\sigma }_i^2{v}_i$ and the orthonormality of the $v_j$. So $\{u_1,\dots ,u_r\}\subseteq {\mathbb{F}}^m$ is an orthonormal set.

Step 3: Completion and the matrix factorisation. Extend $\{u_1,\dots ,u_r\}$ to an orthonormal basis $\{u_1,\dots ,u_m\}$ of ${\mathbb{F}}^m$ by Gram-Schmidt or by completing through the orthogonal complement of the span. Form the matrices $U=[u_1\mid u_2\mid \cdots \mid u_m]\in M_m(\mathbb{F})$ and $V=[v_1\mid v_2\mid \cdots \mid v_n]\in M_n(\mathbb{F})$ with the singular vectors as columns. Both are unitary because their columns are orthonormal bases. Let $\Sigma \in M_{m\times n}(\mathbb{R})$ be the rectangular diagonal matrix with ${\Sigma }_{ii}={\sigma }_i$ for $1\le i\le \min (m,n)$ and zero off the diagonal.

Verify $A=U\Sigma V^{\ast }$ by computing both sides on the orthonormal basis $\{v_1,\dots ,v_n\}$. For $i\le r$, the right-hand side gives $U\Sigma V^{\ast }{v}_i=U\Sigma e_i=U({\sigma }_i{e}_i)={\sigma }_i{u}_i$, and the left-hand side gives $A{v}_i={\sigma }_i{u}_i$ by definition of $u_i$. For $i>r$, the right-hand side gives $U\Sigma e_i=0$ (the $i$-th diagonal entry of $\Sigma$ is zero), and the left-hand side gives $A{v}_i=0$ because $v_i\in \ker A^{\ast }A=\ker A$. The two matrices agree on a basis of ${\mathbb{F}}^n$, so they are equal. $\square$

Corollary (uniqueness of singular values). The singular values ${\sigma }_1\ge {\sigma }_2\ge \cdots$ of $A$ are uniquely determined by $A$: any two singular value decompositions of $A$ have the same diagonal entries of $\Sigma$, listed in the same decreasing order. The proof follows because the ${\sigma }_i^2$ are the eigenvalues of $A^{\ast }A$, which are uniquely determined by $A$, and the singular value convention fixes the order.

Bridge. The singular value decomposition is the universal factorisation of a linear map between inner-product spaces. The existence proof unpacks into four interlocking syntheses, each connecting the SVD to a downstream structure that uses singular values as its building block.

First synthesis: SVD as the orbit decomposition of bi-unitary action. The group $U(m)\times U(n)$ acts on $m\times n$ complex matrices by $(P,Q)\cdot A=PA{Q}^{\ast }$. The orbits of this action are classified by the singular value tuple $({\sigma }_1,{\sigma }_2,\dots )$: two matrices lie in the same orbit iff they have the same singular values. The SVD itself is the statement that the orbit through $A$ contains a unique representative of the form $\mathrm{diag}({\sigma }_1,\dots ,{\sigma }_{\min (m,n)})$. This packages the entire theorem as the orbit-classification result for the bi-unitary action, an instance of the more general invariant-theoretic perspective on linear-algebra canonical forms.

Second synthesis: SVD and the Moore-Penrose pseudoinverse. The least-squares problem ${\min }_x\Vert Ax-b{\Vert }_2$ has a unique minimum-norm solution given by $x=A^{+}b$ with $A^{+}=V{\Sigma }^{+}{U}^{\ast }$ the Moore-Penrose pseudoinverse of $A$, where ${\Sigma }^{+}$ is the rectangular diagonal matrix with reciprocated non-zero singular values. The SVD makes the pseudoinverse explicit and reveals its geometric content: $A^{+}$ inverts the action of $A$ on its row space and annihilates the cokernel. The pseudoinverse underlies linear regression, the normal equations, and the Tikhonov-regularised least-squares used when ${\sigma }_r$ is small.

Third synthesis: SVD as the building block of principal component analysis. For a centred data matrix $X\in M_{N\times d}(\mathbb{R})$ with $N$ samples in ${\mathbb{R}}^d$, the SVD $X=U\Sigma V^T$ identifies the right singular vectors (columns of $V$) as the principal directions of the data: the directions of maximum variance, ordered by decreasing variance ${\sigma }_i^2/N$. Truncating the SVD at $k$ singular values projects the data onto the $k$-dimensional subspace capturing the maximum total variance — the principal component projection. The Eckart-Young theorem (master tier) is the optimality statement: this truncation is the best rank-$k$ approximation in Frobenius norm.

Fourth synthesis: SVD as the finite-dimensional Schmidt decomposition. For compact operators $T$ on a separable Hilbert space, the Schmidt decomposition writes $T={\sum }_n{\sigma }_n{u}_n\otimes v_n^{\ast }$ with ${\sigma }_n\ge 0$ decreasing to zero and $\{u_n\},\{v_n\}$ orthonormal sequences. Finite-rank operators recover the finite-dimensional SVD; general compact operators add the limiting condition ${\sigma }_n\to 0$. The singular values are the singular numbers of $T$ in the infinite-dimensional setting, foundational for the theory of trace-class and Hilbert-Schmidt operators and for the quantum-information notion of entanglement entropy (where the squared singular values of a bipartite density matrix are the Schmidt coefficients).

Exercises [Intermediate+]

Lean formalization [Intermediate+]

Mathlib packages the constituent ingredients of the SVD — the conjugate transpose Matrix.conjTranspose, the Hermitian and positive-semidefinite predicates Matrix.IsHermitian and Matrix.PosSemidef, the lemma isHermitian_transpose_mul_self showing $A^{\ast }A$ is Hermitian for every rectangular matrix $A$, the spectral theorem on Hermitian matrices, and the polar decomposition for square matrices over $\mathbb{C}$. The companion file Codex.Foundations.LinearAlgebra.SVD records the named statements used above.

[object Promise]

This unit is marked lean_status: partial because Mathlib supplies the spectral theorem on $A^{\ast }A$ and the constituent positivity and Hermitian-symmetry lemmas, but does not package the rectangular SVD theorem under a single name with a uniqueness statement and the Eckart-Young companion. The corresponding statement in the companion module is left as a sorry-gated alias pending that packaging.

Advanced results [Master]

Eckart-Young-Mirsky theorem. Let $A$ be an $m\times n$ matrix with SVD $A={\sum }_{i=1}^r{\sigma }_i{u}_i{v}_i^{\ast }$, and let $\Vert \cdot \Vert$ be any unitarily invariant norm on $M_{m\times n}(\mathbb{F})$ — a norm satisfying $\Vert UA{V}^{\ast }\Vert =\Vert A\Vert$ for every pair of unitaries $U,V$. The best rank-$k$ approximation to $A$ in this norm is

$$A_k=\sum _{i=1}^k{\sigma }_i{u}_i{v}_i^{\ast },\underset{\mathrm{rank}B\le k}{\min }\Vert A-B\Vert =\Vert \mathrm{diag}(0,\dots ,0,{\sigma }_{k+1},\dots ,{\sigma }_r)\Vert .$$

For the operator norm (spectral norm) $\Vert \cdot {\Vert }_{\mathrm{op}}$, the minimum is ${\sigma }_{k+1}$. For the Frobenius norm $\Vert \cdot {\Vert }_F$, the minimum is $({\sigma }_{k+1}^2+\cdots +{\sigma }_r^2{)}^{1/2}$. For the nuclear norm (trace norm) $\Vert A{\Vert }_{\ast }={\sum }_i{\sigma }_i$, the minimum is ${\sigma }_{k+1}+\cdots +{\sigma }_r$. The original Eckart-Young (1936) treated the Frobenius case; Mirsky (1960) extended the theorem to arbitrary unitarily invariant norms via von Neumann's theorem that every such norm is a symmetric gauge function applied to the singular values ^{[Eckart, C. & Young, G. — The approximation of one matrix by another of lower rank]}.

Proof sketch. Unitarily-invariant norms depend only on singular values. The matrix $A-B$ has singular values bounded below by the Cauchy interlacing inequalities for compressions: writing ${\sigma }_i(M)$ for the $i$-th singular value of $M$, the Weyl inequality gives ${\sigma }_{k+i}(A)\le {\sigma }_i(A-B)+{\sigma }_k(B)$, and if $\mathrm{rank}B\le k$ then ${\sigma }_{k+1}(B)=0$, forcing ${\sigma }_i(A-B)\ge {\sigma }_{k+i}(A)$ for $i\ge 1$. Hence the singular value sequence of $A-B$ majorises $(0,\dots ,0,{\sigma }_{k+1},{\sigma }_{k+2},\dots ,{\sigma }_r)$, and any monotone unitarily-invariant norm is minimised when this majorisation is tight, which is exactly the truncation $A_k$.

Polar decomposition and the matrix modulus. For a square $n\times n$ matrix $A$ over $\mathbb{F}$, the polar decomposition $A=QP$ has $Q=U{V}^{\ast }$ unitary and $P=V\Sigma V^{\ast }=(A^{\ast }A{)}^{1/2}$ Hermitian positive semidefinite. The factor $P$ is the modulus of $A$, denoted $|A|$ in operator theory, and equals the unique Hermitian positive semidefinite square root of $A^{\ast }A$. The polar decomposition is the matrix generalisation of the polar form $z=r{e}^{i\theta }$ of a complex number, with $r\ge 0$ replaced by $P$ and $e^{i\theta }$ replaced by $Q$. The decomposition is unique when $A$ is invertible; when $A$ is singular the factor $P$ is still unique but $Q$ is only unique on the row space of $A$ and can be chosen arbitrarily unitary on the orthogonal complement. Polar decomposition exists in infinite dimensions for bounded operators on a Hilbert space (with $Q$ a partial isometry in general) and is the starting point for the theory of polar coordinates in operator theory ^{[Horn, R. A. & Johnson, C. R. — Matrix Analysis (2nd ed.)]}.

Moore-Penrose pseudoinverse. For an $m\times n$ matrix $A$ with SVD $A=U\Sigma V^{\ast }$, the Moore-Penrose pseudoinverse is

$$A^{+}=V{\Sigma }^{+}{U}^{\ast }\in M_{n\times m}(\mathbb{F}),$$

where ${\Sigma }^{+}$ is the $n\times m$ rectangular diagonal matrix with $({\Sigma }^{+}{)}_{ii}=1/{\sigma }_i$ for ${\sigma }_i>0$ and zero otherwise. The pseudoinverse satisfies the four Moore-Penrose conditions: (i) $A{A}^{+}A=A$, (ii) $A^{+}A{A}^{+}=A^{+}$, (iii) $(A{A}^{+}{)}^{\ast }=A{A}^{+}$, (iv) $(A^{+}A{)}^{\ast }=A^{+}A$. These four conditions characterise $A^{+}$ uniquely. Application to least-squares: the equation $Ax=b$ may have no exact solution; the vector $x^{\ast }=A^{+}b$ is the unique vector minimising $\Vert Ax-b{\Vert }_2$ among all minimisers $\arg {\min }_x\Vert Ax-b{\Vert }_2$, namely the minimum-norm minimiser. The vector $A{A}^{+}b$ is the orthogonal projection of $b$ onto the column space of $A$, and $I-A^{+}A$ is the orthogonal projection onto the kernel of $A$.

Singular values and unitarily invariant norms. A symmetric gauge function on ${\mathbb{R}}^{\min (m,n)}$ is a norm $\varphi$ that is invariant under permutations and sign changes of its arguments. By von Neumann's theorem (1937), unitarily invariant norms on $M_{m\times n}(\mathbb{F})$ are in bijection with symmetric gauge functions on ${\mathbb{R}}^{\min (m,n)}$: every unitarily invariant norm $\Vert \cdot \Vert$ has the form $\Vert A\Vert =\varphi ({\sigma }_1(A),{\sigma }_2(A),\dots )$ for a unique symmetric gauge function $\varphi$. The three canonical examples — operator, Frobenius, nuclear — correspond to $\varphi (x)={\max }_i|x_i|$, $\varphi (x)=({\sum }_i{x}_i^2{)}^{1/2}$, $\varphi (x)={\sum }_i|x_i|$. This identifies the lattice of unitarily invariant matrix norms with the lattice of symmetric gauge functions, and underlies the Schatten $p$-norms $\Vert A{\Vert }_p=({\sum }_i{\sigma }_i^p{)}^{1/p}$ for $1\le p\le \infty$, interpolating between the nuclear ($p=1$), Frobenius ($p=2$), and operator ($p=\infty$) norms.

Schmidt decomposition for compact operators. A bounded linear operator $T:H\to H$ on a separable Hilbert space $H$ is compact if it maps bounded sets to relatively compact sets, equivalently (in a Hilbert space) if it is the norm limit of a sequence of finite-rank operators. For compact $T$, the operator $T^{\ast }T$ is compact, self-adjoint, positive semidefinite, with a countable orthonormal eigenbasis $\{v_n\}$ and corresponding eigenvalues ${\lambda }_n\ge 0$ accumulating only at zero. Setting ${\sigma }_n=\sqrt{{\lambda }_n}$ and $u_n=T{v}_n/{\sigma }_n$ for ${\sigma }_n>0$, the Schmidt decomposition

$$T=\sum _{n:{\sigma }_n>0}{\sigma }_n{u}_n\otimes v_n^{\ast }$$

converges in the operator norm, with ${\sigma }_n\to 0$. The numbers ${\sigma }_n$ are the singular numbers of $T$. Trace-class operators are those with ${\sum }_n{\sigma }_n<\infty$ (the trace norm is finite); Hilbert-Schmidt operators are those with ${\sum }_n{\sigma }_n^2<\infty$ (the Hilbert-Schmidt norm is finite). The Schatten $p$-class operators have ${\sum }_n{\sigma }_n^p<\infty$; these are the operator-algebraic analogues of ${\ell }^p$ sequences and form the natural framework for non-commutative integration theory ^{[Schmidt, E. — Zur Theorie der linearen und nichtlinearen Integralgleichungen, I. Teil]}.

Numerical computation. The Golub-Kahan SVD algorithm (1965) computes the SVD of an $m\times n$ matrix in $O(\min (m{n}^2,m^{2}n))$ floating-point operations. Stage one reduces $A$ to bidiagonal form by a sequence of Householder reflections from the left and the right, costing $O(mn\min (m,n))$ operations. Stage two diagonalises the bidiagonal form by a variant of the QR algorithm specialised to bidiagonal matrices, costing $O(\min (m,n{)}^2)$ operations per sweep with typically $O(\log (1/ϵ))$ sweeps for accuracy $ϵ$. The algorithm is numerically backward stable: the computed SVD is the exact SVD of a nearby matrix $A+E$ with $\Vert E{\Vert }_F=O({ϵ}_{\mathrm{mach}}\Vert A{\Vert }_F)$. Modern implementations in LAPACK and the BLAS form the computational backbone of every linear-algebra package and underlie the matrix-factorisation methods in numerical optimisation, scientific computing, and machine learning ^{[Golub, G. H. & Kahan, W. — Calculating the singular values and pseudo-inverse of a matrix]}.

Conjugacy classes of the bi-unitary action. The group $U(m)\times U(n)$ acts on $M_{m\times n}(\mathbb{C})$ by $(P,Q)\cdot A=PA{Q}^{\ast }$. The orbits of this action are parametrised by the singular value tuple $({\sigma }_1,{\sigma }_2,\dots ,{\sigma }_{\min (m,n)})$: two matrices lie in the same orbit iff they have the same singular values. The orbit space is the Weyl chamber $\{({\sigma }_1,\dots ,{\sigma }_{\min (m,n)})\in {\mathbb{R}}_{\ge 0}^{\min (m,n)}:{\sigma }_1\ge \cdots \ge {\sigma }_{\min (m,n)}\ge 0\}$. The SVD is the orbit-decomposition: every orbit contains a unique representative of diagonal form with decreasing non-negative entries. This perspective places the SVD in parallel with the Jordan canonical form 01.01.11 (orbit decomposition of ${\mathrm{GL}}_n$ on $M_n$ by conjugation) and the Cartan decomposition $G=KAK$ for a real reductive Lie group $G$ with maximal compact subgroup $K$ and maximal abelian subgroup $A$: in the case $G={\mathrm{GL}}_n(\mathbb{C})$, $K=U(n)\times U(n)$ (acting bi-unitarily), and $A$ the positive diagonal matrices, the Cartan decomposition is the SVD.

Synthesis. First synthesis — SVD as the universal generalisation of diagonalisation: every matrix admits an SVD; only Hermitian positive semidefinite matrices admit a unitary eigendecomposition coinciding with the SVD. The eigenvalue picture from 01.01.08 and the Jordan canonical form from 01.01.11 capture intrinsic features of a square operator under conjugation; the SVD captures extrinsic features of a linear map between inner-product spaces under the bi-unitary action. The two perspectives are complementary: similarity vs unitary equivalence.

Second synthesis — SVD as the geometric content of $A$: the singular values are the semi-axes of the image ellipsoid of the unit ball, and the unitary factors $U,V$ rotate that ellipsoid into standard position. This packages the operator norm (largest singular value), the condition number (${\sigma }_1/{\sigma }_r$ for invertible $A$), the rank (number of non-zero singular values), the kernel (span of right singular vectors with zero singular value), the image (span of left singular vectors with non-zero singular value), and the low-rank structure (truncated SVD) into a single coordinate-free package.

Third synthesis — SVD as the foundation of least-squares and data analysis: the Moore-Penrose pseudoinverse $A^{+}=V{\Sigma }^{+}{U}^{\ast }$ delivers the minimum-norm least-squares solution to $Ax=b$ in one formula. Principal component analysis is the SVD of the centred data matrix. Latent semantic indexing is the truncated SVD of the term-document matrix. Image compression is the truncated SVD of a grayscale pixel matrix. Recommender systems use matrix completion algorithms that fit a low-rank SVD model to sparse observed entries. In each case, the SVD identifies the dominant structure and quantifies the discarded residual.

Fourth synthesis — SVD as the spectral theorem in disguise: the SVD of $A$ encodes two paired spectral theorems, one on $A^{\ast }A$ and one on $A{A}^{\ast }$, with the operator $A$ as the intertwiner. The right singular vectors diagonalise $A^{\ast }A$, the left singular vectors diagonalise $A{A}^{\ast }$, and the singular values appear as the square roots of the shared non-zero eigenvalues. In the infinite-dimensional Hilbert-space setting, this generalises to the Schmidt decomposition for compact operators and to the polar decomposition for arbitrary bounded operators, with the positive part $|A|=(A^{\ast }A{)}^{1/2}$ playing the role of the singular value spectrum and the partial-isometry phase playing the role of $U{V}^{\ast }$. The SVD is therefore the finite-dimensional shadow of the spectral and polar decompositions of operator theory.

Full proof set [Master]

Existence of the SVD via the spectral theorem on $A^ A$.* Let $A\in M_{m\times n}(\mathbb{F})$ with $\mathbb{F}=\mathbb{R}$ or $\mathbb{C}$. The Gram matrix $A^{\ast }A\in M_n(\mathbb{F})$ is Hermitian, since $(A^{\ast }A{)}^{\ast }=A^{\ast }(A^{\ast }{)}^{\ast }=A^{\ast }A$, and positive semidefinite, since $v^{\ast }(A^{\ast }A)v=(Av{)}^{\ast }(Av)=\Vert Av{\Vert }^2\ge 0$ for every $v\in {\mathbb{F}}^n$. By the finite-dimensional spectral theorem on Hermitian operators on an inner-product space, $A^{\ast }A$ admits an orthonormal eigenbasis $\{v_1,v_2,\dots ,v_n\}\subset {\mathbb{F}}^n$ with real eigenvalues ${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_n\ge 0$. Set ${\sigma }_i=\sqrt{{\lambda }_i}$ and let $r$ be the number of non-zero ${\sigma }_i$.

For $1\le i\le r$, define $u_i=A{v}_i/{\sigma }_i\in {\mathbb{F}}^m$. The set $\{u_1,\dots ,u_r\}$ is orthonormal: $u_j^{\ast }{u}_i=(A{v}_j{)}^{\ast }(A{v}_i)/({\sigma }_j{\sigma }_i)=v_j^{\ast }(A^{\ast }A)v_i/({\sigma }_j{\sigma }_i)={\lambda }_i{\delta }_{ij}/({\sigma }_j{\sigma }_i)={\delta }_{ij}$. Extend to an orthonormal basis $\{u_1,\dots ,u_m\}$ of ${\mathbb{F}}^m$ via Gram-Schmidt on the orthogonal complement of the span of the first $r$. Let $U=[u_1\mid \cdots \mid u_m]$, $V=[v_1\mid \cdots \mid v_n]$, and $\Sigma$ the $m\times n$ rectangular diagonal matrix with ${\Sigma }_{ii}={\sigma }_i$ for $i\le \min (m,n)$. Verify on the basis $\{v_i\}$: for $i\le r$, $A{v}_i={\sigma }_i{u}_i=U\Sigma e_i\cdot$ (a column-of-$U$ identity), and for $i>r$, $A{v}_i=0$ because $v_i\in \ker (A^{\ast }A)=\ker A$. Reassembling in matrix form, $AV=U\Sigma$, equivalently $A=U\Sigma V^{\ast }$.

Uniqueness of singular values. The squared singular values ${\sigma }_i^2$ are the eigenvalues of the Hermitian matrix $A^{\ast }A$, which are uniquely determined by $A$ as the roots of the characteristic polynomial ${\chi }_{A^{\ast }A}(t)$, counted with algebraic multiplicity, and the spectral theorem confirms that algebraic multiplicity equals geometric multiplicity for Hermitian matrices. The decreasing-order convention ${\sigma }_1\ge {\sigma }_2\ge \cdots$ then fixes the singular value sequence. The matrices $U$ and $V$ are unique up to a unitary block on each constant-singular-value subspace: replacing $V$ by $VW$ with $W$ a unitary block-diagonal matrix preserving singular-value-equality blocks, and correspondingly $U$ by $U{W}^{\prime }$ on the matched blocks of left singular vectors. For simple non-zero singular values, the freedom reduces to a phase $e^{i\theta }$ (a $\pm 1$ sign over $\mathbb{R}$).

Existence of polar decomposition. Starting from the SVD $A=U\Sigma V^{\ast }$ of a square $n\times n$ matrix, factor as $A=(U{V}^{\ast })(V\Sigma V^{\ast })$. Set $Q=U{V}^{\ast }$ (unitary) and $P=V\Sigma V^{\ast }$. Then $P$ is Hermitian ($P^{\ast }=V{\Sigma }^{\ast }{V}^{\ast }=V\Sigma V^{\ast }=P$, since $\Sigma$ is real diagonal) and positive semidefinite ($w^{\ast }Pw=w^{\ast }V\Sigma V^{\ast }w=\Vert {\Sigma }^{1/2}{V}^{\ast }w{\Vert }^2\ge 0$). The identity $P^2=V{\Sigma }^2{V}^{\ast }=V\Sigma V^{\ast }V\Sigma V^{\ast }=A^{\ast }A$ (after expanding using $V^{\ast }V=I$ and $A^{\ast }A=V{\Sigma }^2{V}^{\ast }$) confirms $P=(A^{\ast }A{)}^{1/2}$, the unique Hermitian positive semidefinite square root.

Uniqueness of polar decomposition for invertible $A$. Suppose $A=Q_1{P}_1=Q_2{P}_2$ with $Q_i$ unitary and $P_i$ Hermitian positive semidefinite. Then $A^{\ast }A=P_1^{\ast }{Q}_1^{\ast }{Q}_1{P}_1=P_1^2=P_2^2$. The Hermitian positive semidefinite square root is unique, so $P_1=P_2=(A^{\ast }A{)}^{1/2}$. When $A$ is invertible, $P_1$ is invertible and $Q_1=A{P}_1^{-1}=A{P}_2^{-1}=Q_2$, so the unitary factor is forced as well.

Moore-Penrose pseudoinverse characterisation. Let $A=U\Sigma V^{\ast }$ and set $A^{+}=V{\Sigma }^{+}{U}^{\ast }$ with ${\Sigma }^{+}$ the rectangular diagonal matrix with ${\Sigma }_{ii}^{+}=1/{\sigma }_i$ for ${\sigma }_i>0$ and zero otherwise. Verify the four Moore-Penrose conditions:

• $A{A}^{+}A=U\Sigma V^{\ast }V{\Sigma }^{+}{U}^{\ast }U\Sigma V^{\ast }=U\Sigma {\Sigma }^{+}\Sigma V^{\ast }=U\Sigma V^{\ast }=A$, using $V^{\ast }V=I$, $U^{\ast }U=I$, and the diagonal identity $\Sigma {\Sigma }^{+}\Sigma =\Sigma$ (the rectangular diagonal matrix satisfies ${\Sigma }_{ii}({\Sigma }^{+}{)}_{ii}{\Sigma }_{ii}={\Sigma }_{ii}$ whether ${\Sigma }_{ii}=0$ or not).
• $A^{+}A{A}^{+}=V{\Sigma }^{+}{U}^{\ast }U\Sigma V^{\ast }V{\Sigma }^{+}{U}^{\ast }=V{\Sigma }^{+}\Sigma {\Sigma }^{+}{U}^{\ast }=V{\Sigma }^{+}{U}^{\ast }=A^{+}$ by the analogous diagonal identity ${\Sigma }^{+}\Sigma {\Sigma }^{+}={\Sigma }^{+}$.
• $(A{A}^{+}{)}^{\ast }=(U\Sigma {\Sigma }^{+}{U}^{\ast }{)}^{\ast }=U(\Sigma {\Sigma }^{+}{)}^{\ast }{U}^{\ast }=U\Sigma {\Sigma }^{+}{U}^{\ast }=A{A}^{+}$, since $\Sigma {\Sigma }^{+}$ is a real diagonal matrix with entries in $\{0,1\}$, hence Hermitian.
• $(A^{+}A{)}^{\ast }=(V{\Sigma }^{+}\Sigma V^{\ast }{)}^{\ast }=V{\Sigma }^{+}\Sigma V^{\ast }=A^{+}A$ by the same argument.

Uniqueness follows from the four conditions: any $X$ satisfying $AXA=A$, $XAX=X$, $(AX{)}^{\ast }=AX$, $(XA{)}^{\ast }=XA$ equals $A^{+}$, because $X=XAX=X(AX{)}^{\ast }=X{X}^{\ast }{A}^{\ast }=X{X}^{\ast }(A{A}^{+}A{)}^{\ast }=X(AX{)}^{\ast }(A^{+}{)}^{\ast }{A}^{\ast }=(XA{)}^{\ast }(A^{+}{)}^{\ast }{A}^{\ast }\cdots$ — a sequence of substitutions using the four conditions forces $X=A^{+}$.

Eckart-Young in Frobenius norm. Let $A={\sum }_{i=1}^r{\sigma }_i{u}_i{v}_i^{\ast }$ and $A_k={\sum }_{i=1}^k{\sigma }_i{u}_i{v}_i^{\ast }$. Suppose $B$ has rank at most $k$; the goal is $\Vert A-B{\Vert }_F\ge \Vert A-A_k{\Vert }_F$. The kernel of $B$ has dimension $\ge n-k$, so there exist orthonormal vectors $w_1,\dots ,w_{n-k}$ in $\ker B\cap \mathrm{span}\{v_1,\dots ,v_{k+1}{\}}^{\ast }$ — by dimension count, $\ker B$ intersects $\mathrm{span}\{v_1,\dots ,v_{k+1}\}$ in a subspace of dimension at least $(n-k)+(k+1)-n=1$. Pick a unit vector $w\in \ker B\cap \mathrm{span}\{v_1,\dots ,v_{k+1}\}$, write $w={\sum }_{i=1}^{k+1}{\alpha }_i{v}_i$ with $\sum |{\alpha }_i{|}^2=1$. Then

$$\Vert A-B{\Vert }_{\mathrm{op}}^2\ge \Vert (A-B)w{\Vert }^2=\Vert Aw{\Vert }^2={\Vert \sum _{i=1}^{k+1}{\alpha }_i{\sigma }_i{u}_i\Vert }^2=\sum _{i=1}^{k+1}|{\alpha }_i{|}^2{\sigma }_i^2\ge {\sigma }_{k+1}^2,$$

using ${\sigma }_i\ge {\sigma }_{k+1}$ for $i\le k+1$ and $\sum |{\alpha }_i{|}^2=1$. So $\Vert A-B{\Vert }_{\mathrm{op}}\ge {\sigma }_{k+1}=\Vert A-A_k{\Vert }_{\mathrm{op}}$, proving Eckart-Young in operator norm.

For the Frobenius norm, the same dimension-count argument generalises. Choose orthonormal vectors $w_1,\dots ,w_{r-k}$ in $\ker B\cap \mathrm{span}\{v_1,\dots ,v_r\}$ (dimension at least $(n-k)+r-n=r-k$ when $r\le n$). Then ${\sum }_j\Vert (A-B)w_j{\Vert }^2={\sum }_j\Vert A{w}_j{\Vert }^2$, and the right-hand side is bounded below by ${\sigma }_{k+1}^2+{\sigma }_{k+2}^2+\cdots +{\sigma }_r^2$ via a Courant-Fischer-type minimax argument on the singular value sum. Hence $\Vert A-B{\Vert }_F^2\ge {\sigma }_{k+1}^2+\cdots +{\sigma }_r^2=\Vert A-A_k{\Vert }_F^2$.

The Mirsky extension to arbitrary unitarily invariant norms uses the Ky Fan dominance theorem: a vector $x$ majorises $y$ iff $\varphi (x)\ge \varphi (y)$ for every symmetric gauge function $\varphi$. Since the singular value sequence of $A-B$ majorises $(0,\dots ,0,{\sigma }_{k+1},\dots ,{\sigma }_r,0,\dots )$, every unitarily invariant norm achieves its minimum on the Eckart-Young truncation $A_k$.

Connections [Master]

• Eigenvalue, eigenvector, characteristic polynomial 01.01.08 — supplies the spectral theorem on Hermitian operators that drives the existence proof of the SVD. The right singular vectors of $A$ are the eigenvectors of $A^{\ast }A$; the singular values are the non-negative square roots of the eigenvalues of $A^{\ast }A$; the left singular vectors are the eigenvectors of $A{A}^{\ast }$. The SVD is therefore the linear-algebra construction that converts the spectral theorem on the Hermitian positive semidefinite Gram matrix into a factorisation of $A$ itself, applicable to arbitrary rectangular matrices. The Cayley-Hamilton identity from 01.01.08 specialises on $A^{\ast }A$ to give finite-dimensional polynomial identities for the singular values, including $\det (tI-A^{\ast }A)={\prod }_i(t-{\sigma }_i^2)$.

• Jordan canonical form and minimal polynomial 01.01.11 — the SVD and the Jordan canonical form are the two universal classifications of matrices, parallel and complementary. The Jordan form is the complete invariant of an $n\times n$ matrix under conjugation by ${\mathrm{GL}}_n(\mathbb{C})$; the SVD is the complete invariant of an $m\times n$ matrix under the bi-unitary action of $U(m)\times U(n)$. The Jordan form respects the algebraic structure of the operator (eigenvalues, generalised eigenspaces, nilpotent part); the SVD respects the geometric structure of the linear map between inner-product spaces (singular values as semi-axes of the image ellipsoid). For Hermitian positive semidefinite matrices, the two coincide: the Jordan form is diagonal, the SVD has $U=V$, and the singular values equal the eigenvalues. For everything else they diverge — and both diverge from the eigendecomposition picture.

• Inner-product space: orthogonality, Gram-Schmidt, spectral theorem 01.01.09 pending — supplies the orthonormalisation machinery and the finite-dimensional spectral theorem on which the SVD existence proof depends. The Gram-Schmidt procedure of 01.01.09 pending is used in the existence proof to extend the partial orthonormal set $\{u_1,\dots ,u_r\}$ of left singular vectors to a full orthonormal basis of ${\mathbb{F}}^m$, completing the unitary $U$. The spectral theorem on the Hermitian operator $A^{\ast }A$ is the load-bearing input. The SVD is the rectangular extension of the spectral theorem and reduces to the spectral theorem when $A$ is itself Hermitian positive semidefinite.

• Bounded and unbounded operators on Hilbert space 02.11.03 — the infinite-dimensional generalisation in which the finite-dimensional SVD becomes the Schmidt decomposition for compact operators on a Hilbert space, the polar decomposition for bounded operators (with $A=Q|A|$ for $Q$ a partial isometry and $|A|=(A^{\ast }A{)}^{1/2}$), and the spectral measure for self-adjoint operators (Stone-von Neumann). The singular numbers of a compact operator $T$ — the non-zero eigenvalues of $|T|$ — generalise singular values and underlie the Schatten $p$-class operators, the trace-class operators ($p=1$), the Hilbert-Schmidt operators ($p=2$), and the operator algebras of non-commutative integration. The finite-dimensional SVD is the rank-finite, no-limiting-condition shadow of this hierarchy.

• Least-squares regression and the normal equations [TODO future unit on regression] — the SVD-based pseudoinverse $A^{+}=V{\Sigma }^{+}{U}^{\ast }$ solves the general least-squares problem ${\min }_x\Vert Ax-b{\Vert }_2$ in one closed-form expression. The classical normal equations $A^{\ast }Ax=A^{\ast }b$ recover the same solution when $A$ has full column rank, but break down when $A$ is rank-deficient or ill-conditioned. The pseudoinverse handles both cases uniformly: it returns the minimum-norm least-squares solution, projecting $b$ onto the column space of $A$ and resolving the kernel degree of freedom by minimising $\Vert x\Vert$. Tikhonov regularisation ${\min }_x\Vert Ax-b{\Vert }^2+{\lambda }^2\Vert x{\Vert }^2$ is equivalent to filtering the singular values, replacing ${\sigma }_i$ by ${\sigma }_i^2/({\sigma }_i^2+{\lambda }^2)$ in the pseudoinverse expansion, and provides numerical stability when ${\sigma }_r$ is small.

• Principal component analysis and data dimensionality reduction [TODO future unit] — for a centred data matrix $X\in M_{N\times d}(\mathbb{R})$ with $N$ samples in ${\mathbb{R}}^d$, the right singular vectors of $X$ (columns of $V$) are the principal directions of the data, ordered by decreasing variance ${\sigma }_i^2/N$. The $k$-truncated SVD $X_k=U_k{\Sigma }_k{V}_k^T$ is the best rank-$k$ approximation in Frobenius norm by Eckart-Young; projecting the data onto the span of $\{v_1,\dots ,v_k\}$ retains the maximum total variance among all $k$-dimensional linear subspaces. PCA is therefore the SVD of the data matrix, with the singular values quantifying the variance explained. Variants — kernel PCA, sparse PCA, robust PCA — replace the Euclidean inner product, the rank-$k$ constraint, or the assumed Gaussian noise model, but all share the SVD-of-a-data-matrix scaffolding.

Historical & philosophical context [Master]

The singular value decomposition has independent origins in the work of Eugenio Beltrami and Camille Jordan in the early 1870s. Beltrami, in Sulle funzioni bilineari (Giornale di Matematiche ad Uso degli Studenti Delle Università 11, 1873, 98–106), proved the existence of a bilinear-form factorisation of a square real matrix into two orthogonal transformations and a diagonal matrix of non-negative entries ^{[Beltrami, E. — Sulle funzioni bilineari]}. One year later, Camille Jordan, in Mémoire sur les formes bilinéaires (Journal de Mathématiques Pures et Appliquées, 2e série, 19, 1874, 35–54), gave an independent treatment for square matrices using a variational approach, characterising the singular values as critical values of the bilinear form $⟨Ax,y⟩$ on the product of two unit spheres ^{[Jordan, C. — Mémoire sur les formes bilinéaires]}. James Joseph Sylvester, in Sur la réduction biorthogonale d'une forme linéo-linéaire à sa forme canonique (Comptes Rendus de l'Académie des Sciences 108, 1889, 651–653), extended the decomposition to rectangular matrices and named the diagonal entries the "canonical multipliers" of the form ^{[Sylvester, J. J. — Sur la réduction biorthogonale d'une forme linéo-linéaire à sa forme canonique]}.

The infinite-dimensional generalisation came from Erhard Schmidt in Zur Theorie der linearen und nichtlinearen Integralgleichungen, I. Teil (Mathematische Annalen 63, 1907, 433–476), in the context of Hilbert's programme on integral equations. Schmidt proved that every compact (then called "completely continuous") integral operator on $L^2$ admits the decomposition $T={\sum }_n{\sigma }_n{u}_n\otimes v_n^{\ast }$ with singular numbers ${\sigma }_n\to 0$ — the Schmidt decomposition of compact operators ^{[Schmidt, E. — Zur Theorie der linearen und nichtlinearen Integralgleichungen, I. Teil]}. Hermann Weyl, in Das asymptotische Verteilungsgesetz der Eigenwerte linearer partieller Differentialgleichungen (Mathematische Annalen 71, 1912, 441–479) and subsequent papers, unified the spectral perspectives on Hermitian and non-Hermitian operators and gave the Weyl inequalities relating singular values to eigenvalues of arbitrary square matrices, central to the modern theory of matrix inequalities and majorisation ^{[Weyl, H. — Das asymptotische Verteilungsgesetz der Eigenwerte linearer partieller Differentialgleichungen]}.

The low-rank approximation theorem was proved by Carl Eckart and Gale Young in The approximation of one matrix by another of lower rank (Psychometrika 1, 1936, 211–218), motivated by the factor analysis of psychological-test data; the theorem identifies the truncated SVD as the optimal rank-$k$ approximation in Frobenius norm ^{[Eckart, C. & Young, G. — The approximation of one matrix by another of lower rank]}. Leon Mirsky, in Symmetric gauge functions and unitarily invariant norms (Quarterly Journal of Mathematics, Oxford Series, 11, 1960, 50–59), extended Eckart-Young to arbitrary unitarily invariant norms using von Neumann's earlier symmetric-gauge-function theorem ^{[Mirsky, L. — Symmetric gauge functions and unitarily invariant norms]}.

The modern numerical algorithm for computing the SVD was developed by Gene Golub and William Kahan in Calculating the singular values and pseudo-inverse of a matrix (Journal of the Society for Industrial and Applied Mathematics, Series B: Numerical Analysis 2, 1965, 205–224) and refined in Golub-Reinsch 1970, with bidiagonalisation by Householder reflections followed by a specialised QR sweep ^{[Golub, G. H. & Kahan, W. — Calculating the singular values and pseudo-inverse of a matrix]}. The Golub-Kahan algorithm and its descendants are implemented in LAPACK and underlie every numerical linear-algebra package. The SVD entered data analysis through factor analysis in psychometrics (Eckart-Young 1936; Tucker 1966), principal component analysis (Hotelling 1933 for the covariance-matrix eigenvalue formulation, recast via SVD in the second half of the twentieth century), and latent semantic indexing (Deerwester et al. 1990 for information retrieval). Recommender systems, image compression (JPEG-2000 uses SVD-adjacent transforms), and large-scale machine learning all rely on truncated-SVD computations.

Bibliography [Master]

Nineteenth-century origins.

• Beltrami, E., "Sulle funzioni bilineari", Giornale di Matematiche ad Uso degli Studenti Delle Università 11 (1873), 98–106.
• Jordan, C., "Mémoire sur les formes bilinéaires", Journal de Mathématiques Pures et Appliquées, 2e série, 19 (1874), 35–54.
• Sylvester, J. J., "Sur la réduction biorthogonale d'une forme linéo-linéaire à sa forme canonique", Comptes Rendus de l'Académie des Sciences 108 (1889), 651–653.

Hilbert-space generalisation.

• Schmidt, E., "Zur Theorie der linearen und nichtlinearen Integralgleichungen, I. Teil", Mathematische Annalen 63 (1907), 433–476.
• Weyl, H., "Das asymptotische Verteilungsgesetz der Eigenwerte linearer partieller Differentialgleichungen", Mathematische Annalen 71 (1912), 441–479.

Low-rank approximation.

• Eckart, C. & Young, G., "The approximation of one matrix by another of lower rank", Psychometrika 1 (1936), 211–218.
• Mirsky, L., "Symmetric gauge functions and unitarily invariant norms", Quarterly Journal of Mathematics, Oxford Series, 11 (1960), 50–59.
• von Neumann, J., "Some matrix inequalities and metrization of matric-space", Tomsk University Review 1 (1937), 286–300.

Numerical computation.

• Golub, G. H. & Kahan, W., "Calculating the singular values and pseudo-inverse of a matrix", Journal of the Society for Industrial and Applied Mathematics, Series B: Numerical Analysis 2 (1965), 205–224.
• Golub, G. H. & Reinsch, C., "Singular value decomposition and least squares solutions", Numerische Mathematik 14 (1970), 403–420.
• Golub, G. H. & Van Loan, C. F., Matrix Computations, 4th ed., Johns Hopkins University Press, 2013.
• Trefethen, L. N. & Bau, D., Numerical Linear Algebra, SIAM, 1997.

Modern textbook treatments.

• Horn, R. A. & Johnson, C. R., Matrix Analysis, 2nd ed., Cambridge University Press, 2013, Ch. 7.
• Strang, G., Introduction to Linear Algebra, 5th ed., Wellesley-Cambridge Press, 2016, Ch. 7.
• Stewart, G. W., "On the early history of the singular value decomposition", SIAM Review 35 (1993), 551–566.
• Bhatia, R., Matrix Analysis, Graduate Texts in Mathematics 169, Springer, 1997.

Applications.

• Hotelling, H., "Analysis of a complex of statistical variables into principal components", Journal of Educational Psychology 24 (1933), 417–441, 498–520.
• Deerwester, S., Dumais, S. T., Furnas, G. W., Landauer, T. K., & Harshman, R., "Indexing by latent semantic analysis", Journal of the American Society for Information Science 41 (1990), 391–407.
• Candès, E. J. & Recht, B., "Exact matrix completion via convex optimization", Foundations of Computational Mathematics 9 (2009), 717–772.

---