The matrix a constraint makes

A condition number sent to infinity

An interior-point method manufactures an ill-conditioned matrix on every iteration, deliberately, because the separating of a diagonal is how it discovers which constraints are active. Written one way the answer keeps fifteen digits at a condition number of 3·10¹⁵. Written the other way — the way almost every code writes it — it has none left.

Worth reading first: The zero that is not a missing entry · The condition number is an amplifier · The road that squares the problem.

Every ill-conditioned matrix on this site so far has been ill conditioned because somebody’s problem was. A Hilbert matrix arrives that way. A Vandermonde on equispaced nodes arrives that way. A covariance with a tiny eigenvalue arrives that way. In each case the condition number is a fact about the question, the site’s spine says the answer will be inaccurate in proportion, and the measurement agrees.

An interior-point method for a constrained optimisation problem is different. It manufactures the ill-conditioning, on every iteration, deliberately, and the manufacturing is the algorithm working rather than failing.

Both forms of an interior-point step: two condition numbers that climb together and two errors that do notOne Newton step of an interior-point method for a quadratic programme with 8 unknowns and 6 constraints, 3 of them active, written twice. The augmented form is [[H, Cᵀ], [C, −D⁻¹]] and the condensed form is H + CᵀDC, and they have the same solution. As the barrier parameter μ falls, the diagonal D separates — z/s runs to 1/μ on the active constraints and to μ on the inactive ones — so both condition numbers climb: 29.04 to 3.044·10¹⁵ for the augmented form and 155.8 to 2.403·10¹⁶ for the condensed one, within a factor of 7.9 of each other. The two forward errors, both measured against the exact rational solution of the system the machine is holding, do not follow: the condensed form's rises to 0.3098 and the augmented form's stays at 1.035·10⁻¹⁵ — fifteen correct digits at a condition number of 3.04·10¹⁵.-14-12-10-8-6-4-2010⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹10¹⁵log₁₀ μ — the barrier parametercondition number, and relative errorκ₂, condensedκ₂, augmentederror, condensederror, augmentedagainst a BigInt answerκ₂ augmented, μ = 10⁻¹⁴3·10¹⁵its relative error10⁻¹⁵κ₂ condensed2.4·10¹⁶its relative error0.31the same step, written two waysand only one of them is solvable
Fig. 1 Two writings of one Newton step. Their condition numbers are within a factor of eight of each other at every point on the axis. Their errors are fourteen orders apart at the right-hand end.

Where the infinity comes from

Minimise ½xᵀHxfx subject to Cxd. Add slacks s = Cxd, require them positive, add multipliers z and require those positive too, and then require sᵢzᵢ = μ on every constraint with μ driven towards zero. That set of requirements has a unique solution for each μ, the solutions trace a curve called the central path, and following it is the method.

At the solution the constraints split in two. An active constraint has sᵢ → 0 with zᵢ bounded away from zero. An inactive one has zᵢ → 0 with sᵢ bounded away. So the diagonal Dᵢ = zᵢ/sᵢ that the Newton step carries runs to 1/μ on one group and to μ on the other, and its spread is 1/μ². Measured, at three decades of μ: 1.84·10⁸, 1.84·10¹⁶, 1.84·10²⁴, which is exactly 1/μ² times a constant the problem fixes.

There is no way to avoid this and remain an interior-point method. Which constraints are active is what the method is finding out, and the diagonal separating is how it finds out. The condition number is the algorithm’s progress bar.

What a factorisation survives across one barrier step, against the step's reduction factorFrom one interior-point iteration to the next, H, C and the sparsity pattern are identical and exactly 6 entries of the matrix change — the diagonal of the (2, 2) block, and nothing off it. That is the sparsity pattern reused for free forever. The factorisation is a different question, and the answer is a threshold: carrying the decomposition of K(μ) to K(σμ) and cleaning up with iterative refinement holds at the rounding level while σ is within a per cent of one, and by σ = 0.9 it buys one step. At the schedules an interior-point method actually uses — σ between 0.5 and 0.1 — the first reuse is already at 1.02 and 2.88·10¹¹. The reason is in the entry count: the 6 entries that moved are the 6 that dominate the matrix, and they moved by a factor of 1/σ, so the relative change in K across one step at σ = 0.1 is 9 rather than the 10⁻³ the collection's reuse essays are about.-1-0.75-0.5-0.25010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹10⁴log₁₀ σ — the barrier's reduction factorresidual after one reused stepconvergedthe pattern free, the factors notentries moved6off-diagonal0survived at σ = 0.996survived at σ = 0.106 steps0 stepsthe few entries that movedare the ones that dominate
Fig. 2 The same separation seen from the sequence field, two essays on: what it does to a factorisation carried from one step to the next.

Why the split has to be discovered rather than guessed

It is worth being clear about what an interior-point method is for, because it explains why the ill-conditioning cannot be designed out.

The hard part of an inequality-constrained problem is combinatorial: which of the p constraints hold with equality at the solution. There are 2^p possibilities, the answer determines everything else — once the active set is known the problem is the equality-constrained one the field’s first essay is about — and no amount of linear algebra finds it.

The barrier’s contribution is to replace that combinatorial question by a continuous one. At μ = 1 every constraint is comfortably inactive and the diagonal D is unremarkable. As μ falls the constraints sort themselves: the ones that will bind have their slacks driven towards zero, the ones that will not have their multipliers driven there instead, and by μ = 10⁻⁸ the two groups are sixteen orders apart in D. Reading off which entries are large is reading off the active set.

So the spread of D is the answer being computed, expressed as a number. A method that kept D bounded would be a method that had not yet distinguished the two groups, and the linear algebra’s difficulty is not a side effect of the search — it is the search, written in the matrix.

