Gram–Schmidt Orthonormalization

Derivation

Gram–Schmidt Orthonormalization

D-075 Home PU-104 Threads symmetry · waves Depends on Invariance of Dimension of a Vector Space, Cauchy–Schwarz and the Triangle Inequality

Statement

Let \((V,\langle\cdot,\cdot\rangle)\) be an inner-product space over \(F=\mathbb{R}\) or \(\mathbb{C}\), and let \(\{v_1,\dots,v_n\}\subset V\) be linearly independent. Then there exists an orthonormal set \(\{e_1,\dots,e_n\}\) — meaning \(\langle e_i,e_j\rangle=\delta_{ij}\) — such that for every \(k\le n\), \(\operatorname{span}\{e_1,\dots,e_k\}=\operatorname{span}\{v_1,\dots,v_k\}\). The set is produced constructively, and is unique up to the phase (sign, in the real case) of each \(e_k\).

Why it matters

An orthonormal basis is the coordinate system in which the inner product becomes trivial: components are read off by projection, \(x=\sum_k\langle e_k,x\rangle\,e_k\), and lengths obey Pythagoras, \(\lVert x\rVert^2=\sum_k|\langle e_k,x\rangle|^2\). Gram–Schmidt is the theorem that says such a basis always exists and, crucially, can be built from any basis without leaving the subspace you started in — the flag \(\operatorname{span}\{v_1,\dots,v_k\}\) is preserved at every stage.

Along the waves thread it is the engine behind orthogonal polynomials (Legendre, Hermite, Chebyshev) and hence Fourier-type expansions; along the symmetry thread it underwrites the \(QR\) decomposition and the existence of orthonormal eigenbases for self-adjoint operators. Nearly every "expand in a basis" argument in physics silently assumes the basis is orthonormal, which is a promise Gram–Schmidt keeps.

Assumptions

The set \(\{v_1,\dots,v_n\}\) is linearly independent.If some \(v_k\) lies in the span of its predecessors then \(u_k=0\) at that stage and the normalisation \(u_k/\lVert u_k\rVert\) divides by zero. The process still works if you simply discard such vectors, yielding an orthonormal basis of the span, but the clean one-to-one correspondence \(v_k\leftrightarrow e_k\) is lost.

\(\langle\cdot,\cdot\rangle\) is a genuine inner product: positive-definite, so \(\langle u,u\rangle>0\) for \(u\neq0\).Positivity is exactly what makes \(\lVert u_k\rVert=\sqrt{\langle u_k,u_k\rangle}\) a real positive number to divide by. Drop it and the "norm" can vanish or go imaginary on nonzero vectors, and the algorithm produces nonsense (see Breaks when).

Convention: \(\langle\cdot,\cdot\rangle\) is conjugate-linear in the first slot and linear in the second (physics convention).The projection coefficient is then \(\langle u_j,v_k\rangle/\langle u_j,u_j\rangle\). With the mathematics convention (linear in the first slot) the numerator becomes \(\langle v_k,u_j\rangle\); using the wrong one over \(\mathbb{C}\) inserts spurious conjugates and destroys orthogonality.

Derivation

\[ u_1=v_1,\qquad e_1=\frac{u_1}{\lVert u_1\rVert},\qquad \lVert u_1\rVert=\sqrt{\langle u_1,u_1\rangle} \]

Base case. Linear independence forces \(v_1\neq0\); positive-definiteness then gives \(\lVert u_1\rVert>0\), so the division is legal and \(\langle e_1,e_1\rangle=1\). A

\[ u_k=v_k-\sum_{j=1}^{k-1}\frac{\langle u_j,v_k\rangle}{\langle u_j,u_j\rangle}\,u_j \qquad(k\ge2) \]

Inductive step: subtract from \(v_k\) its component along each previously built orthogonal direction. Each coefficient is the orthogonal projection of \(v_k\) onto \(u_j\); the sum removes all overlap. B

\[ \langle u_i,u_k\rangle=\langle u_i,v_k\rangle-\sum_{j=1}^{k-1}\frac{\langle u_j,v_k\rangle}{\langle u_j,u_j\rangle}\,\langle u_i,u_j\rangle,\qquad i<k \]

