Superlinear Primal-Dual Affine Scaling Algorithms for LCP

R. D. C. Monteiro^{∗}  S. J. Wright^{†}

PREPRINT MCS-P360-0693, MATHEMATICS AND COMPUTER SCIENCE DIVISION, ARGONNE NATIONAL LABORATORY

Abstract

We describe an interior-point algorithm for monotone linear complementarity problems in which primal-dual affine scaling is used to generate the search directions. The algorithm is shown to have global and superlinear convergence with Q-order up to (but not including) two. The technique is shown to be consistent with a potential-reduction algorithm, yielding the first potential-reduction algorithm that is both globally and superlinearly convergent.

Key words: interior-point methods, primal-dual affine scaling, linear programming, linear complementarity

1 Introduction

During the past three years, we have seen the appearance of several papers dealing with primal-dual interior-point algorithms for linear programs (LP) and monotone linear complementarity problems (LCP) that are superlinearly or quadratically convergent. For LP, these works include McShane [7], Mehrotra [8], Tsuchiya [15], Ye [17], Ye et al. [19], Zhang and Tapia [21], and Zhang, Tapia, and Dennis [22]. For LCP, we mention Ji, Potra, and Huang [1], Ji et al. [2], Kojima, Kurita, and Mizuno [3], Kojima, Meggido, and Noma [4], Ye and Anstreicher [18], and Zhang, Tapia, and Potra [23]. The introductory section of Ye and Anstreicher [18] contains an interesting discussion of the historical development of superlinearly convergent primal-dual interior-point algorithms.

∗School of Industrial and Systems Engineering, Georgia Institute of Technology, Atlanta, GA 30332. 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. This work was done while the first author was a faculty member of the Systems and Industrial Engineering Department at the University of Arizona.

†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.

In this paper, we discuss superlinearly convergent primal-dual affine scaling methods for solving the monotone LCP. This problem consists of finding a vector pair (x, y) ∈ IR^{n} × IR^{n} such that

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

where q ∈ IR^{n}, and M ∈ IR^{n×n} is positive semidefinite. In subsequent discussion, we say that a point (x, y) is feasible if it satisfies the equations (1a) and (1b), and strictly feasible if (1a) and (1b) are satisfied with x > 0, y > 0. We refer to (x, y) as a solution only if all three conditions in (1) hold.

Previous superlinearly convergent algorithms for (1) have required all iterates to belong to a neighborhood of the central path defined by either

N_{2}(β) = {(x, y) feasible | ∥XY e − (x^{T} y/n)e∥_{2} ≤ β(x^{T} y/n)}

or

N_{−∞}(β) = {(x, y) feasible | x_{i}y_{i} ≥ (1 − β)(x^{T} y/n), ∀i = 1, …, n} ,

where

X = diag(x),  Y = diag(y),  e = (1, 1, …, 1)^{T} ,

and β ∈ [0, 1) is a constant. For example, the predictor-corrector algorithms of Ye et al.

[19] for LP and Ji, Potra, and Huang [1] and Ye and Anstreicher [18] for LCP use the neighborhood N_{2}, while the linear programming algorithm of Zhang and Tapia [20] uses N_{−∞}. In this paper, we use a different neighborhood defined with respect to two parameters δ ≥ 0 and η > 0 by

N (δ, η) = {(x, y) strictly feasible | x_{i}y_{i} ≥ η(x^{T} y)^{1+δ}, ∀i = 1, …, n} .  (2)

Clearly, N (δ, η) is equal to the neighborhood N_{−∞}(1 − nη) when δ = 0. The parameters δ and η that define N (δ, η) do not need to be changed as the solution is approached to obtain rapid local convergence. In this respect, the neighborhood (2) differs from N_{2} and N_{−∞}, which need to be expanded during the final stages of the algorithm to achieve superlinear convergence (see Ye [17], Ye and Anstreicher [18]).

Our algorithm uses primal-dual affine scaling search directions. These directions are simply Newton steps for the system of nonlinear equations formed by (1a) and the complementarity condition XY e = 0. Because of the connection to Newton’s method, we would expect such an algorithm to be quadratically convergent if started close to a unique nondegenerate solution and allowed to take full steps. In this paper, we show that by judicious choice of the step size, all iterates will remain in N (δ, η) while simultaneously achieving fast local convergence. Depending on the choice of δ, the Q-order of the local convergence can lie anywhere in the range (1, 2). Moreover, our nondegeneracy assumption requires only that one of the solutions (x^{∗}, y^{∗}) has x^{∗} + y^{∗} > 0, not that the solution is unique.

In Section 2, we define the search directions and find bounds on the components of these directions, in terms of the complementarity gap x^{T} y and the parameters δ and η that define the neighborhood N (δ, η). We pay special attention to the case of M skew-symmetric, which occurs when (1) is derived from a linear programming problem. In this case, the bounds on the search directions are a little tighter and are global; that is, they hold everywhere in the relative interior of the set of feasible points and not just in the neighborhood N (δ, η). In Section 3, we show that a specific choice of step length yields a globally and superlinearly convergent algorithm. For particular choices of the parameter δ, the number of iterates is polynomial in the size of the problem. Finally, in Section 4, we show that the algorithm of Section 3 is consistent with a potential reduction algorithm based on the Tanabe-Todd-Ye potential function

ψ_{q}(x, y) = q log x^{T} y − ∑_{i=1}^{n} log x_{i}y_{i},  (q > n),  (3)

where an Armijo line search with a well-chosen initial trial step length is used. The resulting algorithm is again globally and superlinearly convergent.

Tunçel [16] introduced an algorithm for linear programming that uses affine scaling search directions in conjunction with a penalty function of the form

Ψ_{δ}(x, y) = (δ + 1) log(x^{T} y/n) − log(min_{j}{x_{j}y_{j}}),  (4)

for δ > 0. Step lengths are chosen to keep Ψ_{δ} constant from iteration to iteration. This function is closely related to our neighborhood N (δ, η) since

(x, y) ∈ N (δ, η) ⇔ Ψ_{δ}(x, y) ≤ − log η − (1 + δ) log n.

Tunçel [16] proves global convergence but has no superlinear convergence result. In Section 5, we prove as a consequence of our results that Tunçel’s method is superlinear for δ ∈ (0, 1). Mizuno and Nagasawa [9] describe a method for linear programming which also uses affine scaling search directions and the potential function (3). They prove complexity results for an algorithm that takes a steplength greater than a specified minimum value (defined by a formula not unlike our (35)) which does not increase (3).

The following notational conventions are used in the remainder of the paper: Unless otherwise specified, ∥.∥ denotes the Euclidean norm. For a general vector z ∈ IR^{n} and index set B ⊆ {1, …, n}, z_{B} denotes the vector made up of components z_{i} for i ∈ B. If M ∈ IR^{n×n} and B, N ⊆ {1, …, n}, then M_{B.} refers to the submatrix of M consisting of the rows i ∈ B. Similarly, M_{.N} denotes the submatrix of M corresponding to the columns j ∈ N. If D ∈ IR^{n×n} is diagonal, we write D > 0 when all the diagonal elements are strictly positive, and use the notation D_{B} to denote the diagonal submatrix constructed from D_{ii}, i ∈ B. We say that (B, N) is a partition of {1, …, n} if B ∪ N = {1, …, n} and B ∩ N = ∅.

2 Technical Results

In this section, we state our assumptions and derive bounds on components of the iterates (x^{k}, y^{k}) and the steps (Δx^{k}, Δy^{k}) for k = 0, 1, · · ·. It is assumed throughout that (x^{k}, y^{k}) lies in the neighborhood N (δ, η).

We make use of the following assumptions. The first assumption is implicit throughout the paper; the second is invoked explicitly where needed.

Assumption 1 (Existence of a strictly feasible point) The set of strictly feasible points for (1) is nonempty.

Assumption 2 (Nondegeneracy) There exists a solution (x^{∗}, y^{∗}) of (1) such that x^{∗}+y^{∗} > 0.

