A stopping test is a race

Everything the field has measured so far moves a number in its last digits. This is where the last digits stop being where the answer is.

Iteration counts from 674 to 690, eleven of the thirteen distinct. Every run converged. Every answer is correct to the accuracy that was asked for — the forward errors run from 6.3·10⁻¹¹ to 1.1·10⁻¹⁰, a factor of 1.72, all of them comfortably under the tolerance.

What differs is an integer, it is 2.4% of the work, and it is a property of the machine.

That number is the field’s first result that a scheduler would notice. A disagreement in the eighth digit of a residual is invisible to everything outside the arithmetic; a solve that takes a different number of synchronising steps changes the shape of a parallel program’s execution, and does so for a reason that is not in the program.

Where the reductions are #

The instrument is worth describing because it is unusually honest about what it changes.

The conjugate gradient method performs, per step: one matrix–vector product, three vector updates, and two inner products — rᵀz for the residual update and pᵀAp for the step length. Those two are where a parallel implementation synchronises, they are the all-reduces the communication field counts in the message and the word, and they are the only places in the method where every worker has to agree about a number.

So the experiment replaces the arithmetic context’s inner product with a p-way partitioned reduction and changes nothing else. Addition, multiplication and division stay binary64. The matrix–vector product stays sequential. The method is this site’s own cg, unmodified, and the two ar.dot calls inside it are exactly the two quantities a distributed implementation computes differently.

Those two numbers are also where every scalar in the method comes from: α = rz/pAp is the step length and β = rz’/rz is the direction update. A disagreement of a few hundred ulps in either produces a slightly different iterate, which produces a slightly different residual, which produces a slightly different next step.

Why a few hundred ulps becomes sixteen iterations #

The chain is short and each link is measurable.

The reductions disagree by about κ_sum·u. For pᵀAp the summation condition number is a few hundred — the vector that hides it measures 4.1·10² for a curvature — so the disagreement is around 10⁻¹³ relative. Small.

The iterates diverge geometrically. Conjugate gradients is not a contraction in the way a stationary iteration is; two runs that differ at step k differ by more at step k + 1, because the search directions are built from the residuals and the residuals are built from the directions. By step 600 two runs that started 10⁻¹³ apart are apart by whatever the method’s own error level is.

And the stopping test is a threshold crossing. The residual falls, crosses 10⁻¹⁰ once, and the step at which it does is decided by where the curve happens to be relative to the line. Two curves that differ by a per cent cross a horizontal line at slightly different places, and the difference is quantised to whole steps.

So the mechanism turns a relative 10⁻¹³ into a difference of sixteen steps by two amplifications: a geometric one over hundreds of iterations, and a quantisation at the end.

The part that is not the machine’s fault #

An important qualification, because this figure is easy to over-read.

Conjugate gradients loses its finite-termination property in floating point for reasons that have nothing to do with reductions. The residual vectors stop being orthogonal, the method takes more steps than it should, and how many more is exquisitely sensitive to everything — which is an orthogonalisation nobody calls one’s subject, and the reason the rate the condition number predicts has to be measured against a method that does not terminate when the theory says it will.

So the sensitivity being exploited here is the method’s, not the machine’s. The machine supplies a perturbation of 10⁻¹³ and the method supplies the amplification. A different iterative method with a less chaotic trajectory would show a smaller spread from the same perturbation.

That does not weaken the finding — the perturbation is real, unavoidable and outside the program’s control — but it locates it correctly. The reduction is the cause and the method’s chaos is the amplifier, in exactly the sense the site’s identity separates a backward error from a condition number.

What the spread is a function of #

The drag says it plainly: at κ = 10² every partitioning takes the same number of steps and the field has nothing to show. The spread appears with the conditioning.

That is not a coincidence, it is the same κ twice. A badly conditioned matrix gives a residual whose terms cancel more — so the reductions disagree by more — and gives a method that takes many more steps for the disagreement to be amplified over. Both factors point the same way, and the refusal in this essay’s refutation is the check: on a matrix with κ = 4 the assertion that any problem shows a spread of iteration counts is fed the case and must reject, which it does at 21 steps for every partitioning.

So the operational reading is: a well-conditioned solve is reproducible for the same reason it is accurate, and a badly conditioned one is neither. There is no separate diagnosis to make.

What this costs, and what it does not #

The 2.4% is worth putting in proportion, because it is easy to read the figure as alarming and easier to read it as trivial, and it is neither.

It is not an accuracy problem. Every answer is right. A user who wanted a solution to 10⁻¹⁰ got one, thirteen times.

It is a cost problem, and a small one. Sixteen extra steps out of 674 is 2.4% of a solve. Nobody’s budget turns on that.

It is a scheduling problem, and that one is not small. In a parallel program every worker has to take the same number of iterations, because the reductions are synchronisation points. If a solve is one stage of a pipeline and its length varies by 2.4% between runs, every downstream stage’s timing moves with it, and a system with a deadline has to budget for the worst case rather than the mean.

And it is a reproducibility problem, which is the one that bites. Two runs that took different numbers of steps returned different answers — not just different in the last digits, but arrived at by different trajectories. A test suite comparing an answer with a stored one is comparing across that, which is what a regression test can ask for’s subject.

The residual the test reads is not the residual #

One further mechanism deserves naming, because it interacts with this one and the two are easy to confuse.

