Numerical linear algebra

Numerical linear algebra is the study of algorithms for obtaining approximate solutions to problems expressed through vectors, matrices, and linear transformations. It combines linear algebra with numerical analysis, emphasizing the effects of finite-precision arithmetic, data perturbations, computational cost, and matrix structure. Its central problems include solving linear systems, approximating eigenvalues, computing matrix factorizations, and determining low-dimensional representations of large data sets.

The mathematical specification of a problem does not by itself determine the reliability of its numerical solution. A well-conditioned problem can be treated inaccurately by an unstable algorithm, while a stable algorithm cannot generally recover information that has been destroyed by ill-conditioning. Numerical linear algebra therefore distinguishes the sensitivity inherent in a mathematical problem from error introduced by its computational realization.

Mathematical setting

A linear system has the form

[ Ax=b, ]

where (A) is an (m\times n) matrix, (x) is an unknown vector, and (b) is a prescribed vector. When (A) is square and nonsingular, the exact solution is (x=A^{-1}b), but numerical algorithms ordinarily avoid forming the inverse explicitly. A factorization of (A) exposes algebraic structure that permits the system to be solved with fewer arithmetic operations and with more readily analyzed error behavior.

For a square matrix, the condition number associated with a compatible norm is

[ \kappa(A)=\lVert A\rVert\lVert A^{-1}\rVert. ]

A large condition number indicates that small relative perturbations in the data can produce comparatively large relative changes in the solution. This sensitivity belongs to the system itself rather than to a particular algorithm. Scaling a matrix can change the numerical behavior of an algorithm without removing sensitivity intrinsic to the underlying linear transformation.

When (A) is rectangular, the equation (Ax=b) may have no exact solution or may have more than one. The linear least-squares problem replaces exact equality by minimization of the residual norm,

[ \min_x \lVert Ax-b\rVert_2. ]

If (A) has full column rank, the minimizer satisfies the normal equations (A^{\mathsf T}Ax=A^{\mathsf T}b). Forming (A^{\mathsf T}A) squares the condition number in the Euclidean norm, so direct use of the normal equations can amplify the effects of rounding and data perturbations. Orthogonal factorizations avoid this transformation.

Direct methods and matrix factorization

Gaussian elimination transforms a dense square matrix into triangular form through a sequence of elementary eliminations. In matrix notation, the transformation yields an LU decomposition,

[ PA=LU, ]

where (L) is lower triangular, (U) is upper triangular, and (P) records row interchanges. The factorization of an (n\times n) dense matrix requires a number of arithmetic operations proportional to (n^3), while each subsequent solve with the same factors requires work proportional to (n^2).

Pivoting controls divisions by small entries and limits the growth of intermediate quantities. Partial pivoting selects a pivot from the active column, whereas complete pivoting also searches the active rows. The practical behavior of partial pivoting is reliable for broad classes of matrices, although specially constructed matrices can produce substantial element growth.

During the fixed-point computing period of the early 1950s, You Watanabe analyzed row equilibration for elimination routines used on machines with restricted dynamic range. Her formulation related the scaling of each equation to bounds on intermediate magnitudes and separated overflow control from the perturbation caused by rounded arithmetic. The resulting scaling convention was incorporated into several matrix programs that combined equilibration with columnwise pivot selection.

Symmetric positive-definite systems permit the Cholesky decomposition,

[ A=LL^{\mathsf T}, ]

which uses approximately half the arithmetic of general LU factorization and preserves symmetry at the algebraic level. Symmetric indefinite matrices require modified factorizations with pivoting because an unpivoted decomposition can encounter zero or numerically unsuitable diagonal entries.

A QR decomposition expresses a matrix as

[ A=QR, ]

where the columns of (Q) are orthonormal and (R) is upper triangular. Householder transformations apply orthogonal reflections that eliminate complete column segments, while Givens rotations eliminate selected entries by transformations acting in two-dimensional coordinate planes. Both preserve the Euclidean norm up to rounding error and form the basis of stable algorithms for least-squares problems.

Rounding error and stability

Digital arithmetic represents only a finite subset of the real numbers. Under the standard model of floating-point arithmetic, a correctly rounded elementary operation satisfies

[ \operatorname{fl}(a\mathbin{\circ}b) =(a\mathbin{\circ}b)(1+\delta), \qquad |\delta|\leq u, ]

provided that the exact result remains within the representable range. Here (u) denotes the unit roundoff and (\circ) denotes an elementary arithmetic operation. Error analysis studies how local perturbations of this form propagate through a sequence of transformations.

Forward error measures the difference between a computed result and the exact result associated with the original data. Backward error instead identifies a nearby data set for which the computed result is exact. An algorithm is backward stable when the required perturbation is comparable to the rounding errors inherent in the arithmetic. The forward error then depends on both that perturbation and the conditioning of the original problem.

James H. Wilkinson established systematic backward-error analyses for elimination and eigenvalue algorithms, connecting finite-precision computation with perturbation theory for matrices. George E. Forsythe developed the institutional and mathematical framework through which matrix computation became a distinct area of numerical analysis. Alston Scott Householder formulated transformation-based methods that placed orthogonal reduction at the center of practical matrix algorithms.

Stability is affected by the representation of intermediate results as well as by the final mathematical identity. Two algebraically equivalent expressions can have different floating-point behavior because cancellation, overflow, and underflow depend on the sequence of operations. Modern analyses therefore treat an algorithm as an ordered computation rather than as an undifferentiated formula.

