How to generate uniformly random points on n-spheres and n-balls

For many Monte Carlo methods, such as those in graphical computing, it is critical to uniformly sample from $d$-dimensional spheres and balls. This post describes over twenty different methods to uniformly random sample from the (surface of) a $d$-dimensional sphere or the (interior of) a $d$-dimensional ball.

 

Introduction.

This post gives a comprehensive list of the twenty most frequent and useful methods to uniformly sample from a the surface of a  $d$-sphere, and the interior of the $d$-ball.

Many of them are simple and intuitive; others utilize well-known math methods in sophisticated ways; whilst others can only be described as truly remarkable in their structure. The choice of which one suits your needs will most likely be a balance between code complexity and time efficiency. However, other aspects such as variance reduction may also be of considerable importance.

It is of course, always possible to use the most generalised versions of these formulae as outlined in methods 19-22 for all $d$, but I specifically list some of the specialised methods for smaller $d$ because they may be more efficient, relevant to your application,  and/or amenable to understanding than their  more general expressions. Furthermore, many of the methods applicable to the lower dimensions do not have simple or natural generalization.

 

Definition of the $d-$sphere and $d$-ball.

Let us recap the standard definitions of spheres and balls in $d$-dimensions.

A unit $d-$dimensional sphere is defined such that:

$$ S^d = \left\{   x \in \Bbb{R}^{d+1} : \left \lvert x \right \rvert = 1    \right\} $$

And a unit $d-$dimensional ball is defined such that:

$$ B^d = \left\{   x \in \Bbb{R}^{d} : \left \lvert x \right \rvert  \leq 1    \right\} $$

So the perimeter of a circle is a 1-sphere, and the interior of the circle (a disk) is a 2-ball.

And if the earth was perfectly spherical, than the surface of earth would be a 2-sphere, and the interior of  the earth would be a 3-ball.

The way I remember these definitions, is that for both the $d$-sphere and the $d$-ball, the $d$ signifies how many degrees of freedom it has.

 

Uniform Random Sampling

In this post, we adopt the convention that ‘uniform random sampling‘, ‘(continuous) simple random sampling‘, and even ‘uniform sampling‘ are all fully synonymous. And that for discrete probabilities, this means that all possible elements of $S$ have an equal probability of being selected. For continuous probabilities, this means that the likelihood of an element falling in any subinterval is directly proportional to the length of the subinterval.

Virtually all computer languages have a built-in function for this concept. For example , $\text{rand}(\cdot)$, that produces a random number that is (pseudo)-randomly drawn from uniform distribution $[0,1)$.

 

Note 1.

In contrast to many of my other posts on this blog, we are not talking about quasirandom sampling (aka ‘evenly distributed’), we are just talking about ‘plain old simple vanilla’ uniform random sampling.

 

Note 2.

For the coding examples, I generally use readable python code (built on numpy), except for acceptance/rejection methods where I used pseudo-code. The purpose of the code is aimed for clarity and so does not necessarily include language-specific conventions or optimizations. It uses the notation that $x**y$ is equal to $x^y$.

Although each of them have their strengths and weaknesses, an asterisk (*) indicates methods that I believe are an excellent default option for that section.

Uniformly sampling a 1-sphere (a circle).

*Method 1. Polar

This should be very clear to everyone, and is the foundation for many of the other methods.

theta = 2*PI * random()
x = cos(theta)
y = sin(theta)

 

Method 2. Rejection Method

A lesser-known method (Cook, 1957) that does not directly require trigonometric functions.

You might notice that the expressions have almost identical structure to Pythagorean triplets!

(This method is better known when generalised to the 2-sphere.)

Select u,v ~ U(-1,1)
d2 = u^2+v^2
If d2 >1 then
    reject and go back to step 1
Else
    x = (u^2-v^2)/d2
    y = (2*u*v)/d2
Return (x,y)

 

 

Method 3. Muller , Marsaglia (‘Normalised Gaussians’)

This method was first proposed by Muller and then Marsaglia and the popularised by Knuth, it is very elegant. It naturally generalizes and remains efficient for higher dimensions.

Presumably, it is less well-known due to the less obvious (magical?!) relationship between a circle and the Normal Distribution.

Note that in this version (and its generalizations), the standard deviation does not have to be 1, it just has to be identical for all the random variates. However, in nearly all implementations, a standard error of 1 is picked due to its convenience.

