# Random Sampling a linearly constrained region in n-dimensions…

Hi,

So here is my problem:

Given a nonlinear, discontinous, cost function $f(x_1,x_2,..,x_N)$ along with linear constraints $x_n \ge 0, \forall n$ $x_n \le c_n$ and $\sum_{n=1}^N x_n = 1$ find an optimal (local) solution by randomly sampling the feasible region. $c_n$ are just constants.

The issue I am having is indexing the space, since it is not a simple n-dimensional cube but rather a polytope(i believe its convex). Discretizing and enumerating all possible combinations of points to sample is much too hard since N =20. Approximating the polytope with a n-dimensional cube and sampling from the cube, only about 1% of the samples fall within the feasible region...which is inefficient if I'm trying to generate many samples.

I've tried finding the volume of the space analytically, however the complexity of computing the integral gets overwhelming for many dimensions.

I was wondering if anyone has come across this type of problem and has any recommendation as to different methods I could try to sample this space. Essentially, I need a good way to estimate the volume...am I looking at this the correct way?

any help would be greatly appreciated...

• I don't understand why sampling is the same problem as estimating the volume. Perhaps you can elaborate... I – Igor Rivin Sep 24 '11 at 14:55
• If $1$% of the samples from $[0,c_1]\times ... \times [0,c_n]$ satisfy $\sum x_i \le 1$, that should be usable, even if it seems abstractly inefficient. Also, if the volume helps, you can write out the volume using inclusion-exclusion. You get a sum with about $2^{20}$ terms which should be easy enough to evaluate. – Douglas Zare Sep 24 '11 at 20:16
• Igor, the connection I was making between volume estimation and sampling is that if you can estimate the volume of the polytope, and if you then pick a percentage of the volume randomly, that that percentage can be indexed as a point in the space (a sample) and used to evaluate the cost function. – user1 Sep 25 '11 at 6:41

## 3 Answers

Your constraints $x_n \geq 0$, $\sum_{n=1}^N x_n = 1$, are those for the standard simplex. You could try uniform sampling from the standard simplex, and then reject any sample that doesn't also satisfy the $x_n \leq c_n$ constraints.

An alternative to the procedure described in the linked paper above for uniform sampling from the standard simplex is to generate $n$ exponential(1) random variables $X_1, X_2, \ldots, X_n$ and let $Y_i = X_i/\sum_{i=1}^n X_i$. Then $(Y_1,Y_2,\ldots,Y_n)$ is uniformly distributed on the standard simplex. This can be thought of as generating a random vector from the symmetric Dirichlet distribution. (Also, generating exponential(1) random variables is easy; if $Z \sim U(0,1)$ then $-\ln(Z)$ has an exponential(1) distribution.) Once again, you would then reject any sample that doesn't also satisfy the $x_n \leq c_n$ constraints.

• This is actually working quite well now, since drawing from the dirichlet already ensures one of my constraints is met. As an addition, does anyone know if it is possible to scale the dirichlet samples so that they conform more to the hard constraints given above in each dimension. As an example, if one constraint is a lot tighter than the others then there will be mostly rejects when sampling... – user1 Sep 27 '11 at 2:00
• @user1 So, is it possible to scale the Dirichlet samples so the rejection step is no longer needed? – Tomilov Anatoliy Jul 17 '14 at 9:15

I believe the primary method to generate uniform samples in convex polytopes is via a Markov-chain random walk. For example, the paper by Ravi Kannan and Hariharan Narayanan, "Random walks on polytopes and an affine interior point method for linear programming," STOC 2009, (ACM link) achieves a strongly polynomial mixing time. The basic idea is explained in these older notes by H.E. Romeijn and R.L Smith, "Sampling through Random Walks," (PDF), or, perhaps better, Santosh Vempala's 2008 notes on Algorithmic Convex Geometry (PDF).

The Kannan-Narayanan paper gives a condensed history of algorithms to find (approximately) the volume of a convex polytope, starting from the Dyer, Frieze, and Kannan paper from 1991 whose mixing time was $O(n^{23})$ for a polytope in $\mathbb{R}^n$, through a series of steady improvements reducing the dependency on $n$ down to closer to $n^2$. I am not certain this is the latest word on this topic.

• Is there a simple algorithm for small-dimensional convex polytopes? MCMC seems overkill in something like 5 dimensions. – becko Jan 11 '18 at 17:10

If the $c_n$ are sufficiently small (say, all smaller than $1/\sqrt{n}$) or sufficiently large (say, all at least $1$) the polytope is a simplex, so estimating the volume and sampling are both easy as pointed out in @Mike's answer. On the other hand, trying to find a local optimum by uniform sampling is doomed in $20$-dimensional space, since the distance between the samples will be very large for any reasonable number of samples. (for example, if your polytope were a cube of side one, generating $2^{40}$ points will have inter point distances around $1/4.$ Needless to say, generating $2^{40}$ points might take a while.