A Feasible BFGS Interior Point Algorithm for Solving Convex Minimization Problems.
We propose a BFGS primal-dual interior point method for minimizing a convex function on a convex set defined by equality and inequality constraints. The algorithm generates feasible iterates and consists in computing approximate solutions of the optimality conditions perturbed by a sequence of positive parameters $\mu$ converging to zero. We prove that it converges q-superlinearly for each fixed $\mu$. We also show that it is globally convergent to the analytic center of the primal-dual optimal set when $\mu$ tends to 0 and strict complementarity holds.
Introduction
. We consider the problem of minimizing a smooth convex function
on a convex set defined by inequality constraints. The problem is written as
# R is the function to minimize and c(x) # 0 means that each component
m) of c must be nonnegative at the solution. To simplify
the presentation and to avoid complicated notation, the case when linear equality
constraints are present is discussed at the end of the paper. Since we assume that the
components of c are concave, the feasible set of this problem is convex.
The algorithm proposed in this paper and the convergence analysis require that
f and c are di#erentiable and that at least one of the functions f , -c (1) , . , -c (m) is
strongly convex. The reason for this latter hypothesis will be clarified below. Since the
algorithm belongs to the class of interior point (IP) methods, it may be well suited for
problems with many inequality constraints. It is also more e#cient when the number
of variables remains small or medium, say, fewer than 500, because it updates n - n
matrices by a quasi-Newton (qN) formula. For problems with more variables, limited
memory BFGS updates [39] can be used, but we will not consider this issue in this
paper.
Our motivation is based on practical considerations. During the last 15 years
much progress has been realized on IP methods for solving linear or convex minimization
problems (see the monographs [29, 10, 38, 44, 23, 42, 47, 49]). For nonlinear
convex problems, these algorithms assume that the second derivatives of the functions
used to define the problem are available (see [43, 35, 36, 12, 38, 26]). In practice, how-
# Received by the editors September 15, 1998; accepted for publication (in revised form) January
26, 2000; published electronically August 3, 2000.
http://www.siam.org/journals/siopt/11-1/34472.html
des Sciences, 123 av. A. Thomas, 87060 Limoges Cedex, France (Paul.Armand@
unilim.fr).
# INRIA Rocquencourt, BP 105, 78153 Le Chesnay Cedex, France (Jean-Charles.Gilbert@inria.fr).
- MIP, UFR MIG, Universit-e Paul Sabatier, 118 route de Narbonne, 31062 Toulouse Cedex 4,
France (jan@mip.ups-tlse.fr).
200 P. ARMAND, J. CH. GILBERT, AND S. JAN-J -
EGOU
ever, it is not uncommon to find situations where this requirement cannot be satis-
fied, in particular for large scale engineering problems (see [27] for an example, which
partly motivates this study and deals with the estimation of parameters in a three
phase flow in a porous medium). Despite the possible use of computational di#eren-
tiation techniques [8, 19, 3, 28], the computing time needed to evaluate Hessians or
Hessian-vector products may be so large that IP algorithms using second derivatives
may be unattractive.
This situation is familiar in unconstrained optimization. In that case, qN tech-
niques, which use first derivatives only, have proved to be e#cient, even when there
are millions of variables (see [32, 20] and [9] for an example in meteorology). This fact
motivates the present paper, in which we explore the possibility of combining the IP
approach and qN techniques. Our ambition remains modest, however, since we confine
ourselves to the question of whether the elegant BFGS theory for unconstrained
convex optimization [41, 6] is still valid when inequality constraints are present. For
the applications, it would be desirable to have a qN-IP algorithm in the case when
f and -c are nonlinear and not necessarily convex. We postpone this more di#cult
subject for future research (see [21, 48] for possible approaches).
Provided the constraints satisfy some qualification assumptions, the Karush-
Kuhn-Tucker (KKT) optimality conditions of problem (1.1) can be written (see [17],
for example) as follows: there exists a vector of multipliers # R m such that
where #f(x) is the gradient of f at x (for the Euclidean scalar product), #c(x) is a
matrix whose columns are the gradients #c (i) (x), and is the
diagonal matrix, whose diagonal elements are the components of c. The Lagrangian
function associated with problem (1.1) is defined on R n
Since f is convex and each component c (i) is concave, for any fixed # 0, #) is
a convex function from R n to R. When f and c are twice di#erentiable, the gradient
and Hessian of # with respect to x are given by
Our primal-dual IP approach is rather standard (see [24, 36, 35, 11, 12, 1, 26, 25,
15, 7, 5]). It computes iteratively approximate solutions of the perturbed optimality
system
for a sequence of parameters - > 0 converging to zero. In (1.2),
is the vector of all ones whose dimension will be clear from the context. The last
inequality means that all the components of both c(x) and # must be positive. By
perturbing the complementarity equation of the KKT conditions with the parameter
A BFGS INTERIOR POINT ALGORITHM 201
-, the combinatorial aspect of the problem, inherent in the determination of the active
constraints or the zero multipliers, is avoided. We use the word inner to qualify those
iterations that are used to find an approximate solution of (1.2) for fixed -, while an
outer iteration is the collection of inner iterations corresponding to the same value
of -.
The Newton step for solving the first two equations in (1.2) with fixed - is the
solution
of the linear system
d #f(x) +#c(x)#
in which
xx #(x, #) and This direction is sometimes
called the primal-dual step, since it is obtained by linearizing the primal-dual system
(1.2), while the primal step is the Newton direction for minimizing in the primal
variable x the barrier function
log c (i) (x)
associated with (1.1) (the algorithms in [16, 33, 4] are in this spirit). The two problems
are related since, after elimination of #, (1.2) represents the optimality conditions of
the unconstrained barrier problem
c(x) > 0.
As a result, an approximate solution of (1.2) is also an approximate minimizer of the
barrier problem (1.4). However, algorithms using the primal-dual direction have been
shown to present a better numerical e#ciency (see, for example, [46]).
In our algorithm for solving (1.2) or (1.4) approximately, a search direction d
is computed as a solution of (1.3) in which M is now a positive definite symmetric
matrix approximating # 2
xx #(x, #) and updated by the BFGS formula (see [14, 17] for
material on qN techniques). By eliminating d # from (1.3) we obtain
Since the iterates will be forced to remain strictly feasible, i.e., (c(x), #) > 0, the
positive definiteness of M implies that d x is a descent direction of # - at x. Therefore,
to force convergence of the inner iterates, a possibility could be to force the decrease
of # - at each iteration. However, since the algorithm also generates dual variables #,
we prefer to add to # - the function (see [45, 1, 18])
log # (i) c (i) (x) #
to control the change in #. This function is also used in [30, 31] as a potential function
for nonlinear complementarity problems. Even though the map (x, # -
V(x, #) is not necessarily convex, we will show that it has a unique minimizer, which is
the solution of (1.2), and that it decreases along the direction
this primal-dual merit function can be used to force the convergence of the pairs
(x, #) to the solution of (1.2), using line-searches. It will be shown that the additional
P. ARMAND, J. CH. GILBERT, AND S. JAN-J -
EGOU
function V does not prevent unit step-sizes from being accepted asymptotically, which
is an important point for the e#ciency of the algorithm.
Let us stress the fact that our algorithm is not a standard BFGS algorithm for
solving the barrier problem (1.4), since it is the Hessian of the Lagrangian that is
approximated by the updated matrix M , not the Hessian of # - . This is motivated by
the following arguments. First, the di#erence between # 2
involves first derivatives only. Since these derivatives are considered to be available,
they need not be approximated. Second, the Hessian # 2
xx #, which is approximated by
M , is independent of - and does not become ill-conditioned as - goes to zero. Third,
the approximation of # 2
# obtained at the end of an outer iteration can be used as
the starting matrix for the next outer iteration. If this looks attractive, it has also the
inconvenience of restricting the approach to (strongly) convex functions, as we now
explain.
After the computation of the new iterates
the step-size given by the line-search), the matrix M is updated by the BFGS formula
using two vectors # and #. Since we want the new matrix M+ to be an approximation
of # 2
it satisfies the qN equation M+ # (a property of the
BFGS formula), it makes sense to define # and # by
The formula is well defined and generates stable positive definite matrices provided
these vectors satisfy # > 0. This inequality, known as the curvature condition,
expresses the strict monotonicity of the gradient of the Lagrangian between two successive
iterates. In unconstrained optimization, it can always be satisfied by using
the Wolfe line-search, provided the function to minimize is bounded below. If this is
a reasonable assumption in unconstrained optimization, it is no longer the case when
constraints are present, since the optimization problem may be perfectly well defined
even when # is unbounded below. Now, assuming this hypothesis on the boundedness
of # would have been less restrictive than assuming its strong convexity, but it is not
satisfactory. Indeed, with a bounded below Lagrangian, the curvature condition can
be satisfied by the Wolfe line-search as in unconstrained optimization, but near the
solution the information on # 2
collected in the matrix M could come from a region
far from the optimal point, which would prevent q-superlinear convergence of the it-
erates. Because of this observation, we assume that f or one of the functions -c (i)
is strongly convex, so that the Lagrangian becomes a strongly convex function of x
for any fixed # > 0. With this assumption, the curvature condition will be satisfied
independently of the kind of line-search techniques actually used in the algorithm.
The question whether the present theory can be adapted to convex problems, hence
including linear programming, is puzzling. We will come back to this issue in the
discussion section.
A large part of the paper is devoted to the analysis of the qN algorithm for solving
the perturbed KKT conditions (1.2) with fixed -. The algorithm is detailed in the
next section, while its convergence speed is analyzed in sections 3 and 4. In particular,
it is shown that, for fixed - > 0, the primal-dual pairs (x, #) converge q-superlinearly
toward a solution of (1.2). The tools used to prove convergence are essentially those of
A BFGS INTERIOR POINT ALGORITHM 203
the BFGS theory [6, 13, 40]. In section 5, the overall algorithm is presented and it is
shown that the sequence of outer iterates is globally convergent, in the sense that it is
bounded and that its accumulation points are primal-dual solutions of problem (1.1).
If, in addition, strict complementarity holds, the whole sequence of outer iterates
converges to the analytic center of the primal-dual optimal set.
2. The algorithm for solving the barrier problem. The Euclidean or # 2
norm is denoted by #. We recall that a function # : R n
# R is said to be
strongly convex with modulus # > 0, if for all (x, y) # R n
equivalent definitions, see, for example,
[22, Chapter IV]). Our minimal assumptions are the following.
Assumption 2.1. (i) The functions f and -c (i) (1 # i # m) are convex and
di#erentiable from R n to R and at least one of the functions f , -c (1) , . , -c (m) is
strongly convex. (ii) The set of strictly feasible points for problem (1.1) is nonempty,
i.e., there exists x # R n such that c(x) > 0.
Assumption 2.1(i) was motivated in section 1. Assumption 2.1(ii), also called
the (strong) Slater condition, is necessary for the well-posedness of a feasible interior
point method. With the convexity assumption, it is equivalent to the fact that the
set of multipliers associated with a given solution is nonempty and compact (see
[22, Theorem VII.2.3.2], for example). These assumptions have the following clear
consequence.
Lemma 2.2. Suppose that Assumption 2.1 holds. Then, the solution set of problem
(1.1) is nonempty and bounded.
By Lemma 2.2, the level sets of the logarithmic barrier function # - are compact,
a fact that will be used frequently. It is a consequence of [16, Lemma 12], which we
recall for completeness.
Lemma 2.3. Let f : R n
# R be a convex continuous function and c : R n
be a continuous function having concave components. Suppose that the set {x # R
c(x) > 0} is nonempty and that the solution set of problem (1.1) is nonempty and
bounded. Then, for any # R and - > 0, the set
log c (i) (x) #
is compact (and possibly empty).
Let x 1 be the first iterate of our feasible IP algorithm, hence satisfying c(x 1 ) > 0,
and define the level set
Lemma 2.4. Suppose that Assumption 2.1 holds. Then, the barrier problem (1.4)
has a unique solution, which is denoted by - x - .
Proof. By Assumption 2.1, Lemma 2.2, and Lemma 2.3, L P
1 is nonempty and
compact, so that the barrier problem (1.4) has at least one solution. This solution
is also unique, since # - is strictly convex on {x # R Indeed, by
Assumption 2.1(i), # 2 # - (x) given by (1.6) is positive definite.
To simplify the notation we denote by
z := (x, #)
a typical pair of primal-dual variables and by Z the set of strictly feasible z's:
204 P. ARMAND, J. CH. GILBERT, AND S. JAN-J -
EGOU
The algorithm generates a sequence of pairs (z, M ), where z # Z and M is a
positive definite symmetric matrix. Given a pair (z, M ), the next one (z
obtained as follows. First
is a step-size and is the unique solution of (1.3). The
uniqueness comes from the positivity of c(x) and from the positive definiteness of M
(for the unicity of d x , use (1.5)). Next, the matrix M is updated into M+ by the
where # and # are given by
This formula gives a symmetric positive definite matrix M+ , provided M is symmetric
positive definite and # > 0 (see [14, 17]). This latter condition is satisfied because
of the strong convexity assumption. Indeed, since at least one of the functions f or
-c (i) is strongly convex, for any fixed # > 0, the function x #(x, #) is strongly
convex, that is, there exists a constant # > 0 such that
Since # sizes the displacement in x and #, the merit function used to estimate
the progress to the solution must depend on both x and #. We follow an idea of
Anstreicher and Vial [1] and add to # - a function forcing # to take the value -C(x)
The merit function is defined for z = (x, # Z by
where
log # (i) c (i) (x) # .
Note that
# .
Using # - as a merit function is reasonable provided the problem
z # Z
A BFGS INTERIOR POINT ALGORITHM 205
has for unique solution the solution of (1.2) and the direction descent
direction of # - . This is what we check in Lemmas 2.5 and 2.6 below.
Lemma 2.5. Suppose that Assumption 2.1 holds. Then, problem (2.4) has a
unique solution -
z -x - ), where -
x - is the unique solution of the barrier problem
# - has its ith component defined by ( - ) (i) := -/c (i) (-x - ). Furthermore,
- has no other stationary point than - z - .
Proof. By optimality of the unique solution -
x - of the barrier problem (1.4)
x such that c(x) > 0.
On the other hand, since log t is minimized at
c (i) (- x - ) (i) for all index i, we have
z # Z.
Adding up the preceding two inequalities gives # -z - (z) for all z # Z. Hence
z - is a solution of (2.4).
It remains to show that - z - is the unique stationary point of # - . If z is stationary,
it satisfies
Canceling # from the first equality gives #f(x) -#c(x)C(x) thus
x - is the unique minimizer of the convex function # - . Now, by the
second equation of the system above.
Lemma 2.6. Suppose that z # Z and that M is symmetric positive definite. Let
be the solution of (1.3). Then
so that d is a descent direction of # - at a point z
z - , meaning that # - (z) # d < 0.
Proof. We have # - (z) # d. Using (1.5),
which is nonpositive. On the other hand, when d satisfies the second equation of (1.3),
one has (see [1])
which is also nonpositive. The formula for # - (z) # d given in the statement of the
lemma follows from this calculation. Furthermore, # - (z) # d < 0, if z #= -
z - .
We can now state precisely one iteration of the algorithm used to solve the perturbed
KKT system (1.2). The constants # ]0, 1[ and 0 < # < 1 are given
independently of the iteration index.
206 P. ARMAND, J. CH. GILBERT, AND S. JAN-J -
EGOU
Algorithm A - (for solving (1.2); one iteration).
At the beginning of the iteration, the current iterate
supposed available, as well as a positive definite matrix M approximating
the Hessian of the Lagrangian # 2
xx #(x, #).
1. Compute d := (d x , d # ), the solution of the linear system (1.3).
2. Compute a step-size # by means of a backtracking line search.
2.0.
2.1. Test the su#cient decrease condition:
2.2. If (2.5) is not satisfied, choose a new trial step-size # in [#] and
go to Step 2.1. If (2.5) is satisfied, set z
3. Update M by the BFGS formula (2.1) where # and # are given by (2.2).
By Lemma 2.6, d is a descent direction of # - at z, so that a step-size # > 0 satisfying
(2.5) can be found. In the line-search, it is implicitly assumed that (2.5) is not
satisfied if z holds for the new iterate z + .
We conclude this section with a result that gives the contribution of the line-search
to the convergence of the sequence generated by Algorithm A - . It is in the
spirit of a similar result given by Zoutendijk [50] (for a proof, see [6]). We say that
a function is C 1,1 if it has Lipschitz continuous first derivatives. We denote the level
set of # - determined by the first iterate z
Lemma 2.7. If - is C 1,1 on an open convex neighborhood of the level set L PD
there is a positive constant K such that for any z # L PD
determined by the
line-search in Step 2 of Algorithm A - , one of the following two inequalities holds:
It is important to mention here that this result holds even though - may not be
defined for all positive step-sizes along d, so that the line-search may have to reduce
the step-size in a first stage to enforce feasibility.
3. The global and r-linear convergence of Algorithm A- . In the convergence
analysis of BFGS, the path to q-superlinear convergence traditionally leads
through r-linear convergence (see [41, 6]). In this section, we show that the iterates
generated by Algorithm A - converge to -
z -x - ), the solution of (1.2), with that
convergence speed. We use the notation
Our first result shows that, because the iterates (x, #) remain in the level set L PD
the sequence {(c(x), #)} is bounded and bounded away from zero.
A BFGS INTERIOR POINT ALGORITHM 207
Lemma 3.1. Suppose that Assumption 2.1 holds. Then, the level set L PD
1 is
compact and there exist positive constants K 1 and K 2 such that
1 .
Proof. Since # c(x) -
log(# (i) c (i) (x)) is bounded below by m-(1 - log -),
there is a constant K
1 . By
Assumption 2.1 and Lemma 2.3, the level set L #
1 } is
compact. By continuity, c(L # ) is also compact, so that c(x) is bounded and bounded
away from zero for all z # L PD
1 .
What we have just proven implies that {# -
1 } is bounded
below, so that there is a constant K #
1 . Hence the #-components of the z's in L PD
are bounded and
bounded away from zero.
We have shown that L PD
1 is included in a compact set. Now, it is itself compact
by continuity of # - .
The next proposition is crucial for the technique we use to prove global convergence
(see [6]). It claims that the proximity of a point z to the unique solution of (2.4)
can be measured by the value of # - (z) or the norm of its gradient # - (z). In unconstrained
optimization, the corresponding result is a direct consequence of strong
convexity. Here, # - is not necessarily convex, but the result can still be established
by using Lemma 2.5 and Lemma 3.1. The function # - is nonconvex, for example,
when is minimized on the half-line of nonnegative real numbers.
Proposition 3.2. Suppose that Assumption 2.1 holds. Then, there is a constant
a > 0 such that for any z # L PDa#z - z - # 2
Proof. Let us show that # - is strongly convex in a neighborhood of - z - . Using
(2.3) and the fact that -
# -e, the Hessian of # - at - z - can be written as
# .
From Assumption 2.1, for fixed # > 0, the Lagrangian is a strongly convex function
in the variable x. It follows that its Hessian with respect to x is positive definite at
us show that the above matrix is also positive definite. Multiplying the
matrix on both sides by a vector (u, v) # R n
positive definite and c(-x - ) > 0, this quantity is nonnegative. If
it vanishes, one deduces that next that
positive
definite.
Let us now prove a local version of the proposition: there exist a constant a # > 0
and an open neighborhood N # Z of - z - such that
a #z - z - # 2
208 P. ARMAND, J. CH. GILBERT, AND S. JAN-J -
EGOU
The inequality on the left comes from the fact that # -z - and the strong
convexity of # - near - z - . For the inequality on the right, we first use the local convexity
of for an arbitrary z near -
z -z - (z) z). With the
Cauchy-Schwarz inequality and the inequality on the left of (3.2), one gets
a #
# 1Simplifying and squaring give the inequality on the right of (3.2).
To extend the validity of (3.2) for all z # L PD
su#ces to note that, by virtue
of Lemma 2.5, the ratios
and - (z) -z - )
are well defined and continuous on the compact set L PD
z - is the unique
minimizer of - on L PD
2.5), the ratios are respectively bounded away from
zero and bounded above on L PD
1 \ N , by some positive constants K # 1 and K # 2 . The
conclusion of the proposition now follows by taking
The proof of the r-linear convergence rests on the following lemma, which is part
of the theory of BFGS updates. It can be stated independently of the present context
(see Byrd and Nocedal [6]). We denote by # k the angle between M k # k and
and by # the roundup operator:
Lemma 3.3. Let {M k } be positive definite matrices generated by the BFGS formula
using pairs of vectors {(# k , # k )} k#1 , satisfying for all k # 1
where a 1 > 0 and a 2 > 0 are independent of k. Then, for any r # ]0, 1[, there exist
positive constants b 1 , b 2 , and b 3 , such that for any index k # 1,
for at least #rk# indices j in {1, . , k}.
The assumptions (3.3) made on # k and # k in the above lemma are satisfied in
our context. The first one is due to the strong convexity of one of the functions f ,
-c (1) , . , -c (m) , and to the fact that # is bounded away from zero (Lemma 3.1).
When f and c are C 1,1 , the second one can be deduced from the Lipschitz inequality,
the boundedness of # (Lemma 3.1), and the first inequality in (3.3).
Theorem 3.4. Suppose that Assumption 2.1 holds and that f and c are C 1,1
functions. Then, Algorithm A - generates a sequence {z k } converging to -
z - r-linearly,
meaning that lim sup k#z k -
z - # 1/k < 1. In particular,
z - #.
Proof. We denote by K #
positive constants (independent of the iteration
index). We also use the notation
A BFGS INTERIOR POINT ALGORITHM 209
The bounds on (c(x), #) given by Lemma 3.1 and the fact that f and c are C 1,1
imply that # - is C 1,1 on some open convex neighborhood of the level set L PD
1 , for
example, on
where O is an open bounded convex set containing L PD
1 (this set O is used to have #c
bounded on the given neighborhood).
Therefore, by the line-search and Lemma 2.7, there is a positive constant K #
1 such
that either
or
Let us now apply Lemma 3.3: fix r # ]0, 1[ and denote by J the set of indices j
for which (3.4) holds. Using Lemma 2.6 and the bounds from Lemma 3.1, one has for
# .
Let us denote by K #
4 a positive constant such that #c(x)# K # 4 for all x # L PD
1 . By
using (2.3), (1.5), and the inequality (a
and also, by (1.3),
#d x
P. ARMAND, J. CH. GILBERT, AND S. JAN-J -
EGOU
Combining these inequalities with (3.5) or (3.6) gives for some positive constant K #and for any j # J
The end of the proof is standard (see [41, 6]). Using Proposition 3.2, for j # J ,
# [0, 1[. On the other hand, by the line-search, # - (z k+1
# 1. By Lemma 3.3, |[1, k] # J | #rk# rk, so
that the last inequality gives for any k # 1
where K #
7 is the positive constant using the inequality on
the left in (3.1), one has for all k # 1
from which the r-linear convergence of {z k } follows.
4. The q-superlinear convergence of Algorithm A- . With the r-linear convergence
result of the previous section, we are now ready to establish the q-superlinear
convergence of the sequence {z k } generated by Algorithm A - . By definition, {z k }
converges q-superlinearly to - z - if the following estimate holds:
z
z - #),
which means that #z k+1 - z - #z k - z - # 0 (assuming z k #= -
z - ). To get this result,
f and c have to be a little bit smoother, namely twice continuously di#erentiable near
x - . We use the notation
We start by showing that the unit step-size is accepted asymptotically by the line-search
condition (2.5), provided the updated matrix M k becomes good (or su#ciently
large) in a sense specified by inequality (4.1) below and provided the iterate z k is
su#ciently close to the solution -
z - .
Given two sequences of vectors {u k } and {v k } in some normed spaces and a
positive number #, we write u k # o(#v k # ), if there exists a sequence of {# k } # R
such that # k # 0 and u k # k #v k # for all k.
Proposition 4.1. Suppose that Assumption 2.1 holds and that f and c are twice
continuously di#erentiable near - x - . Suppose also that the sequence {z k } generated by
Algorithm A - converges to - z - and that the positive definite matrices M k satisfy the
estimate
(d x
when k #. Then the su#cient decrease condition (2.5) is satisfied with #
k su#ciently large provided that # < 1
A BFGS INTERIOR POINT ALGORITHM 211
Proof. Observe first that the positive definiteness of -
- with (4.1) implies that
(d x
for some positive constant K # and su#ciently large k. Observe also that d k # 0
(for d x
use (1.5), (4.2), and # - Therefore, for k large enough, z k
and z k are near - z - and one can expand - (z k . A second order
expansion gives for the left-hand side of (2.5)
We want to show that this quantity is negative for k large.
Our first aim is to show that # - (z k smaller than a
term of order o(#d k # 2 ). For this purpose, one computes
. On the other hand, using
one gets from Lemma 2.6
k .
With these estimates, (4.1), and the fact that # 2
with Lemma 3.1 and the boundedness of {#c k }, (4.3) becomes
(d x
(d #
clear that the result will be proven if we show that, for some
positive constant K and k large, # - (z k ) # d k # -K#d k # 2 . To show this, we use the
212 P. ARMAND, J. CH. GILBERT, AND S. JAN-J -
EGOU
last expression of # - (z k ) # d k and an upper bound of |(d x
k |, obtained by the
Cauchy-Schwartz inequality:
(d x
(d #
k .
It follows that
(d x
(d #
k .
Therefore, using (4.2) and Lemma 3.1, one gets
for some positive constant K and k large.
Proposition 4.1 shows in particular that the function V, which was added to # -
to get the merit function - , has the right curvature around - z - , so that the unit
step-size in both x and # is accepted by the line-search.
In the following proposition, we establish a necessary and su#cient condition of
q-superlinear convergence of the Dennis and Mor-e [13] type. The analysis assumes
that the unit step-size is taken and that the updated matrix M k is su#ciently good
asymptotically in a manner given by the estimate (4.5), which is slightly di#erent
from (4.1).
Proposition 4.2. Suppose that Assumption 2.1 holds and that f and c are
twice di#erentiable at - x - . Suppose that the sequence {z k } generated by Algorithm A -
converges to -
z - and that, for k su#ciently large, the unit step-size # accepted
by the line-search. Then {z k } converges q-superlinearly towards -
z - if and only if
Proof. Let us denote by M the nonsingular Jacobian matrix of the perturbed
KKT conditions (1.2) at the solution -
z -x -
# .
A first order expansion of the right-hand side of (1.3) about - z - and the identities
# -e give
Subtracting Md k from both sides and assuming a unit step-size, we obtain
z - #).
A BFGS INTERIOR POINT ALGORITHM 213
Suppose now that {z k } converges q-superlinearly. Then, the right-hand side
of (4.6) is of order o(#z k -
z - #), so that
Then (4.5) follows from the fact that, by the q-superlinear convergence of {z k }, z k -
Let us now prove the converse. By (4.5), the left-hand side of (4.6) is an o(#d k #)
and due to the nonsingularity of M, (4.6) gives z k+1 - z
With a unit step-size, d
z - (z k - z - ), so that we finally get z k+1 - z
For proving the q-superlinear convergence of the sequence {z k }, we need the
following result from the BFGS theory (see [40, Theorem 3] and [6]).
Lemma 4.3. Let {M k } be a sequence of matrices generated by the BFGS formula
from a given symmetric positive definite matrix M 1 and pairs (# k , # k ) of vectors
verifying
< #,
where M is a symmetric positive definite matrix. Then, the sequences {M k } and
k } are bounded and
By using this lemma, we will see that the BFGS formula gives the estimate
#).
Note that the above estimate implies (4.5), from which the q-superlinear convergence
of {z k } will follow.
A function #, twice di#erentiable in a neighborhood of a point x # R n , is said to
have a locally radially Lipschitzian Hessian at x, if there exists a positive constant L
such that for x # near x, one has
Theorem 4.4. Suppose that Assumption 2.1 holds and that f and c are C 1,1
functions, twice continuously di#erentiable near -
x - with locally radially Lipschitzian
Hessians at -
x - . Suppose that the line-search in Algorithm A - uses the constant # <2 . Then the sequence {z k generated by this algorithm converges to
z -x - ) q-superlinearly and, for k su#ciently large, the unit step-size #
accepted by the line-search.
Proof. Let us start by showing that Lemma 4.3 with
M - can be applied.
First, # k # k > 0, as this was already discussed after Lemma 3.3. For the convergence
of the series in (4.7), we use a Taylor expansion, assuming that k is large enough (f
and c are C 2 near -
214 P. ARMAND, J. CH. GILBERT, AND S. JAN-J -
EGOU
With the local radial Lipschitz continuity of # 2 f and # 2 c at - x - and the boundedness
of {# k+1 }, there exist positive constants K # and K # such that
z - # .
Hence the series in (4.7) converges by Theorem 3.4. Therefore, by (4.8) with
and the fact that # k is parallel to d x
#).
By the estimate (4.9) and Proposition 4.1, the unit step-size is accepted when
k is large enough. The q-superlinear convergence of {z k } follows from Proposition
4.2.
5. The overall primal-dual algorithm. In this section, we consider an overall
algorithm for solving problem (1.1). Recall from Lemma 2.2 that the set of primal
solutions of this problem is nonempty and bounded. By the Slater condition (As-
sumption 2.1(ii)), the set of dual solutions is also nonempty and bounded. Let us
denote by -
primal-dual solution of problem (1.1), which is also a solution
of the necessary and su#cient conditions of optimality
Our overall algorithm for solving (1.1) or (5.1), called Algorithm A, consists in
computing approximate solutions of the perturbed optimality conditions (1.2), for a
sequence of -'s converging to zero. For each -, the primal-dual Algorithm A - is used
to find an approximate solution of (1.2). This is done by so-called inner iterations.
Next - is decreased and the process of solving (1.2) for the new value of - is repeated.
We call an outer iteration the collection of inner iterations for solving (1.2) for a fixed
value of -. We index the outer iterations by superscripts j # N\{0}.
Algorithm A (for solving problem (1.1); one outer iteration).
At the beginning of the jth outer iteration, an approximation z j
# Z of the solution -
z of (5.1) is supposed available, as well as a positive
1 approximating the Hessian of the Lagrangian. A value
is given, as well as a precision threshold # j > 0.
1. Starting from z j
use Algorithm A - until z j :=
2. Choose a new starting iterate z j+1
for the next outer iteration, as
well as a positive definite matrix M j+1
1 . Set the new parameters - j+1 > 0
and # j+1 > 0, such that {- j
} and {# j
} converge to zero when j #.
A BFGS INTERIOR POINT ALGORITHM 215
To start the (j+1)th outer iteration, a possibility is to take z j+1
, the updated matrix obtained at the end of the jth outer iteration.
As far as the global convergence is concerned, how z are determined
is not important. Therefore, on that point, Algorithm A leaves the user much freedom
to maneuver, while Theorem 5.1 gives us a global convergence result for such a general
algorithm.
Theorem 5.1. Suppose that Assumption 2.1 holds and that f and c are C 1,1
functions. Then Algorithm A generates a bounded sequence {z j
} and any limit point
of {z j
} is a primal-dual solution of problem (1.1).
Proof. By Theorem 3.4, any outer iteration of Algorithm A terminates with an
iterate z j satisfying the stopping criteria in Step 1. Therefore Algorithm A generates
a sequence {z j
}. Since the sequences {- j
} and {# j
} converge to zero, any limit point
of {z j
} is a solution of problem (1.1). It remains to show that {z j
} is bounded.
Let us first prove the boundedness of {x j
}. The convexity of the Lagrangian
implies that
Using the positivity of # j and c(x 1 ) and next the stopping criteria of Algorithm A, it
follows that
#).
If {x j
} is unbounded, setting t j := #x j
# and y j := x j
one can choose a
subsequence J such that
lim
From the last inequality we deduce that
Moreover, since c(x j ) > 0, we have (-c (i) It follows that
(see, for example, [22, Proposition IV.3.2.5] or
[2, Formula (1)]). Therefore, the solution set of problem (1.1) would be unbounded,
which is in contradiction with what is claimed in Lemma 2.2.
To prove the boundedness of the multipliers, suppose that the algorithm generates
an unbounded sequence of positive vectors {# j
subsequence J # . The
sequence
#)} j#J # is bounded and thus has at least one limit point, say,
Dividing the two inequalities in (5.2) by # j
# and taking limits when j #,
we deduce that # 0, #c(x Using the concavity
of the components c (i) , one has
where the inequality on the right follows from the strict feasibility of the first iterate.
Multiplying by # , we deduce that (# c(x 1 contradiction
with
P. ARMAND, J. CH. GILBERT, AND S. JAN-J -
EGOU
In the rest of this section, we give conditions under which the whole sequence
} converges to a particular point called the analytic center of the primal-dual
optimal set. This actually occurs when the following two conditions hold: strict
complementarity and a proper choice of the forcing sequence # j in Algorithm A,
which has to satisfy the estimate
meaning that # j /- j
Let us first recall the notion of analytic center of the optimal sets, which under
Assumption 2.1 is uniquely defined (see Monteiro and Zhou [37], for related results).
We denote by opt(P ) and opt(D) the sets of primal and dual solutions of problem
(1.1). The analytic center of opt(P ) is defined as follows. If opt(P ) is reduced to
a single point, its analytic center is precisely that point. Otherwise, opt(P ) is a
convex set with more than one point. In that case, f is not strongly convex and,
by Assumption 2.1(i), at least one of the constraint functions, -c (i 0 ) say, is strongly
convex. It follows that the index set
is nonempty (it contains i 0 ). The analytic center of opt(P ) is then defined as the
unique solution of the following problem:
log c (i) (-x) # .
The fact that this problem is well defined and has a unique solution is the matter of
Lemma 5.2 below. Similarly, if opt(D) is reduced to a single point, its analytic center
is that point. In case of multiple dual solutions, the index set
is nonempty (otherwise opt(D) would be reduced to {0}). The analytic center of
opt(D) is then defined as the unique solution of the following problem:
log - # (i)
# .
Lemma 5.2. Suppose that Assumption 2.1 holds. If opt(P ) (resp., opt(D)) is not
reduced to a singleton, then problem (5.3) (resp., (5.4)) has a unique solution.
Proof. Consider first problem (5.3) and suppose that opt(P ) is not a singleton. We
have seen that B is nonempty. By the convexity of the set opt(P ) and the concavity
of the functions c (i) , there exists -
Therefore the
feasible set in (5.3) is nonempty. On the other hand, let -
x 0 be a point satisfying the
constraints in (5.3). Then the set
log c i (-x) # i#B
log c i (-x
is nonempty, bounded (Lemma 2.2), and closed. Therefore, problem (5.3) has a
solution. Finally, by Assumption 2.1(i), we know that there is an index i
A BFGS INTERIOR POINT ALGORITHM 217
that -c (i 0 ) is strongly convex. It follows that the objective in (5.3) is strongly concave
and that problem (5.3) has a unique solution.
By similar arguments and the fact that the objective function in (5.4) is strictly
concave, it follows that problem (5.4) has a unique solution.
By complementarity (i.e., C(-x) - and convexity of problem (1.1), the index
sets B and N do not intersect, but there may be indices that are neither in B nor
in N . It is said that problem (1.1) has the strict complementarity property if
{1, . , n}. This is equivalent to the existence of a primal-dual solution satisfying strict
complementarity.
Theorem 5.3. Suppose that Assumption 2.1 holds and that f and c are C 1,1
functions. Suppose also that problem (1.1) has the strict complementarity property
and that the sequence {# j
} in Algorithm A satisfies the estimate # Then the
sequence {z j
} generated by Algorithm A converges to the point -
z
x 0 is the analytic center of the primal optimal set and -
# 0 is the analytic center of the
dual optimal set.
Proof. Let (- x, -
#) be an arbitrary primal-dual solution of (1.1). Then -
x minimizes
#) and - # so that
Using the convexity of # j ) and the stopping criterion (5.2) of the inner iterations
in Algorithm A, one has
x#,
because
. By Theorem 5.1, there
is a constant C 1 such that m 1
adding the corresponding sides
of the two inequalities above leads to
We pursue this by adapting an idea used by McLinden [34] to give properties of
the limit points of the path -x - ). Let us define #
has for all indices i
c (i)
c (i)
Substituting this in (5.5) and dividing by - j give
c (i) (-x)
c (i)
By assumptions, #
Now supposing that (-x 0 , - # 0 ) is a
limit point of {(x j , # j )} and taking the limit in the preceding estimate provide
c (i) (-x)
P. ARMAND, J. CH. GILBERT, AND S. JAN-J -
EGOU
Necessarily, now that, by strict complementarity,
there are exactly m terms on the left-hand side of the preceding inequality. Hence,
by the arithmetic-geometric mean inequality
c (i) (-x)
or
c (i) (-x) #
One can take -
in this inequality, so that
c (i) (-x) #
c (i) (-x 0 ) and #
This shows that - x 0 is a solution of (5.3) and that - # 0 is a solution of (5.4). Since the
problems in (5.3) and (5.4) have unique solutions, all the sequence {x j
} converges to
x 0 and all the sequence {# j
} converges to - # 0 .
6. Discussion. By way of conclusion, we discuss the results obtained in this
paper, give some remarks, and raise some open questions.
Problems with linear constraints. The algorithm is presented with convex
inequality constraints only, but it can also be used when linear constraints are present.
Consider the problem
min f(x),
obtained by adding linear constraints to problem (1.1). In (6.1), A is a p - n matrix
with p < n and b # R p is given in the range space of A.
Problem (6.1) can be reduced to problem (1.1) by using a basis of the null space
of the matrix A. Indeed, let x 1 be the first iterate, which is supposed to be strictly
feasible in the sense that
Let us denote by Z an n - q matrix whose columns form a basis of the null space
of A. Then, any point satisfying the linear constraints of (6.1) can be written
With this notation, problem (6.1) can be rewritten as the problem in u # R
which has the form (1.1).
Thanks to this transformation, we can deduce from Assumption 2.1 what are
the minimal assumptions under which our algorithm for solving problem (6.2) or,
equivalently, problem (6.1) will converge.
A BFGS INTERIOR POINT ALGORITHM 219
Assumption 6.1. (i) The real-valued functions f and -c (i) (1 # i # m) are convex
and di#erentiable on the a#ne subspace X b} and at least one of the
functions f , -c (1) , . , -c (m) is strongly convex on X. (ii) There exists an x # R n
such that
With these assumptions, all the previous results apply. In particular, Algorithm
A - converges r-linearly (if f and c are also C 1,1 ) and q-superlinearly (if f
and c are also C 1,1 , twice continuously di#erentiable near - x - with locally radially
Lipschitzian Hessian at -
Similarly, the conclusions of Theorem 5.1 apply if f and
c are also C 1,1 .
Feasible algorithms and qN techniques. In the framework of qN methods,
the property of having to generate feasible iterates should not be only viewed as a
restriction limiting the applicability of a feasible algorithm. Indeed, in the case of
problem (6.2), if it is sometimes di#cult to find a strictly feasible initial iterate, the
matrix to update for solving this problem is of order q only, instead of order n for
an infeasible algorithm solving problem (6.1) directly. When q # n, the qN updates
will approach the reduced Hessian of the Lagrangian Z # 2 #)Z more rapidly than
the full Hessian # 2 #, so that a feasible algorithm is likely to converge more rapidly.
About the strong convexity hypothesis. Another issue concerns the extension
of the present theory to convex problems, without the strong convexity assumption
(Assumption 2.1(i)).
this hypothesis, the class of problems to consider encompasses linear
programming (f and c are a#ne). It is clear that for dealing properly with linear
programs, our algorithm needs modifications, since then # and the BFGS formula
is no longer defined. Of course, it would be very ine#ective to solve linear programs
with the qN techniques proposed in this paper (M is the desired matrix), but
problems that are almost linear near the solution may be encountered, so that a
technique for dealing with a situation where # k # k # can be of interest.
To accept # can look at the limit of the BFGS formula (2.1) when
possible update formula could be
The updated matrix satisfies M k+1 # positive semidefinite, provided
is already positive semidefinite. The fact that M k+1 may be singular raises some
di#culties, however. For example, the search direction d x may no longer be defined
(see formula (1.5), in which the matrix M +#c(x)C(x) -1 #c(x) # can be singular).
Therefore, the present theory cannot be extended in a straightforward manner.
On the other hand, the strong convexity assumption may not be viewed as an important
restriction, because a fictive strongly convex constraint can always be added.
An obvious example of fictive constraint is "x # x # K." If the constant K is large
enough, the constraint is inactive at the solution, so that the solution of the original
problem is not altered by this new constraint and the present theory applies.
Better control of the outer iterations. Last but not least, the global convergence
result of section 5 is independent of the update rule of the parameters # j and
In practice, however, the choice of the decreasing values # j and - j is essential for
the e#ciency of the algorithm and would deserve a detailed numerical study.
From a theoretical viewpoint, it would be highly desirable to have an update rule
that would allow the outer iterates of Algorithm A to converge q-superlinearly. Along
P. ARMAND, J. CH. GILBERT, AND S. JAN-J -
EGOU
the same lines, an interesting problem is to design an algorithm in which the barrier
parameter would be updated at every step, while having q-superlinear convergence of
the iterates. Such extensions would involve more di#cult issues.
The global convergence result proved in this paper gives us some reasons to believe
that it is not unreasonable to tackle these open questions.

Acknowledgments

. We would like to thank the referees for their valuable com-
ments. One of them has shown us a direct argument for the last part of the proof of
Proposition 3.2, which is the one we have finally chosen to give in the paper. The other
referee has brought McLinden's paper to our attention, which led us to Theorem 5.3.



--R

On the convergence of an infeasible primal-dual interior-point method for convex programming
Asymptotic analysis for penalty and barrier methods in convex and linear programming
Computational Di
A trust region interior point algorithm for linearly constrained optimization
A. Trust Region Method Based on Interior Point Techniques for
A tool for the analysis of quasi-Newton methods with application to unconstrained minimization
A Primal-Dual Algorithm for Minimizing a Non-convex Function Subject to Bound and Linear Equality Constraints
Principles and Techniques of Algorithmic Di

Interior Point Approach to Linear
A Potential Reduction Method for a Class of Smooth Convex Programming Problems
On the classical logarithmic barrier function method for a class of smooth convex programming problems

Numerical Methods for Unconstrained Optimization and Nonlinear Equations
On the formulation and theory of the Newton interior-point method for nonlinear programming
Sequential Unconstrained Minimization Techniques
Practical Methods of Optimization



Feasible direction interior-point technique for nonlinear optimization

Interior Point Techniques in Optimization-Complementarity
On the method of analytic centers for solving smooth convex problems

A practical interior-point method for convex programming


A Unified Approach to Interior Point Algorithms for Linear Complementarity Problems
A new continuation method for complementarity problems with uniform P
Limiting behavior of trajectories generated by a continuation method for monotone complementarity problems
On the limited memory BFGS method for large scale optimization
The projective SUMT method for convex programming problems
An analogue of Moreau's proximation theorem
An interior point algorithm for solving smooth convex programs based on Newton's method
An extension of Karmarkar-type algorithms to a class of convex separable programming problems with global linear rate of convergence
On the existence and convergence of the central path for convex programming and some duality results

Updating quasi-Newton matrices with limited storage
On the convergence of the variable metric algorithm
Some global convergence properties of a variable metric algorithm for minimization without exact line searches
Theory and Algorithms for Linear Optimization-An Interior Point Approach
"analytical center"
Interior Point Methods of Mathematical Programming
Computational experience with a primal-dual interior-point method for smooth convex programming
Why a pure primal Newton barrier step may be infeasible

Superlinear and quadratic convergence of some primal-dual interior point methods for constrained optimization
Interior Point Algorithms-Theory and Analysis
in Integer Nonlinear Pro- gramming
--TR

--CTR
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
Paul Armand, A Quasi-Newton Penalty Barrier Method for Convex Minimization Problems, Computational Optimization and Applications, v.26 n.1, p.5-34, October
Dingguo Pu , Weiwen Tian, The revised DFP algorithm without exact line search, Journal of Computational and Applied Mathematics, v.154 n.2, p.319-339, 15 May
convex programming;superlinear convergence;constrained optimization;interior point algorithm;analytic center;BFGS quasi-Newton approximations;primal-dual method;line-search