Numerical Linear Algebra
where the algebra stops being the arithmetic
Theme

The thread: Measured, not asserted — page 4

Essays 73 to 96 of the 110 on this theme, in the same order.
0.46-1.11.5-0.470.2-0.181.90.690.42-0.160.028-0.0160-1.50.160.880.099-0.1100.680.762.30.930.0480-0.250.110.952.11.500001.52.3after reflector 2the highlighted entries are the bulge — the only thing that is not Hessenbergentries below the subdiagonal3reflectors used2reflectors in a whole step5the shifts are never formed — only their sum and productand both of those are real Eigenvalues, singular values, rank

Two shifts that are never formed

The double shift is defined as a factorisation of (A − μI)(A − μ̄I), which nobody computes. What is computed is the first column of that product — three numbers — and the bulge those three numbers create, pushed down the subdiagonal by n − 2 reflectors until it falls off the bottom.

 -47-37-27-17-7313233301234log₁₀ of the vector's normfp1611 bitsbfloat168 bitstf3211 bitsbinary3224 bitspale: the format's range · blue: √(Σ(xᵢ/m)²)·m · red: √(Σxᵢ²)fp16 and tf32 have the same eleven significand bitsand their bars do not overlap The arithmetic underneath

A norm that overflows before it is a norm

The vector of sixteen thousands has a Euclidean norm of 4,000, which fp16 represents exactly. Written as the square root of the sum of squares it returns infinity, because squaring doubles the exponent — and the expression costs half the format's range on the one computation every iterative method performs at every step.

 the estimator maximises this quantity over the columns it visitscolumn 1 ‹visited›12column 2 ‹the answer›114column 311.4column 411.4column 511.4column 611.4column 711.4column 811.4column 911.4column 1011.4column 1111.4column 1211.4estimate 12.0a walk that stopped earlythe estimate returned12the true 1-norm114columns visited1products with the matrix5the walk's own stopping test firedand every column it could see was smaller Two errors, and whose fault they are

An estimate that can be fooled

Nobody computes a condition number, because forming an inverse costs more than the solve did. Every library estimates it instead, from four or five products with a factorisation already in hand. The estimate is exactly right on four random matrices out of five — and there is a matrix, three distinct entries wide, on which it returns a twentieth of the truth.

 10⁻²10⁻¹110¹10²10³10⁴10⁻⁴10⁻³10⁻²10⁻¹1backward error, in units of ushare of systems above ituCramereliminationCramer, control24-bit arithmeticworst Cramer, in u395worst elimination, in u1.2control, worst Cramer4.9κ of the worst system3.4·10⁶one derivation, two computationsand only one of them is stable Elimination, and the swap

A rule that is correct and unusable

Cramer's rule gives every component of the solution in closed form, in terms of determinants, and it is a theorem. On two-by-two systems whose rows are nearly parallel it returns an answer with a backward error of 458 units of roundoff where elimination returns 1.3 — on a matrix whose condition number is 32,000 and which elimination solved perfectly.

 10¹10¹.³10¹.⁶10¹.⁹⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰⁰¹10².²00.250.50.751grid points nresidual reduction per stepJacobiGauss–SeidelV-cycleV-cycle spread, 8× in size0.003Jacobi at n = 1270.99work exponent, fitted0.079the dashed curve is cos(πh), Jacobi's closed formthe flat line is the whole method Iterating, instead of factorising

A rate that does not notice the size

The V-cycle reduces the residual by a factor of ten a cycle at fifteen points and at a hundred and twenty-seven. Jacobi on the same four problems goes from 0.981 to 0.9978, climbing towards one. One of those is a constant and the other is an exponent, and that is the whole distinction the field turns on.

 110¹10²10³110¹10²10³off-diagonal entry ccondition number of the eigenvalue√(1 + c²)decoupled: 1measuredthree routes, one number‖A − ZTZᵀ‖/‖A‖1.7·10⁻¹⁵closed form100computed 1/|yᵀx|100worst measured movement46four eigenvalues, two conditioning numbersthe symmetric case has one, and it is 1 Eigenvalues, singular values, rank

A condition number for one eigenvalue

