scientificComputing/regression/lecture/regression.tex

\chapter{\tr{Optimization and Gradient Descent}{Optimierung und Gradientenabstieg}}

Ein sehr h\"aufiges Problem ist, dass die Abh\"angigkeit von
Messwerten von einer Eingangsgr\"o{\ss}e durch ein Modell erkl\"art
werden soll. Das Modell enth\"alt \"ublicherweise einen oder mehrere
Parameter, die den Zusammenhang modifizieren. Wie soll die beste
Parameterisierung des Modells gefunden werden, so dass das Modell die
Daten am besten beschreibt? Dieser Prozess der Parameteranpassung ist
ein Optimierungsproblem, der als Kurvenfit bekannt ist
(\enterm{curve fitting}).

\begin{figure}[t]
  \includegraphics[width=1\textwidth]{lin_regress}\hfill
  \titlecaption{.}{F\"ur eine Reihe von Eingangswerten $x$,
    z.B. Stimulusintensit\"aten, wurden die Antworten $y$ eines
    Systems gemessen (links). Der postulierte lineare Zusammenhang hat
    als freie Parameter die Steigung (mitte) und den
    $y$-Achsenabschnitt (rechts).}\label{linregressiondatafig}
\end{figure}

Die Punktewolke in \figref{linregressiondatafig} legt
zum Beispiel nahe, einen (verrauschten) linearen Zusammenhang zwischen
der Eingangsgr\"o{\ss}e $x$ (\enterm{input}) und der Systemantwort
$y$ (\enterm{output}) zu postulieren.
Wir nehmen also an, dass die Geradengleichung 
\[y = f(x; m, b) = m\cdot x + b \] 
ein gutes Modell f\"ur das zugrundeliegende System sein k\"onnte
(Abbildung \ref{linregressiondatafig}).  Die Geradengleichung hat die
beiden Parameter Steigung $m$ und $y$-Achsenabschnitt $b$ und es wird
die Kombination von $m$ und $b$ gesucht, die die Systemantwort am
besten vorhersagt.

In folgenden Kapitel werden wir anhand dieses Beispiels zeigen, welche
Methoden hinter einem Kurvenfit stecken, wie also numerisch die
optimale Kombination aus Steigung und $y$-Achsen\-abschnitt gefunden
werden kann.


\section{Mittlere quadratischen Abweichung}

Zuerst m\"u{\ss}en wir pr\"azisieren, was wir unter optimalen
Parametern verstehen. Es sollen die Werte der Parameter der
Geradengleichung sein, so dass die entsprechende Gerade am besten die
Daten beschreibt.  Was meinen wir damit? Jeder $y$-Wert der $N$
Datenpaare wird einen Abstand $y_i - y^{est}_i$ zu den durch das
Modell vorhergesagten Werten $y^{est}_i$ (\enterm{estimate}) an den
entsprechenden $x$-Werten haben. In unserem Beispiel mit der
Geradengleichung ist die Modellvorhersage $y^{est}_i=f(x_i;m,b)$
gegeben durch die Geradengleichung
(\figref{leastsquareerrorfig}). F\"ur den besten Fit sollten dieser
Abst\"ande m\"oglichst klein sein.

Wir k\"onnten z.B. fordern, die Summe $\sum_{i=1}^N y_i - y^{est}_i$
m\"oglichst klein zu machen. Das funktioniert aber nicht, da diese
Summe auch dann klein wird, wenn die H\"alfte der $y$-Daten weit
oberhalb der Geraden und die andere H\"alfte weit darunter liegt, da
sich diese positiven und negativen Werte gegenseitig zu Zahlen nahe
Null aufsummieren. Besser w\"are es auf jeden Fall, die Summe des
Betrags der Abst\"ande $\sum_{i=1}^N |y_i - y^{est}_i|$ zu betrachten. Ein
kleiner Wert der Summe kann dann nur erreicht werden, wenn die
Abst\"ande der Datenpunkte von der Kurve tats\"achlich klein sind,
unabh\"angig ob sie \"uber oder unter der Gerade liegen. Statt der
Summe k\"onnen wir genauso gut fordern, dass der \emph{mittlere} Abstand
\begin{equation}
  \label{meanabserror}
  e(\{(x_i, y_i)\}|\{y^{est}_i\}) = \frac{1}{N} \sum_{i=1}^N |y_i - y^{est}_i|
