Singular Value Decomposition

The singular value decomposition (SVD) writes any real matrix as $A = U \Sigma V^\top$: an orthogonal matrix $V^\top$, a diagonal matrix $\Sigma$ of non-negative singular values, and another orthogonal matrix $U$. It exists for every matrix, square or rectangular, of any rank, and it exposes everything about a matrix that does not depend on the choice of orthonormal coordinates: its rank, the subspaces it reaches and ignores, how much it stretches in each direction, and how close it is to a matrix of lower rank.

Computer vision uses it constantly, usually in one line of code: the direct linear transform estimates a homography as the last right singular vector of a stacked system, the eight-point algorithm zeroes the smallest singular value of a fundamental matrix to make it rank 2, noisy rotations are snapped back to the nearest rotation, and the factorization method recovers the shape and motion of an image sequence from its three largest singular values and their singular vectors. This article builds on linear algebra, which covers rank, null spaces and conditioning, and on least squares, which covers overdetermined systems and their solvers.

Definition

The singular value decomposition of a real $m \times n$ matrix $A$ is a factorization

$$A = U \Sigma V^\top,$$

where $U$ is an $m \times m$ orthogonal matrix, $V$ is an $n \times n$ orthogonal matrix, and $\Sigma$ is an $m \times n$ matrix that is zero except for the entries $\sigma_1 \ge \sigma_2 \ge \dots \ge \sigma_{\min(m,n)} \ge 0$ on its main diagonal. The $\sigma_i$ are the singular values of $A$, the columns $\mathbf{u}_i$ of $U$ are its left singular vectors, and the columns $\mathbf{v}_i$ of $V$ are its right singular vectors (Deisenroth, Faisal and Ong, 2020, Theorem 4.22). Every real matrix has one.

Intuition

Read $A\mathbf{x} = U\Sigma V^\top \mathbf{x}$ from right to left. $V^\top$ rotates (or reflects) the input so that the directions $\mathbf{v}_i$ line up with the coordinate axes; $\Sigma$ stretches each axis by $\sigma_i$, dropping or adding dimensions when $m \ne n$; and $U$ rotates the result into place in the output space (Deisenroth, Faisal and Ong, 2020, §4.5.1). Every linear map, however it was written down, is a rotation, an axis-aligned scaling and another rotation.

Two pictures follow from this. Geometrically, the unit sphere of inputs is mapped to an ellipsoid (flattened to a lower-dimensional ellipse if some $\sigma_i$ are zero) whose semi-axes are $\sigma_i \mathbf{u}_i$: the right singular vectors are the input directions that the map sends to the axes of the ellipsoid, and the singular values are the axis lengths. For a $2 \times 2$ image warp, the unit circle becomes an ellipse whose longer axis has length $\sigma_1$ and shorter axis $\sigma_2$, so $\sigma_1 / \sigma_2$ measures how anisotropic the warp is.

Algebraically, the decomposition writes $A$ as a sum of rank-one pieces,

$$A = \sum_{i=1}^{r} \sigma_i \, \mathbf{u}_i \mathbf{v}_i^\top,$$

with $r$ the rank of $A$ (Deisenroth, Faisal and Ong, 2020, eq. 4.91). For a grayscale image stored as a matrix, each term $\mathbf{u}_i \mathbf{v}_i^\top$ is an image whose rows are all multiples of one row pattern $\mathbf{v}_i^\top$, weighted down the image by the column pattern $\mathbf{u}_i$. The singular values say how much of the image each such pattern carries, and they usually fall quickly: a few terms reproduce most of the image, which is the basis of low-rank approximation.

Formal Definition

Notation follows the rest of the site: $A = U \Sigma V^\top$ with singular values in non-increasing order, and $\mathbf{u}_i$, $\mathbf{v}_i$ the columns of $U$ and $V$.

Full and thin forms

The definition above is the full SVD: $U$ and $V$ are square, and $\Sigma$ has the shape of $A$, padded with zero rows (if $m > n$) or zero columns (if $m < n$). The columns of $U$ beyond the first $\min(m, n)$ multiply only zeros, so they can be dropped without changing the product. The thin (or reduced) SVD keeps $k = \min(m, n)$ columns: $U$ is $m \times k$ with orthonormal columns, $\Sigma$ is $k \times k$ and diagonal, and $V$ is $n \times k$ (Deisenroth, Faisal and Ong, 2020, §4.5.3). Keeping only the $r$ non-zero singular values gives the compact form, $U_r \Sigma_r V_r^\top$ with $r$ columns each. The thin form costs less memory and time when $A$ is far from square, as a $2n \times 9$ DLT matrix with hundreds of rows is.

NumPy’s np.linalg.svd(A) returns three arrays U, S, Vh: S is a one-dimensional array of singular values in descending order, and Vh is $V^\top$, not $V$, so the right singular vectors are the rows of Vh (NumPy 2.3 documentation). The default full_matrices=True gives the full form; full_matrices=False gives the thin one, with shapes $(m, k)$, $(k,)$ and $(k, n)$. SciPy’s scipy.linalg.svd follows the same conventions (SciPy 1.16 documentation). A common SVD bug in vision code is taking Vh[:, -1], a column of $V^\top$, instead of Vh[-1], the last right singular vector.

Relation to the eigendecompositions of $A^\top A$ and $AA^\top$

Substituting the decomposition into the two products gives

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

both symmetric eigendecompositions. So the right singular vectors are eigenvectors of $A^\top A$, the left singular vectors are eigenvectors of $AA^\top$, and the non-zero singular values are the square roots of the non-zero eigenvalues of either product (Deisenroth, Faisal and Ong, 2020, §4.5.2–§4.5.3). The two families are linked by

$$A \mathbf{v}_i = \sigma_i \mathbf{u}_i, \qquad A^\top \mathbf{u}_i = \sigma_i \mathbf{v}_i,$$

