A Superquadratic Infeasible-Interior-Point Method for Linear Complementarity Problems

Stephen Wright^{∗}and Yin Zhang^{†}

February 8, 1994 (Revised December, 1994)

RESEARCH REPORT YZ94-03, DEPARTMENT OF MATHEMATICS AND STATISTICS, UMBC PREPRINT MCS-P418-0294, MCS DIVISION, ARGONNE NATIONAL LABORATORY

Abstract

We consider a modification of a path-following infeasible-interior-point algorithm described by Wright. In the new algorithm, we attempt to improve each major iterate by reusing the coefficient matrix factors from the latest step. We show that the modified algorithm has similar theoretical global convergence properties to those of the earlier algorithm, while its asymptotic convergence rate can be made superquadratic by an appropriate parameter choice.

1 Introduction

We describe an algorithm for solving the monotone linear complementarity problem (LCP), in which we aim to find a vector pair (x, y) with

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

where q ∈ IR^{n} and M is an n × n positive semidefinite matrix. The solution set to (1) is denoted by S, while the set S^{c} of strictly complementary solutions is defined as

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

Our algorithm can be viewed as a modified form of Newton’s method applied to the 2n × 2n system y = Mx + q, x_{i}y_{i} = 0, i = 1, 2, ⋯, n, in which all the iterates (x^{k}, y^{k}) are constrained to be strictly positive. For feasible starting points, this primal-dual interior-point approach was first proposed by Kojima, Mizuno, and Yoshise [4]. Superlinearly convergent methods of this type have been described by a number of authors, including McShane [7], Ji, Potra, and Huang [1], Ji et al. [2], and Ye and Anstreicher [19]. The paper of Ye and Anstreicher is particularly interesting because several asymptotic bounds on the steps proved there are used in later works. It also gives a detailed development of superlinearly convergent methods to that time.

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

†Department of Mathematics and Statistics, University of Maryland Baltimore County, Baltimore, Maryland 21228. The work of this author was based on research supported in part by the U.S. Department of Energy under Grant DE-FG02-93ER25171.

More recently, primal-dual algorithms that do not require feasible starting points have been the focus of active research. Zhang [21] described an algorithm with polynomial complexity. In later work, Potra [14], Potra and Sheng [15], Wright [16, 17, 18], and Monteiro and Wright [11] have described superlinearly convergent infeasible-interior-point methods. The algorithm we describe in this paper also falls into this class. It is an extension of the algorithm of Wright [18], which is in turn based on the work of Zhang [21] and Wright [16]. As in [18], the algorithm extends immediately to mixed monotone LCP with few complications.

The algorithm of this paper can be motivated by considering the following locally convergent algorithm for solving the system of nonlinear equations F (z) = 0, where F : IR^{N} → IR^{N} is continuously differentiable.

Choose τ ∈ (0, 1), I ≥ 0, z^{0} ∈ IR^{N}; set k ← 0;

loop:

compute d^{k} = −∇F (z^{k})^{−1}F (z^{k}); z ← z^{k} + d^{k}; for i = 0, 1, ⋯, I (improvement loop) compute d = −∇F (z^{k})^{−1}F (z);

if ∥F (z + d)∥ ≤ τ∥F (z)∥
then z ← z + d
else z^{k+1} ← z; k ← k + 1; go to loop;

end for

z^{k+1} ← z; k ← k + 1; go to loop.

On each iteration, this method takes a single Newton step and follows it up with a number of Newton-like steps calculated with the Jacobian ∇F (z^{k}). Simple analysis shows that if z^{∗} is an isolated solution to the system F (z) = 0 with ∇F (z^{∗}) nonsingular, and if ∥z^{0} − z^{∗}∥ is small enough, then {z^{k}} converges to z^{∗}. Moreover, the inner loop (with iteration index i) eventually executes for all I iterations before control passes back to the main loop and, assuming that ∇F (z) is Lipschitz continuous at z^{∗}, the convergence has Q-order I + 2 (see, for example, [12]). Note that for each value of k, the Jacobian ∇F (z^{k}) is evaluated and factored only once and, in many contexts, the steps d calculated in the improvement loop are not expensive to compute.

Our algorithm, which we describe in Section 2, is identical to that of [18] in that it takes steps of two types—safe steps, which ensure global convergence, and fast steps, which ensure fast local convergence. As in the model algorithm above, each step is followed by an attempt to improve the new iterate without recomputing and refactoring the main coefficient matrix. The inner loop terminates when it fails to make significant progress or after I iterations, whichever comes first. Additional complications arise because of the need to keep the iterates strictly positive and in a wide neighborhood of the central path. These requirements necessitate a certain amount of quite technical analysis, which we present in Sections 4 and 5. The global convergence and complexity analysis is identical to that of the algorithm in [18], which differs from the present algorithm only in the lack of an improve phase. We state the relevant results, omitting most of the details, in Section 3. Some preliminary numerical results appear in Section 6.

The idea of reusing the matrix factors first appeared in Karmarkar, Lagarias, Slutsman and Wang [3]. Mehrotra presented a practical implementation of a predictor-corrector algorithm that reuses the matrix factors in [9] and an asymptotic theoretical analysis in [8]. Mehrotra’s theoretical algorithm differs significantly from the practical method. It requires strict feasibility of all iterates, and it uses the Mizuno-Todd-Ye [10] predictor-corrector framework, which confines the iterates to a narrow Euclidean-norm neighborhood of the central path and requires corrector steps (and hence extra matrix factorizations) to be performed regularly. Our algorithm achieves rapid asymptotic convergence like that described in [8] but for a wider class of problems (LCP). Moreover, our algorithm is closer to computational practice in its use of infeasible iterates and a much wider neighborhood of the central path.

In the remainder of the paper, we use IR_{+}^{n} to denote the nonnegative orthant in IR^{n}. Subscripts on matrices and vectors indicate components, while superscripts on matrices and vectors and subscripts on scalars denote iteration numbers (usually k).

2 The Algorithm

To describe the step between successive iterates, we define for any vector pair (x, y) ∈ IR_{+}^{n} × IR_{+}^{n} the following quantities:

μ = x^{T} y/n,  r = y − Mx − q,  e = (1, 1, ⋯, 1)^{T} ,

and, for any vector x ∈ IR_{+}^{n},

X = diag(x) = diag(x_{1}, x_{2}, ⋯, x_{n}).

When (x, y) = (x^{k}, y^{k}) (that is, the k-th iterate of the algorithm), we use r^{k}, μ_{k}, and X^{k} to denote r, μ, and X, respectively.

During the k-th iteration of the main loop, each search direction (u, v) and step length α̃ is calculated as follows.

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

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

Set

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

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

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

The inequality (4b) ensures that the component-wise products x_{j}y_{j} approach zero at approximately the same rate. They stay in a wide neighborhood of the central path, where x_{j}y_{j} = μ for all j = 1, ⋯, n — hence the term “path-following.” The inequality (4a) ensures that when the current point is infeasible, the decrease in infeasibility ∥r∥ on the current step is at least as great as the decrease in the complementarity gap μ, modulo a factor of (1 − β̃).

The basic form of our algorithm, given below, is the same as the one described in Wright [18] except for the addition of the improve procedure.

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

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

while μ_{k} > 0 solve (2)–(4) with (x, y) = (x^{k}, y^{k}), σ̃ = 0, β̃ = γ̄^{t_{k}}, γ̃ = γ_{min} + γ̄^{t_{k}}(γ_{max} − γ_{min}); if (x^{k} + α̃u)^{T} (y^{k} + α̃v)/n ≤ ρμ_{k} then β_{k} ← β̃, t ← t_{k} + 1, γ ← γ̃; else solve (2)–(4) with (x, y) = (x^{k}, y^{k}), σ̃ ∈ [σ̄, \frac{1}{2}], β̃ = 0, γ̃ = γ_{k};

β_{k} ← 0, t ← t_{k}, γ ← γ_{k};

end if

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

improve ((x, y), t, ν, γ, (x^{k}, y^{k}));

t_{k+1} ← t, ν_{k+1} ← ν, γ_{k+1} ← γ, (x^{k+1}, y^{k+1}) ← (x, y), k ← k + 1; end while.

We refer to the steps that are computed with σ̃ = 0 as fast steps because they lead to rapid local convergence, while the steps with σ̃ ∈ [σ̄, \frac{1}{2}] are safe steps because they ensure global convergence.

The improve procedure, which reuses the coefficient matrix in (2) to improve the new iterate, takes a combination of safe and fast steps, just like the main algorithm. The main difference is that the procedure is terminated if an improvement in μ of at least a factor of τ ∈ (ρ, 1) is not achieved. The user supplies the parameter τ and the nonnegative integer I, where I is the maximum number of steps that can be taken in improve.

improve ((x, y), t, ν, γ, (x^{k}, y^{k}))

Given τ ∈ (ρ, 1), I ≥ 0,

for i = 1, 2, ⋯, I if μ = 0 then return; solve (2)–(4) with σ̃ = 0, β̃ = γ̄^{t}, γ̃ = γ_{min} + γ̄^{t}(γ_{max} − γ_{min}); if (x + α̃u)^{T} (y + α̃v)/n ≤ ρμ then t ← t + 1, γ ← γ̃; else solve (2)–(4) with σ̃ ∈ [σ̄, \frac{1}{2}], β̃ = 0, γ̃ = γ;

if (x + α̃u)^{T} (y + α̃v)/n > τμ then return;

end if

ν ← ν(1 − α̃), (x, y) ← (x, y) + α̃(u, v); end for.

In the special case I = 0, improve is vacuous and the algorithm reduces to the method of [18]. We refer the reader to that paper for the intuitive motivation behind the use of safe and fast steps.

Some of the fundamental properties of the iteration sequence (x^{k}, y^{k}) are not affected by the inclusion of the improve procedure. We still have

r^{k} = ν_{k}r^{0}  (5)

and also the following result, which is similar to Lemma 3.1 of [17].

Lemma 2.1 Suppose that the initial point is infeasible, that is, r^{0} ≠ 0. Then the positive constant β̂ defined by

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

is such that

μ_{k} ≥ β̂ν_{k}μ_{0} = β̂\frac{∥r^{k}∥}{∥r^{0}∥}μ_{0}, ∀k ≥ 0.

We also have the following result, which shows that the algorithm either terminates finitely at a solution of (1) or else generates an infinite sequence {(x^{k}, y^{k})} of strictly positive iterates. The proof is a simple modification of [18, Lemma 3.2] and is omitted.

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

We assume throughout the remainder of the paper that finite termination does not occur, that is, all iterates (x^{k}, y^{k}) and all the intermediate points (x, y) generated in the improve procedure are strictly positive.

3 Global Convergence

The analysis of global convergence and polynomial complexity is nearly identical to that of [18, Section 3]. We need only note that (5) still applies and that all iterates (x^{k}, y^{k}) satisfy x_{j}^{k}y_{j}^{k} ≥ γ_{min}μ_{k}, j = 1, ⋯, n. The intermediate points generated by improve have the same properties. The technical results from [18, Section 3] can therefore be applied to show that nontrivial progress is made at each safe step. The presence of improve and the fast steps cannot hinder (and very often speed) the convergence.

In this section we summarize the main results from [18, Section 3] and state the sole assumption required for global convergence, which is as follows.

Assumption 1 S ≠ ∅.

Theorem 3.1 If a safe step is taken at iteration k, then there is a constant ω > 0 such that the step length α_{k} has

α_{k} ≥ \frac{1}{ω} .

If the initial point (x^{0}, y^{0}) is chosen as

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

where

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

for some (x^{∗}, y^{∗}) ∈ S, then ω = O(n^{2}).

Proof. See [18, Lemma 3.4, Theorem 3.5], where a different definition of ω is used.

The main global convergence result is as follows.

Theorem 3.2 The complementarity gap μ_{k} converges geometrically to zero.

Proof. As in Wright [18, Theorem 3.6], we can show that if a safe step is taken at iteration k, we have

(x^{k} + α_{k}u)^{T} (y^{k} + α_{k}v)/n ≤ (1 − \frac{1}{4ω} ) μ_{k},

while if a fast step is taken, we have

(x^{k} + α_{k}u)^{T} (y^{k} + α_{k}v)/n ≤ ρμ_{k}.

Since the complementarity gap may be decreased further by improve, we have μ_{k+1} ≤ (x^{k} + α_{k}u)^{T} (y^{k} + α_{k}v)/n and therefore

μ_{k+1} ≤ max (1 − \frac{1}{4ω} , ρ) μ_{k},

from which the result follows.

Finally, we state the polynomial complexity result.

Theorem 3.3 [18, Corollary 3.7] Let ϵ > 0 be given. Suppose that the starting point is defined by (6), (7), 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_{ϵ}.

4 Technical Results

In the remainder of the paper, we turn our attention to the latter stages of the algorithm. We show that the algorithm eventually takes only fast steps (that is, the then branch of the main conditional statement is executed). Moreover, the improve routine eventually takes fast steps on all I of its iterations, so that a total of I + 1 fast steps is taken for each factorization of the coefficient matrix in (2).

In this section, we prove some results about the steps generated in this fast phase of the algorithm. In particular, we look at the effects of the inexact coefficient matrix in (2) on the steps calculated within improve.

We start by defining the two assumptions for the local convergence analysis, which will be implicitly assumed to hold throughout the remainder of the paper.

Assumption 2 S^{c} ≠ ∅.

Assumption 3 S is bounded.

For monotone LCP, a sufficient condition for Assumption 3 is the existence of a strictly feasible pair (x̄, ȳ) such that ȳ = Mx̄ + q, (x̄, ȳ) > 0. This can be seen from the fact that for any (x^{∗}, y^{∗}) ∈ S

(x^{∗} − x̄)^{T} (y^{∗} − ȳ) = (x^{∗} − x̄)^{T} M(x^{∗} − x̄) ≥ 0,

implying

x̄^{T} y^{∗} + ȳ^{T} x^{∗} ≤ x̄^{T} ȳ.

By choosing any particular strictly complementary solution (x^{∗}, y^{∗}), we can define index sets B and N by

B = {j | x_{j}^{∗} > 0}, N = {j | y_{j}^{∗} > 0}.

It is well known that the global convergence of the algorithm guarantees that the iteration sequence {(x^{k}, y^{k})} approaches the solution set S (see the error bound result of Mangasarian [6], for example). Therefore, Assumption 3 implies the boundedness of the iteration sequence {(x^{k}, y^{k})}, as given in the following lemma.

Lemma 4.1 There is a positive constant C_{3} such that ∥(x^{k}, y^{k})∥ ≤ C_{3} for all k ≥ 0.

The next two results are simple modifications of results from Wright [18, Section 4]. Since we will apply these results to intermediate points generated by improve as well as to the main iterates (x^{k}, y^{k}), we state them in a more general form than in [18]. The proofs are, however, not affected. Boundedness of the iteration sequence is not necessary for either result, and neither is Assumption 3.

Lemma 4.2 ([18, Lemma 4.1]) Let (x, y) ≥ 0 be such that

r = y − Mx − q = νr^{0} for some ν ∈ [0, \frac{1}{2}],

and μ = x^{T} y/n ≥ β̂νμ_{0} for this value of ν. Then for some constant C_{4} > 0 we have

∥x_{N}∥ ≤ C_{4}μ,  ∥y_{B}∥ ≤ C_{4}μ.  (8)

Lemma 4.3 Let (x, y) be any point with the properties defined in Lemma 4.2, and suppose in addition that x_{j}y_{j} ≥ γ_{min}μ. Let (ū, v̄) be the search direction obtained by solving

[ \begin{matrix}M & −I \\ Y & X \end{matrix} ] [ ū ; v̄ ] = [ \begin{matrix}r \\ −XY e + σ̃μe \end{matrix} ] ,  (9)

where σ̃ ∈ [0, 1). Then there exists a positive constant C_{5} such that

∥ū_{N}∥ ≤ C_{5}μ,  ∥v̄_{B}∥ ≤ C_{5}μ.  (10)

If in addition σ̃ = 0, there is a constant C_{6} > 0 such that

∥ū_{B}∥ ≤ C_{6}μ,  ∥v̄_{N}∥ ≤ C_{6}μ.  (11)

Proof. Follows from Lemma 4.2 and Theorem 4.5 of [18].

We now turn to the “approximate” fast steps computed by (2), where (x, y) is either the current iterate (x^{k}, y^{k}) or some intermediate point generated in the call to improve at iteration k. It is obvious from the algorithm definition that we have

μ = x^{T} y/n ≤ μ_{k}.  (12)

We also assume that the point (x, y) is not too far from (x^{k}, y^{k}) in the sense that there is a constant χ ≥ 1 independent of k such that

∥(x^{k} − x, y^{k} − y)∥ ≤ χμ_{k}.  (13)

Later, in Theorem 5.1, we choose a particular value of χ that ensures that (13) is eventually satisfied by all intermediate points generated by the procedure improve. Hence, the remaining results in this section apply to all iterates visited during the improve phase.

The following result describes some characteristics of the actual search direction (u, v) calculated from (2), partly in terms of the exact search direction (ū, v̄) that satisfies (9).

Lemma 4.4 Let (x, y) be a vector pair satisfying the assumptions of Lemma 4.3 and, in addition, the properties (12) and (13). Then if σ̃ = 0, there are positive constants C_{7}, C_{8}, and C_{9} independent of k and χ such that the following bounds are satisfied:

∥u − ū∥ ≤ C_{7}χμ,  ∥v − v̄∥ ≤ C_{7}χμ,  (14)
∥(u, v)∥ ≤ C_{8}χμ,  (15)
∥ū_{N} − u_{N}∥ ≤ C_{9}χμμ_{k},  ∥v̄_{B} − v_{B}∥ ≤ C_{9}χμμ_{k}.  (16)

Proof. From (2) and (9), respectively, we have that

[ \begin{matrix}M & −I \\ Y^{k} & X^{k} \end{matrix} ] [ \begin{matrix}u \\ v \end{matrix} ] = [ \begin{matrix}r \\ −XY e \end{matrix} ] ,  [ \begin{matrix}M & −I \\ Y & X \end{matrix} ] [ ū ; v̄ ] = [ \begin{matrix}r \\ −XY e \end{matrix} ] .  (17)

By combining these two equations, we obtain

[ \begin{matrix}M & −I \\ Y^{k} & X^{k} \end{matrix} ] [ \begin{matrix}ū − u \\ v̄ − v \end{matrix} ] = [ \begin{matrix}0 \\ (Y^{k} − Y)ū + (X^{k} − X)v̄ \end{matrix} ] .  (18)

Now from (13) and by applying Lemma 4.3 to the second equation in (17), there is a constant C̄_{7} independent of k and χ such that

∥(Y^{k} − Y)ū + (X^{k} − X)v̄∥ ≤ C̄_{7}χμμ_{k}.  (19)

Defining

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

and multiplying the lower block in the system (18) by (X^{k}Y^{k})^{−1/2}, we obtain

D^{k}(ū − u) + (D^{k})^{−1}(v̄ − v) = (X^{k}Y^{k})^{−1/2}[(Y^{k} − Y)ū + (X^{k} − X)v̄].  (20)

Using the upper block of (18), we have (v̄ − v) = M(ū − u), and so it follows from positive semidefiniteness of M that

(ū − u)^{T} (v̄ − v) ≥ 0.  (21)

By taking the Euclidean norm of both sides of (20), and using (21), we have

∥D^{k}(ū − u)∥^{2} + ∥(D^{k})^{−1}(v̄ − v)∥^{2} ≤ ∥(X^{k}Y^{k})^{−1/2}∥^{2}∥(Y^{k} − Y)ū + (X^{k} − X)v̄∥^{2},

Therefore

∥D^{k}(ū − u)∥ ≤ ∥(X^{k}Y^{k})^{−1/2}∥∥(Y^{k} − Y)ū + (X^{k} − X)v̄∥,
∥(D^{k})^{−1}(v̄ − v)∥ ≤ ∥(X^{k}Y^{k})^{−1/2}∥∥(Y^{k} − Y)ū + (X^{k} − X)v̄∥.

Now since x_{j}^{k}y_{j}^{k} ≥ γ_{min}μ_{k}, we have

∥(X^{k}Y^{k})^{−1/2}∥ = max_{j=1,···,n} (x_{j}^{k}y_{j}^{k})^{−1/2} ≤ γ_{min}^{−1/2}μ_{k}^{−1/2}.

Therefore from (19) we have

∥D^{k}(ū − u)∥ ≤ C̄_{7}χγ_{min}^{−1/2}μμ_{k}^{1/2}.

Taking any j = 1, ⋯, n, we find that

|\frac{(y_{j}^{k})^{1/2}}{(x_{j}^{k})^{1/2}} (ū_{j} − u_{j})| ≤ ∥D^{k}(ū − u)∥ ≤ C̄_{7}χγ_{min}^{−1/2}μμ_{k}^{1/2}.

Hence,

|ū_{j} − u_{j}| ≤ max_{j=1,···,n} \frac{(x_{j}^{k})^{1/2}}{(y_{j}^{k})^{1/2}} C̄_{7}χγ_{min}^{−1/2}μμ_{k}^{1/2}
≤ max_{j=1,···,n} \frac{x_{j}^{k}}{(x_{j}^{k}y_{j}^{k})^{1/2}} C̄_{7}χγ_{min}^{−1/2}μμ_{k}^{1/2}
≤ \frac{C_{3}}{γ_{min}^{1/2}μ_{k}^{1/2}} C̄_{7}χγ_{min}^{−1/2}μμ_{k}^{1/2}
≤ \frac{C_{7}}{\sqrt{n}}χμ,

for C_{7} defined in an obvious way. We have proved the first inequality in (14); the proof of the second inequality is similar.

For (16), we repeat the logic above to obtain for i ∈ N that

|ū_{j} − u_{j}| ≤ \frac{(x_{j}^{k})^{1/2}}{(y_{j}^{k})^{1/2}} C̄_{7}χγ_{min}^{−1/2}μμ_{k}^{1/2}
≤ \frac{x_{j}^{k}}{(x_{j}^{k}y_{j}^{k})^{1/2}} C̄_{7}χγ_{min}^{−1/2}μμ_{k}^{1/2}
≤ \frac{C_{4}μ_{k}}{γ_{min}^{1/2}μ_{k}^{1/2}} C̄_{7}χγ_{min}^{−1/2}μμ_{k}^{1/2}
≤ \frac{C_{9}}{\sqrt{n}}χμμ_{k},

where C_{9} is defined appropriately. The bound for ∥v̄_{B} − v_{B}∥ follows similarly. To prove (15), we have from Lemma 4.3 and (14) that

∥(u, v)∥ ≤ ∥(ū, v̄)∥ + ∥(u − ū, v − v̄)∥ ≤ 2(C_{5} + C_{6})μ + 2C_{7}χμ ≤ C_{8}χμ,

where we have defined C_{8} = 2(C_{5} + C_{6} + C_{7}) and used the assumption that χ ≥ 1.

We now state the main result of this section, in which we obtain an estimate for the step length α̃ along a (possibly approximate) fast step direction (u, v). The point (x, y) considered in this theorem represents either the main iterate (x^{k}, y^{k}) itself or one of the intermediate points generated by improve during iteration k. The following positive constants, all of which are independent of k and χ, are used in this result.

C̄_{10} = C_{9}(C_{3} + C_{6}) + C_{7}(C_{4} + C_{5})
C_{10} = 2(C_{5}C_{6} + C̄_{10} + C_{7}C_{9})
C_{12} = \frac{2C_{10}}{(1 − γ̄)(γ_{max} − γ_{min})}  (22)
C_{13} = 2(C_{3}C_{9} + C_{4}C_{7} + C_{8}^{2}),
C_{14} = C_{12} + C_{13}/n.

Theorem 4.5 Let (x, y) be a point that satisfies the assumptions of Lemma 4.4, and in addition

μ ≤ min (1, \frac{n}{C_{13}χ^{2}} ) .  (23)

Let t be a positive integer such that for γ defined by

γ = γ_{min} + γ̄^{t−1}(γ_{max} − γ_{min})

we have x_{j}y_{j} ≥ γμ for j = 1, ⋯, n, and suppose for this value of t that

C_{14}χ^{2} \frac{μ_{k}}{γ̄^{t}} ≤ ρ.  (24)

Then if a fast step is attempted from the point (x, y) with

σ̃ = 0,  β̃ = γ̄^{t},  γ̃ = γ_{min} + γ̄^{t}(γ_{max} − γ_{min}),

and the search direction (u, v) is calculated from (2), the resulting step length α̃ obtained from (2), (3), and (4) satisfies

α̃ ≥ 1 − C_{12}χ^{2} \frac{μ_{k}}{γ̄^{t}} .

Moreover, the fast step is accepted with

(x + α̃u)^{T} (y + α̃v)/n ≤ C_{14}χ^{2} \frac{μ_{k}}{γ̄^{t}} μ.  (25)

Proof. The proof is in three stages. First, we show that the tests (4) are satisfied for all α in the range

[0, 1 − C_{12}χ^{2} \frac{μ_{k}}{γ̄^{t}} ] .  (26)

Second, we show that μ(α) defined by (3) is decreasing on the interval α ∈ [0, 1]. Third, we show that

μ(α̃) ≤ C_{14}χ^{2} \frac{μ_{k}}{γ̄^{t}} μ ≤ ρμ,  (27)

which proves the result.

We first consider the condition (4a). From the left-hand side, we obtain

(x + αu)^{T} (y + αv)
= (x + αū + α(u − ū))^{T} (y + αv̄ + α(v − v̄))
= x^{T} y(1 − α) + α^{2}ū^{T} v̄ + α(x + αū)^{T} (v − v̄) + α(y + αv̄)^{T} (u − ū)  (28)
+α^{2}(u − ū)^{T} (v − v̄).

Now, using Lemma 4.1 and the inequalities (8), (10), (11), (14), (15), and (16), we have

|ū^{T} v̄| ≤ 2C_{5}C_{6}μ^{2} ≤ 2C_{5}C_{6}μμ_{k}

and

|(u − ū)^{T} (y + αv̄)| ≤ ∥u_{N} − ū_{N}∥(∥y_{N}∥ + ∥v̄_{N}∥) + ∥u_{B} − ū_{B}∥(∥y_{B}∥ + ∥v̄_{B}∥)
≤ C_{9}χμμ_{k}(C_{3} + C_{6}μ) + C_{7}χμ^{2}(C_{4} + C_{5})
≤ C̄_{10}χμμ_{k},

where we have used μ ≤ 1 and μ ≤ μ_{k} to derive the last inequality. Similarly, we have

|(v − v̄)^{T} (x + αū)| ≤ C̄_{10}χμμ_{k},

while for the remaining term in (28) we have from Lemma 4.4 that

|(u − ū)^{T} (v − v̄)| ≤ 2C_{7}C_{9}χ^{2}μ^{2}μ_{k} ≤ 2C_{7}C_{9}χ^{2}μμ_{k}.

Hence, since χ ≥ 1, we have from the definition of C_{10} that

|(x + αu)^{T} (y + αv) − (1 − α)x^{T} y| ≤ C_{10}χ^{2}μμ_{k}.  (29)

Since β̃ = γ̄^{t}, we have that (4a) is satisfied provided that

C_{10}χ^{2}μμ_{k} ≤ (1 − α)γ̄^{t}nμ,

which is certainly true provided that

1 − α ≥ \frac{C_{10}χ^{2}μ_{k}}{nγ̄^{t}} .

From the definition of C_{12}, since 1 − γ̄ and γ_{max} − γ_{min} both lie in the range (0, 1), we have

\frac{C_{10}χ^{2}μ_{k}}{nγ̄^{t}} ≤ \frac{C_{12}χ^{2}μ_{k}}{γ̄^{t}} ,

so the inequality (4a) certainly holds for all α in the range (26).

Turning to the second inequality (4b), we have by an argument similar to the one above that

(x_{j} + αu_{j})(y_{j} + αv_{j}) ≥ x_{j}y_{j}(1 − α) − C_{10}χ^{2}μμ_{k} ≥ γμ(1 − α) − C_{10}χ^{2}μμ_{k},

while from (29), we have

(x + αu)^{T} (y + αv)/n ≤ (1 − α)μ + C_{10}χ^{2}μμ_{k}/n.

Hence, the inequality (4b) holds provided that

γ̃μ(1 − α) + C_{10}χ^{2}μμ_{k}(γ̃/n) ≤ γμ(1 − α) − C_{10}χ^{2}μμ_{k},

which is certainly true whenever the inequality

(γ − γ̃)(1 − α) ≥ 2C_{10}χ^{2}μ_{k}  (30)

holds. Since

γ − γ̃ = [γ_{min} + γ̄^{t−1}(γ_{max} − γ_{min})] − [γ_{min} + γ̄^{t}(γ_{max} − γ_{min})] = γ̄^{t−1}(1 − γ̄)(γ_{max} − γ_{min}),

we find that (30) holds whenever

1 − α ≥ \frac{2C_{10}χ^{2}μ_{k}}{γ̄^{t−1}(1 − γ̄)(γ_{max} − γ_{min})},

which, by definition of C_{12}, is true for α in the range (26).

For the second part of the proof, we show that μ(α) defined by (3) is decreasing on the range α ∈ [0, 1]. Taking the derivative, we have

nμ′(α) = (x^{T} v + y^{T} u) + 2αu^{T} v
= (x^{T} v̄ + y^{T} ū) + x^{T} (v − v̄) + y^{T} (u − ū) + 2αu^{T} v
= −x^{T} y + x^{T} (v − v̄) + y^{T} (u − ū) + 2αu^{T} v.  (31)

However, we can use Lemma 4.1 and relations (8), (14), and (16) to obtain

|x^{T} (v − v̄)| ≤ ∥x_{B}∥∥v_{B} − v̄_{B}∥ + ∥x_{N}∥∥v_{N} − v̄_{N}∥
≤ C_{3}C_{9}χμμ_{k} + C_{4}C_{7}χμ^{2}
≤ (C_{3}C_{9} + C_{4}C_{7})χμμ_{k},  (32)

where we have used μ ≤ μ_{k} in the last inequality. A similar bound can be obtained for |y^{T} (u − ū)|. For the final term in (31), we have

|u^{T} v| ≤ C_{8}^{2}χ^{2}μ^{2} ≤ C_{8}^{2}χ^{2}μμ_{k}.  (33)

Substituting these relations in (31) and using the definition of C_{13}, we obtain

nμ′(α) ≤ [−n + C_{13}χ^{2}μ_{k}]μ.

It follows from (23) that the term in brackets is negative, and hence μ′(α) ≤ 0 for all α ∈ [0, 1].

Finally, we observe that the step length α̃ actually selected by the procedure will be at least as long as the upper bound of (26). Hence, using (32), (33), and the definitions of C_{13} and C_{14}, we have

(x + α̃u)^{T} (y + α̃v)
≤ x^{T} y(1 − α̃) + |x^{T} (v − v̄)| + |y^{T} (u − ū)| + |u^{T} v|
≤ (x^{T} y)C_{12}χ^{2} \frac{μ_{k}}{γ̄^{t}} + 2(C_{3}C_{9} + C_{4}C_{7})χμμ_{k} + C_{8}^{2}χ^{2}μμ_{k}
≤ (x^{T} y)\frac{μ_{k}}{γ̄^{t}} [C_{12} + C_{13}/n]χ^{2}
= (x^{T} y)\frac{μ_{k}}{γ̄^{t}} C_{14}χ^{2}.

Therefore (25) holds. Acceptance of this step follows from (24), since we have (x + α̃u)^{T} (y + α̃v)/n ≤ ρμ.

We close this section with a result that is important in defining the onset of the algorithm’s fast phase.

Lemma 4.6 There is a constant η < 1 such that

\frac{μ_{k+1}}{γ̄^{t_{k+1}}} ≤ η \frac{μ_{k}}{γ̄^{t_{k}}} ,  ∀k ≥ 0.  (34)

Proof. When the safe branch of the main algorithm is taken at iteration k, we have from the proof of Theorem 3.2 that

(x^{k} + α̃u)^{T} (y^{k} + α̃v)/n ≤ (1 − \frac{1}{4ω} ) μ_{k},

while the value of t is unaltered. In the subsequent call to improve, the value of t might be incremented. Whenever this happens, we are guaranteed that the complementarity gap μ decreases by a factor of at least ρ, so the ratio μ/γ̄^{t} will also decrease by a factor of at least ρ/γ̄ < 1. Hence, when the safe branch is taken, we have

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

When the fast branch is taken, we have t ← t + 1 and

(x^{k} + α̃u)^{T} (y^{k} + α̃v)/n ≤ ρμ_{k},

so the ratio μ/γ̄^{t} decreases by a factor of at least ρ/γ̄. The comments above ensure that the subsequent call to improve can only accentuate this decrease, so in this case we have

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

The result is obtained by defining

η = max (1 − \frac{1}{4ω} , \frac{ρ}{γ̄} ) .

5 Local Convergence

In this section, we state and prove our two main local convergence results. First, we define a threshold value of μ_{k}/γ̄^{t_{k}} below which both the main algorithm and the procedure improve take only fast steps. Second, we show that the resulting superlinear convergence has Q-order I + 2.

Theorem 5.1 Define

χ = 2(C_{5} + C_{6}) exp (\frac{C_{8}ρ}{1 − ρ}) ,  (35)

and let K_{1} be the smallest value of k for which the following three conditions are satisfied:

ν_{k} ≤ 1/2,  C_{14}χ^{2} \frac{μ_{k}}{γ̄^{t_{k}+I}} ≤ ρ,  μ_{k} ≤ min (1, \frac{n}{C_{13}χ^{2}} ) .  (36)

Then the fast branch is taken in the main algorithm and, moreover, I fast steps are taken in the call to improve.

Proof. Existence of K_{1} is guaranteed by Lemma 4.6. We choose any k ≥ K_{1}. Our proof proceeds by showing first that the step taken from (x^{k}, y^{k}) in the main algorithm is a fast step. We then prove by induction that I fast steps are taken inside the procedure improve. Our main tool in both cases is Theorem 4.5.

For the first part of the proof, we apply Theorem 4.5 with

(x, y) = (x^{k}, y^{k}),  t = t_{k},  γ = γ_{k}.  (37)

Note that the point (x, y) satisfies the assumptions of Lemmas 4.2 and 4.3 (by definition of K_{1}, r^{k}, ν_{k}, etc.) and the conditions (12) and (13) (trivially). Clearly also x_{j}^{k}y_{j}^{k} ≥ γ_{k}μ_{k} for all j = 1, ⋯, n, and the condition (23) also holds. Because

C_{14}χ^{2} \frac{μ_{k}}{γ̄^{t}} = C_{14}χ^{2} \frac{μ_{k}}{γ̄^{t_{k}}} ≤ C_{14}χ^{2} \frac{μ_{K_{1}}}{γ̄^{t_{K_{1}}}} ≤ C_{14}χ^{2} \frac{μ_{K_{1}}}{γ̄^{t_{K_{1}}+I}} ≤ ρ,

the condition (24) also holds. Hence the conditions of Theorem 4.5 are satisfied by the choices (37), and therefore a fast step is taken by the main algorithm.

We turn now to the procedure improve. Our aim is to show inductively that if (x, y) is the current vector pair at the commencement of the i-th iteration of this procedure, then

∥(x^{k}, y^{k}) − (x, y)∥ ≤ [2(C_{5} + C_{6}) ∏_{l=1}^{i−1}(1 + C_{8}ρ^{l})] μ_{k}.  (38)

Moreover, we show that a fast step is taken from this vector (x, y) during the i-th iteration of improve. Note for future reference that

log ∏_{l=1}^{i−1}(1 + C_{8}ρ^{l}) = ∑_{l=1}^{i−1} log(1 + C_{8}ρ^{l}) ≤ ∑_{l=1}^{i−1} C_{8}ρ^{l} ≤ \frac{C_{8}ρ}{1 − ρ},

and so

2(C_{5} + C_{6}) ∏_{l=1}^{i−1}(1 + C_{8}ρ^{l}) ≤ χ, i = 1, ⋯, I.

Consider the case i = 1, that is, the first iteration of improve. We aim to use Theorem 4.5 again, so we start by checking that the point just generated by the main algorithm satisfies the assumptions of this theorem. In other words, the choices

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

must be shown to satisfy these assumptions. It is easy to see that the assumptions of Lemmas 4.2 and 4.3 and the condition (12) are satisfied. To see (13), note that the fast step just taken at iteration k of the main algorithm was computed with an exact coefficient matrix, that is, we have (ū, v̄) = (u, v). Hence we can apply Lemma 4.3 to deduce that

∥(x^{k}, y^{k}) − (x, y)∥ ≤ ∥(u, v)∥ = ∥(ū, v̄)∥ ≤ 2(C_{5} + C_{6})μ_{k}.

Thus the bound (38), and therefore also (13), holds for this point (x, y). The conditions (23) and x_{j}y_{j} ≥ γμ clearly hold, while (24) also holds because

C_{14}χ^{2} \frac{μ_{k}}{γ̄^{t}} = C_{14}χ^{2} \frac{μ_{k}}{γ̄^{t_{k}+1}} ≤ C_{14}χ^{2} \frac{μ_{k}}{γ̄^{t_{k}+I}} ≤ C_{14}χ^{2} \frac{μ_{K_{1}}}{γ̄^{t_{K_{1}}+I}} ≤ ρ.

Hence Theorem 4.5 applies, and we have shown that a fast step is taken on the first iteration of improve.

We now consider the general iteration i of the internal loop of improve. We assume that our assertions hold for iterations 1 through i − 1. Let (x^{−}, y^{−}) denote the value of (x, y) at the start of iteration i − 1, and let (u^{−}, v^{−}) be the search direction calculated during this iteration, while as before (x, y) is the current point at the start of iteration i. To obtain an estimate of ∥(x^{k}, y^{k}) − (x, y)∥, we note by our inductive hypothesis (38) that

∥(x^{k}, y^{k}) − (x^{−}, y^{−})∥ ≤ [2(C_{5} + C_{6}) ∏_{l=1}^{i−2}(1 + C_{8}ρ^{l})] μ_{k}.

We now apply Lemma 4.4 to the step (u^{−}, v^{−}) taken during iteration i − 1, with χ replaced by 2(C_{5} + C_{6}) ∏_{l=1}^{i−2}(1 + C_{8}ρ^{l}), to find that

∥(x^{k}, y^{k}) − (x, y)∥
≤ ∥(x^{k}, y^{k}) − (x^{−}, y^{−})∥ + ∥(x^{−}, y^{−}) − (x, y)∥
≤ 2(C_{5} + C_{6}) ∏_{l=1}^{i−2}(1 + C_{8}ρ^{l})μ_{k} + ∥(u^{−}, v^{−})∥

≤ 2(C_{5} + C_{6}) ∏_{l=1}^{i−2}(1 + C_{8}ρ^{l})μ_{k} + C_{8} [2(C_{5} + C_{6}) ∏_{l=1}^{i−2}(1 + C_{8}ρ^{l})] μ^{−}

≤ [2(C_{5} + C_{6}) ∏_{l=1}^{i−2}(1 + C_{8}ρ^{l})] (1 + C_{8}ρ^{i−1})μ_{k}.

The final inequality follows from the fact that μ^{−} ≤ ρ^{i−1}μ_{k}, since (x^{−}, y^{−}) is arrived at by taking i − 1 fast steps (one step in the main algorithm, followed by i − 2 iterations of the improve loop), at each of which a reduction factor of at least ρ is achieved. We have now shown that the bound (38) continues to hold at iteration i. It is easy to check that the remaining conditions required by Theorem 4.5 hold. We mention only (24), which holds for t = t_{k} + i because

C_{14}χ^{2} \frac{μ_{k}}{γ̄^{t}} = C_{14}χ^{2} \frac{μ_{k}}{γ̄^{t_{k}+i}} ≤ C_{14}χ^{2} \frac{μ_{k}}{γ̄^{t_{k}+I}} ≤ C_{14}χ^{2} \frac{μ_{K_{1}}}{γ̄^{t_{K_{1}}+I}} ≤ ρ.

Hence, we can apply Theorem 4.5 again to deduce that a fast step is taken at iteration i, and our result is proved.

Our final result is to show high-order convergence of the sequence {μ_{k}} to zero. We show that this convergence has a Q-order of at least I + 2, that is, for any ϵ > 0

lim sup_{k→∞} \frac{μ_{k+1}}{μ_{k}^{I+2−ϵ}} = 0.

An equivalent characterization of the Q-order I + 2 convergence is the inequality (39) below (see Potra [13]).

Theorem 5.2 The subsequence {μ_{k}}, k = 0, 1, ⋯, converges to zero with Q-order I + 2, that is,

lim inf_{k→∞} \frac{log μ_{k+1}}{log μ_{k}} ≥ I + 2.  (39)

Proof. Consider k ≥ K_{1}. Since a fast step is taken by the main algorithm and all I iterations of improve, and since Theorem 4.5 applies at all I + 1 steps, we can apply the inequality (25) I + 1 times to bound μ_{k+1} in terms of μ_{k}. The process yields

μ_{k+1} ≤ C_{14}^{I+1}χ^{2(I+1)} \frac{μ_{k}^{I+2}}{γ̄^{(I+1)t_{k}+I(I+1)/2}} .  (40)

It follows from Lemma 4.6 and (40) that

\frac{μ_{k+1}}{μ_{k}} ≤ \frac{C_{14}^{I+1}χ^{2(I+1)}}{γ̄^{I(I+1)/2}} (\frac{μ_{k}}{γ̄^{t_{k}}} )^{I+1} → 0,  (41)

that is, {μ_{k}} converges to zero at least Q-superlinearly.

By taking logarithms, we obtain from (40) that

log μ_{k+1} ≤ log (\frac{C_{14}^{I+1}χ^{2(I+1)}}{γ̄^{I(I+1)/2}} ) + (I + 2) log μ_{k} − (I + 1)t_{k} log γ̄.

We will assume that k is sufficiently large such that μ_{k} < 1. From the above,

\frac{log μ_{k+1}}{log μ_{k}} ≥ (I + 2) + log (\frac{C_{14}^{I+1}χ^{2(I+1)}}{γ̄^{I(I+1)/2}} ) / log μ_{k} − (I + 1) log γ̄ \frac{t_{k}}{log μ_{k}} .  (42)

Obviously, as k → ∞, the second term in the right-hand side vanishes. If we can show that the third term also goes to zero, then the conclusion (39) follows. Since t_{k} ≤ (I + 1)k + 1, it suffices to prove

lim_{k→∞} \frac{k}{log μ_{k}} = 0.  (43)

Suppose otherwise. Then there exist ξ ∈ (0, 1) and a subsequence {μ_{k}}_{K} ⊂ {μ_{k}} such that for all k ∈ K

\frac{k}{log μ_{k}} ≤ \frac{1}{log ξ} , or equivalently ξ^{k} ≤ μ_{k}.

From (41), there exists a positive integer J such that for all k ≥ J, μ_{k+1} ≤ \frac{ξ}{2}μ_{k}. Hence, for all k > J and k ∈ K,

ξ^{k} ≤ μ_{k} ≤ (\frac{ξ}{2})^{k−J} μ_{J}.

That is, for all k > J and k ∈ K, 2^{k} ≤ μ_{J}2^{J}/ξ^{J}. This is clearly a contradiction.

6 Numerical Examples

We include some preliminary numerical results that compare the behavior of our algorithm with the method of [18], in which improve is vacuous (I = 0).

Our test problems have M = AΛA^{T}, where A ∈ IR^{n×n} is dense with elements drawn from a uniform distribution in [−1, 1], and Λ is a diagonal matrix with diagonal elements Λ_{ii} = 10^{4ζ_{i}}, where ζ_{i} is drawn from a uniform distribution in [0, 1]. (We introduce this scaling to aviod making the problems too easy.) A solution (x^{∗}, y^{∗}) is generated so that even-numbered components of x^{∗} and odd-numbered components of y^{∗} are zero, and q is chosen so that the nonzero components of both vectors are uniformly distributed in [0, 1].

Table 1: Average performance of the algorithm for five data sets

n = 20 n = 200
factorizations  36.2  47.2
I = 0 solves  49.8  65.2
improve steps  0  0
factorizations  26.2  36.4
I = 1 solves  71.8  100.6
improve steps 20.4  21.0
factorizations  19.4  31.6
I = 3 solves  95.4  126.4
improve steps 41.4  35.4
factorizations 17.2 30.4
I = 5 solves  114.0 136.4
improve steps 53.8  43.2

The algorithmic constants have the following values:

γ_{min} = 10^{−6},  γ_{max} = 10^{−4},  σ_{min} = 10^{−4},  σ_{max} = .3,
γ̄ = .5,  ρ = min(\sqrt{σ_{max}σ_{min}}, \frac{1}{2} γ̄),  τ = .8.

Experience shows that fast steps usually do not occur until μ becomes quite small. Hence, we modify the algorithm slightly so that the fast step branch of the conditional statements in both the main algorithm and improve is not activated until μ ≤ 1. This modification does not alter any theoretical convergence properties of the algorithm.

The value of σ̃ for the safe step at major iteration k is chosen as

σ_{k} = mid(σ_{min}, μ_{k}/\sqrt{n}, σ_{max}),

where mid() denotes the median of its three arguments. For safe steps in improve, we set σ̃ to the constant σ_{max}. Termination occurs when μ_{k} ≤ 10^{−10}.

Performance of the algorithm for τ = .8 and various values of I is shown in Table 1. For each data set, we noted the number of factorizations (which equals the number of iterations of the main algorithm), the number of back solves, and the total number of steps taken in the improve phase. We averaged these figures over five sets of data to obtain the figures in Table 1.

Clearly, the improve phase decreases the number of factorizations, which is the dominant part of the cost for all reasonable problems. When M is sparse, however, the extra backsolves and other operations required to perform each improve step are not insignificant, so we should avoid taking an excessively large value of I. For our test problems, either I = 3 or I = 5 would seem to be reasonable. Our tests showed that any value of τ reasonably close to 1 gives almost identical results to our particular choice τ = .8.

7 Final Comments

We have analyzed an infeasible-interior-point algorithm that reuses matrix factors to accelerate convergence. In addition to the usual global convergence properties, the new algorithm possesses a local convergence rate of Q-order I + 2. Mehrotra [8] obtained the same convergence rate for his feasible-interior-point predictor-corrector linear programming algorithm. Zhang and Zhang [20] analyzed an infeasible-interior-point algorithm with I = 1 that asymptotically requires only one matrix factorization per iteration. However, they only obtained Q-order 2 convergence instead of Q-order 3.

It is well known that practical behavior of path-following algorithms such as the one analyzed here is much better if a large neighborhood is used; that is, if γ_{min} and γ_{max} are chosen to be small. The values assigned to these constants in the previous section are typical. As a referee has pointed out, however, the number of improve steps take by the algorithm may decrease as these constants are made smaller. The reason lies in the dependence of the constants C_{·} and χ on γ_{min}. The analysis of Wright [18] shows that C_{5} and C_{6} depend inversely on γ_{min}. This inverse dependence spreads to C_{14} and χ through (22) and (35), so that when γ_{min} is smaller, the threshold condition on μ_{k} — the middle condition in (36) — requires μ_{k} to be smaller before the improve steps are fully activated. This observation is not necessarily counter-intuitive, since the improved performance of the algorithm for small γ_{min} comes not so much from the improve phase as from more effective early iterations of the algorithm.

The higher-order convergence rates are mainly of theoretical interest. It is difficult to observe convergence rate higher than cubic in numerical tests, even on very small problems. However, reuse of matrix factors does tend to reduce the total number of factorizations at the cost of increasing the number of back solves, as we show in our numerical results. Since matrix factorizations are usually a good deal more expensive than back substitution, the potential reduction in computational work could be significant for large-scale problems. For linear programming, the practical effectiveness of reusing matrix factors is already well documented [5, 9].

Acknowledgment

We thank two anonymous referees for their constructive comments and suggestions.

References

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

[2] J. Ji, F. A. Potra, R. A. Tapia, and Y. Zhang, An interior-point algorithm with polynomial complexity and superlinear convergence for linear complementarity problems, Report No. 91–23, Department of Mathematical Sciences, Rice University, Houston, July 1991.

[3] N. K. Karmarkar and J. C. Lagarias and L. Slutsman and P. Wang Power– series variants of Karmarkar–type algorithms, AT & T Technical Journal, 68 (1989), pp. 20–36.

[4] M. Kojima, S. Mizuno, and A. Yoshise, An O(\sqrt{n}L) iteration potential reduction algorithm for linear complementarity problems, Mathematical Programming, 50 (1991), pp. 331–342.

[5] I. J. Lustig, R. E. Marsten, and D. F. Shanno, On implementing Mehrotra’s predictor-corrector interior point method for linear programming, SIAM Journal on Optimization, 2 (1992), pp. 435–449.

[6] O. L. Mangasarian, Error bounds for nondegenerate monotone linear complementarity problems, Mathematical Programming, 48 (1990), pp. 437–445.

[7] K. A. McShane, Superlinearly convergent O(\sqrt{n}L)-iteration interior-point algorithms for linear programming and the monotone linear complementarity problem, SIAM Journal on Optimization, 4 (1994), pp. 247–261.

[8] S. Mehrotra, Asymptotic convergence in a generalized predictor-corrector method, Technical Report, Dept. of Industrial Engineering and Management Science, Northwestern University, Evanston, Ill., October 1992.

[9] , On the implementation of a primal-dual interior point method, SIAM Journal of Optimization, 2 (1992), pp. 575–601.

[10] S. Mizuno, M. Todd, and Y. Ye, On adaptive step primal-dual interior-point algorithms for linear programming, Mathematics of Operations Research, 18 (1993), pp. 964– 981.

[11] R. D. C. Monteiro and S. J. Wright, A superlinear infeasible-interior-point affine scaling algorithm for LCP, Preprint MCS–P361–0693, Mathematics and Computer Science Division, Argonne National Laboratory, Argonne, Ill., June 1993. SIAM Journal on Optimization, to appear.

[12] J. M. Ortega and W. C. Rheinboldt, Iterative solution of nonlinear equations in several variables, Academic Press, New York, 1970.

[13] F. A. Potra, On Q-order and R-order of convergence, Journal of Optimization Theory and Applications, 63 (1989), pp. 415–431.

[14] F. A. Potra, An O(nL) infeasible-interior-point algorithm for LCP with quadratic convergence, Reports on Computational Mathematics No. 50, Department of Mathematics, University of Iowa, Iowa City, January 1994.

[15] F. A. Potra and R. Sheng, Predictor-corrector algorithms for solving P_{∗}-matrix LCP from arbitrary positive starting points, Reports on Computational Mathematics No. 58, Department of Mathematics, University of Iowa, Iowa City, August 1994.

[16] S. J. Wright, An infeasible-interior-point algorithm for linear complementarity problems, Mathematical Programming, Series A, 67 (1994), pp. 29–52.

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

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

[19] Y. Ye and K. Anstreicher, On quadratic and O(\sqrt{n}L) convergence of a predictorcorrector algorithm for LCP, Mathematical Programming, 62:3 (1993), pp. 537-552.

[20] D. Zhang and Y. Zhang, A Mehrotra-type predictor-corrector algorithm with polynomial complexity and Q-subquadratic convergence, Technical Report 93–13, Department of Mathematics and Statistics, University of Maryland Baltimore County, September 1993.

[21] Y. Zhang, On the convergence of a class of infeasible-interior-point methods for the horizontal linear complementarity problem, SIAM Journal on Optimization, 4 (1994), pp. 208–227.