In diesem Abschnitt werden "direkte" Verfahren zur numerischen Lösung linearer Gleichungssysteme $Ax=b$ vorgestellt, wobei $A\in \mathbb {R} ^{n\times n}$ eine gegebene Matrix und $b\in \mathbb {R} ^{n}$ ein gegebener Vektor ist.

3.1 Dreieckssysteme

Definition 3.1

Matrizen $L,R\in \mathbb {R} ^{n\times n}$ der Form

L:={\begin{pmatrix}l_{11}&0&...&0\\l_{21}&l_{22}&\ddots &\vdots \\\vdots &&\ddots &0\\l_{n1}&...&...&l_{nn}\end{pmatrix}},\quad R:={\begin{pmatrix}r_{11}&r_{12}&...&r_{1n}\\0&r_{22}&&\vdots \\\vdots &\ddots &\ddots &\vdots \\0&...&0&r_{nn}\end{pmatrix}}

heißen untere bzw. obere Dreiecksmatrizen.

Die Matrizen $L$ und $R$ in der Definition sind offenbar genau dann regulär, wenn gilt:

\det(L)=\prod _{k=1}^{n}l_{kk}\neq 0,\quad \det(R)=\prod _{k=1}^{n}r_{kk}\neq 0.

Für die obere Dreiecksmatrix $R:=(r_{kj})\in \mathbb {R} ^{n\times n}$ mit $r_{kj}=0$ für $k>j$ ist das gestaffelte Gleichungssystem $Rx=z$ für einen gegebenen Vektor $z\in \mathbb {R} ^{n}$ von der Form

{\begin{array}{ccccccccc}r_{11}x_{1}&+&r_{12}x_{2}&+&...&+&r_{1n}x_{n}&=&z_{1},\\&&r_{22}x_{2}&+&...&+&r_{2n}x_{n}&=&z_{2},\\&&&&\ddots &&\vdots &\vdots &\vdots \\&&&&&&r_{nn}x_{n}&=&z_{n}.\end{array}}

Dessen Lösung $x\in \mathbb {R} ^{n}$ kann für reguläres $R$ , d. h. für $r_{ii}\neq 0$ $(i=1,...,n)$ , zeilenweise von unten nach oben durch Auflösen nach der jeweiligen Unbekannten auf der Diagonalen berechnet werden und zwar nach der folgenden Vorschrift:

Für $i=n,n-1,n-2,...,1$ berechne

x_{i}=\left(z_{i}-\sum _{j=i+1}^{n}r_{ij}x_{j}\right)/r_{ii}.

Dabei sind in den Stufen $i=n,n-1,n-2,...,1$ je $n-i$ Multiplikationen und eine Division, d. h. $n-i+1$ wesentliche arithmetische Operationen, durchzuführen. Insgesamt sind dies also

(3.1)

\sum _{i=1}^{n}(n-i+1)=\sum _{m=1}^{n}m={\frac {n(n+1)}{2}}={\frac {n^{2}}{2}}\left(1+{\mathcal {O}}\left({\frac {1}{n}}\right)\right)

Multiplikationen und Divisionen.

Für eine untere Dreiecksmatrix $L:=(l_{ij})\in \mathbb {R} ^{n\times n}$ mit $l_{ij}=0$ für $i<j$ ist das entsprechende gestaffelte Gleichungssystem $Lx=b$ für einen gegebenen Vektor $b\in \mathbb {R} ^{n}$ von der Form

{\begin{array}{ccccccccc}l_{11}x_{1}&&&&&&&=&b_{1},\\l_{21}x_{1}&+&l_{22}x_{2}&&&&&=&b_{2},\\\vdots &&\vdots &&\ddots &&&\vdots &\vdots \\l_{n1}x_{1}&+&l_{n2}x_{2}&+&...&+&_{nn}x_{n}&=&b_{n}.\end{array}}

Seine Lösung $x\in \mathbb {R} ^{n}$ kann für reguläres $L$ , d. h. für $l_{ii}\neq 0$ $(i=1,...,n)$ , zeilenweise von oben nach unten durch Auflösen nach der Unbekannten in der Diagonalen berechnet werden und zwar nach folgender Vorschrift:

Für $i=1,2,3,...,n$ berechne

x_{i}=\left(b_{i}-\sum _{j=1}^{i-1}l_{ij}x_{j}\right)/l_{ii}.

Dabei sind genauso viele wesentliche arithmetische Rechenoperationen durchzuführen wie im Fall eines oberen gestaffelten Gleichungssystems mit $n$ Veränderlichen, nämlich $(n^{2}/2)+(n/2)$ (siehe (3.1)).

3.2 Der Gauß-Algorithmus

3.2.1 Einführung

Seien nun $A:=(a_{ij})\in \mathbb {R} ^{n\times n}$ eine gegebene Matrix mit $\det(A)\neq 0$ und $b\in \mathbb {R} ^{n}$ ein gegebener Vektor. Im Folgenden beschreiben wir den sog. Gauß-Algorithmus. Dieser führt das lineare Gleichungssystem (kurz: LGS) $Ax=b$ in ein äquivalentes oberes gestaffeltes LGS $Rx=z$ über, dessen Lösung $x\in \mathbb {R} ^{n}$ , wie im vorangehenden Abschnitt gezeigt wurde, leicht berechnet werden kann. Erlaubte Operationen, die zu einer äquivalenten Umformung eines LGSs führen, sind offenbar

die Vertauschung der Reihenfolge von Gleichungen,
die skalare Multiplikation von Gleichungen,
die Addition des skalaren Vielfachen einer Gleichung zu einer anderen Gleichung,
die Vertauschung von Spalten der Koeffizientenmatrix.

Dabei entspricht offenbar die Vertauschung von Spalten der Koeffizientenmatrix einer Vertauschung der Reihenfolge, in der die Variablen im LGS auftreten. Wir führen im Folgenden den Gauß-Algorithmus zunächst ohne Zeilen- und Spaltenvertauschungen ein, in welchem Fall er aber auch nur für spezielle Matrizen (wir geben einen solchen Typ an) durchführbar ist.

In der ersten Stufe des Gauß-Algorithmus wird das gegebene LGS

{\begin{array}{ccccccccc}a_{11}x_{1}&+&a_{12}x_{2}&+&...&+&a_{1n}x_{n}&=&b_{1},\\a_{21}x_{1}&+&a_{22}x_{2}&+&...&+&a_{2n}x_{n}&=&b_{2},\\\vdots &&\vdots &&&&\vdots &&\vdots \\a_{n1}x_{1}&+&a_{n2}x_{2}&+&...&+&a_{nn}x_{n}&=&b_{n}\end{array}}

in ein äquivalentes LGS der Form

{\begin{array}{ccccccccc}a_{11}^{(2)}x_{1}&+&a_{12}^{(2)}x_{2}&+&...&+&a_{1n}^{(2)}x_{n}&=&b_{1}^{(2)},\\&&a_{22}^{(2)}x_{2}&+&...&+&a_{2n}^{(2)}x_{n}&=&b_{2}^{(2)},\\&&\vdots &&&&\vdots &&\vdots \\&&a_{n2}^{(2)}x_{2}&+&...&+&a_{nn}^{(2)}x_{n}&=&b_{n}^{(2)}\end{array}}

überführt. Wenn $a_{11}\neq 0$ ist, kann dies für die erste Zeile mit

a_{1j}^{(2)}:=a_{1j}\quad (j=1,...,n),\quad b_{1}^{(2)}:=b_{1}

erreicht werden und für die Zeilen $i=2,3,...,n$ mit Zeilenoperationen der Form

