Numerical Linear Algebra
where the algebra stops being the arithmetic
Concept

Householder reflection — where it appears

A reflection that maps a vector onto a coordinate axis, and the building block of the most stable way to compute a QR factorisation. It is orthogonal to the last bit whatever it is applied to, which is why a factorisation built from reflections has a residual at the working precision by construction.

Named by 25 essays across 5 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 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 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 rounds on the critical pathHouseholder sweep48reduction tree4Cholesky QR4words sentHouseholder sweep1170reduction tree1170Cholesky QR2160two counts, two rankingsrounds, sweep ÷ tree12words, Cholesky ÷ tree1.8arithmetic, tree ÷ sweep1.5the rounds separate the threeand the words do not

The message and the word

Three factorisations of one matrix on sixteen processors: 48 communication rounds, 4, and 4. The words sent are 1,170, 1,170 and 2,160 — so the method with the fewest rounds sends the most words, and the count that separates the three is the one no operation count can see.

cost · Communication -0.16-0.63-0.493.43-0.025-0.632.3-0.69-0.95-1.70.076-0.49-0.694.3-1.6-1.41.93.4-0.95-1.65.40.0141.63-1.7-1.40.0144.60.055-0.0250.0761.91.60.0553.1A, symmetric→-0.164.600004.65.92.500002.51-1.20000-1.24.90.3800000.385.80.2400000.242H = QᵀAQ, tridiagonal‖A − QHQᵀ‖/‖A‖1.1·10⁻¹⁵below the subdiagonal0worst eigenvalue movement7.1·10⁻¹⁵a similarity, so the spectrum is untouched — and every later step is O(n²) rather than O(n³)one reduction, then every iteration is cheapthe eigenvalues did not move

The form that makes it affordable

One Householder reduction, done once, turns every subsequent iteration of the eigenvalue algorithm from cubic to quadratic cost. It changes no answer at all, which is why it is easy to describe as an optimisation and wrong to.

spectra · The QR algorithm 0.46-1.11.5-0.470.2-0.181.90.690.42-0.160.028-0.0160-1.50.160.880.099-0.1100.680.762.30.930.0480-0.250.110.952.11.500001.52.3after reflector 2the highlighted entries are the bulge — the only thing that is not Hessenbergentries below the subdiagonal3reflectors used2reflectors in a whole step5the shifts are never formed — only their sum and productand both of those are real

Two shifts that are never formed

The double shift is defined as a factorisation of (A − μI)(A − μ̄I), which nobody computes. What is computed is the first column of that product — three numbers — and the bulge those three numbers create, pushed down the subdiagonal by n − 2 reflectors until it falls off the bottom.

spectra · Francis 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 1357911131510⁻¹⁶10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1k, where h = 1 − 10⁻ᵏrelative error of 1 − h1 − ‖R⁻ᵀa‖²1 − ‖row of Q₁‖²‖row of Q₂‖²40 × 6k = 12, Cholesky route1.1·10⁻⁴k = 12, thin QR route4.4·10⁻⁴complement, worst k3·10⁻¹⁵one minus a sum of squaresor the sum of the other squares

One minus a leverage is a subtraction

Every deletion diagnostic divides by 1 − h, and computing it as one minus a computed leverage loses digits in proportion to 1/(1 − h), however accurate the leverage. The complementary block of a QR factor gives the same number as a sum of squares and loses κ(A)·u instead: every digit on a well-conditioned design, and half the digits the subtraction loses on a design whose far point is what made 1 − h small.

leastsquares · Leverage 8×8, κ = 10⁸direction at ε = 10⁻², ‖QᵀQ − I‖1.5·10⁻¹⁵scalar at ε = 10⁻², ‖QᵀQ − I‖0.065direction at ε = 10⁻², residual0.00510⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²relative perturbation εdeparture‖QᵀQ − I‖, scalar‖A − QR‖/‖A‖, scalar‖A − QR‖/‖A‖, direction‖QᵀQ − I‖, directionflat line: a different reflection is still a reflectionsloped line: a non-reflection is not

One number that has to be right