Also, in the absence of a native function, a simple and common method to draw from a normal distribution is via the Box-Mueller algorithm.

u = np.random.normal(0,1)
v = np.random.normal(0,1)
d= (u^2+v^2)**0.5  # notation for sqrt(.)
(x,y) = (u,v)/d

To very briefly understand why this works, consider the function

$$ f(u,v) = f(\mathbf{z}) = \left( \frac{1}{\sqrt{2\pi}} e^{-\frac{1}{2}u^2} \right) \left( \frac{1}{\sqrt{2\pi}} e^{-\frac{1}{2}v^2} \right) $$

After some algebra, this simplifies to

$$ f(u,v) = f(\mathbf{z}) = \left( \frac{1}{(\sqrt{2\pi})^{3/2} }\right) e^{-\frac{1}{2}(u^2+v^2)} =  \left( \frac{1}{(\sqrt{2\pi})^{3/2} }\right) e^{-\frac{1}{2}(|\mathbf{z}|^2)}  $$

This illustrates that the probability distribution of $\mathbf{z}$ only depends on its magnitude and not any direction. Thus it must be rotationally invariant (symmetric). And thus, must correspond to a uniform distribution on the circle.

This argument is not specific to $d=2$ and can be generalized to any $d$.

#

 

Uniformly sampling a 2-ball (disk).

 

*Method 4. Rejection Method

This is probably the most intuitive method, and is quite fast for a 2-ball.

However, (as discussed later) it generally gets a bad rap, because if you generalize this algorithm to higher dimensions, it can become catastrophically very inefficient.

The rejection method for $B^2$ has a ~78% acceptance rate.
1. Select x,y ~ U(-1,1)
2. If (x^2+y^2 >1) then reject and go to step 1.
3. Return (x,y)

 

Method 5. Polar + Radial CDF

For balls of any dimension, a very common method of uniformly picking from a $d$-ball, is to first select a random directional unit-vector from a $(d-1)$-sphere, and then multiply this vector by a scalar radial multiplicative factor.

The crucial aspect of this algorithm is the multiplicative factor, $\sqrt{\cdot}$ function. This is intuitively required because area, which corresponds to the cumulative distribution function (CDF), grows according to $r^2$, so we need to apply the inverse of this. See here for a more visual explanation of this.

The following method is notably asymmetrical. Geometrically, one variable is used for radial distance and the other for the angular displacement.

Figure 5. Polar mapping. Horizontal lines map to concentric circles.
u = random()
v = random()
r = u**0.5  # sqrt function
theta = 2* PI *v
x = r*cos(theta)
y = r*sin(theta)

*Method 6. Concentric Map (Shirley 1997)

This variance reduction method  (Shirley, 1997) is fascinating and motivated by the idea that ideal mappings from the square-based random variate $uv$-space to the coordinate $xy$-space (disks) should be area-preserving, two-way continuous, and with low distortion. By comparing the figures for method 5 and method 6, it should be clear why this one generally produces better results in (finite) MC methods. Additional visual explanation and analysis of this method is given here (which also include some other variations including ‘squircle’, elliptical and approximate equal areas).

 

Figure 6. Concentric Mapping maps concentric squares to concentric circles.
u= random()
v = random()
if (u==0 && v==0) return (0,0)
theta=0
r = 1
a=2*u-1
b=2*v-1
if(a*a>b*b)
    {r = a; phi = PI/4*b/a}
else
    {r = b; PI/2- PI/4*a/b}
x= r*cos(theta)
y= r*sin(theta)

 

Method 7. Muller / Marsaglia (‘Normalized Gaussians’)

This method is a hybrid of methods 3 and 5, because it selects a random direction via method 3 and a random radius via method 5.

u = np.random.normal(0,1)
v = np.random.normal(0,1)
norm = (u*u+v*v)**0.5
r = pow( np.random(0,1),1/2.)
(x,y) = r*(u,v)/norm

Method 8. Exponential Distribution

Another stunningly elegant method (Barthe, 2005) that is valid for balls of any dimension, is to draw $d$ random variates from the normal distribution, and a single random variable from the exponential distribution (with $\lambda=1/2$).

(Note that I believe there is a typo in the original paper, as the original paper says $\lambda = 1$, which has then also been copied into Wolfram’s mathworld. )

Recall, that in the absence of a native function, a simple and common method to draw from an exponential distribution with parameter $\lambda$, is via the negative logarithm of a uniform random variate divided by $\lambda$.

