A Path-Following Interior-Point Algorithm for Linear and Quadratic Problems^{∗}

Stephen J. Wright^{†}

December 1993 (revised January, 1995)

PREPRINT MCS-P401-1293, MCS DIVISION, ARGONNE NATIONAL LABORATORY

Abstract

We describe an algorithm for the monotone linear complementarity problem (LCP) that converges from any positive, not necessarily feasible, starting point and exhibits polynomial complexity if some additional assumptions are made on the starting point. If the problem has a strictly complementary solution, the method converges subquadratically. We show that the algorithm and its convergence properties extend readily to the mixed monotone linear complementarity problem and, hence, to all the usual formulations of the linear programming and convex quadratic programming problems.

1 Introduction

The monotone linear complementarity problem (LCP) is to find a vector pair (x, y) ∈ IR^{n}×IR^{n} such that

y = Mx + q,  (x, y) ≥ 0,  x^{T} y = 0,  (1)

where q ∈ IR^{n} and M is an n × n positive semidefinite (p.s.d.) matrix. The mixed monotone linear complementarity problem (MLCP) is to find a vector triple (x, y, z) ∈ IR^{n} × IR^{n} × IR^{m} such that

[ \begin{matrix}y \\ 0 \end{matrix} ] = [ \begin{matrix}M_{11} & M_{12} \\ M_{21} & M_{22} \end{matrix} ] [ \begin{matrix}x \\ z \end{matrix} ] + [ \begin{matrix}q_{1} \\ q_{2} \end{matrix} ]  (2a)
(x, y) ≥ 0,  x^{T} y = 0,  (2b)

where

M = [ \begin{matrix}M_{11} & M_{12} \\ M_{21} & M_{22} \end{matrix} ]

∗This research was supported by the Office of Scientific Computing, U.S. Department of Energy, under Contract W-31-109-Eng-38.

†Mathematics and Computer Science Division, Argonne National Laboratory, 9700 South Cass Avenue, Argonne, Illinois 60439.

is p.s.d. All conventional formulations of the linear programming (LP) and convex quadratic programming (QP) problems can be posed in the form (2) by writing out their conditions for optimality. For instance, consider the QP problem given by

min_{w∈IR^{p}} \frac{1}{2}w^{T} Qw + c^{T} w  (3)

subject to

w_{i} ≥ l_{i} (i ∈ L ⊂ {1, ⋯, p}), w_{i} ≤ u_{i} (i ∈ U ⊂ {1, ⋯, p}),
Cw ≥ d, Aw = b,  (4)

where Q is symmetric p.s.d., C ∈ IR^{m_{I}×p}, A ∈ IR^{m_{E}×p}, and so on. If we define

E_{L} = [e_{i}^{T}]_{i∈L}, E_{U} = [e_{i}^{T}]_{i∈U},

where e_{i} is the i-th unit vector from the standard basis, and

l = [l_{i}]_{i∈L}, u = [u_{i}]_{i∈U},

then we can state the optimality conditions for (3),(4) in the form (2) by defining

n = |L| + |U| + m_{I}, m = p + m_{E},

and

M_{22} = [ \begin{matrix}Q & −A^{T} \\ A & 0 \end{matrix} ] ,  M_{12} = [ \begin{matrix}E_{L} & 0 \\ −E_{U} & 0 \\ C & 0 \end{matrix} ] ,  M_{21} = [ \begin{matrix}−E_{L}^{T} & E_{U}^{T} & −C^{T} \\ 0 & 0 & 0 \end{matrix} ] ,

M_{11} = 0,  q_{1} = [ \begin{matrix}−l \\ u \\ −d \end{matrix} ] ,  q_{2} = [ \begin{matrix}c \\ −b \end{matrix} ] .

In this paper, we focus on an algorithm for (1) and its convergence properties. We then show, using recent work involving the relationship between problems (1) and (2), that this algorithm can be extended painlessly to (2) and, hence, to all the usual LP and QP formulations. Little loss of efficiency is involved in solving LPs and QPs by embedding them in algorithms for (2), provided the linear algebra takes account of the particular structure of each problem. Hence we feel that the linear complementarity formulation is the best one to consider because of its generality, simplicity of notation, and practical efficiency of the algorithms on all its special cases.

In two recent papers [10, 9], we described algorithms for (1) that are globally convergent, have polynomial complexity when the starting point (x^{0}, y^{0}) satisfies certain assumptions, and exhibit subquadratic convergence of the complementarity gap μ_{k} = (x^{k})^{T} y^{k}/n to zero. Neither the starting point nor the iterates are feasible in general. In both algorithms, most of the work at each iteration consists of a single matrix factorization and between one and three triangular solves with the computed factors. These methods were the first interior-point methods with this desirable combination of properties.

The local analysis in both [10] and [9] requires existence of a strictly complementary solution, that is, (x^{∗}, y^{∗}) solving (1) such that x^{∗} + y^{∗} > 0. This assumption is always satisfied if the LCP is a reformulated linear programming problem and, as shown in Monteiro and Wright [7], it is necessary for superlinear convergence of Newton-based primal-dual algorithms. The earlier paper [10] makes an additional assumption: existence of a strictly feasible point, that is, (x̄, ȳ) such that ȳ = Mx̄+q, (x̄, ȳ) > 0. This assumption is undesirable because it is usually not satisfied by large practical problems.

In this paper, we describe an algorithm that is quite similar to the one in [10], except that it does not use the clumsy merit function ϕ(x, y) = x^{T} y +∥y −Mx−q∥ and allows more flexibility in the choice of parameters and starting point. As we mention at the end of Section 3, we can also allow more flexibility in the choice of steplength, bringing the algorithm close to current computational practice. The analysis here is considerably stronger than in [10]. The strict feasibility assumption is no longer required, technical arguments are streamlined, and R-subquadratic convergence of the iterates (x^{k}, y^{k}) to a strictly complementary solution is proved.

The new algorithm is quite different from the one in [9], where the ratio of complementarity gap μ to infeasibility norm ∥r∥ is kept constant until the “fast” phase, when small variations are allowed. To achieve this effect, a correction to the basic path-following step is computed (at a cost of one additional back-substitution) and a messy planar search procedure is performed. The complications in [9] made it possible to prove boundedness of the iteration sequence {(x^{k}, y^{k})}, which was a vital element in the asymptotic convergence analysis. Our use of a stronger technical result (Lemma 4.3) removes the dependence on boundedness. Instead, as mentioned above, we prove convergence of the iteration sequence as a consequence of the main results.

Our algorithm is specified in Section 2. In Section 3 we prove global linear convergence and polynomial complexity. Some technical results are proved in Section 4; these are used to prove superlinear convergence in Section 5. Section 6 shows that the algorithm and its convergence properties can be extended to the mixed problem (2) because (2) can be reformulated as (1). We stress at the outset that this reformulation need not be performed explicitly; it suffices to observe that the (x, y) iterates generated by our extended algorithm are the same as those that would be obtained by reformulating the problem as (1) and applying the algorithm of Section 2 directly, except possibly for some swapping of components to be discussed later.

Unless otherwise specified, ∥ · ∥ denotes the Euclidean norm of a vector. Iteration numbers appear as superscripts on vectors and matrices and as subscripts on scalars. To avoid notational clutter in Sections 3, 4, and 5, we drop the iteration index k from vector and matrix quantities in the proofs. It is retained explicitly in the statement of each result.

We denote the solution set and strictly complementary solution set by

S = {(x^{∗}, y^{∗}) | (x^{∗}, y^{∗}) solves (1)} , S^{c} = {(x^{∗}, y^{∗}) ∈ S | x^{∗} + y^{∗} > 0} ,

respectively. The range space of a matrix is denoted by R(·).

2 The Algorithm

The algorithm generates a sequence of strictly positive iterates (x^{k}, y^{k}). To describe the step between successive iterates, we define

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

We refer to μ_{k} as the complementarity gap and to r^{k} as the residual. Each step is calculated as follows.

Given γ̃ ∈ (0, 1), β̃ ∈ [0, 1), σ̃ ∈ [0, 1), solve

[ \begin{matrix}M & −I \\ Y^{k} & X^{k} \end{matrix} ] [ \begin{matrix}u^{k} \\ v^{k} \end{matrix} ] = [ \begin{matrix}r^{k} \\ −X^{k}Y^{k}e + σ̃μ_{k}e \end{matrix} ] .  (5)

Choose

