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. This site has guarded against it since its first commit 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.
30 min read 24 figures Where the algebra stops being the arithmeticWhose fault is itMeasured, not assertedTwo 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λ_min = -0.1positive: a step existsnegative: a certificate existsone matrix, two questionsstep it turns at6quotient there-0.027share of λ_min recovered0.27λ_min, 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.

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.

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 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 often Cholesky still calls a 8×8 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. 2 The other route to the same verdict, from the elimination field: a Cholesky attempt that fails, which is what “not positive definite” usually means in a library.
The growth factor of a 8×8 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 4 row interchanges, each of which destroys the symmetry that was the reason to use a symmetric factorisation. The bound the general theory allows, 2^7 = 128, is drawn above them both.10¹10³10⁵10⁷10⁹10¹¹110¹10²10³condition number of the matrixgrowth factorbound 2^7partial pivotingCholeskyno pivot to gain fromCholesky growth, every κ1Cholesky interchanges0partial pivoting, worst4the bound, 2^7128both eliminations reach the same growthand only one of them had to swap to get there
Fig. 3 And what that attempt costs when it succeeds, which is the price of finding out.

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.

Matrix–vector products before the curvature test fires, against the size of the negative eigenvalue, n = 40A 40×40 matrix whose spectrum is positive over one decade apart from a single eigenvalue at −λ. At λ = 3 conjugate gradients meets a non-positive curvature after 3 products; at λ = 0.001 it takes 9. The trend is monotone and it runs the reassuring way: a Krylov space finds large eigenvalues first, so the indefiniteness that takes longest to detect is the indefiniteness that matters least. The direction that fires the test recovers between 12 and 42 per cent of λ, so it is a certificate rather than an estimate of the eigenvalue.10⁻³10⁻²10⁻¹10246810size of the negative eigenvalue, −λproducts before the test firesharder to find, and milder8 spectra, n = 40products at the largest λ3products at the smallest9smallest share of λ recovered0.12largest0.42the one that hidesis the one that matters least
Fig. 4 Eight spectra, each positive over one decade apart from a single eigenvalue at −λ. The number of products before the curvature turns, against λ.

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.

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.

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. 5 The reverse arrangement, in the essay that measured it: the smaller the gap, the harder the decision and the more it matters.
Every column's 1-norm on a 24×24 matrix built to stop Hager's walk one column short24 bars, one per column, each the sum of the absolute values in that column — the quantity the estimator maximises, whose largest is the matrix's 1-norm. The walk visits 1 of them and returns 24.00; column 2 has 1-norm 474.24 and is never visited, because its entries alternate in sign and cancel against every sign vector the walk stands on.the estimator maximises this quantity over the columns it visitscolumn 1 ‹visited›24column 2 ‹the answer›474.2column 322.8column 422.8column 522.8column 622.8column 722.8column 822.8column 922.8column 1022.8column 1122.8column 1222.8column 1322.8column 1422.8column 1522.8column 1622.8column 1722.8column 1822.8column 1922.8column 2022.8column 2122.8column 2222.8column 2322.8column 2422.8estimate 24.0a walk that stopped earlythe estimate returned24the true 1-norm474columns visited1products with the matrix5the walk's own stopping test firedand every column it could see was smaller
Fig. 6 And another: a condition estimator whose under-reporting is worst exactly on the matrices where the number would have been worth having.

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.

The condition number of one eigenvalue, at an off-diagonal entry of 1A rising curve against the size of an off-diagonal entry, on logarithmic axes, with a flat line at one below it and a vertical marker at the current value.110¹10²10³110¹10²10³off-diagonal entry ccondition number of the eigenvalue√(1 + c²)decoupled: 1measuredthree routes, one number‖A − ZTZᵀ‖/‖A‖1.4·10⁻¹⁵closed form1.4computed 1/|yᵀx|1.4worst measured movement0.65four eigenvalues, two conditioning numbersthe symmetric case has one, and it is 1
Fig. 7 What computing the eigenvalue properly would involve, and what its own conditioning is.
The plane, and the two vectors inside it, at a gap of 10⁻⁶A circle representing the invariant plane, with two perpendicular radii for the eigenvectors and two more for the perturbed ones, rotated within the same circle.the invariant planesolid: beforedashed: aftersame perturbation, two questionsthe vectors turned, radians0.029the plane turned, radians7.6·10⁻⁸what left the plane5.6·10⁻⁸drawn in the unperturbed plane's own basisa radius is not determined; the circle is
Fig. 8 And the object that is actually well conditioned when the individual eigenvectors are not, which is a subspace rather than a vector.

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.