which is how the decomposition is constructed in textbooks: diagonalize $A^\top A$ with the spectral theorem, then set $\mathbf{u}_i = A\mathbf{v}_i / \sigma_i$. Beltrami’s 1873 derivation, the first published, worked this way through $AA^\top$ and $A^\top A$, for a square, non-singular matrix with distinct singular values; Jordan’s of 1874, independent and complete, characterized the largest singular value as the maximum of $\mathbf{x}^\top A \mathbf{y}$ over unit vectors and deflated the problem one value at a time (Stewart, 1992). Sylvester rediscovered the decomposition in 1889, and Schmidt (1907) gave the infinite-dimensional version for integral equations, with the approximation theorem below. For a symmetric positive semidefinite matrix the SVD and the eigendecomposition coincide.

The construction is a proof, not an algorithm. Forming $A^\top A$ squares the singular values, so information in the small ones is lost to rounding (Deisenroth, Faisal and Ong note the poor numerical behaviour of this route, §4.5.2); libraries work on $A$ directly, as the section on computation describes.

Uniqueness

The singular values are unique. When they are distinct and non-zero, each pair $(\mathbf{u}_i, \mathbf{v}_i)$ is unique up to a simultaneous change of sign, since flipping both leaves $\sigma_i \mathbf{u}_i \mathbf{v}_i^\top$ unchanged; this follows from the eigenvector relation above. When a singular value is repeated, any orthonormal basis of the corresponding subspace will do, and different libraries or versions return different ones. Code that compares singular vectors across runs, or reads geometric meaning into the sign of one, has to fix the sign or the basis itself, as pose recovery from the essential matrix does by requiring $\det U = \det V = +1$.

Properties

Rank, subspaces and numerical rank

The rank of $A$ is the number of non-zero singular values, and the decomposition hands over orthonormal bases of the four fundamental subspaces: the first $r$ columns of $U$ span the column space, the remaining columns of $U$ its orthogonal complement, the first $r$ columns of $V$ the row space, and the last $n - r$ columns of $V$ the null space. This follows from $A\mathbf{v}_i = \sigma_i \mathbf{u}_i$, which sends the first $r$ right singular vectors onto the first $r$ left ones and the rest to zero (Deisenroth, Faisal and Ong, 2020, §4.5.2, state the null-space part).

In floating point, rank becomes a question of thresholds, as the linear algebra article explains, and the SVD is the reason the threshold can be set sensibly. Weyl’s perturbation theorem, in the form Stewart (1992) gives for unsymmetric matrices, says that adding a perturbation $E$ moves each singular value by at most $\lVert E \rVert_2$, the largest singular value of $E$. A matrix of exact rank $k$ measured with noise of spectral norm $\lVert E \rVert_2$ therefore has $k$ singular values at least $\sigma_k - \lVert E \rVert_2$ and the rest at most $\lVert E \rVert_2$: a gap in the singular values that a threshold can be placed in. Tomasi and Kanade (1992) used exactly this test for their measurement matrix, requiring the ratio of the third to the fourth singular value to be large.

Norms and conditioning

The singular values determine the orthogonally invariant matrix norms. The spectral norm is the largest, $\lVert A \rVert_2 = \sigma_1$ (Deisenroth, Faisal and Ong, 2020, Theorem 4.24), and the Frobenius norm is $\lVert A \rVert_F = \sqrt{\sum_i \sigma_i^2}$ (Stewart, 1992), because orthogonal factors do not change the sum of squared entries. For a square invertible matrix, $\kappa(A) = \sigma_1 / \sigma_n$ is the spectral condition number that linear algebra uses to bound error amplification; for a rectangular matrix of full column rank the same ratio $\sigma_1 / \sigma_n$ serves as its condition number. A homography or fundamental-matrix system built from raw pixel coordinates has a huge $\sigma_1 / \sigma_{n-1}$ (Hartley works with its square, a ratio of eigenvalues of $A^\top A$), which is what Hartley’s normalization reduces (normalized eight-point paper).

Pseudoinverse and minimum-norm least squares

With the compact SVD $A = U_r \Sigma_r V_r^\top$, the Moore–Penrose pseudoinverse is

$$A^{+} = V_r \Sigma_r^{-1} U_r^\top = \sum_{i=1}^{r} \frac{1}{\sigma_i} \mathbf{v}_i \mathbf{u}_i^\top,$$

which inverts the non-zero singular values and leaves the zero ones at zero (NumPy 2.3 pinv documentation). Stewart (1992) attributes the pseudoinverse to Moore in 1920 and Penrose in 1955. For a square invertible matrix $A^{+} = A^{-1}$; for a matrix with full column rank, $A^{+} = (A^\top A)^{-1} A^\top$, the matrix of the normal equations.

The vector $\hat{\mathbf{x}} = A^{+}\mathbf{b}$ is the minimum-norm least-squares solution: among all vectors that minimize $\lVert A\mathbf{x} - \mathbf{b} \rVert$, it is the one with the smallest $\lVert \mathbf{x} \rVert$. The reason is visible in the coordinates the SVD provides. Writing $\mathbf{x} = V\mathbf{y}$ and $\mathbf{c} = U^\top \mathbf{b}$, the residual becomes $\lVert \Sigma \mathbf{y} - \mathbf{c} \rVert$, which decouples: $y_i = c_i / \sigma_i$ for $\sigma_i > 0$, and the remaining $y_i$, the null-space components, do not affect the residual at all. Setting them to zero minimizes $\lVert \mathbf{x} \rVert = \lVert \mathbf{y} \rVert$. Only when $A$ has a non-trivial null space is “minimum norm” a real choice; with full column rank the least-squares solution is unique and the pseudoinverse simply computes it. The least squares article compares this SVD solver with QR and the normal equations.

The same coordinates show the danger. A singular value of $10^{-10}$ contributes $c_i / \sigma_i$, amplifying whatever noise is in $c_i$ ten billion times. Solvers therefore treat singular values below a threshold as zero, which is a truncated SVD solution and a form of regularization. The defaults differ: NumPy’s pinv zeroes singular values at or below $10^{-15} \sigma_1$, while lstsq uses machine precision times $\max(m, n)$ times $\sigma_1$ (NumPy 2.3 documentation). For measured data, the threshold should come from the noise level, as with numerical rank.

