Sphere Problem
The question of the optimal arrangements of a number of points on the surface of a sphere returns at regular intervals on USENET. Although the earliest references to this problem date back almost a century, the capability to actually compute configurations with larger numbers of points has become available only recently. The various optimization criteria which may be proposed have been found to have multiple stable solutions and very large numbers of meta-stable solutions. In many cases the polyhedra formed by the convex hull of the solutions are new and of interest. Uniform distributions of points on the sphere are also required by certain computer graphics algorithms.

It is natural to seek regular arrangements of points as the solution of some optimality criterion involving the points such as the scalar distances between the points. Most references to the sphere problem consider the problem of minimizing the electrostatic potential of a number of equal charges on the surface of a unit sphere. This can be shown to be equivalent to minimizing the sum of the inverses of all distances: (F = Sum 1/dkl). I will refer to this problem as the electrostatic problem; it is also sometimes known as the Thomson Problem. It is alternatively possible to maximize the average Euclidean distance for a number of points: (F = Sum dkl). I will refer to this problem as the geometric problem. This problem can alternatively be formulated by using the distances along the surface rather than the Euclidean distances, in which case I refer to it as the arclength problem. Instead of considering the distances of a given point to all others, it is alternatively possible to consider its interactions with only its nearest neighbours by maximizing the minimum distance between points but this is expected to be a much harder problem. This maximum minimum distance problem can also be formulated using either the Euclidean distance or the distance along the surface which seems very natural. Finally, there are other related problems such as the maximum volume of an inscribed polyhedron. This problem is in a different class since its objective function is not expressed in terms of the distances between the points.

My friend Gary G. Ford who has been a sounding board for many of the ideas in this paper has consistently advocated the use of the geometric objective function. He contributed the analytical result that the average distance between infinitely many points distributed uniformly on a sphere is 4/3 times the sphere radius. He has also noted that the geometric problem and the electrostatic problem both maximize certain definitions of the mean of the distances: the geometric problem maximizes the arithmetic mean whereas the electrostatic problem maximizes the harmonic mean.

Literature on the Sphere Problem
I am including a list of literature references. It seems most references concentrate on the electrostatic problem because of its association with molecular chemistry via the classic electron model which is governed by an electrostatic force based on the inverse square law.

Information available on the Internet
The following information on the various sphere problems is available on the internet: Anton Sherwood has a comprehensive overview of the problem of uniform distributions of points on the sphere. Professor Rusin has an FAQ about spheres in his Mathematical Atlas. He also has written a UBASIC for initializing configurations which he makes freely available. Ubasic is a free program for ultra high precision calculations complete with an interpretive BASIC language which was developed by professor Yuji Kida of Rikkyo university in Tokyo, Japan. It has been recommended for number theoretical calculations but I have also found it to be an excellent tool for ultra high precision numerical computations. Neil Sloane provides a list of presumed optimal configurations for the electrostatic problem for order up to 132 and also for some odd cases of higher order. He also lists solutions and data for numerous related problems such as the maximum volume inscribed polyhedron and sphere packings. His interest in the electrostatic problem seems to be limited to finding the global optima and there is no mention of either local optima or meta-stable points. The coordinates of the configurations listed have limited accuracy and the configurations are not normalized. Hugo Pfoertner has models of Sloane's optimum configurations for the electrostatic problem and other interesting material on the sphere problem such as the Optimum Illumination of a Sphere. He also has developed a numerical method of his own. Neubauer, Schilling, Watkins & Zeitlin use a constrained optimization method using Lagrange multipliers to solve the electrostatic problem and report several multiple optima and also mention meta-stable configurations. They also report the asymptotic 4/3 distance law and the asymptotic potential of an infinite covering of electrons. There are pictures of the convex hull for some optimum configurations. Elizabeth D. Dolan and Jorge J. Moré of Argonne National Laboratory use the sphere problem as a benchmark for testing constrained optimization software. Their paper "Benchmarking Optimization Software with COPS" also reports timings. Scott Malloy has a Java applet based animating the optimization which is based on the XYZApp applet for molecule viewing by James Gosling from Sun Microsystems.