{\text{neue Zeile }}i:={\text{alte Zeile }}i-l_{i1}\cdot ({\text{alte Zeile }}1)

für

l_{i1}:={\frac {a_{i1}}{a_{11}}},\quad i=2,3,...,n

bzw. explizit durch

\underbrace {(a_{i1}-l_{i1}a_{11})} _{=0}x_{1}+\underbrace {(a_{i2}-l_{i1}a_{12})} _{=:a_{i2}^{(2)}}x_{2}+...+\underbrace {(a_{in}-l_{i1}a_{1n})} _{=:a_{in}^{(2)}}x_{n}=\underbrace {b_{i}-l_{i1}b_{1}} _{=:b_{i}^{(2)}}.

Diesen Eliminationsschritt wiederholt man nun analog für das um die erste Zeile und Spalte reduzierte LGS für die $(n-1)$ Unbekannten $x_{2},...,x_{n}$ . (Denn Addition eines Vielfachen einer Zeile mit Index $\geq 2$ zu einer anderen Zeile mit Index $\geq 2$ erzeugt wiederum eine Null in der ersten Spalte.) Führt man diesen Eliminationsprozess sukzessive für die jeweils entstehenden Teilsysteme durch, so erhält man also im Fall $a_{kk}^{(k)}\neq 0$ zu $Ax=b$ äquivalente Gleichungssysteme

A^{(k)}x=b^{(k)},\quad k=1,...,n

mit Matrizen der speziellen Gestalt

(3.2)

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

Dabei hat man die folgenden Transformationen zu berechnen:

A:=A^{(1)}\to A^{(2)}\to ...\to A^{(n)}=:R,

b:=b^{(1)}\to b^{(2)}\to ...\to b^{(n)}=:z,

wobei $R$ obere Dreiecksmatrix und $A^{(s)}=A^{(s+1)}=...=A^{(n)}=R$ für $s\leq n$ möglich ist. Den Übergang von $A^{(i)}$ nach $A^{(i+1)}$ bezeichnen wir als $i$ -te Stufe des Gauß-Algorithmus.

Im Zuge des Gauß-Algorithmus sind in der $r$ -ten Stufe $(n-r)(n-r+1)$ Multiplikationen und $n-r$ Divisionen, d. h. $(n-r)^{2}+2(n-r)$ wesentliche arithmetische Rechenoperationen durchzuführen, so dass insgesamt maximal

\sum _{i=1}^{n-1}i^{2}+2\sum _{i=1}^{n-1}i={\frac {n(n+1)(2n+1)}{6}}+n(n-1)={\frac {n^{3}}{3}}\left(1+0\left({\frac {1}{n}}\right)\right)

wesentliche Rechenoperationen anfallen.

Es wurde hier vorausgesetzt, dass die in jeder Stufe auftretenden Diagonalelemente nicht verschwinden, d. h. $a_{kk}^{(k)}\neq 0$ $(k=1,...,n-1)$ ist. Der folgende Satz gibt eine Klasse von Matrizen $A\in \mathbb {R} ^{n\times n}$ an, für die dies der Fall und somit der Gauß-Algorithmus durchführbar ist.

Definition 3.2

Eine Matrix $A=(a_{ij})\in \mathbb {R} ^{n\times n}$ heißt strikt diagonaldominant, wenn gilt:

\sum _{j=1 \atop j\neq i}^{n}|a_{ij}|<|a_{ii}|,\quad i=1,...,n.

Satz 3.3

Wenn $A=(a_{ij})\in \mathbb {R} ^{n\times n}$ strikt diagonaldominant ist, so ist der Gauß-Algorithmus zur Lösung von $Ax=b$ durchführbar.

Beweis.

Es wird mit vollständiger Induktion über $k=1,...,n-1$ gezeigt, dass die Matrizen

B^{(k)}:={\begin{pmatrix}a_{kk}^{(k)}&...&a_{kn}^{(k)}\\\vdots &\ddots &\vdots \\a_{nk}^{(k)}&...&a_{nn}^{(k)}\end{pmatrix}}\in \mathbb {R} ^{(n-k+1)\times (n-k+1)}

strikt diagonaldominant sind und damit insbesondere $|a_{kk}^{(k)}|>0$ $(k=1,...,n-1)$ gilt. Für $B^{(1)}:=A$ ist dies nach Voraussetzung richtig und wir nehmen nun an, dass $B^{(k)}$ für ein beliebiges $k$ strikt diagonaldominant ist. Dann gilt insbesondere $a_{kk}^{(k)}\neq 0$ , so dass der Gauß-Eliminationsschritt auf $B^{(k)}$ anwendbar ist. Er liefert die Matrix $B^{(k+1)}=(a_{ij}^{(k+1)})\in \mathbb {R} ^{(n-k)\times (n-k)}$ mit

a_{ij}^{(k+1)}:=a_{ij}^{(k)}-l_{ik}a_{kj}^{(k)},\quad i,j=k+1,...,n

und den Koeffizienten

l_{ik}=a_{ik}^{(k)}/a_{kk}^{(k)},\quad i=k+1,...,n.

(Man beachte, dass $l_{ik}$ nicht beim nächsten Schritt überschrieben wird und somit nicht mit einem Iterationszähler versehen werden muss.) Für $i=k+1,...,n$ ergibt sich unter Verwendung der Induktionsvoraussetzung

\sum _{j=k+1 \atop j\neq i}^{n}\left|a_{ij}^{(k+1)}\right|\leq \sum _{j=k+1 \atop j\neq i}^{n}\left|a_{ij}^{(k)}\right|+|l_{ik}|\sum _{j=k+1 \atop j\neq i}^{n}\left|a_{kj}^{(k)}\right|=\sum _{j=k \atop j\neq i}^{n}\left|a_{ij}^{(k)}\right|-\left|a_{ik}^{(k)}\right|+|l_{ik}|\left(\sum _{j=k+1}^{n}\left|a_{kj}^{(k)}\right|-\left|a_{ki}^{(k)}\right|\right)

<\left|a_{ii}^{(k)}\right|-\left|a_{ik}^{(k)}\right|+{\frac {\left|a_{ik}^{(k)}\right|}{\left|a_{kk}^{(k)}\right|}}\left(\left|a_{kk}^{(k)}\right|-\left|a_{ki}^{(k)}\right|\right)=\left|a_{ii}^{(k)}\right|-|l_{ik}|\left|a_{ki}^{(k)}\right|\leq \left|a_{ii}^{(k+1)}\right|.

3.2.2 Gauß-Algorithmus mit Pivotsuche

Wir betrachten nun zunächst das folgende Beispiel.

Beispiel 3.4

(1) Für

A:={\begin{pmatrix}0&1\\1&0\end{pmatrix}}

besitzt das Gleichungssystem $Ax=b$ wegen $\det(A)=-1$ für jedes $b\in \mathbb {R} ^{2}$ eine eindeutige Lösung. Jedoch ist der Gauß-Algorithmus in der angegebenen Form wegen $a_{11}=0$ nicht für dessen Lösung durchführbar. Vertauscht man jedoch die Zeilen in dem System und damit in $A$ , so erhält man die Matrix

{\bar {A}}:={\begin{pmatrix}1&0\\0&1\end{pmatrix}}

und kann der Gauß-Algorithmus für die Lösung des zugehörigen Systems erfolgreich angewendet werden.

(2) Die exakte Lösung des Gleichungssystems

(3.3)

{\begin{pmatrix}10^{-4}&1\\1&1\end{pmatrix}}{\begin{pmatrix}x_{1}\\x_{2}\end{pmatrix}}={\begin{pmatrix}1\\2\end{pmatrix}}

