So when I was in the shower earlier and I wondered to myself “Hm. How would I generate random points on an N-sphere?” (true story).

I thought about it for a bit and generated the two obvious solutions. Unfortunately the two obvious solutions are, respectively, wrong and massively inefficient (actually the inefficient one isn’t too bad for low dimensions, but I was interested in the high dimensional case where it rapidly becomes massively inefficient). Once out of the shower I looked the answer up, and found the whole thing interesting enough that I figured someone else might too, so this is a post about that.

Obvious but wrong solution

Generate N random uniform numbers between -1 and 1. Divide the resulting vector by its norm.

Why is this wrong? Well, the initial N random numbers will lie uniformly in the N cube. So the probability density at a given point is proportional to the length of the line that goes from the zero, through that point and out to the N sphere. This line is not of uniform length for the very simple reason that cubes aren’t spheres.

Obvious but slow solution

Generate N uniform random numbers between -1 and 1. Retry this until the resulting norm is <= 1. Divide the resulting vector by its norm. Why does this work? Because we're using a rejection method to make sure we get points inside the unit ball. The result is not uniformly distributed inside the unit ball but it is (evidently) spherically symmetrical. By dividing by the norm we’re then projecting the spherically symmetrical distribution onto the sphere and thus getting a uniform distribution on it.

Why is it slow? Because of the curse of dimensionality. When we do the rejection we’re effectively generating a geometric distribution with probability of success equal to the proportion of the N-cube the N-ball occupies. This means we expect to perform the generation 1 / that probability times, which tends towards infinity quite rapidly as N gets large.

Obvious in retrospect solution

At this point I looked up the answer, though I should have been able to figure it out.

What we want is a spherically symmetric distribution a la the previous solution but which is more efficient to calculate.

It turns out that a good choice for this is an N dimensional normal distribution, which can be obtained by simply generating N Normal(0, 1) random variables. You can then divide by the norm of those variables and get a distribution on the unit sphere.

This is apparently not the best solution (generating normal random numbers is a little expensive), but I couldn’t find a reference to a version which works in N dimensions for any N.

Here’s some C code implementing this: https://gist.github.com/1354700.