A Superlinear Infeasible-Interior-Point Algorithm for Monotone Complementarity Problems

Stephen Wright^{∗} and Daniel Ralph^{†}

PREPRINT MCS-P344-1292, MATHEMATICS AND COMPUTER SCIENCE DIVISION, ARGONNE NATIONAL LABORATORY

Abstract

We use the globally convergent framework proposed by Kojima, Noma, and Yoshise to construct an infeasible-interior-point algorithm for monotone nonlinear complementarity problems. Superlinear convergence is attained when the solution is nondegenerate and also when the problem is linear. Numerical experiments confirm the efficacy of the proposed approach.

1 Introduction

We consider the problem of finding a vector pair (x, y) ∈ IR^{n} × IR^{n} such that

y = f(x),  (x, y) ≥ 0,  x^{T} y = 0,  (1)

where f : IR^{n} → IR^{n} is continuously differentiable in an open set containing the nonnegative orthant of IR^{n} (denoted by IR_{+}^{n}) and monotone, that is,

(x′ − x)^{T} (f(x′) − f(x)) ≥ 0 for all x′, x ∈ IR_{+}^{n}.

Problem (1) is a monotone nonlinear complementarity problem, abbreviated as NCP. We use S to denote the solution set for (1).

Interior-point algorithms for problems of this type have been considered recently by Kojima, Noma, and Yoshise [4], Güler [3], and Potra and Ye [6]. In [4], the authors consider a broad class of infeasible-interior-point algorithms for (1) and show that, assuming continuous differentiability of f, at least one of three scenarios eventually occurs: The algorithm reaches the vicinity of a solution to (1), it reaches the vicinity of a solution to a nearby problem, or it returns an error condition that indicates that no solution of (1) exists in a certain large nonnegative neighborhood of (0, 0). The algorithm we propose in this paper falls into the class considered in [4], and so this global convergence result holds. Under additional assumptions on f and S, our algorithm exhibits superlinear convergence.

∗MCS Division, Argonne National Laboratory, 9700 South Cass Avenue, Argonne, Illinois 60439, U.S.A.

The work of this author was based on research supported by the Office of Scientific Computing, U.S. Department of Energy, under Contract W-31-109-Eng-38.

†Department of Mathematics, The University of Melbourne, Parkville, Victoria 3052, Australia. The work of this author supported by the Australian Research Council.

Our algorithm is based on the one described in [9] for linear complementarity problems (for which f has the form f(x) = Ax + q for some n × n positive semidefinite matrix A). It differs in two main respects. First, the step generated by the algorithm of [9] depends on the entire iteration history through its use of a parameter t_{k}, which counts the number of fast steps taken prior to iteration k. However, in order to fit the framework of [4], the mapping from one iterate to the next must be determined completely by the current iterate (x^{k}, y^{k}), so this explicit dependence on t_{k} must be eliminated. Second, it is no longer practical to choose the step length α to be the largest scalar in (0, 1] that satisfies certain central path/infeasibility conditions, since to do so would require frequent use of expensive root-finding techniques applied to components of f(x). Instead, we use Armijo conditions in conjunction with a backtracking line search technique, in which the initial trial step length is chosen judiciously.

We present some notation and define the algorithm in Section 2. In Section 3, we state the results from Kojima, Noma, and Yoshise [4] and show that the algorithm of Section 2 fits this framework. In Section 4, we prove the rate-of-convergence results under the assumption that the step lengths are of the same order as the current duality gap estimate x^{T} y/n. We show in Section 5 that this assumption holds when the minimizer is unique and nondegenerate. Computational experience is reported in Section 6.

Unless otherwise specified, ∥ · ∥ denotes the Euclidean norm of a vector. We frequently use (x, y) as shorthand for the vector (x^{T} , y^{T})^{T} ∈ IR^{2n}. Iteration indices (usually k) appear as superscripts on vectors and matrices and as subscripts on scalars. Subscripts are used to indicate components of vectors and matrices. The notation B((x̄, ȳ), δ) is used for the closed ball

B((x̄, ȳ), δ) = {(x, y) ∈ IR^{n} × IR^{n} | ∥(x, y) − (x̄, ȳ)∥ ≤ δ}.

The vector (1, 1, ⋯, 1) is denoted by e, while z_{+} is obtained by replacing all negative components in the vector z by zero. If (x^{∗}, y^{∗}) is a solution of (1), we can partition {1, 2, …, n} into two index sets B and N, where

x_{i}^{∗} = 0 ∀i ∈ N,  y_{i}^{∗} = 0 ∀i ∈ B.  (2)

The solution is strictly complementary if x^{∗} + y^{∗} > 0.

2 The Algorithm

Given a starting point with (x^{0}, y^{0}) > (0, 0), the algorithm generates a sequence of iterates (x^{k}, y^{k}) > (0, 0), k = 1, 2, · · ·. With each vector pair (x, y) > 0 we associate the following quantities:

μ = x^{T} y/n,  r = y − f(x),  e = (1, 1, ⋯, 1)^{T} ,
X = diag(x_{1}, x_{2}, ⋯, x_{n}), Y = diag(y_{1}, y_{2}, ⋯, y_{n}).

When (x, y) = (x^{k}, y^{k}), we sometimes attach a subscript or superscript k to the quantities μ, r, X, Y to make the dependence on (x^{k}, y^{k}) explicit.

The main computational operation at each iteration is solution of the 2n-dimensional linear system

[ \begin{matrix}Df(x) & −I \\ Y & X \end{matrix} ] [ \begin{matrix}Δx \\ Δy \end{matrix} ] = [ \begin{matrix}y − f(x) \\ −XY e + σ̃μ_{k}e \end{matrix} ] ,  (3)

where σ̃ ∈ [0, .5]. The steps generated by the algorithm have the form

x(α) = x + αΔx,  (4a)
y(α) = y + αΔy + g(α),  (4b)

where

g(α) = f(x + αΔx) − f(x) − αDf(x)Δx.  (5)

Note that

y(α) − f(x(α)) = (1 − α)(y − f(x)) = (1 − α)r.  (6)

At each iteration, the formulae (3) and (4) are used to calculate a fast step and, if it is unsuccessful, a safe step. Safe steps ensure that desirable global convergence properties hold, while fast steps ensure rapid local convergence. These two types of steps are distinguished by different choices of the centering parameter σ̃ in (3), different choices of the initial trial step size for the Armijo line search, and slightly different acceptance criteria for the step length. While the formal treatment and theoretical utility of fast and safe steps are quite different, the distinction between them need not be so wide in practice. The wide latitude allowed to the user in the choice of σ̃ and initial trial step size for safe steps means that safe steps can be made to perform like fast steps during the later stages of the algorithm.

The overall algorithm is parametrized by a variety of positive scalar constants, which we specify now and explain later, as they arise in subsequent discussions:

χ ∈ (0, 1),  σ̄ ∈ (0, \frac{1}{2}),  ᾱ ∈ (0, 1],
κ ∈ (0, 1),  τ̂ ∈ (0, 1),  β_{0} > 0,
0 < γ_{min} < γ_{max} ≤ \frac{1}{2},  γ̄ ∈ (0, \frac{1}{2}),  ρ ∈ (0, min(γ̄^{1/τ̂}, 1 − κ)).

The starting point (x^{0}, y^{0}) is assumed to satisfy

∥r^{0}∥ ≤ β_{0}μ_{0},  x_{i}^{0}y_{i}^{0} ≥ γ_{max}μ_{0}.  (7)

The main algorithm can now be specified.

for k = 0, 1, 2, ⋯,

if μ_{k} = 0,

then terminate with solution (x^{k}, y^{k}); (x^{k+1}, y^{k+1}) ← fast(x^{k}, y^{k}); if μ_{k+1} > ρμ_{k} then (x^{k+1}, y^{k+1}) ← safe(x^{k}, y^{k}); end if

end for.

Note that the fast step is taken if it produces at least a factor of ρ decrease in the complementarity gap μ. Otherwise, the algorithm reverts to the safe step. The coefficient matrix in (3) is the same for both fast and safe steps, so only one matrix factorization is required per iteration.

The safe step procedure is defined as follows.

safe(x, y):

choose σ̃ ∈ [σ̄, \frac{1}{2}], α^{0} ∈ [ᾱ, 1]; solve (3) to find (Δx, Δy); calculate

γ̃ = min( min_{i=1,···,n} x_{i}y_{i}/μ, γ_{max});  (8)

