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

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.
16 min read 6 figures Where the algebra stops being the arithmeticMeasured, not assertedIdentical algebra, different arithmetic

Look again at the flat line.

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. 1 Loss of orthogonality against condition number, eleven decades, eight seeded matrices at each point. Two of the three methods rise. The third does not move: 1.8·10⁻¹⁵ at κ = 10, and 1.8·10⁻¹⁵ at κ = 10¹¹. Nothing else on this site behaves that way.

Every other quantity here degrades as the problem gets harder. Forward errors rise with κ, loss of orthogonality in Gram–Schmidt rises with κ or κ², approximation errors rise, rank decisions get harder. The Householder line is flat, and the flatness is not a matter of the method being better by a large factor. It is a different kind of thing, and the reason is worth understanding because it generalises far beyond this algorithm.

The construction

Each step of Householder QR takes a column x and reflects it onto a multiple of the first coordinate axis. The reflector is

H = I − 2vvᵀ, where v = (x − αe₁) / ‖x − αe₁‖ and α = ±‖x‖.

v is a unit vector by construction — it is divided by its own norm — and H is the matrix that reflects space in the hyperplane perpendicular to it.

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. 2 One reflection, drawn. The vector, the mirror through the origin, and the image landing on the negative first axis at the same distance from the origin. The badge reports three measurements: the reflection is orthogonal to 10⁻¹⁶, it preserves length to 10⁻¹⁶, and the second component of the reflected vector is zero to 10⁻¹⁶.

The structural argument

Now the point.

Suppose the arithmetic rounds while computing v. What comes out is v̂, a slightly different vector, and — after the normalisation, which is also rounded — a vector of length 1 ± u.

The matrix I − 2v̂v̂ᵀ is exactly a reflection, for whatever v̂ the rounding produced. It reflects in the hyperplane perpendicular to v̂ rather than the one perpendicular to v, which is a slightly different plane, but it is still a reflection: still orthogonal, still norm-preserving, still an involution.

The rounding perturbed which orthogonal transformation was applied. It did not perturb the orthogonality, because orthogonality is not something the arithmetic computes — it is a property of the form I − 2vvᵀ, and that form survives having its v rounded.

Compare Gram–Schmidt. Its Q is assembled by subtraction and division, entry by entry, and nothing about the assembly enforces the property. When the subtractions are inaccurate, the result is a matrix, and there is no structure to fall back on. Its orthogonality is achieved, and it can therefore fail to be achieved.

Achieved against enforced

That distinction is the transferable lesson, and it is worth stating outside the subject.

An invariant maintained by a representation survives perturbation. An invariant maintained by computation does not.

Store a rotation as an angle and it is a rotation no matter what the arithmetic does to the angle. Store it as a 2×2 matrix and repeated multiplication drifts it away from being orthogonal, which is why graphics and robotics code re-orthonormalises rotation matrices periodically, and why quaternions are popular: a normalised quaternion is a rotation by construction.

The same idea explains why the algorithms on this site are made of orthogonal steps rather than inverses. An orthogonal step has κ = 1 and therefore amplifies nothing, and it cannot cease to be orthogonal. Two properties, both structural, and together they are most of what makes an algorithm stable.

The sign, which is the one place it can go wrong

There are two choices of α — plus or minus ‖x‖ — and they differ only in which of the two directions along the axis x is reflected onto. Mathematically it does not matter.

Numerically one of them is a cancellation. If x already points mostly along e₁ with x₁ > 0, choosing α = +‖x‖ makes the first component of v equal to x₁ − ‖x‖, a difference of two nearly equal positive numbers. The direction that defines the reflection is then computed from cancelled digits, and every operation afterwards inherits it — cancellation takes the answer is the mechanism.

Choosing α with the opposite sign to x₁ makes it a sum. Nothing cancels. The figure at the top of this essay reports what the unsafe choice would cost on that particular vector: v would come out at a small fraction of the size of its operands, and the digits would be gone.

The fix is one sign flip and it is in every serious implementation. It is in almost no derivation, because the derivation is correct either way, which is the recurring shape of identical algebra, different arithmetic.

What Householder costs

For an m×n matrix with m ≥ n, Householder QR costs about 2mn² − 2n³/3 operations, which for a square matrix is 4n³/3 — roughly twice Gaussian elimination and about twice Gram–Schmidt.

That factor of two is the price of orthogonality that does not depend on the matrix, and for a dense factorisation it is worth paying without thinking about it. Two situations complicate the choice.

Q is often not needed explicitly. The reflectors are stored as their v vectors, and applying Q or Qᵀ to something costs O(mn) per vector without ever forming the matrix. LAPACK returns the reflectors, not Q, for exactly this reason. If Q is formed explicitly the extra cost is real and frequently unnecessary.

Column-at-a-time access. Krylov methods generate their basis vectors one at a time and cannot factorise a matrix they do not yet have. Householder does not fit that pattern, which is why Gram–Schmidt with re-orthogonalisation survives there.