Low-rank approximation: the Eckart–Young theorem

Truncating the sum of rank-one terms after $k$ terms gives

$$A_k = \sum_{i=1}^{k} \sigma_i \mathbf{u}_i \mathbf{v}_i^\top = U_k \Sigma_k V_k^\top,$$

and $A_k$ is the best approximation of $A$ by a matrix of rank at most $k$:

$$\min_{\operatorname{rank}(B) \le k} \lVert A - B \rVert_F = \lVert A - A_k \rVert_F = \sqrt{\sigma_{k+1}^2 + \dots + \sigma_{r}^2}, \qquad \min_{\operatorname{rank}(B) \le k} \lVert A - B \rVert_2 = \sigma_{k+1}.$$

The error of the best rank-$k$ approximation is exactly what the discarded singular values carry. Eckart and Young (1936) related the problem to the fundamental postulate of factor theory, posed it as a least-squares problem, and simplified its solution by first expressing the matrices in a canonical form, according to their abstract. The theorem is usually called the Eckart–Young theorem, but Stewart (1992) shows that Schmidt had proved it in 1907 for integral operators, calls the name incorrect, and gives Weyl’s 1912 perturbation argument as the elegant proof; Mirsky (1960) showed that $A_k$ minimizes the error in every unitarily invariant norm, which covers the spectral norm as well. Deisenroth, Faisal and Ong (2020, Theorem 4.25) state and sketch the spectral-norm case under Eckart and Young’s name.

The minimizer is unique when $\sigma_k > \sigma_{k+1}$; with a tie, any choice among the tied singular vectors is equally good (our reasoning; Eckart and Young describe the solution as usually unique). A rank-$k$ approximation of an $m \times n$ image needs $k(m + n + 1)$ numbers instead of $mn$, which is why textbooks use it to illustrate image compression (Deisenroth, Faisal and Ong, 2020, §4.6). Its real uses in vision are elsewhere: projecting a noisy matrix onto the rank its model requires (the fundamental matrix, the factorization method below), denoising data that should lie near a low-dimensional subspace, and principal component analysis, whose principal axes are the leading right singular vectors of the centred data matrix.

Homogeneous systems: the last right singular vector

Most linear estimators in multi-view geometry produce a homogeneous system $A\mathbf{h} = \mathbf{0}$, where $\mathbf{h}$ stacks the entries of a matrix defined only up to scale: the nine entries of a homography or a fundamental matrix, the twelve of a camera matrix, the four of a triangulated point. With noisy data there is no exact non-zero solution, and the estimate is

$$\hat{\mathbf{h}} = \arg\min_{\lVert \mathbf{h} \rVert = 1} \lVert A \mathbf{h} \rVert = \mathbf{v}_n,$$

the right singular vector of the smallest singular value, with minimum value $\sigma_n$. The proof takes one line in SVD coordinates: with $\mathbf{h} = V\mathbf{y}$ and $\lVert \mathbf{y} \rVert = 1$, $\lVert A\mathbf{h} \rVert^2 = \lVert \Sigma \mathbf{y} \rVert^2 = \sum_i \sigma_i^2 y_i^2$, which is smallest when all the weight is on the smallest $\sigma_i$. Equivalently, $\mathbf{v}_n$ is the eigenvector of $A^\top A$ for its smallest eigenvalue, $\sigma_n^2$, which is how some papers state it; computing it from $A$ avoids squaring the condition number.

This is the computation behind the direct linear transform for a homography and for the camera matrix in camera calibration, linear triangulation, the eight-point algorithm for the fundamental matrix, and the linear perspective-n-point solution. Two quantities say whether the answer can be trusted. The smallest singular value $\sigma_n$ is the algebraic residual, zero for exact data and growing with the noise. The gap to the next one, $\sigma_{n-1} - \sigma_n$, governs the accuracy of $\mathbf{v}_n$: LAPACK’s error analysis bounds the angle between a computed and a true singular vector by about machine precision times $\sigma_1$ divided by the gap to the nearest other singular value (LAPACK Users’ Guide, §4.9.1), and, by the same argument with the measurement noise in place of rounding error, the gap controls how far noise can rotate $\mathbf{v}_n$ (our reasoning). A degenerate configuration, such as three collinear points among four for a homography or a planar scene for the fundamental matrix, shows up as more than one singular value near zero: the null space has more than one dimension and the data do not determine the answer.

Projecting onto constrained matrices

A linear estimator ignores constraints that the true matrix satisfies, and the SVD is the standard tool for imposing them afterwards: change the singular values to what the constraint requires and keep the singular vectors. By Eckart–Young or by the orthogonal invariance of the Frobenius norm, this gives the nearest constrained matrix in the Frobenius norm in the following cases.

• Rank 2: the fundamental matrix. Setting $\sigma_3 = 0$ gives the nearest singular matrix, $U \operatorname{diag}(\sigma_1, \sigma_2, 0) V^\top$, at Frobenius distance $\sigma_3$. Hartley credits this correction to Tsai and Huang (1984) (normalized eight-point paper). Without it the epipolar lines of a noisy estimate do not meet in a single epipole.
• Two equal singular values: the essential matrix. A $3 \times 3$ matrix is an essential matrix exactly when two singular values are equal and the third is zero, and the nearest one to a linear estimate replaces the singular values with $(1, 1, 0)$ up to scale (essential matrix). The rotation and translation are then read off the same $U$ and $V$.
• Orthogonality: the nearest rotation. The orthogonal matrix nearest to $M = U\Sigma V^\top$ in the Frobenius norm replaces all singular values with 1, giving $UV^\top$. To exclude a reflection, the nearest rotation is $U \operatorname{diag}(1, 1, \det(UV^\top)) V^\top$, which flips the sign of the term with the smallest singular value when needed. Zhang’s method of camera calibration uses it to correct the rotation computed from each view’s homography, and rotation representations uses it to repair drift and to average rotations.

