\svnidlong
{$HeadURL: $}
{$LastChangedDate: $}
{$LastChangedRevision: $}
{$LastChangedBy: $}

\odachapter{Calibration methods available in \oda}
\begin{tabular}{p{4cm}l}
\textbf{Origin:} & CTA memo200802\\
\textbf{Last update:}    & \svnfilemonth-\svnfileyear\\
\end{tabular}

In the calibration or parameter estimation algorithms, the basic idea is to
find the set of model parameters which minimizes the cost function measuring
the distance between the observation and the model prediction. Two different
cost functions are implemented. The first one is similar to equation
(\ref{eq.costfunction1}), while in the second one we add the background
component as follows
\begin{equation}
\label{eq.costfunction2}
  J(x_o) = (x_o-x^b)'(P^b)^{-1}(x_o-x^b)+\sum_{k=1}^N (y^o(k)-Hx(k))R^{-1}(y^o(k)-Hx(k))'
\end{equation}
Where $x^b$ and $P^b$ are the background or initial estimate of $x_o$ and its covariance respectively. This additional component ensures that the solution will not be too far from the initial guess.

The following calibration methods are implemented in \oda :
\begin{itemize}
\item Dud
\item Sparse Dud
\item Simplex
\item Powell
\item Gridded full search
\item Shuffled Complex Evolution (SCE)
\item Generalized Likelihood Uncertainty Estimation (GLUE)
\item (L)BFGS
\item Conjugate Gradient: Fletcher-Reeves, Polak-Ribiere, Steepest Descent
\end{itemize}

\section{Parameter estimation with the Simplex method}
\label{sssec.simplex}
The simplex method that is implemented in COSTA is the one due to Nelder and
Mead (\cite{NelderMead1965}). It is a systematic procedure for generating and
testing candidate vertex solutions to a minimization problem. It begins at an
arbitrary corner of the solution set. At each iteration, the simplex method
selects the variable that will produce the largest change towards the minimum
solution. That variable replaces one of its compatriots that is most severely
restricting it, thus moving the simplex to a different corner of the solution
set and closer to the final solution.

A simplex is the geometrical figure consisting, in \emph{N} dimensions, of
$N+1$ points and all their interconnecting line segments, polygonal faces,
etc. In two dimensions, a simplex is a triangle. In three dimensions it is
a tetrahedron. The simplex method must be started with $N+1$ points of
initial guess, defining the initial simplex. The simplex method now takes a
series of steps. The first step is to move the vertex where the cost is
largest through the opposite face of simplex to a lower point. This step is
called the \emph{reflect} step. When the cost of the new vertex is even
smaller than all the remaining, the method expands the simplex even
further, called the \emph{expand} step. If none of these steps produce a 
better vertex, the method will contract the simplex in the same direction
of the previous step and take a new point. This step is called the
\emph{contract} step. When the cost of the new point is still not better
than the previous one, the method will take the last step called
\emph{shrink}. In this step all of the simplex points, except the one with
the lowest cost, are 'shrinked' toward the best vertex.

As it basically only tries and compares the solution of several different
sets of parameters, the method requires only function evaluations and not
derivatives.

\subsection{Using the Simplex method}
For using the implemented simplex method the user needs to specify in the
input file the \emph{initial-guess} of the set of parameters to calibrate
as well as the \emph{initial-step} for creating the initial simplex. The
initial-step consists of $N$ entries, where $N$ is the number of parameters
to calibrate. An initial vertex is obtained by adding an element of the
initial-step to the corresponding parameter in the initial-guess.
Performing this one by one to all the parameters in the initial-guess,
there are $N+1$ vertices forming the initial simplex. Although it can be
easily extended, at the moment the program only supports the model with
a maximum of 10 parameters.

The stopping criteria for the iteration are the maximum number of iteration
and the maximum cost difference between the worst and the best vertices,
where worst vertex refers to the one with the biggest cost value and best
vertex is the one with the lowest. When the maximum cost difference is very
small we may expect that the (local) minimum of the cost function has been
reached. The optimum solution is the vertex with the lowest cost. However,
as output we display all the final vertices with their respective costs.

