Welchen Algorithmus würden Sie für die schnelle Lösung eines dichten linearen Systems fester Größe (N = 9) empfehlen (Matrix ist symmetrisch, positiv-semidefinit)?
Typen sind 32 und 64 Bit Fließkommazahlen.
Solche Systeme werden millionenfach gelöst, daher sollte der Algorithmus ziemlich schnell in Bezug auf die Dimension sein (n = 9).
P.S. Beispiele für robuste C ++ - Implementierungen für vorgeschlagene Algorithmen werden geschätzt.
1) Was meinst du mit "millionenfach gelöst"? Die gleiche Koeffizientenmatrix mit einer Million von verschiedenen rechten Begriffen oder einer Million verschiedener Matrizen?
Millionen verschiedener Matrizen.
2) Positiv _semi_definite bedeutet, dass die Matrix singulär sein kann (zur Maschinenpräzision). Wie möchten Sie mit diesem Fall umgehen? Einfach einen Fehler melden oder versuchen, eine vernünftige Antwort zurückzugeben?
Der Fehler ist in Ordnung.
Da die Matrix symmetrisch, positiv-semidefinit ist, ist die Cholesky-Zerlegung der LU-Zerlegung streng überlegen. (etwa doppelt so schnell wie LU, unabhängig von der Größe der Matrix. Quelle: "Numerische Lineare Algebra" von Trefethen und Bau)
Es ist auch de facto der Standard für kleine dichte Matrizen (Quelle: Ich promoviere in rechnerischer Mathematik) Iterative Methoden sind weniger effizient als direkte Methoden, wenn das System nicht groß genug wird (schnelle Faustregel, die nichts bedeutet, aber das ist immer schön zu haben: auf jedem modernen Computer ist jede Matrix kleiner als 100 * 100 definitiv eine kleine Matrix, die direkte Methoden benötigt, anstatt iterative)
Jetzt empfehle ich nicht, es selbst zu tun. Es gibt Tonnen von guten Bibliotheken, die gründlich getestet wurden. Aber wenn ich dir einen empfehlen müsste, wäre Eigen :
By the way, hier in der Dokumentation , haben Sie die verschiedenen Vor- und Nachteile ihrer 7 direkte lineare Löser in einem schönen, übersichtlichen Tisch. Es scheint, dass in Ihrem Fall LDLT (eine Variante von Cholesky)
gewinntIm Allgemeinen ist es am besten, wenn Sie eine vorhandene Bibliothek verwenden, anstatt einen Roll-your-Own-Ansatz zu wählen, da es viele langweilige Details gibt, die bei der Suche nach einer schnellen, stabilen numerischen Implementierung berücksichtigt werden müssen.
Hier sind ein paar, um Sie zu beginnen:
Eigene Bibliothek (meine persönliche Vorliebe):
Ссылка
Armadillo: Ссылка
Suche herum und du wirst viele andere finden.
Ich würde LU-Zerlegung empfehlen, besonders wenn "Millionen von Malen gelöst" wirklich "einmal gelöst und auf Millionen von Vektoren angewendet" bedeutet. Sie erstellen die LU-Dekomposition, speichern sie und wenden die Forward-Back-Substitution auf so viele r.h.s. Vektoren wie Sie wünschen.
Es ist stabiler im Hinblick auf die Abrundung, wenn Sie Pivotieren verwenden.
LU für eine symmetrische Semi-Definite-Matrix macht wenig Sinn: Sie zerstören eine nette Eigenschaft Ihrer Eingabedaten und führen unnötige Operationen durch.
Die Wahl zwischen LLT oder LDLT hängt wirklich von der Konditionsnummer Ihrer Matrizen ab und davon, wie Sie Randfälle behandeln wollen. LDLT sollte nur verwendet werden, wenn Sie eine statistisch signifikante Verbesserung der Genauigkeit nachweisen können oder wenn die Robustheit für Ihre App von größter Bedeutung ist.
(Ohne eine Stichprobe Ihrer Matrizen ist es schwer, fundierte Ratschläge zu geben, aber ich vermute, dass bei einer so kleinen Ordnung N = 9 das Verschwenken der kleinen diagonalen Terme zum unteren Teil von D wirklich nicht notwendig ist fange mit dem klassischen Cholesky an und brich einfach die Faktorisierung ab, wenn die Diag-Terme in Bezug auf eine vernünftig gewählte Toleranz zu klein werden.)
Cholesky ist ziemlich einfach zu programmieren, und wenn Sie nach einem wirklich schnellen Code streben, ist es besser, ihn selbst zu implementieren.
Wie die anderen oben, empfehle ich cholesky. Ich habe festgestellt, dass die erhöhte Anzahl von Additionen, Subtraktionen und Speicherzugriffen bedeutet, dass LDLt langsamer ist als Cholesky.
Es gibt in der Tat eine Zahl eine Anzahl von Variationen über cholesky, und welche am schnellsten ist, hängt von der Darstellung ab, die Sie für Ihre Matrizen wählen. Ich benutze im Allgemeinen eine Fortran-Stil-Darstellung, das ist eine Matrix M ist ein double * M mit M (i, j) ist m [i + dim * j]; Ich nehme an, dass ein oberer dreieckiger Cholesky (ein wenig) der schnellste ist, dh man sucht das obere Dreieck U mit U '* U = M.
Für feste, kleine Dimensionen ist es definitiv eine Überlegung wert, eine Version zu schreiben, die keine Schleifen verwendet. Ein relativ einfacher Weg, dies zu tun, ist ein Programm zu schreiben, um dies zu tun. Wie ich mich erinnere, brauchte es eine Routine, die sich mit dem allgemeinen Fall als Vorlage befasste, nur einen Morgen, um ein Programm zu schreiben, das eine bestimmte Version mit fester Dimension schreiben würde. Die Einsparungen können beträchtlich sein. Zum Beispiel dauert meine allgemeine Version 0,47 Sekunden, um eine Million 9x9-Faktorisierungen zu machen, während die loopless-Version 0,17 Sekunden benötigt - diese Timings laufen auf einem 2,6-GHz-PC mit einem einzigen Thread ab.
Um zu zeigen, dass dies keine große Aufgabe ist, habe ich die Quelle eines solchen Programms unten aufgeführt. Es enthält die allgemeine Version der Faktorisierung als Kommentar. Ich habe diesen Code in Situationen verwendet, in denen die Matrizen nicht in der Nähe von Singular sind, und ich denke, es funktioniert dort ok; aber es mag für eine feinere Arbeit zu grob sein.
%Vor%Wegen seiner Benutzerfreundlichkeit können Sie Eigen-Solver nur zum Vergleich verwenden. Für einen spezifischen Anwendungsfall könnte ein bestimmter Solver schneller sein, obwohl ein anderer besser sein soll. Dazu können Sie Laufzeiten für jeden Algorithmus nur für die Auswahl messen. Danach können Sie die gewünschte Option implementieren (oder eine vorhandene finden, die am besten zu Ihren Bedürfnissen passt).
Tags und Links algorithm math c++ linear-algebra equation-solving