Numerical Linear Algebra
where the algebra stops being the arithmetic
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.
31 min read 26 figures Whose fault is itPrecision as the knobMeasured, not assertedAssertions that reject

Worth reading first: The zero you are allowed to write · The zero that means it is finished · Rank is a decision.

The three tolerances of the previous essay were offers. Accept this much error, save this much work; the curve is monotone in both directions, and the only way to be wrong about it is to choose a point nobody wanted.

A detection threshold is a different object, and the difference is not a matter of degree. Both of its directions are failures.

Declare a zero too early and an answer that was available is thrown away. Declare it too late and the next step divides by rounding error. There is no trade because there is nothing being bought — the constant is not purchasing work, it is deciding what the computation concludes.

Where conjugate gradients certifies that a positive definite matrix is indefiniteEvery matrix in this grid is 30×30 and positive definite by construction — its spectrum is κ^(−i/(n−1)), so the smallest eigenvalue is 1/κ and none of them is negative. A filled cell is a run in which pᵀAp came out non-positive and the iteration produced a direction it would report as a proof of indefiniteness. 14 of the 72 runs did. The region is a staircase whose top edge is at twelve significand bits and whose left edge, at eight, is at κ = 10⁵. Nothing rounded incorrectly anywhere: every comparison was performed exactly as written, on a number that was computed as accurately as the format allows.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
Fig. 1 Seventy-two runs, every matrix positive definite by construction. A filled cell is a run in which the iteration produced a direction it would report as a proof that the matrix is not.

One matrix, three answers

Start with the least alarming version. Build a symmetric matrix whose spectrum has three coarse clusters of width 10⁻⁴, each split into three finer ones of width 10⁻¹⁰. The matrix genuinely has a near-invariant subspace of dimension three, another of dimension nine, and an exactly invariant one of dimension twenty-seven. All three are real.

Run an Arnoldi recurrence and vary only the constant it compares its subdiagonal against.

Where the Arnoldi recurrence stops, against the tolerance it calls a zero, n = 27The spectrum has 3 coarse clusters of width 10^-4 and 9 fine ones inside them, so this one matrix has three near-invariant subspaces and all three are real. The tolerance decides which the recurrence stops at: it stops at 3 at the loose end, at 9 in the middle, and at 27 at the tight end. The residual it leaves runs from 1.51·10⁻⁵ to 5.48·10⁻¹⁷. Neither end is a mistake and neither is a trade: stopping too early throws away an answer that was available, and stopping too late means dividing by rounding error.10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²0612182430tolerance called a zerostep the recurrence stops atsubspace of 3subspace of 9subspace of 27too late: dividing by roundingtoo early: an answer thrown awayone matrix, three answerssteps at the tightest tolerance27at the loosest3residual left, tightest5.5·10⁻¹⁷residual left, loosest1.5·10⁻⁵both directions are failuresand the constant decides which
Fig. 2 Where the recurrence stops, against the tolerance it calls a zero. Three steps, and each is one of the three subspaces the matrix has.

It stops at step 3 at the loose end, leaving a relative residual of 1.5·10⁻⁵. It stops at step 9 in the middle, leaving 2.7·10⁻¹⁶. It stops at step 27 at the tight end, leaving 5.5·10⁻¹⁷.

None of those is wrong. Each is the recurrence correctly detecting one of the invariant structures the matrix has, and which one it detects is a property of the constant rather than of the matrix.

The two ends are the two failures, and they are worth naming as failures rather than as choices.

Too loose is an answer thrown away. Stopping at step 3 leaves a residual of 10⁻⁵ that the next six steps would have removed, on a problem where twenty-four steps remain available. The method has finished with 10⁻⁵ of the answer missing and reported success.

Too tight is a division by rounding. Past the point where the subdiagonal has fallen to the roundoff, the recurrence is normalising a vector that is entirely orthogonalisation error, and the “basis vector” it produces points in a direction determined by the arithmetic. Every subsequent step builds on it.

