An enumerative method for the solution of linear complementarity problems

122

European Journal of Operational Research 36 (1988) 122-128 North-Holland

Theory and Methodology

An enumerative method for the solution of linear complementarity problems J.J. J U D I C E

Departamento de Matematica, Universidade de Coimbra, Coimbra, Portugal G. M I T R A

Department of Mathematics and Statistics, Brunel University, Uxbridge and Unicorn, Brunel Science Park, United Kingdom

Abstract: In this paper an enumerative method for the solution of the Linear Complementarity Problem (LCP) is presented. This algorithm either finds all the solutions of any general LCP (that is, no assumption is made concerning the class of the matrix), or it establishes that no such solution exists. The method is extended to the Second Linear Complementarity Problem (SLCP), A new problem which has been introduced for the solution of a general quadratic program. Keywords: Economics, engineering, calculus of variations, numerical analysis, optimization

1. Introduction The Linear Complementarity Problem (LCP) has become one of the important research areas in mathematical programming. This problem is stated as [2, 12]

w=q+Mz,

z>~0,

w>~O, wXz=O,

(1)

where M ~ R "x" and w,q,z ~ R n. There are a number of enumerative methods which are designed to find 'one solution of the LCP'. These methods mainly rely on a set of heuristic rules to prune the tree and are described in [11] and [8]. In [1] AI-Khayyal designed a hybrid method which incorporates a reduced gradient approach to control the tree search. Recently Judice and Faustino [8] have investigated the application of these methods to sparse and structured LCPs. Received October 1985; Revised June 1987

We can identify a range of situations where it is necessary to compute all the solutions of an LCP. Finding all the K u h n - T u c k e r points for a nonconvex quadratic programming problem [2] problem, and finding all the equilibrium points for a bimatrix game and some other problems of economic equilibrium [13] are examples of these. Williams has recently shown [16] that finding all the evolutionary stable strategies for a model representing development of species is equivalent to finding all the solutions of an LCP. In [10] we have shown that a number of nonconvex optimization problems can be cast as a minimum linear complementarity problem (MLCP) [6]. The M L C P is an extension which involves minimising a linear objective function subject to restrictions which take the form of an LCP. Hence one way of solving the M L C P would be to compute all the solutions of the LCP. A number of tree search methods have been proposed for finding all the solutions of the LCPs.

0377-2217/88/$3.50 © 1988, Elsevier Science Publishers B.V. (North-Holland)

J.J. Judice, G. Mitra / Enumerative method for the solutions of LCPs

These are due to Garcia and Lemke [4], Murty [14], Jahanshalou and Mitra [7] and A1-Khayyal [1]. In this paper we describe another enumerative method for this problem. In designing the method we have introduced strategies which speed up solution of subproblems and limit the tree search. Our method can also process the more general problem stated as (2)

A w + Bz + Cu = b, w>~O, z>~O, w V z = O ,

and u,'s are either non-negative or unrestricted. If a quadratic program contains some unrestricted variables and (or) some equality constraints then K u h n - T u c k e r conditions represent a problem of the form w=q+ Mz + Nu,

(3)

0 = p + Rz + Su, Z ~0,

W~0,

zTw:0,

and u is unrestricted. In [10] we have called this the second linear complementarity problem (SLCP). Using the well known method of replacing a free variable by two non-negative variables the SLCP can be easily rest.ated as an LCP. From a computational point of view such a transformation, however, is not advisable and this is particularly true for structured problems. The contents of the paper are organized as follows. In Section 2 we provide an overview of our method and in Section 3 we briefly discuss some relevant aspects of the implementation. To illustrate the different computational steps we provide a worked example in Section 4. A summary of test problems and the experience of our computational investigation are given in Section 5. 2. An overview of the enumerative method

The method is based on the principle that all the solutions of an LCP can be found by generating at most ~ - 0 2i nodes of a tree. depth 1

=,~

.,~

Zt 1

depth 2 depth 3 •

depth n - 1

depth n

123

where (i I . . . . . i n } is a permutation of the set of indices {1,2 . . . . . n }. Prior to tree search a feasible solution which satisfies the linear constraints is first obtained by applying the Phase I of the simplex method. Each node of the tree is derived by solving a linear program which consists of minimising zip or w,p ( p = 1 . . . . . n) subject to the linear constraints of the parent problem. This is achieved by a simple modification of the simplex method whereby the chosen branching variable currently in the basis is nominated as the objective function. We call this the 'Minvar' procedure [9]. The 'Minvar' procedure may terminate in two different ways which may lead to two corresponding conclusions concerning the chosen variable zip

