A Globally and Superlinearly Convergent Potential Reduction Interior Point Method for Convex

Programming^{∗}

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

July 15, 1992

Abstract We consider an interior point algorithm for convex programming in which the steps are generated by using a primal-dual affine scaling technique. A “local” variant of the algorithm is shown to have superlinear convergence with q-order up to (but not including) 2. The technique is embedded in a potential reduction algorithm and the resulting method is shown to be globally and superlinearly convergent. An important feature of the convergence analysis is its use of a novel strict interiority condition, which generalizes the usual conical neighborhood of the central path.

1 Introduction

In the past two years, we have seen the appearance of several papers dealing with the construction of primal-dual interior-point algorithms for linear programs (LP) and linear complementarity problems (LCP) which are superlinearly or quadratically convergent. For LP, these works include McShane [7], Mehrotra [8], Tsuchiya [10], Ye [11], Ye et al. [13] and Zhang et al. [14, 16]. For LCP, we cite Ji et al. [3, 4], Kojima et al. [5, 6], Ye and Anstreicher [12] and Zhang et al. [17]. For more details on the historic development of superlinearly or quadratically convergent primal-dual interior-point algorithms, we refer the reader to the introduction section of Ye and Anstreicher [12].

In this paper, we discuss superlinearly convergent primal-dual affine scaling methods for solving the convex programming problem

min_{x} f(x),  Ax = b,  x ≥ 0,  (1)

∗This research was supported by the National Science Foundation (NSF) under Grant No. DDM-9109404 and the Applied Mathematical Sciences subprogram of the Office of Energy Research, U. S. Department of Energy, under Contract W-31-109-Eng-38.

†Department of Systems and Industrial Engineering, University of Arizona, Tucson, Arizona 85721.

‡Argonne National Laboratory, 9700 South Cass Avenue, Argonne, Illinois 60439.

where x ∈ IR^{n}, A ∈ IR^{m×n}, b ∈ IR^{m} and f : IR^{n} → IR is convex and smooth in a sense to be defined below. We assume that the feasible set {x | Ax = b, x ≥ 0} is nonempty and that m < n. The Wolfe dual problem for (1) can be stated as

max_{x,y} f(x) − x^{T} ∇f(x) + b^{T} y,  ∇f(x) − A^{T} y ≥ 0.  (2)

It is well known that if x̄ solves (1), then there is ȳ ∈ IR^{m} such that (x̄, ȳ) solves (2) and, moreover, the optimal values of the primal and dual objective functions are identical. Introducing a slack vector s by

s = ∇f(x) − A^{T} y,

we deduce that the following conditions must be satisfied by the primal and dual solutions:

s̄ = ∇f(x̄) − A^{T} ȳ,  (3a)
Ax̄ = b,  (3b)
x̄ ≥ 0,  (3c)
s̄ ≥ 0,  (3d)
x̄^{T} s̄ = 0.  (3e)

In the subsequent discussion, we say that a point (x, y, s) is “feasible” if it satisfies the equations (3a)–(3d), and “strictly feasible” if (3a)–(3d) are satisfied with x > 0 and s > 0.

Interior point algorithms for the linear and quadratic versions of (1), and the related complementarity problems, have usually required all iterates to belong to a neighborhood of the central path defined by either

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

or

N_{2}(β) = {(x, y, s) feasible | x_{i}s_{i} ≥ β(x^{T} s/n), ∀i = 1, …, n} ,

where

X = diag(x_{1}, x_{2}, …, x_{n}),  S = diag(s_{1}, s_{2}, …, s_{n}),  e = (1, 1, …, 1)^{T} ,

and β is a constant in (0, 1). For example, the predictor-corrector algorithms of Ye et al. [13] for linear programming and Ji, Potra and Huang [3] and Ye and Anstreicher [12] for linear complementarity problems use neighborhoods of the form N_{1}, while the linear programming algorithm of Zhang and Tapia [15] uses N_{2}. In this paper, we use a different neighborhood defined with respect to two parameters δ > 0 and η > 0 in the following way:

N (δ, η) = {(x, y, s) feasible | x^{T} s ≤ 1, x_{i}s_{i} ≥ η(x^{T} s)^{1+δ}, ∀i = 1, …, n} .  (4)

For analytical purposes, this neighborhood has the advantage over N_{1} and N_{2} that local convergence of q-order arbitrarily close to 2 can be obtained without expanding N as the solution is approached. In the papers by Ye [11] and Ye and Anstreicher [12], it is necessary to expand the neighborhood N_{1}(β) during the final stages in order to obtain rapid local convergence.

The basic algorithm we describe in this paper calculates search directions by using a primal-dual affine scaling technique. This technique is equivalent to finding a Newton direction for the system of nonlinear equations formed by (3a), (3b) and the complementarity condition XSe = 0. Because of this connection to Newton’s method, it is immediately clear that such an algorithm would be quadratically convergent if it is started sufficiently close to a nondegenerate solution and if it is allowed to take full steps. In our local analysis, we show that we can ensure that each iterate lies in the neighborhood N (δ, η) without jeopardizing fast local convergence properties. Moreover, we also show that superlinear convergence can be obtained under certain assumptions which allow the possibility of problem (1) to have multiple solutions! Membership of N (δ, η) is enforced by taking a step length α less than one along the affine scaling direction and by “bending” the search path for s in a way to be described in the next section. We show that it is easy to choose an α to ensure that both continued membership of N (δ, η) and fast local convergence apply. By a suitable choice of parameters, the q-order of the local convergence can lie anywhere in the range (1, 2).

To make the algorithm converge globally from any strictly feasible starting point, we embed the affine scaling technique into the potential reduction algorithm discussed in Monteiro [9]. The search direction is calculated as before, but now we choose the step size so that the Tanabe-Todd-Ye potential function

ϕ_{q}(x, s) = q log x^{T} s − ∑_{i=1}^{n} log x_{i}s_{i}.  (5)

(where q > n) is reduced at every iteration according to the Armijo rule, a requirement that is often used in unconstrained optimization algorithms. The actual amount of reduction must be at least a small multiple of the “predicted” reduction obtained from a first order model of ϕ_{q} around the current point. Rules for making an initial guess of the step size and for reducing the step size when it fails the Armijo test are also given.

Finally, we also show that if q is chosen from the range (n, n + 1/4), we can define an algorithm for which both our global and local convergence theories hold. That is, the method is both globally and superlinearly convergent.

The remainder of the paper is laid out as follows: In Section 2, we start by stating our main assumptions. These are, for the most part, standard. The exception is Assumption 5, and the bulk of the section is taken up with proving that this assumption holds in a variety of familiar circumstances, including the cases of nondegeneracy, weak sharp minima, and quadratic objective functions.

In Section 3, we state the linear system of equations that must be solved to obtain the search directions. We then derive bounds on the components of these search directions, in terms of the duality gap x^{T} s and the parameters δ and η that define the neighborhood N. The local convergence theory for our primal-dual affine scaling algorithm appears in Section 4. We show that if the algorithm is started from a point at which x^{T} s is sufficiently small, and which is not too close to the boundary, superlinear convergence of the duality gap to zero can be obtained. Section 5 contains the global convergence theory. We outline the globally convergent algorithm presented in Monteiro [9] and show that if some iterate K satisfies the initial point conditions for the locally convergent algorithm, superlinear convergence also occurs. Finally, we prove that for appropriate choices of q in ϕ_{q}, namely for every q ∈ (n, n + 1/4), these conditions for superlinear convergence are guaranteed to be satisfied by some iterate. As a consequence, we obtain superlinearly convergent algorithms based on the potential function (5). Since the step sizes of the algorithm are selected by means of Armijo rule applied to the potential function (5), there is no reason to believe that the iterates remain in a neighborhood of the form N_{1}(β) or N_{2}(β). The neighborhood N (δ, η) plays a crucial role: we are able to show that the iterates (generated via Armijo rule) remain within a certain N (δ, η) and that one-step superlinear convergence can be obtained.

The following notational conventions are used in the remainder of the paper: Unless otherwise specified, ∥.∥ denotes the the Euclidean norm. We denote the open ball with center x ∈ IR^{n} and radius ϵ > 0 by

B(x, ϵ) ≡ {u ∈ IR^{n} | ∥u − x∥ < ϵ}.

When x and y are two vectors in IR^{n}, the notation [x, y] denotes the set of vectors on the line joining x and y.

For a general vector z ∈ IR^{n} and index set J ⊆ {1, …, n}, z_{J} denotes the vector made up of components z_{i} for i ∈ J. If M ∈ IR^{n×n} and I, J ⊆ {1, …, n} then M_{IJ} denotes the submatrix consisting of elements M_{ij} for which i ∈ I and j ∈ J. The matrix M_{I.} refers to the submatrix M_{IJ} for which J = {1, …, n}. Similarly, M_{.J} denotes the submatrix M_{IJ} in which I = {1, …, m} If D ∈ IR^{n×n} is diagonal, then D_{B} denotes the diagonal matrix whose diagonal entries are D_{ii} for i ∈ B.

We say that (J, L) is a partition of {1, …, n} if J ∪ L = {1, …, n} and J ∩ L = ∅.

Finally, we define some problem-dependent notation: For ϵ ∈ [0, ∞], let

F_{ϵ} =^{△} {(x, y, s) feasible | x^{T} s ≤ ϵ}.

Clearly, F_{0} is just the set of all primal-dual solutions of (1),(2) and F_{∞} is the set of all primal-dual feasible points. We also define a projection operator Π by Π(x, y, z) =^{△} x and, with a slight abuse of notation,

X_{ϵ} =^{△} Π(F_{ϵ})

for ϵ ∈ [0, ∞]. Note that X_{0} is the solution set for (1). The operator P : IR^{n} → IR^{n} defines projection into X_{0}, that is

P (x) =^{△} arg min_{y∈X_{0}} ∥y − x∥.

2 Assumptions

In this section, we state all the assumptions on problem (1) that will be needed in our development. We then show that certain well-studied cases are special cases of these assumptions.

We start by stating all our major assumptions. In this work, we adopt the convention of explicitly stating the assumptions needed in the statement of each result. The main results of the paper will require all the assumptions below to hold.

Assumption 1 (Differentiability and Convexity of f.) The function f is convex over the set {x | Ax = b, x ≥ 0} and twice continuously differentiable in a neighborhood of this set.

Assumption 2 (Local Lipschitz continuity of the Hessian.) For every x ∈ {x | Ax = b, x ≥ 0}, there exists a neighborhood J (x) of x and a positive constant C_{x} such that for all x^{1}, x^{2} ∈ J (x) we have

∥∇^{2}f(x^{1}) − ∇^{2}f(x^{2})∥ ≤ C_{x}∥x^{1} − x^{2}∥.

Assumption 3 The set of strictly feasible solutions of (3a)–(3d) is non-empty and a strictly feasible solution (x^{0}, y^{0}, s^{0}) is given.

Assumption 4 (Existence of a strictly complementary solution.) There exist a partition (B, N) of {1, …, n} and a solution (x^{∗}, y^{∗}, s^{∗}) of (3) such that x_{B}^{∗} > 0 and s_{N}^{∗} > 0.