Where the Arnoldi recurrence stops, against the tolerance it calls a zero, n = 27The spectrum has 3 coarse clusters of width 10^-8 and 9 fine ones inside them, so this one matrix has three near-invariant subspaces and all three are real. The tolerance decides which the recurrence stops at: it stops at 3 at the loose end, at 9 in the middle, and at 27 at the tight end. The residual it leaves runs from 1.51·10⁻⁹ to 2.92·10⁻¹⁷. Neither end is a mistake and neither is a trade: stopping too early throws away an answer that was available, and stopping too late means dividing by rounding error.10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²0612182430tolerance called a zerostep the recurrence stops atsubspace of 3subspace of 9subspace of 27too late: dividing by roundingtoo early: an answer thrown awayone matrix, three answerssteps at the tightest tolerance27at the loosest3residual left, tightest2.9·10⁻¹⁷residual left, loosest1.5·10⁻⁹both directions are failuresand the constant decides which
Fig. 3 At a coarse width of 10⁻⁸ the first two steps of the staircase have moved close enough together that a single reasonable-looking constant can land on either.
Where the Arnoldi recurrence stops, against the tolerance it calls a zero, n = 27The spectrum has 3 coarse clusters of width 10^-2 and 9 fine ones inside them, so this one matrix has three near-invariant subspaces and all three are real. The tolerance decides which the recurrence stops at: it stops at 3 at the loose end, at 9 in the middle, and at 27 at the tight end. The residual it leaves runs from 0.00152 to 2.52·10⁻¹⁶. Neither end is a mistake and neither is a trade: stopping too early throws away an answer that was available, and stopping too late means dividing by rounding error.10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²0612182430tolerance called a zerostep the recurrence stops atsubspace of 3subspace of 9subspace of 27too late: dividing by roundingtoo early: an answer thrown awayone matrix, three answerssteps at the tightest tolerance27at the loosest3residual left, tightest2.5·10⁻¹⁶residual left, loosest0.0015both directions are failuresand the constant decides which
Fig. 4 And at 10⁻², where they are far apart and the choice is easy — which is the case a default is tuned on.
The Arnoldi subdiagonal and the GMRES residual against the width of an eigenvalue cluster, n = 30, 3 clustersA matrix of size 30 whose spectrum is 3 clusters of width η. At η = 0 the Krylov space closes at dimension 3 and the recurrence breaks down there. As η opens up, h_{4,3} rises with a fitted slope of 1.00 and the relative residual at step 3 with a slope of 1.00 — both exactly proportional to η, at a fixed ratio of 13.4 between them. The residual at step 2, one short of the breakdown, is 0.1264 and stays inside a per cent of that across the whole sweep: the step before the event carries no information about it.10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10²cluster width ηmagnitudeh at step 3residual at 3residual at 2fitted over ten decadesslope of h against η1slope of the residual1h ÷ residual, at every stop13residual at step 20.13the zero arrives in proportionand the step before it learns nothing
Fig. 5 The quantity being compared, from this collection’s first essay on the subject: the subdiagonal at the breakdown, proportional to the cluster width over ten decades.

A count that a constant moves by whole eigenvalues

The second instance is worse, because the quantity it decides is an integer.

The number of infinite eigenvalues of a pencil is n minus the degree of its characteristic polynomial, and on exact data that is a fact. On float data it is the number of singular values of B judged to be zero — a rank decision, and this collection has an essay saying at length that a rank is a decision.

Sweeping the tolerance on one perturbed descriptor pencil whose true count is two:

τ = 10⁻¹⁶ … 10⁻¹² count 0 τ = 10⁻¹⁰ count 1 τ = 10⁻⁸ … 10⁻² count 2 τ = 10⁻¹ count 3 τ = 3·10⁻¹ count 4

One matrix, five different answers, and the differences are not digits. A pencil judged to have two infinite eigenvalues when it has three has one spurious finite eigenvalue in its output, at a value that is the ratio of two quantities which are both rounding error — and there is nothing about its residual that says so, because a small residual is available for it.

