\chapter{Optimization and gradient descent}
\exercisechapter{Optimization and gradient descent}

Optimization problems arise in many different contexts. For example,
to understand the behavior of a given system, the system is probed
with a range of input signals and then the resulting responses are
measured. This input-output relation can be described by a model. Such
a model can be a simple function that maps the input signals to
corresponding responses, it can be a filter, or a system of
differential equations. In any case, the model has same parameter that
specify how input and output signals are related.  Which combination
of parameter values are best suited to describe the input-output
relation? The process of finding the best parameter values is an
optimization problem. For a simple parameterized function that maps
input to output values, this is the special case of a \enterm{curve
  fitting} problem, where the average distance between the curve and
the response values is minimized. One basic numerical method used for
such optimization problems is the so called gradient descent, which is
introduced in this chapter.

%%% Weiteres einfaches verbales Beispiel? Eventuell aus der Populationsoekologie?

\begin{figure}[t]
  \includegraphics[width=1\textwidth]{lin_regress}\hfill
  \titlecaption{Example data suggesting a linear relation.}{A set of
    input signals $x$, e.g. stimulus intensities, were used to probe a
    system. The system's output $y$ to the inputs are noted
    (left). Assuming a linear relation between $x$ and $y$ leaves us
    with 2 parameters, the slope (center) and the intercept with the
    y-axis (right panel).}\label{linregressiondatafig}
\end{figure}

The data plotted in \figref{linregressiondatafig} suggest a linear
relation between input and output of the system. We thus assume that a
straight line
\begin{equation}
  \label{straightline}
  y = f(x; m, b) = m\cdot x + b
\end{equation}
is an appropriate model to describe the system.  The line has two free
parameter, the slope $m$ and the $y$-intercept $b$. We need to find
values for the slope and the intercept that best describe the measured
data.  In this chapter we use this example to illustrate the gradient
descent and how this methods can be used to find a combination of
slope and intercept that best describes the system.


\section{The error function --- mean squared error}

Before the optimization can be done we need to specify what exactly is
considered an optimal fit. In our example we search the parameter
combination that describe the relation of $x$ and $y$ best. What is
meant by this? Each input $x_i$ leads to an measured output $y_i$ and
for each $x_i$ there is a \emph{prediction} or \emph{estimation}
$y^{est}(x_i)$ of the output value by the model. At each $x_i$
estimation and measurement have a distance or error $y_i -
y^{est}(x_i)$. In our example the estimation is given by the equation
$y^{est}(x_i) = f(x_i;m,b)$. The best fitting model with parameters
$m$ and $b$ is the one that minimizes the distances between
observation $y_i$ and estimation $y^{est}(x_i)$
(\figref{leastsquareerrorfig}).

As a first guess we could simply minimize the sum $\sum_{i=1}^N y_i -
y^{est}(x_i)$. This approach, however, will not work since a minimal sum
can also be achieved if half of the measurements is above and the
other half below the predicted line. Positive and negative errors
would cancel out and then sum up to values close to zero. A better
approach is to sum over the absolute values of the distances:
$\sum_{i=1}^N |y_i - y^{est}(x_i)|$. This sum can only be small if all
deviations are indeed small no matter if they are above or below the
predicted line. Instead of the sum we could also take the average
\begin{equation}
  \label{meanabserror}
  f_{dist}(\{(x_i, y_i)\}|\{y^{est}(x_i)\}) = \frac{1}{N} \sum_{i=1}^N |y_i - y^{est}(x_i)|
\end{equation}
Instead of the averaged absolute errors, the \enterm[mean squared
error]{mean squared error} (\determ[quadratischer
Fehler!mittlerer]{mittlerer quadratischer Fehler})
\begin{equation}
  \label{meansquarederror}
  f_{mse}(\{(x_i, y_i)\}|\{y^{est}(x_i)\}) = \frac{1}{N} \sum_{i=1}^N (y_i - y^{est}(x_i))^2
\end{equation}
is commonly used (\figref{leastsquareerrorfig}). Similar to the
absolute distance, the square of the errors, $(y_i - y^{est}(x_i))^2$, is
always positive and thus positive and negative error values do not
cancel each other out. In addition, the square punishes large
deviations over small deviations. In
chapter~\ref{maximumlikelihoodchapter} we show that minimizing the
mean square error is equivalent to maximizing the likelihood that the
observations originate from the model, if the data are normally
distributed around the model prediction.

\begin{exercise}{meanSquaredErrorLine.m}{}\label{mseexercise}%
  Given a vector of observations \varcode{y} and a vector with the
  corresponding predictions \varcode{y\_est}, compute the \emph{mean
    square error} between \varcode{y} and \varcode{y\_est} in a single
  line of code.
\end{exercise}


\section{Objective function}

