AN INTERIOR POINT ALGORITHM FOR LINEARLY CONSTRAINED OPTIMIZATION^{∗}

STEPHEN J. WRIGHT^{†}

SIAM J. OPTIMIZATION. \begin{matrix}c \\ ◯ \end{matrix} 1992 Society for Industrial and Applied Mathematics Vol. 1, No. 4, pp. 000–000, Month 1992 012

Abstract. We describe an algorithm for optimization of a smooth function subject to general linear constraints. An algorithm of the gradient projection class is used, with the important feature that the “projection” at each iteration is performed using a primal-dual interior point method for convex quadratic programming. Convergence properties can be maintained even if the projection is done inexactly in a well-defined way. Higher-order derivative information on the manifold defined by the apparently active constraints can be used to increase the rate of local convergence.

Key words. potential reduction algorithm, gradient porojection algorithm, linearly constrained optimization

AMS(MOS) subject classifications. 65K10, 90C30

1. Introduction. We address the problem

min_{x} f(x)  s.t.  A^{T} x ≤ b, (1)

where x ∈ R^{n} and b ∈ R^{m}, and f is assumed throughout to be twice continuously differentiable on the level set

L = {x | A^{T} x ≤ b, f(x) ≤ f(x^{0})},

where x^{0} is some given initial choice for x. Recent literature on this problem can for the most part be divided into two main classes. On the one hand, there are the “active set” approaches such as sequential quadratic programming, which are most suitable when the constraints A^{T} x ≤ b lack any special structure such as separability. In these algorithms a model of f (for example, the quadratic approximation f(x) + ∇f(x)^{T} d + (1/2)d^{T} ∇^{2}f(x)d) is formed at each “outer” iteration and minimized over some subset of the feasible region. The algorithm tends to move along edges and faces of the boundary of the feasible set, changing its set of currently active constraints by at most one element on each “inner” iteration. A second class of methods, known as “gradient projection” methods, allow more substantial changes to the active set at each iteration by choosing a direction g (for example, ∇f(x) or some scaled version of it) and searching along the piecewise linear path P (x − αg), where α > 0 and P is the projection onto the feasible set. These methods are best suited to the case in which the projection P (.) is easy to perform, for example, when the feasible region is a box whose sides are parallel to the principal coordinate axes.

In this paper, our aim is to describe an algorithm of the gradient projection class, in which we allow the projections to be performed inexactly. We focus on the case of Euclidean norm projections, which can be solved by using interior-point methods for convex quadratic programming problems. In this way, general polyhedral feasible regions can be handled. We thus hope to combine the much-vaunted advantages of

∗ Received by the editors July xx, 1990; accepted for publication (in revised form) September xx, 1991. This research was partly performed at the University of Queensland, Australia, whose support is gratefully acknowledged. The work was also supported by the Applied Mathematical Sciences subprogram of the Office of Energy Research, U. S. Department of Energy, under Contract W-31-109-Eng-38.

† Argonne National Laboratory, 9700 South Cass Avenue, Argonne, Illinois, 60439.

interior-point methods with the desirable properties of gradient projection algorithms — most notably, rapid identification of the final active set. In addition, we allow second-derivative information to be used in the definition of g (as is also done by Dunn [4, 3] and Gafni and Bertsekas [5]) to speed up the asymptotic convergence rate after the correct active set has been identified.

The “inexactness” in the projection is quantified by a duality gap, which is updated at each iteration of the projection subproblem. The global convergence analysis in section 4 is not tied to the use of an interior-point method for the projection; any algorithm (including an active set method) that allows a duality gap to be calculated for each iterate may be used.

The point x^{∗} is a critical point for (1) if there are scalars y_{i} ≥ 0 such that

−∇f(x^{∗}) = ∑_{i∈A} y_{i}a_{i},

where a_{i} are columns of A, and

A = {i = 1, ⋯, m | a_{i}^{T} x^{∗} = b_{i}}.

Equivalently,

− ∇f(x^{∗}) ∈ N(x^{∗}; X), (2)

where X is the feasible set {x | A^{T} x ≤ b}, and N(x; X) is the normal cone to X at x, defined by

N(x; X) = {v | v^{T} (u − x) ≤ 0, for all u ∈ X}.

In the next section, we specify the algorithm. The interior-point method that may be used to perform the projection is discussed in section 3. The global and local convergence properties of the algorithm are analyzed in section 4 and section 5, respectively.

In the remainder of the paper, the following notational conventions will be used: • ∥x∥ = (x^{T} x)^{\frac{1}{2}} (the Euclidean norm), unless otherwise specified.

• P_{Y} (x) denotes the Euclidean projection of the vector x onto the convex set Y ⊂ R^{n}, that is

P_{Y} (x) = arg min_{z∈Y} ∥z − x∥.

If the subscript is omitted from P, projection onto X is assumed.

• intY denotes the interior of Y, and ∂Y denotes its boundary.

• When x is a vector, relations such as x > 0 are meant to apply componentwise.

• Subscripts on vectors and matrices denote components, while superscripts are used to distinguish different iterates. Subscripts on scalars denote iteration numbers.

• When {ξ_{k}} and {ξ̄_{k}} are non-negative sequences, the notation ξ_{k} = O(ξ̄_{k}) means that there is a constant s such that ξ_{k} ≤ sξ̄_{k} for all k sufficiently large. ξ_{k} = o(ξ̄_{k}) means that there is a non-negative sequence {s_{k}} converging to zero such that ξ_{k} ≤ s_{k}ξ̄_{k} for all k sufficiently large.

• The sequence {v^{k}} is said to converge Q-quadratically to v^{∗} if ∥v^{k+1} − v^{∗}∥ = O(∥v^{k} − v^{∗}∥^{2}). It is said to converge R-quadratically if there is a sequence {ξ_{k}} that converges Q-quadratically to zero such that ∥v^{k} − v^{∗}∥ ≤ ξ_{k} for all k.

• If {v^{k}} and {v̄^{k}} are two sequences of vectors, the notation “v^{k} → v̄^{k}” means that lim_{k→∞} ∥v^{k} − v̄^{k}∥ = 0.

• In sections 4 and 5, we introduce constants denoted by C and C̄ with a subscript. In all cases these represent strictly positive constants, even where not stated explicitly.

2. The Algorithm. We start this section by giving an outline of the major operations at each iteration of the basic algorithm. Then we state a formal outline and conclude by mentioning possible variations.

The algorithm first defines an “almost active” set of constraints at each iterate x^{k}. It partitions the gradient into two orthogonal components (which are orthogonal to and tangent to the manifold defined by the almost active set, respectively) and then scales the tangent component by a matrix with suitable positive definiteness properties (possibly an inverse reduced Hessian or a quasi-Newton approximation to it). A projected Armijo-like line search is then performed along the resulting direction.

The activity tolerance at the point x^{k} is ϵ_{k}, where for the moment we require only that ϵ_{k} ≥ 0. The almost active set I^{k} is defined by

I^{k} = {i = 1, ⋯, m | a_{i}^{T} x^{k} ≥ b_{i} − ϵ_{k}∥a_{i}∥}. (3)

We use T^{k} to denote the tangent manifold corresponding to this set:

T^{k} = {z | a_{i}^{T} z = 0, all i ∈ I^{k}}. (4)

The negative gradient is then decomposed using T^{k} by setting

d^{k} = P_{T^{k}}(−∇f(x^{k})), d^{k+} = −[∇f(x^{k}) + d^{k}]. (5)

The tangent component d^{k} is modified by setting

d̃^{k} = D^{k}d^{k}, (6)

where D^{k} is a matrix such that P_{T^{k}} ∘ D^{k} ∘ P_{T^{k}} = D^{k} and

λ_{1}z^{T} z ≤ z^{T} D^{k}z ≤ λ_{2}z^{T} z, all z ∈ T^{k} , (7)

where λ_{1} and λ_{2} are positive constants. The search direction is assembled as

g^{k} = −(d̃^{k} + d^{k+}). (8)

A projected Armijo search is carried out along the path

x^{k}(α) = P (x^{k} − αg^{k}),

where the values α = 1, β, β^{2}, β^{3}, ⋯ (where β ∈ (0, 1) is some constant) are tried. For each such value of α, the projection is calculated with the algorithm described in the next section. This algorithm generates a sequence of feasible approximations to x^{k}(α), which we denote by x^{kj}(α). For each such estimate, the algorithm produces a duality gap γ_{kj}(α). Defining the more convenient quantity

δ_{kj}(α) = \sqrt{2γ_{kj}(α)},

we can obtain upper and lower bounds on the distance from x^{k} − αg^{k} to X, that is

∥x^{kj}(α) − (x^{k} − αg^{k})∥^{2} − δ_{kj}(α)^{2} ≤ ∥x^{k}(α) − (x^{k} − αg^{k})∥^{2} ≤ ∥x^{kj}(α) − (x^{k} − αg^{k})∥^{2}.

These “inner iterations” are stopped at a value of j for which δ_{kj}(α) becomes sufficiently small according to the following criteria:

δ_{kj}(α) ≤ ηα^{τ/2} max (\frac{∥x^{kj}(α) − (x^{k} + αd̃^{k})∥}{α}, ∥d^{k}∥)^{2} (9)

and

δ_{kj}(α) ≤ C_{1}α^{τ/2}. (10)

Here τ, C_{1}, and η are constants that satisfy the conditions

τ > 2, ηC_{1} < 1.

We denote the final computed δ_{kj}(α) by δ_{k}(α), and the corresponding x^{kj}(α) by x^{k}(α; δ_{k}(α)). The step α is then accepted if the following “sufficient decrease” test is satisfied:

f(x^{k}) − f(x^{k}(α; δ_{k}(α))) ≥ σ \{αd^{kT} D^{k}d^{k} + \frac{∥x^{k}(α; δ_{k}(α)) − (x^{k} + αd̃^{k})∥^{2}}{α} \} , (11)

