NumericalMath`NMinimize

A New Standard Package for Global Optimization

Brett Champion
Wolfram Research
brettc@wolfram.com

<< Projects`NumericalMath`NMinimize`;
SeedRandom[1001];
Off[NMinimize::switch];
Off[General::spell];
Off[General::spell1];

Abstract

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 robust enough to allow us to find global minima, and are not easily trapped by local minima.

We will discuss the user interface to NMinimize, as well as detailing the methods and techologies used by NMinimize. We will also demonstrate how NMinimize performs on several examples of varying complexity.

For an example, let us consider Rosenbrock's function in two dimensions, which has a minimum at x=y=1:

NMinimize[{(((y - (x^2)))^2) + (((1 - x))^2),  -2.084 ≤ x ≤ 2.084 && -2.084 ≤ y ≤ 
            2.084},  {x, y}]

Introduction

Mathematica has many functions for finding minima of a function, but none that are suitable for finding a constrained global minimum of a nonlinear objective function. NMinimize aims to fill that need by implementing several algorithms for finding global constrained optima. The available methods are flexible enough to cope with functions that are not differentiable or continuous, and are not easily trapped by local minima. 

NMinimize is intended to fill a gap in the group of optimization functions currently available in Mathematica.

({{ , Objective

            Function, Constraints, Global

            Minimization, Other}, {FindMinimum, Any, None, No, Needs 

              initial points}, {ConstrainedMin, Linear, Linear, Yes, Needs rational

              exact coefficients}, {Minimize, Polynomial, Polynomial, Yes, Needs rational

              coefficients}, {NMinimize, Any, Any, Yes, Needs initial

              region}})

It should, however, be kept in mind that finding a global minimum, even without constraints, can be arbitrarily difficult, and so the methods used may fail. Frequently it may be useful to minimize the function several times, and take the best of the results.

Interface

The general syntax for NMinimize, with objective function obj, constraints cons, and variables vars, is:

NMinimize[obj,  vars]

for an unconstrained problem, and

NMinimize[{obj,  cons},  vars]

for a constrained problem.

The constraints can be either a list or a logical combination of equalities, inequalities, and domain specifications. NMinimize only handles weak inequalities, and automatically converts strong inequalities into weak inequalities. We can also specify in the constraints that certain variables should be integers or reals by using Element.

NMinimize also requires an initial rectangular region in which to create its starting points.  This can be done be giving variables as a list containing the variable followed by the bounds for the initial region, such as {x, a, b}, or by putting inequalities in the constraint, such as a ≤ x ≤ b.

For example if we want to minimize (x^2) + (y^2)subject to the constraint 1 ≤ x ≤ 5 with the initial region where x and y are both between 3 and 5, we would use the following.

NMinimize[{(x^2) + (y^2),  
        1 ≤ x ≤ 5}, {{x,  3,  5},  {y,  3,  
          5}}]

Methods

NMinimize provides five methods for finding global optima: DifferentialEvolution, NelderMead, RandomSearch, SimulatedAnnealing, and Automatic (default.)

DifferentialEvolution

DifferentialEvolution is a genetic method developed by K. Price and R. Storn, and is perhaps the most robust method. This comes at a large cost, since it maintains a large set of points for a gene pool, and is frequently slower than the other methods, since it needs more function evaluations (one per gene in the gene pool per iteration.) This method is the best at discrete optimization problems, although some parameters may need to be tuned for such problems.

NelderMead

NelderMead is the simplex method developed by J. A. Nelder and R. Mead, and is generally the fastest method, requiring only a few function evaluations per iteration. Since NelderMead is deterministic once the initial simplex has been chosen, this method is the most likely to get trapped by a local minimum. This method is good for problems that have relatively large continuous components, with relatively "nice" topologies.

RandomSearch

RandomSearch is a simple method, where we generate a large quantity of points in the initial region, and use each point as the starting point for FindMinimum. This is generally fast, but doesn't scale well, and requires that the objective function be locally continuous.

SimulatedAnnealing

SimulatedAnnealing is an implementation of the simulated annealing algorithm. To improve our minimum we use several starting points for the algorithm, instead of just one. SimulatedAnnealing is generally faster than DifferentialEvolution, but doesn't always perform as well. The method can be used for discrete problems, although like DifferentialEvolution, it may benefit from tuned parameters.

Automatic

