Gibt es ein numpiges/scipy-Punktprodukt, das nur die diagonalen Einträge des Ergebnisses berechnet?

Question

Stellen 2 numpy Arrays mit:

> A, A.shape = (n,p) 
> B, B.shape = (p,p) 

Typischerweise p ist eine kleinere Anzahl (p < = 200), während n beliebig groß sein kann.

Ich tue folgendes:

result = np.diag(A.dot(B).dot(A.T)) 

Wie Sie sehen können, bin ich nur die n Diagonaleinträge zu halten, aber es ist ein Zwischenprodukt (n x n) Array berechnet, aus dem nur die diagonalen Einträge gehalten werden.

Ich wünsche für eine Funktion wie diag_dot(), die nur die diagonalen Einträge des Ergebnisses berechnet und nicht den gesamten Speicher zuweist.

würde ein Ergebnis sein:

> result = diag_dot(A.dot(B), A.T) 

Gibt es eine vorgefertigte Funktionalität wie diese und diese effizient für die Zuteilung des Zwischen (n x n) Array, ohne die Notwendigkeit zu tun?

Answer 1

Ich glaube, ich habe es auf meine eigene, aber dennoch die Lösung teilen:

da nur die Diagonalen einer Matrix-Multiplikation bekommen

> Z = N.diag(X.dot(Y)) 

entspricht der individuellen Summe des Skalarprodukts Reihen von X und Y von Spalten, wird die vorherige Anweisung entspricht:

> Z = (X * Y.T).sum(-1) 

Für die ursprünglichen Variablen bedeutet dies:

> result = (A.dot(B) * A).sum(-1) 

mich bitte korrigieren, wenn ich falsch bin, aber dies sollte es sein ...

Answer 2

Sie können fast alles bekommen Sie jemals mit numpy.einsum geträumt. Bis starten Sie den Dreh raus zu bekommen, scheint es im Grunde wie schwarze Voodoo ...

>>> a = np.arange(15).reshape(5, 3) 
>>> b = np.arange(9).reshape(3, 3) 

>>> np.diag(np.dot(np.dot(a, b), a.T)) 
array([ 60, 672, 1932, 3840, 6396]) 
>>> np.einsum('ij,ji->i', np.dot(a, b), a.T) 
array([ 60, 672, 1932, 3840, 6396]) 
>>> np.einsum('ij,ij->i', np.dot(a, b), a) 
array([ 60, 672, 1932, 3840, 6396]) 

EDIT Sie eigentlich die ganze Sache in einem einzigen Schuss bekommen, es ist lächerlich ...

>>> np.einsum('ij,jk,ki->i', a, b, a.T) 
array([ 60, 672, 1932, 3840, 6396]) 
>>> np.einsum('ij,jk,ik->i', a, b, a) 
array([ 60, 672, 1932, 3840, 6396]) 

BEARBEITEN Sie wollen es nicht zu viel alleine sehen lassen ... Die Antwort des OPs auf seine eigene Frage zum Vergleich hinzugefügt.

n, p = 10000, 200 
a = np.random.rand(n, p) 
b = np.random.rand(p, p) 

In [2]: %timeit np.einsum('ij,jk,ki->i', a, b, a.T) 
1 loops, best of 3: 1.3 s per loop 

In [3]: %timeit np.einsum('ij,ij->i', np.dot(a, b), a) 
10 loops, best of 3: 105 ms per loop 

In [4]: %timeit np.diag(np.dot(np.dot(a, b), a.T)) 
1 loops, best of 3: 5.73 s per loop 

In [5]: %timeit (a.dot(b) * a).sum(-1) 
10 loops, best of 3: 115 ms per loop 

Answer 3

Eine Fußgänger Antwort, die den Bau von großen Zwischenarrays vermeidet, ist:

result=np.empty([n.], dtype=A.dtype) 
for i in xrange(n): 
    result[i]=A[i,:].dot(B).dot(A[i,:])