Numerical Linear Algebra
where the algebra stops being the arithmetic
Theme

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

Essays 25 to 48 of the 54 on this theme, in the same order.
10²020406080100120unknownsiterations2D, no preconditioner2D, preconditioned1D, preconditionedone construction, two dimensions2D steps at 16 unknowns102D steps at 100 unknowns211D steps at 100 unknowns10the same averaging, the same transformand a count that no longer stops growing Structure, and the solver that cannot see it

Two dimensions, and the cluster that thins

The same kernel, the same averaging, the same transform — applied along two axes instead of one. In one dimension the preconditioned step count is 7, 10, 10, 10; on square grids with the same unknown counts it is 10, 18, 20, 21, and the share of the spectrum near one falls from 56% to 17%.

 159131721252910⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹steprelative error in the orthogonal factorNewtonNewton, scaledNewton–Schulzone fixed point, three costsscaled Newton, steps7Newton–Schulz, steps28Newton at step 657scaled Newton at step 64.4·10⁻¹³a Newton step needs an inverseand a Schulz step needs two products Orthogonality, measured

An iteration that only multiplies

Newton's iteration for the polar factor needs an inverse every step. Newton–Schulz needs only matrix products — nothing that reads an entry, nothing that pivots — and it converges if and only if every singular value is below √3. At 1.73205 it converges and at 1.73206 it returns an orthogonal matrix that is not the answer, with a residual of 5·10⁻¹⁶ and nothing to say so.

 the bound211no pivoting127τ = 0.00198τ = 0.003100τ = 0.01106τ = 0.03111τ = 0.1116τ = 0.3117τ = 1138entries in Ua bound, and its slackthe bound, from the graph alone211the worst that occurs138loose by1.5no arithmetic was done to compute the bound — only the pattern of AᵀA and its elimination graphGeorge and Ng: U fits inside chol(AᵀA), whatever the swapsprovable, cheap, and loose Sparsity, and what elimination costs

What the symbolic phase can only bound

Without pivoting, the fill can be computed from the graph and the count is exact — 233 predicted, 233 measured. With pivoting it is 233 predicted and 242 measured, and what survives is a bound that is right at every threshold and loose by 1.7 times at the largest grid drawn.

 3000000120000-21000000-0.5-1.500001.5-0.5000000-2T = ZᵀAZthe highlighted boxes each hold one conjugate pair, and no real rotation removes themthe form, and that it is one‖A − ZTZᵀ‖/‖A‖1.8·10⁻¹⁵‖ZᵀZ − I‖2.5·10⁻¹⁵worst eigenvalue error2.7·10⁻¹⁵surviving subdiagonal26×6, spectrum chosen before the matrix was builtquasi-triangular is as far as the reals go Eigenvalues, singular values, rank

The form a real matrix can reach

A real matrix with complex eigenvalues has no real triangular form, and the reason is one line — a real triangular matrix has a real diagonal, and a similarity does not move the spectrum. What it has instead is triangular except for one two-by-two block per conjugate pair, and the count is decided by the matrix rather than by where the iteration stopped.

 3456789101100.250.50.751bits in the exponent fieldeach curve as a fraction of its own maximumbfloat16fp16rangeprecisionto 617 decadesto 3.9 digitswhat the split buysbfloat16: largest number3.4·10³⁸fp16: largest number6.6·10⁴a bit of exponent doubles the rangea bit of significand adds a third of a digit The arithmetic underneath

The other half of a format

fp16 and tf32 have the same eleven significand bits and their largest numbers are 65,504 and 3.4·10³⁸. For two phases this site simulated the significand alone, so it was obliged to report them as the same format — which is a claim, and a false one.

 D from PAPᵀ = LDLᵀ — the shaded pairs are 2×2 pivots10⁻⁶0.749······0.749·········1.2·10⁻⁶1.4······1.4·········2.1·10⁻⁶0.549······0.549·········3.6·10⁻⁶0.614······0.614·three rules, one matrix‖PAPᵀ − LDLᵀ‖, blocks5.8·10⁻¹⁷‖PAPᵀ − LDLᵀ‖, diagonal3.1·10⁻¹¹growth, blocks1.3growth, diagonal5·10⁵the zero block is what the problem saysand one rule does not need it to be nonzero Elimination, and the swap

When symmetry is not enough

