Lösung linearer Gleichungssysteme

Original post is here: eklausmeier.goip.de

• 1. Konditionszahlen von Matrizen
• 2. Elementare Zeilen- und Spaltenoperationen
• 3. Die $LU$-Zerlegung
• 4. Die Gauß-Elimination

Bei jedem Iterationsschritt eines Newton-Raphson-Verfahrens, bzw. bei jeder Aktualisierung der Iterationsmatrix beim Newton-Kantorovich Iterationsverfahren, fällt die Lösung eines linearen Gleichungssystems an. Es ist daher von Wichtigkeit, für diesen Prozeß ein zuverlässiges und zugleich effizientes Lösungsverfahren parat zu haben.

1. Konditionszahlen von Matrizen #

Wie könnte man erkennen, ob ein berechnetes Ergebnis $x$ für ein lineares Gleichungssystem $Ax=b$, eine gute Näherung darstellt? Naheliegend wäre die Betrachtung des Residuums $Ax-b$. Wie das folgende Beispiel deutlich macht, ist dieses Maß mit Vorsicht zu geniessen.

1. Beispiel: Sei betrachtet das lineare Gleichungssystem $$ \pmatrix{10&7&8&7\cr 7&5&6&5\cr 8&6&1& 9\cr 7&5&9&10\cr} \pmatrix{x_1\cr x_2\cr x_3\cr x_4\cr}= \pmatrix{32\cr 23\cr 33\cr 31\cr}. $$ Einige Vektoren $x$ und ihre Bilder $Ax$ lauten $$ \pmatrix{6\cr -7.2\cr 2.9\cr -0.1\cr}\mapsto \pmatrix{32.1\cr 22.9\cr 32.9\cr 31.1\cr},\qquad \pmatrix{1.50\cr 0.18\cr 1.19\cr 0.89\cr}\mapsto \pmatrix{32.01\cr 22.99\cr 32.99\cr 31.01\cr},\qquad \hbox{richtig:}\quad\pmatrix{1\cr 1\cr 1\cr 1\cr}\mapsto \pmatrix{32\cr 23\cr 33\cr 31\cr}. $$ Die ersten beiden Vektoren könnten den Eindruck erwecken, daß man sich schon sehr nahe der richtigen Lösung befindet, jedoch ist beim ersten Mal keine der angegebenen Nachkommaziffern richtig, und beim zweiten Mal ist nur bei 2 Komponenten eine einzige Stelle hinter dem Komma richtig. Dies alles ist möglich, obwohl die Matrix symmetrisch ist und alle Komponenten der Matrix dem gleichen Größenbereich entstammen.

Ähnliche Verhältnisse liegen vor bei dem nachstehenden Beispiel.

2. Beispiel: Betrachtet sei nun das lediglich zweidimensionale Problem mit den Daten $$ A:=\pmatrix{1.2969&0.8648\cr 0.2161&0.1441\cr},\qquad b:=\pmatrix{0.8642\cr 0.1440\cr},\qquad \overline x:=\pmatrix{\phantom{-}0.9911\cr -0.4870\cr}. $$ Der Residuenvektor ist hier exakt $(-10^{-8},{\mskip 3mu}10^{-8})$, dennoch lautet die exakte Lösung $(2,{\mskip 3mu}-2)$.

Die Empfindlichkeit, mit der ein lineares Gleichungsystem auf Änderung der Eingabedaten reagiert, wird quantitativ erfaßt durch die Konditionszahl.

3. Satz: Voraussetzungen: Die Matrix $A$ sei invertierbar und $x$ sei die exakte Lösung des linearen Gleichungssystems $Ax=b$. Bekannt sei eine Näherungslösung $\hat x$. $r$ sei das Residuum, also $r:=b-A\hat x$.

Behauptung: Es gilt die beidseitige, scharfe Abschätzung $$ {1\over|A|\cdot|A^{-1}|}{|r|\over|b|}\buildrel 1)\over\le {|x-\hat x|\over|x|}\buildrel 2)\over\le \left|A\right| \cdot \left|A^{-1}\right| {\left|r\right|\over\left|b\right|}. $$

Beweis: Zu 1). Wegen $|r|=|b-A\hat x|=|Ax-A\hat x|\le|A|\cdot|x-\hat x|$, ist $|r|/|A|\le|x-\hat x|$. Andererseits ist $|x|=|A^{-1}b|\le|A^{-1}|\cdot|b|$, also $1/(|A^{-1}|\cdot|b|)\le1/|x|$, nach Multiplikation entsprechender Seiten der beiden hergeleiteten Ungleichung die unter 1) gemachte Behauptung $$ {1\over|A|\cdot|A^{-1}|}{|r|\over|b|}\le{|x-\hat x|\over|x|}. $$ Diese Abschätzung ist scharf. Sei zur Abkürzung $e:=x-\hat x$. Die Ungleichung wird zu einer Gleichung, wenn $|Ae|=|A|\cdot|e|$ und $|A^{-1}b|=|A^{-1}|\cdot|b|$. Nach Definition von Matrixnormen gibt es solche Vektoren $e^*$ und $b^*$, sodaß $|Ae^*|=|A|\cdot|e^*|$ und $|A^{-1}b^*|=|A^{-1}|\cdot|b^*|$. Für $b^*$ und $\tilde x:=A^{-1}b^*-e$ wird die Ungleichung zu einer Gleichung.

