\section{Ausgleichsprobleme}\label{KAPITEL-02.4}
%
Bei vielen wissenschaftlichen Versuchen geht es darum, aus Messungen die
Werte von Konstanten $x_1, \dots , x_n$ zu bestimmen. Oft kann man
$x_i$ nicht direkt messen sondern nur eine leichter zug\"angliche Gr\"o{\ss}e
$y$, die als Funktion von $x$ und anderen Parametern $z$ abh\"angt,
$y = f(z, x_1, \dots , x_n).$
%
\begin{example}\label{2.4.1}
Bestimmung von $g$.\\ 
Fallversuche $h = g \frac{t^2}{2}$, $h$ -- Fallh\"ohe, $t$ -- Fallzeit, 
$g$ -- Gravitation. 

Um $x_i$ zu bestimmen, f\"uhre $m > n$ verschiedene Experimente durch und 
erhalte
\[
 y_k = f(z_k, x_1, \dots, x_n),\; k = 1,\dots, m.
\]
%
Dies ergibt \"uberbestimmtes Gleichungssystem. Dies ist  i.a. nicht l\"osbar. 
Daher nur beste Approximation m\"oglich. Zum Beispiel minimiere:
%
\eq{2.4.2}
 \sum^m_{k=1} \left( y_k - f (z_k, x_1, \dots , x_n) \right)^2, \;  
 \mbox  {  kleinste Quadrate (least square), } (\| \bullet \|_2)
\en
%
oder
\[
 \max\limits_{1 \le k \le m} \left| y_k - f (z_k, x_1 , \dots , x_n) \right|,
 \;  (\|\bullet\|_\infty). 
\]
%
f\"ur (\ref{2.4.2}) ergibt sich die notwendige Bedingung 
f\"ur Minimierung durch
%
\eq{2.4.3}
\frac{\partial}{\partial x_i} \left( \sum^m_{k=1} (y_k - f_k (x_1, \dots , x_n))^2
 \right) = 0,  \; i=1:n, 
\en
wobei   $f_k (x_1, \dots , x_n) := f(z_k, x_1, \dots , x_n )$.
%
\end{example}

Dies sind die sogenannten \underline{Normalengleichungen}. 
Ein wichtiger
Spezialfall, das lineare Ausgleichsproblem, liegt vor falls 
$f_k (x_1, \dots , x_n)$ linear in $x_1,\dots,x_n$,
d.h. 
%
\eq{2.4.4}
 \left[ \begin{array}{c}
 f_1 (x_1 , \dots , x_n) \\
 \vdots \\
 f_m (x_1 , \dots , x_n) \end{array}\right] = Ax, \; A \in \K^{m,n}.
\en
%
Dann sind die Normalengleichungen 
\eq{2.4.5}
 \mbox {grad}_x \left( (y-Ax)^* ( y-Ax) \right) = 2A^*Ax - 2A^* y = 0.
\en
%
Also man hat das lineare Gleichungssystem
\eq{2.4.6}
 A^*Ax = A^* y
\en
%
zur L\"osung von
\eq{2.4.7}
 \min\limits_{x \in \K^n} \| Ax-y\|_2.
\en
%
%
\begin{theorem}\label{2.4.8}
Das lineare Ausgleichsproblem (\ref{2.4.7}) hat mindestens eine L\"osung $x_0$.
Ist $x_1$ eine weitere L\"osung, so gilt $Ax_0 = Ax_1$. Das Residuum 
$r = y - Ax_0$ ist eindeutig bestimmt und erf\"ullt $A^* r = 0.$ Jede L\"osung
$x_0$ erf\"ullt auch (\ref{2.4.6}) und umgekehrt.
\end{theorem}
%
\proof
$\K^m =  \,\underbrace{\mathrm{Bild}\,(A)}_L \oplus
\,\underbrace{\mathrm{Bild}\,(A)^\perp}_{L^\perp}.$\\
Daher gibt es eine eindeutige Zerlegung
%
\eq{2.4.9}
 y = s + r,  \mbox {  mit } s \in L, r \in L^\bot
