A matrix that is definite on one machine

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

Twenty-six disagreements is not enough to read a split from either. Over 204 disagreements — two hundred matrices at each of seven gaps — it is 109 to 95, or 53% to 47%, which is a coin flip and is exactly what nine against seventeen is a plausible draw from.

Sweeping the size instead of the gap says the same thing on independent data.

The winner alternates with the size and the totals stay level. At n = 4, 5, 8 and 10 the disagreements number 24, 34, 31 and 26 of two hundred, and the fused build is right 13, 16, 19 and 8 times against the unfused build’s 11, 18, 12 and 18. Fused wins at four and eight, unfused at five and ten, and over all 115 disagreements the split is 56 to 59.

That is a second census on a different axis reaching the same verdict as the first, which is worth more than either alone. The gap sweep gave 109 to 95 over 204; the size sweep gives 56 to 59 over 115. Neither shows a build that is better, and the alternation is what a coin flip looks like when it is sliced four ways — a reader inspecting only n = 8 would conclude the fused form is right three times in five, and a reader inspecting only n = 10 would conclude the opposite with the same confidence.

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 half the time — which is the same conclusion the twenty-six suggested, arrived at without the misleading one-in-three.

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: 7 disagreements at a separation of 10⁻⁸, 2 at 10⁻⁹, 4 at 10⁻¹⁰, 10 at 10⁻¹¹ and 3 at 10⁻¹². Forty matrices at each of five points is not enough to read a shape from, and running two hundred at each of nine says the shape is not there.

It switches on between 10⁻⁷ and 10⁻⁸ and then stays flat at about 14% for seven decades. There is no band and no falling off on the far side. Past the threshold the two builds keep disagreeing at the same rate, because their disagreement is a rounding of fixed size and the pivot whose sign they are deciding stays below it however much more nearly dependent the columns become.

The measurement that explains the flatness is the third one in the table. Where the builds disagree, the last pivot is a tenth of the entry scale’s own rounding — the median of |pivot| ÷ (max|aᵢⱼ|·u) is between 0.09 and 0.19 at every gap. Both builds are deciding a sign from a quantity below their own resolution at every conditioning past the threshold, so nothing about the rate has any reason to move.

What that leaves of the practical reading is a sharper version of it. A test suite sampling at 10⁻⁶ sees perfect agreement and concludes there is nothing to see; the same suite anywhere below 10⁻⁸ sees a seventh of its matrices flip, and moving further down does not make it worse. The threshold is the thing to know, not a band.

Where the threshold is, and what sets it #

A threshold is more useful than a band because it can be computed. It sits where the true smallest pivot falls below the rounding of the entries, and both quantities are available: the pivot is what an exact Cholesky would return and the rounding is max|aᵢⱼ|·u. On the family here the columns are separated by a gap g, the Gram matrix squares it, and the smallest pivot is therefore of order g² times the entry scale — so the threshold is at g² ≈ u, which is g ≈ 1.5·10⁻⁸.

Measured, it is between 10⁻⁷ and 10⁻⁸. The prediction is a square root of the machine’s own precision, and it is the same √u this collection keeps arriving at from other directions — the normal equations’ cliff, the definiteness test’s margin, the weight at which a constraint stops being imposable. Here it is the separation at which a Gram matrix stops having a decidable sign.

That gives the practical rule the census was reaching for. A definiteness verdict is build-dependent exactly when the columns are separated by less than √u, which a caller can estimate without factorising anything: it is a statement about the data, computable from the data, and it says which of two regimes the answer is in before the answer exists.

The rule has a corollary that is worth more than the rule, because it applies to a decision nobody thinks of as numerical. A Gram matrix is formed on purpose, from columns somebody chose, and the √u threshold says how close together those columns may be before the definiteness of their Gram matrix stops being a fact about them. Eight decimal digits of separation is not a demanding requirement in most settings and is an impossible one in a few — a design matrix with two nearly collinear predictors, a kernel matrix at points that nearly coincide, a covariance estimate from samples that are nearly repeats.

In those settings the right response is not a better factorisation. It is to stop asking the question: a matrix whose definiteness is undecidable at the working precision is a matrix whose columns should be merged, dropped or regularised, and every one of those is a decision about the model rather than about the arithmetic. The census is a measurement of when the arithmetic stops being able to answer, and the answer to what then is in a different field.

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:

• a pivot’s sign in a symmetric indefinite factorisation, which decides between a 1 × 1 and a 2 × 2 block — when symmetry is not enough;
• a residual’s comparison with a tolerance, which decides whether to stop — a stopping test is a race;
• a singular value’s comparison with a threshold, which decides a rank — a rank that depends on the thread count;
• a polynomial’s sign near a root, which decides which side a bisection takes — a square that evaluates negative.

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.