On the Superlinear Convergence of an O(n^{3}L) Interior Point Algorithm for Monotone LCP

Kevin A. McShane 2537 Villanova Drive Vienna, VA 22180

27 October 1994, revised 6 July 1995

Abstract

Karmarkar’s partial updating scheme is applied to a polynomial, superlinearly convergent algorithm for the monotone linear complementarity problem. This modification reduces the bound on the number of arithmetic operations required to achieve ϵ-optimality by a factor of \sqrt{n}. The complementarity gap sequence generated by the resulting algorithm converges q-quadratically to zero if the LCP is nondegenerate and has a strictly complementary optimal solution. In general, the convergence rate is q-superlinear if the iterates are assumed to converge to a strictly complementary optimal solution.

1 Introduction

We examine a partial updating variant of a polynomial, superlinearly convergent algorithm for the monotone linear complementarity problem (LCP): find (x, s) ∈ R^{n} × R^{n} which solves

minimize x′s
subject to Mx + q = s

x, s ≥ 0 where M ∈ R^{n×n} is a positive semidefinite (not necessarily symmetric) matrix and q ∈ R^{n}. We assume that the set

F^{++} ≐ {(x, s) ∈ R^{n} |Mx + q = s, x > 0, s > 0}

is nonempty. The optimal objective function value for LCP is zero, hence the optimal values for x and s are complementary, i.e. orthogonal and nonnegative. We call the LCP nondegenerate if it has a unique solution, and degenerate otherwise. LCP is a general class of problems which includes as special cases the linear programming (LP) and convex quadratic programming problems.

The modified LCP algorithm described below has the following properties:

• The bound on the number of arithmetic operations required to achieve ϵ-optimality is a factor of \sqrt{n} smaller than that of the original algorithm. If the data is integer, the computational complexity of the variant is O(n^{3}L), where L is the encoding length of the problem.

• The complementarity gap values (the objective function values) converge to zero at least q-linearly. If the LCP is nondegenerate with a strictly complementary optimal solution (i.e. (x^{∗}, s^{∗}) is feasible for LCP, x^{∗′}s^{∗} = 0, and x^{∗} + s^{∗} > 0), the complementarity gap sequence converges to zero q-quadratically.

• If the problem is degenerate, the convergence rate is q-superlinear under the assumption that the variables converge to a strictly complementary optimal solution.

The use of the partial updating technique to improve the computational complexity of interior point methods originated in Karmarkar’s landmark paper [10]. Karmarkar showed that by updating the current factorization via rank-one updates (based upon approximate values for the variables) rather than performing a complete factorization each iteration, the computational complexity of his algorithm could be reduced by a factor of \sqrt{n}. The technique has since been applied to a variety of other interior point algorithms for LP and LCP to reduce their complexity bounds (for example, [1], [2], [3], [5], [6], [7], [8], [12], [13], [14], [15], [17], [19]).

Our algorithm is based on what Gonzaga [9] refers to as the “longest step path-following” (LSPF) algorithm. The LSPF algorithm was first examined in Mizuno-Yoshise-Kikuchi [16], where it was shown to require at most O(\sqrt{n}L) iterations to solve LCP with integer data. Each iteration requires O(n^{3}) arithmetic operations, so the computational complexity of the algorithm is O(n^{3.5}L). Mizuno and Todd [15] later proposed a partial updating variant of the LSPF algorithm which reduces the computational complexity to O(n^{3}L). The algorithm we propose is different from the Mizuno-Todd algorithm because it uses a different partial updating scheme, does not use the auxiliary sequence of parameters μ_{k} (related to the complementarity gap), and has superlinear convergence properties. To the best of our knowledge, this is the first superlinear convergence result for an O(n^{3}L) algorithm for LP or LCP.

In [11], McShane showed that the complementarity gap sequence generated by the LSPF algorithm converges to zero q-quadratically if the LCP is nondegenerate and has a strictly complementary optimal solution. In the degenerate case, the convergence rate was shown to be q-superlinear, assuming the variables converge to a strictly complementary optimal solution. Bonnans and Gonzaga [4] followed by showing that the iterates do indeed converge in the degenerate case, establishing q-superlinear convergence assuming only the existence of a strictly complementary optimal solution. Recently in [9], Gonzaga has shown that the convergence rate is q-quadratic if a computationally trivial safeguard is applied.

Now, a few notational conventions. The norm ∥ · ∥ refers to the l_{2}-norm unless denoted otherwise by a subscript, and log(·) refers to the natural logarithm. Given a vector v, a subscript of i or j denotes a particular component, while a superscript (usually k) generally refers to an iteration number. A capital letter (say, X or V) version of a vector (x or v) refers to the diagonal matrix with i-th diagonal entry the i-th component of the vector. We may also refer to X as diag(x). The symbol “ ≐ ” assigns a definition to a quantity. Given a sequence of vectors {x^{k}} and scalars {π_{k}}, the notation x = O(π) means that there exists a constant C (dependent on problem data) such that ∥x^{k}∥ ≤ Cπ_{k} for all values of k.

2 Partial Updating

The majority of the work required each iteration in many interior point algorithms is in the formation and factorization of a matrix of the form ADA′, or the form A+D, where A is an n×m (or n×n) matrix which does not change from iteration to iteration, and D is a diagonal matrix which may change from iteration to iteration. These matrix operations have a computational complexity of O(n^{3}). The remainder of the work each iteration is typically O(n^{2}). Karmarkar showed that by using an approximation D̃ for D, one could update the old factorization via rank-one updates–one update for each entry of D̃ which changed that iteration. Each rank-one update requires O(n^{2}) arithmetic operations, so the number of arithmetic operations performed is O(n^{2}K), where K is the number of rank-one updates done, or equivalently, the number of times entries of D̃ change.

Our algorithm maintains approximation vectors x̃ for x and s̃ for s which satisfy

\frac{1}{δ} ≤ \frac{x_{i}}{x̃_{i}} ≤ δ and \frac{1}{δ} ≤ \frac{s_{i}}{s̃_{i}} ≤ δ  (1)

for a fixed constant δ ≥ 1 by using the following generic scheme. We assume D̃_{ii} is a function of x̃_{i} and s̃_{i} alone.

