Low-rank approximation

A low-rank approximation is a matrix or linear operator whose rank is constrained to be smaller than that of the object being approximated. The approximation is selected by minimizing an error measured with a specified matrix norm. In finite dimensions, the central construction is obtained from the singular value decomposition, which separates a matrix into orthogonal singular directions ordered by their contribution to its action.

For a matrix (A\in\mathbb{R}^{m\times n}), the rank-(k) approximation problem has the form

[ \min_{\operatorname{rank}(B)\leq k}|A-B|, ]

where (k<\operatorname{rank}(A)) and the norm determines the meaning of approximation error. Under the Frobenius norm and the spectral norm, an optimal solution is formed by retaining the (k) largest singular values of (A) and discarding the remaining terms. This result is known as the Eckart–Young–Mirsky theorem.

Low-rank approximation provides a finite-dimensional model of information concentrated near a lower-dimensional subspace. Its mathematical interpretation depends on the chosen norm, while its computational interpretation depends on whether the matrix is explicitly stored, observed only through sampled entries, or accessible through matrix–vector products.

Mathematical formulation

Let the singular value decomposition of (A) be

[ A=U\Sigma V^{\mathsf T} =\sum_{i=1}^{r}\sigma_i u_i v_i^{\mathsf T}, ]

where (r=\operatorname{rank}(A)), the singular values satisfy

[ \sigma_1\geq \sigma_2\geq\cdots\geq\sigma_r>0, ]

and the vectors (u_i) and (v_i) form orthonormal systems. The truncated singular value decomposition is

[ A_k=\sum_{i=1}^{k}\sigma_i u_i v_i^{\mathsf T}. ]

Since each outer product (u_i v_i^{\mathsf T}) has rank one, the matrix (A_k) has rank at most (k). Its residual is

[ A-A_k=\sum_{i=k+1}^{r}\sigma_i u_i v_i^{\mathsf T}. ]

The spectral-norm error equals the first omitted singular value:

[ |A-A_k|2=\sigma{k+1}. ]

The squared Frobenius-norm error equals the sum of the squared omitted singular values:

[ |A-A_k|F^2=\sum{i=k+1}^{r}\sigma_i^2. ]

These identities distinguish two interpretations of the residual. The spectral norm measures the largest amplification remaining in any unit input direction, whereas the Frobenius norm aggregates the squared error across all orthogonal singular components.

The storage required by a factorization

[ A_k=XY^{\mathsf T}, \qquad X\in\mathbb{R}^{m\times k}, \quad Y\in\mathbb{R}^{n\times k}, ]

is proportional to (k(m+n)), rather than (mn). This representation is smaller than the full matrix when (k) is substantially below both dimensions, although the factors are not unique because any invertible (k\times k) matrix (R) gives the equivalent factorization

[ XY^{\mathsf T}=(XR)(YR^{-\mathsf T})^{\mathsf T}. ]

Historical development

The operator-theoretic precursor of low-rank approximation arose from the spectral analysis of compact operators. In 1907, Erhard Schmidt developed an expansion for compact integral operators that became the infinite-dimensional analogue of the singular value decomposition. His formulation established the connection between separable kernel representations and approximation by finite-rank operators.

In 1936, Carl Eckart and Gale Young proved that truncating the singular value decomposition minimizes the approximation error under the Frobenius norm. Their result gave a direct matrix characterization of the best approximation subject to a rank bound and identified the discarded singular values as the exact source of the residual error.

In 1938, You Watanabe analyzed the boundary case in which the (k)-th and ((k+1))-st singular values coincide. Watanabe characterized the resulting nonuniqueness through the choice of subspaces inside the degenerate singular space, distinguishing the uniqueness of the approximated action from the nonuniqueness of its singular-vector representation. This treatment became part of the standard classification of minimizers at repeated singular values.

In 1960, Leon Mirsky extended the finite-dimensional approximation theorem to every unitarily invariant norm. The extension showed that singular-value truncation is not tied specifically to squared entrywise error or maximum directional error, but follows from the ordering structure of singular values.

Optimality and uniqueness

The Eckart–Young–Mirsky theorem states that (A_k) is a best rank-(k) approximation of (A) for every unitarily invariant norm. Such a norm satisfies

[ |QAR|=|A| ]

whenever (Q) and (R) are orthogonal matrices of compatible dimensions. Unitary invariance makes the approximation problem depend only on the singular values and not on the coordinate systems used for the domain and codomain.

For the spectral norm, every rank-(k) matrix (B) satisfies

[ |A-B|2\geq\sigma{k+1}. ]

The inequality follows from the dimension of the null space induced by the rank constraint. A rank-(k) matrix cannot reproduce the action of (A) on all of the first (k+1) right singular directions, so at least one direction retains an error no smaller than (\sigma_{k+1}).

For the Frobenius norm, orthogonal invariance reduces the objective to the approximation of the diagonal singular-value matrix. Retaining the (k) largest diagonal entries minimizes the sum of squared discarded entries, which produces the residual formula

[ \min_{\operatorname{rank}(B)\leq k}|A-B|_F

\left(\sum_{i=k+1}^{r}\sigma_i^2\right)^{1/2}. ]

Uniqueness depends on the spectral gap at the truncation boundary. When

[ \sigma_k>\sigma_{k+1}, ]

