3.8 Lineare Gleichungssysteme

Lineare Gleichungssysteme

entstehen in allen Bereichen der Wissenschaft, entweder direkt zur Formulierung eines praktischen Problems oder als Teilschritt innerhalb komplexerer Lösungsstrategien. Dabei sind die Zahlen $a_{i,j}$ und $b_i$ gegeben und die Zahlen $x_1,\ldots,x_n$ gesucht.

Zur Vereinfachung der Darstellung setzt man üblicherweise die Systemmatrix

und den Vektor der rechten Seiten

Das Gleichungssystem bekommt damit die Form

mit dem Lösungsvektor

Es existieren vielfältige numerische Lösungsverfahren mit unterschiedlichen Eigenschaften bzgl. Stabilität, Rechenaufwand, Speicherbedarf und Anwendungsbreite. Schauen uns hier zwei relativ breit anwendbare Verfahren näher an und geben einen kurzen Überblick über weitere Verfahrensideen.

Beschränken uns auf invertierbare quadratische Systemmatrizen, also auf den Fall $m=n$. Insbesondere folgt aus der Forderung der Invertierbarkeit, dass das Gleichungssystem genau eine Lösung besitzen soll. Diese kann dann als

geschrieben werden.

3.8.1Kondition¶

Zur Beurteilung der Kondition beim Lösen eines linearen Gleichungssystems muss zunächst geklärt werden, welche Größen als (potentiell fehlerbehaftete) Eingabe zu betrachten sind. Hier gibt es zwei Ansätze:

• Nur $b$ wird als Eingabe betrachtet und es ist die Kondition der Abbildung $b\mapsto A^{-1}\,b$ gesucht. Dieser Ansatz ist leichter zu handhaben und wird immer verwendet, wenn $A$ sehr genau oder sogar exakt bekannt ist.

• $A$ und $b$ werden beide als Eingaben betrachtet und es ist die Kondition der Abbildung $(A,b)\mapsto A^{-1}\,b$ gesucht.

Weiterhin bestehen gewisse Freiheiten bei der Beantwortung der Frage, wie wir Abweichungen zwischen exakten und fehlerbehafteten Ein- bzw. Ausgaben ausdrücken. Bisher haben wir Konditionsbetrachtungen nur für Abbildungen $f:\bbR^d\to\bbR$ angestellt. Auf der Eingabeseite hatten wir uns für die Summe der komponentenweisen relativen Fehler entschieden. Auf der Ausgabeseite bestand gar keine Wahl. Im Kontext linearer Gleichungssysteme muss nun auch auf der Ausgabeseite entschieden werden, wie die relativen Fehler in den einzelnen Komponenten zu einem Gesamtausgabefehler zusammengeführt werden.

3.8.1.1Vektornormen¶

Für einen Vektor $v\in\bbR^n$, z.B. den Vektor der komponentenweisen relativen Ein- oder Ausgabefehler, können wir verschiedene Normen einführen, also Zahlen $\|v\|$, die die Größe der einzelnen Einträge des Vektors sinnvoll zusammenfassen. Unter “sinnvoll” verstehen wir, dass folgende Eigenschaften gelten sollen:

• $\|v\|\geq 0\quad$ für alle $v\in\bbR^n$,

• $\|v\|=0\;\Rightarrow\; v=0$,

• $\|\alpha\,v\|=|\alpha|\,\|v\|\quad$ für alle $\alpha\in\bbR$ und alle $v\in\bbR^n$,

• $\|v+w\|\leq\|v\|+\|w\|\quad$ für alle $v,w\in\bbR^n$.

Beispiele:

• 1-Norm:

  $$\|v\|_1:=\sum_{i=1}^n|v_i|,$$

  (7)

• $\infty$-Norm:

  $$\|v\|_\infty:=\max_{i\in\{1,\ldots,n\}}|v_i|,$$

  (8)

• euklidische Norm:

  $$\|v\|_2:=\sqrt{\sum_{i=1}^n v_i^2}.$$

  (9)

Alle Normen in $\bbR^n$ sind äquivalent, d.h. für zwei Normen $\|\cdot\|_\alpha$ und $\|\cdot\|_\beta$ gibt es stets Konstanten $c,C>0$, sodass

gilt. Mit Blick auf Fehlerabschätzungen beeinflussen die eingesetzten Normen also höchstens die auftretenden konstanten Faktoren, die ohnehin kaum relevant für die Kernaussagen sind.

3.8.1.2Matrixnormen¶