The partial updating procedure.

STEP 1: Given (x^{k}, s^{k}) ∈ F^{++}and (x̃^{k}, s̃^{k}) which satisfy (1), compute D̃^{k}, x^{k+1}, and s^{k+1};

For i = 1 to n,

if ( \frac{x_{i}^{k+1}}{x̃_{i}^{k}} ≥ δ or \frac{x_{i}^{k+1}}{x̃_{i}^{k}} ≤ \frac{1}{δ} or \frac{s_{i}^{k+1}}{s̃_{i}^{k}} ≥ δ or \frac{s_{i}^{k+1}}{s̃_{i}^{k}} ≤ \frac{1}{δ} )

then set x̃_{i}^{k+1} = x_{i}^{k+1} and s̃_{i}^{k+1} = s_{i}^{k+1};
else set x̃_{i}^{k+1} = x̃_{i}^{k} and s̃_{i}^{k+1} = s̃_{i}^{k};

STEP 2: If the termination criteria are not met, set k = k+1, perform a rank-one update to the old factorization for each changed entry of D̃, go to STEP 1.

The following lemma appears as Lemma 3.3 in Anstreicher [1], and the subsequent theorem is based on ideas contained in the proofs of Lemma 3.1, Corollary 3.2 and Theorem 3.4 in [1]:

Lemma 1 If 0 < a < 1 and b ≥ a, then |log(b)| ≤ \frac{|b − 1| log(a)}{a − 1}.

Proof. The lemma follows from the fact that the function