For everything else — a least-squares fit, a rank decision, an explicit orthonormal basis — Householder is the answer and the choice is not interesting.

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. 3 The two outcomes on the same matrix. Below, Householder’s QᵀQ: the identity to three decimals in every entry. Above, classical Gram–Schmidt’s: unit diagonal, and an off-diagonal entry that has reached one. The factorisations reconstruct the matrix equally well.

Givens rotations, briefly

The other orthogonal building block deserves a mention because it shares the structural property and is used where reflections are clumsy.

A Givens rotation acts on two coordinates at a time, zeroing one entry with a plane rotation built from a cosine and a sine satisfying c² + s² = 1. Like a reflection, it is orthogonal by construction: round c and s, renormalise, and what remains is still exactly a rotation of a slightly different angle.

Rotations cost more than reflections for a dense factorisation — they zero one entry at a time rather than a whole column — and they win where the matrix is sparse or already nearly triangular, because they touch only two rows. Updating a QR factorisation when a row is added, which is the standard move in recursive least squares, is a sequence of Givens rotations.

Both are instances of the same idea, and the idea is the reason this field exists as a separate part of the site: orthogonal transformations are the only ones that can be applied repeatedly without watching the error.

What is asserted here

The figure’s claims are checked as it is drawn, and each is a measurement with a tolerance bracketed rather than a round number chosen.

The reflection preserves length: ‖Hx‖ equals ‖x‖ to 10⁻¹³, checked as two routes to the same number rather than as an inequality. It lands the vector on the axis: the second component of the image is below 10⁻¹³. It is orthogonal: ‖HᵀH − I‖ below 10⁻¹⁵.

And it is an involution: HH must be the identity to 10⁻¹⁴. That is the second route, and it is the one that would catch a sign error in the construction of the reflector — a wrong sign gives a matrix that is still orthogonal and still norm-preserving and is not a reflection, so the first three checks would all pass.

The cancellation claim is a measurement too: the unsafe sign is computed, and the size of the resulting v relative to what it would be with the safe sign is required to show real cancellation. The essay does not warn about it; the build measures it.

Across the loss curve, Householder’s ‖QᵀQ − I‖ is required to stay below a bound proportional to the unit roundoff and the matrix size at every condition number, which is the flat line as an assertion. If it ever rose with κ — which is what a subtle bug in the reflector application would cause — the build would stop, and the essay’s central claim would be the thing that caught it.

The cost accounting, in full

It is worth being precise about the arithmetic, because “Householder costs twice as much” is repeated without the qualifications that make it decision-relevant.

For an m×n matrix with m ≥ n, factorising costs about 2mn² − 2n³/3 flops. Classical or modified Gram–Schmidt costs about 2mn². For a square matrix that is 4n³/3 against 2n³, so Householder is cheaper — the “twice as expensive” folklore comes from the case where Q is formed explicitly, which costs a further 2mn² − 2n³/3 and is usually unnecessary.

The reflectors are stored as their v vectors, in the space the zeroed entries of A vacated, and applying Q or Qᵀ to a vector costs O(mn) by applying the reflectors in sequence. LAPACK returns them in that form for exactly this reason, and code that calls orgqr to materialise Q is frequently spending more than the factorisation cost on something it will only use to multiply.

So the honest summary: Householder is not more expensive than Gram–Schmidt for a dense factorisation, it is more expensive than a Gram–Schmidt that also produces Q explicitly, and it is the wrong shape for algorithms that generate their columns one at a time.

Why the same argument applies to rotations

The structural property is not special to reflections, and the second orthogonal building block shows why the argument is about form rather than about a particular matrix.

A Givens rotation is built from a cosine and a sine with c² + s² = 1. Round c and s, renormalise, and what remains is exactly a rotation through a slightly different angle. Same guarantee, same reason: the property is enforced by the parameterisation rather than achieved by the arithmetic.

Rotations zero one entry at a time rather than a whole column, which makes them more expensive for a dense factorisation and cheaper for a sparse or nearly-triangular one, since they touch only two rows. That is why the QR algorithm for eigenvalues uses rotations on a Hessenberg matrix, and why updating a factorisation after a new row arrives — the standard move in recursive least squares — is a sequence of Givens rotations rather than a refactorisation.

The Jacobi methods used throughout this site are the same object again: each sweep is a sequence of plane rotations, so the accumulated eigenvector matrix cannot drift away from orthogonality. Symmetry is worth more than precision reports ‖VᵀV − I‖ = 10⁻¹⁵ for the Hilbert matrix’s eigenvectors, and that number is flat for the same structural reason the Householder line is flat.

