Kurs:Numerik I/Zerlegung PA = LR

Einführung - Die Zerlegung PA = LR

Häufig ist in der numerischen Mathematik das Gleichungssystem $Ax=b$ für eingegebene (feste) reguläre Matrix $A\in \mathbb {R} ^{n\times n}$ und unterschiedliche rechte Seiten zu lösen.

Effizienzbetrachtung für den Algorithmus

In einem solchen Fall ist es ineffizient, für jede neue rechte Seiten wieder den gesamten Gauß-Algorithmus durchzuführen, da er bei jedem Durchlauf in Bezug auf $A$ wiederholt die gleichen Teiloperationen durchführen würde. Deshalb möchte man die beim Gauß-Algorithmus durchgeführten Teiloperationen von $A$ in irgendeiner Form speichern, die den Rechenaufwand reduziert.

Zerlegung in Teiloperationen

Dieses kann in Form einer Zerlegung von $A$ der Form $PA=LR$ geschehen, wie sie im folgenden Unterabschnitt eingeführt wird, wobei $L\in \mathbb {R} ^{n\times n}$ eine untere Dreiecksmatrix, $R\in \mathbb {R} ^{n\times n}$ eine obere Dreiecksmatrix und $P\in \mathbb {R} ^{n\times n}$ eine Permutationsmatrix ist.

Umformungen als Matrixmultiplikationen

Eine solche Zerlegung bzw. Faktorisierung von $A$ kann man, wie wir zeigen wollen, mittels des Gauß-Algorithmus mit Spaltenpivotsuche gewinnen.

Bemerkung zur PA-LR-Darstellung

Ein Permutationsmatrix ist invertierbar und $P^{-1}$ existiert. Damit lässt sich dann über die Gleichung $P\cdot A=L\cdot R$ Matrix $A$ über folgendes Matrixprodukt darstellen:

$A=P^{-1}\cdot L\cdot R$

Faktorisierung - Lösungsschritte

Liegt eine solche Faktorisierung vor, so kann man das Gleichungssystem $Ax=b$ bzw. $LRx=Pb$ für ein gegebenes $b\in \mathbb {R} ^{n}$ lösen, indem man hintereinander die beiden Dreieckssysteme

$Ly=Pb,\quad Rx=y$

löst, wobei die Lösung $y$ des ersten Systems die rechte Seite des zweiten Systems liefert.

Anwendung für unterschiedliche Vektoren b

Hat man das System $Ax=b$ für mehrere unterschiedliche rechte Seiten $b$ zu lösen, so muss man dann nur einmal die numerisch teuere Zerlegung $PA=LR$ bestimmen, während die Berechnung der Lösung $x$ gemäß $Ly=Pb,\quad Rx=y$ numerisch relativ durch die Verwendung von Dreiecksmatrizen nicht mehr so rechenintensiv ist.

Permutationsmatrizen

In der obigen Beschreibung haben wir die Notwendigkeit identifiziert, dass man für das Lösungsverfahren Zeilen vertauschen müssen. Eine einzelne Zeilvertauschung kann durch ein Matrix $Z_{(i,j)}$ realisiert werden. Durch Matrixmultiplikation kann diese Elementarmatrizen auch zu einer Pemutationsmatrix "zusammenfassen".

Beispiel - Vertauschung 1 und 4 Zeile

$Z_{(1,4)}={\begin{pmatrix}0&0&0&1&0\\0&1&0&0&0\\0&0&1&0&0\\1&0&0&0&0\\0&0&0&0&1\end{pmatrix}}\in \mathbb {R} ^{5\times 5}$.

Maxima - Vertauschung 1 und 4 Zeile

Geben Sie die $Z_{(1,4)}$ als $Z$ und die folgende Matrix in Maxima CAS ein

$A={\begin{pmatrix}1&2&3&4&5\\0&6&0&0&0\\0&0&7&0&0\\0&0&0&8&0\\0&0&0&0&9\end{pmatrix}}\in \mathbb {R} ^{5\times 5}$.

und berechnen Sie das Matrizenprodukt über Z.A mit "." als Operatorsymbol für die Matrixmultiplikation.

Definition - Permutationsmatrix

Jede Matrix $P\in \mathbb {R} ^{n\times n}$ mit genau einer Eins und sonst nur Nullen in jeder Zeile und Spalte heißt Permutationsmatrix.

Beispiel - Permutationsmatrix

Die folgende Matrix stellt eine Permutationsmatrix dar:

$P:={\begin{pmatrix}0&1&0&0\\0&0&1&0\\1&0&0&0\\0&0&0&1\end{pmatrix}}\in \mathbb {R} ^{4\times 4}.$

Definition - Permution

Jede bijektive Abbildung $\pi :\{1,...,n\}\to \{1,...,n\}$ heißt eine Permutation.

Zusammenhang Permutation - Permutationsmatrix

Offensichtlich ist $P\in \mathbb {R} ^{n\times n}$ genau dann eine Permutionsmatrix, wenn es eine Permutation $\pi :\{1,...,n\}\to \{1,...,n\}$ gibt, so dass

$P={\begin{pmatrix}e^{\pi (1)}&...&e^{\pi (n)}\end{pmatrix}}$

gilt, wobei $e^{j}\in \mathbb {R} ^{n}$ die $j$-te Spalte der Einheitsmatrix bezeichnet.

Beispiel - Permutation - Permutationmatrix

Die im obigen Beispiel genannte Permutationmatrix ist die durch folgende Permutation definiert

$\pi (1):=3,\quad \pi (2):=1,\quad \pi (3):=2,\quad \pi (4):=4,$

wobei $\pi :\{1,...,n\}\to \{1,...,n\}$ die folgende Permutationmatrix liefert

$P={\begin{pmatrix}e^{3}&e^{1}&e^{2}&e^{4}\end{pmatrix}}={\begin{pmatrix}e^{\pi (1)}&e^{\pi (2)}&e^{\pi (3)}&e^{\pi (4)}\end{pmatrix}}.$

Mittels Zusammenhang von Permutation und Permutationsmatrix kann man nun die Aussagen des folgenden Satzes erschließen.

Aufgabe - Studierende

Geben Sie in Maxima die Permutationsmatrix an, die die obige Permutation durchführt.

Satz - Permutationsmatrizen

Sei $P\in \mathbb {R} ^{n\times n}$ eine Permutationsmatrix und $\pi :\{1,...,n\}\to \{1,...,n\}$ bezeichne die zu der Matrix $P$ gehörende Permutation. Dann gilt:

• (i) $P$ ist eine orthogonale Matrix, d. h. es ist $P^{-1}=P^{T}$.
• (ii) Es gilt

$P^{-1}=P^{T}={\begin{pmatrix}e^{\pi ^{-1}(1)}&...&e^{\pi ^{-1}(n)}\end{pmatrix}}$.

• (iii) Für $\alpha :=(\alpha _{i})\in \mathbb {R} ^{n}$ gilt

$P\alpha ={\begin{pmatrix}\alpha _{\pi ^{-1}(1)}\\\vdots \\\alpha _{\pi ^{-1}(n)}\end{pmatrix}},\quad \alpha ^{T}P={\begin{pmatrix}\alpha _{\pi (1)}&...&\alpha _{\pi (n)}\end{pmatrix}}.$

• (iv) Für jede Matrix $A\in \mathbb {R} ^{n\times n}$ mit Zeilen $(z^{i})^{T}$ und Spalten $a^{i}$ $(i=1,...,n)$ gilt

(3.6) $PA={\begin{pmatrix}\left(z^{\pi ^{-1}(1)}\right)^{T}\\\vdots \\\left(z^{\pi ^{-1}(n)}\right)^{T}\end{pmatrix}},\quad AP={\begin{pmatrix}a^{\pi (1)}&...&a^{\pi (n)}\end{pmatrix}}.$

Beweis.

