Numerical Linear Algebra
where the algebra stops being the arithmetic
Theme

The thread: Where the algebra stops being the arithmetic — page 6

Essays 121 to 144 of the 259 on this theme, in the same order.
the invariant planesolid: beforedashed: aftersame perturbation, two questionsthe vectors turned, radians0.029the plane turned, radians7.6·10⁻⁸what left the plane5.6·10⁻⁸drawn in the unperturbed plane's own basisa radius is not determined; the circle is Eigenvalues, singular values, rank

The plane survives what its vectors do not

At a gap of 10⁻⁹ a perturbation of 10⁻⁶ turns the two eigenvectors through half a radian and turns the plane they span through 7.6·10⁻⁸ — a ratio of six million. Ask for the subspace instead of the vectors and a hopeless computation becomes a well-conditioned one, with no change to the arithmetic.

 -101-101(x − 1) ÷ 0.001p(x) ÷ 2.68·10⁻¹⁵a square, drawn negativepoints sampled401negative, unfused179negative, fused196signs disagree at98widest gap1.4·10⁻¹⁵unfusedfusedthe true valuethe true curve is the flat oneeverything else is rounding The answer that depends on the machine

A square that evaluates negative

(x − 1)⁶ evaluated near x = 1 comes out negative at 179 of 401 points on one build and 196 on another, and the two disagree about the sign at 98 of them. Neither is nearer the truth: both traces are made entirely of rounding.

 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.

 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.

 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.

 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.

 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.

 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.

 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.

 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.

 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.

 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.

All themes