f(x) = \{ \begin{matrix}log(x)/(x − 1) & for & x > 0, x ≠ 1 \\ 1 & for & x = 1 \end{matrix}

is continuous and monotonically decreasing.

Theorem 2 Suppose an algorithm for LCP generates a sequence of feasible interior points (x^{0}, s^{0}), (x^{1}, s^{1}),…,(x^{K}, s^{K}) via x^{k+1} = x^{k} + d_{x}^{k} and s^{k+1} = s^{k} + d_{s}^{k}. Suppose also that the factorization used to compute the search directions is efficiently updated by using the partial updating procedure with some δ > 1. Let T represent the number of rank-one updates performed. Then

T ≤ \frac{δ}{δ − 1} ∑_{k=0}^{K−1} (∥(X^{k})^{−1}d_{x}^{k}∥_{1} + ∥(S^{k})^{−1}d_{s}^{k}∥_{1}).

Proof. Define T^{i} as the ordered set of iteration indices k such that x̃_{i}^{k} ≠ x̃_{i}^{k−1}. Insert 0 at the beginning of T^{i}. Let a and b, with a < b, be two consecutive elements of T^{i}. Note that x̃_{i}^{b−1} = x̃_{i}^{a} = x_{i}^{a} and s̃_{i}^{b−1} = s̃_{i}^{a} = s_{i}^{a} since x̃_{i} and s̃_{i} are not updated during the iterations between a and b. Since x̃_{i}^{b} ≠ x̃_{i}^{b−1}, we have

\frac{x_{i}^{b}}{x_{i}^{a}} ≥ δ or \frac{x_{i}^{b}}{x_{i}^{a}} ≤ \frac{1}{δ} or \frac{s_{i}^{b}}{s_{i}^{a}} ≥ δ or \frac{s_{i}^{b}}{s_{i}^{a}} ≤ \frac{1}{δ} .

In any case,

log(δ) ≤ |log (\frac{x_{i}^{b}}{x_{i}^{a}} )| + |log (\frac{s_{i}^{b}}{s_{i}^{a}} )| ,

so

log(δ) ≤ |log (\frac{x_{i}^{a+1}}{x_{i}^{a}} \frac{x_{i}^{a+2}}{x_{i}^{a+1}} ⋯ \frac{x_{i}^{b}}{x_{i}^{b−1}} )| + |log (\frac{s_{i}^{a+1}}{s_{i}^{a}} \frac{s_{i}^{a+2}}{s_{i}^{a+1}} ⋯ \frac{s_{i}^{b}}{s_{i}^{b−1}} )|
= |∑_{k=a}^{b−1} log (\frac{x_{i}^{k+1}}{x_{i}^{k}} )| + |∑_{k=a}^{b−1} log (\frac{s_{i}^{k+1}}{s_{i}^{k}} )|

≤ ∑_{k=a}^{b−1} (|log (\frac{x_{i}^{k+1}}{x_{i}^{k}} )| + |log (\frac{s_{i}^{k+1}}{s_{i}^{k}} )|) .

Define

q_{i}^{k} ≐ max \{\frac{x_{i}^{k+1}}{x_{i}^{k}} , \frac{1}{δ}\} ,
r_{i}^{k} ≐ max \{\frac{s_{i}^{k+1}}{s_{i}^{k}} , \frac{1}{δ}\} .

Since log(δ) = |log (\frac{1}{δ})|, we have

log(δ) ≤ ∑_{k=a}^{b−1} (|log(q_{i}^{k})| + |log(r_{i}^{k})|) .

Now, sum over all values for a ∈ T^{i} and over all values for i:

T log(δ) ≤ ∑_{i=1}^{n} ∑_{k=0}^{K−1} (|log(q_{i}^{k})| + |log(r_{i}^{k})|)

where T is the total number of updates done to x̃ (or s̃). Applying Lemma 1 twice, with 1/δ for a and q_{i}^{k} for b, then with 1/δ for a and r_{i}^{k} for b yields

T log(δ) ≤ ∑_{i=1}^{n} ∑_{k=0}^{K−1} \frac{δ log(δ)}{δ − 1} (|1 − q_{i}^{k}| + |1 − r_{i}^{k}|)

≤ ∑_{i=1}^{n} ∑_{k=0}^{K−1} \frac{δ log(δ)}{δ − 1} (|1 − \frac{x_{i}^{k+1}}{x_{i}^{k}} | + |1 − \frac{s_{i}^{k+1}}{s_{i}^{k}} |)

= \frac{δ log(δ)}{δ − 1} ∑_{k=0}^{K−1} (∥(X^{k})^{−1}d_{x}^{k}∥_{1} + ∥(S^{k})^{−1}d_{s}^{k}∥_{1}) ,

which completes the proof.

3 The Algorithm

The search directions we use are based on the LSPF directions and are quite similar to the Mizuno-Todd directions. Given (x, s) ∈ F^{++} and approximation vectors x̃ and s̃, define

Δ ≐ \frac{x′s}{n},

v̄ ≐ −(X̃S̃)^{−1/2}(XSe − Δe),
v̂ ≐ −(X̃S̃)^{−1/2}(XS)e, and
D̃ ≐ X̃^{1/2}S̃^{−1/2}.

Here, e is the n-vector of ones. Our directions are a convex combination of the “centering” and “affine” directions

( d̄_{x} ; d̄_{s} ) ≐ ( \begin{matrix}M & −I \\ D̃^{−1} & D̃ \end{matrix} )^{−1} ( 0 ; v̄ )  (2)
( d̂_{x} ; d̂_{s} ) ≐ ( \begin{matrix}M & −I \\ D̃^{−1} & D̃ \end{matrix} )^{−1} ( 0 ; v̂ )  (3)

that is, our directions are of the form

( \begin{matrix}d_{x}(θ) \\ d_{s}(θ) \end{matrix} ) = θ ( d̄_{x} ; d̄_{s} ) + (1 − θ) ( d̂_{x} ; d̂_{s} )  (4)

where θ ∈ [0, 1]. It is easily verified that the inverse matrix in (2) exists. We also define

D_{x}(θ) ≐ diag(d_{x}(θ)),
D_{s}(θ) ≐ diag(d_{s}(θ)),
x(θ) ≐ x + d_{x}(θ),
s(θ) ≐ s + d_{s}(θ),
Δ(θ) ≐ \frac{x(θ)′s(θ)}{n}, and
θ^{∗} ≐ 1 − \frac{1}{5\sqrt{n}}.

Note that from (3),

Md̂_{x} = d̂_{s},  (5)

and

d̂_{s} = D̃^{−1}v̂ − D̃^{−2}d̂_{x}

so by substitution, d̂_{x} can be formed by solving

(M + D̃^{−2})d̂_{x} = D̃^{−1}v̂,  (6)

then d̂_{s} can be computed from (5). A similar argument holds for d̄_{x} and d̄_{s}, so the partial updating scheme can be efficiently applied to a factorization of the matrix M + D̃^{−2}.

Like the LSPF algorithm, our algorithm is a path-following algorithm. All iterates (x^{k}, s^{k}) are required to satisfy the approximately-centered condition

∥X^{k}S^{k}e − Δ_{k}e∥ ≤ βΔ_{k}  (7)

where β ∈ (0, 1) is a fixed parameter.

We now propose our algorithm:

The LCP algorithm.

STEP 0: Set β = .2, β̂ = β/(2 + β), δ = 1.01, θ^{∗} = 1 − 1/(5\sqrt{n}), and k = 0. Given ε > 0 and (x^{0}, s^{0}) ∈ F^{++} which satisfies (7), set x̃^{0} = x^{0} and s̃^{0} = s^{0};

STEP 1: Given (x^{k}, s^{k}) ∈ F^{++} which satisfies (7), and x̃^{k} and s̃^{k} generated by the partial updating procedure, compute d̄_{x}^{k}, d̄_{s}^{k}, d̂_{x}^{k} and d̂_{s}^{k} using the method in (6). Calculate the roots of the quartic polynomial

g(θ) ≐ ∥X^{k}(θ)S^{k}(θ)e − θΔ_{k}e∥^{2} − β̂^{2}θ^{2}(Δ_{k})^{2}.  (8)

If g has no roots in the interval (0, θ^{∗}], set θ_{k} = 0 and compute the search directions from (4). The updated points will be optimal for LCP. Otherwise, let θ_{k} be the largest root of g in the interval (0, θ^{∗}];

STEP 2: Set x^{k+1} = x^{k}(θ_{k}) and s^{k+1} = s^{k}(θ_{k}). If (x^{k+1})′(s^{k+1}) ≤ ε then terminate, else perform the partial updating procedure, set k = k + 1, then go to STEP 1.

4 Analysis of the Algorithm

The function g(θ) in (8) is similar to the functions used to compute quantities equivalent to θ in the Mizuno-Todd algorithm and the LSPF algorithm as described in Gonzaga [9]. These papers replace Δ_{k} in (7) with a related parameter μ_{k}. The relationship between μ_{k} and Δ_{k} is described in detail in [9]. We choose not to monitor μ, but instead to monitor Δ, in order to exploit the orthogonality of e and XSe − Δe. For our purposes, the function g(θ) is simply an approximation to the updated approximate centering condition

∥X^{k}(θ)S^{k}(θ)e − Δ_{k}(θ)e∥ ≤ βΔ_{k}(θ)

which simplifies the analysis of the algorithm. In the analysis below, we suppress the iteration index to simplify the notation.

Theorem 3 If during an iteration of the LCP algorithm, θ ∈ [0, 1] is chosen so that g(θ) ≤ 0, then ∥X(θ)S(θ)e − Δ(θ)e∥ ≤ βΔ(θ).

Proof. Note that

∥X(θ)S(θ)e − Δ(θ)e∥ ≤ ∥X(θ)S(θ)e − θΔe∥ + ∥(θΔ − Δ(θ))e∥ .

But g(θ) ≤ 0 implies ∥X(θ)S(θ)e − θΔe∥ ≤ β̂θΔ, so

∥X(θ)S(θ)e − Δ(θ)e∥ ≤ β̂θΔ + |θΔ − Δ(θ)| \sqrt{n}.  (9)

Following Gonzaga [9], define

h ≐ X(θ)S(θ)e − θΔe,

hence ∥h∥ ≤ β̂θΔ. If we left-multiply the above expression by (1/n)e′, the result is

Δ(θ) − θΔ = \frac{1}{n}e′h  (10)

thus

|Δ(θ) − θΔ| ≤ \frac{1}{n} ∥e∥ ∥h∥ ≤ \frac{1}{\sqrt{n}}β̂θΔ.  (11)

Also, from (10),

θΔ = Δ(θ) − \frac{1}{n}e′h ≤ Δ(θ) + \frac{1}{\sqrt{n}}β̂θΔ ≤ Δ(θ) + β̂θΔ

so therefore

θΔ ≤ \frac{1}{1 − β̂}Δ(θ).  (12)

If we apply inequalities (11) and (12) to (9), we see that

∥X(θ)S(θ)e − Δ(θ)e∥ ≤ \frac{2β̂}{1 − β̂}Δ(θ) = βΔ(θ)

which completes the proof.

The goal of this section is to show that the LCP algorithm is well-defined. The following three lemmas are similar to results used in the analysis of a variety of path-following O(n^{3}L) interior point algorithms, [8] and [18] for example. The proofs are included for the sake of clarity.

Lemma 4 During each iteration of the LCP algorithm, the following relationships hold:

a) ∥D̃(S − S̃)∥ ≤ δ(δ − 1)\sqrt{1 + β}\sqrt{Δ}

