Różnica w wydajności między numpy a matlab

Question

Obliczam algorytm backpropagation dla rzadkiego autoencodera. Zaimplementowałem go w python przy użyciu numpy i matlab. Kod jest prawie taki sam, ale wydajność jest zupełnie inna. Czas potrzebny na wykonanie zadania matlab to 0.252454 sekundy, podczas gdy numpy 0.973672151566, czyli prawie cztery razy więcej. Wywołasz ten kod kilka razy później w problemie minimalizacji, więc ta różnica prowadzi do kilku minut opóźnienia między implementacjami. Czy to normalne zachowanie? Jak mogę poprawić wydajność w numpy?

realizacja Numpy:

Sparse.rho jest parametrem strojenia sparse.nodes to liczba węzłów w ukrytej warstwie (25) sparse.input (64) liczby węzłów w warstwie wejściowego theta1 i theta2 są macierzami wag odpowiednio dla pierwszej i drugiej warstwy o wymiarach 25x64 i 64x25, m jest równe 10000, rhoest ma wymiar (25,), x ma wymiar 10000x64, a3 10000x64 i a2 10000x25.

UPDATE: Wprowadziłem zmiany w kodzie po niektórych pomysłach odpowiedzi. Wydajność jest teraz numpy: 0.65 vs matlab: 0.25. Realizacja

partial_j1 = np.zeros(sparse.theta1.shape) 
partial_j2 = np.zeros(sparse.theta2.shape) 
partial_b1 = np.zeros(sparse.b1.shape) 
partial_b2 = np.zeros(sparse.b2.shape) 
t = time.time() 

delta3t = (-(x-a3)*a3*(1-a3)).T 

for i in range(m): 

    delta3 = delta3t[:,i:(i+1)] 
    sum1 = np.dot(sparse.theta2.T,delta3) 
    delta2 = (sum1 + sum2) * a2[i:(i+1),:].T* (1 - a2[i:(i+1),:].T) 
    partial_j1 += np.dot(delta2, a1[i:(i+1),:]) 
    partial_j2 += np.dot(delta3, a2[i:(i+1),:]) 
    partial_b1 += delta2 
    partial_b2 += delta3 

print "Backprop time:", time.time() -t 

Matlab:

tic 
for i = 1:m 

    delta3 = -(data(i,:)-a3(i,:)).*a3(i,:).*(1 - a3(i,:)); 
    delta3 = delta3.'; 
    sum1 = W2.'*delta3; 
    sum2 = beta*(-sparsityParam./rhoest + (1 - sparsityParam) ./ (1.0 - rhoest)); 
    delta2 = (sum1 + sum2) .* a2(i,:).' .* (1 - a2(i,:).'); 
    W1grad = W1grad + delta2* a1(i,:); 
    W2grad = W2grad + delta3* a2(i,:); 
    b1grad = b1grad + delta2; 
    b2grad = b2grad + delta3; 
end 
toc 

Answer 1

Byłoby źle powiedzieć "Matlab jest zawsze szybszy niż NumPy" lub vice versa. Często ich wydajność jest porównywalna. Podczas korzystania z NumPy, aby uzyskać dobrą wydajność, należy pamiętać, że szybkość NumPy pochodzi od wywołania funkcji napisanych w języku C/C++/Fortran. Działa to dobrze, gdy stosuje się te funkcje do całych tablic. Ogólnie rzecz biorąc, uzyskuje się gorszą wydajność po wywołaniu funkcji NumPy na mniejszych tablicach lub skalarach w pętli Pythona.

Co jest nie tak z pytaniem o Python? Każda iteracja za pośrednictwem pętli Pythona to wywołanie metody next. Każde użycie indeksowania [] jest wywołaniem metody __getitem__. Każdy numer += to połączenie z numerem __iadd__. Każdy atrybut z kropką (taki jak np. np.dot) obejmuje wywołania funkcji. Te funkcje wywołują , co przyczynia się do znacznego ograniczenia szybkości. Haczyki te dają Pythonowi ekspresywną siłę - indeksowanie ciągów oznacza coś innego niż indeksowanie na przykład dla dyktowania. Ta sama składnia, różne znaczenia. Magię osiąga się, nadając obiektom różne metody: __getitem__.

Ale ta ekspresyjna siła ma koszt w prędkości. Więc jeśli nie potrzebujesz całej dynamicznej ekspresji, aby uzyskać lepszą wydajność, spróbuj ograniczyć się do wywoływania funkcji NumPy na całych tablicach.