u = np.random.normal(0,1) 
v = np.random.normal(0,1)
e = np.random.exponential(0.5)
denom = (u*u + v*v + e)**0.5
(x,y) = (u,v)/denom

 

Method 9. Dropped Coordinates

A somewhat magical method, first noted by (Harman and Lacko 2010)  and then proven ( Voelker 2017) that is not well-known is that if $(x_1,x_2,x_3,x_4)$ is a random vector uniformly distributed on the $3$-sphere, the random vector $(x_1,x_2)$, (that is, with the last two coordinates dropped) is uniformly distributed in the $2$-ball.

See method 20 for a more detailed commentary of this method.

s = np.random.normal(0,1)
t = np.random.normal(0,1)
u = np.random.normal(0,1)
v = np.random.normal(0,1)
norm = (s*s + t*t + u*u + v*v)**(0.5)
(x,y) = (u,v)/norm
Results for n=3000, for methods 8 (radial CDF), 9 (exponential) and 10 (dropped coordinates). It is truly remarkable that such different techniques can produce the same result!

 

 

Uniformly sampling a 2-sphere
(surface of a sphere).

 

Method 10. Polar

This is the natural generalization of method 1. Note that the inverse cosine is required for the same reason that the sqrt is required in method 5.

u = random()
v = random()
theta = 2 * PI * u
phi = np.arccos(2v-1)
x   = np.sin(theta) np.cos(phi)
y   = np.sin(theta) np.sin(phi)
z   = np.cos (theta)

Using $u=\cos \theta$, an equivalent version of this that does not directly involve inverse trig functions is

u   = 2*random()-1
phi = 2* PI * random()
x = cos(phi) * (1-z^2)**0.5
y = sin(phi)* (1-z^2)**0.5
z = u

*Method 11. Muller

This is a direct generalization of method 3.

u = np.random.normal(0,1)
v = np.random.normal(0,1)
w = np.random.normal(0,1)
norm = (u*u + v*v + w*w)**(0.5)
(x,y,z) = (u,v,w)/norm

 

Method 12. Rejection Method #1

This method (Marsaglia, 1972) is a direct generalization of Method 2 .

Select u,v ~ U(-1,1)
d2 = u*u + v*v
If d2>1 then
    reject and go to step 1.
Else
    x = 2*u* sqrt (1-d2)
    y = 2*v* sqrt (1-d2)
    z = 1-2*d2
Return (x,y,z)

Method 13. Rejection Method #2

This method (Cook, 1957) is quite different to any of the previously discussed methods.

For those variables that are accepted, the rules for quaternion transformations, imply that the following points are uniformly distributed.

Select s,t,u,v ~ U(-1,1)
d2 = s^2 + t^2 + u^2 + v^2
If d2>1 then
    reject and go to step 1.
Else
    x = 2*(s*u+t*v)/d2
    y = 2*(u*v-s*t)/d2
    z = (s*s-t*t+u*u-v*v)/d2
Return (x,y,z)

 

Uniformly sampling a 3-ball
(interior of a sphere).

*Method 14. Rejection Method

This method is a direct generalization of Method 2. Again, it is quite efficient as the acceptance rate is still relatively high and the calculations can be done blazingly fast. However, it carries a bad reputation as it is catastrophically inefficient for higher dimensions.

1. Select x,y,z ~ U(-1,1)
2. If (x^2+y^2 +z^2>1) then reject and go to step 1.
3. Return (x,y,z)

 

The rejection method for $B^3$ has a ~52% acceptance rate.

*Method 15. Muller

This is combination of method 7.

u = np.random.normal(0,1)
v = np.random.normal(0,1)
w = np.random.normal(0,1)
r = random()**(1./3)
norm= (u*u+v*v+w*w)**(0.5)
(x,y,z) = r*(u,v,w)/norm

 

*Method 16. Polar

This is the natural generalization of method 5.

u = 2*random()-1
phi = 2* PI * random()
r = random()**(1/3.)
x = r* cos(phi) * (1^2-z^2)**0.5
y = r* sin(phi)* (1^2-z^2)**0.5
z = r*u

Method 17. Exponential Distribution

This is a direct generalisation of method 8. Again, note that I believe there is a typo in the original paper, as the original paper says $\lambda = 1$, which has then also been copied into Wolfram’s mathworld.