Since the model parameters are stored as COSTA state vector, most variable
operations are performed in term of COSTA state vector. The operations are
usually to compute new vertex. This can easily be carried out by using the
combination of functions like cta\_state\_axpy and cta\_state\_scal, with
the aid of functions like cta\_state\_duplicate and cta\_state\_copy for
assigning values to working variables. These COSTA functions reduce
significantly the lines of codes required to perform the operation.

An important function required to implement the simplex method is the sorting
function. This function is used to sort the vertices with respect to their
cost values in descending order. While the vertices are stored as COSTA
state vectors, their cost values are stored as a Fortran array. The sorting
is done by using an external routine, which works only with Fortran
variables. The output of this routine are the sorted array of cost values
and an integer array containing the index of the sorted array.

In the section above  we learnt that the simplex method consists
of several steps which are taken to find a new vertex with lower cost. When
all the steps do not produce a better vertex the program will give a
warning. For some cases this may indicate that the solution does not exist.
This occurs for example with the Lorenz model. Since it is a chaotic model,
small difference in the parameters yields very different cost values. This
makes the Lorenz model not a very good model for calibration tests in the
present setup. Perhaps shortening the interval over which the optimization
is carried out can make the cost function better behaved.
\subsection{The configuration of the Simplex method}
\begin{itemize}
\item {\tt simulation\_span}: the overall timespan to run the model
\item {\tt output}: output specification
\item {\tt iteration}: Number of iterations and convergence tolerance
\item {\tt model<n>}: Model configuration for the initial model state at every
starting vertex
\end{itemize}

\subsection{XML-example}

\begin{verbatim}
  <parameter_calibration>
    <CTA_TIME id="simulation_span" start="0" stop="50.0" step=".1"/>
    <output>
        <filename>results_oscill.m</filename> 
       <CTA_TIME id="times_mean" start="0" stop="50.0" step="1"/>
    </output>
    <iteration maxit="60" tol="0.001">    </iteration>
    <!-- 1st VERTEX:-->
    <model1>
       <modelbuild_sp>
         <xi:include href="models/functions_oscill.xml"/>
          <model>  
             <parameters avg_t_damp="8.95" avg_omega="13.5"> </parameters> 
          </model> -->  
       </modelbuild_sp> 
    </model1>
    <!-- 2nd VERTEX:-->
    <model2>
       <modelbuild_sp>
         <xi:include href="models/functions_oscill.xml"/>
          <model>  
             <parameters avg_t_damp="6.0" avg_omega="15.9"> </parameters> 
          </model> -->  
       </modelbuild_sp> 
    </model2>
    <!-- 3rd VERTEX:-->
    <model3>
       <modelbuild_sp>
         <xi:include href="models/functions_oscill.xml"/>
          <model>  
             <parameters avg_t_damp="9.0" avg_omega="12.5"> </parameters> 
          </model> -->  
       </modelbuild_sp> 
    </model3>
  </parameter_calibration>
\end{verbatim}
 %\caption{Example of input file for simplex method\label{Fig:Simplex_input}}

\section{Parameter estimation with the Conjugate Gradients method}
\label{sssec.conjugrad}
The problem of minimization of multivariable function is usually solved by
determining a \emph{search direction} vector and solve it as a line
minimization problem. If $x$ is a vector containing the variables to be
determined and $h$ is the vector of search direction, at each iteration
step the minimization problem of a function $f$ is formulated as to find
the step size $\lambda$ that minimizes $f(x+\lambda h)$. At the next
iteration, $x$ is replaced by $x+\lambda h$ and a new search direction is
determined. Different methods basically propose different ways of finding
the search direction.

The conjugate gradient method is an algorithm for finding the nearest local
minimum of a function which uses \emph{conjugate directions} for going
downhill. Two vectors $u$ and $v$ are said to be conjugate (with respect to
a matrix $A$) if
\begin{equation}
   u'Av=0
\end{equation}
where in the minimization problem, $A$ is typically the \emph{Hessian}
matrix of the cost function. In the conjugate gradient methods, the search
direction is somehow constructed to be conjugate to the old gradient.

