Testing the "No Free Lunch Theorem"

Circle and Sphere Packings with MathOptimizer Professional

Frank Kampas
WAM Systems, Plymouth Meeting, PA, USA
fkampas@msn.com

János D. Pintér
Pintér Consulting Services Inc. Halifax, NS, Canada
jdpinter@hfx.eastlink.ca
www.pinterconsulting.com

Introduction

In computational optimization, there are circumstances in which the outputs of all procedures solving a particular type of problem are statistically identical. In other words, for any algorithm, any elevated performance over one class of problems is offset by inferior performance over another class (when the statement holds).  

A colorful way of describing such a circumstance, introduced by Wolpert and Macready (1995, 1997), is to say that there is "no free lunch" (NFL) in search and optimization. For a detailed discussion, including downloadable documents, see http://en.wikipedia.org/wiki/No_free_lunch_in_search_and_optimization.

We test the general NFL idea by solving non-trivial instances of circle and sphere packing problems. Specifically, we will seek high-quality (in theory, globally optimized) solutions to the problem of finding a minimal size circle (sphere) that contains a given collection of non-uniform size circles (spheres) in a non-overlapping arrangement. 

First, we perform an initial generic global optimization procedure, and then we use various local search techniques to improve the initial solution.  These local search techniques involve swapping the position of the circles or spheres in order to improve the objective function, the radius of an enclosing circle or sphere.  The different heuristic search methods studied include swapping all pairs of adjacent size circles or spheres, swapping all pairs, swapping adjacent size pairs starting from the smallest pair, swapping all pairs starting from the smallest pairs and performing the previously described swapping techniques either once or repeatedly until the objective function stops improving.  

The various techniques are performed multiple times and finally the best objective function result is plotted vs. time as well as the distribution of results.  The new graphics features in Mathematica version 6, especially for three dimensional plots, are employed to display our results.

The initial global optimization and all local searches are performed using our MathOptimizer Professional third party software package (Pintér and Kampas, 2003), taking advantage of its ability to reuse machine compiled code for model function evaluations. The reuse of machine compiled code enables the circle or sphere swapping to be performed over one thousand times faster than the equivalent calculation performed in Mathematica using FindMinimum.  

MathOptimizer Professional has been described in other works that review its theoretical foundations, key implementation details, and the solution of test problems: the latter also include the solution of various circle packing problems and their potential industrial applications. For additional details not discussed here, consult Kampas and Pintér (2004, 2006), Pintér and Kampas (2005, 2006), Castillo, Kampas and Pintér (2007).

The Basic Swap Step for Circles

Swap the two circles with radii r1 and r2.  The circle radii r, the variables vars, the objective function obj, and the constraints cons could be calculated from the starting solution soln, but are passed as arguments to save time.  Note that the optimizer result is enclosed in Sow[].

swapcircles[r_,vars_,obj_,cons_,soln_,r1_,r2_,max_]:=
Module[
{init},
init = vars /. (soln[[2]]/. {x[r1] -> x[r2], x[r2] -> x[r1],y[r1] -> y[r2], y[r2] -> y[r1]});
Sow[callLGO[obj,cons,vwbLS[r,init,max],Method -> LS,Recompile -> False,TimeStamp -> True]]
]

We shall illustrate the circle swapping function by swapping the two largest circles: we use the global search based solution from the previous section.

slnSwapTwo=AbsoluteTiming[swapcircles[radii,varlist,objtest,constest,slnGlobal[[2]],1,(1/(Sqrt[2])),2.5]]

The plot shows that circles 1 and 2 have exchanged positions, as expected: of course, this exchange could (and typically will) affect also the position of the other circles. As the results show, the objective function (circumscribing radius) 2.00637 became inferior to the previously found value 1.95496. This finding, however, does not need to be valid in general, especially when allowing for more general swapping strategies to be discussed later on.

Show @ GraphicsArray[{p1,showcircles[radii,slnSwapTwo[[2]]]}]

Let us remark that a similar function, implemented with FindMinimum, takes over a thousand times longer and returns a poorer result.

swapcirclesFM[r_,vars_,obj_,cons_,soln_,r1_,r2_,max_]:=
Module[
{init},
init = vars /. (soln[[2]]/. {x[r1] -> x[r2], x[r2] -> x[r1],y[r1] -> y[r2], y[r2] -> y[r1]});
FindMinimum[Join[{obj},cons],vwbLS[r,init,max][[All,{1,3}]]]
]

AbsoluteTiming[swapcirclesFM[radii,varlist,objtest,constest,slnGlobal[[2]],1,(1/(Sqrt[2])),2.5]]

The Basic Swap Step for Spheres

Swap the two spheres with radii r1 and r2.  The sphere radii r, the variables vars, the objective function obj, and the constraints cons could be calculated from the starting solution soln but are passed as arguments to save time.  Note that the optimizer result is enclosed in Sow[].

