How to sample uniformly points around a neighborhood of a point lying on a n-sphere?

Question

Given a point $x$ lying on the surface of a n-sphere $S$, what is an efficient way of randomly sampling points $x_k \in S$ such that their distance from $x$ is at most $r$? ($\|x-x_k\| < r$)

Can this be done uniformly?

I have efficient ways of converting cartesian coordinates to n-spherical coordinates and back.

Answer 1

Sample in a neighborhood of $e_{n+1}=(0,0,\ldots,0,1)$ in $\mathbb{R}^{n+1}$ and then apply any orthogonal transformation (that is, isometry of the sphere) that sends $e_{n+1}$ to $x_k$. This reduces the problem to sampling around $e_{n+1}$.

Geometry shows that the last coordinate of the points, $Z$, will range from $1$ down to $1 - r^2/2$. The distribution of $Z$ is otherwise that of a Beta$(n/2,n/2)$ variate, multiplied by $2$ and shifted down by $1$, as shown at https://stats.stackexchange.com/a/85977/919. Conditional on $Z$, the other $n$ coordinates will lie on a sphere of dimension $n-1$ and radius $\sqrt{1-Z^2}$. Generate these using any of the methods described at How to generate uniformly distributed points on the surface of the 3-d unit sphere? (or otherwise).

To draw values of $Z$, compute the quantile associated with $1-r^2/2$. If $F$ is the Beta$(d/2,d/2)$ distribution function, this will be $F^{-1}(r^2/4)$. For $U$ a uniform variate in $[0, F^{-1}(r^2/4)]$, set $Z=1-2*F(U)$.

Here is R code to illustrate the ideas (and fill in any details I might not have explained sufficiently clearly). Tests of $S^1$ and $S^2$ (which can be directly visualized) and of higher-dimensional spheres with very small distances and the maximum distance ($2$) are consistent with what one would expect, suggesting it is working correctly.

d <- 2       # Dimension of the sphere; one less than dimension of its Euclidean space.
n <- 1e3     # Number of points to generate.
r.max <- 0.2 # Maximum Euclidean distance.
#
# Generate uniformly random points on a `dim`-sphere of radius `radius`.
# Returns (dim+1)-dimensional vectors as rows.
#
rsphere <- function(n, dim, radius=1) {
  x <- matrix(rnorm((dim+1)*n), nrow=n)
  x / (sqrt(rowSums(x^2)) / radius)
}
#
# Generate random heights on the d-sphere.
#
q <- pbeta(r.max^2 / 4, d/2, d/2)           # Limiting quantile
z <- 1 - 2*qbeta(runif(n, 0, q), d/2, d/2)  # Last coordinate
#
# Compute the corresponding radii of the cross-sections at heights `z`.
#
rho <- sqrt(1 - z^2)                        # Radius
#
# Generate the remaining (first `d`) coordinates uniformly.
# Results are in the rows of `x`.
#
x <- cbind(rsphere(n, d-1, radius=rho), z)

Incidentally, it's simple and fast to find an orthogonal transformation of $\mathbb{R}^{n+1}$ that sends $e_{n+1}$ to $x_k$: choose the reflection with this property. Here is sample R code to apply such a reflection to an entire array of coordinates, such as the x produced above.

#
# Reflect points `x` (array of rows) in a way that sends (0,0,..,0,1) to `target`.
#
reflect <- function(x, target) {
  if (!is.matrix(x)) x <- matrix(x, nrow=1)
  n <- length(target)
  v <- c(rep(0, n-1), 1) - target
  v.norm2 <- sum(v^2)
  return(x - outer(2/v.norm2 * c(x%*%v), v))
}

Answer 2

Adding to @whuber's answer here is python code to reproduce the ideas he explained, as well to validate the results.

from scipy.stats import beta
import numpy as np

def rsphere(n,n_dim,rad=1.0):
    X = np.random.normal(size=(n,n_dim+1))
    X_norm=np.divide(np.linalg.norm(X,axis=1),rad)
    X = np.divide(X.T,X_norm).T
    return X

def reflect(X, target):
    n_dim = X.shape[1]
    pole = np.array([0]*(n_dim-1) + [1])
    v = pole - target
    X_new = X - np.outer(2.0/np.dot(v,v) * np.matmul(X,v).T, v)
    return X_new


r_max=0.5 # max radius to sample
n_dim=200 # number of dimension of the space
n_sphere=n_dim-1 # dimension of the n-sphere
n = int(1e3) # numer of points to sample

target = np.zeros(n_dim) # target
target[-2]=1 
pole = np.array([0]*(n_dim-1) + [1]) # north pole

# generate points
X = rsphere(n,n_sphere,rad=1.)
rads = np.random.uniform(0.0,r_max,n)
q = beta.cdf(np.square(rads)/4.0, n_sphere/2.0, n_sphere/2.0)
unif = np.random.uniform(0,q,n)
z = 1 - 2*beta.ppf(unif, n_sphere/2.0, n_sphere/2.0)
rho = np.sqrt(1.0 - np.square(z))
X_pole = np.hstack((rsphere(n, n_sphere-1, rad=rho),z.reshape(-1,1)))
X_nei = reflect(X_pole, target)

To validate results we can look at a histogram of distances of the generated points vs the target:

import matplotlib.pyplot as plt
X_dist = np.linalg.norm(X_nei-target,axis=1)
plt.hist(X_dist)
plt.xlabel('$\|z_i-z_{target} \|$')
plt.show()

And to visualize the 3D case:

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.set_aspect('equal')
# Create a sphere
r = 1
pi = np.pi
cos = np.cos
sin = np.sin
phi, theta = np.mgrid[0.0:pi:40j, 0.0:2.0*pi:40j]
x = r*sin(phi)*cos(theta)
y = r*sin(phi)*sin(theta)
z = r*cos(phi)
# sphere
ax.plot_surface(
    x, y, z,  rstride=1, cstride=1, alpha=0.25, linewidth=0)
# pole
ax.scatter(X_pole[:,0],X_pole[:,1],X_pole[:,2],c='b',alpha=0.5)
ax.scatter(pole[0],pole[1],pole[2],s=40,c='k')
# target
ax.scatter(X_nei[:,0],X_nei[:,1],X_nei[:,2],c='r',alpha=0.5)
ax.scatter(target[0],target[1],target[2],s=40,c='k')
plt.show()