QR Factorization

Gram–Schmidt is a beautiful hand process — take a pile of vectors, strip away the overlaps, and out come tidy, mutually perpendicular unit vectors. But a computer doesn't want a "process" narrated step by step; it wants a matrix equation it can store, reuse, and hand to a solver. The QR factorization is exactly that: Gram–Schmidt, packaged as a single product of two matrices.

Given a matrix A whose columns are linearly independent, we can always write

A = QR,

where Q has orthonormal columns (mutually orthogonal, each of length 1) and R is upper triangular with a positive diagonal. Every column of Q is one of the orthonormal vectors Gram–Schmidt would have built; every entry of R is one of the coefficients it computed along the way. Nothing new is happening — it's the same idea wearing matrix clothes.

Where Q and R come from

Run Gram–Schmidt on the columns \mathbf{a}_1, \mathbf{a}_2, \dots, \mathbf{a}_n of A. It produces orthonormal vectors \mathbf{q}_1, \mathbf{q}_2, \dots, \mathbf{q}_n. Stack those as the columns of Q. Because the \mathbf{q}'s span the same space as the \mathbf{a}'s, every original column can be written back in terms of them:

\mathbf{a}_k = r_{1k}\,\mathbf{q}_1 + r_{2k}\,\mathbf{q}_2 + \dots + r_{kk}\,\mathbf{q}_k.

Collect those coefficients r_{jk} into a matrix R and you have A = QR column by column. And because the columns of Q are orthonormal, reading off each coefficient is just a single dot product — r_{jk} = \mathbf{q}_j \cdot \mathbf{a}_k — which is the same as saying

R = Q^{\mathsf{T}} A.

This last identity is the whole trick in a nutshell: Q is the orthonormal basis, and R = Q^{\mathsf{T}}A records how the original columns sit inside it.

Why R is upper triangular

The triangular shape isn't a lucky accident — it's baked into the order Gram–Schmidt works in. The key fact: the k-th column \mathbf{a}_k only ever gets built from the first k orthonormal vectors, because when Gram–Schmidt reached \mathbf{a}_k, only \mathbf{q}_1, \dots, \mathbf{q}_k existed yet. In symbols,

\mathbf{a}_k \in \operatorname{span}(\mathbf{q}_1, \dots, \mathbf{q}_k).

So the coefficient r_{jk} = \mathbf{q}_j \cdot \mathbf{a}_k is zero whenever j > k: a later orthonormal vector \mathbf{q}_j is, by construction, perpendicular to the whole span that \mathbf{a}_k lives in. All the nonzero entries sit on or above the diagonal — that's precisely what "upper triangular" means:

R = \begin{bmatrix} r_{11} & r_{12} & r_{13} \\ 0 & r_{22} & r_{23} \\ 0 & 0 & r_{33} \end{bmatrix}.

The diagonal entries r_{kk} = \lVert \mathbf{u}_k \rVert are the lengths of the leftover vectors before normalising — all positive when the columns of A are independent, which is why we can insist on a positive diagonal and make the factorization unique.

Worked example 1 — factor a 3×2 matrix

Take the matrix whose columns are the same pair Gram–Schmidt handled, \mathbf{a}_1 = (3, 1) and \mathbf{a}_2 = (1, 2):

A = \begin{bmatrix} 3 & 1 \\ 1 & 2 \end{bmatrix}.

Step 1 — build the orthonormal columns. Keep \mathbf{u}_1 = \mathbf{a}_1 = (3,1), with length \sqrt{10}, so

\mathbf{q}_1 = \tfrac{1}{\sqrt{10}}(3, 1) \approx (0.949,\, 0.316).

Subtract \mathbf{a}_2's shadow on \mathbf{u}_1: \operatorname{proj}_{\mathbf{u}_1}(\mathbf{a}_2) = \tfrac{5}{10}(3,1) = (1.5, 0.5), leaving \mathbf{u}_2 = (1,2) - (1.5,0.5) = (-0.5,\,1.5), with length \sqrt{2.5}, so

