A Fast Algorithm for the Two-Variable Integer ... - Semantic Scholar
A Fast Algorithm for the Two-Variable Integer Programming Problem SIDNIE DRESHER FElT I. T T Programming, Stratford, Connecttcza
Abstract. An algorithm that solves any two-variable integer programming problem is presented. A constant word-length model for the data is assumed. The complexity for a problem with m constraints and word length of L digits iS bounded by the maximum of two values. The first, which is O(mlogm) steps, is a bound on the complexity of finding the convex region bounded by the constraints, each step being an arithmetic orperation or a compare. The second, which Is O(mL) steps, is the complexity of solving m greatest-common-divisor problems. The algorithm finds a minimal binding set of constraints for any given problem, in addition to finding the soluuon set. A new method of solving three constraint problems is introduced. Categories and Subject Descriptors: F.2.0 [Analysis of Algorithms and Problem Complexity]: General; G l 6 [Numerical Analysis]: Optimizatlon--mteger programming; constrained optimization General Terms: Algorithm, Theory 1. Introduction
W e solve the p r o b l e m F i n d integers x a n d y m a x i m i z i n g the objective f u n c t i o n aolx -I- aOEY subjecl to the constraints allX + al2Y >- cl a21x + a22y >- c2
amlX + am2Y > Cm.
T h e au a n d c, are integers. There are no sign restrictions, a n d the region in the plane d e t e r m i n e d b y the constraints m a y be b o u n d e d , u n b o u n d e d , or e m p t y . By a solution p o i n t for the p r o b l e m , we m e a n a pair (x, y ) which satisfies the constraints a n d yields the o p t i m a l (in this case m a x i m a l ) value o f the objective function. T h e set o f all solution points is called the solution set for the p r o b l e m . Part of this research was completed while the author was a visiting fellow at the Department of Computer Soence, Yale Umversity, New Haven, CT. The remainder of the work was supported by I.T.T. Programming. Author's address: I.T.T. Programming, 1000 Oronoque Lane, Stratford, CT 06497. Permission to copy without fee all or part of this material is granted provided that the copies are not made or distributed for direct commercial advantage, the ACM copyright notice and the title of the puhhcation and its date appear, and nonce is given that copying is by permission of the Association for Computing Machinery. To copy otherwise, or to republish, requires a fee and/or specific permission. © 1984 ACM 0004-5411/84/0100-0099 $00.75 Journal of the A ~ t a U o n for Compuung Machinery,Vo|. 3 I, No l, January 1984, pp 99-113
100
SIDNIE DRESHER FElT
We assume a constant word length model for the data; that is, L = word length -- max{log base 2 of datum}, total size of data = m x 3 x L + 2 x L. Complexity will be measured here by counting the number of steps in the calculation. A step is an arithmetic operation or a comparison. Some calculations will involve greatest-common-divisor operations, which require at most O(L) steps. The number of steps to carry out the algorithm is dominated by two processes. The first involves finding the convex region bounded by m constraints. By this is meant finding which constraints actually bound, and finding all vertices. The complexity of this procedure is O(m log m) arithmetic steps. The second process is similar to solving at most m greatest-common-divisor problems. The complexity bound is O ( m L ) steps, that is, a number of steps that is linear in the size of the data. Kannan, in [3], presented an algorithm which was polynomial in the size of data for the special case of nonnegative (a,j) and (c,), and the additional constraints x, y ___ 0. Scarf, in [5], presented a solution that involved solving each threeconstraint subproblem. (There are O(m 3) subproblems.) The results in the present paper grew out of the author's efforts to streamline Scarf's solution to a threeconstraint problem and to reduce the number of three-constraint problems solved. It turns out that it is necessary to solve only one three-constraint problem. Here we give a rough outline of the method that will be used. The algorithms will be fully described in the remaining sections. A pair of integers (x, y) will be called a lattice point. The first step is to transform the data to a standard form' via a change of coordinates given by a 2 x 2 integer matrix of determinant +1 or - 1 (i.e., a unimodular transformation). After the transformation, to maximize the objective function, it will suffice to find a lattice point satisfying all constraints whose x coordinate is minimal. All solution points will lie on a vertical line x = c, which we will call the solution line. Next, the convex region bounded by the constraints will be found, and all constraints not bounding this region will be discarded. 2 At this point we can find out whether the maximum of the objective function is infinite, or if the region is empty and no solution is possible. Some other trivial cases can also be solved immediately. We then shall have a region similar to the one in Figure 1 with constraint set C = (bj . . . . .
bl, hi . . . . , ]'13).
The h bound above and the b bound below. Next some of the constraints will be discarded, leaving a constraint set c ' = (b,, b,+, . . . . .
hj, hj+, . . . . .
hj),
which has the same solution set. Further, b, and h~ will be in some sense essential constraints. J The assumption of uniform word size for the data is needed in this step to prevent the total size of the data from increasing by a factor of m (which could happen if most digits were concentrated in one or two inputs). 2 Eric Denardo has pointed out that a solution could be found by triangulating this region and solving the resulting three-constraint problems.
Two-Variable Integer Programming Problem
101
FIG. I.
A constraint set in standard form.
A typical discard of a constraint will occur when the shaded region in Figure 1 contains no lattice point. (The fact that the region is lattice point free can be discovered by solving a two-constraint problem.) The region will be enlarged by'a discard, but no additional lattice points will be enclosed. It has been pointed out by Bell [ 1] and Scarf [5], that if a solution point exists, there is a choice of at most three constraints which determine the solution line. It will turn out that b, and hj can be chosen as two of the constraints. Furthermore, one of the following will be true. 1. There is no solution point for the problem. 2. All solution points for the original problem lie on the solution line for (b,, hi). 3. A third constraint can be found so that all solution points for the original problem lie on the solution line for (b,, hi, third constraint). A simplified version of a method of Scarf is used to find the solution line for a two-constraint problem. A new method of solving three-constraint problems will be introduced and used here.
2. A Standard Form for the Problem By a unimodular change of coordinates we mean a linear transformation given by a 2 x 2 integer matrix whose determinant is + 1 or - 1 . Let
all A
~
°
a12 •
•
a°la°2 I
\ami am2/
That is, A is an array containing the coefficients for the objective function and all constraints. For a lattice point (x, y), let p = (x, y)' and z = Ap. Then we can restate the problem as follows: Find p so that among all z = Ap with z, ~ c, i = 1. . . . , m, z0 is maximal. If the coordinates in the x, y plane are changed by a unimodular transformation U, the integer lattice is mapped 1-1 onto the integer lattice, and the set of points
z = AUU-~p is not changed. It suffices to solve the problem for any A U, where U may be chosen to give AUa convenient form. Ifq is a solution point for AU, then Uq is a solution point for A.
102
SIDNIE DRESHER FEIT
There is a unimodular change of coordinates which transforms the problem to Find a pair of integers (x, y) with x minimal, subject to (1) x _ d (d possibly -oo), (2) x _< d (d possibly +oo), (3) (x, y) satisfies a set of constraints, each given by a line whose slope is positive. The entries in U will have a worst case bound of 3L + l digits. This can be done quite easily. (The method has been used by Scarf in [5].) Let g be the (positive) greatest c o m m o n divisor of the coefficients aol, ab2 of the objective function. Then there are integers r and s with raoi + sao2 = g,
I rl < l ao2 I,
I s l < l a o , I.
Hence, either of the changes of coordinates,
o, (_:
-aO2g / ao_2 g
transforms the objective function to (-g, 0). 3 At this point, either transformation can be chosen. However, later it may turn out that a particular orientation is needed, and we may have to reverse the choice. After multiplying A by one of these matrices, examine the sign patterns of t h e constraints, and divide them into four types.
I. (+, o) 2. ( - , o ) 3. (*, +) constraints that bound below 4. (*, - ) constraints that bound above. Let LuJ = greatest integer _ u, and ru] = least integer _ u. Any constraint of type 1 is of the form
orequvalen,,y x:f ] All constraints of type 1 can be replaced by the single most binding one. Similarly, if there are constraints of type 2, they can be replaced by the most binding one, in the form
Now a number of trivial cases can be cleared away. Examples of each case are shown in Figure 2. (a) A type l constraint contradicts a type 2 constraint. (b) All constraints are of type l or 2--solution trivial. (c) All constraints are of type 2 and 3, or 2 and 4 only. Then x is not bounded below, and the original objective function is unbounded. 3 Note that m the worst case, t h e n u m b e r o f digits c a n increase to 2 L + I.
103
Two-Variable Integer Programming Problem -.IP,. . q [ . -
(a)
(c ~)
(b)
(c)
(~)
(e)
f)
FIG. 2.
Cases for which a solution is ~mmediate.
(d) There are constraints of type 1 and 3, but none of type 4. Solve easily by finding the highest intersection of a type 3 constraint with the most binding type 1 constraint. (e) There are constraints of type 1 and 4, but none of type 3. Solve easily by finding the highest intersection of a type 4 constraint with the most binding type 1 constraint. Before looking at a final trivial case, it is convenient to perform the next step. Compute tile slope of each line of type 3 (lower bounding) and 4 (upper bounding). Sort type 4 lines in order of decreasing slope and sort the type 3 lines by increasing slope. Now we wish to give all constraint lines of type 3 and 4 a positive slope. To do this, choose a sufficiently large value of T. Subtract T x the y coordinate from the x coordinate, 4 so that every constraint of type 3 has the form ( - , +) and every constraint of type 4 has the form (+, -). Let W be the array obtained from U by subtracting T x column 2 from column 1. Then W is the unimodular transformation can-ring A to the desired form. Now consider the case (f) There is no type 1 constraint and the steepest upper bounding line has slope less than or equal to the slope of the least steep lower bounding line. In this case, either the region is empty or the x values are not bounded below (i.e., the original objective function is unbounded). If none of (a)-(f) hold, there must be constraints of type 3 and type 4. Further, either there is a (+, 0') constraint, or the steepest upper bounding constraint has /
4 T = 1 + max[ I slope I 1 over all constramts will work. Note that when the application o f U is followed by this transformation, 4 L + 3 is an upper b o u n d on the number o f digits in a word. The assumption o f constant word snze will not be essentml m any o f the remaining steps.
104
SIDNIE DRESHER FElT h2
FIG. 3. The convex region R determined by the constraints.
b~
--n
slope greater than the slope of the flattest lower bounding constraint. This implies that there are lower bounds on the x and y coordinates of points in the convex region bounded by the transformed constraint set. If there is a ( - , 0) constraint, set it aside and use it to test the solution for the remaining constraints. If there is a (+, 0) constraint, give it a large d u m m y slope and put it at the top of the type 4 list. A well-known algorithm of complexity O(mlogm) can now be applied to find the convex (possibly unbounded) polyhedron determined by the constraints and the subset of constraints that actually bound. Note that if this polyhedron is empty, we are done. Also, if it~is a line, a ray, or a segment, the problem can be solved immediately. Otherwise, consider Figure 3. Let R be the dosed convex region in the plane determined by the constraints. The algorithm that determines R will return the x and y coordinates of each vertex of R. If a line x --- n meets R in a nonempty segment, the constraints met by the segment and the coordinates of all lattice points in R lying on the segment are easily found. Let b~, ..., bl be the lower lines and h~, ..., hj be the upper lines that together bound the region. Recall that any ( - , 0) constraint has been set aside, and will be used to test the solution of the problem determined by the remaining constraints. From now on, the problem will be in the standard form Find {(x, y)} such that x is minimal and (x, y) lies in a closed convex region bounded by at most one vertical line and a set of lines with positive slope. Further, the region will be such that {yl(x, y) is a solutionl is bounded. The only vertical constraint line that will be included will be a (+, 0) constraint. In order to solve a standard problem, it will be necessary to solve several twoconstraint subproblems and one three-constraint problem.
3. The Discard Algorithm Let C = (b~, ..., bl, hi . . . . . hj) be the constraints for a standard problem, where the b are ordered with increasing slope and the h with decreasing slope. We will delete constraints bl . . . . . b,-i and h~. . . . , hj_l so that the set of constraints c ' = (b, . . . . .
hi, h~ . . . . .
hj)
has the same solution set as the original problem, and b, and hj each block a lattice point accepted by all of the other constraints. The choice of i and j is not unique, in general. Consider the constraint set in Figure 4. If both shaded regions are lattice point free, but there is a lattice point p as shown, then either (b2, hi, h2) or (b~, b2, h2) will satisfy the criteria, but we may not discard both b~ and h~. In the algorithm, a convention will be adopted for the choice. We will first perform all possible discards for lower lines, and then perform discards for upper lines.
Two-Variable Integer Programming Problem
105
b2 FiG. 4. Choice of discards for constraint set C.
lal
If b corresponds to the constraint
- a t x + azy >_ c,
a~ > O, a2 > O,
we define the "objective function parallel to - b " as the linear expression
a ~ x - a2y. Similarly, if h corresponds to the constraint
ajx - a 2 y >- c,
ai > O,
a2 > O,
we define the " objective function parallel to - h " as the linear expression - a l x + a2y.
Algorithm DISCARD (1) t = l , J =
1.
(2) If z = I or slope b,+~ - slope hj, go to (3). Else solve the integer programming problem with constraints b,+l and h~ and objective function parallel to -b,. If all solution points lie on or above line b,, then discard b,, set t -- t + 1, and repeat (2). Else go to (3). (3) I f j - J or slope hj+l _< slope b,, STOP. Else, solve the integer programming problem with constraints b, and h~+~and objective function parallel to -hj. If all solution points lie on or below line hi, then discard hi, set j = j + 1, and repeat (3). Else STOP. The first i - 1 lower constraints and the first j - 1 upper constraints are now discarded, since the remaining constraints will enclose exactly the same set of lattice points, b, and h~ each block a lattice point which is accepted by all other constraints. We call the pair (b,, hi) a frame for constraint set (b, . . . . . bl, hj . . . . . h j); b, is called the lower frame line and hj is called the upper frame line. There can be at most m - 2 discards carried out by the algorithm. Each of these requires the solution of a two-constraint problem. Sections 5 and 7 will show that the complexity of a two-constraint problem is O(L) arithmetic steps. If b,+~ is parallel to h~ or h~+l is parallel to b,, then the parallel lines algorithm of Section 6 should be applied now, to determine whether the parallel lines enclose any lattice points. If they do not, then there is no need to proceed. Therefore, we assume in the remainder of the paper that if one of these pairs is parallel, then lattice points are enclosed.
4. The Four-Constraint Theorem and the Algorithms The arguments that follow are all based on the following theorem which was proved independently by Bell [1] and Scarf [5]. 2 ~ CONSTRAINT THEOREM. Let Q be the closed convex set in R n given by the intersection o f a set of closed half spaces {H,}. I f Q contains no lattice points, then
106
SIDNIE DRESHER FElT
FIG. 5. Lattice point free region to the left of a solution. !
there is a subset of at most 2" of the half spaces, whose intersection also contains no lattice point. In the two-dimensional case, this means that given a closed subset Q of the plane determined by a set of constraint lines, if Q contains no lattice points, then we can choose a set consisting of at most four of the constraints so that there are no lattice points in the set which they determine. Now suppose that we are given a two-variable integer programming problem in standard form, with solution line x = x*. If we add the new constraint X _< X* -- ~,
where ~ is a very small positive number, then a lattice-point-free region is obtained (see Figure 5). The new constraint cannot be deleted without admitting lattice points. However, it must be possible to delete all but three of the original constraints without introducing new lattice points. Thus there is a subset of at most three of the original constraints whose solution points all lie on x = x*. Note that the set of solution points for these three constraints may be larger than the original set. However, all of the solution points will lie on x = x*, and the true solution set can be found by calculating the intersection points of x = x* with the boundary of the region R. We proceed by finding a set of at most three constraints, which will have the following property. I f R contains any lattice points, then the solution line for R is the solution line for the three constraints. To start off, we solve the standard problem with constraints (b,, h~). Let x --- x0 be the solution line, and let S denote the solution set for this problem. Suppose x --- Xo meets the boundary of the region R at constraint lines b' and h'. There are four possibilities: 1. Some point in S is accepted by b' and h' and hence by all constraints. Then this is a solution and we are done. 2. All points in S are rejected, and both b' and h' reject some point in S. Then the region bounded by (b,, hj, b', h') contains no lattice point, for if it did, we would have a region bounded by five constraint lines, none of which could be deleted without introducing a lattice point (see Figure 6a). 3. All points in S are rejected by b' (see Figure 6b). 4. All points in S are rejected by h'. The fourth case can be reduced to the third by reversing the choice of orientation, which turns the picture upside down. The upper and lower bounding lines are
Two-Variable Integer Programming Problem
107
I I
FIG. 6.
(a) Case 2; (b) case 3.
(a)
I q¢
×*
X.=X* (a)
•
FIG. 7.
(b)
Lattice points blocked by hi.
interchanged, and if hj is a vertical line with large dummy slope, it is carried to a lower bounding line with small dummy slope. Thus it suffices to consider the third case. Let us assume that our problem has a nonempty solution set, lying on solution line x --- x*. Under this assumption, we will find three constraints which must have the same solution line. If no solution point for these constraints satisfies all of the original constraints, we will be able to conclude that no lattice point satisfied all of the original constraints. We can observe immediately that b, is an essential, binding constraint. For, as a result of the DISCARD algorithm, we have i = / , or b,+l has slope greater than or equal to the slope of hi, or else lines (b,, b,+l, hi) trap a lattice point to the left of x = x0, the solution line for constraints (b,, hi). (The latter is true because some lower bounding line rejects the solution set S, which implies that b, intersects bz+~ to the left of x0.) In any case, b, cannot be discarded without changing the solution line. The solution line for the problem must be obtainable from b,, some bk that rejects S, and one upper bounding line. We claim that it suffices to use h~ as the upper bounding line. For, from the discard algorithm, either j = J, or slope hj+~ is less than or equal to slope bz, or else at least one lattice point is trapped in the region bounded by (bi, hj, hj+O. In the first two instances, h~ is an essential binding constraint, and cannot be dropped without changing the solution. In the last case, if a trapped lattice point lies to the left o f x --- x*, then hj cannot be discarded without changing the solution. If all trapped points lie on or to the right of x = x*, then hj is the upper boundary line for x = x*, and all of the other upper lines can be ignored (see Figure 7). It remains to choose a third constraint from among the lower bounding lines that reject S. This is done by using a discard algorithm, with objective function parallel to b,. Specifically, let (Ii, ..., Ir) be the lower bounding lines which reject
108
S1DNIE DRESHER FEIT
S. We try to find one of these lines which blocks a lattice point not blocked by the others.
Algorithm DISCARD2 (1) k = I. (2) If k = r, then STOP. (3) If there is no lattice point below 1~in the region bounded by b,, lk, and lk+l, then discard constraint h, set k = k + 1, and go to (2). Else STOP. We claim that it suffices to solve the standard three-constraint problem (hj, b,, lk). The argument showing that/k is an essential constraint is similar to the argument given for hi. From DISCARD2, either k --- r, or else at least one lattice point is trapped in the region bounded by (b,, lk, l~+~). In the first case, lk is an essential binding constraint, and cannot be dropped without changing the solution line. In the second case, if a trapped lattice point lies to the left of x = x*, lk cannot be discarded without changing the solution line. If all trapped points lie on or to the right o f x = x*, then lk is the lower boundary line for x = x*, and (lk+t, . . . , lr) can be discarded without changing the solution. Finally, we use the three-constraint algorithm to solve (b,, lk, h~). If these constraints bound a lattice-free region, then no solution exists for the original problem. If the solution set lies on a line x = n, the solution set for R is contained in the intersection of this line with R. If no lattice points lie on the intersection segment, then R contains no lattice points.
5. The Double Greatest Common Divisor (DGCD) This section and the next contain the technical details for carrying out the computations that lie at the core of the algorithm. The technique for solving two-constraint problems will consist of manipulating the positions of the constraint lines in the plane, using lattice-preserving maps, until a position is reached for which a solution is immediately visible. The only operation used will be the addition of a positive multiple of one column of an array to another. Thus the signs of all subdeterminants will be preserved throughout the algorithm. This operation will be repeated at most O(L) times. An array
A =
(alla2) a2t
a22
a31
a32
(e e)
with sign pattern
+
+
()
and with -a~ l/al2 > -a2t/a22 > -a3t/a32 is to be carried to array
E =
e2~ e22
e31
e32
with sign pattern
+
+
..~
where eli --- el2 and e3t -< e32. The number of steps required for the simple algorithm presented here is bounded by the number of steps required to find the greatest common divisor of ( a l t, at2). The total effect of the algorithm is to multiply A by a 2 x 2 array whose entries are bounded in size by max{a,j].
109
Two-Variable Integer Programming Problem
The algorithm requires that we assign a sign to 0. We will always regard 0 as having a "--" sign. Let us concentrate first on rows 2 and 3. The algorithm below will carry these rows to an array such that sign p a t t e r n = ( +
+),
column s u m = (++).
To simplfy notation, write rows 2 and 3 as
(s) rl
-r2
s2
DGCD Algorithm (1) Let T = I.s~/rd. Set column 2 = column 2 + T x column 1. (2) Let T = tr2/s21. If T = 0, STOP. Elge set column I = column I + T x column 2. (3) Let T : ts,/rd. If T = 0, STOP. Else set column 2 = column 2 + T x column 1, and go to (2)• After this algorithm ends, replace column 1 with column 1 + column 2 to get the desired sign pattern in rows 2 and 3. What pattern will row 1 have now? Because of the sign values of the suNleterminants and the fact that the last replacement occurred in column 1, the pattern is either ( - - ) with e,, __ e,2 and e3, -< e32 and we are done, or else it is still ( - +). In the latter case, we use an algorithm almost identical to the [N3CD. Represent rows 1 and 2 as
r2
-s2
In the earlier case, only negative entries could become 0, which caused no trouble. Now we prevent a " + " element from becoming 0 by defining [IxJJ = greatest integer strictly less than x. Replacing t / by [l J J in the D G C D produces an algorithm that will transform rows 1 and 2, preserving the sign pattern, so that the column sum of the final result is negative. Now replace column 2 by column 1 + column 2, and the desired sign pattern and order relations will be obtained. 6. The Parallel Lines Algorithm
The main result of this section is that a standard problem bounded by a vertical line v and a pair of parallel lines (see Figure 8) ex - f y = c,,
ex - f y = c2,
e, f > 0
can be solved in at most 1 + log2e + logzfsteps. Since v has the form x = d, v can be replaced by x = [d] without changing the solution. The point p in the figure is the first lattice point on v below h'. It is convenient to perform a translation of the coordinates so that p is the origin. If (0, 1) is on or below h, we are done. If not, perform a greatest c o m m o n divisor algorithm on (e, f ) to obtain r and s such that er-fs=g,
__f, e Irl < g ]sl Slope h (Figure 14). Solve the problem with objective function h and two constraints v and bl. If there is no lattice point on or below h, then there are no lattice points in the region. If q is a solution point for this
Two-Variable Integer Programming Problem
113
h
FIG. 12. Case I: slope b~ < slope h.
*q
FIG. 13
Case II' slope bt =slopeh.
I
v
~
h
h'
/ I
FIG. 14. Case 111. slope b~ > s l o p e h
problem, let h' be the line through q parallel to h and solve the parallel lines problem bounded above by h and below by h'. Since q is feasible, there will be a solution and it will satisfy constraint b, and solve the three-constraint problem. ACKNOWLEDGMENTS. I wish to thank Dan Gusfield for many informative conversations, and the referee for several helpful comments and suggestions. REFERENCES (Note:
References [2] and [41 are not cited m the text.)
1. BELt, D.E. A theorem concerning the integer lattice. Stud Appl. Math. 56 (1977), 187-188. 2. HIRSCHBERG, D.S., AND WONG, C.K A polynomml-Ume algorithm for the knapsack problem with two variables. J ACM 23, 1 (Jan. 1976), 147-154. 3. KANNAN,R. A polynomial algorithm for the two-variable integer programming problem. J. ACM 27, 1 (Jan. 1980), !18-122. 4. LENSTRA,H W Integer programming wah a fixed number of variables. Rep. 81-03, Mathematisch lnsutuut, Amsterdam. 5. SCAR~,H.E. ProducUon sets with indlvlsibiliues. Econometrlca (1981). RECEIVED APRIL 1981; REVISEDMARCH 1983; ACCEPTEDMARCH 1983 Journal of the Assocmtlon for Computing
Maehlnery,Vo! 31, No. 1, Janlaartt 1984
Abstract. An algorithm that solves any two-variable integer programming problem is presented. A constant word-length model for the data is assumed. The complexity for a problem with m constraints and word length of L digits iS bounded by the maximum of two values. The first, which is O(mlogm) steps, is a bound on the complexity of finding the convex region bounded by the constraints, each step being an arithmetic orperation or a compare. The second, which Is O(mL) steps, is the complexity of solving m greatest-common-divisor problems. The algorithm finds a minimal binding set of constraints for any given problem, in addition to finding the soluuon set. A new method of solving three constraint problems is introduced. Categories and Subject Descriptors: F.2.0 [Analysis of Algorithms and Problem Complexity]: General; G l 6 [Numerical Analysis]: Optimizatlon--mteger programming; constrained optimization General Terms: Algorithm, Theory 1. Introduction
W e solve the p r o b l e m F i n d integers x a n d y m a x i m i z i n g the objective f u n c t i o n aolx -I- aOEY subjecl to the constraints allX + al2Y >- cl a21x + a22y >- c2
amlX + am2Y > Cm.
T h e au a n d c, are integers. There are no sign restrictions, a n d the region in the plane d e t e r m i n e d b y the constraints m a y be b o u n d e d , u n b o u n d e d , or e m p t y . By a solution p o i n t for the p r o b l e m , we m e a n a pair (x, y ) which satisfies the constraints a n d yields the o p t i m a l (in this case m a x i m a l ) value o f the objective function. T h e set o f all solution points is called the solution set for the p r o b l e m . Part of this research was completed while the author was a visiting fellow at the Department of Computer Soence, Yale Umversity, New Haven, CT. The remainder of the work was supported by I.T.T. Programming. Author's address: I.T.T. Programming, 1000 Oronoque Lane, Stratford, CT 06497. Permission to copy without fee all or part of this material is granted provided that the copies are not made or distributed for direct commercial advantage, the ACM copyright notice and the title of the puhhcation and its date appear, and nonce is given that copying is by permission of the Association for Computing Machinery. To copy otherwise, or to republish, requires a fee and/or specific permission. © 1984 ACM 0004-5411/84/0100-0099 $00.75 Journal of the A ~ t a U o n for Compuung Machinery,Vo|. 3 I, No l, January 1984, pp 99-113
100
SIDNIE DRESHER FElT
We assume a constant word length model for the data; that is, L = word length -- max{log base 2 of datum}, total size of data = m x 3 x L + 2 x L. Complexity will be measured here by counting the number of steps in the calculation. A step is an arithmetic operation or a comparison. Some calculations will involve greatest-common-divisor operations, which require at most O(L) steps. The number of steps to carry out the algorithm is dominated by two processes. The first involves finding the convex region bounded by m constraints. By this is meant finding which constraints actually bound, and finding all vertices. The complexity of this procedure is O(m log m) arithmetic steps. The second process is similar to solving at most m greatest-common-divisor problems. The complexity bound is O ( m L ) steps, that is, a number of steps that is linear in the size of the data. Kannan, in [3], presented an algorithm which was polynomial in the size of data for the special case of nonnegative (a,j) and (c,), and the additional constraints x, y ___ 0. Scarf, in [5], presented a solution that involved solving each threeconstraint subproblem. (There are O(m 3) subproblems.) The results in the present paper grew out of the author's efforts to streamline Scarf's solution to a threeconstraint problem and to reduce the number of three-constraint problems solved. It turns out that it is necessary to solve only one three-constraint problem. Here we give a rough outline of the method that will be used. The algorithms will be fully described in the remaining sections. A pair of integers (x, y) will be called a lattice point. The first step is to transform the data to a standard form' via a change of coordinates given by a 2 x 2 integer matrix of determinant +1 or - 1 (i.e., a unimodular transformation). After the transformation, to maximize the objective function, it will suffice to find a lattice point satisfying all constraints whose x coordinate is minimal. All solution points will lie on a vertical line x = c, which we will call the solution line. Next, the convex region bounded by the constraints will be found, and all constraints not bounding this region will be discarded. 2 At this point we can find out whether the maximum of the objective function is infinite, or if the region is empty and no solution is possible. Some other trivial cases can also be solved immediately. We then shall have a region similar to the one in Figure 1 with constraint set C = (bj . . . . .
bl, hi . . . . , ]'13).
The h bound above and the b bound below. Next some of the constraints will be discarded, leaving a constraint set c ' = (b,, b,+, . . . . .
hj, hj+, . . . . .
hj),
which has the same solution set. Further, b, and h~ will be in some sense essential constraints. J The assumption of uniform word size for the data is needed in this step to prevent the total size of the data from increasing by a factor of m (which could happen if most digits were concentrated in one or two inputs). 2 Eric Denardo has pointed out that a solution could be found by triangulating this region and solving the resulting three-constraint problems.
Two-Variable Integer Programming Problem
101
FIG. I.
A constraint set in standard form.
A typical discard of a constraint will occur when the shaded region in Figure 1 contains no lattice point. (The fact that the region is lattice point free can be discovered by solving a two-constraint problem.) The region will be enlarged by'a discard, but no additional lattice points will be enclosed. It has been pointed out by Bell [ 1] and Scarf [5], that if a solution point exists, there is a choice of at most three constraints which determine the solution line. It will turn out that b, and hj can be chosen as two of the constraints. Furthermore, one of the following will be true. 1. There is no solution point for the problem. 2. All solution points for the original problem lie on the solution line for (b,, hi). 3. A third constraint can be found so that all solution points for the original problem lie on the solution line for (b,, hi, third constraint). A simplified version of a method of Scarf is used to find the solution line for a two-constraint problem. A new method of solving three-constraint problems will be introduced and used here.
2. A Standard Form for the Problem By a unimodular change of coordinates we mean a linear transformation given by a 2 x 2 integer matrix whose determinant is + 1 or - 1 . Let
all A
~
°
a12 •
•
a°la°2 I
\ami am2/
That is, A is an array containing the coefficients for the objective function and all constraints. For a lattice point (x, y), let p = (x, y)' and z = Ap. Then we can restate the problem as follows: Find p so that among all z = Ap with z, ~ c, i = 1. . . . , m, z0 is maximal. If the coordinates in the x, y plane are changed by a unimodular transformation U, the integer lattice is mapped 1-1 onto the integer lattice, and the set of points
z = AUU-~p is not changed. It suffices to solve the problem for any A U, where U may be chosen to give AUa convenient form. Ifq is a solution point for AU, then Uq is a solution point for A.
102
SIDNIE DRESHER FEIT
There is a unimodular change of coordinates which transforms the problem to Find a pair of integers (x, y) with x minimal, subject to (1) x _ d (d possibly -oo), (2) x _< d (d possibly +oo), (3) (x, y) satisfies a set of constraints, each given by a line whose slope is positive. The entries in U will have a worst case bound of 3L + l digits. This can be done quite easily. (The method has been used by Scarf in [5].) Let g be the (positive) greatest c o m m o n divisor of the coefficients aol, ab2 of the objective function. Then there are integers r and s with raoi + sao2 = g,
I rl < l ao2 I,
I s l < l a o , I.
Hence, either of the changes of coordinates,
o, (_:
-aO2g / ao_2 g
transforms the objective function to (-g, 0). 3 At this point, either transformation can be chosen. However, later it may turn out that a particular orientation is needed, and we may have to reverse the choice. After multiplying A by one of these matrices, examine the sign patterns of t h e constraints, and divide them into four types.
I. (+, o) 2. ( - , o ) 3. (*, +) constraints that bound below 4. (*, - ) constraints that bound above. Let LuJ = greatest integer _ u, and ru] = least integer _ u. Any constraint of type 1 is of the form
orequvalen,,y x:f ] All constraints of type 1 can be replaced by the single most binding one. Similarly, if there are constraints of type 2, they can be replaced by the most binding one, in the form
Now a number of trivial cases can be cleared away. Examples of each case are shown in Figure 2. (a) A type l constraint contradicts a type 2 constraint. (b) All constraints are of type l or 2--solution trivial. (c) All constraints are of type 2 and 3, or 2 and 4 only. Then x is not bounded below, and the original objective function is unbounded. 3 Note that m the worst case, t h e n u m b e r o f digits c a n increase to 2 L + I.
103
Two-Variable Integer Programming Problem -.IP,. . q [ . -
(a)
(c ~)
(b)
(c)
(~)
(e)
f)
FIG. 2.
Cases for which a solution is ~mmediate.
(d) There are constraints of type 1 and 3, but none of type 4. Solve easily by finding the highest intersection of a type 3 constraint with the most binding type 1 constraint. (e) There are constraints of type 1 and 4, but none of type 3. Solve easily by finding the highest intersection of a type 4 constraint with the most binding type 1 constraint. Before looking at a final trivial case, it is convenient to perform the next step. Compute tile slope of each line of type 3 (lower bounding) and 4 (upper bounding). Sort type 4 lines in order of decreasing slope and sort the type 3 lines by increasing slope. Now we wish to give all constraint lines of type 3 and 4 a positive slope. To do this, choose a sufficiently large value of T. Subtract T x the y coordinate from the x coordinate, 4 so that every constraint of type 3 has the form ( - , +) and every constraint of type 4 has the form (+, -). Let W be the array obtained from U by subtracting T x column 2 from column 1. Then W is the unimodular transformation can-ring A to the desired form. Now consider the case (f) There is no type 1 constraint and the steepest upper bounding line has slope less than or equal to the slope of the least steep lower bounding line. In this case, either the region is empty or the x values are not bounded below (i.e., the original objective function is unbounded). If none of (a)-(f) hold, there must be constraints of type 3 and type 4. Further, either there is a (+, 0') constraint, or the steepest upper bounding constraint has /
4 T = 1 + max[ I slope I 1 over all constramts will work. Note that when the application o f U is followed by this transformation, 4 L + 3 is an upper b o u n d on the number o f digits in a word. The assumption o f constant word snze will not be essentml m any o f the remaining steps.
104
SIDNIE DRESHER FElT h2
FIG. 3. The convex region R determined by the constraints.
b~
--n
slope greater than the slope of the flattest lower bounding constraint. This implies that there are lower bounds on the x and y coordinates of points in the convex region bounded by the transformed constraint set. If there is a ( - , 0) constraint, set it aside and use it to test the solution for the remaining constraints. If there is a (+, 0) constraint, give it a large d u m m y slope and put it at the top of the type 4 list. A well-known algorithm of complexity O(mlogm) can now be applied to find the convex (possibly unbounded) polyhedron determined by the constraints and the subset of constraints that actually bound. Note that if this polyhedron is empty, we are done. Also, if it~is a line, a ray, or a segment, the problem can be solved immediately. Otherwise, consider Figure 3. Let R be the dosed convex region in the plane determined by the constraints. The algorithm that determines R will return the x and y coordinates of each vertex of R. If a line x --- n meets R in a nonempty segment, the constraints met by the segment and the coordinates of all lattice points in R lying on the segment are easily found. Let b~, ..., bl be the lower lines and h~, ..., hj be the upper lines that together bound the region. Recall that any ( - , 0) constraint has been set aside, and will be used to test the solution of the problem determined by the remaining constraints. From now on, the problem will be in the standard form Find {(x, y)} such that x is minimal and (x, y) lies in a closed convex region bounded by at most one vertical line and a set of lines with positive slope. Further, the region will be such that {yl(x, y) is a solutionl is bounded. The only vertical constraint line that will be included will be a (+, 0) constraint. In order to solve a standard problem, it will be necessary to solve several twoconstraint subproblems and one three-constraint problem.
3. The Discard Algorithm Let C = (b~, ..., bl, hi . . . . . hj) be the constraints for a standard problem, where the b are ordered with increasing slope and the h with decreasing slope. We will delete constraints bl . . . . . b,-i and h~. . . . , hj_l so that the set of constraints c ' = (b, . . . . .
hi, h~ . . . . .
hj)
has the same solution set as the original problem, and b, and hj each block a lattice point accepted by all of the other constraints. The choice of i and j is not unique, in general. Consider the constraint set in Figure 4. If both shaded regions are lattice point free, but there is a lattice point p as shown, then either (b2, hi, h2) or (b~, b2, h2) will satisfy the criteria, but we may not discard both b~ and h~. In the algorithm, a convention will be adopted for the choice. We will first perform all possible discards for lower lines, and then perform discards for upper lines.
Two-Variable Integer Programming Problem
105
b2 FiG. 4. Choice of discards for constraint set C.
lal
If b corresponds to the constraint
- a t x + azy >_ c,
a~ > O, a2 > O,
we define the "objective function parallel to - b " as the linear expression
a ~ x - a2y. Similarly, if h corresponds to the constraint
ajx - a 2 y >- c,
ai > O,
a2 > O,
we define the " objective function parallel to - h " as the linear expression - a l x + a2y.
Algorithm DISCARD (1) t = l , J =
1.
(2) If z = I or slope b,+~ - slope hj, go to (3). Else solve the integer programming problem with constraints b,+l and h~ and objective function parallel to -b,. If all solution points lie on or above line b,, then discard b,, set t -- t + 1, and repeat (2). Else go to (3). (3) I f j - J or slope hj+l _< slope b,, STOP. Else, solve the integer programming problem with constraints b, and h~+~and objective function parallel to -hj. If all solution points lie on or below line hi, then discard hi, set j = j + 1, and repeat (3). Else STOP. The first i - 1 lower constraints and the first j - 1 upper constraints are now discarded, since the remaining constraints will enclose exactly the same set of lattice points, b, and h~ each block a lattice point which is accepted by all other constraints. We call the pair (b,, hi) a frame for constraint set (b, . . . . . bl, hj . . . . . h j); b, is called the lower frame line and hj is called the upper frame line. There can be at most m - 2 discards carried out by the algorithm. Each of these requires the solution of a two-constraint problem. Sections 5 and 7 will show that the complexity of a two-constraint problem is O(L) arithmetic steps. If b,+~ is parallel to h~ or h~+l is parallel to b,, then the parallel lines algorithm of Section 6 should be applied now, to determine whether the parallel lines enclose any lattice points. If they do not, then there is no need to proceed. Therefore, we assume in the remainder of the paper that if one of these pairs is parallel, then lattice points are enclosed.
4. The Four-Constraint Theorem and the Algorithms The arguments that follow are all based on the following theorem which was proved independently by Bell [1] and Scarf [5]. 2 ~ CONSTRAINT THEOREM. Let Q be the closed convex set in R n given by the intersection o f a set of closed half spaces {H,}. I f Q contains no lattice points, then
106
SIDNIE DRESHER FElT
FIG. 5. Lattice point free region to the left of a solution. !
there is a subset of at most 2" of the half spaces, whose intersection also contains no lattice point. In the two-dimensional case, this means that given a closed subset Q of the plane determined by a set of constraint lines, if Q contains no lattice points, then we can choose a set consisting of at most four of the constraints so that there are no lattice points in the set which they determine. Now suppose that we are given a two-variable integer programming problem in standard form, with solution line x = x*. If we add the new constraint X _< X* -- ~,
where ~ is a very small positive number, then a lattice-point-free region is obtained (see Figure 5). The new constraint cannot be deleted without admitting lattice points. However, it must be possible to delete all but three of the original constraints without introducing new lattice points. Thus there is a subset of at most three of the original constraints whose solution points all lie on x = x*. Note that the set of solution points for these three constraints may be larger than the original set. However, all of the solution points will lie on x = x*, and the true solution set can be found by calculating the intersection points of x = x* with the boundary of the region R. We proceed by finding a set of at most three constraints, which will have the following property. I f R contains any lattice points, then the solution line for R is the solution line for the three constraints. To start off, we solve the standard problem with constraints (b,, h~). Let x --- x0 be the solution line, and let S denote the solution set for this problem. Suppose x --- Xo meets the boundary of the region R at constraint lines b' and h'. There are four possibilities: 1. Some point in S is accepted by b' and h' and hence by all constraints. Then this is a solution and we are done. 2. All points in S are rejected, and both b' and h' reject some point in S. Then the region bounded by (b,, hj, b', h') contains no lattice point, for if it did, we would have a region bounded by five constraint lines, none of which could be deleted without introducing a lattice point (see Figure 6a). 3. All points in S are rejected by b' (see Figure 6b). 4. All points in S are rejected by h'. The fourth case can be reduced to the third by reversing the choice of orientation, which turns the picture upside down. The upper and lower bounding lines are
Two-Variable Integer Programming Problem
107
I I
FIG. 6.
(a) Case 2; (b) case 3.
(a)
I q¢
×*
X.=X* (a)
•
FIG. 7.
(b)
Lattice points blocked by hi.
interchanged, and if hj is a vertical line with large dummy slope, it is carried to a lower bounding line with small dummy slope. Thus it suffices to consider the third case. Let us assume that our problem has a nonempty solution set, lying on solution line x --- x*. Under this assumption, we will find three constraints which must have the same solution line. If no solution point for these constraints satisfies all of the original constraints, we will be able to conclude that no lattice point satisfied all of the original constraints. We can observe immediately that b, is an essential, binding constraint. For, as a result of the DISCARD algorithm, we have i = / , or b,+l has slope greater than or equal to the slope of hi, or else lines (b,, b,+l, hi) trap a lattice point to the left of x = x0, the solution line for constraints (b,, hi). (The latter is true because some lower bounding line rejects the solution set S, which implies that b, intersects bz+~ to the left of x0.) In any case, b, cannot be discarded without changing the solution line. The solution line for the problem must be obtainable from b,, some bk that rejects S, and one upper bounding line. We claim that it suffices to use h~ as the upper bounding line. For, from the discard algorithm, either j = J, or slope hj+~ is less than or equal to slope bz, or else at least one lattice point is trapped in the region bounded by (bi, hj, hj+O. In the first two instances, h~ is an essential binding constraint, and cannot be dropped without changing the solution. In the last case, if a trapped lattice point lies to the left o f x --- x*, then hj cannot be discarded without changing the solution. If all trapped points lie on or to the right of x = x*, then hj is the upper boundary line for x = x*, and all of the other upper lines can be ignored (see Figure 7). It remains to choose a third constraint from among the lower bounding lines that reject S. This is done by using a discard algorithm, with objective function parallel to b,. Specifically, let (Ii, ..., Ir) be the lower bounding lines which reject
108
S1DNIE DRESHER FEIT
S. We try to find one of these lines which blocks a lattice point not blocked by the others.
Algorithm DISCARD2 (1) k = I. (2) If k = r, then STOP. (3) If there is no lattice point below 1~in the region bounded by b,, lk, and lk+l, then discard constraint h, set k = k + 1, and go to (2). Else STOP. We claim that it suffices to solve the standard three-constraint problem (hj, b,, lk). The argument showing that/k is an essential constraint is similar to the argument given for hi. From DISCARD2, either k --- r, or else at least one lattice point is trapped in the region bounded by (b,, lk, l~+~). In the first case, lk is an essential binding constraint, and cannot be dropped without changing the solution line. In the second case, if a trapped lattice point lies to the left of x = x*, lk cannot be discarded without changing the solution line. If all trapped points lie on or to the right o f x = x*, then lk is the lower boundary line for x = x*, and (lk+t, . . . , lr) can be discarded without changing the solution. Finally, we use the three-constraint algorithm to solve (b,, lk, h~). If these constraints bound a lattice-free region, then no solution exists for the original problem. If the solution set lies on a line x = n, the solution set for R is contained in the intersection of this line with R. If no lattice points lie on the intersection segment, then R contains no lattice points.
5. The Double Greatest Common Divisor (DGCD) This section and the next contain the technical details for carrying out the computations that lie at the core of the algorithm. The technique for solving two-constraint problems will consist of manipulating the positions of the constraint lines in the plane, using lattice-preserving maps, until a position is reached for which a solution is immediately visible. The only operation used will be the addition of a positive multiple of one column of an array to another. Thus the signs of all subdeterminants will be preserved throughout the algorithm. This operation will be repeated at most O(L) times. An array
A =
(alla2) a2t
a22
a31
a32
(e e)
with sign pattern
+
+
()
and with -a~ l/al2 > -a2t/a22 > -a3t/a32 is to be carried to array
E =
e2~ e22
e31
e32
with sign pattern
+
+
..~
where eli --- el2 and e3t -< e32. The number of steps required for the simple algorithm presented here is bounded by the number of steps required to find the greatest common divisor of ( a l t, at2). The total effect of the algorithm is to multiply A by a 2 x 2 array whose entries are bounded in size by max{a,j].
109
Two-Variable Integer Programming Problem
The algorithm requires that we assign a sign to 0. We will always regard 0 as having a "--" sign. Let us concentrate first on rows 2 and 3. The algorithm below will carry these rows to an array such that sign p a t t e r n = ( +
+),
column s u m = (++).
To simplfy notation, write rows 2 and 3 as
(s) rl
-r2
s2
DGCD Algorithm (1) Let T = I.s~/rd. Set column 2 = column 2 + T x column 1. (2) Let T = tr2/s21. If T = 0, STOP. Elge set column I = column I + T x column 2. (3) Let T : ts,/rd. If T = 0, STOP. Else set column 2 = column 2 + T x column 1, and go to (2)• After this algorithm ends, replace column 1 with column 1 + column 2 to get the desired sign pattern in rows 2 and 3. What pattern will row 1 have now? Because of the sign values of the suNleterminants and the fact that the last replacement occurred in column 1, the pattern is either ( - - ) with e,, __ e,2 and e3, -< e32 and we are done, or else it is still ( - +). In the latter case, we use an algorithm almost identical to the [N3CD. Represent rows 1 and 2 as
r2
-s2
In the earlier case, only negative entries could become 0, which caused no trouble. Now we prevent a " + " element from becoming 0 by defining [IxJJ = greatest integer strictly less than x. Replacing t / by [l J J in the D G C D produces an algorithm that will transform rows 1 and 2, preserving the sign pattern, so that the column sum of the final result is negative. Now replace column 2 by column 1 + column 2, and the desired sign pattern and order relations will be obtained. 6. The Parallel Lines Algorithm
The main result of this section is that a standard problem bounded by a vertical line v and a pair of parallel lines (see Figure 8) ex - f y = c,,
ex - f y = c2,
e, f > 0
can be solved in at most 1 + log2e + logzfsteps. Since v has the form x = d, v can be replaced by x = [d] without changing the solution. The point p in the figure is the first lattice point on v below h'. It is convenient to perform a translation of the coordinates so that p is the origin. If (0, 1) is on or below h, we are done. If not, perform a greatest c o m m o n divisor algorithm on (e, f ) to obtain r and s such that er-fs=g,
__f, e Irl < g ]sl Slope h (Figure 14). Solve the problem with objective function h and two constraints v and bl. If there is no lattice point on or below h, then there are no lattice points in the region. If q is a solution point for this
Two-Variable Integer Programming Problem
113
h
FIG. 12. Case I: slope b~ < slope h.
*q
FIG. 13
Case II' slope bt =slopeh.
I
v
~
h
h'
/ I
FIG. 14. Case 111. slope b~ > s l o p e h
problem, let h' be the line through q parallel to h and solve the parallel lines problem bounded above by h and below by h'. Since q is feasible, there will be a solution and it will satisfy constraint b, and solve the three-constraint problem. ACKNOWLEDGMENTS. I wish to thank Dan Gusfield for many informative conversations, and the referee for several helpful comments and suggestions. REFERENCES (Note:
References [2] and [41 are not cited m the text.)
1. BELt, D.E. A theorem concerning the integer lattice. Stud Appl. Math. 56 (1977), 187-188. 2. HIRSCHBERG, D.S., AND WONG, C.K A polynomml-Ume algorithm for the knapsack problem with two variables. J ACM 23, 1 (Jan. 1976), 147-154. 3. KANNAN,R. A polynomial algorithm for the two-variable integer programming problem. J. ACM 27, 1 (Jan. 1980), !18-122. 4. LENSTRA,H W Integer programming wah a fixed number of variables. Rep. 81-03, Mathematisch lnsutuut, Amsterdam. 5. SCAR~,H.E. ProducUon sets with indlvlsibiliues. Econometrlca (1981). RECEIVED APRIL 1981; REVISEDMARCH 1983; ACCEPTEDMARCH 1983 Journal of the Assocmtlon for Computing
Maehlnery,Vo! 31, No. 1, Janlaartt 1984
