# An Improved Predictor-Corrector Interior-Point Algorithm for Linear Complementarity Problems withO(nL)-Iteration Complexity

Authors: Debin Fang; Qian Yu
Journal: Journal of Applied Mathematics
 (2011)

---

## Abstract

This paper proposes an improved predictor-corrector interior-point algorithm for the linear complementarity problem (LCP) based on the Mizuno-Todd-Ye algorithm. The modified corrector steps in our algorithm cannot only draw the iteration point back to a narrower neighborhood of the center path but also reduce the duality gap. It implies that the improved algorithm can converge faster than the MTY algorithm. The iteration complexity of the improved algorithm is proved to obtainO(nL) which is similar to the classical Mizuno-Todd-Ye algorithm. Finally, the numerical experiments show that our algorithm improved the performance of the classical MTY algorithm.


---

## Body

## 1. Introduction

Since Karmarkar published the first paper on interior point method [1] in 1984, the interior point methodologies have yielded rich theories and algorithms in the fields of linear programming (LP), quadratic programming (QP), and linear complementarity problems (LCP). Among these interior point methods, predictor-corrector interior-point methods play a special role due to their best polynomial complexity and superlinear convergence.In 1993, Mizuno et al. [2] proposed the classical representative of predictor-corrector method for linear programming. The Mizuno-Todd-Ye (MTY) algorithm has O(nL)-iteration complexity which is the best iteration complexity obtained so far for all the interior-point method [3, 4]. Moreover, Ye et al. [5] proved the duality gap of classical MTY algorithm converges to zero quadratically, which dedicated that MTY algorithm has superlinear convergence. The classical MTY algorithm is the first algorithm for LP has both polynomial complexity and superlinear convergence. So the classical MTY algorithm therefore was considered as the most efficient interior point methods for LP.Afterwards, the MTY algorithm was extended to LCP [6, 7] for its excellent performance. In recent papers, Potra [8, 9] presented several predictor-corrector methods for linear complementarity problems acting in a wide neighborhood of the central path, and Stoer et al. [10] proposed a predictor-corrector algorithm with arbitrarily high order of convergence to degenerate sufficient linear complementarity problems.Subsequently, Potra proposed a generalization of the MTY algorithm for infeasible starting points for monotone LCP [11]. Based on the so-called “fast step-safe step” strategy, Wright proposed an infeasible interior point method with polynomial and quadratic convergence for nondegenerate LCP [12]. And Ai and Zhang presented an O(nL)-iteration primal-dual path-following method for monotone LCP based on wide neighborhoods and large updates [13]. These workings improved the MTY algorithm in various aspects. In this paper, the MTY algorithm for LCP will be improved by modifying the corrector of the algorithm to obtain a larger iteration reduce factor than the original, which can guarantee a faster convergence.The typical iterations of the MTY algorithm operate between two neighborhoods of the central path,ρ(1/4) and ρ(1/2). Before starting the MTY algorithm,     an arbitrary initial point ω∈ρ(1/4) will be given.  And then, the predictor produces a point ω̅∈ρ(1/2) by carefully choosing a steplength along the affine-scaling direction at ω and reduces the primal dual gap by a factor of at least (1-χ̂/n) where χ̂=1/84≈0.5946. In the corrector steps, the corrector produces a point ω̂∈ρ(1/4) by taking a unit steplength along the centering direction at ω̅. It is shown that the corrector put the iteration points back to ρ(1/4) which is the narrower neighborhood of the central path and maintains the same duality gap. Thus predictor-corrector reduces the duality gap μk by a factor of at least (1-χ̂/n) and holds the iteration points in the neighborhood ρ(1/4). It follows that the corresponding iterative procedure has O(nL)-iteration complexity.The modified algorithm in this paper has only one corrector step in neighborhood of the central pathρ(1/2) and one predictor step in a narrow neighborhood ρ(1/4) same as the MTY algorithm. Moreover, the modified iteration direction of corrector step can make a reduction for duality gap. This improvement attained a larger iteration reduce factor and results in a faster convergence than the classical MTY algorithm.Section2 describes the procedure of our predictor-corrector algorithm. Section 3 gives the convergence analysis of our algorithm. It shows that the improved algorithm can generate a sequence {(xk,yk)} in the neighborhood ρ(1/4) from an arbitrary initial point (x0,y0)∈ρ(1/4) and converge to optimal solution (xk,yk) with O(nL)-iteration complexity. Finally, the performances of our algorithm are evaluated by numerical experiments in Section 4.


