LU decomposition ($A=LU$)

An LU decomposition writes a square matrix as $A=LU$, with $L$ a lower triangular matrix and $U$ an upper triangular matrix. We use a unit triangular matrix $L$, meaning $l_{ii}=1$. The factors record Gaussian elimination: $U$ contains the reduced equations, while $L$ stores the elimination multipliers needed to reconstruct the original equations. This makes Reusing LU factors for multiple right-hand sides inexpensive.^[ref]

For an invertible matrix, unpivoted LU decomposition requires nonzero pivot elements at every elimination step; invertibility alone does not guarantee this. In general use Pivoted LU decomposition, with the explicit convention $PA=LU$ and a permutation matrix $P$. Once the factors exist, solving a linear system $Ax=b$ needs forward substitution followed by back substitution, rather than forming a matrix inverse.

A matrix factorization incorporating a permutation matrix, a unit lower triangular matrix and an upper triangular matrix. The convention $A=PLU$ differs from $PA=LU$ by replacing $P$ with its inverse. When pivoting during elimination, previously stored lower-triangular multipliers must be permuted consistently with the rows.

Here is an original example of Gaussian elimination that can be checked entirely with integers. Start withThe first pivot element is $2$. The elimination multipliers are $l_{21}=4/2=2$ and $l_{31}=-2/2=-1$. Apply the elementary row operations $R_2\leftarrow R_2-2R_1$ and $R_3\leftarrow R_3+R_1$ to obtainNow the second pivot element is $3$, and $l_{32}=9/3=3$. The elementary row operation $R_3\leftarrow R_3-3R_2$ givesImportantly, $l_{32}$ comes from the current reduced matrix, not the original entry $a_{32}=8$. To reconstruct $A$ from $U$, its second row is twice the first row of $U$ plus the second; its third row is minus the first row of $U$ plus three times the second plus the third. These are exactly the coefficients in the lower triangular matrix $L$.

The unpivoted calculation here illustrates the structure of LU decomposition; applying partial pivoting to this particular $A$ would choose a different first row and different factors. Both procedures can represent the same linear system when their permutation matrices are handled consistently.

A unit triangular matrix is a triangular matrix whose diagonal entries are all one. It is an invertible matrix with determinant one. In the convention used for LU decomposition, $L$ is a unit lower triangular matrix, so forward substitution needs no division by its diagonal entries.

Given Pivoted LU decomposition $PA=LU$, the linear system $Ax=b$ is equivalent to $LUx=Pb$. Introduce the intermediate vector $y=Ux$: first use forward substitution to solve $Ly=Pb$, then back substitution to solve $Ux=y$. With unpivoted LU decomposition, $P=I$ is the identity matrix, so the first right-hand side is simply $b$.^[ref]

For the factors in Computed three-by-three LU decomposition and $b=(7,23,26{)}^T$, forward substitution givesThen back substitution givesDirect matrix multiplication verifies $A(1,2,3{)}^T=(7,23,26{)}^T$. More generally, the back substitution recurrence isThe direction matters: the last row contains only $x_n$, whereas the first row generally involves every unknown. No dense matrix inverse is needed for either triangular solve.

Forward substitution solves a linear system $Ly=c$ with a lower triangular matrix $L$ and nonzero diagonal, from the first row downward:The first sum is empty, so $y_1=c_1/l_{11}$. Each subsequent row uses only components already found. For a unit triangular matrix $L$, omit the division because $l_{ii}=1$. The complementary procedure, back substitution, solves an upper triangular matrix system from the last row upward.^[ref]

We fix the convention $PA=LU$: the permutation matrix $P$ acts on the rows of $A$, and the same row ordering must be applied to $b$ before the Triangular solves for LU decomposition. Equivalently, $A=P^{T}LU$, since $P^{-1}=P^T$. LAPACK documents its factors as $A=PLU$ instead, so check the defining equation rather than just the letter $P$.^[ref] Partial pivoting selects the largest absolute candidate in the current column, among rows $k,\dots ,n$, before forming elimination multipliers.^[ref]

