9
$\begingroup$

I'm trying to sample uniformly on the intersections of faces of several simplicies, with all coordinates being non-negative. That is, given constraints $$A\vec{w}=\vec{b} \ \ and \ \ \vec{w} \geq \vec{0} \ \ and \ \ \sum w_i = 1,$$ I want to sample $\vec{w}$ uniformly. $A$'s dimension is about $100 \times 10000$. A concrete example will be: $$A = \begin{bmatrix} 1 & 1 & 1 \\ 0 & 1 & 2 \end{bmatrix}, \ b=\begin{bmatrix} 1 \\ 0.7 \end{bmatrix}$$, sample $\vec{w}$ uniformly from $Aw=b$ subject to $\vec{w} \geq \vec{0}$ and $\sum w_i = 1$ (This makes the sampling space bounded). Below is a graphical representation of the problem -- to sample uniformly from the red intersection line.

a busy cat

I am well aware that rejection-sampling and MCMC sampling can theoretically solve this problem. However, I have already implemented both approaches in programming, and neither of these two methods performs well enough. This is because the dimension of my sampling space usually goes up to 10000, and rejection sampling simply throws away too many points and MCMC is taking forever to converge. Therefore, I'm desperate to try new methods. Many thanks in advance!! (please do not provide answers using rejection sampling; methods that already have open-source programming implementations are favored)

Improve this question
$\endgroup$

4 Answers 4

Reset to default
2
$\begingroup$

There is some work on billiard walks-- the idea is to start travel within the polytope until a boundary is reached, and then reflect the direction at the boundary and continue. Details are here:

Elena Gryazinaa, Boris Polyak. Random sampling: Billiard Walk algorithm, 2014.

The authors argue that it mixes much more rapidly than Hit-and-Run sampling. I imagine that implementing a billiard walk should be fairly straightforward (since the algorithm's so simple).

Improve this answer
$\endgroup$
0
$\begingroup$

The hit-and-run algorithm and variants are popular choices. These are Monte Carlo methods but should be much better than rejection sampling. Unfortunately I don't know of a canonical reference.

Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ Thanks for answering! I've implemented hit-and-run, but it is not doing well in dimension high as 10000. But can you elaborate on the variant method? Thanks! $\endgroup$
    – Miller Zhu
    Commented Jul 1, 2015 at 23:24
  • $\begingroup$ I don't actually know about them myself, I just know from googling that there are variants with different runtimes and guarantees. But I see what you mean; it seems like the algorithms where anything can be proven are at least on the order of $n^3$ and often higher, which is prohibitive for $n=10000$. I'm not sure what to do about that. $\endgroup$
    – Noah Stein
    Commented Jul 2, 2015 at 1:13
  • $\begingroup$ Yes indeed, the high dimension seems a killer for such MCMC sampling methods, which generally require complexity of $O(n^3)$. But thanks for answering! $\endgroup$
    – Miller Zhu
    Commented Jul 6, 2015 at 17:07
0
$\begingroup$

What about Gibbs sampling ?

At a given point select two dimensions at random, compute the conditional and sample from that ?

Or does the dimension problem kill that too ?

Improve this answer
$\endgroup$
6
  • 1
    $\begingroup$ I'm afraid the dimension problem will kill this too... Since MCMC method requires at least $O(n^3)$ runtime complexity, the dimension will be a big problem. But indeed, the problem itself is very hard, in that the current known methods are all MCMC and requires $O(n^3)$. But now we want to deal with $n \geq 10000$, which simply sounds impossible. Until a totally new method is found, no MCMC scheme might be able to solve this, I guess... $\endgroup$
    – Miller Zhu
    Commented Jul 6, 2015 at 17:14
  • $\begingroup$ Actually, do you have a proof that your method defines a proper probability distribution (int equal to one) ? Now that i loin at it carefully, i dont see why the simplex you define should have finite volume $\endgroup$
    – Guillaume Dehaene
    Commented Jul 6, 2015 at 22:11
  • $\begingroup$ Sorry I forgot to point out that there's one definite linear constraint: $\sum w_i=1$. This should make the intersection of faces of simplices bounded. $\endgroup$
    – Miller Zhu
    Commented Jul 6, 2015 at 23:12
  • $\begingroup$ I see. Everything works with that. Would you be satisfied with an approximation such as vbem or ep ? I m curious whether they would be appropriate here $\endgroup$
    – Guillaume Dehaene
    Commented Jul 7, 2015 at 7:18
  • $\begingroup$ I'm happy to try different methods, even if it's an approximation. How would you use VBEM or EP here? $\endgroup$
    – Miller Zhu
    Commented Jul 7, 2015 at 11:50
0
$\begingroup$

If you'd be satisfied with an approximation, you might try the following method. In large dimensions, an isotropic convex polytope resembles a ball:

enter image description here

(See slides 6 through 9 here for details.) As such, if you uniformly sample from an appropriately stretched ball (to account for any lack of isotropy in your polytope), then your samples will land in the polytope a lot more frequently than standard rejection sampling. In doing so, you will fail to sample the "tentacles" of the polytope, but these make up a vanishing fraction of the polytope by concentration of volume.

Improve this answer
$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .