25
$\begingroup$

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$.

$\endgroup$

7 Answers 7

18
$\begingroup$

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.

$\endgroup$
18
$\begingroup$

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$
$\endgroup$
3
  • $\begingroup$ Can you suggest a reference for the analysis if this algorithm, please? $\endgroup$
    – Elle Najt
    Apr 5, 2018 at 15:16
  • 1
    $\begingroup$ Is this the kind of thing you were looking for? researchgate.net/publication/… $\endgroup$ Apr 6, 2018 at 12:12
  • $\begingroup$ Is this implemented somewhere? $\endgroup$ Dec 6, 2023 at 16:29
4
$\begingroup$

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.

$\endgroup$
5
  • 1
    $\begingroup$ 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. $\endgroup$ Dec 27, 2009 at 16:33
  • $\begingroup$ That idea won't work if you reject 999,999 out of a million samples. There must be some intrinsic way of doing this. $\endgroup$ Dec 27, 2009 at 16:43
  • 1
    $\begingroup$ 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. $\endgroup$ Dec 27, 2009 at 20:53
  • $\begingroup$ 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. $\endgroup$ Dec 29, 2009 at 10:37
  • $\begingroup$ @SteveHuntsman: The idea of covering the polytope with small rectangles is simple, nice and really powerful. Great! $\endgroup$
    – dohmatob
    Nov 2, 2015 at 13:16
2
$\begingroup$

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

$\endgroup$
1
$\begingroup$

Rejection sampling.

$\endgroup$
1
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ Additionally, having to convert to a list of vertices is going to be slow even in principle, since the number of vertices can explode exponentially in the general case. $\endgroup$ Apr 7, 2018 at 16:51
0
$\begingroup$

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]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.