Local Convergence of Interior-Point Algorithms for Degenerate Monotone LCP

Renato D.C. Monteiro^{∗} and Stephen Wright^{†}

April 6, 1993

PREPRINT MCS-P357-0493, MATHEMATICS AND COMPUTER SCIENCE DIVISION, ARGONNE NATIONAL LABORATORY

Abstract Most asymptotic convergence analysis of interior-point algorithms for monotone linear complementarity problems assumes that the problem is nondegenerate, that is, the solution set contains a strictly complementary solution. We investigate the behavior of these algorithms when this assumption is removed.

1 Introduction

In the monotone linear complementarity problem (LCP), we seek a vector pair (x, y) ∈ IR^{n} × IR^{n} that satisfies the conditions

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

where q ∈ IR^{n}, and M ∈ IR^{n×n} is positive semidefinite. We use S to denote the solution set of (1).

An assumption that is frequently made in order to prove superlinear convergence of interior-point algorithms for (1) is the nondegeneracy assumption:

Assumption 1 There is an (x^{∗}, y^{∗}) ∈ S such that x_{i}^{∗} + y_{i}^{∗} > 0 for all i = 1, ⋯, n.

In general, we can define three subsets B, N, and J of the index set {1, ⋯, n} by

B = {i = 1, ⋯, n | x_{i}^{∗} > 0 for at least one (x^{∗}, y^{∗}) ∈ S},
N = {i = 1, ⋯, n | y_{i}^{∗} > 0 for at least one (x^{∗}, y^{∗}) ∈ S},  (2)
J = {i = 1, ⋯, n | x_{i}^{∗} = y_{i}^{∗} = 0 for all (x^{∗}, y^{∗}) ∈ S}.

It is well known that B, N, and J form a partition of {1, ⋯, n}. Another useful result is the following.

∗Systems and Industrial Engineering Department, University of Arizona, Tucson, AZ 85721. The work of this author was based on research supported by the National Science Foundation under grant DDM-9109404 and the Office of Naval Research under grant N00014-93-1-0234.

†Mathematics and Computer Science Division, Argonne National Laboratory, 9700 South Cass Avenue, Argonne, IL 60439. 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.

Lemma 1.1 There is an (x^{∗}, y^{∗}) ∈ S such that x_{i}^{∗} > 0 for all i ∈ B and y_{i}^{∗} > 0 for all i ∈ N.

Proof. Choose |B| + |N| members (x^{i}, y^{i}) of S (where | · | denotes set cardinality) with the property that x_{i}^{i} > 0 for i ∈ B and y_{i}^{i} > 0 for i ∈ N. Define

(x^{∗}, y^{∗}) = \frac{1}{|B| + |N|} ∑_{i∈B∪N}(x^{i}, y^{i}).

Since (x^{i})^{T} y^{j} = (x^{j})^{T} y^{i} = 0 for any two solutions (x^{i}, y^{i}) and (x^{j}, y^{j}) of (1), it is easy to check that y^{∗} = Mx^{∗} + q and (x^{∗})^{T} y^{∗} = 0. Moreover, x_{i}^{∗} > 0 for all i ∈ B and y_{i}^{∗} > 0 for all i ∈ N, giving the result.

Assumption 1 can be restated simply as J = ∅.

An infeasible-interior-point algorithm solves (1) by generating a sequence of strictly positive iterates {(x^{k}, y^{k})}, k = 0, 1, 2, ⋯, while aiming to satisfy the two equality relationships in (1) in the limit as k → ∞. Feasible interior-point algorithms require all iterates to satisfy y^{k} = Mx^{k} + q in addition to strict positivity (x^{k}, y^{k}) > 0.

This paper starts with general results about infeasible-interior-point algorithms. In our presentation, “Q-superlinear convergence” always means Q-superlinear convergence of the complementarity gap (x^{k})^{T} y^{k} to 0. In Section 2, we define the broad class of infeasible algorithms considered in this paper and show that no algorithm of this class can achieve Q-superlinear convergence when J ≠ ∅. Ye and Anstreicher [10] presented a trivial LCP for which J ≠ ∅ and observed that no feasible algorithm whose steps approach the primal-dual affine scaling directions can converge superlinearly for this example. Our result generalizes Ye and Anstreicher’s observation to all instances of problem (1) for which J ≠ ∅ and to a broader class of algorithms that includes many infeasible algorithms.

Section 3 proposes a scheme for estimating the index sets B, N, and J and shows that a finite termination scheme based on these estimates eventually yields an exact solution of (1). The results of this section generalizes the results obtained in Ye [9] for linear programs to the context of degenerate monotone LCPs.

Feasible algorithms in which the step vectors asymptotically converge to the primal-dual affine scaling direction are discussed in Section 4. In this case, we are able to prove stronger results about the asymptotic behavior of the steps (Δx^{k}, Δy^{k}) and to derive formulae for the linear rate of convergence of the complementarity gap in terms of the current iterate and the current stepsize.

In Section 5, we analyze the two-step linear rate of convergence of the predictor-corrector algorithm. We show that this rate is linear with a constant of at most 1 − c/Q^{1/4}, where Q ≤ min{|J|, n−|J|} and c is not too small. This result shows that if |J| contains only a few indices, or it contains all but a few indices, then we can expect a linear rate of convergence that is not too slow.

The following notation is used throughout the paper. Superscripts on matrices and vectors and subscripts on index sets and scalars denote iteration indices (usually k), while subscripts on matrices and vectors define components. The subvector x_{B} denotes [x_{i}]_{i∈B}, while the submatrix M_{BN} is [M_{ij}]_{i∈B,j∈N}. Subvectors and submatrices corresponding to other index sets are defined likewise. We use μ_{k} to denote the normalized complementarity gap μ_{k} = (x^{k})^{T} y^{k}/n. If w ∈ IR^{n} then diag (w) denotes the diagonal matrix having the components of w as diagonal entries. The matrices X^{k} and Y^{k} are defined as X^{k} = diag (x^{k}) and Y^{k} = diag (y^{k}). If u and v are two vectors of the same length and ν ∈ IR, then uv denotes the vector whose i-the component is u_{i}v_{i} and u^{ν} denotes the vector whose i-th component is u_{i}^{ν}. The vector (1, ⋯, 1)^{T}, regardless of its dimension, is denoted by e. For any vector x, the notation x_{+} is used for the vector whose i-th component is max(x_{i}, 0). Unless otherwise specified, ∥ · ∥ denotes ∥ · ∥_{2}. If {δ_{k}} and {ϵ_{k}} are two positive sequences, we say δ_{k} = O(ϵ_{k}) if lim sup_{k} δ_{k}/ϵ_{k} < ∞ and δ_{k} = o(ϵ_{k}) if lim_{k} δ_{k}/ϵ_{k} = 0. The distance function to the solution set S is defined as

dist((x, y), S) = min_{(x̄,ȳ)∈S} ∥(x̄, ȳ) − (x, y)∥.

2 Infeasible Algorithms: Local Convergence

In this section we largely restrict ourselves to discussing interior-point algorithms that fit the following framework, which we refer to as the standard framework.

(a) For all k ≥ 0 we have (x^{k}, y^{k}) > 0 and

∥(x^{k}, y^{k})∥_{∞} ≤ C_{b},  (3)

for some constant C_{b} > 0.

(b) There is a γ̄ ∈ (0, 1) such that

x_{i}^{k}y_{i}^{k} ≥ γ̄μ_{k}  for all i = 1, ⋯, n and k ≥ 0.  (4)

(c) The step has the form (x^{k+1}, y^{k+1}) = (x^{k}, y^{k}) + α_{k}(Δx^{k}, Δy^{k}) with α_{k} ∈ (0, 1]. The search direction (Δx^{k}, Δy^{k}) satisfies the equation

[ \begin{matrix}M & −I \\ Y^{k} & X^{k} \end{matrix} ] [ \begin{matrix}Δx^{k} \\ Δy^{k} \end{matrix} ] = [ \begin{matrix}r^{k} \\ −x^{k}y^{k} + σ_{k}μ_{k}e \end{matrix} ] ,  (5)

where r^{k} = y^{k} − Mx^{k} − q denotes the residual vector and σ_{k} ∈ [0, 1] is the centering parameter.

(d) There are constants \underline{ρ} and ρ̄ such that 0 ≤ \underline{ρ} ≤ ρ̄ and

