\section{Iterative Verfahren}\label{KAPITEL-02.5}
F\"ur viele der Gleichungssysteme, die in der Praxis auftauchen, gilt
$n \gg 10000$ und typischerweise $a_{ij} = 0$ f\"ur $>90\%$ der
Eintr\"age. Wenn man derartige Probleme mit Gau{\ss} oder $QR$ l\"ost, 
kann es passieren, das in den Faktoren $L, Q, R$ alle Eintr\"age $\neq 0$ 
sind.
Extrembeispiel Gau{\ss} f\"ur die Matrix
%
\[
 \left[ \begin{array}{cccc}
  {*}   & * &\cdots& *\\     
  {*}   & * &      & 0\\
  \vdots&   &\ddots&  \\
  {*}   & 0 &      & * \end{array}\right].
\]
%
Die Matrix l\"auft total voll. Es gibt sparse Varianten von 
Gau{\ss}/Cholesky/$QR$.
Dies ist Thema einer Spezialvorlesung.
Eine wichtige Idee, die insbesondere im Zusammenhang mit Vektorrechnern und
Parallelrechnern von gr\"o{\ss}ter Bedeutung ist, 
ist die Konstruktion von Verfahren,
die die Matrix an sich nicht ver\"andern, und wenn m\"oglich auf 
Matrix * Matrix oder Matrix * Vektor--Operationen aufgebaut sind. 
Das f\"uhrt auf iterative Verfahren. Im wesentlichen betrachten wir  hier
2 Ans\"atze.

\subsection{Splitting--Verfahren} 
Zerlege $A$ additiv 
%
\eq{2.5.1}
 A = M - N,
\en
%
wobei $M^{-1}$ eine Approximation an $A^{-1}$ ist, die leicht invertierbar ist, 
z.B. diagonal, blockdiagonal, $\bigtriangleup$--Matrix,\dots ~
Ersetze $Ax = b$ durch
%
\eq{2.5.2}
 Mx = Nx + b \Longleftrightarrow x = M^{-1}Nx + M^{-1}b.
\en
%
Dies ist eine Fixpunktgleichung.

\underline{Kanonisches Fixpunktverfahren:} 
%
\eq{2.5.3}
 x^{(k+1)} = M^{-1} Nx^{(k)} + M^{-1}b, \mbox { mit gegebenen } x^{(0)}.
\en
%
\begin{theorem}\label{2.5.4}
Das Fixpunktverfahren (\ref{2.5.3}) konvergiert f\"ur alle Startwerte $x^{(0)}$
gegen die L\"osung von $Ax=b$ genau dann, wenn der Spektralradius 
$\rho (M^{-1}N) < 1$
 ist, wobei $\rho(T)=\max\{|\lambda|: \;\lambda \mbox{ Eigenwert von } T \}$.
\end{theorem}
%
\proof
Banachscher Fixpunktsatz oder:\\
Sei $e^{(k)} = x- x^{(k)}$ der Fehler in der $k$--ten Iterierten, so gilt:
%
\begin{eqnarray*}
 e^{(k+1)} & = & x - x^{(k+1)}\\
           & = & (M^{-1}Nx + M^{-1}b) - (M^{-1}Nx^{(k)} + M^{-1}b)\\
           & = & M^{-1}N (x-x^{(k)})\\
	   & = & M^{-1}N e^{(k)}\\
           & = & (M^{-1} N)^k e^{(0)}.