The last case has a more important relative. Given two sets of matched 3D points, centred at their centroids and stacked as rows of $P$ and $Q$, the rotation minimizing $\sum_i \lVert \mathbf{q}_i - R\mathbf{p}_i \rVert^2$ maximizes $\operatorname{tr}(R^\top Q^\top P)$, so it is the nearest rotation to the $3 \times 3$ cross-covariance $Q^\top P$. Without the determinant condition this is the orthogonal Procrustes problem, which arose in factor analysis; Stewart (1992) credits Green (1952) and Schönemann (1966) with its SVD solution. The closed-form alignment step of iterative closest point, with the determinant correction that guarantees a rotation, and the trajectory alignment behind absolute trajectory error, which adds a scale factor when the trajectory’s scale is unknown, are this computation.

Projection gives the nearest valid matrix in the Frobenius norm of its entries, an algebraic criterion that does not weigh the image measurements. Hartley noted that the rank-2 correction has little theoretical justification (normalized eight-point paper), and pipelines refine the projected estimate by minimizing a geometric error with nonlinear least squares or bundle adjustment.

Factorization of image sequences

Tomasi and Kanade (1992) applied low-rank approximation to structure from motion. They stacked the image coordinates of $P$ points tracked through $F$ frames into a $2F \times P$ measurement matrix $W$ and subtracted from each row its mean, the centroid of the points in that frame. Under orthographic projection, the registered matrix is the product of a $2F \times 3$ matrix of camera axes and a $3 \times P$ matrix of point coordinates, so their rank theorem states that without noise it has rank at most three. With noise, they took the best rank-3 approximation from the three largest singular values and their singular vectors, $\hat{W} = U_3 \Sigma_3 V_3^\top$, split it as $(U_3 \Sigma_3^{1/2})(\Sigma_3^{1/2} V_3^\top)$, and resolved the remaining invertible $3 \times 3$ ambiguity from the constraint that each frame’s two camera axes are orthogonal unit vectors. These details follow their full technical report on the orthographic case (Cornell TR 92-1270, March 1992); the journal version appeared in the International Journal of Computer Vision later that year. Multi-view geometry covers the extension to affine and projective cameras.

Computing the SVD

LAPACK, which NumPy and SciPy call, proceeds in two stages (LAPACK Users’ Guide, §2.4.6). First, Householder reflections applied from the left and right reduce $A$ to a bidiagonal matrix $B$, nonzero only on the diagonal and the first superdiagonal when $m \ge n$ (the first subdiagonal otherwise), with $A = U_1 B V_1^\top$; when $m$ is much larger than $n$ it can pay to take a QR factorization first and reduce only the $n \times n$ triangle. Second, the SVD of $B$ is computed iteratively, either by a QR iteration (routine xBDSQR) or by divide and conquer (xBDSDC), and the two sets of orthogonal factors are multiplied together. Stewart (1992) credits Golub and Kahan (1965) with introducing the SVD into numerical analysis with a computational algorithm, and the algorithm of Golub and Reinsch (1970) with being the workhorse of the following two decades.

The two drivers built on these stages are xGESVD, the original simple driver, and xGESDD, the divide-and-conquer driver added later, which is much faster for large matrices but needs more workspace (LAPACK Users’ Guide, §2.3.4.3, §3.4.3). NumPy’s np.linalg.svd uses _gesdd; SciPy’s scipy.linalg.svd also defaults to gesdd and offers lapack_driver='gesvd' as an alternative (NumPy 2.3 and SciPy 1.16 documentation). The bidiagonal reduction touches every entry of $A$ once per Householder step, and there are about $\min(m, n)$ steps, so the work grows roughly as $mn \min(m, n)$ for a dense matrix; this is reasoning from the algorithm’s structure rather than a quoted operation count. For the small matrices of geometric vision, a $2n \times 9$ DLT system or a $3 \times 3$ matrix, the cost is negligible; for a large matrix of which only the leading singular vectors are needed, iterative solvers such as SciPy’s scipy.sparse.linalg.svds compute just those.

The algorithms are backward stable: the computed decomposition is the exact SVD of a slightly perturbed matrix $A + E$ with $\lVert E \rVert_2$ of the order of machine precision times $\lVert A \rVert_2$ (LAPACK Users’ Guide, §4.9.1). Each computed singular value is therefore accurate to about $\varepsilon \sigma_1$ in absolute terms, so the large singular values have high relative accuracy and the small ones may not, and computed singular vectors are always orthogonal to working precision but individually accurate only when their singular value is well separated from the others. Computing the same quantities through the eigenvalues of $A^\top A$ loses more: an absolute error of $\varepsilon \sigma_1^2$ in $\sigma_i^2$ becomes a relative error of about $\varepsilon (\sigma_1 / \sigma_i)^2$ in $\sigma_i$, so singular values below $\sqrt{\varepsilon}\,\sigma_1$, about $10^{-8} \sigma_1$ in double precision, are lost entirely (our reasoning from the bound above).

Examples

The examples use NumPy 2.3.5 on synthetic images and synthetic correspondences with exactly known ground truth.

Low-rank approximation of an image

A $120 \times 160$ grayscale image is rendered from a horizontal gradient, a bright rectangle, a dark disk and a patch of diagonal stripes. The code checks the shapes of the full and thin decompositions, the relation to $A^\top A$, and the Eckart–Young errors of the truncations, then compares the truncation with other rank-10 approximations.

import numpy as np

# A 120 x 160 synthetic grayscale image: a horizontal gradient, a bright
# rectangle, a dark disk and a patch of diagonal stripes.
H, W = 120, 160
yy, xx = np.mgrid[0:H, 0:W]
img = 0.2 + 0.5 * xx / (W - 1)                              # gradient
img[20:60, 30:90] = 0.9                                     # rectangle
img[(yy - 80) ** 2 + (xx - 110) ** 2 < 25 ** 2] = 0.05      # disk
stripes = (yy > 70) & (yy < 110) & (xx > 10) & (xx < 60)
img[stripes] = 0.5 + 0.4 * np.sin(2 * np.pi * (xx + yy)[stripes] / 8)

