The Gravitational Method for Linear Programming Katta G. Murty * Department of Industrial & Operations Engineering The University of Michigan Ann Arbor, Michigan 48109-2117, USA Technical Report #86-19 June 1986 *Partially supported by NSF Grant ECS-8521183. Part of this work was completed while visiting Dipartimento di Elettronica of Politechnico di Milano, Milano, Italy on a NATO Grant RG 85-0240; and the Indian Statistical Institute, Delhi, India on the Travel Grant NSF (SFC) - INT - 8521713.

The Gravitational Method for Linear Programming Katta G. Murty ABSTRACT: We discuss a new interior point algorithm for solving linear programs. Geometrically, the method tracks the locus of the center of a drop in the interior of the set of feasible solutions, as it falls under the influence of a powerful gravitational force pulling everything down in the direction of the negative gradient of the objective function. KEY WORDS: Linear programming, interior point, gravitational direction, orthogonal projection.

INTRODUCTION This is a revised version of the algorithm that appeared in the working paper [1]. We consider the linear program (LP) in the form minimize z(x) = cx (1) subject to Ax > b where A is a matrix of order m x n. Sign restrictions on the variables and any other lower or upper bound conditions on the variables, if any, are all included in the above system of constraints. Clearly every LP can be put in this form by well known simple transformations discussed in most LP textbooks (for example, see [2]). If D is any matrix, we denote its ith row by Di., and its jth column by D j. If F is any set, |FI denotes its cardinality. For a real number a, |a| denotes its absolute value. For any vector y, Ilyll denotes its Euclidean norm. NOTE: In practical applications, it usually turns out that the LP model for a practical problem is in standard form min px subject to BX = d (2) X a 0 The dual of this model is directly in form (1) and the gravitational method can be applied to solve the dual of (2) directly. As it will be shown later on, when the gravitational method is applied on the dual of (2), at termination, it will produce an optimum solution for (2), if one exists. ASSUMPTIONS Let K denote the set of feasible solutions of (1). We assume that K $ 0, and that K has a nonempty interior in Rn, and that an initial interior feasible solution x~ (this is a point x~ satisfying AxO > b) of (1) is available. 1

If these assumptions are not satisfied, introduce an artificial variable Xn+1 and modify the problem as follows minimize cx + vxn+1 (3) subject to Ax + exn+1 b where e = (1,..., 1)T Rm and v is a large positive number. For any 2 Rn, let Xn+1 > max{|min{0, Ai R - bi}l: i = 1 to m}, then (, Rn+1) satisfies the constraints in (3) as strict inequalities. Thus the modified problem (3) satisfies all the assumptions made in the above paragraph. We also assume that c + O, as otherwise x~ is optimal to (1), and we can terminate. THE GRAVITATIONAL METHOD The Euclidean distance of x~ from the hyperplane {x: Ai x = bi} is (Ai x~ - bi)/l IAi. I The gravitational approach for solving (1) is the following. Assume that the boundary of K is an impermeable layer separating the inside of K from the outside. Introduce a powerful gravitational force inside K pulling everything down in the direction -cT. Choose e < min{(Ai x0 - bi)/|lAi..: i = 1 to m}. Release a small spherical n-dimensional drop of mercury of diameter 2C with its center at the initial interior feasible solution xO K. The drop will fall under the influence of gravity. During its fall, the drop may touch the boundary, but the center of the drop will always be in the interior of K at a distance > e from the nearest point to it on the boundary. Whenever the drop touches a face of K, it will change direction and will continue to move, if possible, in the gravitational direction that keeps it within K. If the objective function is unbounded below in (1), after changing direction a finite number of times, the drop will continue to fall forever along a half-line in K along which the objective function diverges to -o. If z(x) is bounded below on 2

