diff options
| -rw-r--r-- | reference/sphere.faq | 761 | ||||
| -rw-r--r-- | reference/sphere.rand | 157 |
2 files changed, 918 insertions, 0 deletions
diff --git a/reference/sphere.faq b/reference/sphere.faq new file mode 100644 index 0000000..cc60ca5 --- /dev/null +++ b/reference/sphere.faq @@ -0,0 +1,761 @@ + TOPICS ON SPHERE DISTRIBUTIONS + +[Latest version of 1998/02/26 by Dave Rusin, rusin@math.niu.edu. +This file is available by FTP (or gopher) at + ftp://math.niu.edu/pub/papers/Rusin/known-math/95/sphere.faq +but I recommend instead this URL which gives links to some related areas: + http://www.math.niu.edu/~rusin/known-math/index/spheres.html +If you're reading this, send me your comments and suggestions.] + +I'm writing this FAQ because of the frequency with which topics related +to this subject appear in the USENET newsgroups such as sci.math. This +isn't really my main area of interest so some of the comments are +pretty simple-minded; if you are a power user in need of guidance, you +should look elsewhere for suggestions, while conversely if you can +offer expert guidance I'd like to incorporate it here. Some of the +material already here is lifted from other sources; whenever possible +I try to assign the blame^H^H^H^H^H credit to someone else. + +The overall main topic is the question of how to place points on the +sphere, or how to divide up the sphere into regions, in a regular way. +Along the way we'll talk about map-making, Platonic solids, and so on. +The following questions will be addressed (other suggestions welcome): + +Q0. Notation +Q1. What does it mean to place N points "evenly" on a sphere? +Q2. Are there some values of N which are better than others? +Q3. How can you describe locations on a sphere, anyway? +Q4. Can you give a quick approximation to a good distribution? (*) +Q5. How can we improve an approximate solution? (*) +Q6. What if I want randomly to place points with a uniform distribution? +Q7. Can this be generalized to higher-dimensional spheres? +Q8. How is this related to collections of linear subspaces? +Q9. Where can I get the coordinates of the vertices of the Platonic solids? +Q10. How come this sounds like sphere packing? + +(*) Fetch the accompanying file sphere.bas for a BASIC program demonstrating +these techniques. + +Other topics possibly to include in the future: + What's the formula for the volume of the n-ball? n-sphere? + What's this I hear about exotic differentiable structures on spheres? + Can you describe projective space and its relation to spheres? + What are sources for further information? +If you need some information about these topics, let me know; that will be +my clue that it's time to add something to the FAQ. + +At some point I may separate this into a 2-sphere FAQ and a hi-dimensional- +sphere FAQ. If this seems important, let me know and I'll move it up my +priority list. + +dave + +============================================================================== + ANSWERS TO THE FREQUENTLY-ASKED QUESTIONS +============================================================================== + +Q0. Notation +Just in case there is some confusion I'll informally make the definitions +I'll need later. + EUCLIDEAN SPACE (or n-space, or R^n) is the set of n-tuples of +real numbers. We will use the ordinary metric (distance function). Note +that we commonly speak of ourselves as living in 3-space. + The SPHERE of radius r in R^n is the set of points whose +coordinates satisfy x1^2+...+xn^2=r^2; the unit sphere is the sphere of +radius 1. Note that I am explicitly assuming all spheres are centered +at the origin. The interior of the sphere is called the BALL in R^n; +it is the set of points of distance less than r from the origin. (The +ball and the sphere are disjoint; their union is sometimes called the +n-DISK .) If not otherwise specified we are usually referring to the sphere or +ball in R^3. + The DIMENSION of a (reasonable) subset of R^n is the number of +independent directions one can travel in the set. For example, in the +ordinary sphere motion at most points is a combination of movements +East and North, so the sphere in R^3 is two-dimensional (not three). +The ball in R^n is n-dimensional. + Note that the circle is the sphere (of dimension 1) in R^2; +the sphere in R^1 is two points (+1 and -1). We often write S^n for the +n-dimensional sphere (in R^(n+1).) + A LINE is the set of multiples of a point or vector (other +than the origin). Note in particular that I am referring only to lines +thru the origin. A PLANE is the set of linear combinations of two +non-collinear vectors. (Don't call it a SURFACE ; that's usually +understood to be anything 2-dimensional, including planes, the +2-sphere, paraboloids, and so on.) Likewise we can take any k +linearly independent vectors in R^n and look at all linear +combinations of them; this set is called a LINEAR SUBSPACE of R^n and +is of dimension k. When k=n-1, the set is called a HYPERPLANE. + + The focus of this document is on S^2 but from time to time I +consider greater generality when it is easy to do so. + +============================================================================== + +Q1. What does it mean to place N points "evenly" on a sphere? + +Well, what do you want it to mean? Several valid interpretations exist; +they lead to slightly different distributions of points. + +If you want all distances between distinct points in the set to be equal, +you can have at most 4 points. Indeed, it's easy to see any 3 +equally spaced points mark an equilateral triangle and any 4 must then +mark a regular triangular pyramid. You can have a 5th point equidistant +from the first 3 (glue two such pyramids face to face) but then the +4th and 5th points are too far apart. + +But everyone agrees that the vertices of a regular polygon are equally +spaced in S^1, so there must be another definition we have in +mind. One can ask that every point have some property, such as being +equidistant from its nearest neighbor. Such distributions exist -- for +example, one can space the points equally along the equator of the +sphere. Still, this is not likely to be what we have in mind when we +ask for a uniform distribution of the points. + +One might ask that the points be arranged "symmetrically" on the sphere in +some way; in Q3 we'll discuss how to do this, but true symmetry is available +only for a few values of N. + +It is more helpful to measure uniformity not point-by-point but by considering +some aggregate measure of the dispersion. For example, we could let D be +the minimum distance between any pair of distinct points in the set. This +provides a measurement of how good or poor a distribution we've chosen. +For example, four points equally spaced on the equator have D=sqrt(2)=1.41... +whereas four points on the pyramid have a larger D=(2/3)sqrt(6)=1.63... +So we might ask for an arrangement on the sphere which maximizes D. + +In general this statement of the problem guarantees a solution exists. +Let D be any continuous function of N points on the sphere; then since +the sphere is compact, so is the Cartesian product of N copies of it, +so that D is a continuous function on a compact set, and so achieves a +maximum and a minimum. If D is chosen intelligently, the values of the +variables that lead to an extreme value of D will give an arrangement of +points which is arguably "uniform". + +There are some drawbacks with this perspective. For one, the extreme +values of D will be attained more than once. (Certainly if D depends +symmetrically on the points, any permutation of the N variables will +provide another occurrence of the extreme value. Indeed, it is +reasonable to consider only functions defined on the symmetric product +S^N(S) =(SxSx...xS)/Sym_n rather than on the Cartesian product +(SxSx...xS) of N copies of the sphere S) If D is invariant under +rotations of the sphere, which is reasonable, then in order to hope to +have an isolated extreme configuration, one needs to hold one point +fixed at, say, the North pole, and to allow another point to vary only +along a meridian attached to the pole. + +Another obstacle is that there are several natural choices for D, and an +arrangement of points which is optimal for one such function need not be +optimal for another. Here are some natural choices for D: + D1=min_i min_{j<>i} dist( Pi, Pj ) (the function considered above) + D2=(1/N) Sum_i min_{j<>i} dist( Pi, Pj ) ("average dist to a neighbor") + D3=Sum_i Sum_{i<>j} 1/dist( Pi, Pj ) ("potential energy") +If an electron is placed at each of the points, D3 is the sum of the +electrostatic potentials at these N electrons; they would tend to spread +apart if possible so as to minimize D3. + +Finally, there is the question of how one would compute an optimal +configuration once a decision is made of how optimality would be measured. +If D is a differentiable function of the coordinates of the points, then +the problem can be solved with Lagrange multipliers, but this is hopelessly +cumbersome for any but the smallest values of N and the simplest +functions D. Moreover, some of the most natural functions D to take are +not differentiable everywhere. We shall return to this in Q5. +============================================================================== + +Q2. Are there some values of N which are better than others? + +You bet! As has been pointed out above, values of N <= 4 have very natural +solutions. In addition, there are solutions for some N related to the +Platonic solids. Here's why. + +Suppose we had an arrangement of N points with the following property: +given any two points there is a rotation of the sphere which + (1) carries the one to the other and in the process + (2) carries each of the points in the distribution to some other point + in the set (possibly itself). +These rotations generate a group (a subgroup of SO(3) ) of +permutations of the N points; since the only rotation of the sphere +which fixes two points not antipodal to each other is the identity +(the trivial rotation), the group of all these rotations is finite +(unless N=2). + +Well, as it happens, there are not many finite subgroups of SO(3). + +First, there are all the cyclic subgroups which are sets of rotations +with a common axis, as well as the dihedral groups which adjoin one +rotation about an axis perpendicular to this one. If this is the +symmetry group of your arrangement of points, then all your points +lie on one circle (two circles, if you have a dihedral group and a +non-equatorial arrangement). + +The remaining finite subgroups of SO(3) are the symmetry group of the +pyramid (isomorphic to Alt(4), it has order 12), the symmetry group of +the cube (it's the full symmetric group of order 24, permuting the 4 +diagonals) and the symmetry group of the dodecahedron (the symmetric +group of order 120). The octahedron is dual to the cube and the +icosahedron is dual to the dodecahedron, while the tetrahedron +(pyramid) is self-dual. (Here "dual" means vertices of one are +midpoints of faces of the other, so dual polyhedra have the same +symmetry group.) (For a full analysis of the Platonic solids, see +Coxeter, "Introduction to Geometry", Wiley (NY) 1969.) + +It's not quite true that a set of points permuted by one of these three +groups has to be the set of vertices of a Platonic solid (as a glance +at a soccer ball [non-US: football] will show). However, a firm +insistence that there be sufficient rotations to carry each vertex +to any other limits the number of points to 120. + +So there are only finitely many non-circular arrangements with this high +degree of symmetry. Of course, you could relax the condition that there be +rotations preserving the set of points which permute the points +transitively (any point can be carried to any other one). For example, you +could begin with one of the Platonic solids and mark off, say, the +midpoint of each edge as well as each of the vertices. (The former need +to be scaled of course to make them lie on the sphere.) There is a large +collection of figures that exhibit some partial symmetry. In general, the +more symmetry you require, the fewer the number of values of N which +you can expect to place so symmetrically. + +If you are working with spheres of a higher dimension, there are more +finite groups of rotations; the orbits of a point under these groups will +give sets of points which could reasonably be called "uniform". Indeed, +given any set of points on a sphere (or any other metric space) you can +define the Voronoi region of each point to be the set of elements in the +sphere closer to that point than to any other. Then if the rotation group +acts transitively on the points, it provides enough motions to show that +all the Voronoi regions are congruent. A couple of theorems are relevant +here, I guess: first, given any group G there is a group of rotations +in SO(n) isomorphic to G as long as n is sufficiently large. Second, +given any n there may be infinitely many finite subgroups in SO(n) +(as the cyclic subgroups of SO(3) demonstrate) but there are only finitely +many quotient groups G / Z(G); that is, the G will have an (Abelian) +center of bounded index. + +Departing from these geometric considerations, we will see in Q4 that +some approximate constructions of uniform arrangements are a little +easier for some N's than others, so if you have some freedom in the +selection of N you can make your job a little easier. + +============================================================================== + +Q3. How can you describe locations on a sphere, anyway? + +Clearly any point on a sphere also sits in the ambient Euclidean space, +and so may be specified by its ordinary Cartesian coordinates. But these +coordinates are bound by the equation of the sphere, and it is sometimes +preferable to use coordinates which are unbound. + +Since the sphere in R^n is (n-1)-dimensional, it is natural to try to +parameterize your position with just n-1 variables. This is well known +on the ordinary 2-sphere, where we use latitude and longitude. Such +parameterizations are possible, although there are some nettlesome points +to keep in mind. + +Ideally, a parameterization of the n-sphere would be a function P +defined on a subset of Euclidean space R^n and taking values in the +n-sphere, having the following properties: + 1) P is continuous (i.e. a continuous change in R^n should yield + a continuous change on the sphere. Equivalently, P + should be given by (n+1) continuous coordinate functions + whose squares sum to 1. ) + 2) P is one-to-one: a point on the sphere shouldn't show up more + than once during the parameterization. + 3) P has a continuous inverse F, that is, a coordinate function + which determines for each point in the sphere what values + of the parameters are needed to map to that point. + The existence of such an F implies P would be onto. + 4) P (and/or F) preserves important geometric information like + angles or areas (or the n-dimensional equivalents). + 5) P (and/or F) is differentiable. +Unfortunately, for topological reasons there can be no function P +meeting even just (1) (2) and (3). This is well-known to map-makers +trying the case n=2: any map of the earth's surface has to rip the +surface somewhere. + +However, there are functions which meet some of these conditions; various ones +are appropriate under different conditions. + +1. If you only want to parameterize _part_ of the sphere, you +can get away with a simple projection: use F(x1,...,x{n+1})=(x1,...,xn) +and P(x1,...,xn)=(x1,...,xn, sqrt(1-x1^2-...-xn^2) ). Use the negative +square root to parameterize the bottom half of the sphere. If you need +to have a parameterization valid a neighborhood of some point on the +equator, have F drop out some other coordinate (which P puts +back in with a similar square root formula). + +2. You can get a correspondence P mapping onto all of the sphere +except one point (say, the North pole "NP"); indeed you can define +P : R^n --> S-{NP} and an inverse map F going the other way which are +continuous, differentiable, one-to-one, and onto. Do this with +stereographic projection: place the plane thru the equator of the +sphere and draw rays emanating from the North Pole and pointing +downward some non-zero amount. Each such ray intersects both the +sphere and the plane in one point, and all the points in S-{NP} and +R^n are so met. These rays then set up the one-to-one +correspondence. The formulas are + P(x1,...,xn) = t.(x1,...,xn,0) + (1-t)(0,...,0,1) + F(x1,...,x{n+1})= r . (x1, ..., xn) +(where (0,...,0,1) is the North Pole on the n-sphere in R^{n+1}.) Here +t is chosen to make the distance from the origin be 1 : t=2/(1+x1^2+...+xn^2). +Also r is chosen to invert the first formula: r=1/(1-x{n+1}) + +These formulas are the basis of the polar projection used by map-makers. + +3. If it is necessary to have a parameterization which preserves area +(or its n-dimensional equivalent), this coordinate system F can be +modified by adjusting the function r=r(x{n+1}). + +Here's how we decide what the right function r is. Consider the set of +points in S^n with x{n+1}=a. This is an (n-1) sphere with radius +h=sqrt(1-a^2). Then for small e the set of points in S^n with +x{n+1} lying between a and a+e is a slightly bent version of the +product of an (n-1)-sphere of radius h and an interval of length +e.sqrt[1+(dh/da)^2] = e/sqrt(1-a^2). The function F above carries it +to another such product, but this time the radius is z(a)=r(a)h(a) and +the length of the interval is about e.z'(a). Since F changed the radius +by a factor z(a)/h, the (n-1)-volume will change by a factor +of [z(a)/h]^(n-1). If we want to keep the n-volume unchanged, we need +the ratio of the lengths of the interval, roughly z'(a) sqrt(1-a^2), to +be the reciprocal of this, that is, we need z'(a).z(a)^(n-1) = +h^(n-1)/sqrt(1-a^2) = sqrt(1-a^2)^(n-2). Solving this differential +equation shows we need + r(a) = z(a)/h(a) = [1/sqrt(1-a^2)].[n.integral(1-a^2)^(n/2-1) da]^(1/n) +(As z(-1)=0, the integral should be selected to vanish at a=-1.) + +So r(x)=[1/sqrt(1-x^2)].[n. integral_{-1}^{x} (1-a^2)^(n/2 - 1) da]^(1/n). +For n=2, this gives r(a)=sqrt(2/(1-a)), i.e. z(a)=sqrt(2*(1+a)). In +particular, this will map the whole sphere to the disk of radius 2, +whose area is then 4pi (the area of the sphere) as it should be. + +4. There are also trigonometric parameterizations. +For example, it is easy to parameterize the 1-sphere (the circle) with a +parameter t using P(t) = (cos(t), sin(t)). Here if we really want P +to be one-to-one we need to restrict to an interval such as [0, 2 pi). +This is then indeed one-to-one, onto, and continuous (even differentiable) +and as a bonus preserves lengths. But there can be no continuous inverse +function. F(x,y)=arctan(y/x) will do if branches of the arctan function +are chosen appropriately, but no set of choices can be made consistently +in all quadrants of S^1. + +Still, this parameterization is useful, and admits generalization. The +function + P(u,v) = (cos(u)cos(v), cos(u)sin(v), sin(u)) +is also differentiable, maps onto S^2, and is one-to-one on the set +(-pi/2,pi/2)x[0,2 pi). (The image of P on this set misses the poles, +however.) This P has an inverse in the sense we had on the circle: +(u,v)=F(x,y,z)=(arctan(z/sqrt(x^2+y^2)), arctan(y/x) ). + This is the ordinary parameterization of the sphere in terms of +latitude (u) and longitude (v). By computing the Jacobian for the change of +variable to spherical coordinates, we see that this parameterization +changes areas by a factor of |cos(u)|. + +This may be generalized to the n-sphere: use the function P_n, where +P_n(u1,...,un) = ( cos(un)*P_{n-1}(u1,...,u{n-1}) , sin(un) ) +(for each un in (-pi/2, pi/2), the parameterization traces out an +(n-1)-sphere using P_(n-1). ); n-dimensional volume is reduced by a +factor of | cos(u2) . cos(u3)^2 . ... . cos(un)^(n-1) |. +(I'll discuss another interpretation of this product in Q6.) + +These parameterizations form the basis of various cylindrical maps of +the earth (rectangular maps excluding the poles in which the left and +right edges represent the same points on earth). The Mercator projection +is chosen to remove the distortion of angles; thus (locally) the +Mercator projection gives an accurate representation of shapes. But the +correction used is precisely the inverse of what is needed to correct +the |cos(u)| change in area, so that the Mercator projection is wildly +inaccurate at portraying areas (especially near the poles). I believe +the 'Peters projection' is the name of the map which corrects areas (but +which amplifies the distortions of angles, and thus of shapes). + +5. If the trigonometric expressions are too computation-intensive, +one can use a parameterization requiring only rational functions. Note first +that the circle may be parameterized using a variation of stereographic +projection: + P(t) = ( -2t/(1+t^2) , (1-t^2)/(1+t^2) ) + F(x,y)= (y-1) /x +Thus, we can replace all sine/cosine pairs with these two functions. +For example, the 2-sphere may be parameterized: +P(u,v)=( 4uv/[(1+u^2)(1+v^2)], -2u(1-v^2)/[(1+u^2)(1+v^2)], (1-u^2)/(1+u^2) ) +F(x,y,z)= ( (z-1)/sqrt(x^2+y^2), (y-1)/x ) + +For further details, consult an advanced calculus text; look up the +change-of-variables theorem, particularly spherical coordinates and +the role of the Jacobian (derivative) in determining areas. +[I'll add a reference to the mathematics of map-making as soon as I find one.] + +============================================================================== + +Q4. Can you give a quick approximation to a good distribution? + +Of course. Mark off K equal segments along a meridian using K-1 +points plus the two endpoints. These segments have length pi/K. If you +swing this around the sphere, you'll mark off K-1 circles. Now, in +spherical coordinates, these circles are the points with +ui=-pi/2 +i*(pi/K), where i = 1, 2, ..., K-1. They have radius cos(ui), +hence length 2pi*cos(ui). Thus you can divide them into about +[2pi*cos(ui)]/[pi/K] segments as long as the segments on the meridian. +So take Mi to be an integer near 2*K*cos( -pi/2 + i*(pi/K)) = +2*K*sin( i * (pi/K) ), and mark off points that have u=ui and +v=j*(2*pi/Mi) for j = 1, 2, ..., Mi. This will give you +2+Sum{i=1,...,K-1} M_i points on the sphere. + +To within an error less than K-1, this is the same as +2+Sum 2*K*sin( i* pi/K ) = 2 + 2*K*cot( pi / 2K ). For large values of K, +this is about (4/pi) K^2. (With an area of 4pi on the whole sphere, +this leaves an area of (pi/K)^2 around each point, which is not +surprising since initially the points are about pi/K apart on a +roughly square grid). + +If you have the luxury of placing _approximately_ N points on the +sphere, begin with K an integer roughly equal to sqrt((pi/4)N), and +compute the number of points which would be placed with the method +above; you should get about N points. Here are the values I get +for the first few K using round() to convert to integers when necessary: +for K=1, N=2; for K=2, N=6; then for subsequent K (thru K=22), I get N= +12,22,34,46,82,106,128,156,184,214,248,288,328,368,412,460,510,562,616,... + +Notice that the first layer around the north pole has about 2Ksin(pi/K) +points. For K large, this is about 2*pi = 6.28, so we would take M1=6, +and arrange the first circle of points in the "kissing pennies" pattern. + +This arrangement of points on the sphere has an obvious regularity, +but the points are not spaced as far apart as possible by any reasonable +measure -- for example, roughly sqrt(K/5) of the rows nearest the +equator will have about the same number of points in them, roughly in +a square-looking grid; clearly a honeycomb distribution will place the +same number of points with a greater spacing. This leads to question Q5. + +I have included as a separate file a BASIC program which carries out this +distribution. In view of the previous paragraph, this code can be +improved when N is large, but I'm not really sure of an optimal way to +do so, so I won't. + +============================================================================== + +Q5. How can we improve an approximate solution? + +If the approximate solution in Q4 is insufficient, you'll have to decide why. +This returns you to the discussion on Q1. Let me assume that there is some +function D of the positions of the points which you wish to minimize +or maximize. (Replacing D by -D as necessary we may assume D is to be +maximized.) + +The method outlined in Q4 should give a low value of D. If D is +differentiable, you can use the method of gradients: you adjust the positions +of the points in the direction of the gradient. Professionals in +numerical analysis should scoff at this naive approach -- for one thing, +this method assumes that the critical point is an isolated minimum; also, +it is not clear just how far in the direction of the gradient to +proceed so as not to overshoot the region of convergence; moreover, once +in the region of convergence, the rate of convergence can be quite slow. +However, in practice I have found that if you have a moderate number of +points to place and you are not too demanding in your search for an +absolute minimum, this method works quickly enough if you begin with a +moderately uniform distribution of points. (I would appreciate learning +of improvements from numerical analysts.) + +Here is what the method of gradients means in more detail. If you +have parameterized the sphere and expressed D as a function of N sets +of the parameters, compute the partial derivatives of D with respect +to all these variables, and then evaluate at the present values of the +parameters. Add to each of the parameters a small multiple of this +partial derivative. + +If instead you have calculated D as a function of the Cartesian +coordinates of the N points, you need to take the portion of the +gradient which points in a direction tangent to the sphere. In +practice, this can be accomplished by adjusting all coordinates of the +points as in the preceding paragraph, then rescaling each point's +coordinates to ensure the vector is of length one. + +There is the possibility of creep to be concerned with. That is, it is +possible that with each iteration the function D is not improved; +rather, the points have simply all rotated somewhat around the sphere. +Under the assumption that D is invariant under rotations, there is no +need to allow this kind of motion. One can negate this possibility by +rotating the sphere after each iteration so that one point (say the +North pole) is returned to its previous position and another point +is returned to the meridian it was on previously. + +I have put a BASIC program alongside this file to illustrate these +computations. Since the function D I used polls each pair of points +to find out the repulsion between them, the running time for one +iteration is proportional to N^2. Thus I hardly recommend it for +values of N above a couple of hundred if you are using, say, a +personal computer. + +============================================================================== + +Q6. What if I want randomly to place points with a uniform distribution? + +This is perhaps the easiest way: randomly select the Cartesian +coordinates independently. Toss out any sets of coordinates which take +you outside the sphere (and, should you be so unlucky, any occurrences +of the origin). Scale all the remaining points to push them out to the sphere. + +Another easy way is to use the area-preserving parameterizations of +Q3. If you use the cartographer's polar projection (the third +parameterization in Q3), you need only select random points in the +disk with uniform distribution there, then map each of the selected +points onto the sphere. Of course, you'll never hit the North Pole +this way, but that event was supposed to occur with probability zero +anyway. (The simplest way to select random points in the disk is just +to pick the n coordinates separately, each randomly in [-1,1]. Then +simply discard the n-tuples not lying in the disk.) + +You can also use the trigonometric parameterizations of the sphere +(the cartographer's cylindrical projections): if you are placing +points randomly, the probability that you will have to consider one of +the points where the parameterization fails to be one-to-one, or +continuous, say, is zero. + +However, you need to exercise some care for spheres of dimension +greater than 1 since in that case, the trig parameterizations are not +area-preserving. For example, if you use the spherical coordinates +above to parameterize the 2-sphere, and select u and v randomly in +their respective intervals, you'll find too many points chosen near +the poles. This is because the sets of points with u fixed near ++-pi/2 are only _small_ circles; it is not appropriate to select these +u's as frequently. Indeed, the frequency with which a given u is +selected should be proportional to the circumference of that circle, +which is 2pi |cos(u)|. + +More generally, one can randomly select points in the n-sphere by randomly +selecting the n spherical parameters, but always remembering to select +each u_i with a frequency proportional to the (i-1)-dimensional +measure of the the set of points where u_i is a constant. Fortunately, +we know that this set is just an (i-1) sphere whose radius is cos(u_i) +and whose measure is then proportional to cos(u_i)^(i-1). + +Thus the problem of random point selection on the sphere reduces to this +question: how can one select values of a variable u from an interval +[a,b] so as to achieve a given frequency distribution f? One way is to +let C be the cumulative distribution function, C(u)=integral_a^u f(t)dt. +(This is the probability that the chosen value will be between a and u.) +For then if G is the inverse of this function, then we may select +x randomly with uniform distribution across [0,1] and let y=G(x). This is +a random variable whose distribution function looks like this: + Prob( h < y < h+e ) = Prob( h < G(x) < h+e ) + = Prob(C(h) < x < C(h+e) ) which is roughly + = Prob( C(h) < x < C(h)+C'(h)e ) +using differentials. Well, if x is selected with uniform distribution, +this probability is the length of that interval, namely C'(h)e = f(h).e +since C'=f. But to say that this original probability is roughly f(h)e +is precisely to say that the distribution function for the variable y is f. + +With the frequency distributions at hand, this is not difficult. The +first coordinate u1 is selected uniformly in [0,2 pi). The second +coordinate u2 is to have a probability density function cos(u_2) /2 +on (-pi/2, pi/2) [The extra factor of 1/2 is just a scaling needed to +ensure that the cumulative distribution C(u) ends with C(pi/2) = 1.0]. +Indeed, C(u) = (sin(u)+1)/2 has inverse G(w)= arcsin(2w-1), so one can +get points uniformly on the 2-sphere using the usual trig parameterization +with v selected uniformly in [0, 2pi) and u=arcsin(2w-1) where +w is selected uniformly in [0,1]. The Cartesian coordinates are then +(sqrt(4x-4x^2).cos(v), sqrt(4x-4x^2).sin(v), 2x-1). + +With a slight change of variables, this can be viewed as the equal-area +cylindrical projection in cartography: you can parameterize the surface +of the 2-sphere with coordinates x in [-1, 1] and v in [0,2 pi) +in such a way as to preserve area by mapping (x,v) to +(sqrt(1-x^2) cos(v), sqrt(1-x^2) sin(v), x). This is the projection +which carries points on the sphere (minus the poles) to points on a +cylinder wrapped around the equator by shining a light straight out from +the axis through the sphere to the cylinder. + +On the 3-sphere (and higher) the next variable is selected as u3=G(x3) +where x3 is selected uniformly in [0,1] and G is the inverse of +the function C(u) = (1/2)(u+pi/2+sin(u)cos(u)). The next variable +after that is selected similarly using C(u) = ( cos(u) - sin^3(u)/3 - 1/3). +And so on. + +Observe that if C can be computed, it is not necessary to be able to +solve for its inverse G explicitly. Instead, each of the computations of +y=G(x) may be written x=C(y), an equation to be solved for y when x +is known. This is easily handled with Newton's method. + +I believe a method of generating numbers randomly with a given distribution +is treated in Knuth's Art of Computer Programming. + +============================================================================== +Q7. Can this be generalized to higher-dimensional spheres? + +I have melded this into the previous questions -- anyone reading this +who can think of questions particular to higher dimensions should let me +know what ought to go in here. + +============================================================================== + +Q8. How is this related to collections of linear subspaces? + +I put this question here because someone asked me recently about how +to parameterize the set of planes and other linear subspaces of Euclidean +space. I can't quite remember the question, so I'll not make an +answer at this point. + +One comment to make is that points on the sphere (or, more precisely, +pairs of antipodal points) correspond to lines through the origin in +R^n. The planes in R^n correspond similarly to great circles (circles +of maximal circumference). If n=3, the planes can be parameterized in +a natural way by finding the line perpendicular to them. If you +trace these correspondences through, you'll see that there is a +natural one-to-one correspondence between great circles on the 2-sphere +and pairs of antipodal points. You can get the circle from the points +as the set of points equidistant from the pair; you can get the pair from +the circle as the only points equidistant from all elements of the circle. + +More generally, there is a correspondence between the k-dimensional +and the n-k dimensional subspaces of R^n. The family of all k-dimensional +subspaces in R^n (or n-k dimensional subspaces) is known as the +Grassmannian manifold; it is indeed a manifold and has dimension +k(n-k). When k=1 or k=n-1, this is also known as projective n-1 space; +there is a two-to-one mapping S^{n-1} --> RP^{n-1} which simply sends +a point to the line it's on (this map sends antipodal points of S^{n-1} +to the same element of RP^{n-1} of course. This map turns S^{n-1} +into a two-fold cover of RP^{n-1}; for n>2 this is the universal +covering of projective space.) + +This material is often treated in graduate topology courses, particularly +in a discussion on vector bundles in an algebraic topology course. +You may wish to consult Milnor's book "Characteristic Classes" +(Princeton Univ Press), which is however not elementary. + +============================================================================== + +Q9. Where can I get the coordinates of the vertices of the Platonic solids? + +There are 5 Platonic solids. An interesting but brief discussion (with +pictures) may be found, for example, in Gallian's "Contemporary +Abstract Algebra", D C Heath 1986. + +Solid Number of Vertices # Edges # Faces +Tetrahedron 4 6 4 +Cube 8 12 6 +Octahedron 6 12 8 +Dodecahedron 20 30 12 +Icosahedron 12 30 20 + +Notice the duality between certain pairs. Also note that V-E+F=2 in +all cases. This is deep, and may be thought of as the first fact in +homology theory. + +I will give coordinates of these points (and in indication of the rest +of the structure) in the most compact way I know; if you want the +coordinates so that, for example, the North Pole is one of the vertices, +you'll have to rotate the coordinates (the BASIC program I include has +in it an application of this procedure.) +I have scaled the coordinates so that they lie in the unit sphere, but +this is not always the prettiest form. I have however set it up so that +the coordinates of the dual pairs reflect this duality. + +Tetrahedron +The points are (0,0,1), (1/3)(2sqrt(2), 0, -1), and both points +(1/3)(-sqrt(2), +-sqrt(6), -1). All pairs are joined by an edge; all +triples (sets of three points) lie on a common face. + +Cube +The 8 points are (+-c, +-c, +-c) where c=1/sqrt(3) and the plus/minus signs +may be chosen independently. +The 12 edges join pairs of points which have 2 coordinates equal. +The 6 faces are quadruples of points for which one of the coordinates is + constant. + +Octahedron +The 6 vertices are of three types: (+-1, 0, 0) (0, +-1, 0) and (0, 0, +-1). +The 12 edges join all pairs except the three pairs of antipodal points. +The 8 faces consist of all triples which include one point of each type. + +Icosahedron +(using some information in a post by cjhenrich@delphi.com) +The 12 vertices are of three types: (+-a, +-b, 0), (0, +-a, +-b), and + (+-b, 0, +-a), where a=sqrt((5+sqrt(5))/10) and b=sqrt((5-sqrt(5))/10) +The 30 edges come in two flavors: there are 24 edges joining points of + different types such that their common non-zero coordinate is of + the same sign (e.g., (a, -b, 0) and (0, -a, b) have negative 2nd + coordinates). There are also 6 edges joining pairs of points of + the same type for which the sign in front of the a are equal + (e.g. (a, b, 0) and (a, -b, 0).) +The 20 faces are also of two general sorts: There are 8 triplets consisting + of one point of each type such that whenever two of the points + have a common non-zero coordinate, those coordinates are of the + same sign. Also, there are 12 triplets which contain two points of + the same type having the same sign 'X' in front of the a, and a + third point having 'X'.b in that coordinate (e.g., given the + two points at the end of the last paragraph, the third point on + a face with them is (b,0,+-a). ) + +Dodecahedron: This is dual to the icosahedron +I adapted these coordinates from a post by David Seal. +The 20 vertices include 8 vertices (+-c, +-c, +-c) (c=1/sqrt(3)) and + 12 of the forms (+-x, +-y, 0), (0, +-x, +-y), (+-y, 0, +-x) where + x=(sqrt(5)+1)/(2*sqrt(3)) and y=(sqrt(5)-1)/(2*sqrt(3)). +I will leave it to the reader to dualize the construction of the edges + and vertices of the icosahedron to get the edges and faces of + the dodecahedron. (That means I started it but it took too long. + The edges are the pairs of points of distance sqrt( (2/9)(sqrt(5)-1) ) + from each other.) + +============================================================================== + +Q10. How come this sounds like sphere packing? + +Very observant of you! Actually, there are a lot of questions whose +answers sound like an enumeration of the Platonic solids. Here is one +way to see a connection. + +Suppose you have a packing of Euclidean space with spheres of radius 1. +There is a way to get evenly spaced points on the unit sphere from it. + +For any r>0 let a_r be the number of spheres in the packing which are +contained in the ball of radius r centered at the origin. Note that +if V is the volume of the ball of radius 1 in R^n, then the ball of +radius r has volume V.r^n, so we have an upper bound of a_r <= r^n. +The limit of a_r to r^n as r increases without bound is the 'packing +density' d. Thus for large r we have a_r is roughly d r^n. In +particular, there are roughly d.n.r^(n-1) balls which fit in the +sphere of radius r+1 but not in the sphere of radius r. If we take the +centers of just these balls and scale them back just a bit, we will +get roughly d.n.r^(n-1) points on the sphere of radius r. Since the +sphere centers are of distance at least 2 from each other, the points +on the sphere will be almost 2 units apart, or further. + +Scaling this picture by a factor of r, we get about d.n.r^(n-1) points +on the unit sphere, with distances between them not much less than +2/r. This is not really better than we achieved with the method in Q4, +but it has the advantage that we can also place points on the sphere +corresponding to spheres in the packing which, for example, fit in +the ball of radius r+1.5, or r+2. Possibly these points, when projected +to the sphere of radius r, will lie close to other points already +placed, but then again, maybe not. + +It is not possible, so far as I know, to go backwards and create a +sphere packing from a distribution on a sphere. The difficulty is that +a sphere packing implies a distribution of arbitrarily many points on +the sphere, and gives these distributions in a coherent way, whereas +an individual distribution is done with just one fixed number of points. + +There is abundant literature and a fair amount of controversy about +sphere packings. Certainly the easiest sphere packings to understand +and use are the ones with great regularity: one places a sphere +centered at each point in a crystallographic lattice. An introduction +to these may be found, for example, in Durbin's "Modern algebra" +(Wiley, 1985). One looks at a subgroup of the group of translations in +R^n which is discrete (each group element is far from any other one); +the translates of a single point under this group form a lattice. The +group of symmetries of such a lattice is called a crystallographic +group. Note that there may be some symmetries which leave a point +fixed; the group of these is finite, and indeed there are only a few +possible groups: dropping, if necessary, to a subgroup of index two, +we can assume the group contains only rotations, at which point it +must be cyclic of order 1, 2, 3, 4, or 6; dihedral of order 4, 6, 8, +or 12, or the symmetry groups of the tetrahedron or octahedron. + +This topic comes up from time to time in the popular press, during +which an interview with the inimitable John Conway is de rigueur. +There is also some very interesting work on lattices in higher +dimensions, particularly dimensions which are a multiple of 8, such as +the important Leech lattice. The book by Conway and Sloane is full of +such tidbits. + +============================================================================== + diff --git a/reference/sphere.rand b/reference/sphere.rand new file mode 100644 index 0000000..bf56705 --- /dev/null +++ b/reference/sphere.rand @@ -0,0 +1,157 @@ +From: ags@seaman.cc.purdue.edu (Dave Seaman) +Newsgroups: sci.math.num-analysis +Subject: Re: N-dim spherical random number drawing +Date: 20 Sep 1996 13:04:58 -0500 + +In article <wmEgVDG00VUtMEhUwg@andrew.cmu.edu>, +Shing-Te Li <sl2x+@andrew.cmu.edu> wrote: +>Does anyone know a good algorithm that can draw a vector of random +>numbers uniformly from a N dimensional sphere? + +Here are four methods for generating uniformly-distributed points on a unit +sphere: + + (1) The normal-deviate method. Choose x, y, and z, each + normally-distributed with mean 0 and variance 1. (Actually, the + variance does not matter, provided it remains fixed.) Normalize + the result to obtain a uniform density on the unit sphere. + + This method also generalizes nicely to n dimensions. It works + because the vector chosen (before normalization) has a density + that depends only on the distance from the origin. The method is + mentioned in Knuth, _Seminumerical Algorithms_, 2nd. ed., pp. + 130-131. + + (2) The hypercube rejection method. Choose x, y, and z, uniformly + distributed on [-1,1]. Compute s = x^2 + y^2 + z^2. If s > 1, + reject the triplet and start over. Once you have an acceptable + vector, normalize it to obtain a point on the unit sphere. + + This method also generalizes to n dimensions, but as the dimension + increases the probability of rejection becomes high and many random + numbers are wasted. + + (3) The trig method. This method works only in 3-space, but it is + very fast. It depends on the slightly counterintuitive fact (see + proof below) that each of the three coordinates of a uniformly + distributed point on S^2 is uniformly distributed on [-1,1] (but + the three are not independent, obviously). Therefore, it + suffices to choose one axis (Z, say) and generate a uniformly + distributed value on that axis. This constrains the chosen point + to lie on a circle parallel to the X-Y plane, and the obvious + trig method may be used to obtain the remaining coordinates. + + (a) Choose z uniformly distributed in [-1,1]. + (b) Choose t uniformly distributed on [0, 2*pi). + (c) Let r = sqrt(1-z^2). + (d) Let x = r * cos(t). + (e) Let y = r * sin(t). + + (4) A variation of method (3) that doesn't use trig may be found in + Knuth, _loc. cit._. The method is equivalent to the following: + + (a) Choose u and v uniformly distributed on [-1,1]. + (b) Let s = u^2 + v^2. + (c) If s > 1, return to step (a). + (d) Let a = 2 * sqrt(1-s). + (e) The desired point is (a*u, a*v, 2*s-1). + + This method uses two-dimensional rejection, which gives a higher + probability of acceptance compared to algorithm (2) in three + dimensions. I group this with the trig method because both + depend on the fact that the projection on any axis is uniformly + distributed on [-1,1], and therefore both methods are limited to + use on S^2. + +I have found the trig method to be fastest in tests on an IBM RS/6000 +model 590, using the XLF 3.2 Fortran compiler, with the IMSL +subroutines DRNUN and DRNNOA generating the random numbers. The trig +routines were accelerated by using the MASS library, and sqrt +computations were performed by the hardware sqrt instruction available +on the 590. Timings for 1 million random points on S^2: + + (1) normal-deviate method 9.60 sec. + (2) 3-D rejection 8.84 sec. + (3) trig method 2.71 sec. + (4) polar with 2-D rejection 5.08 sec. + +The fact that algorithm (3) does not involve rejection turns out to be +a huge advantage because it becomes possible to generate all the points +in a single loop with no conditional statements, which optimizes very +effectively and more than offsets the cost of using trig functions. +The alternative with algorithms (2) and (4) is to generate the points +one at a time and test each one for acceptance before proceeding to the +next, which requires nested loops and many more subroutine calls to +generate the random numbers. Algorithm (1) does not explicitly use +rejection, but the IMSL subroutine DRNNOA uses rejection in generating +the normally-distributed numbers that are required. + +In higher dimensions, the normal-deviate method is preferred +because the probability of rejection in the hypercube method +grows very rapidly with the dimension. + +-------------------------------------------------------------------- + +Here is a proof of correctness for the fact that the z-coordinate is +uniformly distributed. The proof also works for the x- and y- +coordinates, but note that this works only in 3-space. + +The area of a surface of revolution in 3-space is given by + + A = 2 * pi * int_a^b f(x) * sqrt(1 + [f'(x}]^2) dx + +Consider a zone of a sphere of radius R. Since we are integrating in +the z direction, we have + + f(z) = sqrt(R^2 - z^2) + f'(z) = -z / sqrt(R^2-z^2) + + 1 + [f'(z)]^2 = r^2 / (R^2-z^2) + + A = 2 * pi * int_a^b sqrt(R^2-z^2) * R/sqrt(R^2-z^2) dz + + = 2 * pi * R int_a^b dz + + = 2 * pi * R * (b-a) + + = 2 * pi * R * h, + +where h = b-a is the vertical height of the zone. Notice how the integrand +reduces to a constant. The density is therefore uniform. + +-- +Dave Seaman dseaman@purdue.edu + ++++ stop the execution of Mumia Abu-Jamal ++++ + ++++ if you agree copy these lines to your sig ++++ +++++ see http://www.xs4all.nl/~tank/spg-l/sigaction.htm ++++ +============================================================================== + +From: George Marsaglia <geo@stat.fsu.edu> +Subject: Random points on a sphere. +Date: Tue, 15 Jun 1999 18:03:07 -0400 +Newsgroups: sci.math.num-analysis,sci.math,sci.stat.math + +Questions on methods for random sampling from the surface +of a sphere keep coming up, and seemingly get lost as +new generations of computer users encounter the problem. + +I described the methods that I had used for years in +"Sampling from the surface of a sphere", Ann Math Stat, v43, 1972, +recommending the two that seemed the fastest: + +1) Choose standard normal variates x1,x2,x3, +put R=1/sqrt(x1^2+x2^2+x3^2) and return the point +(x1*R,x2*R,x3*R). + +2) Generate u and v, uniform in [-1,1] until S=u^2+v^2 < 1, +then return (1-2*S,u*r,v*r), with r=2*sqrt(1-S). + +The second method is fast and easiest to program---unless +one already has a fast normal generator in his library. + +Method 1) seems by far the fastest (and easily applies to higher +dimensions) if my latest normal generator is used. +It produces normal variates at the rate of seven million +per second in a 300MHz Pc. + +George Marsaglia |