When (x^{∗}, y^{∗}) is the vector pair from Assumption 2, we can define a partition (B, N) of {1, …, n} by

B = {i | x_{i}^{∗} > 0}, N = {i | s_{i}^{∗} > 0}.

It is easy to show that all solutions to (1) have x_{N}^{∗} = 0 and y_{B}^{∗} = 0.

We consider in this paper the following class of primal-dual affine scaling algorithms.

Algorithm PDA

initially: Let (x^{0}, y^{0}) be a strictly feasible point;

for k = 0, 1, 2, …

Let (Δx^{k}, Δy^{k}) denote the solution of the linear system

Y Δx + XΔy = −Y Xe,  (5a)
Δy − MΔx = 0,  (5b)

where (x, y) = (x^{k}, y^{k}), X = diag(x), and Y = diag(y);

Choose α_{k} > 0 such that (x^{k} + α_{k}Δx^{k}, y^{k} + α_{k}Δy^{k}) is strictly feasible;

Set (x^{k+1}, y^{k+1}) = (x^{k} + α_{k}Δx^{k}, y^{k} + α_{k}Δy^{k});

end for

We start by stating a well-known result from linear algebra.

Lemma 2.1 Let F ∈ IR^{p×q} be given. Then there exists a nonnegative constant C = C(F) with the following property: for g ∈ IR^{p} such that the system F w = g is feasible, there exists a solution w̄ of F w = g such that

∥w̄∥ ≤ C∥g∥.

The following technical lemma unifies Theorem 2.5 and Lemma A.1 of Monteiro, Tsuchiya and Wang [13], which in turn are based on Theorem 2 of Tseng and Luo [14].

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 D > 0 and any vector h ∈ Range(H), the (unique) optimal solution w̄ = w̄(D, h) of

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

satisfies

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

Proof. We first show that |f^{T} w̄| + ∥h∥_{∞} = 0 if and only if w̄ = 0. The reverse assertion is clear, since from (6) we have h = Hw̄ = 0. For the forward assertion, note that the optimality conditions for (6) are have D^{2}w̄ + f + H^{T} u = 0 for some u ∈ IR^{p}. Since f^{T} w̄ = 0 and h = 0, we have

w̄^{T} D^{2}w̄ = −f^{T} w̄ − w̄^{T} H^{T} u = 0.

Hence, since D > 0, the result follows. Clearly the inequality (7) holds in this case. In the remainder of the proof we show that L can be chosen such that (7) also holds for nonzero solutions to (6).

Assume for contradiction that there exists a sequence of diagonal matrices {D^{k}} with D_{k} > 0 and a sequence {h^{k}} ⊆ Range(H) such that |f^{T} w^{k}| + ∥h^{k}∥_{∞} > 0 for every k and

lim_{k→∞} \frac{∥w^{k}∥_{∞}}{|f^{T} w^{k}| + ∥h^{k}∥_{∞}} = ∞,

where w^{k} is the (unique) optimal solution of (6) with D = D^{k} and h = h^{k}. By taking subsequences of subsequences, if necessary, we can identify a constant L_{1} > 0 and a nonempty index set J ⊆ {1, …, q} such that

\frac{|w_{j}^{k}|}{|f^{T} w^{k}| + ∥h^{k}∥_{∞}} ≤ L_{1},  ∀j ∉ J;  (8)

lim_{k→∞} \frac{|w_{j}^{k}|}{|f^{T} w^{k}| + ∥h^{k}∥_{∞}} = ∞,  ∀j ∈ J .  (9)

Let us define the following linear system

f^{T} w = f^{T} w^{k},  (10a)
Hw = h^{k},  (10b)
w_{j} = w_{j}^{k},  ∀j ∉ J ,  (10c)

and note that w^{k} is a solution of this system. By Lemma 2.1, (10) has a solution ŵ^{k} such that

∥ŵ^{k}∥_{∞} ≤ L_{2} {|f^{T} w^{k}| + ∥h^{k}∥_{∞} + max_{j∉J}( |w_{j}^{k}|)} ,

(the last term being omitted if J = {1, …, q}), where L_{2} ≥ 0 is a constant depending only on f, H, and J. By substituting (8) into this bound, we obtain

∥ŵ^{k}∥_{∞} ≤ L_{2} {|f^{T} w^{k}| + ∥h^{k}∥_{∞} + L_{1} (|f^{T} w^{k}| + ∥h^{k}∥_{∞})}
= L_{3} (|f^{T} w^{k}| + ∥h^{k}∥_{∞}),  (11)

where L_{3} ≡ L_{2} (1 + L_{1}). From (9), there exists K ≥ 0 such that for all k ≥ K we have

|w_{j}^{k}| > L_{3}(|f^{T} w^{k}| + ∥h^{k}∥_{∞}),  ∀j ∈ J .  (12)

Combining (11) and (12), we have

|w_{j}^{k}| > ∥ŵ^{k}∥_{∞}, ∀j ∈ J , ∀k ≥ K.

From this relation and the fact that ŵ^{k} satisfies (10c), we obtain

∥D^{k}ŵ^{k}∥^{2} = ∑_{j∈J}(D_{jj}^{k}ŵ_{j}^{k})^{2} + ∑_{j∉J}(D_{jj}^{k}ŵ_{j}^{k})^{2}
< ∑_{j∈J}(D_{jj}^{k}w_{j}^{k})^{2} + ∑_{j∉J}(D_{jj}^{k}w_{j}^{k})^{2}
= ∥D^{k}w^{k}∥^{2},  ∀k ≥ K.

Hence, since ŵ^{k} satisfies (10a), we have

f^{T} ŵ^{k} + \frac{1}{2}∥D^{k}ŵ^{k}∥^{2} < f^{T} w^{k} + \frac{1}{2}∥D^{k}w^{k}∥^{2},  ∀k ≥ K.  (13)

This relation together with the fact that ŵ^{k} satisfies (10b) contradicts the assertion that w^{k} is an optimal solution of (6) with D = D^{k} and h = h^{k}.

Our main aim now is to derive bounds on the components of the search directions generated by Algorithm PDA. We first deal with the case of M skew-symmetric. This special case is of interest mainly because the linear programming formulation

min_{w} c^{T} w, subject to Aw ≥ b, w ≥ 0,

can, with the introduction of a Lagrange multiplier λ, be posed in the form (1) with

M = [ \begin{matrix}0 & A \\ −A^{T} & 0 \end{matrix} ] , q = [ \begin{matrix}−b \\ c \end{matrix} ] .

The following two lemmas lay the foundations for a global bound on (Δx, Δy) which is found in Theorem 2.5.

Lemma 2.3 Assume that the matrix M in (1) is skew-symmetric, that is, M = −M^{T}, and that Assumption 1 holds. Then for any strictly feasible point (x, y), the solution (Δx, Δy) of (5) solves the quadratic program

min_{(w,z)} q^{T} w + \frac{1}{2}∥D^{−1}w∥^{2} + \frac{1}{2}∥Dz∥^{2},  subject to z − Mw = 0,  (14)

where D ≡ X^{1/2}Y^{−1/2}.

Proof. From (5b), (Δx, Δy) is clearly feasible for (14). It remains to verify that the optimality conditions hold, that is,

q + D^{−2}Δx = M^{T} u,  D^{2}Δy = −u,  (15)

for some u ∈ IR^{n}. By using relations (5a), (5b), and (1a) and the fact that M = −M^{T}, we obtain

q + D^{−2}Δx = q + X^{−1}Y Δx = q − y − Δy = −M(x + Δx) = M^{T} (x + Δx)

and

D^{2}Δy = Y^{−1}XΔy = −x − Δx.

Hence, (Δx, Δy) satisfies (15) with u = x + Δx.

Lemma 2.4 Assume that M in (1) is skew-symmetric. Let a strictly feasible point (x, y) be given, and consider the solution (Δx, Δy) of (5). Then,

q^{T} x = x^{T} y = −q^{T} Δx.  (16)

Proof. We have

q^{T} x = (y − Mx)^{T} x = y^{T} x − x^{T} Mx = y^{T} x,

