An Infeasible-Interior-Point Algorithm for Linear Complementarity Problems^{∗}

Stephen Wright^{†}

October 15, 1992

PREPRINT MCS-P331-1092, MCS DIVISION, ARGONNE NATIONAL LABORATORY

Abstract

In this paper, we discuss a polynomial and Q-subquadratically convergent algorithm for linear complementarity problems that does not require feasibility of the initial point or the subsequent iterates. The algorithm is a modification of the linearly convergent method of Zhang and requires the solution of at most two linear systems with the same coefficient matrix at each iteration.

1 Introduction

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

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

where h ∈ IR^{n} and M is an n × n positive semidefinite matrix. A pair (x, y) is said to be feasible if y = Mx+h and (x, y) ≥ (0, 0), and strictly feasible if the latter inequality is strict.

It is well known that convex quadratic programming problems and linear programming problems can be expressed as linear complementarity problems; the same is true of extended linear-quadratic programming problems (see Rockafellar [9]). Much research has been devoted to interior point methods for (1). Recently, Ji, Potra, and Huang [1] proposed a predictor-corrector algorithm with polynomial complexity and superlinear convergence, while the predictor-corrector algorithm of Ye and Anstreicher [11] is polynomial and Qquadratic. In the latter paper, it was assumed only that a strictly feasible point and a strictly complementary solution (one for which max(x_{i}^{∗}, y_{i}^{∗}) > 0 for i = 1, ⋯, n) exist for (1). Both these algorithms generate a sequence of strictly feasible iterates (x^{k}, y^{k}). A strictly feasible starting point (x^{0}, y^{0}) must therefore be supplied. To find such a point, one often must augment the problem in an artificial way.

∗This research was supported by the Applied Mathematical Sciences subprogram of the Office of Energy Research, U. S. Department of Energy, under Contract W-31-109-Eng-38.

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

More recently, research has focused on algorithms that generate sequences (x^{k}, y^{k}) for which (x^{k}, y^{k}) > 0 but possibly y^{k} ≠ Mx^{k} + h. These infeasible-interior-point methods more nearly reflect computational practice (see, for example, Lustig, Marsten, and Shanno [4]). Also, minor modifications to a solution to a “nearby” problem can produce an excellent starting point for the present problem — an advantage when the underlying problem to be solved is nonlinear. Infeasible algorithms for linear programming have been proposed by Kojima, Meggido, and Mizuno [2], Kojima, Mizuno, and Todd [3], and Potra [7, 8]. Of these, only Potra’s algorithms have both polynomial complexity and superlinear convergence properties. Potra’s methods are of the predictor-corrector type and require three systems of linear equations (two of which have the same coefficient matrix) to be solved at each iteration. The only infeasible-interior-point algorithm for more general problems than linear programs that we are aware of is due to Zhang [12]. He analyzes an algorithm for a class of problems that includes (1) and proves Q-linear convergence of the complementarity gap μ_{k} = (x^{k})^{T} y^{k}/n to zero. Polynomial complexity is obtained for a particular choice of starting point. The algorithm requires the solution of a single system of linear equations at each iteration.

In this paper, we propose modifications of Zhang’s algorithm that retain polynomial complexity and have the added feature that the sequence {μ_{k}} converges superlinearly to zero with Q-order 2. Only the mild assumptions of Ye and Anstreicher [11] are required. Our method requires the solution of at most two linear systems of equations with the same coefficient matrix at each iteration.

When this report was about to be issued, we received a new report by Zhang and Zhang [13] that describes an infeasible-interior-point algorithm that is similar to ours in some respects. They allow relaxed versions of the centering condition and the feasibility dominance condition (cf. below (4e) and (4d), respectively) to be used on some iterations, and they obtain similar convergence properties. However, their algorithm is applicable only to linear programming problems.

Our algorithm is specified in Section 2. Some technical results are proved in Section 3, while in Section 4, we prove Q-linear convergence and polynomial complexity. Results concerning boundedness of the steps are proved in Section 5. Finally, superlinear convergence properties are discussed in Section 6.

Unless otherwise specified, ∥·∥ denotes the Euclidean norm of a vector. Iteration numbers appear as superscripts on vectors and matrices and as subscripts on scalars.

2 The Algorithm

Given a starting point with (x^{0}, y^{0}) > (0, 0), the algorithm generates a sequence of iterates (x^{k}, y^{k}) > (0, 0). The desirability of each point is measured by the merit function

ϕ(x, y) = x^{T} y + ∥y − Mx − h∥

whose two terms measure the complementarity gap and infeasibility. Clearly, a vector pair (x^{∗}, y^{∗}) is a solution of (1) if and only if (x^{∗}, y^{∗}) ≥ (0, 0) and ϕ(x^{∗}, y^{∗}) = 0. We use the shorthand notation ϕ_{k} to denote ϕ(x^{k}, y^{k}).

In order to describe the step between successive iterates, we define

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

The step is calculated as follows.

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

[ \begin{matrix}M & −I \\ Y^{k} & X^{k} \end{matrix} ] [ \begin{matrix}Δx^{k} \\ Δy^{k} \end{matrix} ] = [ \begin{matrix}−h − Mx^{k} + y^{k} \\ −X^{k}Y^{k}e + σ̃μ_{k}e \end{matrix} ] .  (2)

Choose

α̃ = arg min_{α} ϕ(x^{k} + αΔx^{k}, y^{k} + αΔy^{k})  (3)

subject to

α ∈ [0, 1],  (4a)
x^{k} + αΔx^{k} > 0,  (4b)
y^{k} + αΔy^{k} > 0,  (4c)
(x^{k} + αΔx^{k})^{T} (y^{k} + αΔy^{k}) ≥ (1 − β̃)(1 − α)(x^{k})^{T} y^{k},  (4d)
(x_{i}^{k} + αΔx_{i}^{k})(y_{i}^{k} + αΔy_{i}^{k}) ≥ (γ̃/n)(x^{k} + αΔx^{k})^{T} (y^{k} + αΔy^{k}),  i = 1, ⋯, n. (4e)

It has been noted previously that (2) are simply the equations obtained by applying one iteration of Newton’s method to the nonlinear equations

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

starting from the point (x^{k}, y^{k}). The condition (4e), usually referred to as a centering condition, ensures that the iterates do not prematurely approach the edge of the non-negative orthant. The condition (4d) is a relaxation of the condition enforced by Zhang [12, formula (5.7)] to ensure that feasibility is given a higher priority than complementarity (Zhang uses β̃ = 0). Potra’s algorithms replace (4d) with an equality condition in which β̃ = 0. We allow β̃ > 0 to permit superlinear convergence, as will become clear in the analysis that follows.

We can now state our algorithm.

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

t_{0} ← 1, γ_{0} ← 2γ̄;

for k = 0, 1, 2, ⋯

if ϕ_{k} = ϕ(x^{k}, y^{k}) ≤ ϕ̄
then Solve (2)–(4) with σ̃ = μ_{k}, β̃ = γ̄^{t_{k}}, γ̃ = γ̄(1 + γ̄^{t_{k}});
if ϕ(x^{k} + α̃Δx^{k}, y^{k} + α̃Δy^{k}) ≤ ρϕ_{k}
then α_{k} ← α̃, β_{k} ← β̃, σ_{k} ← σ̃, γ_{k+1} ← γ̃;

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

(x^{k+1}, y^{k+1}) ← (x^{k}, y^{k}) + α_{k}(Δx^{k}, Δy^{k}); go to next k;

end if end if

Solve (2)–(4) with σ̃ ∈ [σ, 1/2], β̃ = 0, γ̃ = γ_{k};

