Numerical Linear Algebra
where the algebra stops being the arithmetic
Theme

The thread: Whose fault is it — page 6

Essays 121 to 144 of the 303 on this theme, in the same order.
unpreconditionedbest step20best error0.14edge τ0.01cutoffs at τ = 0.5λdiscrepancy0.032GCV0.016L-curve0.00410⁻⁴10⁻³10⁻²10⁻¹10⁻¹110¹truncation τbest relative errorunpreconditioned floor: 0.1426discrepancy: 0.144GCV: 0.142L-curve: 0.234the cutoff read as 0.5 × each rule's own λthe plateau is wide and the cliff is steep Methods that were designed apart

The rule that is wrong in the right direction

The preconditioner's cutoff is not a new parameter. It is the regularisation parameter this field already knows how to choose, halved — and the rule criticised for choosing λ a factor of two or three too large is the one whose cutoff keeps the floor on every draw, where the rule that chooses λ to within 3% has a worst draw two hundred times off it.

 where each λ sitsoracle's share0.0037corner's share0.1corner's cost25and what a target costsbest fixed target0.0035its cost1estimate ÷ truth at the oracle0.9210⁻⁸10⁻⁶10⁻⁴10⁻²110⁻³10⁻²10⁻¹110¹10²10³10⁴λnoise part ÷ signal parttarget 0.0035oracle λ: share 0.0037corner: share 0.102solid: the share on this draw · dashed: the share a rule can estimatethe rule is a level set of the dashed curve Regularisation, and the answer that is chosen

A rule that has to be told how good its answer will be

The L-curve's corner reads a noise share of 0.10 to 0.21 across fifteen pairings of signal and penalty, and the share the best λ sits at runs from 0.0037 to 0.43 — a factor of a hundred and fifteen. A rule aimed at the right share is within a few per cent of the oracle on every one of them. The right share is about a third to four-fifths of the relative error that λ will achieve, which is the number the answer was wanted for.

 the damagewell scaled, nothing done0.9110^±3, nothing done0.095never certified, of 62at 10^±3rows to unit norm0.64rows by right-hand side0.6start at the rows0.43012300.250.50.751rows rescaled by 10ᵏshare that certifiesnothing donestart at the rowsrows to unit normrows by their right-hand sidethe grey line is a well-scaled programme with nothing doneboth equilibrations are flat, and below it The matrix a constraint makes

Two repairs for one symptom

Rescale a quadratic programme's constraint rows over six decades and the crossover that certified its answer at iterate 1.5 first certifies at 69.8, with two of six programmes never certifying at all. Normalising the rows removes the spread completely — the same numbers at 10¹, 10² and 10³ either way. Starting the method at the magnitudes the rows imply repairs the iteration count completely and the identification only halfway. They are two repairs and they fix different halves.

 8×8 standard normalmatrices60median growth, fully pivoted1.3costliest single decision1.6step 1 — median cost1.624×step 1 — moves a row90%step 2 — median cost1.545×step 2 — moves a row87%step 3 — median cost1.419×step 3 — moves a row88%step 4 — median cost1.333×step 4 — moves a row80%step 5 — median cost1.058×step 5 — moves a row75%step 6 — median cost1.021×step 6 — moves a row80%step 7 — median cost1.000×step 7 — moves a row43%red: growth when that one search is skippedblue: how often that search moves a row Elimination, and the swap

Which of the choices is doing the work

Elimination makes n − 1 decisions and they are not worth the same. On 8×8 standard normal matrices, removing the first pivot search and leaving the other six multiplies the median growth factor by 1.624; removing the last multiplies it by 1.000. The cost falls monotonically along the run, and the worst single matrix in the sweep grows by 2,366 when one early decision goes — so the median is the wrong statistic and the tail is where pivoting earns its reputation.

 total inner iterationscold at η = 1e-107016scaled warm start6837unscaled728910⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²10³10³.⁵10⁴forcing term ηinner iterationsthe previous stepstarted at zerothe previous step, scaledthree curves within eight per cent of each otherand the unscaled guess is the worst of the three When the problem arrives again

A guess worth two per cent