choose α to be the first element in the sequence α^{0}, χα^{0}, χ^{2}α^{0}, ⋯, such that the following conditions are satisfied:

x_{i}(α)y_{i}(α) ≥ γ̃ x(α)^{T} y(α)/n,  (9a)
κα(1 − σ̃)μ ≤ μ − x(α)^{T} y(α)/n ≤ αμ;  (9b)

return (x(α), y(α)).

A nonzero centering term is used, allowing us to move a nontrivial distance along the search direction while staying in the set defined by

{(x, y) | x_{i}y_{i} ≥ γμ}  (10)

where γ = γ̃ (see formula (9a)). In the second acceptance condition (9b), the left inequality ensures a “sufficient decrease” in the objective function. A condition of this kind is present in most optimization algorithms based on line searchs. The purposes of the right inequality in (9b) is to prevent improvement in the complementarity gap μ from outpacing improvement in the infeasibility, measured by r. The relevant result is proved in Lemma (2.2).

Fast-step calculations are a little more complicated. Since it is not permissible in the framework of [4] to maintain a counter t_{k} of the number of fast steps taken prior to iteration k, we form an estimate t̃ of its value by examining the properties of the current iterate. The integer t̃ is in turn used to form γ̂ and β̂, which are used in the acceptance criteria for the step length α. We show in later analysis that γ̂ and β̂ have essentially the same properties as the quantities γ̃ and β̃ in the algorithm of [10]. fast(x, y):

solve (3) with σ̃ = 0 to find (Δx, Δy); calculate γ̃ as in (8); define γ̂ = γ_{min} + γ̄(γ̃ − γ_{min}); calculate

β̃ = \frac{β_{0}μ}{∥r∥},  (11)

and define t̃ to be the smallest positive integer such that

∏_{j=1}^{t̃}(1 − γ̄^{j}) ≤ β̃,

(with t̃= 0 if β̃_{k} ≥ 1);

if t̃= ∞ then return(x, y); define β̂ = γ̄^{t̃+1} (with β̂ = 1 if r = 0); define

α^{0} = 1 − \frac{μτ̂}{min(γ̃ − γ̂, β̂)};  (12)
if  α^{0} ≤ 0  then return(x, y);

choose α to be the first element in the sequence α^{0}, χα^{0}, χ^{2}α^{0}, ⋯, such that the following conditions are satisfied:

x_{i}(α)y_{i}(α) ≥ γ̂ x(α)^{T} y(α)/n,  (13a)
x(α)^{T} y(α) ≥ (1 − α)(1 − β̂)nμ;  (13b)

return(x(α), y(α)).

There is no centering component in this step, since σ̃ = 0. It is therefore necessary to relax the value of γ in (10) from γ̃ to γ̂ to ensure that we can move a nontrivial distance along this direction while staying in this set. The second acceptance criterion (13b) is again motivated by our wish to not allow improvement in μ to outpace improvement in r. There is no Armijo condition for the fast step. Instead, a “sufficient decrease” condition is enforced in the main program, since the fast step is accepted only if μ_{k+1} ≤ ρμ_{k}.

We stress again that we use the subscripted notation γ̃_{k}, γ̂_{k}, β̂_{k}, β̃_{k}, and t̃_{k} to denote the values of γ̃, γ̂, β̂, β̃, and t̃ associated with the k-th iterate (x^{k}, y^{k}). We use α_{k} to denote the value of α used by the step that is actually accepted, whether it comes from fast or safe. For the purposes of subsequent analysis, we also define

ν_{k} = ∏_{j=0}^{k−1}(1 − α_{j}).

From (6), we have that

r^{k} = ν_{k}r^{0}.  (14)

Suppose that successful fast steps are taken at iterations k_{1}, k_{2}, ⋯, with

0 ≤ k_{1} < k_{2} < ⋯,

and that safe steps are taken at all other iterations. If we define t_{k} to be the total number of fast steps accepted prior to iteration k, then clearly

t_{k_{l}} = l − 1,  l = 1, 2, ⋯,  (15)

and

t_{k} = l,  k = k_{l} + 1, ⋯, k_{l+1},  l = 1, 2, · · · .  (16)

The following result relates the settings of γ̂_{k} and β̂_{k} on each fast step to t_{k}. Its proof is rather technical and is relegated to the appendix.

Lemma 2.1 If r^{0} ≠ 0, then for each iteration k_{l}, l = 1, 2, ⋯, at which a fast step is taken, we have for l = 1, 2, ⋯, that the following inequalities hold.

t̃_{k_{l}} ≤ t_{k_{l}} = l − 1,  (17)
γ̃_{k_{l}} − γ_{min} ≥ γ̄^{l−1}(γ_{max} − γ_{min}),  (18)
γ̃_{k_{l}} − γ̂_{k_{l}} ≥ γ̄^{l−1}(1 − γ̄)(γ_{max} − γ_{min}),  (19)
β̂_{k_{l}} ≥ γ̄^{l},  (20)
min(γ̃_{k_{l}} − γ̂_{k_{l}}, β̂_{k_{l}}) ≥ γ̄^{l}(γ_{max} − γ_{min}).  (21)

If r^{0} = 0, the inequalities (18), (19), and (21) are satisfied.

Finally, we show that the improvement in complementarity μ_{k}/μ_{0} cannot exceed the improvement in feasibility ∥r^{k}∥/∥r^{0}∥, modulo a constant factor β_{L} > 0.

Lemma 2.2 If r^{0} ≠ 0, all iterates (x^{k}, y^{k}) satisfy

\frac{μ_{k}}{∥r^{k}∥} ≥ \frac{β_{L}}{β_{0}} ,  (22)

where

1 > β_{L} = ∏_{j=1}^{∞}(1 − γ̄^{j}) > e^{−3/2}.

Proof. The first inequality follows from (9b), (13b), and (6), with similar arguments to those in the proof of Lemma 2.1, which we do not repeat here. The inequality β_{L} > e^{−3/2} is proved in [10, Lemma 3.1].

3 Global Convergence

Kojima, Noma, and Yoshise [4, Section 4.1] analyze an algorithm that would be equivalent to our algorithm if we allowed only safe steps to be taken. They show that it fits into the framework that allows their global convergence result to hold. We show in this section that the use of fast steps does not disqualify our algorithm from the framework of [4], and so the main global convergence result of that paper holds.

Throughout the section we assume only that f is continuously differentiable and monotone.

The model algorithm of Kojima, Noma, and Yoshise consists of three fundamental components:

• An admissible set Ω ⊂ (IR_{++}^{n} × IR_{++}^{n}) ∪ S (where IR_{++}^{n} is the strictly positive orthant in IR^{n}) to which all iterates are confined. In our case,

Ω = {(x, y) ≥ 0 | ∥y − f(x)∥ ≤ (β_{0}/β_{L})(x^{T} y/n), x_{i}y_{i} ≥ γ̄(x^{T} y/n), i = 1, ⋯, n} .

• A merit function ψ(x, y), which in our case is simply ψ(x, y) = x^{T} y = nμ.

• An algorithmic mapping A that produces a new point (x^{+}, y^{+}) ∈ Ω from a given (x, y) ∈ Ω. In our case, A is the fast step calculation whenever it produces a decrease factor of at least ρ in μ; otherwise, A is the safe-step calculation.

The main result of [4] is obtained if the following condition is satisfied. (Note that the items in this condition are obtained by combining Conditions 2.1 and 2.3 of [4].)

Condition 1 The admissible set Ω, merit function ψ, and algorithmic mapping A satisfy the following conditions:

(i) (x^{0}, y^{0}) ∈ Ω;

(ii) Ω = Ω_{++} ∪ Ω_{S} is a subset of (IR_{++}^{n} × IR_{++}^{n}) ∪ S and is closed in (IR^{n} × IR^{n}), where

Ω_{++} = Ω ∩ (IR_{++}^{n} × IR_{++}^{n}), Ω_{S} = Ω ∩ S;

(iii) ψ is a real-valued and continuous function on Ω_{++};

(iv) A is a point-to-set mapping from Ω into the collection of nonempty subsets of Ω. For every (x̄, ȳ) ∈ Ω_{++}, there exist positive numbers δ and ϵ such that if (x, y) ∈ B((x̄, ȳ), δ) ∩ Ω and (x^{+}, y^{+}) ∈ A(x, y), then either