where the last equality is due to fact that M = −M^{T}. Using this fact again, together with (5a) and (1a), we obtain

−x^{T} y = y^{T} Δx+x^{T} Δy = (Mx+q)^{T} Δx+x^{T} MΔx = q^{T} Δx+x^{T} M^{T} Δx−x^{T} M^{T} Δx = q^{T} Δx.

Theorem 2.5 Assume that M = −M^{T}. Then, for every strictly feasible point (x, y), the solution (Δx, Δy) of (5) satisfies

∥(Δx, Δy)∥_{∞} ≤ −C_{1} q^{T} Δx = C_{1} x^{T} y,  (17)

where C_{1} ≥ 0 is a constant independent of (x, y).

Proof. The bound follows directly from Lemmas 2.2, 2.3, and 2.4.

We stress that this result holds when (x, y) is any strictly feasible point, not just when (x, y) lies in the neighborhood N (δ, η). An interesting consequence of this result is that the sequence {(x^{k}, y^{k})} generated by Algorithm PDA when M is skew-symmetric always converges, regardless of the choice of step sizes {α_{k}}, though not necessarily to a solution of (1). We offer a formal proof as our next result.

Corollary 2.6 Assume that M = −M^{T}. Then the sequence {(x^{k}, y^{k})} generated by Algorithm PDA converges.

Proof. From (5) and the skew-symmetry of M, we have

0 ≤ x^{k+1}^{T} y^{k+1} = x^{k}^{T} y^{k}(1 − α_{k}) + α_{k}^{2}Δx^{k}^{T} MΔx^{k} = x^{k}^{T} y^{k}(1 − α_{k}) ≤ x^{k}^{T} y^{k}.

Hence, by Lemma 2.4, the sequence {x^{k}^{T} y^{k}} = {q^{T} x^{k}} is monotonically decreasing and bounded below by zero and therefore convergent. Using Theorem 2.5, we obtain

∥x^{k+1} − x^{k}∥ = α_{k}∥Δx^{k}∥ ≤ −C_{1}α_{k}(q^{T} Δx^{k}) = C_{1}(q^{T} x^{k} − q^{T} x^{k+1}).  (18)

Since the sequence {q^{T} x^{k}} is convergent, the above relation implies that {x^{k}} is a Cauchy sequence and therefore convergent. Therefore {y^{k}} = {Mx^{k} + q} is also convergent, and we have the result.

We now turn to the case of general positive semidefinite matrices M. We start with the following technical result which provides bounds on several quantities involving the direction (Δx, Δy). Its proof is well known and can be found in several papers (see for example Kojima, Mizuno, and Yoshise [6], Mizuno, Todd, and Ye [10], and Monteiro and Adler [12]).

Lemma 2.7 Let (x, y) be a strictly feasible point and (Δx, Δy) be the solution of (5). Then

(a) max_{1≤i≤n} |Δx_{i}Δy_{i}| ≤ x^{T} y/4;

(b) 0 ≤ Δx^{T} Δy ≤ x^{T} y/4;

(c) max {∥D^{−1}Δx∥, ∥DΔy∥} ≤ (x^{T} y)^{1/2}, where D ≡ X^{1/2}Y^{−1/2}.

We next state some simple results concerning boundedness of certain components of (x, y) and (Δx, Δy).

Lemma 2.8 Suppose that Assumption 2 holds. Then there exists a constant r > 0 such that

x_{i} ≤ (x^{T} y)/r,  ∀i ∈ N,  (19a)
y_{i} ≤ (x^{T} y)/r,  ∀i ∈ B,  (19b)

for every feasible point (x, y).

Proof. Consider a solution (x^{∗}, y^{∗}) with x^{∗} + y^{∗} > 0, as in Assumption 2. Because of monotonicity, we obtain

0 ≤ (x − x^{∗})^{T} M(x − x^{∗}) = (x − x^{∗})^{T} (y − y^{∗}) = x^{T} y − x^{∗}^{T} y − x^{T} y^{∗},

where the last equality is due to the fact that x^{∗}^{T} y^{∗} = 0. Hence, x^{T} y^{∗} ≤ x^{T} y, which implies that

x_{i} ≤ (x^{T} y)/y_{i}^{∗}, ∀i ∈ N.

Similarly, we obtain

y_{i} ≤ (x^{T} y)/x_{i}^{∗}, ∀i ∈ B,

and the result follows by choosing r to be the smallest component of (x_{B}^{∗}, y_{N}^{∗}).

Lemma 2.9 Suppose that Assumption 2 holds and that δ ≥ 0 and η > 0 are given. Then for any (x, y) ∈ N (δ, η) with corresponding (Δx, Δy) defined by (5), we have

|Δx_{i}| ≤ \frac{(x^{T} y)^{1− \frac{δ}{2}}}{rη^{1/2}} ,  ∀i ∈ N,  (20a)
|Δy_{i}| ≤ \frac{(x^{T} y)^{1− \frac{δ}{2}}}{rη^{1/2}} ,  ∀i ∈ B,  (20b)

and

x_{i} ≥ rη(x^{T} y)^{δ},  ∀i ∈ B,  (21a)
y_{i} ≥ rη(x^{T} y)^{δ},  ∀i ∈ N,  (21b)

where r is the constant from Lemma 2.8.

Proof. Let D = X^{1/2}Y^{−1/2}, and note from Lemma 2.7(c) that

∥D^{−1}Δx∥ ≤ (x^{T} y)^{1/2}.  (22)

Since (x, y) ∈ N (δ, η), we have

x_{i}y_{i} ≥ η(x^{T} y)^{1+δ},  ∀i = 1, …, n.  (23)

Using relations (22), (19a), and (23), we obtain

|Δx_{i}| ≤ x_{i} (\frac{x^{T} y}{x_{i}y_{i}} )^{1/2} ≤ (\frac{x^{T} y}{r} ) (\frac{1}{η(x^{T} y)^{δ}} )^{1/2} = \frac{(x^{T} y)^{1− \frac{δ}{2}}}{rη^{1/2}} , ∀i ∈ N,

yielding (20a). Relation (20b) follows similarly. To show (21b), we observe that relations (23) and (19a) imply

y_{i} ≥ \frac{η(x^{T} y)^{1+δ}}{x_{i}} ≥ rη(x^{T} y)^{δ}, ∀i ∈ N,

yielding (21b). Similarly, (21a) follows from (23) and (19b).

We now provide upper bounds for the remaining components of the search directions, namely, Δx_{B} and Δy_{N}. This part of the development is based on the approach of Ye and Anstreicher [18] and, in particular, the following lemma.

Lemma 2.10 (Ye and Anstreicher [18, Lemma 3.5]) Suppose that Assumption 2 holds and let (Δx, Δy) be the solution of (5). Then (u, v) = (Δx_{B}, Δy_{N}) solves the problem

min_{(u,v)} \frac{1}{2}∥D_{B}^{−1}u∥^{2} + \frac{1}{2}∥D_{N}v∥^{2},

subject to

I_{.N}v − M_{.B}u = M_{.N}Δx_{N} − I_{.B}Δy_{B}.  (24)

Our bounds for Δx_{B} and Δy_{N} are given in the following result.

Lemma 2.11 Suppose that Assumption 2 holds and that δ ≥ 0 and η > 0 are given. Then there exists a constant C_{2} > 0 independent of δ and η such that for any (x, y) ∈ N (δ, η) with corresponding (Δx, Δy) defined by (5), we have

∥Δx_{B}∥_{∞} ≤ \frac{C_{2}}{η^{1/2}} (x^{T} y)^{1− \frac{δ}{2}} ,  ∥Δy_{N}∥_{∞} ≤ \frac{C_{2}}{η^{1/2}} (x^{T} y)^{1− \frac{δ}{2}} .  (25)

Proof. Applying Lemma 2.2 to problem (24), we conclude that there exists a constant C_{3} > 0 which depends only on the matrix M and the index sets B and N such that

