Is convolution slower in Numpy than in Matlab?

Question

Convolution in Matlab appears to be twice as fast as convolution in Numpy.

Python code (takes 19 seconds on my machine):

import numpy as np
from scipy import ndimage
import time

img = np.ones((512,512,512))
kernel = np.ones((5,5,5))/125

start_time = time.time()
ndimage.convolve(img,kernel,mode='constant')
print "Numpy execution took ", (time.time() - start_time), "seconds"

Matlab code (takes 8.7 seconds on my machine):

img = ones(512,512,512);
kernel = ones(5,5,5) / 125;
tic
convn(img, kernel, 'same');
toc

The two give identical results.

Is there a way to improve Numpy to match or beat Matlab's performance here?

Interestingly, this factor or ~2 difference in runtimes is consistent at many input sizes.

Answer 1

What exact operation are you doing? There are a number of optimizations that ndimage provides if you don't need general N-d convolution.

For example, your current operation:

img = np.ones((512,512,512))
kernel = np.ones((5,5,5))/125
result = ndimage.convolve(img, kernel)

is equivalent to:

img = np.ones((512,512,512))
result = ndimage.uniform_filter(img, 5)

However, there's a large difference in execution speed:

In [22]: %timeit ndimage.convolve(img, kernel)
1 loops, best of 3: 25.9 s per loop

In [23]: %timeit ndimage.uniform_filter(img, 5)
1 loops, best of 3: 8.69 s per loop

The difference is caused by uniform_filter preforming a 1-d convolution along each axis, instead of a generic 3D convolution.

In cases where the kernel is symmetric, you can make these simplifications and get a significant speed up.

I'm not certain about convn, but often matlab's functions make these kind of optimizations behind-the-scenes if your input data matches certain criteria. Scipy more commonly uses one-algorithm-per-function, and expects the user to know which one to pick in which case.

---

You mentioned a "Laplacian of Gaussian" filter. I always get a touch confused with terminology here.

I think what you want in terms of ndimage functions is either scipy.ndimage.gaussian_laplace or scipy.ndimage.gaussian_filter with order=2 (which filters by the second derivative of the gaussian).

At any rate, both simplify the operation down to a 1-d convolution over each axis, which should give a significant (2-3x) speedup.