Global Optimization in Mathematica:  A Comparative Numerical Study

Frank Kampas

WAM Systems, Plymouth Meeting, PA
fkampas@wamsystem.com

János D. Pintér

Pintér Consuling Services, Halifax, NS
jdpinter@hfx.eastlink.ca
www.pinterconsulting.com

Introduction

In this presentation, the built-in global optimization function NMinimize is compared with two third-party tools, Global Optimization (developed by Loehle Enterprises) and MathOptimizer Professional (developed by Pintér Consulting Services and the author).  The comparison was performed for both 71 small problems (1 to 4 variables, 0 to 4 constraints) and for selected larger problems.  We looked at the accuracy of the results and the time required to generate them.  The scope of the comparison was limited to problems with continuous variables only.

Optimizer and Comparison Details

In global optimization we wish to minimize a function f(x) subject to equality constraints g(x) == 0 and inequality constraints which can always be written as h(x) <= 0.  The possible existence of multiple local minima makes the problem a hard one, since a local search from a starting point will only find a local minimum, not is probably not the global minimum. 

The comparison was made between: 

three NMinimize options:  NelderMead, DifferentialEvolution and SimulatedAnnealing, 
two Global Optimization Functions:  GlobalSearch and GlobalPenaltyFn, 
and four MathOptimizer Professional options:  MultiStart Global Search, Global Adaptive Random Search, Branch and Bound and Local Search (for some problems without multiple minima).  

Each optimization function or option was run using its default settings. 

NelderMead examines the values of the objective function at the vertices of a sampling-based polyhedron and then modifies the polyhedron in a way which depends upon the relative values at the vertices, moving the polyhedron through the search space.  DifferentialEvolution constructs a population of solutions which are evolved by "mating" and selection.  SimulatedAnnealing performs randomized trial steps in a population and always accepts those changes for improvements in the objective function and may accept changes which make the objective function worse, depending upon an adaptively chosen Boltzmann factor for the change.  All three methods use a penalty function to handle constraints.

GlobalSearch is described as a multiple-start generalized hill-climbing algorithm that is based on Newton's method but uses a generalized gradient and allows for constraints, which must be analytic.  GlobalPenaltyFn has similar capabilities but can handle the situation in which the constraints cannot be solved for any of the variables or are black box, by using an adaptive penalty function. 

MathOptimizer Professional is a Mathematica interface to Pintér Consulting Services' LGO (Lipschitz Global Optimizer) solver suite, which is available to use in conjunction with various modeling platforms: compilers, spreadhseets, modeling languages, and technical computing systems.  MathOptimizer Professional converts the Mathematica model formulation to C or Fortran (using Mark Sofroniou's Format package) and then compiles the problem using an external compiler.  The compiled code is either called by an LGO executable (as a Windows dll) or linked to an LGO object file to create an executable.  The default solver method, MultiStart Global Search (MSLS) uses a multi-start adaptive global random search.  The Global Adaptive Random Search (GARSLS) method is based on a single-start random search. The Branch and Bound (BBLS) method is an adaptive partition and search method.  All three global search options are automatically followed by a constrained local search (generalized reduced gradient) method (denoted by LS).  The local search method can also be used without a global search phase.  Unlike the other optimizers studied, MathOptimizer Professional only looks inside the region specified by the variable bounds, which it divides and redivides in a systematic way while searching for the global minimum.            

All of the optimizers studied are at least partly stochastic in nature, generating random starting points and/or using random search in the course of the optimization.  Global optimization tools which are rigorous (with deteministic global convergence) or analytical, are – as a rule – prohibitively slow for problems with more than a few variables and constraints.  As a consequence of the embedded stochasticity, multiple runs with different random seeds can sometimes improve the numerical results obtained.  For the 71 small problems, one run was used for each function, except for the Global Optimization functions, for which the default number of 3 was used, unless more were required to obtain a known global minimum.   

Mathematica Version 5.2 was used running under Windows XP on a Dell with a 3GHz Pentium processor.  The Global Optimization version number used was 5.1 (the most recent version).  The external compiler used by MathOptimizer Professional was the Intel Fortran Compiler (Version 9).

Results for Small Problems

