Numerical Linear Algebra
where the algebra stops being the arithmetic
Theme

The thread: Measured, not assumed — page 15

Essays 337 to 360 of the 564 on this theme, in the same order.
10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²-1-0.500.51relative backward error acceptedfraction of the work not doneQR deflation criterionincomplete Cholesky droplow-rank truncationthree fields, no shared arithmeticslope, deflation criterion0.039slope, drop tolerance0.051slope, rank truncation0.23widest apart, as a ratio5.7a tolerance is an offerand the three offers are one offer Two errors, and whose fault they are

A tolerance is priced by the problem

Three tolerances from three fields sit on one pair of axes and agree to within a factor of 5.74. That factor is the ratio of the two curves that cannot move. Change the only problem in the comparison and the third curve's fitted slope swings from 0.188 to 0.040 while the printed spread does not shift by a digit.

 027548110813510⁻⁶10⁻⁴10⁻²110²10⁴10⁶power‖Aᵏ‖Kreiss constant 6760e · n · K‖Aᵏ‖ρᵏtwo routes to one peakspectral radius0.8peak of ‖Aᵏ‖2·10⁴Kreiss constant6757e · n · K1.1·10⁵everything here decays in the endand one of these curves says how much first Eigenvalues, singular values, rank

A spectral radius that grows first

ρ(A) below one guarantees that the powers of A go to zero and says nothing about what they do on the way. Here they rise by a factor of twenty thousand before turning over, and the peak is bracketed above and below by a constant computed from the resolvent norms outside the unit circle — two routes to one number, one through the plane and one through the powers.

 11.31.610⁻¹110¹10²10³log₁₀ ntropical root ÷ actual modulusexactsmallest moduluslargest modulusthree norms, two endslarge root, n = 40.96large root, n = 320.94small root, n = 45.4small root, n = 32230a maximum predicts a maximumand says nothing about a minimum The eigenvalue problem that is not linear

An estimate that does not move

The tropical roots are said to miss the bottom of a spectrum by a factor growing like n². The estimate does not get worse. It changes by four per cent between four masses and sixty-four while the modulus it names falls by a factor of a hundred and sixty, and an estimate that is constant in the variable the answer depends on is a different defect from an inaccurate one.

 110¹10²10³110¹10²10³the accumulating quantity, relative to its first valuethe error, relative to its first valuethe bounds: slope 1what all three do: slope ½three mechanisms, one exponenta left-to-right sum0.49a chain of rotations0.55a residual recurrence0.51every bound's slope1spread of the three0.067a bound is a sum of the roundingsand the roundings have signs The arithmetic underneath

Three walks and one bound

A left-to-right sum, a chain of three thousand rotations and a conjugate gradient residual recurrence share no arithmetic and no vocabulary. Each has a standard bound that is linear in whatever it accumulates against. All three come out at a half — 0.486, 0.554 and 0.507 — and nothing is rescaled.

 n = 512, leaf 16blocks, tightest to loosest592and at the loosest94mean rank, tightest3.3and at the loosest8.9052104156208260the separation constantstored per unknown0.20.30.40.50.951.4stored per unknownmean rank, compressed blocksboth curves move togetherwhich is what the other two knobs never do Neither sparse nor dense

A geometry setting that is a second accuracy

The tolerance raises every rank and moves no block; the leaf moves every block and no rank. The constant that decides how far apart two clusters must be before their block may be compressed moves both — 592 blocks at mean rank 3.25 at one end and 94 at mean rank 8.94 at the other — and the storage falls the whole way while the error against the dense matrix rises by a factor of 3.5, at a tolerance that never changed. It has no best setting, and it is not a third knob on the storage.

 a 28 × 28 saddle-point matrixnatural order, entries113sparsest order, entries72two-by-two pivots, loosest test0and at the strictest3023466992115the constant in the pivot testentries in the factor0.10.30.640.80.95133natural ordersparsest ordergrowth, naturalgrowth, sparsestnumbers above the points: two-by-two pivots takenboth curves favour the sparser order Sparsity, and what elimination costs

The freedom a symmetric factorisation does not have