U, s, Vt = np.linalg.svd(img)                             # full SVD; Vt is V^T
Ut, st, Vtt = np.linalg.svd(img, full_matrices=False)     # thin SVD
print("full:", U.shape, s.shape, Vt.shape, "  thin:", Ut.shape, Vtt.shape)
print("max |U S V^T - A| =", np.abs((Ut * st) @ Vtt - img).max())
lam = np.linalg.eigvalsh(img.T @ img)[::-1]               # eigenvalues of A^T A
print("sqrt of top eigenvalues of A^T A:", np.sqrt(lam[:3]).round(4))
print("top singular values:             ", s[:3].round(4))
print("rank:", np.linalg.matrix_rank(img), "  s_19, s_20 =", s[18].round(4), s[19])

for k in [1, 3, 10, 18]:
    Ak = (Ut[:, :k] * st[:k]) @ Vtt[:k]                   # rank-k truncation
    print(f"k={k:2d}  ||A-A_k||_F = {np.linalg.norm(img - Ak):7.4f}"
          f"  sqrt(sum_(i>k) s_i^2) = {np.sqrt((s[k:] ** 2).sum()):7.4f}"
          f"  ||A-A_k||_2 = {np.linalg.norm(img - Ak, 2):7.4f}  s_(k+1) = {s[k]:7.4f}"
          f"  numbers stored: {k * (H + W + 1) / (H * W):5.1%}")

# Another rank-10 approximation: project onto the span of 10 of the image's columns.
rng = np.random.default_rng(0)
for trial in range(3):
    Q, _ = np.linalg.qr(img[:, rng.choice(W, 10, replace=False)])
    print(f"projection onto 10 random image columns: ||A-B||_F = {np.linalg.norm(img - Q @ (Q.T @ img)):.4f}")

Output:

full: (120, 120) (120,) (160, 160)   thin: (120, 120) (120, 160)
max |U S V^T - A| = 5.9396931817445875e-15
sqrt of top eigenvalues of A^T A: [71.7531 14.8416 12.4584]
top singular values:              [71.7531 14.8416 12.4584]
rank: 19   s_19, s_20 = 0.5393 1.7036756869629975e-14
k= 1  ||A-A_k||_F = 24.2495  sqrt(sum_(i>k) s_i^2) = 24.2495  ||A-A_k||_2 = 14.8416  s_(k+1) = 14.8416  numbers stored:  1.5%
k= 3  ||A-A_k||_F = 14.5791  sqrt(sum_(i>k) s_i^2) = 14.5791  ||A-A_k||_2 =  8.8994  s_(k+1) =  8.8994  numbers stored:  4.4%
k=10  ||A-A_k||_F =  2.5042  sqrt(sum_(i>k) s_i^2) =  2.5042  ||A-A_k||_2 =  1.2133  s_(k+1) =  1.2133  numbers stored: 14.6%
k=18  ||A-A_k||_F =  0.5393  sqrt(sum_(i>k) s_i^2) =  0.5393  ||A-A_k||_2 =  0.5393  s_(k+1) =  0.5393  numbers stored: 26.3%
projection onto 10 random image columns: ||A-B||_F = 6.5533
projection onto 10 random image columns: ||A-B||_F = 10.8806
projection onto 10 random image columns: ||A-B||_F = 12.8019

In NumPy’s output Vt is $V^\top$: the full form returns a $160 \times 160$ orthogonal matrix whose last 40 rows meet only the zero columns of $\Sigma$, and the thin form drops them. The singular values agree with the square roots of the eigenvalues of $A^\top A$ to the printed digits, as they should for the large ones. The rendered image has rank 19, not 120: the gradient, the rectangle and the stripes are each sums of very few row patterns, and most of the rank comes from the disk, whose rows are all different. The twentieth singular value sits at rounding level, $10^{-14}$, the numerical signature of exact rank.

The truncation errors match the discarded singular values in both norms to all printed digits, which is the Eckart–Young theorem: in the Frobenius norm the error is the root of the sum of squares of the discarded values, and in the spectral norm it is the first discarded value. Ten terms, with under 15% of the numbers, leave a Frobenius error of 2.5 against an image norm of about 76. Projecting onto the span of ten of the image’s own columns, a reasonable rank-10 basis, leaves errors between 6.6 and 12.8 for the three draws shown, all worse than the truncated SVD, as the theorem guarantees for any rank-10 matrix.

A homography from the last right singular vector

Fifty points in a $640 \times 480$ image are mapped by a known homography, and the matches are perturbed by 1-pixel Gaussian noise. The DLT stacks two equations per match, after Hartley’s normalization of each point set, and takes the last right singular vector; the error is measured as the largest displacement of the image corners under the estimated homography relative to the true one.

import numpy as np

H_true = np.array([[0.9, 0.05, 30], [0.02, 1.1, 10], [0.0004, 0.0002, 1]])

def project(H, p):
    q = np.c_[p, np.ones(len(p))] @ H.T
    return q[:, :2] / q[:, 2:]

def normalize(p):
    """Similarity T moving the centroid to 0 and the mean distance to sqrt(2)."""
    c = p.mean(0)
    s = np.sqrt(2) / np.linalg.norm(p - c, axis=1).mean()
    return np.array([[s, 0, -s * c[0]], [0, s, -s * c[1]], [0, 0, 1]])

def dlt(p, q):
    """Homography with q ~ H p, from the last right singular vector."""
    T, Tq = normalize(p), normalize(q)
    pn = np.c_[p, np.ones(len(p))] @ T.T
    qn = np.c_[q, np.ones(len(q))] @ Tq.T
    rows = []
    for (x, y, w), (u, v, t) in zip(pn, qn):
        rows.append([0, 0, 0, -t * x, -t * y, -t * w, v * x, v * y, v * w])
        rows.append([t * x, t * y, t * w, 0, 0, 0, -u * x, -u * y, -u * w])
    A = np.array(rows)
    _, s, Vt = np.linalg.svd(A)
    Hn = Vt[-1].reshape(3, 3)                  # last ROW of Vt = last column of V
    H = np.linalg.inv(Tq) @ Hn @ T
    return H / H[2, 2], s, A, Vt[-1]