b) ∥D̃^{−1}(X − X̃)∥ ≤ δ(δ − 1)\sqrt{1 + β}\sqrt{Δ}. Proof. Observe that

∥D̃(S − S̃)∥ = ∥(X̃S̃)^{1/2}(S̃^{−1}S − I)∥
≤ ∥(X̃S̃)^{1/2}∥ ∥S̃^{−1}S − I∥ .  (13)

But (7) implies that

(1 − β)Δ ≤ x_{i}s_{i} ≤ (1 + β)Δ  (14)

so, using (1),

x̃_{i}s̃_{i} ≤ δ^{2}(1 + β)Δ

thus

∥(X̃S̃)^{1/2}∥ ≤ δ\sqrt{1 + β}\sqrt{Δ}.

Since (1) implies that

∥S̃^{−1}S − I∥ ≤ δ − 1,

with (13) we have

∥D̃(S − S̃)∥ ≤ δ(δ − 1)\sqrt{1 + β}\sqrt{Δ}.

Inequality b) can be shown similarly.

Define

ṽ(θ) ≐ θv̄ + (1 − θ)v̂.  (15)

From (2), (3), and (4) we see that

D̃^{−1}d_{x}(θ) + D̃d_{s}(θ) = ṽ(θ).  (16)

Lemma 5 During each iteration of the LCP algorithm, the following relationships hold for any θ ∈ [0, 1]:

a) ∥D̃^{−1}d_{x}∥^{2} + ∥D̃d_{s}∥^{2} ≤ \left∥{ \atop ṽ}\right∥^{2}

b) ∥D̃^{−1}d_{x}∥ + ∥D̃d_{s}∥ ≤ \sqrt{2} \left∥{ \atop ṽ}\right∥

c) ∥D_{x}D_{s}e∥ ≤ \frac{1}{2} \left∥{ \atop ṽ}\right∥^{2}

d) ∥(X − X̃)d_{s} + (S − S̃)d_{x}∥ ≤ \sqrt{2}δ(δ − 1)\sqrt{1 + β}\sqrt{Δ} \left∥{ \atop ṽ}\right∥,

where d_{x} ≐ d_{x}(θ), d_{s} ≐ d_{s}(θ), D_{x} ≐ D_{x}(θ), D_{s} ≐ D_{s}(θ), and ṽ ≐ ṽ(θ). Proof. First, a). From (16),

∥D̃^{−1}d_{x} + D̃d_{s}∥^{2} = \left∥{ \atop ṽ}\right∥^{2} ,

and expanding the leftmost norm and rearranging yields

∥D̃^{−1}d_{x}∥^{2} + ∥D̃d_{s}∥^{2} ≤ \left∥{ \atop ṽ}\right∥^{2} − d′_{x}d_{s}.

But from (2), (3), and (4),

d′_{x}d_{s} = d′_{x}Md_{x} ≥ 0

because M is positive semidefinite. Therefore,

∥D̃^{−1}d_{x}∥^{2} + ∥D̃d_{s}∥^{2} ≤ \left∥{ \atop ṽ}\right∥^{2}

which completes the proof of a).

Inequality b) follows from a) and the fact that

(∥D̃^{−1}d_{x}∥ + ∥D̃d_{s}∥)^{2} ≤ 2 (∥D̃^{−1}d_{x}∥^{2} + ∥D̃d_{s}∥^{2}) .

To show c), note that

∥D_{x}D_{s}e∥^{2} = (e′D_{x}^{2}D̃^{−2})(D̃^{2}D_{s}^{2}e)
≤ ∥D̃^{−2}D_{x}^{2}e∥ ∥D̃^{2}D_{s}^{2}e∥
≤ ∥D̃^{−2}D_{x}^{2}e∥_{1} ∥D̃^{2}D_{s}^{2}e∥_{1}
= ∥D̃^{−1}D_{x}e∥^{2} ∥D̃D_{s}e∥^{2}
≤ \frac{1}{4} (∥D̃^{−1}d_{x}∥^{2} + ∥D̃d_{s}∥^{2})^{2}

by the arithmetic-geometric mean inequality. Inequality c) follows by applying a) to this expression.

Finally,

∥(X − X̃)d_{s} + (S − S̃)d_{x}∥ ≤ ∥(X − X̃)d_{s}∥ + ∥(S − S̃)d_{x}∥

= ∥(X − X̃)D̃^{−1}D̃d_{s}∥ + ∥(S − S̃)D̃D̃^{−1}d_{x}∥
≤ ∥D̃^{−1}(X − X̃)∥ ∥D̃d_{s}∥ + ∥D̃(S − S̃)∥ ∥D̃^{−1}d_{x}∥ .

Inequality d) follows by applying Lemma 4, then Lemma 5 b).

Lemma 6 During each iteration of the LCP algorithm, the following relationships hold for any θ ∈ [0, 1]:

a) ∥ṽ(θ)∥^{2} ≤ δ^{2} (β^{2}/(1 − β) + (1 − θ)^{2}n) Δ

b) ∥ṽ(θ)∥ ≤ δ (β/\sqrt{1 − β} + (1 − θ)\sqrt{n}) \sqrt{Δ}. Proof. From (15),

ṽ(θ) = θv̄ + (1 − θ)v̂
= −(X̃S̃)^{−1/2}(XS)^{1/2} (θ(XS)^{−1/2}(XSe − Δe) + (1 − θ)(XS)^{1/2}e) ,

