Sampling from intersection of sphere and simplex

Question

I need to isotropically sample vectors in $\mathbb{R}^d$ of Euclidean norm $c$ such that components are positive and add up to 1. Any suggestions how to do this in Mathematica?

For $d=3$ this means sampling from one of the circular contours in

SliceContourPlot3D[Sqrt[x^2 + y^2 + z^2], 
 x + y + z == 1, {x, 0, 1}, {y, 0, 1}, {z, 0, 1}]

Things I tried

1. RandomPoint[ RegionIntersection[Sphere[{0, 0, 0}, .9], Simplex[IdentityMatrix[3]]]] -- returns unevaluated
2. Rejection sampling. Inefficient for $d=3$, completely impractical for higher values, I need $d=100$.

Code below is for debugging rejection sampling in $d=3$ case. It takes a sample of 3d points, projects onto $x_1+x_2+x_3=1$ plane, displays them along with the target region marked in gray. Will appreciate any suggestions on making an efficient sampler.

(* simplex visualization code from https://mathematica.stackexchange.com/a/139443/217 *)
visualize[points_, c_] := Module[{},
   d = 3;
   kk = Array[k, 6];
   mat = Partition[kk, 2];
   bb = Array[b, 3];
   source = {{-1/Sqrt[2], 0}, {1/Sqrt[2], 0}, {0, Sqrt[3/2]}};
   target = IdentityMatrix[3];
   eqs = Table[mat . source[[i]] + bb == target[[i]], {i, 3}]; 
   sol = First@Solve[eqs, kk~Join~bb];
   {mat0, bb0} = {mat, bb} /. sol;
   imat0 = PseudoInverse[mat0];
   unmap[point_] := imat0 . (point - bb0);
   expr = mat0 . {x, y} + bb0;
   reg = ImplicitRegion[Reduce[Thread[expr > 0]], {x, y}];
   
   regionPlot = RegionPlot[reg, PlotStyle -> None];
   pointsPlot = ListPlot[unmap /@ points];
   val = Sqrt[Total[expr^2]];
   contourPlot = 
    ContourPlot[val == c, {x, y} \[Element] reg, 
     AspectRatio -> Automatic, PlotPoints -> 10, 
     ContourShading -> None, 
     ContourStyle -> Directive[Opacity[.3], Gray, Thickness[.02]]];
   
   Show[regionPlot, contourPlot, pointsPlot]
   ];

points = RandomVariate[NormalDistribution[], {100, 3}];
visualize[points, .8]

Motivation is to get a diverse sampling of discrete distributions of given complexity. JimB's previous answer achieves some progress for the case when entropy is the complexity measure. This question corresponds to linear entropy as the complexity measure. Related question on math.SE.

Answer 1

Try Ball & InfinitePlane instead of Sphere & Simplex in your first approach:

reg=DiscretizeRegion@  RegionIntersection[Ball[{0, 0, 0}, .9], InfinitePlane[IdentityMatrix[3]]]
p=RandomPoint[reg , 100];
Show[Region[reg], Graphics3D[Point[p]]]

Answer 2

First two remarks:

• The intersection between the simplex given by $x_1 + \ldots + x_d = 1$ and $x_1,\ldots,x_d \geq 0$ and the sphere $(x_1)^2 + \ldots + (x_d)^2 = c^2$ is nonempty if and only if $1/\sqrt{d} \leq c \leq 1$.
• If $1/\sqrt{2} < c < 1$, which seems to be of interest to OP, the intersection is given by $d$ disconnected, congruent "caps". This answer probably makes most sense in that case. I have not actually checked this in detail, let me know if I got this wrong.

An intuitive idea is to randomly pick one corner of the simplex (w below) and then randomly one point on the opposite face (v below) and then pick the convex combination of the two (s*v+(1-s)*w below) that has the right norm c. For the choice of v I now use DirichletDistribution which I learned about from the answer of @JimB but note that I use it for a face of the simplex, not the simplex itself.

I have to warn that the probability measure that I get is not the one OP expects. Intuitively it is not extremely far away either, but this is a low dimensional intuition. It is possible that this code can be improved by picking v better. Here is the code:

randomsimplex[d_]:=RandomVariate[DirichletDistribution[
                      ConstantArray[1,d]]]//Join[#,{1-Total[#]}]&;
random[c_,d_]:=With[{i=RandomInteger[{1,d}]},
  With[{v=Insert[randomsimplex[d-1],0.,i],w=UnitVector[d,i]},
    With[{s=(1-c^2)/(1+Sqrt[c^2+(-1+c^2)*Dot[v,v]])},
      s*v+(1-s)*w]]];

Beware that for $c$ near the minimum it can happen that a complex solution is produced. But for $c \geq 1/\sqrt{2}$ this cannot happen. No resampling is used, and it also works in very high dimensions:

random[0.8,10000]
(* takes about 0.2 seconds *)

To check that it satisfies all conditions, one can use

test[c_,d_]:=With[{r=random[c,d]},
  And[Chop[Total[r]-1]===0,Chop[Norm[r]-c]===0,And@@Thread[r>=0]]];

For example

test[0.8,10000]
(* True *)

Plot. Here is a plot for $d=4$ and the threshold value $c=1/\sqrt{2}$. The plot is internal to $x_1+x_2+x_3+x_4=1$, so the caps are 2-dimensional:

Code for the plot:

With[{X=Orthogonalize[Join[{{1,1,1,1}},IdentityMatrix[4][[1;;3]]]][[2;;4]]},
ListPointPlot3D[Table[random[1/Sqrt[2],4],10000].Transpose[X],BoxRatios->Automatic,Ticks->None]]

Note. The speed of this code is limited by RandomVariate[DirichletDistribution[...]]. To make it faster, generate many variates at once using RandomVariate[DirichletDistribution[...],n].

Answer 3

In my opinion unnecessarily convoluted but still:

RandomPoint[
  ImplicitRegion[
   RegionMember[
    RegionIntersection[Sphere[{0, 0, 0}, .9], 
     Simplex[IdentityMatrix[3]]], {x, y, z}], {x, y, z}], 1000] //
 Graphics3D[
   {Opacity[9/10], Sphere[{0, 0, 0}, .9], Simplex[IdentityMatrix[3]], 
    Point[#]},
   ViewPoint -> {1, 2, 1}] &

Answer 4

Some progress on an approach to address rejection sampling efficiency -- start at the "extreme point" (as much mass concentrated on one coordinate as possible), apply small random rotation, then project back onto the feasible region...repeat until the resulting distribution is close to isotropic (this part needs work)

randomRotate[v_, theta_] := (
   v2 = RandomVariate[NormalDistribution[], d];
   RotationMatrix[theta, {v, v2}] . v);
c = .8;
d = 3;
v = {xx, (1 - xx)/2, (1 - xx)/2};
v0 = v /. First@Solve[{Total[v*v] == c^2, 0 < xx < 1}, xx];
y = ConstantArray[1/d, d];
ii = IdentityMatrix[d];
projector = (ii - PseudoInverse[{y}] . {y});
planeProject[x_] := x . projector + y;
sphereProject[x_] := c*Normalize[x];

theta = Pi/40;
points = {};
v = v0;
For[i = 1, i < 100, i += 1,
  AppendTo[points, v];
  v = randomRotate[v, theta];
  AppendTo[points, v];
  v = Nest[sphereProject[planeProject[#]] &, v, 10];
  AppendTo[points, v];
  ];