[Yt-dev] Fwd: [Numpy-discussion] Summation of large float32/float64 arrays

[Matthew Turk]

Just as a quick note, as I suspect at some point the future this will
be useful to someone else, as well.

Forwarded conversation
Subject: Summation of large float32/float64 arrays
------------------------

From: Matthew Turk <matthewturk at gmail.com>
Date: Fri, May 21, 2010 at 1:13 PM
To: Discussion of Numerical Python <numpy-discussion at scipy.org>


Hi all,

I have a possibly naive question.  I don't really understand this
particular set of output:

In [1]: import numpy

In [2]: a1 = numpy.random.random((512,512,512)).astype("float32")

In [3]: a1.sum(axis=0).sum(axis=0).sum(axis=0)
Out[3]: 67110312.0

In [4]: a1.sum()
Out[4]: 16777216.0

I recognize that the intermediate sums may accumulate error
differently than a single call to .sum(), but I guess my concern is
that it's accumulating a lot faster than I anticipated.  (Interesting
to note that a1.sum() returns 0.5*512^3, down to the decimal; is it
summing up the mean, which should be ~0.5?)  However, with a 256^3
array:

In [1]: import numpy

In [2]: a1 = numpy.random.random((256,256,256)).astype("float32")

In [3]: a1.sum(axis=0).sum(axis=0).sum(axis=0)
Out[3]: 8389703.0

In [4]: a1.sum()
Out[4]: 8389245.0

The errors are much more reasonable.  Is there an overflow or
something that occurs with the 512^3?  These problems all go
completely away with a float64 array, but the issue originally showed
up when trying to normalize an on-disk float32 array of size 512^3,
where the normalization factor was off by a substantial factor (>2x)
depending on the mechanism used to sum.  My suspicion is that perhaps
I have a naive misconception about intermediate steps in summations,
or there is a subtlety I'm missing here.

I placed a sample script I used to test this here:

http://pastebin.com/dGbHwFPK

Thanks for any insight anybody can provide,

Matt

----------
From: Robert Kern <robert.kern at gmail.com>
Date: Fri, May 21, 2010 at 1:26 PM
To: Discussion of Numerical Python <numpy-discussion at scipy.org>


It's not quite an overflow.

In [1]: from numpy import *

In [2]: x = float32(16777216.0)

In [3]: x + float32(0.9)
Out[3]: 16777216.0

You are accumulating your result in a float32. With the a.sum()
approach, you eventually hit a level where the next number to add is
always less than the relative epsilon of float32 precision. So the
result doesn't change. And will never change again as long as you only
add one number at a time. Summing along the other axes creates smaller
intermediate sums such that you are usually adding together numbers
roughly in the same regime as each other, so you don't lose as much
precision.

Use a.sum(dtype=np.float64) to use a float64 accumulator.

--
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth."
 -- Umberto Eco
_______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion at scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

----------
From: Matthew Turk <matthewturk at gmail.com>
Date: Fri, May 21, 2010 at 1:32 PM
To: Discussion of Numerical Python <numpy-discussion at scipy.org>


Hi Robert,
Thank you very much for that explanation; that completely makes sense.
I didn't know about the dtype for accumulators/operators -- that did
just the trick.

Much obliged,

Matt