Numerische lineare Algebra

Die numerische lineare Algebra ist ein zentrales Teilgebiet der numerischen Mathematik. Sie beschäftigt sich mit der Entwicklung und der Analyse von Rechenverfahren (Algorithmen) für Problemstellungen der linearen Algebra, insbesondere der Lösung von linearen Gleichungssystemen und Eigenwertproblemen. Solche Probleme spielen in allen Natur- und Ingenieurwissenschaften, aber auch in der Ökonometrie und in der Statistik eine große Rolle.

Die Algorithmen der numerischen linearen Algebra lassen sich grob in zwei Gruppen einteilen: in die direkten Verfahren, die theoretisch nach endlich vielen Rechenschritten die exakte Lösung eines Problems liefern, und in die iterativen Verfahren, bei denen die exakte Lösung schrittweise immer genauer angenähert wird. Da aber auch die direkten Verfahren wegen der beim Rechnen mit endlicher Genauigkeit entstehenden Rundungsfehler nur Näherungen für die exakte Lösung liefern, ist diese Unterscheidung nur für die Entwicklung und Untersuchung der Verfahren selbst von Bedeutung; für den praktischen Einsatz spielt sie keine große Rolle. Historisch gehen die ersten systematischen Verfahren aus beiden Gruppen – das direkte gaußsche Eliminationsverfahren und das iterative Gauß-Seidel-Verfahren – auf Carl Friedrich Gauß zurück. Beispiele für bedeutende Verfahren des 20. Jahrhunderts, die zahlreiche Verbesserungen und Weiterentwicklungen zur Folge hatten, sind das Zerlegungsverfahren von André-Louis Cholesky, das QR-Verfahren für Eigenwertprobleme von John G. F. Francis und Wera Nikolajewna Kublanowskaja sowie das CG-Verfahren von Eduard Stiefel und Magnus Hestenes als erster Vertreter der wichtigen Krylow-Unterraum-Verfahren.

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

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

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

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