The mean squared error is a so called \enterm{objective function} or
\enterm{cost function} (\determ{Kostenfunktion}). A cost function
assigns to a model prediction $\{y^{est}(x_i)\}$ for a given data set
$\{(x_i, y_i)\}$ a single scalar value that we want to minimize.  Here
we aim to adapt the model parameters to minimize the mean squared
error \eqref{meansquarederror}. In general, the \enterm{cost function}
can be any function that describes the quality of the fit by mapping
the data and the predictions to a single scalar value.

\begin{figure}[t]
  \includegraphics[width=1\textwidth]{linear_least_squares}
  \titlecaption{Estimating the \emph{mean square error}.}  {The
    deviation error, orange) between the prediction (red
    line) and the observations (blue dots) is calculated for each data
    point (left). Then the deviations are squared and the aveage is
    calculated (right).}
  \label{leastsquareerrorfig}
\end{figure}

Replacing $y^{est}$ in the mean squared error \eqref{meansquarederror}
with our model, the straight line \eqref{straightline}, the cost
function reads
\begin{eqnarray}
  f_{cost}(m,b|\{(x_i, y_i)\}) & = & \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}
The optimization process tries to find the slope $m$ and the intercept
$b$ such that the cost function is minimized. With the mean squared
error as the cost function this optimization process is also called
method of the \enterm{least square error} (\determ[quadratischer
Fehler!kleinster]{Methode der kleinsten Quadrate}).

\begin{exercise}{meanSquaredError.m}{}
  Implement the objective function \eqref{mseline} as a function
  \varcode{meanSquaredError()}.  The function takes three
  arguments. The first is a vector of $x$-values and the second
  contains the measurements $y$ for each value of $x$. The third
  argument is a 2-element vector that contains the values of
  parameters \varcode{m} and \varcode{b}. The function returns the
  mean square error.
\end{exercise}


\section{Error surface}
For each combination of the two parameters $m$ and $b$ of the model we
can use \eqnref{mseline} to calculate the corresponding value of the
cost function. The cost function $f_{cost}(m,b|\{(x_i, y_i)\}|)$ is a
function $f_{cost}(m,b)$, that maps the parameter values $m$ and $b$
to a scalar error value.  The error values describe a landscape over the
$m$-$b$ plane, the error surface, that can be illustrated graphically
using a 3-d surface-plot. $m$ and $b$ are plotted on the $x$- and $y$-
axis while the third dimension indicates the error value
(\figref{errorsurfacefig}).

\begin{figure}[t]
  \includegraphics[width=0.75\textwidth]{error_surface}
  \titlecaption{Error surface.}{The two model parameters $m$ and $b$
    define the base area of the surface plot. For each parameter
    combination of slope and intercept the error is calculated. The
    resulting surface has a minimum which indicates the parameter
    combination that best fits the data.}\label{errorsurfacefig}
\end{figure}

\begin{exercise}{errorSurface.m}{}\label{errorsurfaceexercise}
  Generate 20 data pairs $(x_i|y_i)$ that are linearly related with
  slope $m=0.75$ and intercept $b=-40$, using \varcode{rand()} for
  drawing $x$ values between 0 and 120 and \varcode{randn()} for
  jittering the $y$ values with a standard deviation of 15.  Then
  calculate the mean squared error between the data and straight lines
  for a range of slopes and intercepts using the
  \varcode{meanSquaredError()} function from the previous exercise.
  Illustrates the error surface using the \code{surface()} function.
  Consult the documentation to find out how to use \code{surface()}.
\end{exercise}

By looking at the error surface we can directly see the position of
the minimum and thus estimate the optimal parameter combination. How
can we use the error surface to guide an automatic optimization
process?

The obvious approach would be to calculate the error surface for any
combination of slope and intercept values and then find the position
of the minimum using the \code{min} function. This approach, however
has several disadvantages: (i) it is computationally very expensive to
calculate the error for each parameter combination. The number of
combinations increases exponentially with the number of free
parameters (also known as the ``curse of dimensionality''). (ii) the
accuracy with which the best parameters can be estimated is limited by
the resolution used to sample the parameter space. The coarser the
parameters are sampled the less precise is the obtained position of
the minimum.

So we need a different approach. We want a procedure that finds the
minimum of the cost function with a minimal number of computations and
to arbitrary precision.

