Creating random points in the surface of a n-dimensional sphere

I have a point X in the surface of an n-dimensional sphere with center 0.

I want to create random points following a distribution with center X, the points must be in the surface of the n-dimensional sphere, and located very close to X.

With spherical coordinates in 3D, I can put some random noise in the two angles defining the point X.

In the general case, I can create random points with a Gaussian in n-dimensions with center X, and project them into the surface of the sphere, by making the euclidean norm of the random points to be equal to the radius of the sphere (this works because the center of the sphere is 0).

Do you have any better ideas about efficiently creating points like these?

– whuber
Apr 26, 2020 at 12:45
• Why do you object to orthogonal projections? They would seem not only to comply with all your other requirements but also to be extremely efficient.
– whuber
Apr 26, 2020 at 13:45
• Thank you whuber, I just realize making the Euclidean norm of the random points equal to the radius of the sphere is equivalent to orthogonal projections of each point onto the surface. Apr 26, 2020 at 14:10
• Maybe see stats.stackexchange.com/questions/7977/…. The method given there can be generalized to the $n$-sphere Apr 26, 2020 at 14:34
• @Dave Ultimately what we're talking about is how to project a plane to the sphere in such a way that it is nearly an isometry near $x_0.$ It's fruitless to try to re-invent the last 2500 years of geography through ad hoc proposals--nobody has ever been or will be that smart. (There's a famous quotation of Newton about standing on the shoulders of giants...). Use the science to design an appropriate projection. For more about this, see the comments below my answer here.
– whuber
Apr 26, 2020 at 17:22

Using a stereographic projection is attractive.

The stereographic projection relative to a point $$x_0\in S^{n}\subset \mathbb{R}^{n+1}$$ maps any point $$x$$ not diametrically opposite to $$x_0$$ (that is, $$x\ne -x_0$$) onto the point $$y(x;x_0)$$ found by moving directly away from $$-x_0$$ until encountering the tangent plane of $$S^n$$ at $$x_0.$$ Write $$t$$ for the multiple of this direction vector $$x-(-x_0) = x+x_0,$$ so that

$$y = y(x;x_0)= x + t(x+x_0).$$

Points $$y$$ on the tangent plane are those for which $$y,$$ relative to $$x_0,$$ are perpendicular to the Normal direction at $$x_0$$ (which is $$x_0$$ itself). In terms of the Euclidean inner product $$\langle\ \rangle$$ this means

$$0 = \langle y - x_0, x_0 \rangle = \langle x + t(x+x_0) - x_0, x_0\rangle = t\langle x + x_0, x_0\rangle + \langle x-x_0, x_0\rangle.$$

This linear equation in $$t$$ has the unique solution

$$t = -\frac{\langle x-x_0,x_0\rangle}{\langle x + x_0, x_0\rangle}.$$

With a little analysis you can verify that $$|y-x_0|$$ agrees with $$x-x_0$$ to first order in $$x-x_0,$$ indicating that when $$x$$ is close to $$x_0,$$ Stereographic projection doesn't appreciably affect Euclidean distances: that is, up to first order, Stereographic projection is an approximate isometry near $$x_0.$$

Consequently, if we generate points $$y$$ on the tangent plane $$T_{x_0}S^n$$ near its origin at $$x_0$$ and view them as stereographic projections of corresponding points $$x$$ on $$S_n,$$ then the distribution of the points on the sphere will approximate the distribution of the points on the plane.

This leaves us with two subproblems to solve:

1. Generate Normally-distributed points near $$x_0$$ on $$T_{x_0}S^n.$$

2. Invert the Stereographic projection (based at $$x_0$$).

To solve (1), apply the Gram-Schmidt process to the vectors $$x_0, e_1, e_2, \ldots, e_{n+1}$$ where the $$e_i$$ are any basis for $$\mathbb{R}^n+1.$$ The result after $$n+1$$ steps will be an orthonormal sequence of vectors that includes a single zero vector. After removing that zero vector we will obtain an orthonormal basis $$u_0 = x_0, u_1, u_2, \ldots, u_{n}.$$

Generate a random point (according to any distribution whatsoever) on $$T_{x_0}S^n$$ by generating a random vector $$Z = (z_1,z_2,\ldots, z_n) \in \mathbb{R}^n$$ and setting

$$y = x_0 + z_1 u_1 + z_2 u_2 + \cdots + z_n u_n.\tag{1}$$

Because the $$u_i$$ are all orthogonal to $$x_0$$ (by construction), $$y-x_0$$ is obviously orthogonal to $$x_0.$$ That proves all such $$y$$ lie on $$T_{x_0}S^n.$$ When the $$z_i$$ are generated with a Normal distribution, $$y$$ follows a Normal distribution because it is a linear combination of Normal variates. Thus, this method satisfies all the requirements of the question.

To solve (2), find $$x\in S^n$$ on the line segment between $$-x_0$$ and $$y.$$ All such points can be expressed in terms of a unique real number $$0 \lt s \le 1$$ in the form

$$x = (1-s)(-x_0) + s y = s(x_0+y) - x_0.$$

Applying the equation of the sphere $$|x|^2=1$$ gives a quadratic equation for $$s$$

$$1 = |x_0+y|^2\,s^2 - 2\langle x_0,x_0+y\rangle\, s + 1$$

with unique nonzero solution

$$s = \frac{2\langle x_0, x_0+y\rangle}{|x_0+y|^2},$$

whence

$$x = s(x_0+y) - x_0 = \frac{2\langle x_0, x_0+y\rangle}{|x_0+y|^2}\,(x_0+y) - x_0.\tag{2}$$

Formulas $$(1)$$ and $$(2)$$ give an effective and efficient algorithm to generate the points $$x$$ on the sphere near $$x_0$$ with an approximate Normal distribution (or, indeed, to approximate any distribution of points close to $$x_0$$).

Here is a scatterplot matrix of a set of 4,000 such points generated near $$x_0 = (1,1,1)/\sqrt{3}.$$ The standard deviation in the tangent plane is $$1/\sqrt{12} \approx 0.29.$$ This is large in the sense that the points are scattered across a sizable portion of the $$x_0$$ hemisphere, thereby making this a fairly severe test of the algorithm.

It was created with the following R implementation. At the end, this R code plots histograms of the squared distances of the $$y$$ points and the $$z$$ points to the basepoint $$x_0.$$ By construction, the former follows a $$\chi^2(n)$$ distribution. The sphere's curvature contracts the distances the most when they are large, but when $$\sigma$$ is not too large, the contraction is virtually unnoticeable.

#
# Extend any vector x0 to an orthonormal basis.
# The first column of the output will be parallel to x0.
#
gram.schmidt <- function(x0) {
n <- length(x0)
V <- diag(rep(1, n))                 # The usual basis of R^n
if (max(x0) != 0) {
i <- which.max(abs(x0))            # Replace the nearest element with x0
V <- cbind(x0, V[, -i])
}
L <- chol(crossprod(V[, 1:n]))
t(backsolve(L, t(V), transpose=TRUE))
}
#
# Inverse stereographic projection of y relative to the basepoint x0.
# The default for x0 is the usual: (0,0, ..., 0,1).
# Returns a point x on the sphere.
#
iStereographic <- function(y, x0) {
if (missing(x0) || max(abs(x0)) == 0)
x0 = c(1, rep(0, length(y)-1)) else x0 <- x0 / sqrt(sum(x0^2))

if (any(is.infinite(y))) {
-x0
} else {
x0.y <- x0 + y
s <- 2 * sum(x0 * x0.y) / sum(x0.y^2)
x <- s * x0.y - x0
x / sqrt(sum(x^2))                    # (Guarantees output lies on the sphere)
}
}
#------------------------------------------------------------------------------#
library(mvtnorm)                        # Loads rmvnorm
n <- 4e3
x0 <- rep(1, 3)
U <- gram.schmidt(x0)
sigma <- 0.5 / sqrt(length(x0))
#
# Generate the points.
#
Y <- U[, -1] %*% t(sigma * rmvnorm(n, mean=rep(0, ncol(U)-1))) + U[, 1]
colnames(Y) <- paste("Y", 1:ncol(Y), sep=".")

X <- t(apply(Y, 2, iStereographic, x0=x0))
colnames(X) <- paste("X", 1:ncol(X), sep=".")
#
# Plot the points.
#
if(length(x0) <= 8 && n <= 5e3) pairs(X, asp=1, pch=19, , cex=1/2, col="#00000040")
#
# Check the distances.
#
par(mfrow=c(1,2))
y2 <- colSums((Y-U[,1])^2)
hist(y2, freq=FALSE, breaks=30)
curve(dchisq(x / sigma^2, length(x0)-1) / sigma^2, add=TRUE, col="Tan", lwd=2, n=1001)

x0 <- x0 / sqrt(sum(x0^2))
z2 <- colSums((t(X) - x0)^2)
hist(z2, freq=FALSE, breaks=30)
curve(dchisq(x / sigma^2, length(x0)-1) / sigma^2, add=TRUE, col="SkyBlue", lwd=2, n=1001)
par(mfrow=c(1,1))

• Cool! :) Would it be worth mentioning in passing why one might prefer one projection over another? Apr 26, 2020 at 16:43
• @Alexis There are myriad projections (for examples, search Geographic Information Systems). A good one in this circumstance will (a) have little metric distortion in neighborhoods of $x_0$ and $(b)$ be conformal (so that all the distortion it introduces is radial, thereby treating all coordinates on an equal footing). These characteristics, along with the simple formulas involved, lead to the Stereographic projection. For many purposes the Riemannian exponential function (aka "Equidistant projection") would work even better, but is more complicated to implement and isn't conformal except right at $x_0.$
– whuber
Apr 26, 2020 at 16:47
• Instead of Gram-Schmidt, it might be useful to use Stewart's procedure instead. Apr 27, 2020 at 0:16
• @J.M.isnotastatistician Thank you for the suggestion. As the code here shows, I use a Cholesky decomposition. I'm not sure how you propose using Stewart's procedure, which--according to a quick read of the abstract--generates random orthonormal bases. It may be important to have control over this generation (which is achieved in my proposal by means of the original orthonormal basis), so a random procedure will be useful only for isotropic multinormal distributions.
– whuber
Apr 27, 2020 at 12:59
• Is the support of this approximate normal distribution the space of angles $\vec{\theta} \in \mathbb{R}^{n-1}$? Feb 25 at 17:25