The two most important conjugate gradient methods are the
\emph{Fletcher-Reeves} and the \emph{Polak-Ribierre} methods (TODO: Press
et.al., 1989). These algorithms calculate the mutually conjugate directions
of search with respect to the Hessian matrix of the cost function directly
from the function and the gradient evaluations, but without the direct
evaluation of the Hessian matrix. The new search direction $h_{i+1}$ is
determined by using
\begin{equation}
   h_{i+1} = g_{i+1}+\gamma_i h_i
\end{equation}
where $h_i$ is the previous search direction, $g_{i+1}$ is the negative of
local gradient at iteration step $i+1$, while $\gamma_i$ is determined by
using the following equations for \emph{Fletcher-Reeves} and the
\emph{Polak-Ribierre} methods respectively:

\begin{eqnarray}
   \gamma_i=\frac{g_{i+1} \cdot g_{i+1}}{g_{i} \cdot g_{i}} \\
   \gamma_i=\frac{(g_{i+1}-g_{i}) \cdot g_{i+1}}{g_{i} \cdot g_{i}}
\end{eqnarray}

If the vicinity of the minimum has the shape of a long, narrow valley, the
minimum is reached in far fewer steps than would be the case using the
\emph{steepest descent} method, which makes use of minus of the local
gradient as the search direction.

In this study, the line minimization to find the step size $\lambda$ that
minimizes $f(x+\lambda h)$ at every iteration step is done by using the
\emph{golden section search} algorithm. This is an elegant and robust
method of locating a minimum of a line function by bracketing it with three
points: if we can find three points $a,b,$ and $c$ where $f(a)>f(b)<f(c)$
then there must exist at least one minimum point in the interval $(a,c)$.
The points $a,b,$ and $c$ are said to \emph{bracket} the minimum. This
algorithm involves evaluating the function at some point $x$ in the larger
of the two intervals $(a,b)$ or $(b,c)$. If $f(x)<f(b)$ then $x$ replaces
the midpoint $b$ and $b$ becomes an end point. If $f(x)>f(b)$ then $b$
remains the midpoint with $x$ replacing one of the end points. Either way
the width of the bracketing interval will reduce and the position of the
minimum will be better defined. The procedure is then repeated until the
width achieves a desired tolerance. It can be shown that if the new test
point, $x$, is chosen to be a proportion $(3-\sqrt{5})/2$ (hence Golden
Section) along the larger sub-interval, measured from the mid-point $b$,
then the width of the full interval $(a,c)$ will reduce at an optimal rate.

In the absence of an adjoint formalism in COSTA a numerical approximation
was used to compute the gradient of the cost function.


\subsection{Using the Conjugate gradient method}
In this COSTA-method both the Fletcher-Reeves and Polak-Ribiere methods are
implemented. For comparison purpose, we also implement a steepest-descent
method. The user can choose the method to use by specifying the field
method \emph{id} in the input file, where 1 refers to Fletcher-Reeves, 2 to
Polak-Ribiere, and 3 to steepest-descent.

The stopping criteria are the maximum number of iterations, the tolerance
for the step size $\lambda$ as described in subsection
\ref{sssec.conjugrad} and the tolerance for the local gradient. When the
line-minimization does not move the parameters significantly to the new
ones, the cost function may have reached its (local) minimum. On the other
hand, the gradient around the (local) minimum is also close to zero. Hence
it is possible to use the norm of the gradient vector as one of the
stopping criteria. The same tolerance is also used within the golden
section search routine.

Besides specifying the method to use and the stopping criteria parameters,
the user also needs to specify the initial guess of the parameters to
calibrate as well as the initial two points $AX$ and $BX$ for the golden
section search minimization.

Since the method requires the computation of the gradient, while there are no
adjoint models available, the gradient of the cost function is computed
with finite difference. This is fine for small number of parameters like
the ones used in this study. However it will become computationally
expensive for many parameters. The size of the perturbation, $delta$, must
also be specified in the input file.