α̃ = arg min_{α∈[0,α̂]} μ_{k}(α) =^{△} (x^{k} + αu^{k})^{T} (y^{k} + αv^{k})/n,  (6)

where α̂ is the largest number in [0, 1] such that the following inequalities are satisfied for all α ∈ [0, α̂):

(x^{k} + αu^{k})^{T} (y^{k} + αv^{k}) ≥ (1 − β̃)(1 − α)(x^{k})^{T} y^{k},  if r^{k} ≠ 0,  (7a)
(x_{i}^{k} + αu_{i}^{k})(y_{i}^{k} + αv_{i}^{k}) ≥ (γ̃/n)(x^{k} + αu^{k})^{T} (y^{k} + αv^{k}),  i = 1, ⋯, n. (7b)

The search direction obtained from (5) is simply the Newton step for the system of nonlinear equations

F (x, y) = [ \begin{matrix}Mx − y + q \\ XY e \end{matrix} ] = [ \begin{matrix}0 \\ σ̃μ_{k}e \end{matrix} ] ,

from the point (x^{k}, y^{k}).

The inequality (7b), usually referred to as a centering condition, ensures that the iterates do not approach the boundary of the nonnegative orthant too closely. Because we restrict γ̃ to the range [γ_{min}, γ_{max}] for 0 < γ_{min} < γ_{max} ≤ 1/2, (7b) implies that

x_{i}y_{i} ≥ γ_{k}μ_{k} ≥ γ_{min}μ_{k},  ∀k.  (8)

The condition (7a) is used to ensure that improvement in the complementarity gap μ_{k} does not outstrip improvement in the infeasibility ∥r^{k}∥ by too much; a vector pair (x, y) that is complementary but not feasible is of no interest. Note that we need not enforce condition (7a) if the current point is already feasible.

We can now state our algorithm.

Given γ̄ ∈ (0, 1/2), γ_{min} and γ_{max} with 0 < γ_{min} < γ_{max} ≤ 1/2, σ̄ ∈ (0, 1/2), ρ ∈ (0, γ̄), μ̄ ∈ (0, 1], and (x^{0}, y^{0}) with x_{i}^{0}y_{i}^{0} ≥ γ_{max}μ_{0} > 0;

t_{0} ← 1, γ_{0} ← γ_{max};

for k = 0, 1, 2, ⋯

if  μ_{k} = 0  then stop;

(* attempt a fast step *) Solve (5)–(7) with σ̃ = 0, β̃ = γ̄^{t_{k}}, γ̃ = γ_{min} + γ̄^{t_{k}}(γ_{max} − γ_{min}); if (x^{k} + α̃u^{k})^{T} (y^{k} + α̃v^{k})/n ≤ ρμ_{k} then α_{k} ← α̃, β_{k} ← β̃, σ_{k} ← σ̃, γ_{k+1} ← γ̃;

t_{k+1} ← t_{k} + 1;

(x^{k+1}, y^{k+1}) ← (x^{k}, y^{k}) + α_{k}(u^{k}, v^{k});

go to next k; end if

(* revert to a safe step *) Solve (5)–(7) with σ̃ ∈ [σ̄, 1/2], β̃ = 0, γ̃ = γ_{k}; α_{k} ← α̃, β_{k} ← 0, σ_{k} ← σ̃, γ_{k+1} ← γ̃; t_{k+1} ← t_{k}; (x^{k+1}, y^{k+1}) ← (x^{k}, y^{k}) + α_{k}(u^{k}, v^{k}); go to next k;

end for.

The algorithm can be motivated in a few sentences. We begin each major iteration by trying to take a fast step, which uses an affine scaling search direction. To encourage longer steps to be taken we use a strictly positive value of β̃ and a value γ̃ smaller than the current γ_{k}. The fast steps are accepted only if they produce a reduction in μ_{k} of at least a factor of ρ. Otherwise, the algorithm reverts to taking a safe step, whose major distinguishing feature is its use of a strictly positive value σ̃ ≥ σ̄ of the centering parameter.

Safe steps tend to be taken on early iterations, while fast steps are taken toward the tail of the sequence. There may be a gray area in which both safe and fast steps are taken.

The algorithm can be modified to try fast steps only when there is some reasonable hope that they will be accepted. (Earlier versions of the algorithm used a threshold criterion μ_{k} ≤ μ̄, with μ̄ a user-defined parameter, to decide whether to calculate the safe step.) For the sake of simplicity, we do not consider such modifications here, but note simply that the superlinear convergence properties of the algorithm will hold provided that the fast step is eventually tried on every iteration. Besides omitting the threshold μ̄, the algorithm above differs from the one described in [10] in that the duality gap μ is used directly in place of the merit function ϕ, and the particular choices γ_{min} = γ̄ and γ_{max} = 2γ̄ are relaxed.

3 Global Convergence and Polynomial Complexity

In this section, we show that the algorithm converges globally to the solution set of (1) from any starting point (x^{0}, y^{0}) > 0. When the algorithm is initialized in a certain way, the number of iterations is quadratic in the problem dimension n. Throughout the section, we make the following assumption.

Assumption 1 The LCP (1) is feasible; that is, there is a pair (x, y) such that y = Mx + q and (x, y) ≥ 0.

Assumption 1 implies that S ≠ ∅ (see, for example, [1, Theorem 3.1.2]).

If we define the monotonically decreasing sequence {ν_{k}} by

ν_{0} = 1, ν_{k+1} = (1 − α_{k})ν_{k},

it is easy to see that r^{k} = ν_{k}r^{0}. We have the following simple result, whose proof follows that of [9, Lemma 3.1].

Lemma 3.1 The constant β̂ defined by

β̂ =^{△} ∏_{k=0}^{∞}(1 − β_{k}) ≥ ∏_{k=1}^{∞}(1 − γ̄^{k})

is strictly positive, and we have for all k ≥ 0 that

μ_{k} ≥ β̂ν_{k}μ_{0}.  (9)

We can also show that all iterates remain strictly positive, except when finite termination occurs.

Lemma 3.2 For all iterates generated by the algorithm, we have either (x^{k}, y^{k}) > 0 or μ_{k} = 0.

Proof. We prove the result by induction. Note first that the assertion is trivially satisfied by the initial iterate (x^{0}, y^{0}) > 0. If μ_{k} = 0, the algorithm terminates at the k-th iterate. For the remainder of the proof, we assume that (x^{k}, y^{k}) > 0 and prove that either (x^{k+1}, y^{k+1}) > 0 or μ_{k+1} = 0. We consider the cases r^{k} ≠ 0 and r^{k} = 0 separately.

If r^{k} ≠ 0, the constraint (7a) is applied to the choice of α_{k}. Hence, combining (7a) and (7b), we have that

(x_{i}^{k} + αu_{i}^{k})(y_{i}^{k} + αv_{i}^{k}) ≥ (γ̃/n)(1 − β̃)(1 − α)(x^{k}^{T} y^{k}), ∀α ∈ [0, α_{k}).

Since α_{k} ≤ 1, γ̃ > 0, β̃ ∈ [0, 1), and x^{k}^{T} y^{k} > 0, the right-hand side of this expression is strictly positive for all α ∈ [0, α_{k}). Since (x_{i}^{k}, y_{i}^{k}) > 0, it follows that (x_{i}^{k} +αu_{i}^{k}, y_{i}^{k} +αv_{i}^{k}) > 0 for all α ∈ [0, α_{k}). If x_{i}^{k} + α_{k}u_{i}^{k} = 0 or y_{i}^{k} + α_{k}v_{i}^{k} = 0 for some index i, we have from (7b) that

μ_{k+1} = (x^{k} + α_{k}u^{k})^{T} (y^{k} + α_{k}v^{k})/n = 0.

The alternative case is

(x^{k+1}, y^{k+1}) = (x^{k} + α_{k}u^{k}, y^{k} + α_{k}v^{k}) > 0,

so our claim holds.

Consider now the remaining case r^{k} = 0, for which (7a) is not enforced. Suppose first that there exist α values in the range [0, α_{k}] such that

x_{i}^{k} + αu_{i}^{k} = 0 or y_{i}^{k} + αv_{i}^{k} = 0,  (10)

for some index i = 1, ⋯, n, and let ᾱ_{k} be the smallest of these values. It follows from (7b) that

(x^{k} + ᾱ_{k}u^{k})^{T} (y^{k} + ᾱ_{k}v^{k}) = 0.  (11)

