Sampling points in a polytope

[satya123999]

Dear all,

I am defining a polytope in polymake through a set of linear inequalities and could compute the vertices also. In addition, I would like to sample the points in the polytope (similar to a hit-and-run sampling method). Is there a possibility to accomplish this in the polymake ?

Many Thanks & Regards,
Satya

[assarf]

Hi there,

I do not know the "hit-and-run sampling method". At least the name of the method does not ring a bell. But maybe this is what you search?

Code: Select all

rand_inner_points(P, n; Options) -> Polytope Produce a polytope with n random points from the input polytope P. Each generated point is a convex linear combination of the input vertices with uniformly distributed random coefficients. Thus, the output points can't in general be expected to be distributed uniformly within the input polytope; cf. unirand for this. The polytope must be BOUNDED. Arguments: Polytope P the input polytope Int n the number of random points Options: seed => Int controls the outcome of the random number generator; fixing a seed number guarantees the same outcome. Returns Polytope

or this

Code: Select all

unirand(Polytope, n; Options) -> Polytope Produce a polytope with n random points that are uniformly distributed within the given polytope P. P must be bounded and full-dimensional. Arguments: P Polytope Int n the number of random points Options: boundary => Bool forces the points to lie on the boundary of the given polytope seed => Int controls the outcome of the random number generator; fixing a seed number guarantees the same outcome. Returns Polytope

Both functions return a polytope, where the property POINTS contains all the created points.

[satya123999]

Thank you very much ! Is there such a function if the polytope is unbounded ? (i.e. to sample 'n' points)

[assarf]

no, there is nothing for unbounded polyhedra.

one thing you can do is to transform your polyhedron via a projective transformation into a bounded polytope, use one of the two functions above and reverse the transformation afterwards. Or you could produce random convex combinations on your own.

either way, if you are aiming for a uniform distribution: this won't work with unbounded polyhedra.

[satya123999]

Can we try with a small example ?
Suppose I have a set of inequalities in lp format e.g.

Code: Select all

Maximize Subject To c0: - a1 - a2 + a3 = 0 cc2l0: - a1 + a3 <= -2 c1: a2 = -2 cc2l1: a1 + a2 - a3 <= 1 cc3l1: a1 <= -1 c2: - a2 + a3 = -1 cc1l2: - a1 - a2 + a3 <= 0 Bounds a1 free a2 free a3 free End

The set of vertices are

Code: Select all

-1.665334537e-16 -1 3.330669074e-16 -1 1 -1 -2 -3

I guess the bound (P) → Polytope function doesn't work on polytope with negative co-ordinates. So how do you suggest to proceed ?

[assarf]

Hi,

first of all: when you use the lp2poly script you will get a polytope with floating point coordinates. As floats are "evil" you might want to calculate with Rationals instead. So use:

Code: Select all

$p = convert_to<Rational>(lp2poly('test.lp'));

as described here or you input the inequalities directly into polymake, without doing the lp-format conversion.

As for your example: The floating point numbers "-1.665334537e-16" and "3.330669074e-16" are actually zero (rounding issues don't happen with rational coordinates).

But other than that I cannot reproduce your output. When I take your lp file and ask for the vertices (with rational or floating coordinates; with the release and the perpetual beta version of polymake) I get only \( (1, -1, -2, -3) \) as a vertex. And when you look at it more closely you see that the ray \(\lambda(-1, 0, -1)\) violates the equation c2: - a2 + a3 = -1. So the output you posted here, is somewhat strange.
From where I am standing: your case is actually a bounded polytope. It is the convex hull of the point \((-1,-2,-3)\in \mathbb{R}^3\), so everything is fine here

But back to your question:
When you have an unbounded but pointed polyhedron, you want to transform it in a way such that it lies in the positive orthant (meaning every coordinate is positive). That is basic linear algebra (rotating, translating, basis transformation). Then you can use the bound function to do the projective transformation in order to create the bounded version of that polyhedron. Or you do it yourself by transforming the hyperplane \((1,0,\ldots,0)\) to \((1,-1,\ldots,-1)\), for example.
Then you sample everything and reverse every transformation you did before to get the points in your polytope.

[satya123999]

Thank you I was reading the lp file in float but now after changing it to Rational, I also get a single vertex