Both forms of an interior-point step: two condition numbers that climb together and two errors that do notOne Newton step of an interior-point method for a quadratic programme with 8 unknowns and 6 constraints, 4 of them active, written twice. The augmented form is [[H, Cᵀ], [C, −D⁻¹]] and the condensed form is H + CᵀDC, and they have the same solution. As the barrier parameter μ falls, the diagonal D separates — z/s runs to 1/μ on the active constraints and to μ on the inactive ones — so both condition numbers climb: 31.22 to 2.193·10¹⁵ for the augmented form and 102.3 to 2.794·10¹⁶ for the condensed one, within a factor of 12.7 of each other. The two forward errors, both measured against the exact rational solution of the system the machine is holding, do not follow: the condensed form's rises to 0.06652 and the augmented form's stays at 3.49·10⁻¹⁵ — fifteen correct digits at a condition number of 2.19·10¹⁵.-14-12-10-8-6-4-2010⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹10¹⁵log₁₀ μ — the barrier parametercondition number, and relative errorκ₂, condensedκ₂, augmentederror, condensederror, augmentedagainst a BigInt answerκ₂ augmented, μ = 10⁻¹⁴2.2·10¹⁵its relative error3.5·10⁻¹⁵κ₂ condensed2.8·10¹⁶its relative error0.067the same step, written two waysand only one of them is solvable
Fig. 3 Four active constraints of six, where the two groups are as unequal in size as the figure draws.

The two ways to write the step

The Newton step can be written unreduced, as a saddle-point system of the shape this field is about:

[ H Cᵀ ] [Δx] [r₁] [ C −D⁻¹ ] [Δz] = [r₂]

or with the second block eliminated, which is the block elimination the second essay in this field called the range-space method:

(H + CᵀDC) Δx = r₁ + CᵀD r

The second is smaller, it is positive definite, and Cholesky needs no pivoting on it. Almost every implementation solves that one — it is the condensed or normal-equations form, and the name is not a coincidence: the elimination is exactly the move that turns a least-squares problem into its normal equations, and this site has priced that move once.

Both are the same step. They have the same solution, they are related by an exact algebraic identity, and there is nothing to choose between them on paper.

Least squares by QR and by the normal equations in binary32Relative error of the computed coefficients against epsilon, on log axes. The QR route is a flat line near the bottom; the normal-equations route climbs and then stops, at the epsilon where the cross-product matrix becomes exactly singular.10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹ε in Läuchli's matrix (smaller ε, larger κ)relative error in the coefficientsAᵀA exactly singularnormal equationsQRκ from 1.7·10⁸ to 17the cliff is at √u = 2.4·10⁻⁴
Fig. 4 The move, in the field it was first priced in.

The measurement

Both systems, at the same μ, on the same problem, each solved by the routine a code would use: the augmented one by a symmetric indefinite factorisation with 2 × 2 pivots, the condensed one by a Cholesky. Both answers compared against an exact rational solution of the same stored matrix.

μ κ₂ augmented κ₂ condensed error augmented error condensed 10⁰ 2.90·10¹ 1.56·10² 7.6·10⁻¹⁶ 3.7·10⁻¹⁵ 10⁻⁴ 3.02·10⁵ 2.49·10⁶ 1.1·10⁻¹⁵ 9.4·10⁻¹² 10⁻⁸ 3.04·10⁹ 2.51·10¹⁰ 1.2·10⁻¹⁵ 3.0·10⁻⁷ 10⁻¹² 3.04·10¹³ 2.52·10¹⁴ 9.4·10⁻¹⁶ 1.0·10⁻³ 10⁻¹⁴ 3.04·10¹⁵ 2.40·10¹⁶ 1.0·10⁻¹⁵ 3.1·10⁻¹

The two condition numbers climb together, thirteen orders each, and stay within a factor of eight of one another the whole way. The condensed error follows κ·u exactly and arrives at 31 per cent. The augmented error does not move at all: it is 10⁻¹⁵ at μ = 1 and 10⁻¹⁵ at μ = 10⁻¹⁴, with a condition number of 3·10¹⁵.

Fifteen correct digits at a condition number of 3·10¹⁵. That is the finding.

Both forms of an interior-point step: two condition numbers that climb together and two errors that do notOne Newton step of an interior-point method for a quadratic programme with 8 unknowns and 6 constraints, 1 of them active, written twice. The augmented form is [[H, Cᵀ], [C, −D⁻¹]] and the condensed form is H + CᵀDC, and they have the same solution. As the barrier parameter μ falls, the diagonal D separates — z/s runs to 1/μ on the active constraints and to μ on the inactive ones — so both condition numbers climb: 28.23 to 1.439·10¹⁶ for the augmented form and 162.7 to 4.182·10¹⁷ for the condensed one, within a factor of 29.1 of each other. The two forward errors, both measured against the exact rational solution of the system the machine is holding, do not follow: the condensed form's rises to 1 and the augmented form's stays at 1.352·10⁻¹⁵ — fifteen correct digits at a condition number of 1.44·10¹⁶.-14-12-10-8-6-4-2010⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹10¹⁵log₁₀ μ — the barrier parametercondition number, and relative errorκ₂, condensedκ₂, augmentederror, condensederror, augmentedagainst a BigInt answerκ₂ augmented, μ = 10⁻¹⁴1.4·10¹⁶its relative error1.4·10⁻¹⁵κ₂ condensed4.2·10¹⁷its relative error1the same step, written two waysand only one of them is solvable
Fig. 5 With one active constraint rather than three, where the same separation happens and the condensed Cholesky additionally stops working.