therefore (using(1))

∥ṽ(θ)∥^{2} ≤ δ^{2} ∥θ(XS)^{−1/2}(XSe − Δe) + (1 − θ)(XS)^{1/2}e∥^{2} .

Since the two terms above are orthogonal, and since θ ∈ [0, 1], we have

∥ṽ(θ)∥^{2} ≤ δ^{2} (∥(XS)^{−1/2}(XSe − Δe)∥^{2} + (1 − θ)^{2} ∥(XS)^{1/2}e∥^{2})
≤ δ^{2} (∥(XS)^{−1}∥ ∥XSe − Δe∥^{2} + (1 − θ)^{2}nΔ)
≤ δ^{2} (\frac{1}{(1 − β)Δ}β^{2}Δ^{2} + (1 − θ)^{2}nΔ) (by (14) and (7))
= δ^{2} (\frac{β^{2}}{1 − β} + (1 − θ)^{2}n) Δ, which proves a).

Inequality b) follows trivially from a). Theorem 7 During STEP 1 of each iteration of the LCP algorithm, for any θ ∈ [θ^{∗}, 1], (x(θ), s(θ)) ∈ F^{++}. Also, g(θ^{∗}) < 0. Proof. Define d_{x} ≐ d_{x}(θ), d_{s} ≐ d_{s}(θ), D_{x} ≐ D_{x}(θ), D_{s} ≐ D_{s}(θ), and ṽ ≐ ṽ(θ).

If we multiply equation (16) by (X̃S̃)^{\frac{1}{2}} , the result is

X̃d_{s} + S̃d_{x} = (X̃S̃)^{1/2}ṽ = θΔe − XSe.  (17)

Consider the expression in the first term in the definition of g(θ): X(θ)S(θ)e − θΔe = XSe + Xd_{s} + Sd_{x} + D_{x}D_{s}e − θΔe

= XSe + X̃d_{s} + S̃d_{x} + D_{x}D_{s}e − θΔe + (X − X̃)d_{s} + (S − S̃)d_{x}
= D_{x}D_{s}e + (X − X̃)d_{s} + (S − S̃)d_{x} (by (17)).  (18)

Therefore,

∥X(θ)S(θ)e − θΔe∥ ≤ ∥D_{x}D_{s}e∥ + ∥(X − X̃)d_{s} + (S − S̃)d_{x}∥ so applying Lemmas 5 and 6 results in

∥X(θ)S(θ)e − θΔe∥ ≤ \frac{1}{2} \left∥{ \atop ṽ}\right∥^{2} + \sqrt{2}δ(δ − 1)\sqrt{1 + β}\sqrt{Δ} \left∥{ \atop ṽ}\right∥

≤ \frac{δ^{2}}{2} (\frac{β^{2}}{1 − β} + (1 − θ)^{2}n) Δ + α (\frac{β}{\sqrt{1 − β}} + (1 − θ)\sqrt{n}) Δ (19)

where α ≐ \sqrt{2}δ^{2}(δ − 1)\sqrt{1 + β}.

Now, for any θ ∈ [θ^{∗}, 1], and with the choices β = .2, β̂ = β/(2 + β), δ = 1.01, and θ^{∗} = 1 − 1/(5\sqrt{n}), it is easily verified that (19) is bounded above by .06Δ. Therefore, for each index i,

−.06Δ ≤ x_{i}(θ)s_{i}(θ) − θΔ ≤ .06Δ

so

−.06Δ + θΔ ≤ x_{i}(θ)s_{i}(θ).

Moreover, since θ^{∗} ≥ .8,

.7Δ ≤ x_{i}(θ)s_{i}(θ).

Consequently, by continuity, (x(θ), s(θ)) ∈ F^{++} for any θ ∈ [θ^{∗}, 1]. It is easily verified that .06 < β̂θ^{∗}, which implies g(θ^{∗}) < 0.

Corollary 8 The LCP algorithm is well-defined.

Proof. An argument analogous to that in Theorem 7 shows that as long as g(θ) remains less than zero, (x(θ), s(θ)) ∈ F^{++} as θ decreases from the value θ^{∗}. If g has no roots between 0 and θ^{∗}, by continuity, g(0) ≤ 0 and (x(0), s(0)) ≥ 0. But from the definition of g, g(0) ≤ 0 implies that ∥X(0)S(0)e∥ = 0, therefore (x(0), s(0)) solves LCP. Otherwise, if g has a root in the interval [0, θ^{∗}], the choice for θ as the largest root in the interval results in strictly feasible updated iterates which, by Theorem 3, are approximately centered.

5 Complexity of the Algorithm

Theorem 9 Given a value for ε ∈ (0, 1), if x^{0′}s^{0} ≤ 1/ε, then x^{K′}s^{K} ≤ ε after at most K = O (\sqrt{n} log(1/ε)) iterations. The number of arithmetic operations required is O (n^{3} log(1/ε)). If the data is integer (or rational), the algorithm requires O(\sqrt{n}L) iterations and O(n^{3}L) arithmetic operations, where L is the encoding length of the problem.

Proof. For each iteration k, from (10),

Δ_{k} = Δ_{k−1}(θ_{k−1}) = \frac{1}{n}e′h^{k−1} + θ_{k−1}Δ_{k−1}

thus with (11) we have

Δ_{k} ≤ (1 + \frac{β̂}{\sqrt{n}}) θ_{k−1}Δ_{k−1}.

If we multiply this expression by n, the result is

(x^{k})′(s^{k}) ≤ (1 + \frac{β̂}{\sqrt{n}}) θ_{k−1}(x^{k−1})′(s^{k−1})  (20)
≤ (x^{0})′(s^{0}) ∏_{j=0}^{k−1} ((1 + \frac{1}{10\sqrt{n}}) θ_{j}) .  (21)

But

θ_{j} ≤ θ^{∗} = 1 − \frac{1}{5\sqrt{n}}

thus

(x^{k})′(s^{k}) ≤ (x^{0})′(s^{0}) ∏_{j=0}^{k−1} ((1 + \frac{1}{10\sqrt{n}}) (1 − \frac{1}{5\sqrt{n}}))

