Nicht-triviale Summen von äußeren Produkten ohne Provisorien in numpy

Question

Das eigentliche Problem, das ich lösen wollen, einen Satz von N Einheitsvektoren gegeben wird, und einen weiteren Satz von M Vektoren für jede der Einheits berechnen Vektoren den Mittelwert von der absolute Wert des Skalarprodukts davon mit jedem der Vektoren M. Im Wesentlichen berechnet dies das äußere Produkt der zwei Matrizen und summiert und mittelt mit einem absoluten Wert dazwischen.

Für N und M nicht zu groß ist dies nicht schwer und es gibt viele Möglichkeiten (siehe unten), um fortzufahren. Das Problem ist, wenn die erstellten Provisorien groß sind und eine praktische Begrenzung für den vorgesehenen Ansatz bieten. Kann diese Berechnung ohne Erstellung von Provisorien durchgeführt werden? Die Hauptschwierigkeit, die ich habe, ist auf das Vorhandensein des absoluten Wertes zurückzuführen. Gibt es allgemeine Techniken zum "threading" solcher Berechnungen?

als Beispiel den folgenden Code betrachten

N = 7 
M = 5 

# Create the unit vectors, just so we have some examples, 
# this is not meant to be elegant 
phi = np.random.rand(N)*2*np.pi 
ctheta = np.random.rand(N)*2 - 1 
stheta = np.sqrt(1-ctheta**2) 
nhat = np.array([stheta*np.cos(phi), stheta*np.sin(phi), ctheta]).T 

# Create the other vectors 
m = np.random.rand(M,3) 

# Calculate the quantity we desire, here using broadcasting. 
S = np.average(np.abs(np.sum(nhat*m[:,np.newaxis,:], axis=-1)), axis=0) 

Dieser groß ist, ist S nun ein Array mit der Länge N und enthält die gewünschten Ergebnisse. Leider haben wir dabei einige potentiell riesige Arrays erstellt. Das Ergebnis der

np.sum(nhat*m[:,np.newaxis,:], axis=-1) 

ist ein M X N Array. Das Endergebnis ist natürlich nur von der Größe N. Starten Sie die Größen von N und M zu erhöhen und wir schnell in einen Speicherfehler laufen.

Wie oben erwähnt, wenn der absolute Wert dann nicht erforderlich waren, gehen wir konnten wie folgt nun einsum()

T = np.einsum('ik,jk,j', nhat, m, np.ones(M))/M 

mit Das funktioniert und arbeitet schnell sogar für recht große N und M. Für das spezifische Problem muss ich die abs() einschließen, aber eine allgemeinere Lösung (möglicherweise eine allgemeinere ufunc) würde auch von Interesse sein.

Answer 1

auf einige der Kommentare Basierend darauf, dass cython mit scheint, ist der beste Weg zu gehen. Ich habe törichterweise nie versucht, Cython zu benutzen. Es stellt sich heraus, dass es relativ einfach ist, Arbeitscode zu produzieren.

Nach einiger Suche wir folgenden cython Code zusammen. Dies ist nicht der allgemeinste Code, wahrscheinlich nicht der beste Weg, um es zu schreiben, und kann wahrscheinlich effizienter gemacht werden. Trotzdem ist es nur ungefähr 25% langsamer als der einsum() Code in der ursprünglichen Frage, also ist es nicht zu schlecht! Es wurde geschrieben, um explizit mit Arrays zu arbeiten, die wie in der ursprünglichen Frage erstellt wurden (daher die angenommenen Modi der Eingabe-Arrays).
Trotz der Vorbehalte ist es eine einigermaßen effiziente Lösung für das ursprüngliche Problem liefern und als Ausgangspunkt in ähnlichen Situationen dienen kann.

import numpy as np 
cimport numpy as np 
import cython 
DTYPE = np.float64 
ctypedef np.float64_t DTYPE_t 
cdef inline double d_abs (double a) : return a if a >= 0 else -a 