\underline{ρ}\frac{∥r^{0}∥}{μ_{0}} ≤ \frac{∥r^{k}∥}{μ_{k}} ≤ ρ̄\frac{∥r^{0}∥}{μ_{0}}.  (6)

For feasible interior-point algorithms, we have r^{k} = 0, k = 0, 1, ⋯, and so (d) is trivially satisfied.

Infeasible-interior-point algorithms that fit the standard framework include those of Zhang [11] and Wright [7]. This framework is broad enough to include most algorithms that have been proposed to date, including path-following and predictor-corrector algorithms (see, for example, Zhang [11] and Ji, Potra, and Huang [2]). Many algorithms use a central-path neighborhood different from (4); we use (4) partly because these alternative neighborhoods are usually subsets of our neighborhood when γ̄ is chosen sufficiently small.

We note that the analysis of most algorithms does not require Assumption 1 to hold in order to prove global linear convergence or polynomial complexity. It is usually needed only to obtain a local convergence rate that is faster than the global rate (which is usually Q-linear or two-step Q-linear).

The following observation (see Zhang [11, Proposition 3.3]) will be needed in subsequent results. If we define

ν_{0} = 1,  ν_{k} ≡ ∏_{j=0}^{k−1}(1 − α_{j}),  ∀k > 0,

then,

r^{k} = ν_{k}r^{0},  ∀k ≥ 0.  (7)

The first few results of this section give upper and lower bounds on the components x_{i}^{k} and y_{i}^{k} and their ratios. In Lemma 2.1, we show that (3) can be a consequence of other frequently-made assumptions.

Lemma 2.1 If either

(a) \underline{ρ} > 0 or

(b) there is a strictly feasible point for (1), and the sequence {μ_{k}} is bounded,

then (3) holds.

Proof. For case (a), our proof uses techniques similar to those of Mizuno [4, Lemma 3.3], Potra [6, Lemma 4.1], and Wright [8, Lemma 3.2]. Given any point (x̄, ȳ) with (x̄, ȳ) ≥ 0, ȳ = Mx̄ + q, and the starting point (x^{0}, y^{0}), we have

M(ν_{k}x^{0} + (1 − ν_{k})x̄ − x^{k})
= ν_{k}Mx^{0} + (1 − ν_{k})Mx̄ − Mx^{k}
= ν_{k}(y^{0} − q − r^{0}) + (1 − ν_{k})(ȳ − q) − (y^{k} − q − r^{k})
= ν_{k}y^{0} + (1 − ν_{k})ȳ − y^{k}.

Hence, by positive semidefiniteness of M, we have

0 ≤ (ν_{k}x^{0} + (1 − ν_{k})x̄ − x^{k})^{T} (ν_{k}y^{0} + (1 − ν_{k})ȳ − y^{k}).

Rearranging this expression, we obtain

ν_{k}(x^{0}^{T} y^{k} + x^{k}^{T} y^{0})  (8)
≤ ν_{k}^{2}x^{0}^{T} y^{0} + x^{k}^{T} y^{k} + ν_{k}(1 − ν_{k})(x^{0}^{T} ȳ + x̄^{T} y^{0}) + (1 − ν_{k})^{2}x̄^{T} ȳ − (1 − ν_{k})(x̄^{T} y^{k} + x^{k}^{T} ȳ).

Since (x̄, ȳ) ≥ 0, (x^{k}, y^{k}) ≥ 0, and (1 − ν_{k}) ≥ 0, the last term on the right-hand side is nonnegative, so we can drop it without affecting the inequality. Also, we can choose (x̄, ȳ) ∈ S to ensure x̄^{T} ȳ = 0, and after some manipulation (8) becomes

x^{0}^{T} y^{k} + x^{k}^{T} y^{0} ≤ ν_{k}x^{0}^{T} y^{0} + n\frac{μ_{k}}{ν_{k}} + (x^{0}^{T} ȳ + x̄^{T} y^{0})
≤ x^{0}^{T} y^{0} + n\frac{μ_{0}}{\underline{ρ}} + (x^{0}^{T} ȳ + x̄^{T} y^{0}),

where the second inequality follows from (6). If we define

C_{b} = max_{i=1,···,n} \{\frac{1}{min(x_{i}^{0}, y_{i}^{0})} [nμ_{0} (1 + \frac{1}{\underline{ρ}}) + (x^{0}^{T} ȳ + x̄^{T} y^{0})]\} ,

then (3) follows.

For (b), we choose (x̄, ȳ) to be a strictly feasible point, so that (x̄, ȳ) > 0. By rearranging (8), we obtain

(1 − ν_{k})(x̄^{T} y^{k} + x^{k}^{T} ȳ)  (9)
≤ ν_{k}^{2}x^{0}^{T} y^{0} + x^{k}^{T} y^{k} + ν_{k}(1 − ν_{k})(x^{0}^{T} ȳ + x̄^{T} y^{0}) + (1 − ν_{k})^{2}x̄^{T} ȳ − ν_{k}(x^{0}^{T} y^{k} + x^{k}^{T} y^{0}).

Again, the final term is nonnegative, so we can drop it without altering the inequality. We also use ν_{k} ∈ [0, 1] and

1 − ν_{k} ≥ 1 − ν_{1} = α_{0} > 0,  ∀k > 0  (10)

to write

x̄^{T} y^{k} + x^{k}^{T} ȳ ≤ \frac{1}{α_{0}} (x^{0}^{T} y^{0} + x^{k}^{T} y^{k}) + x^{0}^{T} ȳ + x̄^{T} y^{0} + x̄^{T} ȳ.

Boundedness of x^{k}^{T} y^{k} can now be used to obtain (3).

Our next result gives bounds on components of x_{B}^{k} and y_{N}^{k}.

Lemma 2.2 There is a positive constant C_{1} such that for all k ≥ 0,

i ∈ B ⇒ y_{i}^{k} ≤ C_{1}μ_{k},  x_{i}^{k} ≥ γ̄/C_{1};  (11)
i ∈ N ⇒ x_{i}^{k} ≤ C_{1}μ_{k},  y_{i}^{k} ≥ γ̄/C_{1}.  (12)

Proof. Again, we choose (x̄, ȳ) in (8) to be a solution of (1) for which x̄_{i} > 0 for i ∈ B and ȳ_{i} > 0 for i ∈ N. From (9) we have

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

Dropping the (nonnegative) final term on the right-hand side and using ν_{k} ∈ [0, 1], we have

(1 − ν_{k})(x̄^{T} y^{k} + x^{k}^{T} ȳ) ≤ ν_{k}nμ_{0} + nμ_{k} + ν_{k}(1 − ν_{k})(x^{0}^{T} ȳ + x̄^{T} y^{0}).  (13)

Using (10), (13), and (6), we can write

x̄^{T} y^{k} + x^{k}^{T} ȳ ≤ \frac{1}{α_{0}} (ν_{k}nμ_{0} + nμ_{k}) + ν_{k}(x^{0}^{T} ȳ + x̄^{T} y^{0})
≤ \frac{1}{α_{0}} (nρ̄μ_{k} + nμ_{k}) + \frac{ρ̄}{μ_{0}} μ_{k}(x^{0}^{T} ȳ + x̄^{T} y^{0})
≤ C̄_{1}μ_{k},

where C̄_{1} is defined in an obvious way. The upper bounds on x_{N}^{k} and y_{B}^{k} follow when we define

C_{1} = C̄_{1} max (max_{i∈B} \frac{1}{x̄_{i}} , max_{i∈N} \frac{1}{ȳ_{i}} ) .

To prove the remaining inequalities, we use (4). Taking i ∈ B, we have

x_{i}^{k}y_{i}^{k} ≥ γ̄μ_{k} ⇒ x_{i}^{k} ≥ \frac{γ̄μ_{k}}{y_{i}^{k}} ≥ \frac{γ̄μ_{k}}{C_{1}μ_{k}} = \frac{γ̄}{C_{1}} .

The lower bound on y_{i}^{k}, i ∈ N, is proved analogously.

The following consequence of a result of Mangasarian and Shiau [3] bounds the distance to the solution set in terms of μ_{k}.

Lemma 2.3 There is a positive constant C_{2} such that for all k ≥ 0,

dist((x^{k}, y^{k}), S) ≤ C_{2}μ_{k}^{1/2}.  (14)

Proof. From Mangasarian and Shiau [3, Theorem 2.7], there exist positive constants C_{2a} and C_{2b} such that

min_{x̄ | (x̄,Mx̄+q)∈S} ∥x̄ − x^{k}∥_{∞}  (15)

