Division of Research Graduate School of Business Administration The University of Michigan April 1980 NEW METHOD FOR LINEAR PROGRAMMING IS FASTER THAN SIMPLEX FOR A CERTAIN CLASS OF PROBLEMS Working Paper No. 209 Zvi Drezner The University of Michigan Dearborn FOR DISCUSSION PURPOSES ONLY None of this material is to be quoted or reproduced without the express permission of the Division of Research.

New Method for Linear Programming Is Faster than Simplex for a Certain Class of Problems Zvi Drezner School of Management The University of Michigan-Dearborn Abstract We introduce an algorithm which is a combination of Khachian's method and the relaxation method for linear programming. We tested the algorithm on problems of finding a feasible point subject to linear constraints, when we have finite lower and upper bounds for each variable. On a randomly generated set of test problems, our method ran faster than the simplex method in C.P.U. time. The run times for the new method ranged between 4.6 percent and 26 percent of the run times for the simplex method for reasonably sized problems. Introduction In this paper we present the relaxation method for solving a set of linear inequalities [1,5] and relate this method to Khachian's method for linear programming [4]. We show how to perform every iteration efficiently and present two improvements in the basic relaxation method: compound constraints and the nestled ball principle. Computational results are presented. The Basic Relaxation Approach Consider the problem of finding a feasible solution to: n Z a. x. < b. for i=l,...,m. j=l 13 J ' (1) Let: n O.(x) = j=l a..x. - b. JiJ I~ for i=l,...,m. (2)

