Numerical Linear Algebra
where the algebra stops being the arithmetic
Theme

The thread: Measured, not assumed — page 10

Essays 217 to 240 of the 525 on this theme, in the same order.
10¹10²10⁻⁴10⁻³10⁻²10⁻¹1products with Arelative errorHutchinsonHutch++measured at equal costfitted rate, Hutchinson-0.5fitted rate, Hutch++-1.7error at 96, Hutchinson0.017error at 96, Hutch++0.0022both axes count products with Aso the sketch is paid for in the picture Randomised, and the guarantee that changes kind

A rate that belongs to the matrix

Hutchinson's fitted exponent sits near a half on every spectrum measured. Hutch++'s runs from −7.15 to −0.67 across the same four budgets, decided entirely by how fast the singular values fall — so one of the two methods has a convergence rate and the other has a rate per matrix. The ±1 probe's advantage moves the same way, from 1.56× at n = 10 to 1.09× at n = 120.

 110¹670673676679682685688691694pieces the inner products were summed initerations to the tolerance674, the cheapest run690, the dearestthe same solve, pricedpartitionings run13distinct counts11spread, per cent2.4best forward error6.3·10⁻¹¹worst1.1·10⁻¹⁰one matrix, one toleranceand the cost is the machine's The answer that depends on the machine

A stopping test is a race

One matrix, one right-hand side, one tolerance, thirteen partition counts — and eleven different iteration counts between 674 and 690. Every run converged, every answer is right to the accuracy asked for, and what differs is the bill.

 357911110¹grid side meigenvalue of C⁻¹Awithin ½ of onethe cluster grows like the sidem = 4: inside of 169m = 6: inside of 3611m = 8: inside of 6413m = 10: inside of 10017the cluster grows with the sideand the spectrum with the area Structure, and the solver that cannot see it

Four orders of conditioning, and four steps

On a 10×10 grid the two-dimensional kernel's condition number runs from 62 at ρ = 0.5 to 818,561 at ρ = 0.98. The preconditioned step count over the same range runs 18, 21, 21, 22, 21, 19, 18, and the count of eigenvalues the preconditioner actually brings within half a unit of one does not move at all — it is 9, 11, 13, 17 at every correlation the figure will draw.

 1357910⁻⁵10⁻⁴10⁻³10⁻²10⁻¹1eigenvalue λamplificationkeptdiscardeda filter on the starting vectormeasured against the polynomial3.9·10⁻¹¹worst kept direction1best discarded direction0.8the roots are the discarded Ritz valuesand the vertical lines are where they sit Eigenvalues, singular values, rank

Restarting is a filter

A restart throws away the Ritz values it does not want and begins again from a new starting vector. Written in the eigenbasis, that vector's components have been multiplied by a polynomial with its roots at the discarded values — measured component by component, and agreeing with the polynomial to rounding.

 unpreconditionedbest step20rule stops at7its error ÷ best1steps within 10%23truncated at τ = 0.032best step8rule stops at4its error ÷ best1steps within 10%6010203040506010⁻¹110¹steprelative errorrule stops: 7rule stops: 4plain CGLSpreconditionedbars: the steps within 10% of each run's bestthe rule reads the residual, not the window Methods that were designed apart

A stopping rule that follows the run it is given

A preconditioner that reaches the answer four times sooner leaves four steps within 10% of its best instead of sixteen, and a rule that stops by the residual ought to miss so narrow a window more often. Over forty draws of the noise it misses it less: the discrepancy principle stops at 1.030 times the preconditioned run's best against 1.073 times the plain run's. And past the edge it stops within a factor of 1.7 of a run whose own best is 5.5 times Tikhonov's — faithful to a run that has already failed.

 10141822263034384210⁻¹10⁻⁰.⁵110⁰.⁵grid points nrelative error against the continuous signalsampled kernel, hatsintegrated, hatsintegrated, spline0.1% noise, no λsampled kernel, hats, best0.13integrated, hats, best0.12integrated, spline, best0.12truncation of 64 points0.12same nodal unknowns, same databetter at representing the noise too Regularisation, and the answer that is chosen