where σ ∈ (0, 1) is a constant.

The algorithm can be summarized as follows:

Step 1: Choose ϵ_{k}. Compute I^{k} from (3), and g^{k} according to (5)–(8).

Step 2: For α = β^{p}, p = 0, 1, 2, ⋯ (in sequence) approximately calculate x^{k}(α) = P (x^{k} − αg^{k}), terminating when x^{k}(α; δ_{k}(α)) = x^{kj}(α) and its associated δ_{k}(α) = δ_{kj}(α) are found that satisfy (9),(10). If the test (11) is passed for this value of α, set α_{k} = α = β^{p}, x^{k+1} = x^{k}(α; δ_{k}(α)), k ← k + 1, and go to the next iteration. Otherwise, increase p by 1, and try the next α = β^{p}.

In its “exact” form (i.e., δ_{k}(α) ≡ 0), and when D^{k} is defined as the reduced Hessian or a quasi-Newton approximation to it, the step g^{k} is the same as that obtained by specializing the algorithm of Dunn [4] to the linearly constrained case. The calculation of g^{k} is somewhat different in Gafni and Bertsekas [5]. They define an “almost tangent cone” at x^{k} by

C^{k} = {z | a_{i}^{T} z ≤ 0, all i ∈ I^{k}},

and then define d^{k} as the projection of −∇f(x^{k}) onto this cone. Additionally, the conditions on D^{k} are slightly different, and d̃^{k} is the projection of D^{k}d^{k} onto C^{k}. Our reason for following Dunn [4] and using the simpler decomposition relative to T^{k} is our assumption that projection onto the subspace T^{k} can be done exactly and cheaply. This is not unreasonable — the cost would normally be comparable to one iteration of the interior-point algorithm used for the projection onto X. Projection onto C^{k} may, on the other hand, be as expensive as projection onto X. Still, there are intuitive reasons for preferring C^{k} to T^{k}, and it would be of interest to see whether the extra cost per iteration (and the extra algorithmic complexity of doing the projection onto C^{k} inexactly) is justified.

The steplength rule (11) reduces to the one proposed by Gafni and Bertsekas [5] (and also used by Dunn [3]) when δ_{k}(α) ≡ 0. Another obvious possibility, to which we will return briefly in section 5, is

f(x^{k}) − f(x^{k}(α; δ_{k}(α)) ≥ σ {αd^{kT} D^{k}d^{k} + ∇f(x^{k})^{T} [x^{k} + αd̃^{k} − x^{k}(α; δ_{k}(α))]} . (12)

3. Projection onto X. Projection onto the polyhedral set X can be achieved by solving a convex quadratic program or, equivalently, a linear complementarity problem (LCP). In this section, we formulate the problem and outline a primal-dual potential reduction algorithm for solving it. The discussion will be brief, since other papers such as [6, 7, 10, 11] can be consulted for details about motivation, analysis, and implementation issues for this class of interior-point algorithms.

Throughout the remainder of the paper, we use the following assumptions:

(A) The feasible set X has an interior in R^{n}.

(B) At the solution z^{∗} = P (t) of the projection subproblem, the set of vectors

{a_{i} | a_{i}^{T} z^{∗} = b_{i}}

is linearly independent.

The (unique) vector P (t) is obtained by solving

min \frac{1}{2}∥z − t∥^{2} s.t. A^{T} z ≤ b,

or, equivalently,

min \frac{1}{2}∥z − t∥^{2} s.t. A^{T} z + ν = b, ν ≥ 0. (13)

Introducing Lagrange multipliers y for the constraints, we find that (13) is equivalent to the (mixed) LCP

[ \begin{matrix}0 \\ ν \end{matrix} ] = [ \begin{matrix}I & A \\ −A^{T} & 0 \end{matrix} ] [ \begin{matrix}z \\ y \end{matrix} ] + [ \begin{matrix}−t \\ b \end{matrix} ] , ν ≥ 0, y ≥ 0, ν^{T} y = 0. (14)

The coefficient matrix in (14) is clearly positive semi-definite.

The progress of the interior-point algorithm can be gauged by using the potential function defined by

ψ(ν, y) = ρ_{P} log(ν^{T} y) − ∑_{i=1}^{m} log(ν_{i}y_{i}), (15)

where ρ_{P} ≥ m + \sqrt{m}. In Kojima, Mizuno, and Yoshise [7], the step from iterate j to iterate j + 1 is obtained by solving the linear system

[ \begin{matrix}0 \\ Δν \end{matrix} ] = [ \begin{matrix}I & A \\ −A^{T} & 0 \end{matrix} ] [ \begin{matrix}Δz \\ Δy \end{matrix} ] , (16)

together with

ν_{i}^{j}Δy_{i}^{j} + y_{i}^{j}Δν_{i} = (\frac{−ρ_{j}}{1 − ρ_{j}} ) [−ν_{i}^{j}y_{i}^{j} + \frac{γ_{j}}{ρ_{j}} ] , i = 1, ⋯, m, (17)

where γ_{j} = ∑_{i=1}^{m} ν_{i}^{j}y_{i}^{j} and the value of ρ_{j} is discussed below. A steplength θ_{j} is chosen such that

θ_{j} |\frac{Δν_{i}}{ν_{i}^{j}} | ≤ τ, θ_{j} |\frac{Δν_{i}}{ν_{i}^{j}} | ≤ τ, i = 1, ⋯, m, (18)

for some τ ∈ (0, 1). Trivial modifications of the results of Kojima, Mizuno, and Yoshise [7] indicate that for the choices ρ_{j} ≡ ρ_{P} = m+\sqrt{m} and τ = 0.4, we have that

ψ(ν^{j} + θ_{j}Δν, y^{j} + θ_{j}Δy) ≤ ψ(ν^{j}, y^{j}) − 0.2. (19)

When some iterate (z^{j}, ν^{j}, y^{j}) satisfies ψ(ν^{j}, y^{j}) ≤ −O(\sqrt{m}L), it can easily be shown that (ν^{j})^{T} y^{j} ≤ 2^{−O(L)}. This suggests that, provided the initial point (z^{0}, ν^{0}, y^{0}) satisfies ψ(ν^{0}, y^{0}) = O(\sqrt{m}L), convergence to a point with duality gap less than 2^{−O(L)} can be achieved in O(\sqrt{m}L) iterations. (For purposes of the complexity analysis, L is taken to be the “size” of the problem.) Although this choice of ρ_{j} yields the best complexity result to date, it has been observed that, in practice, larger values of ρ_{j} lead to fewer iterations. In Han, Pardalos, and Ye [6], the choice ρ_{j} ≡ m^{2} is made for convex quadratic programs. In the context of linear programs, Zhang, Dennis, and Tapia [11] observe that it is even desirable to let ρ_{j} grow unboundedly large as the solution is approached. (The steps produced by (16),(17) are then very close to being Newton steps for the nonlinear equations formed by the equalities in (14).) Ye et al. [10] have shown such “large” choices of ρ_{j} are not incompatible with obtaining reductions in the potential function. In practical implementations, the line search parameter θ_{j} is also chosen differently. In Han, Pardalos, and Ye [6], the following choice appeared to give good experimental results:

θ_{j} = 0.99 min (  min_{i=1,···,m, Δν_{i}<0} −\frac{ν_{i}^{j}}{Δν_{i}} , min_{i=1,···,m, Δy_{i}<0} −\frac{y_{i}^{j}}{Δy_{i}} ) .

An issue of particular concern in this context is the choice of a feasible initial point at which to start the interior-point iteration. Such a point can be found by augmenting the problem in a simple way. We can reasonably assume that a vector z^{0} that satisfies A^{T} z^{0} < b is available from some previous iteration. If y^{0} is also chosen from a previous iteration, we usually have, from the first equation in (14), that z^{0} + Ay^{0} is similar in magnitude to the primal quantities z and t. We can thus define a (reasonably scaled) vector q by

q = −(z^{0} − t + Ay^{0})

and obtain the following augmented version of (14):

[ \begin{matrix}0 \\ ν \\ ν_{m+1} \end{matrix} ] = [ \begin{matrix}I & A & q \\ −A^{T} & 0 & 0 \\ −q^{T} & 0 & 0 \end{matrix} ] [ \begin{matrix}z \\ y \\ y_{m+1} \end{matrix} ] + [ \begin{matrix}−t \\ b \\ b_{m+1} \end{matrix} ] , (20)

ν ≥ 0, ν_{m+1} ≥ 0, y ≥ 0, y_{m+1} ≥ 0, ν^{T} y + ν_{m+1}y_{m+1} = 0.

The corresponding projection problem is

min \frac{1}{2}∥z − t∥^{2} s.t. A^{T} z ≤ b, q^{T} z ≤ b_{m+1}. (21)

If we choose b_{m+1} to satisfy

b_{m+1} > max(q^{T} z^{0}, q^{T} P (t)),

then we find that a feasible initial point for (20) is

(z, ν, ν_{m+1}, y, y_{m+1}) = (z^{0}, b − A^{T} z^{0}, b_{m+1} − q^{T} z^{0}, y^{0}, 1).

At the optimal solution, ν_{m+1}^{∗} = b_{m+1} − q^{T} P (t) and y_{m+1}^{∗} = 0. A practical choice for b_{m+1} can be made as follows: when t = x − αg with x feasible, note that

q^{T} P (t) = q^{T} [P (x−αg)−(x−αg)]+q^{T} t ≤ ∥q∥∥P (x−αg)−(x−αg)∥+q^{T} t ≤ α∥q∥∥g∥+q^{T} t.

Hence b_{m+1} can be chosen as any number greater than

max(q^{T} z^{0}, α∥q∥∥g∥ + q^{T} t).

We tacitly assume throughout the remainder of the paper that b_{m+1} is chosen large enough that the extra constraint in (21) does not come into play during the projection process (that is, ν_{m+1} stays reasonable large).

Two more points about the computational aspects of the projection should be made since, for many variants of the algorithm described in this paper, it will be the most time-consuming step, apart from the function evaluations. First, note that the cost per interior-point iteration, which is dominated by the cost of solving augmented versions of the linear system (16),(17), is similar to the cost of decomposing the gradient as in (5). (The latter operation may be performed by solving a system whose coefficient matrix is a submatrix of the matrix in (16).) Second, the number of interior-point iterates which will be necessary for a given α should not be too large. A rule of thumb seems to be that around 20–30 iterates are required for an accurate solution when no a priori information about the solution is known. In our case, the situation is better: good starting points will usually be available from previous iterates and from approximate projections for larger values of α. A priori information has been observed to significantly decrease the number of interior-point iterations (see, for example [9]).

In section 5, we assume that the points (z^{j}, ν^{j}, ν_{m+1}^{j}, y^{j}, y_{m+1}^{j}) generated by the interior-point algorithm do not stray too far from the central path defined by

{(z, ν, ν_{m+1}, y, y_{m+1}) feasible in (20) | ν_{i}y_{i} = ∑_{l=1}^{m+1} ν_{l}y_{l}/(m + 1), i = 1, ⋯, m + 1} .

The following assumption is used to prove that unit steps α_{k} = 1 are always eventually used by the method.

(C) There is a constant μ > 1 such that the final iterate (z, ν, ν_{m+1}, y, y_{m+1}) generated by the projection algorithm, each time it is called, satisfies

ν_{i}y_{i} ≥ \frac{∑_{l=1}^{m+1} ν_{l}y_{l}/(m + 1)}{μ}.

Although this assumption conflicts to some extent with the desire for fast asymptotic convergence of the interior-point method, Zhang, Dennis, and Tapia [11, Theorem 3.1] observed that, at least in the case of linear programming that they consider, it appears to hold in practice.

4. Global Convergence. In this section we prove that all accumulation points of the algorithm of section 2 are critical. The result depends crucially on the following lemma, which bounds the distance between x^{k}(α; 0) and x^{k}(α; δ_{k}(α)) in terms of δ_{k}(α).

Lemma 4.1. Suppose that (A) holds and that (B) holds at z^{∗} = x^{k}(α, 0). Then

∥x^{k}(α; 0) − x^{k}(α; δ_{k}(α))∥ ≤ δ_{k}(α).

Proof. Setting t = x^{k} − αg^{k}, we obtain

∥t − x^{k}(α; 0)∥^{2} ≥ ∥t − x^{k}(α; δ_{k}(α))∥^{2} − δ_{k}(α)^{2}
⇒ δ_{k}(α)^{2} ≥ 2[t − x^{k}(α; 0)]^{T} [x^{k}(α; 0) − x^{k}(α; δ_{k}(α))] + ∥x^{k}(α; 0) − x^{k}(α; δ_{k}(α))∥^{2}.

Now since t − x^{k}(α; 0) ∈ N(x^{k}(α; 0) ; X) and x^{k}(α; δ_{k}(α)) ∈ X, the first term on the right-hand side above is non-negative and can be omitted from the inequality. The result follows.

Under appropriate nondegeneracy assumptions, application of the implicit function theorem to a subset of the equalities in (13) (or (20)) would suggest that, locally, a stronger bound of O(γ_{k}(α)) = O(δ_{k}(α)^{2}) might be obtained. In fact, some of the local convergence analysis in section 5 relies on just this observation. In general, however, given a point x^{k} and a search direction g^{k}, there are usually values of α such that the solution of (13) (or (21)) for t = x^{k} − αg^{k} is degenerate. Our result in Lemma 4.1 is similar to, but more specific than, the bound that would be obtained by applying the analysis of Mangasarian and Shiau [8] to (13).

We state without proof the following well-known result, which actually applies for any closed convex X ⊂ R^{n}.

Lemma 4.2. For any x ∈ X and z ∈ R^{n},

a) ∥P (x + αz) − x∥/α is a nonincreasing function of α > 0,