Bei Konditionsuntersuchungen werden wir stets auf Vektoren der Form $A\,x$ und $A^{-1}\,b$ stoßen, die wir in Beziehung zu $x$ bzw. $b$ setzen möchten. Zu diesem Zweck haben sich von einer Vektornorm induzierte Matrixnormen etabliert.

Am häufigsten anzutreffen ist die Spektralnorm, die entsteht, wenn für beide Vektornormen die euklidische Norm genutzt wird:

Die Spektralnorm ist von der gelegentlich auch verwendeten Frobenius-Norm

zu unterscheiden, welche nicht durch Vektornormen induziert wird, aber dennoch die vier Normbedingungen erfüllt.

Die von der 1-Norm oder der $\infty$-Norm induzierten Matrixnormen spielen nur selten eine Rolle, sind aber deutlich einfacher zu berechnen als die Spektralnorm. Für die 1-Norm kann man

zeigen. Für die $\infty$-Norm erhält man

Analog zu den Vektornormen sind auch alle Matrixnormen zu einander äquivalent.

3.8.1.3Fehlerhafte rechte Seite (Konditionszahl)¶

Betrachten wir nur $b$ als Eingabe, so erhalten wir für den relativen Ausgabefehler bei fehlerhafteteter Eingabe $\tilde{b}$ die Abschätzung

wobei die verwendeten Matrix- und Vektornormen kompatibel zu einander sein sollen, ansonsten aber beliebig gewählt werden können (IDV 380).

Beachte, dass beim bisherigen Konditionsbegriff für Abbildungen nach $\bbR$ (statt $\bbR^n$) der Eingabefehler gerade die 1-Norm des Vektors der komponentenweisen relativen Fehler war. Im Kontext linearer Gleichungssysteme ist jedoch die Norm des Vektors der absoluten Fehler geteilt durch die Norm des exakten Eingabevektors als Ausdruck des relativen Eingabefehlers besser handhabbar. Auf der Ausgabeseite wird die gleiche Variante zur Beschreibung des relativen Ausgabefehlers verwendet.

Je nach Wahl der Matrixnorm entstehen unterschiedliche Zahlen $\kappa_1(A)$, $\kappa_2(A)$, $\kappa_\infty(A)$. Wenn im Folgenden kein Index angegeben ist, gelten die Resulte für alle Matrixnormen.

Man kann zeigen, dass stets

gilt.

Für die euklidische Matrixnorm kann man noch

zeigen, wobei $\lambda_{\mathrm{max}}$ und $\lambda_{\mathrm{min}}$ der größte und der kleinste Eigenwert von $A$ sind.

3.8.1.4Fehlerhafte Matrix und rechte Seite¶

Untersuchen noch die Kondition der Abbildung $(A,b)\mapsto A^{-1}\,b$ bei fehlerbehafteter Matrix $\tilde{A}$ und fehlerbehafteter rechter Seite $\tilde{b}$. Hier muss zunächst geklärt werden, unter welchen Bedingungen $\tilde{A}$ überhaupt noch invertierbar ist.

Für den Fehlerzusammenhang zwischen Eingabe und Ausgabe kann man dann

zeigen, wobei Matrix- und Vektornormen wieder kompatibel zu einander sein sollen. Bei kleinem Fehler in $A$ und nicht zu großer Kondition $\kappa(A)$ beschreibt also auch hier $\kappa(A)$ die Fehlerverstärkung.

3.8.1.5Beispiele¶

3.8.1.5.1Einheitsmatrix¶

Für

gilt $A=A^{-1}$ und $\|A\|=1=\|A^{-1}\|$, also $\kappa(A)=1$.

3.8.1.5.2Hilbert-Matrix¶

Für

wächst $\kappa(A)$ exponentiell mit $n$, siehe The Condition of the Finite Segments of the Hilbert Matrix.

3.8.2Wie es nicht geht¶

Für $A\in\bbR^{n\times n}$ kann die Lösung $x$ des Gleichungssystems mit der Cramer’schen Regel ermittelt werden: Zur Berechnung von $x_i$ sei $A_{i}$ wie $A$, aber mit Spalte $i$ ersetzt durch die rechte Seite $b$. Dann gilt

wobei

die Determinate von $A$ ist (analog für $A_i$). Die Summe läuft hier über alle Permutationen $\sigma:\{1,\ldots,n\}\to\{1,\ldots,n\}$, also alle möglichen Anordnungen der Zahlen $1,\ldots,n$. Der Wert $\sgn\sigma$ ist das Vorzeichen der Permutation und ist entweder 1 oder -1.

