INTEGER PROGRAMMING BY THE GROEBNER BASIS METHOD
                                      
                                 
                          by Devendra Kapadia 
                                               Wolfram Research, Inc.

Our aim in this notebook is to define a Mathematica function for solving 
problems in  Integer Programming by the Groebner Basis method [Conti and 
Traverso, AAECC`9, Lecture Notes in Computer Science, Vol. 539, 130-139, 
Springer Verlag, 1991].
      
The general Integer Programming problem can  
be reduced  to that of minimizing a cost vector c.w subject to the conditions 
A.w = b, where the components of the optimal solution w are required  to be 
non-negative integers.
 
The algorithm is divided into three steps:
       
 1. The coefficients of the matrix A  are used  to form the exponents in a 
certain
             binomial ideal.
             
        2. The Groebner 
basis of the binomial ideal with respect to a monomial order
            
given by the vector c is found using the kernel function GroebnerBasis.
     
       
        3. The optimal solution is found by reading off the 
exponents in the answer
            obtained by reducing a certain monomial 
formed from the components of b,
            with respect to the Groebner 
basis and the same monomial order. The kernel
            function 
PolynomialReduce is used for this step.
   
  In contrast to Linear 
Programming, the problem of Integer Programming is  NP-complete. Hence, there 
will be an "expensive step" in any algorithm for the problem. In our case, 
the expensive step is that of finding the Groebner  
  basis.  
 
 We will 
define the function IntegerProgramming below and give several examples of 
problems which can be solved using it. At present, we assume that either all 
the components of A and b are non-negative or that the same condition holds 
for the components of c,  although this restriction may be easy to remove in 
future.  Also, the current version of Mathematica (Version 4.1) is unable to handle weight matrices with rational entries 
correctly and hence the cost matrix is converted into one with integer 
entries in our program.

IntegerProgramming[c_, A_, b_] := 

      Module[{m, n, p, q, CostOrder, CostMatrix, e, t}, 

        p = Dimensions[A][[1]]; q = Dimensions[A][[2]]; 

        CostMatrix = 

          Join[
Table[

              KroneckerDelta[i, j], {i, 1, p + 1}, {j, 1, p + q + 1}], {Join[

                Table[0, {i, 1, p + 1}], 

                If[Apply[And, 

                    Table[c[[j]] ≥ 0, {j, 1, q}]], 

                  Table[c[[j]], {j, 1, q}], 

                  Table[c[[j]] + 

                      Max[Table[-c[[j]]/(Sum[

                                  A[[i, j]], {i, 1, p}]), {j, 1, q}]]*

                        Sum[A[[i, j]], {i, 1, p}], {j, 1, q}]]]}, 

            Table[KroneckerDelta[i, p + q + 4 - j], {i, 3, q + 1}, {j, 1, 

                p + q + 1}]]; 

        CostOrder = Map[LCM[Denominator[#]]*#  &,  CostMatrix]; 

        e = Table[

            Max[0, Max[Table[-A[[i, j]], {i, 1, p}]]], {j, 1, q}]; 

        m = GroebnerBasis[

            Join[{t*Product[z[i], {i, 1, p}] - 1}, 

              Table[t^(e[[j]])*

                    Product[

                      z[i]^(A[[i, j]] + e[[j]]), {i, 1, p}] - 

                  w[j], {j, 1, q}]], 

            Join[{t}, Table[z[i], {i, 1, p}], Table[w[j], {j, 1, q}]], 

            MonomialOrder -> CostOrder]; 

        n = PolynomialReduce[

              t^Max[{0, Max[Table[-b[[i]], {i, 1, p}]]}]*

                Product[

                  z[i]^(b[[i]] + (Max[0, 

                            Max[Table[-b[[i]], {i, 1, 

                                  p}]]])), {i, 1, p}], m, 

              Join[{t}, Table[z[i], {i, 1, p}], Table[w[j], {j, 1, q}]], 

              MonomialOrder -> CostOrder][[2]]; {Sum[

            c[[If [Length[n[[i]]] == 1, n[[i, 1]], 

                    n[[i, 1, 1]]]]]*(If[Length[n[[i]]] == 1, 1, 

                  n[[i, 2]]]), {i, 1, Length[n]}], 

          Table[If[Length[n[[i]]] == 1, n[[i]] -> 1, 

              n[[i, 1]] -> n[[i, 2]]], {i, 1, 

              Length[n]}]}]

In the examples given below, we will simply write the matrices required for 
applying the IntegerProgramming function. The matrices are obtained after 
adding slack variables and/or changing signs in the cost vector  to convert 
the problem from one of maximizing the cost  to one of minimizing the cost.

Example 1:  This  is taken from pg. 372 of the book "Using Algebraic Geometry" by 
Cox, Little and O'Shea,which is the basic reference for our program :

A1 = {{3, -2, 1, 0}, {4, 1, -1, -1}};

b1 = {-1, 5};

c1 = {1, 1000, 1, 100};

IntegerProgramming[c1, A1, b1]

Example 2: This appears on pg. 374 of the above book:

A2 = {{3, 2, 1, 1, 0, 0}, {1, 2, 3, 0, 1, 0}, {2, 1, 1, 0, 0, 

            1}};

b2 = {45, 21, 18};

c2 = {-3, -4, -2, 0, 0, 0};

IntegerProgramming[c2, A2, b2]

Example 3: This is due to Ralph Gomory, who invented the Cutting Plane 
Method for solving Integer Programming problems:

A3 = {{3, 2, 0, 1, 0, 0}, {1, 4, 0, 0, 1, 0}, {3, 3, 1, 0, 0, 

            1}};

b3 = {10, 11, 13};

c3 = {-4, -5, -1, 0, 0, 0};

IntegerProgramming[c3, A3, b3]

Example 4:  I have extracted this example from Adams and Loustaunau, "An 
Introduction to Groebner Bases", pg. 109.

A4 = {{3, -2, 1, -1}, {4, 1, -1, 0}};

b4 = {-1, 5};

c4 = {1000, 1, 1, 100, 0, 0};

IntegerProgramming[c4, A4, b4]

Examples 5 and 6:  The next two are exercises from Adams and Loustaunau, pg. 112 :

c5 = {1000, 1, 1, 1, 0, 0};

IntegerProgramming[c5, A4, b4]

c6 = {1, 1000, 1, 1, 0, 0};

IntegerProgramming[c6, A4, b4]

Example 7 : The example given below is from the original paper by Conti and Traverso 
(pg. 137).  This example
                    is completed in less than 2 
seconds in a development version due to new functionality for computing
     
               Groebner bases of toric ideals added by Daniel Lichtblau.

A7 = {{2, 5, -3, 1, -2}, {1, 7, 2, 3, 

            1}, {4, -2, -1, -5, 3}};

b7 = {214, 712, 331};

c7 = {1, 1, 1, 1, 1};

Timing[IntegerProgramming[c7, A7, b7]]

Example 8: Finally, here is an instance of the knapsack problem from   
    
                                http://www.rpi.edu/~mitchj/matp6640/knapsack/.
 
                   The order of the  variables has been changed to ensure 
that the cost matrix is invertible.

A8 = {{0, 3, 5, 1, 0, 0, 1}, {0, 0, 0, 0, 1, 0, 1}, {0, 1, 0, 0, 0, 1, 

          0}, {1, 0, 1, 0, 0, 0, 0}}

b8 = {7, 1, 1, 1}

c8 = {0, -4, -6, 0, 0, 0, -3}

IntegerProgramming[c8, A8, b8]

Conclusion and Outlook:
Initially, the IntegerProgramming function relied on the use of the 
Groebner package from MathSource. I am deeply grateful to Daniel Lichtblau 
for explaining to me how the kernel functions GroebnerBasis and 
PolynomialReduce can be applied in this situation and for numerous  insights 
into the subject. 

It is very satisyfing to note that, in principle, this 
rather simple function can be used to solve any feasible Integer Programming 
problem subject to the limitations on A, b and c mentioned above. However, it 
is clear that further progress will depend upon choosing the most efficient 
monomial order for each class of problems and also on finding efficient 
methods for calculating the required Groebner bases. At present, we are 
working on refinements of the algorithm for solving special cases of the 
location problem described on 
                                            
http://mat.gsia.cmu.edu/orclass/integer/node8.html.
A truncated version of 
this problem of this type (with 8 neighborhoods) can be solved in the 
development version referred to above in about 8 secs.

Finally, it would 
be interesting to study complexity issues for Integer Programming  along the 
lines suggested by Bayer and Mumford in "What can be computed in Algebraic 
Geometry?".  Of course, this is a highly speculative suggestion !