Linear Algebra

Linear algebra is the mathematics of vectors, matrices and the linear maps between them. Computer vision is written in it: a digital image is an array that can be read as a vector or a matrix, a data set of images or descriptors is a matrix with one sample per row, colour conversions, blurs and the layers of a neural network are matrix products, and the geometry of cameras and views is expressed as $3 \times 3$ and $3 \times 4$ matrices acting on homogeneous coordinates. Estimating a camera, a homography or a fundamental matrix comes down to solving a linear system, and whether the answer can be trusted depends on the rank and conditioning of that system.

This article is an overview of the parts of linear algebra that vision uses, and it links the articles that develop each part. Least squares covers overdetermined systems and their numerical solution. Two factorizations do much of the practical work: the singular value decomposition, which writes any matrix as $A = U \Sigma V^\top$ with orthogonal $U$ and $V$ and non-negative singular values on the diagonal of $\Sigma$, and which underlies rank estimation, low-rank approximation and the solution of homogeneous systems; and the eigendecomposition of symmetric matrices, which underlies principal component analysis.

Definition

Linear algebra studies vector spaces, sets whose elements can be added and scaled, and linear maps between them, functions $f$ with $f(a\mathbf{x} + b\mathbf{y}) = a f(\mathbf{x}) + b f(\mathbf{y})$. In finite dimensions every vector is a list of $n$ numbers, every linear map from $\mathbb{R}^n$ to $\mathbb{R}^m$ is multiplication by an $m \times n$ matrix, and the subject becomes the study of matrices: how they combine, which vectors they send to zero, which outputs they can reach, how they change lengths and angles, and how to undo them.

Intuition

Think of a grayscale image of $H \times W$ pixels as one point in a space with $HW$ axes, one per pixel. Brightening the image scales the point; blending two images is a weighted sum of two points; the difference between two images is an arrow from one point to the other, and its length measures how different they are. Many image operations are linear in this sense: a blur, a colour conversion, a derivative filter or a downsampling replaces each output pixel with a fixed weighted sum of input pixels, so each is a matrix, however large.

Geometry adds a second picture. A $2 \times 2$ or $3 \times 3$ matrix moves points of the plane or of space: it rotates, scales, shears or projects them. The columns of the matrix are where the coordinate axes go, so a matrix can be read directly as a description of what it does to the grid (Boyd and Vandenberghe, 2018, §7.1).

Both pictures lead to the same questions. Can the map be undone, or does it lose information? Which inputs does it ignore entirely? If the outputs are known, can the inputs be recovered, and how sensitive is the answer to noise in the outputs? Rank, null space, determinant and condition number are the answers.

Formal Definition

Notation follows the rest of the site: bold lowercase letters are column vectors, uppercase letters are matrices, $A^\top$ is the transpose, and a data set is $X \in \mathbb{R}^{n \times d}$ with one sample per row, the layout scikit-learn uses. Textbooks often store samples as columns instead (Boyd and Vandenberghe’s feature matrix in §8.2 is $d \times n$), so the transposes in formulas from different sources differ.

Vectors, matrices and images

An $M \times N$ grayscale image is a vector of length $MN$ once its pixels are listed in an agreed order, row by row or column by column; a colour image is a vector of length $3MN$, and a video of $K$ frames one of length $KMN$ (Boyd and Vandenberghe, 2018, §1.1). The same image is also naturally an $M \times N$ matrix, with the row index giving the vertical position and the column index the horizontal one (§6.1). NumPy’s reshape switches between the two views without moving any data, and the digital image article covers the indexing and memory-layout conventions.

Stacking samples gives a data matrix. The pixels of a colour image form an $HW \times 3$ matrix, one pixel per row; a collection of $n$ images of $d$ pixels each forms an $n \times d$ matrix, one image per row; $n$ descriptors such as SIFT vectors form an $n \times 128$ matrix. Statistics of the data are then matrix expressions. With the column means subtracted to give $X_c$, the empirical covariance matrix is

$$C = \frac{1}{n} X_c^\top X_c \in \mathbb{R}^{d \times d},$$

with denominator $n$; the unbiased estimate divides by $n - 1$, which is NumPy’s np.cov default. A covariance matrix is symmetric and positive semidefinite (Deisenroth, Faisal and Ong, 2020, §6.4.2).

Matrix products as linear maps

The product $\mathbf{y} = A\mathbf{x}$ of an $m \times n$ matrix and an $n$-vector is a linear map from $\mathbb{R}^n$ to $\mathbb{R}^m$, and every such map has this form. The $i$-th column of $A$ is the image of the $i$-th unit vector, which gives a quick way to write down the matrix of any linear transformation (Boyd and Vandenberghe, 2018, §7.1). Composing maps multiplies matrices: applying $A$ and then $B$ is the single map $BA$, and the order matters because matrix multiplication is not commutative.