Assumption 5 There exists a constant L > 0 such that

∥∇f(P (x)) − ∇f(x) − ∇^{2}f(x)(P (x) − x)∥ ≤ L(x^{T} s)^{2},  ∀(x, y, s) ∈ F_{1}.  (6)

Since the left hand side of (6) is uniformly zero when f is quadratic, we can immediately note the following result.

Lemma 2.1 If f is linear or quadratic, then Assumption 5 holds.

In Lemmas 2.6 and 2.7 below, we show that Assumption 5 also holds in other well-studied situations.

Boundedness of the set F_{1}, and therefore of the neighborhood N (δ, η), is an immediate consequence of the following lemma:

Lemma 2.2 Suppose that Assumptions 1 and 3 hold. Then there exists a constant C_{1} > 0 such that if (x, y, s) ∈ F_{1} then

x_{i} ≤ C_{1}  for i = 1, …, n,
s_{i} ≤ C_{1}  for i = 1, …, n.

Proof. Consider the strictly feasible point (x^{0}, y^{0}, s^{0}) as in Assumption 3. By the convexity of f, we can easily show that (x − x^{0})^{T} (s − s^{0}) ≥ 0. Hence

x^{T} s^{0} + x^{0}^{T} s ≤ x^{T} s + x^{0}^{T} s^{0} ≤ 1 + x^{0}^{T} s^{0}.

This implies that

x_{i} ≤ (1 + x^{0}^{T} s^{0})/s_{i}^{0},  (7)
s_{i} ≤ (1 + x^{0}^{T} s^{0})/x_{i}^{0},  (8)

for every i = 1, …, n. The result now follows by setting

C_{1} = \frac{1 + x^{0}^{T} s^{0}}{r^{0}} ,

where r^{0} is the smallest component of (x^{0}, s^{0}).

Lemma 2.3 Suppose that Assumptions 1 and 3 hold and assume that there exist constants ϵ > 0 and L_{ϵ} > 0 such that

∥∇f(P (x)) − ∇f(x) − ∇^{2}f(x)(P (x) − x)∥ ≤ L_{ϵ}(x^{T} s)^{2},  ∀(x, y, s) ∈ F_{ϵ}.  (9)

Then Assumption 5 holds.

Proof. If ϵ ≥ 1 then Assumption 5 obviously holds with L = L_{ϵ}. Assume then that ϵ < 1 and consider the following set

U =^{△} {(x, y, s) feasible | x^{T} s ≤ 1 and x^{T} s ≥ ϵ}.

Using Lemma 2.2, it is easy to see that U is compact. Hence,

a ≡ sup_{(x,y,s)∈U} ∥∇f(P (x)) − ∇f(x) − ∇^{2}f(x)(P (x) − x)∥

is finite. It is now easy to see that Assumption 5 is satisfied with L =^{△} max(L_{ϵ}, a/ϵ).

Lemma 2.4 Suppose that Assumptions 1 and 2 hold. Then for every compact subset K ⊆ {x | Ax = b, x ≥ 0}, there exists a constant C = C(K) > 0 such that

∥∇^{2}f(x^{1}) − ∇^{2}f(x^{2})∥ ≤ C∥x^{1} − x^{2}∥,  ∀x^{1}, x^{2} ∈ K.  (10)

Proof. For x ∈ K, consider the constant C_{x} and the neighborhood J (x) of x as in Assumption 2. Since J (x) is a neighborhood of x, there exists ϵ(x) > 0 such that B(x, 2ϵ(x)) ⊆ J (x). By the compactness of K, there exist x̄^{1}, …, x̄^{p} ∈ K such that

K ⊆ B(x̄^{1}, ϵ(x̄^{1})) ∪ … ∪ B(x̄^{p}, ϵ(x̄^{p})).  (11)

Let ω ≡ sup{∥∇^{2}f(x)∥ | x ∈ K} < ∞, ϵ ≡ min{ϵ(x̄^{1}), …, ϵ(x̄^{p})} and C_{max} ≡ max{C_{}_{x̄^{1}}, …, C_{}_{x̄^{p}}}. We will show that C ≡ max{C_{max}, 2w/ϵ} satisfies the requirement of the lemma. Indeed, let x^{1}, x^{2} ∈ K be given. There are two cases to consider. In the first case, assume that ∥x^{1} − x^{2}∥ < ϵ. By (11), there exists an l ∈ {1, …, p} such that x^{1} ∈ B(x̄^{l}, ϵ(x̄^{l})). Hence, x^{2} ∈ B(x̄^{l}, 2ϵ(x̄^{l})) since ∥x^{1} − x^{2}∥ < ϵ ≤ ϵ(x̄^{l}). Hence, it follows that x^{1}, x^{2} ∈ J (x̄^{l}). By Assumption 2, we have

∥∇^{2}f(x^{1}) − ∇^{2}f(x^{2})∥ ≤ C_{}_{x̄^{l}}∥x^{1} − x^{2}∥ ≤ C∥x^{1} − x^{2}∥.

For the second case, assume that ∥x^{1} − x^{2}∥ ≥ ϵ. Then we have

∥∇^{2}f(x^{1}) − ∇^{2}f(x^{2})∥ ≤ 2w ≤ Cϵ ≤ C∥x^{1} − x^{2}∥.

Hence, the lemma follows.

Using Lemmas 2.3 and 2.4, we can now give a general condition which guarantees the validity of Assumption 5.

Lemma 2.5 Suppose that Assumptions 1, 2 and 3 hold. Assume that there exists constants ϵ > 0 and υ > 0 such that the following implication holds:

(x, y, s) ∈ F_{ϵ} ⟹ ∥x − P (x)∥ ≤ υ(x^{T} s).  (12)

Then Assumption 5 holds.

Proof. We may assume without loss of generality that ϵ ≤ 1. Lemma 2.2 implies that X_{1} is bounded, and hence that the sets X_{0} and X_{ϵ} are also bounded. Since these two sets are also closed, they are compact. Hence, the set defined as

K ≡ ⋃ {[x, x̄] | x ∈ X_{ϵ} and x̄ ∈ X_{0}}

is also compact. By lemma 2.4, there exists C > 0 such that (10) holds. Hence, for every (x, y, s) ∈ F_{ϵ}, we obtain

∥∇f(P (x)) − ∇f(x) − ∇^{2}f(x)(P (x) − x)∥
≤ max_{x̃∈[x,P (x)]} ∥[∇^{2}f(x̃) − ∇^{2}f(x)] [P (x) − x]∥
≤ C∥x̃ − x∥ ∥P (x) − x∥
≤ C∥P (x) − x∥^{2}
≤ Cυ^{2}(x^{T} s)^{2},

where the second inequality is due to the fact that both x and x̃ are in K and the last inequality is due to implication (12). Therefore, the assumptions of Lemma 2.3 are satisfied and, as a consequence, Assumption 5 holds.

Two important situations in which the implication (12) holds are stated in the next two lemmas.

Lemma 2.6 Suppose that Assumptions 1, 2 and 3 hold. Then Assumption 5 holds if problem (1) has a weak sharp minimum in the sense of Burke and Ferris [1], that is, if there exists a constant λ > 0 such that

∥x − P (x)∥ ≤ λ [f(x) − f(P (x))] ,

for all x ∈ {x | Ax = b, x ≥ 0}.

Proof. For every feasible (x, y, s), we have

∥x − P (x)∥ ≤ λ [f(x) − f(P (x))] ≤ λ(x^{T} s)

which shows that the implication (12) holds with ϵ = ∞ and υ = λ.

The other important case which guarantees the validity of implication (12), and hence of Assumption 5, requires that the solution of problem (1) be nondegenerate. We now make precise the notion of nondegeneracy of a solution of problem (1).

Definition 1 A solution (x̄, ȳ, s̄) of problem (1) is called nondegenerate if the following conditions are satisfied:

(a) (x̄, ȳ, s̄) is a strictly complementary solution, that is, there exists a partition (B, N) of {1, …, n} such that x̄_{B} > 0 and s̄_{N} > 0;

(b) The matrix Ā defined as

Ā ≡ [ \begin{matrix}A_{B} & A_{N} \\ 0 & I \end{matrix} ]

has linearly independent columns;

(c) ∇^{2}f(x̄) is nonsingular in the null space of Ā.

We observe that the rows of Ā are the gradients of the constraints of problem (1) which are active at the solution x̄. Clearly, a nondegenerate solution must be unique. We are now ready to state the second case under which Assumption 5 holds.

Lemma 2.7 Suppose that Assumptions 1, 2 and 3 hold and assume that problem (1) has the unique nondegenerate solution (x̄, ȳ, s̄). Then Assumption 5 holds.

Proof. Consider the function F : IR^{n} × IR^{m} × IR^{n} → IR^{n} × IR^{m} × IR^{n} defined by

F (x, y, s) = [ \begin{matrix}XSe \\ Ax − b \\ −∇f(x) + A^{T} y + s \end{matrix} ] .  (13)

Using the assumption that (x̄, ȳ, s̄) is a nondegenerate solution, one can easily show that ∇F (x̄, ȳ, s̄) is nonsingular. To simplify notation, let w = (x, y, s) and w̄ = (x̄, ȳ, s̄). The differentiability of F at w̄ implies that we can write

F (w) = F (w̄) + ∇F (w̄)(w − w̄) + ∥w − w̄∥r(w − w̄),  (14)

where lim_{∥h∥→0} r(h) = 0. Define M =^{△} (∥∇F (w̄)∥)^{−1}. Then there exists ξ > 0 such that ∥r(w − w̄)∥ ≤ M/2 for every w ∈ B(w̄, ξ). Hence, for every w ∈ B(w̄, ξ) we obtain using (14) that

∥F (w) − F (w̄)∥ ≥ ∥∇F (w̄)(w − w̄)∥ − ∥r(w − w̄)∥ ∥w − w̄∥
≥ M∥w − w̄∥ − ∥r(w − w̄)∥ ∥w − w̄∥  (15)
≥ \frac{M}{2} ∥w − w̄∥,

where the last inequality is due to the fact that ∥r(w − w̄)∥ ≤ M/2. Using the definition (13) of F and choosing w = (x, y, s) to be feasible, we obtain

∥F (w) − F (w̄)∥ = ∥F (x, y, s) − F (x̄, ȳ, s̄)∥ = ∥XSe∥ ≤ x^{T} s.  (16)

Using (15) and (16), we obtain

∥x − x̄∥ ≤ ∥w − w̄∥ ≤ \frac{2}{M} ∥F (w) − F (w̄)∥ ≤ \frac{2}{M} x^{T} s,  (17)

for every w ∈ B(w̄, ξ). If F_{1} ⊆ B(w̄, ξ), then the assumptions of Lemma 2.5 hold with ϵ = 1 and υ = 2/M, so we are done. Otherwise, consider the set F_{1}\B(w̄, ξ). Clearly, x^{T} s > 0 for all (x, y, s) ∈ F_{1}\B(w̄, ξ). Because of Lemma 2.2, the set F_{1}\B(w̄, ξ) is compact, and hence there exists ϵ > 0 such that

