A factorisation with nothing to pivot for

Everything this site has said about elimination so far has had a pivot in it. The interchange is what makes it stable, the growth factor is what says how much it could still cost, and Wilkinson’s matrix is what shows the bound being attained under the rule that was supposed to prevent it.

For a symmetric positive definite matrix all of that goes away, and it goes away for a reason that takes two lines.

The two lines #

One. In a positive definite matrix the largest entry in absolute value is on the diagonal. If some |a_ij| > max_k a_kk with i ≠ j, take the vector with ±1 in positions i and j and zeros elsewhere; the quadratic form it produces is a_ii + a_jj ∓ 2|a_ij| < 0, contradicting definiteness.

Two. Every Schur complement of a positive definite matrix is positive definite, and its diagonal entries are no larger than the ones it came from: eliminating with pivot a_kk replaces a_ii by a_ii − a_ik²/a_kk, and the subtracted quantity is non-negative.

Put them together. Every entry anywhere in the elimination is bounded by the largest diagonal entry of the current Schur complement, which is bounded by the largest diagonal entry of A, which is the largest entry of A.

So the growth factor is at most 1. And it is at least 1, because the first step has not changed anything yet. It is exactly 1.

What the second curve is for #

Partial pivoting on these matrices also has a growth factor of 1. That is not a coincidence and it is not a reason to shrug: it is the same theorem, since partial pivoting on an SPD matrix is doing the same eliminations in a different order.

What separates the two is underneath the curves and is in the badge: the interchanges. Cholesky makes zero. Partial pivoting makes up to nineteen at n = 24, and every one of them destroys the symmetry that was the reason for using a symmetric factorisation in the first place.

That is the honest cost of the general routine on this class of matrix. Half the storage, half the arithmetic, and no permutation to record — all three lost to a search that finds nothing, at a growth factor identical to the one the symmetric routine gets for free.

The site’s own same arithmetic at a different price is the same shape of argument reached through memory traffic rather than through pivots.

No pivot to gain from, not “no pivot needed” #

The distinction is worth pressing on, because the two sentences suggest different things about risk.

Cholesky does not need pivoting reads as a claim about typical behaviour, of the kind that has an unstated exception somewhere. What is true is stronger and has no exception: there is no interchange that would improve anything, because the quantity pivoting exists to control is already at its minimum possible value. A pivot search on an SPD matrix is a search whose best possible outcome is the answer it already has.

The one thing a symmetric permutation does buy on an SPD matrix is sparsity, which is a different objective entirely and is the order decides the memory. A sparse Cholesky permutes heavily, chooses the permutation from the graph before looking at a single value, and gets the same growth factor of 1 whatever it chooses. That the two objectives do not interfere is exactly why sparse SPD solvers are the easiest kind to write, and it is the property that structure and stability stop being separable shows disappearing the moment the matrix stops being definite.

The failure, which is the test #

Cholesky fails in exactly one way. At step k it asks for the square root of the current a_kk, and if that is not positive there is no real square root and the algorithm stops.

And it fails if and only if the matrix is not positive definite. That is not a coincidence either: the pivots of a Cholesky factorisation are the diagonal of D in LDLᵀ, their product over the first k is the k-th leading principal minor, and Sylvester’s criterion says a symmetric matrix is positive definite exactly when all of those are positive.

So the factorisation is the standard test. It costs n³/3 operations against 4n³/3 for an eigenvalue decomposition, it answers the same yes-or-no question, and it is what every library does when asked isposdef.

In floating point the “if and only if” moves #

assertTheDefinitenessTestMovesWithThePrecision builds QΛQᵀ with Λ = (1, …, 1, δ), so δ is λ_min/λ_max and the true answer is an input rather than a measurement. Then it asks, at each of three precisions, how small δ can be before the test starts saying no.

The edge — the smallest ratio at which every seed succeeds — is at 1.0·u, 1.0·u and 1.8·u at 24, 32 and 53 bits. The absolute thresholds are 6.0·10⁻⁸, 2.3·10⁻¹⁰ and 2.0·10⁻¹⁶: nine orders of magnitude apart, and the same number in each format’s own units.

That is the site’s precision knob applied to a question whose output is a yes or a no rather than a number of digits, and it is the same shape as the mixed-precision threshold at κu ≈ 1 which buying the accuracy back measured.

Why a share and not a verdict #

The figure plots a share of seeds rather than a single edge per precision, and the reason is a finding rather than a presentational choice.

Far below the edge the outcome is not a function of δ at all. The matrix is assembled and stored before the factorisation begins, and assembling it rounds — so the smallest eigenvalue the routine is actually handed is δ plus a rounding error of size about u, whose sign depends on the seed. Three decades below u the share settles at about a half and stays there: a coin flip, at every δ.

A bisection on one seed would have returned a number to fourteen digits, and that number would have been a property of that seed’s rounding. The site has recorded this before — near a threshold the outcome is genuinely erratic, from the mixed-precision phase, and a refusal that waits for a rare event is not a refusal from the expansion. The honest measurement is the share, which is smooth.

The jitter, and what it actually costs #

There is a piece of folklore attached to this threshold that is worth measuring rather than repeating.

