A tolerance is not an accuracy

A square root the residual does not pay took the block upper-triangular preconditioner for a saddle-point system, which with exact solves makes every eigenvalue one in Jordan blocks of size two and lets GMRES finish in two steps, and asked what an inexact solve with the Hessian block does to it. It modelled the inexact solve as an exact solve with a perturbed block, $\tilde H = H + \tau H^{1/2}EH^{1/2}$ for a fixed $E$ of unit norm, so that the outer iteration was still GMRES on one operator. The Jordan pairs split by $\sqrt\tau$; the residual after two steps grew like $\tau$, because GMRES’s polynomial has a double root at one and cancels a split pair to second order; and a rule written in $\tau$ alone, $\lceil 2\log(\mathrm{tol})/\log(c\tau)\rceil$ with $c$ the two-step residual over $\tau$ at one tolerance, predicted the step count within one on 228 of 234 runs.

It said plainly where the model stops. A real inner iteration returns a different error for every vector it is given, so there is no single perturbed operator, and GMRES has to become flexible GMRES, which keeps every preconditioned vector it used rather than one preconditioner. “The prediction with a sign is that with each inner solve stopped at a relative residual of $\tau$, flexible GMRES’s two-step residual still follows $\tau$ rather than $\sqrt\tau$, with a constant at most three times the fixed perturbation’s, and that the step rule above predicts its count within one step on at least nine runs in ten.”

Measured on the same eighteen systems, the first clause holds once “follows $\tau$” is read correctly, and the other two fail. All three turn on the same fact: the tolerance an inner solve is stopped at is not the accuracy it delivers.

The same systems, with a real inner solve #

The eighteen systems are the earlier essay’s: three shapes — twelve unknowns and five constraints, nine and six, fourteen and two — with Hessian blocks of condition number 1, $10^2$ and $10^4$ and constraint blocks of condition number 10 and $10^4$. The Schur complement is still applied exactly. The Hessian block is now solved by conjugate gradients, started from zero for every vector and stopped at the first step whose relative residual is at most $\tau$, for eleven tolerances from $10^{-12}$ to $10^{-2}$. The outer iteration is flexible GMRES to a relative residual of $10^{-10}$.

The fixed model is run through the same flexible GMRES with its one fixed perturbed block, so that both report the same quantity — the unpreconditioned residual $\|b - Kx\|/\|b\|$ — and the only difference between them is the inner solve. Run that way the model behaves as it did under ordinary GMRES: its two-step residual grows with a slope between 0.8 and 1.5 in $\tau$ on every system.

What the inner solve actually delivers #

Before asking what the outer iteration does with an inner solve, it is worth asking what the inner solve does. Every application can be checked: the block is small enough to solve exactly as well, and the error of conjugate gradients’ answer can be measured in $H$’s own norm, the energy norm in which the fixed model’s $\tau$ is stated.

On the six systems with $\kappa(H) = 1$ the Hessian is a multiple of the identity, conjugate gradients’ first step is exact, and every inner solve returns an answer correct to rounding — $6 \cdot 10^{-16}$ — whether it was asked for $10^{-12}$ or $10^{-2}$. Those six systems are solved by flexible GMRES in two steps at every tolerance, which is the exact preconditioner’s count. The fixed model, asked for the same $\tau$, needs seven to nine steps at $10^{-2}$.

On the other twelve the delivered error moves in steps. The blocks have nine to fourteen unknowns, and conjugate gradients on a matrix of $n$ distinct eigenvalues terminates in at most $n$ steps in exact arithmetic; long before that a single step can cut the residual by several orders of magnitude, from just above one tolerance to far below the next. On the nine-by-six system with $\kappa(H) = 10^2$ the delivered residual is $2 \cdot 10^{-9}$ at every tolerance from $10^{-8}$ to $10^{-4}$: nine inner steps reach it, eight do not reach $10^{-4}$. Across all 198 runs, the median residual delivered is ten times smaller than the tolerance or more on 140, and the energy error on 163.

This is not a property of conjugate gradients alone. A stopping test checks a quantity at discrete steps and stops at the first one that passes; the quantity at that step can be anywhere below the threshold. A tolerance that reads its own residual found the same gap in inexact Newton’s inner solves, where a forcing term set the tolerance and the work done depended on how far each inner solve overshot it. Here the overshoot is the main effect rather than a correction to it.

The case a good inner preconditioner makes #

