A constrained system can be reduced by eliminating the multipliers or by eliminating the constrained directions. Both give the same answer in exact arithmetic and inherit different condition numbers — one of them squares the constraint's, and the other does not contain it at all.
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.
Fig. 1 Two condition numbers and two errors, over five decades of the constraint’s conditioning. The
system is the same system at every point on the axis.
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.
Fig. 2 What each route costs, counted rather than described. The crossover is not where the shapes of
the two recipes suggest.
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.
Fig. 3 The reduced problem, on a well-behaved constraint. The three bars per group are three ways of
building Z, which the next essay is about; here they agree.
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.
Fig. 4 With H the identity, so that the only conditioning anywhere in the picture belongs to the
constraint — and the separation is unchanged.
The range-space method solves with S = AH⁻¹Aᵀ. Its condition number satisfies
κ(S) ≤ κ(H) · κ(A)²
and it attains that: the measured κ(S) runs 13.0, 3.99·10⁴, 3.98·10⁸ as κ(A) goes 1, 10², 10⁴,
which is four orders of κ(S) per two decades of κ(A) — the square, to two digits, at every stop.
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.
Fig. 5 The other square in the same identity, from the next essay: κ(Z)² when Z is not orthonormal.
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):
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) once — linearly. What it does not inherit 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.
Fig. 6 The same distinction three essays on, where the elimination that squares is the one every code
performs and the difference reaches fourteen orders.
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.
Fig. 7 At one constraint, where the range-space method’s extra object is a scalar and its penalty is
the whole of what it costs.
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.
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.
Fig. 9 Why the third route can win on a sparse problem: the whole matrix, ordered three ways, with the
fill counted.
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.
Fig. 10 An objective flat in a direction, from the least-squares field. A constraint along that direction
makes the reduced problem definite and leaves H singular.
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.
Fig. 11 The same move in the field it was first priced in: forming a Gram matrix, and the condition
number doubling in exponent.
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.
Fig. 12 At κ(H) = 10⁴, where both lines lift together and the separation between them does not move —
the control that says the gap belongs to the constraint.
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.
Fig. 13 The fourth essay’s picture, which is what the second of those looks like when the preconditioner
declines to know anything about A.
Fig. 14 A mildly conditioned objective, where the null-space line is flat at 21 and the range-space one
has already crossed 10⁶ by the right-hand edge.Fig. 15 And at 10⁶, where κ(H) dominates both and the two errors converge on each other — which is the
case where the choice does not matter.Fig. 16 Eight constraints, where the range-space method’s m × m object is no longer negligible.Fig. 17 The next essay’s comparison at an intermediate stop.Fig. 18 The matrix both routes are eliminating, from the previous essay.Fig. 19 And the reason neither of them is a Cholesky of the whole thing.Fig. 20 A null space met in a different field, where it is the thing that makes an answer non-unique
rather than the thing that makes it computable.Fig. 21 The projection the least-squares field is built on, which is what x = x_p + Zv is doing.Fig. 22 The number that decides whether a computed Z is a basis at all.Fig. 23 What a Schur complement does to a sparsity pattern, from the sparsity field: dense, whatever it
came from.Fig. 24 And the route that eliminates nothing and preconditions instead, from the next essay but one.Fig. 25 The identity underneath the whole comparison, drawn in the field that owns it.Fig. 26 The basis question at an intermediate badness, from the essay that takes it up.Fig. 27 Two constraints rather than three, where the Schur complement is smaller and squares just as hard.Fig. 28 The third route to the inertia, which neither elimination needs.Fig. 29 And what factorising the whole matrix requires, from the sixth essay in this field.Fig. 30 The number that describes the error on the field’s hardest matrix.Fig. 31 And the basis the flat line was quietly choosing.