• I want to create random points following a distribution with center X, the points must be in the surface of the n-dimensional sphere, and located very close to X.

This does not specify the problem in much detail. I will assume that the distribution of the points is spherically symmetric around the point X and that you have some desired distribution for the (Euclidian) distance between the points and X.

You can consider the n-sphere sphere as a sum of (n-1)-spheres, slices/rings/frustrums.

Now we project a point from the n-sphere, onto the n-cylinder around it. Below is a view of the idea in 3 dimensions.

https://en.wikipedia.org/wiki/File:Cylindrical_Projection_basics.svg

The trick is then to sample the height on the cylinder and the direction away from the axis separately.

Without loss of generality we can use the coordinate $$(1,0,0,0,...,0)$$ (solve it for this case and then rotate the solution to your point $$X$$).

Then use the following algorithm:

• Sample the coordinate $$x_1$$ by sampling which slices the points end up in according to some desired distance function.
• Sample the coordinates $$x_2, ..., x_n$$ by determining where the points end up on the (n-1)-spheres (this is like sampling on a (n-1)-dimensional sphere with the regular technique).

Then rotate the solution to the point $$X$$. The rotations should bring the first coordinate $$(1,0,0,0, ..., 0)$$ to the vector $$X$$, the other coordinates should transform to vectors perpendicular to $$X$$, any orthonormal basis for the perpendicular space will do.

