GLOBAL OPTIMIZATION 4.2

GLOBAL NONLINEAR OPTIMIZATION USING MATHEMATICA

Loehle Enterprises
1258 Windemere Ave.
Naperville, IL 60564
USA
630-579-1190
630-579-1195 fax
craigloehl@aol.com


Version 4.2

(C) Copyright 1998-2002 Loehle Enterprises.  No part of this manual or software may be copied or reverse engineered except that purchasers of this package may make back up copies of the software.

*****************************************************************
DISCLAIMER

*****************************************************************

NOTICE: The software product described in this manual is a research tool.  No warranty is offered as to operability of this product.  No guarantee can be made that any particular nonlinear equation can be solved.  No safety critical or financial decisions should be based on the operation of this software.  The manufacturer accepts no liability for losses resulting from its use.
*****************************************************************

i  CONTENTS

TABLE OF CONTENTS

  I	Getting Started                                                                                             	  4

	1  Welcome                                                                                                    	  4
	2  Registering                                                                                                 	  4
	3  Installation                                                                                                    5
	4  A Quick Start and Example							 5

II	Tips for Performance                                                                         		  9

	1  Improving Performance                                                                               	10
	2  The Compile Function                                                                               	11
	3  User Functions: Working with Named Functions				12
	
 III	General Nonlinear Optimization: GlobalSearch, IntervalMin, 
 	GlobalPenaltyFn, and MultiStartMin						13

	1  Introduction								            13
	2  Example									15
	3  Program Operation.								16
		A  Parameters and Default Values					16
		B  Bounds								17
		C  Constraints								17
	4  Output									23
	5  Maximization								23
	6  Limitations									24
	7  Error Messages								26
	8  Performance/Speed								28
	9  Testing and Examples							33

	
 IV	The GlobalMinima Function							41

	1  Introduction							           		 41
	2  The Grid Refinement Approach to Nonlinear Optimization	           		 41
	3  Program Operation.								43
		A  Parameters and Default Values					43
		B  Bounds								45
		C  Constraints								45
		D  Memory Usage							46
	4  Output									46
	5  Maximization								46
	6  Limitations									47
	7  Efficiency									48
	8  Error Messages								49
	9  Step-by-Step Example							50
	10 Examples									56

V	Nonlinear Regression: the NLRegression Function              			67

	1  Introduction									67
	2  Utilizing chi-square fit criteria						74
	3  Multiple Independent Variables Regression					76

VI	Maximum Likelihood Estimation: the MaxLikelihood Function		82

	1  Introduction									82
	2  Examples and Built-In Functions						82

VII	Discrete Variable Problems: 
		InterchangeMethod & TabuSearchMin					89

	1  Introduction									89
	2  Applications									90
		A  Capital Allocation							90
		B  Vehicle Routing/Travelling Salesman				91
		C  Minimum Spanning Tree						95
	
 VIII	The MaxAllocation Function							96

	1  Introduction									96
	2  Applications									99
		A Investment Allocation Problems					99
		B Derivative Hedging with Quadratic Programming			100

IX 	Applications									101
								
	1  Zeros of a Function/Roots of a Polynomial					101
	2  Integer Programming: the Knapsack Problem				103
	3  Differential Equation Models						104
	4  Constraint Equations							106

X	Literature Cited								108

I  GETTING STARTED

I.1  Welcome

Welcome to Global Optimization.  This package provides a suite of tools for solving nonlinear optimization problems, as well as a variety of other applications such as finding the roots or zeros of a nonanalytic function.  The package is easier, simplier, and more robust than most optimization tools, and you will find yourself working more efficiently and more easily than before.

I.2  Registering

Before you continue please register by emailing your name and address with a notice that you have purchased Global Optimization to craigloehl@aol.com.  This will enable us to reach you with free updates and other news.  The developer does not receive customer information when the package is purchased from Wolfram Research or from distributors.  Please include the version of Global Optimization that you purchased and the machine you are using it on.

I.3  Installation

I.3.A  Recommended System Capabilities

This package is designed to work with Mathematica, which must be installed on your computer.  The following are recommended system capabilities for installation:

	•	Personal computer, 200 MHz or faster
	•	16M free RAM or more
	•	Mathematica 3.0 or higher installed

I.3.B  Installation on a Hard Drive

To install the package on your hard drive on any machine, copy the contents of the disk to the Mathematica directory.  The file may need to be unzipped.  On Mac the file may be a self-extracting archive; double-click it to extract the contents.  The go42.mx file contains Mathematica code needed to run the program in Encoded format.  This encoded file can not be read and will not work on other platforms.  It needs to be placed in the 3.0 or 4.0 main directory or in the AddOns/Applications directory.  An example notebook file .nb is included for illustration.  The .nb file(s) may be placed anywhere.

I.4  A Quick Start and Example

Version 4.2 contains ten functions: GlobalSearch, IntervalMin, GlobalPenaltyFn, MultiStartMin, GlobalMinima, NLRegression, MaxLikelihood, InterchangeMethodMin, TabuSearchMin, and MaxAllocation.  The package has been designed for easy use with Mathematica.  It uses all valid Mathematica functions and syntax.  The user must have Mathematica installed.  This manual assumes basic familiarity with Mathematica.  Example input files are in the form of Mathematica notebooks, with a file designation "filename.nb".  These notebook files correspond to examples in this manual.  This manual is executable.  Timing examples in this manual are based on a 600 Mhz Pentium III machine using Mathematica 4.0.

	To begin the execution of a notebook, follow this procedure:

1.	Start Mathematica.

2. 	Use the File pulldown menu to Open one of the example notebook files supplied with the disk.  

3.	Once this notebook is open, pull down the Kernel menu and select Evaluate Notebook.

	In the examples, the package is installed with the Get (Get["go42.mx"] or <<go42.mx) command, with definition for the function and options given as:

Off[Syntax::"spell"];
Off[Syntax::"spell1"];
$HistoryLimit = 0;

<< go42.mx

then several variables are defined, the function to be evaluated is defined, and then the problem is solved.  This process should execute with no errors and should produce output that resembles closely what was in the notebook when you opened it.  If not, some error may exist.  The basic formats are similar to the Mathematica function ConstrainedMin.  The following definitions define the functions in the package:

?GlobalMinima

?GlobalSearch

?GlobalPenaltyFn

?InterchangeMethodMin

?IntervalMin

?MaxAllocation

?MaxLikelihood

?MultiStartMin

?NLRegression

?TabuSearchMin

The following options are available for one or more functions.

?CompileOption

?EvaluateObj

?FastStepping

?MaxIterations

?SensitivityPlots

?ShowProgress

?SimplifyOption

?Starts

?StartsList

?TabuListLength

?UserMemory

?UserResiduals

?Weights

GlobalSearch[expression,inequalities,equalities,{{var1name,lowbound,highbound}..},tolerance,options]

GlobalPenaltyFn[expression,inequalities,equalities,{{var1name,lowbound,highbound}..},tolerance,options]

MultiStartMin[expression,inequalities,equalities,{{var1name,lowbound,highbound}..},tolerance,options]

InterchangeMethodMin[expression,inequalities,{varlist},tolerance,options]

IntervalMin[expression,inequalities,Null,{InitialPoint__Interval},varlist,tolerance,options]

TabuSearchMin[expression,inequalities,{varlist},tolerance,options]

GlobalMinima[expression,inequalities,{{var1name,lowbound,highbound}..},grid,tolerance,
contraction,indifference,options]

MaxAllocation[f[var1name,...],rhs,{varlist}..,tolerance,options]

NLRegression[data,expression,independent_variables,inequalities,{{var1name,lowbound,highbound}..},tolerance,options]

MaxLikelihood[data,expression,independent_variables,inequalities,{{var1name,lowbound,highbound}..},tolerance,options]

Where "expression" is the equation or function name for the function to be minimized, "inequalities" is a list of inequality constraints, and "equalities" is a list of equalities.  Speed of execution is enhanced with Mathematica versions 3.0.1 and higher.  The Mathematica kernel and the executable program need not reside on the same machine.

The functions are each designed for a specific type of problem.  For linear problems, the user should refer to ConstrainedMin (see Mathematica book), although the functions in this package will solve linear problems.  For smooth, unconstrained nonlinear problems, the user should try FindMinimum, which is faster for such problems.  The functions in this package are designed for nonlinear problems with local minima, multiple minima, and/or linear or nonlinear constraints.

The function GlobalSearch is a multiple-start generalized hill-climbing algorithm designed to work with or without constraints.  It is robust to noisy functions and local minima.  It can handle large problems (200+ variables).  It can handle linear and nonlinear equality and inequality constraints, which are assumed to be analytic.   GlobalSearch is the general purpose optimizer in the package, and is the engine used for the NLRegression and MaxLikelihood functions.  

	The function IntervalMin is a general minimizer using Interval methods.  It can handle inequality constraints, which must be analytic.  It is robust to local minima, but is slower than other methods.

The function GlobalPenaltyFn is a multiple-start generalized hill-climbing algorithm designed to work with or without constraints.  It handles the special case of constraints that are not Solvable for any of the variables, or that are black box.  It is robust to noisy functions and local minima.  It can handle large problems (200+ variables).

The function MultiStartMin is a restricted case version of GlobalSearch designed to work for problems with integer or discrete variables or that are highly nonlinear.  To solve this subset of problems, it handles constraints differently and is thus slower than GlobalSearch, particularly as problems get bigger.  Thus for problems with more than 15 variables, one should use GlobalSearch.  It is a multiple-start generalized hill-climbing algorithm designed to work with or without constraints.  It is robust to noisy functions and local minima.  It can handle linear and nonlinear inequality constraints. Equality constraints must be analytic, but inequality constraints need not be and can be black box or logical.  Objective functions may contain any combination of real, integer, and discrete variables.

The function GlobalMinima is designed for problems with many true or local minima or for which a region rather than a point better describes the optimum solution.  It can handle linear and nonlinear inequality constraints.  It is limited to smaller problems (<14 variables) but is very robust to noise and false minima.  For smaller problems (<3 variables), it is faster than MultiStartMin or GlobalSearch.  GlobalMinima can find solution regions, whereas MultiStartMin and GlobalSearch are not efficient at this task.

The function MaxAllocation is designed for allocation problems, such as arise in finance.  In an allocation problem, a fixed quantity, such as an investment sum, is to be distributed to a portfolio.  All variables must be nonnegative and the sum of all investments must equal the total available for investing.  This creates a nonlinear problem with a single equality constraint.

The function InterchangeMethodMin is designed for 0-1 integer problems such as arise in networks, transportation, and scheduling.  The operation of this routine is faster than for problems with continuous variables of the same size.  The function can solve vehicle routing/traveling salesman, minimum spanning tree, and capital allocation problems, among others.

The function TabuSearchMin is designed for 0-1 integer problems such as arise in networks, transportation, and scheduling. The operation of this routine is faster than for problems with continuous variables of the same size. The function can solve vehicle routing/traveling salesman,minimum spanning tree, and capital allocation problems, among others.

The function NLRegression solves nonlinear regression problems.  It is particularly designed for noisy problems or those requiring constraints to achieve a good fit.  L1 and L2 norms are allowed.  Least-squares and chi-square options are also available.  Sensitivity plots and confidence intervals on the parameters are provided as output.

The function MaxLikelihood solves maximum likelihood estimation problems using Log- likelihood estimation.  A library of common distributions is provided and fit statistics are computed.  Summary statistics are computed.

II  TIPS FOR PERFORMANCE

A key factor for numerical optimization is the time taken to compute the objective function.  This is critical because the objective function may be called thousands of times to hundreds of thousands of times.  This section discusses tricks to achieve better performance.

II.1  Improving Performance

Some general tips include the removal of the reading of external files from the user function, computing constant expressions once and then using the result, and avoiding logical operations that will prevent Compile from being used.  A key is to generate a function once as an explicit expression, and then pass this in to the package.  If many steps are involved in generating a user function and these must be performed every time the function is called, this can be very expensive.  It often happens that parts of an expression involve constants.  If these can be evaluated up front, the savings can be enormous.  In the following example, the function "f" evaluates numerically when Evaluate is executed during Function creation.  In function "g" functions like Sin and Tan are evaluated each time the function is called:

ClearAll[f, g]

x = {x1, x2, x3, x4, x5, x6, x7, x8, x9, x10};

f = 
      Function[Evaluate[x], 
        Evaluate[
          Sum[Sin[Pi/2.1]*x[[i]] + Tan[19.4*Pi]/x[[i]] + Cos[ .1] + 
              Sqrt[99], {i, 1, 10}]]]

g[p_] := 
      Module[{d}, {x1, x2, x3, x4, x5, x6, x7, x8, x9, x10} = p; 
        r = Sum[Sin[Pi/2.1]*x[[i]] + Tan[19.4*Pi]/x[[i]] + Cos[ .1] + 
              Sqrt[99], {i, 1, 10}]; Return[r]]

dat = Table[i, {i, 1, 10}]

Apply[f, dat]

g[{1, 2, 3, 4, 5, 6, 7, 8, 9, 10}]

Do[Apply[f, dat], {i, 1, 500}] // Timing

Do[g[{1, 2, 3, 4, 5, 6, 7, 8, 9, 10}], {i, 1, 500}] // 
      Timing

Because the constants in f were Evaluated before execution (during Function definition), f was very fast to execute, even without Compile.  In contrast, g, which computes Sin etc. every time it is called, took 46 times longer to compute.

II.2  The Compile Function

The Mathematica function Compile is used to speed up program execution.  The user input function and constraints are both Compiled.  The advantage of the Compile function is greatest with more complex expressions and larger problems.  For small problems, it can be difficult to even detect a benefit, but there is little cost to the Compile.  The range of speed improvement due to Compile can range from 30% to 30x.  Compilation does not work for all Mathematica expressions.  For example, if "Apply" or "If" are used in the function definition, the function will not Compile.  It is suggested that if a complex expression involving Mathematica special functions is to be used, the user attempt to Compile the expression to test it.  The CompileOption can be set to False (CompileOption->False) in the function call if the user function is not compilable.  It may also be set to False if the user Compiles the expression before passing it in.  It is possible to get errors during the Compile step.  These will be highlighed in blue and will indicate that they are compilation errors.  It is also possible that execution errors can result if non-machine numbers result from some computation of the compiled function.  Mathematica usually reverts to the uncompiled expression in such cases, but the user may wish to run with CompileOption->False in this case to obtain the most accuracy and best speed.  If Compile is attempted on a user function with very complex expressions, it can sometimes hang up during the Compile stage due to a bug in the Mathematica kernel, so testing is again recommended. 

	We can see the benefit of the Compile in the following example:

ClearAll[x1, x2, x3, x4, x5, x6, x7, x8, x9, x10]

x = {x1, x2, x3, x4, x5, x6, x7, x8, x9, x10};

f = Sum[(1 + x[[i]])^2, {i, 1, 10}]

g1 = Function[Evaluate[x], Evaluate[f]];

g2 = Compile[Evaluate[x], Evaluate[f]];

Do[Apply[g1, { .1, 2. , 3. , 4. , 3. ,  .3,  .3,  .3,  .5,  .5}], {i, 
          1, 10000}] // Timing

Do[Apply[g2, { .1, 2. , 3. , 4. , 3. ,  .3,  .3,  .3,  .5,  .5}], {i, 
          1, 10000}] // Timing

We see that Compile speeded up the execution by a factor of more than 4.  On some functions, the speedup can be a factor of 20 or more.  Setting up a function as a Module prevents Compile from working on the objective function, and should be avoided.  Some Mathematica functions, particularly those that are math functions, can be included in a Compiled function.  Examples include Abs, Sqrt, Max, Min, Sin, etc.  Others, such as FindRoot, can not.

II.3  User Functions: Working with Named Functions

The objective (goal) function in all the functions in this package can be defined in several ways.  The function must in some cases be passed in differently, depending on the function format.  In the simplest case, the expression is itself passed in.  Here we find the minimum of the given function, which is zero.

ClearAll[x]

GlobalSearch[x^2 + y^2, , , {{x, 38, 40}, {y, 38, 40}},  .00000001, 
      Starts -> 1]

We can also pass in a function by name:

f := x^2 + y^2

GlobalSearch[f, , , {{x, 38, 40}, {y, 38, 40}},  .00000001, 
      Starts -> 1]

In order to Compile or use the objective function, the programs in this package will attempt to Evaluate the expression.  If there are certain operations involved in the user function, such as obtaining the Inverse or Determinant of a matrix, Mathematica will attempt to Evaluate these symbolically when the function is used in the program.  It will then work with the Evaluated function.  If the user expression is complex, however, this may become impossible.  For example, a symbolic inverse of a matrix with 10 variables may take a very long time.  The user can test this by trying to Evaluate the objective function before passing it in to any of the functions of this package.  If it fails to Evaluate, then use the parameter EvaluateObj->False to prevent evaluation.  Another time when EvaluateObj should be set to False is when any logical operations are performed in the objective function, as in the following example:

f = 
        Function[{x, y}, If[x + y < 0, tem = 0, tem = y + x]; 
          tem];

Apply[f, {-1, -1}]

Apply[f, {1, 1}]

