Numerical Linear Algebra
where the algebra stops being the arithmetic
Theme

The thread: Measured, not assumed — page 12

Essays 265 to 288 of the 545 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.

 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 Orthogonality, measured

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.

 all 2^16 subsets searchedboth minima, one order0.97minimum degree optimal0.82worst fill overshoot1.1worst work overshoot1.111.021.041.0611.041.081.121.16fill ÷ the least filloperations ÷ the leastovershoot squaredon the diagonal: the two overshoots are equalon the upper curve: the arithmetic overshoot is the square Sparsity, and what elimination costs

Two minima that are one minimum

The order that decides the memory found the operation count behaving like the square of the fill, which leaves room for an order with slightly more fill but a shorter heaviest column to do less arithmetic. Searched exactly over every elimination order on forty graphs, that order does not exist: one order attains both minima on thirty-nine of forty, and on the fortieth the least-fill order's arithmetic is 1.0099 times the least. Minimum degree attains both on the same thirty-three graphs and neither on the same seven.

 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.

 01234567110¹10²levelstored entries per rowalgebraicgeometrictwo complexities, one hierarchygrid complexity3.1operator complexity18geometric, for comparison1.5a ring plus random chords — nothing is anywherelevel 6 is 100% dense Iterating, instead of factorising

A hierarchy with no grid behind it

On a graph Laplacian the algebraic V-cycle converges at 0.199 a cycle, its grid complexity is an unremarkable 3.05, and its operator complexity is 17.7 — one level of forty-one unknowns is entirely dense. The number people quote is the one that does not measure the work.

 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.

 5860626466687010⁻¹⁷10⁻¹⁶terms in the dot productmean relative errorthe kernel changes herea constant in a libraryone accumulator3.2·10⁻¹⁷four accumulators2.1·10⁻¹⁷step at the cutoff1.6cutoff64the problem did not changethe loop did The answer that depends on the machine

The length that changes the kernel

A dot product's accuracy steps by a factor of 1.57 between 63 and 64 terms, on vectors drawn identically at both lengths. Nothing about the problem changes there. A library switches from one accumulator to four, at a constant in somebody else's source file.

 05101520253035059118177236295vertices eliminatededges of fill so farminDegree: 71natural: 125reverse: 125random: 160maxDegree: 293fill, by orderingminDegree71natural125reverse125random160maxDegree293edges to start60eliminating a vertex makes a cliqueand the order decides how big The matrix that is a graph

Eliminating a vertex is a graph operation

Gaussian elimination on a Laplacian deletes a vertex and joins its neighbours into a clique with conductances wᵢwⱼ over Σw. The matrix that remains is still a graph — symmetric, zero row sums, nonpositive off the diagonal — and the ordering decides whether the fill is thirty-one edges or four hundred and sixty-five.

 12345678910111210⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹degree of the approximantworst relative error on the target setchosen by the residualdegree reached12error there10⁻¹²even support, same degree1.1·10⁻⁷advantage702nearest support point0.096clustering ratio1.5filled: points chosen by the erroropen: points spread evenly The eigenvalue problem that is not linear

The points the algorithm chose

A rational approximant whose support points are picked by its own residual clusters geometrically at a branch point nobody named — recovering by measurement the rule a hand-built approximant is given. At degree ten it is seven hundred and fifty times more accurate than the same form with its points spread evenly.

 ‖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.

 12345678110¹partial orders kept at each pairingexcess over the exhaustive orderthe exhaustive orderdashes above: each rule's worst casepair by cheapest productpair by smallest resultone knob, two rulesby cost, width 11.5by cost, width 21.1by size, width 11both, width 81by size, worse wider24by cost, worse wider3the width repairs the bad ruleand does not improve the good one Where the flop count stopped predicting the time

The search that got worse as it widened

Keeping two candidate orders instead of one removes a third of the cheapest-product rule's excess, and by eight it has removed all of it — the two greedy rules this field separated become the same rule. On twenty-four of sixty networks a wider search returns a dearer order than a narrower one, and it is the better rule that it more often makes worse.

 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 λ When the index is a tuple

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.

 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.

 unpreconditioneddimension a step1.2best iterate carries24at step20τ = 10⁻³, past the edgedimension a step3.6best iterate carries33its error0.4837111519232731353910⁻¹10⁻⁰.⁵1effective dimensionrelative errorthe arc: dimension 8.0 to 24.0none10⁻⁰.⁵10⁻¹10⁻¹.⁵10⁻²10⁻².⁵10⁻³large dots: each run's own best iteratethe answer sits at one effective dimension Methods that were designed apart

The parameter neither knob is

A preconditioned run has a cutoff and a step count, and neither is the regularisation parameter. The parameter is the effective dimension of the iterate: every cutoff that works puts its own best at 23.7 to 24.3 of it, where the unpreconditioned run's best sits at 23.2, and what the cutoff buys is the rate — 1.27 of it a step with no preconditioner and 3.53 with one. The edge is where a single stride is longer than the distance left.

 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 Regularisation, and the answer that is chosen

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.

 an ordering that is theregreedy at n = 764best ordering at n = 72its largest multiplier1456710¹10²matrix size ngrowth factorgreedy, ties forwardsgreedy, ties backwards, ε = 10⁻¹²greedy, ties backwards, ε = 0the best orderingthe best order is 2 3 4 5 6 1 at n = 6one cyclic shift of the rows Elimination, and the swap

The order the greedy rule cannot choose

Wilkinson's matrix is the standard demonstration that partial pivoting's growth bound of 2^(n−1) is attained. It is attained by the row order the greedy rule picks, and not by the matrix: a single cyclic shift of the rows gives growth 2 at every size, with no multiplier above one. At n = 7 that is 64 against 2. And perturbing one entry by 10⁻¹² leaves the good order exactly where it was while putting every tie-break of the greedy rule back on 64.

 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 Orthogonality, measured

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.

 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.

 40 drawsproducts at target 0.01696deflated277draws inside the target0.810⁻²10⁻¹10⁻²10⁻¹target standard errorrelative error reachedworst draw, no deflationmedian error, no deflationmedian error, deflatedthe grey diagonal is a calibrated rulethe medians are on it and the worst draws are not Randomised, and the guarantee that changes kind

A rule that reads only its own probes

A trace estimator is a mean of independent samples, so its own standard error is estimable from the samples and a stopping rule needs nothing the estimator does not already have. Over forty draws it is calibrated in the middle and not at the edge: at a target relative standard error of 1% the median error reached is 5.3·10⁻³ and the worst of forty is 2.9·10⁻² — three times the target. And the cost of the target is the estimator's own square root: tightening it from 3% to 1% takes the median probe count from 75 to 696.

All themes