The singular values of B, and the number of infinite eigenvalues the exact polynomial says there areThe pencil has 2 algebraic constraints, so 2 of its 6 eigenvalues are infinite. Two routes say so and they are not the same kind of statement. The exact route computes det(A − λB) in BigInt rationals and finds it has degree 4, so n − 4 = 2 eigenvalues are at infinity — an integer, arrived at without rounding. The float route counts the singular values of B below a cut, and the values run 4.27, 2.6, 0.918, 0.416, 1.13·10⁻¹⁶, 4.1·10⁻¹⁵⁷: the gap is a factor of 2.75·10¹⁴⁰ and it falls in the right place, so the decision is comfortable here. It is still a decision, and it is the same decision this site's rank essay is about.σ14.27σ22.6σ30.918σ40.416σ51.13·10⁻¹⁶σ610⁻¹⁸n = 6, and det(A − λB) has degree 4an integer, and a judgementdegree of det(A − λB), exactly4infinite eigenvalues, from the degree2singular values below the cut2largest gap in the spectrum2.7·10¹⁴⁰a degree cannot be nearly threeand a singular value can be nearly zero
Fig. 6 The exact count beside the decision, from the pencil essays: an integer arrived at without rounding, and a spectrum of singular values with a cut somewhere in it.
Singular values of a rank-4 matrix with noise of relative size 10⁻⁶Ten singular values on a logarithmic axis. The first four sit near one; the rest sit at the noise level, and the vertical distance between the two groups is the evidence for the rank.1234567891010⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1indexsingular valuecutoff, σ₁ · 10⁻¹⁰numerical rank 10gap 8.2·10⁴an opiniontrue rank 410×10, built with 4 nonzero valuesrank is a decision
Fig. 7 The decision in its original setting on this site, where the evidence for it is a gap.
A randomised rank-k solve, 4 seeds a rankRelative error against the rank kept, on a logarithmic vertical axis, with a vertical bar at each rank spanning the seeds. The band is 1.68 wide at rank 8, where the method is at its worst, and 1.004 wide at rank 20, where it is at its best. The same computation on the same data returns a different answer each time.0481216202428321rank keptrelative errormedianthe answer movesspread at rank 81.7spread at rank 201best median error0.15widest where the method is worstand the bound does not say so
Fig. 8 And the honest form of such a decision, which is a band rather than a verdict.

There is a band of tolerances that gets it right, and it is wide — six decades. That is the ordinary case and it is why this works most of the time. What the sweep shows is that the width of the band is a property of the matrix, that nothing in the output reports how close to an edge of it the chosen τ was, and that the failure at each edge is silent.

Why the band is not a safety margin

There is a tempting reading of the pencil sweep above — six decades of tolerance give the right answer, so the decision is robust — and it is worth taking apart, because it is the reading that makes this kind of failure surprising when it arrives.

The width of the band is set by the gap in the singular values of B: it runs from the size of the perturbation up to the smallest genuine singular value. Both ends are properties of the problem. On the constructed pencil here the gap is eight orders wide and the band is comfortable; on a pencil assembled from measured data it is whatever the measurement noise and the model leave, and it can be nothing at all.

What makes it dangerous is that the width is not reported. A routine given τ = 10⁻⁸ returns a count, and returns exactly the same count in the comfortable case and in the case where τ was a factor of two from the edge. The number that says which — the ratio of the largest singular value below the cut to the smallest above it — is computed on the way to the answer and discarded.

That is the same complaint as the previous essays’ and it is the reason this collection’s own rank machinery returns a gap and a verdict rather than a verdict.

The diagonal of a pivoted R against the singular values, on a 40×40 Kahan matrix at c = 0.5Two curves down the index. The diagonal of R decays at one constant rate to 0.00366; the singular values follow it and then fall away to 9.07·10⁻¹⁰ at the last one. Column pivoting made no interchange at all, because every trailing column norm is equal at every step — the rule had nothing to choose between.1611162126313610⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²110²indexmagnitude|r_kk|σ_kone factorisation, two verdicts‖AP − QR‖/‖A‖10⁻¹⁶|r_nn|0.0037σ_min9.1·10⁻¹⁰column interchanges0|r_nn| is never below σ_minso the cheap verdict errs one way only
Fig. 9 A cheaper verdict on the same question from the rank-revealing essays, which errs in the flattering direction.
How far |r_nn| sits above σ_min on Kahan's matrix, against the size and the parameter3 curves of |r_nn| ÷ σ_min against n, one per Kahan parameter. Every curve rises without turning over, reaching 2.3·10⁸ at n = 50, c = 0.5. Column pivoting makes no interchange at any point on any of them, so the failure is not a poor choice — there is nothing to choose.81624324048110²10⁴10⁶10⁸size of the matrix|r_nn| ÷ σ_minthe two agreec = 0.2c = 0.35c = 0.5no ceilingratio at n = 1021ratio at n = 502.3·10⁸interchanges, anywhere0the greedy rule never had a choiceand the gap grows with every row
Fig. 10 And the matrix built to make it err, which is what a band with no margin looks like.