Permuting rows and columns together leaves no column to choose, so the conflict between the sparsest pivot and the sound one should be worse rather than better. On a saddle-point matrix whose constraint rows have no diagonal entry at all, it is not there: taking the sparsest available pivot holds 70 entries against the natural order's 113 and a growth of 1.28 against 1.83 — better on both currencies at once, at every setting of the pivot test. The two-by-two blocks that make it legal cost 1.33 entries apiece.

 each method's best, dimensionconjugate gradients24Tikhonov23truncated SVD21and the error thereconjugate gradients0.14Tikhonov0.14truncated SVD0.140510152025301effective dimensionrelative errorconjugate gradientsTikhonovtruncated SVDdashed vertical: the iteration's bestthe iteration drawn while its dimension still rises Methods that were designed apart

One arc, and what each filter pays to be on it

Conjugate gradients and Tikhonov stop at the same effective dimension, and that could have meant two curves crossing once or one curve. It is one curve over a stretch — on six problems, Tikhonov and truncation reach the iteration's error at the iteration's dimension to within 7.2 per cent from 0.7 of the answer to its top — and the two separate on either side. But the curve is shared by a trade, not by an identical answer: at the same dimension Tikhonov carries 22 to 42 per cent more noise than the iteration and up to five per cent less bias.

 worst draw on any grid, % over the oraclediscrepancy13GCV4·10⁶L-curve291median on 96 points, % overdiscrepancy3.4GCV0.18L-curve282030405060708090100110¹grid pointserror ÷ the oracle's, same drawdiscrepancy, mediandiscrepancy, worstGCV, medianGCV, worstL-curve, mediansolid: the median draw; dashed: the worst of sixteenthe axis stops at twenty times the oracle Regularisation, and the answer that is chosen

The data count their dimensions, not the step's

Every grid in the deconvolution essays was chosen with the answer in hand, and so was every λ. From the data alone, the discrepancy principle's worst draw is within 16 per cent of the oracle on every grid from 16 points to 96; generalised cross-validation is better on the median draw and, on grids of thirty points and more, has draws thousands of times worse. And the data can say how many dimensions they carry — about 20, 25 and 29 at three noise levels, one number once the grid exceeds it — but not how many more the step needs: the grid that count chooses is 14 to 19 per cent worse than forty points at the lower two.

 6 numbers, coefficients to 30relations, both agree34relations, they differ1accidents, both agree5accidents, they differ88510152025300123decimal digits the numbers are scaled togap, in digitsrelation, both agreerelation, they differaccident, both agreeaccident, they differhorizontal: a gap of one digitdashed vertical: the digits a double has Exact arithmetic, and what it costs instead

Two precisions guard the other edge

Run a lattice relation search at N digits and again at N/10, and accept its answer only if both runs return the same vector. Among six measured numbers, where the gap between the shortest and next vector vouches for one found relation in 35, the two runs agree on 34. Past a double's sixteen digits they never once agree on an accident, 0 of 561. They do agree on 62 accidents at fifteen digits or fewer — approximate relations that really are the shortest vector there — and the gap, which cannot see a relation among six numbers, can see those.

 no noisesweeps to stop49final error2.5·10⁻¹⁴1% noisesweeps to stop6000largest term, ÷ data0.67110¹10²10³10⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1sweeprelative change in the error, per sweepthe stopping test1% noiseno noise: stops at 10⁻¹⁴the clean run stops; the noisy one crawlsits terms stay the size they were When the index is a tuple

A fit that has an answer and cannot stop

Fit a noisy rank-three tensor with four terms and the prediction was that the spare term would find real structure in the noise and stop pairing off with the others. It does not: the largest cosine between two fitted terms has a median of 0.977 at 1% noise against 0.975 without noise. What changes is the solver. Without noise the four-term fit stops in fifty sweeps; with noise it never stops: its error keeps moving by a part in a million a sweep for six thousand sweeps, while on most tensors its terms stand still. The extra term takes exactly its share of the noise, and the error is settled by sweep 150.

 decay 0.8, 400 draws3% target, c = 10.663% target, c = 1.960.9410% target, c = 10.6611.251.51.7522.252.5556065707580859095100margin c on the standard errordraws inside the target, %3% target10% targeta normal tablethe dashed curve is 2Φ(c) − 1what a margin of c promises if the error is normal Randomised, and the guarantee that changes kind