Clearly, the function is minimized for any set with x + y negative or for x==-y.  Because of the conditional logic in the expression, it should not be symbolically evaluated.  This means that it should not be simplified (SimplifyOption->False) or Evaluated (EvaluateObj->False):

GlobalSearch[f, {}, , {{x, 38, 40}, {y, 38, 40}},  .000001, 
        Starts -> 1, EvaluateObj -> False, 
        SimplifyOption -> False] // Timing

If, in contrast, we set up a function as a Module, the options do not affect the result, but the function can not be Compiled, and thus speed will suffer.

g[x_, y_] := 
        Module[{d}, If[x < 0 && y < 0, tem = 0, tem = 1 + Abs[y] + Abs[x]]; 
          Return[tem]];

GlobalSearch[g, , , {{x, 38, 40}, {y, 38, 40}},  .000001, 
        Starts -> 1, EvaluateObj -> False] // Timing

Thus, if there are any If statements in the objective function or it is a black box model, set EvaluateObj->False or set it up as a Module.

III  GENERAL NONLINEAR OPTIMIZATION: GlobalSearch, GlobalPenaltyFn, IntervalMin AND MultiStartMin

III.1  Introduction

Five general nonlinear solvers are included in this package.  They each are designed for a particular class of problems.  These functions provide tools for global optimization.  Traditional gradient (local) approaches require the user to know how many optima are being sought for the function to be solved, and roughly where the optima are, so that a good initial guess can be made.  The user, however, rarely can make a good initial guess and usually has no information about the existence of multiple solutions.  In the absence of such information, existing algorithms will generally converge to an answer, but this may be only a local solution that may not be globally optimal.  Furthermore, even global optima may not be unique.  The functions in this package can solve black-box objective functions, objective functions that are nondifferentiable, differential equation models, and functions with discrete steps in them.  The initial bounds on the parameters given on input do not need to bound the true solution (except for the GlobalMinima function), but only help get the algorithms started in a good region.

GlobalSearch approaches the difficult problem of finding a global optimum with several techniques.  A generalized hill climbing technique is used that is based on Newton's method but using a generalized gradient rather than a derivative, and allowing for constraints.  Constraints must be analytic, but can be linear or nonlinear.  Multiple starts are used to test for the existence of multiple solutions.  The multiple starts are generated randomly from the region defined by the range of parameter values input by the user.  Feasible starting regions are not needed, but it is assumed that objective function values in this region are Real.  When a step can not be made that improves the solution by at least "tolerance", the program terminates.  If a problem uses only discrete 0-1 variables, the functions InterchangeMethodMin and TabuSearchMin should be used.  GlobalSearch is the solver used by the regression and maximum likelihood functions.  The function is defined by

GlobalSearch[expression,inequalities,equalities,{{var1name,lowbound,highbound}..},tolerance,options]

GlobalPenaltyFn addresses the problem of constraints that are nonanalytic.  Nonanalytic functions can be algorithmic, can have conditional logic, or can be too complex for Solve to separate variables.  If GlobalSearch fails on the constraints, use GlobalPenaltyFn.  A generalized hill climbing technique is used that is based on Newton's method but using a generalized gradient rather than a derivative, and allowing for constraints.  An adaptive penalty method is used that adjusts the penalties for constraint violation depending on the degree of violation of constraints.  Multiple starts are used to test for the existence of multiple solutions.  The multiple starts are generated randomly from the region defined by the range of parameter values input by the user.  Feasible starting regions are not needed, but it is assumed that objective function values in this region are Real.  When a step can not be made that improves the solution by at least "tolerance", the program terminates.  If a problem uses only discrete 0-1 variables, the functions InterchangeMethodMin and TabuSearchMin should be used.  The function is defined by

GlobalPenaltyFn[expression,inequalities,equalities,{{var1name,lowbound,highbound}..},tolerance,options]

IntervalMin uses Interval methods to find a solution.  Intervals are computed in Mathematica based on concepts of limits.  An initial Interval is input by the user.  It is beneficial for this Interval to cover the true solution point, but this is not necessary as long as it is not orders of magnitude distant.  Both the input objective function and the constraints must be analytic.  IntervalMin is robust to local minima, but is not as fast as GlobalSearch because Compile can not be used with Interval variables.  IntervalMin can not solve problems with equality constraints at this time.  We illustrate Intervals below:

x^2 /. x -> Interval[{-1, 2}]

In this case, the true minimum of the function is given by the lower bound of the Interval result.  In simple or separable cases, this will pertain.  In other cases it will not:

1.` - 2.` x + (x^2) + 100.` (x^4) - 200.` (x^2) y + 
        100.` (y^2) /. {x -> Interval[{0, 2}], 
        y -> Interval[{0, 2}]}

The true minimum of this function is 0 at {1,1}.  The true minimum is within the Interval of the result in this case.  Note that when IntervalMin is used, both the estimated parameters and the function value are given in terms of an Interval.  The format for the function is given by:

IntervalMin[expression,inequalities,Null,{InitialPoint__Interval},varlist,tolerance,options]

MultiStartMin addresses the problem of inequality constraints that are nonanalytic. MultiStartMin is effective for highly nonlinear functions with fewer than 15 variables. MultiStartMin may have trouble with multiple constraints if the solution lies in the corner of two constraints.  In such cases, use IntervalMin or GlobalPenaltyFn.  A generalized hill climbing technique is used that is based on Newton's method but using a generalized gradient rather than a derivative, and allowing for constraints.  Inequality constraints are treated as hard boundaries.  Equality constraints must be analytic.  Integer variables should not be used in Equality constraints unless the equality constraints are linear.  All combinations of 2 or more variables at a time (number defined by option SearchDirections) are tested to find the best search direction.  Multiple starts are used to test for the existence of multiple solutions.  The multiple starts are generated randomly from the region defined by the range of parameter values input by the user.  Feasible starting regions are not needed, but it is assumed that objective function values in this region are Real.  Mixtures of Real, Integer, and Discrete (defined by a List) variables are allowed.  If a problem uses only discrete 0-1 variables, the functions InterchangeMethodMin and TabuSearchMin should be used.  The function is defined by

MultiStartMin[expression,inequalities,equalities,{{var1name,lowbound,highbound}..},tolerance,options]

The GlobalMinima function uses an adaptive grid refinement technique.  This technique is robust, and can find multiple solutions in a single run.  Because it's operation is so different, it is documented in it's own section, later.

GlobalMinima[expression,inequalities,{{var1name,lowbound,highbound}..},grid,tolerance,
contraction,indifference,options]

III.2  Example

A simple example is to minimize a square function subject to a constraint, as follows.  The constraint

x > 1 && y > 1

must be converted to standard form as a list:

-x < -1 && -y < -1

-x + 1 < 0 && -y + 1 < 0

c = {-x + 1, -y + 1}

GlobalSearch[x^2 + y^2, 
        c, , {{x, -5, 5}, {y, -5, 5}},  .000001] // Timing

The output shows three solutions from random initial points, the default, with the {x,y} value and the function value in each list element.  Initial feasible points are not needed.  In this case, the function is smooth, so all three solutions are the same, and there is really no need to use multiple starts.

III.3  Program Operation

III.3.A  Parameters and Default Values

The tolerance (T) defines the amount of improvement (absolute, not relative) in the function value required at each iteration.  If at least this much improvement is not found, the program stops.  The tolerance can be set to very small values, such as (10^(-15)), to achieve a highly accurate estimate of the minimum.  If the user sets tolerance = 0, this is an error.

The functions have seven options.  The defaults are MaxIterations->10000, CompileOption->True, ShowProgress->True, StartsList->{}, EvaluateObj->True, SimplifyOption->True, and Starts->3.  MaxIterations prevents the program from running away to infinity when a mistake is made in the function (Min[x] for example).  For very computationally expensive problems (like those with many parameters), it is useful to use Starts->1 to evaluate the time a run takes.  For general usage, Starts->3 (the default) is good. The results of this short run can then be used to define further, more directed runs based on whether all 3 starts found the same solution.  A list of starting values can be input instead of letting the program find them with random search (e.g., StartsList->{{1.,2.},{2.,1.}}).  The CompileOption determines whether the user function and constraints will be compiled.  While Compile reduces execution time, some functions can not be compiled, in which case CompileOption->False should be used.  If the objective function should not be Evaluated, use EvaluateObj->False.  SimplifyOption attempts to simplify the objective function.  If this should not be done or will not improve the solution, this should be set to False.  The MultiStartMin function also has a SearchDirections (default 2) option.  More search directions help solve highly nonlinear problems, but increase the execution time exponentially.

III.3.B  Bounds

Whereas most algorithms require input of an initial guess or starting point (which must be feasible) to initiate the search, GlobalSearch requires only bounds on the variables.  Bounds are generally easier to define than are feasible starting points.  Physical, logical, or economic considerations often provide guidance on realistic bounds for a variable.  If the optimization problem is wing design and one variable is wing thickness, a certain range of thicknesses is sensible to consider.  For a financial problem, the range of values for an investment mix across a product line is clearly bounded between zero and the total budget, and the upper limit can generally be bounded below a much smaller value than this.  If narrower bounds are known, then convergence will require fewer function calls, but even quite wide bounds are acceptable.  If erroneously narrow bounds are input, then initial step size will be too small, and the function will act more like a local solver.  It is better to err on the side of wider bounds.  Upper and lower bounds must be input by the user, as illustrated below.  These bounds do not restrict the searching of GlobalSearch, MultiStartMin, IntervalMin, or GlobalPenaltyFn (as they do for GlobalMinima) but merely provide an initial guide for generating feasible starts.  Initial bounds do not need to produce feasible solutions, the program can find initial feasible points, but the values within these bounds are assumed to be Real.  If hard bounds on variables are necessary (such as positivity restrictions), they can be entered as constraints.

III.3.C  Constraints

Constraints are entered as Mathematica functions, which may be linear or nonlinear.  Equality constraints are entered in the third position in standard form.  For example, positivity restrictions on x and y would be represented by the list:

{-x, -y}

The following problem is a typical LP problem, with the solution at the intersection of two inequality constraints.  The solution is {2/3,10/3,0}).

GlobalSearch[-x1 - 2  x2 + 
          x3, {x1 + x2 + x3 - 4, -x1 + 2  x2 - 2  x3 - 6, 
          2  x1 + x2 - 5, -x1, -x2, -x3}, , {{x1, 0, 1}, {x2, 0, 
            1}, {x3, 0, 1}},  .0000001, Starts -> 2] // 
      Timing

IntervalMin[-x1 - 2  x2 + 
          x3, {x1 + x2 + x3 - 4, -x1 + 2  x2 - 2  x3 - 6, 
          2  x1 + x2 - 
            5, -x1, -x2, -x3}, , {Interval[{0. , 4. }], 
          Interval[{0. , 4. }], Interval[{0. , 1. }]}, {x1, x2, 
          x3},  .0000001] // Timing

We see that IntervalMin is much slower, but there are problems where it performs better.

GlobalPenaltyFn[-x1 - 2  x2 + 
          x3, {x1 + x2 + x3 - 4, -x1 + 2  x2 - 2  x3 - 6, 
          2  x1 + x2 - 5, -x1, -x2, -x3}, , {{x1, 0, 3}, {x2, 0, 
            3}, {x3, 0, 1}},  .000000001, Starts -> 3] // 
      Timing

We see above that GlobalPenaltyFn is slower than GlobalSearch because it can take nonanalytic constraints.  It is also not quite as robust.  MultiStartMin is able to solve this problem only from the ideal starting value.

MultiStartMin[-x1 - 2  x2 + 
          x3, {x1 + x2 + x3 - 4, -x1 + 2  x2 - 2  x3 - 6, 
          2  x1 + x2 - 5, -x1, -x2, -x3}, , {{x1, 0, 3}, {x2, 0, 
            3}, {x3, 0, 1}},  .0000001, Starts -> 2] // 
      Timing

It is useful to illustrate a function in which the constraints are nonanalytic.  Any function for which Solve will not work falls in this category.  We take as the example:

Min[y^2 + x]    s . t .    x == Sin[y]

The above equality can not be Solved for y uniquely.  Thus GlobalSearch  and MultiStartMin fail but GlobalPenaltyFn succeeds.  We first solve the problem by substitution:

res = 
      GlobalSearch[y^2 + Sin[y], , , {{y, 0, 1}},  .00001, Starts -> 1, 
          ShowProgress -> False] // Timing

Where x is then

x = N[Sin[y]] /. res[[2, 1, 1]]

ClearAll[x]

res = 
      GlobalSearch[y^2 + x, , {x - Sin[y]}, {{x, 0, 1}, {y, 0, 1}},  .0001, 
          Starts -> 1, ShowProgress -> False] // Timing

res = 
      GlobalPenaltyFn[
          y^2 + x, , {x - Sin[y]}, {{x, 0, 1}, {y, 0, 1}},  .000001, 
          Starts -> 1, ShowProgress -> False] // Timing

In the next example, we see how a nonlinear constraint can be defined that creates a whole set of solution points, an approximation to which can be obtained with enough starts.  In this example, the solution must lie outside the circle of radius Sqrt[2], but the unconstrained solution is at the origin.  This means that all points on this circle are optimal.

ClearAll[x]

b = 
        GlobalSearch[
          x^2 + y^2, {2 - (x^2 + y^2)}, , {{x, -10, 10}, {y, -10, 
              10}},  .0001, Starts -> 100, 
          CompileOption -> True];

b[[1]]

c = {}; 
    Do[AppendTo[c, {x, y} /. b[[ii, 1]]], {ii, 1, 100}]

ListPlot[c];

In this case, the program finds only a subset of the solutions because of the order in which multiple roots are processed.  MultiStartMin finds a better sample of points around the circle in this case.

b = 
        MultiStartMin[
          x^2 + y^2, {2 - (x^2 + y^2)}, , {{x, -10, 10}, {y, -10, 
              10}},  .0001, Starts -> 100, 
          CompileOption -> True];

b[[1]]

c = {}; 
    Do[AppendTo[c, {x, y} /. b[[ii, 1]]], {ii, 1, 80}]

ListPlot[c];

Equality constraints are input as a list of equality statements in the third parameter position.  In the following problem, the solution is {1/3,1/3,1/3}, with f=1/3.

GlobalSearch[
        x^2 + y^2 + 
          z^2, {}, {x + y + z - 1}, {{x, 0, 1}, {y, 0, 1}, {z, 0, 
            1}},  .00000001, Starts -> 1, ShowProgress -> False] // 
      Timing

MultiStartMin[
        x^2 + y^2 + 
          z^2, {}, {x + y + z - 1}, {{x, 0, 1}, {y, 0, 1}, {z, 0, 
            1}},  .00000001, Starts -> 1, ShowProgress -> False] // 
      Timing

GlobalPenaltyFn[
        x^2 + y^2 + 
          z^2, {}, {x + y + z - 1}, {{x, 0, 1}, {y, 0, 1}, {z, 0, 
            1}},  .00000001, Starts -> 1, ShowProgress -> False] // 
      Timing

In some cases, the solution to a problem occurs at the intersection of constraints.  These problems are difficult but can be solved.  In the following problem, the solution is exactly at the intersection of the line and the circle:

Solve[{x1 + x2 - 1.2 == 0, x1^2 + x2^2 == 1}]

GlobalSearch[(x1 -  .1)^2 + (x2 -  .1)^2 + 
          x3^2, {x1 + x2 - 1.2}, {x1^2 + x2^2 - 1}, {{x1, 0, 1}, {x2, 0, 
            1}, {x3, 5, 10}},  .0000001, Starts -> 4, 
        ShowProgress -> False] // Timing

We see that GlobalSearch found both solutions with 4 starts.  GlobalPenaltyFn also finds both solutions in 4 starts, but 2 of the solutions are not quite dead on.

GlobalPenaltyFn[(x1 -  .1)^2 + (x2 -  .1)^2 + 
          x3^2, {-x1, -x2, 
          x1 + x2 - 1.2}, {x1^2 + x2^2 - 1}, {{x1, 0,  .5}, {x2, 
            0,  .5}, {x3, 5, 10}},  .0000001, Starts -> 4, 
        ShowProgress -> False] // Timing

MultiStartMin finds only one solution, but the solution is correct.

MultiStartMin[(x1 -  .1)^2 + (x2 -  .1)^2 + 
          x3^2, {x1 + x2 - 1.2}, {x1^2 + x2^2 - 1}, {{x1, 0,  .5}, {x2, 
            0,  .5}, {x3, 5, 10}},  .0000001, Starts -> 4, 
        ShowProgress -> False] // Timing

In the following example, the constraint forces a solution outside the unit circle.  The solution is {.707,.707}.  MultiStartMin fails on this problem.  IntervalMin gets pretty close, as does GlobalPenaltyFn.

MultiStartMin[(x -  .1)^2 + (y -  .1)^2, {-(x^2 + y^2 - 
                1)}, {}, {{x, -4, 5}, {y, -4, 5}},  .0000001, 
        Starts -> 3] // Timing

GlobalSearch[(x -  .1)^2 + (y -  .1)^2, {-(x^2 + y^2 - 
                1)}, {}, {{x, -4, 5}, {y, -4, 5}},  .000000001, 
        Starts -> 3] // Timing

IntervalMin[(x -  .1)^2 + (y -  .1)^2, {-(x^2 + y^2 - 
                1)}, , {Interval[{0. , 1. }], Interval[{0. , 1. }]}, {x, 
          y},  .0001] // Timing

GlobalPenaltyFn[(x -  .1)^2 + (y -  .1)^2, {-(x^2 + y^2 - 
                1)}, {}, {{x, -4, 5}, {y, -4, 5}},  .000000001, 
        Starts -> 3] // Timing

WARNING:  Constraints can cause difficulties for GlobalSearch.  Constraints can be mistakenly formulated such that no feasible space exists.  In this case, GlobalSearch will stop with an error message:

GlobalSearch[
        x^2 + y^2 + 
          z^2, {x - 1, -x + 2}, {}, {{x, 0, 1}, {y, 0, 1}, {z, 0, 
            1}},  .00000001, Starts -> 1, ShowProgress -> False] // 
      Timing

WARNING: Equality constraints should not be used with Integer variables.  The result of doing so here is that they are converted into Real variables.

MultiStartMin[(x1 -  .1)^2 + (x2 -  .1)^2 + 
          x3^2, {x1 + x2 - 1.2}, {x1^2 + x2^2 - 1}, {{x1, 0, 1, 
            Integer}, {x2, 0,  .5}, {x3, 5, 10}},  .0000001, 
        Starts -> 1, ShowProgress -> False] // Timing

III.4  Output

Following solution, a list is returned that contains the solution.  For example, two points in a 2 parameter problem would give the list {{{1.1,1.2},-5},{{1.2,1.21},-4.9}}.  Intermediate output can be printed out by setting ShowProgress->True.

GlobalSearch[
      x^2 + y^2, {-x}, , {{x, -5, 5}, {y, -5, 5}},  .000000001, 
      Starts -> 1, ShowProgress -> True]

Following solution, a list is returned that contains the solution as a list of replacement rules with corresponding solutions.

III.5  Maximization

To find a maximum instead of a minimum, simply call GlobalSearch with the negative of the function, keeping in mind that a negative function value will be printed out.

III.6  Limitations

No optimization algorithm can solve all conceivable problems.  Consider a response surface in two dimensions that is flat everywhere except for a very small but very deep pit at one location.  Any hill-climbing algorithm will fail completely to solve this problem, because for almost all starting points there is no gradient.

Some functions have infinitely many solutions.  Consider Sin(1/x) between 0 and 1.  As the function approaches zero, the minima (with z = -1) get closer and closer together, as we can see from the plot following:

Plot [Sin[1/x],  {x,  0,  1}];

Although it is fundamentally impossible to list all solutions when there are infinitely many, GlobalSolve can show that there are very many solutions.  If we run the program with 40 starts on the interval [0, 0.6], we obtain many points with values close to -1.0 (shown below).  This result demonstrates that there are many solutions and that the density of solutions increases with an approach to zero.  This is thus an approximate solution to the true situation.  x is bounded away from 0 to prevent underflows.

GlobalSearch[
      Sin[1/x], {-x +  .00001, x -  .6}, , {{x, 0.0000001,  .6}},  .0001, 
      CompileOption -> False, Starts -> 20]

By restricting the search region, we find more solutions:

GlobalSearch[
      Sin[1/x], {-x +  .000000000001, 
        x -  .2}, , {{x, 0.000001,  .2}},  .00001, Starts -> 20]

These examples illustrate that not all problems can be solved.  Some are not solvable by any algorithm.  Some may be ill-posed.  However, for many problems that cause other algorithms to fail, GlobalSolve either succeeds, has a probability of succeeding, or provides a partial or approximate solution. We may thus say that the set of unsolvable problems is much smaller than it is for other solution methods.  This is particularly so with respect to the inclusion of constraints.

The optimization algorithm fundamentally assumes real valued functions.  It is possible for some parameter values that are tested to return values that are complex or otherwise not Real.  The program assumes that these are illegal values and treats them like values that fail a constraint test.  This means that it can solve a problem where values can be complex.  In some cases, however, complex values or other invalid results are detected on input as a sign of trouble, and the program aborts, as in the following:

GlobalSearch[x^ .5, {}, {}, {{x, -2, 2}},  .00001, 
      Starts -> 1]

When we correct the input range on x, the problem runs:

GlobalSearch[x^ .5, {}, {}, {{x, 0, 2}},  .00001, Starts -> 1, 
      CompileOption -> False]

III.7  Error Messages

An effort has been made to trap as many errors as possible.  Many user errors will prompt Mathematica error messages.  On input, only the order of parameters is evaluated in the passing function call.  If a valid value for tolerance (> 0.0) is not input, the following error message is printed, and the program stops:

	Error: tolerance must be > 0

If the upper bound of the variable definitions is smaller than the lower bound, the program will also stop with an error message.

If a parameter is omitted from the calling sequence, Mathematica cannot continue.  To indicate the problem, it echoes back the names for functions and values for constants in the function call, and then it stops.  This type of echoed-back output means that a format error exists in the calling sequence.  Some syntax errors are trapped within the Global Optimization package.

After input, the program compiles the user input function and the constraints, speeding execution.  If these functions are improperly formatted, this step may produce a compilation error message, which will terminate the run.  If a function is defined such that it does not return numeric values, this will cause error messages from Mathematica, followed by termination of GlobalSearch.

When constraints are used, there may be no feasible region, particularly if the problem has many constraints.  If this occurs, the following error message is printed, and the program terminates:

	Error: no valid initial points found

In this case, check the constraints.  Constraints can be in conflict if care is not taken.

A common mistake is to maximize when the intention was to minimize.  This error can be recognized because the solution will tend to run into one of the parameter bounds and will have an illogical value.  For example, the optimal airplane wing thickness for a design problem might come back as zero.

An unbounded solution may result from an improperly formulated problem (e.g., min(-1/x) over {-1, 1} which becomes -inf at x=0).  Because the function z continues to get larger the closer one gets to zero at an increasing rate, the program will never report convergence.  ShowProgress->True can be used to check for such problems during a run.

The user function passed in to GlobalSearch is Compiled.  This may cause compilation errors.  The result is usually to stop the execution of the program, but it is necessary for the user to realize the source of the error.  Compiled functions can also generate errors at run time if the user function generates non-machine numbers such as high precision numbers or infinity.  Mathematica will usually revert to the uncompiled function and continue running.

If the user defines variables in his program by mistake that should be parameters of the function, this will cause the program to stop or malfunction.  For example, if pp=5 is defined by mistake, this will cause GlobalSearch to detect an error which can be understood from looking at the error output which shows an invalid parameter list:

pp = 5;

\(\(\(GlobalSearch[\((pp - 30)\)^2\ \((y - 60)\)^2, , , {{pp, 0. , 
          100. }, {y, 1, 5}}, 0.001]\)\(
\)\)\)

Note how the program echoes back the parameter input, showing the mistake.

If the function to be solved is by mistake unbounded, then the parameter MaxIterations will stop the program from running forever.

GlobalSearch[
      x + y, {}, {}, {{x, -5, 5}, {y, -5, 5}},  .0000001, 
      Starts -> 1, MaxIterations -> 4]

MaxIterations can also be used to make a preliminary run to test out execution.

GlobalSearch[
      100. *(x^2 - y)^2 +  (1.  - x)^2, {}, {}, {{x, -5, 
          5}, {y, -10. , 10. }},  .0000001, Starts -> 1, 
      MaxIterations -> 5]

We now wish to restart the program with the given ending value.

r = %

rr = r[[1]]

s = {x, y} /. rr[[1]]

GlobalSearch[
        100. *(x^2 - y)^2 +  (1.  - x)^2, {}, {}, {{x, -5, 
            5}, {y, -10. , 10. }},  .0000001, StartsList -> {s}, 
        ShowProgress -> False] // Timing

III.8  Performance/Speed

For larger problems, speed can become a limiting factor.  It is suggested that the objective function and constraints be formulated such that they can be Compiled (see sec. II) which can speed up operation by 2x to 20x.  If a function is very complex, it can be helpful to pull out subexpressions that do not need to be computed each time, and compute them once in advance.  Simplify or FullSimplify may also speed up the function.  Some large problems are run next to illustrate performance.  Either GlobalSearch or GlobalPenaltyFn can be run for large unconstrained problems, for which they use the same basic algorithm.  MultiStartMin can only be used on large problems if SearchDirections->1 is used, and FastStepping->True is recommended for such cases.  Using MultiStartMin will not be advantageous for problems with highly nonlinear functions, but may be the only option for problems with integer or discrete variables.

vars = 
        Table[ToExpression[StringJoin["x", ToString[i]]], {i, 1, 
            100}];

fg = 
        Function[Evaluate[vars], 
          Evaluate[Sum[vars[[i]]^2, {i, 1, 100}]]];

vardefs = 
        Table[{vars[[i]], -10, 10}, {i, 1, 100}];

GlobalSearch[
        fg[x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15, 
          x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x26, x27, x28, 
          x29, x30, x31, x32, x33, x34, x35, x36, x37, x38, x39, x40, x41, 
          x42, x43, x44, x45, x46, x47, x48, x49, x50, x51, x52, x53, x54, 
          x55, x56, x57, x58, x59, x60, x61, x62, x63, x64, x65, x66, x67, 
          x68, x69, x70, x71, x72, x73, x74, x75, x76, x77, x78, x79, x80, 
          x81, x82, x83, x84, x85, x86, x87, x88, x89, x90, x91, x92, x93, 
          x94, x95, x96, x97, x98, x99, x100], , , vardefs,  .000001, 
        Starts -> 1, ShowProgress -> False] // Timing

vars = 
        Table[ToExpression[StringJoin["x", ToString[i]]], {i, 1, 
            200}];

fg = 
        Function[Evaluate[vars], 
          Evaluate[Sum[vars[[i]]^2, {i, 1, 200}]]];

vardefs = 
        Table[{vars[[i]], -10, 10}, {i, 1, 200}];

GlobalSearch[
        fg[x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15, 
          x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x26, x27, x28, 
          x29, x30, x31, x32, x33, x34, x35, x36, x37, x38, x39, x40, x41, 
          x42, x43, x44, x45, x46, x47, x48, x49, x50, x51, x52, x53, x54, 
          x55, x56, x57, x58, x59, x60, x61, x62, x63, x64, x65, x66, x67, 
          x68, x69, x70, x71, x72, x73, x74, x75, x76, x77, x78, x79, x80, 
          x81, x82, x83, x84, x85, x86, x87, x88, x89, x90, x91, x92, x93, 
          x94, x95, x96, x97, x98, x99, x100, x101, x102, x103, x104, x105, 
          x106, x107, x108, x109, x110, x111, x112, x113, x114, x115, x116, 
          x117, x118, x119, x120, x121, x122, x123, x124, x125, x126, x127, 
          x128, x129, x130, x131, x132, x133, x134, x135, x136, x137, x138, 
          x139, x140, x141, x142, x143, x144, x145, x146, x147, x148, x149, 
          x150, x151, x152, x153, x154, x155, x156, x157, x158, x159, x160, 
          x161, x162, x163, x164, x165, x166, x167, x168, x169, x170, x171, 
          x172, x173, x174, x175, x176, x177, x178, x179, x180, x181, x182, 
          x183, x184, x185, x186, x187, x188, x189, x190, x191, x192, x193, 
          x194, x195, x196, x197, x198, x199, x200], , , vardefs,  .000001, 
        Starts -> 1] // Timing

Note that for large problems with constraints or lots of data, execution will be slower.

III.9  Testing and Examples

In this section, the general nonlinear functions are tested.  The test criteria are the ability to obtain accurate answers, time required to execute, size of problem, ability to handle noisy functions, and ability to find multiple solutions.  In the first example, we test for the ability to solve to arbitrary precision:

\(\(\(GlobalSearch[\((x - 30)\)^2\ \((y - 60)\)^2, , , {{x, 0. , 
          100. }, {y, 1, 5}}, 0.000000000000000001, 
      Starts -> 1]\)\(
\)\)\)

The following problem is the difficult Rosenbrock function, but it is easily solved.

GlobalSearch[
        100. *(x^2 - y)^2 +  (1.  - x)^2, {}, {}, {{x, -5, 
            5}, {y, -10. , 10. }},  .0000001, Starts -> 1, 
        ShowProgress -> False] // Timing

MultiStartMin[
        100. *(x^2 - y)^2 +  (1.  - x)^2, {}, {}, {{x, -5, 
            5}, {y, -10. , 10. }},  .0000001, Starts -> 1, 
        ShowProgress -> False] // Timing

GlobalPenaltyFn[
        100. *(x^2 - y)^2 +  (1.  - x)^2, {}, {}, {{x, -5, 
            5}, {y, -10. , 10. }},  .0000001, Starts -> 1, 
        ShowProgress -> False] // Timing

IntervalMin[
        100. *(x^2 - y)^2 +  (1.  - x)^2, , , {Interval[{0. , 3.1}], 
          Interval[{0. , 3.2}]}, {x, y},  .0000001] // Timing

Such extremely flat functions can be slow for IntervalMin to solve.  In this case, IntervalMin did pretty well, but did not get the exact answer.  If we restart with narrower Intervals around the above solution, we get a better result:

IntervalMin[
        100. *(x^2 - y)^2 +  (1.  - x)^2, , , {Interval[{ .98, 
              1.1}], Interval[{ .97, 1.1}]}, {x, y},  .00000001] // 
      Timing

Of particular interest is how the function performs on highly irregular functions.  The following irregular function is tested:

Plot [
        Abs[2  x  +   3  Sin[x]],  {x,  -40,  40}];

GlobalSearch[
        Abs[2  x  +  3 Sin[x]], , , {{x, -40, 40}},  .000001, 
        Starts -> 10] // Timing

All 10 starting values found the solution to this problem.  In the next problem, a function is tested that is step-wise discontinuous:

Plot[Abs[IntegerPart[i]], {i, -10, 10}]

GlobalSearch[Abs[IntegerPart[x]], , , {{x, 0, 100}},  .1, 
      Starts -> 1, CompileOption -> False]

The solution is 0 between -1 and 1, so a valid solution was found.  We can also solve this using an integer solution:

MultiStartMin[Abs[x], , , {{x, 0, 100, Integer}},  .1, 
      Starts -> 1]

Problems can be solved that cause FindMinimum to fail.  In the following, any number below 3 has the minimum function value of 1.

Plot[Max[1, x - 2], {x, 0, 4}]

FindMinimum[Max[1, x - 2], {x, 5, 6}]

GlobalSearch[Max[1, x - 2], , , {{x, 5, 9}},  .1, Starts -> 1, 
      CompileOption -> False]

In the next problem, a more irregular function is tested:

Plot [
        4  Sin[x] + Sin[4  x] + Sin[8  x] + Sin[16  x] + Sin[32  x] + 
          Sin[64  x],  {x, 0, 6.28}];

GlobalSearch[
        4  Sin[x] + Sin[4  x] + Sin[8  x] + Sin[16  x] + Sin[32  x] + 
          Sin[64  x], {-x, x - 6}, , {{x, 0, 6.28}},  .000001, 
        Starts -> 15] // Timing

Eight of the 15 starts found the minimum (-6.98759) and the rest were quite close.  We see next that IntervalMin fails on this problem because the way terms are handled separately does not allow the algorithm to narrow down the range.

IntervalMin[
        4  Sin[x] + Sin[4  x] + Sin[8  x] + Sin[16  x] + Sin[32  x] + 
          Sin[64  x], {}, , {Interval[{0. , 6.28}]}, {x},  .000001] // 
      Timing

In the next problem, the scale of the fluctuations makes solution even more difficult than in the above problem.  Nevertheless, all ten starts find the solution.

Plot [Abs[2  x + x*Sin[x]],  {x, -40, 40}];

GlobalSearch[Abs[2  x + x*Sin[x]], , , {{x, -40, 40}},  .00000001, 
        Starts -> 10] // Timing

All 10 solutions are correct.  With a smaller range for the initial guess, the problem is still solved correctly for most of the starts.

GlobalSearch[Abs[2  x + x*Sin[x]], , , {{x, 36. , 39}},  .000001, 
        Starts -> 10] // Timing

However, with a small enough range for the initial guess, the problem may not be solved because the solver may be converted to a local solver.  We see next that IntervalMin actually solves this problem quickly and correctly brackets zero.

IntervalMin[
        Abs[2  x + 
            x*Sin[x]], {}, , {Interval[{-40. , 40}]}, {x},  .0000001] // 
      Timing

In the next test, a problem with multiple minima, the Branin rcos function, is solved.

b = 
        GlobalSearch[(y - x^2*5.1/(4. *π^2) + 5. *x/π - 
                  6. )^2 + 10. *(1.  - 1. /(8. *π))*Cos[x] + 
            10. , {-40 - x, x - 40, -40 - y, 
            y - 40}, {}, {{x, -40, 40}, {y, -40, 40}},  .000001, 
          Starts -> 30, CompileOption -> True];

c = {}; 
    Do[AppendTo[c, {x, y} /. b[[ii, 1]]], {ii, 1, 30}]

ListPlot[c];

All six solutions to this problem were found with 30 starts.   In the next test, the Csendes function is tested, which has a highly dissected but flat region around the optimum, shaped like an ashtray (see next section for illustrative figures).

GlobalSearch[
          x^6 * (2.0  +  Sin[1.0/x])  +  
            y^6 * (2.0  +  Sin[1.0/y]),  , , {{x,  -1. ,  
              1. },  {y,  -1. ,  1. }}, 0.0000000000000000001] // 
        Timing

We can see that the minimum point {0,0} is being approached with arbitrary accuracy.  In the test of this function with the GlobalMinima function below, we see that GlobalMinima can be used to define the flat region around the minimum.

The following function has many local minima:

Plot3D[
        20  Sin[Pi/2 (x - 2  Pi)] + 
          20  Sin[Pi/2 (y - 2  Pi)] + (x - 2  Pi)^2 + (y - 
                2  Pi)^2, {x, 0, 10}, {y, 0, 10}];

We can see that this could be difficult to fit.  It is known that the solution is -38 at {5.322,5.322}.

GlobalSearch[
        20  Sin[Pi/2 (x - 2  Pi)] + 
          20  Sin[Pi/2 (y - 2  Pi)] + (x - 2  Pi)^2 + (y - 
                2  Pi)^2, {}, , {{x, 0, 10}, {y, 0, 10}},  .000001, 
        Starts -> 4] // Timing

GlobalSearch is able to find the true minimum when several starts are used.  IntervalMin solves this problem reliably:

IntervalMin[
        20  Sin[Pi/2 (x - 2  Pi)] + 
          20  Sin[Pi/2 (y - 2  Pi)] + (x - 2  Pi)^2 + (y - 
                2  Pi)^2, , , {Interval[{0. , 10. }], 
          Interval[{0. , 10. }]}, {x, y},  .0001] // Timing

IV  THE GlobalMinima FUNCTION

IV.1  Introduction

This function represents a robust approach to global optimization.  Traditional gradient (local) approaches require the user to know how many optima are being sought for the function to be solved, and roughly where the optima are, so that a good initial guess can be made.  The user, however, rarely can make a good initial guess and has no information about the existence of multiple minima.  In the absence of such information, existing algorithms will generally converge to an answer, but this is only a local solution that may not be globally optimal.  Furthermore, even global optima may not be unique.  In addition, if the region around the optimum is very flat, most algorithms will report difficulty with convergence.  In reality, such a flat region, for real-world problems, may indicate the existence of an indifference zone of equally good solutions.

This function is designed to overcome these difficulties.  It uses the Adaptive Grid Refinement (AGR) algorithm (an n-dimensional adaptation of Newton's method related to adaptive mesh generation) which finds multiple optima and defines optimal regions in a single run.  Because the user does not need to know much about the function, the AGR algorithm saves considerable time for the user.  It is not necessary to know how many optima exist or where they are; the AGR algorithm will generally find them all.  Feasible starting points do not need to be provided.  The AGR algorithm does this at some cost in computer time by making more function calls, and is thus limited to smaller problems. The AGR algorithm also increases the user's certainty that all useful solutions have been found.  The genetic algorithm method can find multiple optima, but only stochastically, it does not allow constraints, it is difficult to use, and does not define optimal regions. Other global methods that are available are either difficult to use, require special conditions for the user function, or do not allow constraints.  These global methods also do not generally provide information on solution regions, but rather assume that only optimal points are being sought.  Compared to existing algorithms, the AGR algorithm is the easiest to use, is the most general, and provides the most information. It is not feasible for larger problems, however, in which case GlobalSearch, GlobalPenaltyFn, IntervalMin, or MultiStartMin should be used.

IV.2  The Grid Refinement Approach to Nonlinear Optimization

The AGR algorithm works by a grid refinement process.  Consider a function

(x-2)^2

where we wish to minimize z.  In this case, the answer is obviously z= 0 at x= 2, as can be seen from the plot:

Plot[(x - 2)^2,  {x,  -2,  6}];

The Adaptive Grid Refinement  (AGR) algorithm works as follows.  The interval to be searched for a solution (say, -100 to 100) is gridded in this case into n initial grid points (with n = 11 being a useful initial grid density) at a distance d = 18.18 units apart.  At each x point, z is evaluated.  The best points (with the lowest z values) are kept, and the rest are thrown away.  At each x kept, new points (daughters) are evaluated on each side at one-third the distance between the first set of points, and z is again evaluated.  This one-third trisection prevents duplication of points during daughter point generation.  This process of grid refinement continues until a user stopping criterion is met, when all optimal points are displayed.  The same procedure is used for any number of dimensions.  In two dimensions, each point generates eight daughter points.  In three dimensions, 26 points are generated as the corners of a hypercube.  The grid refinement algorithm is in essence a generalized-descent approach.  At each iteration, the population of points in the working set is moving downhill, but over multiple regions and directions.  Consider

\(min \((z)\) = ∑\+\(i = 1\)\%2 x\_i\^\(\(2\)\(\ \ \ \ \ \)\)\)

for a two-dimensional (2D) problem.  For an initial  grid that is symmetric around 0, points in all four quadrants will remain at each iteration and will converge toward 0 from all four directions, very much like four steepest-descent vectors.

	The AGR algorithm is derivative-free and uses no search vectors.  It is therefore more numerically stable than other algorithms.  For example, discontinuities (such as step functions) do not cause problems, nor do evaluations in the vicinity of a constraint boundary.  The algorithm uses very little computational time in operations other than function calls.  The size of the problem (the number of dimensions) is limited only by execution time.    The transparency of the algorithm's operation makes this tool ideal for teaching concepts of optimization to introductory college classes or high school math students.  No advanced mathematical training is required to understand and use it.
	
	The AGR algorithm is similar to a variety of interval methods that have been formally shown to converge (Pinter, 1996).  AGR also has certain similarities to Taboo Search (Cvijovic and Klinowski, 1995) and Jones' algorithm (Jones et al., 1993).
	
	In the worst case, with the optimum falling exactly halfway between two initial grid points, the error estimate of the optimum will decrease as

(1/(2*(3^i))) ((Subscript[x, max]) - (Subscript[x, min]))/n

where i is the number of iterations, the min and max terms reflect the bounds on the search region for the choice parameter x, and n is the number of initial grid points.  Thus, if the initial grid is at 10 unit intervals, after 6 iterations we can estimate the optimal parameter value to within ±0.000686.  After 10 iterations, we can achieve ±0.00000847.  On average, a point will fall away from this worst case scenario, and convergence will be slightly faster.

Grid refinement algorithms tend to increase their computational costs exponentially and can use too much memory space.  In the AGR algorithm, memory constraints are overcome by judicious memory management and by repeated pruning of the working solution set during algorithm operation.  Execution time presents more of a problem.  On a fast PC, small problems run in seconds, problems with up to 7 variables run in minutes, problems from 8 to 11 variables run in up to hours, and problems of 12 to 15 variables run overnight.  This is because millions of function calls may need to be made in higher dimensions.  The AGR algorithm is therefore not suitable for functions whose solution is itself time consuming, such as optimizing  the output of a simulation model.  On a fast Unix machine, up to 20 variables might be feasible.  Although the practical upper limit of 20 variables prevents this algorithm from being universal, its ease of use and ability to find all optima and optimal regions certainly recommend it for the large subset of problems that fall within its domain.

IV.3  Program Operation

IV.3.A  Parameters and Default Values

The user must input several parameters.  The initial grid density defines the initial gridding on each x dimension.  Thus, for a 2D problem, an initial grid density of 11 x 11 will give (11^2)  = 121 initial function calls.  As the dimensionality increases, it is advisable to reduce the initial grid density to reduce execution time.  On the other hand, the more (true or false) minima there are likely to be and the more poorly behaved the function, the larger the initial grid density should be (15 x 15, 20 x 20, or even higher).  Otherwise, minima may be missed.  A high grid density with a high dimensionality may not be feasible, since initial grid exploration requires a function call at the number of grid points raised to the power of the number of dimensions.

n = (grid^d)

The second input parameter, the contraction coefficient (C), defines the degree of contraction of the grid region on each iteration.  For the largest ((Subscript[z, l])) and smallest ((Subscript[z, s])) function values at the current iteration, a point (Subscript[z, i]) will be kept for further refinement only if

(Subscript[z, i]) < (Subscript[z, s]) + ((Subscript[z, l]) - (Subscript[z, s])) C

That is, for C = 0.1, the only points that will be kept are less than 10% of the range larger than the current minimum. If local minima are expected to be a problem, values of 0.25 to 0.5+ will help prevent the true minima from being missed.  The parameter C is cut in half twice as the program proceeds.

The tolerance (T) defines the amount of improvement (absolute, not relative) in the smallest z required at each iteration.  If at least this much improvement is not found, the program stops.  The tolerance can be set to very small values, such as (10^(-15)), to achieve a highly accurate estimate of the minimum.  If the user sets tolerance = 0, the program will stop if the minimum remains the same for two iterations, but if there is a very gradual slope near the minimum and slight improvement is achieved with each iteration, the program will not terminate.  Thus, tolerance values of zero are not advisable.  To select a good tolerance value, the user can make an initial run with a large tolerance to see if multiple optima exist, and then make a second run with a smaller tolerance and bounds that are closer around the rough estimate of the optimum (or around each solution region if there are several).

A second solution criterion is provided by the indifference (I).  When the tolerance criterion causes the program to stop, all points from intermediate iterations with values between (Subscript[z, s]) and (Subscript[z, s]) + I are considered part of the optimal region(s).  This is because the user is often indifferent to extremely small differences in z, but may be interested in the existence of multiple solutions that are (within I) equally good.  Thus, for a steel mill, there may be a range of furnace temperatures within which the output steel quality is approximately the same.  There is no point in trying to control the furnace temperature closer than this target zone; in fact, such control is likely to be costly or impossible.  In general, for real-world optimization problems the existence of zones of indifference is extremely important information, because it is far easier to control a process within some range than to keep it exactly at some "optimal" point.  Only the AGR algorithm provides information on zones of indifference to guide real-world decision making.  Consider the function

z = (x^4)

with I = 0.1, which is minimized at x = 0.  All values of x between -0.562 and +0.562 produce a z less than I, and the user is thus indifferent to these differences.  The "answer" is therefore the range -0.562 < x < 0.562 and not just the point x = 0.  The AGR algorithm provides an approximation to this range that improves as initial grid density improves  If only the absolute minimum point is being sought, I can be set to zero.  If I is set very large, the solution set may define a large region containing many thousands of points.

For higher-dimensional problems, initial grid density may be so coarse that a full picture of the indifference region is not found.  This situation will be evident if there are few points in the indifference region relative to the dimensionality.  The AGR algorithm can be reapplied to a smaller domain around a particular optimum by setting the ranges for (Subscript[x, i]) closer to this zone than in the initial run in order to better define the indifference region.

The AGR algorithm is implemented in Mathematica via the GlobalMinima command, which has eight options.  The defaults are MaxIterations->25, CompileOption->True, ShowProgress->True, StartsList->{}, EvaluateObj->True, SimplifyOption->True, UserMemory->32, and Starts->3.  MaxIterations prevents the program from running away to infinity when a mistake is made in the function (Min[x] for example).  For very computationally expensive problems (like those with many parameters), it is useful to use Starts->1 to evaluate the time a run takes.  For general usage, Starts->3 (the default) is good. The results of this short run can then be used to define further, more directed runs based on whether all 3 starts found the same solution.  A list of starting values can be input instead of letting the program find them with random search (e.g. StartsList->{{1.,2.},{2.,1.}}).  The CompileOption determines whether the user function and constraints will be compiled.  While Compile reduces execution time, some functions can not be compiled, in which case CompileOption->False should be used.  If the objective function should not be Evaluated, use EvaluateObj->False.  SimplifyOption attempts to simplify the objective function.  If this should not be done or will not improve the solution, this should be set to False.  UserMemory, in units of Megabytes, is a user input defining the free memory available on their machine.  GlobalMinima has several ways of executing, depending on how much memory is available.  If more memory is available, it runs faster.

IV.3.B  Bounds

Whereas local algorithms require input of an initial guess or starting point (which must be feasible) to initiate the search, GlobalMinima requires only bounds on the variables.  Bounds are generally easier to define than are feasible starting points.  Physical, logical, or economic considerations often provide guidance on realistic bounds for a variable.  If the optimization problem is wing design and one variable is wing thickness, a certain range of thicknesses is sensible to consider.  For a financial problem, the range of values for an investment mix across a product line is clearly bounded between zero and the total budget, and the upper limit can generally be bounded below a much smaller value than this.  If narrower bounds are known, then convergence will require fewer function calls, but even quite wide bounds are acceptable.  GlobalMinima can not search outside of the bounds, and if the solution found runs right against an input bound, then the true solution may lie outside the bounds and another run may be required if the bounds do not actually represent a constraint.

IV.3.C  Constraints

Constraints are entered as Mathematica functions, which may be linear or nonlinear.  Constraints may be any logical function that defines a region or regions of space.   Constraints are entered as a list.  The terms of the inequalities can be any valid Mathematica inequality expression.  For example, we can define a 3D, doughnut-shaped search region (a sphere with a hollow middle) as follows:

1.  < x^2 + y^2 + z^2 < 10. ;

In contrast to gradient search algorithms, for which constraints are a real problem, GlobalMinima handles constraints easily.  In fact, constraints, by limiting the space to be searched, may greatly reduce the number of function calls.  When constraints are involved, it is advisable to increase the initial grid density somewhat.  Any number of constraints is allowed, but equality constraints are not allowed.   Equality constraints can be approximated by a bounding constraint such as:

0.9 < x1 + x2 + x3 < 1.1;

which becomes in standard form:

{x1 + x2 + x3 - 1.1, -(x1 + x2 + x3) + 0.9};

but in this case, a higher initial grid density is needed to obtain points that fall in the narrow feasible zone.

	Note that simple bounds on the parameters are handled by the variable definition terms {x, (Subscript[x, min]), (Subscript[x, max])}.  Duplicating these bounds as constraints is valid but is unnecessary and will slow program execution.  To be specific, no nonnegativity constraints are needed, as these are handled by the bounds inputs.

WARNING:  Constraints can cause difficulties for GlobalMinima.  Problems are caused by the creation of unreachable regions of the search space by the constraints.  In an unconstrained problem, the gridding is such that the entire search space can be reached by iterating from the initial grid configuration.  This may not be true when constraints are present.  This can happen in several ways.  First, the feasible space may be defined such that no initial grid points fall in the feasible region.  For example, for bounds {-100,100} with a gridding of 5, a constraint {.1<x<.2} will be such a small space that the initial grid won't fall into it except with a very high initial gridding.  If this happens, GlobalMinima will exit with an error message.  Second, a small unreachable space can be created right against the constraint boundary.  For example, for initial grid points {-5,0,5}, the point 0 provides access to the region -2.5<x<2.5.  If a constraint is entered {x>1}, the region {1>x>2.5} is no longer reachable, since only the grid point {5} is available for grid refinement.  Thus in this case, the smallest value of x that can be found is 2.5 if the minimization tends to push the solution up against the constraint.  The solution to this difficulty is to rerun the problem with tighter bounds focused around the solution found, with a higher grid density and/or a shift in the grid.  For example a shift from 4 to 5 grid points will alter the intersection between the gridding and the constraint region, as will a zoom in or an alteration in the left or right bounds.

IV.3.D  Memory Usage

The algorithm used by the program is able to trade memory for speed.  This means that if more memory is available it will run faster.  The default for the option UserMemory is 32 (in units of Megabytes).  If your machine has more memory than this, setting UserMemory to the higher value will likely make it run faster.  If the kernel runs out of memory, try quitting all other applications.  If the kernel still runs out of memory, it may be necessary to restrict UserMemory to smaller values to force GlobalMinima to use the computational path that requires less memory but takes longer.

IV.4  Output

During program execution, output showing the size of the current solution set and the current minimum is optionally printed at each iteration.  ShowProgress->True will print out these intermediate results, which provide information for evaluating convergence and whether the contraction coefficient is too tight or too loose.  Especially for higher dimensional problems, this intermediate output should always be printed.  Following solution, a list is returned that contains the solution.  For example, two points in a 2 parameter problem would give the list {{{1.1,1.2},-5},{{1.2,1.21},-4.9}}.  The total number of function calls is also printed.  Output is illustrated below.

IV.5  Maximization

The algorithm assumes that z is being minimized.  If a function is to be maximized, simply minimize -z.  The results will then be in terms of the correct parameter values, but the z values will be negatives of their true values.

IV.6  Limitations

No optimization algorithm can solve all conceivable problems.  Consider a response surface in two dimensions that is flat everywhere except for a very small but very deep pit at one location.  Steepest-descent algorithms will fail completely to solve this problem, because for almost all starting points there is no gradient.  The algorithms will just stop and say that the gradient is zero.  GlobalMinima will stop after two iterations because no reduction has been made in the minimum, and it will identify points across the entire region as valid.  This result tells the user that the surface is very flat.  The next step that can be taken, if there is reason to believe that an optimum does exist somewhere, is to increase the grid density and run GlobalMinima again.  Initial densities > 1000 are feasible for problems with few dimensions and may by chance find the narrow optimal zone.  Once this zone is identified, convergence is rapid.  There is thus at least a chance that GlobalMinima can solve such an "impossible" problem.

Some functions have infinitely many solutions.  Consider Sin(1/x) between 0 and 1.  As the function approaches zero, the minima (with z = -1) get closer and closer together, as we can see from the plot following:

Plot [Sin[1/x],  {x,  0,  1}];

Although it is fundamentally impossible to list all solutions when there are infinitely many, GlobalMinima can show that there are very many solutions.  If we try an initial grid density of 20 on the interval [0, 1], with contraction = 0.4, we obtain many points with values close to -1.0 (see section IV for details).  This result demonstrates that there are many solutions and that the density of solutions increases with an approach to zero.  This is thus an approximate solution to the true situation.  One can determine that this is an approximate solution by increasing the grid density to find more solutions close to zero.  Making two runs allows us to extrapolate the result that there are very many solutions close to zero, perhaps infinitely many (as in this case).

These examples illustrate that not all problems can be solved.  Some are not solvable by any algorithm.  Some may be ill-posed.  However, for many problems that cause other algorithms to fail, the AGR algorithm either succeeds, has a probability of succeeding, or provides a partial or approximate solution. We may thus say that the set of unsolvable problems is much smaller for the AGR algorithm than for other solution methods.  This is particularly so with respect to the inclusion of constraints.

IV.7  Efficiency

Although the number of function calls is a typical criterion for comparing the speed of optimization algorithms, this measure does not tell the entire story.  Many algorithms do considerable computing for things such as derivatives and gradients and inverting matrices, whereas the AGR algorithm does very little such secondary computing.  The AGR algorithm is thus faster than the number of function calls alone would indicate.  A second important point is that the literature comparing algorithms often indicates that an algorithm had only some probability of finding a solution (for a random initial starting point).  For real problems where the “correct” answer is not known, this means that the problem must be solved many times with different initial conditions to start the algorithm; sometimes hundreds or thousands of runs are necessary.  Genetic algorithms and simulated annealing behave in this way (Ingber and Rosen, 1992).  Thus, the actual efficiency is a function of the total number of runs, not just the single-run function call count or the execution time.  This requirement for making multiple runs creates substantial work for the user.  With GlobalMinima, a single run is sufficient to generate reliable solutions.  In addition, it is easier to define feasible bounds on the domain of a problem than to generate hundreds of feasible starting points for multiple runs, especially if constraints are involved.

It is obvious that a grid refinement approach like that used here will require function calls as an exponential function of the number of dimensions.  The initial gridding can be slow if too great a grid density is input.  For example, for a five-dimensional problem with a grid density of 7, the initial number of grid points is (7^5) = 16, 807.   After the initial grid calculation, execution time is governed by the generation of daughter points in the process of grid refinement.  The formula for the number of daughters d generated by a point in the current solution set is:

d = (3^n) - 1

This gives d = {2, 8, 26, 80, 242 ...} for i = {1, 2, 3, 4, 5 ...}.  With this formulation, large problems take too long to run.  One approach to reducing execution time is to make the contraction coefficient smaller.  For problems with only one minimum, this economy can cut execution time by 0.5 to 0.9.  However, this strategy will cause multiple optima to be missed.  The grid-based approach used here is thus obviously not suitable for very large problems, but does allow for constraints, even nonlinear constraints.  Therefore, a contraction coefficient much below 0.05 is not recommended unless the function is likely to be smooth.  For large problems or problems with complex constraints, the MultiStartMin function is recommended.

IV.8  Error Messages

An effort has been made to trap as many errors as possible.  Much of the setup for running GlobalMinima uses Mathematica formats.  Many user errors will prompt Mathematica error messages.  On input, only the order of parameters is evaluated in the passing function call, not their names.  The user is thus free to use any names, though the names used here and in the example notebooks are the most obvious and will help prevent mistakes.  Therefore using “c” for “contraction” is not an error.  If valid values for input parameters (tolerance >= 0.0, indifference >= 0, 0.0 < contraction < 0.8) are not input, the following error message is printed, and the program stops:

	Error found by GlobalMinima
	Error:  Input parameter outside range

If a parameter is omitted from the calling sequence, Mathematica cannot continue.  To indicate the problem, it echoes back the names for functions and values for constants in the function call, and then it stops.  This type of echoed-back output means that a format error exists in the calling sequence.

After input, the program compiles the user input function and the constraints, speeding execution by a factor of about six.  If these functions are improperly formatted, this step may produce a compilation error message, which will terminate the run.  If a function is defined such that it does not return numeric values, this will cause error messages from Mathematica, followed by termination of GlobalMinima.

If a completely flat function is mistakenly defined, every point that is tried will be kept, and the solution set will keep growing.  The program traps this type of error and exits with a message that the response surface is completely flat after the initial gridding.  If a flat function is input for a large problem with a high initial grid, all available memory could be used up before GlobalMinima can detect that the function is flat.

When constraints are used, there may be no feasible region, particularly if the problem has many constraints.  If this occurs, the following error message is printed, and the program terminates:

	Error found by GlobalMinima
	Error:  No valid grid points generated

In this case, check the constraints and their relation to the bounds on the parameters.  Constraints and bounds can be in conflict if care is not taken.  In some cases, the constraints are valid, but a higher initial grid density is needed to obtain feasible points.

A common mistake is to maximize when the intention was to minimize.  This error can be recognized because the solution will tend to run into one of the parameter bounds and will have an illogical z value.  For example, the optimal airplane wing thickness for a design problem might come back as zero.

An unbounded solution may result from an improperly formulated problem (e.g., min(-1/x) over {-1, 1} which becomes -inf at x=0).  Because the function z continues to get larger the closer one gets to zero at an increasing rate, the program will never report convergence.  If this mistake has been made, GlobalMinima will terminate after MaxIterations and will print out results showing a very large (negative) function value.

The user function passed in to GlobalMinima is Compiled.  This may cause compilation errors.  The result is usually to stop the execution of the program, but it is necessary for the user to realize the source of the error.  Compiled functions can also generate errors at run time if the user function generates non-machine numbers such as high precision numbers or infinity.  Mathematica will usually revert to the uncompiled function and continue running.

If the user defines variables in his program by mistake that should be parameters of the function, this will cause the program to stop or malfunction.  For example, if pp=5 is defined by mistake, this will cause Mathematica to detect an error but GlobalMinima to attempt to execute:

pp = 5;

\(\(\(GlobalMinima[\((pp - 30)\)^2\ \((y - 60)\)^2, , {{pp, 0. , 
          100. }, {y, 1, 5}}, 20, 0.001, 0.3, 0.01]\)\(
\)\)\)

In this example, the parameter list {pp,0.,100.} has become {5,0.,100.} which does not match the function input sequence, causing the program to stop.

IV.9  A STEP-BY-STEP EXAMPLE

Our example is a function with multiple optima, the camelback function (Griewank,  19__):

Plot3D[
        4*x^2 - 2.1*x^4 + x^6/3 + x*y - 4*y^2 + 4*y^4,  {x,  -1,  
          1}, {y, -1, 1}];

This function has two minima, which can be seen in the figure.  GlobalMinima easily finds the two solutions. Even with an initial grid = 7 points and contraction = 0.3, GlobalMinima finds min = -1.293439 at the two solution points {0.4867, -0.74074} and {-0.4867, 0.74074}.  We now illustrate how this solution is obtained.

GlobalMinima is a Mathematica package and must be installed by using the Get command.  The function to be minimized must be defined; in this case, it is the camelback function depicted above.  This is the first parameter for the function call.  The second parameter is the constraints, input as a list of inequality functions separated by commas.  If there are no constraints, an empty list must be entered as {}.  When a set of constraints is entered, Mathematica returns a True or False value.  GolobalMinima tests each potential grid point against the constraint equation set.  The third parameter is the variable names and their bounds, in any order, each defined as a list.  Each variable used in the function must be defined.  The fourth parameter, the indifference parameter, defines the range of values larger than the best solution found that we will consider essentially equivalent and therefore equally good. We must define initial grid density (Integer) and tolerance (Real), where tolerance specifies the improvement in the solution necessary to continue iterating (i.e., if (old min - min) > tolerance, keep going).  Note that this is an absolute rather than a relative tolerance.

	The contraction coefficient is defined for input where the range is typically 0.5 to 0.1. The contraction coefficient is cut in half after the first iteration and in half again after the second.  This is because the initial grid narrows the solution space to regions where faster winnowing of points is possible.

We can now execute GlobalMinima.  Note that parameter names are arbitrary but that order is fixed.  The option ShowProgress causes intermediate output to be printed.  In this first example, no constraints are being used, so an empty list {} is entered.  By setting indifference very small, we only get the best solutions:

toler = 0.001;
grid = 7;
contraction = 0.25;
indif = 0.001;

GlobalMinima[
      4. *a1^2*a1^4 + (1. /3. )*a1^6 + a1*a2 - 4. *a2^2 + 
        4. *a2^4, , {{a1, -1. , 1. }, {a2, -1, 1. }}, grid, toler, 
      contraction, indif, ShowProgress -> True]

GlobalMinima in this example provided stepwise reports on its progress.  After echoing back the input parameters, it reports on the memory needed for computations.  If more memory is available, the program will run faster unless it is already printing out "full efficiency".  At each iteration, it reported the current minimum and maximum and the number of points in the solution set.  For large problems, the execution time can become extremely high if too many points are being kept.  In this case, it can be useful to make the contraction coefficient somewhat smaller.  If a problem is known to be convex and well-behaved (e.g., z = (x^2)), a very small contraction (<0.05) will increase speed dramatically.  A helpful type of printed output is a 2D plot of indifference-region values, for cases with two variables only.  These plots help provide insight into program operation and the process of grid refinement and convergence.
	
	The program can be run with more compact output by omitting the ShowProgress option, which then allows the default False to take precendence.

The camelback problem has two optimal solutions.  A constraint that selects the lower right quadrant as a valid region finds just one of the two solutions, as shown below.

toler = 0.00001;
grid = 7;
contraction = 0.3;
indif = 0.0001;

GlobalMinima[
        4. *a1^2*a1^4 + (1. /3. )*a1^6 + a1*a2 - 4. *a2^2 + 
          4. *a2^4, {-a1, a2}, {{a1, -1. , 1. }, {a2, -1, 1. }}, 
        grid, toler, contraction, indif, ShowProgress -> False] // 
      Timing

We can see that only the solutions in one quadrant were found, as desired.  We can impose nonlinear constraints as well, such as excluding a spherical search region around the origin:

(y1^2) + (y2^2) > 2. ;

When we minimize (y1^2) + (y2^2) with this constraint, we obtain the following solution space (depending on parameters):

GlobalMinima[
        x^2 + y^2, {-(x^2 + y^2) + 2}, {{x, -5, 5. }, {y, -5, 
            5. }}, 30, 0.001, 0.5, 0.0001, 
        ShowProgress -> True];

This ability to represent any number of spaces of any shape with the constraint function is one of the unique and powerful features of the AGR algorithm.

IV.10  Testing

No numerical algorithm can be guaranteed not to fail.  However, the approach taken here is far more robust than other approaches.  The algorithm has been tested with a number of both simple and complex functions and has been shown to perform well for a number of test functions that exhibit false minima, multiple minima, and minimal regions.  Test functions have included positive polynomials with multiple roots, poorly behaved functions with large numbers of false minima, the Rosenbrock function, a number of multiple-minima problems, and many others.  The only difficulty encountered is that for multiple-minima problems, some solutions could be missed if the initial grid density is too low.  The algorithm should be stable for achieving high-precision results, because no derivatives are used, and testing has shown that high accuracy can be achieved easily.  Mathematica uses high precision computations, and stores numbers as exact ratios of whole numbers whenever possible.  Thus very little error is expected from truncation.  Nondifferentiable functions, such as step functions, were tested and caused no numerical difficulties.  The algorithm has no trouble with problems having up to dozens of true optima.  For pathological functions with thousands of optima, the program may be able to find only scores of solutions in a particular run.

This package is not a tool for pure mathematics but rather is designed for applied problems.  In the Csendes function discussed below, a box around the origin contains only solutions less than some value.  Near the origin, values within the box may be (10^(-33)) or less.  GlobalMinima is designed to define this box and to converge near the best minimum.  Convergence to the absolute minimum (in this case, zero at x = 0) is of less interest, because all locations within the defined box region are essentially equally good (to (10^(-33))).  It is thus difficult to compare the AGR algorithm directly to other algorithms, because other approaches seek out the analytical optimum point whereas the AGR algorithm seeks out the practically optimal region(s).  What we can say is that although the AGR algorithm may not always find the analytical minimum of a function, it usually gets very close, and it can also find multiple minima and optimal regions, which other algorithms cannot.

IV.10.A  A Simple Polynomial.

The use of the package is best illustrated with examples.  We begin with some simple cases.  The function (∑_(i = 1)^n) (x^2) is symmetric and for any number of dimensions has a global minimum at the origin with z = 0.  This problem meets the classic criteria of gradient descent algorithms, which converge to the solution from any starting point for functions with no saddles, discontinuities, or false minima.  GlobalMinima converges to the correct answer easily.  Although the AGR algorithm is not as fast as a steepest-descent algorithm, it is quite fast for such well-behaved and steep problems.
	
	A simple modification to this function creates a problem with multiple solutions:

Plot [(x  -  30)^2 (x  -  60)^2,  {x,  0,  
          100}];

Here, z = 0 at x = 30 and at x = 60.  When we run GlobalMinima over the interval {0, 100}, we need an initial grid density of 20 points to find both solutions:

GlobalMinima[(x - 30)^2 (x - 60)^2, , {{x, 0. , 100. }}, 20, 
        0.00001, 0.3, 0.0001] // Timing

Note that for small problems like this, the solution is very fast.  While some of the values look like duplicates, note that they must not be because otherwise the values of the function would be identical.  This results from significant digits that Mathematica does not show on output.  Traditional methods would find one or the other of these solutions, depending on the initial guess given to the algorithm.  In this case, the maximum and minimum coordinates of the indifference region overstate the size of the region, which is really just two widely separated clusters of points.

IV.10.B  The Rosenbrock Function

An interesting problem is the Rosenbrock function (Ingber and Rosen, 1992)  a classic test  of convergence of optimization algorithms.  The plot of this function is:

Plot3D[
        100. *(x^2 - y)^2 + (1 - x)^2, {x, -5, 5}, {y, -10, 
          10}];

This function becomes quite shallow as the optimal region is approached, and some algorithms experience convergence difficulties with it.  We can see that GlobalMinima does not experience difficulties with it and finds the global minimum z = 0 at {1, 1}, even with an initial grid = 5 and contraction = 0.2, so long as the tolerance is small (< (10^(-7))).  The execution time is faster than when using GlobalSearch.

GlobalMinima[
        100. *(x^2 - y)^2 +  (1.  - x)^2, , {{x, -5, 
            5. }, {y, -10. , 10. }}, 8, 0.00000001, 0.2, 0.000000000001] // 
      Timing

IV.10.C  The Branin rcos Function

A related function, the Branin rcos function (Dixon and Szego, 1978), has a more complex surface:

Plot3D[(y - x^2*5.1/(4. *3.14159^2) + 5. *x/3.14159 - 
                  6. )^2 + 10. *(1.  - 1. /(8. *3.14159))*Cos[x] + 
            10. , {x, -40, 40}, {y, -40, 40}];

In the region {-40, 40}, there are six equivalent minima.  With traditional algorithms, it would be very difficult to know what initial guesses to use.  The number of possible minima is certainly not clear from inspection of the equation or the plot.  Zooming in a little on the flat region, we can see a horseshoe-shaped valley:

Plot3D[(y - x^2*5.1/(4. *3.14159^2) + 5. *x/3.14159 - 
                  6. )^2 + 10. *(1.  - 1. /(8. *3.14159))*Cos[x] + 
            10. , {x, -20, 40}, {y, -40, 40}];

However, we still do not have much information about the solutions.  GlobalMinima provides a solution:

ans = 
          GlobalMinima[(y - x^2*5.1/(4. *3.14159^2) + 5. *x/3.14159 - 
                    6. )^2 + 
              10. *(1.  - 1. /(8. *3.14159))*Cos[x] + 
              10. , , {{x, -40, 40}, {y, -40, 40}}, 11, 0.001, 0.3, 
            0.01, ShowProgress -> True];

Here, the semicolon was used to suppress output of the list of points, which are instead stored in "ans".  From this example, we can see that GlobalMinima finds all six optima with the given parameters and that the indifference parameter defines regions around these optima.

IV.10.D  The Csendes Function

A very difficult problem is the Csendes function (Csendes, 1985; Courrieu, 1997):

(Subscript[C, n]) (S) = (∑_(i = 1)^n) ((Subscript[x, i])^6) (2 + 
                    Sin  (1/(Subscript[x, i])))        - 1 < (Subscript[x, i]) < 1

This function has a unique global minimum (0) at x = 0, an oscillation frequency that approaches infinity as one approaches the origin, and a countable infinity of local minima, making local searches impractical.  For a 2D version, we can plot the function:

Plot3D[
        x^6*(2.0 + Sin[1.0/x]) + 
          y^6*(2.0 + Sin[1.0/y]), {x, -1. , 1. }, {y, -1. , 
          1. }];

At this scale, the plot appears smooth, but in cross section and closer to the origin it gets increasingly wavy as we approach zero:

Plot [
        x^6 * (2.   +  
              Sin[1. /x]),  {x,  -  .1,   .1}];

Although the waviness increases close to zero, the absolute magnitude of the oscillations goes to zero as we approach the origin.  We can try to solve this problem as follows, with grid = 10 and tolerance = 0.0000001.

GlobalMinima[
        x^6 * (2.0  +  Sin[1.0/x])  +  
          y^6 * (2.0  +  Sin[1.0/y]),  , {{x,  -1. ,  
            1. },  {y,  -1. ,  1. }}, 10, 0.00001, 0.2, 0.0000001]

The solution to this is not shown here because many pages of printout result from the identification of the indifference region.  It was found that the best solution is (10^(-15)) at {0.0037, 0.0037}, but that the entire region within the box ±0.056 in both dimensions is less than 0.00001 (based on the indifference zone criterion) and that the region within the box -0.167 to 0.21 in both dimensions is less than 0.0001.  Zooming in around the best solution and restarting with a new grid does yield some improvement (z = (10^(-30)) at {9e-06, 9e-06}), but we have to consider whether improvement from (10^(-15)) could possibly have significance.

	We can also turn the problem around and ask how bad the local solution can be in the optimal neighborhood.  In the region within ±0.056, we find that the worst solution found is around 0.0000001.  This enables us to put bounds on the set of points in the optimal region.  If these bounds are wide, then very large Z values might be adjacent to very small ones, in which case it is not advisable to view the region full of optimal solutions as truly optimal, because the region is also full of bad solutions.  From a purely mathematical viewpoint, this is irrelevant; however, from a practical viewpoint, we should not try to implement an optimal policy or design if very bad solutions are very close to it, because our ability to implement a solution is always imperfect.  An obvious example of this difficulty occurs with sin(1/x), whose plot looks like

Plot[Sin [1/x],  {x,  0,  1}];

GlobalMinima with 20 initial grid points yields z = -0.999 at points {0.0033,0.0335,0.0424,0.058,0.09}.  This is not a complete set, but it does show that many equivalent solutions exist.  Solving for the maxima on this same interval, we find that z = -0.999 at {0.022,0.025,0.03,0.037,0.049,0.07}.  We can see that maxima alternate with and are very close to minima over the region.  Unless our implementation is perfect, which is quite unlikely, we cannot in the real world implement the optimal solution for such a problem without great risk of actually implementing a worst possible solution.  Thus, when multiple optima exist, a combined minimization-maximization is useful to define the feasibility of the solutions.

IV.10.E  Wavy Functions, the W Function

Another "difficult" function is the W function (Courrieu,1997):

Wn, 
      k (x) = (1/n) (∑_(i = 1)^n1)-cos ((Subscript[kx, i])) 
              exp (-((Subscript[x, i])^2)/2)

for -π<x<π.  W functions have their unique global minimum (0) at x = 0.  The number of local minima in the search domain is kn (for k odd) or (k + 1)n (for k even).  Here we used only k = 10, which gave 121 local minima for n = 2 and more than 2.59 x (10^10) local minima for n = 10.  The function oscillates between two hulls of constant mean (= 1) whose distance from each other is maximal in the neighborhood of the solution.

	A one-dimensional slice across the W function shows why this is a difficult function for optimization algorithms.

Plot[
          1  -  Cos[10.  *x] Exp[-x^2/2],  {x,  - 3.14159,  
            3.14159}];

The solution of this function with the AGR algorithm depends on having sufficient grid density so that an initial point falls into the inner oscillation that goes to zero.  With an initial grid density of 40 points for a 2D problem, GlobalMinima finds the solution to arbitrary accuracy, depending on Tolerance.  Again for such a problem in real-world applications, it is important to know whether bad solutions are close to good ones.  Here we can zoom in on the best solution and, over the interval {-1, 1} on each dimension, find min(z) = -W.  When we do this, GlobalMinima finds the worst solution z = 1.95 at {-0.311, -0.311} and three other symmetrically located points, all very close to the global minimum (as can also be seen in the one dimensional figure above).  This is a potentially unstable situation.  We must then evaluate whether we can implement a solution to this problem in the real world with sufficient accuracy to ensure that neither noise nor error will lead to the adverse solutions rather than the desired optimal one.

IV.10.F  More Wavy Functions

Some further examples of wavy functions are instructive.  First, we consider the function Sin_1:

z =  |x + 3 Sin(x) |

Plot[Abs[x  +  3 Sin[x]],  {x,  -40,  40}];

When solved by solvers such as IMSL or the Mathematica  FindMinimum function, this function yields local minima unless the initial guess is near the origin.  GlobalMinima solved this to arbitrary accuracy (on {-40, 40}) with 26 initial grid points.  An accuracy of (10^(-12)) was achieved with 150 function calls.  

	The related waxy function Sin_2:

z = | 2x + xSinx |

Plot [
          Abs [2  x  +  x*Sin [x]],  {x,  -40,  40}];

with 12 initial grid points GlobalMinima attained an accuracy of (10^(-12)) with 152 function calls.

IV.10.G  Wavy Fractals

As an example of a very wavy function, a fractal (nondifferentiable) figure on the line was constructed:

Plot[
          4  Sin[x]  +  Sin[4  x]  +  Sin[8  x]  +  Sin[16  x]  +  
            Sin[32  x]  +  Sin[64  x],  {x,  0,  6.28}];

(Here this figure was iterated to just six levels and evaluated between 0 and 2π.) With just 10 initial grid points, GlobalMinima found the solution z = -6.98759 at x = 4.48989.  Even with as few as 5 initial grid points, the optimal solution was found.

GlobalMinima[
          4  Sin[x]  +  Sin[4  x]  +  Sin[8  x]  +  Sin[16  x]  +  
            Sin[32  x]  +  Sin[64  x],  , {{x,  0. , 6.28}},  10, 
          0.00001, 0.2, 0.0000001]//Timing

Note that this is much faster than using MultiStartMin or GlobalSearch and also more reliable on this very wavy function.  GlobalMinima will be faster for small problems (<3 variables without flat regions).

V  NONLINEAR REGRESSION: THE NLRegression FUNCTION

V.1  Introduction

NLRegression is a function that performs nonlinear regression using nonlinear least-squares.  It is able to incorporate constraints which may represent physical limits on the parameters being estimated, such as that a growth rate can not be negative.

The optimizer underlying NLRegression is GlobalSearch.  All options used by GlobalSearch are accepted by NLRegression and passed down to it except StartsList.  See the GlobalSearch section for its options.  The only other option accepted is Weights, which is the weight to be assigned to each data point.  Weighting might be employed when some of the data have a greater reliability than others.  The format for the function call is

NLRegression[data,expression,independent_variables,inequalities,equalities,{{var1name,lowbound,highbound}..},tolerance,options]

The output of the NLRegression procedure is a table of confidence intervals on the parameters, sensitivity plots, the fit parameter values, and the fit statistics.  The sensitivity plot can be very useful for detecting parameters that are redundant.

The typical regression problem involves fitting data to a nonlinear model.  If the model is nonanalytic (black box), the problem can not be solved by NLRegression, but instead must be solved by setting up the least-squares error function which is then passed in to GlobalSearch.  If the function returns values that are below machine accuracy or can return complex results, CompileOption->False should be used.  There are four options: FitStatistic, SensitivityPlots, Weights and Norm:

?FitStatistic

?Weights

?Norm

?SensitivityPlots

?UserResiduals

Let's consider a simple regression problem.

p = Table[i, {i, 0, 20}]

ClearAll[a, b, x]

eqn = a + b*x^2

Where a and b are parameters, and x is the independent variable.

data = {}; 
    Do[AppendTo[data, {x, eqn /. {a -> 1, b ->  .1}}], {x, 0, 
        20}]

data

ListPlot[data]

r = 
      NLRegression[data, eqn, {x}, {}, {}, {{a, 0, 3}, {b, 0, 1}},  .0000001, 
          SensitivityPlots -> False, StartsList -> {{2,  .5}}] // 
        Timing

We recover exactly our original parameters.  It is possible that the best fit to data violates some constraint such as physical feasibility.  In this case, it is necessary to constrain the optimization problem.  We introduce some measurement error into the data.

data[[1, 2]] = data[[1, 2]] - 1

data[[2, 2]] = data[[2, 2]] - 1

data

ListPlot[data]

r = 
      NLRegression[data, 
          eqn, {x}, {}, {}, {{a, 0, 3}, {b, 0, 1}},  .0000001] // 
        Timing

We see above that the error added to the data causes "a" to fall down to .792.  If we know that parameter a can not fall below .9 for some reason, we may add this constraint to the problem.  a>.9 is converted to standard form as {-a+.9}.  We add this to the function call.

r = 
      NLRegression[data, 
          eqn, {x}, {-a +  .9}, {}, {{a, 0, 3}, {b, 0, 1}},  .0000001] // 
        Timing

The addition of the constraint satisfies the restrictions on the problem. Note that with constraints the confidence intervals and sensitivity plots are not defined and can not be printed. Another way to approach this problem is to use the L1 norm, which reduces the influence of outliers.  We try this next, without constraints.

r = 
      NLRegression[data, eqn, {x}, {}, {}, {{a, 0, 3}, {b, 0, 1}},  .0000001, 
          Norm -> L1] // Timing

We note two things immediately.  First, the sensitivity plots are now linear rather than quadratic, because of the different norm.  Second, the L1 norm recovers the true parameters in spite of the error in the data and without the use of constraints.

V.2  Utilizing Chi-Square fit Criteria

In the next example,the fitting of a function to data is illustrated when the data span a very large range of values.  This causes nonlinear estimation to have some difficulties.  Under these circumstances, it is useful to use a chi-square fit statistic:

chisq = 
      Sum[((model - data[[i]])/data[[i]])^2, {i, 1, 
          datalength}]

The problem is illustrated with the same model as used in the previous section, but with a very wide range of values for the independent variable.

<< Graphics`Graphics`;

p = Table[i, {i, 1, 20}]

ClearAll[x, a, b, c, g]

y1 = d*(a*x^g + b) + c*(2*b + a*x^g)

c = 1. ; d = 0. ;
dat1 = {}; curve1 = {}; curve2 = {}; 
    Do[k = y1 /. {x -> p[[i]], a ->  .1, b -> 10. , 
            g -> 2}; AppendTo[dat1, {c, d, p[[i]], k}]; 
      AppendTo[curve1, {p[[i]], k}];, {i, 1, 20}];
c = 0. ; d = 1. ;
Do[
        k = y1 /. {x -> 1000*p[[i]], a ->  .1, b -> 10. , 
              g -> 2}; AppendTo[dat1, {c, d, 1000*p[[i]], k}]; 
        AppendTo[curve2, {1000*p[[i]], k}];, {i, 1, 20}];

p1 = 
      LogLogListPlot[curve1, 
        PlotStyle -> {RGBColor[1, 0, 0], PointSize[ .02]}, 
        GridLines -> Automatic, DisplayFunction -> Identity]; 
    p2 = LogLogListPlot[curve2, 
        PlotStyle -> {RGBColor[0, 0, 1], PointSize[ .02]}, 
        GridLines -> Automatic, DisplayFunction -> Identity]; 
    Show[p1, p2, DisplayFunction -> $DisplayFunction, 
      PlotRange -> All, Frame -> True]

ClearAll[c, d, x]; indepvars = {c, d, x}

If we solve this problem with the default FitStatistic->LeastSquare (the default), the parameter b deviates far from it's true value, and convergence takes a long time.  Even if it is constrained, convergence is very slow.  However, if we use the ChiSquare option, convergence is very fast and the constraints on b are not needed.

r = 
      NLRegression[dat1, 
          y1, {c, d, 
            x}, {}, {}, {{a,  .09,  .12}, {b, 9, 11}, {g, 1.9, 
              2.1}},  .0000001, ShowProgress -> False, 
          CompileOption -> True, FitStatistic -> ChiSquare, 
          SensitivityPlots -> False] // Timing

We see that we have recovered the true parameters {.1,10,2} quite well.

V.3  Multiple Independent Variables in Regression Problems

In the next example,the fitting of a function to data is illustrated when there are multiple independent variables.  A common example of such a problem is when several experiments are done under different conditions and all data from all experiments are to be used to estimate certain parameters.

p = Table[i, {i, 1, 20}]

ClearAll[x, a, b, c, g]

y1 = d (a*x^g + b) + c*(2*b + a*x^g)

c = 1; d = 0;
dat1 = {}; curve1 = {}; curve2 = {}; 
    Do[k = y1 /. {x -> p[[i]], a ->  .1, b -> 10. , 
            g -> 2}; AppendTo[dat1, {c, d, p[[i]], k}]; 
      AppendTo[curve1, {p[[i]], k}];, {i, 1, 20}];
c = 0; d = 1;
Do[k = 
        y1 /. {x -> p[[i]], a ->  .1, b -> 10. , g -> 2}; 
      AppendTo[dat1, {c, d, p[[i]], k}]; 
      AppendTo[curve2, {p[[i]], k}];, {i, 1, 20}]

p1 = ListPlot[curve1]

p2 = ListPlot[curve2]

Show[p1, p2]

ClearAll[c, d, x]

indepvars = {c, d, x}

r = 
      NLRegression[dat1, y1, 
          indepvars, {}, {}, {{a, 0, 1}, {b, 0, 20}, {g, 1, 3}},  .0000001] // 
        Timing

rr = r[[2]]

rrr = rr[[3]]

pars = rrr[[2]]

y1

pars2 = Join[pars, {c -> 1, d -> 0}]

p3 = Plot[y1 /. pars2, {x, 1, 20}]

pars3 = Join[pars, {c -> 0, d -> 1}]

p4 = Plot[y1 /. pars3, {x, 1, 20}]

Show[p1, p2, p3, p4]

We see that we have simultaneously fit both curves.  Problems with multiple independent variables can be fit even if they can not be graphed in the manner above.

VI  MAXIMUM LIKELIHOOD ESTIMATION: THE MaxLikelihood FUNCTION

VI.1  Introduction

In this section, the MaxLikelihood function is described.  This function is an estimation procedure which attempts to maximize the log-likelihood of a function with respect to a data set.  It is a useful alternative to nonlinear regression (the NLRegression function).  The statistics used are based on Crooke et al. (1999).  The function uses GlobalSearch as the minimizing function, which is needed because some functions require constraints to be solved properly.  It thus accepts all GlobalSearch options, except CompileOption and SimplifyOption.  This is because Compile does not work well with Log[], and because a sum of Log terms will generally not simplify.  It is assumed that the function to be solved is algebraic, and black box functions are not allowed.  Weighting is not yet implemented.  Special functions are pre-programmed for efficiency and ease of use.  The pre-programmed functions can be ascertained as follows:

?LikelihoodModels

The syntax for the function is as follows:

MaxLikelihood[data,expression,independent_variables,constraints,{{var1name,lowbound,highbound}..},tolerance,options]

VI.2  Examples and Built-In Functions

We first illustrate use of the MaxLikelihood function for determining the mean μ  and standard deviation σ  from a random sample that has been generated by the Normal distribution:

			\(TraditionalForm\`f(x; μ, σ) = 
        e\^\(-\(\((x + μ)\)\^2\/\(2  σ\^2\)\)\)\/\(σ \@\( 2  π\)\), \ x ∈ ℝ\).
			
In the ContinuousDistribution package,  this probability distribution density is defined as NormalDistribution.

Needs["Statistics`ContinuousDistributions`"];

Needs["Statistics`MultinormalDistribution`"];

Needs["Statistics`DiscreteDistributions`"];

For demonstration purposes, we draw a set of pseudo random numbers from the Normal distribution with μ = 10 and σ = 1.

T = 100;
exact = {μ -> 10, σ -> 1};
data = 
        Table[Random[
            NormalDistribution[μ, σ] /. 
              exact], {T}];

Next, we define the probability density that we would like to fit as a pure function.

Clear[f];
\(\(\(f = PDF[NormalDistribution[μ, σ], #] &;\)\(
\)
    \)\)

g = f[x]

mle = 
      MaxLikelihood[data, 
          g, {x}, {}, {{μ, 8, 12}, {σ,  .8, 1.5}},  .000001, 
          CompileOption -> True] // Timing

We see that we have recovered the true parameters pretty well.  A large sample size would be needed for a better estimate.  We can compare the results using the built-in function NormalModel:

mle = 
      MaxLikelihood[data, 
          NormalModel, {x}, {}, {{μ, 8, 12}, {σ,  .8, 
              1.5}},  .000001, CompileOption -> True] // 
        Timing

We see that the exact same result is obtained, but the run is about 4 times faster.  Note that the order of the parameters, but not their names, is assumed to match the {mean,variance} parameters of the normal model.  In this example, since the data is a sample from the distribution, we do not get back exactly the parameters we used to generate the data unless we draw a very large sample data set.

T = 200;
exact = {μ -> 10, σ -> 1};
data = 
        Table[Random[
            NormalDistribution[μ, σ] /. 
              exact], {T}];

mle = 
      MaxLikelihood[data, 
          NormalModel, {x}, {}, {{μ, 8, 12}, {σ,  .8, 
              1.5}},  .000001] // Timing

We see above that a larger sample size gave a much better estimate of the parameters.  Note also that the execution time is nearly unaffected by the size of the data set for the built in model, although time required goes up linearly with data set size for the standard function.

Next we consider the estimation of the mean λ for a sample data drawn from the Poisson distribution:

		\(TraditionalForm\`\(\(f(
          x; λ) = \(\(e\^\(-λ\)\) λ\^x\)\/\(x!\)\)\(,\)\(x = 0\)\(,\)\(1\)\(,\)\(2\)\(,\)\(…\)\(\ \)\)\).
		
We use the same approach to illustrate the method, generating a sample with a known mean and then estimating the mean from the sample.  We draw sample data from a distribution with λ = 10.

Clear[f];
T = 100;
exact = λ -> 10;
data = 
        Table[Random[
            PoissonDistribution[λ] /. exact], {T}];
f = PDF[PoissonDistribution[λ], #] &;

g = f[x]

mle = 
      MaxLikelihood[data, g, {x}, {}, {{λ, 6, 9}},  .000001] // 
        Timing

mle = 
      MaxLikelihood[data, 
          PoissonModel, {x}, {}, {{λ, 6, 9}},  .000001] // 
        Timing

Once again, the preprogrammed model is faster.  Next we evaluate the Beta function.

Our next example is the Beta Distribution

		\(TraditionalForm\`f(
          x; α, β) = \(\(Γ(α + β)\) \(\((1 - x)\)\^\(q - 1\)\) x\^\(p - 1\)\)\/\(\(Γ(α)\) \(Γ(β)\)\), 0 ≤ x ≤ 1\)
		
where α, β > -1 are parameters and Γ denotes the gamma function.   As above, we draw sample data from the specified distribution.

Clear[f];
T = 200;
exact = {α -> 2, β -> 5};
data = 
        Table[Random[
            BetaDistribution[α, β] /. exact], {T}];
f = PDF[BetaDistribution[α, β], #] &;

g = f[x]

In this case, we must add constraints or the function goes to negative infinity during the solution.

mle = 
      MaxLikelihood[data, 
          g, {x}, {-α + 1, -β}, {{α, 1, 
              5}, {β, 3, 8}},  .000001] // Timing

mle = 
      MaxLikelihood[data, 
          BetaModel, {x}, {-α + 1, -β + 1}, {{α, 
              1, 5}, {β, 3, 8}},  .000001] // Timing

Next we illustrate the Gamma distribution.  Note that the order of parameters is assumed to match that in the Gamma built in function.

Clear[f];
T = 100;
exact={α->1,λ->1};
data=Table[Random[GammaDistribution[α,λ]/.exact],{T}];
f=PDF[GammaDistribution[α,λ],#]&;

g = f[x]

Note that we must constrain the parameters to be positive to prevent the function running away to -Infinity.  In addition, lambda must be constrained away from 0 to prevent underflow during computations, even though the final answer is not near zero.

mle=MaxLikelihood[data,g,{x},{-α +  .0001,α - 5,-λ+.0001,λ - 8},{{α,  .2, 5},{λ,3,8}},.000001]//Timing

mle=MaxLikelihood[data,GammaModel,{x},{-α +  .0001,α - 5,-λ+.0001,λ - 8},{{α,  .2, 5},{λ,3,8}},.000001]//Timing

In this case we get a huge speedup (10 times faster).  Next, we test the LogNormal predefined model.

Clear[f];
T = 100;
exact = {μ -> 20, σ -> 5};
data = 
        Table[Random[
            LogNormalDistribution[μ, σ] /. 
              exact], {T}];
f = PDF[LogNormalDistribution[μ, σ], #] &;

g = f[x]

mle = 
      MaxLikelihood[data, 
          g, {x}, {}, {{μ, 8, 12}, {σ,  .8, 10}},  .000001] // 
        Timing

mle = 
      MaxLikelihood[data, 
          LogNormalModel, {x}, {}, {{μ, 8, 12}, {σ,  .8, 
              10}},  .000001] // Timing

VII  DISCRETE VARIABLE PROBLEMS: THE 
 InterchangeMethodMin & TabuSearchMin FUNCTIONS

V.1I  Introduction

InterchangeMethodMin and TabuSearchMin are functions which maximize a linear or nonlinear function of integer 0-1 variables either with or without constraints.  This type of problem comes up often in operations research (e.g., network, scheduling, and allocation applications).  Integer 0-1 problems are not generally solvable with an analytic approach such as Linear Programming.  The method used here is the Interchange method (Densham and Rushton, 1992; Goldberg and Paz, 1991; Lin and Kernighan, 1973; Teitz and Bart, 1968), which is similar to the author's SWAP algorithm (Loehle, 2000).  It has been shown that the Interchange method is guaranteed to get close to the optimum solution, but for complex problems may not find the exact solution.    On the other hand, it is much faster than Simulated Annealing (Murray and Church, 1995).  The Interchange algorithm begins with a feasible start generated randomly.  It then makes changes in the configuration at each iteration until no more progress is possible.  Note that while it is possible to input this type of problem into MultiStartMin because it accepts discrete variables, MultiStartMin assumes that only a few discrete variables are involved, and may generate a large set of possible moves and run slowly for large problems.  TabuSearchMin operates with a basic Interchange framework, but then adds a tabu list feature.  The tabu list is the list of n previous moves that will not be revisited in looking for the next good move.  This feature adds efficiency because bad moves are not checked repeatedly.

The functions are defined by

InterchangeMethodMin[expr,ineqs,{varlist},tolerance,options]

TabuSearchMin[expr,ineqs,{varlist},tolerance,options]

where expr is the function to be minimized, which need not be linear, and ineqs is the inequalities, which are optional and are in standard form. There are three options for the program.  Starts defines the number of random starting points to test (default 5).  Input starting points can be input with the StartsList option.  CompileOption can be set to False if the function is not compilable.  ShowProgress shows intermediate stages of the computation.  Program operation is shown next in the following examples of applications.  Note that the range of applications is not exhausted by these examples.

VII.2  Applications

VII.2.A Capital Allocation

In many investment problems, the goal is to allocate a fixed amount of money across a series of investments.  When the items to be allocated are discrete, and must be allocated to one or another use in their entirety, then the problem must be solved with integer programming.  In the following example, the variables {x1,x2,x3,x4} are 0-1.  There are several constraints.  The objective function is given by

Max (z)= .2  x1 +  .3  x2 +  .5  x3 +  .1  x4

This problem may be solved as follows.  We minimize the negative of z to maximize, and enter all the constraints, which are less than inequalities:

ClearAll[x1, x2, x3, x4]

InterchangeMethodMin[-( .2  x1 +  .3  x2 +  .5  x3 +  .1  
                x4), { .5  x1 + x2 + 1.5  x3 +  .1  x4 - 
            3.1,  .3  x1 +  .8  x2 + 1.5  x3 +  .4  x4 - 
            2.5,  .2  x1 +  .2  x2 +  .3  x3 +  .1  x4 -  .4}, {x1, x2, x3, 
          x4},  .001, Starts -> 4] // Timing

Seven random starts were used, but only 2 unique starts were used due to duplications.  Two starts found the correct answer {0,0,1,1}.  On larger problems, performance might not be quite this good.  TabuSearch also can solve this problem correctly:

TabuSearchMin[-( .2  x1 +  .3  x2 +  .5  x3 +  .1  x4), { .5  
              x1 + x2 + 1.5  x3 +  .1  x4 - 3.1,  .3  x1 +  .8  x2 + 
            1.5  x3 +  .4  x4 - 
            2.5,  .2  x1 +  .2  x2 +  .3  x3 +  .1  x4 -  .4}, {x1, x2, x3, 
          x4},  .001, Starts -> 4, TabuListLength -> 1] // 
      Timing

Note that Tabu Search may perform worse than Interchange on small problems.  Note also that a small tabu list size is needed for this problem.  If we specify a tabu list too long, a warning results:

TabuSearchMin[-( .2  x1 +  .3  x2 +  .5  x3 +  .1  x4), { .5  
              x1 + x2 + 1.5  x3 +  .1  x4 - 3.1,  .3  x1 +  .8  x2 + 
            1.5  x3 +  .4  x4 - 
            2.5,  .2  x1 +  .2  x2 +  .3  x3 +  .1  x4 -  .4}, {x1, x2, x3, 
          x4},  .001, Starts -> 4, TabuListLength -> 100] // 
      Timing

VII.2.B Vehicle Routing/Travelling Salesman

A very common problem in operations research is vehicle routing.  This problem arises in delivery service, airline, and other transportation contexts.  To illustrate, define a set of 8 cities with coordinates:

c={{1,1},{2,0.1},{3,0.1},{4,1},{1,2},{2,2},{3,2},{4,2}}

Where the objective is to minimize the sum of the links between cities.  We define the objective function as the sum of the distance*existence of the link over all possible links between cities.  Self-loops (e.g., x11) are removed from the lists.  Constraints are added that force every city to have at least 2 links (for entry and exit).  Finally, the sum of all edges must be at least 8.  We really want exactly 8 links, but the algorithm does not accept equality constraints.  The minimization of the sum of link distances will force the solution to exactly 8. Links are defined between cities such that there are no repeated links (e.g., x13 but no x31) to make it faster.

lis = {x12, x13, x23, x14, x24, x34, x15, x25, x35, x45, x16, x26, x36, 
        x46, x56, x17, x27, x37, x47, x57, x67, x18, x28, x38, x48, x58, x68, 
        x78}

InterchangeMethodMin[
        2 (1.3453624047073711` x12 + 2.193171219946131` x13 + 3 x14 + 
                x15 + (Sqrt[2]) x16 + (Sqrt[5]) x17 + (Sqrt[10]) x18 + 1.` x23 + 
                2.193171219946131` x24 + 2.1470910553583886` x25 + 
                1.9` x26 + 2.1470910553583886` x27 + 
                2.758622844826744` x28 + 1.3453624047073711` x34 + 
                2.758622844826744` x35 + 2.1470910553583886` x36 + 
                1.9` x37 + 
                2.1470910553583886` x38 + (Sqrt[10]) x45 + (Sqrt[5]) x46 + (Sqrt[2]) x47 + 
                x48 + x56 + 2 x57 + 3 x58 + x67 + 2 x68 + x78) + 
          Abs[x12 + x13 + x14 + x15 + x16 + x17 + x18 - 2] + 
          Abs[x12 + x23 + x24 + x25 + x26 + x27 + x28 - 2] + 
          Abs[x13 + x23 + x34 + x35 + x36 + x37 + x38 - 2] + 
          Abs[x14 + x24 + x34 + x45 + x46 + x47 + x48 - 2] + 
          Abs[x15 + x25 + x35 + x45 + x56 + x57 + x58 - 2] + 
          Abs[x16 + x26 + x36 + x46 + x56 + x67 + x68 - 2] + 
          Abs[x17 + x27 + x37 + x47 + x57 + x67 + x78 - 2] + 
          Abs[x18 + x28 + x38 + x48 + x58 + x68 + x78 - 
              2], {-(x12 + x13 + x14 + x15 + x16 + x17 + x18 + x23 + 
                  x24 + x25 + x26 + x27 + x28 + x34 + x35 + x36 + x37 + x38 + 
                  x45 + x46 + x47 + x48 + x56 + x57 + x58 + x67 + x68 + 
                  x78) + 8}, lis,  .01, Starts -> 5, 
        CompileOption -> True, ShowProgress -> False] // 
      Timing

