\documentclass[12pt]{article}
\usepackage{graphicx}
\usepackage{url}
\usepackage{subfigure}

\newcommand{\ds}{\displaystyle}
\newcommand{\bt}{\begin{tabbing}}
\newcommand{\et}{\end{tabbing}}
\renewcommand{\baselinestretch}{2}


\begin{document}
\title{BOXCQP: An algorithm for convex quadratic programming subject to bound constraints}
\author{C. Voglis and I.E. Lagaris \\
Dept. of Computer Science, University of Ioannina }

\maketitle

\begin{abstract}
A quadratic programming problem with positive definite Hessian and
bound constraints is solved, using a Lagrange multiplier approach.
The proposed method falls in the category of active set
techniques. The algorithm, at each iteration, modifies both the
minimization parameters in the primal space and the Lagrange
multipliers in the dual space. Comparative results of numerical
experiments are reported.
\end{abstract}

\noindent {\bf PACS:} 02.60.-x; 02.60.Pn; 07.05.Kf; 07.05.-t

\vskip 0.3cm \noindent {\bf Keywords: } Convex quadratic
programming, optimization, active set, Lagrange multipliers.


\section{Introduction}

The problem of minimizing a convex quadratic function  subject to
bound constraints appears quite frequently in applications. For
instance, many problems in computational physics and engineering,
are reduced to quadratic programming problems. Portfolio
management can also be formulated as quadratic programming
problem~\cite{PORT}. In the field of Artificial Intelligence, and
especially in Support Vector Machines (SVM) an efficient quadratic
solver is crucial for the training process~\cite{SVM1, SVM2}. Also
methods for calculating the radiation intensity in oncology
treatment are formulated as quadratic optimization problems \cite{DOSE}.
Finally, many non linear optimization techniques are based on
solving quadratic model subproblems~\cite{ACT1,ACT2,TRUSTQP}.\\
%
The Quadratic Programming problem with simple bounds is
stated as:
%
\begin{equation}
\label{QP} \min_x \  \frac{1}{2}x^TBx +x^Td,
\hbox{ subject to: } a_i \le x_i \le b_i, \forall i \in I=\{1,2,\cdots,N\}
\end{equation}
where $\ds x, d \in R^N$ and $\ds B $ a symmetric, positive
definite $\ds N \times N $ matrix.

For the problem in Eq.~(\ref{QP}) two major strategies exist in the
literature, both of which require feasible steps to be taken. The
first one is the Active Set strategy~\cite{ACT1,ACT2}, which generates
iterates on a face of the feasible box until either a minimizer
of the objective function is found or a point on the boundary
of that face is reached.
%Since the function values are strictly decreasing
%and the iterates are feasible, finite convergence can be proved.
The basic disadvantage of this approach, especially in
the large-scale case, is that constraints are added or removed one
at a time, thus requiring a number of iterations proportional to the
problem size. To overcome this, gradient projection
methods \cite{GP,GP2} were proposed. In that framework the active
set algorithm is allowed to add or remove many constraints per
iteration.

The second strategy consists in treating the inequality
constraints using interior point algorithms. In brief, an interior
point algorithm consists of a series of parametrized barrier
functions which are minimized using Newton's method. The major
computational cost is due to the solution of a linear system,
which provides a feasible search direction.

Our proposed approach to solve the problem of Eq.~(\ref{QP}) is an
active set algorithm, which does not guarantee strictly descent
iterations. In this paper we investigate a series of convex
quadratic test problems. We recognize that bound constraints are a
very special case of linear inequalities, which may in general
have the form $Ax \geq b$, $A$ being an $m x n$ matrix and $b$ is
a vector $\in R^m$. Our investigation is also motivated by the
fact that in the convex case, every problem subject to inequality
constraints can be transformed to a bound constrained one, using
duality, i.e. the problem:
\begin{eqnarray}
\label{PRIMAL}
\min_{x\in R^n} \  \frac{1}{2}x^TBx + x^Td  \\
\textrm{subject to: } Ax \geq b \nonumber
\end{eqnarray}

is equivalent to the dual:

\begin{eqnarray}
\label{DUAL}
\max_{y\in R^m} \  -\frac{1}{2}y^T\tilde{B}y + y^T\tilde{d} \\
\textrm{subject to: } y \geq 0 \nonumber
\end{eqnarray}
where $\tilde{B} = AB^{-1}A^T$  a positive definite matrix and
$\tilde{d} = AB^{-1}d+b$. The dual problem in Eq.~(\ref{DUAL}) is
also a quadratic problem subject to bounds. Let $y^*$ be the
solution of the dual problem. We can then obtain the solution to
the initial problem of Eq.~(\ref{PRIMAL}) as:
\begin{equation}
x^*  = B^{-1}(A^Ty^* - d)
\end{equation}

The paper is organized as follows. The proposed algorithm is
described in detail in Section \ref{SEC:SOLQUAD}. In Section
\ref{SEC:OTHER} we briefly present four quadratic programming
codes that were used against our method on five different problem
types, that are described in Section \ref{SEC:RESULTS}. Finally,
in Section \ref{SEC:TRUST} a new trust region like method is
proposed which takes full advantage of our quadratic programming
algorithm.

\section{Solving the quadratic problem \label{SEC:SOLQUAD}}
For the problem in Eq.~(\ref{QP}), we construct the associated
Lagrangian:
\begin{equation}
L(x,\lambda,\mu)= \frac{1}{2}x^TBx +x^Td
-\lambda^T(x-a)-\mu^T(b-x)
\end{equation}
The KKT necessary conditions at the minimum $\ds x^*, \lambda^*,
\mu^* $ require that:
\begin{eqnarray}
\label{KTP}
Bx^* +d -\lambda^*+\mu^*=0 \nonumber \\
\lambda^*_i \ge 0,\ \ \mu^*_i \ge 0, \ \forall i \in I \nonumber \\
\lambda^*_i(x^*_i-a_i) = 0, \ \forall i \in I \\
\mu^*_i(b_i-x^*_i) = 0, \ \forall i \in I \nonumber \\
x^*_i \in [a_i,b_i],   \ \forall i \in I \nonumber
\end{eqnarray}
A solution to the above system (\ref{KTP}), can be obtained
through an active set strategy described in detail in Algorithm~1:

\vskip 0.4cm \hrule \noindent \vskip 0.2cm {\bf Algorithm 1} \hskip 0.2cm
 BOXCQP \vskip 0.5mm \hrule

\label{ALG:BOXCQP} \vskip 0.4cm \noindent Initially set: $\ds k =
0$, $\ds \lambda^{(0)}=\mu^{(0)}=0$
and $\ds x^{(0)}= -B^{-1}d$. \\
{\bf If } $\ds  x^{(0)} $ is feasible, {\bf Stop},  the solution is: $\ds x^*=x^{(0)}$. \\
At iteration $k$, the quantities $x^{(k)}, \lambda^{(k)},
\mu^{(k)}$ are available.
\begin{enumerate}
\item Define the sets:
\begin{eqnarray*}
  L^{(k)}&=&\{i: x_i^{(k)}<a_i,\hbox{\ or \ }x_i^{(k)}=a_i
     \hbox{\ and \ } \lambda_i^{(k)}\ge 0\}\\
  U^{(k)}&=&\{i: x_i^{(k)}>b_i,\hbox{\ or \ }x_i^{(k)}=b_i
     \hbox{\ and \ } \mu_i^{(k)}\ge 0\}\\
  S^{(k)}&=&\{i: a_i<x_i^{(k)}<b_i,\hbox{\ or \ }x_i^{(k)}=a_i
     \hbox{\ and \ } \lambda_i^{(k)} < 0,\\
    &\hbox{\ or \ }&x_i^{(k)}=b_i \hbox{\ and \ } \mu_i^{(k)} <  0\}
\end{eqnarray*}
Note that $ L^{(k)} \cup U^{(k)} \cup  S^{(k)} = I $ \item Set:
\begin{eqnarray*}
&x_i^{(k+1)}&=\ a_i,\  \mu_i^{(k+1)}=0,\ \ \forall i \in L^{(k)} \\
 &x_i^{(k+1)}&=\ b_i,\  \lambda_i^{(k+1)}=0, \ \  \forall i \in U^{(k)} \\
 &\lambda_i^{(k+1)}&=\ 0,\ \mu_i^{(k+1)}=0, \ \  \forall i \in S^{(k)}
\end{eqnarray*}
\item
 Solve:
\begin{eqnarray*}
Bx^{(k+1)}+d = \lambda^{(k+1)}-\mu^{(k+1)}
\end{eqnarray*}
for the $ \ds N $ unknowns:
\begin{eqnarray*}
x_i^{(k+1)}, \  \forall i\  \in S^{(k)} \\
\mu_i^{(k+1)}, \  \forall i\  \in U^{(k)} \\
\lambda_i^{(k+1)}, \  \forall i\  \in L^{(k)}
\end{eqnarray*}
\item Check if the new point is a solution and decide to either
stop or iterate.
\begin{center}
\bt {\bf If} $\ds  ( x_i^{(k+1)}  \in [a_i,b_i]\ \forall i \in
S^{(k)}$ \= {\bf and}
$\ds  \mu_i^{(k+1)} \ge 0,\  \forall  i   \in U^{(k)} $ \\
\>{\bf and} $\ds  \lambda_i^{(k+1)} \ge 0,\  \forall  i   \in L^{(k)} )$ {\bf Then} \\
\hspace{1 cm} {\bf Stop}, the solution is: $\ds x^*=x^{(k+1)}$. \\
{\bf Else} \\
\hspace{1 cm} set $\ds  k \leftarrow k+1 $ and iterate from {\bf Step 1.} \\
{\bf Endif} \et
\end{center}
\end{enumerate}
 \hrule \vskip 0.4cm

The solution of the linear system in Step~3 of
Algorithm~\ref{ALG:BOXCQP}, needs further consideration. Let us
rewrite the system in a componentwise fashion.

\begin{equation}
\label{DETAIL} \sum_{j \in I}B_{ij}x^{(k+1)}_j +d_i =
\lambda^{(k+1)}_i - \mu^{(k+1)}_i,\ \forall i\in I
\end{equation}
Since $ \ds \forall i \in S^{(k)} $ we have that $ \ds
\lambda^{(k+1)}_i = \mu^{(k+1)}_i = 0$, hence  we can calculate $
\ds x^{(k+1)}_i,\  \forall i \in S^{(k)} $ by splitting the sum in
Eq.~(\ref{DETAIL}) and taking into account Step~2 of the
algorithm, i.e.:
\begin{equation}
\label{SPLITS} \sum_{j \in  S^{(k)}}B_{ij}x^{(k+1)}_j  = -\sum_{j
\in  L^{(k)}}B_{ij}a_j -\sum_{j \in  U^{(k)}}B_{ij}b_j -d_i, \
\forall i\in S^{(k)}
\end{equation}
The submatrix $ \ds B_{ij}, \hbox{ with } i,j \in  S^{(k)}$ is
positive definite as can be readily verified, given that the full
matrix $B$ is. The calculation of $ \ds \lambda^{(k+1)}_i, \
\forall i \in L^{(k)} $ and of $\ds  \mu^{(k+1)}_i, \ \forall i
\in U^{(k)} $ is straightforward and is given by:
\begin{equation}
 \lambda^{(k+1)}_i = \sum_{j \in I}B_{ij}x^{(k+1)}_j +d_i, \ \forall i\in  L^{(k)}
\end{equation}
\begin{equation}
 \mu^{(k+1)}_i = - \sum_{j \in I}B_{ij}x^{(k+1)}_j -d_i, \ \forall i\in  U^{(k)}
\end{equation}

The main computational task of the algorithm above, is the
solution of the linear system in Step~3. The size of the system
may vary according to the size of the active set in each
iteration. In our implementation we solve the linear system using
either the conjugate gradient method (Variant~1) or a direct
solver via $LDL^T$ decomposition (Variant~2). For large scale
problems the conjugate gradient method may be preferable.

\section{Other quadratic codes}
\label{SEC:OTHER}
\subsection{QPBOX}
QPBOX~\cite{QPBOX} is a  Fortran77 package for box constrained
quadratic programs developed in IMM in Technical University of
Denmark.  The bound constrained quadratic program is solved via a
dual problem, which is the minimization of an unbounded, piecewise
quadratic function. The dual problem involves a lower bound of
$\lambda_1$, i.e the smallest eigenvalue of a symmetric, positive
matrix, and is solved by Newton iteration with line search.

\subsection{QLD}
This code~\cite{QLD} is due to K.Schittkowski of the
University of Bayreuth, Germany and is a modification of routines
due to MJD Powell at the University of Cambridge. It is
essentially an active set, interior point method and supports
general linear constraints too.

\subsection{DUALQP}
This method was originally proposed by Goldfarb and Idnani \cite{GOLIDN83} for
positive definite quadratic programming problems. It takes
advantage of the fact that the unconstrained minimum of the
objective function can be used as a starting point (notice that we
choose the same starting point). DUALQP can be applied to problems
with general linear constraints also. The implementation is due to
P. Spellucci of the Technical University of
Darmstadt~\cite{DUALQP}.