Since r^{k} = 0, equation (5) implies that v^{k} = Mu^{k}. Hence, by positive semidefiniteness of M, we have u^{k}^{T} v^{k} ≥ 0. Using the second part of (5), we find that

(x^{k} + αu^{k})^{T} (y^{k} + αv^{k}) = (x^{k}^{T} y^{k})(1 − α(1 − σ̃)) + α^{2}u^{k}^{T} v^{k} ≥ (x^{k}^{T} y^{k})(1 − α(1 − σ̃)). (12)

Now since ᾱ_{k} ∈ [0, α_{k}] ⊂ [0, 1], the relations (11) and (12) can be satisfied simultaneously for α = ᾱ_{k} only if σ̃ = 0 and ᾱ_{k} = α_{k} = 1. Hence we are left with two possibilities. Either ᾱ_{k} = α_{k} = 1 and μ_{k+1} = 0, or there are no α ∈ [0, α_{k}] with the property (10), so (x^{k+1}, y^{k+1}) > 0. Therefore our claim holds again for the case of r^{k} = 0, and we are done.

Finite termination of the algorithm with μ_{k} = 0 and r^{k} = 0 is, of course, the simple case. Throughout the remainder of the paper, we make the implicit assumption that finite termination does not occur and that the algorithm generates an infinite sequence of iterates {(x^{k}, y^{k})}, k = 0, 1, · · ·.

The next lemma contains an inequality that is used in a number of places in the analysis. Similar results appear in Potra [8, Lemma 4.1], Mizuno [5, Lemma 3.3], and Wright [9, Lemma 3.2].

Lemma 3.3 Let (x^{∗}, y^{∗}) ∈ S. Then for all k ≥ 0 there is a constant C_{1} such that

∥y^{k}∥_{1} + ∥x^{k}∥_{1} ≤ C_{1} [μ_{0} + μ_{k}/ν_{k} + ∥y^{∗}∥_{1} + ∥x^{∗}∥_{1}] .  (13)

When the initial point satisfies

x^{0} = ξ_{x}e,  y^{0} = ξ_{y}e,  (14)

for some positive ξ_{x} and ξ_{y}, we have

ξ_{x}∥y^{k}∥_{1} + ξ_{y}∥x^{k}∥_{1} ≤ nμ_{0} + nμ_{k}/ν_{k} + ξ_{x}∥y^{∗}∥_{1} + ξ_{y}∥x^{∗}∥_{1}.  (15)

Proof. As in the proof of [9, Lemma 3.2] we obtain the inequalities

0 ≤ ν_{k}^{2}(x^{0})^{T} y^{0} + (1 − ν_{k})^{2}(x^{∗})^{T} y^{∗} + x^{T} y + ν_{k}(1 − ν_{k})((x^{0})^{T} y^{∗} + (x^{∗})^{T} y^{0})
−ν_{k}((x^{0})^{T} y + (y^{0})^{T} x) − (1 − ν_{k})((x^{∗})^{T} y + x^{T} y^{∗})  (16)

and

ν_{k}((x^{0})^{T} y + (y^{0})^{T} x) ≤ ν_{k}^{2}(x^{0})^{T} y^{0} + x^{T} y + ν_{k}(1 − ν_{k})((x^{0})^{T} y^{∗} + (y^{0})^{T} x^{∗}).  (17)

By defining

C̄_{1} = min_{i=1,···,n} min(x_{i}^{0}, y_{i}^{0}),

and noting that 1 − ν_{k} ≤ 1, we have

∥y∥_{1} + ∥x∥_{1} ≤ C̄_{1}^{−1} [ν_{k}nμ_{0} + nμ_{k}/ν_{k} + ∥x^{0}∥_{∞}∥y^{∗}∥_{1} + ∥y^{0}∥_{∞}∥x^{∗}∥_{1}] ,

which implies (13) for appropriately defined C_{1}. The other inequality (15) follows trivially from (14) and (17).

For purposes of polynomial complexity, we assume that the initial point is defined by (14), where ξ_{x} and ξ_{y} satisfy the following assumptions:

∥x^{∗}∥_{∞} ≤ ξ_{x},  ∥y^{∗}∥_{∞} ≤ ξ_{y},  ξ_{y} ≥ ∥q∥_{∞},  ξ_{y} ≥ ∥Me∥_{∞}ξ_{x} = ∥Mx^{0}∥_{∞}.  (18)

We also find the diagonal matrix D^{k} defined by

D^{k} = (X^{k})^{−1/2}(Y^{k})^{1/2}

useful in the subsequent analysis. The next lemma allows us to bound quantities involving the steps u^{k} and v^{k}.

Lemma 3.4 For all k ≥ 0, there is a constant ω such that

∥D^{k}u^{k}∥^{2} + ∥(D^{k})^{−1}v^{k}∥^{2} ≤ ωμ_{k}.  (19)

If the initial step is chosen according to (14) and (18), then

ω = \frac{2n^{2}}{γ_{min}} (5 + \frac{9}{β̂})^{2} ,  (20)

and so ω = O(n^{2}).

Proof. It is easy to see that the step (u, v) can be partitioned as

(u, v) = (ū, v̄) + (û, v̂),

where

[ \begin{matrix}M & −I \\ Y & X \end{matrix} ] [ ū ; v̄ ] = [ \begin{matrix}0 \\ −XY e + σ_{k}μ_{k}e \end{matrix} ]  (21)

and

[ \begin{matrix}M & −I \\ Y & X \end{matrix} ] [ û ; v̂ ] = [ \begin{matrix}r \\ 0 \end{matrix} ] .  (22)

(Because of nonsingularity of the coefficient matrix in (21) and (22), both (ū, v̄) and (û, v̂) are well defined.) As in [9, Lemma 3.3], we can show that

∥Dū∥ ≤ \frac{2nμ_{k}^{1/2}}{γ_{min}^{1/2}} ,  ∥D^{−1}v̄∥ ≤ \frac{2nμ_{k}^{1/2}}{γ_{min}^{1/2}} .  (23)

For the other component (û, v̂), we have from the second part of (22) that

Dû + D^{−1}v̂ = 0 ⇒ v̂ = −D^{2}û.

Hence from the first part of (22) and positive semidefiniteness of M, we have

Mû − v̂ = r ⇒ (M + D^{2})û = r ⇒ û^{T} D^{2}û ≤ û^{T} r.

Therefore

∥Dû∥^{2} ≤ ∥Dû∥∥D^{−1}r∥.

Using (8) together with Lemmas 3.1 and 3.3, we have that

∥Dû∥ ≤ ∥D^{−1}r∥ ≤ ∑_{i} (\frac{x_{i}}{y_{i}} )^{1/2} ∥r∥_{∞}

= ∑_{i}(x_{i}y_{i})^{−1/2}x_{i}ν_{k}∥r^{0}∥_{∞}
≤ \frac{∥r^{0}∥_{∞}}{γ_{min}^{1/2}μ_{k}^{1/2}} ν_{k}∥x∥_{1}
≤ \frac{∥r^{0}∥_{∞}}{γ_{min}^{1/2}μ_{k}^{1/2}} C_{1} [ν_{k}μ_{0} + μ_{k} + ν_{k}∥y^{∗}∥_{1} + ν_{k}∥x^{∗}∥_{1}]
≤ \frac{∥r^{0}∥_{∞}}{γ_{min}^{1/2}μ_{k}^{1/2}} C_{1}μ_{k} [1 + (μ_{0} + ∥y^{∗}∥_{1} + ∥x^{∗}∥_{1})/(β̂μ_{0})]
= C_{2}μ_{k}^{1/2},  (24)

for some appropriately defined constant C_{2}. Since ∥Dû∥ = ∥D^{−1}v̂∥, we therefore have

∥D^{−1}v̂∥ = ∥Dû∥ ≤ C_{2}μ_{k}^{1/2}.  (25)

By combining (23) and (25), we obtain

∥Du∥^{2} + ∥D^{−1}v∥^{2}
≤ (∥Dū∥ + ∥Dû∥)^{2} + (∥D^{−1}v̄∥ + ∥D^{−1}v̂∥)^{2}
≤ 2 (\frac{2nμ_{k}^{1/2}}{γ_{min}^{1/2}} + C_{2}μ_{k}^{1/2})^{2}
≤ ωμ_{k},

where ω is defined in an obvious way.

The special result (20) is proved by analysis similar to that of Lemma 3.4 in Wright [9].

The key theorem of this section shows that there is a uniform lower bound on the step length α_{k} on each safe step.