Take the inner product of Step 2 with \(u_i\). Linearity in the second slot lets the finite sum and the scalar coefficients pass outside \(\langle u_i,\cdot\rangle\). B

\[ \langle u_i,u_j\rangle=\langle u_i,u_i\rangle\,\delta_{ij}\ (i,j<k)\ \Rightarrow\ \langle u_i,u_k\rangle=\langle u_i,v_k\rangle-\frac{\langle u_i,v_k\rangle}{\langle u_i,u_i\rangle}\langle u_i,u_i\rangle=0 \]

Inductive hypothesis: the earlier \(u_j\) are mutually orthogonal, so only the \(j=i\) term survives the sum. It cancels the first term exactly. Hence \(u_k\perp u_i\) for all \(i<k\); orthogonality propagates. C

\[ u_k=0\ \Rightarrow\ v_k=\sum_{j=1}^{k-1}c_j u_j\in\operatorname{span}\{v_1,\dots,v_{k-1}\}\ \ (\text{contradiction}) \]

Non-vanishing of \(u_k\). If \(u_k=0\), Step 2 writes \(v_k\) as a combination of \(u_1,\dots,u_{k-1}\), each of which lies in \(\operatorname{span}\{v_1,\dots,v_{k-1}\}\) — contradicting linear independence. So \(u_k\neq0\) and \(\lVert u_k\rVert>0\). C

\[ e_k=\frac{u_k}{\lVert u_k\rVert}\ \Rightarrow\ \langle e_i,e_k\rangle=\frac{\langle u_i,u_k\rangle}{\lVert u_i\rVert\,\lVert u_k\rVert}=\delta_{ik} \]

Normalise. Division by the positive number \(\lVert u_k\rVert\) (legal by Step 5) makes \(\langle e_k,e_k\rangle=1\); Step 4 makes the off-diagonal entries vanish. The set \(\{e_1,\dots,e_k\}\) is orthonormal. B

\[ \operatorname{span}\{e_1,\dots,e_k\}\subseteq\operatorname{span}\{v_1,\dots,v_k\},\quad \dim=k=\dim\ \Rightarrow\ \text{equal} \]

Span equality. Each \(e_j\) is a combination of \(v_1,\dots,v_j\) (unwind Steps 1–2), giving the inclusion. The \(e_j\) are orthonormal hence independent, so the left span has dimension \(k\); the \(v_j\) are independent, so the right span has dimension \(k\). A subspace of equal finite dimension inside another must coincide (dimension-invariance). C

Result

\[ u_k=v_k-\sum_{j<k}\frac{\langle u_j,v_k\rangle}{\langle u_j,u_j\rangle}\,u_j,\qquad e_k=\frac{u_k}{\lVert u_k\rVert},\qquad \langle e_i,e_j\rangle=\delta_{ij},\qquad \operatorname{span}\{e_{\le k}\}=\operatorname{span}\{v_{\le k}\} \]

Reading. Process the vectors in order; at each stage strip off everything the new vector already shares with the directions you have fixed, then scale what remains to unit length. What survives is orthogonal to all predecessors by construction, and nonzero precisely because the inputs were independent. The nested spans (the "flag") are preserved exactly — a triangular relationship between the \(v\)'s and the \(e\)'s that is the \(R\) of the \(QR\) decomposition.

Units check. If the \(v_k\) carry a physical dimension \([v]\), then \(u_k\) carries \([v]\) (a difference of like-dimensioned vectors, the projection coefficient \(\langle u_j,v_k\rangle/\langle u_j,u_j\rangle\) being dimensionless since numerator and denominator both scale as \([v]^2\)). Dividing \(u_k\) by \(\lVert u_k\rVert\) (also \([v]\)) makes each \(e_k\) dimensionless, as an orthonormal vector must be: \(\langle e_i,e_j\rangle=\delta_{ij}\) is a pure number.

Limiting cases

