Numerical Linear Algebra
where the algebra stops being the arithmetic
Orthogonality, measured

Orthogonal is a number

"Q is orthogonal" is a claim about a measurable quantity, ‖QᵀQ − I‖, and on the eight-by-eight Hilbert matrix two standard algorithms return 10⁻¹⁵ and 1 for it. The one that returns 1 still reconstructs the matrix perfectly, which is why nothing warns you.
15 min read 6 figures Measured, not assertedWhere the algebra stops being the arithmeticIdentical algebra, different arithmetic

Orthogonality is the load-bearing property of numerical linear algebra. Almost every stable algorithm in the subject works by applying orthogonal transformations, for the reason set out in the condition number is an amplifier: an orthogonal matrix has κ = 1, so it cannot amplify an error that is already present. QR factorisation, the symmetric eigenvalue problem, the SVD, every least-squares method worth using — all of them are sequences of orthogonal steps.

Which makes the word worth examining. “Q is orthogonal” means QᵀQ = I. In the algebra that is an equation. In the arithmetic it is a measurement, and the measurement has a value.

QᵀQ from classical Gram–Schmidt and from Householder on the 8×8 Hilbert matrixTwo eight-by-eight tables of QᵀQ. The upper one has ones on the diagonal and entries as large as one off it; the lower one is the identity to three decimal places everywhere.A = H8 · κ = 1.5·10¹⁰ · both factorisations reconstruct A to 6·10⁻¹⁷the diagonal is 1 in both — every column is a unit vector either way1.000000000001.000000000001.000000000001.000000000001.00000.002-0.002000001.0000.125-0.13300000.0020.1251.000-1.0000000-0.002-0.133-1.0001.000classical Gram–Schmidt1.000000000001.000000000001.000000000001.000000000001.000000000001.000000000001.000000000001.000Householderclassical ‖QᵀQ − I‖1.4Householder ‖QᵀQ − I‖1.4·10⁻¹⁵largest off-diagonal 1 against 3.1·10⁻¹⁶length is not angle
Fig. 1 QᵀQ for the same eight-by-eight matrix, factorised twice. Below, it is the identity to three decimal places in every entry. Above, the diagonal is still exactly one — every column is still a unit vector — and an off-diagonal entry has reached 1, which means two columns that are supposed to be perpendicular have become the same direction.

What the number means

‖QᵀQ − I‖ collects every deviation from orthogonality into one quantity. The diagonal entries measure how far each column is from having unit length; the off-diagonal entries are the cosines of the angles between pairs of columns, scaled by their lengths.

For Householder QR on any matrix, in double precision, it is around 10⁻¹⁵ — a few units of rounding times a modest factor for the size. That number does not depend on the condition number of the matrix, which is the surprising and important part, and it is the subject of a reflection cannot stop being one.

For classical Gram–Schmidt on the eight-by-eight Hilbert matrix it is 1. Not 10⁻³, not 10⁻¹. One. The largest off-diagonal entry of QᵀQ is 1.0, and two unit vectors whose inner product is 1 are the same vector.

Unit length is not perpendicularity

The upper table in the figure has an exact 1 down its whole diagonal, to twelve digits.

That is not a coincidence and it is the reason this failure goes unnoticed. Gram–Schmidt normalises each column as its last act, dividing by the computed norm, so every column comes out with length 1 to within a rounding whatever else has gone wrong. A check that verified the columns were unit vectors — which is the check somebody writes when they are being careful — would pass.

The information is entirely in the off-diagonal, and getting at it costs a matrix product: form QᵀQ, subtract the identity, take a norm. For a thin Q with n columns that is O(mn²), the same order as the factorisation itself, so it is not something to do inside a production loop. It is exactly something to do in a test.

The factorisation is still correct

Here is the part that makes this genuinely difficult rather than merely subtle.

Both factorisations satisfy ‖A − QR‖/‖A‖ ≈ 10⁻¹⁷. The classical Gram–Schmidt factors multiply back to give the original matrix to full double-precision accuracy. The factorisation is of the right matrix, it is accurate, and its Q is not orthogonal.

