6 Lineare Gleichungssysteme

Lineare Gleichungssysteme

entstehen in allen Bereichen der Wissenschaft, entweder direkt zur Formulierung eines praktischen Problems oder als Teilschritt innerhalb komplexerer Lösungsstrategien (vgl. auch Veranstaltung zum wissenschaftlichen Rechnen). 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 hier 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. Erst im nächsten Kapitel gehen wir auf den allgemeinen Fall ein.

6.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.

6.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 sollten:

• $\|v\|\geq 0\quad$ for all $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\|_a$ und $\|\cdot\|_b$ 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.

6.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 Fobenius-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 erhalten wir

(IDVID 613). Die $\infty$-Norm liefert

(IDVID 617).

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

6.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 (IDVID 620).

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 unterschiedlich Zahlen $\kappa_1(A)$, $\kappa_2(A)$, $\kappa_\infty(A)$. Wenn im Folgenden kein Index angegeben ist, gelten die Resulte für alle Matrixnormen.

Es gilt stets

(IDVID 625).

Für die euklidische Matrixnorm kann man noch

zeigen (tun wir aber nicht), wobei $\lambda_{\mathrm{max}}$ und $\lambda_{\mathrm{min}}$ der größte und der kleinste Eigenwert von $A$ sind.

6.1.4Fehlerhafte Matrix und rechte Seite¶

Untersuchen noch die Kondition der 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 erhält man dann

(IDVID 635), 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.

6.1.5Beispiele¶

6.1.5.1Einheitsmatrix¶

Für

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

6.1.5.2Hilbert-Matrix¶

Für

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

6.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 Determinanten. 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 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 (vgl. Veranstaltung zum wissenschaftlichen Rechnen), 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).

6.3Gauß-Algorithmus¶

6.3.1Idee¶

Die LR-Zerlegung von $A$ kann mittels Gauß-Algorithmus (siehe unten) bestimmt werden. Liegt diese vor, so erfolgt das Lösen des Gleichungssystems durch Vorwärtssubstitution und anschließende Rückwärtssubstitution:

1. Löse $L\,y=b$; dazu setze $y_1:=b_1$ und wiederhole für $i=2,\ldots,n$:

   1. $\displaystyle y_i:=b_i-\sum_{j=1}^{i-1} l_{i,j}\,y_j$.

2. Löse $R\,x=y$; dazu setze $x_n:=\frac{y_n}{r_{n,n}}$ und wiederhole für $i=n-1,n-2,\ldots,1$:

   1. $\displaystyle x_i:=\frac{1}{r_{i,i}}\,\left(y_i-\sum_{j=i+1}^n r_{i,j}\,x_j\right)$.

6.3.2LR-Zerlegung¶

Sei $E_{i,j}\in\bbR^{n\times n}$ die Matrix, die in Zeile $i$ und Spalte $j$ eine 1 hat und sonst nur Nullen, und sei $I\in\bbR^{n\times n}$ die Einheitsmatrix. Für $c\in\bbR$ setzen wir

Dann entsteht das Produkt $M_{i,j}(c)\,A$ aus $A$ indem Zeile $j$ von $A$ mit $c$ multipliziert und zu Zeile $i$ von $A$ addiert wird. Umgekehrt wird bei $A\,M_{i,j}(c)$ gerade Spalte $i$ mit $c$ multipliziert und zu Spalte $j$ addiert. Offensichtlich gilt $M_{i,j}(c)^{-1}=M_{i,j}(-c)$.

Mit diesen Matrizen beschreiben wir nun den Algorithmus zum Erhalt der LR-Zerlegung. Dieser erzeugt von der Matrix $A$ ausgehend spaltenweise Nullen unterhalb der Diagonale, indem geeignete Vielfache einer Zeile von den darunter liegenden Zeilen abgezogen werden:

1. Setze $L:=I$ und $R:=A$.

2. Wiederhole für $j=1,\ldots,n-1$:

   1. Wiederhole für $i=j+1,\ldots,n$:

      1. Setze $c:=\frac{r_{i,j}}{r_{j,j}}$.

      2. Setze $L:=L\,M_{i,j}(c)$.

      3. Setze $R:=M_{i,j}(-c)\,R$.

Das Produkt $L\,R$ verändert sich im Laufe der Iteration nicht, bleibt also beim initialen $L\,R=A$. Die Zahlen $c$ sind gerade so gewählt, dass $R$ zu einer oberen Dreiecksmatrix wird. Wegen $i>j$ bleiben die Diagonaleinträge und das obere Dreieck von $L$ unverändert.

