A condition number sent to infinity

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.

Where the infinity comes from #

Minimise ½xᵀHx − fᵀx subject to Cx ≥ d. Add slacks s = Cx − d, 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.

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.

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.

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.

Neither number in that finding belongs to the particular active set, and the way to establish that is to move it. The active set is what the method is discovering, so a claim made at one of them is a claim about one iterate of one problem.

Across active sets of 1, 2, 3, 4 and 5 constraints, at μ = 10⁻¹⁴, the augmented error reads 1.35, 5.02, 1.04, 3.49 and 4.81 ·10⁻¹⁵ while its condition number reads 1.44·10¹⁶, 1.40·10¹⁶, 3.04·10¹⁵, 2.19·10¹⁵ and 2.16·10¹⁵. The condition number moves by a factor of seven across the family and the error does not follow it anywhere — it stays inside one decade, at the unit roundoff, at every one.

The condensed column over the same five reads 1, 1, 0.310, 0.0665 and 0.288. Not one of them has a correct digit, and the two where the active set is smallest have relative errors of exactly one, which is what a solver returns when the answer it produces is no closer to the truth than zero is.

The same finding at two other sizes #

The family has three parameters and the sweep above moves one. The other two are the size of the quadratic and the number of constraints, and both change the shape of the augmented matrix rather than only its diagonal — which is the case where a structural explanation could plausibly stop holding.

So the finding is a property of the two forms and not of the example. Over seven members of the family — five active sets, one larger quadratic, one larger constraint set — the augmented error runs from 2.17·10⁻¹⁶ to 5.02·10⁻¹⁵, a factor of 23, while the augmented condition number runs from 1.18·10¹⁵ to 1.44·10¹⁶ and the condensed error runs from 0.0665 to 2.19. The two forms are not close together at one setting and apart at another; they are twelve to fifteen orders apart everywhere the generator will draw.

How wrong the bound is, in both directions #

“The condensed error follows κ·u exactly” and “the augmented form is twelve orders under it” are two verdicts, and both become more useful as numbers. Divide each measured error by κ₂·u across the whole sweep:

the condensed ratio runs 0.21, 0.13, 0.034, 0.041, 0.11, 0.039, 0.036, 0.12 at μ from 1 to 10⁻¹⁴. The augmented ratio runs 0.24, 8.3·10⁻³, 3.2·10⁻⁵, 3.3·10⁻⁷, 3.6·10⁻⁹, 5.2·10⁻¹¹, 2.8·10⁻¹³ and 3.1·10⁻¹⁵.

For the condensed form the bound is not merely correct, it is right — within a factor of six over fourteen orders of magnitude, which is as well as a worst-case bound ever behaves. Anyone quoting κ·u about the condensed system is not being conservative; they are estimating, accurately.

For the augmented form the same bound over-predicts by a factor that grows exactly as 1/μ², because its numerator is a constant and the bound is not, reaching fifteen orders of magnitude at the end of the sweep.

Which is a sharper statement than “the bound is loose” #

One question the pair raises and this page does not answer: whether the condensed form can be repaired rather than replaced. The two essays before this one in the elimination and least-squares fields both find that iterative refinement rescues a method whose trouble is its route rather than its question — an explicit inverse gets its backward error back, and a Sherman–Morrison update gets everything back. The condensed form’s trouble is squarely a route: the augmented system it is an elimination of has a componentwise condition number of 13.3 at every μ, so the question is well posed and only the form is not.

That makes refinement against the augmented residual the obvious thing to try — solve with the Cholesky, form the residual of the (n + p) × (n + p) system, correct. It would keep the condensed form’s smaller factorisation and its fixed ordering, which are the real reasons codes prefer it, and it is one matrix–vector product a step. Whether the correction converges is not obvious and is not measured here: the corrector is as ill-conditioned as the original solve, so the contraction factor is whatever κ(M)·u is, and at μ = 10⁻¹⁴ that is above one. Somewhere between μ = 10⁻⁸ and 10⁻¹² it would stop converging, which is a boundary worth having and is a phase’s work rather than a paragraph’s.

A bound being loose is ordinary and this collection has measured it a dozen times — an expansion’s rank bound loose by six, a growth bound loose by a square root, a conjugate gradient bound loose by twelve. In every one of those the looseness is a property of the bound, and knowing the factor is enough to use it.

Here it is not a property of the bound. It is tight on one member of a pair and fifteen orders loose on the other, on the same problem, at the same μ, with condition numbers within a factor of eight of each other — 3.04·10¹⁵ against 2.40·10¹⁶ at μ = 10⁻¹⁴.

So κ₂ cannot be calibrated here, even in principle. A calibration is a factor applied to a bound to turn it into an estimate, and the factor would have to be 0.1 for one form and 10⁻¹⁵ for the other, with nothing in κ₂ to say which is in hand. The number does not merely fail to predict the error; it fails to distinguish the two cases at all, which is the stronger failure and is why the repair is a different condition number rather than a corrected version of this one.

That is also why the pair had to be measured together. The condensed row alone would show a bound behaving beautifully and confirm every expectation. The augmented row alone would show a bound loose by fifteen orders, which reads as a badly chosen example. Only the two side by side say that the quantity is not describing the problem, because on this page the problem is held fixed and only the form changes.

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.

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.

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.

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Δx − r₂). 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.

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.

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.

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.