# Uniform Sampling Subject to Linear Equalities and Non-Negativity Constraint

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.

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)

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:

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

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.

• 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! Jul 1 '15 at 23:24
• 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. Jul 2 '15 at 1:13
• 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! Jul 6 '15 at 17:07

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

Or does the dimension problem kill that too ?

• 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... Jul 6 '15 at 17:14
• 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 Jul 6 '15 at 22:11
• 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. Jul 6 '15 at 23:12
• 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 Jul 7 '15 at 7:18
• I'm happy to try different methods, even if it's an approximation. How would you use VBEM or EP here? Jul 7 '15 at 11:50

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

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