One eigenvalue and two steps

Three eigenvalues, and two are the golden ratio preconditioned a saddle-point system

K = [ H   Aᵀ ]
    [ A   0  ]

by the block diagonal of its own two definite pieces, P = blkdiag(H, S) with S = AH⁻¹Aᵀ the Schur complement, and found that P⁻¹K has exactly three distinct eigenvalues: 1, and the two roots of λ² − λ − 1. A minimal polynomial of degree three, so three steps of a Krylov method, at every conditioning of H and A.

The block diagonal is the preconditioner that keeps the symmetry. It throws away the off-diagonal block Aᵀ, and the essay did not ask what that block was worth. There is an obvious experiment: put it back. The block upper-triangular preconditioner

P = [ H   Aᵀ ]
    [ 0   −S ]

costs one more product with Aᵀ per application and keeps everything else the same. And when S is the exact Schur complement, every eigenvalue of P⁻¹K is exactly one. Three distinct values become one.

The reading that suggests itself is that the preconditioned matrix is the identity in disguise and a Krylov method should finish in a single step. The measurement says two steps, on every system tried, and the reason the second step is necessary is the whole content of this essay.

Why every eigenvalue is one #

The triangular factorisation of K makes the claim a two-line computation. Eliminating the constraint rows gives

K = [ I       0 ] [ H   Aᵀ ]
    [ AH⁻¹    I ] [ 0   −S ]

so K = L·P with L unit lower triangular and P exactly the preconditioner. Then P⁻¹K = P⁻¹LP, which is similar to L, and L has only ones on its diagonal. Every eigenvalue of P⁻¹K is one because every eigenvalue of a unit triangular matrix is one.

Similar to L is not equal to the identity, and that is the part that matters. L minus the identity is the single block AH⁻¹ in the lower left; it is nonzero, and its square is zero because a strictly lower block matrix with one nonzero block cannot survive being multiplied by itself. So P⁻¹K = I + N with N² = 0 and N ≠ 0. On the drawn system ‖N‖ is 54.2 and the computed ‖N²‖ is 1.5·10⁻¹³, which is rounding on a matrix whose entries are of order fifty.

A matrix of the form I + N with N nilpotent of index two has (λ − 1)² as its minimal polynomial, not λ − 1. The minimal polynomial is what a Krylov method pays for, as the answer that arrives when the space runs out sets out: GMRES terminates at the degree of the minimal polynomial of the matrix restricted to the starting vector’s Krylov space. Degree two, two steps. The eigenvalue count says one and the Jordan structure says two, and the Jordan structure is the one that is asked.

In the language of that structure, the eigenvalue one has algebraic multiplicity n + m = 17 and geometric multiplicity n = 12: there are m = 5 Jordan blocks of size two hiding under a spectrum that reads as a single point.

The first step is never enough, measured #

Two steps is a prediction, and it has a failure mode worth ruling out: a starting vector that happens to lie in the null space of N would converge in one. So the claim is checked on eighteen systems — three shapes (twelve unknowns with five constraints, nine with six, fourteen with two), three conditionings of H (1, 10², 10⁴) and two of A (10 and 10⁴) — with GMRES run to a relative residual of 10⁻¹⁰.

On every one of the eighteen whose κ(H)κ(A) is at most 10⁶, GMRES takes two steps. The first step’s relative residual ranges from 0.32 to 5.9·10⁻⁹, and the second step’s from 4.9·10⁻¹⁷ to 2.4·10⁻¹³: the first step always leaves at least a thousand times more than the second, which is what a minimal polynomial of degree two looks like from inside the iteration. At the one corner beyond 10⁶ — κ(H) = κ(A) = 10⁴ on the fourteen-by-two shape — it takes three, and only because the second step’s residual lands at 2.2·10⁻¹⁰ against a tolerance of 10⁻¹⁰; a third step then finds the rounding.

The refusal on this essay is the one-step claim. It runs the drawn system on the claim that GMRES finishes in at most one step, and it is required to fail.

The size of the first step’s residual is informative in its own right. It is smallest — 1.0·10⁻⁵, 2.5·10⁻⁵ and 5.9·10⁻⁹ — on the three systems with κ(A) = 10⁴ and κ(H) = 1, and it spans eight orders across the eighteen with no single block’s conditioning ordering it. A nilpotent part does not care how large it is for the count of steps, but the first step’s partial success depends on how much of the starting residual N happens to annihilate, and that depends on the blocks.

A defective eigenvalue does not compute as a point #

The spectrum in the figure is drawn from a Francis decomposition of the explicitly formed P⁻¹K, and it is not a point at one. It is a ring.

That is not an error in the decomposition. A Jordan block of size two moves its eigenvalue by the square root of whatever perturbs it: the matrix [[1, 1], [ε, 1]] has eigenvalues 1 ± √ε. Rounding perturbs I + N by about u‖N‖, so the computed eigenvalues should scatter at a distance of about √(u‖N‖) from one — not at u, which is where a diagonalisable eigenvalue would sit. A condition number for one eigenvalue measures the same phenomenon on an eight-by-eight Jordan block, where the exponent is an eighth rather than a half.