x^{T} s > ϵ, ∀(x, y, s) ∈ F_{1}\B(w̄, ξ).

This implies that F_{ϵ} ⊆ B(w̄, ξ), and hence that implication (12) holds with υ = 2/M. By Lemma 2.5, Assumption 5 must hold.

3 Technical Results

In this section, we briefly review how the primal-dual affine scaling search direction is computed from a given strictly feasible point (x, y, s). The remainder of the section is taken up with finding bounds on the components of the search direction vectors for points (x, y, s) lying in the neighborhood N (δ, η), where δ and η are fixed a priori.

We start by describing how the primal-dual affine scaling search direction is computed. Given a strictly feasible point (x, y, s), we obtain the primal-dual affine scaling search direction (Δx, Δy, Δs) by applying one step of Newton’s method to the nonlinear system defined by (3a), (3b) and the equation XSe = 0. Hence (Δx, Δy, Δs) can be obtained by solving the linear system

SΔx + XΔs = −SXe,  (18a)
AΔx = 0,  (18b)
−∇^{2}f(x)Δx + A^{T} Δy + Δs = 0.  (18c)

The remaining of the section is devoted to finding bounds on the components of the search direction (Δx, Δy, Δs). The main result of this section is stated in Theorem 3.8 which in turn is obtained by combining the results of Lemmas 3.3 and 3.7. We start by stating a simple result that is used in a number of places in this section.

Lemma 3.1 Suppose Assumption 1 holds. Then if (Δx, Δy, Δs) satisfies relations (18b) and (18c) then Δx^{T} Δs ≥ 0.

Proof. The proof follows immediately from the convexity of f(x).

We can now state some simple results concerning bounds on certain components of x, Δx, s, and Δs.

Lemma 3.2 Under Assumptions 1 and 4, if (x, y, s) is feasible then there exists a constant r > 0 such that

x_{i} ≤ (x^{T} s)/r,  ∀i ∈ N,  (19)
s_{i} ≤ (x^{T} s)/r,  ∀i ∈ B.  (20)

Proof. Consider the strictly complementary solution (x^{∗}, y^{∗}, s^{∗}) as in Assumption 4. Using convexity of f, we obtain

0 ≤ (x − x^{∗})^{T} (s − s^{∗}) = x^{T} s − x^{∗}^{T} s − x^{T} s^{∗}

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

s_{i} ≤ (x^{T} s)/x_{i}^{∗}, ∀i ∈ B.

Similarly, we can show that

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

The result now follows by setting r to be the smallest component of (x_{B}^{∗}, s_{N}^{∗}).

Lemma 3.3 Suppose that Assumptions 1 and 4 hold and, for some δ > 0 and η ∈ (0, 1], let (x, y, s) ∈ N (δ, η) be given. Let (Δx, Δy, Δs) denote the solution of (18). Then there exists r > 0 such that

|Δx_{i}| ≤ \frac{(x^{T} s)^{1−(1/2)δ}}{rη^{1/2}} ,  ∀i ∈ N,  (21)
|Δs_{i}| ≤ \frac{(x^{T} s)^{1−(1/2)δ}}{rη^{1/2}} ,  ∀i ∈ B,  (22)
s_{i} ≥ rη(x^{T} s)^{δ},  ∀i ∈ N,  (23)
x_{i} ≥ rη(x^{T} s)^{δ},  ∀i ∈ B.  (24)

Proof. Defining D = X^{1/2}S^{−1/2} and using relation (18a), we obtain

D^{−1}Δx + DΔs = −(SX)^{1/2}e.

Taking the square of the norm of both sides of the above relation, we obtain

∥D^{−1}Δx∥^{2} + ∥DΔs∥^{2} + 2Δx^{T} Δs = x^{T} s,

and since by Lemma 3.1 Δx^{T} Δs ≥ 0, this relation implies

∥D^{−1}Δx∥ ≤ (x^{T} s)^{1/2},  (25)
∥DΔs∥ ≤ (x^{T} s)^{1/2}.  (26)

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

x_{i}s_{i} ≥ η(x^{T} s)^{1+δ},  ∀i = 1, …, n.  (27)

Using relations (25), (19) and (27), we obtain

|Δx_{i}| ≤ x_{i} (\frac{x^{T} s}{x_{i}s_{i}})^{1/2} ≤ (\frac{x^{T} s}{r}) (\frac{1}{η(x^{T} s)^{δ}} )^{1/2} = \frac{(x^{T} s)^{1−(δ/2)}}{rη^{1/2}} ,

for every i ∈ N. This yields (21). Relation (22) follows similarly with the aid of (26), (20) and (27). To show (23), observe that relations (27) and (19) imply

s_{i} ≥ \frac{η(x^{T} s)^{1+δ}}{x_{i}} ≥ rη(x^{T} s)^{δ},

for every i ∈ N. Hence (23) follows. Similarly, (24) follows with the aid of (27) and (20).

Providing upper bounds for the remaining components of the search directions, namely Δx_{B} and Δs_{N}, is more difficult. This part of the development is based on the approach which appears in the paper of Ye and Anstreicher [12].

We start with the following lemma of Ye and Anstreicher [12, Lemma 3.4].

Lemma 3.4 Let M ∈ IR^{p×p} be a positive semi-definite matrix and assume that (J, L) form a partition of {1, …, p}. Then,

R ( \begin{matrix}M_{JJ} & M_{JL} \\ 0 & I \end{matrix} ) = R ( \begin{matrix}M_{JJ}^{T} & M_{LJ}^{T} \\ 0 & −I \end{matrix} ) .  (28)

As a consequence of the above result, we obtain

Lemma 3.5 Let Q ∈ IR^{n×n} be a positive semi-definite matrix and let A ∈ IR^{m×n} be arbitrary. Let (B, N) form an arbitrary partition of {1, …, n}. Define

U_{1} ≡ ( \begin{matrix}Q_{BB} & −A_{.B}^{T} & Q_{BN} \\ A_{.B} & 0 & A_{.N} \\ 0 & 0 & I \end{matrix} ) ,  U_{2} ≡ ( \begin{matrix}Q_{BB}^{T} & A_{.B}^{T} & Q_{NB}^{T} \\ −A_{.B} & 0 & −A_{.N} \\ 0 & 0 & −I \end{matrix} ) .  (29)

Then, R(U_{1}) = R(U_{2}).

Proof. First observe that the matrix

M ≡ ( \begin{matrix}Q_{BB} & −A_{.B}^{T} & Q_{BN} \\ A_{.B} & 0 & A_{.N} \\ Q_{NB} & −A_{.N}^{T} & Q_{NN} \end{matrix} )  (30)

is positive semi-definite since it can be obtained by a symmetric permutation of the positive semi-definite matrix

( \begin{matrix}Q & −A^{T} \\ A & 0 \end{matrix} ) .  (31)

The lemma now follows from Lemma 3.4 applied to the matrix M and the partition (J, L) of {1, …, m + n} for which

M_{JJ} = ( \begin{matrix}Q_{BB} & −A_{.B}^{T} \\ A_{.B} & 0 \end{matrix} ) .  (32)

Using the above result, we next prove a lemma that is similar to Lemma 3.5 of Ye and Anstreicher [12]. In fact, the two results are identical if the function f(x) is quadratic.

Lemma 3.6 Suppose that Assumption 4 holds. Let (x, y, s) be a strictly feasible point and let D ≡ X^{1/2}S^{−1/2}. Let (Δx, Δy, Δs) denote the solution of (18). Then (u, v, w) = (Δx_{B}, Δs_{N}, Δy) solves the problem

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

subject to

A_{.B}u = −A_{.N}Δx_{N},
−Q_{BB}u + A_{.B}^{T}w = Q_{BN}Δx_{N} − Δs_{B},  (34)
−Q_{NB}u + A_{.N}^{T}w + v = Q_{NN}Δx_{N},

where

Q = ∇^{2}f(x), g = −∇f(x̄) + ∇f(x) + ∇^{2}f(x)(x̄ − x),

and x̄ is any primal solution of (1).

Proof. ¿From relations (18b) and (18c) it follows immediately that (u, v, w) = (Δx_{B}, Δs_{N}, Δy) is feasible with respect to the constraints (34). The result follows once we verify that (u, v, w) = (Δx_{B}, Δs_{N}, Δy) satisfies the KKT conditions for the above problem, namely:

( \begin{matrix}g_{B} + D_{B}^{−2}Δx_{B} \\ 0 \\ D_{N}^{2}Δs_{N} \end{matrix} ) ∈ R ( \begin{matrix}A_{.B}^{T} & −Q_{BB}^{T} & −Q_{NB}^{T} \\ 0 & A_{.B} & A_{.N} \\ 0 & 0 & I \end{matrix} ) = R(U_{2}),  (35)

where U_{2} is the matrix defined in Lemma 3.5. Indeed, using equation (18a), it is easily seen that

D^{−2}Δx = −(s + Δs),  (36)
D^{2}Δs = −(x + Δx).  (37)

Let (ȳ, z̄) denote a pair which together with the solution x̄ satisfies (3). Using the fact that both (x, y, s) and (x̄, ȳ, z̄) satisfy (3a)–(3b) and that (Δx, Δy, Δs) satisfy relations (18b) and (18c), one can easily show that

A(x + Δx − x̄) = 0,  (38)
−Q(x + Δx − x̄) + A^{T} (y + Δy − ȳ) + (s + Δs − s̄) = g.  (39)

Since x̄_{N} = 0, s̄_{B} = 0 and using relations (36) and (39), we obtain

g_{B} + D_{B}^{−2}Δx_{B} = g_{B} − (s_{B} + Δs_{B})
= g_{B} − (s_{B} + Δs_{B} − s̄_{B})
= −Q_{B.}(x + Δx − x̄) + A_{.B}^{T}(y + Δy − ȳ)  (40)
= −Q_{BB}(x_{B} + Δx_{B} − x̄_{B}) − Q_{BN}(x_{N} + Δx_{N})
+ A_{.B}^{T}(y + Δy − ȳ).

Also, from (37) and (38) it follows that

D_{N}^{2}Δs_{N} = −(x_{N} + Δx_{N}),  (41)
0 = −A_{.B}(x_{B} + Δx_{B} − x̄_{B}) − A_{.N}(x_{N} + Δx_{N}).  (42)

Using relations (40), (41) and (42), we obtain

( \begin{matrix}g_{B} + D_{B}^{−2}Δx_{B} \\ 0 \\ D_{N}^{2}Δs_{N} \end{matrix} ) ∈ R(U_{1}).  (43)

where U_{1} is the matrix defined in Lemma 3.5. Relation (35) now follows from the fact that R(U_{1}) = R(U_{2}) in view of Lemma 3.5.

Lemma 3.7 Suppose that Assumptions 1, 3, 4, and 5 hold, and that (x, y, s) ∈ N (δ, η) for some δ > 0 and η ∈ (0, 1]. Then there exists a constant C_{2} > 0 independent of δ and η such that

