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ᵢⱼ| > maxₖ aₖₖ with i ≠ j, take the vector with ±1 in positions i and j and zeros elsewhere; the quadratic form it produces is aᵢᵢ + aⱼⱼ ∓ 2|aᵢⱼ| < 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ₖₖ replaces aᵢᵢ by aᵢᵢ − aᵢₖ²/aₖₖ, 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 has a growth factor of 1 as well, or a shade under it. 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 eliminations of the same Schur complements 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.

Across n = 5, 8, 12, 18 and 24 the interchange count runs 2, 4, 9, 13 and 19 — close to three-quarters of n, and every one of them a swap that changes nothing except the symmetry. The bound the general theory allows over the same sizes is 16, 128, 2,048, 1.3·10⁵ and 8.4·10⁶, and Cholesky’s growth is exactly 1 at all five with the argument for it never mentioning n.

Partial pivoting’s growth is not exactly 1, and the sweep is where that shows. It reads 1, 1, 0.9622, 1 and 0.9979 across the five sizes — never above one, and at two of them below it. An interchange on a symmetric matrix destroys the symmetry the two-line argument depends on, so the equality stops being an equality; what survives is the inequality, and on a matrix this well-behaved the largest entry can shrink slightly rather than staying put. That is a smaller claim than the one the section began with and it is the one the figures support.

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ₖₖ, 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·λₘₐₓ, 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 factorisation — and there is a narrow one in which it returns a bad verdict, which the next section measures.

The “only if” fails too, and that half is quieter #

Every measurement in the definiteness section approaches the boundary from above: δ is positive, the matrix is definite, and the question is how small δ can be before the routine starts saying no. The other direction is the one isposdef is actually used for. Taking δ negative, so the matrix is indefinite by construction, and asking how often Cholesky accepts it anyway:

    δ ⁄ u      accepted    worst backward error    smallest stored λmin ⁄ λmax
     +4        24 of 24        1.6·10⁻¹⁶                  +2.5·10⁻¹⁶
     +1        23 of 24        1.5·10⁻¹⁶                  −5.0·10⁻¹⁷
      0        11 of 24        1.5·10⁻¹⁶                  −2.7·10⁻¹⁷
    −0.5       10 of 24        1.5·10⁻¹⁶                  −2.7·10⁻¹⁷
    −0.9        0 of 24             —                          —
     −2         0 of 24             —                          —

Ten of twenty-four seeds at δ = −0.5u are certified as positive definite, and they are not. The third column says the matrix as stored has a genuinely negative smallest eigenvalue in some of the accepted cases, so this is not the construction’s spectrum disagreeing with the assembled matrix — it is the factorisation accepting a matrix an eigensolver rejects.

The band is narrow. Two roundoffs below zero nothing is accepted, so the “if and only if” holds outside a window about one unit roundoff wide, which is the same window and the same coin flip the share-not-a-verdict section describes from the other side. What is different is the consequence. A failure is loud: the routine throws and a caller decides what to do. A false certificate is silent, and the fourth column says why nothing catches it — the factorisation returned in those cases is backward stable to 1.5·10⁻¹⁶, exactly like every other, so ‖A − LL̂ᵀ‖ reports success, a solve with it has a small residual, and every quantity the caller can compute agrees that the computation went well.

So the accurate form of the essay’s claim is two sentences rather than one. Cholesky never returns a bad factorisation. And used as a definiteness test it is right outside a band of width one unit roundoff and a coin flip inside it, in both directions — which matters because the test is the reason the factorisation is being run in a good deal of the code that runs it.

That has a practical shape, and it is the same one the jitter section arrives at from the other direction. A caller who needs the verdict rather than the factors — is this covariance estimate usable, is this Hessian a minimum, is this kernel matrix a valid one — is asking a question whose answer near the boundary is not determined by the matrix. The response is not a better factorisation; it is to decide what the answer should be when the smallest eigenvalue is within a roundoff of zero, say so, and enforce it. Accepting a matrix whose λmin is above +c·u·λmax for a stated c is a decision a reader can audit. Accepting whatever Cholesky happened to accept is a decision made by the seed.

And that is the recurring shape rather than a fact about Cholesky. A yes-or-no question answered by a computation with a rounding in it has a band around its boundary where the answer belongs to the arithmetic, and every essay on this site that has met one has ended in the same place: the band is not removable, so the choice is between naming it and pretending it is not there.

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_\mathrm{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ₖₖ 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.

Where the same pivots come out negative #

A positive definite matrix has positive pivots, and that is what makes the symmetric factorisation need no interchange. In floating point it is a statement about the stored matrix rather than about the mathematical one, and the Hilbert matrix’s fourteenth pivot is negative.