---
title: LU and Cholesky Factorization in Practice
module: Numerical Linear Algebra
moduleNumber: 8
lessonNumber: 2
order: 802
summary: >
  Gaussian elimination, read as a factorization A = LU, turns a linear system
  into two triangular solves. A single near-zero pivot wrecks it, so partial
  pivoting reorders rows to pick the largest available pivot and makes the method
  work for every invertible matrix. For symmetric positive-definite systems,
  Cholesky halves the cost and needs no pivoting.
topics: [Numerical Linear Algebra]
sources:
  - book: Bornemann
    ref: "Part II — Matrix Factorization; §7 Triangular Decomposition"
  - book: Bornemann
    ref: "§8 Cholesky Decomposition"
  - book: Lay
    ref: "§2.5 Matrix Factorizations (LU bridge)"
draft: false
---

Solving $Ax = b$ by hand means Gaussian elimination: clear the first column, clear
the second, and so on until the matrix is triangular, then back-substitute. Von
Neumann and Goldstine observed in 1947 that this is really two separate operations,
first _decompose_ $A$ into a product of triangular matrices, then _solve_ by
substitution, and that separating them is what lets one factorization serve any
number of right-hand sides.[^bor-lu] The decomposition is $A = LU$, and once it is
in hand the system reduces to two [triangular
solves](/linear-algebra/numerical-linear-algebra/numerical-thinking-and-matrix-computation).

## Reading elimination as a factorization

> **Definition (Triangular decomposition).** A **triangular** (or **LU**)
> **decomposition** of $A \in \mathrm{GL}(m;\mathbb{K})$ is a factorization
> $A = LU$ with $L$ lower triangular and $U$ upper triangular. It is
> **normalized** when $L$ is unipotent (unit diagonal).

To solve $Ax = b$, substitute $A = LU$ and split it in two:

$$
L z = b \quad(\text{forward substitution}),
\qquad
U x = z \quad(\text{back substitution}),
$$

so that $L(Ux) = Ax = b$. When a normalized decomposition exists, both factors are
unique: if $L_1 U_1 = L_2 U_2$, then $L_2^{-1} L_1 = U_2 U_1^{-1}$ is at once
unipotent-lower and upper triangular, forcing it to equal $I$.[^bor-lu]

The factors are built one row and column at a time by a recursion that peels off
the first row and column at each step. Partition

$$
A_k = \begin{bmatrix} \alpha_k & u_k^\ast \\ b_k & B_k \end{bmatrix}
    = \begin{bmatrix} 1 & \\ l_k & L_{k+1} \end{bmatrix}
      \begin{bmatrix} \alpha_k & u_k^\ast \\ & U_{k+1} \end{bmatrix}.
$$

Multiplying out the second block row gives $b_k = l_k \alpha_k$ and
$B_k = l_k u_k^\ast + L_{k+1}U_{k+1}$, so provided the **pivot** $\alpha_k \ne 0$,

$$
l_k = \frac{b_k}{\alpha_k},
\qquad
A_{k+1} = B_k - l_k u_k^\ast.
$$

The updated block $A_{k+1}$, the **Schur complement** of $\alpha_k$, is the input
to the next step. The pivots $\alpha_1, \dots, \alpha_m$ end up on the diagonal of
$U$, which yields the exact condition for existence.

> **Lemma (Existence of $LU$).** $A \in \mathrm{GL}(m;\mathbb{K})$ has a normalized
> triangular decomposition if and only if every pivot $\alpha_1, \dots, \alpha_m$
> is nonzero — equivalently, if and only if every leading principal submatrix
> $A_k$ is invertible.