\mathbf{q}_2 = \tfrac{1}{\sqrt{2.5}}(-0.5, 1.5) \approx (-0.316,\, 0.949).

Step 2 — read off R = Q^{\mathsf{T}}A by dot products. Each entry is a \mathbf{q}_j \cdot \mathbf{a}_k:

A = QR: \quad \begin{bmatrix} 3 & 1 \\ 1 & 2 \end{bmatrix} = \begin{bmatrix} 0.949 & -0.316 \\ 0.316 & 0.949 \end{bmatrix}\begin{bmatrix} 3.162 & 1.581 \\ 0 & 1.581 \end{bmatrix}.

The 0 in the bottom-left of R is \mathbf{q}_2\cdot\mathbf{a}_1 — and it had to vanish, because \mathbf{a}_1 was built before \mathbf{q}_2 ever existed, so it lives entirely in \operatorname{span}(\mathbf{q}_1).

Worked example 2 — QR that solves least-squares

Here is why numerical people love QR. Suppose you want the best-fit solution to an overdetermined system A\mathbf{x} = \mathbf{b} — more equations than unknowns, no exact answer. The classic recipe is the normal equations A^{\mathsf{T}}A\,\hat{\mathbf{x}} = A^{\mathsf{T}}\mathbf{b}, but forming A^{\mathsf{T}}A squares the rounding errors and can be badly behaved. Substitute A = QR instead:

A^{\mathsf{T}}A\,\hat{\mathbf{x}} = A^{\mathsf{T}}\mathbf{b} \;\Longrightarrow\; R^{\mathsf{T}}\!\underbrace{Q^{\mathsf{T}}Q}_{=\,I}R\,\hat{\mathbf{x}} = R^{\mathsf{T}}Q^{\mathsf{T}}\mathbf{b} \;\Longrightarrow\; R\,\hat{\mathbf{x}} = Q^{\mathsf{T}}\mathbf{b}.

The Q^{\mathsf{T}}Q = I cancels, and the whole least-squares problem collapses to R\hat{\mathbf{x}} = Q^{\mathsf{T}}\mathbf{b}. Because R is upper triangular, you don't even need to invert anything — just back-substitute from the bottom row up. Using our factored A with \mathbf{b} = (4, 4):

Q^{\mathsf{T}}\mathbf{b} = \begin{bmatrix} 0.949 & 0.316 \\ -0.316 & 0.949 \end{bmatrix}\begin{bmatrix} 4 \\ 4 \end{bmatrix} \approx \begin{bmatrix} 5.06 \\ 2.53 \end{bmatrix}.

Now solve R\hat{\mathbf{x}} = (5.06, 2.53) from the bottom up: the last row gives 1.581\,\hat{x}_2 = 2.53, so \hat{x}_2 = 1.6; the first row 3.162\,\hat{x}_1 + 1.581(1.6) = 5.06 gives \hat{x}_1 = 0.8. (Since this A is square and invertible, that's the exact solution — for a genuinely taller A the very same three lines return the least-squares best fit.)

The other superstar use of QR is the QR algorithm — the workhorse behind almost every eigenvalue routine you'll ever call. The idea is startlingly simple: factor A = Q_1 R_1, then multiply the pieces back in the opposite order to get A_1 = R_1 Q_1; factor that as Q_2 R_2, form A_2 = R_2 Q_2, and repeat. Each A_k is similar to A (same eigenvalues), and under mild conditions the sequence marches toward an upper-triangular matrix — whose diagonal, read straight off, is the list of eigenvalues. No characteristic polynomial, no root-finding: just Gram–Schmidt, over and over.

Why bother turning a clean process into a matrix equation at all? Because once Q and R are stored, they can be reused. Need to solve A\mathbf{x} = \mathbf{b} for a dozen different right-hand sides \mathbf{b}? Factor once, then each new solve is just a cheap Q^{\mathsf{T}}\mathbf{b} followed by a back-substitution. And in practice the factorization is usually computed not by literal Gram–Schmidt (which loses accuracy when columns are nearly parallel) but by Householder reflections — same A = QR at the end, built more stably. The packaging is the point: a reusable object beats a one-shot procedure.