2
$\begingroup$

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]

enter image description here

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.

$\endgroup$
2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – nomaan x Apr 17 '18 at 2:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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