Matrix Decompositions: LU, QR, Cholesky

Splitting one matrix into two triangular pieces turns slow elimination into fast substitution, and the three standard splits cover general, rectangular and symmetric matrices.

Elimination gives A equals L times U, orthonormal columns give A equals Q times R, and symmetric positive definite A gives L times L transpose; triangular diagonals hold the pivots

Why Does This Exist?

Fitting a regression, updating a Kalman filter, taking a Newton step and evaluating a Gaussian likelihood all end in the same place: solve $Ax = b$ for $x$ , given a square matrix $A$ and a vector $b$ . Inverting $A$ directly costs about $2n^3$ flops and squares the condition number, so the computed inverse is less accurate than the factors it came from. Nobody inverts.

Factoring once and substituting twice is the replacement. Elimination on $A$ costs about $2n^3/3$ flops, and each later solve against a new $b$ costs only about $n^2$ . Ten right-hand sides therefore cost little more than one. This page covers the three factorizations and when each applies; linear systems and least squares covers how the system is set up, and determinants covers what the pivots say about it in determinants, inverses and trace.

The one idea you need first

A triangular system needs no elimination. An upper triangular $Ux = y$ solves bottom-up: $x_n = y_n / U_{nn}$ , then each row above uses the values already found. A lower triangular $Ly = b$ solves top-down the same way. Each pass costs about $n^2/2$ flops. The whole point of every decomposition below is to convert the hard problem $Ax = b$ into two easy triangular ones.

Think of It Like This

Disassembling a wardrobe into labelled stacks

You need to move a built wardrobe through a narrow door, so you take it apart. Two stacks result. One stack holds every panel in the order it came off, each labelled with which screws held it and where: that is $L$ , the record of the elimination steps as multipliers. The other stack is the bare frame with a staircase profile, widest at the bottom step: that is $U$ , the staircase shape elimination leaves behind. The screws themselves sit on the diagonal of each stack, one per panel: those are the pivots. If a screw strips (a zero pivot), you reorder the panels before continuing: that is row pivoting. Reassembling in reverse order solves the system.

Where it stops: wardrobe panels are exact while matrix factors round every arithmetic step, and rounding piles up exactly where the condition number says it will. The analogy also hides the rectangular case: QR handles wardrobes with more panels than frame slots, which has no furniture equivalent.

How It Actually Works

LU records elimination. Eliminating $x_1$ from rows $2$ through $n$ uses multipliers $l_{i1} = A_{i1} / A_{11}$ , stored below the diagonal of $L$ with ones on its diagonal. Repeating down the columns gives $A = LU$ : $L$ unit lower triangular holding every multiplier, $U$ upper triangular holding the staircase. Solving $Ax = b$ becomes $Ly = b$ (forward substitution) then $Ux = y$ (backward substitution). The determinant falls out free: $\det(A)$ is the product of the diagonal of $U$ , up to the sign of the row swaps.

Pivoting fixes zero and tiny pivots. If $A_{11}$ is zero, the multiplier divides by zero and elimination halts even when $A$ is perfectly invertible. Partial pivoting swaps in the largest entry below before each column, recorded in a permutation $P$ with $PA = LU$ . The symptom that you skipped it is a division by zero on a matrix you know is nonsingular, or garbage answers from a tiny pivot amplifying rounding. Always pivot; unpivoted LU is a teaching tool, not a solver.

QR handles rectangles and least squares

When $A$ is $m \times n$ with $m \geq n$ , there is no inverse and LU does not apply. QR writes $A = QR$ with $Q$ holding orthonormal columns ( $Q^T Q = I$ ) and $R$ upper triangular. Least squares $\min_x \lVert Ax - b \rVert_2$ becomes $Rx = Q^T b$ , one triangular solve, with no squaring of the condition number the way forming $A^T A$ would. Householder reflections build it in about $2mn^2 - 2n^3/3$ flops. Cost rule of thumb: QR costs roughly twice LU on a square matrix, and buys stability for exactly the problems LU cannot touch.

Cholesky halves the work for symmetric positive definite matrices

If $A$ is symmetric ( $A^T = A$ ) with all eigenvalues positive, it factors as $A = LL^T$ : one lower triangular $L$ with positive diagonal, times its own transpose. Derivation per entry: $L_{jj} = \sqrt{A_{jj} - \sum_{k<j} L_{jk}^2}$ , and each row below uses only values already computed. Cost is about $n^3/3$ flops, half of LU, with half the storage. Covariance matrices and the $X^T X + \lambda I$ of ridge regression are the textbook cases. If $A$ is not positive definite, some radicand goes negative and the algorithm stops: that failure is diagnostic, not a bug to patch over.

Worked example

Take $A = \begin{bmatrix} 4 & 3 \\ 6 & 3 \end{bmatrix}$ with $b = \begin{bmatrix} 7 \\ 9 \end{bmatrix}$ .

Eliminate with multiplier $l_{21} = 6 / 4 = 1.5$ . Row 2 becomes $[6 - 1.5 \times 4,\ 3 - 1.5 \times 3] = [0,\ -1.5]$ . So:

L = \begin{bmatrix} 1 & 0 \\ 1.5 & 1 \end{bmatrix}, \qquad U = \begin{bmatrix} 4 & 3 \\ 0 & -1.5 \end{bmatrix}, \qquad \det(A) = 4 \times (-1.5) = -6

Direct check: $12 - 18 = -6$ . Forward substitution $Ly = b$ : $y_1 = 7$ , then $1.5 \times 7 + y_2 = 9$ gives $y_2 = 9 - 10.5 = -1.5$ . Backward substitution $Ux = y$ : $-1.5\, x_2 = -1.5$ gives $x_2 = 1$ , then $4x_1 + 3 = 7$ gives $x_1 = 1$ . Verify: $4 + 3 = 7$ and $6 + 3 = 9$ .

The companions on small numbers: $B = \begin{bmatrix} 1 & 1 \\ -1 & 1 \end{bmatrix}$ has orthonormal columns up to $1/\sqrt{2} \approx 0.7071$ , so $Q = \begin{bmatrix} 0.7071 & 0.7071 \\ -0.7071 & 0.7071 \end{bmatrix}$ and $R = \begin{bmatrix} 1.4142 & 0 \\ 0 & 1.4142 \end{bmatrix}$ with $\sqrt{2} \approx 1.4142$ . And $C = \begin{bmatrix} 4 & 2 \\ 2 & 3 \end{bmatrix}$ gives $L_{11} = \sqrt{4} = 2$ , $L_{21} = 2/2 = 1$ , $L_{22} = \sqrt{3 - 1} = \sqrt{2} \approx 1.4142$ . Check the product: $2 \times 1 = 2$ and $1 + 2 = 3$ .

Code

The point is the check: multiply the factors back and confirm the solve leaves no residual. Every # -> value below is the printed result.

1import numpy as np2
3A = np.array([[4.0, 3.0], [6.0, 3.0]])4L = np.array([[1.0, 0.0], [1.5, 1.0]])5U = np.array([[4.0, 3.0], [0.0, -1.5]])6b = np.array([7.0, 9.0])7
8print(np.allclose(L @ U, A))            # -> True9print(round(np.linalg.det(A), 6))       # -> -6.010
11y = np.linalg.solve(L, b)               # forward leg, L has unit diagonal12x = np.linalg.solve(U, y)               # backward leg13print(np.round(x, 6))                   # -> [1. 1.]14print(np.round(A @ x - b, 10))          # -> [0. 0.]15
16C = np.array([[4.0, 2.0], [2.0, 3.0]])17Lc = np.linalg.cholesky(C)18print(np.round(Lc, 4))                  # -> [[2.     0.    ] [1.     1.4142]]19print(np.allclose(Lc @ Lc.T, C))        # -> True

Replacing the hand-derived $L$ and $U$ with scipy.linalg.lu returns the same factors up to the row-swap convention; the residual stays zero either way.

Watch Out For

Running LU with no pivoting on a matrix with a zero leading entry

The matrix above with its rows swapped, $\begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix}$ , is perfectly invertible, yet unpivoted elimination divides by $A_{11} = 0$ on the first step and halts. The more insidious version is a tiny but nonzero pivot: no crash, just answers wrong by orders of magnitude because one huge multiplier amplifies every rounding error below it.

The fix is partial pivoting ( $PA = LU$ ), which every serious routine does by default: scipy.linalg.lu returns the permutation, and numpy.linalg.solve pivots internally. If your hand-rolled elimination hits a zero pivot, swap in the largest entry below before dividing, and record the swap. Never special-case the zero away.

Calling Cholesky on a matrix that is not positive definite

Pass a merely symmetric matrix such as $\begin{bmatrix} 1 & 2 \\ 2 & 1 \end{bmatrix}$ (eigenvalues $3$ and $-1$ ) and the algorithm reaches $L_{22} = \sqrt{1 - 4} = \sqrt{-3}$ and raises LinAlgError. The usual cause is forming a covariance from uncentered or rank-deficient data, or rounding that pushed a barely positive eigenvalue negative.

Check before factoring: symmetry up to tolerance, then either an eigenvalue test or simply attempt the factorization and treat failure as the diagnosis. If it fails, the matrix was never Cholesky-eligible: fall back to LU with pivoting for square systems or QR for least squares, and fix the modelling (add the ridge diagonal $\lambda I$ , center the data) rather than clipping the negative radicand to zero.

The Quick Version

Factor once at $\sim 2n^3/3$ flops, then solve each new right-hand side at $\sim n^2$ : that is the whole economy.
$A = LU$ stores elimination multipliers in unit $L$ and the staircase in $U$ ; pivoting ( $PA = LU$ ) is mandatory in practice.
The worked $2 \times 2$ gives $L$ with $1.5$ , $U$ with $-1.5$ , determinant $-6$ , and solution $x = [1, 1]$ against $b = [7, 9]$ .
Determinant is the product of $U$ 's diagonal up to row-swap signs; triangular solves run forward then backward.
$A = QR$ covers rectangular $A$ and least squares without squaring the condition number, at about twice the cost of LU.
$A = LL^T$ halves cost and storage for symmetric positive definite $A$ ; the worked $L$ holds $2$ , $1$ and $\sqrt{2} \approx 1.4142$ .
Never invert explicitly, and never run unpivoted LU or Cholesky outside its eligibility.