I recently ran across a simple way to generate random points uniformly distributed on a sphere: Generate random points in a cube until you get one inside the sphere, then push it to the surface of the sphere by normalizing it [1].
In more detail, to generate random points on the unit sphere, generate triples
(u1, u2, u3)
of independent random values uniformly generated on [−1, 1] and compute the sum of there squares
S² = u1² + u2² + u3²
If S² > 1 then reject the sample and try again. Otherwise return
(u1/S, u2/S, u3/S).
There are a couple things I like about this approach. First, it’s intuitively plausible that it works. Second, it has minimal dependencies. You could code it up quickly in a new environment, say in an embedded system, though you do need a square root function.
Comparison with other methods
It’s not the most common approach. The most common approach to generate points on a sphere in n dimensions is to generate n independent values from a standard normal distribution, then normalize. The advantage of this approach is that it generalizes efficiently to any number of dimensions.
It’s not obvious that the common approach works, unless you’ve studied a far amount of probability, and it raises the question of how to generate samples from normal random variable. The usual method, the Box-Muller algorithm, requires a logarithm and a sine in addition to a square root.
The method starting with points in a cube is more efficient (when n = 3, though not for large n) than the standard approach. About half the samples that you generate from the cube will be thrown out, so on average you’ll need to generate six values from [−1, 1] to generate one point on the sphere. It’s not the most efficient approach—Marsaglia gives a faster algorithm in the paper cited aove—but it’s good enough.
Accept-reject methods
One disadvantage to accept-reject methods in general is that although they often have good average efficiency, their worst-case efficiency is theoretically infinitely bad: it’s conceivable that you’ll reject every sample, never getting one inside the sphere. This is not a problem in practice. However, it may be a problem in practice that the runtime is variable.
When computing in parallel, a task may have to wait on the last thread to finish. The runtime for the task is the runtime of the slowest thread. In that case you may want to use an algorithm that is less efficient on average but that has constant runtime.
It also matters in cryptography whether processes have a constant runtime. An attacker may be able to infer something about the path taken through code by observing how long it took for the code to execute.
Related posts
- Evenly distributing points on a sphere
- Evolution of random number generators
- Testing random number generators
[1] George Marsaglia. Choosing a point from the surface of a sphere. The Annals of Mathematical Statistics. 1972, Vol. 43, No. 2, 645–646.
When you generate a point (x,y,z) that happens to land outside the sphere, consider the point (x-sign(x), y-sign(y), z-sign(z)). This is a point chosen uniformly at random from a spiky shape obtained by rearranging the complement of the unit sphere. The spiky shape contains the cube [-t,t]^3 where t = sqrt(3)-1. So you can get a second bite at the cherry without generating any more random numbers: first check whether (x,y,z) is in the unit sphere; if not, let x’ = (x-sign(x))/(sqrt(3)-1) etc., and check whether (x’,y’,z’) is in the unit sphere.
In principle you could iterate this, but every time you do it you’re magnifying things and therefore losing a few bits of precision. But doing it just once takes you from a fraction pi/6 ~= 52% of the cube to a fraction (1+t^3)pi/6 ~= 73% of the cube; you go from needing about 5.7 random numbers per sphere-surface point to needing about 4.1 random numbers per sphere-surface point.