3
$\begingroup$

Background

Discussion from http://blog.geomblog.org/2005/10/sampling-from-simplex.html and https://cs.stackexchange.com/questions/3227/uniform-sampling-from-a-simplex have shown algorithms of sampling from unit simplex uniformly:

  • Let $X_i \sim Exponential(1)$. Then $X_i/\sum_i X_i$ gives an uniform sample in unit simplex.

  • Let $X_i \sim Uniform([0,1])$. Then the difference of ordered statistic of ${X_i}$ would give an uniform sample in unit simplex.

Problem

Is there an algorithm to sample uniformly at random from the intersection of the unit hypercube and a simplex? i.e.

$\sum_i X_i = c$, where $c \in [0, n]$ and $X_i \in [0,1], \forall i$

I have seen someone asked the exactly same question at https://math.stackexchange.com/questions/4268521/random-sampling-from-the-intersection-of-cube-and-simplex, but it is not solved at the time of asking.

Attempt

Here is my naive approach:

  1. Find out all vertices of the set, denoted $V:=\{V_i\}$.
  2. Draw a uniform sample from unit $|V|$-simplex, denoted $W:=\{W_i\}$.
  3. Since the set is convex, $WV:=\{W_iV_i\}$ would be a sample from the set.

The problem of this naive algorithm is the complexity. Note that it would be $O(2^n)$ in both time and space.

Is there a better algorithm? State the question more precisely.

$\endgroup$
0

1 Answer 1

2
$\begingroup$

I will suggest two algorithms that work if $c$ is not too large, but probably will work poorly for large values of $c$ (which I realize is the important case). As such, I recognize this might not be useful for you in practice. I hope someone else will be able to offer a better answer.

Algorithm 1: rejection sampling

  • Repeat, until $x_1,\dots,x_n$ are all in the hypercube:

    • For $i:=1,2,\dots,n$:

      • Sample $u_i$ from $\text{Exp}(1)$ (the exponential distribution).
    • Set $u := u_1 + \dots + u_n$. Set $x_i := c \cdot u_i/u$ for each $i$.

In other words, in each iteration we sample uniformly at random from the simplex, check if it is the unit hypercube, and if not, repeat until you find a sample that is in the simplex.

Correctness: The correctness of Algorithm 1 follows immediately from the correctness of rejection sampling.

Running time: If $c \le n/\log(n)$, then Algorithm 1 will be efficient. In particular, the expected number of times you repeat will be about some constant. But if $c$ gets larger, the running time of Algorithm 1 will become much worse (probably exponentially rapidly).

Algorithm 2: approximate sampling

  • Set $\alpha := (n+4\sqrt{n})/c$.

  • Repeat, until $x_1,\dots,x_n$ are all in the hypercube:

    • For $i:=1,2,\dots,n$:

      • Repeatedly sample $u_i$ from $\text{Exp}(1)$, until $u_i \le \alpha$ (a truncated exponential distribution).
    • Set $u := u_1 + \dots + u_n$. Set $x_i := c \cdot u_i/u$ for each $i$.

I claim this is efficient if $c$ is not too large, and it samples from approximately the correct distribution if $c$ is not too large.

Running time: Intuitively, the number of iterations of the outer loop is significantly less than Algorithm 1, because Algorithm 2 will "retry" values of $u_i$ that are almost guaranteed to lead to rejection.

Correctness: This algorithm is approximately correct, as long as it doesn't make too many iterations. The rough idea is that if we reject a sample of $u_i$ (because it is larger than $\alpha$), then this sample would almost certainly be rejected as outside the hypercube, so we can immediately redraw it. We'll prove this by showing that the probability that they deviate from the proper distribution is not too large. Let's classify an iteration of the outer loop of either algorithm as good if $u \le n+4\sqrt{n}$, or bad otherwise.

Notice that the probability of an iteration being bad is small (about $3 \times 10^{-5}$ by the Central Limit Theorem, as the exponential distribution has mean 1 and standard deviation 1). As long as Algorithm 2 doesn't do too many iterations, it's likely that all iterations of Algorithm 2 will be good.

Moreover, in a good iteration, the probability distribution of $x_1,\dots,x_n$, conditioned on them being in the hypercube, is the same for both algorithms. Why? In Algorithm 1, if $u_i > \alpha$ in a good iteration, then $$x_i= c \cdot u_i/u > c \alpha/(n+4\sqrt{n}) =1,$$ so $x_1,\dots,x_n$ will surely not be in the hypercube and will be rejected anyway. Hence truncating the exponential distribution in Algorithm 2 does not change the distribution of accepted outputs (conditioned on it being a good iteration).

I don't know how large $c$ can be with Algorithm 2. I think for large enough values of $c$, Algorithm 2 will still exhibit exponential running time, but hopefully Algorithm 2 can support values of $c$ that are a bit larger than Algorithm 1.

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