A square root the residual does not pay

One eigenvalue and two steps put the off-diagonal block back into a saddle-point preconditioner. With $K$ the saddle-point matrix, $H$ its Hessian block, $A$ its constraints and $S = AH^{-1}A^{\mathsf T}$ the Schur complement, the block upper-triangular preconditioner

$$P = \begin{bmatrix} H & A^{\mathsf T} \\ 0 & -S \end{bmatrix}$$

makes $P^{-1}K = I + N$ with $N^2 = 0$ and $N \neq 0$. Every eigenvalue is one, in $m$ Jordan blocks of size two, and the minimal polynomial is $(\lambda - 1)^2$: GMRES takes exactly two steps. The same essay found that the computed spectrum of $I + N$ is not a point but a ring of radius about $\sqrt{u\|N\|}$, the square-root sensitivity of a defective eigenvalue, and measured it against the prediction to a few per cent.

That essay’s exactness had a condition attached. Applying $P^{-1}$ means solving with $H$ and with $S$, and in a code of any size the solve with $H$ is itself iterative and stops at a tolerance. It closed on that: “The nilpotent part is exactly nilpotent only with exact solves. The next measurement on the triangular form is the relation between the inner tolerance and the extra steps it costs, and whether a defective eigenvalue’s square-root sensitivity reappears as a square root in that relation.”

It reappears in the eigenvalues, exactly as predicted, and not in the steps.

An inner solve accurate to τ #

The inexact solve is modelled as an exact solve with a perturbed block. $H$ is replaced by

$$\tilde H = H^{1/2}\,(I + \tau E)\,H^{1/2},$$

with $E$ a fixed symmetric matrix of unit norm. Then $\tilde H^{-1}$ differs from $H^{-1}$ by a relative $\tau$ in $H$’s own norm, which is what an inner conjugate-gradient solve stopped at an energy-norm tolerance $\tau$ leaves behind. Two simplifications are made on purpose. The perturbation is fixed, so the outer iteration is GMRES on one operator rather than a flexible method on a changing one. And it is scaled to $H$, so it stays positive definite at every conditioning — a perturbation of size $\tau\|H\|$ would make $\tilde H$ indefinite once $\tau\kappa(H)$ reached one, which is a property of the model rather than of any solve. The Schur complement stays exact; the earlier essay already swept an approximate one, and where the augmentation puts the cost found a way to make the exact one nearly free.

The systems are that essay’s eighteen: three shapes (twelve unknowns with five constraints, nine with six, fourteen with two), Hessians with $\kappa(H) = 1$, $10^2$ and $10^4$, constraints with $\kappa(A) = 10$ and $10^4$. The inner accuracy runs over thirteen decades, from $10^{-14}$ to $10^{-2}$.

The eigenvalues split by a square root #

The figure at the top of the page is the drawn system: twelve unknowns, five constraints, $\kappa(H) = 100$. It has two populations. The ten eigenvalues that sat in the five Jordan blocks leave one along a line of slope one half: about $0.5\sqrt\tau$, so at $\tau = 10^{-8}$ they are $5\cdot 10^{-5}$ away. The other seven leave along a line of slope one. No eigenvalue is drawn below $\tau = 10^{-8}$, and the reason is a finding of its own, two sections down: below that accuracy the standard decomposition does not finish.

The square root is the arithmetic of a Jordan block. The matrix $\begin{bmatrix} 1 & 1 \\ \varepsilon & 1 \end{bmatrix}$ has eigenvalues $1 \pm \sqrt\varepsilon$: a perturbation in the corner a Jordan block leaves empty moves the eigenvalue by its square root, and a condition number for one eigenvalue measured the same exponent becoming an eighth on an eight-by-eight block. Here each of the $m$ blocks gets its own $\varepsilon$, of either sign, from how the perturbed solve couples the two halves of the block.

Drawn in units of $\sqrt\tau$, the split is a fixed picture. Each pair sits at $1 \pm \sqrt{\tau a_k}$ for a number $a_k$ of either sign, so it lies on the real axis or the imaginary one, symmetric about one. At $\tau = 10^{-6}$ there are two real pairs and three imaginary, the nearest at $0.14\sqrt\tau$ and the furthest at $0.70\sqrt\tau$. Across the dial, from $10^{-8}$ to $10^{-2}$, the positions hold; the seven that move like $\tau$ fall into the centre as $\tau$ shrinks, $0.006\sqrt\tau$ away at $10^{-6}$ and $0.0006\sqrt\tau$ at $10^{-8}$.

