A Large-step Infeasible-interior-point Method for the P_{∗}-matrix LCP

Florian A. Potra^{∗}, Rongqin Sheng^{†}

December, 1994

REPORTS ON COMPUTATIONAL MATHEMATICS, NO. 64/1994, DEPARTMENT OF MATHEMATICS, THE UNIVERSITY OF IOWA

Abstract

A large-step infeasible-interior-point method is proposed for solving P_{∗}(κ)-matrix linear complementarity problems. It is new even for monotone LCP. The algorithm generates points in a large neighborhood of an infeasible central path. Each iteration requires only one matrix factorization. If the problem is solvable, then the algorithm converges from an arbitrary positive starting point (x^{0}, s^{0}). The computational complexity of the algorithm depends on the quality of the starting point. If a well centered starting point is feasible or close to being feasible, then it has O((1 + κ)\sqrt{n}L)-iteration complexity. With appropriate initialization, a modified version of the algorithm terminates in O((1 + κ)^{2}nL) steps either by finding a solution or by determining that the problem is not solvable. High-order local convergence is proved for problems having a strictly complementary solution. We note that while the properties of the algorithm (e.g. computational complexity) depend on κ the algorithm itself does not.

Key Words: linear complementarity problems, P_{∗}-matrices, path-following, interior-point algorithm, polynomiality, superquadratic convergence. Abbreviated Title: A large-step method for LCP.

∗Department of Mathematics, University of Iowa, Iowa City, IA 52242, USA. The work of this author was supported in part by NSF, Grant DMS 9305760.

†Department of Mathematics, University of Iowa, Iowa City, IA 52242, USA.

1 Introduction

The linear complementarity problem (LCP) consists in determining a vector pair (x, s) ∈ R I ^{2n} that satisfies the conditions

s = Mx + q,  x ≥ 0,  s ≥ 0,  x^{T} s = 0,  (1.1)

where q ∈ R I ^{n} and M ∈ R I ^{n×n}. In this paper we consider problem (1.1) with M a P_{∗}- matrix. The class of P_{∗}-matrices was introduced by Kojima, Megiddo, Noma and Yoshise [6] and it contains many types of matrices encountered in practical applications. Let κ be a nonnegative number. A matrix M is called a P_{∗}(κ)−matrix if

(1 + 4κ) ∑_{i∈I_{+}(x)} x_{i}[Mx]_{i} + ∑_{i∈I_{−}(x)} x_{i}[Mx]_{i} ≥ 0 , ∀x ∈ R I ^{n} ,  (1.2)

where

I_{+}(x) = {i : x_{i}[Mx]_{i} > 0} , I_{−}(x) = {i : x_{i}[Mx]_{i} < 0} ,

or equivalently if

x^{T} Mx ≥ −4κ ∑_{i∈I_{+}(x)} x_{i}[Mx]_{i} , ∀x ∈ R I ^{n} .  (1.3)

The class of all P_{∗}(κ)−matrices is denoted by P_{∗}(κ), and the class P_{∗} is defined by

P_{∗} = ⋃_{κ≥0} P_{∗}(κ)

i.e., M is a P_{∗}-matrix if M ∈ P_{∗}(κ) for some κ ≥ 0.

Obviously, P_{∗}(0) = P SD (the class of positive semi-definite matrices) and P_{∗}(κ_{1}) ⊂ P_{∗}(κ_{2}) for 0 ≤ κ_{1} ≤ κ_{2}. Also we have P_{∗} ⊃ P, where P is the the class of all matrices with positive principal minors. This follows from the fact that a P -matrix M is a P_{∗}(κ)-matrix for κ = max \{−\frac{λ_{min}(M)}{4γ(M)}, 0\}, where λ_{min}(M) is the least eigenvalue of (M + M^{T})/2, and γ(M) > 0 is the so-called P -matrix number of M (see [6, Lemma 3.3]).

Most interior-point methods for linear programming problems have been extended for monotone linear complementarity problems, i.e., for the P_{∗}(0)-matrix LCP. Some of them have also been extended to the P_{∗}-matrix LCP (cf. [4], [5], [6], [7]).

In [6] and [7] it is assumed that the starting point (x^{0}, s^{0}) is strictly feasible, in the sense that it belongs to the relative interior of the feasible set

F = {(x, s) ∈ R I _{+}^{2n} : s = Mx + q}.

Such a starting point may be very difficult to find in practical applications. Moreover the existence of a strictly feasible starting point implies that the solution set

F^{∗} = {(x^{∗}, s^{∗}) ∈ F : x^{∗}^{T} s^{∗} = 0}

is nonempty and bounded, which restricts the class of problems to which the methods apply.

We note that all known infeasible-interior-point algorithms for the P_{∗}-matrix LCP depend on the classification number κ for their implementation. However, sometime it is extremely difficult to estimate κ. In this paper, we propose a large-step path-following algorithm independent of κ. If the problem is solvable, the algorithm is globally convergent when starting from an arbitrary positive starting point (x^{0}, s^{0}). The algorithm is defined in a large neighborhood of the central path and requires one matrix factorization per iteration. If a well centered starting point is feasible or close to being feasible, then the algorithm has O((1 + κ)\sqrt{n}L)-iteration complexity. If the starting point is large enough, then the iteration complexity is O((1 + κ)^{2}nL). With appropriate initialization a slightly modified version of the algorithm terminates in O((1 + κ)^{2}nL) steps either by finding a solution or by determining that the problem is not solvable. High-order local convergence is proved for problems having a strictly complementary solution.

We mention that our algorithm is new even for monotone LCP, where it attains O(nL)- iteration complexity for infeasible starting points, compared to O(n^{2}L)-iteration complexity of the algorithms, with similar work per iteration, proposed by Wright [16, 14, 15], Wright and Zhang [17], and Zhang [18]. In this paper, we have adapted the algorithm structure of “fast-safe-improve” of Wright and Zhang [17] by using new search directions and new neighborhood.

2 A Generic Large-step Infeasible-interior-point Algorithm

Our convergence results are proved under the assumption that F^{∗} is not empty. It is easily seen that (x^{∗}, s^{∗}) ∈ F^{∗} if and only if (x^{∗}, s^{∗}) is a nonnegative solution of the following nonlinear system

F (x, s) := ( \begin{matrix}Xs \\ Mx − s + q \end{matrix} ) = 0.  (2.1)

For any given ϵ > 0 we define the set of ϵ-approximate solutions of (1.1) as

F_{ϵ} = {(x, s) ∈ R I _{++}^{2n}, x^{T} s ≤ ϵ, ∥Mx − s + q∥ ≤ ϵ}.

In what follows we present an algorithm that finds a point in this set in a finite number of steps. The starting point of the algorithm can be any pair of strictly positive vectors (x^{0}, s^{0}) > 0. First, let us define

p = \frac{n}{(x^{0})^{T} s^{0}} X^{0}s^{0}, ρ_{0} = min \{[p]_{i} : 1 ≤ i ≤ n\} , τ_{0} = (x^{0})^{T} s^{0}/n.

Obviously, we have

ρ_{0} ≤ 1, ∑_{i=1}^{n}[p]_{i} = n.

The algorithm depends on seven positive parameters:

0 < α_{min} < α_{max} ≤ 1,  (2.2a)
0 < β_{min} < β_{max},  (2.2b)
0 < ρ < γ < 1,  ρ < \overline{ρ},  (2.2c)

and a nonnegative integer I. The complexity of the algorithm will not depend on \overline{ρ} and I. However the asymptotic order of convergence depends on I. Note that these parameters are independent of κ and n. We define a neighborhood

N = {(x, s, τ, α, β) ∈ R I _{++}^{2n+3} : Xs ≥ ατp, ∥Xs − τ p∥ ≤ βτ,
r = (τ/τ_{0})r^{0}, τ ≤ τ_{0}, α_{min} ≤ α ≤ α_{max}, β_{min} ≤ β ≤ β_{max}}

and denote

μ = \frac{1}{n}x^{T} s, r = s − Mx − q.

We note that β ∈ (0, 1) is a typical requirement of short-step algorithms. In our algorithm β does no longer have such a restriction. Therefore N is a large neighborhood of the infeasible central path defined by

Xs = τp,  (2.3a)
s = Mx + q + (τ/τ_{0})r^{0}.  (2.3b)

Also, α, β may change at each step. In our algorithm τ is driven to zero in a specified manner. Since all the points (x, s, τ, α, β) are in N, we have

α_{min}τ ≤ μ ≤ (1 + β_{max})τ,  (2.4)

so that μ will be driven to zero at about the same rate as τ. On the other hand ∥r∥ is driven to zero at exactly the same rate as τ.

