Numerical Linear Algebra
where the algebra stops being the arithmetic
Theme

The thread: Measured, not assumed — page 11

Essays 241 to 264 of the 525 on this theme, in the same order.
-12-10-8-610⁻¹³10⁻¹¹10⁻⁹10⁻⁷log₁₀ of the residual tolerance asked forforward error of the answer1.34×1.48×1.71×1.17×the band does not closetolerances swept4runs at each7accuracy gained1.5·10⁶ratio at 10⁻⁶1.3ratio at 10⁻¹²1.2the bars falland they keep their height The answer that depends on the machine

The tolerance that buys no agreement

Ask for four more orders of accuracy and you get them — the answers improve by a factor of 1.5 million. The ratio between the best and the worst run is 1.34, 1.48, 1.71 and 1.17 across the same sweep. The band falls and it does not close.

 unpreconditionedbest step20best error0.14Tikhonov's best0.14cosine, truncatedbest step5best error0.14directions preconditioned25010203040506010⁻¹110¹steprelative errorTikhonov's best: 0.1405plain CGLSpreconditionedthe blur approximated by a cosine transformtruncation leaves it alone Methods that were designed apart

What a cheap preconditioner has to leave alone

A blur approximated by a matrix the cosine transform diagonalises agrees with the operator everywhere but its first and last seven rows. Made invertible by a shift, as the exact preconditioner was, it never reaches the unpreconditioned run's floor — at α = 10⁻³ its best iterate is 0.749 against 0.143. Made invertible by leaving every eigenvalue below τ alone, it reaches 0.141 in five steps instead of twenty, and the smallest τ that keeps the floor sits at a third to a half of the Tikhonov oracle's λ at three noise levels.

 122028364452606876849210010⁻¹110¹grid points nrelative error against the continuous signalbest with no λ: n = 26no λeach grid's best λTikhonov, 64 points0.1% noise, 16 drawsbest grid with no λ26its error0.1396 points with its λ0.12best λ on fine grids0.0056each point is a median over the same drawsλ belongs to the problem, not to the grid Regularisation, and the answer that is chosen

Where the grid hands over to λ

An unregularised solve on a coarse grid comes within a tenth of the best Tikhonov answer on a fine one, and the pair of a grid and a λ was left unmeasured. Measured, the two do not trade. On every grid up to the best unregularised one no λ helps at all. On every grid of 40 points and more the best λ is the same to within a quarter of a decade — 3.2·10⁻² at 1% noise per sample, 10⁻³ at 0.01% — and the 96-point grid with it beats the best coarse grid by 7, 9 and 13 per cent. The grids between the two, given their own λ, land between them.

 least possible76minimum degree76worst of 300 tie-breaks81reverse Cuthill–McKee85natural order9920 unknownstie-breaks at the least0.44minimum degree ÷ least1every order searched, by the set already eliminatedthe greedy rule found the minimum Sparsity, and what elimination costs

The least fill there is

Finding the elimination order with the least fill is NP-hard, and that is a statement about the hardest graph and the largest size. On a graph of twenty vertices every one of the 20! orders can be searched at once, through the million sets of vertices already eliminated, and the least fill is a number. On the 4×4, 4×5 and 3×7 grids minimum degree finds it exactly. On eighty random sparse graphs of eighteen vertices it finds it on 53 and misses by at most 7.6 per cent, and on every one of the eighty some breaking of its ties finds it.

 groundarcs of the treearcs off the treethe loop worst servedleast-resistanceκ(Z)11κ(ZᵀHZ)143κ, unit diagonal9.3nonzeros in Z426longest loop, arcs26worst path ÷ own arc0.97every tree gives a basis of whole numbersthe resistances decide which one to want Orthogonality, measured

The tree the resistances choose

On a network every basic set is a spanning tree and every null-space basis is a set of loops with entries 0 and ±1, so no tree can make Z badly conditioned. The tree with the best-conditioned Z still gives loop equations 4.8 times worse than the tree of least resistance: the basis has to be chosen against the Hessian, and pivoting finds it only when it pivots on the resistances too.

 -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 Least squares, and the road not to take

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.

 0246810121416110¹10²10³10⁴10⁵10⁶length of the interval Tgrowth factorpartialrookcompleteh = 0.3, 102 unknownspartial, T = 15.01.3·10⁵e^(5T/6)/21.3·10⁵complete2largest κ8.3a boundary-value problem, not a constructionκ single-digit throughout Elimination, and the swap

