NuMa: Lösen des LGS Ax = b

Leitfaden

Beste Möglichkeit Ax = b zu lösen?

In dieser Reihenfolge überprüfen (Kochrezept):

Überprüfe, ob $LR$-Zerlegung existiert:

• Ist $A$ quadratisch? Falls nicht verwende $QR$-Zerlegung!
• Ist $A$ regulär? Falls nicht verwende $QR$-Zerlegung!
• Sind alle $ det(Ak) \ne 0$? Dabei sind die $Ak$ die Hauptabschnittsmatrizen (s.u.). Falls ja verwende eine $LR$-Zerlegung, falls nein, verwende eine $QR$-Zerlegung

($A$ ist genau dann $LR$-zerlegbar, falls alle Determinanten der Hauptabschnittsmatrizen $A_k$ ungleich 0 sind. Hauptabschnittsmatrizen sind dabei die jeweils von $A_{1,1}$ ausgehend quadratischen Matrizen (also die erste ist nur $A_{1,1}$, $A_2$ ist dann eine $2 \times 2$-Matrix) usw. Daraus folgt, dass $A$ auf jeden Fall schonmal regulär sein muss (sonst kann man auf jeden Fall per Gauss-Algorithmus eine Nullzeile generieren und dann ist besonders offensichtlich, dass die Determinante 0 sein muss.) )

• Falls $A$ symmetrisch und positiv definit (d.h. dann ja auf jeden Fall quadratisch) verwende Cholesky-Zerlegung.
• sonst verwende, falls nicht erforderlich oder anders angegeben Gauss-Elimination ohne Spaltenpivotisierung (numerisch instabil).
• sonst verwende Gauss-Elimination mit Spaltenpivotisierung.

Für jede Matrix $A$ existiert eine $QR$-Zerlegung. Auch hier ein Algorithmus zur Wahl des Verfahrens:

• Falls $A$ dünnbesetzt (hauptsächlich Nullen), verwende Givens-Rotation.
• Sonst verwende Householder-Reflexion.
• Theoretisch ist auch das Gram-Schmidt-Verfahren anwendbar, ist so aber erst mal numerisch instabil. Man kann es jedoch auch numerisch stabilisieren.

LR-Zerlegung

Cholesky-Zerlegung

Wie oben bereits angekündigt, muss die Matrix $A$ symmetrisch, positiv definit und quadratisch sein, um die Cholesky-Zerlegung anwenden zu dürfen.

Sind diese Bedingungen erfüllt, dann lässt sich A folgendermaßen zerlegen: $A = LDL^T$, eingesetzt bedeutet dies $ LDL^Tx = b$, wobei $L$ eine untere linke Dreiecksmatrix und $D$ eine Diagonalmatrix mit $d_{i,k} > 0$, falls $i=k$, $ d_{i,k} = 0 $ sonst, ist (Diagonalatrixeinträge sind alle 0, außer auf Hauptdiagonalen, dort sind sie alle echt größer als 0).

Um nun dieses LGS zu lösen, machen wir einige Definitionen zur schrittweisen Lösung.
$$ L\underbrace{D\overbrace{L^Tx}^{=:v}}_{=:y} = b $$

Für $L$ und $D$ finden wir in unserer Formelsammlung folgende Formeln:
$$ d_{k,k} = a_{k,k} – \sum_{j=1}^{k-1} l^{2}_{j,k} k_{j,j} $$

$$
l_{i,k} =
\frac{
a_{i,k} – \sum_{j=1}^{k-1} l_{i,j} k_{j,j} l_{k,j}
}
{d_{k,k}}
$$