In our implementation, the gradient is also stored as COSTA state to make
it easier to assign values between the variables representing model
parameters to calibrate. The norm of the gradient is computed using the
function cta\_state\_nrm2.

Besides the state vector operations, like scaling and adding, in this
method we also need to compute the dot product between two state vectors.
This can easily be done by using the function cta\_state\_dot.


\subsection{The configuration of the CG method}
\begin{itemize}
\item {\tt simulation\_span}: the overall timespan to run the model
\item {\tt maxit}: the maximum number of iterations
\item {\tt tol\_step, tol\_grad}: Stopping criterium  according to the step
  size or the gradient norm
%\item {\tt output}: output specification
\item {\tt AX, BX}: Starting point for golden section search
\item {\tt delta}: Perturbation size 
\item {\tt iteration}: Number of iterations and convergence tolerance
\item {\tt model}: Model configuration for the initial model
\end{itemize}

\subsection {XML-example}
   \begin{verbatim}
  <parameter_calibration>
    <CTA_TIME id="simulation_span" start="0" stop="50.0" step=".1"/>
    <iteration maxit="160" tol_step="0.0001" tol_grad="0.0001">
    </iteration>
    <method id="1" AX="0.0" BX="0.02" delta="1E-8">
    <!-- id=1: Fletcher-Reeves, id=2: Polak-Ribiere, id=3: Steepest-descent-->
    </method>
    <model>
       <modelbuild_sp>
       <functions>
          <create>oscillparam_create</create>
          <covariance>oscillparam_covar</covariance>
          <getobsvals>oscillparam_obs</getobsvals>
          <compute>oscillparam_compute</compute>
       </functions>
       <model>
          <parameters avg_t_damp="8.5" avg_omega="12.9">
          </parameters>
       </model>
       </modelbuild_sp>
    </model>
    <paramstd param1="3" param2="3">
    </paramstd>
  </parameter_calibration>
   \end{verbatim}

\section{Parameter estimation with the LBFGS method}
For the problem of minimizing a multivariable function quasi-Newton methods
are widely employed. These methods involve the approximation of the Hessian
(or its inverse) matrix of the function. The LBFGS (Limited
memory-Broyden-Fletcher-Goldfarb-Shanno) method is basically a method to
approximate the Hessian matrix in the quasi-Newton method of optimization.
It is a variation of the standard BFGS method, which is given by (\cite{Nocedal1980}; \cite{Byrdetall1994})
\begin{equation}
  \label{eq.LBFGS_linemin}
    x_{i+1}=x_i-\lambda_iH_ig_i, \space \space i=0,1,2,\cdots
\end{equation}
where $\lambda_i$ is a steplength, $g_i$ is the local gradient of the cost
function, and $H_i$ is the approximate inverse Hessian matrix which is
updated at every iteration by means of the formula
\begin{equation}
    H_{i+1}=V'_i H_i V_i + \rho_i s_i s'_i