So the site’s rule — no decomposition drawn without its residual printed — is necessary and is not sufficient. ‖A − QR‖ certifies that the factors factorise; ‖QᵀQ − I‖ certifies that one of them is what it claims to be. They are different questions and a figure showing a QR carries both.

This is not a hypothetical concern about a quantity nobody uses. Every downstream use of Q assumes QᵀQ = I. Solving least squares by QR computes Qᵀb and then back-substitutes, which is correct only if Q is orthogonal. An orthogonal-basis method for an eigenvalue problem re-orthogonalises precisely because it does not trust the property to survive. Feed a non-orthogonal Q into any of them and the result is wrong in a way that is not attributable to the factorisation, because the factorisation was fine.

How bad it gets, as a function of the matrix

The loss is not a fixed penalty. It scales with the condition number of the matrix being factorised, and it does so at a rate that differs by method.

Loss of orthogonality against condition number, in binary64A log–log plot of the norm of Q-transpose-Q minus the identity against condition number. Classical Gram–Schmidt rises steeply, modified Gram–Schmidt rises gently, and Householder is flat.110²10⁴10⁶10⁸10¹⁰10¹²10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)‖QᵀQ − I‖classicalmodifiedHouseholderκ²uκu8×8, eight seeds per κ, binary64all three reconstruct A
Fig. 2 ‖QᵀQ − I‖ against condition number for three factorisations of the same matrices, eight seeds each. Classical Gram–Schmidt rises like κ², modified like κ, and Householder is flat. The two dashed lines are the predictions κ²u and κu — they were drawn from theory and the measurements landed on them. Drag the significand width and all three curves slide down together, keeping their slopes.

Reading numbers off that: at κ = 10⁵, classical Gram–Schmidt has ‖QᵀQ − I‖ ≈ 10⁻⁸ and modified has 10⁻¹². At κ = 10⁹ classical is at 0.45, which is to say the columns are no longer usefully a basis at all, and modified is at 10⁻⁸. Householder is at 1.8·10⁻¹⁵ at both, and at every other point on the axis.

The κ² and κ behaviours are classical results — Björck’s analysis of modified Gram–Schmidt is from 1967 — and what the figure adds is that they are this code’s behaviour and not a citation. The predictions were drawn before the measurements and the measurements are on them.

What the flatness of the Householder line means

Every other quantity on this site gets worse as the problem gets harder. The Householder line does not, and the reason is structural rather than a matter of degree.

A Householder reflection is I − 2vvᵀ for a unit vector v. Rounding v gives a slightly different unit vector — and I − 2v’v’ᵀ is still exactly a reflection, still exactly orthogonal, for whatever v’ the rounding produced. The rounding perturbs which reflection is in hand; it cannot make the thing in hand stop being one.

Gram–Schmidt has no such protection. Its Q is assembled from subtractions and divisions, and nothing about the assembly enforces the property. When the subtractions go wrong, the result is a matrix, and there is no structure to fall back on.

That difference — between a property enforced by construction and a property achieved by arithmetic — is one of the recurring lessons in numerical analysis and it is worth taking outside the subject. An invariant maintained by a representation survives rounding. An invariant maintained by careful computation does not.

One Householder reflection, and the sign that keeps it stableA vector, the mirror line through the origin, and its reflection landing on the negative first axis at the same distance from the origin.the mirrorx, length 4.000Hx = (-4.000, 0)v = x − αe₁safe sign: α = −‖x‖, so v is formed from a sum and nothing cancelsunsafe sign: α = +‖x‖ gives ‖v‖ only 35.1% of ‖x‖ + ‖x‖ — the digits go‖HᵀH − I‖5·10⁻¹⁶‖Hx‖ − ‖x‖8.9·10⁻¹⁶second component2.2·10⁻¹⁶built from a unit vectororthogonality is structural
Fig. 3 The structural argument, drawn. A reflection is built from a unit vector and a plane; the badge reports that the resulting matrix is orthogonal to 10⁻¹⁶, that it preserves length to 10⁻¹⁶, and that the reflected vector’s second component is zero to 10⁻¹⁶. Perturbing v moves the mirror. It does not stop the mirror being a mirror.