For the 71 small test problems, the largest number of correct results were found by Global Optimization Global Search (68 correct results), NMinimize Differential Evolution (66 correct results) and MathOptimizer Professional MSLS (65 correct results).  The next runners-up were MathOptimizer Professional GARSLS (61 correct results) and Global Optimization GeneralPenaltyFn (61 correct results).  The smallest number of correct results were found by MathOptimizer Professional BBLS (54 correct results), NMinimize SimulatedAnnealing (53 correct results) and NMinimize NelderMead (50 correct results).  Note that numerical correctness is defined as a precision of (10^(-6))for this study.

For the most difficult problem (the Trefethen problem from the Hundred Dollar, Hundred Digit Challenge) none of the methods found the global minimum using default settings, although 2 of the methods (MSLS and DifferentialEvolution) found the correct result when the number of trials was tripled for each method, as shown below.  The relative order of the elapsed times in this example is typical of the relative order of times required for the small problems.  The larger times required for the MathOptimizer Professional functions are due to the overhead of external compilation of the problem.  The Global Optimization functions typically were among the fastest.  

The syntax of our function compareOpts is compareOpts[ objective function, {constraints}, {variables and their bounds}, {number of runs for each optimizer function} ].

plt @ 
      compareOpts[
        Exp[Sin[50*x]] + Sin[60*Exp[y]] + Sin[70*Sin[x]] + Sin[Sin[80*y]] - 
          Sin[10*(x + y)] + (((x^2) + (y^2))/4), {}, {{x, -3, 
            3}, {y, -3, 3}}, {3, 3, 3, 3, 3, 3, 9, 9}]

Timing Comparison

A simple unconstrained minimization of ∑((Subscript[x, i])^2) for 1, 3, 10, 30, 300, and 1000 variables was performed in order to assess the relative speed of the different optimizers over a wide range of variable number.  The results are shown in the following graph.  Since the problem is a local optimization problem, the MathOptimizer Professional method LS (Local Search) was included.

The four MathOptimizer Professional methods, MSLS, GARSLS and BBLS and LS all have an overhead of around 0.4 seconds due to the external compilation step.  As a result the time required by them remains relatively flat up to 30 or more variables.  Note that the LGO execution time for the LS method for 1000 variables was around 0.1 seconds,  compared to the total time of 1.72 seconds, indicating that compilation of the objective function took longer than the optimization itself, for this problem.

The two Global Optimization methods (GlobalSearch and GlobalPenaltyFn) are the fastest, except for the MathOptimizer Professional BBLS  and LS which are faster when there are several hundred variables or more.  NMinimize NelderMead is faster than the MathOptimizer Professional methods up to about 100 variables and the same speed for more than 100 variables.  NMinimize SimulatedAnnealing and DifferentialEvolution are the slowest for more than 10 or so variables.  Of the three methods that gave the most accurate results on the small test problems, Global Optimization GlobalSearch is the fastest, although the times appear to be converging at 1000 variables.

Hanging Chain Problem

The hanging chain problem involves find the configuration of a chain of rigid links supported at each end.  The objective function to be minimized is the total potential energy of the links.  There are two kinds of constraints.  The boundary condition constraints set the positions of the left and right end of the chain.  The other constraints, one for each link, state that the distance between the ends of a link is less than or equal to the length of the link.  Use of an inequality constraints rather than equality constraints makes the problem easier to solve.

Assume that there are n links in the chain, each with length h and with a total length of 2.  The coordinates of the left and right end of link i are { x[ i - 1 ], y[ i - 1 ] } and { x[ i ], y[ i ] }.  The left end of first link is at { 0, 0 } and the right end of the last link is at  { 1, 0 }.  The statement of the problem is

h := (2/n)

xvars := Table[x[i], {i, 0, n}]      
yvars:=Table[y[i], {i, 0, n}]

allvars := Join[xvars, yvars]

obj :=  (∑_(i = 1)^n)((y[i - 1] + y[i])/2)

bcs := {x[0] == 0, y[0] == 0, x[n] == 1, 
          y[n] == 0};

lincons := 
        Table[(((x[j] - x[j - 1]))^2) + (((y[j] - y[j - 1]))^2) ≤   (h^2), {j, 1, n}];

allcons:=Join[bcs, lincons]

vwb := 
        Transpose[{allvars, Table[-1, {2  n + 2}], 
            Table[1, {2  n + 2}]}];

prob := {obj, allcons, vwb} /.  
        x_[i_] :>   ToExpression[
            ToString[x] <> ToString[i]]