\end{eqnarray*}
%
Wenn Konvergenz f\"ur alle $x^{(0)}$ (also auch f\"ur alle $e^{(0)}$) gelten
soll, dann ist dies \"aquivalent zu
$\lim\limits_{k \to \infty} (M^{-1} N)^k = 0.$\\
Sei $M^{-1} N = PJP^{-1}$ die Transformation auf Jordan--Normalform, d.h.
%
\[
 J = \left[ \begin{array}{ccc}
  J_1 & & \\
  & \ddots & \\
  & & J_k \end{array}\right], \; 
 J_i = \left[ 
       \begin{array}{cccc}
	\lambda_i&  1   &      &  \\
	         &\ddots&\ddots&  \\
                 &      &\ddots& 1\\
                 &      &      & \lambda_i 
       \end{array}\right],
\]
%
so gilt $\lim\limits_{k \to \infty} (M^{-1}N)^k = 0 \Longleftrightarrow
\lim\limits_{k \to \infty} J^k = 0.$\\
$J$ ist obere $\bigtriangleup$--Matrix mit Eigenwerten auf der Diagonale.\\
Also mu{\ss} $(\lambda_i)^k \longrightarrow 0 \Longrightarrow |\lambda_i| < 1
\Longrightarrow \rho (J) < 1 \Longrightarrow \rho (M^{-1}N) < 1.$\\
Umgekehrt, falls $\rho (J) < 1$, so folgt $\lim\limits_{k \to \infty} J^k = 0.$
Man schaue sich dann die Form von $J^k$ an.
\eproof
%
M\"ogliche Wahlen f\"ur $M,N$. Setze 
\[
 \begin{array}{rcccccc}
  A &= &      D            &-&    L              &-& U \\
    &= & \left[\diag\right]&-&\left[\sltri\right]&-&\left[\sutri\right]
 \end{array}
\]
$D$ Diagonale, $L$ unteres $\bigtriangleup$, $U$ oberes $\bigtriangleup$.\\
%
$M = D, \; N = L+U$  ergibt das Verfahren
%
\eq{2.5.5}
 x^{(k+1)} = D^{-1} (L+U) x^{(k)} + D^{-1} b
\en
%
\begin{algorithmus}[Jacobi oder Gesamtschrittverfahren]\label{2.5.6}\quad\\
\begin{tabular}{llll}
\underline{Input:} & \multicolumn{3}{p{13.5cm}}{
			$A \in \K^{n,n}, b \in \K^n, \; a_{ii} \neq 0, 
			i = 1, \dots, n$ und Startvektor 
	                $x^{(0)} = [x^{(0)}_1,\dots,x^{(0)}_n]^T$
			sowie Abbruchgrenze  $\delta$.
		     }\\
\underline{Output:}& \multicolumn{3}{p{13.5cm}}{
			N\"aherung $\tilde{x}$ f\"ur L\"osung $Ax = b$
		     }\\[1.5ex]
&FOR&\multicolumn{2}{l}{$k=0,1,2,\dots$}\\
&   &FOR&$i=1:n$\\
&   &   & $\displaystyle x_i^{(k+1)} =\left(b_i
	   -\sum\limits^{i-1}_{j=1}a_{ij} x^{(k)}_j
           -\sum\limits^n_{j=i+1} a_{ij} x^{(k)}_j\right)/a_{ii}$; \\