α_{k} ← α̃, β_{k} ← 0, σ_{k} ← σ̃, γ_{k+1} ← γ̃;

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

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

go to next k;

end for.

The idea behind this algorithm is simple. While the merit function exceeds ϕ̄, it is identical to Zhang’s algorithm. One linear system of equations is solved at each iteration to yield what we refer to as a safe step. When the merit function falls below the threshold ϕ̄, the algorithm computes a step that may yield rapid convergence by setting σ_{k} equal to the current duality gap and allowing β̃ > 0 in (4d). We refer to such a step as a fast step. The coefficient matrices for the fast and safe steps are the same, and the integer variable t_{k} keeps track of the number of fast steps that have been taken prior to iteration k. If the fast step does not give a significant reduction in ϕ (by at least a factor of ρ), it is discarded, and we take a safe step instead.

This strategy may lead to our rejecting a fast step and taking a safe step that gives a smaller decrease in ϕ at some iterations. Note, however, that we pay a price for taking fast steps in that γ is decreased, and so the Q-linear convergence rate will not be as favorable on subsequent iterations. Therefore, it makes sense to reject fast steps unless they make substantial progress.

In the safe-step calculation, there is considerable scope for user-defined heuristics in the choice of σ_{k}. In a practical implementation, the value of σ_{k} may be adjusted according to the merit function decrease ϕ_{k−1} − ϕ_{k} on the previous step, and the value of α_{k−1}.

In the remainder of the paper, we analyze the convergence properties of this algorithm under the following assumptions:

Assumption 1 M is positive semidefinite.

This assumption implies that if (x^{1}, y^{1}) and (x^{2}, y^{2}) are any two points satisfying y = Mx = h, then

(y^{2} − y^{1})^{T} (x^{2} − x^{1}) = (x^{2} − x^{1})^{T} M(x^{2} − x^{1}) ≥ 0.  (5)

Assumption 2 The problem (1) has a strictly feasible point (x̄, ȳ). That is, (x̄, ȳ) > (0, 0) and ȳ = Mx̄ + h.

Assumption 3 The solution set for (1) is nonempty and, moreover, there is a strictly complementary solution (x^{∗}, y^{∗}).

With respect to (x^{∗}, y^{∗}), we define the partitioning {1, 2, ⋯, n} = N ∪ B, where

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

Assumptions 1, 2, and 3 will be assumed everywhere without being explicitly stated.

3 Technical Results

In this section, we prove a number of results that are needed in the analysis of Sections 4, 5, and 6. In the statement of many of our results, we refer to the pair (x^{k}, y^{k}), which is always understood to be an iterate generated by the algorithm of Section 2.

We start by showing that the dominance of feasibility over complementarity is not completely abandoned by our relaxed condition (4d), but still holds to within a certain constant.

Lemma 3.1

β̂ =^{△} ∏_{k=0}^{∞}(1 − β_{k}) > 0.

Proof. Inspection of the algorithm shows that on each iteration, we have either β_{k} = 0 (safe steps) or β_{k} = γ̄^{t}, t = 1, 2, ⋯ (fast steps). Therefore

β̂ ≥ ∏_{t=1}^{∞}(1 − γ̄^{t}),

and the right-hand side is bounded away from zero since γ̄ ∈ (0, 1).

If we define

ν_{0} = 1,  ν_{k} = ∏_{i=0}^{k−1}(1 − α_{i}),  k = 1, 2, ⋯,

the following result defines upper and lower bounds on μ_{k}.

Lemma 3.2 For all (x^{k}, y^{k}) generated by the algorithm and μ_{k} = (x^{k})^{T} y^{k}/n,

μ_{k} ≥ β̂ν_{k}μ_{0},
μ_{k} ≤ μ_{0} + \frac{1}{n}∥y^{0} − Mx^{0} − h∥.

Proof. By (4d),

μ_{k} ≥ (1 − β_{k−1})(1 − α_{k−1})μ_{k−1} ≥ ∏_{i=0}^{k−1}(1 − β_{i})ν_{k}μ_{0} ≥ β̂ν_{k}μ_{0},

giving the first inequality. The second inequality follows from nμ_{k} ≤ ϕ_{k} ≤ ϕ_{0} = nμ_{0} + ∥y^{0} − Mx^{0} − h∥.

Zhang [12] defines an “auxiliary sequence” (u^{k}, v^{k}) by defining an initial point (u^{0}, v^{0}) such that

(u^{0}, v^{0}) ≤ (x^{0}, y^{0}), v^{0} = Mu^{0} + h,

and subsequent iterates by

u^{k+1} = u^{k} + α_{k}(Δx^{k} + x^{k} − u^{k}), v^{k+1} = v^{k} + α_{k}(Δy^{k} + y^{k} − v^{k}).

He proves the following result. Lemma 3.3 (Zhang [12, Lemma 4.1]) For k ≥ 0,

(i) v^{k} = Mu^{k} + h; (ii) x^{k+1} − u^{k+1} = ν_{k}(x^{0} − u^{0}) ≥ 0 and y^{k+1} − v^{k+1} = ν_{k}(y^{0} − v^{0}) ≥ 0; (iii) If α_{K} = 1 for some K ≥ 0, then (x^{k}, y^{k}) = (u^{k}, v^{k}) and therefore (x^{k}, y^{k}) is strictly

feasible for all k > K. Using Lemma 3.3, we can bound some components of the iterates (x^{k}, y^{k}) Lemma 3.4 There is a constant C_{1} > 0 such that for all iterates (x^{k}, y^{k}),

i ∈ N ⇒ x_{i}^{k} ≤ C_{1}μ_{k},  (6a)
i ∈ B ⇒ y_{i}^{k} ≤ C_{1}μ_{k},  (6b)

and

i ∈ B ⇒ x_{i}^{k} ≥ γ̄/C_{1},  (7a)
i ∈ N ⇒ y_{i}^{k} ≥ γ̄/C_{1}.  (7b)

Proof. From the definition of (u^{k}, v^{k}), we have (x^{k} − x^{∗})^{T} (y^{k} − y^{∗})

= (x^{k} − u^{k} + u^{k} − x^{∗})^{T} (y^{k} − v^{k} + v^{k} − y^{∗})  (8)
= (x^{k} − u^{k})^{T} (y^{k} − v^{k}) + (u^{k} − x^{∗})^{T} (y^{k} − v^{k}) + (x^{k} − u^{k})^{T} (v^{k} − y^{∗}) + (u^{k} − x^{∗})^{T} (v^{k} − y^{∗}).

Now v^{k} = Mu^{k} + h and y^{∗} = Mx^{∗} + h, so by (5),

(u^{k} − x^{∗})^{T} (v^{k} − y^{∗}) ≥ 0.

Therefore, since (x^{∗})^{T} y^{∗} = 0, we have

−(x^{k})^{T} y^{∗} − (x^{∗})^{T} y^{k} + (x^{k})^{T} y^{k}

≥ (x^{k} − u^{k})^{T} (y^{k} − v^{k}) + (u^{k} − x^{∗})^{T} (y^{k} − v^{k}) + (x^{k} − u^{k})^{T} (v^{k} − y^{∗}),

and so

(x^{∗})^{T} y^{k} + (x^{k})^{T} y^{∗}
≤ (x^{k})^{T} y^{k} − (x^{k} − u^{k})^{T} (y^{k} − v^{k}) − (u^{k} − x^{∗})^{T} (y^{k} − v^{k}) − (x^{k} − u^{k})^{T} (v^{k} − y^{∗})
= (x^{k})^{T} y^{k} + (x^{k} − u^{k})^{T} (y^{k} − v^{k}) − (x^{k} − x^{∗})^{T} (y^{k} − v^{k}) − (x^{k} − u^{k})^{T} (y^{k} − y^{∗}). (9)

