Numerical Linear Algebra
where the algebra stops being the arithmetic
Iterating, instead of factorising

The division that cannot be done

Conjugate gradients divides by pᵀAp at every step, and on a matrix that is not positive definite that number can be zero or negative. The guard against it has been here from the first essay and described it as a failure. In the method that made conjugate gradients famous it is the single most valuable object the iteration can produce, and it costs six matrix–vector products.
21 min read 7 figures Where the algebra stops being the arithmeticWhose fault is itMeasured, not assumedTwo routes to a number

Worth reading first: The zero that means it is finished · The rate the condition number predicts · A factorisation with nothing to pivot for.

Two essays about a divisor going to zero, and both of them about what the zero means. This one is about a division that is not the point at all.

Conjugate gradients computes a step length as α = rᵀr / pᵀAp. The denominator is the curvature of the quadratic ½xᵀAx − bᵀx along the search direction p, and on a positive definite matrix it is positive by definition. On anything else it can be zero or negative, at which point the method has no step to take: the quadratic it is minimising has no minimum along that direction, and no length of step is better than any other because arbitrarily long is arbitrarily better.

This site’s own conjugate gradient routine has guarded against that since it was first written, with a comment saying the method “has no step to take” and a flag named breakdown. That is the correct thing for a linear solver to do and it is the whole of what this collection had ever said about it.

It is the wrong thing for the method that made conjugate gradients famous outside linear algebra.

The curvature pᵀAp/pᵀp along conjugate gradients on a matrix with one eigenvalue at −10^-1The matrix is 40×40, positive definite apart from a single eigenvalue at -0.1. Conjugate gradients divides by pᵀAp at every step, and the quantity runs 5.43, 2.72, 1.71, 0.687 before turning negative at step 6, where it is -0.02656. That direction is a proof: one matrix–vector product from outside the iteration confirms it, and it recovers 26.6 per cent of the eigenvalue. MINRES on the same system never forms this quantity, meets nothing, and returns the solution in 37 steps.012345678-7-5-3-11357conjugate gradient steppᵀAp ⁄ pᵀpλₘᵢₙ = -0.1positive: a step existsnegative: a certificate existsone matrix, two questionsstep it turns at6quotient there-0.027share of λₘᵢₙ recovered0.27λₘᵢₙ, by construction-0.1MINRES steps on the same system37the division that cannot be doneis the answer to a different question
Fig. 1 The quantity the method divides by, along the iteration. It changes sign, which is why this axis is linear and not logarithmic.

What the direction is

In a trust-region Newton method the system being solved is the Newton system: A is a Hessian, b is minus a gradient, and the quadratic is a local model of a function somebody is minimising. A direction p with pᵀAp ≤ 0 is a direction along which that model decreases without bound. The quantity being tested is a Rayleigh quotient, and the gap decides the eigenvector is the account of how much of one a matrix will let a computation resolve.

That is not a failure to report. It is a certificate, and it has three properties that make it worth more than the answer the solve was going to produce.

It proves something. The current point is not a minimum. Not “appears not to be”, not “the solver had trouble” — a vector d with dᵀAd < 0 is a proof that A is not positive semidefinite, and a proof that survives being handed to somebody who did not run the iteration — the same standard a bound that is proved holds a bound to.

It is checkable in one product. Anybody can take d, form Ad, take the inner product, and see the sign. The certificate does not require trusting the iteration that produced it, which is the difference between a proof and a report, and the property proving the answer is in the box is built around.

And it costs what has already been spent. On the matrix in the figure above it arrives after six matrix–vector products on a forty-dimensional problem. Establishing the same fact by attempting a Cholesky factorisation costs a third of n³ operations and requires the matrix to exist as entries.

How long it takes to find

The direction is not free — the iteration has to build a Krylov space large enough to contain a direction that sees the negative eigenvalue. How large depends on how negative it is, and the answer runs in the reassuring direction.