\subsection{QUACAN}
This algorithm combines conjugate gradients with gradient
projection techniques, as the algorithm of Mor\'{e} and
Toraldo~\cite{MORE}. A new strategy for the decision of leaving
the current face is introduced, that makes possible to obtain
finite convergence even for a singular Hessian and in the presence
of dual degeneracy. QUACAN~\cite{QUACAN} is specialized for convex
problems subject to simple bounds.


\section{Experimental results}
\label{SEC:RESULTS} To verify the effectiveness of the proposed
approach we used five different problem types, and  measured cpu
times to make a comparison possible. We have implemented BOXCQP in
Fortan 77 and used an Intel Pentium--4 processor with Linux
operating system.

In the subsections below we describe in brief the different test
problems used for the experiments, and report our results.

\subsection{Random problems}
The first set of experiments includes randomly generated problems.
We create random positive--definite $B$ matrices, random bounds
$a$ and $b$ and we choose $d$ so that the unconstrained minimum
falls outside the box. (See Table 1).
%
% Table random
%

\subsection{Circus Tent problem}
The circus tent problem is taken from Matlab's optimization demo
as an example of large-scale quadratic programming with simple
bounds. The problem is to build a \emph{circus tent} to cover a
square lot. The tent is elastic and is to be supported by five
poles. The question is to find the shape of the tent at
equilibrium, that corresponds to the minimum of the energy
function. As we can see in Figure~\ref{FIG:CIRCUS}, the problem
has only lower bounds imposed by the five poles and the ground.

%
% figure circus
%
The surface formed by the elastic tent, is determined by solving
the bound constrained optimization problem:

\begin{eqnarray}
\min_x f(x) = \frac{1}{2}x^THx + x^Tc \\
\textrm{subject to: } l \leq x \nonumber \hskip 0.85cm
\end{eqnarray}

where $f(x)$ corresponds  to  the energy function and $H$ is a
5-point finite difference Laplacian over a square grid.

%
% Table circus
%

\subsection{Biharmonic Equation problem}

We consider the problem of describing small vertical deformations
of an horizontal, elastic membrane clamped on a rectangular
boundary, under the influence of a vertical force.  The membrane
is constrained to remain below an obstacle. For an in depth
discussion of this problem see~\cite{OIALLOI}. The formulation of
the problem is given by Eq.(~\ref{EQN:BIH}).

\begin{eqnarray}
\label{EQN:BIH}
\min_u \frac{1}{2}u^TQu + u^Tf \\
\textrm{subject to: } u \leq \psi \nonumber
\end{eqnarray}

We see  an example in Fig~\ref{FIG:BIHARM} of a membrane under the
influence of a vertical force.

%
% Table biharmo
%

%
% Figure biharmo
%
\subsection{Intensity Modulated Radiation Therapy}

This problem arises in the field of radiotherapy and concerns the
determination of the spatial distribution of the radiation, in a
way that the patient's vital organs are minimally irradiated.
Knowing the beam settings and the intensity profile, one can
calculate the radiation dose. Inversely, when a desired dose is
required, the proper intensity profile for given beam settings can
be retrieved by solving a quadratic problem. The beam settings are
successively modified in an effort to satisfy a set of clinical
constraints, and hence the quadratic subproblem (shown in Eq.
(\ref{EQ:DOSE})), must be solved  a large number of times
\cite{DOSE}.

\begin{eqnarray} \label{EQ:DOSE}
\displaystyle
 \min_f s(f) = \frac{1}{2}f^T A f + f^T b   \\
\nonumber
 \textrm{subject to } f \ge 0 \hskip 0.85cm
\end{eqnarray}

The results reported in Table~\ref{TAB:IMRT}, correspond to real
world data, kindly provided by S. Breedveld \cite{BREED}.
In this example, seven beams are combined resulting to a quadratic problem
with $2342$ parameters.

%
% Table IMRT
%

\subsection{Support Vector Classification}

In this classification problem, the goal is to separate  two
classes using a hyperplane $f(x) = w^Tx + b$,  which is determined
from available examples ($D =\{(x^1, y^1), (x^2, y^2), \ldots
(x^l, y^l)\},\  x \in R^n,\ y \in{-1, 1} $). Furthermore it is
desirable to produce a classifier that will work well on unseen
examples, i.e. it will generalize well. Consider the example in
Fig.~\ref{FIG:SEP}. There are many possible linear classifiers
that can separate the data, but there is only one that maximizes
the distance to the nearest data point of each class. This
classifier is termed the optimal separating hyperplane and
intuitively, one would expect that generalizes optimally.

