Numerical Linear Algebra
where the algebra stops being the arithmetic
Theme

The thread: Whose fault is it — page 3

Essays 49 to 72 of the 252 on this theme, in the same order.
110¹10²10³10⁴10⁵10⁶10⁷10⁻¹⁷10⁻¹⁶10⁻¹⁵10⁻¹⁴10⁻¹³10⁻¹²10⁻¹¹10⁻¹⁰10⁻⁹10⁻⁸1/(1 − h), the leverage of the removed row‖R̄ᵀR̄ − (G − aaᵀ)‖ / ‖G − aaᵀ‖downdatedrefactorisedat h = 1 − 10⁻⁷κ of the downdated matrix4.3κ of the matrix downdated9.3·10⁶rotation's amplification344downdate residual3.5·10⁻¹⁰a hyperbolic rotation is not orthogonaland that is exactly what it is for Least squares, and the road not to take

The observation that cannot be removed

Removing a rank-one term from a Cholesky factor needs a rotation that is not orthogonal, and the number under its square root is 1 − h, where h is the leverage of the row being removed. The algorithm's breakdown condition and the statistician's warning are the same quantity, arrived at from opposite ends, and neither field states it in the other's language.

 567891003691215log₂ nlargest rank in the partitionweak: touching pairs compressedstrong: touching pairs refusedthe test bounds a rank and costs storagestrong, blocks250weak, blocks94strong, numbers6.8·10⁴weak, numbers6.1·10⁴strong ⁄ weak1.1the better partitionis the more expensive one Neither sparse nor dense

The test that costs what it saves

The partition that refuses to compress a touching pair keeps every rank at five while the other lets them climb from nine to thirteen. It also stores more numbers at every size measured — 67,968 against 61,440 at n = 512 — and which of those two facts matters is a question about how large the problem is going to get.

 10¹10²10³03691215what one rebuild is worth, in preconditioned iterationscheapest number of members between rebuildsa dense Cholesky here is worth 6.7 iterationsringed: where the growth rule beats every fixed periodcheapest period, by setup costsetup worth 51setup worth 102setup worth 203setup worth 506setup worth 1008setup worth 20012how long to keep itis a question about what it cost to build When the problem arrives again

What a rebuild is worth

One sequence, one drift, one preconditioner — and six different right answers, because the cheapest rebuild period depends on what a rebuild cost to build and on nothing else. The optimum walks from every member to every twelfth as the setup gets dearer, and the free rule that reads the iteration count beats it in the middle of that range and loses at both ends.

 01234567891011121310⁻¹110¹10²10³iterationinterpolation point σa fixed point that is a conditionorder4iterations13σ against −λ(Aᵣ)2.8·10⁻¹²H₂ error4.3·10⁻⁵it stopped movingand the condition holds there Reduction, and what a model is for

Interpolating at the model’s own poles

One choice of interpolation points is not arbitrary — the mirrored poles of the model about to be built. It is a fixed point rather than a guess, and when it is reached it beats a method costing O(n³) — by 0.4 per cent, which is the honest size of the whole contest.

 14710131619110¹indexeigenvalue of Lan eigenvalue with no vectorλ₂0.35λ₃0.35gap1.7·10⁻¹⁵partitions found3best conductance0.045worst0.045two eigenvalues, one valueand the answer is not a function of the graph The matrix that is a graph

A partition decided in the last digit

On a graph with a symmetry there is no Fiedler vector — there is a plane, and every vector in it is an exact eigenvector. Twenty-four runs with the edge weights nudged by 10⁻¹² return ten different partitions of a cycle and, on a hypercube, two different qualities of answer.

 10⁻¹⁴10⁻¹³10⁻¹²10⁻¹¹10⁻¹⁰10⁻⁹10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1relative size of the entrywise perturbationrelative forward errorκ∞ · εcond(A,x) · εmeasuredboth bounds holdκ∞(A)1.9·10⁸cond(A, x)4.8ratio of the bounds4·10⁷both curves above the data are boundsand only one of them is a measurement Two errors, and whose fault they are

A condition number scaling cannot move