Measured against that prediction:

Within a few per cent at the two gentle conditionings, within 22% at 10⁴, and a factor of three low at 10⁶, where the entries of N span enough orders that the perturbation rounding actually applies is no longer well described by u‖N‖. The check behind the table places each radius in a band around the prediction and requires it to be far above u, which is the statement that matters: an eigenvalue routine run on this matrix reports seventeen eigenvalues scattered about 10⁻⁷ from one — eight or nine orders of magnitude above the rounding a diagonalisable eigenvalue would show — and a reader who took that as a spread would be reading the Jordan structure as a clustering.

That also answers a question the triangular preconditioner raises for anyone who has read the spectrum that predicts nothing. GMRES convergence on a nonnormal matrix is not governed by its eigenvalues in general — any curve is compatible with any spectrum — and P⁻¹K is as nonnormal as a matrix with one eigenvalue can be. What rescues the prediction here is not the spectrum but the minimal polynomial, which is exact.

An approximate Schur complement: one copy of every value against two #

The exact S is not what anyone uses. It costs m solves with H and a dense m × m factorisation, and the point of a block preconditioner is to replace it by something cheap. So the next measurement uses an approximation Ŝ whose quality is a single number: the generalised eigenvalues of the pencil (S, Ŝ) are spread geometrically over [1/c, c], so c = 1 is exact and c = 10 is a Schur approximation wrong by an order of magnitude in both directions.

The two preconditioners then leave recognisably different spectra. With the triangular form, P⁻¹K has eigenvalue one on the n-part and each generalised eigenvalue ν of (S, Ŝ) once. With the diagonal form, each ν produces two eigenvalues, (1 ± √(1 + 4ν))/2 — the golden-ratio formula of the earlier essay with ν in place of 1.

Five distinct values against eleven, and six steps against eleven. The triangular count is one more than its number of values, and the reason is the same as before: the eigenvalue one is still defective, carrying the nilpotent coupling from the n-part, so it contributes a factor (λ − 1)² to the minimal polynomial rather than λ − 1. The diagonal preconditioned matrix is symmetric after the usual congruence, has no Jordan blocks, and its count is exactly its number of values.

While the spectrum has few enough distinct values for a polynomial to interpolate all of them, each method terminates at its count, and the counts do not depend on how wide the spread is — only on how many values there are.

At c = 10 the step counts are exactly those at c = 1.5, 2 and 4. A Schur approximation that is wrong by an order of magnitude costs nothing over one wrong by fifty per cent, on a problem this small, because the count of values is all either method is paying for. That is a finite-dimensional fact and it would be dishonest to generalise it, which is why the next measurement is larger.

When there are too many values to count #

With thirty unknowns and fifteen constraints the approximation leaves fifteen generalised eigenvalues, and the counts start to depend on the spread: 14 triangular steps against 27 diagonal ones at c = 1.5, 15 against 30 at c = 2, 16 against 31 at c = 4, 16 against 37 at c = 10.

Pushed to sixty unknowns and forty constraints at a looser tolerance of 10⁻⁶, where the counts are set by how widely the values spread rather than by how many there are, the ratio of triangular to diagonal steps is 0.54 at c = 1.2, 0.48 at c = 2 and 0.49 at c = 4. The halving survives the change of regime. In the counting regime it is exact — one copy against two — and in the spreading regime it is approximate, and the reason it survives is that the diagonal form’s spectrum is two intervals, one on each side of zero, while the triangular form’s is one interval on the positive side; a polynomial that must be small on two intervals separated by the origin needs about twice the degree of one that must be small on one.

With the approximations a code would actually use #

The spread construction is the right instrument for a measurement, because its error is a single number that can be swept, and it is the wrong one for a verdict, because nobody builds an approximation whose generalised eigenvalues are spread geometrically on purpose. So the two assemblies are run again with the two approximations the golden-ratio essay measured: Ŝ = A diag(H)⁻¹Aᵀ, which replaces H by its diagonal inside the complement, and a scaled AAᵀ, which replaces it by a multiple of the identity. Neither states its own error.

On twelve unknowns and five constraints both give the triangular form 6 GMRES steps and the diagonal form 11 MINRES steps — exactly the counts every spread from 1.5 to 10 gave. On thirty unknowns and fifteen constraints both give 16 against 31, exactly the counts of the spread-4 figure above. The ratio is 0.545 on the small system and 0.516 on the larger one, for both approximations.

That coincidence is the counting regime announcing itself, and it is worth being clear that it is a coincidence of counts rather than of quality. The diagonal approximation and the scaled identity are different matrices with different generalised eigenvalues; on the small system each still leaves five distinct values, and a method that terminates at the number of values cannot tell five good values from five bad ones. The golden-ratio essay’s figure makes the same point from the MINRES side, where the diagonal and identity approximations both cost 11 steps against the exact complement’s 3.

