bottleneck : Wrong mean for float32 array

[matteodefelice]

I think it is better to have this discussion here instead of on the dask page dask/dask#2095

This is the replicable "bug":

ds = xarray.open_dataset('/opt/data/ERAIN/ERAIN-t2m-1983-2012.seasmean.nc')
ds.var167.mean()
Out[14]: 
<xarray.DataArray 'var167' ()>
array(261.6441345214844)
ds.var167.data.mean()
Out[15]: 278.62466

The dataset is ~65 MB, here the file https://www.dropbox.com/s/xtj3fm7ihtbwd5r/ERAIN-t2m-1983-2012.seasmean.nc?dl=0
It is a quite normal NetCDF (no NaN), just processed with CDO as you can see on the dask issue.

[fmaussion]

I can't reproduce this:

In [6]: ds = xr.open_dataset('./Downloads/ERAIN-t2m-1983-2012.seasmean.nc')

In [7]: ds.var167.mean()
Out[7]: 
<xarray.DataArray 'var167' ()>
array(278.6246643066406, dtype=float32)

In [8]: ds.var167.data.mean()
Out[8]: 278.62466

which version of xarray, dask, python are you using?

[matteodefelice]

Ok, I am on MacOS:

• Python 2.7.13 from Macports
• Dask 0.14.1 from Macports
• xarray from GitHub

[andrew-c-ross]

Also on macOS, and I can reproduce.

Using python 2.7.11, xarray 0.9.1, dask 0.14.1 installed through Anaconda. I get the same results with xarray 0.9.1-38-gc0178b7 from GitHub.

In [3]: ds = xarray.open_dataset('ERAIN-t2m-1983-2012.seasmean.nc')

In [4]: ds.var167.mean()
Out[4]:
<xarray.DataArray 'var167' ()>
array(261.6441345214844, dtype=float32)

Curiously, I get the right results with skipna=False...

In [10]: ds.var167.mean(skipna=False)
Out[10]:
<xarray.DataArray 'var167' ()>
array(278.6246643066406, dtype=float32)

... or by specifying coordinates to average over:

In [5]: ds.var167.mean(('time', 'lat', 'lon'))
Out[5]:
<xarray.DataArray 'var167' ()>
array(278.6246643066406, dtype=float32)

[fmaussion]

Does it make a difference if you load the data first? (ds.var167.load().mean()) Or use python 3?

[andrew-c-ross]

I think this might be a problem with bottleneck? My interpretation of _create_nan_agg_method in xarray/core/ops.py is that it may use bottleneck to get the mean unless you pass skipna=False or specify multiple axes. And,

In [2]: import bottleneck
In [3]: bottleneck.__version__
Out[3]: '1.2.0'

In [6]: bottleneck.nanmean(ds.var167.data)
Out[6]: 261.6441345214844

Forgive me if I'm wrong, I'm still a bit new.

[shoyer]

Yes, this is probably related to the fact that .mean() in xarray uses bottleneck if available, and bottleneck has a slightly different mean implementation, quite possibly with a less numerically stable algorithm.

The fact that the dtype is float32 is a sign that this is probably a numerical precision issue. Try casting with .astype(np.float64) and see if the problem goes away.

If you really cared about performance using float32, the other thing to do to improve conditioning is to subtract and add a number close to the mean, e.g., (ds.var167 - 270).mean() + 270.

[matteodefelice]

Thanks all guys for the replies.
@Aegaeon I get the same your results with bottleneck...
@shoyer The point is that I haven't decided the use of float32 and — yes — using .astype(np.float64) solves the issue...the point is that is not an expected behaviour, with such standard dataset I would not expect any problem related to numerical precision...

[shoyer]

@matteodefelice you didn't decide on float32, but your data is stored that way. It's really hard to make choices about numerical precision for computations automatically: if we converted automatically to float64, somebody else would be complaining about unexpected memory usage :).

Looking at our options, we could:

1. Stop using bottleneck on float32 arrays, or provide a flag or option to disable using bottleneck. This is not ideal, because bottleneck is much faster.
2. Automatically convert float32 arrays to float64 before doing aggregations. This is not ideal, because it could significant increase memory requirements.
3. Add a dtype option for aggregations (like NumPy) and consider defaulting to dype=np.float64 when doing aggregations on float32 arrays. I would generally be happy with this, but bottleneck currently doesn't provide the option currently.
4. Write a higher precision algorithm for bottleneck's mean.

[leifdenby]

Sorry to unearth this issue again, but I just got bitten by this quite badly. I'm looking at absolute temperature perturbations and bottleneck's implementation together with my data being loaded as float32 (correctly, as it's stored like that) causes an error on the size of the perturbations I'm looking for.

Example:

In [1]: import numpy as np
   ...: import bottleneck

In [2]: a = 300*np.ones((800**2,), dtype=np.float32)

In [3]: np.mean(a)
Out[3]: 300.0

In [4]: bottleneck.nanmean(a)
Out[4]: 302.6018981933594

Would it be worth adding a warning (until the right solution is found) if someone is doing .mean() on a DataArray which is float32?

Based a little experimentation (https://gist.github.com/leifdenby/8e874d3440a1ac96f96465a418f158ab) bottleneck's mean function builds up significant errors even with moderately sized arrays if they are float32, so I'm going to stop using .mean() as-is from now and always pass in dtype=np.float64.

[shoyer]

Would it be worth adding a warning (until the right solution is found) if someone is doing .mean() on a DataArray which is float32?

I would rather pick option (1) above, that is, "Stop using bottleneck on float32 arrays"

[aquasync]

Is it worth changing bottleneck to use double for single precision reductions? AFAICT this is a matter of changing npy_DTYPE0 to double in the float{64,32} versions of functions in reduce_template.c.