∥(Δx_{B}, Δy_{N})∥_{∞} ≤ C_{3}∥M_{.N}Δx_{N} − I_{.B}Δy_{B}∥_{∞}
≤ C_{3} {∥M_{.N}∥_{∞}∥Δx_{N}∥_{∞} + ∥I_{.B}∥_{∞}∥Δy_{B}∥_{∞}}
≤ C_{3} max\{∥M_{.N}∥_{∞}, ∥I_{.B}∥_{∞}\}\frac{(x^{T} y)^{1− \frac{δ}{2}}}{rη^{1/2}} ,

where the last inequality is due to relations (20a) and (20b). The result now follows by setting

C_{2} =^{△} (C_{3}/r) max{∥M_{.N}∥_{∞}, ∥I_{.B}∥_{∞}}.

We note that the result of Lemma 2.11 with δ = 0 is slightly stronger than the one obtained by Ye and Anstreicher [18, Theorem 3.6] in the sense that the constant C_{2} does not depend on the size of (x, y). Lemma 2.2 plays a crucial role in deriving this stronger version.

We can now merge the results of Theorem 2.5 with Lemmas 2.9 and 2.11 to obtain the following theorem.

Theorem 2.12 Suppose that Assumption 2 holds and that δ > 0 and η > 0 are given. Then there exists a constant C_{4} ≥ 1 independent of δ and η such that for any (x, y) ∈ N (δ, η) with (Δx, Δy) defined by (5), we have

∥Δx∥_{∞} ≤ \frac{C_{4}}{η^{ρ}} (x^{T} y)^{1−ρδ}, ∥Δy∥_{∞} ≤ \frac{C_{4}}{η^{ρ}} (x^{T} y)^{1−ρδ},

where ρ = 0 if M is skew-symmetric and ρ = 1/2 otherwise.

3 The Basic Algorithm

In this section, we develop a special version of a primal-dual affine scaling algorithm for which all iterates lie in a neighborhood N (δ, η), for some δ > 0 and η > 0. We show that when δ ∈ (0, 1/2) (δ ∈ (0, 1), if M is skew-symmetric), the algorithm is superlinearly convergent with q-order equal to 2 − 2δ (respectively, 2 − δ). We also derive a bound on the number of iterations to reduce the duality gap below a specified tolerance. For certain choices of δ, the algorithm is shown to have a polynomial bound on the total number of iterations.

The following notation and definitions will be used in the remainder of the paper. Given the strictly feasible point (x, y) and the search direction (Δx, Δy) from (5), define

(x(α), y(α)) =^{△} (x, y) + α(Δx, Δy);
ϕ =^{△} \frac{Δx^{T} Δy}{x^{T} y};  (26)

χ =^{△} max_{{i=1,···,n \atop Δx_{i}Δy_{i}<0}} (\frac{|Δx_{i}Δy_{i}|}{x_{i}y_{i}}) ,

where we assume that the maximum over the empty set is 0. We use ϕ_{k} and χ_{k} for the values of ϕ and χ at (x, y) = (x^{k}, y^{k}) and (Δx, Δy) = (Δx^{k}, Δy^{k}).

Lemma 3.1 Let (x, y) be a strictly feasible point and (Δx, Δy) be the solution of (5). Then the following statements hold for all α ∈ [0, 1]:

(a) x(α)^{T} y(α) = x^{T} y(1 − α + α^{2}ϕ);

(b) x_{i}(α)y_{i}(α) ≥ x_{i}y_{i}(1 − α − α^{2}χ);

(c) (1 − α)x^{T} y ≤ x(α)^{T} y(α) ≤ (1 − α/2)^{2}x^{T} y;

(d) Δx_{i}Δy_{i} = 0 for all i = 1, ⋯, n if and only if (x(1), y(1)) is a solution of (1) with x(1) + y(1) > 0.

Proof. Statements (a) and (b) follow immediately from (5a) and (26). Lemma 2.7(b) implies that 0 ≤ ϕ ≤ 1/4, which together with (a) implies statement (c).

To prove (d), we assume first that Δx_{i}Δy_{i} = 0 for all i = 1, ⋯, n, so that ϕ = 0. Using (5a), we can partition {1, ⋯, n} into index sets B and N such that

i ∈ B ⇒ Δy_{i} = −y_{i}, Δx_{i} = 0,
i ∈ N ⇒ Δx_{i} = −x_{i}, Δy_{i} = 0.

Therefore

i ∈ B ⇒ y_{i}(1) = 0, x_{i}(1) = x_{i} > 0,
i ∈ N ⇒ x_{i}(1) = 0, y_{i}(1) = y_{i} > 0,

and so (x(1), y(1)) satisfies (1b) and (1c) and x(1) + y(1) > 0, The equations y = Mx + q and (5b) imply that (1a) is also satisfied, so the forward implication is true. For the converse, assume (x(1), y(1)) solution with x(1) + y(1) > 0. Partitioning {1, ⋯, n} into B = {i | x_{i}(1) > 0} and N = {i | y_{i}(1) > 0}, we have using (5a) that

i ∈ B ⇒ Δy_{i} = −y_{i}, Δx_{i} = 0,
i ∈ N ⇒ Δx_{i} = −x_{i}, Δy_{i} = 0.

Hence, the converse is proved, and so (d) holds.

Lemma 3.2 If u ≤ 1 and δ ∈ (0, 1], then (1 − u)^{1+δ} ≤ 1 − (1 + δ)u + δu^{2}.

Proof. By the mean value theorem, we have that

(1 − u)^{δ} = 1 − δu + (δ − 1)δ(1 − ū)^{δ−2} \frac{u^{2}}{2}

for some ū between 0 and u. The last term in the above expression is nonpositive, and therefore (1 − u)^{δ} ≤ 1 − δu. Hence,

(1 − u)^{1+δ} = (1 − u)(1 − u)^{δ} ≤ (1 − u)(1 − δu) = 1 − (1 + δ)u + δu^{2},

giving the result.

Lemma 3.3 Let δ ∈ (0, 1] and η > 0 be given, and assume that (x, y) ∈ N (δ, η). Let (Δx, Δy) be defined by (5), and suppose that Δx_{i}Δy_{i} ≠ 0 for some i = 1, ⋯, n, so that (x(1), y(1)) does not solve (1). Then, for every

α ∈ J =^{△} [0, \frac{δ}{δ + χ + (1 + δ)ϕ}] ,

the following statements hold:

(a) 1 − α − α^{2}χ > 0;

(b) (x(α), y(α)) ∈ N (δ, η).

Proof. Since Δx_{i}Δy_{i} ≠ 0 for some i, we can easily see that either χ > 0 or ϕ > 0, or both. Hence, J ⊆ [0, 1). In view of Lemma 3.1(c) and x^{T} y > 0, this implies that x(α)^{T} y(α) > 0 for all α ∈ J. Using Lemma 3.1(a), Lemma 3.2, and the inclusions ϕ ∈ [0, 1/4], α ∈ [0, 1), we obtain

0 < η [x(α)^{T} y(α)]^{1+δ}
= η(x^{T} y)^{1+δ}[1 − α(1 − αϕ)]^{1+δ}
≤ η(x^{T} y)^{1+δ}[1 − (1 + δ)α(1 − αϕ) + δα^{2}(1 − αϕ)^{2}]
≤ η(x^{T} y)^{1+δ}[1 − (1 + δ)α(1 − αϕ) + δα^{2}]
≤ η(x^{T} y)^{1+δ}[1 − α − α^{2}χ],  (27)

where in the last inequality we used the fact that the interval J is exactly the set of all α ≥ 0 for which

1 − α(1 + δ)(1 − αϕ) + δα^{2} ≤ 1 − α − α^{2}χ.

Clearly, (27) implies statement (a).

Since (x, y) ∈ N (δ, η), we have

x_{i}y_{i} ≥ η(x^{T} y)^{1+δ},  ∀i = 1, ⋯, n.  (28)

Using this relation together with Lemma 3.1(b) and statement (a), we obtain