∥Δx_{B}∥ ≤ \frac{C_{2}}{η^{3/2}} (x^{T} s)^{1−(3/2)δ},  (44)
∥Δs_{N}∥ ≤ \frac{C_{2}}{η^{3/2}} (x^{T} s)^{1−(3/2)δ}.  (45)

Proof. It is well known that since (34) is a consistent set of equalities, there is a triple (ū, v̄, w̄) which satisfies (34) and is bounded above by a constant times the norm of the right hand side of (34). Therefore,

∥(ū, v̄, w̄)∥ ≤ C_{3}∥(Δx_{N}, Δs_{B})∥,  (46)

for some constant C_{3} > 0 independent of δ and η. Define K_{1} ≡ max(∥D_{B}∥, ∥D_{N}^{−1}∥). Then,

\begin{matrix}∥Δx_{B} + D_{B}^{2}g_{B}∥^{2} + ∥Δs_{N}∥^{2} \\ ≤ ∥D_{B}∥^{2} ∥D_{B}^{−1}Δx_{B} + D_{B}g_{B}∥^{2} + ∥D_{N}^{−1}∥^{2} ∥D_{N}Δs_{N}∥^{2} \\ ≤ K_{1}^{2} [ ∥D_{B}^{−1}Δx_{B} + D_{B}g_{B}∥^{2} + ∥D_{N}Δs_{N}∥^{2}] \\ ≤ K_{1}^{2} [ ∥D_{B}^{−1}ū + D_{B}g_{B}∥^{2} + ∥D_{N} v̄∥^{2}] , \end{matrix}  (47)

where the third inequality follows from the fact that (Δx_{B}, Δs_{N}, Δy) is an optimal solution for problem (33)-(34). Relation (47) can be rewritten as

∥( \begin{matrix}Δx_{B} + D_{B}^{2}g_{B} \\ Δs_{N} \end{matrix} )∥ ≤ K_{1} ∥( \begin{matrix}D_{B}^{−1}ū + D_{B}g_{B} \\ D_{N} v̄ \end{matrix} )∥ .

Define K_{2} ≡ max(∥D_{B}^{−1}∥, ∥D_{N}∥). Then, using the triangle inequality twice in the above relation and relations (46), (21) and (22), we obtain

∥(\begin{matrix}Δx_{B} \\ Δs_{N} \end{matrix})∥ ≤ K_{1} \{∥[ \begin{matrix}D_{B}^{−1} & 0 \\ 0 & D_{N} \end{matrix} ]∥ ∥(ū ; v̄)∥ + ∥(\begin{matrix}D_{B}g_{B} \\ 0 \end{matrix})∥\} + ∥(\begin{matrix}D_{B}^{2}g_{B} \\ 0 \end{matrix})∥

≤ K_{1} \{K_{2}C_{3} ∥(\begin{matrix}Δx_{N} \\ Δs_{B} \end{matrix})∥ + ∥D_{B}∥ ∥g_{B}∥\} + ∥D_{B}∥^{2}∥g_{B}∥  (48)

≤ K_{1}K_{2}C_{3}\sqrt{n} \frac{(x^{T} s)^{1−(δ/2)}}{rη^{1/2}} + 2K_{1}^{2}∥g_{B}∥.

It remains to estimate K_{1}, K_{2} and ∥g_{B}∥ in the above relation. First, we estimate K_{1}. Observe that since (x, y, s) ∈ N (δ, η), relation (27) holds. Using (27), the definition of D and Lemma 2.2, we obtain

∥D_{B}^{2}∥ = ∥X_{B}S_{B}^{−1}∥ = max_{i∈B} \{\frac{x_{i}^{2}}{x_{i}s_{i}} \} ≤ \frac{C_{1}^{2}}{η} (x^{T} s)^{−(1+δ)}.

Similarly, we have

∥D_{N}^{−2}∥ = ∥X_{N}^{−1}S_{N}∥ = max_{i∈N} \{\frac{s_{i}^{2}}{x_{i}s_{i}} \} ≤ \frac{C_{1}^{2}}{η} (x^{T} s)^{−(1+δ)}.

Hence,

K_{1}^{2} ≤ \frac{C_{1}^{2}}{η} (x^{T} s)^{−(1+δ)}.  (49)

We now estimate K_{2}. Using (27), the definition of D and Lemma 3.2, we obtain

∥D_{B}^{−2}∥ = ∥X_{B}^{−1}S_{B}∥ = max_{i∈B} \{\frac{s_{i}^{2}}{x_{i}s_{i}} \} ≤ \frac{(x^{T} s/r)^{2}}{η(x^{T} s)^{1+δ}} = \frac{(x^{T} s)^{1−δ}}{ηr^{2}} .

Similarly, we have

∥D_{N}^{2}∥ = ∥X_{N}S_{N}^{−1}∥ = max_{i∈N} \{\frac{x_{i}^{2}}{x_{i}s_{i}} \} ≤ \frac{(x^{T} s/r)^{2}}{η(x^{T} s)^{1+δ}} = \frac{(x^{T} s)^{1−δ}}{ηr^{2}} .

Hence,

K_{2}^{2} ≤ \frac{(x^{T} s)^{1−δ}}{ηr^{2}} .  (50)

Next, we estimate ∥g_{B}∥. We can assume that the primal solution x̄ of Lemma 3.6 is equal to P (x). Using x̄ = P (x) in the definition of g and Assumption 5, we obtain

∥g_{B}∥ ≤ ∥g∥ = ∥ − ∇f(P (x)) + ∇f(x) + ∇^{2}f(x)(P (x) − x)∥ ≤ L(x^{T} s)^{2}.  (51)

Substituting the estimates (49), (50) and (51) into expression (48), we obtain

∥(\begin{matrix}Δx_{B} \\ Δs_{N} \end{matrix})∥ ≤ (\frac{C_{1}^{2}}{η} (x^{T} s)^{−(1+δ)})^{1/2} (\frac{(x^{T} s)^{1−δ}}{ηr^{2}} )^{1/2} \sqrt{n} C_{3} \frac{(x^{T} s)^{1−(1/2)δ}}{rη^{1/2}}

+ 2 (\frac{C_{1}^{2}}{η} (x^{T} s)^{−(1+δ)}) (L(x^{T} s)^{2})

= \frac{(x^{T} s)^{1−(3/2)δ}}{η^{3/2}} [\frac{\sqrt{n} C_{1}C_{3}}{r^{2}} + 2C_{1}^{2}Lη^{1/2}(x^{T} s)^{δ/2}] .

Since x^{T} s ≤ 1 and η ∈ (0, 1], we can bound the bracketed term by a constant C_{2} independent of η and δ. We thus obtain

∥(\begin{matrix}Δx_{B} \\ Δs_{N} \end{matrix})∥ ≤ \frac{C_{2}}{η^{3/2}} (x^{T} s)^{1−(3/2)δ},

which is the desired result.

We can now merge the results of Lemmas 3.3 and 3.7 to obtain the following theorem.

Theorem 3.8 Suppose that Assumptions 1, 3, 4, and 5 hold, and that (x, y, s) ∈ N (δ, η) for some δ > 0 and η ∈ (0, 1]. Then there exists a constant C_{4} independent of δ and η such that

∥Δx∥ ≤ \frac{C_{4}}{η^{3/2}} (x^{T} s)^{1−(3/2)δ},
∥Δs∥ ≤ \frac{C_{4}}{η^{3/2}} (x^{T} s)^{1−(3/2)δ}.

4 A locally convergent algorithm

In this section, we discuss a basic primal-dual affine scaling algorithm based on the step calculation (18) and show that rapid local convergence can be obtained for certain choices of the step lengths. In the next section, this algorithm is embedded in the globally convergent potential reduction algorithm of Monteiro [9].

In describing the algorithm, we first show how the step calculated in (18) can be used to move from one iterate to the next. In general, we cannot simply move along the direction (Δx, Δy, Δs) from the current point (x, y, s), since this may lead to a new point that violates the condition (3a). Instead, we take a “curved” step in the s component and make the transition

(x, y, s) → (x + αΔx, y + αΔy, s(α)),

where

s(α) = ∇f(x + αΔx) − A^{T} (y + αΔy).  (52)

The algorithm proceeds as follows:

Algorithm PDA

initially: choose (x^{0}, y^{0}, s^{0}) strictly feasible

for k = 1, 2, …

Find (Δx^{k}, Δy^{k}, Δs^{k}) by solving (18) with (x, y, s) = (x^{k}, y^{k}, s^{k});
Choose α_{k} > 0 such that (x^{k} + α_{k}Δx^{k}, y^{k} + α_{k}Δy^{k}, s^{k}(α_{k}))

is strictly feasible;

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

end for

The following neighborhood will play an important role in the subsequent development. Let p > 0, δ > 0, η ∈ (0, 1] and γ > 0 be given and define

N_{p}(δ, η, γ) ≡ {(x, y, s) ∈ N (δ, η) | (x^{T} s)^{1−pδ}/(δη^{p}) ≤ 1/γ}.  (53)

The proof of the following lemma is straightforward and is left to the reader.

Lemma 4.1 Assume that p > q > 0, δ > 0 and η ∈ (0, 1]. Then

(x^{T} s)^{1−qδ}/η^{q} ≤ (x^{T} s)^{1−pδ}/η^{p},  ∀x ∈ N (δ, η),  (54)

and

N_{p}(δ, η, γ) ⊆ N_{q}(δ, η, γ),  ∀γ > 0.  (55)

Lemma 4.2 Suppose that all Assumptions 1–5 hold and let (x, y, s) ∈ N_{4}(δ, η, γ) be given, where δ > 0, η ∈ (0, 1], and γ ≥ δ. Then, there exist a constant C_{5} independent of δ, η and γ such that for all α ∈ [0, 1],

|(x_{i} + αΔx_{i})s_{i}(α) − (1 − α)x_{i}s_{i}| ≤ C_{5}α^{2} \frac{(x^{T} s)^{2−3δ}}{η^{3}} ,  ∀i = 1, …, n,  (56)

and

|(x + αΔx)^{T} s(α) − (1 − α)x^{T} s| ≤ nC_{5}α^{2} \frac{(x^{T} s)^{2−3δ}}{η^{3}} .  (57)

Proof. We first note that (57) follows immediately from (56). Using relations (52), (3a) and (18c), we obtain

s(α) = ∇f(x + αΔx) − A^{T} y − αA^{T} Δy
= ∇f(x + αΔx) − (∇f(x) − s) − α(∇^{2}f(x)Δx − Δs)
= s + αΔs + [∇f(x + αΔx) − ∇f(x) − α∇^{2}f(x)Δx]  (58)
= s + αΔs + α ∫_{0}^{1}[∇^{2}f(x + tαΔx) − ∇^{2}f(x)]Δx dt
= s + αΔs + αR(x, αΔx)Δx

where

R(x, Δx) =^{△} ∫_{0}^{1}[∇^{2}f(x + tΔx) − ∇^{2}f(x)] dt.  (59)

Relations (58) and (18a) imply

(x_{i} + αΔx_{i})s_{i}(α) = (x_{i} + αΔx_{i}) [s_{i} + αΔs_{i} + α {R(x, αΔx)Δx}_{i}]
= (1 − α)x_{i}s_{i} + α^{2}Δx_{i}Δs_{i} + α(x_{i} + αΔx_{i}) {R(x, αΔx)Δx}_{i} ,