≤ (1 − \frac{1}{10\sqrt{n}})^{k} (x^{0})′(s^{0}).  (22)

Define K_{1} to be the largest value for k such that (x^{k})′(s^{k}) > ε. Similarly, define K_{2} as the largest k such that (21) > ε, and K_{3} as the largest k such that (22) > ε. A standard argument shows that (22) implies K_{3} = O(\sqrt{n} log(1/ε)). From these definitions it is evident that

K_{1} ≤ K_{2} ≤ K_{3} = O (\sqrt{n} log(1/ε)) .  (23)

Define d_{x}^{k} ≐ d_{x}^{k}(θ_{k}), and d_{s}^{k} ≐ d_{s}^{k}(θ_{k}). To determine the algorithm’s computational complexity, we appeal to Theorem 2. The number of rankone updates required before the last iteration is

T ≤ \frac{δ}{δ − 1} ∑_{k=0}^{K_{1}−1} (∥(X^{k})^{−1}d_{x}^{k}∥_{1} + ∥(S^{k})^{−1}d_{s}^{k}∥_{1})

≤ \frac{δ\sqrt{n}}{δ − 1} ∑_{k=0}^{K_{1}−1} (∥(X^{k})^{−1}d_{x}^{k}∥_{2} + ∥(S^{k})^{−1}d_{s}^{k}∥_{2}) .  (24)

Note that

∥(X^{k})^{−1}d_{x}^{k}∥ = ∥(X^{k})^{−1}X̃^{k}(X̃^{k}S̃^{k})^{− \frac{1}{2}} (D̃^{k})^{−1}d_{x}^{k}∥
≤ ∥(X^{k})^{−1}X̃^{k}∥ ∥(X̃^{k}S̃^{k})^{− \frac{1}{2}} ∥ ∥(D̃^{k})^{−1}d_{x}^{k}∥ ,

so, using (1), (7), and Lemma 5 a),

∥(X^{k})^{−1}d_{x}^{k}∥ ≤ (δ) (\frac{δ}{\sqrt{(1 − β) Δ_{k}}} ) ∥ṽ^{k}(θ)∥

and by Lemma 6

∥(X^{k})^{−1}d_{x}^{k}∥ ≤ (\frac{δ^{3}}{\sqrt{1 − β}} ) (\frac{β}{\sqrt{1 − β}} + (1 − θ_{k})\sqrt{n}) .

A similar argument shows that the same bound applies to ∥(S^{k})^{−1}d_{s}^{k}∥, thus from (24),

T ≤ \frac{2δ^{4}\sqrt{n}}{(δ − 1) \sqrt{1 − β}} ∑_{k=0}^{K_{1}−1} (\frac{β}{\sqrt{1 − β}} + (1 − θ_{k})\sqrt{n}) .

Using the well-known inequality x ≤ − log(1 − x) for all x ∈ [0, 1), we see that (1 − θ_{k}) ≤ − log(θ_{k}), therefore

T ≤ \frac{2δ^{4}\sqrt{n}}{(δ − 1) \sqrt{1 − β}} ∑_{k=0}^{K_{1}−1} (\frac{β}{\sqrt{1 − β}} − log(θ_{k})\sqrt{n})
≤ \frac{2δ^{4}β\sqrt{n}}{(δ − 1)(1 − β)}K_{1} − \frac{2δ^{4}n}{(δ − 1) \sqrt{1 − β}} ∑_{k=0}^{K_{1}−1} log(θ_{k})
≤ O (n log(1/ε)) − \frac{2δ^{4}n}{(δ − 1) \sqrt{1 − β}} ∑_{k=0}^{K_{1}−1} log(θ_{k}) (by (23)). (25)

Now, from the definition of K_{2},

(x^{0})′(s^{0}) ∏_{j=0}^{K_{2}−1} ((1 + \frac{1}{10\sqrt{n}}) θ_{j}) > ε.

Therefore,

log ((x^{0})′(s^{0})) + K_{2} log (1 + \frac{1}{10\sqrt{n}}) + ∑_{j=0}^{K_{2}−1} log(θ_{j}) > log(ε),

hence by rearrangement, and by using the fact that log(1+x) ≤ x for x > −1 we see that

− ∑_{j=0}^{K_{2}−1} log(θ_{j}) < log ((x^{0})′(s^{0})) − log(ε) + \frac{K_{2}}{10\sqrt{n}}

so with the hypotheses of the theorem and (23) we have

− ∑_{j=0}^{K_{2}−1} log(θ_{j}) = O (log(1/ε)) .

Since K_{2} ≥ K_{1}, substituting this into (25) implies that the number of updates required before the last iteration is

T = O (n log(1/ε)) .

Since the last iteration requires O(n^{3}) arithmetic operations, and since each rank-one update requires O(n^{2}) arithmetic operations, the total number of arithmetic operations required is O (n^{3} log(1/ε)).

In the integer case, it is well known that ε = 2^{−2L} suffices to solve LCP, from which the remainder of the theorem follows.

Theorem 10 If the iterates (x^{k}, s^{k}) of the LCP algorithm converge to a strictly complementary optimal solution (x^{∗}, s^{∗}), then {(x^{k})′(s^{k})} → 0 qsuperlinearly. Moreover, after a finite number of iterations, every subsequent iteration will involve a complete factorization.

Proof. For each iteration k, define d_{x}^{k} ≐ d_{x}^{k}(θ_{k}), d_{s}^{k} ≐ d_{s}^{k}(θ_{k}), D_{x}^{k} ≐ D_{x}^{k}(θ_{k}), D_{s}^{k} ≐ D_{s}^{k}(θ_{k}), and ṽ^{k} ≐ ṽ^{k}(θ_{k}). If the algorithm does not terminate finitely, g(θ_{k}) = 0 for each iteration k, so from the definition of g,

∥X^{k}(θ_{k})S^{k}(θ_{k})e − θ_{k}Δ_{k}e∥ = β̂θ_{k}Δ_{k},

that is, by (18),

∥D_{x}^{k}D_{s}^{k}e + (X^{k} − X̃^{k})d_{s}^{k} + (S^{k} − S̃^{k})d_{x}^{k}∥ = β̂θ_{k}Δ_{k}

