Householder QR and Orthogonality Loss in Gram-Schmidt

Orthogonalizing a set of vectors sounds like a solved problem: subtract from each vector its projection onto the ones already processed. Implement that literally and the resulting basis can be orthogonal to about half the digits the arithmetic carries, on matrices that are only moderately ill conditioned. The failure is not a bug in the loop; it is the arithmetic consequence of computing the projection coefficients from the same vectors that are being modified.

The QR factorization computed with Householder reflectors avoids the failure, at roughly twice the flop count. The trade is one of the clearest in numerical linear algebra: pay a constant factor in arithmetic to remove a dependence on the condition number from the orthogonality of the result.

Two versions of the same projection

Classical Gram-Schmidt computes all the coefficients for a new vector against the existing basis, then subtracts all the projections at once. Modified Gram-Schmidt computes one coefficient, subtracts that one projection, then recomputes the next coefficient against the already-modified vector. In exact arithmetic the two produce identical results.

In floating point they differ in the size of the errors that the coefficients inherit. When the vectors are nearly dependent, the projection coefficients are large relative to the components that remain after subtraction, and computing them from the unmodified vector means computing a large number as the difference of two nearly equal quantities. Modified Gram-Schmidt recomputes each coefficient from the partially reduced vector, so the cancellation happens before the coefficient is formed rather than after.

The standard error bounds make the gap explicit. For an m x n matrix with condition number kappa, classical Gram-Schmidt produces a computed basis whose deviation from orthogonality, measured as the norm of Q^T Q - I, is of order kappa^2 * u, while modified Gram-Schmidt produces one of order kappa * u. Both lose orthogonality in proportion to the conditioning, one squared and one not. For a matrix with kappa = 1e8 in double precision, the classical version can lose all of the sixteen digits; the modified version retains about eight.

The reflector, and where it is stored

Householder QR takes a different route. Instead of orthogonalizing column by column against a growing basis, it applies a sequence of reflections that zero out the entries below the diagonal one column at a time. Each reflection is defined by a vector v:

H = I - tau * v * v^T,    tau = 2 / (v^T v)

and the vector is chosen so that H x is a multiple of the first coordinate axis. For a vector x with x_1 not zero, the choice v = x + sign(x_1) * ||x||_2 * e_1 produces H x = -sign(x_1) * ||x||_2 * e_1. The sign is chosen to match x_1 so that the first component of v is a sum rather than a difference, which keeps the reflector away from cancellation.

After n steps, R is the upper triangular result and Q is the product of the reflectors in reverse order. The reflectors are never formed as matrices. Each is stored as the vector v, normalized so that its first entry is one, in the zeroed portion of the column it annihilated, together with the scalar tau. LAPACK’s geqrf uses exactly this layout, which is why the routine’s output overwrites the input array.

#include <math.h>

/* Householder QR of a column-major m x n matrix (m >= n).
   On return the upper triangle holds R; below the diagonal are the
   reflector vectors with implicit leading 1; tau holds the scalars. */
void qr_householder(int m, int n, double *A, double *tau)
{
    for (int k = 0; k < n; k++) {
        double *x = &A[k * m + k];
        double norm = 0.0;
        for (int i = 0; i < m - k; i++)
            norm += x[i] * x[i];
        norm = sqrt(norm);
        if (norm == 0.0) { tau[k] = 0.0; continue; }

        double alpha = (x[0] >= 0.0) ? -norm : norm;  /* avoid cancellation */
        double v0 = x[0] - alpha;                     /* v = x - alpha*e1 */
        double vtv = v0 * v0;
        for (int i = 1; i < m - k; i++)
            vtv += x[i] * x[i];
        tau[k] = 2.0 / vtv;

        /* Apply H = I - tau*v*v^T to the trailing columns. */
        for (int j = k; j < n; j++) {
            double *y = &A[j * m + k];
            double s = v0 * y[0];
            for (int i = 1; i < m - k; i++)
                s += x[i] * y[i];
            s *= tau[k];
            y[0] -= s * v0;
            for (int i = 1; i < m - k; i++)
                y[i] -= s * x[i];
        }
        /* Store v (with its implicit 1) and the diagonal entry of R. */
        x[0] = alpha;
        for (int i = 1; i < m - k; i++)
            x[i] /= v0;
    }
}

Why the error does not contain the condition number

The reason Householder QR keeps its orthogonality is local. Each reflector is orthogonal by construction, to working precision, regardless of the matrix it is applied to: the only rounding in forming H is in tau and in the products, and the resulting matrix differs from an exact reflection by a relative amount of order u. Applying a nearly orthogonal matrix to a vector cannot amplify the error already in that vector by more than a factor close to one.

