A Practical Algorithm for General Large Scale Nonlinear Optimization Problems.
We provide an effective and efficient implementation of a sequential quadratic programming (SQP) algorithm for the general large scale nonlinear programming problem. In this algorithm the quadratic programming subproblems are solved by an interior point method that can be prematurely halted by a trust region constraint.  Numerous computational enhancements to improve the numerical performance are presented.  These include a dynamic procedure for adjusting the merit function parameter and procedures for adjusting the trust region radius.  Numerical results and comparisons are presented.
Introduction
. In a series of recent papers, [3], [6], and [8], the authors have
developed a new algorithmic approach for solving large, nonlinear, constrained optimization
problems. This proposed procedure is, in essence, a sequential quadratic
programming (SQP) method that uses an interior point algorithm for solving the
quadratic subproblems and achieves global convergence through the application of a
special merit function and a trust region strategy. Over the past several years the
theory supporting this approach has been analyzed and strengthened. This theory is
presented in a companion paper [4]. In addition, implementations of the algorithm
have been extensively tested on a variety of large problems, including standard test
problems and problems of engineering and scientific origin, ranging in size from several
hundred to several thousand variables with up to several thousand constraints.
Specific strategies have been developed for handling the parameters utilized by the
algorithm and for dealing with nontrivial pathologies (e. g. , linearly dependent active
constraint gradients or inconsistent linearized constraints in the quadratic subprob-
lem) that often occur in large scale problems. In this paper we present the results of
these efforts.
Based on its theoretical foundation and on our numerical experience we are confident
that this algorithm provides an efficient means for attacking a large, sparse,
nonlinear program with equality and/or inequality constraints. Rigorous comparisons
of algorithms for large nonlinear problems is notoriously difficult, especially given the
extensive set of options typically available in codes for such problems. Nevertheless,
our algorithm, with the (conservative) default parameter settings, has been successful
on problems that have caused difficulties for other algorithms and, consequently, we
are encouraged to believe that it is competitive at the current stage in the development
of methods for solving these large problems.
Below we give an outline of our basic procedure and in the succeeding sections we
provide more specific detail on the component parts of the implemented algorithm,
including the strategies and safeguards that we have used. We also exhibit and comment
on the results of some of our numerical tests. This paper relies heavily on the
results from the paper on the theory for motivation of the basic ideas.
Applied and Computational Mathematics Division, National Institute of Standards and Tech-
nology, Gaithersburg, MD 20899
y Department of Mathematical Sciences, Carnegie Mellon University, Pittsburgh, PA 15213-3890
z Mathematics Department, University of North Carolina, Chapel Hill, NC 27599
We assume the general nonlinear programming problem to be of the form
subject to: g(x)  0
are smooth functions. Nonlinear equality
constraints are not included in our description here in order to avoid distracting
technicalities. The modifications necessary for their insertion can be inferred from [6].
Nonlinear equality constraints are included in our code and in some of the problems
we tested. The sequential quadratic programming method is the backbone of our
algorithm. (See [7] for a review of these techniques.) At the kth step we have an
iterate, x k , denoting the current approximation to the solution of (NLP ). In addition
to the x-iterate we also maintain a non-negative iterate, z k 2 R m , which measures
the infeasibility at x k . At this stage (NLP ) is modeled by a quadratic program of
the form
min
subject to: rg(x k
is taken to be an appropriate approximation to the Hessian of the Lagrangian
for (NLP ), i.e.,
where
and H xx represents the Hessian with respect x of the function to which it is applied.
(See Section 4.5 for a discussion of the choice of B k used in our numerical experiments.)
In this form (QP ) generates a step that provides a search direction for improving the
current iterate.
There are two significant points to be made concerning this phase of our algorithm.
First, we apply an interior point quadratic program solver to (QP ); more specifically,
we use the method found in [1] where solutions are calculated by solving a sequence of
low dimensional quadratic programs. Pertinent details of this solver and its properties
relative to its use in our SQP method can be found in Section 2. Second, we do not
try to solve (QP ) with complete accuracy at each iteration; rather, we often terminate
the interior point method prematurely. In particular, we halt the quadratic program
solver when the steplength exceeds a "trust region radius" that is modified at each
iteration according to how well the improvement in our merit function is predicted.
Thus our algorithm can be said to be a "truncated Newton method" in the sense
of [18] (see also [15]). This particular merit function and a more useful "working
version" are discussed in Section 3 and our strategy for updating the trust region
radius is given in Section 4.2.
The output of the (QP ) solver is a vector that determines the direction of the
step in the x-variable, which in turn yields a step direction for the "slack" variable
z as explained in Section 3. The combined step direction of these two variables is a
descent direction for the working version of the merit function and also for constraint
thus we can choose the steplength in this direction to decrease the merit
function and/or the infeasibility of the iterate. The choice of steplength determines
the new iterate x k+1 and also the new value z k+1 . The strategy for choosing the
steplength and other algorithmic details, including the modifications and safeguards
necessary to make an implementation robust, are given in Section 4.
The results of our numerical tests are contained in Section 5. These results
demonstrate the overall effectiveness of the procedure and highlight the beneficial
effect of our trust region strategy and other procedures. Finally, in Section 6 we
briefly consider weaknesses in the current version of the algorithm and suggest possible
avenues of research to improve its efficiency.
For a discussion of the theoretical and practical questions related to large scale
nonlinear programming see the recent surveys [12], [14] and [21].
2. An Interior Point QP Solver. Interior point methods for linear programming
have been demonstrated to be very successful, especially on large problems,
and recent research has lead to their extension to quadratic programs. A particular
method, the method of optimizing over low-dimensional subspaces, has performed well
on linear programs and has been extended to the quadratic programming case (see
[1], and [2] and the references contained therein). This method, for which good numerical
results for quadratic programs have been reported, has properties that make
it particularly compatible with the SQP algorithm we are describing in this paper. A
brief description of the essential features of this method and their importance for our
purposes follow. The many details of the actual algorithm that are not reported here
may be found in the above references.
The quadratic program that we solve, (QP ), has the form
subject to: A T s
n\Thetan , A 2 R n\Thetam , and b 2 R m . The assumptions on (2.1) that
are necessary to apply the interior point algorithm are that the problem be bounded,
that A have full column rank, and that there exist feasible points (i.e., that the
constraints be consistent). Note that Q can be indefinite and that no assumption of
a full-dimensional interior is required. If equality constraints are present, they are
handled by writing them as two inequalities.
An important prerequisite for solving (2.1) by an interior point method is a feasible
initial point. Our algorithm uses a "Big M " method to construct the Phase I
problem
min
subject to: A T s
where e is a vector of ones and ' is the "artificial" variable. Clearly for '   large
enough the point (s; is feasible for (2.2) and if M is sufficiently large the
algorithm applied to (2.2) will reduce ' until the artificial variable is nonpositive, at
which point the current value of s is feasible and the M' and e' terms are dropped.
If no such value of the artificial variable can be found, then (2.2) is not consistent and
the algorithm stops. As discussed below, we make use of the step obtained from (2.2)
even if it is not feasible for (QP ). Note that when equality constraints are present,
the entire solution procedure takes place in Phase I and ' will always be present.
The defining characteristic of the algorithm is that it proceeds by solving a sequence
of low-dimensional subspace approximations to (2.1). In our application we
follow the reported results in which the dimension of the subspace is taken as three.
The following is an outline of the O3D (for Optimizing over 3-Dimensional subspaces)
version of the algorithm. As the variable ' is treated essentially the same as the
components of s in the O3D algorithm (see, however, Step 6 below) the dependence
on ' is incorporated into the formulation given in (2.1).
O3D Algorithm for Quadratic Programming
1. Given a feasible point, s 0
2. Generate 3 independent search directions be the
matrix whose columns are p i .
3. Form and solve the restricted quadratic program
subject to: A T ~ s
where ~
Call the solution i   .
4. Set s j+1 := s for an appropriate value of the steplength ae 2 (0; 1):
5. If stopping criteria are met, exit.
6. Go to 2. (At this step, if the component of the vector s corresponding to
the artificial variable ' has become nonpositive, it is eliminated from the
problem.)
The search directions in Step 2 are solutions to
\Theta AD 2 A T
where fi is a scalar depending on the current iterate,
with , and the t i are particular values chosen such that one of these
directions is always a descent direction with respect to the objective function. The
steplength ae is set to the lesser of 99% of the distance to the boundary and the distance
to the minimum of the objective function.
The form of the matrix in (2.3) allows for efficient exploitation of the sparsity.
Note that if Q is positive semi-definite, then the matrix in (2.3) is positive definite for
all interior points; otherwise, it may not be. In the latter case, a modification similar
to that in [20] is used. In our application of this algorithm, using this procedure
obviates the need for the matrix B k to be positive definite, which in turn allows us
to use the Hessian of the Lagrangian or a finite difference approximation thereof.
The standard stopping criterion for the algorithm is that at least one of the following
holds: (a) the relative change in two successive values of the objective function
is small; (b) the relative difference between the primal and the dual objective function
values is small; or (c) the difference between two successive iterates is small. For use
in our SQP algorithm we have added: (d) the length of the solution vector exceeds
a specified value. This additional condition has been implemented to allow for trust
region strategies; in particular, this criterion will cause the algorithm to halt if (QP )
is unbounded. In any case, the terminal vector will be a useful direction in the context
of our purposes; this point will be discussed in the next section.
The most recent version of O3D described in [1] contains an option to perform
a special "recentering step" after each subspace optimizing step (Step 4) that has
generally improved the efficiency. This option is not used in the results reported here.
(See Section 6 for a further comment.)
3. Updating the Iterates: the Merit Functions. In this section we review
the definitions and properties of our merit functions and provide formulas for updating
the iterates. The reader is referred to the companion paper for proofs and motivations
of these concepts.
As stated in Section 1, at each iteration our algorithm yields a pair
x k is an approximation to the solution of (NLP ) and z k is the corresponding approximate
slack vector. The step directions for the updated values of these approximations
are based on the (approximate) solution, to the quadratic program
min ffi;'
subject to: rg(x k
obtained as described in the preceding section. The vector ffi k gives the step direction
for x k and we determine the step direction, q k , for the slack vector z k by the formula
\Theta rg(x k
Note that if ffi k is feasible for (QP ) then
\Theta rg(x k
In this case z k is the slack vector for (QP ) corresponding to ffi k and thus is the
slack variable for the linear approximation of g(x k+1 ). Given the step direction we
then update the iterate by means of the formulas
z
for some value of the steplength parameter ff: Observe that if z k  0 then the fact
that means that z k+1 will be non-negative if ff 2 [0; 1]. In
our algorithm the non-negativity of the slack vector iterates is preserved and, in fact,
it sometimes turns out to be useful to maintain the z k at a positive level (see Section
4.8).
It is important to emphasize that the ffi k are determined by (QP ), the quadratic
approximation to (NLP ), and are not dependent on the choice of z k . The z k are
generated solely for use with the merit function described below. That is, we do not
solve the slack variable problem. A comment on the notation is also in order at this
point: We denote the iterate by and the step by (ffi k ; q k ), whereas conventional
notation would be to use ' x k
z k
and
It should be clear from the context what is meant.
In optimization algorithms the value of a steplength parameter is generally chosen
so as to reduce the value of a suitably chosen merit function. Typically, a merit function
for (NLP ) is a scalar-valued function that has an unconstrained minimum at x   ,
a solution to (NLP ). Because a reduction in this function implies that progress is being
made towards the solution, it can be used to determine an appropriate steplength
in a given search direction.
In [5] and [6] a merit function for equality-constrained problems was derived that
has important properties vis-a-vis the steps generated by the SQP algorithm. Using
a slack-variable formulation of (NLP ) a merit function for the inequality constrained
problem can be constructed having the form
d  c(x; z) T
where z is nonnegative, d is a scalar,
and
We use this merit function (and its approximations defined below) for choosing the
value of the steplength parameter ff. As noted above, the approximate slack vectors
generated by our algorithm, z k , always remain non-negative; thus the non-negativity
constraint on the z for / d imposes no theoretical difficulty.
The function  c(x; z) defined above plays an important role in our algorithm as it
is used to measure the feasibility of the pair (x; z). That is, if we define the function
where k\Deltak denotes the standard Euclidean norm and set
then C 0 corresponds to the feasible set of (NLP ) and hence close to feasible
if it is in C j for small j.
For d sufficiently small the merit function / d has the desirable property that a
solution of (NLP ) corresponds to a (constrained) minimum of / d . In addition, if d is
small and ffi k is the exact solution to (QP ) (which implies that ' then the step
descent direction for / d when sufficiently close to feasibility.
Despite these useful properties, / d has two deficiencies that limit its use in an efficient
algorithm. First, (ffi k ; q k ) is a descent direction of / d only near feasibility, and, second,
the evaluation of rf and rg and additional nontrivial computational algebra are
required to assess a prospective point. In order to overcome these difficulties, the
approximate merit function
d
d
where
is developed as a "working" version of / d at As the values of
k and
are fixed, / k
d can be more easily evaluated than / d in a line search algorithm for
choosing an appropriate value of ff. This approximate merit function, / k
d , not only
has essentially the same properties as / d with respect to the step
has the stronger property that the step is a descent direction for / k
d everywhere.
Moreover, for sufficiently small and outside of a ball around the solution
a "sufficient" reduction in / k
d implies a "sufficient" reduction in / d . (We mean by
"sufficient" reduction that a Wolfe condition is satisfied.) Thus we are able to use / k
d
as a surrogate for / d for testing the progress of our iterates towards a minimum.
A further important property of the step ffi k , under the assumption that it is the
exact solution to (QP ), is that it is a descent direction for the function r defined by
(3.4). Thus a basic algorithm for the case where the (QP ) can be solved exactly is as
follows: Given an initial value of j use the steps (ffi k ; q k ) to reduce r until the iterates
are in C j . Once the iterates are contained in C j if a sufficient reduction in / k
d does
not yield a sufficient reduction in / d then reduce j. If, in the course of the algorithm,
remains bounded away from zero, then convergence follows from the fact that the
Wolfe condition is satisfied for / d . If j goes to zero, then convergence follows from
the observation that the radius of the ball in which the Wolfe condition is not satisfied
also goes to zero. This is essentially the algorithm for which global convergence is
proved in the paper on the theory.
In this paper we are primarily interested in enhancements that convert the theoretical
algorithm into one that is practical and efficient. This requires that we make
provisions for situations when the assumptions under which we performed the convergence
analysis are not valid and that we adopt numerical procedures to reduce the
computational effort. As we note below, not all of these modifications have been (or
even can be) theoretically justified, but we believe that the firm foundation of the
underlying algorithm and the evidence accumulated in extensive numerical testing
validate their use.
In the implementation of our algorithm a trust region constraint is used that
possibly truncates the quadratic programming algorithm before an exact solution
is achieved. In this case the theory described above does not apply for the step
obtained from the approximate solution, Although a general
convergence theory based on this step is not yet available, it is shown in the theory
paper that if the approximate solution is obtained from the O3D algorithm and if ' k
is not too large then the resulting step has the appropriate descent properties for the
functions r, / d , and / k
d at In particular, convergence can be achieved if ' k
goes to zero in a suitable manner. These properties justify our use of the truncation
procedure to speed up the algorithm. It is important to note that this approximation
procedure also allows us to handle the difficulty that arises in sequential quadratic
programming methods when the quadratic subproblem is inconsistent.
4. The Truncated SQP Algorithm. In this section we give a somewhat detailed
description of our algorithm. Initially we assume that the Hessian approxima-
are positive definite, the matrices
A k , are nonsingular, and the linearized
constraints in (QP ) are consistent. In real-world applications these assumptions are
not always valid so we have tried to make our algorithm flexible enough to perform
well in situations where these assumptions fail to hold. We describe some of these
adaptations at the end of this section.
The implementation of the algorithm depends upon four important parameters
that need to be either computed or modified throughout the course of the algorithm.
The globalization parameter, j, was introduced in (3.5). It is a measure of the size
of the domain about the feasible region in which the direction (ffi k ; q k ) is a descent
direction for the true merit function / d . A current estimate of j is maintained in the
algorithm. The trust region parameter,  , is an upper bound on the (weighted) norm
of our approximate solution to (QP ),
where D is a positive definite diagonal matrix. The trust region radius  is updated
at every iteration. The parameter, ff, is the steplength parameter. It determines
the length of the step in the variables (x; z) in the direction (ffi k ; q k ). It is chosen
to guarantee progress towards the solution in decreasing either the merit function
or infeasibility. Finally, d, the merit function parameter, must be small enough to
guarantee that the theoretical properties described in the preceding section are valid.
Although the theory allows arbitrarily small values of d, the algorithm becomes very
slow if d is too small, thus it is monitored throughout the algorithm and either increased
or decreased as appropriate.
The outline of the algorithm is followed by specific comments on the procedures
and their justifications. This version contains some of the practical modifications
described above. To simplify the notation we define
Recall that r is given by (3.4).
Basic Truncated SQP Algorithm
1. Initialization: Given x
a. Initialize the slack variable z 0  0;
b. Set k := 0:
2. Calculation of the basic trust region step:
a. While kffik !  , iterate (using O3D) on
min
subject to: rg(x k
to obtain ffi k and ' k .
b. Set
\Theta rg(x k
\Theta rg(x k
otherwise
c. Decrease d if necessary.
3. Computation of the steplength parameter:
a. Choose ff 2 (0; 1] such that / k
d is sufficiently reduced.
b. If
reduce ff if necessary until r is sufficiently reduced.
c. If reduce ff if necessary so that
4. Update of the estimate of the globalization parameter:
a. If
set
5. Update of the variables and check for termination:
a. Set
z
b. If convergence criteria are met, quit.
c. Update B k to B k+1 .
6. Adjustment of the merit function and trust region parameters:
a. Update d if necessary.
b. Adjust the trust region radius  .
7. Return:
a.
b. Go to Step 2.
4.1. The Globalization Parameter. The globalization step is based on work
in [6] and [4]. In Step 3 we require that the approximate merit function be reduced
and, in addition, if the current iterate lies outside the set C j we require that the
constraint infeasibilities also be reduced. This is possible as a result of the descent
properties described in Section 3. If we have a good estimate of j and
then the true merit function can also be reduced; if this is not the case, then our
estimate of j is too large and we reduce its value in Step 4. This procedure will
eventually lead to a sufficiently small value of j. Note that this arrangement allows
steps that may increase the merit function, but only in a controlled way. It also allows
steps that may increase the constraint infeasibilities, but only when inside of C j .
4.2. Updating  . Our procedure for updating  , the trust region radius, in
Step 6b is similar to the standard strategy used in trust region algorithms (see [17] or
[31]) in that we base the decision on how to change  on a comparison of a predicted
relative reduction, pred k ; and an actual relative reduction, ared k , in a function used
to measure the progress toward the solution. (Various formulas for the predicted relative
reduction, pred k , have been suggested for different merit functions, especially for
equality constrained programming problems; see, for example, [19]). What is distinctive
about our procedure is that we use different functions for computing pred k and
ared k depending on the current status of the algorithm. When the linearized constraints
are satisfied we use the approximate merit function to compute the predicted
and actual reductions. When the trust region constraint causes O3D to terminate in
Phase I, i.e., when the linearized constraints are not satisfied, predicted and actual
reductions in infeasibility are used.
In the case when a feasible solution to (QP ) is obtained then / k
d is used to
compute the predicted and actual reductions. Our method for defining pred k differs
from the standard methods used in unconstrained optimization because the step-
finding subproblem is not based solely on the merit function and, moreover, the
trust region constraint does not appear explicitly in the subproblem. Nevertheless
in updating  we want to assess how well an approximation to / k
d agrees with / k
d in
the direction uses a quadratic approximation of the Lagrangian
for the objective function with linearized constraints, we form our approximation to
d based on a quadratic approximation to the function / k
1 given by
and a linear approximation to
Note that / k
d
z). Based on these considerations and the
results of [16] we define the predicted relative reduction by
pred
ae
d
oe
where the derivatives are with respect to x and z and the steplength parameter ff k is
the size of the most recently accepted step. The value of the actual relative reduction,
ared k , is taken to be the difference in the values of / k
d at the points
and divided by the value of / k
valid criticism of the formula for
pred k is its dependence on higher order derivatives. Therefore we use the available
approximation of the Hessian of the Lagrangian for r 2 / k
1 . For example, cell-centered
finite difference approximations to the Hessian of the Lagrangian function were used
in the numerical results presented here, unless analytic second derivative formulas
were readily available.
The above choice for pred k is not used when the step returned by O3D is not
feasible. In these situations the resulting step is dominated by a feasibility improving
component and it makes little sense for the adjustment to  to be determined by / k
rather, a comparison of the predicted and actual improvement in constraint infeasibility
seems more appropriate. Therefore, in this case the function r(x; z) is used for
comparison purposes. The values of pred k and ared k are given as follows for the case
when the O3D algorithm terminates in Phase I:
pred
and
ared
These heuristics for choosing pred k and ared k appear to work well. Specifically,
they allow the trust region radius,  , to be increased even in the event that the step
returned by O3D does not satisfy linearized constraints or it results in an increase in
the true merit function. In our experience, the alternative formulas based solely on
constraint violations never are employed close to the solution. Indeed, the iterates
preceding convergence have always been observed to be well inside C j where satisfying
the linearized constraints and decreasing the merit functions usually pose no problem.
4.3. The steplength ff. The steplength ff is determined in Step 3 of the al-
gorithm. The "sufficient decrease" referred to in 3a and 3b requires that the Wolfe
condition be satisfied. For a given function OE and potential step w from point v this
condition requires that ff satisfy
for some fixed oe 2 (0; 1). In the numerical experiments reported in Section 5 we employed
a simple backtracking procedure (with factor one-half) to find ff to satisfy this
condition for both / k
d and for r. We have also experimented with more sophisticated
line search methods motivated by unconstrained optimization techniques as in [18],
but the observations to date suggest that the more complicated line searches result in
very little improvement of our algorithm, except when the iterates are quite far from
the solution.
4.4. Adjusting d. Choosing an effective value for the merit function parameter
d is essential in our algorithm. While it is clear that (in a compact set) a sufficiently
small value of d will assure that the results given in [4] are valid, there are three
very important practical reasons why the parameter must be adjusted rather than
fixed. First, if the angle between the direction generated by O3D and the gradient of
the approximate merit function becomes nearly orthogonal the steps might become
too small. We adjust d to avoid this possibility. Second, the approximate merit
d , is changing at each iteration and it is possible a previous iterate might
be acceptable to the current / k
d , i.e., cycling might occur. This worry can also be
alleviated by adjusting d. A third reason for changing d is to allow for larger steps.
It is seen from the theory and has been verified by numerical experience that if d
is too small then the form of the merit function forces the path of the iterates to
follow the "nearly active" constraints closely. This causes the algorithm to take very
small steps and, in particular, to be slow in moving away from a nonoptimal active
set. By making it possible to increase d we can significantly improve the algorithm's
performance.
In the implementation of our algorithm there are two opportunities to adjust d:
in Step 2, after solving the quadratic subproblem, and in Step 6, after the step has
been taken. In the first of these adjustments d can only be decreased; in the second,
the parameter may be increased or decreased.
In Step 2, the angle between the gradient of the approximate merit function
d and the step direction (ffi k ; q k ) is computed. If these two vectors become nearly
orthogonal, we conclude that d is not small enough to ensure a good decrease in / k
d ,
and we decrease the parameter. To be more specific, we compute
If w(d)  \Gamma:1 we calculate a value "
d so that w( "
d)  \Gamma:5: We safeguard the procedure
by not allowing more than a certain percentage decrease in d. In the current version
we use 50%.
If d was not decreased in Step 2 we consider modifying it after a step has been
taken (Step 6). Here the primary concern is to avoid cycling. To do so we compute
an interval for the penalty parameter as follows. For a fixed integer  we seek a value
of the parameter,
d, such that
d
Inequality (4.2) implies that none of the past  iterates will be acceptable to the
approximate merit function with the new value of
d. (Thus if cycling would
be possible). To accomplish this, we use the decomposition
1 and / k
are defined in Section 4.2. We then compute the values of / k
and / k
consider the inequalities
d
d
We define d u
i and d l
i to be the upper and lower values of d that ensure that inequality
(4.4) is satisfied. Then letting
d
and
we obtain an interval (d l ; d u ). Assuming that this interval exists it is the case that
if the value of d for the next step is chosen in this interval, the next iterate will
not return to one of the previous  iterates. In practice a value of   5 is usually
more than sufficient to prevent cycling. If the interval doesn't exist, then we make no
change.
Given that we can choose d to avoid cycling, our second objective at this juncture
is to increase d to allow bigger steps. If the d u is larger than the current d then we
can safely increase d without worrying about possible cycling. However, we safeguard
this increase in two ways. First, we require that the predicted reduction based on the
approximate merit function must be greater than the predicted reduction of infeasibility
in the linearized constraints. This restriction prevents d from being increased
prematurely due primarily to a large decrease in constraint infeasibilities. Specifically,
writing the predicted reduction in / k
d (see (4.1)) as
we insist that for a new value of d
Second, we use a maximum allowable change (currently a factor of 2) to limit the
growth of d. Computationally, these simple procedures for updating d appear to be
effective, especially in the presence of highly nonlinear constraints and poorly scaled
problems.
4.5. The Hessian Approximation. In the numerical experimentation reported
here, we have used a finite difference approximation to the Hessian of the Lagrangian
. Although the Hessian of the Lagrangian at a strong solution is positive definite
on the appropriate subspace, it may be indefinite in general. Even if it is positive
definite the finite difference approximation may not be. We experimented with two
approaches for handling this possibility. First, we simply modified the approximate
Hessian matrix by adding non-negative elements to the diagonal ensuring that the
Cholesky factorization of the matrix had positive elements along its diagonal (see
[20]). This modification was easy to implement, but it was observed to slow convergence
on some problems. While this modification guarantees that a positive definite
matrix will be delivered to the (QP ) solver, if it takes place when the iterates get
close to the solution, it generally precludes local q-superlinear convergence.
An alternative to modifying the approximate Hessian of the Lagrangian is simply
to allow O3D to iterate on the indefinite QP subproblem, halting the iterations when
the solution exceeds the trust region radius. We implemented this approach and
it seemed to yield superior results to those obtained by making the approximate
Hessian positive definite (especially when the iterates were close to a solution) even
though, theoretically, we can only prove that we obtain a descent direction when the
approximate Hessian is positive definite.
4.6. Convergence Criteria. The convergence criteria used are standard, and
similar to those in [3]. We first insist that the constraints be satisfied to a close
tolerance; specifically we require
We also require that either
or
The criterion (4.9) is a stronger indication that a KKT point has been reached. The
weaker criterion (4.10) suggests that progress slowed drastically and that iterates may
or may not have drawn close to a solution. For this reason criterion (4.9) is usually
preferable to criterion (4.10). The Lagrange multipliers returned by the quadratic
program are used in (4.9) unless the trust region constraint determines the approximate
solution of the (QP ). In that case, we use the least squares approximation to
the multipliers, replacing all negative multipliers with machine zeros. In all of the
problems solved to date, the trust region never comes into play when the iterates get
close to the solution; therefore the (QP ) multipliers are used for the convergence test
at the solution.
4.7. Inconsistent Quadratic Subproblems. One difficulty that can occur
when making linear approximations to nonlinear constraints is that (QP ) may be
inconsistent. In this case O3D will, even if it runs to completion, not exit Phase
I and will return a positive value of the artificial variable. (Note that this always
occurs if equality constraints are present.) For small ' the resulting direction is a
descent direction for / k
d and for r. As a result, the step taken in this direction
will generally decrease infeasibility, making it less likely that an inconsistent set of
linearized constraints will be encountered during subsequent iterations.
More recent versions of our algorithm include a constraint relaxation procedure
that appears to yield an acceptable step, even in the event that inconsistent
linearizations of constraints are encountered. Because this situation did not surface
during the numerical experiments presented in this paper, we do not include a description
of our perturbation procedure. We do note, however, that we have encountered
important application problems where this procedure was crucial to the performance
of our algorithm (see for example [24]).
4.8. Updating slack variables. One difficulty in our algorithm is the updating
of slacks in the event that the SQP step does not satisfy the linearized constraints
well enough, i.e., ' k is not small enough. This can occur when (QP ) is inconsistent
or when a trust region bound is encountered during the solution of (QP ). In this case
our slack variable updating scheme would ensure that non-negative slacks remain non-
negative, but the direction may not be one of descent. We resolve this dilemma by
opting for descent, i.e., computing q k with ' replacing any negative slacks
using the following rule:
If z k+1
z k+1
ae ffl Mach g i
\Gammag i
where ffl Mach is machine epsilon. This is sometimes referred to as 'closing' the constraints
(see for example [33]).
4.9. Linearly Dependent Constraint Gradients. Linearly dependent constraint
gradients cause many theoretical and computational difficulties in constrained
optimization. In our theoretical algorithm we obtain convergence even when there are
linearly dependent constraint gradients provided the approximate multipliers do not
become unbounded. In practice, even though O3D has no difficulty in dealing with
this problem, evaluating the merit function and computing the least squares approximation
to the Lagrange multipliers become problematical. Computational experience
shows we solve many problems with degeneracy in the constraints. Simply maintaining
slacks to to be positive as described above allows us to factor the crucial matrices
and continue with the algorithm. However, the algorithm failed to solve some problems
that had a large amount of degeneracy in the linearized constraint matrix. This
was, of course, problem dependent but it was observed that the current implementation
can usually solve problems where up to 25 percent of the constraint gradients are
linearly dependent. This degeneracy causes the performance of the merit functions
to deteriorate. In particular, the least squares approximation to Lagrange multipliers
seems to be especially poor, resulting in only very small steps being allowed, even
close to the solution.
5. Numerical Results. The modified algorithm was coded in Fortran and is
installed on a SPARCstation 10 using IEEE floating point arithmetic (64 bit). The
current implementation is being used to solve a wide variety of medium to large
scale problems. In this section we report the results of a set of performance tests
designed specifically to answer questions about the trust region strategy and the
procedure to update the penalty parameter, d. We conclude the section with the
results of our algorithm applied to some test problems that are publicly available. We
emphasize that all of the problems were solved with the same default settings of the
parameters, (see Table 1), i.e., no attempt was made to pick parameter settings to
optimize performance on individual problems.
Although in many of the applications some analytic derivatives were available,
no use of analytic derivative information was used in these numerical experiments.
When possible, first and second derivatives were computed using forward and central
finite differences respectively. A costly one-time calculation provided a zero/non-zero
stencil of the Hessian of the Lagrangian and the Jacobian matrix of the constraint
function. These stencils were then used for the duration of the solution process.
For some problems, these finite difference approximations are not convenient to use.
This can be the case with control problems governed by partial differential equations
(see [29] or [30]). If the partial differential equation is solved using a finite element
method, with piecewise linear elements, then evaluating the derivative of the objective
Parameter Value
z Mach )

Table
Numerical values of default parameters
function with respect to the control variables can be quite cumbersome. In such
cases, which occurred in the control problems in our test suite, one can approximate
the first derivatives of the objective function by solving an adjoint problem with
a computational cost comparable to one function evaluation. (For examples, see
[22].) The objective function portion of the Hessian of the Lagrangian can then be
approximated with forward finite differences.
A set of eight problems was chosen as the first test suite. These problems ranged
in size from 500-1000 variables and from 1000-2000 constraints. The first four are
relatively straightforward nonlinear programming test examples, while the last four
are from actual applications: two discretized control problems, a density estimation
problem from statistics, and a "molecular distance" problem. A more complete description
of these problems is found in the Appendix. The problems all have nonlinear
inequality constraints and exploitable sparsity. Problem 4 (NLP4) was designed to
have a controllable percentage of linear dependency in the constraint gradients to
demonstrate any weaknesses in the algorithm associated with this difficulty. We ran
three versions of our algorithm on each problem; using a positive definite modification
of the Hessian matrix, as discussed in Section 4, with and without the trust region
strategy and using the unmodified Hessian with the trust region. (Using the unmodified
Hessian results in failure in most cases if no trust region strategy is employed.)
In addition, each problem was run from two starting points; one, labeled "c", which
was close to the solution in the sense that each of the variables was of the same order
of magnitude as in the solution and a distant start, labeled "f".
The results of the numerical tests on these problems are summarized in Tables
1-3. The first two columns of each table gives the number of SQP iterations ("nl-i")
and the total number of O3D iterations ("qp-i"). The next two columns contain the
stopping criterion that was met and the value of the gradient of the Lagrangian at the
solution. Unless the algorithm failed, (which is denoted by "Failure" in the tables)
feasibility condition (4.8) was satisfied for all solutions. The stopping criterion is
denoted by either a "1" or a "2" depending on whether (4.9) or (4.10) was satisfied.
If both were conditions were satisfied, a "3" appears in the column. The remaining
columns give information about the values of the parameter d for each run; columns
five through eight giving the initial, maximum, minimum, and final values of this
parameter and the final column giving the last iteration at which d was changed.
The results of the tests illustrate that using the unmodified Hessian with the trust
region was most effective in reducing the number of O3D iterations and the number
of SQP iterations. The trust region strategy prevented long, unprofitable steps from
being generated when far from the solution and the use of the unmodified Hessian
allowed the trust region to become inactive near the solution thus allowing rapid
local convergence. Requiring the Hessian to be positive definite often precluded rapid
local (q-superlinear) convergence and, when used in conjunction with the trust region
strategy, resulted in the trust region's being active close to the solution.
The results also show that the value of the parameter d varied over several orders
of magnitude. The procedures discussed in Section 4 that allowed the value of d to
increase or decrease greatly enhanced the algorithm; earlier tests using either a fixed
value of d or only allowing a reduction in d yielded inferior results.
Another modification in our algorithm, not reflected in the table or included in
the description in the preceding section, was made to force the O3D algorithm to take
a minimum number of steps. We found that when the trust region radius  became
small the algorithm would sometimes exit O3D after only one iteration, resulting in
a poor step direction. This poor step would result in a further decrease in  , and
eventually the algorithm would fail. When we required a minimum number of steps
to be taken in O3D (our choice was 7) this problem disappeared.
Recently a collection of test problems has become available for the testing and
comparing of optimization algorithms, (see [13]). The Constrained and Unconstrained
Testing Enviroment (CUTE), are quickly becoming standards with which researchers
can establish the viability and effectiveness of their numerical algorithms. These problems
are replacing the smaller and well scaled test problems of Hock and Schittkowski
[25] and Schittkowski [32] which were not intended to be used to test large scale al-
gorithms. Our results on the CUTE test problems are summarized in Tables 6, 7
and 8. These problems were solved to the same stopping conditions as the problems
above. Likewise, the same table format was used to present these numerical results.
For detailed description of these problems, structure, motivation, and sources see [9].
While it appears that the CUTE test problem set is rich in both large and small
scale unconstrained and equality constrained test problems, at present there are not
many large scale problems that include inequality constraints (and particularly non-linear
inequality constraints). We chose problems that reflected the class of problems
our algorithm was designed to solve. At least one inequality constraint was present
in each problem. The number of variables and/or constraints was large enough so
that the exploitation of special sparsity structure was important. The problems we
selected from CUTE to report on were "CORKSCREW, MANNE, SVANBERG" and
"ZIGZAG". The associated problem sizes are recorded in Table 5.
It is worth commenting that much of the machinery developed in this paper deals
with effectively handling nonlinear inequality constraints. The performance of our
algorithm on the CUTE test problem set is, therefore, slightly deceiving since many
of the constraints in these problems are simple bounds on the primal variables or
purely linear. (For instance, approximately 83% of the constraints in CORKSCREW,
50% of the constraints in MANNE, and 66% of the constraints in ZIGZAG were linear
and many of them were equality constraints). Although these caused no problem for
our algorithm, the structure of these constraints was not completely exploited and
the extra machinery of our code resulted in an overhead with no performance benefit.
Clearly an algorithm designed specifically to deal with linear equality constraints
should outperform our algorithm on these problems. The problem on which our
algorithm appeared to perform best was SVANBERG, a problem with only inequality
constraints (and a substantial number of them are nonlinear).
We succeeded in solving all four problems with a reasonable number of inner and
outer iterations. However, many of our algorithmic enhancements contributed little
to the solution process. The measure of distance to feasibility (the j-tube strategy),
the nonmonotone updating of penalty parameter d, and the trust region strategy were
essentially dormant during the solution process regardless of the iterates' proximity
to the solution or to feasibility. In fact, the only evidence of our enhancements on the
small number of CUTE test problems that we solved occurred when d was decreased
slightly while solving the problem MANNE employing modified Hessians with a trust
region strategy (see the third and fourth rows of Table 7). It is noteworthy that the
iterates that resulted from solving this problem with the penalty parameter artificially
held fixed at were identical to iterates that resulted for the adjusted d solution.
This appears to illustrate that in this case the adjustment of d was purely superficial.
Truss - c 103 3561 1 4.4e-8 9.87e-1 1.01e00 9.57e-2 9.57e-1 100
Molec-c 37 1376 1 9.8e-9 9.88e-1 1.21e1 1.17e-2 9.81e-1 34
Molec-f

Table
Modified Hessians with no trust region
BCHeat -c 257 2898 1 8.1e-8 9.18e-1 4.15e00 1.38e-1 4.10e00 249
BCHeat -f 289 3071 1 9.9e-8 1.00e00 3.74e00 2.44e-1 3.89e00 281
Molec-c
Molec-f 44 621 1 6.6e-8 1.00e00 1.48e00 7.39e-2 9.52e-2 38

Table
Modified Hessians with trust region
Truss-c 94 2242 1 3.9e-8 9.87e-1 1.03e00 1.26e-1 9.11e-1 87
Molec-c
Molec-f

Table
Un-modified Hessians with trust region
Problem Variables Constraints
CORKSCREW 96 159
ZIGZAG 304 1206

Table
Minimization parameters
Problem nl-i qp-i conv krx lk1 d0 maxd mind final d last d-cha

Table
Modified Hessians with no trust region
Problem nl-i qp-i conv krx lk1 d0 maxd mind final d last d-cha

Table
Modified Hessians with trust region
Problem nl-i qp-i conv krx lk1 d0 maxd mind final d last d-cha

Table
Un-modified Hessians with trust region
6. Future Directions. In this paper we have discussed in some detail an SQP
algorithm for solving large scale nonlinear problems. The numerical results with
default parameter settings indicate that the procedures that we have implemented
are robust, effective, and efficient; the convergence theory in [4] provides a sound
theoretical basis for the procedure. Nevertheless, there are several areas in which the
techniques used here can be improved to allow the solution of larger and more difficult
problems.
Algorithmically, we observe that the current implementation requires the factorization
of both (rg T rg+Z) and (rgrg T ), the latter in O3D. While the sparse matrix
package makes this reasonable for the problems that we have currently considered, it
is clearly expensive to maintain both.
The results reported here use analytic or finite difference Hessian approximations.
An examination of the details of O3D reveals that a limited memory BFGS or limited
memory SR1 could be readily incorporated into the code. We have done some
experimentation with such techniques; the results will be reported elsewhere [26].
Many of the problems that we have seen have been degenerate and this significantly
slows the convergence of the method. The primary culprit is the extremely
poor multiplier estimates provided by the least squares procedure. Improvements in
this area are certainly required.
In some problems (not reported here) that have nonlinear equality constraints,
we have occasionally observed significant difficulty in trying to satisfy the linearized
equality constraints, i.e., in completing Phase I. In these cases we have had some success
in relaxing the constraints [26]. In the context of O3D, this can be accomplished
by simply fixing the artificial variable at some positive value and continuing the O3D
iterations. In this approach, we often find that O3D converges and the "recentering"
procedure mentioned in Section 2 has led to further improvements. The theory in [4]
supports these ideas. The details will, again, be reported elsewhere.



--R

An interior-point method for general large scale quadratic programming problems
An interior point method for linear and quadratic programming problems
A merit function for inequality constrained nonlinear programming problems

A family of descent functions for constrained optimization


A truncated SQP algorithm for large scale nonlinear programming problems
Cute: Constrained and unconstrained testing environment
Functional and numerical solution of a control problem originating from heat transfer
On exact and approximate boundary controlla- bilities for the heat equation
Large scale numerical optimization: Introduction and overview
Lancelot: A Fortran Package for Large-Scale Nonlinear Optimization

SIAM Journal on Numerical Analysis


Globally convergent inexact Newton methods
A robust trust region algorithm with nonmonotonic penalty parameter scheme for constrained optimization
New York
Constrained nonlinear pro- gramming
Numerical Methods for Nonlinear Variational Problems
Molecular conformations from distance matrices
Optimal signal sets for non-gaussian detectors
Test Examples for Nonlinear Programming Codes
The Use of Optimization Techniques in the Solution of Partial Differential Equations from Science and Engineering
The solution of the metric stress and sstress problems in multidimensional scaling using Newtons method
Numerical solution of a nonlinear parabolic control problem by a reduced SQP method
Optimal Control of Systems Governed by Partial Differential Equations

Computing a trust region step
More Test Examples for Nonlinear Programming Codes
On the role of slack variables in quasi-Newton methods for constrained optimiza- tion
Nonparametric Probability Density Estimation
Nonparametric Function Estimation
where it was used to illustrate separability in nonlinear programming.
"distance matrices"
--TR

--CTR
Thomas F. Coleman , Jianguo Liu , Wei Yuan, A New Trust-Region Algorithm for Equality Constrained Optimization, Computational Optimization and Applications, v.21 n.2, p.177-199, February 2002
E. Bradley , A. OGallagher , J. Rogers, Global solutions for nonlinear systems using qualitative  reasoning, Annals of Mathematics and Artificial Intelligence, v.23 n.3-4, p.211-228, 1998
Rubin Gong , Gang Xu, Quadratic surface reconstruction from multiple views using SQP, Integrated image and graphics technologies, Kluwer Academic Publishers, Norwell, MA, 2004
interior point;large scale;trust region;nonlinear programming;SQP;merit function