user_guide:tutorials:optimization

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Both sides previous revision Previous revision
Next revision		 Previous revision
Last revisionBoth sides next revision
tutorial:optimization [2013/08/22 19:56] – [Chvátal-Gomory Closure - Example 1] pfetschuser_guide:tutorials:optimization [2019/01/25 10:56] – ↷ Links adapted because of a move operation oroehrig
Line 11: Line 11:
  
 There are several other tutorials that cover similar topics: There are several other tutorials that cover similar topics:
-  * [[tutorial:ilp_and_hilbertbases|ILP and Hilbert Bases]] +  * [[user_guide:tutorials:ilp_and_hilbertbases|ILP and Hilbert Bases]] 
-  * [[tutorial:michaels_tutorial2|Gomory Cuts]] +  * [[user_guide:tutorials:michaels_tutorial2|Gomory Cuts]] 
-  * [[tutorial:michaels_tutorial|Branch and Bound]] +  * [[user_guide:tutorials:michaels_tutorial|Branch and Bound]] 
-  * [[tutorial:lattice_polytopes_tutorial|Tutorial for Lattice Polytopes]]+  * [[user_guide:tutorials:lattice_polytopes_tutorial|Tutorial for Lattice Polytopes]]
  
 This tutorial is targeted towards the optimization community, since, surprisingly, polymake does not seem to be well known in this community. In particular, the community still tends to use the quite old program ''porta'' to compute convex hulls or inequality descriptions. While ''porta'' still does a decent job here, ''polymake'' offers a much broader feature set. Polymake supports several convex hull algorithms which might be better suited depending on the data. Moreover it offers many visualization tools that can help to better //understand// a given polytope. We think that polymake has many advantages for discrete optimizers and hope that this tutorial will help to spread the usage of polymake. This tutorial is targeted towards the optimization community, since, surprisingly, polymake does not seem to be well known in this community. In particular, the community still tends to use the quite old program ''porta'' to compute convex hulls or inequality descriptions. While ''porta'' still does a decent job here, ''polymake'' offers a much broader feature set. Polymake supports several convex hull algorithms which might be better suited depending on the data. Moreover it offers many visualization tools that can help to better //understand// a given polytope. We think that polymake has many advantages for discrete optimizers and hope that this tutorial will help to spread the usage of polymake.
  
 +You can find files of the example LPs in the folder ''demo/optimization_tut'' in your ''polymake'' directory.
  
 ===== Input: lp2poly ===== ===== Input: lp2poly =====
Line 35: Line 36:
 </code> </code>
  
-Thus, the file describes a 0/1-cube in three dimensions. It should be easy to adapt this format to other cases (If for example ''x1'' does not have any bounds you can write ''x1 free'' instead). +Thus, the file describes a 0/1-cube in three dimensions. It should be easy to adapt this format to other cases (If for example ''x1'' does not have any bounds you can write ''x1 free'' instead). We create a polytope from the file via:
- +
-Now assume that this example is contained in file ''c3t.lp''. We create a polytope from the file via:+
 <code> <code>
-polytope > $f=lp2poly('c3t.lp');+polytope > $f=lp2poly('demo/optimization_tut/c3t.lp');
 </code> </code>
 The polytope ''$f'' is coded via floating point numbers: The polytope ''$f'' is coded via floating point numbers:
Line 76: Line 75:
 1 0 1 1 1 0 1 1
 </code> </code>
-This vertex corresponds to setting ''x1=0, x2=1, x3=3''. The optimal face can also be computed:+This vertex corresponds to setting ''x1=0, x2=1, x3=1''. The optimal face can also be computed:
 <code> <code>
 polytope > print $p->LP->MAXIMAL_FACE; polytope > print $p->LP->MAXIMAL_FACE;
Line 89: Line 88:
 polytope > $p->VISUAL->DIRECTED_GRAPH; polytope > $p->VISUAL->DIRECTED_GRAPH;
 </code> </code>
-{{ :tutorial:c3t_graph.gif?300 }}+{{ user_guide:c3t_graph.gif?300 }}
  
 The minimal and maximal faces can be visualized via The minimal and maximal faces can be visualized via
