mapreduce: use simpler reduction order in some cases

[nlw0]

I'm not sure this could be related to #43501 or not. I just noticed this code for the first time:

julia/base/reduce.jl

Line 258 in 6106e6c

v1 = mapreduce_impl(f, op, A, ifirst, imid, blksize)

Are there any references about this implementation? When could we expect this doubly-recursive approach to be preferable to a linear traversal?

[simeonschaub]

See the discussion in #4039 for example. For linear traversal just use mapfoldl.

[nlw0]

It's about accuracy, TIL...

[nlw0]

I have noticed this seems to be used for Integers as well, where the order of operations should not affect accuracy. sum(::Vector{Int64}) appears to use this method, giving it a penalty over a simple for-loop for smaller arrays.

sum over Int64 values is 50% to 150% slower than a for-loop up to 4k items. I understand the desire to offer the accurate sum for large float arrays, but I find this surprising. Would you consider this a worthy issue, @simeonschaub ?

[nlw0]

I've made an experiment trying to measure the speed and accuracy of sum versus a simple for-loop adding up values from a random vector. I might have made some mistake here, so I'm glad to hear any comments about the code: https://gist.github.com/nlw0/80eca5f6b873e9b75c760b8635c0a593

The main thing I'm not sure is how I'm supposed to measure the accuracy. I've chosen the median of absolute relative errors, but perhaps other measurements are more adequate. Also the distribution of values can probably influence, as well as the data type (I've used Float32), not to mention the machine etc.

And another caveat is that vectorization might have a positive effect in accuracy that should be taken into account. We are not necessarily doing a linear sum here, but a "for-loop with compiler optimizations on".

Here are the results:

Tentative conclusions:

1. The accuracy of the pairwise approach used in sum (i.e. Base.mapreduce_impl) does appear to be limited, or growing very slowly, while the simple linear sum shows a steadily growing inaccuracy.
2. The performance cost does not seem to be trivial until we reach pretty large vector lengths, ~100k values.
3. Performance cost is as high a 5x (!) for smaller arrays.
4. Simple summation resulted in better accuracy for smaller arrays. (insert here all those caveats)

I would suppose this experiment, if there's no major mistake, should incentivize at least trying to tune this block size parameter, or something. It's bigger than that, though because this function appears to work as the "de facto" mapreduce used in Julia, eg #43501 . As a result, anyone using sum for integers is paying a performance price with no accuracy benefit. And that seems to be the case even for small float arrays. I understand the original idea, but I would like to suggest taking a second look at when this function should be employed.

[N5N3]

IIUC, PWSUM has better accuracy when input's length exceeds 2048.
This might be machine dependent, but does show that the default pairwise_blocksize could be tuned more carefully.

As for the high "Performance cost" for small array, it should not be caused by the pairwise reduction, but the initialization.

julia> @btime simpsum($a); @btime pwsum($a, $1024);
  18.337 ns (0 allocations: 0 bytes)
  22.267 ns (0 allocations: 0 bytes)

julia> a = randn(Float32,64);

julia> @btime simpsum($a); @btime pwsum($a, $1024);
  5.400 ns (0 allocations: 0 bytes)
  25.978 ns (0 allocations: 0 bytes)

julia> a = randn(Float32,66);

julia> @btime simpsum($a); @btime pwsum($a, $1024);
  4.500 ns (0 allocations: 0 bytes)
  6.800 ns (0 allocations: 0 bytes)

LLVM doesn't handle unvectorized remainders well, this stongly affects the performance for inputs with small bitsize.