On all eighteen systems the same: from $10^{-8}$ to $10^{-2}$ the median distance of the paired eigenvalues from one is between 0.42 and 0.62 times $\sqrt\tau$, and its slope against $\tau$ is 0.50 on every one of them. The other eigenvalues’ slope is 1.00 on every one.

Which way each pair goes, predicted #

The positions are not just of the right size; they can be written down before anything is decomposed. Write the perturbed solve’s error as $D = \tilde H^{-1} - H^{-1}$, a matrix of size $\tau$. To first order each Jordan pair moves to

$$\lambda = 1 \pm \sqrt{\mu_k}, \qquad \mu_k \in \operatorname{eig}\bigl(S^{-1} A D A^{\mathsf T}\bigr),$$

an $m \times m$ problem built from the constraint block and the inner solve’s error alone. $S^{-1}ADA^{\mathsf T}$ is similar to a symmetric matrix, so every $\mu_k$ is real, and its sign decides the pair’s direction: a positive $\mu_k$ splits along the real axis, a negative one along the imaginary. On the drawn system the five $\mu_k/\tau$ are 0.374, 0.235, $-0.020$, $-0.269$ and $-0.486$ — two real pairs and three imaginary, at 0.61, 0.48, 0.14, 0.52 and 0.70 times $\sqrt\tau$, which are the rings on the figure above.

Measured on every system, the computed pairs sit within 2.0% of $1 \pm \sqrt{\mu_k}$ at $\tau = 10^{-7}$ and 2.2% at $10^{-6}$; the gap grows to 6.7% at $10^{-5}$ and 20% at $10^{-4}$, which is the next term of the expansion, a relative correction that grows like $\sqrt\tau$. The signs vary with the problem and not with $\tau$: two real and three imaginary pairs on four of the six twelve-by-five systems, three and two on the other two, three and three on every nine-by-six system, one and one on every fourteen-by-two. Because $\mu_k/\tau$ is a fixed matrix’s spectrum, the whole pattern in units of $\sqrt\tau$ is fixed, which is why dragging $\tau$ moves nothing.

The residual does not split #

So the eigenvalues have the square-root sensitivity the question asked about. The question was whether the steps do, and the green line in the first figure already answers it. GMRES’s residual after two steps follows slope one: $7.4\cdot 10^{-9}$ at $\tau = 10^{-8}$ — the drawn system’s first residual is 0.25 and its second about $0.74\tau$ — and $7.3\cdot10^{-3}$ at $10^{-2}$. It is linear in $\tau$ across eight decades.

The reason is the polynomial. After two steps GMRES has the freedom of every polynomial of degree two with $p(0) = 1$, and one of them is $(1 - z)^2$, the minimal polynomial of the exact preconditioned matrix. Applied to the perturbed one, $I + N + F$ with $F$ of size $\tau$:

$$\bigl(I - (I + N + F)\bigr)^2 = (N + F)^2 = NF + FN + F^2 ,$$

because $N^2 = 0$. Every surviving term carries at least one factor of $F$, so the residual after two steps is of order $\tau\|N\|$, not $\sqrt\tau$. The polynomial has a double root at one, and a double root cancels a pair of eigenvalues at $1 \pm \sqrt{\tau a}$ to second order: $(1 - z)^2$ at $z = 1 \pm \sqrt{\tau a}$ is $\tau a$. The square root is real; GMRES’s polynomial squares it before it reaches the residual.

Divide each quantity by the power it is claimed to follow. The split over $\sqrt\tau$ is a flat line on every system, bunched between 0.42 and 0.62. The two-step residual over $\tau$ is flat on the nine systems with $\kappa(A) = 10$, to within a factor of 1.4, at levels from 0.013 to 0.74. On the nine with $\kappa(A) = 10^4$ the linear term is far smaller — as little as $3\cdot 10^{-7}\tau$ — and as $\tau$ grows the $F^2$ term overtakes it, so the line rises, by up to a factor of a hundred. Fitted from $10^{-10}$ to $10^{-2}$, the two-step residual’s slope is between 0.86 and 1.27 on every system, and its ratio to $\sqrt\tau$ rises between two hundred and a hundred thousand times. The constant is how much of the starting residual $NF + FN$ can reach, and it depends on the blocks; the power is never one half.

Every pair of steps costs one factor of τ #