Zu 2). Erstens wegen $\left|b\right| = \left|Ax\right| \le \left|A\right| \cdot \left|x\right|$ ist $1 / \left|x\right| \le \left|A\right| / \left|b\right|$. Zweitens ist $\left|x-\hat x\right| = \left|A^{-1}b-A^{-1}A\hat x\right| = \left|A^{-1}r\right| \le \left|A^{-1}\right| \cdot \left|r\right|$. Insgesamt ergibt sich wieder durch Multiplikation entsprechender Seiten $$ {|x-\hat x|\over|x|}\le|A|\cdot|A^{-1}|\cdot{|r|\over|b|}. $$ Auch diese Ungleichung ist scharf. Wähle $x^*$ und $r^*$ mit $|Ax^*|=|A|\cdot|x^*|$ und $\left|A^{-1}r^*\right| = \left|A^{-1}\right| \cdot \left|r^*\right|$. Für $b^*:=A^{-1}x^*$ und $\hat x^*:=x^*-A^{-1}r^*$ gilt Gleichheit.     ☐

4. Definition: Die Zahl $\kappa(A) := \left|A\right| \cdot \left|A^{-1}\right|$ heißt Konditionszahl der quadratischen Matrix $A$ zur Norm $\left|\cdot\right|$.

Für rein theoretische Überlegungen stellt die Konditionszahl eine gute Beschreibung von Störungs- und Empfindlichkeitsphänomenen dar. Jedoch ist man ja gerade an der Inverse, bzw. an der Lösung des linearen Gleichungssystems interessiert. Die hier bei der Analyse auftauchende Konditionszahl ist also bei der praktischen Auflösung nicht bekannt. Es gibt Verfahren, mit denen die Konditionszahl geschätzt werden kann. Diese Verfahren sind stellenweise mit der Gaußschen Eliminationsmethode auf das engste gekoppelt und ergeben sich daher zusammen mit der Rechnung. Direkt aus der Defintion der Konditionszahl ersieht man

5. Eigenschaften: Es gelten $\kappa(A)\ge1$ und $\kappa(AB)\le\kappa(A)\cdot\kappa(B)$.

Mit Hilfe der Konditionszahl lässt sich eine Abschätzung einfach schreiben, welche charakterisiert, inwieweit Störungen in der Matrix $A$ zu Veränderungen in der eigentlich gewünschten Lösung $x$ führen.

6. Satz: Löst man statt des exakten Systems $Ax=b$ das gestörte System $\hat A\hat x=b$, so gilt die scharfe Abschätzung nach oben $$ {\left|x-\hat x\right|\over\left|\hat x\right|} \le \kappa(A) {|A-\hat A|\over\left|A\right|}. $$

Beweis: Zunächst ist $x=A^{-1}b=A^{-1}\hat A\hat x=A^{-1}(A+\hat A-A)\hat x=\hat x+A^{-1}(\hat A-A)\hat x$, also $x-\hat x=A^{-1}(\hat A-A)\hat x$, somit $|x-\hat x|\le|A^{-1}|\cdot|A-\hat A|\cdot|\hat x|$. Durch Erweitern mit $|A|\ne0$ schließlich $$ {\left|x-\hat x\right|\over\left|\hat x\right|} \le \left|A^{-1}\right| \cdot |A-\hat A| %\left|A-\hat A\right| = \left|A^{-1}\right| \cdot \left|A\right| % {\left|A-\hat A\right|\over\left|A\right|}. {|A-\hat A|\over\left|A\right|}. $$ Die Abschätzung ist offensichtlich scharf.     ☐

Die Konditionszahl charakterisiert gleichzeitig auch den Abstand zu allen nicht-invertierbaren Matrizen “in der Nähe” von $A$.

7. Satz: Für alle invertierbaren Matrizen $A$ gilt $$ \min\left{{|A-B|\over|A|}:\hbox{$B$ nicht invertierbar}\right} \ge{1\over\kappa(A)}. $$ Für die Maximumnorm, die 1-Norm und die euklidische Norm gilt Gleichheit.

Beweis: Ist $B$ eine beliebige nicht-invertierbare Matrix, so gibt es also einen Vektor $x\ne0$, mit $Bx=0$. Für diesen Vektor $x$ gilt: $$ \eqalign{ |x| = |A^{-1}Ax| &\le |A^{-1}|\cdot|Ax|=|A^{-1}|\cdot|(A-B)x|\cr &\le |A^{-1}|\cdot|A-B|\cdot|x|.\cr } $$ Daher $1\le|A^{-1}|\cdot|A-B|$, somit $$ {1\over\kappa(A)}={1\over|A|\cdot|A^{-1}|}\le{|A-B|\over|A|}. $$ Für die behauptete Gleichheit bei der Maximumnorm, schätzt man in umgekehrter Richtung ab und zeigt damit indirekt durch Eingabelung, die Gleichheit.

Zur genaueren Unterscheidung von Normen und Betragsstrichen, werde mit $\left|\cdot\right|_\infty$ die Maximumnorm bezeichnet und mit $\left|\cdot\right|$ die gewöhnliche Betragsfunktion für skalare Größen. Sei $v$ ein Einheitsvektor mit den beiden Eigenschaften $\left|v\right|_\infty = 1$ und $\left|A^{-1}v\right|_\infty = \left|A^{-1}\right|_\infty \cdot \left|v\right|_\infty$. Nach Definition der zu einer Vektornorm gehörenden Matrixnorm gibt es solch einen Vektor $v$.