Already orthogonal inputs: every projection coefficient \(\langle u_j,v_k\rangle=0\), so \(u_k=v_k\) and the process only normalises — \(e_k=v_k/\lVert v_k\rVert\).
\(n=1\): the theorem reduces to normalisation, \(e_1=v_1/\lVert v_1\rVert\); no orthogonalisation is needed.
Two vectors in a plane: \(e_2\) is the component of \(v_2\) perpendicular to \(v_1\), normalised — the familiar "subtract the shadow" picture.
Monomials \(\{1,x,x^2,\dots\}\) on \([-1,1]\): Gram–Schmidt reproduces the Legendre polynomials up to normalisation.
Columns of an invertible matrix \(A\): the output assembles into \(A=QR\) with \(Q\) orthonormal and \(R\) upper-triangular with positive diagonal.

Breaks when

The inputs are linearly dependent. At the first dependent \(v_k\) the residual \(u_k\) is exactly zero, and \(e_k=u_k/\lVert u_k\rVert\) is \(0/0\). The algorithm must branch (skip the vector) rather than normalise — the plain formula fails at that step.
Indefinite or degenerate metric (pseudo-inner-product). On Minkowski space a nonzero vector can have \(\langle u,u\rangle=0\) (null) or \(<0\) (timelike). Then \(\lVert u_k\rVert=\sqrt{\langle u_k,u_k\rangle}\) vanishes or is imaginary, and orthonormalisation \(\langle e_i,e_j\rangle=\delta_{ij}\) is impossible; one gets a pseudo-orthonormal set \(\langle e_i,e_j\rangle=\eta_{ij}\) instead, and even that fails whenever a null direction is encountered.
Finite-precision arithmetic with near-dependent inputs. Classical Gram–Schmidt loses orthogonality catastrophically when the \(v_k\) are nearly parallel (small \(\lVert u_k\rVert\) amplifies rounding). The exact theorem holds; the naive computation does not, which is why modified Gram–Schmidt or Householder \(QR\) is used numerically.

Failure modes

Projecting onto the \(v_j\) instead of the \(u_j\). The subtraction must use the already-orthogonalised \(u_j\); projecting onto the original (non-orthogonal) \(v_j\) leaves cross-terms and the output is not orthogonal.
Wrong conjugation slot over \(\mathbb{C}\). Writing \(\langle v_k,u_j\rangle\) when the convention wants \(\langle u_j,v_k\rangle\) conjugates the coefficient and breaks \(\langle u_i,u_k\rangle=0\). The two differ by a complex conjugate.
Normalising too early. If you replace \(u_j\) by \(e_j\) then the correct coefficient is \(\langle e_j,v_k\rangle\) with denominator \(1\); dividing by \(\langle u_j,u_j\rangle\) while projecting onto the already-normalised \(e_j\) double-counts the norm.
Forgetting a term in the sum. Subtracting only the last projection instead of all \(j<k\) leaves \(u_k\) non-orthogonal to the earlier directions.
Assuming the output is unique. Each \(e_k\) is fixed only up to a unit phase \(e^{i\theta}\) (a sign over \(\mathbb{R}\)); pinning the phase requires an extra convention (e.g. positive diagonal of \(R\)).

Discussion

The heart of the proof is a single induction that does double duty: it simultaneously establishes orthogonality (Step 4) and non-vanishing (Step 5). Orthogonality is pure linear algebra — expand and cancel — but non-vanishing is precisely where linear independence is spent. This is the honest accounting the result demands: without independence you still get an orthogonal set, you just cannot guarantee it is a basis of the right size, because some directions collapse to zero.

Geometrically, \(u_k\) is the component of \(v_k\) orthogonal to the subspace \(W_{k-1}=\operatorname{span}\{v_1,\dots,v_{k-1}\}\); the summed projection in Step 2 is the orthogonal projection \(P_{W_{k-1}}v_k\), and \(u_k=(I-P_{W_{k-1}})v_k\). Cauchy–Schwarz is what guarantees this projection is well-behaved: \(|\langle u_j,v_k\rangle|\le\lVert u_j\rVert\,\lVert v_k\rVert\) bounds each coefficient, so no term blows up, and it is the equality condition of Cauchy–Schwarz (parallel vectors) that marks exactly the dependent case where \(u_k\to0\).