What to do about it

Three options, and the choice is usually easy.

Use Householder. It costs about twice the arithmetic of Gram–Schmidt for a square matrix, it is what LAPACK does, and it removes the problem entirely. For a dense factorisation there is no reason to do anything else.

Use modified Gram–Schmidt. One word different from classical, no extra cost, and the loss becomes κu instead of κ²u — which for κ up to about 10⁷ is entirely acceptable in double. Two Gram–Schmidts is about what that one word is and why it works.

Re-orthogonalise. Run the Gram–Schmidt step twice against the already-computed columns. “Twice is enough” is a real theorem (Kahan and Parlett), and it restores orthogonality to rounding level for any Q that was not catastrophically bad to begin with. It costs a second pass and it is what iterative Krylov methods do, because there Q is built one column at a time and Householder does not fit the access pattern.

What is not an option is more precision. Classical Gram–Schmidt at 53 bits is worse than modified at 24 on a mildly ill-conditioned matrix, which the slider on the figure above shows directly: drag it down and the curves move together, and the classical line at full precision remains above the modified line at any precision on the scale.

Measuring in the right precision

One implementation point that keeps the figures on this site honest, because it is easy to get backwards.

The factorisations in the sliding figure are performed at whatever significand width the reader has selected. The measurement of ‖QᵀQ − I‖ is always performed in double.

Measuring a 24-bit factorisation’s orthogonality error in 24-bit arithmetic would report it as very nearly zero, because the error is exactly the size of what 24-bit arithmetic cannot see. The matrix product QᵀQ would round to the identity and the figure would show a flat line at machine zero for every method — a beautiful, entirely false result, and one that nothing else in the pipeline would flag.

That single decision — the measurement is always in double — is one line in lib/matrix.js and it is the most consequential line in the file.

What is asserted here

Four claims, checked as the figures are generated.

Every classical column is a unit vector, to 10⁻¹², which is what makes the failure quiet.

Some pair of them is not perpendicular: the largest off-diagonal entry of QᵀQ must exceed 0.1. If a change ever made classical Gram–Schmidt behave on the Hilbert matrix, this essay’s central figure would stop being true and the build would say so.

Householder’s off-diagonals are below 10⁻¹³, with that tolerance bracketed between the measured noise — around 10⁻¹⁶ — and 10⁻⁶, which would be a genuine failure.

Both factorisations reconstruct A to 10⁻¹⁴, which is the claim that makes the whole essay interesting rather than merely a report of a bug.

And the orthogonality measurement is fed a counterexample: a shear, which has unit columns and is not orthogonal, and which the check must reject. Without that, a scaling error in the routine would make every method on this site look excellent. Assertions that reject is the thread, and this is the instance that guards the most figures.

Where this goes

Two Gram–Schmidts takes the pair apart line by line, and is the clearest instance on the site of two algorithms that a derivation cannot distinguish. A reflection cannot stop being one is the structural argument for why Householder does not degrade. And the projection and the right angle is what orthogonality is for: the least-squares solution is defined by a perpendicularity condition, and the accuracy with which that condition can be enforced is exactly the accuracy of the Q used to enforce it.

The Gram matrix as a diagnostic

QᵀQ has a name outside this context: it is the Gram matrix of the columns, the table of all pairwise inner products. It is worth recognising, because the same object answers several questions that are usually treated separately.

Its diagonal entries are the squared lengths. Its off-diagonal entries, divided by the lengths, are the cosines of the angles between columns. Its determinant is the squared volume of the parallelepiped the columns span, so a Gram determinant near zero means the columns are nearly dependent. Its eigenvalues are the squares of the singular values of Q.