## 2. Description of the Algorithm

In this paper, LCP is to find a vector pair(x,y)∈Rn×Rn such that y=Mx+h, (x,y)≥(0,0) and xTy=0, where h∈Rn and M is a n×n positive semidefinite matrix.Denote the set of all feasible points of LCP byS={(x,y)∣Mx-y+h=0,x≥0,y≥0} and the solution set of LCP by S*={(x*,y*)∣(x*,y*)∈S,x*Ty*=0}. The relative interior of S, S+={(x,y)∣(x,y)∈S,x>0,y>0} is called the set of strictly feasible points, or the set of interior points. We assume S+≠Φ in this paper. It is known that if S+ is nonempty, then for any parameter τ>0 the nonlinear system Mx-y+h=0,xy=τe has a unique positive solution [14]. The set of all such solutions defines the central path C of LCP. So the assumption is sufficient to establish the existence of the central path and the existence of a solution to LCP.The mean value ofxTy is denoted as μ=xTy/n throughout this paper. Then ρ(β)={(x,y)∈S+∣∥Xy-μe∥≤βμ,μ=xTy/n} is considered as the neighborhood of central path, where β>0, e denotes the vector of ones, X=diag(x) and ∥·∥ express the Euclidean norm.The improved predictor-corrector interior-point algorithm (IPCIP) is as follows.


### 2.1. Main Steps of IPCIP Algorithm

Step 0.
Choose an initial pair of interior point(x0,y0) with (x0,y0)∈ρ(1/4), and set the accuracy parameter ε>0. Let k=0.Step 1.
If(xk)Tyk/n≤ε, then stop.Step 2.
Compute the predictor direction(Δxp,Δyp) by solving the system (2.1),
(2.1)MΔxp-Δyp=0,YkΔxp+XkΔyp=23μke-Xkyk.Step 3.
Set(x̂k,ŷk)=(xk,yk)+Qk(Δxp,Δyp) where Qk=min{1/8n,μk/8∥ΔXpΔyp∥}. (We will show later that ∥X̂kŷk-μ̅ke∥≤(1/2)μ̅k, where μ̅k=(1-Qk)μk holds for every (x̂k,ŷk).)Step 4.
Setr=1-1/42n, and compute the corrector direction (Δxc,Δyc) by solving the system (2.2),
(2.2)MΔxc-Δyc=0,ŶkΔxc+X̂kΔyc=rμ̅ke-X̂kŷk.Step 5.
Set(xk+1,yk+1)=(x̂k,ŷk)+(Δxc,Δyc), update k=k+1, and return to Step 1.


## 2.1. Main Steps of IPCIP Algorithm

Step 0.
Choose an initial pair of interior point(x0,y0) with (x0,y0)∈ρ(1/4), and set the accuracy parameter ε>0. Let k=0.Step 1.
If(xk)Tyk/n≤ε, then stop.Step 2.
Compute the predictor direction(Δxp,Δyp) by solving the system (2.1),
(2.1)MΔxp-Δyp=0,YkΔxp+XkΔyp=23μke-Xkyk.Step 3.
Set(x̂k,ŷk)=(xk,yk)+Qk(Δxp,Δyp) where Qk=min{1/8n,μk/8∥ΔXpΔyp∥}. (We will show later that ∥X̂kŷk-μ̅ke∥≤(1/2)μ̅k, where μ̅k=(1-Qk)μk holds for every (x̂k,ŷk).)Step 4.
Setr=1-1/42n, and compute the corrector direction (Δxc,Δyc) by solving the system (2.2),
(2.2)MΔxc-Δyc=0,ŶkΔxc+X̂kΔyc=rμ̅ke-X̂kŷk.Step 5.
Set(xk+1,yk+1)=(x̂k,ŷk)+(Δxc,Δyc), update k=k+1, and return to Step 1.


## 3. Convergence and Complexity Analysis