Householder's orthogonality was called structural: a reflection is built from a unit vector, so rounding the vector names a different reflection rather than a broken one. Tested by breaking it, the claim is narrower and sharper. Perturb every component of the reflector by a relative 10⁻², and ‖QᵀQ − I‖ stays at 1.5·10⁻¹⁵ while the factorisation moves to 5·10⁻³. Perturb the one stored scalar by the same amount and ‖QᵀQ − I‖ is 6.5·10⁻². The structure is one degree of freedom, and the departure is four times its relative error.

orthogonality · Householder what a block storesnumbers in T, block of 11numbers in T, block of 16136‖QᵀQ − I‖ there3.9·10⁻¹⁵10⁻¹⁶10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶block size‖QᵀQ − I‖124816T perturbed by 10⁻⁶T as computed1, 3, 10, 36, 136 computed numbers in the triangleand a factor of five in what they produce

A triangle where the scalar was

Every level-3 QR assembles a block of reflectors into Q = I − Y T Yᵀ, and T is computed by a recurrence whose inputs are its own previous columns. A block of sixteen carries 136 computed numbers where sixteen separate reflections carry sixteen. The orthogonality it produces is 3.9·10⁻¹⁵ against the single reflector's 7.8·10⁻¹⁶ — a factor of five for a hundred and thirty-six times as many things that have to be right.

orthogonality · Householder 81930415263110²10⁴10⁶10⁸10¹⁰10¹²size of the matrix|rₙₙ| ÷ σₘᵢₙthe two agreec = 0.2c = 0.35c = 0.5no ceilingratio at n = 1021ratio at n = 646.8·10¹⁰interchanges, anywhere0the greedy rule never had a choiceand the gap grows with every row

The cheap rank and what it cannot see

Almost nobody computes singular values to decide a rank. The standard substitute is QR with column pivoting, read off the diagonal of R — and there is a triangular matrix on which the greedy rule makes no interchange at all, has no better column available at any step, and reports a matrix eight orders of magnitude further from singular than it is.

spectra · Rank 110¹-101pieces the inner product was summed inqᵢᵀqⱼ ÷ 2.78·10⁻¹⁷the true valuethe norm and its partsdistinct values5of runs10pairs with no fixed sign125of pairs1128κ of this inner product1.1·10¹⁷‖QᵀQ − I‖ moves by1the aggregate is stableand no entry of it is

An inner product with no fixed sign

‖QᵀQ − I‖ is how this site turns "orthogonal" into a number, and across ten partitionings it moves by 2.4%. The entries it is built from are not so lucky: 125 of the 1,128 off-diagonal pairs take both signs, and one of them takes five different values including zero.

machine · Orthogonality 1357911131510⁻¹⁶10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1k, where h = 1 − 10⁻ᵏrelative error of 1 − hsubtractwith a correction stepcorrection kept as a pairfrom the reflectors40 × 6k = 12, subtracting1.1·10⁻⁴k = 12, with a correction1.1·10⁻⁴k = 12, from the reflectors3·10⁻¹⁵reflector operations900dashed: a unit of roundoff over the divisorthree routes sit on it and one does not

The factor a sparse code keeps anyway

Every deletion diagnostic divides by one minus a leverage, and computing it as a subtraction loses a digit for every decade the leverage is from one. The route that does not subtract needs the orthogonal factor, which a sparse factorisation is supposed not to have. Three repairs that avoid it all fail at exactly a unit of roundoff over the divisor — and the fourth, which reaches the orthogonal factor through the Householder vectors a sparse code keeps in order to solve anything at all, returns the same bits as a stored factor in 900 operations.

leastsquares · Leverage heavy row first1 − h at 4³⁰7.4·10⁻¹⁸κ(A) at 4³⁰3.1·10⁸complement's error2.1·10⁻¹⁶110³10⁶10⁹10¹²10¹⁵10¹⁸10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹weight on the rowrelative error of 1 − h1 − ‖R⁻ᵀa‖²1 − ‖row of Q₁‖²‖row of Q₂‖²κ(A)·uthe rows are the same in every positiononly the order the factorisation meets them changes

The weight the factor met first

The route to one minus a leverage through the orthogonal factor was said to lose a digit for every decade of the condition number, whatever else it does. Put a weight on one row and it does not. With the heavy row first, the complement keeps every digit at κ(A) = 2.5·10⁹ while both subtractions return nothing. With the same row last it loses digits as the row's scale grows. And two heavy rows that leave κ(A) at 3.1 still lose six digits when the light rows come first. The law was about the order the factor met the rows, and the condition number had been standing in for it.