def corner_error(H):
    corners = np.array([[0, 0], [640, 0], [640, 480], [0, 480]], float)
    return np.abs(project(H, corners) - project(H_true, corners)).max()

rng = np.random.default_rng(1)
p = rng.uniform([0, 0], [640, 480], size=(50, 2))          # 50 points in a 640 x 480 image
q_exact = project(H_true, p)

H0, s0, _, _ = dlt(p, q_exact)
print("noise-free: s_8 = %.3g, s_9 = %.3g, corner error %.1e px" % (s0[7], s0[8], corner_error(H0)))

q = q_exact + rng.normal(0, 1.0, q_exact.shape)            # 1-pixel Gaussian noise
H1, s1, A, h = dlt(p, q)
print("noisy:      s_8 = %.3g, s_9 = %.3g, corner error %.2f px" % (s1[7], s1[8], corner_error(H1)))
print("||A h|| = %.4g = s_9;  min eigenvalue of A^T A = %.4g = s_9^2"
      % (np.linalg.norm(A @ h), np.linalg.eigvalsh(A.T @ A)[0]))

# Spread over noise draws: the same estimator on 200 noise realizations.
errs = [corner_error(dlt(p, q_exact + np.random.default_rng(k).normal(0, 1.0, q_exact.shape))[0])
        for k in range(200)]
print("corner error over 200 noise draws: median %.2f px, 90th percentile %.2f px"
      % (np.median(errs), np.percentile(errs, 90)))

Output:

noise-free: s_8 = 4.74, s_9 = 1.99e-15, corner error 2.8e-13 px
noisy:      s_8 = 4.74, s_9 = 0.0465, corner error 0.92 px
||A h|| = 0.0465 = s_9;  min eigenvalue of A^T A = 0.002162 = s_9^2
corner error over 200 noise draws: median 1.01 px, 90th percentile 1.49 px

With exact correspondences, the smallest singular value is at rounding level and the estimate is exact: the $100 \times 9$ matrix has a one-dimensional null space and the last right singular vector spans it. Noise lifts $\sigma_9$ off zero but leaves it far below $\sigma_8$, so the minimizer is still well determined, and the gap between them is what keeps the estimate stable. The residual $\lVert A\hat{\mathbf{h}} \rVert$ equals $\sigma_9$, and the smallest eigenvalue of $A^\top A$ equals $\sigma_9^2$, the two equivalent statements of the solution. Over 200 noise draws the corners move by about 1 pixel at the median and 1.5 pixels at the 90th percentile, with the points normalized; the homography article shows how much worse it gets without.

The pseudoinverse of a downsampling operator

Downsampling an 8-pixel image row by averaging pairs of pixels is a wide $4 \times 8$ matrix: four equations in eight unknowns. Recovering the sharp row from the low-resolution one is underdetermined, and the pseudoinverse picks the solution of smallest norm.

import numpy as np
np.set_printoptions(precision=3, suppress=True)

# 2x downsampling of one 8-pixel image row by averaging neighbouring pairs: a wide 4 x 8 matrix.
A = np.kron(np.eye(4), [0.5, 0.5])
x_true = np.array([0.1, 0.3, 0.8, 0.8, 0.2, 0.6, 0.4, 0.4])   # the sharp row
b = A @ x_true                                                   # what the sensor records

U, s, Vt = np.linalg.svd(A)                     # full SVD: Vt is 8 x 8
print("singular values:", s, "  rank:", np.linalg.matrix_rank(A))
x_pinv = Vt[:4].T @ ((U.T @ b) / s)             # V Sigma^+ U^T b, using the 4 non-zero values
print("min-norm solution:  ", x_pinv)
print("np.linalg.pinv agrees:", np.allclose(np.linalg.pinv(A) @ b, x_pinv),
      "  lstsq agrees:", np.allclose(np.linalg.lstsq(A, b)[0], x_pinv))

# Every exact solution is x_pinv plus a null-space vector; the null space is
# spanned by the last 4 right singular vectors: pixel pairs that differ by +d/-d.
print("null-space basis (rows of Vt[4:]):\n", Vt[4:])
print("true row = x_pinv + null-space part:", np.allclose(Vt[4:].T @ (Vt[4:] @ x_true) + x_pinv, x_true))
print("norms: ||x_pinv|| = %.3f, ||x_true|| = %.3f" % (np.linalg.norm(x_pinv), np.linalg.norm(x_true)))

Output:

singular values: [0.707 0.707 0.707 0.707]   rank: 4
min-norm solution:   [0.2 0.2 0.8 0.8 0.4 0.4 0.4 0.4]
np.linalg.pinv agrees: True   lstsq agrees: True
null-space basis (rows of Vt[4:]):
 [[-0.354  0.354 -0.354  0.354  0.5   -0.5    0.     0.   ]
 [-0.354  0.354 -0.354  0.354 -0.5    0.5    0.     0.   ]
 [-0.354  0.354  0.354 -0.354  0.     0.     0.5   -0.5  ]
 [-0.354  0.354  0.354 -0.354  0.     0.    -0.5    0.5  ]]
true row = x_pinv + null-space part: True
norms: ||x_pinv|| = 1.414, ||x_true|| = 1.449

All four singular values equal $1/\sqrt{2}$, so the operator is perfectly conditioned on its row space and loses information only through its null space. The minimum-norm solution replicates each low-resolution pixel into both of its high-resolution pixels, the familiar nearest-neighbour upsampling, and pinv and lstsq return the same vector. The null space is four-dimensional and consists of rows in which the two pixels of each pair differ by $+d$ and $-d$: detail that averaging erases completely. Because the zero singular value is repeated, LAPACK returns one orthonormal basis of that space among many, here mixtures of the four pairs rather than one pair per vector. The true row is the minimum-norm solution plus a null-space component, and it has a larger norm, as every other exact solution must. Choosing the true row requires knowledge that the data do not contain, which is the prior that super-resolution methods supply.

