TUD_MATH_BA/4. Semester/NUMM/TeX_files/Fixpunktiteration.tex

Seien eine reguläre Matrix $A \in \R^{n \times n}$ und $b \in Rn$ gegeben. In diesem Kapitel werden iterative Verfahren zur Lösung des linearen Gleichungssystems
\begin{align}
Ax = b\label{eq_2_2_1}
\end{align}
betrachtet.
%TODO fix counter, this equation should be the 1 and the one in the following section, 2!
\section{Fixpunktiteration}

Grundidee dieser Verfahren ist die geeignete Umformulierung des System $Ax = b$ als Fixpunktaufgabe und die Anwendung des gewöhnlichen Iterationsverfahrens. Die hier betrachtete (zu \cref{eq_2_2_1} äquivalente) Fixpunktaufgabe lautet
\begin{align}
	x = x - B^{-1}(Ax - b),\notag
\end{align}
wobei $B \in \R^{n \times n}$ eine noch zu wählende reguläre Matrix ist. Bei Wahl eines Startpunktes $x^0 \in \Rn$ ergibt sich das gewöhnliche Iterationsverfahren damit zu
\begin{align}
	x^{k+1} := x^{k} - B^{-1}(A x^k -b) = (I - B^{-1}A)x^{k} + B^{-1}b, \qquad k = 0,1,2,\dots \label{eq_2_2_2} %TODO add underbraces for M = I - B^{-1}A and c = B^{-1}b
\end{align}
Mit den Bezeichnung $M := I - B^{-1}A$ und $c:= B^{-1}b$ untersuchen wir deshalb die Iterationsvorschrift
\begin{align}
	x^{k+1} := Mx^k + c. \label{eq_2_2_3}
\end{align}
Die zugehörige Fixpunktabbildung $\Phi: \Rn \to \Rn$ ist damit offenbar gegeben durch
\begin{align}
	\Phi(x) := Mx + c.\notag
\end{align}
\begin{proposition}
	\proplbl{2_2_1}
	Es sei $B \in \Rnn$ regulär und mit $M:= I - B^{-1}A$ gelte
	\begin{align}
		\lambda := \norm{M}_{\ast} < 1 \label{eq_2_1_4}
	\end{align}
	wobei $\norm{\cdot}_{\ast}$ die einer Vektornorm $\norm{\cdot}$ zugeordnete Matrixnorm bezeichnet. Dann gilt:
	\begin{enumerate} %TODO add alph
		\item Die für eine beliebiges $x^0 \in \Rn$ durch \cref{eq_2_2_3} erzeugte Folge $\set{x^k}$ konvergiert gegen die eindeutige Lösung $x*$ des linearen Gleichungssystems \cref{eq_2_2_1}.
		\item Die Abschätzungen \cref{eq_1_2_2} - \cref{eq_1_2_4} sind für alle $k \in \N$ erfüllt.
	\end{enumerate}
\end{proposition}

\begin{proof}
	Direkte Folgerung aus dem Banachschen Fixpunktsatz (\propref{1_1_1})
\end{proof}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3rd lecture %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{remark}
	\proplbl{2_2_2}
	In \propref{2_2_1}a) kann die Folgerung \cref{eq_2_1_4} durch die Bedingung
	\begin{align}
	\rho(M) < 1\label{eq_2_1_5}
	\end{align}
	ersetzt werden. Da
	\begin{align}
		\rho(C) \le \norm{C}_{\ast} \quad \text{ für alle } C \in \Rnn \notag
	\end{align}
	für jede beliebige zugeordnete Matrixnorm $\norm{\cdot}_{\ast}$ gilt (vgl. Übungsaufgabe), ist \cref{eq_2_1_5} eine schwächere Forderung als \cref{eq_2_1_4}. Andererseits gibt es zu jedem Paar $(C,\epsilon) \in \Rnn \times (0,\infty)$ eine zugeordnete Matrixnorm $\norm{\cdot}_{(C,\epsilon)}$, so dass
	\begin{align}
		\norm{C}_{(C,\epsilon)} \le \rho(C) + \epsilon. \notag
	\end{align}
	Dabei ist $\rho(C)$ der \begriff{Spektralradius} der Matrix $C \in \Rnn$, d.h.
	\begin{align}
		\rho(C) := \max_{i = 1,\dots,n}\abs{\lambda_i}, \notag
	\end{align}
	wobei $\lambda_1,\dots,\lambda_n \in \C$ die Eigenwerte der Matrix $C \in \Rnn$ bezeichnen. Man kann weiter zeigen, dass \cref{eq_2_1_5} auch notwendig dafür ist, dass die durch \cref{eq_2_2_2} erzeugte Folge $\set{x^k}$ für jedes $x^0$ gegen $x*$ konvergiert.
