Where the drift lands

A preconditioner is the expensive object in an iterative solve and the obvious one to keep. Every code that solves a sequence keeps one, every one of them has a rule for when to throw it away, and the rule is almost always the same: rebuild when the matrix has changed by more than some fraction of itself.

Written down, that is a threshold on ‖A_t − A₀‖ / ‖A₀‖. It is a perfectly reasonable-looking quantity — it is relative, it is cheap, and it goes to zero when nothing has happened.

It is also the wrong quantity, and the size of the mistake is the condition number.

The experiment #

The preconditioner is the exact Cholesky factorisation of the first matrix, which is not what anybody uses and is exactly what this measurement needs: with it, M⁻¹A₀ is the identity to the last bit and conjugate gradients solves the first member in a single step. So every iteration counted afterwards was bought by the drift and by nothing else — which is the same construction the sequence field uses for its known root, and for the same reason.

The drift is then built in the spectrum rather than added and measured afterwards. A₀ is constructed with a known eigendecomposition at a prescribed condition number, and the perturbation is

E = s · Σ vᵢ vᵢᵀ

summed over the forty smallest eigenvalues or the forty largest, with s chosen so that the Frobenius norm of E is the same fraction of ‖A₀‖ in both cases. The two drifts are the same size in the norm the rule is written in, to fourteen digits, and that equality is asserted rather than assumed.

The measurement #

At a relative drift of 10⁻³ — a tenth of a per cent, a number most rules would not even react to — the two runs cost 19 iterations and 5. At 10⁻² they cost 38 and 9. The rule cannot distinguish them, because in the rule’s own quantity they are identical.

Why, which is one line #

With M the exact factorisation of A₀ and A_t = A₀ + E,

M⁻¹A_t = I + M⁻¹E

and if E is built from A₀’s own eigenvectors with weight wᵢ on the ith, the eigenvalues of the preconditioned matrix are exactly

1 + wᵢ / λᵢ .

The drift is divided by the eigenvalue it lands on. A perturbation sitting on λ_min is amplified by 1/λ_min; the same perturbation on λ_max is divided by λ_max; and ‖E‖ is the same number in both cases because a Frobenius norm does not know where in the spectrum its mass is.

So the quantity that predicts the cost is max(1 + wᵢ/λᵢ) over min(1 + wᵢ/λᵢ), and at a relative drift of 10⁻³ it is 4.47 on the small end and 1.0327 on the large one. Both are available before either solve starts, from the construction rather than from an eigensolve of the product.

The predictor orders and over-predicts #

The predictor gets the ordering right at every drift measured, which is the finding. What it does not get right is the size, and saying so is the second half of an honest measurement.

The conjugate gradient bound applied to a preconditioned condition number of 3.47·10³ asks for 699 iterations. The run takes 58. The bound over-predicts by a factor of twelve, and the fitted slope of iterations against predicted conditioning is 0.176 rather than the bound’s 0.5.

The reason is the one this collection already has an essay about: a stale preconditioner does not spread the spectrum, it splits it. Half of M⁻¹A_t’s eigenvalues are exactly 1, because the drift was placed on the other half, and a method that minimises a polynomial over the spectrum spends nothing on a cluster it has already annihilated. The condition number is a summary of two numbers and the method is reading the whole distribution.

And the over-prediction is not a factor of twelve — it is a factor that grows.

The predicted conditioning grows by two thousand and the actual by three and a half. Across κ = 10², 10³, 10⁴, 10⁵ and 10⁶ the predictor reads 1.48, 35.67, 315.45 and 2,911.41 while the conditioning the method actually meets reads 1.0463, 1.3271, 1.9245 and 3.6667. The over-prediction is therefore 1.4×, 27×, 164× and 794× — not a constant twelve but a factor that itself grows by nearly three orders.

The ordering survives all of it, which is what makes the predictor worth keeping: the small end costs more than the large end at every one of the five conditionings, by 1.7×, 4.2×, 5.3× and 5.1×. So the predictor answers which correctly and is wrong about how much by an amount that depends on the conditioning — and a caller who calibrated the correction at κ = 10² would be out by a factor of five hundred at 10⁶.

So: where the drift lands decides the cost, and the standard bound built on the same quantity is loose by an order of magnitude at the conditioning the essay measured, and by nearly three at the top of the range. Both statements are about the same figure and neither is a correction to the other.

The object a real code holds #

The exact Cholesky above isolates the phenomenon. Nobody keeps one. What a code keeps is an incomplete factorisation, and it ages differently in a way worth drawing.

