CODEWITHSUNDEEP
pythonperformancematlabnumpyconvolution
manimino
manimino

Reputation: 1365

Is convolution slower in Numpy than in Matlab?

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.

Upvotes: 5

Views: 6895

Answers (2)

Joe Kington
Joe Kington

Reputation: 284612

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.

Upvotes: 8

unutbu
unutbu

Reputation: 879411

Not an anwer; just a comment:

Before you compare performance, you need to make sure the two functions are returning the same result:

If Matlab's convn returns the same result as Octave's convn, then convn is different than ndimage.convolve:

octave> convn(ones(3,3), ones(2,2))
ans =

   1   2   2   1
   2   4   4   2
   2   4   4   2
   1   2   2   1

In [99]: ndimage.convolve(np.ones((3,3)), np.ones((2,2)))
Out[99]: 
array([[ 4.,  4.,  4.],
       [ 4.,  4.,  4.],
       [ 4.,  4.,  4.]])

ndimage.convolve has other modes, 'reflect','constant','nearest','mirror', 'wrap', but none of these match convn's default ("full") behavior.

For 2D arrays, scipy.signal.convolve2d is faster than scipy.signal.convolve.

For 3D arrays, scipy.signal.convolve seems to have the same behavior as convn(A,B):

octave> x = convn(ones(3,3,3), ones(2,2,2))
x =

ans(:,:,1) =

   1   2   2   1
   2   4   4   2
   2   4   4   2
   1   2   2   1

ans(:,:,2) =

   2   4   4   2
   4   8   8   4
   4   8   8   4
   2   4   4   2

ans(:,:,3) =

   2   4   4   2
   4   8   8   4
   4   8   8   4
   2   4   4   2

ans(:,:,4) =

   1   2   2   1
   2   4   4   2
   2   4   4   2
   1   2   2   1

In [109]: signal.convolve(np.ones((3,3,3), dtype='uint8'), np.ones((2,2,2), dtype='uint8'))
Out[109]: 
array([[[1, 2, 2, 1],
        [2, 4, 4, 2],
        [2, 4, 4, 2],
        [1, 2, 2, 1]],

       [[2, 4, 4, 2],
        [4, 8, 8, 4],
        [4, 8, 8, 4],
        [2, 4, 4, 2]],

       [[2, 4, 4, 2],
        [4, 8, 8, 4],
        [4, 8, 8, 4],
        [2, 4, 4, 2]],

       [[1, 2, 2, 1],
        [2, 4, 4, 2],
        [2, 4, 4, 2],
        [1, 2, 2, 1]]], dtype=uint8)

Note that np.ones((n,m,p)) creates a float array by default. Matlab ones(n,m,p) appears to create an array of ints. To make a good comparison, you should try to match the dtype of the numpy arrays to the type of the Matlab matrices.

Upvotes: 3

Related Questions