BrevFeed
Sign in

28. Gram-Schmidt orthogonalization process

📖 Reading · 10 min
💡 Every code box below is live — edit it and hit Run.

§20.2 showed an orthogonal set is automatically independent — a one-way street. This lesson builds the road back the other direction: Gram-Schmidt takes any independent set and mechanically produces an orthogonal set with the same span, one vector at a time.

The projection formula

The building block is the projection of \vec u onto \vec v:

\operatorname{proj}_{\vec v}\vec u=\frac{\vec u\cdot\vec v}{\vec v\cdot\vec v}\vec v

— the "shadow" \vec u casts on the line through \vec v: a scalar multiple of \vec v (so it points along \vec v's direction), chosen so that \vec u-\operatorname{proj}_{\vec v}\vec u is orthogonal to \vec v. Verify the orthogonality directly:

(\vec u-\operatorname{proj}_{\vec v}\vec u)\cdot\vec v=\vec u\cdot\vec v-\frac{\vec u\cdot\vec v}{\vec v\cdot\vec v}(\vec v\cdot\vec v)=\vec u\cdot\vec v-\vec u\cdot\vec v=0

— this cancellation is the entire mechanism Gram-Schmidt runs repeatedly: subtract off the part of a vector that points along a direction already handled, and what's left is guaranteed orthogonal to that direction.

The algorithm

Given independent \vec v_1,\dots,\vec v_k, build orthogonal \vec u_1,\dots,\vec u_k with the same span at every stage:

\vec u_1=\vec v_1 \vec u_2=\vec v_2-\operatorname{proj}_{\vec u_1}\vec v_2 \vec u_3=\vec v_3-\operatorname{proj}_{\vec u_1}\vec v_3-\operatorname{proj}_{\vec u_2}\vec v_3 \vdots \vec u_k=\vec v_k-\sum_{i=1}^{k-1}\operatorname{proj}_{\vec u_i}\vec v_k

Each \vec u_j is orthogonal to every \vec u_i with i<j: it's \vec v_j with the projections onto all previous \vec u_i's removed, so by the cancellation identity above (applied to each \vec u_i in turn — and using that \vec u_i\cdot\vec u_m=0 for i\neq m<j, an inductive fact), nothing pointing along any earlier direction survives. \vec u_j\neq\vec0 for every j: if it were zero, \vec v_j would be a linear combination of \vec u_1,\dots,\vec u_{j-1} (hence of \vec v_1,\dots,\vec v_{j-1}, since each \vec u_i is built from earlier \vec v's), contradicting the original set's independence.

\operatorname{span}\{\vec u_1,\dots,\vec u_j\}=\operatorname{span}\{\vec v_1,\dots,\vec v_j\} at every stage — each \vec u_j is \vec v_j minus a combination of earlier \vec u_i's (hence of earlier \vec v_i's), so the spans never change, only the basis describing them.

Doing it in Python

import numpy as np

def gram_schmidt(vectors):
    us = []
    for v in vectors:
        w = v.astype(float).copy()
        for u in us:
            w = w - (v @ u) / (u @ u) * u
        us.append(w)
    return us

v1 = np.array([1., 1., 0.])
v2 = np.array([1., 0., 1.])
v3 = np.array([0., 1., 1.])

us = gram_schmidt([v1, v2, v3])
for i, u in enumerate(us, 1):
    print(f"u{i} =", [round(x, 4) for x in u.tolist()])

print("\npairwise dot products:")
print("u1.u2 =", round(float(us[0] @ us[1]), 10))
print("u1.u3 =", round(float(us[0] @ us[2]), 10))
print("u2.u3 =", round(float(us[1] @ us[2]), 10))

u1 = [1.0, 1.0, 0.0]
u2 = [0.5, -0.5, 1.0]
u3 = [-0.6667, 0.6667, 0.6667]

pairwise dot products:
u1.u2 = 0.0
u1.u3 = 0.0
u2.u3 = 0.0

Confirming the span is preserved at every stage:

import numpy as np

v1 = np.array([1., 1., 0.])
v2 = np.array([1., 0., 1.])

u1 = v1.copy()
u2 = v2 - (v2 @ u1) / (u1 @ u1) * u1

# span{u1, u2} == span{v1, v2}: same rank when stacked together
combined_v = np.column_stack([v1, v2])
combined_u = np.column_stack([u1, u2])
combined_both = np.column_stack([v1, v2, u1, u2])

print("rank{v1,v2} =", np.linalg.matrix_rank(combined_v))
print("rank{u1,u2} =", np.linalg.matrix_rank(combined_u))
print("rank{v1,v2,u1,u2} =", np.linalg.matrix_rank(combined_both),
      " (same as either alone -> identical span)")

rank{v1,v2} = 2
rank{u1,u2} = 2
rank{v1,v2,u1,u2} = 2  (same as either alone -> identical span)

Worked example

Apply Gram-Schmidt to \vec v_1=(1,0), \vec v_2=(2,1).

\vec u_1=\vec v_1=(1,0).

\operatorname{proj}_{\vec u_1}\vec v_2=\frac{(2,1)\cdot(1,0)}{(1,0)\cdot(1,0)}(1,0)=\frac21(1,0)=(2,0)

\vec u_2=(2,1)-(2,0)=(0,1)

\boxed{\vec u_1=(1,0),\quad\vec u_2=(0,1)}

Sanity check. \vec u_1\cdot\vec u_2=1(0)+0(1)=0 ✓ orthogonal. Span check: \{(1,0),(2,1)\} spans \mathbb{R}^2 (independent, since neither is a multiple of the other), and so does \{(1,0),(0,1)\} (the standard basis) — same span, different (now orthogonal) description, exactly as guaranteed.

Your turn

1. Apply Gram-Schmidt to \vec v_1=(1,1), \vec v_2=(0,1).

2. Why must \vec u_1=\vec v_1 exactly, rather than some projection of it?

3. True or false: Gram-Schmidt applied to an already-orthogonal set returns the same vectors unchanged.

Solutions

1. \vec u_1=(1,1). \operatorname{proj}_{\vec u_1}\vec v_2=\frac{(0,1)\cdot(1,1)}{(1,1)\cdot(1,1)}(1,1)=\frac12(1,1)=(0.5,0.5). \vec u_2=(0,1)-(0.5,0.5)=(-0.5,0.5). Check: (1,1)\cdot(-0.5,0.5)=-0.5+0.5=0 ✓.

2. There's nothing to project it onto yet. \vec u_1 is the very first vector in the process — projection removes the component along previously built orthogonal vectors, and there are none before the first step, so \vec v_1 passes through completely unchanged. (This is also why the algorithm's output depends on the order the input vectors are given in: starting from a different \vec v_1 generally produces a different — though equally valid — orthogonal set spanning the same space.)