For 10 links we obtain the following result.  We see that the best results were obtained by the three NMinimize methods (NelderMead, DifferentialEvolution, SimulatedAnnealing), with Global Optimization GlobalPenaltyFn  and MathOptimizerProfessional MSLS also producing good results.

n  =  10;
plt @ 
      compareOpts[Sequence  @@  prob, {1, 1, 1, 1, 1, 1, 3, 3}]

For 20 links, the results are similar, with the result found by MathOptimizer Professional GARSLS fairly close to the results found by the NMinimize methods.

n :=  20

plt @ 
      compareOpts[Sequence  @@  prob, {1, 1, 1, 1, 1, 1, 3, 3}]

Note that a close approximation to the global minimum for this problem can be found by MathOptimizer Professional, using a local search, implying that there is only one minimum.

AbsoluteTiming @ 
      callLGO[Sequence  @@  prob,  Method  ->  LS]

Lennard-Jones Cluster

The Lennard-Jones potential is a model potential for the interaction between two atoms and is proportional to (r^(-12))- (r^(-6)).  To find the configuration of a cluster of "Lennard-Jonesium" one must minimize the potential energy of all the pairs of atoms, so that the objective function has n( n - 1 )/2 terms.  There are no constraints, but the model is highly multimodal.  Known results can be found at http://www-doye.ch.cam.ac.uk/jon/structures/LJ/tables.150.html

The next figure displays the Lennard-Jones potential.

Plot[(1/(r^12)) - (1/(r^6)), {r, 1, 3}]

p[i_]:={x[i], y[i], z[i]}

d2[i_, j_]:=(p[i] - p[j]).(p[i] - p[j]) (* 
            the square of the distance between atom i and atom j *)

v[n_]:=4 Sum[
          d2[i, j]^(-6) - d2[i, j]^(-3), {i, n - 1}, {j, i + 1, n}]

vars[n_]:=Flatten @ Table[p[i], {i, n}]

varswithbounds[n_, b_]:=Transpose[{vars[n], Table[-b + i/100, {i, 3  n}], 
          Table[b + i/100, {i, 3  n}]}]

prob[n_, b_]:={v[n], {}, varswithbounds[n, b]} /.  x_[i_]:>  ToExpression[ToString[x] <> ToString[i]]

All methods found the correct result ( -9.103852 ) for 5 atoms (15 decision variables).  The two Global Optimization Methods (GlobalSearch and GlobalPenaltyFn) were the fastest.

plt @ 
      compareOpts[
        Sequence  @@  prob[5, 1], {1, 1, 1, 1, 1, 1, 3, 3}]

For 20 atoms in the cluster, only MathOptimizer Professional MSLS found the correct result ( -77.177043 ).

plt @ 
      compareOpts[
        Sequence  @@  prob[20, 2], {1, 1, 1, 1, 1, 1, 3, 3}]

If the number of repeated runs is adjusted to give (approximately) equal runtimes to all solvers, then MathOptimizer Professional BBLS also finds the correct result. The GARSLS result also improves in multiple runs, and all MOP options find better solutions than the other solvers; see figure below.

plt@
      compareOpts[
        Sequence @@ prob[20, 2], {1, 12, 16, 2, 1, 1, 5, 6}]

Packing Circles

We wish to pack a given collection of n arbitrary (possibly different) size circles into a circumscribing circle, in a non-overlapping arrangement. The objective is to minimize the radius of the circumscribing circle; see an illustrative arrangement below for 25 circles. Note that this model has n(n-1)/2 non-convex constraints, in addtion to n "easy" convex constraints. To our best knowledge, there are no previosly known results for this general model type, except in our own recent numerical studies.

In the model formulation shown below, circumscribecircles is a function to generate the objective function, constraints and variables with bounds for n circles, where the size of the ith circle can be defined by func[i].  There is only one term in the objective function, but the number of constraints is n + n( n - 1 )/2.