Computational Methods Used
The computational methods used to find uniform distributions of points on a sphere have concentrated on the electrostatic problem rather than the geometric problem. Many references do not even mention the numerical algorithm used to solve the problem. This may be partly due to a certain aloofness of numerical methods but more likely so because the performance of the methods used is too embarrassing. The earlier researchers use a heuristic algorithm whereby the points are moved according to the projections of the resultant force onto the tangent plane at the point. Although this appears to be a natural approach for the electrostatic problem it makes use of first derivative information only, like the method of steepest descent. The later references seem to be moving towards more robust and general methods for constrained optimization. There seem to be no numerical results yet on the geometric problem. None of the references mentions the use of methods using second derivatives of the objective function.

Sphere Program
I regard the sphere problem as an interesting mathematical problem in its own right without the need for a physical justification. It is an excellent and natural benchmark problem for testing methods for constrained optimization and auxiliary methods from numerical linear algebra. The optimization problem is not as well conditioned as one would expect because of the presence of optimum configurations with very small eigenvalues for certain orders. The regular polyhedra formed by the various equilibrium configurations are also of considerable interest.

I have developed a constrained minimization program in the C programming language for the computation of stationary configurations of points on the sphere. The objective function is specified as a subroutine and can be selected as geometric, electrostatic or arclength. This program has seen four major versions since 1987. The first version used spherical polar coordinates but later versions used Cartesian coordinates. The constraints consist of the n unit norm constraints augmented by 3 rotation constraints. The rotation constraints state that the inertia matrix of the configuration remain on main axes. It is interesting that these rotation constraints have been ignored by all researchers. The rotation constraints are not essential but their omission leads to degeneracy of the problem which makes it harder to classify stationary points and to a larger solution space than is necessary. In Cartesian coordinates the problem has 3n unknowns and n+3 constraints for a reduced problem space of 2n-3 dimensions. It would be possible to use spherical polar coordinates [provided absolute differentiation were used for the second derivatives] in which case there are only 3 constraints. However, since there are constraints in any case, and also since the purpose of this program is the development of algorithms for constrained optimization as much as it is an investigative tool for collecting regular polyhedra, I have chosen to use Cartesian coordinates. The presence of two types of constraints is also attractive in this context.

New Constrained Optimization Algorithm
The method of solution involves computing the eigenvalues of the reduced second derivative matrix [Hessian] of the objective function. First, the normals to the constraint surfaces are computed as the first derivatives of the constraint functions. Next, the second derivative matrix (Hessian) of the objective function is computed and reduced by certain factors as required by the theory of the constrained geodesic. This is described in my paper on the constrained geodesic in the chapter on the variation of a scalar function along a constrained geodesic. Next, the block reflector based on the normals to the constraint surfaces is computed. The block reflector is a symmetric orthogonal matrix which is the generalization of the Householder transform to a subspace of vectors. The important theory of the block reflector was developed by Schreiber and Parlett partly based on earlier work by other researchers. The block reflector is used to transform the second derivative matrix of the objective function to a system of coordinates where the last n+3 vectors span the subspace of the normals to the constraint surfaces. Next, the principal minor of order 2n-3 of this transformed second derivative matrix is taken and its eigenvalues and eigenvectors in the reduced space [the orthogonal complement of the normals to the constraint surfaces] computed. Finally, the 2n-3 eigenvectors are transformed back from the reduced space to the original space using the (inverse) block reflector. This novel algorithm computes a unique set of 2n-3 mutually orthogonal search vectors in the tangent plane of the manifold [orthogonal to all individual constraint surfaces], which are the principal directions of the constrained problem and whose corresponding eigenvalues are the curvatures.

The selection and use of a search vector depend on its use. In the search for a meta-stable point, it is not clear which one of the search directions it is best to choose. Once a search direction has been chosen the norm of the first derivative [gradient] of the objective function along the constrained geodesic is minimized. In the search for an optimum, a direction having a positive eigenvalue is selected and the maximum of the function along the constrained geodesic is computed.

Automaton
The program has two modes of operation: the automaton and the bootstrap algorithm. The automaton computes meta-stable points or optima of the objective function starting from a random initial configuration by minimizing the norm of its first derivative (gradient) along a search direction. The search direction is chosen from the positive directions of the reduced second derivative matrix of the objective function whose computation has been described earlier. The automaton has proven very effective in finding all meta-stable configurations for order up to about 20. However, beyond order 20 it becomes increasingly ineffective in finding near optimum configurations.