At the beginning of the k-th iteration of our algorithm, a point (x^{k}, s^{k}, τ_{k}, α_{k}, β_{k}) ∈ N is given. The k-th iteration consists of several steps. Each step has input (x, s, τ, α, β) ∈ N and output (x^{+}, s^{+}, τ_{+}, α_{+}, β_{+}) ∈ N . Such a step is defined as follows:

• Choose λ ∈ [0, 1) (λ = 0 for fast step, λ ∈ (0, 1) for safe step).

• Solve the linear system

S^{k}u + X^{k}v = λτp − Xs,  (2.5a)
Mu − v = (1 − λ)r.  (2.5b)

• Denote

x(θ) = x + θu, s(θ) = s + θv.

• Choose α_{+}, β_{+} such that α_{min} ≤ α_{+}, β ≤ β_{+} ≤ β_{max}.

• Compute \overline{θ} ∈ [0, 1] as the largest step length for which

X(θ)s(θ) ≥ α_{+}(1 − (1 − λ)θ)τp,

and

∥X(θ)s(θ) − (1 − (1 − λ)θ)τ p∥ ≤ β_{+}(1 − (1 − λ)θ)τ,
for all 0 ≤ θ ≤ \overline{θ}.  (2.6)

• Set x^{+} = x + \overline{θ}u, s^{+} = s + \overline{θ}v, τ_{+} = (1 − (1 − λ)\overline{θ})τ.

Let us note that the computation of \overline{θ} above involves the solutions of n quadratic equations and a quartic equation. In the following analysis we assume that these equations are solved exactly. The results are true for appropriate approximate solutions of these equations by use of the bisection method (see also [2]).

Algorithm 2.1 Given α_{max}, α_{min}, β_{max}, β_{min}, γ, ρ satisfying (2.2a) and (x^{0}, s^{0}) > 0; t_{0} ← 0, α_{0} ← α_{max}, β_{0} ← β_{min}, τ_{0} ← μ_{0}; while max{μ_{k}, ∥r^{k}∥} > ϵ * fast branch * solve (2.5)-(2.6) with (x, s) = (x^{k}, s^{k}), λ = 0, α_{+} = α_{min} + γ^{t_{k}+1}(α_{max} − α_{min}), β_{+} = β_{max} − γ^{t_{k}+1}(β_{max} − β_{min}); if 1 − \overline{θ} ≤ ρ then t_{+} ← t_{k} + 1, λ_{k} ← 0;

else * safe branch * solve (2.5)-(2.6) with (x, s) = (x^{k}, s^{k}), λ ∈ (0, 1),

α_{+} = α_{min} + γ^{t_{k}}(α_{max} − α_{min}), β_{+} = β_{max} − γ^{t_{k}}(β_{max} − β_{min});

t_{+} ← t_{k}, λ_{k} ← λ;

end if

* main points * τ_{+} ← τ_{k}(1 − (1 − λ_{k})\overline{θ}), (x^{+}, s^{+}) ← (x^{k}, s^{k}) + \overline{θ}(u, v); improve (x^{k}, y^{k}, (x^{+}, s^{+}, τ_{+}, α_{+}, β_{+}, t_{+}));

* final points * (x^{k+1}, s^{k+1}, τ_{k+1}, α_{k+1}, β_{k+1}) ← (x^{+}, s^{+}, τ_{+}, α_{+}, β_{+}), k ← k + 1;

end while The improve procedure is defined as follows: Given \overline{ρ} ∈ (ρ, 1), I ≥ 0, for i = 1, 2, …, I

(x, s, τ, α, β, t) ← (x^{+}, s^{+}, τ_{+}, α_{+}, β_{+}, t_{+}) if max{μ_{k}, ∥r^{k}∥} ≤ ϵ, then return; solve (2.5)-(2.6) with λ = 0, α_{+} = α_{min} + γ^{t+1}(α_{max} − α_{min}), β_{+} = β_{max} − γ^{t+1}(β_{max} − β_{min}); if 1 − \overline{θ} ≤ ρ then t ← t + 1,

else t_{+} ← t, solve (2.5)-(2.6) with λ ∈ (0, 1), α_{+} = α_{min} + γ^{t}(α_{max} − α_{min}), β_{+} = β_{max} − γ^{t}(β_{max} − β_{min});

if 1 − (1 − λ)\overline{θ} > \overline{ρ} then return;

end if

* intermediary points * τ_{+} ← τ(1 − (1 − λ)\overline{θ}), (x^{+}, s^{+}) ← (x, s) + \overline{θ}(u, v);

end for

The algorithm begins each iteration by trying a fast step, which uses an affine scaling search direction. The fast steps are accepted only if they produce a reduction in τ_{k} or ∥r^{k}∥ of at least a factor of ρ. Otherwise, the algorithm reverts to taking a safe step. Then it goes to improve by reusing the coefficient matrix in (2.5) and taking a combination of safe and fast steps, just like the main algorithm. However, the improve procedure terminates if τ and ∥r∥ are not improved by at least a factor of \overline{ρ} ∈ (ρ, 1). The parameter \overline{ρ} and the nonnegative integer I are supplied by the user, where I is the maximum number of steps that can be taken in improve.

In the next lemma we show that if the main points defined by Algorithm 2.1 are generated in the safe branch then the improvement rate τ_{+}/τ is bounded by a quantity that increases with δ = ∥Uv∥/τ. Later on we will prove global convergence by showing that δ is bounded.

Lemma 2.2 Assume that (x^{k}, s^{k}, τ_{k}, α_{k}, β_{k}) ∈ N and suppose that the main points in Algorithm 2.1 are generated by the safe branch. Then \overline{θ} ∈ [0, 1] defined by (2.6) satisfies

\overline{θ} ≥ θ̂:= min\{1, \frac{λ min((1 − α)ρ_{0} , β)}{δ}\},

where δ = ∥Uv∥/τ. Moreover,

τ^{+}/τ ≤ 1 − min\{1, \frac{λ min((1 − α)ρ_{0} , β)}{δ}\}(1 − λ).  (2.7)

Proof. Since (x, s, τ, α, β) ∈ N, we have

Xs ≥ ατp, ∥Xs − τp∥ ≤ βτ.

By definition we get

X(θ)s(θ) = Xs + θ(Su + Xv) + θ^{2}Uv
= (1 − θ)Xs + θλτp + θ^{2}Uv.

Hence, we obtain

X(θ)s(θ) − α(1 − (1 − λ)θ)τp
= (1 − θ)(Xs − ατp) + θ[(1 − α)λτp + θUv]
≥ θ[(1 − α)λτp + θUv]
≥ θτ[(1 − α)λp − θδe]
≥ 0, for all θ ∈ [0, θ̂].

On the other hand, we have

∥X(θ)s(θ) − (1 − (1 − λ)θ)τp∥ − β(1 − (1 − λ)θ)τ
= ∥(1 − θ)(Xs − τp) + θ^{2}Uv∥ − β(1 − (1 − λ)θ)τ
≤ (1 − θ)∥Xs − τp∥ + θ^{2}∥Uv∥ − β(1 − (1 − λ)θ)τ
≤ θτ(θδ − λβ)
≤ 0, for all θ ∈ [0, θ̂].

Consequently, \overline{θ} ≥ θ̂, and (2.7) follows immediately. □

Theorem 2.3 For any integer k ≥ 0, Algorithm 2.1 defines (x^{k}, s^{k}, τ_{k}, α_{k}, β_{k}) such that

(x^{k}, s^{k}, τ_{k}, α_{k}, β_{k}) ∈ N ,  (2.8)
r^{k} = ψ_{k}r^{0},  (2.9)

where

ψ_{k} = τ_{k}/τ_{0} ≤ ∏_{j=1}^{k}(1 − (1 − λ_{j})\overline{θ}_{j}).  (2.10)

Moreover all the intermediary points produced in the improve procedure also satisfy

(x, s, τ, α, β) ∈ N .  (2.11)

Proof. The proof is by induction. First, we note that

∥X^{0}s^{0} − τ_{0}p∥ = ∥X^{0}s^{0} − μ_{0}p∥ = 0 < β_{0}τ_{0},

and

X^{0}s^{0} ≥ α_{max}X^{0}s^{0} = α_{0}μ_{0}p = α_{0}τ_{0}p,

which shows that (x^{0}, s^{0}, τ_{0}, α_{0}, β_{0}) ∈ N. Hence (2.8) and (2.9) are satisfied for k = 0. Suppose they are satisfied for k = 0, 1, …, l for some l ≥ 0 and denote (x, s) = (x^{l}, s^{l}). Then by Lemma 2.2, \overline{θ} exists, and we have

X(θ)s(θ) ≥ α(1 − (1 − λ)θ)τp > 0, for all 0 ≤ θ ≤ \overline{θ}.  (2.12)

The positivity of x^{+} and s^{+} is proved by contradiction. Suppose for example that [x^{+}]_{i} ≤ 0 for some i. Then there must exist θ′, 0 ≤ θ′ ≤ \overline{θ}, such that [x(θ′)]_{i} = 0. This implies

