Почему этот двойной цикл значительно медленнее в Python по сравнению с Matlab?

Question

Я пытаюсь проанализировать некоторые данные, для которых мне нужно вычислить количество, включающее двойную сумму. Код Python выглядит следующим образом:

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

Для запуска требуется примерно 1.5s. Это в 10 раз больше , чем MATLAB требует для того же вычисления:

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

Выведение:

>> speedTest
Elapsed time is 0.096789 seconds.

Почему мой скрипт на python такой медленный? Как я могу заставить его работать так же быстро, как скрипт MATLAB? (И да, я должен сделать это в Python по ряду других причин)

Помните, что реальные данные соответствуют размеру массива >10,000 элементов, поэтому разность времени становится радикально большой, так как число флопов масштабируется квадратично с количеством элементов.

EDIT:

Я попробовал то же самое без numpy (кроме генерации случайных чисел), используя только списки:

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

В результате:

$ time python3 speedAutoCorr2.py

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

, что еще хуже, чем в случае с numpy.