Krylow-Unterraum-Verfahren
Krylow-Unterraum-Verfahren sind iterative Verfahren zum Lösen großer, dünnbesetzter linearer Gleichungssysteme, wie sie bei der Diskretisierung von partiellen Differentialgleichungen entstehen, oder von Eigenwertproblemen. Sie sind benannt nach dem russischen Schiffbauingenieur und Mathematiker Alexei Nikolajewitsch Krylow, der die Krylow-Unterräume definierte. Mit den hier beschriebenen Verfahren hat das von ihm entwickelte Verfahren zur Berechnung der Koeffizienten des charakteristischen Polynomes nicht mehr viel zu tun. Die Verfahren sind sogenannte Black-Box-Verfahren, die sich durch einfache Implementierung und Robustheit auszeichnen, weshalb sie vielfach verwendet werden. Die Verfahren sind fast alle Projektions-Verfahren.
Die Verfahrensklasse
Gegeben sei das lineare Gleichungssystem <math>\mathbf{A}x=b</math> mit der Matrix <math>\mathbf{A} \in \mathbb{R}^{n \times n}</math>. Zu einer beliebigen Näherungslösung <math>x_0</math> für <math>x</math> und dem Residuum <math>r_0=b-\mathbf{A}x_0</math> ist der <math>m</math>-te Krylow-Unterraum <math>{\mathcal K}_m</math> der von den Vektoren <math>r_0, \mathbf{A}r_0, \dots, \mathbf{A}^{m-1}r_0</math> aufgespannte Untervektorraum.
Fast alle Verfahren finden nun eine bessere Näherungslösung <math>x_{m} \in x_{0}+{\mathcal K}_m</math> mit der Bedingung, dass der Vektor <math>b-\mathbf{A} x_{m}</math> orthogonal zu allen Vektoren eines Unterraumes <math>{\mathcal L}_m</math> steht, eine sogenannte Orthogonalprojektion. Diese Bedingung heißt Galerkin-Bedingung.
Damit ist das Problem auf ein <math>m</math>-dimensionales lineares Gleichungssystem reduziert. Das Ganze wird zu einem iterativen Lösungsverfahren, wenn man die Dimension in jedem Schritt um eins erhöht.
Spezielle Lösungsverfahren ergeben sich durch die konkrete Wahl des Raumes <math>{\mathcal L}_m</math>, sowie durch Ausnutzen von speziellen Eigenschaften der Matrix <math>\mathbf{A}</math>, was das Verfahren beschleunigt, aber die Anwendbarkeit auch einschränkt. Die einfachste Variante ist, für <math>{\mathcal L}_m</math> einfach wieder den Krylow-Unterraum selbst zu wählen. Das besondere an den Verfahren ist, dass sie nur Matrix-Vektor-Multiplikationen sowie Skalarprodukte benötigen. Ist die Matrix dünnbesetzt, so ist das Matrix-Vektor-Produkt schnell ausrechenbar und der Algorithmus praktikabel.
Ein Beispiel ist das Verfahren der Konjugierten Gradienten (CG-Verfahren). Hierbei ist <math>{\mathcal L}_m={\mathcal K}_m</math> und es ist für symmetrische, positiv definite Matrizen gedacht.
Man erhält so eine Vielfalt an Verfahren. Viel wichtiger als die Auswahl der speziellen Krylow-Unterraummethode ist die Wahl des Vorkonditionierers. Dieser formt das lineare Gleichungssystem äquivalent um, so dass die Lösung unverändert bleibt, sich aber günstigere Eigenschaften für die Konvergenz ergeben. Hier sind entscheidende Geschwindigkeitsgewinne zu erzielen, die dazu führen, dass selbst Systeme mit Millionen Unbekannten in wenigen Dutzend Schritten zufriedenstellend gelöst werden können.
Aufwand
Der Aufwand eines Krylow-Unterraum-Verfahrens setzt sich im Wesentlichen aus Matrix-Vektor-Produkten, Skalarprodukten und Vektoraktualisierungen zusammen. Bei vollbesetzten Matrizen ist das GEMV eine BLAS-Level-2-Operation; bei dünnbesetzten Matrizen spricht man entsprechend von Sparse Matrix-Vector Multiplication (SpMV). Skalarprodukte (DOT) und Vektoraktualisierungen (AXPY) sind BLAS-Level-1-Operationen.
Für eine vollbesetzte Matrix <math>\mathbf{A}\in\mathbb{R}^{n\times n}</math> hat das Matrix-Vektor-Produkt <math>y \leftarrow \mathbf{A}x</math> Kosten von <math>2n^2</math> arithmetische Operationen. Skalarprodukte und Vektoraktualisierungen kosten jeweils <math>2n</math> Operationen.
Dünnbesetzte Matrizen
Für <math>\mathbf{A}</math> dünnbesetzt, bezeichnet
<math> m:= \operatorname{nnz}(\mathbf{A}):= \big| \{(i,j)\in \{1, ... ,n\}^2 \mid A_{ij}\neq 0\} \big| \quad \ll n^2 </math>
die Anzahl der Nichtnulleinträge.
Es kostet ein Matrix-Vektor-Produkt der dünnbesetzten Matrix (SpMV) nur <math>2\cdot \operatorname{nnz}(\mathbf{A})</math> arithmetische Operationen, im Gegensatz zum Aufwand <math>O(n^2)</math> des vollbesetzten Matrix-Vektor-Produkts.
Ist die Zahl der Nichtnulleinträge proportional zur Dimension, also <math>m\leq c n</math> mit einer von <math>n</math> unabhängigen Konstante <math>c</math>, so folgt ein Aufwand <math>2m\leq 2 c n</math>. In diesem Fall ist der Aufwand für SpMV von der Größenordnung <math>O(n)</math>.
Lanczos vs. Arnoldi
Beim Aufwand von Krylow-Unterraum-Verfahren ist der Unterschied zwischen Verfahren mit kurzen Rekursionen und Verfahren mit voller Orthogonalisierung entscheidend. Im symmetrischen bzw. hermiteschen Fall führt der Lanczos-Prozess zu kurzen Dreiterm-Rekursionen. Typische Verfahren wie CG, MINRES oder SYMMLQ benötigen daher pro Iteration nur eine feste Anzahl von Skalarprodukten und Vektoraktualisierungen.
Im allgemeinen nichtsymmetrischen Fall verwendet man dagegen häufig den Arnoldi-Prozess, etwa bei FOM oder GMRES. Dabei muss der neue Krylow-Vektor gegen alle bisherigen Basisvektoren orthogonalisiert werden. In der <math>j</math>-ten Iteration entstehen dadurch zusätzlich etwa <math>j</math> Skalarprodukte und <math>j</math> Vektoraktualisierungen.
Ohne Neustart ergibt sich über <math>k</math> Iterationen daher neben den Matrix-Vektor-Produkten ein zusätzlicher Aufwand der Größenordnung <math>O(k^2 n)</math> und ein Speicherbedarf von <math>O(kn)</math>.
Bei CG bleiben die Rekursionen dagegen kurz: Über <math>k</math> Iterationen entsteht neben den Matrix-Vektor-Produkten nur ein zusätzlicher Aufwand der Größenordnung <math>O(kn)</math>, und der Speicherbedarf bleibt <math>O(n)</math>.
CG-Verfahren
Beim CG-Verfahren ohne Vorkonditionierung besteht eine Iteration aus
- einem Matrix-Vektor-Produkt,
- zwei Skalarprodukten,
- drei Vektoraktualisierungen vom AXPY-Typ.
Für eine dünnbesetzte Matrix mit <math>m=\operatorname{nnz}(\mathbf{A})</math> ergibt sich damit pro Iteration im Wesentlichen
<math>2m + 2\cdot 2n + 3\cdot 2n = 2m + 10n</math> arithmetische Operationen.
Für <math>k</math> Iterationen ergibt sich entsprechend näherungsweise der Gesamtaufwand <math>k(2m+10n)</math>, zuzüglich einmaliger Initialisierungskosten.
Ist die Zahl der Nichtnulleinträge proportional zur Dimension, also <math>m\leq c n</math> mit einer von <math>n</math> unabhängigen Konstante <math>c</math>, so folgt ein Aufwand pro Iteration von der Größenordnung <math>O(n)</math>, genau <math>2m+10n \leq (2c+10)n</math>.
Bei einer vollbesetzten Matrix ergibt sich entsprechend näherungsweise
- <math>2n^2 + 10n</math>.
Bei vorkonditionierten Varianten kommt zusätzlich der Aufwand für die Anwendung des Vorkonditionierers hinzu.
Geschichte
Die ersten Krylow-Unterraumverfahren waren das Lanczos-Verfahren, das 1950 und 1952 von Cornelius Lanczos, das Arnoldi-Verfahren, das 1951 von Walter Edwin Arnoldi, und das CG-Verfahren, das 1952 von Magnus Hestenes und Eduard Stiefel veröffentlicht wurde. Die meisten Krylow-Unterraumverfahren sind sogar direkte Löser, die nach spätestens n Schritten bei exakter Arithmetik die exakte Lösung reproduzieren. Allerdings sind die Verfahren als direkte Löser aufgrund Instabilität bei Rundungsfehlern ungeeignet. Der Nutzen als iterativer Löser wurde damals noch nicht erkannt und so verschwanden die Verfahren in der Schublade. Anfang der 1970er wurde der Nutzen des CG-Verfahrens mit Präkonditionierung dann von Reid erkannt und 1971 eine bestechende Fehleranalyse des symmetrischen Lanczos-Verfahrens von Christopher Conway Paige gegeben. Daraufhin entstand eine Vielzahl neuer, besserer und stabilerer Verfahren wie MinRes, SymmLQ, GMRES und QMR, und gänzlich neue Krylow-Unterraumverfahren wie CGS, BiCG, BiCGSTAB und TFQMR wurden entwickelt. Die heutige Klassifizierung der Krylow-Unterraumverfahren nach den Ansätzen QOR und QMR sowie Verallgemeinerungen des CG-Verfahrens nach den Ansätzen Orthores, Orthomin und Orthodir begannen Ende der 1970er.
Inzwischen gibt es angepasste Krylow-Unterraumverfahren für nichtlineare Eigenwertprobleme, wie den rationalen Krylow von Axel Ruhe und den nichtlinearen Arnoldi von Heinrich Voss. Es existieren auch seit geraumer Zeit Verfahren, welche zur Lösung von nichtlinearen Gleichungssystemen in der inneren Schleife eines Newton-Verfahrens Krylow-Unterraumverfahren verwenden, subsumiert unter dem Schlagwort Newton-Krylow-Methoden oder bei Erwähnung des Vorkonditionieres beispielsweise Newton-Krylow-Schwarz-Methoden.
Alphabetische Liste gängiger Krylow-Unterraum-Verfahren
- Arnoldi-Verfahren, zur Eigenwertapproximation
- BiCG, das CG-Verfahren für nicht SPD-Matrizen
- BiCGSTAB, Stabilisierung von CGS
- BiCGSTAB(ell), Stabilisierung von CGS
- BiCGSTABTFQMR, der Ansatz hinter TFQMR angewandt auf BiCGSTAB
- BiOres, eine Variante des BiCG-Verfahrens
- BiOmin, eine Variante des BiCG-Verfahrens
- BiOdir, eine Variante des BiCG-Verfahrens
- CG, zur approximativen Lösung linearer Gleichungssysteme
- CGNE, CG-Verfahren auf den Normalgleichungen, Variante 1
- CGNR, CG-Verfahren auf den Normalgleichungen, Variante 2
- CGS-Verfahren, quadrierte BiCG-Rekursion
- FOM, zur approximativen Lösung linearer Gleichungssysteme
- GMRES, zur approximativen Lösung linearer Gleichungssysteme
- Hessenberg-Verfahren, zur Eigenwertapproximation
- (Stabilized) Induced Dimension Reduction, Kurzterm-Rekursion und Verfahrensklasse für lineare Löser und Eigenwertlöser
- Lanczos-Verfahren, zur Eigenwertapproximation
- MinRes, zur approximativen Lösung linearer Gleichungssysteme
- Orthores, Orthomin und Orthodir, Verallgemeinerungen des CG-Verfahrens
- Ores, eine Variante des CG-Verfahrens
- Omin, eine Variante des CG-Verfahrens
- Odir, eine Variante des CG-Verfahrens
- Potenzmethode, älteste Methode zur Eigenwertapproximation
- QMR, zur approximativen Lösung linearer Gleichungssysteme
- Richardson-Iteration, bei geeigneter Interpretation
- SymmLQ, zur approximativen Lösung linearer Gleichungssysteme
- TFQMR, zur approximativen Lösung linearer Gleichungssysteme
- Tschebyschow-Iteration, semi-iteratives Verfahren auf Basis von Tschebyschow-Polynomen; eng verwandt mit CG
Literatur
- A. Meister: Numerik linearer Gleichungssysteme, 2. Auflage, Vieweg 2005, ISBN 3-528-13135-7
- Y. Saad: Iterative Methods for Sparse Linear Systems, 2nd edition, SIAM Society for Industrial & Applied Mathematics 2003, ISBN 0-89871-534-2