x_{i}(α)y_{i}(α) ≥ x_{i}y_{i}(1 − α − α^{2}χ) ≥ η(x^{T} y)^{1+δ}(1 − α − α^{2}χ),  ∀i = 1, ⋯, n.  (29)

Hence, it follows from (27) and (29) that

min_{i=1,···,n} {x_{i}(α)y_{i}(α)} ≥ η[x(α)^{T} y(α)]^{1+δ} > 0, ∀α ∈ J.

Thus, (x(α), y(α)) ∈ N (δ, η) for all α ∈ J, and the result follows.

Lemma 3.4 Assume that (x, y) ∈ N (δ, η), where δ ∈ (0, 1] and η > 0 are given constants. Then

η(x^{T} y)^{δ} ≤ \frac{1}{n}.  (30)

Proof. Since (x, y) ∈ N (δ, η), we have

x^{T} y = ∑_{i=1}^{n} x_{i}y_{i} ≥ nη(x^{T} y)^{1+δ},

giving the result.

Lemma 3.5 Suppose that δ ∈ (0, 1] and η > 0 are given. There exists a constant C̄ > 0 independent of δ and η such that if (x, y) ∈ N (δ, η), then

max\{ϕ, χ\} ≤ C̄ \frac{(x^{T} y)^{1−δξ}}{3η^{ξ}} ,  (31)

where ξ = 1 if the matrix M is skew-symmetric and ξ = 2 otherwise.

Proof. Consider the constant C_{4} as in the statement of Theorem 2.12, and define C̄ = 3C_{4}^{2}. We will show that C̄ fulfills the requirements of the lemma. For any (x, y) ∈ N (δ, η), Theorem 2.12 and Lemma 3.4 imply that

ϕ = \frac{Δx^{T} Δy}{x^{T} y}
≤ \frac{nC_{4}^{2} [ (x^{T} y)^{2(1−ρδ)} ]/η^{2ρ}}{x^{T} y}

= nη(x^{T} y)^{δ}C_{4}^{2} \frac{(x^{T} y)^{1−(2ρ+1)δ}}{η^{2ρ+1}}
= nη(x^{T} y)^{δ}C_{4}^{2} \frac{(x^{T} y)^{1−δξ}}{η^{ξ}}
≤ C̄ \frac{(x^{T} y)^{1−δξ}}{3η^{ξ}} ,  (32)

where in the last equality we have used the fact that 2ρ + 1 = ξ. For a bound on χ, we use (x, y) ∈ N (δ, η) and Theorem 2.12 to deduce that

χ ≤ max_{1≤i≤n} \{\frac{|Δx_{i}Δy_{i}|}{x_{i}y_{i}}\} ≤ \frac{C_{4}^{2} [ (x^{T} y)^{2(1−ρδ)} ]/η^{2ρ}}{η(x^{T} y)^{1+δ}} = C_{4}^{2} \frac{(x^{T} y)^{1−(2ρ+1)δ}}{η^{2ρ+1}} = C̄ \frac{(x^{T} y)^{1−δξ}}{3η^{ξ}} . (33)

Clearly, (32) and (33) imply the result.

We now state the main theorem concerning convergence of the algorithm. For the purpose of this theorem and the next section we define a measure of centrality of the initial point as follows:

π_{0} = \frac{min_{i=1,···,n}\{x_{i}^{0}y_{i}^{0}\}}{x^{0}^{T} y^{0}/n}.  (34)

Theorem 3.6 Let δ ∈ (0, 1/ξ) and the strictly feasible point (x^{0}, y^{0}) be given, where ξ is the constant defined in Lemma 3.5. Suppose that Algorithm PDA is applied with step lengths

α_{k} = \frac{δ}{δ + χ_{k} + (1 + δ)ϕ_{k}} .  (35)

Then

(a) the sequence {x^{k}^{T} y^{k}} converges monotonically to 0;

(b) given any ϵ ∈ (0, x^{0}^{T} y^{0}), the algorithm achieves x^{k}^{T} y^{k} ≤ ϵ in at most

O ((\frac{x^{0}^{T} y^{0}}{ϵ} )^{δ} (\frac{n}{δπ_{0}} ) log \frac{x^{0}^{T} y^{0}}{ϵ} )  (36)

iterations;

(c) either the algorithm terminates with x^{k}^{T} y^{k} = 0 in a finite number of iterations, or else the sequence {x^{k}^{T} y^{k}} converges to zero superlinearly with Q-order 2 − δξ.

Proof. First, we observe that there is an index K ≥ 0 with Δx_{i}^{K}Δy_{i}^{K} = 0 for all i = 1, ⋯, n if and only if (x^{K+1}, y^{K+1}) is a solution of (1) (see Lemma 3.1(d)). Moreover, using the proof of statement (b) below, one can easily see that if such an index K exists, then it must be of the order in (36). We will henceforth assume that finite termination does not occur and so for all k ≥ 0, there is at least one i = 1, ⋯, n such that Δx_{i}^{k}Δy_{i}^{k} ≠ 0.

(a) Choosing any η > 0 such that (x^{0}, y^{0}) ∈ N (δ, η), we have from relation (35) and Lemma 3.3 that (x^{k}, y^{k}) ∈ N (δ, η) for all k ≥ 0. We now find a lower bound on the step length α_{k}. By Lemma 3.1(c) and the fact that α_{k} ∈ [0, 1] for all k ≥ 0, we have that {x^{k}^{T} y^{k}} is a decreasing sequence. Hence, using (31), we obtain

max\{ϕ_{k}, χ_{k}\} ≤ C̄ \frac{(x^{k}^{T} y^{k})^{1−δξ}}{3η^{ξ}} ≤ C̄ \frac{(x^{0}^{T} y^{0})^{1−δξ}}{3η^{ξ}} , ∀k ≥ 0,

which, in view of (35), implies

α_{k} ≥ ᾱ =^{△} \frac{δ}{δ + C̄(x^{0}^{T} y^{0})^{1−δξ}/η^{ξ}} > 0, ∀k ≥ 0.

Using Lemma 3.1(c) again, we get x^{k+1}^{T} y^{k+1} ≤ (1 − ᾱ/2)^{2}x^{k}^{T} y^{k}, ∀k ≥ 0, which clearly implies statement (a).

(b) Let

η = \frac{π_{0}}{n(x^{0}^{T} y^{0})^{δ}} ,  (37)

where π_{0} is defined in (34). Note that η is the largest possible choice for which (x^{0}, y^{0}) ∈ N (δ, η). By relation (35) and Lemma 3.3, we know that (x^{k}, y^{k}) ∈ N (δ, η) for all k ≥ 0. Assuming that x^{k}^{T} y^{k} ≥ ϵ, we now seek bounds on ϕ_{k} and χ_{k} that are independent of the somewhat murky constant C̄. Using Lemma 2.7(a) and the fact that (x^{k}, y^{k}) ∈ N (δ, η), we obtain

χ_{k} ≤ max_{1≤i≤n} \{\frac{|Δx_{i}^{k}Δy_{i}^{k}|}{x_{i}^{k}y_{i}^{k}}\} ≤ \frac{x^{k}^{T} y^{k}/4}{η(x^{k}^{T} y^{k})^{1+δ}} = \frac{1}{4η(x^{k}^{T} y^{k})^{δ}} .

Hence, from (37) and the assumption that x^{k}^{T} y^{k} ≥ ϵ, we obtain

χ_{k} ≤ \frac{1}{4η(x^{k}^{T} y^{k})^{δ}} ≤ \frac{n}{4π_{0}(ϵ/x^{0}^{T} y^{0})^{δ}} .

Moreover, Lemma 2.7(b) implies that ϕ_{k} ≤ 1/4. These two last relations together with (35) then imply that

α_{k} ≥ \frac{δ}{δ + n/[4π_{0}(ϵ/x^{0}^{T} y^{0})^{δ}] + (1/2)} ≥ \frac{δπ_{0}(ϵ/x^{0}^{T} y^{0})^{δ}}{2n} ,

