Uniformly Sampling from Convex Polytopes

Question

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 \times n$ matrix, $x \in \mathbb{R}^n$, and $b \in \mathbb{R}^m$.) 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$.

Answer 1

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.

Answer 2

Hit and run sampling will perform much better than rejection sampling for higher dimensions.

It converges to a uniform distribution in polynomial time. The basic idea is to start at any point $x_0$ inside the polytope, then follow the following procedure iteratively:

• pick a direction $\alpha$ uniformly at random.
• find the minimum and maximum value of $\theta$ such that $x_n+\theta\cdot\alpha$ is contained in the polytope
• pick a $\theta^*$ uniformly at random from $[\theta_{min},\theta_{max}]$
• your new point is $x_{n+1}=x_n+\theta^*\cdot\alpha$

Answer 3

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.

Answer 4

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

Answer 5

Rejection sampling.

Answer 6

Rejection sampling definitely works if you are able to find a superset $Q$ of the polytope $P$ from which you can sample. If you sample a point from that superset, the probability that it gets accepted is equal to the ratio $\frac{\text{Vol}(P)}{\text{Vol}(Q)}$, so $Q$ should be as small as possible. For instance, it is sample to sample from $Q$ if it is a box or a ball.

In the case where the polytope is specified as a list of inequalities, finding the smallest enclosing ball can be quite hard.

Contrary to what Simon Barthelmé mentions, Boyd and Vandenberghe do not deal with this problem. Actually they deal with the case where the vertices of the polytope are available. Going from the list of inequalities to the set of vertices is also hard (I am actually looking for a MATLAB implementation of that).

One possible approach is to find a small box enclosing the polytope. The box is defined by a set of coordinates $(b_i^{\text{min}},b_i^{\text{max}}), i=1\dots n$, and each coordinate can be found by : $$ b_i^{\text{min}} = \arg \min_x x_i \quad \text{subject to } A x \leq b $$ $$ b_i^{\text{max}}= \arg \max_x x_i \quad \text{subject to } A x \leq b $$ Those are linear programs for which you can use your favourite solver.

Answer 7

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

1. sample in $\sum{x_i} = 1$ by taking i.i.d. exponentials scaled to sum 1
2. 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]