Is there an efficient way to sample uniformly points from the unit n-sphere? Informally, by "uniformly" I mean the probability of picking a point from a region is proportional to the area of that region on the surface of the sphere. Formally, I guess I'm referring to the Haar measure.

I guess "efficient" means the algorithm should take poly(n) time. Of course, it's not clear what I mean by an algorithm since real numbers cannot be represented on a computer to arbitrary precision, so instead we can imagine a model where real numbers can be stored, and arithmetic can be performed on them in constant time. Also, we're given access to a random number generator which outputs a real in [0,1]. In such a model, it's easy to sample from the surface of the n-hypercube in O(n) time, for example.

If you prefer to stick with the standard model of computation, you can consider the approximate version of the problem where you have to sample from a discrete set of vectors that $\epsilon$-approximate the surface of the n-sphere.

  • 3
    $\begingroup$ I would just like to point out that "real numbers cannot be represented on a computer to arbitrary precision" is a common falacy. It is equivalent to "functions on natural numbers cannot be represented on a computer exactly", which is obvious nonsense. $\endgroup$ Commented May 15, 2010 at 6:34
  • $\begingroup$ I don't understand what you mean. Obviously some reals can be represented, but those just form a measure 0 set. Most reals cannot be represented finitely. $\endgroup$
    – Rune
    Commented May 15, 2010 at 12:02

4 Answers 4


Generate $X_1, X_2, \ldots, X_n$ independent, normally distributed random variables See wikipedia for information on how to do this given some standard source of randomness, for example uniform(0,1) random variables.

Then let $Y_i = X_i/\sqrt{X_1^2 + \cdots + X_n^2}$ for $i = 1, \ldots, n$. Then $(Y_1, \ldots, Y_n)$ is uniformly distributed on the surface of the sphere. The time this takes is linear in $n$.

This works because the multivariate normal $(X_1, \ldots, X_n)$ with covariance matrix the identity (that is, $n$ independent unit normals) is rotationally symmetric around the origin.

  • 6
    $\begingroup$ BTW, normal is the only distribution that upon taking $n$-fold product induces a rotationaly invariant distribution. So in some sense, this is the minimal method (i.e. you only need to build a single machine: one for outputting normals). $\endgroup$ Commented May 15, 2010 at 6:28
  • $\begingroup$ For four dimensions, there is a more efficient way. See mathworld.wolfram.com/HyperspherePointPicking.html $\endgroup$ Commented May 14, 2015 at 18:38
  • 1
    $\begingroup$ @Paul Yuryev that's really interesting, do you happen to know of a reference where I can find a proof or at least intuition? $\endgroup$ Commented Aug 1, 2016 at 16:45
  • 3
    $\begingroup$ @Christian Bueno for the product of a distribution with itself to be rotationally symmetric means that f(x)f(y)=f(sqrt(x^2+y^2)) for all x,y in R. Setting g(x)=f(sqrt(x)), this implies g(a)g(b)=g(a+b) for all non-negative real a,b. The only solutions to this are exponentials, so we must have f(x) = e^(cx^2) for some c in R. $\endgroup$ Commented Feb 7, 2022 at 4:41
  • $\begingroup$ @JonathanLove that's a surprisingly concise and elegant argument. Thanks for taking the time to share! $\endgroup$ Commented Feb 7, 2022 at 8:01

You can just take $n$ independent Gaussians and divide the resulting vector by its absolute value. If you have a random number generator that provides you with the standard Gaussian, you can just use it. If all you can generate is a random variable uniformly distributed $[0,1]$, just take two copies $X$ and $Y$ and let $Z=\cos(\pi X)\sqrt{-\log Y}$ (this is a bit simpler than to try to invert erf directly). $Z$ will be a Gaussian with mean $0$ and some variance that doesn't really matter because you normalize in the end anyway.

  • $\begingroup$ Is there a proof/intuitive explanation as to why this works? $\endgroup$ Commented Aug 3, 2016 at 20:15

There is a completely different way of doing this based on Hurwitz's 1897 "factorization of measures". Essentially, you make sequential rotations in $(12)$-plane, then in $(13)$-plane,..., and finally $(n-1,n)$-plane. You can calculate explicitly the $(ij)$-distribution in terms of beta functions. The result is a random rotation in $SO(n)$ which in particular gives you a uniform point in the sphere. See section 2 in this paper by Diaconis and Saloff-Coste for explicit formulas and references.


In case you don't want to worry about dividing by small norms (admittedly, less and less of an issue for higher $n$...)

Generate a point $\mathbf{x}$ on the $(n-1)$-sphere; and generate a number $y\in [-1,1]$ with density $k (1-y^2)^{(n-2)/2}$ for the appropriate constant $k$. Letting $\mathbf{x}'=(y,(1-y^2)^{1/2} \mathbf{x})$ is a uniform point on the $n$-sphere.

Alternatively, in the unlikely event that you don't mind funny correlations between consecutively generated points, you might prefer randomly generating elements of the orthogonal group, for instance letting $A_n$ be sufficiently random${}^1$ antisymmetric matrices, $O_n=(I-A_n)^{-1}(I+A_n)$ is orthonormal, and then the sequence of matrices $U_0=I$, $U_{n+1} = O_n U_n$ is uniformly distributed w.r.t. Haar measure. If later you want to kill consecutive correlations, generate a bunch of them and then sample the sequence at random.

${}^1$: "Sufficiently random": Leveraging the mean-value theorem for all it's worth, we might even take the $A_n$s concentrated on a "generic curve in general position". For sure this isn't efficient, but it gives an idea of how little is needed to get quantitative convergence to uniformity.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.