Lemma 3.1.
Ifb,c∈Rn and b+c=h,bTc≥0,B=diag(b), then(i)
∥Bc∥≤(2/4)∥h∥2;(ii)
bTc≤(1/2)∥h∥2;(iii)
∥b∥2≤∥h∥2.Proof.
The proof of (i) can be seen in Lemma  1 of [2].
From∥h∥2=∥b+c∥2=∥b∥2+∥c∥2+2bTc and bTc≥0, we can obtain bTc≤(1/2)∥h∥2 and ∥b∥2≤∥h∥2. Thus inequalities (ii) and (iii) hold. This completes the proof.Lemma 3.2.
If the point(x̂k,ŷk) was generated by the algorithm, then one has μ̂k=(x̂k)Tŷk/n≥μ̅k, where μ̅k=(1-Qk)μk.Proof.
From Step2 of our algorithm and M is a n×n positive semidefinite matrix, we have (Δxp)TΔyp=(Δxp)TMΔxp≥0. Hence
(3.1)μ̂k=(x̂k)Tŷkn=(xk+QkΔxp)T(yk+QKΔyp)n=(xk)Tyk+Qk[(xk)TΔyp+(Δxp)Tyk]+(Qk)2(Δxp)TΔypn≥(xk)Tykn+Qkn(n⋅23μk-(xk)Tyk)=μk+23Qkμk-Qkμk≥(1-Qk)μk=μ̅k.
This completes the proof.Lemma 3.3.
If(xk,yk)∈ρ(1/4), then the point (x̂k,ŷk) generated by the algorithm satisfies x̂k>0,ŷk>0, ∥X̂kŷk-μ̅ke∥≤(1/2)μ̅k, and (x̂k,ŷk)∈ρ(1/2).Proof.
Let us define(x(Q),y(Q))=(xk+QΔxp,yk+QΔyp), where
(3.2)0≤Q≤Qk=min⁡{18n,μk8‖ΔXpΔyp‖}.
Note that0≤Q≤Qk=min{1/8n,μk/8∥ΔXpΔyp∥}, then we have (2/3)nQ≤1/12 from Q≤1/8n and Q2∥ΔXpΔyp∥≤(1/8)μk from Q≤μk/8∥ΔXpΔzp∥.
Furthermore∥Xkyk-μke∥≤(1/4)μk holds since (xk,yk)∈ρ(1/4). Therefore we have
(3.3)‖X(Q)y(Q)-(1-Q)μke‖=‖Xkyk+Q(XpΔyp+ΔXpyp)+Q2ΔXpΔyp-(1-Q)μke‖=‖Xkyk+Q(23μke-Xkyk)+Q2ΔXpΔyp-(1-Q)μke‖≤(1-Q)‖Xkyk-μke‖+23nQμk+Q2‖ΔXpΔyp‖≤(1-Q)14μk+112μk+18μk=[14+11-Q(112+18)](1-Q)μk≤(14+521)(1-Q)μk≤12(1-Q)μk.
Soxi(Q)yi(Q)≥(1/2)(1-Q)μk>0(i=1,2,…,n) for any Q(0≤Q≤QK).
Since(x(Q),y(Q))=(xk+QΔxp,yk+QΔyp) is continuous with respect to Q,xk>0 and yk>0, it is easily seen that x̂k=x(Qk)>0, ŷk=y(Qk)>0.
Note thatX(Qk)=X̂k,Y(Qk)=Ŷk,(1-Qk)μk=μ̅k, the above argument implies ∥X̂kŷk-μ̅ke∥≤(1/2)μ̅k.
It’s well known that the inequality∥Xz-(xTz/n)e∥≤∥Xz-ce∥ holds for arbitrary vector x,z and arbitrary constant c.
So we have∥X̂kŷk-μ̂ke∥≤∥X̂kŷk-μ̅ke∥≤(1/2)μ̅k≤(1/2)μ̂ from the proof above and Lemma 3.2, thus (x̂k,ŷk)∈ρ(1/2).
This completes the proof.Remark 3.4.
So Lemma3.3 dedicates that the predictors of our algorithm operate in a wide neighborhood of the central path ρ(1/2).Lemma 3.5.
If∥X̂kŷk-μ̅ke∥≤(1/2)μ̅k then ∥(X̂kŶk)-1/2(rμ̅ke-X̂kŷk)∥2≤2(1/4+(1-r)2n)μ̅kProof.
(3.4)‖(X̂kŶk)-1/2(rμ̅ke-X̂kŷk)‖2≤1min⁡1≤i≤nx̂ikŷik‖rμ̅ke-X̂kŷk‖2≤1(1-1/2)μ̅k‖X̂kŷk-μ̅ke+(1-r)μ̅ke‖2≤2μ̅k(‖X̂kŷk-μ̅ke‖2+‖(1-r)μ̅ke‖2)≤2(14+(1-r)2n)μ̅k.
This completes the proof.Lemma 3.6.
Ifn>2 and the point (xk+1,yk+1) was generated by the algorithm, then μk+1≤(1-Qk)μk always hold.Proof.
DenoteD=(X̂k)1/2(Ŷk)-1/2 then
(3.5)D-1Δxc+DΔyc=(X̂kŶk)-1/2(rμ̅ke-X̂kŷk)(D-1Δxc)T(DΔyc)=(Δxc)TΔyc=(Δxc)TMΔxc≥0.
From Lemmas3.1 and 3.5 we have(3.6)(Δxc)TΔyc≤12‖(X̂kŷk)-1/2(rμ̅ke-X̂kŷk)‖2≤(14+(1-r)2n)μ̅k∴μk+1=(xk+1)Tyk+1n=(x̂k+Δxc)T(ŷk+Δyc)n=(x̂k)Tŷk+(x̂k)TΔyc+(Δxc)Tŷk+(Δxc)TΔycn=nrμ̅k+(Δxc)TΔycn≤rμ̅k+(14n+(1-r)2)μ̅k.
Note thatr=1-1/42n, we have
(3.7)(r+14n+(1-r)2)=1-142n+14n+132n=1-42n-932n.
It is obvious that(r+1/4n+(1-r)2)<1 when n>2, so μk+1≤μ̅k=(1-Qk)μk.
This completes the proof.Remark 3.7.
From Lemmas3.2 and 3.6, we can obtain immediately that μk+1≤μ̅k≤μ̂k. That means the corrector reduces μ̂k to μk+1. Because the corrector in the classical MTY algorithm just maintains the same duality gap, the improvement of the corrector in our algorithm will make a faster reduction of μk than MTY algorithm.
Then the polynomial convergence of the algorithm could be established.Theorem 3.8.
Suppose that the sequence{(xk+1,yk+1)} generated by the algorithm, then one has (xk+1,yk+1)∈ρ(1/4).Proof.
Let us define thatxk+1(α)=x̂k+αΔxc,yk+1(α)=ŷk+αΔyc, and denote
(3.8)μk+1(α)=(xk+1(α))Tyk+1(α)n.
As discussed in Lemma 3.6, we also have (D-1Δxc)T(DΔyc)=(Δxc)TΔyc=(Δxc)TMΔxc≥0.
From Lemmas3.1 and 3.5 and r=1-1/42n, then(3.9)‖ΔXcΔyc‖≤24‖(X̂kŶk)-1/2(rμ̅ke-X̂kŷk)‖2≤242(14+(1-r)2n)μ̅k=14(22+22n(1-r)2)μ̅k≤14(22+18)μ̅k≤732μ̅k,μk+1(α)=(xk+1(α))Tyk+1(α)n=(x̂k+αΔxc)T(ŷk+αΔyc)n=(x̂k)Tŷk+α((x̂k)TΔyc+(Δxc)Tŷk)+α2(Δxc)TΔycn≥(x̂k)Tŷk+α((x̂k)TΔyc+(Δxc)Tŷk)n=(x̂k)Tŷk+α(nrμ̅k-(x̂k)Tŷk)n=μ̂k+α(rμ̅k-μ̂k)=(1-α)μ̂k+αrμ̅k≥(1-(1-r)α)μ̅k.
So we can writeμk+1(α)≥(1-(1-r)α)μ̅k, that is, μ̅k≤μk+1(α)/(1-(1-r)α).
Hence(3.10)Xk+1(α)yk+1(α)-μk+1(α)e=(X̂k+αΔXc)(ŷk+αΔyc)-(x̂k+αΔxc)T(ŷk+αΔyc)ne=X̂kŷk+α(X̂kΔyc+ΔXcŷk)+α2ΔXcΔyc-(x̂k)Tŷk+α((x̂k)TΔyc+(Δxc)Tŷk)+α2(Δxc)TΔycne=X̂kŷk+α(rμ̅ke-X̂kŷk)+α2ΔXcΔyc-(x̂k)Tŷk+α(nrμ̅k-(x̂k)Tŷk)+α2(Δxc)TΔycne=(1-α)(X̂kŷk-(x̂k)Tŷkne)+α2(ΔXcΔyc-(Δxc)TΔycne).
Since the inequality∥Xz-(xTz/n)e∥≤∥Xz-ce∥ holds for arbitrary vector x,z and arbitrary constant c, thus(3.11)‖Xk+1(α)yk+1(α)-μk+1(α)e‖≤(1-α)‖X̂kŷk-(x̂k)Tŷkne‖+α2‖ΔXcΔyc-(Δxc)TΔycne‖≤(1-α)‖X̂kŷk-μ̅ke‖+α2‖ΔXcΔyc‖≤(1-α)12μ̅k+α2732μ̅k≤(1/2)(1-α)+(7/32)α21-(1-r)αμk+1(α).
Let us definef(α)=((1/2)(1-α)+(7/32)α2)/(1-(1-r)α), then f′(α)=(-(1/2)r+(7/16)α-(7/32)α2(1-r))/[1-(1-r)α]2.
Whenn>2, we have r≥7/8, this results in f′(α)≤0 and f(α) decreases monotonously, when 0≤α≤1. Because  f(α)≤f(0)=1/2 holds for every α∈[0,1] thus
(3.12)‖Xk+1(α)yk+1(α)-μk+1(α)e‖≤12μk+1(α).
Sinceμk>ε>0 (otherwise the algorithm will be terminated), it follows that
(3.13)Xik+1(α)yik+1(α)≥12μk+1(α)≥12(1-(1-r)α)μ̅k=12(1-(1-r)α)(1-Qk)μk>0.
We have proved thatx̂k=x(Qk)>0,  ŷk=y(Qk)>0 in Lemma 3.3. So we have
(3.14)xk+1(0)=x̂k>0,yk+1(0)=ŷk>0,whenα=0.
Note thatxk+1(α) and yk+1(α) are continuous with respect to α. It implies that xk+1(α)>0 and yk+1(α)>0 hold for any α∈[0,1].
Letα=1, then we have xk+1>0 and yk+1>0.
Furthermore, whenα=1, f(α)=7/32r thus ∥Xk+1yk+1-μk+1e∥≤(7/32r)μk+1≤(1/4)μk+1.
Consider that we can obtainMxk+1-yk+1+h=0 directly from the steps in our algorithm, so (xk+1,yk+1)∈ρ(1/4) holds for every (xk+1,yk+1) generated by the algorithm.
This completes the proof.Theorem 3.9.
The iteration complexity of the algorithm isO(nL).Proof.
LetD′=(Xk)1/2(Yk)-1/2 then
(3.15)(D′)-1Δxp+D′Δyp=(XkYk)-1/2(23μke-Xkyk),((D′)-1Δxp)T(D′Δyp)=(Δxp)TΔyp=(Δxp)TMΔyp≥0.
From Lemma3.1, we have
(3.16)‖ΔXpΔyp‖=‖(D′)-1ΔXpD′Δyp‖≤24⋅‖(Xkyk)-1/2(23μke-Xkyk)‖2≤24⋅(∑i=1n((23μk)2(xik)-1(yik)-1-43μk+xikyik))2.
From∥Xkyk-μke∥≤(1/4)μk, we have xikyik≥(3/4)μk, then (xik)-1(yik)-1≤4/3μk,
(3.17)‖ΔXpΔyp‖≤24⋅∑i=1n(49⋅43⋅μk-43μk)+nμk=24⋅727nμk=72108nμk.
Henceμk/8∥ΔXpΔyp∥≥27/142n, and we can obtain min{1/8n,μk/8∥ΔXpΔyp∥}≥27/142n.
So it implies thatμk+1≤(1-Qk)μk≤(1-27/142n)μk for each k. This  means that μk will decrease at least by a constant factor of (1-27/142n) at each iteration, which guarantees a O(nL)-iteration complexity for the algorithm.
This completes the proof.Remark 3.10.
Note that the reduce factor in the MTY algorithm is(1-χ̂/n) where χ̂=1/84≈0.5946, but the reduce factor in our algorithm is (1-27/142n)=(1-λ/n) where λ=27/142≈1.168, which is larger than χ̂. It implies that our algorithm will converge faster than the MTY algorithm, although both have the same iteration complexity.


