Numerische lineare Algebra

Die Modellierung durch finite Elemente, wie hier zur Spannungsanalyse eines Hubkolbens, führt auf lineare Gleichungssysteme mit sehr vielen Gleichungen und Unbekannten: eine für jeden Knotenpunkt des Gitternetzes

Die numerische lineare Algebra ist ein zentrales Teilgebiet der numerischen Mathematik. Sie beschäftigt sich mit der Entwicklung und der Analyse von Rechenmethoden für Problemstellungen der linearen Algebra, insbesondere der Lösung von linearen Gleichungssystemen und Eigenwertproblemen. Solche Probleme treten häufig in den Ingenieurwissenschaften, der Physik oder der Ökonometrie auf.

Einführung in die Problemstellungen

Ein – auch historisch gesehen – zentraler Anfangspunkt der elementaren linearen Algebra sind lineare Gleichungssysteme. Wir betrachten $n$ Gleichungen der Gestalt

$a_{i1}x_{1}+\ldots +a_{in}x_{n}=b_{i}$

für $n$ Unbekannte $x_{1},\dotsc ,x_{n}$. Die Koeffizienten $a_{ij}$ und $b_{i}$ sind gegebene Zahlen; die gesuchten Werte für $x_{1},\dotsc ,x_{n}$ sollen so bestimmt werden, dass alle Gleichungen erfüllt werden. Die Koeffizienten lassen sich zu einer Matrix $A=(a_{ij})$ zusammenfassen; die Zahlen $b_{i}$ und die Unbekannten $x_{j}$ bilden Spaltenvektoren $\mathbf {b} =(b_{i})$ und $\mathbf {x} =(x_{j})$. Auf diese Weise ergibt sich die Matrix-Vektor-Darstellung

$A\cdot \mathbf {x} =\mathbf {b}$

eines linearen Gleichungssystems: Gesucht ist ein Vektor $\mathbf {x}$, der bei der Matrix-Vektor-Multiplikation mit der gegebenen Matrix $A$ den gegebenen Vektor $\mathbf {b}$ ergibt. Als Teilgebiet der Numerik betrachtet auch die numerische lineare Algebra nur sogenannte korrekt gestellte Probleme, also insbesondere nur solche Probleme, die eine Lösung besitzen und bei denen die Lösung eindeutig bestimmt ist. Insbesondere wird im Folgenden stets angenommen, dass die Matrix $A$ regulär ist, also eine Inverse $A^{-1}$ besitzt. Dann gibt es für jede rechte Seite $\mathbf {b}$ eine eindeutig bestimmte Lösung $\mathbf {x}$ des linearen Gleichungssystems, die formal als $\mathbf {x} =A^{-1}\mathbf {b}$ angegeben werden kann.

Viele wichtige Anwendungen führen allerdings auf lineare Gleichungssysteme mit mehr Gleichungen als Unbekannten. In der Matrix-Vektor-Darstellung $A\cdot \mathbf {x} =\mathbf {b}$ hat in diesem Fall die Matrix $A$ mehr Zeilen als Spalten. Solche überbestimmten Systeme haben im Allgemeinen keine Lösung. Man behilft sich deshalb damit, den Vektor $\mathbf {x}$ so zu wählen, dass die Differenz $\mathbf {r} =A\cdot \mathbf {x} -\mathbf {b}$ in einem noch festzulegenden Sinn „möglichst klein“ wird. Beim mit Abstand wichtigsten Fall, dem sogenannten linearen Ausgleichsproblem, wird dazu die Methode der kleinsten Quadrate verwendet: Hierbei wird $\mathbf {x}$ so gewählt, dass die Quadratsumme $r_{1}^{2}+\ldots +r_{m}^{2}$ minimal wird, wobei $r_{1},\dotsc ,r_{m}$ die Komponenten des Differenzvektors $\mathbf {r}$ bezeichnen. Mithilfe der euklidischen Norm lässt sich das auch so schreiben: Man wähle $\mathbf {x}$ so, dass $\|A\cdot \mathbf {x} -\mathbf {b} \|_{2}^{2}$ minimal wird.

Neben den linearen Gleichungen sind die Eigenwertprobleme ein weiteres zentrales Thema der linearen Algebra. Gegeben ist hierbei eine Matrix $A$ mit $n$ Zeilen und $n$ Spalten; gesucht sind Zahlen $\lambda$ und Vektoren $\mathbf {x} \neq \mathbf {0}$, sodass die Gleichung

$A\cdot \mathbf {x} =\lambda \mathbf {x}$

erfüllt ist. Man nennt dann $\mathbf {x}$ einen Eigenvektor von $A$ zum Eigenwert $\lambda$. Das Problem alle Eigenwerte und Eigenvektoren einer Matrix zu bestimmen ist gleichbedeutend damit sie zu diagonalisieren. Das bedeutet: Man finde eine reguläre Matrix $S$ und eine Diagonalmatrix $D$ mit $S^{-1}\cdot A\cdot S=D$. Die Diagonaleinträge von $D$ sind dann die Eigenwerte von $A$ und die Spalten von $S$ die zugehörigen Eigenvektoren.

Grundprinzipien

“The field of numerical linear algebra is more beautiful, and more fundamental, than its rather dull name may suggest. More beautiful, because it is full of powerful ideas that are quite unlike those normally emphasized in a linear algebra course in a mathematics department. […] More fundamental, because, thanks to a trick of history, ‘numerical’ linear algebra is really applied linear algebra.”

„Das Fachgebiet der numerischen linearen Algebra ist schöner und grundlegender als es sein ziemlich langweiliger Name vermuten lässt. Schöner, weil es voll mächtiger Ideen ist, die ganz anders sind als diejenigen, die normalerweise in einer Vorlesung über lineare Algebra an einem mathematischen Institut als bedeutend herausgestellt werden. […] Grundlegender, weil ‚numerische‘ lineare Algebra dank eines Tricks der Geschichte in Wirklichkeit angewandte lineare Algebra ist.“