where in the last inequality we used the fact that π_{0} ≤ 1 and ϵ < x^{0}^{T} y^{0}. Hence, by Lemma 3.1(c), we conclude that

\frac{x^{k+1}^{T} y^{k+1}}{x^{k}^{T} y^{k}} ≤ (1 − \frac{δπ_{0}(ϵ/x^{0}^{T} y^{0})^{δ}}{4n} )^{2} .  (38)

We have thus shown that (38) holds whenever x^{k}^{T} y^{k} ≥ ϵ. Now, let K be the smallest nonnegative integer k for which x^{k+1}^{T} y^{k+1} ≤ ϵ. Assume for contradiction that

K > (\frac{x^{0}^{T} y^{0}}{ϵ} )^{δ} (\frac{2n}{δπ_{0}} ) log \frac{x^{0}^{T} y^{0}}{ϵ} .

By taking logarithms and using the inequality log(1 − β) ≤ −β for β ∈ (0, 1), we can show that

(1 − \frac{δπ_{0}(ϵ/x^{0}^{T} y^{0})^{δ}}{4n} )^{2K} < \frac{ϵ}{x^{0}^{T} y^{0}} .

Using this relation and the fact that (38) holds for every k ≤ K, we obtain

x^{K}^{T} y^{K} ≤ (1 − \frac{δπ_{0}(ϵ/x^{0}^{T} y^{0})^{δ}}{4n} )^{2K} x^{0}^{T} y^{0} < ϵ,

which contradicts the fact that x^{K}^{T} y^{K} ≥ ϵ. Therefore (36) holds, and we have proved statement (b).

(c) Choose η as in the proof of part (a). From Lemma 3.1(a), (35), and (31), we have

x^{k+1}^{T} y^{k+1} = x^{k}^{T} y^{k}(1 − α_{k} + α_{k}^{2}ϕ_{k})
≤ x^{k}^{T} y^{k} (\frac{χ_{k} + (1 + δ)ϕ_{k}}{δ + χ_{k} + (1 + δ)ϕ_{k}} + ϕ_{k})

≤ x^{k}^{T} y^{k} (\frac{C̄(x^{k}^{T} y^{k})^{1−δξ}/η^{ξ}}{δ} + \frac{C̄(x^{k}^{T} y^{k})^{1−δξ}}{3η^{ξ}} )

≤ \frac{4C̄}{3δη^{ξ}} (x^{k}^{T} y^{k})^{2−δξ} ,

which implies statement (c).

We observe that for certain values of δ, the algorithm of Theorem 3.6 has a polynomial bound on the number of iterations. To see this, assume that L > 0 is such that 2^{L} ≥ max{x^{0}^{T} y^{0}, ϵ^{−1}}. Assume also that the initial point (x^{0}, y^{0}) is chosen in such a way that π_{0} = O(1) is independent of n and L. Then, in terms of n, L, and δ, the algorithm of Theorem 3.6 terminates in O(nL2^{2Lδ}/δ) iterations. If we choose δ = O(1/L), then the number of iterations is O(nL^{2}).

4 A Superlinearly Convergent Potential Reduction Algorithm

Although our algorithm has excellent local convergence properties, preliminary computations have shown that its behavior at points remote from the solution is poor. It is therefore worthwhile to merge our method with other methods with more attractive global convergence properties. In this section, we embed our method in a potential reduction method. The resulting method retains the global convergence behavior of potential reduction methods while exhibiting the fast local convergence associated with the algorithm of Section 3. We are also motivated by a desire to show that potential reduction methods can be superlinearly convergent.

Our potential reduction algorithm is based on the Tanabe-Todd-Ye potential function

ψ_{q}(x, y) = q log x^{T} y − ∑_{i=1}^{n} log x_{i}y_{i},  (q > n).  (39)

We start by specifying the algorithm and stating a global convergence theorem. Then we show that if q lies in the range (n, n + 1/ξ), the algorithm becomes compatible with the method described in the previous section and is superlinearly convergent.

We start by specifying the algorithm.

Algorithm PDPR

initially: Choose fixed constants β > 1 and μ ∈ (0, 1). Let {τ_{k}} be a scalar sequence satisfying τ_{k} ≥ 1/2, and let (x^{0}, y^{0}) be strictly feasible;

for k = 0, 1, 2, …

Find (Δx^{k}, Δy^{k}) by solving (5) with (x, y) = (x^{k}, y^{k});

Define α_{k} = β^{−m_{k}}τ_{k}, where m_{k} is the smallest nonnegative integer for which (x^{k}(β^{−m_{k}}τ_{k}), y^{k}(β^{−m_{k}}τ_{k})) > 0 and

ψ_{q}(x^{k}(β^{−m_{k}}τ_{k}), y^{k}(β^{−m_{k}}τ_{k})) ≤ ψ_{q}(x^{k}, y^{k}) + μ(β^{−m_{k}}τ_{k})∇ψ_{q}(x^{k}, y^{k})^{T} [ \begin{matrix}Δx^{k} \\ Δy^{k} \end{matrix} ];  (40)

Set (x^{k+1}, y^{k+1}) = (x^{k}(α_{k}), y^{k}(α_{k}));

end for

This algorithm is closely related to the algorithm of Monteiro [11] for convex programming. Its global convergence properties can be analyzed by using techniques like those of Monteiro [11] and Kojima et al. [5]. We omit the details of this analysis and simply state the final theorem.

Theorem 4.1 Suppose that Assumption 2 holds. Then the sequence of iterates {(x^{k}, y^{k})} generated by Algorithm PDPR satisfies

lim_{k→∞} x^{k}^{T} y^{k} = 0.

The following assumption ensures that Algorithm PDPR enters a superlinear phase in which the local convergence rate of the algorithm of Theorem 3.6 is attained.

Assumption 3 There is an integer K > 0 such that the iterate (x^{K}, y^{K}) of Algorithm PDPR satisfies

min_{i} x_{i}^{K}y_{i}^{K} ≥ (\frac{C̄}{δ} )^{1/ξ} (x^{K}^{T} y^{K})^{1+1/ξ} ,  (41)

where C̄ and ξ are defined in Lemma 3.5.

Suppose that for some δ ∈ (0, 1/ξ), the sequence of initial step sizes {τ_{k}} is selected as

τ_{k} = max (\frac{1}{2}, \frac{δ}{δ + χ_{k} + (1 + δ)ϕ_{k}} ) ,  ∀k ≥ 0.  (42)

We have the following result.

Lemma 4.2 Suppose that Assumptions 2 and 3 hold, and let {(x^{k}, y^{k})} be the sequence of iterates generated by Algorithm PDPR with the initial step size τ_{k} defined by (42) for all k ≥ 0. Define

η = \frac{min x_{i}^{K}y_{i}^{K}}{(x^{K}^{T} y^{K})^{1+δ}} ,  (43)

where K is the index in Assumption 3. Then, for all k ≥ K, we have

(a) (x^{k}, y^{k}) ∈ N (δ, η); and

(b) τ_{k} = δ/(δ + χ_{k} + (1 + δ)ϕ_{k}).

Proof. Recall that by Lemma 3.1(c), the sequence {x^{k}^{T} y^{k}} decreases. Using this fact together with relations (41) and (43), one can easily verify that

\frac{C̄(x^{k}^{T} y^{k})^{1−ξδ}}{η^{ξ}} ≤ δ,  ∀k ≥ K.  (44)

We show first that (a) implies (b) for each k ≥ K. Assuming (a), we have from (31) and (44) that

\frac{χ_{k} + (1 + δ)ϕ_{k}}{δ + χ_{k} + (1 + δ)ϕ_{k}} ≤ \frac{C̄(x^{k}^{T} y^{k})^{1−ξδ}/η^{ξ}}{δ + C̄(x^{k}^{T} y^{k})^{1−ξδ}/η^{ξ}} ≤ \frac{δ}{2δ} = \frac{1}{2},

and therefore, from (42), statement (b) holds.

