The formula
Every real m×n matrix, no exceptions, factors as
A = U·Σ·Vᵀ
Read it right to left: Vᵀ rotates the input space, Σ stretches along the axes by the singular values, U rotates the result into the output space. U (m×m) and V (n×n) are orthogonal; Σ is diagonal with entries σ₁ ≥ σ₂ ≥ … ≥ 0.
The singular values are square roots of the eigenvalues of AᵀA:
σ_i = √( λ_i(AᵀA) )
That line is where the SVD comes from — but not how you should compute it (see below).
Three diagnostics fall straight out of Σ:
‖A‖₂ = σ₁ rank(A) = number of nonzero σ_i cond(A) = σ₁/σ_n
Largest gain, effective rank, and the amplification factor for relative error, all read off one diagonal.
The pseudoinverse — the machinery underneath every least-squares solver worth trusting:
A⁺ = V·Σ⁺·Uᵀ, where Σ⁺ inverts each nonzero σ_i and zeros the rest
Then x = A⁺·b is the minimum-norm least-squares solution, and small singular values you chose to zero never get inverted into garbage.
And the theorem that makes SVD a compression and filtering tool. Keep only the first k terms:
A_k = σ₁·u₁·v₁ᵀ + σ₂·u₂·v₂ᵀ + … + σ_k·u_k·v_kᵀ
A_k is the closest rank-k matrix to A in both the 2-norm and Frobenius norm, and the error is exactly what you threw away: ‖A − A_k‖₂ = σ_(k+1). Nothing of rank k does better.
Where you meet it
- Test stand calibration review. A six-component strain-gauge balance calibration is a least-squares fit, and the fit matrix has a condition number.
cond = 4×10⁴means roughlylog₁₀(4×10⁴) ≈ 4.6digits of your data's accuracy are gone before the solver starts. When the review board asks why two calibration coefficients trade against each other run to run, the answer is a small singular value — the load schedule never independently exercised that direction. - Multichannel data cleanup. Stack 64 accelerometer channels from a vibration test into a matrix, run the SVD, and the singular values separate coherent structure from uncorrelated noise. A handful of large σ's carry the modal content; the flat tail is the noise floor. Truncate the tail and you have filtered every channel at once without picking a single cutoff frequency.
- Antenna and MIMO work. The SVD of a channel or coupling matrix decomposes an antenna system into independent sub-channels: the singular vectors are the transmit and receive beam pairs, and each σ is that pair's gain. Water-filling power allocation in a MIMO link is literally pouring power onto the largest singular values first.
- Image and pattern compression. A 1024×1024 measurement map stored as its first 30 singular triplets is a 17:1 reduction, and the Eckart–Young theorem tells you the exact worst-case error before you commit — it is σ₃₁.
How it works
Picture the unit sphere in the input space. A maps it to an ellipsoid. The singular values are the ellipsoid's semi-axis lengths; the columns of U point along those axes; the columns of V are the input directions that land there. Every matrix — rectangular, rank-deficient, ugly — does exactly this and nothing more. That is why the SVD works on the matrices eigendecomposition refuses: it never asks for square, symmetric, or diagonalizable.
A concrete one you can check by hand: the 3×2 matrix A = [[3,2],[2,3],[2,−2]] has singular values exactly 5 and 3. Its 2-norm is 5, its condition number is 5/3, and the best rank-1 approximation misses it by a 2-norm error of exactly 3.
Rank is a tolerance decision, not a count. On real hardware data no computed singular value is exactly zero. The honest question is where the gap sits. σ = {2.0, 5×10⁻⁵} (that's the actual spectrum of [[1,1],[1,1.0001]]) is a rank-2 matrix that behaves like rank 1 for any measurement with worse than five-digit accuracy. Numerical rank is "how many σ's stand above your noise," and defending that threshold is your job, not the algorithm's.
The classic mistake: forming AᵀA. The normal-equations shortcut — build AᵀA, eigendecompose or Cholesky it — squares the condition number. A fit matrix at cond = 10⁴ becomes a normal-equations matrix at cond = 10⁸, and in double precision (ε ≈ 2.2×10⁻¹⁶) you just spent half your digit budget on convenience. Golub and Kahan's algorithm exists specifically to avoid this: it bidiagonalizes A directly with Householder reflections and never forms the product [3].
Gotchas worth knowing.
- Singular vectors are not unique. Any
u_i, v_ipair can flip sign together, and when two singular values are equal (or merely close, on noisy data) the corresponding vectors can rotate freely within their subspace. Comparing singular vectors across two test runs without accounting for this generates phantom "changes." - The SVD is unit-blind. Mix millivolts and megapascals in the same data matrix and the singular values rank columns by numerical magnitude, not physical importance. Scale or standardize columns first, and record how — the decomposition changes with the scaling.
- Cost: dense SVD runs on the order of
m·n²flops. Fine at 64 channels, a planning problem at 10⁵. Truncated and randomized variants exist for the big cases; know which one your library called before quoting runtimes. - Truncation is a projection, not magic. Rank truncation removes what is uncorrelated across channels. Noise that is coherent across your channels — a shared ground loop, a common power rail — survives it looking exactly like signal.
History
The SVD was discovered the way a lot of good tools are: repeatedly, by people who didn't read each other. Eugenio Beltrami published it first, in 1873, in a journal written for Italian university students — a teaching exercise on bilinear forms that happened to found a subject [1][2]. Camille Jordan got there independently a year later with a cleaner and more complete treatment; Stewart's history calls the two of them the joint progenitors, Beltrami by priority, Jordan by rigor [1]. James Joseph Sylvester derived it a third time in 1889, apparently unaware of both, and named the singular values the "canonical multipliers" of the matrix [1][2].
The subject then jumped fields. In 1907 Erhard Schmidt — the Schmidt of Gram–Schmidt — built the infinite-dimensional version for integral equations, and proved the low-rank approximation theorem that makes the SVD useful rather than merely true [1][2]. Émile Picard, working in the same integral-equations literature, attached the name "singular values" around 1910 [1][2]. Hermann Weyl added the perturbation theory in 1912 [1], and Léon Autonne extended the decomposition to complex matrices in the 1910s by way of the polar decomposition [1][2]. In 1936 the physicists Carl Eckart and Gale Young proved it for rectangular matrices and re-proved Schmidt's approximation theorem, which has carried their names ever since — Stewart notes, dryly, incorrectly [1][2][5].
None of it was routinely computable until Gene Golub and William Kahan published a stable algorithm in 1965, built on Householder bidiagonalization instead of the condition-squaring AᵀA route [2][3]. The 1970 Golub–Reinsch refinement is, half a century later, still essentially what runs when you type svd [2][4].
Related tools
- /tools/antenna-gain-beamwidth — the per-beam gains an antenna-system SVD hands back, in familiar units
- /tools/strain-gauge-bridge — the channels whose calibration matrix conditioning the SVD diagnoses
- /tools/snr-enob — the noise floor that sets where your singular-value tail flattens out
- /tools/adc-resolution — quantization limits the digits you have before the condition number starts spending them