A preconditioner that need not know the constraint

The previous essay’s preconditioner threw the off-diagonal blocks away and kept the two definite pieces. This one does the opposite: it keeps A exactly, keeps the zero, and replaces H by anything convenient.

P = [ G Aᵀ ] G symmetric, positive definite on the null space of A [ A 0 ]

It is called a constraint preconditioner for the obvious reason: it satisfies the constraint exactly, so a vector it returns is feasible whatever G was. And it has a spectrum with a property the previous one does not — the constraint is not in it.

The statement #

Let Z be any basis for the null space of A. Then P⁻¹K has

eigenvalue 1, with multiplicity 2m the n − m generalised eigenvalues of the pencil (ZᵀHZ, ZᵀGZ)

and that second list does not contain A at all. It does not contain Z either, in the sense that changing Z changes ZᵀHZ and ZᵀGZ by the same congruence and leaves the generalised eigenvalues alone.

The reason is the object the second essay in this field built. A constraint preconditioner is exactly the null-space method used as a preconditioner rather than as a solver: applying P⁻¹ projects onto the feasible set and solves a reduced problem with ZᵀGZ, and applying K measures the reduced problem with ZᵀHZ. So what the iteration sees is the reduced pencil, and the constraint was inverted exactly on the way in.

The 2m eigenvalues at one are what is left over: m from the multiplier block, and m from the directions the projection removes. They are exactly one because P and K agree exactly on those directions — the blocks that produce them are identical between the two matrices.

Why the constraint drops out, said without the pencil #

The pencil statement is the checkable one and it is not the intuitive one. Here is the same fact in a sentence about what the iteration does.

A Krylov method builds its iterate out of vectors P⁻¹K applied repeatedly to a residual. Applying P⁻¹ solves a saddle-point system with the same constraint block, so whatever comes out of it satisfies Ax = 0 exactly for the correction — the correction is feasible. Applying K then measures that feasible correction against the real objective. So every vector the iteration ever sees lives in the null space of A, and the operator it is effectively iterating with is the map that takes a feasible direction, measures its H-curvature, and returns the feasible direction with that G-curvature.

A is used twice per step and never approximated, so it never appears in the answer. It appears in the cost — a solve with P is a solve with a matrix as hard to factorise as K, structurally — and that is the trade this preconditioner makes: pay the constraint’s cost exactly rather than pay for approximating it and then pay again in iterations.

Two routes, and one of them never forms K #

The claim is checked by computing the spectrum twice, and the two computations have nothing in common.

The first forms P⁻¹K explicitly — one solve with P per column — and takes its eigenvalues with a real Schur decomposition, because the product of two symmetric matrices is not symmetric and there is nothing better available. That route is expensive, goes through a non-symmetric algorithm, and produces a list of n + m numbers.

The second builds an orthonormal Z from a QR of Aᵀ, forms ZᵀHZ and ZᵀGZ, and takes the generalised eigenvalues of that symmetric-definite pencil by a Cholesky and a Jacobi sweep — the route the pencil essays established. It is (n − m) × (n − m), it is symmetric throughout, and it never forms K, never forms P, and never multiplies one by the inverse of the other.

They agree to five digits relative at every conditioning tested. That is the site’s two-routes habit in the form it takes when one route is a theorem: the theorem says the second list is a sublist of the first, and the measurement is the agreement.

The measurement, and what it says #

Four condition numbers of A — 1, 10², 10⁴, 10⁶ — with H, G, the shapes and the seed all fixed. The six nontrivial eigenvalues come out as

0.00610895 0.0115277 0.0478812 0.237407 0.552778 1.35628

at every one of the four, to six significant figures, and equal to the generalised eigenvalues of a pencil that never saw A. Six decades of the constraint’s conditioning, and the part of the spectrum that governs the iteration does not move in the sixth digit.

Stated as a design consequence: a preconditioner for a constrained problem may decline to know anything about the constraint. The work of approximating goes entirely into G, which is an approximation to the objective, and the constraint is carried exactly rather than approximately. That is the opposite of what the shape of the matrix suggests — A is half of K’s off-diagonal and all of its structure, and it turns out to be the half that needs no attention.

And the invariance is not a property of the objective’s conditioning or of the problem’s shape either, which is what the other two knobs are for.

Across κ(H) = 1, 10², 10⁴ and 10⁶ the nontrivial spectrum’s movement over six decades of κ(A) reads 1.64·10⁻⁶, 3.02·10⁻⁶, 7.14·10⁻⁶ and 1.17·10⁻⁵. It rises with the objective’s conditioning — by a factor of seven over six decades — and never approaches anything a reader would call motion. The theorem says the pencil does not contain A; the measurement says the arithmetic of not containing it costs about κ(H)^0.15 of relative agreement.

So the invariance holds at five nontrivial eigenvalues, at six and at nine, with the residual movement running 3.84·10⁻⁹, 7.14·10⁻⁶ and 8.83·10⁻⁷ — which is not ordered by size and is ordered by nothing else either. That is the signature of an arithmetic residue rather than of a systematic dependence: the algebra says zero, and what is measured is whatever the two routes’ roundings left behind on that particular problem.