(x^{+}, y^{+}) ∈ Ω_{S}

or

(x^{+}, y^{+}) ∈ Ω_{++} and ψ(x^{+}, y^{+}) ≤ ψ(x, y) − ϵ;

(v) sup{x^{T} y | (x, y) ∈ Ω_{++}, ψ(x, y) ≤ T } < ∞ for every sufficiently large T ∈ IR;

(vi) if (x^{+}, y^{+}) ∈ A(x, y), then

y^{+} − f(x^{+}) = (1 − θ)(y − f(x)),

for some θ ∈ (0, 1].

It is immediately clear from our definition of Ω, ψ, and A that our algorithm satisfies Condition 1(i), (ii), (iii), (v), and (vi). The following two lemmas show that the remaining condition, (iv), is also satisfied. The first of these lemmas essentially shows that the requirements (9) on the step length α_{k} in a safe-step calculation are satisfied for all α_{k} sufficiently small. In keeping with the definition of the generic algorithm of [4], in which the mapping A does not depend explicitly on the iteration history, we state the result without reference to the iteration counter k. The proof of the first result can be found in the appendix.

Lemma 3.1 Suppose that (x̄, ȳ) ∈ Ω_{++}. Then there exist δ̂ > 0 and α̂ ∈ (0, 1] such that for any (x, y) ∈ B((x̄, ȳ), δ̂) ∩ Ω, the calculations (3), (4), and (5) applied to the point (x, y) with σ̃ ∈ [σ̄, 1/2] will yield (x(α), y(α)) satisfying the conditions (9) for all α ∈ [0, α̂].

Lemma 3.2 Given any (x̄, ȳ) ∈ Ω_{++}, there are constants δ > 0 and ϵ > 0 such that if one step of the algorithm is applied to any point (x, y) ∈ B((x̄, ȳ), δ) ∪ Ω, the new point (x^{+}, y^{+}) generated by this process has (x^{+})^{T} y^{+} ≤ x^{T} y − ϵ. Hence Condition 1(iv) is satisfied.

Proof. It is easy to check that the result (x^{+}, y^{+}) of any fast or safe step satisfies the conditions for membership in Ω, so we need only find an ϵ > 0 that satisfies the decrease condition.

Let α̂ and δ̂ be as defined in Lemma 3.1. We can clearly choose a δ ≤ δ̂ and C_{2} > 0 such that

x^{T} y ≥ C_{2} ∀(x, y) ∈ B((x̄, ȳ), δ̂).

If a safe step is calculated from the point (x, y), then from Lemma 3.1 we have that (9) are satisfied for all α ∈ [0, α̂]. Because of the backtracking nature of the step length procedure, it is easy to see that the step length parameter α actually generated by the safe-step procedure satisfies

α ≥ min(χα̂, ᾱ),

where ᾱ is the lower bound on the initial trial step length. Setting ϵ^{(1)} = (κ/2)C_{2} min(χα̂, ᾱ), we have from (9b) that

(x^{+})^{T} y^{+} = x(α)^{T} y(α) ≤ x^{T} y − κ(1 − σ̃)αx^{T} y ≤ x^{T} y − κ(1/2)C_{2} min(χα̂, ᾱ) = x^{T} y − ϵ^{(1)}.

The other algorithmic possibility is that a successful fast step is taken from (x, y). Because of the acceptance criterion associated with such a step, we have that

(x^{+})^{T} y^{+} ≤ ρx^{T} y = x^{T} y − (1 − ρ)x^{T} y ≤ x^{T} y − ϵ^{(2)},

where we have defined ϵ^{(2)} = (1 − ρ)C_{2}.

The result of the lemma follows by taking ϵ = min(ϵ^{(1)}, ϵ^{(2)}).

Having shown that Condition 1 holds in our case, we can state the main results, both of which follow immediately from results in Kojima, Noma, and Yoshise [4]. The first result is essentially [4, Lemma 2.2].

Theorem 3.3 The sequence (x^{k}, y^{k}) behaves in one of the following three ways.

(A) (x^{k}, y^{k}) ∈ S for some k < ∞.

(B) The sequence {(x^{k}, y^{k})} is bounded, and every limiting point of the sequence belongs to S ∩ Ω.

(C) The sequence {(x^{k}, y^{k})} is unbounded.

The second result follows immediately from [4, Theorem 2.5].

Theorem 3.4 Let ϵ be any small positive number and M be any positive number. Then there exists a finite integer p such that at least one of the following three alternatives occurs.

(A′) (x^{p})^{T} y^{p} < ϵ and ν_{p}∥r^{0}∥ < ϵ.

(B′) r^{0} ≠ 0, and ν_{p}∥r^{0}∥ < ϵ.

(C′) r^{0} ≠ 0, ν_{p}∥r^{0}∥ ≥ ϵ, and ν_{p}(r^{0})^{T} x^{p} − (x^{p})^{T} y^{p} ≥ ν_{p}M.

In the case (C′), the set

T (M) = {(x, y) ∈ (IR_{+}^{n} × IR_{+}^{n}) | (r^{0})^{T} x ≤ M}

contains no solutions of (1).

4 Superlinear Local Convergence

We show here that if the algorithm does not terminate finitely or diverge (that is, if alternative (B) of Theorem 3.3 occurs), then the complementarity gap converges superlinearly to zero, under certain assumptions. The Q-order of convergence is at least 1 + τ̂, where τ̂ ∈ (0, 1) is the user-defined parameter in the main algorithm. The analysis in this section is similar to that of Wright [11], which deals with the case of linear f, but there are a number of complications because of the nonlinearity of f and the simpler line search procedure.

We define a near-solution subset of Ω_{++} by

Ω_{++}(δ̃) = {(x, y) ∈ Ω_{++} | dist((x, y), S) ≤ δ̃},

where

dist((x, y), S) =^{△} min{∥(x^{∗}, y^{∗}) − (x, y)∥ | (x^{∗}, y^{∗}) ∈ S}.

Assumption 1 The sequence generated by the algorithm is bounded and does not terminate finitely. Moreover, there is δ̃ > 0 such that

(i) for all (x, y) ∈ Ω_{++}(δ̃) there is a constant C_{3} > 0 such that a fast step (Δx, Δy) calculated via (3) with σ̃ = 0 satisfies ∥(Δx, Δy)∥ ≤ C_{3}μ;

(ii) Df(·) is Lipschitz continuous on the set {x | (x, y) ∈ Ω_{++}(δ̃) for some y}.

We show in Wright [11] that this assumption and therefore the conclusions of this section hold when f is linear and S contains a strictly complementary solution (not necessarily unique).

In our first result, we show that the initial step length in the fast steps (12) approaches 1, giving the possibility of substantial progress. The result is a consequence of Assumption 1 and Theorem 3.3.

Lemma 4.1 If Assumption 1 holds, then

(i) μ_{k} ↓ 0;

(ii) μτ̂ ; k/γ̄^{t_{k}} ↓ 0;

(iii) μτ̂ ; k/ min(γ̃_{k} − γ̂_{k}, β̂_{k}) → 0.

Proof. Assumption 1 implies that alternative (B) of Theorem 3.3 holds, and so all limit points of {(x^{k}, y^{k})} are in S ∩ Ω. Since μ = 0 for all (x, y) ∈ S ∩ Ω and since {μ_{k}} is a decreasing sequence, (i) follows.

For (ii), we consider the effect of iteration k − 1 on the value of t_{k}. If a safe step was taken at iteration k − 1, then we have t_{k} = t_{k−1} and μ_{k} < μ_{k−1}, so certainly

\frac{μτ̂ ; k}{γ̄^{t_{k}}} < \frac{μτ̂ ; k−1}{γ̄^{t_{k−1}}} .

If a fast step is taken, then t_{k} = t_{k−1} + 1 and μ_{k} ≤ ρμ_{k−1}. Hence, by the definition of ρ we have

\frac{μτ̂ ; k}{γ̄^{t_{k}}} ≤ \frac{ρτ̂μτ̂ ; k−1}{γ̄γ̄^{t_{k−1}}} < \frac{γ̄μτ̂ ; k−1}{γ̄γ̄^{t_{k−1}}} = \frac{μτ̂ ; k−1}{γ̄^{t_{k−1}}} .  (23)