(wi). Termination 1. The minimum value of the branching variables is zero. The corresponding branch is generated and the variable is 'starred', that is, zi, = 0 (or wi, = 0). Termination 2. The minimum value of the branching variable is positive. The branch cannot be created and the corresponding part of the tree is terminated.

By using this procedure we guarantee that we are only working with solutions which are feasible for the linear constraints of the LCP. A complementary feasible solution, that is a solution of the LCP, may be reached at any depth of the tree. Since we wish to generate all the solutions of the LCP the search has to be continued until the depth n. This implies that, in case of nondegeneracy, the number of solutions to the LCP is equal to the number of nodes which can be created at depth n. To implement this tree development scheme we need two alternative branching strategies. If the solution (tableau) is complementary then any pair of complementary variables such that neither of them are starred has at least one nonbasic varibale. Therefore the branches for complementary tableaus are of the form:

(D1) °°

star nonbasic variable

apply Minvar

(D2)

124

J.J.

Judice, G. Mitra / Enumerativemethodfor the solutionsof LCPs

If the solution (tableau) is noncomplementary there is at least one pair of complementary nonstarred basic variables (zi, w~) and in order to achieve complementarity each variable is starred in one of the two branches. Hence Minvar is applied to generate both the nodes and the branches are of the form:

the number of starred variables is equal to n go to Step 4. Otherwise go to Step 1. EXIT. If NSOL = 0 the LCP has no solution. Otherwise all the solutions of the LCP have been enumerated and there are NSOL number of solutions.

3. Some aspects of the algorithm implementation

] apply Minvarfor zi / ~

U

apply Minvarfor w~]

(D3)

"o

The enumerative method based on these ideas is summarised below.

Step O. Initiafization. Set NODE = 1 (the node under current investigation), NNODE = 1 (total number of nodes to be investigated), NSOL = 0 (number of solutions of the LCP). Apply the Phase I of the simplex method. If no feasible solution exists go to EXIT. Step 1. Check tableau state. If the tableau is complementary go to Step 2, otherwise go to Step 3. Step 2. Branching for complementary tableaus. Branch by starring a nonbasic variable in any column of the tableau in node (NNODE + 2) and designating its complementary to be starred in node (NNODE + 1). Store the tableau, set NNODE = N N O D E + 2 and N O D E = N N O D E . If the number of starred variables is equal to n go to Step 4. Otherwise reapply Step 2.

Step 3. Branching for noncomplementary tableaus. Consider two basic variables (zip, w~) which constitute a pair of complementary variables. Designate zip w~, to be starred in nodes (NNODE + 1) and (NNODE + 2), store the tableau, set N N O D E = N N O D E + 2 and N O D E = N N O D E and go to Step 6. Step 4. Count and store solutions. Set NSOL = NSOL + 1. Add the solution to the list and go to Step 5. Step 5. Backtrack. If NODE = 1 go to EXIT. Otherwise set NODE = N O D E - 1. If this node (NODE) has already been processed reapply Step 5. Otherwise extract the last tableau stored and go to Step 6. Step 6. Generation of the node. If the chosen variable to be starred has a positive minimum value go to Step 5. Otherwise star this variable. If

In this section we discuss in a summary form four features of the algorithm that we have implemented. These procedures are fully described in our report [9]. (a) For a complementary tableau the nonbasic variable, say zi; to be starred can be chosen such that a basic variable (not complementary to the current one) that is wi, (t 4:p) is forced to have a positive minimum value. As a result the variable zi,, which is currently nonbasic should be starred, which reduces the number of nodes to be developed. (b) V~Tehave designed an algorithm to prepare a noncomplementary tableau for branching. This procedure tries to minimize the pair of chosen complementary variables at the same time. We call this the 'Joint Minvar' procedure. (c) Whenever we reach the depth (n - 1) of the tree it is always possible (by looking at the nonstarred column) to deduce the solutions at depth n. Thus, it is not necessary to generate the tableau at depth n. (d) The second linear complementarity problem can be easily processed by a simple modification in phase 1 which makes basic all the unrestricted variables. Subsequently, the enumerative method is used as before except that these (unrestricted) variables are not candidates to leave the basis.

