Matlab vs Python 2D convolution performance

Question

I like to prototype algorithms in Matlab, but I have the requirement of putting them on a server that also runs quite a bit of Python code. Hence I quickly converted the code to Python and compared the two. The Matlab implementation runs ~1000 times faster (from timing function calls - no profiling). Anyone know off hand why the performance of Python is so slow?

Matlab

% init random data
w = 800;
h = 1200;
hmap = zeros(w,h);

npts = 250;
for i=1:npts
    hmap(randi(w),randi(h)) = hmap(randi(w),randi(h))+1;
end

% Params
disksize = 251;
nBreaks = 25;
saturation = .9;
floorthresh =.05;

fh = fspecial('gaussian', disksize, disksize/7);
hmap = conv2(hmap, fh, 'same');

% Scaling, paritioning etc
hmap = hmap/(max(max(hmap)));
hmap(hmap<floorthresh) = 0;
hmap = round(nBreaks * hmap)/nBreaks;
hmap = hmap * (1/saturation);

% Show the image
imshow(hmap, [0,1])
colormap('jet')

Python

import numpy as np
from scipy.signal import convolve2d as conv2

# Test data parameters
w = 800
h = 1200
npts = 250

# generate data
xvals = np.random.randint(w, size=npts)
yvals = np.random.randint(h, size=npts)

# Heatmap parameters
gaussianSize = 250
nbreaks = 25

# Preliminary function definitions
def populateMat(w, h, xvals, yvals):
    container = np.zeros((w,h))

    for idx in range(0,xvals.size):
        x = xvals[idx]
        y = yvals[idx]
        container[x,y] += 1

    return container

def makeGaussian(size, fwhm):

    x = np.arange(0, size, 1, float)
    y = x[:,np.newaxis]

    x0 = y0 = size // 2

    return np.exp(-4*np.log(2) * ((x-x0)**2 + (y-y0)**2) / fwhm**2)


# Create the data matrix
dmat = populateMat(w,h,xvals,yvals)
h = makeGaussian(gaussianSize, fwhm=gaussianSize/2)

# Convolve
dmat2 = conv2(dmat, h, mode='same')

# Scaling etc
dmat2 = dmat2 / dmat2.max()
dmat2 = np.round(nbreaks*dmat2)/nbreaks
# Show
imshow(dmat2)

Use numpy for the convolution. Have a look at http://stackoverflow.com/questions/16121269/2d-convolution-in-python-similar-to-matlabs-conv2 — kezzos, Nov 12 '14 at 09:48
are you saying that the implementation in `scipy.signal.convolve2d` is the culprit? Will investigate. — Joe, Nov 12 '14 at 09:50
I read the post and it seems to say that the only reason to use numpy is to avoid the scipy dependency. I don't have that restriction. Tried ndimage.convolve but this throws a MemoryError. Will keep trying the differnt options suggested in that thread. — Joe, Nov 12 '14 at 09:59
It is 251x251 in Matlab and 250x250 in Python - Python is massively slower. Sure it is big but Python is being massively outperformed. — Joe, Nov 12 '14 at 10:01
That would mean you need to calculate a kernal of 250x250 pixels, for each pixel in your input image: In pixels this is 60 trillion for your image — kezzos, Nov 12 '14 at 10:07
Matlab: ~0.05 Seconds (using tic-toc) , Python: ~160 (using time.time()). Thats 3200 times slower! I know I am being demanding of the convolve operation but what I am wondering is why Matlab is capable and Python is not? — Joe, Nov 12 '14 at 10:09
My guess is that matlab uses an optimized algorithm for the Gaussian filter, exploiting separability of the the filter and possibly using a recursive approximation. Try using separability under Python. — , Nov 12 '14 at 10:11
Ok, thats interesting. I have seen a few discussions on separability. I guess I need to read up more and figure out how to implement in Python. Unless anyone knows of a lib? — Joe, Nov 12 '14 at 10:12
This is trivial: define two 1D Gaussian, horizontal and vertical and apply them in sequence. — , Nov 12 '14 at 10:13

score 1 · Accepted Answer · answered Nov 12 '14 at 10:29

Ok, problem solved for me thanks to suggestion from @Yves Daust's comments;

The filter scipy.ndimage.filters.gaussian_filter utilises the separability of the kernel and reduces the running time to within a single order of magnitude of the matlab implementation.

import numpy as np
from scipy.ndimage.filters import gaussian_filter as gaussian


# Test data parameters
w = 800
h = 1200
npts = 250

# generate data
xvals = np.random.randint(w, size=npts)
yvals = np.random.randint(h, size=npts)

# Heatmap parameters
gaussianSize = 250
nbreaks = 25

# Preliminary function definitions
def populateMat(w, h, xvals, yvals):
    container = np.zeros((w,h))

    for idx in range(0,xvals.size):
        x = xvals[idx]
        y = yvals[idx]
        container[x,y] += 1

    return container


# Create the data matrix
dmat = populateMat(w,h,xvals,yvals)

# Convolve
dmat2 = gaussian(dmat, gaussianSize/7)

# Scaling etc
dmat2 = dmat2 / dmat2.max()
dmat2 = np.round(nbreaks*dmat2)/nbreaks

# Show
imshow(dmat2)

Matlab vs Python 2D convolution performance

Matlab

Python

1 Answers1