– Lloyd N. Trefethen, David Bau

Ausnutzung von Strukturen

Besetzungsstruktur einer dünnbesetzten Matrix, wie sie bei der Finite-Elemente-Methode auftritt; die kleinen schwarzen Quadrate stehen für die Matrixeinträge ungleich null

Modelle und Fragestellungen in Wissenschaft und Technik können auf Probleme der linearen Algebra mit Millionen von Gleichungen führen. Die Einträge einer Matrix mit einer Million Zeilen und Spalten benötigen im double-precision-Format 8 Terabyte Speicherplatz. Das zeigt, dass bereits die Bereitstellung der Daten eines Problems, geschweige denn seine Lösung, eine Herausforderung darstellen, wenn nicht auch seine spezielle Struktur berücksichtigt wird. Glücklicherweise führen viele wichtige Anwendungen – wie beispielsweise die Diskretisierung partieller Differentialgleichungen mit der Finite-Elemente-Methode – zwar auf sehr viele Gleichungen, in jeder einzelnen Gleichung kommen jedoch nur relativ wenige Unbekannte vor. Für die zugehörige Matrix bedeutet das, dass es in jeder Zeile nur wenige Einträge ungleich null gibt, die Matrix ist wie man sagt dünnbesetzt. Es gibt zahlreiche Methoden, um solche Matrizen effizient abzuspeichern und ihre Struktur auszunutzen. Verfahren, in denen Matrizen nur in Matrix-Vektor-Produkten vorkommen, sind für dünnbesetzte Probleme besonders gut geeignet, da dabei alle Multiplikationen und Additionen mit null, nicht explizit ausgeführt werden müssen. Algorithmen, bei denen die Matrix selbst umgeformt wird, sind hingegen meist nur schwierig zu implementieren, da dann die Dünnbesetztheit meist verloren geht.

Allgemein hat die Besetzungsstruktur, also die Anzahl und die Position der Matrixeinträge ungleich null, einen sehr großen Einfluss auf die theoretischen und numerischen Eigenschaften eines Problems. Das wird am Extremfall von Diagonalmatrizen, also Matrizen, die nur auf der Hauptdiagonale Einträge ungleich null haben, besonders deutlich. Ein lineares Gleichungssystem mit einer Diagonalmatrix kann einfach gelöst werden, indem die Einträge auf der rechten Seite durch die Diagonalelemente dividiert werden, also mittels $n$ Divisionen. Auch lineare Ausgleichsprobleme und Eigenwertprobleme sind für Diagonalmatrizen trivial. Die Eigenwerte einer Diagonalmatrix sind ihre Diagonalelemente und die zugehörigen Eigenvektoren die Standardbasisvektoren $\mathbf {e} _{1},\dotsc ,\mathbf {e} _{n}$.

Ein weiterer wichtiger Spezialfall sind die Dreiecksmatrizen, bei denen alle Einträge oberhalb oder unterhalb der Hauptdiagonale null sind. Gleichungssysteme mit solchen Matrizen können durch Vorwärts- bzw. Rückwärtseinsetzen einfach von oben nach unten bzw. von unten nach oben der Reihe nach aufgelöst werden. Die Eigenwerte von Dreiecksmatrizen sind wiederum trivialerweise die Einträge auf der Hauptdiagonale; zugehörige Eigenvektoren können ebenfalls durch Vorwärts- oder Rückwärtseinsetzen bestimmt werden. Ein weiterer häufiger Spezialfall dünnbesetzter Matrizen sind die Bandmatrizen: Hier sind nur die Hauptdiagonale und einige benachbarte Nebendiagonalen mit Einträgen ungleich null besetzt. Eine Abschwächung der oberen Dreiecksmatrizen sind die oberen Hessenbergmatrizen, bei den auch die Nebendiagonale unter der Hauptdiagonale besetzt ist. Eigenwertprobleme lassen sich mit relativ geringem Aufwand in äquivalente Probleme für Hessenberg- oder Tridiagonalmatrizen transformieren.

Aber nicht nur die Besetzungsstruktur, sondern auch andere Matrixeigenschaften spielen für Entwicklung und Analyse numerischer Verfahren eine wichtige Rolle. Viele Anwendungen führen auf Probleme mit symmetrischen Matrizen. Insbesondere die Eigenwertprobleme sind deutlich einfacher zu handhaben, wenn die gegebene Matrix symmetrisch ist, aber auch bei linearen Gleichungssystemen reduziert sich in diesem Fall der Lösungsaufwand im Allgemeinen um etwa die Hälfte. Weitere Beispiele für Typen von Matrizen, für die spezialisierte Algorithmen existieren, sind die Vandermonde-Matrizen, die Toeplitz-Matrizen und die zirkulanten Matrizen.

Fehleranalyse: Vektor- und Matrixnormen

Als Maße für die „Größe“ eines Vektors $\mathbf {x} =(x_{i})\in \mathbb {R} ^{n}$ werden in der Mathematik unterschiedliche Vektornormen verwendet. Am bekanntesten und verbreitetsten ist die euklidische Norm

$\|\mathbf {x} \|_{2}={\sqrt {\sum _{i=1}^{n}x_{i}^{2}}}$,

also die Wurzel aus der Summe der Quadrate aller Vektorkomponenten. Bei der bekannten geometrischen Veranschaulichung von Vektoren als Pfeile im zwei- oder dreidimensionalen Raum entspricht dies gerade der Pfeillänge. Je nach untersuchter Fragestellung können jedoch auch andere Vektornormen wie etwa die Maximumsnorm $\|\mathbf {x} \|_{\infty }$ oder die 1-Norm $\|\mathbf {x} \|_{1}$ geeigneter sein.

