LU Decomposition

LU decomposition, or LU factorization, stores Gaussian elimination as a product of triangular matrices. In the unpivoted form, $A=LU$, where $L$ is lower triangular and $U$ is upper triangular; with row pivoting, the convention used here is $PA=LU$. The factors support linear-system solves, determinant calculations, and inverse calculations without repeating elimination from the beginning. To solve $A\mathbf x=\mathbf b$, first solve a lower-triangular system for an intermediate vector and then an upper-triangular system for $\mathbf x$. These substitution steps are inexpensive compared with the initial factorization, so retaining the factors is especially useful when the same matrix is paired with several right-hand sides.

Mathematical formulation

Consider a square $n \times n$ matrix $A$:

$$ A = \begin{bmatrix} a_{11} & a_{12} & a_{13} & \cdots & a_{1n} \\ a_{21} & a_{22} & a_{23} & \cdots & a_{2n} \\ a_{31} & a_{32} & a_{33} & \cdots & a_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & a_{n3} & \cdots & a_{nn} \end{bmatrix} $$

The LU decomposition expresses $A$ as:

$$ A = LU $$

where

$$ L = \begin{bmatrix} 1 & 0 & 0 & \cdots & 0 \\ l_{21} & 1 & 0 & \cdots & 0 \\ l_{31} & l_{32} & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ l_{n1} & l_{n2} & l_{n3} & \cdots & 1 \end{bmatrix} $$

$$ U = \begin{bmatrix} u_{11} & u_{12} & u_{13} & \cdots & u_{1n} \\ 0 & u_{22} & u_{23} & \cdots & u_{2n} \\ 0 & 0 & u_{33} & \cdots & u_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & u_{nn} \end{bmatrix} $$

Here, $L$ is lower-triangular with ones on the diagonal, and $U$ is upper-triangular. The factorization is not always guaranteed to exist unless certain conditions are met (e.g., all leading principal minors of $A$ being nonzero when no pivoting is used; with partial pivoting, a factorization $PA = LU$ exists for every nonsingular $A$).

The formulas below use the Doolittle convention, which fixes every diagonal entry of $L$ to one. Specifying that diagonal removes a scaling ambiguity between the two triangular factors and makes the stored elimination multipliers easier to interpret.

Derivation

The derivation of the LU decomposition closely follows the steps of Gaussian elimination. Gaussian elimination transforms the matrix $A$ into an upper-triangular matrix by adding multiples of one row to another. These multipliers can be stored in the entries of a lower-triangular matrix $L$.

Starting from:

$$ A\mathbf{x} = \mathbf{b} $$

we write $A = LU$. Substitute to get:

$$ (LU)\mathbf{x} = \mathbf{b} $$

Introducing an intermediate vector $\mathbf{c}$:

$$ U\mathbf{x} = \mathbf{c} \implies L\mathbf{c} = \mathbf{b} $$

Since $L$ is lower-triangular and nonsingular (with ones on its diagonal), we can quickly solve for $\mathbf{c}$ using forward substitution. Once $\mathbf{c}$ is known, we solve the upper-triangular system $U\mathbf{x} = \mathbf{c}$ via backward substitution.

The process of determining $L$ and $U$ essentially mimics the elimination steps:

I. Use the first row of $A$ to eliminate entries below $a_{11}$.

II. Store these elimination factors in $L$.

III. After the first column is dealt with, the submatrix of $A$ (excluding the first row and column) is similarly factorized.

IV. This process continues recursively until $A$ is fully decomposed into $L$ and $U$.

Recovering the factor formulas

The $(i,j)$ entry of $A=LU$ is $a_{ij}=\sum_k l_{ik}u_{kj}$. Triangularity removes most terms. For $i\leq j$, only $k\leq i$ remain, and $l_{ii}=1$, giving

$$ a_{ij}=\sum_{k=1}^{i-1}l_{ik}u_{kj}+u_{ij}. $$

Solving this identity for $u_{ij}$ gives the upper-row formula below. For $j>i$, the $(j,i)$ entry similarly satisfies

$$ a_{ji}=\sum_{k=1}^{i-1}l_{jk}u_{ki}+l_{ji}u_{ii}, $$

which gives the lower-column formula after division by $u_{ii}$. These identities explain why earlier columns must be available before later columns are computed.

Row swaps change the factorization to $PA=LU$. The same permutation must act on $\mathbf{b}$, and a later swap must also exchange the previously stored multipliers in the relevant rows of $L$. Mixing the conventions $PA=LU$ and $A=PLU$ is a common source of incorrect solves; follow the convention returned by the implementation.

