Numerical Linear Algebra
where the algebra stops being the arithmetic
The answer that depends on the machine

A matrix that is definite on one machine

Two hundred Gram matrices, two conforming builds, and twenty-six of them get different answers to "is this positive definite". The exact verdict, from determinants in BigInt rationals, says the fused build is right nine times and the other one seventeen.
33 min read 28 figures Whose fault is itMeasured, not assertedTwo routes to a numberExact ground truth

Worth reading first: One multiply the compiler removed · A factorisation with nothing to pivot for.

Is this matrix positive definite is the one question in the subject with a two-word answer. There is no tolerance in it, no unit, and no scale: a matrix is positive definite or it is not, and a Cholesky factorisation answers by completing or by stopping at a pivot that is not positive.

Which makes it the sharpest possible test of the previous essay’s claim. A determinant can be a little wrong. A verdict cannot.

Two conforming builds, one matrix, and 26 different verdicts in 200 triesForty 6 × 6 Gram matrices at each of five conditionings, each formed exactly from a tall matrix whose last column is nearly the sum of its first two. Every one of them is asked *are you positive definite* by a Cholesky, once with the inner products accumulated with a single rounding per term and once with each product rounded before it is subtracted. The bar is how many of the forty came back with different verdicts: 7 at 10⁻⁸, 2 at 10⁻⁹, 4 at 10⁻¹⁰, 10 at 10⁻¹¹, 3 at 10⁻¹². The two shaded parts split those by which build was right, settled by the signs of the leading principal minors computed in BigInt rationals — 9 to the fused form and 17 to the rounded one. So neither is the correct implementation. The verdict is a property of the build, and the answer to which build is right is a property of the matrix.-12-11-10-9-8024681012log₁₀ of how nearly dependent the columns areverdicts that disagreed, of 40724103which one is correctmatrices tested200verdicts disagreed26fused was right9unfused was right17lower part: thefused build was righta sign has no last digitso a verdict has nowhere to hide
Fig. 1 Forty 6 × 6 Gram matrices at each of five conditionings. The bar is how many of the forty received different verdicts from two conforming builds; the shaded part is how many of those the fused build got right.

Two hundred matrices, twenty-six disagreements, and the shaded parts are not the whole bars.

The experiment

Each matrix is a Gram matrix AᵀA, formed from a tall matrix whose last column is nearly the sum of its first two. Such a matrix is positive definite whenever the columns are independent, which they are by construction — but only just, and how just is the horizontal axis: the perturbation that separates the last column from the sum of the other two runs from 10⁻⁸ down to 10⁻¹².

The entries of AᵀA are formed exactly, by an error-free product expansion and an exactly rounded sum, so that the matrix handed to both routines is the same matrix and any disagreement afterwards belongs to the factorisation rather than to the assembly.

Then each matrix is asked the question twice.

One build accumulates each inner product with a single rounding per term — the fused form, which is what a compiler emits from s -= L[i][k]*L[j][k] when it contracts.

The other rounds each product before subtracting it — the unfused form, which is what the same source produces on a machine without the instruction, or with the contraction disabled.

Nothing else differs: the same algorithm, the same order, the same precision, the same matrix.

The adjudicator

A disagreement between two routines is only interesting if there is a third thing that settles it, and here there is one that admits no argument.

Sylvester’s criterion: a symmetric matrix is positive definite exactly when every leading principal minor is positive. Every entry of the matrix is a binary64 number and therefore a rational with a power of two underneath it, so its minors are rationals, and their signs are facts rather than measurements. This site’s exact.js computes them in BigInt, which is the same machinery the exact answer, in rationals uses for a Hilbert inverse.

So each of the two hundred matrices has a true verdict, computed by a route that shares no arithmetic with either build.

The result

Twenty-six of the two hundred came back with different verdicts. On each of those, exactly one of the two builds was right — which is a check rather than an observation, and is asserted on every build of this site.

The fused build was the right one nine times. The unfused build was right seventeen.

That is the finding, and it is the opposite of the tidy story. The previous essay’s determinant is exact under the fused form and zero under the other, which invites the conclusion that fusing is correct and contracting is a hazard. It is not: on a computation nobody designed around the instruction, the fusion moves the answer by one rounding in a direction with no reason to be right, and the census says it is right about a third of the time.

The honest statement is therefore not use the fused form. It is:

The verdict is a property of the build, and which build is right is a property of the matrix.

Why a definiteness test is where this shows

Three reasons, and they compound.

The quantity is a difference of nearly equal things. The last pivot of a Cholesky is Aⱼⱼ − Σ L²ⱼₖ, and for a matrix that is barely definite the two are nearly equal, so the summation condition number of that inner product is enormous — which is exactly where the vector that hides it puts every quantity a decision is read from.

The comparison is with zero. A threshold has a margin and a sign does not. Whatever the disagreement between two builds, if it straddles zero the verdicts differ, and the disagreement is guaranteed to straddle zero when the true pivot is smaller than the disagreement.

The matrix is a Gram matrix. Forming AᵀA squares the condition number, which the road that squares the problem is the essay about. So a matrix whose columns are separated by 10⁻¹⁰ arrives at the factorisation with a condition number near 10²⁰, and its smallest pivot is below the rounding of the entries — which is the regime the census is drawn in.

The third one is worth dwelling on because it is not a contrivance. A Gram matrix is what a normal equations solve, a kernel method, a covariance estimate and a least-squares fit all produce, and “nearly dependent columns” is the ordinary situation in all four.

What the disagreement rate is a function of

The bars are not uniform across the horizontal axis and the shape is worth reading rather than averaging away: 7 disagreements at a separation of 10⁻⁸, 2 at 10⁻⁹, 4 at 10⁻¹⁰, 10 at 10⁻¹¹ and 3 at 10⁻¹². The rate does not climb monotonically as the matrices get harder, which is what a first guess would predict.

The reason is that two things move in opposite directions. As the columns become more nearly dependent, the true smallest pivot falls — so it is more often below the rounding, and the two builds are more often deciding a sign from noise. But past a point both builds are deciding from noise, and two coin flips agree half the time. The disagreement rate is highest in the band where the pivot and the rounding are comparable, and falls off on both sides of it: above, because both builds are right; below, because both are guessing and guessing agrees as often as not.

That band is narrow — about two decades wide here — and its position is set by the size of the matrix and the scale of its entries rather than by anything a caller controls. Which is the practical reason the phenomenon is under-reported: a test suite that happens to sample outside the band sees perfect agreement and concludes there is nothing to see, and the same suite on a slightly different problem sees a third of its matrices flip.

What the verdict is actually reporting

Here is the reading that makes the census less alarming and more useful.

When the two builds disagree, the true smallest pivot is smaller than the rounding either build commits. In that regime the question is this matrix positive definite has an exact answer — the minors settle it — and it is a question about the matrix’s last few bits rather than about anything a computation is going to use.

A matrix with a smallest eigenvalue of 10⁻²⁰ relative to its largest is definite, and every computation done with it will behave as though it were singular. Its Cholesky factor, if one is produced, has a last column of noise. The definiteness verdict is true and it is not useful, which is a distinction the subject has a standing answer for: report the pivot, not the verdict.

That is what makes this essay a companion to rank is a decision rather than a scandal. Rank is not a property a floating-point matrix has; it is a decision about a gap, and the honest output of a rank computation is the spectrum with the threshold drawn on it. Definiteness is the same object with the threshold nailed to zero, and nailing a threshold to zero is what makes it a machine-dependent decision rather than a reported one.

The repair, which is not to fix the arithmetic

Three options, and the third is the recommendation.

Compute the pivots exactly. The minors in BigInt settle it, at a cost that is polynomial and enormous — this site can afford it up to about 10 × 10 and not beyond. It is the right instrument for adjudicating a census and the wrong one for a solver.

Make the inner products order-independent. Pre-rounded accumulation removes the disagreement between partitionings, but not between the fused and unfused forms — those differ in how a single product is rounded, which is upstream of the summation policy. So this fixes one of the two causes in the field and not this one.

Report the smallest pivot beside the verdict, and let the caller set a threshold. This is what LAPACK’s info plus the diagonal of the factor already allows and what almost no caller reads. A routine that returns definite, smallest pivot 3.1·10⁻¹⁷, matrix norm 4.2 has said everything, and a caller who wants a machine-independent answer thresholds the pivot at something above the rounding.

The third option is free, and it converts a verdict two builds disagree about into a number they agree about to within a factor.

Where else a sign is a verdict

The census is about definiteness because definiteness is the cleanest case, but the same shape appears wherever a computed sign chooses a branch, and the field’s other essays are the same argument on different objects:

Four verdicts, four different mechanisms, and one common feature: each turns a real number that two builds compute differently into a discrete outcome that has no room to differ a little.

A note on what was held fixed

The census is careful about one thing that is easy to get wrong, and the care is the reason its number means anything.

Both builds are handed the same matrix, and that matrix’s entries are computed exactly. If the Gram matrix were formed in ordinary arithmetic, the fused and unfused builds would each assemble a slightly different AᵀA, and a disagreement afterwards would be a disagreement about two different matrices — which is a much weaker and less interesting statement, and is the one an unguarded experiment would have measured.

Fixing the matrix costs an error-free product expansion per entry and it moves the question from do two builds compute different matrices to do two builds reach different verdicts about one matrix. Only the second is a statement about the factorisation, and only the second has an exact adjudicator, since the exact minors are minors of a matrix that has to be the same one.

The same discipline is why the adjudicator uses rationals rather than a longer float. A quadruple precision Cholesky would agree with one of the two builds most of the time and would be a third opinion rather than a verdict — and this site’s habit is that two routes to a number means two routes that share no arithmetic, not two precisions of one.

The one that cannot be repaired at all

A closing observation about the limit of the recommendation.

Report the number and let the caller threshold it works when there is a number to report. It does not work when the branch is inside a library the caller cannot see — when a factorisation routine tries a 2 × 2 pivot, or a solver switches to a different algorithm on detecting near-singularity, or an eigensolver deflates a converged value. Those decisions are made from computed quantities, taken inside a routine, and never surfaced.

That is the position most callers are in most of the time, and it is why this field’s practical advice ends where the library’s interface does. What can be measured from outside is the spread of the final answer, and what can be asked for is a tolerance that admits it, which is the subject of what a regression test can ask for.

What this costs somebody who has to ship

The uncomfortable version, since a reader who maintains a solver has to do something on Monday.

An algorithm that branches on definiteness — and there are many: a trust-region method choosing between a Newton step and a direction of negative curvature, a preconditioner deciding whether a Cholesky is available, an optimiser testing a Hessian — has a branch whose direction is a property of the build. Two machines then take different code paths, and after a different code path the answers are not close: they are the answers of two different algorithms.

That sounds worse than it is, and the reason is the one this essay’s middle section makes. The branch differs only when the true pivot is below the rounding, and in that regime both branches are defensible — a matrix that is definite by 10⁻²⁰ is a matrix a Cholesky cannot usefully factorise, and the negative-curvature branch is not the wrong thing to do with it. The divergence in the answers is real and the divergence in their quality usually is not.

Which yields the practical rule, and it is the same one negative curvature is the useful direction reaches from the other side: branch on a margin, never on a sign. Test the pivot against a multiple of the rounding — n·u·‖A‖ is the honest scale — rather than against zero. Both builds then take the same branch, because both agree about which side of a threshold that is orders of magnitude above their disagreement a number lies on.

The cost is one constant that has to be chosen and defended, which is the cost this site’s regularisation field spends four essays on. The alternative is a branch that no test can pin down.

At other settings

