Gomory Cuts

In this Tutorial, I describe how to solve an ILP using the gomory cut approach implemented in polymake.

This preamble must be integrated in every script file, so that we can use the polytope functionalities of polymake.

> use application "polytope";

I introduce a subroutine, that returns the non-integer part of some rational number.

> sub f {
>         my $x=shift;
>         my $denom = denominator($x);
>         my $nom = convert_to<Integer>($denom * $x);
>         my $remainder = $nom % $denom;
>         if($remainder < 0){
>           $remainder += $denom;
>         }
>         return new Rational($remainder, $denom);
> }

The main Routine of this script gets invoked with the following parameters:

• The Matrix of equalities as an Matrix
• The Matrix of inequalites as an Matrix
• The Objective function as a list of Rationals.

> sub gomory {
>     
>         my $A = shift;
>         my $ineqs = shift;
>         my @obj_list = @_;
>         my $obj = new Vector<Rational>(@obj_list);
>            
>         my $p = new Polytope<Rational>;
>         $p->EQUATIONS = $A;
>         $p->INEQUALITIES = $ineqs;
> 
> 
> ### Feasibility check
> # Since all the inequalities that we add are valid for the integer hull 
> # of the polytope, we can conclude that the integer hull is empty,
> # when the examined polytope in any iteration cycle is empty.
> 
>         if(!$p->FEASIBLE){
>             print "The integer hull of P is empty.\n";
>             return;
>         }
> 
> 
> ### Solving the relaxation
> 
>         $p->LP = new LinearProgram<Rational>(LINEAR_OBJECTIVE=>$obj);
>     
>         my $x = $p->LP->MAXIMAL_VERTEX;
>         print "Optimal solution of the relaxation: ", $x,"\n";
>     
>     
>         my @B = ();
>         my @N = ();
>         for(my $k = 1; $k < $x->dim(); $k++){
>             if($x->[$k] > 0){
>                 push @B, $k;
>             }else{
>                 push @N, $k;
>             }
>         }
> 
>    
> ### Check for integer solution
> # If we found an integer solution in the linear relaxation,
> # it's the solution of the ILP aswell.
> # If not, we remember the coordinate that is not integer for later steps.
> 
>         my $i = -1;
>         for(my $k=0; $k<@B; $k++){
>             if(f($x->[$B[$k]]) != 0){
>                 $i = $k;
>             }
>         }
>         
>         if($i == -1){
>             print "Found optimal integer solution: ", $x, " -> ", $p->LP->MAXIMAL_VALUE,"\n";
>             return;
>         }
> 
> 
> ### Compute optimal basis
> # In this example, I only implemented the computation of a Basis
> # in the non-degenerated case (that is, where every Basis Variables are strictly positive.
> 
>         my @A_b = ();
>         foreach my $b (@B){
>             push @A_b, $A->col($b);
>         }
>         
>         my $A_b = transpose(new Matrix<Rational>(@A_b));
>         my $A_inv = inv($A_b);
>             
> 
>     
> ### Compute new equality
> # Here is where the remembered non-integer variable comes into play,
> # to compute the Gomory cut itself.
>     
>         my @new_row;
>         my @new_col;
>         
>         for(my $j=0; $j<$A->rows(); $j++){
>             $new_col[$j] = 0;
>         }
>         
>         $new_row[0] = f($A_inv->row($i)*($A->col(0)));
>         foreach my $j (@N){
>             $new_row[$j] = (-1)*f((-1)*$A_inv->row($i)*$A->col($j));
>         }
>         foreach my $j (@B){
>             $new_row[$j] = 0;
>         }
>         push @new_row, 1;
>         my $new_row_vec = new Vector<Rational>(@new_row);
>         my $new_col_vec = new Vector<Rational>(@new_col);
>     
>         print "New equality: $new_row_vec\n";      
> 
>     
> ### Compute new inequality
> # The new inequality demands, that the newly introduced slack variable will be positive.  
>         
>         my @new_ineq_col;    
>         for(my $j=0; $j<$ineqs->rows(); $j++){
>             $new_ineq_col[$j] = 0;
>         }
>         
>         my @new_ineq_row = @new_ineq_col;
>         #push @new_ineq_row, 0; #this leads to dimension mismatches
>         push @new_ineq_row, 1;
>         
>         my $new_ineq_row_vec = new Vector<Rational>(@new_ineq_row);
>         my $new_ineq_col_vec = new Vector<Rational>(@new_ineq_col);    
> 
>     
> ### Construct the new matrices
> # The two matrices constructed are the equality matrix A and the inequality matrix
> # that restricts the variables to be positive.
> # In both cases we add a new column on the right side and a new row at the bottom.
>         
>         my $new_A = new Matrix<Rational>(($A|$new_col_vec) / $new_row_vec); 
>         my $new_ineqs = new Matrix<Rational>(($ineqs|$new_ineq_col_vec) / $new_ineq_row_vec);
>         
>         push @obj_list, 0;
>         
>         gomory($new_A, $new_ineqs, @obj_list);
> }

Finally we create two example matrices and an example objective function, to call the main routine.

> $A = new Matrix<Rational>([
> [-14, 7, -2, 1, 0, 0],
> [-3, 0, 1, 0, 1, 0],
> [-3, 2, -2, 0, 0, 1]]);  
> 
> $ineqs = new Matrix<Rational>([
> [0, 1, 0, 0, 0, 0],
> [0, 0, 1, 0, 0, 0],
> [0, 0, 0, 1, 0, 0],
> [0, 0, 0, 0, 1, 0],
> [0, 0, 0, 0, 0, 1]]);
>     
> @objective = (0,4,-1,0,0,0);
> 
> gomory($A, $ineqs, @objective);
Optimal solution of the relaxation: 1 20/7 3 0 0 23/7
New equality: 5/7 0 0 -2/7 -4/7 0 1
Optimal solution of the relaxation: 1 5/2 7/4 0 5/4 3/2 0
New equality: 1/2 0 0 0 0 0 -1/2 1
Optimal solution of the relaxation: 1 2 1/2 1 5/2 0 1 0
New equality: 1/2 0 0 0 0 -1/2 0 0 1
Optimal solution of the relaxation: 1 2 1 2 2 1 1 0 0
Found optimal integer solution: 1 2 1 2 2 1 1 0 0 -> 7