Sind $\mathbf {x} ,{\tilde {\mathbf {x} }}\in \mathbb {R} ^{n}$ Vektoren, wobei ${\tilde {\mathbf {x} }}$ als eine Näherung für $\mathbf {x}$ aufgefasst werden soll, so lässt sich mithilfe einer Vektornorm $\|\cdot \|$ die Genauigkeit dieser Näherung quantifizieren. Die Norm des Differenzvektors

$\|{\tilde {\mathbf {x} }}-\mathbf {x} \|$

wird als (normweiser) absoluter Fehler bezeichnet. Betrachtet man den absoluten Fehler im Verhältnis zur Norm des „exakten“ Vektors $\mathbf {x} \neq \mathbf {0}$ erhält man den (normweisen) relativen Fehler

${\frac {\|{\tilde {\mathbf {x} }}-\mathbf {x} \|}{\|\mathbf {x} \|}}$.

Da der relative Fehler nicht durch die Skalierung von $\mathbf {x}$ und ${\tilde {\mathbf {x} }}$ beeinflusst wird, ist dieser das Standardmaß für den Unterschied der beiden Vektoren und wird oft auch vereinfacht nur als „Fehler“ bezeichnet.

Auch die „Größe“ von Matrizen wird mit Normen gemessen, den Matrixnormen. Für die Wahl einer Matrixnorm $\|A\|$ ist es wesentlich, dass sie zur verwendeten Vektornorm „passt“, insbesondere soll die Ungleichung $\|A\mathbf {x} \|\leq \|A\|\|\mathbf {x} \|$ für alle $\mathbf {x}$ erfüllt sein. Definiert man $\|A\|$ für eine gegebene Vektornorm als die kleinste Zahl $L$, sodass $\|A\mathbf {x} \|\leq L\|\mathbf {x} \|$ für alle $\mathbf {x}$ gilt, dann erhält man die sogenannte natürliche Matrixnorm. Für jede Vektornorm gibt es also eine davon induzierte natürliche Matrixnorm: Für die euklidische Norm ist das die Spektralnorm $\|A\|_{2}$, für die Maximumsnorm ist es die Zeilensummennorm $\|A\|_{\infty }$ und für die 1-Norm die Spaltensummennorm $\|A\|_{1}$. Analog zu Vektoren kann mithilfe einer Matrixnorm der relative Fehler

${\frac {\|{\tilde {A}}-A\|}{\|A\|}}$

bei einer Näherung einer Matrix $A$ durch eine Matrix ${\tilde {A}}$ quantifiziert werden.

Kondition und Stabilität

Bei Problemen aus der Praxis sind gegebene Größen meist mit Fehlern behaftet, den Datenfehlern. Zum Beispiel kann bei einem linearen Gleichungssystem $A\mathbf {x} =\mathbf {b}$ die gegebene rechte Seite $\mathbf {b}$ aus einer Messung stammen und daher eine Messabweichung aufweisen. Aber auch bei theoretisch beliebig genau bekannten Größen lassen sich Rundungsfehler bei ihrer Darstellung im Computer als Gleitkommazahlen nicht vermeiden. Es muss also davon ausgegangen werden, dass anstelle des exakten Systems $A\mathbf {x} =\mathbf {b}$ in Wirklichkeit ein System $A{\tilde {\mathbf {x} }}={\tilde {\mathbf {b} }}$ mit einer gestörten rechten Seite ${\tilde {\mathbf {b} }}$ und dementsprechend einer „falschen“ Lösung ${\tilde {\mathbf {x} }}$ vorliegt. Die grundlegende Frage ist nun, wie stark sich Störungen der gegebenen Größen auf Störungen der gesuchten Größen auswirken. Wenn der relative Fehler der Lösung nicht wesentlich größer ist als die relativen Fehler der Eingangsgrößen spricht man von einem gut konditionierten, anderenfalls von einem schlecht konditionierten Problem. Für das Beispiel linearer Gleichungssysteme lässt sich hierzu die Abschätzung

${\frac {\|{\tilde {\mathbf {x} }}-\mathbf {x} \|}{\|\mathbf {x} \|}}\leq \|A\|\|A^{-1}\|\cdot {\frac {\|{\tilde {\mathbf {b} }}-\mathbf {b} \|}{\|\mathbf {b} \|}}$

beweisen. Das Problem ist also gut konditioniert, wenn $\|A\|\|A^{-1}\|$, das Produkt der Norm der Koeffizientenmatrix und der Norm ihrer Inversen, klein ist. Diese wichtige Kenngröße heißt Konditionszahl der Matrix $A$ und wird mit $\kappa (A)$ bezeichnet. In realen Problemen wird meist nicht nur, wie hier dargestellt, die rechte Seite $\mathbf {b}$ fehlerbehaftet sein, sondern auch die Matrix $A$. Dann gilt eine ähnliche, kompliziertere Abschätzung, in der aber ebenfalls $\kappa (A)$ die wesentliche Kennzahl zur Bestimmung der Kondition des Problems bei kleinen Datenfehlern ist. Die Definition der Konditionszahl lässt sich auf nicht quadratische Matrizen verallgemeinern und spielt dann auch eine wesentliche Rolle bei Analyse linearer Ausgleichsprobleme. Wie gut ein solches Problem konditioniert ist, hängt allerdings nicht nur wie bei linearen Gleichungssytemen von der Konditionszahl der Koeffizientenmatrix $A$ ab, sondern auch von der rechten Seite $\mathbf {b}$, genauer vom Winkel zwischen den Vektoren $A\mathbf {x}$ und $\mathbf {b}$. Nach dem Satz von Bauer-Fike lässt sich auch die Kondition des Eigenwertproblems mit Konditionszahlen beschreiben. Hier ist es jedoch nicht die Zahl $\kappa (A)$, mit der sich Störungen der Eigenwerte abschätzen lassen, sondern $\kappa (S)$, die Konditionszahl der Matrix $S$, die $A$ via $S^{-1}AS=D$ diagonalisiert.

