Numerical Linear Algebra
where the algebra stops being the arithmetic
The matrix a constraint makes

A preconditioner that need not know the constraint

Keep the constraint block exactly and replace the objective block by anything positive definite on the null space. The preconditioned matrix then has 2m eigenvalues at exactly one, and its remaining n − m are the generalised eigenvalues of a pencil in which the constraint does not appear. Sweep its condition number over six decades and they do not move in six digits.
39 min read 30 figures Exact ground truthTwo routes to a numberWhere the algebra stops being the arithmetic

Worth reading first: Three eigenvalues, and two are the golden ratio · The basis nobody chose on purpose · The zero that is not a missing entry.

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 constraint preconditioner's spectrum, over six decades of the constraint's condition numberP = [[G, Aᵀ], [A, 0]] with G an approximation to H — here a well-conditioned matrix that is not H at all. The theorem says P⁻¹K has eigenvalue 1 with multiplicity 2m = 8 and n − m = 6 others, which are the generalised eigenvalues of the pencil (ZᵀHZ, ZᵀGZ) for any basis Z of the null space of A. A does not appear in that list, and the measurement is the flat lines: κ(A) crosses six decades along the horizontal axis and the 6 eigenvalues move by 7.14·10⁻⁶ relative, which is the arithmetic. A preconditioner for a constrained problem can decline to know anything about the constraint, because the constraint has already been inverted exactly inside it. The one quantity that does move is the drift of the eigenvalues the theorem puts at exactly one: 6.65·10⁻⁹ at κ(A) = 1 and 6.04·10⁻⁵ at 10⁶, which is κ(A) times the unit roundoff.012345610⁻³10⁻²10⁻¹110¹log₁₀ κ(A)eigenvalue of P⁻¹Kthe constraint is not in itnontrivial6at one8movement, six decades7.1·10⁻⁶drift at one6·10⁻⁵2m at onethe marks are a pencil that never saw Aand the lines are the preconditioned matrix
Fig. 1 Six eigenvalues, four condition numbers of A, and six digits of agreement. The horizontal axis is a property of the problem that the vertical axis does not contain.

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.

Three bases for the same null space, at 10 unknowns and 4 constraintsThe same constrained problem solved three times, differing only in which basis Z is used for the null space of A. The orthonormal basis, from a QR of Aᵀ, has κ(Z) = 1 exactly and is dense — 100 per cent of its entries are nonzero. The fundamental basis built on the first 4 columns has κ(Z) = 4.436, and its reduced Hessian comes out at 83.77, which is κ(Z)² to within a factor of 4.26 — the square attained rather than bounded. Its answer is wrong by 1.26·10⁻¹⁵, against 1.29·10⁻¹⁵ for the orthonormal one. Choosing the same kind of basis by pivoting instead gives κ(Z) = 2.258, an error of 1.41·10⁻¹⁵, and the same 50 per cent density: all of the sparsity and none of the loss.κ(A) = 10⁶ throughout · κ(H) = 100 · the answer is the same answer for every basisorthonormal — κ(Z)1κ(ZᵀHZ)32relative error1.29·10⁻¹⁵first m basic — κ(Z)4.44κ(ZᵀHZ)83.8relative error1.26·10⁻¹⁵pivoted basic — κ(Z)2.26κ(ZᵀHZ)43.1relative error1.41·10⁻¹⁵what the choice costsdensity, orthonormal1density, fundamental0.5κ(ZᵀHZ) ÷ κ(Z)², naive4.3error, pivoted choice1.4·10⁻¹⁵every one of them is a basisand one of them loses fourteen digits
Fig. 2 The reduced problem the pencil is built on, from the next essay. All three bars in each group give the same generalised eigenvalues.

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 penalties on the same problem, against the offset in the signalBest relative error for each penalty at four offsets. With no offset the two are within 7% of each other. At an offset of 10 the derivative penalty is 1.85 times better, because a constant lies in its null space and costs it nothing, while the norm penalty pays for the whole offset at every λ.best relative error at each offset‖x‖, offset 00.2035‖L₁x‖, offset 00.1901‖x‖, offset 20.0572‖L₁x‖, offset 20.0524‖x‖, offset 50.0333‖L₁x‖, offset 50.0245‖x‖, offset 100.0240‖L₁x‖, offset 100.0129what the null space buysadvantage at offset 01.1advantage at offset 21.1advantage at offset 51.4advantage at offset 101.9the norm penalty pays for a constantthe derivative penalty does not
Fig. 3 The subspace every vector in the iteration lives in, drawn in the field where it first appeared as an obstacle rather than as a tool.

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 constraint preconditioner's spectrum, over six decades of the constraint's condition numberP = [[G, Aᵀ], [A, 0]] with G an approximation to H — here a well-conditioned matrix that is not H at all. The theorem says P⁻¹K has eigenvalue 1 with multiplicity 2m = 8 and n − m = 6 others, which are the generalised eigenvalues of the pencil (ZᵀHZ, ZᵀGZ) for any basis Z of the null space of A. A does not appear in that list, and the measurement is the flat lines: κ(A) crosses six decades along the horizontal axis and the 6 eigenvalues move by 1.64·10⁻⁶ relative, which is the arithmetic. A preconditioner for a constrained problem can decline to know anything about the constraint, because the constraint has already been inverted exactly inside it. The one quantity that does move is the drift of the eigenvalues the theorem puts at exactly one: 2.57·10⁻⁹ at κ(A) = 1 and 3.3·10⁻⁴ at 10⁶, which is κ(A) times the unit roundoff.012345610⁻¹110¹10²log₁₀ κ(A)eigenvalue of P⁻¹Kthe constraint is not in itnontrivial6at one8movement, six decades1.6·10⁻⁶drift at one3.3·10⁻⁴2m at onethe marks are a pencil that never saw Aand the lines are the preconditioned matrix
Fig. 4 With H the identity, where the pencil is (I, ZᵀGZ) and its eigenvalues are the reciprocal curvatures of G on the null space.

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.

