The reference was a method

A constraint is a weight at infinity measures the method of weighting against a reference, and says exactly how the reference is computed:

Form AᵀA and Aᵀb in floating point, assemble the saddle-point system [[AᵀA, Bᵀ], [B, 0]] with the stored doubles, convert every entry to the dyadic rational it already is, and solve in BigInt. That gives the exact solution of the problem the machine is holding, which is the right reference.

It is the right reference for that page, because the three solvers being compared all see the same rounded object and the comparison is about what each does with it. It is the wrong reference for this one, and noticing why is the whole of the page: the optimality conditions contain AᵀA, so assembling them is the road that squares the problem in a block disguise, and a comparison that starts after the assembly cannot see what the assembly cost.

The two routes #

The saddle-point route writes down the optimality conditions — the residual must be orthogonal to everything the columns can reach except in the directions the constraint spends, which is $A^{\mathsf T}Ax + B^{\mathsf T}\lambda = A^{\mathsf T}b$ alongside $Bx = d$ — assembles them as one $(n + m) \times (n + m)$ system, and eliminates. It is the derivation every textbook gives, it is where the multipliers come from, and it is exactly the object the constraint field is about.

The null-space route eliminates the constraint instead — the move two ways to remove a constraint prices in the constraint field, and the basis it works in is the object a basis that is the same subspace and not the same thing warns can be chosen badly. Factor $B^{\mathsf T} = Q[R; 0]$, so the last $n - m$ columns of Q span the null space of B. Write $x = Q_1 y_1 + Q_2 y_2$; the constraint fixes $y_1$ by a triangular solve with $R^{\mathsf T}$, and $y_2$ is an ordinary unconstrained least-squares problem in $AQ_2$. No cross-product matrix is formed anywhere.

Both are exact in exact arithmetic. Both cost O(mn²) plus a least-squares solve. And the conditioning each of them hands to the arithmetic is completely different.

What each of them loses #

At κ(A) = 10¹¹, on a problem with six unknowns, twelve data rows and two constraints:

The first two lines are the expected ones and the fourth is the finding.

Nearly all of the saddle-point route’s loss happens before any elimination begins. Solving the assembled system exactly, in BigInt rationals, with no rounding of any kind after the assembly, recovers almost nothing: 3.45·10⁻⁴ against the floating-point solve’s 4.64·10⁻⁴. The elimination is responsible for a factor of 1.3 and the assembly for everything else.

That is a different statement from the saddle-point system is ill-conditioned, and it is a stronger one. An ill-conditioned system is one a stable solver handles as well as anything can; here the information is gone by the time the solver is called, and no solver recovers it. The entries of AᵀA are sums of products of numbers spanning eleven decades, each rounded to fifty-three bits, and the small singular directions of A are represented in that sum by digits the sum does not have room for.

And the approximate method beats the exact one #

The third line of the table is worth its own sentence, because it inverts the usual hierarchy.

The method of weighting is an approximation: it replaces the constraint by a weight, its error is 1/τ², and the previous page on this ladder is entirely about how far that limit can be taken. The saddle-point route is exact: it writes down conditions the answer satisfies and solves them.

At κ(A) = 10¹¹ the approximation returns 3.03·10⁻⁹ and the exact route returns 4.64·10⁻⁴ — five orders of magnitude apart, in the approximation’s favour. The reason is the same one: weighting hands its problem to a Householder QR, which never forms a cross-product, so the $\tau^{-2}$ bias it deliberately introduces is far smaller than the rounding the other route’s assembly imposes by accident.

An approximation whose error is stated and bounded is a better object than an exact formulation whose arithmetic is not. That is the general form, and this collection has met it before — it is why the normal equations are the method no library uses despite being the shortest derivation, and the constrained case is the same sentence one block further out.

Why the reference had to change #

The measurement above cannot be made with the previous page’s reference, and it is worth showing why rather than asserting it.

That reference is the exact solution of the assembled system — which is to say, it is the fourth row of the table. Measured against it, the saddle-point route’s error is the elimination’s alone: a factor of 1.3, not five orders of magnitude. And the null-space route, which never forms AᵀA, would be measured against an answer that did, so its error would include the assembly’s rounding as though the route had committed it.

Both numbers would be wrong in the same direction and the comparison would come out backwards. That is not a hypothetical: it is what the first version of this measurement produced, with the two routes within a factor of three of each other at every κ(A) and the null-space route occasionally behind.

