27

I'm working on a data mining algorithm where i want to pick a random direction from a particular point in the feature space.

If I pick a random number for each of the n dimensions from [-1,1] and then normalize the vector to a length of 1 will I get an even distribution across all possible directions?

I'm speaking only theoretically here since computer generated random numbers are not actually random.

Jim Lewis
  • 43,505
  • 7
  • 82
  • 96
Matt
  • 1,513
  • 3
  • 16
  • 32

5 Answers5

50

One simple trick is to select each dimension from a gaussian distribution, then normalize:

from random import gauss

def make_rand_vector(dims):
    vec = [gauss(0, 1) for i in range(dims)]
    mag = sum(x**2 for x in vec) ** .5
    return [x/mag for x in vec]

For example, if you want a 7-dimensional random vector, select 7 random values (from a Gaussian distribution with mean 0 and standard deviation 1). Then, compute the magnitude of the resulting vector using the Pythagorean formula (square each value, add the squares, and take the square root of the result). Finally, divide each value by the magnitude to obtain a normalized random vector.

If your number of dimensions is large then this has the strong benefit of always working immediately, while generating random vectors until you find one which happens to have magnitude less than one will cause your computer to simply hang at more than a dozen dimensions or so, because the probability of any of them qualifying becomes vanishingly small.

Eugene Osovetsky
  • 6,443
  • 2
  • 38
  • 59
Bram Cohen
  • 516
  • 5
  • 3
  • Nice! Thank you for the additional suggestion. – Matt Dec 10 '11 at 06:23
  • 1
    By the way, this is how boost http://www.boost.org/doc/libs/1_47_0/boost/random/uniform_on_sphere.hpp implements it. ;) – Jorge Leitao Nov 13 '12 at 17:36
  • I think this is the best answer! However, there is a small chance that you will get the zero vector (and thus a divide-by-zero error) every once in a while – Kip Feb 27 '13 at 17:01
  • 3
    Here's a reference on why this method is right http://mathworld.wolfram.com/HyperspherePointPicking.html – yoavram Oct 09 '13 at 14:43
  • 9
    A quick explanation of why this works: The probability P of a point being at a given is P(x)*P(y). The gaussian distribution has roughly the form e^(-x^2), so e^(-x^2)*e^(-y^2) is e^(-(x^2+y^2)). That is a function only of the distance of the point from the origin, so the resulting distribution is radially symmetric. This generalizes easily to higher dimensions. – mindoftea Aug 20 '14 at 19:05
  • 2
    Additional note: the Box–Muller transform may be used to generate independent pairs of normally distributed variables from independent pairs of uniformly distributed ones (with no 'waste'). – Rhubbarb Oct 13 '14 at 15:47
  • with `numpy` this would be `vec = numpy.random.randn(dims); return vec / numpy.linalg.norm(vec)` – stav Jun 16 '17 at 22:16
  • You can get a division by 0. What about a safer way? – Powereleven Dec 31 '21 at 07:30
  • @Powereleven just use mag=sqrt(dims), which is the expectation of mag. The result may have a slightly different norm than 1 (though the error moves towards zero with higher dims) – Amit Portnoy Feb 02 '22 at 16:19
  • @Powereleven Another thing is that each x**2 is an estimate of the normal's variance (1) and the sum is very unlikely to be zero... – Amit Portnoy Feb 02 '22 at 16:28
13

You will not get a uniformly distributed ensemble of angles with the algorithm you described. The angles will be biased toward the corners of your n-dimensional hypercube.

This can be fixed by eliminating any points with distance greater than 1 from the origin. Then you're dealing with a spherical rather than a cubical (n-dimensional) volume, and your set of angles should then be uniformly distributed over the sample space.

Pseudocode:

Let n be the number of dimensions, K the desired number of vectors:

vec_count=0
while vec_count < K
   generate n uniformly distributed values a[0..n-1] over [-1, 1]
   r_squared = sum over i=0,n-1 of a[i]^2
   if 0 < r_squared <= 1.0
      b[i] = a[i]/sqrt(r_squared)  ; normalize to length of 1
      add vector b[0..n-1] to output list
      vec_count = vec_count + 1
   else
      reject this sample
end while
Jim Lewis
  • 43,505
  • 7
  • 82
  • 96
  • That's what I was worried about. I just wasn't able to formalize it in my head the way you described. Intuitively I know that I want my possible random vectors to describe a circle. I'm just not seeing how to implement it in code. – Matt Jun 08 '11 at 18:06
  • @Matt: I expanded my answer a bit, hope that helps. – Jim Lewis Jun 08 '11 at 18:19
  • Why would you use an algo with non-deterministic run time AND a branch if you could solve this with a closed-form expression? – John Shedletsky Mar 13 '13 at 00:49
  • In high dimensions, this is extremely inefficient. In six dimensions, for example, only 8% of samples will be accepted. In ten dimensions, this falls to 0.25%. – trbabb Dec 09 '21 at 19:34
2

There is a boost implementation of the algorithm that samples from normal distributions: random::uniform_on_sphere

killogre
  • 1,730
  • 15
  • 26
2

I had the exact same question when also developing a ML algorithm.
I got to the same conclusion as Jim Lewis after drawing samples for the 2-d case and plotting the resulting distribution of the angle.

Furthermore, if you try to derive the density distribution for the direction in 2d when you draw at random from [-1,1] for the x- and y-axis ,you will see that:

f_X(x) = 1/(4*cos²(x)) if 0 < x < 45⁰
and
f_X(x) = 1/(4*sin²(x)) if x > 45⁰

where x is the angle, and f_X is the probability density distribution.

I have written about this here: https://aerodatablog.wordpress.com/2018/01/14/random-hyperplanes/

cmhteixeira
  • 877
  • 11
  • 22
-3
#define SCL1 (M_SQRT2/2)
#define SCL2 (M_SQRT2*2)

// unitrand in [-1,1].
double u = SCL1 * unitrand();
double v = SCL1 * unitrand();
double w = SCL2 * sqrt(1.0 - u*u - v*v);

double x = w * u;
double y = w * v;
double z = 1.0 - 2.0 * (u*u + v*v);
MC Pee Pants
  • 9
  • 1
  • 3
    Code no more is hard to read for non-mechanical guys like me. Any comments on what it does, why it's better than the accepted answer, or something ? – Nikana Reklawyks Oct 27 '12 at 00:45