Der Aufwand für das Berechnen von Determinanten ist enorm! Er liegt bei $n!$ Additionen und $(n-1)\,n!$ Multiplikationen pro Determinante. Für das Lösen eines Gleichungssystems werden $n+1$ Determinanten benötigt. Für $n=20$ würde ein 1-Peta-FLOPS-Rechner mehr als 11 Tage für das Lösen mittels Cramer’scher Regel benötigen.

Zusätzlich sind auf der Cramer’schen Regel beruhende Algorithmen durch die vielen Differenzen als instabil anzusehen.

Auch das explizite Berechnen der Inverse $A^{-1}$, um dann $x=A^{-1}\,b$ zu bekommen, ist im Allgemeinen aufwendiger als das Lösen des Gleichungssystems ohne Verwendung der Inverse. Der Rechenaufwand für den unten behandelten Gauß-Algorithmus entspricht zwar dem des Invertierens (etwa $n^3$ Grundoperationen), der Gauß-Algorithmus ist aber stabiler. Praktisch auftretende Gleichungssysteme haben meist beim Lösen vorteilhaft nutzbare Zusatzeigenschaften (symmetrisch, dünn besetzt, Bandstruktur,...), die den Einsatz schnellerer und speichersparenderer Algorithmen zulassen. So ist beispielsweise die im Allgemeinen vollbesetzte Inverse von großen Tridiagonalmatrizen, wie sie häufig in der Praxis auftreten, viel zu groß um sie überhaupt im Arbeitsspeicher ablegen zu können.

Grundsätzlich gilt: Inverse verweiden. Ausnahmen sind sehr kleine, gut konditionierte Gleichungssysteme, die für viele verschiedene rechte Seiten zu lösen sind (z.B. Koordinatentransformationen in der Computergrafik).

3.8.3Gauß-Algorithmus¶

3.8.3.1Idee¶

Folgende Beobachtung lässt sich zum Lösen linearer Gleichungssysteme einsetzen:

Die LR-Zerlegung von $A$ kann leicht algorithmisch bestimmt werden (verzichten hier auf die Details). Liegt diese vor, so erfolgt das Lösen des Gleichungssystems durch Vorwärtssubstitution und anschließende Rückwärtssubstitution:

3.8.3.2Stabilität¶

Sei $\tilde{x}$ die mittels Gauß-Algorithmus aus exakten Eingaben $A$ und $b$ numerisch ermittelte und somit (rundungs-)fehlerbehaftete Lösung. Liegt der relative Fehler zwischen $b$ und $A\,\tilde{x}$ im Bereich des erwarteten Eingabefehlers, so ist der Algorithmus stabil. Die Differenz $b-A\tilde{x}$ heißt auch Residuum.

Wie groß das Residuum ist, ist allgemein kaum zu bestimmen. Bei schlecht konditioniertem $A$ wird es jedenfalls groß sein, da die Rundungsfehler aus den ersten Rechenschritten analog zu Eingabefehlern verstärkt werden. Der Gauss-Algorithmus ist somit nicht als stabil anzusehen.

Folgende Punkte fassen die Sachlage bzgl. Stabilität zusammen:

• Theoretisch ist der Gauss-Algorithmus instabil.

• Für gewisse Problemklassen (z.B. Tridiagonalmatrizen) ist er garantiert stabil.

• Praktisch tritt Instabilität nur sehr selten auf, sofern der einfache Gauß-Algorithmus durch eine sogenannte Pivotsuche verbessert wird.

• Für nicht zu große Gleichungssysteme ohne spezielle Struktur ist der Gauß-Algorithmus heute das Standardlösungsverfahren.

3.8.3.3Zeit- und Speicherbedarf¶

Rechenaufwand für die LR-Zerlegung: ca. $n^3$ Grundoperationen.

Rechenaufwand für Vorwärts- und Rückwärtssubstitution: ca. $n^2$ Grundoperationen.

Der Rechenaufwand ist für große Gleichungssysteme also relativ hoch. Ist ein System für mehrere rechte Seiten zu lösen, so muss die LR-Zerlegung allerdings nur einmal berechnet werden. Das Lösen der weiteren Systeme benötigt dann nur noch $n^2$ Operationen.

Bei geschickter Implementierung benötigt der Algorithmus im Wesentlichen nur den Speicher für das Ablegen der Matrix $A$, also Größenordnung $n^2$.

3.8.4Cholesky-Verfahren¶

Für positiv definite symmetrische Matrizen $A\in\bbR^{n\times n}$ kann die LR-Zerlegung einfacher ausgedrückt und effizienter berechnet werden.

3.8.4.1Positiv definite Matrizen¶

Eine symmetrische Matrix heißt positiv definit, wenn gilt:

Positiv definite symmetrische Matrizen haben viele vorteilhafte Eigenschaften (siehe Grundlagenliteratur zur linearen Algebra):

