Ausgleichsprobleme

 

 Das lineare Ausgleichsproblem

Gegeben sei eine reelle Funktion $ g$ die neben der Variablen $ t$ noch von einem unbekannten Vektor $ x=(x_1,\hdots,x_n)^T$ abhängt. Dabei habe $ g$ die Gestalt

$\displaystyle g(t,x)=\sum_{i=1}^n x_i g_i(t)$

$ g$ hängt also linear von den Komponenten von $ x$ ab.

Man kann nun z.B. durch ein Experiment versuchen die unbekannten $ x_i$ zu bestimmen. Dafür macht man $ n$ Beobachtungen, d.h. man notiert sich zu jedem getesteten $ t_i$ den Ausgang $ b_i$ des Experiments. Man erhält dann ein lineares Gleichungssystem in den $ n$ Unbekannten:

$\displaystyle g(t_i,x)=b_i$    für $\displaystyle i=1,...,n$
beziehungsweise
\begin{displaymath}\begin{array}{ccccccccc} x_1g_1(t_1)&+&x_2g_2(t_1)&+&...&+&x_... ...g_1(t_n)&+&x_2g_2(t_n)&+&...&+&x_ng_n(t_n)&=&b_n\ \end{array}\end{displaymath}

Schreibt man $ A=(a_{ij})$ mit $ a_{ij}=g_i(t_j)$, so kann man das Problem auch als $ Ax=b$ darstellen. Dieses lineare Gleichungssystem kann genau dann eindeutig gelöst werden, wenn der Rang von $ A$ voll ist, d.h. wenn er $ n$ beträgt.

Hat man dagegen mehr Beobachtungen gemacht als man Parameter hat, so wird das entstandene lineare Gleichungssystem wegen Messfehlern im Allgemeinen keine Lösung mehr besitzen. Liegen also $ m$ Messungen vor, wobei $ m>n$ sei, erhält man wieder ein Gleichungssystem der Form $ Ax=b$, aber mit $ A\in\mathbb{R}^{m,n}$ und $ b\in\mathbb{R}^m$. In diesem Falle kann man aber nach einer Lösung fragen, die die Messdaten in einem gewissen Sinne noch ganz gut trifft (salopp also $ Ax\approx b$ erfüllt). Als geeignetes Maß hat sich erwiesen, die Differenz von $ Ax$ und $ b$ im Sinne der euklidischen Norm auf dem $ \mathbb{R}^n$ möglichst klein zu bekommen. Der Differenzvektor $ Ax-b$ heißt Residuum. Das Ziel ist also eine Lösung zu folgendem Problem zu bestimmen:

minimiere $\displaystyle \vert\vert Ax-b\vert\vert _2$    über alle $\displaystyle x \in \mathbb{R}^n$
Da Quadrieren das Problem nicht verändert, betrachtet man:
minimiere $\displaystyle \vert\vert Ax-b\vert\vert _2^2$    über alle $\displaystyle x \in \mathbb{R}^n$
Dieses Problem nennt man das lineare Ausgleichsproblem. Die Lösung kann man mittels der QR-Zerlegung von $ A$ finden. Man stellt dabei $ A$ als Produkt einer orthogonalen Matrix $ Q\in\mathbb{R}^{m,m}$ und einer $ m\times n$-Matrix
\begin{displaymath}\left(% \begin{array}{c} R \ 0 \ \end{array}% \right)\end{displaymath}
dar, deren oberer Teil aus einer oberen Dreiecksmatrix $ R\in\mathbb{R}^{n,n}$ besteht:
\begin{displaymath}A=Q\left(% \begin{array}{c} R \ 0 \ \end{array}% \righ... ...TA=\left(% \begin{array}{c} R \ 0 \ \end{array}% \right)\end{displaymath}
Weil orthogonale Matrizen die euklidische Norm invariant lassen, gilt
\begin{displaymath}\vert\vert Ax-b\vert\vert _2^2=\vert\vert Q^T(Ax-b)\vert\vert... ...2=\vert\vert Rx-b'\vert\vert _2^2+\vert\vert b''\vert\vert _2^2\end{displaymath}
wenn man den Vektor $ Q^Tb\in\mathbb{R}^m$ mit $ b'\in\mathbb{R}^n$ und $ b''\in\mathbb{R}^{m-n}$ wie folgt partitioniert:
\begin{displaymath}Q^Tb=\left(% \begin{array}{c} b' \ b'' \ \end{array}% \right)\end{displaymath}
Daran sieht man, dass genau das $ x$ den Term $ \vert\vert Ax-b\vert\vert _2^2$ minimiert, welches $ Rx=b'$ löst. Das Ergebnis schreibt man auch als $ x=A^+b$, wobei $ A^+\in \mathbb{R}^{n,m}$ die sog. Pseudoinverse von $ A$ ist. Die minimale Residuumsnorm beträgt in diesem Fall $ \vert\vert b''\vert\vert _2^2$.

 

 

 Das nichtlineare Ausgleichsproblem

Hierbei geht es allgemein um die Minimierung einer Funktion $ g : \mathbb{R}^n\to \mathbb{R}$. Allerdings ist bekannt dass $ g$ die Form $ g(x)=\frac{1}{2}\sum_{i=1}^m(r_i(x))^2$ hat. Die Komponenten der Abbildung $ r(x)=(r_1(x),\hdots,r_m(x))^T$ heißen Residuen. Wieder sollen also die unbekannten $ x_1,\hdots,x_n$ so bestimmt werden, dass dieser Term möglichst klein ist. Dies ist im Allgemeinen nicht mehr direkt zu lösen, sondern man muss mit iterativen Methoden arbeiten, die ausgehend von einer Näherungslösung $ x_0$ eine Folge $ \{x_i\}_{i\geq 0}$ erzeugen, die hoffentlich gegen eine Lösung des Problems konvergiert.

 

 Das Gauß-Newton-Verfahren

Da man die Norm des Residuum möglichst klein bekommen möchte, liegt die Idee nahe ein Verfahren zu benutzen, das sich analog dem Newton-Verfahren zur Lösung der Gleichung $ r(x)=0$ auf sukzessive Linearisierungen von $ \vert\vert r(x)\vert\vert _2$ stützt:
Gegeben sei eine Annäherung $ x_c\in \mathbb{R}^n$ an die gesuchte Lösung $ x$. Mittels einer Taylor-Entwicklung erhält man $ r(x)\approx r(x_c)+J(x_c)(x-x_c)$ wobei $ J$ die Jacobi-Matrix von $ r$ ist. Statt $ \vert\vert r(x)\vert\vert _2$ minimiert man nun in jedem Schritt $ \vert\vert r(x_c)+J(x_c)(x-x_c)\vert\vert _2$. Dies ist ein lineares Ausgleichsproblem. Also erhält man folgendes Verfahren
  1. Wähle $ x_0$ und berechne $ r(x_0)$
  2. $ k=0$
  3. Bestimme $ x_{k+1}$, so dass $ \vert\vert r(x_k)+J(x_k)(x-x_k)\vert\vert _2$ bei $ x=x_{k+1}$ minimal wird.
  4. Berechne $ r(x_{k+1})$
  5. Wenn $ \vert r(x_k)-r(x_{k+1})\vert$ klein ist, dann STOP, sonst weiter bei Schritt 6
  6. Setze $ k=k+1$
  7. Weiter bei Schritt 3

In Schritt 3 wird $ x_{k+1}$ explizit durch $ x_{k+1}=x_k-J^+(x_k)r(x_k)$ berechnet, wobei $ J^+(x_k)$ die Pseudoinverse zu $ J(x_k)$ bezeichnet.

Bei dieser Vorgehensweise ist i.a. nicht gesichert, dass das Residuum $ \vert\vert r(x_{k+1})\vert\vert _2$ kleiner als $ \vert\vert r(x_k)\vert\vert _2$ ist. Die Idee ist nun einen Dämpfungsfaktor $ \lambda_k\leq 1$ einzuführen und in Schritt 3 $ x_{k+1}=x_k-\lambda_k J^+(x_k)r(x_k)$ zu rechnen, so dass auch wirklich $ \vert\vert r(x_{k+1})\vert\vert _2<\vert\vert r(x_k)\vert\vert _2$ erreicht wird. Für den Dämpfungsfaktor werden nacheinander die Werte $ 1, \frac{1}{2}, \frac{1}{4}, \hdots$ ausprobiert, bis sich eine Verbesserung der Residuumsnorm ergibt.

 

 Anpassung von Funktionen

Möchte man eine nichtlineare Funktion $ f(t,x)$ durch gegebene Messpunkte $ (t_i,y_i),\quad i=1,\hdots,m$ laufen lassen, so ist dies sicherlich nicht immer möglich. Aber man kann fordern, dass an jedem Punkt $ t_i$ die Differenz von $ f(t_i,x)$ und $ y_i$ möglichst klein sein soll.

Genauer: Es ist das $ x\in\mathbb{R}^n$ gesucht, für das der Term

$\displaystyle \sum_{i=1}^m\left[y_i-f(t_i,x)\right]^2$
minimal wird. Mit $ r_i(x):=y_i-f(t_i,x)$ ist dies also ein nichtlineares Ausgleichsproblem wie im vorherigen Abschnitt.

Das folgende Programm ermöglicht eine solche Anpassung von einer beliebigen Funktion an beliebige Daten. Als Methoden stehen sowohl das ungedämpfte als auch das gedämpfte Gauß-Newton-Verfahren zur Verfügung. Schafft das Programm keine Anpassung, steckt die Folge $ \{x_i\}_{i\geq 0}$ also z.B. in einem lokalen Minimum fest, so hilft es meist nur die Anfangsvorgabe zu verbessern.