Three families of maps recur in vision.

• Per-pixel maps. A colour conversion applies the same $3 \times 3$ matrix to every pixel’s colour vector. For a data matrix with one pixel per row, the whole image is converted by one product, $Y = X M^\top$. The color spaces article describes the conversions; linear sRGB to CIE XYZ is purely linear, while RGB to YCbCr in JPEG adds an offset and is affine.
• Spatial maps. A blur or a derivative filter computes each output pixel as a weighted sum of nearby input pixels, so with the image written as a vector it is a large, sparse matrix. For 2D convolution of an image with a kernel, $\mathbf{y} = T\mathbf{x}$ for a structured matrix $T$ built from the kernel’s entries (Boyd and Vandenberghe, 2018, §7.4). Libraries never form $T$, but the view explains why blurring is a linear inverse problem and why a convolutional neural network layer, being linear before its activation function, can be computed as one large matrix product.
• Geometric maps. Rotations, scalings, reflections and projections of points are $2 \times 2$ or $3 \times 3$ matrices (Boyd and Vandenberghe, 2018, §7.1). A rotation matrix $R$ is orthogonal, $R^\top R = I$, so its inverse is its transpose and it preserves lengths and angles (Deisenroth, Faisal and Ong, 2020, §3.4); rotation representations covers the ways of parameterizing it.

A translation, $\mathbf{x} \mapsto \mathbf{x} + \mathbf{t}$, is not linear, since it moves the origin. An affine map $\mathbf{x} \mapsto A\mathbf{x} + \mathbf{t}$ becomes linear in homogeneous coordinates, where a 2D point is written $(x, y, 1)^\top$ and the map is a $3 \times 3$ matrix with last row $(0, 0, 1)$. Allowing a general last row gives a projective transformation, the homography, defined only up to scale. A pinhole camera is a $3 \times 4$ matrix $P = K[R \mid \mathbf{t}]$ from homogeneous 3D points to homogeneous image points, with the intrinsic matrix $K$ and the pose $R, \mathbf{t}$ described in camera parameters. The division by the last coordinate that turns a homogeneous vector back into a point is the only nonlinear step, and keeping it outside the matrix is what lets chains of warps and projections collapse into single products. Image warping lists the $3 \times 3$ matrices from translation to homography, and projective geometry places them in a hierarchy by what each preserves.

Column space, null space and rank

Two subspaces describe what a matrix $A \in \mathbb{R}^{m \times n}$ does (Deisenroth, Faisal and Ong, 2020, §2.7.3):

• The column space (image, range) is the set of all outputs $A\mathbf{x}$, the span of the columns of $A$, a subspace of $\mathbb{R}^m$.
• The null space (kernel) is the set of inputs that $A$ sends to zero, $\{\mathbf{x} : A\mathbf{x} = \mathbf{0}\}$, a subspace of $\mathbb{R}^n$. It is the general solution of the homogeneous system $A\mathbf{x} = \mathbf{0}$.

The rank of $A$ is the number of linearly independent columns, which always equals the number of linearly independent rows, and it is the dimension of the column space. A matrix has full rank when its rank is $\min(m, n)$ and is rank deficient otherwise (§2.6.2). The rank–nullity theorem ties the two subspaces together:

$$\dim(\text{null space of } A) + \operatorname{rank}(A) = n$$

(Deisenroth, Faisal and Ong, 2020, Theorem 2.24). Every dimension of the input either reaches the output or is lost in the null space.

In vision, rank is a statement about how much independent information a matrix holds. A stack of images that are all combinations of $k$ underlying patterns has a data matrix of rank at most $k$, however many pixels and images there are; real image collections are often close to such a low-rank matrix, which principal component analysis and low-rank approximation exploit. A fundamental matrix has rank 2 by construction, and its estimators enforce that rank explicitly (fundamental matrix). The camera centre is the null vector of the $3 \times 4$ camera matrix (camera parameters).

Linear systems

A system of $m$ linear equations in $n$ unknowns is $A\mathbf{x} = \mathbf{b}$. It is square if $m = n$, overdetermined if $m > n$ (a tall $A$) and underdetermined if $m < n$ (a wide $A$); with $\mathbf{b} = \mathbf{0}$ it is homogeneous and always has the solution $\mathbf{x} = \mathbf{0}$ (Boyd and Vandenberghe, 2018, §8.3). A solution exists exactly when $\mathbf{b}$ lies in the column space of $A$, that is when appending $\mathbf{b}$ as an extra column does not raise the rank, and if one exists, every vector of the null space can be added to it, so it is unique only when the null space is $\{\mathbf{0}\}$ (Deisenroth, Faisal and Ong, 2020, §2.3.1 and §2.6.2).