$$
% caption: One step of $LU$: the pivot $\alpha_k$ and first row/column are read
% off, then the rank-one update $-l_k u_k^\ast$ overwrites the trailing Schur
% complement, shrinking the active block by one.
\begin{tikzpicture}[scale=1.0, font=\footnotesize]
\definecolor{acc}{HTML}{4A6FA5}
% A_k block
\draw[black] (0,0) rectangle (2.4,2.4);
\draw[black] (0.6,0)--(0.6,2.4); \draw[black] (0,1.8)--(2.4,1.8);
\node[acc] at (0.3,2.1) {$a_{kk}$};
\node at (1.5,2.1) {$u_k^\ast$};
\node at (0.3,0.9) {$b_k$};
\node at (1.5,0.9) {$B_k$};
\node at (1.2,-0.4) {$A_k$};
\draw[->, black, thick] (2.9,1.2) -- (3.9,1.2);
% result: l_k stored, Schur complement
\draw[black] (4.4,0) rectangle (6.8,2.4);
\draw[black] (5.0,0)--(5.0,2.4); \draw[black] (4.4,1.8)--(6.8,1.8);
\node[acc] at (4.7,2.1) {$a_{kk}$};
\node at (5.9,2.1) {$u_k^\ast$};
\node[acc] at (4.7,0.9) {$l_k$};
\node[acc] at (5.9,0.9) {$A_{k+1}$};
\node[acc, align=center] at (5.9,-0.55) {$A_{k+1}=B_k-l_k u_k^\ast$};
\end{tikzpicture}
$$

The cost is dominated by the rank-one update, $2(m-k)^2$ flops at step $k$, so the
full factorization costs

$$
\#\text{flop} \approx 2\sum_{k=1}^{m}(m-k)^2 \approx \frac{2}{3}m^3,
$$

and the two substitutions that follow add only $2m^2$. The factorization is the
expensive part; once $A = LU$ is stored, any number of right-hand sides is cheap.
Because $A$ is read once and overwritten in place by $L$ and $U$, the efficiency
ratio is $q \approx m/3$, so a Level-3 BLAS implementation runs near peak
performance.

## Small pivots and instability

Existence fails the instant a pivot is zero. The invertible matrix

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

has $\alpha_1 = 0$ and no normalized $LU$ at all. The deeper problem is that a
pivot close to zero is nearly as damaging as an exact zero, once arithmetic is
rounded. Take

$$
A = \begin{bmatrix} 10^{-20} & 1 \\ 1 & 1 \end{bmatrix}
  = \begin{bmatrix} 1 & 0 \\ 10^{20} & 1 \end{bmatrix}
    \begin{bmatrix} 10^{-20} & 1 \\ 0 & 1 - 10^{20} \end{bmatrix}.
$$

On a machine carrying sixteen significant digits, $1 - 10^{20}$ rounds to
$-10^{20}$: the subtraction of $1$ falls below the resolution of the
representation. The computed factor is the exact factorization of a _different_
matrix, one with the $(2,2)$ entry silently changed from $1$ to $0$. Solving
$Ax = e_1$ with the rounded factors returns $\hat{x} = (0, 1)^\top$ against the true
$x \approx (-1, 1)^\top$: a one-hundred-percent error in the first component from a
matrix that is perfectly well-behaved.

> **The Zeroth Law of Numerical Analysis.** If theoretical analysis breaks down at
> $\alpha = 0$, numerical analysis breaks down at $\alpha \approx 0$. A tiny pivot
> is a numerical zero.[^bor-lu]

## Partial pivoting

The fix uses a freedom in the problem: the row order of a system is
arbitrary. Before eliminating column $k$, scan it and move the entry of largest
absolute value into the pivot position.

> **Strategy (Partial pivoting).** At step $k$, choose the pivot to be the
> largest-magnitude entry in the first column of the active block $A_k$, and swap
> its row into the pivot row. All later multipliers then satisfy
> $\lvert l_k \rvert \le 1$.

Each row swap is a transposition $\tau_k$, and their accumulated product is a
permutation matrix $P$. The arithmetic of each step is identical to the
unpivoted recursion; only the row selection changes. The result is a factorization
of the row-permuted matrix.

> **Theorem (Pivoted $LU$).** Every $A \in \mathrm{GL}(m;\mathbb{K})$ admits a
> factorization
>
> $$
> P^\ast A = LU, \qquad \lvert L \rvert \le 1,
> $$
>
> with $P$ a permutation matrix, $L$ unipotent lower triangular, and $U$
> invertible upper triangular. Every pivot is nonzero.

The proof is an induction showing each Schur complement $A_{k+1}$ stays invertible,
so its first column is never entirely zero and a nonzero pivot always exists. The
bound $\lvert l_k \rvert \le 1$ is the payoff of choosing the largest pivot: the
multipliers can never blow up, which is what keeps the rounded factorization close
to a true one.