circumscribecircles[n_Integer, func_]:=
Module[
{cons, cons1, cons2, vars, lb, ub, soln, rcirc, pd},
Do[r[i] = func[i], {i, n}];
cons1 = 
              Flatten [
                Table[r[i] + 
                      r[j]  ≤   (Sqrt[(((x[i] - x[j]))^2) + (((y[i] - y[j]))^2)]), {i, n - 1}, {j, i + 1, n}]];   
cons2  =  
              Table[obj  - 
                    r[i]  ≥ (Sqrt[ x[i]^2 + y[i]^2]), {i, n}];    
vars = 
              Join[Flatten[Table[{x[i], y[i]}, {i, n}]], {obj}];   
cons  =  Join[cons1, cons2];   
lb = Join[Table[-2.5, {2*n}], {1}];   
ub = Join[Table[2.5, {2*n}], {3}];   
{obj, cons, 
                Transpose[{vars, lb, ub}]} /.  
              x_[i_] :>   ToExpression[
                  ToString[x] <> ToString[i]] 
]

Packing 5 circles

A comparison for 5 circles (11 variables, 10 non-convex and 5 convex constraints) where the radius of the i'th circle is (i^(-1/2)) is shown below., The best results were achieved by the three MathOptimizer Professional methods, with similar results found by NMiniize SimulatedAnnealing and Global Optimization GlobalPenaltyFn. Notice also the relative runtime differences.

\(plt\ @\ 
      compareOpts[
        Sequence\  @@ \ circumscribecircles[5, #\^\(\(-1\)/2\) &], {1, 1, 1, 
          1, 1, 1, 3, 3}]\)

Packing 10 Circles

For 10 circles, the best results are achieved again by the three MathOptimizer Professional solver options, using also the least time. See the results and figure below.

\(plt\ @\ 
      compareOpts[
        Sequence\  @@ \ circumscribecircles[10, #\^\(\(-1\)/2\) &], {1, 1, 1, 
          1, 1, 1, 3, 3}]\)

Summary of Results and Conclusions

Global optimization covers a very broad class of problems.  Based on the test models used in our numerical experiments, we can draw some (admittedly limited) conclusions.  For small and relatively easy problems all solvers work fairly well, and Global Optimization GlobalSearch seems to be the fastest and most accurate.  For the native function NMinimize, the fastest results are found by the NelderMead method and the most accurate results by the slower DifferentialEvolution method. The methods in MathOptimizer Professional need an external compilation/link operation which in small problems dominates the total runtime for these solvers.

Regarding the more complex test problems, an aggregated comparison of the results found for the Lennard-Jones problem, the hanging chain problem and packing circles is shown in the following figure.  The circle packing results are at the left, curving upwards.  The hanging chain results are along a straight line with a less than 45 degree slope and the Lennard-Jones problem results fall along the x-axis.  For the hanging chain problem (which is a local search problem), the NMinimize NelderMead method gives the best results.  For circle packing, which has a many difficult constraints and many local minima, MathOptimizer Professional MultiStart gives the best results for more than 3 circles.  For the Lennard-Jones problem, which has no constraints but many local minima, the winning optimizer changes from Global Optimization GlobalPenaltyFn to MathOptimizer Professional Branch and Bound to MathOptimizer Professional MultiStart as the number of variables increases, with some winning cases from Global Optimization Global Search.

References

Loehle Enterprises Global Optimization 5.1 www.loehleenterprises.com, www.wolfram.com/products/applications/globalopt

Pintér, J.D. (1996…2005) LGO: A Model Development System for Continuous Global Optimization. User’s Guide. (Current revised edition.) Pintér Consulting Services, Inc., Halifax, NS. 

Pintér, J.D. (2005) Nonlinear optimization in modeling environments: software implementations for compilers, spreadsheets, modeling languages, and integrated computing systems. In: Jeyakumar and Rubinov, Eds., Continuous Optimization: Current Trends and Applications, pp. 147-174. 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 Users with an External Solver Link. User Guide. Pintér Consulting Services, Inc., Halifax, NS, Canada. For a summary, see also www.wolfram.com/products/applications/mathoptpro

Pintér, J.D. and Kampas, F.J. (2005a) Model development and optimization with Mathematica. In: Golden, B., Raghavan, S. and Wasil, E., Eds. The Next Wave in Computing, Optimization, and Decision, pp. 285-302. Springer Science + Business Media, New York.

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

Pintér, J.D. and Kampas, F.J. (2005c) MathOptimizer Professional: key features and illustrative applications. In: Liberti, L., and Macaulan, N., eds., Global Optimization: From Theory to Implementation. Springer Science + Business Media, New York. (To appear.)