Two steps leave a residual of order $\tau$. What happens after that is the same argument again: the next two steps can apply $(1 - z)^2$ once more, and $(N+F)^4 = (NF + FN + F^2)^2$ is of order $\tau^2$.

The histories show the pairing. At $\tau = 10^{-4}$ GMRES’s residuals run 0.25, $7.4\cdot10^{-5}$, $4.1\cdot10^{-6}$, $7.7\cdot10^{-10}$, $8.9\cdot10^{-11}$: drops of 3,400 and 5,300 at the second and fourth steps, of 18 and 9 at the third and fifth. At $10^{-8}$ the second step lands at $7.4\cdot10^{-9}$, the third at $4.1\cdot10^{-10}$, and the fourth at $10^{-16}$. At $10^{-2}$, where $c\tau$ is no longer small, the pattern blurs: the residuals fall by factors between 8 and 56 at every step, and it takes eight. The odd steps are not wasted — they buy the part of the next factor that a degree-three polynomial can already reach — but the large drops come two at a time.

That makes the cost predictable from $\tau$. If each further two steps divides the residual by $c\tau$, with $c$ the two-step residual over $\tau$, then reaching a tolerance tol takes

$$\text{steps} \;\approx\; \left\lceil \frac{2\log \text{tol}}{\log c\tau} \right\rceil,$$

at least two. $c$ is read once per system, at $\tau = 10^{-6}$, from the second residual.

Over 234 runs the rule is exact on 179 and within one step on 228. The misses by one are at the boundaries, where the partial drop at an odd step decides whether a count lands one side or the other. The six misses by two or three are all at the loosest accuracies, $10^{-3}$ and $10^{-2}$, where $c\tau$ is no longer small and “a factor of $c\tau$ per two steps” stops being a clean statement — five of them over-predict and one under-predicts.

Read the other way, the rule is a price list for the inner tolerance. At the drawn system’s $c$ of about 0.7, the two-step answer holds to $\tau = 10^{-10}$; four steps hold to about $10^{-5}$; six to about $10^{-3}$. Loosening the inner solve by five decades costs two outer steps, each of which costs one more inexact solve with $H$. That is the trade the accuracy that is thrown away found for inexact Newton — an inner tolerance far below what the outer method can use buys nothing — with the difference that here the outer method’s appetite is a clean power of $\tau$ and can be stated in advance.

The triangular form keeps its lead #

The other reason to want this number is the comparison the earlier essays set up. The block-diagonal preconditioner keeps symmetry, which allows MINRES and its constant storage; the triangular one needs GMRES, which stores every vector, and repays it with two steps instead of three. Once the inner solve is inexact both lose their clean counts. The question is which loses more.

The triangular form keeps its lead at every $\tau$. At the median GMRES takes 2 steps up to $10^{-9}$, 3 or 4 from $10^{-8}$ to $10^{-4}$, and 5 and 6 at the last two; MINRES takes 3 up to $10^{-10}$, 5 or 6 from $10^{-9}$ to $10^{-5}$, then 8, 9 and 12. On every one of the eighteen systems at every one of the thirteen accuracies, GMRES needs at most as many steps as MINRES. The ratio runs from the exact case’s two to three toward one to two.

The block-diagonal spectrum explains its own pattern. Its three eigenvalues — $1 - \varphi$, $1$ and $\varphi$ from three eigenvalues, and two are the golden ratio — are not defective, so they move like $\tau$, not $\sqrt\tau$, and MINRES’s degree-three polynomial cancels three clusters of width $\tau$ to first order each. On the nine systems with $\kappa(A) = 10$ its residual after three steps is $0.25\tau$ to $0.50\tau$, flat across the eight decades, and on the drawn system at $10^{-4}$ the sixth step reaches $1.2\cdot10^{-9}$, a second factor of about $\tau$: the same law with three in place of two, and the same doubling of the count when $\tau$ crosses the tolerance. A defective preconditioned matrix is more sensitive than a diagonalisable one, and the measurement says GMRES does not pay for that sensitivity, because a polynomial that matches a Jordan block’s multiplicity sees its eigenvalues’ split only squared.

A spectrum the decomposition does not finish #