At λ = 3 the certificate arrives after three products. At λ = 10⁻³ it takes nine. The trend is monotone and it is the trend a Krylov method always has: large eigenvalues in magnitude are found first, small ones last.

The curvature pᵀAp/pᵀp along conjugate gradients on a matrix with one eigenvalue at −10^-3The matrix is 40×40, positive definite apart from a single eigenvalue at -0.001. Conjugate gradients divides by pᵀAp at every step, and the quantity runs 5.44, 2.76, 1.83, 0.884 before turning negative at step 9, where it is -2.811·10⁻⁴. That direction is a proof: one matrix–vector product from outside the iteration confirms it, and it recovers 28.1 per cent of the eigenvalue. MINRES on the same system never forms this quantity, meets nothing, and returns the solution in 39 steps.012345678910-7-5-3-11357conjugate gradient steppᵀAp ⁄ pᵀpλₘᵢₙ = -0.001positive: a step existsnegative: a certificate existsone matrix, two questionsstep it turns at9quotient there-2.8·10⁻⁴share of λₘᵢₙ recovered0.28λₘᵢₙ, by construction-0.001MINRES steps on the same system39the division that cannot be doneis the answer to a different question
Fig. 2 λmin = −10⁻³ on a forty-square Hessian. The quotient turns negative at step 9 of the 120 allowed, and the value it turns at is −2.81·10⁻⁴ — 28.1% of λmin.
The curvature pᵀAp/pᵀp along conjugate gradients on a matrix with one eigenvalue at −10^-2The matrix is 40×40, positive definite apart from a single eigenvalue at -0.01. Conjugate gradients divides by pᵀAp at every step, and the quantity runs 5.44, 2.75, 1.81, 0.865 before turning negative at step 8, where it is -0.004188. That direction is a proof: one matrix–vector product from outside the iteration confirms it, and it recovers 41.9 per cent of the eigenvalue. MINRES on the same system never forms this quantity, meets nothing, and returns the solution in 38 steps.0123456789-7-5-3-11357conjugate gradient steppᵀAp ⁄ pᵀpλₘᵢₙ = -0.01positive: a step existsnegative: a certificate existsone matrix, two questionsstep it turns at8quotient there-0.0042share of λₘᵢₙ recovered0.42λₘᵢₙ, by construction-0.01MINRES steps on the same system38the division that cannot be doneis the answer to a different question
Fig. 3 λmin = −10⁻²: the turn comes at step 8, at −4.19·10⁻³, which is 41.9% of λmin.

Two things are in those numbers and only the first is the section’s claim. The turn step falls 9, 8, 4, 2 as λmin goes from −10⁻³ to −10 — steeply, and in the reassuring direction. But the quotient at the turn is 28.1%, 41.9%, 28.7% and 35.1% of λmin across the same range: the certificate never arrives at the eigenvalue, only at a direction that has seen a third of it.

The curvature pᵀAp/pᵀp along conjugate gradients on a matrix with one eigenvalue at −10^0The matrix is 40×40, positive definite apart from a single eigenvalue at -1. Conjugate gradients divides by pᵀAp at every step, and the quantity runs 5.4, 2.36, 0.566, -0.287 before turning negative at step 4, where it is -0.2868. That direction is a proof: one matrix–vector product from outside the iteration confirms it, and it recovers 28.7 per cent of the eigenvalue. MINRES on the same system never forms this quantity, meets nothing, and returns the solution in 36 steps.012345678-7-5-3-11357conjugate gradient steppᵀAp ⁄ pᵀpλₘᵢₙ = -1positive: a step existsnegative: a certificate existsone matrix, two questionsstep it turns at4quotient there-0.29share of λₘᵢₙ recovered0.29λₘᵢₙ, by construction-1MINRES steps on the same system36the division that cannot be doneis the answer to a different question
Fig. 4 λmin = −1, where the turn is at step 4 and the quotient there is −0.287.
The curvature pᵀAp/pᵀp along conjugate gradients on a matrix with one eigenvalue at −10^1The matrix is 40×40, positive definite apart from a single eigenvalue at -10. Conjugate gradients divides by pᵀAp at every step, and the quantity runs 5.09, -3.51 before turning negative at step 2, where it is -3.507. That direction is a proof: one matrix–vector product from outside the iteration confirms it, and it recovers 35.1 per cent of the eigenvalue. MINRES on the same system never forms this quantity, meets nothing, and returns the solution in 37 steps.012345678-7-5-3-11357conjugate gradient steppᵀAp ⁄ pᵀppositive: a step existsnegative: a certificate existsone matrix, two questionsstep it turns at2quotient there-3.5share of λₘᵢₙ recovered0.35λₘᵢₙ, by construction-10MINRES steps on the same system37the division that cannot be doneis the answer to a different question
Fig. 5 And λmin = −10: the turn at step 2, at −3.51.

