LU and Cholesky Factorization in Practice

Solving $A x = 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.¹ The decomposition is $A = LU$ , and once it is in hand the system reduces to two triangular solves.

Reading elimination as a factorization

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

L z = b (forward substitution), U x = z (back substitution),

so that $L (U x) = A x = 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$ .¹

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} = [α_{k} b_{k} u_{k}^{*} B_{k}] = [1 l_{k} L_{k + 1}] [α_{k} u_{k}^{*} U_{k + 1}] .

Multiplying out the second block row gives $b_{k} = l_{k} α_{k}$ and $B_{k} = l_{k} u_{k}^{*} + L_{k + 1} U_{k + 1}$ , so provided the pivot $α_{k} \neq = 0$ ,

l_{k} = \frac{b _{k}}{α _{k}}, A_{k + 1} = B_{k} - l_{k} u_{k}^{*} .

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

One step of

LU

: the pivot

α_{k}

and first row/column are read off, then the rank-one update

- l_{k} u_{k}^{*}

overwrites the trailing Schur complement, shrinking the active block by one.

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

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

and the two substitutions that follow add only $2 m^{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 = [0111]

has $α_{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 = [1 0^{- 20} 1 11] = [1 1 0^{20} 01] [1 0^{- 20} 0 1 1 - 1 0^{20}] .

On a machine carrying sixteen significant digits, $1 - 1 0^{20}$ rounds to $- 1 0^{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 $A x = e_{1}$ with the rounded factors returns $\overset{x}{^} = (0, 1)^{⊤}$ against the true $x \approx (- 1, 1)^{⊤}$ : a one-hundred-percent error in the first component from a matrix that is perfectly well-behaved.

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.

Each row swap is a transposition $τ_{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.

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 $∣ l_{k} ∣ \leq 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.

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

\leq 1

in size.

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

Algorithm:

\textsc{PivotedLU}(A)

— in-place

P^\ast A = LU

with partial pivoting

1
$p \gets [1, 2, \dots, m]$
permutation record
2
for $k = 1$ to $m-1$ do
3
choose $j \ge k$ maximizing $\lvert A_{jk} \rvert$
pivot search
4
swap rows $k$ and $j$ of $A$ and entries $k, j$ of $p$
5
for $i = k+1$ to $m$ do
6
$A_{ik} \gets A_{ik} / A_{kk}$
multiplier $l$ , stored below the diagonal
7
end for
8
$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
9
end for
10
return $L = \text{unit lower part},\; U = \text{upper part},\; p$

Worked example. Solve $A x = b$ in three-significant-digit arithmetic with

A = [1 0^{- 4} 1 11], b = [12], x_{true} = [1.0001 0.9999] .

Without pivoting, the multiplier is $l = 1/1 0^{- 4} = 1 0^{4}$ . Eliminating the second row,

a_{22} = 1 - 1 0^{4} \cdot 1 = - 9999 \to - 1.00 \times 1 0^{4}, b_{2} = 2 - 1 0^{4} \cdot 1 = - 9998 \to - 1.00 \times 1 0^{4},

where rounding to three digits discards the trailing $1$ and $2$ . Back substitution gives $ξ_{2} = (- 1.00 \times 1 0^{4}) / (- 1.00 \times 1 0^{4}) = 1.00$ and then $ξ_{1} = (1 - 1.00) /1 0^{- 4} = 0$ . The computed $\overset{x}{^} = (0, 1.00)^{⊤}$ is wrong in its first component by $100%$ , on a matrix that is well-behaved.

With partial pivoting, $∣ 1 ∣ > ∣ 1 0^{- 4} ∣$ forces a row swap. Now $l = 1 0^{- 4} /1 = 1 0^{- 4}$ , and

a_{22} = 1 - 1 0^{- 4} \cdot 1 = 0.9999 \to 1.00, b_{2} = 1 - 1 0^{- 4} \cdot 2 = 0.9998 \to 1.00,

so $ξ_{2} = 1.00$ and $ξ_{1} = (2 - 1.00) /1 = 1.00$ . The computed $\overset{x}{^} = (1.00, 1.00)^{⊤}$ matches the true solution to all three digits carried.

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.

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.

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.

The factor is built row by row. Partitioning $A_{k} = L_{k} L_{k}^{*}$ and matching blocks gives, for row $k$ ,

l_{k} = L_{k - 1}^{- 1} a_{k} (forward substitution), λ_{k} = α_{k} - l_{k}^{*} l_{k} .

Positive definiteness of the principal submatrix $A_{k}$ guarantees $α_{k} - l_{k}^{*} 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 $α$ is to positive numbers.²

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

λ_{k}

Algorithm:

\textsc{Cholesky}(A)

— factor s.p.d.

A = L L^\ast

row by row

1
for $k = 1$ to $m$ do
2
$l \gets L_{1:k-1,\,1:k-1}^{-1}\, A_{1:k-1,\,k}$
forward substitution
3
$L_{k,\,1:k-1} \gets l^\ast$
4
$L_{kk} \gets \sqrt{A_{kk} - l^\ast l}$
fails (root of $\le 0$ ) if $A$ is not s.p.d.
5
end for
6
return $L$

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

# flop \approx k = 1 \sum 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 $α_{k} - l_{k}^{*} 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^{*} A = LU$	$\frac{2}{3} m^{3}$	partial (rows)
Symmetric positive-definite $A$	Cholesky, $A = L L^{*}$	$\frac{1}{3} m^{3}$	none needed
Many right-hand sides, fixed $A$	factor once, substitute per $b$	$2 m^{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; whether the algorithm reaches that floor is the question of stability, where pivoted $LU$ turns out to be, in the backward sense, as accurate as the problem allows.

Bornemann, Numerical Linear Algebra, §7 — Triangular Decomposition: the $A = LU$ recursion via Schur complements, uniqueness, the $\frac{2}{3} m^{3}$ cost, the Zeroth Law, and partial pivoting giving $P^{*} A = LU$ with $∣ L ∣ \leq 1$ . ↩ ↩² ↩³
Bornemann, §8 — Cholesky Decomposition: the s.p.d. definition, the row-wise construction $l_{k} = L_{k - 1}^{- 1} a_{k}$ , $λ_{k} = α_{k} - l_{k}^{*} l_{k}$ , and the $\frac{1}{3} m^{3}$ cost. ↩

Reading elimination as a factorization

Small pivots and instability

Partial pivoting

Cholesky for symmetric positive-definite systems

Choosing a solver

Footnotes