Bootstrap Algorithm
The bootstrap algorithm uses the automaton to find a meta-stable point and follows this with a single maximizing step after which it invokes the automaton again to find the next stationary point. In the maximizing step, one of the positive eigenvalues of the reduced second derivative matrix is picked at random and its associated eigenvector computed. The maximum of the objective function along the constrained geodesic starting from this search direction is found using a line search along the constrained geodesic. The idea behind not always choosing the largest positive eigenvalue is to more exhaustively search the solution space in the optimum region.

Geodesic Line Search
The line search of unconstrained optimization must be replaced by the tracing of the constrained geodesic in the case of constrained optimization. It is possible to use a pseudo-geodesic by moving the individual points along great circles in which case the points stay on the unit sphere but the assembly drifts away from the main axes of the inertia matrix. The other possibility is to integrate the constrained geodesic, for example using the fourth order Taylor series expansion for the sphere problem constraints as presented in the geodesic paper. In this case it is necessary to correct for small departures from the constraints of either type due to the inaccuracies of numerical integration. The unit norm constraints are trivially corrected by re-scaling each point. The rotation constraints may be corrected if desired by bringing the configuration to main axes.

Line Optimization
Since both the objective function and its first derivative with respect to arclength are easily computed along the constrained geodesic, two-point Hermite interpolation may be used to find an improved optimum along a constrained geodesic.

Lists of Configurations
The automaton and the bootstrap algorithm run in an infinite loop and maintain lists of configurations storing only ones that are not yet in the list. It is assumed for this purpose that configurations having equal objective functions are the same. This is of course questionable but distinct configurations having identical objective function have not been found except for stereo isomers which are known to exist but which are ignored here.

Classification of Configurations using the Distances Matrix
It is possible to classify configurations by identifying groups of equivalent points whose set of distances to all other points is the same except for a permutation. This is implemented by comparing the rows of the distance matrix of the configuration to see whether two rows are equivalent save for a permutation. For example, for n=8 a stationary solution is the double hexagonal pyramid consisting - after normalization - of six points forming a hexagon in the equatorial plane and two points at the poles. In this case there are two classes of points, one having six members, and another having two and the resulting configuration is denoted as [6,2]. It should also be possible to identify stereo-isomers this way since these will have the same entries in all rows however permuted such that it is not possible to make them identical by exchanging the points [interchanging two rows and the corresponding two columns].

Normalization of Configurations using the Inertia Matrix
The configurations are normalized by performing a rigid rotation such that the inertia matrix of the rotated system is on main axes. The required rotation matrix has the eigenvectors of the inertia matrix as its columns. The eigenvectors are arranged such that the corresponding eigenvalues are in descending order, which makes the X-axis in the rotated system correspond to the largest eigenvalue, the Y-axis to the second largest, and the Z-axis to the smallest eigenvalue. If there are two or three equal eigenvalues, the momentum ellipsoid degenerates into a spheroid or sphere respectively. In the case of two equal eigenvalues, I associate the corresponding eigenvectors with the X-Y plane which makes the Z-axis an axis of rotational symmetry. This is done in the oblate spheroid case but also in the prolate spheroid case when the two equal eigenvalues are smaller than the remaining one - as is sometimes the case. The additional degree of freedom may be used to rotate an arbitrary point in the equatorial plane into [1,0,0]. I have not yet found an approach to normalize configurations with three equal eigenvalues. If the configuration is not perfect, it may be possible to use the intertia matrix of the largest class of points. Finally, if all the points are in a plane, as with the circle solution, the inertia matrix is defective.

Discussion of Results
Most configurations in this paper are for the geometric problem since I find it the most attractive and also since its resulting configurations seem to be somehwat more regular based on certain criteria which I will discuss.

In Table I, I present what I believe to be a fairly complete list of optima of the geometric function for order [number of points] to 96 together with their objective function values and classification. The table also shows the ranking of each optimum in the list of known configurations and the disposition of the eigenvalues of the inertia matrix. The objective function values are normalized by dividing through by the total number of unique distances. The coordinate values for the configurations are in the individual lists of optimum configurations on the download page. The results are in spherical polar coordinates and the coordinate values are fractions of &pi/2 radians for the latitude and fractions of &pi for the longitude.