Hence, the sequence μτ̂ ; k/γ̄^{t_{k}} is decreasing. To see that it decreases to zero, consider two cases. First, if only a finite number of fast steps are taken, we have that t_{k} is constant for all k sufficiently large, while μ_{k} ↓ 0, so we obtain the result in this case. Second, if there are an infinite number of fast steps, we have from (23) that μτ̂ ; k/γ̄^{t_{k}} decreases by a factor of at least ρτ̂/γ̄ < 1 on each such step. Hence, the subsequence corresponding to the fast steps decreases geometrically to zero, so by monotonicity the whole sequence converges to zero, and (ii) holds in this case also.

For (iii), we note from (17) and (21) that

\frac{μτ̂ ; k}{min(γ̃_{k} − γ̂_{k}, β̂_{k})} ≤ \frac{μτ̂ ; k}{γ̄^{t_{k}+2}} ,

and hence the sequence in question is majorized by a sequence that is decreasing monotonically to zero, giving the desired result.

If Assumption 1 holds, then we can define positive constants C_{4} and C_{5} such that for the fast steps (Δx, Δy) we have

|Δx_{i}Δy_{i}| ≤ C_{4}μ^{2},  i = 1, ⋯, n,  (24)

and

∥Δx∥ ≤ C_{5}μ  (25)

for all k sufficiently large. Defining ξ̄ = ξ + C_{5}μ_{0}, where μ_{0} = (x^{0})^{T} y^{0}/n, we have under the prevailing assumption that

∥x + αΔx∥ ≤ ξ̄.  (26)

The next result show that the initial step length (12) eventually satisfies the acceptance criteria (13a) and (13b), making reduction of the step length unnecessary.

Lemma 4.2 Suppose that Assumption 1 holds. Define the constant C_{6} by

C_{6} = 2(C_{4} + ξ̄LC_{5}^{2}).  (27)

Then provided that (x, y) ∈ Ω_{++}(δ̃) and

\frac{C_{6}μ}{min(β̂, γ̃ − γ̂)} ≤ 1,

where β̂, γ̃, and γ̂ are defined in fast, the conditions (13a) and (13b) are satisfied for all

α ∈ [0, 1 − \frac{C_{6}μ}{min(β̂, γ̃ − γ̂)}] .  (28)

Proof. We start with the condition (13a). Using (3), (4), (5), (8), (24), (25), (26), and Assumption 1, we have for α ∈ [0, 1] that

x_{i}(α)y_{i}(α)

= (x_{i} + αΔx_{i})(y_{i} + αΔy_{i} + f_{i}(x + αΔx) − f_{i}(x) − αDf_{i}(x)Δx)
≥ x_{i}y_{i} + α(−x_{i}y_{i}) + α^{2}Δx_{i}Δy_{i} + α(x_{i} + αΔx_{i}) ∫_{0}^{1} [Df_{i}(x + θαΔx) − Df_{i}(x)] Δx dθ
≥ γ̃μ(1 − α) − α^{2}C_{4}μ^{2} − α∥x + αΔx∥ sup_{θ∈(0,1)} ∥Df(x + θαΔx) − Df(x)∥∥Δx∥
≥ γ̃μ(1 − α) − α^{2}C_{4}μ^{2} − α^{2}ξ̄L∥Δx∥^{2}
≥ γ̃μ(1 − α) − C_{4}μ^{2} − ξ̄LC_{5}^{2}μ^{2}.  (29)

Meanwhile, from the right-hand side of (13a), we have

x(α)^{T} y(α)
= (x + αΔx)^{T} (y + αΔy + f(x + αΔx) − f(x) − αDf(x)Δx)
= x^{T} y(1 − α) + α^{2}Δx^{T} Δy + α(x + αΔx)^{T} ∫_{0}^{1} [Df(x + θαΔx) − Df(x)] Δx dθ
≤ x^{T} y(1 − α) + α^{2}C_{4}nμ^{2} + α^{2}ξ̄L∥Δx∥^{2}
≤ x^{T} y(1 − α) + C_{4}nμ^{2} + ξ̄LC_{5}^{2}μ^{2}.  (30)

Combining (29) and (30), we see that (13a) will hold provided that

γ̃μ(1 − α) − C_{4}μ^{2} − ξ̄LC_{5}^{2}μ^{2} ≥ γ̂μ(1 − α) + γ̂C_{4}μ^{2} + (γ̂/n)ξ̄LC_{5}^{2}μ^{2}.  (31)

By replacing γ̂ by 1 in the last term and rearranging, we see that (31) holds if

(γ̃ − γ̂)μ(1 − α) − 2C_{4}μ^{2} − 2ξ̄LC_{5}^{2}μ^{2} ≥ 0 ⟺ 1 − α ≥ \frac{2(C_{4}μ + ξ̄LC_{5}^{2}μ)}{(γ̃ − γ̂)}.

From (27) we see that this inequality is equivalent to

1 − α ≥ \frac{C_{6}μ}{(γ̃ − γ̂)},  (32)

and, since (32) holds for all α in the range of (28), the condition (13a) also holds for all α in this range.

Turning now to (13b), we can perform a similar calculation to (30) (but this time seeking a lower bound) to obtain

x(α)^{T} y(α) ≥ x^{T} y(1 − α) − C_{4}nμ^{2} − ξ̄LC_{5}^{2}μ^{2}.  (33)

From (27) and (33) we deduce that

x(α)^{T} y(α) ≥ x^{T} y (1 − α − C_{4}μ − ξ̄LC_{5}^{2}μ/n) ≥ x^{T} y(1 − α − C_{6}μ).

Therefore a sufficient condition for (13b) is that

x^{T} y(1 − α − C_{6}μ) ≥ (1 − α)(1 − β̂)x^{T} y ⟺ 1 − α ≥ \frac{C_{6}μ}{β̂}.  (34)

Condition (34) is certainly satisfied for all α in the range (28), so we obtain the desired result.

We now give some threshold conditions on μ that ensure that the initial trial step length α^{0} is accepted not only by the conditions (13a), (13b), but also by the outer loop, that is, it yields a reduction factor of at least ρ in the complementarity μ.

Lemma 4.3 Suppose that Assumption 1 holds and that (x, y) ∈ Ω_{++}(δ̃) with

μ ≤ (\frac{1}{C_{6}} )^{1/(1−τ̂)}  (35)

and

\frac{μτ̂}{min(γ̃ − γ̂, β̂)} ≤ \frac{ρ}{2}.  (36)

Then the initial step length α^{0} for the fast step defined by (12) satisfies the acceptance criteria (13), and the fast step is accepted by the main algorithm.

Proof. Note from (35) that

μ^{1−τ̂} ≤ \frac{1}{C_{6}} ⇒ μτ̂ ≥ C_{6}μ.

Using (21) and (36), we have

\frac{C_{6}μ}{min(β̂, γ̃ − γ̂)} ≤ \frac{μτ̂}{min(β̂, γ̃ − γ̂)} ≤ \frac{ρ}{2} ≤ 1.  (37)

We can use Lemma 4.2 and the first inequality in (37) to deduce that α^{0} lies in the range (28) and hence satisfies (13).

To demonstrate acceptance of the fast step by the main algorithm, we need to show that

x(α^{0})^{T} y(α^{0}) ≤ ρx^{T} y.

Using (30) again, we have

x(α)^{T} y(α) ≤ x^{T} y(1 − α) + C_{4}(x^{T} y)μ + ξ̄LC_{5}^{2}(x^{T} y/n)μ ≤ x^{T} y(1 − α + C_{6}μ).  (38)

From this expression, we deduce that the fast step (with α = α^{0}) is accepted by the main algorithm if

1 − α^{0} + C_{6}μ ≤ ρ.  (39)

Note from (36) that

1 − α^{0} = \frac{μτ̂}{min(γ̃ − γ̂, β̂)} ≤ \frac{ρ}{2}.

Using both (35) and (36), we have

C_{6}μ ≤ C_{6}μ^{1−τ̂}μτ̂ ≤ μτ̂ ≤ \frac{μτ̂}{min(γ̃ − γ̂, β̂)} ≤ \frac{ρ}{2}.

The bound (39) clearly follows from these expressions, so the fast step is accepted. We are ready for the main superlinear convergence result.

Theorem 4.4 Suppose that Assumption 1 holds. Then the algorithm eventually always takes fast steps, and

(i) the sequence {μ_{k}} converges superlinearly to zero with Q-order at least 1 + τ̂, and

(ii) the sequence {ν_{k}} converges superlinearly to zero with R-order at least 1 + τ̂.