Skeel's componentwise condition number is invariant under any row scaling — exactly, before any norm is taken, because two diagonal factors cancel entry by entry. It is never larger than the normwise one and can be arbitrarily smaller, and the ratio between them is a diagnostic for which kind of ill-conditioning a matrix has.

 081624324048566410⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹index kmagnitudethe floor: k = 43best truncation: k = 31σₖ|uₖᵀb| exact|uₖᵀb| with noise|uₖᵀe|, smoothedcorrelation ρ = 0.9noise energy, first 28 directions0.98tilt of the noise floor, decades1.1the crossing43the best truncation31the floor tilts, and is still a floorρ = 0.9 Regularisation, and the answer that is chosen

Noise that spares the answer and fools the rules

Make each noise sample remember the last one, keep its size fixed, and the best answer available gets slightly better — 0.1056 to 0.1010 — because slow noise hides in the directions where dividing by σ costs nothing. The Picard crossing still lands two dozen indices past the best truncation. What breaks is the rules. Generalised cross-validation more than doubles the best error on 14 draws of 48 instead of 3, the discrepancy principle's typical cost triples, and the two miss in opposite directions. Whitening by the covariance takes GCV back to 3.

 D from PAPᵀ = LDLᵀ — the shaded pairs are 2×2 pivots10⁻⁶0.749······0.749·········1.2·10⁻⁶1.4······1.4·········2.1·10⁻⁶0.549······0.549·········3.6·10⁻⁶0.614······0.614·three rules, one matrix‖PAPᵀ − LDLᵀ‖, blocks5.8·10⁻¹⁷‖PAPᵀ − LDLᵀ‖, diagonal3.1·10⁻¹¹growth, blocks1.3growth, diagonal5·10⁵the zero block is what the problem saysand one rule does not need it to be nonzero Elimination, and the swap

When symmetry is not enough

The matrix [[0, 1], [1, 0]] is symmetric, nonsingular and perfectly conditioned, and there is no diagonal entry to pivot on. Every factorisation restricted to symmetric interchanges and one-by-one pivots fails on it, at any depth of searching, because every entry it could search is zero. The repair is to take two variables at once.

 κ(A) = 10⁶ throughout · κ(H) = 100 · the answer is the same answer for every basisorthonormal — κ(Z)1κ(ZᵀHZ)25.6relative error1.07·10⁻¹⁵first m basic — κ(Z)1.99·10⁸κ(ZᵀHZ)3.8·10¹⁶relative error0.0518pivoted basic — κ(Z)2.06κ(ZᵀHZ)31.9relative error6.71·10⁻¹⁶what the choice costsdensity, orthonormal1density, fundamental0.5κ(ZᵀHZ) ÷ κ(Z)², naive0.96error, pivoted choice6.7·10⁻¹⁶every one of them is a basisand one of them loses fourteen digits Orthogonality, measured

The basis nobody chose on purpose

A method that eliminates a constraint has to pick a basis for its null space, and every basis is correct. Their condition numbers are eight orders apart, the reduced problem inherits the square, and the choice is usually made by a one-line rule nobody thought of as a numerical decision.

 10⁻²10⁻¹110¹10²10³10⁴10⁻⁴10⁻³10⁻²10⁻¹1backward error, in units of ushare of systems above ituCramereliminationCramer, control24-bit arithmeticworst Cramer, in u395worst elimination, in u1.2control, worst Cramer4.9κ of the worst system3.4·10⁶one derivation, two computationsand only one of them is stable Elimination, and the swap

A rule that is correct and unusable

Cramer's rule gives every component of the solution in closed form, in terms of determinants, and it is a theorem. On two-by-two systems whose rows are nearly parallel it returns an answer with a backward error of 458 units of roundoff where elimination returns 1.3 — on a matrix whose condition number is 32,000 and which elimination solved perfectly.

 02468101210⁻¹⁸10⁻¹⁶10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸member of the sequencebackward error of the solvefresh order · equilibrated kept orderringed: another pivot replaced hereone order, two sets of unitsfresh order, last member9.5·10⁻¹⁷kept order, last member4.8·10⁻⁹equilibrated, last member6.5·10⁻¹⁷pivots replaced, kept2pivots replaced, equilibrated0the perturbation floor √u1.1·10⁻⁸a reused pivot order is safe or notdepending on what the rows are measured in Sparsity, and what elimination costs

The order that was right last time