[x(θ′)]_{i}[s(θ′)]_{i} = 0,

which contradicts (2.12). Hence (2.8) is satisfied for k = l + 1. The fact that (2.9) is verified for k = l + 1 follows from the definition of Algorithm 2.1. Similarly we can prove (2.11). □

In next section, we will show that with appropriate choice of λ and (x^{0}, s^{0}), Algorithm 2.1 achieves O((κ + 1)^{2}nL)-iteration complexity.

3 Global Convergence and Polynomial Complexity

In this section we assume that F^{∗} is nonempty. Under this assumption we will prove that Algorithm 2.1, with ϵ = 0, is globally convergent in the sense that

lim_{k→0} μ_{k} = 0, lim_{k→0} r^{k} = 0,

for any starting point (x^{0}, s^{0}) > 0. We start with the following three technical lemmas.

Lemma 3.1 Assume that F^{∗} is nonempty. Then for any (x^{∗}, s^{∗}) ∈ F^{∗} and (x, s, τ, α, β) ∈ N we have

ψ((x)^{T} s^{0} + (s)^{T} x^{0}) ≤ (1 + β_{max})(1 + 4κ)(2 + ζ)nτ ,  (3.1a)
(1 − ψ)((x)^{T} s^{∗} + (s)^{T} x^{∗}) ≤ (1 + β_{max})(1 + 4κ)((1 + ψ) + (1 − ψ)ζ)nτ , (3.1b)

where

ψ = τ/τ_{0}, ζ = ((x^{0})^{T} s^{∗} + (s^{0})^{T} x^{∗})/((x^{0})^{T} s^{0}).  (3.2)

Proof. According to (2.9) we have

ψs^{0} + (1 − ψ)s^{∗} − s = ψ(r^{0} + M(x^{0} − x^{∗})) − (r + M(x − x^{∗}))
= M(ψx^{0} + (1 − ψ)x^{∗} − x) .

Using the inequalities (x^{∗}, s^{∗}) ≥ 0, (x, s) > 0, and the defining property (1.3) of a P_{∗}(κ)- matrix we can write

[ψx^{0} + (1 − ψ)x^{∗} − x]^{T} [ψs^{0} + (1 − ψ)s^{∗} − s]  (3.3)

≥ −4κ ∑_{i∈I_{+}}[ψx^{0} + (1 − ψ)x^{∗} − x]_{i}[ψs^{0} + (1 − ψ)s^{∗} − s]_{i}

≥ −4κ ∑_{i∈I_{+}} (ψ^{2}[x^{0}]_{i}[s^{0}]_{i} + (1 − ψ)ψ([x^{∗}]_{i}[s^{0}]_{i} + [x^{0}]_{i}[s^{∗}]_{i}) + [x]_{i}[s]_{i})

≥ −4κ(ψ^{2}(x^{0})^{T} s^{0} + (1 − ψ)ψ((x^{∗})^{T} s^{0} + (x^{0})^{T} s^{∗}) + x^{T} s) ,  (3.4)

where

I_{+} = {i : [ψx^{0} + (1 − ψ)x^{∗} − x]_{i}[ψs^{0} + (1 − ψ)s^{∗} − s]_{i} > 0} .

By expanding (3.3), we obtain

[ψx^{0} + (1 − ψ)x^{∗} − x]^{T} [ψs^{0} + (1 − ψ)s^{∗} − s]
= ψ^{2}nμ_{0} + (1 − ψ)ψ((x^{0})^{T} s^{∗} + (s^{0})^{T} x^{∗})
−ψ((x^{0})^{T} s + (s^{0})^{T} x) + x^{T} s − (1 − ψ)(s^{T} x^{∗} + x^{T} s^{∗}) + (1 − ψ)(x^{∗})^{T} s^{∗}. (3.5)

The desired inequalities (3.1) follows from (3.4) and (3.5) by applying Theorem 2.3 and the relations (x^{∗})^{T} s^{∗} = 0, s^{T} x^{∗}+x^{T} s^{∗} ≥ 0, s^{T} x^{0}+x^{T} s^{0} > 0, and μ ≤ (1+β_{max})τ (see (2.4)) . □

In order to prove global convergence we have to show that \overline{θ}’s are bounded away from zero. An essential tool is the following technical result of [4].

Lemma 3.2 Let x, s, a, r be four n-dimensional vectors with x > 0 and s > 0, and let M ∈ R I ^{n×n} be a P_{∗}(κ)-matrix. Then the solution (u, v) of the linear system

Su + Xv = a ,
Mu − v = b

satisfies the following relations:

∥Du∥ ≤ ∥b̃∥ + \sqrt{∥ã∥^{2} + ∥b̃∥^{2} + 2κ∥c̃∥^{2}} ,
∥D^{−1}v∥ ≤ \sqrt{∥ã∥^{2} + ∥b̃∥^{2} + 2κ∥c̃∥^{2}} ,
∥Du∥^{2} + ∥D^{−1}v∥^{2} ≤ ∥ã∥^{2} + 2κ∥c̃∥^{2} + 2∥b̃∥^{2} + 2∥b̃∥\sqrt{∥ã∥^{2} + ∥b̃∥^{2} + 2κ∥c̃∥^{2}}
≡ χ_{1}^{2} ,
∥Uv∥^{2} ≤ \frac{1}{8}∥ã∥^{4} + \frac{1}{4}χ_{1}^{2}(χ_{1}^{2} − ∥ã∥^{2}) ,

where

D = X^{−1/2}S^{1/2} , ã = (XS)^{−1/2}a , b̃ = D^{−1}b , c̃ = ã + b̃ .

We use the above result to find bounds for the quantity δ arising in the estimate given by Lemma 2.2 when the algorithm chooses the safe step.

Lemma 3.3 During the safe branch, δ = ∥Uv∥/τ satisfies

δ ≤ ϕ(κ) \{(1 − λ)[∥p∥/(αρ_{0}) + ∥X^{1/2}S^{−1/2}r∥/\sqrt{τ}] + β/(αρ_{0})\}^{2} ,  (3.6a)

where

ϕ(κ) := [.25(\sqrt{1 + 2κ} + 1)^{4} + .125]^{1/2}.  (3.6b)

Moreover,

∥X^{1/2}S^{−1/2}r∥/\sqrt{τ} ≤ η\sqrt{n(1 + β)},  (3.7a)

where

η = \sqrt{n}(1 + 4κ)(2 + ζ)∥(S^{0})^{−1}r^{0}∥_{∞}\sqrt{\frac{1 + β_{max}}{ρ_{0}α_{min}}} ,  (3.7b)

with ζ given by (3.2). Proof. By applying Lemma 3.2 to the linear system (2.5) with a = λτp−Xs and b = (1−λ)r we have

∥ã∥ = ∥(XS)^{−1/2}(λτp − Xs)∥
≤ (1 − λ)τ∥(XS)^{−1/2}p∥ + ∥(XS)^{−1/2}(Xs − τp)∥
≤ (1 − λ)\sqrt{τ}∥p∥/(αρ_{0}) + \sqrt{τ}β/(αρ_{0}),
∥b̃∥ = ∥D^{−1}r∥ = (1 − λ)∥X^{1/2}S^{−1/2}r∥,
∥c̃∥ = ∥ã + b̃∥ ≤ ∥ã∥ + ∥b̃∥,
δ = ∥Uv∥/τ ≤ ϕ(κ)(∥ã∥ + ∥b̃∥)^{2}/τ,

and (3.6a) follows from the above relations. From Lemma 3.1 we deduce that

∥(XS)^{−1/2}Xr∥/\sqrt{τ} ≤ (ρ_{0}α)^{−1/2}τ^{−1}∥Xr∥

≤ (ρ_{0}α)^{−1/2}τ^{−1}∥Xr∥_{1} = ψ(ρ_{0}α)^{−1/2}τ^{−1}∥Xr^{0}∥_{1}
= ψ(ρ_{0}α)^{−1/2}τ^{−1} ∑_{i=1}^{n} |[x]_{i}[s^{0}]_{i}[r^{0}]_{i}/[s^{0}]_{i}|
≤ (ρ_{0}α)^{−1/2}τ^{−1}∥(S^{0})^{−1}r^{0}∥_{∞}ψ(s^{0})^{T} x
≤ η\sqrt{n(1 + β)}  □

Now let us take

λ = \overline{λ} := \frac{∥p∥/(αρ_{0}) + ∥X^{1/2}S^{−1/2}r∥/\sqrt{τ} + β/(αρ_{0})}{∥p∥/(αρ_{0}) + ∥X^{1/2}S^{−1/2}r∥/\sqrt{τ} + ω β/(αρ_{0})},  (3.8)