Proof. Because of Lemma 4.1, the threshold conditions (35), (36) will be satisfied for all sufficiently large k, so fast steps will eventually always be taken.

For the rate-of-convergence result, note from (12), (30), and τ̂ ∈ (0, 1) that there is a constant C_{10} such that

μ_{k+1} ≤ C_{10} \frac{μ_{k}^{1+τ̂}}{min(γ̃ − γ̂, β̂)} ≤ C_{10} \frac{μ_{k}^{1+τ̂}}{γ̄^{t_{k}+2}} ,  (40)

for all fast iterates k for k sufficiently large. Superlinear convergence of μ_{k} with Q-order at least 1 + τ̂ follows by standard arguments; see Wright [9, Theorem 6.3] and Wright and Zhang [12, Theorem 5.2]. For (ii), we have from (22) that

ν_{k} = \frac{∥r^{k}∥}{∥r^{0}∥} ≤ \frac{β_{0}}{β_{L}∥r^{0}∥}μ_{k},

so {ν_{k}} is majorized by a sequence that converges with Q-order at least 1 + τ̂, giving the result.

5 Nondegenerate Problems

We have already noted that Assumption 1 holds if f is linear and S contains a strictly complementary solution. In this section, we consider nonlinear f. We show that if (1) has a unique, nondegenerate solution, and if there is a vector x̄ > 0 for which f(x̄) > 0, then Assumption 1 is satisfied.

Assumption 2 (i) The solution set S contains the single vector pair (x^{∗}, y^{∗}), where x^{∗} + y^{∗} > 0;

(ii) There is a point x̄ ∈ IR^{n} such that (x̄, f(x̄)) > 0;

(iii) The submatrix

[Df(x^{∗})_{ij}]_{i∈B, j∈B}

(with B defined in (2)) is nonsingular, and Df(·) is Lipschitz continuous in a neighborhood of x^{∗}.

We start by showing boundedness of the iteration sequence.

Lemma 5.1 Suppose that Assumption 2 holds. Then the sequence {(x^{k}, y^{k})} is bounded.

Proof. Define ϵ^{∗} > 0 by

ϵ^{∗} = \frac{1}{2} min_{i=1,···,n} f_{i}(x̄) > 0,

and choose M > 0 large enough that (r^{0})^{T} x^{∗} < M. Then case (C′) of Theorem 3.4 cannot occur, so there is a finite integer K such that we have that ν_{K}∥r^{0}∥ < ϵ^{∗}. Since the sequence {ν_{k}} is decreasing, we have ν_{k}∥r^{0}∥ < ϵ^{∗} for all k ≥ K.

By monotonicity of f, we have

(x^{k} − x̄)^{T} (f(x^{k}) − f(x̄)) ≥ 0.  (41)

Hence, using (14) and the inequalities x̄ > 0, y^{k} > 0, and (x^{k})^{T} y^{k} < (x^{0})^{T} y^{0}, we have

(x^{k})^{T} [f(x̄) − f(x^{k})] ≤ x̄^{T} f(x̄) − x̄^{T} f(x^{k})
⇒  (x^{k})^{T} [f(x̄) + r^{k}] ≤ x̄^{T} f(x̄) + x̄^{T} r^{k} − x̄^{T} y^{k} + (x^{k})^{T} y^{k}
⇒  (x^{k})^{T} [f(x̄) + ν_{k}r^{0}] ≤ x̄^{T} f(x̄) + ν_{k}x̄^{T} r^{0} + (x^{0})^{T} y^{0}.  (42)

For k ≥ K we have from (42) that

ϵ^{∗}e^{T} x^{k} ≤ (x^{k})^{T} [f(x̄) − ϵ^{∗}e] ≤ (x^{k})^{T} [f(x̄) + ν_{k}r^{0}] ≤ x̄^{T} f(x̄) + ϵ^{∗}e^{T} x̄ + (x^{0})^{T} y^{0}.

Hence, since ϵ^{∗} > 0 and x^{k} > 0, it follows that {x^{k}} is bounded.

Boundedness of {y^{k}} follows from

∥y^{k}∥ = ∥f(x^{k}) + r^{k}∥ ≤ ∥f(x^{k})∥ + ∥r^{0}∥.

Because of boundedness, case (C) of Theorem 3.3 does not occur, so the sequence either terminates finitely or converges to (x^{∗}, y^{∗}). An estimate of the size of the step from (3) can also be easily obtained.

Lemma 5.2 Suppose that Assumption 2 holds. Then there is δ̃ > 0 such that for all (x, y) ∈ Ω_{++}(δ̃), the step calculated from (3) satisfies

∥(Δx, Δy)∥ = O(μ).  (43)

Proof. Consider the nonlinear system of equations F (x, y) = (0, 0) defined by

F (x, y) = [ \begin{matrix}f(x) − y \\ XY e \end{matrix} ] ,

and note that F (x^{∗}, y^{∗}) = (0, 0) with nonsingular Jacobian

DF (x^{∗}, y^{∗}) = [ \begin{matrix}Df(x^{∗}) & −I \\ Y^{∗} & X^{∗} \end{matrix} ] .

We can choose constants M > 0 and δ̃ > 0 such that

(x, y) ∈ B((x^{∗}, y^{∗}), δ̃) ⟹ ∥(DF (x^{∗}, y^{∗}))^{−1}∥ ≤ M.

Also, we know from the definition of Ω that

(x, y) ∈ Ω_{++} ⟹ ∥(y − f(x), XY e)∥ = O(μ).

Combining these estimates with (3), we have for (x, y) ∈ Ω_{++}(δ̃) that

∥(Δx, Δy)∥ ≤ M∥(y − f(x), XY e − σ̃e)∥ = O(μ),

giving the result.

We can now state the main result of this section.

Theorem 5.3 Suppose that Assumption 2 holds and that finite termination of the algorithm at the solution point (x^{∗}, y^{∗}) does not occur. Then Assumption 1 is satisfied.

6 Computational Results

The method described here has been implemented and tested on some small NCPs from the collection of Dirkse [2] and some larger extended linear-quadratic programming problems from Zhu and Rockafellar [13]. The following parameter settings are used in our code:

χ = .9,  σ̄ = .01,  ᾱ = .95,  κ = .1,  γ̄ = .5,
γ_{min} = 10^{−4},  γ_{max} = 10^{−2},  τ̂ = .9,  ρ = 0.2.

We modify the algorithm slightly to use a different value for the reduction factor χ for the fast steps. This value, χ_{fast}, is set to 0.98. The value of sigma is chosen at each safe iteration according to the formula

σ = max(σ̄, min(μ, σ_{max})),

where we used σ_{max} = .25. We also avoid calculation of fast steps when they are not likely to succeed. In our implementation, the fast step is not calculated as long as μ_{k} > 0.1; the safe step is always taken when this condition holds.

Successful termination is declared when the criteria

∥r^{k}∥ ≤ n max(TOL, 10^{−9}), μ_{k} ≤ TOL

are both satisfied, where TOL = 10^{−10} in our examples.

We experimented with three nonlinear problems from Dirkse [2]. We omit the details of these problems here, but refer the reader to [2] and the references therein for further details. Briefly, the problems are

nash: Nash equilibrium problem, with n = 10;

Table 1: Performance of the algorithm on small NCP test problems

Problem  nash  josephy  colvncp
Starting Point x^{0} = e x^{0} = 10e x^{0} = e x^{0} = 10e x^{0} = e x^{0} = 10e
μ^{0}  .16(+3) .21(+5) .10(+2) .73(+4) .59(+2) .28(+5)
∥r^{0}∥_{∞}  .30(+4) .18(+5) .12(+2) .75(+3) .65(+3) .31(+5)
Iterations  43  15  9  17  17  24
Solves  47  19  13  22  26  34
Trial steps  702  17  10  17  26  35
Fast steps  2  2  2  3  7  8

josephy: Four-variable problem due to Josephy;

colvncp: A convex programming problem, formulated as an NCP with n = 15.

Two starting points, both of which are distant from the solution, are used for each problem. They are

x^{0} = e, y^{0} = max(1, ∥f(x^{0})∥_{∞})e,

and

x^{0} = 10e, y^{0} = max(1, ∥f(x^{0})∥_{∞})e.

