Haltekriterien für iterative lineare Löser, die auf nahezu singuläre Systeme angewendet werden

16

Betrachte mit nahezu singulär, was bedeutet, dass es einen sehr kleinen Eigenwert von gibt. Das übliche einer iterativen Methode basiert auf dem Residuum und die Iterationen anhalten können, wenn mit der Iterationsnummer. Aber in dem Fall, den wir betrachten, könnte es einen großen Fehler , der in dem mit dem kleinen Eigenwert assoziierten Eigenraum lebt, was einen kleinen Rest ergibt . Angenommen, der anfängliche Rest ist groß, dann kann es vorkommen, dass wir anhaltenAx=bAλ0Arn:=bAxnrn/r0<tolnvλ0Av=λ0vr0rn/r0<tÖl aber der Fehler xn-x ist immer noch groß. Was ist in diesem Fall eine bessere Fehleranzeige? Ist xn-xn-1Ein guter Kandidat?

Hui Zhang
quelle
3
Vielleicht möchten Sie über Ihre Definition von "fast singulär" nachdenken. Die Matrix (mit und der Identitätsmatrix) hat einen sehr kleinen Eigenwert, ist aber weit davon entfernt, singulär zu sein, wie jede Matrix sein könnte. Iϵϵ1ich
David Ketcheson
1
Auchscheint die falsche Schreibweise. ist typischer, nein? ||rn/r0||||rn||/||r0||
Bill Barth
Ja, du hast recht, Bill! Ich werde diesen Fehler korrigieren.
Hui Zhang
1
Was ist mit? und was ist dein Algorithmus genau? bAx/b
Shuhalo
2
Nachtrag: Ich denke, der folgende Aufsatz befasst sich ziemlich genau mit den schlecht konditionierten Systemen, über die Sie sich Sorgen machen, zumindest wenn Sie CG verwenden: Axelson, Kaporin: Fehler-Normschätzung und Stopp-Kriterien in vorkonditionierten konjugierten Gradienteniterationen. DOI: 10.1002 / nla.244
shuhalo

Antworten:

13

Bitte verwenden Sie niemals den Unterschied zwischen aufeinanderfolgenden Iterationen, um ein Stoppkriterium zu definieren. Diese Fehldiagnose führt zu einer Stagnation der Konvergenz. Die meisten nicht symmetrischen Matrixiterationen sind nicht monoton, und selbst GMRES in exakter Arithmetik ohne Neustart kann für eine beliebige Anzahl von Iterationen (bis zur Dimension der Matrix) stagnieren, bevor sie plötzlich konvergieren. Siehe Beispiele in Nachtigal, Reddy und Trefethen (1993) .

Eine bessere Möglichkeit, die Konvergenz zu definieren,

In der Regel interessiert uns die Genauigkeit unserer Lösung mehr als die Größe des Rests. Insbesondere möchten wir garantieren, dass die Differenz zwischen einer Näherungslösung und der exakten Lösung erfüllt für einige Benutzer angegebenen . Es stellt sich heraus, dass dies erreicht werden kann, indem ein so gefunden wird, dass wobei der kleinste singuläre Wert von , bedingt durch x | x n - x | < C c x n | A x n - b | < C & egr; & egr; Axnx

|xn-x|<c
cxn
|EINxn-b|<cϵ
ϵEIN

|xnx|=|A1A(xnx)|1ϵ|AxnAx|=1ϵ|Axnb|<1ϵcϵ=c

wobei wir verwendet haben, dass der größte Singularwert von (zweite Zeile) ist und dass genau (dritte Zeile) löst .1/ϵA1xAx=b

Abschätzen der kleinste Einzelwertϵ

