Gram–Schmidt Builds Orthogonal Directions

Two feature directions can look different while carrying almost the same information. Normalizing each one fixes its length, but it does not remove the overlap. Gram–Schmidt removes the part of each new vector already explained by earlier directions, then normalizes what remains. The length of that remainder becomes a useful diagnostic: zero means dependence, while a very small value warns that the new direction is nearly redundant or numerically fragile.

By the end of this lesson, you will be able to run Gram–Schmidt on two or three vectors, verify the orthonormal output, and interpret residual size as a scale-dependent signal of overlap. You will also know why the order of the inputs and the chosen floating-point tolerance affect the result.

Learning outcome: Execute Gram–Schmidt on two or three vectors and identify instability.

Lesson 046: Orthonormal Bases Simplify Coordinates showed why inner products read coordinates cleanly when the directions are orthonormal. This lesson explains how to construct those directions from a linearly independent list. The official Mathematics for Machine Learning book maps the idea to Section 3.5. The worked vectors, visual, NumPy implementation, AI connection, and assessments here are original.

Each new direction keeps only its unexplained part

Let $a_1,\ldots,a_k$ be vectors in $\mathbb{R}^d$. The input vectors need not be perpendicular or unit length. Gram–Schmidt constructs unit vectors $q_1,\ldots,q_k$ in the same span, provided that each new input contributes a nonzero direction after earlier directions have been removed.

The first vector has nothing to subtract. For later vectors, use the inner product with each earlier unit direction to measure how much of the new vector points along it. Because each $q_j$ has unit length, $\langle q_j,a_i\rangle q_j$ is exactly the projection of $a_i$ onto that direction.

The word residual is doing important work. It is not merely an intermediate vector to discard after normalization. Its norm $\lVert r_i\rVert_2$ is the distance from $a_i$ to the span of the earlier directions. A large residual means that the input contributes a substantial new component relative to the chosen order. A small residual means that the input is close to the earlier span. A zero residual cannot be normalized at all.

The construction is sequential, so the order is part of the calculation. If the inputs are independent, different orders still describe the same span in exact arithmetic, but they can produce different directions, signs, and residual sizes at each step. That matters when the residual is close to the precision or noise floor.

A three-vector calculation, one projection at a time

Use three vectors in $\mathbb{R}^3$:

$$a_1=\begin{bmatrix}1\\1\\0\end{bmatrix}, \qquad a_2=\begin{bmatrix}1\\0\\1\end{bmatrix}, \qquad a_3=\begin{bmatrix}0\\1\\1\end{bmatrix}.$$

The coordinates have no physical unit here; they are three components of a toy feature direction. Keeping the numbers small makes the projection and normalization visible.

The visual below turns that sequence into a reader-controlled walk. Each panel shows the input, the total projection removed, the residual, and the unit direction for one stage. The table is not a decorative duplicate: it is the server-rendered numerical fallback and keeps the complete comparison available to keyboard and assistive-technology users.

Follow the residual through Gram–Schmidt

Use the buttons or the select menu. The selected stage is announced, and the table below remains available without interaction.

Choose a stagestart with a₁remove q₁ from a₂remove q₁ and q₂ from a₃

Previous stage1 of 3: start with a₁Next stage

Input vector

(1, 1, 0)

Projection removed

(0, 0, 0)

Residual before normalization

(1, 1, 0)

Normalized direction

(0.71, 0.71, 0)

All Gram–Schmidt stages and their three-dimensional vectors.

Stage	Input aᵢ	Removed projection	Residual rᵢ	Unit direction qᵢ
start with a₁	(1, 1, 0)	(0, 0, 0)	(1, 1, 0)	(0.71, 0.71, 0)
remove q₁ from a₂	(1, 0, 1)	(0.5, 0.5, 0)	(0.5, −0.5, 1)	(0.41, −0.41, 0.82)
remove q₁ and q₂ from a₃	(0, 1, 1)	(0.67, 0.33, 0.33)	(−0.67, 0.67, 0.67)	(−0.58, 0.58, 0.58)

