veloce Numpy Loops

Question

Come si fa a ottimizzare questo codice (senza vettorializzazione, in quanto ciò conduce ad usare la semantica del calcolo, che è spesso lontano da essere non banale):

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) 

Il punto è che tali loop di tipo corrispondono molto spesso alle operazioni in cui si hanno somme doppie rispetto ad un'operazione vettoriale.

Questo è piuttosto lento:

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

Ok, quindi cerchiamo di cynothize e aggiungere annotazioni di tipo ama non c'è domani:

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 ... che cosa? Numba in soccorso!

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 

Quindi non c'è davvero alcun modo per accelerare questo? Il punto è:

• Se converto il ciclo interno in un vettorizzati versione (la costruzione di una matrice più grande che rappresenta il ciclo interno e quindi chiamando np.outer sulla matrice più grande) ottengo molto codice più veloce.
• se implemento qualcosa di simile in Matlab (R2016a) questo si comporta abbastanza bene a causa di JIT.

Answer 1

Ecco il codice per outer:

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

Così ogni chiamata a outer comporta una serie di chiamate di pitone. Coloro che alla fine chiamano codice compilato per eseguire la moltiplicazione. Ma ognuno di essi comporta un sovraccarico che non ha nulla a che fare con la dimensione dei tuoi array.

Quindi 200 (200 ** 2?) Chiamate a outer avranno tutto questo sovraccarico, mentre una chiamata a outer con tutte le 200 righe ha una serie overhead, seguita da un'operazione rapida compilata.

cython e numba non compila o altrimenti ignora il codice Python in outer. Tutto quello che possono fare è semplificare il codice di iterazione che hai scritto - e questo non sta consumando molto tempo.

Senza entrare nei dettagli, il jolly MATLAB deve essere in grado di sostituire "esterno" con codice più veloce - riscrive l'iterazione. Ma la mia esperienza con MATLAB risale a un tempo prima della sua esca.

Per miglioramenti della velocità reale con cython e numba è necessario utilizzare il codice di Numpy/Python primitivo fino in fondo. O meglio ancora concentrare i tuoi sforzi su pezzi interiori lenti.

Sostituendo il vostro outer con una versione semplificata tagli eseguito tempo circa a metà:

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 

Con il pieno N=200 vostra funzione ha preso 17s per loop. Se sostituisco le due linee interne con pass (nessun calcolo), il tempo scende a 3 ms per loop. In altre parole, il meccanismo di loop esterno non è un grande consumatore, almeno non rispetto a molte chiamate a outer().

Answer 2

memoria permettendo, è possibile utilizzare np.einsum per eseguire questi calcoli pesanti in modo vettorizzati, in questo modo -

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

Si può anche utilizzare np.tensordot -

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

O semplicemente np.dot -

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

Answer 3

Molte esercitazioni e dimostrazioni di Cython, Numba, ecc. Fanno sembrare che questi strumenti possano velocizzare il codice in modo automatico, ma in pratica spesso non è così: dovrai modificare leggermente il tuo codice per estrarre la migliore prestazione. Se hai già implementato un certo grado di vettorizzazione, di solito significa scrivere TUTTI i loop. Ragioni Le operazioni di Numpy array non sono ottimali:

• Un sacco di matrici temporanee vengono create e ripetute;
• Overhead significativo per chiamata se gli array sono piccoli;
• La logica di cortocircuito non può essere implementata, poiché gli array vengono elaborati nel loro complesso;
• A volte l'algoritmo ottimale non può essere espresso utilizzando le espressioni di matrice e ci si accontenta di un algoritmo con una complessità temporale peggiore.

L'utilizzo di Numba o Cython non ottimizza questi problemi! Invece, questi strumenti ti permettono di scrivere codice loopy molto più veloce del semplice Python.

Inoltre, per Numba in particolare, è necessario essere consapevoli della differenza tra "object mode" and "nopython mode". I loop stretti del tuo esempio devono essere eseguiti in modalità nopython per fornire un significativo aumento della velocità. Tuttavia, numpy.outer è not yet supported by Numba, con conseguente funzione di compilazione in modalità oggetto. Decorare con jit(nopython=True) per consentire a tali casi di generare un'eccezione.

esempio per dimostrare un aumento di velocità è infatti possibile:

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

L'esempio ci mostri è una specie di algoritmo inefficiente, in quanto si calcola lo stesso prodotto esterno più volte. La complessità temporale risultante è O (n^4). Può essere ridotto a n^3.

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