Projektionsmethoden und polynomielle

Methoden zur Lösung linearer Gleichungssysteme

B. Kwella und U. Schendel

Preprint No. A-4/95

Forschungsprojekt ”Parallele Numerik” Institut für Mathematik I, FU Berlin

SERIE A MATHEMATIK

Inhaltsverzeichnis

1 Einleitung 2

2 Projektionsmethoden 3 2.1 Balanzierte Projektionsmethoden . . . . . . . . . . . . . . . . . . . . . . . 7 2.2 Beispiele und Algorithmen . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

3 Polynomielle Methoden 10 3.1 Krylov-Projektions-Methoden . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.2 Beispiele und Algorithmen . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.3 Semi-Krylov-Projektions-Methoden . . . . . . . . . . . . . . . . . . . . . . 17 3.4 Beispiele und Algorithmen . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

4 Nicht polynomielle Projektionsmethoden und nicht projektive polynomielle Methoden 19

5 Zusammenfassung und Ausblick 20

Zusammenfassung

Diese Arbeit gibt einen Überblick über iterative Methoden zur Lösung linearer Gleichungssysteme und nimmt eine Klassifizierung vor. Es werden projektive und polynomielle Ansätze untersucht, die Ideen zur Lösung nichtsymmetrischer Systeme liefern. Aufbauend auf der hier vorgestellten allgemeinen Theorie, wie Minimierungs- und Konvergenzeigenschaften von Projektionsmethoden, werden spezielle Krylov-Projektionsmethoden untersucht und konkrete Algorithmen für unterschiedliche Problemstellungen angegeben.

1 Einleitung

In diesem Artikel werden eine Reihe von iterativen Methoden zur Lösung von linearen Gleichungssystemen der Form

Ax = b  (1)

untersucht, wobei A ∈ IR^{N×N} eine gegebene nichtsinguläre N × N-Matrix, b ∈ IR^{N} ein gegebener Vektor und x ∈ IR^{N} der gesuchte Lösungsvektor ist. In Anlehnung an die Arbeit von Joubert und Manteuffel (siehe [JoMa90]) werden effiziente Lösungverfahren für große, schwach besetzte Systeme, deren Koeffizientenmatrix A in der Regel unsymmetrisch ist, vorgestellt.

Iterative Methoden berechnen bekanntlich ausgehend von einem Startvektor x_{0} ∈ IR^{N} eine Folge von Näherungslösungen {x_{i}}_{i∈IN}, x_{i} ∈ IR^{N}, von (1) mit einem Algorithmus der Form

x_{i} = Θ_{i−1}(x_{i−1}), i ∈ IN,  (2)

wobei Θ_{i−1} die Iterationsvorschrift im i-ten Iterationsschritt darstellt. Weiterhin bezeichne x_{∗} ∈ IR^{N} die exakte Lösung des linearen Gleichungssystems Ax = b, e_{i} = x_{∗} − x_{i} den Fehler im i-ten Iterationsschritt und r_{i} = b − Ax_{i} = Ae_{i} das Residuum.

Man sagt, ein iteratives Verfahren zur Lösung von (1) bricht im i-ten Iterationsschritt ab, wenn die Näherungslösung x_{i} nicht existiert oder nicht eindeutig ist. Ein nichtabbrechendes Verfahren konvergiert, falls bezüglich einer Vektornorm

∥e_{i}∥ = ∥x_{i} − x_{∗}∥ ⟶ 0, i → ∞,  (3)

gilt. Erhält man die exakte Lösung nach einer endlichen Anzahl von Iterationsschritten, so spricht man von einem endlichen Iterationsverfahren, d.h. es gibt ein n ∈ IN, so daß für alle i ≥ n x_{i} = x_{∗} gilt. [KwSch93]

Zur weiteren Bearbeitung werden folgende Notationen benötigt: Zwei Vektoren x, y ∈ IR^{N} heißen A-orthogonal, wenn

(x, Ay) = 0

gilt, wobei mit (·, ·) das euklidische Skalarprodukt im IR^{N} bezeichnet wird. Mit A^{T} bezeichnet man die transponierte Matrix von A. Weiterhin heißt eine Matrix A ∈ IR^{N×N} symmetrisch, falls A = A^{T} gilt. Eine symmetrische Matrix heißt positiv bzw. negativ definit, falls für alle x ∈ IR^{N}\{0}

(x^{T} , Ax) > 0 bzw. (x^{T} , Ax) < 0

gilt. Ist A symmetrisch und positiv definit, dann wird durch

∥x∥_{A} := (x, Ax)^{1/2}

die sogenannte A-Norm definiert.

2 Projektionsmethoden

In diesem Abschnitt werden iterative Methoden behandelt, deren Näherungslösung eine Linearkombination der vorher berechneten Lösung und eines Vektors aus einem Unterraum ist. Für ein i bezeichnet man mit L_{i} einen l_{i} dimensionalen und mit Q_{i} einen q_{i} dimensionalen Unterraum des IR^{N}. Hiermit kann folgende Definition eingeführt werden.

Definition 2.1 Ein iteratives Verfahren der Form (2) mit einem beliebigen Startvektor x_{0} heißt Projektionsmethode, wenn für alle i, i ∈ IN, die Näherungslösung x_{i−1} existiert und x_{i} die folgenden Bedingungen erfüllt:

(i) x_{i} − x_{i−1} ∈ Q_{i}

(ii) x_{i} − x_{∗} ⊥ L_{i}

Man nennt x_{i−1} + Q_{i} auch den Lösungsraum, in dem die Näherungslösung x_{i} liegt. Definition 2.1 macht jedoch keine Aussage über Existenz und Eindeutigkeit von x_{i}. Die durch 2.1 (ii) gegebene Orthogonalitätsbeziehung bezeichnet man auch als Petrov-Galerkin Bedingung. Diese gibt an, daß der Fehler im i-ten Iterationsschritt senkrecht zu allen Vektoren des Unterraums L_{i} steht. Die folgende Abbildung mit L_{i} = Q_{i} verdeutlicht die Zusammenhänge.

[figure: x_{i−1}; x^{∗}; α_{i}Q_{i}; x_{i}; e_{i}; Q_{i} = L_{i}; .]

Zur Konstruktion einer neuen Näherungslösung x_{i} wird x_{i−1} parallel zu Q_{i} verschoben und zwar so weit, daß der Fehler im i-ten Iterationsschritt e_{i} senkrecht zu L_{i} steht.

Sei nun L_{i} = (l_{1}, …, l_{l_{i}}) [ Q_{i} = (q_{1}, …, q_{q_{i}}) ] die N ×l_{i} [ N ×q_{i} ] Matrix, deren Spalten eine Basis des L_{i} [ Q_{i} ] bilden. Die Näherungslösung einer Projektionsmethode kann dann nach 2.1 dargestellt werden als

x_{i} = x_{i−1} + Q_{i}α_{i} α_{i} ∈ IR^{q_{i}}.  (4)

Damit erhält man für das Residuum bzw. den Fehler

r_{i} = b − Ax_{i}
= r_{i−1} − AQ_{i}α_{i}  (5)
e_{i} = x_{∗} − x_{i}
= e_{i−1} − Q_{i}α_{i}.

Wendet man die Petrov-Galerkin Bedingung aus 2.1 auf L_{i} an, ergibt sich als Bedingung an α_{i}

L_{i}^{T} e_{i} = L_{i}^{T} e_{i−1} − L_{i}^{T} Q_{i}α_{i} = 0,

wobei L^{T} die transponierte Matrix zu L_{i} ist. Äquivalent dazu erhält man α_{i} als Lösung des linearen Gleichungssystems

L_{i}^{T} Q_{i}α_{i} = L_{i}^{T} e_{i−1}

mit l_{i} Gleichungen und q_{i} Unbekannten. Für l_{i} < q_{i} ergibt sich kein eindeutiges α_{i} und somit auch keine eindeutige Näherungslösung, während für l_{i} > q_{i} gar keine Lösung existiert, bzw. nur eine mit kleinstem Abstand gefunden werden kann. Ist das Gleichungssystem quadratisch, so ergibt sich sofort folgendes Ergebnis.