A better discretisation is a weaker filter

A coarse grid's error has two sources — how well the discrete operator approximates the integral, and how well the grid's function represents the answer — and the grid essay could not separate them. Changed one at a time they separate: integrating the kernel against the hat functions takes a fifth off the 12-point error, reading the answer as a cubic spline takes 15 per cent more, and both roughly double the condition number on every grid. At 0.1% noise the spline discretisation's unregularised solve on 26 points reaches the best truncation of a 64-point grid to 0.3%. At 1% it is worse than the crude grid.

 6101418222610³10⁴10⁵grid side karithmeticND, totalMD, totalMD, critical pathND, critical pathdashed: the time on unbounded processorsthe same time, bought with more work Sparsity, and what elimination costs

An ordering that buys processors, not time

Nested dissection loses to minimum degree on fill and on total work at every grid either measurement could draw. Read along the elimination tree a parallel factorisation works on, it does not win back the time either: its critical path is within 28 per cent of minimum degree's at every size from 8 to 24 points a side, in both directions, and the tree heights and widest columns are nearly the same. What it wins is the ratio. Its total work over its critical path — the most a factorisation on unbounded processors can speed up by — grows from 2.7 to 4.5 while minimum degree's stays between 2.1 and 2.6.

 10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹110²10⁴standard deviation of the noise σmedian growth factormargin 0.00564shootingWilkinsonwhat the growth rests onWilkinson, no noise5.5·10¹¹Wilkinson, σ = 10⁻¹⁴2shooting, σ = 10⁻⁸1.1·10⁴pivot margin0.0056thirty draws at every pointa tie breaks at any noise; a margin needs its own size Elimination, and the swap

A worst case is as fragile as its margin

Wilkinson's matrix grows by 5.5·10¹¹ under partial pivoting, and adding Gaussian noise of 10⁻¹⁴ to every entry takes its median growth to exactly 2. The shooting matrix from a boundary-value problem grows by 1.1·10⁴, and noise ten million times larger leaves it untouched. The difference is what each worst case rests on. Wilkinson's rests on exact ties between candidate pivots, which any noise breaks. The shooting matrix's rests on a choice made by a margin of 5.6·10⁻³, and between 10⁻⁶ and 10⁻³ its median growth is that margin divided by the noise, times a constant between one half and four thirds.

 110¹123base case, columnswords ÷ the best block's√M − 2 = 10M = 144M = 144base of one1.2panels ≤ 10, worst1.2panels of 123the base case is a block size, rounded to a halvingand it has the block's cliff at √M − 2 Where the flop count stopped predicting the time

The block size a recursion still has

A recursive elimination is sold as having no block size, and every real one switches to plain loops below some width. Swept over that width, the traffic is a staircase with its steps at the halvings of n, and its cliff sits where the blocked elimination's does — the first panel wider than √M − 2 moves 1.53 to 4.47 times the words, on six memories of six. On three caches the innermost decides, and a third cache costs every tuned block up to 14 per cent and the recursion nothing.

 110¹10²10³10⁴10⁵10⁶10¹10²10³10⁴10⁵spread of resistances, largest ÷ smallest possiblecondition number, unit diagonalleast-resistance treebreadth-first treegreatest-resistance treenode equationsspread resistances make the loops easyand the nodes hard Orthogonality, measured

Spread resistances make the loops easy

Scaled to a unit diagonal, the loop equations on the least-resistance tree get easier as a network's resistances spread — from 120 to 5.44 over six decades — and stop depending on the grid's size, while the node equations of the same flow get harder, from 538 to 4.6·10⁴. The spread that ruins the range-space formulation rescues the null-space one, though the loops' density means the work saved is a factor of two, not the factor of nine the iteration counts suggest.

 01234567910δ added to Hpositive pivots−λmin(H) = 4.67|μ| = 1minimumsaddlewhat the signs are countingminimum, true count10saddle, true count9saddle read as minimum from1.1guarantee needs δ past5.1the signs count ZᵀHZ + δInot the curvature the problem has The matrix a constraint makes

A shift that certifies a saddle