Weiter sei $y:=A^{-1}v$, $\left|y\right|_\infty =: \left|y_m\right|$. Sei $e_m$ der $m$-te Einheitsvektor; $z:=y_m^{-1}e_m$ und $B:=A-vz^\top$. Wegen $By=AA^{-1}v-y_m^{-1}y_mv=0$ und $y\ne0$, ist $B$ nicht invertierbar. Für beliebige Vektoren $x$ gilt $$ \eqalign { \left|(A-B)x\right|\infty &= \left|y_m^{-1}x_mv\right|\infty = \left|y_m\right|^{-1} \cdot \left|x_m\right| \cdot \left|v\right|\infty\cr &= \left|y_m\right|^{-1} \left|x_m\right|\cr &= \left|y\right|\infty \cdot \left|x_m\right|\cr &= (\left|A^{-1}v\right|\infty)^{-1} \cdot \left|x_m\right|\cr &= \left|A^{-1}\right|\infty^{-1} \cdot \left|v\right|\infty \cdot \left|x_m\right|\cr &= \left|A^{-1}\right|\infty^{-1} \cdot \left|x\right|\infty.\cr } $$ Da $x$ beliebig war gilt somit $$ \left|A-B\right|\infty \le {1\over\left|A^{-1}\right|\infty}, $$ also $$ {\left|A-B\right|\infty\over\left|A\right|\infty} \le {1\over\left|A^{-1}\right|\infty \cdot \left|A\right|\infty} = {1\over\kappa\infty(A)}. $$ Für die $1$-Norm führt man den Beweis ganz ähnlich und für die euklidische Norm benutzt man den Satz über die Existenz einer Singulärwertzerlegung. Der genaue Beweis werde hier nicht ausgeführt.     ☐

In anderer Formulierung des oben schon bewiesenen Satzes, lässt sich schreiben:

8. Satz: Es sei $A$ invertierbar und es sei $A(x+\Delta x)=b+\Delta b$. Dann gilt $$ {\left|\Delta x\right|\over\left|x\right|} \le \kappa(A){\left|\Delta b\right|\over\left|b\right|}. $$

Der nächste Satz zeigt, daß das symmetrisierte Gleichungssystem $A^\top Ax=A^\top b$ eine größere, und damit schlechtere Konditionszahl besitzt, als das ursprüngliche System. Insbesondere gelten diese Überlegungen für das Normalgleichungssystem (Ausgleichung im Sinne kleinster Quadrate), obwohl dort die entsprechenden Matrizen nicht quadratisch, also erst recht nicht invertierbar sind.

9. Satz: Für eine beliebige invertierbare Matrix $A$ gilt $\kappa_s(A)\le\kappa_s(A^\top A)$.

Beweis: Seien $\mu_{\rm max}$ und $\mu_{\rm min}$ entsprechend die größten und kleinsten Eigenwerte von $A^\top A$. Dann ist $\left|A\right|_s=\sqrt{\mu_{\rm max}}$ und $\left|A^{-1}\right|=\sqrt{\mu_{\rm min}^{-1}}$. Weiter ist $\left|A^\top A\right|_s = \mu_{\rm max}$ und $\left|(A^\top A)^{-1}\right| = \mu_{\rm min}^{-1}$. Daher ist $$ \kappa_s = \sqrt{\mu_{\rm max}\over\mu_{\rm min}} \le {\mu_{\rm max}\over\mu_{\rm min}} = \kappa_s(A^\top A). $$     ☐

Welche ist nun die beste Vorkonditionierung mit einer Diagonalmatrix? Es zeigt sich nun, daß dies gerade die normmässige Äquilibrierung aller Zeilenvektoren der Matrix $A$ ist.

10. Satz: Ist die invertierbare $(n\times n)$-Matrix $A=(a_{ik})$ normiert (äquilibriert) gemäß $$ \sum_{k=1}^n \left|a_{ik}\right| = 1,\qquad\hbox{für}\quad i=1,\ldots,n, $$ so gilt für jede Diagonmalmatrix $D$ mit $\det D\ne0$, die Ungleichung $$ \kappa_\infty(DA)\ge\kappa_\infty(A). $$ $\kappa_\infty$ bezeichnet hier die Konditionszahl bezüglich der Zeilensummennorm (verträgliche Norm für die Maximumnorm $\left|\cdot\right|$).

Beweis: siehe Werner (1975)*1972+1, Werner, Helmut. Es sei $\left|\cdot\right|$ die Maximum-Vektornorm bzw. die Zeilensummennorm bei Matrizen. Für jede Diagonalmatrix $D=(d_{ii})$, mit $\det D\ne0$ gilt $$ \left|DA\right| = \max_{i=1}^n \left(\left|d_{ii}\right| \sum_{k=1}^n \left|a_{ik}\right| \right) = \max_{i=1}^n \left|d_{ii}\right| $$ und $$ \left|(DA)^{-1}\right| = \left|A^{-1}D^{-1}\right| = \max_{i=1}^n \left( \sum_{k=1}^n \left|\tilde a_{ik}\right| {1\over d_{kk}}\right) \ge \left|A^{-1}\right| \cdot \min_{i=1}^n {1\over \left|d_{ii}\right|}. $$ Hierbei bezeichnete $\tilde a_{ik}$ die Komponenten der inversen Matrix $A^{-1}$. Aus den beiden obigen Gleichungen folgt $$ \kappa_\infty(DA) = \left|DA\right| \cdot \left|(DA)^{-1}\right| \ge \left|A^{-1}\right| \cdot \underbrace{\max_{i=1}^n \left|d_{ii}\right| \cdot \min_{i=1}^n {1\over\left|d_{ii}\right|}} {\displaystyle{{} = 1 = \left|A\right|}} =\kappa\infty(A). $$     ☐

Versucht man nun hingegen auf beiden Seiten der Matrix eine Äquilibrierung zu erreichen, so kann man sich nach den Worten von Dahlquist und Björck u.U. ebenfalls “in die Nesseln setzen”.