Two conforming builds, one matrix, and 24 different verdicts in 200 triesForty 4 × 4 Gram matrices at each of five conditionings, each formed exactly from a tall matrix whose last column is nearly the sum of its first two. Every one of them is asked *are you positive definite* by a Cholesky, once with the inner products accumulated with a single rounding per term and once with each product rounded before it is subtracted. The bar is how many of the forty came back with different verdicts: 5 at 10⁻⁸, 3 at 10⁻⁹, 6 at 10⁻¹⁰, 2 at 10⁻¹¹, 8 at 10⁻¹². The two shaded parts split those by which build was right, settled by the signs of the leading principal minors computed in BigInt rationals — 13 to the fused form and 11 to the rounded one. So neither is the correct implementation. The verdict is a property of the build, and the answer to which build is right is a property of the matrix.-12-11-10-9-80246810log₁₀ of how nearly dependent the columns areverdicts that disagreed, of 4053628which one is correctmatrices tested200verdicts disagreed24fused was right13unfused was right11lower part: thefused build was righta sign has no last digitso a verdict has nowhere to hide
Fig. 2 Four by four, where the last pivot has fewer inner products to disagree about.
Two conforming builds, one matrix, and 26 different verdicts in 200 triesForty 10 × 10 Gram matrices at each of five conditionings, each formed exactly from a tall matrix whose last column is nearly the sum of its first two. Every one of them is asked *are you positive definite* by a Cholesky, once with the inner products accumulated with a single rounding per term and once with each product rounded before it is subtracted. The bar is how many of the forty came back with different verdicts: 3 at 10⁻⁸, 6 at 10⁻⁹, 6 at 10⁻¹⁰, 6 at 10⁻¹¹, 5 at 10⁻¹². The two shaded parts split those by which build was right, settled by the signs of the leading principal minors computed in BigInt rationals — 8 to the fused form and 18 to the rounded one. So neither is the correct implementation. The verdict is a property of the build, and the answer to which build is right is a property of the matrix.-12-11-10-9-802468log₁₀ of how nearly dependent the columns areverdicts that disagreed, of 4036665which one is correctmatrices tested200verdicts disagreed26fused was right8unfused was right18lower part: thefused build was righta sign has no last digitso a verdict has nowhere to hide
Fig. 3 And ten by ten, where it has more — and the split between which build is right stays near even.
A determinant of one, computed as zero by the expression that is written downThe matrix [[x+2, x+1], [x+3, x+2]] at x = 2ᵏ has determinant exactly 1 at every k. The two products are x² + 4x + 4 and x² + 4x + 3, so they differ by one part in x², and above k = 26 they round to the same binary64 number. Relative error of two evaluations: ad − bc as written, which is exact until k = 27 and then returns zero at every larger size; and Kahan's form, w = bc, e = bc − w recovered by the fusion, (ad − w) − e, which is exact throughout. The difference between them is one rounding of one product — a relative 10⁻¹⁶ — and it is the whole of the answer, because the subtraction that follows removes everything else. Both forms are IEEE-754 conforming and a compiler may emit either from the same source.2022242628303210⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²k, where the entries are near 2ᵏrelative error of the determinantas writtenfused: exactproducts need 54 bitsone rounding, the whole answertrue determinant1naive, k = 300fused, k = 301first wrong at k27sizes returning 06both forms conformand the source does not say which
Fig. 4 The same rounding, on the smallest expression it can decide.
(x − 1)⁶ near x = 1, computed twice, negative in both and not in the same placesThe polynomial x⁶ − 6x⁵ + 15x⁴ − 20x³ + 15x² − 6x + 1 evaluated by Horner's rule at 401 points within 0.001 of x = 1, in units of 2.68·10⁻¹⁵ — the size of the rounding, since the true value here is under 10⁻¹⁸ and is a flat line on this axis. The upper trace rounds the multiply and the add separately; the lower one fuses them. A sixth power cannot be negative and the first is negative at 179 of the 401 points, the second at 196, and the two disagree about the sign at 98. Neither is more nearly right: both are the accumulated rounding of six cancelling terms, and they differ by at most 1.44·10⁻¹⁵ — one rounding of one product, six times. Any test on this value — a bisection, a definiteness check, a sign — gets a different answer from the two builds at a quarter of the points here.-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
Fig. 5 And on a polynomial, where neither build is nearer the truth than the other.
The summation condition number of six inner products this site already computesκ = Σ|xᵢ| ÷ |Σxᵢ| for six reductions taken from computations elsewhere on this site, at 64 terms: an orthogonality check between two columns of a Q, one component of a residual b − Ax on an ill-conditioned system, a nearly degenerate 2 × 2 determinant, a Hutchinson trace probe zᵀAz, a residual norm rᵀr, and a curvature pᵀAp. The range is 1 to 1.01·10¹⁷. The two at the safe end are sums of squares, where no term can cancel another; the three above 10¹⁰ are the orthogonality check, the residual and the determinant — which is to say every quantity a stopping test, an orthogonality test or an orientation test is written in. The reductions that are reproducible are the ones nobody makes a decision from.κ of the suma component of b − Ax1.01·10¹⁷ad − bc, near-degenerate3.6·10¹⁶qᵢᵀqⱼ, an orthogonality check7.39·10¹⁵pᵀAp, a curvature409zᵀAz, a trace probe41.8rᵀr, a residual norm1measured, not assumedhighest10¹⁷lowest1above 10¹⁰3terms64sums of squares are safeand nobody decides anything from one
Fig. 6 Why a barely definite matrix’s last pivot is at the wrong end of the census.
Least squares by QR and by the normal equations in binary32Relative error of the computed coefficients against epsilon, on log axes. The QR route is a flat line near the bottom; the normal-equations route climbs and then stops, at the epsilon where the cross-product matrix becomes exactly singular.10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹ε in Läuchli's matrix (smaller ε, larger κ)relative error in the coefficientsAᵀA exactly singularnormal equationsQRκ from 1.7·10⁸ to 17the cliff is at √u = 2.4·10⁻⁴
Fig. 7 The squaring that puts it there.
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. 8 The same decision with the threshold somewhere other than zero.
How often Cholesky still calls a 12×12 matrix positive definite, against λ_min/λ_max in units of the format's own roundoffThree curves, one per precision, of the share of 24 seeded matrices on which the factorisation succeeds. Measured in units of each format's unit roundoff the three lie almost on top of one another, with the edge — the smallest ratio at which every seed succeeds — at 1.0, 1.0, 1.8 times u. The absolute thresholds are 6·10⁻⁸, 2.3·10⁻¹⁰ and 2·10⁻¹⁶: nine orders of magnitude apart, and the same number in the format's own units.10⁻³10⁻²10⁻¹110¹10²10³10⁴00.250.50.751λ_min / λ_max, in units of the format's own ushare of seeds that succeed24 bits32 bits53 bitsone threshold, three formats24-bit edge, absolute6·10⁻⁸32-bit edge, absolute2.3·10⁻¹⁰53-bit edge, absolute2·10⁻¹⁶in units of u, at 53 bits1.8a yes-or-no question with a precision in itand a coin flip three decades below the edge
Fig. 9 Where a Cholesky stops, in the field that made the argument first.
The growth factor of a 12×12 elimination, against the condition number of the matrixThree curves against κ. Cholesky's growth factor is exactly 1 at every condition number drawn — the elimination never produces an entry larger than the matrix already had. Partial pivoting on the same matrices reaches the same growth, and reaches it by making up to 9 row interchanges, each of which destroys the symmetry that was the reason to use a symmetric factorisation. The bound the general theory allows, 2^11 = 2048, is drawn above them both.10¹10³10⁵10⁷10⁹10¹¹110¹10²10³10⁴condition number of the matrixgrowth factorbound 2^11partial pivotingCholeskyno pivot to gain fromCholesky growth, every κ1Cholesky interchanges0partial pivoting, worst9the bound, 2^112048both eliminations reach the same growthand only one of them had to swap to get there
Fig. 10 And what it does before it gets there.
The exact solution of a 13×13 Hilbert system beside the computed oneTwo columns of numbers: the exact answer, which is the integers one to thirteen, and the answer double-precision elimination returns, with the number of correct digits beside each.H13 x = b, b formed exactly so that x = (1, 2, …, 13)12345678910111213exact1.00002.00003.00133.97865.18864.993010.46720.048921.2657-2.575319.21478.906013.5113computed6.6 correct digits4.8 correct digits3.4 correct digits2.3 correct digits1.4 correct digit0.8 correct digitno correct digitsno correct digitsno correct digitsno correct digitsno correct digits0.6 correct digit1.4 correct digitbackward error2.2·10⁻¹⁷κ = 1.7·10¹⁸The algorithm solved a neighbouring problem perfectly. That problem's answer is this one.right-hand side built in BigInt rationalsthe truth is known
Fig. 11 The kind of standing this essay’s adjudicator provides.
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. 12 A computation that declines to answer, which is what a filtered verdict does.
The same 4096 numbers, added up 26 ways, 21 different answersEach dot is one reduction of one vector of 4096 binary64 numbers whose summation condition number is 1.01·10⁸: the vector is cut into p contiguous pieces, each piece summed left to right, and the pieces combined in a tree. The vertical position is the distance from the exactly rounded sum, relative. Filled dots are the static split a fixed thread count gives; open dots are unequal splits, which is what a work-stealing scheduler produces at the same p. 21 of the 26 runs returned distinct values, spanning 1.55·10⁻¹³ — 2.3·10⁷ ulps of the answer. The dashed line is the classical bound γ₍ₙ₋₁₎Σ|xᵢ|, which every one of them satisfies with 2.58·10⁴ to spare, and which is the same number for all of them: it contains n and Σ|xᵢ| and nothing about the order.110¹10⁻¹¹10⁻⁸10⁻⁵pieces the vector was divided intodistance from the exact sum, relativethe published boundκ · uone vector, one algorithmdistinct answers21runs26spread, in ulps2.3·10⁷κ of the sum10⁸bound ÷ worst error2.6·10⁴nobody chose pand no answer is the answer
Fig. 13 The field’s other cause of two answers.
Five summation policies, 400 permutations each, and how many answers each returnedThe same 1,024 numbers, summation condition number 10⁸, presented in 400 different orders. The bar is the number of distinct binary64 values the policy returned; the figure beside it is the worst relative error it made. A single accumulator returns 303 values and eight pieces return 72; Kahan's compensated loop is the most accurate of the three at 4.96·10⁻¹⁰ and still returns 119. Pre-rounded summation returns one value, at an error of 4.91·10⁻⁶ — four orders worse than the compensated loop and the same on every machine. Exact accumulation returns one value and the right one. The comparison is bitwise because the claim is bitwise: a tolerance here would pass everything.distinct answersone accumulator303eight pieces72compensated119pre-rounded1exact1worst error 1.17·10⁻⁸worst error 4.61·10⁻⁹worst error 4.96·10⁻¹⁰worst error 4.91·10⁻⁶worst error 0bitwise, or not at allpermutations400one accumulator303pre-rounded1its error4.9·10⁻⁶compensated error5·10⁻¹⁰accuracy and agreement are different propertiesand the accurate one is not the agreed one
Fig. 14 And the policy that removes it, which does not remove this one.
One system, one tolerance, and 11 different amounts of workConjugate gradients on a 200 × 200 symmetric positive definite matrix with κ = 10⁴, stopped when the relative residual falls below 10⁻¹⁰, run at each of 13 partition counts. The only thing that differs between the runs is how many pieces the two global inner products were summed in — which on a real machine is the number of workers. Iteration counts run from 674 to 690, 11 of them distinct, a spread of 2.4% of the work. Every run converged and every answer is right to the accuracy asked for; the forward errors span 6.31·10⁻¹¹ to 1.08·10⁻¹⁰, a factor of 1.72. What differs is not a digit at the end of an answer — it is the bill, and it is an integer.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
Fig. 15 A verdict that is an integer rather than a sign.
The numerical rank of one matrix is 10, 12, or somewhere betweenSingular values of one 60 × 14 matrix whose spectrum falls by a fixed factor per index, so there is no gap anywhere and the threshold decides the rank. Each curve is the same matrix, its Gram matrix formed with the inner products summed in a different number of pieces. The curves lie on top of each other for the first several values and separate below about 7.3·10⁻¹³, which is where forming AᵀA has put the rounding. The dashed line is the threshold, σ₁ · 10⁻¹⁴. Counting the values above it gives 12 at p = 1, 12 at p = 2, 12 at p = 4, 10 at p = 8, 10 at p = 16, 11 at p = 32, 11 at p = 60 — the rank of one matrix, as a function of how many workers were available. Not a digit of an answer: the number of columns a model built from this matrix would have.13579111310⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹index ksingular valueσ₁ · 10⁻¹⁴, the thresholdone matrix, three rankspartitionings7lowest rank10highest rank12threshold7.3·10⁻¹⁴σ₁7.3the curves separate in the noiseand the threshold is drawn through it
Fig. 16 And one that is a count of columns.
The smallest defect a regression test can catch, on a machine that will not repeat itselfA 150 × 150 system with κ = 10⁴, solved by conjugate gradients to 10⁻¹⁰ on six partitionings. The shaded band is how far the answer moves across those six when nothing at all is wrong: 3.24·10⁻¹², relative. Each pair of dots is one defect — a single entry of the matrix changed by a relative δ, an assembly slip or a stale coefficient — solved on the same six machines, showing its nearest and furthest run. A defect is catchable only when its nearest run clears the band, because a test whose threshold sits inside the band fails a correct build. Everything at or below 10⁻¹³ overlaps and is invisible; 10⁻¹² is the smallest that separates, at 4.61·10⁻¹². So the tolerance exists, it is bracketed between 3.24·10⁻¹² and 4.61·10⁻¹² — a factor of 1.42 — and it is neither zero nor the 10⁻⁸ that usually gets typed.-15-13-11-910⁻¹³10⁻¹¹10⁻⁹10⁻⁷log₁₀ of the defect, relative to the entry it sits inhow far the answer movedthe machine, on its owna tolerance with two sidesmachine band3.2·10⁻¹²smallest caught10⁻¹²its nearest run4.6·10⁻¹²defects hidden3window, factor1.4below the band nothing is visibleand above it everything is
Fig. 17 What can be asked for once the verdict is known to move.
The 14 eigenvalues of a saddle-point matrix with 10 unknowns and 4 constraints, inside their closed-form bracketsK = [[H, Aᵀ], [A, 0]] with κ(H) = 100 and κ(A) = 100. The shaded bands are the Rusten–Winther brackets, computed from four numbers — the extreme eigenvalues of H and the extreme singular values of A — before the matrix was assembled: negative eigenvalues in [-0.995, -10·10⁻⁵] and positive ones in [0.01, 1.62]. The marks are the computed spectrum. There are exactly 10 above zero and 4 below, which is Sylvester's law of inertia and not a property of this matrix: K is congruent to blkdiag(H, −AH⁻¹Aᵀ), both blocks are definite, and congruence preserves signs. The two groups are separated by a gap containing zero, at a ratio of 39.3 between the innermost positive eigenvalue and the innermost negative one.-1-0.582271-0.1645420.2531860.6709151.088641.506370eigenvalue4 negative10 positivecounted before it was formedpositive10negative4at zero0innermost ratio39the zero block is a theoremand so is the count either side of it
Fig. 18 An inertia that is known before anything runs, which is the opposite situation.
How much a perturbation of the right-hand side is amplified, κ = 10⁶The cumulative distribution of the amplification factor over two hundred random perturbation directions, with the condition number marked as the upper limit.110¹10²10³10⁴10⁵10⁶10⁷00.250.50.751amplification of the input perturbationfraction of directions at or belowκ = 10·10⁵worst found 7.6·10⁵6×6, 200 directionsmedian reaches 0.29 of κ
Fig. 19 The amplifier that puts the last pivot below the rounding.
Backward and forward error against the condition numberA log–log plot over twelve decades of condition number. The backward error is a flat line at ten to the minus sixteen; the forward error rises in proportion to the condition number.110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)relative errorforward errorbackward errorpredicted: κ · u8×8, 20 seeds per κ; dashed is the worstthe problem worsens, not the method
Fig. 20 The instrument that cannot adjudicate this, for the reason the previous essay gives.
What a perturbation of 10⁻⁶ moves, against the gap it is applied acrossThree curves against the gap between two eigenvalues, both axes logarithmic. Two are flat; the third rises steeply as the gap closes and then levels off.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
Fig. 21 A gap being closed, which is what the horizontal axis of the hero figure does.
Forward error of a 6×6 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, …, 6)12 bits3.316 bits0.8620 bits4.324 bits0.09330 bits7.6·10⁻⁴36 bits6.2·10⁻⁶43 bits2·10⁻⁷53 bits4.5·10⁻¹¹dashed: κ · unit roundoffone matrix, eight arithmeticsmeasured against a known answer
Fig. 22 The knob that moves where the disagreements happen and not whether they do.
Loss of orthogonality against condition number, in binary64A log–log plot of the norm of Q-transpose-Q minus the identity against condition number. Classical Gram–Schmidt rises steeply, modified Gram–Schmidt rises gently, and Householder is flat.110²10⁴10⁶10⁸10¹⁰10¹²10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)‖QᵀQ − I‖classicalmodifiedHouseholderκ²uκu8×8, eight seeds per κ, binary64all three reconstruct A
Fig. 23 Orthogonality as a number, computed from an inner product at the wrong end of the census.
QᵀQ from classical Gram–Schmidt and from Householder on the 8×8 Hilbert matrixTwo eight-by-eight tables of QᵀQ. The upper one has ones on the diagonal and entries as large as one off it; the lower one is the identity to three decimal places everywhere.A = H8 · κ = 1.5·10¹⁰ · both factorisations reconstruct A to 6·10⁻¹⁷the diagonal is 1 in both — every column is a unit vector either way1.000000000001.000000000001.000000000001.000000000001.00000.002-0.002000001.0000.125-0.13300000.0020.1251.000-1.0000000-0.002-0.133-1.0001.000classical Gram–Schmidt1.000000000001.000000000001.000000000001.000000000001.000000000001.000000000001.000000000001.000Householderclassical ‖QᵀQ − I‖1.4Householder ‖QᵀQ − I‖1.4·10⁻¹⁵largest off-diagonal 1 against 3.1·10⁻¹⁶length is not angle
Fig. 24 The same quantity entry by entry, where the cancellation that makes it fragile is visible.
Growth factor under partial pivoting to n = 40: the bound, the worst case, and realityGrowth factor against matrix size on a logarithmic vertical axis. The two-to-the-n bound rises as a straight line; Wilkinson's matrix sits exactly on it; random matrices stay near one.0816243240110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴matrix size ngrowth factor max|u| / max|a|the 2ⁿ⁻¹ boundworst of 30 randommedian randomWilkinson's matrix sits on the bound30 Gaussian matrices per sizeat n = 40: bound 5.5·10¹¹, worst 4.8
Fig. 25 A bound whose worst case is never attained, in the field that measured it first.
The backward error of a hierarchical solve, against the error of the representation it usedSolve A_H x = b exactly and the residual against the matrix that was wanted is b − Ax = (A_H − A)x, so ‖b − Ax‖ ⁄ ‖A‖‖x‖ cannot be anything but the representation's own error. Measured across six accuracies spanning ten decades, the two track at a slope of 0.994 and sit a constant 6.7× apart, which is the difference between a norm of a matrix and a norm of that matrix applied to one vector. The reading is the one this field is built on: the accuracy knob is not an accuracy, it is a backward error chosen in advance. Everywhere else on this site a backward error is something an algorithm produced and somebody then measured; here it is a line in the program, and its size is known before the solve starts.10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³‖A − A_H‖ ⁄ ‖A‖, 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
Fig. 26 A backward error chosen rather than incurred, which is the reading this field gives the machine’s own.
Two perturbation bounds and the error that was measured, on a 8×8 matrix spread over 8 decades of unitsThree curves against the size of an entrywise relative perturbation. The normwise bound κ_∞·ε is a valid bound and sits 4·10⁷ times above the componentwise one cond(A, x)·ε, which is also a bound and is nearly attained by the worst of forty random perturbations at each size.10⁻¹⁴10⁻¹³10⁻¹²10⁻¹¹10⁻¹⁰10⁻⁹10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1relative size of the entrywise perturbationrelative forward errorκ_∞ · εcond(A,x) · εmeasuredboth bounds holdκ_∞(A)1.9·10⁸cond(A, x)4.8ratio of the bounds4·10⁷both curves above the data are boundsand only one of them is a measurement
Fig. 27 An error measured entry by entry rather than in a norm, which moves the quantity to the dangerous end of the census.
Backward error with and without the structure, for Levinson and for elimination, n = 12Both solvers are backward stable on this family by the usual measure: the smallest perturbation of any kind that explains either answer stays below 3.1·10⁻¹⁷ at every ρ, and the two curves are indistinguishable. Insisting the perturbation be a symmetric Toeplitz matrix — the kind of object the problem was posed with — lifts both by orders, and at ρ = 0.999 Levinson's rises by a factor of 2.48·10⁶. Neither number is alarming in absolute terms; what is worth carrying is that the reassuring one is the one that is reported, and the two are not measuring the same thing.10⁻³10⁻²10⁻¹110⁻¹⁹10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹1 − ρbackward errorLevinson, Toeplitzelimination, ToeplitzLevinson, any kindelimination, any kindat ρ = 0.999Levinson, any perturbation3.1·10⁻¹⁷Levinson, Toeplitz only7.6·10⁻¹¹the ratio between them2.5·10⁶diagonal defect of the first0.97the number that is reportedand the number that is asked about
Fig. 28 A backward error measured against a structure rather than against any one matrix.

Where this goes next

2022242628303210⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²k, where the entries are near 2ᵏrelative error of the determinantas writtenfused: exactproducts need 54 bitsone rounding, the whole answertrue determinant1naive, k = 300fused, k = 301first wrong at k27sizes returning 06both forms conformand the source does not say which One multiply the compiler removed A determinant whose value is exactly 1, computed as exactly 0 by the expression that is written down, and exactly 1 by the same expression with the multiply and the add fused. Both forms conform to IEEE-754, both are legal compilations of the same source, and nothing in the program says which one you have. 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 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. 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.

What links here

Computed from the collection, not written here: the essays that point at this one.

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.

Cholesky factorisationCondition squaringDefiniteness testExact arithmeticExpression contractionFused multiply-addGram matrixPositive definite