%
% Fig svm1
%
The formulation of the maximum distance linear classifier (if we
omit the constant term $b$ of the hyperplane
equation\footnote{Also known as explicit bias.}) is a convex
quadratic problem with simple bounds on the variables. The
resulting problem has the form:
\begin{eqnarray}
 \min_a \frac{1}{2}a^TQa - a^Te \\
 \textrm{subject to: } 0 \leq a_i \leq C \nonumber
 \label{EQN:SVMOPTIM}
\end{eqnarray}
\noindent
where $e \in R^l$ and with $e_i = 1$, \ $Q_{ij} = y^i y^j
K(x^i,x^j)$ and $K(x, y)$ is the kernel function performing the
non-linear mapping into the feature space. The parameters $a \in
R^l$ are Lagrange multipliers of an original quadratic problem,
that define the separating hyperplane using the relation:
\begin{equation}
{w^*}^Tx = \sum_{i=1}^{l} a^*_i y^i K(x^i, x)
\end{equation}
Hence the separating surface is given by:
\begin{equation}
f(x) =  \textrm{sgn}({w^*}^Tx)
\end{equation}

In our experiments we used the CLOUDS \cite{UCI} data set, which
is a two-dimensional data set with two classes. We have
constructed the problem in Eq. (\ref{EQN:SVMOPTIM}) using an RBF
Kernel function $K(x, y) = exp(\frac{-||x-y||^2} {2^2} )$, and
setting $C = 100$. The experiments conducted follow the procedure:
\begin{itemize}
\item Form the training set by extracting l examples from the
dataset and let the rest examples (5000-l) form the test set.
\item Construct the matrix Q for the problem in Eq. (\ref{EQN:SVMOPTIM})
\item  Apply each solver, obtain the corresponding separating
surface and test-set error. \end{itemize}

In these experiments the large condition number of matrix $Q$
leads to ill conditioned problems. To circumvent this, we added in
the main diagonal of $Q$ a small positive term. The resulting
classification surfaces for $l = 200$, $500$, $1000$ and $2000$
training examples from CLOUDS dataset are shown in
Fig.~\ref{FIG:SVM}.

%
% Table SVMs
%


\section{Bound Constrained Nonlinear Optimization }
\label{SEC:TRUST}
We present here a trust region method for non-linear optimization
with bound constraints, where the trust region is a hyperbox, in
contrast with the usual hypersphere or hyperellipsoid shapes.
The rectangular trust region is natural for problems with bound
constraints, because even when it overlaps with the feasible region,
its geometry is preserved.
Trust region methods fall in the category of sequential quadratic
programming. These algorithms are iterative and the objective function $f(x)$
(assumed to be twice continuously differentiable), is approximated in a proper
neighborhood of the current iterate (the trust region), by a quadratic model.
Namely, at the $k^{th}$ iteration the model is given by:
\begin{equation}
f(x^{k}+s) \approx  m^{(k)}(s) = f(x^{(k)})+ s^Tg^{(k)}
+\frac{1}{2}s^TB^{(k)}s \label{eq:quadmod}
\end{equation}
where $g^{(k)} = \nabla f(x^{(k)})$ and $B^{(k)}$ in the case of Newton's method
is a positive definite modification of the Hessian, while in the case of quasi-Newton
methods is a positive definite matrix produced by the relevant update.

The trust region may be defined by:
\begin{equation}
\mathbf{T}^{(k)} = \{ x \in \Re^n \ | \ ||x-x^{(k)}|| \leq
\Delta^{(k)} \}
\end{equation}
It is obvious that different choices for the norm lead to
different trust region shapes. The Euclidean norm $||\cdot||_2$,
corresponds to a hypershpere, while the $||\cdot||_\infty$ norm
defines a hyperbox.

Given the model and the trust region, we seek a  step
$||s^{(k)}||\leq\Delta^{(k)}$, that minimizes $m^{(k)}(s)$.
We compare the actual reduction $\delta f^{(k)}=f(x^{(k)})- f(x^{(k)}+s^{(k)})$,
to the model reduction  $\delta m^{(k)} = m^{(k)}(0)-m^{(k)}(s^{(k)})$.
If they agree to a certain extend, the step is accepted and the trust
region is either expanded or remains the same. Otherwise the step
is rejected and the trust region is contracted. The basic trust
region algorithm is sketched in Algorithm~2.


