J. SIAM OPTIMIZATION
Vol. 5, No. 3, pp. 590-640, August 1995
{)
1995 Society for Industrial and Applied Mathematics 007
A SEQUENTIAL QUADRATIC PROGRAMMING ALGORITHM USING AN INCOMPLETE SOLUTION OF THE SUBPROBLEM*
WALTER MURRAYt AND FRANCISCO J. PRIETO:
Abstract. We analyze sequential quadratic programming (SQP) methods to solve nonlinear constrained optimization problems that are more flexible in their definition than standard SQP methods. The type of flexibility introduced is motivated by the necessity to deviate from the standard approach when solving large problems. Specifically we no longer require a minimizer of the QP subproblem to be determined or particular Lagrange multiplier estimates to be used. Our main focus is on an SQP algorithm that uses a particular augmented Lagrangian merit function. New results are derived for this algorithm under weaker conditions than previously assumed; in particular, it is not assumed that the iterates lie on a compact set.
Key words, nonlinearly constrained minimization, sequential quadratic programming, quasiNewton method, large-scale optimization
AMS subject classifications. 49D37, 65K05, 90C30
1. Introduction. The problem of interest is the following:
NP
minimize xE
s.t.
F(x) c(x) >_ O,
and c n _+ m. Since we shall not assume that second where F n _+ derivatives are known, computing x*, a point satisfying the first-order Karush-KuhnTucker (KKT) conditions for NP is the best that can be achieved. Such points are feasible and satisfy the following conditions"
(I.i)
VF(x*) Vc(x*)TA*,
)icj(x*)
,
0 j
1,...,m
for some nonnegative multiplier vector A* E m. Whenever the term "KKT point" is used in the following sections, it will mean a point satisfying the first-order KKT conditions for NP. Despite this theoretical limitation, we prefer some KKT points to others to try and satisfy our real purpose of finding a minimizer. For example, if the initial estimate is feasible we do not wish to converge to a nearby KKT point if at that point the objective function is higher. We use the term stationary point to denote a point that is feasible and satisfies (1.1) for some multiplier vector A E m that is not necessarily nonnegative. Typically SQP algorithms generate a sequence of points {xk} converging to a solution, by solving at each point, xk, a quadratic program (QP) of the form
Received by the editors January 1, 1990; accepted for publication (in revised form) March 28, 1994. This research was supported by National Science Foundation grant DDMo9204208, Department of Energy grant DE-FG03-92ER25117, Office of Naval Research grant N00014-90-J-1242, and the Bank of Spain. Systems Optimization Laboratory, Department of Operations Research, Stanford University, Stanford, California 94305-4022 (walter@sol-walter. stanford, edu). Department of Statistics and Econometrics, Universidad Carlos III de Madrid.
590
A SEQUENTIAL QUADRATIC PROGRAMMING ALGORITHM
591
QP
minimize
pE
Soto
VF(xk)Tp + 1/2pTHkp C(Xk) + VC(Xk)p >_ O,
for some positive definite matrix Hk. Let Pk (referred to as the search direction) denote the unique solution to QP. We define Xk+l xk / OkPk, where the steplength ak is chosen to achieve a reduction in a merit function. SQP algorithms are viewed by many as the best approach to the solution of NP when n is small < 200 ). As the size of the problem grows, usually so does the relative importance of the effort to solve QP when compared to the total effort. Indeed, for many large problems the effort to solve QP dominates the total effort. When the minimizer of QP is used to define the search direction, it is not necessary in any theoretical discussion of an SQP algorithm to define how the QP subproblem is solved. Most implementations of SQP methods currently available use an activeset method to solve the QP subproblem. For a comprehensive survey of active-set methods see [18], [13], and [17]. The potential number of iterations to solve a QP using an active-set method grows exponentially with n. In practice the number of iterations grows much more slowly than exponential (if this was not the case active-set methods would be hopelessly inefficient). Nonetheless, the number of iterations required to solve a large QP is usually large. In any implementation of an SQP method it is necessary to limit the number of iterations allowed to solve a given QP subproblem. If the QP solution process is terminated prematurely the SQP algorithm may break down. It is in part for this reason that the development of SQP methods for largescale problems has been inhibited. Even for small problems there are occasions when the number of QP iterations is excessive. Since the definition of "small" continues to increase as computers become more powerful we can expect the cost of solving the subproblems to grow in importance. In the algorithms presented here we have endeavored to improve the efficiency of SQP methods by circumventing the need to determine the minimizer of QP. We show that a suitable search direction may be computed from information available at any stationary point of QP. Stationary points occur as iterates within most active-set methods to solve QP and for such methods the number of iterations to determine a stationary point increases only linearly with the size of the problem. Consequently, the search direction may be found by applying an active-set method to QP and terminating the procedure early. It may be thought that by expending much less effort to compute the search direction, the number of iterations for the outer algorithm may increase. However, it has been observed that large numbers of QP iterations are required only when xk is a poor approximation to x*, that is, when the QP subproblem does not model the nonlinear problem well. We hypothesize that a search direction based on the minimizer of such subproblems is little better than using information at a stationary point. Our preliminary results reported in 6 support this hypothesis. Not solving the QP subproblem also implies that we do not know the QP multipliers, which are often used to estimate the multipliers of NP. In general, SQP methods usually use some specific estimate of the NP multipliers in the definition of the method and hence in the proof of convergence. When solving large problems specific definitions of multiplier estimates are not always computationally attractive. In our analysis we allow for flexibility in how multipliers are defined by requiring only that the multiplier estimates satisfy certain conditions.
=--
592
WALTER MURRAY AND FRANCISCO J. PRIETO
1.1. Incomplete solutions for QP subproblems. There have been other proposals to define the search direction for an SQP algorithm other than as the minimizer of the QP subproblem. In Dembo and Tulowitzki [9] an algorithm is analyzed for which the search direction Pk has the property that
where denotes the minimizer for the kth QP subproblem, (unless stated otherwise all norms in the paper are g2-norms). We follow a different approach and define a search direction for which the effort to compute it has a guaranteed bound. A different algorithm, but using the same approach, was suggested by Gurwitz and Overton [20]. However, no global convergence results were given for their algorithm. In the course of solving a QP an active-set method generates iterates that are stationary points. We show that such points may be used to construct a suitable search direction. The step to the stationary point is not generally an adequate search direction. However, if the stationary point is not a minimizer then there exist nonoptimal multipliers. We show how an auxiliary direction may be constructed using information about the nonoptimal multipliers. This auxiliary direction, when combined with the step to the stationary point, gives a suitable search direction. Terminating the QP algorithm prior to obtaining a solution impacts the SQP algorithm in a number of critical ways. Not only is the search direction different, but also the QP multipliers will not be available. The merit function of principal interest requires the definition of a search direction in the space of the multipliers. In the past, this search direction has been defined using the QP multipliers. The fact that such multipliers are positive was crucial in the analysis of these algorithms. The consequences of terminating the QP solution process early are therefore far reaching. The remainder of this paper is organized as follows. Section 2 describes the form of the general algorithm, and the definition of the search direction. Section 3 studies the convergence properties of the algorithm; it is shown that such an algorithm is globally convergent. In 4 we show that the algorithm converges superlinearly. We also show that the penalty parameter used in the merit function is bounded. Section 5 considers the use of alternative merit functions. Finally, 6 presents numerical results obtained from an implementation that uses the merit function of principal interest.
2. Description of the algorithm. The search direction we propose could be used with most of the merit functions analyzed in the literature. However, our primary interest is the following merit function:
p
(2.1)
LA(x, A, s, p)
F(x) AT(c(x) s)+ 1/2P(C(X) s)T(c(x) s),
where s _> 0 are slack variables, and the scalar p is known as the penalty parameter. This merit function was suggested by Gill et al. [16] and is used in the SQP code NPSOL. It is similar to merit functions proposed by Wright [34] and Schittkowski [32]. Although our primary interest is this specific merit function, we also show (5) how the ideas discussed can be extended to the use of other merit functions. The reason for our focus on this merit function is due to the success in practice of NPSOL. The merit function is also used in a new SQP code, LSSQP [10], designed to solve
large problems.
A SEQUENTIAL QUADRATIC PROGRAMMING ALGORITHM
593
The search is performed on an expanded space, including the Lagrange multiplier A, and the slack variables s. The symbols p, and q will be used to denote the components of the search direction on the corresponding subspaces. In this case, the value of the merit function as a function of the steplength will be denoted by
estimates
,
((; x, p, :k, s, q, p)
,
LA (x + p, )t + c, s + cq, p).
The explicit reference to the parameters will be omitted in what follows. The derivative of with respect to c is denoted by . Also, Ck(C) and (c) will be used to indicate the values of and evaluated at (Xk, Pk, Ak, k, Sk, qk, Pk). The following conventions will be used in the rest of the paper:
gk
VF(xk),
Ak =-- Vc(xk),
Ck
=--C(Xk),
and the symbols and k will be used with the same meaning as Ak and ck, but restricted to the set of active constraints at the given point. The term active constraint will be used to designate a constraint that is satisfied exactly at the current point (cj(x) 0 in NP, or -cj in QP), and the set of all constraints active at a given point will be referred to as the active set at the point.
k
ap
The objective function for the QP subproblem will be denoted by
Ck(P),
(2.3)
Sometimes,
Ck(P)
g[P + PTHkp.
will denote the function of one variable
For any vector v, the notation v- will be used to denote the vector whose jth element is defined as
v-
min(0, vj).
Also, the symbol e denotes the vector (1,..., 1) T, and symbols of the form abc denote fixed scalars related to properties of the problem, or the implementation of the algorithm, where "abc" identifies the specific scalar represented.
Finally, throughout the paper we will use the symbol ]lull to denote the g2-norm of the vector u, unless we explicitly indicate that a different norm is being considered.
2.1. The algorithm. We first present an outline of the algorithm. Given positive definite, x0 and 0, select fig > I1011 and tip > 0.
H0
594
WALTER MURRAY AND FRANCISCO J. PRIETO
ALGORITHM ETSQP
k+-0 repeat Obtain the search direction Pk from the QP subproblem minp Ck(P) g[P + 1/2pTHkp s.t. Akp + ck >_ 0 Compute #k, an estimate of A* such that
if Pk-
0
sk from
Compute
else
(sk)j
(Sk)j
max(0,
Compute
end if
qk +--
sk from
max(0, (ck)j
T if Ck(0) <_ -pkH}p} +- Pk-1 Pk else
AkPk -+-Ck- Sk
Pk end if if Ck (1)
( max.2p_,
_< Ck (0) +
+
liCk Skll 2
/
a (0)
--ce or Ck(&) > Ck(0)+ a&(0) do
else
Select & e
end if while c(xk
end do
ck +---
(0,1) to satisfy (a) < (0)+.a(0), I(a)l <-(0)
+ Pk) a +-- a12
&
Compute 9k+1,
Update Hk to form
Ak+l and c+ Hk+l
k+-k+l
until convergence
The following are some comments on the steps of the algorithm. (i) At each point xk, we form the QP subproblem
(2.4a) (2.4b)
(2.5a) (2.5b)
for some vector rk E
minimize
pE }
gp + 5 lpT
--ck,
subject to Akp >_
and determine a stationary point for QP, that is, a point ihk satisfying
gk
+ Hkk ATkrk, r[ (Akk + Ck) O, Akk + Ck >_ O,
m (the QP multipliers at i5}).
A SEQUENTIAL QUADRATIC PROGRAMMING ALGORITHM
595
From information available at the stationary point we construct a search direction Pk and #k an estimate of )*. The precise conditions that Pk and #k need to satisfy are given later in this section. If Pk -0, we set Ak #k and terminate. Otherwise,
we compute the search direction in the space of the multiplier estimates
k as
(2.6)
(ii) The slack variables sk
k
#k
Ak.
are computed from
max (0, (ck)j)
(Sk)i
max
O, (ck)j
(,kk)j
Pk-1
otherwise.
These values minimize the merit function (2.1) at (Xk,k,Pk-1) with respect to the slack variables. The slack variables sk appear in the merit function (2.1) as part of the term ck Sk. From (2.7), this term takes the value
min (0,
(2.8)
(ck)j -(sk)j
(Ck))) min((ck)j, p-_
if Pk]
0,
otherwise.
The following inequality will be useful in the analysis of the algorithm:
To simplify the notation in the justification of this result, we drop the subscript k. If cj sj cj then clearly Ic sjl-Icl _> 0 <_ Icj sjl. Otherwise, cj sj t cj and If cj sj ?t cj and cj >_ O, then < 0. From (2.8) we get cj- sj < cj < 0, and hence Icj -sjl > Icjl >_ [c- I. We have cj shown Ic-l <_ Icj -sjl under all circumstances, implying (2.9).
c
Ic-I.
(iii) The search direction in the space of the slack variables of slack variables for the QP subproblem, i.e.,
qk
qk
is set to the vector
Akpk -t--ck
Sk.
For a linear constraint this choice keeps the corresponding slack variable at its optimum value.
(iv) The penalty parameter will not be modified if the condition
is satisfied, where
Ck(c)
is defined in
(2.2). Otherwise, we define the penalty param-
eter as
(2.12)
p
max(2p_,,/k, p),
where fp is some positive constant,
596
WALTER MURRAY AND FRANCISCO J. PRIETO
P
and
k(Pk) "+- (2Ak #k)T (Ck St:)
I1 11
Ok was defined in (2.3). It will be shown that the definition (2.12) ensures that (Pk,k, qk) is a sufficient descent direction for the merit function, in the sense that condition (2.11) holds for this value of the penalty parameter. (v) The steplength ak > 0 is computed to reduce Ck(a) while keeping the constraint violation bounded. The termination conditions for the linesearch are as follows:
If
() ,(0) < I(0),
set &
1. Otherwise, find an & E
(0, 1) such that
(2.15a) (2.15b)
where 0
(a) (0) < ai(0), I() > I(0),
<a<r<
1.
If the condition
(2.16)
c(xk + &pk) >_
-ce
holds, we define ak &; otherwise we compute ak by performing a backtracking linesearch from & until (2.15aa) and (2.16) are both satisfied. It will be shown later that such a steplength always exists, and that Algorithm ETSQP is well defined. This
definition of the steplength ensures that c(xk) >_ algorithm could be used to determine ak when anticipate such events will be rare.
--ce for all k. A more sophisticated
(2.16)
does not hold. However, we
(vi) Finally,
xk and
Ak
are updated from
(217)
2.2. The definition of the search direction. At each iteration of Algorithm ETSQP an inner iteration is performed to determine the search direction by solving the QP subproblem (2.4a) using an active-set method. The following is an outline of
a suitable algorithm to determine the search direction. The outer iteration subscript has been omitted, and the subscript refers to the inner iterations.
We assume that positive constants/p,/b,
/M have been defined
(/b _< 1).
A SEQUENTIAL QUADRATIC PROGRAMMING ALGORITHM
597
ALGORITHM SD
Compute p0 satisfying:
Apo + c > 0,
Ilpo
p0,
llc-II,
BTRo
llc-II
A corresponding
Form 0, the active-set matrix at to active QP constraints at p0
i.--0
as the set of all rows in
repeat
Compute/3i from
7i
6--
fi.i
min(1 infj{ -cj
Pi
( )0 +ap
aT-
Pi+l
Set
i+l to be the active-set matrix at
p
+
ii+l
until p is a stationary point. if>0 pp
else
Set vr 1 if #r < b mini #j, otherwise set vr Compute by solving: min{r .i v}
d
"),
d/[[d[[
.--
-
*- 0
+ ai5 min(- (g + H)Td inf,{ c a"d a"d < O} M ) dTHd if [[i5 + q/d[[ > 11i5[I
p+--+Td
p+-p
else end if end if
Some comments on this procedure are presented below. (i) An initial feasible point p0 of the QP subproblem is obtained. When the minimizer of the QP is used as the search direction, then, given the uniqueness of p,
the choice of P0 is irrelevant. If we determine the search direction from a stationary point that is not a minimizer, the sequence of stationary points that we compute depends directly on the value of p0. We wish to define the initial point in such a manner that all stationary points are satisfactory points at which to terminate the solution process. It will be seen that the following conditions on p0 are sufficient to ensure our objective. For some constant p > 0,
(2.18)
Ilpoll
pllc-II
and
gTpo
(ii) A sequence of feasible descent steps are taken, for example, by first computing the unique step i to the minimizer of the QP on the current working set as the least-
598
WALTER MURRAY AND FRANCISCO J. PRIETO
length solution of the system of equations
(2.19)
Ai
0
-#
0
a8
where Pi i8 the current estimate. A step /i i8 taken, where i is obtained one or the step to the nearest constraint,
either
(2.20)
7i
min(1 ipf -cj + aypi
aT_
The QP algorithm may be terminated at any stationary point (Algorithm SD is terminated at the first stationary point.) It will be seen in the proof8 that to alway8 use as the search direction will not in general ensure convergence. (iii) If i8 the minimizer of the QP 8ubproblem the search direction p i8 defined else p
.
,
(.1)
P
+
f 1111 < I1 + dll, otherwise,
where the vector d and the scalar are computed with the following properties: u8eu 0, d I111 X. The rate of descent along d is "suciently" large. By this we mean d satisfies
A
1,
where 0
<d
H + g and d* solve8
min
s.t.
(.a)
Ad o,
There are many procedure8 for computing a suitable vector d. For example, if the are bounded above and below and i ha8 full row rank then a singular values of follows. Define a vector v to 8atis suitable d may be computed
vi
1 0
if i < 0, otherwise.
v and define d
We then compute
the least-length solution of iy
For this direction d we have
(.24)
Td T Tv
1
n
1.
.
it follows I]]] is bounded. We shall now show Under the sumptions made on d is a "sucient" descent direction. Let u* denote the 801ution of the problem
minu
s.t.
0,
Tu IIll
A SEQUENTIAL QUADRATIC PROGRAMMING ALGORITHM
599
We have #Tu* > m mini #j. Define z
(JTd*
since z2
Aid*.
We get
3
T d*
77T( > mllll m.in
> 0. If inequality we get
liAII is bounded it follows IIll is bounded. From .(2.24) and the above
I I ffTd 7?T"id 7?Tv --< ildll rn?j< mll ll Ildll ffTd,
Lemma 2.1 presents some properties of the solutions of (2.23). These properties are based on the observation that the cost vector and the coefficients of each constraint can be normalized without affecting the feasible region or the solutions of the problem. Since we are concerned with sequences we reintroduce the outer subscript. Define k k/llkll and a matrix Bk whose jth row is the normalized jth row of i. The
problem
-
mind
s.t.
(2.25)
AT gkd Bkd > 0,
has the same feasible region and the same solutions as (2.23). We tacitly assume no row of .zii is a zero vector, otherwise it could be omitted from both problems. 0 it implies/hk is the minimizer of the QP. Likewise, if IIk LEMMA 2.1. Given a subsequence of iterates (Xk}, generated by Algorithm ETSQP and such that for all of them k Pk, the directions dk obtained as solutions of -Tz. * 0 along 1. Furthermore, if gkdk (2.23) at each point satisfy y t k < 0 and the subsequence, then either k --* 0 or for any limit (t, B) of the sequence { (k, Bk)}, BT, with > O. defined as in (2.25), it holds that Since Pk --Pk is a feasible descent direction of (2.23) at d 0 it follows Proof. * * that d 0 is not optimal, and the solutions of (2.23) satisfy gkdk < 0 and I]dk 1. Consider now the sequence of problems of the form (2.25) and the problem obtained from a limit of the sequence { (k, Bk) }. The feasible regions of all problems are compact convex polytopes; if we denote the vertices of the polytope corresponding to problem k by {d }, where the index takes a finite number of different values, it holds --. dl, a vertex for the polytope corresponding to the feasible region that for each l, of the limit problem (assume without loss of generality that the convergent subsequence has been chosen so that the number of vertices is the same for all problems in
Ildk I1
-
IIc
d
the subsequence). Any feasible point of the limit problem, d, can be written as a convex combination of the vertices dl, d t dr" We can then construct for any feasible d a sequence {dk}, where each point dk is defined as dk =_ ldk, having the properties that dk is feasible for the kth problem (2.25), and dk d. _-T. ^T * Ifk 0 then g dk ---* 0 implies gk dk 0 and it must hold that d 0 is an optimal solution of the limit problem, implying that there exists a vector > 0 FI BTI. satisfying Note that gk dk O, if and only if kdk --+ O, where dk is a sufficient descent direction. The scalar is given by
Yt
t
(2.26)
min(9,-, 7M),
600
WALTER MURRAY AND FRANCISCO J. PRIETO
where "7M is a specified upper bound on the steplength,
(2.27)
is the largest feasible step
"
in.f{
3
from/5 along d, and
(2.28)
is the step to the minimizer of (/5 / d).
2.3. The multiplier estimates. Equation (2.6) defining the search direction on the ,multiplier space k requires the computation of an estimate #k for the Lagrange multipliers. The estimates {#k } are then used to update { Ak }, the Lagrange multiplier estimate used in the merit function. To allow flexibility in algorithm design we have chosen to specify conditions on the multipliers estimates k rather than give explicit definitions. It will be shown that the following conditions on k are sufficient to ensure that the algorithm is globally convergent. MC1. The estimates #k are uniformly bounded in norm, that is ]]#k]] < MC2. The complementarity condition #(Apk + ck) 0 is satisfied at all iterations.
-
cj
+ @15 ayd
(a + S )rd
drHd
}
We may satisfy these conditions by choosing #k 0. Condition MC2 is made for convenience; condition MC1 and the form in which the multiplier estimates are updated imply that {Ak} are uniformly bounded. LEMMA 2.2. If condition MC1 holds, then ] for all k. The proof is by induction. We select to satisfy A0 Proof. fl. om (2.17),
(2.29)
"k-bl
"k -[-Olk(#k "k), k O.
1, we have
Using norm inequalities and 0 < ak
_ _
.
as required.
2.4. Second-order information. We choose the matrices {Hk} to be positive definite and bounded, with bounded condition number. In practice, such matrices may be generated (see [15]) by updating a quasi-Newton approximation to the Hessian of the Lagrangian function or the Hessian of the augmented Lagrangian function in each iteration together with certain safeguards (for example, if the factors of Hk are updated, by enforcing bounds on the size of the elements, and ensuring sufficiently positive diagonal elements). These conditions can be written as follows: HC1. H < oc is the largest eigenvalue of {H}. HC2. .H > 0 is the smallest eigenvalue of {Hk}. 3. Global convergence results. The results in this section establish global convergence properties for Algorithm ETSQP, under certain assumptions on the problem NP. We first introduce these assumptions, and then, under the condition that they hold, we prove the following results: (i) The iterates {x} lie on a compact set.
A SEQUENTIAL QUADRATIC PROGRAMMING ALGORITHM
601
In Lemma 3.1 we show that the quantities associated with the algorithm
are well defined at all points. In Lemma 3.2 it is shown that if
IlXkll is large then IlPkll cannot be arbitrarily small. In Lemma 3.3 we show that p computed using Algorithm SD satisfies
(p)
where/1 and/2
gTp + pT Hp <_ -lpTHp +/211c- sll
are positive constants.
Lemma 3.4 proves that the sequence {xk} lies on a compact set. Lemma 3.5 shows that the sequence {p} also remains bounded.
(ii) The sequence {llPkll} dominates the sequence {llxk-x*ll}, where x* denotes a
0 KKT point closest to xk. The main implication of this result is that IIP is sufficient to ensure that xk --* x*, a KKT point of NP. It is shown in Lemma 3.6 that the KKT points for problem NP are
isolated.
0 along a subsequence then along the same subsequence 117rk -* 0. Lemma 3.8 introduces another preliminary result, proving that if p 0 along a subsequence then along this subsequence Ilxk x* --* 0, where x* is a KKT point for NP nearest to x. Moreover, for large enough k, Pk is the minimizer of the QP subproblem, and the correct active set at x* is identified. is given in Lemma 3.9. The proof that IlPkll dominates Ilxk (iii) Bounds on the growth of the penalty parameter Pk. We cannot prove that Pk will remain bounded in the algorithm without stronger conditions on the multiplier estimate #k, but we can show that its growth is bounded by certain quantities related with the algorithm, and that is enough to prove convergence. We show in Lemma 3.10 that at all the iterations where the penalty parameter is modified the following bounds hold,
Lemma 3.7 shows that if Ilxk --x*
*
-
x*ll
In Lemmas 3.11 and 3.12 we show that similar inequalities hold at all
iterations.
(iv) The steplength ak is bounded away from zero if we are not close to a solution.
We first need a bound on the second derivatives of (c). In Lemma 3.13 we prove that (ak) _< N for some positive constant N. In Lemma 3.14 we show that, if IlPk is large enough, there exists a value
> 0 independent of the iteration such that ak >_ (. (v) In Theorem 3.15 we show that xk x*. (vi) Finally, we prove that Ak A*.
This result requires stronger conditions on the multiplier estimate #k than just MC1 and MC2. We start by introducing a third condition
-
MC3. Lemma 3.16 strengthens the result in Lemma 3.14 showing that, under the new conditions on the multipliers, ck is uniformly bounded away from zero. In Theorem 3.17 we show that )k A*.
602
WALTER MURRAY AND FRANCISCO J. PRIETO
3.1. Assumptions. Some of the following assumptions make use of the concepts of stationary points and KKT points at infinity. We will say that NP has a stationary point at infinity if there exist sequences (xk} and (k} such that Ilxkll oc and/or IIrkll-- c, and
ck-
As before let Ba denote a matrix whose rows are the normalized rows of Ak and k T denote the normalized gradient vector. Define k so that Ak rlk--gk (B[Pk--k)llgkll. If in addition to the preceding conditions we have p >_ 0, where indicates the limit of some subsequence (k}, we then say there is a KKT point at infinity. Finally, we will say that strict complementarity does not hold at some stationary point at infinity if for the preceding sequences and some constraint j we have
-
0,
AkTrlk
gk --+
0,
rlck
O.
(ck)j
We make the following assumptions. A1. For some constant > 0, the global minimum of the problem
c
- 0 and ()j
0.
minimize
F(x) c(x) >_-13e,
s.t.
is bounded below. A2. There exist no
KKT points at infinity for problem NP. A3. F, cj, and their first and second derivatives are continuous and uniformly bounded in norm on a compact set. A4. The Jacobian corresponding to the active constraints at all KKT points has full rank. A5. A feasible point Pko exists to all the QP subproblems, satisfying
for some constant/p > 0. A6. Strict complementarity holds at all stationary points of NP, including stationary points at infinity, if they exist. A7. The reduced Hessian of the Lagrangian function is nonsingular at all KKT points. The larger the value of 3c, the stronger is assumption A1. There will be problems, for example F(z) f(x)Tf(x), where it is known a priori that Assumption A1 holds with/ oc. Also, if A1 does not hold with 3c 0 then it is possible for any reasonable algorithm to diverge. Assumption A5 imposes conditions, on the initial point for the QP. It is possible that no point satisfies these conditions; this would be the case for example if one of the QP subproblems generated by the algorithm is not feasible. Nevertheless, by introducing an additional variable it is possible to construct a modified problem for which satisfying the conditions on Pko is trivial. Consider the problem
minimize
(3.1)
s.t.
(x, c(x)+e>_O and >_0,
where E and w E [0, 1]. The KKT points for this problem are also KKT points for NP if NP is feasible and w is sufficiently close to one. The modified problem is
A SEQUENTIAL QUADRATIC PROGRAMMING ALGORITHM
always feasible, and the corresponding QP subproblem takes the form
minimize
(p,)en+
(1 -)[++ 1/2
ck
Soto
+ Akp + ke + e
2+p>_O.
_
603
)H
0,
For this QP subproblem the point
(
VkTp0
T (1--W)gk
o
is feasible since we can ensure that 2k _> 0. Therefore there always exists a feasible point that satisfies A5 with p 1 since IIP0[I [[(c + "2ke)-II and
implying that the conditions on Pko in Assumption A5 are trivial to satisfy for (3.1)o 3.2. Existence of the iterates. We start by showing that all the quantities associated with the algorithm are well defined. In particular, we show that the choice of penalty parameter ensures (2.11) is satisfied and that the steplength exists. LEMMA 3.1. Under Assumptions A3, A5 and conditions HC1, HC2, the procedures given in the algorithm to compute the values of the penalty parameter Dk and the steplength ak are well defined. Proof. We drop the subscript k denoting the iteration number, to simplify the
notation.
x,
,
Consider the gradient of the merit function L A, defined in and s,
(2.1),
with respect to
VLA (x, s)
,
g(x)
A(x)TA + pA(x)T(c(x) s)
-(c(x) s)
p((x)
It follows from (2.6), (2.10), and (2.2) that (0) is given by
(3.3)
where g, A, and c are evaluated at x. If IIc- sll 0, from (2.9) and (2.18) we have P0 1/2pTHp _< (Po) 0 it follows that
(O) pT g <_ pT Hp,
implying that p does not need to be modified. If IIc- sll > 0, we obtain from (3.3) that for p- t (defined in
_
)
)
pTg +
0, and since (p)
(2.13))
>o
(o)
gTp + (2A )T(c- ) ,llc- *II
(2.11)
--lpTrr/-/P
which implies the desired descent condition
is satisfied for all p
604
WALTER MURRAY AND FRANCISCO J. PRIETO
An immediate consequence of (2.11) and the properties of bound on the directional derivative:
(.)
(2.12)) that
times.
(0) <-1/2.1111
It follows from the procedure to increase the value of the penalty parameter (see ,ok --o oo if and only if the parameter is increased an infinite number of
.
Hk
is the following
We also need to prove that the value of k introduced in the algorithm is well defined. We show that if condition (2.14) is not satisfied, a steplength &k (0, 1) that satisfies conditions (2.15) always exists (see, for example, Mor and Sorensen [23]). Define the functions
() _= () (0) (0), () () (0), and note that from a < r/and (0) < 0, implied by (2.11), we have X(a) (a) he(0) < (a) (0) (3.5)
for any a. If (2.14) does not hold,
(1)- (0) > he(0) X(1) > O, and we also have x(O) O. om hese wo results and he mean-value heorem, here will be a point & e [0, 1] such that X() > 0, and from (3.5), (5) > 0. From b(0) < 0 we have (0) < 0, and the continuity of (Assumption A3) will imply the existence of a zero of in (0, ). Let & denote the smallest point in (0, &) such that ((c) 0, that is, (0), (c) (3.6) and condition (2.15b) is satisfied at &. From (0) < 0 we must have (a) < 0 Va e [0,&) : (a) < re(0) Va e [0, a), (3.7) implying that condition (2.15b) is not satisfied for any point in [0, Finally, from (3.5) and (3.7), we have x() < 0 w e [0,), and this together with X(0) 0 implies X(8) < 0, that is, ()- (0) < (0), (3.s) showing that & satisfies both conditions (2.15) simultaneously. We still need to consider condition (2.16). For the function h(a) c(x +ap)+e we have from (2.4b) h(O) Ap >_ -c. If-1/2fl _> cj >_ -/3, we have hi(O) >_ 0 and hi(O) >_ 1/2c > 0; if cj >_ -1/23c then hi(O) 1/2c > 0 and in any case there exists a value > 0 such that hi(a) >_ 0 (implying cj(x + ap) >_ -/3c) for all j and all a [0,&], implying that for [0, min(&, &)] both conditions (2.15a) nd (2.16) hold simultaneously.
=
_
This lemma implies that all the quantities associated with the algorithm are well defined.
A SEQUENTIAL QUADRATIC PROGRAMMING ALGORITHM
605
3.3. Boundedness of the iterates. To prove global convergence we show first that if Assumptions A1 and A2 hold, all points in the sequence (xk} generated by the algorithm lie on a compact set. We start by showing that for Ilxkll large enough we cannot have IlPk arbitrary small. LEMMA 3.2. Under Assumptions A2 and A6 and condition HC1, there exist positive constants M and e such that IlXkll >_ M IlPkll e. Assume this result does not hold. Then, for any.positive constants M and Proof. e we can find iterates such that Ilxkll >_ M and IlPkll < e, and we could construct a sequence (xk}, and its associated sequence (Pk}, along which Ilxkll --+ oc and --+ 0. IlPkll O. For this sequence, from IlPkll --+ 0 and (2.4b), we must have and from (2.5a) Also, from the definition of Pk, (2.21), it must hold that Ilihkll 0, and MC1, we must have
=
>--
-
IIc-II
Since IIpll --- 0 and [liSkll 0, using (2.21) and IIdll 1, we also have either k 0 or -k 0 for k large enough. It then follows from (2.26) that either min(k, "k) --+ 0 or k k 0 for k large enough. If k --+ 0 along a subsequence, then (2.27) implies for some constraint j that (rk)i 0 and c(xk) 0, but this would contradict Assumption A6. If ;),k --+ 0 along a subsequence, then from (2.28) and Lemma 2.1 >_ 0 in the limit, where is now defined as a limit point of (k}, where we get
The properties of this sequence,
together with ihk --4 0 and u >_ 0, imply that there exists a KKT point at infinity, Q which violates Assumption A2, so the lemma must hold. Another result we need for the compactness proof is a bound on the value of the QP objective function at the incomplete solution for the QP. LEMMA 3.3. Under Assumption A5 and conditions HC1, HC2, for p computed by Algorithm SD there exist constants 1 > 0 and/2 > 0 such that
(p)
p0,
=_
gTp d- pT H_ < --lpT Hp d- 2 [Ic [[ p_
Proof. The result will be shown by considering first the initial point for the QP, and then the descent achieved in each QP iteration. By definition
-
--Po Hpo + gTpo + pTo Hpo. Since lip011 / pllc-ll and gTpo <_ pllc-ll (Assumption Ah), condition HC1
implies
on
H
(po)
T - po Hpo + pllc-II +
H p211c- I 2
< 0 and c > 0; then for all
Consider the quadratic function b7 + 1/2c72, where b [0,-b/c] (between 0 and the minimizer), we have
(3.10)
7<_--
b
c
=
-y(b+c/)<_0
=
b + 1/2 C"2 < --C72
606
WALTER MURRAY AND FRANCISCO J. PRIETO
The change in the QP objective function at any intermediate QP iteration i can be represented as
(.1)
(p+,)
(p,)
1/2"dHd + ( + Hp)d,
where d is used to denote the QP step obtained from (2.19) or the final step d defined in (2.22), and 7 is a feasible steplength bounded by the steplength to the minimizer along v, as defined in (2.20) or (2.26). We have dTHd > 0 (from condition HC2) and (g + Hp)Td < 0 (from (2.22)), implying that we can apply the bound (3.10) to (3.11) to obtain
<_
If we have taken N iterations to compute p (the search direction), by adding the inequalities (3.12) for i= 0,..., N and using (3.9) we obtain
N
(P)
(3.13)
<
)(Po) "- E(2(Pi) --)(Pi--1))
i--1
-1
Hpo +
E 72dTi Hd
i=1
+ PllC-II + fltHfl2P IIc-112"
We can use the convexity of the function pTHp (implied by property HC2) to
write
PHP+E 7dT Hdi >- N + 1
i=1
o
+
E 7idi
i=1
H
o
+
E 7idi
i=1
N+1
pTHp"
This result together with
(3.13)
implies
(3.14)
(P)
-< -2N(1+ 1) pTHp +
pllc-II +
Z H p211c-
2
[3 Since c- _> flee the desired result follows from this inequality and (2.9). We can now prove the main result of this section. LEMMA 3.4. Under Assumptions A1, A2, A3, A5, and A6, and conditions MC1, HC1 and HC2, the sequence {xk} generated by the algorithm lies on a compact set. Proof. First we show the set of points at which the penalty parameter is modified lies on a compact set. If Pk remains bounded it follows from the manner the penalty parameter is nodified, (2.12), that there is only a finite set of such points. Therefore x. Consider the iterations k where the we need only study the case when pk penalty parameter is modified. From condition MC1 and the boundedness of the multiplier estimates )k (Lemma 2.2), we have
This result, together with Lemma 3.3 and the definition of the penalty parameter (2.13), gives
PkllCk 8kll 2
-Pk H "3fltt)llCk <: (ill +
T gkPk
I_T
kPk
+ (2Ak k)T(ck T kll fllPk HkPk.
A SEQUENTIAL QUADRATIC PROGRAMMING ALGORITHM
607
As
we have assumed pk -* o,
(3.16) implies [[Ck
we have
Skll
--*
O, and from (2.9) also
From Lemma 3.3 and (3.15)
(3.17a) (3.175)
[ + (e ,)r( ) T 1..T <-- --tk HkPl llPk HkPk + (/31 + 3/3)11ck sk]l.
If IIpll > 0 along an infinite subsequence, then it follows from I1 11 0 and MC2 that there exists an index K such that for all k >_ K in the subsequence,
>-
From (3.17b) we obtain the following bound on wk,
(3.18)
wk
T <_ --Pk Hkpl,
for k >_ K. From (3.17a) and the bounds (3.18) and (3.3), we have for sufficiently large k
This last inequality implies that Pk is not modified for all k >_ K, which contradicts our assumption that the penalty parameter was modified an infinite number of times.
We have shown that IlPkll "+ 0 along the subsequence at which the penalty parameter is modified. The boundedness of Ilxkll along this subsequence follows from Lemma 3.2. We now consider those points corresponding to iterations where the penalty parameter is not modified. From condition (2.16) on the linesearch and Assumption A1, we have F(xk) _> /3F > -oc for all k. Also, from Lemma 2.2 IIAkll is bounded, implying that
(3.19)
LA(Xk,.k, Sk, Pk) >_ F
max
(p, m/c)
>
--oo.
Since Ilxk[I is bounded when Pk P- and Ln(xk, Ak, sk, pk) is reduced when Pk Pk-1 it follows that L(xk, Ak, sk, Pk) is bounded. Moreover, for a sequence of iterations for which Pk is not changed the reduction in L(Xk, A,Sk,pk) is bounded. Let I denote the index at which Pk is modified and let I _< k _< K denote the iterates for which Pk remains fixed. It follows from the above reasoning that there exists N such that
(3.20)
)l
CK
k--I
Z (qk
K
K
Ck-I
N,
where to simplify the notation we have used Ck Ck(0). From the termination condition for the linesearch (2.15a), also have
K
(3.2)
1/2Z=H
k--I
(llpll <_ -(- Ck-I-1)
k--I
_
(3.4) and (3.20),
we
N.
608
WALTER MURRAY AND FRANCISCO J. PRIETO
This result implies that akllPll is bounded. Hence if IIxll is not bounded there must exist sets of iterates with indices, say st < k < rt for 1, 2,..., such that M, Ilxk > M for M large enough, limt_ rt c, and Ilxst It follows that if M is chosen so that M 3, max{llxiII } then Pk is constant in the interval st < k < rt. The existence of an index such that Ilxs < M is assured since we have Ilxxll _< M and at least one index in the interval for which [Ixk[I > M. From these assumptions and definitions it follows that
<-
rt --1
k--st
It follows from Lemma 3.2 that ]]pa > e for st
/1
_< k <_ rt. From (3.22)
we get
j-st
j--stq-1
but this contradicts (3.21), implying that the points generated by the algorithm must lie on a compact set. To complete this section, we show that the search direction computed from the QP subproblem is bounded. LEMMA 3.5. Under the assumptions of Lemma 3.4, the sequence {Pk } is bounded. Proof. We drop the subscript k in the proof. As all the steps taken in the solution of the QP subproblem are descent steps, we have from (2.3),
(P0)
p >-- (P) gTp + lpT.. 1/211H 1/2 p + H- 1/2 gll 2 1/2gTH-lg,
implying from HC2 and
Ilall
Ila + bll + Ilbll,
V/,HIlpll < IIn1/2Pll < IIH-1/2gll / IIH1/2p/ H-1/2gll < IIH-1/2gll + V/2(p0) + gTH-lg
The boundedness of Ilpll follows from this result, Lemma 3.4, conditions HC1 and [] HC2 and the bound (3.9). It is tacitly assumed in the remaining proofs that the Assumptions A1-A7 and conditions MC1, MC2, HC1, and HC2 hold.
3.4. The sequence of search directions {pk}. In this section we relate the behavior of the sequence {x x*}, where x* denotes a KKT point closest to xk, to 0 implies xk --+ x*, that of the sequence {pk}. In particular, we show that IlPkll and so it is enough to prove that IlPk 0 to establish global convergence. Although the KKT point x* introduced above may nog be unique, the assumptions made on the problem, and more specifically Assumption A7, imply that if IlXk- x* is sufficiently small then x* is unique, as the following lemma shows. This result allows us to work with a well-defined sequence {Xk --x*}, at least close to a KKT point; it will also imply that the limit point of the sequence generated by the algorithm is unique. LEMMA 3.6. The KKT points .for problem NP are isolated. Proof. Assume that the result does not hold, and let x* denote a KKT point for NP that is not isolated, that is, for any e > 0 there exists a KKT point y x*
-
-
-
A SEQUENTIAL QUADRATIC PROGRAMMING ALGORITHM
609
IIx* --YII < e. Consequently, there exists a sequence {Yk} such that Yk is a KKT point for all k, yk x* and Yk x*. For sufficiently small IIx* -Ykll the active sets at yk and x* must be the same; otherwise we would have for some constraint j that cj(x*) 0 with both cj(yk) > 0 and (A)j 0 along some subsequence, where A is the multiplier vector at Yk. From
satisfying
:
Assumptions A3 and A4 and (1.1) we have Ak A*, the multiplier vector at x*, but this would imply cj(x*) ,Xj 0, contradicting Assumption A6. Let Zk denote a basis for the null-space of 7((yk), the Jacobian of the active constraints at Yk, and Z* denote the corresponding basis at x*. Among all possible bases, Z is selected to have continuous first derivatives in a ball around x*. It follows from Assumption A4 and the fact the active set is constant that such bases exist. For any element of the sequence Yk and for x* we have from (1.1)
,
ZT VF(y)
The Taylor series expansion of
0
0
and
around
TVF(x*)
x*
gives
Z[VF(yk)
zTkVF(yt)= Z[(VF(yk)- Vc(yk)T,x*)
Z *T(VF(x*)
(3.23)
Vc(x*)TA*) + (VZ(x*)(VF(x*) Vc(x*)T,x*) -t-- z*T72L(x*, A*))(y x*) + o(llyk
(1.1)
in
where L(x,,) is the Lagrangian function of NP. Using by [[Yk- x*]l gives
(3.23), and dividing
(3.24)
then
z*Tv2L(x*,,X*)hk
satisfies
o(1),
where
5
Yk
x*
[lY
x*ll"
Let denote the subset of constraints active at x* and Yk. If e is sufficiently small
6k
(3.25) e(y )
we have from
0
x*) / o(lly
*11)
o(1).
{5}, with limit
Finally, for any convergent subsequence of the bounded sequence (3.24) and (3.25),
re(x*)5
0,
[] contradicting Assumption A7. This result, together with Assumption A2, implies that the number of KKT points lying on any compact region is finite. The distinctness and finiteness of the KKT and points implies the existence of e* > 0 such that for any two KKT points, say we have I]x < *, where x* is a KKT > 2*. It follows that if Ilxk point nearest to xk, then x* is unique. The next result presents some properties of the QP multipliers that will be useful for the analysis of the convergence and rate of convergence of the algorithm. LEMMA 3.7. Given a sequence of iterates {x} and the associated sequence of search directions {Pk } such that Xk x*, a KKT point for NP with multiplier vector and Pk O, then
x,
xl
x*ll
x
*
-
-
610
WALTER MURRAY AND FRANCISCO J. PRIETO
where rk are the QP multipliers at the stationary point
.
Furthermore,
j such that cj(x*) 61 > 0 we must have (r)j 0 for large enough k. If Pk 0 it follows from (2.21) that 0. Consequently, it follows from Assumption A3 that for k sufficiently large iSk <_ ough > 0. I1 11
if IlXk -x*ll <_ Kllkll for some constant K. Proof. We first show that for any constraint
I1 11
-
+
>_
> 0,
implying that the multiplier for this constraint is zero. Let .* and denote the corresponding J acobian matrices restricted to the active and #k denote their respective multipliers. From (1.1) and (2.5a) set at x* and let we have
k
*
implying
(3.26)
,T(* #k) g*
gk
HkPk
(.* .k)T#k.
From Assumption A4 that will also have full rank for large enough k, implying that #k is bounded in norm, and these results together with (2.21), Pk -+ 0 and HC1 yield rk ---* A*. Using Taylor series expansions in (3.26), we obtain
(3.27)
t*T(7k
where L(x, ) denotes the Lagrangian function for NP. The required result follows from (3.27), the condition we have imposed on the sequences {ihk} and {Xk- x*}, the boundedness of [[#kll, Assumptions A3 and ha and condition HC1. We now analyze the sequence of search directions {Pk }. The following result shows that as Pk ---+ 0 we get close to KKT points of NP and we only need to consider values Pk obtained as the minimizers for the corresponding subproblems. We complete this result by showing that a small value of IlPk also implies that the correct active set at x* is identified, in the sense that the active QP constraints at Pk correspond to the active NP constraints at x*. LEMMA 3.8. If along a subsequence pk --, 0 then along this subsequence IlXk O, where x* is a KKT point nearest to xk. For k large enough, x* is unique, Pk is the QP minimizer and the correct active set at x* is identified. Proof. A subsequence such that Pk --* 0 exists if and only if a subsequence exists 0 and the active set at Pk is constant. Let {r} denote the sequence of such that Pk indices for such a subsequence. From the definition (2.21) of Pr it follows immediately that Arp + c >_ O. From --, 0 and Assumption A3 it must hold that --, 0 and i5 0. p From (2.5) we have
x*ll
-
-
c-
-
(3.28)
ATrr
g,.
H,./?,.
0 and
rT(AD / Or)
O.
Since 15r
(3.29)
-
A SEQUENTIAL QUADRATIC PROGRAMMING ALGORITHM
611
0 it follows that
ATrr
g
We now show that for large enough r that Pr must have been computed as the minimizer for the QP. It follows from pr -* 0 and IIdll 1 that either there exists K such that for all r > K we have r 0 or r 0 (see (2.26)). If we assume the
latter it follows that
- 0,
rrCr
T
0 and
c- --, 0.
-
min(r, /r)
O.
(i) If -, 0 along a subsequence, then from will have for some constraint j
(2.27) along
and
this subsequence we
Vcj(x)T(r + /dr) + cj(x)
0
(r)j
0,
where (r)j 0 follows from the fact that the QP constraint j is limiting the step, and so it cannot be active at i5. These equations imply
contradicting Assumption A6. 0 along a subsequence, then from (ii) If
dT H.d. which implies from condition HC1 and Ildrll 1 that (0) (H,.,. / g.)Tdr ---, O. If the condition number of -i along the subsequence is bounded, condition (2.24) will
hold and for some constraint j we have (rr)i < 0, (r)j 0, giving
0 and Vcj(xr)Tr/c(x,.)
-
c(x,.)---O
and
(r)j
0,
(2.28),
cj(Xr)
0
and
(rr)/ 0,
again contradicting Assumption A6. Otherwise, from Lemma 2.1 in the limit we have that Vc(x*)T,* VF(x*) with >_ 0, implying that x* is a KKT point with a rank-deficient J acobian matrix for the active constraints, violating Assumption A4. 0 for r > K and this together with (3.28) implies We conclude therefore that >_ 0, which together pr is the minimizer of the QP subproblem. For r large enough with (3.29) and Assumption A3 implies 0, where x* is the nearest KKT point to Xr. For r large enough x* is unique. Finally, we prove that for r large enough the active set of the QP coincides with the active set of NP at x*. First note that for r large enough the active set of the QP must be a subset of the constraints active at x*, otherwise Pr is a step to a nonactive constraint implying IlPrll > > 0. Assume that for the subsequence we have Vcj(x)p + cj(x,.) > 0 and cj(x*) 0. From (2.5b) we must have (rr)j 0, implying from Lemma 3.7 that Aj 0, but this violates Assumption A6, and for r large enough the correct active set is known. This result shows that there is an e > 0 such that if IlPkll < e, then Pk is the solution of the QP subproblem, and the correct active set is known. 0 along a subsequence, then Xk x*. To We have just shown that if Pk show Pk 0, we need a stronger result, giving a relationship between the rates of convergence of the sequences {Xk x*} and {Pk}.
*
"
IIx x*ll
,
612
WALTEI MUI:ttY AND FRANCISCO J. PRIETO
LEMMA 3.9. If x* denotes a KKT point closest to xk, then there exists a constant M such that
Proof. If IIpkll > e for all k then the result holds trivially since IIxk and ]Ix* are both bounded. Again let {r} denote the indices of a subsequence such that Pr --* 0 and the active set at Pr is constant. From Lemma 3.8, for this subsequence we have Ilxr x* I 0. We assume for the rest of this proof that r is large enough so that x* is unique, p is the minimizer of the QP and the correct active set has been identified. Let A, and # denote the corresponding quantities restricted to the constraints has full row rank, and we in the active set. From Assumption A4 we know that also has full rank. assume that r is large enough so that Let Z denote a basis for the null space of with uniformly bounded norm and continuous first derivatives. From the optimality conditions for p, (2.5), we get
,
*
r,
(3.30)
Since
h(x)
Ar
Pr
(
e
)
x*),
AT*)
Cr
).
h(x*)
0, we have from the Taylor series expansion that
hj(x)
where Sj((O)j)
Sj((O)j)(xr
Vhj(x* + (O)j(x x*))
and 0 < (O)j <_ 1. We have therefore
(3.31)
From (3.23) we get
(ZTgr)=-S(O)(x-x*)
S(O)--( Z*TV2L(x*A*)
and Assumptions A4 and A7 imply that S(0) is nonsingular. It follows that for sufficiently large values of r, S(8) is also nonsingular. It then follows from (3.31) that for some positive constant
(3.32)
lix,.
x*lJ < MI(IIzT g,.II +
From Assumption A3, property HC1 and (3.30) it. follows that
for some positive constant M2. Since the subsequence {p}} such that p} --. 0 is composed of a finite number of subsequences for which pr 0 and the active set at Pr is constant, the required result [] follows from (3.32) and (3.33).
3.5. Bounds on the penalty parameter. The conditions we have imposed on the algorithm (and more specifically on the multiplier estimate) are not sufficient to ensure that the penalty parameter is bounded. However, bounds on Pk are related to the behavior of different quantities in the algorithm, and in particular to IlPkl[ and
-
A SEQUENTIAL QUADRATIC PROGRAMMING ALGORITHM
613
lick Sk I1" The following lemmas introduce bounds on the size of Pk in terms of these quantities. We start by presenting the results for those iterations where the penalty parameter is modified, and then we extend the results to general iterations. The notation kl is used in all that follows to indicate iterations at which the value of the penalty parameter needs to be modified. LEMMA 3.10. For any iteration kl in which the value of p is modified,
for some constant N. Proof. All quantities
dropped.
in the proof refer to iteration k, and so this subscript is
From the definition of/5, (2.13), and Lemma 3.3 we get
llc_ sl]2
gTp + 5PT Hp + (2 #)T(c- s) <_ --DlpTHp+ 21]c-- sll + (2- #)T(c-- S) <_ (/2 + ]12- #ll)l]C-- S.I],
(3.15)
and the above result we obtain
where/31 and/32
are positive constants. From the first bound in the Lemma,
(3.34)
If the penalty parameter needs to be modified, condition Pkz- and (3.3) implies
,
(2.11)
cannot hold for
(0)
It follows that
gTp + (2A- #)T(c-- s) llc- s]l 2 > --1/2pTHp.
gTp + lpT p + (2A #) T (c 8) > 0. Replacing in (3.35) the bound for gTp 4- 1/2pTHp given in Lemma 3.3 we obtain
(3.35)
(2A- #)T(c- 8) 4- 2]1c- 81] > pTHp,
which together with Lemma 2.2 implies
(3.36)
3/3 +
jl
lie-
> pTHp.
From condition HC2 we have I[pll <_ (1/fl.H)pTHp. If we multiply both sides of this inequality by/5 and use (3.36) to bound pTHp, we obtain
lsvH where the last inequality follows from (3.34). The second desired bound then follows
from
at
IIPlI2
We now extend these results to all iterations. To simplify notation, we shall use I and K to denote kl and k+ respectively. Thus, the penalty parameter is increased
-
1
SVH
2 pTHp <_ 3fl +f12 llc_ sl _< (33 + 2)
xx and XK in order to satisfy condition (2.11), and remains fixed at px for iterations
I,...,K- 1.
614
WALTER MURRAY AND FRANCISCO J. PRIETO
LEMMA 3.11. There exists
a constant
M such that for all l,
kt+l --I
(3.37)
Pk,
k--k
IlakPkll2 <c M.
Proof. For I _< k <_ K- 1, property (2.15a) imposed by the choice of ck, and the fact that the penalty parameter is not increased, imply that
Summing these inequalities for k
I to K- 1, 0 <_ a
K-1
1 together with
(3.4)
gives
(3.38)
lsvg
k:I
{kPk] 2 <__ 1- CK.
Consider the term p,(,-
). From (2.1) and (2.2),
F- T(o_ )+ ]_ This equation, together with the boundedness of p[[c- sl and p[cK -s
(3.39)
p(- ) M1
(implied by p > p and Lemma 3.10), and that of the multiplier estimates (Lemma 2.2), implies that for some M1 > 0,
.
+ p(F- F).
Consider now iterations for which []p][ e, so that Lemma 3.8 applies and the minimizer for the subproblem (for all other iterations p has been obtained Lemma 3.10 implies that p is bounded, and the result follows from Assumption A3, (3.39), and (3.38)). Expanding F and c about x, we get
(.40) (.40b) (3.41)
As p
F F,
,
(z x,), + O([lz, ), d,( x,) + O(]x, x[).
Mp[[p[ and [[x-
From Lemma 3.9 we have
]x,-
Mp[[p[I.
was obtained as the solution of the QP, condition (2.5a) must hold with mul0. This condition together with (3.40ha), (3.40ab), and (3.41) tiplier vector
implies
(.4)
Using again
F,
F
(i
c),, + O (m(,Ill
, ll)).
(2.11) cannot hold at
(2.5),
cIT
(0)
--pTA T I
_gTp_ pTHp"
Since p is increed at iteration I, we must have that condition that iteration, implying
gTp + (2A, )T(c S) p_,I[C S, 2 > --pTHp.
A SEQUENTIAL QUADRATIC PROGRAMMING ALGORITHM
615
The previous two results imply
pTrTcx < -px-p Hp + p
s
and this, together with the positive-definiteness of H (condition HC2), the boundedness of the multipliers (condition MC1 and Lemma 2.2) and Lemma 3.10, gives
(3.43)
flcTi < p(2A,- )T(c- 81)
M2,
0 we must have
for some M2 > 0. Consider now the term
cT in (3.42). om
--pCg
T
pcT
and from (2.9) we have conclude that there exists a constant
c;
Using c -s]. that p < p and Lemma 3.10, we such
M3
(3.44)
pcT, < M3.
om
Finally, consider the third term on the right-hand side of (3.42). It follows Lemma 3.10 and the relation p < p that there exists Ma and M5 such that
p[p]2 < M4
and hence for some constant
and
p]p]2 < Ms,
M6
(3.45)
Combining
p,O(max(llp, IIp l )) < M6.
(3.43), (3.44), and (3.45), we obtain the bound
p(F
Fr) < M2 + M3 + M6, which, together with (3.39) and (3.38) implies the desired result. LEMMA 3.12. There exists a constant M such that, for all k, M. (3.46) Pk ]]c Sk Proof. As in the preceding Lemma, let I k and K kt+. om Lemma 3.10, (3.46) is immediate for k I and k K.
To verify a bound for k I + 1,..., K- 1 we analyze some intermediate iterations k nd k + 1. We drop the iteration subscript; also let quantities evaluated at Xk+ be denoted with a tilde. om (2.8), p( ) min(p, ). Consider the following two cses: then (i) If p
-,
(3.47)
xk gives
p, le
-
IX I.
(ii) Assume now that p,j < -Ij I Expanding the jth constraint function around
cj
+ aafp + o(ll pll ).
Rewriting the previous expression, we obtain:
(3.48)
j
(1 a)cj + a(ap + cj) +
616
WALTER MURRAY AND FRANCISCO J. PRIETO
Adding and subtracting (1 -a)sy on the right-hand side of (3.48) gives
(3.49)
Since sy,
(x
)(aj
) + (x c)y + (ap + ) + O(llpl12).
(1 a)sy + a(ayp + cy) >_ O,
ap + cy, a and 1 -a are all nonnegative, we get
(3.49)
we obtain
and using this bound in
Since we assume pxSj < -lAy[ we have 5y y -gy 1- a < 1 in (3.50) we get the following inequality:
< 0. Using this bound and
Icj- jl + O(IIPlI2)
-y Iyl Iy
jl
-(1
pl
a)(cy
sy) + O(llPll 2)
Multiplying both sides by
gives
For a given iteration k < K- 1 and constraint j we have one of the following two
situations.
(i) For some iteration l, I < r iterations =/,..., k- 1, and use
Pl(ck)j
< k, pi(ct)y >_ -[(At)j]. If (3.47), we get
we add
(3.51) for
(sk)Jl
<- Pz[(Cl)j
(Sl)Jl W PzO
(
The boundedness of p,l(ck)y -(sk)yl then follows from Lemmas 2.2 and 3.11. (ii) For all iterations l, I < <_ k we have p(ct)y < -I(At)yl. We add (3.51) for r I to k- 1, to obtain
Pz[(Ck)j
(Sk)jl
<-- Pz[(Ct)j
(8t)J[ q- PxO
(
and now the desired result follows from Lemmas 3.10 and 3.11.
3.6. Boundedness of ak. Given the result of Lemma 3.11, all that is left to establish the global convergence of the algorithm is to show that the steplength is bounded away from zero. As a consequence of the weak assumptions imposed on the multiplier estimate #k, it is not possible to show that such a bound exists. However, it can be proved that the bound does exist if there is no subsequence along which IlPk O. This is enough to prove convergence. We first derive a bound on the norm of the second derivative along the linesearch. there exists a positive constant N such that LEMMA 3.13. For 0 <_ 0 <_
II--
,
(0) _< N.
A SEQUENTIAL QUADRATIC PROGRAMMING ALGORITHM
617
Proof. We again drop the subscript
k. From
(3.2),
V2LA
V2F- -j(Aj
From the definition of
,
p(cj si))V2cy -A -pA
+ pATA -AT
0
--pAT
I pI
I
given in
(2.2),
we get
"(0)
(3.52)
w + E((0)- (0)) p v (o)p
+ p(A(O)p- q)T(A(O)p- q) 2T(A(O)p q),
s(O)
=_ s
where the argument 0 denotes quantities evaluated at x + Op, except for and
+ Oq
w vF(0)- E( + 0)v(0).
We now derive bounds on the terms on the right-hand side of (3.52). For the first term we can write
(3.53)
pTWp < NIlIP
M1,
for some constant M1, using Assumption A3, the boundedness of I111 and dition MC1 and Lemma 2.2), and the boundedness of IlPll (Lemma 3.5). Expanding cj in a Taylor series about x gives
I111 (con-
(o)
where 0
() + o,(x)p +
p(1
< Oj < 0. Using (2.10) and multiplying both sides by p gives
p(cj(O) (sj(O))
is bounded for 0
O)(cj(x) sj) + p 1/202pTV2cj(Oj)p.
Lemma 3.12 implies that plcj(x)- sjl is bounded, Lemma 3.11 implies that pllOpll 2 < c, and Assumption A3 implies that IIV2cj(Oj)ll is also bounded. Consequently,
p l(cj(O)
sj(O)) < N,
(3.52)
where N is a constant. This result and Lemma 3.5 imply the second term in is also bounded, that is,
(3.5a)
Ip(cy(O) sy(O))pTVcy(O)pl <_ N:llPll _< M,
(3.52).
where N2 and M2 are constants. Consider now pllA(O)p- qll 2, the third term on the right-hand side of Using Taylor series, we have
(.)
where 0 < 0j < 0. From
( + op)p
+ OpV()p,
we obtain
(2.10) and Lemmas 3.11 and 3.12,
plIA(O)p- qll < Ma,
(3.56)
618
WALTER MURRAY AND FRANCISCO J. P RIETO
where M3 is a constant.
From (3.55), (2.10), Assumption A3, and the boundedness of the final term on the right-hand side of (3.52) is also bounded,
I111 (Lemma 2.2),
2T(A(O)p q)
(3.57)
--2T(Ap q) +
iOpTV2ci(j)P
< 2T(c- s) + X4llpll z < M4,
where N4 and M4 are constants. [3 The desired bound follows from (3.52), (3.53), (3.54), (3.56), and (3.57). LEMMA 3.14. For any > O, if IlPk[I > there exists a value () such that is the steplength computed by the algorithm. ck > (e) > O, where We drop the subscript k corresponding to the iteration number. We start Proof. by proving that & (as defined in (2.14) and (2.15)) is bounded away from zero if IlPll > e. If condition (2..14)is satisfied at a given iteration, then & 1, trivially bounded away from zero. We assume therefore that & is chosen to satisfy (2.15). In the proof of Lemma 3.1 it was shown that the linesearch procedure was well defined, and in particular, that there exists a value C) E (0, 1] satisfying (2.15) and such that condition (2.15b) is not satisfied for any value of a e [0, C)); see (3.6), (3.8), and (3.7). From the Taylor series expansion of i at C) we have
()
where 0 obtain
(0) + "(0),
r/
<0<
C).
Therefore, using (3.6) and noting that
< 1 and b(0) < 0,
we
(z.s)
,,(0)
( ) ,,(0---S
(Since C) > 0, 0 must be such that "(0) > 0.) If IlPll > e, then from (3.4) we have that I(0)1 > we also have 1/(0) < N, implying
1/2svHe2, and from Lemma 3.13
If condition (2.16) is satisfied for C), then the previous bound holds for c. Otherwise, for some constraint j we must have hi(c)) ci(x + c)p) + tic < 0 (using the notation introduced in Lemma 3.1). If hi(0 > 1/2tic > 0, from the continuity of h 0 and hi(a > 0 for all c [0, 6]. From there exists a value C) < & such that hi(c) the mean-value theorem
hi(&)-hi(O)
hi(O)
h(O)
Ih(0)l
for some /9 e [0,6]. But as Ih(0)l laj(x + Op)Tpl K for some K > 0 (from Assumption A3 and the boundedness of IIPlI, Lemma 3.5), we have
(3.59)
c)
> fie 2K"
A SEQUENTIAL QUADRATIC PROGRAMMING ALGORITHM
619
If hi(0) <
1/2c, we must have from (2.4b),
hj(O)
ap> -cj c
hy(O) >
2
-.
From hy(0) > 0 and hi(d) < 0 there must exist a value & < d such that hi(&) < 0, implying the existence of & < & such that h() 0 and h(c) _> 0 for all c E [0, &] (also, h() > 0 for all e [0, &]). From the mean-value theorem,
for some 0 e [0, 0]. But h(0) _> 1/2/c, and Ih(0)l IpTV2c(x + OP)Pl < [( for some K > 0, from Assumption A3 and the boundedness of IIPlI, Lemma 3.5, implying again
(3.60)
0
>
2R"
The procedure to construct a will ensure that a >_ and so the result presented in the lemma will hold. We can now prove the global convergence theorem for the algorithm.
1/2,
THEOREM 3.15. The sequence {Xk} generated by the algorithm converges to a unique KKT point for NP. O, it is sufficient to Proof. It follows from Lemma 3.9 that to prove Ilxk x* show
(3.61)
lim
IlPkll -+ O.
If (3.61) is true then there exists K such that IlXk < (*/2 and IlPkll < for all k > K, where 2(* is the minimum distance between two KKT points. It follows that x* is unique for k > K (the sequence converges to the unique KKT point nearest to xK), otherwise it implies that for some k > K that either Ilxk > (*/2 or (*- Consequently, to prove the theorem it is sufficient to show (3.61) is true. IlPkll > 0 for any k, the algorithm terminates and the theorem is true. Hence If IlPk we assume that IlPk # 0 for any k. If Pk L 0, there must exist a subsequence {Pl}, and a positive constant (, such that IlPlll > ( for all 1. In this case, from Lemma 3.14 there will exist a uniform lower bound on hi, al > > 0, but then
--x*ll
*
--x*ll
contradicting the fact that pllapll is bounded (Lemma 3.11). In the bounded case, we know that there exists a value t5 and an iteration index K such that p- for all k > K. Again, the proof is by_ contradiction. Consider only indices such that >/. Every such iteration after K must yield a strict decrease in the merit function because the termination condition for the linesearch (2.15a), together with the boundedness of the steplength (from Lemma 3.14 and IlPtll > () and (3.4) imply
_<
_<
--1/2aaZ .HIIP ll = < 0.
The adjustment of the slack variables s in (2.7) can only lead to a further reduction in the merit function, as L A is quadratic in s and the minimizer with respect to sj
620
WALTER MURRAY AND FRANCISCO J. PI:tIETO
is given by cj Aj/p. From the fact that the penalty parameter is not modified, for iterations from the subsequence we have
(x+) () <
-1/2o.
Therefore, since the merit function with p t5 decreases by at least a fixed quantity at every step in the subsequence, it must be unbounded below, contradicting (3.19). It follows that (3.61) must hold. Having established the global convergence of the algorithm, the next step is to A*. In order to prove this result, we need to show that the multiplier estimate Ak strengthen our conditions on the multiplier estimate #k (if #k does not converge then
A will not converge either). Following is the additional condition. MC3. II#k--,k*ll- O(llxk--x*ll), where denotes any multiplier vector associated
with a KKT point closest to xk. This condition requires that
*
in condition MC1 be chosen so that
(3.62)
Estimates satisfying MC1, MC2, and MC3 may be obtained by computing a multiplier for the "active" constraints (say, least-squares estimates of least length), and expanding to the full multiplier space by augmenting this vector with zeros corresponding to the inactive constraints. If such an estimate does not satisfy MC1, then a suitable estimate may be determined by appropriate scaling. The multipliers at the stationary point of the QP also satisfy the three conditions. Note that if x* is not unique then from Lemma 3.6, Ilxk -x* > e for some e > 0, and MC3 holds for any vector #k that is bounded. We first show that under the stronger conditions on # the steplength ak is uniformly bounded away from zero. LEMMn 3.16. Under MC3 and all earlier assumptions and conditions,
O.
Proof. We again drop the subscript k. We first tighten the bound on b"() given Lemma 3.13. From (3.53) and (3.54), we have that the first two terms on the right-hand side of (3.52) are bounded by a multiple of Ilpll 2. For the remaining terms,
in
from
(3.55) and (2.10)we obtain
(p(A(t?)p-q)-2)T(A(O)p-q)
p(cj
E(OpTVcy(y)p--cy+sy--2j)(0pTVcy(y)p--cy+sy).
J
Expanding this expression, and using Lemmas 3.11 and 3.12 to bound the terms sj)OpTVcj(j)p and p82(pTVc(j)p)2, we obtain
(.6)
IIA() 11
2(A(0)p
(3.3) and MC2,
for some constant M. Observe that from
p (c- )r(_ )
+ 2(c_
-b (0) + pT(g_ AT.) #Ts.
(3.64)
Using Taylor expansions and Lemma 3.9 it follows that
pT(g AT#) pT(g*
A.T#)
A SEQUENTIAL QUADRATIC PROGRAMMING ALGORITHM
621
From this result and MC3 there exists a constant
(3.65)
pT(g_ AT#)<_
From #k --* A*, strict complementarity at a KKT point (Assumption A6), and the fact that the correct active set is identified for ]JpJJ small enough (Lemma 3.8), we eventually have # >_ 0 and #Ts >_ 0. Consequently, it follows from (3.52), (3.53), (3.54), (3.63), (3.64), and (3.65)that
(9)
_
such that
-k (0) + Yllpk
2
for some constant N > 0. This result and (3.4) can be used with there exists a value & satisfying (2.15) such that
(3.58) to imply that
> 0.
>_ (1
+ i)ll .ll
+
The desired result then follows from an argument identical to that given in the final part of Lemma 3.14. This lemma also implies that the effort needed to compute the value for the steplength is uniformly bounded in the algorithm. We now establish the convergence of the multiplier estimate. THEOREM 3.17. Under MC3 and all other assumptions and conditions,
lim
Ak
A*.
Proof. From (2.29),
k
(3.66)
j=O
where
(3.67)
with 1 and condition that A0
/kk
a,
1.
tk
a H (1
r--l+l
k
a),
< k,
a
a
at,
(This convention is used because of the speciM initial
we observe that
0.) From Lemma 3.16 and (3.67),
0
(3.68a)
(3.68b)
<
,
k
/=0
a
tk
1
for all l,
1,
(3.68c)
From condition MC3 we have
tk
(1- )k-t,
< k.
(3.69)
k
* + Mkhtk,
]]tk]] K,
1.
with ]Mk M, 5k ]]Xk X*] and can choose a value K1 so that, for k
(3.70)
M5k]
.
om Theorem 3.15, for any e > 0 we
622
WALTER MURRAY AND FRANCISCO J. PRIETO
Given any e > 0, we can also define an iteration index K2 with the following property:
(3.71)
for k _> K2 k_> 2K,
+ I.
Let K
2(k + i)(i + 2) max(Kt,K2). Then, from (3.66) and (3.69),
K
/=0
(i- ()k _<
we have for
k
l=K+l
Hence it follows from (3.68b) that
K
/=0
k
l=K+l
From the bounds on
IIll (condition MC1), IIll, and (3.62), we obtain
K
l=O
k
l=K+l
Since we assume k
K
_> 2K,
K
it follows from
K
(3o68a) and (3.68c) that
0k
_<
(1- )k- -(1- )2K-I <_ (K + 1)(1- )g.
_<
/=o
/=o
Using (3.71), we thus obtain the following bound for the first term on the right-hand side of (3.72)"
K
(3.73)
2/g
/=o
To bound the second term in (3.72), we use (3.6Sb) and (3.70)"
k
k
(3.74)
/=K+I
klM 61 <_
/=K+I
<_
given any e
Combining
(3.72)-(3.74), we obtain the following result:
<_
> 0, we can find
K such that
for k _> 2K + 1,
[:3
which implies the desired result.
4. Rate of convergence. In this section we shall show under additional assumptions on the multiplier estimate that the algorithm converges at a superlinear rate, independently of the asymptotic behavior of the penalty parameter. Since Pk -+ 0, we may assume without loss of generality that Pk has been obtained as the minimizer for the QP subproblem, and that the correct active set has been
identified. We again start by presenting an outline of the steps taken. (i) Bounds on the rate of growth of the penalty parameter introduced in Lemmas 3.10, 3.11, and 3.12 are tightened.
A SEQUENTIAL QUADRATIC PROGRAMMING ALGORITHM
623
In Lemma 4.1 we prove that at all iterations at which Pk is increased (if we
have an infinite sequence of such iterations)
Pllck skll
0 and
pkllPkll 2
--,
O.
In Lemmas 4.2 and 4.3 these results are extended to all iterations.
(ii) In Lemma 4.4 it is shown that #Tsk 0 for sufficiently large k. (iii) Lemma 4.5 proves the superlinear convergence of the sequence {xk +Pk-x*}, under certain assumptions on Hk. (iv) For k sufficiently large, ak 1. Lemma 4.6 gives the relationship between the descent in one iteration Ck(1)Ck(0) and the initial derivative in the linesearch (0).
Theorem 4.7 shows that ak 1 for all sufficiently large k, implying superlinear convergence. (v) Finally, Theorem 4.8 shows that under an additional condition on the multipliers, the penalty parameter remains bounded. The first two lemmas introduce refinements on the results presented in Lemmas 3.10, 3.11, and 3.12, and their proofs are based on the corresponding proofs for these lemmas. LEMMA 4.1. If kl --* x3, where kl denotes an iteration at which the penalty parameter is increased, then
lim p Ilck
Sk
0
and
lim p
IIP
2
0.
Proof. We drop the subscript kt in what follows. Since p is the minimizer of QP, condition (2.5a) holds for a nonnegative vector From (2.4b) and (2.5a) we have gTp / 1/2pTHp _rTc and using this result in the definition of 5, (2.13),
From (2.12) we have p _< 2, and using Theorem 3.17, MC3, and Lemma 3.7 we obtain
l From (3.36) and (4.1) we have lim__, Pk IIPk 2-. O, completing the proof. 4.2. For general iterations k, limk-.o pkllPkll 2 -O. LEMMA Proof. Define I _-- kt and K kg+l. If p is bounded, the result follows from Theorem 3.15. If p is increased in an infinite number of iterations, from (3.38) and Lemma 3.14 we only need to show that
,
--+
O.
p1
From the boundedness of ]lAkll (Lemma 2.2), Lemma 4.1 and the fact that
PK we have
<
We also have from Lemma 4.1,
624
WALTER MURRAY AND FRANCISCO J. PRIETO
These results and the definition of
,
(2.2), imply
p,(,- CK)- px(Fx--FK) 40.
We now analyze the asymptotic behavior of the term px(Fx
F). We have
F F
(c, c)T, +
O(max([Ip, 2,
(3.43) also holds
Using the same arguments as in the proof of Lemma 3.11, inequality in this case, and from (3.15),
A second bound for this term can be obtained from r, >_ 0 and s, >_ 0, implying
Since
I1,11
is bounded, it follows from applying
pTr
T
Lemma 4.1 to (4.3) and (4.4) that
(4.5)
From (2.9), the boundedness of
cx
----+
O.
I111 and Lemma 4.1,
(4.6)
We can again use Lemma 4.1 to obtain
(4.7)
p,O(m x(llp, IIp, ll ))
0.
From (3.42), (4.5), (4.6), and (4.7) we have that the sequence {p,(F F)} is bounded above by a sequence that converges to zero. It then follows from 0 and (4.2) that p(x- ) --+ 0 and the desired result follows from (3.38) and Lemma 3.16. LEMMA 4.3. For general iterations k, limk_. Pk lick Sk O. Proof. If p is bounded the result follows from c* >_ 0, A* _> 0, c* 0, Theorems 3.15 and 3.17 and (2.8). We assume therefore that p is increased an infinite number of times. Consider two cases. Case 1. If constraint j is such that > 0, then j 0 and from (2.8),
AT
c.
but from Theorem 3.17 and Assumptions A3 and A6, eventually ,j < pcj, implying
Case 2. For those j such that -0, implying )j > 0, consider iteration indices large enough that the correct active set is identified (Lemma 3.8), implying + cj 0. From the Taylor series expansion for cj and the boundedness of the steplengh,
c.
ap
C:i(Xk + okPk)
Cj(Xk) + ak(ak)pk + o(ll p, ll
(1
Ok)Cj(Xk) + O(llp ll ).
A SEQUENTIAL QUADRATIC PROGRAMMING ALGORITHM
625
Recurring this relationship for k, I < k < K, we get
but as 0 < at _< 1 we must have
(4.8)
From
P I(c ) I
<- pI(a) I / P
IpI )
pl(cx)j -(s)j
c.
0, Assumptions A3 and A6, and (2.8), eventually it must hold that pxlc()jl, and using Lemma 4.1, (4.8), and Lemma 4.2,
From this result, definition (2.8), Assumptions A3 and A6, and Theorem 3.17, for k large enough
pal(oa)
This completes the proof.
]min(pk(ck), (Ak))l
O.
Ip (oa) l
0.
Proof. If constraint j is such that > 0, then for k large enough (ck)j _> e > 0, and (ak)Pk + (Ck) >_ 1/2e > 0. It therefore follows from MC2 that (#k) 0. If j is such that 0, then from Assumption A6, > 0. Also, from Lemma 4.3,
pk((Ck)j --(Sk)j) min(pa(ck)j, (Ak)j) O, and for large enough k Theorem 3.17 will imply pk(Ck)d <_ (Ak)j; these two results and definition (2.7) imply
LEMMA 4.4. For k large enough #Sk
c.
c.
-
.
(sk)d
completing the result.
max(0, (ck)j
To prove that the algorithm converges superlinearly it is necessary to assume that Hk converges to an approximation of V2xL(x*, *) in some sense, where L(x, denotes the Lagrangian function for problem NP. Define Wk as
(4.9)
Wk
2 =-- V2=L(xk, Ak) V2xF(xk)- y(Ak)VC(Xk).
J
We impose the following additional condition on Hk. HC3. Following Boggs, Tolle, and Wang [3], we assume
Z[(H
where
o( IIp II),
Zk
active at
x*,
is a basis for the null space of k, the Jacobian of xk of those constraints that is bounded in norm and has its smallest singular value bounded
away from 0. The proof proceeds by first showing that the sequence (xk + Pk x* } converges superlinearly, and then proving that a steplength of one is eventually attained.
626
WALTER MURRAY AND FRANCISCO J. PRIETO
The following lemma corresponds to Theorem 3.1 in [3]. LEMMA 4.5. Under Assumptions A1-A7, and conditions MC1-MC3, HC1-HC3,
(4.10)
The results presented on bounds for the growth rate of the penalty parameter allow us to obtain an asymptotic expansion for the quantities involved in the linesearch termination criterion. We want to prove that condition (2.14) is satisfied for k sufficiently large. It is shown in the following lemma that the satisfaction of (2.14) is directly related to the asymptotic properties of T =_ P(gk ATk#k) + pkTwkp. LEMMA 4.6. The following relationship holds:
Ck(1) Ck(0)
with Xk
1/2(0) + 1/2Tk / o(llPkll2).
Proof. In the proof we drop the subscript k, and we denote quantities associated + Pk by a tilde, that is,/ F(xk + pk) while F =_ F(xk). From the definition of the merit function (2.2) and (2.1) we have f #T(5- s q) + AT(c- s) (1) (0) p ,o (4.11) (,_
+
From the Taylor series expansion of c around x and (2.10) we have
5j
sj qj
j
cj
ayp
and using this result with the Taylor expansions for c and F in
(1) (0)
(4.12)
gTp + pTV 2 Fp -j#j pTV2cjp + ,T(c 8) P + -(,pTV2cjp, 2- P (c- s)r(c- s)+ o(llpll =) 8,...,j
#
From (2.6), condition MC3 and Theorem 3.17 we have
(4.13)
, + , + o(1).
Also, from Lemma 4.2 and Assumption A3 we have
-
pT,-,2
(4.11)
we obtain
ppTV2cjp 0(1) and p(pTV2cjp)2 o(llpl12). Replacing these results in (4.12) and reordering the terms we obtain
(i)
(o)
9 + v erV + P ( ,)T( ,) / o(llPll ) + 1/2,(- ,)to simplify this expression,
Using
(4.9) and (3.3)
(1)
(4.14)
(0)
1/2(0) + 1/2 (gTp + pTWp + #T(c
s)) +
Finally, from condition MC2 we have #Tc --#TAp, and from Lemma 4.4 we know that eventually #Ts 0, implying in particular that #Ts o(llpl12), and replacing these bounds in (4.14) we have
(1) 0(0)
1/2(0) + 1/2 (pTWp + pT(g
AT,)) + o(llPl12),
A SEQUENTIAL QUADRATIC PROGRAMMING ALGORITHM
627
completing the result.
The main result of this section is given in the next theorem. It is shown that, if condition MC3 is replaced by a stronger condition, then after a finite number of iterations a steplength of one is taken for all iterations thereafter, implying that the algorithm achieves superlinear convergence. The new condition is MC3. I1 o(11 *11). It is possible to prove superlinear convergence without the need to strengthen the conditions on the multipliers. It is shown in [29] that there exists a constant M such that if Pk > M, condition MC3 is sufficient. THEOREM 4.7. If MC3 and all other assumptions and conditions hold then eventually a unit step is always taken and the algorithm converges superlinearly. Proof. As in Powell and Yuan [28], observe that the continuity of second derivatives gives the following relationships:
*11
(4.15)
F(xk + Pk)
C(Xk + Pk)
1(g(xk) + g(xk + Pk) Pk + o(llPkll2), C(Xk) + 1/2 (Z(xk) + Z(xk + Pk))Pk + o(llPkll2).
F(xk) +
From the Taylor series expansions we have
-
F(z + p) (x + )
and since
ITV2F x )p + o(llr,ll), F(x) + (x)rp + -,, (x) + (x)r + lpTV2 () + o(1111),
(4.10) and Lemma 3.9 imply g(xk + Pk) g* + o(llPkll), ay(xk + Pk) a.i + o(llPkll), we get from (4.15) and (4.16) that (we drop the subscript k)
(4.17a) (4.17b)
pTV2Fp (g* g)Tp + o(llPl12), pTV2cjp (a aj )Tp -t- O(llPl12),
Condition MC3, Theorem 3.17, and (4.13) give -:j Aj pTW2cjp # pTW2cjp + o(llpl12), and if we apply this bound to the result of adding (4.17a) to (4.17b) multiplied
by ,j, we have
(4.18)
Condition MC3,
pTWp pT(g*
(1.1),
A,T#)
pT(g AT#) + o(llplle).
o(llplle),
and Lemma 3.9 imply
pT(g*
and from
A,T#) pTA*T (A*
-
#)
(4.18),
T
(4.19)
pTWp + pT(g AT#) pT(g*
(1) (0)
A,T#) + o(llpil 2)
o(llpl12).
From Lemma 4.6 and (4.19) we get
1/2(0) + o(llpll).
Since (0)
< 0, the above relationship and Theorem 3.15 imply that condition (2.14)
is eventually satisfied for k sufficiently large.
628
WALTER MURRAY AND FRANCISCO J. PRIETO
Regarding condition
(4.20)
for some Oj E
[0, 1], and aj(xk OjPk)Tpk (4.21) o e [0, ].
Using Theorem 3.15 and the boundedness of IIV2cj(xk /jPk)ll A3 and Lemma 3.4) in (4.21), for k large enough
-
(2.16), we can use Taylor series expansions for cj cj(xk q- Pk) Cj(Xk) q- aj(x} q- OjPk)Tpk
aj(xk)Tpk + pV2cj(xk -}- jPk)Pk,
(from Assumption
aj(xk q-Ojpk)Tpk
and from
_
to write
aj(xk)Tpk 1/2c,
(2.4b),
aj(xk + Oypk)Tpk > aj(xk)Tpk
1 > _()_ 1/2. 2
Replacing this bound in (4.20), we obtain for all k large enough c(xk + Pk) >_ --1/2ce, and condition (2.16) will also be satisfied, giving Xk+l xk + Pk. The required result then follows from Lemma 4.5.
4.1. Boundedness of the penalty parameter. The last result in this section shows that, if condition MC3 is replaced by a slightly stronger condition, the penalty parameter needs to be modified in at most a finite number of iterations (and consequently it remains bounded). The criterion presented will be satisfied, for example, by the least-squares multipliers computed at xk + Pk. THEOREM 4.8. If the multiplier estimates #k in the algorithm satisfy
(a.)
and all other assumptions and conditions hold then there exists a constant M such that Pk M for all k. Proof. We may assume k large enough so that ck 1. From (2.5), (2.4b), and ,-Ts k_ 0, we have rk >
gkPk
_
I1 *11 o(llx +
+ PkHkPk pAkTrk --cTrk <_ --(ci 8k)TTrk,
where rk denotes the QP multipliers at iteration k. From that a unit steplength is accepted, it follows that
(4.24)
(0) <_ --pkgkpkT / tl2, I1 11 pll 11 From (4.22), HC2, and Lemmas 3.9, 3.8, and 3.7 we must have
112#k-1- #k rkll <_ M1 IlPkll
I]2#-1
#
-
(3.3), (4.23), and the fact
, ,
<-- M2pHkpk
for some positive constants M1, M2. It then follows using a 2 -F b2
_
.
2ab that
I_.T llll 11 < M2 V/pHpII 11 < v Hw / M2 I1 11
,
implying from
(4.24) that _/_/ I(0) < -,
, + (M
1/2M22, condition (2.11) will be satisfied, and the penalty parameter will not be increased. Given that we are using the rule (2.12) for updating Pk, it must hold that Pk
From this inequality it follows that if p >_
A SEQUENTIAL QUADRATIC PROGRAMMING ALGORITHM
629
5. Other merit functions. Several merit functions have been proposed and analyzed in the literature (a review can be found in Powell [27]). The question arises if the convergence results using early termination in the solution of the QP subproblem depend on our specific merit function, or if they are fairly independent of this choice.
We shall show in this section that the choice of merit function is not critical. What we present is how to adapt our SQP algorithm to the use of other merit functions rather than examine other methods explicitly to see if the particular QP subproblem posed
and the manner the search is performed can be adapted to the use of an incomplete solution. For example, we still perform a search in the x and A spaces. Slack variables do not appear in the merit functions we shall consider, consequently the search in the space of the slack variables is no longer required. We have selected as examples the study of two particular merit functions. The first one corresponds to a class of merit functions that includes among others the merit function analyzed in Han [21], Byrd and Nocedal [5], and Burke and Han [4]. This general merit function takes the form:
(5.1)
=_ max(0,-cj(x)). Again, we where an gp-norm will omit the subscript if we refer to the g2-norm, and we will explicitly include it whenever we refer to a general gp-norm. The second merit function we consider is
(x, A) F(x) (1 _< p _< cx)) is used, and c-(x)
(5.2)
(x,
where we use the 2-norm. This merit function has been studied among others by Powell and Yuan [28] (applied to the equality-constrained problems only) and Schittkowski [32]. Unlike either of these algorithms, where the multiplier estimate was treated as a function of the iterate A.(x), we do not explicitly define the form of the multiplier estimates although the ones used in both methods satisfy the criteria MC1, MC2, and MC3. Indeed the one used in [28] also satisfies MC3. We still assume A1-A7 hold for the problem. However, when the merit function (5.1) is used, the multiplier estimate #k is only required to satisfy MC1. This condition is trivial to satisfy. For example, we may choose 0 0 and #k 0 making the search in the multiplier space void. Such a choice reduces (5.1) to the well-known gl merit function and our algorithm becomes very similar to that analyzed in [21]. When (5.2) is used, we assume conditions MC1 and MC2 hold. We have also assumed in the proofs that A0 _> 0 and #k _> 0. We omit the proofs that the iterates lie on a compact set. For the first merit function (5.1) this proof is relatively straightforward, since it will be shown that the penalty parameter is bounded. The proof for the second merit function (5.2) is very similar to that for the Augmented Lagrangian merit function. The criteria (2.15) for the choice of steplength ak assume the merit function has continuous first derivatives. This property does not necessarily hold for the merit functions under consideration. Therefore we use the following criteria for determining a value Define
(5.3)
Ak =-- [Pk + (k Xk)Tc-(Xk) PkllC-(Xk)llp.
Ck(&k) --= (Xk + &kPk, )k +
We start by selecting a value &k satisfying
(5.4)
630
WALTER MURRAY AND FRANCISCO J. PRIETO
and either
(5.5)
or
& >_ 7 > 0 & > 7k
and
(5.6)
a < 1 and k > 0. For a discussion of these criteria where 0 < 7t < 7u < 1, 0 < and the estence of see Calamai and Mor [6]. In addition to these conditions, we also also want to limit the size of the infeibilities. If 5k satisfies condition (2.16), then we let ak k. Otherwise, we compute ak by performing a backtracking linesearch om k until conditions (5.4) and (2.16) are both satisfied. Our preference for the criteria given in 2 is bed on our belief that in practice they lead to a better choice of a. In the definition of our algorithm we could have used other steplength criteria without impacting the convergence properties. The following basic relationships will be used to establish the convergence results,
(5.7a) (5.7b) (5.7c) (5.7d)
c (x + ap)
c(x + p) c(x)
aap
min(0, cj(x) + aap)
nin(O, c(x) + aap) (1 a)c (x),
--wTAp --]]c- (x)],
-Ap
-c- (x).
Vc(x). Also, is a diagonal matrix such that -Ap is an element of the subdifferential of c-(x + take values in [0, 1], are zero whenever cj(x) > 0 and take the value one whenever cj(z) < 0. Finally, wTAp represents an element of 0(0), the subdifferential of c-(x + ap)p at 0. The elements of w are given by
In these inequalities A
and have the property that w Tc(x) -]c-(x)]p. Consider now the case when has been defined from that Ak 0 and (2.4b),
(5.1). om our assumption
Ak(Akpk + ck)
(0)
(5.8)
0
for all k. It follows from this inequality and the relationships given in
(5.7) that
g[Pk + [c-(Xk) AkAp pkw[AkPk
T
We select Pk such that
This rule is analogous to the ones used in Byrd and Nocedal
[5], and Burke and Han
The first step is to establish that such a value of p exists.
we have
om (3.14) and (5.3)
g
+
+
A SEQUENTIAL QUADRATIC PROGRAMMING ALGORITHM
631
If we now use (2.6), property MC1, and Lemma 2.2 to bound the multiplier term
h
h ud
Defining Pu k. This result also shows that the value of p will remain bounded in the algorithm. THEOREM 5.1. The algoHthm modified to use the merit function (5.1) converges globally. O. Proof. Given the bound in Lemma 3.9, it suffices to show that ]Pk As p cannot grow without bound, any strategy for increing p by a finite quantity whenever it is required to increase p implies that there exists an iteration value K such that Pk PK for all k K. We consider only iterations of this form. For k K, from (5.4), (5.8) and condition MC2,
I1i1 1111,, w obtain in (5.9) )PkHkPk + (2 + 3, p)IIZiI. (2 + 3), for any value p p condition (5.8) is satisfied for any
(a)- (a-l)
om the boundedness of
(5.10)
If []p
(Assumption A3), it follows that
0, convergence follows from Lemma 3.9. Otherwise, if for a subsequence 0 along the subsequence, and from the we must have ak termination conditions for the linesearch (5.4), (5.5), and (5.6), k 0, as the step required to satis condition (2.16) is uniformly bounded away from zero (see (3.59)
]]Pk]] > e, from (5.10)
-
av
1I
0.
-uZ..allpll
.
and
(3.60)).
Finally, from
(5.6)
we must also have
0.
In the following relationships we drop the subscript k corresponding to the iteration number, and we denote by a tilde the value of functions evaluated at x + p (i.e.,
e (z + p)). om the definition of the merit function (5.1),
() (0)
a + r(e+(
c-) + llc-I, F ) + (lle- I1 (1 )llc- II,).
re-
For the last term, from (5.7a) and (5.7b), it follows that
(,) (o) < (gTp + AT(e- c-) + C T- cwlIc-
If we use again (5.7a) and (5.7b) on the terms associated with the multiplier estimates (given that by assumption A 4- >_ 0), and the Taylor series expansions for F and c,
we obtain
() (0)
_
+(
F
CgTp) + pile
c
AplI,,.
gTp / Ej(j / bj)lj cj ap 4- (1 ()ATc ATc + (1 C)Tc PlI- II,, + O(IIPlI).
632
WALTER MURRAY AND FRANCISCO J. PRIETO
After simplifying this expression we have
() (0) < (r + ( )%- llc-iI) + vllllll- ii + Replacing this bound in (5.6) implies 0 < (1 a)A + 2l[[c-[ p + O(l[p2).
Since from (5.8) and condition HC2, A --svU]p]l 2, and we have sumed that lipS] > e, it follows by taking limits along the subsequence that
0
-(1 a)He2.
0 for the whole sequence.
However, this is not possible, implying ]Pk
Consider now the second merit function direction at (Xk, Ak) is given by
g}Pk
(5.2).
The subgradient along the search
+ [c-(x}) ADkA}pk pkc-(xk)TAkpk
g[p +
where
=_
A )Tc-
p ll -(x )ll
Note that Ak >_ 0 implies
(kik + pkc-)T(Akpk + C)
O.
In this case it is not immediately evident that Pk remains bounded. The convergence proof we give is similar to the one introduced in 3. The definition of p given in that section will be preserved, except c- s is replaced by c-. THEOREM 5.2. The algothm modified to use the merit function (5.2) converges
globally.
O. Proof. Again, from Lemma 3.9 it is enough to show that Pk First sume that p is bounded. The argument used is similar to the one in Theorem 5.1. om (5.4), (5.8), condition MC2 and the boundedness of (5.10) must hold also for this case. If ]]pk O, convergence follows om Lemma 3.9. Otherwise, if for a subsequence 0, and from condition (5.6) and the > e, from (5.10) we must have ak Pk boundedness of the step to satis (2.16), Gk 0. om (5.2), (5.7a) and (5.7b), we also have (we again drop the index k in the following relationships, and use a tilde to indicate values at x + p)
,
+
and again using (5.7a) and (5.7b) on the terms associated with the multiplier estimates, we obtain
(a)
(0)
a(gTp + (
)Tc- pllc-II 2) + 211c-II (11511 + pll -II) + O(ll pll )
A SEQUENTIAL QUADRATIC PROGRAMMING ALGORITHM
633
Replacing this bound in
0
(5.6)
implies
< (1 --o)&A 4" =ll-II
(ll :ll + 1/2pll -II) + o(ll p[l=).
Since from (5.8) and condition HC2, z _< --/.gllpll and we have assumed that IlPll > e and p is bounded, it follows by taking limits along the subsequence that
,
0 <_
-(1 a)/l,He.2.
However, this is not possible, which implies Ilpa I 0 for the whole sequence. Assume now that pk grows without bound. In this case we have that for all
iterations where the value of the penalty parameter is increased
The proof of this result is basically that of Lemma 3.10. From these bounds it is possible to show that we must also have
for all k (the proof is similar to the one for Lemma D convergence of the algorithm.
3.11), implying p
--+
0 and the
6. Numerical results. In this section we present numerical results obtained from an implementation of our algorithm. As a first step we have modified the code NPSOLo We have called the modified routine INPSOL. Apart from the definition of the search direction all other aspects of INPSOL are identical to those of NPSOL. A detailed description of NPSOL is given in Gill et al. [15]. It should be noted that NPSOL does not incorporate linear constraints into the merit function. An initial point is obtained that is feasible with respect to the linear constraints and thereafter feasibility is retained (by incorporating the linear constraints in the QP subproblem). On many practical problems the feasible region with respect to the linear constraints is compact. On such problems this approach ensures Assumption A2 is s...