Improving Numpy Performance Improving Numpy Performance numpy numpy

Improving Numpy Performance


The code in scipy for doing 2d convolutions is a bit messy and unoptimized. See http://svn.scipy.org/svn/scipy/trunk/scipy/signal/firfilter.c if you want a glimpse into the low-level functioning of scipy.

If all you want is to process with a small, constant kernel like the one you showed, a function like this might work:

def specialconvolve(a):    # sorry, you must pad the input yourself    rowconvol = a[1:-1,:] + a[:-2,:] + a[2:,:]    colconvol = rowconvol[:,1:-1] + rowconvol[:,:-2] + rowconvol[:,2:] - 9*a[1:-1,1:-1]    return colconvol

This function takes advantage of the separability of the kernel like DarenW suggested above, as well as taking advantage of the more optimized numpy arithmetic routines. It's over 1000 times faster than the convolve2d function by my measurements.


For the particular example 3x3 kernel, I'd observe that

1  1  11 -8  11  1  1  1  1  1     0  0  0= 1  1  1  +  0 -9  0  1  1  1     0  0  0

and that the first of these is factorable - it can be convoluted by convolving (1 1 1) for each row, and then again for each column. Then subtract nine times the original data. This may or may not be faster, depending on whether the scipy programmers made it smart enough to automatically do this. (I haven't checked in a while.)

You probably want to do more interesting convolutions, where factoring may or may not be possible.


Before going to say C with ctypes, I'd suggest running a standalone convolve in C, to see where the limit is.
Similarly for CUDA, cython, scipy.weave ...

Added 7feb: convolve33 8-bit data with clipping takes ~ 20 clock cycles per point,2 clock cycles per mem access, on my mac g4 pcc with gcc 4.2. Your mileage will vary.

A couple of subtleties:

  • do you care about correct clipping to 0..255 ? np.clip() is slow,cython etc. don't know.
  • Numpy/scipy may need memory for temps the size of A (so keep 2*sizeof(A) < cache size).
    If your C code, though, does a running update inplace, that's half the mem but a different algorithm.

By the way, google theano convolve =>"A convolution op that should mimic scipy.signal.convolve2d, but faster! In development"