Numerical Linear Algebra
where the algebra stops being the arithmetic
Concept

Normal equations — where it appears

The square system AᵀAx = Aᵀb, whose condition number is the square of A's and which is exactly singular below ε = √u. It is taught first and used by nobody, because squaring the condition number loses half the digits before any arithmetic is done.

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

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

leastsquares · Least-squares 10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹ε in Läuchli's matrix (smaller ε, larger κ)relative error in the coefficientsAᵀA exactly singularnormal equationsQRκ from 1.7·10⁸ to 17the cliff is at √u = 2.4·10⁻⁴

The road that squares the problem

The normal equations are the first method every course teaches and the method no library uses. Forming AᵀA squares the condition number, and below ε = √u it does not degrade — it produces a matrix that is exactly singular, from data that was perfectly usable.

leastsquares · Normal equations 01234510⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹log₁₀ κ(A)condition number, and relative errorκ(S)κ(ZᵀHZ)range-space errornull-space erroragainst a BigInt answerκ(S) at κ(A) = 10⁵4·10¹⁰κ(ZᵀHZ), all stops21range-space forward error5.3·10⁻⁶null-space forward error5.8·10⁻¹²both are the same algebraand only one squares

Two ways to remove a constraint

A constrained system can be reduced by eliminating the multipliers or by eliminating the constrained directions. Both give the same answer in exact arithmetic and inherit different condition numbers — one of them squares the constraint's, and the other does not contain it at all.

constraint · Saddle-point systems 1357911131510⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1index kσₖ ÷ σ₁σ₁√utwo factorstheir productthe normal equations, againroutes agree to k =6product floors at9.7·10⁻¹⁰σ₁√u2.3·10⁻⁹κ(P)κ(Q)1.8·10³⁶√ of it1.3·10¹⁸do not form the productthe σ below the line are the bound

The product nobody had to form

The Hankel singular values are the square roots of the eigenvalues of PQ. Form that product and half of them stop existing, at a floor this site can predict from one number — and the fix is the one the least-squares field has had since its first essay, arriving in a place with no least-squares problem in it.

reduction · Balanced truncation 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 00.250.50.75110⁻¹110¹share of the noise placed in the matrixleast-squares error ÷ total least-squares errorequally accuratetotal leastsquares aheadordinary leastsquares aheadthe model, not the methodadvantage, all noise in b0.28advantage, all noise in A2.5seeds at each share40the same total noise at every pointand only where it sits changes

When the matrix is wrong too

Every least-squares problem here has assumed A is exact and b is not, and moved b onto the column space of A. Where both were measured, the smallest correction that makes the system consistent moves the matrix as well — and on the problems where that answer is more accurate, it has the larger residual, by construction rather than by luck.

leastsquares · Total least-squares 10²10³10⁴10⁵10⁶10⁷10⁸10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²condition number‖QᵀQ − I‖one passsweeptwicetreewhat the second pass removesfitted slope, one pass2fitted slope, two passes0.96rounds, two passes6rounds, the sweep24one pass squares the condition numberand two do not

Doing it twice

Cholesky QR squares the condition number — a fitted slope of 1.95 in κ against the Householder sweep's 1.00. Run the identical routine a second time on the Q it returned and the slope is 0.93, the orthogonality is at or below the sweep's at every κ, and the price is one more all-reduce.

cost · Communication -14-12-10-8-6-4-2010⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹10¹⁵log₁₀ μ — the barrier parametercondition number, and relative errorκ₂, condensedκ₂, augmentederror, condensederror, augmentedagainst a BigInt answerκ₂ augmented, μ = 10⁻¹⁴3·10¹⁵its relative error10⁻¹⁵κ₂ condensed2.4·10¹⁶its relative error0.31the same step, written two waysand only one of them is solvable

A condition number sent to infinity

An interior-point method manufactures an ill-conditioned matrix on every iteration, deliberately, because the separating of a diagonal is how it discovers which constraints are active. Written one way the answer keeps fifteen digits at a condition number of 3·10¹⁵. Written the other way — the way almost every code writes it — it has none left.

constraint · Interior-point conditioning 0246810121410⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹10⁴10⁷10¹⁰log₁₀ τ — the weight on the constraintrelative error against the exact answerτ = 1/√uGram–Schmidtnormal equationsHouseholder QRthe ceiling is the method'sHouseholder, τ = 10¹⁴4.8·10⁻¹⁵normal equations14Gram–Schmidt4.9·10¹⁰1/√u6.7·10⁷a constraint is a weight at infinityand the solver decides how far infinity is

A constraint is a weight at infinity

Stack an equality constraint on top of a least-squares problem with a large weight and the answer approaches the constrained one like 1/τ². The limit is takeable to any accuracy — and how far it can be taken is a property of the solver, not of the problem. One of them stops at the square root of the precision, and one of them does not stop.

leastsquares · Constrained least-squares the diagonal of the hat matrix, hᵢ = aᵢᵀ(AᵀA)⁻¹aᵢaverage p/m = 0.20010the leave-one-out residual: eᵢ/(1 − hᵢ), and forty refitsbars: closed form · dots: refitted without that pointone number, two fieldsΣ hᵢ, exactly p10largest leverage0.5closed form against refits8.9·10⁻¹³1 − h of the first row0.5y appears in the residualand nowhere in the leverage

Influence is decided before the data

The diagonal of the hat matrix sums to the number of columns and the response appears nowhere in it, so a fit has exactly p units of influence to hand out among m observations. The same row at h = 0.5 is a ten-fold outlier on one design and a boundary case on another, and which of those it is was settled before a single measurement was taken.