b) ∥P (x + αz) − x∥/α ≤ ∥z∥.

Before proving the main result (Theorem 4.5), we show that the conditions (9),(10) ensure that the projection is computed exactly when x^{k} is critical (Lemma 4.3) and, in a technical result, show that the algorithm produces descent at a noncritical point (Lemma 4.4).

Lemma 4.3. Suppose that (A) holds and that (B) holds at z^{∗} = x^{k}. When x^{k} is critical, then δ_{k}(α) = 0 for all α ∈ [0, 1], and x^{k}(α; δ_{k}(α)) = x^{k} for all α ∈ [0, 1].

Proof. Clearly the result is true for α = 0. For the remainder of the proof, we assume that α ∈ (0, 1].

All vectors in the subspace T^{k} are orthogonal to N(x^{k}; X). Hence by (2) and (5), d^{k} = d̃^{k} = 0 and d^{k+} = −∇f(x^{k}). Also, by (2),

x^{k}(α, 0) = P (x^{k} − α∇f(x^{k})) = x^{k},

and so

∥x^{k}(α; δ_{k}(α)) − (x^{k} + αd̃^{k})∥ ≤ ∥x^{k}(α; 0) − x^{k}∥ + ∥x^{k}(α; δ_{k}(α)) − x^{k}(α; 0)∥ ≤ δ_{k}(α).

Substituting this expression in (9), we have

δ_{k}(α) ≤ ηα^{τ/2} \frac{δ_{k}(α)^{2}}{α^{2}}

and hence

δ_{k}(α)[1 − ηα^{τ/2−2}δ_{k}(α)] ≤ 0. (22)

From (10) and the fact that ηC_{1} < 1,

1 − ηα^{τ/2−2}δ_{k}(α) ≥ 1 − ηC_{1}α^{τ−2} > 1 − α^{τ−2} ≥ 0,

since α ∈ [0, 1] and τ > 2. Since δ_{k}(α) ≥ 0, the inequality (22) can hold only if δ_{k}(α) = 0. Thus, the first statement is proved. Proof of the second statement follows immediately.

Lemma 4.4. Suppose I^{k} is defined by (3), where ϵ_{k} is any positive number. Suppose that (A) holds and that (B) holds for all z^{∗} = x^{k}(α, 0) for α ∈ [0, 1]. Then, given any σ̄ ∈ (0, 1) there exists an ᾱ(σ̄) ∈ (0, ϵ_{k}/∥g^{k}∥) such that

∇f(x^{k})^{T} [x^{k} − x^{k}(α; δ_{k}(α))]
(23)  ≥ σ̄ [αd^{kT} D^{k}d^{k} + \frac{1}{α}∥x^{k}(α; δ_{k}(α)) − (x^{k} + αd̃^{k})∥^{2}]