So the reference is part of the method being measured, and choosing it is a decision about what question is being asked. Against the assembled system: how well does each solver solve what it was given? Against the data: how well does each route answer the question? The first is what a backward-error analysis asks and the second is what a caller wants, and they are different by the whole of an assembly step.

What the exact reference is, here #

The reference is built by accumulating AᵀA and Aᵀb in exact rationals from the stored doubles of A and b — every entry of A is a dyadic rational exactly, the products and sums are exact, and the resulting saddle-point system is solved in BigInt.

So it is the exact answer to the problem the machine holds as data, rather than to the cross-product the machine would have formed from it. Nothing is rounded between the input and the answer. The distinction is the same one an exact answer to a measured problem draws one field over, where an exact solve of a system defined by data accurate to fourteen places is wrong in its first digit — there the gap is between the data and the world, and here it is between the data and a matrix assembled from it.

The conditioning that decides it #

One column of the figure is not an error and is the reason the errors are what they are.

κ(A) is 10¹¹ at the right-hand end and κ(AQ₂) is 4.9·10⁶ — the conditioning of what is actually solved once the constraint has been eliminated. The null-space route’s error is governed by that number and not by κ(A), which is why it is 1.89·10⁻⁹ rather than the 10⁻⁵ that κ(A)·u would suggest.

The saddle-point route has no such quantity. Its matrix contains AᵀA whole, weak directions and all, and the constraint appears as a border rather than as an elimination — so the difficulty the constraint could have removed is still in the matrix when it is assembled, and the assembly rounds it away rather than removing it cleanly.

That is the mechanism in one sentence, and it is the same ordering question the projection and the right angle settles for the unconstrained problem. One route removes the constraint and then squares nothing; the other squares everything and then imposes the constraint. The order of the two operations is the whole of the difference, and neither derivation mentions it.

Where the curve bends, and what it says about the rounding #

The shape of the two error curves is worth reading rather than just their endpoints, because it says which quantity each route is paying in.

The null-space route’s error runs 1.21·10⁻¹⁵, 6.68·10⁻¹⁵, 3.06·10⁻¹³, 4.56·10⁻¹², 4.01·10⁻¹¹ and 1.89·10⁻⁹ across κ(A) = 10, 10³, 10⁵, 10⁷, 10⁹ and 10¹¹ — six orders of error across ten of conditioning, which is a slope of about 0.6 rather than 1. It is not paying κ(A)·u; it is paying κ(AQ₂)·u, and κ(AQ₂) runs 4.5, 77, 1.2·10³, 1.9·10⁴, 3.1·10⁵ and 4.9·10⁶ — six orders, matching the error’s six.

The saddle-point route’s error runs 1.73·10⁻¹⁶, 1.61·10⁻¹³, 1.65·10⁻¹¹, 1.31·10⁻⁹, 7.08·10⁻⁷ and 4.64·10⁻⁴: twelve orders across ten of κ(A), a slope of about 1.2. That is the signature of a squared conditioning being only partly squared — the (1, 1) block is AᵀA and carries κ(A)², while the border carries κ(B) and the right-hand side is not all cross-product — and it is why the gap between the two routes widens with κ(A) rather than being a constant factor.

So the two curves are not the same curve at different heights. One tracks a conditioning the constraint has already reduced and the other tracks something close to its square, and the widening gap is the measurement that distinguishes those two explanations from a fixed penalty.

That is worth stating because a fixed penalty would be tolerable. A route that is always a hundred times worse is a route with a known price; a route whose penalty grows as κ(A)^0.6 has no price at all until the problem is in front of it, which is exactly the situation a small residual is not a small error describes from the other side.

What the two routes cost #

Accuracy is half a comparison and the other half is worth having, because the null-space route is not free.

The saddle-point route forms AᵀA at mn²/2 multiplications, borders it, and eliminates an (n + m) × (n + m) system at (n + m)³/3. For six unknowns, twelve rows and two constraints that is a few hundred operations and the assembly dominates at large m.

The null-space route factors $B^{\mathsf T}$ — an $n \times m$ Householder QR at $2nm^2 - 2m^3/3$ — applies Q to A at $2mn$ per column, and then solves a least-squares problem in an (rows) × (n − m) matrix. The last of those is the expensive part and it is an ordinary QR, so the route costs roughly what an unconstrained fit of the same size costs, plus the constraint’s own factorisation.

