Creation of MxM matrix: Python version much slower than Matlab version

Question

I have written the same code in Python (NumPy) and in Matlab, and I tried to use the same structure for both language and follow the same procedure. Now my problem is that when I run my code in Python it's very very slow but it runs fast in Matlab.

How can I improve my Python code? I want to work on the Python version of my code.

In the code I am trying to create an MxM matrix (M can vary between 100 to 1000) and trying to fill this matrix 16 times, each time I add the previous value of M[i,j] with the current value. The operation inside the code is simple just some add and divide, and now I wonder why Matlab performs better than Python.

Python:

import numpy as np
import matplotlib.pyplot as plt

from datetime import datetime
startTime = datetime.now()

bb = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]
bb = ['{0:04b}'.format(x) for x in bb]


step_prob=0.01

step = int(1/step_prob)

out_t = np.zeros((step,step), dtype=float)
out_f = np.zeros((step,step), dtype=float)
s = np.zeros((step,step), dtype=float)
improve = np.empty((step,step))

for ii in range(1, step+1):
    print ii
    pti=ii*step_prob
    for  jj in range(1, step+1):
        pfi=jj*step_prob
        improve[ii-1,jj-1]=0
        si=(1-pfi+pti)/2
        for k in range(len(bb)):
            f1 = int(bb[k][0])
            f2 = int(bb[k][1])  
            f3 = int(bb[k][2])
            f4 = int(bb[k][3])
            for i in range(1, step+1):
                for j in range(1, step+1):
                    ptj=i*step_prob
                    pfj=j*step_prob
                    out_t[i-1,j-1] = f1*pti*ptj+f2*pti*(1-ptj)+f3*(1-pti)*ptj+f4*(1-pti)*(1-ptj)
                    out_f[i-1,j-1] = f1*pfi*pfj+f2*pfi*(1-pfj)+f3*(1-pfi)*pfj+f4*(1-pfi)*(1-pfj)
                    sj=(1-pfj+ptj)/2;
                    s[i-1,j-1] = ((1-out_f[i-1,j-1]+out_t[i-1,j-1])/2)-max(si,sj);

#         temp=s*(s>0)*tril(ones(size(s)));
            temp = s*(s>0)*np.tril(np.ones((len(s),len(s))).astype(int))
            improve[ii-1,jj-1]=improve[ii-1,jj-1]+sum(sum(temp))

im = plt.imshow(improve*np.tril(np.ones((len(improve),len(improve))).astype(int)),interpolation='bilinear', aspect='auto') # s*[s>0]
# im = plt.imshow(improve, interpolation='bilinear', aspect='auto') # s*[s>0]
plt.colorbar(im, orientation='vertical')
plt.show()

Matlab:

close all
clear all
clc


bb = de2bi([0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15],4,'left-msb');
% bb = de2bi([1],4);

step_prob=0.01;
for ii=1:1*(1/step_prob)
    ii
    pti=ii*step_prob;
    for jj=1:1*(1/step_prob)
        pfi=jj*step_prob;
        improve(ii,jj)=0;
        si=(1-pfi+pti)/2;
        for k=1:size(bb,1)
            f1=bb(k,1);
            f2=bb(k,2);
            f3=bb(k,3);
            f4=bb(k,4);
            for i=1:1*(1/step_prob)
                for j=1:1*(1/step_prob)
                    ptj=i*step_prob;
                    pfj=j*step_prob;
                    out_t(i,j)=f1*pti*ptj+f2*pti*(1-ptj)+f3*(1-pti)*ptj+f4*(1-pti)*(1-ptj);
                    out_f(i,j)=f1*pfi*pfj+f2*pfi*(1-pfj)+f3*(1-pfi)*pfj+f4*(1-pfi)*(1-pfj);
                    sj=(1-pfj+ptj)/2;
                    s(i,j)=((1-out_f(i,j)+out_t(i,j))/2)-max(si,sj);
                end
            end
            temp=s.*(s>0).*tril(ones(size(s)));
            improve(ii,jj)=improve(ii,jj)+sum(sum(temp));
            clear f1 f2 f3 f4
        end
    end