The six systems with $\kappa(H) = 1$ look like a degenerate corner of the sweep, a block so easy that the question of inner accuracy does not arise. In a real code they are the case the inner solve is built to approach. The Hessian block of a discretised problem is not solved by unpreconditioned conjugate gradients; it is solved with a multigrid cycle or an incomplete factorisation as the inner preconditioner, and the whole point of that preconditioner is to make the block’s preconditioned spectrum cluster near one. Changing the condition number on purpose is the general statement, and a rate that does not notice the size is what a V-cycle does to it: a tenfold reduction per cycle at every grid size.

A clustered spectrum is exactly where conjugate gradients takes the largest steps. The rate the condition number predicts found the textbook bound in $\sqrt{\kappa}$ to be a bound and not an estimate, and the gap between them is largest when the eigenvalues bunch: the iteration removes a cluster’s error almost at once and the next residual is far below the last. So the better the inner preconditioner, the larger the jump past any given tolerance on the step that crosses it, and the further the delivered accuracy sits from the nominal one. The extreme of that is the identity block here, where one step delivers rounding error at every tolerance.

That turns the overshoot from a curiosity of small blocks into something a solver designer has to plan for. The tolerance passed to an inner solve is a ceiling on its error, not a forecast of it, and how far below the ceiling the answer lands depends on the inner preconditioner’s quality — which is the one thing that is usually tuned independently of the outer iteration.

The residual follows the error, not its square root #

The first clause of the prediction was that the outer residual would follow $\tau$ and not $\sqrt\tau$ — that the defective eigenvalue’s square-root sensitivity, visible in the split of the Jordan pairs, would again not reach the residual. With the tolerance replaced by the error delivered, that holds.

On eleven of the twelve systems the two-step residual climbs with a slope between 0.69 and 1.06 in the delivered error; the twelfth, where conjugate gradients delivers $2 \cdot 10^{-9}$ at most tolerances and flexible GMRES finishes in four steps regardless, has too little range to fit a slope to and gives 0.34. Nothing sits on the line of slope one half. The mechanism the earlier essay gave survives the change of iteration: a residual polynomial with a double root at one cancels each pair split by $\sqrt{\varepsilon}$ to order $\varepsilon$, and flexible GMRES, which minimises the same residual over a space spanned by the preconditioned vectors it actually used, has that polynomial available whenever the first two inner solves are close enough to each other.

So the square root still stays in the spectrum. That part of the argument is about GMRES’s polynomial, and it does not care whether the perturbation came from a model or from an iteration.

A constant the model does not supply #

The second clause fails, and by more than the prediction’s factor of three. The fixed model’s two-step residual is a constant times its $\tau$ — between 0.5 and 360 depending on the system, 44 at the median. Divide each flexible run’s two-step residual by its own system’s model constant times the energy error the inner solves delivered, and the ratio should be about one if a real inner solve of a given accuracy costs what a fixed perturbation of that accuracy costs.

It has a median of 12 on the systems whose constraint block has $\kappa(A) = 10$, and 99 on those with $\kappa(A) = 10^4$, with a tenth of the runs above 160 and 2,000 respectively. Twenty-four runs of 66 sit within three of the model on the well-conditioned constraints, thirteen of 66 on the ill-conditioned ones. The figure above shows it as a band of points lying one to three decades over the model’s line, roughly parallel to it.

A real inner solve’s error is not a random perturbation. Conjugate gradients removes error fastest in the directions of $H$’s extreme eigenvalues and leaves it last in the interior and lower part of the spectrum, so the error it hands back is concentrated in a few eigenvectors of $H$, and it is a different few each time it is applied to a different vector. The fixed $E$ spreads its perturbation evenly. Which of the two the constraint block amplifies more depends on how $A$ couples to those eigenvectors, and the ratio’s dependence on $\kappa(A)$ says that the coupling matters; why it comes out at ten to a hundred, rather than at one or at a thousand, is not established here.

The rule was a fixed operator’s rule #

The third clause fails too. Read in the delivered error, the rule predicts the flexible step count within one on 158 runs of 198; read in the tolerance, as the prediction stated it, on 156. Both are about four in five, where nine in ten was predicted and the fixed model managed 228 of 234.

The misses gather where they would be expected. Thirty-one of the forty are at tolerances of $10^{-5}$ and looser, and twenty-six are on the systems with $\kappa(H) = 10^4$, where an inner solve’s error varies most from one application to the next. Twenty-seven of them predict too many steps. The rule assumes each pair of outer steps divides the residual by the same factor $c\tau$, which is true of GMRES on a fixed operator whose spectrum sits in clusters of radius $\sqrt\tau$. Flexible GMRES on a sequence of different inner solves has no such geometric rate. A late application that happens to be accurate can finish the solve early; a loose one can stall it. The constant read off one tolerance does not carry to another, because the errors at different tolerances are not multiples of one another — they are different vectors.