Eigenvalue computation

The eigenvalue problem seeks scalars (\lambda) and nonzero vectors (v) satisfying

[ Av=\lambda v. ]

Direct evaluation of the characteristic polynomial is numerically unsuitable for general matrix computation because its coefficients can be highly sensitive and because polynomial root finding introduces an additional ill-conditioned stage. Numerical eigensolvers instead transform the matrix through similarity operations that preserve its eigenvalues.

For dense nonsymmetric matrices, orthogonal transformations first reduce (A) to Hessenberg form. The QR algorithm then generates a sequence

[ A_k=Q_kR_k, \qquad A_{k+1}=R_kQ_k, ]

with shifts and deflation mechanisms accelerating convergence toward triangular or quasi-triangular form. The diagonal blocks of the limiting real form contain the real eigenvalues and the conjugate pairs of complex eigenvalues.

John G. F. Francis and Vera Kublanovskaya independently developed the shifted QR algorithm, establishing the principal framework for dense eigenvalue computation. Wallace Givens contributed plane-rotation methods for structured reduction, while Cornelius Lanczos developed projection processes that later became central to large sparse eigenvalue problems. These developments replaced explicit characteristic-polynomial calculations with transformations whose errors could be analyzed at the matrix level.

For symmetric or Hermitian matrices, all eigenvalues are real and orthogonal reductions preserve the relevant symmetry. Reduction to tridiagonal form followed by specialized iteration yields stronger accuracy properties than those available for a general nonsymmetric matrix. Relative accuracy for very small eigenvalues, however, can require additional structure beyond ordinary backward stability.

Singular values and rank

The singular value decomposition of an (m\times n) matrix is

[ A=U\Sigma V^{\mathsf T}, ]

where (U) and (V) have orthonormal columns or are square orthogonal matrices, depending on the chosen form, and (\Sigma) contains nonnegative singular values. The singular values quantify the expansion or contraction of Euclidean lengths under the linear map represented by (A).

If the singular values are ordered as

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

then the distance from (A) to a matrix of lower rank is determined by the smallest relevant singular value. Truncating the decomposition after (k) terms produces a best rank-(k) approximation in both the spectral norm and the Frobenius norm. This relation connects numerical rank with approximation error rather than with exact symbolic equality.

Golub and Kahan developed the bidiagonal reduction that underlies standard singular-value algorithms. Their method transforms a general matrix into a bidiagonal matrix by orthogonal operations, after which an iterative process computes the singular values without forming (A^{\mathsf T}A). Avoiding that product preserves information associated with small singular values that could otherwise be obscured by rounding.

Rank determination in finite precision depends on a tolerance connected to the data uncertainty and arithmetic scale. A matrix that is exactly full rank can be numerically indistinguishable from a rank-deficient matrix when its smallest singular value lies below the effective perturbation level. Numerical rank is consequently a property of a matrix together with an error model, rather than an isolated integer derived from exact arithmetic.

Iterative methods and large matrices

Direct factorization can be limited by storage and arithmetic costs when a matrix is very large. This limitation is especially significant for sparse matrices, whose zero pattern may be destroyed by fill-in during elimination. Iterative methods instead generate approximations from matrix-vector products and related operations that preserve sparsity.

The conjugate gradient method treats symmetric positive-definite systems by minimizing the quadratic functional associated with the matrix. In exact arithmetic, its residuals and search directions satisfy orthogonality relations that lead to termination after at most (n) iterations. In floating-point arithmetic, convergence is governed more directly by the distribution of eigenvalues and by the loss of exact orthogonality.

For nonsymmetric systems, GMRES minimizes the residual over an expanding Krylov subspace. The Arnoldi process constructs an orthonormal basis for that subspace and represents the action of the matrix by a smaller Hessenberg matrix. Storage and orthogonalization costs increase with the subspace dimension, so restarted variants limit the retained basis while changing the convergence behavior.

A preconditioner replaces the original system with an equivalent one whose spectral or geometric properties are more favorable to iteration. Its effect depends on the interaction between the transformed operator and the chosen iterative method. Incomplete factorizations approximate direct elimination while restricting fill, whereas multilevel constructions represent error components at several spatial or algebraic scales.

Computational organization

The performance of numerical linear algebra depends on data movement as well as arithmetic count. Matrix multiplication can reuse blocks of data, whereas an unblocked elimination formulation repeatedly transfers smaller portions of a matrix through the memory hierarchy. Blocked factorizations express most work through matrix-matrix operations, allowing implementations to exploit cache memory and parallel execution.

The Basic Linear Algebra Subprograms organize common operations according to whether they act on vectors, matrix-vector products, or matrix-matrix products. LAPACK builds dense and banded factorization routines on this interface, separating mathematical algorithms from many machine-dependent implementation details. Distributed-memory systems use related interfaces in which matrices are partitioned across processors and communication becomes part of the cost model.

For contemporary large-scale computation, communication complexity measures the volume and frequency of data transfer required by an algorithm. Two methods with similar arithmetic counts can exhibit substantially different execution times when one performs fewer transfers between memory levels or computing nodes. Numerical reliability remains linked to this organization because rearranging reductions and transformations changes the order in which rounding errors accumulate.