Usunąć pętlę for; jeśli to możliwe, stosuj równania "wektoryzowane". Na przykład, zamiast

for i in range(m): 
    delta3 = -(x[i,:]-a3[i,:])*a3[i,:]* (1 - a3[i,:])  

można obliczyć delta3 dla każdego i wszystko na raz:

delta3 = -(x-a3)*a3*(1-a3) 

Podczas gdy w for-loopdelta3 jest wektorem, używając vectorized równanie delta3 jest macierzą.

---

Niektóre obliczenia w for-loop nie zależą i i dlatego powinny być podnoszone na zewnątrz pętli. Na przykład, sum2 wygląda stałej:

sum2 = sparse.beta*(-float(sparse.rho)/rhoest + float(1.0 - sparse.rho)/(1.0 - rhoest)) 

---

Oto runnable z alternatywnym przykładem realizacji (alt) o kodzie (orig).

Moja odniesienia timeit pokazuje 6.8x poprawę szybkości:

In [52]: %timeit orig() 
1 loops, best of 3: 495 ms per loop 

In [53]: %timeit alt() 
10 loops, best of 3: 72.6 ms per loop 

---

import numpy as np 


class Bunch(object): 
    """ http://code.activestate.com/recipes/52308 """ 
    def __init__(self, **kwds): 
     self.__dict__.update(kwds) 

m, n, p = 10 ** 4, 64, 25 

sparse = Bunch(
    theta1=np.random.random((p, n)), 
    theta2=np.random.random((n, p)), 
    b1=np.random.random((p, 1)), 
    b2=np.random.random((n, 1)), 
) 

x = np.random.random((m, n)) 
a3 = np.random.random((m, n)) 
a2 = np.random.random((m, p)) 
a1 = np.random.random((m, n)) 
sum2 = np.random.random((p,)) 
sum2 = sum2[:, np.newaxis] 

def orig(): 
    partial_j1 = np.zeros(sparse.theta1.shape) 
    partial_j2 = np.zeros(sparse.theta2.shape) 
    partial_b1 = np.zeros(sparse.b1.shape) 
    partial_b2 = np.zeros(sparse.b2.shape) 
    delta3t = (-(x - a3) * a3 * (1 - a3)).T 
    for i in range(m): 
     delta3 = delta3t[:, i:(i + 1)] 
     sum1 = np.dot(sparse.theta2.T, delta3) 
     delta2 = (sum1 + sum2) * a2[i:(i + 1), :].T * (1 - a2[i:(i + 1), :].T) 
     partial_j1 += np.dot(delta2, a1[i:(i + 1), :]) 
     partial_j2 += np.dot(delta3, a2[i:(i + 1), :]) 
     partial_b1 += delta2 
     partial_b2 += delta3 
     # delta3: (64, 1) 
     # sum1: (25, 1) 
     # delta2: (25, 1) 
     # a1[i:(i+1),:]: (1, 64) 
     # partial_j1: (25, 64) 
     # partial_j2: (64, 25) 
     # partial_b1: (25, 1) 
     # partial_b2: (64, 1) 
     # a2[i:(i+1),:]: (1, 25) 
    return partial_j1, partial_j2, partial_b1, partial_b2 


def alt(): 
    delta3 = (-(x - a3) * a3 * (1 - a3)).T 
    sum1 = np.dot(sparse.theta2.T, delta3) 
    delta2 = (sum1 + sum2) * a2.T * (1 - a2.T) 
    # delta3: (64, 10000) 
    # sum1: (25, 10000) 
    # delta2: (25, 10000) 
    # a1: (10000, 64) 
    # a2: (10000, 25) 
    partial_j1 = np.dot(delta2, a1) 
    partial_j2 = np.dot(delta3, a2) 
    partial_b1 = delta2.sum(axis=1) 
    partial_b2 = delta3.sum(axis=1) 
    return partial_j1, partial_j2, partial_b1, partial_b2 

answer = orig() 
result = alt() 
for a, r in zip(answer, result): 
    try: 
     assert np.allclose(np.squeeze(a), r) 
    except AssertionError: 
     print(a.shape) 
     print(r.shape) 
     raise 

---

Wskazówka: Zauważ, że zostawiłem w komentarzach kształt wszystkim pośrednia tablice. Znajomość kształtu tablic pomogła mi zrozumieć, co robił twój kod. Kształt macierzy może pomóc w wyborze odpowiednich funkcji NumPy.A przynajmniej zwracanie uwagi na kształty może pomóc ci wiedzieć, czy operacja jest rozsądna. Na przykład, jeśli obliczyć

