LU Decomposition · aiwiki.page

LU decomposition is a matrix factorization that expresses a matrix as the product of lower and upper triangular matrices. In its simplest form, $A=LU$ , where $L$ is lower triangular and $U$ is upper triangular. Practical implementations usually incorporate row permutations, giving $PA=LU$ . The factorization records the operations of Gaussian elimination in reusable form and is a central tool in numerical linear algebra, particularly for solving systems of linear equations. (cs.cornell.edu)

Definition and existence

For a square $n\times n$ matrix, the triangular structure means

l_{ij}=0\quad(i<j),\qquad u_{ij}=0\quad(i>j).

The usual normalization sets $l_{ii}=1$ ; such an $L$ is called unit lower triangular. This convention removes the freedom to redistribute nonzero diagonal scaling between the factors. For a nonsingular matrix, a normalized unpermuted factorization exists and is unique precisely when every leading principal submatrix is nonsingular. Equivalently, all its leading principal minors are nonzero. These are the determinants of the upper-left $k\times k$ submatrices, for $k=1,\ldots,n$ . (cs.utexas.edu)

Nonsingularity alone does not guarantee an unpermuted LU factorization. For example,

A=\begin{pmatrix}0&1\\1&0\end{pmatrix}

is invertible, but its first pivot is zero. Exchanging its rows produces the identity matrix, which has trivial triangular factors. Every nonsingular square matrix admits a row-permuted factorization $PA=LU$ , with $P$ a permutation matrix. Some sources instead write $A=PLU$ , using the inverse permutation; these conventions must be distinguished when interpreting software outputs. (cs.cornell.edu)

Rectangular matrices also admit pivoted LU factorizations. For $A$ of size $m\times n$ , setting $r=\min(m,n)$ gives factors of sizes $m\times r$ and $r\times n$ , with triangular or trapezoidal structure. A singular matrix can still be factored, although zero diagonal entries in $U$ prevent ordinary triangular solution of an arbitrary system. (netlib.org)

Construction by elimination

At elimination step $k$ , assuming a nonzero pivot, the algorithm forms multipliers

l_{ik}=\frac{a^{(k)}_{ik}}{a^{(k)}_{kk}},\qquad i>k,

and updates the remaining entries:

a^{(k+1)}_{ij} =a^{(k)}_{ij}-l_{ik}a^{(k)}_{kj}, \qquad i,j>k.

Here $a^{(k)}$ denotes the matrix at the beginning of the step. The eliminated entries supply the multipliers stored in $L$ , while the remaining upper triangle becomes $U$ . The trailing update is an outer-product subtraction; in block form, the remaining matrix is a Schur complement. (cs.utexas.edu)

As a direct example,

\begin{pmatrix}2&1\\4&3\end{pmatrix} = \begin{pmatrix}1&0\\2&1\end{pmatrix} \begin{pmatrix}2&1\\0&1\end{pmatrix}.

The multiplier $2$ subtracts twice the first row from the second, leaving the second pivot equal to $1$ .

For a dense square matrix, conventional elimination requires approximately $\tfrac23n^3$ floating-point operations, counting multiplication and addition separately. Its computational complexity is therefore $O(n^3)$ . Implementations commonly store both factors in one array: $U$ occupies the diagonal and upper triangle, and the strict lower triangle stores $L$ , whose unit diagonal is implicit. (cs.utexas.edu)

Pivoting and numerical accuracy

In floating-point arithmetic, a very small pivot can produce large multipliers and magnify rounding errors. Partial pivoting chooses an entry of largest absolute value in the active column and exchanges its row with the pivot row. This bounds the elimination multipliers in magnitude by one. Complete pivoting searches the entire active submatrix and exchanges both rows and columns, producing $PAQ=LU$ . (cs.cornell.edu)

Partial pivoting generally provides good numerical stability, but it is not an unconditional guarantee. Exceptional matrices can exhibit substantial growth in intermediate entries. Backward-error analysis relates the computed factors to a slightly perturbed input matrix, with error bounds depending on the factor magnitudes. Even a small backward error need not imply a small solution error: sensitivity also depends on the matrix’s condition number. (cs.cornell.edu)

Solving systems and computing determinants

Once $PA=LU$ is available, solving $Ax=b$ reduces to

Ly=Pb,\qquad Ux=y.

Forward substitution solves the first triangular system, and backward substitution solves the second. Because $L$ has unit diagonal,

y_i=(Pb)_i-\sum_{j<i}l_{ij}y_j, \qquad x_i=\frac{y_i-\sum_{j>i}u_{ij}x_j}{u_{ii}}.

Each new right-hand side costs $O(n^2)$ operations, so the expensive factorization can be reused. (cs.cornell.edu)

For unit-diagonal $L$ ,

\det(A)=(-1)^s\prod_{i=1}^{n}u_{ii},

where $s$ is the number of row exchanges. The permutation sign is essential. Thus, after factorization, the determinant requires only a diagonal product and parity calculation. (cs.cornell.edu)

Structured matrices and related factorizations

For a sparse matrix, elimination may create nonzero entries where the original matrix contained zeros, a phenomenon called fill-in. Sparse solvers use ordering strategies to limit storage and arithmetic. They often employ $PAQ=LU$ , with column permutations primarily reducing fill and row permutations supporting stability. (cs.cornell.edu)

Dense implementations instead exploit block updates and matrix multiplication to improve memory reuse and support parallel computing. Pivot searches and row exchanges introduce communication costs that can limit parallel performance. (cs.cornell.edu)

Related factorizations exploit different properties. For symmetric or Hermitian positive-definite matrices, Cholesky decomposition uses the symmetry and permits stable factorization without pivoting. Rectangular least-squares problems commonly use QR decomposition, while singular value decomposition supports least-squares solutions when rank deficiency must be accommodated. (cs.cornell.edu)