Lösung eines spärlichen und stark schlecht konditionierten Systems

Ich beabsichtige, Ax = b zu lösen, wobei A eine komplexe, spärliche, unsymmetrische und stark schlecht konditionierte (Bedingungsnummer ~ 1E + 20) quadratische oder rechteckige Matrix ist. Ich konnte das System mit ZGELSS in LAPACK genau lösen. Mit zunehmenden Freiheitsgraden in meinem System dauert es jedoch lange, das System auf einem PC mit ZGELSS zu lösen, da die Sparsity nicht ausgenutzt wird. Kürzlich habe ich SuperLU (unter Verwendung von Harwell-Boeing-Speicher) für dasselbe System ausprobiert, aber die Ergebnisse waren für die Bedingungsnummer> 1E + 12 ungenau (ich bin nicht sicher, ob dies ein numerisches Problem beim Schwenken ist).

Ich bin eher geneigt, bereits entwickelte Löser zu verwenden. Gibt es einen robusten Löser, der das erwähnte System schnell (dh unter Ausnutzung der Sparsity) und zuverlässig (angesichts der Bedingungsnummern) lösen kann?

linear-solver sparse-matrix lapack condition-number user1234
quelle

Können Sie es vorkonditionieren? In diesem Fall könnten Krylov-Subraummethoden effektiv sein. Selbst wenn Sie auf direkten Methoden bestehen, hilft die Vorkonditionierung bei der Kontrolle numerischer Fehler.

Geoff Oxberry

Ich habe auch gute Erfahrungen mit der Vorkonditionierung gemacht, so wie es hier beschrieben wird: en.wikipedia.org/wiki/… Sie können die Vorkonditionierung in exakter Arithmetik durchführen. Meine Matrizen sind jedoch alle dicht, sodass Sie hier nicht auf spezifischere Methoden / Routinen verweisen können.

AlexE

Warum ist die Zustandsnummer so groß? Vielleicht kann die Formulierung verbessert werden, um das System besser konditioniert zu machen? Im Allgemeinen können Sie nicht erwarten, dass ein Residuum genauer als

ausgewertet werden kann , was Krylov von geringem Wert macht, wenn Ihnen die Bits ausgehen. Wenn die Bedingungsnummer wirklich

, sollten Sie die Quad-Genauigkeit verwenden ( mit GCC, unterstützt von einigen Paketen, einschließlich PETSc).

(machine precision) \cdot (condition number)

$(\text{machine precision})\cdot(\text{condition number})$

10^{20}

$10^{20}$ __float128

Jed Brown

Woher erhalten Sie diese Schätzung der Bedingungsnummer? Wenn Sie Matlab bitten, die Bedingungsnummer einer Matrix mit einem Nullraum zu schätzen, erhalten Sie möglicherweise eine Unendlichkeit oder manchmal nur eine wirklich große Zahl (wie das, was Sie haben). Wenn das System, das Sie betrachten, einen Nullraum hat und Sie wissen, was es ist, können Sie es projizieren und was Ihnen übrig bleibt, hat möglicherweise eine bessere Bedingungsnummer. Dann können Sie PETSc oder Trilinos verwenden oder was haben Sie.

Daniel Shapero

min ‖ A x - b ‖

$\min \| Ax-b \|$

p e r p (N (A))

$perp(N(A))$

Antworten:

Wenn Sie ZGELSS verwenden, um dieses Problem zu lösen, verwenden Sie die abgeschnittene Singularwertzerlegung, um dieses äußerst schlecht konditionierte Problem zu regulieren. Es ist wichtig zu verstehen, dass diese Bibliotheksroutine nicht versucht, eine Lösung der kleinsten Quadrate für , sondern versucht, eine Lösung zu finden, die minimiert gegen die Minimierung von. $Ax=b$ $\| x \|$ $\| Ax-b \|$

Beachten Sie, dass der an ZGELSS übergebene Parameter RCOND verwendet werden kann, um anzugeben, welche Singularwerte in die Berechnung der Lösung einbezogen und von dieser ausgeschlossen werden sollen. Jeder Singularwert kleiner als RCOND * S (1) (S (1) ist der größte Singularwert) wird ignoriert. Sie haben uns nicht mitgeteilt, wie Sie den RCOND-Parameter in ZGELSS eingestellt haben, und wir wissen nichts über den Rauschpegel der Koeffizienten in Ihrer Matrix oder auf der rechten Seite Daher ist es schwer zu sagen, ob Sie ihn verwendet haben eine angemessene Menge an Regularisierung. $A$ $b$