What the solver reported, and what it achievedFinal relative residual and final relative error against the working precision, on a logarithmic vertical axis. The residual sits at the stopping tolerance at every precision — the method converged, by every test available to it. The error runs from 9.06·10⁻¹³ at 53 bits to 0.0114 at 8, tracking the unit roundoff, which is drawn beside it.614223038465410⁻¹³10⁻¹⁰10⁻⁷10⁻⁴working significand bitsrelative sizeerrorresidualuthe gap the solver cannot seeerror ÷ residual at 24 bits5.9·10⁵error ÷ residual at 16 bits1.2·10⁸error ÷ residual at 8 bits2.7·10¹⁰every convergence test passesand the answer is wrong in proportion to u
Fig. 9 The same shape from the mixed-precision essays: a quantity that decides the outcome and is not among the ones reported.

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⁻¹⁴.

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.

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: κ · u40×40, 20 seeds per κ; dashed is the worstthe problem worsens, not the method
Fig. 10 The gap between a small residual and a correct answer, which is what a method returning a saddle point would be sitting in.

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

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 -10. At Δ = 0.25 the truncated conjugate gradient gets 39.6 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 35.2 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 Δ = 0.25.110¹00.30.60.91.2trust-region radius Δshare of the exact model decreasethe exact subproblem7 radii, n = 40share at the smallest radius0.4share at the largest0.35products, at most2radii stopped by the curvature7a few products against an eigendecompositionand most of the decrease
Fig. 12 At a large negative eigenvalue the curvature stops the step at every radius, and the share is what a single direction is worth.
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.001. At Δ = 0.25 the truncated conjugate gradient gets 89.8 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 68.3 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 Δ = —.110¹00.30.60.91.2trust-region radius Δshare of the exact model decreasethe exact subproblem7 radii, n = 40share at the smallest radius0.9share at the largest0.68products, at most8radii stopped by the curvature0a few products against an eigendecompositionand most of the decrease
Fig. 13 And at a tiny one the boundary stops it first almost everywhere, so the share is what the ordinary conjugate-gradient path is worth.
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 λ_min recovered0.35λ_min, by construction-10MINRES steps on the same system37the division that cannot be doneis the answer to a different question
Fig. 14 The same crossover from the other side: at λ = −10 the quotient turns almost immediately.
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λ_min = -0.001positive: a step existsnegative: a certificate existsone matrix, two questionsstep it turns at9quotient there-2.8·10⁻⁴share of λ_min recovered0.28λ_min, by construction-0.001MINRES steps on the same system39the division that cannot be doneis the answer to a different question
Fig. 15 And at λ = −10⁻³ it takes nine products to get there, by which time the iteration has done most of a solve.

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. 16 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 exactly the condition for a Rayleigh quotient to be lost: the quotient is about λ_min, the rounding in the dot product is about u·‖A‖, so the sign is unreliable once λ_min/λ_max falls below about u·n. At double precision that requires a condition number past 10¹⁵, which is contrived. At the precisions the hardware is moving towards it is not, and the last essay of this collection’s current run measures exactly where.

Relative error of two algebraically identical expressions for (1 − cos x)/x², in binary32A log–log plot of relative error against x. The expression written as it reads loses accuracy as x falls and is eventually wrong in every digit; the rearranged form stays at rounding level.10⁻⁸10⁻⁶10⁻⁴10⁻²110⁻⁹10⁻⁶10⁻³1xrelative error of the computed value(1 − cos x)/x², as written2 sin²(x/2)/x²no digits left at allbinary32 throughoutone function, two spellings · zero below 3.5·10⁻⁴
Fig. 17 The mechanism, from the arithmetic field: what a subtraction of two nearly equal numbers leaves behind, and how much of it is the answer.
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. 18 And the same computation at eight precisions, which is the knob this whole site is built around.

