Verarbeitung von oberen dreieckigen Elementen nur mit NumPy einsum

Question

Ich benutze numpy einsum, um die Skalarprodukte eines Arrays von Spaltenvektoren pts der Form (3, N) mit sich selbst zu berechnen, die auf einer Matrix dotps der Form (N, N) mit allen Punktprodukten. Dies ist der Code, den ich verwende:

dotps = np.einsum('ij,ik->jk', pts, pts) 

Das funktioniert, aber ich brauche nur die Werte über der Hauptdiagonale. dh. der obere dreieckige Teil des Ergebnisses ohne die Diagonale. Ist es möglich, nur diese Werte mit einsum zu berechnen? oder auf eine andere Weise, die schneller ist als einsum, um die ganze Matrix zu berechnen?

Mein pts-Array kann ziemlich groß sein, also wenn ich nur die Werte berechnen könnte, die ich brauche, würde das meine Rechengeschwindigkeit verdoppeln.

Answer 1

können Sie relevanten Spalten schneiden und dann mit np.einsum -

R,C = np.triu_indices(N,1) 
out = np.einsum('ij,ij->j',pts[:,R],pts[:,C]) 

Probelauf -

In [109]: N = 5 
    ...: pts = np.random.rand(3,N) 
    ...: dotps = np.einsum('ij,ik->jk', pts, pts) 
    ...: 

In [110]: dotps 
Out[110]: 
array([[ 0.26529103, 0.30626052, 0.18373867, 0.13602931, 0.51162729], 
     [ 0.30626052, 0.56132272, 0.5938057 , 0.28750708, 0.9876753 ], 
     [ 0.18373867, 0.5938057 , 0.84699103, 0.35788749, 1.04483158], 
     [ 0.13602931, 0.28750708, 0.35788749, 0.18274288, 0.4612556 ], 
     [ 0.51162729, 0.9876753 , 1.04483158, 0.4612556 , 1.82723949]]) 

In [111]: R,C = np.triu_indices(N,1) 
    ...: out = np.einsum('ij,ij->j',pts[:,R],pts[:,C]) 
    ...: 

In [112]: out 
Out[112]: 
array([ 0.30626052, 0.18373867, 0.13602931, 0.51162729, 0.5938057 , 
     0.28750708, 0.9876753 , 0.35788749, 1.04483158, 0.4612556 ]) 

---

weiter optimieren -

Lassen Sie uns Zeit, unsere Vorgehensweise und sehen, ob es Möglichkeit zur Verbesserung der Leistung.

In [126]: N = 5000 

In [127]: pts = np.random.rand(3,N) 

In [128]: %timeit np.triu_indices(N,1) 
1 loops, best of 3: 413 ms per loop 

In [129]: R,C = np.triu_indices(N,1) 

In [130]: %timeit np.einsum('ij,ij->j',pts[:,R],pts[:,C]) 
1 loops, best of 3: 1.47 s per loop 

innerhalb der Speichereinschränkungen Der Aufenthalt, es sieht nicht wie wir viel über die Optimierung np.einsum tun können. Also, lassen Sie uns den Fokus auf np.triu_indices verschieben.

Für N = 4 haben wir:

In [131]: N = 4 

In [132]: np.triu_indices(N,1) 
Out[132]: (array([0, 0, 0, 1, 1, 2]), array([1, 2, 3, 2, 3, 3])) 

Es scheint ein regelmäßiges Muster zu erstellen, wie eine Art Verschiebung eines though. Dies könnte mit einer kumulativen Summe geschrieben werden, die sich an den Positionen 3 und 5 verschiebt. generisch zu denken, würden wir bis Codierung es so etwas wie dieses Ende -

def triu_indices_cumsum(N): 

    # Length of R and C index arrays 
    L = (N*(N-1))/2 

    # Positions along the R and C arrays that indicate 
    # shifting to the next row of the full array 
    shifts_idx = np.arange(2,N)[::-1].cumsum() 

    # Initialize "shift" arrays for finally leading to R and C 
    shifts1_arr = np.zeros(L,dtype=int) 
    shifts2_arr = np.ones(L,dtype=int) 

    # At shift positions along the shifts array set appropriate values, 
    # such that when cumulative summed would lead to desired R and C arrays. 
    shifts1_arr[shifts_idx] = 1 
    shifts2_arr[shifts_idx] = -np.arange(N-2)[::-1] 

    # Finall cumsum to give R, C 
    R_arr = shifts1_arr.cumsum() 
    C_arr = shifts2_arr.cumsum() 
    return R_arr, C_arr 

Lassen Sie sich Zeit, die für verschiedene N's!

In [133]: N = 100 

In [134]: %timeit np.triu_indices(N,1) 
10000 loops, best of 3: 122 µs per loop 

In [135]: %timeit triu_indices_cumsum(N) 
10000 loops, best of 3: 61.7 µs per loop 

In [136]: N = 1000 

In [137]: %timeit np.triu_indices(N,1) 
100 loops, best of 3: 17 ms per loop 

In [138]: %timeit triu_indices_cumsum(N) 
100 loops, best of 3: 16.3 ms per loop 

So sieht es aus wie für anständige N's, die maßgeschneiderte cumsum basierte triu_indices vielleicht einen Blick wert sein!