Where the augmentation puts the cost

Two earlier essays have measured what an approximate Schur complement costs a block preconditioner. Three eigenvalues, and two are the golden ratio found that the exact complement S = AH⁻¹Aᵀ gives three eigenvalues and three steps, and that every affordable replacement moves the spectrum away from them. One eigenvalue and two steps found that the way the blocks are assembled decides whether each approximate value is paid for once or twice. Both treated S as given and asked how well something cheap could imitate it.

There is a third move, and it is the one with a literature of its own. Do not approximate the Schur complement: change the problem so that its Schur complement becomes something cheap. The saddle-point system

[ H   Aᵀ ] [ x ]   [ f ]
[ A   0  ] [ y ] = [ g ]

has the same solution as

[ H + γAᵀA   Aᵀ ] [ x ]   [ f + γAᵀg ]
[ A          0  ] [ y ] = [ g        ]

for any γ, because the added term γAᵀ(Ax − g) is zero at the solution. The augmented system’s Schur complement is A(H + γAᵀA)⁻¹Aᵀ, which by a Woodbury identity equals (S⁻¹ + γI)⁻¹. As γ grows that tends to I/γ — a multiple of the identity, which costs nothing to apply. The augmented Lagrangian preconditioner blkdiag(H + γAᵀA, I/γ) uses exactly that.

The appeal is that a knob replaces an approximation. The cost is that the knob acts on the other block, and this essay measures where that cost goes.

A spectrum with a closed form at every γ #

The block-diagonal preconditioner with the exact Schur complement gives the eigenvalues 1 and (1 ± √5)/2. With an approximation Ŝ it gives 1 and (1 ± √(1 + 4ν))/2 for every generalised eigenvalue ν of the pencil (S, Ŝ) — the formula the triangular essay used. Here the approximation is I/γ for the augmented complement (S⁻¹ + γI)⁻¹, and its generalised eigenvalues can be written down: for each eigenvalue s of the original S,

ν = γs / (1 + γs)

which is between zero and one, tends to one as γs grows, and tends to zero as it shrinks. So the whole preconditioned spectrum is a function of γ and the eigenvalues of S, and the figure draws it both ways: the dashed curves from this formula, the marks from a Francis decomposition of the formed matrix. They agree to seven digits at every γ drawn, which is the check that the formula has been applied to the right S.

The formula says what the figure shows. The eigenvalues approach the golden-ratio pair from the inside — ν below one gives (1 + √(1 + 4ν))/2 below φ and (1 − √(1 + 4ν))/2 above 1 − φ — and how close they get is set by γ times the smallest eigenvalue of S. The largest distance from {1 − φ, φ}:

At γ = 0.01 the furthest eigenvalue is 0.617 from the pair, which is φ − 1 itself: a pair of eigenvalues sitting near one and near zero, as if the preconditioner had done nothing on that direction. At γ = 10⁶ it is 7.6·10⁻⁶ away, and MINRES takes five steps. The exact complement’s three steps are approached and not reached, because ν never equals one at any finite γ.

The γ that is large enough is set by the constraint’s conditioning #

“Large enough γ” has a scale, and the formula names it: γ has to be large against 1/s for the smallest eigenvalue s of S. S = AH⁻¹Aᵀ has eigenvalues that shrink with the square of the smallest singular value of A, so a worse-conditioned constraint pushes the required γ up by the square of the change.

Measured, on the same H with κ(A) raised from 10 to 1,000: at γ = 1 the smallest ν falls from 0.0558 to 6.2·10⁻⁶; at γ = 100, from 0.855 to 6.2·10⁻⁴; and at γ = 10⁴ it is 0.058, against 0.998 on the gentler constraint. A hundredfold worse κ(A) needs a ten-thousandfold larger γ for the same ν — γ = 10⁴ on the hard constraint gives 0.058, which is what γ = 1 gave on the easy one.