Theorem 3.5 If a safe step is taken at iteration k, then

α_{k} ≥ \frac{σ̄(1 − γ_{max})}{2ω},

where ω is as defined in Lemma 3.4.

Proof. We prove the result by showing that the conditions (7a) and (7b) hold for all α in the range

[0, \frac{σ̄(1 − γ_{max})}{2ω}] .  (26)

We further show that the complementarity gap μ_{k}(α) defined in (6) is decreasing on the interval (26). These observations are sufficient to prove the result.

Because of (5), we have that (7a) is satisfied if

(x + αu)^{T} (y + αv) ≥ (1 − α)x^{T} y
⇔ x^{T} y(1 − α + ασ_{k}) + α^{2}u^{T} v ≥ (1 − α)x^{T} y
⇔ σ_{k}x^{T} y + αu^{T} v ≥ 0.  (27)

From Lemma 3.4 we have

|u^{T} v| ≤ ∥Du∥∥D^{−1}v∥ ≤ ωμ_{k}.  (28)

Since σ_{k} ≥ σ̄ and α is in the range (26), we have

σ_{k}x^{T} y + αu^{T} v ≥ σ̄x^{T} y − α|u^{T} v|
≥ σ̄x^{T} y − \frac{σ̄(1 − γ_{max})}{2ω}ωμ_{k}
= σ̄x^{T} y [1 − \frac{1 − γ_{max}}{2n} ]
> 0.

Therefore (27) and hence (7a) are satisfied for α in the range (26).

By using (5), we note that (7b) is satisfied if

(x_{i} + αu_{i})(y_{i} + αv_{i}) ≥ γ_{k}(x + αu)^{T} (y + αv)/n
⇔ x_{i}y_{i}(1 − α) + ασ_{k}μ_{k} + α^{2}u_{i}v_{i} ≥ (γ_{k}/n) [x^{T} y(1 − α + σ_{k}α) + α^{2}u^{T} v] . (29)

Because of (8), (29) is true provided that

γ_{k}(1 − α)μ_{k} + ασ_{k}μ_{k} + α^{2}u_{i}v_{i} ≥ γ_{k} [(1 − α + σ_{k}α)μ_{k} + α^{2}u^{T} v/n]
⇔ σ_{k}(1 − γ_{k})μ_{k} + α(u_{i}v_{i} − γ_{k}u^{T} v/n) ≥ 0.

Now from (28) and

|u_{i}v_{i}| = |D_{ii}u_{i}||D_{ii}^{−1}v_{i}| ≤ ∥Du∥∥D^{−1}v∥ ≤ ωμ_{k},

we have

σ_{k}(1 − γ_{k})μ_{k} + α(u_{i}v_{i} − γ_{k}u^{T} v/n) ≥ σ_{k}(1 − γ_{k})μ_{k} − α(ωμ_{k} + ωμ_{k}/n) ≥ σ̄(1 − γ_{k})μ_{k} − 2αωμ_{k}.

Since α lies in the range (26) and γ_{k} ∈ [γ_{min}, γ_{max}], we have

σ̄(1 − γ_{k})μ_{k} − 2αωμ_{k} ≥ σ̄(1 − γ_{k})μ_{k} − 2ω \frac{σ̄(1 − γ_{max})}{2ω}μ_{k} ≥ 0.

Therefore (29) and hence (7b) hold for the interval in question.

Finally, from (6), (28), and (26), we have

nμ′_{k}(α) = (x^{T} v + y^{T} u) + 2αu^{T} v
≤ (σ_{k} − 1)x^{T} y + 2αωμ_{k}
≤ (σ_{k} − 1)nμ_{k} + σ̄(1 − γ_{max})μ_{k}
≤ (σ_{k} + σ̄/n − 1)nμ_{k}
< 0,

since σ_{k} ∈ [σ̄, 1/2]. Hence the minimizer of μ_{k}(α) subject to the conditions (7a) and (7b) lies beyond the interval (26), and we have the result.

We can now show that μ decreases by a factor strictly less than one on each safe iteration.

Theorem 3.6 If a safe step is taken at iteration k, then

μ_{k+1} ≤ μ_{k} [1 − \frac{σ̄(1 − γ_{max})}{8ω}] ,  (30)

where ω is as defined in Lemma 3.4.

Proof. By Theorem 3.5, μ_{k}(α) is decreasing for

α ∈ [0, \frac{σ̄(1 − γ_{max})}{2ω}] ,

and α_{k} lies beyond this interval. Because of (28), we have

μ_{k}(α) = μ_{k}[1 − α(1 − σ_{k})] + α^{2}u^{T} v/n
≤ μ_{k} [1 − \frac{σ̄(1 − γ_{max})}{2ω}(1 − σ_{k})] + \frac{σ̄^{2}(1 − γ_{max})^{2}}{4nω^{2}} ωμ_{k}.

Since 1 − σ_{k} ≥ 1/2, σ̄ ∈ (0, 1/2), 1 − γ_{max} < 1, and n ≥ 1, we have

μ_{k}(α) ≤ μ_{k} [1 − \frac{σ̄(1 − γ_{max})}{4ω}] + \frac{σ̄(1 − γ_{max})}{8ω}μ_{k} = μ_{k} [1 − \frac{σ̄(1 − γ_{max})}{8ω}] ,

giving the desired result.

If we take fast iterations into account, we find that

μ_{k+1} ≤ μ_{k} max (1 − \frac{σ̄(1 − γ_{max})}{8ω}, ρ) ,  (31)

so we have geometric convergence of {μ_{k}} to zero from any starting point (x^{0}, y^{0}) > 0. For the special choice of starting point (14), (18), we have the following polynomial complexity result.

Corollary 3.7 Let ϵ > 0 be given. Suppose that the starting point is defined by (14), (18) where μ_{0} = ξ_{x}ξ_{y} ≤ 1/ϵ^{τ} for some constant τ ≥ 0 independent of n. Then there is an integer K_{ϵ} with

K_{ϵ} = O(n^{2} log(1/ϵ))

such that μ_{k} ≤ ϵ for all k ≥ K_{ϵ}.

Proof. From (30) and the fact that ω = O(n^{2}) (Lemma 3.4), we find that there is a constant δ independent of n such that

μ_{k+1} ≤ (1 − δ/n^{2})μ_{k}  (32)

when a safe step is taken on iteration k. By adjusting δ if necessary, the inequality (32) also holds for fast steps. The result follows from this inequality by a standard argument (see, for example, Zhang [14, Theorem 7.2]).

We note in passing that the same global convergence and complexity results can be obtained even if the exact minimizing α̃ is not found in (6). Instead, we can merely require α̃ to satisfy conditions like those often seen in line search methods for unconstrained minimization. One such pair of conditions is

μ_{k}(α̃) ≤ μ_{k}(0) + ϵα̃μ′_{k}(0),  (33a)
α̃ ≥ τα̃_{max},  (33b)

where ϵ ∈ (0, \frac{1}{2}), τ ∈ (0, 1), and α̃_{max} is the largest value in [0, α̂] for which (33a) holds (where α̂ is defined in (7)). It is easy to show, using the techniques in this section, that the exact minimum from (6) satisfies (33); that Theorem 3.5 is still true if we scale the lower bound by τ; and that Corollary 3.7 continues to hold.

If the relaxed conditions (33) are adopted, we can use techniques more like those in practical codes to choose α̃. In particular, we can use an Armijo-like procedure in which we calculate the step length α_{top} along (u^{k}, v^{k}) to the boundary of the positive orthant, then back off from this length until (33) and (7) are all satisfied. For appropriate choices of the parameters, the popular choice α̃ = .9995 min(1, α_{top}) is usually accepted.

4 Bounds for the Fast-Step Components

In this section we show that when σ_{k} = 0, the step norms ∥u^{k}∥ and ∥v^{k}∥ are both O(μ_{k}) for all sufficiently large k. This result is essential to the local convergence analysis of the next section. For notational convenience, we use u_{B}^{k} to denote the vector whose components are u_{i}^{k} for i ∈ B, u_{N}^{k} as the subvector made up of u_{i}^{k}, i ∈ N, and so on.

We make the following assumption throughout the remainder of this section.

Assumption 2 S^{c} ≠ ∅.

Given any strictly complementary solution (x^{∗}, y^{∗}), we can define the partition

{1, 2, ⋯, n} = N ∪ B,

where

B = {i | x_{i}^{∗} > 0},  N = {i | y_{i}^{∗} > 0}.  (34)