lautet $x_{1}=1.000\ 100\ 01,x_{2}=0.999\ 899\ 99$ . Bei 3-stelliger Rechnung ergibt sich mit dem Gauß-Algorithmus der Umrechnungsfaktor $l_{21}=10^{4}$ , damit das Tableau

\left({\begin{array}{cc|c}10^{-4}&1&1\\0&-10^{4}&-10^{4}\end{array}}\right)

und die mit großen Fehlern behaftete "Lösung" $x_{1}=0,x_{2}=1$ . Vertauscht man aber die Zeilen in (3.3), so gelangt man mit ${\tilde {l}}_{21}=10^{-4}$ zu dem Tableau

\left({\begin{array}{cc|c}1&1&2\\0&1&1\end{array}}\right)

und damit zu der guten Näherungslösung $x_{1}=1,x_{2}=1$ .

Insbesondere können also große Umrechnungsfaktoren zu numerischen Instabilitäten führen. Zur Vermeidung solcher etwaigen Instabilitäten bietet sich die folgende Modifikation des Gauß-Algorithmus mit einer Spaltenpivotsuche an. Dabei nimmt man, sofern dies erforderlich ist, im $k$ -ten Schritt eine Zeilenvertauschung derart vor, dass man das auf bzw. unterhalb der Diagonale befindliche betragsmäßig größte Element der aktuellen $k$ -ten Spalte an die Position des $k$ -ten Diagonalelementes bringt.

Algorithmus 1 (k-te Stufe des Gauß-Algorithmus mit Spaltenpivotsuche)

(0) Gib die Matrix $A^{(k)}$ der Gestalt (3.2).

(1) Bestimme ein $p\in \{k,...,n\}$ mit

\left|a_{pk}^{(k)}\right|\geq \left|a_{ik}^{(k)}\right|,\quad i=k,...,n.

(2) Erzeuge ${\hat {A}}^{(k)}:=({\hat {a}}_{ij}^{(k)})\in \mathbb {R} ^{n\times n}$ aus $A^{(k)}$ sowie ${\hat {b}}^{(k)}:=({\hat {b}}_{i}^{(k)})\in \mathbb {R} ^{n}$ aus $b^{(k)}$ durch Vertauschung der $p$ -ten und der $k$ -ten Zeile von $A^{(k)}$ bzw. $b^{(k)}$ , d. h. setze

\left.{\begin{matrix}{\hat {a}}_{pj}^{(k)}:=a_{kj}^{(k)},\quad {\hat {a}}_{kj}^{(k)}:=a_{pj}^{(k)},\\{\hat {a}}_{ij}^{(k)}:=a_{ij}^{(k)}{\mbox{ für }}i\notin \{p,k\}\end{matrix}}\right\}j=1,...,n

sowie

{\begin{matrix}{\hat {b}}_{k}^{(k)}:=b_{p}^{(k)},\quad {\hat {b}}_{p}^{(k)}:=b_{k}^{(k)},\\{\hat {b}}_{i}^{(k)}:=b_{i}^{(k)}{\mbox{ für }}i\notin \{p,k\}.\end{matrix}}

(3) Führe den Eliminationsschritt ${\hat {A}}^{(k)}\to A^{(k+1)},{\hat {b}}^{(k)}\to b^{(k+1)}$ wie oben für $A^{(k)}\to A^{(k+1)}$ und $b^{(k)}\to b^{(k+1)}$ beschrieben aus.

(4) Setze $k:=k+1$ .

Das Element $a_{pk}^{(k)}$ wird als Pivotelement bezeichnet. Für $\det(A)\neq 0$ ist $a_{pk}^{(k)}\neq 0$ für jedes $k$ , d. h. ist der Gauß-Algorithmus mit Spaltenpivotsuche durchführbar. Denn in der $k$ -ten Stufe des Verfahrens muss $a_{ik}^{(k)}\neq 0$ für mindestens ein $i\in \{k,k+1,...,n\}$ gelten, da man anderenfalls zur nächsten Stufe übergehen könnte und damit ${\hat {A}}^{(n)}$ eine Dreiecksmatrix wäre, für die mindestens ein Diagonalelement identisch Null und somit $\det({\hat {A}}^{(n)})=0$ wäre. Letzteres ist jedoch wegen $\left|\det({\hat {A}}^{(n)})\right|=|\det(A)|$ nicht möglich, da die beim Gauß-Algorithmus mit Pivotsuche in jeder Stufe durchgeführten Operationen (Zeilenvertauschung und Addition eines Vielfachen einer Zeile zu einer anderen Zeile) höchstens einen Vorzeichenwechsel der Determinante zur Folge haben.

Es sei darauf hingewiesen, dass man offenbar durch Multiplikation der entsprechenden Gleichung mit einem geeigneten Skalar jedes Element $\neq 0$ in der Pivotspalte zum Pivotelement machen kann und somit eine geeignete Skalierung des zugrunde liegenden linearen Gleichungssystems den Verlauf des Gauß-Algorithmus wesentlich beeinflussen kann. Mehr dazu findet man z. B. bei Deuflhard/Hohmann. In dem Buch dieser Autoren findet man auch Aussagen zur Stabilität des Gauß-Algorithmus. So wird festgestellt, dass Gauß-Elimination mit Spaltenpivotsuche über der Menge aller invertierbaren Matrizen betrachtet nicht stabil, für wichtige Unterklassen, wie beispielsweise die der positiv definiten Matrizen (vgl. Abschnitt 3.3.6) und die der strikt diagonaldominanten Matrizen, aber stabil ist.

3.3 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. In einem solchen Fall ist es ineffizient, für jede neue rechte Seiten wieder den Gauß-Algorithmus durchzuführen, da er bei jedem Durchlauf in Bezug auf $A$ dieselben Resultate liefern würde. Deshalb möchte man die beim Gauß-Algorithmus durchgeführten Transformationen von $A$ in irgendeiner Form speichern. 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. Eine solche Zerlegung bzw. Faktorisierung von $A$ kann man, wie wir zeigen wollen, mittels des Gauß-Algorithmus mit Spaltenpivotsuche gewinnen. 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

(3.4)

Ly=Pb,\quad Rx=y

löst, wobei die Lösung $y$ des ersten Systems die rechte Seite des zweiten Systems liefert. 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äß (3.4) numerisch relativ "billig" ist.

3.3.1 Permutationsmatrizen

Definition 3.5

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 3.6

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 3.7

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

Bemerkung 3.8

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

(3.5)

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 3.9

Der Matrix aus Beispiel 3.6 liegt die durch

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

definierte Permutation $\pi :\{1,...,n\}\to \{1,...,n\}$ zugrunde, denn es gilt

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 Bemerkung 3.8 kann man nun die Aussagen des folgenden Satzes erschließen.

Satz 3.10

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 3.12

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$ .

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.

3.3.6 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 3.18

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:

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

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.

Lemma 3.19

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.

Den Beweis des Lemmas findet man in Büchern der Linearen Algebra. Wir beweisen hier:

Lemma 3.20

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.