3. True. If \vec v_j is already orthogonal to \vec v_1,\dots, \vec v_{j-1} (hence to \vec u_1,\dots,\vec u_{j-1}, which span the same space), every projection term \operatorname{proj}_{\vec u_i}\vec v_j=\frac{\vec v_j\cdot\vec u_i}{\vec u_i\cdot\vec u_i}\vec u_i has numerator \vec v_j\cdot\vec u_i=0, so it contributes nothing — \vec u_j=\vec v_j at every step. Gram-Schmidt is idempotent on already-orthogonal input, a good sanity check to run on the algorithm itself.

Check yourself in code

Apply Gram-Schmidt to \vec v_1=(1,1,1), \vec v_2=(1,0,2), rounding each result to 4 decimal places, and confirm orthogonality.

Print exactly this:

u1 = [1.0, 1.0, 1.0]
u2 = [0.0, -1.0, 1.0]
u1 . u2 = 0.0

import numpy as np

v1 = np.array([1., 1., 1.])
v2 = np.array([1., 0., 2.])

u1 = v1.copy()
u2 = v2 - (v2 @ u1) / (u1 @ u1) * u1

print("u1 =", [round(x, 4) for x in u1.tolist()])
# print u2 the same way, then u1 . u2

import numpy as np

v1 = np.array([1., 1., 1.])
v2 = np.array([1., 0., 2.])

u1 = v1.copy()
u2 = v2 - (v2 @ u1) / (u1 @ u1) * u1

print("u1 =", [round(x, 4) for x in u1.tolist()])
print("u2 =", [round(x, 4) for x in u2.tolist()])
print("u1 . u2 =", round(float(u1 @ u2), 10))

Gram-Schmidt turns any independent set into an orthogonal one with identical span, one vector at a time: subtract off every projection onto vectors already built, and what remains is guaranteed orthogonal to all of them, and guaranteed nonzero by independence.

Next: normalizing a Gram-Schmidt output into orthonormal vectors — orthogonal and unit length — the single most convenient kind of basis this course ever works with.

← Previous Next →