(It is well known that B and N are independent of the particular choice of (x^{∗}, y^{∗}).) We start by showing that x_{N}^{k} and y_{B}^{k} can be bounded in terms of μ_{k}.

Lemma 4.1 Let K_{1/2} be the smallest integer such that ν_{k} ≤ 1/2 for all k ≥ K_{1/2}. Then there is a constant C_{4} > 0 such that

0 < x_{i}^{k} ≤ C_{4}μ_{k}, ∀i ∈ N;  0 < y_{i}^{k} ≤ C_{4}μ_{k}, ∀i ∈ B.  (35)

Proof. Note that K_{1/2} is well defined, by the results of Section 3 and the fact that {ν_{k}} is a decreasing sequence. Let (x^{∗}, y^{∗}) be a strictly complementary solution. By rearranging (16) and noting that

(x^{0})^{T} y + x^{T} y^{0} > 0, (x^{∗})^{T} y^{∗} = 0,

we have

x^{T} y^{∗} + y^{T} x^{∗} ≤ \frac{ν_{k}^{2}}{1 − ν_{k}} nμ_{0} + \frac{nμ_{k}}{1 − ν_{k}} + ν_{k} ((x^{0})^{T} y^{∗} + (x^{∗})^{T} y^{0})
≤ 2ν_{k}^{2}nμ_{0} + 2nμ_{k} + ν_{k} ((x^{0})^{T} y^{∗} + (x^{∗})^{T} y^{0}) .

By (9), we can therefore define a constant C̄_{4} > 0 such that

x^{T} y^{∗} + (x^{∗})^{T} y ≤ C̄_{4}μ_{k}.

Since (x, y) > 0 and (x^{∗}, y^{∗}) ≥ 0 we have

i ∈ N ⇒ x_{i}y_{i}^{∗} ≤ C̄_{4}μ_{k} ⇒ x_{i} ≤ \frac{C̄_{4}μ_{k}}{y_{i}^{∗}}
i ∈ B ⇒ x_{i}^{∗}y_{i} ≤ C̄_{4}μ_{k} ⇒ y_{i} ≤ \frac{C̄_{4}μ_{k}}{x_{i}^{∗}} .

The result (35) follows when we define

C_{4} = C̄_{4} max (max_{i∈B} \frac{1}{x_{i}^{∗}} , max_{i∈N} \frac{1}{y_{i}^{∗}} ) .

The next result gives the required bounds on half the components of the vector pair (u^{k}, v^{k}).

Lemma 4.2 Let K_{1/2} be as defined in Lemma 4.1. Then there is a constant C_{5} > 0 such that for all k ≥ K_{1/2} we have

∥u_{N}^{k}∥ ≤ C_{5}μ_{k},  ∥v_{B}^{k}∥ ≤ C_{5}μ_{k}.  (36)

Proof. For i ∈ N we have from (19) that

\frac{y_{i}}{x_{i}} (u_{i})^{2} = |D_{ii}u_{i}|^{2} ≤ ∥Du∥^{2} ≤ ωμ_{k}.

Hence from (35) and (8), we have

(u_{i})^{2} ≤ \frac{ωμ_{k}x_{i}}{y_{i}} = \frac{ωμ_{k}x_{i}^{2}}{x_{i}y_{i}} ≤ \frac{ωμ_{k}C_{4}^{2}μ_{k}^{2}}{γ_{min}μ_{k}} = \frac{ωC_{4}^{2}}{γ_{min}} μ_{k}^{2}.

Therefore

∥u_{N}∥ ≤ (\frac{nω}{γ_{min}} )^{1/2} C_{4}μ_{k},

giving the first inequality in (36). The proof of the second inequality is similar.

To obtain bounds on the remaining components of (u^{k}, v^{k}), we need the following two technical lemmas.

Lemma 4.3 (Monteiro and Wright [6, Lemma 2.2].) Let f ∈ IR^{q} and H ∈ IR^{p×q} be given. Then there exists a nonnegative constant L = L(f, H) with the property that for any diagonal matrix S > 0 and any vector h ∈ Range(H), the (unique) optimal solution w̄ = w̄(S, h) of

min_{w} f^{T} w + \frac{1}{2} ∥Sw∥^{2} ,  subject to Hw = h,  (37)

satisfies

\left∥{ \atop w̄}\right∥_{∞} ≤ L \{|f^{T} w̄| + ∥h∥_{∞}\} .

The second technical lemma identifies the components u_{B}^{k} and v_{N}^{k} as the solution of a quadratic program. It is an extension of a result of Ye and Anstreicher [12, Lemma 3.5] and is proved in Wright [10]. We use D_{B}^{k} to denote the diagonal submatrix composed of the elements D_{ii}^{k} for i ∈ B, and so on.

Lemma 4.4 (Wright [10, Lemma 5.2]) The vector pair (u_{B}^{k}, v_{N}^{k}) solves the convex quadratic program

min_{(w,z)} \frac{1}{2}∥D_{B}^{k}w∥^{2} − σ_{k}μ_{k}e_{B}^{T}(X_{B}^{k})^{−1}w + \frac{1}{2}∥(D_{N}^{k})^{−1}z∥^{2} − σ_{k}μ_{k}e_{N}^{T}(Y_{N}^{k})^{−1}z,  (38)

subject to

M_{BB}w = r_{B}^{k} − M_{BN}u_{N}^{k} + v_{B}^{k},  (39a)
M_{NB}w − z = r_{N}^{k} − M_{NN}u_{N}^{k}.  (39b)

The main result of this section is as follows.

Theorem 4.5 Let K_{1/2} be as defined in Lemma 4.1. Then if σ_{k} = 0, there is a constant C_{6} > 0 such that

∥u^{k}∥ ≤ C_{6}μ_{k}, ∥v^{k}∥ ≤ C_{6}μ_{k}.

Proof. It follows from Lemmas 4.3 and 4.4 and the inequality (9) that there is a constant L > 0 such that

∥(u_{B}^{k}, v_{N}^{k})∥ ≤ L(∥r^{k}∥ + ∥u_{N}^{k}∥ + ∥v_{B}^{k}∥) ≤ L(ν_{k}∥r^{0}∥ + 2C_{5}μ_{k}) ≤ L (∥r^{0}∥/(β̂μ_{0}) + 2C_{5}) μ_{k}.

By combining this inequality with (36), we obtain the result.

Note that the bound on the solution of (37) is independent of the positive diagonal matrix S. This feature is not present in the analysis of [10, Section 5] and [9, Section 5], where estimates of the elements of D_{B}^{k} and D_{N}^{k} are used in the bound for (u_{B}^{k}, v_{N}^{k}). In [10], these estimates follow from boundedness of the iterates (x^{k}, y^{k}) which, in turn, follow from existence of a strictly feasible point. Since the estimates are no longer required, the strict feasibility assumption is not required either. In [9], the estimates are obtained by synchronizing reductions in μ_{k} and ∥r_{k}∥, which complicates the algorithm considerably. Awkward technical devices such as auxiliary sequences are used in both paper, but are not needed here.

5 Local Superlinear Convergence

In this section, we show that for all k sufficiently large, step k is a fast step and that consequently the sequence {μ_{k}} converges subquadratically to zero. The treatment in this section follows that of Wright [10, Section 6] and [9, Section 6].

Throughout the analysis, we will make use of the constant C_{7} defined by

C_{7} =^{△} max (1, 2C_{6}^{2}) .  (40)

We start with a simple technical result.

Lemma 5.1 For all k ≥ 0, we have

\frac{μ_{k+1}}{γ̄^{t_{k+1}}} ≤ η \frac{μ_{k}}{γ̄^{t_{k}}} ,

where

η = max (1 − \frac{σ̄(1 − γ_{max})}{8ω}, \frac{ρ}{γ̄} ) < 1,

and so the sequence

\{\frac{μ_{k}}{γ̄^{t_{k}}} \}

converges monotonically and geometrically to zero.

Proof. When a safe step is taken, we have from (30) and t_{k} = t_{k+1} that

\frac{μ_{k+1}}{γ̄^{t_{k+1}}} ≤ [1 − \frac{σ̄(1 − γ_{max})}{8ω}] \frac{μ_{k}}{γ̄^{t_{k}}} .

For a fast step, we have μ_{k+1} ≤ ρμ_{k} and t_{k+1} = t_{k} + 1 and so

\frac{μ_{k+1}}{γ̄^{t_{k+1}}} ≤ (\frac{ρ}{γ̄} ) \frac{μ_{k}}{γ̄^{t_{k}}} .