The miss a normal table already priced

A trace estimator that stops when its own standard error reaches a target misses the target on about a third of draws, and the essay that measured it read the loose targets as the worst calibrated. Over 400 draws the loose target is the better covered — 76% at 10% against 64% at 3% — because the warm-up stops most of its runs with probes to spare. Where the criterion decides, the misses are a normal distribution's: a margin of c on the standard error buys what a normal table says, 94.8% at 1.96, and costs c² in probes, 3.86 times.

 heavy row first1 − h at 4³⁰7.4·10⁻¹⁸κ(A) at 4³⁰3.1·10⁸complement's error2.1·10⁻¹⁶110³10⁶10⁹10¹²10¹⁵10¹⁸10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹weight on the rowrelative error of 1 − h1 − ‖R⁻ᵀa‖²1 − ‖row of Q₁‖²‖row of Q₂‖²κ(A)·uthe rows are the same in every positiononly the order the factorisation meets them changes Least squares, and the road not to take

The weight the factor met first

The route to one minus a leverage through the orthogonal factor was said to lose a digit for every decade of the condition number, whatever else it does. Put a weight on one row and it does not. With the heavy row first, the complement keeps every digit at κ(A) = 2.5·10⁹ while both subtractions return nothing. With the same row last it loses digits as the row's scale grows. And two heavy rows that leave κ(A) at 3.1 still lose six digits when the light rows come first. The law was about the order the factor met the rows, and the condition number had been standing in for it.

 noise none, mediansthe answer's step0fresh start, corrected6.2·10⁻⁹carried start, corrected1.3·10⁻¹¹0200400600800100010⁻¹²10⁻¹⁰10⁻⁸steps along the streamlargest relative coefficient errorfresh solve, one correctioncarried answer, one correctionfresh solve, two correctionsHouseholder on the rowseach correction contracts its start by the same factorthe nearer start wins When the problem arrives again

The answer the last window left

A sliding window that corrects its least-squares answer at every step could start each correction from the previous step's corrected answer instead of from a fresh solve: the two windows share all but one row. On a stream with any noise in it, that start is three orders worse. The window's exact answer moves by 0.79 of itself in one step at κ(A) = 3·10⁶, a fresh seminormal solve is wrong by only 1.8·10⁻⁴, and one correction contracts either start by the same factor — so the fresh start ends at 2.2·10⁻⁸ and the carried one at 3.7·10⁻⁵. Only on data that agree exactly does carrying win.

 Péclet 100, orders 1, 2, 3, 6nodal edge at 1,000, outflow8.9weighted edge at 1,000-0.097unstable nodal models7110¹10²10³-20246810largest cell over smallestright edge of the numerical rangenodal, small cells at the outflownodal, small cells at the inflowcell-size weightedabove the dashed line the certificate is goneit is a sufficient condition, not a failure Reduction, and what a model is for

The inner product the mesh already computed

A Galerkin reduced model is certified stable when the operator's numerical range sits in the left half plane, and the certificate belongs to the coordinates. The coordinates that broke it before were a constructed rescaling. Grade a convection–diffusion mesh towards its outflow boundary layer — the grading anyone resolving the layer would choose — and nodal values break it too: the range's edge is past zero at a ratio of ten and reaches +8.9 at a thousand, and at Péclet 100 seven reduced models come back unstable. Weight the projection by the cell sizes the discretisation already computed and every model at every grading is stable, with the range's edge back at −0.098.

 three seeds, medianblock MGS, Qᵀb, κ 10⁸0.23block MGS, b as a block1.1·10⁻⁹Householder2.8·10⁻⁹10²10⁴10⁶10⁸10¹⁰10¹²10⁻¹⁶10⁻¹²10⁻⁸10⁻⁴110⁴10⁸κ(A)relative error of xblock MGS, Qᵀbblock CGS, either routeblock MGS, b as a blockHouseholderthe same Q in both block MGS routesonly the order in which b meets it differs Orthogonality, measured