In the symmetric case every eigenvalue has condition number exactly one. In this four-by-four matrix two of them have condition number 100.005 and the other two have exactly 1, and the number belongs to the eigenvalue rather than to the matrix.

 10⁻⁷10⁻⁵10⁻³10⁻¹⁰10⁻⁸10⁻⁶magnitudespacing to the next numberthe smallest normalgradualflush to zerosmallest normal6.1·10⁻⁵smallest subnormal6·10⁻⁸octaves of subnormals10pairs that lie under FTZ10the spacing stops halving and stays putwhich is what makes x − y = 0 mean x = y The arithmetic underneath

The numbers below the smallest one

Below the smallest normal number the spacing stops halving and stays put, all the way to zero. That is what gradual underflow is, and the thing it buys is the sentence every algorithm assumes without being told — x minus y is zero only when x equals y.

 does this matrix look nearly singular?green: the test agrees with the truth · red: it does not · the bar under each number is its magnitude, over sixty-two decades|det A||det A|^(1/n)σ_min1/κ = σ_min/σ_max0.1·I at n = 40perfectly conditioned10⁻⁴⁰0.10.11κ = 10¹⁰, |det| = 1nearly singular1110·10⁻⁶10·10⁻¹¹Hilbert at n = 8nearly singular2.7·10⁻³³8.5·10⁻⁵1.1·10⁻¹⁰6.6·10⁻¹¹the two counterexamplesκ of the scaled identity1its determinant10⁻⁴⁰κ of the normalised matrix10¹⁰its determinant1det(cA) = cⁿ det(A)so a determinant carries the units n times over Two errors, and whose fault they are

The number that decides nothing

The determinant is the first scalar anybody attaches to a matrix and the last one worth consulting. A tenth of the identity has a determinant of 10⁻⁶⁰ and a condition number of exactly one. The Hilbert matrix's determinant stops being right at n = 13 and stops being a number at n = 29, and nothing in between reports either.

 10²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1κ₂(A)relative errorforward, A⁻¹bforward, A\bbackward, A⁻¹bbackward, A\bat κ = 10¹⁴η, LU solve2.2·10⁻¹⁷η, via the inverse4.5·10⁻⁵forward, LU solve2.8·10⁻⁴forward, via the inverse0.015one factorisation, two ways to use itand one of them forfeits the backward error Elimination, and the swap

The inverse that is never formed

x = A⁻¹b is how the solution of a linear system is written and it is not how it is computed. The usual reason given is cost — three times the arithmetic. The real reason is that one of the two routes is backward stable and the other is not, and at κ = 10¹⁴ they differ by twelve orders of magnitude in the number that says whose fault a wrong answer is.

 -0.0833-0.167-0.0833-0.1671-0.167-0.0833-0.167-0.0833R · A · P, normalised to a unit centrestored entries per rowlevel 0 · 31×314.87/rowlevel 1 · 15×158.22/rowlevel 2 · 7×77.37/rowlevel 3 · 3×35.44/rowlevel 4 · 1×11.00/rowstill a stencil, still annihilates a constantentries in an interior row9weight outside the 3×30row sum0the isotropic model problemnine, at every level below the first Iterating, instead of factorising

The coarse problem is a different problem

In one dimension the Galerkin coarse operator is the coarse discretisation, entry for entry — this site asserted it. In two dimensions a five-point operator produces a nine-point coarse one, so the recursion solves a different discretisation at every level below the first, and converges at 0.20 a cycle regardless.

 10⁻³10⁻¹10¹10³01magnitude448NaN — no ∞0.0156 — the smallest normaldrawn from the format's own rulespositive finite values126largest finite value448worst round-trip error0the subnormals are the evenly spaced ticks at the lefteverything a byte can be The arithmetic underneath

Eight bits, and a format that breaks the rules

E4M3 reuses the exponent code IEEE reserves for infinities, so it reaches 448 where the same bits under IEEE's rules would reach 240 — and has no infinity left to signal an overflow with. The same computation is a NaN on one conforming device and 448 on another.

 10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻⁸10⁻⁵10⁻²gap between the two eigenvalueshow far it movedthe eigenvectorsthe eigenvaluestheir plane‖E‖ / gapone perturbation, three answerseigenvalue shift, spread over the sweep1plane angle, spread over the sweep1eigenvector angle, spread1.6·10⁵the dashed line is Davis–Kahan's ‖E‖/gaptwo of the three never noticed Eigenvalues, singular values, rank