(i) folgt mit (3.5) aus

$\left(P^{T}P\right)_{ij}=\left[{\begin{pmatrix}\left(e^{\pi (1)}\right)^{T}\\\vdots \\\left(e^{\pi (n)}\right)^{T}\end{pmatrix}}{\begin{pmatrix}e^{\pi (1)}&...&e^{\pi (n)}\end{pmatrix}}\right]_{ij}=\left[\left(e^{\pi (i)}\right)^{T}e^{\pi (j)}\right]_{ij}=(\delta _{\pi (i),\pi (j)})_{ij}=(\delta _{ij})_{ij}=I.$

Die Behauptung in (ii) ergibt sich aus

$\left[{\begin{pmatrix}e^{\pi ^{-1}(1)}&...&e^{\pi ^{-1}(n)}\end{pmatrix}}{\begin{pmatrix}e^{\pi (1)}&...&e^{\pi (n)}\end{pmatrix}}\right]_{ij}=\sum _{\ell =1}^{n}\left(e^{\pi ^{-1}(\ell )}\right)_{i}\left(e^{\pi (j)}\right)_{\ell }=\sum _{\ell =1}^{n}\delta _{i,\pi ^{-1}(\ell )}\delta _{\pi (j),\ell }$

$=\sum _{k=1}^{n}\delta _{i,k}\delta _{\pi (j),\pi (k)}=\delta _{\pi (j),\pi (i)}=\delta _{ji}=\delta _{ij}.$

Die erste Behauptung in (iii) folgt aus

$P\alpha =\sum _{j=1}^{n}\alpha _{j}e^{\pi (j)}=\sum _{\ell =1}^{n}\alpha _{\pi ^{-1}(\ell )}e^{\ell }=(\alpha _{\pi ^{-1}(1)},...,\alpha _{\pi ^{-1}(n)})^{T},$

die zweite folgt mit (ii) und letzterer Identität aus

$\alpha ^{T}P=\left(P^{T}\alpha \right)^{T}=\left(\sum _{j=1}^{n}\alpha _{j}e^{\pi ^{-1}(j)}\right)^{T}=(\alpha _{\pi (1)},...,\alpha _{\pi (n)}).$

Die entsprechenden Matrixversionen in (iv) folgen analog aus

