Sherman-Morrison-Woodbury-Formel ihre Gültigkeit und (259) ist daher tunlichst als Al-gorithmus zur Lösung des gestörten Gleichungssystemes
(A+V SWT)y=b zu lesen, wenn ein Löser „A−1“ für das ungestörte System
Ax=b zur Verfügung steht177
Eine solche algorithmische Lesart der Sherman-Morrison-Woodbury-Formel könnte etwa wie folgt aussehen
y= z }| {x
A−1b− z }| {X
A−1V (S−1+WT z }| {X
A−1V
| {z }
Z
)−1
z }| {d
WT z }| {x
A−1b
| {z }
r
und zur folgenden Rechenvorschrift führen
Algorithmus 4.13 (Sherman-Morrison-Woodbury-Algorithmus; Variante 1) Ax=b; ein Gleichungssystem lösen!
AX =V; k Gleichungssysteme lösen!
Z :=WTX;
d:=WTx;
(S−1 +Z)r =d; ein Gleichungssystem lösen!
y:=x−Xr.
Andere Varianten, die von speziellen Gegebenheiten Gebrauch machen (vgl. das einleitende Beispiel) sind denkbar. Liegt fürA z.B. eineLR-Zerlegung vor, so kann die Sherman-Mor-rison-Formel z.B. wie folgt verwirklicht werden:
(A=LR vorgelegt)
Algorithmus 4.14 (Sherman-Morrison-Woodbury-Algorithmus, Variante 2) Ax=b; ein Gleichungssystem lösen!
LX(1) =V; k Dreieckssysteme lösen!
RTX(2) =W; k Dreieckssysteme lösen!
Z :=X(2)TX(1); d:=WTx;
(S−1 +Z)r =d; ein Gleichungssystem lösen!
c:=X(1)r;
Rw=c; ein Dreieckssystem lösen!
y:=x−w;
177Die Penetranz, mit der hier wieder auf dem Verbot der inversen Matrix herumgetrampelt wird, erklärt sich aus einer ebensolchen Penetranz, mit der eine doch recht große Anzahl (sonst eigentlich ganz intelligent erscheinender) Anwender numerischer Verfahren (und das schließt auch (sogenannte) Mathematikerleider
5 Iterative Verfahren zur Lösung großer Linearer Syste-me
Iterative Verfahren zur Lösung linearer Gleichungssysteme
Ax=b, A∈R(n,n), b ∈Rn gegeben, x∈Rn gesucht. (261) hatten wir schon im Abschnitt 3.1.3 als Anwendung des Banachschen Fixpunktsatzes dis-kutiert. Diese Iterationen waren speziell Iterationen der Form
xn+1 = Φ(xn), n ∈N
gewesen, bei denen der nächste Iterationsvektor mit einer festen Verfahrensfunktion Φ : Rn −→ Rn nur aus dem aktuellen Iterationsvektor xn berechnet wird, und bei dem die Lösungen von (261) Fixpunkte vonΦsind. Diese Verfahren nennt man „stationäre Verfah-ren‘ “.
Es gibt auch „nichtstationäre Iterationen“. Hier kann sich die Iterationsfunktion von Schritt zu Schritt ändern
xn+1 = Φn(xn), n ∈N.
Es kann die Iteration hier etwa nicht nur vom letzten Iterationswert abhängen sondern von mehreren. So kann die Iteration aus ihrem konkreten Verlauf „lernen“.
Wir können in diesem Abschnitt nur einige weitere Ideen der Verfahren vorstellen, denn mit der Schilderung des aktuellen Forschungsstandes zu iterativen Verfahren178 für lineare Gleichungssysteme ließen sich ohne Probleme fünf und mehr schöne dicke Bücher füllen, wie z.B. [DMY], [BEA],[YS], [AG], [OA], [HAVDV].
Um die Frage, wann solche Iterationsverfahren direkten Methoden vorzuziehen sind, haben wir uns oben elegant herumgedrückt, indem wir sie lieber gar nicht erst aufgeworfen haben.
So verfahren auch die meisten anderen Autoren. Einige der „Verfechter iterativer Methoden“
versteigen sich zwar zu der Behauptung, iterative Methoden seien stets anzuwenden, wenn die Systemmatrix groß und dünnbesetzt sei. Dies ist so aber nicht haltbar.
Leider kann man dem Anwender der numerischen Mathematik wohl (noch?) keine einfache Anleitung an die Hand geben.
Überhuber fasst in seinem Buch [UEB] seine Erfahrungen in der folgenden Tabelle zusam-men.
direkte Verfahren iterative Verfahren Genauigkeit nicht beeinflussbar wählbar
Rechenaufwand vorhersagbar meist nicht vorhersagbar aber oft kleiner neue rechte Seite rasch keine Zeitersparnis
Speicherbedarf größer kleiner
Startwert-Vorgabe nicht erforderlich meist vorteilhaft Algorithmus Parameter nicht erforderlich müssen gesetzt werden Black-box-Verwendung möglich oft nicht möglich
Rubustheit ja nein
Bethke und Voss (vgl [BV]) haben eine MATLAB-Umgebung bereitgestellt, mit der der potentielle Benutzer erstens einsehen kann, dass ein klare (überall geltende) Aussage nicht machbar ist, mit der er aber auch testen kann, welche der bis 2003 verfügbaren Verfahren
178Die Forschung kann hier auf gar keinen Fall als irgendwie abgeschlossen bezeichnet werden. Die Lei-stung iterativer Verfahren ist bislang noch alles andere als wirklich zufriedenstellend.
für seine Problemklasse angemessen erscheint.
Dabei können (siehe [BV]) manchmal schon überraschende Ergebnisse herauskommen.
Wenn man die Ausführungen von Überhuber noch ein wenig ergänzen und vertiefen will, so sind vermutlich die folgenden Aussagen als normative Ideen geeignet.
Iterative Verfahrungen werden angewandt, wenn
1. eine direkte Lösung aus Speicherplatzgründen nicht durchführbar ist, 2. die Matrix eine für eine iterative Lösung besonders geeignete Struktur hat, 3. eine (sehr) gute Näherung für die Lösung existiert
4. ein direkter Löser für eine in einer Norm von der aktuellen Systemmatrix nur leicht abweichende Matrix vorhanden ist
5. iterative Löser schneller sind als direkte Löser
Hierzu sind zunächst einige allgemeine Kommentare angebracht
Ad 1. Es kann natürlich leicht vorkommen, dass trotz raffinierter Speichertechniken die Elemente der LR-Zerlegung einer Matrix nicht mehr speicherbar sind, die Multipli-kation eines Vektors mit der Systemmatrix aber ganz unproblematisch ist179. Wenn dann direkte Verfahren, die Teile der Matrix und ihrer Zerlegungen auf externe Zwi-schenspeicher legen, auch nicht mehr anwendbar sind, wird man alternativ versuchen, mit iterativen Verfahren Erfolg zu haben.
Wenn man sich auf iterative Verfahren einlässt, weiß man i.a nicht, welches der zig existierenden Verfahren man nehmen soll, man weiß i.a. nicht, wie schnell es konver-giert, ja, man weiß i.a. noch nicht einmal, ob es überhaupt konvergiert.
Ad 2. Man hat heute schon eine Reihe von Anwendungsgebieten identifiziert, deren Glei-chungssysteme für iterative Verfahren gut geeignet sind. Insbesondere sind dies die elliptischen partiellen Differentialgleichungen, deren Diskretisierungen auf Matrizen führen, die mit speziell angepassten Iterationsverfahren viel schneller lösbar sind als mit direkten Lösern. Die Matrix A von Seite 29 braucht man bei iterativer Anwen-dung z.B. überhaupt nicht zu speichern, weil die Multiplikation mit A vollständig mit dem Differenzenstern (51) zu bewältigen ist. Außerdem werden hier problemspe-zifische Iterationsverfahren (siehe unten z.B. Mehrgitter-Verfahren) eingesetzt.
Ad 3. Wenn schon eine sehr gute Näherung für die Lösung bekannt ist180, liegt es na-he, die „kleinen Fehler“, die die Näherung noch hat, mit einem Iterationsverfahren auszumerzen. Da direkte Verfahren keine Lösungsnäherung arbeitssparend einsetzen können, sind sie in solch einer Situation erst einmal im Nachteil.
Ad 4. Wenn ein Löser A˜−1 für eine Matrix A˜ vorhanden ist181 , die A recht nahe kommt
∥A−A˜∥ ≤ε ( klein ).
Dann ist
xn+1 :=xn−A˜−1(Axn−b)
179Man wird sehen, dass die meisten heute verwendeten iterativen Verfahren die Systemmatrix im We-sentlichen nur durch ihre multiplikative Anwendung auf einen Vektor in’s Spiel bringen.
180Gelegenheiten hierfür sprechen wir noch an.
− −
eine Iteration, deren Iterationsmatrix182
M :=I−A˜−1A= ˜A−1(A−A)˜ wegen
∥M∥ ≤ ∥A˜−1∥ · ∥A−A˜∥ ≤ ∥A˜−1∥ ·ε bei moderat großem ∥A˜−1∥ einen kleinen Spektralradius hat.
Ad 5. Hierzu ist nichts zu sagen. Wenn man diese Information hat, verwendet man na-türlich iterative Verfahren.
Eine Situation, die die unter 3. und 4. beschriebenen Gegebenheiten vereint, findet man bei der sogenannten „Nachiteration‘ “.
Algorithmus 5.1 (Nachiteration) Ein lineares Gleichungssystem
Ax=b
werde durch eine Gauss-Elimination gbearbeitet, die am Ende eine recht ordentliche Lö-sungsapproximation x0 liefert sowie eine approximative LR-Zerlegung
A ≈LR
mit normierter unterer Dreiecksmatrix L und rechter oberer Dreiecksmatrix R.
Bei exakter Rechnung wäre x0 die Lösung und A das Produkt von L und R. Tatsächlich ist bei der Ausführung auf dem Rechner zu erwarten, dass x0 ̸=x und LR ̸= A. Es wird aber andererseits auch zu erwarten sein, dassx0 eine gute Näherung der Lösungx ist und der Lösungsprozess R−1L−1 eine gute Approximation von A−1. MitA˜−1 =R−1L−1 lautet die Nachiteration also
xn+1 :=xn−R−1L−1(Ax0−b), n≥0.
Zu sehr guten Ergebnisse gelangt man, wenn man das sogenannte „Residuum“ Ax0−b mit erhöhter Rechnergenauigkeit berechnet.
Beispiel 5.2 (Nachiteration)
Führt man den Gauss-Algorithmus (mit LR-Zerlegung) für das System Ax=b für
A=
2 −1 0 0 0 0 0
−1 2 −1 0 0 0 0
0 −1 2 −1 0 0 0
0 0 −1 2 −1 0 0
0 0 0 −1 2 −1 0
0 0 0 0 −1 2 −1
0 0 0 0 0 −1 2
und b= (1,1,1,1,1,1,1)T
mit zweistelliger Dezimalgleitpunktarithmetik aus, so ergibt sich als Approximation für die korrekte Lösung
x= (3.5000,6.0000,7.5000,8.0000,7.5000,6.0000,3.5000)T
182Vgl. Seite 43
die ziemlich falsche Näherung
x0 = (4.0000,7.0000,8.0000,8.0000,7.0000,5.0000,2.5000)T mit den approximativen L- und R-Matrizen
Lˆ =
1.0000 0 0 0 0 0 0
−0.5000 1.0000 0 0 0 0 0
0 −0.7000 1.0000 0 0 0 0
0 0 −0.8000 1.0000 0 0 0
0 0 0 −0.8000 1.0000 0 0
0 0 0 0 −0.8000 1.0000 0
0 0 0 0 0 −0.8000 1.0000
und
Rˆ =
2.0000 −1.0000 0 0 0 0 0
0 1.5000 −1.0000 0 0 0 0
0 0 1.3000 −1.0000 0 0 0
0 0 0 1.2000 −1.0000 0 0
0 0 0 0 1.2000 −1.0000 0
0 0 0 0 0 1.2000 −1.0000
0 0 0 0 0 0 1.2000
.
Ganz offensichtlich ist
Lˆ·Rˆ =
2.0000 −1.0000 0 0 0 0 0
−1.0000 2.0000 −1.0000 0 0 0 0
0 −1.0500 2.0000 −1.0000 0 0 0
0 0 −1.0400 2.0000 −1.0000 0 0
0 0 0 −0.9600 2.0000 −1.0000 0
0 0 0 0 −0.9600 2.0000 −1.0000
0 0 0 0 0 −0.9600 2.0000
nicht gleichA. Führen wir mit diesen fehlerhaften MatrizenLˆundRˆ nun die Nachiteration aus
xk+1 =xk−Rˆ−1Lˆ−1(Axk−b), k ≥0, so ergeben sich folgende Iterationsvektoren
x1 =
3.4464 5.8928 7.3393 7.8410 7.3693 5.8911 3.4277
, x2 =
3.5027 6.0053 7.5080 8.0050 7.4953 5.9922 3.4941
, x3 =
3.4996 5.9993 7.4989 7.9989 7.4992 5.9993 3.4995
mit konsekutiven Fehlern E(k) = ∥xk−x∥2 in der 2-Norm
E = [1.936,0.3156,0.1559 10−1,0.2080 10−2,0.1013 10−3,0.1375 10−4,0.6654 10−6,0.9051 10−7, . . . und einer Folge von Verbesserungsfaktoren Vn+1; =En+1/En