The constraint preconditioner's spectrum, over six decades of the constraint's condition numberP = [[G, Aᵀ], [A, 0]] with G an approximation to H — here a well-conditioned matrix that is not H at all. The theorem says P⁻¹K has eigenvalue 1 with multiplicity 2m = 8 and n − m = 6 others, which are the generalised eigenvalues of the pencil (ZᵀHZ, ZᵀGZ) for any basis Z of the null space of A. A does not appear in that list, and the measurement is the flat lines: κ(A) crosses six decades along the horizontal axis and the 6 eigenvalues move by 1.17·10⁻⁵ relative, which is the arithmetic. A preconditioner for a constrained problem can decline to know anything about the constraint, because the constraint has already been inverted exactly inside it. The one quantity that does move is the drift of the eigenvalues the theorem puts at exactly one: 5.85·10⁻⁹ at κ(A) = 1 and 1.94·10⁻⁵ at 10⁶, which is κ(A) times the unit roundoff.012345610⁻⁴10⁻³10⁻²10⁻¹110¹log₁₀ κ(A)eigenvalue of P⁻¹Kthe constraint is not in itnontrivial6at one8movement, six decades1.2·10⁻⁵drift at one1.9·10⁻⁵2m at onethe marks are a pencil that never saw Aand the lines are the preconditioned matrix
Fig. 5 At κ(H) = 10⁶, where the spectrum does spread — because the axis it depends on is the objective’s conditioning and not the constraint’s.

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 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.

The share of shifts inside a pair of eigenvalues that count it wrongly, against the pair's separationTwo eigenvalues at 1 and 1 + gap, with six others spread around them, and 100 shifts placed strictly between the pair — where the count must read 3. Down to a separation of 10⁻¹² every shift reads it correctly. At 10⁻¹³ one of 100 does not, at 10⁻¹⁴ 19 do not, and at 10⁻¹⁵ none of them reads it correctly. The boundary sits where n‖A‖u puts it — -4.405 for this matrix — because the floating-point count is the exact count of a matrix within that distance of A. This is the only place in the method where the answer can be wrong**, and it is wrong by a whole eigenvalue when it is: the failure is a miscount, not a small error.-15-13-11-9-7-5-300.250.50.751log₁₀ separation of the pairshare of shifts counted wronglyn‖A‖uthe only place it failswrong at 10⁻¹²0wrong at 10⁻¹⁴19wrong at 10⁻¹⁵100n‖A‖u-4.4wrong by a whole eigenvalueor not wrong at all
Fig. 6 The same distinction in the spectra field’s version, where a count rather than a value is what survives the arithmetic.

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). The method converges in the same number of steps at κ(A) = 1 and at 10⁶.

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 constraint preconditioner's spectrum, over six decades of the constraint's condition numberP = [[G, Aᵀ], [A, 0]] with G an approximation to H — here a well-conditioned matrix that is not H at all. The theorem says P⁻¹K has eigenvalue 1 with multiplicity 2m = 8 and n − m = 6 others, which are the generalised eigenvalues of the pencil (ZᵀHZ, ZᵀGZ) for any basis Z of the null space of A. A does not appear in that list, and the measurement is the flat lines: κ(A) crosses six decades along the horizontal axis and the 6 eigenvalues move by 7.14·10⁻⁶ relative, which is the arithmetic. A preconditioner for a constrained problem can decline to know anything about the constraint, because the constraint has already been inverted exactly inside it. The one quantity that does move is the drift of the eigenvalues the theorem puts at exactly one: 6.65·10⁻⁹ at κ(A) = 1 and 6.04·10⁻⁵ at 10⁶, which is κ(A) times the unit roundoff.012345610⁻³10⁻²10⁻¹110¹log₁₀ κ(A)eigenvalue of P⁻¹Kthe constraint is not in itnontrivial6at one8movement, six decades7.1·10⁻⁶drift at one6·10⁻⁵2m at onethe marks are a pencil that never saw Aand the lines are the preconditioned matrix
Fig. 7 The claim again, at the hero’s conditioning, with the drift of the unit eigenvalues printed — which is the half that does depend on κ(A).

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.

