9

There is a surprising number of posts asking about the uniform sampling of convex polytopes (e.g. 1, 2, 3), and an equally surprising number of non-answers (e.g. a, b).

From what I have read, this seems to be a difficult problem; and known solutions consist in partitioning the polytope into a collection of simplices, select one randomly, and sample a uniformly random point from it. That is assuming that we know the vertices of the polytope in the first place (otherwise the problem is apparently NP-hard?).

From what I know, sampling uniformly random points from a d-dimensional simplex with d+1 known vertices is done in two steps:

  1. Sample d+1 weights from a Dirichlet distribution with parameters $(\alpha_k = 1)_{k=0..d}$
  2. Compute the corresponding weighted sum of vertices

Assuming this is correct, my question is: why can we not apply the same method to sample any convex polytope uniformly randomly? I.e. by computing a weighted sum of its vertices, with weights sampled from a Dirichlet distribution.

Is it because it has not been proven that this would result in a uniform distribution, or because there are known cases where it does not?

Jonathan H
  • 2,393

2 Answers2

9

Well, it seems the short answer is: "because it doesn't work".

It is still unclear to me why that is; and I should probably first try to understand why it works in the case of a simplex. But trying the same method on regular polygons yields a sampling density that is not uniform:

Sampling density using Dirichlet weights

Jonathan H
  • 2,393
  • 1
    Intuitively, I suppose this has to do with intersections. Take e.g. the pentagon, and label the vertices 1 to 5 in adjacent order. Draw a line from 1 to 3 and one from 2 to 4. The intersection is a point that you can express with at least two different Dirichlet samples. However, a point on the line from 1 to 2 is only expressible by one sample. This is not a problem in the simplex. – cangrejo Sep 17 '20 at 10:09
  • @cangrejo That's a nice way of looking at it! And it almost seems like there should be a way to "correct" for such duplicate combinations a posteriori, without having to partition the polytope with simplices... But I am sure this is more difficult than it seems as well :) – Jonathan H Sep 18 '20 at 00:05
  • This picture is really worth a thousand words. – William Bell Dec 02 '21 at 04:45
1

One might think that one can fix the above phenomenon by first sampling three random vertices (forming a triangle), then computing a weighted sum of its three vertices from a Dirichlet distribution. Unfortunately, this also fails. It does work much better than the above approach on an almost-regular pentagon.

almost-regular pentagon

But, on a not-at-all-regular pentagon, this approach also fails. Essentially, we have greater density if the triangle has a small area. not-regular pentagon

Perhaps one might try to fix this by giving lesser weight to triplets that form small-area triangles. But, if you have P points, you'll have to compute the area of P-choose-3 triangles, which is expensive.

Code below:

import numpy as np
import matplotlib.pyplot as plt
from sklearn.utils import check_random_state


rng = check_random_state(42)
N = 50000;


plt.figure();
pentagon = np.array([[0, 0], [0.8, 0], [1.2, 0.8], [0.4, 1.4], [-0.3, 0.8]])
plt.scatter(pentagon[:, 0], pentagon[:, 1]);
samples = np.zeros((N, 2));
for n in range(N):
triplet = pentagon[rng.choice(5, size=(3,), replace=False), :]
samples[n, :] = rng.dirichlet(np.ones(3)) @ triplet
plt.scatter(samples[:, 0], samples[:, 1], s=0.05);
plt.axis('equal');


plt.figure();
pentagon = np.array([[0, 0], [2, 0], [1.2, 0.8], [0.4, 1.4], [-0.3, 0.8]])
plt.scatter(pentagon[:, 0], pentagon[:, 1]);
samples = np.zeros((N, 2));
for n in range(N):
triplet = pentagon[rng.choice(5, size=(3,), replace=False), :]
samples[n, :] = rng.dirichlet(np.ones(3)) @ triplet
plt.scatter(samples[:, 0], samples[:, 1], s=0.05);
plt.axis('equal');
calmcc
  • 271