A condition number sent to infinity
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.
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.
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.
At other settings
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.
- A constraint is a weight at infinity — both name condition number, exact ground truth, normal equations, saddle-point systems
- A tensor that cannot be decomposed — both name backward error, condition number, exact ground truth, forward error
- The units the matrix is measured in — both name componentwise condition number, condition number, exact ground truth, forward error
- A small residual is not a small error — both name backward error, condition number, forward error
- An accuracy that is a backward error — both name backward error, condition number, forward error
- An eigenvalue count that cannot be slightly wrong — both name condition number, exact ground truth, saddle-point systems
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