The Gram–Schmidt Process: Orthogonality from Scratch
Given any basis for a subspace, Gram–Schmidt produces an orthonormal basis spanning exactly the same subspace. The algorithm is constructive, finite, and geometrically transparent: at each step, you subtract off the components that lie in directions you’ve already handled, and normalize what remains. The residual points in a genuinely new direction—perpendicular to everything before it.
This article develops the process from the ground up: the geometry of projection, the algorithm itself with full proofs, worked examples in $\mathbb{R}^2$ and $\mathbb{R}^3$ with visual diagrams, the connection to QR factorization, the numerical instability of the classical algorithm and its fix, and applications across mathematics, physics, and computation.
1. Why Orthogonality Matters
An orthonormal basis $\{e_1, \ldots, e_n\}$ satisfies $\langle e_i, e_j \rangle = \delta_{ij}$—each vector has unit length, and distinct vectors are perpendicular. This is a strong condition, and it buys you enormous computational and theoretical leverage.
Coordinates are free. If $\{e_1, \ldots, e_n\}$ is orthonormal, then the coordinates of any vector $v$ are simply inner products:
$$v = \sum_{i=1}^{n} \langle v, e_i \rangle \, e_i$$No matrix inversion. No solving systems of equations. You just project.
Lengths decompose cleanly. By the Pythagorean theorem in $n$ dimensions:
$$\|v\|^2 = \sum_{i=1}^{n} |\langle v, e_i \rangle|^2$$This is Parseval’s identity, and it says that the total “energy” in a vector equals the sum of squared energies in each orthogonal component. No cross terms. No interference.
Matrices become simple. The matrix representing a linear transformation with respect to an orthonormal basis has the property that its inverse equals its transpose: $Q^{-1} = Q^T$ (or $Q^* = Q^{-1}$ in the complex case). These orthogonal (or unitary) matrices preserve lengths and angles—they are the rigid motions of linear algebra.
The core promise
Orthogonal bases turn hard problems (solving systems, minimizing errors, decomposing signals) into easy ones (taking inner products). Gram–Schmidt is the machine that manufactures these bases from arbitrary inputs.
2. Orthogonal Projection: The Geometric Engine
The entire Gram–Schmidt process rests on one operation: orthogonal projection. Before we state the algorithm, we need to understand this operation geometrically and algebraically.
Projection onto a single vector
Given a nonzero vector $u$ and any vector $v$, the orthogonal projection of $v$ onto $u$ is the vector in the direction of $u$ that is closest to $v$:
$$\operatorname{proj}_u v = \frac{\langle v, u \rangle}{\langle u, u \rangle} \, u$$
The geometric interpretation: drop a perpendicular from the tip of $v$ down to the line spanned by $u$. The foot of that perpendicular is $\operatorname{proj}_u v$. The vector from the projection to $v$—namely $v - \operatorname{proj}_u v$—is orthogonal to $u$.
The key algebraic fact: the residual is orthogonal to $u$.
Verification
Compute:
$$\left\langle v - \frac{\langle v, u \rangle}{\langle u, u \rangle} u, \; u \right\rangle = \langle v, u \rangle - \frac{\langle v, u \rangle}{\langle u, u \rangle} \langle u, u \rangle = \langle v, u \rangle - \langle v, u \rangle = 0 \quad \checkmark$$$\blacksquare$
Projection onto a subspace
If $\{u_1, \ldots, u_k\}$ is an orthogonal set (pairwise perpendicular, but not necessarily unit length), then the projection of $v$ onto $W = \operatorname{span}\{u_1, \ldots, u_k\}$ decomposes into independent projections onto each direction:
$$\operatorname{proj}_W v = \sum_{i=1}^{k} \frac{\langle v, u_i \rangle}{\langle u_i, u_i \rangle} \, u_i$$
This is the key simplification that orthogonality provides. For a non-orthogonal basis, projecting onto a subspace requires solving a system of equations (the normal equations). For an orthogonal basis, you just sum the individual projections—they don’t interfere with each other because the cross terms vanish.
3. The Gram–Schmidt Algorithm
Now we have the tool. Here is the machine that uses it.
Input: A linearly independent set $\{v_1, v_2, \ldots, v_n\}$ in an inner product space.
Output: An orthogonal set $\{u_1, u_2, \ldots, u_n\}$ such that $\operatorname{span}\{u_1, \ldots, u_k\} = \operatorname{span}\{v_1, \ldots, v_k\}$ for each $k$.
Classical Gram–Schmidt (CGS)
- Set $u_1 = v_1$
- For $k = 2, 3, \ldots, n$:
$$u_k = v_k - \sum_{j=1}^{k-1} \frac{\langle v_k, u_j \rangle}{\langle u_j, u_j \rangle} \, u_j$$ - To normalize: set $e_k = u_k / \|u_k\|$ for each $k$.
The idea at each step is the same: take the next input vector $v_k$, subtract off its projections onto all previously computed orthogonal directions, and keep the residual. The residual is perpendicular to everything before it by construction—we subtracted away every component that wasn’t.
The geometric intuition
Think of building a coordinate system one axis at a time. The first axis is easy—just take $v_1$. For the second axis, you need a direction perpendicular to the first. Take $v_2$ and remove its shadow on $u_1$; what remains points in a genuinely new direction. For the third axis, take $v_3$ and remove its shadows on both $u_1$ and $u_2$; the residual is perpendicular to the entire plane spanned by the first two. And so on.
4. The Span Preservation Property
A subtle but crucial property: at every stage $k$, the span of the first $k$ orthogonal vectors equals the span of the first $k$ input vectors. This is not obvious, and it matters.
Theorem (Span Preservation)
For each $k = 1, \ldots, n$: $\;\operatorname{span}\{u_1, \ldots, u_k\} = \operatorname{span}\{v_1, \ldots, v_k\}$.
Proof by induction on $k$.
Base case ($k = 1$): $u_1 = v_1$, so the spans are equal.
Inductive step: Assume $\operatorname{span}\{u_1, \ldots, u_{k-1}\} = \operatorname{span}\{v_1, \ldots, v_{k-1}\}$. By the formula:
$$u_k = v_k - \sum_{j=1}^{k-1} \frac{\langle v_k, u_j \rangle}{\langle u_j, u_j \rangle} u_j$$Since $u_1, \ldots, u_{k-1} \in \operatorname{span}\{v_1, \ldots, v_{k-1}\}$, the sum is in $\operatorname{span}\{v_1, \ldots, v_{k-1}\}$. So $u_k = v_k - (\text{something in } \operatorname{span}\{v_1, \ldots, v_{k-1}\})$, which means $u_k \in \operatorname{span}\{v_1, \ldots, v_k\}$.
Conversely, rearranging: $v_k = u_k + \sum_{j=1}^{k-1} c_j u_j$, so $v_k \in \operatorname{span}\{u_1, \ldots, u_k\}$. Combined with the inductive hypothesis, both spans are equal.
Moreover, $u_k \neq 0$ because $v_k \notin \operatorname{span}\{v_1, \ldots, v_{k-1}\}$ (the $v_i$ are linearly independent).
$\blacksquare$
This property is what makes Gram–Schmidt useful beyond abstract orthogonalization. It means the process doesn’t “drift” into a different subspace. The orthonormal basis spans exactly the same space as the original basis. This is essential for the QR factorization, for least-squares problems, and for any application where you need an orthonormal basis for a specific subspace.
5. Worked Example: $\mathbb{R}^2$
Let $v_1 = \begin{pmatrix} 3 \\ 1 \end{pmatrix}$ and $v_2 = \begin{pmatrix} 2 \\ 4 \end{pmatrix}$. These are linearly independent (neither is a scalar multiple of the other), but they are not orthogonal: $\langle v_1, v_2 \rangle = 6 + 4 = 10 \neq 0$.
Step 1: Set $u_1 = v_1$
$$u_1 = \begin{pmatrix} 3 \\ 1 \end{pmatrix}$$Step 2: Compute $u_2$
Subtract the projection of $v_2$ onto $u_1$:
$$\operatorname{proj}_{u_1} v_2 = \frac{\langle v_2, u_1 \rangle}{\langle u_1, u_1 \rangle} u_1 = \frac{10}{10} \begin{pmatrix} 3 \\ 1 \end{pmatrix} = \begin{pmatrix} 3 \\ 1 \end{pmatrix}$$ $$u_2 = v_2 - \operatorname{proj}_{u_1} v_2 = \begin{pmatrix} 2 \\ 4 \end{pmatrix} - \begin{pmatrix} 3 \\ 1 \end{pmatrix} = \begin{pmatrix} -1 \\ 3 \end{pmatrix}$$Verification: $\langle u_1, u_2 \rangle = (3)(-1) + (1)(3) = -3 + 3 = 0$. $\checkmark$
Step 3: Normalize
$$e_1 = \frac{u_1}{\|u_1\|} = \frac{1}{\sqrt{10}} \begin{pmatrix} 3 \\ 1 \end{pmatrix}, \qquad e_2 = \frac{u_2}{\|u_2\|} = \frac{1}{\sqrt{10}} \begin{pmatrix} -1 \\ 3 \end{pmatrix}$$6. Worked Example: $\mathbb{R}^3$
Let $v_1 = \begin{pmatrix} 1 \\ 1 \\ 0 \end{pmatrix}$, $v_2 = \begin{pmatrix} 1 \\ 0 \\ 1 \end{pmatrix}$, $v_3 = \begin{pmatrix} 0 \\ 1 \\ 1 \end{pmatrix}$. These span all of $\mathbb{R}^3$ but are not orthogonal.
Check: $\langle u_1, u_2 \rangle = \frac{1}{2} - \frac{1}{2} + 0 = 0$. $\checkmark$
$\|u_2\|^2 = \frac{1}{4} + \frac{1}{4} + 1 = \frac{3}{2}$
Verification:
- $\langle u_1, u_3 \rangle = -\frac{2}{3} + \frac{2}{3} + 0 = 0$ $\checkmark$
- $\langle u_2, u_3 \rangle = \frac{1}{2}(-\frac{2}{3}) + (-\frac{1}{2})(\frac{2}{3}) + 1 \cdot \frac{2}{3} = -\frac{1}{3} - \frac{1}{3} + \frac{2}{3} = 0$ $\checkmark$
The orthonormal basis
Normalizing:
$$e_1 = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 \\ 1 \\ 0 \end{pmatrix}, \quad e_2 = \sqrt{\frac{2}{3}} \begin{pmatrix} 1/2 \\ -1/2 \\ 1 \end{pmatrix}, \quad e_3 = \frac{1}{\sqrt{4/3}} \begin{pmatrix} -2/3 \\ 2/3 \\ 2/3 \end{pmatrix} = \frac{1}{\sqrt{3}} \begin{pmatrix} -1 \\ 1 \\ 1 \end{pmatrix}$$7. The QR Factorization
Gram–Schmidt does more than produce an orthonormal basis. It implicitly computes a matrix factorization that is one of the most important in numerical linear algebra.
Let $A = [v_1 \mid v_2 \mid \cdots \mid v_n]$ be the matrix whose columns are the input vectors, and let $Q = [e_1 \mid e_2 \mid \cdots \mid e_n]$ be the matrix whose columns are the orthonormal vectors. Then:
$$A = QR$$
where $Q$ is orthogonal ($Q^TQ = I$) and $R$ is upper triangular with positive diagonal entries:
$$R = \begin{pmatrix} \langle v_1, e_1 \rangle & \langle v_2, e_1 \rangle & \cdots & \langle v_n, e_1 \rangle \\ 0 & \langle v_2, e_2 \rangle & \cdots & \langle v_n, e_2 \rangle \\ \vdots & & \ddots & \vdots \\ 0 & 0 & \cdots & \langle v_n, e_n \rangle \end{pmatrix}$$Why is $R$ upper triangular? Because Gram–Schmidt preserves spans. The vector $e_k$ is orthogonal to $v_1, \ldots, v_{k-1}$’s span, so $\langle v_j, e_k \rangle = 0$ whenever $j < k$. The entries below the diagonal are all zero.
Why are the diagonal entries positive? Because $r_{kk} = \langle v_k, e_k \rangle = \|u_k\|$, the norm of the orthogonal vector before normalization, which is positive for linearly independent inputs.
QR of the $\mathbb{R}^3$ example
From our earlier computation:
$$\underbrace{\begin{pmatrix} 1 & 1 & 0 \\ 1 & 0 & 1 \\ 0 & 1 & 1 \end{pmatrix}}_A = \underbrace{\begin{pmatrix} 1/\sqrt{2} & 1/\sqrt{6} & -1/\sqrt{3} \\ 1/\sqrt{2} & -1/\sqrt{6} & 1/\sqrt{3} \\ 0 & 2/\sqrt{6} & 1/\sqrt{3} \end{pmatrix}}_Q \underbrace{\begin{pmatrix} \sqrt{2} & 1/\sqrt{2} & 1/\sqrt{2} \\ 0 & \sqrt{3/2} & 1/\sqrt{6} \\ 0 & 0 & 2/\sqrt{3} \end{pmatrix}}_R$$The QR factorization is used everywhere: solving least-squares problems, computing eigenvalues (the QR algorithm), solving linear systems more stably than LU decomposition, and as a building block for SVD computation.
8. Numerical Stability: Classical vs. Modified Gram–Schmidt
The classical algorithm (CGS) described above is mathematically correct but numerically dangerous. In floating-point arithmetic, the computed vectors $u_k$ gradually lose orthogonality due to rounding errors. For ill-conditioned matrices, the loss can be catastrophic.
The problem is that CGS computes all projections using the original $v_k$. If the early projections introduce rounding errors into partial results, those errors compound. The fix is simple: project sequentially.
Modified Gram–Schmidt (MGS)
- Set $u_1 = v_1$, $\;e_1 = u_1 / \|u_1\|$
- For $k = 2, 3, \ldots, n$:
- Set $u_k^{(0)} = v_k$
- For $j = 1, 2, \ldots, k-1$:
- $u_k^{(j)} = u_k^{(j-1)} - \langle u_k^{(j-1)}, e_j \rangle \, e_j$
- Set $u_k = u_k^{(k-1)}$, $\;e_k = u_k / \|u_k\|$
The difference: instead of computing all projections from the original $v_k$, MGS subtracts one projection at a time, updating the intermediate vector after each subtraction. Mathematically, CGS and MGS produce identical results. Numerically, MGS is dramatically more stable.
The difference, intuitively
Classical GS: “Compute all the shadows at once and subtract them all from the original vector.”
Modified GS: “Subtract the first shadow. Then, from what remains, subtract the second shadow. Then from that, subtract the third. And so on.”
Modified GS is more accurate because each subsequent projection operates on a vector that has already been cleaned of its earlier components. Rounding errors in the first projection don’t contaminate the second projection’s input.
| Property | Classical GS | Modified GS | Householder QR |
|---|---|---|---|
| FLOPS for $m \times n$ matrix | $2mn^2$ | $2mn^2$ | $2mn^2 - \frac{2}{3}n^3$ |
| Orthogonality of computed $Q$ | $\|Q^TQ - I\| \sim \kappa(A) \cdot \epsilon$ | $\|Q^TQ - I\| \sim \kappa(A) \cdot \epsilon$ | $\|Q^TQ - I\| \sim \epsilon$ |
| Residual $\|A - QR\|$ | $\sim \epsilon \|A\|$ | $\sim \epsilon \|A\|$ | $\sim \epsilon \|A\|$ |
| Produces $Q$ explicitly? | Yes | Yes | Only as a product of reflectors |
Here $\kappa(A)$ is the condition number and $\epsilon$ is machine epsilon ($\approx 10^{-16}$ in double precision). For well-conditioned matrices ($\kappa \approx 1$), CGS and MGS are equally fine. For ill-conditioned matrices ($\kappa \gg 1$), CGS can produce vectors that are far from orthogonal, while MGS degrades more gracefully. For the best numerical orthogonality, Householder reflections are preferred, but MGS is often sufficient and has the advantage of producing the columns of $Q$ one at a time.
9. Gram–Schmidt in General Inner Product Spaces
Nothing in the algorithm requires $\mathbb{R}^n$. The Gram–Schmidt process works in any inner product space: $\mathbb{C}^n$ with the Hermitian inner product, function spaces with $\langle f, g \rangle = \int_a^b f(x) \overline{g(x)} \, dx$, or any space where you have a notion of length and angle.
Example: Orthogonal polynomials
Start with the monomials $\{1, x, x^2, x^3, \ldots\}$ in the space of polynomials on $[-1, 1]$ with the inner product:
$$\langle f, g \rangle = \int_{-1}^{1} f(x) g(x) \, dx$$Applying Gram–Schmidt produces the Legendre polynomials:
$$P_0(x) = 1, \quad P_1(x) = x, \quad P_2(x) = \tfrac{1}{2}(3x^2 - 1), \quad P_3(x) = \tfrac{1}{2}(5x^3 - 3x), \;\ldots$$These are orthogonal: $\int_{-1}^{1} P_m(x) P_n(x) \, dx = 0$ for $m \neq n$.
Different inner products on the same underlying space of polynomials produce different orthogonal families:
| Inner product | Resulting orthogonal polynomials |
|---|---|
| $\int_{-1}^{1} f \, g \, dx$ | Legendre polynomials |
| $\int_{-1}^{1} \frac{f \, g}{\sqrt{1 - x^2}} \, dx$ | Chebyshev polynomials |
| $\int_{0}^{\infty} f \, g \, e^{-x} \, dx$ | Laguerre polynomials |
| $\int_{-\infty}^{\infty} f \, g \, e^{-x^2} \, dx$ | Hermite polynomials |
All of these are Gram–Schmidt applied to $\{1, x, x^2, \ldots\}$ with different weight functions. The same algorithm, different geometries, different results.
10. Applications
Least-squares approximation
Given an overdetermined system $Ax = b$ (more equations than unknowns), the least-squares solution minimizes $\|Ax - b\|^2$. Via the QR factorization $A = QR$:
$$Ax = b \;\implies\; QRx = b \;\implies\; Rx = Q^T b$$Since $R$ is upper triangular, this is trivially solved by back-substitution. No need to form $A^T A$ (which squares the condition number) or invert anything. This is why QR-based least squares is the standard in practice.
Signal processing and Fourier analysis
The Fourier series is an orthogonal expansion. The functions $\{1, \cos x, \sin x, \cos 2x, \sin 2x, \ldots\}$ are orthogonal on $[-\pi, \pi]$ with respect to $\langle f, g \rangle = \int_{-\pi}^{\pi} f(x)g(x) \, dx$. The Fourier coefficients are just inner products—projections onto orthogonal directions—exactly because the basis is orthogonal. This is Gram–Schmidt’s payoff: once you have an orthogonal basis, decomposition is projection.
Quantum mechanics
In quantum mechanics, states live in a Hilbert space (a complete inner product space). Observable quantities correspond to self-adjoint operators, and their eigenstates form an orthonormal basis for the space. Gram–Schmidt is used to orthonormalize degenerate eigenstates (eigenstates with the same eigenvalue that are not automatically orthogonal). The Born rule—the probability of measuring a state $|\psi\rangle$ in eigenstate $|e_k\rangle$ is $|\langle e_k | \psi \rangle|^2$—works precisely because the eigenstates are orthonormal.
Krylov methods and iterative solvers
The Arnoldi iteration (for non-symmetric matrices) and the Lanczos algorithm (for symmetric matrices) are variants of Gram–Schmidt applied to Krylov subspaces $\{b, Ab, A^2b, \ldots\}$. These generate orthonormal bases for the subspaces in which iterative methods like GMRES and conjugate gradients search for approximate solutions to $Ax = b$. Without Gram–Schmidt (or Householder) orthogonalization, these methods would not be numerically stable.
Data science and statistics
In regression, orthogonalizing the predictor variables via Gram–Schmidt reveals the sequential contribution of each variable. The residual at step $k$ shows what the $k$-th predictor adds that the previous $k-1$ predictors cannot explain. This is the geometric foundation of variable importance and sequential ANOVA.
11. The View from Above
Gram–Schmidt is a theorem about decomposition. It says that any finite-dimensional inner product space has an orthonormal basis, and provides an explicit construction. More precisely:
Every finite-dimensional inner product space has an orthonormal basis.
Every linearly independent set can be orthogonalized while preserving its span.
Every matrix with linearly independent columns has a unique QR factorization with $R$ having positive diagonal.
These are three phrasings of the same theorem, and Gram–Schmidt is the constructive proof of all three.
The process is named after Jørgen Pedersen Gram (1879) and Erhard Schmidt (1907), though the essential idea appeared earlier in the work of Laplace and Cauchy. Schmidt’s contribution was to formalize it for infinite-dimensional function spaces, where the convergence questions are more delicate.
At its heart, the algorithm is almost embarrassingly simple: subtract shadows, keep what’s left. But this simple operation—applied iteratively—manufactures orthogonality from arbitrary inputs. It turns messy coordinates into clean ones, complicated projections into trivial ones, and unstable computations into stable ones. It is the bridge between the bases you are given and the bases you want.