rolling().sum() is numerically unstable

[ilemhadri]

What happened?

On an input array like

array([0.        , 0.        , 0.        , 0.57392103, 0.57392103,
       0.57392103, 0.57392103, 0.57392103, 0.57392103, 0.57392103,
       0.57392103, 0.57392103, 0.        , 0.57392103, 0.57392103,
       0.57392103, 0.57392103, 0.57392103, 0.57392103, 2.29551022,
       2.29551022, 2.29551022, 2.29551022, 2.29551022, 2.29551022,
       2.29551022, 2.29551022, 2.29551022, 2.29551022, 2.29551022,
       2.29551022, 2.29551022, 2.29551022, 0.57383408, 0.57383408,
       0.57383408, 0.57383408, 0.57383408, 0.57383408, 0.        ,
       0.        , 0.        , 0.        , 0.        , 0.        ,
       0.        , 0.        , 0.        , 0.        , 0.        ,
       0.        , 0.        , 0.        , 0.        , 0.        ,
       0.        , 0.        , 0.        , 0.        , 0.        ,
       0.        , 0.        , 0.        , 0.        , 0.        ,
       0.        , 0.        , 0.        , 0.        , 0.        ,
       0.        , 0.        , 0.        , 0.        , 0.        ,
       0.        , 0.        , 0.        , 0.        , 0.        ,
       0.        , 0.        , 0.        , 0.        , 0.        ,
       0.        , 0.        , 0.        , 0.        , 0.        ,
       0.        , 0.        , 0.        , 0.        , 0.        ,
       0.        , 0.        , 0.        , 0.        , 0.        ])

that has all positive values, and then zeros, computing the rolling sum (or mean) leads to negative values

mydata.rolling(t=3, min_periods=1).sum().values[0,0,:]
array([-3.33066907e-14, -3.33066907e-14, -3.33066907e-14,  5.73921029e-01,                                                                                                                                                                                                                     
        2.14784206e+00,  1.72176309e+00,  1.72176309e+00,  1.72176309e+00,                                                                                                                                                                                                                     
        1.72176309e+00,  1.72176309e+00,  1.72176309e+00,  1.72176309e+00,                                                                                                                                                                                                                     
        1.14784206e+00,  1.14784206e+00,  1.14784206e+00,  1.72176309e+00,                                                                                                                                                                                                                     
        1.72176309e+00,  1.72176309e+00,  1.72176309e+00,  3.44335228e+00,
        5.16494146e+00,  6.88653065e+00,  6.88653065e+00,  6.88653065e+00,
        6.88653065e+00,  6.88653065e+00,  6.88653065e+00,  6.88653065e+00,
        6.88653065e+00,  6.88653065e+00,  6.88653065e+00,  6.88653065e+00,
        6.88653065e+00,  5.16485452e+00,  3.44317838e+00,  1.72150224e+00,
        1.72150224e+00,  1.72150224e+00,  1.72150224e+00,  1.14766816e+00,
        5.73834081e-01, -3.35287353e-14, -3.35287353e-14, -3.35287353e-14,
       -3.35287353e-14, -3.35287353e-14, -3.35287353e-14, -3.35287353e-14,
       -3.35287353e-14, -3.35287353e-14, -3.35287353e-14, -3.35287353e-14,
       -3.35287353e-14, -3.35287353e-14, -3.35287353e-14, -3.35287353e-14,
       -3.35287353e-14, -3.35287353e-14, -3.35287353e-14, -3.35287353e-14,
       -3.35287353e-14, -3.35287353e-14, -3.35287353e-14, -3.35287353e-14,
       -3.35287353e-14, -3.35287353e-14, -3.35287353e-14, -3.35287353e-14,
       -3.35287353e-14, -3.35287353e-14, -3.35287353e-14, -3.35287353e-14,
       -3.35287353e-14, -3.35287353e-14, -3.35287353e-14, -3.35287353e-14,
       -3.35287353e-14, -3.35287353e-14, -3.35287353e-14, -3.35287353e-14,
       -3.35287353e-14, -3.35287353e-14, -3.35287353e-14, -3.35287353e-14,
       -3.35287353e-14, -3.35287353e-14, -3.35287353e-14, -3.35287353e-14,
       -3.35287353e-14, -3.35287353e-14, -3.35287353e-14, -3.35287353e-14,
       -3.35287353e-14, -3.35287353e-14, -3.35287353e-14, -3.35287353e-14,
       -3.35287353e-14, -3.35287353e-14, -3.35287353e-14, -3.35287353e-14])

Both arrays have dtype = float64.

The issue aggravates as the rolling window increases.

What did you expect to happen?

the rolling calculation could be more numerically precise by keeping track for instance of the Kahan compensation term.
https://en.wikipedia.org/wiki/Kahan_summation_algorithm

MVCE confirmation

• Minimal example — the example is as focused as reasonably possible to demonstrate the underlying issue in xarray.
• Complete example — the example is self-contained, including all data and the text of any traceback.
• Verifiable example — the example copy & pastes into an IPython prompt or Binder notebook, returning the result.
• [X ] New issue — a search of GitHub Issues suggests this is not a duplicate.

Environment

this is reproducible across xarray versions, but mine is 2022.09.0.

[mathause]

Do you have bottleneck installed? Could you try xr.set_options(use_bottleneck=False)?

See pydata/bottleneck#379

[ilemhadri]

Do you have bottleneck installed? Could you try xr.set_options(use_bottleneck=False)?

use_bottleneck=False seems to solve it on the latest version of xarray, in the sense that the negative values are now exactly 0. However, I don't think it's compatible with previous versions of xarray like 0.11.3:

ValueError: argument name 'use_bottleneck' is not in the set of valid options set(['keep_attrs', '
enable_cftimeindex', 'cmap_sequential', 'arithmetic_join', 'warn_for_unclosed_files', 'file_cache_maxsize', 'cmap_
divergent', 'display_width'])                             

In our setting, we would love to preserve backwards-compatibility and be able to do this in v0.11.3 as well!

[mathause]

You'll have to uninstall bottleneck then - the use_bottleneck option was added more recently.

[ilemhadri]

You'll have to uninstall bottleneck then - the use_bottleneck option was added more recently.

@mathause Does xarray use bottleneck in other places than rolling operations?

That'd help me assess the scope of this change!