The practical form of the caution is short: a certificate should be verified in the precision the claim is made in. 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 Krawczyk operator on a box of half-width 0.4Two rectangles in the plane. The outer one is the box handed in; the inner one is K(X), the image the operator returns; and the marked point is the root (√2, √2), known in closed form. The image is narrower than the box — 0.2263 against 0.8 — and the verdict is: exactly one root.11.522.511.52xyexactly one roota verdict, not a bound‖I − C F′(X)‖0.28width of X0.8width of K(X)0.23strictly inside is a proofand overlapping is nothing at all
Fig. 19 The strongest version of that habit on this site: a computation that returns an interval containing the answer, which is a claim anybody can check without repeating the work.
The widths a proof is available at, against the precisionTwo boundaries against the number of significand bits, on a logarithmic vertical axis, with the band between them shaded. The narrow boundary runs from 0.00162 at 12 bits to 3.33·10⁻¹⁶ at 53, a small multiple of u at each. The wide boundary is 1.414214 at every precision — which is √2, the half-width at which the box first touches the line where the Jacobian is singular.9172533414910⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹significand bitsbox half-widthwidest provablenarrowest provableutwo boundaries, two causesnarrowest, in units of u3widest, at 53 bits1.4√21.4one end is the arithmeticthe other is the function
Fig. 20 And what verifying rather than estimating costs, measured across precisions.

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

How often Cholesky still calls a 10×10 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. 21 The modification’s own knob, from the elimination field: how far a diagonal has to be pushed before the factorisation succeeds, and what that push has done to the problem.

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 block diagonal of a Bunch–Kaufman factorisation of an 8×8 saddle-point matrix at τ = 10⁻⁶The matrix D from PAPᵀ = LDLᵀ, drawn as a matrix. 4 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.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
Fig. 22 The nearest thing to it in the elimination field: a symmetric matrix with no diagonal entry to pivot on, where the repair is to take two variables at once. The event there is also a zero that names what has to be done differently rather than what has gone wrong.

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.

Matrix–vector products before the curvature test fires, against the size of the negative eigenvalue, n = 20A 20×20 matrix whose spectrum is positive over one decade apart from a single eigenvalue at −λ. At λ = 3 conjugate gradients meets a non-positive curvature after 3 products; at λ = 0.001 it takes 11. The trend is monotone and it runs the reassuring way: a Krylov space finds large eigenvalues first, so the indefiniteness that takes longest to detect is the indefiniteness that matters least. The direction that fires the test recovers between 6 and 33 per cent of λ, so it is a certificate rather than an estimate of the eigenvalue.10⁻³10⁻²10⁻¹1024681012size of the negative eigenvalue, −λproducts before the test firesharder to find, and milder8 spectra, n = 20products at the largest λ3products at the smallest11smallest share of λ recovered0.064largest0.33the one that hidesis the one that matters least
Fig. 23 At twenty unknowns the same trend, with the whole curve a little lower: a smaller Krylov space reaches the extremes of the spectrum in fewer steps because there are fewer eigenvalues between them.
Matrix–vector products before the curvature test fires, against the size of the negative eigenvalue, n = 60A 60×60 matrix whose spectrum is positive over one decade apart from a single eigenvalue at −λ. At λ = 3 conjugate gradients meets a non-positive curvature after 4 products; at λ = 0.001 it takes 11. The trend is monotone and it runs the reassuring way: a Krylov space finds large eigenvalues first, so the indefiniteness that takes longest to detect is the indefiniteness that matters least. The direction that fires the test recovers between 0 and 29 per cent of λ, so it is a certificate rather than an estimate of the eigenvalue.10⁻³10⁻²10⁻¹1024681012size of the negative eigenvalue, −λproducts before the test firesharder to find, and milder8 spectra, n = 60products at the largest λ4products at the smallest11smallest share of λ recovered8.4·10⁻⁴largest0.29the one that hidesis the one that matters least
Fig. 24 And at sixty, where the whole curve has barely moved: a Krylov space sees the extremes of a spectrum whatever the dimension of the space it sits in.

What is next

The three essays so far have all been about one matrix. The next three are about problems with two — Ax = λBx — where the zero moves into the second matrix and takes an eigenvalue with it, and where 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.

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