Dense factorization costs $O(n^3)$ once, followed by $O(n^2)$ for each new right-hand side. Reusing the factors therefore avoids repeating elimination and usually avoids the extra work of forming an explicit inverse. In the example below, exact arithmetic gives $(x_1,x_2,x_3)=(-5/23,29/69,-1/23)$, which can be substituted into the original equations to verify the rounded results.

Algorithm steps

Given an $n \times n$ matrix $A$, the LU decomposition algorithm without pivoting can be described as follows:

Step 1: Initialization.

Set $L = I$ (the $n \times n$ identity matrix) and $U = 0$ (the $n \times n$ zero matrix).

Step 2: Main loop. (for $i = 1$ to $n$):

Compute the diagonal and upper elements of $U$:

For $j = i$ to $n$:

$$ u_{ij} = a_{ij} - \sum_{k=1}^{i-1} l_{ik} u_{kj}. $$

Compute the lower elements of $L$:

For $j = i+1$ to $n$:

$$ l_{j i} = \frac{1}{u_{ii}}\left(a_{j i} - \sum_{k=1}^{i-1} l_{jk} u_{k i}\right). $$

III. After these loops complete, $A = LU$ is obtained.

Step 4: Solving $A\mathbf{x} = \mathbf{b}$.

Forward substitution for $L\mathbf{c} = \mathbf{b}$:

For $i = 1$ to $n$:

$$ c_i = b_i - \sum_{k=1}^{i-1} l_{ik} c_{k}. $$

Backward substitution for $U\mathbf{x} = \mathbf{c}$:

For $i = n$ down to $1$:

$$ x_i = \frac{c_i - \sum_{k=i+1}^{n} u_{ik}x_{k}}{u_{ii}}. $$

Step 5: Partial pivoting (what the implementation does).

The steps above assume every pivot $u_{ii}$ is nonzero. The implementation instead performs partial pivoting: at step $i$ it swaps row $i$ with the row $p \geq i$ whose entry in column $i$ has the largest absolute value, recording the swap in a permutation matrix $P$ (and in the already computed columns of $L$), and then eliminates as before. The result is the factorization

$$ PA = LU, $$

with $L$ unit lower-triangular, $U$ upper-triangular and $P$ a permutation matrix, returned as (P, L, U). If the largest available pivot is not larger than $n\,\varepsilon_{\text{mach}}\max_{ij}|a_{ij}|$ the matrix is treated as singular and an error is raised.

Solving $A\mathbf{x}=\mathbf{b}$ then reads $PA\mathbf{x} = P\mathbf{b}$, so forward substitution is applied to $P\mathbf{b}$:

$$ L\mathbf{c} = P\mathbf{b}, \qquad U\mathbf{x} = \mathbf{c}. $$

The determinant follows from $\det(P)\det(A) = \det(L)\det(U)$. Since $\det(L)=1$ and $\det(P) = (-1)^s$ where $s$ is the number of row swaps,

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

(Without pivoting, $s=0$ and this reduces to $\det(A)=\prod u_{ii}$.)

Example

The example below is worked without pivoting (so $P = I$) to keep the arithmetic simple. The implementation, with partial pivoting, would first swap rows 1 and 2 (the entry $3$ is the largest in the first column) and so produce a different, equally valid factorization $PA = LU$; the solution $\mathbf{x}$ is the same.

Consider the system of equations:

$$ \begin{aligned} 2x + 3y - 4z &= 1, \\ 3x - 3y + 2z &= -2, \\ -2x + 6y - z &= 3. \end{aligned} $$

In matrix form:

$$ A = \begin{bmatrix} 2 & 3 & -4 \\ 3 & -3 & 2 \\ -2 & 6 & -1 \end{bmatrix} $$

$$ \mathbf{b} = \begin{bmatrix} 1 \\ -2 \\ 3 \end{bmatrix} $$

Step-by-Step LU decomposition

Step 1: Initialize

Set:

$$ L = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} $$

$$ U = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} $$

Compute First Row of $U$:

$$ u_{11} = a_{11} = 2, \quad u_{12} = a_{12} = 3, \quad u_{13} = a_{13} = -4. $$

Thus:

$$ U = \begin{bmatrix} 2 & 3 & -4 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} $$

Compute First Column of $L$ (below diagonal):

For $i = 2$ to 3:

$$ l_{21} = \frac{a_{21}}{u_{11}} = \frac{3}{2} = 1.5, \quad l_{31} = \frac{a_{31}}{u_{11}} = \frac{-2}{2} = -1. $$

Now:

$$ L = \begin{bmatrix} 1 & 0 & 0 \\ 1.5 & 1 & 0 \\ -1 & 0 & 1 \end{bmatrix} $$