And a proof of something false

The third instance is the one worth the essay’s title.

Conjugate gradients decides a matrix is indefinite by testing the sign of pᵀAp. That is not a tolerance at all — the constant is zero, the comparison is exact, and there is no judgement in it. It is the cleanest possible detection: a sign.

A sign of a computed quantity is unreliable when the quantity is smaller than its own rounding error. pᵀAp/pᵀp is a Rayleigh quotient, so it is bounded below by λ_min; the rounding in the inner product is about u·‖A‖. So the sign means nothing once λ_min/λ_max falls below roughly u·n — and past that point, a positive definite matrix produces a direction that certifies it is not.

The hero is that region. Every matrix in it has spectrum κ^(−i/(n−1)), so its smallest eigenvalue is 1/κ and none of them is negative. The filled cells are runs where the iteration produced a certificate, and the number in each is the step at which it did.

The boundary sits at twelve significand bits. At fourteen and above, over eight decades of conditioning, the test is never wrong. At twelve and below it is wrong somewhere, and at eight it is wrong from κ = 10⁵ — which is not an extreme condition number by any standard.

Where conjugate gradients certifies that a positive definite matrix is indefiniteEvery matrix in this grid is 16×16 and positive definite by construction — its spectrum is κ^(−i/(n−1)), so the smallest eigenvalue is 1/κ and none of them is negative. A filled cell is a run in which pᵀAp came out non-positive and the iteration produced a direction it would report as a proof of indefiniteness. 18 of the 72 runs did. The region is a staircase whose top edge is at twelve significand bits and whose left edge, at eight, is at κ = 10⁶. Nothing rounded incorrectly anywhere: every comparison was performed exactly as written, on a number that was computed as accurately as the format allows.significand bits10³10⁴10⁵10⁶10⁷10⁸10⁹10¹⁰κ(A)8136460199426314512345105113211154244912624014162024every matrix positive definiteruns producing a false certificate18of runs in total72never above, in significand bits12first κ at eight bits10⁶the comparison was correctand what it proved was not true
Fig. 11 At sixteen unknowns the region is smaller: a shorter iteration is fewer chances for a search direction to become one the rounding cannot keep positive.
Where conjugate gradients certifies that a positive definite matrix is indefiniteEvery matrix in this grid is 48×48 and positive definite by construction — its spectrum is κ^(−i/(n−1)), so the smallest eigenvalue is 1/κ and none of them is negative. A filled cell is a run in which pᵀAp came out non-positive and the iteration produced a direction it would report as a proof of indefiniteness. 18 of the 72 runs did. The region is a staircase whose top edge is at twelve significand bits and whose left edge, at eight, is at κ = 10⁶. Nothing rounded incorrectly anywhere: every comparison was performed exactly as written, on a number that was computed as accurately as the format allows.significand bits10³10⁴10⁵10⁶10⁷10⁸10⁹10¹⁰κ(A)8674344381691001014868201017173138412711170611210014162024every matrix positive definiteruns producing a false certificate18of runs in total72never above, in significand bits12first κ at eight bits10⁶the comparison was correctand what it proved was not true
Fig. 12 And at forty-eight it is larger. The certificate gets less trustworthy as the problem gets bigger, which is the direction that matters, because bigger is why anybody is using an iterative method.

Twelve bits is not an exotic place to be

Eight significand bits is the E4M3 format that this collection has an essay about, and half the hardware sold for numerical work now implements it. Eleven is half precision. Twelve is between them.

None of this is a warning about a museum piece. It is a statement that the arithmetic the field is moving towards has a unit roundoff at which a test everybody treats as exact stops being exact, and that the failure mode is a proof rather than a wrong number.