where ω > 1 is a fixed number. Clearly,

\overline{λ} ∈ (1/ω, 1).  (3.9)

From (3.6a) we have

δ ≤ ω^{2}ϕ(κ)β^{2}/(αρ_{0})^{2}.  (3.10)

Then by virtue of (2.7), (3.8), (3.9) and (3.10), we deduce that

τ^{+}/τ ≤ 1 − min \{1, λ\frac{min((1 − α)ρ_{0} , β)}{δ}\} (1 − λ)
≤ 1 − min \{1, \frac{min((1 − α)ρ_{0} , β)}{ω^{3}ϕ(κ)β^{2}/(αρ_{0})^{2}}\} (1 − \overline{λ})
≤ 1 − min \{1, \frac{min\{(1 − α)ρ_{0} , β\}}{ω^{3}ϕ(κ)β^{2}/(αρ_{0})^{2}}\} \frac{(ω − 1)β}{∥p∥ + ηαρ_{0}\sqrt{n(1 + β)} + ωβ}
≤ 1 − θ^{∗},  (3.11)

where

θ^{∗} := min \{1, \frac{min\{(1 − α_{max})ρ_{0} , β_{min}\}}{ω^{3}ϕ(κ)β_{max}^{2}/(α_{min}ρ_{0})^{2}}\} \frac{(ω − 1)β_{min}}{∥p∥ + ηα_{max}ρ_{0}\sqrt{n(1 + β_{max})} + ωβ_{max}} (3.12)

Suppose the algorithm takes the fast branch at iteration k. Since the improve procedure never increases the value of ψ, we must have ψ_{k+1}/ψ_{k} ≤ ρ; if a safe branch is taken we must have ψ_{k+1}/ψ_{k} ≤ 1 − θ^{∗}. Therefore we obtain

ψ_{k+1}/ψ_{k} ≤ max{1 − θ^{∗}, ρ}.  (3.13)

This observation shows that our algorithm converges at a global linear rate.

Theorem 3.4 Suppose that the optimal set F^{∗} is nonempty, and in the safe branch λ_{k} is defined by (3.8).

(i) If ϵ = 0 then Algorithm 2.1 either finds an optimal solution z^{∗} ∈ F^{∗} in a finite number of steps or produces an infinite sequence z^{k} = (x^{k}, s^{k}) such that

lim_{k→∞}(x^{k})^{T} s^{k} = 0, lim_{k→∞}(r^{k}) = 0.

(ii) If ϵ > 0 then Algorithm 2.1 terminates with a z ∈ F_{ϵ} in at most

K_{ϵ} = ⌈\frac{| ln(\frac{ϵ}{ϵ_{0}} )|}{| ln(max\{1 − θ^{∗}, ρ\})|}⌉

iterations, where ϵ_{0} = max{(1 + β)(x^{0})^{T} s^{0}, ∥r^{0}∥}, and ⌈χ⌉ denotes the smallest integer greater than or equal to χ.

For starting points that are feasible or close to being feasible we obtain the following complexity result.

Corollary 3.5 Under the hypothesis of Theorem 3.4 suppose that ρ_{0} = Ω(1), and the starting point is chosen such that there is a constant C independent of n satisfying the inequality

(2 + ζ)∥(S^{0})^{−1}r^{0}∥_{∞} ≤ \frac{C}{(1 + κ)\sqrt{n}}.

Then Algorithm 2.1 terminates in at most

K̃_{ϵ} = O((1 + κ)\sqrt{n} ln(ϵ_{0}/ϵ))

iterations.

Most of the complexity results on infeasible interior-point methods are obtained for starting points of the form

x^{0} = ρ_{p}e, s^{0} = ρ_{d}e ,  (3.14)

where ρ_{p} and ρ_{d} are sufficiently large positive constants (“big M initialization”). For such starting points we have clearly ρ_{0} = 1 and

ζ = ∥x^{∗}∥_{1}/(nρ_{p}) + ∥s^{∗}∥_{1}/(nρ_{d}), for some (x^{∗}, s^{∗}) ∈ F^{∗} ,
∥(S^{0})^{−1}r^{0}∥_{∞} ≤ 1 + (ρ_{p}/ρ_{d})∥Me∥_{∞} + (1/ρ_{d})∥q∥_{∞}.

Therefore if ρ_{p} and ρ_{d} satisfy the inequalities

ρ_{p} ≥ n^{−1}∥x^{∗}∥_{1} , ρ_{d} ≥ max{ρ_{p}∥Me∥_{∞} , ∥q∥_{∞}, n^{−1}∥s^{∗}∥_{1}} ,  (3.15)

for some (x^{∗}, s^{∗}) ∈ F^{∗} then η ≤ O((1 + κ)\sqrt{n}). Hence, from (3.12) and Theorem 3.4 we deduce the following computational complexity bound.

Corollary 3.6 Under the hypothesis of Theorem 3.4 suppose that the starting point satisfies (3.15) for some (x^{∗}, s^{∗}) ∈ F^{∗}. Then Algorithm 2.1 terminates in at most

K̃_{ϵ} = O((1 + κ)^{2}n ln(ϵ_{0}/ϵ))  (3.16)

iterations.

Let us end this section by noting that while our algorithm does not use any information on the classification number κ, its computational complexity depends on κ.

4 Infeasibility Detection

All the above results have been proved under the assumption that F^{∗} is nonempty. It turns out that Algorithm 2.1 can be modified in such a way that it can detect within polynomial time whether F^{∗} contains points of norm less than a quantity chosen in advance provided the parameter κ is known. Let \overline{ρ}_{p}, \overline{ρ}_{d} be such quantities. Also let us define

\overline{ζ} = (∥x^{0}∥\overline{ρ}_{d} + ∥s^{0}∥\overline{ρ}_{p})/((x^{0})^{T} s^{0}),

\overline{η} = \sqrt{n}(1 + 4κ)(2 + \overline{ζ})∥(S^{0})^{−1}r^{0}∥_{∞}\sqrt{\frac{1 + β_{max}}{α_{min}ρ_{0}}} ,  (4.1)

\overline{θ}^{∗} := min \{1, \frac{min\{(1 − α_{max})ρ_{0} , β_{min}\}}{ω^{3}ϕ(κ)β_{max}^{2}/(α_{min}ρ_{0})^{2}}\} \frac{(ω − 1)β_{min}}{∥p∥ + \overline{η}α_{max}ρ_{0}\sqrt{n(1 + β_{max})} + ωβ_{max}} .

In what follows, we always assume that λ is defined by (3.8). Then we can prove the following theorem.

Theorem 4.1 Suppose that the following instruction “If (x^{0})^{T} s^{k}+(s^{0})^{T} x^{k} > (1+4κ)[(τ_{k}/τ_{0})(x^{0})^{T} s^{0}+(τ_{0}/τ_{k})(x^{k})^{T} s^{k}+(1−τ_{k}/τ_{0})(\overline{ρ}_{d}∥x^{0}∥+\overline{ρ}_{p}∥s^{0}∥)] then terminate” is inserted just before the while loop of Algorithm 2.1. Then the new algorithm terminates at most

\overline{K}_{ϵ} = ⌈\frac{| ln(\frac{ϵ}{ϵ_{0}} )|}{| ln(max\{1 − \overline{θ^{∗}}, ρ\})|}⌉

iterations either by finding a ϵ − approximate solution or by determining either there is no z^{∗} = (x^{∗}, s^{∗}) ∈ F^{∗} with ∥x^{∗}∥ ≤ \overline{ρ}_{p}, ∥s^{∗}∥ ≤ \overline{ρ}_{d}, or M is not a P_{∗}(κ)-matrix.

Proof. Suppose that M is a P_{∗}(κ)-matrix and that the inequality

(x^{0})^{T} s^{k} + (s^{0})^{T} x^{k} ≤ (1 + 4κ)[(τ_{k}/τ_{0})(x^{0})^{T} s^{0} + (τ^{0}/τ_{k})(x^{k})^{T} s^{k} + (1 − τ_{k}/τ_{0})(\overline{ρ}_{d}∥x^{0}∥ + \overline{ρ}_{p}∥s^{0}∥)]

(4.2)

holds for all 0 ≤ k ≤ \overline{K}_{ϵ}. Then we have

ψ_{k}((x^{0})^{T} s^{k} + (s^{0})^{T} x^{k}) ≤ (1 + 4κ)(1 + β_{max})(2 + \overline{ζ})nτ_{k}, 0 ≤ k ≤ \overline{K}_{ϵ}.  (4.3)

