Numerische lineare Algebra

Kernaufgaben und Bedeutung

Die numerische lineare Algebra ist ein zentrales Teilgebiet der numerischen Mathematik. Sie entwickelt und untersucht Algorithmen für Probleme der linearen Algebra, vor allem für lineare Gleichungssysteme und Eigenwertprobleme. Solche Aufgaben treten in Natur- und Ingenieurwissenschaften, Ökonometrie und Statistik häufig auf, zum Beispiel bei Statik, elektrischen Netzwerken, volkswirtschaftlichen Verflechtungen, Stabilitätsuntersuchungen, Resonanzphänomenen, PageRank, Hauptkomponentenanalyse sowie bei der Diskretisierung partieller Differentialgleichungen durch Differenzen- oder Finite-Elemente-Verfahren.

Ein lineares Gleichungssystem mit n Gleichungen und n Unbekannten hat die Form a_i1 x_1 + ... + a_in x_n = b_i. Die Koeffizienten werden zu einer Matrix A=(a_ij) zusammengefasst, die rechten Seiten und Unbekannten zu Vektoren b und x. In Matrix-Vektor-Schreibweise lautet das System A · x = b. Betrachtet werden hier korrekt gestellte Probleme, insbesondere solche mit einer regulären Matrix A, also einer Matrix mit Inverser A^-1. Dann existiert für jede rechte Seite b genau eine Lösung, formal x = A^-1 b.

Viele Anwendungen führen zu überbestimmten Systemen mit mehr Gleichungen als Unbekannten. Diese haben meist keine exakte Lösung. Dann sucht man x so, dass das Residuum r = A · x - b möglichst klein wird. Beim linearen Ausgleichsproblem wird die Methode der kleinsten Quadrate verwendet: Man minimiert r_1^2 + ... + r_m^2, also ||A · x - b||_2^2.

Beim Eigenwertproblem ist eine quadratische Matrix A gegeben. Gesucht sind Zahlen λ und Vektoren x ≠ 0 mit A · x = λx. Dann heißt x Eigenvektor zum Eigenwert λ. Alle Eigenwerte und Eigenvektoren zu bestimmen, entspricht der Diagonalisierung: Man sucht eine reguläre Matrix S und eine Diagonalmatrix D mit S^-1 · A · S = D. Die Diagonaleinträge von D sind die Eigenwerte, die Spalten von S die Eigenvektoren. Der Artikel betrachtet der Einfachheit halber reelle Matrizen und Vektoren; viele Verfahren lassen sich auf komplexe Zahlen übertragen.

Direkte und iterative Verfahren

Die Algorithmen der numerischen linearen Algebra werden grob in direkte und iterative Verfahren eingeteilt. Direkte Verfahren liefern theoretisch nach endlich vielen Rechenschritten die exakte Lösung. Iterative Verfahren erzeugen schrittweise immer bessere Näherungen. In der Praxis liefern aber auch direkte Verfahren wegen Rundungsfehlern beim Rechnen mit endlicher Genauigkeit nur Näherungen. Deshalb ist die Unterscheidung vor allem für Entwicklung und Analyse wichtig.

Historisch gehen wichtige Grundverfahren auf Carl Friedrich Gauß zurück: das direkte gaußsche Eliminationsverfahren und das iterative Gauß-Seidel-Verfahren. Weitere zentrale Entwicklungen sind die Cholesky-Zerlegung 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 Krylow-Unterraum-Verfahren.

Im 20. Jahrhundert wurden Matrixzerlegungen zu Standardmethoden für mäßig große Systeme. Die Cholesky-Zerlegung wurde 1923 veröffentlicht. Für Eigenwertprobleme entstanden ab 1929 Vektoriterationsverfahren wie die Potenzmethode von Richard von Mises und 1944 die inverse Iteration von Helmut Wielandt. Der Durchbruch für die Berechnung aller Eigenwerte nicht zu großer Matrizen war 1961–1962 das QR-Verfahren. Für sehr große dünnbesetzte Matrizen wurde 1952 das CG-Verfahren entwickelt. Daraus entstand die wichtige Klasse der Krylow-Unterraum-Verfahren, darunter BiCG, MINRES, GMRES, QMR und TFQMR.