end
scale=step_prob:step_prob:1;
imagesc(scale,scale,improve.*tril(ones(size(improve))))
xlabel('Pfj')
ylabel('Ptj')
axis xy

Answer 1

You use very short variable names, so its hard to follow what your code is doing.

For speed in numpy, you need to use vector operations. So instead of writing

for ii in range(step):
    out[ii] = ii * x + b

You should write:

out = numpy.arange(step)*x + b

Individually accessing the elements of numpy array is actually quite slow. You need to write your code to do everything in vector operations.

Here is my rewrite of your inner loop to be vectorized. You could vectorize the outer loop as well, but this should make the biggest difference.

        f1 = int(bb[k][0])
        f2 = int(bb[k][1])  
        f3 = int(bb[k][2])
        f4 = int(bb[k][3])

        ptj = np.arange(step_prob, 1.0 + step_prob, step_prob)[:, None]
        pfj = np.arange(step_prob, 1.0 + step_prob, step_prob)[None, :]

        out_t = f1*pti*ptj+f2*pti*(1-ptj)+f3*(1-pti)*ptj+f4*(1-pti)*(1-ptj)
        out_f = f1*pfi*pfj+f2*pfi*(1-pfj)+f3*(1-pfi)*pfj+f4*(1-pfi)*(1-pfj)
        sj = (1 - pfj + ptj) / 2

        the_max = np.maximum(sj, [si])
        s = ((1-out_f+out_t)/2)-the_max

        temp = s*(s>0)*np.tril(np.ones((len(s),len(s))).astype(int))
        improve[ii-1,jj-1]=improve[ii-1,jj-1]+temp.sum()

Answer 2

My observations, in terms of implementing this in Python (without regards to improve the algorithm, that will be later):

1. You are doing a LOT of unnecessary type conversions in Python that you aren't doing in MATLAB. Type conversions are always slow, so this is going to hurt your performance a lot.
2. You shouldn't do for k in range(len(bb)):, always iterate over the sequence directly. You can use unpacking to get the values without needing to index in this case.
3. You are iterating over range(1, foo+1), then subtracting one to put it back into Python 0-based indexing. This is an unnecessary mathematical operation, better to do range(foo), or in your case enumerate(range(1, foo+1)). Better yet, pre-compute your values, then enumerate over those.
4. numpy.unpackbits is the equivalent of de2bi and will be much faster.
5. Python has in-place operations, which will be faster. So a += 1 instead of a = a + 1.
6. Numpy sums work of flattened arrays by default. Further, summing is a method of numpy arrays, which is simpler and probably faster.
7. Numpy arrays are floats by default if created using ones or zeros, but there is a dtype argument that can set it to whatever you want.
8. Numpy has a ones_like and zeros_like to create an array of the same shape as another. By default it also creates one with the same dtype, but you can override that.
9. PEP8 is the recommended style for python code.

Now, in terms of improving the algorithm:

1. si and sj are not very big (on the order of tens of megabytes, tops). You can vectorize those and re-use the values. In fact, they are identical, so you only need one.
2. You not only never use the other values of out_t and out_f besides the current one, you overwrite them repeatedly. So these are better off as scalars.
3. You can vectorize many of the values for the two innermost loops.
4. If you put the third loop as the outer loop, you can vectorize out_t and out_f.

So here is my partially-vectorized version:

bb0 = np.arange(16, dtype='uint8')
bb = np.unpackbits(bb0[None], axis=0)[-4:, :].T

step_prob = 0.01
step = int(1/step_prob)

pxx = np.linspace(step_prob, 1, step)
pxxt = pxx[:, None]

sx = (1-pxx+pxxt)/2

improve = np.zeros((step,step))
for f4, f3, f2, f1 in bb:
    out_x = f1*pxx*pxxt + f2*pxx*(1-pxxt) + f3*(1-pxx)*pxxt + f4*(1-pxx)*(1-pxxt)
    for ii, (out_t, sxi) in enumerate(izip(out_x, sx)):
        for jj, (out_f, si) in enumerate(izip(out_x, sxi)):
            s = (1-out_f[:, None]+out_t) - np.maximum(si, sxi)
            temp = s*(s>0)*np.tril(np.ones_like(s))
            improve[ii, jj] = temp.sum()

It is orders of magnitude faster than your version.