Efficiently sampling points uniformly from the surface of an n-sphere

Question

Is there an efficient way to sample uniformly points from the unit n-sphere? Informally, by "uniformly" I mean the probability of picking a point from a region is proportional to the area of that region on the surface of the sphere. Formally, I guess I'm referring to the Haar measure.

I guess "efficient" means the algorithm should take poly(n) time. Of course, it's not clear what I mean by an algorithm since real numbers cannot be represented on a computer to arbitrary precision, so instead we can imagine a model where real numbers can be stored, and arithmetic can be performed on them in constant time. Also, we're given access to a random number generator which outputs a real in [0,1]. In such a model, it's easy to sample from the surface of the n-hypercube in O(n) time, for example.

If you prefer to stick with the standard model of computation, you can consider the approximate version of the problem where you have to sample from a discrete set of vectors that $\epsilon$-approximate the surface of the n-sphere.

Answer 1

Generate $X_1, X_2, \ldots, X_n$ independent, normally distributed random variables See wikipedia for information on how to do this given some standard source of randomness, for example uniform(0,1) random variables.

Then let $Y_i = X_i/\sqrt{X_1^2 + \cdots + X_n^2}$ for $i = 1, \ldots, n$. Then $(Y_1, \ldots, Y_n)$ is uniformly distributed on the surface of the sphere. The time this takes is linear in $n$.

This works because the multivariate normal $(X_1, \ldots, X_n)$ with covariance matrix the identity (that is, $n$ independent unit normals) is rotationally symmetric around the origin.

Answer 2

You can just take $n$ independent Gaussians and divide the resulting vector by its absolute value. If you have a random number generator that provides you with the standard Gaussian, you can just use it. If all you can generate is a random variable uniformly distributed $[0,1]$, just take two copies $X$ and $Y$ and let $Z=\cos(\pi X)\sqrt{-\log Y}$ (this is a bit simpler than to try to invert erf directly). $Z$ will be a Gaussian with mean $0$ and some variance that doesn't really matter because you normalize in the end anyway.

Answer 3

There is a completely different way of doing this based on Hurwitz's 1897 "factorization of measures". Essentially, you make sequential rotations in $(12)$-plane, then in $(13)$-plane,..., and finally $(n-1,n)$-plane. You can calculate explicitly the $(ij)$-distribution in terms of beta functions. The result is a random rotation in $SO(n)$ which in particular gives you a uniform point in the sphere. See section 2 in this paper by Diaconis and Saloff-Coste for explicit formulas and references.

Answer 4

In case you don't want to worry about dividing by small norms (admittedly, less and less of an issue for higher $n$...)

Generate a point $\mathbf{x}$ on the $(n-1)$-sphere; and generate a number $y\in [-1,1]$ with density $k (1-y^2)^{(n-2)/2}$ for the appropriate constant $k$. Letting $\mathbf{x}'=(y,(1-y^2)^{1/2} \mathbf{x})$ is a uniform point on the $n$-sphere.

---

Alternatively, in the unlikely event that you don't mind funny correlations between consecutively generated points, you might prefer randomly generating elements of the orthogonal group, for instance letting $A_n$ be sufficiently random${}^1$ antisymmetric matrices, $O_n=(I-A_n)^{-1}(I+A_n)$ is orthonormal, and then the sequence of matrices $U_0=I$, $U_{n+1} = O_n U_n$ is uniformly distributed w.r.t. Haar measure. If later you want to kill consecutive correlations, generate a bunch of them and then sample the sequence at random.

---

${}^1$: "Sufficiently random": Leveraging the mean-value theorem for all it's worth, we might even take the $A_n$s concentrated on a "generic curve in general position". For sure this isn't efficient, but it gives an idea of how little is needed to get quantitative convergence to uniformity.