31. QR decomposition

§20.5's normal equations need A^TA — and forming that product squares A's condition number (§22.2 makes this precise), amplifying numerical error. This closing lesson of Module 20 packages §20.3–20.4's Gram-Schmidt process into a matrix factorization that solves the same least-squares problem without ever computing A^TA.

The factorization

For an m\times n matrix A with independent columns, QR decomposition writes

A=QR

where Q is m\times n with orthonormal columns, and R is n\times n upper triangular. Built via Gram-Schmidt as below, R's diagonal comes out positive automatically (each entry is a norm); this is the standard mathematical convention, and the one the worked example below follows. Library routines are free to return either sign per column — numpy.linalg.qr uses Householder reflections internally and does not guarantee a positive diagonal, so a code-computed Q,R may differ from the Gram-Schmidt version by a sign flip in matching columns of Q and rows of R (harmless: QR still reconstructs A either way).

Construction, directly from Gram-Schmidt: run §20.3–20.4's process on A's columns \vec a_1,\dots,\vec a_n to get orthonormal \vec q_1,\dots,\vec q_n — these are Q's columns. Each \vec a_j is then a linear combination of \vec q_1,\dots,\vec q_j only (never later \vec q_i's, since Gram-Schmidt builds \vec q_j from \vec a_j and earlier vectors alone) — precisely §20.4's coordinate formula, \vec a_j=\sum_{i=1}^j(\vec a_j\cdot\vec q_i)\vec q_i. Stacking these coefficients as columns produces exactly the upper-triangular R, with R_{ij}=\vec a_j\cdot\vec q_i for i\le j and zero below the diagonal.

Solving least squares via QR

Substitute A=QR into §20.5's normal equations:

A^TA\hat{\vec x}=A^T\vec b\ \Longrightarrow\ (QR)^T(QR)\hat{\vec x}=(QR)^T\vec b\ \Longrightarrow\ R^TQ^TQR\hat{\vec x}=R^TQ^T\vec b

Q^TQ=I exactly because Q's columns are orthonormal — §20.4's P^{-1}=P^T identity, applied here even though Q itself may not be square. This cancels Q^TQ entirely:

R^TR\hat{\vec x}=R^TQ^T\vec b\ \Longrightarrow\ R\hat{\vec x}=Q^T\vec b