\begin{ibox}[t]{\label{differentialquotientbox}Difference quotient and derivative}
  \includegraphics[width=0.33\textwidth]{derivative}
  \hfill
  \begin{minipage}[b]{0.63\textwidth}
    The difference quotient
    \begin{equation}
      \label{difffrac}
      m = \frac{f(x + \Delta x) - f(x)}{\Delta x}
    \end{equation}
    of a function $y = f(x)$ is the slope of the secant (red) defined
    by the points $(x,f(x))$ and $(x+\Delta x,f(x+\Delta x))$ with the
    distance $\Delta x$.

    The slope of the function $y=f(x)$ at the position $x$ (yellow) is
    given by the derivative $f'(x)$ of the function at that position.
    It is defined by the difference quotient in the limit of
    infinitesimally (orange) small distances $\Delta x$:
    \begin{equation}
      \label{derivative}
      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{equation}
  \end{minipage}\vspace{2ex} 

  It is not possible to calculate the derivative, \eqnref{derivative},
  numerically. The derivative can only be estimated using the
  difference quotient, \eqnref{difffrac}, by using sufficiently small
  $\Delta x$.
\end{ibox}

\begin{ibox}[tp]{\label{partialderivativebox}Partial derivative and gradient}
  Some functions that depend on more than a single variable:
  \[ z = f(x,y) \]
  for example depends on $x$ and $y$. Using the partial derivative
  \[ \frac{\partial f(x,y)}{\partial x} = \lim\limits_{\Delta x \to 0} \frac{f(x + \Delta x,y) - f(x,y)}{\Delta x} \]
  and
  \[ \frac{\partial f(x,y)}{\partial y} = \lim\limits_{\Delta y \to 0} \frac{f(x, y + \Delta y) - f(x,y)}{\Delta y} \]
  one can estimate the slope in the direction of the variables
  individually by using the respective difference quotient
  (Box~\ref{differentialquotientbox}).  \vspace{1ex}

  \begin{minipage}[t]{0.5\textwidth}
    \mbox{}\\[-2ex]
    \includegraphics[width=1\textwidth]{gradient}
  \end{minipage}
  \hfill
  \begin{minipage}[t]{0.46\textwidth}
    For example, the partial derivatives of
    \[ f(x,y) = x^2+y^2 \] are 
    \[ \frac{\partial f(x,y)}{\partial x} = 2x \; , \quad \frac{\partial f(x,y)}{\partial y} = 2y \; .\]

    The gradient is a vector that is constructed from the partial derivatives:
    \[ \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) \]
    This vector points into the direction of the strongest ascend of
    $f(x,y)$.
  \end{minipage}

  \vspace{0.5ex} The figure shows the contour lines of a bi-variate
  Gaussian $f(x,y) = \exp(-(x^2+y^2)/2)$ and the gradient (thick
  arrows) and the corresponding two partial derivatives (thin arrows)
  for three different locations.
\end{ibox}


\section{Gradient}
Imagine to place a small ball at some point on the error surface
\figref{errorsurfacefig}. Naturally, it would roll down the steepest
slope and eventually stop at the minimum of the error surface (if it had no
inertia). We will use this picture to develop an algorithm to find our
way to the minimum of the objective function. The ball will always
follow the steepest slope. Thus we need to figure out the direction of
the steepest slope at the position of the ball.

The \entermde{Gradient}{gradient} (Box~\ref{partialderivativebox}) of the
objective function is the vector
\begin{equation}
  \label{gradient}
  \nabla f_{cost}(m,b) = \left( \frac{\partial f(m,b)}{\partial m},
    \frac{\partial f(m,b)}{\partial b} \right)
\end{equation}
that points to the strongest ascend of the objective function.  The
gradient is given by partial derivatives
(Box~\ref{partialderivativebox}) of the mean squared error with
respect to the parameters $m$ and $b$ of the straight line. There is
no need to calculate it analytically because it can be estimated from
the partial derivatives using the difference quotient
(Box~\ref{differentialquotientbox}) for small steps $\Delta m$ and
$\Delta b$. For example, the partial derivative with respect to $m$
can be computed as
\[\frac{\partial f_{cost}(m,b)}{\partial m} = \lim\limits_{\Delta m \to
  0} \frac{f_{cost}(m + \Delta m, b) - f_{cost}(m,b)}{\Delta m}
\approx \frac{f_{cost}(m + \Delta m, b) - f_{cost}(m,b)}{\Delta m} \;
. \]
The length of the gradient indicates the steepness of the slope
(\figref{gradientquiverfig}). Since want to go down the hill, we
choose the opposite direction.


\begin{figure}[t]
  \includegraphics[width=0.75\textwidth]{error_gradient}
  \titlecaption{Gradient of the error surface.}  {Each arrow points
    into the direction of the greatest ascend at different positions
    of the error surface shown in \figref{errorsurfacefig}. The
    contour lines in the background illustrate the error surface. Warm
    colors indicate high errors, colder colors low error values. Each
    contour line connects points of equal
    error.}\label{gradientquiverfig}
\end{figure}