\end{remark}

Um eine Matrix $B$ zu finden, so dass einerseits der Aufwand pro Iteration \cref{eq_2_2_2} niedrig und andererseits die Bedingung \cref{eq_2_1_4} bzw. \cref{eq_2_1_5} erfüllt ist, betrachten wir die folgende Zerlegung
\begin{align}
	A = L + D + R \notag
\end{align}

der Matrix $A$, wobei $D:= \diag(a_{11}, \dots, a_{nn})$ die aus den Diagonalelementen von $A$ bestehende Diagonalmatrix bezeichnet und $L$ bzw. $R$ eine untere bzw. obere Dreiecksmatrix ist mit

\begin{align}
	L = 
	\begin{pmatrix}
		0      &        &        &           & \\
		a_{21} & 0      &        &           & \\
		a_{31} & a_{32} & 0      &           & \\
		\vdots &        & \ddots & \ddots    & \\
		a_{n1} & \cdots & \cdots & a_{n,n-1} & 0
	\end{pmatrix}
	 \bzw R = 
	 \begin{pmatrix}
	 0      & a_{12} & a_{13} & \dots     & a_{1n}\\
	        & 0      & a_{23} & \dots     & a_{2n}\\
	        &        & \ddots & \ddots    & \vdots\\
	        &        &        & 0         & a_{n-1,n}\\
	        &        &        & a_{n,n-1} & 0
	 \end{pmatrix}.\notag
\end{align}

\subsection{Das \person{Jacobi}-Verfahren}
Wir setzen hier voraus, dass $D$ regulär ist und wählen
\begin{align}
	B:= D \label{eq_1_2_6}
\end{align}
Damit ergibt sich die Iterationsvorschrift
\begin{align}
	x^{k+1} = x^k - D^{-1}(Ax^k - b) = -D^{-1}(L+R)x^k+D^{-1}b. \label{eq_1_2_7}
\end{align}
In \cref{eq_2_2_3} ist entsprechend
\begin{align}
	M:= M_J := -D^{-1}(L+R) \text{ und } c:= c_J := D^{-1}b\notag 
\end{align}
zu wählen. Dieses Verfahren heißt \begriff{Gesamtschrittverfahren} oder \begriff{Jacobi-Verfahren}. Der Aufwand pro Schritt (Berechnung von $x^{k+1}$ aus $x^k$) beträgt $\Landau(n^2)$ bei voll besetzter Matrix $A$ und mindestens $\Landau(n)$, falls $A$ schwach besetzt ist.

\begin{proposition}
	Die Matrix $A$ sei streng diagonaldominant (vgl. Definition 3.1 der Vorlesung ENM). Dann ist die Matrix $B$ aus \cref{eq_1_2_6} regulär und es gilt
	\begin{align}
	\norm{M_J}_{\infty} \le \lambda_{SD} := \max_{i = 1,\dots,n} \frac{1}{\abs{a_{ii}}} \sum_{\substack{j =1 \\ j\neq i}}^n \abs{a_{ij}} < 1. \notag
	\end{align}
\end{proposition}

\begin{proof}
	Die Regularität von $B$ ergibt sich sofort aus der strengen Diagonaldominanz von $A$. Nutzt man die Definition der Zeilensummennorm $\norm{\cdot}_{\infty}$ erhält man sofort
	\begin{align}
		\norm{M_J}_{\infty} = \norm{D^{-1}(L+R)}_{\infty} =\max_{i = 1,\dots,n} \frac{1}{\abs{a_{ii}}} \sum_{\substack{j =1 \\ j\neq i}}^n \abs{a_{ij}}  = \lambda_{SD}.\notag
	\end{align}
	Die vorrausgesetzte strenge Diagonaldominanz von $A$ sichert $\lambda_{SD} < 1$.