Jacobi sweeps on the 6×6 Hilbert matrixA plot of the off-diagonal norm against sweep number, falling from about one to ten to the minus seventeen in five sweeps, with the matrix shown at three stages beneath it.01234567810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹sweep‖off-diagonal‖eigenvalues, largest first1.6190.24240.016326.157·10⁻⁴1.257·10⁻⁵1.083·10⁻⁷κ = 1.5·10⁷‖VᵀV − I‖ = 1.9·10⁻¹⁵1.000.500.330.250.200.170.500.330.250.200.170.140.330.250.200.170.140.130.250.200.170.140.130.110.200.170.140.130.110.100.170.140.130.110.100.09as given1.620.050.00-0.040.00-0.000.050.240.000.010.00-0.000.000.000.000.00-0.00-0.00-0.040.010.000.020.000.000.000.00-0.000.000.00·-0.00-0.00-0.000.00·0.00after one sweep1.62······0.24······0.00······0.02······0.00······0.00after 8 sweepseach sweep is n(n−1)/2 rotationsthe trace is conserved
Fig. 4 Rotations at work. Jacobi on the six-by-six Hilbert matrix reaches an off-diagonal norm below 10⁻¹⁷ in four sweeps, and the eigenvector matrix it accumulates is orthonormal to 10⁻¹⁵ — because every one of the sixty rotations it performed was exactly a rotation, whatever its cosine rounded to.

What “structural” buys elsewhere

The general lesson is worth taking out of the subject, because it applies wherever a computation must preserve a property.

Represent a rotation as three angles and it is a rotation no matter what happens to the angles. Represent it as a 3×3 matrix and repeated composition drifts, which is why graphics and robotics code re-orthonormalises periodically and why unit quaternions are popular — a normalised quaternion is a rotation by construction, and normalising is one division.

Represent a probability distribution as unnormalised weights plus a normalisation step and it sums to one exactly, every time. Represent it as numbers that are supposed to sum to one and it drifts.

Represent a symmetric matrix by its upper triangle and it is symmetric. Store all n² entries and compute with them, and the asymmetry that rounding introduces is enough to make a symmetric eigensolver refuse — which is exactly the precondition this site’s Jacobi routine enforces, and exactly why (A + Aᵀ)/2 appears in careful code.

In each case the pattern is the same: choose a representation in which the invariant cannot be violated, and the arithmetic can be as sloppy as it likes.

Backward and forward error against the condition numberA log–log plot over twelve decades of condition number. The backward error is a flat line at ten to the minus sixteen; the forward error rises in proportion to the condition number.110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)relative errorforward errorbackward errorpredicted: κ · u8×8, 20 seeds per κ; dashed is the worstthe problem worsens, not the method
Fig. 5 And the limit of what structure buys. Orthogonal steps keep the algorithm’s contribution flat — that is the lower line, unmoved across twelve orders of magnitude of conditioning. They do nothing about the upper line, because that one belongs to the problem. No amount of good construction makes an ill-conditioned question well posed.

Where this leaves the field

Three essays on orthogonality, and a single practical conclusion.

Orthogonal is a number established that the property is a measurement, that it can be badly violated by a factorisation that reconstructs its matrix perfectly, and that unit columns are not perpendicular columns. Two Gram–Schmidts established that one word of difference produces eight orders of magnitude, and that the R factors agree while the Q factors do not. This essay establishes why one method escapes the problem entirely.

The conclusion: use Householder for dense factorisations, use modified Gram–Schmidt with re-orthogonalisation where the access pattern requires it, and measure ‖QᵀQ − I‖ either way. The first two are free. The third costs a matrix product and is the only thing that would reveal whether either had gone wrong.

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 contrast worth ending on. Gaussian elimination’s stability rests on a quantity that is measured to be small and bounded only by something useless — see the bound that is never attained. Householder’s rests on a property that cannot be violated. Both algorithms work; only one of them has a theorem behind it, and it is the more expensive one.

Which is the honest summary of the trade. Elimination is a sequence of choices is faster and relies on an observation; orthogonal is a number measures what the alternative buys; and the projection and the right angle is where the guarantee is actually spent.

What this does not protect against

Being precise about the guarantee, because “orthogonality is structural” can be over-read.

Householder QR is backward stable: the computed Q and R satisfy QR = A + E with ‖E‖ at the level of rounding, and the computed Q is within rounding of an exactly orthogonal matrix. Both of those are theorems.

Neither says the factorisation is accurate in the sense of being close to the exact QR of A. For an ill-conditioned matrix the exact R and the computed R can differ substantially, because R inherits the conditioning of A — and it must, since R and A have the same singular values. What is guaranteed is that the computed pair factorises something very close to A, which is the strongest statement available and is not the same statement.

The distinction is exactly the one the exact answer to a nearby problem sets up, applied to a factorisation rather than a solution. The structural argument in this essay is about one factor’s orthogonality; the conditioning of the problem still governs how much the factors themselves can be trusted, and orthogonal steps guarantee only that the algorithm adds nothing to it.

That is worth saying because the flat line in the first figure invites the opposite conclusion. It says the method does not degrade. It does not say the answer improves.

Where this goes next

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 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. 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. 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⁻¹⁶.