Struktur, Fehler und Stabilität

Numerische lineare Algebra nutzt gezielt Strukturen von Matrizen. Bei sehr großen Problemen kann allein das Speichern schwierig sein: Eine Matrix mit einer Million Zeilen und Spalten benötigt im double-precision-Format 8 Terabyte Speicherplatz. Viele Anwendungsprobleme führen aber auf dünnbesetzte Matrizen, bei denen pro Zeile nur wenige Einträge ungleich null sind. Verfahren, die Matrizen nur in Matrix-Vektor-Produkten verwenden, sind dafür besonders geeignet. Verfahren, die die Matrix selbst umformen, können die Dünnbesetztheit dagegen verlieren.

Besonders einfache Strukturen sind Diagonalmatrizen und Dreiecksmatrizen. Bei Diagonalmatrizen löst man lineare Gleichungssysteme durch n Divisionen; ihre Eigenwerte sind die Diagonalelemente. Dreieckssysteme werden durch Vorwärts- oder Rückwärtseinsetzen gelöst. Weitere wichtige Formen sind Bandmatrizen, Hessenbergmatrizen und Tridiagonalmatrizen. Symmetrische Matrizen sind ebenfalls bedeutsam, weil sie viele Probleme vereinfachen und den Lösungsaufwand oft etwa halbieren.

Fehler werden mit Normen gemessen. Für Vektoren ist die euklidische Norm ||x||2 = sqrt(sum{i=1}^n x_i^2) besonders verbreitet. Der absolute Fehler einer Näherung x~ zu x ist ||x~ - x||. Der relative Fehler ist ||x~ - x|| / ||x|| und ist das Standardmaß, weil er sich bei Skalierung nicht ändert. Matrixnormen messen entsprechend die Größe von Matrizen; natürliche Matrixnormen passen zu Vektornormen und erfüllen ||Ax|| ≤ ||A|| ||x||.

Die Kondition beschreibt, wie stark Datenfehler die Lösung beeinflussen. Für A x = b gilt bei Störung der rechten Seite die Abschätzung ||x~ - x|| / ||x|| ≤ ||A|| ||A^-1|| · ||b~ - b|| / ||b||. Die Zahl κ(A)=||A|| ||A^-1|| heißt Konditionszahl. Kleine Konditionszahlen bedeuten gute Kondition, große schlechte Kondition. Stabilität ist dagegen eine Eigenschaft eines Algorithmus. Vorwärtsstabilität bedeutet, dass das Ergebnis nicht wesentlich stärker vom exakten Ergebnis abweicht, als aufgrund von Datenfehlern und Kondition zu erwarten ist. Rückwärtsstabilität bedeutet, dass das berechnete Ergebnis als exakte Lösung eines nur wenig gestörten Eingabeproblems verstanden werden kann; ein rückwärtsstabiler Algorithmus ist auch vorwärtsstabil.

Orthogonalität und Transformationen

Orthogonale Matrizen sind Matrizen, deren Spalten eine Orthonormalbasis bilden, also paarweise senkrecht stehen und euklidische Länge 1 haben. Für eine orthogonale Matrix Q gilt Q^T Q = I und Q^-1 = Q^T. Dadurch lässt sich Qx = b einfach lösen: x = Q^T b. Außerdem erhalten orthogonale Matrizen die euklidische Norm, also ||Qx||_2 = ||x||_2. Daher gilt ||Q||_2 = 1 und κ(Q)=1. Multiplikation mit orthogonalen Matrizen vergrößert relative Fehler nicht.

Orthogonale Matrizen sind besonders wichtig bei Eigenwertproblemen. Symmetrische Matrizen lassen sich nach dem Spektralsatz orthogonal diagonalisieren: Für A^T = A existieren eine orthogonale Matrix Q und eine Diagonalmatrix D mit Q^T A Q = D. Die Diagonaleinträge von D sind die Eigenwerte, die Spalten von Q bilden eine Orthonormalbasis aus Eigenvektoren.