The previous Newton step looks like a free guess at the next one, and it is worth nothing. Started from it unscaled, the inner solve costs 4 to 61 per cent more than starting from zero, because the guess is 15 to 209 times too large. Scaled by the ratio of the two residual norms it is the right size and halves the starting residual — which buys a constant handful of inner iterations, not a share, because conjugate gradients costs the logarithm of its tolerance.

 48 products, 24 drawsbest share, decay 0.70.45best share, decay 0.950.2worst cost of a third2.800.10.20.30.410⁻³10⁻²10⁻¹share spent on the sketchmedian relative errorthe published thirddecay 0.7decay 0.85decay 0.95large dots: the best split on each curveand the dashed line is the one a library picks Randomised, and the guarantee that changes kind

The split nobody is in a position to choose

Hutch++ spends two thirds of its budget on a sketch and a third on probes, and the third is published as a constant. Swept across six rates of spectral decay at a fixed budget of 48 products, the best share is 0.45 on the fastest and 0.00 on the slowest — sketch nothing at all — and the published third costs between 1.09 and 4.11 times the best error. The decay that decides it is readable from the sketch's own singular values, for products the estimator was going to spend anyway.

 an eigenvalue repeated 2 times, in a 200-dimensional problemblock of 1never returns them allblock of 232 products, 16 stepsblock of 351 products, 17 stepsblock of 460 products, 15 stepsblock of 570 products, 14 stepsblock of 684 products, 14 stepsthe narrowest that workscheapest width2its products with A32the widest drawn84narrower than the multiplicity is not slowit is absent Eigenvalues, singular values, rank

How wide the block should be

A block narrower than the multiplicity does not converge slowly — it never returns the missing copy at all. Above the multiplicity every extra column buys iterations at about ten products with A each. And the mechanism that is supposed to make the choice unimportant never fires from a random start.

 ‖H − Hₙ₋₁‖∞ ÷ σₙpositive residues2.0000mixed residues2.0000clustered poles2.0000four poles2.0000wide spectrum2.0000‖ΔH‖ 3.89·10⁻⁴ σ 1.94·10⁻⁴‖ΔH‖ 3.92·10⁻⁴ σ 1.96·10⁻⁴‖ΔH‖ 0.00502 σ 0.00251‖ΔH‖ 8.41·10⁻⁴ σ 4.21·10⁻⁴‖ΔH‖ 4.6·10⁻⁴ σ 2.3·10⁻⁴2four systems, four ratios, one number: 2.0000one state, one equalityand the sum is that, repeated Reduction, and what a model is for

The state that is removed is not a mode

Balanced truncation removes one state and pays exactly twice one Hankel singular value. The natural reading is that the state removed is the model's least important mode and that the σ is that mode's own size — and on five systems that reading over-estimates by between 1.21 and 7.46, never once under.

 57911131517192110⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²degreerelative error of the worst rootmeasuredκ × ua polynomial is its coefficientsworst root, measured0.0076predicted, κ × u0.0023root condition number10¹³largest coefficient1.4·10¹⁹the roots are integersand the coefficients are not the roots Two errors, and whose fault they are

The roots are not the coefficients

A polynomial whose roots are the integers one to twenty, expanded exactly, handed to the routine every library uses. The computed roots are wrong in the third digit, the computation is backward stable for the matrix it factorised, and above degree eighteen the coefficients are not double-precision numbers at all.

 10³10⁶10⁹10¹²10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹κ(A), identical for both problemsrelative error in the answerthe weak directions in the constraint's null spacethe weak directions in the constraint's own rowsone κ(A), two problemsκ(A), both10¹²κ left, constrained1κ left, free10¹²error, constrained4.7·10⁻¹⁶error, free3·10⁻⁴between them6.4·10¹¹a constraint is informationand κ(A) does not know it arrived Least squares, and the road not to take

The condition number that does not know

Two constrained fits with the same size, the same number of constraints and the same κ(A) to twelve figures. One returns 4.7·10⁻¹⁶ and the other 3.0·10⁻⁴. What separates them is the conditioning of A restricted to the constraint's null space — 1.00 against 10¹² — which every solver computes on the way and none reports.

 total operationswork, no dissection10⁵work, full dissection1.5·10⁵work, one level1.3·10⁵on unbounded processorscritical path, no dissection4.1·10⁴least, at depth 43.4·10⁴critical path, one level5.4·10⁴10⁵switching depthoperations01234568total operationscritical pathdepth 0 is minimum degree, depth 8 is nested dissectiondepth 1 is worse than both Sparsity, and what elimination costs

