bugGNU Octave - Bugs: bug #63848, Mean: some single class inputs can...

 
 

bug #63848: Mean: some single class inputs can hit precision limit with sum

Submitter:  Clément Bret <clem_bret>
Submitted:  Tue 28 Feb 2023 02:49:55 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  nrjank
Originator Name:  Open/Closed:  * Closed
Release:  * 7.2.0 Operating System:  * Any
Fixed Release:  8.1.0 Planned Release:  8.1.0
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Mon 04 Mar 2024 04:00:27 PM UTC, comment #29: 

just to add to what markus said, that test was added as a known failure condition for simple mean calculation with double precision numbers.  Unless we implement an alternative mean or summation calculation to avoid such limits, or move to standard operation in quad precision arithmetic, that test is expected to fail. 

Note in the test report it says
PASSES 79 out of 80 tests (1 known bug)

the presence of the bug number indicates to the test routine that it's a known failure.  if you ran the full test suite it would show up as an XFAIL (expected fail).

the comments below do link to some of the other discussions about implementing block based or other summation routines that could preserve precision for some cases, where most effort has stalled out. if you're able to contribute to those it could be of value to octave.

Nicholas Jankowski <nrjank>
Group Member
Mon 04 Mar 2024 03:47:32 PM UTC, comment #28: 

The expected value for that average is not 35184372088832 like you wrote.

`flintmax("double")` is 2^53 (or 9007199254740992). Take that to Wolfram Alfa (or another "calculator" of your choice that has the required precision) and you'll see that the expected result is actually what the test expects currently:
https://www.wolframalpha.com/input?i=%289007199254740992+%2B+2%5E8+-1%29%2F2%5E8

Markus Mützel <mmuetzel>
Group administrator
Mon 04 Mar 2024 12:28:12 PM UTC, comment #27: 

Host OS: Windows 11 23H2 22631.3155
VM: VM VirtualBox 7.0.14 r161095 (Qt5.15.2)
Guest OS: Ubuntu 23.10x64
Octave: 10.0.0

Test mean fails in double test case.

>> test mean

*** assert <63848> (mean ([flintmax("double"), ones(1, 2^8-1, "double")]), 35184372088833-1/(2^8), eps(35184372088833))

  Location  |  Observed  |  Expected  |  Reason
     ()      35184372088832 35184372088833   Abs err 1 exceeds tol 0.0078125 by 1
PASSES 79 out of 80 tests (1 known bug)

The initial investigation concludes that the test case failed due to the wrong assertion condition. The second and third arguments should be 35184372088832 and eps(35184372088832) instead of 35184372088833-1/(2^8) and eps(35184372088833), respectively. Please verify the case described above and suggested solution.

Krzysztof Kijanowski <kkija>
Mon 06 Mar 2023 04:14:19 PM UTC, comment #26: 

Verified fix on the stable branch.  Marking Bug as fixed and closing report.

Rik <rik5>
Group administrator
Fri 03 Mar 2023 03:06:55 AM UTC, comment #25: 

pushed the same changes to stable as
https://hg.savannah.gnu.org/hgweb/octave/rev/942ca5f2c71c

only difference was the old version relied on the sums preserving singles so had to edit the ending output conversion block.

Nicholas Jankowski <nrjank>
Group Member
Thu 02 Mar 2023 08:02:00 PM UTC, comment #24: 

got it. By empty i meant something like 'git merge -s ours' which should do a merge and then completely ignore any incoming changes . it sort of would work out like Rik's option 2.

the mean in default is now quite different from that in stable, so even with such simple changes we'll see what it thinks.

Nicholas Jankowski <nrjank>
Group Member
Thu 02 Mar 2023 07:54:49 PM UTC, comment #23: 