The one thing that does move #

The 2m eigenvalues the theorem places at exactly one come out at 1 ± 6.65·10⁻⁹ at κ(A) = 1 and 1 ± 6.04·10⁻⁵ at κ(A) = 10⁶. The drift rises by four orders as κ(A) rises by six, and it is the only quantity anywhere on the figure that knows about A.

It is the arithmetic rather than the algebra, and it is a useful reminder about what an exact multiplicity is. The theorem gives an eigenvalue of multiplicity 2m; a computation gives 2m numbers near one, and how near depends on how well the routine could invert a badly conditioned block. A code that clustered the computed eigenvalues by rounding them would report the right answer at κ(A) = 1 and the wrong number of distinct values at 10⁶.

The size of that drift is worth reading across the figures above, because it is the one number on the page that behaves like a computation rather than like a theorem. At κ(A) = 10⁶ it is 1.22·10⁻⁷ on the eight-variable problem, 6.04·10⁻⁵ on the ten-variable one and 1.13·10⁻⁴ at fourteen — three orders of magnitude of spread across a factor of under two in the size, on eigenvalues the theorem places at exactly one in every case. And it is not ordered by κ(H) at all: at κ(H) = 1 the drift reaches 3.3·10⁻⁴, five times what it reaches at κ(H) = 10⁶. Whatever governs it, it is not the conditioning of either block.

The refusal the library publishes covers exactly that: the assertion is fed the claim that the unit eigenvalues are exactly one in floating point and required to reject it. An exact multiplicity is not an exactly computed one, and the difference is measurable.

What a flat line is evidence of, and what it is not #

Six digits of agreement across six decades is a strong measurement and it is worth being precise about what it rules out.

It does not say that a badly conditioned constraint is harmless. K’s own condition number still contains σₘ(A), the system is still hard to solve accurately, and the range-space method on the same problem still loses seven orders. What it says is narrower and more useful: the iteration count under this preconditioner does not depend on κ(A) — for as long as there is an iteration count, which is a qualification the next section had to be written to earn.

The accuracy it converges to is a different quantity and does depend on κ(A), through the solve with P that every step performs. So the two halves separate cleanly — the constraint decides how accurately each step can be taken and the objective decides how many steps there are — and a figure of iteration counts and a figure of forward errors would look nothing alike.

This is the same separation the site’s spine always makes, in an unusual place: the condition number governs the accuracy, and something else entirely governs the work.

The count is flat because of which residual is being counted #

The sentence above is the kind this collection has learned to distrust: a claim about what an iteration does, argued from a spectrum, with no iteration run. So it was run. Left-preconditioned GMRES on K, with P applied by an exact solve, n = 16, m = 6, κ(H) = 10⁴, and κ(A) swept across eight decades. Two residuals are recorded at each step — the preconditioned one the method minimises and prints, and the true one, ‖b − Kx‖ ⁄ ‖b‖.

The left column is the flat line the theorem predicts, and for the first three rows the two columns agree and the claim is exactly right: fourteen steps at κ(A) = 1 and fourteen at 10⁴, with the constraint’s conditioning crossing four decades in between. That is the measurement the section above was asserting without having made.

The last two rows are the part the spectrum could not have shown. At κ(A) = 10⁶ the reported residual still crosses 10⁻⁹ at step fourteen and the true residual never gets closer than 3.7·10⁻⁶. At 10⁸ the reported residual crosses it at step seven — sooner than at any easier problem, which read on its own is a preconditioner that improves as the constraint degrades — while the true relative residual never leaves the third digit. The iteration has not converged faster. It has stopped converging, and the number being counted has quietly become a different number.

The reason is mechanical and has nothing to do with the algebra. Left preconditioning gives GMRES the system P⁻¹Kx = P⁻¹b, and the residual it reports is relative to ‖P⁻¹b‖. That norm grows with κ(A), because P inherits A and inverting a badly conditioned constraint block amplifies whatever it is handed. A fixed relative threshold on a growing denominator is crossed early. The count is honest about the quantity it counts and the quantity is no longer convergence.

The drift already measured is the other half of it. The 2m eigenvalues the algebra places exactly at one arrive at 1 ± 6·10⁻⁵ by κ(A) = 10⁶, and a Krylov method finishes early because a cluster is a single point to a polynomial of low degree. A cluster whose width is larger than the tolerance is not a cluster, and the minimal-polynomial argument that gives the step count stops applying at exactly the κ(A) where the true residual stalls.

So the clean separation stands, with a boundary on it. The constraint decides the accuracy each step can be taken to; the objective decides how many steps there are; and beyond about κ(A) = 10⁴ on this problem the first has eaten the second, because a step that cannot be taken accurately cannot be counted. assertTheFlatIterationCountIsTheReportedResidual holds all of it — that the true count is flat while it exists, that it ceases to exist at the two worst conditionings, that the reported count falls there, and that the forward error rises by seven orders across the eight. It is the same lesson as a stopping test that is a race, arriving from the side: a residual is a measurement of something, and which something is a question worth asking before counting steps of it.

