The Gram-Schmidt Procedure

Q: What exactly does the procedure guarantee about its output?

Three things: the output vectors are pairwise orthogonal, none of them is the zero vector, and their span equals the span of the inputs. The second conclusion is the only one that requires the input set to be linearly independent, and the third is what allows the output to be substituted for the input in any question about the subspace.

Q: What happens if I run it on a linearly dependent set?

The first vector that lies in the span of its predecessors produces an exactly zero remainder, since everything is projected away. That is a valid dependence detector in exact arithmetic. In floating point the remainder is small rather than zero, and deciding whether it counts as zero is a tolerance judgement that determines the computed rank.

Q: Why does the modified version behave so much better numerically?

Because each coefficient is computed against a vector that has already had earlier components removed, so the coefficient itself is smaller and carries less of the accumulated error. The bound on loss of orthogonality improves from $\varepsilon\kappa(A)^{2}$ to $\varepsilon\kappa(A)$. With a condition number of $10^{8}$ in binary64 that is the difference between a useless and a usable basis.

Q: How does the procedure relate to the QR factorisation?

It is one algorithm for computing it. Collect the inputs as the columns of $A$ and the normalised outputs as the columns of $Q$; the projection coefficients and the norms are exactly the entries of the upper triangular $R$, and $A = QR$. Householder reflections compute the same factorisation by a different and more stable route.

Q: Does the result depend on the order in which I process the vectors?

Yes. The first vector survives untouched and determines the first direction, so any permutation of the inputs generally gives a different orthogonal set. The span of the output is always the same, so every version is equally valid mathematically, but they differ in conditioning and in interpretability.

Q: Should I use Gram-Schmidt or Householder in production code?

Householder, unless the vectors arrive one at a time. Householder QR is backward stable without qualification and needs no reorthogonalisation, at comparable cost. Gram-Schmidt keeps its place in Arnoldi, Lanczos and GMRES iterations, where the next vector does not exist until the previous one has been processed.

Q: Can the procedure be applied to something other than column vectors?

Yes. Nothing in the argument uses the entries of the vectors — only the axioms of the inner product. It therefore runs unchanged in any inner product space, including spaces of polynomials, matrices or functions, which is how families such as the Legendre and Hermite polynomials are constructed.

Overview

An arbitrary basis is awkward to compute with. Coordinates relative to it require solving a linear system, the coefficients of nearby vectors can be wildly different, and nothing decouples. An orthogonal basis has none of these problems: coordinates come from single inner products, energies add without cross terms, and the associated matrix has a trivially computable inverse. The Gram-Schmidt procedure is the constructive answer to the obvious question — given any independent set, can an orthogonal one spanning the same subspace always be produced? It can, and the construction is explicit.

The idea is one step repeated. Keep the first vector unchanged. For each subsequent vector, subtract its projection onto every direction already produced, so that what remains is orthogonal to all of them. Because only multiples of earlier vectors are subtracted, nothing leaves the subspace spanned by the inputs processed so far, and because the input set is independent, the remainder is never zero. Those two observations are the whole content of the theorem: the output is orthogonal, non-zero, and spans exactly what the input spanned.

The span-preservation clause is what makes the procedure useful rather than merely tidy. It means an orthogonal set can be substituted for the original wherever the subspace is what matters — in a column space, a null space basis, a Krylov subspace or a set of regression predictors — without changing the answer to any question about that subspace. In proofs, it licenses the standing assumption that any finite-dimensional inner product space has an orthonormal basis.

In numerical practice the procedure comes with a serious caveat. Written literally as stated, it loses orthogonality catastrophically when the input vectors are nearly dependent: the computed set can drift so far that later vectors are no longer usefully orthogonal to earlier ones. A rearrangement known as modified Gram-Schmidt is mathematically identical and numerically far better, and for production work Householder QR is better still. The classical formulation remains the right way to state and prove the theorem, and the wrong way to implement it.

Definition

Gram-Schmidt Procedure, Column Vectors

GSPCV

Let $S = \{\vec{v}_1, \vec{v}_2, \ldots, \vec{v}_p\}$ be a linearly independent set in $\mathbb{C}^m$ . Define vectors $\vec{u}_1, \ldots, \vec{u}_p$ recursively by

$\vec{u}_i = \vec{v}_i - \frac{\langle \vec{v}_i, \vec{u}_1 \rangle}{\langle \vec{u}_1, \vec{u}_1 \rangle}\vec{u}_1 - \frac{\langle \vec{v}_i, \vec{u}_2 \rangle}{\langle \vec{u}_2, \vec{u}_2 \rangle}\vec{u}_2 - \cdots - \frac{\langle \vec{v}_i, \vec{u}_{i-1} \rangle}{\langle \vec{u}_{i-1}, \vec{u}_{i-1} \rangle}\vec{u}_{i-1}$

