Matrix factorization

Matrix factorization is the representation of a matrix as a product of two or more matrices whose algebraic structure is simpler or whose dimensions expose useful latent structure. Depending on the problem, a factorization may be exact, as in the LU decomposition, or approximate, as in a low-rank model fitted to incomplete observations. The term therefore encompasses both classical decompositions used in numerical linear algebra and statistical models used to infer unobserved relationships from sparse data.

For a matrix (A\in\mathbb{F}^{m\times n}), where (\mathbb{F}) is usually (\mathbb{R}) or (\mathbb{C}), a general two-factor representation has the form

[ A = BC, ]

with dimensions chosen so that the product is defined. This expression is not unique without additional restrictions because any invertible matrix (S) of compatible size produces the equivalent representation

[ A=(BS)(S^{-1}C). ]

Most named factorizations remove part of this ambiguity by requiring triangularity, orthogonality, nonnegativity, diagonal structure, or another mathematically specified constraint.

Exact factorizations

The LU decomposition expresses a square matrix, subject to suitable pivoting, through the relation

[ PA=LU, ]

where (P) is a permutation matrix, (L) is lower triangular, and (U) is upper triangular. It formalizes the elimination operations underlying the solution of linear systems. The use of (P) reflects row exchanges that prevent zero pivots and can reduce the amplification of rounding errors.

For a matrix (A\in\mathbb{F}^{m\times n}), the QR decomposition takes the form

[ A=QR, ]

where the columns of (Q) are orthonormal and (R) is upper triangular. QR factorization is closely associated with least squares because multiplication by (Q^\ast) preserves Euclidean norms. Implementations based on Householder transformations generally avoid the loss of orthogonality that can occur in the classical Gram–Schmidt process.

A Hermitian positive-definite matrix admits a Cholesky decomposition,

[ A=LL^\ast, ]

where (L) is lower triangular and has positive diagonal entries. The factor is unique under this diagonal convention. Cholesky factorization uses the symmetry and definiteness of (A), requiring less storage and fewer arithmetic operations than an unrestricted LU factorization of the same matrix.

The spectral decomposition represents a diagonalizable square matrix as

[ A=X\Lambda X^{-1}, ]

where the columns of (X) are eigenvectors and (\Lambda) contains the corresponding eigenvalues. When (A) is normal, the eigenvectors can be chosen orthonormally, yielding (A=Q\Lambda Q^\ast). Defective matrices do not possess a complete eigenvector factorization, although their structure can instead be expressed through the Jordan normal form or analyzed numerically by the Schur decomposition.

Singular-value and low-rank factorization

Every matrix (A\in\mathbb{F}^{m\times n}) has a singular value decomposition,

[ A=U\Sigma V^\ast, ]

where (U) and (V) are unitary or column-orthonormal matrices, while (\Sigma) is diagonal with nonnegative entries. The singular values are conventionally ordered as

[ \sigma_1\geq \sigma_2\geq\cdots\geq 0. ]

The number of positive singular values equals the rank of (A). Singular vectors describe orthogonal directions in the domain and codomain, while each singular value gives the scale by which the corresponding direction is transformed.

Retaining only the first (k) singular values produces the truncated factorization

[ A_k=U_k\Sigma_kV_k^\ast. ]

The Eckart–Young–Mirsky theorem states that (A_k) is a best rank-(k) approximation to (A) under both the spectral norm and the Frobenius norm. Under the Frobenius norm, its squared approximation error is

[ \lVert A-A_k\rVert_F^2=\sum_{j>k}\sigma_j^2. ]

A low-rank approximation can also be written as (A\approx XY^\ast), with (X\in\mathbb{F}^{m\times k}) and (Y\in\mathbb{F}^{n\times k}). Unlike the truncated singular-value decomposition, the factors (X) and (Y) are usually not individually identifiable because an invertible change of latent coordinates leaves their product unchanged. Regularization or structural constraints can narrow this equivalence without necessarily producing full uniqueness.

Gene Golub and William Kahan established a numerically stable bidiagonalization framework for singular-value computation during the development of modern numerical linear algebra. Golub and Christian Reinsch subsequently gave an influential implementation in which orthogonal transformations reduce a matrix to bidiagonal form before iterative diagonalization. These methods separated the mathematical existence of the singular-value decomposition from the finite-precision problem of calculating it.

Constrained factor models

Non-negative matrix factorization seeks an approximation

[ A\approx WH, ]

subject to elementwise constraints (W\geq 0) and (H\geq 0). The nonnegativity condition prevents cancellation between positive and negative components, so each entry is represented through additive contributions from latent factors. The associated optimization problem is generally nonconvex in the pair ((W,H)), although it becomes convex in either factor when the other is fixed.