The eigenvalues above stop at $\tau = 10^{-8}$ because below it the Francis iteration — the decomposition every library runs, which the algorithm the libraries actually run walks through — does not converge on these matrices. With a budget of 400 double steps per row it converges on every system from $10^{-8}$ up. At $10^{-9}$ it fails on nine of the eighteen, at $10^{-10}$ on sixteen, at $10^{-11}$ and $10^{-12}$ on seventeen, and at $10^{-14}$ on eight. The exact case does not have the problem: with $\tau = 0$ the earlier essay’s matrices converge in between 0 and 50 double steps. What the stalled runs have in common is ten eigenvalues within a few times $10^{-5}$ of one another and of one, arranged as pairs that are almost Jordan blocks; the exact case, where the same ten are an exact Jordan structure, converges at once. Why the near case is so much harder than the exact one is not traced here.

The failure is quiet if the convergence flag is not read. An unconverged quasi-triangular matrix still has a diagonal, and read as a spectrum it is a plausible one: on the drawn system at $10^{-10}$ it reports five real pairs at 0.58 to 1.10 times $\sqrt\tau$ — the right size, the wrong directions, three of them real where the pairs are imaginary. A second trap sits beside it. The routine that reads the two-by-two blocks off the quasi-triangular form calls a subdiagonal zero below $10^{-12}$ of the matrix’s norm, and an imaginary pair split by $\sqrt\tau$ for $\tau$ below about $10^{-7}$ has a block whose subdiagonal sits under that line, so even a converged decomposition read at the default reports it as two real eigenvalues. The form a real matrix can reach is the essay on what those blocks are; here they are read at the iteration’s own deflation tolerance, $10^{-14}$, and only where it converged.

What this says about a ring #

The earlier essay drew a ring of radius $\sqrt{u\|N\|}$ around one and called it the computed spectrum. That ring is this essay’s law at $\tau = u$: rounding is an inexact solve with an accuracy of one unit of roundoff. And the earlier essay found GMRES finishing in two steps anyway, with second residuals between $4.9\cdot10^{-17}$ and $2.4\cdot10^{-13}$. Both observations were the same fact. The eigenvalues it could compute had moved $10^{-8}$ to $10^{-6}$ from one; the residual GMRES reached was $10^{-17}$ to $10^{-13}$. A spectrum drawn from a decomposition shows the square root; an iteration that builds the right polynomial does not.

That distinction is worth keeping in view whenever a preconditioned spectrum is shown as the argument for an iteration count. The answer that arrives when the space runs out gives the count as the degree of the minimal polynomial on the starting vector’s Krylov space. A perturbation destroys that degree at once — the exact minimal polynomial of $I + N + F$ has degree $n + m$ — and an iteration count read from the spectrum would jump from two to seventeen. What survives is that a low-degree polynomial is nearly a minimal one, and how nearly is set by the polynomial’s behaviour at the perturbed eigenvalues: here $\tau$, through a double root, rather than $\sqrt\tau$.

What eighteen small systems do not show #

The systems are dense and small, so the preconditioned matrix can be formed and decomposed; the largest has twenty unknowns. The inner solve is a fixed perturbation, which a real inner iteration is not: conjugate gradients stopped at a tolerance returns a different error for every right-hand side, the preconditioner then changes from step to step, and the outer method has to be flexible GMRES. The theory of that case — an inner accuracy that may be relaxed as the outer residual falls — is well developed, and nothing here tests it. The Schur complement is exact; with an approximate one the Jordan blocks are already split by the approximation, the eigenvalues are distinct, and this essay’s distinction between $\sqrt\tau$ and $\tau$ is replaced by the earlier essay’s spread. And every count is to a tolerance of $10^{-10}$; a looser outer tolerance moves the rule’s thresholds by the same logarithm.

Still open: a changing inner solve, an approximate Schur complement, and restarting #

An inner solve that changes. A real inner iteration returns a different error for each vector it is applied to. Under flexible GMRES the argument above no longer holds as stated, because there is no single $F$. 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.

Both blocks inexact. With the Schur complement approximated to a spread $c$ as well, the pairs are split by $c$ before $\tau$ arrives. The prediction is that the two-step residual is then the larger of the two effects — the spread’s term, independent of $\tau$, and the $c\tau$ term — so that tightening the inner solve below the spread’s level buys nothing, and the useful inner accuracy is set by the Schur approximation.

Restarting. GMRES stores every vector, and the triangular form’s count of six at $\tau = 10^{-2}$ is still small. A restart length of two — the minimal polynomial’s degree — would keep the storage of MINRES. The prediction is that GMRES(2) with the triangular preconditioner converges at every $\tau$ up to $10^{-3}$, at a rate of one factor of $c\tau$ per cycle, and fails to converge at all on the systems with the largest $c$ at $10^{-2}$.