Now, by Lemma 3.3,

x^{k} − u^{k} ≥ 0, y^{k} − v^{k} ≥ 0, y^{k} ≥ 0, x^{k} ≥ 0,

and so

(y^{k})^{T} (x^{k} − u^{k}) ≥ 0, (x^{k})^{T} (y^{k} − v^{k}) ≥ 0.

By substitution in (9), we have

(x^{∗})^{T} y^{k} + (x^{k})^{T} y^{∗}

≤ (x^{k})^{T} y^{k} + (x^{k} − u^{k})^{T} (y^{k} − v^{k}) + (x^{∗})^{T} (y^{k} − v^{k}) + (y^{∗})^{T} (x^{k} − u^{k})
= (x^{k})^{T} y^{k} [1 + \frac{(x^{k} − u^{k})^{T} (y^{k} − v^{k})}{(x^{k})^{T} y^{k}} + \frac{(x^{∗})^{T} (y^{k} − v^{k})}{(x^{k})^{T} y^{k}} + \frac{(y^{∗})^{T} (x^{k} − u^{k})}{(x^{k})^{T} y^{k}}] .

Now, by Lemmas 3.2 and 3.3, we can bound the term in the square brackets by C̄_{1}, where

C̄_{1} =^{△} [1 + \frac{(x^{0} − u^{0})^{T} (y^{0} − v^{0})}{β̂(x^{0})^{T} y^{0}} + \frac{(x^{∗})^{T} (y^{0} − v^{0})}{β̂(x^{0})^{T} y^{0}} + \frac{(y^{∗})^{T} (x^{0} − u^{0})}{β̂(x^{0})^{T} y^{0}}] ,

and so

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

Hence, for i ∈ N,

x_{i}^{k}y_{i}^{∗} ≤ nC̄_{1}μ_{k} ⇒ x_{i}^{k} ≤ \frac{nC̄_{1}}{y_{i}^{∗}} μ_{k},

while for i ∈ B,

y_{i}^{k} ≤ \frac{nC̄_{1}}{x_{i}^{∗}} μ_{k}.

Hence (6) is obtained by taking

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

For (7a), we simply note that

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

The proof of (7b) is similar.

Assumption 2 can be used to show that the iterates remain bounded.

Lemma 3.5 There is a constant C_{2} > 0 such that for k ≥ 0 and i = 1, 2, ⋯, n,

0 < x_{i}^{k} ≤ C_{2}, 0 < y_{i}^{k} ≤ C_{2}.

Proof. Because of Assumption 2, there is a strictly feasible point (x̄, ȳ). Now

(x^{k} − x̄)^{T} (y^{k} − ȳ)
= (x^{k} − u^{k} + u^{k} − x̄)^{T} (y^{k} − v^{k} + v^{k} − ȳ)
= (x^{k} − u^{k})^{T} (y^{k} − v^{k}) + (u^{k} − x̄)^{T} (y^{k} − v^{k}) + (x^{k} − u^{k})^{T} (v^{k} − ȳ)
+(u^{k} − x̄)^{T} (v^{k} − ȳ)
= (x^{k} − u^{k})^{T} (y^{k} − v^{k}) + (u^{k} − x^{k})^{T} (y^{k} − v^{k}) + (x^{k} − x̄)^{T} (y^{k} − v^{k})

+(x^{k} − u^{k})^{T} (v^{k} − y^{k}) + (x^{k} − u^{k})^{T} (y^{k} − ȳ) + (u^{k} − x̄)^{T} (v^{k} − ȳ)

= −(x^{k} − u^{k})^{T} (y^{k} − v^{k}) + (x^{k})^{T} (y^{k} − v^{k}) − x̄^{T} (y^{k} − v^{k})
+(x^{k} − u^{k})^{T} y^{k} − (x^{k} − u^{k})^{T} ȳ + (u^{k} − x̄)^{T} (v^{k} − ȳ).  (10)

Now, since (u^{k}, v^{k}) and (x̄, ȳ) are both feasible, the last term in (10) is non-negative. Moreover,

x^{k} > 0, y^{k} − v^{k} ≥ 0 ⇒ (x^{k})^{T} (y^{k} − v^{k}) ≥ 0,
y^{k} > 0, x^{k} − u^{k} ≥ 0 ⇒ (y^{k})^{T} (x^{k} − u^{k}) ≥ 0,
(x^{k} − u^{k}) = ν_{k}(x^{0} − u^{0}),  (y^{k} − v^{k}) = ν_{k}(y^{0} − v^{0}).

Hence, from (10),

(x^{k} − x̄)^{T} (y^{k} − ȳ)
≥ −(x^{k} − u^{k})^{T} (y^{k} − v^{k}) − x̄^{T} (y^{k} − v^{k}) − ȳ^{T} (x^{k} − u^{k})
= −ν_{k}^{2}(x^{0} − u^{0})^{T} (y^{0} − v^{0}) − ν_{k}x̄^{T} (y^{0} − v^{0}) − ν_{k}ȳ^{T} (x^{0} − u^{0}),

and so, using Lemma 3.2 to bound (x^{k})^{T} y^{k}, we obtain

x̄^{T} y^{k} + ȳ^{T} x^{k}
≤ (x^{k})^{T} y^{k} + x̄^{T} ȳ + (x^{0} − u^{0})^{T} (y^{0} − v^{0}) + x̄^{T} (y^{0} − v^{0}) + ȳ^{T} (x^{0} − u^{0})
≤ ϕ_{0} + x̄^{T} ȳ + (x^{0} − u^{0})^{T} (y^{0} − v^{0}) + x̄^{T} (y^{0} − v^{0}) + ȳ^{T} (x^{0} − u^{0})
=^{△} C̄_{2}.

Hence

0 < y_{i}^{k} ≤ \frac{C̄_{2}}{x̄_{i}} ,  0 < x_{i}^{k} ≤ \frac{C̄_{2}}{ȳ_{i}} ,  i = 1, 2, ⋯, n.

The result is obtained by setting

C_{2} = C̄_{2} max ( sup_{i=1,···,n} \frac{1}{x̄_{i}} , sup_{i=1,···,n} \frac{1}{ȳ_{i}} ) .

We can use Lemma 3.5 to define lower bounds on some other components of (x^{k}, y^{k}).

Lemma 3.6 For all k ≥ 0,

i ∈ B ⇒ y_{i}^{k} ≥ \frac{1}{C_{2}} γ̄μ_{k},
i ∈ N ⇒ x_{i}^{k} ≥ \frac{1}{C_{2}} γ̄μ_{k},

where C_{2} is as defined in Lemma 3.5.

Proof. We have from Lemma 3.5 that

i ∈ B, x_{i}^{k}y_{i}^{k} ≥ γ_{k}μ_{k} ⇒ y_{i}^{k} ≥ \frac{γ_{k}μ_{k}}{x_{i}^{k}} ≥ \frac{γ_{k}μ_{k}}{C_{2}} ≥ \frac{γ̄μ_{k}}{C_{2}} .

4 Linear Convergence and Polynomial Complexity

In this section, we modify some results of Zhang [12] to show that the algorithm of Section 2 produces a sequence {ϕ_{k}} that converges Q-linearly. When the starting point is chosen appropriately, the method has polynomial complexity.

We start with a result that can be used to derive global bounds on the step (Δx^{k}, Δy^{k}). We define