• Square and full rank. There is exactly one solution, $\mathbf{x} = A^{-1}\mathbf{b}$, for every $\mathbf{b}$ (Boyd and Vandenberghe, 2018, §11.2). Small calibrated problems, such as converting a pixel to a ray with $K^{-1}$, are of this kind.
• Overdetermined. With noisy measurements, more equations than unknowns are generally inconsistent, and the usual answer minimizes $\lVert A\mathbf{x} - \mathbf{b} \rVert^2$; the least squares article covers the normal equations, the geometric view as a projection, and QR- and SVD-based solvers. Fitting a line to edge points, a translation to a window of image gradients as in Lucas–Kanade, or a camera to calibration points all take this form.
• Underdetermined. Fewer equations than unknowns leave a null space of solutions, and a choice must be made, typically the solution of smallest norm or one favoured by a regularizer. Recovering a high-resolution image from a low-resolution one is underdetermined in this sense, since the downsampling matrix is wide, which is why super-resolution needs a prior on what images look like.
• Homogeneous with a constraint. Many geometric estimators produce equations that are linear in the unknown entries of a matrix and have zero right-hand side. The direct linear transform (DLT) for a homography stacks two equations per correspondence into $A\mathbf{h} = \mathbf{0}$, with $\mathbf{h}$ the nine entries of $H$; the eight-point algorithm does the same for the fundamental matrix, one equation per match. The trivial solution is excluded by requiring $\lVert \mathbf{h} \rVert = 1$, which matches the fact that $H$ is defined only up to scale. With exact data the answer spans the one-dimensional null space of $A$; with noisy data $A$ has full column rank, and the unit vector minimizing $\lVert A\mathbf{h} \rVert$ is the right singular vector of $A$ for its smallest singular value, as the homography and least squares articles show.

Norms and inner products

The inner product (dot product) of two vectors is $\mathbf{x}^\top\mathbf{y} = \sum_i x_i y_i$, and the Euclidean norm is $\lVert \mathbf{x} \rVert = \sqrt{\mathbf{x}^\top \mathbf{x}}$. They define distance, $\lVert \mathbf{x} - \mathbf{y} \rVert$, and angle, $\cos\theta = \mathbf{x}^\top\mathbf{y} / (\lVert \mathbf{x} \rVert \lVert \mathbf{y} \rVert)$ (Deisenroth, Faisal and Ong, 2020, §3.1–§3.4). Other norms weight the coordinates differently: the $\ell_1$ norm sums absolute values and the $\ell_\infty$ norm takes the largest one. For matrices, the Frobenius norm $\lVert A \rVert_F$ treats the entries as one long vector, and the spectral norm $\lVert A \rVert_2$ is the largest factor by which $A$ stretches any vector, its largest singular value.

On images these are familiar quantities. The squared Euclidean distance between two images, divided by the number of pixels, is the mean squared error. The inner product of two mean-subtracted patches, divided by their norms, is normalized cross-correlation, the cosine of the angle between them. A linear classifier scores an input by an inner product with a weight vector. The similarity and distance measures article compares these measures and others used for descriptors and embeddings.

Orthogonality and projections

Two vectors are orthogonal when their inner product is zero. A set of mutually orthogonal unit vectors is orthonormal, and a square matrix with orthonormal columns is an orthogonal matrix, $Q^\top Q = I$, whose inverse is its transpose (Deisenroth, Faisal and Ong, 2020, §3.4). Coordinates in an orthonormal basis are just inner products, $c_i = \mathbf{q}_i^\top \mathbf{x}$, which is why orthonormal bases are the convenient ones: the columns of a rotation matrix and the principal axes of a covariance matrix are both of this kind.

A projection is a linear map $P$ with $P^2 = P$: applying it twice changes nothing. The orthogonal projection onto the column space of a matrix $B$ with independent columns is

$$P = B (B^\top B)^{-1} B^\top,$$

which sends $\mathbf{x}$ to the closest point of that subspace, leaving a residual orthogonal to every column of $B$ (Deisenroth, Faisal and Ong, 2020, §3.8). If the columns of $B$ are orthonormal it simplifies to $P = BB^\top$. Least squares is exactly this projection of the measurements onto the column space of the design matrix, as the least squares article shows. Projecting an image onto the span of a few basis images, keeping what they can express and discarding the rest, is the idea behind subspace methods for faces and appearance.

