Warum ist B = numpy.dot (A, x) so viel langsamer durch B [i,:,:] = numpy.dot (A [i,:,:], x))?

Question

Ich bekomme einige Effizienz Testergebnisse, die ich nicht erklären kann.

Ich möchte eine Matrix B zusammenstellen, deren i-te Einträge B [i,:,:] = A [i,:,:] Punkt (x), wo jedes A [i,:,:] ist eine 2D-Matrix, und so ist x.

Ich kann dies drei Möglichkeiten, um die Leistung zu testen Ich mache zufällige (numpy.random.randn) Matrizen A = (10,1000,1000), x = (1000,1200). Ich erhalte die folgenden Zeitergebnisse:

(1) Einzel mehrdimensionalen Skalarprodukt

B = A.dot(x) 

total time: 102.361 s 

(2) bis i Looping und Durchführung 2D-Punktprodukte

# initialize B = np.zeros([dim1, dim2, dim3]) 
    for i in range(A.shape[0]): 
     B[i,:,:] = A[i,:,:].dot(x) 

total time: 0.826 s 

(3) numpy.einsum

So ist Option (2) bei weitem die schnellste. Betrachtet man nur (1) und (2), sehe ich nicht den großen Unterschied zwischen ihnen. Wie kann das Durchschleifen und Ausführen von 2D-Punktprodukten ~ 124 mal schneller sein? Beide benutzen numpy.dot. Irgendwelche Einsichten?

I enthalten den Code für die obigen Ergebnisse verwendet knapp unter:

import numpy as np 
import numpy.random as npr 
import time 

dim1, dim2, dim3 = 10, 1000, 1200 
A = npr.randn(dim1, dim2, dim2) 
x = npr.randn(dim2, dim3) 

# consider three ways of assembling the same matrix B: B1, B2, B3 

t = time.time() 
B1 = np.dot(A,x) 
td1 = time.time() - t 
print "a single dot product of A [shape = (%d, %d, %d)] with x [shape = (%d, %d)] completes in %.3f s" \ 
    % (A.shape[0], A.shape[1], A.shape[2], x.shape[0], x.shape[1], td1) 


B2 = np.zeros([A.shape[0], x.shape[0], x.shape[1]]) 
t = time.time() 
for i in range(A.shape[0]): 
    B2[i,:,:] = np.dot(A[i,:,:], x) 
td2 = time.time() - t 
print "taking %d dot products of 2D dot products A[i,:,:] [shape = (%d, %d)] with x [shape = (%d, %d)] completes in %.3f s" \ 
    % (A.shape[0], A.shape[1], A.shape[2], x.shape[0], x.shape[1], td2) 

t = time.time() 
B3 = np.einsum("ijk, kl -> ijl", A, x) 
td3 = time.time() - t 
print "using np.einsum, it completes in %.3f s" % td3 

Answer 1

Bei kleineren dims 10,100,200, erhalte ich eine ähnliche Ranking

In [355]: %%timeit 
    .....: B=np.zeros((N,M,L)) 
    .....: for i in range(N): 
       B[i,:,:]=np.dot(A[i,:,:],x) 
    .....: 
10 loops, best of 3: 22.5 ms per loop 
In [356]: timeit np.dot(A,x) 
10 loops, best of 3: 44.2 ms per loop 
In [357]: timeit np.einsum('ijk,km->ijm',A,x) 
10 loops, best of 3: 29 ms per loop 

In [367]: timeit np.dot(A.reshape(-1,M),x).reshape(N,M,L) 
10 loops, best of 3: 22.1 ms per loop 

In [375]: timeit np.tensordot(A,x,(2,0)) 
10 loops, best of 3: 22.2 ms per loop 

die itererative schneller ist, wenn auch nicht so stark wie in Ihrem Fall.

Dies ist wahrscheinlich wahr, solange diese iterierende Dimension im Vergleich zu den anderen klein ist. In diesem Fall ist der Overhead der Iteration (Funktionsaufrufe usw.) im Vergleich zur Berechnungszeit gering. Und alle Werte auf einmal zu verwenden, verbraucht mehr Speicher.

Ich versuchte eine dot Variante, wo ich A in 2d umgeformt, denken, dass dot macht diese Art der Umgestaltung intern. Ich bin überrascht, dass es tatsächlich am schnellsten ist. tensordot macht wahrscheinlich die gleiche Umformung (dieser Code, wenn Python lesbar).