so that

\begin{matrix}|(x_{i} + αΔx_{i})s_{i}(α) − (1 − α)x_{i}s_{i}| \\ ≤ α^{2}|Δx_{i}Δs_{i}| + α|(x_{i} + αΔx_{i})| | \{R(x, αΔx)Δx\}_{i} |. \end{matrix}  (60)

We will now estimate each of the quantities which appears in the right hand side of (60). Regardless of whether i ∈ B or i ∈ N, we have from Lemmas 3.3 and 3.7 that

|Δx_{i}Δs_{i}| ≤ \frac{(x^{T} s)^{1−(1/2)δ}}{rη^{1/2}} \frac{C_{2}}{η^{3/2}} (x^{T} s)^{1−(3/2)δ} ≤ \frac{C_{6}}{η^{2}} (x^{T} s)^{2−2δ},  (61)

where

C_{6} ≡ C_{2}/r.  (62)

Relation (55) with p = 4 and q = 3/2 and the assumption that (x, y, s) ∈ N_{4}(δ, η, γ) imply that (x, y, s) ∈ N_{3/2}(δ, η, γ). Hence, (x^{T} s)^{1−(3δ/2)}/η^{3/2} ≤ δ/γ ≤ 1 in view of (53) and the assumption that γ ≥ δ. Using this inequality, Lemma 2.2 and Theorem 3.8, we obtain

|x_{i} + αΔx_{i}| ≤ |x_{i}| + |Δx_{i}| ≤ C_{1} + C_{4} \frac{(x^{T} s)^{1−(3δ/2)}}{η^{3/2}} ≤ C_{1} + C_{4}.  (63)

Consider now the set

K =^{△} {x ∈ IR^{n} | Ax = b, 0 ≤ x_{i} ≤ C_{1} + C_{4}, ∀i = 1, …, n}.

Clearly, K is a compact subset of {x | Ax = b, x ≥ 0}. Hence, by Lemma 2.4 there exists a constant C > 0 such that

∥∇^{2}f(x^{1}) − ∇^{2}f(x^{2})∥ ≤ C∥x^{1} − x^{2}∥, ∀x^{1}, x^{2} ∈ K.

Using this last relation and Theorem 3.8, we obtain

\begin{aligned}| \{R(x, αΔx)Δx\}_{i} | & ≤ & ∥R(x, αΔx)∥ ∥Δx∥ \\ & ≤ & (∫_{0}^{1} ∥∇^{2}f(x + tαΔx) − ∇^{2}f(x)∥ dt) ∥Δx∥ \\ & ≤ & α\frac{C}{2} ∥Δx∥^{2} \\ & ≤ & α\frac{CC_{4}^{2}}{2} \frac{(x^{T} s)^{2−3δ}}{η^{3}} . \end{aligned}  (64)

Substituting relations (61), (63) and (64) into (60), we obtain

\begin{matrix}|(x_{i} + αΔx_{i})s_{i}(α) − (1 − α)x_{i}s_{i}| \\ ≤ α^{2} \frac{C_{6}}{η^{2}} (x^{T} s)^{2−2δ} + α (C_{1} + C_{4}) (α\frac{CC_{4}^{2}}{2} \frac{(x^{T} s)^{2−3δ}}{η^{3}} ) \\ ≤ α^{2} \frac{(x^{T} s)^{2−3δ}}{η^{3}} \{C_{6}η(x^{T} s)^{δ} + (C_{1} + C_{4})(CC_{4}^{2}/2)\} \\ ≤ α^{2} \frac{(x^{T} s)^{2−3δ}}{η^{3}} \{C_{6} + (C_{1} + C_{4})(CC_{4}^{2}/2)\} . \end{matrix}  (65)

Relation (56) now follows if we let

C_{5} ≡ C_{6} + (C_{1} + C_{4})(CC_{4}^{2}/2).  (66)

In Lemma 4.3, we show that the property of belonging to the set N_{4}(δ, η, γ) is “inheritable” when γ is sufficiently large; that is, we can take a substantial step with parameter α_{k} from the current iterate (x^{k}, y^{k}, s^{k}) and remain within N_{4}(δ, η, γ).

Lemma 4.3 Suppose that δ ∈ (0, 1/4) and η ∈ (0, 1] are given and that Assumptions 1–5 hold. Then there exists a constant γ̄ ≥ 1 independent of δ and η that satisfies the following implication. If γ ≥ γ̄ and (x, y, s) ∈ N_{4}(δ, η, γ) then (x + αΔx, y + αΔy, s(α)) ∈ N_{4}(δ, η, γ) for all

α ∈ [0, 1 − \frac{γ}{δη^{4}} (x^{T} s)^{1−4δ}] .  (67)

Moreover, we have that

(x + αΔx)^{T} s(α) ≤ (1 − α (1 − \frac{δ}{2})) x^{T} s,  ∀α ∈ [0, 1].  (68)

Proof. Consider the constant C_{5} as in the statement of Lemma 4.2 and define

γ̄ ≡ max{1, (2n + 1)C_{5})}.  (69)

We will show that γ̄ fulfills the requirements of the lemma. Indeed, let γ ≥ γ̄ and (x, y, s) ∈ N_{4}(δ, η, γ) be given. Note that (69) implies that γ ≥ 1. Our first task is to show that the range of α values such that

(x_{i} + αΔx_{i})s_{i}(α) ≥ η [(x + αΔx)^{T} s(α)]^{1+δ} ,  (70)

for all i = 1, 2, …, n, includes the interval in (67). By relations (56) and (27), we have

(x_{i} + αΔx_{i})s_{i}(α) ≥ (1 − α)x_{i}s_{i} − C_{5}α^{2} \frac{(x^{T} s)^{2−3δ}}{η^{3}}
≥ (1 − α)η(x^{T} s)^{1+δ} − C_{5}α^{2} \frac{(x^{T} s)^{2−3δ}}{η^{3}}  (71)
≥ η(x^{T} s)^{1+δ} [1 − α (1 + C_{5} \frac{(x^{T} s)^{1−4δ}}{η^{4}} )] .

for all α ∈ [0, 1]. Similarly, from relation (57), we obtain

(x + αΔx)^{T} s(α) ≤ x^{T} s [1 − α (1 − nC_{5} \frac{(x^{T} s)^{1−3δ}}{η^{3}} )] ,  ∀α ∈ [0, 1].  (72)

Relation (55) with p = 4 and q = 3 and the assumption that (x, y, s) ∈ N_{4}(δ, η, γ) imply that (x, y, s) ∈ N_{3}(δ, η, γ). Hence, using (53) we obtain

(x^{T} s)^{1−3δ}/η^{3} ≤ δ/γ.  (73)

We now use the fact that for all δ ∈ (0, 1) and u ≥ −1, (1 + u)^{δ} ≤ 1 + δu. Using this inequality, (72), (73) and the fact that nC_{5}δ/γ ≤ 1, we find that

\begin{matrix}η [(x + αΔx)^{T} s(α)]^{1+δ} \\ ≤ η(x^{T} s)^{1+δ} [1 − α (1 − nC_{5} \frac{(x^{T} s)^{1−3δ}}{η^{3}} )] [1 − δα (1 − nC_{5} \frac{(x^{T} s)^{1−3δ}}{η^{3}} )] \\ ≤ η(x^{T} s)^{1+δ} [1 − α(1 + δ) (1 − nC_{5} \frac{(x^{T} s)^{1−3δ}}{η^{3}} ) + δα^{2} (1 − nC_{5} \frac{(x^{T} s)^{1−3δ}}{η^{3}} )^{2}] \\ ≤ η(x^{T} s)^{1+δ} [1 − α(1 + δ) (1 − nC_{5} \frac{(x^{T} s)^{1−3δ}}{η^{3}} ) + δα^{2}] . \end{matrix}  (74)

The inequality (70) will hold if the right-hand side of (74) is smaller than or equal to the right-hand side of (71). In other words, we want α to be chosen so that

η(x^{T} s)^{1+δ} [1 − α(1 + δ) (1 − nC_{5} \frac{(x^{T} s)^{1−3δ}}{η^{3}} ) + δα^{2}]
≤ η(x^{T} s)^{1+δ} [1 − α (1 + C_{5} \frac{(x^{T} s)^{1−4δ}}{η^{4}} )] .

This is equivalent to

(1 + δ) (1 − nC_{5} \frac{(x^{T} s)^{1−3δ}}{η^{3}} ) − δα ≥ (1 + C_{5} \frac{(x^{T} s)^{1−4δ}}{η^{4}} ) ,

that is,

α ≤ \frac{1}{δ} \{(1 + δ) [1 − nC_{5} \frac{(x^{T} s)^{1−3δ}}{η^{3}} ] − [1 + C_{5} \frac{(x^{T} s)^{1−4δ}}{η^{4}} ]\}  (75)

To show that (x + αΔx, y + αΔy, s(α)) ∈ N (δ, η) for every α satisfying (67), it is sufficient to show that the right hand side of (75) is greater than or equal to 1 − (γ/δη^{4})(x^{T} s)^{1−4δ}. Indeed, we have

\frac{1}{δ} \{(1 + δ) [1 − nC_{5} \frac{(x^{T} s)^{1−3δ}}{η^{3}} ] − [1 + C_{5} \frac{(x^{T} s)^{1−4δ}}{η^{4}} ]\}
≥ \frac{1}{δ} \{δ − (1 + δ)nC_{5} \frac{(x^{T} s)^{1−3δ}}{η^{3}} − C_{5} \frac{(x^{T} s)^{1−4δ}}{η^{4}} \}
≥ \frac{1}{δ} \{δ − \frac{(x^{T} s)^{1−4δ}}{η^{4}} [C_{5} + 2nC_{5}η(x^{T} s)^{δ}]\}
≥ 1 − \frac{γ}{δη^{4}} (x^{T} s)^{1−4δ},

where the last inequality follows from the fact that η(x^{T} s)^{δ} ≤ 1, γ ≥ γ̄ and relation (69). We have thus shown that (70) holds for every α satisfying (67). The validity of relation (68) guarantees that (x + αΔx, y + αΔy, s(α)) ∈ N_{4}(δ, η, γ) for all α satisfying (67). Therefore, to complete the proof of the lemma, it is sufficient to show that (68) holds. Indeed, using (72), (73) and (69), we obtain

(x + αΔx)^{T} s(α) ≤ x^{T} s [1 − α + αnC_{5} \frac{(x^{T} s)^{1−3δ}}{η^{3}} ]
≤ x^{T} s [1 − α + α\frac{nC_{5}δ}{γ} ]
≤ x^{T} s (1 − α (1 − \frac{δ}{2})) ,

for every α ∈ [0, 1].

The following lemma will be useful later in the selection of the step sizes α_{k} for algorithm PDA.

Lemma 4.4 Suppose that Assumptions 1–5 hold and let δ ∈ (0, 1/4), η ∈ (0, 1] and γ ≥ γ̄ be given, where γ̄ is the constant specified in Lemma 4.3. Let (x, y, s) ∈ N_{4}(δ, η, γ) be given. Consider the solution (Δx, Δy, Δs) of (18) and define α̂ as