≤ C_{2a} \{∥[ \begin{matrix}(x^{k})^{T} (Mx^{k} + q) \\ −Mx^{k} − q \\ −x^{k} \end{matrix} ]_{+}∥_{2} + [(x^{k})^{T} (Mx^{k} + q) + C_{2b} ∥[ \begin{matrix}−Mx^{k} − q \\ −x^{k} \end{matrix} ]_{+}∥_{2}]^{1/2}\} .

In our case, Mx^{k} + q = y^{k} − r^{k} and (x^{k}, y^{k}) > 0, so by (6),

\{(x^{k})^{T} (Mx^{k} + q)i\}_{+} ≤ |(x^{k})^{T} (Mx^{k} + q)| ≤ (x^{k})^{T} y^{k} + ∥x^{k}∥∥r^{k}∥ ≤ (n + C_{b}ρ̄\frac{∥r^{0}∥}{μ_{0}}) μ_{k},

∥(−Mx^{k} − q)_{+}∥ = ∥(r^{k} − y^{k})_{+}∥ ≤ ∥r_{+}^{k}∥ ≤ ∥r^{k}∥ ≤ C_{b}ρ̄\frac{∥r^{0}∥}{μ_{0}}μ_{k},

∥(−x^{k})_{+}∥ = 0.

Substitution of these inequalities into (15) indicates that

min_{x̄ | (x̄,Mx̄+q)∈S} ∥x̄ − x^{k}∥_{∞} = O(μ_{k}^{1/2} + μ_{k}).  (16)

Now

dist((x^{k}, y^{k}), S) = min_{(x̄,ȳ)∈S} ∥[ \begin{matrix}x̄ − x^{k} \\ ȳ − y^{k} \end{matrix} ]∥

= min_{x̄ | (x̄,Mx̄+q)∈S} ∥[ \begin{matrix}x̄ − x^{k} \\ (Mx̄ + q) − (Mx^{k} + q + r^{k}) \end{matrix} ]∥
≤ min_{x̄ | (x̄,Mx̄+q)∈S}(1 + ∥M∥)∥x̄ − x^{k}∥ + ∥r^{k}∥
≤ (1 + ∥M∥)\sqrt{n} min_{x̄ | (x̄,Mx̄+q)∈S} ∥x̄ − x^{k}∥_{∞} + ρ̄∥r^{0}∥μ_{k}/μ_{0}
= O(μ_{k}^{1/2} + μ_{k}),

giving the result.

The next two results give bounds on the ratios of components of x^{k} to components of y^{k}.

Lemma 2.4 If C_{2} is the constant from Lemma 2.3, then for all k ≥ 0,

\frac{γ̄}{C_{2}} μ_{k}^{1/2} ≤ x_{i}^{k} ≤ C_{2}μ_{k}^{1/2},  \frac{γ̄}{C_{2}} μ_{k}^{1/2} ≤ y_{i}^{k} ≤ C_{2}μ_{k}^{1/2},  ∀i ∈ J,  (17)

and

\frac{γ̄}{C_{2}^{2}} ≤ \frac{x_{i}^{k}}{y_{i}^{k}} ≤ \frac{C_{2}^{2}}{γ̄}  ∀i ∈ J.  (18)

Proof. We prove the bounds (17) only for x_{i}^{k}, i ∈ J, since the results for y_{i}^{k} are similar. Since x_{i}^{∗} = y_{i}^{∗} = 0 for all i ∈ J, we have

x_{i}^{k} = x_{i}^{k} − x_{i}^{∗} ≤ dist((x^{k}, y^{k}), S) ≤ C_{2}μ_{k}^{1/2},

which gives the upper bounds. For the lower bounds, we use (4) to write

x_{i}^{k} ≥ \frac{γ̄μ_{k}}{y_{i}^{k}} ≥ \frac{γ̄μ_{k}}{C_{2}μ_{k}^{1/2}} = \frac{γ̄}{C_{2}} μ_{k}^{1/2}.

The bounds (18) follow immediately from (17).

Lemma 2.5 If C_{1} is the constant from Lemma 2.2, then for all k ≥ 0,

\frac{x_{i}^{k}}{y_{i}^{k}} ≤ \frac{C_{1}^{2}}{γ̄} μ_{k},  ∀i ∈ N,  (19)
\frac{y_{i}^{k}}{x_{i}^{k}} ≤ \frac{C_{1}^{2}}{γ̄} μ_{k},  ∀i ∈ B.  (20)

(21)

Proof. The proof follows immediately from (11) and (12).

It is possible to prove global linear convergence and polynomial complexity for a number of algorithms that fit the standard framework without assuming nondegeneracy (see, for example, Zhang [11], Wright [7], and Ji, Potra, and Huang [2]). However, Assumption 1 is used to prove Q-superlinear convergence in Ye and Anstreicher [10], Ji, Potra, and Huang [2], and Wright [7]. We are therefore led to pose the question, Is nondegeneracy necessary for superlinear convergence? The following two results resolve this question in the affirmative.

Theorem 2.6 Suppose that J ≠ ∅. Then there is a constant ϵ > 0 such that μ_{k+1}/μ_{k} ≥ ϵ for all k sufficiently large.

Proof. Suppose for contradiction that there is an infinite subsequence K such that μ_{k+1} = o(μ_{k}) for k ∈ K. Taking any index i ∈ J, we have from (17) that

0 < x_{i}^{k+1} ≤ C_{2}μ_{k+1}^{1/2} = o(μ_{k}^{1/2}),
0 < y_{i}^{k+1} ≤ C_{2}μ_{k+1}^{1/2} = o(μ_{k}^{1/2}).

Now

x_{i}^{k+1} = x_{i}^{k} + α_{k}Δx_{i}^{k} ⇒ α_{k}Δx_{i}^{k} = −x_{i}^{k} + x_{i}^{k+1} = −x_{i}^{k} + o(μ_{k}^{1/2})

and, similarly,

α_{k}Δy_{i}^{k} = −y_{i}^{k} + o(μ_{k}^{1/2}).

From (5), we have

y_{i}^{k}Δx_{i}^{k} + x_{i}^{k}Δy_{i}^{k} = −x_{i}^{k}y_{i}^{k} + σ_{k}μ_{k}
⇒ y_{i}^{k}(α_{k}Δx_{i}^{k}) + x_{i}^{k}(α_{k}Δy_{i}^{k}) = −α_{k}x_{i}^{k}y_{i}^{k} + σ_{k}α_{k}μ_{k}
⇒ −2x_{i}^{k}y_{i}^{k} + o(μ_{k}^{1/2})x_{i}^{k} + o(μ_{k}^{1/2})y_{i}^{k} = −α_{k}x_{i}^{k}y_{i}^{k} + σ_{k}α_{k}μ_{k}.

If we divide this last expression by x_{i}^{k}y_{i}^{k}, we obtain

−2 + \frac{o(μ_{k}^{1/2})}{x_{i}^{k}} + \frac{o(μ_{k}^{1/2})}{y_{i}^{k}} = −α_{k} + σ_{k}α_{k} \frac{μ_{k}}{x_{i}^{k}y_{i}^{k}} ,

and using the lower bound in (17) for x_{i}^{k} and y_{i}^{k}, we conclude that

2 − α_{k} = −σ_{k}α_{k} \frac{μ_{k}}{x_{i}^{k}y_{i}^{k}} + o(1) ≤ −σ_{k}α_{k}/n + o(1).

If we take the limit for k ∈ K, k → ∞, 2 − α_{k} is bounded below by 1, while the right-hand side of the above inequality will eventually be less than any positive constant. This gives a contradiction, and therefore we cannot have μ_{k+1} = o(μ_{k}).

Corollary 2.7 When the assumptions of Theorem 2.6 hold, the sequence {μ_{k}} cannot exhibit ℓ-step superlinear convergence to zero, for any integer ℓ ≥ 1.

Proof. The proof follows from μ_{k+ℓ}/μ_{k} ≥ ϵ^{ℓ} for all sufficiently large k.

As well as precluding Q-superlinear convergence of infeasible-interior-point algorithms such as the one in Wright [7], Corollary 2.7 shows that 2-step superlinear convergence is not possible for predictor-corrector algorithms such as the one in Ji, Potra, and Huang [2] applied to degenerate problems.

3 Infeasible Algorithms: Finite Termination

In this section, we propose a technique for estimating the index sets B, N, and J and for performing a projection from the current iterate onto the solution set S. We prove that finite termination can be achieved from any sufficiently advanced iterate (that is, when k is sufficiently large). The projection scheme is similar to that of Ye [9] for linear programming.

In the interests of generality, we make a number of assumptions on the iterates (x^{k}, y^{k}) themselves, rather than on the algorithm used to generate them. This strategy allows the analysis of this section to be applied to algorithms such as the one in [8], in addition to algorithms that fall within the standard framework of Section 2.

We require the following bounds to be satisfied for all sufficiently large k:

i ∈ J ⇒ \frac{1}{C_{3}} ≤ \frac{x_{i}^{k}}{y_{i}^{k}} ≤ C_{3},  0 < x_{i}^{k} ≤ C_{2}μ_{k}^{1/2},  0 < y_{i}^{k} ≤ C_{2}μ_{k}^{1/2},  (22a)

i ∈ B ⇒ 0 < \frac{y_{i}^{k}}{x_{i}^{k}} ≤ C_{4}μ_{k},  x_{i}^{k} ≥ C_{5},  (22b)

i ∈ N ⇒ 0 < \frac{x_{i}^{k}}{y_{i}^{k}} ≤ C_{4}μ_{k},  y_{i}^{k} ≥ C_{5},  (22c)
∥(x^{k}, y^{k})∥ ≤ C_{b},  (22d)
\frac{∥r^{k}∥}{μ_{k}} ≤ ρ̄\frac{∥r^{0}∥}{μ_{0}},  (22e)
lim_{k} μ_{k} = 0,  (22f)

where C_{3}, C_{4}, C_{5}, C_{b}, and ρ̄ are positive constants with C_{3} > 1 and C_{4} > 1.

The inequalities (22a), (22b), and (22c) immediately suggest a scheme for estimating B, N, and J. We define estimates B_{k}, N_{k}, and J_{k} by

B_{k} = \{i | \frac{y_{i}^{k}}{x_{i}^{k}} ≤ min(1/2, μ_{k}^{1/2})\} ,

N_{k} = \{i | \frac{x_{i}^{k}}{y_{i}^{k}} ≤ min(1/2, μ_{k}^{1/2})\} ,
J_{k} = {1, 2, ⋯, n}\(B_{k} ∪ N_{k}).

Lemma 3.1 For all k sufficiently large, we have

B_{k} = B,  N_{k} = N,  J_{k} = J.

Proof. Choose a positive integer K_{1} such that for k ≥ K_{1} the bounds (22) hold and, in addition,

μ_{k} < min(1/(2C_{4}^{2}), 1/C_{3}^{2}).

Then

i ∈ B ⇒ \frac{y_{i}^{k}}{x_{i}^{k}} ≤ C_{4}μ_{k} ≤ (C_{4}μ_{k}^{1/2})μ_{k}^{1/2} ≤ μ_{k}^{1/2},

and, since C_{4} > 1,

i ∈ B ⇒ \frac{y_{i}^{k}}{x_{i}^{k}} ≤ C_{4}μ_{k} ≤ 1/(2C_{4}) < 1/2.

Therefore i ∈ B_{k}. Similar logic shows that i ∈ N ⇒ i ∈ N_{k} for k ≥ K_{1}. For i ∈ J, we have

\frac{x_{i}^{k}}{y_{i}^{k}} ≥ \frac{1}{C_{3}} > μ_{k}^{1/2} ⇒ i ∉ N_{k},

while

\frac{y_{i}^{k}}{x_{i}^{k}} ≥ \frac{1}{C_{3}} > μ_{k}^{1/2} ⇒ i ∉ B_{k}.

Therefore i ∈ J_{k}, and the proof is complete.

On later iterations, when the index set estimates B_{k}, N_{k}, and J_{k} have stopped fluctuating from one iteration to the next, the following projection subproblem can be solved in an attempt to find an exact solution to (1). The problem is an equality-constrained quadratic program and hence is easier to solve than the original LCP.

Problem P (k):  min_{x_{B_{k}} ,y_{N_{k}}} \frac{1}{2}∥x_{B_{k}} − x_{B_{k}}^{k}∥^{2} + \frac{1}{2}∥y_{N_{k}} − y_{N_{k}}^{k}∥^{2}

subject to

0 = M_{B_{k},B_{k}}x_{B_{k}} + q_{B_{k}},
0 = M_{J_{k},B_{k}}x_{B_{k}} + q_{J_{k}},
y_{N_{k}} = M_{N_{k},B_{k}}x_{B_{k}} + q_{N_{k}}.

Let (x̄_{B_{k}}^{k}, ȳ_{N_{k}}^{k}) denote the optimal solution of this problem. A candidate solution (x̄^{k}, ȳ^{k}) of (1) is obtained by setting

x̄_{i}^{k} = \{ \begin{matrix}(x̄_{B_{k}})_{i} & i ∈ B_{k} \\ 0 & otherwise \end{matrix} ,  ȳ_{i}^{k} = \{ \begin{matrix}(ȳ_{N_{k}})_{i} & i ∈ N_{k} \\ 0 & otherwise \end{matrix} .  (23)

Theorem 3.2 For all k sufficiently large, the solution of Problem P (k) and the construction (23) yield a solution of (1).

Proof. Assume that k is large enough that the bounds (22) and the result of Lemma 3.1 hold. (We can subsequently refer to B_{k}, N_{k}, and J_{k} as B, N, and J without confusion.) For any (x^{∗}, y^{∗}) ∈ S, we have that (x_{B}^{∗}, y_{N}^{∗}) satisfies the equality constraints in Problem P (k), so the feasible set for P (k) is nontrivial. Using a change of variable

x̃_{B}^{k} = x_{B} − x_{B}^{k}, ỹ_{N}^{k} = y_{N} − y_{N}^{k},

we can reformulate P (k) as

min_{x̃_{B}^{k},ỹ_{N}^{k}} \frac{1}{2}∥x̃_{B}∥^{2} + \frac{1}{2}∥ỹ_{N}∥^{2},  (24)

subject to

−M_{B,B}x̃_{B}^{k} = q_{B} + M_{B,B}x_{B}^{k},
− M_{J,B}x̃_{B}^{k} = q_{J} + M_{J,B}x_{B}^{k},  (25)
ỹ_{N}^{k} − M_{N,B}x̃_{B}^{k} = q_{N} + M_{N,B}x_{B}^{k} − y_{N}^{k}.

Since the constraints (25) are consistent, Hoffman’s lemma [1] implies that there is a sequence of vector pairs {(x̃_{B}^{k}, ỹ_{N}^{k})} such that (x̃_{B}^{k}, ỹ_{N}^{k}) is feasible for (25) and

∥[ \begin{matrix} & x̃_{B}^{k} \\  & ỹ_{N}^{k} \end{matrix} ]∥ = O (∥[ \begin{matrix}q_{B} + M_{B,B}x_{B}^{k} \\ q_{J} + M_{J,B}x_{B}^{k} \\ q_{N} + M_{N,B}x_{B}^{k} − y_{N}^{k} \end{matrix} ]∥) .  (26)

Since (x̄_{B}^{k} − x_{B}^{k}, ȳ_{N}^{k} − y_{N}^{k}) is optimal for problem (24)–(25), we have

∥[ \begin{matrix}x̄_{B}^{k} − x_{B}^{k} \\ ȳ_{N}^{k} − y_{N}^{k} \end{matrix} ]∥ ≤ ∥[ \begin{matrix} & x̃_{B}^{k} \\  & ỹ_{N}^{k} \end{matrix} ]∥ .  (27)

Now, using the bounds (22), we obtain

∥[ \begin{matrix}q_{B} + M_{B,B}x_{B}^{k} \\ q_{J} + M_{J,B}x_{J}^{k} \\ q_{N} + M_{N,B}x_{B}^{k} − y_{N}^{k} \end{matrix} ]∥

≤ ∥q + Mx^{k} − y^{k}∥ + ∥[ \begin{matrix}M_{B,N}x_{N}^{k} + M_{B,J}x_{J}^{k} − y_{B}^{k} \\ M_{J,N}x_{N}^{k} + M_{J,J}x_{J}^{k} − y_{J}^{k} \\ M_{N,N}x_{N}^{k} + M_{N,J}x_{J}^{k} \end{matrix} ]∥
= ∥r^{k}∥ + O(∥x_{N}^{k}∥ + ∥x_{J}^{k}∥ + ∥y_{B}^{k}∥ + ∥y_{J}^{k}∥)
= O(μ_{k}^{1/2}).  (28)

It follows from (26)–(28) that there is a constant C_{8} > 0 such that

∥x̄_{B}^{k} − x_{B}^{k}∥ ≤ C_{8}μ_{k}^{1/2}, ∥ȳ_{N}^{k} − y_{N}^{k}∥ ≤ C_{8}μ_{k}^{1/2}.

Hence, we have

x̄_{B}^{k} ≥ x_{B}^{k} − C_{8}μ_{k}^{1/2}e ≥ C_{5}e − C_{8}μ_{k}^{1/2}e,
ȳ_{N}^{k} ≥ y_{N}^{k} − C_{8}μ_{k}^{1/2}e ≥ C_{5}e − C_{8}μ_{k}^{1/2}e,

and for k sufficiently large, we conclude that (x̄_{B}^{k}, ȳ_{N}^{k}) > 0. The construction (23) then yields a vector pair (x̄^{k}, ȳ^{k}) that solves (1).

4 Feasible Algorithms

In this section we consider only feasible algorithms (r^{k} ≡ 0) in which all iterates lie in the relative interior of the feasible region. Condition (6) of the standard framework of Section 2 becomes redundant. Theorem 2.6 shows that no algorithm from this framework can converge superlinearly when J ≠ ∅. Hence, the best we can hope to show is that the algorithm has a linear rate of convergence that is not too slow. Specifically, if we define

Λ =^{△} lim sup_{k→∞} \frac{x^{k+1}^{T} y^{k+1}}{x^{k}^{T} y^{k}} = lim sup_{k→∞} \frac{μ_{k+1}}{μ_{k}} ,  (29)

then we hope that the value Λ ∈ (0, 1) is not too close to 1.

The main result of this section — Theorem 4.1 — assumes that σ_{k} converges to 0 and gives a formula for the linear rate of convergence in terms of just the current iterate and the current step size.

Theorem 4.1 Assume that r^{k} = 0 for all k ≥ 0 and lim_{k→∞} σ_{k} = 0. Then Λ defined by (29) satisfies

Λ = lim sup_{k→∞} [1 − α_{k} + α_{k}^{2} \frac{x_{J}^{k}^{T} y_{J}^{k}}{4 x^{k}^{T} y^{k}} ] .  (30)

The proof of this theorem is postponed until we have proved several preliminary results. Note that if α_{k} converges to 1, then the rate of convergence becomes simply

Λ = lim sup_{k→∞} [\frac{x_{J}^{k}^{T} y_{J}^{k}}{4 x^{k}^{T} y^{k}} ] ≤ \frac{1}{4}.

This expression indicates a reasonably fast linear convergence rate, but we have not been able to design an algorithm that achieves this rate since it is difficult to enforce the condition lim_{k→∞} α_{k} = 1.

We now derive lower and upper bounds on the linear rate of convergence. The following result is a consequence of Theorem 4.1.

Corollary 4.2 Assume that r^{k} = 0 for all k ≥ 0 and lim_{k→∞} σ_{k} = 0. If lim inf_{k→∞} α_{k} = ᾱ, then

1 − ᾱ + \frac{1}{4}ᾱ^{2}L ≤ lim sup_{k→∞} \frac{x^{k+1}^{T} y^{k+1}}{x^{k}^{T} y^{k}} ≤ 1 − ᾱ + \frac{1}{4}ᾱ^{2}U,  (31)

where

L = γ̄ \frac{|J|}{n}, U = γ̄ \frac{|J|}{n} + (1 − γ̄) < 1.

Proof. Using relation (4), we can easily show that

L ≤ \frac{x_{J}^{k}^{T} y_{J}^{k}}{x^{k}^{T} y^{k}} ≤ U.  (32)

Hence,

lim sup_{k→∞} \frac{x^{k+1}^{T} y^{k+1}}{x^{k}^{T} y^{k}} ≤ lim sup_{k→∞} [1 − α_{k} + \frac{1}{4}α_{k}^{2}U] = 1 − ᾱ + \frac{1}{4}ᾱ^{2}U,

where the last equality follows because the quadratic 1 − α + (U/4)α^{2} is decreasing on [0, 1]. We have therefore proved the second inequality of (31); the first inequality follows by a similar argument.

We now concentrate our efforts on the proof of Theorem 4.1.

Lemma 4.3 Let a ∈ IR_{++}^{p}, b ∈ IR_{++}^{p}, and z ∈ IR^{p} be given. Let H ⊆ IR^{2p} be a subspace with the property that

(u, v) ∈ H ⇒ u^{T} v ≥ 0.  (33)

Then, the relations

au + bv = z,  (u, v) ∈ H,  (34)

have at most one solution (u, v).

Proof. Assume for contradiction that (u^{1}, v^{1}) ∈ IR^{2p} and (u^{2}, v^{2}) ∈ IR^{2p} are two solutions of (34). We have (u^{1} − u^{2}, v^{1} − v^{2}) ∈ H since H is a subspace. Therefore, in view of (33), we have

(u^{1} − u^{2})^{T} (v^{1} − v^{2}) ≥ 0.  (35)

On the other hand, we obtain

a (u^{1} − u^{2}) + b (v^{1} − v^{2}) = 0,  (36)

which in turn implies

v^{1} − v^{2} = −b^{−1}a (u^{1} − u^{2}), u^{1} − u^{2} = −a^{−1}b (v^{1} − v^{2}).

Combining these relations with relation (35), we obtain

−∥(b^{−1}a)^{1/2}(u^{1} − u^{2})∥ ≥ 0, −∥(a^{−1}b)^{1/2}(v^{1} − v^{2})∥ ≥ 0.

It follows from these last two inequalities that u^{1} = u^{2} and v^{1} = v^{2}, so our result is proved.

Lemma 4.4 There hold

Δx_{N}^{k} = O(μ_{k}),  Δy_{B}^{k} = O(μ_{k});  (37)
Δx_{J}^{k} = O(μ_{k}^{1/2}),  Δy_{J}^{k} = O(μ_{k}^{1/2});  (38)
Δx_{B}^{k} = O(μ_{k}^{1/2} + σ_{k}),  Δy_{N}^{k} = O(μ_{k}^{1/2} + σ_{k}).  (39)

Proof. The proof is a modification of earlier results of Ye and Anstreicher [10] and Wright [7, 8]. For completeness, we include it in the appendix.

Lemma 4.5 Assume that r^{k} = 0 for all k ≥ 0 and lim_{k→∞} σ_{k} = 0. Then

i ∈ J ⇒ lim_{k→∞} \frac{Δx_{i}^{k}}{x_{i}^{k}} = −\frac{1}{2},  lim_{k→∞} \frac{Δy_{i}^{k}}{y_{i}^{k}} = −\frac{1}{2}.

Proof. It follows from Lemmas 2.4 and 4.4 that the sequence {((x_{J}^{k})^{−1}Δx_{J}^{k}, (y_{J}^{k})^{−1}Δy_{J}^{k})} is bounded. Let (w^{x}, w^{y}) be an accumulation point of this sequence, so that there is an infinite subsequence K with

lim_{k∈K}((x_{J}^{k})^{−1}Δx_{J}^{k}, (y_{J}^{k})^{−1}Δy_{J}^{k}) = (w^{x}, w^{y}).  (40)

By further restriction of K if necessary, we use Lemmas 2.4 and 4.4 again to deduce that there are vector pairs (l^{x}, l^{y}) and (δ^{x}, δ^{y}) such that

lim_{k∈K} \frac{1}{μ_{k}^{1/2}} (x_{J}^{k}, y_{J}^{k}) = (l^{x}, l^{y}) > 0,  (41)

and

lim_{k∈K} \frac{1}{μ_{k}^{1/2}} (Δx_{J}^{k}, Δy_{J}^{k}) = (δ^{x}, δ^{y}).  (42)

By using Lemma 4.4, the fact that lim_{k→∞} σ_{k} = 0, and the relations

x_{J}^{k}Δy_{J}^{k} + y_{J}^{k}Δx_{J}^{k} = σ_{k}μ_{k}e − x_{J}^{k}y_{J}^{k},
Δy^{k} − MΔx^{k} = 0,

we can easily verify that

l^{x}δ^{y} + l^{y}δ^{x} = −l^{x}l^{y},  (43)
I_{.J}δ^{y} − M_{.J}δ^{x} ∈ Range[I_{.N}, −M_{.B}].  (44)

Let (x̄, ȳ) be an arbitrary solution of (1). Then, we have

0 = y^{k} − Mx^{k} − q
= y^{k} − Mx^{k} − (ȳ − Mx̄)
= (y^{k} − ȳ) − M(x^{k} − x̄).

This relation and the fact that x̄_{J∪N} = 0 and ȳ_{B∪J} = 0 then imply

I_{.B} \frac{y_{B}^{k}}{μ_{k}^{1/2}} + I_{.J} \frac{y_{J}^{k}}{μ_{k}^{1/2}} − M_{.J} \frac{x_{J}^{k}}{μ_{k}^{1/2}} − M_{.N} \frac{x_{N}^{k}}{μ_{k}^{1/2}} ∈ Range[I_{.N}, −M_{.B}].  (45)

It follows from Lemma 2.2 that the first and fourth terms in the left-hand side of relation (45) converge to 0, so we obtain

I_{.J}l^{y} − M_{.J}l^{x} ∈ Range[I_{.N}, −M_{.B}].  (46)

It is easy to verify that the subspace H of the vectors (u, v) ∈ IR^{2|J|} satisfying relation (46) (with (l^{x}, l^{y}) replaced by (u, v)) has the property expressed by relation (33). Hence, by Lemma 4.3, the system defined by (43) and (44) has at most one solution. From (46), it is clear that

(δ^{x}, δ^{y}) = −\frac{1}{2}(l^{x}, l^{y})

solves (43) and (44) and is therefore the unique solution. Hence, from (40), (41), and (42), we obtain

(w^{x}, w^{y}) = ((l^{x})^{−1}δ^{x}, (l^{y})^{−1}δ^{y}) = −\frac{1}{2}(e, e).

Since every accumulation point of \{((x_{J}^{k})^{−1}Δx_{J}^{k}, (y_{J}^{k})^{−1}Δy_{J}^{k})\} is equal to −\frac{1}{2}(e, e), the result follows.

We are now in a position to give the proof of Theorem 4.1.

Proof of Theorem 4.1. From the standard framework of Section 2, it is straightforward to see that

\frac{x^{k+1}^{T} y^{k+1}}{x^{k}^{T} y^{k}} = 1 − α_{k} + σ_{k}α_{k} + α_{k}^{2} \frac{Δx^{k}^{T} Δy^{k}}{x^{k}^{T} y^{k}} .  (47)

Observe that relation (30), and therefore Theorem 4.1, follows immediately from (47) once we show that

\frac{Δx^{k}^{T} Δy^{k}}{x^{k}^{T} y^{k}} = \frac{x_{J}^{k}^{T} y_{J}^{k}}{4 x^{k}^{T} y^{k}} + o(1).  (48)

To show (48), define

t_{x}^{k} = (x_{J}^{k})^{−1}Δx_{J}^{k} − \frac{1}{2}e,  t_{y}^{k} = (y_{J}^{k})^{−1}Δy_{J}^{k} − \frac{1}{2}e.  (49)

By Lemma 4.5, we know that t_{x}^{k} = o(1) and t_{y}^{k} = o(1). Using (49), we have

Δx_{i}^{k}Δy_{i}^{k} = x_{i}^{k}y_{i}^{k} [ \frac{1}{2} + (t_{x}^{k})_{i}] [ \frac{1}{2} + (t_{y}^{k})_{i}] = \frac{x_{i}^{k}y_{i}^{k}}{4} [1 + q_{i}^{k}],  (50)

where

q_{i}^{k} ≡ 2 [ (t_{x}^{k})_{i} + (t_{y}^{k})_{i} + 2(t_{x}^{k})_{i}(t_{y}^{k})_{i}] = o(1).  (51)

From (50), we obtain

\frac{Δx_{J}^{k}^{T} Δy_{J}^{k}}{x^{k}^{T} y^{k}} = \frac{x_{J}^{k}^{T} y_{J}^{k}}{4 x^{k}^{T} y^{k}} + ∑_{i∈J} q_{i}^{k} \frac{x_{i}^{k}y_{i}^{k}}{4x^{k}^{T} y^{k}} = \frac{x_{J}^{k}^{T} y_{J}^{k}}{4 x^{k}^{T} y^{k}} + o(1).  (52)

By using Lemma 4.4, we can easily verify that

\frac{Δx_{B}^{k}^{T} Δy_{B}^{k}}{x^{k}^{T} y^{k}} = o(1),  \frac{Δx_{N}^{k}^{T} Δy_{N}^{k}}{x^{k}^{T} y^{k}} = o(1).  (53)

Relation (48) now follows by combining (52) and (53).

5 Convergence of the Predictor-Corrector Algorithm

In this section we use the results of the preceding section to analyze the two-step linear rate of convergence of the predictor corrector algorithm (see, for example, [2], [5], and [10]) when applied to degenerate LCPs. We show that this rate is less than or equal to 1−c/Q^{1/4} where Q ≤ min{|J|, n − |J|} and c is a positive constant that is not too small. This result shows that if |J| contains only a few indices, or all but a few indices, then we can attain a linear rate of convergence that is not too slow.

Although the predictor-corrector algorithm fits the standard framework of Section 2 (and can be described accordingly by letting predictor steps be taken at even values of k and corrector steps at odd values), it is more convenient to use the description below, which closely follows [2]. For a given constant β ∈ (0, 1), define

N (β) ≡ {(x, y) ≥ 0 | y = Mx + q, ∥Xy − μe∥ ≤ βμ},

where, as earlier, μ ≡ x^{T} y/n. For (x, y) > 0 with y = Mx + q, the following system defines a direction (Δx, Δy) which is used in the description of the predictor-corrector algorithm:

MΔx − Δy = 0,  (54a)
yΔx + xΔy = σμe − xy,  (54b)

where σ ∈ [0, 1].

Predictor-Corrector Algorithm: Let the constants β ∈ (0, 1/4] and τ ∈ (0, β], and a strictly feasible solution (x^{0}, y^{0}) ∈ N (β) be given. Set k = 0, and go to step 1.

(1) Compute the predictor step (Δx^{k}, Δy^{k}) by solving system (54) with (x, y) = (x^{k}, y^{k}) and σ_{k} = 0.

(2) Compute the step size α_{k} > 0 by

α_{k} ≡ max{α ∈ [0, 1] | (x^{k} + α̃Δx^{k}, y^{k} + α̃Δy^{k}) ∈ N (β + τ), ∀α̃ ∈ [0, α]},

and set (x̂^{k}, ŷ^{k}) ≡ (x^{k}, y^{k}) + α_{k}(Δx^{k}, Δy^{k}).

(3) Compute the corrector step (Δx̂^{k}, Δŷ^{k}) by solving system (54) with (x, y) = (x̂^{k}, ŷ^{k}) and σ_{k} = 1.

(4) Compute the new iterate as

(x^{k+1}, y^{k+1}) = (x̂^{k}, ŷ^{k}) + (Δx̂^{k}, Δŷ^{k}).

Set k = k + 1, and go to step 1. We now state several properties of the predictor-corrector algorithm. The main properties

of the centering steps are given by the following result. Lemma 5.1 Let μ̂_{k} ≡ (x̂^{k})^{T} ŷ^{k}/n for all k ≥ 0. The following statements hold for every k ≥ 0:

(a) ∥Δx̂^{k}Δŷ^{k}∥ ≤ βμ̂^{k};

(b) (Δx̂^{k})^{T} Δŷ^{k} ≤ μ̂^{k}/8;

(c) (x^{k+1}, y^{k+1}) = (x̂^{k}, ŷ^{k}) + (Δx̂^{k}, Δŷ^{k}) ∈ N (β);

(d) (x^{k+1})^{T} y^{k+1} = (x̂^{k} + Δx̂^{k})^{T} (ŷ^{k} + Δŷ^{k}) ≤ [1 + 1/(8n)] (x̂^{k})^{T} ŷ^{k}. Proof. Statements (a) and (c) follows from arguments similar to the ones used in Lemma

2.3 of Ji, Potra, and Huang [2], while (b) and (d) are based on Lemma 3.1 of [2].

As a consequence of Lemma 5.1, we have the following result. Corollary 5.2 For all k ≥ 0, we have

\frac{(x^{k+1})^{T} y^{k+1}}{x^{k}^{T} y^{k}} ≤ (1 + \frac{1}{8n}) (1 − α_{k} + α_{k}^{2} \frac{(Δx^{k})^{T} Δy^{k}}{x^{k}^{T} y^{k}} ) .  (55)

Proof. The result follows immediately from Lemma 5.1(d) and the fact that (x̂^{k})^{T} ŷ^{k} = (x^{k} + α_{k}Δx^{k})^{T} (y^{k} + α_{k}Δy^{k}) = (1 − α_{k})x^{k}^{T} y^{k} + α_{k}^{2}(Δx^{k})^{T} Δy^{k}.

For the predictor steps, we have the following result. Lemma 5.3 The following statements hold for every k ≥ 0:

(a) ∥Δx^{k}Δy^{k}∥ ≤ (x^{k})^{T} y^{k}/2;

(b) (Δx^{k})^{T} Δy^{k} ≤ (x^{k})^{T} y^{k}/4; (c) The step size α_{k} satisfies

α_{k} ≥ (\frac{τ}{n})^{1/2}

and

α_{k} ≥ \frac{2}{1 + \sqrt{1 + 4∥z^{k}∥/τ}} ,  (56)

where

z^{k} ≡ \frac{1}{μ_{k}} (Δx^{k}Δy^{k} − \frac{(Δx^{k})^{T} Δy^{k}}{n} e);  (57)

(d) If τ ≥ 1/(16n), then

\frac{x^{k+1}^{T} y^{k+1}}{x^{k}^{T} y^{k}} ≤ (1 − \frac{\sqrt{τ}}{2\sqrt{n}}) = 1 − O(1/\sqrt{n}).  (58)

For a proof of Lemma 5.3, we refer the reader to Lemma 3.1, Lemma 3.2, and Theorem 3.1 of Ji, Potra and Huang [2] and Lemma 2.2 of Ye and Anstreicher [10].

It is well known that Lemma 5.3(d) implies polynomial convergence of the predictorcorrector algorithm. Our final goal — stated in the next theorem — is to provide an upper bound on the ratio of successive elements of the sequence {(x^{k})^{T} y^{k}} which is sharper than the one implied by relation (58).

Theorem 5.4 Assume that τ ≥ 1/(8n) and 0 < |J| < n. Then,

lim sup_{k→∞} \frac{x^{k+1}^{T} y^{k+1}}{x^{k}^{T} y^{k}} ≤ 1 − \frac{(τ/2)^{1/2}}{2Q^{1/4}} ,  (59)

where

Q ≡ \frac{(n − |J|)|J|}{n} ≤ min\{|J|, n − |J|\}.

The proof of Theorem 5.4 is postponed until we have proved the following preliminary result.

Lemma 5.5 Assume that n ≥ 2. For all k sufficiently large, we have

∥z^{k}∥ ≤ Q^{1/2}/2.  (60)

Proof. Define w^{k} ∈ IR^{n} as

w^{k} = ( \begin{matrix}w_{B}^{k} \\ w_{N}^{k} \\ w_{J}^{k} \end{matrix} ) ≡ ( \begin{matrix}0 \\ 0 \\ x_{J}^{k}s_{J}^{k} \end{matrix} ) .

Using Lemmas 4.4 and 4.5, we can easily verify that

∥μ_{k}z^{k}∥ = ∥Δx^{k}Δy^{k} − \frac{(Δx^{k})^{T} Δy^{k}}{n} e∥ = \frac{1}{4} ∥w^{k} − \frac{(x_{J}^{k})^{T} y_{J}^{k}}{n} e∥ + o(μ_{k}).  (61)

Since (x^{k}, y^{k}) ∈ N (β), a simple argument shows that

∥x_{J}^{k}y_{J}^{k} − \frac{(x_{J}^{k})^{T} y_{J}^{k}}{|J|} e∥ ≤ ∥x_{J}^{k}y_{J}^{k} − μ_{k}e∥ ≤ βμ_{k}.

Also,

(x_{J}^{k})^{T} y_{J}^{k} ≤ (1 + β)|J|μ_{k},

so we obtain

∥w^{k} − \frac{(x_{J}^{k})^{T} y_{J}^{k}}{n} e∥^{2}
= ∥x_{J}^{k}y_{J}^{k} − \frac{(x_{J}^{k})^{T} y_{J}^{k}}{n} e∥^{2} + (n − |J|) (\frac{(x_{J}^{k})^{T} y_{J}^{k}}{n} )^{2}
= ∥x_{J}^{k}y_{J}^{k} − \frac{(x_{J}^{k})^{T} y_{J}^{k}}{|J|} e∥^{2} + ∥\frac{(x_{J}^{k})^{T} y_{J}^{k}}{|J|} e − \frac{(x_{J}^{k})^{T} y_{J}^{k}}{n} e∥^{2} + (n − |J|) (\frac{(x_{J}^{k})^{T} y_{J}^{k}}{n} )^{2}
≤ (βμ_{k})^{2} + \frac{(n − |J|)^{2}}{n^{2}|J|} (x_{J}^{k}^{T} y_{J}^{k})^{2} + (n − |J|) (\frac{(x_{J}^{k})^{T} y_{J}^{k}}{n} )^{2}
= (βμ_{k})^{2} + \frac{n − |J|}{n|J|}(x_{J}^{k}^{T} y_{J}^{k})^{2}
≤ (βμ_{k})^{2} + \frac{n − |J|}{n|J|}|J|^{2}(1 + β)^{2}μ_{k}^{2}
≤ μ_{k}^{2}(1 + β)^{2} [β^{2} + \frac{(n − |J|)|J|}{n}]
≤ μ_{k}^{2}(1 + β)^{2} [β^{2} + Q] .

Hence, for all k sufficiently large, we have from (61) that

∥μ_{k}z^{k}∥ ≤ μ_{k} \frac{(1 + β)}{4} [β^{2} + Q]^{1/2} + o(μ_{k}) ≤ μ_{k} \frac{(1 + β)}{4} [β + Q^{1/2}] ≤ \frac{μ_{k}}{2} Q^{1/2},

where the last inequality follows from the fact that β ≤ 1/4 and Q^{1/2} ≥ 1/\sqrt{2}. We are now in a position to give the proof of Theorem 5.4.

Proof of Theorem 5.4. Using relations (56) and (60), we obtain

α_{k} ≥ \frac{2}{1 + \sqrt{1 + 2Q^{1/2}/τ}} ≥ \frac{2}{2\sqrt{2Q^{1/2}/τ}} = \frac{\sqrt{τ/2}}{Q^{1/4}},  (62)

where the second inequality can easily be verified by using τ ≤ 1/4 and Q^{1/2} ≥ 1/\sqrt{2}. Using relations (55), (62), and Lemma (5.3)(b), we conclude that for all k sufficiently large, we have

\frac{x^{k+1}^{T} y^{k+1}}{x^{k}^{T} y^{k}} ≤ (1 + \frac{1}{8n}) (1 − \frac{α_{k}}{2} )^{2}
≤ (1 + \frac{1}{8n} − \frac{α_{k}}{2} ) (1 − \frac{α_{k}}{2} )
≤ (1 − \frac{α_{k}}{2} )
≤ (1 − \frac{(τ/2)^{1/2}}{2Q^{1/4}} )

where the third inequality follows from (62) and the fact that τ ≥ 1/(8n) and Q^{1/2} ≤ n.

A Appendix

We prove Lemma 4.4 for the infeasible case r^{k} ≠ 0, although a proof for the feasible case would suffice for the purposes of Section 4. Before doing do, we prove some useful auxiliary results. First, we recall a result due to Ye and Anstreicher [10].

Lemma A.1 Let M be a positive semi-definite matrix, and partition M as

M = ( \begin{matrix}M_{PP} & M_{PQ} \\ M_{QP} & M_{QQ} \end{matrix} ) ,

where the pair of index sets P ⊆ {1, ⋯, n} and Q ⊆ {1, ⋯, n} forms a nontrivial partition of {1, ⋯, n}. Then

Range ( \begin{matrix}M_{PP} & M_{PQ} \\ 0 & I \end{matrix} ) = Range ( \begin{matrix}M_{PP}^{T} & M_{QP}^{T} \\ 0 & −I \end{matrix} )

As a consequence of Lemma A.1, we obtain the following result which will be explicitly used in the proof of Lemma 4.4.

Lemma A.2 There holds

Range ( \begin{matrix}M_{BB} & M_{BJ} & M_{BN} \\ 0 & 0 & I \end{matrix} ) = Range ( \begin{matrix}M_{BB}^{T} & M_{JB}^{T} & M_{NB}^{T} \\ 0 & 0 & −I \end{matrix} ) .  (63)

Proof. Applying Lemma A.1 with P = B ∪ J and Q = N, we obtain

Range ( \begin{matrix}M_{BB} & M_{BJ} & M_{BN} \\ M_{JB} & M_{JJ} & M_{JN} \\ 0 & 0 & I \end{matrix} ) = Range ( \begin{matrix}M_{BB}^{T} & M_{JB}^{T} & M_{NB}^{T} \\ M_{BJ}^{T} & M_{JJ}^{T} & M_{NJ}^{T} \\ 0 & 0 & −I \end{matrix} ) ,

which in turn immediately implies (63).

The following lemma is similar to Lemma 3.5 of Ye and Anstreicher [10] and Lemma 5.2 of Wright [7]. For this proof and the proof of Lemma 4.4, we drop the iteration index k on matrices and vectors for clarity, and define the diagonal matrix D = X^{−1/2}Y^{1/2}. (The principal submatrices D_{N} and D_{B} are defined in an obvious way.)

Lemma A.3 The vector pair (Δx_{B}, Δy_{N}) is the unique solution of the convex quadratic programming problem

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

subject to

M_{BB}w  = r_{B} − M_{BJ}Δx_{J} − M_{BN}Δx_{N} + Δy_{B},
M_{JB}w = r_{J} − M_{JJ}Δx_{J} − M_{JN}Δx_{N} + Δy_{J}, (65)
M_{NB}w − z = r_{N} − M_{NJ}Δx_{J} − M_{NN}Δx_{N}.

Proof. By the Karush-Kuhn-Tucker conditions, a candidate solution (w, z) is optimal if it is feasible with respect to (65) and, in addition,

( \begin{matrix}D_{B}^{2}w − σ_{k}μ_{k}X_{B}^{−1}e \\ D_{N}^{−2}z − σ_{k}μ_{k}Y_{N}^{−1}e \end{matrix} ) ∈ Range ( \begin{matrix}M_{BB}^{T} & M_{JB}^{T} & M_{NB}^{T} \\ 0 & 0 & −I \end{matrix} ) .  (66)

We prove the result by showing that (Δx_{B}, Δy_{N}) satisfies this condition. Clearly, (Δx_{B}, Δy_{N}) is feasible with respect to (65). Using (5), we have

D_{B}^{2}Δx_{B} − σ_{k}μ_{k}X_{B}^{−1}e
= −y_{B} − Δy_{B}
= −y_{B} − (M_{BB}Δx_{B} − r_{B} + M_{BN}Δx_{N} + M_{BJ}Δx_{J})
= −M_{BB}(x_{B} + Δx_{B}) − M_{BN}(x_{N} + Δx_{N}) − M_{BJ}(x_{J} + Δx_{J}),

and

D_{N}^{−2}Δy_{N} − σ_{k}μ_{k}Y_{N}^{−1}e = −(x_{N} + Δx_{N}).

Therefore,

( \begin{matrix}D_{B}^{2}w − σ_{k}μ_{k}X_{B}^{−1}e \\ D_{N}^{−2}z − σ_{k}μ_{k}Y_{N}^{−1}e \end{matrix} ) ∈ Range ( \begin{matrix}M_{BB} & M_{BJ} & M_{BN} \\ 0 & 0 & I \end{matrix} ) .  (67)

From relations (63) and (67), it follows that (Δx_{B}, Δy_{N}) satisfies (66).

We are now in a position to give the proof of Lemma 4.4.

Proof of Lemma 4.4. First, we show that

∥DΔx∥ = O(μ_{k}^{1/2}),  ∥D^{−1}Δy∥ = O(μ_{k}^{1/2}).  (68)

We introduce vector pairs (û, v̂) and (ū, v̄) that satisfy

[ \begin{matrix}M & −I \\ Y & X \end{matrix} ] [ û ; v̂ ] = [ \begin{matrix}r \\ 0 \end{matrix} ] ,  (69a)
[ \begin{matrix}M & −I \\ Y & X \end{matrix} ] [ ū ; v̄ ] = [ \begin{matrix}0 \\ −xy + σ_{k}μ_{k}e \end{matrix} ] ,  (69b)

so that

(Δx, Δy) = (û, v̂) + (ū, v̄).

It follows exactly as in [8, Lemma 3.3] that

∥Dū∥ = O(μ_{k}^{1/2}),  ∥D^{−1}v̄∥ = O(μ_{k}^{1/2}).  (70)

If we use (6) and the boundedness assumption (3), minor modifications to the proof of [8, Lemma 3.4] can be used to show

∥Dû∥ = O(μ_{k}^{1/2}),  ∥D^{−1}v̂∥ = O(μ_{k}^{1/2}).  (71)

The relations (68) follow from (70) and (71).

We now prove (37) for Δx_{N} (the proof for Δy_{B} is similar). Taking i ∈ N, we have from (68) and (19) that

|(\frac{y_{i}}{x_{i}} )^{1/2} Δx_{i}| ≤ ∥DΔx∥ = O(μ_{k}^{1/2})

⇒ |Δx_{i}| ≤ (\frac{x_{i}}{y_{i}} )^{1/2} O(μ_{k}^{1/2}) ≤ (\frac{C_{1}^{2}}{γ̄} μ_{k})^{1/2} O(μ_{k}^{1/2}) = O(μ_{k}),

as required.

To prove (38), we observe from (68) and (18) that for i ∈ J, we have

|(\frac{y_{i}}{x_{i}} )^{1/2} Δx_{i}| = O(μ_{k}^{1/2})

⇒ |Δx_{i}| = (\frac{x_{i}}{y_{i}} )^{1/2} O(μ_{k}^{1/2}) = O(μ_{k}^{1/2}),

as required. The proof for Δy_{J} is identical.

For the remaining inequality (39), we use the result of Lemma A.3. Since the feasible set for (65) is nonempty, there is a feasible vector pair (w̄, z̄) such that

∥(w̄, z̄)∥ = O(∥r∥) + O(∥Δx_{J}∥ + ∥Δx_{N}∥ + ∥Δy_{J}∥ + ∥Δy_{B}∥)
= O(μ_{k}^{1/2}),

where the last equality follows from (6), (37), and (38). The remainder of the proof follows by using identical logic to that of Lemma 5.3 of Wright [7], so we omit the details.

References

[1] A. J. Hoffman, On approximate solutions of systems of linear inequalities, J. Res. Nat. Bur. Standards, 49 (1952), pp. 263–265.

[2] J. Ji, F. Potra, and S. Huang, A predictor-corrector method for linear complementarity problems with polynomial complexity and superlinear convergence, Tech. Rep. 18, Department of Mathematics, University of Iowa, Iowa City, Iowa, August 1991.

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

[4] S. Mizuno, Polynomiality of Kojima-Meggido-Mizuno infeasible-interior-point algorithm for linear programming, Tech. Rep. 1006, School of Operations Research and Industrial Engineering, Cornell University, 1992.

[5] S. Mizuno, M. Todd, and Y. Ye, On adaptive-step primal-dual interior point algorithms for linear programming, Tech. Rep. 944, School of Operations Research and Industrial Engineering, Cornell University, Ithaca, N.Y., 1990. To appear in Mathematics of Operations Research.

[6] F. A. Potra, An infeasible interior-point predictor-corrector algorithm for linear programming, Tech. Rep. 26, Department of Mathematics, University of Iowa, Iowa City, Iowa, June 1992.

[7] 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, Illinois, October 1992.

[8] , A path-following infeasible-interior-point algorithm for linear complementarity problems, Preprint MCS–P334–1192, Mathematics and Computer Science Division, Argonne National Laboratory, Argonne, Illinois, November 1992. To appear in Optimization Methods and Software.

[9] Y. Ye, On the finite convergence of interior-point algorithms for linear programming, Tech. Rep. 91–5, Department of Management Sciences, University of Iowa, Iowa City, February 1991.

[10] Y. Ye and K. Anstreicher, On quadratic and O(\sqrt{n}L) convergence of a predictorcorrector algorithm for LCP, Tech. Rep. 91–20, Department of Management Sciences, University of Iowa, Iowa City, Iowa, November 1991.

[11] Y. Zhang, On the convergence of an infeasible interior-point algorithm for linear programming and other problems, Tech. Rep. 92–07, Department of Mathematics and Statistics, University of Maryland, Baltimore County, Maryland, April 1992.