0

I'm trying to analyze some data, for which I need to calculate a quantity involving a double sum. The python code looks like this:

import numpy as np

tmax = 1000
array = np.random.rand(tmax)

tstart = 500

meanA = np.mean(array[tstart:])

quantity = np.zeros(tmax-tstart)
for f in range(1,tmax-tstart,1):
    count = 0
    integrand = 0

    for ff in range(tstart+1,tmax-f):
            count += 1
            dAt = array[ff] - meanA
            dAtt = array[ff:ff+f+1] - meanA
            integrand += np.sum(dAt * dAtt)

    if count != 0:
            integrand /= f*count
            quantity[f] = integrand

This takes roughly 1.5s to run. This is 10 times more time than MATLAB takes to do the same computation:

tic;
tmax = 1000;
tstart = 500;

array = rand(1,tmax);
meanA = mean(array(tstart:tmax));
quantity = zeros(1,tmax);

for f=1:tmax-tstart
    integrand = 0;
    count = 0;
    for ff=tstart:tmax-f
    count = count + 1;
    dAt = array(ff)-meanA;
    dAtt = array(ff:ff+f)-meanA;
    integrand = integrand + sum(dAt*dAtt);
end
integrand = integrand/(f*count);
autocorr(f) = integrand;
end

toc

Outputting:

>> speedTest
Elapsed time is 0.096789 seconds.

Why is my python script so slow? How can I make it run as fast as the MATLAB script? (And yes, I have to do this in python for a number of other reasons)

Mind that the real data corresponds to an array size of >10,000 elements, so the time different becomes radically large, as the number of flops scale quadratically with the number of elements.

EDIT:

I tried the same without numpy (except for random number generation), using solely lists:

import numpy as np

tmax = 1000
array = np.random.rand(tmax)

array = list(array)

tstart = 500

meanA = sum((array[tstart:]))/len(array[tstart:])

quantity = [0] * (tmax-tstart)
for f in range(1,tmax-tstart,1):
    count = 0
    integrand = 0

    for ff in range(tstart+1,tmax-f):
            count += 1
            dAt = array[ff] - meanA
            dAtt = array[ff:ff+f+1] - meanA
            try:
                    integrand += sum([dAt * i for i in dAtt])
            except:
                    integrand += dAt * dAtt
    if count != 0:
            integrand /= f*count
            quantity[f] = integrand

This results in:

$ time python3 speedAutoCorr2.py

real    0m6.510s
user    0m6.731s
sys     0m0.123s

which is even worse than the case with numpy.

sodiumnitrate
  • 2,899
  • 6
  • 30
  • 49
  • Your python script calculates that `count !=0` part on every iteration of inner loop, while MATLAB script does it in outer loop. – Alexander Dmitriev Mar 17 '19 at 19:41
  • 2
    If you want to get good performance out of for-loop heavy `numpy` code, consider the `numba` JIT for `numpy`. `numpy` will be slow like this. Another possibility is Cython. Your other option is to translate your numpy approach into something that takes advantages of the very fast vectorized operations on `numpy` arrays, not using for-loops and indexing like this, or the routines exposed in `scipy`. – juanpa.arrivillaga Mar 17 '19 at 19:46
  • @AlexanderDmitriev that's a typo, the original code doesn't have it. I can verify that moving it to the outer for loop has no effect on speed. – sodiumnitrate Mar 17 '19 at 19:49
  • @juanpa.arrivillaga Thanks for the suggestions. It wouldn't surprise me if `numpy` was slower compared to `C` but `matlab`? I'll look into `cython`, but it would also be nice to know why `numpy` is slow with loops. – sodiumnitrate Mar 17 '19 at 19:52
  • 2
    MATLAB has a JIT, Python doesn’t (unless you explicitly use numba). Yes, MATLAB is much faster than Python. – Cris Luengo Mar 18 '19 at 02:44
  • Possible duplicate of [for loop in python is 10x slower than matlab](https://stackoverflow.com/questions/17242684/for-loop-in-python-is-10x-slower-than-matlab) – sodiumnitrate Mar 18 '19 at 14:36
  • Thanks everyone, I've decided to solve this by writing a small C program and running it from the jupyter notebook within my python code. I also realized that this question was asked before and discussed at length, so I marked it as a duplicate. – sodiumnitrate Mar 18 '19 at 14:37

0 Answers0