Welches Prinzip steckt hinter der Konvergenz der Krylov-Subraummethoden zur Lösung linearer Gleichungssysteme?

24

Nach meinem Verständnis gibt es zwei Hauptkategorien iterativer Methoden zum Lösen linearer Gleichungssysteme:

  1. Stationäre Methoden (Jacobi, Gauß-Seidel, SOR, Multigrid)
  2. Krylov-Subraum-Methoden (Conjugate Gradient, GMRES usw.)

Ich verstehe, dass die meisten stationären Methoden durch iteratives Relaxieren (Glätten) der Fourier-Modi des Fehlers funktionieren. Wie ich es verstehe, funktioniert die Methode des konjugierten Gradienten (Krylov-Subraum-Methode), indem sie einen optimalen Satz von Suchrichtungen anhand der Potenzen der Matrix "durchläuft", die auf das te Residuum angewendet werden . Ist dieses Prinzip allen Krylov-Subraummethoden gemeinsam? Wenn nicht, wie charakterisieren wir das Prinzip der Konvergenz von Krylov-Subraummethoden im Allgemeinen?n

Paul
quelle
2
Ihre Analyse stationärer Methoden ist durch einfache Modellprobleme verzerrt, da diese anhand von Fourier-Moden analysiert werden können. Außerdem werden ADI (Alternating Direction Implicit) und viele andere Methoden ignoriert. Bei den meisten "stationären Methoden" geht es darum, viele einfache "approximative partielle" Löser zu einem iterativen Löser zu kombinieren . Der Sinn der Krylov-Methoden ist es, die Konvergenz einer gegebenen stationären linearen Iteration zu beschleunigen (oder sogar zu erzwingen).
Thomas Klimpel
4
Ein Artikel, von dem ich glaube, dass er verfasst wurde, um Ihre Fragen zu beantworten, ist Ipsen und Meyer, Die Idee hinter den Krylov-Methoden, Amer. Mathematik. Monthly 105 (1998), S. 889-899. Es ist ein wunderbar gut geschriebenes und klarstellendes Papier, das hier erhältlich ist .
Andrew T. Barker
@ AndrewT.Barker: Großartig! Vielen Dank Andrew! :)
Paul

Antworten:

21

Im Allgemeinen suchen alle Krylov-Methoden im Wesentlichen ein Polynom, das klein ist, wenn es im Spektrum der Matrix ausgewertet wird. Insbesondere kann der te Rest einer Krylov-Methode (mit Null-Anfangsschätzung) in der Form geschrieben werdenn

rn=Pn(A)b

wobei ein monisches Polynom vom Grad n istPnn .

Wenn diagonalisierbar ist, haben wir mit A = V Λ V - 1AA=VΛV1

rnVPn(Λ)V1b=κ(V)Pn(Λ)b.

Für den Fall, dass normal ist (z. B. symmetrisch oder einheitlich), wissen wir, dass κ ( V ) = 1. GMRES konstruiert ein solches Polynom durch Arnoldi-Iteration, während CG das Polynom unter Verwendung eines anderen inneren Produkts konstruiert (siehe diese Antwort)Aκ(V)=1. für Details). . In ähnlicher Weise konstruiert BiCG sein Polynom durch den nicht symmetrischen Lanczos-Prozess, während die Chebyshev-Iteration vorherige Informationen über das Spektrum verwendet (normalerweise Schätzungen der größten und kleinsten Eigenwerte für symmetrische bestimmte Matrizen).

Betrachten Sie als cooles Beispiel (motiviert von Trefethen + Bau) eine Matrix, deren Spektrum wie folgt lautet:

Spektrum der Matrix

In MATLAB habe ich dies konstruiert mit:

A = rand(200,200);
[Q R] = qr(A);
A = (1/2)*Q + eye(200,200);

Wenn wir GMRES betrachten, das Polynome konstruiert, die tatsächlich den Residuum über alle monischen Polynome des Grades n minimieren , können wir die Residuumshistorie leicht vorhersagen, indem wir das Kandidatenpolynom betrachtenn

Pn(z)=(1z)n

was in unserem Fall gibt

|Pn(z)|=12n