That matters for what the certificate may be used for. It is a proof that the Hessian is indefinite and it is not an estimate of how indefinite — a third of the answer, with the fraction itself moving between 28% and 42% for no reason the caller can see. A method that turned the quotient into a step length proportional to the curvature would be using a number that is wrong by a factor of three and inconsistently so.

And the cost of the detection is flat, which the trend above hides. MINRES takes 39, 38, 36 and 37 steps to reach the rounding level across the same four Hessians — a spread of 8% while the turn step moves by a factor of 4.5. So the expensive part of the computation does not know or care how negative the eigenvalue is; only the early exit does. A caller who budgeted from the turn step would be budgeting from the one quantity that varies.

So the indefiniteness that hides is the indefiniteness that matters least. A Hessian with an eigenvalue at −3 is a point the optimiser should leave immediately and the iteration says so at once; a Hessian with an eigenvalue at −10⁻³ is a point that is very nearly a minimum, and the extra six products spent finding that out are six products spent on a distinction that changes little.

That is not a general law about detection — it is a fact about the interaction between which eigenvalues a Krylov space finds and which eigenvalues matter here, and the two happen to be the same ones. It is worth naming precisely because the reverse arrangement is what one usually gets. On the rank decision this collection has already drawn, on the condition estimate, on the breakdown of the previous essay, the quantity that is hardest to detect is the one whose absence does the most damage.

And the direction is not the eigenvector

The certificate proves indefiniteness and it does not deliver λmin. Its Rayleigh quotient recovers between eleven and forty-two per cent of the negative eigenvalue across the eight spectra above, with no obvious pattern in which.

That is the honest reading and it is enough. A trust-region method does not need λmin; it needs a direction of descent that the model does not bound, and any direction with negative curvature is one. Getting λmin as well would cost the full eigenvalue problem, which is the thing the whole approach exists to avoid.

What the site’s own guard was doing

The guard in this collection’s conjugate gradient routine is worth quoting against itself, because it is a small, correct, and complete example of the mistake this essay is about.

It tests whether pᵀAp is positive, and if it is not it sets a flag and returns. The flag is called breakdown; the comment beside it says the method has no step to take. Every word of that is true. What it does not do is keep the vector.

The vector is right there — it is p, the direction the loop was about to divide by — and by the time the flag has been set and the function has returned, it is gone. Recovering it costs another run of the whole iteration. Keeping it costs one array reference.

That asymmetry is why this essay exists rather than a note in a source comment. The routine was written to solve linear systems, it does that correctly, and the information it discards is worthless for that purpose and is the entire answer to a different question. Nothing about the code is wrong. What is missing is a return value, and nothing in the way the routine is tested could have said so: every assertion on it is about the solution it produces on matrices where it produces one.

The same matrix, asked a different question

Put the certificate beside what a different method does with the identical system and the point of this essay is a comparison rather than a claim.