\end{equation}
where
\begin{eqnarray}
    \rho_i=\frac{1}{y'_i s_i} \\
    V_i=I-\rho_i y_i s'_i
\end{eqnarray}
and
\begin{eqnarray}
    s_i=x_{i+1}-x_i \label{eq.LBFGS_s} \\
    y_i=g_{i+1}-g_i \label{eq.LBFGS_y}
\end{eqnarray}

Using this method, instead of storing the matrices $H_i$, one stores a
certain number of pairs, say $m$, of pairs $\{s_i,y_i\}$ that define them
implicitly. The product of $H_i g_i$ is obtained by performing a sequence
of inner products involving $g_i$ and the $m$ most recent vector pairs
$\{s_i,y_i\}$ to define the iteration matrix.

Like in conjugate gradient methods, the line minimization for determining
$\lambda$ in equation (\ref{eq.LBFGS_linemin}) is implemented by using the
Golden Section Search algorithm.

\subsection{Using the LBFGS method}
In the implementation of LBFGS method we use COSTA matrices to store the
$s_i$ and $y_i$ vectors in equation (\ref{eq.LBFGS_s}) and
(\ref{eq.LBFGS_y}). These vectors are stored as the column of the COSTA
matrices. The elimination of the oldest vectors, however, requires the
combination of a function for getting a column of a COSTA matrix and
another function for setting the column. The latter function is already
available in the present COSTA, i.e. cta\_matrix\_setcol. However, the
former one is not yet available. We resorted to worked around this by
creating a loop consisting a combination of cta\_matrix\_getval and
cta\_vector\_setvals for performing this operation. 

Vectors $s_i$ and $y_i$ are stored as a COSTA state for the same reason as in
the implementation of the conjugate gradient and simplex methods. The
operations with a COSTA matrix, however, require the variable to be stored as a 
COSTA vector. Therefore at some points in the program we were required to
use cta\_state\_getvec and cta\_state\_setvec to get the values of a COSTA
state and assign it to a COSTA vector and vice versa.

\subsection{The configuration of the LBFGS method}

Like in the conjugate gradient method, the implemented LBFGS method
requires maximum number of iterations, step-size tolerance, tolerance for
gradient, and the perturbation size, $delta$, for computing the gradient.
However for using this method the user also needs to specify $nstore$,
which is the number of vectors $s_i$ and $y_i$ to store.

\begin{itemize}
\item {\tt simulation\_span}: the overall timespan to run the model
\item {\tt output}: output specification
\item {\tt iteration}: Number of iterations and convergence tolerance, {\tt
  maxln} is the maximum number of iteration steps in the line search part.
\item {\tt model}: Model configuration for the initial model
\item {\tt method}: Method configuration.
    {\tt nstore} specifies the storage size, i.e. max number of vector s and y to store.
    {\tt c1} and {\tt c2} are the parameters used in the 1st and 2nd Wolfe
    conditions (these are used in the line search part). 
    {\tt delta} specifies the perturbation size in computing the gradient.
\end{itemize}


\subsection{XML-example}
\begin{verbatim}
  <parameter_calibration>
    <!-- Simulatie timespan en stapgrootte via CTA_Time -->
    <CTA_TIME id="simulation_span" start="0" stop="50.0" step=".1"/>
    <output>
        <filename>results_oscill.m</filename> 
       <CTA_TIME id="times_mean" start="0" stop="50.0" step="1"/>
    </output>
    <iteration maxit="10" maxln="20" tol_step="0.0001" tol_grad="0.0002">
    </iteration>
    <method nstore="3" c1="1E-4" c2="0.5" delta="1E-8">

    </method>
       <model>
         <xi:include href="models/lorenz_sp.xml"/>
       </model>
  </parameter_calibration>
\end{verbatim}

\section{Parameter estimation with Dud}

Dud (Doesn't use derivative) is one of the optimization algorithms, which do
not use any derivative of the function being evaluated. It can be seen as a
Gauss-Newton method, in the sense that it transforms the nonlinear least square
problem into the well-known linear square problem. The difference is that
instead of approximating the nonlinear function by its tangent function, the
Dud uses an affine function for the linearization. For $N$ calibration
parameters, Dud requires $(N+1)$ set of parameters estimates. The affine
function for the linearization is formed through all these $(N+1)$ guesses.
Note that the affine function gives exact value at each of the $(N+1)$ points.
The resulting least square problem is then solved along the affine function to
get a new estimate, whose cost is smaller than those of all other previous
estimates. If it does not produce a better estimate, the Dud will perform
different steps, like searching in opposite direction and/or decreasing
searching-step, until a better estimate is found. Afterwards, the estimate with
the largest cost is replaced with the new one and the procedure is repeated for
the new set of $(N+1)$ estimates. The procedure is stopped when one of the
stopping criteria is fulfilled.

\subsection{Dud Algorithm}

Suppose we have a numerical model with a vector of $p$ uncertain parameters,
say $\mathbf{x} \in \mathbb{R}^p$ . Let $\mathbf{y} \in \mathbb{R}^n$ be a set
of n data points and let $f_i : \mathbb{R}^p \rightarrow \mathbb{R}$ be the
model prediction corresponding to the $i^{th}$ data point. The goal is to find
a vector of parameters $\mathbf{x}^{opt} \in \mathbb{R}^p$ that minimizes the
following costfunction:
\begin{equation}\label{eq:dudcost}
(\mathbf{x}) = [\mathbf{y} - \mathbf{f(x)}]^T [\mathbf{y} - \mathbf{f(x)}].
\end{equation}
The idea of the Dud algorithm is similar to that of the Gauss-Newton algorithm.
The key difference is that the Dud algorithm uses an approximation of gradients
instead of the exact gradients. In short the Dud algorithm works as follows.

The memory of the algorithm always contains $p + 1$ estimations of the optimal
parameter set and their corresponding model predictions. First the model
functions $f_i(\mathbf{x})$ are linearized based on interpolation between the
model predictions that are stored in memory. Then the linearized version of
equation~\ref{eq:dudcost} is solved, yielding a new estimation of the optimal
parameter set. Assuming the new estimation is better than all previous ones,
it replaces the worst estimation in the memory of the algorithm. This process
is repeated until the best estimation is ”sufficiently” close to the optimal
parameter set.

Let us now formulate the algorithm in more detail. Let $\mathbf{x}^1 , \ldots
\mathbf{x}^{p+1}$ and
$\mathbf{f}(\mathbf{x}^1), \ldots, \mathbf{f}(\mathbf{x}^{p+1})$ be the
parameter estimations and corresponding model predictions that are stored in
memory. The first step is to fill this memory with $p+1$ initial guesses. These
guesses are sorted according to the value of their respective costfunctions
$S(\mathbf{x}^i)$, meaning $\mathbf{x}^1$ is the parameter set with the highest
costfunction and $\mathbf{x}^{p+1}$ the one with the lowest.

Then the model functions $f_i(\mathbf{x})$ are linearized by interpolation
between the model predictions that are stored in memory. Let us denote the
linearized model functions by $l_i(\mathbf{x})$. Then
\begin{equation}
l_i (\Delta \mathbf{x}) = f_i(\mathbf{x}^{p+1}) + M \Delta \mathbf{x}
\end{equation}
with $\Delta \mathbf{x} = \mathbf{x} - \mathbf{x}^{p+1}$ and $M$ the
approximation of the Jacobian matrix by interpolation. Note that $M = F P^{−1}$
where $F$ is the matrix with $j^{th}$ column equal to $\mathbf{f}(\mathbf{x}^j)
- \mathbf{f}(\mathbf{x}^{p+1})$ and $P$ the matrix with $j^{th}$ column equal
to $\mathbf{x}^i - \mathbf{x}^{p+1}$.

Then $\mathbf{f}(\mathbf{x}) = \mathbf{l}(\mathbf{x})$ is substituted into cost
equation~\ref{eq:dudcost}. Let us call the linearized costfunction
$Q(\mathbf{x})$. By solving $Q'(\Delta \mathbf{x}) = 0$ we find that for the
optimal value of $\Delta \mathbf{x}$ holds
\begin{equation}
M^T M \Delta \mathbf{x} = M^T \left[\mathbf{y} - \mathbf{f}(\mathbf{x}^{p+1})\right]
\end{equation}
yielding a new parameter estimation $\mathbf{x}^\ast = \mathbf{x}^{p+1}
+ \Delta \mathbf{x}$. If $\mathbf{x}^\ast$ has a lower costfunction than
$\mathbf{x}^{p+1}$, the worst estimate $\mathbf{x}^1$ is tossed out of the
memory, making place for the new estimation. The elements in the memory are
again sorted according to their costfunctions, so $\mathbf{x}^{p+1}
= \mathbf{x}^\ast$ , $\mathbf{x}^p = \mathbf{x}^{p+1}$ etc.

It may however happen that the new estimation $\mathbf{x}^\ast$ is not better
than one or more of the previous estimations. In this case a line search is
done. A better estimation is searched in the direction from $\mathbf{x}^{p+1}$
to $\mathbf{x}^\ast$ , so on the line
\begin{equation}
\mathbf{r}(\epsilon) = \mathbf{x}^{p+1} + \epsilon (\mathbf{x}^\ast - \mathbf{x}^{p+1} )
\end{equation}
The step size $\epsilon \in \mathbb{R}$ is iteratively being reduced until a
step size $\epsilon^\ast$ is found for which $\mathbf{r}(\epsilon^\ast)$ has a
lower costfunction than $\mathbf{x}^{p+1}$. Note that $\epsilon^\ast$ may very
well be negative. The new estimation $\mathbf{r}(\epsilon^\ast)$ is then stored
into memory as before.



\section{Dud with constraints}

The parameters $x_i, i \in {1, \ldots, p}$ in the calibration process sometimes
represent phyisical parameters that should for instance be positive or, based
on measurements, need to lie in a specific interval. Due to the linearization
process in the Dud algorithm it may very well happen that at some point the Dud
algorithm proposes a parameter $x_j$ which does not belong to an expected
interval $[a_j , b_j]$. This could blow up the model computation, but it may
also be the case that the model automatically corrects such an unexpected value
to a value within the expected interval. In the latter case, the Dud algorithm
may be ruined. To prevent such unpredicted behaviour, the possibilty to add a
set of constraints to the minimization problem in the Dud algorithm was
created.

The method of quadratic programming was used for this. First note that
minimizing $Q(\Delta \mathbf{x})$ is equivalent to minimizing
\begin{eqnarray}
q(\Delta \mathbf{x}) &=& \Delta \mathbf{x}^T M^TM \Delta \mathbf{x} − 2\Delta \mathbf{x}^T M^T \left[\mathbf{y} - \mathbf{f}(\mathbf{x}^{p+1})\right]\nonumber\\
 &=& \Delta \mathbf{x}^T G\Delta \mathbf{x} - 2\Delta \mathbf{x}^T \mathbf{c}
\end{eqnarray}
where $G = M^T M$ and $\mathbf{c} = M^T \left[\mathbf{y}
- \mathbf{f}(\mathbf{x}^{p+1})\right]$.

To solve this minimization problem with inequality constraints, we need to know
how to solve it with equality constraints first. Suppose for now that we have a
set of $m$ equality constraints such that $D\mathbf{x} = d$, with
$D \in \mathbb{R}^{m \times p}$ and d$ \in \mathbb{R}^m$ . We want to solve the following minimization problem
\begin{eqnarray}
\min_{x} & & q(\mathbf{x})\nonumber \\
\textrm{subject to} & & D\mathbf{x} = d
\end{eqnarray}
A vector $\lambda^\ast \in \mathbb{R}^m$ of Lagrange multipliers exists, such
that for the solution $\mathbf{x}^\ast$ holds
\begin{equation}\label{eq:lagrangianmultipliers}
\left[ \begin{array}{cc} G & -D^T \\ D & 0 \end{array} \right]
\left[ \begin{array}{c} \mathbf{x}^\ast \\ \lambda^\ast \end{array}\right]
= \left[ \begin{array}{c} -\mathbf{c} \\ \mathbf{d} \end{array}\right]
\end{equation}
This system can be rewritten in a form that is useful for computations by
expressing $\mathbf{x}^\ast$ as $\mathbf{x}^\ast = \mathbf{x} + \mathbf{s}$,
where $\mathbf{x}$ is some estimate of the solution and $\mathbf{s}$ is the
desired step. Rearranging equation~\ref{eq:lagrangianmultipliers} yields
\begin{equation}
\left[ \begin{array}{cc} G & D^T \\ D & 0 \end{array} \right]
\left[ \begin{array}{c} - \mathbf{s} \\ \lambda^\ast \end{array}\right]
= \left[ \begin{array}{c} \mathbf{g} \\ 0 \end{array}\right]
\end{equation}
where $\mathbf{g} = \mathbf{c} + G\mathbf{x}$ and $\mathbf{s} =
\mathbf{x}^\ast - \mathbf{x}$.

Now let us look at how to handle inequality constraints. We want to solve the
following problem
\begin{eqnarray}\label{eq:constrainedproblem}
\min_{\Delta x} & & q(\Delta \mathbf{x})\nonumber \\
\textrm{subject to} & & A \Delta \mathbf{x} \geq \Delta \mathbf{b}
\end{eqnarray}

where $\Delta \mathbf{b} = \left[ \begin{array}{c}
\mathbf{a} - \mathbf{x}^{p+1} \\ \mathbf{x}^{p+1} - \mathbf{b}
\end{array}\right]$
and $A = \left[ \begin{array}{c} I \\ -I \end{array}\right]$ with
$I \in \mathbb{R}^{p \times p}$ the identity matrix.

We find the solution $\Delta \mathbf{x}\ast$ iteratively with initial guess
$\Delta \mathbf{x} = 0$.

{\bf First Iteration}\\
First we ignore all constraints given by equation~\ref{eq:constrainedproblem},
yielding an optimal step $\mathbf{s}$ which satisfies $−G\mathbf{s}
= \mathbf{g}^0$ where $\mathbf{g}^0 = \mathbf{c} + G \Delta \mathbf{x}^0$.

If $A(\mathbf{s} + \Delta \mathbf{x}^0) \geq \Delta \mathbf{b}$, the
constraints are satisfied and we are done with $\Delta \mathbf{x}^\ast
= \mathbf{s} + \Delta \mathbf{x}^0$. If not, we find the constraint $i$ for
which $A_i(\mathbf{s} + \Delta \mathbf{x}^0) - \Delta b_i$ is smallest. Here
$A_i$ denotes the $i_{th}$ row of the matrix $A$. By choosing an $\alpha^0 \in
[0,1]$ such that $A_i (\alpha^0 \mathbf{s} + \Delta \mathbf{x}^0 ) = \Delta
b_i$ and setting $\Delta \mathbf{x}^1 = \alpha^0 \mathbf{s}
+ \Delta \mathbf{x}^0$ we ensure that $\Delta \mathbf{x}^1$ satisfies all
constraints. The $i^{th}$ constraint is then added to the so-called {\it
working set} $\mathbf{W}^1$. The working set $\mathbf{W}^k$ contains all
constraints which are active (i.e.\ act as equality constraints) at the
$k^{th}$ iteration.

{\bf $k^{th}$ Iteration}\\
We solve the following problem
\begin{eqnarray}
\min_{\mathbf{s}} && \mathbf{s}^TG\mathbf{s} + \mathbf{s}^T \mathbf{g}^k \nonumber \\
\textrm{subject to} && A_i \mathbf{s}= 0, i \in \mathbf{W}_k
\end{eqnarray}

If $\mathbf{s} = 0$, we check if all Lagrange multipliers $\lambda_i$ are
positive. If so, we are done with $\Delta \mathbf{x}^\ast
= \Delta \mathbf{x}^{k-1}$. If not, we remove the constraint corresponding to
the smallest Lagrange multiplier from working set $\mathbf{W}_{k-1}$, yielding
a new working set $W_k$. Setting $\Delta \mathbf{x}^k
= \Delta \mathbf{x}^{k-1}$ we start a new iteration.

If $\mathbf{s} \neq 0$, we check wether $\mathbf{s} + \Delta \mathbf{x}^{k-1}$
satisfies the inequality constraints that are not in working set
$\mathbf{W}^{k-1}$. If so, we set $\Delta \mathbf{x}^k = \mathbf{s}
+ \Delta \mathbf{x}^{k-1}$ and start the next iteration. If not we find the
constraint $j \notin \mathbf{W}^{k-1}$ for which $A_j (\mathbf{s}
+ \Delta \mathbf{x}^{k-1}) - \Delta b_j$ is smallest and choose $\alpha^k \in
[0,1]$ such that $A_j(\alpha^k \mathbf{s} + \Delta \mathbf{x}^{k-1}) = \Delta
b_j$. We add this constraint to the working set $\mathbf{W}_{k-1}$, yielding an
updated working set $\mathbf{W}_k$ . We set $\Delta \mathbf{x}^k
= \alpha^k \mathbf{s} + \Delta \mathbf{x}^{k-1}$ and start the next iteration.