or by rearrangement,

θ_{k} = \frac{∥D_{x}^{k}D_{s}^{k}e + (X^{k} − X̃^{k})d_{s}^{k} + (S^{k} − S̃^{k})d_{x}^{k}∥}{β̂Δ_{k}}.  (26)

In order to show q-superlinear convergence of the complementarity gap sequence to zero, (20) implies that it suffices to show that θ_{k} → 0 as k → ∞. Due to (26), it is enough to show that

\frac{∥D_{x}^{k}D_{s}^{k}e∥}{Δ_{k}} → 0  (27)

and

\frac{∥(X^{k} − X̃^{k})d_{s}^{k} + (S^{k} − S̃^{k})d_{x}^{k}∥}{Δ_{k}} → 0.  (28)

We first show (27). The proof is similar to that of Theorem 4.4 in [11]. Consider component i of x^{k} and s^{k}. Assume without loss of generality

that x_{i}^{k} → x_{i}^{∗} > 0 and s_{i}^{k} → 0 as k → ∞. Since x_{i}^{k} + (d_{x}^{k})_{i} → x_{i}^{∗} implies x_{i}^{k}/x_{i}^{∗} + (d_{x}^{k})_{i}/x_{i}^{∗} → 1 as k → ∞, it follows that

(d_{x}^{k})_{i} → 0 as k → ∞.  (29)

Let K be such that for k ≥ K,

\frac{1}{2}x_{i}^{∗} ≤ x_{i}^{k} ≤ 2x_{i}^{∗}.  (30)

For k ≥ K, (7) and (30), followed by (1) imply

\frac{1 − β}{2δx_{i}^{∗}} ≤ \frac{s̃_{i}^{k}}{Δ_{k}} ≤ \frac{2(1 + β)δ}{x_{i}^{∗}} .  (31)

Now, Lemmas 5 a) and 6 b) imply

∥D̃^{k}d_{s}^{k}∥ ≤ γ\sqrt{Δ_{k}},

where γ ≐ δ (β/\sqrt{1 − β} + \sqrt{n}), therefore

|\begin{matrix}√ & \overline{x̃_{i}^{k}} \\ √ & \overline{s̃_{i}^{k}} \end{matrix} (d_{s}^{k})_{i}| ≤ γ\sqrt{Δ_{k}}

hence (using (30), (31), and (1)),

|\frac{(d_{s}^{k})_{i}}{Δ_{k}} | ≤ \frac{γ\sqrt{s̃_{i}^{k}}}{\sqrt{x̃_{i}^{k}}\sqrt{Δ_{k}}} ≤ \frac{2δγ\sqrt{1 + β}}{x_{i}^{∗}} .  (32)

Since a similar argument holds for each i, (29) and (32) imply (27).

Next, we prove (28). Again, consider component i, and assume that x_{i}^{k} → x_{i}^{∗} > 0 and s_{i}^{k} → 0 as k → ∞. It is evident that

s̃_{i}^{k} → 0 as k → ∞

because for any fixed index m, since s_{i}^{k} → 0, there is a smallest index p > m such that

\frac{s_{i}^{p}}{s̃_{i}^{m}} < \frac{1}{δ},

hence s̃_{i}^{m} will be updated no later than iteration p. Consequently s̃_{i}^{k} is updated to s_{i}^{k} infinitely often, and since x̃_{i}^{k} is updated to x_{i}^{k} whenever s̃_{i}^{k} is updated, and since x_{i}^{k} → x_{i}^{∗}, we have

x̃_{i}^{k} → x_{i}^{∗} as k → ∞.  (33)

Expression (32) shows that |(d_{s}^{k})_{i}| = O(Δ_{k}) and (33) implies |x_{i}^{k} − x̃_{i}^{k}| → 0, therefore

\frac{|(x_{i}^{k} − x̃_{i}^{k}) (d_{s}^{k})_{i}|}{Δ_{k}} → 0.  (34)

Note also that

|(s_{i}^{k} − s̃_{i}^{k}) (d_{x}^{k})_{i}| = |((s_{i}^{k} − s̃_{i}^{k}) \begin{matrix}√ & \overline{x̃_{i}^{k}} \\ √ & \overline{s̃_{i}^{k}} \end{matrix} ) (\begin{matrix}√ & \overline{s̃_{i}^{k}} \\ √ & \overline{x̃_{i}^{k}} \end{matrix} ) (d_{x}^{k})_{i}| .

But from Lemma 4 a),

|((s_{i}^{k} − s̃_{i}^{k}) \begin{matrix}√ & \overline{x̃_{i}^{k}} \\ √ & \overline{s̃_{i}^{k}} \end{matrix} )| = O (\sqrt{Δ_{k}}) ,

and since (32) implies that

|(\begin{matrix}√ & \overline{s̃_{i}^{k}} \\ √ & \overline{x̃_{i}^{k}} \end{matrix} )| = O (\sqrt{Δ_{k}}) ,

with (29) we conclude that

\frac{|(s_{i}^{k} − s̃_{i}^{k}) (d_{x}^{k})_{i}|}{Δ_{k}} → 0.  (35)

Since a similar argument holds for each i, expressions (34) and (35) imply (28), so θ_{k} → 0 and the algorithm is q-superlinearly convergent.

Recall from (16) that for each iteration,

D̃^{−1}d_{x} + D̃d_{s} = ṽ.

If we multiply this equation by (X̃S̃)^{−1/2}, we see that

X̃^{−1}d_{x} + S̃^{−1}d_{s} = (X̃S̃)^{−1/2}ṽ.  (36)

Once again, consider component i, and suppose x_{i}^{k} → x_{i}^{∗} > 0 and s_{i}^{k} → 0 as k → ∞. Since θ_{k} → 0 and ((d_{x}^{k})_{i}/x_{i}^{k}) → 0, let M̄ be such that k ≥ M̄ implies

\frac{θ_{k}}{1 − β} + |\frac{(d_{x}^{k})_{i}}{x_{i}^{k}} | < \frac{1}{δ}.  (37)

From the argument above, there exists an index m > M̄ such that s̃_{i}^{m} = s_{i}^{m} and x̃_{i}^{m} = x_{i}^{m}. This implies that, for iteration m, the i-th row of (36) simplifies to