np.dot(A, B) 

i A.shape = (n, m) i B.shape = (m, p), następnie np.dot(A, B) będzie tablicą kształcie (n, p).

---

To może pomóc zbudować tablice w C_CONTIGUOUS zamówienie (przynajmniej, jeśli za pomocą np.dot). Nie może być tak jak z prędkością 3x się w ten sposób:

Poniżej x jest taka sama jak xf oprócz tego x jest C_CONTIGUOUS i xf jest F_CONTIGUOUS - i ta sama zależność dla y i yf.

import numpy as np 

m, n, p = 10 ** 4, 64, 25 
x = np.random.random((n, m)) 
xf = np.asarray(x, order='F') 

y = np.random.random((m, n)) 
yf = np.asarray(y, order='F') 

assert np.allclose(x, xf) 
assert np.allclose(y, yf) 
assert np.allclose(np.dot(x, y), np.dot(xf, y)) 
assert np.allclose(np.dot(x, y), np.dot(xf, yf)) 

%timeit benchmarki pokazują różnicę w prędkości:

In [50]: %timeit np.dot(x, y) 
100 loops, best of 3: 12.9 ms per loop 

In [51]: %timeit np.dot(xf, y) 
10 loops, best of 3: 27.7 ms per loop 

In [56]: %timeit np.dot(x, yf) 
10 loops, best of 3: 21.8 ms per loop 

In [53]: %timeit np.dot(xf, yf) 
10 loops, best of 3: 33.3 ms per loop 

---

chodzi benchmarkingu w Pythonie:

It can be misleading korzystania z różnicy w parach time.time() wywołań odniesienia prędkości kodu w języku Python. Musisz powtórzyć pomiar wiele razy. Lepiej wyłączyć automatyczny garbage collector. Ważne jest również mierzenie dużych rozpiętości czasu (na przykład co najmniej 10 sekund powtórzeń) w celu uniknięcia błędów spowodowanych słabą rozdzielczością w taktowaniu zegara i zmniejszenia znaczenia wezwań na wywołania time.time. Zamiast pisać cały ten kod samodzielnie, Python dostarcza ci timeit module. Zasadniczo używam tego na czas kawałków kodu, z wyjątkiem tego, że dzwonię do niego przez IPython terminal dla wygody.

Nie jestem pewien, czy ma to wpływ na testy porównawcze, ale należy pamiętać, że może to mieć znaczenie. W question I linked to, zgodnie z time.time, dwa kawałki kodu różniły się o czynnik 1,7x, podczas gdy testy porównawcze z użyciem timeit wykazały, że kawałki kodu przebiegały w zasadniczo identycznych ilościach czasu.

Answer 2

chciałbym zacząć inplace operacji, aby uniknąć przeznaczyć nowych macierzach za każdym razem:

partial_j1 += np.dot(delta2, a1[i,:].reshape(1,a1.shape[1])) 
partial_j2 += np.dot(delta3, a2[i,:].reshape(1,a2.shape[1])) 
partial_b1 += delta2 
partial_b2 += delta3 

można zastąpić to wyrażenie:

a1[i,:].reshape(1,a1.shape[1]) 

z prostszy i szybszy (dzięki Bi Rico):

a1[i:i+1] 

Również ta linia:

sum2 = sparse.beta*(-float(sparse.rho)/rhoest + float(1.0 - sparse.rho)/(1.0 - rhoest)) 

wydaje się być taka sama w każdej pętli, nie należy go przeliczyć.

I, prawdopodobnie niewielka optymalizacja, można zastąpić wszystkie wystąpienia x[i,:] z x[i].

Wreszcie, jeśli jesteś w stanie przeznaczyć m razy więcej pamięci, można śledzić unutbu sugestię i wektoryzacji pętlę:

for m in range(m): 
    delta3 = -(x[i]-a3[i])*a3[i]* (1 - a3[i]) 

z:

delta3 = -(x-a3)*a3*(1-a3) 

i można zawsze używaj Numby i zwiększaj prędkość znacznie bez wektoryzacji (i bez użycia większej ilości pamięci).

Answer 3

Różnice w wydajności między numpy i matlab zawsze mnie denerwowały. Często w końcu sprowadzają się do leżących pod nimi bibliotek. O ile mi wiadomo, Matlab używa domyślnie pełnego atlasu, podczas gdy numpy używa światła na klapie. Matlab uważa, że ​​ludzie nie dbają o przestrzeń i masę, podczas gdy numpy uważa, że ​​ludzie to robią. Similar question z dobrą odpowiedzią.