Numerical Linear Algebra
where the algebra stops being the arithmetic
Series

Householder — the series

7 essays on one idea, from the one that introduces it to the one that assumes the rest.
  1. 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.

    part 2 · orthogonality
  2. 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.

    part 3 · orthogonality
  3. 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.

    part 4 · orthogonality
  4. 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.

    part 5 · orthogonality
  5. 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.

    part 6 · orthogonality
  6. 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.

    part 7 · orthogonality
  7. 10¹10²10³rowsdeparture after 64 reflectors, over u100200400700toward zeroto nearestgrowing like mgrowing like √mdotted: like m · dashed: like √m, both from the first pointnot like the length

    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.

    part 8 · orthogonality

All series