QR and Cholesky: Factorizations That Solve

By Better Wise

Module 1 turned KL divergences into objectives. Every one of those objectives, once you differentiate it, hands you a linear system: a Gauss-Newton step, a ridge-regularized readout, a Gaussian process posterior. This lesson is about the two factorizations you will reach for most often, and about the one-line shortcut that silently destroys your answer while looking like the textbook derivation.

By the end you should be able to:

• Write the least-squares problem and its normal equations, and state exactly how much accuracy the shortcut throws away.
• Factor a matrix with Cholesky and with Householder QR, and quote the flop cost of each.
• Predict when Cholesky fails, and choose between jitter, a ridge term, and a different factorization.
• Read a condition number and pick the route before you run the code.

The normal equations and their exact price

For $\min_x \lVert Ax - b \rVert_2$ with $A \in \mathbb{R}^{m \times n}$, $m \ge n$, setting the gradient to zero gives the shortest derivation in numerical linear algebra:

$$ A^\top A x = A^\top b $$

The matrix $A^\top A$ is symmetric and positive semidefinite, so Cholesky applies. The catch is the condition number:

$$ \kappa(A^\top A) = \kappa(A)^2 $$

With double-precision unit roundoff $\varepsilon \approx 1.1\times10^{-16}$, a backward-stable QR solve carries relative error near $\varepsilon,\kappa(A)$, while the normal-equations route carries roughly $\varepsilon,\kappa(A)^2$. Squaring a condition number means squaring your lost digits.

Worked example (why this is not academic). Take

$$ A = \begin{bmatrix} 1 & 1 \ 1 & 1+\delta \end{bmatrix}, \qquad \delta = 10^{-8} $$

Exact arithmetic gives $\det(A^\top A) = \delta^2 = 10^{-16}$, while the entries of $A^\top A$ are all near $4$. The computed determinant is a difference of two numbers of size $4$, whose rounding error is about $4\varepsilon \approx 4.4\times10^{-16}$, which exceeds the true value. The computed determinant is noise and can come out negative, at which point Cholesky dies on a negative pivot. Meanwhile $\kappa(A) \approx 4/\delta = 4\times10^{8}$: QR keeps about seven digits, the normal equations keep none. Same matrix, same machine, two different answers.

Cholesky: fast when the matrix is honestly SPD

Cholesky computes $A = LL^\top$ with $L$ lower triangular and positive diagonal. It exists if and only if $A$ is symmetric positive definite, costs about $n^3/3$ flops (half of LU), and each subsequent solve costs $2n^2$ for two triangular substitutions. For SPD input it is backward stable: the computed factor satisfies $LL^\top = A + E$ with $\lVert E\rVert \approx \varepsilon \lVert A \rVert$.

Two readings worth keeping. First, $L$ is a matrix square root of $A$: to draw $z \sim \mathcal{N}(0, \Sigma)$, factor $\Sigma = LL^\top$ and set $z = L\xi$ with $\xi \sim \mathcal{N}(0, I)$. Second, $L$ can be updated in $O(n^2)$ when a rank-one term is added to $A$, which is why recursive least squares and Kalman filters keep a running Cholesky factor rather than refactoring. Downdates are the fragile direction: removing a rank-one term can push a pivot through zero.

The hard limit: symmetry and positive definiteness are properties $A^\top A$ has in exact arithmetic only. Form it in floating point and the property can be false. A negative pivot mid-factorization is your signal, not a bug in LAPACK.

Warning: Never "fix" a negative pivot by taking the absolute value. A negative pivot means the matrix you passed in is not positive definite, so any factor you force out describes a different problem than the one you meant to solve.

Householder QR: paying flops for digits

QR computes $A = QR$ with $Q$ having orthonormal columns and $R$ upper triangular. Householder reflectors build it in about $2mn^2 - \tfrac{2}{3}n^3$ flops, roughly twice the Cholesky-on-normal-equations cost when $m \gg n$, and then $Rx = Q^\top b$ costs $n^2$. It never forms $A^\top A$, so it never squares the condition number.

1. Factor A = Q R with Householder reflectors (never form Q explicitly; keep the reflectors).
2. Apply Q^T to b: c = Q^T b, in place.
3. Solve R x = c[0:n] by back substitution.
4. Check the residual ||A x - b|| and the diagonal of R.

A note on Gram-Schmidt, because it is the version people write by hand: classical Gram-Schmidt effectively squares the condition number too (error near $\varepsilon\kappa^2$), and modified Gram-Schmidt improves this to about $\varepsilon\kappa$ at the cost of slowly losing orthogonality in $Q$. Householder QR is the default for a reason.

If the diagonal of $R$ decays toward zero, the problem is rank deficient or nearly so. QR with column pivoting reorders columns to push the small diagonal entries to the bottom right, giving a usable numerical rank. It is a diagnostic, not a cure; the minimum-norm solution and a defensible truncation threshold come from the SVD, which is the next lesson.

Cost and accuracy side by side

flowchart TD
    Start["Linear system or least squares"] --> SPD{"A square and numerically SPD?"}
    SPD -- yes --> Chol["Cholesky A = L L^T<br/>n^3/3 flops"]
    SPD -- no --> QR["Householder QR A = Q R<br/>error ~ eps * kappa(A)"]
    QR --> Piv{"Diagonal of R collapsing?"}
    Piv -- yes --> PQR["Column pivoting<br/>numerical rank"]
    Piv -- no --> Tri["Triangular solve"]
    Chol --> Tri
    PQR --> Tri
    Tri --> Out["x"]