Satz 2.1 Ist l_{i} = q_{i}, dann bricht das Verfahren genau dann ab, wenn L_{i}^{T} Q_{i} singulär ist, wenn also Rang L_{i}^{T} Q_{i} < l_{i} bzw. q_{i} ist.

In diesem Zusammenhang ergibt sich eine weitere wichtige Aussage.

Korollar 2.1 Für ein e_{0} sei das Projektionsverfahren nichtabbrechend und aus x_{i−1} ≠ x_{∗} für alle i folgt dim Q_{i−1} < dim Q_{i}. Dann erhält man für ein N × N-Problem die exakte Lösung (bei exakter Arithmetik) nach spätestens N Schritten, das Verfahren ist also endlich.

Beweis : Die Näherungslösungen eines nicht abbrechenden Projektionsverfahren existieren und sind nach Definition eindeutig. Somit ist

x_{i} = x_{∗}
⟺ e_{i} = x_{∗} − x_{i} = 0
⟺ e_{i} = x_{∗} − x_{i−1} − α_{i}Q_{i}
= e_{i−1} − α_{i}Q_{i} = 0
⟺ e_{i−1} ∈ Q_{i}

Da dim Q_{i} streng monoton wächst, existiert ein n ≤ N mit e_{n−1} ∈ Q_{n}. Somit ist x_{n} = x_{∗} mit n ≤ N die exakte Lösung. □

Existiert nun die Inverse von L_{i}^{T} Q_{i}, ist also insbesondere l_{i} = q_{i}, dann gilt

α_{i} = (L_{i}^{T} Q_{i})^{−1}L_{i}^{T} e_{i−1}.  (6)

Somit lassen sich (4) und (5) schreiben als

x_{i} = x_{i−1} + Q_{i}(L_{i}^{T} Q_{i})^{−1}L_{i}^{T} e_{i−1}  (7)

und

r_{i} = (I − Q_{i}(L_{i}^{T} Q_{i})^{−1}L_{i}^{T})r_{i−1}
e_{i} = (I − Q_{i}(L_{i}^{T} Q_{i})^{−1}L_{i}^{T})e_{i−1}.  (8)

Definiert man nun

P_{i} := I − Q_{i}(L_{i}^{T} Q_{i})^{−1}L_{i}^{T}  (9)

so ist \underline{P_{i} eine Projektion}, da

P_{i}^{2} = I − 2Q_{i}(L_{i}^{T} Q_{i})^{−1}L_{i}^{T} + Q_{i}(L_{i}^{T} Q_{i})^{−1}(L_{i}^{T} Q_{i})(L_{i}^{T} Q_{i})^{−1}L_{i}^{T}
= P_{i}

gilt. Damit ergibt sich die für Projektionsmethoden charakteristische Darstellung des Fehlers

e_{i} = P_{i} e_{i−1}

bzw. des Residuums

r_{i} = P_{i} r_{i−1}.

Im Folgenden wird nun untersucht, unter welchen Bedingungen ein Projektionsverfahren nicht abbrechend ist. Es kann immer eine Matrix B_{i} ∈ IR^{N×N} konstruiert werden mit L_{i} = B_{i}^{T} Q_{i}, wenn L_{i} und Q_{i} die gleiche Dimension ≠ 0 besitzen. Dabei ist B_{i}^{T} := L_{i}(Q_{i}^{T} Q_{i})^{−1}Q_{i}^{T}. Mit Hilfe von B_{i} = Q_{i}(Q_{i}^{T} Q_{i})^{−1}L_{i}^{T} kann (9) geschrieben werden als

P_{i} = I − Q_{i}(Q_{i}^{T} B_{i}Q_{i})^{−1}Q_{i}^{T} B_{i}.

Die Matrix B_{i} beschreibt die Abhängigkeit zwischen den beiden Unterräumen Q_{i} und L_{i}. Ist B_{i} definit bzgl. Q_{i}, also v^{T} B_{i}v ≠ 0 für alle v ∈ Q_{i}\{0} gilt, so ist auch Q_{i}^{T} B_{i}Q_{i} = L_{i}^{T} Q_{i} definit und somit nichtsingulär, das Verfahren bricht nach Satz 2.1 im i-ten Iterationsschritt nicht ab. Dieses Ergebnis motiviert den folgenden Satz:

Satz 2.2 (Beschränktheit) Sei l_{i} = q_{i} ≠ 0 und L_{i} = B_{i}^{T} Q_{i} für alle i, wobei Q_{i}, L_{i} und B_{i} von e_{0} abhängig sein können. Weiterhin nimmt man an, daß ein c_{i}, welches unabhängig von e_{0} ist, existiert mit der Eigenschaft

0 < ∥B_{i}∥∥y∥^{2} ≤ c_{i}|y^{T} B_{i}y|

für alle e_{0} und alle y ∈ Q_{i}\{0}. Dann bricht das Verfahren in keinem Schritt ab und es gilt die Abschätzung

∥x_{i}∥ ≤ ∥x_{i−1}∥ + c_{i}∥e_{i−1}∥.

Beweis : Sei y = Q_{i}v, dann ist y^{T} B_{i}y = v^{T} L_{i}^{T} Q_{i}v. Nach Voraussetzung gilt nun :

\frac{v^{T} L_{i}^{T} Q_{i}v}{v^{T} v} \frac{∥v∥^{2}}{∥Q_{i}v∥^{2}} = \frac{y^{T} B_{i}y}{∥y∥^{2}} ≥ \frac{∥B_{i}∥}{c_{i}} > 0.

Daraus folgt, daß L_{i}^{T} Q_{i} definit und somit nichtsingulär ist, das Verfahren bricht in diesem Schritt also nicht ab. Um die Abschätzung zu zeigen, betrachtet man eine Näherungslösung

x_{i} = x_{i−1} + Q_{i}(Q_{i}^{T} B_{i}Q_{i})^{−1}Q_{i}^{T} B_{i}e_{i−1}.

Setzt man hier Q̃_{i} = Q_{i}R mit R quadratisch und Q̃_{i}^{T} Q̃_{i} = I, dann kann man x_{i} schreiben als

x_{i} = x_{i−1} + Q̃_{i}(Q̃_{i}^{T} B_{i}Q̃_{i})^{−1}Q̃_{i}^{T} B_{i}e_{i−1}.

Dann gilt die Abschätzung

∥x_{i}∥ ≤ ∥x_{i−1}∥ + ∥Q̃_{i}(Q̃_{i}^{T} B_{i}Q̃_{i})^{−1}∥∥Q̃_{i}^{T} B_{i}e_{i−1}∥.

Um eine untere Schranke für ∥x_{i}∥ zu finden, betrachtet man weiterhin