The depth that is worse than both ends

Nested dissection to a chosen depth and minimum degree below it is the ordering codes ship, and sweeping the depth was supposed to find a setting that keeps most of dissection's parallelism for most of minimum degree's work. It does not exist: total work rises with the depth at every grid size, and one depth — the first — is worse than both extremes on work and on the critical path at all four sizes measured. One bisection buys nothing because there is no recursion under it to amortise the separator.

 what one line doesinner work saved, worst case0.09best case0.27forward error cost, worst263610⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴outer toleranceforward errorerror, with the floorerror, without itinner work, rescaledthe grey diagonal is the residual the test asked forboth runs meet it and one is far more accurate When the problem arrives again

One line that buys a quarter of the run

The adaptive forcing rule has a floor on it that no published statement of the rule carries: do not solve a step to an accuracy the outer loop will not use. Removing it costs 9 to 27 per cent of the whole inner run. Keeping it costs between 23 and 2,600 times the forward error — accuracy the residual test never asked for and both runs satisfy the test either way. The line is a trade between a residual and an error, and which of the two the caller meant decides whether it is a saving.

 10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹10³10⁵σ — how far the periodic wrap is from singularrelative forward error10⁻¹10⁻²10⁻³10⁻⁴10⁻⁵10⁻⁶10⁻⁷10⁻⁸10⁻⁹10⁻¹⁰10⁻¹¹10⁻¹²periodic correctionκ(C)·uelimination on Tn = 64, median of five answersperiodic correction, σ = 10⁻⁸1.2·10⁻⁴elimination on T, σ = 10⁻⁸1.1·10⁻¹⁴κ(C) at σ = 10⁻⁸4·10⁸κ(T) at σ = 10⁻⁸1712the same matrix T in every columnonly the wrap it is solved through changes Structure, and the solver that cannot see it

The circulant the problem did not contain

A matrix that differs from a circulant in two corner entries can be solved through the circulant, by a transform and a two-by-two correction, and the cost claim is exact. The accuracy claim is not. On tridiag(−1, 2 + σ, −1), whose condition number stops at 1,712, the correction is wrong by 1.2·10⁻⁴ at σ = 10⁻⁸ while elimination is right to 1.1·10⁻¹⁴ — because the periodic neighbour is singular at σ = 0 and the two-by-two system inherits that. Solved by Cramer's rule, as here, the two amplifications multiply; a later measurement found that a pivoted solve of the same two-by-two system removes the second.

 081624324000.250.50.7511.25angle of the strong direction (degrees)ratio, and convergence factorenergy ratio of oneenergy ratioclassicalminimiserwith the constraintwhere the formula is optimalenergy ratio at 0°1energy ratio at 45°1.2classical rate at 45°0.35the minimiser's rate there0.28a ratio of exactly one while the assumption holdsand a diagnostic when it stops Iterating, instead of factorising

The formula that was already optimal

Ask for the interpolation that minimises the energy of its own columns and the answer is the classical AMG formula — to zero at every row of the one-dimensional Laplacian, and to four digits in two dimensions. On the operator rotated to 45° the two part company, and the gap between them is a diagnostic that needs no reference solution.

 the problem you posedA = H10b = A·(1, 2, …, 10)the problem it answered exactlyA + δA, b + δb‖δ‖ / ‖A‖ = 2.3·10⁻¹⁷the answer you wantedx = (1, 2, …, 10), exactlythe answer you gotx̂, wrong by 2.7·10⁻⁴ relativebackward error 2.3·10⁻¹⁷forward error 2.7·10⁻⁴κ = 1.6·10¹³κ · η = 3.6·10⁻⁴, and the measured forward error is 2.7·10⁻⁴.The algorithm is not at fault. The problem is.H10, LU with partial pivotingresidual and error differ Two errors, and whose fault they are

Three errors and one number