Da wir meistens $3\times 3$ Matrizen lösen müssen, hier noch einmal zur Veranschaulichung, bzw. zur Überprüfung, ob Sie die Formeln richtig angewendet haben (Danke an Sebastian an dieser Stelle):
$$
\begin{aligned}
A = LDL^T & = &
\begin{pmatrix}
1 & 0 & 0\\
l_{2,1} & 1 & 0\\
l_{3,1} & l_{3,2} & 1
\end{pmatrix}
\begin{pmatrix}
d_{1,1} & 0 & 0 \\
0 & d_{2,2} & 0 \\
0 & 0 & d_{3,3}
\end{pmatrix}
\begin{pmatrix}
1 & l_{2,1} & l_{3,1} \\
0 & 1 & l_{3,2} \\
0 & 0 & 1
\end{pmatrix}
\\
&
=
&
\begin{pmatrix}
d_{1,1} & 0 & 0 \\
l_{2,1} d_{1,1} &d_{2,2} & 0 \\
l_{3,1} d_{1,1} &l_{3,2} d_{2,2} & d_{3,3}
\end{pmatrix}
\begin{pmatrix}
1 & l_{2,1} & l_{3,1} \\
0 & 1 & l_{3,2} \\
0 & 0 & 1
\end{pmatrix}
\\
&
=
&
\begin{pmatrix}
d_{1,1} & d_{1,1} \cdot l_{2,1} & d_{1,1} \cdot l_{3,1} \\
d_{1,1} \cdot l_{2,1} & l_{2,1}^2 \cdot d_{1,1} + d_{2,2} & l_{2,1} \cdot d_{1,1} \cdot l_{2,1} + l_{3,2} \cdot d_{2,2} \\
d_{1,1} \cdot l_{3,1} & l_{2,1} \cdot d_{1,1} \cdot l_{2,1} + l_{3,2} \cdot d_{2,2} & l_{3,1}^2 \cdot d_{1,1} + l_{3,2}^2 \cdot d_{2,2} + d_{3,3}
\end{pmatrix}
\end{aligned}
$$

Jetzt lösen wir einfach schrittweise folgende GLS:
$$
Ly = b
$$
$$
Dv = y
$$
$$
L^Tx = v
$$

Beispiel folgt noch nach Fertigstellung der anderen Sektionen dieses Artikels

Gausselimination ohne Spaltenpivotisierung

Sei $A$ eine Matrix mit $i$ Zeilen und $j$ Spalten.Im Prinzip funktioniert das Ganze wie in der Schule. Wir eliminieren spaltenweise alle Einträge unterhalb der Hauptdiagonalen. Betrachten wir den ersten Iterationsschritt. Wir subtrahieren die entsprechenden Vielfachen der ersten Zeile von der $i$-ten Zeile.

Da uns die Nullen unterhalb der Hauptdiagonalen nicht interessieren, notieren wir uns an deren Stellen ebene jene gerade erwähnte Faktoren (Bezogen auf die 1. Iteration:
Sie müssen sich dabei folgende Frage stellen: Was ist der Faktor mit dem ich die erste Zeile mutiplizieren, sodass das Ergebnis verschwindet, wenn ich diese Zeile z.B. von der zweiten abziehe).

Beispiel: Bestimme die $LR$-Zerlegung von $A$, wobei
$$
A =
\begin{pmatrix}
8 & 4 & 12 \\
2 & 4 & 16 \\
4 & 8 & 2
\end{pmatrix}
$$

Außerdem seien die in Boxen dargestellten Einträge, die o.a. Faktoren, aus denen sich nachher unsere $L$-Matrix zusammensetzt.

Hinweis: wir könnten nebenbei auch das LGS $Ax = b$ lösen, indem wir $b$ als zusätzliche Spalte in $A$ eintragen.

Nach der ersten Iteration:
$$
A_{{LR}_1} =
\begin{pmatrix}
8 & 4 & 12 \\
\boxed{\frac{1}{4}} & 3 & 13 \\
\boxed{\frac{1}{2}} & 6 & -4
\end{pmatrix}
$$
$
l_{2,1} = \frac{a_{2,1}}{a_{1,1}} = \frac{1}{4} \\
l_{3,1} = \frac{a_{2,1}}{a_{1,1}} = \frac{1}{2}
$

Nach der zweiten Iteration:
$$
A_{{LR}_2} =
\begin{pmatrix}
8 & 4 & 12 \\
\boxed{\frac{1}{4}} & 3 & 13 \\
\boxed{\frac{1}{2}} & \boxed{2} & -30
\end{pmatrix}
$$
$
l_{3,2} = \frac{a_{3,2}}{a_{2,2}} = 2
$

Jetzt erhalten wir unsere Matrizen $L$ und $R$ als:

$$
L =
\\
\begin{pmatrix}
1 & 0 & 0 \\
\frac{1}{4} & 1 & 0 \\
\frac{1}{2} & 2 & 1
\end{pmatrix}
$$

$$
R =
\begin{pmatrix}
8 & 4 & 12 \\
0 & 3 & 13 \\
0 & 0 & -30
\end{pmatrix}
$$

Was, wann, warum Spaltenpivotisierung?

