\section{Der Cholesky--Algorithmus, Bandmatrizen}\label{KAPITEL-02.2}
%
Ein guter numerischer Algorithmus sollte spezielle Eigenschaften 
des mathematischen Problems ber\"ucksichtigen. Der Gau{\ss}--Algorithmus ist 
auf allgemeine lineare Gleichungssysteme zugeschnitten. 
Nun sind aber viele der Anwendungsproblem
aus Ingenierwissenschaft oder Physik, Chemie so, da{\ss} die Gleichungssysteme
eine spezielle Struktur haben, wie zum Beispiel symmetrische, Hermitesche
Matrizen, Bandstruktur etc. Wir betrachten nun zuerst Hermitesch, postitiv
definite Systeme
%
\eq{2.2.1.}
 Ax = b \; \mbox { mit } \; A = A^*, \; \mbox { d.h. } 
 a_{ij} = \overline{a}_{ij},
\en
%
\[ 
  \mbox { und } x^* Ax > 0, \; \forall x \in \C^n\setminus\{0\}
\]
%
und untersuchen, wie wir den Gau{\ss}--Algorithmus verbessern k\"onnen.
%
\begin{theorem}[Cholesky Zerlegung]\label{2.2.2}\quad\\
Sei $A\in \K^{n,n}$ ($\K = \C$ oder $\K = \R$), Hermitesch positiv definit.
Dann existiert eine untere $\bigtriangleup$--Matrix $G\in\K^{n,n}$ mit
positiven Diagonalelementen, so da{\ss}
\[
 A = GG^*.
\]
%
\end{theorem}
%
\proof
Da $A$ positiv definit ist, sind alle Hauptabschnittsmatrizen
\[
 A(1:k,1:k), \;  k = 1, \dots, n
\]
%
positiv definit, haben also Determinante $\neq 0.$
Also existiert eine eindeutige $LR$--Zerlegung $A=LR$ mit $L$ untere
$\bigtriangleup$--Matrix mit 1--Diagonale, und $R$ hat Diagonalelemente
$>0$ (gleich den Pivots). 
Also kann man $R$ schreiben als $D\tilde{R}$,
wobei $\tilde{R}$ 1--Diagonale hat. Sei $D=\dgl(d_1, \dots, d_n)$. 
Die Pivots sind gerade die Determinanten der Hauptabschnittsmatrizen, 
also reell positiv (alle Eigenwerte sind reell, positiv), also ist 
$d_i>0, \;\forall i=1,\dots,n.$ Wir haben 
$A=LD \tilde{R} = LD^{\frac{1}{2}}D^{\frac{1}{2}} \tilde{R}$. 
$D^{\frac{1}{2}} = \dgl (d_1^{\frac{1}{2}},\dots, d_n^{\frac{1}{2}})$.\\
Also folgt
\[
  D^{-\frac{1}{2}} L^{-1} A L^{-*} D^{-\frac{1}{2}} =
  D^{\frac{1}{2}} \tilde{R} L^{-*} D^{-\frac{1}{2}}.
\]
%
Links steht eine Hermitesche Matrix, rechts eine obere $\bigtriangleup$--Matrix
mit 1--Diagonale. Also mu{\ss} rechts die Einheitsmatrix stehen. Also gilt
\[
 D^{\frac{1}{2}} \tilde{R} = D^{\frac{1}{2}} L^* \Longleftrightarrow
 \tilde{R} = L^*.
\]
%
Setze $G:= LD^{\frac{1}{2}}.$
\eproof
%
Damit k\"onnen wir den Gau{\ss}--Algorithmus verbessern und erhalten den 
Cholesky--Algorithmus.
%
\begin{algorithmus}[Cholesky Zerlegung GAXPY--Version]\label{2.2.3}
\begin{tabular}{llll}
\underline{Input:} & \multicolumn{3}{l}{Hermitesch positiv definite Matrix 
			$A \in \K^{n,n}$.}\\
\underline{Output:}& \multicolumn{3}{p{13.5cm}}{ Untere 
			$\bigtriangleup$--Matrix $G \in \K^{n,n}$ mit
			positiver Diagonale, so da{\ss} $A = GG^*.$ F\"ur 
			$i\le j$ \"uberschreibt $G(i,j)$ den Wert von $A(i,j)$.
		     }\\[4ex]