We now prove (a) by induction. Clearly, by the definition (43), (a) holds for k = K. Assume now that (a) is satisfied at some arbitrary k ≥ K. Then (b) also holds at the index k, and we can use Lemma 3.3(b) and the fact that α_{k} ∈ (0, τ_{k}] to conclude that (x^{k+1}, y^{k+1}) ∈ N (δ, η), giving the result.

A result like Theorem 3.6(c) ensures superlinear convergence provided we show that m_{k} = 0 for all k sufficiently large, that is, α_{k} = τ_{k}. We do so in the following theorem.

Theorem 4.3 Suppose that Assumptions 2 and 3 hold, and let {(x^{k}, y^{k})} be the sequence of iterates generated by Algorithm PDPR with the initial step size τ_{k} defined by (42) for some δ ∈ (0, 1/ξ). Assume also that τ_{k} < 1 for all k. Then the sequence {x^{k}^{T} y^{k}} converges to zero superlinearly with Q-order 2 − δξ.

Proof. Global convergence of the sequence {x^{k}^{T} y^{k}} to zero follows from Theorem 4.1. We next show that α_{k} = τ_{k} for all k ≥ K sufficiently large, from which superlinear convergence of {x^{k}^{T} y^{k}} to zero follows as a consequence of Lemma 4.2(b) and Theorem 3.6(c). To show that α_{k} = τ_{k}, we need only show that

(x^{k}(τ_{k}), y^{k}(τ_{k})) > 0,  (45)

and

ψ_{q}(x^{k}(τ_{k}), y^{k}(τ_{k})) − ψ_{q}(x^{k}, y^{k}) ≤ μτ_{k}∇ψ_{q}(x^{k}, y^{k})^{T} [ \begin{matrix}Δx^{k} \\ Δy^{k} \end{matrix} ] .  (46)

Lemma 4.2(a) implies that (45) holds for all k ≥ K. To prove that (46) holds for all k ≥ K sufficiently large, we start by finding a lower bound on its right-hand side. Using (39) and (5a), we have

μτ_{k}∇ψ_{q}(x^{k}, y^{k})^{T} [ \begin{matrix}Δx^{k} \\ Δy^{k} \end{matrix} ]

= μτ_{k} \{[\frac{q}{x^{k}^{T} y^{k}} y^{k} − (X^{k})^{−1}e]^{T} Δx^{k} + [\frac{q}{x^{k}^{T} y^{k}} x^{k} − (Y^{k})^{−1}e]^{T} Δy^{k}\}

= μτ_{k} \{\frac{q}{x^{k}^{T} y^{k}} (y^{k}^{T} Δx^{k} + x^{k}^{T} Δy^{k}) − ∑_{i=1}^{n} [\frac{Δx_{i}^{k}}{x_{i}^{k}} + \frac{Δy_{i}^{k}}{y_{i}^{k}} ]\}
= −μτ_{k}(q − n) ≥ −μ(q − n).  (47)

To find an upper bound for the left-hand side of (46), note that for all k ≥ K, we have

\frac{τ_{k}^{2} max(ϕ_{k}, χ_{k})}{1 − τ_{k}} = (\frac{δ}{δ + χ_{k} + (1 + δ)ϕ_{k}} )^{2} (\frac{δ + χ_{k} + (1 + δ)ϕ_{k}}{χ_{k} + (1 + δ)ϕ_{k}} ) max(ϕ_{k}, χ_{k})

= (\frac{δ^{2}}{δ + χ_{k} + (1 + δ)ϕ_{k}} ) \frac{max(ϕ_{k}, χ_{k})}{χ_{k} + (1 + δ)ϕ_{k}}
≤ δ.

Using this bound together with statements (a) and (b) of Lemma 3.1, we obtain

ψ_{q}(x^{k}(τ_{k}), y^{k}(τ_{k})) − ψ_{q}(x^{k}, y^{k})
= q log [\frac{x^{k}(τ_{k})^{T} y^{k}(τ_{k})}{x^{k}^{T} y^{k}}] − ∑_{i=1}^{n} log [\frac{x_{i}^{k}(τ_{k})y_{i}^{k}(τ_{k})}{x_{i}^{k}y_{i}^{k}}]

≤ q log(1 − τ_{k} + τ_{k}^{2}ϕ_{k}) − ∑_{i=1}^{n} log(1 − τ_{k} − τ_{k}^{2}χ_{k})
= (q − n) log(1 − τ_{k}) + q log (1 + \frac{τ_{k}^{2}ϕ_{k}}{1 − τ_{k}} ) − n log (1 − \frac{τ_{k}^{2}χ_{k}}{1 − τ_{k}} )
≤ (q − n) log(1 − τ_{k}) + q log(1 + δ) − n log(1 − δ).  (48)

Since x^{k}^{T} y^{k} → 0, it follows from (31) and (42) that 1 − τ_{k} → 0. Therefore, from (48), the left-hand side of (46) approaches −∞ as k → ∞. Hence, in view of (47), the inequality (46) will hold for all sufficiently large k ≥ K.

The proof of superlinearity now follows as in Theorem 3.6(c).

Note that if τ_{k} = 1 for some k, we have Δx_{i}^{k}Δy_{i}^{k} = 0 for all i; therefore, from Lemma 3.1(d), (x^{k+1}, y^{k+1}) is an exact solution of (1). In this case, Algorithm PDPR will reject the step α_{k} = τ_{k} because (45) is not satisfied! Any implementation would surely recognize this special case, so we have avoided the complication that it introduces into the analysis above.

Obviously, we cannot explicitly identify the index K required by Assumption 3, since we do not know the value of C̄ in general. However, we can be sure that Assumption 3 holds if we choose q to be in the interval (n, n + 1/ξ), as we show in the following theorem.

Theorem 4.4 Suppose that Assumptions 1 and 2 hold and that q ∈ (n, n+1/ξ). Then there is an index K such that (41) holds at the iterate (x^{K}, y^{K}) generated by Algorithm PDPR.

Proof. Let us define Q and Q̄ as

Q =^{△} ψ_{q}(x^{0}, y^{0}),  Q̄ =^{△} Q − (n − 1) log(n − 1).

Since ψ_{q} is reduced at each iteration of Algorithm PDPR, we have ψ_{q}(x^{k}, y^{k}) ≤ Q for all k ≥ 0. Hence, for all k ≥ 0 and i = 1, ⋯, n, we have

(q − n + 1) log x^{k}^{T} y^{k} − log x_{i}^{k}y_{i}^{k}
≤ Q − (n − 1) log x^{k}^{T} y^{k} + ∑_{j≠i} log x_{j}^{k}y_{j}^{k}
≤ Q − (n − 1) log (x^{k}^{T} y^{k} − x_{i}^{k}y_{i}^{k}) + ∑_{j≠i} log x_{j}^{k}y_{j}^{k}
≤ Q − (n − 1) log(n − 1) = Q̄,

where the third inequality follows from the arithmetic-geometric mean inequality, namely: (∑_{i=1}^{p} a_{i})/p ≥ (Π_{i=1}^{p}a_{i})^{1/p} for any positive scalars a_{1}, …, a_{p}. Hence,

log \frac{x_{i}^{k}y_{i}^{k}}{(x^{k}^{T} y^{k})^{q−n+1}} ≥ −Q̄ ⇒ x_{i}^{k}y_{i}^{k} ≥ e^{−Q̄} (x^{k}^{T} y^{k})^{q−n+1} ,

and so

x_{i}^{k}y_{i}^{k} ≥ (\frac{C̄}{δ} )^{1/ξ} (x^{k}^{T} y^{k})^{1+1/ξ} \{e^{−Q̄} (\frac{δ}{C̄} )^{1/ξ} (x^{k}^{T} y^{k})^{q−n−1/ξ}\} .

Since q − n − 1/ξ < 0 and x^{k}^{T} y^{k} → 0, there is a K > 0 such that the last bracketed term in the expression is greater than 1 for all k ≥ K, giving the result.

Finally, we combine the results of the last two theorems to obtain the following result.