leastsquares · Leverage 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 heavy row firstresidual at 4³⁰-2.8·10⁻¹⁰complement's error9.1·10⁻¹⁶110³10⁶10⁹10¹²10¹⁵10¹⁸10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹weight on the rowrelative error of its residualb − Axb − Ax, exactlyQ₂Q₂ᵀbu·wthe same rows, the same data, in every positionthe solution misfits the heavy row by its own rounding

The residual the solution cannot hold

Sorting a weighted fit's rows heaviest first gave every digit of one minus the heavy row's leverage back. It gives nothing back to the heavy row's residual, if that residual is computed the way every textbook computes it — as the datum minus the fitted value. The fitted value is a double, and a double cannot resolve a misfit smaller than its own last digit times the weight: at a weight of 4²⁴ the residual formed from the solution is wrong in its second digit in every order, and forming the subtraction exactly changes nothing. Taken from the same orthogonal factor as the divisor, the residual keeps fifteen digits, and so does Cook's distance at 3.4·10¹⁷.

leastsquares · Leverage 96 × 64, κ = 10⁸blocks of eight, final1.4·10⁻¹⁴their root-sum-square1.2·10⁻¹⁴one at a time, final2.9·10⁻¹⁴081624324048566410⁻¹⁵10⁻¹⁴10⁻¹³columns factorised‖QᵀQ − I‖sum of the blocks' ownone reflector at a timeblocks of eightroot-sum-squareeach block is one factor, each reflector anotherthe departure grows with the square root of the count

Eight blocks and sixty-four reflections

One block of sixteen reflectors, assembled as I − Y T Yᵀ, departed from orthogonality five times as far as a single reflection, and the question was whether a factorisation of many blocks multiplies that factor. It does not. On 96 × 64 matrices eight blocks of eight end at 1.42·10⁻¹⁴ — within a fifth of the root-sum-square of their own departures, and less than half their sum — while the same factorisation taken one reflector at a time ends at 2.91·10⁻¹⁴. A block departs more than a reflector, and there are an eighth as many of them. And a nearly dependent column that swells the triangle's entries to 10²⁶ costs the product nothing.

orthogonality · Householder κ = 10⁸, acrossQᵀb ÷ appended, ρ = 01.3·10⁷Qᵀb ÷ appended, ρ = 10.6710⁻¹⁶10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1relative residual ρforward error010⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1block MGS, Qᵀbblock MGS, appendedHouseholderdotted: κ squared times u times ρthe residual brings every route to the same term

The residual the appended block cannot remove

Appending b as one more block made block modified Gram–Schmidt solve least squares as well as Householder, ten million times better than the same Q through Qᵀb at κ = 10⁸ — on problems with no residual. Give b a component outside the range and every stable route's error rises with it, while Qᵀb's, already at κ²u, does not move. The appended block's advantage then falls as one over the residual: 5,400 at a relative residual of 10⁻⁴, 54 at a per cent, none at one. It never falls behind Householder by more than a factor of four. What the residual decides is whether the extra block is worth its synchronisations, and the answer is yes up to a residual of about a per cent.

orthogonality · Gram–Schmidt median departureone reflector at a time, formed7.5·10⁻¹⁵one reflector at a time, applied3.4·10⁻¹⁵blocks of 64, applied3.7·10⁻¹⁵10⁻¹⁵10⁻¹⁴block sizedeparture from orthogonality12481632642·10⁻¹⁵4·10⁻¹⁵6·10⁻¹⁵8·10⁻¹⁵Q formedapplied through the blocksopen dots: each of the five seedsblocking helps the formed matrix and not the operator

The factor nobody forms

A blocked Householder factorisation's orthogonal factor, multiplied out, departs from orthogonality half as far in blocks of sixteen as one reflector at a time, and that was read as blocking buying a factor of two. Libraries do not multiply it out. Applied to vectors through its stored blocks — which is how every caller uses it — the same factor departs by 3.1 to 4.0·10⁻¹⁵ at every block size from one to sixty-four, and stops growing after about twenty reflectors instead of adding them up. The factor of two was the price of forming the product, and a factor that is never formed never pays it.

