Random point on hypersphere surface at a uniform distance of another point on the surface

Question

Suppose I have an n-sphere of radius 1 centered in $(0,0,...,0)$, where each point on the surface represents a multinomial distribution:

• Given coordinates of a point $S=(x_1, x_2,...,x_n)$ on the sphere surface, with $-1 \leq x_i \leq 1$, the associated distribution would be $P$ with parameters $(x_1^2, x_2^2, ..., x_n^2)=(p_1,p_2,...,p_n)$, where $p_i$ represents the probability of value $i$.

I would like to obtain, given $S$ on the sphere surface, a random point $S'$ on the sphere surface such that the euclidean distance $d(S,S')$ follows a uniform distribution in the range $(0, D)$, which would be equivalent of a method of adding noise to the multinomial distribution $P$ associated to $S$ such that $P'$ is also a correct multinomial distribution.

The final goal is to use this method in a simulated annealing program to optimize the parameters of a multinomial distribution, as a way to find a neighbor of the current estimation of the probability distribution.

Is there a feasible, perhaps approximate way to find $S'$ in a reasonable computation time for an $n$-sphere with $n \ge 3$?

Answer 1

I ended up finding what I think is an exact solution.

Let $\sigma$ a plane that passes through the origin $O$, where the normal vector is described by the point $S$. The neighbors of $S$ are in the surface of a spherical cap with a base parallel to $\sigma$.

A random unit vector $V$ parralel to $\sigma$ may be obtained by projecting a random point in the sphere surface onto $\sigma$ and normalizing.

If we take $d$ as the minimum distance from $S$ to the center of the sperical cap base $B$, we can determine $||OB||=1-d$ and the maximum distance to a neighbor $S'$ of $S$ (in the spherical cap) as $||BS'||=\sqrt{1-(1-d)^2}$ with Pythagoras theorem. The distance of $S$ to $S'$ can also be determined with Pythagoras theorem: $||SS'||=\sqrt{2d}$.

Taking last two paragraphs into account, a random neighbor $S'$ of $S$ at a distance $||SS'||=D$ can be obtained as $OB+V||BS'||$, as all the unknowns can be obtained with $d$ which can be obtained by $d=D^2/2$. If a value $D'$ is chosen from a random uniform in $[0,D]$, and then $S'$ is obtained for that $D'$, then $d(S,S')$ will follow a uniform distribution in $[0,D]$. Here goes a graphic representation, (you have to imagine a sphere).

Here goes some Python code for those interested:

import numpy as np
from numpy.linalg import norm
from mayavi import mlab

# https://math.stackexchange.com/questions/1585975/how-to-generate-random-points-on-a-sphere
def random_point_in_sphere(n):
    vec = np.random.randn(n)
    vec /= norm(vec)
    return vec

# https://stackoverflow.com/questions/9605556/how-to-project-a-point-onto-a-plane-in-3d
def random_vector_on_plane(normal_vec):
    # Get a random point in sphere and project it to the plane
    P = random_point_in_sphere(normal_vec.shape[0])
    dist = np.dot(P, normal_vec)  # distance from P to plane along the normal
    vec = P - dist*normal_vec
    vec /= norm(vec)
    return vec

# Generate neighbor S' of S in the same sphere surface such that the euclidean 
# distance d(S,S') follows a uniform distribution in [0, d]
def generate_neighbor(S, d):
    # a value of d describes a spherical cap, its height is d**2/2
    cap_height = np.random.uniform(0, d)**2 / 2
    # Generate a direction on plane parallel to spherical cap base
    D = random_vector_on_plane(S)
    # Point on shperical cap perimeter built as sum of two vectors:
    # 1. Vector from origin to S that stops at cap base
    # 2. Random vector parallel to cap base. Its length can be computed with
    # the use of pythagoras theorem.
    neighbor = S * (1-cap_height) + D*np.sqrt(1-(1-cap_height)**2)
    return neighbor


# Test that the generate neighbor function works correctly
n = 10  # 10-sphere
n_exp = 50000
n_points = 50
d = 0.5
test_points = np.array([random_point_in_sphere(n) for _ in range(n_points)])

means_distances = []
std_devs_distances = []
max_distances = []
for S in test_points:
    distances = []
    for _ in range(n_exp):
        neighbor = generate_neighbor(S, d)
        distances.append(np.linalg.norm(S-neighbor, ord=2))  # eucl. distance
    means_distances.append(np.mean(distances))
    std_devs_distances.append(np.std(distances))
    max_distances.append(np.max(distances))


print(f'Mean of euclidean distance: {np.mean(means_distances)}'
      f' expected: {d / 2}')
print(f'Std.dev. of mean across experiments (expected near 0): '
      f'{np.std(means_distances)}')
print(f'Std.dev. of euclidean distance: {np.mean(std_devs_distances)}'
      f' expected: {np.sqrt(d**2 / 12)}')
print(f'Std.dev. of std.dev across experiments (expected near 0): '
      f'{np.std(std_devs_distances)}')
print(f'Maximum distance encountered: {np.max(max_distances)}')

It outputs:

Mean of euclidean distance: 0.25000092324558143 expected: 0.25
Std.dev. of mean across experiments (expected near 0): 0.0004582988807651125
Std.dev. of euclidean distance: 0.1443138166035262 expected: 0.14433756729740643
Std.dev. of std.dev across experiments (expected near 0): 0.00028820418816544
Maximum distance encountered: 0.4999995726161469

A visualization for $n=3$ was also made for many neighbors (black) of $S$(white):

# https://stackoverflow.com/questions/31768031/plotting-points-on-the-surface-of-a-sphere-in-pythons-matplotlib
r = 1
phi, theta = np.mgrid[0:np.pi:101j, 0:2 *np.pi:101j]
x = r*np.sin(phi)*np.cos(theta)
y = r*np.sin(phi)*np.sin(theta)
z = r*np.cos(phi)

mlab.figure(1, bgcolor=(1, 1, 1), fgcolor=(0, 0, 0), size=(400, 300))
mlab.clf()
mlab.mesh(x , y , z, color=(0.0,0.5,0.5), mode='axes')

# Generate S and some S'
S = random_point_in_sphere(3)
d = 0.3
neighbors = np.array([generate_neighbor(S, d) for _ in range(10000)])
# Add S and its neighbors to plot
xx, yy, zz = neighbors[:, 0], neighbors[:,1], neighbors[:,2]
mlab.points3d(xx, yy, zz, scale_factor=0.02)
mlab.points3d(S[0], S[1], S[2], scale_factor=0.05, color=(1, 1, 1))

# Show x,y,z axes (there should be a better way to do this)
mlab.points3d(np.linspace(1, 2, 100, dtype=float), [0]*100, [0]*100,
              scale_factor=0.01, color=(1, 0, 0))
mlab.points3d([0]*100, np.linspace(1, 2, 100, dtype=float), [0]*100,
              scale_factor=0.01, color=(0, 1, 0))
mlab.points3d([0]*100, [0]*100, np.linspace(1, 2, 100, dtype=float),
              scale_factor=0.01, color=(0, 0, 1))

mlab.show()