The real solve costs no more #

All three clauses of the prediction were about how well a model describes the real thing. The practical question is what the real thing costs, and the answer is better than the model.

Over all 198 runs, flexible GMRES with a real inner solve takes no more outer steps than the fixed model on 181. At the median it takes the same number up to $10^{-6}$ — three to four — and fewer beyond it: four against five at $10^{-4}$, five against seven at $10^{-3}$, seven or eight against nine at $10^{-2}$. The seventeen runs where it takes more are scattered through the tolerances with a cluster at the loosest. At $10^{-2}$ it takes more on seven of the eighteen systems, up to fourteen steps against nine on the nine-by-six system with both blocks ill conditioned. That is the one setting where a single bad inner solve can cost several outer steps, and it is the setting where the fixed model looked worst anyway.

The reason is the overshoot. A solve asked for $10^{-4}$ that delivers $10^{-9}$ is, for the outer iteration, a solve at $10^{-9}$, and the outer iteration’s cost is set by what it is given.

Loosening the inner solve bought nothing #

That leaves the question an inner tolerance is chosen to answer: how loose can it be before the extra outer steps cost more than the inner work saved?

On these systems, never. With $\kappa(H) = 10^2$ the median total is 34 inner steps at $10^{-12}$ and climbs to 85 at $10^{-2}$; with $\kappa(H) = 10^4$ it is 70 at $10^{-12}$, between 65 and 74 up to $10^{-6}$, and 118 at $10^{-2}$. The tightest tolerance is within a fifth of the cheapest on both, and the loosest is the dearest on both.

The arithmetic is visible in the two counts. One inner solve on a nine-to-fourteen-unknown block takes 11 to 26 steps at $10^{-12}$ and 6 to 18 at $10^{-2}$ — a saving of about a third — while the outer iteration goes from two to four steps to eight or nine, and every outer step needs one more inner solve. A loose inner tolerance trades a third of each solve for twice as many solves.

That arithmetic belongs to small blocks. On a block of a million unknowns conjugate gradients’ steps grow like the logarithm of $1/\tau$ times the square root of $\kappa(H)$, there is no finite termination to fall into, and the trade looks different; it is the question a tolerance that reads its own residual and two forcing terms, one floor answer for Newton’s method, where the inner tolerance is chosen step by step. What the small blocks do show is that the saving a loose tolerance offers has to be measured in delivered accuracy and in total work, and that the nominal tolerance predicts neither.

What eighteen small blocks do not show #

Every inner solve here starts from zero, and the Schur complement is exact. A real code would warm-start each inner solve from the previous one and approximate the Schur complement too, as a preconditioner that need not know the constraint does with the Hessian block; both would change what each application delivers. The blocks are small enough that conjugate gradients terminates, which is what makes the overshoot so large; on larger blocks the delivered error would sit closer to the tolerance, and the gap between the real solve and the fixed model would be expected to narrow. The explanation offered for the constant — that conjugate gradients’ error is concentrated in a few eigenvectors and the constraints amplify those — is consistent with the dependence on $\kappa(A)$ and is not tested by any measurement here. One eigenvalue and two steps and three eigenvalues, and two are the golden ratio set the exact counts all of this is measured against.

Still open: an inner error in chosen directions, and a test on delivered accuracy #

Where the error lies. If the constant comes from conjugate gradients concentrating its error in a few eigenvectors of $H$, a fixed model built with its perturbation in those directions should reproduce it. The prediction with a sign is that a fixed perturbation $\tau H^{1/2}EH^{1/2}$ with $E$ supported on the eigenvectors of $H$’s three smallest eigenvalues gives a two-step residual within a factor of three of flexible GMRES’s at matched energy error on at least nine of the twelve systems, where the evenly spread $E$ manages it on fewer than half.

Stopping on the error delivered. An inner solve can estimate its own energy error cheaply from the quantities conjugate gradients already carries. The prediction is that stopping each inner solve when that estimate falls below a target $\eta$, rather than when the residual falls below $\tau$, makes the step rule — written in $\eta$ — count the outer steps within one on at least 180 of 198 runs, because the error each application delivers then sits within a factor of a few of the target instead of anywhere below it.