orthogonality · Householder a generic tall matrixHouseholder: aimed over median8b appended: aimed over median11column MGS: aimed over median1110⁻⁴10⁻³10⁻²10⁻¹1error ÷ κ²uρthe boundHouseholderb appendedcolumn MGSdots: random directions · bar: aimed · ring: the fixed directionaiming closes a factor of seven

The worst residual belongs to the route

Every stable least-squares route's error rises with the residual at about a fiftieth of the bound κ²uρ, and the question left open was whether a residual aimed at the weak directions closes the gap. It can be aimed exactly: the map from residual to error is, on these problems, one direction and rounding. Aimed, Householder's constant rises from a median of 0.014 to 0.11 — √(m − n) = 6.9 times a typical direction, at 32, 64 and 128 rows to three figures — and stops a factor of nine short of the bound. The direction is each route's own: aimed at one, another route draws a fifth of its worst. And aimed at the appended block, the one-per-cent rule that made it worth its extra block falls to half a per cent, with the block twelve times behind Householder on the same data.

orthogonality · Gram–Schmidt 0.000.050.100.15aimed error ÷ κ²uρill-conditioning spread through the columnsill-conditioning between column blocksHouseholder, sequential0.113tree of two leaves0.090tree of four leaves0.130Householder, sequential0.014tree of two leaves0.042tree of four leaves0.031bars: aimed constantsame order, different directions

A tree leaks along its leaves

A stable least-squares route's error, for a fixed factorisation, is a linear map from the residual with one dominant direction, and each route leaks along its own. A tree of Householder factorisations — the shape a distributed code factorises in — was expected to have a direction of its own too, with b appended at every leaf inheriting Householder's accuracy there. Measured at κ = 10⁸ on two placements of the ill-conditioning, a tree's map has one direction like every route's, and its aimed constant is of Householder's order: within a fifth of it when the ill-conditioning is spread, two to three times it when it sits between column blocks. Its direction is its own — at a cosine of 0.40 to 0.58 from the sequential route's, 0.13 to 0.18 between blocks — and the split of the rows decides it. Appending b at the leaves changes nothing: the same constant to five figures, the same direction to ten.

orthogonality · Gram–Schmidt 0816243240485664reflectors appliedworst ‖Qᵀ(Qx) − x‖02·10⁻¹⁵4·10⁻¹⁵6·10⁻¹⁵8·10⁻¹⁵96 rows192 rows384 rowsmedian of five matrices, 64 columns eachlonger vectors, a higher plateau, reached as soon

The plateau is set by the length

A Householder factor applied through its reflectors, never formed, loses orthogonality over its first twenty or so reflectors and then stops. The saturation was predicted to be a property of the vector length: on 384-long vectors the growth should continue longer before it stops. It does not. On 96-, 192- and 384-long vectors the departure reaches nine-tenths of its plateau within twenty reflectors, and on the longest it gets there soonest. What the length sets is the plateau itself — 3.6, 5.0 and 6.7·10⁻¹⁵ — which is the unit roundoff times the square root of the length, times a constant between 3.1 and 3.4. Twice the reflectors on the same length move it by nine per cent.

orthogonality · Householder

Every trip comes back shorter

An applied Householder factor stops losing orthogonality after about twenty reflectors, at a height that grows like the square root of the vectors' length, because rounding to nearest makes errors of random sign. Rounding every operation toward zero makes errors of one sign, and the plateau was predicted to grow like the length itself — doubling for each doubling. In single precision it does nothing of the kind: from 96 rows to 768, eight times the length, the departure after 64 reflectors grows by a factor of 1.6. An inner product does pile up its errors toward zero, a third of a unit a term, but only when its terms share a sign; Householder's inner products do not, and theirs grow like the square root as before. What rounding toward zero adds is something nearest never does: every vector sent through the factor and back comes home about a unit shorter per reflector, at every length. That shrink is a floor under the departure and has no plateau; rescaling every reflection to its input's length removes it and leaves the departure still climbing.

orthogonality · Householder

Named alongside it

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

OrthogonalityQR factorisationUnit roundoffCondition numberExact ground truthGram–SchmidtLoss of orthogonalityResidualCatastrophic cancellationLeast squaresBackward errorBackward stability

All concepts