ll = %;

ll3 = ll[[2, 1]]

We see that the function found the solution in 4 out of 5 starts, which is a loop around the points in a rough circle.  Even the starts that failed got quite close to the best solution and had the proper number of links defined.  The timing is quite good.  The problem of subgraphs that are not linked to form a complete vehicle route can occur in this problem but is unlikely unless the subgraphs produce exactly the same total travel distance as the linked route.

TabuSearchMin[
        2 (1.3453624047073711` x12 + 2.193171219946131` x13 + 3 x14 + 
                x15 + (Sqrt[2]) x16 + (Sqrt[5]) x17 + (Sqrt[10]) x18 + 1.` x23 + 
                2.193171219946131` x24 + 2.1470910553583886` x25 + 
                1.9` x26 + 2.1470910553583886` x27 + 
                2.758622844826744` x28 + 1.3453624047073711` x34 + 
                2.758622844826744` x35 + 2.1470910553583886` x36 + 
                1.9` x37 + 
                2.1470910553583886` x38 + (Sqrt[10]) x45 + (Sqrt[5]) x46 + (Sqrt[2]) x47 + 
                x48 + x56 + 2 x57 + 3 x58 + x67 + 2 x68 + x78) + 
          Abs[x12 + x13 + x14 + x15 + x16 + x17 + x18 - 2] + 
          Abs[x12 + x23 + x24 + x25 + x26 + x27 + x28 - 2] + 
          Abs[x13 + x23 + x34 + x35 + x36 + x37 + x38 - 2] + 
          Abs[x14 + x24 + x34 + x45 + x46 + x47 + x48 - 2] + 
          Abs[x15 + x25 + x35 + x45 + x56 + x57 + x58 - 2] + 
          Abs[x16 + x26 + x36 + x46 + x56 + x67 + x68 - 2] + 
          Abs[x17 + x27 + x37 + x47 + x57 + x67 + x78 - 2] + 
          Abs[x18 + x28 + x38 + x48 + x58 + x68 + x78 - 
              2], {-(x12 + x13 + x14 + x15 + x16 + x17 + x18 + x23 + 
                  x24 + x25 + x26 + x27 + x28 + x34 + x35 + x36 + x37 + x38 + 
                  x45 + x46 + x47 + x48 + x56 + x57 + x58 + x67 + x68 + 
                  x78) + 8}, lis,  .01, Starts -> 5, 
        TabuListLength -> 60] // Timing

We see that the TabuSearchMin function found the solution in 1 out of 5 starts, which is a loop around the points in a rough circle, but in 1/3 the time of the InterchangeMethod.  Even the starts that failed got quite close to the best solution and had the proper number of links defined.  The timing is quite good.  The problem of subgraphs that are not linked to form a complete vehicle route can occur in this problem but is unlikely unless the subgraphs produce exactly the same total travel distance as the linked route.

VII.2.C Minimum Spanning Tree

The construction of a minimum spanning tree is a problem that comes up in problems of delivery from fixed warehouses, utility network design, and other applications.  We begin with the weighted graph, with nodes (a,b,c,d,e,f), and with edges and weights defined by:

{xab,xac,xae,xbc,xbd,xcd,xce,xcf,xdf,xef}

{7, 2, 2, 1, 3, 3, 4, 4, 5, 6}

This can be solved with the interchange method.  The objective function is to minimize the sum of the weighted edges.  This produces the minimum tree we desire.  We know that we want 5 edges as the final solution (#nodes -1), so we add an inequality that the sum of the edges must be greater than 5.  This allows more than 5 initially, but then edges are trimmed as the algorithm proceeds.  If an equality constraint is put in, the program will fail.  Additional constraints force every node to have at least one edge touching it.

InterchangeMethodMin[
        7  xab + 2  xac + 2  xae + 1  xbc + 3  xbd + 3  xcd + 4  xce + 
          4  xcf + 5  xdf + 
          6  xef, {-(xab + xac + xae + xbc + xbd + xcd + xce + xcf + 
                  xdf + xef) + 5, -(xab + xac + xae) + 
            1, -(xab + xbc + xbd) + 
            1, -(xac + xbc + xcd + xce + xcf) + 
            1, -(xbd + xcd + xdf) + 1, -(xae + xce + xef) + 
            1, -(xcf + xdf + xef) + 1}, {xab, xac, xae, xbc, xbd, 
          xcd, xce, xcf, xdf, xef},  .1, Starts -> 1, 
        CompileOption -> True] // Timing

We see that the function returns the correct solution {ac,ae,bc,bd,cf} with a very fast execution time.  This reflects the fast computational time for binary variables.  Large problems can thus be solved.

TabuSearchMin[
        7  xab + 2  xac + 2  xae + 1  xbc + 3  xbd + 3  xcd + 4  xce + 
          4  xcf + 5  xdf + 
          6  xef, {-(xab + xac + xae + xbc + xbd + xcd + xce + xcf + 
                  xdf + xef) + 5, -(xab + xac + xae) + 
            1, -(xab + xbc + xbd) + 
            1, -(xac + xbc + xcd + xce + xcf) + 
            1, -(xbd + xcd + xdf) + 1, -(xae + xce + xef) + 
            1, -(xcf + xdf + xef) + 1}, {xab, xac, xae, xbc, xbd, 
          xcd, xce, xcf, xdf, xef},  .1, Starts -> 1, 
        TabuListLength -> 10, CompileOption -> True] // 
      Timing

We see that the function returns the correct solution {ac,ae,bc,bd,cf} with a very fast execution time. For such smaller problems, the Tabu Search approach does not provide superior speed.

VIII  THE MaxAllocation FUNCTION

VIII.1  Introduction

MaxAllocation is a function to maximize a nonlinear function subject to a linear constraint.

\(Max \((f \((x\_i)\))\)\ s . t  .0 < x\_i < 
        rhs, \ \(\(∑x\_i\)\+\(i = 1\)\%n\) \(\(=\)\(rhs\)\)\)

This type of problem comes up often in finance, economics, and investment.  Traditional algorithms have a great deal of trouble with such a problem, though apparently simple.  A new algorithm enables a very efficient solution of this type of problem.  The algorithm is based on the idea of path following.  If one can stay on the constraint line, then the solution stays feasible at all times.  An approach to path following is based on the SWAP algorithm (Loehle, 2001). We first discretize the interval of the rhs.  This discretization is made proportional to the number of variables, N, for convenience.  The initial 2N pieces are placed onto all the variables equally (each variable is given a value 2*rhs/N).  At iteration 1, each variable is tested to see if swapping a single piece from that variable would improve the solution.  This performs a rectangular move in N space.  After no more pieces can be moved from variable 1 (either all have been moved or no improvement is obtained), testing proceeds to variable 2 to see if any pieces can be moved, and so on.  Note that some of the pieces moved from 1 to 2 might be redistributed at this point.  After the algorithm has tested all the variables in order for the possibility of making a swap, the number of pieces on each variable is doubled (with half the length) for the next iteration.  This enables a finer resolution of the function.

The function is defined by

MaxAllocation[function, varlist, rhs, tolerance, options]

where rhs is the right hand side value (from eqn (1)) and tolerance defines the stopping criterion.  There are three options for the program.  MaxIterations will stop the program at the specified number of iterations if desired.  CompileOption is useful if the function is not compilable.  ShowProgress shows intermediate stages of the computation.

Operation of the function is shown in the following example.  This example is a sum of diminishing return terms.  This is a function with no nonlinear interaction terms.  It can be shown that the algorithm converges to the analytic optimum in this case.

ClearAll[x2]

varlist = {x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, 
          x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x26, 
          x27, x28, x29, x30, x31, x32, x33, x34, x35, x36, x37, x38, x39, 
          x40, x41, x42, x43, x44, x45, x46, x47, x48, x49, x50, x51, x52, 
          x53, x54, x55, x56, x57, x58, x59, x60};

v[x_] := 
      Sum[i*(1.0 - E^(-3. *x[[i]])), {i, 1, 60}]

Note the technique used here to take advantage of indexed variables.  The function v is defined with a tensor as an argument.  Then function vvv is defined to accept a list and turn it into a tensor by passing it down to v.  This function can be compiled, and is in the following example.

vvv = 
        Function[{x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, 
            x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x26, 
            x27, x28, x29, x30, x31, x32, x33, x34, x35, x36, x37, x38, x39, 
            x40, x41, x42, x43, x44, x45, x46, x47, x48, x49, x50, x51, x52, 
            x53, x54, x55, x56, x57, x58, x59, x60}, 
          v[{x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, 
              x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x26, 
              x27, x28, x29, x30, x31, x32, x33, x34, x35, x36, x37, x38, 
              x39, x40, x41, x42, x43, x44, x45, x46, x47, x48, x49, x50, 
              x51, x52, x53, x54, x55, x56, x57, x58, x59, 
              x60}]];

t = Table[1, {i, 1, 60}]

v[t]

Apply[vvv, t]

MaxAllocation[
        vvv[x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15, 
          x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x26, x27, x28, 
          x29, x30, x31, x32, x33, x34, x35, x36, x37, x38, x39, x40, x41, 
          x42, x43, x44, x45, x46, x47, x48, x49, x50, x51, x52, x53, x54, 
          x55, x56, x57, x58, x59, x60], varlist, 1.0,  .1, 
        CompileOption -> True, ShowProgress -> False] // 
      Timing

Next, we can compare the value of compiling the function in terms of timing.

MaxAllocation[
        vvv[x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15, 
          x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x26, x27, x28, 
          x29, x30, x31, x32, x33, x34, x35, x36, x37, x38, x39, x40, x41, 
          x42, x43, x44, x45, x46, x47, x48, x49, x50, x51, x52, x53, x54, 
          x55, x56, x57, x58, x59, x60], varlist, 1.0,  .1, 
        CompileOption -> False] // Timing

We can see that the compiled version runs 5 times faster.

It would be useful to have a proof that this algorithm works. At this time, no proof is available.  The algorithm has been tested on a suite of problems with known solutions, and gives good results.   The algorithm has been designed to guard against infinite loops.   It is possible that a problem could be devised with local minima that would prevent a global optimum from being reached, but no such problem has yet been encountered. The initial equitable distribution of pieces across the variables guards against interaction (e.g., x1*x2) type terms which can not be improved by rectangular moves if they are initialized to zero.

The above problem with 250 variables took 10 minutes to run on a Pentium III 600 MHz machine.  Because the execution time goes up nonlinearly with the number of variables, it is likely that the upper limit for an overnight run is 350-1000 variables, depending on the machine used.

VIII.2  Applications

VIII.2.A Investment Allocation Problems

In many investment problems, the goal is to allocate a fixed amount of money across a series of investments.  We need not invest in all the options available.  A typical problem to solve is the following, (from A.K. Dixit, Optimization in Economic Theory):  A capital sum C is available for allocation among n investment projects.  The expected return from a portfolio of (Subscript[x, j]) projects is

(∑_(j = 1)^n)[(Subscript[α, j]) (Subscript[x, j]) - 
                0.5  (Subscript[β, j]) (Subscript[x, j]^2)]

subject to:

(∑_(j = 1)^((Subscript[nx, j])))=C

This problem may be solved as follows:

a = {1000, 2000, 3000};

rhs = 10.0;

b = {100, 500, 1000. };

h[x_] := 
        Sum[a[[i]]*x[[i]] -  .5*b[[i]]*x[[i]]^2, {i, 1, 
            3}];

fff = Function[{a1, a2, a3}, h[{a1, a2, a3}]];

\(\(\(MaxAllocation[fff[a1, a2, a3], {a1, a2, a3}, rhs,  .1, 
      ShowProgress -> True]\)\(
\)\)\)

Here we find the solution which is an interesting mix of investments.  This problem can be extended to 350 to 1000 variable cases for overnight runs, depending on machine speed.

VIII.2.B Derivative Hedging with Quadratic Programming

A common problem in finance is derivative construction to hedge an investment.  One technique for derivative construction is to construct a target payoff value (T) for the instrument, and then to allocate investment to a series of instruments with return R to get as close to the payoff value as possible using a quadratic objective function (for 3 instruments here):

(((∑_(i = 1)^3)[(Subscript[T, i]) - (∑_(j = 1)^3) 
                            (Subscript[x, j]) (Subscript[R, i, j])])^2)

(∑_(j = 1)^3) (Subscript[x, j]) (Subscript[p, j]) ≤ C

This can be solved with MaxAllocation.  First, for a problem with 3 possible investments, we make a dummy variable (index 4 here) to take up the slack in the investment (to represent the part of C NOT invested).  Also, the equation is rearranged so that x represents the dollar amount invested in each instrument, so x is divided by the price p in eqn. (3).  Then C represents our rhs and the problem is given by:

p = {1, 1, 1};

targetpayoff = {1000, 2000, 3000};

payoff = {{1000, 0, 0}, {0, 2000, 0}, {0, 0, 3000}};

\(\(fy[
          x_] := \(-Sum[\((targetpayoff[\([i]\)] - 
                  Sum[\((x[\([j]\)]/p[\([j]\)])\)*
                      payoff[\([i, j]\)], 
						{j, 1, 3}])\)^2, {i, 
              1, 3}]\);\)\)

fff = Function[{a1, a2, a3, a4}, fy[{a1, a2, a3}]];

\(\(\(MaxAllocation[fff[a1, a2, a3, a4], {a1, a2, a3, a4}, 
      6. ,  .1]\)\(
\)\)\)

In this example, we know the exact answer, and the algorithm converged to the exact answer quickly.  Such derivative construction problems can be solved with hundreds of variables.

IX  APPLICATIONS

IX.1  Zeros of a Function/Roots of a Polynomial

Although the Mathematica function RootFind can find the zeros of a function, it can occur that a function is not analytic, in which case RootFind will not work.  For example, a function might need to be calculated recursively, algorithmically, or by simulation.  FindRoot also will only find a single root at a time from a given starting point.  In addition, we may invoke the need for identifying the indifference zone: we may want all the solutions that are within some indifference region of the zero or root.  It is a simple matter to use either GlobalMinima or MultiStartMin to find roots of a function.  Consider the function

(x^2) - 25

where the roots or zeros are obviously +5  and -5.  As formulated, the function looks like

Plot[x^2 - 25,  {x,  -10, 10}];

and the minimum is at -25.  If, however, we evaluate the absolute value of the function

Plot[Abs[x^2 - 25],  {x,  -10, 10}];

then we can simply minimize the function as in all the above examples and we will find the roots or zeros of the function, which are now the minima.  In cases where there is a region of solutions with values near zero, this procedure will likewise allow us to define these regions.

GlobalMinima[Abs[x^2 - 25],  , {{x, -10, 10}},  10, 
          0.00001, 0.1, 0.0000001]//Timing

MultiStartMin[Abs[x^2 - 25], {}, {}, { {x,  -10, 10}},  .00001, 
        Starts -> 4] // Timing

GlobalSearch[Abs[x^2 - 25], {}, {}, { {x,  -10, 10}},  .00001, 
        Starts -> 4] // Timing

We can see that all three functions find the two solutions but on this small problem GlobalMinima is faster.  The other 2 functions also need multiple starts to find both solutions.

IX.2  Integer Programming: the Knapsack Problem

As discussed above, MultiStartMin allows integer variables.  This implies that it can solve integer programming problems.  Optimization problems with integer variables are difficult because they violate the assumptions of continuous methods such as LP or gradient descent.  Various methods have been applied to integer programming problems, including heuristic search.  The method used here for integer variables is a generalized descent with discrete step sizes, combined with limited interchange.  This method is illustrated for the knapsack problem.  In the knapsack problem, the optimum return on packing of discrete variables is desired.  This could occur when packing a truck, where different packages have different shipping values.  On the other hand, the size of each package may differ, and the total space available is constrained.  We next test such a problem.  Because we wish to maximize, we put a negative sign in front of the objective function.

MultiStartMin[-(x1 + 2  x2 + 4  x3 + 15  x4 + 80  x5 + 
              100  x6), {-x1, -x2, -x3, -x4, -x5, -x6, 1  x1 + 3  x2 + 15  x3 + 29  x4 + 70  x5 + 59  x6 - 300}, , {{x1, 0, 
            2, Integer}, {x2, 0, 2, Integer}, {x3, 0, 2, Integer}, {x4, 0, 2, 
            Integer}, {x5, 0, 2, Integer}, {x6, 0, 4, Integer}},  .00001, 
        ShowProgress -> False, Starts -> 10] // Timing

We see that the correct solution {{5,0,0,0,0,5},-505.} was not found  out of 10 starts.  A better approach is to first solve the continuous problem, and then pass the result back to the discrete problem:

GlobalSearch[-(x1 + 2  x2 + 4  x3 + 15  x4 + 80  x5 + 
              100  x6), {-x1, -x2, -x3, -x4, -x5, -x6, 1  x1 + 3  x2 + 15  x3 + 29  x4 + 70  x5 + 59  x6 - 300}, , {{x1, 0, 
            2}, {x2, 0, 2}, {x3, 0, 2}, {x4, 0, 2}, {x5, 0, 2}, {x6, 0, 
            4}},  .00001, ShowProgress -> False, Starts -> 1] // 
      Timing

We next pass this result, after we make it Integer, back to MultiStartMin, obtaining the correct solution:

MultiStartMin[-(x1 + 2  x2 + 4  x3 + 15  x4 + 80  x5 + 
              100  x6), {-x1, -x2, -x3, -x4, -x5, -x6, 1  x1 + 3  x2 + 15  x3 + 29  x4 + 70  x5 + 59  x6 - 300}, , {{x1, 0, 
            2, Integer}, {x2, 0, 2, Integer}, {x3, 0, 2, Integer}, {x4, 0, 2, 
            Integer}, {x5, 0, 2, Integer}, {x6, 0, 4, Integer}},  .00001, 
        ShowProgress -> False, StartsList -> {{0, 0, 0, 0, 0, 5}}] // 
      Timing

IX.3  Differential Equation Models

The following example shows how to fit parameters in a differential equation model.

sol = 
      NDSolve[{y'[x] == y[x], y[0] == 1}, 
        y, {x, 0, 2}]

y[1.5] /. sol

Plot[y[t] /. sol, {t, 0, 2}]

fn = 
      Function[{a}, 
        so = NDSolve[{y'[x] == y[x], y[0] == a}, 
            y, {x, 0, 2}]; j = y[1.5] /. so; j[[1]]]

fn[1]

gg[x_] := Abs[fn[x] - fn[1]]

GlobalSearch[gg, {}, , {{x, 0, 2}},  .00001, 
        CompileOption -> False, SimplifyOption -> False, 
        Starts -> 1] // Timing

In this example we see that we recovered the original parameter a=1.  Note that CompileOption->False was necessary because a Mathematica function was passed in, which prohibits compilation.  SimplifyOption->False was also needed to prevent error messages (although the correct answer is still obtained without it).  It is also possible to fit models in cases where we want to fit the entire time trajectory to the data, in which case the function is evaluated at many points and some measure of fit computed.  FindMinimum is unable to solve this problem by any solution method:

FindMinimum[gg[x], {x, 0, 2}]

In this example we see that we recovered the original parameter a=1.  Note that CompileOption->False was necessary because a Mathematica function was passed in, which prohibits compilation.  SimplifyOption->False was also needed to prevent error messages (although the correct answer is still obtained without it).  It is also possible to fit models in cases where we want to fit the entire time trajectory to the data, in which case the function is evaluated at many points and some measure of fit computed.  FindMinimum is unable to solve this problem by any solution method.

IX.4  Constraint Equations

One of the hurdles to solving nonlinear optimization problems is finding a feasible starting point.  This is particularly true when there are a large number of constraints.  In the following example, it is clear that the feasible region is the positive quarter of a sphere.

ClearAll[z]

constraints = {-x, -y, -z, x^2 + y^2 + z^2 - 2}

It is possible to use optimization to find the feasible region.  We set up the problem to solve for a constant objective function, so that the program stops when it finds feasible initial points.

MultiStartMin[2, 
        constraints, , {{x, -5, 5}, {y, -5, 5}, {z, -5, 
            5}},  .001, Starts -> 10, ShowProgress -> False] // 
      Timing

GlobalSearch[2, 
        constraints, , {{x, -5, 5}, {y, -5, 5}, {z, -5, 
            5}},  .001, Starts -> 10, ShowProgress -> False] // 
      Timing

The above result gives a set of feasible points, which gives a basis for defining the feasible region.

X.  Literature Cited

Cadzow, J. A. and Martens, H.R., 1970. Discrete-Time and Computer Control Systems, Prentice-Hall, Inc.

Courrieu, P.  1997.  The hyperbell algorithm for global optimization: A random walk using Cauchy densities.  J. Global Optimization 10:37-55.

Crooke, P., L. Froeb, and S. Tschantz.  1999.  Maximum likelihood estimation.  Mathematica in Education and Research 8:17-23.

Csendes, T.  1985.  A simple but hard-to-solve global optimization test problem.  IIASA Workshop on Global Optimization, Sopron, Hungary.

Cvijovic, D. and J. Klinowski.  1995.  Taboo search: an approach to the multiple minima problem.  Science 267:664-666.

Densham, P., and G. Rushton.  1992.  A more efficient heuristic for solving large p-nedian problems.  Pap. Reg. Sci. 71:307-329.

Dixon, L.C.W., and G.P. Szego.  1978.  The global optimization problem: An introduction.  Towards Global Optimization 2, North-Holland, Amsterdam.

Goldberg, J. and L. Paz.  1991.  Locating emergency vehicle bases when service time depends on call location.  Transportation Science  25:264-280.

Griewank, A.O.,  19 ___.  Generalized descent for global optimization.

Ingber, L., and B. Rosen.  1992.  Genetic algorithms and very fast simulated annealing: A comparison.  Mathematical and Computer Modelling 16:87-100.

Jones, D.R., C.D. Perttunen, and B.E. Stuckman.  1993.  Lipschitzian optimization without the Lipschitz constant.  J. Optimization Theory and Applications 79:157-181.

Lin, S. and B. Kernighan.  1973.  An effective heuristic algorithm for the traveling salesman problem.  Operations Research 21:498-516.

Loehle, C. 2000.  Optimal control of spatially distributed process models.  Ecological Modelling 131:79-95. 

Murray, A.T. and R.L. Church.  1995.  Heuristic solution approaches to operational forest planning problems.  OR Spektrum 17:193-203.

Nelson, J. and J.D. Brodie.  1990.  Comparison of random search algorithm and mixed integer programming for solving area-based forest plans.  Canadian J. Forest Research 20:934-942.

Pinter, J.D.  1996.  Global Optimization in Action.  Kluwer Academic Publishers, Boston.

Phillips, C. L. and Nagle, H. T., 1995. Digital Control System Analysis and Design, 3rd edition, Prentice-Hall, Inc.

Shahian, B. and Hassul, M., 1993. Control System Design Using MATLAB, Prentice-Hall, Inc.

Teitz, M. and P. Bart.  1968.  Heuristic methods for estimating the generalized vertex median of a weighted graph.  Operations Research 16:955-961.