Zwei häufig verwendete orthogonale Transformationen sind Householder-Matrizen und Givens-Rotationen. Householder-Matrizen haben die Form H = I - 2vv^T mit ||v||_2 = 1 und beschreiben Spiegelungen. Sie können einen gegebenen Vektor a auf ein Vielfaches von e_1 transformieren, also H a = σ e_1 mit σ = ±||a||_2, und damit alle Einträge außer dem ersten zu null machen. Givens-Rotationen drehen nur eine zweidimensionale Ebene und lassen die übrigen n-2 Dimensionen unverändert; sie können gezielt einzelne Matrixeinträge auf null setzen.

Für Eigenwertberechnungen sind Ähnlichkeitstransformationen zentral. Zwei Matrizen A und B heißen ähnlich, wenn B = S^-1 A S mit regulärer Matrix S gilt. Ähnliche Matrizen haben dieselben Eigenwerte, und Eigenvektoren lassen sich ineinander umrechnen. Das charakteristische Polynom χ_A(λ)=det(λI-A) ist theoretisch wichtig, aber numerisch oft ungeeignet, weil die Nullstellenberechnung schlecht konditioniert sein kann. Daher formen viele Verfahren A durch Ähnlichkeitstransformationen in einfachere Matrizen um, meist mit orthogonalen Transformationsmatrizen.

Faktorisierungen und klassische Iterationen

Das gaußsche Eliminationsverfahren löst lineare Gleichungssysteme, indem es Variablen systematisch eliminiert. Durch Subtraktion geeigneter Vielfacher von Gleichungen entsteht eine Stufenform, die durch Rückwärtseinsetzen gelöst wird. Numerisch wichtig ist die Pivotisierung: Zeilen werden vertauscht, damit die Quotienten l_ik = a_ik / a_kk klein bleiben und Instabilitäten durch Stellenauslöschung vermieden werden.

Faktorisierungsverfahren zerlegen die Koeffizientenmatrix A in Produkte einfacher Matrizen, etwa A = BC. Dann löst man nacheinander B y = b und C x = y. Der Vorteil ist besonders groß, wenn mehrere Systeme mit derselben Matrix A, aber verschiedenen rechten Seiten gelöst werden: Die aufwendige Faktorisierung muss nur einmal berechnet werden.

Bei der LR-Zerlegung wird A ohne Zeilenvertauschungen als A = LR geschrieben, mit unterer Dreiecksmatrix L und oberer Dreiecksmatrix R. Mit Pivotisierung lautet die Form PA = LR, wobei P eine Permutationsmatrix ist. Danach löst man L y = P b durch Vorwärtseinsetzen und R x = y durch Rückwärtseinsetzen. Mit geeigneter Pivotisierung ist die Methode in praktischen Anwendungen fast immer stabil, obwohl pathologische Beispiele mit exponentiellem Fehlerwachstum existieren.

Die Cholesky-Zerlegung gilt für symmetrische positiv definite Matrizen, also Matrizen mit A^T = A und nur positiven Eigenwerten. Dann gibt es eine untere Dreiecksmatrix L mit A = LL^T. Wegen der Symmetrie ist der Aufwand etwa halb so groß wie bei der LR-Zerlegung. Sie kann bei Normalgleichungen A^T A x = A^T b für lineare Ausgleichsprobleme eingesetzt werden, ist aber nur für gut konditionierte Probleme mit wenigen Unbekannten empfehlenswert.

Die QR-Zerlegung schreibt A = QR mit orthogonaler Matrix Q und oberer Dreiecksmatrix R. Dann wird R x = y mit y = Q^T b gelöst. Wegen der guten Kondition orthogonaler Matrizen vermeidet sie Instabilitäten der LR-Zerlegung, kostet aber meist etwa doppelt so viel Rechenaufwand. Sie ist ein gängiges Verfahren für nicht zu große, gut konditionierte lineare Ausgleichsprobleme.

Fixpunktiteration mit Splitting-Verfahren beginnt mit einem Startvektor x_0 und berechnet x_{k+1}=M x_k + c. Solche Verfahren konvergieren genau dann, wenn alle Eigenwerte von M Betrag kleiner als 1 haben. Bei Splitting-Verfahren zerlegt man A = B + C mit leicht invertierbarem B und erhält x = -B^-1 C x + B^-1 b. Das Jacobi-Verfahren nutzt als B die Diagonalmatrix von A, das Gauß-Seidel-Verfahren den unteren Dreiecksanteil. Mit Relaxation, also x_{k+1}=x_k+ω Δx_k, entsteht zum Beispiel aus Gauß-Seidel das SOR-Verfahren.

