Numerical Optimization in Mathematica: an insider's view of NMinimize

Brett Champion

NumericalMath`NMinimize is a new standard add-on package for finding solutions to constrained and unconstrained global optimization problems. The package provides several derivative-free methods for finding global minima. This allows us to solve problems where the objective function is not differentiable (or even continuous). The methods are sufficiently robust to find global minima insofar as they are not readily trapped by local minima. We will discuss the user interface to NMinimize and also provide details of  the methods and technologies utilized. Technical issues such as how initial regions are chosen, integer variables handled, and the like will also be covered.

NMinimize is a new function in Mathematica 4.2 for solving a wide range of global optimization problems. In order to use this function, we need first to load one of the Standard AddOn Packages:

<< NumericalMath`NMinimize`

The basic syntax for unconstrained problems is

NMinimize[f,  vars]

where f is the function to minimized, and vars is a list of variables, possibly specifying initial search regions. For constrained problems, we have the extended syntax

NMinimize[{f,  cons},  vars]

where cons is the set of constraints represented as either a list or a logical combination of inequalities, equalities, and variable domains. We can also specify different algorithms for minimization, most of which treat the objective function as a black box. Each of the methods has various suboptions that allow fine tuning of the algorithms. The syntax for specifying methods and method options is, respectively

NMinimize[f,  vars,  Method -> method]

NMinimize[f,  vars,  
      Method -> {method,  (Subscript[opt, 1]),  (Subscript[opt, 2]) ... }]

where method is an acceptable method, and the (Subscript[opt, i])are method specific options (those in Options[method]). Normal NMinimize options (those in Options[NMinimize]) are specified at the same level as the Method option. For more information on NMinimize syntax, and examples using NMinimize, please see the NMinimize entry in the Mathematica Help Browser.

It should be noted that this paper discusses very technical aspects of the implementation of NMinimize in Mathematica 4.2, and some of the implementation details or interfaces may change in future versions of Mathematica.

Methods

DifferentialEvolution

DifferentialEvolution is a genetic algorithm developed by K. Price and R. Storn [1]. This method maintains a population of specimens (Subscript[x, 1]), …, (Subscript[x, n]), represented as vectors of d-dimensional real numbers (“genes”), where d is the number of variables. Every iteration, for each vector (Subscript[x, i]), it chooses random integers 1 ≤ a, b, c ≤ n and constructs the mate (Subscript[y, i]) = 
        (Subscript[x, i]) + γ((Subscript[x, a]) + ((Subscript[x, b]) - (Subscript[x, c]))), where γ is the value of ScalingFactor. Then (Subscript[x, i]) is mated with (Subscript[y, i]) according to the value of CrossProbability, giving us the child (Subscript[z, i]). At this point (Subscript[x, i]) competes against (Subscript[z, i]) for position i in the population. The default value of SearchPoints is Automatic, which is Min[10*d, 50].

NelderMead

NelderMead is an implementation of J. A. Nelder and R. Mead's simplex algorithm [2]. For simplicity, we assume here that minimization is being done. We start with the highest vertex (Subscript[v, h]) of a simplex, and reflect it across the centroid (Subscript[v, c]) of the remaining points to a new point (Subscript[v, a]) such that ‖(Subscript[v, a]) - 
            (Subscript[v, c])‖/‖(Subscript[v, h]) - 
            (Subscript[v, c])‖ equals the ReflectRatio. If (Subscript[v, a]) is lower than all other vertices, we expand reflection by finding the point (Subscript[v, b]) such that ‖(Subscript[v, b]) - 
            (Subscript[v, a])‖/‖(Subscript[v, a]) - 
            (Subscript[v, c])‖ equals the ExpandRatio. Then the lower of (Subscript[v, a]) and (Subscript[v, b]) replaces (Subscript[v, h]) in the simplex and the process starts over. If (Subscript[v, a]) is lower than the highest vertex but not lower than the lowest vertex, (Subscript[v, a]) replaces (Subscript[v, h]) and the process starts over. Finally, if (Subscript[v, a]) is the highest vertex, we construct (Subscript[v, b]) so that ‖(Subscript[v, c]) - 
            (Subscript[v, b])‖/‖(Subscript[v, c]) - 
            (Subscript[v, h])‖ equals the ContractRatio. If (Subscript[v, b]) is lower than (Subscript[v, h]), we replace (Subscript[v, h]) with (Subscript[v, b]) and start over; otherwise, we move each of the simplex vertices (Subscript[v, i]) toward the lowest vertex (Subscript[v, l]), giving new vertices (Subscript[v, i]^′) where \(TraditionalForm\`‖v\_i - 
            v\_i\%′‖/‖v\_i - v\_l‖\) matches the value of ShrinkRatio, and start over.

NelderMead also has a mechanism for restarting if it converges to a point that doesn't satisfy the constraints. In this case, it will attempt to follow constraint gradients to find a point that does satisfy the constraints, generate a new simplex around this point, and start again. Each restart has the full MaxIterations available for use. The maximum number of restarts is three if NelderMead is being used as a result of Method->Automatic, and five if NelderMead was specified by the Method option.

RandomSearch

RandomSearch is a simple method that generates a large set of points and uses each one as the starting point for FindMinimum. Since FindMinimum doesn't allow constraints, a modified objective function is used. The default value of SearchPoints is Automatic, which is Min[100*d, 1024], where d is the number of variables.

SimulatedAnnealing

SimulatedAnnealing is an implementation of a biased random-walk search method [3]. It generates a random set of starting points, and for each starting point picks a random direction in which to move. If the move is a better point then it is accepted. If the move is not a better point, we calculate a probability ρ and compare it against a random value r ∈ (0, 1), and if r < ρ we accept the new point despite it not being an improvement. The probability ρ is given by \(TraditionalForm\`\(\(ρ = 
        e\^\(b \(\(\)\([i \(\(,\)\(\ \)\) Δ  f \(\(,\)\(\ \)\) f\_0]\)\)\)\)\(,\)\)\) where b is the function defined by BoltzmannExponent, i is the current iteration, Δ  f is the change in objective function value, and (Subscript[f, 0]) is the value of the objective function from the previous iteration. This is repeated for each starting point until the maximum number of iterations is reached, the method converges to a point based on the accuracy or precision goals, or the method stays at the same point consecutively for the number of iterations given by LevelIterations. When BoltzmannExponent -> Automatic the function used is Function[{i, df, f0}, -df*Log[i+1]/10]. The default value of SearchPoints is Automatic, which is Min[2*d, 50], where d is the number of variables.

Automatic

With Method -> Automatic, NMinimize picks which method to use based on the type of problem. If the objective function and constraints are linear then a pure linear programming routine is called. If there are integer variables, or if the head of the objective function isn't a numeric function, DifferentialEvolution is used. For everything else, it uses NelderMead, and if that does poorly it tries DifferentialEvolution.

Handling of Constraints

NMinimize handles constraints in a number of ways depending on the their type. The first step is to expand the constraints to the Form Or[And[..], And[..], ..] so that we can solve the problem over each of the Or branches. Inequalities are changed to the form g≤0, equalities to h==0, and domain constraints are reduced to a list of integer variables.

If linear equalities are present, NMinimize solves the equations exactly, and plugs the solution into the other constraints, reducing the dimension of the problem. It will not solve for integer variables, although solutions for non-integer variables may contain the integer variables.

Linear inequalities are currently handled differently depending on whether they give a bound on a variable. If so, then after a new point is generated, the point is checked against the bound constraints, and violations are truncated to satisfy the constraints. Other linear inequalities are treated the same as nonlinear inequality constraints, as next described.

Nonlinear constraints, including certain linear inequalities as just noted, are handled with penalty functions applied outside the feasible region. Penalty functions are functions of two variables, with the first variable being the positive constraint violation, and the second variable being the current iteration number. The default penalty function doubles the penalty term multiplier every ten iterations and uses quadratic penalty terms. This can be changed using the PenaltyFunction option. In the following example we double the multiplier every iteration, with quadratic penalty terms:

NMinimize[{obj,  cons},  vars,  
      PenaltyFunction -> 
        Function[{d, i}, (d^2) (2^i)]]

Domains are given by adding Element[x, Integers] or Element[x, Reals] to the constraints. All variables must be either integers or real numbers at this point.

NMinimize handles integer variables (Subscript[x, i]) by substituting Round[(Subscript[x, i])] for (Subscript[x, i]) in the objective function and all constraints. This is made possible by using methods that don't rely on derivatives. It is also possible for users to handle integer variables on their own, e.g. by adding terms like (((x-Round[x]))^2)to the objective function or by using a mapping from real numbers to integers in the objective function.

Choosing an Initial Region

Since NMinimize methods generally need several initial points, NMinimize takes a starting region as input, analogous to the way FindMinimum and FindRoot take starting points. The initial regions are calculated from the variables and constraints passed into NMinimize. Each variable can be given as either the variable by itself, or in a list with the upper and lower bounds of the initial region for that variable. If no initial bound is given in the variable list, we look to the constraints for information, and if nothing there can be used, the default range for the variable is the interval [-1,1]. Here are several examples:

NMinimize[f[x],  {x}]

Since there are no constraints, and no explicit initial region, the default initial region [-1,1] is used.

NMinimize[f[x],  {{x,  -5,  5}}]

Here there still are no constraints but there is an explicit initial region, [-5,5], that is used.

NMinimize[{f[x],  
        0 ≤ x ≤ 5},  {x}]

This example has upper and lower bounds on the variable, but no explicit initial region, so the region defined by the constraints,  [0,5], is used.

NMinimize[{f[x],  
        0 ≤ x ≤ 5},  {{x,  1,  
          6}}]

This example has an explicit initial region, as well as bounds on the variable. We use the intersection of the explicit initial region and the constraint region, giving [1,5] as the actual initial region.

NMinimize[{f[x],  x ≤ 3},  {x}]

Here there is an upper bound, but no lower bound, and no explicit initial region. We use Min[-1, UpperBound-2] as the lower bound, resulting in [-1,3] as the initial region.

NMinimize[{f[x],  x ≥ 0},  {x}]

This example is the opposite of the previous one, with a lower bound, and no upper bound, but still no explicit initial region.  In this case we use Max[1, LowerBound+2] as the upper bound, which gives us [0,2] as the initial region.

The initial region is just that—a place to start. The minimization process is able to move outside the region, so long as constraints aren't violated. It is also possible to specify an explicit set of starting points via the InitialPoints option. This takes a list of starting points, whose ordering is the same as the ordering of the variables given to NMinimize. If a method needs more starting points than given by InitialPoints, it randomly generates points inside the initial region as specified above. If fewer points are needed, it takes points sequentially from the list until enough points are selected.

NMinimize will generate up to ten times the number of needed points in an attempt to find points that additionally satisfy the non-rectangular inequality constraints. If it cannot find enough initial points that satisfy the constraints it issues a message and uses those generated points that have the smallest constraint violations. Equality constraints are ignored at this point.

When to Say When

NMinimize uses several criteria in determining when to stop. MaxIterations gives the maximum number of iterations to be used. How this is interpreted may vary from method to method. For example, RandomSearch will use the MaxIterations value as the maximum number of iterations for each of the individual FindMinimum calls. AccuracyGoal and PrecisionGoal behave much as they do in other functions—they provide a way to stop when either the absolute or relative change in function evaluations becomes smaller than a given amount. Most NMinimize methods also look at the sum of squares of the constraint violations, and compare this to the value of the Tolerance option.

NMinimize will run until it reaches the maximum number of iterations, unless the sum of squares of the constraint violations is below Tolerance and either the absolute change in function values drops below (10^(-AccuracyGoal))or the relative change drops below (10^(-PrecisionGoal)). For problems with integer variables, NMinimize runs until the maximum number of iterations are used. One can make NMinimize use the full MaxIterations by setting both AccuracyGoal->∞ and PrecisionGoal->∞. For NelderMead and RandomSearch, AccuracyGoal and PrecisionGoal are checked every iteration, and for DifferentialEvolution and SimulatedAnnealing, every ten iterations.

Post-Processing

NMinimize generally does some post-processing of the results using FindMinimum. This frequently adds a good deal of local refinement once we are in the general area of the global minimum. This is configurable in all methods except Automatic by changing the value of the PostProcess method-option. Possible values for PostProcess are True, False, Automatic, "FindMinimum", or a list of the form {"FindMinimum", <options>}. At present the values Automatic and FindMinimum have the same effect as True. If {"FindMinimum", <options>} is used, then the options listed are passed to FindMinimum in the post-processing stage in place of the default values used by NMinimize.

Since FindMinimum does not accept constraints we need to use penalty functions to enforce the constraints in post-processing. For constraints that are not satisfied when post-processing begins, we use penalty terms based on the function given by the PenaltyFunction option to encourage the solution to move inside the feasible region. For the constraints that are already satisfied, we use inverse barrier terms to keep the solution inside the feasible region.

Since FindMinimum uses descent methods which do not work well with integer variables, we do not try to improve integer components of solutions. Instead we fix the integer variables at the rounded value, and proceed to post-process with the remaining variables (if any). This, of course, may give less than optimal results.

And the Winner Is...

Once we reach the end, we typically have several potential solutions, from which we need to choose the "best" solution. The solutions are first ranked according to the sum of squares of the constraint violations, and we discard any that violate the constraints by more than the value of Tolerance plus the smallest constraint violation. The remaining solutions are then sorted by objective function value, and the best is returned. If this solution is not the one with the smallest constraint violation, a message is issued that notes the candidate solution with smallest constraint violation.

Examples

Many elementary examples may be found in the on-line documentation for NMinimize, which is available at

by following the links Standard Packages->NumericalMath->NMinimize. More examples may be found in talks from the 2001 Mathematica Developer Conference at the following URLs:

Permutation Example

Here is a simple example that demonstrates a technique that works well for permutation problems. We use NMinimize to find the permutation that minimizes a dot product. We define a function on a list of real numbers that turns the list numbers into a permutation using Ordering, and takes the dot product of the numbers one through four with the permutation. The answer from NMinimize is returned in terms of real numbers, and so a bit of work is done to get the answer back in terms of integers.

f[x : {_Real .. }] := 
        Range[4] . Ordering[x];
sol = 
        NMinimize[f[{a, b, c, d}],  {a, b, c, d},  
          Compiled -> False];
sol[[2]] = 
        Thread[{a, b, c, d} -> 
            Ordering[{a, b, c, d} /. 
                Last[sol]]];
sol

SIAM $100 Challenge Example

Here is a solution to the fourth problem in the SIAM $100 Challenge [4]. We use DifferentialEvolution since it tends to be a bit more robust, albeit slower, than the default method.

f = (E^Sin)[50 x] + Sin[60 (E^y)] + 
          Sin[70 Sin[x]] + Sin[Sin[80 y]] - Sin[10 (x + y)] + 
          (1/4) ((x^2) + (y^2));
NMinimize[f,  {x, y},  Method -> "DifferentialEvolution",  
      WorkingPrecision -> 32]

Polynomial Example with Constraints

Here is an example with nonlinear constraints. This problem is tough to solve numerically because multiple constraints are active at the solution [5]. The exact solution is -44 when at {0, 1, 2, -1}.

f[w_Real, x_Real, y_Real,  z_Real] := 
        (w^2) + (x^2) + 2 (y^2) + (z^2) - 5 w - 5 x - 21 y + 7 z + 
          Random[NormalDistribution[]];
c1 = 
        2 (w^2) + (x^2) + (y^2) + 2 w - x - z - 5 ≤ 0;
c2 = 
        (w^2) + (x^2) + (y^2) + (z^2) + w - x + y - z - 8 ≤ 
          0;
c3 = 
        (w^2) + 2 (x^2) + (y^2) + 2 (z^2) - w - x - 10 ≤ 
          0;

NMinimize[{f[w, x, y, z], c1 ∧ c2 ∧ c3}, {w, x, y, z}, 
      Method -> DifferentialEvolution,  
      PrecisionGoal -> ∞,  
      AccuracyGoal -> ∞,  
      MaxIterations -> 500]

Stochastic Objective Function Example

Next we will add some Gaussian noise to the objective function from the previous example.

<< "Statistics`NormalDistribution`";
fr[w_Real, x_Real, y_Real,  z_Real] := 
        (w^2) + (x^2) + 2 (y^2) + (z^2) - 5 w - 5 x - 21 y + 7 z + 
          Random[NormalDistribution[]];

We will use multiple objective function evaluations for improved performance, although this will be slower. We still get solutions near the noise-free minimum, but not as accurately. Our results appear to be comparable in quality to those of dedicated stochastic solvers [6].

(Subscript[f, 5])[w_Real, x_Real, y_Real,  z_Real] := 
        Sum[fr[w, x, y, z], {i, 5}]/5;

NMinimize[{(Subscript[f, 5])[w, x, y, z], c1 ∧ c2 ∧ c3}, {w, x, y, z}, 
      Method -> DifferentialEvolution,  MaxIterations -> 500,  
      PrecisionGoal -> ∞,  
      AccuracyGoal -> ∞]

References

[1]  K. Price and R. Storn, “Differential Evolution,” Dr. Dobb's Journal, April 1997, pp. 18–24, 78. 
[2]  R. Mead and J. A. Nelder, “A Simplex Method for Function Minimization,” Computer Journal, 7 (4), 1965 pp. 308–312. 
[3]  L. Ingber, “Simulated annealing: Practice versus theory,” Mathematical Computer Modelling, 18 (11), 1993 pp. 29–57. 
[4]  N. Trefethen, "A Hundred-dollar, Hundred-digit Challenge,"SIAM News, 35 (1), 2002 p. ??. 
[5]  H.-P. Shwefel, Evolution and Optimum Seeking, New York: John Wiley and Sons, 1995. 
[6]  I.-J. Wang and J. C. Spall, “A Constrained Simultaneous Perturbation Stochastic Approximation Algorithm Based on Penalty Functions,” in Proceedings of the American Control Conference, San Diego, CA.1999 pp. 393–399.