Permutiere eine Matrix in numpy

27

Ich möchte eine dichte quadratische Übergangsmatrix direkt ändern, indem ich die Reihenfolge mehrerer Zeilen und Spalten mithilfe der Numpy-Bibliothek von Python ändere. Mathematisch entspricht dies einer Vormultiplikation der Matrix mit der Permutationsmatrix P und einer Nachmultiplikation mit P ^ -1 = P ^ T, dies ist jedoch keine rechnerisch sinnvolle Lösung.

Im Moment tausche ich manuell Zeilen und Spalten aus, aber ich hätte erwartet, dass numpy eine nette Funktion f (M, v) hat, wobei M n Zeilen und Spalten und v n Einträge hat, so dass f (M, v) aktualisiert wird M entsprechend der Indexpermutation v. Möglicherweise scheitere ich gerade am Suchen des Internets.

So etwas könnte mit Numpys "Advanced Indexing" möglich sein, aber ich verstehe, dass eine solche Lösung nicht vorhanden wäre. Auch für einige einfache Situationen kann es ausreichend sein, eine Indexpermutation nur separat zu verfolgen, aber dies ist in meinem Fall nicht zweckmäßig.

Hinzugefügt:
Wenn von Permutationen die Rede ist, bedeutet dies manchmal nur das Abtasten von zufälligen Permutationen, beispielsweise als Teil einer Prozedur zum Erhalten von p-Werten in Statistiken. Oder sie bedeuten das Zählen oder Aufzählen aller möglichen Permutationen. Ich spreche nicht über diese Dinge.

Hinzugefügt:
Die Matrix ist klein genug, um in den Arbeitsspeicher des Desktops zu passen, aber groß genug, um sie nicht unnötig zu kopieren. Eigentlich würde ich gerne Matrizen verwenden, die so groß wie möglich sind, aber ich möchte nicht mit der Unannehmlichkeit fertig werden, sie nicht im RAM zu halten, und ich führe O (N ^ 3) LAPACK-Operationen auf der Matrix aus, was auch der Fall wäre Begrenzen Sie die praktische Matrixgröße. Ich kopiere derzeit Matrizen dieser Größe unnötig, aber ich hoffe, dass dies für die Permutation leicht vermieden werden kann.

keiner
quelle
3
Es wäre gut, wenn Sie die Frage aktualisieren könnten, um die Größe Ihrer Matrizen anzugeben. "Gigantisch" bedeutet nicht für alle Menschen dasselbe.
Bill Barth
2
Sie haben Recht, dass die erweiterte (oder so genannte ausgefallene) Indizierung eine Kopie erstellt. Aber wenn Sie akzeptieren, mit dieser Tatsache zu leben, besteht Ihr Code nur darin M[v], die Zeilen zu permutieren.
Daniel Velkov
@ Daniel: Und es wäre M [v,:] [:, v], die ganze Permutation zu machen? Wäre dies der beste Weg, um die Permutation mit ausgefallener Indizierung zu erhalten? Und würde es das Dreifache des Matrixspeichers verwenden, einschließlich der Größe der ursprünglichen Matrix, der zeilen- und spaltenpermutierten Matrix und der temporären zeilenpermutierten Matrix?
Keine
Das ist richtig, Sie hätten Ihre Originalmatrix und 2 Kopien. Übrigens, warum müssen Sie sowohl Zeilen als auch Spalten gleichzeitig permutieren?
Daniel Velkov
4
Was machst du mit der permutierten Matrix? Es ist möglicherweise besser, den Vektor beim Anwenden des Operators einfach zu permutieren.
Jed Brown

Antworten:

9

Laut der Dokumentation gibt es in numpy keine direkte Permutationsmethode wie ndarray.sort .

Ihre Optionen lauten also (vorausgesetzt, es Mhandelt sich um eine Matrix und den Permutationsvektor).N×Np

  1. Implementierung eines eigenen Algorithmus in C als Erweiterungsmodul (aber In-Place-Algorithmen sind zumindest für mich schwierig!)
  2. SpeicheraufwandN

    for i in range(N):
        M[:,i] = M[p,i]
    for i in range(N):
        M[i,:] = M[i,p]
  3. SpeicheraufwandN2

    M[:,:] = M[p,:]
    M[:,:] = M[:,p]

Hoffe, dass diese suboptimalen Hacks nützlich sind.

Stefano M
quelle
@keins ist Hack 2. was nennst du "manuelles Vertauschen von Zeilen und Spalten"?
Stefano M
1
Ich würde die Optionen 1 und 2 kombinieren: Schreibe C-Code, der einen Puffer der Ordnung N verwendet, um jede permutierte Spalte zu schreiben, und schreibe ihn dann dorthin zurück, wo er herkommt. dann machen Sie dasselbe für Zeilen. Wie @Stefano schreibt, benötigt dies nur zusätzlichen Speicher, den Sie bereits verwenden, um die Permutation p zu speichern . O(N)p
Erik P.
O(N)O(N)
2
Dies ist ein wirklich guter Kandidat für eine Cython-Funktion. Es sollten nicht mehr als 10 Zeilen sein. . . willst du, dass ich es knacke?
Meawoppl
Lol. Ich habe damit angefangen, Cython zu verwenden, und dann die richtige Antwort in einer Funktion gefunden, die ich ständig benutze. Doh. Siehe meine gepostete Antwort.
Meawoppl
6