On a constrained problem whose Hessian has four negative eigenvalues, a saddle-point matrix is quasi-definite only once H + δI is positive definite — past δ = 5.08 here. Its pivot signs then count the curvature of ZᵀHZ + δI rather than of ZᵀHZ, so a saddle with a negative curvature of −1 is certified a minimum from δ = 1.05 on, and every saddle shallower than δ goes the same way. Iterative refinement against the unregularised matrix keeps the second-order test the count gave up: it contracts on the minimum at δ/(μ + δ), 0.980 a step at δ = 5, and on the saddle it grows at exactly 1.25.

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

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.

 4914192429343900.20.40.60.81decimal digits of the numbers the lattice is givenshare of trials recovering the relationa double has about this many digitsnumbers that are exactthe same numbers, measureda window with two edgescoefficients up to30exact, at 12 digits1exact, at 401doubles, at 121doubles, at 250doubles, at 400.13one edge is the relationand the other is the data Exact arithmetic, and what it costs instead

A relation among digits that were not there

A lattice finds the exact integer combination of several numbers that vanishes, and it needs enough digits of them to do it. Give it more digits than the numbers have and it finds a relation anyway — among the rounding, with the same confidence and no warning. Recovery runs 8 of 8 at twelve digits and 1 of 8 at forty.

 10⁻²10⁻¹10¹10²drift between consecutive membersinner steps over the whole runthe previous answer — a warm startthe parabola through the last threethe line through the last twodegree zero is a choicemembers per factorisation4warm, at drift 0.02160linear, at 0.0212quadratic, at 0.0226the factor13drifts where warm fails2one more term in the extrapolationand a factor of thirteen When the problem arrives again

A warm start is degree zero

The previous answer used as the next member's starting point costs 160 inner steps over twenty members. The line through the last two answers costs 12 — a factor of thirteen, for three vector operations and no extra storage. The parabola through the last three costs 26, which is worse than the line and better than the point.

 12345610⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹Ritz value, largest firstdistance from 10one vectora block of twoa space, not a ratecopies found, one vector1copies found, block of two2Krylov dimension, one vector16the second copy is not in the spaceat any number of steps Eigenvalues, singular values, rank

An eigenvalue one vector cannot see

A matrix with an exactly doubled eigenvalue at 10. Twelve Lanczos steps find it once; twenty-four find it once, on a Krylov space of dimension 23 in a 24-dimensional problem. A block of two vectors finds it twice. This is not slow convergence — the second copy is not in the space.

 half a stepthe 32 values of one block, in the order they arrivestep 0.0625one scale, thirty-two valuesoctaves inside the block1.9entries rounded to zero0worst error over its bound131 levels either side of zerothe largest entry chose the step The arithmetic underneath

One exponent for thirty-two numbers

Share the exponent across a block and the cost per value drops from eight bits to 6.25, and the accuracy improves — up to about three octaves of spread inside a block. Past that a single outlier deletes the thirty-one values beside it, and the 2-norm barely notices.

 1357911131510⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1vertexrelative error against the exact answerthree routes, one numberworst, by elimination4.8·10⁻¹⁶worst, by spectrum10·10⁻¹⁶solve residual1.5·10⁻¹⁵Foster's sum15n − 115difference3.6·10⁻¹⁵the answer is a ratio of integersso the error is known The matrix that is a graph

A distance computed by a solve

Effective resistance is the one quantity in this field with no combinatorial route to it — it is defined by a linear system. On a small unweighted graph the answer is a ratio of two integers, so for once the error is known rather than estimated, and every resistance in a graph has to add up to a number fixed in advance.

 048121610⁻⁸10⁻⁷10⁻⁶10⁻⁵extra columns in the sample, prelative compression errorthe best rank-8 representationfive draws, mean and rangea band, not a lineexcess at p = 012excess at p = 83spread at p = 00.44spread at p = 80.12the optimum of this rank7.8·10⁻⁸one seed shows the meanand five show the risk Randomised, and the guarantee that changes kind

What a single draw cannot report