$$
% caption: Partial pivoting scans column $k$ for the largest-magnitude entry
% (here $8$, below the diagonal) and swaps its row up before eliminating, keeping
% every multiplier $\le 1$ in size.
\begin{tikzpicture}[scale=1.0, font=\footnotesize]
\definecolor{acc}{HTML}{4A6FA5}
% column of candidates
\draw[black] (0,0) rectangle (1.0,3.0);
\node at (0.5,2.6) {$0.3$};
\node at (0.5,1.8) {$2$};
\node[acc] at (0.5,1.0) {$8$};
\node at (0.5,0.2) {$-1$};
\draw[acc, thick] (0,0.6) rectangle (1.0,1.4);
\node[acc, anchor=west] at (1.2,1.0) {largest $\lvert\cdot\rvert$ in column};
% swap arrow
\draw[->, acc, thick] (4.7,2.6) .. controls (5.4,1.9) and (5.4,1.7) .. (4.7,1.0);
\draw[->, acc, thick] (4.7,1.0) .. controls (4.1,1.7) and (4.1,1.9) .. (4.7,2.6);
% after swap
\draw[black] (6.0,0) rectangle (7.0,3.0);
\node[acc] at (6.5,2.6) {$8$};
\node at (6.5,1.8) {$2$};
\node at (6.5,1.0) {$0.3$};
\node at (6.5,0.2) {$-1$};
\draw[acc, thick] (6.0,2.2) rectangle (7.0,3.0);
\node[anchor=north, align=center] at (0.5,-0.15) {column $k$};
\node[anchor=north, align=center] at (6.5,-0.15) {pivot on top};
\end{tikzpicture}
$$

Nothing about the solve changes: $Ax = b$ becomes $P^\ast A x = P^\ast b$, so one
permutes the right-hand side and runs the same forward and back substitution.

```algorithm
caption: $\textsc{PivotedLU}(A)$ — in-place $P^\ast A = LU$ with partial pivoting
$p \gets [1, 2, \dots, m]$ // permutation record
for $k = 1$ to $m-1$ do
  choose $j \ge k$ maximizing $\lvert A_{jk} \rvert$ // pivot search
  swap rows $k$ and $j$ of $A$ and entries $k, j$ of $p$
  for $i = k+1$ to $m$ do
    $A_{ik} \gets A_{ik} / A_{kk}$ // multiplier $l$, stored below the diagonal
  end for
  $A_{k+1:m,\, k+1:m} \gets A_{k+1:m,\, k+1:m} - A_{k+1:m,\, k}\, A_{k,\, k+1:m}$ // Schur update
end for
return $L = \text{unit lower part},\; U = \text{upper part},\; p$
```

> **Worked example.** Solve $Ax = b$ in three-significant-digit arithmetic with
>
> $$
> A = \begin{bmatrix} 10^{-4} & 1 \\ 1 & 1 \end{bmatrix},
> \qquad
> b = \begin{bmatrix} 1 \\ 2 \end{bmatrix},
> \qquad
> x_{\text{true}} = \begin{bmatrix} 1.0001 \\ 0.9999 \end{bmatrix}.
> $$
>
> _Without pivoting_, the multiplier is $l = 1 / 10^{-4} = 10^{4}$. Eliminating
> the second row,
>
> $$
> a_{22} = 1 - 10^{4}\cdot 1 = -9999 \to -1.00\times 10^{4},
> \qquad
> b_2 = 2 - 10^{4}\cdot 1 = -9998 \to -1.00\times 10^{4},
> $$
>
> where rounding to three digits discards the trailing $1$ and $2$. Back
> substitution gives $\xi_2 = (-1.00\times 10^4)/(-1.00\times 10^4) = 1.00$ and
> then $\xi_1 = (1 - 1.00)/10^{-4} = 0$. The computed $\hat x = (0,\, 1.00)^\top$
> is wrong in its first component by $100\%$, on a matrix that is well-behaved.
>
> _With partial pivoting_, $\lvert 1\rvert > \lvert 10^{-4}\rvert$ forces a row
> swap. Now $l = 10^{-4}/1 = 10^{-4}$, and
>
> $$
> a_{22} = 1 - 10^{-4}\cdot 1 = 0.9999 \to 1.00,
> \qquad
> b_2 = 1 - 10^{-4}\cdot 2 = 0.9998 \to 1.00,
> $$
>
> so $\xi_2 = 1.00$ and $\xi_1 = (2 - 1.00)/1 = 1.00$. The computed $\hat x =
> (1.00,\, 1.00)^\top$ matches the true solution to all three digits carried.

