The Gram-Schmidt Process

A procedure that turns any independent set of vectors into an orthonormal set spanning the same space.

The idea

Gram-Schmidt builds an orthonormal basis one vector at a time. Start with the first vector and normalize it. For each subsequent vector, subtract its projections onto all the vectors already chosen, leaving a residual orthogonal to them, then normalize that residual. The result is an orthonormal set that spans exactly the same subspace as the original vectors.

The steps

• u1 = v1, then q1 = u1 / |u1|
• u2 = v2 minus its projection onto q1, then q2 = u2 / |u2|
• u3 = v3 minus its projections onto q1 and q2, then normalize
• continue until every input vector is processed

Connection to QR

Gram-Schmidt applied to the columns of a matrix A produces the QR factorization: the orthonormal vectors form Q, and the projection coefficients accumulate into an upper triangular R, so that A = QR. This is the conceptual origin of QR, though practical software uses more stable variants.

Numerical caution

The classical version loses orthogonality when vectors are nearly parallel, because subtracting large nearly-equal quantities amplifies roundoff. The modified Gram-Schmidt process, which subtracts each projection immediately rather than all at once, is more stable. For high accuracy, Householder reflections are preferred over either form.

import numpy as np
def mgs(A):
    A = A.astype(float).copy(); m, n = A.shape
    Q = np.zeros((m, n)); R = np.zeros((n, n))
    for j in range(n):
        v = A[:, j].copy()
        for i in range(j):
            R[i, j] = Q[:, i] @ v; v -= R[i, j] * Q[:, i]
        R[j, j] = np.linalg.norm(v); Q[:, j] = v / R[j, j]
    return Q, R

Orthogonalizing basis functions this way keeps spectral and modal representations of physical fields well conditioned, preventing redundant directions from corrupting a simulation.