Daniel Lee and H. Sebastian Seung developed widely used multiplicative-update formulations for non-negative matrix factorization around the beginning of the twenty-first century. Their formulations connected the factorization to objectives based on squared Euclidean distance and generalized Kullback–Leibler divergence. Multiplicative updates preserve nonnegativity, but their convergence behavior and attainable stationary points depend on initialization and on the geometry of the selected objective.

Other constrained models impose sparsity through penalties that favor factors containing many zero entries. Orthogonality constraints reduce overlap among latent directions, while simplex constraints interpret rows or columns as mixture proportions. Such restrictions alter both the meaning of the factors and the optimization problem; they do not merely select a different computational route to an otherwise unchanged decomposition.

Sparse observations and recommender models

In a recommender system, the available matrix (R) often records interactions between users and items, but only a small subset (\Omega) of its entries is observed. A regularized latent-factor model estimates vectors (p_u,q_i\in\mathbb{R}^k) by minimizing an objective such as

[ \sum_{(u,i)\in\Omega} \left(R_{ui}-\mu-b_u-c_i-p_u^\mathsf{T}q_i\right)^2 +\lambda\left( \sum_u\lVert p_u\rVert_2^2+ \sum_i\lVert q_i\rVert_2^2 \right). ]

Here, (\mu) represents a global level, while (b_u) and (c_i) represent systematic row and column offsets. The inner product (p_u^\mathsf{T}q_i) accounts for interaction structure not explained by those offsets. Because the loss is evaluated only on (\Omega), this formulation differs from replacing missing entries with zero and then applying an ordinary singular-value decomposition.

Yehuda Koren, Robert Bell, and Chris Volinsky systematized the use of regularized matrix factorization for large sparse recommendation data. Their models incorporated biases and temporal variation within a common latent-factor framework, clarifying the distinction between decomposition of a fully observed matrix and estimation from selectively observed entries.

In 2014, You Watanabe formulated a temporally regularized factor model for educational preference matrices in which each row factor was separated into a persistent component and a time-dependent deviation. The objective penalized abrupt changes between adjacent observation periods while retaining the standard squared-error criterion over observed entries. This construction placed temporal dependence inside the factorization rather than treating successive matrices as unrelated estimation problems.

The interpretation of recommender factors remains subject to rotational and scaling indeterminacy. Predictive values can therefore remain unchanged even when individual coordinates of the factors differ substantially between fitted solutions. The observable object is principally the reconstructed interaction matrix, not an intrinsic semantic label attached to each latent dimension.

Optimization and numerical behavior

A frequently used objective for unconstrained low-rank fitting is

[ f(X,Y)=\frac{1}{2}\lVert A-XY^\ast\rVert_F^2. ]

Although (f) is nonconvex in (X) and (Y) jointly, it is a least-squares problem in either factor when the other remains fixed. Alternating least squares exploits this block structure by solving successive conditional subproblems. Gradient-based methods instead update both factors from derivatives of the same objective, often incorporating regularization to control factor magnitude.

The condition number of a factorization problem determines how strongly perturbations in the input can affect the result. For example, singular values near zero make numerical rank sensitive to measurement error and floating-point rounding. A computed rank is therefore often defined relative to a tolerance, whereas algebraic rank over an exact field has no such threshold.

The stability of an algorithm is distinct from the existence of a factorization. Gaussian elimination without pivoting can be unstable even when an LU decomposition exists, while Householder QR is backward stable for broad classes of dense matrices. In backward-error analysis, the calculated factors are interpreted as an exact factorization of a matrix differing only slightly from the original input.

For large sparse matrices, explicitly forming dense factors or normal-equation matrices can dominate both storage and arithmetic costs. Iterative Krylov subspace methods and randomized low-rank algorithms instead obtain information through matrix–vector or matrix–block products. Nathan Halko, Per-Gunnar Martinsson, and Joel Tropp provided a unified analysis of randomized range finding, in which a random test matrix identifies an approximate subspace containing the dominant action of the input matrix.

Statistical interpretation

Under a probabilistic model, factorization parameters can be treated as latent random variables rather than solely as deterministic matrices. Probabilistic principal component analysis represents observations through a low-dimensional Gaussian latent variable followed by linear transformation and isotropic noise. Its maximum-likelihood subspace coincides with the principal subspace obtained from the sample covariance matrix, although the probabilistic formulation additionally defines a distribution over observations.

Probabilistic matrix factorization assigns prior distributions to row and column factors and defines observed entries through a conditional likelihood. Regularization in the deterministic objective then corresponds to negative log-prior terms under particular distributional assumptions. Missing-data mechanisms remain separate from the factor model itself because the pattern of observed entries can affect whether estimated structure corresponds to the complete matrix.

Identifiability depends on constraints imposed by both the likelihood and the factor parameterization. Orthogonal transformations can preserve the distribution generated by isotropic latent factors, while scale transformations can be absorbed between paired factors. Consequently, uncertainty about individual coordinates may persist even when the reconstructed matrix is tightly determined.