ξ_{k} = (\frac{n}{γ_{k}} )^{1/2} \frac{(x^{k} − u^{k})^{T} y^{k} + (y^{k} − v^{k})^{T} x^{k}}{(x^{k})^{T} y^{k}}
η_{k} = (1 − 2σ_{k} + \frac{σ_{k}^{2}}{γ_{k}} ) + 2ν_{k} \frac{(x^{0} − u^{0})^{T} (y^{0} − v^{0})}{β̂(x^{0})^{T} y^{0}}
ω_{k} = (ξ_{k} + \sqrt{ξ_{k}^{2} + η_{k}})^{2} ,

and the diagonal matrix

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

which also appears in much of the subsequent analysis.

Lemma 4.1 For all k,

∥D^{k}Δx^{k}∥^{2} + ∥(D^{k})^{−1}Δy^{k}∥^{2} ≤ ω_{k}(x^{k})^{T} y^{k}.  (12)

Moreover, there is a constant ω > 0

ω_{k} ≤ ω.  (13)

Proof. Minor modifications of the proofs of Lemma 6.2 and Theorem 7.1 in Zhang [12] yield (12). Since, from Lemma 3.2, we have

(x^{k})^{T} y^{k} ≥ β̂ν_{k}(x^{0})^{T} y^{0},

we can modify the proof of Zhang [12, Lemma 6.1] to show that

ξ_{k} ≤ (\frac{n}{γ_{k}} )^{1/2} [1 + \frac{(x^{0} − u^{0})^{T} y^{∗} + (y^{0} − v^{0})^{T} x^{∗} + (x^{0} − u^{0})^{T} (y^{0} − v^{0})}{β̂(x^{0})^{T} y^{0}}] .

Since σ_{k} ∈ [0, 1) and γ_{k} ≥ γ̄ > 0, it is easy to see that {ξ_{k}} and {η_{k}} are bounded sequences. Hence {ω_{k}} is also bounded, and so we have (13).

We can now prove linear convergence.

Theorem 4.2 There is a constant δ ∈ (0, 1) such that

ϕ_{k+1} ≤ (1 − δ)ϕ_{k},  k = 0, 1, 2, ⋯,  (14)

that is, the algorithm converges globally and Q-linearly.

Proof. Consider a safe step with 0 < σ ≤ σ_{k} ≤ 1/2. As in the proof of Theorem 7.1 of Zhang [12], we can show that

ϕ_{k+1} ≤ (1 − δ_{k})ϕ_{k},

where

δ_{k} ≥ (1 − \frac{2(1 − γ_{k})σ}{n} ) \frac{(1 − γ_{k})σ}{nω} ≥ (1 − \frac{2(1 − γ̄)σ}{n} ) \frac{(1 − 2γ̄)σ}{nω} > 0,

where the second inequality follows from γ_{k} ∈ (γ̄, 2γ̄].

When a successful fast step is taken, we have by definition that

ϕ_{k+1} ≤ ρϕ_{k}.

The result follows by setting

δ = min ((1 − \frac{2(1 − γ̄)σ}{n} ) \frac{(1 − 2γ̄)σ}{nω} , 1 − ρ) .  (15)

The complexity result depends on a particular choice of starting point. Zhang [12] suggests the following choice. First, define

(u^{0}, v^{0}) = arg min_{(u,v)} ∥u∥^{2} + ∥v∥^{2}, subject to v = Mu + h,  (16)

and

r^{∗} = inf {∥(x^{∗}, y^{∗})∥ | (x^{∗}, y^{∗}) is a solution of (1) } .  (17)

Now, we choose r ≥ ∥(u^{0}, v^{0})∥ and set

(x^{0}, y^{0}) = r(e, e).  (18)

Theorem 4.3 Suppose that (x^{0}, y^{0}) is defined by (16)–(18) and that there is a constant β independent of n such that

r ≥ r^{∗}/(β\sqrt{n}).

Suppose that, for a given ϵ > 0, we have

ϕ_{0} ≤ 1/ϵ^{τ},

where τ is a constant independent of n. Then there is an integer K_{ϵ} such that

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

and ϕ_{k} ≤ ϵ for k ≥ K_{ϵ}.

Proof. As in Zhang [12, Lemma 7.1], with minor modifications, it can be shown that if we choose

ω = lim sup_{k→∞} ω_{k},

then ω = O(n). Equation (15) then implies that δ ≥ δ̄/n^{2}, for some δ̄ > 0 independent of n. Therefore ϕ_{k} ≤ ϵ when

k ≥ \frac{log(ϵ/ϕ_{0})}{log(1 − δ)} = \frac{O(log(1/ϵ))}{δ}.

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

5 Boundedness of the Steps

In this section, we obtain bounds for components of the steps Δx^{k} and Δy^{k} in terms of γ_{k}, σ_{k}, and μ_{k}. These bounds are finer than the global bounds implied by (12) and are important for the analysis of superlinear convergence. We treat components in the basic and nonbasic index sets B and N differently and define

Δx_{B}^{k} = [Δx_{i}^{k}]_{i∈B}, Δx_{N}^{k} = [Δx_{i}^{k}]_{i∈N},

and similarly for Δy_{B}^{k}, Δy_{N}^{k}, etc. Without loss of generality, we can assume that the components of all vectors (and the rows and columns of M) are rearranged so that the indices in B are listed first. We can then partition M and D^{k} in an obvious way as

M = [ \begin{matrix}M_{BB} & M_{BN} \\ M_{NB} & M_{NN} \end{matrix} ] , D^{k} = [ \begin{matrix}D_{B}^{k} \\ & D_{N}^{k} \end{matrix} ] .

(X^{k} and Y^{k} may be partitioned similarly to D^{k}.) Here, e_{B} and e_{N} denote vectors of the appropriate length whose components are all 1.

We start with a simple result that bounds components of Δx_{N}^{k} and Δy_{B}^{k}.

Lemma 5.1 There is a constant C_{4} such that

i ∈ N ⇒ |Δx_{i}^{k}| ≤ C_{4}μ_{k},  (19a)
i ∈ B ⇒ |Δy_{i}^{k}| ≤ C_{4}μ_{k}.  (19b)

Proof. From (12) and (13), we have for i ∈ N that

|Δx_{i}^{k} (\frac{y_{i}^{k}}{x_{i}^{k}} )^{1/2}| ≤ (nμ_{k}ω)^{1/2},

and so

|Δx_{i}^{k}| ≤ (nω)^{1/2} \frac{μ_{k}^{1/2}}{(x_{i}^{k}y_{i}^{k})^{1/2}} x_{i}^{k} ≤ \frac{(nω)^{1/2}}{γ_{k}^{1/2}} x_{i}^{k}.

If we define C_{4} = C_{1}(nω/γ̄)^{1/2}, the inequality (19a) follows from Lemma 3.4. The proof of (19b) is similar.

Bounds on the remaining components Δx_{B}^{k} and Δy_{N}^{k} are more difficult to obtain. We need the following lemma, which is similar to results obtained by Ye and Anstreicher [11, Lemma 3.5] and Monteiro and Wright [5, Lemma 3.6]. For clarity, we drop the superscript k from vectors and matrices in the following two results.

Lemma 5.2 The vector pair (w, z) = (Δx_{B}, Δy_{N}) solves the convex quadratic programming problem

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

subject to

M_{BB}w = (−h − Mx + y)_{B} − M_{BN}Δx_{N} + Δy_{B},  (21a)
M_{NB}w − z = (−h − Mx + y)_{N} − M_{NN}Δx_{N}.  (21b)

Proof. The first-order necessary conditions for (w, z) to solve (20),(21) are also sufficient, since the problem is convex. These conditions are that (w, z) is feasible with respect to (21) and