Then

T = \{\vec{u}_1, \vec{u}_2, \ldots, \vec{u}_p\}

is an orthogonal set of non-zero vectors and

\operatorname{span}(T) = \operatorname{span}(S)

The first step is degenerate: the sum is empty, so $\vec{u}_1 = \vec{v}_1$ . Every denominator is non-zero because each $\vec{u}_k$ produced so far has been shown to be non-zero.

Orthogonal Projection onto a Vector

OPV

The component of $\vec{v}$ lying along a non-zero $\vec{u}$ is $\operatorname{proj}_{\vec{u}}(\vec{v}) = \frac{\langle \vec{v}, \vec{u} \rangle}{\langle \vec{u}, \vec{u} \rangle}\vec{u}$ . Subtracting it leaves a remainder orthogonal to $\vec{u}$ , which is a one-line verification and is the single operation the whole procedure repeats.

The projection is invariant to the scale of $\vec{u}$ : replacing $\vec{u}$ by $\alpha\vec{u}$ leaves the projection unchanged, so intermediate vectors may be rescaled freely for arithmetic convenience.

Orthonormalisation

The two-stage process of running Gram-Schmidt to obtain an orthogonal set $T$ and then dividing each member by its norm to obtain $\vec{w}_i = \vec{u}_i / \|\vec{u}_i\|$ . The result is an orthonormal set with the same span, and the matrix with these columns satisfies $Q^{*}Q = I_p$ .

Concepts

Why the remainder is orthogonal to everything before it

Fix $i$ and take the inner product of $\vec{u}_i$ with some earlier $\vec{u}_j$ , $j < i$ . Additivity distributes the product across the subtraction, giving $\langle \vec{v}_i, \vec{u}_j \rangle$ minus a sum of terms $\frac{\langle \vec{v}_i, \vec{u}_k \rangle}{\langle \vec{u}_k, \vec{u}_k \rangle}\langle \vec{u}_k, \vec{u}_j \rangle$ . By the inductive hypothesis the earlier vectors are already mutually orthogonal, so every term with $k \ne j$ vanishes. The single surviving term is $\frac{\langle \vec{v}_i, \vec{u}_j \rangle}{\langle \vec{u}_j, \vec{u}_j \rangle}\langle \vec{u}_j, \vec{u}_j \rangle = \langle \vec{v}_i, \vec{u}_j \rangle$ , which cancels the leading term exactly. The coefficient was chosen for precisely this cancellation.

Span preservation runs in both directions

Each $\vec{u}_i$ is $\vec{v}_i$ plus a combination of $\vec{u}_1, \ldots, \vec{u}_{i-1}$ , so by induction every output lies in the span of the inputs and $\operatorname{span}(T) \subseteq \operatorname{span}(S)$ . Rearranging the same equation for $\vec{v}_i$ expresses each input as a combination of outputs, giving the reverse inclusion. The two together force equality. Note the triangular structure: $\vec{u}_i$ involves only $\vec{v}_1$ through $\vec{v}_i$ and vice versa, which is exactly why the matrix $R$ of the associated factorisation is upper triangular.

Non-vanishing is where independence is used

Independence of the input is not decorative. If some $\vec{u}_i$ came out as $\vec{0}$ , rearranging its defining equation would express $\vec{v}_i$ as a combination of $\vec{u}_1, \ldots, \vec{u}_{i-1}$ , hence — by span preservation on the earlier vectors — as a combination of $\vec{v}_1, \ldots, \vec{v}_{i-1}$ . That contradicts independence. Conversely, feeding a dependent set to the procedure produces a zero vector at the first redundant input, which is exactly how the process detects dependence, though in floating point it produces a very small vector rather than an exactly zero one.

From orthogonal to orthonormal, and to QR

Setting $\vec{w}_i = \vec{u}_i/\|\vec{u}_i\|$ normalises without disturbing orthogonality. Collecting the inputs as the columns of $A$ and the normalised outputs as the columns of $Q$ , the rearranged recurrence reads $\vec{v}_i = \sum_{k<i}\langle \vec{v}_i, \vec{w}_k \rangle \vec{w}_k + \|\vec{u}_i\|\vec{w}_i$ , which is precisely $A = QR$ with $[R]_{ki} = \langle \vec{v}_i, \vec{w}_k \rangle$ above the diagonal and $[R]_{ii} = \|\vec{u}_i\|$ on it. Gram-Schmidt is therefore not merely related to the QR factorisation; it is one algorithm for computing it.

The order of the inputs matters

