Numerical Linear Algebra
where the algebra stops being the arithmetic
Theme

The thread: Measured, not assumed — page 9

Essays 193 to 216 of the 505 on this theme, in the same order.
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.

 10⁻¹⁶10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²110⁻⁷10⁻⁵10⁻³10⁻¹10¹relative error in the datarelative error in the answerthe answer's errorthe data's errornothing here is the arithmetic'sdata error10⁻¹²answer error10·10⁻⁴amplification10·10⁸rounding committed0residual0the residual is exactly zeroand the answer is wrong anyway Exact arithmetic, and what it costs instead

An exact answer to a measured problem

The residual is the zero vector, nothing was rounded at any step, and the answer is wrong in its first digit. Data accurate to fourteen places, an exact solve of the system it defines, and an error of 10⁻⁵ — because conditioning was never a statement about arithmetic and removing the arithmetic error removes none of it.

 10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²ε the blocks were compressed at‖A − LLᵀ‖ ⁄ ‖A‖the factorisationthe knob, on a factorisationtruncations34residual at 10⁻⁴2.2·10⁻⁵residual at 10⁻¹⁰9.9·10⁻¹²slope1worst ratio to the representation1.8ten decadesand a slope of one Neither sparse nor dense

A knob calibrated in residuals

A formatted Cholesky has two numbers in it and only one of them is an accuracy. Across twelve trees — three sizes by four leaf sizes — the leaf moves the truncation count from 0 to 258 and moves the ranks of the blocks not at all, while the residual follows the tolerance at slopes between 1.022 and 1.046 and sits at about a tenth of it throughout.

 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.

 15913172110¹10⁶10¹¹10¹⁶10²¹10²⁶number of indicesnumbers storedthe tensor: n^dthe core: r^d + dnrthe train: (d − 2)nr² + 2nrthe curse, movedentries at d = 2010²⁶core1.1·10¹²train5920core ⁄ train1.9·10⁸entries ⁄ core9.5·10¹³the definition is repairedthe size is not When the index is a tuple

A compression of 10¹⁴ that still does not fit

A Tucker core of a twenty-index array at rank four is 1.1·10¹² numbers against the tensor's 1.05·10²⁶ — a compression by a factor of 9.5·10¹³ that is still nearly nine terabytes. The ratio is not the verdict. The verdict is a ceiling, and the ceiling is a number of indices.

 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.

 1110¹10²10³10⁴10⁵the noise level it is told ÷ the true oneerror ÷ the oracle's0.40.50.71.523twice the oraclecliff, ρ = 0.69told the truthworst of 400middle 80%median draweach draw against its own oraclecliff, median ρ0.69told the truth, median1.1told a half, median3201shaded: the middle 80% of drawsbelow the cliff the error has doubled Regularisation, and the answer that is chosen

Thirty-two coefficients instead of a noise level

The discrepancy principle has to be told the noise, and told too little it does not degrade — it falls off a cliff, at 0.80 of the truth when the noise is 10% and at 0.58 when it is 0.001%, exactly where the understatement forces the filter past its best truncation. The missing number is in the data. The root mean square of the last thirty-two coefficients never sends the rule over the cliff at or below 1% noise in four hundred draws, where eight coefficients with the same median do so thirty-five times.

 12345610⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹digits of measurement accuracybest relative errorTikhonovtruncationCGLS stepslope 0.89, optimalslope 2/3, Tikhonov's limithow fast the error fallsTikhonov slope0.7truncation slope0.94CGLS step slope0.89optimal 2ν/(2ν + 1)0.89Tikhonov's limit0.67better data, smaller errorat a rate the method may not be able to keep Methods that were designed apart

The method that cannot use a smooth answer

