Modified Gram-Schmidt#
Householder reflections and Givens rotations compute QR by introducing zeros in a matrix. Gram-Schmidt takes a different approach: it builds an orthonormal basis for the column space, one vector at a time. The coefficients used to express the original columns in this basis form the triangular factor \(R\).
This approach is useful when vectors arrive sequentially. For example, an iterative method may generate its next vector using the basis vectors already computed. Gram-Schmidt can extend that basis without changing the earlier vectors. We will use this construction later in Arnoldi iteration. Householder QR can also be extended, but it stores the basis implicitly as a product of transformations.
Throughout this page, \(A\in\mathbb{R}^{m\times n}\) with \(m\ge n\) has linearly independent columns. We seek the reduced QR factorization
where \(R\in\mathbb{R}^{n\times n}\) is upper triangular with positive diagonal entries. Since \(R\) is invertible, \(Q=AR^{-1}\). This viewpoint is sometimes called triangular orthogonalization, in contrast to applying orthogonal transformations to triangularize \(A\). The algorithm does not form \(R^{-1}\).
Removing One Direction at a Time#
Let \(a_j\) denote the original \(j\)-th column of \(A\). The first column of \(A=QR\) satisfies \(a_1=r_{11}q_1\), so
For every remaining column, remove its component along \(q_1\):
Here \(a_j^{(1)}\) is the remaining part of \(a_j\) after the first projection. It is perpendicular to \(q_1\), because
Normalize \(a_2^{(1)}\) to obtain \(q_2\), then remove the \(q_2\) component from each remaining column. Continuing in this way gives modified Gram-Schmidt (MGS). Each projection coefficient is computed from the current residual column, after the preceding components have been removed.
Fig. 10 Remove the component of \(a_2\) along \(q_1\), then normalize what remains. On the left, the purple arrow represents \(a_2^{(1)}\), drawn from the tip of the projection to the tip of \(a_2\). On the right, the same remainder is drawn from the origin; normalizing it changes its length but preserves its perpendicular direction. Both \(q_1\) and \(q_2\) have length \(1\).#
At step \(k\), write \(a_j^{(k-1)}\) for this residual, with \(a_j^{(0)}=a_j\). The formulas are
The update is an orthogonal projection onto the space perpendicular to \(q_k\):
There is also an outer-product interpretation. If \(r_k^T\) denotes row \(k\) of \(R\), then
Each step extracts one rank-one term. Since \(r_{kj}=0\) for \(j<k\), only the remaining columns need to be updated.
The In-Place Algorithm#
The superscripts above explain the derivation. In an implementation, one working array is enough. Its column \(k\) is normalized and stored as \(q_k\); its later columns hold the current residuals. Initialize \(R\) to zero.
for k = 1, ..., n:
r[k,k] = norm(A[:,k], 2)
if r[k,k] == 0:
stop: the current column adds no independent direction
A[:,k] = A[:,k] / r[k,k]
for j = k+1, ..., n:
r[k,j] = A[:,k]^T A[:,j]
A[:,j] = A[:,j] - A[:,k] r[k,j]
On completion, the working array \(A\) contains \(Q\), and \(R\) is stored separately. The input matrix is overwritten; retain a copy if it is needed later.
The zero test describes breakdown in exact arithmetic. In floating-point arithmetic, a tiny residual may contain mostly rounding error even when its norm is nonzero. Detecting numerical dependence requires a tolerance relative to the input scale and desired accuracy. A small \(r_{kk}\) alone does not establish that the original columns are dependent.
Why the Algorithm Works#
The proof needs to establish three facts: the vectors are orthonormal, normalization never divides by zero under our rank assumption, and the stored coefficients reconstruct \(A\).
Proof. Correctness of Modified Gram-Schmidt in Exact Arithmetic. Suppose that \(q_1,\ldots,q_{k-1}\) are orthonormal and span the first \(k-1\) original columns. At the start of step \(k\), each remaining column satisfies
For \(k=1\), the sum is empty and these statements hold with \(a_j^{(0)}=a_j\).
The residual \(a_k^{(k-1)}\) cannot be zero. Otherwise, the displayed identity would place \(a_k\) in the span of its preceding columns, contradicting linear independence. Thus \(r_{kk}>0\), and normalizing gives a unit vector \(q_k\) perpendicular to the previous vectors. Also,
so the first \(k\) basis vectors span the first \(k\) original columns.
For \(j>k\), the update removes the component along \(q_k\):
It also preserves perpendicularity to every earlier \(q_i\), since
Finally, substituting \(a_j^{(k-1)}=q_kr_{kj}+a_j^{(k)}\) into the reconstruction identity extends its sum through \(i=k\). This proves the induction step.
At completion, every original column satisfies
These column identities give \(A=QR\), with \(Q^TQ=I_n\) and \(R\) upper triangular with positive diagonal.
If the columns are dependent, the basic algorithm stops at the first zero residual. QR still exists, but completing an orthonormal basis requires an additional choice. This is one reason to distinguish existence of QR from successful completion of this particular algorithm.
Comparison with Classical Gram-Schmidt#
Classical Gram-Schmidt (CGS) computes all projection coefficients for a column from the original column:
It then sets \(r_{jj}=\|v\|_2\) and \(q_j=v/r_{jj}\). MGS instead uses the successive residuals. For one column, its computation can be written as
This column-by-column ordering computes the same MGS projections as the in-place algorithm above. It is convenient when new vectors are generated one at a time.
In exact arithmetic, CGS and MGS agree. Removing a component along \(q_\ell\) does not change the inner product with another orthogonal vector \(q_i\):
In floating-point arithmetic, the computed basis vectors are not exactly orthogonal, so this equality no longer holds exactly. The choice of which vector enters each dot product then matters.
Numerical Accuracy: The Condition-Number Factor#
Use the two backward-error measures introduced for Householder QR. For MGS, \(\widehat Q\) is the computed reduced factor, so its orthogonality equation uses \(I_n\).
For a full-column-rank matrix, let
where \(\sigma_{\max}\) and \(\sigma_{\min}\) are the largest and smallest singular values. Under the same rounding assumptions, MGS gives the bounds
The \(C_i(m,n)\) are dimension-dependent factors independent of \(A\). The orthogonality bound additionally requires \(C_0(m,n)u\kappa_2(A)\) to be sufficiently small, so rounding errors remain below the scale of the smallest singular value.
The new factor is \(\kappa_2(A)\) in the orthogonality bound. MGS still reconstructs \(A\) accurately, but nearly dependent columns can cause substantial loss of orthogonality. Householder and Givens have no condition-number factor in either bound.
The reason is cancellation in the projection steps. A nearly dependent column leaves a residual much smaller than the quantities subtracted from it. Errors in those quantities can be large relative to the residual, and normalization turns them into errors in the new basis direction. MGS usually reduces this effect compared with one-pass CGS, but does not eliminate it. This motivates a second orthogonalization pass.
See Higham’s QR overview for the comparison, and Chapter 19 of his Accuracy and Stability of Numerical Algorithms for the analysis.
Reorthogonalization#
If accurate orthogonality is needed, project each residual against the previous basis vectors a second time before normalizing. This is called reorthogonalization. The second pass removes components left by rounding errors in the first pass.
For column \(j\), first perform the column-by-column MGS pass above. Then, using its remaining vector \(v\), compute
Finally set \(r_{jj}=\|v\|_2\) and \(q_j=v/r_{jj}\). Adding the correction to \(r_{ij}\) is essential: both passes remove components that must be included in the identity \(a_j=\sum_{i<j}q_ir_{ij}+q_jr_{jj}\).
For numerically independent columns, two passes can restore orthogonality to a dimension-dependent multiple of \(u\). This guarantee requires rounding errors to be small relative to the independent part of each new column; a second pass cannot recover a direction that has already been lost in rounding. See Giraud, Langou, and Rozloznik’s analysis of reorthogonalization for the precise assumptions and bounds.
A Numerical Example#
Consider three almost parallel columns:
We have \(A_\varepsilon^TA_\varepsilon=\boldsymbol{1}\boldsymbol{1}^T+\varepsilon^2I_3\), where \(\boldsymbol{1}=(1,1,1)^T\). The rank-one matrix \(\boldsymbol{1}\boldsymbol{1}^T\) has eigenvalue \(3\) along \(\boldsymbol{1}\) and eigenvalue \(0\) on its perpendicular complement. Adding \(\varepsilon^2I_3\) and taking square roots gives the singular values \(\sqrt{3+\varepsilon^2},\varepsilon,\varepsilon\). Thus the columns are independent, but \(\kappa_2(A_\varepsilon)=\sqrt{3+\varepsilon^2}/\varepsilon\) becomes large as \(\varepsilon\) decreases.
The following experiment compares one-pass CGS, one-pass MGS, MGS with a second pass, and Householder QR. It uses the column-by-column ordering so the second pass is easy to identify.
import numpy as np
def gram_schmidt(A, classical=False, passes=1):
m, n = A.shape
Q = np.zeros((m, n))
R = np.zeros((n, n))
for j in range(n):
v = A[:, j].copy()
for _ in range(passes):
if classical:
# All coefficients use the same vector.
coefficients = Q[:, :j].T @ v
R[:j, j] += coefficients
v -= Q[:, :j] @ coefficients
else:
# Each coefficient uses the updated residual.
for i in range(j):
coefficient = Q[:, i] @ v
R[i, j] += coefficient
v -= coefficient * Q[:, i]
R[j, j] = np.linalg.norm(v)
if R[j, j] == 0:
raise ValueError("Zero residual: cannot normalize this column.")
Q[:, j] = v / R[j, j]
return Q, R
epsilon = 1e-8
A = np.vstack((np.ones((1, 3)), epsilon * np.eye(3)))
factorizations = [
("CGS", gram_schmidt(A, classical=True)),
("MGS", gram_schmidt(A)),
("MGS, two passes", gram_schmidt(A, passes=2)),
("Householder", np.linalg.qr(A, mode="reduced")),
]
print(f"{'Method':<18} {'Factorization error':>20} {'Orthogonality error':>20}")
for name, (Q, R) in factorizations:
residual = np.linalg.norm(A - Q @ R, "fro") / np.linalg.norm(A, "fro")
orthogonality = np.linalg.norm(np.eye(3) - Q.T @ Q, 2)
print(f"{name:<18} {residual:20.1e} {orthogonality:20.1e}")
Method Factorization error Orthogonality error
CGS 8.3e-26 5.0e-01
MGS 1.4e-25 8.2e-09
MGS, two passes 6.4e-17 3.7e-17
Householder 3.2e-24 6.9e-16
All four methods reconstruct this input accurately. Their orthogonality errors differ substantially: MGS improves on CGS, while the second pass brings MGS close to Householder accuracy. The exact printed values can vary slightly with the numerical library.
Computational Cost and Choice of Method#
For each pair \(k<j\), the dot product costs about \(2m\) flops, and the residual update costs another \(2m\). There are \(n(n-1)/2\) such pairs. Computing a norm and normalizing each column costs about \(3m\) flops, so the total is
CGS has the same leading arithmetic cost. A second MGS pass doubles the projection work, giving \(4mn^2+O(mn)\) flops when used for every column.
Householder QR costs \(2mn^2-\tfrac23n^3+O(mn)\) flops to compute \(R\) and store \(Q\) implicitly; explicitly forming the reduced \(Q\) requires additional work. MGS produces the reduced \(Q\) directly. Thus being tall and thin does not by itself make MGS the preferred method. Householder QR is the usual choice for a dense matrix available in advance; MGS with reorthogonalization is useful when an explicit basis must be built sequentially.
For complex matrices, replace transposes in the inner products by conjugate transposes. The projection formulas and exact-arithmetic proof then apply with \(Q^HQ=I_n\).