optimizationapldyalog

Dyalog APL does sum `+/` from the wrong end, violating the specification


Of course, many functions f are not associative in the sense that a f(b f c) (which is often just written a f b f c in APL) is not always the same as (a f b)f c. This is also the case when f is addition + of floating-point numbers.

Fortunately, in APL there is a simple convention to evaluate from right to left. As an example, ÷/2 1 2 gives 4 (two divided by one half), not 1 (two wholes divided by two).

But to my horror, some optimization (they tell me "SIMD") in the implementation of the Dyalog interpreter violates this!

As a first example, try:

(a b c)←1E¯6 ¯1E13 1E13

a+b+c    ⍝ gives 0.000001

+/a b c  ⍝ gives 0 but should be equivalent to the previous one

The +/ adds up in a bad order.

As a second example:

x←?10000000⍴0     ⍝ 10 million random floats

(+/x)-({⍺+⍵}/x)   ⍝ the two expressions "+/x" and "{⍺+⍵}/x" should be the same
                  ⍝ but this difference is nonzero with a very high probability

Again, +/ sums in a bad order.

Try it online!

In every case I have tried where +/x was wrong, it produced exactly the same value as {⍺+⍵}/⊖x, i.e. what you would get if you summed the vector from the wrong end.

My question here is: Is there any reason why the optimization used behind the scene could not take the numbers from the correct end? Is this bug already reported/documented?


Addition: I am realizing there are lots of more "natural" examples. You can try:

(a b c)←0.6 0.7 0.4

⎕PP←17   ⍝ be sure we can detect different numbers

a+b+c    ⍝ gives 1.7000000000000002

+/a b c  ⍝ gives 1.6999999999999997 but should be equivalent to the previous one

Further addition: Product ×/ is also affected; for example there is a nonzero difference between 0.6×0.7×0.4 and ×/0.6 0.7 0.4.


Least-Common-Multiple Reduce ∧/ does not seem to be affected. (Note that least common multiple is so implemented that already for two numbers, the order matters, for example 0.5∧0.7 is different from 0.7∧0.5; but operator reduce / does what it should on as far as I can see.)


Solution

  • As Silas writes, you are not the first person to notice this, I am sure there are people who "discover" this at regular intervals. I am going to ask the documentation team to add some words about it in the next revision of the documentation.

    Dyalog APL has done +/ from left to right since the beginning of time (1980's). The original implementors are no longer with us, but I'll wager it was done because sequential memory access performed better, because they were probably implementing +\ at the same time, and perhaps had a pragmatic thought along the lines of "anyone who has actual floating-point data that depends on the order that the sums are done in is living dangerously anyway".

    Of course, there are people who are forced to live dangerously, using 64-bit binary floats to represent large financial quantities (which is why we implemented support for 128-bit decimal floats in 2011). People living on the edge devise strategies to do things like sorting the data and rounding regularly to improve accuracy, and from that point they don't care about the order that Dyalog APL computes the sums, so long as every release of Dyalog APL does it the same way as the previous one.

    We did once put out a pre-release where we had optimized some summations and also did them according to the spec, but were immediately smacked down by a customer who said this new APL system failed their internal QA. They wouldn't even like us to produce "better" results; once you are producing financial reports you want the same numbers every year.

    The issue regularly pops up in planning meetings, but so far it hasn't made it to the top of the list. For the reason explained above, we cannot change the default behaviour of +/, and having some kind of global switch might complicate building applications from components that might need different behaviours. It has been suggested that 1⊥, which is equivalent to +/ or +⌿ on some arrays, can be allowed to use SIMD or other optimisations – but I find it distasteful to use this "trick" as the official way to get a different implementation. The obvious alternative would be to use the variant operator, so you could define something along the lines of:

     fastsum←+/⍠'Method' 'Fast'
    

    ... to allow SIMD and other optimisations that might not be stable. Or

     isosum←+/⍠'Method' 'ISO'
    

    ... if you want the right-to-left iso standard implementation.

    We welcome ideas on how to go forward with this, it's going to happen at some point!

    Disclaimer: I am the CTO of Dyalog.