leastsquares · Leverage 110¹10²10³10⁴10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹sweeprelative error, and the largest term's sizerising: the swamp's largest rank-one termfalling, slowly: its errorfalling, once: a fit with an answera plateau with a rising floorswamp error0.0014swamp term10term growth2.2benign error9.7·10⁻¹⁵benign term growth1the error alone cannot tellthe size of the terms can

An iteration that walks out of the set

Every sweep of alternating least squares is the exact minimiser of its own subproblem, so the objective can only fall. What it cannot do is converge, when the target's nearest rank-r point is not in the rank-r set — and a plateau at a small residual looks identical to slow convergence unless the size of the terms is plotted beside it.

tensor · Alternating least-squares 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 -1012345612345678xyall 32without one twinwithout boththe pair at x = 6leverage of each0.44pair's joint share0.89deleted alone, mean0.62deleted together, mean3each twin covers for the otheronly removing both shows the shift

Two observations that hide each other

Two observations at the same place, wrong by the same amount, each look harmless when deleted alone, because a fit without one still has the other. Single deletion sees the shared error cut by (1 − 2h)/(1 − h) — measured at 261 times at the far end — and only the pair's two-by-two block of the hat matrix says what the two of them hold.

leastsquares · Leverage 10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1λ, the ridge on every subproblemerror, and the largest rank-one termthe terms a swamp was growingthe error it costs a fit that has an answerthe error it costs one that has nota bound, at its own priceterms, no ridge6.7terms, heaviest1.2sweeps, no ridge3000sweeps, heaviest69benign error at λ0.089terms still the right terms1the plateau is boundedand the answer is biased by exactly λ

The repair that costs exactly itself

A ridge on every subproblem is the standard cure for a swamp and it works — the terms stop at 1.61 instead of climbing past 6.7, and a run that never finished finishes in 798 sweeps. On a tensor that does have an answer the error it costs is the ridge itself, to within a factor of two, at every setting from 10⁻¹⁰ to 10⁻¹. And the terms it returns are still the right terms.

tensor · Alternating least-squares what each reachesbest pair0.0034‖x‖ alone0.0048‖L₁x‖ alone0.0035λ₁, on ‖x‖λ₂, on ‖L₁x‖ — log₁₀0-6-5-4-3-2-100-6-5-4-3-2-10outlined row and column: one penalty switched offred cell: the best pair, 0.00343darker is worse, on a logarithmic scalethe optimum sits on or beside an edge

A second penalty is not a second parameter

Penalise ‖x‖ and ‖L₁x‖ at once and there are two λ to choose. Over ten draws on five signals the best pair beats the better single penalty by between 0.00% and 3.1%, and one of its two parameters is exactly zero on 30 to 70 per cent of draws. Choosing the wrong one of the two costs up to 54%. The surface is a choice between two curves with a knob nobody needs.

regularisation · Parameter choice 10¹10³10⁵10⁷10⁹10¹¹10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1κ(A), the conditioning of the fitrelative error in the answerthe saddle-point route, in floating pointdashes: forming AᵀA, then solving exactlythe method of weighting, τ = 10⁸the null-space routethe reference was a methodκ(A)10¹¹κ of what is left4.9·10⁶null space1.9·10⁻⁹saddle point4.6·10⁻⁴forming AᵀA alone3.4·10⁻⁴weighting3·10⁻⁹the damage is in the formingand not in the solving

The reference was a method

The optimality conditions of a constrained fit contain AᵀA, so solving them is the road that squares the problem wearing a block structure. At κ(A) = 10¹¹ the route that never forms a cross-product returns 1.89·10⁻⁹ and the route that does returns 4.64·10⁻⁴ — and forming AᵀA and then solving it in exact rationals returns 3.45·10⁻⁴, so nearly all of the loss happens before any elimination begins.

leastsquares · Constrained least-squares 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 10²10⁴10⁶10⁸10¹⁰0285684112140168196224condition numberstepsnormal equationsbidiagonalisationone sequence, two costssteps at κ = 10², both16at κ = 10⁶, ratio1.1at κ = 10¹⁰, ratio1.9the same iterates in the algebraand twice the work at κ = 10¹⁰

One sequence and two recurrences

CGLS and LSQR compute the same iterates — the minimiser over a space is unique, so there is nothing to choose between them in the algebra. At κ = 10⁶ they cost 42 steps and 47. At κ = 10¹⁰ they cost 110 and 209, across four seeds, and the quantity that separates them is the orthogonality of a basis neither of them keeps.

iterative · Krylov ‖QᵀQ − I‖ of the implied Qκ = 10⁸, Cholesky0.37κ = 10⁸, sweep8.5·10⁻⁹κ = 10¹⁰, Cholesky1.3κ = 10¹⁰, sweep3·10⁻⁷κ = 3·10¹⁰, Choleskyrefusedκ = 3·10¹⁰, sweep6.7·10⁻⁶κ = 10¹², Cholesky1.5κ = 10¹², sweep1.4·10⁻⁴the safe run is the one that failsrefusals in the range1wholly non-orthogonal returns2the pivot it refused on-6.5·10⁻¹⁷one of these outcomes is safeand it is the refusal

The licence is not the boundary

Cholesky QR is licensed by κ²u ≪ 1, which reaches equality at κ = 9.5·10⁷ in double precision. At 10⁸ the factor it returns is already 0.37 away from orthogonal, and it goes on returning factors as far as 10¹³ — refusing at scattered condition numbers in between, at different ones for eight columns and for six.

machine · Algorithm selection

Named alongside it

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

Condition numberLeast-squaresResidualUnit roundoffQR factorisationCondition squaringForward errorGram matrixLeverageOrthogonal projectionSaddle-point systemsExact ground truth

All concepts