That is the first half of where the cost goes. The augmentation does not make a badly conditioned constraint easy; it converts the constraint’s conditioning into a demand on γ, and the demand grows like κ(A)².

The block the knob acts on #

The preconditioner applies (H + γAᵀA)⁻¹ at every step, and nothing about that block got easier. Its condition number is the one number that decides what the second half of the cost is, and the added term has rank m and norm γ‖A‖², so for large γ it grows like γ:

The sweep that follows runs the full solve at each γ, with MINRES taken to a relative residual of 10⁻¹⁴ and each application of (H + γAᵀA)⁻¹ counted as the conjugate gradient steps it would need, and it compares the answer against a solution computed exactly in rationals.

Read down the columns. The outer count falls, 21 to 6, which is the spectrum approaching the golden-ratio pair. The inner count rises, 14 to 43, which is the augmented block’s conditioning growing like γ; the count rises far more slowly than the square root of κ that the rate the condition number predicts would allow — κ rises by nearly eight orders and the count triples — which is consistent with the augmentation adding only m = 5 large eigenvalues, since conjugate gradients pays a few steps for each isolated outlier rather than the square root of the whole ratio. That reading is not measured separately. The product does not fall. Its minimum, 143, is at γ = 1, and from there it rises, to 258 at 10⁸. Taken at face value — outer steps times inner steps — the work of the whole solve is flat to within a factor of two across ten decades of the parameter that was supposed to make it easy.

On the hard constraint the same sweep shows the trade even more plainly: outer steps fall from 25 at γ = 1 to 8 at 10⁸, inner steps rise from 14 to 40, and the product goes from 350 to 320. The knob moves two thirds of the outer steps into inner steps and leaves the total where it was.

The digits go with it #

The last column is the part a step count cannot see. With the inner solve done exactly — no conjugate gradient tolerance to blame — the answer’s forward error rises from 6.6·10⁻¹⁶ at γ = 0.01 to 2.4·10⁻⁸ at γ = 10⁸. Seven and a half digits, from a change that leaves the solution of the system mathematically untouched.

The mechanism the numbers are consistent with is the condition number as an amplifier, applied to a block the augmentation made ill-conditioned on purpose. Every application of the preconditioner is a Cholesky solve with H + γAᵀA, and the augmented right-hand side f + γAᵀg is γ times larger than the original in the directions of the constraint; the solution x is recovered by cancelling γAᵀAx against γAᵀg to leave something of the original size, and a cancellation of that kind loses about u times the ratio of what is cancelled to what is left. The measurement does not separate the error of the Cholesky factor from the error of the cancellation. What it establishes is that the loss tracks the block’s condition number. Across the table the forward error tracks that condition number: forward error divided by κ is 7.6·10⁻¹⁸ at γ = 0.01, 7.6·10⁻¹⁸ at 10², 7.4·10⁻¹⁸ at 10⁴, 2.6·10⁻¹⁸ at 10⁶ and 6.1·10⁻¹⁸ at 10⁸ — a constant to within a factor of three over eight orders of κ, and below κu by a factor of fifteen to forty. The exception is γ = 1, where the error is 5.0·10⁻¹⁵ against a κ of 84 — a ratio eight times the others’, which is recorded and not explained.

On the hard constraint the answer starts worse and ends in the same place: 8.9·10⁻¹³ at γ = 1, where the constraint’s own conditioning is already visible, and 2.2·10⁻⁸ at 10⁸. At large γ the augmentation caps the achievable accuracy at a few times 10⁻¹⁸ times κ(H + γAᵀA), whatever the constraint started at.

This is the pattern a condition number sent to infinity found in interior-point methods and the regularisation that legalises every order found in quasi-definite factorisation: a parameter introduced to make one part of a computation easy, whose price is paid in the accuracy of another part, and which therefore has an optimum rather than a direction to push in. A constraint is a weight at infinity is the nearest relative. The method of weighting approximates a constraint by adding τ²AᵀA to a least-squares problem and takes τ to infinity; the augmented Lagrangian adds γAᵀA and keeps the constraint exactly, so it needs no limit — but its block H + γAᵀA is the same kind of ill-conditioned matrix the weighting method factorises, and it costs digits for the same reason.

