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.