Determinants and invertibility

The determinant of a square matrix is the signed volume of the parallelepiped spanned by its columns: the factor by which the map scales areas (in 2D) or volumes, with a negative sign when it reverses orientation (Deisenroth, Faisal and Ong, 2020, §4.1). A square matrix is invertible exactly when its determinant is non-zero, equivalently when it has full rank (Theorems 4.1 and 4.3). The determinant is multiplicative, $\det(AB) = \det A \det B$, and equals the product of the eigenvalues (Theorem 4.16). A rotation has determinant $+1$, a reflection $-1$, an affine image warp scales every area by $|\det A|$ for its $2 \times 2$ block, and the rank-2 constraint on a fundamental matrix is $\det F = 0$.

As a numerical test of invertibility the determinant is nearly useless, and Deisenroth, Faisal and Ong note that direct numerical methods have replaced its explicit computation (§4.1). The reason is scale: for the $100 \times 100$ matrix $0.1 I$, which is perfectly well behaved, the determinant is $10^{-100}$, while a matrix with determinant 1 can be close to singular. Whether a system can be solved reliably is measured by the condition number below, and whether a matrix is singular in practice by its rank with a tolerance.

Eigenvalues and eigenvectors

An eigenvector of a square matrix $A$ is a non-zero vector that $A$ only scales, $A\mathbf{v} = \lambda\mathbf{v}$; the scale $\lambda$ is its eigenvalue. Eigenvectors are the directions a linear map leaves unturned. A general matrix may have complex eigenvalues or too few independent eigenvectors, but symmetric matrices, which are the ones vision meets most, are well behaved: by the spectral theorem, a real symmetric matrix has real eigenvalues and an orthonormal basis of eigenvectors, so $A = Q \Lambda Q^\top$ with $Q$ orthogonal and $\Lambda$ diagonal (Deisenroth, Faisal and Ong, 2020, Theorem 4.15). Any matrix of the form $A^\top A$ is symmetric and positive semidefinite, with non-negative eigenvalues (Theorem 4.14). This factorization is the eigendecomposition.

Two symmetric positive semidefinite matrices recur throughout vision:

• Covariance matrices. The eigenvectors of the covariance of a data matrix are the principal axes of the data and the eigenvalues are the variances along them; principal component analysis keeps the axes with the largest eigenvalues. For image patches or face images this gives a low-dimensional basis, and for the pixel colours of a region it gives the dominant colour direction.
• The structure tensor. The $2 \times 2$ matrix of summed gradient products over a window has eigenvalues that classify the window as flat, edge or corner. Both are small in a flat region, one is large on an edge, and both are large at a corner; the Harris corner detector scores corners from the determinant and trace, and the aperture problem is the case of one zero eigenvalue.

The singular value decomposition extends the spectral theorem to every rectangular matrix: the singular values of $A$ are the square roots of the eigenvalues of $A^\top A$, and the right singular vectors are its eigenvectors (Deisenroth, Faisal and Ong, 2020, §4.5).

Properties

Conditioning

A system can be exactly solvable and still useless if a small change in the data causes a large change in the answer. For a square invertible $A$, perturb the right-hand side of $A\mathbf{x} = \mathbf{b}$ to $\mathbf{b} + \delta\mathbf{b}$. The solution moves by $\delta\mathbf{x} = A^{-1}\delta\mathbf{b}$, so $\lVert \delta\mathbf{x} \rVert \le \lVert A^{-1} \rVert \lVert \delta\mathbf{b} \rVert$, and $\lVert \mathbf{b} \rVert \le \lVert A \rVert \lVert \mathbf{x} \rVert$. Dividing gives the standard bound

$$\frac{\lVert \delta\mathbf{x} \rVert}{\lVert \mathbf{x} \rVert} \le \kappa(A) \, \frac{\lVert \delta\mathbf{b} \rVert}{\lVert \mathbf{b} \rVert}, \qquad \kappa(A) = \lVert A \rVert \, \lVert A^{-1} \rVert = \frac{\sigma_{\max}}{\sigma_{\min}},$$

where the last form holds for the spectral norm and uses the largest and smallest singular values; np.linalg.cond computes it this way by default. The condition number is at least 1, equals 1 for orthogonal matrices, and is infinite for singular ones. It is a property of the problem, not of the algorithm: for the worst perturbations the bound is attained, so no method can guarantee better relative accuracy than $\kappa(A)$ times the relative error in the data, and rounding errors alone, about $10^{-16}$ in double precision, are data errors of this kind.