4. A worked example illustrating the enumerative method In this section an example taken from [7] is solved by the enumerative method. Consider the LCP given by the following tableau:

1 w1= w2= w4 =

-6

z 1

-2 -10

3

-1

-4

w3=

-

2

-20

--

z 2

--

1 -1 2 -3

z 3

--

3 1

-4 -1

-1

z 4

(T1) 2

1

3

125

J.J. Judice, G. Mitra / Enumerative method for the solutions of LCPs

The Phase I of the simplex m e t h o d requires one pivot t r a n s f o r m a t i o n to get an initial feasible solution c o n t a i n e d in the following tableau: w1= z1 = w3 = W4=

1

-- w 2

-- Z 2

-- z 3

-- 2 4

2.8 0.4 2.6 2.0

-0.2 -0.1 0.1 --2.0

1.2 0.1 1.9 --1.0

2.8 -0.1 -0.9 --1.0

-3.8 0.1 1.9 5.0

(T2)

z1=10.99, w 2=97.99, z3=7.99, =205.99, wl = z 2 = w 3 = z 4 = 0 .

Step 5 (Backtrack) follows a n d node 2 is generated by starring the variable w 4. To do this, tableau (T4) is considered and M i n v a r is applied. O n e iteration of M i n v a r is sufficient to make w 4 n o n b a s i c a n d the following c o m p l e m e n t a r y tableau is obtained: 1

This tableau is n o n c o m p l e m e n t a r y . So Step 3 is applied a n d following the Joint M i n v a r procedure a pivot t r a n s f o r m a t i o n with ~32 as the pivot is performed which leads to the following tableau: 1

w1= z~ = z2 = w4 =

-

1.15

w

2

-0.26 - 0.1 0.05 - 1.94

0.26

1.37 3.36

--

w 3

--

-0.63 - 0.05 0.52 0.52

2 3

3.36 - 0.05 - 0.47 1.47

--

1

-

0.39 0.28 1.53 3.87

w2

-0.07 -0.1 0.01 - 2.06

- w3

0.18 -0.06 0.43 0.25

-

wI

-0.29 -0.01 0.14 0.43

-5 0 1 6

(T3)

- za

-1.48 -0.07 0.29 3.81

(T4)

This is a c o m p l e m e n t a r y tableau. Since w x is starred (w 1 = 0) the c o r r e s p o n d i n g c o l u m n is disregarded whereby z 3 attains its m i n i m u m value. As this is positive, the c o m p l e m e n t a r y variable w 3 is starred. By applying Step 2 two following b r a n c h e s are created. g

g 4

i

n

starred / / / ~

W 4

starred

I n node 3 we have a c o m p l e m e n t a r y tableau (T4) with three starred variables ( 2 4 = W 1 = 1'9 3 = 0 ) . As explained in section 3 this leads to the two solutions set out below: Solution

1. z I = 0.28, z 2 3.87,

=

z3 = z1 = z2 = z4 =

1.85 0.36 1.22 1.01

--

w 2

-0.88 -0.15 0.17 -- 0.54

-- W 3

-- w 1

I [ ] ]

] I I /

-- w 4

0.38 0.02 -0.07 0.26

(T5)

g 4

I n (T3) zl attains its m i n i m u m value. Since this value is positive w a has to be starred. F o r this the m i n v a r procedure is applied a n d a pivot transform a t i o n with ~ 3 as the pivot is performed yielding the tableau (T4). z3 = zl = z2 = w4 =

w4

1.53, z 3 = 0.39, w 4 =

W1 = W 2 : W 3 ~ z 4 = O .

Solution 2. ( O b t a i n e d from solution 1 b y a pivot t r a n s f o r m a t i o n with d31 as the pivot.)

We note that it is not necessary to u p d a t e the second a n d the third c o l u m n s since the variables w 3 and w a are starred. At this stage the n u m b e r of starred variables is equal to 3 a n d there is exactly one n o n s t a r r e d n o n b a s i c variable. So node 2 provides two more solutions which are 3. z 1 = 0 . 3 6 , z 2 = 1.22, z 3 = 1.85, z 4 = 1.01, w l = w 2 = w 3 = w 4 = 0 . Solution 4. z 1 = 1 . 4 1 , w 2 = 6 . 9 7 , z 3 = 7 . 9 9 , w 4 = 4.79, w l = z 2 = w 3 = w 4 = 0 .

