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

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.

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[], y = X[]}, {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][];

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.

• can you elaborate a little on the algorithm? I can't make out what this pull back is. If this is something standard then please suggest a reference. – nomaan x Apr 17 '18 at 2:27