(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.

(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.

(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 3.21

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.

3.4 Orthogonalisierungsverfahren

In diesem Abschnitt soll für eine gegebene Matrix $A\in \mathbb {R} ^{n\times k}$ mit $1\leq k\leq n$ eine Zerlegung der Form

(3.37)

A=QS

bestimmt werden, wobei $Q\in \mathbb {R} ^{n\times n}$ eine orthogonale Matrix, d. h. eine reguläre Matrix mit

(3.38)

Q^{-1}=Q^{T}

ist und $S$ mit einer oberen Dreiecksmatrix $R$ wie folgt gebildet wird:

(3.39)

S={\begin{pmatrix}R\\\mathbf {0} \end{pmatrix}}\in \mathbb {R} ^{n\times k},\quad R:={\begin{pmatrix}\ast &...&*\\&\ddots &\vdots \\&&*\end{pmatrix}}\in \mathbb {R} ^{k\times k},\quad \mathbf {0} :=(0)\in \mathbb {R} ^{(n-k)\times k}.

Im Fall $n=k$ ist demnach insbesondere $S=R$ . Wie wir zeigen wollen, ermöglicht eine solche QR- bzw. QS-Zerlegung von $A$ sowohl die stabile Lösung linearer Gleichungssysteme $Ax=b$ als auch die stabile Lösung von Ausgleichsproblemen

\min _{x\in \mathbb {R} ^{k}}\|b-Ax\|_{2}.

Dafür benötigen wir einige Eigenschaften orthogonaler Matrizen.

3.4.1 Elementare Eigenschaften orthogonaler Matrizen

Lemma 3.23

Sei $Q\in \mathbb {R} ^{n\times n}$ eine orthogonale Matrix. Dann ist auch $Q^{T}$ eine orthogonale Matrix und es gilt

(3.40)

\|Qx\|_{2}=\|x\|_{2}=\left\|Q^{T}x\right\|_{2},\quad x\in \mathbb {R} ^{n}.

Beweis.

Wegen

(Q^{T})^{-1}=(Q^{-1})^{-1}=Q=(Q^{T})^{T}

ist auch $Q^{T}$ eine orthogonale Matrix. Des Weiteren hat man für $Q$

\|Qx\|_{2}^{2}=(Qx)^{T}Qx=x^{T}Q^{T}Qx=x^{T}Q^{-1}Qx=x^{T}x=\|x\|_{2}^{2}

und somit auch $\left\|Q^{T}x\right\|_{2}=\|x\|_{2}$ für die orthogonale Matrix $Q^{T}$ .

q.e.d.

Die Beziehungen (3.40) besagen, dass die durch $Q$ und $Q^{T}$ definierten linearen Abbildungen isometrisch bezüglich der Euklidischen Norm sind.

Lemma 3.24

Es gilt:

(i) $\|Q\|_{2}=\left\|Q^{T}\right\|_{2}=1$ ,

(ii) $|\det(Q)|=\left|\det(Q^{T})\right|=1$ .

Beweis.

Lemma 3.23, angewandt auf $Q$ und $Q^{T}$ , impliziert Aussage (i) wegen

\max _{\|x\|_{2}=1}\|Qx\|_{2}=\max _{\|x\|_{2}=1}\|x\|_{2}=1=\max _{\|x\|_{2}=1}\left\|Q^{T}x\right\|_{2}.

Aussage (ii) folgt aus

1=\det(I)=\det(QQ^{T})=\det(Q)\det(Q^{T})=[\det(Q)]^{2}=\left[\det(Q^{T})\right]^{2}.

q.e.d.

Korollar 3.25

Seien $A\in \mathbb {R} ^{n\times n}$ eine reguläre und $Q\in \mathbb {R} ^{n\times n}$ eine orthogonale Matrix. Dann gilt

\operatorname {cond} _{2}(QA)=\operatorname {cond} _{2}(A).

Beweis.

Nach Lemma 3.23 gilt $\|QAx\|_{2}=\|Ax\|_{2}$ für $x\in \mathbb {R} ^{n}$ und somit

\|A\|_{2}=\max _{\|x\|_{2}=1}\|Ax\|_{2}=\max _{\|x\|_{2}=1}\|QAx\|2=\|QA\|_{2}.

Weiter gilt mit Lemma 3.24 (i)

\left\|A^{-1}\right\|_{2}=\left\|A^{-1}Q^{T}Q\right\|_{2}\leq \left\|A^{-1}Q^{T}\right\|_{2}\|Q\|_{2}=\left\|A^{-1}Q^{T}\right\|_{2}\leq \left\|A^{-1}\right\|_{2}\left\|Q^{T}\right\|_{2}=\left\|A^{-1}\right\|_{2}

und folglich $\|A^{-1}\|_{2}=\left\|A^{-1}Q^{T}\right\|_{2}$ . Damit ergibt sich

\operatorname {cond} _{2}(QA)=\|QA\|_{2}\left\|A^{-1}Q^{-1}\right\|_{2}=\|QA\|_{2}\left\|A^{-1}Q^{T}\right\|_{2}=\|A\|_{2}\left\|A^{-1}\right\|_{2}=\operatorname {cond} _{2}(A).

q.e.d.

Multiplikation einer quadratischen Matrix $A$ von links mit einer orthogonalen Matrix führt also auf eine Matrix mit derselben Kondition wie $A$ . Weiter gilt:

Lemma 3.26

Seien $Q_{1},Q_{2}\in \mathbb {R} ^{n\times n}$ orthogonale Matrizen. Dann ist auch $Q_{1}Q_{2}$ eine orthogonale Matrix.

Beweis.

Die Aussage des Lemma folgt aus

(Q_{1}Q_{2})^{-1}=Q_{2}^{-1}Q_{1}^{-1}=Q_{2}^{T}Q_{1}^{T}=(Q_{1}Q_{2})^{T}.

q.e.d.

3.4.2 Lösung linearer Gleichungssysteme mittels QR-Zerlegung

Wir betrachten nun wieder die Lösung eines linearen Gleichungssystems $Ax=b$ für eine reguläre Matrix $A\in \mathbb {R} ^{n\times n}$ und beliebige rechte Seite $b\in \mathbb {R} ^{n}$ und gehen davon aus, dass eine Zerlegung der Form $A=QR$ von $A$ mit einer orthogonalen Matrix $Q\in \mathbb {R} ^{n\times n}$ und einer oberen Dreiecksmatrix $R\in \mathbb {R} ^{n\times n}$ gegeben ist (vgl. (3.37) - (3.39)). Wegen Lemma 3.24 gilt offenbar

\det(A)\neq 0\Leftrightarrow \det(R)\neq 0.

Weiter hat man die Äquivalenzen

Ax=b\Leftrightarrow QRx=b\Leftrightarrow Rx=Q^{T}b

und mit Korollar 3.25 die Beziehungen

\operatorname {cond} _{2}(R)=\operatorname {cond} _{2}(Q^{T}A)=\operatorname {cond} _{2}(A).

D. h., das Gleichungssystem $Ax=b$ ist äquivalent zu dem gestaffelten Gleichungssystem $Rx=Q^{T}b$ , wobei die Matrix $R$ bezüglich der Spektralnorm keine schlechtere Kondition als $A$ aufweist und gemäß 3.40 $\left\|Q^{T}b\right\|_{2}=\|b\|_{2}$ gilt.

3.4.3 QR-Zerlegung mittels Gram-Schmidt-Orthogonalisierung

Es sei nun $A\in \mathbb {R} ^{n\times n}$ eine quadratische Matrix mit $\det(A)\neq 0$ und es seien eine orthogonale Matrix $Q\in \mathbb {R} ^{n\times n}$ und eine obere Dreiecksmatrix $R\in \mathbb {R} ^{n\times n}$ gesucht, so dass

(3.41)

A=QR

gilt. Schreiben wir $A,Q$ und $R$ mit Vektoren $a^{j},q^{j}\in \mathbb {R} ^{n}$ in der Form

A={\begin{pmatrix}a^{1}&...&a^{n}\end{pmatrix}},\quad Q={\begin{pmatrix}q^{1}&...&q^{n}\end{pmatrix}},\quad R={\begin{pmatrix}r_{11}&...&r_{1n}\\&\ddots &\vdots \\&&r_{nn}\end{pmatrix}},

so ist die Gleichung (3.41) offenbar äquivalent mit den Gleichungen

(3.42)

a^{j}=\sum _{i=1}^{j}r_{ij}q^{i},\quad j=1,...,n,

wobei $r_{ij}=0$ für $i>j$ berücksichtigt wurde. Dabei sind die $a^{j}$ $(j=1,...,n)$ linear unabhängig und die $q^{j}$ $(j=1,...,n)$ wegen der Orthogonalität von $Q$ paarweise orthonormal und damit auch linear unabhängig. Die Gleichungen (3.42) implizieren somit für $1\leq s\leq n$

(3.43)

\operatorname {span} \{a^{1},...,a^{s}\}=\operatorname {span} \{q^{1},...,q^{s}\}.

Folglich suchen wir eine orthogonale Basis $\{q^{1},...,q^{n}\}$ des $\mathbb {R} ^{n}$ , für welche die Gleichungen (3.42) erfüllt sind. Wir wollen im Folgenden zeigen, dass man eine solche durch Orthogonalisierung der $a^{j}$ $(j=1,...,n)$ mittels des aus der Linearen Algebra bekannten Gram-Schmidt-Orthonormierungsverfahrens erhält.

Sind $a^{j}$ $(j=1,...,n)$ die gegebenen Spalten von $A$ , welche nach Voraussetzung eines Basis des $\mathbb {R} ^{n}$ bilden, so lautet für diese das Gram-Schmidt-Orthonormierungsverfahren, wobei wir die durch das Verfahren erzeugten Vektoren gleich mit $q^{j}$ bezeichnen:

(3.44)

{\hat {q}}^{j}:=aj-\sum _{i=1}^{j-1}\left[(a^{j})^{T}q^{i}\right]q^{i},\quad q^{j}:={\frac {{\hat {q}}^{j}}{\|{\hat {q}}^{j}\|_{2}}},\quad j=1,2,...,n.

Bekanntlich gilt für die so erzeugten Vektoren $q^{j}$ $(j=1,...,n)$ die Identität (3.43) und sind diese paarweise orthonormal.

Aus den Gleichungen (3.44) leitet man unmittelbar die folgenden Gleichungen her:

a^{j}=\left\|{\hat {q}}^{j}\right\|_{2}q^{j}+\sum _{i=1}^{j-1}\left[(a^{j})^{T}q^{i}\right]q^{i},\quad j=1,2,...,n.

Wie ein Vergleich mit den Gleichungen (3.42) zeigt, erhält man demnach die gewünschte QR-Faktorisierung von $A$ für die Matrix $Q$ mit den durch das Gram-Schmidt-Verfahren erzeugten Vektoren $q^{j}$ $(j=1,...,n)$ und die obere Dreiecksmatrix $R$ , welche folgende Elemente hat:

r_{ij}:=(a^{j})^{T}q^{i},\quad i=1,...,j-1,\quad r_{jj}:=\left\|{\hat {q}}^{j}\right\|_{2}\quad (j=1,...,n).

Der hier beschriebene Orthogonalisierungsprozess ist aber unter Umständen nicht gutartig, etwa für $\|{\hat {q}}^{j}\|_{2}\ll 1$ und damit $r_{jj}\ll 1$ (s. Beispiel 3.4 in Band 1 von Quarteroni/Sacco/Saleri und S. 177 bei Stoer). Mit wachsendem $j$ kann die Orthogonalität immer stärker verloren gehen. Deshalb werden zumeist andere Methoden, wie die im folgenden Abschnitt dargestellte, zur QR-Faktorisierung von $A$ herangezogen.

3.4.4 QS-Zerlegung mittels Householder-Transformationen

Sei nun allgemeiner $A\in \mathbb {R} ^{n\times k}$ mit $1\leq k\leq n$ gegeben. Gegenstand dieses Abschnitts ist die Bestimmung einer Zerlegung der Form $A=QS$ mit einer orthogonalen Matrix $Q\in \mathbb {R} ^{n\times n}$ und einer im Fall $n>k$ nichtquadratischen Matrix $S\in \mathbb {R} ^{n\times k}$ wie in (3.39) mittels sog. Householder-Transformationen. In diesem Zusammenhang zeigen wir zunächst:

Lemma 3.27

Es sei $\omega \in \mathbb {R} ^{n}$ ein Vektor mit $\|\omega \|_{2}=1$ und ${\mathcal {H}}\in \mathbb {R} ^{n\times n}$ die Matrix

(3.45)

{\mathcal {H}}:=I-2\omega \omega ^{T}.

Dann gilt

(3.46) ${\mathcal {H}}^{T}={\mathcal {H}}$ ( ${\mathcal {H}}$ ist symmetrisch),

(3.47) ${\mathcal {H}}^{2}=I$ ( ${\mathcal {H}}$ ist involutorisch),

(3.48) ${\mathcal {H}}^{T}{\mathcal {H}}=I$ ( ${\mathcal {H}}$ ist orthogonal).

Beweis.

Die Beziehungen (3.46) und (3.47) ergeben sich aus

{\mathcal {H}}^{T}=I-2(\omega \omega ^{T})^{T}=I-2\omega \omega ^{T}={\mathcal {H}},

{\mathcal {H}}^{2}=(I-2\omega \omega ^{T})(I-2\omega \omega ^{T})=I-2\omega \omega ^{T}-2\omega \omega ^{T}+4\omega \underbrace {(\omega ^{T}\omega )} _{=1}\omega ^{T}=I.

Die Identität (3.48) folgt aus (3.46) und (3.47).

q.e.d.

Eine Matrix vom Typ ${\mathcal {H}}$ (3.45) für ein $\omega \in \mathbb {R} ^{n}$ mit $\|\omega \|_{2}=1$ nennen wir Householder-Matrix und eine lineare Abbildung

h:\mathbb {R} ^{n}\to \mathbb {R} ^{n},\quad h(x):={\mathcal {H}}x

mit einer Householder-Matrix ${\mathcal {H}}\in \mathbb {R} ^{n\times n}$ bezeichnen wir als Householder-Transformation.

Householder-Matrizen kann man dazu verwenden, einen Block von Komponenten eines gegebenen Vektors $x$ (durch Multiplikation mit einer solchen Matrix von links) zu Null zu setzen. Will man insbesondere alle Komponenten von $x$ außer der $m$ -ten Komponente Null setzen, so muss man den Vektor $\omega$ zur Konstruktion von ${\mathcal {H}}$ so wählen, wie er im folgenden Lemma angegeben wird. Dabei ist wieder mit $e^{m}\in \mathbb {R} ^{n}$ die $m$ -te Spalte von $I\in \mathbb {R} ^{n\times n}$ gemeint.

Lemma 3.28

Gegeben sei $x\in \mathbb {R} ^{n}\setminus \{0\}$ mit $x\notin \operatorname {span} \{e^{m}\}$ . Weiter sei

(3.49)

\omega :={\frac {x+\sigma e^{m}}{\|x+\sigma e^{m}\|_{2}}}

mit

(3.50)

\sigma :=\pm \|x\|_{2}.

Dann hat man $\|\omega \|_{2}=1$ und für ${\mathcal {H}}:=I-2\omega \omega ^{T}$

(3.51)

{\mathcal {H}}x=-\sigma e^{m}.

Beweis.

Wegen $x\notin \operatorname {span} \{e^{m}\}$ ist der Nenner in (3.49) ungleich Null und damit $\omega$ in (3.49) wohldefiniert. Offenbar ist $\|\omega \|_{2}=1$ . Ferner gilt

\|x+\sigma e^{m}\|_{2}^{2}=\|x\|_{2}^{2}+2\sigma (e^{m})^{T}x+\sigma ^{2}=2\|x\|_{2}^{2}+2\sigma (e^{m})^{T}x=2(x+\sigma e^{m})^{T}x

und damit

2\omega ^{T}x={\frac {2(x+\sigma e^{m})^{T}x}{\|x+\sigma e^{m}\|_{2}}}=\|x+\sigma e^{m}\|_{2}.

Zusammen mit (3.49) gelangt man zu der Darstellung

{\mathcal {H}}x=\left(I-2\omega \omega ^{T}\right)x=x-\|x+\sigma e^{m}\|_{2}\omega =x-(x+\sigma e^{m})=-\sigma e^{m}.

Bemerkung 3.29

Zur Vermeidung von Stellenauslöschungen bei der Berechnung von (3.50) ist es offenbar sinnvoll, $\sigma :=\operatorname {sgn} (x_{m})\|x\|_{2}$ zu wählen, wobei $x_{m}$ die $m$ -te Komponente von $x$ bezeichnet.

Wir geben ein Beispiel.

Beispiel 3.30

Sei $x:=(1,1,1,1)^{T}$ und $m:=3$ . Dann errechnen wir $\|x\|_{2}=2$ und setzen wir somit $\sigma :=\operatorname {sgn} (x_{3})\|x\|_{2}=2$ sowie

\omega :={\frac {x+\sigma e^{m}}{\|x+\sigma e^{m}\|_{2}}}={\frac {1}{2{\sqrt {3}}}}(1,1,3,1)^{T}.

Es ergibt sich

2\omega \omega ^{T}={\frac {1}{6}}{\begin{pmatrix}1\\1\\3\\1\end{pmatrix}}(1,1,3,1)={\frac {1}{6}}{\begin{pmatrix}1&1&3&1\\1&1&3&1\\3&3&9&3\\1&1&3&1\end{pmatrix}}

und demnach

{\mathcal {H}}:=I-2\omega \omega ^{T}={\frac {1}{6}}{\begin{pmatrix}5&-1&-3&-1\\-1&5&-3&-1\\-3&-3&-3&-3\\-1&-1&-3&5\end{pmatrix}},\quad {\mathcal {H}}x={\begin{pmatrix}0\\0\\-2\\0\end{pmatrix}}.

Will man für ein $k\geq 1$ die Komponenten $x_{i}$ $(i=1,...,k-1)$ von $x\in \mathbb {R} ^{n}$ unverändert lassen und gleichzeitig $x_{i}=0$ $(i=k+1,...,n)$ erreichen, bekommt man dies, indem man $x$ von links mit der Transformationsmatrix

P_{k}:={\begin{pmatrix}I_{k-1}&0\\0&{\mathcal {H}}_{k}\end{pmatrix}}

multipliziert. Dabei ist $I_{\ell }$ die $(\ell \times \ell )$ -Einheitsmatrix und ${\mathcal {H}}_{k}$ die $(n-k+1)\times (n-k+1)$ -Householder-Matrix der Form

{\mathcal {H}}_{k}:=I_{n-(k-1)}-2\omega ^{k}(\omega ^{k})^{T},\quad \omega ^{k}:={\frac {x^{n-k+1}+\sigma _{k}e^{n-k+1,1}}{\|x^{n-k+1}+\sigma _{k}e^{n-k+1,1}\|_{2}}},\quad \sigma _{k}:=\pm \left\|x^{n-k+1}\right\|_{2},

welche mit dem aus den letzten $n-k+1$ Komponenten von $x$ bestehenden Vektor $x^{n-k+1}\in \mathbb {R} ^{n-k+1}$ und der ersten Spalte $e^{n-k+1,1}\in \mathbb {R} ^{n-k+1}$ der Einheitsmatrix $I_{n-k+1}$ gebildet wird. Denn ist $x^{k-1}\in \mathbb {R} ^{k-1}$ der aus den ersten $k-1$ Komponenten von $x$ bestehende Vektor, so hat man nach Lemma 3.28

P_{k}x={\begin{pmatrix}I_{k}&0\\0&{\mathcal {H}}_{k}\end{pmatrix}}{\begin{pmatrix}x^{k-1}\\x^{n-k+1}\end{pmatrix}}={\begin{pmatrix}x^{k-1}\\{\mathcal {H}}_{k}x^{n-k+1}\end{pmatrix}}={\begin{pmatrix}x^{k-1}\\\sigma _{k}e^{n-k+1,1}\end{pmatrix}}.

Nun ist klar, wie man eine Matrix $A\in \mathbb {R} ^{n\times k}$ in die Form $QS$ mit einer orthogonalen Matrix $Q\in R^{n\times n}$ und eine Matrix $S\in \mathbb {R} ^{n\times k}$ der Gestalt (3.39) zerlegen kann. Ausgehend von $A^{(1)}:=A$ werden sukzessive Matrizen der Form

(3.52)

A^{(j)}:={\begin{pmatrix}a_{11}^{(j)}&...&...&...&...&a_{1k}^{(j)}\\&\ddots &&&&\vdots \\&&a_{j-1,j-1}^{(j)}&...&...&a_{j-1,k}^{(j)}\\&&&a_{jj}^{(j)}&...&a_{jk}^{(j)}\\&&&\vdots &\ddots &\vdots \\&&&a_{nj}^{(j)}&...&a_{nk}^{(j)}\end{pmatrix}}\in \mathbb {R} ^{n\times k},\quad j=2,...,k+1

bestimmt, so dass dann schließlich $S:=A^{(k+1)}$ die gewünschte Gestalt hat. Die Matrizen in (3.52) werden dabei sukzessive durch Transformationen der Form

(3.53)

A^{(j+1)}:={\mathcal {P}}^{(j)}A^{(j)},\quad {\mathcal {P}}^{(j)}:={\begin{pmatrix}I_{j-1}&0\\0&{\mathcal {H}}^{(j)}\end{pmatrix}}

mit Householder-Matrizen

{\mathcal {H}}^{(j)}:=I_{n-(j-1)}-2\omega ^{(j)}(\omega ^{(j)})^{T}

erzeugt, wobei $\omega ^{(j)}\in \mathbb {R} ^{n-(j-1)}$ mit

{\hat {a}}^{(j)}:=(a_{jj}^{(j)},...,a_{nj}^{(j)})^{T}

so zu wählen ist, dass mit einem $\sigma _{j}\in \mathbb {R}$

{\mathcal {H}}^{(j)}{\hat {a}}(j)=-\sigma _{j}e^{n-(j-1),1}

gilt. Im Fall ${\hat {a}}^{(j)}\notin \operatorname {span} \{e^{n-(j-1),1}\}$ erreicht man dies gemäß Lemma 3.28 mit

\omega ^{(j)}:={\frac {{\hat {a}}^{(j)}+\sigma _{j}e^{n-(j-1),1}}{\|{\hat {a}}^{(j)}+\sigma _{j}e^{n-(j-1),1}\|_{2}}},\quad \sigma _{j}:=\operatorname {sgn} ({\hat {a}}_{1}^{(j)})\left\|{\hat {a}}^{(j)}\right\|_{2}.

Im selteneren Fall ${\hat {a}}^{(j)}\in \operatorname {span} \{e^{n-(j-1),1}\}$ kann man offenbar ${\mathcal {H}}^{(j)}:=I_{n-(j-1)}$ bzw. ${\mathcal {P}}^{(j)}:=I$ in (3.53) wählen. Nach Lemma 3.27 sind dann die Matrizen ${\mathcal {H}}^{(1)},...,{\mathcal {H}}^{(k)}$ und damit die Matrizen ${\mathcal {P}}^{(1)},...,{\mathcal {P}}^{(k)}$ orthogonal und symmetrisch, so dass man wegen

A^{(1)}=A,

A^{(2)}={\mathcal {P}}^{(1)}A^{(1)}={\mathcal {P}}^{(1)}A,

A^{(3)}={\mathcal {P}}^{(2)}A^{(2)}={\mathcal {P}}^{(2)}{\mathcal {P}}^{(1)}A,

\vdots

A^{(k+1)}={\mathcal {P}}^{(k)}A^{(k)}={\mathcal {P}}^{(k)}{\mathcal {P}}^{(k-1)}\cdots {\mathcal {P}}^{(1)}A

mit

S:={\mathcal {P}}^{(k)}{\mathcal {P}}^{(k-1)}\cdots {\mathcal {P}}^{(1)}A,\quad Q:={\mathcal {P}}^{(1)}{\mathcal {P}}^{(2)}\cdots {\mathcal {P}}^{(k)}

wegen (3.47) die gewünschte Zerlegung $A=QS$ erhält. Gemäß Lemma 3.26 ist $Q$ tatsächlich eine orthogonale Matrix.

Man beachte, dass man $Q$ für die Lösung des Gleichungssystems $Ax=b$ mit $n=k$ (sowie für die des Ausgleichsproblems in Abschnitt 4.3) gar nicht benötigt, sondern nur den Vektor $Q^{T}b$ . Wegen $\left({\mathcal {P}}^{(j)}\right)^{T}={\mathcal {P}}^{(j)}$ ist dabei

Q^{T}=\left({\mathcal {P}}^{(k)}\right)^{T}\left({\mathcal {P}}^{(k-1)}\right)^{T}\cdots \left({\mathcal {P}}^{(1)}\right)^{T}={\mathcal {P}}^{(k)}{\mathcal {P}}^{(k-1)}\cdots {\mathcal {P}}^{(1)},

so dass man mit $b$ wie mit $A$ verfahren kann:

b^{(1)}:=b,

b^{(2)}:={\mathcal {P}}^{(1)}b^{(1)}={\mathcal {P}}^{(1)}b,

b^{(3)}:={\mathcal {P}}^{(2)}b^{(2)}={\mathcal {P}}^{(2)}{\mathcal {P}}^{(1)}b,

\vdots

b^{(k+1)}:={\mathcal {P}}^{(k)}b^{(k)}{\mathcal {P}}^{(k)}{\mathcal {P}}^{(k-1)}\cdots {\mathcal {P}}^{(1)}b=Q^{T}b.

Man beachte, dass mit

B^{(j)}:={\begin{pmatrix}a_{jj}^{(j)}&...&a_{jk}^{(j)}\\\vdots &\ddots &\vdots \\a_{nj}^{(j)}&...&a_{nk}^{(j)}\end{pmatrix}}

gemäß (3.52) und (3.53)

A^{(j+1)}={\begin{pmatrix}I_{j-1}&0\\0&{\mathcal {H}}^{(j)}\end{pmatrix}}{\begin{pmatrix}\#&*\\0&B^{(j)}\end{pmatrix}}={\begin{pmatrix}\#&*\\0&{\mathcal {H}}^{(j)}B^{(j)}\end{pmatrix}}

gilt, also wie beim Gauß-Algorithmus nur die verbleibende Restmatrix $B^{(j)}$ in $A^{(j)}$ und analog der verbleibende Anteil der rechten Seite zu transformieren ist. Weiter sei gesagt, dass Householder-Transformationsmatrizen niemals explizit durch Matrixmultiplikation gebildet werden sollten. Denn eine Matrixmultiplikation ${\mathcal {H}}B$ mit ${\mathcal {H}}\in \mathbb {R} ^{m\times m}$ und $B\in \mathbb {R} ^{m\times \ell }$ kostet ${\mathcal {O}}(m^{2}\ell )$ Multiplikationen. Die benötigten Matrixmultiplikationen für ${\mathcal {H}}:=I-2\omega \omega ^{T}$ berechne man daher wie folgt:

(3.54)

{\mathcal {H}}B=(I-2\omega \omega ^{T})B=B-\omega v^{T},\quad v^{T}:=2\omega ^{T}B.

Zunächst ermittele man also den Vektor $v$ und dann "aufdatiere" man die Matrix $B$ . Für ein solches Vorgehen benötigt man nur $2m\ell$ Multiplikationen.

Bemerkung 3.31

Will man die Vektoren $\omega ^{(j)}$ aufbewahren, weil man z. B. ein Gleichungssystem $Ax=b$ mit derselben Matrix und verschiedenen rechten Seiten zu lösen hat, so kann man für jedes $j\in \{1,...,k\}$ das Diagonalelement $a_{jj}^{(j+1)}$ zunächst gesondert speichern und anschließend den Vektor $\omega ^{(j)}$ in der $j$ -ten Spalte der Matrix $A^{(j+1)}$ eintragen (s. Beispiel 3.32). Auf diese Weise spart man Speicherplatz, was bei sehr großen Matrizen relevant ist.

Beispiel 3.32

Man berechne eine QS-Zerlegung von $A$ und $Q^{T}b$ für

A:={\begin{pmatrix}1&0\\1&3\\1&4\\1&7\end{pmatrix}},\quad b:={\begin{pmatrix}1\\2\\6\\4\end{pmatrix}}.

Wir setzen $A^{(1)}:=A$ und $b^{(1)}:=b$ und bekommen zunächst

{\hat {a}}^{(1)}:=(1,1,1,1)^{T},\quad \left\|{\hat {a}}^{(1)}\right\|_{2}=2,\quad \sigma _{1}=\operatorname {sgn} (1)\cdot 2=2.

Damit folgt

{\tilde {\omega }}^{(1)}:={\hat {a}}^{(1)}+\sigma _{1}e^{4,1}=(1,1,1,1)^{T}+2(1,0,0,0)^{T}=(3,1,1,1)^{T},

\left\|{\tilde {\omega }}^{(1)}\right\|_{2}=2{\sqrt {3}},

\omega ^{(1)}:={\frac {{\tilde {\omega }}^{(1)}}{\left\|{\tilde {\omega }}^{(1)}\right\|_{2}}}={\frac {1}{2{\sqrt {3}}}}(3,1,1,1)^{T}.

Für ${\mathcal {H}}^{(1)}:=I-2\omega ^{(1)}(\omega ^{(1)})^{T}$ benötigt man nun die Produkte ${\mathcal {H}}^{(1)}A^{(1)}$ und ${\mathcal {H}}^{(1)}b^{(1)}$ , die man analog zu (3.54) berechnet. Es ist

\left(v^{(1)}\right)^{T}:=2\left(\omega ^{(1)}\right)^{T}A^{(1)}={\frac {1}{\sqrt {3}}}{\begin{pmatrix}3&1&1&1\end{pmatrix}}{\begin{pmatrix}1&0\\1&3\\1&4\\1&7\end{pmatrix}}={\frac {1}{\sqrt {3}}}{\begin{pmatrix}6&14\end{pmatrix}},

\left(u^{(1)}\right)^{T}:=2\left(\omega ^{(1)}\right)^{T}b^{(1)}={\frac {1}{\sqrt {3}}}{\begin{pmatrix}3&1&1&1\end{pmatrix}}{\begin{pmatrix}1\\2\\6\\4\end{pmatrix}}={\frac {15}{\sqrt {3}}}

und demnach

\omega ^{(1)}\left(v^{(1)}\right)^{T}={\frac {1}{6}}{\begin{pmatrix}3\\1\\1\\1\end{pmatrix}}{\begin{pmatrix}6&14\end{pmatrix}}={\begin{pmatrix}3&7\\1&7/3\\1&7/3\\1&7/3\end{pmatrix}},

\omega ^{(1)}\left(u^{(1)}\right)^{T}={\frac {15}{6}}{\begin{pmatrix}3\\1\\1\\1\end{pmatrix}}={\begin{pmatrix}15/2\\5/2\\5/2\\5/2\end{pmatrix}}.

Demzufolge ergibt sich

A^{(2)}:={\mathcal {H}}^{(1)}A^{(1)}=A^{(1)}-\omega ^{(1)}\left(v^{(1)}\right)^{T}={\begin{pmatrix}-2&-7\\0&2/3\\0&5/3\\0&14/3\end{pmatrix}},

b^{(2)}:={\mathcal {H}}^{(1)}b^{(1)}=b^{(1)}-\omega ^{(1)}\left(u^{(1)}\right)^{T}={\begin{pmatrix}-13/2\\-1/2\\7/2\\3/2\end{pmatrix}}.

Wollte man den Vektor $\omega ^{(1)}$ aufbewahren, um damit zu einem späteren Zeitpunkt die Zerlegung von $A$ wieder erzeugen zu können, so legt man einen Arbeitsvektor $d\in \mathbb {R} ^{2}$ an, setzt man $d_{1}:=a_{11}^{(2)}=-2$ und trägt man $\omega ^{(1)}$ in die erste Spalte von $A^{(2)}$ ein, so dass sich ergibt:

{\hat {A}}^{(2)}:={\begin{pmatrix}{\frac {1}{2}}{\sqrt {3}}&-7\\{\frac {1}{6}}{\sqrt {3}}&2/3\\{\frac {1}{6}}{\sqrt {3}}&5/3\\{\frac {1}{6}}{\sqrt {3}}&14/3\end{pmatrix}}

In der Praxis kann man $A^{(2)}$ und auch ${\hat {A}}^{(2)}$ wieder $A$ nennen, da man $A$ selbst nicht mehr benötigt.

Nun verfährt man analog mit der Restmatrix bzw. dem Restvektor

{\tilde {A}}^{(2)}:={\begin{pmatrix}2/3\\5/3\\14/3\end{pmatrix}},\quad {\tilde {b}}^{(2)}:={\begin{pmatrix}-1/2\\7/2\\3/2\end{pmatrix}}.

Man bekommt

{\hat {a}}^{(2)}:=\left({\frac {2}{3}},{\frac {5}{3}},{\frac {14}{3}}\right)^{T},\quad \left\|{\hat {a}}^{(2)}\right\|_{2}=5,\quad \sigma _{2}=\operatorname {sgn} \left({\frac {2}{3}}\right)\cdot 5=5,

{\tilde {\omega }}^{(2)}:={\hat {a}}^{(2)}+\sigma _{2}e^{3,1}=\left({\frac {2}{3}},{\frac {5}{3}},{\frac {14}{3}}\right)^{T}+5(1,0,0)^{T}=\left({\frac {17}{3}},{\frac {5}{3}},{\frac {14}{3}}\right)^{T},

\left\|{\tilde {\omega }}^{(2)}\right\|_{2}={\frac {\sqrt {510}}{3}},

\omega ^{(2)}:={\frac {{\tilde {\omega }}^{(2)}}{\left\|{\tilde {\omega }}^{(2)}\right\|_{2}}}={\frac {1}{\sqrt {510}}}(17,5,14)^{T},

\left(v^{(2)}\right)^{T}:=2\left(\omega ^{(2)}\right)^{T}{\tilde {A}}^{(2)}={\frac {2}{\sqrt {510}}}{\begin{pmatrix}17&5&14\end{pmatrix}}{\begin{pmatrix}2/3\\5/3\\14/3\end{pmatrix}}={\frac {510}{3{\sqrt {510}}}},

\left(u^{(2)}\right)^{T}:=2\left(\omega ^{(2)}\right)^{T}{\tilde {b}}^{(1)}={\frac {1}{\sqrt {510}}}{\begin{pmatrix}17&5&14\end{pmatrix}}{\begin{pmatrix}-1\\7\\3\end{pmatrix}}={\frac {60}{\sqrt {510}}},

{\mathcal {H}}^{(2)}{\tilde {A}}^{(2)}={\tilde {A}}^{(2)}-\omega ^{(2)}\left(v^{(2)}\right)^{T}={\begin{pmatrix}2/3\\5/3\\14/3\end{pmatrix}}-{\frac {1}{3}}{\begin{pmatrix}17\\5\\14\end{pmatrix}}={\begin{pmatrix}-5\\0\\0\end{pmatrix}},

{\mathcal {H}}^{(2)}{\tilde {b}}^{(2)}={\tilde {b}}^{(2)}-\omega ^{(2)}\left(u^{(2)}\right)^{T}={\begin{pmatrix}-1/2\\7/2\\3/2\end{pmatrix}}-{\frac {60}{510}}{\begin{pmatrix}17\\5\\14\end{pmatrix}}={\begin{pmatrix}-5/2\\99/34\\-5/34\end{pmatrix}}.

Damit ergibt sich nun

S:=A^{(3)}={\begin{pmatrix}-2&-7\\0&-5\\0&0\\0&0\end{pmatrix}},\quad R:={\begin{pmatrix}-2&-7\\0&-5\end{pmatrix}},\quad Q^{T}b=b^{(3)}:={\begin{pmatrix}-13/2\\-5/2\\99/34\\-5/34\end{pmatrix}}.

Will man die $\omega ^{(j)}$ wiederverwenden, so bilde man mit ${\hat {A}}^{(2)}$ und dem oben definierten $d_{1}$

{\hat {A}}^{(3)}={\begin{pmatrix}{\frac {1}{2}}{\sqrt {3}}&-7\\{\frac {1}{6}}{\sqrt {3}}&{\frac {1}{\sqrt {510}}}17\\{\frac {1}{6}}{\sqrt {3}}&{\frac {1}{\sqrt {510}}}5\\{\frac {1}{6}}{\sqrt {3}}&{\frac {1}{\sqrt {510}}}14\end{pmatrix}},\quad d:={\begin{pmatrix}-2\\-5\end{pmatrix}}.

Aus ${\hat {A}}^{(3)}$ und $d$ könnte man $R$ in offensichtlicher Weise gewinnen.

Man kann sich überlegen, dass eine QS-Zerlegung von $A\in \mathbb {R} ^{n\times k}$ mittels Householder-Transformationen etwa

2\left[{\frac {k^{3}}{3}}+(n-k){\frac {k^{2}}{2}}\right]

Multiplikationen benötigt, also

(3.55)

\sim nk^{2}

für

n\gg k

und

\sim {\frac {2}{3}}n^{3}

für

n\approx k

und damit bei quadratischen Matrizen etwa doppelt so viele wesentliche Rechenoperationen wie der Gauß-Algorithmus. Neben dem Einsatz für die bereits besprochene stabile Lösung von linearen Gleichungssystemen werden wir im nächsten Abschnitt eine weitere wichtige Anwendung einer solchen QS-Zerlegung geben.

This article is issued from Wikiversity. The text is licensed under Creative Commons - Attribution - Sharealike. Additional terms may apply for the media files.