Ich bekomme eine Matrix , die symmetrisch, invertierbar, positiv definitiv und dicht ist. Ich muss testen, ob wobei J die All-One-Matrix ist.Q det ( Q ) = det ( 12 I - Q - J )
Ich mache das gerade mit der Gürteltierbibliothek, aber es stellt sich als zu langsam heraus. Die Sache ist, dass ich dies für eine Billion Matrizen tun muss und es stellt sich heraus, dass die Berechnung der beiden Determinanten der Engpass meines Programms ist. Daher habe ich zwei Fragen
Gibt es einen Trick, mit dem ich die Determinante schneller berechnen könnte, da ich ihre Größe kenne? Ist vielleicht eine unordentliche Erweiterung für Matrizen, die in diesem Fall funktionieren könnte?
Gibt es eine andere effiziente Möglichkeit, die Gleichheit zu testen?
Bearbeiten. Um die Kommentare zu beantworten. Ich muss alle verbundenen nicht selbstkomplementären Graphen der Ordnung so berechnen, dass und die gleiche Anzahl von Spannbäumen haben. Die Motivation dafür finden Sie in diesem Mathoverflow- Beitrag. Die Maschine wird parallel auf einer 8-Kern-Maschine mit 3,4 GHh ausgeführt.
Bearbeiten. Ich konnte die erwartete Laufzeit um 50% reduzieren, indem ich ein C-Programm zur spezifischen Berechnung der Determinante einer Matrix erstellte. Vorschläge sind weiterhin willkommen.
quelle
Antworten:
Da Sie bereits C ++ verwenden und Ihre Matrizen symmetrisch positiv definit sind, würde ich eine nicht drehbare Faktorisierung von und auch von . Hier gehe ich davon aus, dass auch positiv definit ist, andernfalls muss das Gründen der numerischen Stabilität geschwenkt werden (es ist auch möglich, dass das Schwenken nicht erforderlich ist, obwohl es nicht positiv eindeutig ist, aber Sie müssen es versuchen). Q 12 I - Q - J 12 I - Q - J L D L T.LDLT Q 12I−Q−J 12I−Q−J LDLT
Dies ist schneller als eine LU-Faktorisierung und auch schneller als Cholesky, da Quadratwurzeln vermieden werden. Die Determinante ist einfach das Produkt der Elemente der diagonalen Matrix. Der Code ein LDL - Faktorisierung durchzuführen ist so einfach , dass man es in weniger als 50 Zeilen C schreiben können Wikipedia - Seite auf den Algorithmus beschreibt, und ich habe einige einfache Templat Code Cholesky zu tun hier . Sie können dies erheblich vereinfachen und ändern, um die Quadratwurzel zur Implementierung der Faktorisierung zu vermeiden .L D L T.D LDLT
Da Sie auch das Speicherformat steuern können, können Sie die Routine weiter optimieren, um nur die Hälfte der Matrix zu speichern und in ein lineares Array zu packen, um die Speicherlokalität zu maximieren. Ich würde auch einfache benutzerdefinierte Punktprodukt- und Rang-1-Aktualisierungsroutinen schreiben, da die Problemgrößen so klein sind, dass Sie den Compiler die Routinen einbinden lassen sollten, um den Anrufaufwand zu verringern. Da es sich um eine Schleife mit fester Größe handelt, sollte der Compiler in der Lage sein, Dinge bei Bedarf automatisch ein- und auszurollen.
Ich möchte vermeiden , zu spielen Tricks versuchen , sich die Tatsache zunutze zu tragen , dass enthält innerhalb des Ausdrucks. Es ist wahrscheinlich, dass diese Tricks bei solch kleinen Problemgrößen langsamer sind als nur zwei separate Determinantenberechnungen durchzuführen. Die einzige Möglichkeit, diese Behauptungen zu überprüfen, besteht natürlich darin, es zu versuchen.Q.12I−Q−J Q
quelle
Ohne einige Informationen über die Konstruktion dieser positiv definierten reellen symmetrischen Matrizen sind die zu machenden Vorschläge notwendigerweise ziemlich begrenzt.12×12
Ich habe das Armadillo-Paket von Sourceforge heruntergeladen und mir die Dokumentation angesehen. Versuchen Sie, die Leistung der getrennten Berechnung von und zu verbessern , wobei die Rang-1-Matrix aller ist, indem Sie z . In der Dokumentation wird darauf hingewiesen, dass dies die Standardeinstellung für Matrizen bis zu einer Größe von ist. Daher gehe ich davon aus, dass die Option eine Standardeinstellung für den Fall ist.det(Q) det(12I−Q−J) J 4×4 12×12
det(Q,slow=false)
slow=true
WasQ 12×12 Die Determinante wird durchgeführt, weil die Berechnung zu diesem Zeitpunkt an LAPACK (oder ATLAS) "über eine Mauer geworfen" wird, ohne dass darauf hingewiesen wird, dass kein Schwenken erforderlich ist. siehe
slow=true
vermutlich tut, ist ein teilweises oder vollständiges Schwenken, um eine Reihenebenenform zu erhalten, aus der die Determinante leicht gefunden werden kann. Sie wissen jedoch im Voraus, dass die Matrix eindeutig positiv ist, sodass ein Schwenken für die Stabilität nicht erforderlich ist (zumindest vermutlich für den Großteil Ihrer Berechnungen. Es ist unklar, ob das Armadillo-Paket eine Ausnahme auslöst, wenn die Drehpunkte übermäßig klein werden, dies sollte jedoch eine sein Funktion eines vernünftigen numerischen linearen Algebra-Pakets. BEARBEITEN: Ich habe den Armadillo-Code gefunden, der in der Header-Datei implementiert ist , wobei C ++ - Vorlagen für wesentliche Funktionen verwendet wurden. Die Einstellung scheint keinen Einfluss darauf zu haben, wie eindet
include\armadillo_bits\auxlib_meat.hpp
slow=false
det_lapack
und seine Aufrufe in dieser Datei.Der andere Punkt wäre, der Empfehlung zu folgen, das Armadillo-Paket zu erstellen, das mit Hochgeschwindigkeitsersatz für BLAS und LAPACK verknüpft ist, wenn Sie diese tatsächlich verwenden. siehe Kap. 5 der Armadillo README.TXT-Datei für Details. [Die Verwendung einer dedizierten 64-Bit-Version von BLAS oder LAPACK wird auch für die Geschwindigkeit auf aktuellen 64-Bit-Computern empfohlen.]
Die Zeilenreduktion auf die Staffelform ist im Wesentlichen eine Gaußsche Eliminierung und hat eine arithmetische Komplexität . Für beide Matrizen beträgt dies dann das Doppelte dieser Arbeit oder . Diese Operationen mögen zwar der "Engpass" in Ihrer Verarbeitung sein, aber es gibt wenig Hoffnung, dass ohne eine spezielle Struktur in (oder einige bekannte Beziehungen zwischen den Billionen Testfällen, die eine Amortisation ermöglichen) die Arbeit auf reduziert werden könnte .23n3+O(n2) 43n3+O(n2) Q O(n2)
Zum Vergleich beinhaltet die Expansion einer allgemeinen Matrix durch CofaktorenMultiplikationsoperationen (und ungefähr ebenso viele Additionen / Subtraktionen), so dass für der Vergleich ( vs. ) eindeutig die Eliminierung gegenüber Cofaktoren begünstigt.n×n n! n=12 12!=479001600 23n3=1152
Ein anderer Ansatz, der Arbeit erfordert, wäre die Reduzierung von auf die tridiagonale Form mit Householder-Transformationen, wodurch auch in die tridiagonale Form gebracht wird. Das Berechnen von und kann danach in Operationen durchgeführt werden. [Der Effekt der Rang-1-Aktualisierung in der zweiten Determinante kann als ausgedrückt werden, der durch Lösen eines tridiagonalen Systems gegeben ist.]43n3+O(n2) Q 12I−Q det(Q) det(12I−Q−J) O(n) −J
Die Implementierung einer solchen unabhängigen Berechnung kann sich lohnen, um die Ergebnisse erfolgreicher (oder fehlgeschlagener) Aufrufe von Armadillos
det
Funktion zu überprüfen .Sonderfall: Nehmen wir an, dass wobei wie zuvor die (Rang 1) Matrix aller ist und ist nicht singuläre (positive) diagonale Matrix. In der Tat wären dies für die vorgeschlagene Anwendung in der Graphentheorie ganzzahlige Matrizen. Dann lautet eine explizite Formel für :Q=D−J J D=diag(d1,…,dn) det(Q)
Eine Skizze seines Beweises bietet die Möglichkeit, eine breitere Anwendbarkeit zu veranschaulichen, dh wenn eine bekannte Determinante hat und das System schnell gelöst ist. Beginnen Sie mit dem Ausklammern:D Dv=(1…1)T
Jetzt ist wieder Rang 1, nämlich . Beachten Sie, dass die zweite Determinante einfach ist:D−1J (d−11…d−1n)T(1…1)
wobei das charakteristische Polynom von . Als Rang-1-Matrix muss (mindestens) Faktoren von , um seinen Nullraum zu berücksichtigen. Der "fehlende" Eigenwert ist , wie aus der Berechnung ersichtlich ist:f(x) D−1J f(x) n−1 x ∑d−1i
Daraus folgt, dass das charakteristische Polynom und wie oben für , .f(x)=xn−1(x−∑d−1i) f(1) det(I−D−1J) 1−∑d−1i
Beachten Sie auch, dass wenn , dann , eine diagonale Matrix, deren Determinante einfach das Produkt ihrer diagonalen Einträge ist.Q=D−J 12I−Q−J=12I−D+J−J=12I−D
quelle
Wenn Sie eine strukturierte Methode zur Aufzählung der Diagramme haben, für die Sie die Determinanten berechnen möchten, können Sie möglicherweise Aktualisierungen mit niedrigem Rang finden, die Sie von einem Diagramm zum anderen übertragen.
Wenn ja, können Sie das Matrix-Determinanten-Lemma verwenden , um die Determinante des nachfolgenden zu zählenden Graphen unter Verwendung Ihrer Kenntnis der Determinante des aktuellen Graphen kostengünstig zu berechnen.
Das heißt, für eine Matrix und Vektoren : Dies kann verallgemeinert werden, wenn U und V sind Matrizen und ist :A u,v det(A+uvT)=(1+vTA−1u)det(A) n×m A n×n det(A+UVT)=det(Im+VTA−1U)det(A)
Um die Inverse effizient zu berechnen, können Sie die Sherman-Morrison-Formel verwenden , um die Inverse der nachfolgenden Matrix aus der aktuellen zu erhalten:(A+uvT)−1=A−1−A−1uvTA−11+vTA−1u
quelle