[ \begin{matrix}D_{B}^{2}w − σ_{k}μ_{k}X_{B}^{−1}e_{B} \\ D_{N}^{−2}z − σ_{k}μ_{k}Y_{N}^{−1}e_{N} \end{matrix} ] ∈ R ([ \begin{matrix}M_{BB}^{T} & M_{NB}^{T} \\ 0 & −I \end{matrix} ]) ,  (22)

where R(.) denotes the range space of a matrix. Ye and Anstreicher [11, Lemma 3.4] prove that

R ([ \begin{matrix}M_{BB}^{T} & M_{NB}^{T} \\ 0 & −I \end{matrix} ]) = R ([ \begin{matrix}M_{BB} & M_{BN} \\ 0 & I \end{matrix} ]) ,

and so (22) becomes

[ \begin{matrix}D_{B}^{2}w − σ_{k}μ_{k}X_{B}^{−1}e_{B} \\ D_{N}^{−2}z − σ_{k}μ_{k}Y_{N}^{−1}e_{N} \end{matrix} ] ∈ R ([ \begin{matrix}M_{BB} & M_{BN} \\ 0 & I \end{matrix} ]) .  (23)

We now show that (Δx_{B}, Δy_{N}) satisfies (23). Equation (2), appropriately partitioned and scaled, yields

M_{BB}Δx_{B} = (−h − Mx + y)_{B} − M_{BN}Δx_{N} + Δy_{B},  (24a)
M_{NB}Δx_{B} − Δy_{N} = (−h − Mx + y)_{N} − M_{NN}Δx_{N},  (24b)
D_{B}^{2}Δx_{B} + Δy_{B} = −y_{B} + σ_{k}μ_{k}X_{B}^{−1}e_{B},  (24c)
D_{N}^{−2}Δy_{N} + Δx_{N} = −x_{N} + σ_{k}μ_{k}Y_{N}^{−1}e_{N}.  (24d)

It follows from (24a) and (24b) that (w, z) satisfies the constraints (21). Now, since

y^{∗} = Mx^{∗} + h ⇒ 0 = M_{BB}x_{B}^{∗} + h_{B},

we have by substitution of (24c) in (24a) that

D_{B}^{2}Δx_{B} − σ_{k}μ_{k}X_{B}^{−1}e_{B}
= −y_{B} − [M_{BB}Δx_{B} + h_{B} + M_{BB}x_{B} + M_{BN}x_{N} − y_{B} + M_{BN}Δx_{N}]
= −M_{BB}(x_{B} + Δx_{B} − x_{B}^{∗}) − M_{BN}(x_{N} + Δx_{N}).  (25)

From (24d),

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

Together, (25) and (26) imply (23).

Lemma 5.3 There is a constant C_{5} > 0 such that

∥Δx_{B}∥ ≤ \frac{C_{5}}{2} (μ_{k} + σ_{k}).  (27a)
∥Δy_{N}∥ ≤ \frac{C_{5}}{2} (μ_{k} + σ_{k}).  (27b)

Proof. Since the feasible set for (21) is nonempty, there is a feasible vector pair (w̄, z̄) such that

∥[ w̄ ; z̄ ]∥ = O (∥ − h − Mx + y∥) + O (∥Δx_{N}∥) + O (∥Δy_{B}∥) = O(μ_{k}),

where the estimate for the infeasibility term follows from Lemma 3.2, since

∥ − h − Mx^{k} + y^{k}∥ = ν_{k}∥ − h − Mx^{0} + y^{0}∥ ≤ \frac{μ_{k}}{β̂μ_{0}} ∥ − h − Mx^{0} + y^{0}∥.

Note that the objective function (20) may be written as

\frac{1}{2}∥D_{B}w − σ_{k}μ_{k}D_{B}^{−1}X_{B}^{−1}e_{B}∥^{2} + \frac{1}{2}∥D_{N}^{−1}z − σ_{k}μ_{k}D_{N}Y_{N}^{−1}e_{N}∥^{2} + constant.

Hence

∥Δx_{B} − σ_{k}μ_{k}D_{B}^{−2}X_{B}^{−1}e_{B}∥^{2} + ∥Δy_{N} − σ_{k}μ_{k}D_{N}^{2}Y_{N}^{−1}e_{N}∥^{2}
≤ ∥D_{B}^{−1}∥^{2}∥D_{B}Δx_{B} − σ_{k}μ_{k}D_{B}^{−1}X_{B}^{−1}e_{B}∥^{2} + ∥D_{N}∥^{2}∥D_{N}^{−1}Δy_{N} − σ_{k}μ_{k}D_{N}Y_{N}^{−1}e_{N}∥^{2}
≤ C̄_{5}^{2} {∥D_{B}Δx_{B} − σ_{k}μ_{k}D_{B}^{−1}X_{B}^{−1}e_{B}∥^{2} + ∥D_{N}^{−1}Δy_{N} − σ_{k}μ_{k}D_{N}Y_{N}^{−1}e_{N}∥^{2}}
≤ C̄_{5}^{2} {∥D_{B}w̄ − σ_{k}μ_{k}D_{B}^{−1}X_{B}^{−1}e_{B}∥^{2} + ∥D_{N}^{−1}z̄ − σ_{k}μ_{k}D_{N}Y_{N}^{−1}e_{N}∥^{2}} ,

where C̄_{5} = max(∥D_{B}^{−1}∥, ∥D_{N}∥). The last inequality follows from optimality of (Δx_{B}, Δy_{N}) in (20),(21). Hence

∥[ \begin{matrix}Δx_{B} − σ_{k}μ_{k}D_{B}^{−2}X_{B}^{−1}e_{B} \\ Δy_{N} − σ_{k}μ_{k}D_{N}^{2}Y_{N}^{−1}e_{N} \end{matrix} ]∥ ≤ C̄_{5} ∥[ \begin{matrix}D_{B}w̄ − σ_{k}μ_{k}D_{B}^{−1}X_{B}^{−1}e_{B} \\ D_{N}^{−1}z̄ − σ_{k}μ_{k}D_{N}Y_{N}^{−1}e_{N} \end{matrix} ]∥ ,

and so

∥[ \begin{matrix}Δx_{B} \\ Δy_{N} \end{matrix} ]∥ ≤ C̄_{5} \{∥[ \begin{matrix}D_{B} \\ & D_{N}^{−1} \end{matrix} ]∥ ∥[ w̄ ; z̄ ]∥ + σ_{k}μ_{k} ∥[ \begin{matrix}D_{B}^{−1} \\ & D_{N} \end{matrix} ]∥ ∥[ \begin{matrix}X_{B}^{−1}e_{B} \\ Y_{N}^{−1}e_{N} \end{matrix} ]∥\}

+σ_{k}μ_{k} ∥[ \begin{matrix}D_{B}^{−2} \\ & D_{N}^{2} \end{matrix} ]∥ ∥[ \begin{matrix}X_{B}^{−1}e_{B} \\ Y_{N}^{−1}e_{N} \end{matrix} ]∥ .

Now

∥D_{B}^{−1}∥ = sup_{i∈B} (\frac{x_{i}^{k}}{y_{i}^{k}} )^{1/2} ≤ C_{2} \frac{1}{(γ̄μ_{k})^{1/2}} ,  ∥D_{N}∥ = C_{2} \frac{1}{(γ̄μ_{k})^{1/2}} ,  (Lemmas 3.5, 3.6),

