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: 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`](http://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.filters.gaussian_laplace.html#scipy.ndimage.filters.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.

