Hier werden verschiedene Methoden behandelt, um eine Matrix
geeignet für die Verarbeitung vorzubereiten, und zwar
so, daß sich dabei möglichst Vorteile ergeben -
sozusagen Vereinfachungsmethoden für schwach besetzte
Matrizen.:)
Voraussetzungen
Einleitung, vor allem der Überblick über die verschiedenen Speichermethoden, da dort die einzelnen Begriffe und Matrixarten vorgestellt werden.
Vielleicht sollten wir an dieser Stelle erst einmal definieren, was wir überhaupt mit Vorkonditionierung meinen...:
Definition Eine Matrix M nennt man vorkonditionierend,
wenn sie die Transformation eines linearen Systems in ein äquivalentes
System (also mit gleicher Lösung) mit vorteilhafteren
Spektraleigenschaften bewirkt:
Für M bieten sich sinnvollerweise solche Matrizen an,
die A oder A-1 ähnelich sind (also
approximieren).
Bei der Verwendung dieser Methoden muß immer auch eine Abwägung
zwischen den Vorteilen (Steigerung der Konvergenzgeschwindigkeit) und den
Nachteilen (zusätzlicher Aufwand) getroffen werden.
Der Aufwand der meisten Methoden ist proportional
zur Variablenanzahl n. Der Aufwand vergrößert sich
also pro Iteration um einen konstanten multiplikativen Faktor n.
In der Praxis wird die Matrix M = M1M2
faktorisiert und das Gleichungssystem transformiert:
Wäre es nicht einfacher, statt dieser umständlichen Transformation
einfach M-1A zu verwenden ?
Die Matrizen M1 und M2
werden übrigens links- bzw. rechtskonditionierend genannt.
Das "Kochrezept" für die Vorkonditionierung
eines iterativen Verfahrens mit M1 und M2:
Bei dieser Methode (auch Punkt-Vorkonditionierung genannt)
besteht die Matrix M nur aus der Diagonale von A:
Theoretisch wäre hier keine zusätzliche Speicherung nötig,
aber in der Praxis werden die Reziprokwerte der Diagonaleinträge gespeichert,
um überflüssige, wiederholte Divisionen zu vermeiden.
Zerlegt man die Indexmenge J := {1,...,n} in paarweise
disjunkte Teilmengen Ji, dann erhält man eine
vorkonditionierende Matrix in Blockdiagonalform:
Diese Methode wird dann Jacobi-Block-Vorkonditionierung genannt.
Baut man die symmetrische Systemmatrix auf der Zerlegung
A = D + L + LT auf, dann ist die
SSOR-Matrix definiert als
oder parametrisiert als
Durch eine optimale Wahl von w (Omega) wird die Anzahl der Iterationsschritte
stark reduziert - in der Praxis ist diese "optimale" Wahl aber schwer zu treffen.
Dieses Verfahren hat gewisse Ähnlichkeiten mit dem folgenden, da hier
M und dort A in faktorisierter Form vorliegt.
Es ist übrigens immer möglich, auf diese Vorkonditionierung zu kommen,
da die Zerlegung von vornherein gegeben ist - sie kann also immer
angewendet werden.
Definition Eine Matrix heißt unvollständig,
wenn im Zerlegungsprozeß bestimme Auffüllelemente ignoriert werden.
Die Matrix M ist dann durch
gegeben, und ihre Wirksamkeit für die Konvergenzbeschleunigung ist von ihrer
"Approximationsqualität"
abhängig.
Die Berechnung einer unvollständigen Faktorisierung kann eventuell zu
einem Abbruch führen (Division durch Null als Pivotelement)
oder sie ergibt indefinite Matrizen (negatives Pivotelement).
Ihre Existenz ist aber in vielen Fällen gesichert, wenn die Ausgangsmatrix
A bestimmte Eigenschaften besitzt.
Bei iterativen Verfahren muß man immer auf den Rechenaufwand des
Faktorisierungsprozesses achten - dieser zahlt sich nur dann aus, wenn
das iterative Verfahren ohne Vorkonditionierung sehr viele
Schritte benötigen würde oder die gleiche Matrix M
für mehrere lineare Systeme verwendet werden kann.
(z.B. beim Newton-Verfahren für große nichtlineare Systeme mit schwach
besetzter Funktionalmatrix)
Hier geht man von einer Partitionierung der Systemmatrix aus.
Der wichtigste Unterschied zu den vorigen Methoden liegt in der
Inversion der Pivotblöcke - hier kommen zwei neue Probleme
dazu, die uns bisher erspart blieben:
Daraus folgt die Notwendigkeit zur...
Da die Berechnung einfach und schnell bleiben soll, fallen die Berechnung
der vollen Inversen und das Herausnehmen der entsprechenden Bandmatrix
einmal weg.
Die einfachste Approximation (Annäherung) der Inversen A-1
einer Bandmatrix A ist die Diagonalmatrix, die aus den
Reziprokwerten der Elemente der Hauptdiagonale von A besteht.
Anwendung vor allem bei partiellen Differentialgleichungen, weil dort
die Diagonalblöcke der Koeffizientenmatrix gewöhnlich stark diagonaldominant
sind.
Diese treten in der Praxis sehr häufig auf.
Wir haben also:
Daraus ergibt sich folgender Algorithmus zur Berechnung der approximativen
Pivotinversen:
Blockmethoden sind besonders interessant, weil sie sich sehr gut für Vektor-
und Parallelrechner eignen.
Wenn man ausgeht von:
kann man auf zwei Wegen zu einer unvollständigen Faktorisierung gelangen:
Die meisten der folgenden Methoden arbeiten auf diese Weise.
Dagegen wird die Anzahl der Iterationen durch die Vorkonditionierung
gewöhnlich nur um einen konstanten, von der Größe des
Systems unabhängigen Wert verbessert.
Sicherlich schon, aber diesem Ansatz fehlen zwei wichtige Eigenschaften,
die für den Erfolg mancher iterativer Methoden entscheidend
sind: Symmetrie und positive Definitheit.
Unsere Transformation M1-1AM2-1
weist diese Eigenschaften für M1 = M2T
auf.
mit der Koeffizientenmatrix M1-1AM2-1.
Vorteile:
Nachteile:
Auffüllelemente sind Nichtnullelemente der Faktormatrizen L und
U an Positionen, an denen A Nullen hat.Durchführung einer unvollständigen Faktorisierung
Unvollständige Blockfaktorisierung
Es wird mit A eine normale unvollständige Faktorisierung duchgeführt,
aber es werden Teilblöcke der Matrix als Basiseinheiten verwendet.
Approximation der Inversen
(Für andere Möglichleiten schlag' nach bei Axelsson, Eijkhout sowie Axelsson, Polman.)
Block-Tridiagonalmatrizen
Sie erfüllt alle gewünschten Eigenschaften, da außerhalb der Diagonalblöcke Aii
keine Auffüllung auftreten kann.
Parallele unvollständige Blockfaktorisierung
(D ist dabei die Block-Diagonalmatrix aus den Pivotblöcken)
:
Die Lösung des Gesamtsystems erreicht man hier über die Lösung kleinerer Teilsysteme mit den Xi-Matrizen.
Hier ergeben sich durch die dominante Multiplikation mit den Yi-Blöcken, was große Vorteile auf Vektorrechnern bringt.
Unvollständige LQ-Faktorisierung
Diese Methode ist eine Alternative zur unvollständigen LU-Faktorisierung und ebenso zur Zerlegung schwach besetzter Matrizen geeignet.
Ausgangspunkt ist hier die Orthogonalisierung der Zeilen mit dem
Schmidtschen Verfahren.
Mit dem Faktor L der unvollständigen LQ-Zerlegung der Matrix A
erhält man schließlich die unvollständige Cholesky-Zerlegung
von AAT.
Aus der Praxis weiß man, daß diese Methode bei einigen relevanten Problemen zufriedenstellende Effizienz ermöglicht.
Polynomiale Vorkonditionierung
Wie in der Einleitung erwähnt, gibt es auch vorkonditionierende
Matrizen, mit denen A-1 approximiert wird.
Zu diesen gehören auch die polynomialen Ansätze.
Nehmen wir an, A kann in der Form A = I - B dargestellt
werden, wobei der Spektralradius von B kleiner als 1 ist.
Durch Anwendung der Neumann-Reihe kann man die Inverse von A
als unendliche Reihe anschreiben:
Dadurch kann man eine Näherung für A-1 durch Abbrechen dieser Reihe erreichen.
Beispiel: Dubois, Greenbaum und Rodrigue untersuchten das Verhältnuis zwischen einer einfachen Vorkonditionierung mittels der Zerlegung A = M - N und einer plynomial vorkonditionierten Methode mit:
Dabei ergab sich hauptsächlich, daß für "klassische"
Iterationsverfahren k Schritte der polynomialen Vorkonditionierung
exakt äquivalent zu kp Schritten der ursprünglichen Methode
sind.
Für beschleunigte Verfahren, speziell für die Tschebyscheff- Iteration,
kann man also durch die Vorkonditionierung die Anzahl der Iterationen maximal
um den Faktor p verringern.
Die Anzahl der Matrix-Vektor-Produkte kann aber nicht verringert werden.
Da aber ein großer Teil der inneren Produkte und der Aktualisierungsoperationen eliminiert wird, kann durch das zuletzt genannte Verfahren trotzdem oft eine große Effizienzsteigerung erreicht werden.
Noch ein Beispiel: Man kann sich M auch als ein
Matrixpolynom M = Pn(A) mit der Normalisierung P(0) = 1
vorstellen. Das günstigste Polynom ist dann jenes, das || I - M-1A ||
minimiert.
Bei Verwendung der Maximumnorm ergeben sich dadurch Tschebyscheff- Polynome.
Die erforderliche Spektrumsabschätzung kann man dann leicht aus einer CG-Iteration
berechnen.