Numerische Lineare Algebra: LR-Zerlegung, Cholesky, Jacobi-Verfahren

Auf der Seite Gauß-Verfahren löst du Gleichungssysteme mit Zeilenumformungen. Für den Computer ist dasselbe Verfahren die Grundlage fast aller Rechnungen in Technik und Wissenschaft: Strukturberechnung, Wettermodelle, Regression, Bildverarbeitung. Dort geht es um Systeme mit $n={10}^3$ bis ${10}^9$ Unbekannten, und zwei Fragen entscheiden: Wie viel Rechenarbeit braucht das Verfahren, und wie verlässlich ist das Ergebnis trotz Rundungsfehlern?

Die Antwort ist eine geschickte Umformulierung: Man zerlegt die Matrix einmal in zwei Dreiecksmatrizen und löst danach beliebig viele rechte Seiten billig. Für symmetrische positiv definite Matrizen geht es mit der Cholesky-Zerlegung doppelt so schnell. Für sehr große dünnbesetzte Systeme nimmt man stattdessen Iterationsverfahren. Dazu brauchst du die Begriffe aus Numerik-Grundlagen (Rundungsfehler, Kondition), das Lösen linearer Systeme und für die Konvergenz den Banachschen Fixpunktsatz aus Metrische Räume und Topologie.

Dreieckssysteme

Ein System mit Dreiecksmatrix lässt sich ohne Elimination lösen. Bei der oberen Dreiecksmatrix $R$ löst man von unten nach oben (Rückwärtseinsetzen): $x_n=\frac{y_n}{r_{nn}}$ und $x_i=\frac{1}{r_{ii}}(y_i-\sum _{k=i+1}^n{r}_{ik}{x}_k)$. Bei der unteren Dreiecksmatrix $L$ geht es von oben nach unten (Vorwärtseinsetzen): $y_i=\frac{1}{l_{ii}}(b_i-\sum _{k=1}^{i-1}{l}_{ik}{y}_k)$. Beides kostet nur etwa $n^2$ Operationen, denn für die $i$-te Zeile braucht man $n-i$ beziehungsweise $i-1$ Multiplikationen.

Die LR-Zerlegung

Das Gauß-Verfahren erzeugt aus $A$ eine obere Dreiecksmatrix $R$. Merkt man sich die Eliminationsfaktoren $l_{ik}=\frac{a_{ik}}{a_{kk}}$, mit denen man das $l_{ik}$-Fache von Zeile $k$ von Zeile $i$ abgezogen hat, in einer Matrix $L$, so gilt $A=LR$.

Beweisidee (Induktion). Schreibe $A=\left(\begin{array}{cc}{a}_{11}& r^{\top }\\ c& A^{\prime }\end{array}\right)$ mit $a_{11}\ne 0$. Dann ist

Die verbleibende Matrix $S=A^{\prime }-\frac{c{r}^{\top }}{a_{11}}$ (das Schur-Komplement) ist eine $(n-1)\times (n-1)$-Matrix, auf die man dasselbe anwendet. Die Voraussetzung über die Hauptminoren stellt sicher, dass jedes Pivotelement ungleich 0 ist. Die Eindeutigkeit folgt so: Aus $L_1{R}_1=L_2{R}_2$ ergibt sich $L_2^{-1}{L}_1=R_2{R}_1^{-1}$. Links steht eine untere Dreiecksmatrix mit Einsen auf der Diagonale, rechts eine obere Dreiecksmatrix. Beides zugleich ist nur die Einheitsmatrix, also ist $L_1=L_2$ und $R_1=R_2$.

Aufwand. Im Schritt $k$ werden für $n-k$ Zeilen je $n-k$ Multiplikationen und Subtraktionen gebraucht. Das sind $\sum _{k=1}^{n-1}{(n-k)}^2\approx \frac{n^3}{3}$ Multiplikationen und ebenso viele Subtraktionen, zusammen etwa $\frac{2}{3}{n}^3$ Operationen. Für $n=1000$ sind das gut $\mathrm{6,7}·{10}^8$ Operationen, auf einem Rechner Sekundenbruchteile. Das Berechnen der Inversen kostet etwa das Dreifache, und das Auswerten der Determinante nach der Leibniz-Formel mit $n!$ Summanden wäre für $n=20$ schon unmöglich.