Während die Kondition eine Eigenschaft des zu lösenden Problems ist, ist Stabilität eine Eigenschaft des dafür verwendeten Verfahrens. Ein numerischer Algorithmus liefert – auch bei exakt gedachten Eingangsdaten – im Allgemeinen nicht die exakte Lösung des Problems. Zum Beispiel muss ein iteratives Verfahren, das eine wahre Lösung schrittweise immer genauer annähert, nach endlich vielen Schritten mit der bis dahin erreichten Näherungslösung abbrechen. Aber auch bei direkten Verfahren, die theoretisch in endlich vielen Rechenschritten die exakte Lösung ergeben, kommt es bei der Umsetzung auf dem Computer bei jeder Rechenoperation zu Rundungsfehlern. In der numerischen Mathematik werden zwei unterschiedliche Stabilitätsbegriffe verwendet, die Vorwärtsstabilität und Rückwärtsstabilität. Sei dazu allgemein $u$ eine Eingabegröße eines Problems und $v=f(u)$ seine exakte Lösung, aufgefasst als Wert einer Funktion $f$ angewendet auf $u$. Auch wenn man die Eingabegröße als exakt vorgegeben betrachtet, wird die Berechnung mit einem Algorithmus ein anderes, „falsches“ Ergebnis ${\tilde {v}}=\operatorname {alg} (u)$ liefern, aufgefasst als Wert einer anderen, „falschen“ Funktion $\operatorname {alg}$ ebenfalls angewendet auf $u$. Ein Algorithmus heißt vorwärtsstabil, wenn sich ${\tilde {v}}$ nicht wesentlich stärker von $v$ unterscheidet, als es aufgrund der Fehler in der Eingangsgröße $u$ und der Kondition des Problems sowieso zu erwarten wäre.ref>Peter Deuflhard, Andreas Hohmann: Numerische Mathematik 1. 4. Auflage. Walter de Gruyter, Berlin 2008, ISBN 978-3-11-020354-7, S. 44. </ref> Mit einer formalen Definition dieses Begriffs erhält man zwar ein naheliegendes und relativ anschauliches Maß für die Stabilität, aber bei komplizierten Algorithmen ist es oft schwierig, ihre Vorwärtsstabilität zu untersuchen. Daher wird im Allgemeinen zunächst eine sogenannte Rückwärtsanalyse betrachtet: Dazu wird ein ${\tilde {u}}$ bestimmt mit $\operatorname {alg} (u)=f({\tilde {u}})$, das heißt: Der durch das Verfahren berechnete „falsche“ Wert wird aufgefasst als „richtiger“ Wert, der aber mit einem anderen Wert der Eingabegröße berechnet wurde. Ein Algorithmus heißt rückwärtsstabil, wenn sich ${\tilde {u}}$ nicht wesentlich stärker von $u$ unterscheidet, als es aufgrund der Fehler in dieser Eingangsgröße sowieso zu erwarten wäre. Es lässt sich beweisen, dass ein rückwärtsstabiler Algorithmus auch vorwärtsstabil ist.

Orthogonalität und orthogonale Matrizen

Wie die lineare Algebra zeigt, besteht ein enger Zusammenhang zwischen Matrizen und Basen des Vektorraums $\mathbb {R} ^{n}$. Sind $n$ linear unabhängige Vektoren $\mathbf {b} _{1},\dotsc ,\mathbf {b} _{n}$ im $\mathbb {R} ^{n}$ gegeben, so sind diese eine Basis des Raums und jeder andere Vektor kann eindeutig als Linearkombination der Basisvektoren dargestellt werden. Ein Basiswechsel entspricht dabei der Multiplikation gegebener Vektoren und Matrizen mit einer Transformationsmatrix. Einen wichtigen Spezialfall bilden die Orthonormalbasen. Hierbei sind die Basisvektoren paarweise orthogonal zueinander („stehen senkrecht aufeinander“) und sind zudem alle auf euklidische Länge 1 normiert, so wie die Standardbasis $(\mathbf {e} _{1},\mathbf {e} _{2},\mathbf {e} _{3})$ im dreidimensionalen Raum. Fasst man die Basisvektoren spaltenweise zu einer Matrix

$(\mathbf {b} _{1}|\cdots |\mathbf {b} _{n})$

zusammen, so erhält man im Fall einer Orthonormalbasis eine sogenannte orthogonale Matrix.

Orthonormalbasen und orthogonale Matrizen besitzen zahlreiche bemerkenswerte Eigenschaften, auf denen die wichtigsten Verfahren der modernen numerischen linearen Algebra basieren. Die Tatsache, dass bei einer orthogonalen Matrix $Q$ die Spalten eine Orthonormalbasis bilden, lässt sich in Matrixschreibweise durch die Gleichung $Q^{T}Q=I$ ausdrücken, wobei $Q^{T}$ die transponierte Matrix und $I$ die Einheitsmatrix bezeichnen. Das zeigt wiederum, dass eine orthogonale Matrix regulär ist und ihre Inverse gleich ihrer Transponierten ist: $Q^{-1}=Q^{T}$. Die Lösung eines linearen Gleichungssystems $Q\mathbf {x} =\mathbf {b}$ lässt sich daher sehr einfach bestimmen, es gilt $\mathbf {x} =Q^{T}\mathbf {b}$. Ein andere grundlegende Eigenschaft ist es, dass eine Multiplikation eines Vektors mit einer orthogonalen Matrix seine euklidische Norm unverändert lässt

$\|Q\mathbf {x} \|_{2}=\|\mathbf {x} \|_{2}$.