The drift is the one this collection already has a field about: the five-point operator’s two directions stop being equally weighted, ε running from 1 down to 0.033 over nineteen members, with the sparsity pattern unchanged — so the incomplete factorisation from the first member stays applicable for the whole run, which is exactly the situation in which it gets kept.

The kept factorisation goes

18 21 24 26 29 32 34 36 37 38 38 38 40 42 44 46 48 50 52

iterations. A factorisation rebuilt at every member goes

18 18 18 18 18 17 17 17 17 16 16 16 15 15 15 15 14 14 13

in the opposite direction. The problem is getting easier — an anisotropic operator is a friendlier problem for a factorisation that knows about the anisotropy — and the kept preconditioner is getting worse at it. The two start at the same 18 by construction and never meet again.

The failure mode nobody sees #

The incomplete case is the one worth worrying about in practice, and not because it is worse. It is milder, and that is the problem.

A stale factorisation stops converging — a cliff, an error, something in a log. A stale preconditioner never stops converging. It costs 18 iterations, then 30, then 52, and every member returns a correct answer to the tolerance it was asked for. Nothing in the output says that the run is now taking three times as long as it needs to, because nothing in the run knows what it would have taken with a fresh preconditioner.

That is a class of defect this collection keeps meeting: the symptom is absence. There is no wrong number to find, no assertion to trip, and no gate that could be written on the output of a single member. The only way to see it is to run the comparison, and running the comparison means building the fresh preconditioner, which is the thing being avoided.

Two routes to the same number #

The predicted conditioning is computed here from the construction — the weights wᵢ were chosen, the eigenvalues λᵢ are known, and 1 + wᵢ/λᵢ is arithmetic on two lists. That is the cheap route and it is only available because the experiment built the drift on purpose.

It would be worth nothing if it did not agree with the expensive route, so the agreement is the check. Forming M⁻¹A_t explicitly and taking the ratio of its extreme eigenvalues gives the same number, and it must: the two are the same quantity computed from two different objects, one of them never assembled. The collection’s standing habit is that a number arrived at by one route has been wrong every time it has been published, and this is the cheapest available second route.

What the agreement licenses is the sentence in the figure’s badge: the conditioning of the preconditioned matrix is knowable before the solve, from quantities a code that built the preconditioner already has. It is not knowable cheaply for a general drift — nobody has A₀’s eigenvectors — but the experiment is not about what is cheap. It is about which quantity the answer depends on, and the answer depends on that one and not on the norm.

What the rule should be written in #

If not ‖E‖/‖A₀‖, then what?

The honest answer from this measurement is ‖A₀⁻¹E‖, or anything that behaves like it: the drift measured in the units the preconditioner works in rather than the units the matrix is written in. That quantity is 4.47 and 0.033 for the two drifts above, which is a factor of 135 between two things the naive rule calls equal.

Estimating it is not free — it is a norm of a product with an inverse, which is the object the condition-estimation essays are about — but it is not expensive either, and this collection already has the machinery: an estimate of ‖A₀⁻¹E‖ costs a handful of solves with a factorisation that is already in memory.

The cheaper answer, and the one the essay on what a rebuild is worth is about, is to stop trying to predict and read the iteration count instead. That has its own failure — the count moves with the right-hand side too — and the two essays are best read together.

A minority of the drift sets the price #

The table above places the whole drift at one end or the other, which makes the two rows the extremes. The interesting question is what a mixture does, because no real drift is confined to half a spectrum, and the answer decides whether the norm rule is imprecise or wrong.

Hold the total relative drift at 10⁻² and move energy from the large half to the small one. Every row below has ‖E‖/‖A₀‖ equal to the same 10⁻² to fourteen digits, so the rule under test calls all of them the same problem:

with none of the energy on the small end the solve takes 9 iterations, and with a millionth of it, 9. With a ten-thousandth, still 9. With a thousandth, 14. With one per cent, 21. With a tenth, 31. With half, 40 — and with all of it, 40.

One per cent of the drift’s energy costs 21 of the 40 iterations, and half of it costs the whole penalty. The relationship saturates, and it saturates early.

The mechanism is one word. The quantity that decides the cost is max(1 + wᵢ/λᵢ), and a maximum is set by the worst-placed component. A norm is an average over all of them. So the two disagree whenever one component is much worse placed than the rest, which is not an edge case — it is what “structured” means, and a drift with any structure at all has its worst component somewhere in particular.