Every value E4M3 can hold, and the one it cannotA logarithmic axis with a vertical tick at every representable magnitude, crowding together towards the left and thinning to the right, with the top of the range marked.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
Fig. 13 The eight-bit format in question, from the arithmetic field, where its exponent and significand are drawn out.
Precision against range, with the 4-bit-exponent formats markedA scatter of formats with unit roundoff across and largest finite value up, both logarithmic. The two eight-bit formats sit at the lower left, joined to each other and to a marked position below one of them.10⁻¹⁶10⁻¹²10⁻⁸10⁻⁴10²10⁸10¹⁴10²⁰10²⁶10³²10³⁸unit roundofflargest finite valueE4M3E5M2bfloat16fp16tf32binary32binary64240 — what IEEE's rules would givethe octave the rule-breaking boughtE4M3 largest finite448under IEEE's rules240E5M2 largest finite5.7·10⁴seven bits and a sign, spent two waysno infinities in one of them
Fig. 14 And the trade its designers made between range and precision, which is where the twelve bits above come from.
Forward error of a 8×8 Hilbert solve at eight precisionsA bar for each significand width from 12 to 53 bits showing the relative error in the computed solution, with the condition number times the unit roundoff marked as a prediction.κ = 1.5·10¹⁰ · the exact answer is (1, 2, …, 8)12 bits4.116 bits1.620 bits1024 bits0.530 bits0.5636 bits0.005943 bits5.5·10⁻⁶53 bits3.2·10⁻⁸dashed: κ · unit roundoffone matrix, eight arithmeticsmeasured against a known answer
Fig. 15 The site’s own knob, which is why this measurement was available at all: the same computation at every significand width from eight bits to double.

The mechanism, in one line

It is worth writing out why twelve rather than sixteen or eight, because the boundary is not arbitrary and the argument gives a rule.

The quotient being tested is q = pᵀAp/pᵀp, whose true value is at least λ_min. The rounding in a dot product of length n is about n·u·‖p‖‖Ap‖, and ‖Ap‖ ≤ λ_max‖p‖, so the computed q is q plus something of size n·u·λ_max. The sign is unreliable once

λ_min ⁄ λ_max < n · u

which is 1/κ < n·u, or κ > 1/(n·u).

At fifty-three bits and n = 30 that is κ > 5·10¹⁴. At twenty-four bits, κ > 2·10⁶. At eight bits, κ > 1.7·10². The measured boundary is looser than that at every precision, because the conjugate gradient directions do not reach the smallest eigenvalue immediately and a shorter run never gets near it — which is the same effect that made the certificate hard to find in the essay on negative curvature, arriving here as protection rather than as cost.

The rule that falls out is worth carrying: a sign test on a Rayleigh quotient is trustworthy when κ·n·u is well below one, and is a coin flip when it is not. It costs nothing to check and it is not checked anywhere.

Relative error of two algebraically identical expressions for (1 − cos x)/x², in 11-bitA log–log plot of relative error against x. The expression written as it reads loses accuracy as x falls and is eventually wrong in every digit; the rearranged form stays at rounding level.10⁻⁶10⁻⁴10⁻²110⁻⁵10⁻²10¹xrelative error of the computed value(1 − cos x)/x², as written2 sin²(x/2)/x²no digits left at all11-bit throughoutone function, two spellings · zero below 0.031
Fig. 16 The mechanism at eleven significand bits, from the arithmetic field, where a subtraction of nearly equal quantities leaves nothing behind.
A 6.25-bit block format against 8-bit E4M3, read two waysFour curves of relative error against the spread of the data. Two, for the block format, rise steeply when read as a median and gently when read as a norm; the other two are nearly flat.048121610⁻²10⁻¹1octaves of spread within a blockrelative errorblock, medianblock, 2-normE4M3, medianE4M3, 2-normthe same data, two readingsblock format, bits a value6.3norm across 16 octaves1.8median across 16 octaves62entries deleted at the wide end383640 values, blocks of 32a norm is dominated by what a block format keeps
Fig. 17 And what else changes at that width, from the low-precision essays.

What is actually being decided

Put the three instances beside each other and the shape they share is clear.

In each of them, a quantity is compared against a constant and the conclusion changes. Not the accuracy — the conclusion. A recurrence stops at three or at nine or at twenty-seven; a pencil has one or two or four infinite eigenvalues; a matrix is positive definite or it is not.

Those are discrete outcomes reached by comparing a continuous quantity against a constant, which is the definition of a decision. And a decision has a property that a trade does not: it can be wrong, and the wrongness is not measured in digits.