Solution

After this we backtrack to n o d e 1 a n d the tree search is completed. Whereas the complete enum e r a t i o n consists of 31 nodes, only 3 nodes are e x a m i n e d b y our search procedure. T o c o m p u t e these four solutions of the L C P only 6 pivot t r a n s f o r m a t i o n s were required.

5.

Computational

experience

and

discussion

of

re-

stilts

I n this section a set of test LCPs a n d SLCPs either constructed or taken from k n o w n sources are presented in T a b l e 1 together with experimental results of applying our e n u m e r a t i v e m e t h o d to these problems. T o carry out these experiments the enumerative m e t h o d was p r o g r a m m e d using a contracted tableau simplex m e t h o d with reinversion procedure for the purpose of accuracy. The tree was developed using a last in first out ( L I F O ) scheme a n d tableaus were stored explicitly in the disk backing store. The investigation was carried out using a Honeywell DP68 Computer. I n Table 1 the p r o b l e m n u m b e r is presented in the first column. The source of the respective

J.J. Judice, G. Mitra / Enumerative methodfor the solutions of LCPs

126 Table 1 Problems and results Problem

P1 P2 P3 P4 P5 P6 P7 P8 P9 PI0 Pll P12 P13 P14 P15 P16 P17 PI8 P19 P20 P21 P22 P23 P24 P25

Source

$1 $2 $3 $4 $4 $4 $4 $4 $4 $5 $5 $5 $5 $5 $6 $6 $6 $6 $6 $6 $6 $6 $6 $6 $6

Dimensions

Matrix

V

U

I

E

order

4 4 4 9 6 10 9 10 16 8 8 9 11 12

0 0 0 2 0 2 0 0 0 2 1 0 2 0

3 3 3 6 12 9 15 15 10 7 12 12 9 8

0 0 0 0 2 0 0 0 1 1 2 3 3 0

7 7 7 17 20 21 24 25 27 18 23 24 25 20 20 25 30 33 35 37 40 42 45 48 50

problem is given in the second column by a parameter which is explained at the end of the table. The first fourteen SLCPs and LCPs are the Kuhn-Tucker conditions of different nonconvex quadratic programming problems. The dimensions of these different quadratic programs are set out in the third column under the following headings. V U I E

-

number number number number

of of of of

nonnegative variables, unrestricted variables, inequality constraints, equality constraints.

The order of the matrix of the SLCP or LCP is presented in the fourth column. The number of solutions, pivot steps and nodes and the CPU time in seconds of the application of the enumerative method to the test problems are presented in the remaining columns. List of souces for the test problems Sl, $2, $3 - Hansen and Mathiesen problems 1.1, 1.4 and 1.5 [5];

Solutions

Iterations

Nodes

3 3 11 15 1 4 3 3 3 7 3 37 376 589 7 62 29 44 66 40 117 29 50 33 27

23 17 43 663 204 610 2090 3534 4548 49 314 720 6509 6761 293 482 513 429 1568 1407 2282 1904 2479 1686 1640

11 15 25 239 79 163 479 785 971 15 117 335 2603 2835 69 169 117 101 503 407 711 465 735 605 521

$4 - Hansen and Mathiesen random problems [5]; $5 - Let the quadratic program be of the form Minimise

pVz + ½zXDz,

subjectto

A z ( >~=)b,

z~R',

b ~ R m,

where the variables may be nonegative or unrestricted. Then these test problems are obtained by using coefficient values pi = --1~

-2

d,j=

if

i=j,

1