K, after changing direction a finite number of times, the drop will come to a halt. The algorithm tracks the path of the center of the drop as it falls in free fall under the influence of gravity. Let J denote this path of the center of this drop in its fall. THE GRAVITATIONAL DIRECTION AT AN INTERIOR POINT x6 K Suppose a drop of radius e, with its center at x is inside K. So (Ai.x - bi)/I IAi. | I e, i = 1 to m. (4) At every point x on the locus of the center of the drop in the gravitational method, (4) will always be satisfied. Given a point x on, define J(x) {i: (Ai.x - bi)/|IAi. || e (5) The hyperplane {x: Ai x = bi} is touching the drop of radius c when its center is at the interior point xE K only if i J(x). Now, define y= - -cT/ icll (6) if J(x) = 0 (i.e., if (A.x - bi)/I Ai, | > e for all i = 1 to m), when the drop is in a position with its center at x, it will move in the gravitational direction y~. The distance that it will move in this direction is C (Ai x - b) - ||A I|A| 0 = minimum -- ---- b ----- -Ai.: 1 i i < m and i such _Ai.Y0 (7) that Ai. < 0 where we adopt the convention that the minimum in the empty set is +. If 0 = + W in (7), then the drop continues to move indefinitely along the half-line {x + Xy0: X. 0}, and z(x) is unbounded below on this feasible half-line, terminate. If 0 is finite in (7), at the end of this move, the drop will be in a position with its center at x + Oy0, touching the boundary of K, and it will 3

either halt (see the conditions for this, discussed later on) or change direction into the gravitational direction at x + Oy~ and move in that direction. When x is such that J(x) f 0, that is, min{(Ai x - bi)/llAi.l: i = 1 to m} = E (8) the direction that the drop will move next, called the gravitational direction at x, can be defined using many different principles. One principle to define the gravitational direction at x, where x is an interior point of K satisfying (8) is by the following procedure, which may take several steps. STEP 1: If the drop moves in the direction y~ from x, the position of its center will be x + Xy for some X > 0. Since (4) holds, the ith constraint will block the movement of the drop in the direction y0, only if i J(x) and Ai. y 0. Define J {i: i- J(x), and Ai y~ < O} CASE 1: J1 = 0: If J1 = 0, y is the gravitational direction at x, and the distance it can move in this direction is determined as in (7). CASE 2: J1 f 0: If J1 + 0, each of the constraints Ai.x 2 bi for i Ji, is currently blocking the movement of the drop in the direction yO. Define T1 = J1, and let D1 be the matrix of order IT11 x n whose rows are Ai. for iE T1. Let E1 be the submatrix of D1 of order (rank of D1) x n, whose set of rows is a maximal linearly independent subset of row vectors of D1. Let P1 = {i: Ai. is a row vector of El}. So P1CT1. Let F1 be the subspace {x: D1X 0 } = {x: Elx = 0}, F1 is the subspace corresponding to the set of all constraints which are blocking the movement of the drop in the direction y. Let E1 be the orthogonal projection of yO in the subspace F1, that is = (I - ET(E1ET E)yO. (9 4

SUBCASE 2.1: E1 0: If 1 f O, let y1 = E1/l 111, go to Step 2. SUBCASE 2.2: 1 = 0: If 1 = 0, let the row vector i = (pi: ie P1) = -IIc I((E1ET)- ElyO)T. Then lE = c. SUBCASE 2.2.1: =1 = 0 and p > 0: If p >a O, define the row vector i = (ir) by 1ri 0, if i 0 P1 i= i, if if P1 Then i is a basic feasible solution to the dual of (1). In this case, as will be shown later on, the drop halts in the current position, it cannot roll any further, under the gravitational force. SUBCASE 2.2.2: 01 = 0, p 4 0: If 1 = 0 and u 4 O, delete the i corresponding to the most negative pi from the set P1 (any other commonly used rule for deleting one or more of the i associated with negative i from P1 can be applied in this case). Redefine the matrix E1 to be the one whose rows are A. for i in the new set P1, compute the new orthogonal projection E' as in (9) using the new E1 and repeat Subcase 2.1 or 2.2 as appropriate with the new 1. GENERAL STEP r: Let yr 1 be the direction determined in the previous step. Define Jr {i: iE J(x) and Aiyr-1 < 0}. CASE 1: Jr = 0: If Jr = 0, yr1 is the gravitational direction at x, and the distance the drop can move in this direction is determined as in (3) with yr-1 replacing yO r CASE 2: Jr + 0: Define Tr = - Js and let Dr be the matrix of order ITrl s=1 x n whose rows are Ai for iE Tr. Let Er be the submatrix of Dr of order (rank of Dr) x n, whose set of rows is a maximal linearly independent subset of row 5

vectors of Dr. Let Pr = {i: Ai. is a row vector of Er}. Let Fr be the supspace {x: Drx = 0} = (x: Erx = 0}. Let Er be the orthogonal projection of y0 in the subspcae Fr, that is r = (I - ET(ErET)-1ErJy SUBCASE 2.1: r + 0: Let yr = |r/llr'1, to go Step r+1. SUBCASE 2.2: r = 0: Let p = (i: i Pr) | | IcI((ErE T)-1 rT SUBCASE 2.2.1: Er = 0, and p > O: Define i = (7i) by 0i = 0 for i Pr = i' ifor i Pr ir is a basic feasible solution to the dual of (1). In this case the drop halts, it cannot roll any further under the gravitational force. SUBCASE 2.2.2: Er = 0 and p 4 O: If Er = 0 and p 4 O, proceed exactly as under Subcase 2.2.2 described under Step 1, with Pr replacing P1. It can be shown that this procedure does produce the gravitational direction at x, finitely, if the drop can move at all. We are currently working on developing efficient methods for choosing the index set Pr of maximal linearly independent subset of row vectors of Dr, in Case 2, and on the best strategies for deleting a subset of constraints associated with negative Pi in Subcase 2.2.2. We are also looking at other principles for defining the gravitational direction at the interior point x of K. CONDITIONS FOR THE HALTING OF THE DROP Let e be the radius of the drop and x~ K satisfy (4). We have the following theorem. THEOREM 1: When the center of the drop is at x, it halts iff J(x) defined in (5) is f 0, and there exists a dual feasible solution = (-i.) for the dual of (1) satisfying 6

iri = 0, for all i J(x) (10) PROOF: The drop will halt when its center is at x, iff there exists no direction at x along which the drop could move within the interior of K, that will slide its center on a line of decreasing objective value for some positive length. That is, iff there exists no y satisfying cy < 0 (Ai.(x + Xy) - bi)/| Ai. || a, i = 1 to m for 0 AX < a, for some a > 0. Since x satisfies (4), and from the definition of J(x) in (5), this implies that the drop will halt when its center is at x iff the system Aiy a 0, for all iE J(x) cy < 0 has no solution y. By the well known Farkas' lemma (see, for example [2]) this holds iff there exists a r - (ir-: i = 1 to m) feasible to the dual of (1) satisfying (10). WHAT TO DO WHEN THE DROP HALTS? THEOREM 2: Suppose the drop of radius e halts with its center at xE K. Then the LP (1) has a finite optimum solution. Let z be the optimum objective value in (1). Let ir - (ri) be the dual feasible solution satisfying (10) guaranteed to exist by Theorem 1. Then cx b + ij(x) i (11) icJ(x) 1 and cx < z+ + - IT (12) i(J(x)i 7

PROOF: If the drop halts, by Theorem 1, the dual of (1) is feasible. So, the LP (1) has a finite optimum solution by the duality theory of LP. Consider the perturbed LP minimize z(x) = cx subject to A. b, for i J(x) Ai x > (13) bi + ~, for if J(x) The hypothesis in the theorem implies that x, iT, together satisfy the primal, dual feasibility and the complementary slackness optimality conditions for (13) and its dual. Hence, by the duality theorem of LP, (11) holds. Also, by the weak duality theorem of LP, (12) holds. Hence, if the drop halts with its center at position x, and a ir satisfying (10) is found, and e Z -I). is small, then x can be taken as a near optimum ieJ() 1 solution to (1) and the algorithm terminated. Also, in this case wr is an optimum solution for the dual of (1), and the true optimum solution of (1) can be obtained by well known pivotal methods that move from x to an extreme point without increasing the objective value (see [2]). THEOREM 3: Suppose the drop of radius E halts with its center at x~ K. If the system of equations Ai.x = bi, i6 J(x) (14) has a solution R which is feasible to (1), then R is an optimum feasible solution of (1). PROOF: Let 1r be the dual feasible solution satisfying (10) guaranteed by Theorem 1. It can be verified that X, w together satisfy the complementary slackness optimality conditions for (1) and its dual, so R is an optimum solution for (1). In this case it is optimum to the dual of (1). U 8

If the drop of radius c halts with its center at xE K, and there exists no solution to the system of equations (14) which is feasible to (1), then this drop is unable to move any further down in K under the gravitational force, eventhough it is not close to an optimum solution for (1). See Figure 1. 9

x Figure 1: The set K is on the side of the arrow marked on each constraint. The gravitational force is pulling the drop straight down, but it cannot move any further, because it is squeezed between hyperplanes 1 and 2. 10

Suppose the drop of radius e halts with its center at x. If the system Ai. x bi, i c J(x) (15) has no feasible solution, the gravitational method reduces the radius of the drop, see below, keeping the center at x, and continues. On the other hand, suppose the drop of radius e halts with its center at x, and the system (15) is feasible. Let E be the matrix whose rows form a maximal linearly independent subset of rows of {Ai.: i E J(x). Then the nearest point to x in the flat {x: Ai. x bi, i 6 J(x)} is x - x + ET(EET) 1 (d-Ex) where d is the column vector of bi for i such that Ai is a row of E. If x is feasible to (1), then by Theorem 3, x is an optimum feasible solution for (1) and the method terminates. Otherwise, at this stage the gravitational method reduces the radius of the drop (for example, replace e by c/2), keeping the center at x, and traces the locus of the center of the new drop as it now begins to fall under the influence of gravity again. The same process is repeated when the new drop halts. See Figure 2 for an illustration of the path of the drop in a convex polyhedron in R3. 11

-c direction 0. Figure 11.9 Path of the drop in the gravitational method in a convex polyhedron in R3 1~2

References 1. K. G. Murty, "A new interior variant of the gradient projection method for linear programming" working paper 85-18, Department of Industrial and Operations Engineering, The University of Michigan, May 1985. 2. K. G. Murty, Linear Programming, 1983, Wiley, New York. 13