eines linearen Gleichungssystems: Gesucht ist ein Vektor ${\displaystyle \mathbf {x} }$, der bei der Matrix-Vektor-Multiplikation mit der gegebenen Matrix ${\displaystyle A}$ den gegebenen Vektor ${\displaystyle \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 ${\displaystyle A}$ regulär ist, also eine Inverse ${\displaystyle A^{-1}}$ besitzt. Dann gibt es für jede rechte Seite ${\displaystyle \mathbf {b} }$ eine eindeutig bestimmte Lösung ${\displaystyle \mathbf {x} }$ des linearen Gleichungssystems, die formal als ${\displaystyle \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 ${\displaystyle A\cdot \mathbf {x} =\mathbf {b} }$ hat in diesem Fall die Matrix ${\displaystyle A}$ mehr Zeilen als Spalten. Solche überbestimmten Systeme haben im Allgemeinen keine Lösung. Man behilft sich deshalb damit, den Vektor ${\displaystyle \mathbf {x} }$ so zu wählen, dass die Differenz ${\displaystyle \mathbf {r} =A\cdot \mathbf {x} -\mathbf {b} }$, das Residuum, 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 ${\displaystyle \mathbf {x} }$ so gewählt, dass die Quadratsumme ${\displaystyle r_{1}^{2}+\ldots +r_{m}^{2}}$ minimal wird, wobei ${\displaystyle r_{1},\dotsc ,r_{m}}$ die Komponenten des Differenzvektors ${\displaystyle \mathbf {r} }$ bezeichnen. Mithilfe der euklidischen Norm lässt sich das auch so schreiben: Man wähle ${\displaystyle \mathbf {x} }$ so, dass ${\displaystyle \|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 ${\displaystyle A}$ mit ${\displaystyle n}$ Zeilen und ${\displaystyle n}$ Spalten; gesucht sind Zahlen ${\displaystyle \lambda }$ und Vektoren ${\displaystyle \mathbf {x} \neq \mathbf {0} }$, sodass die Gleichung

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

erfüllt ist. Man nennt dann ${\displaystyle \mathbf {x} }$ einen Eigenvektor von ${\displaystyle A}$ zum Eigenwert ${\displaystyle \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 ${\displaystyle S}$ und eine Diagonalmatrix ${\displaystyle D}$ mit ${\displaystyle S^{-1}\cdot A\cdot S=D}$. Die Diagonaleinträge von ${\displaystyle D}$ sind dann die Eigenwerte von ${\displaystyle A}$ und die Spalten von ${\displaystyle S}$ die zugehörigen Eigenvektoren.

Diese Probleme treten in allen Natur- und Ingenieurwissenschaften auf. Sie spielen aber auch in den Wirtschaftswissenschaften sowie in der Statistik – und damit in allen Gebieten, die sich statistischer Methoden bedienen – eine große Rolle. Lineare Gleichungssysteme beschreiben beispielsweise Modelle in der Statik, elektrische Netzwerke oder volkswirtschaftliche Verflechtungen. So scheinbar unterschiedliche Aufgabenstellungen wie die Stabilitätsuntersuchung dynamischer Systeme, Resonanzphänomene bei Schwingungen, die Bestimmung eines PageRanks oder die Hauptkomponentenanalyse in der Statistik führen alle auf Eigenwertprobleme. Lineare Gleichungen entstehen auch durch Linearisierung und Diskretisierung innerhalb anderer numerischer Verfahren. So lassen sich beispielsweise zahlreiche Modelle in Naturwissenschaft und Technik durch partielle Differentialgleichungen beschreiben. Ihre numerische Lösung durch Differenzen- oder Finite-Elemente-Verfahren führt auf Systeme mit sehr vielen Unbekannten.

In diesem Übersichtsartikel wird der Einfachheit halber angenommen, dass alle gegebenen Matrizen und Vektoren reell sind, das heißt, dass alle ihre Einträge reelle Zahlen sind. Meist lassen sich die angesprochenen Verfahren direkt auf komplexe Zahlen verallgemeinern; an die Stelle der orthogonalen Matrizen tritt dann beispielsweise ihr komplexes Pendant, die unitären Matrizen. Mitunter ist es auch vorteilhaft, ein gegebenes komplexes Problem – etwa durch Betrachtung von Real- und Imaginärteil – auf ein reelles zurückzuführen. Zusatzüberlegungen treten allerdings bei Eigenwertproblemen mit nichtsymmetrischen reellen Matrizen auf, denn diese können auch nichtreelle Eigenwerte und Eigenvektoren haben.

Bereits seit der Antike sind Lösungen konkreter Problemstellungen überliefert, die aus heutiger Sicht als lineare Gleichungssysteme angesehen werden können. Die Neun Kapitel der Rechenkunst, in denen der Stand der chinesischen Mathematik des 1. Jahrhunderts n. Chr. zusammengefasst ist, enthielten dabei bereits tabellarische Rechenvorschriften, die Eliminationsverfahren in Matrixdarstellung entsprachen.^[1] Die systematische Untersuchung linearer Gleichungssysteme setzte gegen Ende des 17. Jahrhunderts mit ihrer Formulierung mithilfe allgemeiner Koeffizienten ein. Nach Vorarbeiten von Gottfried Wilhelm Leibniz und Colin Maclaurin veröffentlichte Gabriel Cramer 1750 eine explizite Lösungsformel für beliebig viele Unbekannte mithilfe von Determinanten. Mit dieser cramerschen Regel war das Problem theoretisch vollständig gelöst, auch in Hinblick auf Existenz und Eindeutigkeit von Lösungen. Für deren praktische Berechnung erwies sich die Formel jedoch als völlig ungeeignet, weil der Rechenaufwand dabei mit der Anzahl der Unbekannten astronomisch schnell anwächst (siehe auch Cramersche Regel#Rechenaufwand).^[2]

Die erste Verwendung und Beschreibung systematischer Rechenverfahren für lineare Gleichungen geht auf Carl Friedrich Gauß (1777–1855) zurück. Zu Beginn des 19. Jahrhunderts waren die Bestimmung der Bahndaten astronomischer Objekte und die Landesvermessung durch Triangulation die wichtigsten Anwendungsaufgaben der mathematischen Praxis. 1801 erregte Gauß großes Aufsehen, als es ihm gelang, die Bahn des neu entdeckten Kleinplaneten Ceres aus wenigen Beobachtungen so genau zu bestimmen, dass Ceres Ende des Jahres wiedergefunden werden konnte. Für das zugehörige überbestimmte Gleichungssystem verwendete er die von ihm entdeckte Methode der kleinsten Quadrate. Das von ihm zur Berechnung der Lösung eingesetzte Eliminationsverfahren beschrieb Gauß systematisch ab 1809 im Rahmen der Bahnbestimmung des Asteroiden Pallas, allerdings noch direkt angewendet auf Quadratsummen.^[3]

Auch das erste Iterationsverfahren zur Lösung linearer Gleichungssysteme – das Gauß-Seidel-Verfahren – stammt von Gauß. In einem Brief an Christian Ludwig Gerling berichtete er 1823 von einem neuen einfachen Verfahren, mit dem die Lösung Schritt für Schritt immer besser angenähert werden könne. Gauß, der inzwischen mit der Triangulation des Königreichs Hannover beschäftigt war, schreibt darin, er rechne fast jeden Abend noch einen Iterationsschritt; das sei eine angenehme Abwechslung zur einförmigen Aufnahme der Messdaten. Das Verfahren sei so wenig anfällig für Fehler, dass es sich sogar „halb im Schlaf“ ausführen lasse.^[4] 1845 veröffentlichte Carl Gustav Jacob Jacobi ein anderes, ähnliches Iterationsverfahren, das Jacobi-Verfahren. Als Philipp Ludwig von Seidel, ein Schüler Jacobis, 1874 ein System mit 72 Unbekannten lösen musste, entwickelte er eine modifizierte, verbesserte Version dieser Methode. Wie sich im Nachhinein herausstellte, ist dieses Verfahren äquivalent zum Iterationsverfahren von Gauß, von dem Seidel jedoch vermutlich nichts wusste.^[5] Jacobi veröffentlichte 1846 auch ein iteratives Verfahren zur Transformation von Matrizen, das sich zur Lösung des Eigenwertproblems für symmetrische Matrizen eignet und heute ebenfalls als Jacobi-Verfahren bezeichnet wird. Er selbst verwendete es jedoch nur als Vorbereitungsschritt, um die Diagonaleinträge der Matrix stärker dominant zu machen.^[6]

Im Jahr 1923 wurde ein von André-Louis Cholesky entwickeltes Verfahren veröffentlicht, das bestimmte lineare Gleichungssysteme löst, indem die Koeffizientenmatrix in ein Produkt zweier einfacherer Matrizen zerlegt wird, die Cholesky-Zerlegung. Auch das gaußsche Eliminationsverfahren stellte sich im Nachhinein als ein Spezialfall solcher Matrixzerlegungsverfahren heraus. Algorithmen aus dieser Verfahrensgruppe sind auch heute noch die Standardverfahren zur Lösung mäßig großer Systeme.^[7]

Ab Ende der 1920er Jahre kamen auch neue Ideen zur iterativen Lösung von Eigenwertproblemen auf, beginnend 1929 mit der Vorstellung der Potenzmethode durch Richard von Mises. Wie auch bei der Weiterentwicklung zur inversen Iteration durch Helmut Wielandt 1944, können mit diesen einfachen Vektoriterationsverfahren immer nur Eigenvektoren zu einem einzelnen Eigenwert berechnet werden.^[8] Eine vollständige Lösung des Eigenwertproblems für beliebige Matrizen blieb aufwändig. Der Durchbruch kam hier 1961–1962 mit der Entwicklung des QR-Verfahrens durch den britischen Informatiker John G. F. Francis und unabhängig davon durch die russische Mathematikerin Wera Nikolajewna Kublanowskaja. Während Kublanowskaja in ihrer Arbeit von Anfang an ein tiefes Verständnis der Konvergenzeigenschaften der Methode aufzeigte, arbeitete Francis vor allem an Implementierungsdetails, die das Verfahren schnell und stabil machten. Das QR-Verfahren ist bis heute das Standardverfahren zur Berechnung aller Eigenwerte und Eigenvektoren nicht allzu großer Matrizen.^[9]

Die Lösung linearer Gleichungssysteme mit sehr großen Matrizen, wie sie bei der Diskretisierung von partiellen Differentialgleichungen auftreten, blieb weiterhin schwierig. Diese Matrizen haben nur relativ wenige Einträge, die ungleich null sind, und es ist von entscheidender Bedeutung, dass ein numerisches Verfahren diese Eigenschaft ausnutzt. Ein neuer Ansatz dazu, der sich als Startpunkt zahlreicher Weiterentwicklungen herausstellen sollte, war das 1952 von Eduard Stiefel und Magnus Hestenes entwickelte CG-Verfahren. Dabei wird das lineare Gleichungssystem in dem Spezialfall, dass die Matrix symmetrisch und zusätzlich positiv definit ist, durch ein äquivalentes Optimierungsproblem ersetzt. Als noch fruchtbarer erwies sich ein anderer Zugang zum CG-Verfahren, der gleichzeitig von Cornelius Lanczos entdeckt wurde: Die durch das CG-Verfahren berechneten Näherungen befinden sich in einer aufsteigenden Kette von Unterräumen, den Krylow-Räumen.^[10]

Trotz der Entdeckung dieser Zusammenhänge dauerte es relativ lange, bis konkrete Verallgemeinerungen des CG-Verfahrens entwickelt wurden. Das 1974 von Roger Fletcher veröffentlichte BiCG-Verfahren ist zwar theoretisch für beliebige Matrizen anwendbar, erwies sich jedoch in der Praxis in vielen Fällen als instabil. Das 1975 erschienene MINRES-Verfahren ist ein Krylow-Unterraum-Verfahren, für das die Matrix zwar symmetrisch sein muss, aber nicht unbedingt positiv definit wie beim CG-Verfahren.^[11] In der Folgezeit wurden zahlreiche Weiterentwicklungen untersucht, insbesondere Stabilisierungsversuche für das BiCG-Verfahren. Ein Beispiel für ein weit verbreitetes Krylow-Unterraum-Verfahren für beliebige lineare Gleichungssysteme ist eine Verallgemeinerung des MINRES-Verfahrens, das 1986 von Yousef Saad und Martin H. Schultz vorgestellte GMRES-Verfahren. Weitere Verfahren verwenden Synthesen aus Ideen der BiCG-Gruppe und GMRES, so das QMR-Verfahren (Roland W. Freund und Noel M. Nachtigal, 1991) sowie das TFQMR-Verfahren (Freund, 1993).^[12] Von Anfang an wurden Krylow-Unterraum-Verfahren auch zur Berechnung von Eigenwerten verwendet, Ausgangspunkte waren hier ein Verfahren von Lanczos 1950 und das Arnoldi-Verfahren von Walter Edwin Arnoldi 1951.^[13]

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 im Allgemeinen verloren geht.^[15]

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 ${\displaystyle 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 ${\displaystyle \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 die 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,^[16] 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.^[17]

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

${\displaystyle \|\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 ${\displaystyle \|\mathbf {x} \|_{\infty }}$ oder die 1-Norm ${\displaystyle \|\mathbf {x} \|_{1}}$ geeigneter sein.

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

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

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

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

Da der relative Fehler nicht durch die Skalierung von ${\displaystyle \mathbf {x} }$ und ${\displaystyle {\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.^[18]

Auch die „Größe“ von Matrizen wird mit Normen gemessen, den Matrixnormen. Für die Wahl einer Matrixnorm ${\displaystyle \|A\|}$ ist es wesentlich, dass sie zur verwendeten Vektornorm „passt“, insbesondere soll die Ungleichung ${\displaystyle \|A\mathbf {x} \|\leq \|A\|\|\mathbf {x} \|}$ für alle ${\displaystyle \mathbf {x} }$ erfüllt sein. Definiert man ${\displaystyle \|A\|}$ für eine gegebene Vektornorm als die kleinste Zahl ${\displaystyle L}$, sodass ${\displaystyle \|A\mathbf {x} \|\leq L\|\mathbf {x} \|}$ für alle ${\displaystyle \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 ${\displaystyle \|A\|_{2}}$, für die Maximumsnorm ist es die Zeilensummennorm ${\displaystyle \|A\|_{\infty }}$ und für die 1-Norm die Spaltensummennorm ${\displaystyle \|A\|_{1}}$. Analog zu Vektoren kann mithilfe einer Matrixnorm der relative Fehler

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

bei einer Näherung einer Matrix ${\displaystyle A}$ durch eine Matrix ${\displaystyle {\tilde {A}}}$ quantifiziert werden.^[19]

Bei Problemen aus der Praxis sind gegebene Größen meist mit Fehlern behaftet, den Datenfehlern. Zum Beispiel kann bei einem linearen Gleichungssystem ${\displaystyle A\mathbf {x} =\mathbf {b} }$ die gegebene rechte Seite ${\displaystyle \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 ${\displaystyle A\mathbf {x} =\mathbf {b} }$ in Wirklichkeit ein System ${\displaystyle A{\tilde {\mathbf {x} }}={\tilde {\mathbf {b} }}}$ mit einer gestörten rechten Seite ${\displaystyle {\tilde {\mathbf {b} }}}$ und dementsprechend einer „falschen“ Lösung ${\displaystyle {\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

${\displaystyle {\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 ${\displaystyle \|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 ${\displaystyle A}$ und wird mit ${\displaystyle \kappa (A)}$ bezeichnet. In realen Problemen wird meist nicht nur, wie hier dargestellt, die rechte Seite ${\displaystyle \mathbf {b} }$ fehlerbehaftet sein, sondern auch die Matrix ${\displaystyle A}$. Dann gilt eine ähnliche, kompliziertere Abschätzung, in der aber ebenfalls ${\displaystyle \kappa (A)}$ die wesentliche Kennzahl zur Bestimmung der Kondition des Problems bei kleinen Datenfehlern ist.^[20] Die Definition der Konditionszahl lässt sich auf nicht quadratische Matrizen verallgemeinern und spielt dann auch eine wesentliche Rolle bei der Analyse linearer Ausgleichsprobleme. Wie gut ein solches Problem konditioniert ist, hängt allerdings nicht nur wie bei linearen Gleichungssystemen von der Konditionszahl der Koeffizientenmatrix ${\displaystyle A}$ ab, sondern auch von der rechten Seite ${\displaystyle \mathbf {b} }$, genauer vom Winkel zwischen den Vektoren ${\displaystyle A\mathbf {x} }$ und ${\displaystyle \mathbf {b} }$.^[21] Nach dem Satz von Bauer-Fike lässt sich auch die Kondition des Eigenwertproblems mit Konditionszahlen beschreiben. Hier ist es jedoch nicht die Zahl ${\displaystyle \kappa (A)}$, mit der sich Störungen der Eigenwerte abschätzen lassen, sondern ${\displaystyle \kappa (S)}$, die Konditionszahl der Matrix ${\displaystyle S}$, die ${\displaystyle A}$ via ${\displaystyle S^{-1}AS=D}$ diagonalisiert.^[22]

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 ${\displaystyle u}$ eine Eingabegröße eines Problems und ${\displaystyle v=f(u)}$ seine exakte Lösung, aufgefasst als Wert einer Funktion ${\displaystyle f}$ angewendet auf ${\displaystyle u}$. Auch wenn man die Eingabegröße als exakt vorgegeben betrachtet, wird die Berechnung mit einem Algorithmus ein anderes, „falsches“ Ergebnis ${\displaystyle {\tilde {v}}=\operatorname {alg} (u)}$ liefern, aufgefasst als Wert einer anderen, „falschen“ Funktion ${\displaystyle \operatorname {alg} }$ ebenfalls angewendet auf ${\displaystyle u}$. Ein Algorithmus heißt vorwärtsstabil, wenn sich ${\displaystyle {\tilde {v}}}$ nicht wesentlich stärker von ${\displaystyle v}$ unterscheidet, als es aufgrund der Fehler in der Eingangsgröße ${\displaystyle u}$ und der Kondition des Problems sowieso zu erwarten wäre.^[23] 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 nach einer Idee von James H. Wilkinson zunächst eine sogenannte Rückwärtsanalyse betrachtet: Dazu wird ein ${\displaystyle {\tilde {u}}}$ bestimmt mit ${\displaystyle \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.^[24] Ein Algorithmus heißt rückwärtsstabil, wenn sich ${\displaystyle {\tilde {u}}}$ nicht wesentlich stärker von ${\displaystyle 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.^[25]

Wie die lineare Algebra zeigt, besteht ein enger Zusammenhang zwischen Matrizen und Basen des Vektorraums ${\displaystyle \mathbb {R} ^{n}}$. Sind ${\displaystyle n}$ linear unabhängige Vektoren ${\displaystyle \mathbf {b} _{1},\dotsc ,\mathbf {b} _{n}}$ im ${\displaystyle \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 ${\displaystyle (\mathbf {e} _{1},\mathbf {e} _{2},\mathbf {e} _{3})}$ im dreidimensionalen Raum. Fasst man die Basisvektoren spaltenweise zu einer Matrix

${\displaystyle (\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.^[26] Die Tatsache, dass bei einer orthogonalen Matrix ${\displaystyle Q}$ die Spalten eine Orthonormalbasis bilden, lässt sich in Matrixschreibweise durch die Gleichung ${\displaystyle Q^{T}Q=I}$ ausdrücken, wobei ${\displaystyle Q^{T}}$ die transponierte Matrix und ${\displaystyle I}$ die Einheitsmatrix bezeichnen. Das zeigt wiederum, dass eine orthogonale Matrix regulär ist und ihre Inverse gleich ihrer Transponierten ist: ${\displaystyle Q^{-1}=Q^{T}}$. Die Lösung eines linearen Gleichungssystems ${\displaystyle Q\mathbf {x} =\mathbf {b} }$ lässt sich daher sehr einfach bestimmen, es gilt ${\displaystyle \mathbf {x} =Q^{T}\mathbf {b} }$. Eine andere grundlegende Eigenschaft ist es, dass eine Multiplikation eines Vektors mit einer orthogonalen Matrix seine euklidische Norm unverändert lässt

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

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

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

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

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 ${\displaystyle A}$, für die ${\displaystyle A^{T}=A}$ gilt, existiert eine orthogonale Matrix ${\displaystyle Q}$ und eine Diagonalmatrix ${\displaystyle D}$ mit

${\displaystyle Q^{T}AQ=D}$.

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

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

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

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

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

Householder-Transformationen und Givens-Rotationen können also dazu benutzt werden, eine gegebene Matrix ${\displaystyle A}$ auf eine obere Dreiecksmatrix zu transformieren, oder anders ausgedrückt, eine QR-Zerlegung ${\displaystyle 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.^[30]

In der linearen Algebra wird zur Untersuchung des Eigenwertproblems ${\displaystyle A\mathbf {x} =\lambda \mathbf {x} }$ einer Matrix ${\displaystyle A}$ mit ${\displaystyle n}$ Zeilen und ${\displaystyle n}$ Spalten das charakteristische Polynom ${\displaystyle \chi _{A}(\lambda )=\det(\lambda I-A)}$ verwendet, ein Polynom vom Grad ${\displaystyle n}$. Die Eigenwerte von ${\displaystyle A}$ sind genau die Nullstellen von ${\displaystyle \chi _{A}}$. Mit dem Fundamentalsatz der Algebra ergibt sich daraus direkt, dass ${\displaystyle A}$ genau ${\displaystyle 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 ${\displaystyle 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 Rundungsfehler 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.^[31] Viele numerische Verfahren zur Berechnung von Eigenwerten und Eigenvektoren beruhen daher auf einer anderen Grundidee, den Ähnlichkeitstransformationen: Zwei quadratische Matrizen ${\displaystyle A}$ und ${\displaystyle B}$ werden ähnlich genannt, wenn es eine reguläre Matrix ${\displaystyle S}$ mit

${\displaystyle B=S^{-1}AS}$

gibt. Es kann gezeigt werden, dass zueinander ähnliche Matrizen die gleichen Eigenwerte haben, bei einer Ähnlichkeitstransformation der Matrix ${\displaystyle A}$ auf die Matrix ${\displaystyle B}$ ändern sich also die gesuchten Eigenwerte nicht. Auch die zugehörigen Eigenvektoren lassen sich leicht ineinander umrechnen: Ist ${\displaystyle \mathbf {x} }$ ein Eigenvektor von ${\displaystyle B}$, dann ist ${\displaystyle S\mathbf {x} }$ ein Eigenvektor von ${\displaystyle A}$ zum gleichen Eigenwert. Das führt zu Grundideen, die in zahlreichen Algorithmen zum Einsatz kommen. Die Matrix ${\displaystyle 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 ${\displaystyle S}$ meist orthogonale Matrizen verwendet.^[32]

Das klassische Eliminationsverfahren von Gauß zur Lösung linearer Gleichungssysteme ist ein Ausgangspunkt und Vergleichsmaßstab 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 ${\displaystyle k}$-ten Diagonalelement ${\displaystyle a_{kk}}$ der Matrix ${\displaystyle A}$ ein Element ${\displaystyle a_{ik}}$ in derselben Spalte eliminiert werden, dann muss mit dem Quotienten

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

das ${\displaystyle l_{ik}}$-fache der ${\displaystyle k}$-ten Zeile von der ${\displaystyle i}$-Zeile subtrahiert werden. Dazu muss zumindest ${\displaystyle a_{kk}\neq 0}$ gelten, was sich durch geeignete Zeilenvertauschungen für eine reguläre Matrix ${\displaystyle A}$ stets erreichen lässt. Aber mehr noch: Ist ${\displaystyle |a_{kk}|}$ sehr klein im Vergleich zu ${\displaystyle |a_{ik}|}$, dann ergäbe sich ein sehr großer Betrag von ${\displaystyle 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 Beträge ${\displaystyle |l_{ik}|}$ möglichst klein bleiben.^[33]

Die wichtigsten direkten Verfahren zur Lösung linearer Gleichungssysteme lassen sich als Faktorisierungsverfahren darstellen. Deren Grundidee ist es, die Koeffizientenmatrix ${\displaystyle A}$ des Systems ${\displaystyle A\mathbf {x} =\mathbf {b} }$ in ein Produkt aus zwei oder mehr Matrizen zu zerlegen, allgemein etwa ${\displaystyle A=BC}$. Das lineare Gleichungssystem lautet damit ${\displaystyle BC\mathbf {x} =\mathbf {b} }$ und wird in zwei Schritten gelöst: Zuerst wird die Lösung ${\displaystyle \mathbf {y} }$ des Systems ${\displaystyle B\mathbf {y} =\mathbf {b} }$ berechnet und anschließend die Lösung ${\displaystyle \mathbf {x} }$ des Systems ${\displaystyle C\mathbf {x} =\mathbf {y} }$. Es gilt dann ${\displaystyle A\mathbf {x} =BC\mathbf {x} =B\mathbf {y} =\mathbf {b} }$, also ist ${\displaystyle \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 ${\displaystyle B}$ und ${\displaystyle 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 ${\displaystyle A}$, aber unterschiedlichen rechten Seiten gelöst werden sollen: Die Faktorisierung von ${\displaystyle A}$, im Allgemeinen der aufwändigste Verfahrensschritt, muss dann nur einmal berechnet werden.

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

${\displaystyle PA=LR}$.

Nach dem Grundprinzip der Faktorisierungsverfahren werden zur Lösung von ${\displaystyle A\mathbf {x} =\mathbf {b} }$ also zunächst wie beschrieben die Dreiecksmatrizen ${\displaystyle L}$ und ${\displaystyle R}$ sowie gegebenenfalls die zugehörige Permutation bestimmt. In nächsten Schritt wird ${\displaystyle L\mathbf {y} =P\mathbf {b} }$ mit der zeilenpermutierten rechten Seite durch Vorwärtseinsetzen und schließlich ${\displaystyle 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.^[34]

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

${\displaystyle A=LL^{T}}$.

Ein allgemeiner Ansatz für die Matrixeinträge von ${\displaystyle 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 ${\displaystyle A}$ reduziert sich der Rechenaufwand gegenüber der LR-Zerlegung auf etwa die Hälfte.^[35]

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, ${\displaystyle \|A\mathbf {x} -\mathbf {b} \|_{2}^{2}}$ zu minimieren, äquivalent damit ist, das lineare Gleichungssystem

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

zu lösen. Die Koeffizientenmatrix ${\displaystyle A^{T}A}$ dieser Normalgleichungen ist symmetrisch und, wenn die Spalten von ${\displaystyle A}$ linear unabhängig sind, auch positiv definit. Es kann also mit dem Cholesky-Verfahren gelöst werden.^[36] 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 ${\displaystyle A}$ zu verwenden.

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

${\displaystyle A=QR}$

direkt nach dem allgemeinen Prinzip der Faktorisierungsverfahren gelöst werden; es ist nur noch ${\displaystyle R\mathbf {x} =\mathbf {y} }$ mit ${\displaystyle \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.^[37] Allerdings ist der Rechenaufwand im Allgemeinen etwa doppelt so groß, sodass unter Umständen eine Abwägung der Verfahren getroffen werden muss.^[38]

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 ${\displaystyle \|A\mathbf {x} -\mathbf {b} \|_{2}^{2}}$

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

${\displaystyle \|\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 ${\displaystyle Q^{T}}$ orthogonal ist, also die euklidische Norm erhält, und dass ${\displaystyle Q^{T}Q=I}$ gilt. Der letzte Ausdruck lässt sich einfach durch Rückwärtseinsetzen der ersten ${\displaystyle n}$ Zeilen von ${\displaystyle R\mathbf {x} =Q^{T}\mathbf {b} }$ minimieren.^[39]

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

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

mit einer geeigneten Matrix ${\displaystyle M}$ und einem geeigneten Vektor ${\displaystyle \mathbf {c} }$. Es lässt sich beweisen, dass solche Verfahren genau dann konvergieren, wenn alle Eigenwerte von ${\displaystyle M}$ einen Betrag echt kleiner als 1 haben. In diesem Fall konvergieren die Iterierten ${\displaystyle \mathbf {x} _{k}}$ gegen eine Lösung der Gleichung ${\displaystyle \mathbf {x} =M\mathbf {x} +\mathbf {c} }$, also gegen einen Fixpunkt der Iterationsfunktion ${\displaystyle 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 ${\displaystyle A}$ in eine Summe

${\displaystyle A=B+C}$

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

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

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

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

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

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

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

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

Ein einfaches, aber zuverlässiges iteratives Verfahren zur Lösung des Eigenwertproblems für symmetrische Matrizen ist das Jacobi-Verfahren.^[43] Es erzeugt durch sukzessive Ähnlichkeitstransformationen mit Givens-Rotationen eine Folge von symmetrischen Matrizen, die alle ähnlich zu der gegebenen symmetrischen Matrix ${\displaystyle A}$ sind und gegen eine Diagonalmatrix ${\displaystyle D}$ konvergieren. Bricht man das Verfahren nach einer geeigneten Anzahl von Schritten ab, erhält man deshalb mit den Diagonaleinträgen von ${\displaystyle D}$ Näherungen für die gesuchten Eigenwerte von ${\displaystyle A}$.

In jedem Schritt wird die Givens-Rotation, in diesem Zusammenhang auch als Jacobi-Rotation bezeichnet, so gewählt, dass der Eintrag an der Matrixposition ${\displaystyle (i,j)}$ und der symmetrisch dazu liegende bei ${\displaystyle (j,i)}$ zu null transformiert werden. Dabei ist jedoch zu beachten, dass bei dieser Transformation die ganze ${\displaystyle i}$-te und ${\displaystyle j}$-te Zeile sowie die ganze ${\displaystyle i}$-te und ${\displaystyle j}$-te Spalte der Matrix geändert wird. Deshalb werden die in einem Schritt erzeugten Nullen im Allgemeinen in den folgenden Schritten wieder zunichtegemacht. Dennoch konvergieren bei geeigneter Wahl der Positionen für die Jacobi-Rotationen alle Nichtdiagonalelemente gegen null. Das klassische Jacobi-Verfahren wählt dazu in jedem Iterationsschritt diejenige Position ${\displaystyle (i,j)}$, an der sich das Nichtdiagonalelement mit dem größten Absolutbetrag befindet. Bei einer Handrechnung war diese Position normalerweise schnell zu erkennen, bei der Umsetzung als Computerprogramm ist der Aufwand für die Suche danach im Vergleich zu den übrigen Rechenoperationen jedoch erheblich. Daher wird heute meist das zyklische Jacobi-Verfahren verwendet. Dabei werden die Positionen in einer vorher fest gewählten Reihenfolge zyklisch durchlaufen, etwa einfach spaltenweise.^[44] Es lässt sich beweisen, dass sowohl das klassische als auch das zyklische Jacobi-Verfahren stets konvergieren. Im Vergleich zu moderneren Algorithmen ist die Konvergenz allerdings relativ langsam. Für dünnbesetzte Matrizen ist das Jacobi-Verfahren nicht geeignet, da im Laufe der Iteration die Matrix mit immer mehr Nichtnulleinträgen aufgefüllt wird.^[45]

Eine einfache Ausgangsidee zur Berechnung von Eigenvektoren einer Matrix ${\displaystyle A}$ ist die Potenzmethode. Ein Startvektor ${\displaystyle \mathbf {x} _{0}\neq \mathbf {0} }$ wird iterativ immer wieder mit ${\displaystyle A}$ multipliziert

${\displaystyle \mathbf {x} _{k+1}=A\mathbf {x} _{k}}$

oder, ausgedrückt mit der ${\displaystyle k}$-ten Matrixpotenz, es wird ${\displaystyle \mathbf {x} _{k}=A^{k}\mathbf {x} _{0}}$ berechnet. Dahinter steckt die geometrische Anschauung, dass der Vektor durch ${\displaystyle A}$ in jedem Schritt am stärksten in die Richtung des Eigenvektors mit dem größten Eigenwert gestreckt wird. In dieser einfachen Form ist die Vektoriteration jedoch für die Praxis ungeeignet, da im Allgemeinen die Einträge von ${\displaystyle \mathbf {x} _{k}}$ schnell sehr klein oder sehr groß werden. Daher wird der Vektor in jedem Schritt zusätzlich zur Multiplikation mit ${\displaystyle A}$ noch mit einer Vektornorm auf ${\displaystyle \|\mathbf {x} _{k}\|=1}$ normiert. Man kann dann unter gewissen Voraussetzungen an die Lage der Eigenwerte beweisen, dass dieses Verfahren bis auf möglicherweise einen skalaren Vorfaktor tatsächlich gegen einen Eigenvektor zum betragsgrößten Eigenwert konvergiert.

Wendet man diese Idee formal auf die inverse Matrix ${\displaystyle A^{-1}}$ an, so erhält man einen Eigenvektor zum betragskleinsten Eigenwert von ${\displaystyle A}$. Hierzu wird freilich nicht die Inverse selbst berechnet, sondern es wird in jedem Schritt das lineare Gleichungssystem

${\displaystyle A\mathbf {x} _{k+1}=\mathbf {x} _{k}}$

gelöst. Eine weitere Verallgemeinerung der Idee erhält man mithilfe eines sogenannten Shiftparameters ${\displaystyle \sigma }$. Ein Eigenvektor von ${\displaystyle A}$ zu dem am nächsten bei ${\displaystyle \sigma }$ liegenden Eigenwert ist nämlich auch ein Eigenvektor zum betragskleinsten Eigenwert der „geshifteten“ Matrix ${\displaystyle A-\sigma I}$. Mit der zugehörigen Iteration

${\displaystyle (A-\sigma I)\mathbf {x} _{k+1}=\mathbf {x} _{k}}$

und Normierung von ${\displaystyle \mathbf {x} _{k}}$ in jedem Schritt ergibt sich das Verfahren der inversen Vektoriteration.

Vektoriterationsverfahren berechnen also zunächst einen bestimmten Eigenvektor von ${\displaystyle A}$, der zugehörige Eigenwert kann mithilfe des Rayleigh-Quotienten erhalten werden. Sie sind offenbar dann gut geeignet, wenn – wie häufig in bestimmten Anwendungsfällen – nur der größte, nur der kleinste oder allgemeiner nur ein einzelner Eigenwert mitsamt seinem Eigenvektor gesucht ist.^[46]

Das QR-Verfahren ist zurzeit der wichtigste Algorithmus zur Berechnung aller Eigenwerte und Eigenvektoren von nicht zu großen vollbesetzten Matrizen ${\displaystyle A}$.^[47] Es ist ein Iterationsverfahren, das in jedem Schritt eine QR-Zerlegung verwendet, um durch wiederholte Ähnlichkeitstransformationen eine Matrixfolge zu erzeugen, die schnell gegen eine obere Dreiecksmatrix konvergiert. Startend mit der Ausgangsmatrix ${\displaystyle A_{0}:=A}$ wird in seiner Grundidee im ${\displaystyle k}$-ten Schritt die Matrix QR-zerlegt,

${\displaystyle A_{k}=Q_{k}R_{k}}$,

und anschließend werden die beiden Faktoren in umgekehrter Reihenfolge wieder zusammenmultipliziert:

${\displaystyle A_{k+1}:=R_{k}Q_{k}}$,

um die neue Näherungsmatrix zu erhalten. Wegen ${\displaystyle Q_{k}^{-1}=Q_{k}^{T}}$ ergibt sich ${\displaystyle R_{k}=Q_{k}^{T}A_{k}}$ und daraus ${\displaystyle A_{k+1}=Q_{k}^{T}A_{k}Q_{k}}$; es handelt sich bei dieser Umformung also tatsächlich um eine Ähnlichkeitstransformation mit einer orthogonalen Matrix. Wie eine genauere Analyse zeigt, besteht ein enger Zusammenhang zur Potenzmethode: Die QR-Iteration lässt sich auffassen als eine Potenzmethode, die simultan auf alle Vektoren einer Orthonormalbasis angewendet wird; durch die QR-Zerlegung in jedem Schritt wird dabei sichergestellt, dass diese Vektoren im Laufe der Iteration auch numerisch stabil orthonormiert bleiben (siehe auch Unterraumiteration). Aus dieser Darstellung ergibt sich auch ein Beweis, dass das Verfahren unter geringen Voraussetzungen an ${\displaystyle A}$ gegen eine obere Dreiecksmatrix konvergiert.^[48]

In dieser einfachen Form ist das QR-Verfahren aus zwei Gründen noch nicht für die Praxis geeignet. Zum einen ist der Rechenaufwand für die QR-Zerlegung, die in jedem Schritt bestimmt werden muss, sehr groß. Zum anderen findet die Konvergenz im Allgemeinen nur langsam statt, es müssen also viele Schritte durchgeführt werden, um eine gewünschte Genauigkeit zu erhalten. Dem ersten Punkt lässt dich dadurch begegnen, dass in einem Vorbereitungsschritt die Matrix ${\displaystyle A}$ durch Ähnlichkeitstransformationen auf Hessenberg-Gestalt gebracht wird. Das lässt sich durch ${\displaystyle n-2}$ Transformationen mit geeigneten Householder-Matrizen erreichen. Da eine Hessenberg-Matrix nur noch ${\displaystyle n-1}$ Nichtnulleinträge unter der Hauptdiagonale hat, lässt sie sich schnell mit den entsprechenden ${\displaystyle n-1}$ Givens-Rotationen QR-zerlegen. Wie sich leicht zeigen lässt, erhält ein Schritt des QR-Verfahrens Symmetrie und Hessenberg-Gestalt. Da eine symmetrische Hessenberg-Matrix eine Tridiagonalmatrix ist, vereinfacht sich das Verfahren im symmetrischen Fall nochmals erheblich. Die Konvergenzgeschwindigkeit kann ähnlich wie bei der inversen Vektoriteration deutlich erhöht werden, wenn in jedem Schritt anstelle der Matrix ${\displaystyle A_{k}}$ die Matrix ${\displaystyle A_{k}-\sigma _{k}I}$ mit einem geschickt gewählten Shiftparameter ${\displaystyle \sigma _{k}}$ transformiert wird. Für die Wahl von ${\displaystyle \sigma _{k}}$, der Wert sollte eine Näherung an einen Eigenwert von ${\displaystyle A}$ sein, existieren verschiedene sogenannte Shiftstrategien.^[49]

Mit einer Variante des QR-Verfahrens kann auch die sogenannte Singulärwertzerlegung einer Matrix berechnet werden.^[50] Diese Verallgemeinerung der Diagonalisierung auf beliebige – sogar nicht quadratische – Matrizen wird in einigen Anwendungen, wie etwa in der Bildkompression, direkt verwendet. Mithilfe der Singulärwertzerlegung können auch große, schlecht konditionierte lineare Ausgleichsprobleme gelöst werden.^[51]

Die Krylow-Unterraum-Verfahren mit ihren zahlreichen Varianten und Spezialisierungen sind die wichtigste Verfahrensgruppe zur Lösung sowohl von linearen Gleichungssystemen als auch von Eigenwertproblemen, wenn die gegebene Matrix ${\displaystyle A}$ sehr groß und dünnbesetzt ist. Der historisch erste Algorithmus aus dieser Gruppe ist das Verfahren der konjugierten Gradienten, kurz CG-Verfahren (von englisch conjugate gradients) zur Lösung linearer Gleichungssysteme mit symmetrischen und positiv definiten Koeffizientenmatrizen.

Der fruchtbare Zusammenhang des CG-Verfahrens mit Krylow-Unterräumen wurde erst später erkannt, seine Grundidee ist eine andere: Es löst anstelle des Gleichungssystems ein dazu äquivalentes Optimierungsproblem. Ist nämlich ${\displaystyle A}$ symmetrisch und positiv definit, so ist die Lösung von ${\displaystyle A\mathbf {x} =\mathbf {b} }$ die eindeutig bestimmte Minimalstelle der Funktion

${\displaystyle E(\mathbf {x} )={\frac {1}{2}}\mathbf {x} ^{T}A\mathbf {x} -\mathbf {b} ^{T}\mathbf {x} }$.

Damit stehen grundsätzlich alle numerischen Verfahren zur Lösung von Optimierungsproblemen auch für das lineare Gleichungssystem zur Verfügung, insbesondere die sogenannten Abstiegsverfahren. Das sind iterative Verfahren, die ausgehend von der aktuellen Näherung ${\displaystyle \mathbf {x} _{k}}$ im ${\displaystyle k}$-ten Schritt entlang einer geeigneten Suchrichtung ${\displaystyle \mathbf {s} _{k}}$ das eindimensionale Optimierungsproblem

„Suche ${\displaystyle \alpha _{k}\in \mathbb {R} }$, sodass ${\displaystyle E(\mathbf {x} _{k}+\alpha _{k}\mathbf {s} _{k})}$ minimal wird.“

lösen. Die dabei gefundene Stelle ${\displaystyle \mathbf {x} _{k+1}=\mathbf {x} _{k}+\alpha _{k}\mathbf {s} _{k}}$ wird die neue Näherung für den nächsten Schritt. Eine zunächst naheliegende Wahl für die Suchrichtung ${\displaystyle \mathbf {s} _{k}}$ ist die Richtung des steilsten Abstiegs, was auf das Gradientenverfahren zur Bestimmung der Minimalstelle führt. Allerdings zeigt sich, dass die so berechneten Näherungen ${\displaystyle \mathbf {x} _{k}}$ sich im Allgemeinen nur sehr langsam und in einem „Zickzackkurs“ der wahren Lösung annähern.^[52] Wesentlich besser geeignet sind Suchrichtungen, die die spezielle Gestalt der zu minimierenden Funktion ${\displaystyle E}$ berücksichtigen. Die Niveaumengen von ${\displaystyle E}$ sind ${\displaystyle (n-1)}$-dimensionale Ellipsoide (im anschaulichen, zweidimensionalen Fall Ellipsen), daher ist es günstig, die Suchrichtungen zueinander konjugiert zu wählen (im Anschauungsfall entspricht das den konjugierten Durchmessern). Dabei heißen zwei Richtungen ${\displaystyle \mathbf {s} _{k}}$ und ${\displaystyle \mathbf {s} _{l}}$ konjugiert, wenn ${\displaystyle \mathbf {s} _{k}^{T}A\mathbf {s} _{l}=0}$ gilt. Das CG-Verfahren wählt daher für die erste Suchrichtung die Richtung des steilsten Abstiegs, aber die folgenden so, dass alle Suchrichtungen zueinander konjugiert sind. Es lässt sich zeigen, dass dann nach ${\displaystyle n}$ Abstiegen die wahre Lösung erreicht wird. Meist ist aber eine ausreichend genaue Näherungslösung schon nach deutlich weniger Schritten erreicht und das Verfahren kann vorzeitig abgebrochen werden.^[53]

In den Rechenschritten des CG-Verfahrens geht die Matrix ${\displaystyle A}$ nur in der Form von Matrix-Vektor-Produkten ein. Sie selbst wird nicht zerlegt oder umgeformt – ein großer Vorteil, wenn sie dünnbesetzt ist. Nimmt man zur Vereinfachung (aber ohne Beschränkung der Allgemeinheit) an, dass als Startvektor ${\displaystyle \mathbf {x} _{0}}$ der Nullvektor gewählt wird, so zeigt eine genauere Analyse, dass jede Näherung ${\displaystyle \mathbf {x} _{k}}$ eine Linearkombination der Vektoren ${\displaystyle \mathbf {b} ,A\mathbf {b} ,A^{2}\mathbf {b} ,\dotsc ,A^{k-1}\mathbf {b} }$ ist, also aus wiederholten Multiplikationen der rechten Seite ${\displaystyle \mathbf {b} }$ mit ${\displaystyle A}$ aufgebaut wird. Anders ausgedrückt: Jedes ${\displaystyle \mathbf {x} _{k}}$ liegt in einem Krylow-Unterraum

${\displaystyle {\mathcal {K}}_{k}(A,\mathbf {b} )=\operatorname {span} (\mathbf {b} ,A\mathbf {b} ,A^{2}\mathbf {b} ,\dotsc ,A^{k-1}\mathbf {b} )}$.

Diese Eigenschaft ist das Kennzeichen der Krylow-Unterraum-Verfahren: Sie erzeugen iterativ für ${\displaystyle k=0,1,2,\dotsc }$ Näherungen ${\displaystyle \mathbf {x} _{k}}$ mit ${\displaystyle \mathbf {x} _{k}\in {\mathcal {K}}_{k}(A,\mathbf {b} )}$. Dabei wird ${\displaystyle \mathbf {x} _{k}}$ zusätzlich so gewählt, dass das Residuum ${\displaystyle \mathbf {r} _{k}=A\mathbf {x} _{k}-\mathbf {b} }$ in einem noch festzulegenden Sinne möglichst klein ist. Beim CG-Verfahren ist die Bedingung nicht unbedingt naheliegend, aber für die spezielle Struktur des Problems gut geeignet: Mit der durch ${\displaystyle A^{-1}}$ gewichteten Vektornorm ${\displaystyle \|\mathbf {v} \|_{A^{-1}}={\sqrt {\mathbf {v} ^{T}A^{-1}\mathbf {v} }}}$ ist in jedem Schritt ${\displaystyle \|\mathbf {r} _{k}\|_{A^{-1}}}$ minimal.^[54] Der Nachteil liegt dabei darin, dass dies nur funktioniert, wenn ${\displaystyle A}$ tatsächlich symmetrisch und positiv definit ist, anderenfalls ist ${\displaystyle \|\cdot \|_{A^{-1}}}$ gar keine Norm. Im Allgemeinen werden die Zusatzbedingungen, die Krylow-Unterraum-Verfahren an die Wahl von ${\displaystyle \mathbf {x} _{k}\in {\mathcal {K}}_{k}}$ stellen, als sogenannte Projektionsbedingung formuliert. Man verlangt dabei, dass das Residuum ${\displaystyle \mathbf {r} _{k}}$ orthogonal zu allen Vektoren aus einem ${\displaystyle k}$-dimensionalen Unterraum ${\displaystyle {\mathcal {L}}_{k}}$ ist, in Symbolen

${\displaystyle \mathbf {r} _{k}\perp {\mathcal {L}}_{k}}$.

Die ${\displaystyle {\mathcal {L}}_{k}}$ sind normalerweise selbst Krylow-Unterräume, im einfachsten Fall, wie auch beim CG-Verfahren, zum Beispiel ${\displaystyle {\mathcal {L}}_{k}={\mathcal {K}}_{k}}$.^[55] Für die konkrete Berechnung der Näherungen werden sukzessive Orthonormalbasen der beteiligten Krylow-Unterräume aufgebaut. Das bekannte Gram-Schmidt-Verfahren zur Orthonormalisierung ist in seiner Standardform leider numerisch instabil. Es lässt sich jedoch mit einer kleinen Modifikation stabilisieren.^[56]

Aus den genannten Grundideen ergeben sich zahlreiche Variationen, Anpassungen und Verbesserungen innerhalb dieser Verfahrensklasse, von denen nur einige exemplarisch genannt werden sollen. Eine direkte Verallgemeinerung des CG-Verfahrens ist das BiCG-Verfahren. Es hebt die Einschränkung auf symmetrische Matrizen dadurch auf, dass es zusätzlich zu dem mit ${\displaystyle A}$ gebildeten Krylow-Unterräumen, auch die zur transponierten Matrix ${\displaystyle A^{T}}$ gehörigen verwendet. Eine Optimierung, die die zusätzlichen Multiplikationen mit ${\displaystyle A^{T}}$ vermeidet, ist das CGS-Verfahren. Beide Verfahrenstypen sind in vielen praktischen Fällen instabil, bilden aber die Grundlage für verschiedene Stabilisierungsversuche, etwa in der Gruppe der BiCGSTAB-Verfahren. Wichtige und im Allgemeinen stabile Verfahren sind GMRES und seine Spezialisierung für symmetrische Matrizen, MINRES. Sie setzen direkt bei den Residuen an und bestimmen ${\displaystyle \mathbf {x} _{k}}$ im Krylow-Unterraum so, dass ${\displaystyle \|A\mathbf {x} _{k}-\mathbf {b} \|_{2}}$ minimal ist. Weitere Verbesserungen dieses Grundprinzips sind etwa das QMR- und das TFQMR-Verfahren.^[57]

Krylow-Unterraum-Verfahren können nicht nur für sehr große dünnbesetzte lineare Gleichungssysteme verwendet werden, sondern auch zur Lösung ebensolcher Eigenwertprobleme – ein weiterer Grund für ihre große Bedeutung in der modernen numerischen linearen Algebra. Natürlich kann in Eigenwertproblemen nicht mit ${\displaystyle \mathbf {x} _{0}=\mathbf {0} }$ gestartet werden (${\displaystyle \mathbf {0} }$ ist ja per Definition kein Eigenvektor). Es werden hier die Näherungen ${\displaystyle \mathbf {x} _{k}}$ und zugehörige ${\displaystyle \lambda _{k}}$ so bestimmt, dass

${\displaystyle A\mathbf {x} _{k}-\lambda _{k}\mathbf {x} _{k}\perp {\mathcal {K}}_{k}(A,\mathbf {x} _{0})}$

mit ${\displaystyle \mathbf {x} _{0}\neq \mathbf {0} }$ gilt. Dieses Vorgehen führt auf ein nur ${\displaystyle k}$-dimensionales Eigenwertproblem, das sich für kleine ${\displaystyle k}$ leicht lösen lässt und Näherungen an einige Eigenwerte von ${\displaystyle A}$ liefert.^[58] Der zugehörige Grundalgorithmus ist das Arnoldi-Verfahren. Wie stets bei Eigenwertproblemen ergeben sich für symmetrische Matrizen deutliche Vereinfachungen; diese führen auf das Lanczos-Verfahren.^[59]

• Steffen Börm, Christian Mehl: Numerical Methods for Eigenvalue Problems. Walter de Gruyter, Berlin/Boston 2012, ISBN 978-3-11-025033-6. 
• Folkmar Bornemann: Numerische lineare Algebra – Eine konzise Einführung mit MATLAB und Julia. Springer Spektrum, Wiesbaden 2016, ISBN 978-3-658-12883-8. 
• 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. 

• Jean-Luc Chabert u. a. (Hrsg.): A History of Algorithms. Springer, Berlin/Heidelberg 1999, ISBN 978-3-540-63369-3. 
• Yousef Saad, Henk A. van der Vorst: Iterative solution of linear systems in the 20th century. In: Journal of Computational and Applied Mathematics. Band 123, 2000, S. 1–33. 
• Gene H. Golub, Henk A. van der Vorst: Eigenvalue computation in the 20th century. In: Journal of Computational and Applied Mathematics. Band 123, 2000, S. 35–65. 

• Numerische Lineare Algebra. Institut Computational Mathematics der TU Braunschweig, abgerufen am 20. August 2018 (Überblick über das Fachgebiet). 
• Peter Spellucci: Numerische Lineare Algebra. (PDF) 16. Juli 2009, abgerufen am 20. August 2018 (Vorlesungsskript der TU Darmstadt). 
• Frank Wübbeling: Numerische Lineare Algebra im WS 2012/13. (PDF) 29. März 2013, abgerufen am 20. August 2018 (Vorlesungsskript der Universität Münster). 
• Netlib. Abgerufen am 20. August 2018 (Sammlung freier Software zur numerischen linearen Algebra).