(cancel R^T from both sides, valid since R is invertible — triangular with nonzero diagonal). R\hat{\vec x}=Q^T\vec b solves by simple back substitution (§16.5's triangular-solve trick, reused directly) — A^TA never gets formed at all, sidestepping the numerical amplification it causes, which is exactly why production least-squares solvers (including numpy.linalg.lstsq) use QR (or the closely related SVD, §21.4) rather than the normal equations directly.

Doing it in Python

import numpy as np

A = np.array([[1., 1.], [1., 2.], [1., 3.]])
Q, R = np.linalg.qr(A)

print("Q =")
for row in Q.tolist():
    print([round(v, 4) for v in row])
print("R =")
for row in R.tolist():
    print([round(v, 4) for v in row])

print("\nQ^T Q == I:", bool(np.allclose(Q.T @ Q, np.eye(2))))
print("Q @ R == A:", bool(np.allclose(Q @ R, A)))

Q =
[-0.5774, 0.7071]
[-0.5774, 0.0]
[-0.5774, -0.7071]
R =
[-1.7321, -3.4641]
[0.0, -1.4142]

Q^T Q == I: True
Q @ R == A: True

Solving least squares via R\hat{\vec x}=Q^T\vec b, and confirming it matches §20.5's normal-equations answer:

import numpy as np
import scipy.linalg as la

A = np.array([[1., 1.], [1., 2.], [1., 3.]])
b = np.array([2., 3., 5.])

Q, R = np.linalg.qr(A)
x_qr = la.solve_triangular(R, Q.T @ b)

x_normal = np.linalg.solve(A.T @ A, A.T @ b)

print("via QR:             ", [round(v, 4) for v in x_qr.tolist()])
print("via normal equations:", [round(v, 4) for v in x_normal.tolist()])
print("match:", bool(np.allclose(x_qr, x_normal)))

via QR:              [0.3333, 1.5]
via normal equations: [0.3333, 1.5]
match: True

Worked example

Find the QR decomposition of A=\begin{pmatrix}1&0\\1&1\end{pmatrix} using Gram-Schmidt directly.

\vec a_1=(1,1), \vec a_2=(0,1). \|\vec a_1\|=\sqrt2, so \vec q_1=(1,1)/\sqrt2.

\operatorname{proj}_{\vec q_1}\vec a_2=(\vec a_2\cdot\vec q_1)\vec q_1=\left(\frac1{\sqrt2}\right)\frac{(1,1)}{\sqrt2}=\left(\frac12,\frac12\right). \vec u_2=(0,1)-(\tfrac12,\tfrac12)=(-\tfrac12,\tfrac12), \|\vec u_2\|=\tfrac1{\sqrt2}, so \vec q_2=(-1,1)/\sqrt2.

Q=\frac1{\sqrt2}\begin{pmatrix}1&-1\\1&1\end{pmatrix}

R_{11}=\vec a_1\cdot\vec q_1=\sqrt2 (equivalently just \|\vec a_1\|). R_{12}=\vec a_2\cdot\vec q_1=\tfrac1{\sqrt2}. R_{22}=\|\vec u_2\|=\tfrac1{\sqrt2}. R_{21}=0 (upper triangular, by construction).

\boxed{R=\begin{pmatrix}\sqrt2&\tfrac1{\sqrt2}\\0&\tfrac1{\sqrt2}\end{pmatrix}}

Sanity check. QR=\frac1{\sqrt2}\begin{pmatrix}1&-1\\1&1\end{pmatrix}\begin{pmatrix}\sqrt2&1/\sqrt2\\0&1/\sqrt2\end{pmatrix}=\frac1{\sqrt2}\begin{pmatrix}\sqrt2&1/\sqrt2-1/\sqrt2\\\sqrt2&1/\sqrt2+1/\sqrt2\end{pmatrix}=\frac1{\sqrt2}\begin{pmatrix}\sqrt2&0\\\sqrt2&\sqrt2\end{pmatrix}=\begin{pmatrix}1&0\\1&1\end{pmatrix}=A ✓ — reconstructs A exactly, and Q^TQ=\frac12\begin{pmatrix}1&1\\-1&1\end{pmatrix}\begin{pmatrix}1&-1\\1&1\end{pmatrix}=\frac12\begin{pmatrix}2&0\\0&2\end{pmatrix}=I ✓, confirming Q's columns really are orthonormal.

Your turn

1. Why must R's diagonal entries be positive (in the standard convention)?

2. If A is already n\times n with orthonormal columns, what are Q and R?

3. True or false: QR decomposition requires A to be square.

Solutions

1. R_{jj}=\|\vec u_j\| (§20.3's Gram-Schmidt vector, before normalizing) — a norm, always nonnegative, and strictly positive whenever \vec u_j\neq\vec0, which holds for every j precisely because A's columns are independent (§20.3's nonzero-output guarantee). A zero diagonal entry would mean a dependent column, contradicting the independence assumption QR relies on.

2. Q=A itself (its columns are already orthonormal, so no Gram-Schmidt adjustment is needed) and R=I — the trivial factorization A=AI, matching the fact that Gram-Schmidt applied to an already-orthonormal set changes nothing (§20.3's "your turn" question 3, one level further: orthonormal, not just orthogonal).

3. False. QR is defined for any m\times n matrix with independent columns, m\ge n — tall or square, never wide with m<n (which can't have n independent columns in only m<n dimensions to begin with, by §17.1). This is precisely the shape §20.5's least squares needs: more equations than unknowns.

Check yourself in code

Compute the QR decomposition of A=\begin{pmatrix}0&1\\1&1\\1&0\end{pmatrix} and confirm QR=A and Q^TQ=I.

Print exactly this:

Q @ R == A: True
Q^T Q == I: True

import numpy as np

A = np.array([[0., 1.], [1., 1.], [1., 0.]])
Q, R = np.linalg.qr(A)

print("Q @ R == A:", bool(np.allclose(Q @ R, A)))
# print whether Q^T Q equals the 2x2 identity

import numpy as np

A = np.array([[0., 1.], [1., 1.], [1., 0.]])
Q, R = np.linalg.qr(A)

print("Q @ R == A:", bool(np.allclose(Q @ R, A)))
print("Q^T Q == I:", bool(np.allclose(Q.T @ Q, np.eye(2))))

QR decomposition, A=QR with orthonormal Q and upper-triangular R, is Gram-Schmidt packaged as a matrix factorization. It solves least-squares via R\hat{\vec x}=Q^T\vec b (simple back substitution) without ever forming A^TA, which is why real solvers use it instead of §20.5's normal equations directly.

This closes Module 20. Every idea here — norm, orthogonality, projection — required only the dot product's three axioms. Module 21 asks what happens when a matrix's own structure (specifically, symmetry) interacts with this geometry, starting with a genuinely beautiful fact: a symmetric matrix's eigenvectors are always orthogonal.