# Sampling uniformly from unit ball around a point on a Riemannian manifold

I'm interested in a mechanical procedure (i.e., algorithm) for sampling uniformly from $$\{y \in \mathcal{M} : d_\mathcal{M}(x, y) \le 1\}$$ for some fixed $$x \in \mathcal{M}$$, using well-known tools. The measure is the one induced by the corresponding Riemannian metric $$\{g_x\}_{x \in \mathcal{M}}$$.

When does such an algorithm exist and when is it non-trivial to find one?

We could start from the following schema:

• Sample a tangent vector uniformly from $$\{ u \in \mathcal{T}_x \mathcal{M} : \lVert u \rVert_x \le 1 \}$$.

• $$y = \textrm{exp}_x(u)$$, where $$\textrm{exp}(\cdot)$$ is the exponential map.

Clearly the first step is non-trivial, but that aside, would this approach yield uniform samples in the unit ball around $$x$$?

And then, is it possible to specify a class of manifolds for which it is possible to automate the first step? Perhaps for embedded submanifolds where we could sample uniformly from the ambient space and then somehow restrict to $$\mathcal{T}_x \mathcal{M}$$ (and apply the inverse Riemannian metric if different than the inherited one)?

Take for instance the Lorentz model of hyperbolic geometry: $$\mathcal{H}^n = \{ x \in \mathbb{R}^{n+1} : \langle x, x \rangle_\mathcal{L} = -1, x_0 > 0 \},\quad\textrm{with}\ \langle x, y \rangle_\mathcal{L} = -x_0 y_0 + \sum_{i=1}^n x_i y_i.$$ The tangent space is $$\mathcal{T}_x \mathcal{H}^n = \{ v \in \mathbb{R}^{n+1} : \langle x, v \rangle_\mathcal{L} = 0 \}$$. How could one sample from the unit ball of tangent vectors uniformly?

Your proposed scheme will be ruined by curvature. For example, consider a 2-sphere of radius $$1/\pi,$$ so that the disc $$\exp(B_1) = \{\exp(u) : \|u\| \le 1\}$$ perfectly covers it. The outer region $$0.5 \le \|u\| \le 1$$ of the unit disc in the tangent space contains $$3/4$$ of the measure, so samples drawn using the uniform distribution on $$B_1$$ would hit this region $$3/4$$ of the time; but this corresponds to only half of the sphere in the Riemannian measure. (In your hyperbolic example you would find the opposite - the outer regions have much more Riemannian measure than "exponential" measure.)

If we use the exponential map on $$B_1$$ as a coordinate system $$x$$, this is reflected in the equation $$d\mu_g = \sqrt{\det g}\ d^nx. \tag 1$$ Here $$d\mu_g$$ is the Riemannian measure associated to $$g,$$ while $$d^nx$$ is the standard/Euclidean/uniform measure on the tangent/coordinate space. Thus, unless you happen to have $$\det g=1$$ in these coordinates, the measures are not the same.

However, equation $$(1)$$ should tell you how you can sample from $$d\mu_g:$$ if you have an expression for $$g$$ in normal coordinates, then $$w=\sqrt{\det g}$$ is a scalar function you can compute, so we have reduced the problem to that of sampling a weighted measure $$d\mu_g = w\ d^nx,$$ which you should be able to do programmatically (e.g. via the rejection method).

(One caveat: this only works if the injectivity radius at $$x$$ is at least $$1$$; i.e. if the exponential map is injective on $$B_1.$$ Otherwise, different points in $$B_1$$ will map to overlapping points on the manifold, so you'd have to restrict to a subset of $$B_1$$ making the exponential map cover each point in the unit ball exactly once.)

• Thanks, that's very helpful! I was wondering about your brief comment on the hyperbolic space: the Lorentz model is a semi-Riemannian manifold with constant metric $g = \textrm{diag}(-1, 1, \ldots, 1)$, so $|\det g| = 1$. Doesn't this mean that the above sampling will work as it is? Sep 9 '19 at 14:46
• Actually, I'm not sure if that makes sense because that metric tensor is the $(n + 1)$-dimensional one used in the Minkowski space, not the $n$-dimensional one at $x \in \mathbb{H}^n$ once a coordinate system is fixed. Sep 9 '19 at 15:21
• @Călin: your second comment is exactly correct - this method works only when using the "intrinsic" coordinates on the $n$-manifold of interest, not the "extrinsic" coordinates of any $>n$-dimensional space it's embedded in. (Note that it does, however, work for any intrinsic coordinates, not just geodesic normal coordinates - so if you a coordinate system in which you can calculate $\det g,$ you should be able to do some "automated" sampling.) Sep 9 '19 at 22:55
• Could you please just expand a little on how rejection sampling would be used here? If I understand correctly, the same relationship applies between the Riemannian measure and any other measure (not just the uniform one, $\mathbb{1}_\mathcal{X} (x) / \mathcal{Vol}(\mathcal{X})$). So if I want to sample a pdf $p(x)$ defined on the manifold, I can sample instead $p(x) \sqrt{\det g(x)}$ in some exponential chart (assuming injectivity). But how exactly would I use rejection sampling here? I first need to find a function that upper bounds it, right? Sep 10 '19 at 12:30
• @Călin: To sample with the weight $f(x)=p(x)\sqrt{\det g(x)},$ first normalize $\tilde f = C f$ so that its maximum is less than $1$ (ideally equal to $1$), then generate random points $x$ in the chart coordinates and reject them with probability $1-\tilde f(x).$ Sep 11 '19 at 0:43