That is the whole difference from the previous essay. A tolerance that buys work costs accuracy in proportion, so a bad choice is a suboptimal point on a curve. A tolerance that decides costs a conclusion, so a bad choice is a different answer to a question with only one right one.

Smoothed aggregation against the angle of the anisotropy (ε = 0.01)Convergence factor against the rotation angle, with the stencil's largest axis coupling and its diagonal coupling on the same axis. Aligned, the factor is 0.1927. At 45° it is 0.7774, and the reason is beside it: the axis couplings are 0.5050 while the coupling along the direction the anisotropy runs in is 0.2475. The matrix does not contain the anisotropy.-213284300.250.50.751rotation of the anisotropy (degrees)factor / couplingusableconvergence factoraxis couplingdiagonal couplingthe standard answer, and the anglefactor at 0°0.19factor at 45°0.78axis ÷ diagonal coupling at 45°2the hierarchy reads the matrixand the matrix lost the direction
Fig. 18 A decision of the same kind in the multigrid field, where a threshold decides which directions the method calls strongly connected.
Two subdivision rules on roots 0.2 apart, one of which sits on the first cutSix bars. Cutting a box at its midpoint puts the root at x = 1 on a boundary of every box the search ever asks about, and the operator cannot verify a root it never sees the interior of: 0 verified after 147 evaluations. Cutting at 0.485 of the width verifies both roots in 57.both rules see the same problem and the same operatormidpoint: verified0midpoint: undecided6midpoint: evaluations147off-centre: verified2off-centre: undecided0off-centre: evaluations57where the box is cutroots verified, midpoint0roots verified, off-centre2evaluations saved90a root on a boundary is never in an interiorand 0.485 is the whole repair
Fig. 19 And another, from the branch-and-bound essay, where the rule decides rather than computes.

What can be done about them

Three things, in increasing order of how rarely they are done.

Report which side of the decision was close. Every one of these thresholds has a quantity beside it, and the ratio of that quantity to the threshold is a free measure of how confident the decision was. A rank decision backed by a gap of 10¹⁴ and one backed by a gap of 3 are different objects; this collection’s own rank machinery returns the gap alongside the count for exactly that reason, and almost nothing else does.

Report what the decision was about. A breakdown detected at step k should say whether the vectors were small or the inner product was — which the previous essays in this run showed is the difference between finishing and failing, costs two norms, and cannot be recovered afterwards.

And verify a conclusion in the precision the claim is made in. A certificate of negative curvature costs one matrix–vector product to check from outside, and the check is independent of the iteration that produced it. At eight bits it is the only thing standing between a proof and a rounding error.

The Krawczyk operator on a box of half-width 0.4Two rectangles in the plane. The outer one is the box handed in; the inner one is K(X), the image the operator returns; and the marked point is the root (√2, √2), known in closed form. The image is narrower than the box — 0.2265 against 0.8 — and the verdict is: exactly one root.11.522.511.52xyexactly one roota verdict, not a bound‖I − C F′(X)‖0.28width of X0.8width of K(X)0.23strictly inside is a proofand overlapping is nothing at all
Fig. 20 The strongest form of that habit on this site: a computation that returns an interval containing the answer, which anybody can check without repeating the work.
Past the cliff: what each method returnsLoss of orthogonality at four condition numbers past κ²u = 1. Cholesky QR returns a factor whose implied Q is 0.37 away from orthogonal at κ = 10⁸ and 1.49 at 10¹² — and at 3·10¹⁰ it refuses outright, on a pivot of -6.5·10⁻¹⁷. The sweep and the tree return a usable factorisation of every one of these matrices.‖QᵀQ − I‖ of the implied Qκ = 10⁸, Cholesky0.37κ = 10⁸, sweep8.5·10⁻⁹κ = 10¹⁰, Cholesky1.3κ = 10¹⁰, sweep3·10⁻⁷κ = 3·10¹⁰, Choleskyrefusedκ = 3·10¹⁰, sweep6.7·10⁻⁶κ = 10¹², Cholesky1.5κ = 10¹², sweep1.4·10⁻⁴the safe run is the one that failsrefusals in the range1wholly non-orthogonal returns2the pivot it refused on-6.5·10⁻¹⁷one of these outcomes is safeand it is the refusal
Fig. 21 And the general rule this collection keeps arriving at: a routine that cannot do what was asked has two options, and returning an answer is the worse one.
Where a bound can be proved, against κuEvery combination of 4 working precisions and 5 Hilbert sizes, placed by the product of the condition number and the unit roundoff. Filled marks are the cases where a bound was proved; open ones are the cases the method refused. The two sets separate at κu ≈ 1 — the largest verified is 0.24 and the smallest declined is 0.89 — which is the threshold iterative refinement's convergence sits at.10⁻⁵10⁻²10¹10⁴10⁷10¹⁰132129κ · usignificand bitsκu = 1a bound was provedthe method refusedit never returns a wrong boundlargest κu with a proof0.24smallest κu without one0.89cases refused, of the grid15a refusal is not a wide bound — it is no bound at alland it is the only failure mode here
Fig. 22 Where declining is the right answer, from the verification essays: the boundary past which the method refuses rather than reports.