11. Beispiel: Man betrachte für $0 < \left|\varepsilon\right| < 1$ die Matrix $$ A := \pmatrix{\varepsilon&-1&1\cr -1&1&1\cr 1&1&1\cr},\qquad A^{-1} = {1\over4}\pmatrix{0&-2&2\cr -2&1-\varepsilon&1+\varepsilon\cr 2&1+\varepsilon&1-\varepsilon\cr}. $$ Sein nun $D_1 := \mathop{\rm diag}(1,\varepsilon,\varepsilon)$ und $D_2 := \mathop{\rm diag}(\varepsilon^{-1},1,1)$. Die skalierte Matrix $$ B:=D_2AD_1=\pmatrix{1&-1&1\cr -1&\varepsilon&\varepsilon\cr 1&\varepsilon&\varepsilon\cr} $$ ist zwar jetzt zeilenäquilibriert, jedoch beträgt die Konditionszahl jetzt $\kappa(B)\approx3/\varepsilon$, während hingegen $\kappa(A)=3$.

2. Elementare Zeilen- und Spaltenoperationen #

Wenn auch die Cramersche Regel eine elegante Darstellung der Lösung eines linearen Gleichungssystems liefert, so ist doch selbst bei bestmöglicher Ausrechnung aller Determinanten der Aufwand höher, als derjenige von Verfahren, die im folgenden vorgestellt werden. Würde man die $(n+1)$ Determinanten ($n$ Zählerdeterminanten und eine Nennerdeterminante) als Summe von jeweils $n$ Faktoren berechnen, so gelänge man zu Rechengrößenordnungen der Form $(n!)$. Schon $n=50$ würde mehr als $10^{64}$ Gleitkommamultiplikationen erfordern, was selbst für Größtrechenanlagen unvertretbar lange Rechenzeiten heraufbeschwören würde. Aber, wie eben erwähnt, selbst bei bestmöglicher und effizientester Auswertung von Determinanten, wären immer noch größenordnungsmässig $n^4$ Operationen nötig, während hingegen die im weiteren Verlaufe dargestellten Verfahren in der Größenordnung $n^3$ liegen.

1. Ein lineares Gleichungssystem mit $p$ Gleichungen und $n$ Unbekannten $x_1$, $x_2$ bis $x_n$ hat die Form $$ \eqalign{ a_{11}x_1+a_{12}x_2+\cdots+a_{1n}x_n &= b_1\cr a_{21}x_1+a_{22}x_2+\cdots+a_{2n}x_n &= b_2\cr &>\vdots\cr a_{p1}x_1+a_{p2}x_2+\cdots+a_{pn}x_n &= b_p\cr } $$ Hierbei sind $a_{ij}$ und $b_k$ feste gegebene Zahlen. Die $a_{ij}$ heißen die Koeffizienten. Der erste Index (hier $i$) heißt Zeilenindex, der zweite Index (hier $j$) heißt Spaltenindex. Der Vektor $(b_1,\ldots,b_p)$ heißt Vektor der rechten Seite [^{right hand side (RHS)}]. Das lineare Gleichungssystem heißt homogen, falls der Vektor der rechten Seite gleich dem Nullvektor ist, also $b_k=0$, für alle $k=1,\ldots,p$.

2. Daß man beim Operieren mit Gleichungssystemen vorsichtig sein muß, zeigt das folgende Beispielgleichungssystem $$ \eqalign{ x+y+z &= 1\cr x-y+z &= 0\cr -x+y-z &= 0 } $$ Zieht man nun die drei Folgerungen, erstens die 1.te Gleichung beizubehalten, zweitens alle drei Gleichungen zusammenzuaddieren und drittens die 2.te zur 3.ten Gleichung zu addieren, so erhält man $$ \eqalign{ x+y+z &= 1\cr x+y+z &= 1\cr 0 &= 0\cr } $$ Aufgrund der Konstruktion ist jedes Tripel $(x,y,z)$, welches das ursprüngliche Gleichungssystem löst auch gleichzeitig Lösung des neuen umgeformten Systems. Die Umkehrung gilt jedoch nicht! Das Tripel mit $x=y=z=1/3$ löst zwar das neue, umgeformte System, nicht aber das ursprüngliche. Durch Ziehen von Folgerungen können also Lösungen hinzukommen!

3. Gibt es Umformungen, die die Lösungsmenge nicht verändern? Ja, es gibt bei linearen Gleichungssystemen unendlich viele Umformungsmöglichkeiten, die die Lösungsgesamtheit nicht verändern. Drei besonders wichtige sind die nachstehenden Umformungen:

1. Vertauschen zweier Gleichungen, und
2. Multiplikation einer der Gleichungen mit einer Zahl${}\ne0$,
3. Addition einer mit einer beliebigen Zahl multiplizierten Gleichung zu einer weiteren Gleichung.

Die nicht betroffenen Gleichungen des Systems werden beibehalten, so wie sie sind. Es gilt nun, daß die obigen drei Vertreter von Umformungen, die Lösungsgesamtheit nicht verändern. Diese drei ausgezeichneten Umformungen, heißen {\it ^{elementare Umformungen}}. Mit den obigen drei Umformungen, wäre das obige malheur nicht passiert.

4. Satz: Bei elementaren Umformungen ändert sich die Lösungsmenge nicht.

Beweis: Das oben angeschriebene Gleichungssystem lautet nach Anwendung der dritten Umformungsregel $$ \eqalign{ a_{11}x_1+a_{12}x_2+\cdots+a_{1n}x_n &= b_1\cr (a_{21}+\lambda a_{11})x_1+(a_{22}+\lambda a_{12})x_2+\cdots +(a_{2n}+\lambda a_{1n})x_n &= b_2+\lambda b_1\cr &>\vdots\cr a_{p1}x_1+a_{p2}x_2+\cdots+a_{pn}x_n &= b_p\cr } $$ Die erste Zeile des Gleichungssytems wurde mit $\lambda$ multipliziert und zur zweiten Gleichung hinzuaddiert. Die restlichen Zeilen wurden völlig unverändert übernommen. Die Lösung des alten Gleichungssystems ist zugleich auch Lösung des neuen umgeformten Lösungssystems, da aber diese Umformung rückgängig gemacht werden kann, hat das neue System genau die gleichen Lösungen. Die Rückgängigmachung geschähe dadurch, daß man die erste Gleichung mit $(-\lambda)$ multipliziert und zur zweiten Gleichung hinzuaddiert. Der Rest der Gleichungen wird wieder belassen.     &#9744;

