# bug #56884: mean can overflow

 Submitted by: Marco Caliari Submitted on: Thu 12 Sep 2019 09:07:22 AM UTC Category: Octave Function Severity: 3 - Normal Priority: 5 - Normal Item Group: Inaccurate Result Status: Patch Submitted Assigned to: None Originator Name: Open/Closed: Open Release: dev Operating System: GNU/Linux

## Post a Comment

Add a New Comment

You are not logged in

Please log in, so followups can be emailed to you.

## Discussion

Tue 17 Sep 2019 09:25:36 AM UTC, comment #16:

Before replacing the current implementation, the correct handling of all options and input classes has to be taken care of. And I do not know whether this is worth the effort: the best solution would definitely be to do pairwise_mean on the compiled level. After all, my implementation is still much slower than optimal, because with vectorized code you lose cache efficiency, and once this is implemented correctly, there would be no need to touch this ever again.

Michael Leitner <mleitner>
Tue 17 Sep 2019 08:27:33 AM UTC, comment #15:

@Michael: you are right, pairwise_mean is the best code, on average. I think it should replace the actual m-file mean. And probably also the harmonic mean (that is, mean(x,'h') should become 1./pairwise_mean(1./x)). The example il

octave:43> v=1./[1.7e308,1.5e308]
v =

5.882352941176472e-309   6.666666666666668e-309

octave:44> mean(v,'h')
ans =    0.000000000000000e+00
octave:45> 1./pairwise_mean(1./v)
ans =   6.250000000000003e-309

Marco Caliari <caliari>
Tue 17 Sep 2019 07:35:19 AM UTC, comment #14:

Thanks for testing. But please note that when we talk about performance, we have to consider the following: Let's say we have x=randn(N,M), and we want to take the mean over the first dimension. If both N and M are small, the function invocation will take much longer than the actual evaluations, so this is no issue. So the point is whether N or M (or both) are large. For pairwise_mean, the computational effort should scale as N*M*log(N), for iter_mean like N*(M+const), where I have written the time per loop iteration explicitly. At my computer, const is about 2000, that is, only when the M becomes larger than 2000, the cost of the looping starts to become insignificant. Your example is such a case, you have M=201*202. So indeed, in this case your code is faster. However, I would argue that most often (at least for my use cases I can claim that) averages are taken over one-dimensional arrays. And in this case, pairwise_mean is always (much) faster. For instance, with N=100000 and M=1, I have 5 milliseconds for pairwise_mean and 2.7 seconds for iter_mean.

Further, it seems to me that it is also more accurate. As neither mine nor your algorithm as yet has any notion of object classes, a quick test is given by

x=randn(10000,1);
double([iter_mean(single(x)) pairwise_mean(single(x))])-mean(x)

Here pairwise_mean seems to be more accurate by about an order of magnitude on average, which I would expect to even increase for larger vectors.

Michael Leitner <mleitner>
Tue 17 Sep 2019 07:06:06 AM UTC, comment #13:

I see the following

octave:29> x=randn(200,201,202);
octave:30> tic,mean(x);,toc
Elapsed time is 0.00880623 seconds.
octave:31> tic,mean(x);,toc
Elapsed time is 0.00801992 seconds.
octave:32> tic,iter_mean(x);,toc
Elapsed time is 0.055238 seconds.
octave:33> tic,iter_mean(x);,toc
Elapsed time is 0.0553889 seconds.
octave:34> tic,pairwise_mean(x);,toc
Elapsed time is 0.352085 seconds.
octave:35> tic,pairwise_mean(x);,toc
Elapsed time is 0.358778 seconds.
octave:36> ver
----------------------------------------------------------------------
GNU Octave Version: 5.1.1 (hg id: 232102b1e920)
GNU Octave License: GNU General Public License
Operating System: Linux 4.18.0-25-generic #26~18.04.1-Ubuntu SMP

So, mean is the fastest, then iter_mean, finally pairwise_mean. I attach a new version of iter_mean which fixes a problem with simple vectors.

Marco Caliari <caliari>
Tue 17 Sep 2019 06:52:41 AM UTC, comment #12:

No, I did not. I was just using Marco's file as template and did not look at the tests. Indeed there was an error to fix. Now the tests pass the same as Marco's: all pass deterministically but the last, which can fail stochastically. I would say, about one time in five tries, and this holds for both algorithms (but for different vectors x). But this is no problem, this test is about subnormal numbers at the extreme end, where the relative eps becomes of the order of one in any case.

As to Kahan summation: no, I was talking about native sum, not m-code mean. But this was probably premature: What I was testing was code like

a=ones(100,1);a(1)=flintmax();sum(a)-flintmax()

putting the flintmax either to the beginning or the end and varying the vector length.

Indeed, at my office computer (Octave 4.0.3 on amd64) these behave exactly like naive summation (that is, not so bad when flintmax is at the end, but very bad when it is at the beginning), while using sum(a,"extra") is more accurate, which looks like it was using pairwise summation, but definitely not Kahan.

In contrast, at my home computer (older Octave on x86, I think) these lines of code evaluated much more accurate even without "extra", that's why I thought that this was Kahan summation at work. But perhaps it was rather 80-bit extended precision (which would probably the best bet for "extra", being perhaps even faster)? I will have to look again.

I did not look beyond libinterp/corefcn/data.cc, where the array_value ().sum etc. would be defined. But anyway, this bug is about mean, not sum.

Michael Leitner <mleitner>
Tue 17 Sep 2019 12:23:41 AM UTC, comment #11:

Michael Leitner: You did not really intend the tests
for your pairwise_mean to test iter_mean did you?

Also, I did not see anything that looked like Kahan summation
in scripts/statistics/mean.m. Am I missing something?

Michael Godfrey <godfrey>
Mon 16 Sep 2019 08:12:56 PM UTC, comment #10:

Now with attachment.

Michael Leitner <mleitner>
Mon 16 Sep 2019 08:11:35 PM UTC, comment #9:

The problem with Marco's implementation is that it loops over the entries to be averaged, and thus is unsuited for large dimensions as long as there is not JIT. Attached is another proposal based on pairwise summation (implemented in vectorized form), including the idea of the iterative mean algorithm that evades over- and underflows. At first glance, I would say that it should not be worse than iterative mean (I would hope that it is better, as it uses only division by two), and it definitely is faster when written as m-code. I would say the best solution would be to implement my proposal as compiled code, where you can use locality of reference much better. Thus, I would expect it to be even faster than the native sum, which seems to be using Kahan summation. At the cost of increased runtime, this is practically perfectly accurate (that is, the returned value is the nearest floating-point value to the exact solution), however for mean the accuracy is not needed (I do not see the point in the returned value being more accurate on an absolute scale than the inputs).

Michael Leitner <mleitner>
Mon 16 Sep 2019 03:04:51 PM UTC, comment #8:

Marco,

Your code looks good. But, careful testing would be a
good idea.

There are choices about how to handle the various
algorithms for evaluating mean(x).

First, it would be good to decide on what algorithms should
be implemented, and it should be decided which should be
the default. A suggestion would be to provide at
least your iterative and also Kahan's algorithm which is
given in pseudocode at:
https://en.wikipedia.org/wiki/Kahan_summation_algorithm
and the simple code as is.

The choice of the default deserves some thought and testing.
I would be tempted to suggest the Kahan algorithm based on
its robustness (this would tend to help people who have not
thought much about it and just expect "right" answers no matter
what input values they use). It may be that your iter_mean is
faster, but I do not know..

Anyhow, your code (with suitable text for the Manual) should
be installed.

Thanks!

Michael Godfrey <godfrey>
Mon 16 Sep 2019 06:54:53 AM UTC, comment #7:

I attach a simple implementation of the iterative mean. I think I used subsref for the first time in my life. Comments welcome.

Marco Caliari <caliari>
Sat 14 Sep 2019 06:51:54 PM UTC, comment #6:

>> x = rand(10,1) * 1e-323;
>> n = 10; dim = 1;
>> y = sum (x / n, dim)

y = 0

Oh no! Now my edge case is broken :(

Ceral Paquet <octavebugs>
Thu 12 Sep 2019 03:24:32 PM UTC, comment #5:

I think the proposed simple solution can underflow for small values of 'x' or large vector size 'n'. This is not better, anyway.

y = sum (x / n, dim);

Armin MÃ¼ller <arm_in>
Thu 12 Sep 2019 12:23:44 PM UTC, comment #4:
Michael Godfrey <godfrey>
Thu 12 Sep 2019 12:01:53 PM UTC, comment #3:

The iterative mean was proposed sometime in the early 1960's.
It is definitely better than what is currently in Octave.
I have the impression that Kahan is better, but a bit more
review would be good. Checking on how R does it would
make sense.

Michael Godfrey <godfrey>
Thu 12 Sep 2019 11:44:42 AM UTC, comment #2:

I do not think that Kahan summation can avoid the overflow. Maybe the iterative algorithm http://www.heikohoffmann.de/htmlthesis/node134.html

Marco Caliari <caliari>
Thu 12 Sep 2019 11:36:11 AM UTC, comment #1:

There is a large literature on accurate summation.
A place to start is:

Certainly, the current implementation is not appropriate.

Kahan's algorithm is described at:
https://en.wikipedia.org/wiki/Kahan_summation_algorithm

Michael Godfrey <godfrey>
Thu 12 Sep 2019 09:07:22 AM UTC, original submission:

The function mean can overflow

> mean([1.7e308,1.5e308])
ans =  Inf

even if the result is clearly 1.6e308. It appear to me that the simplest solution is to replace

y = sum (x, dim) / n;

with

y = sum (x / n, dim);

around line 153 of mean.m. By the way, matlab overflows too.

Marco Caliari <caliari>

## Attached Files

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:

Comment:

Attached Files
file #47537:  pairwise_mean.m added by mleitner (1KiB - text/x-matlab)
file #47536:  pairwise_mean.m added by mleitner (1KiB - text/x-objcsrc)
file #47524:  iter_mean.m added by caliari (844B - text/x-objcsrc)

## Dependencies

Depends on the following items: None found

Items that depend on this one: None found

## Mail Notification Carbon-Copy List

Carbon-Copy List
• -email is unavailable- added by mleitner (Posted a comment)
• -email is unavailable- added by octavebugs (Posted a comment)
• -email is unavailable- added by arm_in (Posted a comment)
• -email is unavailable- added by godfrey (Posted a comment)
• -email is unavailable- added by caliari (Submitted the item)
•

## Votes

Do you think this task is very important?
If so, you can add your encouragement to it.
This task has 0 encouragements so far.

Only project members can vote.

## History

Follow 4 latest changes.

Date Changed by Updated Field Previous Value => Replaced by
2019-09-17 mleitner Attached File- Added pairwise_mean.m, #47537
2019-09-16 mleitner Attached File- Added pairwise_mean.m, #47536
2019-09-16 caliari Attached File- Added iter_mean.m, #47524
StatusNone Patch Submitted

Powered by Savane 3.5