Creation of MxM matrix: Python version much slower than Matlab version 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&gt;0)*tril(ones(size(s)));
            temp = s*(s&gt;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&gt;0]
# im = plt.imshow(improve, interpolation='bilinear', aspect='auto') # s*[s&gt;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&gt;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
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&gt;0)*np.tril(np.ones((len(s),len(s))).astype(int))
        improve[ii-1,jj-1]=improve[ii-1,jj-1]+temp.sum()