As an easy improvement, we can turn off pw-reduction if we know pw-reduction and simp-reduction are exactly the same.
e.g. +/* for HWInt, & / | / min / max for HWReal.
(It should be safe as the inference accuracy won't inference the result)

[nalimilan]

@nlw0 Do you have evidence showing that integer sum could be significantly faster if we switched to a simple loop instead of using pairwise summation? I imagine this kind of change could be accepted (contrary to dropping pairwise summation for floating point).

[nlw0]

#45000 (comment)

@nalimilan Do you have any power to e.g. revert this closed issue back to open? It's a pretty simple experiment, but I'm not going to spend more time with this issue unless someone with power to do something concrete is listening. I've provided enough evidence and arguments, and it's easy to see from the theory why this piecewise reduction is bound to introduce obstacles to compiler optimizations. As much as I like discussing these things, I don't really see that I should spend more time with this issue. I have provided evidence and arguments, you could try disproving them already. As for not doing piecewise summation for floats, it's in the end a design decision, I happen to think it's a bad one, trying to save users from themselves. In the end it's up to the core designers. By the way, I have shown better precision and time for the non-piecewise summation, so even if we believe it's important to save users from themselves, piecewise summation is not doing it for some input lengths. A case of "victim of your own hubris" in my humble opinion. If anybody really does care either about floating point summation accuracy and/or running time, this is like a sore thumb.

[StefanKarpinski]

I think there's a number of reasons why this issue hasn't gotten any traction. The title suggests that you want to always do linear traversal in mapreduce (that's the title, after all), which sounds like a guarantee, whereas having the option to do other reduction orders is pretty crucial for both speed and accuracy of floating-point reduction, and other cases. Even in the integer case, you really don't want to guarantee strict left-to-right reduction since that precludes SIMD reordering which is very important for performance. So let me assume that the actual intention of this issue is much milder: to point out that there are some specific cases where it might be better to use a different algorithm that we currently are using for mapreduce.

Another reason that I think you're having trouble getting traction is that there are too many different things being discussed here and some of them are convincing and some of them aren't. One fairly simple thing that seems to be very plausible is that it would be faster specifically for integers to use a straightforward SIMD reduction with no pairwise accumulation logic. Since integer + and * are associative, this has no effect on accuracy and is strictly about performance. It's an easy argument to make: change the implementation, show that it's faster. Should be an easy PR to get merged if someone wants to do it.

The floating-point accuracy thing is going to be a harder case to make. What I do know is that simple left-to-right accumulation (as the title and experiments suggest you're advocating for) is the worst possible algorithm in terms of both speed and accuracy, because it has the most skewed tree of operations. The more balanced the tree, the better the error behavior and that tends to be better for speed as well. Whether the current implementation is optimal or not is an entirely different matter, but I think much of the rejection of the issue comes from people (imo correctly) disbelieving that left-to-right reduction would be better.

So, in summary:

• It's definitely possible to make improvements to mapreduce, but a guarantee of strict left-to-right reduction is not going to happen—we need to keep the latitude to reduce in a better order.
• Neither integers nor floats will be faster with a strict left-to-right reduction which the title suggest—you need to at least take advantage of SIMD, which isn't going to be strictly left-to-right. It seems like what you're actually advocating for is a "simple" left-to-right modulo SIMD (I'll call this LTRS) reduction, but one really needs to get into the weeds of the issue to understand that.
• The simplest case to be made kind of gets lost in the noise, which is that integers should use LTRS. I believe this entirely and think it was probably a mistake that they got included in the pairwise reduction logic—pairwise reduction is totally unnecessary overhead for associative numerical types like integers.
• I also fully believe that lift-to-right modulo SIMD reduction for floats is also somewhat faster, but I'm skeptical of claims that it's more accurate; indeed, if I'm reading your experiments right, they show linear error growth for left-to-right and asymptotically constant error for pairwise reduction, which is what is theoretically expected.
• I'm very curious about the LTRS error being lower for arrays smaller than 2^11; I can't understand why that would be the case. I wouldn't be surprised by the error being the same, but why would pairwise have worse error for smaller arrays? If we can understand why that's the case, it certainly suggests that the base case size for pairwise summation should be around 2^11 elements.

A suggestion for improving the summation comparison: use the Xsum package—it computes the exact sum of a collection of floats efficiently and more precisely, even, than using BigFloats.