Large-Scale Active-Set Box-Constrained Optimization Method with Spectral Projected Gradients.
A new active-set method for smooth box-constrained minimization is introduced. The algorithm combines an unconstrained method, including a new line-search which aims to add many constraints to the working set at a single iteration, with a recently introduced technique (spectral projected gradient) for dropping constraints from the working set. Global convergence is proved. A computer implementation is fully described and a numerical comparison assesses the reliability of the new algorithm.
Introduction
The problem considered in this paper consists in the minimization of a
smooth with bounds on the variables. The feasi-
Department of Computer Science IME-USP, University of S~ao Paulo, Rua do Mat~ao
1010, Cidade Universitaria, 05508-090, S~ao Paulo SP, Brazil. This author was supported
by PRONEX-Optimization 76.79.1008-00, FAPESP (Grants 99/08029-9 and 01/04597-4)
and CNPq (Grant 300151/00-4). e-mail: egbirgin@ime.usp.br
y Department of Applied Mathematics IMECC-UNICAMP, University of Campinas,
This author was supported by PRONEX-
Optimization 76.79.1008-00, FAPESP (Grant 01/04597-4), CNPq and FAEP-UNICAMP.
e-mail: martinez@ime.unicamp.br
ble
set
is dened
by
Box-constrained minimization algorithms are used as subalgorithms for solving
the subproblems that appear in many augmented Lagrangian and penalty
methods for general constrained optimization. See [11, 12, 16, 17, 18, 19,
20, 21, 26, 28, 31]. A very promising novel application is the reformulation
of equilibrium problems. See [1] and references therein. The methods introduced
in [11] and [26] are of trust-region type. For each iterate x k 2
a quadratic approximation of f is minimized in a trust-region box. If the
objective function value at the trial point is su-ciently smaller than f(x k ),
the trial point is accepted. Otherwise, the trust region is reduced. The
dierence between [11] and [26] is that, in [11], the trial point is in the face
dened by a \Cauchy point", whereas in [26] the trial point is obtained
by means of a specic box-constrained quadratic solver, called QUACAN.
See [2, 15, 23, 25, 31] and [13] (p. 459). Other trust-region methods for
box-constrained optimization have been introduced in [3, 29].
QUACAN is an active-set method that uses conjugate gradients within
the faces, approximate internal-face minimizations, projections to add constraints
to the active set and an \orthogonal-to-the-face" direction to leave
the current face when an approximate minimizer in the face is met. In [17]
a clever physical interpretation for this direction was given.
Numerical experiments in [16] suggested that the e-ciency of the algorithm
[26] relies, not in the trust-region strategy, but in the strategy of
QUACAN for dealing with constraints. This motivated us to adapt the
strategy of QUACAN to general box-constrained problems. Such adaptation
involves two main decisions. On one hand, one needs to choose an
unconstrained minimization algorithm to deal with the objective function
within the faces. On the other hand, it is necessary to dene robust and
e-cient strategies to leave faces and to add active constraints. Attempts for
the rst decision have been made in [6] and [10]. In [10] a secant multipoint
minimization algorithm is used and in [6] the authors use the second-order
minimization algorithm of Zhang and Xu [36].
In this paper we adopt the leaving-face criterion of [6], that employs the
spectral projected gradients dened in [7, 8]. See, also, [4, 5, 32, 33, 34].
For the internal minimization in the faces we introduce a new general algorithm
with a line search that combines backtracking and extrapolation.
The compromise in every line-search algorithm is between accuracy in the
localization of the one-dimensional minimizer and economy in terms of functional
evaluations. Backtracking-like line-search algorithms are cheap but,
sometimes, tend to generate excessively small steps. For this reason, back-tracking
is complemented with a simple extrapolation procedure here. The
direction chosen at each step is arbitrary, provided that an angle condition
is satised.
In the implementation described in this paper, we suggest to choose
the direction using the truncated-Newton approach. This means that the
search vector is an approximate minimizer of the quadratic approximation
of the function in the current face. We use conjugate gradients to nd
this direction, so the rst iterate is obviously a descent direction, and this
property is easily monitorized through successive conjugate gradient steps.
The present research is organized as follows. In Section 2 we describe an
\unconstrained" minimization algorithm that deals with the minimization
of a function on a box. The algorithm uses the new line-search technique.
Due to this technique it is possible to prove that either the method nishes
at a point on the boundary where, perhaps, many constraints are added,
or it converges to a point in the box where the gradient vanishes. The
box-constrained algorithm is described in Section 3. Essentially, we use
the algorithm of Section 2 to work within the \current face" and spectral
projected gradients [7] to leave constraints. The spectral projected gradient
technique also allows one to leave many bounds and to add many others to
the working set at a single iteration. This feature can be very important
for large-scale calculations. In this section we prove the global convergence
of the box-constrained algorithm. The computational description of the
code (GENCAN) is given in Section 4. In Section 5 we show numerical
experiments using the CUTE collection. In Section 6 we report experiments
using some very large problems (up to 10 7 variables). Finally, in Section 6
we make nal comments and suggest some lines for future research.
In this section we assume that f : IR
ug:
The set B will represent each of the closed faces
of
in Section 3. The
dimension  n in this section is the dimension of the reduced subspace of the
Section 3 and the gradient of this section is composed by the derivatives
with respect to free variables in Section 3. We hope that using the notation
rf in this section will not lead to confusion.
From now on, we denote
Our objective is to dene a general iterative algorithm that starts in the
interior of B and, either converges to an unconstrained stationary point, or
nishes in the boundary of B having decreased the functional value. This
will be the algorithm used \within the faces" in the box-constrained method.
Algorithm 2.1 is based on line searches with Armijo-like conditions and
extrapolation. Given the current point x k and a descent direction d k , we
nish the line search if x k +d k satises a su-cient descent criterion and if the
directional derivative is su-ciently larger than hg(x k ); d k i. If the su-cient
descent criterion does not hold, we do backtracking. If we obtained su-cient
descent but the increase of the directional derivative is not enough, we try
extrapolation.
Let us explain why we think that this philosophy is adequate for large-scale
box-constrained optimization.
1. Pure backtracking is enough for proving global convergence of many
optimization algorithms. However, to accept the rst trial point when
it satises an Armijo condition can lead to very small steps in critical
situations. Therefore, steps larger than the unity must be tried when
some indicator says that this is worthwhile.
2. If the directional derivative su-ciently larger than
we consider that there is not much to decrease increasing
the steplength in the direction of d k and, so, we accept the unit
steplength provided it satises the Armijo condition. This is reasonable
since, usually, the search direction contains some amount of
second-order information that makes the unitary steplength desirable
from the point of view of preserving a satisfactory order of convergence.
3. If the unitary steplength does not satisfy the Armijo condition, we
do backtracking. In this case we judge that it is not worthwhile to
compute gradients of the new trial points, which would be discarded
if the point is not accepted.
4. Extrapolation is especially useful in large-scale problems, where it is
important to try to add as many constraints as possible to the working
set. So, we extrapolate in a rather greedy way, multiplying the
steplength by a xed factor while the function value decreases.
5. We think that the algorithm presented here is the most simple way
in which extrapolation devices can be introduced with a reasonable
balance between cost and e-ciency. It is important to stress that this
line search can be coupled with virtually any minimization procedure
that computes descent directions.
For all z 2 IR n , the Euclidean projection of z onto a convex set S will be
denoted P S (z). In this section, we denote P (y). The symbol k  k
represents the Euclidean norm throughout the paper.
Algorithm 2.1: Line-search based algorithm
The algorithm starts with x 0 2 Int(B). The non-dimensional parameters
are given.
We also use the small tolerances  abs ;  rel > 0. Initially, we set k 0.
Step 1. Computing the search direction
Step 1.1 If kg k
Step 1.2 Compute such that
Step 2. Line-search decisions
Step 2.1 Compute
set  minf
then go to Step 2.2
else go to Step 2.3.
Step 2.2 (At this point we have x k
If
take and go to Step 5
else go to Step 3 (Extrapolation)
else go to Step 4 (Backtracking).
Step 2.3 (At this point we have x k
If
take  k   max and x
such that f(x k+1 )  and go to Step 5
(In practice, such a point is obtained performing Step 3
of this algorithm (Extrapolation).)
else go to Step 4 (Backtracking).
Step 3. Extrapolation
Step 3.1 If ( <  max and N >  max ) then set  trial  max
else set  trial N.
Step 3.2 If (   max and kP
take
the execution of Algorithm 2.1.
Step 3.3 If (f(P
take to Step 5
else set   trial and go to Step 3.1.
Step 4. Backtracking
Step 4.1 Compute new .
Step 4.2 If (f(x k
take and go to Step 5
else go to Step 4.1.
Step 5. If  k   max terminate the execution of Algorithm 2.1
else set to Step 1.
Remarks. Let us explain here the main steps of Algorithm 2.1 and their
motivations. The algorithm perform line-searches along directions that satisfy
the angle-cosine condition (2). In general, this line search will be used
with directions that possess some second-order information, so that the \nat-
ural" step must be initially tested and accepted if su-cient-descent
and directional-derivative conditions ((3) and (4)) are satised.
The rst test, at Step 2.1, asks whether x k is interior to the box. If
this is not the case, but f(x k we try to obtain smaller
functional values multiplying the step by a xed factor and projecting onto
the box. This procedure is called \Extrapolation". If x k is not interior
and backtracking.
is interior but the Armijo condition (3) does not hold, we
also do backtracking. Backtracking stops when the Armijo condition (6)
is fullled. If (3) holds, we test the directional derivative condition (4).
As we mentioned above, if (4) is satised too, we accept x k new
point. However, if (3) holds and (4) does not, we judge that, very likely,
taking larger steps along the direction d k will produce further decrease of
the objective function. So, in this case we also do Extrapolation.
In the Extrapolation procedure we try successive projections of x k +d k
onto the box, with increasing values of . If the entry point
interior but x k +Nd k is not, we make sure that the point x will be
tested rst. The extrapolation nishes when decrease of the function is not
obtained anymore or when the distance between two consecutive projected
trial points is negligible.
The iteration of Algorithm 2.1 nishes at Step 5. If the corresponding
iterate x k+1 is on the boundary of B, the algorithm stops, having encountered
a boundary point where the functional value decreased with respect
to all previous ones. If x k+1 is in the interior of B the execution continues
increasing the iteration number.
The
ux-diagrams in Figures 1 and 2 help to understand the structure
of the line-search procedure.
x+d -Int aa-amax
Line Search
bb-condition
xnew-x+d
End
Armijo
Backtracking
Extrapolation

Figure

1: Line Search procedure.
Extrapolation
a - aamax and
dist < ee
aa-aatrial
aatrial-NNaa
a < aamax and
NNa > aamax
f(P(x+aatrial d))
xnew-P(x+aad)
End

Figure

2: Extrapolation strategy.
In the following theorem we prove that any sequence generated by Algorithm
2.1, either stops at an unconstrained stationary point, or stops in the
boundary of B, or generates, in the limit, unconstrained stationary points.
Theorem 2.1. Algorithm 2.1 is well dened and generates points with
strictly decreasing functional values. If fx k g is a sequence generated by
Algorithm 2.1, one of the following possibilities holds.
(i) The sequence stops at x k , with g(x k
(ii) The sequence stops at x
(iii) The sequence is innite, it has at least one limit point, and every limit
point x  satises g(x
Proof. Let us prove rst that the algorithm is well dened and that it
generates a sequence with strictly decreasing function values. To see that
it is well dened we prove that the loops of Steps 3 and 4 necessarily nish
in nite time. In fact, at Step 3 we multiply the nonnull direction d k by a
number greater than one, or we take the maximum allowable feasible step.
Therefore, eventually, the boundary is reached or the increase condition (5)
is met. The loop of Step 4 is a classical backtracking loop and nishes
because of well-known directional derivative arguments. See [14]. On exit,
the algorithm always requires that f(x k so the sequence
strictly decreasing.
It remains to prove that, if neither (i) nor (ii) hold, then any cluster
point x  of the generated sequence satises g(x  be an innite
subset of IN such that
lim
Suppose rst that ks k k is bounded
away from zero for k 2 K 1 . Therefore, there exists  > 0 such that ks k k
for all k 2 K 1 . By (3) or (6) we have that
for all k 2 K 1 . Therefore, by (2),
By the continuity of f this implies that lim k2K 1
Suppose now that ks k k is not bounded away from zero for k 2 K 1 . So,
there exists K 2 , an innite subset of K 1 , such that lim k2K 2
ks
be the set of indices such that  k is computed at Step 2.2 for
Analogously, let K 4  K 2 be the set of indices such that  k is
computed at Step 3 for all k 2 K 4 and let K 5  K 2 be the set of indices such
that  k is computed at Step 4 for all k 2 K 5 . We consider three possibilities:
(i) K 3 is innite.
(ii) K 4 is innite.
(iii) K 5 is innite.
Consider, rst, the case (i). By (4) we have that
and
ks
ks
for all k 2 K 3 . Since K 3 is innite, taking an convergent subsequence
taking limits in (7) and using continuity, we obtain that
Since  2 (0; 1), this implies that hg(x  ); di  0. But, by (2) and continuity,
Consider, now, Case (iii). In this case, K 5 is innite. For all k 2 K 5
there exists s 0
k such that
and
ks 0
ks k
By (10), lim ks 0
by (9), we have for all k 2 K 5 ,
So, by the Mean-Value theorem, there exists  k 2 [0; 1] such that
for all k 2 K 5 . Dividing by ks 0
k k, taking limits for a convergent subsequence
d) we obtain that
This inequality is similar to (8). So, g(x  follows from the same
arguments.
Consider, now, Case (ii). Since we are considering cases where an innite
sequence is generated it turns out that, in (5), P
Moreover, by Step 3.1,  trial  N and P
Therefore, for all k 2 K 4 , writing  0
we have that  0
and
Therefore, by the Mean-Value theorem, for all k 2 K 4 there exists  k 2
k ] such that
Thus, for all k 2 K 4 , since  0
we have that
dividing by kd k k and taking a convergent subsequence of d k =kd k k, we obtain:
hg(x  ); di  0:
But, by (2), taking limits we get hg(x  ); di  kg(x  )k. This implies that
This completes the proof. 2
3 The box-constrained algorithm
The problem considered in this section is
Minimize f(x) subject to x
where
is given by (1).
As in [23], let us divide the feasible
set
into disjoint open faces, as
follows. For all I
We I the smallest a-ne subspace that contains F I and S I the
parallel linear subspace to V I . The (continuous) projected gradient at xis dened as
For all x 2 F I , we dene
I
[g P (x)]:
The main algorithm considered in this paper is described below.
Algorithm 3.1: GENCAN
Assume that x 0is an arbitrary initial point,  2 (0; 1) and 0 <  min
I be the face that contains the current iterate x k . Assume
that g P (otherwise the algorithm terminates). At the main iteration
of the algorithm we perform the test
kg I
If (13) takes place, we judge that it is convenient that the new iterate belongs
to
F I (the closure of F I ) and, so, we compute x k+1 doing one iteration of
Algorithm 2.1, with the set of variables restricted to the free variables in F I .
So, the set B of the previous section corresponds to
F I here.
If (13) does not hold, we decide that some constraints should be abandoned
and, so, the new iterate x k+1 is computed doing one iteration of the SPG
method described by Algorithm 3.2. In this case, before the computation of
x k+1 we compute the spectral gradient coe-cient  k in the following way.
Otherwise, dene
and
Algorithm 3.2 is the algorithm used when it is necessary to leave the
current face, according to the test (13).
Algorithm 3.2: SPG
Compute as the next iterate of a monotone SPG iteration [7, 8] with
the spectral step  k . Namely, we dene the search direction d k as
and we compute x in such a way that
trying rst perhaps, reducing this coe-cient by means of a
safeguarded quadratic interpolation procedure.
Remark. Observe that x
F I if x k 2 F I and x k+1 is computed by Algorithm
3.2. In this case, (13) does not hold, so kg P
the components corresponding to the free variables of g I
are the same, this means that g P components corresponding
to xed variables. Therefore,
F I for all  > 0. So,
F I for all  > 0. But, according to the SPG iteration,
for some  > 0,  0 > 0. This implies that x
F I .
We nish this section giving some theoretical results. Roughly speaking,
we prove that the algorithm is well dened and that a Karush-Kuhn-Tucker
point is computed up to an arbitrary precision. Moreover, under dual-
nondegeneracy, the (innite) algorithm identies the face to which the limit
belongs in a nite number of iterations.
Theorem 3.1. Algorithm 3.1 is well dened.
Proof. This is a trivial consequence of the fact that Algorithm 2.1 and Algorithm
3.2 (the SPG algorithm [7]) are well dened. 2
Theorem 3.2. Assume that fx k g is generated by Algorithm 3.1. Suppose
that there exists  k 2 f0; I for all k
k. Then,
every limit point of fx k g is rst-order stationary.
Proof. In this case, x k+1 is computed by Algorithm 2.1 for all k   k. Thus,
by Theorem 2.1, the gradient with respect to the free variables tends to zero.
By a straightforward projection argument, it follows that kg I
Since (13) holds, this implies that kg P every limit point is
rst-order stationary. 2
Theorem 3.3. Suppose that for all k 2 f0; I , there exists
such that x k
I . Then, there exists a limit point of fx k g that is
rst-order stationary.
Proof. See Theorem 3.3 of [6]. 2
Theorem 3.4. Suppose that all the stationary points of (12) are nondegen-
erate. ( @f
the hypothesis of Theorem
3.2 (and, hence, its thesis) must hold.
Proof. See Theorem 3.4 of [6]. 2
Theorem 3.5. Suppose that fx k g is a sequence generated by Algorithm 3.1
and let " be an arbitrary positive number. Then, there exists k 2 f0;
such that kg P
Proof. This result is a direct consequence of Theorems 3.2 and 3.3. 2
Implementation
At iteration k of Algorithm 2.1 the current iterate is x k and we are looking
for a direction d k satisfying condition (2). We use a truncated-Newton
approach to compute this direction. To solve the Newtonian system we call
Algorithm 4.1 (described below) with A  r 2 f(x k
The following algorithm applies to the problem
s:
The initial approximation to the solution of (14) is s The algorithm
nds a point s  which is a solution or satises q(s  ) < q(s 0 ). Perhaps, the -
nal point is on the boundary of the region dened by ksk   and  l  s   u.
Algorithm 4.1: Conjugate gradients
The parameters   << 1 and k max 2 IN are given. The algorithm starts with
Step 1. Test stopping criteria
set s
Step 2. Compute conjugate gradient direction
Step 2.1 If
else compute
Step 2.2 If (p T
Step 3. Compute step
Step 3.1 Compute
ug.
Step 3.2 Compute
Step 3.3 If (
If (
If (
Step 4. Compute new iterate
Step 4.1 Compute
Step 4.2 If (b T s k+1 > kbkks k+1
set s
Step 4.3 If (
set s
Step 5. Compute
set to Step 1.
This algorithm is a modication of the one presented in [27] (p. 529)
for symmetric positive denite matrices A and without constraints. The
modications are the following:
At Step 2.2 we test if p k is a descent direction at s k , i.e., if hp k ; rq(s k )i <
To force this condition we multiply p k by 1 if necessary. If the
matrix-vector products are computed exactly, this safeguard is not
necessary. However, in many cases the matrix-vector product Ap k is
replaced by a nite-dierence approximation. For this reason, we perform
the test in order to guarantee that the quadratic decreases along
the direction p k .
At Step 3.3 we test if p T
inequality holds, the step  k
in the direction p k is computed as the minimum among the conjugate-gradient
step and the maximum positive step preserving feasibility. If
and we are at the rst iteration of CG, we set  k  max .
In this way CG will stop with s
the angle condition of Step 4.2. If we are not at iteration zero of CG,
we keep the current approximation to the solution of (14) obtained so
far.
At Step 4.2 we test whether the angle condition (2) is satised by
the new iterate or not. If this condition is not fullled, we stop the
algorithm with the previous iterate. We also stop the algorithm if the
boundary of the feasible set is achieved (Step 4.3).
The convergence criterion  for the conjugate-gradient algorithm is dy-
namic. It varies linearly with the logarithm of the norm of the continuous
projected gradient, beginning with the value   i and nishing with   f . We
dene
where
a
log
log
and  is used in the stopping criterion kg P (x)k 1 <  of Algorithm 3.1.
The parameter k max is the maximum number of CG-iterations for each
call of the conjugate-gradient algorithm. It also varies dynamically in such
a way that more iterations are allowed at the end of the process than at the
beginning. The reason is that we want to invest a larger eort in solving
quadratic subproblems when we are close to the solution than when we are
far from it. In fact,
where
log
In the Incremental-quotient version of GENCAN, r 2 f(x k ) is not computed
and the matrix-vector products r 2 f(x k )y are approximated by
with In fact, only the components correspondent
to free variables are computed and the existence of xed variables
is conveniently exploited in (15).
5 Numerical experiments with the CUTE collec-
tion
In order to assess the reliability GENCAN, we tested this method against
some well-known alternative algorithms using all the non-quadratic bound-
constrained problems with more than 50 variables from the CUTE [12] col-
lection. The algorithms that we used for comparing GENCAN are BOX-
QUACAN [26] (see, also, [28]), LANCELOT [11, 12] and the Spectral Projected
Gradient method (SPG) (described as SPG2 in [7]; see also [8]). All
the methods used the convergence criterion kg P
stopping criteria were inhibited.
In GENCAN we used
Algorithm 3.1), and
Algorithm 4.1). In all algorithms we used  . The
parameters of (the line search of) Algorithm 3.2 were the default parameters
mentioned in [7] and the same used in (the line search of) Algorithm 2.1,
i.e.,
In LANCELOT we used exact second derivatives and we did not use
preconditioning in the conjugate-gradient method. The reason for this is
that, in GENCAN, the conjugate gradient method for computing directions
is also used without preconditioning. The other options of LANCELOT
were the default ones. A small number of modications were made in
BOX-QUACAN to provide a fair comparison. These modications were:
(i) the initial trust-region radius of GENCAN was adopted; (ii) the maximum
number of conjugate-gradient iterations was xed in
the accuracy for solving the quadratic subproblems was dynamic in BOX-
QUACAN varying from 0:1 to 10 5 , as done in GENCAN, (iv) the minimum
trust-region radius  min was xed in 10 3 to be equal to the corresponding
parameter in GENCAN.
The codes are written in Fortran 77. The tests were done using an ultra-SPARC
from SUN, with 4 processors at 167 MHz, 1280 mega bytes of main
memory, and Solaris 2.5.1 operating system. The compiler was WorkShop
Compilers 4.2 Oct 1996 FORTRAN 77 4.2. Finally we have used the
ag -O4 to optimize the code.
In the rst four tables we report the full performance of LANCELOT,
SPG, BOX-QUACAN, GENCAN (true Hessian) and GENCAN (Incremental-
quotients). The usual denition of iteration in LANCELOT involves only
one function evaluation. However, in order to unify the comparison we call
\iteration" to the whole process that computes a new iterate with lower functional
value, starting from the current one. Therefore, a single LANCELOT-
iteration involves one gradient evaluation but, perhaps, several functional
evaluations. At each iteration several trust-region problems are solved approximately
and each of them uses a number of CG-iterations. Problems
HADAMALS and SCON1LS have bounds where the lower limit is equal to
the upper limit. BOX-QUACAN does not run under this circumstances,
so the performance of this method in that situation is not reported in the
corresponding table. In these tables, we report, for each method:
IT: Number of iterations;
FE: Functional evaluations;
GE: Gradient evaluations;
CG: Conjugate gradient iterations, except in the case of SPG, where CG
iterations are not computed;
Time: CPU time in seconds;
nal functional value obtained;
of the projected gradient at the nal point.
The next 3 tables repeat the information of the rst ones in a more
compact and readable form. In Table 5 we report the nal functional value
obtained for each method, in the cases where there was at least one dierence
between them, when computed with four signicant digits.
In

Table

7 we report, for each method, the numbers FE and (GE+CG).
The idea is that a CG iteration is sometimes as costly as a gradient-evaluation.
The cost is certainly the same when we use the incremental-quotient version
of GENCAN. Roughly speaking, GE+CG represents the amount of work
used for solving subproblems and FE represents the work done on the true
problem trying to reduce the objective function.

Table

8 reports the computer times for the problems where at least one
of the methods used more than 1 second. The computer time used by
LANCELOT must be considered under the warning made in [9] page 136,
\LANCELOT [.] does not require an interface using the CUTE tools. It
is worth noting that LANCELOT exploits much more structure than that
Problem n IT FE GE CG Time f(x) kg P (x)k1
6.467D 06
QRTQUAD 120 144 178 145 570 1.39 3.625D+06 3.505D 06
CHEBYQAD 50 22 28 23 463 2.22 5.386D 7.229D 06
LINVERSE 1999 22 28 23 2049 47.22 6.810D+02 3.003D 06

Table

1: Performance of LANCELOT.
provided by the interface tools". As a consequence, although GENCAN used
less iterations, less functional evaluations, less gradient evaluations and less
conjugate-gradient iterations than LANCELOT in SCON1LS, its computer
time is greater than the one spent by LANCELOT. In some problems, like
QR3DLS and CHEBYQAD, the way in which LANCELOT takes advantage
of the SIF structure is also impressive.
Now we include an additional table that was motivated by the observation
of Table 7. It can be observed that the number of functional evaluations
per iteration is larger in GENCAN than in LANCELOT and BOX-
QUACAN. The possible reasons are three:
Many SPG-iterations with, perhaps, many functional evaluations per
iteration.
Many TN-iterations with backtracking.
Many TN-iterations with extrapolations.
We classify the iterations with extrapolation in successful and unsuccessful
ones. A successful extrapolation is an iteration where the extrapolation
produced a functional value smaller than the one corresponding to the rst
Problem n IT FE GE Time f(x) kg P (x)k1
7.896D 06
QRTQUAD 120 598 1025 599 0.20 3.624D+06 8.049D 06
HADAMALS 1024
CHEBYQAD 50 841 1340 842 33.75 5.386D 9.549D 06
NONSCOMP 10000 43 44 44 2.81 3.419D 10 7.191D 06

Table

2: Performance of SPG.
trial point. An unsuccessful extrapolation corresponds to a failure in the
rst attempt to \double" the steplength. Therefore, in an unsuccessful ex-
trapolation, an additional \unnecessary" functional evaluation is done and
the \next iterate" corresponds to the rst trial point. According to this, we
report, in Table 9, the following features of GENCAN (incremental-quotient
SPG-IT: SPG iterations, used for leaving the current face.
SPG-FE: functional evaluations in SPG-iterations.
TN-IT: TN iterations.
TN-FE: functional evaluations in TN-iterations.
TN-(Step 1)-IT: TN-iterations where the unitary step was accepted.
TN-(Step 1)-FE: functional evaluations in TN-iterations where the
unitary step was accepted. This is necessarily equal to TN-(Step 1)-
IT.
TN-(Backtracking)-IT: TN-iterations where backtracking was necessary

Problem n IT FE GE CG Time f(x) kg P (x)k1
5.742D 06
EXPQUAD 120
28 5.23 9.133D+03 6.388D 07
QRTQUAD 120 22 28 23 214 0.10 3.625D+06 5.706D 07
CHEBYQAD 50 52 66 53 960 45.70 5.387D 9.535D 06
6.559D 06

Table

3: Performance of BOX-QUACAN.
TN-(Backtracking)-FE: functional evaluations at iterations with back-tracking

TN-(Extrap(+))-IT: successful iterations with extrapolation.
TN-(Extrap(+))-FE: functional evaluations at successful iterations with
extrapolation.
TN-(Extrap( ))-IT: unsuccessful iterations with extrapolation.
TN-(Extrap( ))-FE: functional evaluations at unsuccessful iterations
with extrapolation. This number is necessarily equal to twice the
corresponding number of iterations.
Problem n IT FE GE CG Time f(x) kg P (x)k1
EXPQUAD 120
3.813D 06
NONSCOMP 10000 17 43 19
9.053D 06

Table

4: Performance of GENCAN (true-hessian version).
Problem n IT FE GE CG Time f(x) kg P (x)k1
EXPLIN 120 17 43 19
EXPQUAD 120 21 51 23 53 0.03 3.626D+06 2.236D 06
CHEBYQAD 50 31 43 2.929D 06
NONSCOMP 10000

Table

5: Performance of GENCAN (incremental-quotient version).
BDEXP 1.969D 2.744D 1.967D
QRTQUAD 3.625D+06 3.624D+06 3.625D+06 3.625D+06 3.625D+06
CHEBYQAD 5.386D 5.386D 5.387D 5.386D 5.386D
DECONVB 6.393D 09 4.826D 08 5.664D 6.043D
QR3DLS 2.245D 08 1.973D 05 1.450D
SCON1LS 5.981D 04 1.224D+00 | 1.269D 4.549D 04

Table

Final functional values.
Problem LANCELOT SPG BOX-QUACAN GENCAN-QUOT
FE GE+CG FE GE FE GE+CG FE GE+CG
43 58
43
EXPQUAD
QRTQUAD 178 715 1025 599 28 236 75 101
CHEBYQAD 28 486 1340 842 66 1012 43 918
28 2072 1853 1023 19 415 34 87
SCON1LS 9340 5750468 7673022 5000002 | | 8565 4995260

Table

7: Functional and equivalent-gradient (GE+CG) evaluations.
MCCORMCK 4.24 2.27 5.23 4.57 3.56
HADAMALS 4.40 1.63 | 1.80 1.19
CHEBYQAD 2.22 33.75 45.70 13.86 22.26
QR3DLS 439.31 2203.97 2286.09 976.13 523.50

Table

8: Computer time.
Type of GENCAN iterations Details of Truncated Newton iterations
SPG Iteration TN iterations Step=1 Backtracking Extrap(+) Extrap( )
Problem IT FE IT FE IT FE IT FE IT FE IT FE

Table

9: GENCAN features.
Observing Table 9 we realise that:
1. The number of SPG-iterations is surprisingly small. Therefore, only
in few iterations the mechanism to \leave the face" is activated. So,
in most iterations, the number of active constraints remains the same
or is increased. Clearly, SPG-iterations are not responsible for the
relatively high number of functional evaluations.
2. The number of iterations where backtracking was necessary is, also,
surprisingly small. Therefore, extrapolations are responsible for the
functional-evaluations phenomenon. Since an unsuccessful extrapolation
uses only one additional (unnecessary) functional evaluations, its
contribution to increasing FE is also moderate. In fact, unsuccessful
extrapolations are responsible for 116 functional evaluations considering
all the problems, this means less than 8 evaluations per problem.
It turns out that many functional evaluations are used in successful
extrapolations. Considering the overall performance of the method
this seems to be a really good feature. An extreme case is BDEXP,
where only one TN-iteration was performed, giving a successful extrapolation
that used 11 functional evaluations and gave the solution
of the problem.
Further remarks
Convergence was obtained for all the problems with all the methods
tested, with the exception of SPG that did not solve SCON1LS after more
than thirty hours of computer time. The method that, in most cases, obtained
the lowest functional values was GENCAN-QUOT, but the dier-
ences do not seem to be large enough to reveal a clear tendency.
As was already mentioned in [7], the behavior of SPG is surprisingly
good. Although it is the only method that fails to solve a problem in reasonable
time, its behavior in the problems where it works is quite e-cient.
This indicates the existence of families of problems where SPG is, probably,
the best possible alternative. This observation has already been made in [8].
BOX-QUACAN has been the less successful method in this set of ex-
periments. This is not surprising, since the authors of [16] had observed
that this method outperformed LANCELOT in quadratic problems but is
not so good when the function is far from being quadratic. In fact, it was
this observation what motivated the present work. Nevertheless, there is
still a large scope for improvements of BOX-QUACAN, as far as we take
into account that improvements in the solution of quadratic subproblems are
possible and that sophisticated strategies for updating trust-region radius
can be incorporated.
6 Experiments with very large problems
We wish to place q circles of radius r in the rectangle [0; d 1 in such
a way that for all qg, the intersection
between the circle i and the circle j is at most one point. Therefore, given
I i  the goal is to determine c
solving the problem:
Minimize
subject to
r  c i
r  c i
The points c are the centers of the desired circles. If the objective
function value at the solution of this minimization problem is zero, then the
original problem is solved.
When I the problem above is
known as the Cylinder Packing problem [22]. The present generalization is
directed to Sociometry applications.

Table

describes the main features of some medium and large scale
problems of this type. In the problems 9-15 the sets I i were randomly
generated with the Schrage's random number generator [35] and seed = 1.
In all cases we used 0:5. Observe that n, the number of variables, is
equal to 2 q.

Tables

11 and 12 show the performances of GENCAN and LANCELOT.
Internal limitations of the big-problems installation of CUTE forbid the
solution of larger instances of this problems using SIF. We show the CPU
Problem n #I i box

Table

10: Medium- and large-scale classical and modied cylinder packing
problems.
times of GENCAN both using SIF (SIF-Time) and Fortran subroutines (FS-
for computing function and gradient. We used a random initial point
(generated inside the box with the Schrage's algorithm and seed equal to 1).
Both methods found a global solution in all the cases. In Table 12 we also
report the number of free variables at the solution so far found by GENCAN.
GENCAN
using Fortran subroutines using SIF LANCELOT
IT FE GE CG Time IT FE GE CG Time IT FE GE CG Time

Table

11: GENCAN and LANCELOT with cylinder packing problems.
9
28 10649.79 9914682

Table

12: GENCAN with very large problems.
7 Final remarks
Numerical algorithms must be analyzed not only from the point of view of its
present state but also from considerations related to their possibility of im-
provement. The chances of improvement of active-set methods like the one
presented in this paper come from the development of new unconstrained algorithms
and from the adaptation of known unconstrained algorithms to the
specic characteristics of our problem. In our algorithm, the computation
of the search direction is open to many possibilities. As we mentioned in
the introduction, a secant multipoint scheme (with a dierent procedure for
leaving the faces) was considered in [10] and a negative-curvature Newtonian
direction for small problems was used in [6], where leaving faces is also associated
to SPG. A particularly interesting alternative is the preconditioned
spectral projected gradient method introduced in [30].
The extension of the technique of this paper to general linearly constrained
optimization is another interesting subject of possible research.
From the theoretical point of view, the extension is straightforward, and the
convergence proofs do not oer technical di-culties. The only real di-culty
is that we need to project onto the feasible set, both in the extrapolation
steps and in the SPG iterations. In theory, extrapolation can be avoided
without aecting global convergence, but projections are essential in SPG
iterations. Sometimes, the feasible polytope is such that projections are
easy to compute. See [8]. In those cases, the extension of GENCAN would
probably be quite e-cient.

Acknowledgements

.
The authors are very grateful to Nick Gould, who helped them in the
use of the SIF language. We are also indebted to two anonymous referees
whose comments helped us a lot to improve the nal version of the paper.



--R

On the resolution of the generalized nonlinear complementarity problem.

A limited memory algorithm for bound constrained minimization.
Restricted opti- mization: a clue to a fast and accurate implementation of the Common Re ection Surface stack method




CUTE: constrained and unconstrained testing environment.

Global convergence of a class of trust region algorithms for optimization with simple bounds.
A globally convergent augmented Lagrangean algorithm for optimization with general constraints and simple bounds

Numerical methods for unconstrained optimization and nonlinear equations.

Comparing the numerical performance of two trust-region algorithms for large-scale bound-constrained minimization





Optimising the palletisation of cylinders in cases




Matrix Computations.


Preconditioned spectral gradient method for unconstrained optimization problems


On the Barzilai and Borwein choice of steplength for the gradient method
The Barzilai and Borwein gradient method for the large scale unconstrained minimization problem.
A more portable Fortran random number generator.
A class of inde
--TR
Global convergence of a class of trust region algorithms for optimization with simple bounds
A globally convergent augmented Lagrangian algorithm for optimization with general constraints and simple bounds
A limited memory algorithm for bound constrained optimization
Matrix computations (3rd ed.)
Gradient Method with Retards and Generalizations
Estimation of the optical constants and the thickness of thin films using unconstrained optimization
A More Portable Fortran Random Number Generator
Trust-region methods
Validation of an Augmented Lagrangian Algorithm with a Gauss-Newton Hessian Approximation Using a Set of Hard-Spheres Problems
Duality-based domain decomposition with natural coarse-space for variational inequalities0
Algorithm 813
On the Resolution of the Generalized Nonlinear Complementarity Problem
A Class of Indefinite Dogleg Path Methods for Unconstrained Minimization
The Barzilai and Borwein Gradient Method for the Large Scale Unconstrained Minimization Problem
Nonmonotone Spectral Projected Gradient Methods on Convex Sets
Newton''s Method for Large Bound-Constrained Optimization Problems
Constrained Quadratic Programming with Proportioning and Projections
Augmented Lagrangians with Adaptive Precision Control for Quadratic Programming with Equality Constraints

--CTR
G. Birgin , J. M. Martnez, Structured minimal-memory inexact quasi-Newton method and secant preconditioners for augmented Lagrangian optimization, Computational Optimization and Applications, v.39 n.1, p.1-16, January   2008
G. Birgin , R. A. Castillo , J. M. Martnez, Numerical Comparison of Augmented Lagrangian Algorithms for Nonconvex Problems, Computational Optimization and Applications, v.31 n.1, p.31-55, May       2005
G. Birgin , J. M. Martnez , F. H. Nishihara , D. P. Ronconi, Orthogonal packing of rectangular items within arbitrary convex regions by nonlinear optimization, Computers and Operations Research, v.33 n.12, p.3535-3548, December 2006
active-set strategies;numerical methods;box-constrained minimization;Spectral Projected Gradient