\vskip 0.4cm \hrule \noindent \vskip 0.2cm  {\bf Algorithm~2} \hskip 0.2cm
 Basic trust region  \vskip 0.5mm \hrule
%\begin{algorithm}{}

\begin{itemize}
 \item [1.] Pick the initial point and trust region parameter $x^{(0)}$
and $\Delta^{(0)}$, and set $k=0$.
 \item [2.] Construct a quadratic model:
$$
f(x^{k}+s) \approx  m^{(k)}(s) = f(x^{(k)})+ s^Tg^{(k)}
+\frac{1}{2}s^TB^{(k)}s
$$
 \item [3.]  Minimize $m^{(k)}(s)$ and hence determine  $||s^{(k)}|| \leq \Delta^{(k)} $
 \item [4.] Compute the ratio of actual to expected reduction:
$r^{(k)} = \frac{\delta f^{(k)}}{\delta m^{(k)}}$,
and  update the trust region, following the strategy of  Dennis \& Schnabel
\cite{DENSCH} (Appendix A, page 338).
 \item [5.] Increment $k \leftarrow k+1$ and repeat from 1.
\end{itemize}
\hrule \vskip 0.2cm
Consider the bound constrained problem:
\begin{equation}
 \hspace{1cm} \min_x f(x),  \hbox{\ \ }
 \hbox{subject to: }  l_i \leq x_i \leq u_i \nonumber
\end{equation}
(The unconstrained case is obtained by letting  $u_i=-l_i \rightarrow \infty$.)

\noindent
Let $x^{(k)}$ be the $k$-th iterate of the trust region algorithm. \\
Hence step 3 of Algorithm~2 becomes:
\begin{eqnarray}
 & & \hspace{2cm} \min_s \:m^{(k)}(s) =  s^Tg{^{(k)}} +  \frac{1}{2}s^T B^{(k)} s  \\
& & \hbox{ subject to: } \max(l_i -x^{(k)}_i ,  -\Delta^{(k)}) \leq s_i \leq \min(u_i - x^{(k)}_i,  \Delta^{(k)}) \nonumber
\end{eqnarray}

In the unconstrained case, our experiments (that used a BFGS update),
showed similar performance to spherical trust region implementations.
For the bound constrained case our method is obviously superior, since
it maintains the simplicity of the rectangular trust region, where
our efficient quadratic solver is applicable \cite{TRUSTQP}.
Note that the trust region formed by the intersection of a sphere
with the rectangular box defined by the parameter bounds, is not easy to treat.
Concluding, further numerical tests showed that our method performs similarly
to active set methods with line--search.
\section{Conclusions}
We have presented an active set algorithm for solving bound
constrained convex quadratic problems, using a Lagrange multiplier
approach that modifies, at each iteration, both the primal
and the dual variables. Extensive experimental testing showed significant
improvement over several known quadratic programming codes.
In addition a trust region method for nonlinear objectives has emerged,
that takes advantage of the proposed algorithm in order to efficiently
solve unconstrained and bound constrained problems.

\begin{thebibliography}{25}
\bibitem{ACT1} R. Fletcher, {\it Practical Methods of Optimization}, 2nd ed., John Wiley \&Sons, Inc., New York
(1987).
\bibitem{ACT2} P. E. Gill, W. Murray, and M. H. Wright, {\it Practical Optimization}, Academic Press, New York
(1981).
\bibitem{TR1} A. R. Conn, N. I. M. Gould and Ph. L. Toint, {\it Testing a class
of methods for solving minimization problems with simple bounds on
the variables}, Mathematics of Computation, {\bf 50}, pp. 399-430,
(1988).
\bibitem{MORE} J.J. More',  G. Toraldo, On the solution of quadratic programming
problems with bound constraints - SIAM J. on Opt. 1,  (1991),
pp.93-113.
\bibitem{GP} A. R. Conn, N. I. M. Gould, and Ph. L. Toint, {\it LANCELOT: A Fortran package for large-scale nonlinear
 optimization (Release A)}, Springer Series in Computational Mathematics{\bf 17},