The result follows immediately from these two bounds.

We now show that fast steps are taken for all k sufficiently large. Theorem 5.2 Define K to be the smallest index such that

\frac{μ_{k}}{γ̄^{t_{k}}(γ_{max} − γ_{min})} ≤ \frac{ρ}{2C_{7}}  (41)

for all k ≥ K. Then a fast step is taken on iteration k, with step length α_{k} satisfying

1 ≥ α_{k} ≥ 1 − C_{7} \frac{μ_{k}}{γ̄^{t_{k}}(γ_{max} − γ_{min})}.  (42)

Moreover, we have

μ_{k+1} ≤ \frac{2C_{7}}{γ̄^{t_{k}}(γ_{max} − γ_{min})}μ_{k}^{2} ≤ ρμ_{k}.  (43)

Proof. The proof is structured like that of Theorem 3.5, in that we show that the conditions (7a) and (7b) hold for all α satisfying

α ∈ [0, 1 − C_{7} \frac{μ_{k}}{γ̄^{t_{k}}(γ_{max} − γ_{min})}]  (44)

and then show that μ_{k}(α) is decreasing on this interval.

We start with condition (7a). Note from Theorem 4.5, (40), and α ∈ (0, 1], we have

(x + αu)^{T} (y + αv)/n = (1 − α)μ_{k} + α^{2}u^{T} v/n ≥ (1 − α)μ_{k} − \frac{C_{7}}{2} μ_{k}^{2}.  (45)

Now for α in the interval (44), we have

α ≤ 1 − C_{7} \frac{μ_{k}}{γ̄^{t_{k}}(γ_{max} − γ_{min})} ⇒ \frac{C_{7}μ_{k}}{2} ≤ (1 − α)γ̄^{t_{k}} \frac{γ_{max} − γ_{min}}{2} ≤ (1 − α)γ̄^{t_{k}}.

Therefore

(1−α)μ_{k} − \frac{C_{7}μ_{k}}{2} μ_{k} ≥ (1−α)μ_{k} −γ̄^{t_{k}}(1−α)μ_{k} = (1−γ̄^{t_{k}})(1−α)μ_{k} = (1−β_{k})(1−α)μ_{k}, (46)

since β_{k} = γ̄^{t_{k}} for a fast step. Relation (7a) follows from (45) and (46).

For (7b), we again use Theorem 4.5, (5), (40), and (8) to derive

(x_{i}+αu_{i})(y_{i}+αv_{i}) = x_{i}y_{i}(1−α)+α^{2}u_{i}v_{i} ≥ γ_{k}(1−α)μ_{k}−∥u∥∥v∥ ≥ γ_{k}(1−α)μ_{k}− \frac{C_{7}}{2} μ_{k}^{2}. (47)

Meanwhile, assuming that a fast step is computed, we have

γ_{k+1}(x + αu)^{T} (y + αv)/n = γ_{k+1}(1 − α)μ_{k} + γ_{k+1}α^{2}u^{T} v/n ≤ γ_{k+1}(1 − α)μ_{k} + \frac{C_{7}}{2} μ_{k}^{2}. (48)

Combining (47) and (48), we find that (7b) is satisfied if

(γ_{k} − γ_{k+1})(1 − α) ≥ C_{7}μ_{k}.  (49)

Now, since γ̄ ∈ (0, 1/2], we have

γ_{k}−γ_{k+1} = γ̄^{t_{k}−1}(γ_{max}−γ_{min})−γ̄^{t_{k}}(γ_{max}−γ_{min}) = γ̄^{t_{k}−1}(1−γ̄)(γ_{max}−γ_{min}) ≥ γ̄^{t_{k}}(γ_{max}−γ_{min}).

Therefore (49) is satisfied if

γ̄^{t_{k}}(γ_{max} − γ_{min})(1 − α) ≥ C_{7}μ_{k}.  (50)

But (50) is clearly satisfied for all α in the interval (44), so we deduce that (7b) holds. Finally, we examine μ_{k}(α) from (6). For α ∈ [0, 1] we have

μ′_{k}(α) = −μ_{k} + 2αu^{T} v/n ≤ −μ_{k} + C_{7}μ_{k}^{2} = −(1 − C_{7}μ_{k})μ_{k}.

Now from (41), we clearly have C_{7}μ_{k} < 1, and therefore μ′_{k}(α) < 0. Hence the complementarity gap is certainly decreasing on the interval (44). We deduce that the step length α_{k} lies above the upper bound of the interval (44), so the proof of (42) is complete.

For (43), we use (42) to obtain

μ_{k+1} = (x + α_{k}u)^{T} (y + α_{k}v)/n
= (1 − α_{k})μ_{k} + α_{k}^{2}u^{T} v/n
≤ \frac{C_{7}}{γ̄^{t_{k}}(γ_{max} − γ_{min})}μ_{k}^{2} + \frac{C_{7}}{2} μ_{k}^{2}
≤ \frac{2C_{7}}{γ̄^{t_{k}}(γ_{max} − γ_{min})}μ_{k}^{2},  (51)

giving the first inequality in (43). The second inequality is an immediate consequence of (41).

We can now state our asymptotic rate-of-convergence result.

Theorem 5.3 The sequence {μ_{k}} converges superlinearly to zero with Q-order 2.

Proof. See [10, Theorem 6.3 (ii)].

So far, we have used the term “convergence” to denote convergence of μ_{k} and ∥r^{k}∥ to zero. Fast convergence of the actual iterates (x^{k}, y^{k}) to a solution was not proved in Wright [10], [9], but it is easy to show. A key result is a bound on the distance of (x^{k}, y^{k}) to the solution set S.

Lemma 5.4 Suppose that Assumption 2 holds. Then there is a constant C_{3} such that

dist ((x^{k}, y^{k}), S) =^{△} min_{(x,y)∈S} ∥(x^{k}, y^{k}) − (x, y)∥_{∞} ≤ C_{3}μ_{k}.

Proof. Note first that for any (x^{∗}, y^{∗}) ∈ S, we have

y − y^{∗} = (Mx + q + r) − (Mx^{∗} + q) = M(x − x^{∗}) + r,

and so from Lemma 3.1

∥(x, y) − (x^{∗}, y^{∗})∥_{∞} ≤ (1 + ∥M∥_{∞})∥x − x^{∗}∥_{∞} + ν_{k}∥r^{0}∥_{∞}
≤ (1 + ∥M∥_{∞})∥x − x^{∗}∥_{∞} + μ_{k}∥r^{0}∥_{∞}/(β̂μ_{0}).  (52)

Now from Mangasarian [3, Theorem 2.6], we have that there is a constant C̄_{3} such that

min_{(x^{∗},Mx^{∗}+q)∈S} ∥x − x^{∗}∥_{∞} ≤ C̄_{3} ∥(−Mx − q, −x, x^{T} (Mx + q))_{+}∥_{2}
≤ C̄_{3} ∥(r − y, x^{T} (y − r))_{+}∥_{2}
≤ C̄_{3}(∥r∥_{2} + x^{T} y + ∥x∥_{1}∥r∥_{∞})
≤ C̄_{3}(ν_{k}∥r^{0}∥_{2} + nμ_{k} + ν_{k}∥x∥_{1}∥r^{0}∥_{∞}).

If we substitute from (13) and (9), we obtain from the last inequality that there is a constant Ĉ_{3} such that

min_{(x^{∗},Mx^{∗}+q)∈S} ∥x − x^{∗}∥_{∞} ≤ Ĉ_{3}μ_{k}.  (53)

The result follows by combining (52) with (53).

Theorem 5.5 If Assumption 2 holds, the iterates (x^{k}, y^{k}) converge R-subquadratically to a point in S.

Proof. We show first that the sequence {(x^{k}, y^{k})} is Cauchy. Since fast steps are taken for all k ≥ K (Theorem 5.2), then for any k_{1}, k_{2} with K ≤ k_{1} < k_{2}, we have from Theorem 4.5 that

∥(x^{k_{2}}, y^{k_{2}}) − (x^{k_{1}}, y^{k_{1}})∥ ≤ ∑_{k=k_{1}}^{k_{2}−1} ∥(u^{k}, v^{k})∥ ≤ \sqrt{2}C_{6} ∑_{k=k_{1}}^{k_{2}−1} μ_{k}.

From (43), we find that {μ_{k}} is bounded by a geometric sequence and therefore