\frac{(d_{x}^{m})_{i}}{x_{i}^{m}} + \frac{(d_{s}^{m})_{i}}{s_{i}^{m}} = \frac{1}{\sqrt{x_{i}^{m}s_{i}^{m}}} ṽ_{i}^{m},

which, after solving for (d_{s}^{m})_{i}/s_{i}^{m}, yields

\frac{(d_{s}^{m})_{i}}{s_{i}^{m}} = \frac{θ_{m}Δ_{m}}{x_{i}^{m}s_{i}^{m}} − 1 − \frac{(d_{x}^{m})_{i}}{x_{i}^{m}} .  (38)

Now,

\frac{s_{i}^{m+1}}{s̃_{i}^{m}} = \frac{s_{i}^{m} + (d_{s}^{m})_{i}}{s_{i}^{m}}
= \frac{θ_{m}Δ_{m}}{x_{i}^{m}s_{i}^{m}} − \frac{(d_{x}^{m})_{i}}{x_{i}^{m}} (by (38))
≤ \frac{θ_{m}}{1 − β} + |\frac{(d_{x}^{m})_{i}}{x_{i}^{m}} | (by (7))
< \frac{1}{δ} (by (37)),

so the partial updating scheme will update s̃_{i}^{m+1} = s_{i}^{m+1} and x̃_{i}^{m+1} = x_{i}^{m+1}. Repeating the argument by replacing m with m + 1 and so on shows that s̃_{i}^{k} = s_{i}^{k} and x̃_{i}^{k} = x_{i}^{k} for all k ≥ m. A similar argument holds for each i, therefore from some point on, complete factorizations will be performed every successive iteration.

Corollary 11 If the LCP is nondegenerate and has a strictly complementary solution (x^{∗}, s^{∗}), the complementarity gap sequence generated by the LCP algorithm converges to zero q-quadratically.

Proof. The fact that the complementarity gap converges to zero is enough to guarantee that (x^{k}, s^{k}) → (x^{∗}, s^{∗}), so we do not need to assume that the iterates converge in this case. The previous theorem shows that from some point on complete factorizations will be performed, so the algorithm becomes in essence the LSPF algorithm, which is known to converge q-quadratically in the nondegenerate case with a strictly complementary optimal solution (see [11]).

6 Concluding Remarks

We have shown that, in the nondegenerate case with a strictly complementary optimal solution, a particular partial updating variant of the LSPF algorithm is q-quadratically convergent and has a computational complexity bound for achieving ε-optimality which is a factor of \sqrt{n} smaller than that of the original algorithm. A similar q-superlinear convergence result is obtained for the degenerate case under the assumption that the iterates converge to a strictly complementary optimal solution.

References

[1] K. M. Anstreicher. A standard form variant, and safeguarded linesearch, for the modified Karmarkar algorithm. Math. Programming, 47:337–351, 1990.

[2] K. M. Anstreicher and R. A. Bosch. Long steps in a O(n^{3}L) algorithm for linear programming. Math. Programming, 54:251–265, 1992.

[3] K. M. Anstreicher and R. A. Bosch. A new infinity-norm path following algorithm for linear programming. SIAM J. Optim., 5(2):236–246, 1995.

[4] J. F. Bonnans and C. C. Gonzaga. Convergence of interior point algorithms for the monotone linear complementarity problem. Technical Report 2074, INRIA, Rocquencourt, France, 1993.

[5] R. A. Bosch. On Mizuno’s rank-one updating algorithm for linear programming. SIAM J. Optim., 3(4):861–867, 1993.

[6] R. A. Bosch and K. M. Anstreicher. On partial updating in a potential reduction linear programming algorithm of Kojima, Mizuno, and Yoshise. Algorithmica, 9:184–197, 1993.

[7] D. den Hertog, C. Roos, and J. P. Vial. A complexity reduction for the long-step path following algorithm for linear programming. SIAM J. Optim., 2:71–87, 1992.

[8] C. C. Gonzaga. An algorithm for solving linear programming problems in O(n^{3}L) operations. In N. Megiddo, editor, Progress in Mathematical Programming: Interior-point and Related Methods, pages 1–28. Springer-Verlag, 1989.

[9] C. C. Gonzaga. The largest step path following algorithm for monotone linear complementarity problems. Report of the Faculty of Technical Mathematics and Informatics 94-07, Delft University of Technology, Delft, the Netherlands, 1994.

[10] N. Karmarkar. A new polynomial-time algorithm for linear programming. Combinatorica, 4:373–395, 1984.

[11] K. A. McShane. Superlinearly convergent O(\sqrt{n}L)-iteration interior point algorithms for linear programming and the monotone linear complementarity problem. SIAM J. Optim., 4(2):247–261, May 1994.

[12] S. Mehrotra. Deferred rank one updates in O(n^{3}L) interior point algorithms. Technical Report 90-11, Department of Industrial Engineering and Management Sciences, Northwestern University, Evanston, IL, 1990.

[13] S. Mizuno. A rank one updating interior algorithm for linear programming. The Arabian J. Sci. Engrg., 15:671–677, 1990.

[14] S. Mizuno. O(n^{ρ}L) iteration O(n^{3}L) potential reduction algorithms for linear programming. Linear Algebra Appl., 152:155–168, 1991.

[15] S. Mizuno and M. J. Todd. An O(n^{3}L) adaptive path following algorithm for a linear complementarity problem. Math. Programming, 52:587–595, 1991.

[16] S. Mizuno, A. Yoshise, and T. Kikuchi. Practical polynomial time algorithms for linear complementarity problems. J. Oper. Res. Soc. of Japan, 32:75–92, 1989.

[17] R. D. C. Monteiro and I. Adler. Interior path following primal-dual algorithms, Part I: Linear programming. Math. Programming, 44:27– 42, 1989.

[18] R. D. C. Monteiro and I. Adler. Interior path following primal-dual algorithms, Part II: Convex quadratic programming. Math. Programming, 44:43–66, 1989.

[19] Y. Ye. An O(n^{3}L) potential reduction algorithm for linear programming. Math. Programming, 50:239–258, 1991.