Damit folgt für die Spektralnorm $\|Q\|_{2}=1$ und für die Konditionszahl ebenfalls

$\kappa (Q)=\|Q\|_{2}\|Q^{-1}\|_{2}=1$,

denn $Q^{-1}=Q^{T}$ ist ebenfalls eine orthogonale Matrix. Multiplikationen mit orthogonalen Matrizen bewirken also keine Vergrößerung des relativen Fehlers.

Orthogonale Matrizen spielen auch eine wichtige Rolle in der Theorie und der numerischen Behandlung von Eigenwertproblemen. Nach der einfachsten Version des Spektralsatzes lassen sich symmetrische Matrizen orthogonal diagonalisieren. Damit ist gemeint: Zu einer Matrix $A$, für die $A^{T}=A$ gilt, existiert eine orthogonale Matrix $Q$ und eine Diagonalmatrix $D$ mit

$Q^{T}AQ=D$.

Auf der Diagonale von $D$ stehen die Eigenwerte von $A$ und die Spalten von $Q$ bilden eine Orthonormalbasis aus Eigenvektoren. Mit der sogenannten schurschen Normalform existiert eine Verallgemeinerung dieser orthogonalen Transformation für nicht symmetrische Matrizen. Insbesondere ist nach dem oben erwähnten Satz von Bauer-Fike das Eigenwertproblem für symmetrische Matrizen stets gut konditioniert.

Es gibt zwei spezielle, leicht handhabbare Arten orthogonaler Matrizen, die in zahllosen konkreten Verfahren der numerischen linearen Algebra zum Einsatz kommen: die Housholder-Matrizen und die Givens-Rotationen. Householder-Matrizen haben die Gestalt

$H=I-2\mathbf {v} \mathbf {v} ^{T}$

mit einem Vektor $\mathbf {v} \in \mathbb {R} ^{n}$ mit $\|\mathbf {v} \|_{2}=1$. Geometrisch beschreiben sie Spiegelungen des $n$-dimensionalen Raums $\mathbb {R} ^{n}$ an der $(n-1)$-dimensionalen Hyperebene durch den Nullpunkt, die orthogonal zu $\mathbf {v}$ ist. Ihre wesentliche Eigenschaft ist die folgende: Zu einem gegebenen Vektor $\mathbf {a}$ lässt sich leicht ein Vektor $\mathbf {v}$ bestimmen, sodass die zugehörige Householder-Matrix $H=I-2\mathbf {v} \mathbf {v} ^{T}$ den Vektor $a$ auf ein Vielfaches von $\mathbf {e} _{1}$ transformiert: $H\mathbf {a} =\sigma \mathbf {e} _{1}$ mit $\sigma =\pm \|a\|_{2}$. Dieses $H$ transformiert also alle Einträge von $\mathbf {a}$ bis auf den ersten zu null. Wendet man auf diese Weise geeignete Householder-Transformationen Spalte für Spalte nacheinander auf eine Matrix $A$ an, so können alle Einträge von $A$ unterhalb der Hauptdiagonale zu null transformiert werden.

Givens-Rotationen sind spezielle Drehungen des $\mathbb {R} ^{n}$, die eine zweidimensionale Ebene drehen und die anderen $n-2$ Dimensionen fest lassen. Die Transformation eines Vektors $\mathbf {a}$ mit einer Givens-Rotions verändert daher nur zwei Einträge von $\mathbf {a}$. Durch geeignete Wahl des Drehwinkels kann dabei einer der beiden Einträge auf null gesetzt wird. Während Householder-Transformation, angewendet auf Matrizen, ganze Teilspalten transformieren, können Givens-Rotationen dazu verwendet werden, gezielt einzelne Matrixeinträge zu ändern.

Mit Householder-Transformationen und Givens-Rotationen können also dazu verwendet werden eine gegebene Matrix $A$ auf eine obere Dreiecksmatrix zu transformieren, oder anders ausgedrückt, eine QR-Zerlegung $A=QR$ in eine orthogonale Matrix und eine obere Dreiecksmatrix zu berechnen. Die QR-Zerlegung ist ein wichtiges und vielseitiges Werkzeug, das in zahlreichen Verfahren aus allen Bereichen der numerischen linearen Algebra zum Einsatz kommt.

Ähnlichkeitstransformationen

In der linearen Algebra wird zur Untersuchung des Eigenwertproblems $A\mathbf {x} =\lambda \mathbf {x}$ einer Matrix $A$ mit $n$ Zeilen und $n$ Spalten das charakteristische Polynom $\chi _{A}(\lambda )=\det(\lambda I-A)$ verwendet, ein Polynom vom Grad $n$. Die Eigenwerte von $A$ sind genau die Nullstellen von $\chi _{A}$. Mit dem Fundamentalsatz der Algebra ergibt sich daraus direkt, dass $A$ genau $n$ Eigenwerte besitzt, wenn sie mit ihrer Vielfachheit gezählt werden. Allerdings können diese Eigenwerte, auch bei reellen Matrizen, komplexe Zahlen sein. Ist jedoch $A$ eine reelle symmetrische Matrix, dann sind ihre Eigenwerte alle reell.

Das charakteristische Polynom hat zwar eine große theoretische Bedeutung für das Eigenwertproblem, zur numerischen Berechnung ist es jedoch nicht geeignet. Das liegt vor allem daran, dass das Problem, aus gegebenen Koeffizienten die Nullstellen des zugehörigen Polynoms zu berechnen, im Allgemeinen sehr schlecht konditioniert ist: Kleine Störungen wie Rundefehler an Koeffizienten eines Polynoms können zu einer starken Verschiebung seiner Nullstellen führen. Damit würde ein gegebenenfalls gut konditioniertes Problem – die Berechnung der Eigenwerte – durch ein zwar mathematisch äquivalentes, aber schlecht konditioniertes Problem – die Berechnung der Nullstellen des charakteristischen Polynoms – ersetzt. Viele numerische Verfahren zur Berechnung von Eigenwerten und Eigenvektoren beruhen daher auf einer anderen Grundidee, den Ähnlichkeitstransformationen: Zwei quadratische Matrizen $A$ und $B$ werden ähnlich genannt, wenn es eine reguläre Matrix $S$ mit