• $A$ ist stets invertierbar und $A^{-1}$ ist wieder symmetrisch und positiv definit.

• Alle Eigenwerte von $A$ sind reell und positiv.

• $\det A>0$.

• Alle Diagonaleinträge sind positiv.

• Der betragsgrößte Eintrag von $A$ liegt auf der Diagonale.

Positiv definite symmetrische Matrizen treten in zahlreichen Anwendungen auf. Beispielsweise kann man leicht zeigen, dass Matrizen der Form $A=B^\rmT\,B$ für invertierbares $B\in\bbR^{n,n}$ stets symmetrisch und positiv definit sind. Solche Matrizen treten unter anderem beim Berechnen von Skalarprodukten in transformierten Koordinaten auf. $B$ ist dabei die Koordinatentransformation.

3.8.4.2Cholesky-Zerlegung¶

Man folgende spezielle Form der LR-Zerlegung herleiten:

Diese sogenannte Cholesky-Zerlegung erhält man direkt aus der LR-Zerlegung, wenn man die Symmetrie $A=A^\rmT$ ausnutzt. Insbesondere gilt $R=D\,L^\rmT$.

Allerdings geht es auf direktem Wege etwas schneller (verzichten hier auf die Details).

3.8.4.3Zeit- und Speicherbedarf¶

Der Rechenaufwand liegt wie bei der LR-Zerlegung in der Größenordnung $n^3$, ist aber trotzdem nur etwa eine Viertel so groß.

Der Speicherbedarf ist etwa doppelt so hoch wie beim Gauß-Algorithmus, da gewissen Zwischenergebnisse zusätzlich zur Matrix $A$ im Speicher abgelegt werden müssen.

3.8.4.4Gesamtalgorithmus¶

Ist die Cholesky-Zerlegung bekannt, kann das Gleichungssystem $A\,x=b$ analog zur LR-Zerlegung mittels Vorwärts- und anschließender Rückwärtssubstitution gelöst werden:

3.8.4.5Stabilität¶

Man kann zeigen, dass das Cholesky-Verfahren rückwärtsstabil ist (also auch vorwärtsstabil). Insbesondere gilt das Cholesky-Verfahren als stabiler als der Gauss-Algorithmus.

3.8.5Weitere Verfahren¶

Wir erwähnen kurz weitere Verfahrensideen um bei Bedarf mit den Begrifflichkeiten vertraut zu sein.

3.8.5.1QR-Zerlegung¶

Die LR-Zerlegung steht nur für quadratische Matrizen zur Verfügung. Für allgemeine rechteckige Matrizen $A\in\bbR^{m\times n}$, $m\geq n$ (überbestimmtes Gleichungssystem) kann man jedoch stets die sogenannte QR-Zerlegung $A=Q\,R$ finden. Dabei ist $R\in\bbR^{m\times n}$ wieder eine obere Dreiecksmatrix, wobei die unteren $m-n$ Zeilen mit Nullen aufgefüllt werden, und $Q\in\bbR^{m\times m}$ eine orthogonale Matrix. Die Spalten von $Q$ sind also orthogonal zu einander und haben Norm 1. Für orthogonale Matrizen gilt stets $Q^{-1}=Q^\rmT$, sodass sich das Lösen des Gleichungssystems $A\,x=b$ auf Rückwärtseinsetzen und Multiplikation mit $Q^\rmT$ reduziert.

Für quadratische Matrizen ($m=n$) kann man die QR-Zerlegung mit etwa $n^3$ Grundoperationen berechnen.

3.8.5.2Iterative Verfahren¶

Iterative Verfahren liefern im Gegensatz zum Gauß-Algorithmus und zum Cholesky-Verfahren nicht die (bis auf Rundungsfehler) exakte Lösung des Gleichungssystems $A\,x=b$, sondern nur Näherungslösungen. Die Verfahren berechnen zunächst eine sehr grobe Näherung, die dann Schritt für Schritt verfeinert wird. Dieser Ansatz erlaubt deutlich geringere Rechenzeiten bei sehr großen Gleichungssystemen.

3.8.5.2.1Richardson-Iteration¶

Für die Lösung $x$ des Gleichungssystems und eine beliebige reelle Zahle $\tau>0$ gilt

Dies ist eine sogenannte Fixpunktgleichung, da das Auswerten der rechten Seite wieder auf $x$ führt. Daraus kann man die Iterationsvorschrift

ableiten. Man kann zeigen, dass für hinreichend kleines $\tau$ die Folge $x_0:=0,x_1,x_2,\ldots$ gegen die Lösung des Gleichungssystems konvergiert. Je nach gewünschter Genauigkeit kann die Iteration früher oder später abgebrochen werden.

Die Iterationsvorschrift kann zu

verallgemeinert werden, wobei $M\in\bbR^{n\times n}$ und $c\in\bbR^n$ von $A$ und $b$ abhängen. So erhält man eine Vielzahl verschiedener Lösungsverfahren, die je nach Eigenschaften von $A$ vorteilhafte Eigenschaften wie beispielsweise besonders schnelle Konvergenz haben können. Die Konvergenz der Iterierten zur exakten Lösung ist dabei keinesfalls garantiert, sondern kann nur unter recht engen Voraussetzungen an $M$ und $c$ abgesichert werden.

3.8.5.2.2CG-Verfahren¶

Für das Verfahren der konjugierten Gradienten (kurz: CG-Verfahren, CG = conjugate gradients) wird das Gleichungssystem durch das Minimierungsproblem

ersetzt. Man kann zeigen, dass die Lösung des Minimierungsproblems mit der Lösung des Gleichungssystems übereinstimmt.

Zur Lösung des Minimierungsproblems stehen eine Vielzahl numerischer Optimierungsverfahren zur Verfügung. Das einfachste ist das so genannte Gradienten-Verfahren, welches, beginnend an einem Punkt $x_0$, den Gradient $\nabla J$ der Zielfunktion $J$ berechnet und dann die aktuelle Position in Richtung des negativen Gradienten (Richtung des steilsten Abstiegs!) verändert:

Dabei ist $s_k$ die Schrittweite, welche auf verschiedene Arten gewählt werden kann, z.B. konstant. Für den Gradient erhält man

Wählt man $s_k$ stets so, dass $J(x_{k-1}-s\,\nabla J(x_{k-1}))$ für $s=s_k$ minimal wird, so erhält man das CG-Verfahren, welches besonders schnell konvergiert. Für konstantes $s_k$ erhält man hingegen das Verfahren des steilsten Abstiegs, welches verhältnismäßig langsam konvergiert.

3.8.5.2.3Verfahren für dünn besetzte Matrizen¶

In zahlreichen Anwendungen treten dünn besetzte (englisch: sparse) sehr große Matrizen $A\in\bbR^{n\times n}$ auf, also Matrizen, die nur wenige von Null verschiedene Einträge haben. Die Anzahl der Nicht-Null-Einträge ist üblicherweise ein kleines Vielfaches von $n$. Diese Matrizen speichert man nicht als Block aus $n^2$ Zahlen, sondern als Liste der Nicht-Null-Einträge und der zugehörigen Zeilen- und Spalten-Indizes.

Beispielsweise sind Adjazenzmatrizen von Graphen oft dünnbesetzt (und groß). Auch beim numerischen Lösen vieler naturwissenschaftlicher Probleme treten große dünnbesetzte Matrizen auf.

Der Zugriff auf einen durch Zeilen- und Spaltenindex gegebenen Matrixeintrag ist bei dünn besetzten Matrizen mit hohem Aufwand verbunden (Liste durchsuchen!). Matrix-Vektor-Produkte können jedoch schnell berechnet werden (IDV 385). Somit sind iterative Lösungsverfahren gegenüber Gauß-Algorithmus und Cholesky-Verfahren hier klar im Vorteil, da iterative Verfahren meist nur Matrix-Vektor-Produkte auswerten und keine anderweitigen Zugriffe auf die Matrixeinträge tätigen.

Heute übliche Verfahren für große dünn besetzte Gleichungssysteme sind Krylow-Unterraum-Verfahren (Verallgemeinerung des CG-Verfahrens) und Mehrgitterverfahren.

3.8.5.3Vorkonditionierung¶

Die Konvergenz iterativer Lösungsverfahren kann deutlich beschleunigt werden, wenn das Gleichungssystem im Vorfeld geeignet äquivalent umgeformt wird. Auch für nicht iterative Verfahren kann eine geeignete Umformung zur Verbesserung der Matrixkondition führen. Eine übliche Umformung ist die Multiplikation mit einer geeignet gewählten Matrix $M\in\bbR^{n\times n}$:

Diese Matrix soll leicht zu beschaffen sein und dabei die Inverse $A^{-1}$ möglichst gut annähern. Wie $M$ konkret zu wählen ist, hängt auch vom eingesetzten Lösungsverfahren ab.

Eine einfache Technik zur Vorkonditionierung ist die Äquilibrierung: Jede Zeile des Gleichungssystems wird durch die Norm der entsprechenden Zeile der Systemmatrix als Vektor geteilt. Bei Verwendung der euklidischen Norm also