Automatic is the default method for NMinimize, and we try to pick a method that is good for the problem we're given. If we're given a linear programming problem, we make a call to ConstrainedMin. If any of the variables are integers, we use DifferentialEvolution. Otherwise we try NelderMead, and if that does poorly, we switch to DifferentialEvolution. This gives a generally fast method, but also a way to drop back to a slower, but more robust method if needed.

Examples

Basic Example

Suppose we want to find the minimum of (x^2) + (y^2) + 10 cos(x y). As we can see below, this function has many local minima, which can confuse a gradient based method such as FindMinimum. We've marked the global minimum below in black.

p3d = 
        Plot3D[(x^2) + (y^2) + 10 Cos[x y],  {x,  -5,  
            5},  {y,  -5,  5},  PlotPoints -> 100,  
          Mesh -> False,  
          DisplayFunction -> Identity];
pt = 
        Graphics3D[{PointSize[0.05], 
            Point[{-1.7147112177123554`, -1.7147113257417117`, -3.9174895055337746`  +  .1}]}];
Show[{p3d,  pt},  
      DisplayFunction -> $DisplayFunction]

NMinimize handles something like this just fine.

NMinimize[
      (x^2) + (y^2) + 
        10 Cos[x y],   {{x,  -5,  5},  {y,  -5,  
          5}}]

Example with Constraints

In this example, let's do something with constraints. Note that we don't need to restrict ourselves to linear constraints.

NMinimize[{1000 - (x^2) - 2  (y^2) - (z^2) - x y - x z,  
        (x^2) + (y^2) + (z^2) == 25  && 8  x + 14  y + 7  z == 56 && 
          x ≥ 0 && y ≥ 0 && 
          z ≥ 0},  {{x,  0,  10},  {y,  0,  10},  {z,  
          0,  10}}]

{(x^2) + (y^2) + (z^2), 8  x + 14  y + 7  z} /. %[[2]] // 
      InputForm

Example with an Integer Variable

Here is an example where we specify domain information. Note that in this case, the nearest integer is 0, but that violates one of the other constraints, so the solution is at x->1

NMinimize[{(((x - (1/3)))^2),  
        x ∈ Integers && x ≥ (1/4)},  {{x,  0,  
          10}}]

A Discrete Optimization Example

A common mathematical puzzle is to arrange 8 queens on a chessboard such that no queen is attacking any other queen.  We'll set this up as an optimization problem, and let NMinimize find a solution for us. To make it a bit easier, we'll set up the problem as a permutation problem, which will mean we don't have to worry about two queens being in the same row or column.

Initialization

(*attackQ[pos1,  
          pos2] is True iff pos1 is attacking pos2*) 
attackQ[{i_Integer, xi_Integer},  {j_Integer, 
          xj_Integer}] := 
      Or[i - j == 0,  xi - xj == 0,  i - j == xi - xj,  
        i - j == xj - xi]

countAttacks[q : {_Real .. }] := 
      Block[{qns = Ordering[q], n = Length[q]}, 
Count[
          Flatten[Table[
              attackQ[{i, qns[[i]]}, {j, 
                  qns[[j]]}],  {i,  
                1,  n - 1}, {j,  i + 1,  n}]],  True]]

grid = 
      Graphics[Join[
          Table[Line[{{1, i}, {9, i}}],  {i,  9}], 
Table[
            Line[{{i, 1}, {i, 9}}],  {i, 9}]]]

QueenPlot[vec_] := 
      Show[grid, 
        Graphics[ 
          Table[Disk[{i +  .5,  
                vec[[i]] +  .5},  
              0.25],  {i,  Length[vec]}]],  
        AspectRatio -> 1]

Optimize

\(\(\( (*create\ some\ variables*) \)\(
\)\(\(vars\  = \ \(Q\_# &\) /@ \ 
          Range[8];\)

 (*create\ some\ vars\ with\ bounds*) 

    \(bvars = \({#, 0, 1} &\) /@ vars;\)\)\)\)

{f, sol} = 
        NMinimize[countAttacks[vars], bvars, 
          Method -> DifferentialEvolution, 
          PostProcess -> False, 
Compiled -> 
            False];
 (*clean up the solution a bit*) 
{f, 
      Thread[vars -> (sol = 
              Ordering[
                sol[[All, 
                  2]]])]}

QueenPlot[sol]