Sie scheinen mit den regulierten Lösungen, die Sie mit ZGELSS erhalten, zufrieden zu sein. Es scheint also, dass die Regularisierung durch die abgeschnittene SVD-Methode erfolgt (die eine minimale -Lösung unter den Lösungen mit den kleinsten Quadraten findet, die minimieren über den Raum von Lösungen, die von den Singularvektoren überspannt werden, die den Singularwerten zugeordnet sind, die größer als RCOND * S (1) sind, ist für Sie zufriedenstellend. $\| x \|$ $\| Ax-b \|$

Ihre Frage könnte wie folgt umformuliert werden: "Wie kann ich effizient regulierte Lösungen für kleinste Quadrate für dieses große, spärliche und sehr schlecht konditionierte lineare Problem der kleinsten Quadrate erhalten?"

Meine Empfehlung wäre, eine iterative Methode (wie CGLS oder LSQR) zu verwenden, um das explizit regulierte Problem der kleinsten Quadrate zu minimieren

$\min \| Ax-b \|^{2} + \alpha^{2} \| x \|^{2}$

Dabei wird der Regularisierungsparameter so angepasst, dass das Problem der gedämpften kleinsten Quadrate gut konditioniert ist und Sie mit den resultierenden regulierten Lösungen zufrieden sind. $\alpha$

Brian Borchers
quelle

Ich entschuldige mich dafür, dass ich dies zu Beginn nicht erwähnt habe. Das zu lösende Problem ist die Helmholtz-Akustikgleichung mit FEM. Das System ist aufgrund der zur Annäherung der Lösung verwendeten ebenen Wellenbasis schlecht konditioniert.

user1234

Woher kommen die Koeffizienten in und ? Sind sie Messdaten? "exakte" Werte aus der Konstruktion eines Objekts (die in der Praxis nicht mit Toleranzen von 15 Stellen bearbeitet werden können ...)?

A

$A$

b

$b$

Brian Borchers

Die Matrizen A und b werden unter Verwendung der schwachen Formulierung von Helmholtz PDE gebildet, siehe: asadl.org/jasa/resource/1/jasman/v119/i3/…

user1234

Jed Brown hat bereits in den Kommentaren zu der Frage darauf hingewiesen, aber es gibt wirklich nicht viel, was Sie mit üblicher doppelter Genauigkeit tun können, wenn Ihre Bedingungsnummer groß ist: In den meisten Fällen erhalten Sie wahrscheinlich keine einzige Ziffer der Genauigkeit in Ihre Lösung und, schlimmer noch, Sie können es nicht einmal sagen, weil Sie den Rest, der Ihrem Lösungsvektor entspricht, nicht genau auswerten können. Mit anderen Worten: Es ist keine Frage, welchen linearen Löser Sie wählen sollten - kein linearer Löser kann für solche Matrizen etwas Nützliches tun.

Solche Situationen treten normalerweise auf, weil Sie eine ungeeignete Basis wählen. Zum Beispiel erhalten Sie solche schlecht konditionierten Matrizen, wenn Sie die Funktionen als Grundlage für eine Galerkin-Methode wählen. (Dies führt zu der Hilbert-Matrix, die notorisch schlecht konditioniert ist.) Die Lösung besteht in solchen Fällen nicht darin, zu fragen, welcher Löser das lineare System lösen kann, sondern zu fragen, ob es bessere Basen gibt, die verwendet werden können. Ich möchte Sie ermutigen, dasselbe zu tun: Denken Sie darüber nach, Ihr Problem neu zu formulieren, damit Sie nicht mit solchen Matrizen enden. $1,x,x^2,x^3,...$

Wolfgang Bangerth
quelle

Wenn wir ein schlecht gestelltes Problem für eine PDE diskretisieren, z. B. eine Rückwärtswärmegleichung, erhalten wir definitiv eine schlecht konditionierte Matrixgleichung. Dies ist nicht der Fall, den wir lösen können, indem wir die Gleichung neu formulieren oder einen effizienten Matrixlöser auswählen oder die Genauigkeit der Gleitkommazahl verbessern. In diesem Fall [dh akustische inverse Probleme] ist ein Regularisierungsverfahren erforderlich.

17.

Der einfachste / schnellste Weg, um schlecht konditionierte Probleme zu lösen, besteht darin, die Genauigkeit der Berechnungen (durch rohe Gewalt) zu erhöhen. Eine andere (aber nicht immer mögliche) Möglichkeit besteht darin, Ihr Problem neu zu formulieren.

Möglicherweise müssen Sie die vierfache Genauigkeit (34 Dezimalstellen) verwenden. Obwohl in einem Kurs 20 Ziffern verloren gehen (aufgrund der Bedingungsnummer), erhalten Sie immer noch 14 korrekte Ziffern.

Wenn es von Interesse ist, sind jetzt auch spärliche Löser mit vierfacher Genauigkeit in MATLAB verfügbar.

(Ich bin der Autor der genannten Toolbox).

Pavel Holoborodko
quelle