The triangular structure is not incidental — it is the content of the \(QR\) decomposition. Collecting \(v_k=\sum_{j\le k}r_{jk}e_j\) with \(r_{jk}=\langle e_j,v_k\rangle\) and \(r_{kk}=\lVert u_k\rVert>0\) gives \(A=QR\): the algorithm is a constructive proof that every full-rank matrix factors into an isometry times an upper-triangular matrix with positive diagonal. That factorisation is the basis of the \(QR\) eigenvalue algorithm, least-squares solving, and the numerical construction of orthogonal polynomials.

In infinite dimensions the same recursion runs verbatim on any countable independent set and, applied to \(\{1,x,x^2,\dots\}\) in \(L^2[a,b]\) with a weight \(w(x)\), generates the classical orthogonal polynomials; completeness of the resulting system (that the \(e_k\) span a dense subspace, not merely each finite flag) is a separate analytic theorem, not part of Gram–Schmidt itself. The abstract statement also has a coordinate-free reading: Gram–Schmidt is the Iwasawa/\(QR\) factorisation \(GL(n)=O(n)\cdot T_+\) of the general linear group into orthogonal and upper-triangular-positive parts, and the uniqueness-up-to-phase is exactly the intersection \(O(n)\cap T_+=\{\pm1\}^n\) being discrete.

Common misconceptions. Gram–Schmidt does not require finite dimension (only a countable independent set and a genuine inner product), and it does not by itself deliver a basis of the whole space — only of the span of the inputs. It is also not the same as "diagonalising the Gram matrix": that gives an orthonormal set but scrambles the nested spans, whereas Gram–Schmidt is the unique-up-to-phase procedure that keeps every partial flag intact.

Worked examples

\(\mathbb{R}^3\) with the dot product. Orthonormalise \(v_1=(1,1,0),\ v_2=(1,0,1),\ v_3=(0,1,1)\).

Setup: standard inner product \(\langle a,b\rangle=\sum_i a_ib_i\), process in the given order. A

\[ u_1=(1,1,0),\quad \lVert u_1\rVert=\sqrt{2},\quad e_1=\tfrac{1}{\sqrt2}(1,1,0) \]

Base case: normalise the first vector. A

\[ \frac{\langle u_1,v_2\rangle}{\langle u_1,u_1\rangle}=\frac{1}{2},\quad u_2=(1,0,1)-\tfrac12(1,1,0)=\left(\tfrac12,-\tfrac12,1\right),\quad e_2=\tfrac{1}{\sqrt6}(1,-1,2) \]

Subtract the projection onto \(u_1\); \(\lVert u_2\rVert^2=\tfrac14+\tfrac14+1=\tfrac32\), so scaling \((\tfrac12,-\tfrac12,1)\) to unit length gives \((1,-1,2)/\sqrt6\). B

\[ \frac{\langle u_1,v_3\rangle}{\langle u_1,u_1\rangle}=\tfrac12,\quad \frac{\langle u_2,v_3\rangle}{\langle u_2,u_2\rangle}=\frac{1/2}{3/2}=\tfrac13,\quad u_3=(0,1,1)-\tfrac12(1,1,0)-\tfrac13\left(\tfrac12,-\tfrac12,1\right)=\left(-\tfrac23,\tfrac23,\tfrac23\right) \]

Subtract projections onto both earlier orthogonal directions. \(\langle u_1,v_3\rangle=1,\ \langle u_2,v_3\rangle=-\tfrac12+1=\tfrac12\). B

\[ \lVert u_3\rVert^2=3\cdot\tfrac49=\tfrac43,\quad e_3=\frac{(-\tfrac23,\tfrac23,\tfrac23)}{2/\sqrt3}=\tfrac{1}{\sqrt3}(-1,1,1) \]

Normalise. Check: \(e_1\!\cdot\!e_2=\tfrac{1-1+0}{\sqrt{12}}=0\), \(e_1\!\cdot\!e_3=\tfrac{-1+1+0}{\sqrt6}=0\), \(e_2\!\cdot\!e_3=\tfrac{-1-1+2}{\sqrt{18}}=0\). C