The matrix [[0, 1], [1, 0]] is symmetric, nonsingular and perfectly conditioned, and there is no diagonal entry to pivot on. Every factorisation restricted to symmetric interchanges and one-by-one pivots fails on it, at any depth of searching, because every entry it could search is zero. The repair is to take two variables at once.

 01234567810⁻⁵10⁻⁴10⁻³10⁻²10⁻¹110¹μ, the entry above the diagonalsep, and the eigenvalue gapmin |λᵢ + μⱼ|sep(A, B)the spectra never moveeigenvalue gap, throughout2sep at μ = 02sep at μ = 89.5·10⁻⁴amplification there1056solvability is the eigenvaluesand conditioning is not Structure, and the solver that cannot see it

An equation whose unknown is a matrix

AX + XB = C is linear in X, so it has a coefficient matrix, and writing it down is the obvious thing to do. At n = 100 that matrix has a hundred million entries for a problem with ten thousand unknowns, and the algorithm everybody uses instead never forms it. Its conditioning is not the eigenvalue gap either, which is the number a reader is invited to consult.

 15731pointsfinest gridone unknown — the recursion bottoms out in a divisionthe coarse operator, two ways‖RA_hP − A_2h‖/‖A_2h‖10⁻¹⁸unknowns / finest grid1.7cycle cost, in fine sweeps12each coarse point reaches three fine ones½, 1, ½ — and the restriction is its transpose Iterating, instead of factorising

The same problem on a coarser grid

Restriction, the coarse operator and interpolation are three matrices with nine distinct entries between them. Two of the three are each other's transpose, and their product with the fine operator is the coarse discretisation exactly — not approximately, entry for entry, at every level.

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

 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.

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

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

 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.

 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.

 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.

 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.

 171319253137434910⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²110²indexmagnitude|r_kk|σ_kone factorisation, two verdicts‖AP − QR‖/‖A‖1.1·10⁻¹⁶|r_nn|8.7·10⁻⁴σ_min3.7·10⁻¹²column interchanges0|r_nn| is never below σ_minso the cheap verdict errs one way only Eigenvalues, singular values, rank

The cheap rank and what it cannot see

Almost nobody computes singular values to decide a rank. The standard substitute is QR with column pivoting, read off the diagonal of R — and there is a triangular matrix on which the greedy rule makes no interchange at all, has no better column available at any step, and reports a matrix eight orders of magnitude further from singular than it is.

 00.250.50.75100.250.50.751xuexactcentral differencesupwindthe oscillation is exact‖Ax − b‖/‖b‖ for the central answer4.6·10⁻¹⁸values outside [0, 1]16worst excursion0.52the dashed lines are 0 and 1, which the equation guaranteesno solver was involved Iterating, instead of factorising

The stencil that is not symmetric

Past a cell Péclet number of exactly one — measured by bisection at 1.0000000000000002 — the central-difference solution of a convection–diffusion problem oscillates from point to point and leaves the interval the equation guarantees, at 16 of 31 grid points. It is the exact solution of its own linear system, to 4.6·10⁻¹⁸. No solver was involved.

 darker is smaller: 10⁻¹⁰, 10⁻⁶, 10⁻³, 10⁻¹one spectrum, two matricesspectral radius0.8reach of the 10⁻³ level1.3eigenvalues, all at0.8the circle is |z| = 1and every eigenvalue is well inside it Eigenvalues, singular values, rank

The eigenvalues that are not there

For a normal matrix the resolvent norm is exactly one over the distance to the nearest eigenvalue, so a picture of it carries nothing the spectrum did not. Move one entry above the diagonal and the region a perturbation of 10⁻⁸ can put an eigenvalue into stops being a disc and reaches out past the unit circle, while every eigenvalue stays at 0.8.

 00.250.50.75100.51xuexact & tunedupwindcentralone of these is exacttuned, worst nodal error2.4·10⁻¹⁷upwind, worst nodal error0.14central, points outside [0, 1]16exact at every nodeand only at the nodes Iterating, instead of factorising

The diffusion that makes the answer exact

Upwinding adds h/2 of artificial diffusion. Central differencing adds none. Add ε·ξ·Pe with ξ = coth(Pe) − 1/Pe and the computed solution is the exact one at every grid point, to 2.4·10⁻¹⁷ — at every Péclet number, on the problem it was derived from and on no other.

All themes