Eine genaue Schätzung des kleinsten Singularwerts ist normalerweise nicht direkt aus dem Problem verfügbar, kann jedoch als Nebenprodukt eines konjugierten Gradienten oder einer GMRES-Iteration geschätzt werden. Es ist zu beachten, dass, obwohl Schätzungen der größten Eigenwerte und Singularwerte normalerweise nach nur wenigen Iterationen recht gut sind, eine genaue Schätzung des kleinsten Eigen- / Singularwerts normalerweise erst erhalten wird, wenn die Konvergenz erreicht ist. Vor der Konvergenz wird die Schätzung im Allgemeinen deutlich größer sein als der wahre Wert. Dies deutet darauf hin , dass Sie müssen tatsächlich die Gleichungen lösen , bevor Sie die richtigen Toleranz definieren . Eine automatische Konvergenztoleranz , die eine vom Benutzer bereitgestellte Genauigkeit nimmtϵcϵcFür die Lösung und Schätzungen könnte der kleinste singuläre Wert mit dem aktuellen Stand der Krylov-Methode zu früh konvergieren, da die Schätzung von viel größer war als der wahre Wert.ϵϵ

Anmerkungen

  1. Die obige Diskussion funktioniert auch, wenn durch den links-vorkonditionierten Operator und den vorkonditionierten Rest oder wenn der rechts-vorkonditionierte Operator und der Fehler P ( x n - x )AP1AP1(Axnb)AP1P(xn-x) . Wenn ein guter Vorkonditionierer ist, wird der vorkonditionierten Operator gut konditionierte. Für links Präkonditionierung, bedeutet dies den vorkonditionierten Rest klein gemacht werden kann, aber der wahre Rest möglicherweise nicht. Für die richtige Vorkonditionierung giltwird leicht klein gemacht, aber der wahre FehlerP-1|P(xnx)||xnx|möglicherweise nicht. Dies erklärt, warum die linke Vorkonditionierung besser ist, um Fehler klein zu machen, während die rechte Vorkonditionierung besser ist, um den Rest klein zu machen (und um instabile Vorkonditionierer zu debuggen).
  2. In dieser Antwort finden Sie weitere Informationen zu Normen, die durch GMRES und CG minimiert wurden.
  3. Die Schätzungen der extremalen Einzelwerte können unter Verwendung überwacht werden -ksp_monitor_singular_valuemit jedem PETSc Programm. Siehe KSPComputeExtremeSingularValues () singuläre Werte von Code zu berechnen.
  4. Wenn GMRES mit singulären Werten zu schätzen, ist es entscheidend, dass Neustarts nicht verwendet werden (zB -ksp_gmres_restart 1000 , in PETSc).
Jed Brown
quelle
1
'' funktioniert auch, wenn A durch einen vorkonditionierten Operator ersetzt wird '' - Dies gilt jedoch nur für den vorkonditionierten Rest wenn P - 1 A verwendet wird. zum vorkonditionierten Fehler P - 1 δ x , wenn A P - 1 verwendet wird. P1rP1EINP1δxAP-1
Arnold Neumaier
1
Guter Punkt, ich habe meine Antwort bearbeitet. Beachten Sie, dass Sie mit dem rechtskonditionierten Fall steuern können , indem Sie den Vorkonditionierer abwickeln (indem Sie P - 1 anwenden ). PδxP1
Jed Brown
6

Ein anderer Weg , um dieses Problem zu betrachten , ist , die Werkzeuge zu prüfen , aus diskreten inversen Problemen, das heißt, die Probleme beinhalten die Lösung oder min | | A x - b | | 2 , wo A sehr schlecht konditioniert ist (dh das Verhältnis zwischen dem ersten und dem letzten Singulärwert σ 1 / σ nEINx=bMindest||EINx-b||2EINσ1/σn groß ist ).

Hier haben wir verschiedene Methoden zur Auswahl des Abbruchkriteriums, und für ein iteratives Verfahren, würde ich das L-Kurve Kriterium empfehlen, da es nur Mengen handelt, die bereits verfügbar sind (HAFTUNGSAUSSCHLUSS: Mein Berater diese Methode Pionierarbeit geleistet, so dass ich auf jeden Fall bin voreingenommen gegen es). Ich habe dies mit Erfolg in einem iterativen Verfahren verwendet.

Die Idee ist es, die Restnorm zu überwachen und die Lösung Norm η k = | | x k | | 2 , wobei x k das k istρk=||EINxk-b||2ηk=||xk||2xkk -te Iteration ist. Wie Sie Iterierte, diese die Form eines L in einem loglog (rho, eta) Grundstück und an der Stelle , an der Ecke des L zu zeichnen beginnen , ist die optimale Wahl.