MINRES solves the same Ax = b on the same indefinite A. It builds the same Krylov space by the same products. It never forms pᵀAp at all — it minimises the residual over the space rather than the A-norm of the error, and nothing anywhere in it requires A to be definite. So it never meets the event, and on the matrix at the top of this essay it returns the solution in thirty-seven steps at a relative residual of 2.4·10⁻¹⁴ — a number worth reading with a small residual is not a small error in mind, since a residual is what the method minimises and not what it was asked for.

Conjugate gradients on the same matrix stops after six products with a residual nowhere near zero and a vector.

Both are right. They were asked different questions. One was asked for the solution of a linear system, which exists — A is nonsingular — and it produced it. The other was asked to minimise a quadratic, which has no minimum, and it produced the proof that there is none. The breakdown belongs to the question rather than to the matrix.

That reading has a consequence worth stating for anybody choosing between them. If the system is genuinely a linear system and A merely happens to be indefinite, conjugate gradients is the wrong method and MINRES is the right one, and the “breakdown” is the method correctly declining a problem it does not solve. If the system is a Newton system, MINRES is the wrong method: it will happily return the stationary point of a model that has no minimum, which is a saddle, and the caller will step to it.

What following it is worth

The direction on its own is unbounded, so a trust-region method bounds it: minimise the model inside a ball of radius Δ, and follow the negative-curvature direction to the boundary of that ball. Steihaug’s method is conjugate gradients with two extra tests — one for leaving the ball and one for the curvature — and it stops at whichever happens first.

The question is what that cheap step gets compared with solving the trust-region subproblem exactly, which needs a full symmetric eigendecomposition and a secular equation solved in the shift.

Steihaug's model decrease as a share of the exact trust-region answer, against the radiusThe quadratic model is ½xᵀAx − bᵀx on the same 40×40 matrix, positive definite apart from one eigenvalue at -0.1. At Δ = 0.25 the truncated conjugate gradient gets 89.2 per cent of what the exact subproblem solution gets, for 2 matrix–vector products against a full symmetric eigendecomposition and a secular equation. The share falls monotonically to 36.9 per cent at Δ = 16. The open marks are the radii at which the step stopped because the curvature turned rather than because it reached the boundary; the crossover is at Δ = 4.110¹00.30.60.91.2trust-region radius Δshare of the exact model decreasethe exact subproblem7 radii, n = 40share at the smallest radius0.89share at the largest0.37products, at most6radii stopped by the curvature3a few products against an eigendecompositionand most of the decrease
Fig. 6 Steihaug’s model decrease as a share of the exact one, against the radius. The open marks are radii at which the curvature stopped the step rather than the boundary.

At Δ = 0.25 the truncated iteration gets 89.2 per cent of the exact decrease for two matrix–vector products. The share falls monotonically as the region grows — 84.9, 76.1, 53.6, 48.1, 42.9 — and at Δ = 16 it is 36.9 per cent. Every stop costs at most six products.

The shape of that curve is the argument for the method and it is worth reading carefully, because both ends matter.

At a small radius it is nearly exact. That is where a trust-region method spends its time when it is in trouble: the radius shrinks precisely when the model has been a poor predictor, and it is exactly then that the cheap step loses almost nothing.

At a large radius it is well short. When the region is generous the truncated path is a conjugate-gradient path with a stopping rule bolted on, and it does not curve the way the exact solution curves. A method that spent most of its time at large radii would be leaving half the available decrease on the table.

And the reason for stopping changes across the sweep. Below Δ = 4 the step is stopped by the ball; above it, by the curvature. That crossover moves with the negative eigenvalue: a very negative one is met immediately and the curvature stops every step, while a tiny one is met after the boundary has already been reached.

What it does not certify

Two things the direction does not prove, and both are worth being explicit about because the word certificate invites over-reading.

It does not say how many negative eigenvalues there are. One direction with negative curvature proves at least one, and that is all. A matrix with three negative eigenvalues and a matrix with one produce the same kind of object, and distinguishing them needs an inertia count — which is a factorisation, and is the thing the whole approach was avoiding.

