Differences
This shows you the differences between two versions of the page.
Both sides previous revision Previous revision Next revision | Previous revision | ||
tutorial:optimization [2013/08/22 19:54] – pfetsch | user_guide:tutorials:optimization [2019/02/04 22:55] (current) – external edit 127.0.0.1 | ||
---|---|---|---|
Line 1: | Line 1: | ||
- | ====== polymake and Optimization ====== | + | {{page>.:latest:@FILEID@}} |
- | + | ||
- | By Sebastian Pokutta and Marc E. Pfetsch. | + | |
- | + | ||
- | ===== Introduction ===== | + | |
- | + | ||
- | Polymake offers many interesting features that can help (discrete) optimizers to analyze optimization problems. For example | + | |
- | * linear optimization can be performed exactly and visualized in small dimensions | + | |
- | * the convex hull of feasible points of an integer program can be computed and analyzed | + | |
- | * Hilbert bases can be computed | + | |
- | + | ||
- | There are several other tutorials that cover similar topics: | + | |
- | * [[tutorial: | + | |
- | * [[tutorial: | + | |
- | * [[tutorial: | + | |
- | * [[tutorial: | + | |
- | + | ||
- | This tutorial is targeted towards the optimization community, since, surprisingly, | + | |
- | + | ||
- | + | ||
- | ===== Input: lp2poly ===== | + | |
- | + | ||
- | The first important step is to get the desired input into '' | + | |
- | + | ||
- | < | + | |
- | Minimize | + | |
- | | + | |
- | Subject to | + | |
- | C1: x1 + x2 + x3 <= 2 | + | |
- | Bounds | + | |
- | 0 <= x1 <= 1 | + | |
- | 0 <= x2 <= 1 | + | |
- | 0 <= x3 <= 1 | + | |
- | End | + | |
- | </ | + | |
- | + | ||
- | 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 '' | + | |
- | + | ||
- | Now assume that this example is contained in file '' | + | |
- | < | + | |
- | polytope > $f=lp2poly(' | + | |
- | </ | + | |
- | The polytope '' | + | |
- | < | + | |
- | polytope > print $f-> | + | |
- | Polytope< | + | |
- | </ | + | |
- | + | ||
- | We convert it to a rational polytope via: | + | |
- | < | + | |
- | polytope > $p = new Polytope< | + | |
- | </ | + | |
- | + | ||
- | Now, '' | + | |
- | + | ||
- | + | ||
- | ===== Linear Optimization ===== | + | |
- | + | ||
- | Polymake can be used to perform several actions related to linear optimization (linear programming - LP). For instance, one can exactly solve a linear program (via lrs or cdd). Before we explain the corresponding usage, we first need to have a linear optimization problem at hand. | + | |
- | + | ||
- | Assuming that we are given the above example in variable '' | + | |
- | < | + | |
- | polytope > print $p-> | + | |
- | 0 1 1 1 | + | |
- | </ | + | |
- | Thus - as described in the file - the objective function coefficients are 1 for all three variables (and there is an offset of 0). | + | |
- | + | ||
- | Now, we can solve the corresponding linear program via | + | |
- | < | + | |
- | polytope > print $p-> | + | |
- | 2 | + | |
- | </ | + | |
- | Thus, the maximal value that we can obtain via the above linear objective function is 2. We can also get an optimal vertex via | + | |
- | < | + | |
- | polytope > print $p-> | + | |
- | 1 0 1 1 | + | |
- | </ | + | |
- | This vertex corresponds to setting '' | + | |
- | < | + | |
- | polytope > print $p-> | + | |
- | {4 5 6} | + | |
- | </ | + | |
- | This means that the optimal face is the convex hull of three vertices (with indices 4, 5, 6). | + | |
- | + | ||
- | Of course, by replacing '' | + | |
- | + | ||
- | The directed graph obtained by directing the graph of the polytope in the direction of increasing objective function can be obtained via | + | |
- | < | + | |
- | polytope > $p-> | + | |
- | </ | + | |
- | {{ : | + | |
- | + | ||
- | The minimal and maximal faces can be visualized via | + | |
- | <code> | + | |
- | polytope > $p-> | + | |
- | </ | + | |
- | {{ : | + | |
- | + | ||
- | + | ||
- | ===== Computing Facets ===== | + | |
- | + | ||
- | An important action that is often needed to come up with new facet describing inequalities for combinatorial optimization problems is the computation of convex hulls for small examples. | + | |
- | + | ||
- | ==== Pure Integer Case ==== | + | |
- | + | ||
- | We begin with the case in which all variables are required to be integral, i.e., the //pure integer case//. Moreover, the approach depends on whether the polyhedron is bounded or not. | + | |
- | + | ||
- | === 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. | + | |
- | + | ||
- | For our example consider the 5-cycle, i.e., the graph C< | + | |
- | < | + | |
- | Maximize | + | |
- | obj: x#1 + x#2 + x#3 + x#4 + x#5 | + | |
- | Subject to | + | |
- | | + | |
- | | + | |
- | | + | |
- | | + | |
- | | + | |
- | Bounds | + | |
- | 0 <= x#1 <= 1 | + | |
- | 0 <= x#2 <= 1 | + | |
- | 0 <= x#3 <= 1 | + | |
- | 0 <= x#4 <= 1 | + | |
- | 0 <= x#5 <= 1 | + | |
- | General | + | |
- | x#1 x#2 x#3 x#4 x#5 | + | |
- | End | + | |
- | </ | + | |
- | Here, '' | + | |
- | + | ||
- | We assume that the above information is contained in the file '' | + | |
- | < | + | |
- | polytope > $f=lp2poly(' | + | |
- | polytope > $p = new Polytope< | + | |
- | </ | + | |
- | + | ||
- | We are now interested in all feasible solutions to the above problem, i.e., all assignments of 0 or 1 to the variables such that the above inequalities are satisfied. These feasible points can be computed via: | + | |
- | < | + | |
- | polytope > $p-> | + | |
- | </ | + | |
- | To understand these points and make computational use of this information, | + | |
- | < | + | |
- | polytope > $s=new Polytope(POINTS=> | + | |
- | </ | + | |
- | Here, the coordinate labels, i.e., the variable names, are copied to the new polytope. | + | |
- | + | ||
- | Now, the facets of the new polytope can be computed and listed via: | + | |
- | < | + | |
- | polytope > print_constraints($s); | + | |
- | Facets: | + | |
- | 0: x#3 >= 0 | + | |
- | 1: -x#4 -x#5 >= -1 | + | |
- | 2: x#1 >= 0 | + | |
- | 3: -x#2 - x#3 >= -1 | + | |
- | 4: -x#1 - x#2 - x#3 - x#4 - x#5 >=-2 | + | |
- | 5: -x#1 - x#2 >=-1 | + | |
- | 6: x#4 >= 0 | + | |
- | 7: -x#1 - x#5 >= -1 | + | |
- | 8: -x#3 - x#4 >= -1 | + | |
- | 9: x#2 >= 0 | + | |
- | 10: x#5 >= 0 | + | |
- | </ | + | |
- | The facet defining inequalities can be interpreted as follows: | + | |
- | * There are five trivial inequalities '' | + | |
- | * The five original ' | + | |
- | * We have the so-called //odd-cycle inequality// | + | |
- | + | ||
- | Of course, one can also use the usual polymake output, e.g., '' | + | |
- | + | ||
- | This example showed one of the routine actions often performed by discrete optimizers. Of course, this action can also be performed by a script, which makes the computation a one-line command. | + | |
- | + | ||
- | Note that the size of instances that can be handled will probably be small. Usually, things become difficult from dimension 15 on, but it depends on the particular structure of your instances, i.e., on the number of facets and lattice points. | + | |
- | + | ||
- | + | ||
- | === Unbounded Polyhedra === | + | |
- | + | ||
- | If the underlying polyhedron is unbounded, the approach above does not work anymore, since there are infinitely many lattice points. Arguably, | + | |
- | + | ||
- | The following mathematical insights are important to treat the unbounded case. First, we have to assume that the data, i.e., the inequality description of the polyhedron '' | + | |
- | + | ||
- | To illustrate the construction, | + | |
- | < | + | |
- | Minimize | + | |
- | | + | |
- | Subject to | + | |
- | C1: x1 + x2 >= 0.5 | + | |
- | C2: x1 - 2 x2 <= 1.5 | + | |
- | C3: x2 - 2 x1 <= 1.5 | + | |
- | General | + | |
- | x1 x2 | + | |
- | End | + | |
- | </ | + | |
- | {{ : | + | |
- | + | ||
- | We now assume that the example is contained in the file '' | + | |
- | < | + | |
- | polytope > $f = lp2poly(' | + | |
- | polytope > $pin = new Polytope< | + | |
- | </ | + | |
- | The visualization in the picture can be generated with '' | + | |
- | + | ||
- | We now extract the rays of the recession cone | + | |
- | < | + | |
- | polytope > $rays = $pin-> | + | |
- | </ | + | |
- | This command first computes all vertices of the polyhedron (this includes unbounded vertices); note that is involves a convex hull computation. The set '' | + | |
- | < | + | |
- | polytope > print $rays; | + | |
- | 0 1 1/2 | + | |
- | 0 1 2 | + | |
- | </ | + | |
- | Thus, there are two rays that are generators of the recession cone. | + | |
- | + | ||
- | We now have to construct the Minkowski hull of all intervals '' | + | |
- | < | + | |
- | my $zero = unit_vector< | + | |
- | my $B = new Polytope< | + | |
- | foreach my $r (@$rays) | + | |
- | { | + | |
- | my $M = new Matrix< | + | |
- | | + | |
- | $M = $M / $zero; | + | |
- | my $ptemp = new Polytope< | + | |
- | $B = minkowski_sum($B, | + | |
- | } | + | |
- | </ | + | |
- | The code first generates a polytope '' | + | |
- | + | ||
- | The next step is to obtain the bounded part '' | + | |
- | < | + | |
- | polytope > $Qpoints = $pin-> | + | |
- | polytope > $Q = new Polytope< | + | |
- | </ | + | |
- | + | ||
- | The two polytopes are now combined: | + | |
- | < | + | |
- | polytope > $p = minkowski_sum($Q, | + | |
- | </ | + | |
- | + | ||
- | We now generate the lattice points (as in the bounded part) and add the rays from above: | + | |
- | < | + | |
- | polytope > $latticemat = new Matrix< | + | |
- | polytope > $newpoints = new Matrix< | + | |
- | </ | + | |
- | Here, '' | + | |
- | + | ||
- | Finally, the polytope we are interested in is: | + | |
- | < | + | |
- | polytope > $q = new Polytope(POINTS=> | + | |
- | </ | + | |
- | + | ||
- | The facets can be viewed as usual: | + | |
- | < | + | |
- | polytope > print_constraints($q); | + | |
- | Facets: | + | |
- | 0: 2x1 - x2 >= -1 | + | |
- | 1: 0 >= -1 | + | |
- | 2: -x1 + 2x2 >= -1 | + | |
- | 3: x1 + x2 >= 1 | + | |
- | </ | + | |
- | + | ||
- | {{ : | + | |
- | + | ||
- | Note that the upper right part (including the red vertices) arises from truncation of the polyhedron for visualization. | + | |
- | + | ||
- | + | ||
- | ==== Mixed-Integer Case ==== | + | |
- | + | ||
- | Let us now briefly discuss how to proceed if there are variables that are allowed to be integral. In this case there are several different types of information that one might be interested in. Let us first consider the question of how to compute the convex hull of all feasible integral variables, i.e., we consider the projection to the integral variables and then consider the convex hull of all feasible solutions. We only consider the bounded case, i.e., the original polyhedron is bounded. | + | |
- | + | ||
- | Consider the following example: | + | |
- | < | + | |
- | Minimize | + | |
- | | + | |
- | Subject to | + | |
- | C1: s1 - 10 x1 <= 0 | + | |
- | C2: s2 - 10 x2 <= 0 | + | |
- | C3: s1 + s2 <= 1.5 | + | |
- | C4: s1 + s2 >= 0.5 | + | |
- | Bounds | + | |
- | 0 <= s1 | + | |
- | 0 <= s2 | + | |
- | 0 <= x1 <= 1 | + | |
- | 0 <= x2 <= 1 | + | |
- | General | + | |
- | x1 x2 | + | |
- | End | + | |
- | </ | + | |
- | In this example there are two integral variables '' | + | |
- | < | + | |
- | polytope > $m=lp2poly(' | + | |
- | polytope > $p = new Polytope< | + | |
- | </ | + | |
- | We project the polyhedron in '' | + | |
- | < | + | |
- | polytope > $q=projection($p, | + | |
- | </ | + | |
- | We now construct the convex hull of all feasible points as above: | + | |
- | < | + | |
- | polytope > $s=new Polytope(POINTS=> | + | |
- | polytope > print_constraints($s); | + | |
- | Facets: | + | |
- | 0: | + | |
- | 1: | + | |
- | 2: | + | |
- | </ | + | |
- | Thus, as expected, the convex hull equals the triangle with vertices '' | + | |
- | + | ||
- | ===== 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: | + | |
- | + | ||
- | 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// | + | |
- | + | ||
- | Using the polymake extension [[https:// | + | |
- | + | ||
- | ==== Example: Explicit Matrix ==== | + | |
- | + | ||
- | In a first example, we directly create an integral matrix | + | |
- | < | + | |
- | polytope > $M=new Matrix< | + | |
- | </ | + | |
- | The total unimodularity of this matrix can be checked as follows: | + | |
- | < | + | |
- | polytope > print is_totally_unimodular($M); | + | |
- | 1 | + | |
- | </ | + | |
- | Thus, the given matrix is totally unimodular. | + | |
- | + | ||
- | ==== Example: Matrix from Input ==== | + | |
- | + | ||
- | In the second example, we reuse the file '' | + | |
- | < | + | |
- | polytope > $f=lp2poly(' | + | |
- | polytope > $p = new Polytope< | + | |
- | </ | + | |
- | We now want to check whether the constraint matrix defined by the inequalities is totally unimodular (note that there are no equations in this example). Thus we first extract the inequality matrix without the first column (as an integer matrix) and then perform the test: | + | |
- | < | + | |
- | polytope > $A = new Matrix< | + | |
- | polytope > print is_totally_unimodular($A); | + | |
- | 1 | + | |
- | </ | + | |
- | Thus, this matrix is totally unimodular as well. | + | |
- | + | ||
- | ===== Total dual integrality ===== | + | |
- | + | ||
- | Computations with respect to total dual integrality (TDI) can also be performed in polymake. Currently (August 2013), you need the perpetual beta version of polymake to access this functionality. | + | |
- | + | ||
- | The main functions are: | + | |
- | * The function '' | + | |
- | * The function '' | + | |
- | + | ||
- | Note that the input has to be full-dimensional in order to use these functions (this can be handled via '' | + | |
- | + | ||
- | To demonstrate the behavior of these functions, consider the 5-cycle example from above again: | + | |
- | < | + | |
- | Maximize | + | |
- | obj: x#1 + x#2 + x#3 + x#4 + x#5 | + | |
- | Subject to | + | |
- | | + | |
- | | + | |
- | | + | |
- | | + | |
- | | + | |
- | Bounds | + | |
- | 0 <= x#1 <= 1 | + | |
- | 0 <= x#2 <= 1 | + | |
- | 0 <= x#3 <= 1 | + | |
- | 0 <= x#4 <= 1 | + | |
- | 0 <= x#5 <= 1 | + | |
- | General | + | |
- | x#1 x#2 x#3 x#4 x#5 | + | |
- | End | + | |
- | </ | + | |
- | + | ||
- | Let us test whether the inequality system of this example is TDI. Thus, we first load the data as usual: | + | |
- | < | + | |
- | polytope > $f = lp2poly(' | + | |
- | polytope > $p = new Polytope< | + | |
- | </ | + | |
- | We now extract the corresponding inequality system and check it for TDIness: | + | |
- | < | + | |
- | polytope> | + | |
- | polytope > totally_dual_integral($M); | + | |
- | 0 | + | |
- | </ | + | |
- | 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, | + | |
- | < | + | |
- | polytope > $q = make_totally_dual_integral($p); | + | |
- | polytope > print_constraints($q); | + | |
- | Inequalities: | + | |
- | 0: x5 >= 0 | + | |
- | 1: x4 >= 0 | + | |
- | 2: x3 >= 0 | + | |
- | 3: x2 >= 0 | + | |
- | 4: x1 >= 0 | + | |
- | 5: -x1 - x2 >= -1 | + | |
- | 6: -x1 - x5 >= -1 | + | |
- | 7: -x2 - x3 >= -1 | + | |
- | 8: -x3 - x4 >= -1 | + | |
- | 9: -x4 - x5 >= -1 | + | |
- | 10: -x1 - x2 - x3 - x4 - x5 >= -5/2 | + | |
- | 11: 0 >= -1 | + | |
- | </ | + | |
- | As expected, the right hand side is non integral (otherwise, we know from general theory that the polytope would be integral as well). | + | |
- | + | ||
- | ===== 0/ | + | |
- | + | ||
- | ===== Chvátal-Gomory Closure ===== | + | |
- | + | ||
- | In the following we want to briefly show how closures of polytopes with respect to certain cutting-plane operators can be computed. We consider the two well-known cutting-plane operators here. The first one is the Chvátal-Gomory generator and the second one is the Lift-and-project operator as defined by Balas. For simplicity we will assume that the considered polytope is full-dimensional. | + | |
- | + | ||
- | ==== Chvátal-Gomory Closure - Example 1 ==== | + | |
- | + | ||
- | We first consider the polytope from the stable set problem from above: | + | |
- | < | + | |
- | Maximize | + | |
- | obj: x#1 + x#2 + x#3 + x#4 + x#5 | + | |
- | Subject to | + | |
- | | + | |
- | | + | |
- | | + | |
- | | + | |
- | | + | |
- | Bounds | + | |
- | 0 <= x#1 <= 1 | + | |
- | 0 <= x#2 <= 1 | + | |
- | 0 <= x#3 <= 1 | + | |
- | 0 <= x#4 <= 1 | + | |
- | 0 <= x#5 <= 1 | + | |
- | General | + | |
- | x#1 x#2 x#3 x#4 x#5 | + | |
- | End | + | |
- | </ | + | |
- | As before we read in the file using '' | + | |
- | < | + | |
- | polytope > $f = lp2poly(' | + | |
- | polytope > $p = new Polytope< | + | |
- | </ | + | |
- | + | ||
- | The Chvatal-Gomory closure of a polytope can be computed with the function '' | + | |
- | < | + | |
- | polytope > $g=gc_closure($p); | + | |
- | polytope > print print_constraints($g); | + | |
- | Inequalities: | + | |
- | 0: x5 >= 0 | + | |
- | 1: x4 >= 0 | + | |
- | 2: x3 >= 0 | + | |
- | 3: x2 >= 0 | + | |
- | 4: x1 >= 0 | + | |
- | 5: -x1 - x2 >= -1 | + | |
- | 6: -x1 - x5 >= -1 | + | |
- | 7: -x2 - x3 >= -1 | + | |
- | 8: -x3 - x4 >= -1 | + | |
- | 9: -x4 - x5 >= -1 | + | |
- | 10: -x1 - x2 - x3 - x4 - x5 >= -2 | + | |
- | 11: 0 >= -1 | + | |
- | </ | + | |
- | Let us check whether the resulting polytope is integral by examining its vertices: | + | |
- | < | + | |
- | polytope > print $g-> | + | |
- | 1 0 0 0 1 0 | + | |
- | 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 | + | |
- | </ | + | |
- | Thus, in this case, we have obtained the integer hull by one step of the Chvatal-Gomory-closure. | + | |
- | + | ||
- | ==== Example 2 for Chvátal-Gomory Closure ==== | + | |
- | + | ||
- | 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: | + | |
- | < | + | |
- | polytope > $M = new Matrix< | + | |
- | polytope > $t = new Polytope< | + | |
- | polytope > $t1 = gc_closure($t); | + | |
- | polytope > $t1-> | + | |
- | polytope > print_constraints($t1); | + | |
- | Facets: | + | |
- | 0: x4 >= 0 | + | |
- | 1: x3 >= 0 | + | |
- | 2: x2 >= 0 | + | |
- | 3: x1 >= 0 | + | |
- | 4: -x1 - x2 - x3 >= -1 | + | |
- | 5: -x1 - x2 - x4 >= -1 | + | |
- | 6: -x1 - x3 - x4 >= -1 | + | |
- | 7: -x2 - x3 - x4 >= -1 | + | |
- | polytope > $t1-> | + | |
- | polytope > print $t1-> | + | |
- | 1 1 0 0 0 | + | |
- | 1 0 0 0 0 | + | |
- | 1 1/3 1/3 1/3 1/3 | + | |
- | 1 0 1 0 0 | + | |
- | 1 0 0 1 0 | + | |
- | 1 0 0 0 1 | + | |
- | </ | + | |
- | Thus, one round was not enough to produce an integral polytope. However, one more round is enough: | + | |
- | < | + | |
- | polytope > $t2 = gc_closure($t1); | + | |
- | polytope > $t2-> | + | |
- | polytope > print_constraints($t2); | + | |
- | Facets: | + | |
- | 0: x4 >= 0 | + | |
- | 1: x3 >= 0 | + | |
- | 2: x2 >= 0 | + | |
- | 3: x1 >= 0 | + | |
- | 4: -x1 - x2 - x3 - x4 >= -1 | + | |
- | </ | + | |
- | + | ||
- | + | ||
- | ===== Lift-and-project closure ===== | + | |
- | + | ||
- | ===== Other Software ===== | + | |
- | + | ||
- | + | ||