&FOR&\multicolumn{2}{l}{$j = 1:n$}\\
&   &IF & $j > 1$\\
&   &   & $A(j:n,j) = A(j:n,j)-A(j:n,1:j-1)*A(j,1:j-1)^*$;\\
&   &END&\\ 
&   &\multicolumn{2}{l}{$A(j:n,j) = A(j:n,j) / \sqrt{A(j,j)}$;}\\
&END&   & 
\end{tabular}
\end{algorithmus}
\underline{Kosten:}  $n^3/3$ \flops.\\
\underline{Fehleranalyse:} Analog wie bei Gau{\ss}--Zerlegung.

Wie im Fall der $LR$--Zerlegung kann man auch die Cholesky--Zerlegung anders
aufbauen:
%
\[
 A = \left[ \begin{array}{cc}
     \alpha & v^*\\
     v & B \end{array}\right] = \left[ \begin{array}{cc}
     \beta & 0\\
     v/\beta & I_{n-1} \end{array}\right] 
     \left[ \begin{array}{cc}
     1 & 0\\
     0 & B-vv^* / \alpha \end{array}\right]
     \left[ \begin{array}{cc}
     \beta & v^*/\beta\\
     0 & I_{n-1} \end{array}\right]
\]
%
$\beta = \sqrt{\alpha}$. Wenn die Cholesky Zerlegung von
$B-vv^*/\alpha = G_1G_1^*$ gegeben ist, dann erh\"alt man die Zerlegung,
$A = GG^*$ mit
\[
 G = \left[ \begin{array}{cc}
     \beta & 0\\
     v/\beta & G_1 \end{array}\right].
\]
%
%
\begin{algorithmus}[Cholesky--Zerlegung: \"Au{\ss}ere Produkte Version]\label{2.2.4}
\begin{tabular}{llll}
\underline{Input:} & \multicolumn{3}{l}{Hermitesch positiv definite Matrix 
		     $A \in \K^{n,n}$.}\\
\underline{Output:}& \multicolumn{3}{p{13.5cm}}{
                	Untere $\bigtriangleup$--Matrix $G \in \K^{n,n}$ mit 
			positiver Diagonale, so da{\ss} $A = GG^*$. 
			F\"ur $i \ge j$ \"uberschreibt $G(i,j)$
			jeweils $A(i,j)$.
		     }\\[4ex]
&FOR&\multicolumn{2}{l}{$k = 1:n$}\\
&   &\multicolumn{2}{l}{$A(k,k) = \sqrt{A(k,k)}$;}\\
&   &\multicolumn{2}{l}{$A(k+1:n,k) = A(k+1:n,k)/A(k,k)$;}\\
&   &FOR& $j=k+1:n$\\
&   &   & $A(j:n,j) = A(j:n,j) - A(j:n,k)*A(j,k)$;\\
&   &END&\\ 
&END&   & 
\end{tabular}
\end{algorithmus}
\underline{Kosten:}  $n^3/3$ \flops.\\
\underline{Fehleranalyse} Analog wie bei Gau{\ss}--Zerlegung.

Beachte: Die $j$--Schleife berechnet den unteren $\bigtriangleup$--Anteil
von
\[
 A(k+1:n, k+1:n) = A(k+1:n, k+1:n)-A(k+1:n,k) A(k+1:n,k)^*
\]
%
\underline{Pivotisierung:} Um Symmetrie zu erhalten, sollte man nur auf der 
Diagonalen pivotisieren. Also Permutationen $PAP^*$ verwenden, die die 
Diagonalelemente vertauschen. 

\subsection{Verallgemeinerung f\"ur nicht positiv definite Systeme.}
$LDL^*$--Zerlegung $D$ blockdiagonal mit $1 \times 1$ oder $2 \times 2$
Bl\"ocken. Hier \underline{mu{\ss}} man i.a. permutieren, d.h. man erzeugt 
$P^*AP=LDL^*$.\\
$\longrightarrow$ \"Ubung.

L\"osung des Gleichungssystems dann wie bei Gau{\ss} mit Vorw\"arts-- und 
R\"uckw\"artseinsetzen.

\subsection{Bandsysteme} 
Eine weitere h\"aufig vorkommende spezielle Struktur ist die Bandstruktur. Eine
Matrix hei{\ss}t \underline{$p-q$ Bandmatrix}
mit \underline{unterer Bandbreite $p$} und \underline{oberer Bandbreite $q$}
falls $a_{ij} = 0$ f\"ur $i>j+p$ und $j>i+q$.
%
\begin{example}\label{2.2.5}\quad\\
$A = 
\left[ 
 \begin{array}{ccccc}
  {*}& * & * & 0 & 0\\
  {*}& * & * & * & 0\\
   0 & * & * & * & *\\
   0 & 0 & * & * & *\\
   0 & 0 & 0 & * & *\\
   0 & 0 & 0 & 0 & 0 
 \end{array}
\right]$ 
\begin{tabular}{l}
{ untere Bandbreite } 1\\
{ obere Bandbreite } 2\\
{ 1--2 Bandmatrix. } 
\end{tabular}
\end{example}
%
\begin{theorem}\label{2.2.6}
$A\in \K^{n,n}$ habe eine $LR$--Zerlegung $A=LR$. 
Falls $A$ eine $p-q$--Bandmatrix ist, so hat $L$ untere Bandbreite $p$ und 
$R$ obere Bandbreite $q$.
\end{theorem}
%
\proof
Mit vollst\"andiger Induktion:
\[
 A = \left[ \begin{array}{cc}
      \alpha & w^*\\
      v & B
     \end{array}\right] =
%
  \left[ \begin{array}{cc}
     1 & 0\\
     v/\alpha & I_{n-1}
     \end{array}\right]
%
  \left[ \begin{array}{cc}
     1 & 0\\
     0 & B-vw^*/\alpha
     \end{array}\right]
%
  \left[ \begin{array}{cc}
     \alpha & w^*\\
     0 & I_{n-1}
     \end{array}\right]
\]
%
$B-vw^*/\alpha$ ist $p-q$ Bandmatrix da nur die ersten $q$ Komponenten von $w$
und die ersten $p$ Komponenten von $v$ ungleich $0$ sind. Sei $L_1R_1$ die
 $LR$--Zerlegung von $B-vw^*/\alpha.$ Nach I.V. und
der Sparsit\"at von $v,w$ folgt da{\ss}
\[
 L = \left[ \begin{array}{cc}
     1 & 0\\
     v/\alpha & L_1 \end{array}\right] \quad \mbox { und } \quad
 R = \left[ \begin{array}{cc}
     \alpha & w^*\\
     0 & R_1 \end{array}\right]
\]
%
die gew\"unschten Brandbreiten haben und es gilt $A = LR.$
\eproof
%
Es folgt also, da{\ss} man die Bandstruktur optimal ausnutzen kann 
vorausgesetzt, die $LR$--Zerlegung existiert!\reversemarginpar\marginpar{\huge \hfill !}

%
\begin{algorithmus}[Band Gau{\ss}--Elimination, \"au{\ss}ere
Produkte Version]\label{2.2.7} 
\begin{tabular}{lllll}
\underline{Input:} & \multicolumn{4}{l}{$A\in \K^{n,n} p-q$--Bandmatrix, 
			welches eine $LR$--Zerlegung besitzt.
}\\
\underline{Output:}& \multicolumn{4}{p{13.5cm}}{
			Zerlegung $A=LR$, $A(i,j)$ \"uberschrieben durch
			$L(i,j)$ f\"ur $i>j$ und $R(i,j)$ sonst.
		     }\\[1.5ex]
&FOR&\multicolumn{3}{l}{$k = 1:n-1$}\\
&   &FOR&\multicolumn{2}{l}{$i=k+1:\min (k+p,n)$}\\
&   &   &\multicolumn{2}{l}{$A(i,k) = A(i,k)/A(k,k)$;}\\
&   &END&  & \\ 
&   &FOR&\multicolumn{2}{l}{$j=k+1:\min (k+q,n)$}\\
&   &   &FOR&$i=k+1:\min (k+p,n)$\\
&   &   &   &$A(i,j) = A(i,j) - A(i,k) *A(k,j)$;\\
&   &   &END&\\
&   &END&   &\\
&END&   &   &
\end{tabular}
\end{algorithmus}
\underline{Kosten:} f\"ur $n \gg p,q$:  $\approx 2npq$ \flops\\
\underline{Fehleranalyse:} wie normaler Gau{\ss}.\\
 Cholesky f\"ur Band analog.\\
$\longrightarrow$ \"Ubung .

Spezielle Speichertechniken m\"oglich. Dann kann man auch den 
$\bigtriangleup$--L\"oser verbessern.
%
\begin{algorithmus}[Band--Vorw\"arts--Aufl\"osen]\label{2.2.8}\quad\\
\begin{tabular}{llll}
\underline{Input:} & \multicolumn{3}{l}{$L \in \K^{n,n}$ untere 
			$\bigtriangleup$--Matrix mit 1--Diagonale und unterer 
			Bandbreite $p, b \in \K^n$.}\\
\underline{Output:}& \multicolumn{3}{p{13.5cm}}{
			$b$ \"uberschrieben mit der L\"osung von $Lx=b$.
		     }\\[1.5ex]
&FOR&\multicolumn{2}{l}{$j = 1:n$}\\
&   &FOR&$i=j+1:\min (j+p,n)$\\
&   &   &$b(i) = b(i) - L(i,j)*b(j)$;\\
&   &END& \\ 
&END&   &  
\end{tabular}
\end{algorithmus}
\underline{Kosten:} f\"ur $n \gg p$: $\approx 2np$ \flops.

%
\begin{algorithmus}[Band--R\"uckw\"arts--Aufl\"osen]\label{2.2.9} \quad\\
\begin{tabular}{llll}
\underline{Input:} & \multicolumn{3}{l}{$R \in \K^{n,n}$ obere 
			$\bigtriangleup$--Matrix
			mit oberer Bandbreite $q, b \in \K^n.$}\\
\underline{Output:}& \multicolumn{3}{p{13.5cm}}{
			$b$ \"uberschrieben mit der L\"osung von $Rx=b$.
		     }\\[1.5ex]
&FOR&\multicolumn{2}{l}{$j = n:-1:1$}\\
&   &\multicolumn{2}{l}{$b(j)=b(j)/u(j,j)$;}\\
&   &FOR&$i = \max (1,j-q):j-1$\\
&   &   &$b(i) = b(i) - L(i,j)*b(j)$;\\
&   &END& \\ 
&END&   &  
\end{tabular}
\end{algorithmus}






Tridiagonalmatrizen werden als 3 Vektoren gespeichert, oder bei Hermitesch positiv definiten
Matrizen als 2 Vektoren.
%
\begin{algorithmus}[Symmetrisch, positiv definiter tridiagonal L\"oser]\label{2.2.10}
\begin{tabular}{lll}
\underline{Input:} & \multicolumn{2}{p{13.5cm}}{$A\in\R^{n,n}$ symmetrisch, 
			positiv definit, tridiagonal, $b \in \R^n$. $A$ 
			gespeichert als $d(1:n), e(1:n-1)$.}\\
\underline{Output:}& \multicolumn{2}{p{13.5cm}}{
			L\"osung von $Ax = b$ \"uberschrieben auf $b$.
		     }\\[1.5ex]
&FOR&{$k = 2:n$}\\
&   & $t = e(k-1)$;\\
&   & $e(k-1) = t/d(k-1)$;\\
&   & $d(k) = d(k) - t*e(k-1)$;\\
&END& \\
&FOR&{$k = 2:n$}\\
&   & $b(k)=b(k) - e(k-1)*b(k-1)$;\\
&END& \\
&\multicolumn{2}{l}{$b(n) = b(n) /d(n)$;}\\
&FOR&{$k = n-1:-1:1$}\\
&   & $b(k) = b(k) /d(k) - e(k)*b(k+1)$;\\
&END&     
\end{tabular}
\end{algorithmus}
\underline{Kosten:} $8n$ \flops.\\
$\longrightarrow$ \"Ubung \ref{numI/uebung026}
%