\[ e_1=\tfrac{1}{\sqrt2}(1,1,0),\quad e_2=\tfrac{1}{\sqrt6}(1,-1,2),\quad e_3=\tfrac{1}{\sqrt3}(-1,1,1) \]

Reading. Three independent but non-orthogonal vectors are replaced by a mutually perpendicular unit triad spanning all of \(\mathbb{R}^3\), with each partial span preserved: \(\operatorname{span}\{e_1\}=\operatorname{span}\{v_1\}\), \(\operatorname{span}\{e_1,e_2\}=\operatorname{span}\{v_1,v_2\}\).

Units check. Inputs and outputs are pure-number vectors; all pairwise inner products come out as \(\delta_{ij}\), dimensionless.

Function space \(L^2[-1,1]\) with \(\langle f,g\rangle=\int_{-1}^{1}f(x)g(x)\,dx\). Orthonormalise the monomials \(v_1=1,\ v_2=x,\ v_3=x^2\).

Setup: same recursion, inner product now an integral (real, so no conjugation). B

\[ u_1=1,\quad \langle u_1,u_1\rangle=\int_{-1}^{1}dx=2,\quad e_1=\tfrac{1}{\sqrt2} \]

Base case: normalise the constant. A

\[ \langle u_1,v_2\rangle=\int_{-1}^{1}x\,dx=0\ \Rightarrow\ u_2=x,\quad \langle u_2,u_2\rangle=\int_{-1}^{1}x^2dx=\tfrac23,\quad e_2=\sqrt{\tfrac32}\,x \]

\(x\) is already orthogonal to the constant (odd function integrates to zero), so only normalisation is needed. B

\[ \langle u_1,v_3\rangle=\int_{-1}^{1}x^2dx=\tfrac23,\quad \langle u_2,v_3\rangle=\int_{-1}^{1}x^3dx=0,\quad u_3=x^2-\frac{2/3}{2}\cdot1-0=x^2-\tfrac13 \]

Subtract the projection onto \(u_1\); the projection onto \(u_2\) vanishes by parity. B

\[ \langle u_3,u_3\rangle=\int_{-1}^{1}\!\left(x^2-\tfrac13\right)^2dx=\tfrac25-\tfrac29=\tfrac{8}{45},\quad e_3=\sqrt{\tfrac{45}{8}}\left(x^2-\tfrac13\right)=\tfrac{\sqrt{10}}{4}\,(3x^2-1) \]

Normalise: \(\int_{-1}^{1}(x^4-\tfrac23x^2+\tfrac19)dx=\tfrac25-\tfrac23\cdot\tfrac23+\tfrac19\cdot2=\tfrac25-\tfrac29=\tfrac{8}{45}\). C

\[ e_1=\tfrac{1}{\sqrt2},\quad e_2=\sqrt{\tfrac32}\,x,\quad e_3=\tfrac{\sqrt{10}}{4}\,(3x^2-1) \]

Reading. Gram–Schmidt on \(\{1,x,x^2\}\) reproduces the first three Legendre polynomials \(P_0=1,\ P_1=x,\ P_2=\tfrac12(3x^2-1)\), each rescaled to unit \(L^2\) norm. This is the concrete link to the waves thread: orthogonal polynomials are Gram–Schmidt applied to monomials under a weighted integral inner product.

Units check. With \(x\) dimensionless, each \(e_k\) has dimension \([f]=[\text{length}]^{-1/2}\) inherited from \(\int f^2\,dx=1\); the normalisation constants carry it consistently.

Problems

