Schnelle Numpy Loops

Question

Wie optimieren Sie diesen Code (ohne Vektorisierung, da dies die Semantik der Berechnung mit bis führt, die oft weit entfernt, nicht-triviale):

slow_lib.py: 
import numpy as np 

def foo(): 
    size = 200 
    np.random.seed(1000031212) 
    bar = np.random.rand(size, size) 
    moo = np.zeros((size,size), dtype = np.float) 
    for i in range(0,size): 
     for j in range(0,size): 
      val = bar[j] 
      moo += np.outer(val, val) 

Der Punkt ist, dass solche Art Schleifen ziemlich häufig Operationen entsprechen, wo Sie Doppelsummen über einige Vektoroperation haben. Diese

ist ziemlich langsam:

>>t = timeit.timeit('foo()', 'from slow_lib import foo', number = 10) 
>>print ("took: "+str(t)) 
took: 41.165681839 

Ok, also dann ist es lassen cynothize und Typ Anmerkungen hinzufügen mag es kein Morgen:

c_slow_lib.pyx: 
import numpy as np 
cimport numpy as np 
import cython 
@cython.boundscheck(False) 
@cython.wraparound(False) 

def foo(): 
    cdef int size = 200 
    cdef int i,j 
    np.random.seed(1000031212) 
    cdef np.ndarray[np.double_t, ndim=2] bar = np.random.rand(size, size) 
    cdef np.ndarray[np.double_t, ndim=2] moo = np.zeros((size,size), dtype = np.float) 
    cdef np.ndarray[np.double_t, ndim=1] val 
    for i in xrange(0,size): 
     for j in xrange(0,size): 
      val = bar[j] 
      moo += np.outer(val, val) 


>>t = timeit.timeit('foo()', 'from c_slow_lib import foo', number = 10) 
>>print ("took: "+str(t)) 
took: 42.3104710579 

... ehr ... was? Numba zur Rettung!

numba_slow_lib.py: 
import numpy as np 
from numba import jit 

size = 200 
np.random.seed(1000031212) 

bar = np.random.rand(size, size) 

@jit 
def foo(): 
    bar = np.random.rand(size, size) 
    moo = np.zeros((size,size), dtype = np.float) 
    for i in range(0,size): 
     for j in range(0,size): 
      val = bar[j] 
      moo += np.outer(val, val) 

>>t = timeit.timeit('foo()', 'from numba_slow_lib import foo', number = 10) 
>>print("took: "+str(t)) 
took: 40.7327859402 

Also gibt es wirklich keine Möglichkeit, dies zu beschleunigen? Der Punkt ist:

• wenn ich die innere Schleife in eine vektorisierte Version (Gebäude eine größere Matrix, die die innere Schleife und dann ruft np.outer auf der größeren Matrix) umwandeln ich viel schnelleren Code.
• Wenn ich etwas Ähnliches in Matlab (R2016a) implementiere, wird dies durch JIT recht gut.

Answer 1

Hier ist der Code für outer:

def outer(a, b, out=None):  
    a = asarray(a) 
    b = asarray(b) 
    return multiply(a.ravel()[:, newaxis], b.ravel()[newaxis,:], out) 

Also jeder Aufruf outer beinhaltet eine Reihe von Python-Anrufe. Diese rufen schließlich den kompilierten Code auf, um die Multiplikation durchzuführen. Aber jeder hat einen Overhead, der nichts mit der Größe Ihrer Arrays zu tun hat.

Also 200 (200 ** 2?) Anrufe zu outer haben all diesen Aufwand, während ein Anruf an outer mit allen 200 Zeilen hat einen Overhead-Satz, gefolgt von einem schnellen kompilierten Vorgang.

cython und numba nicht kompilieren oder sonst umgehen den Python-Code in outer. Alles, was sie tun können, ist die Straffung des von Ihnen geschriebenen Iterationscodes - und das kostet nicht viel Zeit.

Ohne ins Detail zu gehen, muss der MATLAB jit die 'äußere' mit schnellerem Code ersetzen können - er schreibt die Iteration neu. Aber meine Erfahrung mit MATLAB stammt aus einer Zeit vor ihrem JIT.

Für echte Geschwindigkeit Verbesserungen mit cython und numba müssen Sie primitiven numpy/Python-Code den ganzen Weg hinunter verwenden. Oder besser konzentrieren Sie sich auf langsame innere Teile.