In vision, ill-conditioning has recognizable sources. Undoing a blur is ill conditioned because blurring nearly erases fine detail, so the matrix has tiny singular values for high-frequency patterns. A homography or fundamental-matrix system built from raw pixel coordinates is badly conditioned because its columns differ in scale by orders of magnitude; translating and scaling the coordinates first, as the normalized eight-point paper argued, fixes most of it, and the homography article measures the effect. Forming $A^\top A$ squares the condition number, the reason least-squares solvers avoid the normal equations on hard problems.

Solve, do not invert

Mathematically $\mathbf{x} = A^{-1}\mathbf{b}$, but numerical code almost never forms $A^{-1}$. A library solver factors $A$ (LU with pivoting, QR or Cholesky) and solves the triangular systems that result. Boyd and Vandenberghe count about $2n^3$ floating-point operations to solve a square system through a QR factorization and about $3n^3$ to compute the inverse through the same factorization, and a factorization can be reused for many right-hand sides at little extra cost (2018, §11.3). Deisenroth, Faisal and Ong add that computing the inverse is generally not recommended for reasons of numerical precision (2020, §2.3.4).

The accuracy argument is subtler than often stated. Druinsky and Toledo (2012) showed that if the computed inverse $V$ is a good left inverse, with $\lVert VA - I \rVert$ of order $\kappa(A)$ times machine precision, then $V\mathbf{b}$ is as accurate as the solution of a backward-stable solver; but the residual $A\mathbf{x} - \mathbf{b}$ of $V\mathbf{b}$ is not small in general, so the computed answer does not carry the usual certificate of a small backward error, and an inverse computed one way can be a good right inverse without being a good left inverse. Solving directly with a backward-stable factorization gives both the accuracy and the small residual without depending on how an inverse was computed, costs less, and is what np.linalg.solve, scipy.linalg.solve and the backslash operator of other systems do.

Rank in floating point

Exact rank deficiency does not survive rounding or noise: the columns of a measured matrix are almost never exactly dependent, so in floating point nearly every matrix has full rank. Numerical rank is therefore defined with a tolerance, as the number of singular values above a threshold. NumPy’s matrix_rank uses $\sigma_{\max} \cdot \max(m, n) \cdot \varepsilon$ by default, with $\varepsilon$ the machine precision, and its documentation advises a tolerance near the measurement uncertainty when the data are noisier than that (NumPy 2.3 documentation). The same issue appears whenever an algorithm tests for singularity: a matrix whose smallest singular value is a millionth of its largest is singular for any practical purpose, whatever its determinant.

Examples

The examples use NumPy 2.3.5 and SciPy 1.16.3 on synthetic images.

An image as a data matrix, and colour conversion as a matrix

The JPEG File Interchange Format (JFIF) converts 8-bit RGB to YCbCr by fixed weights plus an offset of 128 on the two chroma channels, and gives the inverse coefficients as well (Hamilton, 1992). Written with one pixel per row, the conversion of a whole image is one matrix product.

import numpy as np
np.set_printoptions(suppress=True)

# A 4 x 6 synthetic RGB image: left half orange, right half sky blue, one white pixel.
img = np.zeros((4, 6, 3))
img[:, :3] = [230, 120, 30]
img[:, 3:] = [90, 160, 230]
img[0, 0] = [255, 255, 255]

H, W, _ = img.shape
X = img.reshape(-1, 3)                 # data matrix: one pixel per row, n = H*W
print("image", img.shape, "-> data matrix", X.shape)

# JFIF RGB -> YCbCr: a 3x3 matrix plus an offset, i.e. an affine map.
M = np.array([[ 0.299,   0.587,   0.114 ],
              [-0.1687, -0.3313,  0.5   ],
              [ 0.5,    -0.4187, -0.0813]])
c = np.array([0.0, 128.0, 128.0])

Y = X @ M.T + c                        # every pixel converted by one matrix product
ycc = Y.reshape(H, W, 3)
print("YCbCr of the white, orange and blue pixels:")
print(np.round(ycc[0, 0], 1), np.round(ycc[1, 0], 1), np.round(ycc[1, 5], 1))

# The same map as one 3x4 matrix acting on homogeneous pixels [R, G, B, 1].
A = np.hstack([M, c[:, None]])
Xh = np.hstack([X, np.ones((X.shape[0], 1))])
print("homogeneous form agrees:", np.allclose(Xh @ A.T, Y))

# The map is invertible: det(M) != 0, and inv(M) reproduces JFIF's inverse coefficients.
print("det(M) =", round(np.linalg.det(M), 4))
print("inv(M) =\n", np.round(np.linalg.inv(M), 4))
X_back = np.linalg.solve(M, (Y - c).T).T
print("max round-trip error:", np.abs(X_back - X).max())