One thing this does not indict, and the distinction matters for the section after next. It is a measurement of the operator P⁻¹K, which is the object whose spectrum this essay computes. The projected conjugate gradient the field ends on measures its residual in the reduced space, where the constraint has already been inverted exactly and never appears in a norm. Same algebra, different arithmetic — which is the more precise version of the claim that the two techniques are one technique.

And the price: the preconditioner is indefinite #

P has the same zero block as K and the same A, so by the inertia argument it has the same inertia: n positive eigenvalues and m negative. It is not positive definite, and the library’s assertion refuses the claim that it is.

That costs the method the iteration the previous essay used. MINRES needs a positive definite preconditioner, because the symmetric preconditioning L⁻¹KL⁻ᵀ needs a Cholesky of P to exist. With an indefinite P there is no such L, and preconditioned MINRES is not available.

What is available is the projected conjugate gradient: run CG on the reduced problem ZᵀHZ, preconditioned by ZᵀGZ, which the theorem above says is the same iteration. That works, it is what codes actually do, and it makes the constraint preconditioner exactly the null-space method with an inner iteration. Which is a satisfying place for the field to arrive at: the two apparently different techniques of the second essay and this one are the same technique, and the difference is whether the reduced problem is solved directly or approximately.

Which G, and the one constraint on it #

G has to be positive definite on the null space of A, not everywhere. That is a weaker requirement than it looks and it is what makes the method practical: a G that is indefinite in the constrained directions is fine, because the projection never sees those directions.

Three choices are usual. G = diag(H) is the cheapest thing that is usually definite enough. G = I gives a pencil (ZᵀHZ, ZᵀZ), which for an orthonormal Z is just the reduced Hessian’s spectrum — so the unpreconditioned reduced problem, and the iteration count is whatever κ(ZᵀHZ) says. And G = an incomplete factorisation of H is the version that behaves like a real preconditioner, with all the usual questions about how much fill to allow.

What none of them has to do is look at A. The figure’s G is a well-conditioned matrix that is not H and was not built from it, and the spectrum it produces is the same list at every κ(A) on the axis.

The Z that is not there #

Everything above was stated for “any basis Z”, and a reader who has met the next essay’s measurements will notice that this is a strong claim: that essay’s whole subject is that different bases for one null space have condition numbers decades apart and produce reduced Hessians whose condition numbers differ by their square.

Both are true, and the reconciliation is that a generalised eigenvalue of a pencil is invariant under congruence. Replacing Z by ZM for any nonsingular M sends ZᵀHZ to MᵀZᵀHZM and ZᵀGZ to MᵀZᵀGZM, and det(MᵀXM − λMᵀYM) = det(M)² det(X − λY): the characteristic polynomial is scaled and its roots are unchanged. Any two bases for the same null space differ by such an M.

So the spectrum is basis-free even though every matrix in it is not, and the badly conditioned basis of the next essay gives the same six numbers as the orthonormal one. What it does not give is the same computed six numbers, which is the arithmetic again — and it is why the routine that produces this figure uses the orthonormal basis and says so.

Where the two preconditioners meet #

Both of this field’s preconditioners have exact spectra and they are exact in different directions, which is worth putting side by side.

The block-diagonal one is exact when S = AH⁻¹Aᵀ is exact, and then the spectrum is three numbers that contain neither H nor A — the purest statement available, and unaffordable. Approximating S spreads all three.

The constraint one is exact whenever A is kept, whatever G is, and then the spectrum is 2m ones and a list that contains H and G and not A. It is affordable, because keeping A exactly costs a factorisation of P and that factorisation is a saddle-point solve of the same shape as the problem — which sounds circular and is not, because P’s (1, 1) block was chosen to be factorisable and K’s was not.

So the trade is legible: one preconditioner is exact in the objective and approximate in the constraint, the other is exact in the constraint and approximate in the objective, and the second is the one whose approximation the problem leaves to be chosen.

What this costs to build, honestly #

A constraint preconditioner needs P factorised, and P is a saddle-point matrix. That sounds like the method has replaced one hard problem by an identical hard problem, and the reason it has not is worth writing out, because it is the whole economic case.

K’s (1, 1) block is H, which in the problems this field is about is either large and dense — a Hessian assembled from an objective — or expensive to factorise for reasons of fill. P’s (1, 1) block is G, and G is chosen. Choose it diagonal and the factorisation of P collapses: eliminating the diagonal block leaves A G⁻¹ Aᵀ, an m × m matrix, and the whole solve is one small factorisation and two triple products. Choose it as an incomplete factorisation of H and P costs what that factorisation costs plus the same m × m object.

So the structure of the constraint is paid for exactly, once, and the part that was expensive is the part that was approximated. That is the reverse of the usual arrangement, where a preconditioner approximates everything a bit; and it is available only because the exact part is what the theorem needed kept.

The cost figure in the second essay counts the two eliminations; a constraint preconditioner sits between them, at the cost of the cheap one and the conditioning of the expensive one.