Two ways to remove a constraint

The previous essay left the matrix

K = [ H Aᵀ ] [ A 0 ]

with an inertia that is known, a Cholesky that stops at a predictable row and a spectrum in two pieces. None of that is a way of solving anything. This one is about the two ways that are, and about the fact that they are not two orderings of one computation.

Eliminate the multipliers #

The first block row says Hx + Aᵀy = f, so x = H⁻¹(f − Aᵀy), and substituting into the second gives

A H⁻¹ Aᵀ y = A H⁻¹ f − g

The matrix on the left is S = AH⁻¹Aᵀ, the Schur complement of H in K, and the previous essay already met it: it is positive definite whenever A has full row rank, and its first diagonal entry is the pivot Cholesky stops on. It is m × m, which is small when the constraints are few.

So the recipe is: factorise H once, solve m systems with it to build S, factorise S, solve for y, and back-substitute for x. This is the range-space method, named for the fact that the multiplier lives in the range of A. It is the obvious thing to do when m is a handful and n is large — one dense m × m factorisation on top of one factorisation of H.

Eliminate the constrained directions #

The second block row says Ax = g, which pins m of the n degrees of freedom and leaves n − m. Write x = x_p + Zv with Ax_p = g and AZ = 0: any x of that form is feasible, and every feasible x is of that form.

Substituting into the first row and multiplying by Zᵀ makes the multiplier term vanish, because ZᵀAᵀ = (AZ)ᵀ = 0. What is left is

(ZᵀHZ) v = Zᵀ(f − Hx_p)

an unconstrained system of size n − m with a matrix that is symmetric positive definite: vᵀ ZᵀHZ v = ‖H¹ᐟ²Zv‖² > 0 because Z has full column rank. So the null-space method turns a constrained problem into a smaller unconstrained one, on which every tool the rest of this site has applies — a Cholesky, conjugate gradients, a preconditioner.

ZᵀHZ is the reduced Hessian, and it is the object an optimisation code reports curvature about. Its eigenvalues are the second derivatives of the objective along the feasible directions, which is what “the problem is well conditioned on the manifold” means when anybody says it.

Both are correct, and that is the control #

Before any comparison is worth making, the two have to be shown to be solving the same problem. They are: the library runs both on one system, and the difference between the two answers is at the rounding level whenever nothing in the problem is ill conditioned — 3.9·10⁻¹⁶, 6.9·10⁻¹⁶ and 5.5·10⁻¹⁶ on the three bases the next essay compares.

That control is the thing that makes the rest of the page a measurement rather than a demonstration of two different algorithms. There is one answer. Both routes reach it in exact arithmetic. Everything below is about the arithmetic.

What each one inherits #

The range-space method solves with S = AH⁻¹Aᵀ. Its condition number satisfies

κ(S) ≤ κ(H) · κ(A)²

and the measured κ(S) runs 13.0, 3.99·10⁴, 3.98·10⁸ as κ(A) goes 1, 10², 10⁴ — four orders of κ(S) per two decades of κ(A), which is the square to two digits at every stop. What that squaring is a statement about is worth a section of its own, because it is not the bound.

The reason is the same reason the normal equations square the condition number, and it is the same algebra: AH⁻¹Aᵀ is (H⁻¹ᐟ²Aᵀ)ᵀ(H⁻¹ᐟ²Aᵀ), a Gram matrix, and a Gram matrix’s singular values are the squares of its factor’s. The site has priced that move once, in the least-squares field, where the conclusion was that nobody should form AᵀA. Here it is being formed on purpose, by a method that is otherwise the sensible one, because the alternative is a factorisation of something n × n.

The null-space method solves with ZᵀHZ. Its condition number satisfies

κ(ZᵀHZ) ≤ κ(H) · κ(Z)²