On the collection's own signal four regularisers reach the same floor to a few per cent. Score them instead against answers of increasing smoothness and one stops improving. Across six decades of noise Tikhonov's error falls with a fitted slope of 0.70 whether the answer is twice or four times as smooth, while truncation's rises to 0.94 — and at 10⁻⁶ noise Tikhonov's best is 38 times truncation's.

 -1-0.500.511.520eigenvalue of P⁻¹K1 − φ1φS never formedγ1smallest ν0.056furthest from φ, 1 − φ0.57MINRES steps11an exact Schur approximation, bought by changing Hthe golden ratio without S The matrix a constraint makes

Where the augmentation puts the cost

Add γAᵀA to the objective block of a saddle-point system and its Schur complement tends to I/γ, so the cheapest possible approximation becomes the right one and the golden-ratio spectrum arrives — within 7.6·10⁻⁶ at γ = 10⁶. MINRES falls from 21 steps to 6. The inner solve with the augmented block rises from 14 conjugate gradient steps to 43, their product does not fall at all, and the answer loses seven and a half digits on the way.

 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.

 0510152025303540012345matrix size ngrowth factorpartial pivotingrook pivotingcomplete pivotingsolid: median of 30 · dashed: worst of themat n = 40partial pivoting, median3.3rook pivoting, median2complete pivoting, median1.6rook, worst of the draw2.8the same matrices under every rulerook 3.3× partial's search Elimination, and the swap

A pivot that searches one row and one column

Rook pivoting looks down a column for its largest entry, along that entry's row for a larger one, and back down that entry's column, until it finds an entry largest in both. On Gaussian matrices of size 64 it keeps the median growth factor at 2.53 against partial pivoting's 4.06 and complete pivoting's 1.88, and it compares 6,987 entries against 2,080 and 89,440. On Wilkinson's matrix it holds the growth at exactly 2 where partial pivoting reaches 9.2·10¹⁸. And on a matrix built to make it walk it compares 113,376 entries — more than complete pivoting.

 08162432404810⁻¹10²10⁵10⁸10¹¹10¹⁴10¹⁷10²⁰degreeκ of the design matrix1/u: past this a double holds nothingmonomialsChebyshevArnoldi on the pointsκ at degree 48monomials1.5·10¹⁷Chebyshev1.6·10⁷Arnoldi on the points1200 two intervals with a gap between themthree bases for one space of polynomials Least squares, and the road not to take

A basis built from the points

A polynomial fit computed in monomials and in an orthogonal basis gives the same curve on exact data, and the valley essay drew the two lying on top of each other. Add 0.1% noise and they separate — by 1.7·10⁻⁵ at degree 40 and 0.004 at degree 48 — because the fitted curve moves with the basis by its condition number times the rounding times the residual. Chebyshev polynomials keep that small only on points spread like their weight; on a sample with a hole in it they reach κ = 1.55·10⁷. A basis orthogonalised against the sample points themselves stays at 1 on every set.

 10²10³10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹steps takenlargest relative error in a coefficientcarried factorrecomputedcarried + one correction+ a second correctionHouseholder QR, same rowsdrift of the factorafter 2875 stepsκ(A) of the window2·10⁴carried factor2.8·10⁻⁸recomputed5.8·10⁻⁹carried + one correction7.1·10⁻¹³+ a second correction3.7·10⁻¹²Householder QR, same rows1.1·10⁻¹²measured against the coefficients in exact rationalsthe refresh repairs the factor, not the answer When the problem arrives again

The repair the drift did not need

A sliding window's carried Cholesky factor drifts 3.9·10⁻¹⁴ from its data, and multiplying by κ(AᵀA) predicts eight lost digits in the coefficients, a stream conditioned at 10¹² losing the answer, and a periodic refresh of the factor as the default repair. Measured against coefficients computed exactly in rationals, all three come out differently. On a stream made ill-conditioned by scaling, the conditioning never reaches the coefficients. On a collinear stream, a freshly recomputed factor is as wrong as the drifted one. And one correction from the window's own rows reaches Householder's accuracy for a fraction of a refresh's cost.

 significand bits10³10⁴10⁵10⁶10⁷10⁸10⁹10¹⁰κ(A)87892634456329625911109435761199127514162024every matrix positive definiteruns producing a false certificate14of runs in total72never above, in significand bits12first κ at eight bits10⁵the comparison was correctand what it proved was not true Two errors, and whose fault they are