α̂ = sup{α ≥ 0 | x + αΔx ≥ 0, s + αΔs ≥ 0}.  (76)

Then,

|α̂ − 1| ≤ \frac{γ}{δη^{4}} (x^{T} s)^{1−4δ}.

Proof. ¿From the definition of α̂, it is easy to show that

α̂^{−1} = max(−X^{−1}Δx, −S^{−1}Δs).  (77)

It follows from relations (23), (24), (62) and (54) with p = 4 and q = 5/2, and Lemma 3.7 that

max \{max_{i∈B} (\frac{|Δx_{i}|}{x_{i}}) , max_{i∈N} (\frac{|Δs_{i}|}{s_{i}})\} ≤ \frac{C_{2}(x^{T} s)^{1−(3/2)δ}/η^{3/2}}{rη(x^{T} s)^{δ}}
≤ \frac{C_{6}(x^{T} s)^{1−(5/2)δ}}{η^{5/2}}  (78)
≤ \frac{C_{6}(x^{T} s)^{1−4δ}}{η^{4}} .

By (69) and (66), we have C_{6} ≤ C_{5} ≤ γ̄/(2n + 1) ≤ γ/2. Using this fact in (78), together with δ < 1 and (x, y, s) ∈ N_{4}(δ, η, γ), we find that

max \{max_{i∈B} (\frac{|Δx_{i}|}{x_{i}}) , max_{i∈N} (\frac{|Δs_{i}|}{s_{i}})\} ≤ \frac{C_{6}(x^{T} s)^{1−4δ}}{η^{4}}
≤ \frac{γ}{2δη^{4}} (x^{T} s)^{1−4δ} ≤ \frac{1}{2}.  (79)

¿From (18a), we have that

\frac{Δx_{i}}{x_{i}} + \frac{Δs_{i}}{s_{i}} = −1, ∀i = 1, …, n.

Using this relation and (79), we obtain

\begin{matrix}|\frac{Δx_{i}}{x_{i}} + 1| ≤ \frac{|Δs_{i}|}{s_{i}} ≤ \frac{γ}{2δη^{4}} (x^{T} s)^{1−4δ} ≤ \frac{1}{2},  ∀i ∈ N, \\ |\frac{Δs_{i}}{s_{i}} + 1| ≤ \frac{|Δx_{i}|}{x_{i}} ≤ \frac{γ}{2δη^{4}} (x^{T} s)^{1−4δ} ≤ \frac{1}{2},  ∀i ∈ B. \end{matrix}  (80)

¿From (79) and (80) it follows that

max(−X^{−1}Δx, −S^{−1}Δs) = max(−X_{N}^{−1}Δx_{N}, −S_{B}^{−1}Δs_{B}).  (81)

Using (80) again and this last relation, we obtain

|α̂^{−1} − 1| ≤ \frac{γ}{2δη^{4}} (x^{T} s)^{1−4δ} ≤ \frac{1}{2}.  (82)

Hence,

|α̂ − 1| ≤ α̂\frac{γ}{2δη^{4}} (x^{T} s)^{1−4δ} ≤ \frac{γ}{δη^{4}} (x^{T} s)^{1−4δ},

since α̂ ≤ 2 due to (82).

Finally we show that, given (strong) assumptions on the starting point for Algorithm PDA, a sequence of iterates for which the duality gap x^{k}^{T} s^{k} converges superlinearly to zero can be generated.

Theorem 4.5 Suppose that the Assumptions 1–5 hold and let δ ∈ (0, 1/4) be given. Let γ̄ be the constant as in the statement of Lemma 4.3. Assume that a constant γ ≥ γ̄ and a strictly feasible solution (x^{0}, y^{0}, s^{0}) satisfy

min_{i} x_{i}^{0}s_{i}^{0} ≥ (\frac{2γ}{δ} )^{1/4} (x^{0}^{T} s^{0})^{5/4} ,  (83)

and

η =^{△} (\frac{2γ}{δ} )^{1/4} (x^{0}^{T} s^{0})^{(1−4δ)/4} ≤ 1.  (84)

Then the following statements hold.

(i) If (x^{0}, y^{0}, s^{0}) is used as the starting point for Algorithm PDA and the sequence of step sizes {α_{k}} satisfies

α_{k} ∈ J_{k} =^{△} [0, 1 − \frac{γ}{δη^{4}} (x^{k}^{T} s^{k})^{1−4δ}] ,  ∀k = 1, 2, …  (85)

then

(x^{k}, y^{k}, s^{k}) ∈ N_{4}(δ, η, 2γ), ∀k ≥ 0.

(ii) Let δ̄ be an arbitrary constant satisfying δ̄ ∈ (0, 1 − 4δ). If the sequence of step sizes {α_{k}} is defined as

α_{k} = max (\frac{1}{3}, 1 − σ (x^{k}^{T} s^{k})δ̄) α̂_{k},  ∀k ≥ 0,  (86)

where

σ = (x^{0}^{T} s^{0})^{−(1−4δ)} = \frac{2γ}{δη^{4}} ,  (87)

and α̂_{k} is any number satisfying

|α̂_{k} − 1| ≤ \frac{γ}{δη^{4}} (x^{k}^{T} s^{k})^{1−4δ} ,  (88)

then, α_{k} ∈ J_{k} for all k ≥ 0 and the sequence of duality gaps corresponding to {x_{k}, y_{k}, s_{k}} eventually converges q-superlinearly to zero according to

x^{k+1}^{T} s^{k+1} ≤ \frac{31}{9}σ (x^{k}^{T} s^{k})^{(1+δ̄)} .  (89)

Proof. We first show (i). In view of Lemma 4.3, it is sufficient to show that {x_{0}, y_{0}, s_{0}} ∈ N_{4}(δ, η, 2γ). Indeed, relation (84) implies that x^{0}^{T} s^{0} ≤ 1 and

(x^{0}^{T} s^{0})^{1−4δ}/(δη^{4}) = 1/(2γ).

Moreover, using (83) and (84), we obtain

min_{i} x_{i}^{0}s_{i}^{0} ≥ (\frac{2γ}{δ} )^{1/4} (x^{0}^{T} s^{0})^{5/4}
= (\frac{2γ}{δ} )^{1/4} (x^{0}^{T} s^{0})^{(1−4δ)/4}(x^{0}^{T} s^{0})^{1+δ}
= η(x^{0}^{T} s^{0})^{1+δ}.

Hence, {x_{0}, y_{0}, s_{0}} ∈ N_{4}(δ, η, 2γ).

We now show (ii), starting with a proof that α_{k} ∈ J_{k} for all k ≥ 0. Indeed, using relation (88) and the fact that (x^{k}, y^{k}, s^{k}) ∈ N_{4}(δ, η, 2γ), and in particular that

(x^{k}^{T} s^{k})^{1−4δ}/(δη^{4}) ≤ 1/(2γ),  (90)

we obtain

\frac{1}{3}α̂_{k} ≤ \frac{1}{3} (1 + \frac{γ}{δη^{4}} (x^{k}^{T} s^{k})^{1−4δ})
≤ \frac{1}{3} (1 + \frac{1}{2}) = \frac{1}{2} ≤ 1 − \frac{γ}{δη^{4}} (x^{k}^{T} s^{k})^{1−4δ} .  (91)

Moreover, from (88), (87) and the fact that δ̄ ∈ (0, 1 − 4δ), it follows that

(1 − σ (x^{k}^{T} s^{k})δ̄) α̂_{k}
≤ (1 − \frac{2γ}{δη^{4}} (x^{k}^{T} s^{k})δ̄) (1 + \frac{γ}{δη^{4}} (x^{k}^{T} s^{k})^{1−4δ})
≤ 1 − \frac{2γ}{δη^{4}} (x^{k}^{T} s^{k})δ̄ + \frac{γ}{δη^{4}} (x^{k}^{T} s^{k})^{1−4δ}
≤ 1 − \frac{γ}{δη^{4}} (x^{k}^{T} s^{k})^{1−4δ} .  (92)

Together, (91) and (92) imply that α_{k} ∈ J_{k}. Before showing that (89) holds, we show that x^{k}^{T} s^{k} converges at least q-linearly to zero. Indeed, by (90) and (88), we have

α_{k} ≥ \frac{1}{3}α̂_{k} ≥ \frac{1}{3} [1 − \frac{γ}{δη^{4}} (x^{k}^{T} s^{k})^{1−4δ}] ≥ \frac{1}{6}.

Using this inequality together with (68) and δ ∈ (0, 1/4), we obtain

x^{k+1}^{T} s^{k+1} ≤ (1 − α_{k}(1 − δ/2)) x^{k}^{T} s^{k} ≤ (1 − \frac{1}{6}(1 − δ/2)) x^{k}^{T} s^{k} < \frac{41}{48}x^{k}^{T} s^{k}, (93)

as required. We now show that (89) holds. Because of the q-superlinear convergence property, we can choose an integer K such that

σ(x^{k}^{T} s^{k})δ̄ ≤ 2/3,  ∀k ≥ K.  (94)

In view of (86), this implies that

α_{k} = (1 − σ(x^{k}^{T} s^{k})δ̄) α̂_{k},  ∀k ≥ K.  (95)

¿From relations (88), (95) and (72), it follows that

x^{k+1}^{T} s^{k+1} ≤ x^{k}^{T} s^{k} [1 − α_{k} (1 − nC_{5} \frac{(x^{k}^{T} s^{k})^{1−3δ}}{η^{3}} )] ≤

x^{k}^{T} s^{k} [1 − (1 − σ(x^{k}^{T} s^{k})δ̄) (1 − \frac{γ}{δη^{4}} (x^{k}^{T} s^{k})^{1−4δ}) (1 − nC_{5} \frac{(x^{k}^{T} s^{k})^{1−3δ}}{η^{3}} )] ,

for all k ≥ K. Using (69), (87), relation (54) with p = 4 and q = 3, and (94), we have that

nC_{5} \frac{(x^{k}^{T} s^{k})^{1−3δ}}{η^{3}} ≤ \frac{γ}{δ} \frac{(x^{k}^{T} s^{k})^{1−4δ}}{η^{4}} ≤ \frac{σ}{2} (x^{k}^{T} s^{k})δ̄ ≤ \frac{1}{3}.

Hence,

x^{k+1}^{T} s^{k+1} ≤ (x^{k}^{T} s^{k}) [1 − (1 − σ(x^{k}^{T} s^{k})δ̄)^{3}]
≤ (x^{k}^{T} s^{k}) [3σ(x^{k}^{T} s^{k})δ̄ + σ^{3}(x^{k}^{T} s^{k})^{3δ̄}]
≤ (3 + \frac{4}{9}) σ(x^{k}^{T} s^{k})^{1+δ̄},

for every k ≥ K. Hence, (89) holds.

A number of choices for α̂_{k} that insures the validity of (88) are possible. The simplest is to set α̂_{k} ≡ 1 but, from Lemma 4.4, we see that choosing α̂_{k} according to (76) is another possibility. A third possibility is to choose

