# How to address numerical non-associativity for parallel reduction?

A parallel reduction assumes that the corresponding operation is associative. This assumption is violated for addition of floating point numbers. You might ask why I care about this. Well, it makes results less reproducible. And it gets worse when simulated annealing is used to optimize (or fit parameters) over subroutines producing such non-reproducible results.

What are the common ways to deal with this problem? What can be said about the following strategies?

• Don't care about the non-reproducibility.
• Don't use parallel reduction with floating point numbers and addition.
• Create suitably sized work-packages in a reproducible way, and do the final reduction by hand.
• Use higher precision for the addition (but not all compilers offer higher precision floating point types).
• Are you concerned about reproducibility on the same number of processes or reproducibility on a different number of processors? How much of a performance hit are you willing to take for bitwise reproducibility? Are you only interested in simulated annealing? – Jed Brown May 29 '12 at 21:31
• @JedBrown I'm concerned about the possibility to get reproducible results, for example in order to debug potential problems. It's fine for me if there is a way to reproduce the results, for example by using the same number of processors (or by just "knowing" the number of processors used originally). I would be willing to take the performance hit associated with using higher precision floating point types for the addition itself. My concrete problems were mostly related to simulated annealing and unexpected differences, but these all turned out to be real bugs . – Thomas Klimpel May 29 '12 at 21:49

A reduction implemented using MPI_Allreduce() is reproducible as long as you use the same number of processors, provided the implementation observed the following note appearing in Section 5.9.1 of the MPI-2.2 standard.

Advice to implementors. It is strongly recommended that MPI_REDUCE be implemented so that the same result be obtained whenever the function is applied on the same arguments, appearing in the same order. Note that this may prevent optimizations that take advantage of the physical location of processors. (End of advice to implementors.)

If you need to guarantee reproducibility at all costs, you can follow the guidelines in the next paragraph:

Advice to users. Some applications may not be able to ignore the non-associative nature of floating-point operations or may use user-defined operations (see Section 5.9.5) that require a special reduction order and cannot be treated as associative. Such applications should enforce the order of evaluation explicitly. For example, in the case of operations that require a strict left-to-right (or right-to-left) evaluation order, this could be done by gathering all operands at a single process (e.g., with MPI_GATHER), applying the reduction operation in the desired order (e.g., with MPI_REDUCE_LOCAL), and if needed, broadcast or scatter the result to the other processes (e.g., with MPI_BCAST). (End of advice to users.)

In the broader scheme of things, efficient algorithms for most applications take advantage of locality. Since the algorithm is really different when run on a different number of processes, it is just not practical to exactly reproduce results when run on a different number of processes. A possible exception is multigrid with damped Jacobi or polynomial (e.g. Chebyshev) smoothers, where it is possible for this simple method to perform very well.

With the same number of processes, it is often beneficial to performance to process messages in the order that they are received (e.g. using MPI_Waitany()), which introduces non-determinism. In such cases, you can implement two variants, the fast one that receives in any order and a "debug" one that receives in a static order. This requires that all underlying libraries are also written to offer this behavior.

For debugging in some cases, you can isolate part of a computation that does not offer this reproducible behavior and perform that redundantly. Depending on how the components were designed, that change may be a tiny amount of code or very intrusive.

For the most part I ditto Jed's answer. However, there is a different way out: Given the size of normal floating point numbers, you can store every number in a 4000-or-so bit fixed point number. So if you do a reduction on floating point numbers thus embedded, you get an exact calculation, no matter the associativity. (Sorry, I don't have the reference to who came up with this idea.)

You can implement a numerically stable reduction algorithm in MPI the same as you can do in serial. There may be a performance hit, of course. If you can afford to replicate the vector, just use MPI_Gather and do the numerically stable reduction in serial on the root. In some cases, you may find the performance hit is not a big deal.

Another solution is to use wide accumulators as described here. You can do this with MPI as a user-defines reduction, although it will use a lot more bandwidth.

A compromise for the above is to use compensated summation. See references “Kahan summation” for details. Higham’s “Accuracy and Stability of Numerical Algorithms” is an excellent resource on this topic.