\en
%
und es gibt ein $x_0$ mit $A^*x_0 = s.$ Da $A^*r = 0$, so folgt
$A^*y = A^* s = A^*Ax_0.$ Also l\"ost $x_0$ (\ref{2.4.6}) die Normalengleichungen.\\
Umgekehrt sei $x_1$ L\"osung von (\ref{2.4.6}). Dann hat $y$ die Zerlegung 
$y = s+r$ mit $s = Ax_1, r = y-Ax_1, \;  s \in L, r \in L^\bot.$ Wegen der 
Eindeutigkeit der Zerlegung folgt $Ax_1 = Ax_0$. \\
Weiter erf\"ullt jede  L\"osung $x_0$ von (\ref{2.4.6})
das Ausgleichsproblem (\ref{2.4.7}). Denn
f\"ur beliebiges $x$ setze $z = Ax - Ax_0, \;r = y-Ax_0$, dann gilt wegen
$r^* z = 0$:
%
\[
 \|y - Ax\|^2_2 = \|r - z\|^2_2 = \|r\|^2_2 + \|z\|^2_2 \ge
 \|r\|^2_2 = \|y-Ax_0\|^2_2.
\]
\eproof
%
Falls $A$ vollen Rang hat, so ist $A^*A$ positiv definit. Nun k\"onnte man einfach
folgenden Algorithmus verwenden:
%
\begin{algorithmus}[Normalengleichung]\label{2.4.10}\quad\newline
\begin{tabular}{ll}
\underline{Input:} & \multicolumn{1}{p{13.5cm}}{
			$A \in \K^{m,n}$ mit $\rnk(A) = n$, und $b \in \K^m.$
		     }\\
\underline{Output:}& \multicolumn{1}{p{13.5cm}}{
			L\"osung $x$ des Minimierungsproblems (\ref{2.4.7}).
		     }\\[1.5ex]
& Berechne das untere Dreieck von $C = A^*A.$\\
& $d = A^*b.$\\
& Berechne Cholesky Zerlegung von $C = GG^*.$ \\
& L\"ose $Gy = d$ und $G^*x = y.$
\end{tabular}
\end{algorithmus}
\underline{Kosten:} $(m+n/3)n^2$ \flops.\\
\underline{Fehleranalyse:} Angenommen, keine Fehler bei $C = A^*A$ und 
			   $d = A^*b.$
Dann gilt wegen der Analyse der Cholesky Zerlegung:
\[
 (A^*A+E) \tilde{x} = A^* b, 
\]
mit $\|E\|_2 \approx \eps \|A^*\|_2 \|A \|_2 \approx \eps \|A^*A\|_2.$,
d.h. wir erwarten
%
\[
  \frac{\|x-\tilde{x}\|_2}{\|x\|_2} \approx \eps  \kappa_2 (A^*A)
  = \eps \kappa_2 (A)^2,
\]
%
d.h. das Quadrat der Kondition geht ein. {\bf Schlecht!}\reversemarginpar\marginpar{\huge \hfill !}
%
\begin{example}\label{2.4.11}
\[
 A = \left[ \begin{array}{cc}
  1 & 1\\
  10^{-3} & 0\\
  0 & 10^{-3} \end{array}\right],\; b = \left[ \begin{array}{c}
  2 \\ 10^{-3} \\ 10^{-3} \end{array}\right],
  \; x = \left[ \begin{array}{c}
 1 \\ 1 \end{array}\right], \; \kappa_2 (A) = 1.4 \cdot 10^3.
\]
%
Mit 6--stelliger Arithmetik ist $A^*A = \left[ \begin{array}{ll}
  1 & 1\\
  1 & 1 \end{array}\right]$ singul\"ar. Mit 7--stelliger Arithmetik folgt