A fourth threshold, which this essay’s own machinery contains

Every essay in this run is drawn by a generator, and several of those generators contain a threshold of exactly the kind being described. It is worth admitting one.

The Arnoldi recurrence in this collection’s own library declares a breakdown when the subdiagonal is below 10⁻¹⁰ times the Frobenius norm of the matrix. That constant was not derived. It was raised from 10⁻¹³ while these figures were being drawn, because at 10⁻¹³ a spectrum with twelve distinct eigenvalues did not break down at all — the measured subdiagonal there is 5·10⁻¹², which is above the threshold — and the figure that was meant to show the recurrence stopping at twelve showed it running past.

The value it was raised to is nine orders below the smallest genuine subdiagonal in any of the spectra drawn, so the band is wide and the choice is safe. It is still a choice, it is still a constant somebody typed, and the essay it enables is about exactly that.

Recording it here rather than in a comment is the point. A tolerance defended in the prose of the argument it supports is a tolerance a reader can disagree with; one buried in a source file is a number nobody will ever look at again.

The step at which the Arnoldi recurrence breaks down, against the number of distinct eigenvalues, n = 30Every matrix here is 30×30 and every one of them breaks down before step 18. The step is the number of distinct eigenvalues in the spectrum, exactly, at every one of the 8 spectra drawn — the points lie on the diagonal and the horizontal line at n = 30, which is the bound every course states, is off the top of the picture. The relative residual at the breakdown step is at most 5.64·10⁻¹⁶; one step earlier it is at least 3.5·10¹⁰ times larger.024681012141618024681012141618distinct eigenvalues in the spectrumstep the recurrence stops atthe step is m, not nn = 30 throughoutspectra drawn8every one breaking at m8worst residual at the breakdown5.6·10⁻¹⁶smallest gain over the step before3.5·10¹⁰an invariant subspace contains the answerand its dimension is what the method costs
Fig. 23 The figure that needed the constant raised, from the first essay of this run: the breakdown step against the number of distinct eigenvalues.

The refusal

The assertion is fed the claim that a false certificate appears at double precision.

It is the refusal that keeps the finding honest in the flattering direction. Everything above says the curvature test can be wrong; a routine that reported it wrong at fifty-three significand bits over this range of conditioning would be over-stating the result, and the whole point of drawing a boundary is that there is something on the other side of it.

Over the three conditionings tested at double precision — 10⁶, 10⁸ and 10¹⁰ — there is not one false certificate. The boundary is where the grid says it is, and the assertion fails if the machinery ever starts finding one there.

Where this run of essays ends

Twelve essays, and one sentence underneath them.

A quantity going to zero is one arithmetic event with two opposite meanings — an object was found, or an object was lost — and nothing in the number says which. An Arnoldi subdiagonal of 10⁻¹⁴ is the answer arriving; a two-sided divisor of 10⁻¹⁴ is a method with nowhere to go. A determinant that is identically zero is a question with no content; a determinant that has merely lost degree is a model with constraints in it. A curvature that changes sign is the most useful object an iteration can produce, or a rounding error wearing a proof’s clothes, and which of those it is depends on the significand width.

The last two essays are the reason the sentence is about programs rather than about mathematics. In floating point no quantity ever is zero, so every one of those events is somebody’s comparison against a constant — and the constant decides which of the two meanings the computation acts on.