Two tolerances, and which count is a count #

The spectrum figures report MINRES steps to a relative residual of 10⁻¹⁰ and the sweep reports them to 10⁻¹⁴, and the difference between the two columns is informative rather than a nuisance. At γ = 0.01 the count goes from 13 to 21 when four more digits are asked for. At γ = 1 and γ = 100 it goes from 11 to 11. At γ = 10⁴ it goes from 7 to 10, and at 10⁶ from 5 to 8.

Eleven is not an accident. With five constraints and seven unconstrained directions the preconditioned matrix has at most 1 + 2·5 = 11 distinct eigenvalues — one, and a pair for each generalised eigenvalue — so MINRES terminates in eleven steps whatever the tolerance. At γ = 1 and 100 the eleven values are spread widely enough that nothing short of all of them will do, and the count is a count of values: tightening the tolerance costs nothing, because the twelfth step does not exist.

At γ = 10⁴ and 10⁶ the values have gathered near the golden-ratio pair and a polynomial of low degree is small on all of them long before it has to vanish on each one. The count is then a rate, and it pays for digits: about three steps for four more of them, on both settings. The augmentation has changed the kind of convergence as well as its speed, from finite termination on a scattered spectrum to geometric convergence on a clustered one — which is why the outer count keeps falling slowly past the point where the spectrum looks finished.

And at γ = 0.01 the eight extra steps are consistent with one eigenvalue sitting at about −5.9·10⁻⁴, very near zero: a residual polynomial must equal one at the origin and be small at that eigenvalue, and the closer the eigenvalue is to the origin the more degree each digit costs. That reading is not separately measured; the eight steps are.

Two harder systems, and where the least work falls #

The sweep above is one system, and on it the least work is at γ = 1, near the small end of the range. Two harder variants move that minimum, and they move the accuracy optimum with it.

On the constraint conditioned at κ(A) = 1,000 the outer count starts much higher — 27 steps at γ = 0.01, 25 at 1, 21 at 100 — because the smallest generalised eigenvalue is 6·10⁻⁸, 6·10⁻⁶ and 6·10⁻⁴ there and the preconditioner has done almost nothing on those directions. It falls to 12 at γ = 10⁴, where the smallest ν reaches 0.058, and to 8 at 10⁸. The inner count rises from 14 to 40 as before. So the product falls, from 378 to 228 at γ = 10⁴, and then rises again, to 320 at 10⁸. On this system the augmentation genuinely saves work — forty per cent at its best — and the saving sits at a γ ten thousand times larger than the easy system’s.

The answer’s accuracy has an interior optimum here too, which the easy system did not show. The forward error is 6.6·10⁻¹³ at γ = 0.01 and 8.9·10⁻¹³ at γ = 1 — the constraint’s own conditioning already costing digits before any augmentation, the effect two ways to remove a constraint measured in the eliminations that inherit it — and then 2.7·10⁻¹⁴ at γ = 100, thirty times better, before rising to 1.2·10⁻¹², 3.3·10⁻¹¹ and 2.2·10⁻⁸. A moderate augmentation improves the answer on a badly conditioned constraint; a large one takes the improvement back and more.

On a Hessian conditioned at κ(H) = 10⁴ with the gentle constraint the same shape appears for a different reason. The augmented block’s condition number first falls, from 4,490 at γ = 0.01 to 1,420 at γ = 1, because the added term lifts the directions in which H alone is smallest; only after that does it grow like γ, to 6.2·10¹⁰ at γ = 10⁸. The least work, 176, is at γ = 1, and so is the best forward error, 5.2·10⁻¹⁵ against 4.6·10⁻¹⁴ at γ = 0.01 and 2.3·10⁻⁷ at 10⁸. The inner count at 10⁸ is 56, the largest in any of the three sweeps.