I have used the automaton to compute all stable and meta-stable configurations for the geometric objective for order up to 20-25 when the number of solutions becomes prohibitive. I have only a partial list for n=24 and n=25 since in the latter case there are about 29,000 stationary solutions. For order upwards of n=20 it becomes increasingly hard to find the more optimal solutions using the automaton.

I have used the bootstrap algorithm to find most optimum configurations for order to 96. It has proven quite reliable in locating optima and usually requires only 5-8 stationary points to reach a local or global optimum, largely independent of the order. Iteration of the bootstrap algorithm seems to exhaust the optimum solutions in somewhere from 500-2000 tries. The result are very consistent although some optima are much more commonly found than others and may take many tries. The results are very reproducible and it is possible that all optima for order up to 60 have been found and most of the optima for order up to 96.
The bootstrapping algorithm also saves the stationary points visited during successive overall iterations in lists of near optimal solutions. These lists which are not presented here because of their size contain most solutions in the optimum region regardless of classification.

Number of Local Maxima
In contrast to the large number of stationary points which we will discuss next, local maxima are relatively rare as shown in Table II. The lowest order for which a second maximum is known is n=16 where it is the third best solution, next come n=22 and n=32 (both again the third in ranking). From n=32 secondary maxima become more frequent. The largest number of maxima found is 19 for n=92 and n=94. In many cases, secondary maxima are within the most optimum 2-6 objective function values, although for higher order this is not always the case. For higher order starting with 51 there are secondary maxima which are irregular configurations.

The number of maximum solutions found for the geometric objective function for order up to 96 is listed in Table II. Figure I (LARGE!) shows a graph of the global maxima of the geometric objective function for order less than 32. The graph also shows the objective function value of the conjectured minimum stationary solution, the circle solution. In the circle solution the points are arranged at uniform angles within a circle. The optimum objective function for the circle configuration can be computed analytically. The asymptotic value for infinite number of points is &pi/4. The maximum and minimum objective function values (circle) are shown in graph as function of the order in Figure 1. The range of the objective function lies between the lower line which represents the maximum distance for all points on a circle which has a limit of 4/&pi and the upper line for the optimal solution which has a limit of 4/3. The stationary function values are clustered much more densely near the optimum value and much more sparsely near the circle solution.

I have not performed a systematic search for the solutions to the electrostatic problem but have verified that Sloane's solutions in the range from 40-90 are indeed optimum although they are not very accurate and the configurations are not normalized nor are the points classified. In some cases such as n=72 and n=80 the lack of accuracy is due to the presence of very small or even zero eigenvalues in the reduced Hessian which makes finding the optimum solution very difficult.

Number of Meta-Stable Configurations
A remarkable aspect of the sphere problem is the astonishing large number of meta-stable configurations which grows exponentially with the order. The combined number of stable and meta-stable solutions for the geometric objective function is shown in Table III.

I have used the automaton to compile lists of stationary solutions for order up to 25. The numbers presented here are probably accurate for order up to 16 and quite likely higher. The computation of the set of 86 solutions for n=13 requires about 5000 tries and is complete in 3 minutes. The computation of the set of 2000 solutions for n=20 requires about 10 million tries and requires several days of CPU time. The computation of the approximately 29,000 solutions for n=25 requires several months of CPU time and is not yet complete.
The meta-stable configurations having the least optimal objective function value are usually the circular configuration [n points arranged uniformly in a circle], next the pyramid [n-1 points arranged uniformly in a circle and a pole], a non-descript figure [with two points in one hemisphere and the remaining points arranged in some kind of an asteroid belt] and the double pyramid [n-2 points arranged uniformly in a circle in the equatorial plane and two poles]. The circular solution is always the lowest ranking meta-stable solution. It is possible to compute its objective function analytically. The asymptotic difference between the circle solution for a given order and its limit [for n going to infinity] 4/&pi can be analyzed using the Euler-MacLaurin summation formula. The more interesting question is whether a sharp analytical bound can be found for the departure of the optimal solutions from its limit 4/3.

