It should be intuitively clear that rotating the sphere (or the cube) won’t change the distribution of the random points before projection, hence the distribution of the projected points must be independent of the orientation of the sphere, and hence independent of any particular location on the sphere.
Or in other words, if you take the “dotted” sphere out of the cube afterwards, you won’t be able to tell by the dots which way it was originally oriented within the cube.
The part that makes this work is the rejection aspect.
What would be biased is if you inscribed a cube in the unit sphere. This would require additional transformations to create a uniform distribution. If you simply "throw away" the extra corner bits that aren't used, it won't affect the distribution.
> The advantage of this approach is that it generalizes efficiently to any number of dimensions.
I am unsure about whether this is true. The ratio of a ball’s volume to its enclosing hypercube’s volume should decrease to 0 as dimensionality increases. Thus, the approach should actually generalize very poorly.
Let S = {S_i} be any set of cubes that covers a d-sphere. Choose a point in a cube and an integer i in [0, |S|). Now you have a random point in S. With a judicious choice of S you obtain a uniformly random point in the unit sphere with high probability.
Please forgive me my naivete, but won't generating two random polar coordinates do? I'm bad at math, so I might as well be very very wrong here, but I'd like to know.
Edit: see @srean's excellent explanation why that won't do.
At least when trying to end up with stereographically projected coordinates, in general it seems to be faster to uniformly generate a point in the disk by rejection sampling and then transform it by a radially symmetric function to lie on the sphere, rather than uniformly generating a point in the ball and then projecting outward. For one thing, fewer of the points get rejected because the disk fills more of the square than the ball fills of the cube.
My favorite way to generate random points on a n-dimensional sphere is to just sample n times from a Gaussian distribution to get a n-dimensional vector, and then normalizing that vector to the radius of the desired sphere.
This is exactly the method the article describes as the most common method (though the article uses the more specific “standard normal” rather than the more general “gaussian” when describing the distribution), and notes generalizes efficiently to higher dimensions unlike accept/reject, but, as the article notes, the accept/reject method is more efficient for n=3.
Hmm, I don't buy it. The simplicity of just normalizing some Gaussian random deviates (especially since you generate them two at a time using Box-Muller) seems better than accept-reject. Especially considering that the ratio of the volume of the n-dimensional ball to the volume of [-1, 1]^n tends to zero as n tends to infinity exponentially fast...
Accept-reject methods are nonstarters when the architecture makes branching excessively expensive, specifically SIMD and GPU, which is one of the domains where generating random points on a sphere is particularly useful.
The Box-Muller transform is slow because it requires log, sqrt, sin, and cos. Depending on your needs, you can approximate all of these.
log2 can be easily approximated using fast inverse square root tricks:
constexpr float fast_approx_log2(float x) {
x = std::bit_cast<int, float>(x);
constexpr float a = 1.0f / (1 << 23);
x *= a;
x -= 127.0f;
return x;
}
(conveniently, this also negates the need to ensure your input is not zero)
sqrt is pretty fast; turn `-ffast-math` on. (this is already the default on GPUs) (remember that you're normalizing the resultant vector, so add this to the mag_sqr before square rooting it)
The slow part of sin/cos is precise range reduction. We don't need that. The input to sin/cos Box-Muller is by construction in the range [0,2pi]. Range reduction is a no-op.
For my particular niche, these approximations and the resulting biases are justified. YMMV. When I last looked at it, the fast log2 gave a bunch of linearities where you wanted it to be smooth, however across multiple dimensions these linearities seemed to cancel out.
The range reduction of sin/cos is slow only because of the stupid leftover from the 19th century that is the measuring of the phases and plane angles in radians.
A much better unit for phase and plane angle is the cycle, where range reduction becomes exact and very fast.
The radian had advantages for symbolic computation done with pen and paper, by omitting a constant multiplicative factor in the derivation and integration formulae.
However, even for numeric computations done with pen and paper, the radian was inconvenient so in the 19th century the sexagesimal degrees and the cycles continued to be used in parallel with the radians.
Since the development of automatic computers there has remained no reason whatsoever to use the radian for anything. Radians are never used in input sensors or output actuators, because that can be done only with low accuracy, but in physical inputs and outputs angles are always measured in fractions of a cycle. Computing derivatives or integrals happens much more seldom than other trigonometric function evaluations. Moreover, the functions that are derived or integrated almost always have another multiplicative factor in their argument, e.g. frequency or wavenumber, so that factor will also appear in the derivation/integration formula and the extra multiplication with 2Pi that appears in the derivative (or with 1/2Pi that appears in the integral) can usually be absorbed in that factor, or in any case it can be done only once for a great number of function evaluations. Therefore switching the phase unit of measurement from radian to cycle greatly reduces the amount of required computation, while also increasing the accuracy of the results.
A mathematical library should implement only the trigonometric functions of 2Pi*x, and their inverses, which suffice for all applications. There is no need for the cos and sin functions with arguments in radians, which are provided by all standard libraries now.
For reasons that are completely mysterious for me, the C standard and the IEEE standard of FP arithmetic have been updated to include the trigonometric functions of Pi*x, not the functions of 2Pi*x.
It is completely beyond my power of comprehension why this has happened. All the existing applications want to measure the phases in cycles, not in half-cycles. At most there are some applications where measuring the phases or plane angles in right angles could be more convenient, i.e. those might like trigonometric functions of (Pi/2)*x, but again not even in those cases there is any use for half-cycles.
So with the functions of the updated math libraries, e.g. cospi and sinpi, you hopefully can avoid the slow range reduction, but you still have to add a superfluous scaling by 2, due to the bad function definitions.
Similar considerations apply to the use of binary logarithms and exponentials instead of the slower and less accurate hyperbolic (a.k.a. natural) logarithms and exponentials.
Isn't the density distribution of values going to be higher along the directions pointing in the cube's corners? There's more volume between the sphere and nearby the corners than between the sphere and nearby the faces' centres.
An easy way to make this more efficient is to proceed as normal, but if the point is outside the sphere, run the algorithm again using cyclic xor's of the coordinates. This gives you a free second try without generating new random deviates.
You can't do the XOR in floating point representation, but you can if you do the entire algorithm in fixed point or if you retain the random bits before converting to a floating point value.
This decreases the number of random numbers that need to be generated from ~5.72 to ~3.88.
OP is able to create random points to infinite precision or the spherical cow has as many points on it as you like. Many here don't have the luxury of a Hilbert monitor.
If you convert (non mathematician here!) your sphere into an n-agon with an arbitrarily fine mesh of triangular faces, is the method described by OP still valid. ie generate ...
... now I come to think of it, you now have finite faces which are triangular and that leads to a mad fan of tetrahedrons and I am trying to use a cubic lattice to "simplify" finding a series of random faces. Well that's bollocks!
Number the faces algorithmically. Now you have a linear model of the "sphere". Generating random points is trivial. For a simple example take a D20 and roll another D20! Now, without toddling off to the limit, surely the expensive part becomes mapping the faces to points on a screen. However, the easy bit is now the random points on the model of a sphere.
When does a triangular mesh of faces as a model of a sphere become unworkable and treating a sphere instead as a series of points at a distance from a single point - with an arbitrary accuracy - become more or less useful?
I don't think that will be an issue for IT - its triangles all the way and the more the merrier. For the rest of the world I suspect you normally stick with geometry and hope your slide rule can cope.
What about generating 4 random numbers, each uniformly distributed. Make it a vector and normalize it. Now use it as a quaternion rotation, apply it to to vector (1,0,0)
Can someone tell me how to do the opposite? Given a collection of points on a sphere, how does one calculate the probability those points are randomly distributed? I vaguely recall a (1970's?) paper by John E. Westfall that used "nearest neighbor analysis" to determine if the craterlets on the floor of the Moon's Plato crater were randomly distributed, but I don't recall the method. Not sure if that method would generalize to a sphere.
25 comments
[ 4.9 ms ] story [ 47.7 ms ] threadMaybe; my first instinct is that there'll be some bias somewhere.
Maybe I'll have some time tonight to play with this in p5js.
Or in other words, if you take the “dotted” sphere out of the cube afterwards, you won’t be able to tell by the dots which way it was originally oriented within the cube.
What would be biased is if you inscribed a cube in the unit sphere. This would require additional transformations to create a uniform distribution. If you simply "throw away" the extra corner bits that aren't used, it won't affect the distribution.
I am unsure about whether this is true. The ratio of a ball’s volume to its enclosing hypercube’s volume should decrease to 0 as dimensionality increases. Thus, the approach should actually generalize very poorly.
Edit: see @srean's excellent explanation why that won't do.
https://observablehq.com/@jrus/stereorandom
At least when trying to end up with stereographically projected coordinates, in general it seems to be faster to uniformly generate a point in the disk by rejection sampling and then transform it by a radially symmetric function to lie on the sphere, rather than uniformly generating a point in the ball and then projecting outward. For one thing, fewer of the points get rejected because the disk fills more of the square than the ball fills of the cube.
The Box-Muller transform is slow because it requires log, sqrt, sin, and cos. Depending on your needs, you can approximate all of these.
log2 can be easily approximated using fast inverse square root tricks:
(conveniently, this also negates the need to ensure your input is not zero)sqrt is pretty fast; turn `-ffast-math` on. (this is already the default on GPUs) (remember that you're normalizing the resultant vector, so add this to the mag_sqr before square rooting it)
The slow part of sin/cos is precise range reduction. We don't need that. The input to sin/cos Box-Muller is by construction in the range [0,2pi]. Range reduction is a no-op.
For my particular niche, these approximations and the resulting biases are justified. YMMV. When I last looked at it, the fast log2 gave a bunch of linearities where you wanted it to be smooth, however across multiple dimensions these linearities seemed to cancel out.
A much better unit for phase and plane angle is the cycle, where range reduction becomes exact and very fast.
The radian had advantages for symbolic computation done with pen and paper, by omitting a constant multiplicative factor in the derivation and integration formulae.
However, even for numeric computations done with pen and paper, the radian was inconvenient so in the 19th century the sexagesimal degrees and the cycles continued to be used in parallel with the radians.
Since the development of automatic computers there has remained no reason whatsoever to use the radian for anything. Radians are never used in input sensors or output actuators, because that can be done only with low accuracy, but in physical inputs and outputs angles are always measured in fractions of a cycle. Computing derivatives or integrals happens much more seldom than other trigonometric function evaluations. Moreover, the functions that are derived or integrated almost always have another multiplicative factor in their argument, e.g. frequency or wavenumber, so that factor will also appear in the derivation/integration formula and the extra multiplication with 2Pi that appears in the derivative (or with 1/2Pi that appears in the integral) can usually be absorbed in that factor, or in any case it can be done only once for a great number of function evaluations. Therefore switching the phase unit of measurement from radian to cycle greatly reduces the amount of required computation, while also increasing the accuracy of the results.
A mathematical library should implement only the trigonometric functions of 2Pi*x, and their inverses, which suffice for all applications. There is no need for the cos and sin functions with arguments in radians, which are provided by all standard libraries now.
For reasons that are completely mysterious for me, the C standard and the IEEE standard of FP arithmetic have been updated to include the trigonometric functions of Pi*x, not the functions of 2Pi*x.
It is completely beyond my power of comprehension why this has happened. All the existing applications want to measure the phases in cycles, not in half-cycles. At most there are some applications where measuring the phases or plane angles in right angles could be more convenient, i.e. those might like trigonometric functions of (Pi/2)*x, but again not even in those cases there is any use for half-cycles.
So with the functions of the updated math libraries, e.g. cospi and sinpi, you hopefully can avoid the slow range reduction, but you still have to add a superfluous scaling by 2, due to the bad function definitions.
Similar considerations apply to the use of binary logarithms and exponentials instead of the slower and less accurate hyperbolic (a.k.a. natural) logarithms and exponentials.
You can't do the XOR in floating point representation, but you can if you do the entire algorithm in fixed point or if you retain the random bits before converting to a floating point value.
This decreases the number of random numbers that need to be generated from ~5.72 to ~3.88.
If you convert (non mathematician here!) your sphere into an n-agon with an arbitrarily fine mesh of triangular faces, is the method described by OP still valid. ie generate ...
... now I come to think of it, you now have finite faces which are triangular and that leads to a mad fan of tetrahedrons and I am trying to use a cubic lattice to "simplify" finding a series of random faces. Well that's bollocks!
Number the faces algorithmically. Now you have a linear model of the "sphere". Generating random points is trivial. For a simple example take a D20 and roll another D20! Now, without toddling off to the limit, surely the expensive part becomes mapping the faces to points on a screen. However, the easy bit is now the random points on the model of a sphere.
When does a triangular mesh of faces as a model of a sphere become unworkable and treating a sphere instead as a series of points at a distance from a single point - with an arbitrary accuracy - become more or less useful?
I don't think that will be an issue for IT - its triangles all the way and the more the merrier. For the rest of the world I suspect you normally stick with geometry and hope your slide rule can cope.
Essentially:
1. Generate uniformly random values u,v between 0 and 1.
2. Then find their spherical coordinates as:
theta = 2 pi u; phi = acos(2 v -1)
Honestly I’m unsure why you’d choose a method any more complicated.
(Unreasonable Effectiveness of Quasirandom Sequences)
has a section dedicated to spheres.
TL;DR: Generate two orthogonal sequences, map with sin/cos.
I think this method is currently state of the art. It provably doesn't clump in higher dimensions.
Also, it is blue enough for Monte Carlo.