In \(\mathbb{R}^2\) orthonormalise \(v_1=(3,0),\ v_2=(1,2)\) (dot product).
Solution
\(u_1=(3,0),\ \lVert u_1\rVert=3,\ e_1=(1,0)\). Coefficient \(\langle u_1,v_2\rangle/\langle u_1,u_1\rangle=3/9=\tfrac13\), so \(u_2=(1,2)-\tfrac13(3,0)=(0,2)\), \(\lVert u_2\rVert=2\), \(e_2=(0,1)\). Result \(e_1=(1,0),\ e_2=(0,1)\) — these inputs happen to determine the standard basis.
Orthonormalise \(v_1=(1,1,1),\ v_2=(1,1,0)\) in \(\mathbb{R}^3\), in that order.
Solution
\(u_1=(1,1,1),\ \langle u_1,u_1\rangle=3,\ e_1=\tfrac{1}{\sqrt3}(1,1,1)\). Coefficient \(\langle u_1,v_2\rangle/3=2/3\), so \(u_2=(1,1,0)-\tfrac23(1,1,1)=(\tfrac13,\tfrac13,-\tfrac23)\), \(\lVert u_2\rVert^2=\tfrac19+\tfrac19+\tfrac49=\tfrac23\), \(e_2=\tfrac{1}{\sqrt6}(1,1,-2)\). Check \(e_1\!\cdot\!e_2=\tfrac{1+1-2}{\sqrt{18}}=0\).
In \(\mathbb{C}^2\) with \(\langle a,b\rangle=\bar a_1 b_1+\bar a_2 b_2\), orthonormalise \(v_1=(1,i),\ v_2=(1,0)\).
Solution
\(u_1=(1,i),\ \langle u_1,u_1\rangle=1+1=2,\ e_1=\tfrac{1}{\sqrt2}(1,i)\). Coefficient \(\langle u_1,v_2\rangle/2\) with \(\langle u_1,v_2\rangle=\bar1\cdot1+\bar i\cdot0=1\), so \(u_2=(1,0)-\tfrac12(1,i)=(\tfrac12,-\tfrac{i}{2})\), \(\langle u_2,u_2\rangle=\tfrac14+\tfrac14=\tfrac12\), \(e_2=\tfrac{1}{\sqrt2}(1,-i)\). Check \(\langle e_1,e_2\rangle=\tfrac12(\bar1\cdot1+\bar i\cdot(-i))=\tfrac12(1+(-i)(-i))=\tfrac12(1-1)=0\); the conjugate in the first slot is essential.
Show that applying Gram–Schmidt to the columns of a real invertible \(n\times n\) matrix \(A\) yields a factorisation \(A=QR\) with \(Q\) orthogonal and \(R\) upper-triangular with positive diagonal.
Solution
Let \(v_1,\dots,v_n\) be the columns (independent since \(A\) is invertible) and \(e_1,\dots,e_n\) the Gram–Schmidt output; set \(Q=[e_1\,\cdots\,e_n]\), which is orthogonal because \(Q^{\!\top}Q=[\langle e_i,e_j\rangle]=I\). By the span property each \(v_k\in\operatorname{span}\{e_1,\dots,e_k\}\), so \(v_k=\sum_{j\le k}r_{jk}e_j\) with \(r_{jk}=\langle e_j,v_k\rangle\) and \(r_{jk}=0\) for \(j>k\) — i.e. \(R=[r_{jk}]\) is upper-triangular. Its diagonal is \(r_{kk}=\langle e_k,v_k\rangle=\lVert u_k\rVert>0\) (nonzero by independence, Step 5). Stacking the columns gives \(A=QR\).
In \(L^2[0,1]\) with \(\langle f,g\rangle=\int_0^1 fg\,dx\), orthonormalise \(v_1=1,\ v_2=x\).
Solution
\(u_1=1,\ \langle u_1,u_1\rangle=\int_0^1 dx=1,\ e_1=1\). Coefficient \(\langle u_1,v_2\rangle=\int_0^1 x\,dx=\tfrac12\), so \(u_2=x-\tfrac12\). Norm \(\int_0^1(x-\tfrac12)^2dx=\int_0^1(x^2-x+\tfrac14)dx=\tfrac13-\tfrac12+\tfrac14=\tfrac1{12}\), so \(e_2=\sqrt{12}\,(x-\tfrac12)=\sqrt3\,(2x-1)\). These are the first two shifted (normalised) Legendre polynomials on \([0,1]\); check \(\int_0^1 1\cdot\sqrt3(2x-1)\,dx=\sqrt3\int_0^1(2x-1)dx=0\).