für in dem Spektrum von A .zA

Wenn wir nun GMRES mit einer zufälligen RHS ausführen und die Residuenhistorie mit diesem Polynom vergleichen, sollten sie ziemlich ähnlich sein (die Kandidatenpolynomwerte sind kleiner als die GMRES-Residuen, weil ):b2>1

Restliche Geschichte

Reid. Atcheson
quelle
Können Sie klarstellen, was Sie unter "klein im Spektrum der Matrix" verstehen?
Paul
2
Als komplexes Polynom genommen hat das Polynom einen kleinen Modul in einem Bereich der komplexen Ebene, der das Spektrum von A enthält . Stellen Sie sich ein Konturdiagramm vor, das einem Streudiagramm der Eigenwerte überlagert ist. Wie klein ist klein? Es hängt vom Problem ab, ob A normal ist und die rechte Seite b . Die Grundidee ist jedoch, dass die Folge von Polynomen ( P n ) im Spektrum immer kleiner werden soll, so dass die Restschätzung in meiner Antwort gegen 0 geht . PnAAb.(Pn)0
Reid.Atcheson
@ Reid.Atcheson: Sehr gut ausgedrückt. Könnte ich empfehlen, als κ ( V ) zu schreiben und zu erwähnen, dass es eine für normale Matrizen ist? VV1κ(V)
Jack Poulson
Der mit optimaler SOR vorkonditionierte Laplace-Wert hat ein Spektrum, das dieser Beispielmatrix sehr ähnlich ist. Details hier: scicomp.stackexchange.com/a/852/119
Jed Brown
Streng genommen ist CGNE unabhängig vom Spektrum, da es nur von singulären Werten abhängt.
Jed Brown
17

On norms

nth iteration, GMRES finds the polynomial Pn that minimizes the 2-norm of the residual

rn=Axnb=(Pn(A)1)bb=Pn(A)b.

Suppose A is SPD, so A induces a norm and so does A1. Then

rnA1=rnTA1rn=(Aen)TA1Aen=enTAen=enA

where we have used the error

en=xnx=xnA1b=A1rn

Thus the A-norm of the error is equivalent to the A1 norm of the residual. Conjugate gradients minimizes the A-norm of the error which makes it relatively more accurate at resolving low energy modes. The 2-norm of the residual, which GMRES minimizes, is like the ATA-norm of the error, and thus is weaker in the sense that low-energy modes are less well-resolved. Note that the A-norm of the residual is essentially worthless because it is even weaker on low-energy modes.

Sharpness of convergence bounds

Finally, there is interesting literature regarding different Krylov methods and subtleties of GMRES convergence, especially for non-normal operators.

Jed Brown
quelle
You left off the excellent book by Olavi Nevanlinna: books.google.com/…
Matt Knepley
11

Iterative methods in a nutshell:

  1. Stationary methods are in essence fixed point iterations: To solve Ax=b, you pick an invertible matrix C and find a fixed point of

    x=x+CbCAx
    This converges by Banach's fixed point theorem if ICA<1. The various methods then correspond to a specific choice of C (e.g., for Jacobi iteration, C=D1, where D is a diagonal matrix containing the diagonal elements of A).
  2. Krylov methods subspace methods are in essence projection methods: You pick subspaces U,VCn and look for a x~U such that the residual bAx~ is orthogonal to V. For Krylov methods, U of course is the space spanned by powers of A applied to an initial residual. The various methods then correspond to specific choices of V (e.g., V=U for CG and V=AU for GMRES).

    The convergence properties of these methods (and projection methods in general) follow from the fact that due to the respective choice of V, the x~ are optimal over U (e.g., they minimize the error in the energy norm for CG or the residual for GMRES). If you increase the dimension of U in every iteration, you are guaranteed (in exact arithmetic) to find the solution after finitely many steps.

    As pointed out by Reid Atcheson, using Krylov spaces for U allows you to prove rates of convergence in terms of the eigenvalues (and thus the condition number) of A. In addition, they are crucial for deriving efficient algorithms for computing the projection x~.

    This is nicely explained in Youcef Saad's book on iterative methods.

Christian Clason
quelle