The output is not a function of the input set but of the input sequence. Permuting the inputs generally produces a different orthogonal set, though always with the same span. The first vector is always kept intact, so it determines the first direction entirely. In practice this is exploited: ordering the most reliable or most significant vector first, or applying column pivoting to process the largest remaining vector next, improves both numerical behaviour and interpretability.

Classical form, modified form, and when to abandon both

The classical form computes every coefficient $\langle \vec{v}_i, \vec{u}_k \rangle$ against the original $\vec{v}_i$ and subtracts them all at the end. The modified form subtracts each projection immediately and computes the next coefficient against the partially updated vector. In exact arithmetic the two are identical. In floating point the modified form is dramatically better: its loss of orthogonality grows like $\varepsilon \kappa(A)$ rather than $\varepsilon \kappa(A)^{2}$ , where $\kappa$ is the condition number. The change costs nothing, so there is no reason to implement the classical form.

Even the modified form degrades on badly conditioned input. The standard remedy is to run the orthogonalisation step a second time against the already-computed vectors — the rule of thumb that "twice is enough" is supported by rigorous analysis. Where the full orthonormal basis is the deliverable and accuracy matters more than the ability to produce columns one at a time, Householder QR is superior: it is backward stable unconditionally and produces $Q$ orthogonal to machine precision regardless of conditioning. Gram-Schmidt keeps its place where vectors arrive one at a time, as in Arnoldi and GMRES iterations.

Procedure: orthogonalising a set of $p$ vectors

Confirm independence, or accept detectionThe theorem requires an independent input. If independence is unknown, run the procedure and treat a vanishing

\vec{u}_i

as evidence that

\vec{v}_i

lies in the span of its predecessors.

InitialiseSet

\vec{u}_1 = \vec{v}_1

and record

\langle \vec{u}_1, \vec{u}_1 \rangle

. Nothing is subtracted at the first step.

Project and subtractFor each new

\vec{v}_i

, subtract

\frac{\langle \vec{v}_i, \vec{u}_k \rangle}{\langle \vec{u}_k, \vec{u}_k \rangle}\vec{u}_k

for every

k < i

. In the modified form, update the working vector after each subtraction and compute the next coefficient from the updated vector.

Check the remainderIf

\|\vec{u}_i\|

is negligible relative to

\|\vec{v}_i\|

, the input was numerically dependent. Either stop, or drop that vector and continue with a reduced rank.

Reorthogonalise if requiredRepeat the subtraction pass against the vectors already produced when the norm dropped by more than an order of magnitude. One repetition is normally sufficient.

Normalise and record RDivide each

\vec{u}_i

by its norm to obtain

\vec{w}_i

. The coefficients used, together with the norms, populate the upper triangular factor

R

A = QR

Equations

Gram-Schmidt recurrence

EQ-GSP-01

Gram-Schmidt recurrence

The whole procedure in one line. For $i = 1$ the sum is empty, so $\vec{u}_1 = \vec{v}_1$ .

Projection form of a single step

EQ-GSP-02

{proj}_{\vec{u}} (\vec{v}) = \frac{⟨ \vec{v}, \vec{u} ⟩}{⟨ \vec{u}, \vec{u} ⟩} \vec{u}, ⟨ \vec{v} - {proj}_{\vec{u}} (\vec{v}), \vec{u} ⟩ = 0

Subtracting the projection leaves a remainder orthogonal to the direction projected out. The coefficient is chosen to make exactly this happen.

Conclusions of the procedure

EQ-GSP-03

⟨ {\vec{u}}_{i}, {\vec{u}}_{j} ⟩ = 0 (i \neq j), {\vec{u}}_{i} \neq \vec{0}, span (T) = span (S)

Orthogonality, non-vanishing and span preservation. The second requires the input set to be linearly independent; the other two do not.

Normalisation to an orthonormal set

EQ-GSP-04