What the appended block inherits

Modified Gram–Schmidt on [A b] solves least squares as well as Householder, although its Q is not orthogonal. A block code appends b as one more block. Block modified Gram–Schmidt inherits the rescue at every placement of the ill-conditioning: at κ = 10⁸ the appended block gives 6.9·10⁻¹⁰ where the same Q through Qᵀb gives 8.9·10⁻³. Block classical Gram–Schmidt gets the same wrong answer both ways, to the last bit. And the variant whose Q is orthogonal to 10⁻¹⁵ — two passes with Cholesky QR inside — is a hundred thousand times worse than Householder when the ill-conditioning is inside the blocks, because its R is wrong.

 noise that halves the growthshooting, partial pivoting0.0018shooting, τ = 0.10.56Wilkinson, τ = 0.10.1800.10.20.30.40.50.60.70.80.9110⁻¹⁶10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1one minus the thresholdstored noise that halves the growthshooting, T = 12, h = 0.3Wilkinson, n = 24the same growth at every thresholdheld by a choice noise must push further to reverse Elimination, and the swap

A threshold that holds the growth still

Partial pivoting's large growth on the shooting matrix rests on a margin of 5.6·10⁻³, and Wilkinson's rests on exact ties, which a perturbation at rounding breaks. A sparse code pivots with a threshold instead, taking the sparsest row among candidates within a factor τ of the largest. At every threshold from 1 to 0.1 both matrices grow by exactly as much as under partial pivoting. But the growth is now held by row counts, which no perturbation of the values can reverse: at τ = 0.1 it survives noise on the stored entries up to 0.56 on the shooting matrix, three hundred times more, and 0.18 on Wilkinson's, where partial pivoting's goes at 10⁻¹⁶. The margin that predicts it is the chosen candidate's height above the threshold line, and it needs no division by the growth.

 -12-10-8-6-4-2024681011.522.533.5leaf width minus (√M − 2)words ÷ pure recursionM = 64M = 144M = 256the edge is the same column at every memoryand the cheapest leaf sits on it Where the flop count stopped predicting the time

The leaf that sits on the edge

A recursive elimination's base case was found to be a block size in disguise, with a cliff where the blocked elimination's is, and the choice read as a trade: the processor wants wide leaves, the cache wants narrow ones. Measured at every width rather than at the halvings of 96, there is no trade inside the edge. Leaves exactly √M − 2 wide are the cheapest the recursion can have in words as well as calls — 0.81, 0.80 and 0.84 of the pure recursion's traffic at 64, 144 and 256 words — and one column wider moves 1.77 to 3.15 times it. And a matrix of 100 columns, which halves unevenly, meets the cliff in two steps rather than one.

 n = 64, twist of half a stepnearest sample, Laplacian0.0024nearest sample, squared5.8·10⁻⁶κ(wrap), squared2.8·10⁶10⁻³10⁻²10⁻¹110⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹θsymbolπ/nLaplacian, θ² near 0its square, θ⁴ near 0no twist puts a sample further than π/n from θ = 0and the symbol's order decides what it sees there Structure, and the solver that cannot see it

A zero no twist can step around

Solved through its circulant wrap with the small capacitance system factorised stably, a banded matrix was found to lose only what the cancellation in the correction costs, and a twist of half a step kept that cancellation near one. That holds only while the cancellation is large. With it kept small, the corrected solve's error is a quarter of the wrap's condition number times u on every case measured — and on a symbol with a zero of fourth order, the square of the Laplacian, no twist can keep the wrap well conditioned, because the nearest sample any twist can reach sees (π/n)⁴. At n = 64 the corrected solve is 34 times less accurate than elimination, where the Laplacian's is four.

 κ(A) = 10⁸, 152 saddlespassed, ordinary63passed, on wrong counts57passed, always510⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹110¹10⁻⁴10⁻²110²10⁴10⁶10⁸curvature needed, μshift applied ÷ μordinary schedulenull space on wrong countsnull space alwaysbelow the dashed line a saddle was certifiedthe arbitrated shift sits at twice the need The matrix a constraint makes