$B=S^{-1}AS$

gibt. Es kann gezeigt werden, dass zueinander ähnliche Matrizen die gleichen Eigenwerte haben, bei einer Ähnlichkeitstransformation der Matrix $A$ auf die Matrix $B$ ändern sich also die gesuchten Eigenwerte nicht. Auch die zugehörigen Eigenvektoren lassen sich leicht ineinander umrechnen: Ist $\mathbf {x}$ ein Eigenvektor von $B$, dann ist $S\mathbf {x}$ ein Eigenvektor von $A$ zum gleichen Eigenwert. Das führt zu Grundideen, die in zahlreichen Algorithmen zum Einsatz kommen. Die Matrix $A$ wird durch Ähnlichkeitstransformation in eine Matrix überführt, für die das Eigenwertproblem effizienter zu lösen ist, oder es wird eine Folge von Ähnlichkeitstransformationen konstruiert, bei denen sich die Matrix einer Diagonal- oder Dreiecksmatrix immer weiter annähert. Aus den oben genannten Gründen werden dabei für die Transformationsmatrizen $S$ meist orthogonale Matrizen verwendet.

Verfahren und Verfahrensklassen

Gaußsches Eliminationsverfahren

Das klassische Eliminationsverfahren von Gauß zur Lösung linearer Gleichungssysteme ist ein Ausgangspunkt und Vergleichsmaßstand für weiterentwickelte Verfahren. Es wird aber auch immer noch als einfaches und zuverlässiges Verfahren – insbesondere in seiner Modifikation als LR-Zerlegung (siehe unten) – für nicht zu große, gut konditionierte Systeme in der Praxis verbreitet eingesetzt. Das Verfahren eliminiert systematisch Variablen aus den gegebenen Gleichungen, indem geeignete Vielfache einer Gleichung von einer anderen Gleichung subtrahiert werden, bis ein System in Stufenform entsteht, das der Reihe nach von unten nach oben aufgelöst werden kann.

Numerische Überlegungen kommen ins Spiel, wenn die Stabilität des Verfahrens betrachtet wird. Soll mit dem $k$-ten Diagonalelement $a_{kk}$ der Matrix $A$ ein Element $a_{ik}$ in derselben Spalte eliminiert werden, dann muss mit dem Quotienten

$l_{ik}:={\frac {a_{ik}}{a_{kk}}}$

das $l_{ik}$-fache der $k$-ten Zeile von der $i$-Zeile subtrahiert werden. Dazu muss zumindest $a_{kk}\neq 0$ gelten, was sich durch geeignete Zeilenvertauschungen für eine reguläre Matrix $A$ stets erreichen lässt. Aber mehr noch: Ist $|a_{kk}|$ sehr klein im Vergleich zu $|a_{ik}|$, dann ergäbe sich ein sehr großer Betrag von $l_{ik}$. In den nachfolgenden Schritten bestünde dann die Gefahr von Stellenauslöschungen durch Subtraktionen großer Zahlen und das Verfahren wäre instabil. Daher ist es wichtig, durch Zeilenvertauschungen, sogenannte Pivotisierung dafür zu sorgen, dass die Quotienten $l_{ik}$ möglichst klein bleiben.

Faktorisierungsverfahren

Die wichtigsten direkten Verfahren zur Lösung linearer Gleichungssysteme lassen sich als Faktorisierungsverfahren darstellen. Deren Grundidee ist es, die Koeffizientenmatrix $A$ des Systems $A\mathbf {x} =\mathbf {b}$ in ein Produkt aus zwei oder mehr Matrizen zu zerlegen, allgemein etwa $A=BC$. Das lineare Gleichungssystem lautet damit $BC\mathbf {x} =\mathbf {b}$ und wird in zwei Schritten gelöst: Zuerst wird die Lösung $\mathbf {y}$ des Systems $B\mathbf {y} =\mathbf {b}$ berechnet und anschließend die Lösung $\mathbf {x}$ des System $C\mathbf {x} =\mathbf {y}$. Es gilt dann $A\mathbf {x} =BC\mathbf {x} =B\mathbf {y} =\mathbf {b}$, also ist $\mathbf {x}$ die Lösung des ursprünglichen Problems. Auf den ersten Blick scheint dabei nur die Aufgabe, ein lineares Gleichungssystem zu lösen, durch die Aufgabe, zwei lineare Gleichungssysteme zu lösen, ersetzt zu werden. Die Idee dahinter ist es jedoch, die Faktoren $B$ und $C$ so zu wählen, dass die beiden Teilsysteme wesentlich einfacher zu lösen sind als das Ausgangssystem. Ein offensichtlicher Vorteil der Verfahrensklasse ergibt sich im Fall, dass mehrere lineare Gleichungssysteme mit derselben Koeffizientenmatrix $A$ aber unterschiedlichen rechten Seiten gelöst werden sollen: Die Faktorisierung von $A$, im Allgemeinen der aufwändigste Verfahrensschritt, muss dann nur einmal berechnet werden.

LR-Zerlegung