That changes the shape of the criticism. The section above says the norm rule cannot tell which end a drift is on. What the mixture says is stronger: a drift can be 99 per cent well-placed and still cost most of the badly-placed price. A rule watching ‖E‖/‖A₀‖ is not looking at a blurred version of the right quantity; it is looking at a statistic that a one-per-cent minority of the data can move by a factor of twenty-five while leaving the statistic itself unchanged to fourteen digits.

It also explains why the failure is so hard to notice in a real code. A drift that is mostly harmless still costs, so the rule’s predictions are not wildly wrong on most members — they are wrong by the amount that the small unlucky component contributes, which varies from member to member with no pattern anybody would spot. The run is slower than it should be, by a different factor each time, and nothing in the output says so.

What would have to be true for the norm rule to work #

It is worth stating the condition under which the standard rule of thumb is right, because it is not empty and because naming it is more useful than saying the rule is wrong.

A threshold on ‖E‖/‖A₀‖ predicts the cost correctly when the drift is spread evenly over the spectrum — when the perturbation has about as much weight on the small eigenvalues as anywhere else. Then max(1 + wᵢ/λᵢ) is proportional to ‖E‖/λ_min, the condition number enters once as a constant, and a threshold calibrated on one problem transfers to another of similar conditioning.

That is a real situation. A drift that is genuinely random with respect to the operator’s eigenvectors — a perturbation from measurement noise, say, or from a rounding — looks like that, and for it the rule is a reasonable proxy.

What is not like that is any drift with structure, which is most of them: a time step that changes one coefficient, a load that acts on part of a domain, an operator whose anisotropy grows. Those have their weight somewhere in particular, and where they have it is exactly the question the norm discards. So the rule is not a bad approximation to the right quantity; it is the right quantity multiplied by an unknown that runs from 1 to the condition number, and nothing in the norm says which.

What the right rule would cost #

Naming a better quantity is easy and is worth nothing on its own, so it is worth being concrete about what a code would have to do to use ‖A₀⁻¹E‖ and whether it could afford to.

The factorisation of A₀ is in memory — that is the whole premise of the situation — so applying A₀⁻¹ to a vector is a pair of triangular solves. A one-norm estimator of ‖A₀⁻¹E‖ needs a handful of products with A₀⁻¹E and its transpose, which is four or five of those pairs plus four or five products with E. On the 80 × 80 problem here that is under a tenth of one preconditioned iteration, and on a sparse operator it is a few matrix–vector products.

So the right quantity is affordable, and it is affordable because the preconditioner is being kept: the object that makes the estimate cheap is the object whose lifetime is in question. That is a pleasant arrangement and it is worth stating, because the usual objection to a better rule is that it needs something the code does not have.

Two caveats, both real. E has to be formed, or at least applied, and a code that assembles a new matrix each member has it while one that only ever applies an operator does not. And the estimator is an estimator: the condition-estimation essays measure how far a one-norm estimate can sit from the truth, and a factor of two or three there is harmless for a threshold and would not be for a prediction. What is being decided is whether to rebuild, which is a comparison against a number somebody chose anyway — so an estimate good to a factor of three is enough, against a norm rule that is wrong by a factor of the condition number.

The refusal #

The claim under test is that a relative drift threshold is a usable rule. It is fed the case where it happens to be right, and the case is chosen to be one where being right means nothing.

At a relative drift of 10⁻⁶, the small end costs 4 iterations and the large end costs 2. The assertion that the small end costs more than two and a half times the large one is fed those numbers and fails — correctly, because at that drift neither placement costs anything worth measuring and there is no distinction to make. The threshold rule is right there, in the sense that it makes no error; it is right the way a rule that always says “do nothing” is right on a problem where nothing needs doing.

That is what the refusal is for: to stop the main measurement from being read as a claim that the norm is always wrong. It is wrong where the answer matters and correct where it does not, which is the least useful pattern a rule can have.

One more thing the norm cannot see #

There is a second quantity hiding in the incomplete-factorisation figure, and it is worth naming because it is the same mistake in a different disguise.

The drift there is an anisotropy: the operator’s two directions stop being equally weighted. In the Frobenius norm that is a large change — by the last member the matrix is a substantial distance from the first — and a norm-based rule would have rebuilt many members earlier. It would have been right to, but not for the reason it thought: the rule would have fired because the matrix moved, and the correct reason to fire is that the smoothing directions the incomplete factorisation encodes have stopped being the directions the operator is stiff in.

This collection has a whole field about that distinction. Line relaxation along the strong direction and coarsening only along it are both statements about where in the problem the difficulty is, and neither of them is a statement about how far the matrix has moved. A preconditioner is a claim about the structure of the difficulty, and it goes stale when the structure moves rather than when the entries do.