The two are within a small factor of each other whenever m is small, which is the case a constraint usually is: a curve required to pass through a point, a set of coefficients required to sum to one, a flow required to balance at a node. So the choice is not a trade at all in the regime that matters — one route is more accurate and costs about the same, and the derivation everybody writes down is the other one.

Where the saddle-point route genuinely wins is when there are many constraints, since n − m is then small and the null space is cheap to work in but the QR of $B^{\mathsf T}$ is not. That crossover is not measured here and it is a counting exercise rather than an arithmetic one.

What a code should actually do #

Three lines, and they follow from the table rather than from the derivation.

Use the null-space route when the constraint is small and exact. It is the most accurate thing measured here, it costs about what an unconstrained fit costs, and its error is governed by the conditioning of the problem that is left rather than by the conditioning of the fit.

Use the method of weighting when a least-squares routine is all there is, which is the case the previous page is about — and then put the weighted rows first and use a solver that does not form a cross-product, which is what that page’s whole second half is about. At τ = 10⁸ the measurements here put it within a factor of two of the null-space route at every κ(A) drawn.

Assemble the optimality conditions when the multipliers are wanted, and know that the assembly costs what the table says. If the multipliers are wanted and the fit is ill-conditioned, the answer is the augmented form rather than either of these.

None of that is available from the derivations, which are all correct and all silent about arithmetic. It is available from four numbers, and the four numbers need a reference that is outside every method being compared.

The refusals #

The claim worth breaking is the one that makes the saddle-point route sound safe: these are the exact optimality conditions, so an error in the answer is the elimination’s and a stable solver removes it. The assertion that solving the assembled system exactly recovers the answer is fed the pair at κ(A) = 10¹¹ — 4.64·10⁻⁴ in floating point, 3.45·10⁻⁴ in exact rationals — and required to fail.

The second refusal covers the reading that this is an argument against the formulation. Fed the measurements at κ(A) = 10¹ and 10³, the claim that the saddle-point route is unusable has to fail: 1.73·10⁻¹⁶ and 1.61·10⁻¹³, which is every digit the machine has.

And the third is a control on the reference, since the whole page turns on it. The exact answer computed from the data is fed the claim that it misses its own constraints, and has to refuse — a page comparing constrained solvers needs its yardstick to be feasible, and feasibility here is an equality rather than a tolerance.

What this does not say about the saddle-point system #

Three qualifications, because the constraint field is built on that system and none of this is an argument against it.

Not that the formulation is wrong. The optimality conditions are the definition of the answer, the multipliers are quantities a caller often wants, and the constraint field prices what it costs to solve them well. What is measured here is one particular way of assembling them — with AᵀA in the (1, 1) block — and that is a choice rather than a necessity.

Not that it cannot be assembled better. The augmented form $[[I, A], [A^{\mathsf T}, 0]]$ with the residual as an unknown has no cross-product in it at all, is larger, and is what a code that wants multipliers without squaring should use. It is not measured here and it is the obvious next comparison.

And not that κ(A) = 10¹¹ is a reasonable place to work. At κ(A) = 10³ the two routes are within a factor of thirty of each other and both are at 10⁻¹³, which is where most problems sit. The measurement is about which route degrades and how, not about a regime anybody should be in.

Still open: the augmented form, and the multipliers #

Two things this page leaves.

The augmented system. $[[I, A], [A^{\mathsf T}, 0]]$ for the unconstrained problem, and its bordered version for the constrained one, replace the cross-product with an identity block and a copy of A. The matrix is larger by the number of rows, it is symmetric indefinite rather than a saddle point with a definite corner, and its conditioning is κ(A) rather than κ(A)². Whether it recovers the null-space route’s accuracy at the null-space route’s cost, or sits somewhere between, is one sweep and is not run here.

The multipliers. Everything above measures x and nothing measures λ. The null-space route does not compute the multipliers at all — it eliminates them — so a caller who wants them has either to recover them by a second solve or to use a route that produces them. What that recovery costs in accuracy, and whether a multiplier recovered from an accurate x is better than one from a route that computed both badly, is the measurement that would decide which route a code with sensitivity analysis behind it should use.

One line #

The optimality conditions of a constrained fit contain a cross-product matrix, so assembling them loses the digits and solving them exactly does not bring the digits back.