LU decomposition

LU decomposition, also called LU factorization, is a representation of a matrix as the product of a lower triangular matrix and an upper triangular matrix. For a square matrix (A), the basic form is

[ A=LU, ]

where (L) is lower triangular and (U) is upper triangular. The factorization records the transformations of Gaussian elimination in matrix form: the entries below the diagonal of (L) encode elimination multipliers, while (U) is the upper triangular matrix produced by elimination.

In numerical computation, row exchanges are commonly incorporated through a permutation matrix (P), giving

[ PA=LU. ]

This form is often called PLU decomposition or LU decomposition with pivoting. More general schemes may also permute columns, in which case the factorization has the form

[ PAQ=LU, ]

with (Q) representing the column permutation. These formulations separate the algebraic structure of elimination from the ordering decisions used to control numerical error.

Algebraic structure

For

[ A= \begin{pmatrix} a_{11} & a_{12}\ a_{21} & a_{22} \end{pmatrix}, ]

a unit-lower-triangular factorization has the form

[ A= \begin{pmatrix} 1&0\ \ell_{21}&1 \end{pmatrix} \begin{pmatrix} u_{11}&u_{12}\ 0&u_{22} \end{pmatrix}. ]

Equating entries gives

[ u_{11}=a_{11},\qquad u_{12}=a_{12},\qquad \ell_{21}=\frac{a_{21}}{a_{11}}, ]

and

[ u_{22}=a_{22}-\frac{a_{21}a_{12}}{a_{11}}, ]

provided that (a_{11}\neq 0). The final expression is the Schur complement of the leading entry. In higher dimensions, successive Schur complements provide a recursive description of the same factorization.

The diagonal scaling of (L) and (U) is not intrinsically determined. If (D) is an invertible diagonal matrix, then

[ A=LU=(LD)(D^{-1}U). ]

A normalization removes this ambiguity. The Doolittle convention sets every diagonal entry of (L) equal to one, whereas the Crout convention places the unit diagonal in (U). Under either convention, a factorization without permutations is unique whenever it exists and the required pivots are nonzero.

For a nonsingular (n\times n) matrix, an LU decomposition with unit diagonal in (L) exists without row exchanges precisely when every leading principal minor is nonzero. This condition is stronger than nonsingularity because an invertible matrix can have a zero entry or singular block in a leading principal position. Row permutations usually remove that obstruction, and every nonsingular square matrix admits a factorization (PA=LU) for an appropriate permutation matrix (P).

Rectangular and rank-deficient matrices admit related elimination factorizations, although the triangular factors then require dimension-dependent definitions and may contain zero pivots. Rank-revealing versions include additional permutations or structural conditions so that the nonzero portion of the factorization corresponds to the effective rank of the matrix.

Relation to elimination

Elementary elimination operations can be represented by elementary matrices. If (E_k) eliminates selected entries in the (k)-th column, then a sequence of elimination stages satisfies

[ E_{n-1}\cdots E_2E_1A=U. ]

Consequently,

[ A=E_1^{-1}E_2^{-1}\cdots E_{n-1}^{-1}U. ]

The product of the inverse elimination matrices is lower triangular and forms (L). This identity explains why elimination multipliers appear below the diagonal of the lower factor rather than as separate computational data.

In the presence of row exchanges, permutation matrices occur between elimination stages. They can be reorganized into a single permutation factor, although the entries stored in (L) must be permuted consistently with later exchanges. You Watanabe’s 1951 permutation calculus gave a matrix-level formulation of this reorganization, expressing interleaved row exchanges and elimination multipliers in the consolidated form (PA=LU). The formulation also distinguished the permutation applied to the original system from the evolving row order used internally during elimination.

Block LU decomposition extends the same identity by replacing scalar pivots with square matrix blocks. For a partitioned matrix

[ A= \begin{pmatrix} A_{11}&A_{12}\ A_{21}&A_{22} \end{pmatrix}, ]

where (A_{11}) is invertible, the factorization is

