Generating random vectors with constraints

I need to create random vectors of real numbers a_i satisfying the following constraints:

abs(a_i) < c_i;
sum(a_i)< A;        # sum of elements smaller than A
sum(b_i * a_i) < B; # weighted sum is smaller than B
aT*A*a < D          # quadratic multiplication with A smaller than D

where c_i, b_i, A, B, D are constants.


What would be the typical algorithm to generate efficiently this kind of vector?

• What do you mean by the fourth constraint, "The magnitude of a is .." – M. Tibbits Mar 11 '11 at 15:53
• My mistake. Finished description. Thanks for the feedback. – LouisChiffre Mar 11 '11 at 16:18
• How is that a_i follows distribution p_i and is also less that c? That's because the distribution p_i is also less that c? In which distribution are you thinking? – deps_stats Mar 12 '11 at 0:31
• @deps_stats. Very good points. Pseudo code was not very clear. The distribution I have in mind is the poisson distribution. Each element follows a poisson distribution with a different lambda. With that in mind I guess the first condition (a_i < c) is not necessary since I can just rescale a_i at the end of the generation to satisfy it. – LouisChiffre Mar 12 '11 at 1:43
• Let me ask something else... Are the c,A,B and lambdas fixed? – deps_stats Mar 15 '11 at 0:05

If I understand you correctly, only points in some small volume of n-dimensional space meet your constraints.

Your first constraint constrains it to the interior of a hypersphere, which reminds me of the comp.graphics.algorithms FAQ "Uniform random points on sphere" and How to generate uniformly distributed points in the 3-d unit ball? The second constraint slices a bit out of the hypersphere, and the other constraints further whittle away at the volume that meets your constraints.

I think the simplest thing to do is one of the approaches suggested by the FAQ:

• choose some arbitrary Axis-aligned bounding box that we are sure contains the entire volume. In this case, -c < a_1 < c, -c < a_2 < c, ... -c < a_n < c contains the entire constrained volume, since it contains the hypersphere described by the first constraint, and the other constraints keep whittling away at that volume.
• The algorithm uniformly picks points throughout that bounding box. In this case, the algorithm independently sets each coordinate of a candidate vector to some independent uniformly distributed random number from -c to +c. (I am assuming you want points distributed with equal density throughout this volume. I suppose you could make the algorithm select some or all coordinates with a Poisson distribution or some other non-uniform distribution, if you had some reason to do that).
• Once you have a candidate vector, check each constraint. If it fails any of them, go back and pick another point.
• Once you have a candidate vector, store it somewhere for later use.
• If you don't have enough stored vectors, go back and try to generate another one.

With sufficiently high-quality random number generator, this gives you a set of stored coordinates that meet your criteria with (expected) uniform density.

Alas, if your have a relatively high dimensionality n (i.e., if you construct each vectors out of a relatively long list of coordinates), the inscribed sphere (much less your whittled-down volume) has a surprisingly small part of the total volume of the total bounding box, so it might need to execute many iterations, most of them generating rejected points outside your constrained area, before finding a point inside your constrained area. Since computers these days are pretty fast, will that be fast enough?

• So what you're suggesting is effectively to sample the space. I've got a similar problem except finding a bounding box cannot be done statically (IE, cannot be hard-coded). From experience, this breaks down if your constraints are of the form f1(x1) + f2(x2) == C Any suggestions here? – Groostav May 31 '15 at 4:02
• Yes, the sampling method doesn't work if you have equality ("==") constraints. Constraints such as points that are on the surface of a sphere, or on the surface of a cylinder, etc. (radius == R), rather than inside the sphere (radius <= R). Uniformly selecting points over the whole volume will "never" (probability close to 0) hit your desired surface. So you need to somehow only pick points that are on that surface -- i.e., to find points [x1, x2, x3] such that f1(x1) + f2(x2) == C, you could randomly pick x1 and then force x2 = inverse_f2( C - f1(x1) ). – David Cary May 31 '15 at 15:16
• For the special case of uniformly distributed points on the surface of a sphere, see "Uniform random points on sphere". – David Cary May 31 '15 at 15:18
• @Groostav: Perhaps your question is different enough from the original question that you could post it as a new top-level question? "I've just been told I have to post a follow-up question, why and how?" – David Cary May 31 '15 at 15:49