Deciding that a zero has arrived

The previous tolerances were offers — accept this much error, save this much work. A detection threshold is not an offer, because both directions are failures. One matrix here has three genuinely near-invariant subspaces, and the constant somebody typed decides which of them the recurrence stops at; at eight significand bits the same kind of constant produces a proof of something false.

 263850627410⁻⁵10⁻³10⁻¹rows in the linearisation‖g − r‖ on the target setpolynomialrationalequal cost, two basesreach0.06largest size drawn78rational there2·10⁻⁵polynomial there0.0015the ratio73no winner on an easy targetand two orders on a hard one The eigenvalue problem that is not linear

Two approximants and one matrix size

A polynomial approximant linearises to nd rows and a rational one to n(m+1), so the fair contest fixes the matrix and varies the basis. On an easy target set the two are indistinguishable and the ordering flips with the noise; on one that reaches a branch point the rational pulls away by two orders.

 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.

 02505007501000250275300325350375additionsrunning totalround to nearest: nothing arrivesexactstochasticnearesta thousand additionshalf an ulp at 2561moves, round to nearest0moves, stochastic46relative error, nearest0.28relative error, stochastic0.0228 significand bits, unbounded exponenta flat line is not a small error The arithmetic underneath

A coin flip that fixes the average

Add 0.1 to 256 a thousand times at eight significand bits and the answer is 256. Not approximately — the total never moves, not once, and no error bound says so. Round up one time in twenty instead of never, and it arrives at 348 against a true 356.

 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.

 -10123456789-101λbranch point at −0.4the closed formwhich of these is an answereigenvalues returned36wanted6past the branch point6complex24worst against the closed form7.7·10⁻⁴all of them exactfor a problem nobody asked The eigenvalue problem that is not linear

The eigenvalues that are answers to nothing

A rational approximant of degree five turns a six-by-six problem into a thirty-six-by-thirty-six one, and thirty-six numbers come back. Six are the answer. The rest are exact eigenvalues of the approximant, lying where the function it approximates is not a real number at all.

 places settled, of 14the arithmetic14 of 14one missing link3 of 14α from 0.80 to 0.9011 of 14entries move by 4.05·10⁻¹¹entries move by 0.00443entries move by 0.00123three kinds of doubtcomparisons14smallest gap6.3·10⁻⁵largest gap0.0024arithmetic band4.1·10⁻¹¹one-link band0.0044α band0.0012the arithmetic settles nearly all of themthe data settles three The matrix that is a graph

A ranking whose order is not determined

Of the fourteen adjacent comparisons in a top-fifteen, fourteen survive perturbing the arithmetic at the rounding level, eleven survive moving the teleportation parameter across its usual range, and three survive removing one link. The computation is the strongest part of the answer.

 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.

 02040608010010⁻¹¹10⁻¹⁰10⁻⁹10⁻⁸truncations performed by the factorisationrelative errorthe representation's own error‖A − LLᵀ‖ ⁄ ‖A‖no decomposition without its residualresidual, 0 truncations2.1·10⁻¹⁰residual, 98 truncations1.1·10⁻⁹representation, deepest1.4·10⁻⁹residual ⁄ representation0.81levels, deepest5a hundred approximate stepsand an exact-looking factorisation Neither sparse nor dense

The count that is not the budget

A Cholesky performed inside a low-rank format truncates 0, 2, 10, 34 and 98 times as the leaf falls from 128 to 8, and those five integers are the same at every accuracy from 10⁻¹² to 10⁻². Across all ten decades the factorisation's residual stays below the representation's own error at a ratio between 0.81 and 1.00 — with two entries that read 1.83 and 1.78, and neither of them is accumulation.

 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.

All themes