Man bezeichnet dasjenige Diagonalelement als Pivotelement, das im aktuellen Iterationsschritt betrachtet wird. Spaltenpivotisierung bedeutet nun, dass Zeilen der Matrix derart
vertauscht werden, dass das betragsmäßig größte Element Pivotelement ist. Keine Sorge, weiter unten ist ein Beispiel.

Eine solche Vertauschung nennt man Permutation. Um nun die Vertauschung festzuhalten, muss von links an die Matrix noch die sogenannte Permutationsmatrix multipliziert werden.
Die Permutationsmatrix $P$ ist fast die Einheitsmatrix $I$, außer dass zwei Zeilen entsprechend der im Gaussalgorithmus vertauschten Zeilen vertauscht werden.

Beispiel: werden im ersten Iterationsschritt die ersten zwei Zeilen vertauscht ($P_1$), dann im zweiten Iterationsschritt die erste und dritte Zeile vertauscht ($P_2$), dann berechnet sich $P$ folgendermaßen:
$$
P = P_2 \cdot P_1 =
\begin{pmatrix}
0 & 0 & 1 \\
0 & 1 & 0 \\
1 & 0 & 0
\end{pmatrix}
\begin{pmatrix}
0 & 1 & 0 \\
1 & 0 & 0 \\
0 & 0 & 1
\end{pmatrix}
\\
=
\begin{pmatrix}
0 & 0 & 1 \\
1 & 0 & 0 \\
0 & 1 & 0
\end{pmatrix}
$$

Für die Spaltenpivotisierung gibt es zwei Gründe:

• Erforderlich, wenn Pivotelement 0,
  da sonst Gaussalgorithmus undurchführbar. $(I)$
• Stabilisierung des Gaussalgorithmus. $(II)$

Kurz wir machen immer dann eine Spaltenpivotisierung, wenn nötig ($I$) oder gefordert ($II$). Ein Beispiel folgt im nächsten Abschnitt.

Gausselimination mit Spaltenpivotisierung

Dieses Beispiel stammt aus Dahmen/Reusken S.72. Lösen wir ein Gleichungssystem der Form $Ax=b$ mit 4-stelliger Gleitpunktarithmetik (GPA), wobei
$$
A =
\begin{pmatrix}
0.00031 & 1 \\
1 & 1
\end{pmatrix}
\quad
\\
b=
\begin{pmatrix}
-3 \\
-7 \\
\end{pmatrix}
$$

Dann vergleichen wir das Ergebnnis mit einer exakten Rechnung. Ohne Pivotisierung erhalten wir
$$
Rc =
\begin{pmatrix}
0.00031 &
1 &
-3
\\
\boxed{\frac{1}{0.00031}}&
1- \frac{1}{0.00031} &
-7 – \frac{-3}{0.00031}
\end{pmatrix}
$$
Bei vierstelliger GPA kommen wir da auf:
$
\tilde{x}_1 = -6.452 \quad
\tilde{x}_2 = -2.998
$

Bei genauer Rechnung kommen wir da auf:
$
x_1=-4.00124… \quad
x_2=-2.998759…
$

$\tilde{x}_1$ ist also auf keiner Stelle richtig. Dabei beträgt die Kondition der Matrix gerade einmal $\kappa_\infty (A)$ = 4. Dies bedeutet, dass der Gaußalgorithmus ohne Pivotisierung extrem instabil ist. Mit Pivotisierung können wir ihn dann stabilisieren und erhalten bei vierstelliger GPA:
$$
\tilde{x}_1 = -4.001 \quad
\tilde{x}_2 = -2.999
$$

Dies ist völlig akzeptabel.

QR-Zerlegung

Eigenschaften der QR-Zerlegung

$Q$ ist nicht irgendeine beliebige Matrix, sondern sie ist orthogonal.Die wichtigste Eigenschaft einer orthogonalen Matrix, die wir hier ausnutzen ist: $Q^T = Q^{-1}$. Wie jetzt vielleicht auch leicht zu sehen ist:
$$
Ax = b
\overset{\exists \text{QR-Zerlegung}}{\Leftrightarrow}
QRx = b
\Leftrightarrow
Rx = Q^{-1}b
\Leftrightarrow
Rx = Q^{T}b
$$
Rückwärtseinsetzen liefert uns dann in $O(n^2)$ $x$.

Givens-Rotation

Bevor ich das Verfahren vorstelle erstmal die Vor- und Nachteile dieser Methode vorausgeschickt:

Vorteile:

• Auf beliebige $m\times n$ Matrizen anwendbar.
• Sehr stabil, Pivotisierung weder nötig noch sinnvoll.
• Flexibel auf dünnbesetzte Matrix anwendbar, da
  Einträge einzeln eliminiert werden.

Nachteile:

• Aufwendig im vergleich zum Gaußalgorithmus:
  $\frac{4}{3} n^3$ Operationen falls $ m \approx n $
  und $2mn^2$. Durch schnelle Givens-Rotation
  kann der Aufwand halbiert werden.
  Bei dünnbesetzten Matrizen ist der Aufwand jedoch
  wesentlich geringer.

Zum Verfahren:
wie oben bereits erwähnt, kann in jedem Iterationsschritt immer nur eine Stelle eliminiert werden, anstatt einer kompletten Spalte wie beim Gaussverfahren. Die Eliminiation dieser Stellen erfolgt durch linksseitige Multiplikation mit Givensmatrizen $G_{i,j}$. Anders als im Buch Dahmen/Reusken bezeichnet $i$ die entsprechende Zeile von $A$, $j$ die Spalte. Wir werden nachher sehen, wie sich die Indizees auf die $G$-Matrizen auswirken. Bezeichne die Notation $G_k$ diejenige Matrix, die im $k$-ten Schritt links an $A$ multipliziert wird und $A_k$ diejenige Matrix für die gilt (Induktion, wobei $A_0 = A$): $G_k \cdot A_{k-1} = A_k$ Dann ergeben sich:
$$
Q = G_n^T \cdot G_{n-1}^T \cdot \dots \cdot G_2^T \cdot G_1^T \\
R = A_n
$$
Keine Sorge, das wird alles an einem Beispiel noch gezeigt.

Zusätzlich muss man darauf achten, eine geeignete Reihenfolge bei der Elimination zu der Einträge wählen. Sie sind auf der sicheren Seite, wenn Sie entweder spaltenweise
oder zeilenweise unterhalb der entsprechenden Diagonalelementen eliminieren.

(Der Grund dafür ist, dass die Multiplikation mit einer $G_{i,j}$-Matrix jeweils (nur) die $i$-te und $j$-te Zeile von $A$ beeinflusst und man vermeiden möchte, dass bereits
generierte Nulleinträge im nächsten Schritt wieder verloren gehen. Probieren sie ruhig mal eine falsche Reihenfolge aus, um zu sehen was passiert.)

Nun wie sind $G$-Matrizen überhaupt aufgebaut? In einer solchen Form:
$$
G_{i,j} =
\begin{pmatrix}
1 & 0 & \cdots
\\
0 & 1 & 0 & \cdots
\\
\vdots & 0 & \ddots
\\
& & & c & \cdots & -s
\\
& & & \vdots & & \vdots
\\
& & & s & \cdots & c
\\
& & & & & & \ddots
\\
& & & & & & & 1 \\
& & & & & & & & 1 \\
\end{pmatrix}
$$

Ohne lange Erklärungen und ohne geometrische Interpretationen hier das Kochrezept wie man fortfährt:

$$
c = \frac{a_{jj}}{\sqrt{a_{jj}^2 + a_{ij}^2}} \\
s = \frac{-a_{ij}}{\sqrt{a_{jj}^2 + a_{ij}^2}}
$$

Nun wie ist die Position der der $c$, $s$, $-s$ in $G$? Nun da $G_k$ nur die $i$-te und $j$-te Zeile von $A_{k-1}$ beeinflussen kann, müssen die $c$, $s$, $-s$ sich in den
$i$-ten und $j$-ten Spalten und Zeilen von $G_k$ befinden.

Beispiel für $G_{2,4}$ oder für $G_{4,2}$, für den fall $G \in \mathbb{R}^{4 \times 4}$:
$$
G =
\begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & c & 0 & -s \\
0 & 0 & 1 & 0 \\
0 & s & 0 & c
\end{pmatrix}
$$

Nun das ganze Verfahren mal an einem Beispiel durchexerziert:

Sei
$$
A =
\begin{pmatrix}
3 & 4 \\
0 & 2 \\
0 & 1 \\
0 & 1
\end{pmatrix}
$$
wir sehen: $A$ ist dünnbesetzt $\Rightarrow$ Givens-Rotation sinnvoll!

