Numerical Linear Algebra
where the algebra stops being the arithmetic
Concept

Gram–Schmidt — where it appears

Orthogonalising a set of vectors one at a time by subtracting their projections onto the ones already done. The order of the subtractions is algebraically irrelevant and numerically decisive: the classical and modified variants differ by fifteen decades of orthogonality.

Named by 11 essays across 4 fields — each of them below, with the objects they name alongside it.

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.

orthogonality · Orthogonality 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.

orthogonality · Gram–Schmidt 05101520253010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹iteration‖RᵀR − I‖ and ‖r‖/‖b‖step n‖RᵀR − I‖residual30×30, run for exactly n stepsexact arithmetic would end here

An orthogonalisation nobody calls one

Conjugate gradients are derived as a minimisation and behave as an orthogonalisation, which is why the finite-termination property in every textbook is not a property the method has in floating point.

iterative · Krylov 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.

orthogonality · Householder 10¹10²10³10⁴10⁵10⁶10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹number of terms addedrelative error against the exact sumin orderin a treecompensatedbinary32 · terms are 1/icompensated: 3·10⁻⁸

The order they are added in

Addition is associative in the algebra and is not associative in the arithmetic. The same million numbers, added in a different order, give answers that differ in the third significant figure — and the fix is not a wider float, it is a different order.

arithmetic · Summation classical Gram–Schmidt4.62·10⁻¹⁰modified Gram–Schmidt1.49·10⁻¹²Householder, one sweep2.03·10⁻¹⁴reduction tree, 16 leaves1.48·10⁻¹⁵departure from orthogonality, logarithmicthe tree, at four depths‖AᵀA − RᵀR‖/‖AᵀA‖, depth 13.4·10⁻¹⁵‖AᵀA − RᵀR‖/‖AᵀA‖, depth 24.3·10⁻¹⁵‖AᵀA − RᵀR‖/‖AᵀA‖, depth 31.7·10⁻¹⁵‖AᵀA − RᵀR‖/‖AᵀA‖, depth 41.5·10⁻¹⁵the same algebra, four timestwo of them are products of reflections

A reduction that changes the order

A tall-skinny QR computed as a tree of independent block factorisations touches a 512×12 matrix once instead of twelve times, computes a completely different sequence of roundings from the sweep it replaces, and returns ‖AᵀA − RᵀR‖/‖AᵀA‖ = 1.65·10⁻¹⁵ against the sweep's 9.95·10⁻¹⁵. On the same matrix classical Gram–Schmidt returns 4.6·10⁻¹⁰.

cost · Communication how far is it from A to an orthogonal matrix?smaller is nearer · the polar factor minimises this in every unitarily invariant normpolar factor U1.8554QR, signs fixed2.1265QR as returned3.8226200 drawn at randomκ = 10polar factor1.9QR, signs fixed2.1QR as returned3.8best of 200 random2.7‖A − QR‖ is the same either wayand ‖A − Q‖ is not

The nearest orthogonal matrix

Every field that has to clean up a drifted rotation reaches for QR, and QR does not answer the question. The nearest orthogonal matrix is the orthogonal factor of the polar decomposition — nearer by about a tenth, and, more to the point, the same matrix whatever order the columns were written in. QR's answer changes completely.

orthogonality · Polar decomposition polar: U₁P − U₂QR: Q₁P − Q₂-3.3·10⁻¹⁶-3.3·10⁻¹⁶2.2·10⁻¹⁶-1.9·10⁻¹⁶5.6·10⁻¹⁷-1.7·10⁻¹⁶6.7·10⁻¹⁶-8.9·10⁻¹⁶2.8·10⁻¹⁶-5.6·10⁻¹⁶-3.6·10⁻¹⁶6.1·10⁻¹⁶-8.3·10⁻¹⁶8.3·10⁻¹⁶-7.2·10⁻¹⁶-10·10⁻¹⁶2.2·10⁻¹⁶4.4·10⁻¹⁶1.7·10⁻¹⁶-4.4·10⁻¹⁶-5.6·10⁻¹⁷5.6·10⁻¹⁷0-7.2·10⁻¹⁶-4.4·10⁻¹⁶3.9·10⁻¹⁶5.6·10⁻¹⁷2.8·10⁻¹⁶-1.1·10⁻¹⁶-3.3·10⁻¹⁶2.2·10⁻¹⁶1.1·10⁻¹⁶-2.8·10⁻¹⁶5.6·10⁻¹⁶-2.8·10⁻¹⁶5.6·10⁻¹⁶0.67-0.35-0.23-0.620.19-0.72-0.770.170.20.290.190.93-0.180.55-0.0290.06-0.460.11-0.17-0.61-0.230.38-0.55-0.380.750.870.661.1-0.36-0.530.42-0.670.180.180.52-0.3Frobenius norms, columns reordered‖U₁P − U₂‖2.8·10⁻¹⁵‖Q₁P − Q₂‖3‖Q‖, for scale2.4κ of the matrix100Frobenius distance between the two answers, on one scale0 to 4polar2.8·10⁻¹⁵QR3.048‖Q‖ = 2.449the column space did not moveand one of the two answers did

