A polynomial-time algorithm for a class of linear complementarity ...
Mathematical Programming 44 (1989) 1-26 North-Holland
A POLYNOMIAL-TIME ALGORITHM FOR LINEAR COMPLEMENTARITY PROBLEMS
1
A CLASS
OF
Masakazu KOJIMA Department of Information Sciences, Tokyo Institute of Technology, Oh-Okayama, Meguro-ku, Tokyo 152, Japan Shinji M I Z U N O a n d Akiko Y O S H I S E Department of Industrial Engineering and Management, Tokyo Institute of Technology, Oh-Okayama, Meguro-ku, Tokyo 152, Japan
Received 27 July 1987 Revised manuscript received 8 February 1988
Given an n x n matrix M and an n-dimensional vector q, the problem of finding n-dimensional vectors x and y satisfying y=Mx+q,
x~>0, y~>0, x~yi=O ( i = 1 , 2 , . . . , n )
is known as a linear complementarity problem. Under the assumption that M is positive semidefinite, this paper presents an algorithm that solves the problem in O(n 3L) arithmetic operations by tracing the path of centers, {(x, y) E S: x~y~ = I.* (i = 1, 2 , . . . , n) for some/~ > 0} of the feasible region S = {(x, y) >~0: y = Mx + q}, where L denotes the size of the input data of the problem. Key words: Linear complementarity problem, polynomial-time algorithm, path of centers, Karmarkar's algorithm.
1. I n t r o d u c t i o n
Let M be a n n × n matrix, a n d q c R ' . The p r o b l e m of finding a n ( x , y ) c R 2n satisfying y--Mx+q,
( x , y ) > ~ O , x~yi=O ( i = l , 2 , . . . , n )
(1.1)
is k n o w n as a linear c o m p l e m e n t a r i t y p r o b l e m ( a b b r e v i a t e d by LCP), w h i c h has various i m p o r t a n t a p p l i c a t i o n s in linear a n d convex q u a d r a t i c p r o g r a m s , b i m a t r i x games a n d some other areas of e n g i n e e r i n g [3, 12, 18, etc.]. Here R" d e n o t e s the n - d i m e n s i o n a l E u c l i d e a n space. Several c o m p u t a t i o n a l m e t h o d s have b e e n d e v e l o p e d for solving LCPs [3, 12, 24, etc.]. These m e t h o d s apply a s e q u e n c e of pivoting o p e r a t i o n s to the system of linear e q u a t i o n s y = M x + q or a certain artificial system of e q u a t i o n s associated with the LCP. I n some worst cases, they require an e x p o n e n t i a l n u m b e r of p i v o t i n g o p e r a t i o n s [1, 4 a n d 16]. It is w e l l - k n o w n that a n L C P with a n arbitrary m a t r i x M is N P - c o m p l e t e [2].
2
M. Kofima et al. / A polynomial-time algorithm for LCPs
In the field of linear programs, there have been developed many algorithms with a polynomially bounded computational complexity [5, 7, 9, 10, 11, 19, 25, 26, etc.]. For convex quadratic programs, K a p o o r and Vaidya [8] and Ye and Tse [27] have recently proposed polynomially bounded algorithms. These algorithms may be roughly classified into three groups: 1. Ellipsoid algorithms which originated from the first polynomially bounded linear programming algorithm by Khachiyan [10] in 1979. 2. Projective rescaling algorithm by K a r m a r k a r [9], its variations and extensions [5, 8, 26, 27, etc.]. 3. Algorithms using the idea of tracing the path of centers of a polytope [7, 11, 19, 25, etc.]. The algorithms in the last group are relatively new and are closely related with the second group (See [6]). Among these algorithms, the ones given by Gonzaga [7] and Vaidya [25] in the last group for solving linear programs have attained the O(n3L) computational complexity in terms of the number of arithmetic operations. The idea on which the algorithms of the last group are based have been studied in more general framework including convex minimization problems and linear complementarity problems by Megiddo [15], Sonnevend [20] and Tanabe [23]. In the previous p a p e r [11], the authors have proposed an algorithm that solves linear programs in O ( n 4 L ) arithmetic operations by using Newton's method as a numerical tool to trace the path of centers simultaneously in the primal and dual feasible regions. In the present paper we modify and extend their algorithm to a class of LCPs with positive semi-definite matrices. The main emphasis will be laid on the theoretical computational complexity of the algorithm. We do not refer to a practically efficient implementation of the algorithm, which should be studied in the future though. Let R~ and R 2 . denote the nonnegative orthant {x c Rn : x ~>0} of R n and the positive orthant { x c Rn: x > 0 } of R n, respectively. We employ the symbol S for the set of all the feasible solutions of the LCP, Sint its interior and Sop the set of all the solutions of the LCP; S = { ( x , y ) ~ R + 2~.. y = M x + q } ,
Sint=S ~ R++ 2n = {(x, y) c R++. 2,~ . y = M x + q } , Scp={(x,y) ~ S: xiy, = 0 ( i = 1 , 2 , . . . , n)}. Throughout the paper, we impose the following assumptions on the LCP: Assumptions. (i) n I> 2. (If n = 1 then the LCP could be solved trivially.) (ii) All the elements of the n x n matrix M and the vector q are integers. (iii) The matrix M is positive semi-definite, i.e., x T M x >~0 for every x c R n. We may further assume without loss of generality that (iv) each row of the matrix M has at least one nonzero element. To see this, assume on the contrary that all the elements in the ith row of the
M. Kojima et al. / A polynomial-time algorithm for LCPs
3
matrix M are zero. I f q~ is negative then the LCP has no solution. Otherwise we can reduce the size of the LCP to be solved by eliminating the variables xi and yi. We define the size L of the LCP (1.1) by n
L= [2 i~l
n+l
Y~ l o g ( [ a 0 ] + l ) + l o g ( n 2 ) J + l ,
(1.2)
j=l
where aii denotes the (i,j)th element of the n × (n + 1) matrix A = [ M q] consisting of the coefficient matrix M and the constant vector q of the system of equations of the LCP (1.1) and [~:J the largest integer not greater than ~:~ R+. The assumption (iv) ensures the inequality n + l o g ( n 2) ~ 0 } = {(x, y) c Sint : XYe =/ze for s o m e / x > 0}. When Sint ~ 0 and Assumption (iii) is satisfied, the system (1.6) has a unique solution for each/x > 0 and the path SoCnof centers forms a smooth curve in the set Sint. We can also characterize solutions of the system (1.6) in terms of the logarithmic barrier function method. These facts have been indicated and studied partially by Megiddo [15]. We will give a complete proof for these facts in Appendix A. Geometrically the path Seen of centers runs through the interior Sin t of the feasible region S to a solution of the LCP which lies in the boundary of S. Starting from a known initial point in a neighborhood of the path Soen, we trace the path Scan until we attain a sufficiently small parameter ~. The idea of this approach has been suggested by Megiddo [15]. Generally, we are not able to trace the path S ~ accurately because it runs nonlinearly through the interior S~nt of the feasible region S. We are forced to stray from the path Sc~ even if a given initial point lies on the path. As a measure for the deviation of each (x, y) ~ S~t from the path S ~ , we employ the quantity min IIH(~,x,y)ll = min I[XYe-l.~el] ~R+
/~ER+ = IIXYe-(x~y/n)ell
.
It should be noted that the first equality follows from the definition (1.7) of the mapping H and (x,y)e S~,t, and the second because the point (xVy/n)e is an orthogonal projection of the point XYe onto the line {~e:/z c R}. Obviously an (X~ y)C Sin t lies on the path S ~ if and only if IIXYe-(xTy/n)ell = 0. We want to control our approximation of the path so that the deviation IIXYe--(xTy/n)ell converges zero at least linearly as the error xTy for the complementarity slackness tends to zero. This leads to the definition of the c~-center neighborhood S ~ (c~) of the path Scan : S ~ (c~) = {(x, y) ~ Si~t: IlXYe - (xTy/n)ell 0, (x, y ) ~ Sin t and (x,~ y) ~ c R++ 2, as inputs and generates the new point (~, y) as an output, which always satisfies the system of equations = Mff + q
for any (x, y) c Si,t and a n y / ~ > 0.
(3.5)
Concerning the modified N e w t o n procedure we have the following result: Theorem 2. Let a , ~ and 6 be constants such that O < a < ~ O . 1 , 6 = a / ( 1 - a ) and 0 /x - 0.828~'/n 1/2
>~[1-6/n'/2-O.82a/n1/2]~
(by (3.7))
~>0.85~" (since 0 < ~ < 0 . 1 , 3 = a / ( 1 - c ~ ) < 0 . 1 2
and n>~2).
On the other hand, we see,
I I ~ e - ( e l l ~ I I ~ ? e - ~ell ~~O. We then see that each coordinate of :~ is bounded by 2C/n 2, i.e., ~ < (2C/nZ)e. For every coordinate of the basic feasible solution (~,)7) can be represented as the ratio A 1 / A 2 of the determinants dl and A 2 of some n × n submatrices of [E - M q] such that 1 ~< ]32[ and O0 and )7o>0, we obtain ~o=0. []
7. Conclusions
We have presented two algorithms, Algorithms 1 and 2 that solve an LCP satisfying Assumptions (i), (ii) and (iii) in O ( n ° S L ) iterations, where L denotes the size of the input data of the LCP (see (1.2)). An essential idea behind the algorithms is "tracing the path of centers of the feasible region by using Newton's method", which has been successfully utilized by several authors to develop polynomial-time algorithms for linear programs. Algorithm 1 requires O(n 3) arithmetic operations per iteration; hence it has the O(n3SL) computational complexity in terms of arithmetic operations. Algorithm 2 is rather complicated mainly because it incorporates the rank-one update procedure to save O(n °s) arithmetic operations on the average per iteration, but it attains O(n3L) computational complexity.
Acknowledgement
The authors are grateful to the referees for helpful comments which have improved the paper. The remarks on bit computational complexity in Section 5 and in Appendix B have been added according to a suggestion by one of the referees.
M. Kojima et al. / A polynomial-time algorithm for LCPs
21
Appendix A. A characterization of centers in terms of the logarithmic barrier function method
Consider the following quadratic programming problem (QP): Minimize
x Xy
subject to
(x, y) c S,
where S = { ( x , y ) : y - M x = q, x~>0, y>~0}. The QP is equivalent to the LCP in the sense that (x, y) is a solution of the LCP if and only if it is a minimal solution of the QP with the objective value 0. We now apply the logarithmic barrier function method to the QP to replace the nonnegativity condition (x, y)~> 0 by additional logarithmic barrier function terms to the objective function. L(/~):
Minimize
xVy - I~ ~ log xj - ~ ~ log yj j~l
subject to
j
1
(x, y) c Si.t.
The following two theorems characterize a point in the path See. of centers, i.e., a solution of the system (1.6) with ~ > 0, in terms of the problem L(~). Theorem A.1. Let ~ > O. I f (x, y) satisfies the system (1.6) of equations then it is a
minimal solution of L(l~ ). Proof. The objective function of the problem L(/~) can be rewritten as
[xjyj - i~ log(xyj )]. j=l
Since each term x j y j - l x log(xjyj) in brackets [. ] attains the minimum under the condition (x~, yj ) > 0 if and only if XJYJ--/x, the desired result follows. [] Theorem A.2. Let I~ > O. Let the LCP satisfy Assumption (iii) in Section 1. Suppose
that (x, y) is a minimal solution of L(ix ). Then it is a solution of the system (1.6). Proof. We first observe that (x,y) satisfies the K a r u s h - K u h n - T u c k e r optimality
condition (see, for example, [13]) with a Lagrangian multiplier vector u c R n : y-~X-le+MXu=O,
x-/xy-le-u=O
and
y-Mx=q.
We shall show that u -- 0. Then the condition above is equivalent to the system (1.6). Multiplying the diagonal matrices X=diag(xl,xz,...,x,)
and
Y=diag(y~,ya,...,y,)
to the first and the second equalities above, respectively, we obtain Xy-p~e+XM Tu=O
and
Xy-/~e-
Yu--O,
M. Kojima et al. / A polynomial-time algorithmfor LCPs
22 which imply
X ( M T + X -l Y ) u = 0. By A s s u m p t i o n (iii) and (x, y ) > 0 , we see that the n x n matrix X ( M T + X -1 Y) is nonsingular, so that we can conclude that the Lagrange multiplier vector u is zero. [] Theorem A.3. Let the LCP satisfy Assumption (iii) and Sin t ¢ •. Then the problem
L(I~ ) has a unique minimal solution for ever), i~ > O. Proof. Let tx > 0 be fixed. The p r o b l e m L(tx) can be rewritten as
Minimize
x T ( M x + q ) - - t z ~ l o g x i - - / z ~ logyj j=l
subject to
j=l
(x, y) e Sint.
Since the objective function to be minimized is a strictly convex function on a n o n e m p t y convex constraint set S~nt, a minimal solution is unique if it exists. In order to see the existence of a solution to this p r o b l e m , we need only prove that the set S~nt(w) of all the feasible solutions of the p r o b l e m L ( # ) with the objective value not greater than ~o is n o n e m p t y , closed and b o u n d e d for some sufficiently large w. The closedness of the set S~,lt(w) can be verified easily. By the assumption Si,t # 0, the set of all the solutions to the L C P is b o u n d e d . (See [14].) Hence so is the set {(x, y) c S: x-r(Mx + q) x V ( M x + q) is convex by A s s u m p t i o n (iii), the set {(x, y) c S: xT ( M x + q) O. This implies that the m i n i m u m of the quadratic function xT (Mx + q) over the subset
S(r) : {(x, y) c S: II(x, Y)II = r} of S grows at least linearly in r for sufficiently large r. H e n c e the m i n i m u m of the objective function of the p r o b l e m L(/x), which differs only the logarithmic barrier terms --/z ~ l o g x j - - / z ~ l o g y i j--I
j=l
from the quadratic function x T ( M x + q ) , over the set S(r) diverges as r tends to the infinity. This ensures that the set Si.t(w) is b o u n d e d for any w. Finally, we see by the a s s u m p t i o n Sint ¢ 0 that the set Sint (~o) is n o n e m p t y for sufficiently large w. []
M. Kojima et aL / A polynomial-time algorithm ,for LCPs
23
Theorem A.4. Let the L C P satisfy Assumptions (iii) in Section 1, and Sint ~ ~0. Let (x(/~),y(/~)) denote the solution o f the system (1.6). Then (x(/~),y(/~)) is C 1 (continuously differentiable) at every i~ > O. Proof. It is easily verified under Assumptions (iii) that the Jacobian matrix of the mapping H with respect to (x, y) is nonsingular and C 1 at every solution (x, y) of the system (1.6) with g > 0. Thus, the desired result follows from the well-known implicit function theorem (see, for example, [17]). []
Appendix B. Computing an exact solution of the LCP
We shall assume that an approximate solution (£,)~) ~ S of the LCP (1.1) satisfying the relation (1.3) and show how we compute an exact solution by using this information. We employ the notation z for the 2n-dimensional variable vector (x, y) and A for the n x 2n matrix [E - M ] , where E denotes the n x n identity matrix. Then the feasible region S of the LCP (1.1) can be represented as the set of solutions of the system of equations Az=q,z=(x,y)c
R+2, .
(B.1)
Lemma B. Suppose ~ c S. L e t / ~ = {k: ~k < 2-L}. Then there is a vertex z* o f S such that z* = 0
(B.2)
for all k ~ I£.
Proof. By Cramer's rule, we know that each element of a basic feasible solution of the system (B.1) can be represented as the ratio ~11/32 of the determinants of some n × n subrnatrix of [A q]. From the definition (1.2) of L, we see the absolute value of the determinant a of any n × n submatrix of [A q] is bounded by 2L/n 2. See, for example, [21]. This implies that, for any basic feasible solution v, v~=0
ifv~ln22 -L
ifvi>0
(B.3)
holds. Since v is a vertex of S if and only if it is a basic feasible solution of (B.1), (B.3) is true for any vertex of S. On the other hand, the point ~ can be represented as the sum of the convex combination of some vertices v 1, v 2, . . . , v p of S and some unbounded direction u of S such that p
g= ~ cjv;+u, j=l
p
~ cj=l, j~l
cj~>0 ( j = 1 , 2 , . . . , p ) .
24
M. Kojima et aL / A polynomial-time algorithm Jot LCPs
See, for example, [22, Theorem 2.12.6]. In view of Caratheodory's theorem [20, Theorem 2.2.12], we may assume that p ~< 1 + n. Hence we can find an index r such that c~/> 1/(1 + n). We shall show that z* - v r satisfies (B.2). Assume on the contrary that v~>0
for s o m e k ~ / £ .
v~n22
-L
Then
because v r is a vertex of S so that (B.3) holds with v = v r. Since all the components of the vectors v 1, v 2, . . . , v p and u are nonnegative, we obtain P
Zk = ~ CjVJ+Uk>~CrVrk>~(1/(l+n))n22-L>2 -L
(since 2~2-L}. As we will see below, we can move from 2 to a point ~ e S in O(n 3) arithmetic operations such that K ( E ) c K(~) and that the set of columns of the matrix A with indices in K ( ~ ) c is linearly independent. We now consider the system of equations Az=q,
zcR+
2n
,
(B.4) zk=0
for e v e r y k e K ( g ) .
Applying L e m m a B, we see that this system of equations has a solution z* = (x*, y* ), and by the assumption (1.3) and K ( ~ ) c K ( 5 ) that (x*, y*) is an exact solution of the LCP satisfying (1.4). On the other hand, since the columns of the matrix A associated with indices in K(ff) c is linearly independent, tile solution z * = (x*, y*) is unique and can be computed in O(n 3) arithmetic operations. Now we shall show how to move from the point ~ e S to a point g ~ S such that K ( f ) c K ( ~ ) and that the set of columns of the matrix A with indices in K(~) ° is linearly independent. Let ff = ~. We consider the polyhedral set P(ff) consisting of the solutions z = (x, y) of the system Az=q,
z ~ R 2n,
(B.5) zk=~k
for e v e r y k e K ( ~ ) .
By applying the Gaussian elimination to the homogeneous system A u = 0,
(B.6)
M. Kojima et al. / A polynomial-time algorithm Jbr LCPs
25
we c o m p u t e a s o l u t i o n u satisfying uj = 0 ( j c K ( ~ ) ) a n d uk > 0 for s o m e k c K ( ~ ) c if it exists, a n d a m a x i m a l i n d e x subset K o f K (~)~ such t h a t the set o f the c o l u m n s o f the m a t r i x A with the indices in K is l i n e a r l y i n d e p e n d e n t . This requires O ( n 3) a r i t h m e t i c o p e r a t i o n s . I f K = K ( ~ ) ° t h e n ~ itself is a d e s i r e d p o i n t in S. O t h e r w i s e we have a s o l u t i o n u o f the system (B.6) satisfying UJ = 0 ( j c K ( ~ ) ) a n d uk > 0 for s o m e k c K ( ~ ) c. In this case we can m o v e f r o m ~ t o w a r d the d i r e c t i o n - u to o b t a i n a p o i n t z c P ( ~ ) such t h a t ]K(z)] < [K(~)] b y a p p l y i n g a ratio test to the s o l u t i o n o f the n o n h o m o g e n e o u s system (B.5) with the s o l u t i o n u o f the h o m o g e n e o u s system (B.6). H e r e ]K] d e n o t e s the n u m b e r o f e l e m e n t s in an i n d e x set K. T h e n , r e p l a c i n g b y z, we p e r f o r m p i v o t o p e r a t i o n s to the h o m o g e n e o u s system (B.6) to g e n e r a t e a n e w s o l u t i o n u such that uj = 0 ( j c K ( ~ ) ) a n d Uk > 0 for s o m e k c K ( ~ ) c if it exists, a n d r e p e a t the same p r o c e d u r e until we find a p o i n t 2 c S satisfying the d e s i r e d p r o p e r t y . This iteration t e r m i n a t e s in at m o s t n steps since 0
A POLYNOMIAL-TIME ALGORITHM FOR LINEAR COMPLEMENTARITY PROBLEMS
1
A CLASS
OF
Masakazu KOJIMA Department of Information Sciences, Tokyo Institute of Technology, Oh-Okayama, Meguro-ku, Tokyo 152, Japan Shinji M I Z U N O a n d Akiko Y O S H I S E Department of Industrial Engineering and Management, Tokyo Institute of Technology, Oh-Okayama, Meguro-ku, Tokyo 152, Japan
Received 27 July 1987 Revised manuscript received 8 February 1988
Given an n x n matrix M and an n-dimensional vector q, the problem of finding n-dimensional vectors x and y satisfying y=Mx+q,
x~>0, y~>0, x~yi=O ( i = 1 , 2 , . . . , n )
is known as a linear complementarity problem. Under the assumption that M is positive semidefinite, this paper presents an algorithm that solves the problem in O(n 3L) arithmetic operations by tracing the path of centers, {(x, y) E S: x~y~ = I.* (i = 1, 2 , . . . , n) for some/~ > 0} of the feasible region S = {(x, y) >~0: y = Mx + q}, where L denotes the size of the input data of the problem. Key words: Linear complementarity problem, polynomial-time algorithm, path of centers, Karmarkar's algorithm.
1. I n t r o d u c t i o n
Let M be a n n × n matrix, a n d q c R ' . The p r o b l e m of finding a n ( x , y ) c R 2n satisfying y--Mx+q,
( x , y ) > ~ O , x~yi=O ( i = l , 2 , . . . , n )
(1.1)
is k n o w n as a linear c o m p l e m e n t a r i t y p r o b l e m ( a b b r e v i a t e d by LCP), w h i c h has various i m p o r t a n t a p p l i c a t i o n s in linear a n d convex q u a d r a t i c p r o g r a m s , b i m a t r i x games a n d some other areas of e n g i n e e r i n g [3, 12, 18, etc.]. Here R" d e n o t e s the n - d i m e n s i o n a l E u c l i d e a n space. Several c o m p u t a t i o n a l m e t h o d s have b e e n d e v e l o p e d for solving LCPs [3, 12, 24, etc.]. These m e t h o d s apply a s e q u e n c e of pivoting o p e r a t i o n s to the system of linear e q u a t i o n s y = M x + q or a certain artificial system of e q u a t i o n s associated with the LCP. I n some worst cases, they require an e x p o n e n t i a l n u m b e r of p i v o t i n g o p e r a t i o n s [1, 4 a n d 16]. It is w e l l - k n o w n that a n L C P with a n arbitrary m a t r i x M is N P - c o m p l e t e [2].
2
M. Kofima et al. / A polynomial-time algorithm for LCPs
In the field of linear programs, there have been developed many algorithms with a polynomially bounded computational complexity [5, 7, 9, 10, 11, 19, 25, 26, etc.]. For convex quadratic programs, K a p o o r and Vaidya [8] and Ye and Tse [27] have recently proposed polynomially bounded algorithms. These algorithms may be roughly classified into three groups: 1. Ellipsoid algorithms which originated from the first polynomially bounded linear programming algorithm by Khachiyan [10] in 1979. 2. Projective rescaling algorithm by K a r m a r k a r [9], its variations and extensions [5, 8, 26, 27, etc.]. 3. Algorithms using the idea of tracing the path of centers of a polytope [7, 11, 19, 25, etc.]. The algorithms in the last group are relatively new and are closely related with the second group (See [6]). Among these algorithms, the ones given by Gonzaga [7] and Vaidya [25] in the last group for solving linear programs have attained the O(n3L) computational complexity in terms of the number of arithmetic operations. The idea on which the algorithms of the last group are based have been studied in more general framework including convex minimization problems and linear complementarity problems by Megiddo [15], Sonnevend [20] and Tanabe [23]. In the previous p a p e r [11], the authors have proposed an algorithm that solves linear programs in O ( n 4 L ) arithmetic operations by using Newton's method as a numerical tool to trace the path of centers simultaneously in the primal and dual feasible regions. In the present paper we modify and extend their algorithm to a class of LCPs with positive semi-definite matrices. The main emphasis will be laid on the theoretical computational complexity of the algorithm. We do not refer to a practically efficient implementation of the algorithm, which should be studied in the future though. Let R~ and R 2 . denote the nonnegative orthant {x c Rn : x ~>0} of R n and the positive orthant { x c Rn: x > 0 } of R n, respectively. We employ the symbol S for the set of all the feasible solutions of the LCP, Sint its interior and Sop the set of all the solutions of the LCP; S = { ( x , y ) ~ R + 2~.. y = M x + q } ,
Sint=S ~ R++ 2n = {(x, y) c R++. 2,~ . y = M x + q } , Scp={(x,y) ~ S: xiy, = 0 ( i = 1 , 2 , . . . , n)}. Throughout the paper, we impose the following assumptions on the LCP: Assumptions. (i) n I> 2. (If n = 1 then the LCP could be solved trivially.) (ii) All the elements of the n x n matrix M and the vector q are integers. (iii) The matrix M is positive semi-definite, i.e., x T M x >~0 for every x c R n. We may further assume without loss of generality that (iv) each row of the matrix M has at least one nonzero element. To see this, assume on the contrary that all the elements in the ith row of the
M. Kojima et al. / A polynomial-time algorithm for LCPs
3
matrix M are zero. I f q~ is negative then the LCP has no solution. Otherwise we can reduce the size of the LCP to be solved by eliminating the variables xi and yi. We define the size L of the LCP (1.1) by n
L= [2 i~l
n+l
Y~ l o g ( [ a 0 ] + l ) + l o g ( n 2 ) J + l ,
(1.2)
j=l
where aii denotes the (i,j)th element of the n × (n + 1) matrix A = [ M q] consisting of the coefficient matrix M and the constant vector q of the system of equations of the LCP (1.1) and [~:J the largest integer not greater than ~:~ R+. The assumption (iv) ensures the inequality n + l o g ( n 2) ~ 0 } = {(x, y) c Sint : XYe =/ze for s o m e / x > 0}. When Sint ~ 0 and Assumption (iii) is satisfied, the system (1.6) has a unique solution for each/x > 0 and the path SoCnof centers forms a smooth curve in the set Sint. We can also characterize solutions of the system (1.6) in terms of the logarithmic barrier function method. These facts have been indicated and studied partially by Megiddo [15]. We will give a complete proof for these facts in Appendix A. Geometrically the path Seen of centers runs through the interior Sin t of the feasible region S to a solution of the LCP which lies in the boundary of S. Starting from a known initial point in a neighborhood of the path Soen, we trace the path Scan until we attain a sufficiently small parameter ~. The idea of this approach has been suggested by Megiddo [15]. Generally, we are not able to trace the path S ~ accurately because it runs nonlinearly through the interior S~nt of the feasible region S. We are forced to stray from the path Sc~ even if a given initial point lies on the path. As a measure for the deviation of each (x, y) ~ S~t from the path S ~ , we employ the quantity min IIH(~,x,y)ll = min I[XYe-l.~el] ~R+
/~ER+ = IIXYe-(x~y/n)ell
.
It should be noted that the first equality follows from the definition (1.7) of the mapping H and (x,y)e S~,t, and the second because the point (xVy/n)e is an orthogonal projection of the point XYe onto the line {~e:/z c R}. Obviously an (X~ y)C Sin t lies on the path S ~ if and only if IIXYe-(xTy/n)ell = 0. We want to control our approximation of the path so that the deviation IIXYe--(xTy/n)ell converges zero at least linearly as the error xTy for the complementarity slackness tends to zero. This leads to the definition of the c~-center neighborhood S ~ (c~) of the path Scan : S ~ (c~) = {(x, y) ~ Si~t: IlXYe - (xTy/n)ell 0, (x, y ) ~ Sin t and (x,~ y) ~ c R++ 2, as inputs and generates the new point (~, y) as an output, which always satisfies the system of equations = Mff + q
for any (x, y) c Si,t and a n y / ~ > 0.
(3.5)
Concerning the modified N e w t o n procedure we have the following result: Theorem 2. Let a , ~ and 6 be constants such that O < a < ~ O . 1 , 6 = a / ( 1 - a ) and 0 /x - 0.828~'/n 1/2
>~[1-6/n'/2-O.82a/n1/2]~
(by (3.7))
~>0.85~" (since 0 < ~ < 0 . 1 , 3 = a / ( 1 - c ~ ) < 0 . 1 2
and n>~2).
On the other hand, we see,
I I ~ e - ( e l l ~ I I ~ ? e - ~ell ~~O. We then see that each coordinate of :~ is bounded by 2C/n 2, i.e., ~ < (2C/nZ)e. For every coordinate of the basic feasible solution (~,)7) can be represented as the ratio A 1 / A 2 of the determinants dl and A 2 of some n × n submatrices of [E - M q] such that 1 ~< ]32[ and O0 and )7o>0, we obtain ~o=0. []
7. Conclusions
We have presented two algorithms, Algorithms 1 and 2 that solve an LCP satisfying Assumptions (i), (ii) and (iii) in O ( n ° S L ) iterations, where L denotes the size of the input data of the LCP (see (1.2)). An essential idea behind the algorithms is "tracing the path of centers of the feasible region by using Newton's method", which has been successfully utilized by several authors to develop polynomial-time algorithms for linear programs. Algorithm 1 requires O(n 3) arithmetic operations per iteration; hence it has the O(n3SL) computational complexity in terms of arithmetic operations. Algorithm 2 is rather complicated mainly because it incorporates the rank-one update procedure to save O(n °s) arithmetic operations on the average per iteration, but it attains O(n3L) computational complexity.
Acknowledgement
The authors are grateful to the referees for helpful comments which have improved the paper. The remarks on bit computational complexity in Section 5 and in Appendix B have been added according to a suggestion by one of the referees.
M. Kojima et al. / A polynomial-time algorithm for LCPs
21
Appendix A. A characterization of centers in terms of the logarithmic barrier function method
Consider the following quadratic programming problem (QP): Minimize
x Xy
subject to
(x, y) c S,
where S = { ( x , y ) : y - M x = q, x~>0, y>~0}. The QP is equivalent to the LCP in the sense that (x, y) is a solution of the LCP if and only if it is a minimal solution of the QP with the objective value 0. We now apply the logarithmic barrier function method to the QP to replace the nonnegativity condition (x, y)~> 0 by additional logarithmic barrier function terms to the objective function. L(/~):
Minimize
xVy - I~ ~ log xj - ~ ~ log yj j~l
subject to
j
1
(x, y) c Si.t.
The following two theorems characterize a point in the path See. of centers, i.e., a solution of the system (1.6) with ~ > 0, in terms of the problem L(~). Theorem A.1. Let ~ > O. I f (x, y) satisfies the system (1.6) of equations then it is a
minimal solution of L(l~ ). Proof. The objective function of the problem L(/~) can be rewritten as
[xjyj - i~ log(xyj )]. j=l
Since each term x j y j - l x log(xjyj) in brackets [. ] attains the minimum under the condition (x~, yj ) > 0 if and only if XJYJ--/x, the desired result follows. [] Theorem A.2. Let I~ > O. Let the LCP satisfy Assumption (iii) in Section 1. Suppose
that (x, y) is a minimal solution of L(ix ). Then it is a solution of the system (1.6). Proof. We first observe that (x,y) satisfies the K a r u s h - K u h n - T u c k e r optimality
condition (see, for example, [13]) with a Lagrangian multiplier vector u c R n : y-~X-le+MXu=O,
x-/xy-le-u=O
and
y-Mx=q.
We shall show that u -- 0. Then the condition above is equivalent to the system (1.6). Multiplying the diagonal matrices X=diag(xl,xz,...,x,)
and
Y=diag(y~,ya,...,y,)
to the first and the second equalities above, respectively, we obtain Xy-p~e+XM Tu=O
and
Xy-/~e-
Yu--O,
M. Kojima et al. / A polynomial-time algorithmfor LCPs
22 which imply
X ( M T + X -l Y ) u = 0. By A s s u m p t i o n (iii) and (x, y ) > 0 , we see that the n x n matrix X ( M T + X -1 Y) is nonsingular, so that we can conclude that the Lagrange multiplier vector u is zero. [] Theorem A.3. Let the LCP satisfy Assumption (iii) and Sin t ¢ •. Then the problem
L(I~ ) has a unique minimal solution for ever), i~ > O. Proof. Let tx > 0 be fixed. The p r o b l e m L(tx) can be rewritten as
Minimize
x T ( M x + q ) - - t z ~ l o g x i - - / z ~ logyj j=l
subject to
j=l
(x, y) e Sint.
Since the objective function to be minimized is a strictly convex function on a n o n e m p t y convex constraint set S~nt, a minimal solution is unique if it exists. In order to see the existence of a solution to this p r o b l e m , we need only prove that the set S~nt(w) of all the feasible solutions of the p r o b l e m L ( # ) with the objective value not greater than ~o is n o n e m p t y , closed and b o u n d e d for some sufficiently large w. The closedness of the set S~,lt(w) can be verified easily. By the assumption Si,t # 0, the set of all the solutions to the L C P is b o u n d e d . (See [14].) Hence so is the set {(x, y) c S: x-r(Mx + q) x V ( M x + q) is convex by A s s u m p t i o n (iii), the set {(x, y) c S: xT ( M x + q) O. This implies that the m i n i m u m of the quadratic function xT (Mx + q) over the subset
S(r) : {(x, y) c S: II(x, Y)II = r} of S grows at least linearly in r for sufficiently large r. H e n c e the m i n i m u m of the objective function of the p r o b l e m L(/x), which differs only the logarithmic barrier terms --/z ~ l o g x j - - / z ~ l o g y i j--I
j=l
from the quadratic function x T ( M x + q ) , over the set S(r) diverges as r tends to the infinity. This ensures that the set Si.t(w) is b o u n d e d for any w. Finally, we see by the a s s u m p t i o n Sint ¢ 0 that the set Sint (~o) is n o n e m p t y for sufficiently large w. []
M. Kojima et aL / A polynomial-time algorithm ,for LCPs
23
Theorem A.4. Let the L C P satisfy Assumptions (iii) in Section 1, and Sint ~ ~0. Let (x(/~),y(/~)) denote the solution o f the system (1.6). Then (x(/~),y(/~)) is C 1 (continuously differentiable) at every i~ > O. Proof. It is easily verified under Assumptions (iii) that the Jacobian matrix of the mapping H with respect to (x, y) is nonsingular and C 1 at every solution (x, y) of the system (1.6) with g > 0. Thus, the desired result follows from the well-known implicit function theorem (see, for example, [17]). []
Appendix B. Computing an exact solution of the LCP
We shall assume that an approximate solution (£,)~) ~ S of the LCP (1.1) satisfying the relation (1.3) and show how we compute an exact solution by using this information. We employ the notation z for the 2n-dimensional variable vector (x, y) and A for the n x 2n matrix [E - M ] , where E denotes the n x n identity matrix. Then the feasible region S of the LCP (1.1) can be represented as the set of solutions of the system of equations Az=q,z=(x,y)c
R+2, .
(B.1)
Lemma B. Suppose ~ c S. L e t / ~ = {k: ~k < 2-L}. Then there is a vertex z* o f S such that z* = 0
(B.2)
for all k ~ I£.
Proof. By Cramer's rule, we know that each element of a basic feasible solution of the system (B.1) can be represented as the ratio ~11/32 of the determinants of some n × n subrnatrix of [A q]. From the definition (1.2) of L, we see the absolute value of the determinant a of any n × n submatrix of [A q] is bounded by 2L/n 2. See, for example, [21]. This implies that, for any basic feasible solution v, v~=0
ifv~ln22 -L
ifvi>0
(B.3)
holds. Since v is a vertex of S if and only if it is a basic feasible solution of (B.1), (B.3) is true for any vertex of S. On the other hand, the point ~ can be represented as the sum of the convex combination of some vertices v 1, v 2, . . . , v p of S and some unbounded direction u of S such that p
g= ~ cjv;+u, j=l
p
~ cj=l, j~l
cj~>0 ( j = 1 , 2 , . . . , p ) .
24
M. Kojima et aL / A polynomial-time algorithm Jot LCPs
See, for example, [22, Theorem 2.12.6]. In view of Caratheodory's theorem [20, Theorem 2.2.12], we may assume that p ~< 1 + n. Hence we can find an index r such that c~/> 1/(1 + n). We shall show that z* - v r satisfies (B.2). Assume on the contrary that v~>0
for s o m e k ~ / £ .
v~n22
-L
Then
because v r is a vertex of S so that (B.3) holds with v = v r. Since all the components of the vectors v 1, v 2, . . . , v p and u are nonnegative, we obtain P
Zk = ~ CjVJ+Uk>~CrVrk>~(1/(l+n))n22-L>2 -L
(since 2~2-L}. As we will see below, we can move from 2 to a point ~ e S in O(n 3) arithmetic operations such that K ( E ) c K(~) and that the set of columns of the matrix A with indices in K ( ~ ) c is linearly independent. We now consider the system of equations Az=q,
zcR+
2n
,
(B.4) zk=0
for e v e r y k e K ( g ) .
Applying L e m m a B, we see that this system of equations has a solution z* = (x*, y* ), and by the assumption (1.3) and K ( ~ ) c K ( 5 ) that (x*, y*) is an exact solution of the LCP satisfying (1.4). On the other hand, since the columns of the matrix A associated with indices in K(ff) c is linearly independent, tile solution z * = (x*, y*) is unique and can be computed in O(n 3) arithmetic operations. Now we shall show how to move from the point ~ e S to a point g ~ S such that K ( f ) c K ( ~ ) and that the set of columns of the matrix A with indices in K(~) ° is linearly independent. Let ff = ~. We consider the polyhedral set P(ff) consisting of the solutions z = (x, y) of the system Az=q,
z ~ R 2n,
(B.5) zk=~k
for e v e r y k e K ( ~ ) .
By applying the Gaussian elimination to the homogeneous system A u = 0,
(B.6)
M. Kojima et al. / A polynomial-time algorithm Jbr LCPs
25
we c o m p u t e a s o l u t i o n u satisfying uj = 0 ( j c K ( ~ ) ) a n d uk > 0 for s o m e k c K ( ~ ) c if it exists, a n d a m a x i m a l i n d e x subset K o f K (~)~ such t h a t the set o f the c o l u m n s o f the m a t r i x A with the indices in K is l i n e a r l y i n d e p e n d e n t . This requires O ( n 3) a r i t h m e t i c o p e r a t i o n s . I f K = K ( ~ ) ° t h e n ~ itself is a d e s i r e d p o i n t in S. O t h e r w i s e we have a s o l u t i o n u o f the system (B.6) satisfying UJ = 0 ( j c K ( ~ ) ) a n d uk > 0 for s o m e k c K ( ~ ) c. In this case we can m o v e f r o m ~ t o w a r d the d i r e c t i o n - u to o b t a i n a p o i n t z c P ( ~ ) such t h a t ]K(z)] < [K(~)] b y a p p l y i n g a ratio test to the s o l u t i o n o f the n o n h o m o g e n e o u s system (B.5) with the s o l u t i o n u o f the h o m o g e n e o u s system (B.6). H e r e ]K] d e n o t e s the n u m b e r o f e l e m e n t s in an i n d e x set K. T h e n , r e p l a c i n g b y z, we p e r f o r m p i v o t o p e r a t i o n s to the h o m o g e n e o u s system (B.6) to g e n e r a t e a n e w s o l u t i o n u such that uj = 0 ( j c K ( ~ ) ) a n d Uk > 0 for s o m e k c K ( ~ ) c if it exists, a n d r e p e a t the same p r o c e d u r e until we find a p o i n t 2 c S satisfying the d e s i r e d p r o p e r t y . This iteration t e r m i n a t e s in at m o s t n steps since 0
