Beste Wahl des Lösers für ein großes, spärlich symmetrisches (aber nicht positiv definiertes) System

Ich arbeite derzeit an der Lösung sehr großer symmetrischer (aber nicht positiv definierter) Systeme, die von bestimmten Algorithmen generiert werden. Diese Matrizen haben eine schöne Blocksparsity, die zum parallelen Lösen verwendet werden kann. Ich kann mich jedoch nicht entscheiden, ob ich einen direkten Ansatz (wie Multi-Frontal) oder einen iterativen Ansatz (vorkonditioniertes GMRES oder MINRES) verwenden soll. Alle meine Studien zeigen, dass iterative Löser (selbst bei ziemlich schneller Konvergenz von 7 inneren Iterationen) den direkten '\' - Operator in MATLAB nicht schlagen können. Theoretisch sollten direkte Methoden jedoch teurer sein. Wie passiert das? Gibt es für einen solchen Fall ein aktuelles Dokument oder Papier? Kann ich Blocksparsity in parallelen Systemen mit direkten Methoden genauso effizient einsetzen wie flexible iterative Löser wie GMRES?

Der MUMPS Sparse Direct Solver kann symmetrische unbestimmte Systeme verarbeiten und ist frei verfügbar ( http://graal.ens-lyon.fr/MUMPS/ ). Ian Duff war einer der Autoren von MUMPS und MA57, daher weisen die Algorithmen viele Ähnlichkeiten auf.

MUMPS wurde für Parallelcomputer mit verteiltem Speicher entwickelt, funktioniert aber auch auf Einzelprozessor-Computern. Wenn Sie es mit einer Multithread-BLAS-Bibliothek verknüpfen, können Sie auf Multicore-Prozessoren mit gemeinsamem Speicher angemessene Beschleunigungen erzielen.

Sie haben nicht gesagt, wie groß "sehr groß" ist, aber MUMPS verfügt auch über eine Out-of-Core-Funktion, um Probleme zu lösen, bei denen die faktorisierte Matrix nicht in den verfügbaren Speicher passt.

Das allgemeine Problem, unter dem Direktlöser leiden, ist das Füllphänomen, was bedeutet, dass die Inverse einer dünnen Matrix dicht sein kann. Dies führt zu einem enormen Speicherbedarf, wenn die Struktur der Matrix nicht "geeignet" ist.

Es gibt Versuche, diese Probleme zu luumgehen , und die Standardfunktion von MATLAB verwendet einige davon. Unter http://www.mathworks.de/de/help/matlab/ref/lu.html finden Sie eine Übersicht über Permutationen, diagonale Skalierung und dergleichen. Gleiches gilt natürlich auch für fortgeschrittenere Pakete wie MUMPS ( http://graal.ens-lyon.fr/MUMPS/ ).

Wenn Ihr Problem jedoch groß genug ist, stoßen Sie im Allgemeinen an diese Speichergrenze und müssen iterative (Krylov) Methoden verwenden.

Wenn Ihr Problem symmetrisch und unbestimmt ist, ist MINRES möglicherweise die naheliegende Wahl. Beachten Sie jedoch, dass GMRES und MINRES in exakter Arithmetik dasselbe tun, wenn das Problem symmetrisch ist, GMRES jedoch tendenziell weniger unter dem Verlust der Orthogonalität leidet. Daher betrachten einige Leute GMRES als "die beste Implementierung von MINRES".

In jedem Fall werden Sie wahrscheinlich von der Vorkonditionierung Ihres Systems profitieren. Wenn Sie keine Zeit mit dem Anpassen eines Vorkonditionierers verbringen möchten, bringt Sie ein unvollständiger LU-Faktorisierungs-Vorkonditionierer (ILU) möglicherweise bereits irgendwohin.

Jeder iterative Löser kann direkte Methoden nur schlagen, wenn das Problem ausreichend groß ist (groß, hängt von mehreren Faktoren ab, wie z. B. Speicherbedarf, Effizienz der Implementierung). Und auch jede Krylov-Methode (zum Beispiel GMRES) ist nur dann gut, wenn Sie eine geeignete Vorkonditionierung verwenden (in der Praxis). Wenn Sie keine Vorkonditionierung verwenden, sind Krylov-Methoden im Allgemeinen nicht hilfreich. Auch ich arbeite mit blockarmen Matrizen (diese stammen aus PDE-Anwendungen) und habe festgestellt, dass blockkonditionierte (überlappende additive schwarz) Krylov-Löser (entweder GMRES oder BiCG-Stab neu gestartet) in Verbindung mit einer Grobgitterkorrektur (oder einem Multigrid) für diese skalierbar sind Art der Matrizen.

Der Matlab '\' - Operator ist dank erstklassiger Programmierung hocheffizient. Vielleicht sehen Sie einige der Ideen und wie sie alle möglichen Abkürzungen in Timothy Davis 'Buch verwendet haben.

In Matlab können Sie gmres verwenden, das gut entwickelt und stabil ist. Wahrscheinlich sind Minres, die theoretisch für Ihren Fall ideal sein sollten, möglicherweise nicht so zuverlässig (zumindest meiner Erfahrung nach). Sie sollten die gleiche oder eine höhere Effizienz von matlab gmres erhalten, wenn

1. Ihr System ist groß genug (mindestens einige Tausend mal einige Tausend).
2. Wenn Sie die richtigen Parameter verwenden (RESTART, MAXIT, X0). Hierfür liegen keine klaren Richtlinien vor. Nutzen Sie Ihre Erfahrung.
3. Verwenden Sie einen guten Vorkonditionierer. Sie können ILU oder sogar billiger einen Block Jacobi verwenden. Dies wird den Aufwand erheblich reduzieren.
4. WICHTIGSTES: Wenn Ihre Matrix dünn ist, verwenden Sie das Matlab-Sparse-Format. Dafür ist Matlab gmres ideal gebaut. Dies wird die Kosten erheblich senken.

Verwenden Sie für noch größere Systeme ein Werkzeug wie PETSc.

Wenn Ihre Matrix eine Dimension in der Mitte von Zehntausenden oder weniger hat, verwenden Sie eine direkte Methode, obwohl es nicht viele frei verfügbare direkte Methoden für unbestimmte symmetrische Systeme gibt (eigentlich sind keine, von denen ich weiß, Open Source). Es gibt MA57 von der HSL, aber das ist nur für den akademischen Gebrauch kostenlos. Sie können die Symmetrie sicherlich ignorieren und UMFPACK verwenden .

Bei einer Dimension von etwa zehn Hundert Hundert beginnt die Speichernutzung einer direkten Methode über das hinauszugehen, was ein typischer Desktop-Computer vernünftigerweise verarbeiten kann. Wenn Sie also nicht über eine leistungsfähige Shared-Memory-Maschine verfügen, müssen Sie auf iterative Methoden umsteigen. Für unbestimmte Probleme können Sie BiCG (Bikonjugatgradient) für symmetrische Systeme spezialisieren, obwohl ein Zusammenbruch möglich ist. Es gibt ein kürzlich veröffentlichtes MINRES-QLP für symmetrische Systeme, das mehr numerische Stabilität bietet.

Die beiden Methoden erfordern ziemlich unterschiedliche Implementierungen, da Sie für direkte Methoden tatsächlich die Matrix bilden müssen, während Sie im iterativen Ansatz die Matrix normalerweise nicht explizit bilden.

Es gibt eine Reihe von Gründen, warum ein Ansatz schneller sein kann als der andere, insbesondere in Abhängigkeit von der Matrixdimension. Bei stark schlecht konditionierten Systemen können iterative Methoden ziemlich stark stagnieren. Bei nicht so spärlichen Matrizen erzeugen direkte Methoden am Ende viel Füllmaterial, was die Dinge sehr verlangsamt. Außerdem sind die direkten Methoden in Matlab stark optimiert (es wird intern UMFPACK oder MA57 verwendet), während iterative Methoden normalerweise direkt in Matlab codiert werden und es weniger Möglichkeiten gibt, BLAS der Stufe 3 auszunutzen, da die Engpässe die Anwendung von Matvecs und Punktprodukten sind.