@cython.boundscheck(False) 
@cython.wraparound(False) 
def process_vectors (np.ndarray[DTYPE_t, ndim=2, mode="fortran"] nhat not None, 
        np.ndarray[DTYPE_t, ndim=2, mode="c"] m not None) : 
    if nhat.shape[1] != m.shape[1] : 
     raise ValueError ("Arrays must contain vectors of the same dimension") 
    cdef Py_ssize_t imax = nhat.shape[0] 
    cdef Py_ssize_t jmax = m.shape[0] 
    cdef Py_ssize_t kmax = nhat.shape[1] # same as m.shape[1] 
    cdef np.ndarray[DTYPE_t, ndim=1] S = np.zeros(imax, dtype=DTYPE) 
    cdef Py_ssize_t i, j, k 
    cdef DTYPE_t val, tmp 
    for i in range(imax) : 
     val = 0 
     for j in range(jmax) : 
      tmp = 0 
      for k in range(kmax) : 
       tmp += nhat[i,k] * m[j,k] 
      val += d_abs(tmp) 
     S[i] = val/jmax 
    return S 

Answer 2

Dies ist ein wenig langsamer, aber schafft nicht die große Zwischenmatrix.

vals = np.zeros(N) 
for i in xrange(N): 
    u = nhat[i] 
    for v in m: 
     vals[i]+=abs(np.dot(u,v)) 
    vals[i]=vals[i]/M 

edit: verschoben Division durch M außerhalb der for-Schleife.

edit2: neue Idee, alte für Nachwelt und relevanten Kommentar beibehalten.

m2 = np.average(m,0) 
vals = np.zeros(N) 
for i in xrange(N): 
    u=nhat[i] 
    vals[i]=abs(np.dot(u,m2)) 

Das ist schnell, aber manchmal gibt unterschiedliche Werte, arbeite ich an, warum, aber vielleicht kann es in der Zwischenzeit helfen.

edit 3: Ah, es ist das absolute Wert Ding.hmm

>>> S 
array([ 0.28620962, 0.65337876, 0.37470707, 0.46500913, 0.49579837, 
     0.29348924, 0.27444208, 0.74586928, 0.35789315, 0.3079964 , 
     0.298353 , 0.42571445, 0.32535728, 0.87505053, 0.25547394, 
     0.23964505, 0.44773271, 0.25235646, 0.4722281 , 0.33003338]) 
>>> vals 
array([ 0.2099343 , 0.6532155 , 0.33039334, 0.45366889, 0.48921527, 
     0.20467291, 0.16585856, 0.74586928, 0.31234917, 0.22198642, 
     0.21013519, 0.41422894, 0.26020981, 0.87505053, 0.1199069 , 
     0.06542492, 0.44145805, 0.08455833, 0.46824704, 0.28483342]) 

time to compute S: 0.000342130661011 seconds 
time to compute vals: 7.29560852051e-05 seconds 

bearbeiten 4: Nun, wenn Sie überwiegend positive Werte für die Einheitsvektoren haben diese schneller ausgeführt werden soll, um die Vektoren in m unter der Annahme, sind immer positiv, wie sie in Ihren Dummy-Daten sind.

Answer 3

Ich glaube nicht, dass es einen einfachen Weg gibt (außerhalb von Cython und dergleichen), um Ihre genaue Operation zu beschleunigen. Vielleicht möchten Sie aber überlegen, ob Sie wirklich berechnen müssen, was Sie berechnen. Denn wenn statt dem Mittelwert der absoluten Werte Sie die root mean square nutzen könnten, würden Sie noch irgendwie Größen der inneren Produkte von durchschnittlich, aber man es in einem einzigen Schuss wie bekommen konnte: als

rms = np.sqrt(np.einsum('ij,il,kj,kl,k->i', nhat, nhat, m, m, np.ones(M)/M)) 

Dies ist die gleiche tun:

rms_2 = np.sqrt(np.average(np.einsum('ij,kj->ik', nhat, m)**2, axis=-1)) 

Ja, es ist nicht genau das, was Sie gefragt, aber ich fürchte, es ist so nah, wie Sie mit einem vektorisiert Ansatz zu bekommen. Wenn Sie diesen Weg zu gehen entscheiden, wie gut np.einsum führt für große N und M: es eine Tendenz, wenn zu viele Parameter und Indizes übergeben zu versinken hat.