Multiply hohe Ordnung Matrizen mit numpy

Question

Ich habe dieses Spielzeug Problem, dass mein viel größeres Problem widerspiegelt:

import numpy as np 

ind = np.ones((3,2,4)) # shape=(3L, 2L, 4L) 
dist = np.array([[0.1,0.3],[1,2],[0,1]]) # shape=(3L, 2L) 

ans = np.array([np.dot(dist[i],ind[i]) for i in xrange(dist.shape[0])]) # shape=(3L, 4L) 
print ans 

""" prints: 
    [[ 0.4 0.4 0.4 0.4] 
    [ 3. 3. 3. 3. ] 
    [ 1. 1. 1. 1. ]] 
""" 

ich es so schnell wie möglich tun will, so numpy Funktionen mit ans berechnen sollte der beste Ansatz sein , da diese Operation schwer ist und meine Matrizen ziemlich groß sind.

Ich sah this post, aber die Formen sind unterschiedlich und ich kann nicht verstehen, welche axes ich für dieses Problem verwenden sollte. Ich bin mir aber sicher, dass tensordot die Antwort haben sollte. Irgendwelche Vorschläge?

EDIT: Ich nahm @ajcr's answer, aber bitte auch meine eigene Antwort lesen, kann es andere helfen ...

Answer 1

Sie np.einsum nutzen, um die Operation zu tun, da es für eine sehr sorgfältige Kontrolle ermöglicht über die multiplizierten Achsen und die summiert werden:

>>> np.einsum('ijk,ij->ik', ind, dist) 
array([[ 0.4, 0.4, 0.4, 0.4], 
     [ 3. , 3. , 3. , 3. ], 
     [ 1. , 1. , 1. , 1. ]]) 

die Funktion multipliziert die Einträge in der ersten Achse ind mit den Einträgen in der ersten Achse von dist (Subskript 'i'). Dito für die zweite Achse jedes Arrays (Index 'j'). Anstatt ein 3D-Array zurückzugeben, teilen wir einsum mit, die Werte entlang der Achse 'j' zu summieren, indem wir es aus den Ausgabeskripten auslassen, wodurch ein 2D-Array zurückgegeben wird.

---

np.tensordot ist schwieriger für dieses Problem anzuwenden. Es summiert automatisch die Produkte von Achsen. Wir wollen aber zwei Sätze von Produkten, sondern nur ein von ihnen zu summieren.

Schreiben np.tensordot(ind, dist, axes=[1, 1]) (wie in der Antwort, die Sie verknüpft) berechnet die richtigen Werte für Sie, aber gibt ein 3D-Array mit Form (3, 4, 3) zurück. Wenn Sie die Speicherkosten eines größeren Arrays leisten können, könnten Sie verwenden:

np.tensordot(ind, dist, axes=[1, 1])[0].T 

Dies gibt Ihnen das richtige Ergebnis, sondern weil tensordot zuerst eine viel größere als erforderlich Array schafft, scheint einsum besser zu sein Möglichkeit.

Answer 2

Nach @ajcr's great answer wollte ich, um zu bestimmen, welche Methode die schnellste ist, so habe ich timeit:

import timeit 

setup_code = """ 
import numpy as np 
i,j,k = (300,200,400) 
ind = np.ones((i,j,k)) #shape=(3L, 2L, 4L) 
dist = np.random.rand(i,j) #shape=(3L, 2L) 
""" 

basic ="np.array([np.dot(dist[l],ind[l]) for l in xrange(dist.shape[0])])" 
einsum = "np.einsum('ijk,ij->ik', ind, dist)" 
tensor= "np.tensordot(ind, dist, axes=[1, 1])[0].T" 

print "tensor - total time:", min(timeit.repeat(stmt=tensor,setup=setup_code,number=10,repeat=3)) 
print "basic - total time:", min(timeit.repeat(stmt=basic,setup=setup_code,number=10,repeat=3)) 
print "einsum - total time:", min(timeit.repeat(stmt=einsum,setup=setup_code,number=10,repeat=3)) 

Die überraschenden Ergebnisse waren:

tensor - total time: 6.59519493952 
basic - total time: 0.159871203461 
einsum - total time: 0.263569731028 

So offensichtlich tensordot wurde unter Verwendung der falsche Weg zu Mach es (ganz zu schweigen von memory error in größeren Beispielen, wie @ajcr angegeben).

Da dieses Beispiel klein war, änderte ich die Matrizen Größe i,j,k = (3000,200,400), nur die Reihenfolge umgedreht zu sein, sicherzustellen, dass es einen weiteren Test mit einer höheren Anzahl von Wiederholungen nicht bewirken und eingerichtet hat:

print "einsum - total time:", min(timeit.repeat(stmt=einsum,setup=setup_code,number=50,repeat=3)) 
print "basic - total time:", min(timeit.repeat(stmt=basic,setup=setup_code,number=50,repeat=3)) 

Die Ergebnisse im Einklang mit dem ersten Lauf waren:

einsum - total time: 13.3184077671 
basic - total time: 8.44810031351 

, eine andere Art von Größenwachstum Prüfung jedoch - i,j,k = (30000,20,40) folgenden Ergebnissen geführt:

einsum - total time: 0.325594117768 
basic - total time: 0.926416766397 

Siehe die Kommentare nach Erklärungen für diese Ergebnisse.

Die Moral ist, wenn für die schnellste Lösung für ein bestimmtes Problem suchen, versucht, Daten zu erzeugen, die ähnlich wie bei den Originaldaten wie möglich, in der Bezeichnung der Arten und Formen sind. In meinem Fall i ist viel kleiner als j,k und so blieb ich mit der hässlichen Version, die in diesem Fall auch die schnellste ist.