\end{proof}

\subsection{Das \person{Gauss-Seidel}-Verfahren}
Wir setzen hier voraus, dass $L + D$ regulär ist und wählen

\begin{align}
	B := L + D \label{1_2_8}
\end{align}

Damit ergibt sich die Iterationsvorschrift

\begin{align}
	x^{k+1} = x^k - (L+D)^{-1}(Ax^k - b) = - (L+D)^{-1}R x^k + (L+D)^{-1}b. \label{1_2_9}.
\end{align}

In \cref{eq_2_2_3} ist entsprechend

\begin{align}
	M:= M{GS} := - (L+D)^{-1}R \text{ und } c:= c_{GS} := (L+D)^{-1}b \notag
\end{align}

zu wählen. Dieses Verfahren heißt \begriff{Einzelschrittverfahren} oder \begriff{Gauß-Seidel-Verfahren}. Der Aufwand pro Schritt beträgt im ungünstigsten Fall $\Landau(n^2)$. Verbesserungen sind möglich, wenn eine Sparse-Struktur in $A$ ausgenutzt werden kann.

\begin{proposition}
	Die Matrix $A$ sei streng diagonaldominant ($\nearrow$ Definition 3.1 der Vorlesung ENM). Dann ist die Matrix $B$ aus \cref{1_2_8} regulär und es gilt
	\begin{align}
		\norm{M_{GS}}_{\infty} \le \lambda_{SD} <1. \notag
	\end{align}
\end{proposition}

\begin{proof}
	Die Regularität von $B$ folgt sofort aus der strengen Diagonaldominanz von $A$. Weiter ergibt sich
	\begin{align}
		\norm{M_{GS}}_{\infty} = \norm{(L+D)^{-1}R}_{\infty} = \sup_{\norm{y}_{\infty}=1} \norm{(L+D)^{-1}Ry}_{\infty}. \notag
	\end{align}
	Um für einen festen Vektor $y$ mit $\norm{y}_{\infty} = 1$ eine Abschätzung für die rechte Seite zu erhalten, setzen wir $z:= (L+D)^{-1}Ry$. Damit gilt
	\begin{align}
		(D+L)z = Ry \label{1_2_10}
	\end{align}
	und
	\begin{align}
	z_1 = \frac{1}{a_{11}} \sum_{j=1}^{n} a_{1j}y_j. \notag
	\end{align}
	Daraus folgt (da $\lambda_{SD} < 1$ wegen der strengen Diagonaldominanz von $A$)
	\begin{align}
		\abs{z_1} 
		\le \frac{1}{\abs{a_{11}}} \sum_{j=2}^{n} \abs{a_{1j}}\abs{y_j} 
		\le \sum_{j=2}^{n} \abs{a_{1j}} \le \lambda_{SD} < 1.\notag
	\end{align}
	Nehmen wir nun an, dass
	\begin{align}
		\abs{z_1} \le \text{ für } i = 1, \dots, k-1, \notag
	\end{align}
	für ein $k \in \set{2,\dots,n}$ gilt. Dann folgt wegen \cref{1_2_10} und $\norm{y}_{\infty} = 1$
	\begin{align}
		\abs{z_k} = \frac{1}{\abs{a_{kk}}} \abs{-\sum_{i=1}^{k-1} a_{ki}z_i + \sum_{i=k+1}^{n} a_{ki}y_i} 
		\le \frac{1}{\abs{a_{kk}}} \brackets{\sum_{i=1}^{k-1} \abs{a_{ki}} + \sum_{i=k+1}^{n} \abs{a_{ki}}} \le \lambda_{SD}. \notag
	\end{align}
	Somit hat man induktiv $\abs{z_k} \le \lambda_{SD}$ für $k = 1, \dots, n$ und damit
	\begin{align}
		\norm{(L+D)^{-1}Ry}_{\infty} = \norm{z}_{\infty} \le \lambda_{SD} \notag
	\end{align}
	für beliebige $y$ mit $\norm{y}_{\infty} = 1$.
\end{proof}

\subsection{SOR-Verfahren}

Um dieses verfahren zu beschreiben, nehmen wir an, dass für ein $\omega \neq 0$ die Matrix

\begin{align}
	B:=L + \frac{1}{\omega}D \label{eq_2_2_11}