The refusal this page publishes

The reading that has to be closed is not a careless one. It is the site’s own spine, quoted correctly: forward error ⪅ condition number × backward error, so a backward-stable solve of a matrix with κ = 10¹³ has a relative forward error of about 10¹³ · 10⁻¹⁶ = 10⁻³, and three digits survive.

The bound is true. It is an upper bound, it is attained by the condensed form on this very page, and the augmented form is twelve orders under it.

The assertion is fed the pair at μ = 10⁻¹² — κ = 3.04·10¹³ and a measured relative error of 9.4·10⁻¹⁶ — and required to reject the claim that the error is near κ·u. It does. What that leaves is a question the rest of the field has to answer: if κ₂ is not the number describing this matrix, what is?

The next essay in the error field has it, and it is a number that has been in this collection since the scaling essays: the componentwise condition number of the augmented matrix is 13.3, and it does not move while κ₂ crosses thirteen decades.

The augmented step's two condition numbers, and the one the error obeysκ₂ is a worst case over all perturbations of the same NORM, and a normwise perturbation is allowed to put its whole budget on the smallest entry of a matrix. The componentwise number — Skeel's ‖ |A⁻¹||A||x| ‖ / ‖x‖ — is a worst case over perturbations proportional to the entries, which is what a backward-stable factorisation actually makes. For the augmented interior-point matrix the two are 3.044·10¹³ and 13.25 at μ = 10⁻¹², a ratio of 2.296·10¹², and the measured error is 9.434·10⁻¹⁶ — which is the second number times the unit roundoff and has nothing to do with the first. The componentwise line is still finding its level down to about μ = 10⁻⁴, where the two groups of constraints have not yet separated, and from there it does not move in the fourth digit; where it settles — 13.25 here — is set by how many constraints are active rather than by μ. Eliminating the second block to reach the condensed form destroys the distinction: there the componentwise number is 1.324·10¹⁴, tracking the normwise one, because the large entries have been summed into CᵀDC and are no longer separately identifiable.-12-10-8-6-4-2010⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹10¹⁵log₁₀ μcondition number, and relative errorκ₂, augmentedcomponentwise, condensedcomponentwise, augmentederror, augmentedone matrix, two numbersκ₂ at μ = 10⁻¹²3·10¹³componentwise, same matrix13their ratio2.3·10¹²measured error9.4·10⁻¹⁶every library prints the top lineand the error obeys the third
Fig. 6 The answer, which the next essay is about: two condition numbers of one matrix, thirteen orders apart, and the error obeying the smaller one.

Why the structure survives one form and not the other

The mechanism is the one the structured backward error essays established, applied to a matrix whose structure is unusually explicit.

A backward-stable factorisation perturbs each entry by an amount relative to that entry. The huge entries of −D⁻¹ are perturbed by huge absolute amounts that are tiny relative to themselves; the ordinary entries of H and C are perturbed by ordinary tiny amounts. So the perturbation the arithmetic actually makes is structured: it respects the scale of every entry.

A normwise condition number cannot see that. κ₂ is a worst case over perturbations of a given norm, and a normwise perturbation is permitted to put its entire budget on the smallest entry of the matrix — a perturbation no backward-stable factorisation would ever make. That worst case is what κ₂ = 3·10¹⁵ is reporting, and it is a fact about a perturbation nobody is applying.

Eliminating the second block destroys the distinction, and that is why the condensed form gets no protection. The huge entries are summed into CᵀDC and stop being separately identifiable: an entry of the condensed matrix is a sum of a small number and a huge one, and perturbing it relatively perturbs the small contribution by an absolute amount that is huge. The structure was in the sparsity of the scale, and the elimination smeared it.

The backward error of a hierarchical solve, against the error of the representation it usedSolve A_H x = b exactly and the residual against the matrix that was wanted is b − Ax = (A_H − A)x, so ‖b − Ax‖ ⁄ ‖A‖‖x‖ cannot be anything but the representation's own error. Measured across six accuracies spanning ten decades, the two track at a slope of 0.994 and sit a constant 6.7× apart, which is the difference between a norm of a matrix and a norm of that matrix applied to one vector. The reading is the one this field is built on: the accuracy knob is not an accuracy, it is a backward error chosen in advance. Everywhere else on this site a backward error is something an algorithm produced and somebody then measured; here it is a line in the program, and its size is known before the solve starts.10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³‖A − A_H‖ ⁄ ‖A‖, the representation‖b − Ax‖ ⁄ ‖A‖‖x‖, the solveequala backward error, chosenslope0.99representation at 10⁻⁸1.2·10⁻⁹backward error there1.3·10⁻¹⁰representation at 10⁻¹²1.2·10⁻¹³backward error there2.5·10⁻¹⁴the compression is not an approximationit is a perturbation of the problem
Fig. 7 The identity being used correctly, in the field that owns it — where the backward error is what the algorithm did and the condition number is what the problem did with it.

And the condensed form does something worse than lose digits

H + CᵀDC is positive definite for every positive D — H is, CᵀDC is positive semidefinite, and a sum of the two is. It is a theorem with no hypotheses left to check.

At μ = 10⁻¹⁴ with one or two active constraints the Cholesky refuses it. A pivot comes out non-positive and the routine this collection uses as the test for definiteness answers no about a matrix whose smallest eigenvalue is 0.0107.

So at the far end of the barrier the condensed form has not merely lost its digits; it has no answer to return. The library asserts both halves — every eigenvalue positive, and the factorisation refusing — because the interesting claim is the pair rather than either alone. The algebra and the arithmetic disagree about a property, not about a value.