Hence, provided x^{k} is not critical, there is an α̂(σ̄) ∈ (0, ᾱ] such that f(x^{k}) > f(x^{k}(α, δ_{k}(α)) for all α ∈ (0, α̂].

Proof.

∇f(x^{k})^{T} [x^{k} − x^{k}(α; δ_{k}(α))]
(24)  = ∇f(x^{k})^{T} [x^{k} − x^{k}(α; 0)] + ∇f(x^{k})^{T} [x^{k}(α; 0) − x^{k}(α; δ_{k}(α))]

and for α ∈ (0, ϵ_{k}/∥g^{k}∥), it can be proved by using a similar argument to that in [5, Proposition 1 (b)] that

∇f(x^{k})^{T} [x^{k} − x^{k}(α, 0)]
≥ αd^{kT} D^{k}d^{k} + \frac{1}{α}∥x^{k}(α, 0) − (x^{k} + αd̃^{k})∥^{2}.

By the smoothness assumptions on f, there is a constant B such that

∥∇f(x)∥ ≤ B  for all  x ∈ L.

Since all x^{k} ∈ L, we have, using Lemma 4.1, that

|∇f(x^{k})^{T} [x^{k}(α; 0) − x^{k}(α; δ_{k}(α))]| ≤ Bδ_{k}(α). (25)

Now

∥x^{k}(α; 0) − (x^{k} + αd̃^{k})∥^{2} = ∥x^{k}(α; δ_{k}(α)) − (x^{k} + αd̃^{k})∥^{2} (26)

+2[x^{k}(α; 0) − x^{k}(α; δ_{k}(α))]^{T} [x^{k}(α; δ_{k}(α)) − (x^{k} + αd̃^{k})] + ∥x^{k}(α; 0) − x^{k}(α; δ_{k}(α))∥^{2}

≥  ∥x^{k}(α; δ_{k}(α)) − (x^{k} + αd̃^{k})∥^{2} − 2δ_{k}(α)∥x^{k}(α; δ_{k}(α)) − (x^{k} + αd̃^{k})∥,

and so from (24)–(26)

∇f(x^{k})^{T} [x^{k} − x^{k}(α; δ^{k}(α))]
(27)  ≥ αd^{kT} D^{k}d^{k} + \frac{1}{α}∥x^{k}(α; δ_{k}(α)) − (x^{k} + αd̃^{k})∥^{2}
−\frac{2}{α}δ_{k}(α)∥x^{k}(α; δ_{k}(α)) − (x^{k} + αd̃^{k})∥ − Bδ_{k}(α).

Now,

\frac{1}{α}∥x^{k}(α; δ_{k}(α)) − (x^{k} + αd̃^{k})∥
(28)  ≤ \frac{1}{α}∥x^{k}(α; δ_{k}(α)) − x^{k}(α; 0)∥ + \frac{1}{α}∥x^{k}(α; 0) − (x^{k} + αd̃^{k})∥
≤  \frac{1}{α}δ_{k}(α) + \frac{1}{α}∥P (x^{k} + αd^{k+} + αd̃^{k}) − (x^{k} + αd̃^{k})∥.

The following simple argument shows that x^{k} + αd̃^{k} ∈ X for α ∈ (0, ϵ_{k}/∥g^{k}∥):

i ∉ I^{k} ⇒ a_{i}^{T} [x^{k} + αd̃^{k}] ≤ b_{i} − ϵ_{k}∥a_{i}∥ + α∥g^{k}∥∥a_{i}∥ < b_{i}
i ∈ I^{k} ⇒ a_{i}^{T} [x^{k} + αd̃^{k}] = a_{i}^{T} x^{k} ≤ b_{i}.

Hence by Lemma 4.2(b), (28) becomes

\frac{1}{α}∥x^{k}(α; δ_{k}(α)) − (x^{k} + αd̃^{k})∥ ≤ \frac{1}{α}δ_{k}(α) + ∥d^{k+}∥ ≤ \frac{1}{α}δ_{k}(α) + B. (29)

Hence (27) becomes

∇f(x^{k})^{T} [x^{k} − x^{k}(α; δ_{k}(α))] (30)

≥ αd^{kT} D^{k}d^{k} + \frac{1}{α}∥x^{k}(α; δ_{k}(α)) − (x^{k} + αd̃^{k})∥^{2} − \frac{2}{α}δ_{k}(α)^{2} − 3Bδ_{k}(α).

We now consider two cases. First, suppose that

\frac{1}{α}∥x^{k}(α; δ_{k}(α)) − (x^{k} + αd̃^{k})∥ ≥ ∥d^{k}∥.

Then from (9) it follows that

δ_{k}(α) ≤ ηα^{τ/2} \frac{∥x^{k}(α; δ_{k}(α)) − (x^{k} + αd̃^{k})∥^{2}}{α^{2}} .

Using this, together with (10) and the fact that ηC_{1} < 1, we have from (30) that

∇f(x^{k})^{T} [x^{k} − x^{k}(α; δ_{k}(α))]
(31) ≥  αd^{kT} D^{k}d^{k} + \frac{1}{α}∥x^{k}(α; δ_{k}(α)) − (x^{k} + αd̃^{k})∥^{2}
− [\frac{2}{α}(C_{1}α^{τ/2})ηα^{τ/2} + 3Bηα^{τ/2}] (\frac{1}{α^{2}} )∥x^{k}(α; δ_{k}(α)) − (x^{k} + αd̃^{k})∥^{2}
≥ αd^{kT} D^{k}d^{k} + \frac{1}{α} [1 − 2α^{τ−2} − 3Bηα^{\frac{τ−2}{2}} ] ∥x^{k}(α; δ_{k}(α)) − (x^{k} + αd̃^{k})∥^{2}.

The inequality (23) will be satisfied provided

1 − 2α^{τ−2} − 3Bηα^{\frac{τ−2}{2}} ≥ σ̄. (32)

Setting β = α^{\frac{τ−2}{2}} , we find that the quadratic 2β^{2} + (3Bη)β + (σ̄ − 1) has one positive root. Hence we can find an ᾱ_{1} > 0 such that the required inequality will be satisfied for all α ∈ (0, ᾱ_{1}].

For the second case, assume that

\frac{1}{α}∥x^{k}(α; δ_{k}(α)) − (x^{k} + αd̃^{k})∥ ≤ ∥d^{k}∥.

Then from (9),

δ_{k}(α) ≤ ηα^{τ/2}∥d^{k}∥^{2},

and so from (30),

∇f(x^{k})^{T} [x^{k} − x^{k}(α; δ_{k}(α))] (33)

≥ αd^{kT} D^{k}d^{k} + \frac{1}{α}∥x^{k}(α; δ_{k}(α)) − (x^{k} + αd̃^{k})∥^{2} − \frac{2}{α}ηC_{1}α^{τ} ∥d^{k}∥^{2} − 3Bηα^{τ/2}∥d^{k}∥^{2}.

From (7) it follows that

∥d^{k}∥^{2} ≤ \frac{1}{λ_{1}} d^{kT} D^{k}d^{k},

and so using ηC_{1} < 1, we have

∇f(x^{k})^{T} [x^{k} − x^{k}(α; δ_{k}(α))]
(34)≥ α [1 − \frac{2}{λ_{1}} α^{τ−2} − \frac{3Bη}{λ_{1}} α^{\frac{τ−2}{2}} ] d^{kT} D^{k}d^{k} + \frac{1}{α}∥x^{k}(α; δ_{k}(α)) − (x^{k} + αd̃^{k})∥^{2}.

For (23), it is sufficient that

1 − \frac{2}{λ_{1}} α^{τ−2} − \frac{3Bη}{λ_{1}} α^{\frac{τ−2}{2}} ≥ σ̄.

A similar argument to that above shows that a positive value ᾱ_{2} can be found so that this inequality is satisfied for α ∈ (0, ᾱ_{2}]. Hence the first part of the result follows by setting

ᾱ(σ̄) = min(1, ϵ_{k}/(2∥g^{k}∥), ᾱ_{1}, ᾱ_{2}).

The second part of the result (i.e., that f(x^{k}) > f(x^{k}(α, δ_{k}(α))) for sufficiently small α) is obtained by modifying the argument of Gafni and Bertsekas [5, Proposition 1 (b)]. By the mean value theorem, we can find a point ζ^{k}(α) on the line joining x^{k} to x^{k}(α, δ_{k}(α)) such that

f(x^{k}) − f(x^{k}(α; δ_{k}(α))) = ∇f(ζ^{k}(α))^{T} [x^{k} − x^{k}(α; δ_{k}(α))].

Hence from (23), for α ∈ (0, ᾱ),

\frac{1}{α}[f(x^{k}) − f(x^{k}(α, δ_{k}(α)))]
(35)  ≥ σ̄ [d^{kT} D^{k}d^{k} + \frac{1}{α^{2}} ∥x^{k}(α; δ_{k}(α)) − (x^{k} + αd̃^{k})∥^{2}]
+\frac{1}{α}[∇f(ζ^{k}(α)) − ∇f(x^{k})]^{T} [x^{k} − x^{k}(α; δ_{k}(α))].

Again, writing

x^{k} − x^{k}(α; δ_{k}(α)) = x^{k} − x^{k}(α; 0) + x^{k}(α; 0) − x^{k}(α; δ_{k}(α))

and using

\frac{∥x^{k} − x^{k}(α; 0)∥^{2}}{α^{2}} ≤ ∥g^{k}∥^{2} = ∥d̃^{k}∥^{2} + ∥d^{k+}∥^{2} ≤ (λ_{2} + 1)∥∇f(x^{k})∥^{2} ≤ (λ_{2} + 1)B^{2}, (36)

we have

\frac{1}{α}[∇f(ζ^{k}(α)) − ∇f(x^{k})]^{T} [x^{k} − x^{k}(α; δ_{k}(α))]
≥  −∥∇f(ζ^{k}(α)) − ∇f(x^{k})∥ [B\sqrt{λ_{2} + 1} + \frac{1}{α}δ_{k}(α)]
≥ −∥∇f(ζ^{k}(α)) − ∇f(x^{k})∥ [B\sqrt{λ_{2} + 1} + C_{1}α^{τ/2−1}] = O(α).

When d^{k} ≠ 0, it follows from (35) that

lim_{α→0} \frac{f(x^{k}) − f(x^{k}(α; δ_{k}(α)))}{α} ≥ σ̄d^{kT} D^{k}d^{k} > 0.

On the other hand, when d^{k} = 0,

\frac{f(x^{k}) − f(x^{k}(α; δ_{k}(α)))}{α} ≥ \frac{σ̄}{α^{2}} ∥x^{k}(α; δ_{k}(α)) − x^{k}∥^{2} + O(α)
(37)  ≥ \frac{σ̄}{α^{2}} ∥x^{k}(α; 0) − x^{k}∥^{2} − \frac{2σ̄}{α^{2}} δ_{k}(α)∥x^{k}(α; 0) − x^{k}∥ + O(α).

A straightforward application of Lemma 4.2(a) shows that

\frac{1}{α}∥x^{k}(α; 0) − x^{k}∥ ≥ ∥x^{k}(1, 0) − x^{k}∥.

Also, from Lemma 4.2(b), we have for α ∈ (0, ᾱ) that

∥x^{k}(α; 0) − x^{k}∥ ≤ α∥d^{k+}∥ ≤ αB.

Using these inequalities, together with (10), we have from (37) that

\frac{f(x^{k}) − f(x^{k}(α; δ_{k}(α)))}{α} ≥ σ̄∥x^{k}(1; 0) − x^{k}∥^{2} − \frac{2σ̄B}{α} δ_{k}(α) + O(α)
≥ σ̄∥x^{k}(1; 0) − x^{k}∥^{2} − 2σ̄BC_{1}α^{τ/2−1} + O(α).

Taking the limit, we have

lim_{α→0} \frac{f(x^{k}) − f(x^{k}(α; δ_{k}(α)))}{α} ≥ σ̄∥x^{k}(1; 0) − x^{k}∥^{2} > 0.

In either case, there is an α̂ ≤ ᾱ with the desired property.

For the main result of this section, we need to be more specific about the choice of ϵ_{k}. We now assume that

ϵ_{k} = min(ϵ, ĉ_{k}ϵ̂(x_{k})), (38)

where there is a constant B̂ such that

ĉ_{k} ∈ [1, B̂],

and

ϵ̂(x) is a continuous function of x that is zero only when x is critical.

Theorem 4.5. Suppose that ϵ_{k} satisfies condition (38), that (A) holds, and that (B) holds for x^{k}(α, 0), for all α ∈ [0, 1] and all k sufficiently large. Then every accumulation point x^{k} generated by the algorithm is critical.

Proof. The proof is quite similar to the proof of Proposition 2 of Gafni and Bertsekas [5]. Some modifications are necessary because of the inexactness in x^{k}(α) and because of the need for the quantity σ̄ in Lemma 4.4. We include most of the details here, and refer the reader to [5] for the remainder.

Suppose for contradiction that there is a noncritical point x^{∗} and a subsequence K such that lim_{k∈K} x^{k} = x^{∗}. If α_{k} denotes the steplength used in the step from x^{k} to x^{k+1}, (11) implies that

lim_{k∈K} α_{k}d^{kT} D^{k}d^{k} = 0 (39)

lim_{k∈K} \frac{1}{α_{k}} ∥x^{k}(α_{k}; δ_{k}(α_{k})) − (x^{k} + α_{k}d̃^{k})∥^{2} = 0. (40)

Taking a subsequence if necessary, assume that

lim_{k∈K} α_{k} = α^{∗}

for some α^{∗} ∈ [0, 1].

Two cases arise. First assume that α^{∗} > 0. Then from (39), d^{k} →^{k∈K} 0 and so d̃^{k} →^{k∈K} 0 and d^{k+} →^{k∈K} −∇f(x^{k}). Also from (40),

lim_{k∈K} ∥x^{k}(α_{k}; δ_{k}(α_{k})) − (x^{k} + α_{k}d̃^{k})∥ = 0, (41)

and so from (9),

lim_{k∈K} δ_{k}(α_{k}) = 0.

Using this limit together with (41), we get

x^{∗}(α^{∗}, 0) = P (x^{∗} − α^{∗}∇f(x^{∗})) = x^{∗},

which implies that x^{∗} is critical.

For the second case, take α^{∗} = 0. Then for k ∈ K sufficiently large, the test (11) will fail at least once, thus, using the notation

α_{k}^{−} = \frac{α_{k}}{β} ,

we have that

f(x^{k}) − f(x^{k}(α_{k}^{−}; δ_{k}(α_{k}^{−})) (42)

< σ \{α_{k}^{−}d^{kT} D^{k}d^{k} + _{\frac{1}{α_{k}^{−}}} ∥x^{k}(α_{k}^{−}; δ_{k}(α_{k}^{−})) − (x^{k} + α_{k}^{−}d̃^{k})∥^{2}\} .

Since, by (38), ϵ_{k} is bounded away from zero, and since it follows from (36) that ∥g^{k}∥ is bounded above, we have

lim inf_{k∈K} ϵ_{k}/∥g^{k}∥ > 0. (43)

Hence, setting σ̄ = (σ + 1)/2, Lemma 4.4 can be applied to find an ᾱ > 0 such that (23) holds for α ∈ (0, ᾱ]. Moreover, closer examination of the proof of Lemma 4.4 shows that, because of (43), the value of ᾱ can be chosen independently of x^{k}, for k sufficiently large. Now since lim_{k∈K} α_{k}^{−} = 0, we have for k sufficiently large that

∇f(x^{k})^{T} [x^{k} − x^{k}(α_{k}^{−}; δ_{k}(α_{k}^{−}))]
(44)  ≥ \frac{σ+1}{2} \{α_{k}^{−}d^{kT} D^{k}d^{k} + _{\frac{1}{α_{k}^{−}}} ∥x^{k}(α_{k}^{−}; δ_{k}(α_{k}^{−})) − (x^{k} + α_{k}^{−}d̃^{k})∥^{2}\} .

Using the mean value theorem, and combining (42) and (44), we have

\frac{1−σ}{2} \{α_{k}^{−}d^{kT} D^{k}d^{k} + _{\frac{1}{α_{k}^{−}}} ∥x^{k}(α_{k}^{−}; δ_{k}(α_{k}^{−})) − (x^{k} + α_{k}^{−}d̃^{k})∥^{2}\}
(45)  ≤ ∇f(x^{k})^{T} [x^{k} − x^{k}(α_{k}^{−}; δ_{k}(α_{k}^{−}))] − f(x^{k}) + f(x^{k}(α_{k}^{−}; δ_{k}(α_{k}^{−}))
=  [∇f(x^{k}) − ∇f(ζ^{k})]^{T} [x^{k} − x^{k}(α_{k}^{−}; δ_{k}(α_{k}^{−}))]

for some ζ^{k} on the line joining x^{k} to x^{k}(α_{k}^{−}, δ_{k}(α_{k}^{−})). Note that

\frac{1}{α_{k}^{−}} ∥x^{k} − x^{k}(α_{k}^{−}; δ_{k}(α_{k}^{−}))∥ ≤ \frac{1}{α_{k}^{−}} ∥x^{k} − x^{k}(α_{k}^{−}; 0)∥ + \frac{δ(α_{k}^{−})}{α_{k}^{−}} ≤ ∥g^{k}∥ + C_{1}(α_{k}^{−})^{τ/2−1},

which is bounded because of (36). Hence the right-hand side of (45) is o(α_{k}^{−}), and dividing both sides of (45) by α_{k}^{−} we have that

lim_{k∈K} d^{kT} D^{k}d^{k} = 0, (46)

lim_{k∈K} \frac{1}{(α_{k}^{−})^{2}} ∥x^{k}(α_{k}^{−}; δ_{k}(α_{k}^{−})) − (x^{k} + α_{k}^{−}d̃^{k})∥^{2} = 0. (47)

From (46), lim_{k∈K} d^{k} = 0 and so lim_{k∈K} d̃^{k} = 0. Since in addition d̃^{k} ∈ T^{k}, we have that x^{k} + α_{k}^{−}d̃^{k} ∈ X for k sufficiently large. Lemma 4.2(a) can be applied to show that

\frac{1}{α_{k}^{−}} ∥x^{k}(α_{k}^{−}; 0) − (x^{k} + α_{k}^{−}d̃^{k})∥
(48)  = _{\frac{1}{α_{k}^{−}}} ∥P ((x^{k} + α_{k}^{−}d̃^{k}) + α_{k}^{−}d^{k+}) − (x^{k} + α_{k}^{−}d̃^{k})∥
≥  ∥P ((x^{k} + α_{k}^{−}d̃^{k}) + d^{k+}) − (x^{k} + α_{k}^{−}d̃^{k})∥.

Meanwhile Lemma 4.2(b) implies that

\frac{1}{α_{k}^{−}} ∥x^{k}(α_{k}^{−}; 0) − (x^{k} + α_{k}^{−}d̃^{k})∥ ≤ ∥d^{k+}∥ ≤ B. (49)

Taking the sequence in (47), and using Lemma 4.1, (10), (48), and (49), we have

\frac{1}{(α_{k}^{−})^{2}} ∥x^{k}(α_{k}^{−}; δ_{k}(α_{k}^{−})) − (x^{k} + α_{k}^{−}d̃^{k})∥^{2}
≥ _{\frac{1}{(α_{k}^{−})^{2}}} ∥x^{k}(α_{k}^{−}; 0) − (x^{k} + α_{k}^{−}d̃^{k})∥^{2} − \frac{2δ_{k}(α_{k}^{−})}{α_{k}^{−}} (_{\frac{1}{α_{k}^{−}}} ) ∥x^{k}(α_{k}^{−}; 0) − (x^{k} + α_{k}^{−}d̃^{k})∥ (50)
≥  ∥P ((x^{k} + α_{k}^{−}d̃^{k}) + d^{k+}) − (x^{k} + α_{k}^{−}d̃^{k})∥^{2} − 2BC_{1}(α_{k}^{−})^{τ/2−1}.

Since the second term in this expression approaches zero, it follows from (50) that in the limit,

P (x^{∗} − ∇f(x^{∗})) = x^{∗},

and so x^{∗} is critical, again giving a contradiction.

5. Local Convergence. For the exact algorithm, the local convergence analysis is quite simple because when convergence occurs to a local minimum that satisfies the “standard” assumptions, the iterates eventually all lie on the manifold defined by the constraints which are active at the solution. This does not occur in our case, where the iterates remain in the interior of X. We thus need to ensure that the distance of the iterates to the active manifold is decreasing sufficiently quickly so as not to interfere with the (rapid) convergence in the tangent direction. Fortunately, some inherent properties of the path following projection algorithm prove to be useful here.

In this section we prove R-quadratic convergence of an algorithm in which D^{k} is a reduced Hessian. Much of the analysis is devoted to showing that steplengths of α_{k} = 1 are used for all sufficiently large k. We start by defining a scheme for choosing ϵ_{k}, then state an active set identification result. Eventual unit steplength is established in a sequence of lemmas and Theorem 5.6. We conclude with the main rate-of-convergence result in Theorem 5.7.

In addition to the assumptions made in the preceding sections, we use the following:

(D) x^{∗} is a strict local minimum that is nondegenerate, that is,

−∇f(x^{∗}) ∈ ri N(x^{∗}; X),

where ri N(x^{∗}; X) is the interior of N(x^{∗}; X) relative to the affine hull of N(x^{∗}; X).

(E) ϵ_{k} is defined as

ϵ_{k} = min(ϵ, ĉ_{k}ϵ(x_{k})),

where ϵ > 0 is a positive constant,

ϵ(x) = ∥x − P (x − ∇f(x))∥,

and ĉ_{k} ∈ [1, B̂] for some B̂ < ∞. (ĉ_{k} is a “random” quantity and need not be a function of x^{k}.)

If Assumption (B) also holds at x^{∗}, then Assumption (D) implies that there are unique scalars y_{i}^{∗} > 0 such that

− ∇f(x^{∗}) = ∑_{i∈A} y_{i}^{∗}a_{i}, (51)

where A is as defined in section 1. For later reference we introduce the notation

Ā = [a_{i}]_{i∈A},  Ā ∈ R^{n×r}, r ≤ n.

Orthonormal matrices Z ∈ R^{n×(n−r)} and Y ∈ R^{n×r} can be defined such that Z^{T} Ā = 0 and Z^{T} Y = 0.

Our relaxed definition of ϵ_{k} is motivated by the fact that calculation of x − P (x − ∇f(x)) involves a projection onto X and hence will be carried out inexactly by the algorithm of section 3. The following scheme can be used:

Algorithm to calculate ϵ_{k}:

Step 1: Given some constant Ĉ ∈ (0, 1), apply the algorithm of section 2 to find P (x^{k} − ∇f(x^{k})), terminating when the duality gap ϵ_{P}^{2} /2 satisfies the inequality

ϵ_{P} ≤ (1 − Ĉ)∥x̂^{k} − x^{k}∥,

where x̂^{k} is the latest estimate of the solution.

Step 2: Set ϵ_{k} = min(ϵ, 2∥x̂^{k} − x^{k}∥).

With the notation x̂^{k∗} = P (x^{k} − ∇f(x^{k})), Lemma 4.1 and the conditions on ϵ_{P} can be used to show that

Ĉ ≤ 1 − \frac{ϵ_{P}}{∥x̂^{k} − x^{k}∥} ≤ \frac{∥x̂^{k∗} − x^{k}∥}{∥x̂^{k} − x^{k}∥} ≤ 1 + \frac{ϵ_{P}}{∥x̂^{k} − x^{k}∥} ≤ 2 − Ĉ.

Hence

2∥x̂^{k} − x^{k}∥ = ĉ_{k}∥x̂^{k∗} − x^{k}∥,  where  ĉ_{k} ∈ [\frac{2}{2 − Ĉ} , \frac{2}{Ĉ} ] ,

and so the requirements of Assumption (D) are satisfied.

From this definition of ϵ_{k}, the following active set identification result can be proved:

Lemma 5.1. Suppose that assumptions (A), (D), and (E) hold and that (B) holds for x^{k}(α, 0), for α ∈ [0, 1] and all k sufficiently large. Assuming that x^{∗} is a limit point of the sequence {x^{k}}, we have lim_{k→∞} x^{k} = x^{∗} and I^{k} = A for all k sufficiently large.

Proof. The result follows from Lemma B.1 of Gafni and Bertsekas [5]; trivial modifications are required because of our relaxed definition of ϵ_{k}. The Assumption (B) in [5] corresponds to our Assumption (C) (see Theorem 2.8 in Burke and Moré [1]).

We next show that the steplengths do not vanish as k → ∞.

Lemma 5.2. Under the assumptions of Lemma 5.1, there is α̂ > 0 such that

α_{k} ≥ α̂

for all k sufficiently large.

Proof. From Lemma 5.1 we have that for k sufficiently large, I^{k} = A. Since d^{k} = P_{T^{k}}(−∇f(x^{k})) → 0, it follows that ∥g^{k}∥ → ∥∇f(x^{∗})∥ ≤ B. Now in Lemma 4.4, we are free to set ϵ_{k} uniformly equal to a constant ϵ̃ > 0 which is chosen so that

i ∉ I^{k} = A ⇒ a_{i}^{T} x^{k} ≤ b_{i} − 2ϵ̃∥a_{i}∥

for all sufficiently large k. Hence ϵ_{k}/∥g^{k}∥ is bounded away from zero. Now, given any σ̄ ∈ (0, 1), we can apply Lemma 4.4 to find an ᾱ(σ̄) > 0 such that for all α ∈ (0, ᾱ(σ̄)],

∇f(x^{k})^{T} [x^{k} − x^{k}(α; δ_{k}(α))] ≥ σ̄ [αd^{kT} D^{k}d^{k} + \frac{1}{α}∥x^{k}(α; δ_{k}(α)) − (x^{k} + αd̃^{k})∥^{2}] .

If we use L as an upper bound on ∇^{2}f(x) for x in some neighborhood of x^{∗}, it follows exactly as in Gafni and Bertsekas [5] that

f(x^{k}) − f(x^{k}(α; δ_{k}(α)))
≥ α(σ̄ − Lλ_{2}α)d^{kT} D^{k}d^{k} + (\frac{σ̄}{α} − L)∥x^{k}(α; δ_{k}(α)) − (x^{k} + αd̃^{k})∥^{2}.

If we choose

α̃ = sup_{σ̄∈[σ,1)} min (ᾱ(σ̄), \frac{σ̄ − σ}{Lλ_{2}} , \frac{σ̄ − σ}{L} ) ,

it follows from the line search mechanism (11) that

α_{k} ≥ α̂ =^{def} \frac{α̃}{β},

and the result follows since clearly α̂ > 0.

The next result follows easily from Lemma 5.2, Lemma B.2 of [5], and the analysis of Dunn [2]

Lemma 5.3. Under the assumptions of Lemma 5.1, we have for all k sufficiently large that

x^{k}(α; 0) = x^{∗} + Zv^{k}(α) ∈ x^{∗} + T^{k},

for all α ∈ [α_{k}, 1], where v^{k}(α) ∈ R^{n−r}. Also,

(x^{k} + α(d̃^{k} + d^{k+})) − x^{k}(α; 0) = Āȳ^{k}(α) ∈ N(x^{∗}; X), (52)

for all α ∈ [α_{k}, 1], where ȳ^{k}(α) ∈ R^{r} has ȳ_{i}^{k}(α) > C_{2}α for i = 1, ⋯, r and some constant C_{2} > 0.

Proof. We prove only the last statement concerning the lower bound on ȳ^{k}(α). Since d^{k} → 0 and d^{k+} → −∇f(x^{k}), we can combine (51) and (52) to obtain

x^{k}(α; 0) − (x^{k} + αĀy^{∗}) + Āȳ^{k}(α) → 0, (53)

where y^{∗} = {y_{i}^{∗}}_{i∈A}. Since x^{k+1} − x^{k} → 0 we have from (9) that δ_{k}(α_{k}) → 0. Hence

0 ≤ ∥x^{k}(α_{k}; 0)−x^{k}∥ ≤ ∥x^{k}(α_{k}; δ_{k}(α_{k}))−x^{k}(α_{k}; 0)∥+∥x^{k+1}−x^{k}∥ ≤ δ_{k}(α_{k})+∥x^{k+1}−x^{k}∥ → 0.

Now by Lemma 4.2, and since α ∈ [α_{k}, 1],

\frac{1}{α}∥x^{k}(α; 0) − x^{k}∥ ≤ \frac{1}{α_{k}} ∥x^{k}(α_{k}; 0) − x^{k}∥.

Since α_{k} ≥ α̂, it follows from this inequality that x^{k}(α; 0) − x^{k} → 0. Hence from (53), using the full rank of Ā, we have that

ȳ^{k}(α) → αy^{∗}.

Since y^{∗} > 0, the result follows.

Corollary 5.4. Under the assumptions of Lemma 5.1, we have for all k sufficiently large that

Z^{T} s^{k+} = 0  where  s^{k+} = x^{k}(1; 0) − (x^{k} + d̃^{k}).

In addition, a_{i}^{T} (x^{k} + s^{k+}) = b_{i} for i ∈ A.

Proof. The statement Z^{T} s^{k+} = 0 follows from the second expression in Lemma 5.3 by setting α = 1 and noting that Z^{T} d^{k+} = 0. For the second part, note that

x^{k} + s^{k+} = x^{k}(1; 0) − d̃^{k} ∈ x^{∗} + T^{k},

from the first expression in Lemma 5.3 and the fact that d̃^{k} ∈ T^{k}. Hence a_{i}^{T} (x^{k} + s^{k+}) = a_{i}^{T} x^{∗} = b_{i}, as required.

For the remainder of this section we use the following notational conventions: • γ_{k} = δ_{k}(α_{k})^{2}/2 is the final duality gap for the step from x^{k} to x^{k+1};

• The error in the approximate unit step is separated into two components:

x^{k}(1; δ_{k}(1)) − x^{k}(1; 0) = e^{k} = ẽ^{k} + e^{k+},

where ẽ^{k} = P_{T^{k}}(e^{k}) = ZZ^{T} e^{k} and e^{k+} = Y Y^{T} e^{k};

• δ_{k} denotes δ_{k}(1).

A technical result is needed before establishing eventual unit steps.

Lemma 5.5. Suppose that Assumption (C) and the assumptions of Lemma 5.1 hold and that the “special case” in the projection algorithm (i.e., x^{k} − α_{k}g^{k} ∈ X) occurs only finitely often. Then for k sufficiently large, there are positive constants C_{3}, C_{4}, and C_{5} such that

(d^{k+})^{T} s^{k+} ≥ C_{3}γ_{k−1}∥d^{k+}∥,
∥s^{k+}∥ ≤ C_{4}γ_{k−1}.

Further, if α_{k} = 1, we have that

∥e^{k+}∥ ≤ C_{5}γ_{k},
|d^{kT} ẽ^{k}| ≤ δ_{k}∥d^{k}∥.

Proof. Assume k is large enough that the “special case” never occurs after iteration k − 1. Assume further that α_{k−1} and all subsequent steplengths are bounded below by α̂, as in Lemma 5.2. Recall that x^{k} = x^{k−1}(α_{k−1}; δ_{k−1}(α_{k−1})). Let ν^{k−1} and y^{k−1} be the final values of the ν and y variables in the projection algorithm of section 3 which was used to compute x^{k}. We start by finding bounds on elements of ν^{k−1} in terms of γ_{k−1}; these are needed for the first three inequalities.

As discussed in the proof of Lemma 5.3, δ_{k}(α_{k}) → 0; that is, the projection subproblem is solved more and more accurately. Recall that the matrix equation in (20) holds at every iteration of the projection algorithm. The first part of this equation yields that

x^{k} − (x^{k−1} + α_{k−1}(d̃^{k−1} + d^{(k−1)+})) + Ay^{k−1} + qy_{m+1}^{k−1} = 0. (54)

From the second part of the equation and the choice of ϵ̃ in the proof of Lemma 5.2,

ν_{i}^{k−1} = b_{i} − a_{i}^{T} x^{k} ≥ 2ϵ̃∥a_{i}∥ > 0  for  i ∉ A.

Since γ_{k−1} = δ_{k−1}(α_{k−1})^{2}/2 = ∑_{i=1}^{m+1} ν_{i}^{k−1}y_{i}^{k−1}, we have for i ∉ A that

0 < y_{i}^{k−1} ≤ \frac{γ_{k−1}}{ν_{i}^{k−1}} ≤ \frac{γ_{k−1}}{2ϵ̃∥a_{i}∥} → 0. (55)

Since we have assumed that ν_{m+1}^{∗} is bounded away from zero for all projection subproblems,

y_{m+1}^{k−1} → 0. (56)

Now, using k − 1 instead of k in (52) and setting α = α_{k−1},

x^{k−1}(α_{k−1}; 0) − (x^{k−1} + α_{k−1}(d̃^{k−1} + d^{(k−1)+})) + Āȳ^{k−1}(α_{k−1}) = 0. (57)

Comparing (54) with (57), we have

x^{k−1}(α_{k−1}; 0) − x^{k} = Ay^{k−1} + qy_{m+1}^{k−1} − Āȳ^{k−1}(α_{k−1}).

Using (55)–(56), the full rank of Ā, and noting that ∥x^{k}−x^{k−1}(α_{k−1}; 0)∥ ≤ δ_{k−1}(α_{k−1}) → 0, we have that for some constant C_{2} > 0,

y_{i}^{k−1} → ȳ_{i}^{k−1}(α_{k−1}) ≥ C_{2}α_{k−1} for i ∈ A.

Hence for k sufficiently large, with α_{k−1} ≥ α̂, there is a constant C_{2L} > 0 such that

y_{i}^{k−1} ≥ C_{2L} for i ∈ A.

Also, by full rank of Ā and boundedness of ∇f, there is a C_{2U} > 0 such that

y_{i}^{k−1} ≤ C_{2U}.

By Assumption (C), we have for i ∈ A that

ν_{i}^{k−1} ≥ \frac{γ_{k−1}}{μy_{i}^{k−1}} ≥ \frac{γ_{k−1}}{C_{2U}μ}.

Also,

ν_{i}^{k−1} ≤ \frac{γ_{k−1}}{y_{i}^{k−1}} ≤ \frac{γ_{k−1}}{C_{2L}} .

From these last two expressions, we can define positive constants C_{7L} and C_{7U} such that

C_{7L}γ_{k−1} ≤ ν_{i}^{k−1} ≤ C_{7U}γ_{k−1}  for i ∈ A  . (58)

For the first result, note from d^{k+} → −∇f(x^{∗}) and Z^{T} d^{k+} = 0 that d^{k+} = Āt^{k}, where t^{k} → y^{∗}. Hence t^{k} > 0 for k sufficiently large. From Corollary 5.4, (20), and (58), we have for i ∈ A that

a_{i}^{T} s_{i}^{k+} = b_{i} − a_{i}^{T} x^{k} = ν_{i}^{k−1} ≥ C_{7L}γ_{k−1}.

Hence, noting that ∥t^{k}∥ ≥ ∥d^{k+}∥/∥Ā∥, we have

(d^{k+})^{T} s^{k+} = (Āt^{k})^{T} s^{k+} = t^{kT} (Ā^{T} s^{k+}) ≥ C_{7L}γ_{k−1}∥t^{k}∥ ≥ C_{3}γ_{k−1}∥d^{k+}∥,

for C_{3} = C_{7L}/∥Ā∥, giving the first result.

For the second result, we have from Lemma 5.3, Corollary 5.4, and (20) that for some u^{k} ∈ R l ^{r},

s^{k+} = Āu^{k}  and  Ā^{T} s^{k+} = [ν_{i}^{k−1}]_{i∈A}.

Hence

Ā^{T} Āu^{k} = [ν_{i}^{k−1}]_{i∈A},

which, by full rank of Ā, boundedness of ρ_{i}, and (58), gives that

∥u^{k}∥ ≤ C_{8}γ_{k−1}

for some constant C_{8} > 0. Since

∥s^{k+}∥ ≤ ∥Ā∥∥u^{k}∥,

the result follows by setting C_{4} = C_{8}∥Ā∥.

For the third inequality, we again use (20) and Lemma 5.3 to deduce that for i ∈ A,

ν_{i}^{k} = b_{i} − a_{i}^{T} x^{k}(1; δ_{k}) = a_{i}^{T} [x^{k}(1; 0) − x^{k}(1; δ_{k})] = −a_{i}^{T} e^{k+}.

Now e^{k+} = Āv^{k} for some v^{k}, so an identical argument to that of the preceding paragraph can be used to give the result.

The fourth inequality follows simply from

|d^{kT} ẽ^{k}| ≤ ∥d^{k}∥∥ẽ^{k}∥ ≤ ∥d^{k}∥∥e^{k}∥ ≤ δ_{k}∥d^{k}∥.

Theorem 5.6. Suppose that Assumptions (A), (C), (D), and (E) hold and that Assumption (B) holds in a neighborhood of x^{∗}. Suppose that Z^{T} ∇^{2}f(x^{∗})Z is positive definite and that for k sufficiently large, the tangent component of the step is given by

d̃^{k} = Z(Z^{T} ∇^{2}f(x^{k})Z)^{−1}Z^{T} d^{k}.

Suppose there is a non-negative sequence {ξ_{k}} such that lim_{k→∞} ξ_{k} = 0 and that, in addition to (9),(10), the sequence {δ_{k}} satisfies

(59)  ∥d^{k}∥δ_{k} ≤ ξ_{k}γ_{k−1}
δ_{k}^{2} ≤ ξ_{k}γ_{k−1}.

Assume that σ < 0.5 in (11). Then α_{k} = 1 for all sufficiently large k.

Proof. First, we consider the special case of x^{∗} ∈ intX, for which we have ∇f(x^{∗}) = 0. By Lemma 5.1, the two-metric gradient projection method reduces to Newton’s method when k is sufficiently large. Consequently, d^{k} = −∇f(x^{k}), and d̃^{k} = −(∇^{2}f(x^{k}))^{−1}∇f(x^{k}). By Lemma 5.3, x^{k}(1; 0) = x^{k} + d̃^{k} ∈ intX for k sufficiently large. Correspondingly, exactness of the projection yields γ_{k} = δ_{k} ≡ 0 for all such k. Now

f(x^{k}) − f(x^{k}(1; δ_{k})) = f(x^{k}) − f(x^{k} + d̃^{k})
= −∇f(x^{k})^{T} d̃^{k} − \frac{1}{2}(d̃^{k})^{T} ∇^{2}f(x^{k})d̃^{k} + o(∥d̃^{k}∥^{2})
= \frac{1}{2}d^{kT} (∇^{2}f(x^{k}))^{−1}d^{k} + o(∥d^{k}∥^{2}).

The second term on the right-hand side of (11) is zero, so for k sufficiently large, (11) is satisfied for α_{k} = 1.

In the remaining case, x^{∗} ∈ ∂X; thus, by Assumption (C), ∇f(x^{∗}) ≠ 0. Moreover, the “special case” does not occur in the projection algorithm for sufficiently large k (this follows directly from (52), which states in particular that (x^{k} − α_{k}g^{k}) − P (x^{k} − α_{k}g^{k}) ≠ 0). Since

∇f(x^{k}) = −d^{k} − d^{k+}

and

x^{k}(1; δ_{k}) = x^{k} + d̃^{k} + s^{k+} + e^{k},

we have that

f(x^{k}) − f(x^{k}(1; δ_{k})) = ∇f(x^{k})^{T} [x^{k} − x^{k}(1; δ_{k})]

−(1/2)[x^{k} − x^{k}(1; δ_{k})]^{T} ∇^{2}f(x^{k})[x^{k} − x^{k}(1; δ_{k})] + o(∥x^{k} − x^{k}(1; δ_{k})∥^{2})

= [−d^{k} − d^{k+}]^{T} [−d̃^{k} − s^{k+} − e^{k}] − (1/2)[−d̃^{k} − (s^{k+} + e^{k})]^{T} ∇^{2}f(x^{k})[−d̃^{k} − (s^{k+} + e^{k})]
+o(∥x^{k} − x^{k}(1; δ_{k})∥^{2}).

It can be easily shown that d̃^{kT} ∇^{2}f(x^{k})d̃^{k} = d^{kT} d̃^{k}, and so, after some rearrangement,

f(x^{k}) − f(x^{k}(1; δ_{k}))
= \{\frac{1}{2}d^{kT} d̃^{k} + (s^{k+} + e^{k})^{T} (s^{k+} + e^{k})\} + (d^{k+})^{T} s^{k+} − ∇f(x^{k})^{T} e^{k}

−\frac{1}{2}[s^{k+} + e^{k}]^{T} [∇^{2}f(x^{k}) + 2I][s^{k+} + e^{k}] − [s^{k+} + e^{k}]^{T} ∇^{2}f(x^{k})d̃^{k} (60)

+o(∥x^{k} − x^{k}(1; δ_{k})∥^{2}).

Now Lemma 5.5 can be used to deduce the following inequalities:

(d^{k+})^{T} s^{k+} ≥ C_{3}γ_{k−1}∥d^{k+}∥ ≥ C̄_{3}γ_{k−1},

for some C̄_{3} > 0, since ∥d^{k+}∥ → ∥∇f(x^{∗})∥ ≠ 0;

|∇f(x^{k})^{T} e^{k}| ≤ |d^{kT} ẽ^{k}| + |(d^{k+})^{T} e^{k+}| ≤ δ_{k}∥d^{k}∥ + C_{5}Bγ_{k};

\frac{1}{2} |[s^{k+} + e^{k}]^{T} [∇^{2}f(x^{k}) + 2I][s^{k+} + e^{k}]| ≤ C_{11}∥s^{k+} + e^{k}∥^{2}
≤ C_{12}γ_{k−1}^{2} + C_{13}γ_{k−1}δ_{k} + C_{14}δ_{k}^{2};

|[s^{k+} + e^{k}]^{T} ∇^{2}f(x^{k})d̃^{k}| ≤ C_{15}∥d^{k}∥(γ_{k−1} + δ_{k}),

By substituting in (60), we obtain

f(x^{k}) − f(x^{k}(1; δ_{k})) ≥ \{\frac{1}{2}d^{kT} d̃^{k} + (s^{k+} + e^{k})^{T} (s^{k+} + e^{k})\}
+C̄_{3}γ_{k−1} − δ_{k}∥d^{k}∥ − C_{5}Bγ_{k} − C_{12}γ_{k−1}^{2} − C_{13}γ_{k−1}δ_{k} − C_{14}δ_{k}^{2}
−C_{15}∥d^{k}∥γ_{k−1} − C_{15}∥d^{k}∥δ_{k} + o(∥d̃^{k} + s^{k+} + e^{k}∥^{2})
=  \{\frac{1}{2}d^{kT} d̃^{k} + (s^{k+} + e^{k})^{T} (s^{k+} + e^{k})\}
+γ_{k−1} [C̄_{3} − ξ_{k} − (C_{5}Bξ_{k}/2) − C_{12}γ_{k−1} − C_{13}δ_{k} − C_{14}ξ_{k}
−C_{15}∥d^{k}∥ − C_{15}ξ_{k}] + o(∥d̃^{k} + s^{k+} + e^{k}∥^{2}).

As k → ∞, the term in square brackets approaches C̄_{3} > 0; that is, it is positive for sufficiently large k. It is easy to see that the final o(∥d̃^{k} + s^{k+} + e^{k}∥^{2}) term is eventually dominated by the term in curly brackets. Hence, since σ ∈ (0, \frac{1}{2}), we have for k sufficiently large that

f(x^{k}) − f(x^{k}(1; δ_{k})) ≥ σ {d^{kT} d̃^{k} + ∥x^{k}(1, δ_{k}) − (x^{k} + d̃^{k})∥^{2}} ,

and so α_{k} = 1 passes the acceptance test (11) and x^{k}(1; δ_{k}) will be accepted as the new iterate.

The conditions (59) should be imposed only in the final stages of the algorithm, when there is a suspicion that the active manifold has been identified. Otherwise, it could happen that at some early iterate, x^{k} − g^{k} ∈ intX, in which case the projection is performed exactly (γ_{k} = δ_{k} = 0) and, because of (59), exact projections would be demanded at all subsequent iterations.

A similar result to Theorem 5.6 can be stated for the alternative acceptance test (12), and it can be proved in almost identical fashion.

We can now prove the final result.

Theorem 5.7. Suppose that the assumptions of Theorem 5.6 hold and that the sequence {γ_{k}} converges Q-quadratically to zero, that is, there is a constant C_{10} such that

δ_{k} ≤ C_{10}γ_{k−1}. (61)

Then the rate of local convergence of the algorithm is R-quadratic.

Proof. In the case x^{∗} ∈ intX, we actually obtain Q-quadratic convergence, since the algorithm eventually reduces to Newton’s method. We therefore focus on the case of x^{∗} ∈ ∂X.

By setting ξ_{k} = max(∥d^{k}∥, δ_{k}), it is easy to see that (61) implies (59), and so Theorem 5.6 applies. By the definition of s^{k+},

x^{k} + d̃^{k} − x^{k}(1; 0) = −s^{k+}.

Multiplying through by Z^{T}, and using the definition of d̃^{k}, we obtain

Z^{T} (x^{k} − x^{k}(1; 0)) − (Z^{T} ∇^{2}f(x^{k})Z)^{−1}Z^{T} ∇f(x^{k}) = 0. (62)

By optimality of x^{∗}, Z^{T} ∇f(x^{∗}) = 0, so by Taylor series expansion, and since ZZ^{T} + Y Y^{T} = I,

(63) Z^{T} (∇f(x) + ∇^{2}f(x)(x^{∗} − x)) = O(∥x − x^{∗}∥^{2})
⇒ Z^{T} ∇f(x^{k}) − Z^{T} ∇^{2}f(x^{k})ZZ^{T} (x^{k} − x^{∗}) = Z^{T} ∇^{2}f(x^{k})Y Y^{T} (x^{k} − x^{∗}) + O(∥x^{k} − x^{∗}∥^{2}).

Multiplying (63) by (Z^{T} ∇f(x^{k})Z)^{−1}, and adding to (62), we have

∥Z^{T} (x^{∗} − x^{k}(1; 0))∥ = O(∥Y^{T} (x^{k} − x^{∗})∥) + O(∥x^{k} − x^{∗}∥^{2}). (64)

Recall that x^{k} = x^{k−1}(1; δ_{k−1}) and that by Lemma 5.3

Y^{T} (x^{k−1}(1; 0) − x^{∗}) = 0 (65)

for all sufficiently large k. Hence, using the third inequality in Lemma 5.5, we have

∥Y^{T} (x^{k} − x^{∗})∥ = ∥Y^{T} (x^{k−1}(1; δ_{k−1}) − x^{k−1}(1, 0))∥
(66)  = ∥e^{(k−1)+}∥
≤ C_{20}γ_{k−1}.

From (64)–(66),

∥x^{∗} − x^{k+1}∥ ≤ ∥x^{∗} − x^{k}(1; 0)∥ + ∥e^{k}∥
(67)  ≤ C_{20}γ_{k−1} + δ_{k} + O(∥x^{k} − x^{∗}∥^{2}).

Now δ_{k} ≤ C_{10}γ_{k−1}, and so we can choose a constant C_{21} ≥ max(1, C_{20} + C_{10}) such that

∥x^{∗} − x^{k+1}∥ ≤ C_{21} max(γ_{k−1}, ∥x^{∗} − x^{k}∥^{2}). (68)

Given any τ < C_{21}^{−1}, we can choose an integer k̄ sufficiently large that

γ_{k−1} ≤ τ^{2}  and  ∥x^{k} − x^{∗}∥ ≤ τ for all k ≥ k̄.

An inductive argument based on (68) then shows that

∥x^{k̄+j} − x^{∗}∥ ≤ τ̄_{j},

where

τ̄_{0} = τ
τ̄_{j+1} = C_{21}τ̄_{j}^{2}.

Clearly the sequence {τ̄_{j}} is Q-quadratically convergent, so the result follows.

Results similar to Theorems 5.6 and 5.7 could be proved for other choices of d̃^{k} — for example, where d̃^{k} is a quasi-Newton or inexact Newton method step rather than the reduced Newton step. These would be of practical importance in applications in which it is difficult to compute or factor the reduced Hessian.

Finally, we note that it may be efficient to include a second “local” phase in the basic algorithm of section 2. When it appears that the active constraint set has been identified, the current iterate could be projected onto the appropriate manifold (placing it on ∂X). Standard methods for equality-constrained nonlinear programming could then be applied to identify the minimum on this manifold. However, it is likely that the basic algorithm would also be quite efficient in this situation because, as the final few iterates are close together, a good starting point for the projection would be readily available.

Acknowledgements. I thank David Stewart for helpful discussions during the preparation of this paper, and the referees and editor for their comments on the first draft.

REFERENCES

[1] J. V. Burke and J. J. Moré, On the identification of active constraints, SIAM Journal on Numerical Analysis, 25 (1988), pp. 1197–1211.

[2] J. C. Dunn, On the convergence of projected gradient processes to singular critical points, Journal of Optimization Theory and Applications, 55 (1987), pp. 203–216.

[3] , Gradient projection methods for systems optimization problems, Control and Dynamic Systems, 29 (1988), pp. 135–195.

[4] , A projected Newton method for minimization problems with nonlinear inequality constraints, Numerische Mathematik, 53 (1988), pp. 377–409.

[5] E. M. Gafni and D. P. Bertsekas, Two-metric projection methods for constrained optimization, SIAM Journal on Control and Optimization, 22 (1984), pp. 936–964.

[6] C.-G. Han, P. Pardalos, and Y. Ye, Computational aspects of an interior point algorithm for quadratic programming problems with box constraints, in Large-Scale Numerical Optimization, T. F. Coleman and Y. Li, eds., SIAM Publications, Philadelphia, PA, 1990, pp. 92–112.

[7] M. Kojima, S. Mizuno, and A. Yoshise, An o(\sqrt{n}l) iteration potential reduction algorithm for linear complementarity problems, Tech. Rep. Research Report B–217, Department of Information Sciences, Tokyo Institute of Technology, 1988.

[8] O. L. Mangasarian and T.-H. Shiau, Error bounds for monotone linear complementarity problems, Mathematical Programming, 36 (1986), pp. 81–89.

[9] S. J. Wright, Interior point methods for optimal control of discrete-time systems, Tech. Rep. MCS–P229–0491, Mathematics and Computer Science Division, Argonne National Laboratory, Argonne, IL, April 1991.

[10] Y. Ye, K. O. Kortanek, J. A. Kaliski, and S. Huang, Near-boundary behavior of primaldual potential reduction algorithms for linear programming, Tech. Rep. Working Paper Number 90–9, College of Business Administration, The University of Iowa, 1990.

[11] Y. Zhang, J. E. Dennis, and R. A. Tapia, On the superlinear and quadratic convergence of primal-dual interior point linear programming algorithms, Tech. Rep. TR90–6, Department of Mathematical Sciences, Rice University, 1990.