I would try Markus's suggestion in comment #22.  Make the addition of "double" to the calls to sum() on the stable branch, commit it, then merge stable to default and see if there are any conflicts.  I would hope that the merging is actually painless.  If it is not painless then I would try a different trick.  Do 'hg update -C' on the development branch to forget about the changes you just attempted.  Change directory to scripts/statistics and copy mean.m to meantmp.m.  Then do 'hg merge stable' and just accept the changes whatever they are.  Now copy meantmp.m to mean.m overwriting whatever the conflicts were.  Now checkin with 'hg ci -m "maint: merge stable to default".

Rik <rik5>
Group administrator
Thu 02 Mar 2023 05:11:48 PM UTC, comment #22: 

You could:
- commit your changes to the stable branch,
- merge to the default branch,
- resolve any potential conflicts on the default branch,
- commit the merge,
- push to the remote repository.

(I don't know what an "empty merge" is.)

Markus Mützel <mmuetzel>
Group administrator
Thu 02 Mar 2023 05:00:13 PM UTC, comment #21: 

since mean has deviated quite a bit between stable and default, is there any issue with just pushing a patch to stable to make the same effective change, then doing an 'empty' merge to remove the extra head without actually merging anything back to default?  or should the default changeset be backed out and work from that direction?

Nicholas Jankowski <nrjank>
Group Member
Thu 02 Mar 2023 04:47:45 PM UTC, comment #20: 

@rik, sorry didn't reload before the push and missed your comment. will look at getting this into stable.

Nicholas Jankowski <nrjank>
Group Member
Thu 02 Mar 2023 04:13:28 PM UTC, comment #19: 

yes, turns out it was as simple as turning each call into sum (x, dim, "double") except for the large int handling path. Pushed
https://hg.savannah.gnu.org/hgweb/octave/rev/cfeda68b01ad
along with BISTs to capture the single precision issue and a single xfail test to show a similar limit with doubles. marking as ready for test

Nicholas Jankowski <nrjank>
Group Member
Thu 02 Mar 2023 06:20:51 AM UTC, comment #18: 

Given that this produces the wrong answer, I think you should make the changeset to apply against the stable branch of Octave.  This shouldn't go another year before programmer's see a fixed release.

Rik <rik5>
Group administrator
Wed 01 Mar 2023 05:57:26 PM UTC, comment #17: 

i'll look at patching mean now. for the protection of octave i don't work in anything but the simple compiled code by choice, so someone else is going to have to play with bug #61143

Nicholas Jankowski <nrjank>
Group Member
Wed 01 Mar 2023 05:55:15 PM UTC, comment #16: 

OK w/ me

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Wed 01 Mar 2023 05:35:39 PM UTC, comment #15: 

sorry, that's patch #10314

Nicholas Jankowski <nrjank>
Group Member
Wed 01 Mar 2023 05:35:02 PM UTC, comment #14: 

since blocking moves but doesn't completely eliminate where you hit the precision limit, i'd suggest we keep this one just about mean, implement the 'process as double' approach for singles now that patch 10341 has pushed, and then refocus 61143 on sum routines.

Nicholas Jankowski <nrjank>
Group Member
Wed 01 Mar 2023 05:08:11 PM UTC, comment #13: 

Why not to implement Michael Leitner <mleitner> proposal?

In any case should we close this bug ad dup of https://savannah.gnu.org/bugs/?61143

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Wed 01 Mar 2023 04:58:52 PM UTC, comment #12: 

ok, if not consensus, then the result of inertia :)

maybe something good to add to the short(?) projects wiki list.

in any case, mean can still try to workaround the problem in m-code at least for single. we can leave double as a to-do until we have a native quadruple precision type or something happens in bug #61143 to address it direcetly in sum.

Nicholas Jankowski <nrjank>
Group Member
Wed 01 Mar 2023 04:53:55 PM UTC, comment #11: 

I do not think there was any consensus.

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Wed 01 Mar 2023 04:37:41 PM UTC, comment #10: 