Irregular Configurations
Examining the list of stationary solutions for example for n=24 shows that there are large numbers of solutions lacking any obvious regularity whose distance matrix classification therefore has all points in separate classes. These irregular solutions seem to have no apparent order yet are completely reproducible and there is no doubt as to their genuineness. Nor are they necessarily ill-conditioned like some of the lower ranking solutions where clouds of points are in a kind of asteroid belt and where the optima are weak. The nature of the irregular solutions is yet unknown but it is possible that they possess some kind of quasi regularity which we cannot yet quantify like some kind of a helix. In almost all cases the local and global optimum solution are regular solutions. There is a startling exception for n=61 where the optimum solution for the geometric objective is irregular. For this order, there are three maximum solutions which is not unusual but the globally largest [and also the third largest] are irregular, the second largest being [2,2,2...2]. It seems that n=61 is truly a troubled number since the most regular solution found so far is [3,3,3...3,1].

Perfect Configurations
The number of classes of distinct points in a solution may be considered a measure of its regularity. For certain order there exist perfect solutions all of whose points are within a single class. The convex hull of these perfect solutions does not necessarily have all faces equal. In the cases n=4,6,8,12 the optimum solution is also perfect but for n=16 it is not. Thus - aside from the trivial case of the circle - perfect but non-optimal solutions have been found for various even n. These are for n=6 (1), n=8 (1(cube)) and n=10 (2) n=12 (4), n=14(2), n=16 (2). Perfect solutions often consist of an even number of parallel planes having equal numbers of points in each set of parallel planes and the planes may either be at the same horizontal angles or staggered. The recently discovered crystal form of carbon having sixty atoms [C60] is also a perfect solution but it is not an optimum solution by any means in terms of the objective functions used here.

Center of Gravity
For many configurations the center of gravity [the first moment of the points] is in the origin. There are however amongst the lower ranking solutions ones where this is not the case.

Inertia Matrix Eigenvalues
For some configurations the inertia matrix [the second moment of the points] has three equal eigenvalues and the momentum ellipsoid is a sphere. Such solutions have been found for n=4, (the regular tetrahedron), for n=8 (the cube) but not the square anti-prism, for n=16, for n=24 for the optimum solution (both stereo-isomers) but also for a second non-trivial solution and finally for n=60 for the soccer ball.

Comparison of Solutions to the Geometric and Electrostatic Problems
The ordered lists of solutions to these two problems from the global optimum downwards to the circle solution are pairwise very similar with pairs of corresponding solutions having the same general shape and classification however with missing entries in either list. The number of solutions for the geometric and the electrostatic problem is approximately equal with the number of solutions for the electrostatic problem sometimes marginally higher. For order less than 10 the number of solutions is the same for both objective functions but from n=10 onwards the ranking of a few similar solutions changes positions in the lists and in some cases an extra solution exists for one problem and not for the other. For example for n=10 we have (28,29) solutions respectively and for n=12 (58,57) solutions. The number of optimum solutions for the geometric objective function is tabulated in Table II but I have not generated the equivalent table for the electrostatic objective function. In most cases, similar shapes are found in each list at almost the same ranking order although the numerical values of the corresponding coordinates may differ appreciably. In a much smaller number of cases, the solutions to both problems are the same. The similarity of corresponding solutions to both problems is usually such that a solution for one problem converges very rapidly to the corresponding solution for the other problem if there is one.

For certain number of points the geometric problem has as its optimum solution a more regular figure than the electrostatic problem: the optimum solutions for both n=66 and n=78 have classification [6,6,6...6] for the geometric problem but [2,2,2...2] for the electrostatic problem. Finally, the range of the eigenvalues of the inertia matrix for similar solutions for the three objective functions seems to be smallest for the geometric objective function, larger for the electrostatic objective function and largest for the arclength objective function. This suggests that the solutions to the geometric problem are the most spherical. I have collected some more results for reference on a separate page for those who are interested.