With the help of (4.3) we can show that τ_{k+1}/τ_{k} ≤ max\{1 − \overline{θ^{∗}}, ρ\} for all 0 ≤ k ≤ \overline{K}_{ϵ}, which implies (x^{\overline{K}^{ϵ}}, s^{\overline{K}^{ϵ}}) ∈ F_{ϵ}. On the other hand, if there exists (x^{∗}, s^{∗}) ∈ F^{∗} such that ∥x^{∗}∥ ≤ \overline{ρ}_{p}, ∥s^{∗}∥ ≤ \overline{ρ}_{d}, then (3.4) and (3.5) imply that ( 4.2) must hold for all k ≥ 0. Therefore, if (4.2) is violated for some k ≤ \overline{K}_{ϵ}, then there is no (x^{∗}, s^{∗}) ∈ F^{∗} such that ∥x^{∗}∥ ≤ \overline{ρ}_{p}, ∥s^{∗}∥ ≤ \overline{ρ}_{d}. This completes the proof of our theorem. □ From the above theorem we can obtain polynomial complexity for our new algorithm under certain assumptions on the starting point.

Corollary 4.2 Under the hypothesis of Theorem 4.1 suppose ρ_{0} = Ω(1) and the starting point satisfies

(2 + \overline{ζ})∥(S^{0})^{−1}r^{0}∥_{∞} ≤ \frac{C}{(1 + κ)\sqrt{n}},

for some constant C independent of n. Then the new algorithm terminates in at most

K̃_{ϵ} = O((1 + κ)\sqrt{n} ln(ϵ_{0}/ϵ))

iterations either by finding an ϵ − approximate solution or by determining either there is no z^{∗} = (x^{∗}, s^{∗}) ∈ F^{∗} with ∥x^{∗}∥ ≤ \overline{ρ}_{p}, ∥s^{∗}∥ ≤ \overline{ρ}_{d}, or M is not a P_{∗}(κ)-matrix.

Suppose now that the starting point satisfies

x^{0} = ρ̂_{p}e, s^{0} = ρ̂_{d}e,  (4.4)

where

ρ̂_{p} ≥ \overline{ρ}_{p}/\sqrt{n},  (4.5)
ρ̂_{d} ≥ max\{ρ̂_{p}∥Me∥_{∞}, ∥q∥_{∞}, \overline{ρ}_{d}/\sqrt{n}\}.  (4.6)

Then we obtain

\overline{ζ} ≤ 2,  (4.7)
∥(S^{0})^{−1}r^{0}∥_{∞} ≤ 1 + (ρ̂_{p}/ρ̂_{d})∥Me∥_{∞} + (1/ρ̂_{d})∥q∥_{∞} ≤ 3.  (4.8)

From (4.1), (4.7) and (4.8), it follows that \overline{η} < O((1 + κ)\sqrt{n}). Hence we have the following complexity result.

Corollary 4.3 Under the assumption of Theorem 4.1 suppose that the starting point satisfies (4.4)-(4.6). Then the new algorithm terminates in at most

K̂_{ϵ} = O((1 + κ)^{2}n ln(ϵ_{0}/ϵ))

iterations either by finding an ϵ − approximate solution or by determining either there is no z^{∗} = (x^{∗}, s^{∗}) ∈ F^{∗} with ∥x^{∗}∥ ≤ \overline{ρ}_{p}, ∥s^{∗}∥ ≤ \overline{ρ}_{d}, or M is not a P_{∗}(κ)-matrix.

5 Technical Results

In what follows we study the asymptotic convergence properties of Algorithm 2.1 under the further assumption that (1.1) has a strictly complementary solution. We prove some results about the steps generated in the fast branch of the algorithm. In particular, we analyze the effect of using an approximate Jacobian for solving (2.5) as is done during the improve procedure. Let us denote by F^{c} the set of all strictly complementary solutions, i.e.,

F^{c} = {(x, s) ∈ F^{∗} : [x]_{i} + [s]_{i} > 0, i = 1, 2, ⋯, n} .

It is well known that there is a unique partition

B ∪ N = {1, 2, ⋯, n} , B ∩ N = ∅ ,

such that for any (x, s) ∈ F^{c} we have ([x]_{i} > 0, [s]_{i} = 0, ∀i ∈ B) and ([x]_{i} = 0, [s]_{i} > 0, ∀i ∈ N). Let us denote the corresponding partition of M by

M = ( \begin{matrix}M_{BB} & M_{BN} \\ M_{NB} & M_{NN} \end{matrix} ) .

Also, for any vector y ∈ R I ^{n} we denote by y_{B} the vector of components [y]_{i}, i ∈ B, and by y_{N} the vector of components [y]_{i}, i ∈ N.

From the next lemma it follows that the points generated by Algorithm 2.1 are bounded.

Lemma 5.1 There is a constant C_{1} > 0 such that

∥(x, s)∥ ≤ C_{1},  (5.1)

for all (x, s, τ, α, β) ∈ N .

Proof. The proof is straitforward by observing that (3.1a) implies

x^{T} s^{0} + s^{T} x^{0} ≤ (1 + β)(1 + 4κ)(2 + ζ)nτ_{0}. □

Lemma 5.2 There is a constant C_{2} > 0 such that for any (x, s, τ, α, β) ∈ N and for any z^{∗} := (x^{∗}, s^{∗}) ∈ F^{∗} we have

∥(x + u^{a} − x^{∗}, s + v^{a} − s^{∗})∥ ≤ C_{2} \frac{∥z^{k} − z^{∗}∥^{2}}{τ} ,  (5.2)

where z := (x, s) and (u^{a}, v^{a}) satisfy

Su^{a} + Xv^{a} = −Xs,  (5.3a)
Mu^{a} − v^{a} = r.  (5.3b)

Proof. Let z^{∗} = (x^{∗}, s^{∗}) ∈ F^{∗}. By using (5.3) we have:

( \begin{matrix}S & X \\ M & −I \end{matrix} ) ( \begin{matrix}x + u^{a} − x^{∗} \\ s + v^{a} − s^{∗} \end{matrix} ) = ( \begin{matrix}(X − X^{∗})(s − s^{∗}) \\ 0 \end{matrix} ) .

Applying Lemma 3.2 to the above linear system, we obtain

∥D(x + u^{a} − x^{∗})∥ ≤ \sqrt{1 + 2κ}∥(XS)^{− \frac{1}{2}} (X − X^{∗})(s − s^{∗})∥
≤ \sqrt{1 + 2κ}\begin{matrix}∥x − x^{∗}∥ & ∥s − s^{∗}∥ \\ √ & \overline{α_{min}τ} \end{matrix}
≤ \frac{\sqrt{1 + 2κ}}{2 \sqrt{α_{min}}} \begin{matrix}∥z − z^{∗}∥ & 2 \\ √ & \overline{τ} \end{matrix} ,  (5.4)

where D = X^{− \frac{1}{2}} S^{\frac{1}{2}} . Thus

∥x + u^{a} − x^{∗}∥ = ∥D^{−1}D (x + u^{a} − x^{∗})∥
≤ ∥X∥∥(XS)^{−1/2}∥ ∥D(x + u^{a} − x^{∗})∥
≤ \frac{\sqrt{1 + 2κ} C_{1}}{2α_{min}} \frac{∥z − z^{∗}∥^{2}}{τ}
= C_{2} \frac{∥z − z^{∗}∥^{2}}{τ} , with C_{2} := \frac{\sqrt{1 + 2κ}C_{1}}{2α_{min}} .

The corresponding inequality for s + v^{a} − s^{∗} can be obtained similarly. □

Lemma 5.3 If F^{c} is nonempty, then there is a constant C_{3} > 0 such that for any (x, s, τ, α, β) ∈ N with τ ≤ τ_{1},

∥x_{N}∥ ≤ C_{3} τ, ∥s_{B}∥ ≤ C_{3}τ,

for all k ≥ 1.

Proof. Let (x^{∗}, s^{∗}) ∈ F^{c}. If ψ = τ/τ_{0}, then ψ ≤ ψ_{1} ≤ 1 − θ_{0} so that we can write

(x)^{T} s^{∗} + (s)^{T} x^{∗} ≤ (1 + β_{max})(1 + 4κ)((1 + ψ_{1})/(1 − ψ_{1}) + ζ)nτ_{k}
≤ (1 + β_{max})(1 + 4κ)((2 − \overline{θ}_{0})/\overline{θ}_{0} + ζ)nτ_{k} .

Therefore,

[x]_{i} ≤ \frac{(1 + β_{max})(1 + 4κ)((2 − \overline{θ}_{0})/\overline{θ}_{0} + ζ)}{[s^{∗}]_{i}}nτ_{k} , ∀ i ∈ N ,

and

[s]_{i} ≤ \frac{(1 + β_{max})(1 + 4κ)((2 − \overline{θ}_{0})/\overline{θ}_{0} + ζ)}{[x^{∗}]_{i}}nτ_{k} , ∀ i ∈ B .

Hence the desired result holds with

