Welcome to MathOverflow

Question

2

How to choose a point uniformly from a convex polytope $P \subset [0,1]^n$ defined by some inequalities, Ax < b ? (Here A is an m-by-n matrix, x is n-by-1 and b is m-by-1.) I imagine that you could start with a uniformly chosen point in the cube and do some process to get a point with Ax < b.

flag	asked Dec 27 2009 at 4:57 John Mangual 6,150●4●40

David Eppstein · Answer 1 · 2009-12-27 05:07:47Z UTC

See the answers to this question — uniform sampling from polytopes forms the basis for the known algorithms for calculating their volumes. The methods from the papers mentioned in those answers mostly take the form of a random walk inside the polytope; they differ in the details of the walk and in the analysis of its mixing time.

Simon Barthelmé · Answer 2 · 2009-12-27 10:46:01Z UTC

2

Rejection sampling will definitely work. Take a hypercube that you know contains the polytope, sample from the hypercube, and accept only those samples that belong in the polytope. However if the relative volume of the polytope is small you'll end up rejecting most samples, and the method might get painfully slow. Depending on your needs you might want to find, e.g., the smallest enclosing ball first, so that you can draw your uniform samples from that instead of the hypercube.
See Boyd & Vanderberghe's book on Convex Optimisation (it's online) for finding smallest enclosing sets.

answered Dec 27 2009 at 10:46

Simon Barthelmé
31●1

	For practicality it might be better to cover the polytope with boxes. Doing this optimally is hard/ill-posed/etc, but I would imagine it wouldn't be too challenging to fix the number of elements in a cover (start with a uniform "layering" of the hypercube) and minimize volume of the cover under the constraints of fixed cardinality and no overlap. Then rejection sample from that. – Steve Huntsman Dec 27 2009 at 16:33
	That idea won't work if you reject 999,999 out of a million samples. There must be some intrinsic way of doing this. – John Mangual Dec 27 2009 at 16:43
	Also, you might set up the box covering by taking a uniform mesh on the hypercube, then discarding boxes that don't intersect the polytope at all. Keep indices for the remaining boxes as well as a supplementary Boolean variable to indicate whether or not the box is entirely contained in the polytope. Sample uniformly on indices, then rejection sample or uniformly sample on the resulting boxes depending on the Boolean variable you've already set. This middle road is probably close to the best you can do in silico. – Steve Huntsman Dec 27 2009 at 20:53
	John: if you want an intrinsic mechanism, you can use MCMC techniques to sample from the density asymptotically. For example, it would be be fairly easy to implement Gibbs sampling, but it might take a long time for the chain to converge. In general the problem you want to solve is very hard, computationally speaking. To find an efficient algorithm, you are going to make some assumptions about the structure of your polytope, or to use some preprocessing to find out more about it. That's what Steve is suggesting. – Simon Barthelmé Dec 29 2009 at 10:37

Steve Huntsman · Answer 3 · 2009-12-27 05:03:47Z UTC

1

Rejection sampling.

answered Dec 27 2009 at 5:03

Steve Huntsman
8,139●1●22●56

Tim Benham · Answer 4 · 2012-01-07 15:25:34Z UTC

1

If you use MATLAB cprnd on their File Exchange solves the problem.

answered Jan 7 2012 at 15:25

Tim Benham
11●1

denis · Answer 5 · 2012-06-10 13:00:04Z UTC

A case with a fast simple method: to sample the "right simplex" $\ \sum{x_i} \le 1,\ x_i \ge 0$:

sample in $\sum{x_i} = 1$ by taking i.i.d. exponentials scaled to sum 1
scale by random-uniform$^\frac{1}{dim}$.

(I have no idea how to generalize this.)

In Python with NumPy, this is

def random_simplex_sum1( N, dim ):
    """ N uniform-random points >= 0, sum x_i == 1 """
    X = np.random.exponential( size=(N,dim) )
    X /= X.sum(axis=1)[:,np.newaxis]
    return X

def random_simplex_le1( N, dim ):
    """ N uniform-random points >= 0, sum x_i <= 1 """ 
    return random_simplex_sum1( N, dim ) \
        * (np.random.uniform( size=N ) ** (1/dim)) [:,np.newaxis]

Uniformly Sampling from Convex Polytopes

Remember to vote up questions/answers you find interesting or helpful (requires 15 reputation points)

5 Answers

You can accept an answer to one of your own questions by clicking the check mark next to it. This awards 15 reputation points to the person who answered and 2 reputation points to you.