α̂_{k} = sup {α | (x^{k} + θΔx^{k}, s^{k}(θ)) ≥ 0 for all θ ∈ [0, α] } ,

which can be shown, by a technique similar to that used in Lemma 4.4, to satisfy (88). However, this last choice is not attractive since α̂_{k} is not easily computable.

5 A globally and superlinearly convergent potential reduction algorithm

In this section we embed the locally convergent algorithm of the previous section in the potential reduction algorithm of Monteiro [9].

We start by outlining a version of the algorithm in [9] in which the centering parameter (namely, the parameter σ_{k} in the notation of [9]) is chosen to be zero. The directions so obtained are referred to as primal-dual affine scaling directions. For simplicity of exposition, we fix certain parameters in his algorithm to numerical values. The algorithm makes use of the Tanabe-Todd-Ye potential function:

ϕ_{q}(x, s) = q log x^{T} s − ∑_{i=1}^{n} log x_{i}s_{i}  (96)

where q > n.

Algorithm PDPR

initially: choose fixed constants β > 1 and μ ∈ (0, 1) and let (x^{0}, y^{0}, s^{0}) be a strictly feasible point;

for k = 1, 2, …

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

Set τ_{k} = θ_{k}α̂_{k} where θ_{k} ≥ 1/3 and

α̂_{k} ≡ max{α | x_{k} + αΔx_{k} ≥ 0, s_{k} + αΔs_{k} ≥ 0};  (97)

Define α_{k} = β^{−m_{k}}τ_{k}, where m_{k} is the smallest non-negative integer for which

x^{k} + β^{−m_{k}}τ_{k}Δx^{k} > 0, s^{k}(β^{−m_{k}}τ_{k}) > 0,

and

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

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

end for

Some comments are in order. First, there is nothing special with the constant 1/3 which gives a lower bound for θ_{k}. In fact, any positive constant could be used in place of 1/3. Second, α_{k} defined in (97) is well-defined in the sense that the maximum of the set exists. This can be shown as follows. By relation (18a), it follows that x^{k}^{T} Δs^{k}+s^{k}^{T} Δx^{k} = −x^{k}^{T} s^{k} < 0. Hence, some component of either Δx^{k} or Δs^{k} must be negative. Thus, the maximum in (97) is finite.

The following global convergence theorem can be proved as a consequence of the results presented in Monteiro [9].

Theorem 5.1 In Algorithm PDPR,

lim_{k→∞} x^{k}^{T} s^{k} = 0,

and the sequence {(x^{k}, y^{k}, s^{k})} has at least one limit point. If (x̄, ȳ, s̄) is one such limit point then x̄ solves the primal problem (1), while (x̄, ȳ) solves the corresponding dual problem (2).

The theorem above guarantees that Algorithm PDPR converges globally from any feasible starting point. To prove local superlinear convergence of Algorithm PDPR, we need the assumption stated below that some iterate (x^{K}, y^{K}, s^{K}) of the algorithm satisfies a condition similar to (83) and (84). In the final result of this section, we show that this assumption can be dispensed with for certain values of q.

Assumption 6 For fixed constants γ ≥ γ̄ and δ ∈ (0, 1/4), assume that an integer K > 0 is known for which the iterate (x^{K}, y^{K}, s^{K}) of Algorithm PDPR satisfies

min_{i} x_{i}^{K}s_{i}^{K} ≥ (\frac{2γ}{δ} )^{1/4} (x^{K}^{T} s^{K})^{5/4} , and  (99)
η = [\frac{2γ}{δ} ]^{1/4} (x^{K}^{T} s^{K})^{(1−4δ)/4} ≤ 1.  (100)

If Assumption 6 holds then the following rule for selecting the sequence of initial step sizes can be employed:

τ_{k} =^{△} max (\frac{1}{3}, 1 − σ_{k} (x^{k}^{T} s^{k})δ̄) α̂_{k},  (101)

where δ̄ ∈ (0, 1 − 4δ) and