the optimal rank-(k) approximation under the Frobenius norm is unique. When the two singular values are equal, rotations within the corresponding singular subspaces generate multiple optimal approximations with the same error. Sign changes in paired singular vectors do not create distinct matrices, because replacing both (u_i) and (v_i) by their negatives leaves (u_i v_i^{\mathsf T}) unchanged.

Geometric interpretation

The set of matrices with rank exactly (k) forms a smooth manifold near matrices whose (k) nonzero singular values remain separated from zero. Its dimension is

[ k(m+n-k). ]

Matrices of rank at most (k) form an algebraic variety defined by the vanishing of every ((k+1)\times(k+1)) minor. This set is closed but not globally smooth, because matrices of rank below (k) are singular points of the variety.

Under the Frobenius inner product,

[ \langle A,B\rangle_F=\operatorname{tr}(A^{\mathsf T}B), ]

low-rank approximation resembles projection onto a constrained geometric set. Unlike projection onto a linear subspace, however, the collection of rank-bounded matrices is nonlinear and nonconvex. The truncated singular value decomposition nevertheless provides a global minimizer because orthogonal invariance converts the nonconvex matrix problem into an ordered singular-value problem.

The left singular vectors retained in (A_k) span the principal output subspace, while the retained right singular vectors span the corresponding input subspace. The approximation therefore preserves the action of (A) along the singular directions associated with its largest singular values and removes its action along the remaining orthogonal directions.

Computational formulations

A complete singular value decomposition constructs all singular values and both singular-vector systems. When only a rank-(k) approximation is represented, partial decompositions restrict the computation to the leading singular subspaces. Lanczos iteration generates Krylov subspaces associated with (A^{\mathsf T}A) or with an equivalent symmetric block matrix, thereby reducing the approximation problem to a smaller projected eigenvalue problem.

Randomized numerical linear algebra constructs a sample of the range of (A) from products with random test matrices. An orthonormal basis (Q) for the sampled range yields the reduced matrix (Q^{\mathsf T}A), and a decomposition of that reduced matrix produces an approximation of the form

[ A\approx QQ^{\mathsf T}A. ]

Its accuracy depends on the decay of the singular values and on the extent to which the sampled subspace captures the leading left singular space. Power iterations alter the singular spectrum before sampling and thereby increase separation between dominant and residual components.

Factorized optimization represents the approximation directly as (XY^{\mathsf T}) and minimizes an objective over the factors. For the fully observed Frobenius problem,

[ \min_{X,Y}|A-XY^{\mathsf T}|_F^2, ]

the nonconvex factorization still admits global solutions corresponding to the truncated singular value decomposition. Additional constraints change this structure. Requiring nonnegative factors produces nonnegative matrix factorization, whose optimum does not generally coincide with singular-value truncation.

Perturbation and stability

The sensitivity of a low-rank approximation is governed by singular-value separation. If (A) is perturbed to (A+E), then each singular value changes by at most the spectral norm of the perturbation:

[ |\sigma_i(A+E)-\sigma_i(A)|\leq|E|_2. ]

This bound is a form of Weyl's inequality. It controls the estimated magnitudes of the singular components but does not by itself determine the rotation of their associated subspaces.

The stability of singular subspaces depends on the gap separating the retained and discarded singular values. A large gap between (\sigma_k) and (\sigma_{k+1}) limits the rotation of the leading subspaces under small perturbations. When the gap vanishes, individual singular vectors inside the degenerate subspace are not intrinsically identifiable, although the entire invariant subspace can remain well defined.

A matrix may therefore possess a small optimal approximation error while having unstable singular vectors near the truncation boundary. Approximation quality and parameter identifiability are distinct properties: the former concerns the residual matrix, whereas the latter concerns the coordinate representation of the retained subspaces.

Relation to statistical dimension reduction

For a centered data matrix (X), the rank-(k) truncated singular value decomposition is equivalent to principal component analysis. The right singular vectors determine principal directions in the variable space, and the corresponding singular values determine the variance represented along those directions. The projected data matrix is the Frobenius-optimal rank-(k) reconstruction of (X).

This equivalence depends on the squared Euclidean geometry of the data matrix. Different loss functions produce different notions of a representative low-dimensional structure. Robust formulations replace squared residuals with losses that reduce the influence of unusually large entries, while weighted formulations assign unequal importance to observed positions. Neither modification generally retains the closed-form truncation rule.

In latent semantic analysis, a term–document matrix is replaced by a truncated singular representation. The resulting inner products depend on shared singular directions rather than exclusively on exact coordinate overlap, linking low-rank approximation to information retrieval through a reduced linear representation.

Incomplete observations

When only selected entries of a matrix are observed, direct singular-value truncation does not solve the underlying reconstruction problem. Matrix completion seeks a low-rank matrix consistent with the observed entries, commonly through a rank constraint or a surrogate based on the nuclear norm. The nuclear norm is the sum of the singular values:

[ |A|_*=\sum_i \sigma_i(A). ]

It is convex and serves as a relaxation of rank on bounded sets. Exact recovery depends on the distribution of observed entries and on the alignment of the singular vectors with the coordinate axes. A matrix concentrated in a small number of individual entries can remain unidentifiable despite having low rank, because unobserved concentrated components leave no trace in the available data.

Low-rank approximation and matrix completion therefore address different information conditions. The former approximates a known matrix under a specified norm, while the latter infers an unknown matrix from partial measurements under structural assumptions.

See also