---

einsum richtet eine 'Produktsumme' Iteration Einbeziehung 4 Variablen, die i,j,k,m - das ist dim1*dim2*dim2*dim3 Schritte mit der C-Ebene nditer. Je mehr Indizes Sie haben, desto größer ist der Iterationsraum.

Answer 2

Ich bin nicht allzu vertraut mit numpy der C-API und die numpy.dot ist eine solche eingebaute Funktion, die unter _dotblas in früher früher Versionen.

Trotzdem sind hier meine Gedanken.

1)numpy.dot verwendet unterschiedliche Pfade für 2-dimensionale Arrays und n-dimensionale Arrays. Aus den numpy.dot ‚s online documentation:

für 2-D-Arrays es äquivalent zu Matrizenmultiplikation ist, und für 1-D-Arrays Skalarprodukt von Vektoren (ohne komplexe Konjugation). Für N Dimensionen ist es ein Summenprodukt über der letzten Achse von a und dem vorletzten von b

Punkt (a, b) [i, j, k, m] = Summe (a [i, j ,:] * b [k,:, m])

Also für 2-D-Arrays Sie immer garantiert werden, einen Anruf zu BLAS der dgemm jedoch für ND-Arrays haben numpy die Multiplikation wählen könnten Achsen für Arrays, die entspricht vielleicht nicht der schnellsten sich ändernden Achse (wie Sie aus dem von mir geposteten Auszug sehen können), und als Ergebnis könnte die volle Leistung von dgemm verpasst werden.

2) Ihr Array A ist zu groß, um in den CPU-Cache geladen zu werden. In Ihrem Beispiel verwenden Sie A mit Abmessungen (10,1000,1000) die

In [1]: A.nbytes 
80000000 
In [2]: 80000000/1024 
78125 

dass fast 80MB ist gibt, viel größer als Ihre Cache-Größe. Also wieder verlieren Sie die meisten dgemm's Macht genau dort.

3) Auch die Funktionen sind etwas ungenau. Die time-Funktion in Python ist bekanntermaßen nicht präzise. Verwenden Sie stattdessen timeit.

Also alle, die die oben genannten Punkte im Auge, wollen wir versuchen, mit Arrays zu experimentieren, die auf den Cache

dim1, dim2, dim3 = 20, 20, 20 
A = np.random.rand(dim1, dim2, dim2) 
x = np.random.rand(dim2, dim3) 

def for_dot1(A,x): 
    for i in range(A.shape[0]): 
     np.dot(A[i,:,:], x) 

def for_dot2(A,x): 
    for i in range(A.shape[0]): 
     np.dot(A[:,i,:], x)  

def for_dot3(A,x): 
    for i in range(A.shape[0]): 
     np.dot(A[:,:,i], x) 

und hier sind die Zeiten, die ich (mit numpy 1.9.2 gegen OpenBLAS 0.2.14 gebaut) erhalten geladen werden kann:

In [3]: %timeit np.dot(A,x) 
10000 loops, best of 3: 174 µs per loop 
In [4]: %timeit np.einsum("ijk, kl -> ijl", A, x) 
10000 loops, best of 3: 108 µs per loop 
In [5]: %timeit np.einsum("ijk, lk -> ijl", A, x) 
10000 loops, best of 3: 97.1 µs per loop 
In [6]: %timeit np.einsum("ikj, kl -> ijl", A, x) 
1000 loops, best of 3: 238 µs per loop 
In [7]: %timeit np.einsum("kij, kl -> ijl", A, x) 
10000 loops, best of 3: 113 µs per loop 
In [8]: %timeit for_dot1(A,x) 
10000 loops, best of 3: 101 µs per loop 
In [9]: %timeit for_dot2(A,x) 
10000 loops, best of 3: 131 µs per loop 
In [10]: %timeit for_dot3(A,x) 
10000 loops, best of 3: 133 µs per loop 

Beachten Sie, dass es immer noch einen Zeitunterschied gibt, aber nicht in Größenordnungen. Beachten Sie auch die Wichtigkeit von choosing the axis of multiplication. Nun kann vielleicht ein nüchterner Entwickler etwas Licht auf das werfen, was numpy.dot tatsächlich unter der Haube für N-D-Arrays tut.