5. Zu diesen drei elementaren Umformungen korrespondieren die folgenden sogenannten Elementarmatrizen. Die Vertauschung zweier Zeilen $$ \begin{pmatrix} & & & & & & & \cr & 1 & & & & & & 0\cr & & \ddots & & & & \unicode{x22F0} & \cr i\to & & & 0 & \ldots & 1 & \ldots & \cr & \vdots & & \vdots & & \vdots & & \vdots\cr j\to & & \ldots & 1 & \ldots & 0 & & \cr & \vdots & \unicode{x22F0} & & & & \ddots & \vdots\cr & 0 & & & & & & 1\cr \end{pmatrix} $$ Addiert man zur $i$-ten Zeile von $L(\lambda)$ die $j$-te Zeile multipliziert mit $f(\lambda)$, so ist die äquivaltent mit der Linksmultiplikation mit $$ \begin{pmatrix} & & & & & j\downarrow & & \cr & 1 & & & & & & \cr & & \ddots & & & & & \cr i\to & & & 1 & \ldots & f(\lambda) & & \cr & & & & \ddots & \vdots & & \cr & & & & & 1 & & \cr & & & & & & \ddots & \cr & & & & & & & 1\cr \end{pmatrix} $$ Die gleiche Operation für die Spalten ist äquivalent mit der Multiplikation von rechts mit der transponierten Matrix $$ \begin{pmatrix} & & & & & & & \cr & 1 & & & & & & \cr & & \ddots & & & & & \cr i\to & & & 1 & & & & \cr & & & \vdots & \ddots & & & \cr j\to & & & f(\lambda) & \ldots & 1 & & \cr & & & & & & \ddots & \cr & & & & & & & 1\cr \end{pmatrix} $$ Schließlich die Multiplikation der $i$-ten Zeile (Spalte) von $L(\lambda)$ mit einer Zahl $a\ne0$ ist äquivalent mit der Multiplikation von links (rechts) mit der Matrix $$ \begin{pmatrix} & & & & j\downarrow & & & \cr & 1 & & & & & & \cr & & \ddots & & & & & \cr & & & 1 & & & & \cr i\to & & & & a & & & \cr & & & & & 1 & & \cr & & & & & & \ddots & \cr & & & & & & & 1\cr \end{pmatrix} $$

3. Die $LU$-Zerlegung #

1. Die Produktzerlegung der Matrix $A\in\mathbb{C}^{n\times n}$ in $A=LU$, mit Subdiagonalmatrix $L$ und Superdiagonalmatrix $U$ heißt eine $LU$-Zerlegung. Wie das Beispiel der Matrix $A=({0\atop1}{1\atop0})$ zeigt, braucht es nicht unbedingt immer eine $LU$-Zerlegung zu geben.

Gibt es jedoch eine $LU$-Zerlegung, so ist diese Zerlegung unter gewissen Normierungsbedingungen auch eindeutig.

2. Satz: Die Zerlegung einer invertierbaren Matrix $A$ in das Produkt $A=L\cdot U$, einer Superdiagonalmatrix $U$ und einer normierten Subdiagonalmatrix $L$ (Subdiagonalmatrix mit lauter Einsen in der Diagonalen), ist eindeutig, wenn sie existiert.

Beweis: $A$ sei zerlegbar in die beiden Produkte $A=LU=\hat L\hat U$, mit Superdiagonalmatrizen $U$, $\hat U$ und normierten Subdiagonalmatrizen $L$, $\hat L$. Das Produkt von Superdiagonalmatrizen ist wieder eine Superdiagonalmatrix. Durch Transponierung erhält man entsprechend, daß das Produkt von Subdiagonalmatrizen wieder eine Subdiagonalmatrix ist und, daß weiter das Produkt von zwei normierten Subdiagonalmatrizen wieder eine normierte Subdiagonalmatrix ist. Nun folgt aus der Gleichung $U\hat U^{-1}=L^{-1}\hat L=I$, dann $\hat L^{-1}=L^{-1}$ und $\hat U^{-1}=U^{-1}$ (Eindeutigkeit der Inversen).     &#9744;

Die folgende Tabelle gibt an: die wesentliche Anzahl der durchzuführenden Operationen und der benötigte Speicherplatzbedarf bei vollbesetzten Matrizen, bei symmetrischen $(m,m)$-Bandmatrizen und dies in Abhängigkeit davon, ob eine Pivotsuche durchgeführt wird oder nicht.

3. Definition: $A\in\mathbb{C}^{n\times n}$ heißt ^{$(m,k)$-Bandmatrix} ($m,k\in{0,\ldots,n-1}$) mit linksseitiger Bandbreite $m$ und rechtsseitiger Bandbreite $k$, wenn $A$ insgesamt $m$ ^{Unterdiagonalen} und $k$ ^{Oberdiagonalen} enthält, welche Nichtnullelemente besitzen und sonst $A$ nur aus Nullen besteht. Eine ^{Tridiagonalmatrix} ist somit eine $(1,1)$-Bandmatrix, eine ^{obere Hessenbergmatrix} ist eine $(1,n-1)$-Bandmatrix und eine ^{untere Hessenbergmatrix} ist eine $(n-1,1)$-Bandmatrix. Hermitesche $(m,k)$-Bandmatrizen sind stets $(m,m)$-Bandmatrizen. Gelegentlich ist es nützlich zu sagen, daß $A$ eine $(m,k)$-Bandmatrix ist, wenn die entsprechenden Unter- oder Überdiagonalen Nichtnullelemente enthalten können. So wäre dann die Nullmatrix eine $(m,k)$-Bandmatrix für jedes $m,k\in{0,\ldots,n-1}$.