Output:

image (4, 6, 3) -> data matrix (24, 3)
YCbCr of the white, orange and blue pixels:
[255. 128. 128.] [142.6  64.4 190.3] [147.  174.8  87.3]
homogeneous form agrees: True
det(M) = 0.2363
inv(M) =
 [[ 1.     -0.      1.402 ]
 [ 1.     -0.3441 -0.7141]
 [ 1.      1.772  -0.0001]]
max round-trip error: 4.263256414560601e-14

White maps to chroma $(128, 128)$, no colour, because each chroma row of $M$ sums to zero, so any grey $(g, g, g)$ has zero chroma before the offset. The orange pixel has a red-difference value above 128 and a blue-difference value below it, and the sky blue pixel the reverse. Appending a 1 to each pixel turns the affine map into a single $3 \times 4$ matrix, the same device that homogeneous coordinates use for geometry. The determinant is non-zero, so the conversion loses nothing, and the inverse matrix reproduces JFIF’s published inverse coefficients 1.402, 0.34414, 0.71414 and 1.772 to four decimals; the $-0.0001$ where JFIF has 0 comes from the forward coefficients being rounded to four decimals in the standard. The round trip uses solve rather than the inverse, and returns the pixels to within $10^{-13}$.

Rank of an image collection

Two hundred $32 \times 32$ images are generated as random mixtures of three patterns: a horizontal ramp, a vertical ramp and a bright disk. As a $200 \times 1024$ data matrix, one image per row, the collection has the rank of the patterns, not of its size.

import numpy as np

rng = np.random.default_rng(0)
s = 32
yy, xx = np.mgrid[0:s, 0:s] / (s - 1)

# Three basis patterns: horizontal ramp, vertical ramp, a bright disk.
ramp_x = xx
ramp_y = yy
disk = ((xx - 0.5) ** 2 + (yy - 0.5) ** 2 < 0.25 ** 2).astype(float)
B = np.stack([ramp_x, ramp_y, disk]).reshape(3, -1)       # 3 x 1024

# 200 images, each a random mix of the three patterns: X = C B, one image per row.
C = rng.uniform(0, 1, size=(200, 3))
X = C @ B                                                  # 200 x 1024
print("data matrix:", X.shape, " rank:", np.linalg.matrix_rank(X))

# Null space: any image orthogonal to all three patterns is invisible to X.
null_dim = X.shape[1] - np.linalg.matrix_rank(X)
print("null space dimension:", null_dim)

# A new mixed image lies in the row space; an unrelated pattern does not.
def residual_outside(img, basis):
    coef, *_ = np.linalg.lstsq(basis.T, img, rcond=None)   # project onto span of rows
    return np.linalg.norm(img - basis.T @ coef) / np.linalg.norm(img)

new_mix = (0.3 * ramp_x + 0.9 * disk).ravel()
checker = ((np.floor(xx * 4) + np.floor(yy * 4)) % 2).ravel()
print(f"relative residual outside span: mix {residual_outside(new_mix, B):.1e}, "
      f"checkerboard {residual_outside(checker, B):.2f}")

# Sensor noise makes the matrix full rank in floating point, but not "really".
Xn = X + rng.normal(0, 0.01, size=X.shape)
print("noisy rank, default tolerance:", np.linalg.matrix_rank(Xn))
print("largest singular values:", np.round(np.linalg.svd(Xn, compute_uv=False)[:5], 2))
print("noisy rank, tolerance 1.0:", np.linalg.matrix_rank(Xn, tol=1.0))

Output:

data matrix: (200, 1024)  rank: 3
null space dimension: 1021
relative residual outside span: mix 1.0e-15, checkerboard 0.75
noisy rank, default tolerance: 200
largest singular values: [329.43  53.05  38.36   0.46   0.45]
noisy rank, tolerance 1.0: 3

The rank is 3 because every row is a combination of three patterns, and rank–nullity leaves $1024 - 3 = 1021$ dimensions of images that no combination of the rows can produce. A new mixture lies in the three-dimensional row space to rounding error, while three quarters of the checkerboard’s norm lies outside it: projection onto a subspace separates what a model can express from what it cannot. Adding Gaussian noise with a standard deviation of 0.01 makes the matrix full rank, 200, under NumPy’s default tolerance, yet the singular values show three large values and then a floor near 0.46 set by the noise. A tolerance placed in that gap recovers the rank that the data really have, which is the practical meaning of numerical rank.

Geometric maps and a homogeneous system

The corners of a $640 \times 480$ image are mapped by an affine matrix and by a homography in homogeneous coordinates. The homography is then recovered from the four correspondences as the null space of the DLT system.

