An Interior-Point Algorithm for Nonconvex Nonlinear Programming.
The paper describes an interior-point algorithm for nonconvex nonlinear programming which is a direct extension of interior-point methods for linear and quadratic programming. Major modifications include a merit function and an altered search direction to ensure that a descent direction for the merit function is obtained. Preliminary numerical testing indicates that the method is robust. Further, numerical comparisons with MINOS and LANCELOT show that the method is efficient, and has the promise of greatly reducing solution times on at least some classes of models.
INTRODUCTION
In this paper, we describe modifications that we used to convert the quadratic programming
(QP) solver LOQO into a general nonconvex nonlinear programming solver (of the
same name). As a code for quadratic programming, LOQO implements an interior-point
method. The complete details of the QP implementation can be found in [25]. In view of
the dramatic success of modern interior-point methods for linear and quadratic program-
ming, our guiding philosophy was to modify the QP interior-point method in LOQO as little
as possible to make a robust and efficient general nonconvex nonlinear optimizer.
For notational simplicity, we begin by considering the following nonlinear programming
problem:
subject to h i (x) # 0,
(1)
where x is a vector of dimension n and f (x) and the h i (x) are assumed to be twice continuously
differentiable. This is a simplification of the general nonlinear programming
problem, which can include equality constraints and bounds on the variables. In fact, the
method developed in this paper and its implementation in LOQO handles all of these cases
efficiently, and the exact way in which this is accomplished is discussed in detail later. For
1991 Mathematics Subject Classification. Primary 90C30, Secondary 49M37, 65K05.
Key words and phrases. Nonlinear programming, interior-point methods, nonconvex optimization.
Research of the first author supported by NSF grant CCR-9403789 and by ONR grant N00014-98-1-0036.
Research of the second author supported by AFOSR grant F49620-95-1-0110.
the present, we consider this version of the problem, as it greatly simplifies the terminol-
ogy, and the extension to the more general case is quite straightforward.
The interior-point approach taken in this paper is described as follows. First, add slack
variables w i to each of the constraints (1), reformulating the problem as
subject to h(x) -w= 0,
(2)
where h(x) and w represent the vectors with elements respectively. We then
eliminate the inequality constraints in (2) by placing them in a barrier term, resulting in the
problem
subject to h(x) -
where the objective function
is the classical Fiacco-McCormick [8] logarithmic barrier function. The Lagrangian for
this problem is
y,
and the first-order conditions for a minimum are
#w
where W is the diagonal matrix with elements is the vector of all ones, and #h(x)
is the Jacobian matrix of the vector h(x). We now modify (5) by multiplying the second
equation by W , producing the standard primal-dual system
where again Y is the diagonal matrix with elements y i . Note that the second equation
implies that y is nonnegative, which is consistent with the fact that y is the vector of
Lagrange multipliers associated with what were originally inequality constraints.
The basis of the numerical algorithm for finding a solution to the primal-dual system
(6) is Newton's method, which is well known to be very efficient for linear and convex
quadratic programming. In order to simplify notation and at the same time highlight connections
with linear and quadratic programming, we introduce the following definitions:
and
AN INTERIOR-POINT ALGORITHM FOR NONCONVEX NONLINEAR PROGRAMMING 3
The Newton system for (6) is then
#x
#w
#y
3 5 .
The system (7) is not symmetric, but is easily symmetrized by multiplyingthe first equation
by -1 and the second equation by -W -1 , yielding
#x
#w
#y
-#
where
y,
Note that # measures primal infeasibility. By analogy with linear programming, we refer
to # as dual infeasibility. Also, note that # , and # depend on x , y, and w even though
we don't show this dependence explicitly.
It is (8), or a small modification of it, that LOQO solves at each iteration to find the
search directions #x , #w, #y. Since the second equation can be used to eliminate #w
without producing any off-diagonal fill-in in the remaining system, one normally does this
elimination first. Hence, #w is given by
and the resulting reduced KKT system is given by
-H (x, y) A T (x)
#y
.
The algorithm then proceeds iteratively from an initial point x (0) , w (0) , y (0) through a
sequence of points determined from the search directions described above:
For linear programming different steplengths, # (k)
d , are used with the primal and
dual directions, whereas for nonlinear programming a common steplength is employed.
1.1. Formulas for search directions. While LOQO solves the reduced KKT system (12),
it is nonetheless useful for subsequent analysis to give explicit formulas for #x and #w.
Let
y,
denote the dual normal matrix. In the following theorem, we drop the explicit indication
of the dependence of A, N , and # f on x , y, and w.
Theorem 1. If N is nonsingular, then (8) has a unique solution. In particular,
4 ROBERT J. VANDERBEI AND DAVID F. SHANNO
Remark. The formulas for #x and #w involve three terms each. In both cases, the first
term can be viewed as an optimality direction, the second term as a centrality direction,
and the third term as a feasibility direction.
Proof. Solving the second block of equations in (12) for #y and eliminating #y from the
first block of equations yields a system involving only #x whose solution is
Using this formula, we can then solve for #y and finally for #w. The resulting formula
for #w is:
From the definitions of # and # , it follows that
Finally, we use this formula to eliminate # and # from (16) and (17). The formulas (14)
and (15) for #x and #w then follow easily.
1.2. Critical features. The critical features of the algorithm are the choices of # (k) at
each iteration and the modification of the system (8) in order to find a local minimizer of
(1). For linear and convex quadratic programming, modification of the linear system is
never required except possibly to deal with problems of numerical accuracy, and the step
length at each iteration is determined by a simple ratio test. (See, for example, Lustig,
Marsten and Shanno [18] and Vanderbei [25].) For convex nonlinear programming, again
the linear system need not be modified, but the method for choosing # (k) at each iteration
becomes more complex, as it is well known that for general convex nonlinear problems,
with a poor initial estimate Newton's method may diverge. In order to achieve convergence
to a solution to the system (6), El Bakry et al. [7] introduced the merit function
and showed that for a proper choice of -, there exists a sequence of step lengths {# (k)
such that # 0 decreases monotonically and the iterates converge to a point at which #
provided the Jacobian of the system (6) remains nonsingular. Shanno and Simantiraki [21]
tested a variant of this algorithm on the Hock and Schittkowski set of test problems [16]
and found that while the algorithm was often efficient, it also often converged to local
maxima or saddle points, all of which satisfy the first-order necessary conditions. Further,
the Jacobian of the system (6) sometimes becomes singular, causing the algorithm to fail.
Thus it became apparent that a better merit function, and a way of ensuring a nonsingular
Jacobian, would be necessary for a successful general algorithm. Our approach to this is
discussed in detail in the next section.
1.3. Related work. It has recently come to our attention that Lasdon et al. [12] have
also been studying this basic algorithm. While their algorithm differs from ours in many
particulars, the general approach is similar. For related work on primal-dual interior-point
methods see [10], [2], [13], [5]. An alternative barrier approach to inequality-constrained
problems is discussed in [1].
AN INTERIOR-POINT ALGORITHM FOR NONCONVEX NONLINEAR PROGRAMMING 5
2. ALGORITHM MODIFICATIONS FOR CONVEX OPTIMIZATION
As mentioned in the previous section, the search directions given by (8) together with a
steplength selected simply to ensure that the vectors w and y remain component-wise non-negative
yield an efficient and robust algorithm for convex quadratic programming. How-
ever, for nonquadratic convex optimization problems, further reduction in the steplength
may be necessary to guarantee convergence, as is trivially illustrated by the unconstrained
minimization of the univariate function f using an initial estimate x (0)
with |x (0)
| > 1. Merit functions are used to guide one in deciding how much to shorten
the steplength. In this section we describe the specific merit function that we have implemented

2.1. The merit function. Merit functions for equality constrained nonlinear programming
have been the subject of a great deal of research over the past twenty-odd years. The idea of
a merit function is to ensure that joint progress is made both toward a local minimizer and
toward feasibility. This progress is achieved by shortening the steplength along the search
directions defined by (8) as necessary so that sufficient reduction in the merit function is
made. One possible merit function is
This merit function is exact, which means that there exists a # 0 such that, for all # 0 , a
minimizer of (19) is guaranteed to be feasible and, under general conditions, a local minimizer
of the original problem. While exactness is a useful property, the nondifferentiability
of # 1 -norms can cause difficulties with numerical algorithms. The merit function defined
by
is the penalty function for equality constrained nonlinear programming studied by Fiacco
and McCormick [8]. It is differentiable everywhere. However, it has the theoretical disadvantage
of requiring # to tend to infinity to ensure convergence to a feasible point, which
again is hoped to be a local minimizer of the original problem. In spite of this apparent
disadvantage, in practice we have had no such difficulty, and hence have chosen to use it.
The l 2 merit function (20) when applied to problems of the form (3) is
#,-
log(w
(recall that #(x, h(x)). The following theorem shows, among other things, that
for large enough #'s the search directions defined by (8) are descent directions for #,-
whenever the problem is strictly convex.
Theorem 2. Suppose that the dual normal matrix N is positive definite. Then, the search
directions have the following properties:
#w b  T
#x
#w
(2) There exists # min # 0 such that, for every # min ,
#w#,-
#x
#w
In both cases, equality holds if and only if (x, w) satisfies (6) for some y.
6 ROBERT J. VANDERBEI AND DAVID F. SHANNO
Proof. It is easy to see that
From the expressions for #x and #w given in Theorem 1, we then get
#w b  T
#x
#w
Assuming that and that N is positive definite, we see that
#w b  T
#x
#w
which completes the proof of the first property.
For the second property, first note that the merit function is the barrier function plus a
constant times a measure of infeasibility:
#,-
We now address the infeasibility term. It is easy to check that
#w #(x, w)# 2
I  # .
From this explicit expression for the gradient of the infeasibility term, it follows that
#x
#w
Combining the inner products using the barrier function and the infeasibility term, we get
#w#,-
#x
#w
.
Suppose that (x, w) is not feasible (#= 0). For (#x, #w) to fail to be a descent
direction for the merit function, it is necessary that
-e
When this is the case, setting
ensures that (#x, #w) is a descent direction for #,- for every # min .
AN INTERIOR-POINT ALGORITHM FOR NONCONVEX NONLINEAR PROGRAMMING 7
Now suppose that (x, w) is feasible (i.e., 0). In this case, the barrier function's
inner product and the merit function's inner product agree. Furthermore, either (#x, #w)
is a descent direction for both or
In the latter case, we use the positive definiteness of N to conclude that
Introducing a dual variable e, we see that this last equation together with
the assumed feasibility comprise precisely the first-order optimality conditions (5) for the
barrier problem.
In LOQO, # is initialized to 0 and is unchanged as long as (#x, #w) is a descent direction
for #,- . When (#x, #w) fails as a search direction, then # is calculated from
using This algorithm has proved quite satisfactory in practice, as # is
changed very few times during the calculations on problems we have tried so far, and on
many problems remains 0. While this may seem surprising, recall that the search direction
is the Newton direction to try to find a feasible first-order point, and generally moves
jointly toward a minimizer and a feasible point. To ensure that we will not simply approach
an infeasible optimum very slowly, if the step along the search vector becomes very small,
# is increased by a factor of 10. The current version of the code starts the next iteration
from the new point with the increased value of # . In future work, we will test the efficiency
gained by simply using the new # with the old search directions.
Thus at each step an # max is chosen by the standard ratio test to ensure that nonnegative
variables remain nonnegative. In the context of the discussion of this section, this means
that, for must remain strictly positive so
Then the interval [0, # max ] is searched by successive halving, using merit-function evaluations
only, for a value of # that produces a reduction in the merit function satisfying an
Armijo condition.
We conclude this section with two comments. First, the above discussion assumed that
N is positive definite. Subsequent sections will demonstrate how the algorithm can be
modified to relax this assumption.
Second, the method listed here is admittedly ad hoc. A convergence proof for a method
of this kind will certainly require a more carefully developed algorithm in terms of selection
of the parameter # . Yet in preliminary testing this method works surprisingly well,
while attempts at more rigorous algorithms were less successful. Thus it can be viewed
similarly to the recent method of Fletcher and Leyffer [9], which essentially accepts any
point that improves either optimality or infeasibility. While much remains for further study,
our experience to date is that Newton's method should be hindered as little as possible.
3. ALGORITHM MODIFICATIONS FOR NONCONVEX OPTIMIZATION
Using the merit function described in the previous section to guide the selection of a
steplength (together with the stabilization technique of the next section) yields an efficient
and robust algorithm for convex optimization. For nonconvex optimization problems, there
is another important issue to consider, namely, that the matrix N(x, y,
8 ROBERT J. VANDERBEI AND DAVID F. SHANNO
A T (x)W may fail to be positive semidefinite. In this section, we discuss the
impact of indefiniteness and describe the algorithmic changes that we made to address this
issue.
3.1. The search direction. In Theorem 2 we showed that the search directions (#x, #w)
have desirable descent properties for both the barrier function and the merit function provided
that N is positive definite. When N is indefinite, the algorithm might still converge
to something, but it might not be a local minimum for the problem one wishes to solve.
Consider, for example, the problem of minimizing the concave function 4x(1 - x) subject
to the bound constraint 0 # x # 1. The algorithm presented so far when applied to
this problem and initialized with 0.4 # x (0)
# 0.6 has the very undesirable property of
converging to the global maximum at 0.5.
These results suggest replacing H (x, y) in the definition of the search directions with a
diagonal perturbation thereof:
If # is chosen so that N is positive definite, then Theorem 2 applies and ensures descent
properties for the barrier and merit functions. Of course, the following question naturally
arises: what is the impact of such a perturbation on other important properties such as
reduction in primal and dual infeasibility? The following theorem addresses this question.
Henceforth, we assume that the search directions are computed using -
H instead of H .
Theorem 3. Let
y).
Then
y.
Proof. We begin with -
#:
where the last equality follows from the third block of equations in (8). For -
# , the analysis
is similar but one must use the first block of equations in (8) with -
y
AN INTERIOR-POINT ALGORITHM FOR NONCONVEX NONLINEAR PROGRAMMING 9
Finally, the analysis of w T y follows closely the analogous development for linear pro-
gramming. It starts as follows:
Then, we rewrite the linear term in t as follows:
y.
Substituting (24) into (23) we get the desired result.
As for linear programming, - is usually chosen so that 0 < # < 1 and therefore w T y
decreases provided that the step length is short enough. However, if # > 0, Theorem 3
shows that dual infeasibility may fail to decrease even with arbitrarily small steps. This
seems a small price to pay for the other desirable properties and, even though we don't
have a proof of convergence for the perturbed method, empirical evidence suggests that #
is zero most of the time (using the rules given shortly) and that the dual infeasibility does
eventually decrease to zero once the iterates are in a neighborhood of the optimal solution
(where N should be positive definite).
In our nonlinear LOQO, we pick # > 0 whenever necessary to keep N positive defi-
nite. By doing this, we guarantee that the algorithm, if convergent, converges to a local
minimum.
The value of # is computed as follows. We first do an L DL T factorization of a symmetric
permutation of the reduced KKT matrix in (12). The permutation is computed at the
beginning based solely on structure/fill-in considerations; see [25] for details. If H (x, y)
is positive definite, then the reduced KKT matrix is quasidefinite as defined in [24] and
hence, as shown in [24], the factorization is guaranteed to exist. Furthermore, it has the
property that every diagonal element of D associated with an original diagonal element of
must be positive and every element associated with an original diagonal element
of -WY -1 must be negative. Hence, after factoring, we scan the diagonal elements of D
looking for elements of the wrong sign and noting the one with largest magnitude. Let # 0
denote the largest magnitude. If all of the diagonal elements have the correct sign, then
must be positive definite and no perturbation is required. Otherwise, we do an
initial perturbation with . For the trivial case of a 1 - 1 matrix, such a perturbation
is guaranteed to produce a positive definite matrix. However, for larger matrices, there
is no such guarantee. Therefore, we factor the perturbed matrix and check as before. If
the perturbation proves to be insufficient, we keep doubling it until a perturbation is found
that gives a positive definite matrix. If, on the other hand, the initial perturbation is sufficiently
large, then we do successive halving until we find a perturbation that is too small.
Finally, we double once. It is not clear whether this halving process is necessary. We plan
to investigate this issue in the future. The method described here is very simple, and can
undoubtedly be improved, but it has in general worked very well in practice. Note that as
w and y are always positive, only H (x, y) ever needs to be modified.
4. BOUNDS AND RANGES
As noted in the introduction, the general nonlinear programming problem may have
equality constraints as well as inequality constraints. Furthermore, simple bounds on the
variables, while they may be considered as inequality constraints, are generally handled
separately for efficiency. All of these cases were treated specifically in the quadratic solver
LOQO and are described in [25]. Our nonlinear modification continues to accept all such
formulations. We describe the algorithmic issues only briefly here.
A more general form of inequality constraint is a range constraint of the form
This is easily converted to a system of constraints
Equality constraints are treated by simply declaring them to be
range constraints with r Note that this does not increase the number of nonlinear
constraints, but does add a linear constraint and an extra variable. In [23], Vanderbei shows
how to reduce the resulting KKT system to one identical in size and character to (12), but
with modifications to the diagonal of matrix W -1 Y .
Bounds on the variables are handled similarly. The algorithm assumes that each variable
x j has (possibly infinite) upper and lower bounds:
where
-#< u j #.
As in interior-point methods for linear programming, the finite inequalities are converted to
equalities by adding slack variables, and the resulting system of equations is again reduced
to one of the exact form and size of (12) with the only change being the addition of a
diagonal matrix to the H (x, y) block of the reduced KKT matrix. One should note that,
in contrast to the lower-bound shift employed by others, we treat lower bounds exactly
the same way that we treat upper bounds; such a symmetric treatment of upper and lower
bounds was first suggested by Gill, et al., in [14].
The only addition to the discussion of this section is that all nonnegative variables added
through either range or bound constraints must be included in the logarithmic barrier term
in the merit function and any extra linear equalities are included in the penalty term. Thus
the final form of the reduced KKT system solved is
#y
,
is an n - n nonnegative diagonal matrix, Em is an positive diagonal
matrix and # 1 and # 2 are appropriately modified right-hand sides; see [25] for explicit
expressions. If every variable has either a finite upper bound or a finite lower bound, i.e. if
there are no free variables, then E n too is a positive diagonal matrix.
AN INTERIOR-POINT ALGORITHM FOR NONCONVEX NONLINEAR PROGRAMMING 11
4.1. Free variables and stabilization. A diagonal element of the nonnegative matrix E n
introduced in the previous section is positive if the corresponding variable has a finite upper
or lower bound and it is zero if the variable is a free variable. Therefore as things stand,
for linear programming problems having free variables, the reduced KKT matrix is not
quasidefinite. To get a quasidefinite matrix even for this case, we write each free variable
x j as a difference of two nonnegative variables:
Such a splitting of free variables is a common trick but the next step is usually to eliminate
them from the problem. We do not do this. Instead, we retain each free x j and add the
constraint shown above. We then do the algebra to reduce the system again to the usual
reduced KKT system. The result is exactly the same as before except that now E n is strictly
positive in every component, even those associated with free variables. This improves the
conditioning of the reduced KKT system and thereby helps stabilize the accuracy of the
factorization. Of course, the new nonnegative variables, t j and g j , must be incorporated
in the logarithmic barrier term in the merit function and the linear equations in the penalty
term.
This implicit handling of free variables ensures that the matrix in (25) is quasidefinite
whenever the problem is a convex optimization problem (in particular, whenever it is a
linear programming problem). Furthermore, by adding something positive to H (x, y) we
increase the likelihood that the system is quasidefinite even when the problem is nonconvex
and hence we greatly reduce the number of times that we need to add diagonal perturbations
to H (x, y). In fact, on the Hock and Schittkowski test set on which we report in
Section 6 LOQO took a total of 2174 Newton iterations with 3298 matrix factorizations for
an average of 1.52 factorizations per iteration. As previously stated, our search algorithm
for # is very simplistic and could lead to significantly more factorizations than a more
sophisticated version. However, these results seem to verify that the addition of the positive
diagonal matrix E n and the semidefinite matrix A T (x)W greatly
reduces the number of perturbations required.
For full details of the computation of E n and the modification of the right-hand side, as
well as full details from the previous section, the reader is referred to [23] and [25].
5. IMPLEMENTATION DETAILS
This section provides implementation details for the nonlinear version of LOQO that we
have coded and tested. In particular, we deal with some specific algorithmic choices. In the
next section we give the results of computational testing, with comparison to MINOS [20]
and LANCELOT [6] on a test set that includes all of the well-known Hock and Schittkowski
problems [16] as well as several other large-scale real-world problems.
5.1. Choice of the barrier parameter -. So far, we have not specified any particular
choice for the barrier parameter -. Most theoretical analyses of interior-point algorithms
choose - as
In [7], El Bakry et al. show that for this choice of -, their algorithm
converges to a zero of the merit function (18). Computational experience has shown that
the algorithm performs best when the complementary products w i y i approach zero at a
uniform rate. We measure the distance from uniformity by computing
Clearly, only if w i y i is a constant over all values of i. When far
from uniformity, a larger - promotes uniformity for the next iteration. Consequently, we
use the following heuristic for our choice of -, which has proven very effective in practice:
denotes the steplength parameter, which by default is 0.95, and # is a
settable scale factor, which defaults to 0.1.
5.2. The initial point. Nonlinear programming traditionally requires that a starting point
be given as part of the problem data, and comparative numerical testing is done using these
traditional starting points. For an interior-point code, more is required, however, as both
slack variables and split variables are added to the problem. Thus even with initial values
for the x j 's, initial values for the w i 's, y i 's etc. must be determined by the program. For
slack variables, for instance the w i 's, given an x (0) we can compute h i
There are two difficulties with (26). First, if x (0) is not feasible, (26) gives an initial negative
value to some of the w i 's, which is not allowed for interior-point methods. Second,
even if x (0) is feasible, it may lie on the boundary of the feasible region or so close to it
that some of the initial w i 's are very close to zero and progress is hindered. Much as in
interior-point methods for linear programming, it is necessary to specify a # > 0 so that
all initial values of variables constrained to be nonnegative are at least as large as # . Hence
w i is initially set as follows:
An analogous method can be used to set the slacks on the range constraints used to handle
equality constraints, and also the slack variables that transform simple bounds into
equalities.
The latter brings up another interesting point. Some of the standard test problems have
simple bounds on the variables and give an initial point that violates these bounds. Such an
initial point seems unrealistic and can be computationally dangerous, as sometimes bounds
on variables are used to ensure that the function value can actually be calculated, such as
enforcing appears somewhere in the calculations. In view of this, we have
added an option to the code that can force bounds to be honored. In this case, any initial
point that lies outside the bounds is altered to lie inside the bounds and is set exactly #
from the bound if there is only one bound and, if there are both upper and lower bounds,
is set to a 90%-10% mixture of the two bounds with the higher value placed on the nearer
bound. Finally, the split parts associated with a free variable x are set so that
the difference is equal to x j and the smaller of the two is equal to # .
The remaining point of interest is the choice of # . For linear programming,
has proved computationally efficient. For the nonlinear problems we have tested to date,
however, proved by far the best choice. A good problem-dependent choice for
# remains a topic for further study.
AN INTERIOR-POINT ALGORITHM FOR NONCONVEX NONLINEAR PROGRAMMING 13
5.3. Other issues. Issues such as stopping rules and matrix reordering heuristics are handled
in exactly the same manner as documented in [23].
6. COMPUTATIONAL RESULTS
The LOQO implementation described in this paper, in addition to its ability to read
industry-standard MPS files for expressing linear and quadratic programming problems, is
interfaced with two of the most popular mathematical programming languages, AMPL [11]
and GAMS [4]. Currently, only AMPL can provide a solver with second-order information,
i.e. Hessians, and therefore all of our testing was performed using AMPL together with
LOQO. Many other solvers are also interfaced with AMPL, including MINOS [20] and
LANCELOT [6]. For the tests described here, we used LOQO version 3.10 (19971027),
MINOS version 5.4 (19940910), and LANCELOT version 20/03/1997. Both MINOS and
LANCELOT were run using default options although we did have to increase the number
of superbasics and the maximum number of major/minor iterations on some of the larger
problems.
All tests were performed on an R4600 SGI workstation with 160 MBytes of real mem-
Kbytes of data cache, and a 133 MHz clock. There are three components to the
stopping rule for LOQO: (i) primal feasibility, (ii) dual feasibility, and (iii) lack of duality
gap. The default rule declares a solution primal/dual feasible if the relative infeasibility
is less than 1.0e - 6 and declares the problem optimal if, in addition, there are 8 or more
digits of agreement between the primal and dual objective function values.
6.1. Hock and Schittkowski suite. Table 1 shows the solution times in seconds for the
Hock and Schittkowski set of test problems [16]. Table 2 shows how many interior-point
iterations LOQO used to solve each problem. Each solver was run with default parameters.
Since LOQO is still under development some of the defaults might change in the future.
Therefore, we list here those parameters that might change and give the current default
settings (for nonlinear problems):
. mufactor, called above, is set to 0.1.
. bndpush, called above, is set to 100.
. honor bnds, a boolean flag indicating whether bounds on variables are to be enforced
throughout the iteration process as opposed to only at optimality, is set to 1
(i.e., true).
All of the problems in the Hock and Schittkowski set are very small. For linear pro-
gramming, the simplex method is usually more efficient than interior-point methods on
small problems. One would expect the same for nonlinear programming. As expected,
MINOS, which is simplex-based, was faster on a majority of the problems. Nonetheless,
we were pleasantly surprised to see how well LOQO did compared with MINOS. Since all
of the times are small, direct time comparisons have limited value and so we have summarized
in Table 3 how many times each solver came in first, second, and third place.
Finally, we computed the total time for LOQO and MINOS on problems that both codes
solved to optimality. The total time for LOQO was 21.9 seconds whereas for MINOS it was
28.1 seconds. It should be noted here that MINOS uses only first derivatives. It is possible
that a true Newton variant of MINOS might improve speed. However, our experience to
date is that for small dense problems coded in AMPL the cost of evaluating second derivatives
is significant and may well offset any other improved algorithmic efficiency. It is our
understanding that the default LANCELOT uses exact second derivatives.
14 ROBERT J. VANDERBEI AND DAVID F. SHANNO
Time in Seconds Time in Seconds Time in Seconds
Name Minos Lancelot Loqo Name Minos Lancelot Loqo Name Minos Lancelot Loqo


1. Solution times on Hock-Schittkowski problems. Legend: (1)
could not find a feasible solution, (2) erf() not available, (3) step got too
small, (4) too many iterations, (5) could not code model in AMPL, (6)
unbounded or badly scaled, (7) core dump.
In order to interpret the results of Table 1, we note first that hs067 was not run with
any algorithm. This problem contains an internal fixed-point calculation which is difficult
to code in AMPL. LOQO failed to converge for only one case of the remaining test set,
hs013. This is a classic problem where the optimal point does not satisfy the KKT condi-
tions. As LOQO currently only looks for KKT points, it was unable to solve this problem.
Identifying optima that do not satisfy the KKT conditions as well as correctly determining
unboundedness and infeasibility remain topics for further study.
Another point of interest is the relatively high iteration count and solution time for
LOQO on hs027. These can be dramatically reduced by increasing # early in the iteration
sequence, but all methods we have tried to date that significantly help with hs027 hurt
enough other problems so as to be less efficient overall. Thus an optimal means of choosing
# , which one hopes would be part of a provably globally convergent algorithm, also
remains for further study.
AN INTERIOR-POINT ALGORITHM FOR NONCONVEX NONLINEAR PROGRAMMING 15
Name Iters Name Iters Name Iters Name Iters Name Iters


2. Iteration counts on Hock-Schittkowski problems.
Solver 1st 2nd 3rd
LOQO 44


3. Rankings
Nonconvex problems often have alternative optima, which is certainly the case with
some problems in the Hock and Schittkowski test set. In our experiments, LOQO found
alternative optima that were worse than the reported optima on 8 problems (hs002, hs016,
and better on 7 problems (hs047, hs088, hs089,
While finding optima that are worse is generally to be ex-
pected, finding optima that are better is quite a surprise because those reported by Hock and
Schittkowski represent the best optima known from years of experimentation with the test
set. For the suboptimal problems the desired optimal solution could generally be obtained
by altering # , thereby changing the initial solution.
In addition to the default value for # , we tried two other values. For
on two problems and is noticeably less efficient on others. For it fails on seven
problems and is again noticeably less efficient on others.
Name
markowitz2 201 1200 201200 200 4m43s 1m56s 13m35s
polygon2 195 42 766 880 0.7s *60m0s 1m1s
sawpath 198 5 784 25 2s 5s *8s
structure4 720 1536 5724 20356 2m40s *43m6.89s *
trafequil2 628 1194 5512 76 5.1s 5.3s 2m39s


4. Preliminary computational results for several application ar-
eas. An asterisk indicates that the solution obtained was either infeasible
or suboptimal. The double asterisk indicates not enough memory to
solve the problem.
We also did a run in which we did not force variables to remain within their bounds;
i.e., honor bnds=0. On this run, three problems in addition to hs013 failed to converge.
We did not find this very surprising for interior-point methods, where we feel that bounds
should be honored from the beginning.
As a final note, the iteration counts in Table 2 represent the total number of times first
and second partial derivatives were computed. A few extra function evaluations were used
in reducing the merit function, but as the time taken to calculate them is very small compared
to other components of the algorithm, we chose to report execution times rather than
function counts.
6.2. Large-scale real-world problems. We are assembling a collection of real-world
problems encoded in AMPL. The collection is available from the first author's web site
[22].

Table

4 gives a brief summary of problem statistics and computational results for
some of the larger problems that we have encountered. Most of these problems are quite
difficult. As the table shows, LOQO can be a very efficient and robust code on large difficult
problems when compared to the other solvers. We made naive attempts to set parameters
to appropriate non-default values in all codes. Clearly, we were better able to adjust parameters
in LOQO than in the other codes. Note that we set the time limit for each run at one
hour. Our inability to run LANCELOT due to memory limitations on structure4, is the
result of not knowing how to adjust the memory estimates made by AMPL when invoking
LANCELOT.
As these problems have not been previously part of any standard test set, we give here
a brief description of some of the interesting ones. Non-default LOQO parameters, when
used, are noted below with the problem descriptions.
6.2.1. Antenna Array Synthesis (antenna). An important problem in electrical engineering
is to determine how to combine the signals from an array of antennas so as to reinforce the
signal from a desired direction and suppress the signal from undesired ones. There are
various formulations of this problem, some of which fall under convex optimization. The
details of one such formulation are given in [17]. To solve this problem with LOQO, we set
parameters as follows: sigfig=5, bndpush=2, convex, inftol=1.0e-1.
6.2.2. Electrons on a Sphere (fekete2). Given n electrons placed on a conducting sphere,
the problem is to find a distribution of these electrons that minimizes the total Coulomb
AN INTERIOR-POINT ALGORITHM FOR NONCONVEX NONLINEAR PROGRAMMING 17
potential. This problem is nonconvex and has a large number of local minima. Problem
was solved with
6.2.3. optimization (markowitz2). The Markowitz model in portfolio optimization
seeks the portfolio of investments that optimizes a linear combination of expected return
and expected risk. Risk is modeled as the variance of the return. Problem markowitz2
is a separable convex quadratic formulation of the problem with linear constraints. Expected
reward and risk are computed using historical data.
6.2.4. Minimal Surfaces (minsurf). Given a domain in R 2 and boundary data, the minimal
surface problem is to find an interpolation of the boundary data into the interior of the domain
so that the surface so generated has minimal surface area. It is a convex optimization
problem.
6.2.5. Largest Small Polygon (polygon2). Given n, the problem is to find the n-sided polygon
of maximal area whose diameter does not exceed 1. This problem sounds trivial, but
Graham showed in [15] that the optimal hexagon is not the regular hexagon. To solve
polygon2 with LOQO, we set parameters as follows: convex.
6.2.6. Saw Path Tracking (sawpath). Given a list of points describing the center line of a
wood piece, the problem is to fit the best polynomial that minimizes the sum of the squares
of errors subject to three sets of side constraints: (1) the polynomial must go through
the first point, (2) the polynomial must have a specified initial slope, and (3) the radius
of curvature must never exceed a given value. This problem comes from F. Grondin of
6.2.7. Structural Optimization (structure4). In structural optimization, the problem is to
decide how to build a structure to minimize weight or compliance, subject to the constraint
that the structure can support one or more possible loadings. There are various
models, some of which are convex while others are not. Also, some of these models are
for pin-jointed truss-like structures while others are based on finite-element models. Problem
structure4 is to find an optimal bracket design obtained using a convex finite-element
model. The corresponding truss-like model is described in [3].
6.2.8. Traffic Equilibrium (trafequil2). The problem is to find flows through a network
that minimize a nonlinear function, called the Beckman objective.
6.3. Mittelmann's quadratic programming set. As a service to the Operations Research
community, Hans Mittelmann has done extensive testing of optimization software. We
summarize here some of his results [19]. He used an AMPL model to generate a number
of random, feasible, quadratic programming problems. These problems are determined by
specifying values for five parameters, n, m, p, s, and pf; see [19] for their definitions.
He used two versions of the random model generator; one that generates convex problems
and one that does not. The results for the convex problems are shown in Table 5 and
for the nonconvex problems the results are shown in Table 6.
The first table shows how efficient interior-point methods are compared to other methods
on convex quadratic programming problems. Regarding the second table in which
problems were nonconvex, inevitably the different algorithms found different optima. We
were also encouraged by the fact that LOQO successfully solved six of the nine problems
efficiently. However, LOQO did fail on three of the problems. The reason is that the algorithm
eventually took small steps and therefore failed to make sufficient progress before
the iteration limit was reached. No merit reductions to the steplength were incurred and
100 20 200 1.1 34 28 4
.3 293 98 14
500 100 1000 .1 .1 4304 811 85
1000 200 500 .02 .1 5247 1434 147


5. Mittelmann's results for convex quadratic programming.
Numbers shown under the solvers are solution times in seconds. For
LOQO, the following parameter settings were used: bndpush=100
honor bnds=0 pred corr=1 mufactor=0. (These are the LP defaults. In
the next release of LOQO, these values will be the defaults for QPs too.)
.3 37 * 19
500 100 1000 .1 .1 542 * 172
500 200 1000 .1 .1 462 1389 *
1000 200 2000 .02 .1 7106 * 1311


6. Mittelmann's results for nonconvex quadratic programming.
Numbers shown under the solvers are solution times in seconds. An
asterisk indicates that the solution obtained was either infeasible or sub-
optimal. For LOQO, the following parameter settings were used: bnd-
push=100 honor bnds=0 pred corr=1 mufactor=0.
so the short steps must have arisen from some variables being very close to their bounds.
Future research will address identifying and remedying this phenomenon.
7. CONCLUSIONS
The aim of our work has been to take an existing interior-point code for quadratic
programming, namely LOQO, and modify it as little as possible to develop an efficient
code for nonconvex nonlinear programming. We found that essentially only two changes
were needed:
(1) A merit function to ensure proper steplength control.
(2) Diagonal pertubation to the Hessian matrix to ensure that the search directions are
descent directions when the problem is not convex.
AN INTERIOR-POINT ALGORITHM FOR NONCONVEX NONLINEAR PROGRAMMING 19
As noted throughout the paper, the algorithm documented herein represents our first attempt
at an interior-point code for nonconvex nonlinear programming problems. We believe
that significant improvements can and will be made over time. However, we personally
were both gratified and somewhat surprised at how robust and how efficient this initial
code has proved to be.
8.

ACKNOWLEDGEMENTS

We'd like to thank David Gay for doing much of the hard work by providing the automatic
differentiation algorithms in AMPL to compute gradients and Hessians.
9.
The second author of this paper has known Olvi Mangasarian for approximately thirty
years. I believe we first met at the first of the four Madison conferences on mathematical
programming held in May, 1970. All four conferences, organized first by Olvi, Ben
Rosen and Klaus Ritter and later by Olvi, Bob Meyer and Steve Robinson, in my opinion
contributed greatly to the development of mathematical programming during a period of
considerable excitement and great progress. The excitement was perhaps greater than Olvi
might have wished at the first conference, where tear gas pervaded the campus, the National
Guard was omnipresent, and we had to move the sessions off campus to avoid the
antiwar protests. Since these early days, Olvi and I have remained close friends, meeting
regularly at conferences, where our wives have shared many happy days visiting the art
museums while we attended sessions (and often wished we were with our wives). I have
been around long enough to have taught from Olvi's book when it was new, and have followed
all of his seminal work from the theorems of the alternative through the work of his
student Shi Ping Han on SQP to his work on proximal point algorithms, complementarity,
and most recently on detecting breast cancer, to name just a few of his accomplishments.
Through all of this he has been a great friend, congenial host, and valued companion. The
first author is young, and new to nonlinear programming. This is a wonderful way for him
to make the acquaintance of one of the icons of our field.



--R

A globally convergent Lagrangian barrier algorithm for optimization with general inequality constraints and simple bounds.
A primal-dual algorithm for minimizing a non-convexfunction subject to bound and linear equality constraints
Optimization methods for truss geometry and topology design.
GAMS: A User's Guide.
An interior point algorithm for large scale nonlinear programming.
LANCELOT: a Fortran Package for Large-Scale Nonlinear Optimization (Release
On the formulation and theory of the Newton interior-point method for nonlinear programming
Nonlinear Programming: Sequential Unconstrainted Minimization Tech- niques
Nonlinear programming without a penalty function.

AMPL: A Modeling Language for Mathematical Programming.
Computational experience with a safeguarded barrier algorithm for sparse nonlinear programming.
A primal-dual interior method for nonconvex nonlinear pro- gramming
Solving reduced KKT systems in barrier methods for linear and quadratic programming.
The largest small hexagon.
Test examples for nonlinear programmingcodes.
Antenna array pattern synthesis via convex optimization.
Interior point methods for linear programming: computational state of the art.
Benchmarks for optimization software.
MINOS 5.4 user's guide.


LOQO: An interior point code for quadratic programming.
Symmetric quasi-definite matrices
Linear Programming: Foundations and Extensions.
--TR
On the formulation and theory of the Newton interior-point method for nonlinear programming
An interior point potential reduction method for constrained equations
A globally convergent Lagrangian barrier algorithm for optimization with general inequality constraints and simple bounds
Lancelot
Test Examples for Nonlinear Programming Codes

--CTR
Alicia Troncoso Lora, Advances in optimization and prediction techniques: Real-world applications: Thesis, AI Communications, v.19 n.3, p.295-297, August 2006
V. Adetola , M. Guay, Brief paper: Parameter convergence in adaptive extremum-seeking control, Automatica (Journal of IFAC), v.43 n.1, p.105-110, January, 2007
Hande Y. Benson , Robert J. Vanderbei , David F. Shanno, Interior-Point Methods for Nonconvex Nonlinear Programming: Filter Methods and Merit Functions, Computational Optimization and Applications, v.23 n.2, p.257-272, November 2002
Sasan Bakhtiari , Andr L. Tits, A Simple Primal-Dual Feasible Interior-Point Method for Nonlinear Programming with Monotone Descent, Computational Optimization and Applications, v.25 n.1-3, p.17-38
Wei-Peng Chen , Lui Sha, An energy-aware data-centric generic utility based approach in wireless sensor networks, Proceedings of the third international symposium on Information processing in sensor networks, April 26-27, 2004, Berkeley, California, USA
Hande Y. Benson , Arun Sen , David F. Shanno , Robert J. Vanderbei, Interior-Point Algorithms, Penalty Methods and Equilibrium Problems, Computational Optimization and Applications, v.34 n.2, p.155-182, June      2006
Igor Griva, Numerical Experiments with an Interior-Exterior Point Method for Nonlinear Programming, Computational Optimization and Applications, v.29 n.2, p.173-195, November 2004
Richard H. Byrd , Jorge Nocedal , Richard A. Waltz, Feasible Interior Methods Using Slacks for Nonlinear Optimization, Computational Optimization and Applications, v.26 n.1, p.35-61, October
M. D'Apuzzo , M. Marino, Parallel computational issues of an interior point method for solving large bound-constrained quadratic programming problems, Parallel Computing, v.29 n.4, p.467-483, 01 April
Hande Y. Benson , David F. Shanno, An exact primal---dual penalty method approach to warmstarting interior-point methods for linear programming, Computational Optimization and Applications, v.38 n.3, p.371-399, December  2007
Helmut Maurer , Hans D. Mittelmann, Optimization Techniques for Solving Elliptic Control Problems with Control and State Constraints: Part 1. Boundary Control, Computational Optimization and Applications, v.16 n.1, p.29-55, April 2000
Hans D. Mittelmann, Verification of Second-Order Sufficient Optimality Conditions for Semilinear Elliptic and Parabolic Control Problems, Computational Optimization and Applications, v.20 n.1, p.93-110, October 2001
L. N. Vicente, Local Convergence of the Affine-Scaling Interior-Point Algorithm for Nonlinear Programming, Computational Optimization and Applications, v.17 n.1, p.23-35, Oct. 2000
Silvia Bonettini , Valeria Ruggiero, Some iterative methods for the solution of a symmetric indefinite KKT system, Computational Optimization and Applications, v.38 n.1, p.3-25, September 2007
Stefania Bellavia , Benedetta Morini, Global convergence enhancement of classical linesearch interior point methods for MCPs, Journal of Computational and Applied Mathematics, v.151 n.1, p.171-199, 1 February
P.-A. Absil , Andr L. Tits, Newton-KKT interior-point methods for indefinite quadratic programming, Computational Optimization and Applications, v.36 n.1, p.5-41, January   2007
Helmut Maurer , Hans D. Mittelmann, Optimization Techniques for Solving Elliptic Control Problems with Control and State Constraints. Part 2: Distributed Control, Computational Optimization and Applications, v.18 n.2, p.141-160, Feb. 1, 2001
Roman A. Polyak, Nonlinear Rescaling as Interior Quadratic Prox Method in Convex Optimization, Computational Optimization and Applications, v.35 n.3, p.347-373, November  2006
Nicholas I. M. Gould , Dominique Orban , Philippe L. Toint, GALAHAD, a library of thread-safe Fortran 90 packages for large-scale nonlinear optimization, ACM Transactions on Mathematical Software (TOMS), v.29 n.4, p.353-372, December
Francisco Facchinei , Giampaolo Liuzzi , Stefano Lucidi, A Truncated Newton Method for the Solution of Large-Scale Inequality Constrained Minimization Problems, Computational Optimization and Applications, v.25 n.1-3, p.85-122
Silvia Bonettini , Emanuele Galligani , Valeria Ruggiero, Inner solvers for interior point methods for large scale nonlinear programming, Computational Optimization and Applications, v.37 n.1, p.1-34, May       2007
A. Migdalas , G. Toraldo , V. Kumar, Nonlinear optimization and parallel computing, Parallel Computing, v.29 n.4, p.375-391, 01 April
nonconvex optimization;nonlinear programming;interior-point methods