Second Pivot (i=2):

Compute $u_{22}$:

$$ u_{22} = a_{22} - l_{21} u_{12} = (-3) - (1.5)(3) = -3 - 4.5 = -7.5. $$

Compute $u_{23}$:

$$ u_{23} = a_{23} - l_{21}u_{13} = 2 - (1.5)(-4) = 2 + 6 = 8. $$

Thus:

$$ U = \begin{bmatrix} 2 & 3 & -4 \\ 0 & -7.5 & 8 \\ 0 & 0 & 0 \end{bmatrix} $$

For $l_{32}$:

$$ l_{32} = \frac{a_{32} - l_{31}u_{12}}{u_{22}} = \frac{6 - (-1)(3)}{-7.5} = \frac{6 + 3}{-7.5} = \frac{9}{-7.5} = -1.2. $$

Update $L$:

$$ L = \begin{bmatrix} 1 & 0 & 0 \\ 1.5 & 1 & 0 \\ -1 & -1.2 & 1 \end{bmatrix} $$

Third Pivot (i=3):

Compute $u_{33}$:

$$ u_{33} = a_{33} - l_{31}u_{13} - l_{32}u_{23} = (-1) - (-1)(-4) - (-1.2)(8). $$

Carefully evaluate:

$$ (-1) - (-1 \times - 4) - (-1.2 \times 8) = (-1) - (4) - (-9.6) = -5 + 9.6 = 4.6. $$

Thus:

$$ U = \begin{bmatrix} 2 & 3 & -4 \\ 0 & -7.5 & 8 \\ 0 & 0 & 4.6 \end{bmatrix} $$

So finally, we have:

$$ L = \begin{bmatrix} 1 & 0 & 0 \\ 1.5 & 1 & 0 \\ -1 & -1.2 & 1 \end{bmatrix} $$

$$ U = \begin{bmatrix} 2 & 3 & -4 \\ 0 & -7.5 & 8 \\ 0 & 0 & 4.6 \end{bmatrix} $$

Forward Substitution ($L\mathbf{c} = \mathbf{b}$):

$$ c_1 = b_1 = 1 $$

$$ c_2 = b_2 - l_{21} c_1 = -2 - (1.5)(1) = -3.5 $$

$$ c_3 = b_3 - l_{31}c_1 - l_{32}c_2 = 3 - (-1)(1) - (-1.2)(-3.5) = 3 + 1 - 4.2 = -0.2 $$

Backward Substitution ($U\mathbf{x} = \mathbf{c}$):

$$ x_3 = \frac{c_3}{u_{33}} = \frac{-0.2}{4.6} \approx - 0.0434783 $$

$$ x_2 = \frac{c_2 - u_{23}x_3}{u_{22}} = \frac{-3.5 - (8)(-0.0434783)}{-7.5} = \frac{-3.5 + 0.3478264}{-7.5} = \frac{-3.1521736}{-7.5} \approx 0.42029 $$

$$ x_1 = \frac{c_1 - u_{12}x_2 - u_{13}x_3}{u_{11}} = \frac{1 - 3(0.42029) - (-4)(-0.0434783)}{2} = \frac{1 - 1.26087 - 0.1739132}{2} = \frac{-0.4347832}{2} = -0.2173916 $$

Advantages

• Once $A = LU$ is computed, solving multiple systems $A\mathbf{x} = \mathbf{b}$ becomes efficient, as only forward and backward substitution are required for each new right-hand-side vector $\mathbf{b}$. This is particularly beneficial in applications requiring repeated solves with the same matrix $A$.
• LU decomposition organizes the elimination steps into the matrices $L$ (lower triangular) and $U$ (upper triangular), simplifying the process and providing a structured representation of the system. Partial pivoting can be incorporated, enhancing numerical stability for a wide range of problems.
• It allows for the efficient computation of matrix determinants (via $\det(A) = (-1)^{s}\prod u_{ii}$, where $s$ is the number of row swaps in $PA = LU$; $s=0$ without pivoting) and matrix inverses, and serves as a building block for advanced numerical techniques, such as eigenvalue computations and solving partial differential equations.

Limitations

• Not all matrices are directly LU decomposable without row interchanges. For many practical cases, partial pivoting is required, resulting in a decomposition of the form $PA = LU$, where $P$ is a permutation matrix.
• If the matrix $A$ has zero diagonal elements or does not meet certain structural conditions, direct LU decomposition without permutations may fail or lead to numerical instability.
• For large sparse matrices, naive LU decomposition can cause significant fill-in, where new nonzero elements appear in $L$ and $U$. This increases both memory usage and computational complexity, potentially rendering the method impractical without specialized sparse matrix techniques.