Kernel and covariance matrices are positive semi-definite by construction and positive definite only if the data says so. In practice they arrive with a smallest eigenvalue near zero — duplicated observations, a kernel width larger than the spread of the points, a covariance estimated from fewer samples than dimensions — and chol throws. The standard response is to add jitter·I and try again, doubling the jitter until it works.

What the figure above says is that the loop terminates when the jitter reaches about u·λ_max, and what it does not say is what has been bought. Adding εI to a matrix whose smallest eigenvalue is zero does not recover the missing direction; it replaces it with a direction of size ε, which is a statement about arithmetic and not about the data. The factorisation then succeeds, the residual ‖A + εI − LLᵀ‖ is at rounding, and every quantity a solver can see says the computation went perfectly.

The honest summary is that jitter converts a factorisation that fails into a factorisation of a different matrix, and the difference is exactly the thing the failure was reporting. That is the same move as a modified Cholesky, performed by the caller rather than by the library, usually without the record of how far the matrix was moved.

What Cholesky’s backward error actually promises #

Growth factor 1 is a statement about the intermediate quantities. The statement about the answer is Wilkinson’s, and it has a condition in it that is worth reading.

The computed factors satisfy L̂L̂ᵀ = A + ΔA with ‖ΔA‖ ≤ c(n)·u·‖A‖ — backward stable, with a modest polynomial in n and no growth factor anywhere, which is the payoff of the two lines at the top of this essay. And the theorem’s hypothesis is that the algorithm runs to completion: if c(n)·u·κ₂(A) < 1 it does, and past that it may not.

So the failure mode and the stability statement are the same boundary seen from two sides. Below it, Cholesky is as stable as an elimination gets; above it, it does not return an answer at all. There is no region in which it returns a bad one.

That is unusual and it is worth naming, because almost nothing else on this site behaves that way. The road that squares the problem degrades smoothly and silently; the swap that is not optional returns an answer of the right shape that is wrong; the normal equations lose digits without comment.

The square root, and the variant without it #

LDLᵀ with D diagonal and L unit lower triangular is the same factorisation with the square roots pulled out: L_chol = L·D^{1/2}. Everything above applies to it unchanged — same growth, same pivots, same failure condition, with the failure now a non-positive d_kk rather than a square root of one.

Two reasons it gets used. Square roots are slow, historically much more so than they are now, and there are n of them. And it extends: LDLᵀ makes sense for an indefinite symmetric matrix, where D is allowed negative entries, and Cholesky does not. That extension is the whole of the next essay, and the D that comes out of it is not diagonal.

There is a third reason, and it is the one that survives modern hardware. In a Kalman filter or a square-root information filter the quantity being propagated is a factorisation rather than a matrix, and a step that requires the covariance to be re-formed and re-factorised loses the guarantee that it stays definite at all. LDLᵀ can be updated in place, and its definiteness is then a property of the signs in D rather than of an arithmetic that might have drifted.

What a library does when the test says no #

Two things, and they answer different questions.

A modified Cholesky perturbs the matrix until it is definite and factorises that — replacing a non-positive pivot by something small and positive, or adding a multiple of the identity — and returns the factors of a different matrix, along with how far away it is. That is what an optimisation code does inside a Newton step, where the Hessian may be indefinite and a descent direction is wanted rather than the exact Newton one.

A symmetric indefinite factorisation does not perturb anything and does not need definiteness at all, at the cost of a pivot rule that can take two variables at once. That is when symmetry is not enough, and it is the essay after this one.

The distinction matters more than it looks. The first changes the problem and says so; the second solves the problem given. A code that silently did the first would be returning an answer to a question its caller did not ask, and the reason chol throws rather than repairing is that there is no repair that is right for every caller.

What “positive definite” is, computationally #

A last observation about the word, because it is used in two ways that come apart exactly at the threshold this essay measures.

Mathematically, positive definite means xᵀAx > 0 for every nonzero x, which is a statement about infinitely many vectors and is equivalent to every eigenvalue being positive.

Computationally, it means the matrix a routine was handed factorises. Those agree above the edge and not below it, and below it the second is not even a property of the matrix — it is a property of the matrix, the format, and the seed that generated it.

The consequence for a caller is small and specific: isposdef is not a predicate. It is a test whose answer is a fact about the arithmetic as much as about the argument, and two libraries can disagree about the same matrix without either being wrong.

That is unusual and it is worth being clear that it is the matrix rather than the test that is at fault. A matrix whose smallest eigenvalue is 10⁻¹⁸ of its largest is, for every computational purpose, indistinguishable from one whose smallest is −10⁻¹⁸, and no test that reads the stored entries can distinguish them because the stored entries are the same.

What is worth carrying #

Growth factor exactly one, at every size and every conditioning. Not bounded — equal, with a two-line proof and nine measurements. That is the strongest stability statement anywhere on this site, and it belongs to the algorithm with no pivot rule in it.

The generality of partial pivoting is not free on this class, and what it costs is not accuracy. It is the symmetry, and with it half the work, half the storage, and the permutation nobody has to record.

The failure is the test, and it is the cheapest definiteness test there is.

And in floating point the test’s boundary is a small multiple of the unit roundoff, at every precision — so “is this matrix positive definite” is a question whose answer depends on the format it is asked in, and three decades below the boundary the answer is a coin flip.