σ_{k} = \{ \begin{matrix}any value, & if k < K; \\ σ & if k ≥ K, \end{matrix}  (102)

with

σ =^{△} \frac{2γ}{δη^{4}} = (x^{K}^{T} s^{K})^{−(1−4δ)}.

Superlinear convergence follows from the results of Section 4 if we can show that α_{k} = τ_{k} for all k sufficiently large. That is, the relation (98) must be satisfied with m_{k} = 0 for all sufficiently large k. We have the following result:

Theorem 5.2 Suppose that Assumptions 1–6 hold. Assume that the sequence of initial trial step sizes is chosen as described in the previous paragraph. Then Algorithm PDPR converges superlinearly with q-order at least 1 + δ̄.

Proof. First, we observe that if (x^{k}, y^{k}, s^{k}) ∈ N (δ, η, 2γ) then Lemma 4.4 and Theorem 4.5(ii) imply that the initial k-th step size satisfies τ_{k} ∈ J_{k}. Since α_{k} ≤ τ_{k}, we also have that α_{k} ∈ J_{k}. Lemma 4.3 then implies that (x^{k+1}, y^{k+1}, s^{k+1}) ∈ N (δ, η, 2γ). Since (x^{K}, y^{K}, s^{K}) ∈ N (δ, η, 2γ), the above argument implies that (x^{k}, y^{k}, s^{k}) ∈ N (δ, η, 2γ), for all k ≥ K. As mentioned in the discussion above, superlinear convergence with q-order at least 1+δ̄ follows as a consequence of Theorem 4.5(ii) once we show that τ_{k} satisfies the following relation

ϕ_{q}(x^{k} + τ_{k}Δx^{k}, s^{k}(τ_{k})) − ϕ_{q}(x^{k}, s^{k}) ≤ μτ_{k}∇ϕ_{q}(x^{k}, s^{k})^{T} [ \begin{matrix}Δx^{k} \\ Δs^{k} \end{matrix} ] ,  (103)

for every k sufficiently large. Indeed, from Theorem 5.1, we know that x^{k}^{T} s^{k} → 0 as k tends to ∞, that is Algorithm PDPR converges globally. Therefore, we can choose K_{1} ≥ K such that

σ (x^{k}^{T} s^{k})δ̄ ≤ \frac{2}{3}, ∀k ≥ K_{1},

and so

τ_{k} = max (\frac{1}{3}, 1 − σ (x^{k}^{T} s^{k})δ̄) α̂_{k} = (1 − σ(x^{k}^{T} s^{k})δ̄) α̂_{k},  ∀k ≥ K_{1}.  (104)

Moreover, it follows from (104) and Lemma 4.4 that

|τ_{k} − (1 − σ(x^{k}^{T} s^{k})δ̄)| = (1 − σ(x^{k}^{T} s^{k})δ̄) |α̂_{k} − 1|
≤ |α̂_{k} − 1|
≤ \frac{γ}{δη^{4}} (x^{k}^{T} s^{k})^{1−4δ}
≤ \frac{σ}{2} (x^{k}^{T} s^{k})δ̄.

Hence, we obtain that

(1 − \frac{3}{2}σ(x^{k}^{T} s^{k})δ̄) ≤ τ_{k} ≤ (1 − \frac{1}{2}σ(x^{k}^{T} s^{k})δ̄) .  (105)

In the remainder of the proof, we drop the superscript and subscript k for the sake of clarity. The right hand side of (103) is easy to evaluate. Indeed, from (18a) and the fact that by (105), τ ≤ 1, we have

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

= μτ \{[\frac{q}{x^{T} s}s − X^{−1}e]^{T} Δx + [\frac{q}{x^{T} s}x − S^{−1}e]^{T} Δs\}

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

We now evaluate the left hand side of (103). ¿From (57) it follows that

\frac{(x + αΔx)^{T} s(α)}{x^{T} s} ≤ 1 − α + nC_{5} \frac{(x^{T} s)^{1−3δ}}{η^{3}} ,  ∀α ∈ [0, 1].  (106)

Moreover, it follows from (56) that for i = 1, …, n

\begin{aligned}\frac{(x_{i}+αΔx_{i})s_{i}(α)}{x_{i}s_{i}} & ≥ & 1 − α − C_{5} \frac{(x^{T} s)^{2−3δ}}{η^{3}(x_{i}s_{i})} \\ & ≥ & 1 − α − C_{5} \frac{(x^{T} s)^{1−4δ}}{η^{4}} ,  ∀α ∈ [0, 1], \end{aligned}  (107)

where the last inequality is due to the fact that x_{i}s_{i} ≥ η(x^{T} s)^{1+δ}. Using (106), (107) and the definition of ϕ_{q} given in (96), we obtain

ϕ_{q}(x + αΔx, s(α)) − ϕ_{q}(x, s)
= q log [\frac{(x + αΔx)^{T} s(α)}{x^{T} s}] − ∑_{i=1}^{n} log [\frac{(x_{i} + αΔx_{i})s_{i}(α)}{x_{i}s_{i}}]
≤ q log [1 − α + nC_{5} \frac{(x^{T} s)^{1−3δ}}{η^{3}} ] − ∑_{i=1}^{n} log [1 − α − C_{5} \frac{(x^{T} s)^{1−4δ}}{η^{4}} ]
= (q − n) log(1 − α) + q log [1 + \frac{nC_{5}(x^{T} s)^{1−3δ}}{η^{3}(1 − α)} ]

− n log [1 − \frac{C_{5}(x^{T} s)^{1−4δ}}{η^{4}(1 − α)} ] ,

for all α ∈ [0, 1]. Setting α = τ_{k} in the above relation and using (105), we obtain

ϕ_{q}(x + τΔx, s(τ)) − ϕ_{q}(x, s)
≤ (q − n) log [\frac{3}{2}σ(x^{T} s)δ̄] + q log [1 + \frac{nC_{5}(x^{T} s)^{1−3δ−δ̄}}{η^{3}(σ/2)}]

− n log [1 − \frac{C_{5}(x^{T} s)^{1−4δ−δ̄}}{η^{4}(σ/2)}] .

Consider each of the three terms in this expression. Since x^{k}^{T} s^{k} → 0 and δ̄ ∈ (0, 1 − 4δ), we have for k sufficiently large that

(q − n) log [\frac{3σ}{2} (x^{T} s)δ̄] ≤ −2(q − n)μ,

q log [1 + \frac{nC_{5}(x^{T} s)^{1−3δ−δ̄}}{η^{3}(σ/2)}] ≤ \frac{q − n}{2}μ,
−n log [1 − \frac{C_{5}(x^{T} s)^{1−4δ−δ̄}}{η^{4}(σ/2)}] ≤ \frac{q − n}{2}μ.

Hence, for k sufficiently large,

ϕ_{q}(x^{k} + τ_{k}Δx^{k}, s^{k}(τ_{k})) − ϕ_{q}(x^{k}, s^{k}) ≤ −(q − n)μ ≤ μτ_{k}∇ϕ_{q}(x^{k}, s^{k})^{T} [ \begin{matrix}Δx^{k} \\ Δs^{k} \end{matrix} ] .

We now address two important issues in the implementation of Algorithm PDPR, namely, the choice of γ greater than the (unknown) constant γ̄, and the availability of an index K such that the conditions (99) and (100) are satisfied. We address the latter issue first. Given any γ > 0, we show that a K satisfying the two conditions of Assumption 6 is always available provided that the parameter q that defines the potential function ϕ_{q} lies in the range (n, n + 1/4).

Theorem 5.3 Let γ > 0 and δ ∈ (0, 1/4) be given and assume that q ∈ (n, n + 1/4). Then there exists an integer K such that both (99) and (100) hold.

Proof. For convenience, we define the quantity ν = q − n so that ν ∈ (0, 1/4). Since the potential function ϕ_{q} decreases at each iteration, we have

ϕ_{q}(x^{k}, s^{k}) ≤ ϕ_{q}(x^{0}, s^{0}) =^{△} M.

Now, for any i ∈ {1, 2, …, n}, we have

(1 + ν) log x^{k}^{T} s^{k} − log x_{i}^{k}s_{i}^{k}
≤ M − (n − 1) log x^{k}^{T} s^{k} + ∑_{j≠i} log x_{j}^{k}s_{j}^{k}
≤ M − (n − 1) log (x^{k}^{T} s^{k} − x_{i}^{k}s_{i}^{k}) + ∑_{j≠i} log x_{j}^{k}s_{j}^{k}
≤ M − (n − 1) log(n − 1)
=^{△} M̄.

Hence

log \frac{x_{i}^{k}s_{i}^{k}}{(x^{k}^{T} s^{k})^{1+ν}} ≥ −M̄

⟹  x_{i}^{k}s_{i}^{k} ≥ e^{−M̄} (x^{k}^{T} s^{k})^{1+ν}
⟹  x_{i}^{k}s_{i}^{k} ≥ [\frac{2γ}{δ} ]^{1/4} (x^{k}^{T} s^{k})^{5/4} \{e^{−M̄} [\frac{δ}{2γ} ]^{1/4} (x^{k}^{T} s^{k})^{ν−1/4}\} .

Hence (99) will be satisfied for k = K if we have

e^{−M̄} [\frac{δ}{2γ} ]^{1/4} (x^{K}^{T} s^{K})^{ν−1/4} ≥ 1.

Since x^{k}^{T} s^{k} → 0 by Theorem 5.1, this last condition, and hence (99), is guaranteed to hold for any sufficiently large choice of K. By increasing K if necessary, we can also ensure that (100) is also satisfied. Hence, the result follows.

In the remainder of the paper, we assume that q ∈ (n, n + 1/4).

The remaining issue is the choice of a constant γ. In order for the theory of sections 4 and 5 to hold, we require that γ ≥ γ̄, where the constant γ̄, which depends on the problem data but not on δ, is not known. Our strategy is similar to one that has been suggested for the choice of penalty parameter in the nonsmooth penalty function approach for constrained optimization (see Fletcher [2]); namely, to embed Algorithm PDPR in an outer loop in which γ is successively increased until a suitable value is found. We formalize this procedure as follows:

Algorithm G-PDPR

initially: Choose fixed constants β > 1, δ̄ ∈ (0, 1), δ ∈ (0, (1 − δ̄)/4), μ ∈ (0, 1) and γ_{0} ≥ δ;

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

for ℓ = 0, 1, …

Phase 1:

Apply Algorithm PDPR with γ = γ_{ℓ}, defining the initial steps τ_{k} by (101), with σ_{k} arbitrary;

terminate Phase 1 at the first iterate k = K_{ℓ} for which (99) and (100) hold with K = K_{ℓ}.

Phase 2:

Set σ = (x^{K_{ℓ}}^{T} s^{K_{ℓ}})^{−(1−4δ)};

Set η = (2γ/δ)^{1/4}(x^{K_{ℓ}}^{T} s^{K_{ℓ}})^{(1−4δ)/4};

Continue Algorithm PDPR with σ_{k} ≡ σ;

terminate Phase 2 if any of the following three conditions is violated by the current iterate (x^{k}, y^{k}, s^{k}):

(x^{k} + τ_{k}Δx^{k}, y^{k} + τ_{k}Δy^{k}, s^{k}(τ_{k})) ∈ N (δ, η, 2γ);  (108)
|α̂^{k} − 1| ≤ \frac{γ}{δη^{4}} (x^{k}^{T} s^{k})^{1−4δ};  (109)
(x^{k} + τ_{k}Δx^{k})^{T} s^{k}(τ_{k}) ≤ \frac{31}{9}σ (x^{k}^{T} s^{k})^{(1+δ̄)} .  (110)
Set γ_{ℓ+1} ← 10γ_{ℓ} and ℓ ← ℓ + 1;

end for

The increase factor of 10 for γ can of course be replaced by any constant that exceeds 1. Note that Algorithm G-PDPR is a special case of Algorithm PDPR, in which we define the choice of initial step size in a somewhat elaborate way. Note also that the choices of δ̄ and δ guarantee that δ ∈ (0, 1/4) and that δ̄ ∈ (0, 1 − 4δ).

The following result guarantees global and superlinear convergence of Algorithm G-PDPR.

Theorem 5.4 Suppose that q ∈ (n, n+1/4) and that Assumptions 1–5 hold. Then Algorithm G-PDPR converges converges globally and superlinearly with q-order at least 1 + δ̄

Proof. As noted above, Algorithm G-PDPR is simply a special case of Algorithm PDPR, and so global convergence follows as a consequence of Theorem 5.1.

To prove local convergence properties, we show first that the number of “outer loops” in Algorithm G-PDPR is finite. Let us assume for contradiction that this is not true. Then there must be a loop index ℓ for which γ_{ℓ} > γ̄. Theorem 5.3 implies that Phase 1 of loop ℓ terminates, while Lemma 4.4 and Theorem 4.5 ensure that conditions (108), (109), and (110) are not violated on any subsequent iteration. This means that Phase 2 of loop ℓ does not terminate, giving a contradiction.

Now consider the final loop of Algorithm G-PDPR. Since all the conditions (108), (109) and (110) are satisfied for every k sufficiently large, it is clear from (110) that superlinear convergence with Q-order at least (1 + δ̄) is obtained if we can show that α_{k} = τ_{k}, for all k sufficiently large. The proof that this holds is identical to that of Theorem 5.2. We need to observe that (105) can be obtained by using relation (109), while the relations (106) and (107) follow from Lemma 4.2 and (108).

6 Conclusions and discussion

We have described an algorithm for general convex programming problems that converges globally and superlinearly under fairly weak assumptions. An interesting feature of the global convergence framework is that the choice of q for the potential function does not depend on δ. Since smaller values of q−n are typically (and, for some algorithms, rigorously) associated with slower convergence, the choices of q we have described in Section 5 are preferable to the more obvious choice q = n + δ, which would have made for simpler analysis.

The analysis can be generalized in at least two ways. First, a centering parameter can be introduced into the step calculation (18). Equation (18a) can be replaced by

S^{k}Δx^{k} + X^{k}Δs^{k} = −S^{k}X^{k}e + σ_{k}(x^{k}^{T} s^{k})/n,

where σ_{k} ∈ [0, σ], σ ∈ [0, 1). The second term on the right hand side serves to bias the step towards the central path, and is a common feature of most interior point algorithms. The global convergence analysis in Monteiro [9] still holds when σ_{k} is chosen in this way, while our local convergence analysis still goes through provided that σ_{k} approaches zero sufficiently rapidly (that is, like some power of (x^{k}^{T} s^{k})) during the latter stages of the algorithm.

A second possible generalization is to functions whose Hessian satisfies a Hölder continuity condition

∥∇^{2}f(x^{1}) − ∇^{2}f(x^{2})∥ ≤ C_{x}∥x^{1} − x^{2}∥^{γ}

(where γ ∈ (0, 1]) in place of the Lipschitz continuity condition in Assumption 2. In this case, it is still possible to design an algorithm with global and superlinear convergence. The q-order for the convergence of the duality gap to zero may lie in the range [1, 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 thanks Argonne for its financial support and for the congenial scientific atmosphere that it provided.

References

[1] J. V. Burke and M. C. Ferris, Weak sharp minima in mathematical programming, Research Report 1050, iComputer Sciences Department, University of Wisconsin, Wisconsin, Madison, 1991. SIAM Journal on Control and Optimization, to appear.

[2] R. Fletcher, Practical Methods of Optimization, John Wiley & Sons, New York, second ed., 1987.

[3] J. Ji, F. Potra, and S. Huang, A predictor-corrector method for linear complementarity problems with polynomial complexity and superlinear convergence, Technical Report, Department of Mathematics, The University of Iowa, Iowa City, Iowa, 52240, USA, August 1991.

[4] J. Ji, F. Potra, R. A. Tapia, and Y. Zhang, An interior–point method with polynomial complexity and superlinear convergence for linear complementarity problems, Technical Report TR–91–23, Dept. of Mathematical Sciences, Rice University, Houston, TX 77251, USA, 1991.

[5] M. Kojima, Y. Kurita, and S. Mizuno, Large–step interior point algorithms for linear complementarity problems, Research Reports on Information Sciences, Ser. B : Operations Research B–243, Dept. of Information Sciences, Tokyo Institute of Technology, Oh–Okayama, Meguro–ku, Tokyo 152, Japan, 1991.

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

[7] K. A. McShane, A superlinearly convergent O(\sqrt{n}L) iteration primal–dual linear programming algorithm, manuscript, U. S. Government, 2537 Villanova Drive, Vienna, Virginia, USA, 1991.

[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] 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, Tucson, AZ 85721, USA, July 1991.

[10] 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.

[11] Y. Ye, On the q-order of convergence of interior-point algorithms for linear programming, Working Paper 91–17, The College of Business Administration, The University of Iowa, Iowa City, IA 52242, USA, 1991.

[12] Y. Ye and K. Anstreicher, On quadratic and O(\sqrt{n}L) convergence of a predictorcorrector algorithm for LCP, Technical Report, Department of Management Sciences, The University of Iowa, Iowa City, IA 52242, USA, November 1991.

[13] Y. Ye, O. Güler, R. A. Tapia, and Y. Zhang, A quadratically convergent O(\sqrt{n}L)–iteration algorithm for linear programming, Technical Report TR–91–26, Dept. of Mathematical Sciences, Rice University, Houston, TX 77251, USA, August 1991.

[14] Y. Zhang and R. A. Tapia, A quadratically convergent polynomial primal–dual interior–point algorithm for linear programming, Technical Report TR 90–40, Dept. of Mathematical Sciences, Rice University, Houston, TX 77251, USA, 1990.

[15] , 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.

[16] Y. Zhang, R. A. Tapia, and J. E. Dennis, On the superlinear and quadratic convergence of primal–dual interior point linear programming algorithms, Technical Report TR 90–6, Dept. of Mathematical Sciences, Rice University, Houston, TX 77251, USA, January 1990.

[17] Y. Zhang, R. A. Tapia, and F. Potra, On the superlinear convergence of interior point algorithms for a general class of problems, Technical Report TR 90–9, Dept. of Mathematical Sciences, Rice University, Houston, TX 77251, USA, March 1990.