MINRES on a saddle-point system under three Schur-complement approximations, and under noneThe same system at 12 unknowns and 5 constraints, solved four ways. With the exact Schur complement the preconditioned matrix has three distinct eigenvalues and the residual falls to 2.947·10⁻¹⁴ in three steps, after which nothing is left to remove. Replacing S by A diag(H)⁻¹Aᵀ costs 13 steps and replacing it by a scaled AAᵀ costs 13; the unpreconditioned system takes 25. The exact preconditioner is unaffordable — forming S costs 5 solves with H and a decomposition — so its value is as the statement the cheap ones are measured against, and the measurement needs no reference solution: the distance from {1 − φ, 1, φ} is a property of the approximation alone.0369121518212410⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1MINRES steprelative residualthree stepssteps to a residual of 10⁻¹¹exact S3diag(H)13scaled AAᵀ13none25three eigenvalues, three stepsand the approximations pay for the difference
Fig. 8 The iteration the previous essay’s definite preconditioner allows, for contrast with the one this essay’s does not.

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.

What a hierarchical preconditioner costs, against how accurate it is, at κ = 1.75e+2The iteration is conjugate gradients on the true matrix and the preconditioner is a solve with the representation, so what is measured is how much of the compression the outer loop repairs and what the repair costs. The count falls monotonically — 16, 14, 7, 4, 3, 2 steps against 72 with no preconditioner — and it keeps falling all the way to the tightest accuracy, which is the number a solver prints and the one everybody watches. The total work, one dense product with A plus one solve with A_H per step, is the other curve, and its minimum is at ε = 0.5 where the preconditioner has rank 1. Past that point every extra digit in the preconditioner is bought and not needed: at the tightest setting the run costs 1.59 times the cheapest, having taken 14 fewer steps to get there.-9-7-5-3-1110⁵10⁶10⁷log₁₀ ε of the preconditionermultiplications for the whole solveno preconditioner: 72 stepsiteration count, rescaledmultiplications spenttwo curves, opposite waysκ175steps, no preconditioner72cheapest ε0.5its rank1tightest ⁄ cheapest1.6the count is what is printedand the cost is what is spent
Fig. 9 The usual question about a G, asked in the field where it was first asked: how much of a factorisation to keep.

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.

Three bases for the same null space, at 10 unknowns and 4 constraintsThe same constrained problem solved three times, differing only in which basis Z is used for the null space of A. The orthonormal basis, from a QR of Aᵀ, has κ(Z) = 1 exactly and is dense — 100 per cent of its entries are nonzero. The fundamental basis built on the first 4 columns has κ(Z) = 1.994·10⁸, and its reduced Hessian comes out at 3.801·10¹⁶, which is κ(Z)² to within a factor of 0.956 — the square attained rather than bounded. Its answer is wrong by 0.05179, against 1.07·10⁻¹⁵ for the orthonormal one. Choosing the same kind of basis by pivoting instead gives κ(Z) = 2.06, an error of 6.71·10⁻¹⁶, and the same 50 per cent density: all of the sparsity and none of the loss.κ(A) = 10⁶ throughout · κ(H) = 100 · the answer is the same answer for every basisorthonormal — κ(Z)1κ(ZᵀHZ)25.6relative error1.07·10⁻¹⁵first m basic — κ(Z)1.99·10⁸κ(ZᵀHZ)3.8·10¹⁶relative error0.0518pivoted basic — κ(Z)2.06κ(ZᵀHZ)31.9relative error6.71·10⁻¹⁶what the choice costsdensity, orthonormal1density, fundamental0.5κ(ZᵀHZ) ÷ κ(Z)², naive0.96error, pivoted choice6.7·10⁻¹⁶every one of them is a basisand one of them loses fourteen digits
Fig. 10 The bases the invariance is invariant to, from the next essay: κ(Z) = 1, 2.06 and 1.99·10⁸ for the same null space.

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.

The spectrum of P⁻¹K with S ≈ A diag(H)⁻¹Aᵀ, at 10 unknowns and 4 constraintsP = blkdiag(H, S). With S the exact Schur complement AH⁻¹Aᵀ the preconditioned matrix has exactly three distinct eigenvalues — 1 with multiplicity n − m = 6, and (1 ± √5)/2 with multiplicity 4 each. Those are 1 − φ = -0.618034 and φ = 1.61803, the golden ratio, which arrives from λ² − λ − 1 = 0 rather than from anything anybody chose. The dashed lines are that closed form and the marks are the computed spectrum; here they are 9 distinct values and the largest distance from the closed form anywhere is 2.955. A minimal polynomial of degree three means a Krylov method finishes in three steps, which is what the next figure measures.-3-1.89435-0.7887010.3169491.42262.528253.63394.739550eigenvalue of P⁻¹Kwritten down, then computeddistinct9at 16φ computed4.6off the closed form31 − φ1φthe preconditioner's effect is a theoremand the golden ratio is in it
Fig. 11 What approximating the other half looks like: the three points of the previous essay, spread.

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.