• How do you achieve a close approximation to a Normal distribution?
– whuber
Apr 27, 2020 at 13:05
• @whuber I did not write down explicitly how to achieve this, since this is not clear to me what 'a close approximation to a Normal distribution' is supposed to mean. Also, I do not recognise this in the OP's question. I just gave a general formalism. For any distribution function of the distance (which is obviously limited by the half the circumference), you can use this algorithm. One only has to sample from this distance-distribution and after that compute the coordinate $x_1$ (or the height on the cylinder) based on that distance. Apr 27, 2020 at 13:37
• The problem is that "compute the coordinate" can introduce extreme distortion. As far as "close approximation" goes, that seems intuitively clear, although it is open to various different mathematical interpretations: but at a minimum, you should offer one reasonable interpretation and show why your proposal satisfies that criterion.
– whuber
Apr 27, 2020 at 13:40
• @whuber where does this criterion come from? I have just given an algorithm to compute a distribution for any given distance distribution. Also, not that it is not an approximation, it is an exact solution. Apr 27, 2020 at 13:42
• @whuber the OP's question is very thin on these details. I have just provided a simple method for the case when we can describe the distribution as an isotropic/uniform distribution with the only density variation based on distance from the point. Yes, if anisotropy is involved (but only a single point $X$ is given so this does hijack the OP's question) then the method for rotation needs to be adjusted and requires additional axes/points to determine the rotation (this is the same for any method). If such details are implicitly considered, then the question becomes too broad and unclear. Apr 27, 2020 at 15:48

First, it is not possible to have the positions be exactly Gaussian since restriction to the surface of a sphere imposes a bound on the range of the coordinates.

You could look at using truncated, to $$(-\pi, \pi)$$, normals for each component. To be clear, for a 2-sphere (in 3-space) you have fixed the radius, and must choose 2 angles. I am suggesting you put truncated normal distributions on the angles.

• Thank you. I just edited the question to make it more clear. About your first point, Gaussian points have component standard deviations of 0.01 the radius of the sphere so every random point is really near X and quite near its projection in the surface. The point of the angles I still dont get it, in 3d i can create random points by modifying randomly the 2 angles of polar coordinates, but in n-D i find more convenient to create random points with Gaussian distributions and then project onto sphere surface. Apr 26, 2020 at 10:00
• It depends whether your goal is to get 'some distribution on the sphere' or for that distribution to be some particular (truncated) gaussian. Your projection idea would produce a distribution, but this is unlikely to be gaussian. Consider the circle. If you stand at the origin and draw a standard bivariate normal centered on the origin, the distribution of the projected points should be uniform on the circle. Apr 26, 2020 at 17:42
• That example doesn't sound relevant to the question, because the origin is not an example of a point $X$ "in the surface of the ... sphere" around which the distribution should be focused.
– whuber
Apr 26, 2020 at 19:48
• I was responding to, "i find more convenient to create random points with Gaussian distributions and then project onto sphere surface." Apr 27, 2020 at 6:51

Your center $$X$$ is the projection (normalization) of a vector which has not a normalized norm. You could define a vector $$x$$ which would serve as your unnormalized center and then select data points around $$x$$ (typically using a Gaussian distribution).
There is a free parameter: the norm of $$x$$. As a matter of fact what will matter is the ratio between the standard deviation of $$X$$ and that norm. You will get a value similar to the $$\kappa$$ value of the multi-dimensional Von Mises destribution.