and A does not appear. With an orthonormal Z the bound is κ(H) and the measurement is flat: 21.13 at κ(A) = 1, 21.13 at 10², 21.13 at 10⁴. The reduced Hessian does not know how badly conditioned the constraint was.

The other factor in the bound, which turns out not to matter #

Both bounds carry κ(H), and everything above was measured with it held at 100. Turning that knob instead is the check the identity asks for, and it answers a question the κ(A) sweep cannot: of the two factors in κ(H)·κ(A)², which one is doing the damage?

Two more decades of it, which is the range over which a Hessian in an optimisation actually varies — early iterates are well scaled and late ones, near an active set, are not:

The reduced Hessian is better conditioned than the Hessian it came from, at every stop, and by a margin that widens:

That κ(ZᵀHZ) ≤ κ(H) is not luck, it is Cauchy interlacing. ZᵀHZ is a section of H onto a subspace, so its eigenvalues interlace H’s: the largest cannot exceed H’s largest and the smallest cannot fall below H’s smallest. The bound κ(ZᵀHZ) ≤ κ(H)·κ(Z)² is therefore loose in the same direction at every stop, and with an orthonormal Z the true statement is an inequality rather than the equality the bound’s shape suggests.

And neither error tracks κ(H) at all. Across six orders of magnitude of it the range-space error moves from 7.34·10⁻⁷ to 2.47·10⁻⁶ — one and a half orders, and not monotonically — and the null-space error from 2.03·10⁻¹² to 4.42·10⁻¹², a factor of two. Compare that with the κ(A) sweep, where the range-space error grows by six orders or more.

So the bound κ(S) ≤ κ(H)·κ(A)² has two factors and only one of them is a live quantity on this family. A reader who reads it as a warning about ill-conditioned Hessians has taken the wrong half: the Hessian’s conditioning is the factor that enters linearly and is then cut down by the projection, and the constraint’s is the factor that enters squared and is not cut down by anything.

The square is in the exponent, and the bound is a parallel line #

Two things are being said at once by “κ(S) ≤ κ(H)·κ(A)², and the measurement squares”. The first is that κ(S) is proportional to κ(A)². The second is that the constant of proportionality is κ(H). Extended a decade and read as a ratio, the first is exact and the second is out by a factor of twenty-five.

The third column settles at 3.98 and does not move again — a fitted slope of 2.000 in the log, across four decades. The exponent is not approximately two, it is two. The fourth column settles at 25.1 and does not move either, so the bound and the measurement are parallel lines a decade and a half apart, and 25.1 is κ(H)/3.98.

Which is the good case, and the reason to separate the two claims rather than let one carry the other. A bound whose slack is a constant is worth having: it predicts, it can be divided by 25 and used, and the prediction survives a change of conditioning. A bound that is attained is a different and stronger thing, and the range-space method never gets within an order of this one. The essay’s argument — that the range-space route squares what the null-space route does not — needs only the exponent, and the exponent is the part that is exact.

The κ(H)/4 in the slack is not mysterious either. κ(S) = κ(B)² for B = L⁻¹Aᵀ, and κ(B) sits anywhere between κ(A)/√κ(H) and κ(A)√κ(H); here it lands at almost exactly 2κ(A), which is nearly the middle of that interval on a log scale. So the bound is loose because it charges for the worst possible interaction between H and A, and a random H and a random A do not interact that way.

There is a practical reading of the constant, and it is the reason to measure a slack rather than quote a bound. A code deciding between the two routes wants to know how many digits the range-space method will cost on this problem, and the bound answers with κ(H)κ(A)², which on the system here overstates the loss by 1.4 decades — enough to reject a route that would have been fine. The measured relation κ(S) ≈ 4κ(A)² answers the same question and is right, and the only thing needed to obtain it was to run the sweep one decade further and divide.