swapspheres[r_,vars_,obj_,cons_,soln_,r1_,r2_,max_]:=
Module[
{init},
init = vars /. (soln[[2]]/. {x[r1] -> x[r2], x[r2] -> x[r1],y[r1] -> y[r2], y[r2] -> y[r1],z[r1] -> z[r2], z[r2] -> z[r1]});
Sow[callLGO[obj,cons,vwbLS[r,init,max],Method -> LS,Recompile -> False,TimeStamp -> True]]
]

Illustrate the sphere swapping function by swapping the two largest spheres

slnSwapTwo=AbsoluteTiming[swapspheres[radii,varlist,objtest,constest,slnGlobal[[2]],1,(1/(2^(1/3))),2.5]]

The plot shows that spheres 1 and 2 have exchanged positions, as expected.

Show @ GraphicsArray[{p1,showspheres[radii,slnSwapTwo[[2]]]}]

A similar function, implemented with FindMinimum (below) takes a few thousand times longer, and often returns a poorer result than MOP.

swapspheresFM[r_,vars_,obj_,cons_,soln_,r1_,r2_,max_]:=
Module[
{init},
init = vars /. (soln[[2]]/. {x[r1] -> x[r2], x[r2] -> x[r1],y[r1] -> y[r2], y[r2] -> y[r1],z[r1] -> z[r2], z[r2] -> z[r1]});
FindMinimum[Join[{obj},cons],vwbLS[r,init,max][[All,{1,3}]]]
]

AbsoluteTiming[swapspheresFM[radii,varlist,objtest,constest,slnGlobal[[2]],1,(1/(2^(1/3))),2.5]]

Search Strategies Included in the Comparative Study

Global:  Global search (MathOptimizer Professional used in its default search mode) 

Random Search:  Local search starting from randomly chosen starting values (MathOptimizer Professional used in its local search mode) 

SwapAllAdj:  Swap all adjacent size circles (spheres), and keep the best result (i.e., the configuration that needs the smallest enclosing circle or sphere):
Best[ 10⇌9, 9⇌8, 8⇌7, 7⇌6, 6⇌5, 5⇌4, 4⇌3, 3⇌2, 2⇌1 ]

SwapAllAdjRpt:  Repeat the previous strategy as long as the objective function value keeps improving

SeqSwapAdj:  Swap adjacent size circles (spheres), starting with the smallest pair, move up to larger pairs, and return the first result that is better than the input:  
First Improved Result[ 10⇌9, 9⇌8, 8⇌7, 7⇌6, 6⇌5, 5⇌4, 4⇌3, 3⇌2, 2⇌1 ]

SeqSwapAdjRpt:  Repeat the previous strategy, but keep going as long as the objective function value keeps improving

SwapAll:  Swap all pairs of circles (spheres) and return the best result:
Best[ 10⇌9 ... 2⇌1, 10⇌8 ... 3⇌1, 10⇌7 ... 4⇌1, 10⇌6 ... 5⇌1, 10⇌5 ... 6⇌1, ... , ... ]

SwapAllRpt:  Repeat the previous strategy as long as the objective function value keeps improving

SeqSwapAll:  Swap all pairs of circles (spheres) starting with the smallest pair, move up to larger pairs and return the first result that is better than the input:
First Improved Result[ 10⇌9 ... 2⇌1, 10⇌8 ... 3⇌1, 10⇌7 ... 4⇌1, 10⇌6 ... 5⇌1, 10⇌5 ... 6⇌1, ... , ... ]

SeqSwapAllRpt:  Repeat the previous strategy as long as the objective function value keeps improving

Summary of Circle Packing Results

The best result :

This section compares and sums up our numerical experiments related to the 10-circle packing example. We collect the results from all previous sections, and present the best results in table format. This is followed by graphical summaries.

The next table displays the best objective function value (radius of circumscribing circle) by solution method.  The best results were obtained by SwapAllRpt, SeqSwapAdjRpt and SeqSwapAllRpt.  For the record, we print the results here with higher precision than earlier.

({{"Global", 1.9223094120209718}, {"RandomSearch", 1.9240487093536316}, {"SwapAllAdj", 1.9274601565517050}, {"SwapAllAdjRpt", 1.9174260008224506}, {"SeqSwapAdj", 1.9221391373214323}, {"SeqSwapAdjRpt", 1.9158205231840948}, {"SwapAll", 1.9263515495346892}, {"SwapAllRpt", 1.9148313247862758}, {"SeqSwapAll", 1.9221391373214323}, {\, 1.9148293930115357}})

The following plot displays the best objective function value by method and elapsed time.  Note that intermediate results are also shown: in other words, the best results from intermediate circle swappings are also displayed, together with the final results for a technique.  The number of trials for each method was adjusted so that the total time required was approximately 200 seconds.  Average times for each method were

