1
$\begingroup$

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.

Improve this question
$\endgroup$
1
  • $\begingroup$ Maybe try converting the intersection of all half-spaces into a single ConicHullRegion first? $\endgroup$
    – Silvia
    Jan 19 at 16:33

1 Answer 1

Reset to default
0
$\begingroup$

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.

Improve this answer
$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – hilberts_drinking_problem
    Aug 27, 2021 at 4:58
  • $\begingroup$ 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. $\endgroup$
    – David G. Stork
    Aug 27, 2021 at 4:59
  • $\begingroup$ 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. $\endgroup$
    – hilberts_drinking_problem
    Aug 27, 2021 at 5:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.