A loop that asks the null space why

An inertia-correction loop sees only an integer, and two different faults produce the same wrong one: curvature that needs a shift, and a constraint too weak for the count to see. One QR of the constraint matrix on a wrong count tells them apart — it shifts none of the 22 weak-constraint minima the ordinary loop shifted by up to 2,621 — and its reduced eigenvalue gives the shift a saddle needs in one step, twice the need exactly, where the schedule overshoots by up to 17,783 times. But at κ(A) = 10⁸ the loop still certifies 57 saddles of 152, because a false certificate is a count that read right, and a check made only on wrong counts never sees it. Asking every time leaves five, all shallower than 2·10⁻⁸.

 75 aggregates over 225 unknownsthe matrix chose thismean extent along y3mean extent along x1points adopted by pass two15no coordinate enters the methodand the shape follows the coupling Iterating, instead of factorising

Aggregating what the matrix calls strong

The depth phase measured every method it had on the 45°-rotated anisotropic operator — 0.784, 0.883, 0.844 — and diagnosed the failure as being in the discretisation rather than in the hierarchy. Smoothed aggregation is the standard answer to anisotropy. It returns 0.789.

 κ of the suma component of b − Ax8.09·10¹⁷ad − bc, near-degenerate3.6·10¹⁶qᵢᵀqⱼ, an orthogonality check3.26·10¹⁶zᵀAz, a trace probe95.7pᵀAp, a curvature73.7rᵀr, a residual norm1measured, not assumedhighest8.1·10¹⁷lowest1above 10¹⁰3terms128sums of squares are safeand nobody decides anything from one The answer that depends on the machine

Two machines, one certificate

Nothing a solver returns says which of its answers you got. Four things could be reported instead — the summation condition number, the partition count, an exactly accumulated residual and a directed-rounding interval — and each costs about one pass over data the routine already has in hand.

 finite eigenvalues (degree of det Q)13at infinity (2n − degree)1at infinity, by the rank of M12n, if M were nonsingular14a degree, not a decisiondegree of det Q13at infinity1by the rank of M1singular-value gap∞the count is a degreeand the other route is a judgement The eigenvalue problem that is not linear

One mass removed, and one eigenvalue gone

A coordinate with no inertia reads like a coordinate that has been deleted, and a chain of eight masses with one of them removed would then be a chain of seven, with fourteen eigenvalues. It has fifteen. The massless coordinate is still there, still carrying a damper, and it contributes a first-order equation rather than none.

 each bar is a percentage — of the sample, or of the true condition numberexactly right86.0%inside 10%91.0%inside a factor of 291.0%worst in the sample, ×10035.7%the constructed matrix, ×1007.7%usually exactexact share0.86worst of the sample0.36the constructed matrix0.077a routine that is right most of the timeand never wrong in the safe direction Two errors, and whose fault they are

The tail a sample never reaches

Hager's estimator is exactly right on four random matrices in five, and that share is stable — between 80.5 and 87.5 per cent across nine sizes. The worst underestimate is not stable at all: it falls every time more matrices are drawn, from 0.746 at sixty to 0.377 at four hundred, and the matrix built to defeat the estimator sits five times below anything four hundred draws found.

 02468101201change of units, by exponent10⁰10¹10²10³10⁴10⁶10¹⁰10¹⁶10¹⁹10²⁰10⁴⁰10¹⁵⁰as writtenafter scalinga range questionlargest finite value3.4·10³⁸predicted boundary γ1.8·10¹⁹last γ that forms10¹⁹stops where scaling fails0not a poor answerno answer at all The arithmetic underneath

The units that overflow before the answer does

A change of variable that is exact in the algebra requires γ² times a matrix to be a number the format can hold. In binary64 that is a bound nobody meets by accident. In binary32 it arrives at 10¹⁹ and in fp16 at 256, and past it there is no answer rather than a poor one.

All themes