4. Lemma: $A\in\mathbb{C}^{n\times n}$ sei eine $(m,k)$-Bandmatrix und besitze eine $LU$-Zerlegung. Dann ist $L$ eine $(m,0)$-Bandmatrix und $U$ ist eine $(0,k)$-Bandmatrix.

Beweis: Es ist $a_{ij}=\sum_{\nu=1}^{\min(i,j)} \ell_{i\nu} u_{\nu j}$, also $$ \eqalignno{ &u_{1j} = a_{1j} \quad (j=1,\ldots,n), \qquad \ell_{i1} = {a_{i1} \over u_{11}}, \cr &i=2,\ldots,n: \cr &\qquad\qquad u_{ij} = a_{ij} - \sum_{\nu=1}^{i-1} \ell_{i\nu} u_{\nu j} \quad(j=i,i+1,\ldots,n),\cr &\qquad\qquad \ell_{ji} = {1\over u_{ii}} \left(a_{ji} - \sum_{\nu=1}^{i-1} \ell_{i\nu} u_{\nu j}\right) \quad(j=i+1,\ldots,n).\cr } $$     &#9744;

Umgekehrt ist das Produkt einer $(m,0)$-Bandmatrix mit einer $(0,k)$-Bandmatrix immer mindestens eine $(m,k)$-Bandmatrix.

5. Beispiel: Die Inversen von Bandmatrizen können einen hohen Auffüllungsgrad aufweisen. Dies zeigt $$ A = \pmatrix{ 1 & -1 & & & \cr -1 & 2 & -1 & & \cr & -1 & 2 & -1 & \cr & & -1 & 2 & -1\cr & & & -1 & 2\cr}, \qquad L^\top = U = \pmatrix{ 1 & -1 & & & \cr & 1 & -1 & & \cr & & 1 & -1 & \cr & & & 1 & -1\cr & & & & 1\cr}. $$ Die Matrix $A$ ist diagonal-dominant mit positiven Diagonaleinträgen, also positiv definit. $L$ und $U$ sind Dreiecksbandmatrizen, dennoch ist die Inverse von $A$ vollbesetzt, nämlich $$ A^{-1} = \pmatrix{ 5 & 4 & 3 & 2 & 1\cr 4 & 4 & 3 & 2 & 1\cr 3 & 3 & 3 & 2 & 1\cr 2 & 2 & 2 & 2 & 1\cr 1 & 1 & 1 & 1 & 1\cr}. $$

4. Die Gauß-Elimination #