A pivot order computed once and reused across a sequence saves the symbolic phase, and the price is that a pivot which was large may now be small. Replacing it with √u·‖A‖ costs eight orders of backward error and iterative refinement recovers a factor of 8.8 of them. Divide each row by its largest entry first and the same reuse costs nothing at all.

 02468101210⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²binary64 operations per elementrelative errorruntime orderfixed partitionpre-rounded, one passpre-rounded, two passescompensated (Kahan)exact accumulationwhat determinism costsreproducible policies3cheapest, ops/element3its error0.0024two passes10⁻¹⁷ordinary, 8 pieces3.6·10⁻⁷nothing is cheap and exactthe middle of the figure is the answer The answer that depends on the machine

Accuracy and agreement are different properties

The most accurate policy on this site's summation figure returns 119 different answers, and the one that returns a single answer is four orders less accurate. Neither property implies the other, and the vocabulary has one word for both.

 26101418222630343810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹stepsizeleast error: 20‖QᵀQ − I‖filter disagreementan identity with an expiry datedisagreement at step 1210⁻¹³disagreement at step 400.014least error at step20exact while the basis is orthogonaland false where the method is best Methods that were designed apart

An expiry date the noise does not move

The polynomial description of conjugate gradients leaves the level of rounding at step 17 or 18 on this operator, at every noise level from 10% to 0.1%. The step worth stopping at moves from 3 to 44 across the same range. They coincide at about 1% noise, which is where the coincidence was first read, and it is a fact about the noise rather than about the method.

 110¹10²10³10⁴10⁵10⁶error ÷ the oracle's, on the same drawten times the oracleover tendiscrepancy principletold ‖e‖0 of 1000unbiased risktold σ²15 of 1000cross-validationtold nothing57 of 1000quasi-optimalitytold nothing0 of 1000L-curvetold nothing0 of 1000blue: median · bar to the 90th percentile · line to the 99th · red: worstthe counts are the tail Regularisation, and the answer that is chosen

One draw in twenty

Sixteen draws gave generalised cross-validation a worst case of 12%. A thousand draws at each of five noise levels give it a second answer on four to six in every hundred, ten to seven million times worse than the oracle, while its median stays among the best of five rules. The quasi-optimality criterion, told nothing either, never costs more than 1.41 in five thousand draws. The share settles by a thousand draws, and letting the search look further down more than triples it.

 the smallest perturbation of any kind — 1.89·10⁻¹⁷the smallest Toeplitz one — 2.65·10⁻¹³both exact for the same x̂smallest of any kind1.9·10⁻¹⁷smallest Toeplitz one2.7·10⁻¹³the price of the constraint1.4·10⁴diagonal defect, unconstrained0.97an exact answer to a nearby problemof a kind nobody posed Structure, and the solver that cannot see it

A nearby problem of the wrong kind

A good algorithm returns the exact answer to a nearby problem. A hundred and eighteen essays have measured the distance and not one has asked what the nearby problem looks like. On a Toeplitz system it is a rank-one matrix that is constant along none of its diagonals — and the smallest one that is Toeplitz is two and a half million times larger.

 the estimator maximises this quantity over the columns it visitscolumn 1 ‹visited›12column 2 ‹the answer›114column 311.4column 411.4column 511.4column 611.4column 711.4column 811.4column 911.4column 1011.4column 1111.4column 1211.4estimate 12.0a walk that stopped earlythe estimate returned12the true 1-norm114columns visited1products with the matrix5the walk's own stopping test firedand every column it could see was smaller Two errors, and whose fault they are

An estimate that can be fooled

Nobody computes a condition number, because forming an inverse costs more than the solve did. Every library estimates it instead, from four or five products with a factorisation already in hand. The estimate is exactly right on four random matrices out of five — and there is a matrix, three distinct entries wide, on which it returns a twentieth of the truth.

 110¹10²10³110¹10²10³off-diagonal entry ccondition number of the eigenvalue√(1 + c²)decoupled: 1measuredthree routes, one number‖A − ZTZᵀ‖/‖A‖1.7·10⁻¹⁵closed form100computed 1/|yᵀx|100worst measured movement46four eigenvalues, two conditioning numbersthe symmetric case has one, and it is 1 Eigenvalues, singular values, rank

A condition number for one eigenvalue

In the symmetric case every eigenvalue has condition number exactly one. In this four-by-four matrix two of them have condition number 100.005 and the other two have exactly 1, and the number belongs to the eigenvalue rather than to the matrix.

 10¹10¹.³10¹.⁶10¹.⁹⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰¹10².²00.250.50.751grid points nresidual reduction per stepJacobiGauss–SeidelV-cycleV-cycle spread, 8× in size0.003Jacobi at n = 1270.99work exponent, fitted0.079the dashed curve is cos(πh), Jacobi's closed formthe flat line is the whole method Iterating, instead of factorising

