schnelle Weg von

Question

Diagonalen von XMX-^T in Python Berechnung muss ich ohne die Diagonalen von XMX-^T berechnen for-Schleife, oder in anderen Worten, für Schleife die folgenden ersetzt:

X = nump.random.randn(10000, 100) 
M = numpy.random.rand(100, 100) 
out = numpy.zeros(10000) 
for n in range(10000): 
    out[n] = np.dot(np.dot(X[n, :], M), X[n, :]) 

Ich weiß, irgendwie sollte ich numpy.einsum verwenden, aber ich konnte nicht herausfinden, wie?

Vielen Dank! Hier

Answer 1

ist ein einfacheres Beispiel:

M = array([[ 0, 4, 8], 
      [ 1, 5, 9], 
      [ 2, 6, 10], 
      [ 3, 7, 11]]) 

X = array([[ 0, 4, 8], 
      [ 1, 5, 9], 
      [ 2, 6, 10], 
      [ 3, 7, 11]]) 

Was Sie suchen - die Summe der Diagonalelemente - häufiger als Spur in der Mathematik bekannt ist. Sie können die Spuren Ihres Matrixproduktes, ohne Schleife erhalten, durch:

In [102]: np.dot(X, np.dot(M,X.T)).trace() 
Out[102]: 692 

Answer 2

Sicher gibt es eine np.einsum Art und Weise, wie so -

np.einsum('ij,ij->i',X.dot(M),X) 

This missbraucht die schnelle Matrix-Multiplikation auf der ersten Ebene mit X.dot(M) und verwendet dann np.einsum, um die erste Achse zu halten und die Summe reduziert die zweite Achse.

Runtime Test -

In diesem Abschnitt werden alle bisher geschrieben Ansätze, das Problem zu lösen.

In [132]: # Setup input arrays 
    ...: X = np.random.randn(10000, 100) 
    ...: M = np.random.rand(100, 100) 
    ...: 
    ...: def original_app(X,M): 
    ...:  out = np.zeros(10000) 
    ...:  for n in range(10000): 
    ...:  out[n] = np.dot(np.dot(X[n, :], M), X[n, :]) 
    ...:  return out 
    ...: 

In [133]: np.allclose(original_app(X,M),np.einsum('ij,ij->i',X.dot(M),X)) 
Out[133]: True 

In [134]: %timeit original_app(X,M) # Original solution 
10 loops, best of 3: 97.8 ms per loop 

In [135]: %timeit np.dot(X, np.dot(M,X.T)).trace()# @Colonel Beauvel's solution 
1 loops, best of 3: 2.24 s per loop 

In [136]: %timeit np.einsum('ij,jk,ik->i', X, M, X) # @hpaulj's solution 
1 loops, best of 3: 442 ms per loop 

In [137]: %timeit np.einsum('ij,ij->i',X.dot(M),X) # Proposed in this post 
10 loops, best of 3: 28.1 ms per loop 

Answer 3

In [210]: X=np.arange(12).reshape(4,3) 
In [211]: M=np.ones((3,3)) 

In [212]: out=np.zeros(4) 
In [213]: for n in range(4): 
    out[n]= np.dot(np.dot(X[n,:],M), X[n,:]) 
    .....:  

In [214]: out 
Out[214]: array([ 9., 144., 441., 900.]) 

One einsum Ansatz:

In [215]: np.einsum('ij,jk,ik->i', X, M, X) 
Out[215]: array([ 9., 144., 441., 900.]) 

anderen einsum Vergleich:

In [218]: timeit np.einsum('ij,jk,ik->i', X, M, X) 
100000 loops, best of 3: 8.98 µs per loop 

In [219]: timeit np.einsum('ij,ij->i',X.dot(M),X) 
100000 loops, best of 3: 11.9 µs per loop 

Das ist ein bisschen schneller, aber die Ergebnisse mit größerer Größe diff können.

einsum spart die Berechnung vieler unnötiger Werte (vgl. Die diagonal oder trace Ansätze).

ähnliche Verwendung von einsum - Combine Einsum Expresions