Chapter 75. LU Decomposition

LU decomposition is a factorization of a matrix into a lower triangular part and an upper triangular part. It is one of the standard matrix decompositions used to solve linear systems, compute determinants, and organize Gaussian elimination in matrix form.

For a square matrix $A$ , an LU decomposition has the form

A = LU,

where $L$ is lower triangular and $U$ is upper triangular. A lower triangular matrix has zero entries above the main diagonal. An upper triangular matrix has zero entries below the main diagonal. This is the standard lower-upper factorization described in numerical linear algebra references.

75.1 Triangular Matrices

A lower triangular matrix has the form

L = \begin{bmatrix} \ell_{11} & 0 & 0 & \cdots & 0 \\ \ell_{21} & \ell_{22} & 0 & \cdots & 0 \\ \ell_{31} & \ell_{32} & \ell_{33} & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \ell_{n1} & \ell_{n2} & \ell_{n3} & \cdots & \ell_{nn} \end{bmatrix}.

An upper triangular matrix has the form

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}.

Triangular systems are easy to solve because each equation contains only a controlled set of unknowns. A lower triangular system is solved from top to bottom. An upper triangular system is solved from bottom to top.

This is the reason LU decomposition is useful. It replaces one general linear system by two triangular systems.

75.2 The Basic Factorization

Let

A = \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{bmatrix}.

An LU decomposition writes $A$ as

A = LU.

In many treatments, $L$ is chosen to have ones on the diagonal:

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

Such a matrix is called unit lower triangular. With this convention, the diagonal scaling is placed in $U$ . This convention is often called the Doolittle form of LU decomposition. A related convention, called the Crout form, places the unit diagonal in $U$ instead.

75.3 LU and Gaussian Elimination

LU decomposition is Gaussian elimination recorded as a matrix factorization.

Suppose Gaussian elimination transforms $A$ into an upper triangular matrix $U$ . At each step, a multiple of one row is subtracted from a lower row to eliminate an entry below the pivot. The multipliers used in these eliminations become the entries of $L$ .

For example, to eliminate $a_{21}$ using $a_{11}$ , one uses the multiplier

\ell_{21} = \frac{a_{21}}{a_{11}}.

To eliminate $a_{31}$ , one uses

\ell_{31} = \frac{a_{31}}{a_{11}}.

These multipliers are stored below the diagonal in $L$ . The remaining matrix after elimination is $U$ .

Thus $U$ records the result of elimination, while $L$ records how the elimination was performed.

75.4 A Two by Two Example

Let

A = \begin{bmatrix} 2 & 3 \\ 4 & 7 \end{bmatrix}.

We seek

A = LU,

with

L = \begin{bmatrix} 1 & 0 \\ \ell_{21} & 1 \end{bmatrix}, \qquad U = \begin{bmatrix} u_{11} & u_{12} \\ 0 & u_{22} \end{bmatrix}.

Multiplying gives

LU = \begin{bmatrix} u_{11} & u_{12} \\ \ell_{21}u_{11} & \ell_{21}u_{12} + u_{22} \end{bmatrix}.

Equating this with $A$ , we obtain

u_{11} = 2, \qquad u_{12} = 3,

\ell_{21}u_{11} = 4,

\ell_{21} = 2.

The final entry gives

\ell_{21}u_{12} + u_{22} = 7.

Therefore

2 \cdot 3 + u_{22} = 7,

u_{22} = 1.

Hence

L = \begin{bmatrix} 1 & 0 \\ 2 & 1 \end{bmatrix}, \qquad U = \begin{bmatrix} 2 & 3 \\ 0 & 1 \end{bmatrix}.

Indeed,

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

75.5 Solving Linear Systems with LU

Suppose

Ax = b

and

A = LU.

Then

LUx = b.

Introduce an intermediate vector

y = Ux.

Then the original system becomes two triangular systems:

Ly = b

and

Ux = y.

The first system is solved by forward substitution. The second system is solved by backward substitution.

