The Problem: Building an Orthogonal Basis
In the previous lesson we saw how powerful orthogonal sets of vectors are — projections simplify, computations stabilize, and formulas become elegant. But in practice, the vectors you start with are rarely orthogonal. You might have three directions that span a subspace, but they point at oblique angles to each other.
The Gram-Schmidt process solves this problem: given any basis {v₁, v₂, …, vₙ} for a subspace, it produces an orthonormal basis {q₁, q₂, …, qₙ} for the same subspace. The key idea is sequential projection — each new basis vector is obtained by taking the next input vector and subtracting off all components in the directions already found.
An orthogonal basis has mutually perpendicular vectors (dot products are zero). An orthonormal basis goes further — each vector also has unit length. Orthonormality is the gold standard: when columns of a matrix Q are orthonormal, QTQ = I, which makes every formula involving Q fast and numerically stable.
The Algorithm, Step by Step
Start with linearly independent vectors v₁, v₂, …, vₙ. Gram-Schmidt produces orthonormal vectors q₁, q₂, …, qₙ in three conceptual steps per iteration:
- Start fresh — take the next input vector vₖ.
- Subtract projections — remove the components of vₖ in each direction already established (q₁ through qₖ₋₁).
- Normalize — divide the result by its length to get a unit vector.
The projection of vₖ onto qⱼ is simply (qⱼᵀvₖ)qⱼ, because qⱼ is already a unit vector. The scalar qⱼᵀvₖ is the inner product (or correlation) between vₖ and the direction qⱼ — it measures how much of vₖ lies along qⱼ. By subtracting all such components, the residual uₖ is guaranteed to be orthogonal to all previous q vectors.
The First Two Steps in Detail
Step 1 (k = 1): The first vector is trivial — just normalize v₁. Set u₁ = v₁ and q₁ = u₁ / ‖u₁‖. There is nothing to subtract yet.
Step 2 (k = 2): Take v₂ and subtract its projection onto q₁. The component of v₂ along q₁ is (q₁ᵀv₂)q₁, so:
Check: q₁ᵀu₂ = q₁ᵀv₂ − (q₁ᵀv₂)(q₁ᵀq₁) = q₁ᵀv₂ − q₁ᵀv₂ = 0. The subtraction perfectly cancels the q₁ component — u₂ is orthogonal to q₁ by construction.
A Worked Example
Let v₁ = [1, 1, 0]ᵀ and v₂ = [1, 0, 1]ᵀ and v₃ = [0, 1, 1]ᵀ in ℝ³. These are linearly independent (check: their determinant is nonzero). Apply Gram-Schmidt:
Step 1: u₁ = v₁ = [1, 1, 0]ᵀ. Length: ‖u₁‖ = √2. So q₁ = [1/√2, 1/√2, 0]ᵀ.
Step 2: q₁ᵀv₂ = (1/√2)(1) + (1/√2)(0) + 0 = 1/√2. Subtract: u₂ = v₂ − (1/√2)q₁ = [1, 0, 1]ᵀ − (1/√2)[1/√2, 1/√2, 0]ᵀ = [1, 0, 1]ᵀ − [1/2, 1/2, 0]ᵀ = [1/2, −1/2, 1]ᵀ. Length: ‖u₂‖ = √(1/4 + 1/4 + 1) = √(3/2) = √6/2. So q₂ = [1/√6, −1/√6, 2/√6]ᵀ.
Step 3: q₁ᵀv₃ = (1/√2)(0) + (1/√2)(1) + 0 = 1/√2. q₂ᵀv₃ = (1/√6)(0) + (−1/√6)(1) + (2/√6)(1) = 1/√6. Subtract: u₃ = v₃ − (1/√2)q₁ − (1/√6)q₂ = [0,1,1]ᵀ − [1/2,1/2,0]ᵀ − [1/6,−1/6,2/6]ᵀ = [−2/3, 2/3, 2/3]ᵀ. After normalizing, q₃ = [−1/√3, 1/√3, 1/√3]ᵀ.
Verify: q₁·q₂ = (1/√2)(1/√6) + (1/√2)(−1/√6) + 0 = 0 ✓. And q₁·q₃ = (1/√2)(−1/√3) + (1/√2)(1/√3) + 0 = 0 ✓. Each computed qⱼ has unit length and each pair is orthogonal — exactly what Gram-Schmidt guarantees.
The QR Decomposition
Gram-Schmidt does more than produce an orthonormal basis — it secretly factors the original matrix. Arrange the input vectors as columns of A = [v₁ v₂ … vₙ]. The Gram-Schmidt process yields Q = [q₁ q₂ … qₙ] with orthonormal columns. What happened to A?
Since each vₖ is a linear combination of q₁, …, qₖ (by construction), the matrix A = QR where R is upper triangular with positive diagonal entries rₖₖ = ‖uₖ‖. This is the QR decomposition:
The entry rᵢⱼ = qᵢᵀvⱼ is the inner product computed during Gram-Schmidt — the amount of vⱼ that projects onto qᵢ. Because Gram-Schmidt only subtracts projections onto earlier basis vectors (q₁, …, qₖ₋₁), each vₖ only involves q₁ through qₖ — making R upper triangular rather than full. This is not a coincidence; it is the algebraic signature of the sequential subtraction process.
Why QR Is Fundamental
The QR decomposition is one of the most important matrix factorizations in numerical linear algebra. It is used to:
- Solve least squares problems stably (via back-substitution on R)
- Compute eigenvalues (the QR algorithm iterates QR factorizations)
- Orthogonalize vectors in Krylov subspace methods (GMRES, Lanczos)
- Implement the Fourier transform efficiently
Numerical Stability: Classical vs. Modified Gram-Schmidt
The classical Gram-Schmidt algorithm described above is mathematically correct, but in floating-point arithmetic it can accumulate errors. The issue: when you compute u₂ = v₂ − (q₁ᵀv₂)q₁ and then u₃ = v₃ − (q₁ᵀv₃)q₁ − (q₂ᵀv₃)q₂ using a pre-computed q₂, rounding errors in q₂ propagate into u₃.
Modified Gram-Schmidt (MGS) reorders the same computation to reduce error growth. Instead of subtracting all projections at once using the original vₖ, MGS updates a running vector after each projection subtraction:
Both algorithms produce the same result in exact arithmetic. The difference is purely about floating-point behavior. MGS keeps the intermediate vectors closer to orthogonal because it re-orthogonalizes against the already-cleaned-up directions rather than the original (potentially correlated) vᵢ. For nearly linearly dependent vectors, the difference can be dramatic.
Orthonormal Columns and the Projection Formula
One of the big payoffs of Gram-Schmidt is the simplification of the projection formula. Recall from Lesson 5.1: the projection onto col(A) uses P = A(AᵀA)⁻¹Aᵀ, which requires computing and inverting AᵀA. When A = Q has orthonormal columns, QᵀQ = I (the identity), so the formula collapses:
This is the computational power of Gram-Schmidt: by investing the effort upfront to orthonormalize a basis, every subsequent projection, least squares solve, or coordinate extraction becomes a simple dot product — no solving systems, no inverting matrices.
Connection to the Fourier Series
The Fourier series is precisely this formula applied to an infinite orthonormal basis of sines and cosines. The Fourier coefficient aₙ = ⟨f, cos(nθ)⟩ is the inner product of f with the n-th basis function — the analog of qⱼᵀb. The Fourier series then reconstructs f as a sum of these projections. Gram-Schmidt, in its most abstract form, is the finite-dimensional version of Fourier analysis.
Gram-Schmidt converts any linearly independent set {v₁, …, vₙ} into an orthonormal set {q₁, …, qₙ} spanning the same subspace. Each new vector subtracts projections onto all previous directions: uₖ = vₖ − Σⱼ (qⱼᵀvₖ)qⱼ, then qₖ = uₖ/‖uₖ‖. Applied to the columns of A, Gram-Schmidt produces the QR decomposition A = QR with QᵀQ = I and R upper triangular. The Modified Gram-Schmidt algorithm is numerically preferred. Once orthonormal, projections reduce to QQᵀb — no matrix inversion needed, just dot products.