The block diagonal of a Bunch–Kaufman factorisation of an 10×10 saddle-point matrix at τ = 10⁻⁶The matrix D from PAPᵀ = LDLᵀ, drawn as a matrix. 5 of its blocks are 2×2 — the shaded pairs — and the rest are single entries. A 2×2 block is taken exactly where no single diagonal entry is large enough to divide by safely, which on a matrix whose lower-right block is zero by construction is most of it. The badge carries the residual of all three pivot rules.D from PAPᵀ = LDLᵀ — the shaded pairs are 2×2 pivots10⁻⁶-0.8········-0.8···········2·10⁻⁶0.9········0.9···········2·10⁻⁵1········1············-0.7········-0.73·10⁻⁶··········2·10⁻⁶0.7········0.7·three rules, one matrix‖PAPᵀ − LDLᵀ‖, blocks1.2·10⁻¹⁶‖PAPᵀ − LDLᵀ‖, diagonal1.9·10⁻¹¹growth, blocks1growth, diagonal2.9·10⁵the zero block is what the problem saysand one rule does not need it to be nonzero
Fig. 7 Where the count does come from: a symmetric indefinite factorisation, whose diagonal blocks carry the inertia and which costs a factorisation to obtain.

And it does not certify a direction of descent for the function, only for the model. The quadratic is a local approximation, the direction is unbounded below in the approximation, and how far along it the real function keeps falling is a question the model cannot answer. That is exactly what the trust region is for, and it is why the two halves of this essay are one essay: the certificate is what makes the step possible and the radius is what makes it meaningful.

Where the test can be wrong

Everything above assumes the sign of pᵀAp is the sign of pᵀAp. In floating point it is the sign of a computed inner product, and an inner product whose true value is smaller than its own rounding error has a sign that means nothing.

The condition for that is the condition for a Rayleigh quotient to be lost: the quotient is about λmin, the rounding in the dot product is about u·‖A‖, and the standard backward-error bound for a dot product of length n carries a factor of n. So the sign should be unreliable once λmin/λmax falls below about u·n.

It is worth measuring where it actually goes, because the factor of n in that bound is the same factor eight bits and a format that breaks the rules finds to be absent from the error itself. Taking the exact eigenvector of the negative eigenvalue, rounding it into the format and evaluating the quotient there, over twelve spectra at each of a fine grid of ratios:

bits u u·n median first failure ÷ u u·n ÷ measured
11 4.9·10⁻⁴ 2.0·10⁻² 3.2·10⁻⁵ 0.06 620
24 6.0·10⁻⁸ 2.4·10⁻⁶ 3.2·10⁻⁹ 0.05 750
53 1.1·10⁻¹⁶ 4.4·10⁻¹⁵ 1.0·10⁻¹⁷ 0.09 440

The threshold is about u/15, and it is u/15 at all three. Thirteen orders of magnitude in u, and the measured threshold is 0.06, 0.05 and 0.09 of it — consistent to within a factor of two. So the bound is not the answer scaled by a bit; it is the answer times a constant of about six hundred, and n does not appear.

That changes the practical reading in the direction that matters.

In double the regime is unreachable rather than contrived. The threshold sits at a ratio of 10⁻¹⁷, which is below what a double can hold as a distinct eigenvalue at all. There is no double-precision matrix on which the curvature sign is lost this way, and a failed certificate in double is a bug rather than a rounding limit.

At the reduced precisions the caution stands and moves out by the same factor. Binary32 needs κ ≈ 3·10⁸ rather than the 4·10⁵ the n·u rule implies, and eleven significand bits need κ ≈ 3·10⁴ rather than 50 — which is the difference between “the test is unusable at fp16” and “the test is usable at fp16 on any matrix a solver would not already have refused”.

Rounding the matrix into the format as well as accumulating in it changes nothing on this experiment, which confirms the premise the reasoning started from: the dominant rounding is in the accumulation. assertTheSignSurvivesToAFifteenthOfU measures the whole table and requires the threshold to be the same fraction of u at every precision — because that is the claim that says n is not in it.