Two-level convergence factor for three interpolations, at ε = 0.001Three bars per operator: the classical formula, the unconstrained energy minimiser, and the minimiser constrained to reproduce a constant. On the isotropic and aligned operators the constraint costs a factor of 4.8; on the rotated one it buys a factor of 1.33.lower is better; each bar is the factor the error falls by per two-level cycleisotropic, classical0.0609isotropic, minimiser0.0609isotropic, constrained0.0848aligned, classical0.0612aligned, minimiser0.0615aligned, constrained0.2933rotated 45°, classical0.3464rotated 45°, minimiser0.2843rotated 45°, constrained0.2131what the constraint is worthconstraint at isotropic1.4constraint at aligned4.8constraint at rotated 45°0.75the constraint costs where the method worksand buys where it does not
Fig. 12 The two routes counted, with the preconditioner’s arithmetic sitting on the cheaper line and its accuracy on the other.

At other settings

The constraint preconditioner's spectrum, over six decades of the constraint's condition numberP = [[G, Aᵀ], [A, 0]] with G an approximation to H — here a well-conditioned matrix that is not H at all. The theorem says P⁻¹K has eigenvalue 1 with multiplicity 2m = 8 and n − m = 6 others, which are the generalised eigenvalues of the pencil (ZᵀHZ, ZᵀGZ) for any basis Z of the null space of A. A does not appear in that list, and the measurement is the flat lines: κ(A) crosses six decades along the horizontal axis and the 6 eigenvalues move by 3.02·10⁻⁶ relative, which is the arithmetic. A preconditioner for a constrained problem can decline to know anything about the constraint, because the constraint has already been inverted exactly inside it. The one quantity that does move is the drift of the eigenvalues the theorem puts at exactly one: 1.93·10⁻⁸ at κ(A) = 1 and 3.11·10⁻⁵ at 10⁶, which is κ(A) times the unit roundoff.012345610⁻²10⁻¹110¹log₁₀ κ(A)eigenvalue of P⁻¹Kthe constraint is not in itnontrivial6at one8movement, six decades3·10⁻⁶drift at one3.1·10⁻⁵2m at onethe marks are a pencil that never saw Aand the lines are the preconditioned matrix
Fig. 13 At κ(H) = 100, where the six nontrivial eigenvalues span two orders instead of three.
The constraint preconditioner's spectrum, over six decades of the constraint's condition numberP = [[G, Aᵀ], [A, 0]] with G an approximation to H — here a well-conditioned matrix that is not H at all. The theorem says P⁻¹K has eigenvalue 1 with multiplicity 2m = 4 and n − m = 8 others, which are the generalised eigenvalues of the pencil (ZᵀHZ, ZᵀGZ) for any basis Z of the null space of A. A does not appear in that list, and the measurement is the flat lines: κ(A) crosses six decades along the horizontal axis and the 8 eigenvalues move by 5.01·10⁻⁷ relative, which is the arithmetic. A preconditioner for a constrained problem can decline to know anything about the constraint, because the constraint has already been inverted exactly inside it. The one quantity that does move is the drift of the eigenvalues the theorem puts at exactly one: 1.31·10⁻⁸ at κ(A) = 1 and 5.88·10⁻⁶ at 10⁶, which is κ(A) times the unit roundoff.012345610⁻³10⁻²10⁻¹110¹log₁₀ κ(A)eigenvalue of P⁻¹Kthe constraint is not in itnontrivial8at one4movement, six decades5·10⁻⁷drift at one5.9·10⁻⁶2m at onethe marks are a pencil that never saw Aand the lines are the preconditioned matrix
Fig. 14 Two constraints, so eight nontrivial eigenvalues and four at one.
The spectrum of P⁻¹K with S = AH⁻¹Aᵀ, exactly, at 10 unknowns and 4 constraintsP = blkdiag(H, S). With S the exact Schur complement AH⁻¹Aᵀ the preconditioned matrix has exactly three distinct eigenvalues — 1 with multiplicity n − m = 6, and (1 ± √5)/2 with multiplicity 4 each. Those are 1 − φ = -0.618034 and φ = 1.61803, the golden ratio, which arrives from λ² − λ − 1 = 0 rather than from anything anybody chose. The dashed lines are that closed form and the marks are the computed spectrum; here they are 3 distinct values and the largest distance from the closed form anywhere is 2.909·10⁻¹⁴. A minimal polynomial of degree three means a Krylov method finishes in three steps, which is what the next figure measures.00.3670080.7340171.101031.468031.835040eigenvalue of P⁻¹Kwritten down, then computeddistinct3at 16φ computed1.6off the closed form2.9·10⁻¹⁴1 − φ1φthe preconditioner's effect is a theoremand the golden ratio is in it
Fig. 15 The previous essay’s three points, for reading against this page’s flat lines.
The reduced Hessian and the answer, against how badly the basic columns were chosenγ is the condition number of the first 4 columns of A, and the naive rule calls exactly those columns basic. Its reduced Hessian's condition number climbs from 83.77 to 4.871·10¹⁸ across ten decades of γ — the square, because κ(ZᵀHZ) ≤ κ(H)κ(Z)² and the bound is attained. Choosing the basic columns by a column-pivoted QR instead holds it between 31.72 and 43.09 at every stop. The two forward errors, measured against a BigInt answer, follow their own condition numbers: 0.3173 against 8.07·10⁻¹⁶ at γ = 10¹⁰. κ(A) is 10⁶ throughout and never moves — nothing about the problem's difficulty changes across this axis, only which columns a one-line rule happens to pick.024681010⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹10¹⁵10¹⁹log₁₀ γ — conditioning of the first m columnscondition number, and relative errorκ(ZᵀHZ), naiveκ(ZᵀHZ), pivotederror, naiveerror, orthonormalone problem, two rulesκ(ZᵀHZ) naive, γ = 10¹⁰4.9·10¹⁸κ(ZᵀHZ) pivoted, worst43error, naive0.32error, pivoted8.1·10⁻¹⁶κ(A) does not move across this axisand the answer moves by fourteen digits
Fig. 16 The next essay’s sweep, which is what happens to ZᵀHZ when Z is not orthonormal — and which this page’s claim is invariant to.
What each way of eliminating a constraint inherits, over five decades of κ(A)The same system solved twice, at 8 unknowns and 3 constraints with κ(H) = 100. The range-space method forms S = AH⁻¹Aᵀ and inherits κ(S), which rises from 12.99 to 3.981·10¹⁰ — the square of κ(A), for the reason the normal equations square it. The null-space method solves with the reduced Hessian ZᵀHZ, whose condition number is 21.13 at the start of the sweep and 21.13 at the end: it does not contain κ(A) at all. The two forward errors, measured against a solution computed in BigInt rationals, follow their own condition numbers: 5.314·10⁻⁶ against 5.788·10⁻¹² at the far end. Both methods are correct and one of them is usable.01234510⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹log₁₀ κ(A)condition number, and relative errorκ(S)κ(ZᵀHZ)range-space errornull-space erroragainst a BigInt answerκ(S) at κ(A) = 10⁵4·10¹⁰κ(ZᵀHZ), all stops21range-space forward error5.3·10⁻⁶null-space forward error5.8·10⁻¹²both are the same algebraand only one squares
Fig. 17 The direct route the preconditioner is an inner iteration inside.
The 14 eigenvalues of a saddle-point matrix with 10 unknowns and 4 constraints, inside their closed-form bracketsK = [[H, Aᵀ], [A, 0]] with κ(H) = 100 and κ(A) = 100. The shaded bands are the Rusten–Winther brackets, computed from four numbers — the extreme eigenvalues of H and the extreme singular values of A — before the matrix was assembled: negative eigenvalues in [-0.995, -10·10⁻⁵] and positive ones in [0.01, 1.62]. The marks are the computed spectrum. There are exactly 10 above zero and 4 below, which is Sylvester's law of inertia and not a property of this matrix: K is congruent to blkdiag(H, −AH⁻¹Aᵀ), both blocks are definite, and congruence preserves signs. The two groups are separated by a gap containing zero, at a ratio of 39.3 between the innermost positive eigenvalue and the innermost negative one.-1-0.582271-0.1645420.2531860.6709151.088641.506370eigenvalue4 negative10 positivecounted before it was formedpositive10negative4at zero0innermost ratio39the zero block is a theoremand so is the count either side of it
Fig. 18 P’s inertia, which is K’s, and is why MINRES is unavailable here.
Relative eigenvalue error of two reductions of Ax = λBx, against κ(B), n = 12The pencil is built as B = LLᵀ and A = LΛLᵀ, so its generalised eigenvalues are exactly Λ and the error is a measurement rather than a comparison. Forming B⁻¹A gives a matrix whose departure from symmetry reaches 1.08 — about its own size — and reducing through the Cholesky factor gives one that is nearer symmetric by a factor of at least 6078. That difference is real and it does not appear in the accuracy: the two error curves have fitted slopes of 0.92 and 0.98 against κ(B) and stay within a factor of 2.3 of each other over fifteen decades. The conditioning belongs to the pencil, and no reduction of it escapes.10¹10⁴10⁷10¹⁰10¹³10¹⁶10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹κ(B)relative error, and asymmetryvia B⁻¹Avia Choleskyasymmetry of B⁻¹Au · κ(B)against a spectrum known exactlyslope, via B⁻¹A0.92slope, via Cholesky0.98worst ratio between them2.3asymmetry of B⁻¹A1.1the symmetry claim is trueand it is not about the accuracy
Fig. 19 Generalised eigenvalues taken two ways, in the field that owns the pencil.
What a perturbation of 10⁻⁶ moves, against the gap it is applied acrossThree curves against the gap between two eigenvalues, both axes logarithmic. Two are flat; the third rises steeply as the gap closes and then levels off.10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻⁸10⁻⁵10⁻²gap between the two eigenvalueshow far it movedthe eigenvectorsthe eigenvaluestheir plane‖E‖ / gapone perturbation, three answerseigenvalue shift, spread over the sweep1plane angle, spread over the sweep1eigenvector angle, spread1.6·10⁵the dashed line is Davis–Kahan's ‖E‖/gaptwo of the three never noticed
Fig. 20 What decides whether a computed eigenvalue of a pencil is worth its digits.
LSQR steps to a normal residual of 10⁻¹², with and without a sketched preconditionerTwo curves against κ(A). The unpreconditioned solve takes 15, 22, 32, 59, 150 steps as the conditioning worsens, and at the worst one it does not reach the tolerance inside 150. The preconditioned solve takes 10, 10, 10, 10, 10 — flat across eight decades. Both return the same vector, to the accuracy the problem allows: the preconditioned answer and the direct one differ by about κ(A)·u at every κ, 5.7·10⁻⁶ at the worst. The SAME sketch used as an answer instead of as a preconditioner is 4.5 away.10²10⁴10⁶10⁸10¹⁰110¹10²10³condition number of the matrixLSQR steps to 10⁻¹²no preconditionersketchedthe cost, not the answersteps, preconditioned10steps, without150answer gap, preconditioned5.7·10⁻⁶κ(A)·u at the worst κ1.1·10⁻⁶answer gap, sketch-and-solve4.5one sketch, two things to do with itand only one of them moves the answer
Fig. 21 Iteration counts under a preconditioner, counted the ordinary way.
Conjugate gradients at κ = 104, against the bound κ permitsA semi-logarithmic plot of the relative A-norm error against iteration count. The measured curve falls below a smooth dashed curve showing the classical condition-number bound.0408012016020024010⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹iteration‖e‖_A / ‖e₀‖_Ameasuredκ bound119 steps40×40, spectrum spread evenly in logbound permits 1417
Fig. 22 The iteration a projected conjugate gradient is, once the projection has been applied.
The augmented step's two condition numbers, and the one the error obeysκ₂ is a worst case over all perturbations of the same NORM, and a normwise perturbation is allowed to put its whole budget on the smallest entry of a matrix. The componentwise number — Skeel's ‖ |A⁻¹||A||x| ‖ / ‖x‖ — is a worst case over perturbations proportional to the entries, which is what a backward-stable factorisation actually makes. For the augmented interior-point matrix the two are 3.044·10¹³ and 13.25 at μ = 10⁻¹², a ratio of 2.296·10¹², and the measured error is 9.434·10⁻¹⁶ — which is the second number times the unit roundoff and has nothing to do with the first. The componentwise line is still finding its level down to about μ = 10⁻⁴, where the two groups of constraints have not yet separated, and from there it does not move in the fourth digit; where it settles — 13.25 here — is set by how many constraints are active rather than by μ. Eliminating the second block to reach the condensed form destroys the distinction: there the componentwise number is 1.324·10¹⁴, tracking the normwise one, because the large entries have been summed into CᵀDC and are no longer separately identifiable.-12-10-8-6-4-2010⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹10¹⁵log₁₀ μcondition number, and relative errorκ₂, augmentedcomponentwise, condensedcomponentwise, augmentederror, augmentedone matrix, two numbersκ₂ at μ = 10⁻¹²3·10¹³componentwise, same matrix13their ratio2.3·10¹²measured error9.4·10⁻¹⁶every library prints the top lineand the error obeys the third
Fig. 23 And the field’s central measurement, from the essay after next.
QᵀQ from classical Gram–Schmidt and from Householder on the 8×8 Hilbert matrixTwo eight-by-eight tables of QᵀQ. The upper one has ones on the diagonal and entries as large as one off it; the lower one is the identity to three decimal places everywhere.A = H8 · κ = 1.5·10¹⁰ · both factorisations reconstruct A to 6·10⁻¹⁷the diagonal is 1 in both — every column is a unit vector either way1.000000000001.000000000001.000000000001.000000000001.00000.002-0.002000001.0000.125-0.13300000.0020.1251.000-1.0000000-0.002-0.133-1.0001.000classical Gram–Schmidt1.000000000001.000000000001.000000000001.000000000001.000000000001.000000000001.000000000001.000Householderclassical ‖QᵀQ − I‖1.4Householder ‖QᵀQ − I‖1.4·10⁻¹⁵largest off-diagonal 1 against 3.1·10⁻¹⁶length is not angle
Fig. 24 The number that says whether the Z in the second route was a basis at all.
The constraint preconditioner's spectrum, over six decades of the constraint's condition numberP = [[G, Aᵀ], [A, 0]] with G an approximation to H — here a well-conditioned matrix that is not H at all. The theorem says P⁻¹K has eigenvalue 1 with multiplicity 2m = 6 and n − m = 7 others, which are the generalised eigenvalues of the pencil (ZᵀHZ, ZᵀGZ) for any basis Z of the null space of A. A does not appear in that list, and the measurement is the flat lines: κ(A) crosses six decades along the horizontal axis and the 7 eigenvalues move by 2.92·10⁻⁶ relative, which is the arithmetic. A preconditioner for a constrained problem can decline to know anything about the constraint, because the constraint has already been inverted exactly inside it. The one quantity that does move is the drift of the eigenvalues the theorem puts at exactly one: 2.73·10⁻⁸ at κ(A) = 1 and 3.14·10⁻⁵ at 10⁶, which is κ(A) times the unit roundoff.012345610⁻³10⁻²10⁻¹110¹log₁₀ κ(A)eigenvalue of P⁻¹Kthe constraint is not in itnontrivial7at one6movement, six decades2.9·10⁻⁶drift at one3.1·10⁻⁵2m at onethe marks are a pencil that never saw Aand the lines are the preconditioned matrix
Fig. 25 Three constraints, so seven nontrivial eigenvalues and six at one.
MINRES on a saddle-point system under three Schur-complement approximations, and under noneThe same system at 12 unknowns and 5 constraints, solved four ways. With the exact Schur complement the preconditioned matrix has three distinct eigenvalues and the residual falls to 6.534·10⁻¹⁶ in three steps, after which nothing is left to remove. Replacing S by A diag(H)⁻¹Aᵀ costs 3 steps and replacing it by a scaled AAᵀ costs 3; the unpreconditioned system takes 13. The exact preconditioner is unaffordable — forming S costs 5 solves with H and a decomposition — so its value is as the statement the cheap ones are measured against, and the measurement needs no reference solution: the distance from {1 − φ, 1, φ} is a property of the approximation alone.0246810121410⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1MINRES steprelative residualthree stepssteps to a residual of 10⁻¹¹exact S3diag(H)3scaled AAᵀ3none13three eigenvalues, three stepsand the approximations pay for the difference
Fig. 26 The other preconditioner’s control case.
The reduced Hessian and the answer, against how badly the basic columns were chosenγ is the condition number of the first 2 columns of A, and the naive rule calls exactly those columns basic. Its reduced Hessian's condition number climbs from 148.6 to 2.711·10¹⁹ across ten decades of γ — the square, because κ(ZᵀHZ) ≤ κ(H)κ(Z)² and the bound is attained. Choosing the basic columns by a column-pivoted QR instead holds it between 64.71 and 75.33 at every stop. The two forward errors, measured against a BigInt answer, follow their own condition numbers: 0.8965 against 2.01·10⁻¹⁵ at γ = 10¹⁰. κ(A) is 10⁶ throughout and never moves — nothing about the problem's difficulty changes across this axis, only which columns a one-line rule happens to pick.024681010⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹10¹⁵10¹⁹log₁₀ γ — conditioning of the first m columnscondition number, and relative errorκ(ZᵀHZ), naiveκ(ZᵀHZ), pivotederror, naiveerror, orthonormalone problem, two rulesκ(ZᵀHZ) naive, γ = 10¹⁰2.7·10¹⁹κ(ZᵀHZ) pivoted, worst75error, naive0.9error, pivoted2·10⁻¹⁵κ(A) does not move across this axisand the answer moves by fourteen digits
Fig. 27 The bases the pencil is invariant to, at two constraints.
The share of shifts inside a pair of eigenvalues that count it wrongly, against the pair's separationTwo eigenvalues at 1 and 1 + gap, with six others spread around them, and 100 shifts placed strictly between the pair — where the count must read 3. Down to a separation of 10⁻¹² every shift reads it correctly. At 10⁻¹³ one of 100 does not, at 10⁻¹⁴ 0 do not, and at 10⁻¹⁵ none of them reads it correctly. The boundary sits where n‖A‖u puts it — 4807 for this matrix — because the floating-point count is the exact count of a matrix within that distance of A. This is the only place in the method where the answer can be wrong**, and it is wrong by a whole eigenvalue when it is: the failure is a miscount, not a small error.-15-13-11-9-7-5-300.250.50.751log₁₀ separation of the pairshare of shifts counted wronglyn‖A‖uthe only place it failswrong at 10⁻¹²0wrong at 10⁻¹⁴0wrong at 10⁻¹⁵100n‖A‖u4807wrong by a whole eigenvalueor not wrong at all
Fig. 28 And the band in which a computed multiplicity stops being a multiplicity.
The pivots Cholesky computes on a saddle-point matrix, and the one it stops at — row 11 of 14Every pivot Cholesky forms on K = [[H, Aᵀ], [A, 0]], in order, with the 10 taken inside H drawn as positive and the 11-th, which is not, drawn as the failure. It is not a near miss and it is not a rounding accident: the leading n × n block is positive definite, so the first 10 pivots are positive at any conditioning, and the next one is 0 − a₁ᵀH⁻¹a₁ where a₁ᵀ is the first row of A. That equals −(AH⁻¹Aᵀ)₁₁ = -3.7312, and the pivot the routine actually reports is -3.7312 — the same number to 15 digits. The failure has an address and a value, and both were available before the factorisation started.02468101210⁻³10⁻²10⁻¹110¹10²pivot, in order|pivot|negative, at row 11|S₁₁| from AH⁻¹Aᵀa failure with an addressrows before it10the pivot-3.7−(AH⁻¹Aᵀ)₁₁-3.7disagreement3.6·10⁻¹⁶it does not fail somewhereit fails at the constraint
Fig. 29 Why neither preconditioner can be a Cholesky.
Minimise ‖Ax − b‖ subject to Bx = d, solved as a weighted least-squares problem, three waysStack the constraint on top of the objective with a weight τ and solve the ordinary least-squares problem that results. In exact arithmetic the answer approaches the constrained one like 1/τ² — measured here as exactly four orders of error per two decades of τ, against a solution computed in BigInt rationals from the problem's own optimality conditions. What stops the limit is the solver and not the problem. The normal equations on the weighted problem form entries of size τ², so the constraint block is lost once τ² passes 1/u: the last weight at which they are within an order of the right answer is 10⁴ here, against 1/√u = 6.71·10⁷, and by τ = 10¹⁴ they are wrong by 14.3. Householder QR has no such ceiling and is at 4.8·10⁻¹⁵ at the same weight. Classical Gram–Schmidt is worse than either, at 4.9·10¹⁰.0246810121410⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹10⁴10⁷10¹⁰log₁₀ τ — the weight on the constraintrelative error against the exact answerτ = 1/√uGram–Schmidtnormal equationsHouseholder QRthe ceiling is the method'sHouseholder, τ = 10¹⁴4.8·10⁻¹⁵normal equations14Gram–Schmidt4.9·10¹⁰1/√u6.7·10⁷a constraint is a weight at infinityand the solver decides how far infinity is
Fig. 30 And a constrained problem imposed by a weight instead, from the least-squares field.

