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

[Michael Leitner]

Follow-up Comment #30, bug #61143 (project octave):

I have attached an implementation of pair-wise summation. In order to allow
simple testing, I did not write it as a patch against the internal sum, but as
an external compiled oct-file. You just have to call

++
mkoctfile sum_pair3.cc -o sum_pair3.oct
--

to use it as sum_pair3(). Again for simplicity, it assumes the input to be in
class double, it doesn't do any input validation, and it sums over all
elements (has no notion of multi-dimensional arrays).

Some remarks: on the architecture where I implemented that (a few years old
Intel notebook, running 32-bit octave on 32-bit debian), the internal sum
seems to use 80-bit extended precision internally. You can test that for your
own architecture by doing

++
delta=1e-5;sum([1 delta -1])-delta
--

for different delta. The returned discrepancy is indicative of the used
precision -- on my computer it is always on the order of about 5e-20, where
the precision of the 63-bit mantissa of 80-bit extended precision would lead
to 2^-63=1.1e-19. 

So internal sum has a bit more than three decimal digits headroom to guard
against cancellation. For random input values, the accumulated error will show
Brownian motion, thus scale proportional to eps*sqrt(N). This means that you
have to go to about N=1e7 until you will see that. But again, this is only for
non-systematic inputs, while significant deficiencies can be exposed for
special inputs, see below.

My implementation of pair-wise summation conceptually corresponds to padding
the input to a power of two, and then having

++
function Sum=pairwise_sum(input)
if length(input)>1
  Sum=pairwise_sum(input(1:end/2))+pairwise_sum(input(end/2+1:end));
else
  Sum=input;
end
--

Of course it is implemented more efficiently. It emits also the time taken for
the summation, and this is because seemingly calling oct-files has an overhead
proportional to the length of input -- I suspect this is because the input is
copied. Thus, a plain tic/toc will give you longer timings than what is really
due to the summing (and if this should eventually replace the internal sum,
there will be no copying). In any case, contrary to expectation, the time
spent on summing in my implementation is less than for the naive sum, about
60%. And just as for the internal sum, I used extended precision for the
running totals, which does not have any effect on the timings compared to
using doubles.

And to the accuracy: summing always the same number (with low-order bits set)
is most critical for naive summing, see

++
N=10000000;b=pi*ones(N,1);sum(b)-pi*N
--

Characteristically, the error goes with the square of N in this case and
easily reaches 1000 eps. On the other hand, sum_pair3() is exact. 

And this is not only for such a constructed example: For instance, summing
sin() over a full period evaluated at equi-distant points should sum exactly
to zero. 

++
N=1000000;x=(0.5:N)/N*2*pi;d=sin(x);[sum(d) sum_pair3(d)]
--

shows here an error of 1.9e-12. This could be thought of to be due to
inaccuracies in sin() itself, but the pair-wise sum indeed gives zero, so this
is the error of the naive summation.

Concludingly: this shall show that pair-wise summation can practically be
implemented in compiled code, ist as efficient (actually, in my implementation
more efficient) than the naive sum as it is implemented now, and solves
accuracy issues. This is for summing double inputs, while according to my
point of view, summing single-precision inputs (which this bug report was
initially about) should be done by the identical algorithm (also in extended
precision) -- there is no argument against doing that, as far as I can see.

(file #51950)
    _______________________________________________________

Additional Item Attachment:

File name: sum_pair3.cc                   Size:1 KB
    <https://file.savannah.gnu.org/file/sum_pair3.cc?file_id=51950>



    _______________________________________________________

Reply to this item at:

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

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