That last identity is why forming a Gram matrix squares a condition number, and it connects this essay directly to the road that squares the problem: the normal equations are the Gram matrix of A, and the reason they lose accuracy is the reason QᵀQ is a sensitive thing to compute.

For a Q that is supposed to be orthonormal, the Gram matrix should be the identity, and the extent to which it is not is the single number this whole field is about. Printing it, as the figure at the top of this essay does, is more informative than printing the norm — the norm says how bad, the table says which pair.

What a bad Q does to a Krylov method

The most common place loss of orthogonality causes real damage is not in a dense factorisation, where Householder is available, but in iterative methods, where it is not.

Arnoldi and Lanczos build an orthonormal basis for a Krylov subspace one vector at a time, each new vector obtained by applying A to the previous one and orthogonalising against everything so far. The access pattern is inherently Gram–Schmidt-shaped: the vectors do not exist in advance, so there is no matrix to reflect.

When orthogonality is lost, two things happen and both are famous. In Lanczos, eigenvalues appear multiple times in the computed spectrum — “ghost” eigenvalues, which are not eigenvalues of anything, and which arise because the basis has silently started to repeat directions. In GMRES, the computed residual stops tracking the true residual, so the method reports convergence it has not achieved.

The standard defences are re-orthogonalisation, which is expensive but reliable, and selective re-orthogonalisation, which is a genuinely subtle piece of engineering. Both exist because the property this essay measures does not survive on its own.

How far a perturbation of size ε moves an eigenvalue, 8×8A log–log plot of eigenvalue movement against perturbation size. The symmetric case lies on a line of slope one; the non-symmetric case lies on a line of slope one eighth, and at a perturbation of ten to the minus sixteen it has already moved by a hundredth.10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹size of the perturbation ‖δA‖how far the eigenvalues moveJordan block, ε^(1/8)symmetric, ≤ ‖δA‖rounding error alone moves it to 10⁻²six seeds per symmetric point; Jordan is closed formsymmetry beats precision
Fig. 4 Why the eigenvalue case is where it hurts most. A basis that has lost orthogonality is a perturbation of the problem, and for a non-symmetric matrix a perturbation of size ε moves the spectrum by a fractional power of ε. A loss of 10⁻⁸ in the basis can move eigenvalues by 10⁻¹.

Orthogonality in more than one norm

A brief technical point that the figures here elide, since it occasionally matters.

‖QᵀQ − I‖ can be measured in the Frobenius norm, which is the root sum of squares of all entries, or in the 2-norm, which is the largest singular value. This site uses Frobenius throughout, because it is cheap and because it aggregates every deviation rather than only the worst direction.

The two differ by at most a factor of √n, so for the matrix sizes here they tell the same story. The 2-norm is the right choice when the question is about a specific direction — “how badly does this basis fail for one particular vector” — and Frobenius is the right choice for the question this site asks, which is “how far is this from being an orthonormal basis at all”.

Stating which norm is being used is not pedantry. A great many quoted error figures in this subject are incomparable because their norms differ, and the difference can be a factor of the matrix dimension.

Fitting the same degree-11 polynomial in two basesOn the left, the data and two fitted curves that lie on top of each other. On the right, the condition numbers of the two design matrices, ten orders of magnitude apart.00.10.20.30.40.50.60.70.80.910.20.40.60.81xyboth fits, drawn on top of one anothercondition number of the design matrixmonomial1.2·10⁸Chebyshev2.5largest fitted coefficientmonomial113Chebyshev0.51rms residual: 3.9·10⁻⁶ and 2.9·10⁻⁷ — the data is fitted either way.30 points, degree 11, single precisionthe basis is part of the problem
Fig. 5 And the reason the choice of basis matters before any factorisation happens. The same fit in two bases: monomials, whose columns are nearly parallel and whose design matrix has κ = 1.2·10⁸, and Chebyshev polynomials, which are nearly orthogonal by design and give κ = 2.5. Orthogonality is something a basis can be chosen to have, not only something an algorithm can try to achieve.