Das Gaußsche Eliminationsverfahren kann als Faktorisierungsverfahren aufgefasst werden. Trägt man die Koeffizienten $l_{ik}$ für $1\leq k\leq i\leq n$ in eine Matrix ein, ergibt sich ohne Zeilenvertauschungen $A=LR$ mit einer unteren Dreiecksmatrix $L$ und einer oberen Dreiecksmatrix $R$. Zusätzlich ist $L$ unipotent, das heißt alle Einträge auf der Hauptdiagonale von $L$ sind gleich 1. Wie gesehen müssen im Allgemeinen bei der Gauß-Elimination Zeilen von $A$ vertauscht werden. Das lässt sich formal mit Hilfe einer Permutationsmatrix $P$ darstellen, indem anstelle von $A$ die zeilenpermutierte Matrix $PA$ faktorisiert wird:

$PA=LR$.

Nach dem Grundprinzip der Faktorisierungsverfahren werden zur Lösung von $A\mathbf {x} =\mathbf {b}$ also zunächst wie beschrieben die Dreiecksmatrizen $L$ und $R$ sowie gegebenenfalls die zugehörige Permutation bestimmt. In nächsten Schritt wird $L\mathbf {y} =P\mathbf {b}$ mit der zeilenpermutierten rechten Seite durch Vorwärtseinsetzen und schließlich $R\mathbf {x} =\mathbf {y}$ durch Rückwärtseinsetzen gelöst.

Die LR-Zerlegung und damit das gaußsche Eliminationsverfahren ist mit geeigneter Pivotisierung „fast immer stabil“, das heißt in den meisten praktischen Anwendungsaufgaben tritt keine große Fehlerverstärkung auf. Es lassen sich jedoch pathologische Beispiele angeben, bei denen die Verfahrensfehler exponentiell mit der Anzahl der Unbekannten anwachsen.

Cholesky-Zerlegung

Die Cholesky-Zerlegung ist wie die LR-Zerlegung eine Faktorisierung der Matrix $A$ in zwei Dreiecksmatrizen für den in vielen Anwendungen auftretenden Fall, dass $A$ symmetrisch und positiv definit ist, also $A^{T}=A$ erfüllt und nur positive Eigenwerte besitzt. Unter diesen Voraussetzungen gibt es eine untere Dreiecksmatrix $L$ mit

$A=LL^{T}$.

Ein allgemeiner Ansatz für die Matrixeinträge von $L$ führt auf ein explizites Verfahren, mit dem diese spaltenweise oder zeilenweise nacheinander berechnet werden können, das Cholesky-Verfahren. Durch diese Ausnutzung der Symmetrie von $A$ reduziert sich der Rechenaufwand gegenüber der LU-Zerlegung auf etwa die Hälfte.

Symmetrische und positiv definite Koeffizientenmatrizen treten klassisch bei der Formulierung der sogenannten Normalgleichungen zur Lösung linearer Ausgleichsprobleme auf. Man kann zeigen, dass das Problem, $\|A\mathbf {x} -\mathbf {b} \|_{2}^{2}$ zu minimieren, äquivalent zur Lösung des linearen Gleichungssystems

$A^{T}A\mathbf {x} =A^{T}\mathbf {b}$

ist. Die Koeffizientenmatrix $A^{T}A$ dieser Normalgleichungen ist symmetrisch und, wenn die Spalten von $A$ linear unabhängig sind, auch positiv definit. Es kann also mit dem Cholesky-Verfahren gelöst werden. Dieses Vorgehen empfiehlt sich jedoch nur für gut konditionierte Probleme mit wenigen Unbekannten. Im Allgemeinen ist nämlich das System der Normalgleichungen deutlich schlechter konditioniert als das ursprünglich gegebene lineare Ausgleichsproblem. Es ist dann besser, nicht den Umweg über die Normalgleichungen zu gehen, sondern direkt eine QR-Zerlegung von $A$ zu verwenden.

QR-Zerlegung

Das lineare Gleichungssystem $A\mathbf {x} =\mathbf {b}$ kann nach der Berechnung einer QR-Zerlegung

$A=QR$

direkt nach dem allgemeinen Prinzip der Faktorisierungsverfahren gelöst werden; es ist nur noch $R\mathbf {x} =\mathbf {y}$ mit $\mathbf {y} =Q^{T}\mathbf {b}$ durch Rückwärtseinsetzen zu bestimmen. Aufgrund der guten Kondition orthogonaler Matrizen treten dabei die möglichen Instabilitäten der LR-Zerlegung nicht ein. Allerdings ist der Rechenaufwand im Allgemeinen etwa doppelt so groß, sodass unter Umständen eine Abwägung der Verfahren getroffen werden muss.

Die QR-Zerlegung ist auch das gängige Verfahren zur Lösung nicht zu großer, gut konditionierter linearer Ausgleichsprobleme. Für das Problem

Minimiere $\|A\mathbf {x} -\mathbf {b} \|_{2}^{2}$

gilt mit $A=QR$ und $\mathbf {r} =A\mathbf {x} -\mathbf {b}$

$\|\mathbf {r} \|_{2}^{2}=\|Q^{T}\mathbf {r} \|_{2}^{2}=\|Q^{T}QR\mathbf {x} -Q^{T}\mathbf {b} \|_{2}^{2}=\|R\mathbf {x} -Q^{T}\mathbf {b} \|_{2}^{2}$.

Dabei wurde verwendet, dass $Q^{T}$ orthogonal ist, also die euklidische Norm erhält, und dass $Q^{T}Q=I$ gilt. Der letzte Ausdruck lässt sich einfach durch Rückwärtseinsetzen der ersten $n$ Zeilen von $R\mathbf {x} =Q^{T}\mathbf {b}$ minimieren.ref>Peter Deuflhard, Andreas Hohmann: Numerische Mathematik 1. 4. Auflage. Walter de Gruyter, Berlin 2008, ISBN 978-3-11-020354-7, S. 76 f. </ref>

Fixpunktiteration mit Splitting-Verfahren