$$
% caption: The two computed solutions of the near-singular system against the true
% solution. Without pivoting the first component collapses to zero; partial
% pivoting lands on the true point.
\begin{tikzpicture}[scale=1.0, font=\footnotesize]
\definecolor{acc}{HTML}{4A6FA5}
\draw[->, black] (-0.3,0) -- (6.2,0) node[anchor=west] {$x_1$};
\draw[->, black] (0,-0.3) -- (0,2.9) node[anchor=south] {$x_2$};
\foreach \i/\x in {0/0, 1/3, 2/6} {\draw[black] (\x,-0.08)--(\x,0.08); \node[anchor=north] at (\x,-0.1) {$\i$};}
\draw[black] (-0.08,2)--(0.08,2); \node[anchor=east] at (-0.1,2) {$1$};
\draw[black, dashed] (0,2) -- (3,2);
\node[anchor=north] at (1.5,1.78) {component 1 lost};
\fill[acc] (3,2) circle (2.4pt);
\node[acc, anchor=south west] at (3.1,2.05) {true and pivoted};
\fill[black] (0,2) circle (2.4pt);
\node[anchor=west] at (0.15,2.5) {no pivot};
\end{tikzpicture}
$$

Partial pivoting is the default dense solver everywhere; it is what a call to
`A\b` or LAPACK's `xGETRF` runs underneath.

## Cholesky for symmetric positive-definite systems

A large class of systems — least-squares normal equations, covariance matrices,
stiffness matrices in mechanics — carry extra structure that a general solver
wastes. These are the symmetric (Hermitian) positive-definite matrices.

> **Definition (s.p.d.).** $A \in \mathbb{K}^{m \times m}$ is **self-adjoint** if
> $A^\ast = A$ (symmetric for $\mathbb{K}=\mathbb{R}$, Hermitian for
> $\mathbb{K}=\mathbb{C}$), and **positive definite** if $x^\ast A x > 0$ for all
> $x \ne 0$. A self-adjoint positive-definite matrix is **s.p.d.** For such $A$
> the quantity $x^\ast A x$ is always real, and $A$ is invertible.

For s.p.d. matrices the factorization can be made symmetric, and no pivoting is
ever needed because every leading principal submatrix inherits the s.p.d.
property and so is invertible.

> **Theorem (Cholesky).** Every s.p.d. matrix $A$ has a unique factorization
>
> $$
> A = L L^\ast,
> $$
>
> with $L$ lower triangular and a strictly positive diagonal. (Equivalently
> $A = R^\ast R$ with $R = L^\ast$ upper triangular.)

The factor is built row by row. Partitioning $A_k = L_k L_k^\ast$ and matching
blocks gives, for row $k$,

$$
l_k = L_{k-1}^{-1} a_k \quad(\text{forward substitution}),
\qquad
\lambda_k = \sqrt{\alpha_k - l_k^\ast l_k}.
$$

Positive definiteness of the principal submatrix $A_k$ guarantees
$\alpha_k - l_k^\ast l_k > 0$, so the square root is real and the diagonal stays
positive. The name reads as a matrix square root: Cholesky is to s.p.d. matrices
what $\sqrt{\alpha}$ is to positive numbers.[^bor-chol]

$$
% caption: Cholesky reads only the lower half of the symmetric $A$ and produces
% one triangular factor $L$; row $k$ needs a forward solve against $L_{k-1}$ plus a
% square root for the diagonal entry $\lambda_k$.
\begin{tikzpicture}[scale=1.0, font=\footnotesize]
\definecolor{acc}{HTML}{4A6FA5}
% A symmetric, lower half shaded
\draw[black] (0,0) rectangle (2.4,2.4);
\draw[black] (0,2.4)--(2.4,0);
\fill[acc!10] (0,2.4)--(2.4,0)--(0,0)--cycle;
\node[black] at (1.5,1.9) {$A^\ast{=}A$};
\node[acc] at (0.7,0.9) {read};
\node at (1.2,-0.4) {$A$};
\draw[->, black, thick] (2.9,1.2) -- (3.9,1.2);
% L lower triangle with row k highlighted
\draw[black] (4.4,0) rectangle (6.8,2.4);
\draw[black] (4.4,2.4)--(6.8,0);
\fill[acc!10] (4.4,2.4)--(6.8,0)--(4.4,0)--cycle;
\draw[acc, thick] (4.4,0.8) rectangle (6.8,1.2);
\node[acc, anchor=west] at (7.0,1.0) {row $k$};
\node at (5.6,-0.4) {$L$};
\node[acc, align=center] at (5.6,-1.15) {$A = L L^\ast$};
\end{tikzpicture}
$$