1. Definition: Es bezeichnet $B[k]:=(b_{ij})_{i,j=1,\ldots,k}$ die Matrix gebildet aus den führenden $k$ Zeilen und Spalten. Entsprechend entgegengesetzt bezeichnet $B\!\left]\ell\right[:=(b_{ij})_{i,j=n-\ell,\ldots,n}$ die Matrix zusammengesetzt aus den letzten $\ell$ Zeilen und Spalten.

2. Das Gaußsche Eliminationsverfahren mit totaler Pivotsuche (GECP: Gaussian Elimination with complete pivoting) geht wie folgt vor sich. Es sei $A$ eine komplexe $n\times n$ Matrix.

1. Durchsuche die Matrix $A$ nach dem betragsmässig größten Element, das sogenannte Pivotelement. Dieses sei $a_{ij}$. Vertausche die $i$-te Zeile mit der ersten Zeile und vertausche die $j$-te Spalte mit der ersten Spalte. Die so neu gebildete Matrix heiße $A^{(0)}$.
2. Pivotiere bzgl. $a_{11}^{(0)}$, d.h. addiere ein Vielfaches der ersten Zeile zu allen darunterliegenden Zeilen, sodaß in der ersten Spalte, nur noch Nullen stehen (außer in dem ersten Element, dem Pivot, stehen dann nur noch Nullen im ersten Spaltenvektor). Die jetzt neu erhaltene Matrix heiße $A^{(1)}$.
3. In gleicher Weise führe den zweiten Schritt durch für $k=2,\ldots,n-1$, d.h. finde in $A^{(k-1)}\!\left]n-k+1\right[$ das betragsmässig größte Element, nenne es Pivotelement, und vertausche dann die entsprechende Zeile und Spalte, sodaß dieses Element in die Position $(k,k)$ kommt. Wenn das Pivotelement nicht Null ist, so pivotiere, andernfalls breche das Verfahren ganz ab: die Matrix $A$ ist nicht invertierbar. Die neu erhaltene Matrix nenne $A^{(k)}$.

Die Gauß-Elimination liefert entweder eine obere Dreiecksmatrix $A^{(n-1}$, oder die Information, daß $A$ nicht invertierbar ist. Wie erwähnt, heißen die Diagonalelemente $a_{kk}^{(k-1)}$ Pivotelemente. Eine Matrix $A$ ist genau dann invertierbar, wenn alle Pivotelemente nicht verschwinden (Produkt der Diagonalelemente ist bis auf das Vorzeichen der Wert einer Determinante). Wenn man das betragsmässig größte Element der Restmatrix $A\!\left]n-k\right[$ nur in der ersten Spalte von $A\!\left]n-k\right[$ sucht, also das betragsmässig größte Element der jeweils ersten Spalte sucht, anstatt in der ganzen Matrix $A\!\left]n-k\right[$ danach zu suchen (w.o.), so spricht man von Gauß-Elimination mit partieller Pivotsuche. Beschränkt man sich sogar bei der Suche auf ein beliebiges nicht verschwindendes Element, so spricht man von gewöhnlicher Gauß-Elimination. Bei Handrechnung verwendet man häufig die gewöhnliche Gauß-Elimination und wählt als Pivot möglichst kleine ganze Zahlen, falls vorhanden, z.B. 1. Programmiert wird i.d.R. die Gauß-Elimination mit partieller Pivotwahl, während GECP eher selten angewendet wird.

Nach dem $k$-ten Eliminationsschritt sieht die umgeformte Matrix $A$ dann wie folgt aus $$ A^{(k-1)} = \left( \begin{array}{cccc|cccc} a_{11}^{(0)} & a_{12}^{(0)} & \ldots & a_{1,k}^{(0)} && a_{1,k+1}^{(0)} & \ldots & a_{1,n}^{(0)}\cr & a_{22}^{(1)} & \ldots & a_{2,k}^{(1)} && a_{2,k+1}^{(1)} & \ldots & a_{2,n}^{(1)}\cr & & \ddots & \vdots && \vdots & \ddots & \vdots\cr 0 & & & a_{k,k}^{(k-1)}{\mskip 5mu} && a_{k,k+1}^{(k-1)} & \ldots & a_{k,n}^{(k-1)}\cr \hline &&& * && * & \ldots & *\cr &&& * && * & \ldots & *\cr &0&& \vdots && \vdots & \ddots & \vdots\cr &&& * && * & \ldots & *\cr \end{array} \right) $$ Wüßte man im voraus, welche Zeilen und Spalten jeweils zu vertauschen wären, so könnte man diese gleich im voraus durchführen. Mit dieser so vorpräparierten Matrix bräuchte man dann keinerlei Zeilen- und Spaltenvertauschungen durchzuführen. Alle diejenigen Matrizen, bei denen diese Vorvertauschungen schon durchgeführt sind, sollen CP heißen (^{completely pivoted}).

Es sei $$ g(A) = {\displaystyle{\max_{i,j,k} |a_{ij}^{(k)}|} \over \displaystyle{\max_{i,j} \left|a_{ij}\right|} } . $$ Diese Größe heißt Wachstum der Pivotstrategie.

Bei partieller Pivotwahl kann gelten $g(A)=2^n$. Bei totaler Pivotsuche sagt die 1992 falsifizierte Wilkinsonsche Vermutung (local copy), nach James Hardy Wilkinson (1919--1986), $g(A)\le n$. Diese Abschätzung ist scharf, wie man anhand von sogenanten Hadamard-Matrizen, Jacques S. Hadamard, (1865--1963), zeigen kann. Für komplexe Matrizen sind diese Schranken zu erhöhen. Hierzu Jane Day und Brian Peterson in Day/Peterson (1988):

Wilkinson's conjecture is very intriguing---easy to state, soon believed, and apparently very difficult to resolve.

3. Proposition: $A\in\mathbb{C}^{n\times n}$ sei invertierbar und CP. Es werde der $k$-te Eliminationsschritt der GECP durchgeführt, $k<i,j\le n$. Dann gilt $$ a_{ij}^{(k)} = {A_{1,\ldots,k,i}^{1,\ldots,k,j} \over A_{1\ldots k}^{1\ldots k}}. $$

Beweis: Da $A$ CP und da Pivotierung (nämlich Linearkombination von Zeilen) die Determinante und die Hauptminoren von $A$ nicht ändern, ergibt der Laplacesche Entwicklungssatz nach der $k$-ten Spalte (oder Zeile) sofort $\displaystyle{ A_{1,\ldots,k,j}^{1,\ldots,k,i} = a_{ij}^{(k)} A_{1\ldots k}^{1\ldots k}. }$ Aufgrund der Invertierbarkeit von $A$ folgt die gemachte Aussage.     &#9744;

4. Corollar: Das Pivotelement $a_{kk}^{(k-1)}$ bei GECP ist $a_{kk}^{(k-1)}={1\over x}$; $x$ ist das $(k,k)$-Element der Inverse der Matrix $A[1\ldots k]$, also von $(A[1\ldots k])^{-1}$. Insbesondere ist $a_{nn}^{(n-1)}$ das Reziproke eines Elementes von $A^{-1}$.

Beweis: Sei $B:=A[1\ldots k]$. Nach der Proposition ist dann $$ (B^{-1}){k,k} = {B{1,\ldots,k-1}^{1,\ldots,k-1} \over B_{1\ldots k}^{1\ldots k}} = {A_{1,\ldots,k-1}^{1,\ldots,k-1} \over A_{1\ldots k}^{1\ldots k}} = {1 \over a_{kk}^{(k-1)}} . $$     &#9744;

5. Bemerkungen: (1) GECP kann wie folgt interpretiert werden: Hat man die ersten $k-1$ Zeilen und Spalten gewählt, so wählt man die $k$-te Zeile und Spalte deswegen aus, weil dann die führende $k\times k$ Determinante maximalen Betrag aufweist.

(2) Alternativ kann man argumentieren: Hat man die ersten $k-1$ Zeilen und Spalten gewählt, so wählt man die $k$-te Zeile und Spalte deswegen aus, weil dann $\det A\!\left]n-k\right[$ minimalen Betrag aufweist.

(3) Geometrisch gesprochen: man wählt die $k$-te Zeile und Spalte deswegen aus, weil dann das $k$-Volumen des Spates gebildet aus den ersten $k$ Zeilenvektoren und Spaltenvektoren aus $A[k]$ maximal wird.

(4) Umgekehrt: man wählt die $k$-te Zeile und Spalte deswegen aus, weil damit das $(n-k)$-Volumen des projizierten Spates minimal wird, denn Pivotieren bzgl. $a_{kk}^{(k-1)}$ heißt Projizieren des Spates aufgespannt durch die Zeilen von $A^{(k-1)}\!\left]n-k+1\right[$ in den Spat aufgespannt durch die Zeilen von $A^{(k)}\!\left]n-k\right[$.

6. Satz: Darstellungssatz für die Gauß-Elimination nach Day/Peterson (1988). Es sei $A\in\mathbb{C}^{n\times n}$ invertierbar und $k<n$. Die Restmatrix nach dem $k$-ten Eliminationsschritt ist gegeben durch $$ A^{(k)}!\left]n-k\right[ = \left(A^{-1}!\left]n-k\right[\right)^{-1}. $$

D.h. die nach dem $k$-ten Eliminationsschritt noch nicht in Dreiecksform vorliegende Restmatrix $A^{(k)}\!\left]n-k\right[$ ist nichts anderes als von der eigentlichen Inversen $A^{-1}$ die invertierte Restmatrix $\bigl(A^{-1}\!\left]n-k\right[\bigr)^{-1}$. Damit stellt die Restmatrix nicht nur irgendeine Hilfsmatrix dar, sondern steht im Gegenteil mit der Inversen schon in engster Verbindung.

Beweis: Für $i,j>k$ sei $a_{ij}^{(k)}$ ein beliebiges Element von $A^{(k)}\!\left]n-k\right[$. Nach der vorhergehenden Proposition und dem Satz über Minoren Inverser gilt $$ a_{ij}^{(k)} = {A_{1,\ldots,k,i}^{1,\ldots,k,j} \over A_{1\ldots k}^{1\ldots k}} = {(-1)^{i+j} \left|A\right| (A^{-1}){k+1,\ldots,\hat\imath,\ldots,n} ^{k+1,\ldots,\hat\jmath,\ldots,n}\over A{1\ldots k}^{1\ldots k}} $$ Erneut wegen dem Satz über die Minoren Inverser gilt $$ \left|A^{-1}!\left]n-k\right[\right| = (A^{-1}){k+1,\ldots,n}^{k+1,\ldots,n} = {\alpha{k+1,\ldots,n}^{k+1,\ldots,n} \over \left|A\right|} = {A_{1\ldots k}^{1\ldots k} \over \left|A\right|} $$ Durch Einsetzen $$ a_{ij}^{(k)} = (-1)^{i+j} {(A^{-1}){k+1,\ldots,\hat\imath,\ldots,n}^{k+1,\ldots,\hat\jmath,\ldots,n} \over |A^{-1}!\left]n-k\right[|} . $$ Damit ist $a{ij}^{(k)}$ genau der $(i,j)$-Eintrag von $(A^{-1}\!\left]n-k\right[)^{-1}$.     &#9744;

7. Corollar: Führt man vor der eigentlichen Elimination sämtliche Zeilen- und Spaltenvertauschungen im voraus durch (also Matrix ist CP), so hat dies die Bedeutung, daß für $k=1,\ldots,n-1$ der Betrag der Determinante von $A^{-1}\!\left]n-k\right[$ nicht vergrößert werden kann durch irgendwelche Zeilen- und Spaltenvertauschungen der letzten $n-k+1$ Zeilen und Spalten von $A$.

Beweis: Nach den obigen Bemerkungen (3) und (4) ist GECP gleichbedeutend mit der Minimierung von $\left|\det A^{(k)}\!\left]n-k\right[\right|$ in jedem Schritt $k$, also der Maximierung von $\left|\det A^{-1}\!\left]n-k\right[\right|$.     &#9744;

Für positiv definite Matrizen $A$ ist GE immer anwendbar und zugleich liefert GE ein einfaches Kriterium zur Überprüfung auf positve Definitheit. Es gilt nämlich

8. Satz: Voraussetzung: $A\in\mathbb{C}^{n\times n}$ sei hermitesch. (Hermite, Charles (1822--1901))

Behauptung: GE durchführbar $\land$ $a_{ii}^{(k)}>0$ $\iff$ $A\succ0$.

Beweis: $A$ positiv definit $\iff$ $A_{1\ldots r}^{1\ldots r}>0$ $\iff$ $a_{ii}^{(k)}>0$, da Determinanten oder Hauptminoren sich nicht ändern bei Addition von Vielfachen von Zeilen zueinander.     &#9744;

9. Corollar: Jede positiv definite (hermitesche) Matrix $A$ besitzt genau eine $LU$-Zerlegung der Form $A=LU=LDL^\top$, mit einer normierten Subdiagonalmatrix $L$ und einer Diagonalmatrix $D$ mit lauter positiven Diagonalelementen.

Die Gauß-Elimination mit Diagonalstrategie mit positiven (Diagonal-) Pivots ist genau dann ausführbar, wenn die Matrix positiv definit ist. Also bei positiv definiten Matrizen sind Zeilen- und/oder Spaltenvertauschungen prinzipiell nicht erforderlich. Dies ist insofern von besonderem Interesse, als daß bei sehr großdimensionalen Matrizen ($n>1000$ beispielsweise) man besonders Wert legt auf einen geringen Auffüllungsgrad, welcher mit einer Pivotstrategie i.d.R. in einem Zielkonflikt steht. Konzentriert man sich daher bei positiv definiten Matrizen allein darauf, den Auffüllungsgrad gering zu halten, so bleibt dennoch die Gauß-Elimination immer durchführbar.

Genauso zeigt man: GE durchführbar $\iff$ $A_{1\ldots r}^{1\ldots r}\ne0$, da auch hier wieder sich die Hauptminoren nicht ändern bei Linearkombination von Zeilen. Damit hat man: Eine Matrix $A$ besitzt genau dann eine $LU$-Zerlegung, wenn alle führenden Hauptminoren nicht verschwinden. Dies deswegen, weil die Existenz einer $LU$-Zerlegung äquivalent ist mit der Durchführbarkeit der Gauß-Elimination ohne irgendwelche Zeilen- oder Spaltenvertauschungen.