Where this goes next

00.3670080.7340171.101031.468031.835040eigenvalue of P⁻¹Kwritten down, then computeddistinct3at 16φ computed1.6off the closed form2.9·10⁻¹⁴1 − φ1φthe preconditioner's effect is a theoremand the golden ratio is in it Three eigenvalues, and two are the golden ratio Precondition a saddle-point system by the block diagonal of its own two definite pieces and the preconditioned matrix has exactly three distinct eigenvalues — 1, and the two roots of λ² − λ − 1. A minimal polynomial of degree three means three steps, at every conditioning, and the preconditioner nobody can afford turns out to be the statement the affordable ones are measured against. 01234510⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹log₁₀ κ(A)condition number, and relative errorκ(S)κ(ZᵀHZ)range-space errornull-space erroragainst a BigInt answerκ(S) at κ(A) = 10⁵4·10¹⁰κ(ZᵀHZ), all stops21range-space forward error5.3·10⁻⁶null-space forward error5.8·10⁻¹²both are the same algebraand only one squares Two ways to remove a constraint 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. κ(A) = 10⁶ throughout · κ(H) = 100 · the answer is the same answer for every basisorthonormal — κ(Z)1κ(ZᵀHZ)25.6relative error1.07·10⁻¹⁵first m basic — κ(Z)1.99·10⁸κ(ZᵀHZ)3.8·10¹⁶relative error0.0518pivoted basic — κ(Z)2.06κ(ZᵀHZ)31.9relative error6.71·10⁻¹⁶what the choice costsdensity, orthonormal1density, fundamental0.5κ(ZᵀHZ) ÷ κ(Z)², naive0.96error, pivoted choice6.7·10⁻¹⁶every one of them is a basisand one of them loses fourteen digits The basis nobody chose on purpose A method that eliminates a constraint has to pick a basis for its null space, and every basis is correct. Their condition numbers are eight orders apart, the reduced problem inherits the square, and the choice is usually made by a one-line rule nobody thought of as a numerical decision.

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.

Named objects

A flat tag is an object no other essay names yet.

Block preconditionerCondition numberConstraint preconditionerInertiaMatrix pencilMINRESNull spacePreconditioningReduced hessianSaddle-point systems