A rate that does not notice the size

The V-cycle reduces the residual by a factor of ten a cycle at fifteen points and at a hundred and twenty-seven. Jacobi on the same four problems goes from 0.981 to 0.9978, climbing towards one. One of those is a constant and the other is an exponent, and that is the whole distinction the field turns on.

 23456710⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²poles in the approximantresiduals and errorforward error‖T(λ)x‖‖T̃(λ)x‖one extra evaluationagainst the approximant1.5·10⁻¹²against the problem asked1.9·10⁻⁶forward error4.8·10⁻⁵‖g − r‖ there8.5·10⁻⁵the free residual is flatand the answer is not The eigenvalue problem that is not linear

The problem the solver was actually given

A linearisation is exact — it has the polynomial's eigenvalues, with their multiplicities, and the whole loss is arithmetic. A nonlinear eigenvalue problem does not offer that. Every algorithm replaces the function first, and the term that replacement contributes is committed before any number is rounded and appears in no residual.

 2022242628303210⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²k, where the entries are near 2ᵏrelative error of the determinantas writtenfused: exactproducts need 54 bitsone rounding, the whole answertrue determinant1naive, k = 300fused, k = 301first wrong at k27sizes returning 06both forms conformand the source does not say which The answer that depends on the machine

One multiply the compiler removed

A determinant whose value is exactly 1, computed as exactly 0 by the expression that is written down, and exactly 1 by the same expression with the multiply and the add fused. Both forms conform to IEEE-754, both are legal compilations of the same source, and nothing in the program says which one you have.

 0173451688510211910⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹1steps of the walkdistance from stationarya rate read as a timeλ₂ of the walk0.9steps measured122steps predicted131ratio0.93the dashed line is the eigenvaluethe curve is the walk The matrix that is a graph

The rate is the second eigenvalue

A walk forgets where it started at a rate the graph's second eigenvalue names exactly. Across three orders of magnitude in the step count the prediction is five per cent high — and the published rate for PageRank is right for a reason nobody states, which is that a link graph is in pieces.

 the diagonal of the hat matrix, hᵢ = aᵢᵀ(AᵀA)⁻¹aᵢ · dashed: its average p/m = 0.200p/m10the 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 Least squares, and the road not to take

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.

 10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹110⁻⁴10⁻³10⁻²10⁻¹110¹10²10³10⁴λnoise part ÷ signal part, in ‖x‖where every corner sitsthe corner, ×2.29the oraclethe corner reads ‖x‖noise share at the corner0.23noise share at the oracle0.026corner ÷ oracle, this draw2.3corner ÷ oracle, 60-draw median1.5noise share: ‖L·(noise part)‖ ÷ ‖L·(signal part)‖the corner reads that share Regularisation, and the answer that is chosen

The corner reads the norm it is drawn in

The L-curve was the costliest rule this field scored, and the cost was not the rule's. On the same sixty draws, with the same best achievable error, the corner of ‖x‖ against the residual costs 1.53 times the oracle and the corner of ‖L₁x‖ costs 1.003. Across five signals and three penalties the corner lands wherever amplified noise is between a tenth and a fifth of the norm being plotted, and it finds the oracle only when the oracle happens to sit there — twenty-nine times too costly on a smooth signal under ‖x‖, within half a per cent on four spikes.

 unpreconditionedbest step20best error0.14Tikhonov's best0.14α = 0.001best step1best error0.14eigenvalues sent near one22051015202530354010⁻¹110¹steprelative errorTikhonov's best: 0.1405plain CGLSpreconditioneda better preconditionerarrives at the noise sooner Methods that were designed apart

A preconditioner that arrives past the answer

On a system that is solved to convergence a preconditioner changes how fast the answer arrives and not what it is. On a problem regularised by stopping it changes where every step lands. Conjugate gradients preconditioned by AᵀA + αI reaches its best answer in one step at α = 10⁻³, and at α = 10⁻⁶ its best answer is its first step, with an error of 1.35 against the unpreconditioned run's 0.1426 — while the count of eigenvalues it has clustered at one rises from 22 to 32.

All themes