Wir wollen als erstes das Element $a_{3,2}$ eliminieren.
Daher:
$c_{3,2} = \frac{a_{jj}}{\sqrt{a_{jj}^2+a_{ij}^2}} = \frac{2}{\sqrt{5}}$
und
$s_{3,2} = \frac{-a_{ij}}{\sqrt{a_{jj}^2+a_{ij}^2}} = \frac{-1}{\sqrt{5}}$

Für $G_1$ ergibt sich:
$$
G_1 =
\begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & \frac{2}{\sqrt{5}} & \frac{1}{\sqrt{5}} & 0 \\
0 & \frac{-1}{\sqrt{5}} & \frac{2}{\sqrt{5}} & 0 \\
0 & 0 & 0 & 1
\end{pmatrix}
$$

Dann ist:
$$
A_1 = G_1 \cdot A =
\begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & \frac{2}{\sqrt{5}} & \frac{1}{\sqrt{5}} & 0 \\
0 & \frac{-1}{\sqrt{5}} & \frac{2}{\sqrt{5}} & 0 \\
0 & 0 & 0 & 1
\end{pmatrix}
\cdot
\begin{pmatrix}
3 & 4 \\
0 & 2 \\
0 & 1 \\
0 & 1
\end{pmatrix}
=
\begin{pmatrix}
3 & 4 \\
0 & \sqrt{5} \\
0 & 0 \\
0 & 1
\end{pmatrix}
$$

Wir wollen als nächtes das Element $a_{4,2}$ eliminieren. Daher:
$c_{4,2} = \frac{a_{jj}}{\sqrt{a_{jj}^2+a_{ij}^2}} = \frac{\sqrt{5}}{\sqrt{6}}$
und
$s_{4,2} = \frac{-a_{ij}}{\sqrt{a_{jj}^2+a_{ij}^2}} = \frac{-1}{\sqrt{6}}$

Für $G_2$ ergibt sich:
$$
G_2 =
\begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & \frac{\sqrt{5}}{\sqrt{6}} & 0 & \frac{1}{\sqrt{6}} \\
0 & 0 & 1 & 0 \\
0 & \frac{-1}{\sqrt{6}} & 0 & \frac{\sqrt{5}}{\sqrt{6}}
\end{pmatrix}
$$

Dann ist:
$$
R = A_2 = G_2 \cdot G_1 \cdot A = G_2 \cdot A_1 =
\\
\begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & \frac{\sqrt{5}}{\sqrt{6}} & 0 & \frac{1}{\sqrt{6}} \\
0 & 0 & 1 & 0 \\
0 & \frac{-1}{\sqrt{6}} & 0 & \frac{\sqrt{5}}{\sqrt{6}}
\end{pmatrix}
\cdot
\begin{pmatrix}
3 & 4 \\
0 & \sqrt{5} \\
0 & 0 \\
0 & 1
\end{pmatrix}
\\
=
\begin{pmatrix}
3 & 4 \\
0 & \sqrt{6} \\
0 & 0 \\
0 & 0
\end{pmatrix}
$$

Nun brauchen wir noch $Q$:
$$
Q = G_2^T \cdot G_1^T =
\begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & \frac{\sqrt{5}}{\sqrt{6}} & 0 & \frac{-1}{\sqrt{6}} \\
0 & 0 & 1 & 0 \\
0 & \frac{1}{\sqrt{6}} & 0 & \frac{\sqrt{5}}{\sqrt{6}}
\end{pmatrix}
\cdot
\begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & \frac{2}{\sqrt{5}} & \frac{-1}{\sqrt{5}} & 0 \\
0 & \frac{1}{\sqrt{5}} & \frac{2}{\sqrt{5}} & 0 \\
0 & 0 & 0 & 1
\end{pmatrix}
\\
=
\begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & \frac{2}{\sqrt{6}} & \frac{-1}{\sqrt{6}} & \frac{-1}{\sqrt{6}} \\
0 & \frac{1}{\sqrt{5}} & \frac{2}{\sqrt{5}} & 0 \\
0 & \frac{2}{\sqrt{30}} & \frac{-1}{\sqrt{30}} & \frac{\sqrt{5}}{\sqrt{6}}
\end{pmatrix}
$$

$Q$ ist auch orthogonal, wie leicht nachzuprüfen ist: die $\norm{\cdot}_2$ aller Zeilen und Spalten ist 1.

Householder-Transformation