- 2 - The basic relaxation method is: Step 1: Choose an initial solution x(0); set k=0. Normalize each n 2 constraint by dividing aij and b. by ( Z a.i). ijj 1 13 Step 2: Find the constraint which is most violated, i.e., find (0k) = max ((x(k)) i 1 (k) and let r be the constraint for which 0(k) is obtained. Step 3: If 0(k) <, stop. x(k) is feasible. Step 4: Otherwise, update x(k): (k+l) (k) (k) kl) = xk) - (a. for j=l,...,n. 3 3 rj Note that e (x (k+l) 0. Go to step 2. r It is not clear whether the relaxation method will lead to a feasible solution if there is one, or how to recognize a nonfeasible situation. The next two theorems are concerned with these questions. Theorem 1: If a.. and b. are integers, then the relaxation method determines in a 13 1 finite number of iterations if a feasible solution exists. Proof: The proof is based on the following observation: According to Khachian [4], there exists a radius R such that if a feasible solution exists, it must exist in the ball centered at x(0) with radius R(). Let us assume that we have proven that in the k'th iteration a feasible solution must exist inside a ball centered at (k) with radius R(k). The hyperplane passing through x(k) cuts (k) off more than one half of the ball. In fact, it is distant by (k) from the center of the k'th ball, which is x(k). Therefore, the k+l'th ball, which is centered at x(k+l) has a radius of (R(k) - (k) )i As shown by Khachian [4],

- 3 - (k)>2-L(where L is defined there) if there is no feasible solution (k) -L and if () 0 2- a feasible solution exists. Therefore, if (k)< 2-L, we know that there is a feasible solution and if 0(k)> 2-L for every k, then in a finite number of steps we get R(k < 0, and there is no feasible solution. The following theorem is trivial by the proof of theorem 1: Theorem 2: N fkr2 (0)2 If for some N Z 0(k)> R(0) 2 then there is no feasible solution k=0 to problem (1). The relaxation method is probably not polynomial since the number of iterations is not polynomial. The main difficulty in applying Khachian's method [4] is the great accuracy needed in order to assure the finding of a feasible solution if one exists. In [2] a practical method for a special class of problems has been presented. That is; there are given bounds on the variables j.<xj.<j. We can therefore, replace R( by R(0)2 n R = Z (j-gj) 2/4 j=l with center at x.() = (I+Qj)/2 for j=l,...,n. In addition, we assume a given accuracy of e>0 It is preferable to perform a linear transformation yielding xj = (x-j)/(Pj-j), (k)0) x < x!, R(O) 2(n/4) x. = 0.5. Going Beyond the Constraints One can multiply 0(k) by 1+a for a given -1<a<l. There is no sense in using <0. c=0 yields the greatest decrease in R2 since

-4 - R(k+l)2 (k)22 ( 2 (k)2 We have checked the possibility of usinga >0 and have found that for our test problems U=0.8 gave the best results. A More Efficient Computational Procedure Every iteration in the basic relaxation method is of complexity O(mn) We can perform the computation of each interation in O(m)+O(n) if a set=up procedure of complexity O(m2n) and a space for a matrix of size m(m-l)/2 are permitted. Let: T A.. = a.a. ij a a for i=l,...,m. Note that and i..=l. Let (k) = 0 (x (k)) We keep the values of O9k) and replace the updating formula in step 4 of the algorithm to: xk+l) = (k) - (ka for j=l,...,n (4a) j J ' J (k+l) =0 - 0 (k) for i=l,..m (4b) 1 i lr R(k+l) 2 R(k)2 ( (k)2 (4c) Updating the vector x is of complexity 0(n), updating the vector e(k) is of complexity 0(m), and finding the maximum violated constraint is also of complexity 0(m) (it can be even calculated in the same loop in which 0(k) is updated). Therefore, every interation is simple and fast and is of complexity 0(m)+O(n). Calculating all A.. is of complexity O(m2n). 1J

-5 - The Compound Constraint We can replace two or more constraints by a convex combination of constraints, yielding a "better" cut. This idea was presented by Goldfarb and Todd [3] as a surrogate cut. Since an important feature of the relaxation method is the low complexity of each iteration, we would like to retain this feature, and therefore we restrict ourselves to checking only some pairs of constraints, as will be explained later. Let us consider two constraints: alx < b1 T a2x < b2 with 91(x)=alx-b1 and 02(x)=a2x-b2. For 0 < B < 1 it must be true that: T (Pa1 + (1-p)a2)x < lb1 b+ (l-p)b2 (5) We would like to maximize the residual 0 of the compound constraint (5) Since the constraint (5) must be first normalized, we find that: 0M() = (p01 + (1-p)02)/ 02 Ja1 + (1-p)a2. T T T T Since alal=a2a2=l and ala2=a2al2 then 0 (v) = (Pe1 + (l-i)92)/[p2 + (1-V)2 + 2X12 (1-p)]~. (6) Solving d0(P)/dp = 0 leads to the constraint: T T (aa1 + 2a2)x < abl + 8b2 (7) where: a = (01 1282)/(1- 12) S = (2-X1201e)/(1-X22)

- 6 - This combination is valid only if a>O and >,0. The updating formulas (4) are now: (k+l) (k) - - x j = x. -cia. 3a. f x( )x - C lj - a2j f i i i - 2i R(k+1)2 R(k)2 (02 + 02 t cn e 1 = It can easily be verified that 0(k+) (k+l) 1 2 or j=l,...,n (8a) for i=l,...,m (8b) - 2- 0 2)/(-X,).(8c) 12 102)/(1-2). (8c) = 0. We apply the compound constraint formula by the following strategy. Let 01 be a maximum violated constraint. We check combinations of another constraint combined with the maximum violated one. Since 01 = max ( k) 1 i 1 we have a>0 for every i. Therefore, we check only if 0. - Ail1 00 Note that: 2 2 2 2 (0 + 02- - 2X 10)/(-X12) = 01 + (02- A120 ) / (112 so we look for: =0 m ax (0.i- li01)2/(1-i) * (10) i i 1i,0 0 If 9 exists we choose constraint number 2 as the constraint that maximizes 0 Calculations performed by this scheme retain the complexity 0(m)+O(n) for each iteration, but the computational effort is increased. The Principle of the Nestled Ball In this section we will prove that if after some interations the ball is inside the initial ball, then there is no feasible solution. Let R(0) = R R(k) r, an (k) -(0))T(x(kx() = d2 Let R, R r, and (x -x (x X Theorem 3: If R>r+d then there is no feasible solution to the problem (1)

- 7 - Proof: If a feasible solution exists, it must be in the ball centered at x(k) with radius r. Therefore, there must be a feasible solution in the ball centered at x(0) with radius r+d. Since R>r+d, we could have started the iterations with a ball with radius of r+d as R(). Since the value of R() does not affect the values of O., xetc., we would have passed through exactly the same points and would have reached x(k) in k iterations. Since R2r2 remains unchanged, we would have ended with a smaller r (let it be r', where r'2 = (r+d)2 - (R2-r2) ) This smaller r yields yet a smaller R), namely R()= t'+d(, and so on. We get a sequence of r's; let the sequence be ro, rl,... where: rg = r (11) r 2 = (rk+d)2 - (R2-r2) There are two possibilities, either r2<O for some k, or r2>O for k k every k. If r2<0 for some k, then there is no feasible solution, since we k I have proven that a feasible solution must lie in a nonexisting ball. If rk>0 for every k, then a limit to the sequence rk exists, since rk is monotonically decreasing and bounded by zero. Let this limit be z. For the limit point we must have by equation (11): z2 = (z+d)2 - (R2-r2) z = (R2-r2-d2)/2d Since R>r+d: z > ((r+d)2-r2-d2)/2d = r But z>r is impossible, since r0=r and the sequence is monotonically decreasing. The theorem is proven.

-8 - Theorem 3 provides us with a better stopping criterion in case of infeasible solution. The geometric interpretation of Theorem 3 is quite interesting. There is a growing ball centered at x(0) such that if x(k) enters this ball there is no feasible solution. The radius of this ball is R-r, and R-r = (Rr2)/ (R+r) = BO(k) /(R+r) It can be shown that there exists a ball centered at x(0) inside which there is no feasible solution. The condition of Theorem 3 holds if the radius of that ball is greater than R(O) We have not yet found a "good use" for this result, even though it gives us more information about the set of possible locations of feasible points. More research is yet to be done. Further Suggestions Equality Constraints If there are equality constraints in the problem, we can, of course, change each of them to a pair of inequalities. However, I believe that the following is a better approach. Suppose we have an equality constraint: n j aljxj = b1 There must exist some alj- 0; or else, if bl=0, the constraint can be ignored, and if bI#O, then the system is inconsistent. Choose any a.lj0 (or choose the maximal in absolute value); let it be all. We have: n 1 = (bl - ax )/all We can substitute xl in all other inequalities (and equalities), thus reducing the number of variables by one. We must add two constraints for the bounds of xl, namely, n.1 (b1 - a xi)/all) Pi. j=2

-9 - In substituting x1 for all equalities, we replace each equality in turn by two inequalities but reduce the number of variables. We can also employ the following strategy. Let every constraint be defined as: n b -ei< Z aijx bi (12) j=l j i where e. is a big number if the left constraint is not applicable. Every 1 two inequalities of type (12) are, for practical purposes, one constraint only. The solution procedure is almost unchanged. We still have to calculate only 0(k) but we must check for max ()(, -e.-e(k) i' ' I i instead of max (0k) 3, which requires almost no additional computational effort. We believe that this approach is superior to that of handling an equality as a pair of inequalities, since the presence of an equality gives a feasible region of zero volume but transforming into a lower dimension space probably yields a nonzero volume for the feasible region. Large Problems The basic relaxation method is well suited for working with direct-access secondary storage. We keep in core the vectors x(, k) Pi and e. If we use the transformation (3), we no longer need x() p, and S. in the memory core, so the in-core storage is only of size 2m+n. a.. is stored on a disk by rows, and so is A.. For every iteration we need only one 13 row of a.. and one of A.. Therefore, only two vectors, one of length n and the other of length m, must be read into the core for every iteration. Sparse Matrices If every row of aij is given by a list of all nonzero elements as (j,aij), then updating x(k) is trivial and very fast, and calculating A.. is also very 13convenient. It might be even more economical to calculate a row of convenient. It might be even more economical to calculate a row of A.. whenever 1j

- 10 - it is needed rather than calculating all X's and storing them. Every iteration will require one pass of the entire matrix a.. If all nonzero a.. are ones, a bit presentation may be very efficient. We calculate a vector holding the number of ones in each row vector of a.. 13 instead of normalization. Calculating %.. simply involves taking a "logical and" operation between words and counting the number of ones in the resulting words. If the number of ones in a word is small, we can efficiently count the number of ones by the following observation. If we arithmetically subtract "one" in the rightmost position of the word from a nonzero number and perform a "logical and" between the result and the word, then the rightmost one is wiped out while the rest of the bits in the word remain the same. By this we do the following: Step 1: Set count to 0. Step 2: If word = O, go to step 7. Step 3: Set count=count+l. Step 4: Calculate word'=word-l. Step 5: Perform a "logical and" between word and word' and put result in word. Step 6: Go to Step 2. Step 7: Exit with count as the number of "ones" in word. Computational Comparison We have generated two types of problems, feasible problems and infeasible problems. The coefficients a.. for the feasible problems were g1,1).3 n generated uniformly on the segment (-1,1). We set b. = E a. i/4 1 j=l1

- 11 - so that x=0.25 is a feasible point. The feasible simplex is not necessarily big. In fact, if m>>n, the point x=0.25 is practically the only feasible point. In addition, we assumed 0 4 x. < 1. The constraints for the infeasible problems were identically generated, but the last constraint was replaced by: m-1 amj =- a.. for j=l,...,n i=l1 m-l b = - (Z b. + O.ln((m-1)/3)). i=l The expected value of the term added to b after normalization of the m constraint is equal to 0.1. We use single precision variables on AMDAHL 470/v7 at the University of Michigan, Ann Arbor, Michigan. We use c=10 4(i.e., 8(k) < means feasible solution). All run times exclude input and are expressed in terms of seconds of C.P.U. In Tables 1 and 2, the basic approach is compared with the basic approach with compound constraints. We have used a=0.8. The number of iterations decreases when the compound constraints are used, but run times on feasible problems remain almost the same. We have decided not to include the compound constraints in the new method for further comparisons.

- 12 - Table 1: Run Times for Compound Constraints: Feasible Problems m n Single Constraints. Compound Constraints.Iterations J Time (sec) Iterations Time (sec). S..~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ I 10 10 10 50 50 50 50 100 100 100 100 200 200 200 200 10 20 50 10 50 100 200 10 50 100 200 10 20 100 200 5 5 6 114 47 36 30 86 299 100 60 75 138 10000* 167 j m 0.004 0.004 0.006 0.032 0.080 0.148 0.290 0.083 0.370 0.582 1.115 0.274 0.497 9.061 4.415 4 3 4 71 34 20 16 55 905 70 31 48 79 10000* 117 0.004 0.004 0.006 0.040 0.090 0.160 0.314 0.096 0.848 0.619 1.179 0.313 0.565 16.095 4.693 I * Run terminated due to iteration limit of 10000.

- 13 - Table 2: Run Times for Compound Constraints: Infeasible Problems m n Single Constraints Compound Constraints I__Iterations Time (sec) Iterations Time (sec) 5 10 11 0.005 4 0.005 10 10 7 0.005 2 0.004 10 100 25 0.018 9 0.015 20 20 2 0.009 2 0.010 20 50 7 0.016 6 0.018 20 100 27 0.035 15 0.036 50 50 37 0.079 15 0.084 50 100 74 0.166 33 0.173 100 100 88 0.585 40 0.613 _~ - ~iii.... nn - l

- 14 - In Tables 3 and 4 compare three methods for linear programming: the simplex method, Khachian's method using deep cuts [2], and our new relaxation method. Since there are so many codes for the simplex method, we must define our method exactly. We first tried the MPS code. Run times were surprisingly high (excluding input-output). This is probably because MPS works with lists of coefficients, a method which is not suitable to our dense matrix problems. Double precision is probably used with other sophisticated techniques which are time=consuming; however, in order to give the simplex a "fair chance," I have coded the "good old" simplex method in single precision with e=104 (i.e., if the coefficients of the objective function are less than c, I assumed optimality). The run time of this unsophisticated simplex was only a fraction of the run time on MPS for all the problems that were tested (not all problems in Tables 3 and 4 were tested on MPS). For example, the feasible problem of 50 by 50 was run by MPS in 7.0 seconds'of C.P.U.. thd initialization phase (input, adding slacks and artificials, etc.) took 3.6 seconds which leaves 3.4 C.P.U. seconds for the 87 iterations needed. In my simplex program the same problem was solved in 37 iterations and only 0.343 seconds which is 10 percent of the time of MPS! The results presented in Tables 3 and 4 speak for themselves and need no additional commentary.

a I - 15 - Table 3: Feasible Problems m n Simplex Khachian Relaxation ____ Iter. Time Iter. Time Iter. Time 5 5 5 0.006 8 0.003 6 0.002 10 10 7 0.010 4 0.003 5 0.004 10 20 7 0.014 6 0.009 5 0.004 20 10 21 0.035 61 0.022 21 0.006 20 20 19 0.046 17 0.021 14 0.008 20 30 15 0.049 10 0.026 8 0.010 20 100 9 0.141 36 0.830 14 0.030 30 50 25 0.171 43 0.261 25 0.033 30 80 22 0.280 32 0.496 21 0.047 40 20 38 0.144 109 0.129 28 0.025 40 60 32 0.314 58 0.518 37 0.065 40 80 28 0.411 50 0.788 27 0.082 50 50 37 0.343 119 0.765 47 0.080 50 100 48 1.059 106 2.511 36 0.148

c.. a - 16 - Table 4: Infeasible Problems m n Simplex Khachian Relaxation Iter. Time Iter. Time Iter. Time 5 10 4 0.006 2 0.003 11 0.005 10 10 5 0.007 5 0.004 7 0.005 10 100 10 0.144 12 0.254 25 0.018 20 20 14 0.035 18 0.021 2 0.009 20 50 16 0.095 27 0.153 7 0.016 20 100 21 0.328 33 0.723 27 0.035 50 50 66 0.610 116 0.728 37 0.079 50 100 160 3.591 107 2.454 74 0.166 100 100 216 7.451 385 9.530 88 0.585

- 17 - References 1. Agmon S. "The Relaxation Method for Linear Inequalities." Canadian Journal of Mathematics 6 (1954), pp. 382-392. 2. Drezner Zvi "Improved Formulas for the Russian Method for Linear Programming." Submitted for publication. 3. Goldfarb, D. and Todd, M. "Combining Inequalities and Improving Stability in the Russian Method for Linear Programming" Polynomial Time Algorithm Workshop, New York, February 1980. 4. Khachian L. G. "A Polynomial Algorithm in Linear Programming" (English Translation). Soviet Mathematics Doklady 20 (1979), pp. 191-194. 5. Motzkin, T., and Schoenberg, I. J. "The Relaxation Method for Linear Inequalities." Canadian Journal of Mathematics 6 (1954) pp. 393-404., I