The practical form of the caution is short: a certificate should be verified in the precision the claim is made in — which is what two machines, one certificate is about, at a different scale. The vector costs one product to check, the check is independent of the iteration that produced it, and a certificate that fails its own check is a rounding error rather than a proof. That is why the direction is worth returning as a vector rather than as a boolean, and it is the strongest argument for this whole framing: a report cannot be re-checked and an object can.

The cost of not looking

The alternative to all of this is the one most codes take, and it is worth pricing rather than disparaging.

Modify the Hessian until it is positive definite — add a multiple of the identity, which is the regularisation when the answer is a choice prices in another setting, or replace the negative eigenvalues by their absolute values, or run a modified Cholesky that adds to the diagonal whenever a pivot comes out non-positive, which is a factorisation with nothing to pivot for given something to pivot for after all — and then solve. The result is a descent direction, the solver never meets the event, and nothing in the code has to know what indefiniteness is.

What it costs is the direction. A modified Hessian’s Newton step points somewhere sensible and it does not point along the negative curvature, so the method walks away from a saddle at the rate a gradient method walks rather than at the rate the curvature offers. On a problem whose saddles are the interesting feature — which is most non-convex problems of any size — that is the difference between escaping in a few iterations and escaping in many.

The measured version of the same trade is the trust-region curve above: 89 per cent of the exact decrease at a small radius, for two products. A modified-Cholesky step is not on that curve at all, because it is answering the question the model has no answer to by changing the model.

The two meanings, once more

Three essays and three zeros. The first was an object arriving: the Krylov space closed and the answer was inside it. The second was an object failing to arrive: two sequences lost contact and there was no next term. This one is neither — the quantity that changes sign is not a divisor whose vanishing ends a recurrence, it is a measurement of the problem, and what it says is that the problem is not the one the method was written for.

What the three have in common is the shape this collection keeps meeting: the number that stops the method carries information that the method’s own output does not, and whether that information is good news, bad news or a different question entirely is not visible in the number.

The refusal

The assertion behind this essay is fed the claim that a certificate of negative curvature exists on a matrix built to be positive definite.

It is the failure that would make everything above worthless, and it is not far-fetched: it is what happens when the sign of a rounded inner product is read as a fact about the matrix. A routine that accepted it would be manufacturing proofs, and the whole value of the direction as a certificate rests on the claim that a routine which reports one has found one.

The second refusal beside it is smaller and belongs to the trust-region half: a subproblem posed with a region of radius zero, where the model decrease is exactly zero and any share of it is a division by zero. Both are cases where a number comes back that has the right type and no meaning.

Where this goes next

Everything above is about one matrix. The generalised problem Ax = λBx moves the zero into the second matrix and takes an eigenvalue with it, and there a problem can turn out to have no answer at all while every routine that is asked returns one.

Where this goes next

ran to the end3173lucky — a subspace closed295serious, cured by a block of two495serious, cured by a longer block29serious, incurable at any length8counted, not estimatedserious, as a fraction0.13of those, cured at two0.93incurable8matrices tried4000measure zero on the realsand an eighth of the integers The same zero, and nothing was found Change the recurrence by two lines and the divisor stops being a norm. It becomes an inner product of two vectors from two different sequences, and an inner product of two different vectors is zero on a whole hyperplane — with neither vector anywhere near zero, nothing invariant, and nothing converged. The arithmetic event is identical and the meaning is opposite. 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. 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 A factorisation with nothing to pivot for Cholesky's growth factor is not bounded by one. It is equal to one, at every size and every condition number, and the two-line reason is why the algorithm needs no pivoting at all — not "usually gets away without it". Its only failure is the square root of a non-positive number, which is exactly the test for definiteness, and in floating point that test moves with the precision.

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.

CertificateCholeskyConjugate gradientsEigenvaluesIndefinite matrixKrylov subspaceMINRESNegative curvatureTrust-region