Springer--Verlag (1992).
\bibitem{GP2} D. P. Bertsekas, {\it Nonlinear Programming}, Athena Scientific (1996)
\bibitem{SVM1} S. Gunn, {\it Support vector machines for classification and
regression}, ISIS technical report, Image Speech \& Intelligent
Systems Group, University of Southampton (1997).
\bibitem{SVM2}E. Osuna, R. Freund, and F. Girosi, {\it Support vector machines: Training and
applications},  A.I. Memo (in press) 1602, MIT A. I. Lab. (1997).
\bibitem{QPOPT} P.E. Gill, W. Murray, M.A. Saunders, and M.H. Wright,
{\it Inertia-controlling methods for general quadratic
programming} , SIAM Rev. 33 (1991), pp. 1--36
\bibitem{QLD} K. Schittkowski, {\it QLD: A FORTRAN Code for Quadratic
Programming, User’s Guide}, Mathematisches Institut,
Universit\"{a}t Bayreuth, Germany (1986).
\bibitem{QPBOX}K. Madsen, H. B. Nielsen,  and M.C. Pinar, {\it Bound Constrained
Quadratic Programming via Piecewise Quadratic Functions},
Mathematical Programming {\bf 85} (1) (1999)
\bibitem{GOLIDN83} D. Goldfarb and A. Idnani,
                {\it A Numerically Stable Dual Method for Solving Strictly
                Convex Quadratic Programs},
                Mathematical Programming, {\bf 27} (1983) 1--33.
