Every real matrix has a singular value decomposition. Not most matrices, not the square ones, not the invertible ones. A result proved by Beltrami in 1873 and made computable by the Golub–Kahan algorithm in 1965 guarantees that any m × n real matrix factors as A = UΣVᵀ, with U and V orthogonal and Σ nonnegative diagonal.
What is singular value decomposition?
Singular value decomposition splits a matrix into three factors with fixed geometric jobs.
V selects a set of orthonormal input directions. Σ scales each of those directions by a nonnegative number. U places the scaled results into an orthonormal set of output directions. Nothing else happens. Any linear map, no matter how ugly the matrix looks, is a rotation followed by an axis-aligned stretch followed by another rotation.
The nonnegative numbers on Σ are the singular values, written σ₁ ≥ σ₂ ≥ … ≥ σₚ ≥ 0 where p = min(m, n).
This is what separates SVD from eigendecomposition. Eigendecomposition needs a square matrix and still fails for defective ones. SVD handles a 3×7 matrix, a rank-1 matrix, and a matrix of zeros without a special case. Free reference sheets for the surrounding machinery sit in the linear algebra cheatsheet.
What does the SVD formula mean?
Write it once and read the shapes.
A = UΣVᵀ, where A is m × n, U is m × m orthogonal, Σ is m × n with σ₁ … σₚ on the main diagonal and zeros elsewhere, and V is n × n orthogonal. For complex matrices the transpose becomes a conjugate transpose, A = UΣV*, which is where polar notation for complex numbers starts paying off.
The shapes are the part students get wrong on exam scripts.
6×3
6×6
6×3
3×3
6×3
6×3
3×3
3×3
The hatched block in U holds m − n columns spanning the cokernel of A. They are genuine basis vectors, and for a tall data matrix they are almost always wasted work.
How do you compute an SVD?
By hand, the route runs through AᵀA, which is symmetric and positive semidefinite, so its eigenvalues are real and nonnegative.
Take A = [[3, 2, 2], [2, 3, −2]]. Because A is 2×3, working with AAᵀ = [[17, 8], [8, 17]] is faster than AᵀA. Its eigenvalues are 25 and 9, so σ₁ = 5 and σ₂ = 3.
The eigenvectors give u₁ = (1, 1)/√2 and u₂ = (1, −1)/√2. Then v₁ = Aᵀu₁/5 = (1, 1, 0)/√2 and v₂ = Aᵀu₂/3 = (1, −1, 4)/(3√2).
V still needs a third column, and this is the step most worked examples skip. Because A has rank 2 and three columns, its kernel is one-dimensional. Solving Ax = 0 gives v₃ = (−2, 2, 1)/3, which extends V to a full orthonormal basis of ℝ³.
Verify with numpy.linalg.svd and the singular values return as exactly 5 and 3. Signs on individual singular vectors may flip, which is not an error.
Does every matrix have an SVD?
Yes. Without exception, for every real or complex matrix of any shape.
The proof is short. AᵀA is symmetric positive semidefinite, so the spectral theorem hands over an orthonormal eigenbasis with nonnegative eigenvalues. Take square roots and the construction above always terminates. Nothing about A needs to be square, invertible, or diagonalisable.
Uniqueness is weaker than existence, and exams test the gap.
The singular values are always unique. The singular vectors are unique only up to sign when the singular values are distinct. Where two singular values coincide, any orthonormal basis of that shared subspace works, so infinitely many valid U and V exist for the same A.
unit circle, ‖x‖ = 1
rotate, circle unchanged
stretch by 4 and by 2
rotate, long axis on y = x
Which SVD variant should you compute?
Four forms appear in coursework and in library documentation. They differ only in how much of U and V gets stored.
The gap widens fast. At 5,000×500 the full SVD costs 3.3 times the economy version. At 50,000×100 it costs 143 times as much, because the 4m²n term dominates everything else.
NumPy defaults to full_matrices=True. On a tall data matrix that default is the single most expensive line in most student scripts.
What do singular values actually tell you?
Four readings, all direct.
The count of nonzero singular values is the rank. σ₁ is the spectral norm ‖A‖₂. The ratio σ₁/σₙ is the 2-norm condition number, which bounds how badly a solve amplifies input error. And the squares σᵢ² partition the total Frobenius energy of A.
The 8×8 Hilbert matrix, Hᵢⱼ = 1/(i + j − 1), shows what a spectrum looks like when it collapses.
Read the last row carefully. σ₈ is roughly 1.1×10⁻¹⁰, which is tiny but far above the float64 rounding floor, so numpy.linalg.matrix_rank still reports rank 8. Its default threshold is σ₁ × max(m, n) × eps, about 3.0×10⁻¹⁵ here.
Mathematical rank and numerical rank are different objects. This matrix is full rank and effectively rank 2 at the same time.
Where does SVD show up in coursework?
Least squares first. The Moore–Penrose pseudoinverse A⁺ = VΣ⁺Uᵀ, where Σ⁺ inverts each nonzero singular value, returns the minimum-norm least-squares solution x = A⁺b even when AᵀA is singular and the normal equations fail.
Principal component analysis is the same computation wearing different notation. Centre the data matrix, take its SVD, and the right singular vectors are the principal directions while σᵢ²/(n − 1) gives the variance each one explains.
Compression follows from the Eckart–Young–Mirsky theorem, proved by Eckart and Young in Psychometrika in 1936 and extended by Mirsky to every unitarily invariant norm. Truncating at rank k is optimal, and the spectral-norm error equals σₖ₊₁ exactly.
The saving is concrete. A 1,920×1,080 matrix holds 2,073,600 entries; a rank-60 truncation holds 180,060, an 11.5-fold reduction. Related worked material sits in the AskSia Explore library and across the linear algebra course pages.
Frequently Asked Questions
What is the formula for singular value decomposition?
The formula is A = UΣVᵀ for real matrices and A = UΣV* for complex ones. For an m × n matrix A, U is m × m orthogonal, Σ is m × n with singular values σ₁ ≥ σ₂ ≥ … ≥ σₚ ≥ 0 on the main diagonal where p = min(m, n), and V is n × n orthogonal. An equivalent and often more useful form is the singular value expansion A = σ₁u₁v₁ᵀ + … + σᵣuᵣvᵣᵀ, a sum of exactly r rank-1 matrices where r = rank(A). The expansion form makes truncation obvious: stop the sum early and you have the best rank-k approximation available under any unitarily invariant norm. The singular values relate to eigenvalues by σᵢ = √λᵢ, where λᵢ are the eigenvalues of AᵀA. Practise the expansion form first, since it is the version that appears in PCA, compression and least squares derivations.
How does singular value decomposition work?
Geometrically, A maps the unit sphere in ℝⁿ to an ellipsoid in ℝᵐ. Vᵀ rotates the sphere so its axes align with the directions A treats specially. Σ stretches axis i by σᵢ, turning the sphere into an axis-aligned ellipsoid. U rotates that ellipsoid into its final orientation. For A = [[3,1],[1,3]], the unit circle becomes an ellipse with semi-axes 4 and 2, long axis along the line y = x. Numerically the process does not follow the AᵀA hand method, because squaring a matrix squares its condition number and destroys small singular values. LAPACK instead bidiagonalises A with Householder reflections, then runs an implicit-shift QR iteration on the bidiagonal form, an approach introduced by Golub and Kahan in 1965. Confirm your hand result against numpy.linalg.svd, and treat sign differences on singular vectors as expected rather than as errors.
What are the properties of singular value decomposition?
Six properties carry most of the exam weight. Singular values are always real and nonnegative, even for complex matrices. The number of strictly positive singular values equals rank(A). σ₁ equals the spectral norm ‖A‖₂, and the Frobenius norm equals the square root of the sum of all σᵢ². The condition number κ₂(A) = σ₁/σₙ measures error amplification; the 8×8 Hilbert matrix reaches 1.53×10¹⁰, meaning roughly ten decimal digits of accuracy are lost in a solve. Singular vectors supply orthonormal bases for all four fundamental subspaces at once: column space, row space, kernel and cokernel. Finally, singular values are unique while singular vectors are unique only up to sign, and not even that when singular values repeat. Build a Flashcards deck on these six and let FSRS scheduling place the reviews before your final rather than the night before it.
What is a multilinear singular value decomposition?
A multilinear SVD, usually called the higher-order SVD or HOSVD, generalises the factorisation from matrices to tensors with three or more indices. It was formalised by De Lathauwer, De Moor and Vandewalle in 2000 and writes a tensor as a core tensor multiplied by an orthogonal factor matrix along each mode. The important caveat is that the Eckart–Young guarantee does not carry over. Truncating a HOSVD gives a quasi-optimal approximation, typically within a factor of √N of the best possible for an order-N tensor, but not the exact optimum that matrix truncation delivers. Best rank-k tensor approximations can fail to exist at all, a result established by de Silva and Lim in 2008. Treat HOSVD as a separate topic rather than an extension, and check whether your unit covers tensors before spending revision time on it.
What is the difference between SVD and eigendecomposition?
Eigendecomposition writes A = PDP⁻¹ using one basis for both input and output, which forces A to be square and diagonalisable. SVD uses two different orthonormal bases, V for the input space and U for the output space, which is why it exists for every matrix of every shape. The two coincide in one case: when A is symmetric positive semidefinite, its eigenvalues and singular values are identical and U = V. For a symmetric matrix with negative eigenvalues, the singular values are the absolute values, so A = [[−4, 0], [0, 2]] has eigenvalues −4 and 2 but singular values 4 and 2. Eigendecomposition also uses a non-orthogonal P in general, which can be badly conditioned; U and V are always perfectly conditioned. Reach for SVD whenever the matrix is rectangular, rank-deficient, or numerically suspect.
When is SVD the wrong tool?
Sparsity is the first limit. A dense SVD produces dense U and V, so factorising a sparse 10⁶ × 10⁶ matrix would demand more memory than the machine has. Iterative methods such as Lanczos bidiagonalisation or randomised range finders recover only the leading k singular triplets and leave the sparsity intact.
Stability is the second. Singular values are always well conditioned, but individual singular vectors are not. When σᵢ and σᵢ₊₁ are close, small perturbations rotate the corresponding vectors substantially, so any interpretation attached to a single principal component in that region is fragile.
Cost is the third. At roughly 21m³ flops for a square matrix, computing a full SVD to answer a question that rank, QR or a Cholesky factorisation could settle is the wrong trade.
And SVD describes linear structure only. If the relationship in your data curves, no orthogonal factorisation will find it, and a low-rank fit will simply be wrong in an orderly way.