∑_{k=k_{1}}^{k_{2}−1} μ_{k} ≤ ∑_{k=k_{1}}^{k_{2}−1} ρ^{k−k_{1}}μ_{k_{1}} ≤ \frac{μ_{k_{1}}}{1 − ρ}.

Hence,

∥(x^{k_{2}}, y^{k_{2}}) − (x^{k_{1}}, y^{k_{1}})∥ ≤ \frac{\sqrt{2}C_{6}}{1 − ρ} μ_{k_{1}},  (54)

which approaches zero as k_{1}, k_{2} → ∞. Hence, the sequence is Cauchy and therefore convergent to a limit point (x^{∗}, y^{∗}) (say). Because of Lemma 5.4 and the fact that S is closed, we have (x^{∗}, y^{∗}) ∈ S.

By using the same arguments that led to (54), we have for k ≥ K that

∥(x^{k}, y^{k}) − (x^{∗}, y^{∗})∥ ≤ \frac{\sqrt{2}C_{6}}{1 − ρ} μ_{k}.

The R-subquadratic convergence of (x^{k}, y^{k}) to (x^{∗}, y^{∗}) follows from this bound and Theorem 5.3.

6 Extension to MLCP, LP, and QP

We conclude by showing that the algorithm of this paper can be extended to the mixed LCP (2) and, hence, to all the usual formulations of linear and convex quadratic programming problems. The crucial result here is due to Güler [2], who showed that any generalized linear complementarity problem involving a maximal monotone operator can be reformulated as a standard LCP (1). The analysis in this section shows that the operator represented by (2a) is in fact maximal monotone and, hence, satisfies the assumptions of Theorem 3.2 in [2]. The following extension of Assumption 1 is essential to our analysis.

Assumption 3 The MLCP (2) is feasible; that is, there is a vector triple (x, y, z) with (x, y) ≥ 0 satisfying (2a).

We define T to be a multivalued mapping from IR^{n} to subsets of IR^{n} such that y ∈ T (x) whenever there exists a z ∈ IR^{m} such that (x, y, z) satisfies (2a). The graph of T is given by

G(T) = {(x, y) | y = M_{11}x + M_{12}z + q_{1}, 0 = M_{21}x + M_{22}z + q_{2}, some z ∈ IR^{m}}. (55)

Because of Assumption 3, we have G(T) ≠ ∅. It is easy to check that T is monotone, that is, for all (x, y) ∈ G(T) and (x̄, ȳ) ∈ G(T), we have (x − x̄)^{T} (y − ȳ) ≥ 0. In Theorem 6.3 below, we show that T is in fact maximal monotone, that is, there is no other monotone operator T̄ : IR^{n} → IR^{n} such that G(T) is strictly contained in G(T̄). The following two technical lemmas lay the foundation for this theorem.

Lemma 6.1 Suppose M is positive semidefinite. Then Mu = 0 ⇔ M^{T} u = 0.

Proof.

Mu = 0 ⇒ u^{T} Mu = 0 ⇒ u^{T} (M + M^{T})u = 0,

and so u is a minimizer of the convex quadratic function f(u) = u^{T} (M + M^{T})u. Therefore

∇f(u) = 0 ⇒ (M + M^{T})u = 0 ⇒ M^{T} u = 0,

and the forward implication is proved. The converse is similar.

Lemma 6.2 Let the p.s.d. matrix M be partitioned as

M = [ \begin{matrix}M_{11} & M_{12} \\ M_{21} & M_{22} \end{matrix} ]

where M_{11} and M_{22} are square and

M_{·2} =^{△} [ \begin{matrix}M_{12} \\ M_{22} \end{matrix} ]

has full column rank. Then if D is a diagonal matrix with strictly positive diagonal entries, then

[ \begin{matrix}M_{11} + D & M_{12} \\ M_{21} & M_{22} \end{matrix} ]

is nonsingular.

Proof. Let (x, z) be a vector pair such that

[ \begin{matrix}M_{11} + D & M_{12} \\ M_{21} & M_{22} \end{matrix} ] [ \begin{matrix}x \\ z \end{matrix} ] = [ \begin{matrix}0 \\ 0 \end{matrix} ] .  (56)

Then, using the p.s.d. property of M, we have

[ x^{T} z^{T} ] [ \begin{matrix}M_{11} + D & M_{12} \\ M_{21} & M_{22} \end{matrix} ] [ \begin{matrix}x \\ z \end{matrix} ] = 0

⇒ [ x^{T} z^{T} ] M [ \begin{matrix}x \\ z \end{matrix} ] + x^{T} Dx = 0 ⇒ x^{T} Dx ≤ 0 ⇒ x = 0.

From (56) we have M_{·2}z = 0, and thus z = 0 by the full rank assumption. Hence (56) is satisfied only by (x, z) = (0, 0), and the result is proved.

Our main result is the following.

Theorem 6.3 Suppose that Assumption 3 holds. Then the mapping T whose graph is given by (55) is maximal monotone.

Proof. We prove the result by appealing to a theorem of Minty [4], which states that T is maximal monotone if and only if R(I + T) = IR^{n}. In the remainder of the proof, we show that for any w ∈ IR^{n}, there is an (x, y) ∈ G(T) such that x+y = w, and hence w ∈ R(I +T). In other words, for any w the following linear system must have a solution triple (x, y, z):

[ \begin{matrix}M_{11} & M_{12} \\ M_{21} & M_{22} \end{matrix} ] [ \begin{matrix}x \\ z \end{matrix} ] − [ \begin{matrix}y \\ 0 \end{matrix} ] = [ \begin{matrix}−q_{1} \\ −q_{2} \end{matrix} ] ,  x + y = w.  (57)

We show first that solutions to (57) can be obtained from solutions to the following system:

[ \begin{matrix}M_{11} & M̄_{12} & −I \\ M̄_{21} & M̄_{22} & 0 \\ I & 0 & I \end{matrix} ] [ x ; z̄ ; y ] = [ \begin{matrix}−q_{1} \\ −q̄_{2} \\ w \end{matrix} ] ,  (58)

where z̄ ∈ IRm̄, and q̄_{2} ∈ IRm̄ and the matrices M̄_{12}, M̄_{21}, M̄_{22} are quantities to be defined, with the crucial property that the submatrix

M̄_{·2} = [ M̄_{12} ; M̄_{22} ] ∈ IR^{(n+m̄)×m̄}  (59)

has full column rank.

Let v ∈ IR^{m} be any vector such that M_{·2}v = 0, where M_{·2} is defined as in Lemma 6.2. Then

M [ \begin{matrix}0 \\ v \end{matrix} ] = 0

and hence, by Lemma 6.1,

M^{T} [ \begin{matrix}0 \\ v \end{matrix} ] = 0 ⇒ [ \begin{matrix}M_{21}^{T} \\ M_{22}^{T} \end{matrix} ] v = 0.  (60)

Let V_{2} ∈ IR^{m×d} (0 ≤ d ≤ m) be the matrix whose columns form a basis for the subspace with the property M_{·2}v = 0, and define V_{1} ∈ IR^{m×(m−d)} so that the columns of V_{1} are a basis for R(V_{2})^{⊥}. Then, defining the nonsingular matrix V = [V_{1} | V_{2}], we have

[ \begin{matrix}M_{11} & M_{12} \\ M_{21} & M_{22} \end{matrix} ] [ \begin{matrix}x \\ z \end{matrix} ] = [ \begin{matrix}y − q_{1} \\ −q_{2} \end{matrix} ]

⇔ [ \begin{matrix}I & 0 \\ 0 & V^{T} \end{matrix} ] [ \begin{matrix}M_{11} & M_{12} \\ M_{21} & M_{22} \end{matrix} ] [ \begin{matrix}I & 0 \\ 0 & V \end{matrix} ] [ \begin{matrix}I & 0 \\ 0 & V^{−1} \end{matrix} ] [ \begin{matrix}x \\ z \end{matrix} ] = [ \begin{matrix}y − q_{1} \\ −V^{T} q_{2} \end{matrix} ]

⇔ [ \begin{matrix}M_{11} & M̄_{12} & 0 \\ M̄_{21} & M̄_{22} & 0 \\ 0 & 0 & 0 \end{matrix} ] [ \begin{matrix}x \\ (V^{−1}z)_{1} \\ (V^{−1}z)_{2} \end{matrix} ] = [ \begin{matrix}y − q_{1} \\ −V_{1}^{T} q_{2} \\ −V_{2}^{T} q_{2} \end{matrix} ] ,  (61)