Conjugate gradients does not recompute b − Ax at every step. It updates the residual by a recurrence, r ← r − αAp, which costs nothing extra and is exact in exact arithmetic. In floating point the recurrence’s residual drifts away from the true one, and the stopping test is written in the recurrence’s — which is the residual the method reports, and the gap it measures is a quantity no stopping test can see.

So the number the test compares with a tolerance is doubly removed from what anyone wants: it is the residual of the recurrence rather than of the matrix, and it is computed by a reduction whose value depends on the partition. The two are independent and they compose. The drift is a systematic quantity bounded by a multiple of u·max‖xⱼ‖; the disagreement is a random quantity of size κ_sum·u. A run that stops early because its recurrence has drifted low is a different failure from a run that stops early because its reduction landed low, and only the first has a repair inside the method — residual replacement, one extra product every few steps.

Measured together they are why two runs of the same solve return answers a factor of 1.72 apart rather than a factor of 1.0000001. Neither mechanism alone would do it.

Two repairs, and one that looks like a repair #

A fixed iteration count. Deciding in advance to run 700 steps and not to test at all removes the race completely: every machine does the same work in the same order and the only variation left is the last digits of the answer. It costs the difference between 700 and however many the problem needs, and it requires knowing the number in advance, which a rate that is known in advance says is available for the model problem and for very little else.

An order-independent reduction. Pre-rounded accumulation on the two inner products makes both of them functions of the vectors rather than of the partition, so every machine follows the same trajectory and stops at the same step. What determinism costs prices it: on a dense matrix, 5% of a step for exact accumulation and 1% for pre-rounding; on a five-point stencil, 122% and 22%. The sparse case is where the trade is real, and it is the case where the reductions are also the smallest part of the answer’s accuracy.

And the one that looks like a repair: tightening the tolerance. It does not work, and the measurement is the next essay’s: the tolerance that buys no agreement runs the same sweep at four tolerances and finds the ratio between the best and worst run at 1.34, 1.48, 1.71 and 1.17 while the accuracy improves by a factor of 1.5 million.

What a caller sees, and what they conclude #

Worth walking through, because the observed symptom is not the one this essay describes and the wrong conclusion is easy to reach.

A caller does not see thirteen partition counts. They see a job that took 4 minutes 12 on Tuesday and 4 minutes 18 on Wednesday, on the same input, with the same binary. If they look at the iteration count they see 674 and 690. Nothing in either run reports a partitioning, because the partitioning is chosen inside a runtime that has no reason to mention it.

The available explanations, in the order people reach for them:

The machine was busier. Reasonable, and testable by looking at the wall-clock per iteration rather than the total — which will be flat, because the iteration count moved rather than the speed.

The input changed. Testable, and usually the first thing checked.

The library was updated. Sometimes true, and the one explanation that is both plausible and verifiable — which is why it absorbs a great deal of investigation that belongs elsewhere.

Floating point is non-deterministic. Reached last, stated as folklore, and treated as a non-explanation because it comes without a number.

The contribution of this field is the number. The two reductions in this solve have summation condition numbers of about 4·10² and 1, so their answers differ between partitionings by about 10⁻¹³, which the method amplifies over 674 steps into a residual curve that crosses the tolerance up to sixteen steps away from where it did yesterday. That is a statement somebody can act on: it says which quantity to make order-independent, what it will cost, and what to expect if they do nothing.

What a stopping test is, restated #

The essay’s title is not a metaphor and it is worth ending on the literal reading.

A stopping test compares a computed number with a constant. The computed number is a reduction whose value depends on how the work was divided. The constant does not. So the comparison’s outcome depends on a division of work that no part of the program chose, and the first step at which the comparison comes out true is the outcome of a race between the residual’s fall and the arithmetic’s disagreement.

Every iterative method on this site has one. The essay that measures what a tolerance actually obtains — a tolerance that reads its own residual — asks what accuracy a given stopping value buys, and answers with a number that this field now qualifies: the accuracy obtained is what that essay says, and which step obtains it is what this one says.

The same race in three other methods #

Conjugate gradients is the clean case because it has exactly two reductions and both matter. The shape recurs wherever a method’s control flow reads a computed number.

GMRES. Every step orthogonalises against every previous vector, so the reductions are a Gram–Schmidt’s worth rather than two, and the residual is available from the Hessenberg least-squares problem rather than from a recurrence. More reductions, at inner products whose κ is the orthogonality check’s — which is at the top of the census — and the same threshold crossing at the end.

Lanczos with a convergence test on the Ritz values. The test asks whether a Ritz value has settled, which is a comparison of two computed numbers rather than of one number with a constant, and the arithmetic that produces them is what an orthogonalisation nobody calls one shows is already fragile. Two builds converge different eigenvalues at different steps, which changes what the method returns rather than only when.

An inexact Newton method. The inner solve’s stopping test decides how accurate a step is, and the outer iteration reads the step. A difference of one inner iteration is a different step, which is a different linearisation point, which is a different problem for the next inner solve — the accuracy that is thrown away is the essay about how much of that inner accuracy survives, and the answer there is not much, which is the reassuring half: an outer loop that discards its inner accuracy also discards the machine’s disagreement about it.

That last point generalises, and it is the field’s one genuinely comforting result. A loop that recomputes its residual from the matrix at every step forgets what the previous step got wrong, so the disagreement does not accumulate across an outer iteration. It accumulates within one, which is what the sixteen steps here are.

At other settings #