Worked example (ridge makes Cholesky honest again). Suppose $\sigma_{\max} = 1$ and $\sigma_{\min} = 10^{-6}$, so $\kappa(A) = 10^{6}$ and $\kappa(A^\top A) = 10^{12}$. Normal equations give relative error near $10^{-4}$: four digits gone before you look at the data. Add a ridge term $\lambda = 10^{-6}$ and the effective condition number becomes

$$ \kappa_\lambda = \frac{\sigma_{\max}^2 + \lambda}{\sigma_{\min}^2 + \lambda} = \frac{1 + 10^{-6}}{10^{-12} + 10^{-6}} \approx 10^{6} $$

which restores about six more digits and guarantees a positive pivot. The regularizer fixes the statistics and the arithmetic at the same time, a coincidence worth remembering when the next lesson asks you to justify $\lambda$.

Check yourself. You have $m = 10^5$, $n = 10^3$, and $\kappa(A) = 10^3$ in double precision. Which route, and why?

Click to reveal the answer

Normal equations plus Cholesky is defensible here: $\varepsilon\kappa^2 \approx 10^{-10}$, far below any tolerance you care about, and it costs about half the flops of QR ($mn^2 = 10^{11}$ versus $2mn^2$). If that condition number were $10^6$ instead, the same reasoning flips: $\varepsilon\kappa^2 \approx 10^{-4}$ and you switch to Householder QR.

Choosing a solver without regret

• Threshold rule: the normal-equations route is safe while $\kappa(A) \lesssim \sqrt{1/\varepsilon}$, about $10^8$ in double and about $4\times10^3$ in float32. If you cannot bound $\kappa$, assume the worst and use QR.
• Give Cholesky a reason to succeed. For kernel and Gram matrices, add jitter proportional to the mean diagonal, and check the smallest pivot rather than catching an exception and retrying blindly.
• Never invert $L$, $R$, or $A$. Triangular solves cost $n^2$; explicit inverses cost $n^3$ and lose accuracy.
• Compute the residual $b - Ax$ with care (or in higher precision) when you want to verify a solve; forming it naively in the working precision loses exactly the digits you are trying to measure.
• Common mistake: solve(A.T @ A, A.T @ b) in production code. It is one character shorter than the stable version and one condition number worse.

What to carry forward

Cholesky is the fast, honest route for symmetric positive definite systems, and it stops being honest the moment positive definiteness is only approximate. Householder QR costs about twice as much and never squares the condition number, which is why it is the default for least squares. Column pivoting tells you when the rank is collapsing but does not tell you what to do about it.

That last gap is exactly where the next lesson starts: the full SVD gives the pseudoinverse, a principled truncation threshold, and the rank control that pivoted QR can only gesture at.

Forming the normal equations replaces the condition number κ(A) with _____, which is why their accuracy degrades quadratically. Cholesky factors a matrix as A = LLᵀ and exists if and only if A is _____. Householder QR costs about _____ flops for an m-by-n matrix with m ≥ n. A mid-factorization _____ in Cholesky means the input matrix is not numerically positive definite. Adding a ridge term λI improves the effective condition number to _____.

Front: What are the normal equations for least squares? Back: AᵀA x = Aᵀb, obtained by setting the gradient of ||Ax − b||² to zero.

Front: How does forming AᵀA change the condition number? Back: It squares it: κ(AᵀA) = κ(A)².

Front: Relative error in x for a backward-stable QR least-squares solve? Back: About ε·κ(A) in the working precision.

Front: Relative error in x when solving via the normal equations? Back: About ε·κ(A)², which is quadratically worse.

Front: What factorization does Cholesky compute, and for which matrices? Back: A = LLᵀ with positive diagonal L, for symmetric positive definite A only.

Front: Flop cost of Cholesky factorization versus LU? Back: About n³/3 versus 2n³/3, so Cholesky is roughly twice as cheap.

Front: Cost of Householder QR for m ≥ n? Back: About 2mn² − (2/3)n³ flops, plus n² for the triangular solve.

Front: What does a negative pivot during Cholesky tell you? Back: The input is not numerically positive definite, so the problem or the matrix formation is wrong.

Front: Why is classical Gram-Schmidt avoided in production? Back: Its error grows like ε·κ(A)², so it squares the condition number like the normal equations.

Front: What does QR with column pivoting buy you? Back: A usable numerical rank estimate from the decaying diagonal of R, at roughly twice the cost of plain QR.

Front: How does adding λI change the condition number of AᵀA? Back: κ becomes (σ_max² + λ)/(σ_min² + λ), which shrinks toward 1 as λ grows.

Front: When is the normal-equations route safe in double precision? Back: While κ(A) is below roughly √(1/ε) ≈ 10⁸, and far below that if you want margin.

Front: Why keep a Cholesky factor in a Kalman filter instead of refactoring? Back: Rank-one updates refresh L in O(n²) instead of O(n³).

Front: How do you sample from N(0, Σ) using Cholesky? Back: Factor Σ = LLᵀ, draw ξ ~ N(0, I), and return z = Lξ.