C_{3} = n\sqrt{n}(1 + β_{max})(1 + 4κ)((2 − \overline{θ}_{0})/\overline{θ}_{0} + ζ) max \{max_{i∈B} \frac{1}{[x^{∗}]_{i}} , max_{i∈N} \frac{1}{[s^{∗}]_{i}} \} . □

Lemma 5.4 If F^{c} is nonempty, then there exists a constant C_{4} > 0 such that for each (x, s, τ, α, β) ∈ N there is a solution z^{∗} = (x^{∗}, s^{∗}) ∈ F^{∗} such that

∥z − z^{∗}∥ ≤ C_{4}τ.  (5.5)

Proof. Consider the following problem:

\begin{aligned}M_{BB}x_{B} & & & & = −q_{B} \\ M_{NB}x_{B} & & & −s_{N} & = −q_{N} \\ & x_{N} & & & = 0 \\ & & s_{B} & & = 0 \\ & & x_{B}, & s_{N} & ≥ 0 . \end{aligned}  (5.6)

Under the assumption that F^{c} ≠ ∅, (5.6) has a solution and the solution set of (5.6) is F^{∗}. By Hoffman’s lemma [3] , there is a constant C′_{4} such that there is a z^{∗} ∈ F^{∗} satisfying

∥z − z^{∗}∥ ≤ C′_{4} ∥(M_{BB}x_{B} + q_{B}, M_{NB}x_{B} + q_{N} − s_{N}, x_{N}, s_{B})∥
≤ C′_{4} ∥(−M_{BN}x_{N} + s_{B}, −M_{NN}x_{N}^{k}, x_{N}, s_{B})∥ + C′_{4}∥r∥ .  (5.7)

We have

∥r∥ = ∥ψr^{0}∥ = \frac{∥r^{0}∥}{τ_{0}}(ψτ_{0}) = \frac{∥r^{0}∥}{τ_{0}}τ.  (5.8)

Finally, (5.5) follows from Lemma 5.3 and (5.7)-(5.8) . □

Lemma 5.5 There exists a constant C_{5} such that for any (x, s, τ, α, β) ∈ N we have

∥(u^{a}, v^{a})∥ ≤ C_{5}τ,  (5.9)

where (u^{a}, v^{a}) is the solution of (5.3).

Proof. Let z^{∗} = (x^{∗}, s^{∗}) ∈ F^{∗} satisfy (5.5). From Lemma 5.2 and Lemma 5.4 we deduce that

∥u^{a}∥ ≤ ∥x + u^{a} − x^{∗}∥ + ∥x − x^{∗}∥
≤ (C_{2}C_{4}^{2} + C_{4})τ := C_{5}τ .

The corresponding inequality for ∥v^{a}∥ is obtained in a similar manner. □

We now turn to the approximate fast steps computed by (2.5)-(2.6), where (x, s) is either the current iterate (x^{k}, s^{k}) or some intermediate point generated in the call to improve at the k-th iteration. By the definition of the algorithm, we have

τ ≤ τ_{k}.  (5.10)

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

∥(x^{k} − x, s^{k} − s)∥ ≤ χτ_{k}.  (5.11)

We will show later on that such a χ exists (see Lemma 6.1). The following lemma describes some properties of the actual search direction (u, v) calculated from (2.6), in terms of the exact search direction (u^{a}, v^{a}) that satisfies (5.3).

Lemma 5.6 Suppose that F^{c} is nonempty and consider the notation of Algorithm 2.1. Under the assumption (5.11) suppose that (u, v) is the solution of (2.5) with λ = 0, then there exists a constant C_{6} independent of k and χ such that the following bounds are satisfied:

∥u − u^{a}∥ ≤ C_{6}χτ, ∥v − v^{a}∥ ≤ C_{6}χτ,  (5.12a)
∥(u, v)∥ ≤ C_{6}χτ,  (5.12b)
∥u_{N} − u_{N}^{a}∥ ≤ C_{6}χττ_{k}, ∥v_{B} − v_{B}^{a}∥ ≤ C_{6}χττ_{k}.  (5.12c)

Proof. From (2.5), we have that

( \begin{matrix}S^{k} & X^{k} \\ M & −I \end{matrix} ) ( \begin{matrix}u \\ v \end{matrix} ) = ( \begin{matrix}−Xs \\ r \end{matrix} ) ,  (5.13)

while from (5.3), we get

( \begin{matrix}S & X \\ M & −I \end{matrix} ) ( \begin{matrix}u^{a} \\ v^{a} \end{matrix} ) = ( \begin{matrix}−Xs \\ r \end{matrix} ) ,  (5.14)

and therefore

( \begin{matrix}S^{k} & X^{k} \\ M & −I \end{matrix} ) ( \begin{matrix}u^{a} \\ v^{a} \end{matrix} ) = ( \begin{matrix}−Xs + (S^{k} − S)u^{a} + (X^{k} − X)v^{a} \\ r \end{matrix} ) .  (5.15)

By (5.13) and (5.15) we obtain

( \begin{matrix}S^{k} & X^{k} \\ M & −I \end{matrix} ) ( \begin{matrix}u^{a} − u \\ v^{a} − v \end{matrix} ) = ( \begin{matrix}(S^{k} − S)u^{a} + (X^{k} − X)v^{a} \\ 0 \end{matrix} ) .  (5.16)

Then from (5.11) and Lemma 5.5, there is a constant C′_{6} independent of k and χ such that

∥(S^{k} − S)u^{a} + (X^{k} − X)v^{a}∥ ≤ C′_{6}χττ_{k}.  (5.17)

Applying Lemma 3.2 to the linear system (5.16) with a = (S^{k} − S)u^{a} + (X^{k} − X)v^{a}, b = 0, we deduce that

∥D^{k}(u^{a} − u)∥ ≤ \sqrt{1 + 2κ}∥(X^{k}S^{k})^{−1/2}((S^{k} − S)u^{a} + (X^{k} − X)v^{a})∥
≤ \frac{\sqrt{1 + 2κ}C′_{6}}{\sqrt{α_{min}ρ_{0}}} χτ\sqrt{τ_{k}},  (5.18)

where D^{k} = (X^{k})^{−1/2}(S^{k})^{1/2}. Therefore

∥u^{a} − u∥ ≤ ∥(D^{k})^{−1}∥∥D^{k}(u^{a} − u)∥ ≤ \frac{\sqrt{1 + 2κ}C_{1}C′_{6}}{α_{min}ρ_{0}} χτ.  (5.19)

Then from (5.18) and Lemma 5.3 we deduce that

∥u_{N}^{a} − u_{N}∥ ≤ ∥(D_{N}^{k})^{−1}∥∥D^{k}(u^{a} − u)∥
= ∥(X_{N}^{k})(X_{N}^{k}S_{N}^{k})^{−1/2}∥∥D^{k}(u^{a} − u)∥
≤ \frac{\sqrt{1 + 2κ}C_{3}C′_{6}}{α_{min}ρ_{0}} χττ_{k}.  (5.20)

The bounds for (v^{a} − v) and (v_{B}^{a} − v_{B}) are obtained similarly. Since

∥(u, v)∥ ≤ ∥(u^{a}, v^{a})∥ + ∥(u − u^{a}, v − v^{a})∥,

relation (5.12) follows from Lemma 5.5 with an appropriate C_{6} > 0. □

The next lemma gives an estimate for the step length \overline{θ} along a (possibly approximate) fast step direction (u, v). The point (x, s) considered in the lemma represents either the main iterate (x^{k}, y^{k}) itself or one of the intermediate points generated by the improve procedure during the k-th iteration. For the purpose of this result, let us define

C′_{7} = ρ_{0}(1 − γ) min{α_{max} − α_{min}, β_{max} − β_{min}},
C′′_{7} = max{2C_{5}^{2}, 4C_{6}(C_{1} + C_{3})},
C_{7} = C′′_{7}/C′_{7}.

Lemma 5.7 During the k-th iteration of Algorithm 2.1 suppose (x, s, τ, α, β) ∈ N, where

α = α_{min} + γ^{t}(α_{max} − α_{min}), β = β_{max} − γ^{t}(β_{max} − β_{min}),

and the quantity t satisfies

C_{7}χ^{2} \frac{τ_{k}}{γ^{t}} ≤ ρ.  (5.21)

If \overline{θ} is computed from (2.5)-(2.6) with λ = 0 and

α_{+} = α_{min} + γ^{t+1}(α_{max} − α_{min}), β_{+} = β_{max} − γ^{t+1}(β_{max} − β_{min}),

then

1 − \overline{θ} ≤ C_{7}χ^{2} \frac{τ_{k}}{γ^{t}} ≤ ρ.  (5.22)

Proof. By definition, we have

X(θ)s(θ) = Xs + θ(Su + Xv) + θ^{2}Uv
= Xs + θ(Su^{a} + Xv^{a}) + θ[S(u − u^{a}) + X(v − v^{a})] + θ^{2}Uv
= (1 − θ)Xs + θ[S(u − u^{a}) + X(v − v^{a})] + θ^{2}Uv.  (5.23)