This site's identity has two factors and a division of blame between them. Two fields have now added a third party and a fourth, and only one of the four is a property of anything — the others are decisions, made before the arithmetic, reported by nothing.

 -15-13-11-910⁻¹³10⁻¹¹10⁻⁹10⁻⁷log₁₀ of the defect, relative to the entry it sits inhow far the answer movedthe machine, on its owna tolerance with two sidesmachine band3.2·10⁻¹²smallest caught10⁻¹²its nearest run4.6·10⁻¹²defects hidden3window, factor1.4below the band nothing is visibleand above it everything is The answer that depends on the machine

What a regression test can ask for

The machine's own variation on one solve is 3.2·10⁻¹², and the smallest defect whose answers clear it is one part in 10¹². The tolerance exists, it is bracketed on both sides by a factor of 1.42, and it is neither zero nor the 10⁻⁸ that usually gets typed.

 10⁻²10⁻¹110¹10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹|λ|backward error of the computed pairone residual, four denominatorspairs measured16worst, all three1.4·10⁻¹⁵worst, K alone1.8·10⁻¹⁴worst, M alone9·10⁻¹²K-only factor, low1.1high46solid: every coefficient may movedashed: only one may The eigenvalue problem that is not linear

A perturbation that moves every coefficient

The backward error of a polynomial eigenpair is measured against perturbations of all three coefficients at once. Restrict it to the one coefficient anybody is willing to move and the same computed answers are stable at one eigenvalue and unstable at another, by a factor that runs from 1.06 to 6,370 across a single spectrum.

 10⁻¹110¹10²10³1234eigenvalue of the preconditioned Laplacian, and the stretchbreadth-firstκ = 34.62randomκ = 55.17random, secondκ = 37.11the bar is the spectrumthe tick is a combinatorial bound The matrix that is a graph

A preconditioner that is a tree

Every eigenvalue of a tree-preconditioned Laplacian is at least one and at most the total stretch — a combinatorial integer with no arithmetic in it. Measured, the bound is two to four times loose, and on a grid the preconditioner makes the conditioning worse by a factor of 1.85 at every size.

 n = 512, ε = 10⁻⁸least storage, at leaf16numbers per unknown there133at the largest leaf211050100150200leaf sizestored per unknown48163264stored per unknownshare in the dense blockslarge dot: the least storage availableand the dashed curve is why there is one Neither sparse nor dense

Two knobs on one number

The tolerance raises every block's rank and leaves the partition alone. The leaf size does the opposite — it changes how many blocks there are and does not move a single rank. Both act on the same storage, and only one of them has a best setting: at eight digits the least storage is at a leaf of 16 and the largest leaf tried costs 59 per cent more. The best leaf rises with the accuracy, from 4 at two digits to 16 at twelve, and the influence runs only that way.

 -7-6-5-4-3-210⁻¹⁶10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²110²10⁴log₁₀ of the twist's distance from landing, in turnsrelative forward errorsimple zero, landed on at 0simple zero, landed on at 0.3double zero, d = 2n = 64, twist 10⁻⁵ of a turn from landingpair: error, κ(A) = 4.1·10⁶5.4·10⁻¹²single: error, κ(A) = 3.8·10⁶3·10⁻⁸double: error, κ(A) = 4.2·10¹²7515the same condition number of the wrapa different number of modes near zero Structure, and the solver that cannot see it

Two near-zeros cost less than one

Solve a well-conditioned tridiagonal matrix through a nearly singular wrap and the correction's accuracy is not set by how singular the wrap is. At κ = 4·10⁷ one wrap returns the answer to 8.8·10⁻¹¹ — better than κ·u — and another, at κ = 3.8·10⁷, returns it to 3.3·10⁻⁶. The difference is how many of its samples sit near the symbol's zeros. A real wrap lands on a conjugate pair, a rank-two correction absorbs the pair exactly, and its two-by-two system has condition number 1.00 — which mattered because that system was solved by Cramer's rule; solved with pivoting, the single landing costs what the pair costs.

 10⁻³10⁻²10⁻¹110¹10²10³10⁴10⁵the shift the reduced Hessian neededshift taken ÷ shift needed10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹110¹median of eightexactly enoughκ(A) = 10, eight draws per curvatureworst ratio, curvature 10⁻⁸·²⁵1.8·10⁴worst ratio above 10⁻²7.3trials stopped below the curvature0factorisations a solve, mean3.2above one: more convex than the problemon the floor: the count passed a saddle The matrix a constraint makes