A zero first pivot element does not imply a singular matrix. For example,Unpivoted Gaussian elimination would divide $2$ by zero. Indeed, no $A_0=LU$ with a unit lower triangular matrix $L$ exists: $u_{11}=0$ would force $(LU{)}_{21}=l_{21}{u}_{11}=0$, contradicting $(A_0{)}_{21}=2$. One row interchange resolves the obstruction. For $b_0=(1,5{)}^T$, use $P{b}_0=(5,1{)}^T$; back substitution then gives $x=(1,1{)}^T$. Using the original $b_0$ with these factors would solve a different linear system.

Partial pivoting bounds each elimination multiplier by one in absolute value in exact arithmetic. This helps control roundoff error, but intermediate entries can still grow; it is not an unconditional guarantee of numerical stability for every matrix.^[ref] It also does not remove the original problem's condition number. For the diagonal matrix $D=\mathrm{diag}(1,\epsilon )$ with $0<\epsilon <1$, the condition number in the Euclidean norm is ${\kappa }_2(D)=1/\epsilon$. Row permutations preserve this condition number, because they preserve Euclidean norms. Even an exact solve changes from $(1,1{)}^T$ to $(1,1+\eta /\epsilon {)}^T$ when the right-hand side changes from $(1,\epsilon {)}^T$ to $(1,\epsilon +\eta {)}^T$. This sensitivity is a property of the linear system, not a zero-pivot element failure.^[ref]

For a permutation matrix $P$ and elimination matrices $E_i$, put $F_i=P{E}_i{P}^{-1}$. Cancellation of adjacent $P^{-1}P$ gives $P{E}_k\cdots E_{1}M=F_k\cdots F_1(PM)$. If the $F_i$ and their inverses are lower triangular, then $PM=LU$ with $U=P{E}_k\cdots E_{1}M$ and $L=F_1^{-1}\cdots F_k^{-1}$. Lower triangularity needs checking: an arbitrary permutation does not preserve it.

When a fixed matrix $A$ is used with several forcing or measurement vectors, compute its Pivoted LU decomposition once. Collect $r$ right-hand sides into the columns of $B=(b^{(1)}\cdots b^{(r)})$, and the corresponding solutions into $X$. ThenEach column undergoes the same row permutation, forward substitution and back substitution. Neither $L$, $U$ nor $P$ depends on $B$. This is the useful separation between preparing a linear system and applying its solution procedure repeatedly. LAPACK implements it by taking previously computed factors, permuting all right-hand sides, and performing the two triangular solves.^[ref]

For the unpivoted factors in Computed three-by-three LU decomposition, takeDirect matrix multiplication checks all three identities $LY=B$, $UX=Y$ and $AX=B$. The first column reproduces the previous solve; the other two use exactly the same factors. New right-hand sides can also arrive one at a time. If $A$ changes, these stored factors generally need to be recomputed.

For classical dense Gaussian elimination over the real numbers, count each scalar addition, subtraction, multiplication or division as one arithmetic operation. At step $k$, there are $n-k$ elimination multipliers to divide out and $(n-k{)}^2$ remaining entries, each requiring one multiplication and one subtraction. Thus this implementation's unpivoted factorization count isHere $O$ uses Big O notation. Partial pivoting adds $O(n^2)$ comparisons and row movement, leaving the leading arithmetic count unchanged. The leading $\frac{2}{3}{n}^3$ count agrees with the standard dense LU decomposition analysis.^[ref]

With unit $L$, forward substitution uses $n(n-1)$ arithmetic operations. Back substitution uses another $n(n-1)$ plus $n$ diagonal divisions. Together these cost $2n^2-n$ per right-hand side; applying a stored row permutation moves $n$ entries and requires no scalar arithmetic. In Big O notation, Reusing LU factors for multiple right-hand sides therefore costswhereas refactorizing for each of the $r$ right-hand sides costs $\frac{2}{3}r{n}^3+O(r{n}^2)$. For $n=1000$ and $r=10$, the leading counts are about $6.87\times 10^8$ versus $6.69\times 10^9$ arithmetic operations, nearly a tenfold saving. These are operation counts, not measured runtimes: memory access and optimized matrix multiplication also affect speed.^[ref]