6.3.3Stabilität¶

Ist $\tilde{x}$ die mittels Gauß-Algorithmus aus exakten Eingaben $A$ und $b$ numerisch ermittelte und somit (rundungs-)fehlerbehaftete Lösung, so haben wir im Sinne der Rückwärtsanalyse

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.

Ein alternative Analyse liefert, die eine Darstellung $\tilde{x}=\hat{A}^{-1}\,b$ nutzt, findet man in Fehleranalyse für die Gauß-Elimination zur Berechnung der Lösung minimaler Länge (Satz 4.5). Auch die dort erzielte Abschätzung für den Unterschied zwischen $A$ und $\hat{A}$ liefert nicht in jedem Fall Stabilität.

Folgende Punkte fassen die Sachlage zusammen:

• Theoretisch ist der Gauss-Algorithmus instabil.

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

• Praktisch tritt Instabilität nur sehr selten auf, sofern Pivotsuche eingesetzt wird (siehe unten).

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

6.3.4Zeit- und Speicherbedarf¶

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

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.

Beachte: Die LR-Zerlegung kann in einer Matrix der Größe $n\times n$ gespeichert werden, da von $L$ nur das untere Dreieck ohne Diagonalelemente (alle 1) benötigt wird. Diese Matrix kann den Speicherbereich von $A$ nutzen, da $R$ schrittweise aus $A$ entsteht. Die Nicht-Null-Einträge von $L$ sind in jedem Schritt gerade dort, wo $R$ Nullen enthält. Sofern die Matrix $A$ nicht mehr benötigt wird, also überschrieben werden kann, benötigt der Gauß-Algorithmus also keinen zusätzlichen Speicher.

6.3.5Pivotsuche¶

Der Gauß-Algorithmus in der bisher vorgestellten Form hat zwei Schwächen:

• Alle Diagonalelemente von $A$ müssen von Null verschieden sein.

• Die Division durch $r_{j,j}$ kann bei kleinem $r_{j,j}$ zu Instabilität führen.

Beide Probleme kann man relativ leicht beheben. Zunächst ein Beispiel für die Instabilität: Rechen mit dezimalen Gleitkommazahlen mit zweistelliger Mantisse und lösen

Bei exaktem Rechnen erhält man $x_1=1.\overline{001}$ und $x_2=0.\overline{998}$, gerundet auf das genutzte Zahlenformat also $x_1=1.0$ und $x_2=1.0$. Führt man den Gauß-Algorithmus mit Runden nach jeder Operations aus, so liefert dieser $x_1=0.0$ und $x_2=1.0$ (IDVID 650), also eine grob falsche Lösung. Führt man den Algorithmus jedoch mit vertauschten Zeilen aus, so stimmt die numerische Lösung mit der gerundeten exakten Lösung überein.

Bei Gleichungssystemen ist die Reihenfolge der Gleichungen egal. Man kann also die Zeilen der Systemmatrix $A$ beliebig anordnen, solange man die Reihenfolge der Einträge in der rechten Seite $b$ entsprechend anpasst. Diese Freiheit nutzen wir zur Lösung der beiden genannten Probleme. Vor jeder neuen Iteration der äußeren Schleife bei der LR-Zerlegung sortieren wir die verbleibenden Zeilen so, dass $r_{j,j}$ (das sogenannte Pivotelement) so groß wie möglich wird. Insbesondere ist das Pivotelement dann nie Null (sonst wäre das Gleichungssystem nicht lösbar). Aus Sicht der Stabilität ist dieses Vorgehen vorteilhaft, weil dann die Faktoren $\frac{r_{i,j}}{r_{j,j}}$ klein sind, also Rundungsfehler aus den vorhergehenden Schritten so wenig wie möglich verstärkt werden.

Dieses Vorgehen wird Spaltenpivotsuche genannt, da nur innerhalb einer Spalte gesucht wird. Der zusätzliche Rechenaufwand liegt in der Größenordnung $n^2$ ($n$-mal suchen in höchstens $n$ Einträgen). Zur weiteren Erhöhung der Stabilität kann auch spaltenübergreifend gesucht werden. Dann müssen ggf. Spalten von $A$ vertauscht werden, was einem Verändern der Reihenfolge der Einträge in $x$ entspricht. Der zusätzliche Rechenaufwand steigt durch den größeren Suchbereich jedoch auf ca. $n^3$, weshalb diese vollständige Pivotsuche kaum eingesetzt wird.

6.3.6Implementierung¶