Line 95: Line 94:
 polytope > $p->VISUAL->MIN_MAX_FACE; polytope > $p->VISUAL->MIN_MAX_FACE;
 </code> </code>
-{{ :tutorial:c3t_maxface.gif?300 |}}+{{ user_guide:c3t_maxface.gif?300 |}}
  
  
Line 108: Line 107:
 === Bounded Polyhedra === === Bounded Polyhedra ===
  
-Let us illustrate the approach via the example of the //stable set problem//: Here one is given an (undirected) Graph G = (V,E) with node set V and edges E. The goal is to find a largest subset of node V' such that any two nodes in V' are not connected by an edge.+Let us illustrate the approach via the example of the //stable set problem//: Here one is given an (undirected) Graph G = (V,E) with node set V and edges E. The goal is to find a largest subset of nodes V' such that any two nodes in V' are not connected by an edge.
  
 For our example consider the 5-cycle, i.e., the graph C<sub>5</sub> with five nodes {1, 2, 3, 4, 5} and edges {1,2}, {2,3}, {3,4}, {4,5}, {5,1}. A formulation of the stable set problem for this graph looks as follows: For our example consider the 5-cycle, i.e., the graph C<sub>5</sub> with five nodes {1, 2, 3, 4, 5} and edges {1,2}, {2,3}, {3,4}, {4,5}, {5,1}. A formulation of the stable set problem for this graph looks as follows:
Line 134: Line 133:
 We assume that the above information is contained in the file ''stab.lp''. We now read it into polymake and convert it to rational form, as explained above: We assume that the above information is contained in the file ''stab.lp''. We now read it into polymake and convert it to rational form, as explained above:
 <code> <code>
-polytope > $f=lp2poly('stab.lp');+polytope > $f=lp2poly('demo/optimization_tut/stab.lp');
 polytope > $p = new Polytope<Rational>($f); polytope > $p = new Polytope<Rational>($f);
 </code> </code>
Line 194: Line 193:
 End End
 </code> </code>
-{{ :tutorial:ip-unb.gif?300 | Picture of unbounded polyhedron (truncated at upper right)}}+{{ user_guide:ip-unb.gif?300 | Picture of unbounded polyhedron (truncated at upper right)}}
  
 We now assume that the example is contained in the file ''unbounded.lp'' and proceed as above We now assume that the example is contained in the file ''unbounded.lp'' and proceed as above
 <code> <code>
-polytope > $f = lp2poly('unbounded.lp');+polytope > $f = lp2poly('demo/optimization_tut/unbounded.lp');
 polytope > $pin = new Polytope<Rational>($f); polytope > $pin = new Polytope<Rational>($f);
 </code> </code>
Line 263: Line 262:
 </code> </code>
  
-{{ :tutorial:ip-unb-integerhull.gif?300 |}}+{{ user_guide:ip-unb-integerhull.gif?300 |}}
  
 Note that the upper right part (including the red vertices) arises from truncation of the polyhedron for visualization. Note that the upper right part (including the red vertices) arises from truncation of the polyhedron for visualization.
Line 292: Line 291:
 In this example there are two integral variables ''x1'' and ''x2'', while ''s1'' and ''s2'' are continuous variables. Assuming the data is contained in the file ''mip.lp'', we proceed as follows: In this example there are two integral variables ''x1'' and ''x2'', while ''s1'' and ''s2'' are continuous variables. Assuming the data is contained in the file ''mip.lp'', we proceed as follows:
 <code> <code>
-polytope > $m=lp2poly('mip.lp');+polytope > $m=lp2poly('demo/optimization_tut/mip.lp');
 polytope > $p = new Polytope<Rational>($m); polytope > $p = new Polytope<Rational>($m);
 </code> </code>
Line 312: Line 311:
 ===== Integral Polytopes and Total Unimodularity ===== ===== Integral Polytopes and Total Unimodularity =====
  
-As explained in the previous example, the integral points in a polytope are of particular interest in discrete optimization. These points are called //lattice points// in polymake and the corresponding convex hull //lattice polytope//. The handling of such polytopes is explained in more detail in the [[tutorial:lattice_polytopes_tutorial|Tutorial for Lattice Polytopes]].+As explained in the previous example, the integral points in a polytope are of particular interest in discrete optimization. These points are called //lattice points// in polymake and the corresponding convex hull //lattice polytope//. The handling of such polytopes is explained in more detail in the [[user_guide:tutorials:lattice_polytopes_tutorial|Tutorial for Lattice Polytopes]].
  