Projecting onto essential matrices and rotations

The first part estimates an essential matrix from 100 noisy matches between two calibrated views by the linear eight-point method, then projects it to rank 2 and onto the essential matrices. A matrix with singular values $(\sigma, \sigma, 0)$ satisfies $2EE^\top E - \operatorname{tr}(EE^\top)E = 0$, since both terms equal $2\sigma^3$ times the same rank-2 matrix; the code uses this cubic identity as a test. The second part aligns 30 points on an object with a rotated, noisy observation of them, once by fitting an unconstrained $3 \times 3$ matrix and projecting it to a rotation, and once by solving the Procrustes problem directly, and repeats the comparison over 200 noise draws.

import numpy as np
np.set_printoptions(precision=4, suppress=True)
rng = np.random.default_rng(3)

def rot(axis, deg):
    a = np.asarray(axis, float) / np.linalg.norm(axis); t = np.deg2rad(deg)
    K = np.array([[0, -a[2], a[1]], [a[2], 0, -a[0]], [-a[1], a[0], 0]])
    return np.eye(3) + np.sin(t) * K + (1 - np.cos(t)) * K @ K

def angle_deg(R1, R2):
    return np.degrees(np.arccos(np.clip((np.trace(R1.T @ R2) - 1) / 2, -1, 1)))

# --- Part 1: a linear eight-point estimate of E, then the nearest valid matrices.
R, t = rot([0.2, 1, 0.1], 12), np.array([1.0, 0.1, 0.05])
X = rng.uniform([-2, -2, 4], [2, 2, 8], size=(100, 3))          # scene points
x1 = X / X[:, 2:]                                                 # normalized image coords
X2 = X @ R.T + t
x2 = X2 / X2[:, 2:]
f = 800.0                                                         # 1-pixel noise at f = 800
x1[:, :2] += rng.normal(0, 1 / f, (100, 2)); x2[:, :2] += rng.normal(0, 1 / f, (100, 2))
A = np.einsum("ni,nj->nij", x2, x1).reshape(-1, 9)              # rows of x2^T E x1 = 0
E_lin = np.linalg.svd(A)[2][-1].reshape(3, 3)
U, s, Vt = np.linalg.svd(E_lin)
print("singular values of the linear estimate (scaled):", s / s[0])

E_rank2 = U @ np.diag([s[0], s[1], 0]) @ Vt                     # nearest rank-2 matrix
E_ess = U @ np.diag([1, 1, 0]) @ Vt                              # nearest essential, up to scale
for name, E in [("linear", E_lin), ("rank-2", E_rank2), ("essential", E_ess)]:
    En = E / np.linalg.norm(E)
    sv = np.linalg.svd(En, compute_uv=False)
    cubic = np.linalg.norm(2 * En @ En.T @ En - np.trace(En @ En.T) * En)
    print(f"{name:9s} sigma = {sv}  |2EE^T E - tr(EE^T)E| = {cubic:.1e}")
print("Frobenius distance to the rank-2 matrix = s_3:", np.isclose(np.linalg.norm(E_lin - E_rank2), s[2]))

# --- Part 2: nearest rotation (orthogonal Procrustes) for 3D points seen twice.
P = rng.uniform(-0.5, 0.5, size=(30, 3))                        # points on an object, metres
R_true = rot([1, 2, 3], 40)
Q = P @ R_true.T + rng.normal(0, 0.02, P.shape)                 # 2 cm noise; centred, no translation

M, *_ = np.linalg.lstsq(P, Q, rcond=None)                       # unconstrained linear fit Q ~ P M
M = M.T
print("linear fit: det = %.3f, max |M^T M - I| = %.3f" % (np.linalg.det(M), np.abs(M.T @ M - np.eye(3)).max()))

def nearest_rotation(M):
    U, _, Vt = np.linalg.svd(M)
    D = np.diag([1, 1, np.sign(np.linalg.det(U @ Vt))])       # forbid a reflection
    return U @ D @ Vt

R_proj = nearest_rotation(M)                                    # fit, then project
R_proc = nearest_rotation(Q.T @ P)                              # Procrustes: maximize tr(R^T Q^T P)
for name, Rh in [("projected linear fit", R_proj), ("Procrustes", R_proc)]:
    print(f"{name:21s} det = {np.linalg.det(Rh):.6f}  error {angle_deg(Rh, R_true):.3f} deg"
          f"  sum of squared residuals {np.sum((Q - P @ Rh.T) ** 2):.4f}")

# The same comparison over 200 noise draws (angle error in degrees).
errs = []
for k in range(200):
    g = np.random.default_rng(100 + k)
    Qk = P @ R_true.T + g.normal(0, 0.02, P.shape)
    Mk = np.linalg.lstsq(P, Qk, rcond=None)[0].T
    errs.append([angle_deg(nearest_rotation(Mk), R_true), angle_deg(nearest_rotation(Qk.T @ P), R_true)])
errs = np.array(errs)
print("median error over 200 draws: projected fit %.3f deg, Procrustes %.3f deg; Procrustes smaller in %d of 200"
      % (*np.median(errs, 0), (errs[:, 1] < errs[:, 0]).sum()))

Output:

singular values of the linear estimate (scaled): [1.     0.9952 0.0073]
linear    sigma = [0.7088 0.7054 0.0052]  |2EE^T E - tr(EE^T)E| = 7.1e-03
rank-2    sigma = [0.7088 0.7054 0.    ]  |2EE^T E - tr(EE^T)E| = 4.8e-03
essential sigma = [0.7071 0.7071 0.    ]  |2EE^T E - tr(EE^T)E| = 2.5e-16
Frobenius distance to the rank-2 matrix = s_3: True
linear fit: det = 0.994, max |M^T M - I| = 0.033
projected linear fit  det = 1.000000  error 0.591 deg  sum of squared residuals 0.0384
Procrustes            det = 1.000000  error 0.710 deg  sum of squared residuals 0.0383
median error over 200 draws: projected fit 0.788 deg, Procrustes 0.769 deg; Procrustes smaller in 118 of 200