Den Gauß-Algorithmus sollte man (heute) nicht selbst implementieren, da hochwertige Implementierung in allen gängigen Software-Biliotheken verfügbar sind. Meistens basieren diese auf LAPACK, einer sehr ausgereiften Bibliothek für lineare Algebra, die in für verschiedene Prozessorarchitekturen optimierten Versionen vorliegt.

6.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.

6.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 (vgl. Veranstaltung zum wissenschaftlichen Rechnen). Beispielsweise sind Matrizen der Form $A=B^\rmT\,B$ für invertierbares $B\in\bbR^{n,n}$ stets symmetrisch und positiv definit (IDVID 655). Solche Matrizen treten unter anderem beim Berechnen von Skalarprodukten in transformierten Koordinaten auf. $B$ ist dabei die Koordinatentransformation (vgl. Hauptkomponentenanalyse (“PCA”) in den Datenwissenschaften).

6.4.2Cholesky-Zerlegung¶

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$ (IDVID 660).

Aus der LR-Zerlegung von $A$ bekommt man ohne relevante Zusatzkosten stets die Cholesky-Zerlegung. Allerdings geht es auf direktem Wege etwas schneller. Betrachten dazu die Einträge von

und von

als Unbekannte und lösen das nichtlineare Gleichungssystem

Dieses kann explizit gelöst werden (IDVID 665). Man erhält den folgenden Lösungsalgorithmus:

6.4.3Zeit- und Speicherbedarf¶

Der Rechenaufwand liegt wie bei der LR-Zerlegung in der Größenordnung $n^3$, ist aber trotzdem nur etwa halb so groß oder, bei Zwischenspeichern der Produkte $l_{j,k}\,d_{k,k}$, sogar nur ein Viertel so groß (IDVID 670).

Der Speicherplatz des unteren Dreiecks von $A$ kann beim berechnen der Zerlegung mit den Einträgen von $L$ überschrieben werden, da alle Einträge von $A$ nochmal im oberen Dreieck vorhanden sind. Die Einträge von $D$ können auf die entsprechenden Diagonaleinträge von $A$ geschrieben werden, da die Diagonalelemente von $A$ im weiteren Verlauf nicht mehr benötigt werden. Somit benötigt das Berechnen der Cholesky-Zerlegung keinen zusätzlichen Speicherplatz. Lediglich für das eventuelle Zwischenspeichern der Produkte $l_{j,k}\,d_{k,k}$ wird Speicherplatz in der Größenordnung $n^2$ benötigt.

6.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:

6.4.5Stabilität¶

Das Lösen von linearen Gleichungssystemen mittels Cholesky-Zerlegung ist rückwärtsstabil. Ist $\tilde{x}$ die aus dem Cholesky-Verfahren resultierende (rundungsfehlerbehaftete) Lösung, so findet man eine Matrix $\hat{A}$ mit $\tilde{x}=\hat{A}^{-1}\,b$. Für diese Matrix kann man

zeigen, wobei $\varepsilon$ das Maschinenepsilon ist. Die Konstante $c_n$ ist meist klein und liegt nur in Ausnahmefällen in der Nähe der garantierten oberen Schranke $c_n\leq 4\,n\,(3\,n+1)$ (siehe Accuracy and Stability of Numerical Algorithms (Abschnitt 10.1.1) für Herleitung und Details).

Insbesondere gilt das Cholesky-Verfahren als stabiler als der Gauss-Algorithmus. Ein günstiges Vertauschen von Zeilen im Vorfeld zur Verbesserung der Stabilität wie beim Gauss-Algorithmus ist beim Cholesky-Verfahren nicht zweckmäßig, da dadurch die Symmetrie der Matrix verloren gehen würde.

6.5Weitere Verfahren¶

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

6.5.1QR-Zerlegung¶

Die LU-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 Länge 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.

6.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.

6.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.

6.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.

Nur Lösung des Minimierungsproblems stehen eine vielzahl numerischer Optimierungsverfahren zur Verfügung. Das einfachste ist das sogenannte 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.

6.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-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 (vgl. Veranstaltung zum wissenschaftlichen Rechnen).

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 (IDVID 680). 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.

6.5.3Vorkonditionierung¶

Die Konvergenz iterative Lösungsverfahren kann deutlich beschleunigt werden, wenn das Gleichungssystem im Vorfeld geeignet äquivalent umgeformt wird. Auch für nicht iterative Verfahren kann eine geeignet Umformung zur Verbesserung der Matrixkondition führen. Ein ü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 Länge der entsprechenden Zeile der Systemmatrix als Vektor geteilt. Bei Verwendung der euklidischen Norm also