Eigenwerte und große dünnbesetzte Systeme

Für symmetrische Matrizen ist das Jacobi-Verfahren ein einfaches, zuverlässiges iteratives Eigenwertverfahren. Es verwendet sukzessive Ähnlichkeitstransformationen mit Givens-Rotationen, auch Jacobi-Rotationen genannt, sodass eine Folge ähnlicher symmetrischer Matrizen gegen eine Diagonalmatrix D konvergiert. Die Diagonaleinträge von D liefern Näherungen für die Eigenwerte. Das klassische Verfahren wählt jeweils das größte Nichtdiagonalelement, das zyklische Verfahren durchläuft die Positionen in fester Reihenfolge. Beide konvergieren, aber relativ langsam; für dünnbesetzte Matrizen ist das Verfahren ungeeignet, weil die Matrix während der Iteration mit Nichtnulleinträgen aufgefüllt wird.

Vektoriterationsverfahren berechnen einzelne Eigenwerte und Eigenvektoren. Bei der Potenzmethode wird ein Startvektor x_0 ≠ 0 wiederholt mit A multipliziert, also x_k = A^k x_0, und zwischendurch normiert. Unter bestimmten Voraussetzungen konvergiert dies gegen einen Eigenvektor zum betragsgrößten Eigenwert. Die inverse Vektoriteration löst in jedem Schritt A x_{k+1}=x_k und zielt auf den betragskleinsten Eigenwert. Mit Shiftparameter σ verwendet man (A - σI) x_{k+1}=x_k, um einen Eigenvektor zu dem Eigenwert zu finden, der σ am nächsten liegt. Der zugehörige Eigenwert kann mit dem Rayleigh-Quotienten berechnet werden.

Das QR-Verfahren ist der wichtigste Algorithmus zur Berechnung aller Eigenwerte und Eigenvektoren nicht zu großer vollbesetzter Matrizen. Es startet mit A_0=A, berechnet A_k=Q_k R_k und setzt A_{k+1}=R_k Q_k. Wegen Q_k^-1=Q_k^T ist dies eine orthogonale Ähnlichkeitstransformation: A_{k+1}=Q_k^T A_k Q_k. Praktisch wird A vorher auf Hessenberg-Gestalt gebracht; bei symmetrischen Matrizen entsteht eine Tridiagonalmatrix. Shiftstrategien mit A_k - σ_k I beschleunigen die Konvergenz. Eine Variante berechnet auch die Singulärwertzerlegung, die etwa in der Bildkompression und bei schlecht konditionierten Ausgleichsproblemen genutzt wird.

Krylow-Unterraum-Verfahren sind besonders wichtig für sehr große dünnbesetzte lineare Gleichungssysteme und Eigenwertprobleme. Beim CG-Verfahren für symmetrische positiv definite Matrizen wird A x = b als Minimierungsproblem für E(x)=1/2 x^T A x - b^T x verstanden. Statt steilsten Abstiegsrichtungen nutzt CG konjugierte Suchrichtungen mit s_k^T A s_l = 0 und erreicht theoretisch nach n Abstiegen die exakte Lösung, oft aber schon früher eine gute Näherung.

Die Näherungen liegen in Krylow-Unterräumen K_k(A,b)=span(b, Ab, A^2b, ..., A^{k-1}b). Allgemeine Krylow-Verfahren wählen x_k in diesem Raum so, dass das Residuum r_k=A x_k-b nach einer Projektionsbedingung klein wird. Weitere Verfahren sind BiCG, CGS und BiCGSTAB, außerdem die meist stabilen Verfahren GMRES und MINRES, die ||A x_k - b||_2 minimieren, sowie QMR und TFQMR. Für Eigenwertprobleme führen Krylow-Ideen auf das Arnoldi-Verfahren und im symmetrischen Fall auf das Lanczos-Verfahren.