\end{equation}
der Menge der $N$ Datenpaare $(x_i, y_i)$ gegeben die Modellvorhersagen
$y_i^{est}$ klein sein soll.

Am h\"aufigsten wird jedoch bei einem Kurvenfit der \determ{mittlere
  quadratische Abstand} (\enterm{mean squared distance} oder
\enterm{mean squared error})
\begin{equation}
  \label{meansquarederror}
  e(\{(x_i, y_i)\}|\{y^{est}_i\}) = \frac{1}{N} \sum_{i=1}^N (y_i - y^{est}_i)^2
\end{equation}
verwendet (\figref{leastsquareerrorfig}). Wie beim Betrag sind die
quadratischen Abst\"ande immer positiv, unabh\"angig ob die Datenwerte
\"uber oder unter der Kurve liegen. Durch das Quadrat werden
zus\"atzlich gro{\ss}e Abst\"ande st\"arker gewichtet.
      
\begin{exercise}{meanSquareError.m}{}\label{mseexercise}%
  Schreibe eine Funktion \code{meanSquareError}, die die mittlere
  quadratische Abweichung zwischen einem Vektor mit den beobachteten
  Werten $y$ und einem Vektor mit den entsprechenden Vorhersagen
  $y^{est}$ berechnet.\newpage
\end{exercise}


\section{Zielfunktion}

$e(\{(x_i, y_i)\}|\{y^{est}_i\})$ ist eine sogenannte
\determ{Zielfunktion}, oder \determ{Kostenfunktion} (\enterm{objective
  function}, \enterm{cost function}), da wir die Modellvorhersage so
anpassen wollen, dass der mittlere quadratische Abstand, also die
Zielfunktion, minimiert wird. In
Kapitel~\ref{maximumlikelihoodchapter} werden wir sehen, dass die
Minimierung des mittleren quadratischen Abstands \"aquivalent zur
Maximierung der Wahrscheinlichkeit ist, dass die Daten aus der
Modellfunktion stammen, unter der Vorraussetzung, dass die Daten
um die Modellfunktion normalverteilt streuen.

\begin{figure}[t]
  \includegraphics[width=1\textwidth]{linear_least_squares}
  \titlecaption{Ermittlung des mittleren quadratischen Abstands.}
  {Der Abstand (\enterm{error}, orange) zwischen der Vorhersage (rote
    Gerade) und den Messdaten (blaue Punkte) wird f\"ur jeden
    gemessenen Datenpunkt ermittelt (links). Anschlie{\ss}end werden
    die Differenzen zwischen Messwerten und Vorhersage quadriert
    (\enterm{squared error}) und der Mittelwert berechnet (rechts).}
  \label{leastsquareerrorfig}
\end{figure}

Die Kostenfunktion mu{\ss} nicht immer der mittlere quadratische
Abstand sein. Je nach Problemstellung kann die Kostenfunktion eine
beliebige Funktion sein, die die Parameter eines Modells auf einen
Wert abbildet, der in irgendeiner Weise die Qualit\"at des Modells
quantifiziert. Ziel ist es dann, diejenigen Parameterwerte zu finden,
bei der die Kostenfunktion --- oder eben ``Zielfunktion'' ---
minimiert wird.
%%% Einfaches verbales Beispiel?