if [ i - j [ = l ,

0

otherwise,

bi= -1,

aij=

i=l,...,m;

i,j=l

. . . . . l,

- a O,

j=l,...,l,

where a 0 is a random number drawn from a uniform distribution in the interval [0,1], $6 -

These are LCPs whose matrix M ~ and vector q ~ R n are given by

qi = --5Cto,

R "×"

J.J. Judice, G. Mitra / Enumerative method for the solutions of LCPs

i +j-

References

ifi+j-2
i+j-2,

mij =

2-

n,

otherwise,

i,j=l,...,n, w h e r e a 0 is a r a n d o m n u m b e r d r a w n u n i f o r m d i s t r i b u t i o n i n t h e i n t e r v a l [0,1].

from a

Brief discussion of the results To our knowledge there is no extensive computational experience with the other enumerative methods (see Section 1) for finding all solutions of the LCP. For the purpose of comparison we can only use the results reported by Hansen and Mathiesen [5]. As stated earlier we solved three exact problems and few randomly generated problems taken from [5]. For these 3 (small) problems the results are set out in Table 2. For the randomly generated problems, although constructed in the same way, similar comparisons are not as meaningful. An examination of the number of pivot transformations for problems of some dimension shows the same trend between the two methods. If the matrix of an LCP is nonnegative then for any i such that qi > 0 the variable z i can be starred before starting the algorithm. Because of this property only nonpositive vectors q are considered in the LCPs of source 6. Note that Lemke's method [12] cannot be applied to these LCPs since there is no vector d > 0 such that M ~ L ( d ) [3]. The results presented in Table 1 for these LCPs are encouraging and show that the enumerative method does not perform as poorly as expected of a tree search approach with an increase of the dimension of the LCP. In a sense this confirms Kaneko and Hallman's experience of finding one solution for LCPs with non-negative matrices [11].

Table 2 Problem

127

Number of pivot transformations Enumerative method

Hansen and Mathiesen approach

1.1

23

124

1.4

17

85

1.5

43

232

[l] Al-Khayval, F.A., "An implicit enumeration procedure for the general linear complementarity problem", Mathematical Programming Studies 31 (1987) 1-20. [2] Cottle, R.W. and Dantzig, G.B., "Complementary pivot theory of mathematical programming", in: G.B. Dantzig and A.F. Veinott Jr. (eds.), Mathematics of Decision Sciences, American Mathematical Society, Providence, 1968, 115-135. [3] Garcia, G.B., "Some classes of matrices in linear complementarity theory", Mathematical Programming 5(1973) 299-310. [4] Garcia, G.B., and Lemke, C.E., "All solutions to linear complementarity problems by implicit search", paper presented at the 39th National meeting of Operations Research Society of America, 1971. [5] Hansen, J. and Mathiesen, L., "Generating stationary points for a nonconcave quadratic program by Lemke's almost complementarity pivot algorithm", Discussion paper 11/73, Norwegian School of Economics and Business Administration, Bergen, 1973. [6] Ibaraki, T., "Complementary programming", Operations Research 19 (1971) 1523-1529. [7] Jahanshalou, G.R. and Mitra, G., "Linear complementarity problem and a tree search algorithm for its solution", in: A. Prekopa, (ed.), Survey of Mathematical Programming, Budapest, 1979, 35-56. [8] Judice, J.J. and Faustino, A., "An experimental investigation of enumerative methods for the linear complementarity problem", Technical Report, Department of Mathematics, University of Coimbra, 1986, presented at Euro VIII, Lisbon, September 1986. [9] Judice, J.J. and Mitra, G., "An enumerative method for the solution of linear complementarity problems", Technical Report, TR/04/83, Department of Mathematics and Statisics, Brunel University, 1983. [10] Judice, J.J. and Mitra, G., "Reformulation of mathematical programming problems as linear complementary problems and an investigation of their solution methods", to appear in JOTA, 1987. [11] Kaneko, I. and Hallman, W.P., "An enumeration algorithm for a general linear complementarity problem", Technical Report WP 78-11, Department of Industrial Engineering, University of Wisconsin, Madison, 1978. [12] Lemke, C.E., "On complementary pivot theory", in: G.B. Dantzig and A.F. Veinott Jr. (eds.), Mathematics of Decision Sciences, American Mathematical Society, Providence, 1968, 95-114. [13] Lemke, C.E., "A survey of complementarity theory", in: F. Giannessi and R.W. Cottle (eds.), Variational Inequalities and Complementarity Problems: Theory and Applications, Wiley, New York, 1980. [14] Murty, K., "An algorithm for finding all the feasible complementary bases for a linear complementary problem", Technical Report 75-2, Department of Industrial" Engineering, University of Michigan, 1972.

128 [15] Pardalos, P.M. and Rosen, J.B., "Global optimisation approach to the linear complementarity problem", Technical Report, Computer Science Department, University of Minnesota, 1985.

[16] Williams, H.P., "Evolution, games theory and polyhedra' to appear in the Journal of Mathematical Biology.