Table 1 summarizes the performance of our code on these problems. The entries for μ^{0} and ∥r^{0}∥_{∞} are self-explanatory; they indicate that our starting point was distant from the solution for each of the problems. The number of “iterations” equals the number of evaluations of the Jacobian Df(·) and also the number of matrix factorizations. The number of “solves” indicates the number of times the factors were used to compute a (safe or fast) step. The number of solves typically exceeds the number of iterations because both a safe step and a fast step are computed on some iterations. The “trial steps” entry is the total number of candidate step lengths α_{k} that were tried during the entire algorithm. This number is equal to the total number of evaluations of the function vector f.

The results in the table represent good performance of the method, with the exception of the Nash problem from the first starting point. The large number of unsuccessful trial values of α_{k} is due to the lack of sophistication of our Armijo line search, which would certainly be replaced by a safeguarded polynomial interpolation scheme in a more practical implementation of the algorithm. When χ is changed from .9 to a more conservative .5, the number of trial steps for this case decreases to 139, without affecting the number of factorizations.

Our second set of test problems is quite different from the first. They are large, mixed linear complementarity problems, in which, given an index set I ⊂ {1, ⋯, n} and its complement Ī = {1, ⋯, n}\I, we aim to find a vector pair (x, y) such that

y = Mx + q,  (x_{I}, y_{I}) ≥ 0,  x_{I}^{T} y_{I} = 0,  yĪ = 0,  (44)

where the coefficient matrix M is positive semidefinite. A few obvious modifications to the algorithm are required to take account of the mixed nature of the problem, and we do not discuss these here. We obtain the problems by reformulating the stagewise extended linear-quadratic programming (ELQP) problems considered by Rockafellar [7] and Zhu and Rockafellar [13]. The motivation for these stagewise problems comes from optimal control, and we refer the reader to the cited references for more details. Stagewise ELQPs are quite complicated to state. They are defined in terms of the quadratic function

J(u_{0}, ⋯, u_{N}; v_{1}, ⋯, v_{N+1}) = ∑_{i=1}^{N} [p_{i}^{T} u_{i} + q_{i}^{T} v_{i} + \frac{1}{2}u_{i}^{T} P_{i}u_{i} − \frac{1}{2}v_{i}^{T} Q_{i}v_{i} − v_{i}^{T} D_{i}u_{i}]

+p_{0}^{T} u_{0} + q_{N+1}^{T}v_{N+1} + \frac{1}{2}u_{0}^{T} P_{0}u_{0} − \frac{1}{2}v_{N+1}^{T}Q_{N+1}v_{N+1} − ∑_{i=1}^{N+1}[C_{i}^{T} v_{i} + c_{i}]^{T} x_{i−1},

where the intermediate variables (x_{0}, ⋯, x_{N}) satisfy

x_{0} = B_{0}u_{0} + b_{0},  x_{i} = A_{i}x_{i−1} + B_{i}u_{i} + b_{i},  i = 1, ⋯, N.

The matrices P_{i} and Q_{i} are all positive semidefinite. The problem is to find a saddle point of J(u_{0}, ⋯, u_{N}; v_{1}, ⋯, v_{N+1}) subject to u_{i} ∈ U_{i}, i = 0, ⋯, N, and v_{i} ∈ V_{i}, i = 1, ⋯, N + 1, where each U_{i} and V_{i} is a polyhedral subset of IR^{n_{i}} and IR^{m_{i}}, respectively. Associated with the ELQP is a primal problem

min_{(u_{0},···,u_{N})∈U_{0}×···×U_{N}} f(u_{0}, ⋯, u_{N}),

where

f(u_{0}, ⋯, u_{N}) = max_{(v_{1},···,v_{N+1})∈V_{1}×···×V_{N+1}} J(u_{0}, ⋯, u_{N}; v_{1}, ⋯, v_{N+1}),

and a dual problem

max_{(v_{1},···,v_{N+1})∈V_{1}×···×V_{N+1}} g(v_{1}, ⋯, v_{N+1}),

where

g(v_{1}, ⋯, v_{N+1}) = min_{(u_{0},···,u_{N})∈U_{0}×···×U_{N}} J(u_{0}, ⋯, u_{N}; v_{1}, ⋯, v_{N+1}).

By introducing explicit representations of the polyhedral sets U_{i} and V_{i}, and introducing intermediate and slack variables, we can formulate the ELQP as a mixed monotone linear complementarity problem of the form (44). Moreover, by a “stagewise” ordering of the variables, we can ensure that the coefficient matrix M is banded, where the bandwidth is independent of N. It is the bandedness that makes the complementarity formulation practical, since the time taken to factor and solve the linear system at each iteration of our algorithm is O(N), rather than the O(N^{3}) that would be obtained by a more naive formulation. In our code, the LAPACK band solve routines DGBTRF and DGBTRS [1] are used to solve this linear system.

We tested our algorithm on modifications of the random stagewise ELQPs that are described in Section 6 of Zhu and Rockafellar [13]. The problems are obtained by discretizing Table 2: Performance of the algorithm on smaller ELQP problems

n_{i} = 10, m_{i} = 10, N = 64, dimension = 2600

Data Set 1 3 5 7 9
μ_{0}  .19(+4) .11(+6) .15(+4) .43(+4) .18(+4)
Iterations  15  19  15  23  18
Solves  20  24  18  35  24
Trial steps  16  24  19  35  24
Fast steps  3  4  5  3  5
CPU time (sec) 26.1  33.3  25.4  41.2  31.5

continuous generalized optimal control problems, where the data for the continuous problems is generated randomly. The matrices P_{i} and Q_{i} are chosen to be diagonal and positive semidefinite, while the remaining matrices and vectors in the problem are dense. The polyhedra U_{i} and V_{i} are rectangles, that is, Cartesian products of intervals on the real line. In [13], the matrices Q_{i} and P_{i} are all strictly positive definite, but we modify them here by setting Q_{i} ≡ 0. As discussed in Rockafellar [7], this choice corresponds to “hard” constraints on the primal variables u_{i}; in fact, the primal problem above reduces to a quadratic program. The algorithms described by Rockafellar [8] and Zhu and Rockafellar [13] tend not to perform well in this important case, however. They lend themselves better to fully quadratic problems in which P_{i} and Q_{i} are all positive definite. On the other hand, our interior-point algorithm seems indifferent to this property.

The code for our algorithm does not take advantage of the linear nature of these problems, even though the Armijo line search could be replaced easily by an exact determination of the maximal α_{k}. Our aim is to demonstrate that the algorithm that we analyze in this paper is quite effective for both linear and nonlinear problems.

We use the code discussed in [13] to generate the data and starting points for our test problems. Besides resetting each Q_{i} to 0, we modify the starting points slightly to ensure strict interiority. We use the odd-numbered data sets from the problem generator, which accounts for our numbering scheme for the test problems. Tables 2 and 3 contain a summary of our results. The first line of each table contains the dimensions of the problem according to the notation above. The total dimension of the linear system to be solved at each iteration depends not only on the total number of components in (u_{0}, ⋯, u_{N}) and (v_{1}, ⋯, v_{N+1}) but also on the number of intermediate and slack variables. Some components of the linear system can be eliminated conveniently, leaving a subproblem whose size is the “dimension” indicated in Tables 2 and 3. The numbers of iterations, solves, and trial step lengths are as in Table 1. The last row contains CPU times on a Sun SPARCstation IPX.

The results indicate good performance, with between 15 and 28 iterations performed for each problem. These results could have been improved by finer tuning of the user-defined parameters, but they substantiate our claim that a single set of parameters can give good performance on very different problems (large and small, linear and nonlinear). Note that the average number of iterations is slightly higher in Table 3, reflecting the oft-made observation Table 3: Performance of the algorithm on larger ELQP problems

n_{i} = 20, m_{i} = 20, N = 128, dimension = 10320

Data Set 1 3 5 7 9
μ_{0}  .28(+4) .11(+6) .14(+5) .18(+5) .10(+5)
Iterations  21  26  20  28  22
Solves  29  38  25  45  32
Trial steps  34  48  31  52  36
Fast steps  6  4  6  1  4
CPU time (sec) 435.  548.  406.  592.  452.

that this number grows slowly with problem size in most interior-point methods. Note, too, that a number of fast steps were taken at the tail of the iteration sequence, and rapid local convergence was observed.

7 Conclusions