Where the Arnoldi recurrence stops, against the tolerance it calls a zero, n = 27The spectrum has 3 coarse clusters of width 10^-6 and 9 fine ones inside them, so this one matrix has three near-invariant subspaces and all three are real. The tolerance decides which the recurrence stops at: it stops at 3 at the loose end, at 9 in the middle, and at 27 at the tight end. The residual it leaves runs from 1.51·10⁻⁷ to 2.85·10⁻¹⁷. Neither end is a mistake and neither is a trade: stopping too early throws away an answer that was available, and stopping too late means dividing by rounding error.10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²0612182430tolerance called a zerostep the recurrence stops atsubspace of 3subspace of 9subspace of 27too late: dividing by roundingtoo early: an answer thrown awayone matrix, three answerssteps at the tightest tolerance27at the loosest3residual left, tightest2.8·10⁻¹⁷residual left, loosest1.5·10⁻⁷both directions are failuresand the constant decides which
Fig. 24 At a coarse width of 10⁻⁶ the first two steps of the staircase are a decade apart, which is the comfortable case and the one a default is chosen in.
Where conjugate gradients certifies that a positive definite matrix is indefiniteEvery matrix in this grid is 40×40 and positive definite by construction — its spectrum is κ^(−i/(n−1)), so the smallest eigenvalue is 1/κ and none of them is negative. A filled cell is a run in which pᵀAp came out non-positive and the iteration produced a direction it would report as a proof of indefiniteness. 11 of the 72 runs did. The region is a staircase whose top edge is at twelve significand bits and whose left edge, at eight, is at κ = 10⁵. Nothing rounded incorrectly anywhere: every comparison was performed exactly as written, on a number that was computed as accurately as the format allows.significand bits10³10⁴10⁵10⁶10⁷10⁸10⁹10¹⁰κ(A)8103107722995729138136951015386111214162024every matrix positive definiteruns producing a false certificate11of runs in total72never above, in significand bits12first κ at eight bits10⁵the comparison was correctand what it proved was not true
Fig. 25 And at forty unknowns, where the region has grown into the precisions a working code would use.
Where the Arnoldi recurrence stops, against the tolerance it calls a zero, n = 27The spectrum has 3 coarse clusters of width 10^-3 and 9 fine ones inside them, so this one matrix has three near-invariant subspaces and all three are real. The tolerance decides which the recurrence stops at: it stops at 3 at the loose end, at 9 in the middle, and at 27 at the tight end. The residual it leaves runs from 1.51·10⁻⁴ to 5.26·10⁻¹⁷. Neither end is a mistake and neither is a trade: stopping too early throws away an answer that was available, and stopping too late means dividing by rounding error.10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²0612182430tolerance called a zerostep the recurrence stops atsubspace of 3subspace of 9subspace of 27too late: dividing by roundingtoo early: an answer thrown awayone matrix, three answerssteps at the tightest tolerance27at the loosest3residual left, tightest5.3·10⁻¹⁷residual left, loosest1.5·10⁻⁴both directions are failuresand the constant decides which
Fig. 26 At 10⁻³ the three steps are as far apart as this family puts them.

Where this goes next

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 The zero you are allowed to write A deflation criterion sets a subdiagonal entry to zero because it is small. A drop tolerance discards an entry of a factor because it is small. A truncation discards a singular value because it is small. Three fields, three vocabularies, no shared arithmetic — and plotted as work saved against error accepted, one curve. 1234567891010⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1indexsingular valuecutoff, σ₁ · 10⁻¹⁰numerical rank 10gap 8.2·10⁶an opiniontrue rank 410×10, built with 4 nonzero valuesrank is a decision Rank is a decision A floating-point matrix does not have a rank. It has a spectrum of singular values, and somewhere in that spectrum is a place where the values stop being signal and start being noise. Deciding where is a judgement, and the evidence for it is a gap. 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 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.

Shares its objects with

Essays that name at least two of the same things, and that neither author linked.

Named objects

A flat tag is an object no other essay names yet.

Backward errorCertificateCounterexampleHalf precisionInfinite eigenvalueInvariant subspaceKrylov subspaceMachine epsilonNegative curvatureNumerical rankTolerance