A test with no answer in it

A caller with no reference answer can still ask whether a routine answered the right question: reverse the columns, run it again, compare. The polar factor's two answers agree to 10⁻¹⁵ at every conditioning drawn; a QR's differ by 2.353 on matrices whose own norm is 2.449. The test has a floor, and the floor is measurable too.

orthogonality · Polar decomposition 10¹10³10⁵10⁷10⁹10¹¹10¹³10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)‖QᵀQ − I‖classical once, Householdermodified once, Householderclassical twice, Householderclassical twice, Cholesky QRκ²uκuat κ = 10⁸classical once, Householder0.0042modified once, Householder5.7·10⁻⁹classical twice, Householder3.1·10⁻¹⁵classical twice, Cholesky QR10⁻¹⁵64×16 in blocks of 4, three seedsHouseholder inside does not help between

A stable block is not a stable basis

Block Gram–Schmidt orthogonalises twice over — between blocks, and inside each one. Householder inside the blocks does not stop the classical between-block step losing orthogonality like κ², 4.2·10⁻³ at κ = 4.3·10⁷, and a second pass does not stop Cholesky QR inside the blocks breaking down at κ = 10⁸. Each level fails only on ill-conditioning placed at its own level, and one variant holds 3·10⁻¹⁵ on every placement.

orthogonality · Gram–Schmidt 10²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)‖x̂ − x‖ / ‖x‖modified, through Qᵀbclassical, either route‖QᵀQ − I‖, modifiedHouseholdermodified, on [A b]κ²uκuforward error at κ = 10⁸Householder4·10⁻⁹modified, through Qᵀb0.13modified, on [A b]2.7·10⁻¹⁰classical840×8, three seeds, against an exact rational solveQ is not orthogonal; x is right

The right-hand side as one more column

Modified Gram–Schmidt's Q is 4.3·10⁻⁹ from orthogonal at κ = 10⁸, and a least-squares solve that multiplies b by it is wrong by 0.13. Hand the same routine b as an extra column instead and the answer is right to 2.7·10⁻¹⁰ — closer than Householder's 4.0·10⁻⁹. Classical Gram–Schmidt gains nothing from the same trick, to the last bit.

orthogonality · Gram–Schmidt three seeds, medianblock MGS, Qᵀb, κ 10⁸0.23block MGS, b as a block1.1·10⁻⁹Householder2.8·10⁻⁹10²10⁴10⁶10⁸10¹⁰10¹²10⁻¹⁶10⁻¹²10⁻⁸10⁻⁴110⁴10⁸κ(A)relative error of xblock MGS, Qᵀbblock CGS, either routeblock MGS, b as a blockHouseholderthe same Q in both block MGS routesonly the order in which b meets it differs

What the appended block inherits

Modified Gram–Schmidt on [A b] solves least squares as well as Householder, although its Q is not orthogonal. A block code appends b as one more block. Block modified Gram–Schmidt inherits the rescue at every placement of the ill-conditioning: at κ = 10⁸ the appended block gives 6.9·10⁻¹⁰ where the same Q through Qᵀb gives 8.9·10⁻³. Block classical Gram–Schmidt gets the same wrong answer both ways, to the last bit. And the variant whose Q is orthogonal to 10⁻¹⁵ — two passes with Cholesky QR inside — is a hundred thousand times worse than Householder when the ill-conditioning is inside the blocks, because its R is wrong.

orthogonality · Gram–Schmidt

Named alongside it

The objects these essays reach for when they reach for this one.

OrthogonalityHouseholder reflectionLoss of orthogonalityQR factorisationCondition numberBackward stabilityReorthogonalisationBlock methodsCholesky QRExact ground truthGaussian eliminationHouseholder

All concepts