# RandomPoint does not return with many region intersections

I am trying to generate random points on an n-dimensional unit sphere that satisfy a number of linear inequalities (a point $$v \in \mathbb{S}^N$$ has to satisfy $$v \cdot x \geq v \cdot x_k$$ for all $$x_k$$ in some list of vectors; the set of such v's is non-empty). The distribution has to be uniform on the region defined by the linear constraints. The following code defines a region to sample as an intersection of a bunch of half-spaces and a sphere. It works as expected when the number of half-spaces (numPts - 1) is small, and stops returning when the number of half-spaces gets large. Since there is no error message, I am having trouble diagnosing the problem. I am on version 12.1.0.0.

Clear["Global*"];
SeedRandom[0];
numDims = 15;
numPts = 200; (*works*)
(*numPts = 700; (*does not work*) *)
y = RandomReal[1, numDims];
xs = RandomPoint[Simplex[numDims], numPts];
vs = Dot[xs, y];
j = Ordering[vs, -1][[1]];
x = xs[[j]];
rest = Delete[xs, j];
hs = Table[HalfSpace[rest[[k]] - x, 0], {k, 1, Length[rest]}];
reg = Fold[RegionIntersection, Region[Sphere[numDims]],  hs];
samplePts = RandomPoint[reg, 10]


I saw a recent related post that suggested using DiscretizeRegion, but I could not get that to work. I also tried using ImplicitRegion without success.

I am trying to get the sampler to work reliably for around 1000 half-spaces.

• Maybe try converting the intersection of all half-spaces into a single ConicHullRegion first? Jan 19 at 16:33

Try this (which works in 0.1 second):

Clear[numNeeded, dimension, thevecs, thevertices, theWeights,
numConstraints];
numNeeded = 50;
dimension = 15;
numConstraints = 1000;
thevecs =
Normalize /@ RandomReal[{-1, 1}, {numConstraints, dimension}];
thevertices =
Flatten[Table[Drop[thevecs, {i}], {i, Length[thevecs]}], 1];
theWeights =
Normalize /@ Table[RandomReal[1, Length[thevertices]], numNeeded];
finalPoints =
Normalize /@ Table[theWeights[[i]].thevertices, {i, numNeeded}]


How I (tried) to make this work:

• numNeeded is the number of random points you seek
• dimension is the dimension of the sphere
• thevecs are the many constraint vectors you have (which might be 1000)
• theWeights are a set of relative weightings of the vectors at the vertices of the simplex defined by the constraint vectors missing one other vector
• The weights multiply the vertices to find a random point within the simplex. There are numNeeded of these.

I may have slipped up, but this is the general approach.

• Thanks. My understanding is that RandomPoint` is supposed to result in a uniform distribution in the region. If I understand correctly, this method does not guarantee a uniform distribution. Aug 27, 2021 at 4:58
• Indeed: This gives a point in the center of the region so defined. By the way, if you have 1000 constraints, you'll need $2^{1000}$ sample points in order to get just a single random point that satisfies your constraints. That's a lot. Aug 27, 2021 at 4:59
• To give some background, I am able to sample from a uniform distribution of points meeting the constraints when the points lie in a unit simplex (instead of a sphere). This can be done with a hit-and-run sampler and does not require computing $2^{1000}$ vertices of the region. My main interest in using Mathematica is to see if there is a way to do the same on a sphere instead of a polytope. Aug 27, 2021 at 5:23