Warnung: Das folgende Beispiel funktioniert ordnungsgemäß, aber die Verwendung des vollständigen Parametersatzes, der am Ende des Dokuments vorgeschlagen wird, macht einen Fehler oder zumindest eine "undokumentierte Funktion" in der Funktion numpy.take () sichtbar. Einzelheiten finden Sie in den Kommentaren unten. Fehlerbericht eingereicht .

Sie können dies direkt mit der take () - Funktion von numpy tun , aber es erfordert ein wenig Reifenspringen.

Hier ist ein Beispiel für eine zufällige Permutation der Zeilen einer Identitätsmatrix:

import numpy as np
i = np.identity(10)
rr = range(10)
np.random.shuffle(rr)
np.take(i, rr, axis=0)
array([[ 0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.],
       [ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.]])

Um dies zu tun, müssen Sie lediglich den "out" -Parameter so festlegen, dass er mit dem Eingabearray identisch ist, UND Sie müssen mode = "clip" oder mode = "wrap" einstellen. Wenn Sie den Modus nicht festlegen, wird eine Kopie erstellt, um den Array-Status in einer Python-Ausnahme wiederherzustellen (siehe hier) .

Abschließend scheint take eine Array-Methode zu sein, also nicht

np.take(i, rr, axis=0)

du könntest anrufen

i.take(rr, axis=0)

wenn das mehr nach deinem geschmack ist. Insgesamt sollten Sie also ungefähr so ​​aussehen:

#Inplace Rearrange
arr = makeMyBixMatrix()
pVec0, pVec1 = calcMyPermutationVectors()
arr.take(pVec0, axis=0, out=arr, mode="clip")
arr.take(pVec1, axis=1, out=arr, mode="clip")

Um sowohl Zeilen als auch Spalten zu permutieren, muss man entweder zwei Mal damit arbeiten oder ein paar hässliche Spielereien mit numpy.unravel_index ziehen , die mir Kopfschmerzen bereiten .

meawoppl
quelle
Wie gesagt, die vorhandenen Algorithmen sind schwierig. Ihre Lösung funktioniert nicht mit NumPy 1.6.2. und 1.7.1 (doppelte Zeilen / Spalten). Hatte keine Zeit zu überprüfen, ob 1.8.x dieses Problem behebt
Stefano M
Hmmm. Können Sie irgendwo Testcode posten? In meinem Kopf habe ich das Gefühl, dass es eine Sortieroperation für die Indizes geben muss, die zuerst vor dem Zupfen erfolgt. Ich werde diese PM genauer untersuchen.
Meawoppl
1
Wenn ich laufe diesen Code ich 1.6.2, test take, not overwriting: True, test not-in-place take: True , test in-place take: False, rr [3, 7, 8, 1, 4, 5, 9, 0, 2, 6], arr [30 70 80 70 40 50 90 30 80 90], ref [30 70 80 10 40 50 90 0 20 60]. So np.takezumindest für numpy 1.6.2 ist nicht bekannt , auf eine in-Place - Permutation und Verwirrungen , Dinge zu tun.
Stefano M
Yeouch. Gut demonstriert. Dies ist wahrscheinlich ein Bug IMHO. Zumindest sollte in den Dokumenten angegeben werden, dass Eingabe und Ausgabe nicht dasselbe Array sein können. Überprüfen Sie dies wahrscheinlich, und es sei denn, dies ist der Fall.
Meawoppl
Einig in Bezug auf den Fehler: Vielleicht sollten Sie Ihrem Beitrag eine Notiz hinzufügen, um die Leser zu warnen, dass Ihre Lösung falsche Ergebnisse liefern kann.
Stefano M
2

Wenn Sie eine dünne Matrix im COOFormat gespeichert haben , kann Folgendes hilfreich sein

    A.row = perm[A.row];
    A.col = perm[A.col];

ACOOpermnumpy.arraymm

Vincent Traag
quelle
aber was ist der Speicheraufwand für das Speichern einer Matrix mit voller Dichte als C00Matrix mit geringer Dichte an erster Stelle?
Federico Poloni
intfloatfloatn2 in einem dichten Fall , so dass man sich wahrscheinlich besser an reguläre numpy.ndarrays halten könnte.
Vincent Traag
1

Ich habe nicht genug Reputation, um einen Kommentar abzugeben, aber ich denke, die folgende SO-Frage könnte hilfreich sein: https://stackoverflow.com/questions/4370745/view-onto-a-numpy-array

Die grundlegenden Punkte sind, dass Sie Basic Slicing verwenden können und eine Ansicht auf das Array erstellen, ohne es zu kopieren. Wenn Sie jedoch Advanced Slicing / Indexing ausführen, wird eine Kopie erstellt.

hadsed
quelle
Das OP bittet um eine Permutation, und dies ist mit dem Basic Slicing nicht möglich.
Stefano M
Sie haben natürlich Recht. Ich dachte, es wäre nützlich für das OP, zu verstehen, was mit dem Schneiden passiert (falls sie es nicht wussten), da sie besorgt waren, wann Kopien passieren würden. Wenn er etwas aus Ihrer Antwort verwenden würde, wäre das meines Erachtens gut zu wissen, da Sie sie in Ihren Loops verwenden.
Hadded
-1

Wie wäre es mit

my_array [:, [0, 1]] = my_array [:, [1, 0]]

Johnsankey
quelle
1
Damit baut er ein Provisorium auf, das genau das ist, was er vermeiden will.
Michael Grant