Auf diese Weise können Sie ein Kriterium implementieren, bei dem Sie nach dem Überqueren der Ecke ein Auge auf die Steigung von (ρk,ηk) ), und dann die Iterierte wählen , die an der Ecke lag.

Die Art und Weise, wie ich es gemacht habe, beinhaltete das Speichern der letzten 20 Iterationen und wenn der Gradient abs(log(ηk)log(ηk1)log(ρk)log(ρk1)) war größer als ein gewisser Schwellwert 20 für aufeinanderfolgende Iterationen, wusste ich , dass ich auf dem vertikalen Teil der Kurve war und daß ich die Ecke passiert hatte. Ich habe dann die erste Iteration in meinem Array (dh die vor 20 Iterationen) als meine Lösung genommen.

Darüber hinaus gibt es weitere Einzelheiten der Verfahren zur Ecke zu finden, und diese Arbeit besser, aber eine beträchtliche Anzahl von Iterierten erfordern speichern. Spielen Sie mit ihm ein bisschen herum. Wenn Sie in Matlab sind, können Sie die Toolbox Regularisierung Werkzeuge verwenden, die einige dieser implementiert (insbesondere die „Ecke“ Funktion anwendbar ist).

Beachten Sie, dass dieser Ansatz besonders für große Probleme geeignet ist, da die zusätzliche Rechenzeit winzig ist.

OscarB
quelle
1
Danke vielmals! Im loglog (rho, eta) -Diagramm beginnen wir also rechts von der L-Kurve und enden oben auf L, oder? Ich weiß nur nicht , das Prinzip hinter diesem Kriterium. Können Sie erklären , warum es immer wie eine L - Kurve verhalten und warum wir um die Ecke wählen?
Hui Zhang
||EINx-b||2=||e||2ebexeinct=b+e. Für weitere Analysen siehe Hansen, PC & O'Leary, DP (1993). Die Verwendung der L-Kurve bei der Regularisierung diskreter, schlecht gestellter Probleme. SIAM Journal on Scientific Computing, 14. Beachten Sie, dass ich den Beitrag gerade leicht aktualisiert habe.
OscarB
4
@HuiZhang: es ist nicht immer ein L. Wenn die Regularisierung nicht eindeutig ist , kann es ein Doppel - L sein, um zwei Kandidaten für die Lösung führt, eine mit Brutto featurse besser gelöst, das andere mit bestimmten Details besser gelöst. (Und natürlich können weitere komplexe Formen auftreten.)
Arnold Neumaier
Gilt die L-Kurve für schlecht konditionierte Probleme, bei denen es eine eindeutige Lösung geben sollte? Das heißt, ich bin in Problemen interessiert Ax = b , wobei b bekannt ist „genau“ und A fast singulär , aber immer noch technisch umkehrbar. Mir scheint, wenn Sie so etwas wie GMRES verwenden, ändert sich die Norm Ihrer aktuellen Schätzung von x im Laufe der Zeit nicht allzu sehr, insbesondere nach den ersten, jedoch vielen Iterationen. Es scheint mir , dass der vertikale Teil der L-Kurve tritt auf, weil in einem schlecht gestellten Problem keine eindeutige / gültige Lösung ist; Wäre dieses vertikale Merkmal bei allen schlecht konditionierten Problemen vorhanden?
Nukeguy
An einem Punkt, erreichen Sie eine solche vertikale Linie, in der Regel , weil die numerische Fehler in Ihrer Lösungsmethode Ergebnis in || Ax-b || nicht abnehmend. Sie haben jedoch Recht, dass bei solchen rauschfreien Problemen die Kurve nicht immer wie ein L aussieht, was bedeutet, dass Sie in der Regel nur wenige Ecken zur Auswahl haben und eine über die andere wählen kann schwierig sein. Ich glaube, dass das Papier, auf das ich oben in meinem Kommentar verwiesen habe, kurz rauschfreie Szenarien behandelt.
OscarB