∥D_{B}∥ = sup_{i∈B} (\frac{y_{i}^{k}}{x_{i}^{k}} )^{1/2} ≤ C_{1} (\frac{μ_{k}}{γ̄} )^{1/2} ,  ∥D_{N}^{−1}∥ ≤ C_{1} (\frac{μ_{k}}{γ̄} )^{1/2} ,  (Lemma 3.4),
∥X_{B}^{−1}e_{B}∥ = sup_{i∈B} \frac{1}{x_{i}^{k}} ≤ \frac{C_{1}}{γ̄} ,  ∥Y_{N}^{−1}e_{N}∥ ≤ \frac{C_{1}}{γ̄} ,  (Lemma 3.4),
∥w̄∥ = O(μ_{k}),  ∥z̄∥ = O(μ_{k}),  C̄_{5} = max (∥D_{B}^{−1}∥, ∥D_{N}∥) ≤ \frac{C_{2}}{(γ̄μ_{k})^{1/2}} .

Hence

∥[ \begin{matrix}Δx_{B} \\ Δy_{N} \end{matrix} ]∥ ≤ O(μ_{k}^{−1/2}) \{O(μ_{k}^{3/2}) + σ_{k}O(μ_{k}^{1/2})\} + σ_{k}μ_{k}O(μ_{k}^{−1})
≤ \frac{C_{5}}{2} (μ_{k} + σ_{k}),

for an appropriately defined constant C_{5}.

For “fast” steps, the results of Lemmas 5.1 and 5.3 can be combined to give the following result.

Theorem 5.4 If σ_{k} = μ_{k}, then

|Δx_{i}^{k}Δy_{i}^{k}| ≤ C_{4}C_{5}μ_{k}^{2}, i = 1, ⋯, n.

6 Asymptotic Subquadratic Convergence

In this section, we show that all steps after a certain point in the algorithm are fast steps and that the sequences {μ_{k}} and {ϕ_{k}} converge Q-subquadratically to zero.

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

C_{6} =^{△} max (1, 2C_{4}C_{5}) .  (28)

We start by defining a threshold condition involving μ_{k} and γ_{k}, and finding a lower bound on the step length when this condition is satisfied.

Lemma 6.1 Suppose at iteration k that

\frac{μ_{k}}{(γ_{k} − γ̄)(1 − γ̄)} ≤ \frac{ρ}{3C_{6}}  (29)

and that a fast step is calculated. Then the step length α̃ will satisfy

1 ≥ α̃ ≥ 1 − C_{6} \frac{μ_{k}}{(γ_{k} − γ̄)(1 − γ̄)}.

Proof. Before proceeding, note that if the fast step is successful, the algorithm sets γ_{k+1} to γ̄(1 + γ̄^{t_{k}}), and so

γ_{k} − γ_{k+1} = γ̄(1 + γ̄^{t_{k}−1}) − γ̄(1 + γ̄^{t_{k}}) = γ̄^{t_{k}}(1 − γ̄) = (γ_{k} − γ̄)(1 − γ̄).  (30)

Under these circumstances, condition (29) is equivalent to

\frac{μ_{k}}{γ_{k} − γ_{k+1}} ≤ \frac{ρ}{3C_{6}} .  (31)

We use (29) and (31) interchangeably for the remainder of the proof.

The proof is in three stages. First, we show that the relaxed centrality condition (4e) holds for all α satisfying

α ∈ [0, 1 − C_{6} \frac{μ_{k}}{γ_{k} − γ_{k+1}} ] .  (32)

Second, we show that the other main condition on α, namely, (4d), also holds for α satisfying (32). Together, (4d) and (4e) imply that (4b) and (4c) also hold. Third, we show that ϕ(x^{k} + αΔx^{k}, y^{k} + αΔy^{k}) is decreasing for α in the range (32), which, by (3), implies the result.

We start with (4e). Theorem 5.4, (2), (4e), (28), and the fact that σ̃ = σ_{k} = μ_{k} imply that

|Δx_{i}^{k}Δy_{i}^{k}| ≤ C_{4}C_{5}μ_{k}^{2} ≤ (C_{6}/2)μ_{k}^{2},
x_{i}^{k}y_{i}^{k} ≥ γ_{k}μ_{k},
x_{i}^{k}Δy_{i}^{k} + y_{i}^{k}Δx_{i}^{k} = −x_{i}^{k}y_{i}^{k} + σ_{k}μ_{k} = −x_{i}^{k}y_{i}^{k} + μ_{k}^{2}.

Hence

(x_{i}^{k} + αΔx_{i}^{k})(y_{i}^{k} + αΔy_{i}^{k})
= x_{i}^{k}y_{i}^{k} + α(x_{i}^{k}Δy_{i}^{k} + y_{i}^{k}Δx_{i}^{k}) + α^{2}(Δx_{i}^{k})^{T} Δy_{i}^{k}
≥ γ_{k}μ_{k}(1 − α) + αμ_{k}^{2} − α^{2}(C_{6}/2)μ_{k}^{2}.

On the other hand,

\frac{1}{n}γ_{k+1}(x^{k} + αΔx_{k})^{T} (y^{k} + αΔy_{k})

≤ γ_{k+1} [μ_{k}(1 − α + αμ_{k}) + α^{2}(C_{6}/2)μ_{k}^{2}] .

Therefore condition (4e) will be satisfied provided that

γ_{k}μ_{k}(1 − α) + αμ_{k}^{2} − α^{2}(C_{6}/2)μ_{k}^{2}
≥ γ_{k+1} [μ_{k}(1 − α + αμ_{k}) + α^{2}(C_{6}/2)μ_{k}^{2}] ,

or, equivalently,

(γ_{k} − γ_{k+1})μ_{k}(1 − α) + αμ_{k}^{2}(1 − γ_{k+1}) − α^{2}μ_{k}^{2}(C_{6}/2)(1 + γ_{k+1}) ≥ 0.

Since γ_{k+1} ∈ (0, 1) and α ∈ (0, 1], this last inequality will hold if

(γ_{k} − γ_{k+1})μ_{k}(1 − α) − αC_{6}μ_{k}^{2} ≥ 0 ⇒ 0 ≤ α ≤ ᾱ_{I} =^{△} \frac{1}{1 + C_{6} \frac{μ_{k}}{γ_{k}−γ_{k+1}}} .

By assumption (31), the second term in the denominator is less than one, so we can write

ᾱ_{I} ≥ 1 − C_{6} \frac{μ_{k}}{γ_{k} − γ_{k+1}} .

Therefore, (4e) is satisfied for α in the range (32).

Turning now to condition (4d), we note that

\frac{1}{n}(x^{k} + αΔx^{k})^{T} (y^{k} + αΔy^{k}) = μ_{k}(1 − α + αμ_{k}) + \frac{1}{n}α^{2}(Δx^{k})^{T} Δy^{k}
≥ μ_{k}(1 − α + αμ_{k}) − α(C_{6}/2)μ_{k}^{2},

where we have used Theorem 5.4, (28), and α ∈ (0, 1] to derive the inequality. Therefore, (4d) will hold if

μ_{k}(1 − α + αμ_{k}) − α(C_{6}/2)μ_{k}^{2} ≥ (1 − α)(1 − γ̄^{t_{k}})μ_{k}
⇔ γ̄^{t_{k}} − α [γ̄^{t_{k}} − μ_{k} + (C_{6}/2)μ_{k}] ≥ 0.  (33)

From (28), (29), and (30), we have that

\frac{μ_{k}}{γ̄^{t_{k}}} ≤ \frac{ρ}{3C_{6}} (1 − γ̄) < 1.

Hence the term in square brackets in (33) is positive, and so (4d) is satisfied if

0 ≤ α ≤ ᾱ_{II} =^{△} \frac{γ̄^{t_{k}}}{γ̄^{t_{k}} − μ_{k} + (C_{6}/2)μ_{k}} .

Now

ᾱ_{II} ≥ \frac{1}{1 + (C_{6}/2)μ_{k}/γ̄^{t_{k}}} ≥ 1 − (C_{6}/2)\frac{μ_{k}}{γ̄^{t_{k}}} .

