Numerical Linear Algebra
where the algebra stops being the arithmetic
Theme

The thread: Two routes to a number — page 5

Essays 97 to 120 of the 216 on this theme, in the same order.
10²10³10⁴10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴terms summedrelative error+∞, −∞ 1.01zero 1.00stochastic 0.50nearest 0.47√n against n, fittednearest, fitted exponent0.47stochastic, fitted exponent0.5toward +∞, fitted exponent1twelve seeds averaged at each sizethe slope is the bias, not the precision The arithmetic underneath

The direction the error leans

The size of one rounding error is set by the precision. How ten thousand of them combine is set by something else entirely — the rounding mode — and the fitted exponents are 0.47 for round-to-nearest and 1.01 for round-toward-infinity, on identical data at identical precision.

 11.31.61.92.210⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10²log₁₀ quadrature pointsdistance from the true counthalf an eigenvaluean integer, eventuallytrue count2at 4 points2at 128 points2finest error1.4·10⁻¹³the integral is an integerand a rounding hides how far it was Randomised, and the guarantee that changes kind

Counting what is inside a circle

A trace of a matrix nobody wants to form, integrated around a contour, gives an integer — how many eigenvalues are inside. It converges exponentially, it is estimated with random probes, and the probe block is a ceiling that the answer does not mention.

 12345678910111210⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹order of the reduced modelworst relative errora model made of measurementssamples used24degree read6gap at the cut2·10⁸best model, order12its error2.7·10⁻¹¹states in the original40filled: at its own samplesopen: everywhere else Reduction, and what a model is for

A model with no matrices behind it

Twenty-four numbers — values of a transfer function at twenty-four points — produce a sixth-order model of a forty-state system that passes through every sample to 10⁻¹² and matches the function it never saw to 10⁻⁹. Noise of 10⁻¹⁰ on those numbers takes the rank decision's gap from 10⁸ to 4.

 02468101210⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹singular value, in orderσ ⁄ σ₁eight digitsrenumberedin the separator's orderthe same cliff, from the other endσ₂ ⁄ σ₁0.062σ₄ ⁄ σ₁1.3·10⁻⁴σ₆ ⁄ σ₁1.5·10⁻⁸renumbered σ₄ ⁄ σ₁0.63rank at 10⁻⁸6the same entriesin two orders Sparsity, and what elimination costs

The cliff behind the count

The fill's rank is an integer between three and six across every separator two dense half-eliminations can afford, and this field has already recorded that a handful of such integers cannot carry a law. The singular values underneath are real numbers. They say the cliff's first step is 23.0 at a separator of eleven, 19.1 at fifteen and 16.2 at twenty-three — and that a control with no differential operator behind it gives 14,672.

 10²02040size niterationsno preconditionerwrapped (Strang)averaged (T. Chan)both are circulant approximations‖C − T‖/‖T‖, averaged0.22‖C − T‖/‖T‖, wrapped0.23smallest eigenvalue, wrapped, n = 160.23one of them is positive definiteand it is the one that is nearer Structure, and the solver that cannot see it

A speedup with a ceiling of its own

At ρ = 0.5 the averaged circulant takes 5 conjugate gradient steps at n = 512 against an unpreconditioned 30 — and that 30 is where the unpreconditioned count stops. It reads 29, 28, 30, 30 at n = 64 to 512 and then 29, 28, 26, 27, 25 at every doubling out to 16,384, because κ has reached 99.9% of Szegő's limit and the count has nothing left to grow with.

 10²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)‖x̂ − x‖ / ‖x‖modified, through Qᵀbclassical, either route‖QᵀQ − I‖, modifiedHouseholdermodified, on [A b]κ²uκuforward error at κ = 10⁸Householder4·10⁻⁹modified, through Qᵀb0.13modified, on [A b]2.7·10⁻¹⁰classical840×8, three seeds, against an exact rational solveQ is not orthogonal; x is right Orthogonality, measured

The right-hand side as one more column

Modified Gram–Schmidt's Q is 4.3·10⁻⁹ from orthogonal at κ = 10⁸, and a least-squares solve that multiplies b by it is wrong by 0.13. Hand the same routine b as an extra column instead and the answer is right to 2.7·10⁻¹⁰ — closer than Householder's 4.0·10⁻⁹. Classical Gram–Schmidt gains nothing from the same trick, to the last bit.

 110¹10²10⁵block size bwords movedbest block: b = 10the recursion: no block size126,742 wordsbest block, scanned1.1·10⁵the recursion1.3·10⁵recursion ÷ best1.2M = 144 words; the recursion never reads it1.16× the best of 24 blocks Where the flop count stopped predicting the time

