I would like to uniformly sample points from the d-spherical cap. The problem is as follow:

Given a point X generated on S^D and the maximum angular seperation theta, generate a point Y such as angle(X,Y) <= theta. The point Y should be generate uniformly from the d-spherical cap.

To generate a random point on S^D we can use the following code:

def get_point(N, D):
    X = np.random.randn(N, D + 1)
    X /= np.linalg.norm(X, axis=1)[:, np.newaxis]
    return X

I tried to apply the rotation matrix to a given input point X

def get_sample_around_point(point, theta_max):
    D = point.shape[0] - 1
    theta = np.random.uniform(0, theta_max)
 
    rotation_axis = np.random.normal(size=(D + 1,))
    rotation_axis /= np.linalg.norm(rotation_axis)
    cos_theta = np.cos(theta)
    sin_theta = np.sin(theta)
    rotation_matrix = (
        cos_theta * np.eye(D + 1)
        + sin_theta * np.outer(rotation_axis, rotation_axis)
        + (1 - cos_theta) * np.eye(D + 1)
    )
    
    new_point = np.dot(rotation_matrix, point.T).reshape(1, -1)
    new_point /= np.linalg.norm(new_point, axis=1)[:, np.newaxis]
    return new_point.flatten()

Example:

N = 1
D = 2
theta = np.pi/4
x = get_point(N, D)[0]
get_sample_around_point(x, theta)

However, the points are not sampled uniformly.

N = 1
D = 2
x = get_point(N, D)[0]
theta = np.pi/3
Y = np.array([get_sample_around_point(x, theta) for _ in range(500)])

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x[0], x[1], x[2], label='X', s=500, color='black')
ax.scatter(Y[:, 0], Y[:, 1], Y[:, 2], color='red', label='')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()

To confirm it we can plot the histogram of angular distances from the sampled point x. The distribution is not uniform.