Using (30) again, we have that γ̄^{t_{k}} ≥ (γ_{k} − γ_{k+1}), and so

ᾱ_{II} ≥ 1 − C_{6} \frac{μ_{k}}{γ_{k} − γ_{k+1}} .

Therefore (4d) holds for α in the range (32).

Finally, we examine ϕ(x^{k} + αΔx^{k}, y^{k} + αΔy^{k}) for α ∈ [0, 1]. Using the notation r_{k} = ∥y^{k} − Mx^{k} − h∥, we have that

ϕ(x^{k} + αΔx^{k}, y^{k} + αΔy^{k}) = (x^{k} + αΔx^{k})^{T} (y^{k} + αΔy^{k}) + (1 − α)r_{k}
= nμ_{k} + α[(x^{k})^{T} Δy^{k} + (y^{k})^{T} Δx^{k}] + α^{2}(Δx^{k})^{T} Δy^{k} + (1 − α)r_{k}
= nμ_{k} − αn(1 − μ_{k})μ_{k} + α^{2}(Δx^{k})^{T} Δy^{k} + (1 − α)r_{k}.

If (Δx^{k})^{T} Δy^{k} ≤ 0, then, since μ_{k} < 1 by (31), ϕ(x^{k} + αΔx^{k}, y^{k} + αΔy^{k}) is decreasing on [0, 1]. Otherwise, the unconstrained minimum occurs at

α_{min} = \frac{n(1 − μ_{k})μ_{k} + r_{k}}{2(Δx^{k})^{T} Δy^{k}} ≥ \frac{(2/3)nμ_{k}}{2(Δx^{k})^{T} Δy^{k}} ≥ \frac{(2/3)nμ_{k}}{2nC_{4}C_{5}μ_{k}^{2}} ≥ \frac{(2/3)}{C_{6}μ_{k}}.  (34)

But from (31),

(3/2)C_{6}μ_{k} ≤ \frac{ρ}{2}(γ_{k} − γ_{k+1}) ≤ 1,  (35)

and so α_{min} ≥ 1. Again, we deduce that ϕ(x^{k} + αΔx^{k}, y^{k} + αΔy^{k}) is decreasing on [0, 1]. This implies that the value of α̃ chosen by the algorithm is the largest value that satisfies (4), and so

1 ≥ α̃ ≥ 1 − C_{6} \frac{μ_{k}}{γ_{k} − γ_{k+1}} ,

as required.

The next two results show that once the threshold condition (29) is satisfied at some iteration, then fast steps may be taken on this and all subsequent iterations, and subquadratic convergence can be attained.

Lemma 6.2 If (29) is satisfied at iteration k and a fast step is taken, then

μ_{k+1} ≤ ρμ_{k},  (36a)
ϕ_{k+1} ≤ ρϕ_{k},  (36b)

μ_{k+1} ≤ (\frac{3C_{6}}{γ_{k} − γ_{k+1}} ) μ_{k}^{2},  (37a)
ϕ_{k+1} ≤ (\frac{3C_{6}}{n(γ_{k} − γ_{k+1})}) ϕ_{k}^{2},  (37b)

and

\frac{μ_{k+1}}{(γ_{k+1} − γ̄)(1 − γ̄)} ≤ \frac{μ_{k}}{(γ_{k} − γ̄)(1 − γ̄)}.  (38)

Proof. Consider

\frac{1}{n}(x^{k} + αΔx^{k})^{T} (y^{k} + αΔy^{k}) = μ_{k}[1 − α(1 − μ_{k})] + α^{2} \frac{(Δx^{k})^{T} Δy^{k}}{n} .  (39)

The argument in the last part of Lemma 6.1, in particular, formulae (34) and (35), can be applied to show that the quadratic function (39) is decreasing on [0, 1]. Therefore

α_{k} ≥ 1 − C_{6} \frac{μ_{k}}{γ_{k} − γ_{k+1}} ,

and so

μ_{k+1} = \frac{1}{n}(x^{k} + α_{k}Δx^{k})^{T} (y^{k} + α_{k}Δy^{k})
≤ μ_{k} [1 − (1 − C_{6} \frac{μ_{k}}{γ_{k} − γ_{k+1}} ) (1 − μ_{k})] + \frac{1}{n} |(Δx^{k})^{T} Δy^{k}|
≤ μ_{k} [C_{6} \frac{μ_{k}}{γ_{k} − γ_{k+1}} + μ_{k}] + (C_{6}/2)μ_{k}^{2}  (40)
≤ μ_{k} [\frac{ρ}{3} + \frac{ρ}{3}] + \frac{ρ}{3}μ_{k}
≤ ρμ_{k},

yielding (36a). For (36b), we again use the notation r_{k} = ∥y^{k} − Mx^{k} − h∥ and note that

r_{k+1} = (1 − α_{k})r_{k} ≤ C_{6} \frac{μ_{k}}{γ_{k} − γ_{k+1}} r_{k} ≤ ρr_{k}.

Therefore,

ϕ_{k+1} = nμ_{k+1} + r_{k+1} ≤ nρμ_{k} + ρr_{k} = ρϕ_{k}.

For (37a), using (40) and the fact that

1 ≤ C_{6} ≤ C_{6} \frac{1}{γ_{k} − γ_{k+1}} ,

we have that

μ_{k+1} ≤ μ_{k} [C_{6} \frac{μ_{k}}{γ_{k} − γ_{k+1}} + μ_{k} + C_{6}μ_{k}] ≤ \frac{3C_{6}}{γ_{k} − γ_{k+1}} μ_{k}^{2},

as required. Also,

ϕ_{k+1} = nμ_{k+1} + r_{k+1}
≤ \frac{3C_{6}}{γ_{k} − γ_{k+1}} nμ_{k}^{2} + \frac{C_{6}}{γ_{k} − γ_{k+1}} μ_{k}r_{k}
≤ \frac{3C_{6}}{γ_{k} − γ_{k+1}} μ_{k}[nμ_{k} + r_{k}]
≤ \frac{3C_{6}}{n(γ_{k} − γ_{k+1})}[nμ_{k} + r_{k}]^{2}
= \frac{3C_{6}}{n(γ_{k} − γ_{k+1})}ϕ_{k}^{2},

giving (37b).

From (37a) and (30), we have

μ_{k+1} ≤ \frac{3C_{6}}{γ̄^{t_{k}}(1 − γ̄)}μ_{k}^{2}.

Therefore,

\frac{μ_{k+1}}{(γ_{k+1} − γ̄)(1 − γ̄)} = \frac{μ_{k+1}}{γ̄^{t_{k}+1}(1 − γ̄)} ≤ \frac{1}{γ̄^{t_{k}+1}(1 − γ̄)} \frac{3C_{6}}{γ̄^{t_{k}}(1 − γ̄)}μ_{k}^{2} = \frac{3C_{6}}{γ̄} [\frac{μ_{k}}{γ̄^{t_{k}}(1 − γ̄)}]^{2} . (41)

From (29) and (30), we have

\frac{3C_{6}}{γ̄} \frac{μ_{k}}{γ̄^{t_{k}}(1 − γ̄)} = \frac{3C_{6}}{γ̄} \frac{μ_{k}}{(γ_{k} − γ̄)(1 − γ̄)} ≤ \frac{ρ}{γ̄} < 1,

where the last inequality follows from the definition of ρ and γ̄ in the algorithm. Substitution of this inequality into (41) gives (38).

Theorem 6.3 Suppose that the condition (29) is satisfied at iteration K and that

ϕ_{K} ≤ ϕ̄.

Then (i) the algorithm will take fast steps at iteration K and at all subsequent iterations, and