\begin{exercise}{meanSquaredGradient.m}{}\label{gradientexercise}%
  Implement a function \varcode{meanSquaredGradient()}, that takes the
  $x$- and $y$-data and the set of parameters $(m, b)$ of a straight
  line as a two-element vector as input arguments. The function should
  return the gradient at the position $(m, b)$ as a vector with two
  elements.
\end{exercise}

\begin{exercise}{errorGradient.m}{}
  Extend the script of exercises~\ref{errorsurfaceexercise} to plot
  both the error surface and gradients using the
  \varcode{meanSquaredGradient()} function from
  exercise~\ref{gradientexercise}. Vectors in space can be easily
  plotted using the function \code{quiver()}. Use \code{contour()}
  instead of \code{surface()} to plot the error surface.
\end{exercise}


\section{Gradient descent}
Finally, we are able to implement the optimization itself. By now it
should be obvious why it is called the gradient descent method. All
ingredients are already there. We need: (i) the cost function
(\varcode{meanSquaredError()}), and (ii) the gradient
(\varcode{meanSquaredGradient()}). The algorithm of the gradient
descent works as follows:
\begin{enumerate}
\item Start with some given combination of the parameters $m$ and $b$
  ($p_0 = (m_0, b_0)$).
\item \label{computegradient} Calculate the gradient at the current
  position $p_i$.
\item If the length of the gradient falls below a certain value, we
  assume to have reached the minimum and stop the search. We are
  actually looking for the point at which the length of the gradient
  is zero, but finding zero is impossible because of numerical
  imprecision. We thus apply a threshold below which we are
  sufficiently close to zero (e.g. \varcode{norm(gradient) < 0.1}).
\item \label{gradientstep} If the length of the gradient exceeds the
  threshold we take a small step into the opposite direction:
  \[p_{i+1} = p_i - \epsilon \cdot \nabla f_{cost}(m_i, b_i)\]
  where $\epsilon = 0.01$ is a factor linking the gradient to
  appropriate steps in the parameter space.
\item Repeat steps \ref{computegradient} -- \ref{gradientstep}.
\end{enumerate}

\Figref{gradientdescentfig} illustrates the gradient descent --- the
path the imaginary ball has chosen to reach the minimum. Starting at
an arbitrary position on the error surface we change the position as
long as the gradient at that position is larger than a certain
threshold. If the slope is very steep, the change in the position (the
distance between the red dots in \figref{gradientdescentfig}) is
large.

\begin{figure}[t]
  \includegraphics[width=0.45\textwidth]{gradient_descent}
  \titlecaption{Gradient descent.}{The algorithm starts at an
    arbitrary position. At each point the gradient is estimated and
    the position is updated as long as the length of the gradient is
    sufficiently large.The dots show the positions after each
    iteration of the algorithm.} \label{gradientdescentfig}
\end{figure}

\begin{exercise}{gradientDescent.m}{}
  Implement the gradient descent for the problem of fitting a straight
  line to some measured data. Reuse the data generated in
  exercise~\ref{errorsurfaceexercise}.
  \begin{enumerate}
  \item Store for each iteration the error value.
  \item Plot the error values as a function of the iterations, the
    number of optimization steps.
  \item Plot the measured data together with the best fitting straight line.
  \end{enumerate}\vspace{-4.5ex}
\end{exercise}


\section{Summary}

The gradient descent is an important numerical method for solving
optimization problems. It is used to find the global minimum of an
objective function.

Curve fitting is a common application for the gradient descent method.
For the case of fitting straight lines to data pairs, the error
surface (using the mean squared error) has exactly one clearly defined
global minimum. In fact, the position of the minimum can be analytically
calculated as shown in the next chapter.

Problems that involve nonlinear computations on parameters, e.g. the
rate $\lambda$ in an exponential function $f(x;\lambda) = e^{\lambda
  x}$, do not have an analytical solution for the least squares. To
find the least squares for such functions numerical methods such as
the gradient descent have to be applied.

The suggested gradient descent algorithm can be improved in multiple
ways to converge faster.  For example one could adapt the step size to
the length of the gradient. These numerical tricks have already been
implemented in pre-defined functions. Generic optimization functions
such as \matlabfun{fminsearch()} have been implemented for arbitrary
objective functions, while the more specialized function
\matlabfun{lsqcurvefit()} i specifically designed for optimizations in
the least square error sense.

%\newpage
\begin{important}[Beware of secondary minima!]
  Finding the absolute minimum is not always as easy as in the case of
  fitting a straight line. Often, the error surface has secondary or
  local minima in which the gradient descent stops even though there
  is a more optimal solution, a global minimum that is lower than the
  local minimum. Starting from good initial positions is a good
  approach to avoid getting stuck in local minima. Also keep in mind
  that error surfaces tend to be simpler (less local minima) the fewer
  parameters are fitted from the data. Each additional parameter
  increases complexity and is computationally more expensive.
\end{important}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\printsolutions