Then from Lemmas 5.1, 5.3 and 5.6, we deduce that

∥S(u − u^{a}) + X(v − v^{a})∥
≤ ∥s_{B}∥∥u_{B} − u_{B}^{a}∥ + ∥x_{B}∥∥v_{B} − v_{B}^{a}∥ + ∥s_{N}∥∥u_{N} − u_{N}^{a}∥ + ∥x_{N}∥∥v_{N} − v_{N}^{a}∥
≤ 2C_{6}(C_{1} + C_{3})χττ_{k} ≤ \frac{C′′_{7}}{2} χ^{2}ττ_{k}.  (5.24)

Let θ′′ be the unique positive root of the quadratic equation

C′_{7}γ^{t}(1 − θ) − \frac{C′′_{7}}{2} χ^{2}τ_{k}θ − \frac{C′′_{7}}{2} χ^{2}τ_{k}θ^{2} = 0.

Then it is easily seen that θ′′ ∈ (0, 1) and

C′_{7}γ^{t}(1 − θ) + \frac{C′′_{7}}{2} χ^{2}τ_{k}θ − \frac{C′′_{7}}{2} χ^{2}τ_{k}θ^{2} ≥ 0, for all θ ∈ [0, θ′′].

Hence, if θ ∈ [0, θ′′], then from (5.23), (5.24) and Lemma 5.6, we deduce that

X(θ)s(θ) − α_{+}(1 − θ)τp
= (1 − θ)(Xs − α_{+}τp) + θ[S(u − u^{a}) + X(v − v^{a})] + θ^{2}Uv
= (1 − θ)(Xs − ατp) + (1 − θ)(α − α_{+})τp + θ[S(u − u^{a}) + X(v − v^{a})] + θ^{2}Uv
≥ (1 − γ)γ^{t}(α_{max} − α_{min})(1 − θ)τp − θ∥S(u − u^{k}) + X(v − v^{k})∥e − θ^{2}∥Uv∥e
≥ C′_{7}γ^{t}(1 − θ)τe − \frac{C′′_{7}}{2} χττ_{k}θe − \frac{C′′_{7}}{2} χ^{2}τ^{2}θ^{2}e
≥ τ[C′_{7}γ^{t}(1 − θ) − \frac{C′′_{7}}{2} χ^{2}τ_{k}θ − \frac{C′′_{7}}{2} χ^{2}τ_{k}θ^{2}]e ≥ 0,

and

∥X(θ)s(θ) − (1 − θ)τp∥ − β_{+}(1 − θ)τ
≤ (1 − θ)∥Xs − τp∥ + θ∥S(u − u^{a}) + X(v − v^{a})∥ + θ^{2}∥Uv∥ − β_{+}(1 − θ)τ
≤ −(1 − γ)γ^{t}(β_{max} − β_{min})(1 − θ)τ + θ∥S(u − u^{a}) + X(v − v^{a})∥ + θ^{2}∥Uv∥
≤ −C′_{7}γ^{t}(1 − θ)τ + \frac{C′′_{7}}{2} χ^{2}ττ_{k}θ + \frac{C′′_{7}}{2} χ^{2}τ^{2}
≤ τ[−C′_{7}γ^{t}(1 − θ) + \frac{C′′_{7}}{2} χ^{2}τ_{k} + \frac{C′′_{7}}{2} χ^{2}τ_{k}] ≤ 0.

Therefore, \overline{θ} ≥ θ′′ which implies

1 − \overline{θ} ≤ 1 − θ′′ = \frac{1}{C′_{7}γ^{t}} (\frac{C′′_{7}}{2} χ^{2}τ_{k}θ′′ + \frac{C′′_{7}}{2} χ^{2}τ_{k}(θ′′)^{2})
≤ C_{7}(χ^{2}τ_{k}/γ^{t}) ≤ ρ. □

Lemma 5.8 There exists a constant C_{8} ∈ (0, 1) such that

\frac{τ_{k+1}}{γ^{t_{k+1}}} ≤ C_{8} \frac{τ_{k}}{γ^{t_{k}}} , for all k ≥ 0.  (5.25)

Proof. When the fast step is taken with a certain t, then t_{+} = t + 1 and

τ_{+} ≤ ρτ.

Hence

\frac{τ_{+}}{γ^{t_{+}}} ≤ ρ\frac{τ}{γ^{t+1}} ≤ \frac{ρ}{γ} \frac{τ}{γ^{t}} .

On the other hand, if a safe step is taken, we have from (3.13) and the definition of the algorithm that

τ_{+} ≤ max\{1 − θ^{∗}, \overline{ρ}\}τ,

while the value of t keeps unchanged. So it follows that

\frac{τ_{+}}{γ^{t_{+}}} = \frac{τ_{+}}{γ^{t}} ≤ max\{1 − θ^{∗}, \overline{ρ}\}\frac{τ}{γ^{t}} .

Therefore, (5.25) holds with C_{8} = max\{ρ/γ, 1 − θ^{∗}, \overline{ρ}\}. □

6 Local Convergence

We are ready to show that there exists a threshold value of τ_{k}/γ^{t_{k}} below which both the main algorithm and the procedure improve take only fast steps.

Lemma 6.1 Define

χ = max \{1, C_{5} exp (\frac{C_{6}ρ}{1 − ρ})\} ,  (6.1)

and let K be the smallest index such that

C_{7}χ^{2} \frac{τ_{K}}{γ^{t_{K}+I}} ≤ ρ.  (6.2)

Then at the k-th iteration with k ≥ K, the fast branch is taken in the main algorithm and, moreover, I fast steps are taken in the call to improve.

Proof. The existence of K is guaranteed by Lemma 5.8 because

\frac{τ_{k}}{γ^{t_{k}+I}} = \frac{1}{γ^{I}} \frac{τ_{k}}{γ^{t_{k}}} ≤ \frac{τ_{0}}{γ^{t_{0}+I}} C_{8}^{k} → 0, as k → ∞.

Now let us choose any k ≥ K. Our proof proceeds by showing first that the step taken from (x^{k}, s^{k}) in the main algorithm is a fast step. We then prove by induction that I fast steps are taken inside the procedure improve. Our main tool is Lemma 5.7. For the first part of the proof, we apply Lemma 5.7 with

(x, s) = (x^{k}, s^{k}), t = t_{k}.  (6.3)

Note that the point (x, s) trivially satisfies the assumption (5.11). Since

C_{7}χ^{2} \frac{τ_{k}}{γ^{t}} = C_{7}χ^{2} \frac{τ_{k}}{γ^{t_{k}}} ≤ C_{7}χ^{2} \frac{τ_{K}}{γ^{t_{K}}} ≤ C_{7}χ^{2} \frac{τ_{K}}{γ^{t_{K}+I}} ≤ ρ,

the condition (5.21) of Lemma 5.7 is satisfied by the choices (6.3), and therefore a fast step is taken by the main algorithm. We turn now to the procedure improve. Our aim is to show inductively that if (x, s) is the current vector pair at the commencement of the i-th iteration of this procedure, then

∥(x^{k}, s^{k}) − (x, s)∥ ≤ [C_{5} ∏_{l=0}^{i−1}(1 + C_{6}ρ^{l})] τ_{k},  (6.4)

and that a fast step is taken from this vector pair (x, s) during the i-th iteration of improve. Note for future reference that

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

and thus

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

Consider the case i = 1, that is the first iteration of improve. Then we have

(x, s) = (x^{k}, s^{k}) + \overline{θ}(u^{k}, v^{k}),

where (u^{k}, v^{k}) is the solution of (5.3) with (x, s) = (x^{k}, s^{k}) since a fast step was just taken at iteration k in the main algorithm. Hence by Lemma 5.5, we deduce that

∥(x^{k}, s^{k}) − (x, s)∥ ≤ ∥(u^{k}, v^{k})∥ ≤ C_{5}τ_{k}.

Thus the bound (6.4), and therefore also (5.11), holds for the point (6). Because

C_{7}χ^{2} \frac{τ_{k}}{γ^{t}} = C_{7}χ^{2} \frac{τ_{k}}{γ^{t_{k}}} ≤ C_{7}χ^{2} \frac{τ_{K}}{γ^{t_{K}}} ≤ C_{7}χ^{2} \frac{τ_{K}}{γ^{t_{K}+I}} ≤ ρ,

Lemma 5.7 applies, which implies that a fast step is taken at the first iteration of improve. We now consider the i-th iteration of the internal loop of improve. Let us assume that our assertions hold for iterations 1 through i − 1. Let (x^{−}, s^{−}) denote the value of (x, s) at the start of iteration i − 1, and let (u^{−}, v^{−}) be the search direction calculated during this iteration. We reserve the notation (x, s) for the current point at the start of iteration i. To obtain an estimate of ∥(x^{k}, s^{k}) − (x, s)∥, we note that our inductive hypothesis (6.4) implies that

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

