Rolling window (moving average, moving std, and more)
Hi, Implementing moving average, moving std and other functions working over rolling windows using python for loops are slow. This is a effective stride trick I learned from Keith Goodman's <kwgoodman@gmail.com> Bottleneck code but generalized into arrays of any dimension. This trick allows the loop to be performed in C code and in the future hopefully using multiple cores. import numpy as np def rolling_window(a, window): """ Make an ndarray with a rolling window of the last dimension Parameters ---------- a : array_like Array to add rolling window to window : int Size of rolling window Returns ------- Array that is a view of the original array with a added dimension of size w. Examples -------- >>> x=np.arange(10).reshape((2,5)) >>> rolling_window(x, 3) array([[[0, 1, 2], [1, 2, 3], [2, 3, 4]], [[5, 6, 7], [6, 7, 8], [7, 8, 9]]]) Calculate rolling mean of last dimension: >>> np.mean(rolling_window(x, 3), -1) array([[ 1., 2., 3.], [ 6., 7., 8.]]) """ if window < 1: raise ValueError, "`window` must be at least 1." if window > a.shape[-1]: raise ValueError, "`window` is too long." shape = a.shape[:-1] + (a.shape[-1] - window + 1, window) strides = a.strides + (a.strides[-1],) return np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides) Using np.swapaxes(-1, axis) rolling aggregations over any axis can be computed. I submitted a pull request to add this to the stride_tricks module. Erik
Hi Erik, This is really neat ! Do I understand correctly, that you mean by "stride tricks", that your rolling_window is _not_ allocating any new memory ? IOW, If I have a large array using 500MB of memory, say of float32 of shape 125,1000,1000 and I want the last axis rolling of window size 11, what would the peak memory usage of that operation be ? How about renaming the option `window` to `window_size` (first I was thinking of things like hamming and hanning windows...)... ? Thanks, Sebastian Haase On Sat, Jan 1, 2011 at 5:29 AM, Erik Rigtorp <erik@rigtorp.com> wrote:
Hi,
Implementing moving average, moving std and other functions working over rolling windows using python for loops are slow. This is a effective stride trick I learned from Keith Goodman's <kwgoodman@gmail.com> Bottleneck code but generalized into arrays of any dimension. This trick allows the loop to be performed in C code and in the future hopefully using multiple cores.
import numpy as np
def rolling_window(a, window): """ Make an ndarray with a rolling window of the last dimension
Parameters ---------- a : array_like Array to add rolling window to window : int Size of rolling window
Returns ------- Array that is a view of the original array with a added dimension of size w.
Examples -------- >>> x=np.arange(10).reshape((2,5)) >>> rolling_window(x, 3) array([[[0, 1, 2], [1, 2, 3], [2, 3, 4]], [[5, 6, 7], [6, 7, 8], [7, 8, 9]]])
Calculate rolling mean of last dimension: >>> np.mean(rolling_window(x, 3), -1) array([[ 1., 2., 3.], [ 6., 7., 8.]])
""" if window < 1: raise ValueError, "`window` must be at least 1." if window > a.shape[-1]: raise ValueError, "`window` is too long." shape = a.shape[:-1] + (a.shape[-1] - window + 1, window) strides = a.strides + (a.strides[-1],) return np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides)
Using np.swapaxes(-1, axis) rolling aggregations over any axis can be computed.
I submitted a pull request to add this to the stride_tricks module.
Erik
On Mon, Jan 3, 2011 at 05:13, Sebastian Haase <seb.haase@gmail.com> wrote:
Hi Erik, This is really neat ! Do I understand correctly, that you mean by "stride tricks", that your rolling_window is _not_ allocating any new memory ? Yes, it's only a view.
IOW, If I have a large array using 500MB of memory, say of float32 of shape 125,1000,1000 and I want the last axis rolling of window size 11, what would the peak memory usage of that operation be ? It's only a view of the array, no copying is done. Though some operations like np.std() will copy the array, but that's more of a bug. In general It's hard to imagine any speedup gains by copying a 10GB array.
How about renaming the option `window` to `window_size` (first I was thinking of things like hamming and hanning windows...)... ? Sounds fare.
On Mon, Jan 3, 2011 at 5:37 AM, Erik Rigtorp <erik@rigtorp.com> wrote:
It's only a view of the array, no copying is done. Though some operations like np.std() will copy the array, but that's more of a bug. In general It's hard to imagine any speedup gains by copying a 10GB array.
I don't think that np.std makes a copy of the input data if the input is an array. If the input is, for example, a list, then an array is created.
On Mon, Jan 3, 2011 at 10:36, Keith Goodman <kwgoodman@gmail.com> wrote:
On Mon, Jan 3, 2011 at 5:37 AM, Erik Rigtorp <erik@rigtorp.com> wrote:
It's only a view of the array, no copying is done. Though some operations like np.std() will copy the array, but that's more of a bug. In general It's hard to imagine any speedup gains by copying a 10GB array.
I don't think that np.std makes a copy of the input data if the input is an array. If the input is, for example, a list, then an array is created.
When I tried it on a big array, it tried to allocate a huge amount of memory. As I said it's probably a bug.
On Mon, Jan 3, 2011 at 7:41 AM, Erik Rigtorp <erik@rigtorp.com> wrote:
On Mon, Jan 3, 2011 at 10:36, Keith Goodman <kwgoodman@gmail.com> wrote:
On Mon, Jan 3, 2011 at 5:37 AM, Erik Rigtorp <erik@rigtorp.com> wrote:
It's only a view of the array, no copying is done. Though some operations like np.std() will copy the array, but that's more of a bug. In general It's hard to imagine any speedup gains by copying a 10GB array.
I don't think that np.std makes a copy of the input data if the input is an array. If the input is, for example, a list, then an array is created.
When I tried it on a big array, it tried to allocate a huge amount of memory. As I said it's probably a bug.
Yes, that would be a big bug. np.std does have to initialize the output array. If the window size is small compared to arr.shape[axis] then the memory taken by the output array is of the same order as that of the input array. Could that be what you are seeing?
a = np.arange(10)
Small window, output array shape (8,):
rolling_window(a, 2) array([[0, 1], [1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9]])
Big window, output array shape (2,):
rolling_window(a, 9) array([[0, 1, 2, 3, 4, 5, 6, 7, 8], [1, 2, 3, 4, 5, 6, 7, 8, 9]])
On Mon, Jan 3, 2011 at 10:52, Keith Goodman <kwgoodman@gmail.com> wrote:
On Mon, Jan 3, 2011 at 7:41 AM, Erik Rigtorp <erik@rigtorp.com> wrote:
On Mon, Jan 3, 2011 at 10:36, Keith Goodman <kwgoodman@gmail.com> wrote:
On Mon, Jan 3, 2011 at 5:37 AM, Erik Rigtorp <erik@rigtorp.com> wrote:
It's only a view of the array, no copying is done. Though some operations like np.std() will copy the array, but that's more of a bug. In general It's hard to imagine any speedup gains by copying a 10GB array.
I don't think that np.std makes a copy of the input data if the input is an array. If the input is, for example, a list, then an array is created.
When I tried it on a big array, it tried to allocate a huge amount of memory. As I said it's probably a bug.
Yes, that would be a big bug.
np.std does have to initialize the output array. If the window size is small compared to arr.shape[axis] then the memory taken by the output array is of the same order as that of the input array. Could that be what you are seeing?
a = np.arange(10)
Small window, output array shape (8,):
rolling_window(a, 2) array([[0, 1], [1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9]])
Big window, output array shape (2,):
rolling_window(a, 9) array([[0, 1, 2, 3, 4, 5, 6, 7, 8], [1, 2, 3, 4, 5, 6, 7, 8, 9]])
No the array was (500,2000) and i did np.std(rolling_window(a,252),-1) and it started to allocate > 2GB.
On Fri, Dec 31, 2010 at 8:29 PM, Erik Rigtorp <erik@rigtorp.com> wrote:
Implementing moving average, moving std and other functions working over rolling windows using python for loops are slow. This is a effective stride trick I learned from Keith Goodman's <kwgoodman@gmail.com> Bottleneck code but generalized into arrays of any dimension. This trick allows the loop to be performed in C code and in the future hopefully using multiple cores.
I like using strides for moving window functions. The one downside I found is that it is slow when window * (arr.shape[axis] - window) is large:
a = np.random.rand(1000000) b = rolling_window(a, 5000)
import bottleneck as bn timeit bn.nanmean(b, axis=1) 1 loops, best of 3: 7.1 s per loop timeit bn.move_nanmean(a, window=5000, axis=0) 100 loops, best of 3: 7.99 ms per loop
On 12/31/2010 06:29 PM, Erik Rigtorp wrote:
Hi,
Implementing moving average, moving std and other functions working over rolling windows using python for loops are slow. This is a effective stride trick I learned from Keith Goodman's <kwgoodman@gmail.com> Bottleneck code but generalized into arrays of any dimension. This trick allows the loop to be performed in C code and in the future hopefully using multiple cores.
An alternative is to go straight to C, with a cython interface. If you look in the num/src subdirectory of http://currents.soest.hawaii.edu/hgstage/hgwebdir.cgi/pycurrents/ you will find this approach, labeled "ringbuf" and "runstats". See pycurrents/setup.py, and its driver, pycurrents/runsetup.py, to see how runstats is presently being built. Instead of calculating statistics independently each time the window is advanced one data point, the statistics are updated. I have not done any benchmarking, but I expect this approach to be quick. The code is old; I have not tried to update it to take advantage of cython's advances over pyrex. If I were writing it now, I might not bother with the C level at all; it could all be done in cython, probably with no speed penalty, and maybe even with reduced overhead. Eric
On Mon, Jan 3, 2011 at 11:26, Eric Firing <efiring@hawaii.edu> wrote:
Instead of calculating statistics independently each time the window is advanced one data point, the statistics are updated. I have not done any benchmarking, but I expect this approach to be quick.
This might accumulate numerical errors. But could be fine for many applications.
The code is old; I have not tried to update it to take advantage of cython's advances over pyrex. If I were writing it now, I might not bother with the C level at all; it could all be done in cython, probably with no speed penalty, and maybe even with reduced overhead.
No doubt this would be faster, I just wanted to offer a general way to this in NumPy.
On Mon, Jan 3, 2011 at 5:32 PM, Erik Rigtorp <erik@rigtorp.com> wrote:
On Mon, Jan 3, 2011 at 11:26, Eric Firing <efiring@hawaii.edu> wrote:
Instead of calculating statistics independently each time the window is advanced one data point, the statistics are updated. I have not done any benchmarking, but I expect this approach to be quick.
This might accumulate numerical errors. But could be fine for many applications.
The code is old; I have not tried to update it to take advantage of cython's advances over pyrex. If I were writing it now, I might not bother with the C level at all; it could all be done in cython, probably with no speed penalty, and maybe even with reduced overhead.
No doubt this would be faster, I just wanted to offer a general way to this in NumPy. _______________________________________________
BTW, some of these operations can be done using scipy's ndimage - right ? Any comments ? How does the performance compare ? ndimage might have more options regarding edge handling, or ? Cheers, Sebastian Haase
On Tue, Jan 4, 2011 at 8:06 AM, Sebastian Haase <seb.haase@gmail.com> wrote:
On Mon, Jan 3, 2011 at 5:32 PM, Erik Rigtorp <erik@rigtorp.com> wrote:
On Mon, Jan 3, 2011 at 11:26, Eric Firing <efiring@hawaii.edu> wrote:
Instead of calculating statistics independently each time the window is advanced one data point, the statistics are updated. I have not done any benchmarking, but I expect this approach to be quick.
This might accumulate numerical errors. But could be fine for many applications.
The code is old; I have not tried to update it to take advantage of cython's advances over pyrex. If I were writing it now, I might not bother with the C level at all; it could all be done in cython, probably with no speed penalty, and maybe even with reduced overhead.
No doubt this would be faster, I just wanted to offer a general way to this in NumPy. _______________________________________________
BTW, some of these operations can be done using scipy's ndimage - right ? Any comments ? How does the performance compare ? ndimage might have more options regarding edge handling, or ?
Take a look at the moving window function in the development version of the la package: https://github.com/kwgoodman/la/blob/master/la/farray/mov.py Many of the moving window functions offer three calculation methods: filter (ndimage), strides (the strides trick discussed in this thread), and loop (a simple python loop). For example:
a = np.random.rand(500,2000) timeit la.farray.mov_max(a, window=252, axis=-1, method='filter') 1 loops, best of 3: 336 ms per loop timeit la.farray.mov_max(a, window=252, axis=-1, method='strides') 1 loops, best of 3: 609 ms per loop timeit la.farray.mov_max(a, window=252, axis=-1, method='loop') 1 loops, best of 3: 638 ms per loop
No one method is best for all situations. That is one of the reasons I started the Bottleneck package. I figured Cython could beat them all.
On Tue, Jan 4, 2011 at 8:14 AM, Keith Goodman <kwgoodman@gmail.com> wrote:
On Tue, Jan 4, 2011 at 8:06 AM, Sebastian Haase <seb.haase@gmail.com> wrote:
On Mon, Jan 3, 2011 at 5:32 PM, Erik Rigtorp <erik@rigtorp.com> wrote:
On Mon, Jan 3, 2011 at 11:26, Eric Firing <efiring@hawaii.edu> wrote:
Instead of calculating statistics independently each time the window is advanced one data point, the statistics are updated. I have not done any benchmarking, but I expect this approach to be quick.
This might accumulate numerical errors. But could be fine for many applications.
The code is old; I have not tried to update it to take advantage of cython's advances over pyrex. If I were writing it now, I might not bother with the C level at all; it could all be done in cython, probably with no speed penalty, and maybe even with reduced overhead.
No doubt this would be faster, I just wanted to offer a general way to this in NumPy. _______________________________________________
BTW, some of these operations can be done using scipy's ndimage - right ? Any comments ? How does the performance compare ? ndimage might have more options regarding edge handling, or ?
Take a look at the moving window function in the development version of the la package:
https://github.com/kwgoodman/la/blob/master/la/farray/mov.py
Many of the moving window functions offer three calculation methods: filter (ndimage), strides (the strides trick discussed in this thread), and loop (a simple python loop).
For example:
a = np.random.rand(500,2000) timeit la.farray.mov_max(a, window=252, axis=-1, method='filter') 1 loops, best of 3: 336 ms per loop timeit la.farray.mov_max(a, window=252, axis=-1, method='strides') 1 loops, best of 3: 609 ms per loop timeit la.farray.mov_max(a, window=252, axis=-1, method='loop') 1 loops, best of 3: 638 ms per loop
No one method is best for all situations. That is one of the reasons I started the Bottleneck package. I figured Cython could beat them all.
I added four new function to Bottleneck: move_min, move_max, move_nanmin, move_nanmax. They are much faster than using SciPy's ndimage.maximum_filter1d or the strides trick:
a = np.random.rand(500,2000) timeit la.farray.mov_max(a, window=252, axis=-1, method='filter') # ndimage 1 loops, best of 3: 336 ms per loop timeit bn.move_max(a, window=252, axis=-1) # bottleneck 100 loops, best of 3: 14.1 ms per loop
That looks too good to be true. Are the outputs the same?
a1 = la.farray.mov_max(a, window=252, axis=-1, method='filter') a2 = bn.move_max(a, window=252, axis=-1) np.testing.assert_array_almost_equal(a1, a2)
Yes.
participants (4)
-
Eric Firing -
Erik Rigtorp -
Keith Goodman -
Sebastian Haase