Laid side by side, the three sweeps say that the knob has an optimum in both currencies and that the optimum moves with the problem: γ = 1, 10⁴ and 1 for the work, and γ = 0.01, 100 and 1 for the digits. Two of the three least-work settings sit where the smallest generalised eigenvalue is about 0.06, and the third where it is 0.40 — suggestive of a rule written in the smallest ν rather than in γ, and three systems are not enough to call it one. At the top of the range all three agree on the price: forward error divided by κ(H + γAᵀA) is 6.1·10⁻¹⁸, 5.3·10⁻¹⁸ and 3.7·10⁻¹⁸ at γ = 10⁸.

What the knob is worth, stated fairly #

Three things the measurement supports, and one it does not.

The spectrum is exactly what the theory says, at every γ. That is not a small thing: it means the step count can be predicted from γ and the smallest eigenvalue of S without running anything, and the prediction has seven-digit agreement behind it.

The outer iteration does get easier, and the gain saturates. From γ = 10⁴ upward the spectrum is within 10⁻³ of the golden-ratio pair and MINRES takes between five and ten steps, against the exact complement’s three. The last factor of two is the most expensive to buy.

The inner solve pays most of it back. Counted in conjugate gradient steps with an unpreconditioned inner solve, the total is flat on the easy system and falls by at most forty per cent on the hard constraint. That is the part of the measurement most specific to its setting. A production augmented Lagrangian code does not solve with H + γAᵀA by unpreconditioned conjugate gradients; it uses a preconditioner or a multigrid method designed for the augmented block, whose whole difficulty is that the added term has rank m and norm γ‖A‖², so the block has m eigenvalues that grow with γ and n − m that do not, and the inner method must handle both. The flat product above is the cost when that design problem has not been solved, and it is what the parameter buys on its own.

What it does not support is a preference between this preconditioner and the diagonal Schur approximations of the golden-ratio essay. On this problem at γ = 1 the augmented preconditioner takes 11 outer steps, the same as a diagonal approximation A diag(H)⁻¹Aᵀ; which is cheaper in total depends on the cost of the inner solves in each, and the comparison would need the same inner preconditioner for both. Nor does it rank the augmentation against a preconditioner that need not know the constraint, which keeps A exactly and whose spectrum does not move with κ(A) at all — the one property the hard-constraint sweep above shows the augmentation lacks.

The answer’s digits are the part that transfers. Whatever inner solver is used, a forward error that tracks κ(H + γAᵀA) is a property of the augmented system, and the only way to keep the digits is to keep γ below the point where γ‖A‖² overwhelms H — which is, roughly, the point past which the outer iteration stopped improving much anyway.

Where this goes from here #

An inner solve designed for the augmented block. The flat product above assumes the inner solve pays the full square root of κ(H + γAᵀA). A preconditioner for the augmented block that handles the rank-m, norm-γ term — a block solver that treats the constraint directions separately — would change the inner column, and whether the product then falls with γ, and what that does to the digits, is the measurement that decides whether the augmentation is worth doing at all.

An inexact inner solve. Every count above uses an exact application of (H + γAᵀA)⁻¹ inside the preconditioner. Solving it to a tolerance changes the preconditioner at every outer step, and the triangular essay’s nilpotent structure and this essay’s closed-form spectrum are both exact only with exact solves. The relation between the inner tolerance, the outer count and the forward error is unmeasured, and the collection’s measurement of how much accuracy an inner solve is worth, the accuracy that is thrown away, is about Newton steps rather than preconditioners and would have to be made again here.

And choosing γ. The three sweeps have optima for work at γ = 1, 10⁴ and 1 and for digits at γ = 0.01, 100 and 1, so the choice is a balance between two measurable quantities rather than a push in one direction. A rule that picks γ from the smallest eigenvalue of an approximation of S — available at the cost of a few Lanczos steps — and a statement of how close it lands to that balance would turn this essay’s table into a method.