yes, that matches with what I had to do to work around the int64/uint64 precision issues that couldn't just be cast to double for bug #54567

the consensus over on 61143 seemed to be to leave naïve summation as is rather than implement a pairwise or similar algorithm, although it's been a few years since that discussion. there's a linked matlab blog article that indicates they implemented a more accurate block based summation across all data types in 2020 after doing so for single in 2017, although as expected they guard their blocksize determination routine.  I'm wondering if that changes in any way the thoughts about a better summation algorithm in octave (with the notion that octave may produce inaccurate answers for certain summations/means/etc where matlab does not)

Nicholas Jankowski <nrjank>
Group Member
Wed 01 Mar 2023 03:43:54 PM UTC, comment #9: 

It is a saturation at flintmax('single') = 16777216.
With double it is 9007199254740992

There are better summations algorithms.
See lots of discussions and proposals in

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

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Wed 01 Mar 2023 01:01:30 PM UTC, comment #8: 

Hmmm, thanks. Didn’t think hard enough on what was actually am happening. The output being ~half the true mean made me think it was the same as with, e.g, (intmax+intmax)/2 overflow, where the fix had to occur in mean, not sum, because the summation overflow was expected and the mean operation  ordering just needed to be fixed (or type changed where possible). As I read bug #61143 Dmitri linked I see it is the small addend size running into the precision limit of the growing summation. Retitling report.

I assume we can hit the same issues for “double” through judicious choice of inputs?

@dmitri, for the most part octave/Matlab does default to double output by except in a few cases (and there is an outtype option to control this) so trying to match expected output type is a compatibility issue that we currently meet and I wouldn’t want to discard.

Nicholas Jankowski <nrjank>
Group Member
Wed 01 Mar 2023 07:44:50 AM UTC, comment #7: 

Is this really an overflow issue with `sum`? The expected approximate result of the summation (80e6*0.5) is still far from `realmax('single')` (approx. 3.4e38).
I'd guess it is an issue caused by the limited precision of floating point numbers...

Markus Mützel <mmuetzel>
Group administrator
Tue 28 Feb 2023 08:37:15 PM UTC, comment #6: 

ok. I was planning on starting the mean/median/std/var port from statistics 1.5.4 shortly over at patch #10314. after starting to bring that in i will look at adding single overflow fixes without breaking any large int64/uint64 handling i fixed in there from bug #54567

Nicholas Jankowski <nrjank>
Group Member
Tue 28 Feb 2023 08:23:13 PM UTC, comment #5: 

If it was up to me, I would leave the output "double" and let
user truncate/cast to whatever they want.

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Tue 28 Feb 2023 07:31:50 PM UTC, comment #4: 

Matlab seems to be following a strategy where they

1. Record class of input
2. Convert input to double
3. Perform calculation
4. Convert result back to input class

We could do that in Octave as well.  Right now we rely on the underlying functions like sum() which automatically perform operations in the class of the input for floating point values.

In this particular case, it would advantageous not to create a duplicate of the original matrix, but merely of type double, which would use a lot of additional memory.  Instead, as Dmitri suggested, add "double" the function calls to sum() and then convert back to result class at the end of the function with a call to cast().

Rik <rik5>
Group administrator
Tue 28 Feb 2023 05:02:03 PM UTC, comment #3: 

just went through fixing something similar for large ints in the statistics pkg version of mean, which we're planning on porting to octave core now that statistics v1.5.4 has been released. trying to maintain output type compatibility with matlab and the 'outtype' option wasn't exactly trivial, but once that's ported over I can look at whether things with single can be handled similarly to how they avoided overflow with ints.  there we also had to be careful of large ints causing precision issues with double, and that shouldn't be an issue here.

Nicholas Jankowski <nrjank>
Group Member
Tue 28 Feb 2023 04:12:00 PM UTC, comment #2: 

A quick-n-dirty fix is toadd "double" qualifier to all
calls to sum() in the mean.m file:

+verbose+
diff -u /usr/local/share/octave/8.0.91/m/statistics/mean.m  ../scripts/statistics/mean.m
--- /usr/local/share/octave/8.0.91/m/statistics/mean.m 2023-02-27 13:04:16.784440819 -0500
+++ ../scripts/statistics/mean.m 2023-02-28 10:45:49.227600518 -0500
@@ -136,7 +136,7 @@
       n = numel (x(:));
       if (omitnan)
         idx = isnan (x);
-        n -= sum (idx(:));
+        n -= sum (idx(:),'double');
         x(idx) = 0;
       endif
       y = sum (x(:), 1) ./ n;
@@ -147,10 +147,10 @@
       n = size (x, dim);
       if (omitnan)
         idx = isnan (x);
-        n = sum (! idx, dim);
+        n = sum (! idx, dim,'double');
         x(idx) = 0;
       endif
-      y = sum (x, dim) ./ n;
+      y = sum (x, dim,'double') ./ n;
     endif
 
   else
@@ -166,10 +166,10 @@
       n = size (x, dim);
       if (omitnan)
         idx = isnan (x);
-        n = sum (! idx, dim);
+        n = sum (! idx, dim,'double');
         x(idx) = 0;
       endif
-      y = sum (x, dim) ./ n;
+      y = sum (x, dim,'double') ./ n;
 
     else
 
@@ -186,10 +186,10 @@
           n = numel (x(:));
           if (omitnan)
             idx = isnan (x);
-            n -= sum (idx(:));
+            n -= sum (idx(:),'double');
             x(idx) = 0;
           endif
-          y = sum (x(:), 1) ./ n;
+          y = sum (x(:), 1,'double') ./ n;
 
         ## for 1 dimension left, return column vector
         case 1
@@ -200,7 +200,7 @@
             if (omitnan)
               x_vec = x_vec(! isnan (x_vec));
             endif
-            y(i) = sum (x_vec, 1) ./ numel (x_vec);
+            y(i) = sum (x_vec, 1,'double') ./ numel (x_vec);
           endfor
           y = ipermute (y, [misdim, dim]);
 
-verbose-

Be aware that this will make mean() to output "double" even for "single" input.

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Tue 28 Feb 2023 03:26:07 PM UTC, comment #1: 

This looks like a duplicate of
https://savannah.gnu.org/bugs/?61143

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Tue 28 Feb 2023 02:49:55 PM UTC, original submission:  

I expect the following two lines to give both 0.5:


>> mean (rand (1,80e6, "single"))
ans = 0.2097
>> mean (rand (1,80e6))
ans = 0.5000


The two results seem to differ once vectors are larger than about 20e6 elements.

Clément Bret <clem_bret>

 

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

Attach Files:
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by kkija (Posted a comment)
  • -email is unavailable- added by mmuetzel (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by nrjank (Posted a comment)
  • -email is unavailable- added by dasergatskov (Posted a comment)
  • -email is unavailable- added by clem_bret (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 12 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2023-03-06 rik5 Open/ClosedOpen Closed
    2023-03-06 rik5 StatusReady For Test Fixed
    2023-03-03 nrjank Fixed ReleaseNone 8.1.0
        Planned Release9.1.0 8.1.0
    2023-03-02 nrjank StatusIn Progress Ready For Test
        Planned ReleaseNone 9.1.0
    2023-03-01 nrjank StatusConfirmed In Progress
        Assigned toNone nrjank
    2023-03-01 nrjank SummaryMean: single class inputs can overflow with sum producing wrong mean Mean: some single class inputs can hit precision limit with sum
    2023-02-28 nrjank StatusNone Confirmed
        Operating SystemGNU/Linux Any
        SummaryMean: wrong result for large single vectors Mean: single class inputs can overflow with sum producing wrong mean

    Back to the top

    Powered by Savane 3.13-02a9.
    Corresponding source code