The habit worth keeping

Any time code produces a matrix that is supposed to be orthogonal — a QR factor, an eigenvector matrix, a rotation, a basis — compute ‖QᵀQ − I‖ once and look at it.

It costs a matrix product, it is a single number, and it distinguishes the working case from every common way of getting it wrong. Most of the time the answer is 10⁻¹⁵ and the check was unnecessary. The cases where it is not are the cases where every downstream result is wrong, and nothing else would have said so, because the columns will still be unit vectors and the factorisation will still reconstruct.

That habit, generalised, is the site’s rule. No decomposition without its residual is the same instinct applied to a different property, and between them the two numbers cover most of what can go wrong with a factorisation before anything downstream sees it.

Growth factor under partial pivoting: the bound, the worst case, and realityGrowth factor against matrix size on a logarithmic vertical axis. The two-to-the-n bound rises as a straight line; Wilkinson's matrix sits exactly on it; random matrices stay near one.0816243240110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴matrix size ngrowth factor max|u| / max|a|the 2ⁿ⁻¹ boundworst of 30 randommedian randomWilkinson's matrix sits on the bound30 Gaussian matrices per sizeat n = 40: bound 5.5·10¹¹, worst 4.8
Fig. 6 The comparison that puts this field in context. Elimination’s stability is empirical — a measured quantity that stays near three against a bound of 5.5·10¹¹, as the bound that is never attained sets out. Orthogonality under Householder is structural, and needs no such argument.

Where this goes next

for each previous column i, subtract the projection of column j onto qᵢclassicalr[i][j] = qᵢ · a[j] ↑ the ORIGINAL columnv = v − r[i][j] · qᵢthe three worst |qᵢ · qⱼ|:columns 7 and 8: 1columns 6 and 8: 0.13columns 6 and 7: 0.13modifiedr[i][j] = qᵢ · v ↑ what is LEFT of itv = v − r[i][j] · qᵢthe three worst |qᵢ · qⱼ|:columns 1 and 8: 4.4·10⁻⁷columns 2 and 8: 2.7·10⁻⁷columns 3 and 8: 2.4·10⁻⁸The two R factors agree to 1.2·10⁻⁶ relative. The two Q factors do not.the 8×8 Hilbert matrixone word, eight orders Two Gram–Schmidts One argument changes. Classical Gram–Schmidt projects the original column onto each previous direction; modified projects what is left of it. In exact arithmetic the coefficients are identical. In floating point they differ by eight orders of magnitude in the thing that matters. the mirrorx, length 4.000Hx = (-4.000, 0)v = x − αe₁safe sign: α = −‖x‖, so v is formed from a sum and nothing cancelsunsafe sign: α = +‖x‖ gives ‖v‖ only 35.1% of ‖x‖ + ‖x‖ — the digits go‖HᵀH − I‖5·10⁻¹⁶‖Hx‖ − ‖x‖8.9·10⁻¹⁶second component2.2·10⁻¹⁶built from a unit vectororthogonality is structural A reflection cannot stop being one Householder QR holds orthogonality at 10⁻¹⁵ whatever the condition number of the matrix, and Gram–Schmidt does not. The reason is not that it is more careful. It is that its Q is built from unit vectors, and rounding a unit vector gives a different reflection rather than a broken one. everything Ax can reachb = (1.1, 0.4, 1.5)Ax, the closest reachable pointr = b − Ax‖Aᵀr‖ / (‖A‖‖r‖)1.7·10⁻¹⁶‖b‖² − ‖Ax‖² − ‖r‖²1.3·10⁻¹⁵‖r‖1.3200 random nearby points of the plane were tried; none is closer.a 3×2 system, Householder QRperpendicularity is checked The projection and the right angle The least-squares solution is the one whose residual is perpendicular to everything the columns can reach. That is not a mnemonic — it is an equation, Aᵀr = 0, and the computed answer satisfies it to 10⁻¹⁶.