Global 			0.44 seconds
Random Search	0.20 seconds
SwapAllAdj		2.01 seconds
SwapAllAdjRpt	8.52 seconds
SeqSwapAdj		0.76 seconds
SeqSwapAdjRpt	5.11 seconds
SwapAll		9.50 seconds
SwapAllRpt	          83.5  seconds
SeqSwapAll		 0.78 seconds
SeqSwapAllRpt          16.8 seconds

Key to the two plots:

Global Optimization With Different Random Seeds
Local Optimization Starting From Random Points
Sequential Adjacent Circle Swapping
Repeated Sequential Adjacent Circle Swapping
Adjacent Circle Swapping
Repeated Adjacent Circle Swapping
Swapping All Pairs of Circles
Repeated Swapping of All Pairs of Circles
Sequential Swapping of All Circles
Repeated Sequential Swapping of All Circles

The second plot shows the cumulative objective function distributions by method.  This plot shows that the best three methods are the three repeated methods with swapping all circles, superior to swapping all adjacent size circles, which is slightly superior to sequential swapping of adjacent size circles.

Repeated sequential swapping of all pairs of circles, moving on with the first improvement of the objective function, appears to be the best technique, achieving the best value found for the objective function in the shortest time.  The overall results seem to indicate that the global optimum should be close to 1.91: accepting this postulated value, all strategies seem to be able to generate at their best a solution which is within ~1.5% of the putative optimum. If the optimum is in fact close to the best value found ~1.91483, then the best results with all strategies are all within 0.66% of the putative optimum value.  Each of these results were generated in a few seconds or minutes.

Summary of Sphere Packing Results

The best result :

This section compares and sums up our numerical experiments related to the 10-sphere packing example. We collect the results from all previous sections, and present the best ones in a table; this is followed by graphical summaries.

The next table shows the best objective function value (radius of enclosing sphere) obtained by each strategy used.  The best results were obtained by SwapAllRpt, SeqSwapAdjRpt and SwapAllAdjRpt, the last one being the best.  Again, for the record, we write out the results here with higher precision than earlier.

({{"Global", 1.8513359075427580}, {"RandomSearch", 1.8487216814810763}, {"SwapAllAdj", 1.8382279699258777}, {"SwapAllAdjRpt", 1.8259562513496725}, {"SeqSwapAdj", 1.8417251648541804}, {"SeqSwapAdjRpt", 1.8412968938489418}, {"SwapAll", 1.8468200654372806}, {"SwapAllRpt", 1.8301947846570883}, {"SeqSwapAll", 1.8417251648541804}, {"SeqSwapAllRpt", 1.8330810637200052}})

The overall quality of the best results is fairly similar: the ratio of the worst and best result is ~ 1.013899378; it.e. the relative difference between these is less than 1.4 %.

The number of trials in the runs for each method was adjusted so that the total time was approximately 200 seconds.  The average time for each trial for the different methods was as follows:

Global			0.46 seconds
Random Search	0.24 seconds
SwapAllAdj		2.29 seconds
SwapAllAdjRpt          13.9  seconds
SeqSwapAdj		0.78 seconds
SeqSwapAdjRpt	7.12 seconds
SwapAll	          11.7 seconds
SwapAllRpt		69  seconds
SeqSwapAll		 0.79 seconds
SeqSwapAllRpt	21 seconds

We show next the best objective function value by method and elapsed time.  Again, intermediate results are also shown: in other words, the best results from intermediate sphere swappings are shown, not just the final results for a technique.  The best result at longer times is obtained with repeated swapping of all adjacent size spheres.

Key to the two plots (above and below):

Global Optimization With Different Random Seeds
Local Optimization Starting From Random Points
Sequential Adjacent Sphere Swapping
Repeated Sequential Adjacent Sphere Swapping
Adjacent Sphere Swapping
Repeated Adjacent Sphere Swapping
Swapping All Pairs of Spheres
Repeated Swapping of All Pairs of Spheres
Sequential Swapping of All Spheres
Repeated Sequential Swapping of All Spheres

Finally, we show the cumulative objective function value distributions, for each method.  Repeated Adjacent Size Sphere Swapping gives the best result.

Conclusions

In this article, we present a fully scalable model development procedure and propose strategies to solve, in principle, arbitrary size instances of circle (and sphere) packing problems.  For illustration, we use a 10-circle and another 10-sphere model instance. These models are solved (approximately) by using various solution approaches, with the objective of comparing these approaches, as well as of finding high-quality global optimum estimates.  (As far as we know, nobody has published the proven global solution to these problems.)