Corollary 4.5 Suppose that Assumption 2 holds, that q ∈ (n, n + 1/ξ), that the sequence of initial step sizes {τ_{k}} is defined by (42) for some δ ∈ (0, 1/ξ), and that τ_{k} < 1 for all k. Then the sequence {x^{k}^{T} y^{k}} generated by Algorithm PDPR converges to zero superlinearly with Q-order 2 − δξ.

5 Concluding Remarks

In this concluding section, we discuss the close relationship that exists between the algorithm of Theorem 3.6 and the one presented by Tunçel [16]. We also show that Tunçel’s algorithm is superlinearly convergent whenever the parameter δ that appears in the potential function (4) lies in the interval (0, 1).

For simplicity, we assume that the affine scaling direction (Δx, Δy) at any strictly feasible point satisfies Δx_{i}Δy_{i} ≠ 0, for some i ∈ {1, ⋯, n}, so that finite termination of the algorithm never occurs. For the purpose of this section, we also assume that the matrix M in (1) is skew-symmetric. Since Tunçel’s algorithm is for linear programs, it fits into this framework.

We describe Tunçel’s algorithm in an equivalent but slightly different way from [16] in order to better point out its connection to the neighborhood N (δ, η). Let the strictly feasible point (x^{k}, y^{k}) be the k-th iterate of the algorithm. Tunçel’s algorithm sets (x^{k+1}, y^{k+1}) = (x^{k}(α_{k}), y^{k}(α_{k})), where the stepsize α_{k} ∈ (0, 1) is such that the point (x^{k}(α_{k}), y^{k}(α_{k})) satisfies the equation in (x, y) defined by

\frac{min_{i=1,···,n}\{x_{i}y_{i}\}}{(x^{T} y)^{1+δ}} = \frac{min_{i=1,···,n}\{x_{i}^{k}y_{i}^{k}\}}{(x^{k}^{T} y^{k})^{1+δ}}.

Tunçel shows that such a stepsize always exists and is unique. In the proof of this fact, he also shows that (x^{k}(α), y^{k}(α)) satisfies the strict inequality

\frac{min_{i=1,···,n}\{x_{i}^{k}(α)y_{i}^{k}(α)\}}{[x^{k}(α)^{T} y^{k}(α)]^{1+δ}} > \frac{min_{i=1,···,n}\{x_{i}^{k}y_{i}^{k}\}}{(x^{k}^{T} y^{k})^{1+δ}}

for all α ∈ (0, α_{k}). Observe that the iterates of this algorithm satisfy

\frac{min_{i=1,···,n}\{x_{i}^{k}y_{i}^{k}\}}{(x^{k}^{T} y^{k})^{1+δ}} = \frac{min_{i=1,···,n}\{x_{i}^{0}y_{i}^{0}\}}{(x^{0}^{T} y^{0})^{1+δ}} =^{△} η, ∀k ≥ 0,

which, in terms of the potential function (4), is equivalent to Ψ_{δ}(x^{k}, y^{k}) = Ψ_{δ}(x^{0}, y^{0}) for all k ≥ 0. It is now easy to see that Tunçel’s step size is the largest α > 0 for which (x^{k}(α), y^{k}(α)) ∈ N (δ, η). In view of Lemma 3.3(b), it follows that our step size (35) is less than or equal to Tunçel’s step size. Therefore Tunçel’s algorithm achieves larger or equal reduction of the duality gap at each iteration while generating all points within the neighborhood N (δ, η). Hence, the same proof given for Theorem 3.6(c) can be used to show that Tunçel’s algorithm is superlinearly convergent with q-order of convergence equal to 2−δ, whenever δ ∈ (0, 1).

Acknowledgment

Part of this work was done while the first author was visiting the Mathematics and Computer Science Division of Argonne National Laboratory. He is grateful to Argonne for its financial support and its congenial scientific atmosphere.

References

[1] J. Ji, F. A. Potra, and S. Huang, A predictor-corrector method for linear complementarity problems with polynomial complexity and superlinear convergence, Technical Report 18, Department of Mathematics, University of Iowa, Iowa City, Iowa, 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, Technical Report 91–23, Department of Mathematical Sciences, Rice University, Houston, TX, July 1991.

[3] M. Kojima, Y. Kurita, and S. Mizuno, Large–step interior point algorithms for linear complementarity problems, SIAM Journal on Optimization, 3 (1993), pp. 398–412.

[4] M. Kojima, N. Megiddo, and T. Noma, Homotopy continuation methods for nonlinear complementarity problems, Mathematics of Operations Research, 16 (1991), pp. 754–774.

[5] M. Kojima, N. Megiddo, T. Noma, and A. Yoshise, A unified approach to interior-point algorithms for linear complementarity problems, no. 538 in Lecture Notes in Computer Science, Springer-Verlag, Berlin, 1991.

[6] M. Kojima, S. Mizuno, and A. Yoshise, A polynomial-time algorithm for a class of linear complementarity problems, Mathematical Programming, 44 (1989), pp. 1–26.

[7] K. A. McShane, A superlinearly convergent O(\sqrt{n}L) iteration primal–dual linear programming algorithm, SIAM Journal on Optimization, 4 (1994), pp. 247–261.

[8] S. Mehrotra, Quadratic convergence in a primal-dual method, Technical Report 91– 15, Dept. of Industrial Engineering and Management Science, Northwestern University, Evanston, IL 60208, USA, September 1991.

[9] S. Mizuno and A. Nagasawa, A primal-dual affine scaling potential reduction algorithm for linear programming, Research Memorandum 427, The Institute of Statistical Mathematics, Minami-Azabu, Minato-ku, Tokyo, Japan, February 1992.

[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, A globally convergent interior point algorithm for convex programming, Technical Report 91–021, Dept. of Systems and Industrial Engineering, University of Arizona, July 1991.

[12] R. D. C. Monteiro and I. Adler, Interior path-following primal-dual algorithms. Part I: Linear programming, Mathematical Programming, 44 (1989), pp. 27–41.

[13] R. D. C. Monteiro, T. Tsuchiya, and Y. Wang, A simplified global convergence proof of the affine scaling algorithm, Annals of Operations Research, 47 (1993), pp. 443– 482.

[14] P. Tseng and Z. Q. Luo, On the convergence of the affine–scaling algorithm, Mathematical Programming, 56 (1992), pp. 301–319.

[15] T. Tsuchiya, Quadratic convergence of Iri and Imai’s method for degenerate linear programming problems, Research Memorandum 412, The Institute of Statistical Mathematics, 4–6–7 Minami–Azabu, Minato–ku, Tokyo 106, Japan, 1991.

[16] L. Tuncel, Constant potential primal-dual algorithms: A framework, Tech. Rep. 1018, School of Operations Research and Industrial Engineering, Cornell University, Ithaca, NY, August 1992.

[17] Y. Ye, On the Q-order of convergence of interior-point algorithms for linear programming, in Proceedings of the 1992 Symposium on Applied Mathematics, W. Fang, ed., Institute of Applied Mathematics, Chinese Academy of Sciences, 1992, pp. 534–543.

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

[19] Y. Ye, O. Güler, R. A. Tapia, and Y. Zhang, A quadratically convergent O(\sqrt{n}L)-iteration algorithm for linear programming, Mathematical Programming, Series A, 59 (1993), pp. 151–162.

[20] Y. Zhang and R. A. Tapia, Superlinear and quadratic convergence of primal–dual interior–point methods for linear programming revisited, Journal of Optimization Theory and Applications, 73 (1992), pp. 229–242.

[21] , A superlinearly convergent polynomial primal–dual interior–point algorithm for linear programming, SIAM Journal on Optimization, 3 (1993), pp. 118–133.

[22] Y. Zhang, R. A. Tapia, and J. E. Dennis, On the superlinear and quadratic convergence of primal-dual interior point linear programming algorithms, SIAM Journal on Optimization, 2 (1992), pp. 304–324.

[23] Y. Zhang, R. A. Tapia, and F. A. Potra, On the superlinear convergence of interior point algorithms for a general class of problems, SIAM Journal on Optimization, 3 (1993), pp. 413–422.