\end{align}
regulär ist. Damit ergibt sich die Iterationsvorschrift

\begin{align}
	x^{k+1} := x^k - \brackets{L+\frac{1}{\omega}D}^{-1}(Ax^k - b) = M(\omega)x^k + c(\omega)\notag
\end{align}

\begin{align}
	M(\omega) := I -  \brackets{L + \frac{1}{\omega} D^{-1} A} 
	= \brackets{L + \frac{1}{\omega} D}^{-1}
%	\brackts{\frac{1-\omega}{\omega} D - R} \notag
\end{align}

und

\begin{align}
	c(\omega) :=  \brackets{L+\frac{1}{\omega}D}^{-1}b.
\end{align}

Für $\omega = 1$ erhält man offenbar als Spezialfall das Gauß-Seidel-Verfahren, so dass der folgende Satz auch dafür Anwendung finden kann. Man beachte dazu \propref{2_2_2}.

\begin{proposition}
	Die Matrix $A$ sei symmetrisch und positiv definit. Dann ist die Matrix $B$ aus \cref{eq_2_2_11} regulär (für jedes $\omega \neq 0$). Falls $\omega \in (0,2)$, dann gilt
	\begin{align}
		\rho(M(\omega)) < 1\notag
	\end{align}
	und umgekehrt.
\end{proposition}

\begin{proof}
	Da $A$ positiv definit ist, gilt $e_i^TA e_i = a_{ii} > 0$ für $i = 1, \dots, n$. Also ist $D$ positiv definit und damit $B$ regulär für alle $\omega \neq 0$.\\
	Sei $\lambda \in \C$ ein Eigenwert von $M(\omega)$ und $z \in \Cn$ ein zugehöriger Eigenvektor. Mit
	\begin{align}
		A = A - M(\omega)^T AM(\omega) + M(\omega)^T AM(\omega) \notag
	\end{align}
	sowie (unter Berücksichtigung der Definition von $M$ und von $A = A^T$ und $R = L^T$)
	\begin{align}
	A - M(\omega)^T AM(\omega) 
	&= A - (I - B^{-1}A)^T A(I-B^{-1}A)\notag\\
	&= AB^T A + AB^{-1}A - AB^{T-1}AB^{-1}A\notag\\
	&= (B^{-1}A)^T (B + B^T - A)(B^{-1}A)\notag\\
	&= (B^{-1}A)^T \brackets{L + \frac{1}{\omega}D + L^T + \frac{1}{\omega}D - L - D - L^T})(B^{-1}A)\notag\\
	&= (B^{-1}A)^T\brackets{\frac{2 -\omega}{\omega}D}(B^{-1}A)\notag
	\end{align}
	ergibt sich daher
	\begin{align}
		z^H Az = (AB^{T-1}z)^H\brackets{\frac{2 -\omega}{\omega}D}(B^{-1}Az) + z^H M(\omega)^T AM(\omega)z.\notag
	\end{align}
	Da die Diagonalmatrix $D$ positiv-definit ist, besitzt $\frac{2-\omega}{\omega}D$ dieselbe Eigenschaft für $\omega \in (0,2)$. Es folgt
	\begin{align}
		(AB^{T-1}z)^H \brackets{\frac{2 -\omega}{\omega}D} (B^{-1 Az}) > 0 \notag
	\end{align}
	und damit
	\begin{align}
		\abs{\lambda} < 1.
	\end{align}
	Also gilt $\rho(M(\omega)) = \max_{i = 1,\dots,n}\abs{\lambda_i} < 1$, sofern $\omega \in (0,2)$. Die Umkehrung der Aussage ergibt sich aus dem Satz von \person{Kahan} ($\nearrow$ Übungsaufgabe). %TODO add later.
\end{proof}

Es ist nun naheliegend, dass man $\omega \in(0,2)$ so wählen möchte, dass $\rho(\omega)$ möglichst klein ist. Dies ist in bestimmten Fällen näherungsweise möglich, ansonsten beschränkt man sich auf geeignete Heuristiken zur Wahl von $\omega$. Auf der Fixpunktiteration \cref{eq_2_2_3} beruhende Verfahren werden häufig auch \begriff{Splitting-Methoden} genannt. Es gibt noch weitere solche Verfahren, auf die hier nicht eingegangen wird.