Three Different Algorithms for GeneratingUniformly Distribut.docx

Three Different Algorithms for Generating
Uniformly Distributed Random Points on the
N-Sphere
Jan Poland
Oct 24, 2000
Abstract
We present and compare three different approaches to generate
random points on the N -sphere: A simple Monte Carlo
algorithm, a
coordinate-by-coordinate strategy and a method based on the
rotation
invariance the normal distribution. The latter algorithm is the
fastest.
1 Introduction
Computer scientists sometimes encounter the problem of
generating random
points on the N -dimensional unit sphere SN = {x ∈ RN +1 :
‖x‖2 = 1}.
In most cases, the random points should be distributed
uniformly. For this
task, there exists a nice little algorithm already mentioned in
Knuth ([3]),
that is based on the fact that the N -dimensional normal
distribution is in-
variant under rotation. This algorithm is indubitably the most

efficient, in
particular in high dimensions. But if a special non-uniform
distribution is re-
quested, it may be profitable to know other approaches such as
a coordinate-
by-coordinate strategy or a simple rejection method.
When N is small, each algorithm is sufficiently fast, and the
problem
is not very interesting. For the 1-sphere, i.e. the unit circle in
R2, U =
(cos(α), sin(α)) with α uniformly distributed in [0, 2π] does the
job. Here, we
will concentrate on the high dimensional case. Applications are
for example
in the initialization of a counterpropagation neural network (see
[7]) or a
variant of generalized continuous recombination in an evolution
strategy (cf.
[4]).
Other discussions of our problem can be found in [6] or [5].
2 A Monte-Carlo algorithm
An immediate algorithm is the following. Take a random point
in [−1, 1]N +1
and project it onto the unit sphere . This results in a non-
uniform distri-
bution, therefore one extra step is necessary. Generate a random
point in
X ∈ [−1, 1]N +1 and reject it if it is outside the unit ball, i.e.
‖X‖ > 1. If
not, U = ‖X‖−1 ·X gives a random point on SN (neglecting the

unlikely case
that X = 0), and produces a uniform distribution on the sphere.
This algorithm is often referred to as the rejection method (see
[6]). The
major drawback of the algorithm is it’s expected running time:
It grows
worse than exponentially in N . This is due to the fact that the
volume of
the N -dimesional unit ball
V (N ) =
π
N
2
Γ(1 + N
2
)
more than exponentially decreases as N → ∞, where Γ(·) is the
gamma
function
Γ(x) =
∫ ∞
0
e−t · tx−1 dt.
3 The coordinate approach
In high dimensions, the distribution of each single coordinate of
a uniformly

distributed random point U on the N -sphere becomes quite
complicated.
Suppose we are given a random vector U uniformly distributed
on SN .
Consider the projection U1 of U to its first coordinate, lets say
the x-coordinate.
Then the density function f1(x) of U1 is a bounded function
from [−1, 1] into
R+. We now want to state the f1 explicitely.
The following computations will need the incomplete beta
function
Bx(a, b) =
∫ x
0
ta−1 · (1 − t)b−1 dt
and the beta function B(a, b) = B1(a, b).
To determine f1, take a point u ∈ SN and consider a small move
along
the first coordinate to another point ũ ∈ SN . We observe that
the move
approximately covers the distance
‖ũ − u‖ ≈ | arcsin(x̃) − arcsin(x)|,
and as ũ → u we obtain exact equality, where x and x̃ denote the
first
coordinates of u and ũ, respectively. Hence, we have
∂u

∂x
= arcsin(x)′ =
1√
1 − x2
.
On the other hand, whith a fixed first coordinate x ∈ [−1, 1],
the point u
lays on a (N − 1)-sphere with radius
√
1 − x2 which has the volume
R(x) = V (N − 1) · (
√
1 − x2)N−1.
Since the density of the distribution of U1 has to be linear in
both
∂u
∂x
and
R(x), we obtain
f1(x) = s · arcsin(x)′ · (
√
1 − x2)N−1
= B

(
1
2
,
N
2
)−1
· (
√
1 − x2)N−2
for x ∈ [−1, 1]. The factor s is the appropriate scaling factor
assuring that
the integral over f1 is 1 and thus evaluates to
s =
(∫ 1
−1
(
√
1 − x2)N−2dx
)−1
= B
(
1

2
,
N
2
)−1
.
Now it is straightforward to generate random numbers with
density f1.
We have to find F1 by integrating f1 and finally get F
−1
1 by inverting F1.
Thus, U1 = F
−1
1 (Z) will have the desired properties, if Z is uniformly dis-
tributed on [0, 1]. Performing the integration yields
F1(x) =
∫ x
−1
f1(y)dy =
1
2
+ sign(x) · Bx2
(
1
2

, N
2
)
2 · B
(
1
2
, N
2
).
The last step of inverting F1 cannot be done in a closed form.
So one has
to employ an approximation algorithm such as Newton’s method
to perform
this task numerically.
Once constructed a properly distributed U1, the rest is done by
recursion:
The remaining task is to generate a random point distributed
uniformly on
the (N − 1)-sphere with radius
√
1 − x2. The recursion will finally stop at
N = 0, where it remains to randomly choose a point out of {−1,
1}. For
the reason of efficiency and accuracy, one can also treat the
cases N = 1

