Numpy summation with sliding window is really slow

Question

Code:

shape = np.array([6, 6])
grid = np.array([x.ravel() for x in np.meshgrid(*[np.arange(x) for i, x in enumerate(shape)], indexing='ij')]).T
slices = [tuple(slice(box[i], box[i] + 2) for i in range(len(box))) for box in grid]
score = np.zeros((7,7,3))
column = np.random.randn(36, 12) #just for example

column
>> array([[  0,   1,   2,   3, ... 425, 426, 427, 428, 429, 430, 431]])

column = column.reshape((16, 3, 3, 3))
for i, window in enumerate(slices):
        score[window] += column[i]
score
>> array([[[0.000e+00, 1.000e+00, 2.000e+00],
    [3.000e+01, 3.200e+01, 3.400e+01],
    [9.000e+01, 9.300e+01, 9.600e+01], ... 
    [8.280e+02, 8.300e+02, 8.320e+02],
    [4.290e+02, 4.300e+02, 4.310e+02]]])

It works but last 2 lines take really much time as they will be in loop. The problem is that 'grid' variable contains an array of windows. And I don't now how to speed up the process.

Last 2 lines seem to be easily vectorizable... pure python for loops are really slow for iterating and summing.. vectorized operations fix exactly that — rafaelc, Aug 15 '18 at 21:27
Nope.. I meant summing np arrays directly.. Post a sample data, will make it easier to demonstrate — rafaelc, Aug 15 '18 at 21:30
by vectorization, he means you could directly add the arrays using something like np.sum(). — Anurag Reddy, Aug 15 '18 at 21:37
After some tweaks to get the code to work, I tried a `%%timeit on two parts. Setup averaged 5ms, the for loop 6ms. Not a drastic difference. Iteration over a range of 600+ will take time. — hpaulj, Aug 15 '18 at 21:46
You are summing (3,3) windows. Is this some sort of convolution? Or something that can solved with `as_strided`? — hpaulj, Aug 15 '18 at 21:52

hpaulj · Accepted Answer · 2018-08-17T22:48:19.667

Let's simplify the problem at bit - reduce the dimensions, and drop the final size 3 dimension:

In [265]: shape = np.array([4,4])
In [266]: grid = np.array([x.ravel() for x in np.meshgrid(*[np.arange(x) for i
     ...: , x in enumerate(shape)], indexing='ij')]).T
     ...: grid = [tuple(slice(box[i], box[i] + 3) for i in range(len(box))) fo
     ...: r box in grid]
     ...: 
     ...: 
In [267]: len(grid)
Out[267]: 16
In [268]: score = np.arange(36).reshape(6,6)
In [269]: X = np.array([score[x] for x in grid]).reshape(4,4,3,3)
In [270]: X
Out[270]: 
array([[[[ 0,  1,  2],
         [ 6,  7,  8],
         [12, 13, 14]],

        [[ 1,  2,  3],
         [ 7,  8,  9],
         [13, 14, 15]],

        [[ 2,  3,  4],
         [ 8,  9, 10],
         [14, 15, 16]],

        ....
        [[21, 22, 23],
         [27, 28, 29],
         [33, 34, 35]]]])

This is a moving window - one (3,3) array, shift over 1,..., shift down 1, etc

With as_strided is is possible to construct a view of the array, that consists of all these windows, but without actually copying values. Having worked with as_strided before I was able construct the equivalent strides as:

In [271]: score.shape
Out[271]: (6, 6)
In [272]: score.strides
Out[272]: (48, 8)
In [273]: ast = np.lib.stride_tricks.as_strided
In [274]: x=ast(score, shape=(4,4,3,3), strides=(48,8,48,8))
In [275]: np.allclose(X,x)
Out[275]: True

This could be extended to your (28,28,3) dimensions, and turned into the summation.

Generating such moving windows has been covered in previous SO questions. And it's also implemented in one of the image processing packages.

Adaptation for a 3 channel image,

In [45]: arr.shape
Out[45]: (6, 6, 3)
In [46]: arr.strides
Out[46]: (144, 24, 8)
In [47]: arr[:3,:3,0]
Out[47]: 
array([[ 0.,  1.,  2.],
       [ 6.,  7.,  8.],
       [12., 13., 14.]])

In [48]: x = ast(arr, shape=(4,4,3,3,3), strides=(144,24,144,24,8))
In [49]: x[0,0,:,:,0]
Out[49]: 
array([[ 0.,  1.,  2.],
       [ 6.,  7.,  8.],
       [12., 13., 14.]])

Since we are moving the window by one element at a time, the strides for x are easily derived form the source strides.

For 4x4 windows, just change the shape

x = ast(arr, shape=(3,3,4,4,3), strides=(144,24,144,24,8))

In Efficiently Using Multiple Numpy Slices for Random Image Cropping

@Divikar suggests using skimage

With the default step=1, the result is compatible:

In [55]: from skimage.util.shape import view_as_windows
In [63]: y = view_as_windows(arr,(4,4,3))
In [64]: y.shape
Out[64]: (3, 3, 1, 4, 4, 3)
In [69]: np.allclose(x,y[:,:,0])
Out[69]: True

Ok. I'm not familiar with as_strided function yet. What does the 'strides' parameter mean? — Nurislam Fazulzyanov, Aug 16 '18 at 05:57
`strides` is a basic attribute of a numpy array (along with `shape` and `dtype`). It's part of what makes this multidimensional code so powerful. But its use in `as_strided` is a rather advanced topic. — hpaulj, Aug 16 '18 at 06:20
I have updated my code a bit. I understand that it can your solution could be turned into summation but in my particular situation I want the result to be of 'score' shape, because I want to iteratively add 'columns' to 'score' variable. — Nurislam Fazulzyanov, Aug 16 '18 at 06:37

Numpy summation with sliding window is really slow

1 Answers1

Linked