import numpy as np
from scipy.linalg import null_space
np.set_printoptions(suppress=True, precision=4)

# Corners of a 640 x 480 image, one homogeneous point (x, y, 1) per column.
P = np.array([[0, 640, 640, 0],
              [0, 0, 480, 480],
              [1, 1, 1, 1]], dtype=float)

def apply(T, P):
    Q = T @ P
    return Q[:2] / Q[2]                       # back to Cartesian: divide by the last row

def side_directions(Q):
    top, bottom = Q[:, 1] - Q[:, 0], Q[:, 2] - Q[:, 3]
    return top / np.linalg.norm(top), bottom / np.linalg.norm(bottom)

th = np.deg2rad(10)
A = np.array([[1.2 * np.cos(th), -np.sin(th), 50],     # affine: last row (0, 0, 1)
              [0.9 * np.sin(th), np.cos(th), 20],
              [0, 0, 1]])
H = np.array([[0.9, 0.05, 30],                           # homography: general last row
              [0.02, 1.1, 10],
              [0.0004, 0.0002, 1]])

for name, T in [("affine", A), ("homography", H)]:
    t, b = side_directions(apply(T, P))
    print(f"{name:10s} top/bottom edges parallel: {np.isclose(abs(t @ b), 1)}"
          f"  det = {np.linalg.det(T):.3f}")

# Composition and inversion are matrix products and inverses.
print("H^-1 H P == P:", np.allclose(apply(np.linalg.inv(H) @ H, P), P[:2]))

# Recover H from the four correspondences: each gives two rows of A h = 0 (the DLT).
Q = apply(H, P)
rows = []
for (x, y), (u, v) in zip(P[:2].T, Q.T):
    rows.append([-x, -y, -1, 0, 0, 0, u * x, u * y, u])
    rows.append([0, 0, 0, -x, -y, -1, v * x, v * y, v])
M = np.array(rows)                                       # 8 x 9
N = null_space(M)                                        # orthonormal basis of {h : M h = 0}
print("M:", M.shape, " rank:", np.linalg.matrix_rank(M), " null space dim:", N.shape[1])
H_est = N[:, 0].reshape(3, 3)
print("recovered H (scaled so H[2,2] = 1):\n", H_est / H_est[2, 2])

Output:

affine     top/bottom edges parallel: True  det = 1.191
homography top/bottom edges parallel: False  det = 0.974
H^-1 H P == P: True
M: (8, 9)  rank: 8  null space dim: 1
recovered H (scaled so H[2,2] = 1):
 [[ 0.9     0.05   30.    ]
 [ 0.02    1.1    10.    ]
 [ 0.0004  0.0002  1.    ]]

The affine map keeps the image’s top and bottom edges parallel, as every affine map does, while the homography’s non-zero bottom row makes them converge, which is how a photographed rectangle becomes a trapezoid. Both determinants are non-zero, so both maps can be undone by their inverse matrices. Each correspondence contributes two equations that are linear in the nine entries of $H$, because the unknown scale of $H\mathbf{x}$ cancels when $u$ and $v$ are cross-multiplied. Four correspondences give eight equations in nine unknowns, a rank-8 matrix whose null space is one-dimensional, and that null space is $H$ up to scale. SciPy’s null_space computes it from the singular value decomposition. With more than four noisy correspondences the matrix has full column rank and null_space returns no vectors; the estimate is then the singular vector of the smallest singular value, as described for the homography.

Ill-conditioning: undoing a blur

A row of 64 pixels with two bright bars is blurred by a Gaussian of standard deviation 1.5 pixels, written as a $64 \times 64$ matrix whose rows are normalized weights. Recovering the sharp row is solving $A\mathbf{x} = \mathbf{b}$. The smallest errors and residuals depend on the BLAS and LAPACK library that NumPy is linked to: the output below came from NumPy 2.3.5 with OpenBLAS 0.3.30 on an Apple arm64 machine, and other builds can print different values in those lines.

import numpy as np

n = 64
# One row of a synthetic image: a dark background with two bright bars.
x_true = np.full(n, 0.2)
x_true[12:20] = 0.9
x_true[36:52] = 0.6

# Gaussian blur (sigma = 1.5 pixels) of the row, written as an n x n matrix.
i = np.arange(n)
A = np.exp(-0.5 * ((i[:, None] - i[None, :]) / 1.5) ** 2)
A /= A.sum(axis=1, keepdims=True)                 # each output pixel is a weighted mean
b = A @ x_true                                    # the blurred row

print(f"cond(A) = {np.linalg.cond(A):.1e}")

