Python - Numpy: 3D-Matrix * 2D-Vektor schnelle Berechnung

Question

Hallo an alle,

Hier ist, was ich tun möchte. Ich habe zwei Arrays bekam:

• rotation_matrices die 50 zweidimensionale Drehmatrizen enthält. Jede Rotationsmatrix ist als (2,2) geformt. So wird rotation_matrices als (2,2,50) geformt.
• Vektoren, die 50 zweidimensionale Vektoren enthält. Somit ist es geformt als (2,50).

Ich möchte (falls vorhanden) eine einzeilige numpy Operation, die mir die (2,50) Array gibt, die die gedrehten Vektoren enthält, nennen wir es rotated_vectors. Was ich meine ist, dass das k-ith-Element von rotierte_vektoren das Produkt der k-ith-Rotationsmatrix mit dem k-i-ten Vektor enthält.

Im Moment habe ich mit der Schleife kommen, der folgt:

for ind,elt in enumerate(np.arange(nb_of_vectors)): 
      rotated_vector[ind] = np.dot(rotation_matrices[:,:,ind], vectors[:,ind]) 

ich Verbesserungs denke, es gibt Raum. Wenn Sie einen Vorschlag haben, sind Sie willkommen.

Vielen Dank für Ihre Zeit.

Jagaral

Answer 1

Für solche Reduzierungen, die Ausrichtung erfordern entlang einer oder mehrerer Achsen, kann man verwenden np.einsum -

rotated_vector = np.einsum('ijk,jk->ki',rotation_matrices,vectors) 

Bitte beachten Sie, dass die Ausgabe der Form sein würde (N,2), wo N ist die Anzahl der Vektoren. Wenn Sie statt dessen nach einer Ausgabe der Form (2,N) suchen und den ursprünglichen Code stattdessen wie folgt benötigen: rotated_vector[:,ind] = np.dot(...), bearbeiten Sie einfach die Ausgabezeichenfolgenschreibweise zu ik anstelle von ki.

Runtime Test -

In [24]: def org_app(rotation_matrices,vectors): 
    ...:  nb_of_vectors = vectors.shape[1] 
    ...:  r = np.zeros((nb_of_vectors,2)) 
    ...:  for ind,elt in enumerate(np.arange(nb_of_vectors)): 
    ...:   r[ind] = np.dot(rotation_matrices[:,:,ind], vectors[:,ind]) 
    ...:  return r 
    ...: 

In [25]: # Input arrays 
    ...: rotation_matrices = np.random.rand(2,2,50) 
    ...: vectors = np.random.rand(2,50) 
    ...: 

In [26]: out1 = org_app(rotation_matrices,vectors) 

In [27]: out2 = np.einsum('ijk,jk->ki',rotation_matrices,vectors) 

In [28]: np.allclose(out1,out2) # Verify results 
Out[28]: True 

In [29]: %timeit org_app(rotation_matrices,vectors) 
10000 loops, best of 3: 196 µs per loop 

In [30]: %timeit np.einsum('ijk,jk->ki',rotation_matrices,vectors) 
100000 loops, best of 3: 5.12 µs per loop 

Das wiederum beweist, warum in NumPy Iterieren ist im Grunde nur schlecht!

Answer 2

Ihre Achsen sind in einer ungewöhnlichen Reihenfolge. Zuerst werden Sie die Matrix setzen wollen Achsen zuletzt:

rotation_matrices = np.rollaxis(rotation_matrices, -1) # shape (50, 2, 2) 
vectors = np.rollaxis(vectors, -1)      # shape (50, 2) 

die es erlauben würden Sie Ihre bestehende Schleife besser lesbar zu machen:

for ind in np.arange(nb_of_vectors): 
    rotated_vector[ind] = np.dot(rotation_matrices[ind], vectors[ind]) 

Aber stattdessen Sie die Matrixmultiplikation Operator (oder np.matmul in Python < 3.5)

rotated_vectors = (a @ vectors[...,None])[...,0] 
# rotated_vectors = np.matmul(a, vectors[...,None])[...,0] 

Die [...,None] wandelt die Anordnung von Vektoren (Form (n,) in eine Anordnung von Spaltenmatrizen (Form (n, 1)) und das nachlauf [...,0] wandelt die Spaltenmatrizen zurück in Vektoren

Answer 3

Dies ist eine explizite Formel

Version

result = np.array([vectors[0,:]*rotation_matrices[0,0,:] + 
        vectors[1,:]*rotation_matrices[0,1,:], 
        vectors[0,:]*rotation_matrices[1,0,:] + 
        vectors[1,:]*rotation_matrices[1,1,:]]).transpose() 

viel (14x) schneller als die Original-Code, aber langsamer (2.6x) als einsum auf meiner Maschine