diff --git a/programming/lecture/programming.tex b/programming/lecture/programming.tex index 975ebeb..2aa6219 100644 --- a/programming/lecture/programming.tex +++ b/programming/lecture/programming.tex @@ -119,7 +119,6 @@ wie eine einzelne Variable gel\"oscht wird. Der Datentyp bestimmt, wie die im Speicher abgelegten Bitmuster interpretiert werden. Die wichtigsten Datentpyen sind folgende: - \begin{itemize} \item \codetterm{integer} - Ganze Zahlen. Hier gibt es mehrere Unterarten, die wir in \matlab{} (meist) ignorieren k\"onnen. @@ -148,7 +147,7 @@ unterschiedlichem Speicherbedarf und Wertebreich. Daten gespeichert werden. Dennoch lohnt es sich, sich ein wenig mit den Datentypen auseinanderzusetzen (Box \ref{daqbox}). -\begin{ibox}[tp]{\label{daqbox}Digitalisierung von Messdaten} +\begin{ibox}[t]{\label{daqbox}Digitalisierung von Messdaten} Szenario: Die elektrische Aktivit\"at (z.B. die Membranspannung) einer Nervenzelle wird gemessen. Die gemessene Spannung wird mittels Messkarte digitalisiert und auf dem Computer f\"ur weitere Analysen diff --git a/regression/code/errorSurface.m b/regression/code/errorSurface.m index 636569f..ba3b7fe 100644 --- a/regression/code/errorSurface.m +++ b/regression/code/errorSurface.m @@ -1,6 +1,3 @@ -clear -close all - load('lin_regression.mat'); % compute mean squared error for a range of sloopes and intercepts: @@ -16,14 +13,10 @@ end % plot the error surface: figure() [N,M] = meshgrid(intercepts, slopes); -s = surface(M,N,error_surf); +s = surface(M, N, error_surf); xlabel('slope', 'rotation', 7.5) ylabel('intercept', 'rotation', -22.5) zlabel('error') set(gca,'xtick', (-5:2.5:5)) grid on view(3) - -set(gcf, 'paperunits', 'centimeters', 'papersize', [15, 15], ... - 'paperposition', [0., 0., 15, 15]) -saveas(gcf, 'error_surface', 'pdf') diff --git a/regression/code/lsqGradient.m b/regression/code/lsqGradient.m index c848b81..2736370 100644 --- a/regression/code/lsqGradient.m +++ b/regression/code/lsqGradient.m @@ -1,7 +1,17 @@ -function gradient = lsqGradient(parameter, x, y) -h = 1e-6; +function gradient = lsqGradient(x, y, parameter) +% The gradient of the least square error +% +% Arguments: x, the input values +% y, the corresponding measured output values +% parameter, vector containing slope and intercept +% as the 1st and 2nd element +% +% Returns: the gradient as a vector with two elements -partial_m = (lsqError([parameter(1)+h, parameter(2)],x,y) - lsqError(parameter,x,y))/ h; -partial_n = (lsqError([parameter(1), parameter(2)+h],x,y) - lsqError(parameter,x,y))/ h; + h = 1e-6; % stepsize for derivatives -gradient = [partial_m, partial_n]; \ No newline at end of file + partial_m = (lsqError(x, y, [parameter(1)+h, parameter(2)]) - lsqError(x, y, parameter))/ h; +partial_n = (lsqError(x, y, [parameter(1), parameter(2)+h]) - lsqError(x, y, parameter))/ h; + + gradient = [partial_m, partial_n]; +end diff --git a/regression/lecture/derivative.py b/regression/lecture/derivative.py new file mode 100644 index 0000000..7d85716 --- /dev/null +++ b/regression/lecture/derivative.py @@ -0,0 +1,53 @@ +import numpy as np +import matplotlib.pyplot as plt + +#plt.xkcd() +fig = plt.figure( figsize=(2.5,3.4) ) +ax = fig.add_subplot(1, 1, 1) + +# parabula: +x1 = -0.2 +x2 = 0.9 +x = np.linspace(x1, x2, 200) +y = x*x +ax.set_xlim(x1, x2) +ax.set_ylim(-0.2, 0.7) +ax.plot(x, y, c='b', lw=4, zorder=0) +# secant: +x = np.asarray([0.1, 0.7]) +y = x*x +ax.set_xticks(x) +ax.set_yticks(y) +ax.set_xticklabels(['$x$','$x+\Delta x$']) +ax.set_yticklabels(['','']) +ax.scatter(x, y, c='#CC0000', edgecolor='none', s=150, zorder=10) +# function values: +ax.plot([x[0], x[0], x1],[-0.2, y[0], y[0]], '--k', lw=1, zorder=6) +ax.plot([x[1], x[1], x1],[-0.2, y[1], y[1]], '--k', lw=1, zorder=6) +ax.text(x1+0.05, y[0]+0.05, '$f(x)$', zorder=6) +ax.text(x1+0.05, y[1]+0.05, '$f(x+\Delta x)$', zorder=6) +# slope triangle: +ax.plot([x[0], x[1], x[1]],[y[0], y[0], y[1]], '-k', lw=2, zorder=7) +ax.text(np.mean(x), y[0]-0.08, '$\Delta x$', ha='center', zorder=7) +ax.text(x[1]+0.05, np.mean(y), '$f(x+\Delta x)-f(x)$', va='center', rotation='vertical', zorder=7) +# secant line: +m = np.diff(y)/np.diff(x) +xl = [x1, x2] +yl = m*(xl-x[0])+y[0] +ax.plot(xl, yl, c='#CC0000', lw=3, zorder=7) + +# derivative: +md = 2.0*x[0] +yl = md*(xl-x[0])+y[0] +ax.plot(xl, yl, c='#FFFF66', lw=3, zorder=5) + +# limit: +for ml in np.linspace(md, m, 5)[1:] : + yl = ml*(xl-x[0])+y[0] + xs = 0.5*(ml+np.sqrt(ml*ml-4.0*(ml*x[0]-y[0]))) + ax.scatter([xs], [xs*xs], c='#FFCC00', edgecolor='none', s=80, zorder=3) + ax.plot(xl, yl, c='#FFCC00', lw=2, zorder=3) + +plt.tight_layout(); +plt.savefig('derivative.pdf') +#plt.show(); diff --git a/regression/lecture/error_surface.py b/regression/lecture/error_surface.py index 124f13d..30ca2f6 100644 --- a/regression/lecture/error_surface.py +++ b/regression/lecture/error_surface.py @@ -41,8 +41,7 @@ def plot_error_plane(ax, x, y, m, n): linewidth=0, shade=True) # Minimum: mini = np.unravel_index(np.argmin(error_surf), error_surf.shape) - #ax.scatter([m], [n], [0.0], color='#cc0000') - ax.scatter(slopes[mini[1]], intercepts[mini[0]], [0.0], color='#cc0000') + ax.scatter(slopes[mini[1]], intercepts[mini[0]], [0.0], color='#cc0000', s=60) if __name__ == "__main__": diff --git a/regression/lecture/gradient.py b/regression/lecture/gradient.py new file mode 100644 index 0000000..b762217 --- /dev/null +++ b/regression/lecture/gradient.py @@ -0,0 +1,42 @@ +import numpy as np +import matplotlib.pyplot as plt + +#plt.xkcd() +fig = plt.figure( figsize=(3.5,3.7) ) +ax = fig.add_subplot(1, 1, 1) + +def gaussian(x, y) : + return np.exp(-0.5*(x*x+y*y)) + +def gaussiangradient(x, y) : + return np.asarray([-x*np.exp(-0.5*(x*x+y*y)), -y*np.exp(-0.5*(x*x+y*y))]) + +# Gauss contour lines: +x1 = -2.0 +x2 = 2.0 +x = np.linspace(x1, x2, 200) +y = np.linspace(x1, x2, 200) +X, Y = np.meshgrid(x, y) +Z = gaussian(X, Y) +ax.set_xlabel('x') +ax.set_ylabel('y', rotation='horizontal') +ax.set_xticks(np.arange(x1, x2+0.1, 1.0)) +ax.set_yticks(np.arange(x1, x2+0.1, 1.0)) +ax.contour(X, Y, Z, zorder=0) + +# gradients: +xxg = [-1.1, 1.4, -1.0] +yyg = [-1.1, 0.6, 0.75] +for xg, yg in zip(xxg, yyg) : + zg = gaussian(xg, yg) + g = gaussiangradient(xg, yg) + ax.scatter(xg, yg, c='k', s=100, zorder=10) + ax.quiver([xg, xg], [yg, yg], [g[0], 0.0], [0.0, g[1]], units='xy', angles='uv', scale_units='x', scale=0.5, zorder=10) + ax.quiver([xg], [yg], [g[0]], [g[1]], units='xy', angles='uv', scale_units='x', scale=0.5, zorder=10, lw=2) +ax.text(-0.8, -1.5, '$\partial f/\partial x$', ha='center', zorder=20) +ax.text(-1.55, -0.8, '$\partial f/\partial y$', va='center', rotation='vertical', zorder=20) +ax.text(-0.4, -0.8, r'$\nabla f$', ha='left', zorder=20) + +plt.tight_layout(); +plt.savefig('gradient.pdf') +#plt.show(); diff --git a/regression/lecture/regression.tex b/regression/lecture/regression.tex index f0ac5d8..4985637 100644 --- a/regression/lecture/regression.tex +++ b/regression/lecture/regression.tex @@ -201,83 +201,124 @@ 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. +\begin{ibox}[t]{\label{differentialquotientbox}Differenzenquotient und Ableitung} + \includegraphics[width=0.33\textwidth]{derivative} + \hfill + \begin{minipage}[b]{0.63\textwidth} + Der Differenzenquotient + \[ m = \frac{f(x + \Delta x) - f(x)}{\Delta x} \] einer Funktion $y = f(x)$ ist + die Steigung der Sekante (rot) durch die beiden Punkte $(x,f(x))$ und + $(x+\Delta x,f(x+\Delta x))$ mit dem Abstand $\Delta x$. + + Die Steigung einer Funktion $y=f(x)$ an einer Stelle $x$ (gelb) wird durch + die Ableitung $f'(x)$ der Funktion an dieser Stelle berechnet. Die + Ableitung ist \"uber den Grenzwert (orange) des Differenzenquotienten f\"ur + unendlich kleine Abst\"ande $\Delta x$ definiert: + \[f'(x) = \frac{{\rm d} f(x)}{{\rm d}x} = \lim\limits_{\Delta x \to 0} \frac{f(x + \Delta x) - f(x)}{\Delta x} \] + \end{minipage} +\end{ibox} + +\begin{ibox}[t]{\label{partialderivativebox}Partielle Ableitungen und Gradient} + Bei Funktionen + \[ z = f(x,y) \] + die von mehreren Variablen, z.B. $x$ und $y$ abh\"angen, + kann die Steigung in Richtung jeder dieser Variablen + mit den partiellen Ableitungen + \[ \frac{\partial f(x,y)}{\partial x} = \lim\limits_{\Delta x \to 0} \frac{f(x + \Delta x,y) - f(x,y)}{\Delta x} \] + und + \[ \frac{\partial f(x,y)}{\partial y} = \lim\limits_{\Delta y \to 0} \frac{f(x, y + \Delta y) - f(x,y)}{\Delta y} \] + definiert \"uber den jeweiligen Differenzenquotienten + (Box~\ref{differentialquotientbox}) berechnet werden. \vspace{1ex} + + \begin{minipage}[t]{0.44\textwidth} + \mbox{}\\[-2ex] + \includegraphics[width=1\textwidth]{gradient} + \end{minipage} + \hfill + \begin{minipage}[t]{0.52\textwidth} + Z.B. lauten die partiellen Ableitungen von + \[ f(x,y) = x^2+y^2 \] + \[ \frac{\partial f(x,y)}{\partial x} = 2x \; , \quad \frac{\partial f(x,y)}{\partial y} = 2y \; .\] + + Der Gradient ist der aus den partiellen Ableitungen gebildete Vektor + \[ \nabla f(x,y) = \left( \begin{array}{c} \frac{\partial f(x,y)}{\partial x} \\[1ex] \frac{\partial f(x,y)}{\partial y} \end{array} \right) \] + und zeigt in Richtung des st\"arksten Anstiegs der Funktion $f(x,y)$. + \end{minipage} + + \vspace{1ex} Die Abbildung zeigt die Kontourlinien einer bivariaten + Gau{\ss}glocke $f(x,y) = \exp(-(x^2+y^2)/2)$ und den Gradienten mit + seinen partiellen Ableitungen an drei verschiedenen Stellen. +\end{ibox} + \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}\] +Wenn eine Kugel an einem beliebigen Startpunkt auf der Fehlerfl\"ache +\figref{errorsurfacefig} losgelassen werden w\"urde, dann w\"urde sie +entlang des steilsten Gef\"alles zum Minimum der Fehlerfl\"ache rollen +und dort zum Stehen kommen. Um den Weg der Kugel mit einem +Computerprogramm nachvollziehen zu k\"onnen, ben\"otigen wir +Information \"uber die Richtung des steilsten Gef\"alles an der +jeweils aktuellen Position. + +Der Gradient (Box~\ref{partialderivativebox}) der Kostenfunktion +\[ \nabla e(m,b) = \left( \begin{array}{c} \frac{\partial + e(m,b)}{\partial m} \\[1ex] \frac{\partial f(m,b)}{\partial + b} \end{array} \right) \] +bzgl. der beiden Parameter $m$ und $b$ der Geradengleichung ist ein +Vektor, der in Richtung des steilsten Anstiegs der Kostenfunktion +$e(m,b)$ zeigt. Die L\"ange des Gradienten gibt die St\"arke des +Anstiegs an (\figref{gradientquiverfig})). Da wir aber abw\"arts zum +Minimum laufen wollen, m\"ussen wir die dem Gradienten +entgegengesetzte Richtung einschlagen. -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. +Die partiellen Ableitungen m\"ussen nicht analytisch berechnet werden +sondern k\"onnen numerisch entsprechend dem Differenzenquotienten +(Box~\ref{differentialquotientbox}) mit kleinen Schrittweiten $\Delta +m$ und $\Delta b$ angen\"ahert werden: +\[\frac{\partial e(m,b)}{\partial m} = \lim\limits_{\Delta m \to + 0} \frac{e(m + \Delta m, b) - e(m,b)}{\Delta m} \approx \frac{e(m + \Delta m, b) - + e(m,b)}{\Delta m} \; . \] \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 + \includegraphics[width=0.75\columnwidth]{error_gradient} + \titlecaption{Der Gradienten der Fehlerfl\"ache.} + {Jeder Pfeil zeigt die Richtung und die Steigung f\"ur verschiedene Parameterkombination aus Steigung und - y-Achsenabschnitt an. Die Kontourlinien im Hintergrund + $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. +\begin{exercise}{lsqGradient.m}{}\label{gradientexercise}% + Implementiere eine Funktion \code{lsqGradient}, die die $x$- und + $y$-Werte der Messdaten sowie den Parametersatz $(m, b)$ der + Geradengleichung als 2-elementigen Vektor als Argumente + entgegennimmt und den Gradienten an dieser Stelle 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 + Benutzt die Funktion aus der vorherigen \"Ubung (\ref{gradientexercise}), + um f\"ur 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 +Zu guter Letzt muss nur noch der Gradientenabstieg implementiert +werden. Die daf\"ur ben\"otigten Zutaten haben wir aus den +vorangegangenen \"Ubungen bereits vorbereitet. 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)$. + b_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 @@ -287,14 +328,14 @@ f\"ur den Abstieg lautet: 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)\] + \[p_{i+1} = p_i - \epsilon \cdot \nabla e(m_i, b_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 +\"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}. @@ -318,3 +359,37 @@ Punkte in Abbildung \ref{gradientdescentfig} gro{\ss}. \item Erstelle einen Plot, der den besten Fit in die Daten plottet. \end{enumerate} \end{exercise} + + +\section{Fazit} +Mit dem Gradientenabstieg haben wir eine wichtige Methode zur +Bestimmung eines globalen Minimums einer Kostenfunktion +kennengelernt. F\"ur den Fall des Kurvenfits mit einer +Geradengleichung zeigt der mittlere quadratische Abstand als +Kostenfunktion in der Tat ein einziges klar definiertes Minimum. + +Wie wir im n\"achsten Kapitel sehen werden, kann die Position des +Minimums bei Geradengleichungen sogar analytisch bestimmt werden, der +Gradientenabstieg w\"are also gar nicht n\"otig +\matlabfun{polyfit}. Bei allen nichtlinearen Funktionen, wie z.B. der +Exponentialfunktion $f(x;\lambda) = \exp(\lambda x)$ gibt es aber +keine analytische L\"osung, und das Minimum der Kostenfunktion muss +numerisch, z.B. mit dem Gradientenabstiegsverfhren bestimmt werden. + +Um noch schneller das Minimum zu finden, kann der Gradientenabstieg +auf vielf\"altige Wei{\ss}e verfeinert werden. z.B. kann die +Schrittweite an die St\"arke des Gradienten angepasst werden. Diese +numerischen Tricks sind in bereits vorhandenen Funktionen +implementiert. Allgemeine Funktionen sind f\"ur beliebige +Kostenfunktionen gemacht \matlabfun{fminsearch}, w\"ahrend spezielle +Funktionen z.B. f\"ur die Minimierung des quadratischen Abstands bei +einem Kurvenfit angeboten werden \matlabfun{lsqcurvefit}. + +\begin{important} + Das Finden des globalen Minimums ist leider nur selten so leicht wie + bei einem Geradenfit. Oft gibt es viele Nebenminima, in denen der + Gradientenabstieg enden kann, obwohl das gesuchte globale Minimum + noch weit entfernt ist. Darum ist es meist sehr wichtig, wirklich + gute Startwerte f\"ur die zu bestimmenden Parameter der + Kostenfunktion zu haben. +\end{important}