## 4. Examples and Numerical Results

Finally, the numerical experiments were carried out to evaluate the performance and practical efficiency of the algorithm,Test Problem 1
Find a vector pair(x,y)∈R3×R3 such that y=M1x+h1,(x,y)≥(0,0),xTy=0, where h1∈R3 and M1 is a 3×3 positive semidefinite matrix,
(4.1)M1=(213320115),h1=(-10-2).Test Problem 2
Find a vector pair(x,y)∈R5×R5 such that y=M2x+h2,(x,y)≥(0,0),xTy=0, where h2∈R5 and M2 is a 5×5 positive semidefinite matrix,
(4.2)M2=(7002028359003030146180025),h2=(-9.5-36.5-5-14-18.5).Test Problem 3
This example is a general test problem used by Noor et al. [15]. The problem is to find a vector pair (x,y)∈Rn×Rn such that y=M3x+h3,(x,y)≥(0,0),xTy=0, where h3∈Rn and M3 is a n×n positive semidefinite matrix. To test the efficiency of our algorithm by large-scale problem, we set the dimension of test problem 3 as 100, that is, n=100,
(4.3)M3=(4-20⋯014-2⋯0014⋯0⋮⋮⋮⋮⋮000⋯4),h3=(11⋮11).These test problems are solved by our IPCIP algorithm and the classical MTY algorithm. The experiments are run on a PC (2.6 GHz CPU, 2 G DDR RAM) using MATLAB 7. The results including the corresponding numbers of iterations,μk, the approximate solutions, and computing times are shown in Tables 1, 2, and 3.Table 1
Computational results of test problem 1.
MTYIPCIPx0(1, 1, 1)T(1, 1, 1)Ty0(5, 5, 5)T(5, 5, 5)Tε1e-061e-06Iterations1211μk3.76e-075.73e-07xk(0.0, 0.0004, 0.3999)T(0.0, 0.0009, 0.3998)Tyk(0.2002, 0.0009,9.29e-07)T(0.2004, 0.0018,1.81e-07)TRun time (s)0.0782200.049975Table 2
Computational results of test problem 2.
MTYIPCIPx0(1.5, 1.5, 1.5, 1.5, 1.5)T(1.5, 1.5, 1.5, 1.5, 1.5)Ty0(4, 4, 4, 4, 4)T(4, 4, 4, 4, 4)Tε1e-061e-06Iterations109μk2.61e-082.18e-07xk(0.8402,1.3726,0.0347,1.8094,1.6320)T(0.8402,1.3726,0.0347,1.8094,1.6320)Tyk(0, 0, 0, 0, 0)T(0, 0, 0, 0, 0)TRun time (s)0.0713350.037116Table 3
Computational results of test problem 3.
MTYIPCIPx0(1, 1, 1,…, 1, 1)T(1, 1, 1,…, 1, 1)Ty0(3, 4, 4,…, 4, 6)T(3, 4, 4,…, 4, 6)Tε1e-051e-05Iterations87μk1.00e-081.92e-08xk(0, 0, 0,…, 0, 0)T(0, 0, 0,…, 0, 0)Tyk(1, 1, 1,…, 1, 1)T(1, 1, 1,…, 1, 1)TRun time (s)0.2173750.177515


## 5. Conclusion

This paper modified the MTY algorithm for solving monotone LCPs to strengthen the convergence results. Although the iteration complexity of the improved algorithm was proved to beO(nL) which is similar to the classical MTY algorithm, but the reduce factor of duality gap was enhanced to 1.168 and results in a faster convergence than the classical MTY algorithm.The numerical results show that our IPCIP algorithm is more efficient than the MTY algorithm. The number of iterations was decreased, and the computing times could be reduced by nearly 20% to 50%. That indicated IPCIP algorithm has a better performance than the MTY algorithm.



---

*Source: 340192-2011-12-29.xml*