Floating point: why two runs of one model disagree

Two engineers run the same model, from the same input file, with the same version of the same program, and get factors of safety of 1.4173 and 1.4172. Nothing is broken. The difference is real, it is reproducible on each machine, and it comes from the one property of floating-point arithmetic that most numerical code is written as though it did not have.

Every number in this article is recomputed from its inputs by a script that runs on each build of this site, in IEEE 754 double precision — the same arithmetic being described.

Addition is not associativePermalink to “Addition is not associative”#

The familiar demonstration is that

$$0.1 + 0.2 = 0.30000000000000004$$

which is not 0.3. The residual is $5.5511151231257827 \times 10^{-17}$, and it is there because neither 0.1 nor 0.2 is representable in binary; the sum is the correctly rounded result of adding two values that were already approximations.

That is a representation problem and it is well known. The associativity problem is separate and worse. Take four numbers and add them left to right:

$$10^{16} + 1 - 10^{16} + 1$$

The first partial sum is $10^{16} + 1$, and the spacing between representable doubles near $10^{16}$ is 2, so that value cannot be stored. It rounds back to $10^{16}$, the 1 is gone, and the running total finishes at 1. Reorder the same four terms so the two large ones cancel first and the total is 2, which is the exact answer.

No rounding rule was violated. Every individual operation returned the correctly rounded result of its two arguments. The sum is order-dependent because rounding is, and there is no implementation that avoids it.

What this looks like in an assembly loopPermalink to “What this looks like in an assembly loop”#

The four-term example is contrived. Here is one that is not.

A nodal force is assembled by adding contributions from every element that touches the node. Take one large contribution of 45000000 and a thousand small ones of 0.37, whose exact total is 45000370:

The two naive orders differ by $2.682209 \times 10^{-6}$, which is a relative error of $5.960415 \times 10^{-14}$ and, more usefully stated, exactly 360 units in the last place of the result. Each small addition lost part of itself to rounding against a running total far larger than it was, and the losses accumulated in one direction.

Smallest-first is better because each addition is between operands of comparable size. Kahan summation is better still: it keeps a running compensation term for what each addition lost and adds it back, recovering the exact value here at the cost of three extra operations per element.

Where the order comes from, and why it changesPermalink to “Where the order comes from, and why it changes”#

An assembly loop does not usually choose its summation order deliberately. The order falls out of things that are not thought of as numerical decisions at all:

• Thread count. A parallel reduction sums partial results in whatever order the threads finish, unless it is written not to. Two runs on the same machine with the same binary can differ.
• Element numbering. Renumbering a mesh for bandwidth reduction changes the order in which contributions arrive at each node.
• Compiler flags. -ffast-math and its relatives explicitly permit the compiler to reassociate arithmetic, which is to say they permit it to change the answer.
• Fused multiply-add. An FMA computes $a \times b + c$ with a single rounding instead of two. It is more accurate, and it is a different result. Whether the compiler emits one depends on the target architecture.
• Library versions. A BLAS that blocks a dot product differently between releases sums in a different order.

None of these appear in a model description, and all of them move the last digits.

Why it matters more than the size suggestsPermalink to “Why it matters more than the size suggests”#

360 units in the last place is a relative error of $6 \times 10^{-14}$, which is negligible against every source of uncertainty in a geotechnical problem. If the question is “will this slope stand”, the answer is unaffected.

But two other questions are affected, and both are questions about the software rather than about the ground.

Verification. A regression test compares today’s output against a stored reference. If the program is not deterministic, the test must be given a tolerance, and a tolerance wide enough to absorb non-determinism is also wide enough to hide a genuine defect. The value of a benchmark suite depends on the difference between “passes” and “passes to within something”.

Nonlinear amplification. In a plastic analysis a stress point either yields or does not, and that is a branch on a comparison. A perturbation of $10^{-14}$ in a stress can flip a yield check, which changes which region of the yield surface the return targets, which changes the tangent, which changes how many Newton iterations are taken, which changes the increment. Near a limit load, where the solution is genuinely sensitive, a difference in the last place can become a difference in the third.

What to doPermalink to “What to do”#

Determinism in floating point is achievable, and it is achieved by construction rather than by care:

• Fix the reduction order. Sum into per-thread buckets in a fixed layout and combine them in a fixed sequence, rather than letting completion order decide.
• Reject reassociation. Do not compile numerical kernels with fast-math. The performance it buys is small next to a result that cannot be reproduced.
• Pin the FMA decision. Contract or do not, explicitly, and record which.
• Compensate where it is cheap. Kahan summation on assembly loops costs a few percent of a phase that is not usually the bottleneck.
• Test bit-for-bit. A regression suite that asserts exact equality is asserting something much stronger than “close enough”, and it is the only kind of test that catches a reassociation introduced by an innocuous refactor.

The underlying position is the one that runs through everything on this site: a result that cannot be reproduced is not yet a result. Floating-point arithmetic is deterministic — every operation has one correctly rounded answer. What is not deterministic is a program that declines to say what order it does things in.