%
\[
 \tilde{x} = \left[ \begin{array}{c}
  2.00001\\
  0 \end{array}\right] \Longrightarrow \frac{\| \tilde{x} - x \|_2}{\|x\|_2}
  \approx \underbrace{\eps}_{10^{-6}} \underbrace{\kappa_2 (A)^2}_{10^6}.
\]
%
\end{example}

Normalengleichung kann also zu sehr schlechten Ergebnissen f\"uhren. 
Verwende daher $QR$--Zerlegung. Bilde
%
\[
\begin{array}{rcl}
 Q^*A &= &R = 
          \left[ \begin{array}{c}
           R_1 \\ 0 \end{array}\right] 
          \begin{array}{l}
           n \\ m-n \end{array}, \; 
           R_1 \mbox { obere } \bigtriangleup \mbox {--Matrix}\\
%
 Q^*b &=& \left[ \begin{array}{c}
  c \\d \end{array}\right] \begin{array}{l}
  n \\ m-n \end{array}
\end{array}
\]
%
Dann gilt
%
\[
 \|Ax - b\|^2_2 = \|Q^*Ax - Q^*b\|^2_2 =
 \|R_1x - c\|^2_2 + \|d\|^2_2,
\]
%
da $Q$ unit\"ar, f\"ur alle $x\in\K^n.$ Falls $\rnk(A)$ = $\rnk(R_1)=n$,
so gilt f\"ur die L\"osung $R_1x = c.$
%
Also hat man: 
%
\begin{algorithmus}[$QR$ L\"osung des Ausgleichsproblem]\label{2.4.12}

\begin{tabular}{lll}
\underline{Input:} & \multicolumn{2}{p{13.5cm}}{
			$A \in \K^{m,n}$ mit $\rnk(A) = n$, und $b \in \K^m.$
		     }\\
\underline{Output:}& \multicolumn{2}{p{13.5cm}}{
			$x \in \K^n$, so da{\ss} 
			$\|Ax-b\|_2 \stackrel{!}{=} \min$.
		     }\\[1.5ex]
& \multicolumn{2}{p{13.5cm}}{
   Verwende Algorithmus \ref{2.3.6}, um $A$ mit seiner $QR$--Zerlegung zu 
   \"uberschreiben.}\\
&FOR& $j=1:n$\\ 
&   & $v(j) = 1;$\\
&   & $v(j+1:m) = A(j+1:m,j);$\\
&   & $b(j:m) = ROWHOUSE(b(j:m), v(j:m));$\\
&END&\\
& \multicolumn{2}{p{13.5cm}}{
   L\"ose $R(1:n,1:n)x=b(1:n)$ mit R\"uckw\"arts--Einsetzen.}
\end{tabular}
\end{algorithmus}
\underline{Kosten:} $2n^2 (m - \frac{n}{3})$ \flops.\\
\underline{Fehleranalyse:} Das berechnete $\tilde{x}$ l\"ost 
$\min  \|(A+ \delta A)x - (b+ \delta b)\|_2$ wobei
%
\begin{eqnarray*}
 \| \delta A\|_F & \le & (6m - 3n +41) n \eps  \|A\|_F + {\cal O}
 (\mbox {eps}^2)\\
 \| \delta b \|_2 & \le & (6m - 3n + 40) n \eps  \|b\|_2 +
 {\cal O} (\eps^2)
\end{eqnarray*}
%
{\bf Vorsicht!} Probleme tauchen auf, wenn 
$\kappa_2(A) = \kappa_2(R) \approx \frac{1}{\eps}$, oder nat\"urlich
wenn Rang$(A) < n$.\reversemarginpar\marginpar{\huge \hfill !}

Abhilfe schafft hier die sogenannte Singul\"arwertzerlegung $A=U\Sigma V^*$
mit $U,V$ unit\"ar, $\Sigma= \left[ 
        \begin{array}{ccc}
         \sigma_1&        &        0\\
	         & \ddots &         \\ 
             0   &        & \sigma_n\\ \hline
      & 0 & \end{array}\right].$ 
Dies kann zu diesem Zeitpunkt hier nicht gemacht werden.