The shift that stops at the first right count

A nonconvex solver that finds the wrong inertia adds δI to H and tries again, and the δ it settles on is used as though it measured the curvature it corrects. It does not. On a well-conditioned constraint it is the schedule's number — 1.8·10⁴ times the need at a curvature of 5.6·10⁻⁹, between one and 7.3 times above 10⁻² — and on an ill-conditioned one the loop stops wherever the count first reads right: 63 of 152 saddles at κ(A) = 10⁸, the deepest with curvature 56. Refining the shift by bisection removes the first error and adds to the second.

 10⁻⁹10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²110¹10²10³10⁴standard deviation of the noise σmedian growth factorthe marginevery entrynonzeros onlyrelative, per entryone error in Estep 0.3, margin 0.00564, thirty drawsno noise1.1·10⁴every entry, σ = 10⁻⁵536nonzeros only, σ = 10⁻³1.1·10⁴one error in E, σ = 10⁻², draws kept21every curve is the same matrixonly what the noise touches changes Elimination, and the swap

Noise the growth amplifies

The shooting matrix's growth of 1.1·10⁴ fell under noise as its pivot margin divided by the noise, and the explanation offered was noise reversing partial pivoting's choices. But noise of 10⁻⁵ is five hundred times smaller than that margin. Add the same noise only to the entries that are not zero and every draw keeps the whole growth up to 10⁻³. The dense noise was not reversing the comparisons by itself: it sat in the zeros, the elimination multiplied it by the growth already made, and a comparison flips when that product reaches about four tenths of the margin — at every noise level from 10⁻⁶ to 10⁻³.

 10⁻¹101020304050z, the effective noise over thicknesstrials mirrored, %5 points, equal20 points, equal80 points, equal● five precise, 1/σ²■ five weighted ×100per cent mirrored at z = 0.45, 1,000 trialstwenty points1.8eighty points1.3five points9.5five precise, weighted 1/σ²7.9equal noise, five weighted ×1006.8both weighted schemes have five effective pointstwenty and eighty points share one curve Orthogonality, measured

Five precise points are five points

Weighting each sighting by its reliability is the standard form of an attitude or registration fit, and it changes how often the nearest orthogonal matrix comes back as a mirror. Measured, the rate is a function of two numbers: the weighted noise over thickness, and the effective count (Σw)²/Σw². Five points with a tenth of the noise, weighted by 1/σ², carry the information of 515 equal points and mirror like five — 7.9 per cent at a noise ratio where twenty points mirror 1.8 and five mirror 9.5. The √m the earlier measurement left unchecked is right, and it counts what carries the thin direction.

 10⁻¹⁴10⁻¹³10⁻¹²10⁻¹¹10⁻¹⁰10⁻⁹10⁻⁸10⁻⁷10⁻⁶050100150200250relative residual tolerance of every solveinner steps over twenty membersprevious answerline, last twoparabola, last threea straight path, twenty members, drift 0.02tolerance 10⁻⁶: previous answer62tolerance 10⁻⁶: line11tolerance 10⁻⁶: parabola17tolerance 10⁻¹⁰: previous answer160tolerance 10⁻¹⁰: line12tolerance 10⁻¹⁰: parabola14tolerance 10⁻¹⁴: previous answer261tolerance 10⁻¹⁴: line17tolerance 10⁻¹⁴: parabola23tighter to the lefteach point is a whole run of twenty When the problem arrives again

A straight path has nothing for a parabola to fit

The line through the last two answers beat the parabola through the last three on a drifting sequence, and the reason offered was that the parabola amplifies the stored answers' error. Tightening the solve tolerance from 10⁻⁶ to 10⁻¹⁴ should have reversed that, and it does not: the line needs 7 to 18 inner steps over twenty members at every tolerance and the parabola 14 to 23. The sequence's roots move along a straight line, so the line is exact and there is nothing else to fit. Bend the path by a part in ten thousand and the parabola wins below 10⁻⁹, by 105 steps to 75 at 10⁻¹⁴.

All themes