{\vec{w}}_{i} = \frac{1}{‖ {\vec{u}}_{i} ‖} {\vec{u}}_{i}, ⟨ {\vec{w}}_{i}, {\vec{w}}_{j} ⟩ = {\begin{cases} 1 & i = j \\ 0 & i \neq j \end{cases}

Scaling by a positive real cannot destroy orthogonality, so normalisation is always safe once the orthogonal set is in hand.

Matrix form: the QR factorisation

EQ-GSP-05

A = Q R, [R]_{k i} = ⟨ {\vec{v}}_{i}, {\vec{w}}_{k} ⟩ (k < i), [R]_{i i} = ‖ {\vec{u}}_{i} ‖

The triangular structure of $R$ is the algebraic shadow of the fact that $\vec{u}_i$ depends only on $\vec{v}_1$ through $\vec{v}_i$ .

Modified Gram-Schmidt update

EQ-GSP-06

Modified Gram-Schmidt update

Identical in exact arithmetic to the classical form, but each coefficient is computed from the partially updated vector, which is what improves the numerical behaviour.

Loss of orthogonality bound

EQ-GSP-07

‖ Q^{*} Q - I ‖ \sim ε κ (A) (modified), ‖ Q^{*} Q - I ‖ \sim ε κ (A)^{2} (classical)

The practical reason to prefer the modified form. With $\kappa(A) = 10^{8}$ in binary64 the classical form returns columns that are not orthogonal at all.

Variable Definitions

Symbols used on this page
Symbol	Name	Meaning	Domain / type
$S$	Input set	The linearly independent set $\{\vec{v}_1, \ldots, \vec{v}_p\}$ to be orthogonalised	subset of C^m
$T$	Output set	The orthogonal set $\{\vec{u}_1, \ldots, \vec{u}_p\}$ with the same span	subset of C^m
${\vec{v}}_{i}$	Input vector	The $i$ -th vector of the original set	element of C^m
${\vec{u}}_{i}$	Orthogonal vector	The $i$ -th output, orthogonal to all its predecessors	element of C^m
${\vec{w}}_{i}$	Orthonormal vector	$\vec{u}_i$ scaled to unit norm; column $i$ of $Q$	element of C^m
$p$	Set size	Number of vectors being orthogonalised, at most $m$	positive integer
$m$	Vector size	Number of entries in each vector	positive integer
$R$	Triangular factor	Upper triangular matrix of projection coefficients and norms	p x p matrix
$κ (A)$	Condition number	Ratio of largest to smallest singular value of the input matrix; controls loss of orthogonality	real, at least 1

Worked Numerical Example

Problem statement

A test rig samples a response at four equally spaced instants $t = 0, 1, 2, 3$ . The natural model basis is $\{1,\ t,\ t^{2}\}$ , which sampled at those instants gives three column vectors in $\mathbb{R}^4$ . These columns are badly correlated, so orthogonalise them by Gram-Schmidt, normalise the result, and use the resulting orthogonal basis to fit measured data without solving a linear system.

Sample the model basis
Evaluating $1$ , $t$ and $t^{2}$ at $t = 0, 1, 2, 3$ gives three columns. Their pairwise inner products are large and positive, so the design matrix is strongly correlated — the classic ill-conditioning of a raw polynomial fit.
${\vec{v}}_{1} = [\begin{matrix} 1 \\ 1 \\ 1 \\ 1 \end{matrix}], {\vec{v}}_{2} = [\begin{matrix} 0 \\ 1 \\ 2 \\ 3 \end{matrix}], {\vec{v}}_{3} = [\begin{matrix} 0 \\ 1 \\ 4 \\ 9 \end{matrix}]$
First vector passes through unchanged
There is nothing to project against, so $\vec{u}_1 = \vec{v}_1$ . Record $\langle \vec{u}_1, \vec{u}_1 \rangle = 4$ for reuse; every later step divides by it.
${\vec{u}}_{1} = [\begin{matrix} 1 \\ 1 \\ 1 \\ 1 \end{matrix}], ⟨ {\vec{u}}_{1}, {\vec{u}}_{1} ⟩ = 4$
Remove the mean from the linear term
The coefficient is $\langle \vec{v}_2, \vec{u}_1 \rangle / \langle \vec{u}_1, \vec{u}_1 \rangle = 6/4 = 3/2$ , which is the arithmetic mean of the sample times. Subtracting it centres the linear column, and the result is orthogonal to the constant column by construction.
${\vec{u}}_{2} = {\vec{v}}_{2} - \frac{3}{2} {\vec{u}}_{1} = [\begin{matrix} - 3 / 2 \\ - 1 / 2 \\ 1 / 2 \\ 3 / 2 \end{matrix}], ⟨ {\vec{u}}_{2}, {\vec{u}}_{2} ⟩ = 5$
Project the quadratic term out of both earlier directions
Two coefficients are needed: $\langle \vec{v}_3, \vec{u}_1 \rangle / 4 = 14/4 = 7/2$ and $\langle \vec{v}_3, \vec{u}_2 \rangle / 5 = 15/5 = 3$ . Subtracting both leaves the pure quadratic direction, with all constant and linear content removed.
${\vec{u}}_{3} = {\vec{v}}_{3} - \frac{7}{2} {\vec{u}}_{1} - 3 {\vec{u}}_{2} = [\begin{matrix} 1 \\ - 1 \\ - 1 \\ 1 \end{matrix}], ⟨ {\vec{u}}_{3}, {\vec{u}}_{3} ⟩ = 4$
Verify orthogonality explicitly
Never trust the construction without checking. $\langle \vec{u}_2, \vec{u}_1 \rangle = -\tfrac{3}{2} - \tfrac{1}{2} + \tfrac{1}{2} + \tfrac{3}{2} = 0$ ; $\langle \vec{u}_3, \vec{u}_1 \rangle = 1 - 1 - 1 + 1 = 0$ ; and $\langle \vec{u}_3, \vec{u}_2 \rangle = -\tfrac{3}{2} + \tfrac{1}{2} - \tfrac{1}{2} + \tfrac{3}{2} = 0$ . All three distinct pairs vanish, so $T$ is orthogonal and therefore linearly independent.
Normalise and read off the QR factors
The norms are $2$ , $\sqrt{5}$ and $2$ . Dividing gives the orthonormal columns of $Q$ , and the coefficients already computed populate $R$ : the diagonal holds the norms and the off-diagonal entries hold $\langle \vec{v}_i, \vec{w}_k \rangle$ .
$Q = [\begin{matrix} 1 / 2 & - 3 / (2 \sqrt{5}) & 1 / 2 \\ 1 / 2 & - 1 / (2 \sqrt{5}) & - 1 / 2 \\ 1 / 2 & 1 / (2 \sqrt{5}) & - 1 / 2 \\ 1 / 2 & 3 / (2 \sqrt{5}) & 1 / 2 \end{matrix}], R = [\begin{matrix} 2 & 3 & 7 \\ 0 & \sqrt{5} & 3 \sqrt{5} \\ 0 & 0 & 2 \end{matrix}]$
Check the factorisation on one column
Reconstruct the third column of the input from the third column of $R$ : $7\vec{w}_1 + 3\sqrt{5}\,\vec{w}_2 + 2\vec{w}_3$ evaluates entry by entry to $\begin{bmatrix} 0 & 1 & 4 & 9 \end{bmatrix}^{t}$ , which is $\vec{v}_3$ . A single column check of this kind catches most transcription and sign errors.
Fit measured data with decoupled coefficients
Suppose the rig returns $\vec{y} = \begin{bmatrix} 1 & 2 & 5 & 10 \end{bmatrix}^{t}$ . Because the basis is orthogonal, each coefficient is an independent quotient of inner products, with no normal equations to solve: $c_1 = 18/4 = 9/2$ , $c_2 = 15/5 = 3$ , $c_3 = 4/4 = 1$ .
$\vec{y} = \frac{9}{2} {\vec{u}}_{1} + 3 {\vec{u}}_{2} + 1 {\vec{u}}_{3}$

Result

The correlated sampled monomials have been replaced by three mutually orthogonal columns — a constant, a centred ramp and an alternating quadratic contrast — that span exactly the same space of quadratic responses. Fitting the data required three inner products instead of a $3 \times 3$ solve, and the coefficients are stable: adding a cubic term to the model would leave all three unchanged, which is never true of a raw monomial basis. The leading coefficient $c_3 = 1$ recovers the fact that the data is exactly $1 + t^{2}$ .

Applications & Industry Use

Statistics and regression

Orthogonal polynomial fitting

Fitting a polynomial in raw powers produces a design matrix with a condition number that grows explosively with degree. Orthogonalising the sampled powers gives discrete orthogonal polynomials for which the normal equations are diagonal, coefficients are independent of the model order, and the fit is numerically stable to far higher degree.

Iterative linear solvers

Arnoldi and GMRES

Krylov subspace methods generate basis vectors one at a time by repeated matrix-vector products, then orthogonalise each new vector against all its predecessors. Gram-Schmidt is used rather than Householder precisely because the vectors are not available in advance, and reorthogonalisation is standard practice to stop the basis degenerating.

Wireless communications

Interference cancellation and beamforming

A receiver can null a known interferer by projecting the array snapshot onto the orthogonal complement of the interference direction — a single Gram-Schmidt step. Successive interference cancellation applies the same operation repeatedly, orthogonalising each user's channel against those already decoded.

Computational chemistry and physics

Orthogonalising basis functions

Atomic orbital basis sets overlap, so the associated overlap matrix is not the identity. Orthogonalisation, whether by a Gram-Schmidt sweep or by symmetric methods, produces a basis in which the eigenvalue problem takes standard form, which is a prerequisite for most electronic structure algorithms.

Computer graphics and robotics

Re-orthonormalising drifting frames

A rotation matrix integrated over many time steps slowly loses orthonormality, so that lengths and angles start to distort. A short Gram-Schmidt pass over its three columns restores the frame cheaply, and is the standard correction applied in attitude-tracking and physics-integration loops.

Geodesy and surveying

Least-squares network adjustment

Adjustment of a measurement network is a large least-squares problem solved through a QR factorisation rather than by forming the normal equations, because squaring the design matrix doubles the condition number. The orthogonalisation stage of that factorisation is the procedure described here.

Design Considerations

Never implement the classical form

The modified formulation costs the same number of operations, requires only a reordering of loops, and reduces the loss of orthogonality from $\varepsilon\kappa(A)^{2}$ to $\varepsilon\kappa(A)$ . There is no scenario in which the classical arrangement is preferable in floating point, and it should appear only in proofs and textbook statements.

Choose Householder when the full basis is the deliverable

Householder QR is unconditionally backward stable and produces $Q$ orthogonal to machine precision no matter how ill-conditioned the input, at essentially the same cost. Prefer Gram-Schmidt only when vectors arrive incrementally, as in Krylov methods, or when the projection coefficients themselves have interpretive value.

Budget for reorthogonalisation

A cheap and effective safeguard is to monitor the norm of the working vector during the subtraction pass. If it drops by more than about an order of magnitude, catastrophic cancellation has occurred and the orthogonalisation should be repeated against the same vectors. Doing this at most once is sufficient in practice and roughly doubles the worst-case cost.

Order the input deliberately

The procedure preserves the first vector exactly and progressively degrades later ones, so the sequence should be chosen rather than inherited. Placing the largest or most trusted vector first, or applying column pivoting to select the vector with the largest remaining norm at each step, improves conditioning and yields a rank-revealing factorisation as a by-product.

Decide what a near-zero remainder means

In exact arithmetic a zero remainder proves dependence. In floating point the remainder is merely small, and the threshold at which it is declared zero determines the computed rank. That threshold must be scaled — typically against $\|\vec{v}_i\|$ or the largest diagonal of $R$ — and stated explicitly wherever a rank is reported.

Keep the unnormalised form if square roots are expensive

Normalisation introduces $p$ square roots and $p$ divisions. Where the orthogonal set is used only for projections, the unnormalised form with explicit $\langle \vec{u}_k, \vec{u}_k \rangle$ denominators is exact-arithmetic friendly and avoids irrational entries entirely, which matters for symbolic and fixed-point work.

Standards & Codes

Notation, interchange and numerical standards that govern how this material is written down, stored and computed in production systems.

Applicable standards, conventions and reference implementations
Reference	Title	Relevance to this topic
`LAPACK / BLAS reference`	Linear Algebra PACKage reference implementation	Provides `xGEQRF` (Householder QR), `xGEQP3` (pivoted QR) and `xORGQR` / `xUNGQR` for forming $Q$ explicitly. LAPACK deliberately does not ship a classical Gram-Schmidt routine, which is itself a standards-level judgement about the algorithm.
`BLAS Level 2`	Basic Linear Algebra Subprograms, matrix-vector operations	Modified Gram-Schmidt is naturally expressed as repeated `xGEMV` and `xAXPY` calls; block variants promote the inner loop to Level 3 `xGEMM` for cache efficiency at the cost of extra reorthogonalisation.
`IEEE 754-2019`	IEEE Standard for Floating-Point Arithmetic	Defines the unit roundoff $\varepsilon$ that appears in every loss-of-orthogonality bound, and the cancellation behaviour that makes subtraction of nearly equal vectors the critical step of the algorithm.
`ISO 80000-2`	Quantities and units — Part 2: Mathematics	Fixes the notation for inner products, norms and the span of a set used throughout the statement of the procedure and its conclusions.
`ISO/IEC 14882`	Programming languages — C++	Specifies `std::inner_product` and the numeric requirements that a hand-rolled orthogonalisation loop relies on; template libraries such as Eigen expose `HouseholderQR` and `ColPivHouseholderQR` in preference to a Gram-Schmidt implementation.

Material Selection

For a mathematical topic, "material" is the numeric representation: the scalar field, storage format and precision the computation is built from.

Representation and precision selection
Representation	Select when	Trade-off
Exact rational arithmetic	Small symbolic problems, textbook verification, and construction of exact orthogonal bases for integer data.	Orthogonality is exact and dependence is detected with certainty, but fractions proliferate rapidly: denominators multiply at every step and can grow beyond practical size by around ten vectors.
IEEE 754 binary64 with modified Gram-Schmidt	General numerical use where vectors arrive incrementally.	The standard working choice; orthogonality holds to about $\varepsilon\kappa(A)$ , which is acceptable up to condition numbers near $10^{8}$ and inadequate beyond.
IEEE 754 binary64 with reorthogonalisation	Condition numbers above $10^{8}$ , or long Krylov recurrences where the basis must stay orthogonal over many iterations.	Restores near machine-precision orthogonality at up to twice the cost; selective reorthogonalisation triggered by a norm-drop test recovers most of the benefit for a small fraction of the cost.
Householder reflectors, binary64	The whole matrix is available and a numerically defensible orthonormal basis is required.	Backward stable regardless of conditioning and no reorthogonalisation needed, but the basis is produced only at the end and $Q$ is stored implicitly unless explicitly generated.
IEEE 754 binary32	Real-time attitude correction and graphics frames, where only three or four short vectors are involved.	Fast and adequate for well-conditioned small frames, but the unit roundoff of about $10^{-7}$ leaves little margin, so periodic re-orthonormalisation must be frequent.
Complex scalars	Array processing, spectral bases and quantum states.	Every coefficient is complex and the conjugate must fall on the correct argument; the algorithm is otherwise unchanged, at roughly four times the real arithmetic cost per inner product.

Manufacturing Notes

Implementation notes — how the result is actually produced by hand, by algorithm and by library, including cost and numerical behaviour.

Operation count

Orthogonalising $p$ vectors of length $m$ requires about $2mp^{2}$ floating-point operations, the same leading order as Householder QR for the factorisation itself. Step $i$ costs $O(mi)$ , so the work grows linearly through the sweep — a useful property when the number of vectors is not known in advance, as in an iterative solver that may terminate early.

Hand procedure

Tabulate $\langle \vec{u}_k, \vec{u}_k \rangle$ as each vector is completed; those denominators are reused at every later step and recomputing them is the commonest source of arithmetic slips. Rescaling an intermediate $\vec{u}_k$ by any convenient factor to clear fractions is legitimate and does not change the subsequent results, because the projection is invariant to the scale of the vector projected onto.

Library behaviour

numpy.linalg.qr and scipy.linalg.qr both use Householder reflections via LAPACK, not Gram-Schmidt; scipy.linalg.orth uses the singular value decomposition. SymPy's GramSchmidt works in exact arithmetic and is intended for symbolic use only. A hand-written Gram-Schmidt in a production numerical pipeline should be treated as a red flag unless the incremental structure genuinely requires it.

Verification technique

Three checks close the loop. Compute the full Gram matrix and inspect $\|Q^{*}Q - I_p\|$ ; reconstruct $Q R$ and compare with the original $A$ in the Frobenius norm; and confirm that the diagonal of $R$ is positive and decreasing in magnitude if pivoting was applied. The factorisation residual can be small even when orthogonality has been lost, so both tests are needed — neither implies the other.

Cancellation is the failure mechanism

When $\vec{v}_i$ is nearly inside the span of its predecessors, the subtraction removes almost all of it and the surviving remainder consists largely of rounding error. The leading digits cancel, the relative error explodes, and the resulting direction is essentially arbitrary. Watching the ratio $\|\vec{u}_i\| / \|\vec{v}_i\|$ gives direct visibility of this and is the cheapest available diagnostic.

Failure Modes & Common Mistakes

Failure modes, root causes and prevention
Failure mode / mistake	Impact	Root cause	Prevention & detection
Implementing the classical rather than modified form	high	Coding the recurrence exactly as written, computing all coefficients against the original vector before subtracting.	Subtract each projection immediately and compute the next coefficient from the updated vector; the loop reordering is trivial and the accuracy gain is orders of magnitude.
Feeding a linearly dependent set	high	Running the procedure on columns that are not independent, producing a zero or near-zero vector whose normalisation then divides by almost nothing.	Test $\\|\vec{u}_i\\|$ against a scaled tolerance before normalising, and either terminate or drop the offending vector and continue with reduced rank.
Silent loss of orthogonality on ill-conditioned input	high	Cancellation during subtraction leaves a remainder dominated by rounding error, so later vectors are not orthogonal to earlier ones.	Monitor the norm drop at each step and reorthogonalise when it exceeds an order of magnitude; report $\\|Q^{*}Q - I\\|$ with the result.
Reusing stale denominators	medium	Caching $\langle \vec{u}_k, \vec{u}_k \rangle$ and failing to update it after rescaling or reorthogonalising $\vec{u}_k$ .	Recompute the denominator whenever the vector it belongs to is modified, or normalise immediately so that every denominator is $1$ .
Assuming the output is independent of input order	medium	Treating the procedure as a function of a set rather than a sequence, then being surprised when a permuted input gives a different basis.	Fix and document the ordering; where a canonical result is required, use pivoting or the singular value decomposition, which does not depend on column order.
Putting the conjugate on the wrong argument	high	Complex data with the coefficient computed as $\langle \vec{u}_k, \vec{v}_i \rangle$ instead of $\langle \vec{v}_i, \vec{u}_k \rangle$ .	Verify orthogonality of the first two output vectors explicitly on a complex test case; a swapped conjugate leaves a purely imaginary residual that is easy to spot.
Normalising before the whole sweep is complete	low	Scaling $\vec{u}_i$ to unit norm partway through and then continuing with denominators that assume the unnormalised vector.	Either normalise every vector immediately as it is produced and drop all denominators, or normalise nothing until the end. Do not mix the two conventions.
Forming normal equations instead of using QR	medium	Solving least squares via $A^{*}A$ , which squares the condition number and discards the accuracy that orthogonalisation was meant to provide.	Solve the triangular system $R\vec{x} = Q^{}\vec{b}$ directly; never form $A^{}A$ for an ill-conditioned design matrix.

FAQs

What exactly does the procedure guarantee about its output?

Three things: the output vectors are pairwise orthogonal, none of them is the zero vector, and their span equals the span of the inputs. The second conclusion is the only one that requires the input set to be linearly independent, and the third is what allows the output to be substituted for the input in any question about the subspace.

What happens if I run it on a linearly dependent set?

The first vector that lies in the span of its predecessors produces an exactly zero remainder, since everything is projected away. That is a valid dependence detector in exact arithmetic. In floating point the remainder is small rather than zero, and deciding whether it counts as zero is a tolerance judgement that determines the computed rank.

Why does the modified version behave so much better numerically?

Because each coefficient is computed against a vector that has already had earlier components removed, so the coefficient itself is smaller and carries less of the accumulated error. The bound on loss of orthogonality improves from $\varepsilon\kappa(A)^{2}$ to $\varepsilon\kappa(A)$ . With a condition number of $10^{8}$ in binary64 that is the difference between a useless and a usable basis.

How does the procedure relate to the QR factorisation?

It is one algorithm for computing it. Collect the inputs as the columns of $A$ and the normalised outputs as the columns of $Q$ ; the projection coefficients and the norms are exactly the entries of the upper triangular $R$ , and $A = QR$ . Householder reflections compute the same factorisation by a different and more stable route.

Does the result depend on the order in which I process the vectors?

Yes. The first vector survives untouched and determines the first direction, so any permutation of the inputs generally gives a different orthogonal set. The span of the output is always the same, so every version is equally valid mathematically, but they differ in conditioning and in interpretability.

Should I use Gram-Schmidt or Householder in production code?

Householder, unless the vectors arrive one at a time. Householder QR is backward stable without qualification and needs no reorthogonalisation, at comparable cost. Gram-Schmidt keeps its place in Arnoldi, Lanczos and GMRES iterations, where the next vector does not exist until the previous one has been processed.

Can the procedure be applied to something other than column vectors?

Yes. Nothing in the argument uses the entries of the vectors — only the axioms of the inner product. It therefore runs unchanged in any inner product space, including spaces of polynomials, matrices or functions, which is how families such as the Legendre and Hermite polynomials are constructed.

References

Beezer, R. A. A First Course in Linear Algebra, Version 0.70. University of Puget Sound, 2006. Section O, subsection GSP. Licensed under the GNU Free Documentation License v1.2.
Björck, Å. Numerical Methods for Least Squares Problems. Society for Industrial and Applied Mathematics, 1996.
Golub, G. H. and Van Loan, C. F. Matrix Computations, 4th edition. Johns Hopkins University Press, 2013.
Trefethen, L. N. and Bau, D. Numerical Linear Algebra. Society for Industrial and Applied Mathematics, 1997.
Higham, N. J. Accuracy and Stability of Numerical Algorithms, 2nd edition. Society for Industrial and Applied Mathematics, 2002.
ISO 80000-2:2019, Quantities and units — Part 2: Mathematics. International Organization for Standardization.

AI Suggested Questions

Show me a 3-column matrix with condition number near $10^{8}$ and compare the orthogonality loss of classical against modified Gram-Schmidt in binary64.
Derive the entries of $R$ in the QR factorisation directly from the Gram-Schmidt recurrence and explain why the matrix is upper triangular.
Apply the procedure to the monomials $1, t, t^{2}, t^{3}$ under the integral inner product on $[-1,1]$ and identify the polynomials it produces.
How does column pivoting change the Gram-Schmidt sweep, and why does it make the factorisation rank-revealing?
Explain the norm-drop criterion for triggering reorthogonalisation and justify the rule that twice is enough.
Why does forming the normal equations $A^{*}A$ destroy the accuracy that a QR-based least-squares solve preserves?

Related Calculators

Gram-Schmidt Calculator

Orthogonalise a set of real or complex vectors step by step, showing every projection coefficient and the resulting orthonormal set.

QR Factorisation Calculator

Factor a matrix as $A = QR$ by modified Gram-Schmidt or Householder reflections and report the reconstruction residual.

Orthogonality Loss Monitor

Measure $\|Q^{*}Q - I\|$ against the condition number of the input and recommend when reorthogonalisation is required.