This procedure is especially useful when many systems have the same coefficient matrix $A$ but different right-hand sides $b$ . The factorization $A = LU$ is computed once. Each new right-hand side then requires only two triangular solves.

75.6 Forward Substitution

Let

Ly = b,

where $L$ is lower triangular. Written componentwise,

\ell_{11}y_1 = b_1,

\ell_{21}y_1 + \ell_{22}y_2 = b_2,

\ell_{31}y_1 + \ell_{32}y_2 + \ell_{33}y_3 = b_3,

and so on.

The first equation gives $y_1$ . The second equation then gives $y_2$ . The third equation gives $y_3$ . Continuing in this order gives all entries of $y$ .

For a general lower triangular matrix,

y_i = \frac{ b_i - \sum_{j=1}^{i-1} \ell_{ij}y_j }{ \ell_{ii} }.

If $L$ is unit lower triangular, then $\ell_{ii} = 1$ , and the formula becomes

y_i = b_i - \sum_{j=1}^{i-1} \ell_{ij}y_j.

75.7 Backward Substitution

Let

Ux = y,

where $U$ is upper triangular. Written componentwise, the last equation is

u_{nn}x_n = y_n.

This gives $x_n$ . The previous equation gives $x_{n-1}$ , since $x_n$ is already known. The process continues upward.

For a general upper triangular matrix,

x_i = \frac{ y_i - \sum_{j=i+1}^{n} u_{ij}x_j }{ u_{ii} }.

Backward substitution works only when the diagonal entries needed for division are nonzero. In the LU setting, this means the pivots must be nonzero.

75.8 Example: Solving a System

Let

A = \begin{bmatrix} 2 & 3 \\ 4 & 7 \end{bmatrix}, \qquad b = \begin{bmatrix} 8 \\ 18 \end{bmatrix}.

From the previous example,

A = LU,

where

L = \begin{bmatrix} 1 & 0 \\ 2 & 1 \end{bmatrix}, \qquad U = \begin{bmatrix} 2 & 3 \\ 0 & 1 \end{bmatrix}.

First solve

Ly = b.

That is,

\begin{bmatrix} 1 & 0 \\ 2 & 1 \end{bmatrix} \begin{bmatrix} y_1 \\ y_2 \end{bmatrix} = \begin{bmatrix} 8 \\ 18 \end{bmatrix}.

The equations are

y_1 = 8,

2y_1 + y_2 = 18.

Therefore

y_2 = 18 - 16 = 2.

y = \begin{bmatrix} 8 \\ 2 \end{bmatrix}.

Now solve

Ux = y.

That is,

\begin{bmatrix} 2 & 3 \\ 0 & 1 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} 8 \\ 2 \end{bmatrix}.

The second equation gives

x_2 = 2.

The first equation gives

2x_1 + 3x_2 = 8.

Thus

2x_1 + 6 = 8,

x_1 = 1.

The solution is

x = \begin{bmatrix} 1 \\ 2 \end{bmatrix}.

75.9 Existence of LU Decomposition

An LU decomposition without row exchanges does not always exist. The elimination process requires nonzero pivots. If a pivot is zero, division by that pivot cannot be performed.

For example,

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

cannot be decomposed as $A = LU$ with $L$ unit lower triangular and $U$ upper triangular. The first pivot would be zero.

However, row exchanges fix this problem. Instead of factoring $A$ directly, one factors a row-permuted version of $A$ :

PA = LU.

Here $P$ is a permutation matrix. This form is called LU decomposition with pivoting, or LUP decomposition. Standard numerical algorithms usually include pivoting because it avoids zero pivots and improves numerical stability.

75.10 Pivoting

Pivoting means exchanging rows before or during elimination.

The most common form is partial pivoting. At each step, the algorithm chooses a pivot row by looking down the current column and selecting an entry with large absolute value. That row is swapped into the pivot position.

The factorization then has the form

PA = LU.

Equivalently,

A = P^{-1}LU.

Since a permutation matrix satisfies

P^{-1} = P^T,

one may also write

A = P^TLU.