The growth a boundary-value problem supplies

Large growth under partial pivoting is usually said to need a matrix built for it. A two-point boundary-value problem solved by multiple shooting supplies one without being asked: its growth factor is e^(5T/6)/2 to four figures — 1.1·10⁴ at an interval of 12, 1.3·10⁵ at 15 — on a matrix whose condition number never exceeds 8.3, while rook and complete pivoting keep it below 2. And it is the finer shooting grid that grows: below a step of 0.3397 the choice partial pivoting makes turns on one entry against one, and above it there is no growth at all.

 -18-15-12-9-6-300log₁₀ ‖PKPᵀ − LDLᵀ‖ / ‖K‖share of orderingsbestworstexistence and stabilityfactorise1unregularised0.67worst growth5·10⁷growth × γ0.5the ordering is free to chooseand not free of consequence The matrix a constraint makes

The perturbation that does the work

A saddle-point matrix made quasi-definite is perturbed in both blocks, and the laws measured for it moved both together. Moved apart, the laws all belong to one block. The zero block's perturbation γ decides whether every ordering factorises, sets the worst ordering's growth at 0.51/γ, and costs the answer 1,451 per unit — the reciprocal of the smallest eigenvalue of AH⁻¹Aᵀ to three figures. The perturbation of H moves none of the first two and costs 19 per unit. Refinement removes each block's perturbation at the rate its own Schur complement sets, so γ's limit sits fifty times nearer than δ's.

 10⁴1ceiling on the largest intermediate, numbers heldarithmetic ÷ the unconstrained cheapest orderno evaluation order fits below thisminimising the peak insteadthe unconstrained cheapest ordera ceiling, not a targetfree order, peak9504leanest possible2106cheapest that fits1leanest order costs1.6the difference1.6largest input2106the same memoryat two prices Where the flop count stopped predicting the time

A ceiling is not a target

Asked to hold the smallest possible intermediate, an evaluation order costs a median of 1.55 times the cheapest order's arithmetic. Asked to hold no more than that same amount, it costs 1.15. The two answers hold exactly the same number of numbers, and on one network they are 5.17 times apart in work.

 10¹110¹skew of the starting basisdistance ÷ the distance to the nearest pointthe nearest lattice pointdashes: the skewed basis's orthogonality defectrounding in the skewed basisrounding in the reduced basisthe same three lines of arithmeticskew40defect there40worst, skewed40mean, skewed26mean, reduced1exact share, skewed0one lattice, one targetand two descriptions Exact arithmetic, and what it costs instead

Rounding a coordinate in the wrong basis

The nearest lattice point is found by writing the target in the basis and rounding each coordinate, which is three lines of arithmetic and is wrong. On a basis skewed by forty it lands a mean of 26 times too far and a worst of 40.02 — the basis's own orthogonality defect, to three figures — and the identical three lines on the reduced basis are exact at every target.

 0204060801001share of shufflings, per centarithmetic ÷ the sorted order'sthe sorted order, and the nearest-neighbour paththe order is a free variabledrift between neighbours0.04median penalty1.5worst penalty1.8factorisations, sorted4worst shuffle's16greedy ÷ sorted1the same sixteen problemsat two prices When the problem arrives again

The order a batch arrives in

Sixteen problems over a parameter, solved in the order the loop produced them, cost a median of 1.54 times what the same sixteen cost sorted, and 2.80 times at the worst shuffling. A nearest-neighbour path computed from the parameter values alone recovers the sorted cost exactly, at every drift and every shuffle.

 the grid the operator came froman edge is a coupling the matrix calls strongkeptinterpolatedwhat the entries decidedcoupling ratio, x against y0.001strong couplings across x0strong couplings along y110rows kept or dropped whole1111×11 grid, θ = 0.25the coarse grid, from the matrix alone Iterating, instead of factorising

The coarse grid the matrix chooses