With no oversampling the construction's error is 11.6 to 50 times the best representation of the same rank, and the spread across five seeds runs from 16 to 146 per cent of the mean. Sixteen extra columns bring the excess to between 2.0 and 2.9 at every rank measured and the spread to between 4 and 20 per cent — and the second number is the one a single run cannot report and the one that decides whether the first is a measurement.

 A6×6B6×6I ⊗ A + Bᵀ ⊗ I36×36one equation, two objectsentries in A and B72entries in the coefficient matrix1296two routes, relative gap1.7·10⁻¹⁶‖AX + XB − C‖/‖C‖1.4·10⁻¹⁶the small squares are the problemand the large one is the notation Structure, and the solver that cannot see it

The elimination the matrix does not need

The Kronecker form of AX + XB = C is dismissed with a hundred million entries and (2/3)n⁶ operations. Both price an elimination, and after the reduction both routes take, the matrix has exactly n³ nonzeros, none of them above the block diagonal, and nothing left to eliminate.

 234567891010⁻⁴10⁻²1ADI steps krelative error, and max|r|²max|r|²the measured errora bound that is an identitysteps drawn9ratio, worst0.77spread of the ratio1error at k = 103.9·10⁻⁴its residual1.1·10⁻⁴four orders of errorand one constant Reduction, and what a model is for

Bracketing an error nobody can measure

The error of a low-rank Gramian factor is the one quantity a caller cannot compute, because computing it needs the Gramian the factor exists to avoid forming. Two numbers that can be computed sit either side of it — a rational factor known before the run, and a residual known after — and they stay a factor of four apart across a fourfold change of size.

 123456024681012moments Krank of the block Hankel12 eigenvalues insidethe ceiling K·ℓthe ranka ceiling with a knob on itprobes2eigenvalues inside12rank at K = 12rank at K = 612solves, at every K512one ceiling per probeand K of them per moment The eigenvalue problem that is not linear

A ceiling with a knob on it

A contour method returns at most as many eigenvalues as its probe block has columns, and the object that comes back does not distinguish that from having found everything. One line of the derivation multiplies the ceiling by a number the caller chooses, and it costs no extra solves at all.

 110¹10²10⁻¹²10⁻⁸10⁻⁴110⁴ncondition number, and residual reachedmarks above: κ of the step · dashes: its closed formmiddle: a fit from a random startbelow: the same fit started at the answera decomposition that is ill-conditionedκ at n = 1283.3·10⁴its closed form3.3·10⁴cosine of the terms1from a random start0.012from the answer1.4·10⁻¹⁰the answer existsand cannot be found Two errors, and whose fault they are

A tensor that cannot be decomposed

Every member of a certain sequence is exactly a sum of two rank-one terms, and both terms are written down in closed form. A three-hundred-sweep fit from a random start does not find them, and stalls at the same one per cent however far the sequence goes — while a fit started at the answer loses digits exactly as 2n² says it should.

 -13-11-9-7-5-3-102468log₁₀ of the accuracy asked forlargest train rank0.30 of rank per decadea cost that is typed inentries216stored at 10⁻²90stored at 10⁻¹²288rank per decade0.3worst error ⁄ tolerance0.83the storage is chosena constant of rank a decade When the index is a tuple

The digit that costs more than the tensor

Ask a three-index reciprocal tensor on six points a side for seven digits and its train is 288 numbers against 216 entries. The break-even rank is n − 1 at all four grids measured, and a train that reaches it fits with exactly n numbers to spare.

 10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹relative compression error, chosen‖x − x*‖ ⁄ ‖x*‖κ × the chosen backward errorthe error in the answerboth factors known firstκ21chosen at 10⁻⁸1.4·10⁻⁹predicted forward2.9·10⁻⁸measured forward1.6·10⁻⁹bound ⁄ measured26the amplifier, used forwardsfor once Neither sparse nor dense

The knob that moved two things

Decide how many digits the answer needs, divide by the condition number, and compress to that. It is the one rule licensed in advance here, and its two factors are not the independent inputs it reads as: the partition's leaf moves neither of them and moves the answer by nearly a factor of three, and the only knob here that raises κ halves the ranks while it does so.

All themes