The recursion that was never told the memory

A blocked elimination has to be tuned to its fast memory, and tuned to one memory it costs up to 2.9 times the best at another. A recursive elimination splits the columns in half down to one and reads no memory size at all. On eight fast memories from 36 to 576 words it moves between 0.94 and 1.28 times the words of the best tuned block, with the same 585,200 operations and the same pivots — and on a machine with two caches it beats the block tuned to either cache on six machines of seven.

 1611162126313610⁻²²10⁻¹⁹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹vertex, clique first then tail|entry| relative to the largestone rounding of the largest entrya bound that is provedPerron root11bracket, low11bracket, high11bracket width1.2·10⁻⁴smallest entry-1.5·10⁻²⁰entries below zero4every entry is positiveand the picture disagrees The matrix that is a graph

An eigenvector that must not change sign

Perron's theorem says the leading eigenvector of a connected nonnegative matrix is strictly positive. On a clique with a long tail, four of its thirty-six entries come back negative — and beside them is the only two-sided bound in these essays that is proved rather than estimated.

 10⁻⁴10⁻³10⁻²10⁻¹100.250.50.751anisotropy εsmoothing factor μpoint + fully-line + fullpoint + semi-y⅓, the one-dimensional answertwo routes, three curvesgap between the repairs2.2·10⁻¹⁶scan against closed form2.2·10⁻¹⁶the dashed curve lies on the solid one beneath itone repair, written two ways Iterating, instead of factorising

Smoothing a whole line at once

Solve every grid line in the strong direction exactly rather than sweeping over it, and the smoothing factor goes from 0.9993 back to 0.3340 — which is the one-dimensional answer, on a problem that is not one-dimensional. The repair replaces one ε in the closed form by a one.

 110¹10²10⁴10⁶10⁸10¹⁰10¹²n, points along one axismultiplicationsa dense factorisationthrough the eigenbasis6 decompositions of an n × nunknowns1.6·10⁴multiplications9.4·10⁵dense factorisation2.5·10¹²fitted exponent7‖Ax − b‖ ⁄ ‖b‖1.8·10⁻¹⁵nothing of size n^dis ever factorised When the index is a tuple

Five indices are cheaper than two

The same 4,096 unknowns cost 1.049·10⁶ multiplications indexed as a 64 × 64 grid and 1.966·10⁵ indexed as six axes of four. The dense factorisation that ignores the indexing costs 4.581·10¹⁰ at every one of them, and the residual improves in the same direction as the cost.

 condition numbers, on a logarithmic scaleκ(P)5.66·10¹⁸κ(Q)5.71·10¹⁸√(κ(P)κ(Q)) — the SVD route5.68·10¹⁸κ(P)κ(Q) — the product route3.23·10³⁷1/u4.5·10¹⁵the mean, or the productand only one of them fits Reduction, and what a model is for

A condition number that is not the model's

κ(P)κ(Q) is quoted as the reason one route to the Hankel singular values fails, and it is. It is also read as a measure of how reducible a model is, and it is not — the norm of the Gramian does not move at all as the McMillan degree runs from four to fourteen, and the condition number wanders over a factor of thirty-five with no trend.

 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.

 10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³relative compression error, the representation‖b − Ax‖ ⁄ ‖A‖‖x‖, the solveequala backward error, chosenslope0.99representation at 10⁻⁸1.2·10⁻⁹backward error there1.3·10⁻¹⁰representation at 10⁻¹²1.2·10⁻¹³backward error there2.5·10⁻¹⁴the compression is not an approximationit is a perturbation of the problem Two errors, and whose fault they are

An accuracy that is a backward error

Every backward error on this site is something an algorithm produced and somebody then measured. This one is a line in the program. Solving with a compressed matrix gives a residual that is the compression's own error, at a slope of 1.000 over ten decades, so the knob that sets the storage sets the backward error directly.

 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.

 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.

 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.

 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.

 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.

 point · strong in y0.951point · strong in x0.956y-line · strong in y0.037y-line · strong in x0.967semi-y · strong in y0.107semi-y · strong in x0.9661.00 — no convergenceone problem, seen from two sidesy-line, aimed0.037y-line, turned sideways0.97the ratio2631×31 grid, twelve V-cyclesthe direction is in the code, not in the problem Iterating, instead of factorising

Coarsening in one direction only

Leave the smoother alone and halve only the strong direction, and the smoothing factor is 0.3340 — identical to line relaxation's, at every anisotropy and every weight, to twelve digits. The convergence factors are then a factor of three apart, and at 45° both repairs fail outright.

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

All themes