Given a tridiagonal matrix and no information about a grid, the coarsening keeps every other point and derives the weights ½, 1, ½ — the operators the geometric method was handed. Given the anisotropic operator, it discovers semi-coarsening, in the right direction, without a coordinate.

 13579111310⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹index ksingular valueσ₁ · 10⁻¹⁴, the thresholdone matrix, three rankspartitionings7lowest rank10highest rank12threshold7.3·10⁻¹⁴σ₁7.3the curves separate in the noiseand the threshold is drawn through it The answer that depends on the machine

A rank that depends on the thread count

One 60 × 14 matrix, one threshold, seven partitionings of the inner products that build its Gram matrix — and numerical ranks of 12, 12, 12, 10, 10, 11 and 11. Not a digit of an answer: the number of columns a model built from this matrix would have.

 19172533414911.251.51.752indexsample ÷ originalkept, and not keptedges kept344of1225eigenvalue ratio, low0.43high1.7degree ratio, low0.56high1.5diameter, after3filled dots are eigenvaluesopen dots are degrees The matrix that is a graph

A graph with a tenth of the edges

Keeping 344 of 1,225 edges, sampled by effective resistance and reweighted, preserves every eigenvalue of the Laplacian to within a factor of 1.7. It preserves no degree — half of them are wrong by more than a third — and it takes the diameter from one to three.

 56789100108216324432log₂ nstored per unknowndenseweak: every off-diagonal blockstrong: only the admissible onesa constant per doublingstrong, n = 5126.8·10⁴weak, n = 5126.1·10⁴dense, n = 5122.6·10⁵per doubling26relative compression error3.8·10⁻¹⁰the dense line doublesand the other two add a constant Neither sparse nor dense

The offset that moved the slope

The accuracy is supposed to lift a storage curve and leave its growth alone. Measured at seven tolerances, the strong partition adds 10.9 numbers per unknown per doubling at two digits and 36.5 at twelve — the growth rate more than triples, so ten decades of accuracy cost 33 per cent more storage at 64 unknowns and 118 per cent at 512.

 10⁻³10⁻²10⁻¹110⁻¹⁹10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹1 − ρbackward errorLevinson, Toeplitzelimination, ToeplitzLevinson, any kindelimination, any kindat ρ = 0.999Levinson, any perturbation2.9·10⁻¹⁷Levinson, Toeplitz only6·10⁻¹²the ratio between them2·10⁵diagonal defect of the first0.97the number that is reportedand the number that is asked about Structure, and the solver that cannot see it

The number that cannot rank them

Levinson and Gaussian elimination are indistinguishable on the backward error a library reports — every one of ninety-six measurements between 1.16·10⁻¹⁷ and 5.73·10⁻¹⁷. The structured backward error separates them by up to a hundredfold, in whichever direction the point happens to give. Only the forward error ranks them, and only because this family's exact answer is known.

 0246810121410¹⁵10¹⁹10²³10²⁷10³¹10³⁵10³⁹10⁴³k|CAᵏB|6 exact zerosa definition that will not computestates16exact zeros6‖A‖₂, measured1146last step-to-step ratio1330range across the run2·10²⁹zero for the travel timethen ‖A‖₂ᵏ Reduction, and what a model is for

The definition asks for more of what defeats it

The rank of a p × p Hankel matrix of Markov parameters resolves a degree of p, and it needs 2p parameters to do it. Those parameters grow like the norm of the state matrix raised to their index, so the count that buys resolution is the same count that buys dynamic range. One model, four run lengths, and a spread that runs from 10¹¹ to 10¹⁰².

 entries rounded to zero, of 320block, as given310block, sorted by size22E4M3, either order0median entry's relative errorblock, as given1block, sorted0.00321E4M30.022the same numbers, three waysdeleted, as given310deleted, sorted22deleted, per-element08-bit significands, blocks of 32sorting is free and changes no value The arithmetic underneath

A bit buys an octave

The outlier a block survives is exactly two raised to its significand width — 8 at three bits, 32 at five, 128 at seven, 512 at nine. Each extra bit doubles the range the block tolerates and halves the ordinary entry's error. Reordering the same numbers buys every octave at once and costs nothing.

 1234567891010⁻⁵⁵10⁻⁴⁹10⁻⁴³10⁻³⁷10⁻³¹10⁻²⁵10⁻¹⁹10⁻¹³10⁻⁷10⁻¹cyclesizethe residualrepairedthe reported boundwhat the stopping rule readsreported at the last cycle9.4·10⁻⁴¹the residual there5.7·10⁻⁵with the border recomputed3.8·10⁻¹⁴products, cheap and repaired9a bound with nothing under itand one product a cycle to fix it Eigenvalues, singular values, rank

