numpy: arrays addieren mit Besonderheit

[bords0]

Ich habe ein 2-dimensionales quadratisches array, aus dem ich ein neues array der gleichen Größe erzeugen will, indem ich bestimmte Elemente addiere. Am einfachsten als Code:

Code: Alles auswählen

import numpy as np
n = 9
# irgendein n x n - array
A = np.array([np.random.randint(-n, n) for _ in range(n * n)]).reshape(n, n)
# eine Permutation
P = np.random.permutation(n)

B = ???  # was muss hier hin...

B[i, j] == -(A[P[i], i] + A[P[j], j]) + A[P[i], j] + A[P[j], i]  # ...damit das hier gilt?

Nur, dass ich jetzt gerne wüsste, was ich statt der ??? hinschreiben muss, damit die Gleichheit zwei Zeilen weiter unten erfüllt ist
Als Schleife ist das kein Problem, aber langsam. Und es ist so langsam, dass es ein Problem ist. Und bestimmt geht es geschickt mit numpy, aber ich kenne mich nicht genug aus.

Eigentlich muss ich noch ein Matrix dazuaddieren, mit folgender Eigenschaft:

Code: Alles auswählen

C[i, j] == 2 * D[i, j] * E[P[i], P[j]]

Aber das habe ich hinbekommen:

Code: Alles auswählen

C = 2 * D * E[P][:, P]

[MagBen]

Am einfachsten geht es mit einer doppelten Schleife:

Code: Alles auswählen

B = np.empty((n,n))
for i in range(n):
    for j in range(n):
        B[i, j] = -(A[P[i], i] + A[P[j], j]) + A[P[i], j] + A[P[j], i]  

Wobei es mit Numpy oft ohne solche Schleifen geht.

[bords0]

Wobei es mit Numpy oft ohne solche Schleifen geht.

Danke, aber genau nach der Lösung ohne Schleife suche ich. Weshalb ich ja auch schrieb:

Als Schleife ist das kein Problem, aber langsam.

[MagBen]

Eine Schleife kann ich eliminieren:

Code: Alles auswählen

import numpy as np

n = 9
# irgendein n x n - array
A = np.array([np.random.randint(-n, n) for _ in range(n * n)]).reshape(n, n)
# eine Permutation
P = np.random.permutation(n)
 
B = np.empty((n,n))
for i in range(n):
    for j in range(n):
        B[i, j] = -(A[P[i], i] + A[P[j], j]) + A[P[i], j] + A[P[j], i]  
        
B2 = np.zeros((n,n))
for i in range(n):
    B2[i, :] -= A[P[i], i]
    B2[:, i] -= A[P[i], i]
    B2[i, :] += A[P[i], :]
    B2[:, i] += A[P[i], :]

np.testing.assert_array_equal(B,B2)
print("B==B2")

[Sirius3]

ohne Schleife:

Code: Alles auswählen

I = np.arange(n).reshape(-1,1).repeat(n,1)
B = -(A[P[I], I] + A[P[I.T], I.T]) + A[P[I], I.T] + A[P[I.T], I]

[MagBen]

Beeindruckend!

[Sirius3]

Noch zur Vollständigkeit, ein Timing-Test:
2 Schleifen: 0.23s
1 Schleife: 0.32s
0 Schleifen: 0.05s

Durch das viele hin- und herkopieren der Matrix wird der 1-Schleifen-Fall deutlich langsamer.

[MagBen]

für n=9 mag die 1-Schleifenlösungen am langsamsten sein, trotzdem glaube ich die Zeit von 0.3s nicht. Ich hab's nämlich mit n=2000 getestet und habe folgende Zeiten gemessen:
2 Schleifen: 12.4s
1 Schleife: 0.4s
0 Schleifen: 1.9s
Gemessen mit time.time(), die Zeiten sind mit 0.1s Genauigkeit reproduzierbar.

Durch das viele hin- und herkopieren der Matrix wird der 1-Schleifen-Fall deutlich langsamer.

Das Kopieren ist nicht das Problem, sondern das Neuanlegen von temporären Arrays:

1. temporäre nxn Array: A[P, I] + A[P[I.T], I.T]
2. temporäre nxn Array: -(...)
3. temporäre nxn Array: ... + A[P, I.T]
4. temporäre nxn Array: ... + A[P[I.T], I]

Code: Alles auswählen

import numpy as np
import time
 
n = 2000
# irgendein n x n - array
A = np.array([np.random.randint(-n, n) for _ in range(n * n)]).reshape(n, n)
# eine Permutation
P = np.random.permutation(n)

a1 = time.time()
B = np.empty((n,n))
for i in range(n):
    for j in range(n):
        B[i, j] = -(A[P[i], i] + A[P[j], j]) + A[P[i], j] + A[P[j], i]  
a2 = time.time()
print("%.1f" % (a2-a1))

a1 = time.time()
B2 = np.zeros((n,n))
for i in range(n):
    B2[i, :] -= A[P[i], i]
    B2[:, i] -= A[P[i], i]
    B2[i, :] += A[P[i], :]
    B2[:, i] += A[P[i], :]
a2 = time.time()
print("%.1f" % (a2-a1))

#np.testing.assert_array_equal(B,B2)
#print("B==B2")

a1 = time.time()
I = np.arange(n).reshape(n,1).repeat(n,1)
B3 = -(A[P[I], I] + A[P[I.T], I.T]) + A[P[I], I.T] + A[P[I.T], I]
a2 = time.time()
print("%.1f" % (a2-a1))

#np.testing.assert_array_equal(B,B3)
#print("B==B3")

[bords0]

@MagBen, @Sirius3: Danke für die beiden Lösungen!

Bei meinen timing-Tests ist es noch extremer:
Die Ergebnisse sind 90.9, 0.6, 3.7 (Timing-Code von MagBen).

Das war für mich etwas überraschend, dass für große n die Lösung ohne Schleifen deutlich schlechter ist.
Ich habe letztendlich folgende Lösung verwendet, die ohne Schleifen und ohne (große) zusätzliche arrays auskommt (sondern das den broadcast-Mechanismus machen lässt, wenn ich das richtig verstanden habe).

Code: Alles auswählen

AP0 = A[P]
AP = np.diag(AP0)
B4 = AP0 + A.T[:,P] - AP - AP[:, None]

Damit bekomme ich konsistent leicht bessere Werte als die anderen Lösungen; für n=2000 sind es 0.5.

(Eigentlich brauche ich nur n unterhalb 1000, dort sind bei mir die drei nicht-2-Schleifen-Lösungen fast gleich.)

Nochmals Danke für die Tipps!