Pivoting has two purposes. First, it prevents division by zero when a possible nonzero pivot is available. Second, it can reduce the growth of rounding errors in floating point arithmetic.

75.11 Determinants from LU

LU decomposition also simplifies determinant computation.

A = LU,

then

\det(A) = \det(L)\det(U).

The determinant of a triangular matrix is the product of its diagonal entries. Therefore,

\det(A) = \left(\prod_{i=1}^{n} \ell_{ii}\right) \left(\prod_{i=1}^{n} u_{ii}\right).

If $L$ is unit lower triangular, then

\det(L) = 1.

Hence

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

If pivoting is used and

PA = LU,

then

\det(P)\det(A) = \det(L)\det(U).

Thus

\det(A) = \det(P)\det(L)\det(U),

because $\det(P) = \pm 1$ . Each row swap changes the sign of the determinant.

75.12 Inverses from LU

LU decomposition can be used to compute the inverse of a matrix, although direct inversion is often avoided in numerical work when only a solution of $Ax=b$ is needed.

To find $A^{-1}$ , solve

AX = I.

A = LU,

then

LUX = I.

This is equivalent to solving

LY = I

and then

UX = Y.

Each column of $X$ is found by solving two triangular systems. The columns of $X$ are the columns of $A^{-1}$ .

This method is systematic, but it is usually better to solve $Ax=b$ directly by LU rather than compute $A^{-1}$ explicitly.

75.13 Algorithm

The following algorithm computes an LU decomposition without pivoting, using the convention that $L$ has ones on the diagonal.

lu_decomposition(A):
    n = number of rows of A

    L = identity matrix of size n
    U = zero matrix of size n by n

    for k = 1 to n:
        for j = k to n:
            U[k,j] = A[k,j]
            for s = 1 to k-1:
                U[k,j] = U[k,j] - L[k,s] * U[s,j]

        for i = k+1 to n:
            L[i,k] = A[i,k]
            for s = 1 to k-1:
                L[i,k] = L[i,k] - L[i,s] * U[s,k]
            L[i,k] = L[i,k] / U[k,k]

    return L, U

The division by $U[k,k]$ requires a nonzero pivot. If the pivot is zero or very small, the algorithm should use pivoting.

75.14 Operation Count

For a dense $n \times n$ matrix, LU decomposition requires on the order of

\frac{2}{3}n^3

arithmetic operations.

Solving the two triangular systems requires on the order of

2n^2

operations for each right-hand side.

Thus the expensive step is the factorization. Once the factorization is available, solving for additional right-hand sides is much cheaper.

This explains why LU decomposition is useful in repeated-solve problems.

75.15 Interpretation

LU decomposition separates elimination into two pieces.

The matrix $U$ is the upper triangular matrix obtained after elimination. The matrix $L$ stores the multipliers that produced the eliminations. The equation

A = LU

says that the original matrix can be reconstructed from these two pieces.

In computational terms, LU decomposition is a reusable form of Gaussian elimination. In algebraic terms, it is a factorization into simpler matrices. In geometric terms, it decomposes a linear transformation into transformations with triangular structure.

75.16 Common Variants

Variant	Form	Description
LU	$A = LU$	Factorization without pivoting
LUP	$PA = LU$	LU with row pivoting
PLU	$A = PLU$	Equivalent convention using a permutation factor
LDU	$A = LDU$	Separates diagonal scaling into $D$
Crout	$A = LU$	$U$ has unit diagonal
Doolittle	$A = LU$	$L$ has unit diagonal

Different conventions place diagonal entries in different factors. The essential idea remains the same: a general matrix is represented by triangular factors.

75.17 Summary

LU decomposition factors a matrix into lower and upper triangular matrices:

A = LU.

With pivoting, the more stable and general form is

PA = LU.

The decomposition is Gaussian elimination written as a matrix product. It is used to solve linear systems, compute determinants, and perform repeated solves efficiently. Its practical value comes from replacing one general system by two triangular systems:

Ly = b, \qquad Ux = y.

The triangular systems are solved by forward and backward substitution. This makes LU decomposition one of the central tools of computational linear algebra.