How often Cholesky still calls a 12×12 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. 8 Where that boundary sits in general, from the elimination field: a matrix definite by a few ulps and a Cholesky that fails anyway.

The one place the two forms are genuinely different problems

There is a difference between the two systems that is not about accuracy at all, and it is worth separating from everything above.

The augmented system returns Δx and Δz together. The condensed system returns Δx, and Δz has to be recovered by back-substitution: Δz = D(CΔxr₂). That recovery multiplies by D, whose entries reach 1/μ, so whatever relative error Δx carries is inherited by Δz — and the multipliers are what the method uses to decide the next step length and to test optimality.

So the condensed form’s loss is not confined to the primal step. It arrives in the dual variables amplified, and a method that is losing digits in z near the solution is a method whose stopping test is being computed from the least accurate quantity it has. The augmented form computes both blocks in one solve and neither is downstream of the other.

That is the sort of asymmetry a table of forward errors on Δx alone would hide, and it is one reason the measurement above compares the full stacked solution rather than its first block.

‖AX̂ − I‖ and ‖X̂A − I‖ for a computed Hilbert inverse, to n = 12Two residuals of the same computed inverse, each divided by ‖A‖‖X̂‖ so both are dimensionless. The exact inverse satisfies both to zero. The computed one satisfies ‖AX̂ − I‖ at 3.5·10⁻¹⁸ — the side its columns were solved along — and ‖X̂A − I‖ at 6.8·10⁻¹⁶, 195 times larger, at n = 12. Each column of X̂ is the exact solution of a slightly perturbed system, but a different perturbation for each column, so there is no single nearby matrix whose inverse X̂ is.468101210⁻²⁰10⁻¹⁹10⁻¹⁸10⁻¹⁷10⁻¹⁶10⁻¹⁵10⁻¹⁴10⁻¹³nresidual / (‖A‖ ‖X̂‖)‖X̂A − I‖‖AX̂ − I‖the two sides, comparedn = 6, ratio6.4n = 8, ratio8.3n = 10, ratio82n = 12, ratio195the exact inverse satisfies bothand the computed one satisfies the side it was computed along
Fig. 9 The same shape of error in the field that owns it: a quantity computed through an intermediate, inheriting the intermediate’s conditioning rather than its own.

What a code should do, and what most do

The augmented form is (n + p) × (n + p) rather than n × n, needs a symmetric indefinite factorisation with 2 × 2 pivots rather than a Cholesky, and on a sparse problem has a fill that depends on an ordering that has to cope with pivoting. The condensed form is smaller, definite, and its ordering can be computed once. Every one of those is a real advantage and they are why the condensed form is the default in most implementations.

What the measurement says is that the advantage has a price, that the price is invisible in the condition number, and that it arrives late — at μ = 10⁻⁴ the condensed error is 10⁻¹¹ and nobody would notice. The digits go where the method is trying to finish.

There is a repair, and it is the sixth essay in this field: regularise the augmented form so that its factorisation exists under any ordering, and the sparse objection to it disappears. That is why interior-point codes that do use the augmented form regularise, and it is why the regularisation is not a hack.

What the regularisation costs, and what 6 steps of refinement take backSolving [[H + δI, Aᵀ], [A, −δI]] instead of K gives the exact answer to a different problem, so its error is proportional to δ: measured at 1489·δ across six decades, which is a slope of one and not a trend. Refining against the unregularised matrix — the residual formed with K and the correction solved with the regularised factorisation — removes that term entirely, because the perturbation was never in the residual. It works while δ is below σ_min(K) = 6.797·10⁻⁴, marked on the axis, and stops working above it: the iteration's contraction factor is δ/σ_min and a fixed point needs that under one. So the trade-off curve every regularisation essay on this site has drawn — a term falling in δ against a term rising in it — has, here, a left branch that can simply be removed.-14-12-10-8-6-4-210⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹log₁₀ δrelative error against the exact answerδ = σ_min(K)no refinement6 stepsa left branch that can be removederror ÷ δ, unrefined1489σ_min(K)6.8·10⁻⁴refined at δ = 10⁻⁶3.6·10⁻¹⁶refined at δ = 10⁻²0.62the perturbation is known exactlybecause the code chose it
Fig. 10 The repair that makes the accurate form affordable, from the essay that measures what it costs.

The measurement’s own foundations

A page claiming fifteen correct digits at κ = 10¹⁵ has to say what it compared against, because the obvious method of comparison would be circular: solving the same system more carefully and calling that the truth is exactly the practice this collection was built to avoid.

Every finite double is a dyadic rational. So the matrix the machine is holding — after the rounding that happened when H, C and D were formed — is a matrix of exact rationals, and exactSolve solves it in BigInt with no rounding anywhere. The comparison is against the answer to the problem the machine actually has, not against the answer to the problem somebody meant, and the difference between the two is not a numerical question at all.

That distinction matters more here than anywhere else on the site. The problem somebody meant has D exactly equal to z/s; the problem the machine has has D rounded, and at μ = 10⁻¹⁴ the rounding of a quantity of size 10¹⁴ is a change of size 10⁻². Measuring against the intended problem would report an error of 10⁻² for both forms and conclude nothing. Measuring against the stored problem separates the arithmetic of the solve from the arithmetic of the assembly, and it is the second that the whole page depends on being able to ignore.