The gap decides the eigenvector

A symmetric matrix's eigenvalues move by at most the size of the perturbation, whatever the spectrum looks like. Its eigenvectors are governed by a completely different quantity — the distance to the neighbouring eigenvalue — and at a gap of 10⁻⁹ the same perturbation turns them through 27°.

 θx (frequency across x)θy0π/2π0π/2πthe coarse grid'sunder 0.2under 0.40under 0.60under 0.80under 0.95under 1.01damping per sweeptwo routessmoothing factor, scanned1closed form130×30 frequency cellsthe marker is the mode nothing removes Iterating, instead of factorising

A direction the smoother cannot see

Give the Laplacian a strong direction and multigrid stops working — from 0.2016 a cycle to 0.9565 — with every component unchanged and the condition number identical to twelve digits. The problem did not get harder. The link between the method's two halves broke.

 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.

 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.

 λ = 105 timesλ = 9.55 timesλ = 95 timesλ = 8.55 timesλ = 2.952 timesλ = 2.92 timeseigenvalues that arrived more than once — the matrix has 40 distinct onesa spectrum with the wrong multiplicitiesextra copies, no reorthogonalisation25extra copies, full reorthogonalisation0worst relative error among the copies1.9·10⁻⁸steps taken of 80 asked for, full40no arithmetic error was madeevery one of these is right to eight digits Eigenvalues, singular values, rank

An eigenvalue that arrives twice

A matrix with forty distinct eigenvalues, handed to Lanczos for eighty steps, returns twenty-five extra copies of thirteen of them — the largest arriving five times. Every copy is accurate to 1.9·10⁻⁸ relative. No arithmetic error was made, nothing overflowed, and a caller counting eigenvalues gets the wrong multiplicity from a computation in which no individual number is wrong.

 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.

 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.

 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.

 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.

 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.

 the grid the operator came froman edge is a coupling the matrix calls strongkeptinterpolatedwhat the entries decidedcoupling ratio, x against y0.001strong couplings across x0strong couplings along y110rows kept or dropped whole1111×11 grid, θ = 0.25the coarse grid, from the matrix alone Iterating, instead of factorising

The coarse grid the matrix chooses

Given a tridiagonal matrix and no information about a grid, the coarsening keeps every other point and derives the weights ½, 1, ½ — the operators the geometric method was handed. Given the anisotropic operator, it discovers semi-coarsening, in the right direction, without a coordinate.

 1234567891010⁻⁵⁵10⁻⁴⁹10⁻⁴³10⁻³⁷10⁻³¹10⁻²⁵10⁻¹⁹10⁻¹³10⁻⁷10⁻¹cyclesizethe residualrepairedthe reported boundwhat the stopping rule readsreported at the last cycle9.4·10⁻⁴¹the residual there5.7·10⁻⁵with the border recomputed3.8·10⁻¹⁴products, cheap and repaired9a bound with nothing under itand one product a cycle to fix it Eigenvalues, singular values, rank

Keeping the vectors, and losing the bound

Thick restarting keeps the Ritz vectors instead of filtering the starting vector — the same eigenvalues for a third of the products with A. Its residual bound reaches 9.4·10⁻⁴¹ while the residual it bounds sits at 5.7·10⁻⁵, and the eigenvalues are correct to 4.3·10⁻¹⁴ the whole time, so nothing reports it.

 01234567110¹10²levelstored entries per rowalgebraicgeometrictwo complexities, one hierarchygrid complexity3.1operator complexity18geometric, for comparison1.5a ring plus random chords — nothing is anywherelevel 6 is 100% dense Iterating, instead of factorising

A hierarchy with no grid behind it

On a graph Laplacian the algebraic V-cycle converges at 0.199 a cycle, its grid complexity is an unremarkable 3.05, and its operator complexity is 17.7 — one level of forty-one unknowns is entirely dense. The number people quote is the one that does not measure the work.

All themes