u = np.random.normal(0,1)
v = np.random.normal(0,1)
w = np.random.normal(0,1)
e = np.random.exponential(0.5)
denom = (e + u*u + v*v + w*w)**0.5
(x,y,z) = (u,v,w)/denom

 

Method 18. Dropped Coordinates

A natural extension of method X.

Based on the fact that if $(x_1,x_2,x_3,x_4,x_5)$ is a random vector uniformly distributed on the $4$-sphere, the random vector $(x_1,x_2,x_3)$, (that is, with the last two coordinates dropped) is uniformly distributed in the $3$-ball.

See method 20 for a more detailed commentary of this method.

s = np.random.normal(0,1)
t = np.random.normal(0,1)
u = np.random.normal(0,1)
v = np.random.normal(0,1)
w = np.random.normal(0,1)
norm = (s*s + t*t + u*u + v*v + w*w)**(0.5)
(x,y,z) = (u,v,w)/norm

 

Uniformly sampling a $d$-sphere.

*Method 19. Muller

For higher dimensions where $d>3$, most of the above methods do not neatly simplify. Therefore, the one that is almost always used for a $d$-sphere is the generalised version of methods 3 and 9.

u = np.random.normal(0,1,d)  # an array of d normally distributed random variables
d=np.sum(u**2) **(0.5)
(x,y,z) = (u,v,w)/d

 

Uniformly sampling a $d$-ball.

Although the rejection methods 4 and 13 are extremely intuitive and trivially generalise to higher dimensions, one quickly encounters ‘the curse of dimensionality’. As mentioned, for $d=1$ the acceptance rate is ~71%, and for $d=2$ the acceptance rate is 52%. Although this might appear inefficient, due to not requiring to calculate any complex functions this is often one of the fastest ways to pick points in a $2-$ball or $3$-ball.

By comparing the hypervolume of a $d$-ball with the enclosing cube, one can deduce that for $d=10$, the acceptance rate is less than 0.25% and it gets exponentially worse for higher $d$.

*Method 20. Muller

In practice, the almost exclusive method to calculate points inside a $d-$ball is the generalized version of method 14. Namely,

u = np.random.normal(0,1,d)  # an array of d normally distributed random variables
norm=np.sum(u**2) **(0.5)
r = random()**(1.0/d)
x= r*u/norm

 

Method 21. Exponential Distribution

This is a direct generalisation of method 8 and 16.

x = np.random.normal(0,1,d)
e = np.random.exponential(0.5)
denom = (e + np.sum(u**2))**0.5 
(x,y,z) = (u,v,w)/denom

 

*Method 22. Dropped Coordinates

A very elegant method, first noted by (Harman and Lacko 2010)  and then proven ( Voelker 2017) that is not well-known is that if $ (x_1,x_2,x_3,…,x_d,x_{d+1},x_{d+2})$ is a random vector uniformly distributed on the $(d+1)$-sphere, the random vector $ (x_1,x_2,x_3,…,x_d)$, (that is, with the last two coordinates dropped) is uniformly distributed in the $d$-ball.

This is non-trivial result where even developing an intuition is difficult, (even the base case of $d=1$ is not trivial!)

However, it may be interesting/useful to note that an $\exp(1/2)$ variable $Y$ in method 19, has the same distribution as $(Y_1^2+Y_2^2)$  where $Y_1,Y_2$ are standard normal – thus connecting this methods 19 and 20.

Thus,

u = np.random.normal(0,1,d+2)  # an array of (d+2) normally distributed random variables
norm=np.sum(u**2) **(0.5)
u = u/norm
x = u[0:d] #take the first d coordinates

 

Conclusion

This post gives a comprehensive list of the most frequent and useful methods to uniformly sample from a the surface of a  $d$-sphere, an the interior of the $d$-ball.

Many of them are simple and intuitive, whilst others are truly remarkable in their structure.

It is of course, possible to always use the most generalised versions of these formulae as outlined in methods 19-22, but sometimes specialised methods for smaller $d$ may be more efficient and/or amenable to understanding.

 

 

***

My name is Martin Roberts. I have a PhD  in theoretical physics. I love maths and computing. I’m open to new opportunities – consulting, contract or full-time – so let’s have a chat on how we can work together!

Come follow me on Twitter: @Techsparx!
My other contact details can be found here. 

Other Posts you may like

Leave a Reply

Your email address will not be published. Required fields are marked *