(ii) the sequences {μ_{k}} and {ϕ_{k}} converge superlinearly to zero with Q-order 2.

Proof. The first part follows inductively from Lemma 6.2. Since (29) is satisfied at iteration K and since (36b) holds, the fast step will be accepted at this iteration. Formula (38) implies that (29) will again hold at iteration K + 1. Since ϕ_{K+1} ≤ ρϕ_{K} < ϕ̄, a fast step will be accepted at iteration K + 1, and so on.

For (ii), we use (37a). Since

μ_{k+1} ≤ \frac{3C_{6}}{γ̄^{t_{k}}(1 − γ̄)}μ_{k}^{2},

we have

log μ_{k+1} ≤ 2 log μ_{k} + log(3C_{6}) − t_{k} log γ̄ − log(1 − γ̄).  (42)

We can now use an argument like that of Ye [10]. From (36a) and

log ρ < log γ̄ < 0  and  1 ≤ t_{k} ≤ k + 1,

we have for sufficiently large choice of k that

log μ_{k} ≤ log μ_{K} + (k − K) log ρ ≪ (k + 1) log γ̄ ≤ t_{k} log γ̄,

that is, the first term on the right hand side of (42) will eventually dominate the third term and, in fact,

lim_{k→∞} \frac{t_{k}}{log μ_{k}} = 0.

Dividing (42) by log μ_{k}, we have

lim inf_{k→∞} \frac{log μ_{k+1}}{log μ_{k}} = 2.  (43)

From Potra [6], (43) implies that the Q-order of convergence for {μ_{k}} is 2.

The argument for {ϕ_{k}} is similar.

Finally, we show that the threshold condition (29) will eventually be met and, hence, that subquadratic convergence will be obtained.

Theorem 6.4 Define

f =^{△} \frac{log γ̄}{log ρ} ∈ (0, 1).

Define a constant ϵ̂ as follows:

ϵ̂ = \{ \begin{matrix}ϕ̄ & if \frac{ϕ̄}{nγ̄(1−γ̄)} ≤ \frac{ρ}{3C_{6}} , \\ [\frac{ρ}{3C_{6}} \frac{nγ̄^{2}(1−γ̄)}{ϕ̄}]1/(1−f) ϕ̄ & otherwise. \end{matrix}  (44)

Then if K is the smallest positive integer such that

ϕ_{K} ≤ ϵ̂,  (45)

we have that

\frac{μ_{K}}{γ̄^{t_{K}}(1 − γ̄)} ≤ \frac{ρ}{3C_{6}}  (46)

and

ϕ_{K} ≤ ϕ̄,  (47)

and hence the conditions of Theorem 6.3 are satisfied.

Proof. First, consider the case of

\frac{ϕ̄}{nγ̄(1 − γ̄)} ≤ \frac{ρ}{3C_{6}} .  (48)

Then ϵ̂ = ϕ̄, and so (45) immediately implies (47). Since K is the first iterate below ϕ̄, the algorithm cannot have taken any fast steps yet, so t_{K} = 1. Hence, by (45), (48), and the fact that μ_{K} ≤ ϕ_{K}/n,

\frac{μ_{K}}{γ̄^{t_{K}}(1 − γ̄)} = \frac{μ_{K}}{γ̄(1 − γ̄)} ≤ \frac{ϕ_{K}/n}{γ̄(1 − γ̄)} ≤ \frac{ϕ̄/n}{γ̄(1 − γ̄)} ≤ \frac{ρ}{3C_{6}} ,

which gives (46).

In the remaining case, note that

\frac{ρ}{3C_{6}} \frac{nγ̄^{2}(1 − γ̄)}{ϕ̄} ≤ \frac{ρ}{3C_{6}} \frac{nγ̄(1 − γ̄)}{ϕ̄} < 1,

and so ϵ̂ < ϕ̄ in (44), so again (45) implies (47). This inequality also implies that the algorithm may have taken some fast steps prior to iteration K but, since a reduction of at least ρ occurs on each fast step, the number of such steps is bounded by

⌈\frac{log(ϵ̂/ϕ̄)}{log ρ}⌉ .

Hence

t_{K} ≤ ⌈\frac{log(ϵ̂/ϕ̄)}{log ρ}⌉ + 1 ≤ \frac{log ϵ̂− log ϕ̄}{log ρ} + 2.

Now,

γ̄^{t_{K}} ≥ γ̄^{2} exp [\frac{log ϵ̂− log ϕ̄}{log ρ} log γ̄]
= γ̄^{2} exp [f log ϵ̂− f log ϕ̄]
= γ̄^{2}ϵ̂^{f} ϕ̄^{−f}.

Hence

\frac{μ_{K}}{γ̄^{t_{K}}(1 − γ̄)} ≤ \frac{μ_{K}}{γ̄^{2}(1 − γ̄)ϵ̂^{f} ϕ̄^{−f}}

≤ \frac{ϕ_{K}/n}{γ̄^{2}(1 − γ̄)ϵ̂^{f} ϕ̄^{−f}} ≤ \frac{ϵ̂}{nγ̄^{2}(1 − γ̄)ϵ̂^{f} ϕ̄^{−f}} = ϵ̂^{1−f} \frac{ϕ̄^{f}}{nγ̄^{2}(1 − γ̄)}.

Substituting in this expression from (44), we have

\frac{μ_{K}}{γ̄^{t_{K}}(1 − γ̄)} ≤ \frac{ρ}{3C_{6}} \frac{nγ̄^{2}(1 − γ̄)}{ϕ̄}ϕ̄^{1−f} \frac{ϕ̄^{f}}{nγ̄^{2}(1 − γ̄)} = \frac{ρ}{3C_{6}} ,

and so (46) is satisfied.

The following result is immediate from Theorems 6.3 and 6.4.

Corollary 6.5 The algorithm converges Q-subquadratically.

References

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

[2] M. Kojima, N. Meggido, and S. Mizuno, A primal-dual exterior point algorithm for linear programming, Tech. Rep. RJ 8500, IBM Alamaden Research Center, San Jose, Calif., 1991.

[3] M. Kojima, S. Mizuno, and M. Todd, Infeasible interior-point primal-dual potential-reduction algorithms for linear programming, Tech. Rep. 1023, School of Operations Research and Industrial Engineering, Cornell University, Ithaca, New York, August 1992.

[4] I. J. Lustig, R. E. Marsten, and D. F. Shanno, Computational experience with a globally convergent primal-dual predictor-corrector algorithm for linear programming, Tech. Rep. SOR 92–10, Program in Statistics and Operations Research, Department of Civil Engineering and Operations Research, Princeton University, Princeton, New Jersey, September 1992.

[5] R. D. C. Monteiro and S. J. Wright, A globally and superlinearly convergent potential reduction interior point method for convex programming, Tech. Rep. MCS– P316–0792, Mathematics and Computer Science Division, Argonne National Laboratory, Argonne, Illinois, August 1992.

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

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

[8] , A quadratically convergent infeasible interior-point algorithm for linear programming, Tech. Rep. 28, Department of Mathematics, University of Iowa, Iowa City, Iowa, July 1992.

[9] R. T. Rockafellar, Computational schemes for large-scale problems in extended linear-quadratic programming, Mathematical Programming, 48 (1990), pp. 447–474.

[10] Y. Ye, On the Q-order of convergence of interior-point algorithms for linear programming, Tech. Rep. 91–17, Department of Management Sciences, University of Iowa, October 1991.

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

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

[13] Y. Zhang and D. Zhang, Superlinear convergence of infeasible interior-point methods for linear programming, Tech. Rep. 92–15, Department of Mathematics and Statistics, University of Maryland, Baltimore County, September 1992.