```algorithm
caption: $\textsc{Cholesky}(A)$ — factor s.p.d. $A = L L^\ast$ row by row
for $k = 1$ to $m$ do
  $l \gets L_{1:k-1,\,1:k-1}^{-1}\, A_{1:k-1,\,k}$ // forward substitution
  $L_{k,\,1:k-1} \gets l^\ast$
  $L_{kk} \gets \sqrt{A_{kk} - l^\ast l}$ // fails (root of $\le 0$) if $A$ is not s.p.d.
end for
return $L$
```

> **Worked example.** Cholesky-factor the s.p.d. matrix
>
> $$
> A = \begin{bmatrix} 4 & 2 & 2 \\ 2 & 5 & 3 \\ 2 & 3 & 6 \end{bmatrix}.
> $$
>
> Column 1: $\lambda_{11} = \sqrt{4} = 2$, then $l_{21} = 2/2 = 1$ and $l_{31} =
> 2/2 = 1$. Column 2: $\lambda_{22} = \sqrt{5 - 1^2} = 2$, then $l_{32} = (3 -
> l_{31}l_{21})/\lambda_{22} = (3 - 1)/2 = 1$. Column 3: $\lambda_{33} = \sqrt{6 -
> l_{31}^2 - l_{32}^2} = \sqrt{6 - 1 - 1} = 2$. Hence
>
> $$
> L = \begin{bmatrix} 2 & 0 & 0 \\ 1 & 2 & 0 \\ 1 & 1 & 2 \end{bmatrix},
> \qquad
> L L^\top = A.
> $$
>
> Every radicand stayed positive, which certifies $A$ positive definite; had one
> been $\le 0$, the square root would have failed and the algorithm would have
> reported the matrix as indefinite.

Only the lower half of $A$ is touched, and only one triangular factor is stored, so
Cholesky costs

$$
\#\text{flop} \approx \sum_{k=1}^{m} k^2 \approx \frac{1}{3}m^3,
$$

half the cost of a general $LU$, with the same efficiency ratio $q \approx m/3$ and
the same near-peak Level-3 BLAS performance. The square root also doubles as a
test: if the argument $\alpha_k - l_k^\ast l_k$ is ever nonpositive, $A$ was not
positive definite, and the algorithm reports it.

## Choosing a solver

| Situation | Method | Cost (flop) | Pivoting |
| --- | --- | --- | --- |
| General invertible $A$ | pivoted $LU$, $P^\ast A = LU$ | $\tfrac{2}{3}m^3$ | partial (rows) |
| Symmetric positive-definite $A$ | Cholesky, $A = LL^\ast$ | $\tfrac{1}{3}m^3$ | none needed |
| Many right-hand sides, fixed $A$ | factor once, substitute per $b$ | $2m^2$ per solve | reuse the factors |

The recurring pattern is to factor once and then answer each new right-hand side
with two triangular solves. How much error a solve inherits from the matrix is set
by its [conditioning](/linear-algebra/numerical-linear-algebra/conditioning-and-floating-point);
whether the algorithm reaches that floor is the question of
[stability](/linear-algebra/numerical-linear-algebra/stability-and-error-analysis),
where pivoted $LU$ turns out to be, in the backward sense, as accurate as the
problem allows.

[^bor-lu]: **Bornemann**, _Numerical Linear Algebra_, §7 — Triangular Decomposition: the $A = LU$ recursion via Schur complements, uniqueness, the $\tfrac{2}{3}m^3$ cost, the Zeroth Law, and partial pivoting giving $P^\ast A = LU$ with $\lvert L \rvert \le 1$.
[^bor-chol]: **Bornemann**, §8 — Cholesky Decomposition: the s.p.d. definition, the row-wise construction $l_k = L_{k-1}^{-1}a_k$, $\lambda_k = \sqrt{\alpha_k - l_k^\ast l_k}$, and the $\tfrac{1}{3}m^3$ cost.