% Except otherwise noted,
% this work is © 2009 Dmitri Fedorov, under CC-BY-SA license, version 3+.

\chapter{Minimization and optimization}
\section{Introduction}
{\em Minimization} is a problem of finding the minimum (or the maximum)
of a given, generally non-linear, real valued function $f(\ve{x})$ (often called the {\em objective function})
of an $n$-dimensional argument $\ve{x}\doteq\{x_1,\dots,x_n\}$.

Minimization is a simple case of a more general poblem --- {\em
optimization} --- which includes finding best available values of some
objective function within a given domain and subject to given constrains.

Minimization is related to root-finding as at the minimum all partial
derivatives of the objective function vanish,
	\begin{equation}\label{}
\frac{\partial f}{\partial x_i} = 0 \:.
	\end{equation}
One can alternatively solve this system of non-linear equations.

\section{Local minimization}
\subsection{Newton's methods}
Newton's method is based on the quadratic approxiamtion of the
objective function $f$ in the vicinity of the suspected minimum,
	\begin{equation}\label{}
f(\ve{x}+\Delta\ve{x}) \approx f(\ve{x})+\nabla
f(\ve{x})^T\Delta\ve{x} +
\frac12 \Delta\ve{x}^T \ma{H}(\ve{x}) \Delta\ve{x}
\:,
	\end{equation}
where the vector $\nabla f(\ve{x})$ is the gradient of the objective
function at the point $\ve{x}$,
	\begin{equation}\label{}
\nabla f(\ve{x}) \doteq
\left[ \frac{\partial f(\ve{x})}{\partial x_i} \right]_{i=1,\dots,n} \:,
	\end{equation}
and $\ma{H}(\ve{x})$ is the {\em Hessian matrix}\index{Hessian matrix} --
a square matrix made of second-order partial derivatives of the objective
function at the point $\ve{x}$,
	\begin{equation}\label{}
\ma{H}(\ve{x}) \doteq \left[ \frac{\partial^2 f(\ve{x})}{\partial
x_i\partial x_j} \right]_{i,j\in 1,\dots,n} \:.
	\end{equation}
The minimum of this quadratic form as function of $\Delta\ve{x}$ is found
at the point where its gradient with respect to $\Delta\ve{x}$ vanishes,
	\begin{equation}\label{o:eq:nf}
\nabla f(\ve{x}) + \ma{H}(\ve{x}) \Delta\ve{x} = 0 \:.
	\end{equation}
This gives an approximate step towards the minimum, called the {\em
Newton's step}\index{Newton's step},
	\begin{equation}\label{}
\Delta\ve{x} = -\ma{H}(\ve{x})^{-1}\nabla f(\ve{x}) \:.
	\end{equation}

The original Newton's method is simply the iteration,
	\begin{equation}\label{}
\ve{x}_{k+1} = \ve{x}_{k} - \ma{H}(\ve{x}_k)^{-1}\nabla f(\ve{x}_k) \:,
	\end{equation}
where at each iteration the full Newton's step is taken and the Hessian
matrix is recalculated.  In practice, instead of calculating $\ma{H}^{-1}$
one rather solves the linear equation~(\ref{o:eq:nf}).

Usually Newton's method is modified to include a smaller step size,
$\lambda\Delta\ve{x}$, with $0<\lambda<1$.  The step-size $\lambda$
can be found by a backtracking algoritm similar to that in the Newton's
method for root-finding. One starts with $\lambda=1$ and than
backtracks, $\lambda\leftarrow\lambda/2$,
until the {\em Armijo condition},
	\begin{equation}\label{}
f(\ve{x}+\lambda\Delta\ve{x}) < f(\ve{x}) +
\alpha \nabla f(\ve{x})\lambda\Delta\ve{x} \:,
	\end{equation}
is satisfied.  The parameter $\alpha$ can be chosen as small as
$10^{-4}$.

\subsection{Quasi-Newton methods}
Quasi-Newton methods are variations of the Newton's method which attempt
to avoid recalculation of the Hessian matrix at each iteration, trying
instead certain updates based on the analysis of the gradient vectors.
The update $\delta \ma{H}$ is usually chosen to satisfy the condition
	\begin{equation}\label{o:eq:up}
\nabla f(\ve{x}+\ve{s}) = \nabla f(\ve{x}) + (\ma{H}+\delta
\ma{H})\ve{s} \:,
	\end{equation}
called {\em secant equation}, which is the Taylor expansion of the
gradient.  This secant condition is under-determined in more than one
dimension as it consists of only $n$ equations for the $n^2$ unknown
elements of the update $\delta \ma{H}$.  Various quasi-Newton methods use
different choices for the form of the solution of the secant equation.

In quasi-Newton methods one often starts with a unity matrix as the
zeroth approximation for the Hessian matrix and then applies the updates.
In this case one can actually use the inverse Hessian matrix---thus
avoiding the need to solve the linear equation at each iteration---and
apply the updates directly to the inverse matrix.

\subsubsection{Broydens's update}
The Broyden's update is chosen in the form
	\begin{equation}\label{}
\delta \ma{H} = \ve{c}\ve{s}^T \;.
	\end{equation}
where the vector $\ve{c}$ is found from the condition~(\ref{o:eq:up}).
The resulting update is
	\begin{equation}\label{}
\ma{H} \rightarrow \ma{H}+
\frac{\ve{y}-H\ve{s}}
{\ve{s}^T\ve{s}}\ve{s}^T \:,
	\end{equation}
where $\ve{y} \doteq \nabla f(\ve{x}+\ve{s})-\nabla f(\ve{x})$.
The update for the inverse matrix is given as\footnote{using the
Sherman-Morrison formula,
$$
(A+uv^T)^{-1} = A^{-1} - \frac{A^{-1}uv^T A^{-1} }{ 1 + v^T A^{-1}u}.
$$
}
	\begin{equation}\label{}
\ma{H}^{-1} \rightarrow \ma{H}^{-1}+
\frac{(\ve{s}-\ma{H}^{-1}\ve{y})\ve{s}^T\ma{H}^{-1}}
{\ve{y}^T \ma{H}^{-1}\ve{s}} \:.
	\end{equation}

\subsubsection{SR1 update}
The symmetric-rank-1 update (SR1) in chosen in the form
	\begin{equation}\label{}
\delta \ma{H} = \ve{u}\ve{u}^T \;,
	\end{equation}
where the vector $\ve{u}$ is found from the condition~(\ref{o:eq:up}). The
resulting update is
	\begin{equation}\label{}
\ma{H} \rightarrow \ma{H}+
\frac{(\ve{y}-\ma{H}\ve{s})(\ve{y}-\ma{H}\ve{s})^T}
{(\ve{y}-\ma{H}\ve{s})^T\ve{s}} \:.
	\end{equation}
And for the inverse matrix
	\begin{equation}\label{}
\ma{H}^{-1} \rightarrow \ma{H}^{-1}+
\frac{(\ve{s}-\ma{H}^{-1}\ve{y})(\ve{s}-\ma{H}^{-1}\ve{y})^T}
{(\ve{s}-\ma{H}^{-1}\ve{y})^T\ve{y}} \:.
	\end{equation}

\subsubsection{Davidon–Fletcher–Powell (DFP) update}
	\begin{equation}
\delta\ma{H}^{-1}=\frac{\ve{y}\ve{y}^T}{\ve{y}^Ts}-
\frac{\ma{H}^{-1}\ve{s}\ve{s}^T\ma{H}^{-1}}{\ve{s}^T\ma{H}^{-1}\ve{s}}
	\end{equation}

%\begin{table}
%\renewcommand{\arraystretch}{2}
%\begin{tabular}{|c|c|c|}
%\hline
%Method & $\delta H$ & $\delta H^{-1}$ \\
%\hline
%DFP &
%$\frac{\ma{H}\ve{y}\ve{y}^T\ma{H}}{\ve{y}^T\ma{H}\ve{y}}
%+\frac{\ve{s}\ve{s}^T}{\ve{y}^T\ve{s}}$
%& \\
%\hline
%BFGS & & \\
%\hline
%Broyden &
%$
%\frac{\ve{u}-H\ve{s}}
%{\ve{s}^T\ve{s}}\ve{s}^T
%$
%&
%$\left.\frac{(\ve{s}-H^{-1}\ve{y})\ve{y}^TH^{-1}}
%{\ve{y}^T H^{-1}\ve{s}}\right.$ \\
%\hline
%SR1 &
%$\frac{(\ve{y}-H\ve{s})(\ve{y}-H\ve{s})^T}
%{(\ve{y}-H\ve{s})^T\ve{s}}$
%&
%$
%\frac{(\delta\ve{x}-H^{-1}\ve{d})(\delta\ve{x}-H^{-1}\ve{d})^T}
%{(\delta\ve{x}-H^{-1}\ve{d})^T\ve{d}}
%$\\
%\hline
%\end{tabular}
%\end{table}

\subsection{Downhill simplex method}
The {\em downhill simplex method} (also called Nelder-Mead method or
amoeba method) is a commonnly used minimization algorithm, where the
minimum of a function in an $n$-dimensional space is found by transforming
a simplex---a polytope with $n$+1 vertexes---according to the function
values at the vertexes, moving it downhill until it converges towards
the minimum.

The advantages of the downhill simplex method is its stability and the
lack of use of derivatives. However, the convergence is realtively slow
as compared to Newton's methods.

In order to introduce the algorithm we need the following definitions:
\begin{itemize}

\item Simplex: a figure (polytope) represented by $n$+1 points, called
vertexes, $\{\ve{p}_1,\dots,\ve{p}_{n+1}\}$ (where each point
$\ve{p}_{k}$ is an $n$-dimensional vector).

\item Highest point: the vertex, $\ve{p}_\mathrm{hi}$,
with the highest value of the function:
$f(\ve{p}_\mathrm{hi})=\max _k f(\ve{p}_{k})$.

\item Lowest point: the vertex, $\ve{p}_\mathrm{lo}$, with the
lowest value of the function:
$f(\ve{p}_\mathrm{lo}) = \min _k f(\ve{p}_{k})$.

\item Centroid: the center of gravity of all points, except for the
highest:
$\ve{p}_\mathrm{ce} = \frac{1}{n}\sum_{(k\ne\mathrm{hi})}\ve{p}_k$

\end{itemize}

The simplex is moved downhill by a combination of the following elementary
operations:
\begin{enumerate}

\item Reflection: the highest point is reflected against the centroid,
$\ve{p}_\mathrm{hi} \rightarrow  \ve{p}_\mathrm{re}
= \ve{p}_\mathrm{ce} + (\ve{p}_\mathrm{ce} -
\ve{p}_\mathrm{hi})$.

\item Expansion: the highest point reflects and then doubles its
distance from the centroid, $\ve{p}_\mathrm{hi} \rightarrow
\ve{p}_\mathrm{ex} = \ve{p}_\mathrm{ce} +
2(\ve{p}_\mathrm{ce} - \ve{p}_\mathrm{hi})$.

\item Contraction: the highest point halves its distance from the
centroid,
$\ve{p}_\mathrm{hi} \rightarrow
\ve{p}_\mathrm{co} = \ve{p}_\mathrm{ce} +
\frac{1}{2}(\ve{p}_\mathrm{hi} - \ve{p}_\mathrm{ce})$.

\item Reduction: all points, except for the lowest, move towards the
lowest points halving the distance.
$\ve{p}_{k\ne\mathrm{lo}} \rightarrow
\frac{1}{2}(\ve{p}_{k}+\ve{p}_\mathrm{lo})$.

\end{enumerate}
A C-implementation of the simplex operations is listed in
Table~\ref{o:t:simpo}.

Finally, Table~\ref{o:t:dsa} shows one possible algorithm for the downhill
simplex method.
\begin{table}
\caption{Downhill simplex (Nelder-Mead) algorithm}
\label{o:t:dsa}
\lstset{frame=single,basicstyle=\footnotesize}
\begin{lstlisting}[frame=single,mathescape=true,tabsize=2]
REPEAT :
	find highest, lowest, and centroid points of the simplex
	try reflection
	IF $f(\mathrm{reflected})< f(\mathrm{lowest})$ :
		try expansion
		IF $f(\mathrm{expanded})< f(\mathrm{reflected})$ :
			accept expansion
		ELSE :
			accept reflection
	ELSE :
		IF $f(\mathrm{reflected})< f(\mathrm{highest})$ :
			accept reflection
		ELSE :
			try contraction
			IF $f(\mathrm{contracted})< f(\mathrm{highest})$ :
				accept contraction
			ELSE :
				do reduction
UNTIL converged (e.g. size(simplex)<tolerance)
\end{lstlisting}
\end{table}

\section{Global optimization}

Global optimization is the problem of locating (a good approximation to)
the global minimum of a given objective function in a search space large
enough to prohibit exhaustive enumeration.

When only a small sub-space of the search space can be realistically
sampled the stochastic methods usually come to the fore.

A good local minimizer converges to the nearest local miminum relatively
fast, so one possible global minimizer can be constructed by simply
starting the local miminizer from different random starting points.

In the following several popular global minimization algorithms are
shortly described.

\subsection{Simulated annealing}

Simulated annealing is a stochastic metaheuristic algorithm for global
minimization. The name and inspiration come from annealing --- heating
up and cooling slowly --- in material science.  The slow cooling allows
a piece of material to reach a state with "lowest energy".

The objective function in the space of states is interpreted as some
sort of potential energy and the states --- the points in the search
space --- are interpeted as physical states of some physical system.
The system attempts to make transitions from its current state to some
randomly sampled nearest states with the goal to eventually reach the
state with minimal energy -- the global minimum.

The system is attached to a thermal resevoir with certain
temperature. Each time the energy of the system is measured the
reservoir supplies it with a random amount of thermal energy
sampled from the Boltzmann distribution,
	\begin{equation}
P(E) = Te^{-E/T} \:.
	\end{equation}

If the temperature equals zero the system can only make transitions
to the neighboring states with lower potential energy.  In this case
the algorithm turns merely into a local minimizer with random sampling.

If temperature is finite the system is able to climb up the ridges of
the potential energy --- about as high as the current temperature ---
and thus escape from local minima and hopefully eventually reach the
global minimum.

One typically starts the simulation with some finite temperature on the
order of the height of the typical hills on the potential energy surface,
letting the system to wander almost unhindered around the landscape with a
good chance to locale if not the best then at least a good enough minimum.
The temperature is then slowly reduced following some annealing schedule
which may be supplied by the user but must end with $T=0$ towards the
end of the alloted time budget.

Table~\ref{o:t:sa} lists one possible variant of the algorithm.

\begin{table}
\caption{Simulated annealing algorithm}
\label{o:t:sa}
\lstset{frame=single,basicstyle=\footnotesize}
\begin{lstlisting}[frame=single,mathescape=true,tabsize=2]
state $\leftarrow$ start_state
T $\leftarrow$ start_temperature
energy $\leftarrow$ E(state) + random_thermal_energy(T)
REPEAT :
	new_state $\leftarrow$ neighbour(state)
	new_energy $\leftarrow$ E(new_state) + random_thermal_energy(T)
	IF new_energy < energy :
		state $\leftarrow$ new_state
		energy $\leftarrow$ new_energy
	reduce_temperature_according_to_schedule(T)
UNTIL terminated
\end{lstlisting}
\end{table}

The function {\tt neigbour(state)} should return a randomly chosen neighbour
of the given state.

Downhill simplex method can incorporate simulated annealing by
adding the stochastic thermal energy to the values of the objective
function at the verices.

\subsection{Quantum annealing}

Quantum annealing is a general global minimization algorithm which ---
like simulated annealing --- also allows the search path to escape from
local minima.  However instead of the thermal jumps over the potential
barriers quantum annealing allows the system to tunnel through the
barriers.

In its simple incarnation the quantum annealing algorithm allows the
system to attempt transitions not only to the nearest states but also
to distant states within certain "tunneling distance" from the current
state.  The transition is accepted only if it reduces the potential
energy of the system.

At the beginning of the minimization prcedure the tunnelling distance
is large --- on the order of the size of the region where the global
minimum is suspected to be localed --- allowing the system to explore
the region.  The tunneling distance is then slowly reduced according
to a schedule such that by the end of the alloted time the tunnelling
distance reduces to zero at which point the system hopefully is in the
state with minimal energy.
\begin{table}
\caption[tbh]{Quantum annealing algorithm}
\label{o:t:qa}
\lstset{frame=single,basicstyle=\footnotesize}
\begin{lstlisting}[frame=single,mathescape=true,tabsize=2]
state $\leftarrow$ start_state
energy $\leftarrow$ E(state)
R $\leftarrow$ start_radius
REPEAT :
	new_state $\leftarrow$ neighbour_within_radius(state,R)
	new_energy $\leftarrow$ E(new_state)
	IF new_energy < energy :
		state $\leftarrow$ new_state
		energy $\leftarrow$ new_energy
	reduce_radius_according_to_schedule(R)
UNTIL terminated
\end{lstlisting}
\end{table}

\subsection{Evolutionary algorithms}

Unlike annealing algorithms, which follow the motion of only one point
in the search space, the evolutionary algorithms typically follow a set
of points called a population of individuals. A bit like the downhill
simplex method which follows the motion of a set of points -- the simplex.

The population evolves towards more fit individuals where fitness
is understood in the sence of minimizing the objective functions. The
parameters of the individuals (for example, the coordinates of the points
in multi-dimentional minimization of a continuous objective function)
are called genes.

The algorithm proceeds iteratively with the population in each iteration
called a generation.  In each generation the fitness of each individual
--- typically, the value of the objective function --- is evaluated and
the new generation is generated stochastically from the gene pool of
the current generation through crossovers and mutations such that the
genes of more fit individuals have a better chance of propagating into
the next generation.

Each new individual in the next generation is produced from a pair of
"parent" individuals of the current generation.  The use of two "parents"
is biologically inspired, in practice more than two "parents" can be
used as well.  The parents for a new individual are selected from the
individuals of the current generation through a fitness based stochastic
process where fitter individuals are more likely to be selected.

The "child" individual shares many characteristics of its "parents".
In the simplest case the "child" may get its genes by simply averaging
the genes of its parents.  Then a certain amount of mutations --- random
changes in the genes --- are added to the "child's" genes.

Generation of "children" continues until the population of the new
generation reaches the appropriate size after which the iteration
repeats itself.

The algorithm is terminated when the fitness level of the population is
deemed sufficient or when the allocated budget is exhausted.

\section{Implementation in C}

\begin{table}
\caption{C implementation of simplex operations}
\label{o:t:simpo}
\lstset{frame=single,basicstyle=\footnotesize,tabsize=3}
\lstinputlisting[language=c]
{c/downhill.simplex/simplex_operations.c} 
\lstinputlisting[language=c]
{c/downhill.simplex/simplex_update.c} 
\end{table}

\begin{table}
\label{o:t:a}
\caption{C implementation of downhill simplex algorithm}
\lstset{frame=single,basicstyle=\footnotesize,tabsize=2}
\lstinputlisting[language=c]
{c/downhill.simplex/downhill_simplex.c} 
\end{table}