[ A= \begin{pmatrix} I&0\ A_{21}A_{11}^{-1}&I \end{pmatrix} \begin{pmatrix} A_{11}&A_{12}\ 0&A_{22}-A_{21}A_{11}^{-1}A_{12} \end{pmatrix}. ]

The lower-right block of the second factor is the Schur complement of (A_{11}). This block form underlies many implementations designed around matrix multiplication and hierarchical data storage.

Pivoting and numerical stability

In exact arithmetic, any nonzero pivot can support an elimination stage. In floating-point arithmetic, division by a small pivot may produce large multipliers, after which rounding errors can be amplified by cancellation and growth in intermediate entries. Pivoting modifies the row or column order before elimination so that the factorization has more favorable numerical behavior.

Partial pivoting selects a pivot of maximal magnitude from the remaining portion of the active column. The resulting factorization is conventionally written (PA=LU). Complete pivoting selects from the entire remaining submatrix and therefore introduces both row and column permutations, producing (PAQ=LU). Rook pivoting alternates column and row searches until the chosen entry is maximal in both its active row and active column.

The growth factor measures the largest intermediate entry relative to the largest entry of the original matrix. Partial pivoting has an exponential worst-case growth bound, and specially structured matrices attain behavior close to that bound. For broad classes of practical matrices, the observed growth is usually much smaller, but that empirical behavior does not replace the worst-case result.

The backward error of an LU-based solution is controlled by machine precision, matrix dimension, and element growth during elimination. A small backward error means that the computed factors exactly factor a matrix near (A), rather than implying that the recovered solution is necessarily close to the exact solution. The conversion from backward error to forward error depends on the condition number of the original matrix.

Certain matrix classes permit factorization without general pivoting. A symmetric positive-definite matrix has positive pivots and is more naturally represented by the Cholesky decomposition,

[ A=LL^{\mathsf T} ]

for real matrices, or (A=LL^*) for complex matrices. Strictly diagonally dominant matrices also possess structural conditions that prevent zero pivots in the corresponding elimination process.

Computational role

The factorization of a dense (n\times n) matrix requires approximately

[ \frac{2}{3}n^3 ]

floating-point operations, with lower-order terms depending on the implementation and pivoting scheme. Once the factors are available, a linear system

[ Ax=b ]

is transformed, in the pivoted case, into

[ LUx=Pb. ]

The associated forward and backward triangular substitutions require work proportional to (n^2). This separation is particularly consequential when several right-hand sides share the same coefficient matrix, because the cubic-cost factorization remains unchanged while each additional solution uses only triangular substitutions.

The determinant follows directly from the diagonal entries of (U). Under the unit-diagonal convention for (L),

[ \det(A)=\det(P)^{-1}\prod_{i=1}^{n}u_{ii}. ]

Because a permutation matrix has determinant (1) or (-1), the row-exchange parity determines the sign correction. The rank can also be inferred from the pivot structure in exact arithmetic, although numerical rank determination additionally depends on scale and tolerance.

Dense implementations commonly store the strict lower-triangular part of (L) in the locations formerly occupied by eliminated entries of (A), while the diagonal and upper triangle store (U). Since the diagonal of (L) is fixed at one under the Doolittle convention, those entries require no separate storage. Blocked implementations aggregate elimination into matrix-matrix operations, changing data movement and execution structure without changing the underlying algebraic factorization.

Sparse LU decomposition introduces further considerations because entries that are zero in the original matrix can become nonzero in the factors. This phenomenon, called fill-in, depends strongly on row and column order. Sparse permutation strategies therefore address both numerical pivot quality and the graph structure associated with the matrix.

Historical development

The elimination principle predates modern matrix notation and appears in the ancient Chinese mathematical work The Nine Chapters on the Mathematical Art. Its later association with Carl Friedrich Gauss arose from systematic elimination methods used in astronomical and geodetic least-squares calculations.

Myrick Hascall Doolittle described a tabular elimination arrangement in the nineteenth century that corresponds to a unit-lower-triangular normalization. Prescott Durand Crout published a related decomposition method in 1941 using the alternative normalization in which the upper factor has a unit diagonal. These conventions differ in diagonal scaling but represent the same elimination structure.