-Of particular interest for discrete optimization are properties of the orginal inequality system to define a lattice polytope, i.e., a polytope such that all of its vertices are integral. One particularly interesting case occurs if the matrix defining the polytope is //totally unimodular// and the right hand side is integral.+Of particular interest for discrete optimization are properties of the orginal inequality system to define a lattice polytope, i.e., a polytope such that all of its vertices are integral (this can be tested by checking the property ''LATTICE''). One particularly interesting case occurs if the matrix defining the polytope is //totally unimodular// and the right hand side is integral.
  
 Using the polymake extension [[https://github.com/xammy/unimodularity-test/wiki/Polymake-Extension|Unimodularity]] by Matthias Walter, this can be checked as illustrated in the following examples. Using the polymake extension [[https://github.com/xammy/unimodularity-test/wiki/Polymake-Extension|Unimodularity]] by Matthias Walter, this can be checked as illustrated in the following examples.
Line 335: Line 334:
 In the second example, we reuse the file ''c3t'' from the example above. We read it into polymake: In the second example, we reuse the file ''c3t'' from the example above. We read it into polymake:
 <code> <code>
-polytope > $f=lp2poly('c3t.lp');+polytope > $f=lp2poly('demo/optimization_tut/c3t.lp');
 polytope > $p = new Polytope<Rational>($f); polytope > $p = new Polytope<Rational>($f);
 </code> </code>
Line 354: Line 353:
   * The function ''make_totally_dual_integral'' takes a polytope and returns a new polytope with inequalities that are TDI.   * The function ''make_totally_dual_integral'' takes a polytope and returns a new polytope with inequalities that are TDI.
  
-Note that the input has to be full-dimensional in order to use these functions (this can be handled via ''ambient_lattice_normalization'').+Note that the input has to be full-dimensional in order to use these functions.
  
 To demonstrate the behavior of these functions, consider the 5-cycle example from above again: To demonstrate the behavior of these functions, consider the 5-cycle example from above again:
Line 379: Line 378:
 Let us test whether the inequality system of this example is TDI. Thus, we first load the data as usual: Let us test whether the inequality system of this example is TDI. Thus, we first load the data as usual:
 <code> <code>
-polytope > $f = lp2poly('stab.lp');+polytope > $f = lp2poly('demo/optimization_tut/stab.lp');
 polytope > $p = new Polytope<Rational>($f); polytope > $p = new Polytope<Rational>($f);
 </code> </code>
 We now extract the corresponding inequality system and check it for TDIness: We now extract the corresponding inequality system and check it for TDIness:
 <code> <code>
-polytope> $M = new Matrix<Rational>($p->INEQUALITIES);+polytope > $M = new Matrix<Rational>($p->INEQUALITIES);
 polytope > totally_dual_integral($M); polytope > totally_dual_integral($M);
-0+       
 </code> </code>
-Hence, the system is not TDI, which we expected from general theory, since we know that the polytope is not integral, but the system is integral. Consequently, let us construct a TDI-system for this polytope:+Since the last line is empty, the system is not TDI, which we expected from general theory, since we know that the polytope is not integral, but the system has integral coefficients. Consequently, let us construct a TDI-system for this polytope:
 <code> <code>
 polytope > $q = make_totally_dual_integral($p); polytope > $q = make_totally_dual_integral($p);
Line 406: Line 405:
 11: 0 >= -1 11: 0 >= -1
 </code> </code>
-As expected, the right hand side is non integral (otherwise, we know from general theory that the polytope would be integral as well).+As expected, the right hand side is non integral (otherwise, we know from general theory that the polytope would be integral as well). The result is now TDI: 
 +<code> 
 +polytope > $N = new Matrix<Rational>($q->INEQUALITIES); 
 +polytope > print totally_dual_integral($N); 
 +1 
 +</code> 
 +*Note* that we need to take the inequalities instead of facets here, since facets are irredundant and thus might not be TDI, although the complete set of inequalities is TDI.
  
-===== 0/1-Polytopes ===== 
  
 ===== Chvátal-Gomory Closure ===== ===== Chvátal-Gomory Closure =====
Line 438: Line 442:
 As before we read in the file using ''lp2poly'': As before we read in the file using ''lp2poly'':
 <code> <code>
-polytope > $f = lp2poly('stab.lp');+polytope > $f = lp2poly('demo/optimization_tut/stab.lp');
 polytope > $p = new Polytope<Rational>($f); polytope > $p = new Polytope<Rational>($f);
 </code> </code>
  
-The Chvátal-Gomory closure of a polytope can be computed with the function ''gc_closure''. The function takes a full-dimensional polytope and returns a new polytope. This contains the system of inequlities defining the closure in the property ''INEQUALITIES''. For our example, we obtain:+The Chvátal-Gomory closure of a polytope can be computed with the function ''gc_closure''. The function takes a full-dimensional polytope and returns a new polytope. This contains the system of inequalities defining the closure in the property ''INEQUALITIES''. For our example, we obtain:
 <code> <code>
 polytope > $g = gc_closure($p); polytope > $g = gc_closure($p);
Line 460: Line 464:
 11: 0 >= -1 11: 0 >= -1
 </code> </code>
-Let us check whether the resulting polytope is integral by examining its vertices:+Let us check whether the resulting polytope is integral:
 <code> <code>
-polytope > print $g->VERTICES; +polytope > print $g->LATTICE; 
-1 0 0 0 1 0 +1
-1 0 0 0 0 1 +
-1 0 1 0 0 0 +
-1 0 0 0 0 0 +
-1 1 0 0 0 0 +
-1 0 0 1 0 0 +
-1 1 0 1 0 0 +
-1 0 1 0 0 1 +
-1 1 0 0 1 0 +
-1 0 1 0 1 0 +
-1 0 0 1 0 1+
 </code> </code>
 Thus, in this case, we have obtained the integer hull by one step of the Chvatal-Gomory-closure. Thus, in this case, we have obtained the integer hull by one step of the Chvatal-Gomory-closure.
Line 479: Line 473:
 ==== Chvátal-Gomory Closure - Example 2 ==== ==== Chvátal-Gomory Closure - Example 2 ====
  
-Let us now consider the classical example of a polytope with the vertices of simplex in d dimensions and the point 1/2 times (1, ..., 1). It can be shown that such a polytope has rank at least log(d) - 1, see [Pokutta, 2011]. In our example, we use d = 4:+Let us now consider the classical example of a polytope with the vertices of simplex in d dimensions and the point 1/2 times (1, ..., 1). It can be shown that such a polytope has rank at least log(d) - 1, see [[http://www.box.net/shared/at1y8i3pq434bxt6m9xm|Pokutta, 2011]]]. In our example, we use d = 4:
 <code> <code>
 polytope > $M = new Matrix<Rational>([[1,0,0,0,0],[1,1,0,0,0],[1,0,1,0,0],[1,0,0,1,0],[1,0,0,0,1],[1,1/2,1/2,1/2,1/2]]); polytope > $M = new Matrix<Rational>([[1,0,0,0,0],[1,1,0,0,0],[1,0,1,0,0],[1,0,0,1,0],[1,0,0,0,1],[1,1/2,1/2,1/2,1/2]]);
Line 495: Line 489:
 6: -x1 - x3 - x4 >= -1 6: -x1 - x3 - x4 >= -1
 7: -x2 - x3 - x4 >= -1 7: -x2 - x3 - x4 >= -1
 +polytope > print $t1->LATTICE;
 +0
 +</code>
 +Thus, one round was not enough to produce an integral polytope. Indeed, the vertices are
 +<code>
 polytope > $t1->VERTICES; polytope > $t1->VERTICES;
 polytope > print $t1->VERTICES; polytope > print $t1->VERTICES;
Line 504: Line 503:
 1 0 0 0 1 1 0 0 0 1
 </code> </code>
-Thus, one round was not enough to produce an integral polytope. However, one more round is enough:+ 
 +However, one more round is enough:
 <code> <code>
 polytope > $t2 = gc_closure($t1); polytope > $t2 = gc_closure($t1);
Line 515: Line 515:
 3: x1 >= 0 3: x1 >= 0
 4: -x1 - x2 - x3 - x4 >= -1 4: -x1 - x2 - x3 - x4 >= -1
 +polytope > print $t2->LATTICE;
 +1
 </code> </code>
  
 ===== Lift-and-project closure ===== ===== Lift-and-project closure =====
  
-===== Other Software =====+The lift-and-project closure of a 0/1-polytope P can be generated as follows: for each coordinate compute the intersection of P with the pair of opposite cube faces and compute the convex hull. Then intersect the result with P. The following subroutine performs this operation - the code is somewhat complicated throught the fact that we need to check whether parts are empty. 
 +<code> 
 +sub lpclosure 
 +{ 
 +   my $p = shift; 
 +   my $d = $p->AMBIENT_DIM; 
 +   my $q = new Polytope<Rational>($p); 
 +   for (my $k = 0; $k < $d; $k = $k+1) 
 +   { 
 +      if ( $q->DIM == -1 )         # can stop as soon as $q is empty 
 +      { 
 + return $q; 
 +      } 
 + 
 +      # create reversed opposite inequalities of 0/1-cube and corresponding polyhedra 
 +      my $v1 = new Vector<Rational>(0 | -unit_vector($d, $k)); 
 +      my $v2 = new Vector<Rational>(-1 | unit_vector($d, $k)); 
 + 
 +      # create intersection of corresponding polyhedra with iterated polyhedron $q 
 +      my $b1 = new Polytope<Rational>(INEQUALITIES => $v1 / $q->FACETS); 
 +      my $b2 = new Polytope<Rational>(INEQUALITIES => $v2 / $q->FACETS); 
 + 
 +      if ( ($b1->DIM > -1) && ($b2->DIM > -1) ) 
 +      { 
 + my $c = conv($b1, $b2); 
 + $q = intersection($q, $c); 
 +      } 
 +      elsif ( ($b1->DIM > -1) && ($b2->DIM == -1) ) 
 +      { 
 + $q = intersection($q, $b1); 
 +      } 
 +      elsif ( ($b1->DIM == -1) && ($b2->DIM > -1) ) 
 +      { 
 + $q = intersection($q, $b2); 
 +      } 
 +   } 
 +   return $q; 
 +} 
 +</code> 
 + 
 +==== Lift-and-Projet Closure - Example 1 ==== 
 + 
 +For our well known stable set example, we get the following: 
 +<code> 
 +polytope > $q = lpclosure($p); 
 +polytope > $q->FACETS; 
 +polytope > print_constraints($q); 
 +Facets: 
 +0: -x4 - x5 >= -1 
 +1: -x2 - x3 >= -1 
 +2: -x1 - x2 - x3 - x4 - x5 >= -2 
 +3: -x1 - x5 >= -1 
 +4: -x3 - x4 >= -1 
 +5: x2 >= 0 
 +6: x1 >= 0 
 +7: x3 >= 0 
 +8: -x1 - x2 >= -1 
 +9: x4 >= 0 
 +10: x5 >= 0 
 +</code> 
 +Thus, the lift-and-project closure in this case gives the integral hull (as we have seen above). 
 + 
 +==== Lift-and-Projet Closure - Example 2 ==== 
 + 
 +Let us now consider the same example as for CG-closures: 
 +<code> 
 +polytope > $M = new Matrix<Rational>([[1,0,0,0,0],[1,1,0,0,0],[1,0,1,0,0],[1,0,0,1,0],[1,0,0,0,1],[1,1/2,1/2,1/2,1/2]]); 
 +polytope > $p = new Polytope<Rational>(POINTS => $M); 
 +polytope > $q = lpclosure($p); 
 +polytope > $q->FACETS; 
 +polytope > print_constraints($q); 
 +Facets: 
 +0: x2 >= 0 
 +1: x4 >= 0 
 +2: -x1 - x2 - x3 - x4 >= -1 
 +3: x3 >= 0 
 +4: x1 >= 0 
 +</code> 
 +Thus, we have obtained the integral hull in a single step of the lift-and-project closure as opposed to two steps in the CG-closure. 
  
  
  
  
  • user_guide/tutorials/optimization.txt
  • Last modified: 2019/02/04 22:55
  • by 127.0.0.1