matrici Moltiplica di ordine superiore con NumPy

Question

ho creato questo problema giocattolo che riflette il mio problema molto più grande:

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. ]] 
""" 

voglio farlo il più velocemente possibile, in modo da utilizzare le funzioni di NumPy per calcolare ans dovrebbe essere l'approccio migliore , poiché questa operazione è pesante e le mie matrici sono piuttosto grandi.

Ho visto this post, ma le forme sono diverse e non riesco a capire quale axes dovrei usare per questo problema. Tuttavia, sono sicuro che tensordot dovrebbe avere la risposta. Eventuali suggerimenti?

EDIT: ho accettato @ajcr's answer, ma si prega di leggere la mia risposta pure, può aiutare gli altri ...

Answer 1

si potrebbe usare np.einsum per fare l'operazione in quanto consente per molto attento controllo su quali assi sono moltiplicati e che sono riassunti:

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

la funzione moltiplica le voci nel primo asse di ind con le voci del primo asse di dist (pedice 'i'). Idem per il secondo asse di ciascun array (pedice 'j'). Invece di restituire un array 3D, si dice ad einsum di sommare i valori lungo l'asse 'j' omettendolo dagli indici di output, restituendo in tal modo un array 2D.

---

np.tensordot è più difficile da applicare a questo problema. Riassume automaticamente i prodotti degli assi. Tuttavia, vogliamo due serie di prodotti ma per sommare solo uno di essi.

Scrivere np.tensordot(ind, dist, axes=[1, 1]) (come nella risposta a cui è collegato) calcola i valori corretti per te, ma restituisce un array 3D con forma (3, 4, 3). Se si può permettersi il costo di memoria di un array più grande, è possibile utilizzare:

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

Questo ti dà il risultato corretto, ma perché tensordot crea una matrice molto più grande del necessario prima, einsum sembra essere una migliore opzione.

Answer 2

seguito @ajcr's great answer, ho voluto per determinare quale metodo è il più veloce, quindi ho usato 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)) 

I sorprendenti risultati sono stati:

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

Così, ovviamente, utilizzando tensordot è stato il modo sbagliato di fallo (per non parlare dello memory error in esempi più grandi, proprio come ha affermato @ajcr).

Da questo esempio era piccolo, ho cambiato la dimensione matrici essere i,j,k = (3000,200,400), capovolto l'ordine solo per essere sicuri che non ha effetto e impostare un altro test con alto numero di ripetizioni:

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)) 

I risultati erano coerenti con la prima esecuzione:

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

Tuttavia, testando un altro tipo di crescita dimensionale - i,j,k = (30000,20,40) ha portato i seguenti risultati:

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

Vedere i commenti per le spiegazioni di questi risultati.

La morale è, quando si cerca la soluzione più veloce per un problema specifico, provare a generare dati che siano il più possibile simili ai dati originali, in termini di tipi e forme. Nel mio caso i è molto più piccolo di j,k e così sono rimasto con la versione brutta, che è anche la più veloce in questo caso.