Perché l'einsum di Numpy è più lento delle funzioni incorporate di numpy?

Question

Di solito ho ottenuto una buona performance dalla funzione di einsum di Numpy (e mi piace che sia sintassi). @ La risposta di Ophion a this question mostra che - per i casi testati - einsum supera costantemente le funzioni "built-in" (a volte un po ', a volte molto spesso). Ma ho appena incontrato un caso in cui l'einsum è molto più lento. Considerate le seguenti funzioni equivalenti:

(M, K) = (1000000, 20) 
C = np.random.rand(K, K) 
X = np.random.rand(M, K) 

def func_dot(C, X): 
    Y = X.dot(C) 
    return np.sum(Y * X, axis=1) 

def func_einsum(C, X): 
    return np.einsum('ik,km,im->i', X, C, X) 

def func_einsum2(C, X): 
    # Like func_einsum but break it into two steps. 
    A = np.einsum('ik,km', X, C) 
    return np.einsum('ik,ik->i', A, X) 

mi aspettavo func_einsum a correre più veloce, ma non è questo che ho incontrato. Eseguito su una CPU quad-core con hyperthreading, versione NumPy 1.9.0.dev-7ae0206, e multithreading con OpenBLAS, ottengo i seguenti risultati:

In [2]: %time y1 = func_dot(C, X) 
CPU times: user 320 ms, sys: 312 ms, total: 632 ms 
Wall time: 209 ms 
In [3]: %time y2 = func_einsum(C, X) 
CPU times: user 844 ms, sys: 0 ns, total: 844 ms 
Wall time: 842 ms 
In [4]: %time y3 = func_einsum2(C, X) 
CPU times: user 292 ms, sys: 44 ms, total: 336 ms 
Wall time: 334 ms 

Quando ho aumentare K a 200, le differenze sono più estreme :

In [2]: %time y1= func_dot(C, X) 
CPU times: user 4.5 s, sys: 1.02 s, total: 5.52 s 
Wall time: 2.3 s 
In [3]: %time y2= func_einsum(C, X) 
CPU times: user 1min 16s, sys: 44 ms, total: 1min 16s 
Wall time: 1min 16s 
In [4]: %time y3 = func_einsum2(C, X) 
CPU times: user 15.3 s, sys: 312 ms, total: 15.6 s 
Wall time: 15.6 s 

Qualcuno può spiegare perché l'einsum è molto più lento qui?

Se è importante, qui è il mio config NumPy:

In [6]: np.show_config() 
lapack_info: 
    libraries = ['openblas'] 
    library_dirs = ['/usr/local/lib'] 
    language = f77 
atlas_threads_info: 
    libraries = ['openblas'] 
    library_dirs = ['/usr/local/lib'] 
    define_macros = [('ATLAS_WITHOUT_LAPACK', None)] 
    language = c 
    include_dirs = ['/usr/local/include'] 
blas_opt_info: 
    libraries = ['openblas'] 
    library_dirs = ['/usr/local/lib'] 
    define_macros = [('ATLAS_INFO', '"\\"None\\""')] 
    language = c 
    include_dirs = ['/usr/local/include'] 
atlas_blas_threads_info: 
    libraries = ['openblas'] 
    library_dirs = ['/usr/local/lib'] 
    define_macros = [('ATLAS_INFO', '"\\"None\\""')] 
    language = c 
    include_dirs = ['/usr/local/include'] 
lapack_opt_info: 
    libraries = ['openblas', 'openblas'] 
    library_dirs = ['/usr/local/lib'] 
    define_macros = [('ATLAS_WITHOUT_LAPACK', None)] 
    language = f77 
    include_dirs = ['/usr/local/include'] 
lapack_mkl_info: 
    NOT AVAILABLE 
blas_mkl_info: 
    NOT AVAILABLE 
mkl_info: 
    NOT AVAILABLE 

Answer 1

si può avere il meglio dei due mondi:

def func_dot_einsum(C, X): 
    Y = X.dot(C) 
    return np.einsum('ij,ij->i', Y, X) 

Sul mio sistema:

In [7]: %timeit func_dot(C, X) 
10 loops, best of 3: 31.1 ms per loop 

In [8]: %timeit func_einsum(C, X) 
10 loops, best of 3: 105 ms per loop 

In [9]: %timeit func_einsum2(C, X) 
10 loops, best of 3: 43.5 ms per loop 

In [10]: %timeit func_dot_einsum(C, X) 
10 loops, best of 3: 21 ms per loop 

Quando disponibili , np.dot utilizza BLAS, MKL o qualsiasi altra libreria tu abbia. Quindi la chiamata a np.dot è quasi certamente in fase di multithreading. np.einsum ha i propri loop, quindi non usa nessuna di queste ottimizzazioni, a parte il proprio uso di SIMD per velocizzare le cose su un'implementazione di vaniglia C.

---

Poi c'è la chiamata einsum multi-input che viene eseguito molto più lento ... La fonte NumPy per einsum è molto complesso e non capiscono fino in fondo di esso. Quindi, essere informati che la segue è speculativa nella migliore delle ipotesi, ma ecco quello che penso sta succedendo ...

Quando si esegue qualcosa come np.einsum('ij,ij->i', a, b), del vantaggio nel fare np.sum(a*b, axis=1) viene da evitare di dover istanziare la matrice intermedia con tutte le prodotti e due volte su di esso. Quindi a livello basso ciò che accade è qualcosa di simile:

for i in range(I): 
    out[i] = 0 
    for j in range(J): 
     out[i] += a[i, j] * b[i, j] 

dire ora che siete alla ricerca di qualcosa di simile:

np.einsum('ij,jk,ik->i', a, b, c) 

Si potrebbe fare la stessa operazione

np.sum(a[:, :, None] * b[None, :, :] * c[:, None, :], axis=(1, 2)) 

E quello che penso einsum è di eseguire questo ultimo codice senza dover istanziare l'enorme array intermedio, che certamente rende le cose più veloci:

In [29]: a, b, c = np.random.rand(3, 100, 100) 

In [30]: %timeit np.einsum('ij,jk,ik->i', a, b, c) 
100 loops, best of 3: 2.41 ms per loop 

In [31]: %timeit np.sum(a[:, :, None] * b[None, :, :] * c[:, None, :], axis=(1, 2)) 
100 loops, best of 3: 12.3 ms per loop 

Ma se si guarda attentamente, liberarsi della memoria intermedia può essere una cosa terribile.Questo è quello che penso einsum sta facendo a livello basso:

for i in range(I): 
    out[i] = 0 
    for j in range(J): 
     for k in range(K): 
      out[i] += a[i, j] * b[j, k] * c[i, k] 

ma si stanno ripetendo una tonnellata di operazioni! Se invece fatto:

for i in range(I): 
    out[i] = 0 
    for j in range(J): 
     temp = 0 
     for k in range(K): 
      temp += b[j, k] * c[i, k] 
     out[i] += a[i, j] * temp 

si sarebbe fare I * J * (K-1) meno moltiplicazioni (e I * J aggiunte extra), e risparmiare un sacco di tempo. La mia ipotesi è che einsum non sia abbastanza intelligente da ottimizzare le cose a questo livello. Nel source code c'è un suggerimento che ottimizza solo operazioni di 1 o 2 operandi, non 3. In ogni caso automatizzare questo per ingressi generali sembra tutt'altro che semplice ...

Answer 2

einsum ha un caso specializzato per '2 operandi , ndim = 2 '. In questo caso ci sono 3 operandi e un totale di 3 dimensioni. Quindi deve usare un generale nditer.

Durante il tentativo di capire come l'ingresso stringa viene analizzata, ho scritto un simulatore puro Python einsum, https://github.com/hpaulj/numpy-einsum/blob/master/einsum_py.py

L'(ridotta) einsum e somma di prodotti funzioni sono:

def myeinsum(subscripts, *ops, **kwargs): 
    # dropin preplacement for np.einsum (more or less) 
    <parse subscript strings> 
    <prepare op_axes> 
    x = sum_of_prod(ops, op_axes, **kwargs) 
    return x 

def sum_of_prod(ops, op_axes,...): 
    ... 
    it = np.nditer(ops, flags, op_flags, op_axes) 
    it.operands[nop][...] = 0 
    it.reset() 
    for (x,y,z,w) in it: 
     w[...] += x*y*z 
    return it.operands[nop] 

Debug output per myeinsum('ik,km,im->i',X,C,X,debug=True) con (M,K)=(10,5)

{'max_label': 109, 
'min_label': 105, 
'nop': 3, 
'shapes': [(10, 5), (5, 5), (10, 5)], 
....}} 
... 
iter labels: [105, 107, 109],'ikm' 

op_axes [[0, 1, -1], [-1, 0, 1], [0, -1, 1], [0, -1, -1]] 

Se si scrive una funzione sum-of-prod come questo mi n cython si dovrebbe ottenere qualcosa di simile al einsum generalizzato.

Con l'intero (M,K), questo einsum simulato è 6-7 volte più lento.

---

qualche edificio tempi sulle altre risposte:

In [84]: timeit np.dot(X,C) 
1 loops, best of 3: 781 ms per loop 

In [85]: timeit np.einsum('ik,km->im',X,C) 
1 loops, best of 3: 1.28 s per loop 

In [86]: timeit np.einsum('im,im->i',A,X) 
10 loops, best of 3: 163 ms per loop 

Questo 'im,im->i' step is substantially faster than the other. The sum dimension, m is only 20. I suspect einsum` sta trattando questo come un caso speciale.

In [87]: timeit np.einsum('im,im->i',np.dot(X,C),X) 
1 loops, best of 3: 950 ms per loop 

In [88]: timeit np.einsum('im,im->i',np.einsum('ik,km->im',X,C),X) 
1 loops, best of 3: 1.45 s per loop 

I tempi per questi calcoli compositi sono semplicemente le somme dei pezzi corrispondenti.