where

M̄_{12} = M_{12}V_{1},  M̄_{22} = V_{1}^{T} M_{22}V_{1},  M̄_{21} = V_{1}^{T} M_{21}.

To verify that M̄_{·2} has full column rank, let w be a vector for which M̄_{·2}w = 0. Then

[ \begin{matrix}M_{12} \\ V_{1}^{T} M_{22} \end{matrix} ] V_{1}w = 0 ⇒ [ \begin{matrix}M_{12} \\ V^{T} M_{22} \end{matrix} ] V_{1}w = 0 ⇒ M_{·2}V_{1}w = 0 ⇒ V_{1}w ∈ R(V_{2}),

and hence w = 0 by definition of V_{1} and V_{2}, proving the assertion.

For (61) to be consistent, it is necessary that V_{2}^{T} q_{2} = 0. This follows, however, from Assumption 3 and (60), since

q_{2} ∈ R([M_{21} | M_{22}]) ⇒ V_{2}^{T} q_{2} = 0.

Defining m̄ = m − d and

z̄ = (V^{−1}z)_{1} (first m̄ components of V^{−1}z),
ẑ = (V^{−1}z)_{2} (last m − m̄ components of V^{−1}z),
q̄_{2} = V_{1}^{T} q_{2},

we find that any solution (x, y, z̄) of (58) can be transformed into a solution (x, y, z) of (57) by setting

z = V [ z̄ ; ẑ ] ,

where ẑ is chosen arbitrarily.

It remains to show that (58) does in fact have a solution. By adding the third block row in (58) to the first block row, we obtain the equivalent system

[ \begin{matrix}M_{11} + I & M̄_{12} & 0 \\ M̄_{21} & M̄_{22} & 0 \\ I & 0 & I \end{matrix} ] [ x ; z̄ ; y ] = [ \begin{matrix}−q_{1} + w \\ −q̄_{2} \\ w \end{matrix} ] .  (62)

Since, by our choice of transformations V_{1} and V_{2}, we clearly have that

M̄ = [ \begin{matrix}M_{11} & M̄_{12} \\ M̄_{21} & M̄_{22} \end{matrix} ]

is p.s.d., while M̄_{·2} has full rank, it follows from Lemma 6.2 that the upper left 2 × 2 block in the coefficient matrix of (62) is nonsingular and, hence, that (62) has a unique solution triple (x, y, z̄). Hence we have identified (x, y) ∈ G(T) with x + y = w, and our proof is complete.

Because of (61), the graph of T can be restated as

G(T) = {(x, y) | y = M_{11}x + M̄_{12}z̄ + q_{1}, 0 = M̄_{21}x + M̄_{22}z̄ + q̄_{2}, some z̄ ∈ IRm̄}, with M̄_{·2} full rank. If we define a matrix W ∈ IR^{(n+m̄)×n} such that the columns of W form a basis for R(M̄_{·2})^{⊥}, we can eliminate z̄ from the definition of G(T) altogether and write

G(T) = {(x, y) | F x − Gy = a},  (63)

where

F = W^{T} [ M_{11} ; M̄_{21} ] ,  G = W^{T} [ \begin{matrix}I \\ 0 \end{matrix} ] ,  a = W^{T} [ \begin{matrix}−q_{1} \\ −q̄_{2} \end{matrix} ] .

The form (63) is the canonical form used by Güler. It is not difficult to verify that a condition used by Güler to prove maximal monotonicity of T —namely, nonsingularity of F + G—is satisfied by (63). As shown in [2, Theorem 3.2] conversion of (63) to the form (1) can now be achieved by premultiplying F x − Gy = a by a nonsingular operator and possibly swapping some components of x and y. As Güler notes, this reformulation process need not actually be carried out to apply the algorithm of Section 2 to the problem (2). Instead, we can extend our algorithm to (2) and note that the sequence of (x^{k}, y^{k}) iterates generated by the extended algorithm is the same as the sequence that would be generated by the basic algorithm applied to the LCP reformulation, subject possibly to the swapping of components between x^{k} and y^{k} just mentioned.

The extension of our algorithm to (2) is fairly obvious, so we omit the details, noting simply that the linear system to be solved at each iteration is

[ \begin{matrix}M_{11} & M_{12} & −I \\ M_{21} & M_{22} & 0 \\ Y^{k} & 0 & X^{k} \end{matrix} ] [ \begin{matrix}u^{k} \\ w^{k} \\ v^{k} \end{matrix} ] = [ \begin{matrix}y^{k} − M_{11}x^{k} − M_{12}z^{k} − q_{1} \\ −M_{21}x^{k} − M_{22}z^{k} − q_{2} \\ −X^{k}Y^{k}e + σ̃μ_{k}e \end{matrix} ] ,  (64)

which is a generalization of (5). Existence of a solution to (64) follows from the same reasoning as in the proof of Theorem 6.3. We need to note that when (2) is feasible, we have

−M_{21}x^{k} − M_{22}z^{k} − q_{2} ∈ R([M_{21} | M_{22}])

and, hence,

V_{2}^{T} [−M_{21}x^{k} − M_{22}z^{k} − q_{2}] = 0.

Also, we need to apply Lemma (6.2) with D = (X^{k})^{−1}Y^{k}. Although the step components u^{k} and v^{k} are uniquely determined by (64), the w^{k} components are not, unless the submatrix M_{·2} has full rank.

Finally, we note that our analysis depends crucially on the existence of feasible points for (1) and (2). When (2) arises from an LP or QP, a feasible primal-dual point must exist. This requirement is a little troubling, since in practice many LPs are either primal or dual infeasible. Ye, Todd, and Mizuno [13] and Xu, Hung, and Ye [11] have alleviated this difficulty in the case of linear programming by describing augmentation/reformulation of an LP in standard form, which has the property that the resulting mixed LCP possesses both a feasible point and a strictly complementary solution. Since all the assumptions of this paper are satisfied by these reformulations, our algorithm can be applied with confidence.

References

[1] R. W. Cottle, J.-S. Pang, and R. E. Stone, The Linear Complementarity Problem (Academic Press, 1992).

[2] O. Güler, Generalized linear complementarity problems, Manuscript (1993).

[3] O. L. Mangasarian, Error bounds for nondegenerate monotone linear complementarity problems, Math. Progr. 48(1990) 437–445.

[4] G. J. Minty, Monotone (nonlinear) operators in Hilbert space, Duke J. Math. 29(1992) 341–346.

[5] S. Mizuno, Polynomiality of Kojima-Meggido-Mizuno infeasible-interior-point algorithm for linear programming, Tech. Report 1006, School of Operations Research and Industrial Engineering, Cornell University, Ithaca, N.Y. (1993).

[6] R. D. C. Monteiro and S. J. Wright, Superlinear primal-dual affine scaling algorithms for LCP, Preprint MCS–P360–0693, Mathematics and Computer Science Division, Argonne National Laboratory, Argonne, Ill. (1993). To appear in Math. Progr., Series A.

[7] R. D. C. Monteiro and S. J. Wright, Local convergence of interior-point algorithms for degenerate monotone LCP, Comput. Optim. Appls 3(1994) 131–155.

[8] F. A. Potra, An infeasible interior-point predictor-corrector algorithm for linear programming, Tech. Report 26, Department of Mathematics, University of Iowa, Iowa City, Iowa (1992). To appear in SIAM J. Optim.

[9] S. J. Wright, A path-following infeasible-interior-point algorithm for linear complementarity problems, Optim. Methods Softw. 2(1993) 79–106.

[10] S. J. Wright, An infeasible-interior-point algorithm for linear complementarity problems, Math. Progr. 67(1994) 29–52.

[11] X. Xu, P. Hung, and Y. Ye, A simplified homogeneous and self-dual linear programming algorithm and its implementation, Manuscript (1993).

[12] Y. Ye and K. Anstreicher, On quadratic and O(\sqrt{n}L) convergence of a predictorcorrector algorithm for LCP, Math. Progr. Series A 62(1993) 537–551.

[13] Y. Ye, M. J. Todd, and S. Mizuno, An O(\sqrt{n}L)-iteration homogeneous and self-dual linear programming algorithm, Working Paper 92-11, College of Business Administration, The University of Iowa (1992). To appear in Math. Oper. Res.

[14] Y. Zhang, On the convergence of a class of infeasible-interior-point methods for the horizontal linear complementarity problem, SIAM J. Optim. 4(1994) 208–227.