Ersetzen Sie Ihre outer mit eine vereinfachte Version Zeit etwa in der Hälfte laufen Cuts:

def foo1(N): 
     size = N 
     np.random.seed(1000031212) 
     bar = np.random.rand(size, size) 
     moo = np.zeros((size,size), dtype = np.float) 
     for i in range(0,size): 
       for j in range(0,size): 
         val = bar[j] 
         moo += val[:,None]*val 
     return moo 

Mit der vollen N=200 Ihre Funktion nahm 17s pro Schleife. Wenn ich die inneren zwei Zeilen durch pass (keine Berechnung) ersetze, fällt die Zeit auf 3ms pro Schleife. Mit anderen Worten, der äußere Schleifenmechanismus ist kein großer Zeitverbraucher, zumindest nicht verglichen mit vielen Aufrufen an outer().

Answer 2

Speicher erlaubt, können Sie np.einsum verwenden, um diese schweren Berechnungen in einer vektorisierten Weise durchzuführen, wie so -

moo = size*np.einsum('ij,ik->jk',bar,bar) 

One auch np.tensordot verwenden -

moo = size*np.tensordot(bar,bar,axes=(0,0)) 

Oder einfach np.dot -

moo = size*bar.T.dot(bar) 

Answer 3

Viele Tutorials und Demonstrationen von Cython, Numba, etc. lassen es so aussehen, als könnten diese Tools Ihren Code automatisch beschleunigen, aber in der Praxis ist dies oft nicht der Fall: Sie müssen Ihren Code ein wenig modifizieren um die beste Leistung zu erzielen. Wenn Sie bereits ein gewisses Maß an Vektorisierung implementiert haben, bedeutet das normalerweise, ALLE Schleifen zu schreiben. Gründe, die Numpy-Array-Operationen nicht optimal sind:

• Viele temporäre Arrays werden erstellt und durchlaufen;
• Signifikanter Overhead pro Anruf, wenn die Arrays klein sind;
• Kurzschlusslogik kann nicht implementiert werden, da Arrays als Ganzes verarbeitet werden;
• Manchmal kann der optimale Algorithmus nicht mit Array-Ausdrücken ausgedrückt werden und Sie entscheiden sich für einen Algorithmus mit einer schlechteren Zeitkomplexität.

Mit Numba oder Cython werden diese Probleme nicht optimiert! Stattdessen ermöglichen diese Tools das Schreiben von Code, der viel schneller ist als der einfache Python.

Auch für Numba sollten Sie den Unterschied zwischen "object mode" and "nopython mode" beachten. Die engen Schleifen aus Ihrem Beispiel müssen im nopython-Modus ausgeführt werden, um eine signifikante Beschleunigung zu erreichen. numpy.outer ist jedoch not yet supported by Numba, was dazu führt, dass die Funktion im Objektmodus kompiliert wird. Dekorieren Sie mit jit(nopython=True), um solche Fälle eine Ausnahme auslösen zu lassen.

Beispiel eine Beschleunigung zu zeigen, ist in der Tat möglich:

import numpy as np 
from numba import jit 

@jit 
def foo_nb(bar): 
    size = bar.shape[0] 
    moo = np.zeros((size, size)) 
    for i in range(0,size): 
     for j in range(0,size): 
      val = bar[j] 
      moo += np.outer(val, val) 
    return moo 

@jit 
def foo_nb2(bar): 
    size = bar.shape[0] 
    moo = np.zeros((size, size)) 
    for i in range(size): 
     for j in range(size): 
      for k in range(0,size): 
       for l in range(0,size): 
        moo[k,l] += bar[j,k] * bar[j,l] 
    return moo 

size = 100 
bar = np.random.rand(size, size) 

np.allclose(foo_nb(bar), foo_nb2(bar)) 
# True 

%timeit foo_nb(bar) 
# 1 loop, best of 3: 816 ms per loop 
%timeit foo_nb2(bar) 
# 10 loops, best of 3: 176 ms per loop 

Answer 4

Das Beispiel, das Sie uns zeigen, Art ineffizienter Algorithmus ist, da Sie die gleichen äußeren Produkt mehrfach berechnen. Die resultierende Zeitkomplexität ist O (n^4). Es kann auf n^3 reduziert werden.

for i in range(0,size): 
    val = bar[i] 
    moo += size * np.outer(val, val)