The local convergence results represent a natural extension of those presented in [11] for the case of linear f. In this case, we require only existence of a strictly complementary solution for ∥(Δx, Δy)∥ = O(μ) to hold. Monteiro and Wright [5] show that the strict complementarity assumption is necessary for superlinear convergence of methods that behave like Newton’s method near the solution. There is a considerable gap between the weak assumption of the linear case and the stronger assumptions of Section 5 for nonlinear f, and current research is aimed at bridging this gap as far as possible.

Our results can be extended easily to the case of Hölder continuous Df(·), that is, ∥Df(x^{1}) − Df(x^{2})∥ ≤ L∥x^{1} − x^{2}∥^{τ} for some τ ∈ (0, 1] and x^{1}, x^{2} in a neighborhood of x^{∗}. When τ > (\sqrt{5} − 1)/2, we can show local convergence with Q-order 1 + τ̂, provided that τ̂ lies in the range (0, τ^{2} + τ − 1). The algorithm can be modified to allow reduction of τ̂ after failed fast iterations, so that that the inclusion τ̂ ∈ (0, τ^{2} + τ − 1), and therefore superlinear convergence, eventually occurs. Knowledge of τ is not necessary. In the case of Df(·) merely continuous near x^{∗}, we still have global convergence, from the analysis of Section 3.

Acknowledgments

We are grateful to Ciyou Zhu for numerous discussions of ELQPs and for providing the code for the algorithm described in [13].

Appendix

Proof of Lemma 2.1. We deal only with the case of r^{0} ≠ 0, since the proofs for the case of r^{0} = 0 are simple specializations.

We prove (17) and (18) by induction on l. Taking l = 1, we know that safe steps are taken for k = 0, 1, ⋯, k_{1} − 1. Therefore, from (9b), we have

μ_{k_{1}}/μ_{0} ≥ ∏_{j=1}^{k_{1}−1}(1 − α_{j}) = ν_{k_{1}},

while, from (14),

r^{k_{1}} = ν_{k_{1}}r^{0}.

Therefore, by (11), we have β̃_{k_{1}} ≥ 1 and t̃_{k_{1}} = 0 = t_{k_{1}}, yielding (17). For (18), we note that by choice of (x^{0}, y^{0}) and (8), we have

γ̃_{0} = γ_{max}.

Since safe steps are taken for k = 0, 1, ⋯, k_{1} − 1, we have by repeated application of (9a) and (8) that

γ_{max} = γ̃_{0} = γ̃_{1} = ⋯ = γ̃_{k_{1}}.

Hence,

γ̃_{k_{1}} − γ_{min} = γ_{max} − γ_{min},

so (18) is satisfied for l = 1.

We now assume that (17) and (18) are satisfied for k_{1}, k_{2}, ⋯, k_{l} and prove that they continue to be true at k_{l+1}. From (11) and (14), we have that

μ_{k_{l}}/μ_{0} ≥ ν_{k_{l}} β̃_{k_{l}} ≥ ν_{k_{l}} ∏_{j=1}^{t̃_{k_{l}}}(1 − γ̄^{j}).

Applying (13b) to the fast step that is taken at iteration k_{l}, we have that

μ_{k_{l}+1}/μ_{0} ≥ (1 − α_{k_{l}})(1 − β̂_{k_{l}})(μ_{k_{l}}/μ_{0})
≥ (1 − α_{k_{l}})(1 − β̂_{k_{l}})ν_{k_{l}} β̃_{k_{l}}
= ν_{k_{l}+1}(1 − β̂_{k_{l}})β̃_{k_{l}}
= ν_{k_{l}+1}(1 − γ̄^{t̃_{k_{l}}+1})β̃_{k_{l}}
≥ ν_{k_{l}+1} ∏_{j=1}^{t̃_{k_{l}}+1}(1 − γ̄^{j}).

Since safe steps are taken at iterations k_{l} + 1, ⋯, k_{l+1} − 1, we have from (9b) that the ratio (μ_{k}/ν_{k}) is nondecreasing for k = k_{l} + 1, ⋯, k_{l+1}. In particular,

μ_{k_{l+1}}/μ_{0} ≥ ν_{k_{l+1}} ∏_{j=1}^{t̃_{k_{l}}+1}(1 − γ̄^{j}),

and therefore

β̃_{k_{l+1}} ≥ ∏_{j=1}^{t̃_{k_{l}}+1}(1 − γ̄^{j}).  (45)

From (45), the definition of t̃_{k_{l+1}}, and the inductive hypothesis for (17), we have

t̃_{k_{l+1}} ≤ t̃_{k_{l}} + 1 ≤ t_{k_{l}} + 1 = l = t_{k_{l+1}},

so (17) continues to hold at k_{l+1}.

Combining (13a) at k = k_{l} with (8) at k = k_{l} + 1, we have that γ̃_{k_{l}+1} ≥ γ̂_{k_{l}}. In fact, as observed earlier, γ̃_{k} is nondecreasing over the safe iterations k = k_{l} + 1, ⋯, k_{l+1} − 1, and so

γ̃_{k_{l+1}} ≥ ⋯ ≥ γ̃_{k_{l}+1} ≥ γ̂_{k_{l}}.

Hence, by using the definition of γ̂_{k} and the inductive hypothesis for (18), we have that

γ̃_{k_{l+1}} − γ_{min} ≥ γ̂_{k_{l}} − γ_{min} = γ̄(γ̃_{k_{l}} − γ_{min}) ≥ γ̄^{l}(γ_{max} − γ_{min}).

Therefore (18) continues to hold at iteration k_{l+1}. We conclude that (17) and (18) hold for all l = 1, 2, · · ·.

To prove (19), we note by (18) and the definition of γ̂_{k_{l}} that

γ̃_{k_{l}} − γ̂_{k_{l}} = γ̃_{k_{l}} − [γ_{min} + γ̄(γ̃_{k_{l}} − γ_{min})] = (1 − γ̄)(γ̃_{k_{l}} − γ_{min}) ≥ (1 − γ̄)γ̄^{l−1}(γ_{max} − γ_{min}),

as required.

For (20), we need only note that

β̂_{k_{l}} = γ̄^{t̃_{k_{l}}+1} ≥ γ̄^{(l−1)+1} = γ̄^{l}.

The final inequality (21) follows immediately from (19) and (20), if we note that 1−γ̄ ≥ γ̄.

Proof of Lemma 3.1. Define

δ̂ = \frac{1}{2} min_{i=1,···,n} (min(x̄_{i}, ȳ_{i})) > 0,

and choose positive constants C_{3} and C_{2} such that

x^{T} y ≥ C_{2} and ∥(x, y)∥ ≤ C_{3} for all (x, y) ∈ B((x̄, ȳ), δ̂).

Note that the coefficient matrix

[ \begin{matrix}Df(x) & −I \\ Y & X \end{matrix} ]

is nonsingular at (x̄, ȳ) and continuous in an open set containing B((x̄, ȳ), δ̂) ⊂ (IR_{+}^{n} × IR_{+}^{n}). The right-hand side of (3) is also continuous with respect to both (x, y) and σ̃. Therefore we can define a constant C_{1} > 0 such that

∥(Δx, Δy)∥ ≤ C_{1}

for all (x, y) ∈ B((x̄, ȳ), δ̂) and σ̃ ∈ [σ̄, 1/2].

Taking α ∈ [0, 1] and (x, y) ∈ B((x̄, ȳ), δ̂), we have

x(α)^{T} y(α)

= (x + αΔx)^{T} (y + αΔy + f(x + αΔx) − f(x) − αDf(x)Δx)
= x^{T} y − α(1 − σ̃)x^{T} y + α^{2}Δx^{T} Δy + α ∫_{0}^{1}(x + αΔx)^{T} [Df(x + θαΔx) − Df(x)]Δx dθ
≤ x^{T} y[1 − α(1 − σ̃)] + α^{2}C_{1}^{2} + αC_{1}(C_{3} + C_{1}) sup_{θ∈(0,1)} ∥Df(x + θαΔx) − Df(x)∥. (46)

Hence the left inequality in (9b) is certainly satisfied if

α^{2}C_{1}^{2} + αC_{1}(C_{3} + C_{1}) sup_{θ∈(0,1)} ∥Df(x + θαΔx) − Df(x)∥ ≤ (1 − κ)α(1 − σ̃)x^{T} y.  (47)

By continuity of Df, we can make the following assertion: for each (x, y) ∈ B((x̄, ȳ), δ̂), there is η_{1}(x) > 0 such that

∥Df(x + d) − Df(x)∥ ≤ \frac{(1 − κ)C_{2}}{4C_{1}(C_{3} + C_{1})}

