# Efficient Sampling of Points from the Unit Simplex with Convex Hull Containing a Given Set

I have a set of points $$P = \{p_1, \ldots, p_k\}$$, where each $$p_i \in \mathbb{R}^n$$ and these points reside on the unit simplex, such that for each $$i$$, we have $$p_i > 0$$ and $$\sum_{j=1}^{n} p_{ij} = 1$$, where $$p_{ij}$$ denotes the $$j$$-th component of point $$p_i$$.

I want to sample $$n$$ points $$\{s_1, \ldots, s_n\}$$ from the unit simplex such that the all points in $$P$$ is contains in the the convex hull formed by these sampled, $$P \subset \text{conv}(s_1, \ldots, s_n)$$.

A straightforward approach is to uniformly sample points from the unit simplex and discard any set of points whose convex hull does not contain $$P$$. However, this method results in a high rejection rate, which seems inefficient.

Could anyone suggest a more efficient algorithm or method to achieve this sampling? Any insights on how to directly sample a set of points that satisfy the convex hull condition would be greatly appreciated.

Thank!

• How big can the dimension $n$ be in the problem that you need such a simulation?
– Amir
Commented Feb 11 at 15:09
• Hi Amir! The dimension n in our problem is quite big, roughly in the lower hundreds. The size of the set P is quite small, between 2 and 10. I'll update my question with this information. Commented Feb 13 at 1:19

The one quick improvement that you should try first is to immediately (i.e. upon sampling, not upon completing the set) discard points that lie inside the convex hull of $$P$$, because those can never be part of the result. I wouldn’t be surprised if this alone would lead to an acceptable rejection rate (but that of course depends on your application).

If that’s not good enough, you could try the following approach. (Everything in the following lives on the $$(n-1)$$-dimensional hyperplane that contains the unit simplex, so e.g. when I talk about a hyperplane, it’s an $$(n-2)$$-dimensional hyperplane within that $$(n-1)$$-dimensional hyperplane.)

For any hyperplane that touches the convex hull of $$P$$, i.e. that contains at least one $$p_j$$ and has all the $$p_i$$ on one side, at least one of the $$s_\ell$$ must be on the other side. Conversely, if the $$s_\ell$$ satisfy this constraint for all such hyperplanes, then their convex hull contains $$P$$.

Select some set of such hyperplanes, including at least the faces of the convex hull of $$P$$. Whether it makes sense to include further hyperplanes might depend on the sorts of sets $$P$$ you have. Each hyperplane in the set represents a constraint that at least one of the $$s_\ell$$ must lie on the other side of the hyperplane than $$P$$.

Now sample candidates for the $$s_\ell$$ and keep track of which of the constraints they satisfy, i.e. form a bitset where the bit corresponding to a constraint is set iff the candidate satisfies the constraint. Again, discard any candidates that don’t satisfy any of the constraints (since in particular they lie on the inside of all faces of the convex hull of $$P$$, and thus inside that convex hull, and thus can’t be part of the result). Keep going until at least one subset of $$n$$ of the candidates satisfies all constraints; for all such subsets, check whether their convex hull actually contains $$P$$ (it might not since you didn’t include all possible constraints) and uniformly randomly select one of the subsets that pass the check (or if none does, keep going).

There are some implementational details to be filled in for efficiently determining which $$n$$-point subsets satisfy all the constraints. The most straightforward approach would be to just form all $$n$$-point subsets and OR their bitsets to see whether a full bitset results, but you could improve on that, e.g. keep track of which constraints are satisfied by only one point so you know those points need to be included; or keep a count of the constraints satisfied by each point and avoid generating subsets whose count total is less than the number of constraints.

You don’t specify the distribution you’re aiming for, so I guess either you don’t particularly care about the details of the distribution or you’re aiming for the distribution that your rejection scheme would yield. If the latter, I don’t think this approach would replicate that distribution; I suspect it favours “tighter” solutions because those are less likely to have formed solutions with other points at earlier stages. If you want to avoid that bias and exactly replicate the distribution from your rejection scheme, you could determine how many candidates you typically need in order to get an admissible subset of $$n$$, and then always generate that many candidates. But that might sometimes generate too many $$n$$-point subsets with full bitsets whose convex hulls you need to construct.