Krystal Program
The Krystal program is a Visual Basic application designed for viewing lists of stationary configurations. The configurations may be viewed from a user defined viewpoint or they may be animated by rotating them around the vertical axis. The viewpoint may be changed using latitude and longitude slider controls. The configurations are retrieved from a Microsoft Access database using SQL calls. An auxiliary program is required to populate the database from the flat files containing the configurations. The configurations are displayed by tracing the polygons spanning their convex hull which are computed using the Qhull program. QHull is a program available from the Geometry Center of the University of Minnesota which computes the convex hull of a set of points [in an arbitrary number of dimensions] and which has many other related capabilities. It is written in the C language and has been compiled into a DLL to interface with the Krystal program which is written in Visual Basic. Only the convex hull functionality of the program is used although there is an option to calculate its volume. The configurations are retrieved one at a time but the program can move forward and backward with considerable speed. The classification and stability information are also displayed. There is currently no provision to edit configurations or write them back to the database. It is interesting to examine the lists of solutions for larger order using the Krystal program. The better configurations with the higher objective function values look like regular polyhedra (or less regular for the irregular configurations). The lower objective function values become increasingly more rambling the lower the objective function value, having shapes reminiscent of circus tents with many guy wires or intricately curved saddles much like the saddle dome, a certain stadium in my former home town Calgary. There are many such configurations consisting of three or four lobes where the points seem to be arranged in curved lines like a giant tent. The variety in these configurations is truly amazing. The lowest configuration is always the circle for which it is possible to give an analytical expression for the stationary solution. The third worst solution is usually the single pyramid with one point at the north pole and n-1 points arranged in a circle at a latitude below the equator. There is usually an interesting solution in between the circle and the single pyramid which consists of two points at equal latitude in the northern hemisphere and the remaining n-2 points arranged in some kind of an asteroid belt below the equator. This solution is invariably very indeterminate and hard to converge, indicating that its stability almost extends beyond a single point. The double pyramid is always a solution which is slightly more optimal than the single pyramid.

Initial Configuration
The initial configurations are generated by using a random generator to generate points within a cube of base 2 and projecting points inside the unit sphere onto its surface while discarding the points outside the unit sphere. I use a high quality random generator R250 by Kirkpatrick and Stoll, an implementation of which may be found in Berndt M. Gammel's Matpack library. The random generator is also used by the bootstrapping algorithm to choose from the available positive directions. The quality of the random generator is expected to have a positive influence on the results but this has not been investigated.

Bootstrapping Algorithm Efficiency
The bootstrapping algorithm requires typically 6 but almost invariably between 5-8 stationary points to move from a random initial solution [counted as the first] to an optimum solution. Typical performance for the bootstrapping algorithm on a 350 MHz Pentium II for n=60 is about 600 optimal solutions per day, 37% of which are the global maximum, 53% the secondary maximum and 10% the third local maximum. Like the automaton, the bootstrapping algorithm saves intermediate stationary solutions in a file. This file is read in again when the program is restarted and over time becomes a repository of solutions in the optimum region with some lesser solutions appended which may be discarded. I usually retain all solutions whose objective function exceeds the lowest optimum found and any additional regular configurations having one or more classes with three or more points. The random selection of an eigenvector from the ones having positive eigenvalues also causes the solution space in the optimum region to be exhaustively searched. The algorithm reliably generates the meta-stable solutions with the highest objective function values between runs. For example, for order 24, the total number of solutions is about 20,000 but when using the Automaton, members of the optimum 1000 solutions are found only very rarely, maybe a few out of 100,000 tries. However, using the bootstrapping algorithm 10,000 times to find a maximum solution while cataloguing the intermediate solutions, the resulting set of about 60,000 solutions contains most of the 1000 most optimal solutions.

Discussion
The constrained optimization algorithm presented here whereby the block reflector is used to transform the reduced Hessian matrix to the orthogonal complement of the normals to the constraint surfaces and to transform the computed eigenvectors back is entirely new. The use of the reduced Hessian as specified by the theory of the constrained geodesic is entirely new.
The bootstrapping algorithm where the meta-stable points of the objective function are used as stepping stones in a deterministic search for optima is a new idea.

Several of the ideas used in the sphere problem such as the classification of configurations using the distance matrix, the normalization of configurations using the inertia matrix and the normalization of the objective function are also new.

There are still lots of open questions: The most intriguing is whether continuous families of optimum configurations exist and their discovery would be very exciting. It is now known that there are many optimum configurations which seem to be surrounded by a region which is nearly optimum (one or more very small eigenvalues). There is also the issue of the nature of irregular configurations. Finally, there is the question of the asymptotic behaviour of the optimum values of the objective function.

There are also unresolved algorithmic issues such as whether to use the norm of the gradient of the objective function as a new objective function. This is quite complicated and has not yet been tried.