What the practical approximations do establish is the part a reader would want before choosing: the halving is not an artefact of a construction built to have it. With the approximations a saddle-point code actually reaches for, the triangular assembly still needs about half the steps of the diagonal one, on both sizes, and the check behind it requires the ratio to stay under three quarters on both.

What they do not establish is how the two assemblies compare once the step counts are set by spread rather than by count with a practical approximation — the regime the sixty-by-forty measurement reached with a stated spread. A diagonal approximation to the Schur complement of a discretised flow problem has an error that grows with the mesh, and there the ratio could move; that is a measurement on a larger family of problems than the ones measured here.

The price is the recurrence #

Halving the step count is not halving the work, and the cost ledger is where the two preconditioners genuinely part company.

Per application, the triangular preconditioner costs one solve with H, one solve with Ŝ, and one product with Aᵀ; the diagonal one costs the two solves. The extra product is cheap — A is the constraint matrix, usually sparse — so per step the two are close.

Per step of the Krylov method, they are not close. The triangular P is not symmetric, so P⁻¹K is not symmetric in any inner product the preconditioner supplies, and MINRES — which needs a symmetric operator and a symmetric positive definite preconditioner — cannot be run. GMRES can, and GMRES keeps every basis vector it builds and orthogonalises each new one against all of them: step k costs k inner products and step k stores k vectors. An orthogonalisation nobody calls one is the collection’s account of what that orthogonalisation costs and what goes wrong in it. MINRES, on the symmetric diagonal form, runs on a three-term recurrence and stores three vectors however long it runs — the structure one sequence and two recurrences is about.

So at sixteen steps against thirty-one, GMRES has done about 16·17/2 = 136 inner products and stores seventeen vectors; MINRES has done a fixed handful per step, about a hundred in total, and stores three. In orthogonalisation work the two are level at this size, and in memory MINRES is five times cheaper. On a problem needing a hundred triangular steps GMRES would store a hundred vectors and do five thousand inner products against MINRES’s two hundred steps and constant storage; restarted GMRES would recover the memory and give back some of the steps, and how many is a property of the spectrum’s spread that has not been measured here.

The clean way to say it: the triangular preconditioner buys half the steps with the recurrence, and whether that is a bargain depends on how many steps there were to halve. When the Schur approximation is good enough that both counts are small, the halving is nearly free. When it is poor enough that both are large, the recurrence is worth more than the steps.

What this adds to the two essays before it #

The golden-ratio essay established that a Schur complement preconditioner is a statement about which values the preconditioned matrix can have, and that the affordable approximations are measured by how far they move those values. A preconditioner that need not know the constraint kept the constraint exactly and approximated the objective instead, and paid for it with an indefinite preconditioner.

This essay adds the structural axis neither of them had. The same two blocks, H and an approximation of S, can be assembled into a preconditioner two ways, and the assembly decides three things at once: whether each approximate eigenvalue appears once or twice, whether the eigenvalue one is defective, and whether the Krylov method has a short recurrence. The triangular assembly wins the first, pays for it in the second — one extra step that no eigenvalue count predicts — and loses the third.

It is also the collection’s clearest case of an eigenvalue count misleading in the direction nobody expects. The eigenvalues that are not there and the non-normality field more broadly warn that eigenvalues can make a matrix look better conditioned than it is. Here a single eigenvalue makes a matrix look like the identity, and the matrix is two steps away from it, with a nilpotent part of norm 4,860 at κ(H) = 10⁴ that no spectral plot shows.

What has not been measured #

The step counts use exact solves with H and Ŝ inside each application. A practical block preconditioner replaces those by a multigrid cycle or an incomplete factorisation, and an inexact inner solve breaks the exact nilpotency: N² is no longer zero, only small, and the second step’s residual is no longer at rounding. How the two-step structure degrades under an inexact solve with H — whether gracefully, like a perturbed Jordan block, or by adding steps in proportion to the inner tolerance — is the obvious next measurement and is not made.

The problems are small and dense. The counting regime, where step counts are set by the number of distinct values, is a property of a problem with few constraints; the spreading regime at sixty unknowns is closer to practice and was measured at one looser tolerance only.

And the comparison is GMRES against MINRES in steps and in stored vectors. A lower-triangular variant, a right-preconditioned form, and the constraint preconditioner of the previous essay assembled triangularly are each a different matrix with the same eigenvalues and different Jordan structure; none is measured.

Where this goes from here #

Making the Schur complement easy instead of approximating it. Every approximation above treats S as fixed and asks how well Ŝ matches it. The augmented Lagrangian approach changes the problem instead: add γAᵀA to H, which does not move the solution, and the new Schur complement tends to I/γ, which is trivial to apply. Where the augmentation puts the cost measures how close to the golden-ratio spectrum that brings the diagonal preconditioner, and what it does to the solve with the augmented H — which is where the cost, and the digits, go.

An inexact inner solve. 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.

And restarting. The triangular form’s advantage is in steps and its disadvantage is in storage, and restarted GMRES trades one for the other. The restart length at which the triangular preconditioner stops being worth its recurrence is a single number per problem, and it has not been measured.