for all d with ∥d∥ ≤ η_{1}(x). Since B((x̄, ȳ), δ̂) is compact, we have that

η̄_{1} =^{△} inf_{(x,y)∈B((x̄,ȳ),δ̂)∩Ω} η_{1}(x) > 0.

Therefore, setting

α̂^{(1)} = min (\frac{(1 − κ)C_{2}}{4C_{1}^{2}} , \frac{η̄_{1}}{C_{1}} ) ,

we have for (x, y) ∈ B((x̄, ȳ), δ̂), σ̃ ∈ [σ̄, 1/2], and α ∈ [0, α̂^{(1)}] that

α^{2}C_{1}^{2} + αC_{1}(C_{3} + C_{1}) sup_{θ∈(0,1)} ∥Df(x + θαΔx) − Df(x)∥

≤ αα̂^{(1)}C_{1}^{2} + αC_{1}(C_{3} + C_{1}) sup_{θ∈(0,1)} ∥Df(x + θαΔx) − Df(x)∥
≤ α(1 − κ)(C_{2}/4) + αC_{1}(C_{3} + C_{1})(1 − κ)(C_{2}/4)/(C_{1}(C_{3} + C_{1}))
≤ (1 − κ)α(C_{2}/2)
≤ (1 − κ)α(1 − σ̃)x^{T} y.  (48)

Hence (47) is satisfied, yielding the left inequality in (9b).

We now show that (9a) holds. For any i = 1, ⋯, n, (x, y) ∈ B((x̄, ȳ), δ̂), and α ∈ [0, α̂^{(1)}], we have

x_{i}(α)y_{i}(α)

= (x_{i} + αΔx_{i})(y_{i} + αΔy_{i} + f_{i}(x + αΔx) − f_{i}(x) − αDf_{i}(x)Δx)
≥ x_{i}y_{i} + α[−x_{i}y_{i} + σ̃(x^{T} y)/n] − α^{2}C_{1}^{2} − αC_{1}(C_{3} + C_{1}) sup_{θ∈(0,1)} ∥Df_{i}(x + θαΔx) − Df_{i}(x)∥
≥ x_{i}y_{i}(1 − α) + σ̃α(x^{T} y)/n − α^{2}C_{1}^{2} − αC_{1}(C_{3} + C_{1}) sup_{θ∈(0,1)} ∥Df(x + θαΔx) − Df(x)∥
≥ γ̃(x^{T} y/n)(1 − α) + σ̃α(x^{T} y/n) − α^{2}C_{1}^{2}  (49)

−αC_{1}(C_{3} + C_{1}) sup_{θ∈(0,1)} ∥Df(x + θαΔx) − Df(x)∥.

Meanwhile, we have as in (46) that

x(α)^{T} y(α) ≤ x^{T} y[1 − α(1 − σ̃)] + α^{2}C_{1}^{2} + αC_{1}(C_{3} + C_{1}) sup_{θ∈(0,1)} ∥Df(x + θαΔx) − Df(x)∥.

By combining the last two expressions, we find that (9a) will be satisfied if

(x^{T} y/n)ασ̃(1 − γ̃) − 2α^{2}C_{1}^{2} − 2αC_{1}(C_{3} + C_{1}) sup_{θ∈(0,1)} ∥Df(x + θαΔx) − Df(x)∥ ≥ 0,

which in turn is true if

(C_{2}/n)σ̄(1 − γ_{max}) − 2αC_{1}^{2} − 2C_{1}(C_{3} + C_{1}) sup_{θ∈(0,1)} ∥Df(x + θαΔx) − Df(x)∥ ≥ 0. (50)

Just as we chose η̄_{1} above, we can find η̄_{2} such that

∥Df(x + d) − Df(x)∥ ≤ \frac{(C_{2}/n)σ̄(1 − γ_{max})}{4C_{1}(C_{3} + C_{1})}

for all (x, y) ∈ B((x̄, ȳ), δ̂) and ∥d∥ ≤ η̄_{2}. Hence, if we choose

α̂^{(2)} = min (α̂^{(1)}, \frac{η̄_{2}}{C_{1}} , \frac{(C_{2}/n)σ̄(1 − 2γ̄)}{4C_{1}^{2}}) ,

a calculation similar to (48) shows that (50) (and hence (9a)) is satisfied for all α ∈ [0, α̂^{(2)}] and (x, y) ∈ B((x̄, ȳ), δ̂).

For the remaining inequality, we have as in (46) that

x(α)^{T} y(α) ≥ x^{T} y[1 − α(1 − σ̃)] − α^{2}C_{1}^{2} − αC_{1}(C_{3} + C_{1}) sup_{θ∈(0,1)} ∥Df(x + θαΔx) − Df(x)∥.

Therefore, the right inequality in (9b) holds if

ασ̃x^{T} y − α^{2}C_{1}^{2} − αC_{1}(C_{3} + C_{1}) sup_{θ∈(0,1)} ∥Df(x + θαΔx) − Df(x)∥ ≥ 0,

which in turn is satisfied if

σ̄C_{2} − αC_{1}^{2} − C_{1}(C_{3} + C_{1}) sup_{θ∈(0,1)} ∥Df(x + θαΔx) − Df(x)∥ ≥ 0.  (51)

As before, we can choose η̄_{3} > 0 such that

∥Df(x + d) − Df(x)∥ ≤ \frac{σ̄C_{2}}{2C_{1}(C_{3} + C_{1})}

for all (x, y) ∈ B((x̄, ȳ), δ̂) and ∥d∥ ≤ η̄_{3}. Another calculation like (48) shows that the inequality (51) is satisfied if α ∈ [0, α̂^{(3)}], where

α̂^{(3)} = min (α̂^{(2)}, \frac{σ̄C_{2}}{2C_{1}^{2}} , \frac{η̄_{3}}{C_{1}} ) .

The conclusion of the lemma is obtained by setting α̂ = α̂^{(3)}.

References

[1] E. Anderson, Z. Bai, C. Bischof, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, S. Ostrouchov, and D. Sorensen, LAPACK User’s Guide, SIAM, Philadelphia, 1992.

[2] S. P. Dirkse, Solving complementarity problems via GAMS: A model library. Manuscript, May 1993.

[3] O. Güler, Existence of interior points and interior paths in nonlinear monotone complementarity problems, Mathematics of Operations Research, 18 (1993), pp. 128–147.

[4] M. Kojima, T. Noma, and A. Yoshise, Global convergence and detecting infeasibility in interior-point algorithms, Research Reports on Information Sciences, Series B: Operations Research B–257, Department of Information Sciences, Tokyo Institute of Technology, Tokyo, Japan, September 1992.

[5] R. D. C. Monteiro and S. J. Wright, Local convergence of interior-point algorithms for degenerate monotone LCP, Computational Optimization and Applications, 3 (1994), pp. 131–155.

[6] F. A. Potra and Y. Ye, Interior point methods for nonlinear complementarity problems, Technical Report 15, Department of Mathematics, University of Iowa, Iowa City, July 1991.

[7] R. T. Rockafellar, Linear-quadratic programming and optimal control, SIAM Journal on Control and Optimization, 25 (1987), pp. 781–814.

[8] , Computational schemes for large-scale problems in extended linear-quadratic programming, Mathematical Programming, 48 (1990), pp. 447–474.

[9] S. J. Wright, An infeasible-interior-point algorithm for linear complementarity problems, Preprint MCS–P331–1092, Mathematics and Computer Science Division, Argonne National Laboratory, Argonne, Ill., October 1992. To appear in Mathematical Programming, Series A.

[10] , A path-following infeasible-interior-point algorithm for linear complementarity problems, Optimization Methods and Software, 2 (1993), pp. 79–106.

[11] , A path-following interior-point algorithm for linear and quadratic optimization problems, Preprint MCS–P401–1293, Mathematics and Computer Science Division, Argonne National Laboratory, Argonne, Ill., December 1993.

[12] S. J. Wright and Y. Zhang, A superquadratic infeasible-interior-point method for linear complementarity problems, Preprint MCS–P418–0294, Mathematics and Computer Science Division, Argonne National Laboratory, Argonne, Ill., February 1994.

[13] C. Zhu and R. T. Rockafellar, Primal-dual projected gradient algorithms for extended linear-quadratic programming, SIAM Journal on Optimization, 3 (1993), pp. 751– 783.