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)
 A: 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).
A: 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.
A: 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 ?
A: 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.