Pivotisierung

Ist ein Pivotelement $a_{kk}$ klein, ist der Faktor $l_{ik}=\frac{a_{ik}}{a_{kk}}$ riesig, und die Rundungsfehler wachsen ungebremst.

Die Spaltenpivotsuche vertauscht in jedem Schritt die Zeile so, dass das betragsgrößte Element der Spalte Pivotelement wird. Dann gilt $\left|l_{ik}\right|\le 1$. Das Ergebnis ist $PA=LR$. Mit dieser Strategie ist das Gaußsche Verfahren in der Praxis rückwärtsstabil: Die berechnete Lösung ist die exakte Lösung eines nur wenig gestörten Systems. Das hat James Hardy Wilkinson in den 1960er-Jahren gezeigt. Die Genauigkeit des Ergebnisses hängt danach nur noch von der Konditionszahl ab, siehe Numerik-Grundlagen. Das Residuum $r=b-A\tilde{x}$ ist ein guter Indikator: Es gilt $\frac{\Vert x-\tilde{x}\Vert }{\Vert x\Vert }\le \kappa (A)\frac{\Vert r\Vert }{\Vert b\Vert }$.

Die Cholesky-Zerlegung

Für symmetrische, positiv definite Matrizen (siehe Definitheit), wie sie bei Normalgleichungen und Energieminimierung entstehen, geht es schneller und ohne Pivotisierung.

Beweisidee. Eine positiv definite Matrix hat nur positive Hauptminoren, also existiert $A=LR$ ohne Pivotisierung, und die Pivotelemente $r_{ii}$ sind positiv. Aus der Symmetrie folgt $R=D{L}^{\top }$ mit der Diagonalmatrix $D=\mathrm{diag}(r_{ii})$. Dann ist $A=LD{L}^{\top }=(L{D}^{1/2}){(L{D}^{1/2})}^{\top }$, also $G=L{D}^{1/2}$. Die Wurzeln sind wegen $r_{ii}>0$ reell. Umgekehrt ist das Verfahren ein Test auf positive Definitheit: Scheitert es an einer negativen Zahl unter der Wurzel, ist $A$ nicht positiv definit.

Der Aufwand ist etwa $\frac{1}{3}{n}^3$ Operationen, halb so viel wie bei der LR-Zerlegung, weil nur die halbe Matrix gespeichert und berechnet wird. Das Verfahren ist ohne Pivotisierung stabil.

Iterative Verfahren

Bei sehr großen Systemen, die zudem dünnbesetzt sind (die meisten Einträge sind 0, wie bei der Diskretisierung von Differentialgleichungen), füllt die LR-Zerlegung die Nullen mit Zahlen auf und sprengt Speicher und Rechenzeit. Iterationsverfahren erzeugen stattdessen eine Folge $x^{(0)},x^{(1)},\dots$ von Näherungen, bei der jeder Schritt nur eine Matrix-Vektor-Multiplikation mit der dünnen Matrix braucht.

Man zerlegt $A=L+D+U$ in den strikt unteren Teil $L$, die Diagonale $D$ und den strikt oberen Teil $U$ (diese Matrizen haben nichts mit den Faktoren der LR-Zerlegung zu tun) und löst die $i$-te Gleichung nach $x_i$ auf:

Beim Gauß-Seidel-Verfahren benutzt man in dieser Formel für $j<i$ bereits die neuen Werte $x_j^{(k+1)}$ dieses Schritts, die schon berechnet sind. Das konvergiert meist doppelt so schnell und braucht nur einen Vektor Speicher.

Beweis für Jacobi. Die Iteration ist $x^{(k+1)}=B{x}^{(k)}+c$ mit $B=-D^{-1}(L+U)$. Die Zeilensummennorm ist $\Vert B{\Vert }_{\infty }=\underset{i}{\max }\sum _{j\ne i}\frac{\left|a_{ij}\right|}{\left|a_{ii}\right|}=:q<1$. Damit ist $x\mapsto Bx+c$ eine Kontraktion im vollständigen Raum $({ℝ}^n,\Vert ·{\Vert }_{\infty })$, und der Banachsche Fixpunktsatz liefert Existenz und Konvergenz, mit der Fehlerschranke $\Vert x^{(k)}-x{\Vert }_{\infty }\le \frac{q^k}{1-q}\Vert x^{(1)}-x^{(0)}{\Vert }_{\infty }$.

Moderne Verfahren wie das CG-Verfahren (konjugierte Gradienten, für symmetrisch positiv definite Matrizen) und GMRES sind viel schneller, bauen aber auf derselben Idee auf: Näherungen aus Matrix-Vektor-Produkten. Für Ausgleichsprobleme $\min \Vert Ax-b{\Vert }_2$ nutzt man die QR-Zerlegung (in der Praxis mit Householder-Spiegelungen, im Prinzip mit dem Gram-Schmidt-Verfahren der Seite Skalarprodukt und Orthogonalität), weil sie die Kondition nicht quadriert wie die Normalgleichungen $A^{\top }Ax=A^{\top }b$.

Typische Fehler

Übungen

Geschichte: von den Neun Kapiteln zum Computer

Das Jiu Zhang Suanshu (Neun Kapitel der Rechenkunst), entstanden zwischen etwa 200 v. Chr. und 100 n. Chr., enthält eine klare Demonstration des Eliminationsverfahrens, und Liu Hui kommentierte es 263. In Europa beschrieb Isaac Newton die Elimination in seinen Algebra-Vorlesungen, die 1707 gedruckt wurden. Joseph-Louis Lagrange führte 1759 für quadratische Formen eine Umformung ein, die einer Zerlegung $A=U^{\top }DU$ mit einer oberen Dreiecksmatrix $U$ und einer Diagonalmatrix $D$ entspricht. Gauß erwähnte 1809 die „gewöhnliche Elimination“ für die Normalgleichungen der Ausgleichsrechnung und entwickelte in seiner Arbeit über die Bahn des Asteroiden Pallas von 1810 eine Rechenschreibweise dafür. Diese Schreibweise übernahmen später berufsmäßige Rechner, deshalb trägt das Verfahren heute seinen Namen. Matrizen benutzte man dabei noch nicht. Die Deutung als Zerlegung in Dreiecksmatrizen kam erst im 20. Jahrhundert. Das iterative Gauß-Seidel-Verfahren erwähnte Gauß 1823 in einem Brief an Gerling. Veröffentlicht wurde es von Seidel 1874. Tadeusz Banachiewicz beschrieb die Zerlegung 1938 als Matrixprodukt. John von Neumann und Alan Turing untersuchten die Rundungsfehler, Turing 1948, und gaben dem Verfahren seine heutige Form. James Hardy Wilkinson zeigte in den 1960er-Jahren, dass das Verfahren mit Pivotisierung rückwärtsstabil ist.

Häufige Fragen

Quellen und weiterführende Literatur

• Wikipedia: LR-Zerlegung
• Wikipedia: Cholesky-Zerlegung
• Wikipedia: Gauß-Seidel-Verfahren
• Peter Deuflhard, Andreas Hohmann: Numerische Mathematik 1, De Gruyter, Berlin.
• Lloyd N. Trefethen, David Bau: Numerical Linear Algebra, SIAM, Philadelphia.
• Alan M. Turing: Rounding-off errors in matrix processes, Quarterly Journal of Mechanics and Applied Mathematics 1 (1948).
• Joseph F. Grcar: How ordinary elimination became Gaussian elimination, Historia Mathematica 38 (2011), S. 163 bis 218 (englisch)