generating vertices of LP optimal solutions

[sdementen]

I have successfully solved an LP with pycddlib of the form

max c.x
st A.x <= b

I need to generate the set of vertices for the optimal solutions of both this primal problem and its dual. If I am not mistaken, it is exactly what pycddlib is good at ! But I can't figure out from the documentation how to do it. Could you give me some hints ?

If it can be of any help, here is the problem

begin
 19 13 rational
 0 1 0 0 0 0 0 0 0 0 0 0 0
 0 0 1 0 0 0 0 0 0 0 0 0 0
 0 0 0 1 0 0 0 0 0 0 0 0 0
 0 0 0 0 1 0 0 0 0 0 0 0 0
 0 0 0 0 0 1 0 0 0 0 0 0 0
 0 0 0 0 0 0 1 0 0 0 0 0 0
 0 0 0 0 0 0 0 1 0 0 0 0 0
 0 0 0 0 0 0 0 0 1 0 0 0 0
 0 0 0 0 0 0 0 0 0 1 0 0 0
 0 0 0 0 0 0 0 0 0 0 1 0 0
 0 0 0 0 0 0 0 0 0 0 0 1 0
 0 0 0 0 0 0 0 0 0 0 0 0 1
 2 -1 -1 -1 0 0 0 0 0 0 0 0 0
 1 0 0 0 -1 -1 -1 0 0 0 0 0 0
 2 0 0 0 0 0 0 -1 -1 -1 0 0 0
 2 0 0 0 0 0 0 0 0 0 -1 -1 -1
 1 -1 0 0 -1 0 0 -1 0 0 -1 0 0
 2 0 -1 0 0 -1 0 0 -1 0 0 -1 0
 3 0 0 -1 0 0 -1 0 0 -1 0 0 -1
end
maximize
 0 8 6 3 7 5 2 6 4 1 6 4 1

with the python code

import cdd
mat = cdd.Matrix([[ 0,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
       [ 0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
       [ 0,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  0],
       [ 0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1,  0,  0],
       [ 0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1,  0],
       [ 0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1],
       [ 2, -1, -1, -1,  0,  0,  0,  0,  0,  0,  0,  0,  0],
       [ 1,  0,  0,  0, -1, -1, -1,  0,  0,  0,  0,  0,  0],
       [ 2,  0,  0,  0,  0,  0,  0, -1, -1, -1,  0,  0,  0],
       [ 2,  0,  0,  0,  0,  0,  0,  0,  0,  0, -1, -1, -1],
       [ 1, -1,  0,  0, -1,  0,  0, -1,  0,  0, -1,  0,  0],
       [ 2,  0, -1,  0,  0, -1,  0,  0, -1,  0,  0, -1,  0],
       [ 3,  0,  0, -1,  0,  0, -1,  0,  0, -1,  0,  0, -1]], number_type='fraction')
mat.obj_type = cdd.LPObjType.MAX
mat.obj_func = [0, 8, 6, 3, 7, 5, 2, 6, 4, 1, 6, 4, 1]
print(mat)
print(mat.obj_func)
lp = cdd.LinProg(mat)
lp.solve()
lp.status == cdd.LPStatusType.OPTIMAL

print(lp.obj_value)
print(" ".join("{0}".format(val) for val in lp.primal_solution))
print(" ".join("{0}".format(val) for val in lp.dual_solution))

[stephane-caron]

I'm not sure what you mean by "the set of vertices for the optimal solutions" of an LP, as the optimal solution of the LP is a single point.

Maybe what you want, and what is commonly computed using cdd, is the set of vertices of your linear constraint A.x <= b. If that's so, you can use e.g. the function compute_polytope_vertices() from pypoman. (This function is just a thin wrapper around cdd, you can take a look at the code.)

Here is sample code for your example:

import numpy
import pypoman

A = numpy.array([
    [-1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
    [0, -1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
    [0,  0, -1,  0,  0,  0,  0,  0,  0,  0,  0,  0],
    [0,  0,  0, -1,  0,  0,  0,  0,  0,  0,  0,  0],
    [0,  0,  0,  0, -1,  0,  0,  0,  0,  0,  0,  0],
    [0,  0,  0,  0,  0, -1,  0,  0,  0,  0,  0,  0],
    [0,  0,  0,  0,  0,  0, -1,  0,  0,  0,  0,  0],
    [0,  0,  0,  0,  0,  0,  0, -1,  0,  0,  0,  0],
    [0,  0,  0,  0,  0,  0,  0,  0, -1,  0,  0,  0],
    [0,  0,  0,  0,  0,  0,  0,  0,  0, -1,  0,  0],
    [0,  0,  0,  0,  0,  0,  0,  0,  0,  0, -1,  0],
    [0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0, -1],
    [1,  1,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0],
    [0,  0,  0,  1,  1,  1,  0,  0,  0,  0,  0,  0],
    [0,  0,  0,  0,  0,  0,  1,  1,  1,  0,  0,  0],
    [0,  0,  0,  0,  0,  0,  0,  0,  0,  1,  1,  1],
    [1,  0,  0,  1,  0,  0,  1,  0,  0,  1,  0,  0],
    [0,  1,  0,  0,  1,  0,  0,  1,  0,  0,  1,  0],
    [0,  0,  1,  0,  0,  1,  0,  0,  1,  0,  0,  1]])
b = numpy.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 1, 2, 2, 1, 2, 3])

vertices = pypoman.compute_polytope_vertices(A, b)
print "Found %d vertices" % len(vertices)

Which finds 388 vertices on my machine.

[sdementen]

An LP can have as solution x a polyhedron no?

[stephane-caron]

Aah, that's right! The set of optimal solutions is indeed a face of the initial polyhedron, and what you want is to enumerate vertices of the face corresponding to your optimal solution.

I don't know how to do this directly with cdd.

What I would do off the top of my head is:

• Take the optimal solution x*
• Enumerate its active set, i.e. the set act of constraint indexes s.t. A_act x* == b_act
• Call vertex enumeration (function above) on the face {A x <= b, A_act x == b_act}; or, use the function call above and just discard vertices v for which A_act v != b_act
• Apply the same to the dual

That's definitely not efficient since the LP solver already knows the active set internally.

Looking forward to a more elegant solution :wink:

[mcmtroffaes]

One way is to add a constraint to fix the value of the optimal solution, and then to enumerate the vertices of the resulting extended problem. More elegantly, the simplex method can (in principle) also be used to enumerate all optimal solutions, without having to resolve an extended problem from scratch. I'm not sure how to do that with cdd. Regardless, none of this has any bearing on the Python interface of cddlib, so I will close this. I suggest that you address your maths questions to mathoverflow instead, to find out the best method for doing what you want to do.