The diagram rounds values to two decimal places; the worked calculation in the lesson keeps exact radicals and fractions.

NumPy mirrors the derivation and exposes the checks

The implementation can keep directions as rows. For each input, it starts with a copy of the original vector, subtracts one projection at a time, measures the residual norm, and normalizes only when that norm clears the tolerance. numpy.linalg.norm provides the Euclidean lengths used here. This one-at-a-time form is often called modified Gram–Schmidt; in exact arithmetic it gives the same subtractions as the displayed sum, while the update order makes the current residual explicit.

import numpy as np


def gram_schmidt(vectors, rtol=1e-12):
    directions = []
    residual_norms = []

    for index, vector in enumerate(vectors, start=1):
        original = np.asarray(vector, dtype=float)
        if original.ndim != 1:
            raise ValueError(f"vector {index} must be one-dimensional")

        residual = original.copy()
        for q in directions:
            coefficient = q @ residual
            residual -= coefficient * q

        residual_norm = np.linalg.norm(residual)
        residual_norms.append(residual_norm)
        scale = max(1.0, np.linalg.norm(original))
        if residual_norm <= rtol * scale:
            raise ValueError(
                f"vector {index} is dependent at rtol={rtol:g}"
            )

        directions.append(residual / residual_norm)

    return np.array(directions), np.array(residual_norms)


A = np.array([
    [1.0, 1.0, 0.0],
    [1.0, 0.0, 1.0],
    [0.0, 1.0, 1.0],
])

Q, residual_norms = gram_schmidt(A)
np.set_printoptions(precision=8, suppress=True)
print(Q)
print(residual_norms)
print(Q @ Q.T)
print(np.allclose(Q @ Q.T, np.eye(3), rtol=1e-12, atol=1e-12))

The verified output is:

[[ 0.70710678  0.70710678  0.        ]
 [ 0.40824829 -0.40824829  0.81649658]
 [-0.57735027  0.57735027  0.57735027]]
[1.41421356 1.22474487 1.15470054]
[[ 1.  0.  0.]
 [ 0.  1.  0.]
 [ 0.  0.  1.]]
True

Because the directions are stored as rows, Q @ Q.T is the Gram matrix of the output. The identity matrix confirms unit lengths on the diagonal and zero cross-inner-products elsewhere. If an implementation stores directions as columns instead, the corresponding check is Q.T @ Q; the shape convention must be decided before the multiplication is written.

np.allclose compares floating-point values with explicit relative and absolute tolerances rather than demanding exact decimal equality. That is a test of the computed result, not a license to treat every small number as zero. The tolerance has to match the scale and noise of the feature data.

Residual size is a useful overlap warning, not a verdict

Consider a nearly duplicated pair:

$$a_1=(1,1,0), \qquad a_2=(1,1+10^{-10},0).$$

After removing the $q_1$ component, the second residual is approximately

$$r_2\approx(-5\times10^{-11},\,5\times10^{-11},\,0), \qquad \lVert r_2\rVert_2\approx7.07\times10^{-11}.$$

With rtol=1e-12, the residual clears the threshold for this input scale and the code accepts a very small new direction. With rtol=1e-9, it is flagged as dependent. Neither setting discovers a universal truth about the features: the first is asking for finer separation, while the second treats that separation as below the application’s useful resolution.

An exactly dependent vector is clearer. If $a_2=2a_1$, subtracting its projection leaves $r_2=0$. Division by $\lVert r_2\rVert_2$ would be division by zero, so a robust implementation must stop, drop the vector, or report the redundancy. It must not emit nan and continue as if a new direction existed.

The scale-aware check in the code uses $\lVert r_i\rVert_2\le\texttt{rtol}\,\max(1,\lVert a_i\rVert_2)$ as its decision rule. In a real feature pipeline, choose the threshold from expected measurement noise, feature scaling, and the downstream decision. Record that choice with the diagnostic. A tolerance without a scale is hard to interpret.

Order changes the basis and the stability story