Applying Lemma 5.6 to the step (u^{−}, v^{−}) taken during iteration i − 1, with χ replaced by C_{5} ∏_{l=0}^{i−2}(1 + C_{6}ρ^{l}), we deduce that

∥(x^{k}, s^{k}) − (x, s)∥
≤ ∥(x^{k}, s^{k}) − (x^{−}, s^{−})∥ + ∥(x^{−}, s^{−}) − (x, s)∥
≤ [C_{5} ∏_{l=0}^{i−2}(1 + C_{6}ρ^{l})] τ_{k} + ∥(u^{−}, v^{−})∥

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

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

The final inequality follows from the fact that τ^{−} ≤ ρ^{i−1}τ_{k}, since (x^{−}, s^{−}) is obtained by taking i − 1 fast steps ( one step in the main algorithm, followed by i − 2 iterations of the improve loop). We have shown that the bound (6.4) holds at iteration i as well. Since

C_{7}χ^{2} \frac{τ_{k}}{γ^{t}} = C_{7}χ^{2} \frac{τ_{k}}{γ^{t_{k}+i}} ≤ C_{7}χ^{2} \frac{τ_{K}}{γ^{t_{k}+I}} ≤ C_{7}χ^{2} \frac{τ_{K}}{γ^{t_{K}+I}} ≤ ρ,

Lemma 5.7 applies again, and therefore a fast step is accepted at iteration i. This completes our induction argument and proves our lemma. □

We end this section by showing that the sequence {μ_{k}} converges to zero with Q-order at least I + 2. Since γ_{min}ρ_{0}τ ≤ μ ≤ (1 + β_{max})τ, this is equivalent to show that for any ϵ > 0,

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

An equivalent characterization of the Q-order I + 2 convergence is

lim inf_{k→∞} \frac{ln τ_{k+1}}{ln τ_{k}} ≥ I + 2  (6.5)

(see Ortega and Rheinboldt [9] or Potra [10]).

Theorem 6.2 Suppose F^{c} is not empty, and λ during the safe branch is defined by (3.8), then the main sequence {μ_{k}}, k = 0, 1, …, converges to zero with Q-order at least I + 2.

Proof. We will show (6.5). Consider k ≥ K. By Lemma 6.1, a fast step is taken by the main algorithm and all I iterations of improve, so we can apply inequality (5.22) I + 1 times to bound τ_{k+1}:

τ_{k+1} ≤ \frac{C_{7}^{I+1}χ^{2(I+1)}τ_{k}^{I+1}}{γ^{(I+1)t_{k}+I(I+1)/2}} .  (6.6)

From Lemma 5.7 and (6.6) we deduce that

\frac{τ_{k+1}}{τ_{k}} ≤ C_{9} (\frac{τ_{k}}{γ^{t_{k}}} )^{I+1} → 0,  (6.7)

where

C_{9} := \frac{C_{7}^{I+1}χ^{2(I+1)}}{γ^{I(I+1)/2}} .

This means that {τ_{k}} (and therefore {μ_{k}}) converges to zero at least Q-superlinearly. By taking logarithms in (6.6), we get

ln τ_{k+1} ≤ ln C_{9} + (I + 2) ln τ_{k} − (I + 1)t_{k} ln γ.  (6.8)

We may assume that k is sufficiently large so that τ_{k} < 1. Then from (6.8) we have

\frac{ln τ_{k+1}}{ln τ_{k}} ≥ I + 2 + \frac{C_{9}}{ln τ_{k}} − ((I + 1) ln γ)\frac{t_{k}}{ln τ_{k}} .  (6.9)

Obviously, as k → ∞, the second term in the right-hand side vanishes. Since t_{k} ≤ (I + 1)k + 1, it suffices to prove

lim_{k→∞} \frac{k}{ln τ_{k}} = 0.  (6.10)

Suppose, by contradiction, that (6.10) does not hold. Then there exists a subsequence, also denoted by τ_{k}, such that there exists ν > 0,

\frac{k}{− ln τ_{k}} ≥ ν, for all k > 0,

which implies (1/e)^{k} ≤ τ_{k}^{ν}. From (6.7), there exists a positive integer J such that for all k ≥ J, τ_{k+1} ≤ (.5/e)^{1/ν}τ_{k}. Hence, for all k > J,

(1/e)^{k} ≤ τ_{k}^{ν} ≤ ((.5/e)^{(k−J)/ν}τ_{J})^{ν} = (.5/e)^{(k−J)}τ_{J}^{ν}.

It follows that for all k > J we must have 2^{k} ≤ (.5/e)^{−J}τ_{J}^{ν}, which is a contradiction. This completes the proof. □

7 Concluding Remarks

• Our algorithm can solve all P_{∗}-matrix linear complementarity problems , without knowing the classification number κ in advance. The computational complexity of the algorithm depends on κ.

• The algorithm is new even for monotone LCP where it attains O(nL)-iteration complexity, which improves the O(n^{2}L)-iteration complexity of the algorithm of [17].

• The algorithm requires only one matrix factorization at each iteration, compared to two matrix factorizations of predictor-corrector algorithms (cf. [8], [11], [12], [13]).

• Asymptotically, the algorithm requires I + 1 backsolves at each iteration to get fast local convergence of Q-order at least I+2 for problems having a strictly complementary solution.

• The algorithm can be easily extended to horizontal LCP and mixed LCP (see [1]).

Acknowledgment

The authors would like to thank Mihai Anitescu for carefully reading the manuscript.

References

[1] M. Anitescu, G. Lesaja, and F. A. Potra. On the equivalence between different formulations and algorithms for the P_{∗}(κ) − LCP. Reports on Computational Mathematics, Department of Mathematics, The University of Iowa, Iowa City, IA 52242, USA, October 1994.

[2] C. C. Gonzaga. The largest step path following algorithm for monotone linear complementary problems. Technical Reports, January 1994.

[3] A. J. Hoffman. On approximate solutions of systems of linear inequalities. Journal of Research of the National Bureau of Standards, 49:263–265, 1952.

[4] J. Ji and F. A. Potra. An infeasible-interior-point method for the P_{∗}-matrix LCP. Reports on Computational Mathematics 52, Department of Mathematics, University of Iowa, Iowa City, IA 52242, USA, February 1994.

[5] J. Ji, F. A. Potra, and R. Sheng. A predictor–corrector method for solving the P_{∗}- matrix LCP from infeasible starting points. Reports on Computational Mathematics 55, Department of Mathematics, University of Iowa, Iowa City, IA 52242, USA, June 1994.

[6] M. Kojima, N. Megiddo, T. Noma, and A. Yoshise. A unified approach to interior point algorithms for linear complementarity problems. Lecture Notes in Comput. Sci., 538, 1991.

[7] J. Miao. A quadratically convergent O((1 + κ)\sqrt{n}L)-iteration algorithm for the P_{∗}(κ)- matrix linear complementarity problem. Research Report RRR 93, RUTCOR-Rutgers Center for Operations Research, Rutgers University, P.O Box 5063, New Brunswick, NJ, USA, 1993.

[8] S. Mizuno, F. Jarre, and J. Stoer. A unified approach to infeasible-interior-point algorithms via

geometrical linear complementarity problems. Preprint Nr. 213, Mathematische Institute der Universität Würzburg, 97074 Würzburg, Germany, April 1994.

[9] J. M. Ortega and W. C. Rheinboldt. Iterative Solution of Nonlinear Equations in Several Variables. Academic Press, New York, 1970.

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

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

[12] F. A. Potra and R. Sheng. A modified O(nL) infeasible-interior-point algorithm for LCP with quadratic convergence. Reports on Computational Mathematics 54, Department of Mathematics, The University of Iowa, Iowa City, IA 52242, USA, April 1994.

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

[14] S. J. Wright. A path–following infeasible–interior–point algorithm for linear complementarity problems. Optimization Methods and Software, 2:79–106, 1993.

[15] S. J. Wright. A path–following interior–point algorithm for linear and quadratic problems. Preprint MCS–P401–1293, Mathematics and Computer Science Division, Argonne National Laboratory, Argonne, IL 60439, USA, December 1993.

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

[17] S. J. Wright and Y. Zhang. A superquadratic infeasible–interior–point algorithm for linear complementarity problems. Preprint MCS–P418–0294, Mathematics and Computer Science Division, Argonne National Laboratory, Argonne, IL 60439, USA, February 1994.

[18] Y. Zhang. On the convergence of a class of infeasible interior–point methods for the horizontal linear complementarity problem. SIAM J. Optimization, 4(1):208–227, 1994.