The product of n such reflectors therefore deviates from orthogonality by an amount that accumulates additively, of order n * u, with no factor of kappa. The computed factor R satisfies A + E = Q R with ||E|| of order n * u * ||A||, which is the same backward error statement as for LU with a good pivot sequence. The condition number still governs the forward error of a least-squares solution, but it does not govern whether the computed Q is orthogonal.

That distinction matters when the basis is used for something other than solving a single system. In subspace iteration, in Krylov methods that maintain an explicit basis, and in any algorithm that projects onto the computed range, a basis that has lost orthogonality silently changes the subspace being represented. The projection then converges to something, but not necessarily to the thing that was asked for.

Modified Gram-Schmidt with reorthogonalization is the alternative when the basis vectors are needed explicitly rather than as a product of reflectors. Running the projection twice restores orthogonality to working precision for most matrices, a result usually summarized as “twice is enough”: after the second pass the components along the existing basis are at the level of the rounding error, so a third pass has nothing left to remove. The cost is a second pass over the same data, which doubles the flops of the orthogonalization but keeps the storage pattern of a column-oriented algorithm.

What the extra flops buy

For an m x n matrix with m >= n, Householder QR costs about 2mn^2 - 2n^3/3 flops, compared with about 2mn^2 for either version of Gram-Schmidt. The difference is the 2n^3/3 term, which is the work of applying each reflector to the columns that follow it rather than only to the current one. For a tall matrix, where m is much larger than n, the two are within a small constant of each other; for a square matrix the reflector version does roughly twice the work.

The solve phase is cheap in either case. Given A = QR and a right-hand side b, the least-squares solution is obtained by forming Q^T b, which costs about 2mn flops, discarding the entries below position n, and solving the n x n triangular system R x = c, which costs n^2/2. The triangular solve is where the conditioning of the problem enters, and it is the step that benefits from the accurate factorization.

Storage is the other axis. The reflectors occupy the lower triangle of the m x n array, plus n scalars, so no additional allocation is needed beyond the input. Gram-Schmidt also overwrites the input with Q and needs a separate n x n array for R, which for a tall matrix is smaller than the reflector storage but not free.

The normal equations and the square of the condition number

The alternative to QR for least squares is to form A^T A and solve the normal equations. That halves the flop count and, for a tall matrix, turns a rectangular problem into a small square one. It also squares the condition number: kappa_2(A^T A) = kappa_2(A)^2, so a matrix with kappa = 1e8 produces a normal-equations system with kappa = 1e16, at which point the solution has no correct digits in double precision.

The loss is not a defect of the solver; it is a property of the formulation, and it is the clearest illustration of the split between the condition number of a problem and the backward error of an algorithm. Forming A^T A discards the information in the small singular values by squaring them, and no amount of care in the subsequent factorization recovers it. The same squaring appears when the singular values are computed as the eigenvalues of A^T A, which is why the singular value decomposition is computed by bidiagonalizing A directly rather than by forming the cross-product.

QR avoids the squaring because it never forms the cross-product. The condition number of the least-squares problem is kappa(A) plus a term that depends on the residual, and the QR-based solution achieves an error proportional to that quantity. The caveat is that the relevant condition number for a least-squares problem with a large residual is not kappa(A) alone; when the residual is comparable to ||A|| ||x||, the sensitivity of the solution is governed by the augmented matrix, and a QR of A alone is not sufficient to characterize it. The details belong with the sensitivity analysis in Condition Number, Residual, and Backward Error. What follows from the squaring is that a singular value computation must not go through A^T A at all; the bidiagonalization route that avoids it is described in Lanczos Iteration and Shift-Invert for Extreme Eigenvalues.

When the diagonal of R is not a rank test

Plain Householder QR processes the columns in the given order. If a column is a combination of the ones before it, the corresponding diagonal entry of R is zero in exact arithmetic and tiny in floating point, but a tiny entry is not a reliable rank estimate: it can be small because the matrix is rank deficient, or because the column is nearly dependent, or because the matrix is simply badly scaled, and the three cases call for different responses.

Column-pivoted QR reorders the remaining columns at each step to bring the one with the largest remaining norm to the front. The result is a factorization A P = Q R in which the diagonal of R decreases in magnitude, so the rank can be estimated by counting the entries above a threshold relative to the largest. This is what LAPACK’s geqp3 computes, and it is the standard tool for the rank-revealing problem. The pivoting costs an extra pass over the trailing column norms per step, which is a small fraction of the total for a tall matrix and a noticeable one for a square matrix.

The rank estimate is still a heuristic, because the threshold has to be chosen and because the diagonal entries of a pivoted QR are not guaranteed to be within a bounded factor of the singular values. The guarantee requires a stronger factorization, such as a complete orthogonal decomposition or a rank-revealing QR with a specific column selection rule. For most applications the pivoted factorization with a threshold at a small multiple of u * ||A|| is sufficient, and it has the practical advantage of being a single pass over the matrix.