$PA=\sum _{j=1}^{n}e^{\pi (j)}(z^{j})^{T}=\sum _{\ell =1}^{n}e^{\ell }(z^{\pi ^{-1}(\ell )}={\begin{pmatrix}\left(z^{\pi ^{-1}(1)}\right)^{T}\\\vdots \\\left(z^{\pi ^{-1}(n)}\right)^{T}\end{pmatrix}}$

sowie unter Verwendung dieser Beziehung und (ii) aus

$AP=\left(P^{T}A^{T}\right)^{T}=\left[P^{T}{\begin{pmatrix}(a^{1})^{T}\\\vdots \\(a^{n})^{T}\end{pmatrix}}\right]^{T}={\begin{pmatrix}(a^{\pi (1)})^{T}\\\vdots \\(a^{\pi (n)})^{T}\end{pmatrix}}^{T}.$

q.e.d.

Man beachte, dass die Indizes $\pi ^{-1}(1),...,\pi ^{-1}(n)$ gemäß Aussage (ii) von Satz 3.10 gerade die Indizes der Einheitsvektoren, welche die Spalten von $P^{T}$ bzw. die Zeilen von $P$ bilden, sind. Somit bewirkt also die Multiplikation einer Matrix $A$ mit einer Permutationsmatrix von links bzw. rechts eine Permutation der Zeilen bzw. Spalten von $A$, die der Permutation der Zeilen bzw. Spalten von $P$ im Vergleich mit der Einheitsmatrix entspricht.

Beispiel 3.11

Seien

$P:={\begin{pmatrix}0&0&1\\1&0&0\\0&1&0\end{pmatrix}},\quad A:={\begin{pmatrix}1&2&3\\4&5&6\\7&8&9\end{pmatrix}}.$

Dann folgt

$PA={\begin{pmatrix}0&0&1\\1&0&0\\0&1&0\end{pmatrix}}{\begin{pmatrix}1&2&3\\4&5&6\\7&8&9\end{pmatrix}}={\begin{pmatrix}7&8&9\\1&2&3\\4&5&6\end{pmatrix}},$

$AP={\begin{pmatrix}1&2&3\\4&5&6\\7&8&9\end{pmatrix}}{\begin{pmatrix}0&0&1\\1&0&0\\0&1&0\end{pmatrix}}={\begin{pmatrix}2&3&1\\5&6&4\\8&9&7\end{pmatrix}}.$

---

In numerischen Implementierungen erfolgt die Abspeicherung einer Permutationsmatrix mit der zugehörigen Permutation $\pi$ in Form eines Vektors

$(\pi ^{-1}(1),...,\pi ^{-1}(n))^{T}\in \mathbb {R} ^{n}$ oder $(\pi (1),...,\pi (n))^{T}\in \mathbb {R} ^{n}$.

Eine besondere Rolle spielen Elementarpermutationen $\pi :\{1,...,n\}\to \{1,...,n\}$, die zwei Zahlen vertauschen und die restlichen Zahlen unverändert lassen. Im Fall einer Elementarpermutation gibt es also zwei Zahlen $i,r\in \{1,...,n\}$ mit

(3.7) $\pi (i)=r,\quad \pi (r)=i,\quad \pi (j)=j\quad (j\notin \{i,r\}).$

Hier gilt wegen

$\pi ^{-1}(r)=i=\pi (r),\quad \pi ^{-1}(i)=r=\pi (i),\quad \pi ^{-1}(j)=j=\pi (j)\quad (j\notin \{i,r\})$

die Identität $\pi ^{-1}=\pi$, so dass sich für die zu $\pi$ gehörende Permutationsmatrix $P\in \mathbb {R} ^{n\times n}$ (vgl. (3.5)) ergibt:

(3.8) $P^{-1}=P^{T}=P.$

3.3.2 Frobenius-Matrizen

Wir betrachten eine weitere wichtige Klasse von Matrizen.

Definition - Frobeniusmatrix

Sei $k\in \{1,...,n\}$. Jede Matrix der Form

(3.9) ${\begin{pmatrix}1&&&&&\\&\ddots &&&&\\&&1&&&\\&&-l_{k+1,k}&1&&\\&&\vdots &&\ddots &\\&&-l_{n,k}&&&1\end{pmatrix}}\in \mathbb {R} ^{n\times n}$

heißt Frobenius-Matrix vom Index $k$.

Bemerkung - Frobeniusmatrix

Eine Frobenius-Matrix vom Index $k$ unterscheidet sich von der Einheitsmatrix gleicher Größe also nur in der $k$-ten Spalte und dort auch nur unterhalb der Diagonalen. Insbesondere lässt sich die prinzipielle Vorgehensweise bei den Zeilenoperationen der $k$-ten Stufe des Gauß-Algorithmus durch Multiplikation mit einer Frobenius-Matrix vom Index $k$ beschreiben. So gilt für Vektoren $z^{j}\in \mathbb {R} ^{n}$ $(j=1,...,n)$

(3.10) ${\begin{pmatrix}1&&&&&\\&\ddots &&&&\\&&1&&&\\&&-l_{k+1,k}&1&&\\&&\vdots &&\ddots &\\&&-l_{n,k}&&&1\end{pmatrix}}{\begin{pmatrix}\left(z^{1}\right)^{T}\\\vdots \\\left(z^{n}\right)^{T}\end{pmatrix}}={\begin{pmatrix}\left(z^{1}\right)^{T}\\\vdots \\\left(z^{k}\right)^{T}\\\left(z^{k+1}\right)^{T}-l_{k+1,k}\left(z^{k}\right)^{T}\\\vdots \\\left(z^{n}\right)^{T}-l_{n,k}\left(z^{k}\right)^{T}\end{pmatrix}}$

Offenbar lässt sich die Frobenius-Matrix in (3.9) mit

(3.11) $f^{k}:=(0,...,0,l_{k+1,k},...,l_{n,k})^{T}\in \mathbb {R} ^{n}$

wegen

$f^{k}(e^{k})^{T}={\begin{pmatrix}0\\\vdots \\0\\0\\l_{k+1,k}\\\vdots \\l_{n,k}\end{pmatrix}}(0,...,0,\overbrace {1} ^{k{\text{-te Stelle}}},0,...,0)={\begin{pmatrix}0&...&0&0&0&...&0\\\vdots &\ddots &\vdots &\vdots &\vdots &\ddots &\vdots \\0&...&0&0&0&...&0\\0&...&0&0&0&...&0\\0&...&0&l_{k+1,k}&0&...&0\\\vdots &\ddots &\vdots &\vdots &\vdots &\ddots &\vdots \\0&...&0&l_{n,k}&0&...&0\end{pmatrix}}$

in der Form

(3.12) $F_{k}=I-f^{k}(e^{k})^{T}$

darstellen, wobei $I\in \mathbb {R} ^{n\times n}$ die Einheitsmatrix und $e^{k}\in \mathbb {R} ^{n}$ wieder die $k$-te Spalte von $I$ bezeichnet.

Lemma 3.13

Für $k=1,...,n-1$ sind die Frobenius-Matrizen $F_{k}\in \mathbb {R} ^{n\times n}$ vom Index $k$ regulär und es gilt für $k=1,...,n-1$

(3.13) $F_{k}^{-1}={\begin{pmatrix}1&&&&&\\&\ddots &&&&\\&&1&&&\\&&l_{k+1,k}&1&&\\&&\vdots &&\ddots &\\&&l_{n,k}&&&1\end{pmatrix}}=I+f^{k}(e^{k})^{T}$

sowie

(3.14) $F_{1}^{-1}\cdots F_{n-1}^{-1}={\begin{pmatrix}1&&&&\\l_{21}&1&&&\\\vdots &l_{32}&1&&\\\vdots &\vdots &\ddots &\ddots &\\l_{n1}&l_{n2}&...&l_{n,n-1}&1\end{pmatrix}}=I+\sum _{j=1}^{n-1}f^{j}(e^{j})^{T}.$

Beweis.

Für $F_{k}$ hat man mit (3.11) die Darstellung (3.12). Wegen

$\left[I+f^{k}(e^{k})^{T}\right]\underbrace {\left[I-f^{k}(e^{k})^{T}\right]} _{=F_{k}}=I+f^{k}(e^{k})^{T}-f^{k}(e^{k})^{T}-f^{k}\underbrace {\left[(e^{k})^{T}f^{k}\right]} _{=0}(e^{k})^{T}=I$

folgt

$\det \left(I+f^{k}(e^{k})^{T}\right)\det \left(I+f^{k}(e^{k})^{T}\right)=1,$

also die Regularität von $F_{k}$ sowie die behauptete Darstellung von $F_{k}^{-1}$. Im Folgenden soll nun mittels vollständiger Induktion die Identität

(3.15) $F_{1}^{-1}\cdots F_{k}^{-1}=I+\sum _{j=1}^{k}f^{j}(e^{j})^{T},\quad k=1,...,n-1$

nachgewiesen werden, welche im Fall $k:=n-1$ der Formel (3.14) entspricht.

Die Darstellung in (3.15) ist sicher richtig für $k:=1$. Wir nehmen nun weiter an, dass sie für ein beliebiges $k\in \{1,...,n-2\}$ richtig ist. Dann gilt die Darstellung in (3.15) auch für $k+1$, denn

$F_{1}^{-1}\cdots F_{k}^{-1}F_{k+1}^{-1}=\left[I+\sum _{j=1}^{k}f^{j}(e^{j})^{T}\right]\left[I+f^{k+1}(e^{k+1})^{T}\right]$

$=I+f^{k+1}(e^{k+1})^{T}+\sum _{j=1}^{k}f^{j}(e^{j})^{T}+\sum _{j=1}^{k}f^{j}\underbrace {\left[(e^{j})^{T}f^{k+1}\right]} _{=0}(e^{k+1})^{T}=I+\sum _{j=1}^{k+1}f^{j}(e^{j})^{T}.$

q.e.d.

Lemma 3.14

Sei $F_{k}$ eine wie in (3.12) mit (3.11) dargestellte Frobenius-Matrix vom Index $k$ und $P$ eine Permutationsmatrix mit zugehöriger Elementarpermutation $\pi$ von der Form (3.7) mit $i,r\in \{k+1,...,n\}$. Dann entsteht die Matrix $PF_{k}P$ aus $F_{k}$ durch Vertauschen der Einträge $i$ und $r$ in der $k$-ten Spalte, d. h.

$PF_{k}P=I-\left(Pf^{k}\right)(e^{k})^{T}.$

Beweis.

Die Aussage ergibt sich unmittelbar aus

$PF_{k}P=P\left[I-f^{k}(e^{k})^{T}\right]P=\underbrace {P^{2}} _{=I}-[Pf^{k}][\underbrace {(e^{k})^{T}P} _{=(e^{k})^{T}}],$

wobei (3.8) und Satz 3.10 (iii) sowie die Forderungen für $i$ und $r$ eingehen.

q.e.d.

3.3.3 Die LR-Zerlegung mittels Gauß-Algorithmus

Im Folgenden wird die allgemeine Vorgehensweise beim Gauß-Algorithmus zur sukzessiven Erzeugung von Matrizen $A^{(k)}\in \mathbb {R} ^{n\times n}$ der Form

(3.16) $A^{(k)}={\begin{pmatrix}a_{11}^{(k)}&a_{12}^{(k)}&...&...&...&a_{1n}^{(k)}\\&a_{22}^{(k)}&...&...&...&a_{1n}^{(k)}\\&&\ddots &&&\vdots \\&&&a_{kk}^{(k)}&...&a_{kn}^{(k)}\\&&&\vdots &\ddots &\vdots \\&&&a_{nk}^{(k)}&...&a_{nn}^{(k)}\end{pmatrix}}\in \mathbb {R} ^{n\times n}$

mit Spaltenpivotsuche als Folge spezieller Matrix-Operationen beschrieben. Und zwar wird im $k$-ten Schritt entsprechend Algorithmus 1 eine Zeilenvertauschung $A^{(k)}\to P_{k}A^{(k)}$ vorgenommen. Hierbei ist $P_{k}\in \mathbb {R} ^{n\times n}$ eine elementare Permutationsmatrix, die nur eine Vertauschung der Zeilen $k$ und $p_{k}\geq k$ von $A^{(k)}$ bewirkt. ($p_{k}=k$ und somit $P_{k}=I$ ist möglich.) Damit ergibt sich gemäß (3.6) und (3.10)

(3.17) $A^{(k+1)}=F_{k}P_{k}A^{(k)}$

mit

(3.18) $F_{k}:={\begin{pmatrix}1&&&&&\\&\ddots &&&&\\&&1&&&\\&&-l_{k+1,k}&1&&\\&&\vdots &&\ddots &\\&&-l_{n,k}&&&1\end{pmatrix}}$

für

(3.19) $l_{p_{k},k}:={\frac {a_{kk}^{(k)}}{a_{p_{k},k}^{(k)}}},\quad l_{ik}={\frac {a_{ik}^{(k)}}{a_{p_{k},k}^{(k)}}}\quad (i=k+1,...,n,i\neq p_{k})$

sowie

(3.20) $P_{k}:={\begin{pmatrix}1&&&&&&&&&&\\&\ddots &&&&&&&&&\\&&1&&&&&&&&\\&&&0&&&&1&&&\\&&&&1&&&&&&\\&&&&&\ddots &&&&&\\&&&&&&1&&&&\\&&&1&&&&0&&&\\&&&&&&&&1&&\\&&&&&&&&&\ddots &\\&&&&&&&&&&1\end{pmatrix}}$ ${\begin{array}{ll}\leftarrow &{\text{Zeile }}k\\&\\&\\&\\\leftarrow &{\text{Zeile }}p_{k}\end{array}}$

Der Index $p_{k}$ bezeichnet dabei die Position der Zeile aus $A^{(k)}$, welche das Pivotelement enthält.

Satz 3.15

Mit den Definitionen (3.16) - (3.20) gilt die Identität $PA=LR$ für

$P:=P_{n-1}\cdots P_{1},\quad R:=A^{(n)}$

und

(3.21) $L:={\begin{pmatrix}1&&&&\\{\hat {l}}_{21}&1&&&\\\vdots &{\hat {l}}_{32}&1&&\\\vdots &\vdots &\ddots &\ddots &\\{\hat {l}}_{n1}&{\hat {l}}_{n2}&...&{\hat {l}}_{n,n-1}&1\end{pmatrix}}$

mit

(3.22) ${\begin{pmatrix}0\\\vdots \\0\\1\\{\hat {l}}_{k+1,k}\\\vdots \\{\hat {l}}_{n,k}\end{pmatrix}}:=P_{n-1}\cdots P_{k+1}{\begin{pmatrix}0\\\vdots \\0\\1\\l_{k+1,k}\\\vdots \\l_{n,k}\end{pmatrix}}$

Beweis.

Für $k=1,2,...$ gilt mit (3.17) sowie (3.8)

$A^{(2)}=F_{1}P_{1}A=F_{1}(P_{1}A),$

$A^{(3)}=F_{2}P_{2}A^{(2)}=F_{2}(P_{2}F_{1}\underbrace {P_{2})(P_{2}} _{=I}P_{1}A),$

$A^{(4)}=F_{3}P_{3}A^{(3)}=F_{3}(P_{3}F_{2}\underbrace {P_{3})(P_{3}} _{=I}P_{2}F_{1}P_{2}\underbrace {P_{3})(P_{3}} _{=I}P_{2}P_{1}A)$

und so weiter, was schließlich

(3.23) $R=A^{(n)}={\hat {F}}_{n-1}\cdots {\hat {F}}_{1}PA$

ergibt mit

$P:=P_{n-1}\cdots P_{1}$

und den Frobenius-Matrizen ${\hat {F}}_{n-1}=F_{n-1}$ und

${\hat {F}}_{k}:=P_{n-1}\cdots P_{k+1}F_{k}P_{k+1}\cdots P_{n-1}={\begin{pmatrix}1&&&&&\\&\ddots &&&&\\&&1&&&\\&&-{\hat {l}}_{k+1,k}&\ddots &&\\&&\vdots &&\ddots &\\&&-{\hat {l}}_{nk}&&&1\end{pmatrix}}$

für $k=1,...,n-2$, wobei in die letzte Identität Lemma 3.14 eingeht. Eine Umformung von (3.23) liefert dann

$PA={\begin{pmatrix}{\hat {F}}_{1}^{-1}&...&{\hat {F}}_{n-1}^{-1}\end{pmatrix}}R=LR,$

wobei die letzte Gleichheit mit Lemma 3.13 folgt. Damit ist alles bewiesen.

q.e.d.

Man beachte, dass die Matrix $L$ (3.21) also gerade aus den aktuellen Umrechnungsfaktoren ${\hat {l}}_{ik}$ $(i>k)$ gebildet wird, nachdem (sofern erforderlich) die Zeilenvertauschung, welche die Zeile mit dem jeweils gewählten Pivotelement an die richtige Position bringt, erfolgt ist. In Implementierungen werden die frei werdenden Anteile des linken Dreiecks der Matrix $A$ sukzessive mit den Einträgen der unteren Dreiecksmatrix $L$ überschrieben, während sich in dem rechten Dreieck der Matrix $A$ die Einträge der Dreiecksmatrix $R$ ergeben. Die Permutationsmatrix $P$, deren Zeilen $(e^{r_{i}})^{T},i=1,...,n$ genannt seien, lässt sich einfach in Form eines Buchhaltungsvektors $r\in \mathbb {N} ^{n}$ angeben, und es gilt, wie man mit Satz 3.10 erschließt,

${\begin{pmatrix}r_{1}\\\vdots \\r_{n}\end{pmatrix}}=P{\begin{pmatrix}1\\\vdots \\n\end{pmatrix}}$

Wir wollen die Vorgehensweise an einem Beispiel vorführen.

Beispiel 3.16

Gegeben sei die Matrix

$A:={\begin{pmatrix}0&0&1&1\\2&2&2&2\\1&2&2&2\\1&2&3&6\end{pmatrix}}.$

Nach Anhängen des für die Speicherung der Zeilenpermutationen zuständigen Buchhaltervektors liefert der Algorithmus folgendes (unterhalb der Treppe ergeben sich sukzessive die Einträge von $L$ aus (3.21), (3.22)):

${\begin{pmatrix}0&0&1&1\\2&2&2&2\\1&2&2&2\\1&2&3&6\end{pmatrix}},{\begin{pmatrix}1\\2\\3\\4\end{pmatrix}}{\underset {\longrightarrow }{\text{ Zeilentausch }}}{\begin{pmatrix}2&2&2&2\\0&0&1&1\\1&2&2&2\\1&2&3&6\end{pmatrix}},{\begin{pmatrix}2\\1\\3\\4\end{pmatrix}}$

${\underset {\longrightarrow }{\text{Elimination }}}{\begin{pmatrix}2&2&2&2\\0&0&1&1\\{\frac {1}{2}}&{\underline {1}}&1&1\\{\frac {1}{2}}&1&2&5\end{pmatrix}},{\begin{pmatrix}2\\1\\3\\4\end{pmatrix}}{\underset {\longrightarrow }{\text{ Zeilentausch }}}{\begin{pmatrix}2&2&2&2\\{\frac {1}{2}}&1&1&1\\0&0&1&1\\{\frac {1}{2}}&1&2&5\end{pmatrix}},{\begin{pmatrix}2\\3\\1\\4\end{pmatrix}}$

${\underset {\longrightarrow }{\text{Elimination }}}{\begin{pmatrix}2&2&2&2\\{\frac {1}{2}}&1&1&1\\0&0&{\underline {1}}&1\\{\frac {1}{2}}&1&1&4\end{pmatrix}},{\begin{pmatrix}2\\3\\1\\4\end{pmatrix}}{\underset {\longrightarrow }{\text{ Elimination }}}{\begin{pmatrix}2&2&2&2\\{\frac {1}{2}}&1&1&1\\0&0&1&1\\{\frac {1}{2}}&1&1&3\end{pmatrix}},{\begin{pmatrix}2\\3\\1\\4\end{pmatrix}}$

Dabei ist das jeweils gewählte Pivotelement unterstrichen. Der letzte Permutationsvektor $(2,3,1,4)^{T}$ besagt, dass

$\pi (1)=3,\quad \pi (2)=1,\quad \pi (3)=2,\quad \pi (4)=4$

für die zu $P$ gehörende Permutation $\pi$ gilt und dass also $PA$ aus $A$ hervorgeht, indem man die erste Zeile von $A$ in die dritte Position bringt, die zweite in die erste, usw. Es ergibt sich somit die Faktorisierung

$PA={\begin{pmatrix}2&2&2&2\\1&2&2&2\\0&0&1&1\\1&2&3&6\end{pmatrix}}={\begin{pmatrix}1&&&\\{\frac {1}{2}}&1&&\\0&0&1&\\{\frac {1}{2}}&1&1&1\end{pmatrix}}{\begin{pmatrix}2&2&2&2\\&1&1&1\\&&1&1\\&&&3\end{pmatrix}}=LR.$

---

Man beachte, dass die Zerlegung $PA=LR$ einer Matrix nicht eindeutig ist. Man könnte ja beispielsweise die Matrix $L$ mit einem Skalar $\alpha \neq 0$ und die Matrix $R$ mit dem Skalar $1/\alpha$ multiplizieren.

Mit dem Gauß-Algorithmus kann man bekanntlich auch die Determinante von $A$ berechnen. So gilt im Fall, dass eine Zerlegung $PA=LR$ vorliegt,

$\det(P)=(-1)^{\sigma }$

und

$\det(PA)=\det(P)\det(A)=\det(L)\det(R)=\prod _{i=1}^{n}r_{ii},$

wobei $\sigma$ die Anzahl von paarweisen Zeilenvertauschungen ist, die die Überführung von $I$ in $P$ erfordert bzw. welche beim Gauß-Algorithmus vorgenommen wird und die $r_{ii}$ die Diagonalelemente von $R$ sind. Demnach hat man

$\det(A)=(-1)^{\sigma }\prod _{i=1}^{n}r_{ii}.$

3.3.4 Nachiteration

Aufgrund von Rundungsfehlern errechnet man in der Praxis nicht eine Zerlegung $LR$, sondern eine Zerlegung ${\tilde {L}}{\tilde {R}}$ von $PA$, so dass

$PA\approx {\tilde {L}}{\tilde {R}}$

gilt. Statt der (hier als eindeutig vorausgesetzten Lösung) $x$ von $Ax=b$ bzw. $PAx=Pb$ berechnet man demnach unter Verwendung von ${\tilde {L}}$ und ${\tilde {R}}$ eine Näherungslösung $x^{0}$ und den durch sie erzeugten Defekt

$d^{0}:=Pb-PAx^{0}\neq 0.$

Daher ist die folgende Nachiteration sinnvoll: wiederum unter Verwendung der vorliegenden Zerlegung ${\tilde {L}}{\tilde {R}}$ von $PA$ bestimmt man die Lösung $\delta x^{0}$ der Defektgleichung

(3.24) $PA\delta x=d^{0}$

und setzt man anschließend

$x^{1}:=x^{0}+\delta x^{0}.$

Bei exakter Lösung des Systems (3.24) (mit $PA$ und nicht ${\tilde {L}}{\tilde {R}}$) hätte man dann

$PAx^{1}=PAx^{0}+PA\delta x^{0}=Pb-d^{0}+d^{0}=Pb.$

Da man i. a. jedoch nicht exakt rechnet, könnte man diesen Prozess wiederholen. Normalerweise genügt es jedoch, $d^{0}$ und $\delta x^{0}$ mit doppelter Genauigkeit zu berechnen und die beschriebene Nachiteration nur einmal durchzuführen (vgl. Deuflhard/Hohmann).

Beispiel 3.17

Wir betrachten das Gleichungssystem $Ax=b$ mit

$A:={\begin{pmatrix}1.05&1.02\\1.04&1.02\end{pmatrix}},\quad b:={\begin{pmatrix}1\\2\end{pmatrix}}.$

Die Lösung $x$ des Systems lautet

$x\approx (-100,103.921\ 569)^{T}.$

Gauß-Elimination ohne Zeilenvertauschung liefert bei 3-stelliger Rechnung

${\tilde {L}}:={\begin{pmatrix}1&0\\.990&1\end{pmatrix}},\quad {\tilde {R}}:={\begin{pmatrix}1.05&1.02\\0&0.01\end{pmatrix}}$

sowie

${\tilde {L}}{\tilde {R}}-A={\begin{pmatrix}0&0\\5_{10}-4&2_{10}-4\end{pmatrix}}.$

Man errechnet

$x^{0}:=(-97.1,101)^{T}$

mit dem Defekt

$d^{0}:=b-Ax^{0}={\begin{cases}(0,0)^{T}&(3{\text{-stellig}}),\\(0.065,0.035)^{T}&(6{\text{-stellig}}).\end{cases}}$

Nachiteration mit 6-stelliger Rechnung ergibt

$x^{1}:=x^{0}+\delta x^{0}=(-99.9,104)^{T}.$

3.3.5 Direkte LR-Zerlegung

In gewissen Situationen ist es möglich und zwecks Ausnutzung vorhandener Strukturen der Matrix $A:=(a_{ij})\in \mathbb {R} ^{n\times n}$ auch sinnvoll, auf eine Pivotstrategie zu verzichten und mittels einer unteren Dreiecksmatrix $L:=(l_{ij})\in \mathbb {R} ^{n\times n}$ und einer oberen Dreiecksmatrix $R:=(r_{ij})\in \mathbb {R} ^{n\times n}$ eine LR-Zerlegung der Form

(3.25) ${\begin{pmatrix}a_{11}&a_{12}&...&a_{1n}\\a_{21}&a_{22}&...&a_{2n}\\\vdots &\vdots &\ddots &\vdots \\a_{n1}&a_{n2}&...&a_{nn}\end{pmatrix}}={\begin{pmatrix}1&&&\\l_{21}&1&&\\\vdots &\ddots &\ddots &\\l_{n1}&...&l_{n,n-1}&1\end{pmatrix}}{\begin{pmatrix}r_{11}&r_{12}&...&r_{1n}\\&r_{22}&...&r_{2n}\\&&\ddots &\vdots \\&&&r_{nn}\end{pmatrix}}$

auf direkte Weise zu bestimmen. Eine solche Zerlegung einer regulären Matrix $A$ existiert genau dann, wenn $\det(H_{k})\neq 0$ $(k=1,...,n)$ für die Hauptuntermatrizen

(3.26) $H_{k}:={\begin{pmatrix}a_{11}&...&a_{1k}\\\vdots &\ddots &\vdots \\a_{k1}&...&a_{kk}\end{pmatrix}}\in \mathbb {R} ^{k\times k},\quad k=1,...,n$

von $A$ gilt (z. B. Sätze 2.14 und 2.17 bei Kanzow). Wegen $\det(A)\neq 0$ ist dann auch $\det(R)\neq 0$, also $r_{ii}\neq 0$ $(i=1,...,n)$ und damit das folgende Vorgehen möglich.

Existiert eine LR-Zerlegung von $A$ wie in (3.25), so verwendet den Ansatz in (3.25), um für die $n^{2}$ gesuchten Größen $r_{ik}$ $(i\leq k)$ und $l_{ik}$ $(i>k)$ die $n^{2}$ Bestimmungsgleichungen

$a_{ik}=\sum _{j=1}^{n}l_{ij}r_{jk},\quad i,k=1,...,n$

zu erhalten, welche wegen $l_{ij}=0$ für $i<j$ und $r_{jk}=0$ für $j>k$ mit

(3.27) $a_{ik}=\sum _{j=1}^{\min(i,k)}l_{ij}r_{jk},\quad i,k=1,...,n$

identisch sind. Dabei gibt es verschiedene Möglichkeiten, wie man aus den Gleichungen (3.27) die Einträge von $L$ und $R$ bestimmen kann. Zum Beispiel führt eine Berechnung der Zeilen von $R$ und der Spalten von $L$ entsprechend der Parkettierung auf den folgenden Algorithmus:

Algorithmus 2 (Direkte LR-Zerlegung)

(0) Gib $A:=(a_{ij})\in \mathbb {R} ^{n\times n}$ mit $\det(A)\neq 0$ und setze $i:=1$.

(1) Berechne

$r_{ik}:=a_{ik}-\sum _{j=1}^{i-1}l_{ij}r_{jk}\quad (k=i,...,n),$

$l_{ki}:=\left(a_{ki}-\sum _{j=1}^{i-1}l_{kj}r_{ji}\right)/r_{ii}\quad (k=i+1,...,n).$

(2) Falls $i=n$, stop!

(3) Setze $i:=i+1$ und gehe nach (1).

Für eine solche direkte LR-Zerlegung sind insgesamt $(n^{3}/3)(1+0(1/n))$, d. h. die gleiche Größenordnung von Multiplikationen und Divisionen wie für den Gauß-Algorithmus erforderlich.

Cholesky-Zerlegung positiv definiter Matrizen

In diesem Abschnitt zeigen wir, dass die eben vorgestellte LR-Zerlegungen für positiv definite Matrizen besonders attraktiv ist.

Definition - positiv definite Matrix

Eine Matrix $A\in \mathbb {R} ^{n\times n}$ heißt positiv definit, falls $A$ symmetrisch, d. h. $A=A^{T}$ ist und falls gilt:

$\langle x,Ax\rangle =x^{T}Ax>0,\quad x\in \mathbb {R} ^{n}\setminus \{0\}.$

Bemerkung - positiv definite Matrix

Mit dieser Definition selbst lässt sich, wie es ähnlich häufig in der Mathematik der Fall ist, nur schwer arbeiten. Deshalb sind äquivalente Bedingungen von Interesse.

Hesse-Matrix und ein lokales Minimum einer Fehlerfunktion

Sei $f\colon U\subset \mathbb {R} ^{n}\to \mathbb {R}$ eine zweimal stetig differenzierbare Funktion. Dann ist die Hesse-Matrix von $f$ am Punkt $x=(x_{1},\ldots ,x_{n})\in U$ definiert durch

$\operatorname {H} _{f}(x):=\left({\frac {\partial ^{2}f}{\partial x_{i}\partial x_{j}}}(x)\right)_{i,j=1,\dots ,n}={\begin{pmatrix}{\frac {\partial ^{2}f}{\partial x_{1}\partial x_{1}}}(x)&{\frac {\partial ^{2}f}{\partial x_{1}\partial x_{2}}}(x)&\cdots &{\frac {\partial ^{2}f}{\partial x_{1}\partial x_{n}}}(x)\\[0.5em]{\frac {\partial ^{2}f}{\partial x_{2}\partial x_{1}}}(x)&{\frac {\partial ^{2}f}{\partial x_{2}\partial x_{2}}}(x)&\cdots &{\frac {\partial ^{2}f}{\partial x_{2}\partial x_{n}}}(x)\\\vdots &\vdots &\ddots &\vdots \\{\frac {\partial ^{2}f}{\partial x_{n}\partial x_{1}}}(x)&{\frac {\partial ^{2}f}{\partial x_{n}\partial x_{2}}}(x)&\cdots &{\frac {\partial ^{2}f}{\partial x_{n}\partial x_{n}}}(x)\end{pmatrix}}.$

Ist die Hesse-Matrix an einer Stelle $x=(x_{1},\ldots ,x_{n})\in U$ positiv definit, wobei gleichzeitig der Gradient der Nullvektor ist, dann besitzt $f$ ein lokales Minimum.

Lemma - Kriterien für positive Definitheit

Folgende Aussagen sind für eine Matrix $A:=(a_{ij})\in \mathbb {R} ^{n\times n}$ mit $A=A^{T}$ äquivalent:

• (i) $A$ ist positiv definit.
• (ii) Alle Eigenwerte $\lambda _{i}$ $(i=1,...,n)$ von $A$ sind reell und positiv.
• (iii) Die Determinanten der $n$ Hauptuntermatrizen $H_{k}$ $(k=1,...,n)$ von $A$ in (3.26) sind alle positiv.

Bemerkung - Kriterien für positive Definitheit

Den Beweis des Lemmas findet man in Büchern der Linearen Algebra. Im Folgenden beweisen wir Eigenschaften positiv definite Matrizen $A$ die positiv Definitheit von allen Untermatrizen und der Invertierbarkeit der Matrix.

Lemma - positiv Definitheit von Untermatrizen

Die Matrix $A:=(a_{ij})\in \mathbb {R} ^{n\times n}$ sei positiv definit. Dann gilt:

• (i) $A$ ist regulär,
• (ii) alle Untermatrizen von $A$ der Form

(3.29) ${\begin{pmatrix}a_{kk}&...&a_{k,k+\ell }\\\vdots &\ddots &\vdots \\a_{k+\ell ,k}&...&a_{k+\ell ,k+\ell }\end{pmatrix}}\in \mathbb {R} ^{(\ell +1)\times (\ell +1)}$

sind positiv definit,

• (iii) $\det(A)>0$.

Beweis

Der Beweis erfolgt in der Reihenfolge der drei Eigenschaften von positiv definiten Matrizen aus dem Lemma.

Beweis (i)

(i) Wäre $A$ singulär, so gäbe es einen Vektor $x\in \mathbb {R} ^{n}\setminus \{0\}$ mit $Ax=0$. Damit wäre auch $x^{T}Ax=0$, was einen Widerspruch zur positiven Definitheit der Matrix $A$ darstellte.

Beweis (ii)

(ii) Sei nun $B\in \mathbb {R} ^{(\ell +1)\times (\ell +1)}$ eine Untermatrix der Form (3.29) und $x:=(x_{i})\in \mathbb {R} ^{\ell +1}\setminus \{0\}$. Wegen der vorausgesetzten Symmetrie von $A$ ist auch $B$ symmetrisch. Für $z:=(z_{i})\in \mathbb {R} ^{n}$ mit

$z_{i}:={\begin{cases}x_{i+1-k},&k\leq i\leq k+\ell \\0,&{\text{sonst}}\end{cases}}$

gilt dann $z\neq 0$ sowie

$x^{T}Bx=\sum _{i,j=k}^{k+\ell }a_{ij}x_{i+1-k}x_{j+1-k}=\sum _{i,j=k}^{k+\ell }a_{ij}z_{i}z_{j}=\sum _{i,j=1}^{n}a_{ij}z_{i}z_{j}=z^{T}Az>0.$

Beweis (iii)

(iii) Die Eigenwerte $\lambda _{1},...,\lambda _{n}$ von $A$ sind (reell und) positiv, denn für $x^{(k)}\in \mathbb {R} ^{n}\setminus \{0\}$ mit $Ax^{(k)}=\lambda _{k}x^{(k)}$ gilt ja

$0<(x^{(k)})^{T}Ax^{(k)}=\lambda _{k}\left[(x^{(k)})^{T}x{(k)}\right],\quad k=1,...,n$

(siehe auch Lemma 3.19). Weiter ist die Matrix $A$ nach einem Ergebnis aus der Linearen Algebra diagonalisierbar, d. h., es gibt eine reguläre Matrix $S\in \mathbb {R} ^{n\times n}$ mit

$SAS^{-1}=D,$

wobei $D$ die Matrix

$D:=\operatorname {diag} (\lambda _{1},...,\lambda _{n})\in \mathbb {R} ^{n\times n}$

ist. Somit folgt

$\det(A)=\det(S^{-1})\det(D)\det(S)=\det(D)=\prod _{k=1}^{n}\lambda _{k}>0.$

q.e.d.

Satz - Produktzerlegung mit unterer Dreiecksmatrix

Die Matrix $A\in \mathbb {R} ^{n\times n}$ sei positiv definit. Dann gibt es genau eine untere Dreiecksmatrix $L:=(l_{kj})\in \mathbb {R} ^{n\times n}$ mit $l_{kk}>0$ $(k=1,...,n)$ und

(3.30) $A=LL^{T}.$

Beweis.

Der Beweis wird mit vollständiger Induktion geführt. Für $n:=1$ ist eine positiv definite Matrix $A=(\alpha )\in \mathbb {R} ^{1\times 1}$ eine positive Zahl $\alpha >0$. Eine solche kann eindeutig in der Form

$\alpha =l\cdot l,\quad l={\sqrt {\alpha }}$

geschrieben werden. Wir nehmen nun an, dass die Behauptung für positiv definite Matrizen bis zur Dimension $n-1$ richtig ist und betrachten jetzt eine positiv definite Matrix $A\in \mathbb {R} ^{n\times n}$. Diese lässt sich mit der nach Lemma 3.20 positiv definiten Matrix $A_{n-1}\in \mathbb {R} ^{(n-1)\times (n-1)}$ und einem Vektor $s\in \mathbb {R} ^{n}$ in der Form

$A={\begin{pmatrix}A_{n-1}&s\\s^{T}&a_{nn}\end{pmatrix}}$

partitionieren, wobei $A_{n-1}$ nach der Induktionsvoraussetzung mittels einer eindeutig bestimmten unteren Dreiecksmatrix $L_{n-1}:=(l_{kj})\in \mathbb {R} ^{(n-1)\times (n-1)}$ mit $l_{kk}>0$ $(k=1,...,n-1)$ zerlegt werden kann in

(3.31) $A_{n-1}=L_{n-1}L_{n-1}^{T}.$

Für die gesuchte Matrix $L\in \mathbb {R} ^{n\times n}$ machen wir nun einen Ansatz der Form

$L:={\begin{pmatrix}L_{n-1}&0\\c^{T}&\alpha \end{pmatrix}}$

und versuchen wir $c\in \mathbb {R} ^{n-1}$ und $\alpha \in \mathbb {R}$ so zu bestimmen, dass

(3.32) $A={\begin{pmatrix}A_{n-1}&s\\s^{T}a_{nn}\end{pmatrix}}={\begin{pmatrix}L_{n-1}&0\\c^{T}&\alpha \end{pmatrix}}{\begin{pmatrix}L_{n-1}^{T}&c\\0&\alpha \end{pmatrix}}$

gilt. Gleichheit in (3.32) hat man nun wegen (3.31) genau dann, wenn

(3.33) $L_{n-1}c=s,$

(3.34) $c^{T}c+\alpha ^{2}=a_{nn}$

gilt. Die Gleichung (3.33) besitzt sicher eine eindeutige Lösung $c=L_{n-1}^{-1}s$, da $L_{n-1}\in \mathbb {R} ^{(n-1)\times (n-1)}$ als untere Dreiecksmatrix mit positiven Diagonalelementen regulär ist. Auch die zweite Gleichung (3.34) besitzt offenbar eine Lösung $\alpha \in \mathbb {C}$. Aufgrund von (3.32) gilt außerdem

$\det(A)=\det {\begin{pmatrix}L_{n-1}&0\\c^{T}&\alpha \end{pmatrix}}\det {\begin{pmatrix}L_{n-1}^{T}&c\\0&\alpha \end{pmatrix}}=\alpha ^{2}[\det(L_{n-1})]^{2},$

so dass wegen $\det(A)>0$ (vgl. Lemma 3.20) und $\det(L_{n-1})=\prod _{k=1}^{n}l_{kk}>0$ auch $\alpha ^{2}>0$ ist und somit die Gleichung (3.34) eine eindeutige Lösung $\alpha >0$ hat. Damit ist alles gezeigt.

q.e.d.

Die Zerlegung $A=LL^{T}$ einer positiv definiten Matrix $A:=(a_{ij})\in \mathbb {R} ^{n\times n}$ bezeichnet man als Cholesky-Zerlegung von $A$. Ein direkter Ansatz zu ihrer Bestimmung ist es, die Gleichungen $A=LL^{T}$ bzw. die Gleichungen

${\begin{pmatrix}a_{11}&a_{12}&...&a_{1n}\\a_{21}&a_{22}&...&a_{2n}\\\vdots &\vdots &\ddots &\vdots \\a_{n1}&a_{n2}&...&a_{nn}\end{pmatrix}}={\begin{pmatrix}l_{11}&&&\\l_{21}&l_{22}&&\\\vdots &\ddots &\ddots &\\l_{n1}&...&l_{n,n-1}&l_{nn}\end{pmatrix}}{\begin{pmatrix}l_{11}&l_{21}&...&l_{n1}\\&l_{22}&...&l_{n2}\\&&\ddots &\vdots \\&&&l_{nn}\end{pmatrix}}$

als $n(n+1)/2$ Bestimmungsgleichungen für die $n(n+1)/2$ gesuchten Größen $l_{ik}$ $(i\geq k)$ aufzufassen:

$a_{ik}=\sum _{j=1}^{k}l_{ij}l_{kj},\quad 1\leq k\leq i\leq n\quad (i=1,...,n).$

Spaltenweise Berechnung der Einträge der unteren Dreiecksmatrix $L\in \mathbb {R} ^{n\times n}$ aus diesen Gleichungen führt auf den folgenden Algorithmus:

Algorithmus 3 (Cholesky-Zerlegung)

(0) Gib eine positiv definite Matrix $A:=(a_{ij})\in \mathbb {R} ^{n\times n}$ und setze $k:=1$.

(1) Berechne

$l_{kk}:=\left(a_{kk}-\sum _{j=1}^{k-1}l_{kj}^{2}\right)^{1/2},$

$l_{ik}:=\left(a_{ik}-\sum _{j=1}^{k-1}l_{ij}l_{kj}\right)/l_{kk}\quad (i=k+1,...,n).$

(2) Falls $k=n$, stop!

(3) Setze $k:=k+1$ und gehe nach (1).

Beispiel 3.22

Gegeben sei die positiv definite Matrix

$A:=(a_{ij})={\begin{pmatrix}1&2&1\\2&5&2\\1&2&10\end{pmatrix}}.$

Dann errechnet man für den Spaltenindex $k:=1$ die Einträge

$l_{11}:={\sqrt {a_{11}}}={\sqrt {1}}=1,$

$l_{21}:=a_{21}/l_{11}=2/1=2,$

$l_{31}:=a_{31}/l_{11}=1/1=1,$

für den Spaltenindex $k:=2$ die Einträge

$l_{22}:={\sqrt {a_{22}-l_{21}^{2}}}={\sqrt {5-4}}=1,$

$l_{32}:=(a_{32}-l_{31}l_{21})/l_{22}=(2-1\cdot 2)/1=0$

und schließlich für den Spaltenindex $k:=3$ den Eintrag

$l_{33}:={\sqrt {a_{33}-l_{31}^{2}-l_{32}^{2}}}={\sqrt {10-1^{2}-0^{2}}}=3.$

Somit erhält man für $A$ die Cholesky-Zerlegung

$A={\begin{pmatrix}1&0&0\\2&1&0\\1&0&3\end{pmatrix}}{\begin{pmatrix}1&2&1\\0&1&0\\0&0&3\end{pmatrix}}.$

---

Eine Cholesky-Zerlegung erfordert insgesamt die folgende Anzahl von Multiplikationen, Divisionen und Berechnungen von Quadratwurzeln:

$\sum _{k=1}^{n}\left[k+\sum _{i=k+1}^{n}k\right]=\sum _{k=1}^{n}[k+(n-k)k]=\sum _{k=1}^{n}k+n\sum _{k=1}^{n}k-\sum _{k=1}^{n}k^{2}$

$={\frac {n(n+1)}{2}}+n{\frac {n(n+1)}{2}}-{\frac {n(n+1)(2n+1)}{6}}={\frac {n^{3}}{6}}\left(1+0\left({\frac {1}{n}}\right)\right).$

Dies sind etwa halb so viele wesentliche Rechenoperationen, wie sie der Gauß-Algorithmus bzw. eine direkte LR-Zerlegung für eine beliebige reguläre Matrix erfordern.

3.3.7 Bandmatrizen

Bei der Diskretisierung von gewöhnlichen und partiellen Differentialgleichungen oder auch der Berechnung der Momente kubischer Splines (vgl. Abschnitt 7) ergeben sich lineare Gleichungssysteme $Ax=b$, bei denen $A=(a_{ij})\in \mathbb {R} ^{n\times n}$ mit $\det(A)\neq 0$ eine Bandmatrix der Bandbreite $p+q+1$ ist, d. h. bei denen $A$ die Gestalt

(3.35) $A:={\begin{pmatrix}a_{11}&...&a_{1,q+1}&&&\\\vdots &\ddots &&\ddots &&\\a_{p+1,1}&&\ddots &&\ddots &\\&\ddots &&\ddots &&a_{n-q,n}\\&&\ddots &&\ddots &\vdots \\&&&a_{n,n-p}&...&a_{nn}\end{pmatrix}}$

hat und somit

$a_{ik}:=0,\quad i=1,...,n,\quad 1\leq k<i-p\leq n,\quad 1\leq i+q<k\leq n$

mit gewissen $p,q\in \{1,...,n-1\}$ gilt. Insbesondere spricht man im Fall $p=q=1$, d. h. im Fall

$A:={\begin{pmatrix}a_{11}&a_{12}&&&&\\a_{21}&a_{22}&a_{23}&&&\\&a_{32}&a_{33}&\ddots &&\\&&\ddots &\ddots &a_{n-2,n-1}&\\&&&a_{n-1,n-2}&a_{n-1,n-1}&a_{n-1,n}\\&&&&a_{n,n-1}&a_{nn}\end{pmatrix}}$

von einer Tridiagonalmatrix.

Bei Gleichungssystemen mit Bandmatrizen lässt sich der zu betreibende Aufwand bei allen in diesem Kapitel angesprochenen Methoden verringern, außer bei denen mit Pivotstrategien, da diese die Bandstruktur im Allgemeinen zerstören. Exemplarisch soll das Vorgehen für Bandmatrizen am Beispiel der direkten LR-Zerlegung demonstriert werden. Wenn eine LR-Zerlegung für $A$ in (3.35) möglich ist (und $\det(A)\neq 0$ ist), so ist diese im Fall, dass man die Diagonaleinträge von $L$ als 1 wählt, eindeutig und von der Gestalt

$A={\begin{pmatrix}1&&&&&\\l_{21}&\ddots &&&&\\\vdots &\ddots &\ddots &&&\\l_{p+1,1}&&\ddots &\ddots &&\\&\ddots &&\ddots &\ddots &\\&&l_{n,n-p}&...&l_{n,n-1}&1\end{pmatrix}}{\begin{pmatrix}r_{11}&...&r_{1,q+1}&&\\&\ddots &&\ddots &\\&&\ddots &&r_{n-q,n}\\&&&\ddots &\vdots \\&&&&r_{nn}\end{pmatrix}}$

(siehe z. B. Satz 2.29 bei Kanzow). Komponentenschreibweise geschrieben heißt dies

$a_{ik}=\sum _{j=j_{0}}^{\min\{i,k\}}l_{ij}r_{jk},\quad i=1,...,n,\quad k=\max\{1,i-p\},...,\min\{n,i+q\},\quad j_{0}:=\max\{1,i-p,k-q\},$

was bei einer Parkettierung wie in (3.28) auf den folgenden Algorithmus zur Bestimmung der LR-Zerlegung der Bandmatrix $A$ führt:

Algorithmus 4 (LR-Zerlegung für Bandmatrizen)

(0) Gib eine Matrix $A:=(a_{ij})\in \mathbb {R} ^{n\times n}$ mit

$a_{ik}:=0,\quad i=1,...,n,\quad 1\leq k<i-p\leq n,\quad 1\leq i+q<k\leq n$

für gegebene $p,q\in \{1,...,n-1\}$ und setze $i:=1$.

(1) Für $k=i,...,\min\{i+q,n\}$ berechne $j_{0}=\max\{1,i-p,k-q\}$ und

$r_{ik}:=a_{ik}-\sum _{j=j_{0}}^{i-1}l_{ij}r_{jk}.$

(2) Für $k=i+1,...,\min\{i+p,n\}$ berechne $j_{0}=\max\{1,i-p,k-q\}$ und

$l_{ki}:=\left(a_{ki}-\sum _{j=j_{0}}^{i-1}l_{kj}r_{ji}\right)/r_{ii}.$

(3) Falls $i=n$, stop!

(4) Setze $i:=i+1$ und gehe nach (1).

Ist $A$ eine Tridiagonalmatrix und schreibt man

(3.36) ${\begin{pmatrix}\alpha _{1}&\beta _{2}&&\\\gamma _{2}&\alpha _{2}&\ddots &\\&\ddots &\ddots &\beta _{n}\\&&\gamma _{n}&\alpha _{n}\end{pmatrix}}={\begin{pmatrix}1&&&\\l_{2}&1&&\\&\ddots &\ddots &\\&&l_{n}&1\end{pmatrix}}{\begin{pmatrix}d_{1}&r_{2}&&\\&\ddots &\ddots &\\&&\ddots &r_{n}\\&&&d_{n}\end{pmatrix}},$

so vereinfacht sich Algorithmus 4 zu

Algorithmus 4* (LR-Zerlegung Tridiagonalmatrizen)

(0) Gib eine Tridiagonalmatrix $A\in \mathbb {R} ^{n\times n}$ und schreibe $A,L$ und $R$ wie in (3.36). Setze

$d_{1}:=\alpha _{1},\quad r_{2}:=\beta _{2}$

und $i:=2$.

(1) Berechne

$l_{i}:=\gamma _{i}/d_{i-1},\quad d_{i}:=\alpha _{i}-l_{i}r_{i},\quad r_{i+1}:=\beta _{i+1}.$

(2) Falls $i=n-1$, berechne

$l_{n}:=\gamma _{n}/d_{n-1},\quad d_{n}:=\alpha _{n}-l_{n}r_{n}$

und stoppe!

(3) Setze $i:=i+1$ und gehe nach (1).

Man kann zeigen, dass im Fall einer Tridiagonalmatrix eine LR-Zerlegung wie in (3.36) möglich ist, wenn gilt (Lemma 2.28 bei Kanzow):

$|\alpha _{1}|>|\beta _{2}|,\quad |\alpha _{i}|\geq |\gamma _{i}|+|\beta _{i+1}|,\quad i=2,...,n-1,$

$|\alpha _{n}|\geq |\gamma _{n}|,\quad \gamma _{i}\neq 0,\quad i=2,...,n.$

Diese Bedingungen besagen offenbar, dass man für die erste Zeile strikte und für die anderen Zeilen nur normale Diagonaldominanz fordern muss. Die Forderung $\gamma _{i}\neq 0$ $(i=2,...,n)$ macht Sinn, da im anderen Fall eine LR-Zerlegung mit $L:=I$ und $R:=A$ existiert und folglich nicht berechnet werden muss. Für die LR-Zerlegung einer Tridiagonalmatrix sind offenbar nur

$(n-2)\cdot 2+2=2n-2$

wesentliche Rechenoperationen erforderlich.

Siehe auch

• Lösung linearer Gleichungssysteme - (Foliensatz)
• Orthonormalisierungsverfahren - (Foliensatz)

Seiteninformation

Diese Lernresource können Sie als Wiki2Reveal-Foliensatz darstellen.

Wiki2Reveal

Dieser Wiki2Reveal Foliensatz wurde für den Lerneinheit Kurs:Numerik I' erstellt der Link für die Wiki2Reveal-Folien wurde mit dem Wiki2Reveal-Linkgenerator erstellt.

• Die Seite wurde als Dokumententyp PanDocElectron-SLIDE erstellt.
• Link zur Quelle in Wikiversity: https://de.wikiversity.org/wiki/Kurs:Numerik%20I/Zerlegung%20PA%20=%20LR
• siehe auch weitere Informationen zu Wiki2Reveal und unter Wiki2Reveal-Linkgenerator.