The same sweep says the other half more strongly than the section above did. κ(ZᵀHZ) is 21.13 at every one of the six conditionings, unchanged in the fourth digit while κ(A) crosses five decades. That is the claim the whole null-space case rests on, and it is not merely flat, it is constant. assertTheSquareIsInTheExponentNotInTheBound holds the whole shape: the bound holds, the slack is constant and never under ten, the exponent is two to two decimal places, and the reduced Hessian does not move.

And the errors follow the condition numbers rather than the problem #

The two forward errors are measured against an answer computed in BigInt rationals from the same stored doubles, so neither of them is being compared against a better float. Over four decades of κ(A):

range-space 6.6·10⁻¹⁶ 1.4·10⁻¹² 3.0·10⁻⁸ null-space 8.4·10⁻¹⁶ 5.7·10⁻¹⁵ 8.8·10⁻¹³

The range-space error grows by a factor of 4.6·10⁷ and the null-space one by 1.1·10³. Seven orders against three, on one system, with one answer.

The null-space error is not zero and does not stay at the rounding level either, and saying why is the honest half. The particular solution x_p = A⁺g is computed from A, so it inherits κ(A). It does not inherit it linearly, though — fitted over five decades the slope is 0.863, and the sequence is not even monotone: the error falls from 8.36·10⁻¹⁶ at κ(A) = 1 to 4.16·10⁻¹⁶ at 10 before it starts rising. So the growth factor of 1.1·10³ quoted above is a ratio of two endpoints, and one of those endpoints is where the rounding happened to land on the easiest problem in the sweep. A slope over six points is the honest version and it says less than linear. What the null-space route does not inherit at all is the square. The claim the assertion makes is therefore three-part: the range-space error grows by six orders or more, the null-space error grows by fewer than five, and the two separate by at least three. All three are measured rather than bounded.

Which one to use, which is not settled by the above #

If the argument stopped here the answer would be “always the null-space method”, and it is not, for two reasons that the cost figure shows and the conditioning figure cannot.

Forming Z costs a factorisation of Aᵀ, and the orthonormal Z from a QR is dense: n × (n − m) entries, where A itself may have been sparse. For m small and n large that is the larger object in the problem, and the reduced Hessian ZᵀHZ is (n − m) × (n − m) — nearly as large as H, and dense whatever H was. The range-space method’s extra object is m × m.

And the reduced Hessian has to be formed to be used. ZᵀHZ costs n²(n − m) to build explicitly. A code that only needs matrix–vector products can apply it as three multiplications without ever assembling it, which is what a projected conjugate gradient does, and that route keeps the conditioning and drops the cost. It also gives up the ability to factorise, which is the trade the iterative field has made once already.

So the choice is a shape question — is m small or is n − m small — with a conditioning penalty attached to one branch, and the penalty is worth knowing before the shape decides.

Where the square comes from, in one line #

It is worth having the mechanism rather than the citation, because it is the same mechanism in three fields and the citation is different in each.

Write H = LLᵀ and put B = L⁻¹Aᵀ, which is n × m. Then S = AH⁻¹Aᵀ = BᵀB. The singular values of B are σᵢ(B), and the eigenvalues of BᵀB are σᵢ(B)². So

κ(S) = σ₁(B)² / σₘ(B)² = κ(B)²

and κ(B) is between κ(A)/√κ(H) and κ(A)√κ(H). Forming S is forming a Gram matrix, and a Gram matrix’s condition number is the square of its factor’s. There is no arithmetic in that derivation at all: the squaring happens in exact arithmetic, and what floating point adds is only that the squared quantity is the one the solve is conditioned on.

The same three lines are the whole of why the normal equations lose twice as many digits as a QR, why the Gram matrix in a sketching argument has to be handled by a factorisation rather than formed, and why CholeskyQR needs a second pass. Four fields, one identity, and this is the fourth time the site has met it.

The third route, which is neither #

There is a way of solving Kz = b that eliminates nothing: factorise the whole matrix, with a symmetric indefinite factorisation that allows 2 × 2 pivots. The essay that introduced it is about the pivot rule; what matters here is that it inherits κ(K) and nothing squared, so it is the accurate route as well as the general one.