The exact solution of a 13×13 Hilbert system beside the computed oneTwo columns of numbers: the exact answer, which is the integers one to thirteen, and the answer double-precision elimination returns, with the number of correct digits beside each.H13 x = b, b formed exactly so that x = (1, 2, …, 13)12345678910111213exact1.00002.00003.00133.97865.18864.993010.46720.048921.2657-2.575319.21478.906013.5113computed6.6 correct digits4.8 correct digits3.4 correct digits2.3 correct digits1.4 correct digit0.8 correct digitno correct digitsno correct digitsno correct digitsno correct digitsno correct digits0.6 correct digit1.4 correct digitbackward error2.2·10⁻¹⁷κ = 1.7·10¹⁸The algorithm solved a neighbouring problem perfectly. That problem's answer is this one.right-hand side built in BigInt rationalsthe truth is known
Fig. 11 The older form of the same discipline, where the entries were rationals before they were stored.

What is not claimed

Three things this page does not say, and separating them from what it does say is most of the value.

It does not say the augmented form is always accurate. It says that on this family, at these μ, with this factorisation, it is — and the reason it is has been isolated to a property that is measurable, so it can be checked elsewhere rather than assumed.

It does not say κ₂ is useless. It is the right number for a perturbation of a given norm, which is what an error in the data is: a measurement error in H or C is not structured, and against that kind of error the augmented system really is as sensitive as κ₂ says.

And it does not say the barrier is a bad idea. The ill-conditioning is what the method uses to locate the active set, and a method that avoided it would be a method that never found out which constraints bind.

How much a perturbation of the right-hand side is amplified, κ = 10⁶The cumulative distribution of the amplification factor over two hundred random perturbation directions, with the condition number marked as the upper limit.110¹10²10³10⁴10⁵10⁶10⁷00.250.50.751amplification of the input perturbationfraction of directions at or belowκ = 10·10⁵worst found 7.6·10⁵6×6, 200 directionsmedian reaches 0.29 of κ
Fig. 12 The bound doing its job, in the case it was written for.

At other settings