The linear estimate is close to an essential matrix but not one: its third singular value is not zero and its first two differ by about 0.5%, so it fails the cubic test. Zeroing $\sigma_3$ gives a rank-2 matrix at Frobenius distance exactly $\sigma_3$ that still fails the test, because rank 2 is necessary but not sufficient; replacing the singular values with $(1, 1, 0)$ satisfies it to rounding error. In the second part, the unconstrained fit is nearly but not exactly orthogonal, and both SVD routes produce proper rotations with determinant 1. Procrustes has the smaller sum of squared residuals in the shown draw, as it must, even though its rotation error happens to be the larger of the two there, since it minimizes that sum over all rotations, while fitting first and projecting afterwards optimizes a different criterion. The median rotation errors of the two over 200 draws differ by about two hundredths of a degree, and Procrustes is the more accurate in 118 of them: the practical difference is small at this noise level, but only Procrustes is the least-squares rotation.

The rank of a tracked image sequence

Sixty points on the faces of a box are tracked through 20 orthographic frames while the box turns by 1.5 degrees per frame, 28.5 degrees in all, and the tracks are perturbed by 0.5-pixel noise. After the centroid of each frame is subtracted, the $40 \times 60$ measurement matrix should have rank 3.

import numpy as np
np.set_printoptions(precision=2, suppress=True)
rng = np.random.default_rng(0)

def rot(axis, deg):
    a = np.asarray(axis, float) / np.linalg.norm(axis); t = np.deg2rad(deg)
    K = np.array([[0, -a[2], a[1]], [a[2], 0, -a[0]], [-a[1], a[0], 0]])
    return np.eye(3) + np.sin(t) * K + (1 - np.cos(t)) * K @ K

# 60 feature points on the faces of a box, tracked through 20 orthographic frames
# while the box turns by 1.5 degrees per frame; image coordinates in pixels.
P = rng.uniform(-1, 1, size=(60, 3)) * [100, 60, 40]
face = rng.integers(0, 3, 60)
P[np.arange(60), face] = np.sign(P[np.arange(60), face]) * np.array([100, 60, 40])[face]
F = 20
W = np.vstack([(P @ rot([0.3, 1, 0.2], 1.5 * f).T)[:, :2].T + rng.uniform(-50, 50, (2, 1))
               for f in range(F)])                                # 2F x P: x and y rows per frame
W_noisy = W + rng.normal(0, 0.5, W.shape)                         # 0.5-pixel tracking noise

for name, M in [("noise-free", W), ("noisy", W_noisy)]:
    Wc = M - M.mean(axis=1, keepdims=True)                        # register: subtract centroids
    U, s, Vt = np.linalg.svd(Wc, full_matrices=False)
    W3 = (U[:, :3] * s[:3]) @ Vt[:3]                              # best rank-3 approximation
    rms = np.sqrt(np.mean((W3 - Wc) ** 2))
    print(f"{name:10s} s_1..s_6 = {s[:6]}  RMS residual of rank 3 = {rms:.3f} px")

Output:

noise-free s_1..s_6 = [2556.02 1408.34  169.9     0.      0.      0.  ]  RMS residual of rank 3 = 0.000 px
noisy      s_1..s_6 = [2555.89 1408.51  169.54    6.52    6.19    5.97]  RMS residual of rank 3 = 0.462 px

Without noise, exactly three singular values are non-zero, Tomasi and Kanade’s rank theorem. With noise, all forty become non-zero, but the fourth and later ones form a floor about 26 times below the third, set by the tracking noise, and the gap identifies the rank unambiguously. The third singular value is much smaller than the first two because the box turns through less than 30 degrees, so depth is seen through small changes in the projections; a narrower range of motion would shrink $\sigma_3$ toward the noise floor (our reasoning), which is why the authors proposed testing the ratio of the third to the fourth singular value. The best rank-3 approximation leaves an RMS residual a little below the 0.5-pixel noise, because the three retained components absorb a small fraction of the noise; the factorization then splits this rank-3 matrix into camera axes and shape.

Common Misconceptions

• “Vh holds the right singular vectors in its columns.” NumPy and SciPy return $V^\top$. The right singular vectors are the rows of Vh, and the last right singular vector, the DLT solution, is Vh[-1], not Vh[:, -1].
• “The SVD is the eigendecomposition.” Only for symmetric positive semidefinite matrices. In general the eigendecomposition exists only for some square matrices, its eigenvalues can be negative or complex and its eigenvectors need not be orthogonal, whereas the SVD exists for every matrix, uses two different orthonormal bases and has non-negative singular values (Deisenroth, Faisal and Ong, 2020, §4.5.3).
• “The Eckart–Young theorem holds only in the Frobenius norm.” Truncation is optimal in the spectral norm and in every unitarily invariant norm (Mirsky, as reported by Stewart, 1992). What the theorem does not cover is a weighted or masked error, such as one that ignores the missing entries of a measurement matrix in which tracked points are occluded; Tomasi and Kanade (1992) handled that case with a separate iterative procedure.
• “Projection by SVD gives the best constrained estimate.” It gives the nearest valid matrix in the Frobenius norm of the matrix entries, not the most likely one given the image measurements; it is an initialization for refinement of a geometric error.

Where It Is Used

• Two-view and multi-view geometry. The fundamental matrix (eight-point algorithm and rank-2 correction), the essential matrix (projection and pose recovery), the homography (DLT), triangulation, perspective-n-point, and the factorization methods of multi-view geometry and structure from motion.
• Calibration and pose. Camera calibration solves homogeneous systems for homographies and intrinsics and projects estimated rotations onto rotations; rotation representations, iterative closest point and absolute trajectory error use the Procrustes solution.
• Numerical linear algebra. Rank and null-space computation, condition numbers, pseudoinverses, and the SVD-based least squares solver behind np.linalg.lstsq.
• Data analysis and learning. Principal component analysis is the SVD of a centred data matrix, applied to face images, image patches and descriptors; low-rank approximation compresses and denoises image and feature matrices.