def report(name, x):
    err = np.linalg.norm(x - x_true) / np.linalg.norm(x_true)
    res = np.linalg.norm(A @ x - b) / np.linalg.norm(b)
    print(f"{name:16s} relative error {err:.1e}   relative residual {res:.1e}")

report("solve(A, b)", np.linalg.solve(A, b))
V = np.linalg.inv(A)
report("inv(A) @ b", V @ b)
I = np.eye(n)
print(f"||inv(A) A - I|| = {np.linalg.norm(V @ A - I, 2):.1e}   ||A inv(A) - I|| = {np.linalg.norm(A @ V - I, 2):.1e}")

# A tiny perturbation of the data, like 8-bit quantization of the blurred row.
b_q = np.round(b * 255) / 255
print(f"relative change in b from quantization: {np.linalg.norm(b_q - b) / np.linalg.norm(b):.1e}")
x_q = np.linalg.solve(A, b_q)
print(f"relative error after quantization: {np.linalg.norm(x_q - x_true) / np.linalg.norm(x_true):.1e}")

Output:

cond(A) = 3.1e+04
solve(A, b)      relative error 8.1e-13   relative residual 2.3e-16
inv(A) @ b       relative error 8.9e-11   relative residual 4.8e-13
||inv(A) A - I|| = 1.7e-10   ||A inv(A) - I|| = 1.2e-12
relative change in b from quantization: 2.2e-03
relative error after quantization: 9.0e+00

A blur of only 1.5 pixels already gives a condition number of about $3 \times 10^4$. With exact data, solve returns a residual at machine precision and an error within $\kappa(A)$ times machine precision, about $7 \times 10^{-12}$, as the conditioning bound allows. In this run, multiplying by the computed inverse leaves a residual three orders of magnitude larger and an error about a hundred times larger. The two residuals of the inverse itself suggest why: this computed inverse is a good right inverse, $A V \approx I$, but a much poorer left inverse, and the accuracy guarantee of Druinsky and Toledo needs a good left inverse. Rounding the blurred row to 8-bit grey levels, a relative change of 0.2%, turns the solution into noise nine times larger than the image itself, an amplification of several thousand that stays below the bound $\kappa(A)$. No solver can do better on this data; deblurring needs a different problem, typically regularized least squares, which Boyd and Vandenberghe (2018, §15.3) apply to image deblurring.

Common Misconceptions

• “A matrix is either invertible or not.” Exactly, yes; numerically, a matrix with a condition number of $10^{12}$ is invertible in name only. Check the condition number or the singular values, not whether the solver raised an error.
• “A small determinant means nearly singular.” The determinant scales with the size of the entries and with the dimension: $0.1 I$ in 100 dimensions has determinant $10^{-100}$ and condition number 1. Conditioning is measured by $\kappa(A)$.
• “Compute the inverse, then multiply.” It costs more and, as the blur example shows, can lose accuracy and the small residual a direct solver guarantees. Use solve, lstsq or a factorization; form $A^{-1}$ only when its entries are themselves needed, as when composing a homography’s inverse with other transforms.
• “Rank is a fixed number.” For measured data it depends on a tolerance. Choose it from the noise level, not from machine precision, when the noise is far above machine precision.
• “Affine and linear mean the same thing.” A linear map fixes the origin; an affine map adds a translation. Homogeneous coordinates turn affine and projective maps into linear ones, at the cost of one extra coordinate and, for projective maps, a final division.
• “Rows and columns are interchangeable conventions.” They are, but formulas change with them. With one sample per row, a per-sample linear map is $X M^\top$ and the covariance is $X_c^\top X_c / n$; with one sample per column they are $M X$ and $X_c X_c^\top / n$. Mixing sources without checking their layout silently transposes results.

Where It Is Used

• Camera models and geometry. The pinhole camera model, camera parameters, homographies, the fundamental matrix and essential matrix, triangulation and camera calibration are matrices acting on homogeneous coordinates and are estimated from homogeneous linear systems, then refined by nonlinear least squares and bundle adjustment.
• Image processing. Colour conversion (color spaces), linear filtering and convolution, image warping and resampling are linear maps on pixel vectors or pixel coordinates.
• Local image analysis. The structure tensor’s eigenvalues drive the Harris corner detector, feature tracking and the aperture problem; Lucas–Kanade solves a small least-squares system per window.
• Estimation. Least squares, the Kalman filter and RANSAC’s minimal solvers are linear algebra on noisy measurements.
• Learning. Data matrices, covariance and its eigenvectors for principal component analysis, inner products in similarity measures and embeddings, and the matrix products inside every neural network layer, whose gradients backpropagation computes with transposed matrices.