&   &END&\\
&   &\multicolumn{2}{l}{IF $ \|x^{(k+1)} - x^{(k)} \| < \delta$,  STOP}\\
&END&
\end{tabular}
\end{algorithmus}
\underline{Kosten:} im wesentlichen ein Matrix--Vektor Produkt (unter 
		    Ausnutzung der Sparsit\"at), 1 Vektoraddition 
		    pro Iterationsschritt.\\
\underline{Fehleranalyse:} Wie bei Matrix * Vektor, Vektor+Vektor.
Verfahrensfehler abh\"angig von $\delta, \rho (D^{-1}(L+U))$.

$M = D-L, \; N = U$  ergibt das Verfahren
\eq{2.5.7}
 x^{(k+1)} = (D-L)^{-1} U x^{(k)} + (D-L)^{-1} b
\en

\begin{algorithmus}[Gau{\ss}--Seidel oder Einzelschrittverfahren]\label{2.5.8}\quad\\
\begin{tabular}{llll}
\underline{Input:} & \multicolumn{3}{p{13.5cm}}{
			$A \in \K^{n,n}, b \in \K^n, \; a_{ii} \neq 0, 
			i = 1, \dots, n$ und Startvektor 
	                $x^{(0)} = [x^{(0)}_1,\dots,x^{(0)}_n]^T$
			sowie Abbruchgrenze  $\delta$.
		     }\\
\underline{Output:}& \multicolumn{3}{p{13.5cm}}{
			N\"aherung $\tilde{x}$ f\"ur L\"osung $Ax = b$
		     }\\[1.5ex]
&FOR&\multicolumn{2}{l}{$k=0,1,2\dots$}\\
&   &FOR&$i=1:n$\\
&   &   & $\displaystyle x_i^{(k+1)} =\left(b_i
	   -\sum\limits^{i-1}_{j=1}a_{ij} x^{(k+1)}_j
           -\sum\limits^n_{j=i+1} a_{ij} x^{(k)}_j\right)/a_{ii}$; \\
&   &END&\\
&   &\multicolumn{2}{l}{IF $ \|x^{(k+1)} - x^{(k)} \| < \delta$,  STOP}\\
&END&
\end{tabular}
\end{algorithmus}

\underline{Kosten:} im wesentlichen wie Jacobi, nur schwerer als Vektoroperation
zu implementieren.

\underline{Fehleranalyse.} Fehler abh\"angig von Fehler in L\"osung von $(D-L)x = c$ sowie
Verfahrensfehler abh\"angig von $\rho ((D-L)^{-1} U), \delta.$
%

Beide Verfahren sind einfach, aber konvergieren i.a. sehr langsam, 
insbesondere bei den Problemen die in der Praxis auftauchen. In 
Spezialf\"allen gilt: Falls Jacobi konvergiert, so konvergiert
Gau{\ss}--Seidel schneller, falls Jacobi divergiert, so divergiert
Gau{\ss}--Seidel schneller (Satz von Stein/Rosenberg).
%
\begin{theorem}[Ostrowski/Reich]\label{2.5.9} 
Sei $A \in \R^{n,n}$ symmetrisch positiv definit, so
konvergiert das Gau{\ss}--Seidel Verfahren f\"ur alle $x^{(0)}$.
\end{theorem}
%
\proof
Siehe Stoer 2, Golub/Van Loan.
\eproof
%

\subsection{Konvergenzbeschleunigung}
Eine M\"oglichkeit, die Konvergenz zu beschleunigen ist die sukzessive 
\"Uberrelaxation (SOR).

$M = D-\omega L, N = (1 - \omega) D + \omega U, \omega \in(0,2)$, 
$(M-N) = \omega(D-L-U)$ ergibt Verfahren
%
\eq{2.5.10}
x^{(k+1)} = (D-\omega L)^{-1}((1 - \omega) D + \omega U)x^{(k)}+
	    (D-\omega L)^{-1}\omega b.
\en

\begin{algorithmus}[SOR]\label{2.5.11} \quad\\
\begin{tabular}{llll}
\underline{Input:} & \multicolumn{3}{p{13.5cm}}{
			$A \in \K^{n,n}, b \in \K^n, \; a_{ii} \neq 0, 
			i = 1, \dots, n$ und Startvektor 
	                $x^{(0)} = [x^{(0)}_1,\dots,x^{(0)}_n]^T$
			sowie Abbruchgrenze  $\delta$.
		     }\\
\underline{Output:}& \multicolumn{3}{p{13.5cm}}{
			N\"aherung $\tilde{x}$ f\"ur L\"osung $Ax = b$
		     }\\[1.5ex]
&FOR&\multicolumn{2}{l}{$k=0,1,2\dots$}\\
&   &FOR&$i=1:n$\\
&   &   & $\displaystyle x_i^{(k+1)} =
           \omega \left( b_i - \sum\limits^{i-1}_{j=1} a_{ij} x_j^{(k+1)}
           -\sum\limits^n_{j=i+1} a_{ij} x^{(k)}_j \right) / a_{ii} 
           + (1-\omega) x^{(k)}_i$;\\
&   &END&\\
&   &\multicolumn{2}{l}{IF $ \|x^{(k+1)} - x^{(k)} \| < \delta$,  STOP}\\
&END&
\end{tabular}
\end{algorithmus}
Kosten und Fehler wie oben ($\omega=1$: Gau{\ss}--Seidel).
 
%

In einigen Spezialf\"allen von strukturierten Matrizen aus partiellen
Differentialgleichungen gibt es Konvergenzbeweise, optimales $\omega$, etc.
Weitere Beschleunigungsmethoden: Tschebyschow--Beschleunigung, 
semiiterative Verfahren, unvollst\"andige Dreieckszerlegungen.

\subsection{Konjugierte Gradienten}
\underline{Grundidee:} Betrachte quadratisches Funktional
\[ 
\Phi (x) = \frac{1}{2} x^T Ax - x^T b, \;A \in \R^{n,n} 
\mbox{ symmetrisch, positiv definit}, \;b \in \R^n.
\]
Minimum von $\Phi$ ist der Wert $-b^TA^{-1}b/2$. Dieser
wird erreicht bei $x = A^{-1} b.$ Also ist Minimierung von $\Phi(x)$
\"aquivalent zur L\"osung von $Ax=b.$

\underline{Methode des Steilsten Abstiegs:} In jedem Schritt lege 
die Richtung des steilsten Abstiegs fest und gehe in diese Richtung
so weit wie m\"oglich. Beim Punkt $x_c$ angelangt, ist die Richtung des 
steilsten Abstiegs die Richtung des negativen Gradienten.
%
\eq{2.5.12}
 - \left.\nabla \Phi \right|_{x_c} = b - Ax_c = r_c, \;
   \mbox { Residium von } x_c.
\en
%
Falls $r_c \neq 0$, dann gibt es ein $\alpha$, 
so da{\ss} $\Phi (x_c + \alpha r_c) < \Phi (x_c).$ W\"ahle
%
\eq{2.5.13}
 \alpha = \frac{ r^T_c r_c}{r^T_c Ar_c}.
\en
%
Das ergibt Minimum von $\Phi (c_x + \alpha r_c)$ \"uber
$\alpha.$ Man h\"atte damit folgende Methode des steilsten Abstiegs.

\begin{tabular}{ll}
\multicolumn{2}{l}{$k=0; x_0 = 0; r_0 = b;$}\\
WHILE& $r_k \neq 0$\\
     & $k=k+1$;\\
     & $\alpha_k  =  r^T_{k-1} r_{k-1} / r^T_{k-1} Ar_{k-1}$;\\
     & $x_k       =  x_{k-1} + \alpha_k r_{k-1}$;\\
     & $r_k       =  b - Ax_k;$\\
END&
\end{tabular}

Sehr langsame Konvergenz wenn $\kappa_2(A)$ gro{\ss}. 
Besser ist es, die Suchrichtungen
so zu w\"ahlen, da{\ss} die Residuen senkrecht aufeinanderstehen, so da{\ss} man in exakter 
Rechnung in $\le n$ Schritten fertig ist.
%
\begin{algorithmus}[Konjugierte Gradienten--Verfahren]\label{2.5.14}\quad\\
\begin{tabular}{llll}
\underline{Input:} & \multicolumn{3}{p{13.5cm}}{
			$A \in \R^{n,n}$, symmetrisch, positiv definit, 
                        $b \in \R^n$, $k_{\max}$ (Maximum an 
			Iterationsschritten), $\varepsilon>0$ Toleranz.
		     }\\
\underline{Output:}& \multicolumn{3}{p{13.5cm}}{
	              L\"osung von $Ax = b$.
		     }\\[1.5ex]
&\multicolumn{3}{l}{$k=0;\; x_0 = 0;\; r_0 = b;\; \rho_0 = \|r_0\|^2_2;$}\\
&WHILE& \multicolumn{2}{l}{$\sqrt{\rho_k}>\varepsilon\sqrt{\rho_0}$ AND 
	$k<k_{\max}$}\\
&     & \multicolumn{2}{l}{$k=k+1$;}\\
&     & IF & $k=1$\\
&     &    & $p=r;$\\
&     &ELSE&\\
&     &    & $\beta  =  \rho_{k-1} / \rho_{k-2};$  
             \% $( = r_{k-1}^T r_{k-1} /r^T_{k-2} r_{k-2})$\\
&     &    & $p = r + \beta p;$ \% $(= p_k)$\\
&     &END &\\
&     & \multicolumn{2}{l}{$w =Ap$;}\\
&     & \multicolumn{2}{l}{$\alpha = \rho_{k-1}/p^T w$;  
	\% $(=r^T_{k-1} r_{k-1} /p_k^T Ap_k)$}\\
&     & \multicolumn{2}{l}{$x  =  x + \alpha p$;  \% $(= x_k)$}\\
&     & \multicolumn{2}{l}{$r = r - \alpha w$; \% $(=r_k)$}\\
&     & \multicolumn{2}{l}{$\rho_k  =  \|r\|^2_2$; \% $(=r^T_k r_k)$}.\\
& END  
\end{tabular}
\end{algorithmus}
%
\underline{Kosten:} 1 Matrix * Vektor--Multiplikation  + 2 Skalarprodukte
+ 2 GAXPY pro Schritt.\\
\underline{Fehleranalyse:} Rundungsfehler zerst\"oren Konvergenz in $n$
Schritten und Orthogonalit\"at der Residuen.
%
Bemerkung: 
\begin{itemize}
\item
 F\"ur die Vektoren $p\equiv p_k$, $r\equiv r_k$ aus Algorithmus \ref{2.5.14} 
 gilt:
 $p_k^TAp_l = 0, r_k^Tr_l=0$ f\"ur alle $k>l\geqslant 0$.\\
\item 
 F\"ur die Funktion $\Phi$ gilt mit $p_1,\dots,p_k$ aus 
 Algorithmus \ref{2.5.14}:
\[
\min_{\alpha_1,\dots,\alpha_k\in\R} \Phi(x+\sum_{l=1}^k\alpha_l p_l)=
\min_{\alpha_k\in\R} \left(\min_{\alpha_1,\dots,\alpha_{k-1}\in\R} \Phi(x+\sum_{l=1}^k\alpha_l p_l)\right).
\]
\end{itemize}
 $\rightarrow$ \"Ubung \ref{numI/uebung042}.

\begin{theorem}\label{2.5.15} 
Sei $A \in \R^{n,n}$ symmetrisch positiv definit und $b \in \R.$ 
F\"ur die Iterierten in Algorithmus \ref{2.5.14} gilt:
%
\eqa
\nonumber \| x - x_k\|_A & \equiv & \sqrt{(x-x_k)^T A(x-x_k)}\\
\label{2.5.16}           &\le& 2\|x-x_0\|_A 
			  \left( \frac{\sqrt{\kappa_2(A)}-1}
                                      {\sqrt{\kappa_2(A)}+1} \right)^k.
\ena
\end{theorem}
%
\proof
Siehe Golub/Van Loan.
\eproof
%
Genauigkeit ist oft beser als (\ref{2.5.16}) vorhersagt. Je
n\"aher $\kappa_2 (A)$ an $1$, je schneller die Konvergenz. 
Konvergenzbeschleunigung
durch Pr\"akonditionierung. \\
L\"ose
\eq{2.5.17}
 (C^{-1} A C^{-T})(C^Tx) = C^{-1}b
\en
%
anstatt $Ax=b$, wobei $C$ nicht singul\"ar und so gew\"ahlt,
da{\ss} $\kappa_2 (C^{-1}AC^{-T})$ m\"oglichst nahe an $1$ liegt.
Best m\"ogliche Wahl:
\[
 C^{-1} = A^{-\frac{1}{2}} \Longrightarrow \kappa_2 (C^{-1}AC^{-T})=1.
\]
%
Das ist nat\"urlich genau so schwierig wie urspr\"ungliches Problem.\\
In vielen Anwendungen wie etwa Partielle Differentialgleichungen,
Finite Element Methoden hat man nat\"urliche Pr\"akonditionierer. 
Pr\"akonditioniertes Kongugiertes Gradientenverfahren 
f\"ur symmetrisch positiv definite Systeme i.a. das Beste. 




 