Wenn wir nun in unsere Gleichung \eqref{meansquarederror} f\"ur die
Modellvorhersage $y^{est}$ die Geradengleichung einsetzen, erhalten wir
f\"ur die Zielfunktion
\begin{eqnarray}
  e(\{(x_i, y_i)\}|m,b) & = & \frac{1}{N} \sum_{i=1}^N (y_i - f(x_i;m,b)^2 \label{msefunc} \\
  & = & \frac{1}{N} \sum_{i=1}^N (y_i - m x_i - b)^2 \label{mseline}
\end{eqnarray}
den mittleren quadratischen Abstand der Datenpaare $(x_i, y_i)$
gegeben die Parameterwerte $m$ und $b$ der Geradengleichung. Ziel des
Kurvenfits ist es, die Werte f\"ur $m$ und $b$ so zu optimieren, dass
der Fehler \eqnref{mseline} minimal wird.

\begin{exercise}{lsqError.m}{}
  Implementiere die Zielfunktion f\"ur die Optimierung mit der
  linearen Geradengleichung als Funktion \code{lsqError}.
  \begin{itemize}
  \item Die Funktion \"ubernimmt drei Argumente: Das erste ist ein
    Vektor mit den $x$-Werten, an denen gemessen wurde, und das zweite
    ein Vektor mit den zugeh\"origen $y$-Werten.  Das dritte Argument
    ist ein 2-elementiger Vektor, der die Parameter \code{m} und
    \code{b} enth\"alt.
  \item Die Funktion gibt als Ergebniss den Fehler als mittleren
    quadratischen Abstand \eqnref{mseline} zur\"uck.
  \item Die Funktion soll die Funktion \code{meanSquareError} der
    vorherigen \"Ubung benutzen.
  \end{itemize}
\end{exercise}


\section{Fehlerfl\"ache}

Die beiden Parameter $m$ und $b$ der Geradengleichung spannen eine
F\"ache auf. F\"ur jede Kombination aus $m$ und $b$ k\"onnen wir den
Wert der Zielfunktion, hier der mittlere quadratische Abstand
\eqnref{meansquarederror}, berechnen.  Es gibt also f\"ur jeden Punkt
in der sogenannten \emph{Fehlerfl\"ache} einen Fehler. In diesem
Beispiel eines 2-dimensionalen Problems (zwei freie Parameter) kann
die Fehlerfl\"ache graphisch durch einen 3-d ``surface-plot''
dargestellt werden. Dabei werden auf der $x$- und der $y$-Achse die
beiden Parameter und auf der $z$-Achse der Fehlerwert aufgetragen
(\figref{errorsurfacefig}).

\begin{figure}[t]
  \includegraphics[width=0.75\columnwidth]{error_surface.pdf}
  \caption{\textbf{Fehlerfl\"ache.} Die beiden freien Parameter
    unseres Modells $m$ und $b$ spannen die Grundfl\"ache des Plots
    auf. F\"ur jede Kombination von Steigung $m$ und
    $y$-Achsenabschnitt $b$ wird die errechnete Vorhersage des Modells
    mit den Messwerten verglichen und der Fehlerwert geplottet. Die
    sich ergebende Fehlerfl\"ache hat ein Minimum (roter Punkt) bei
    den Werten von $m$ und $b$, f\"ur die die Gerade die Daten am
    besten beschreibt.}\label{errorsurfacefig}
\end{figure}

\begin{exercise}{errorSurface.m}{}\label{errorsurfaceexercise}%
  Lade den Datensatz \textit{lin\_regression.mat} in den Workspace (20
  Datenpaare in den Vektoren \code{x} und \code{y}). Schreibe ein Skript
  \file{errorSurface.m}, dass den Fehler, berechnet als mittleren
  quadratischen Abstand zwischen den Daten und einer Geraden mit
  Steigung $m$ und $y$-Achsenabschnitt $b$, in Abh\"angigkeit von $m$
  und $b$ als surface plot darstellt (siehe Hilfe f\"ur die
  \code{surf} Funktion).
\end{exercise}

An der Fehlerfl\"ache kann direkt erkannt werden, bei welcher
Parameterkombination der Fehler minimal, beziehungsweise die
Parameterisierung optimal an die Daten angepasst ist. Wie kann die
Fehlerfunktion und die durch sie definierte Fehlerfl\"ache nun benutzt
werden, um den Optimierungsprozess zu leiten?

Die naheliegenste Variante ist, von der Fehlerfl\"ache einfach den Ort
des globalen Minimums zu bestimmen. Das ist im Allgemeinen jedoch zu
rechenintensiv, da f\"ur jede m\"ogliche Kombination der Parameter der
Fehler berechnet werden muss. Die Anzahl der n\"otigen Berechnungen
steigt exponentiell mit der Anzahl der Parameter (``Fluch der
Dimension''). Auch eine bessere Genauigkeit, mit der das Minimum
bestimmt werden soll erh\"oht die Anzahl der n\"otigen
Berechnungen. Wir suchen also ein Verfahren, dass das Minimum der
Kostenfunktion mit m\"oglichst wenigen Berechnungen findet.


\section{Gradient}

Man kann sich den Optimierungsprozess veranschaulichen wenn man sich
vorstellt, dass eine Parameterkombination einen Punkt auf der
Fehlerfl\"ache definiert. Wenn von diesem Punkt aus eine Kugel
losgelassen w\"urde, dann w\"urde sie, dem steilsten Gef\"alle
folgend, zum Minimum der Fehlerfl\"ache rollen und dort zum Stehen
kommen. Um dem Computer zu sagen, in welche Richtung die Position
ver\"andert werden soll muss also die Steigung an der aktellen
Position berechnet werden.

Die Steigung einer Funktion an einer Stelle ist die Ableitung der
Funktion an dieser Stelle, die durch den Differenzquotienten f\"ur
unendlich kleine Schritte $h$ bestimmt wird.
\[f'(x) = \lim\limits_{h \rightarrow 0} \frac{f(x + h) - f(x)}{h} \]
Bei unserem Fittingproblem h\"angt die Fehlerfunktion von zwei
Parametern ab und wir bestimmen die Steigung partiell f\"ur jeden
Parameter einzeln. Die partielle Ableitung nach $m$ sieht so aus:
\[\frac{\partial g(m,n)}{\partial m} = \lim\limits_{h \rightarrow 0} \frac{g(m + h, n) - g(m,n)}{h}\]
Da wir die Ver\"anderung des Fehlers in Abh\"angigkeit der beiden
Parameter bewerten, ist der Gradient an der Stelle $(m,n)$ ein
zweielementigen Vektor, der aus den partiellen Ableitungen nach $m$ und
nach $n$ besteht. Die Richtung des Gradienten zeigt die Richtung der
gr\"o{\ss}ten Steigung an, seine L\"ange repr\"asentiert die St\"arke
des Gef\"alles.

\[\bigtriangledown g(m,n) = \left( \frac{\partial g(m,n)}{\partial m}, \frac{\partial g(m,n)}{\partial n}\right)\]
Die partiellen Ableitungen m\"ussen nicht analytisch berechnet werden
sondern wir n\"ahern sie numerisch durch einen sehr kleinen Schritt
an.
\[\frac{\partial g(m,n)}{\partial m} = \lim\limits_{h \rightarrow
  0} \frac{g(m + h, n) - g(m,n)}{h} \approx \frac{g(m + h, n) -
  g(m,n)}{h}\]

Die graphische Dargestellung der Gradienten (Abbildung
\ref{gradientquiverfig}) zeigt, dass diese in die Richtung der
gr\"o{\ss}ten Steigung zeigen. Um zum Minimum der Fehlerfunktion zu
gelangen sollte man also die entgegengesetzte Richtung einschlagen.

\begin{figure}
  \includegraphics[width=0.75\columnwidth]{figures/error_gradient}
  \caption{\textbf{Darstellung der Gradienten f\"ur jede
      Parameterkombination.} Jeder Pfeil zeigt die Richtung und die
    Steigung f\"ur verschiedene Parameterkombination aus Steigung und
    y-Achsenabschnitt an. Die Kontourlinien im Hintergrund
    illustrieren die Fehlerfl\"ache. Warme Farben stehen f\"ur
    gro{\ss}e Fehlerwerte, kalte Farben f\"ur kleine. Jede
    Kontourlinie steht f\"ur eine Linie gleichen
    Fehlers.}\label{gradientquiverfig}
\end{figure}

\begin{exercise}{lsqGradient.m}{}\label{gradientexercise}
  Implementiert eine Funktion \code{lsqGradient}, die den aktuellen
  Parametersatz als 2-elementigen Vektor, die x-Werte und die y-Werte
  als Argumente entgegennimmt und den Gradienten zur\"uckgibt.
\end{exercise}

\begin{exercise}{errorGradient.m}{}
  Benutzt die Funktion aus vorheriger \"Ubung (\ref{gradientexercise}),
  um f\"ur die jede Parameterkombination aus der Fehlerfl\"ache
  (\"Ubung \ref{errorsurfaceexercise}) auch den Gradienten zu
  berechnen und darzustellen. Vektoren im Raum k\"onnen mithilfe der
  Funktion \code{quiver} geplottet werden.
\end{exercise}

\section{Gradientenabstieg}

Zu guter Letzt muss ``nur'' noch der Gradientenabstieg implementiert
werden. Die daf\"ur ben\"otigten Zutaten sollten wir aus den
vorangegangenen \"Ubungen haben. Wir brauchen: 1. Die Fehlerfunktion
(\code{meanSquareError.m}), 2. die Zielfunktion (\code{lsqError.m})
und 3. den Gradienten (\code{lsqGradient.m}).  Der Algorithmus
f\"ur den Abstieg lautet:
\begin{enumerate}
\item Starte mit einer beliebigen Parameterkombination $p_0 = (m_0,
  n_0)$.
\item \label{computegradient} Berechne den Gradienten an der akutellen Position $p_i$.
\item Wenn die L\"ange des Gradienten einen bestimmten Wert
  unterschreitet, haben wir das Minum gefunden und k\"onnen die Suche
  abbrechen.  Wir suchen ja das Minimum, bei dem der Gradient gleich
  Null ist. Da aus numerischen Gr\"unden der Gradient nie exakt Null
  werden wird, k\"onnen wir nur fordern, dass er hinreichend klein
  wird (z.B. \code{norm(gradient) < 0.1}).
\item \label{gradientstep} Gehe einen kleinen Schritt ($\epsilon =
  0.01$) in die entgegensetzte Richtung des Gradienten:
  \[p_{i+1} = p_i - \epsilon \cdot \bigtriangledown g(m_i, n_i)\]
\item Wiederhole die Schritte \ref{computegradient} -- \ref{gradientstep}.
\end{enumerate}

Abbildung \ref{gradientdescentfig} zeigt den Verlauf des
Gradientenabstiegs. Von einer Startposition aus wird die Position
solange ver\"andert, wie der Gradient eine bestimmte Gr\"o{\ss}e
\"uberschreitet. An den Stellen, an denen der Gradient sehr stark ist
ist auch die Ver\"anderung der Position gro{\ss} und der Abstand der
Punkte in Abbildung \ref{gradientdescentfig} gro{\ss}.

\begin{figure}[hb]
  \includegraphics[width=0.6\columnwidth]{figures/gradient_descent}
  \caption{\textbf{Gradientenabstieg.} Es wird von einer beliebigen
    Position aus gestartet und der Gradient berechnet und die Position
    ver\"andert. Jeder Punkt zeigt die Position nach jedem
    Optimierungsschritt an.} \label{gradientdescentfig}
\end{figure}

\begin{exercise}{gradientDescent.m}{}
  Implementiere den Gradientenabstieg f\"ur das Problem der
  Parameteranpassung der linearen Geradengleichung an die Messdaten in
  der Datei \code{lin\_regression.mat}.
  \begin{enumerate}
  \item Merke Dir f\"ur jeden Schritt den Fehler zwischen
    Modellvorhersage und Daten.
  \item Erstelle eine Plot, der die Entwicklung des Fehlers als
    Funktion der Optimierungsschritte zeigt.
  \item Erstelle einen Plot, der den besten Fit in die Daten plottet.
  \end{enumerate}
\end{exercise}