and N = 2 seperately: For N = 1 one has to construct one
coordinate of
a point on the unit circle, while for N = 2 the formulas yield f1
≡ 12 and
F1(x) =
1
2
(x + 1).
The time complexity of this algorithm is clearly linear, as
solving the
equation can be done in constant time. The main drawback of a
practical
implementation of this algorithm is the approximation that has
to be done
to solve the equation. Even a rapidly convergent approximation
such as
Newton’s method restrains the performance considerably.
Note:
var(randball(d)) =
√
πΓ( d+1
2
)
2Γ( d+4
2
)B( 1

2
, d+1
2
)
=
Γ(1 + d
2
)
2Γ(2 + d
2
)
4 Using normally distributed random vectors
This is the algorithm mentioned in Knuth ([3]). It makes use of
a strong
property of a normally distributed random vector. We briefly
recall the
definition.
Definition. A random variable has distribution N (0, 1) if it has
the density
function
f (x) =
1√
2π
e−
1
2
x2 .

A d-dimensional random vector X has distribution N (0, I) if the
components
are independant and have distribution N (0, 1) each. In this
case, the density
of X is given by
f (x) =
1
(
√
2π)d
e−
1
2
<x,x>.
In fact, the latter condition is also sufficient for X having
independant
and normally distributed components. This follows from the
uniqueness of
the Fourier transform and the particular form of the Fourier
transform of the
normal distribution (see [1]).
Theorem. Let X be a d-dimensional random vector with
distribution
N (0, I) and U ∈ Rd×d an orthogonal matrix (i.e. U U t = U tU
= I). Then
Y = U X has distribution N (0, I), too.
Proof. For any measurable set A ⊂ Rd we have

P (Y ∈ A) = P (X ∈ U tA)
=
∫
U tA
1
(
√
2π)d
e−
1
2
<x,x>
=
∫
A
1
(
√
2π)d
e−
1
2
<U x,U x>

=
∫
A
1
(
√
2π)d
e−
1
2
<x,x>
0 5 10 15 20 25 30 35 40
0
5
10
15
20
25
N

tim
e
[
s]
Monte Carlo algorithm
Coordinate algorithm
Normal distribution
Figure 1: The performance of the three algorithms
by orthogonalty of U , hence the conclusion follows. 2
As any rotation is in fact just a multiplication with an
appropriate or-
thogonal matrix, we conclude from this theorem that normally
distributed
random vectors are invariant under rotation. Thus, generating X
∈ RN +1
with distribution N (0, I) and then projecting it onto the sphere
SN produces
random vectors U = ‖X‖−1 ·X that are uniformly distributed on
the sphere.
For generating the random vector X, we can make use of the
Box-Muller
method (see [2]). Thus, the time time comlexity of the
algorithm is clearly
linear.
5 Comparison of the algorithms
As the latter of our algorithms neither rejects points nor
involves slow ap-
proximation methods, we expect it to be the most efficient one.

This is
perfectly affirmed by Fig. 1. The figure displayes the time (in
seconds) it
took to generate 1000 random points on the N -sphere with
varying N . As
the experiments were carried out in Matlab, the data is a little
biased: While
the normal distribution algorithm profits from the very fast built
in generator
for normally distributed random numbers, the code for the
incomplete beta
function is very slow. Nontheless, the relations are correct. The
time used
by the Monte Carlo algorithm drastically increases at about N =
10, the
algorithm is not usable for larger N . The coordinate algorithm
takes linear
time. The same does the normal distribution algorithm, but with
a much
smaller ascent.
6 Conclusions
For solving the problem of generating uniformly distributed
random points
on a high dimensional unit sphere in practice, the last algorithm
is obviously
the most efficient and therefore to be preferred. However, if a
specific non-
uniform distribution is needed, a modification of the coordinate
alogrithm
may be appropriate. Even a modification of the Monte Carlo
algorithm can
be suitable, even though it will probably result in an inefficient

algorithm for
high dimensions.
References
[1] H. Bauer. Wahrscheinlichkeitstheorie. deGruyter, 4th
edition, 1991.
[2] G. Box and M. Muller. A note on the generation of random
normal
deviates. Annals of Mathematical Statistics, 29:610–611, 1958.
[3] D. E. Knuth. The Art of Computer Programming, vol. 2:
Seminumerical
Algorithms. Addison-Wesley, 1969.
[4] I. Rechenberg. Evolutionsstrategie ’94. frommann-holzboog,
Stuttgart,
1994.
[5] D. Rusin. Topics on sphere distributions.
http://www.math.niu.edu/ rusin/known-math/95/sphere.faq.
[6] D. Seaman. Topics on sphere distributions.
http://www.math.niu.edu/ rusin/known-math/96/sph.rand.
[7] A. Zell. Simulation neuronaler Netze. Addison-Wesley,
Bonn, 1994.
C2
0.85
C1
0.95

Three Different Algorithms for GeneratingUniformly Distribut.docx

Recommended

More Related Content

Similar to Three Different Algorithms for GeneratingUniformly Distribut.docx (20)

More from herthalearmont (20)

Recently uploaded (20)

Three Different Algorithms for GeneratingUniformly Distribut.docx