Gram–Schmidt is a span-preserving construction, not an order-free sorting operation. Reordering the independent inputs changes which direction is built first, so it changes the projections and the residual norms that follow. The final lists still span the same subspace in exact arithmetic, but the basis vectors can be rotated within that subspace or negated. Near dependence makes the difference more consequential: putting a nearly redundant vector early can make later residuals look small, while a better-conditioned direction first can leave more numerical headroom.

For a small teaching calculation, the explicit loop is valuable because every projection is visible. For a larger or ill-conditioned matrix, use a tested QR factorization such as numpy.linalg.qr and inspect its rank or reconstruction diagnostics. That is a change of numerical tool, not a change in the geometric goal.

Feature directions are not statistically independent variables

Suppose rows of a matrix are feature directions: perhaps hand-designed measurement patterns, learned probes, or directions used to inspect a model’s representation. Gram–Schmidt can produce an orthogonal coordinate frame for their span. The residual norm tells the practitioner how much new linear direction each row adds relative to the rows processed earlier. This can guide feature review or reveal a near-duplicate before a downstream linear analysis.

It does not prove that the original random variables are statistically independent, and it does not promise better predictions. Orthogonality compares fixed vectors under a chosen inner product. Statistical independence is a statement about a joint probability distribution. Even zero covariance only rules out a particular linear relationship.

For a concrete counterexample, take centered samples

$$x=(-2,-1,0,1,2), \qquad y=(2,-1,-2,-1,2)=x^2-2.$$

Their dot product is

$$x^{\mathsf T}y=(-2)(2)+(-1)(-1)+0(-2)+1(-1)+2(2)=0.$$

The sample columns have zero linear overlap, yet $y$ is completely determined by $x$ through a nonlinear rule. In an AI report, say that directions are orthogonal or that centered sample covariance is zero under the stated convention. Reserve “independent” for a conclusion supported by a separate probabilistic argument.

Gram–Schmidt therefore gives a practical bridge from the orthonormal-coordinate idea in the previous lesson to feature diagnostics: remove the explained part, measure what is left, normalize only when the remainder is meaningful, and report the order and tolerance that made the decision. The next lesson will use the same residual idea to describe directions orthogonal to an entire subspace.

Check your understanding

Question 1

For input vectors a₁,…,aᵢ and existing unit directions q₁,…,qᵢ₋₁, write the residual rᵢ and explain why its norm is informative.

Show the step-by-step solution

Start by subtracting the projection onto each earlier direction:

$$r_i=a_i-\sum_{j=1}^{i-1}\langle q_j,a_i\rangle q_j.$$

The coefficient $\langle q_j,a_i\rangle$ measures the component along $q_j$ because $q_j$ has unit length. After all those components are removed, $\lVert r_i\rVert_2$ is the distance from $a_i$ to the span of the earlier directions. A zero value means no new direction remains; a small value means the input is nearly in that span at the chosen numerical scale.

Question 2

Run the first two Gram–Schmidt steps for a₁=(1, 1) and a₂=(1, 0). Give q₁, the projection removed from a₂, r₂, and q₂.

Show the step-by-step solution

The first norm is

$$\lVert a_1\rVert_2=\sqrt{1^2+1^2}=\sqrt2, \qquad q_1=\begin{bmatrix}1/\sqrt2\\1/\sqrt2\end{bmatrix}.$$

The projection coefficient for $a_2$ is

$$q_1^{\mathsf T}a_2=\frac1{\sqrt2}.$$

Therefore the projection vector is $(1/2,1/2)$ and

$$r_2=\begin{bmatrix}1\\0\end{bmatrix} -\begin{bmatrix}1/2\\1/2\end{bmatrix} =\begin{bmatrix}1/2\\-1/2\end{bmatrix}.$$

Its norm is $1/\sqrt2$, so normalizing gives

$$q_2=\begin{bmatrix}1/\sqrt2\\-1/\sqrt2\end{bmatrix}.$$

The cross-inner-product is $q_1^{\mathsf T}q_2=1/2-1/2=0$, which checks the new direction.

Question 3

In the lesson’s three-vector example, why is the total projection removed from a₃ equal to (2/3, 1/3, 1/3)? Show both projection vectors before adding them.

Show the step-by-step solution