\bibitem{UCI}C. L. Blake, and C. J. Merz, {\it  UCI Repository of machine learning
databases},
\url{http://www.ics.uci.edu/~mlearn/MLRepository.html}, Irvine,
University of California, Department of Information and Computer
Science (1998).
\bibitem{PORT} A. F. Perold, {\it Large-Scale Portfolio Optimization}, Management Science {\bf 30}, 1143-1160 (1984).
\bibitem{TRUSTQP} C. Voglis and I. E. Lagaris, {\it A Rectangular Trust-Region
Approach for Unconstrained and Box-Constrained Optimization
Problems}, International Conference of Computational Methods in
Sciences and Engineering (ICCMSE 2004), Athens, (2004)
\bibitem{QUACAN} A. Friedlander and M. Martinez, On the maximization of a concave
quadratic function with box constraints, SIAM J. Optimization 4,
pp. 177-192 (1994).
\bibitem{DUALQP} \url{http://www.mathematik.tu-darmstadt.de:8080/ags/ag8/Mitglieder/spellucci_de.html}
\bibitem{DOSE} S.~Breedveld, P.~R.~M.~Storchi, M.~Keijzer and B.~J.~M.~Heijmen {\it Fast, multiple optimizations of quadratic dose
objective functions in IMRT}
\bibitem{BREED} S.~Breedveld, private communication.
\bibitem{DENSCH} J. E. Dennis, R. B. Schnabel {\it Numerical Methods
for Unconstrained Optimization and Nonlinear Equations}, SIAM (1996).
\bibitem{OIALLOI} K. Kunisch, and F. Rendl, {\it An Infeasible Active Set Method For Quadratic Problems With Simple Bounds},
SIAM J. Optimization 14(1), pp. 35-52 (2003)
\end{thebibliography}

\vfill \eject
\begin{table}[!h]
\centering
\begin{tabular}{|c|c|c|c|c|c|}
\hline
 Prob. Size & BOXCQP & QPBOX & QLD & DUALQP & QUACAN \\ \hline \hline
 100 & 0.41& 1.22& 0.61& 0.91& 1.33 \\ \hline
 200 & 0.48& 1.53& 0.83& 1.06& 1.63 \\ \hline
 300 & 0.67& 1.74& 0.91& 1.12& 2.05 \\ \hline
 500 & 0.8& 1.93& 1.32& 1.67&  2.22\\ \hline
 1000 & 1.1& 2.76& 2.12& 2.33&  3.11\\ \hline
 2000 & 1.8& 3.17& 2.37& 3.13&  4.12\\ \hline
\end{tabular}
\caption{CPU times (secs) for the Random test problems.}
\end{table}

\vfill \eject

\begin{table}[!h]
\centering
\begin{tabular}{|c|c|c|c|c|c|}
\hline
 Prob. Size & BOXCQP & QPBOX & QLD & DUALQP & QUACAN \\ \hline \hline
 400 & 0.02& 0.48& 0.20& 0.38& 0.05 \\ \hline
 625 & 0.03& 0.69& 0.79& 0.66& 0.12 \\ \hline
 900 & 0.04 & 0.77 & 2.30& 1.11&0.35  \\ \hline
 1225 & 0.05& 0.93& 5.57& 1.61& 0.69  \\ \hline
 1600 & 0.07 & 1.21& 6.32 & 1.87 & 0.74 \\ \hline
 2025 & 0.10& 1.56 & 8.11 & 2.16 & 0.89 \\ \hline
 2500 & 0.12 & 1.78 & 9.12 & 2.51 & 1.02 \\ \hline
 3600 & 0.17 & 2.01 & 11.89 & 3.24 &  1.23\\ \hline
\end{tabular}
\caption{CPU times (secs) for the Circus Tent problem.}
\end{table}

\vfill \eject

\begin{table}[!h]
\centering
\begin{tabular}{|c|c|c|c|c|c|}
\hline
 Prob. Size & BOXCQP & QPBOX & QLD & DUALQP & QUACAN \\ \hline \hline
  400 & 0.05& 0.07 & 0.56 & 0.38 & 0.18 \\ \hline
  625 & 0.42 & 0.65 & 1.13& 1.71 & 1.12 \\ \hline
  900 & 0.81 & 1.07 & 1.29  & 2.23 & 1.65  \\ \hline
  1225 & 1.37 & 2.47 & 3.01 & 3.11 & 2.10  \\ \hline
 1600 & 2.18 & 2.98 & 3.87 & 4.21 & 2.76  \\ \hline
 2025 & 2.67 & 3.34 & 4.42 & 4.88 & 3.64 \\ \hline
 2500 & 3.52& 4.07 & 5.10 & 5.43 &  4.41 \\ \hline
 3600 & 4.67 & 5.33 & 6.02 & 6.23 & 5.53 \\ \hline
\end{tabular}
\caption{CPU times (secs) for the Biharmonic Equation problem.}
\end{table}

\vfill \eject
\begin{table}[!h]
\centering
\begin{tabular}{|c|c|c|c|c|c|}
\hline
 Prob. Size & BOXCQP & QPBOX & QLD & DUALQP & QUACAN \\ \hline \hline
  2342 & 1.91 & 2.12 & 1.78 & 1.56 & 1.67 \\ \hline
\end{tabular}
\caption{CPU times (secs) for the IMRT problem.} \label{TAB:IMRT}
\end{table}

\vfill \eject
\begin{table}[!h]
\centering
\begin{tabular}{|c|c|c|c|c|c|}
\hline
 Prob. Size & BOXCQP & QPBOX & QLD & DUALQP & QUACAN \\ \hline \hline
  200 & 0.34 &  0.51  &  0.04 &  0.08&  0.12\\ \hline
  500 & 2.23 &  2.55 &  4.22 &   3.21&  7.86\\ \hline
  1000 & 3.45 &  4.12 &  10.81 &  9.76 & 10.12 \\ \hline
  2000 & 15.66 &  18.63 &  20.51 &  29.01&  25.11\\ \hline
  3000 & 24.23 &  29.91 &   33.62 &  52.11 & 31.49 \\ \hline
  4000 & 36.77 &   43.63 &  50.04  &  69.94 & 45.22 \\ \hline
\end{tabular}
\caption{CPU times (secs) for the SVM training problem.}
\label{TAB:SVM}
\end{table}

\vfill \eject
\begin{figure}[!h]
\centering \includegraphics[width=5.5in,height=3in]{circus}
\caption{Circus tent problem}\label{FIG:CIRCUS}
\end{figure}

\vfill \eject
\begin{figure}[!h]
\centering
\includegraphics[scale=2.3]{biharm} \caption{The force is presented on the left, and
the deformation on the right} \label{FIG:BIHARM}
\end{figure}

\vfill \eject
\begin{figure}[!h]
\centering \includegraphics[scale=0.7]{optimal} \caption{Maximum
distance classifier} \label{FIG:SEP}
\end{figure}

\vfill \eject
\begin{figure}[!h]
{\centering \subfigure[200 Training
Examples]{\includegraphics[scale=0.45]{svm1}}\subfigure[500
Training Examples]{\includegraphics[scale=0.45]{svm2}} }
{\centering \subfigure[1000 Training
Examples]{\includegraphics[scale=0.45]{svm3}}\subfigure[2000
Training Examples]{\includegraphics[scale=0.45]{svm4}} }
\caption{SVM classification surfaces}\label{FIG:SVM}
\end{figure}

\end{document}
