# The way to generate uniformly random factors on n-spheres and in n-balls

*by*Phil Tadros

For a lot of Monte Carlo strategies, corresponding to these in graphical computing, it’s essential to uniformly pattern from $d$-dimensional spheres and balls. This put up describes over twenty totally different strategies to uniformly random pattern from the (floor of) a $d$-dimensional sphere or the (inside of) a $d$-dimensional ball.

First printed: 10*th March, 2019* Final up to date: *16 April, 2020*

This put up provides a complete record of the twenty most frequent and helpful strategies to uniformly pattern from a the floor of a $d$-sphere, and the inside of the $d$-ball.

A lot of them are easy and intuitive; others make the most of well-known math strategies in subtle methods; while others can solely be described as actually outstanding of their construction. The selection of which one fits your wants will probably be a stability between code complexity and time effectivity. Nonetheless, different points corresponding to variance discount might also be of appreciable significance.

It’s after all, all the time potential to make use of probably the most generalised variations of those formulae as outlined in strategies 19-22 for all $d$, however I particularly record a few of the specialised strategies for smaller $d$ as a result of they might be extra environment friendly, related to your utility, and/or amenable to understanding than their extra basic expressions. Moreover, most of the strategies relevant to the decrease dimensions should not have easy or pure generalization.

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

Allow us to recap the usual definitions of spheres and balls in $d$-dimensions.

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

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

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

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

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

And if the earth was completely spherical, than the floor of earth can be a 2-sphere, and the inside of the earth can be a 3-ball.

The best way I keep in mind these definitions, is that for each the $d$-sphere and the $d$-ball, the $d$ signifies what number of levels of freedom it has.

### Uniform Random Sampling

On this put up, we undertake the conference that ‘** uniform random sampling**‘, ‘(

*steady) easy random sampling*‘, and even ‘

*uniform sampling*‘ are all totally synonymous. And that for discrete possibilities, because of this all potential parts of $S$ have an equal chance of being chosen. For steady possibilities, because of this the probability of a component falling in any subinterval is instantly proportional to the size of the subinterval.

Just about all laptop languages have a built-in operate for this idea. For instance , $textual content{rand}(cdot)$, that produces a random quantity that’s (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.

That is, uniformly sampling points on the circumferende of 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{2pi}} e^{-frac{1}{2}u^2} right) left( frac{1}{sqrt{2pi}} e^{-frac{1}{2}v^2} right) $$

After some algebra, this simplifies to

$$ f(u,v) = f(mathbf{z}) = left( frac{1}{(sqrt{2pi})^{3/2} }right) e^{-frac{1}{2}(u^2+v^2)} = left( frac{1}{(sqrt{2pi})^{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$.

That is, uniformly sampling points inside a circle.

### *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.

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.

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).

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

That is, uniformly sampling points on the 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)

That is, uniformly sampling points on the inside 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)

### *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-u^2)**0.5 y = r* sin(phi)* (1-u^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

### *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

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 $sim 0.71$, and for $d=2$ the acceptance rate is $sim 0.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.0025 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 primary d coordinates

This put up provides a complete record of probably the most frequent and helpful strategies to uniformly pattern from a the floor of a $d$-sphere, an the inside of the $d$-ball.

A lot of them are easy and intuitive, while others are actually outstanding of their construction.

It’s after all, potential to all the time use probably the most generalised variations of those formulae as outlined in strategies 19-22, however generally specialised strategies for smaller $d$ could also be extra environment friendly and/or amenable to understanding.

***

My title is Dr Martin Roberts, and I’m a contract Principal Information Science guide, who loves working on the intersection of maths and computing.

“I remodel and modernize organizations by means of progressive knowledge methods options.”

You’ll be able to contact me by means of any of those channels.LinkedIn: https://www.linkedin.com/in/martinroberts/

Twitter:

@TechSparxhttps://twitter.com/TechSparxe mail: Martin (at) RobertsAnalytics (dot) com

Extra particulars about me might be discovered

here.