The two coefficients are

$$q_1^{\mathsf T}a_3=\frac1{\sqrt2}, \qquad q_2^{\mathsf T}a_3=\frac1{\sqrt6}.$$

Multiply each coefficient by its direction:

$$\frac1{\sqrt2}q_1 =\begin{bmatrix}1/2\\1/2\\0\end{bmatrix}, \qquad \frac1{\sqrt6}q_2 =\begin{bmatrix}1/6\\-1/6\\1/3\end{bmatrix}.$$

Adding corresponding coordinates gives

$$\begin{bmatrix}1/2+1/6\\1/2-1/6\\0+1/3\end{bmatrix} =\begin{bmatrix}2/3\\1/3\\1/3\end{bmatrix}.$$

That vector is the part of $a_3$ already explained by the first two directions; subtracting it leaves the residual used to build $q_3$.

Question 4

What should a Gram–Schmidt implementation do when a new input is exactly a multiple of an earlier input, such as a₂=2a₁?

Show the step-by-step solution

If $a_2=2a_1$, the first direction $q_1$ spans $a_2$. The projection is

$$(q_1^{\mathsf T}a_2)q_1=2\lVert a_1\rVert_2q_1=2a_1=a_2,$$

so the residual is $r_2=0$. Its norm is zero, and the expression $r_2/\lVert r_2\rVert_2$ is undefined. The implementation must report the dependence, drop the duplicate, or apply an explicitly documented policy; it must not normalize the zero vector or silently keep nan values.

Question 5

Why can two different input orders produce different q-vectors even when they span the same subspace? What remains invariant in exact arithmetic?

Show the step-by-step solution

Each new direction is made orthogonal to the directions that came before it. Changing the order changes that earlier span at each step, so the projection, residual, and normalized direction can all change. A direction can also flip sign without changing the one-dimensional subspace it represents.

For linearly independent inputs, exact arithmetic still preserves the span of the complete list: every original vector can be reconstructed from the output directions and the output directions lie in the original span. The individual basis vectors and the intermediate residual sizes are order-dependent, which is why order is part of a numerical diagnostic.

Question 6

The code stores q-vectors as rows. Which matrix check verifies orthonormality, and what row-versus-column mistake would change the check?

Show the step-by-step solution

With rows, the $(i,j)$ entry of $Q Q^{\mathsf T}$ is the inner product of row $i$ and row $j$. Therefore the correct check is

$$Q Q^{\mathsf T}\approx I.$$

In NumPy, use np.allclose(Q @ Q.T, np.eye(Q.shape[0]), atol=...) with a stated tolerance. If the directions are stored as columns, the corresponding Gram matrix is $Q^{\mathsf T}Q$ instead. Using the column check on row-stored data can inspect the wrong pairwise relationships or even produce a shape that does not match the intended identity.

Question 7

A new feature direction has residual norm 0.00000000007 after earlier directions are removed. Does that automatically mean the feature is useless?

Show the step-by-step solution

No. The value says that the new direction contributes very little relative to the earlier span at the current input scale and ordering. Whether that is negligible depends on measurement noise, feature units, downstream sensitivity, and the chosen tolerance.

For example, rtol=1e-12 may accept a residual near $7\times10^{-11}$ for a unit-scale vector, while rtol=1e-9 may classify it as dependent. The correct interpretation is a documented numerical decision—“below this resolution” or “retained at this resolution”—not a universal claim about predictive value.

Question 8

Why does zero centered dot product between feature samples fail to prove statistical independence? Use x=(−2,−1,0,1,2) and y=(2,−1,−2,−1,2).

Show the step-by-step solution

First compute the linear overlap:

$$x^{\mathsf T}y=-4+1+0-1+4=0.$$

Because both sample vectors are centered, this also gives zero empirical covariance under the usual equal-weight convention. But the relation $y=x^2-2$ determines every value of $y$ from $x$. The pair therefore has a nonlinear dependence even though its linear inner product is zero.

Orthogonality and zero covariance describe one linear comparison. Independence requires a statement about the full joint distribution, so it needs separate evidence.