Its cost is a factorisation of an (n + m) × (n + m) matrix, which is more arithmetic than either elimination and less than either elimination plus the object it forms — and for a sparse K it can be very much less, because the fill of the whole matrix under a good ordering can be smaller than the fill of S, which is dense whenever any two constraints share a variable. The essay on what an ordering can be chosen for takes that up, and it needs one more ingredient first.

The other thing the range-space method needs, and does not have #

There is a hypothesis in the range-space derivation that is easy to walk past: H must be invertible. The null-space method does not need it — ZᵀHZ can be positive definite while H itself is only positive semidefinite, and in a great many real problems it is.

A least-squares fit with more unknowns than data has a singular AᵀA. A structure with a mechanism has a singular stiffness matrix until the supports are applied. An optimisation problem whose objective is flat in some direction has a singular Hessian in that direction, and the constraint is precisely what removes the flatness. In every one of those the whole system K is nonsingular, the problem has a unique answer, and the range-space method cannot start.

That is a genuine asymmetry rather than a technicality, and it is the reason optimisation codes lean towards the null-space route while flow codes lean towards the Schur one: in a flow problem the (1, 1) block is a discretised diffusion and is definite, and in an optimisation problem it is whatever the objective’s curvature happens to be. The measurement above compares the two where both are legal; the choice in practice is often made where only one is.

What the measurement says about the identity #

The site’s spine is forward error ⪅ condition number × backward error, and this page is an unusually clean instance of it. Both methods are backward stable in the ordinary sense: each returns the exact answer to a nearby version of the system it actually solved. The trouble is that the systems they actually solve are different, and one of them was manufactured with a condition number that is the square of anything in the original problem.

So the forward error is large for the reason the identity says, and the blame is not with the arithmetic and not with the problem. It belongs to a step of the method — a step whose whole purpose was to make the problem smaller. That is a third author, and the site’s usual pair does not have a slot for it.

The nearest thing already written is the least-squares field’s verdict on the normal equations, and the shape is identical: a reformulation that is algebraically exact and numerically a choice. What is new here is that both reformulations are exact, both are standard, and which one squares depends on which block is eliminated.

The refusal #

The reading this page has to close is the one that sounds like caution rather than carelessness: two block eliminations of the same nonsingular system are two orderings of the same arithmetic, so they lose the same accuracy. It is what “both are backward stable” would mean if backward stability were a property of the answer rather than of the system solved.

The assertion is fed the pair at κ(A) = 10⁴ — 3.0·10⁻⁸ and 8.8·10⁻¹³ — and required to reject the claim that they agree within two orders. It does, by four.

The same file’s other two refusals guard the parts of the argument that are easiest to over-read. One is fed the claim that AH⁻¹Aᵀ is indefinite and required to refuse it, because the whole range-space method rests on that matrix being factorisable by a Cholesky. The other is fed a Z whose product with A is not zero and required to refuse it as a null-space basis, because everything in the second half of the page assumes the term ZᵀAᵀy vanishes exactly and not approximately.

What is left over #

Two things this page does not settle, and both are taken up later in the field.

The first is that the null-space method’s advantage was measured with one particular Z — the orthonormal basis from a QR of Aᵀ, which has κ(Z) = 1 by construction. Nothing in the derivation says the basis has to be that one, and the identity κ(ZᵀHZ) ≤ κ(H)κ(Z)² has a second square in it that the flat line above was quietly setting to one. The next essay puts a different basis in and the flat line stops being flat.

The second is that neither route was given a preconditioner, and the whole question changes when one is available: a method that never eliminates anything, and instead moves the spectrum of K itself, is the third essay in this field. It has an unusual property for a preconditioner — its effect is a theorem rather than a measurement, and the theorem has the golden ratio in it.