A simple way to generate random points on a sphere(johndcook.com) |
A simple way to generate random points on a sphere(johndcook.com) |
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.
The method presented in the parent article is derived exactly from this method of Marsaglia, and it should also work for any dimension.
Here's an outline of where the θ = 2πu, φ = acos(2v-1) where u, v are uniform random values from 0 to 1 comes from.
If you just picked a uniform random latitude and longitude for each point the distribution would not be uniform, because lines of latitude vary in length. They are smaller the farther they are away from the equator. You need to pick the longer lines of latitude more often than you pick the shorter lines of latitude. The probability density function (PDF) for picking latitude needs to be proportional to latitude length.
If you have a non-uniform distribution and you need to pick an x from it but only have a uniform random number generator there is a trick. Figure out the cumulative density function (CDF) for your PDF. Then to pick a random value use your uniform random number generator to pick a value y from 0 to 1, and find the x such that CDF(x) = y. That's your x. I.e., pick x = CDF^(-1)(y).
Latitude length is proportional to sin(φ) or cos(φ) of latitude depending on whether you are using 0 for the equator or 0 for one of the poles (I think the Mathworld article is using 0 for one of the poles), which makes the PDF proportional to sin(φ) or cos(φ). The CDF is the integral of PDF so ends up proportional to cos(φ) or sin(φ). Using the inverse of that on your uniform random numbers then gives the right distribution for you latitude lengths. Thus we have the acos() in the formula for φ. The argument is 2v-1 rather than v to cover both hemispheres.
You could do the same things in higher dimensions. For a 4D sphere, you could slice it into parallel "lines" of "hyper-latitude". Those would be 2D instead of the 1D of the lines of latitude on a 3D sphere. Then to pick a random point you would chose a line of hyper-latitude at random, and then randomly place a point on that line of hyper-latitude.
Like with the 3D sphere you would need to weigh the lines of hyper-latitude differently depending on where they are on the 4D sphere's surface. Do something similar with the PDF and CDF to use the CDF^(-1) to pick one the lines of hyper-latitude with your uniform random number generator.
You then need to place a point uniformly on that slice of hyper-latitude. That's a 2D thing in this case rather than the 1D we had to deal with before so we will need two more random numbers to place our point. I have no idea what the hell it looks like, but I'd guess it is not going to be "square", so we can't simply use two uniform random numbers directly. I suspect 2D thing would also have to sliced up in parallel lines of "sub-latitude" and we'd have to do the whole PDF and CDF and inverse CDF things there too.
I think in general for an N dimensional sphere you would end up with placing your random point involves picking 1 coordinate directly with your uniform random number generator, and the other N-2 would involve inverse CDFs to get the right distributions from the uniform random number generator.
I have no idea whatsoever if those CDFs would all be simple trig functions like we have in the 3D case, or would be more complicated and possibly not have a reasonably efficient way to compute their inverses.
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;
}
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.
Edit: see @srean's excellent explanation why that won't do.
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.
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.
Maybe; my first instinct is that there'll be some bias somewhere.
Maybe I'll have some time tonight to play with this in p5js.
https://en.m.wikipedia.org/wiki/Moran%27s_I
Create the spatial weights matrix with great circle distances?
(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.
[1] https://docs.scipy.org/doc/scipy/reference/generated/scipy.s...
I’d also point out that the usual way Gaussians are generated (sample uniformly from the unit interval and remap via the Gaussian percentile) can be expressed as sampling from a d-dimensional cube and then post-processing as well, with the advantage that it works in one shot. (Edit: typo)
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.
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.
If you want fast sqrt (or more generally, if you care at all about not getting garbage), I would recommend using an explicit approx sqrt function in your programming language rather than turning on fastmath.
And what's the mathematical basis? (that is, is this technique formalized anywhere?)
It seems insane to me that you run Newton's algorithm straight on the IEEE 754 format bits and it works, what with the exponent in excess coding and so on
It turns out you do have an ultra fast way to compute log_2: you bitcast a float to an integer, and then twiddle some bits. The first 8 bits (after the sign bit, which is obviously zero because we're assuming our input is positive) or whatever are the exponent, and the trailing 23 bits are a linear interpolation between 2^n and 2^(n+1) or whatever. exp_2 is the same but in reverse.
You can simply convert the integer to floating point, multiply by -.5, then convert back to integer. But multiplying -.5 by x can be applied to a floating point operating directly on its bits, but it's more complicated. You'll need to do some arithmetic, and some magic numbers.
So you're bitcasting to an integer, twiddling some bits, twiddling some bits, twiddling some bits, twiddling some bits, and bitcasting to a float. It turns out that all the bit twiddling simplifies if you do all the legwork, but that's beyond the scope of this post.
So there you go. You've computed exp_2(-.5 log_2 x). You're done. Now you need to figure out how to apply that knowledge to the inverse square root.
It just so happens that 1/sqrt(x) and exp(-.5 log x) are the same function. exp(-.5 log x) = exp(log(x^-.5)) = x^-.5 = 1/sqrt(x).
Any function where the hard parts are computing log_2 or exp_2 can be accelerated this way. For instance, x^y is just exp_2(y log_2 x).
Note that in fast inverse square root, you're not doing Newton's method on the integer part, you're doing it on the floating point part. Newton's method doesn't need to be done at all, it just makes the final result more accurate.
Here's a blog here that gets into the nitty gritty of how and why it works, and a formula to compute the magic numbers: https://h14s.p5r.org/2012/09/0x5f3759df.html
A lot of old-school algorithms like CORDIC went the same way.
There's a related technique to compute exponentials with FMA that's somewhat more useful in ML (e.g. softmax), but it has similar precision issues and activation functions are so fast relative to matmul that it's not usually worth it.
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.
To appreciate why, consider strips along two constant latitudes. One along the Equator and the other very close to the pole. The uniformly random polar coordinates method will assign roughly the same number points to both. However the equatorial strip is spread over a large area but the polar strip over a tiny area. So the points will not be uniformly distributed over the surface.
What one needs to keep track of is the ratio between the infinitesimal volume in polar coordinates dphi * dtheta to the infinitesimal of the surface area. In other words the amount of dilation or contraction. Then one has apply the reciprocal to even it out.
This tracking is done by the determinant of the Jacobian.
This gives an algorithm for sampling from a sphere: choose randomly from a cylinder and then project onto a sphere. In polar coordinates:
sample theta uniformly in (0,2pi)
sample y uniformly in (-1,1)
project phi = arcsin(y) in (-pi,pi)
polar coordinates (theta, phi) define describe random point on sphere
One way to fix the problem is to sample uniformly not on the latitude x longitude rectangle but the sin (latitude) x longitude rectangle.
The reason this works is because the area of a infinitesimal lat long patch on the sphere is dlong x lat x cosine (lat). Now, if we sample on the long x sin(lat) rectangle, an infinitesimal rectangle also has area dlong x dlat x d/dlat sin(lat) = dlong x dlat cos (lat).
Unfortunately, these simple fixes do not generalize to arbitrary dimensions. For that those that exploit rotational symmetry of L2 norm works best.
-First of all, it's intuitive to me that the "candidate" points generated in the cube are randomly distributed without bias throughout the volume of the cube. That's almost by definition.
-Once you discard all of the points outside the sphere, you're left with points that are randomly distributed throughout the volume of the sphere. I think that would be true for any shape that you cut out of the cube. So this "discard" method can be used to create randomly distributed points in any 3d volume of arbitrary shape (other than maybe one of those weird pathological topologies.)
-Once the evenly distributed points are projected to the surface of the sphere, you're essentially collapsing each radial line of points down to a single point on the sphere. And since each radial line has complete rotational symmetry with every other radial line, each point on the surface of the sphere is equally likely to be chosen via this process.
That's not a rigorous proof by any means, but I've satisfied myself that it's true and would be surprised if it turned out not to be.
It looks reasonably random to my eye: https://editor.p5js.org/fancybone/sketches/DUFhlJvOZ
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.
You're correctly pointing out that the values of r won't be uniformly distributed. There will be many more points where the value of r is close to 1 then there will be where the value of r is close to 0. This is a natural consequence of the fact that the points are uniformly distributed throughout the volume, but there's more volume near the surface than there is near the center. That's all true.
But now look at the final step. By projecting every point to the surface of the sphere, you've just overwritten every single point's r-coordinate with r=1. Any bias in the distribution of r has been discarded. This step is essentially saying "ignore r, all we care about are the values of theta and phi."
[0]https://en.wikipedia.org/wiki/Spherical_coordinate_system
What I meant is that there's a specific use case for the rejection sampling algorithm, namely computer graphics, and in that use case the asymptotic behavior is irrelevant because n will never not be 3. What is relevant is that the algorithm is more obvious to non-statisticians than Box-Muller, and Marsaglia's variant in particular is also more efficient. Sqrt, ln, and sincos aren't particularly fast operations even on modern hardware (and the latter two aren't SIMDable either), whereas generating uniform variates is almost free in comparison.
EDIT: Plotting it out as a point cloud seems to confirm your suspicion.
Without scaling: https://editor.p5js.org/spyrja/sketches/7IK_RssLI
Scaling fixed: https://editor.p5js.org/spyrja/sketches/kMxQMG0dj