Keeping the vectors, and losing the bound

Thick restarting keeps the Ritz vectors instead of filtering the starting vector — the same eigenvalues for a third of the products with A. Its residual bound reaches 9.4·10⁻⁴¹ while the residual it bounds sits at 5.7·10⁻⁵, and the eigenvalues are correct to 4.3·10⁻¹⁴ the whole time, so nothing reports it.

 123456110²10⁴10⁶moments Kκ of the block Hankelintegrated in zin (z − c)/ρone division per pointradius6unscaled, K = 62.3·10⁶scaled, K = 61750the factor between1291the ceiling rises by Kand so did the conditioning The eigenvalue problem that is not linear

The conditioning that rises with the ceiling

Higher moments multiply a contour method's ceiling by K and grade its block Hankel over ρ to the 2K, so the two knobs are the same knob. One division per quadrature point separates them, and the measurement of what it is worth grows from twenty to twenty thousand.

 -12-10-8-6-4-2010⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹10¹⁵log₁₀ μcondition number, and relative errorκ₂, augmentedcomponentwise, condensedcomponentwise, augmentederror, augmentedone matrix, two numbersκ₂ at μ = 10⁻¹²3·10¹³componentwise, same matrix13their ratio2.3·10¹²measured error9.4·10⁻¹⁶every library prints the top lineand the error obeys the third Two errors, and whose fault they are

Two condition numbers of one matrix

κ₂ is a worst case over perturbations of a given norm, and a normwise perturbation may put its whole budget on the smallest entry. The componentwise number is a worst case over perturbations proportional to the entries, which is what a backward-stable factorisation actually makes. On one matrix they are 3·10¹³ and 13.3, and the error obeys the second.

 10¹10²110¹rank of the trains at the moment of evaluationcost ÷ the cost of this rank's own best ordereach rank's own best orderthe order compiled at rank 4, carriedthe smallest-result rule, recomputedthe cheapest-product rule, recomputedone plan, eight rankscompiled at rank4worst carried8carried at rank 2568fresh rule at 2562.9carried at rank 21fresh rule at 21.1the search was exhaustiveand the dimensions moved Where the flop count stopped predicting the time

The plan that was right at rank four

An evaluation order is chosen once and paid for thousands of times, and the dimensions it was chosen at are not the dimensions it runs at. Compiled at rank four and run at rank 256 it costs 8.01 times the order that rank deserves; compiled at 256 and run at 2 it costs 301 times. A one-line rule recomputed on arrival costs 2.92 and 1.11.

 01210⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹terms fitted1 − the largest cosine between two fitted terms3 — the tensor's rank45below this line, two terms are the same termthe residual says nothing3 terms, collinearity0.524 terms, collinearity15 terms, collinearity13 terms, error2.1·10⁻¹⁴4 terms, error10⁻¹³4 terms, sweeps76one term too manyand two of them cancel When the index is a tuple

One term too many

A tensor built from three rank-one terms, fitted with four, reaches a relative error of 1.0·10⁻¹³ in seventy-six sweeps. Two of the four terms come back with a cosine of 0.99927 between them, subtracting each other, and the smallest weight is a fifth of the largest rather than a rounding. Nothing in the residual says any of it.

 10¹10³10⁵10⁷10⁹10¹¹10¹³10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹κ(B), the conditioning of the constraint blockrelative error, and feasibilitydashes: the method of weightingthe saddle-point routethe null-space routealong the bottom: how nearly every answer satisfies the constraintsfeasible and wrongκ(B)4.6·10¹²κ(B)·u0.001null space1.2·10⁻⁴saddle point0.0082weighting2.6feasibility10⁻¹⁵every constraint is satisfiedand the answer has no digits Least squares, and the road not to take

Feasible and wrong

A third constraint that nearly repeats the first takes the best route's answer from 2.96·10⁻¹⁵ to 1.16·10⁻⁴, and the other two routes to no correct digit at all. Every one of those answers satisfies every constraint to 10⁻¹⁵. The quantity a caller checks after a constrained solve is the one quantity here that says nothing.

All themes