∥Q̃_{i}(Q̃_{i}^{T} B_{i}Q̃_{i})^{−1}∥^{2} = 1/ min_{∥w∥=1} ∥Q̃_{i}(Q̃_{i}^{T} B_{i}Q̃_{i}w∥^{2}
= 1/ min_{∥w∥=1} w^{T} Q̃_{i}(Q̃_{i}^{T} B_{i}Q̃_{i}(ww^{T} + (I − ww^{T}))Q̃_{i}(Q̃_{i}^{T} B_{i}Q̃_{i}w
= 1/ min_{∥w∥=1}(|w^{T} Q̃_{i}(Q̃_{i}^{T} B_{i}Q̃_{i}w|^{2} + ∥(I − ww^{T}))Q̃_{i}(Q̃_{i}^{T} B_{i}Q̃_{i}w∥^{2})
≤ 1/ min_{∥w∥=1} |w^{T} Q̃_{i}(Q̃_{i}^{T} B_{i}Q̃_{i}|^{2}
≤ 1/ min_{y∈Q_{i},∥y∥=1} |y^{T} B_{i}y|^{2}.

Dabei folgt die letzte Ungleichung mit y = Q̃_{i}w aus der Tatsache, daß ∥Q̃_{i}w∥ = ∥w∥ gilt. Beachtet man die Voraussetzungen an B, so erhält man schließlich

∥x_{i}∥ ≤ ∥x_{i−1}∥ + sup_{y∈Q_{i}\\{0\}} \frac{∥y∥^{2}}{|y^{T} B_{i}y∥}∥B_{i}∥∥e_{i−1}∥
≤ ∥x_{i−1}∥ + c_{i}∥e_{i−1}∥.

Zusammenfassend kann festgestellt werden, daß ein Projektionsverfahren im i-ten Iterationschritt nicht abbricht, wenn die Matrix B_{i}, die die Beziehung zwischen den zwei dimensionsgleichen Unterräumen L_{i} und Q_{i} beschreibt, definit ist bzw. gewissen Schranken unterliegt. Wächst außerdem die Dimension der Unterrräume in jedem Schritt, so ist das Verfahren endlich.

2.1 Balanzierte Projektionsmethoden

In diesem Abschnitt werden Projektionsmethoden untersucht, wo es möglich ist, die Abhängigkeit der Unterräume durch eine feste Matrix B darzustellen. Ihre besonderen Eigenschaften werden vorgestellt.

Definition 2.2 Man nennt eine Projektionsmethode balanziert, falls eine feste Matrix B ∈ IR^{N×N} existiert, so daß für alle i und e_{0} gilt:

L_{i} = B^{T} Q_{i}.

Die Petrov-Galerkin Bedingung e_{i} ⊥ L_{i} aus Definition 2.1 (ii) ist somit äquivalent zu e_{i} ⊥ B^{T} Q_{i}, also gilt

(Be_{i}, z) = 0 ∀z ∈ Q_{i}.

Ist die Matrix B definit, also y^{T} By ≠ 0 für alle y ≠ 0, so sind die Voraussetzungen von Satz 2.2 erfüllt. Beachtet man, daß eine reelle Matrix B genau dann definit ist, wenn ihr symmetrischer Anteil H(B) := \frac{B^{T} +B}{2} positiv bzw. negativ definit ist (siehe [JoYo87]), so ergibt sich folgendes Ergebnis sofort aus Satz 2.2.

Korollar 2.2 (Beschränktheit) Eine balanzierte Projektionsmethode bricht in keinem Schritt ab, falls B definit ist. Gilt außerdem

ρ_{min}(B) = inf_{x≠0} \frac{|y^{T} By|}{y^{T} y},

dann ist

∥x_{i}∥ ≤ ∥x_{i−1}∥ + \frac{∥B∥}{ρ_{min}(B)}∥e_{i−1}∥.

Beweis : Der Beweis ist eine direkte Anwendung von Satz 2.2. □

Wiederum erhält man ein wichtiges Resultat, wenn die Dimension der Unterräume L_{i} und Q_{i} in jedem Iterationsschritt wächst.

Satz 2.3 (gleichmäßige Beschränktheit) Seien die Vorausetzungen von Korollar 2.2 erfüllt und gilt zusätzlich Q_{i} ⊆ Q_{i+1} für alle i, dann existiert für alle i die Schranke

∥x_{i}∥ ≤ ∥x_{0}∥ + \frac{∥B∥}{ρ_{min}(B)}∥e_{0}∥.

Beweis : Der Beweis kann analog zum Beweis von Satz 2.2 geführt werden unter Beachtung, daß eine Näherungslösung x_{i} unter der Voraussetzung Q_{i} ⊆ Q_{i+1} folgende Form hat:

x_{i} = x_{0} + Q_{i}(Q_{i}^{T} B_{i}Q_{i})^{−1}Q_{i}^{T} B_{i}e_{0}.
□

Ist die Matrix B selbst symmetrisch und positiv definit, dann ist der Projektor P_{i} = Q_{i}(Q_{i}^{T} B_{i}Q_{i})^{−1}Q_{i}^{T} B_{i} orthogonal bzgl. des inneren Produktes (·, ·)_{B} := (B·, ·). In diesem Fall wird durch ∥ · ∥_{B} eine Norm definiert und man findet eine Näherungslösung, die den Fehler bzgl. dieser Norm minimiert.

Satz 2.4 (Minimierung) Basiert eine balanzierte Projektionsmethode auf einer Matrix B, die symmetrisch und positiv definit ist, dann ist x_{i} das eindeutige Minimum des Fehlers ∥e_{i}∥_{B} über dem affinen Lösungsraum x_{i−1} + Q_{i} und es gilt ∥e_{i}∥_{B} ≤ ∥e_{i−1}∥_{B}.

Beweis : Unter obigen Voraussetzungen und Beachtung der Darstellung e_{i} = P_{i}e_{i−1} erhält man nun, daß e_{i} die B-orthogonale Projektion von e_{i−1} auf das Komplement von Q_{i} ist und es gilt:

e_{i}^{T} Be_{i} = e_{i−1}^{T}Be_{i−1} − e_{i−1}BQ_{i}(Q_{i}^{T} BQ_{i})^{−1}Q_{i}^{T} Be_{i−1}
≤ e_{i−1}^{T}Be_{i−1}.

Man nimmt nun an, daß eine weitere Näherungslösung x′ aus dem Lösungsraum existiert. Sei e′ = x_{∗} − x′ der dazugehörige Fehler, dann ist für ein w

e′ = e_{i} + Q_{i}w.

Wendet man die Eigenschaften des Projektors an, so erhält man

e^{′ T} Be′ = e_{i}^{T} Be_{i} + w^{T} Q_{i}^{T} BQ_{i}w
> e_{i}^{T} Be_{i} für w ≠ 0,

was zeigt, daß x_{i} das eindeutige Minimum ist. □

Eine Aussage über die Konvergenz von Projektionsverfahren macht der folgende Satz, der ein Theorem von [EiES83] verallgemeinert.

Satz 2.5 (Konvergenz) Für eine balanzierte Projektionsmethode, die auf einer symmetrisch positiv definiten Matrix B basiert, für die Z := BA^{−1} definit und außerdem r_{i−1} ∈ Q_{i} ist, gilt: Es existiert ein ϵ > 0, unabhängig von e_{0}, so daß für alle i die Abschätzung

∥e_{i}∥_{B} ≤ (1 − ϵ)∥e_{i−1}∥_{B}

gültig ist.

Beweis: Ist B symmetrisch und positiv definit, dann wird der Fehler e_{i} über dem affinen Unterraum e_{i−1} + Q_{i} bzgl. der B-Norm minimiert. Sei nun

ẽ_{i} = e_{i−1} − λAe_{i−1},

wobei λ durch die Beziehung Bẽ_{i}⊥Ae_{i−1} = r_{i−1} definiert ist. Falls dies eine Minimierung über einem kleineren Unterraum nach sich zieht, muß gelten:

∥e_{i}∥_{B} ≤ ∥ẽ_{i}∥_{B}.

Man erhält

λ = \frac{(Ae_{i−1})^{T} B(e_{i−1})}{(Ae_{i−1})^{T} B(Ae_{i−1})}.

Mit den projektiven Eigenschaften ergibt sich

∥ẽ_{i}∥_{B}^{2} = ∥e_{i−1}∥_{B}^{2} − \frac{|e_{i−1}^{T}BAe_{i−1}|^{2}}{(Ae_{i−1})^{T} B(Ae_{i−1})}.

Unter Vernachlässigung der Indizes und Beachtung von Ae = r erhält man schließlich

\frac{∥ẽ∥_{B}^{2}}{∥e∥_{B}^{2}} = 1 − \frac{|e^{T} BAe|^{2}}{(e^{T} Be)(e^{T} A^{T} BAe)}
= 1 − \frac{|r^{T} A^{−1 T} Br|^{2}}{(r^{T} A^{−1 T} BA^{−1}r)(r^{T} Br)}
= 1 − \frac{|r^{T} Zr|^{2}}{(r^{T} A^{−1 T} BA^{−1}r)(r^{T} Br)}
≤ 1 − \frac{ρ_{min}(Z)^{2}}{∥A^{−1}∥^{2} ∥B∥^{2}}

mit

ρ_{min}(M) = min_{v≠0} \frac{v^{T} Mv}{v^{T} v}

für jede Matrix M, die definit ist genau dann, wenn ρ_{min}(M) ≠ 0 ist. □

Zusammenfassend kann nun für balanzierte Projektionsmethoden festgestellt werden, daß sie, basierend auf einer definiten Matrix B, nicht abbrechend sind. Ist die Matrix B weiterhin sogar symmetrisch und positiv definit, so erhält man Verfahren, die eine minimale Näherungslösung erzeugen und konvergent sind.

2.2 Beispiele und Algorithmen

In diesem Abschnitt wird die Konstruktion eines speziellen Projektionsverfahren vorgestellt, welches in [GrPfSch91] genau analysiert wird.

Man definiert sich zunächst zu einer Näherungslösung x_{i−1} eine Hyperebene E_{i} ∈ IR^{N−1} der Form

E_{i} := {v ∈ IR^{N}|p_{i−1}^{T}Av − p_{i−1}^{T}b = 0},

mit einem geeigneten p_{i−1}. Wählt man Q_{i} := A^{T} p_{i−1}, so ist Q_{i} orthogonal zu E_{i}. Charakteristisch für Projektionsmethoden ist eine Näherungslösung x_{i} somit darstellbar als

x_{i} = x_{i−1} + λ_{i−1}A^{T} p_{i−1}, λ_{i−1} ∈ IR,

sie wird also parallel zu A^{T} p_{i−1} um den Betrag λ_{i−1} verschoben und somit orthogonal auf E_{i} projiziert. Da aber auch x_{i} ∈ E_{i} ist, gilt

λ_{i−1} = −\frac{p_{i−1}^{T}r_{i−1}}{p_{i−1}^{T}AA^{T} p_{i−1}} .

Beim Projektionsverfahren von Kaczmarz unterteilt man jeden Iterationsschritt in N Teilschritte, d.h. man wählt p_{i−1} := e_{j}, j = 1, …, N, den j-ten Einheitsvektor und erhält so

x_{i} = x_{i−1} − \frac{e_{j}r_{i−1}}{∥A^{T} e_{j}∥}A^{T} e_{j}, j = 1, …, N.

Man kann dieses Verfahren verallgemeinern, in dem man den Vektor p_{i−1} durch eine Matrix V_{i} ersetzt, es wird also nicht mehr unbedingt auf eine Hyperebene, sondern auf einen Unterraum beliebiger Dimension projiziert.

Seine Anwendung findet das Projektionsverfahren von Kaczmarz im Zusammenhang mit blockiterativen Methoden, so daß aus ihnen parallele Algorithmen gewonnen werden können.

3 Polynomielle Methoden

In diesem Kapitel werden iterative Methoden vorgestellt, die ebenso wie Projektionsmethoden eine Näherungslösung als Linearkombination der vorher berechneten Lösung und eines Vektors aus einem Unterraum bestimmen. Hier wird jedoch a priori keine Orthogonalitätseigenschaft des Fehlers gefordert. Zur Konstruktion der Lösungsräume werden zunächst einige Begriffe eingeführt.

Mit Π_{k}^{n} bezeichne man die Menge aller Polynome P_{k} vom Grad ≤ k, für die P_{k}(0) = n gilt. Damit ergibt sich folgende Definition.

Definition 3.1 Sei r_{i} = b−Ax_{i} ≠ 0 das Residuum eines Iterationsverfahrens zur Lösung von (1) mit einem beliebigen Startvektor x_{0}. Dann heißt ein Polynom R_{ik} ∈ Π_{ik}^{1}, k ∈ IN, Residuenpolynom von A, wenn für das Residuum

r_{i} = R_{ik}(A)r_{0},

gilt.

Besitzt das Residuum eine Residuenpolynomdarstellung, dann gilt auch für den Fehler

e_{i} = R_{ik}(A)e_{0}.

Mit diesem Begriff kann nun folgende Definition eingeführt werden.

Definition 3.2 Ein iteratives Verfahren der Form (2) mit einem beliebigen Startvektor x_{0} heißt polynomielle (polynomial-based) Methode k-ten Grades, wenn für das Residuum eine Residuenpolynomdarstellung existiert und das Residuenpolynom im i-ten Iterationsschritt vom Grad ik ist.

In jedem Iterationsschritt ist daher das Residuenpolynom vom Grad k.

Ebenso wie die Definition des Residuenpolynoms verallgemeinert auch Definition 3.2 die sonst in der Literatur gebräuchlichen Begriffe. Einige Autoren, wie beispielsweise Elman und Saad ([Elma82], [Saad87]), definieren ein Residuenpolynom nur für den Fall, daß k = 1 gilt, woraus die in Abschnitt 3.1 gezeigten speziellen Eigenschaften folgen.

Zuvor jedoch wird noch eine generelle Aussage bezüglich der Konvergenz dieser Verfahren angegeben.

Lemma 3.1 Ein Iterationsverfahren der Form (2), das außerdem eine polynomielle Methode k-ten Grades ist, konvergiert, falls

lim_{i→∞} ∥R_{ik}(A)∥ = 0  (10)

gilt.

Beispielsweise bestimmen Tschebyscheff-Methoden, die zu den polynomiellen Methoden gehören, ein Residuenpolynom in Abhängigkeit vom Spektrum der Matrix A so, daß dessen Norm möglichst klein ist. Die im nächsten Abschnitt untersuchten Iterationsverfahren verwenden andere Ansätze zur Konstruktion eines geeigneten Residuenpolynoms.

3.1 Krylov-Projektions-Methoden

Polynomielle Methoden ersten Grades besitzen also im i-ten Iterationsschritt ein Residuenpolynom vom Grade ≤ i, also r_{i} = R_{i}(A)r_{0}. Somit gilt

b − Ax_{i} = R_{i}(A)r_{0}
b − r_{0} − (c_{1}A + … + c_{i}A^{i})r_{0} = Ax_{i}
Ax_{0} − (c_{1}A + … + c_{i}A^{i})r_{0} = Ax_{i}.

Da A regulär ist, liegt eine Näherungslösung x_{i} in

x_{i} ∈ x_{0} + span {r_{0}, Ar_{0}, A^{2}r_{0}, …, A^{i−1}r_{0}}
x_{i} ∈ x_{0} + K_{i}(r_{0}, A).  (11)

K_{k}(x, A) = span {x, Ax, A^{2}x, …, A^{k−1}x} nennt man Krylov-Raum. Deshalb bezeichnet man polynomielle Verfahren ersten Grades auch als Krylov-Methoden. Der Krylov-Raum hat die Dimension k, wenn die durch A ∈ IR^{N×N} und x ∈ IR^{N}\{0} erzeugten Vektoren {x, Ax, A^{2}x, …, A^{k−1}x}, k ∈ IN, k ≤ N linear unabhängig sind. Hierzu betrachtet man das Minimalpolynom eines Vektors x ∈ IR^{N}\{0} bezüglich der Matrix A, welches das Polynom kleinsten Grades P_{k} ∈ Π_{k} mit k ≤ N ist, für das

P_{k}(A)x = 0

gilt. Das Minimalpolynom P_{k}(A) zu x ≠ 0 ist genau dann vom Grad k, wenn für jedes j < k die Vektoren {x, Ax, A^{2}x, …, A^{j}x} linear unabhängig und die Vektoren {x, Ax, A^{2}x, …, A^{k}x} linear abhängig sind. Somit gilt K_{i}(x, A) ⊂ K_{i+1}(x, A) für alle i ≤ k.

Eine Näherungslösung x_{i} einer Krylov-Methode ist nach (11)

x_{i} = x_{i−1} + α_{i−1}p_{i−1} α_{i−1} ∈ IR, p_{i−1} ∈ K_{i}(r_{0}, A).  (12)

Somit ist

e_{i} = x_{∗} − x_{i} = x_{∗} − x_{i−1} − α_{i−1}p_{i−1}
= e_{i−1} − α_{i−1}p_{i−1},
r_{i} = r_{i−1} − α_{i−1}Ap_{i−1}.

Mit p_{0} := r_{0} liegt r_{i} ∈ K_{i}(r_{0}, A), besitzt also eine Residuenpolynomdarstellung r_{i} = R_{i}(A)r_{0} und erfüllt so Definition 3.2. Viele polynomielle Verfahren zur Lösung des linearen Gleichungssystem Ax = b basieren auf dem Prinzip anstelle von (1) ein äquivalentes Minimierungsproblem zu formulieren. Sie lösen in jedem Iterationsschritt ein Optimierungsproblem, d.h. sie bestimmen α_{i} derart, daß der Fehler bzw. das Residuum bzgl. irgendeiner Vektornorm über einem affinen Unterraum minimal ist. Die bekanntesten Vertreter dieser Methoden sind Abstiegsverfahren. Charakteristisch für diese Klasse von Iterationsverfahren ist der Abstieg in Gradientenrichtung der zu minimierenden Funktion. In Abschnitt 3.2 werden konkrete Algorithmen angegeben, die ihre Richtungsvektoren in Abhängigkeit von Eigenschaften der Koeffizientenmatrix A bzw. der durch die Unterräume bestimmten Matrix B wählen. Beispielsweise konstruiert man A-orthogonale bzw. A^{T} A-orthogonale Richtungsvektoren p_{i}, die zugehörigen Residuen ergeben sich dann als orthogonal bzw. A-orthogonal zueinander.

Zur Konstruktion polynomieller Methoden, die außerdem projektive Eigenschaften besitzen, wählt man Q_{i} = K_{i}(r_{0}, A) und bestimmt die Richtungsvektoren p_{i} derart, daß die Petrov-Galerkin Bedingung

Be_{i}⊥K_{i}(r_{0}, A)

erfüllt ist. Diese Krylov-Methoden sind somit auch balanzierte Projektionsmethoden; man nennt sie Krylov-Projektionsmethoden. Aus der Orthogonalitätsbeziehung e_{i}⊥B^{T} K_{i}(r_{0}, A) und der Existenz von α_{i}p_{i} mit α_{i} ≠ 0, erhält man sofort, daß

p_{i} ∈ K_{i+1}(r_{0}, A)
p_{i} ⊥ B^{T} K_{i}(r_{0}, A)  (13)

gilt.

Polynomielle Projektionsverfahren ersten Grades sind endlich. Dies folgt einerseits aus Satz 2.1, da sie balanzierte Projektionsmethoden sind und die Dimension der Unterräume in jedem Schritt wächst. Anderseits erhält man dieses Ergebnis für polynomielle Iterationsverfahren, die Näherungslösungen mit orthogonalen Residuen in einem Krylov-Raum bestimmen, durch Betrachtungen des Minimalpolynoms.

Satz 3.1 Die Näherungslösung x_{n} einer polynomielle Methode ersten Grades mit orthogonalen Residuen ist gleich der exakten Lösung x_{∗} des Gleichungssystems Ax = b (bei exakter Arithmetik) genau dann, wenn der Grad des Minimalpolynoms von r_{0} bezüglich der Matrix A gleich n ist.

Beweis : Sei nun x_{n} = x_{∗}, dann gilt

r_{n} = R_{n}(A)r_{0} = 0.

Dies bedeutet aber, daß die Vektoren {r_{0}, Ar_{0}, A^{2}r_{0}, …, A^{n}r_{0}} linear abhängig sind, da sie eine nichttriviale Darstellung der Null bilden.

Weiterhin gilt für alle j < n

r_{j} = R_{j}(A)r_{0} ≠ 0,

und nach obigen Betrachtungen folgt die lineare Unabhängigkeit der Vektoren {r_{0}, Ar_{0}, A^{2}r_{0}, …, A^{j}r_{0}}. Weiterhin gilt dies aber genau dann, wenn das Minimalpolynom von r_{0} bezüglich A den Grad n hat.

Umgekehrt impliziert aber die lineare Abhängigkeit der {r_{0}, …, A^{n}r_{0}} und die lineare Unabhängigkeit der {r_{0}, …, A^{n−1}r_{0}} die Gleichheit der durch sie aufgespannten Krylov-Räume K_{n}(r_{0}, A) = K_{n−1}(r_{0}, A). Nach Konstruktion ist r_{n} ∈ K_{n}(r_{0}, A) = K_{n−1}(r_{0}, A) und außerdem orthogonal zu K_{n−1}(r_{0}, A). Dies kann aber nur erfüllt sein, wenn r_{n} = 0 ist. Weiterhin gibt es aber auch kein j < n mit r_{j} = 0, da dies im Widerspruch zur linearen Unabhängigkeit steht. □

Da im IR^{N} aber N + 1-Vektoren linear abhängig sind und deshalb

r_{N} = R_{N}(A)r_{0} = c_{0}r_{0} + c_{1}Ar_{0} + ⋯ + c_{N}A^{N}r_{0} = 0

gilt, erhält man daraus sofort Korollar 3.1 Für ein N×N-Problem erhält man nach spätestens N Schritten die Lösung.

Weitere wichtige Aussagen über Krylov-Projektionsmethoden ergeben sich analog zu Satz 2.2 bis 2.5. Auch hier erhält man Resultate, die abhängig von Eigenschaften der Matrix B sind. So ist eine polynomielle Projektionsmethode mit einer festen Matrix B genau dann nicht abbrechend, wenn B definit ist. Dies bedeutet aber auch, ist B indefinit, so bricht das Verfahren im ersten Iterationsschritt ab, da für ein r_{0} ≠ 0 gilt: r_{0}^{T} Br_{0} = 0. Allgemeiner kann festgestellt werden:

Satz 3.2 Ein Krylov-Projektionsmethode mit einem beliebigen Startresiduum r_{0}, basierend auf der Matrix A und einer festen Matrix B, sei nichtabbrechend im i-ten Schritt und es gelte x_{i} ≠ x_{∗}. Dann bricht das Verfahren auch im i + 1-ten Iterationsschritt genau dann nicht ab, wenn für den, bis auf Vielfache, eindeutigen Richtungvektor p_{i}, der die Bedingungen (13) erfüllt, gilt:

p_{i}^{T} Bp_{i} ≠ 0.

Beweis: Ein Induktionsargument liefert: In jedem Iterationsschritt existiert der eindeutige Richtungsvektor p_{i} ∈ K_{i+1}(r_{0}, A) und erfüllt p_{i}⊥B^{T} K_{i}(r_{0}, A) genau dann, wenn p_{j}^{T} Bp_{j} ≠ 0 für j < i. Ein eindeutiges α_{i} existiert nun genau dann, wenn p_{i}^{T} Bp_{i} ≠ 0 gilt. □

Wendet man nun die Orthogonalitätsbeziehung e_{i+1}⊥B^{T} K_{i+1}(r_{0}, A) auf p_{i} ∈ K_{i+1}(r_{0}, A) an, so erhält man sofort (e_{i}, Bp_{i}) − α_{i}(p_{i}, Bp_{i}) = 0, woraus sich unter Beachtung von (p_{i}, Bp_{i}) ≠ 0

α_{i} = \frac{(e_{i}, Bp_{i})}{(p_{i}, Bp_{i})}  (14)

ergibt. Weiterhin existiert ein Vektor p_{i} ≠ 0 nur, wenn i kleiner als der Grad des Minimalpolynoms von r_{0} bzgl. der Matrix A ist. Aus (13) und (14) erhält man, daß die Richtungsvektoren B orthogonal sind, d.h. p_{j}^{T} Bp_{i} = 0 für alle j < i. Somit sind sie insbesondere linear unabhängig und bilden eine Basis für den Raum span {p_{0}, …, p_{i−1}} = span {r_{0}, …, r_{i−1}} = span {r_{0}, Ar_{0}, …, A^{i−1}r_{0}} = K_{i}(r_{0}, A).

Wählt man für eine Krylov-Projektionsmethode einen geeigneten Startvektor x_{0}, so daß der Fehler e_{0} und x_{0} gewissen Schranken unterliegen, und ist wiederum die Matrix B definit, so ergibt sich analog zu Satz (2.2) auch hier die Beschränktheit der Folge von Näherungslösungen. Auch die Eigenschaft der Minimierung überträgt sich für eine symmetrisch und positiv definite Matrix B. Man erhält also eine eindeutig bestimmte Näherungslösung x_{i}, die den Fehler e_{i} bzgl. der B-Norm über dem Lösungsraum x_{0} +K_{i}(r_{0}, A) minimiert. Für weitere Diskussionen sei auf [JoMa90] verwiesen.

Wie schon oben erwähnt, ergeben sich durch die unterschiedliche Wahl von B, die natürlich auch abhängig ist von Eigenschaften der Koeffizientenmatrix A, verschiedene Iterationsverfahren, auf die im nächsten Abschnitt näher eingegangen wird.

3.2 Beispiele und Algorithmen

Krylov-Projektionsmethoden können in drei wichtige Unterklassen, abhängig von Eigenschaften, die die Matrix B erfüllt, eingeteilt werden. Verfahren, die durch eine symmetrisch positiv definite Matrix B charakterisiert sind, heißen konjugierte Gradientenmethoden. Das klassische konjugierte Gradienten-Verfahren (cg-Verfahren) liefert für eine symmetrisch positiv definite Matrix A eine eindeutige Näherungslösung x_{i}, die den Fehler über dem Unterraum minimiert. Genaue Untersuchungen findet man beispielsweise bei [Hest80].

Für einen beliebigen Startvektor x_{0} ∈ IR^{N} und r_{0} = b − Ax_{0} führt man folgende Iterationen mit p_{0} = r_{0} durch:

x_{i} = x_{i−1} + α_{i−1}p_{i−1}
mit α_{i−1} = \frac{(r_{i−1}, p_{i−1})}{(p_{i−1}, Ap_{i−1})}
r_{i} = r_{i−1} − α_{i−1}p_{i−1}
p_{i} = r_{i} − β_{i−1}p_{i−1}
mit β_{i−1} = \frac{(r_{i}, Ap_{i−1})}{(p_{i−1}, Ap_{i−1})}

Die durch das cg-Verfahren erzeugten Richtungsvektoren p_{i} sind für eine symmetrische positiv definite Koeffizientenmatrix A-orthogonal, d.h. (p_{i}, Ap_{j}) = 0, für i ≠ j. Insbesondere sind die Vektoren {p_{0}, …, p_{i−1}} linear unabhängig und spannen einen Unterraum des IR^{N} auf. Die Residuen r_{i} sind orthogonal und somit insbesondere linear unabhängig. Das cg-Verfahren ist eine polynomielle Methode ersten Grades, da die Menge der Richtungsvektoren bzw. der Residuen einen Krylov-Raum der Dimension i aufspannen, d.h. r_{i} ∈ K_{i}(r_{0}, A), es existiert also eine Residuenpolynomdarstellung r_{i} = R_{i}(A)r_{0}. Eine Näherungslösung x_{i} liegt in x_{0} + K_{i}(r_{0}, A) und minimiert die Funktion E(x) = ∥e∥_{A}^{2} über dem affinen Krylovraum. Ebenso gehört das cg-Verfahren mit B = A zu den balanzierten Projektionsmethoden mit einer symmetrisch positiv definiten Matrix B. Der Fehler e_{i} steht somit senkrecht zu AK_{i}.

Das konjugierte Residuen-Verfahren (cr-Verfahren) ist eine Variante des cg- Verfahrens für eine symmetrisch positiv definite Matrix A, bei der nicht der Fehler, sondern das Residuum minimiert wird. Der Grundalgorithmus des cg-Verfahrens bleibt bei diesem Verfahren erhalten. Da jedoch hier das Optimum der Funktion F (x) := ∥r∥_{2}^{2} zu bestimmen ist, unterscheiden sich die α_{i} und β_{i} von denen des cg-Verfahrens. Wählt man α_{i} = \frac{(r_{i},Ap_{i})}{(Ap_{i},Ap_{i})} und β_{i} = \frac{(Ar_{i−1},Ap_{i})}{(Ap_{i},Ap_{i})}, so konstruiert man A^{T} A-orthogonale Richtungsvektoren p_{i} und die Residuen sind A-orthogonal.

Das cr-Verfahren ist auch eine Krylov-Projektionsmethode mit B = A^{T} A symmetrisch und positiv definit. Man beachte, daß die Hesse-Matrix HF (x) := \frac{δ^{2}F}{δx_{i}δx_{j}} hier für beliebiges nichtsinguläres A wieder symmetrisch und positiv definit ist, das Minimierungsproblem also für eine größere Klasse von Aufgaben definiert ist.

Das cr-Verfahren bzw. die Minimierung der Funktion F (x) bildet nun die Grundlage für eine Klasse von Lösungsverfahren mit einer nichtsymmetrischen N × N-Matrix A, deren symmetrischer Anteil H(A) positiv definit ist.

Die vom cr-Verfahren erzeugten Richtungsvektoren sind A^{T} A-orthogonal, falls A symmetrisch und positiv definit ist. Ist die Symmetrie von A nicht mehr gegeben, so verliert man im allgemeinen diese Orthogonalität. Dies führt auf die zweite Unterklasse von polynomiellen Projektionsmethoden, den Orthogonalen Fehlermethoden, charakterisiert durch eine unsymmetrische, aber definite Matrix B sind.

Beispielsweise erzeugt das allgemeine cr-Verfahren (gcr-Verfahren) ([EiES83]) eine Menge A^{T} A-orthogonaler Richtungsvektoren, indem alle vorher berechneten Richtungsvektoren benutzt werden, um p_{i+1} zu berechnen. Es ergibt sich der folgende Algorithmus, der auf eine definite Matrix A, deren symmetrischer Anteil H(B) zusätzlich positiv definit sein muß, angewandt werden kann.

x_{i} = x_{i−1} + α_{i−1}p_{i−1}
mit α_{i−1} = \frac{(r_{i−1}, Ap_{i−1})}{(Ap_{i−1}, Ap_{i−1})}
r_{i} = r_{i−1} − α_{i−1}Ap_{i−1}
p_{i} = r_{i} − ∑_{j=0}^{i−1} β_{j}^{(i)}p_{j}
mit β_{j}^{(i)} = \frac{(Ar_{i}, Ap_{j})}{(Ap_{j}, Ap_{j})}

Weiterhin sind die Residuen des gcr-Verfahrens A-orthogonal und besitzen eine Residuenpolynomdarstellung. Mit B = A definit ist es eine Krylov-Projektionsmethode, die die Residuennorm über dem affinen Krylov-Raum x_{0} + K_{i}(r_{0}, A) minimiert. Allgemein kann jedoch über die Minimierungseigenschaft von orthogonalen Fehlermethoden keine Aussage gemacht werden. Der Nachteil des gcr-Verfahrens, der es praktisch undurchführbar macht, liegt in dem großen Arbeitsaufwand pro Iterationsschritt und vor allem in dem hohen Speicherbedarf, da alle bisher berechneten Richtungsvektoren gespeichert werden müssen. Es bildet jedoch die Grundlage für iterative Projektionsmethoden, die in Abschnitt 3.4 vorgestellt werden.

Die dritte Unterklasse polynomieller und projektiver Iterationsverfahren basieren auf einer symmetrischen, nicht notwendig definiten Matrix B; man bezeichnet sie als Symmetrische Krylov-Projektionsmethoden. Beachtet man jedoch, daß wesentliche Eigenschaften nur für definite Matrizen B gelten, können diese Verfahren abbrechend und ihre Näherungslösungen unbeschränkt sein. Sie werden deshalb in diesem Zusammenhang nicht näher untersucht. Ein klassisches Beispiel ist die SYMMLQ-Methode, die in [PaSa75] genauer beschrieben ist.

3.3 Semi-Krylov-Projektions-Methoden

Eine Möglichkeit um beispielsweise unökonomische orthogonale Fehlermethoden praktisch durchführbar zu machen, bzw. für bestimmte A überhaupt nichtabbrechende Verfahren zu entwickeln, bieten polynomielle Methoden k-ten Grades. Sie besitzen nach Definition 3.2 in jedem Iterationsschritt ein Residuenpolynom vom Grad ≤ k, d.h. es gilt: r_{i} = R_{ik}(A)r_{0}. Zu einem fest vorgegeben k ∈ IN mit k << N wird hier in jedem Iterationsschritt die Fehler- bzw. Residuennorm über einem k-dimensionalen Unterraum minimiert. Das ursprüngliche Problem der Optimierung über dem ganzen Krylov-Raum wird dadurch verlagert auf das Lösen einer k-dimensionalen Minimierungsaufgabe in jedem Iterationsschritt. Mit r_{i} = R_{ik}(A)r_{0} ergibt sich

b − Ax_{i} = R_{k}(A)r_{i−1}
b − r_{i−1} − (c_{1}A + … + c_{k}A^{k})r_{i−1} = Ax_{i}
Ax_{i−1} − (c_{1}A + … + c_{k}A^{k})r_{i−1} = Ax_{i}

und daraus folgt für eine Näherungslösung

x_{i} ∈ x_{i−1}K_{k}(r_{i−1}, A).

In jedem Iterationsschritt spannen die Richtungsvektoren bzw. die Residuen einer polynomiellen Methode k-ten Grades einen Krylov-Raum der Dimension k auf. Man nennt sie daher Semi-Krylov-Methoden. Eine Näherungslösung x_{i} liegt in x_{0}+K_{ik}(r_{0}, A), jedoch wird in der Regel nicht der ganzen Raum genutzt. Für eine Semi-Krylov-Projektionsmethode ist daher Q_{ik} ⊆ K_{ik} und L_{ik} ⊆ B_{ik}K_{ik}. In der Regel gilt dim Q_{ik}, dim L_{ik} ≤ ik. Existiert zusätzlich eine feste Matrix B mit L_{ik} = B^{T} Q_{ik}, so nennt man das Verfahren balanziert.

Semi-Krylov-Projektionsmethoden sind ökonomischer als ”volle” Krylov-Verfahren, besitzen aber keine Eigenschaften wie Endlichkeit oder Minimierung über dem ganzen Lösungsraum. Jedoch minimiert x_{i} die Fehler- bzw. Residuennorm über dem affinen Unterraum x_{i−1} +K_{k}(r_{i−1}, A). Da diese Optimierung in jedem Schritt über dem zugehörigen Krylov-Raum geschieht, erhofft man sich, daß die Norm des Fehlers bzw. des Residuums über dem ganzen Lösungsraum nicht anwächst. Über die Endlichkeit der Methoden kann im Allgemeinen auch keine Aussage gemacht werden, da die Dimension des aufspannenden Unterraums nicht notwendig wächst.

Im nächsten Abschnitt werden einige Beispiele von Semi-Krylov-Projektionsmethoden vorgestellt, die sich wiederum in der Konstruktion der Richtungsvektoren, die abhängig von Eigenschaften der Matrix A bzw. B sind, unterscheiden.

3.4 Beispiele und Algorithmen

Typische Beispiele für balanzierte Semi-Krylov-Projektionsmethoden sind sog. Truncated und Restarted Methoden. Das in Abschnitt 3.2 vorgestellte gcr-Verfahren liefert für eine definite Matrix A mit positiv definitem symmetrischen Anteil eine Lösung. Jedoch ist der Arbeits- und Speicheraufwand so erheblich, daß es praktisch nicht angewandt wird. Die erste Idee ist nun, die Summe zur Berechnung der p_{i} nach k Schritten abzubrechen (truncated gcr), d.h. man benutzt anstelle aller Richtungsvektoren nur die letzten k und erhält somit eine Menge von k + 1 A^{T} A- orthogonalen Richtungsvektoren {p_{i−k}, …, p_{i−1}, p_{i}}. Diesem Prinzip folgt das von Vinsome (siehe [Vins76]) vorgestellte Verfahren Orthomin(k) mit folgendem Algorithmus:

x_{i} = x_{i−1} + α_{i−1}p_{i−1}
mit α_{i−1} = \frac{(r_{i−1}, Ap_{i−1})}{(Ap_{i−1}, Ap_{i−1})}
r_{i} = b − Ax_{i}
p_{i} = r_{i} − ∑_{j=max(0,i−k)}^{i−1} β_{j}^{(i)}p_{j}

mit β_{j}^{(i)} = \frac{(Ar_{i}, Ap_{j})}{(Ap_{j}, Ap_{j})}

Dann gilt auch hier für eine Matrix A, deren symmetrischer Anteil positiv definit ist: Innerhalb des Iterationsschrittes sind die letzten k durch den Algorithmus erzeugten Richtungsvektoren A^{T} A-orthogonal, d.h. (Ap_{i}, Ap_{j}) = 0 i ≥ k, j = i − k, …, i − 1. Zwei aufeinanderfolgende Residuen sind A-orthogonal und x_{i} minimiert die Residuennorm über dem affinen Unterraum x_{i−k−1} + span{p_{i−k−1}, …, p_{i−1}}.

Ein zweiter Ansatz, das gcr-Verfahren durchfürbar zu machen, beruht darauf, den gcr-Algorithmus nach jeweils k Iterationen neu zu starten (restarted gcr, gcr(k-1)). Dabei konstruiert man sukzessive stets einen k-dimensionalen Unterraum, aus dem dann eine Näherungslösung bestimmt wird. Der Algorithmus berechnet Zwischennäherungen x̃_{j}, so daß für die x_{i−1} =: x̃_{0}, …, x̃_{k} =: x_{i} gilt, mit:

x̃_{j} = x̃_{j−1} + α̃_{j−1}p̃_{j−1}
mit α̃_{j−1} = \frac{(r̃_{j−1}, Ap̃_{j−1})}{(Ap̃_{j−1}, Ap̃_{j−1})}
r̃_{j} = r̃_{j−1} − α̃_{j−1}Ap̃_{j−1}
p̃_{j} = r̃_{j} − ∑_{l=0}^{j−1} β̃_{l}^{(j)}p̃_{l}
mit β̃_{l}^{(j)} = \frac{(Ar̃_{j}, Ap̃_{l})}{(Ap̃_{l}, Ap̃_{l})}.

Für eine Matrix A mit positiv definitem symmetrischen Anteil sind auch hier innerhalb jedes Iterationsschrittes die Richtungsvektoren p̃_{i} A^{T} A-orthogonal und die Residuen r̃_{i} sind A-orthogonal. Die Norm des Residuum wird über dem affinen Krylovraum x_{i−1} + K(r_{i−1}, A, k + 1) minimiert. Das gcr(k)-Verfahren ist eine polynomielle Methode k + 1-ten Grades, die mit B = A^{T} A definit in keinem Iterationsschritt abbricht.

Das Verfahren gmres(k), welches bei [SaSc86] genau analysiert wird, kann nun auf eine beliebige nichtsinguläre Matrix A angewandt werden. Im ersten Teil des gmres(k)- Algorithmus erzeugt man orthonormale Richtungsvektoren mit dem Gram-Schmidt Verfahren, die eine Basis für einen Krylov-Raum bilden.

p̃_{j} = Ap_{j−1} − ∑_{l=1}^{j−1} h_{lj}p_{l}
mit h_{lj} = (Ap_{j−1}, p_{l})
p_{j} = p̃_{j}/h_{j,j−1}
mit h_{j,j−1} = ∥p̃_{j}∥_{2}

Man gibt sich also k Suchrichtungen p_{j} vor und löst nun im zweiten Teil ein explizites Minimierungsproblems im IR^{k} mit der Methode der kleinsten Quadrate, um eine Näherungslösung x_{i} = x_{i−1} + ∑_{j=1}^{k} c_{j}p_{j} zu bestimmen.

gmres(k) ist ein polynomielles Iterationsverfahren k-ten Grades, welches in der Praxis benutzt wird, um lineare Gleichungssysteme mit beliebiger regulärer Koeffizientenmatrix zu lösen. Das Verfahren bricht nicht ab, aber es stagniert möglicherweise. Ist der symmetrische Anteil von A positiv definit, dann erzeugt gmres(k), bei exakter Arithmetik, die gleiche Lösung wie gcr(k-1).

4 Nicht polynomielle Projektionsmethoden und nicht projektive polynomielle Methoden

Schließlich sollen in diesem Kapitel noch Methoden erwähnt werden, die entweder nur projektiv bzw. nur polynomiell sind. Zu den nicht polynomiellen Projektionsmethoden gehören beispielsweise blockorientierte konjugierte Gradientenverfahren, sowie die in [SaSiYip88] untersuchten Methoden USYMLQ und USYMQR. Diese Verfahren sind definiert durch

x_{i} ∈ x_{0} + K̃_{i}(r_{0}, A^{T}), e_{i}⊥A^{T} K̃_{i}(r_{0}, A)

bzw.

x_{i} ∈ x_{0} + K̃_{i}(r_{0}, A^{T}), e_{i}⊥A^{T} AK̃_{i}(r_{0}, A^{T})

über dem sogenannten Quasi-Krylov-Raum

K̃_{i}(v, A) = span {v, Av, AA^{T} v, AA^{T} Av, …}.

Diese Methoden sind in der Regel nicht polynomiell, nach Definition zählen sie aber zu den Projektionsverfahren. Die USYMLQ-Methode kann jedoch trotzdem abbrechend sein, während das USYMQR-Verfahren balanziert ist mit einer Matrix B = A^{T} A, die für eine beliebige nichtsinguläre Matrix A symmetrisch und positiv definit ist, und somit die in Abschnitt 2.1 vorgestellten Eigenschaften besitzt.

Ein Beispiel für eine nicht projektive polynomielle Methode ist das quadratische konjugierte Gradientenverfahren (CGS), welches von [Sonn89] vorgestellt wurde. Ein durch diese Methode erzeugtes Residuum kann dargestellt werden als

r_{i} = (P_{i}(A))^{2}r_{0},

ist somit nicht mehr projektiv, erzeugt jedoch wiederum eine Näherungslösung über einem Krylov-Raum.

5 Zusammenfassung und Ausblick

In dieser Arbeit wurden iterative Methoden zur Lösung des linearen Gleichungssystems Ax = b untersucht. Eine Klassifizierung und Analyse wurde anhand spezieller Eigenschaften vorgenommen. So kann eine Aussage über Konvergenz der Iterationsverfahren, Minimierungseigenschaften sowie Fehlerschranken getroffen werden. Charakteristisch ist, daß Verfahren mit ”schönen” Eigenschaften nur auf eine eingeschränkte Klasse von Koeffizientenmatrizen A anwendbar ist. Die folgende Abbildung verdeutlicht die Zusammenhänge.

[figure] [figure] [figure] [figure] [figure] [figure] [figure] [figure] [figure] [figure] [figure] [figure] [figure] [figure] [figure] [figure] [figure] [figure] [figure] [figure]

Projektions-Methoden Polynomielle Methoden Hermitsche Krylov-Projektions-Methoden Konjugierte Gradienten Methoden Orthogonale Fehler- Methoden

Krylov-Projektions-Methoden

Semi-Krylov-Projektions-Methoden

Der abschließende Überblick liefert eine Auswahl typischer Verfahren, deren Verwendung von Eigenschaften der Matrix A abhängig ist.

Konjugierte Gradienten-Methoden cg-verfahren, cr-Verfahren

· Koeffizientenmatrix A: symmetrisch, positiv definit; · endliche Verfahren mit n ≤ N, eindeutig bestimmte Näherungslösung; Minimierung über Lösungsraum.

Orthogonale Fehlermethoden gcr-Verfahren

· A nicht notwendig symmetrisch, der symmetrische Anteil H(A) muß jedoch positiv definit sein; · allgemein kann keine Aussage über die Minimierungseigenschaften von Näherungslösungen gemacht werden.

Hermitsche/Symmetrische Krylov-Projektionsmethoden SYMMLQ-Verfahren

· A symmetrisch, nicht notwendig definit; · abbrechende verfahren, unbeschränkte Näherungslösungen.

Semi-Krylov-Projektionsmethoden orthomin(k)-Verfahren (truncated Verfahren), gcr(k-1)-Verfahren (restarted Verfahren)

· A besitzt positiv definiten symmetrischen Anteil; · Näherungslösungen mit Minimierungseigenschaften über einem k-dimensionalen Unterraum.

gmres(k) (restarted Verfahren)

· A beliebig nichtsingulär; · Näherungslösungen mit Minimierungseigenschaften über einem k-dimensionalen Unterraum; nicht abbrechendes Verfahren, möglicherweis stagnierend.

Literatur

[EiES83] S.C.Eisenstat, H.C.Elman und M.H.Schultz, Variational iterative methods for nonsymmetric systems of linear equations, SIAM Journal on Numerical Analysis, Vol.20, S.345-357, 1983.

[Elma82] H.C.Elman, Iterative Methods for Large, Sparse, Nonsymmetric Systems of Linear Equations, Yale University Doctoral Dissertation, New Haven, 1982.

[GrPfSch91] E.Groß, M.Pfister, U.Schendel, Projektionsmethoden zur Lösung linearer Gleichungssysteme, Preprint No. A-91-24, Fachbereich Mathematik, Freie Universität Berlin, 1991.

[Hest80] M.R.Hestenes, Conjugate Direction Methods in Optimization, Springer-Verlag, New York, 1980.

[JoMa90] W.D.Joubert, Th.A.Manteuffel, Iterative Methods for Nonsymmetric Linear Systems in... Academic Press, 1990.

[JoYo87] W.D.Joubert, D.M.Young, Necessary and sufficient conditions for the simplification of generalized conjugate gradient algorithms, Linear Algebra and its Applications, Vol. 88/89, S.449-485, 1987.

[KwSch93] B.Kwella, U.Schendel, Operator-Koeffizienten-Methoden, Preprint No. A-5/93, Serie A Mathematik, FU Berlin

[PaSa75] C.C.Paige, M.A.Saunders, Solution of sparse indefinite systems of linear equations, SIAM Journal on Numerical Analysis, Vol.12, S.617-629, 1975.

[SaSc86] Y.Saad und M.H.Schultz, GMRES: A Generalized Minimum Residual Algorithm for Solving Nonsymmetric Linear Systems, SIAM Journal on Scientific and Statistical Computing, Vol.7, S.856-869, 1986.

[Saad87] Y.Saad, Least Squares Polynomials in the Complex Plane and their Use for Solving Nonsymmetric Linear Systems, SIAM Journal on Numerical Analysis, Vol.24, S.155-169, 1987.

[SaSiYip88] M.A.Saunders, H.D.Simon, E.L.Yip, Two conjugate-gradient-type methods for sparse unsymmetric linear equations, SIAM Journal on Numerical Analysis, Vol.25, S.927-940, 1988.

[Sonn89] P.Sonneveld, CGS, a Fast Lanczos-type Solver for Nonsymmetric Linear Systems, SIAM J. Sci. Stat. Comput., Vol. 10, S.35-52, 1989.

[Vins76] P.K.W.Vinsome, Orthomin, an iterative method for solving sparse sets of simultaneous linear equations, in Proc. Fourth Symposium on Reservoir Simultation, Society of Petroleum Engineers of AIME, S.149-159, 1976.