38. LU, Cholesky, QR decompositions (unified review)

This course has now built four matrix factorizations across three modules: LU (§16.5), QR (§20.6), the Spectral/SVD pair (§21.0, §21.4). This short opening lesson of Module 22 puts the first three side by side, adds one new one — Cholesky, LU's specialization for positive-definite matrices — and asks the practical question each was motivated by but never stated outright: which decomposition for which job?

The four, side by side

Each row removes a restriction the row above it needed: LU wants a square matrix; QR drops that entirely; SVD drops even the independent-columns requirement. This isn't a coincidence — it's the order that matches how much structure each proof needs to lean on.

Cholesky decomposition

For A symmetric positive definite (§21.2), Cholesky is LU specialized to exploit the symmetry: A=LL^T for a unique lower-triangular L with positive diagonal entries. It's derivable directly from LU: A=LU (§16.5) combined with A=A^T=U^TL^T and uniqueness of LU forces U=DL^T for a diagonal D of positive entries (PD guarantees this), and absorbing \sqrt D into L gives A=(L\sqrt D)(L\sqrt D)^T, exactly Cholesky's form.

Why it's worth having as its own algorithm: it costs half the arithmetic of general LU (only one triangular factor to compute, not two), and it doubles as a positive-definiteness test — Cholesky succeeds if and only if A is PD (an attempt to take \sqrt{} of a negative number during the algorithm signals failure), making it the practical alternative to §21.2's eigenvalue or Sylvester tests.

Choosing a decomposition

• Solving A\vec x=\vec b for square, well-behaved A: LU — the cheapest option, and §16.5's reusable-across-many-\vec b property still applies.
• A symmetric positive definite (common for covariance matrices, and normal-equation matrices A^TA from §20.5): Cholesky — half the cost of LU, plus a free PD check.
• Least squares, or anything needing an orthonormal basis for \operatorname{Col}(A): QR — §20.6's numerically stable route, avoiding A^TA's error amplification entirely.
• Rank, best low-rank approximation, or genuinely singular/rectangular A: SVD — the most expensive to compute, but the most information it reveals, and the only one of the four that never fails to exist.

Doing it in Python

import numpy as np
import scipy.linalg as la

A = np.array([[4., 2., 2.], [2., 5., 3.], [2., 3., 6.]])   # symmetric PD

L_chol = la.cholesky(A, lower=True)
print("Cholesky L =")
for row in L_chol.tolist():
    print([round(v, 4) for v in row])
print("L L^T == A:", bool(np.allclose(L_chol @ L_chol.T, A)))

P, L_lu, U_lu = la.lu(A)
n_lu = np.count_nonzero(L_lu) + np.count_nonzero(U_lu)
n_chol = np.count_nonzero(L_chol)
print(f"\nLU stores {n_lu} nonzeros total (L+U); Cholesky stores {n_chol} (L alone)")

Cholesky L =
[2.0, 0.0, 0.0]
[1.0, 2.0, 0.0]
[1.0, 1.0, 2.0]
L L^T == A: True

LU stores 12 nonzeros total (L+U); Cholesky stores 6 (L alone)

Cholesky as a PD test — succeeding for a PD matrix, failing for an indefinite one:

import numpy as np
import scipy.linalg as la

A_pd = np.array([[2., 1.], [1., 2.]])
A_indef = np.array([[1., 2.], [2., 1.]])   # from Section 21.2/21.3

for name, A in [("A_pd", A_pd), ("A_indef", A_indef)]:
    try:
        la.cholesky(A)
        print(f"{name}: Cholesky succeeded -> positive definite")
    except np.linalg.LinAlgError:
        print(f"{name}: Cholesky failed -> NOT positive definite")

A_pd: Cholesky succeeded -> positive definite
A_indef: Cholesky failed -> NOT positive definite

Worked example

Find the Cholesky decomposition of A=\begin{pmatrix}9&3\\3&5\end{pmatrix} by hand.

L=\begin{pmatrix}l_{11}&0\\l_{21}&l_{22}\end{pmatrix}. From LL^T=A: l_{11}^2=9\Rightarrow l_{11}=3 (positive root, by convention). l_{11}l_{21}=3\Rightarrow l_{21}=1. l_{21}^2+l_{22}^2=5\Rightarrow1+l_{22}^2=5\Rightarrow l_{22}=2.

\boxed{L=\begin{pmatrix}3&0\\1&2\end{pmatrix}}

Sanity check. LL^T=\begin{pmatrix}3&0\\1&2\end{pmatrix}\begin{pmatrix}3&1\\0&2\end{pmatrix}=\begin{pmatrix}9&3\\3&5\end{pmatrix}=A ✓. Every entry under the square root was positive throughout (9 and then 4), confirming A is genuinely PD — matching §21.2's Sylvester criterion independently: leading minors 9>0 and \det A=45-9=36>0.

Your turn

1. Why does Cholesky fail immediately (rather than producing a wrong answer) when applied to a non-PD matrix?

2. For a 1000\times1000 system that needs solving for 50 different right-hand sides, which decomposition minimizes total work, and why?

3. True or false: Cholesky decomposition is unique.

Solutions

1. The algorithm computes l_{ii}=\sqrt{a_{ii}-\sum_{k<i}l_{ik}^2} at each diagonal step — if A isn't PD, this quantity can come out negative, requiring a square root of a negative number (undefined for real L). The algorithm doesn't silently produce a nonsensical answer; it halts exactly where positive-definiteness is violated, which is precisely why it doubles as a PD test rather than needing a separate eigenvalue check.

2. LU (or Cholesky, if additionally symmetric PD): factor once (O(n^3)), then each of the 50 solves is a cheap O(n^2) triangular pair (§16.5) — total cost dominated by the single factorization, not by the 50 right-hand sides. Recomputing a fresh elimination for each \vec b separately (no factorization reuse) would cost 50\times O(n^3) instead — the entire point of §16.5's reusability argument, now with concrete numbers attached.

3. True, with the stated convention (positive diagonal). If L_1L_1^T=L_2L_2^T=A with both diagonals positive, then L_1=L_2 — this is what makes Cholesky useful as a canonical form (unlike LU with pivoting, whose P can vary) and is part of why it's preferred whenever applicable.

Check yourself in code

Find the Cholesky factor of A=\begin{pmatrix}16&4\\4&5\end{pmatrix} and verify LL^T=A.

Print exactly this:

L =
[4.0, 0.0]
[1.0, 2.0]
L L^T == A: True

import numpy as np
import scipy.linalg as la

A = np.array([[16., 4.], [4., 5.]])
L = la.cholesky(A, lower=True)
print("L =")
for row in L.tolist():
    print(row)
# print whether L @ L.T equals A

import numpy as np
import scipy.linalg as la

A = np.array([[16., 4.], [4., 5.]])
L = la.cholesky(A, lower=True)
print("L =")
for row in L.tolist():
    print(row)
print("L L^T == A:", bool(np.allclose(L @ L.T, A)))

LU, Cholesky, QR, and SVD each trade computational cost against generality: LU is cheapest for square systems, Cholesky halves that cost when symmetry and positive-definiteness apply, QR handles any shape stably for least squares, and SVD — the most expensive — never fails to exist and reveals the most structure. Reusing a factorization across many right-hand sides is the standard reason to compute one at all.

Next: a decomposition that works for any square matrix, even defective ones §19.4 could only handle with Jordan form — the Schur decomposition.