Both forms of an interior-point step: two condition numbers that climb together and two errors that do notOne Newton step of an interior-point method for a quadratic programme with 8 unknowns and 6 constraints, 2 of them active, written twice. The augmented form is [[H, Cᵀ], [C, −D⁻¹]] and the condensed form is H + CᵀDC, and they have the same solution. As the barrier parameter μ falls, the diagonal D separates — z/s runs to 1/μ on the active constraints and to μ on the inactive ones — so both condition numbers climb: 28.34 to 1.401·10¹⁶ for the augmented form and 160.6 to 1.673·10¹⁷ for the condensed one, within a factor of 11.9 of each other. The two forward errors, both measured against the exact rational solution of the system the machine is holding, do not follow: the condensed form's rises to 1 and the augmented form's stays at 5.024·10⁻¹⁵ — fifteen correct digits at a condition number of 1.4·10¹⁶.-14-12-10-8-6-4-2010⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹10¹⁵log₁₀ μ — the barrier parametercondition number, and relative errorκ₂, condensedκ₂, augmentederror, condensederror, augmentedagainst a BigInt answerκ₂ augmented, μ = 10⁻¹⁴1.4·10¹⁶its relative error5·10⁻¹⁵κ₂ condensed1.7·10¹⁷its relative error1the same step, written two waysand only one of them is solvable
Fig. 13 Two active constraints, where the ill-conditioning lives in two directions rather than three.
Both forms of an interior-point step: two condition numbers that climb together and two errors that do notOne Newton step of an interior-point method for a quadratic programme with 8 unknowns and 6 constraints, 5 of them active, written twice. The augmented form is [[H, Cᵀ], [C, −D⁻¹]] and the condensed form is H + CᵀDC, and they have the same solution. As the barrier parameter μ falls, the diagonal D separates — z/s runs to 1/μ on the active constraints and to μ on the inactive ones — so both condition numbers climb: 31.53 to 2.161·10¹⁵ for the augmented form and 97.8 to 1.482·10¹⁶ for the condensed one, within a factor of 6.86 of each other. The two forward errors, both measured against the exact rational solution of the system the machine is holding, do not follow: the condensed form's rises to 0.2882 and the augmented form's stays at 4.805·10⁻¹⁵ — fifteen correct digits at a condition number of 2.16·10¹⁵.-14-12-10-8-6-4-2010⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹10¹⁵log₁₀ μ — the barrier parametercondition number, and relative errorκ₂, condensedκ₂, augmentederror, condensederror, augmentedagainst a BigInt answerκ₂ augmented, μ = 10⁻¹⁴2.2·10¹⁵its relative error4.8·10⁻¹⁵κ₂ condensed1.5·10¹⁶its relative error0.29the same step, written two waysand only one of them is solvable
Fig. 14 Five, where nearly every constraint binds and the inactive group is a single row.
The augmented step's two condition numbers, and the one the error obeysκ₂ is a worst case over all perturbations of the same NORM, and a normwise perturbation is allowed to put its whole budget on the smallest entry of a matrix. The componentwise number — Skeel's ‖ |A⁻¹||A||x| ‖ / ‖x‖ — is a worst case over perturbations proportional to the entries, which is what a backward-stable factorisation actually makes. For the augmented interior-point matrix the two are 2.52·10¹³ and 11.19 at μ = 10⁻¹², a ratio of 2.252·10¹², and the measured error is 6.091·10⁻¹⁶ — which is the second number times the unit roundoff and has nothing to do with the first. The componentwise line is still finding its level down to about μ = 10⁻⁴, where the two groups of constraints have not yet separated, and from there it does not move in the fourth digit; where it settles — 11.19 here — is set by how many constraints are active rather than by μ. Eliminating the second block to reach the condensed form destroys the distinction: there the componentwise number is 2.759·10¹⁴, tracking the normwise one, because the large entries have been summed into CᵀDC and are no longer separately identifiable.-12-10-8-6-4-2010⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹10¹⁵log₁₀ μcondition number, and relative errorκ₂, augmentedcomponentwise, condensedcomponentwise, augmentederror, augmentedone matrix, two numbersκ₂ at μ = 10⁻¹²2.5·10¹³componentwise, same matrix11their ratio2.3·10¹²measured error6.1·10⁻¹⁶every library prints the top lineand the error obeys the third
Fig. 15 The componentwise reading at ten constraints, where the flat line rises with the count and not with μ.
The 14 eigenvalues of a saddle-point matrix with 10 unknowns and 4 constraints, inside their closed-form bracketsK = [[H, Aᵀ], [A, 0]] with κ(H) = 100 and κ(A) = 100. The shaded bands are the Rusten–Winther brackets, computed from four numbers — the extreme eigenvalues of H and the extreme singular values of A — before the matrix was assembled: negative eigenvalues in [-0.995, -10·10⁻⁵] and positive ones in [0.01, 1.62]. The marks are the computed spectrum. There are exactly 10 above zero and 4 below, which is Sylvester's law of inertia and not a property of this matrix: K is congruent to blkdiag(H, −AH⁻¹Aᵀ), both blocks are definite, and congruence preserves signs. The two groups are separated by a gap containing zero, at a ratio of 39.3 between the innermost positive eigenvalue and the innermost negative one.-1-0.582271-0.1645420.2531860.6709151.088641.506370eigenvalue4 negative10 positivecounted before it was formedpositive10negative4at zero0innermost ratio39the zero block is a theoremand so is the count either side of it
Fig. 16 The augmented form’s inertia, which is the field’s first essay applied to a matrix whose (2, 2) block is not quite zero.
What each way of eliminating a constraint inherits, over five decades of κ(A)The same system solved twice, at 8 unknowns and 3 constraints with κ(H) = 100. The range-space method forms S = AH⁻¹Aᵀ and inherits κ(S), which rises from 12.99 to 3.981·10¹⁰ — the square of κ(A), for the reason the normal equations square it. The null-space method solves with the reduced Hessian ZᵀHZ, whose condition number is 21.13 at the start of the sweep and 21.13 at the end: it does not contain κ(A) at all. The two forward errors, measured against a solution computed in BigInt rationals, follow their own condition numbers: 5.314·10⁻⁶ against 5.788·10⁻¹² at the far end. Both methods are correct and one of them is usable.01234510⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹log₁₀ κ(A)condition number, and relative errorκ(S)κ(ZᵀHZ)range-space errornull-space erroragainst a BigInt answerκ(S) at κ(A) = 10⁵4·10¹⁰κ(ZᵀHZ), all stops21range-space forward error5.3·10⁻⁶null-space forward error5.8·10⁻¹²both are the same algebraand only one squares
Fig. 17 The elimination that produces the condensed form, priced on a problem with no barrier in it.
Two perturbation bounds and the error that was measured, on a 8×8 matrix spread over 8 decades of unitsThree curves against the size of an entrywise relative perturbation. The normwise bound κ_∞·ε is a valid bound and sits 4·10⁷ times above the componentwise one cond(A, x)·ε, which is also a bound and is nearly attained by the worst of forty random perturbations at each size.10⁻¹⁴10⁻¹³10⁻¹²10⁻¹¹10⁻¹⁰10⁻⁹10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1relative size of the entrywise perturbationrelative forward errorκ_∞ · εcond(A,x) · εmeasuredboth bounds holdκ_∞(A)1.9·10⁸cond(A, x)4.8ratio of the bounds4·10⁷both curves above the data are boundsand only one of them is a measurement
Fig. 18 The componentwise bound in the field where it was introduced.
κ(A) and κ of the row-equilibrated matrix, along both families, n = 8Write A = D·X with D diagonal and every row of X of unit norm. Every theorem about relative accuracy is a hypothesis on κ(X), and κ(A) appears in none of them — which is easy to read past and is the whole difference between the two families here. Along the graded family κ(A) climbs from 271 to 7.04·10⁵⁰ and κ(X) is 4.892 at every one of the six matrices — the grading is exactly what the diagonal factor absorbs. Along the uniform family κ(X) climbs to 3.87·10²⁵. That is the number that says which question has an answer, and it is not the number anybody prints.01020304050110¹⁰10²⁰10³⁰10⁴⁰10⁵⁰log₁₀ κ(A)log₁₀ of the condition numberκ(X) = κ(A)graded familyuniform familyone number decides, and it is not κκ(X), graded, at every grading4.9κ(A), graded, at the widest7·10⁵⁰κ(X), uniform, at the widest3.9·10²⁵κ(A) there8·10²⁸κ is a fact about the matrixand the hypothesis is about a factor of it
Fig. 19 And what a change of units does to a normwise condition number, which is the same fact from the other side.
Backward error with and without the structure, for Levinson and for elimination, n = 12Both solvers are backward stable on this family by the usual measure: the smallest perturbation of any kind that explains either answer stays below 3.1·10⁻¹⁷ at every ρ, and the two curves are indistinguishable. Insisting the perturbation be a symmetric Toeplitz matrix — the kind of object the problem was posed with — lifts both by orders, and at ρ = 0.999 Levinson's rises by a factor of 2.48·10⁶. Neither number is alarming in absolute terms; what is worth carrying is that the reassuring one is the one that is reported, and the two are not measuring the same thing.10⁻³10⁻²10⁻¹110⁻¹⁹10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹1 − ρbackward errorLevinson, Toeplitzelimination, ToeplitzLevinson, any kindelimination, any kindat ρ = 0.999Levinson, any perturbation3.1·10⁻¹⁷Levinson, Toeplitz only7.6·10⁻¹¹the ratio between them2.5·10⁶diagonal defect of the first0.97the number that is reportedand the number that is asked about
Fig. 20 Backward errors of several solvers on one matrix, from the error field.
Iterative refinement from a 24-bit factorisation, κ = 10⁴A semi-logarithmic plot of forward error against refinement step. One curve falls steeply to the level of a double-precision solve; the other is nearly flat.012345610⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹refinement step‖x − x*‖ / ‖x*‖a full double-precision solveresidual in24-bitresidual indoubleone argument apartκ·u of the factorisation6·10⁻⁴double residual, final3.2·10⁻¹³same-precision, final1.3·10⁻⁴30×30, κ = 10⁴, same factors in both runsidentical cost
Fig. 21 The repair a condensed solve would need, and which cannot recover digits the factorisation never had.
What a factorisation survives across one barrier step, against the step's reduction factorFrom one interior-point iteration to the next, H, C and the sparsity pattern are identical and exactly 4 entries of the matrix change — the diagonal of the (2, 2) block, and nothing off it. That is the sparsity pattern reused for free forever. The factorisation is a different question, and the answer is a threshold: carrying the decomposition of K(μ) to K(σμ) and cleaning up with iterative refinement holds at the rounding level while σ is within a per cent of one, and by σ = 0.9 it buys one step. At the schedules an interior-point method actually uses — σ between 0.5 and 0.1 — the first reuse is already at 0.0527 and 1.49·10¹⁰. The reason is in the entry count: the 4 entries that moved are the 4 that dominate the matrix, and they moved by a factor of 1/σ, so the relative change in K across one step at σ = 0.1 is 9 rather than the 10⁻³ the collection's reuse essays are about.-1-0.75-0.5-0.25010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹10⁴log₁₀ σ — the barrier's reduction factorresidual after one reused stepconvergedthe pattern free, the factors notentries moved4off-diagonal0survived at σ = 0.996survived at σ = 0.106 steps0 stepsthe few entries that movedare the ones that dominate
Fig. 22 What the same separation does to reuse, from two essays on.
The growth factor of a 12×12 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 9 row interchanges, each of which destroys the symmetry that was the reason to use a symmetric factorisation. The bound the general theory allows, 2^11 = 2048, is drawn above them both.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
Fig. 23 Cholesky’s growth factor, which is exactly one — and says nothing about whether it will complete.
The spectrum of P⁻¹K with S = AH⁻¹Aᵀ, exactly, at 10 unknowns and 4 constraintsP = blkdiag(H, S). With S the exact Schur complement AH⁻¹Aᵀ the preconditioned matrix has exactly three distinct eigenvalues — 1 with multiplicity n − m = 6, and (1 ± √5)/2 with multiplicity 4 each. Those are 1 − φ = -0.618034 and φ = 1.61803, the golden ratio, which arrives from λ² − λ − 1 = 0 rather than from anything anybody chose. The dashed lines are that closed form and the marks are the computed spectrum; here they are 3 distinct values and the largest distance from the closed form anywhere is 2.909·10⁻¹⁴. A minimal polynomial of degree three means a Krylov method finishes in three steps, which is what the next figure measures.00.3670080.7340171.101031.468031.835040eigenvalue of P⁻¹Kwritten down, then computeddistinct3at 16φ computed1.6off the closed form2.9·10⁻¹⁴1 − φ1φthe preconditioner's effect is a theoremand the golden ratio is in it
Fig. 24 And the preconditioner two essays back, whose spectrum this matrix’s barrier would move off its closed form.
The augmented step's two condition numbers, and the one the error obeysκ₂ is a worst case over all perturbations of the same NORM, and a normwise perturbation is allowed to put its whole budget on the smallest entry of a matrix. The componentwise number — Skeel's ‖ |A⁻¹||A||x| ‖ / ‖x‖ — is a worst case over perturbations proportional to the entries, which is what a backward-stable factorisation actually makes. For the augmented interior-point matrix the two are 1.822·10¹³ and 15.57 at μ = 10⁻¹², a ratio of 1.17·10¹², and the measured error is 6.124·10⁻¹⁶ — which is the second number times the unit roundoff and has nothing to do with the first. The componentwise line is still finding its level down to about μ = 10⁻⁴, where the two groups of constraints have not yet separated, and from there it does not move in the fourth digit; where it settles — 15.57 here — is set by how many constraints are active rather than by μ. Eliminating the second block to reach the condensed form destroys the distinction: there the componentwise number is 4.19·10¹⁴, tracking the normwise one, because the large entries have been summed into CᵀDC and are no longer separately identifiable.-12-10-8-6-4-2010⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹10¹⁵log₁₀ μcondition number, and relative errorκ₂, augmentedcomponentwise, condensedcomponentwise, augmentederror, augmentedone matrix, two numbersκ₂ at μ = 10⁻¹²1.8·10¹³componentwise, same matrix16their ratio1.2·10¹²measured error6.1·10⁻¹⁶every library prints the top lineand the error obeys the third
Fig. 25 The explanation at five constraints.
Both forms of an interior-point step: two condition numbers that climb together and two errors that do notOne Newton step of an interior-point method for a quadratic programme with 8 unknowns and 8 constraints, 3 of them active, written twice. The augmented form is [[H, Cᵀ], [C, −D⁻¹]] and the condensed form is H + CᵀDC, and they have the same solution. As the barrier parameter μ falls, the diagonal D separates — z/s runs to 1/μ on the active constraints and to μ on the inactive ones — so both condition numbers climb: 28.05 to 2.822·10¹⁵ for the augmented form and 96.2 to 2.275·10¹⁶ for the condensed one, within a factor of 8.06 of each other. The two forward errors, both measured against the exact rational solution of the system the machine is holding, do not follow: the condensed form's rises to 0.1919 and the augmented form's stays at 5.753·10⁻¹⁶ — fifteen correct digits at a condition number of 2.82·10¹⁵.-14-12-10-8-6-4-2010⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹10¹⁵log₁₀ μ — the barrier parametercondition number, and relative errorκ₂, condensedκ₂, augmentederror, condensederror, augmentedagainst a BigInt answerκ₂ augmented, μ = 10⁻¹⁴2.8·10¹⁵its relative error5.8·10⁻¹⁶κ₂ condensed2.3·10¹⁶its relative error0.19the same step, written two waysand only one of them is solvable
Fig. 26 Eight constraints, where the inactive group is larger and the picture is unchanged.
500 random symmetric orderings of a regularised saddle-point matrix, and what each factorisation reproducesA quasi-definite matrix — [[H + δI, Aᵀ], [A, −δI]] with δ = 10⁻⁸ — has an LDLᵀ factorisation with diagonal D for every symmetric permutation, so the ordering can be chosen for fill with no numerical veto at all. That is the count: 500 of 500 orderings factorise here, against 69.0 per cent for the same matrix with the zero block left where the problem put it. Each bar counts the orderings whose factorisation left a relative residual ‖PKPᵀ − LDLᵀ‖/‖K‖ in that decade, with the tallest holding 130 of the 500. They span 1.65·10⁻¹⁶ to 6.17·10⁻⁴ — 13 orders — and the growth factor across them runs from 1 to 6.43·10⁷ — which is 0.643/δ, an inverse law that holds at every δ this figure is drawn at. Existence is not stability, and the theorem says only the first.-18-15-12-9-6-300log₁₀ ‖PKPᵀ − LDLᵀ‖ / ‖K‖share of orderingsbestworstexistence and stabilityfactorise1unregularised0.69worst growth6.4·10⁷growth × δ0.64the ordering is free to chooseand not free of consequence
Fig. 27 The regularisation that makes the accurate form affordable.
What each way of eliminating a constraint inherits, over five decades of κ(A)The same system solved twice, at 8 unknowns and 2 constraints with κ(H) = 100. The range-space method forms S = AH⁻¹Aᵀ and inherits κ(S), which rises from 11.43 to 1.127·10¹¹ — the square of κ(A), for the reason the normal equations square it. The null-space method solves with the reduced Hessian ZᵀHZ, whose condition number is 55.28 at the start of the sweep and 55.28 at the end: it does not contain κ(A) at all. The two forward errors, measured against a solution computed in BigInt rationals, follow their own condition numbers: 2.137·10⁻⁶ against 1.053·10⁻¹¹ at the far end. Both methods are correct and one of them is usable.01234510⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹log₁₀ κ(A)condition number, and relative errorκ(S)κ(ZᵀHZ)range-space errornull-space erroragainst a BigInt answerκ(S) at κ(A) = 10⁵1.1·10¹¹κ(ZᵀHZ), all stops55range-space forward error2.1·10⁻⁶null-space forward error1.1·10⁻¹¹both are the same algebraand only one squares
Fig. 28 And the elimination that squares, on a problem with no barrier in it.
ν(σ), the number of eigenvalues below σ, counted from the signs of an unpivoted LDLᵀThe staircase is the number of negative pivots in an LDLᵀ factorisation of A − σI, taken with no pivoting at all. Congruence preserves the signs, so that count is the number of eigenvalues below σ — Sylvester's law of inertia, used as an algorithm. The vertical marks are the 14 eigenvalues a Jacobi decomposition returns, which is a completely different computation; every step of the staircase is at one of them and every one of them has a step. Over 2000 shifts placed at random across the interval the two answers disagree 0 times. That is not a statement about a tolerance: the output is a count, so it is exactly right or wrong by a whole eigenvalue, and there is nothing in between for a rounding error to land in.-4-2.35371-0.7074290.9388572.585144.231435.8777102468101214shift σν(σ)an answer that is an integershifts2000disagreements0eigenvalues14steps14the marks are a Jacobi decompositionand the staircase never saw one
Fig. 29 An inertia counted rather than computed, from the spectra field.
The reduced Hessian and the answer, against how badly the basic columns were chosenγ is the condition number of the first 4 columns of A, and the naive rule calls exactly those columns basic. Its reduced Hessian's condition number climbs from 83.77 to 4.871·10¹⁸ across ten decades of γ — the square, because κ(ZᵀHZ) ≤ κ(H)κ(Z)² and the bound is attained. Choosing the basic columns by a column-pivoted QR instead holds it between 31.72 and 43.09 at every stop. The two forward errors, measured against a BigInt answer, follow their own condition numbers: 0.3173 against 8.07·10⁻¹⁶ at γ = 10¹⁰. κ(A) is 10⁶ throughout and never moves — nothing about the problem's difficulty changes across this axis, only which columns a one-line rule happens to pick.024681010⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹10¹⁵10¹⁹log₁₀ γ — conditioning of the first m columnscondition number, and relative errorκ(ZᵀHZ), naiveκ(ZᵀHZ), pivotederror, naiveerror, orthonormalone problem, two rulesκ(ZᵀHZ) naive, γ = 10¹⁰4.9·10¹⁸κ(ZᵀHZ) pivoted, worst43error, naive0.32error, pivoted8.1·10⁻¹⁶κ(A) does not move across this axisand the answer moves by fourteen digits
Fig. 30 And the other square this collection measures, from the orthogonality field.

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.

Active setBackward errorBarrier parameterComponentwise condition numberCondition numberExact ground truthForward errorInterior point methodNormal equationsSaddle-point systems