For both models, we first perform an initial generic global optimization procedure, and then use various relatively quick heuristic strategies combined with a local search procedure to improve the initial solution.  These local search techniques are all based on swapping the position of several circles (spheres) and then repeating the local search, in order to improve the solution estimate.  The initial global optimization and all subsequent local searches are performed using our MathOptimizer Professional third party software package.  For both circles and spheres, the various techniques are performed multiple times and finally the best results are plotted vs. time, together with a plot of the empirical distribution of the optimum estimates.

In our extensive, but admittedly still limited numerical experience, our study in itself does not support the unqualified use of the "No Free Lunch" meta-theorem. Specifically, our results show that more thorough swapping strategies rather consistently produce somewhat better results; at the same time, these strategies also take more time.  This would give a different meaning to the NFL paradigm.

The overall results also indicate that the solution strategies proposed are able to produce solutions that are fairly close to the putative global optimum. Our numerical results can also be used to apply statistical techniques to estimate the global optimum value. Another deterministic approach could be to use verified optimum estimates based e.g. on using interval arithmetic.

Finally, let us remark that larger model instances (with more circles or spheres to pack) can lead to much more difficult models: the number of decision variables grows only linearly,  but the number of non-convex constraints grows (approximately) as the square of the number of circles. The current shipment version of MathOptimizer Professional enables users to solve nonlinear optimization models with up to a few thousand variables and constraints. In other numerical tests not reported here, we have generated global search based high-quality feasible solutions to model instances with up to 60 circles in a few minutes, using average desktop computers. The above facts show the overall numerical viability of the generic approach based on global optimization, the heuristics suggested, and the use of our MOP package to handle these and similar problems. For additional numerical results and further details, consult e.g. Kampas and Pintér (2004, 2006), Pintér and Kampas (2005, 2006), Castillo, Kampas and Pintér (2007).

References

Castillo, I., Kampas, F.J. and Pintér, J.D. (2007) Solving Circle Packing Problems by Global Optimization: Numerical Results and Industrial Applications. European Journal of Operations Research (To appear.)

Horst, R., and Pardalos, P.M., Eds. (1995) Handbook of Global Optimization. Vol. 1. Kluwer Academic Publishers, Dordrecht.

Kampas, F.J. and Pintér, J.D. (2004) Generalized Circle Packings: Model Formulations and Numerical Results. In: P. Mitic, C. J. Jacob, and J. Carne, Eds. New Ideas in Symbolic Computation: Proceedings of the 6th International Mathematica Symposium. Positive Corporation, Hampshire, 2004.

Kampas, F.J. and Pintér, J.D. (2006) Configuration Analysis and Design by Using Optimization Tools in Mathematica. The Mathematica Journal 10 (1), 128-154.

Pardalos, P.M., and Romeijn, H.E., Eds. (2002) Handbook of Global Optimization. Vol. 2. Kluwer Academic Publishers, Dordrecht.

Pintér, J.D (1996) Global Optimization in Action – Continuous and Lipschitz Optimization: Algorithms, Implementations and Applications. Kluwer Academic Publishers, Dordrecht

Pintér, J.D (2002) Global Optimization: Software, Test Problems, and Applications. In: Pardalos, P. M. and Romeijn, H. E., Eds. Handbook of Global Optimization, Volume 2, pp. 515-569. Kluwer Academic Publishers, Dordrecht.

Pintér, J. D.  (2005) Nonlinear Optimization in Modeling Environments: Software Implementations for Compilers, Spreadsheets, Modeling Languages, and Integrated Computing Systems. In: V. Jeyakumar and A. M. Rubinov, Eds., Continuous Optimization: Current Trends and Applications, pp. 147-173. Springer Science + Business Media, New York.

Pintér, J.D. and Kampas, F.J. (2003) MathOptimizer Professional: An Advanced Modeling and Optimization System for Mathematica, Using the LGO Solver Engine. User Guide. Published and distributed by Pintér Consulting Services, Inc., Halifax, NS. For summary information and other references, see also http://www.wolfram.com/products/applications/mathoptpro/.

Pintér, J.D. and Kampas, F.J. (2005) Nonlinear Optimization in Mathematica with MathOptimizer Professional. Mathematica in Education and Research 10 (2), 1-18.

Pintér, J.D. and Kampas, F.J. (2006) MathOptimizer Professional: Key Features and Illustrative Applications. In: Liberti, L., and Maculan, N., Eds. Global Optimization: From Theory to Implementation, pp. 263-279. Springer Science + Business Media, New York, 2006.

Wolpert, D.H., Macready, W.G. (1995), No Free Lunch Theorems for Search, Technical Report SFI-TR-95-02-010, Santa Fe Institute.

Wolpert, D.H., Macready, W.G. (1997), No Free Lunch Theorems for Optimization, IEEE Transactions on Evolutionary Computation 1, 67. http://ic.arc.nasa.gov/people/dhw/papers/78.pdf