How to pick a random points in higher dimensional parametric regions?

Question

I am using a simple code to pick points randomly (uniformly) from a parametric region. This seems to work perfectly with 2D and 3D regions. For instance

a = ParametricRegion[{Cos[x],Sin[x]Cos[y],Sin[x]Sin[y]}, {{x, 0, \[Pi]}, {y, 0, 2 \[Pi]}}];
pts = RandomPoint[a, 5000];
Graphics3D[{Darker[Blue], PointSize[Small], Point[pts]},Boxed -> False]

I would like to do the same with other regions in higher dimensions. In 4D this process works but is slower (slightly). But as soon as I go to 5D, it seems to fail. With even as low as 50 points it takes too long (I terminated it after about 15 minutes).

Is there a way around this? Alternately, is there another way to pick points uniformly from a parametric region that will be faster/more efficient?

Edit: I have looked at the first answer given here but this uses DiscretizeRegion which will not work for dimensions greater than 3.

Answer 1

Let's consider a parameterization $f \colon \varOmega \to \mathbb{R}^d$ with domain $\varOmega \subset \mathbb{R}^n$. Let's denote the surface measure on $\varSigma = f(\varOmega)$ by $\omega$ so that $\mu = \omega / \int_\varSigma \operatorname{d} \omega$ is the uniform probability measure on $\varSigma$.

Assuming that random points in the domain $\varOmega$ of the parameterization can be obtained quite cheaply, we can pull the uniform density on the surface back to the domain. More precisely, we are looking for a function $\varrho \colon \varOmega \to \mathbb{R}$ so that

$$\int_A \operatorname{d} \omega = \int_{f^{-1}(A)} \, \varrho \, \operatorname{d} \lambda \quad \text{for all measurable $A \subset \varSigma$,}$$

where $\lambda$ denotes the Lebesgue measure on $\varOmega \subset \mathbb{R}^n$

One says $f^* \omega(B) = \int_B \varrho \, \operatorname{d} \lambda$ is the pullback measure of $\omega$. Unfortunately, Wikipedia has only an entry for the pushforward measure. But since $f$ is one-to-one between $\varOmega$ and $\varSigma$, pushforward and pullback are inverse of each other. Likewise, $(f^* \omega) / \int_\varOmega \operatorname{d} (f^* \omega)= f^* \mu$ is the pullback measure of the probability measure $\mu$ that we are looking for. So instead of sampling from $\mu$ directly, we can also sample from $f^* \mu$ and map the sampled point cloud back to $\varSigma$ by applying $f$.

All this would be rather useless if we hadn't a simple formula for $\varrho$:

$$\varrho(x) = \sqrt{\det ( Df(x)^\intercal Df(x))} \quad \text{for $x \in \varOmega$.}$$

But how do we sample from $f^* \mu$? Well, we are going to use a Monte-Carlo approach. Let's assume that $\varrho$ is bounded and define $m = \max_{x \in \varOmega} \varrho(x)$. We sample uniformly from $\varOmega \times [0,m]$ and discard points $(x ,t) \in \varOmega \times [0,m]$ with $t > \varrho(x)$. Then the chance for an accepted point $(x,t)$ to have in $x \in B$ is proportional to the volume of the set $\{(x,t) | \text{$x \in B$ and $\varrho(x) \leq t$}\}$ and this volume is equal to $\int_B \varrho \, \operatorname{d} \lambda$. Written down with conditional probablities, this reads as

$$ \begin{align*} P(x \in B \mid \text{$(x,t)$ accepted}) &= \frac{\text{$P(x \in B$ and $(x,t)$ accepted})}{P(\text{$(x,t)$ accepted})} = \frac{\int_B \varrho \, \operatorname{d} \lambda}{\int_\varOmega \varrho \, \operatorname{d} \lambda} = \int_B \operatorname{d} (f^* \mu) \end{align*} $$

I'll write down the algorithm for the sphere example. f is the parameterization, ρ the pullback of the density (it is not a probability density!)

f = X \[Function] With[{x = X[[1]], y = X[[2]]}, {Cos[x], Sin[x] Cos[y], Sin[x] Sin[y]}];
dim = 2;
ρ = X \[Function] 
   With[{Df = D[f[Table[X[[i]], {i, 1, dim}]], {Table[X[[i]], {i, 1, dim}], 1}]},
    Simplify[Sqrt[Det[Df\[Transpose].Df]]]
    ];

Specify the domain and compute the maximum of ρ on the domain.

domain = Rectangle[{0, 0}, {Pi, 2 Pi}];
m = NMaximize[ρ[Array[x, dim]], Array[x, dim] ∈ domain][[1]];

The idea of the algorithm is to compute random points (uniformly distributed) in the domain and for each point a uniformly distributed random number between $0$ and $m$, where $m$ is the maximum of $\varrho$ on the domain. We will accept a point $x$ if and only if its number $t$ satisfies $t \leq \varrho(x)$. In the end we map the accepted $x$ onto the surface by $f$. In order to speed this up, we can compile the density $\varrho$ and the parameterization $f$.

Block[{XX, X},
  XX = Table[Compile`GetElement[X, i], {i, 1, dim}];
  cρ = With[{density = ρ[XX]},
    Compile[{{X, _Real, 1}},
     density,
     CompilationTarget -> "C",
     RuntimeAttributes -> {Listable},
     Parallelization -> True,
     RuntimeOptions -> "Speed"
     ]];

  cf = With[{code = f[XX]},
    Compile[{{X, _Real, 1}},
     code,
     CompilationTarget -> "C",
     RuntimeAttributes -> {Listable},
     Parallelization -> True,
     RuntimeOptions -> "Speed"
     ]];
  ];

Finally, this produces the random points on the sphere:

pointcount = 2500000;
rand = cf@With[{pts = RandomPoint[domain, pointcount]},
      Pick[pts, 
       UnitStep[
        Subtract[cρ[pts], RandomReal[{0, m}, pointcount]]], 1]
      ]; // RepeatedTiming // First
rand // Dimensions

0.306

{1592292, 3}

So this returns about 1.6 million random points on the sphere in the same time that RandomPoint[a, 5000] needed to compute 5000 random points.