Eine völlig andere Idee, um $A\mathbf {x} =\mathbf {b}$ zu lösen, besteht darin, einen Startvektor $\mathbf {x} _{0}$ zu wählen und daraus schrittweise $\mathbf {x} _{1}=F(\mathbf {x} _{0})$, $\mathbf {x} _{2}=F(\mathbf {x} _{1})$ immer neue Näherungen an die gesuchte Lösung zu berechnen. Im Fall der Konvergenz der Folge $(\mathbf {x} _{0},\mathbf {x} _{1},\mathbf {x} _{2},\dotsc )$ gegen $\mathbf {x}$ wird dann diese Iteration nach einer geeigneten Anzahl $N$ von Schritten mit einer ausreichend genauen Näherung $\mathbf {x} _{N}$ für $\mathbf {x}$ abgebrochen. Die einfachsten und wichtigsten Verfahren dieser Art verwenden eine Iteration der Gestalt

$\mathbf {x} _{k+1}=M\mathbf {x} _{k}+\mathbf {c}$

mit einer geeigneten Matrix $M$ und einem geeigneten Vektor $\mathbf {c}$. Es lässt sich beweisen, dass solche Verfahren genau dann konvergieren, wenn alle Eigenwerte von $M$ einen Betrag echt kleiner als 1 haben. In diesem Fall konvergieren die Iterierten $\mathbf {x} _{k}$ gegen eine Lösung der Gleichung $\mathbf {x} =M\mathbf {x} +\mathbf {c}$, also gegen einen Fixpunkt der Iterationsfunktion $F(\mathbf {y} )=M\mathbf {y} +\mathbf {c}$.

Ein systematisches Vorgehen bei der Suche nach geeigneten Algorithmen dieser Gestalt ermöglicht die Idee der Splitting-Verfahren. Dabei wird die Matrix $A$ in eine Summe

$A=B+C$

zerlegt mit einer leicht zu invertierenden Matrix $B$ und dem Rest $C=A-B$. Durch Einsetzen und Umstellen ergibt sich damit aus $A\mathbf {x} =\mathbf {b}$ die Fixpunktgleichung

$\mathbf {x} =-B^{-1}C\mathbf {x} +B^{-1}\mathbf {b}$.

Mit $M=-B^{-1}C$ und $\mathbf {c} =B^{-1}\mathbf {b}$ erhält man so ein Iterationsverfahren der Gestalt $\mathbf {x} _{k+1}=M\mathbf {x} _{k}+\mathbf {c}$, das im Falle der Konvergenz die Lösung von $A\mathbf {x} =\mathbf {b}$ liefert. Die Konvergenzgeschwindigkeit ist umso größer, je kleiner der betragsgrößte Eigenwert der Iterationsmatrix $M$ ist. Dieser lässt sich auch durch beliebige Matrixnormen von $M$ abschätzen.

Als klassische Beispiele für Splitting-Verfahren verwendet das Jacobi-Verfahren für $B$ die Diagonalmatrix mit der Hauptdiagonale von $A$, das Gauß-Seidel-Verfahren den unteren Dreiecksanteil von $A$. Zur Konvergenzbeschleunigung der Fixpunktverfahren lässt sich die Idee der Relaxation nutzen. Denkt man sich die Iteration in der Form

$\mathbf {x} _{k+1}=\mathbf {x} _{k}+\Delta \mathbf {x} _{k}$

mit der Korrektur $\Delta \mathbf {x} _{k}$ im $k$-ten Schritt dargestellt, geht man mit einem geeignet gewählten Relaxationsparameter $\omega$ zu

$\mathbf {x} _{k+1}=\mathbf {x} _{k}+\omega \Delta \mathbf {x} _{k}$

über. Zum Beispiel erhält man auf diese Weise aus dem Gauß-Seidel-Verfahren das SOR-Verfahren.

Literatur

• Steffen Börm, Christian Mehl: Numerical Methods for Eigenvalue Problems. Walter de Gruyter, Berlin/Boston 2012, ISBN 978-3-11-025033-6. 
• Wolfgang Bunse, Angelika Bunse-Gerstner: Numerische lineare Algebra. B. G. Teubner, Stuttgart 1985, ISBN 978-3-519-02067-7. 
• Biswa Nath Datta: Numerical Linear Algebra and Applications. 2. Auflage. SIAM, Philadelphia 2010, ISBN 978-0-89871-685-6. 
• James W. Demmel: Applied Numerical Linear Algebra. SIAM, Philadelphia 1997, ISBN 978-0-89871-389-3. 
• Gene H. Golub, Charles F. Van Loan: Matrix Computations. 3. Auflage. The Johns Hopkins University Press, Baltimore 1996, ISBN 0-8018-5413-X. 
• Nicholas J. Higham: Accuracy and Stability of Numerical Algorithms. 2. Auflage. SIAM, Philadelphia 2002, ISBN 0-89871-521-0. 
• Andreas Meister: Numerik linearer Gleichungssysteme. 5. Auflage. Springer Spektrum, Wiesbaden 2015, ISBN 978-3-658-07199-8. 
• Granville Sewell: Computational Methods of Linear Algebra. 3. Auflage. World Scientific, Singapur 2014, ISBN 978-981-4603-85-0. 
• Lloyd N. Trefethen, David Bau, III: Numerical Linear Algebra. SIAM, Philadelphia 1997, ISBN 978-0-89871-361-9. 

Weblinks

• Numerische Lineare Algebra. Institut Computational Mathematics der TU Braunschweig, abgerufen am 2. März 2016 (Überblick über das Fachgebiet). 
• Peter Spellucci: Numerische Lineare Algebra. (PDF) 16. Juli 2009, abgerufen am 2. März 2016 (Vorlesungsskript der TU Darmstadt). 
• Frank Wübbeling: Numerische Lineare Algebra im WS 2012/13. (PDF) 29. März 2013, abgerufen am 2. März 2016 (Vorlesungsskript der Universität Münster).