[Octave-bug-tracker] [bug #61143] Functions sum and mean returns wrong answer for single precision input

[Kai Torben Ohlhus]

Update of bug #61143 (project octave):

                  Status:                    None => Need Info              

    _______________________________________________________

Follow-up Comment #20:

Interestingly, the topic of having a more precise summation pops up every few
years (last time in 2019 for bug #56884).  However, due to the rare
occurrences in practice with users aware of numerical computations it gets
lost in the sand after a hot discussion 😇

Personally, I tend to agree to comment #18.  Leave things as they are, improve
and better document the algorithm used for the "extra" option.  And close this
either as duplicate of bug #56884 or as "Won't fix".

https://octave.org/doc/v6.3.0/XREFsum.html

Exceptional (large / ill-conditioned input) cases require special compensated
summation algorithms (which there exist a lot in research).  To my knowledge a
"one size fits all" numerical algorithm is not available today 😇

Like jwe said in comment #7, mostly the bar to when things "break" is only
pushed a little higher with higher computational costs.

For example Matlab R2021a seems to accumulate single precision in double
precision (nobody knows what Matlab does 🤫) as the OPs example works. 
However, just changing to another pathological example


>> x=single(0.1*ones(1e+8,1));
>> sum(x)

ans =

  single

     9999961


and you can argue you paid hundreds of dollars for a wrong result (even
without getting a simple warning) or you got what you can get from numerical
computations. The interpretation depends on your point of view.  Is it a bug
or is it a feature?

Regarding comment #19: if nothing is said about the summation algorithm, like
in Octave (perhaps Matlab too), most likely good-old recursive (or naive) sum
is used.  With this default, I know what I get and can expect.  This is speed
by processor and compiler internals and IEEE-754 arithmetic which has well
studied and known limitations, a user aware of numerical computations can deal
with 🙂

However, a "heuristic black box sum" (perhaps Matlab) that covers some
exceptional rare input cases "good enough" (see example above), still fails
other pathological (mostly constructed) examples (as above in Matlab), by
users who don't know or care about their input data.

A users who needs 100% summation precision must look at symbolic computation
tools.  Numerical tools like Octave/Matlab offer speed for known precision
limitations.

Finally, some pointers to the relevant code in Octave:

https://hg.savannah.gnu.org/hgweb/octave/file/aedebbc6b765/libinterp/corefcn/data.cc#l2904

https://hg.savannah.gnu.org/hgweb/octave/annotate/aedebbc6b765/liboctave/operators/mx-inlines.cc#l1688

The currently 2009 (https://hg.savannah.gnu.org/hgweb/octave/rev/192d94cff6c1)
implemented sum algorithms of "extra" for double is

Ogita, Rump, Oishi: "Accurate Sum and Dot Product"
SIAM J. Sci. Comput., 26(6), 1955–1988, 2005.
https://doi.org/10.1137/030601818

If you make changes to the default (recursive) summation algorithm: please
study this algorithm pretty well, perform lots of tests and document for
future programmers any assumption (for optimization, correctness, etc.) you
made.


    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?61143>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/