neldermead is a R port of a module originally developed for Scilab version 5.2.1 by Michael Baudin (INRIA - DIGITEO). Information about this software can be found at www.scilab.org. The following documentation as well as the content of the functions .Rd files are adaptations of the documentation provided with the original Scilab neldermead module.

neldermead currently does not include any adaptation of the Scilab ‘nmplot’ function series that is available in the original neldermead module.

Overview

Description

The goal of this toolbox is to provide several direct search optimization algorithms based on the simplex method. The optimization problem to solve is the minimization of a cost function, with bounds and nonlinear constraints.

\[ \begin{array}{l l} min f(x)\\ l_i \le{} x_i \le{} h_i, & i = 1,n \\ g_j(x) \ge{} 0, & j = 0,nb_{ineq} \\\\ \end{array} \]

where \(f\) is the cost function, \(x\) is the vector of parameter estimates, \(l\) and \(h\) are vectors of lower and upper bounds for the parameter estimates, \(n\) is the number of parameters and \(nb_{ineq}\) the number of inequality constraints \(g(x)\).

The provided algorithms are direct search algorithms, i.e. algorithms which do not use the derivative of the cost function. They are based on the update of a simplex, which is a set of \(k \ge n+1\) vertices, where each vertex is associated with one point and one function value.

The following algorithms are available:

  • The fixed shape simplex method of Spendley, Hext and Himsworth: this algorithm solves an unconstrained optimization problem with a fixed shape simplex made of \(k = n+1\) vertices.
  • The variable shape simplex method of Nelder and Mead: this algorithm solves an unconstrained optimization problem with a variable shape simplex made of \(k = n+1\) vertices [1].
  • Box’s complex method: this algorithm solves an constrained optimization problem with a variable shape simplex made of an arbitrary k number of vertices (\(k = 2n\) is recommended by Box).

Basic object

The basic object used by the neldermead package to store the configuration settings and the history of an optimization is a ‘neldermead’ object, i.e. a list typically created by neldermead and having a strictly defined structure (see ?neldermead for more details).

The cost function

The function element of the neldermead object allows to configure the cost function. The cost function is used to compute the objective function value \(f\). If the nbineqconst element of the neldermead object is configured to a non-zero value, the cost function must also compute the value of the nonlinear, positive, inequality constraints \(c\). The cost function can also take as input/output an additional argument, if the costfargument element is configured. The function should be defined as described in vignette('manual', package = 'optimbase'):

  costf <- function(x, index, fmsfundata){
    # Define f and c here #
    return(list(f, g = NULL, c, gc = NULL, index = index,
                this = list(costfargument = fmsfundata)))
  }

where:

  • x: is the current point, as a column vector,
  • index: (optional), an integer representing the value to compute,
  • fmsfundata: an user-provided input/output argument,
  • f: the value of the objective function (a scalar),}
  • g: typically the gradient of the objective function in the context of the optimbase functions; must be set to NULL as the Nelder-Mead is not gradient-based,
  • c: the vector of values of non-linear, positive, inequality constraints,
  • gc: typically the gradient of the constraints in the context of the optimbase functions; must be set to NULL as the Nelder-Mead is not gradient-based, and
  • this: must be set to list(costfargument = fmsfundata).

The index input parameter tells the cost function what to return as output arguments (as described in vignette('manual', package = 'optimbase'). It has the following meaning:

  • index = 2: compute f,
  • index = 5: compute c,
  • index = 6: compute f and c

The fmsdata argument is both input and output. This feature may be used in the situation where the cost function has to update its environment from call to call. Its simplest use is to count the number of calls to the cost function, but this feature is already available directly. Consider the more practical situation where the optimization requires the execution of an underlying Newton method (a chemical solver for example). This Newton method requires an initial guess \(x_0\). If the initial guess for this underlying Newton method is kept constant, the Newton method may have problems to converge when the current optimization point get far away from the its initial point. If a costfargument element is defined in the neldermead object, it can be passed to the cost function as the fmsdata argument. In this case, the initial guess for the Newton method can be updated so that it gets the value of the previous call. This way, the Newton method will have less problems to converge and the cost function evaluation may be faster.

We now present how the feature works. Every time the cost function is called back, the costfargument element is passed to the cost function as an input argument. If the cost function modifies its content in the output argument, the content of the costfargument element is updated accordingly. Once the optimization is performed, the user may call the neldermead.get function and get back an updated costfargument content.

The output function

The outputcommand element of the neldermead object allows to configure a command which is called back at the start of the optimization, at each iteration and at the end of the optimization. The output function must be defined as follows:

outputcmd <- function(state, data, myobj)

where:

  • state: is a string representing the current state of the algorithm. Available values are ‘init,’ ‘iter,’ and ‘done,’
  • data: a list containing at least the following entries:
    • x: the current optimum,
    • fval: the current function value,
    • iteration: the current iteration index,
    • funccount: the number of function evaluations,
    • simplex: the current simplex,
    • step: the previous step in the algorithm. The following values are available: ‘init,’ ‘done,’ ‘reflection,’ ‘expansion,’ ‘insidecontraction,’ ‘outsidecontraction,’ ‘reflectionnext,’ and ‘shrink,’
  • myobj: a user-defined parameter. This input parameter is defined with the outputcommandarg element of the neldermead object.

The output function may be used when debugging the specialized optimization algorithm, so that a verbose logging is produced. It may also be used to write one or several report files in a specialized format (ASCII, LaTeX, Excel, etc…). The user-defined parameter may be used in that case to store file names or logging options.

The data list argument may contain more fields than the current presented ones. These additional fields may contain values which are specific to the specialized algorithm, such as the simplex in a Nelder-Mead method, the gradient of the cost function in a BFGS method, etc…

Termination

The current package takes into account several generic termination criteria. The following termination criteria are enabled by default:

  • maxiter,
  • maxfunevals,
  • tolxmethod,
  • tolsimplexizemethod.

The neldermead.termination function uses a set of rules to compute if the termination occurs and sets optimization status to one of the following: ‘continue,’ ‘maxiter,’ ‘maxfunevals,’ ‘tolf,’ ‘tolx,’ ‘tolsize,’ ‘tolsizedeltafv,’ ‘kelleystagnation,’ ‘tolboxf’ or ‘tolvariance.’ The value of the status may also be a user-defined string, in the case where a user-defined termination function has been set.

The following set of rules is examined in this order.

  • By default, the status is ‘continue’ and the terminate flag is FALSE.
  • The number of iterations is examined and compared to the maxiter element of the neldermead object: if iterations \(\ge\) maxiter, then the status is set to ‘maxiter’ and terminate is set to TRUE.
  • The number of function evaluations is examined and compared to the maxfunevals elements: if funevals \(\ge\) maxfunevals, then the status is set to ‘maxfuneval’ and terminate is set to TRUE.
  • The tolerance on function value is examined depending on the value of the tolfunmethod.
    • FALSE: then the criteria is just ignored,
    • TRUE: if \(|currentfopt| < tolfunrelative \cdot |previousfopt| + tolfunabsolute\), then the status is set to ‘tolf’ and terminate is set to TRUE.
    The relative termination criteria on the function value works well if the function value at optimum is near zero. In that case, the function value at initial guess \(fx0\) may be used as \(previousfopt\). This criteria is sensitive to the \(tolfunrelative\) and \(tolfunabsolute\) elements. The absolute termination criteria on the function value works if the user has an accurate idea of the optimum function value.
  • The tolerance on x is examined depending on the value of the tolxmethod element.
    • FALSE: then the criteria is just ignored,
    • TRUE: if \(norm(currentxopt - previousxopt) < tolxrelative \cdot norm(currentxopt) + tolxabsolute\), then the status is set to ‘tolx’ and terminate is set to TRUE.
    This criteria is sensitive to the tolxrelative and tolxabsolute elements. The relative termination criteria on x works well if x at optimum is different from zero. In that case, the condition measures the distance between two iterates. The absolute termination criteria on x works if the user has an accurate idea of the scale of the optimum x. If the optimum x is near 0, the relative tolerance will not work and the absolute tolerance is more appropriate.
  • The tolerance on simplex size is examined depending on the value of the tolsimplexizemethod element.
    • FALSE: then the criteria is just ignored,
    • TRUE: if \(ssize < tolsimplexizerelative \cdot simplexsize0 + tolsimplexizeabsolute\), where \(simplexsize0\) is the size of the simplex at iteration 0, then the status is set to ‘tolsize’ and terminate is set to TRUE.
    The size of the simplex is computed from the ‘sigmaplus’ method of the optimsimplex package. This criteria is sensitive to the tolsimplexizeabsolute and the tolsimplexizerelative elements.
  • The absolute tolerance on simplex size and absolute difference of function value is examined depending on the value of the tolssizedeltafvmethod element.
    • FALSE: then the criteria is just ignored,
    • TRUE: if both the following conditions \(ssize < tolsimplexizeabsolute\) and \(shiftfv < toldeltafv\) are true where \(ssize\) is the current simplex size and \(shiftfv\) is the absolute value of the difference of function value between the highest and lowest vertices, then the status is set to ‘tolsizedeltafv’ and terminate is set to TRUE.
  • The stagnation condition based on Kelley sufficient decrease condition is examined depending on the value of the kelleystagnationflag element.
    • FALSE: then the criteria is just ignored,
    • TRUE: if \(newfvmean \le oldfvmean - alpha \cdot t(sg) \cdot sg\) where \(newfvmean\) (resp. \(oldfvmean\)) is the function value average in the current iteration (resp. in the previous iteration), then the status is set to ‘kelleystagnation’ and terminate is set to TRUE. Here, \(alpha\) is a non-dimensional coefficient and \(sg\) is the simplex gradient.
  • The termination condition suggested by Box is examined depending on the value of the boxtermination element.
    • FALSE: then the criteria is just ignored,
    • TRUE: if both the following conditions \(shiftfv < boxtolf\) and \(boxkount == boxnbmatch\) are true, where \(shiftfv\) is the difference of function value between the best and worst vertices, and boxkount is the number of consecutive iterations where this criteria is met, then the status is set to ‘tolboxf’ and terminate is set to TRUE. Here, the boxtolf parameter is the value associated with the boxtolf element of the neldermead object and is a user-defined absolute tolerance on the function value. The boxnbmatch parameter is the value associated with the boxnbmatch element and is the user-defined number of consecutive match.
  • The termination condition based on the variance of the function values in the simplex is examined depending on the value of the tolvarianceflag element.
    • FALSE: then the criteria is just ignored,
    • TRUE: if \(var < tolrelativevariance \cdot variancesimplex0 + tolabsolutevariance\), where \(var\) is the variance of the function values in the simplex, then the status is set to ‘tolvariance’ and terminate is set to TRUE. Here, the tolrelativevariance parameter is the value associated with the tolrelativevariance element of the neldermead object and is a user-defined relative tolerance on the variance of the function values. The tolabsolutevariance parameter is the value associated with the tolabsolutevariance element and is the user-defined absolute tolerance of the variance of the function values.
  • The user-defined termination condition is examined depending on the value of the myterminateflag element.
    • FALSE: then the criteria is just ignored,
    • TRUE: if the term boolean output argument returned by the termination function is TRUE, then the status is set to the user-defined status and terminate is set to TRUE.

Notes about optimization methods

Kelley’s stagnation detection

The stagnation detection criteria suggested by Kelley is based on a sufficient decrease condition, which requires a parameter alpha > 0 to be defined [2]. The kelleynormalizationflag element of the neldermead object allows to configure the method to use to compute this alpha parameter. Two methods are available, where each method corresponds to a different paper by Kelley:

  • constant: in ‘Detection and Remediation of Stagnation in the Nelder-Mead Algorithm Using a Sufficient Decrease Condition,’ Kelley uses a constant alpha, with the suggested value 1.e-4, which is the typical choice for line search method.
  • normalized: in ‘Iterative Methods for Optimization,’ Kelley uses a normalized alpha, computed from the following formula: \(alpha = alpha0 \cdot sigma0 / nsg\), where \(sigma0\) is the size of the initial simplex and \(nsg\) is the norm of the simplex gradient for the initial guess point.

O’Neill’s factorial optimality test

In ‘Algorithm AS47 - Function minimization using a simplex procedure,’ O’Neill presents a fortran 77 implementation of the simplex method [3]. A factorial test is used to check if the computed optimum point is a local minimum. If the restartdetection element of the neldermead object is set to ‘oneill,’ that factorial test is used to see if a restart should be performed.

Method of Spendley et al.

The original paper may be implemented with several variations, which might lead to different results [4]. This section defines what algorithmic choices have been used in the present package.

The paper states the following rules. * ‘Rule 1. Ascertain the lowest reading y, of yi … yk+1 Complete a new simplex Sp by excluding the point Vp corresponding to y, and replacing it by V* defined as above.’ * ‘Rule 2. If a result has occurred in (k + 1) successive simplexes, and is not then eliminated by application of Rule 1, do not move in the direction indicated by Rule 1, or at all, but discard the result and replace it by a new observation at the same point.’ * ‘Rule 3. If y is the lowest reading in So , and if the next observation made, y* , is the lowest reading in the new simplex S , do not apply Rule 1 and return to So from Sp . Move out of S, by rejecting the second lowest reading (which is also the second lowest reading in So).’

We implement the following ‘rules’ of the Spendley et al. method:

  • Rule 1 is strictly applied, but the reflection is done by reflection of the high point, since we minimize a function instead of maximizing it, like Spendley.
  • Rule 2 is NOT implemented, as we expect that the function evaluation is not subject to errors.
  • Rule 3 is applied, i.e. reflection with respect to next to the high point. The original paper does not mention any shrink step. When the original algorithm cannot improve the function value with reflection steps, the basic algorithm stops. In order to make the current implementation of practical value, a shrink step is included, with shrinkage factor sigma. This perfectly fits into to the spirit of the original paper. Notice that the shrink step makes the rule #3 (reflection with respect to next-to-worst vertex) unnecessary. Indeed, the minimum required steps are the reflection and shrinkage. Nevertheless, the rule #3 has been kept in order to make the algorithm as close as it can be to the original.

Method of Nelder and Mead

The purpose of this section is to analyse the current implementation of Nelder-Mead’s algorithm. The algorithm that we use is described in ‘Iterative Methods for Optimization’ by Kelley.

The original paper uses a ‘greedy’ expansion, in which the expansion point is accepted whatever its function value. The current implementation, as most implementations, uses the expansion point only if it improves over the reflection point, that is,

  • if fe<fr, then the expansion point is accepted,
  • if not, the reflection point is accepted.

The termination criteria suggested by Nelder and Mead is based on an absolute tolerance on the standard deviation of the function values in the simplex. We provide this original termination criteria with the tolvarianceflag element of the neldermead object, which is disabled by default.

Box’s complex algorithm

In this section, we analyse the current implementation of Box’s complex method [5]. The initial simplex can be computed as in Box’s paper, but this may not be safe. In his paper, Box suggests that if a vertex of the initial simplex does not satisfy the non linear constraints, then it should be ‘moved halfway toward the centroid of those points already selected.’ This behaviour is available when the scalingsimplex0 element of the neldermead object is set to ‘tocenter.’ It may happen, as suggested by Guin [6], that the centroid is not feasible if the constraints are not convex. In this case, the initial simplex cannot be computed. This is why we provide the ‘tox0’ option, which allows to compute the initial simplex by scaling toward the initial guess, which is always feasible.

In Box’s paper, the scaling into the non linear constraints is performed ‘toward’ the centroid, that is, by using a scaling factor equal to 0.5. This default scaling factor might be sub-optimal in certain situations. This is why we provide the boxineqscaling element, which allows to configure the scaling factor.

In Box’s paper, whether we are concerned with the initial simplex or with the simplex at a given iteration, the scaling for the non linear constraints is performed without end. This is because Box’s hypothesis is that ‘ultimately, a satisfactory point will be found.’ As suggested by Guin, if the process fails, the algorithm goes into an infinite loop. In order to avoid this, we perform the scaling until a minimum scaling value is reached, as defined by the guinalphamin element.

We have taken into account the comments by Guin, but it should be emphasized that the current implementation is still as close as possible to Box’s algorithm and is not Guin’s algorithm. More precisely, during the iterations, the scaling for the non linear constraints is still performed toward the centroid, be it feasible or not.

User-defined algorithm

The mymethod element of the neldermead object allows to configure a user-defined simplex-based algorithm. The reason for this option is that many simplex-based variants of Nelder-Mead’s algorithm have been developed over the years, with specific goals. While it is not possible to provide them all, it is very convenient to use the current structure without being forced to make many developments.

The value of the mymethod element is expected to be a R function with the following structure:

  myalgorithm <- function( this ){
    ...
    return(this)
  }

where this is the current neldermead object.

In order to use the user-defined algorithm, the method element must be set to ‘mine.’ In this case, the component performs the optimization exactly as if the user-defined algorithm was provided by the component.

The user interested in that feature may use the internal scripts provided in the distribution as templates and tune his own algorithm from that point. There is of course no warranty that the user-defined algorithm improves on the standard algorithm, so that users use this feature at their own risks.

User-defined termination

Many termination criteria are found in the literature. Users who aim at reproducing the results exhibited in a particular paper may find that that none of the provided termination criteria match the one which is used in the paper. It may also happen that the provided termination criteria are not suitable for the specific test case. In those situation the myterminate element of the neldermead object allows to configure a user-defined termination function. The value of the myterminate element is expected to be a R function with the following structure:

  mystoppingrule <- function( this , simplex ){
    ...
    return(list(this = this, terminate = terminate, status = status))
  }

where this is the current neldermead object and simplex is the current simplex. The terminate output argument is a logical flag which is FALSE if the algorithm must continue and TRUE if the algorithm must stop. The status output argument is a string which is associated with the current termination criteria.

In order to enable the use of the user-defined termination function, the value of the myterminateflag element must be set to TRUE in the neldermead object. At each iteration, if the myterminateflag element has been set to TRUE, the user-defined termination is called. If the terminate output argument is TRUE, then the algorithm is stopped. In that case, the value of the status element of the neldermead.get function output is the value of the `status} output argument of the user-defined termination function.

Specialized functions

fminsearch

The fminsearch function is based on a specialized use of the more general neldermead function bundle and searches for the unconstrained minimum of a given cost function. This function corresponds to the Matlab (or Scilab) fminsearch function. In the context of fminsearch, the function to be minimized is not a cost function as described above but an objective function (returning a numeric scalar). Additional information and examples are available in ?fminsearch from a R environment.

Examples

We present in this section basic examples illustrating the use of neldermead functions to optimize unconstrained or constrained systems. More complex examples are described in a Scilab-based document written by Michael Baudin and available at http://forge.scilab.org/index.php/p/docneldermead/. Because the R port of the Scilab neldermead module is almost literal, the user should be able to reproduce the described examples in R with minimal adaptations.

Example 1: Basic use

In the following example, we solve a simple quadratic test case. We begin by defining the cost function, which takes 3 input arguments and returns the value of the objective function as the \(f\) element of a list. The standard starting point [-1.2 1.0] is used. neldermead creates a new neldermead object. Then we use neldermead.set to configure the parameters of the problem. We use all default settings and perform the search for the optimum. neldermead.get is finally used to retrieve the optimum parameters.

  quadratic <- function(x = NULL,index = NULL,fmsfundata = NULL){
    return(list(f = x[1]^2 + x[2]^2,
                g = c(),
                c = c(), 
                gc = c(),
                index = index,
                this = list(costfargument = fmsfundata)))
  }

  x0 <- transpose( c(1.0,1.0) )
  nm <- neldermead()
  nm <- neldermead.set(nm,'numberofvariables',2)
  nm <- neldermead.set(nm,'function',quadratic)
  nm <- neldermead.set(nm,'x0',x0)
  nm <- neldermead.search(nm)
  summary(nm)
## 
## Number of Estimated Variable(s): 2
## 
## Estimated Variable(s):
##   Initial         Final
## 1       1 -1.010582e-08
## 2       1 -1.768891e-07
## 
## Cost Function:
## function(x = NULL,index = NULL,fmsfundata = NULL){
##     return(list(f = x[1]^2 + x[2]^2,
##                 g = c(),
##                 c = c(), 
##                 gc = c(),
##                 index = index,
##                 this = list(costfargument = fmsfundata)))
##   }
## <bytecode: 0x558e34ed3e48>
## 
## Cost Function Argument(s):
## [1] ""
## 
## Optimization:
## - Status: "maxfuneval"
## - Initial Cost Function Value: 2.000000
## - Final Cost Function Value: 0.000000
## - Number of Iterations (max): 52 (100)
## - Number of Function Evaluations (max): 100 (100)
## 
## Simplex Information:
## 
## - Simplex at Initial Point:
## Dimension: n=2
## Number of vertices: nbve=3
##   Vertex #1/3 : fv=2.000000e+00, x=1.000000e+00 1.000000e+00
##   Vertex #2/3 : fv=5.000000e+00, x=2.000000e+00 1.000000e+00
##   Vertex #3/3 : fv=5.000000e+00, x=1.000000e+00 2.000000e+00
## 
## - Simplex at Optimal Point:
## Dimension: n=2
## Number of vertices: nbve=3
##   Vertex #1/3 : fv=3.139189e-14, x=-1.010582e-08 -1.768891e-07
##   Vertex #2/3 : fv=1.290894e-13, x=-3.557479e-07 5.032676e-08
##   Vertex #3/3 : fv=1.601186e-13, x=2.637847e-07 3.008924e-07

Example 2: Customized use

In the following example, we solve the Rosenbrock test case. We begin by defining the Rosenbrock function, which takes 3 input arguments and returns the value of the objective function. The standard starting point [-1.2 1.0] is used. neldermead creates a new neldermead object. Then we use neldermead.set to configure the parameters of the problem. The initial simplex is computed from the axes and the single length 1.0 (this is the default, but is explicitely written here as an example). The variable simplex algorithm by Nelder and Mead is used, which corresponds to the -method ‘variable’ option. neldermead.search performs the search for the minimum. Once the minimum is found, we represent part of the search space using the contour} function (this is possible since our problem involves only 2 parameters) and we superimpose the starting point (in red), the optimisation path (in blue), and the optimum (in green) to the plot. The history of the optimisation can be retrieved (usingneldermead.get`) because the ‘storehistory’ option was set to TRUE.

  rosenbrock <- function(x = NULL,index = NULL,fmsfundata = NULL){
    return(list(f = 100*(x[2]-x[1]^2)^2+(1-x[1])^2,
                g = c(),
                c = c(),
                gc = c(),
                index = index,
                this = list(costfargument = fmsfundata)))
  }
  x0 <- transpose(c(-1.2,1.0))
  nm <- neldermead()
  nm <- neldermead.set(nm,'numberofvariables',2)
  nm <- neldermead.set(nm,'function',rosenbrock)
  nm <- neldermead.set(nm,'x0',x0)
  nm <- neldermead.set(nm,'maxiter',200)
  nm <- neldermead.set(nm,'maxfunevals',300)
  nm <- neldermead.set(nm,'tolfunrelative',10*.Machine$double.eps)
  nm <- neldermead.set(nm,'tolxrelative',10*.Machine$double.eps)
  nm <- neldermead.set(nm,'simplex0method','axes')
  nm <- neldermead.set(nm,'simplex0length',1.0)
  nm <- neldermead.set(nm,'method','variable')
  nm <- neldermead.set(nm,'verbose',FALSE)
  nm <- neldermead.set(nm,'storehistory',TRUE)
  nm <- neldermead.set(nm,'verbosetermination',FALSE)
  nm <- neldermead.search(nm)
  
  xmin <- ymin <- -2.0 
  xmax <- ymax <- 2.0 
  nx <- ny <- 100
  stepy <- stepx <- (xmax - xmin)/nx
  ydata <- xdata <- seq(xmin,xmax,stepx)
  zdata <- apply(expand.grid(xdata,ydata),1,
                 function(x) neldermead.function(nm,transpose(x)))
  zdata <- matrix(zdata,ncol = length(ydata))
  optimpath <- matrix(unlist((neldermead.get(nm,'historyxopt'))),
                      nrow = 2)
  optimpath <- data.frame(x = optimpath[1,],y = optimpath[2,])

  contour(xdata,ydata,zdata,levels = c(1,10,100,500,1000,2000))
  par(new = TRUE,ann = TRUE)
  plot(c(x0[1],optimpath$x[158]), c(x0[2],optimpath$y[158]),
       col = c('red','green'),pch = 16,xlab = 'x[1]',ylab = 'x[2]',
       xlim = c(xmin,xmax),ylim = c(ymin,ymax))
  par(new = TRUE,ann = FALSE)  
  plot(optimpath$x,optimpath$y,col = 'blue',type = 'l',
       xlim = c(xmin,xmax),ylim = c(ymin,ymax))

Setting the ‘verbose’ element of the neldermead object to 1 allows to get detailed information about the current optimization process. The following is a sample output for an optimization based on the Nelder and Mead variable-shape simplex algorithm. Only the output corresponding to the iteration #156 is displayed. In order to display specific outputs (or to create specific output files and graphics), the ‘outputcommand’ option should be used.

 == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == = 
Iteration \#156 (total = 156)
Function Eval \#298                                                             
Xopt: 0.99999999999991 0.999999999999816                                        
Fopt: 8.997809e-27                                                              
DeltaFv: 4.492261e-26                                                           
Center: 1.00000000000003 1.00000000000007                                       
Size: 4.814034e-13                                                              
Vertex \#2/3 : fv = 2.649074e-26, x = 1.000000e+00 1.000000e+00               
Vertex \#3/3 : fv = 5.392042e-26, x = 1.000000e+00 1.000000e+00
Reflect                     
xbar = 1.00000000000001 1.00000000000003                                          
Function Evaluation \#299 at [0.99999999999996  ]
Function Evaluation \#299 at [0.999999999999907  ]    
xr = [0.99999999999996 0.999999999999907], f(xr) = 0.000000                         
  > Perform reflection                                                          
Sort 

Example 3: Optimization with bound constraints

In the following example, we solve a simple quadratic test case used in Example 1 but in the case where bounds are set for parameter estimates. We begin by defining the cost function, which takes 3 input arguments and returns the value of the objective function as the f element of a list. The starting point [1.2 1.9] is used. neldermead creates a new neldermead object. Then we use neldermead.set to configure the parameters of the problem including the lower -boundsmin and upper -boundsmax bounds. The initial simplex is computed from boxnbpoints random points within the bounds. The variable simplex algorithm by Box is used, which corresponds to the -method ‘box’ option. neldermead.search finally performs the search for the minimum.

  quadratic <- function(x = NULL,index = NULL,fmsfundata = NULL){
    return(list(f = x[1]^2 + x[2]^2,
                g = c(),
                c = c(), 
                gc = c(),
                index = index,
                this = list(costfargument = fmsfundata)))
  }
  set.seed(0)
  x0 <- transpose(c(1.2,1.9))
  nm <- neldermead()
  nm <- neldermead.set(nm,'numberofvariables',2)
  nm <- neldermead.set(nm,'function',quadratic)
  nm <- neldermead.set(nm,'x0',x0)
  nm <- neldermead.set(nm,'verbose',FALSE)
  nm <- neldermead.set(nm,'storehistory',TRUE)
  nm <- neldermead.set(nm,'verbosetermination',FALSE)
  nm <- neldermead.set(nm,'method','box')
  nm <- neldermead.set(nm,'boundsmin',c(1,1))
  nm <- neldermead.set(nm,'boundsmax',c(2,2))
  nm <- neldermead.search(nm)
  summary(nm)
## 
## Number of Estimated Variable(s): 2
## 
## Estimated Variable(s):
##   Initial    Final Lower bound Upper bound
## 1     1.2 1.000001           1           2
## 2     1.9 1.000001           1           2
## 
## Cost Function:
## function(x = NULL,index = NULL,fmsfundata = NULL){
##     return(list(f = x[1]^2 + x[2]^2,
##                 g = c(),
##                 c = c(), 
##                 gc = c(),
##                 index = index,
##                 this = list(costfargument = fmsfundata)))
##   }
## <bytecode: 0x558e35b24908>
## 
## Cost Function Argument(s):
## [1] ""
## 
## Optimization:
## - Status: "maxfuneval"
## - Initial Cost Function Value: 5.050000
## - Final Cost Function Value: 2.000004
## - Number of Iterations (max): 90 (100)
## - Number of Function Evaluations (max): 100 (100)
## 
## Simplex Information:
## 
## - Simplex at Initial Point:
## Dimension: n=2
## Number of vertices: nbve=3
##   Vertex #1/3 : fv=5.050000e+00, x=1.200000e+00 1.900000e+00
##   Vertex #2/3 : fv=7.610000e+00, x=2.000000e+00 1.900000e+00
##   Vertex #3/3 : fv=5.440000e+00, x=1.200000e+00 2.000000e+00
## 
## - Simplex at Optimal Point:
## Dimension: n=2
## Number of vertices: nbve=3
##   Vertex #1/3 : fv=2.000004e+00, x=1.000001e+00 1.000001e+00
##   Vertex #2/3 : fv=2.000004e+00, x=1.000001e+00 1.000001e+00
##   Vertex #3/3 : fv=2.000004e+00, x=1.000001e+00 1.000001e+00

Example 4: Optimization with nonlinear inequality constraints

In the following example, we solve Michalewicz’s \(G_6\) test problem using Box’s methods [7] (example suggested by Pascal Grandeau}. This problem consists in minimizing: \(G_{6}(x) = (x_{1}-10)^3+(x_{2}-20)^3\), given the nonlinear constraints:

\[ \begin{array}{l l} c1: & (x_{1}-5)^2+(x_{2}-5)^2-100 \ge{} 0 \\ c2: & -(x_{1}-6)^2-(x_{2}-5)^2+82.81 \ge{} 0 \\ \end{array} \]

and bounds: \(13 \le{} x_{1} \le{} 100,\: 0 \le{} x_{2} \le{} 100\).

We begin by defining the michalewicz function, which takes 3 input arguments and return the value of the objective function and the constraint evaluations as the f and c elements of a list. neldermead creates a new neldermead object. Then we use neldermead.set to configure the parameters of the problem, including the lower -boundsmin and upper -boundsmax bounds. The initial simplex is computed from boxnbpoints random points within the bounds. The variable simplex algorithm by Box is used, which corresponds to the -method ‘box’ option. neldermead.search finally performs the search for the minimum. The starting point ([15 4.99]) like all the vertices of the optimization simplex must be feasible, i.e. they must satisfy all constraints and bounds. Constraints are enforced by ensuring that all arguments of c in the cost function output are positive or null. Note that the boundaries were set to stricter ranges to limit the sensitivity of the solution to the initial guesses.

  michalewicz <- function(x = NULL,index = NULL,fmsfundata = NULL){
    f <- c()
    c <- c()
    if (index == 2 | index == 6) 
      f <- (x[1]-10)^3+(x[2]-20)^3
    
    if (index == 5 | index == 6)
      c <- c((x[1]-5)^2+(x[2]-5)^2 -100, 
          82.81-((x[1]-6)^2+(x[2]-5)^2))
    varargout <- list(f = f,
        g = c(),
        c = c, 
        gc = c(),
        index = index,
        this = list(costfargument = fmsfundata))
    return(varargout)
  }
  set.seed(0)
  x0 <- transpose(c(15,4.99))
  nm <- neldermead()
  nm <- neldermead.set(nm,'numberofvariables',2)
  nm <- neldermead.set(nm,'nbineqconst',2)
  nm <- neldermead.set(nm,'function',michalewicz)
  nm <- neldermead.set(nm,'x0',x0)
  nm <- neldermead.set(nm,'maxiter',300)
  nm <- neldermead.set(nm,'maxfunevals',1000)
  nm <- neldermead.set(nm,'simplex0method','randbounds')
  nm <- neldermead.set(nm,'boxnbpoints',3)
  nm <- neldermead.set(nm,'storehistory',TRUE)
  nm <- neldermead.set(nm,'method','box')
  nm <- neldermead.set(nm,'boundsmin',c(13,0))
  nm <- neldermead.set(nm,'boundsmax',c(20,10))
  nm <- neldermead.search(nm)
  summary(nm)
## 
## Number of Estimated Variable(s): 2
## 
## Estimated Variable(s):
##   Initial      Final Lower bound Upper bound
## 1   15.00 14.0950000          13          20
## 2    4.99  0.8429608           0          10
## 
## Number of Inequality Contraints: 2
## 
## Cost Function:
## function(x = NULL,index = NULL,fmsfundata = NULL){
##     f <- c()
##     c <- c()
##     if (index == 2 | index == 6) 
##       f <- (x[1]-10)^3+(x[2]-20)^3
##     
##     if (index == 5 | index == 6)
##       c <- c((x[1]-5)^2+(x[2]-5)^2 -100, 
##           82.81-((x[1]-6)^2+(x[2]-5)^2))
##     varargout <- list(f = f,
##         g = c(),
##         c = c, 
##         gc = c(),
##         index = index,
##         this = list(costfargument = fmsfundata))
##     return(varargout)
##   }
## <bytecode: 0x558e33d3bb20>
## 
## Cost Function Argument(s):
## [1] ""
## 
## Optimization:
## - Status: "impossibleimprovement"
## - Initial Cost Function Value: -3256.754501
## - Final Cost Function Value: -6961.813876
## - Number of Iterations (max): 236 (300)
## - Number of Function Evaluations (max): 794 (1000)
## 
## Simplex Information:
## 
## - Simplex at Initial Point:
## Dimension: n=2
## Number of vertices: nbve=3
##   Vertex #1/3 : fv=-3.256755e+03, x=1.500000e+01 4.990000e+00
##   Vertex #2/3 : fv=-3.276394e+03, x=1.506683e+01 4.953517e+00
##   Vertex #3/3 : fv=-3.188984e+03, x=1.507561e+01 5.082317e+00
## 
## - Simplex at Optimal Point:
## Dimension: n=2
## Number of vertices: nbve=3
##   Vertex #1/3 : fv=-6.961814e+03, x=1.409500e+01 8.429608e-01
##   Vertex #2/3 : fv=-6.961814e+03, x=1.409500e+01 8.429608e-01
##   Vertex #3/3 : fv=-6.961814e+03, x=1.409500e+01 8.429608e-01

Example 5: Passing data to the cost function

In the following example, we use a simple example to illustrate how to pass user-defined arguments to a user-defined cost function. We try to find the mean and standard deviation of some normally distributed data using maximum likelihood (actually a modified negative log-likelihood approach) (example suggested by Mark Taper).

We begin by defining the negLL function, which takes 3 input arguments and return the value of the objective function. The random dataset is then generated and stored in the list fmsdundata. neldermead creates a new neldermead object. Then we use neldermead.set to configure the parameters of the problem, including costfargument, set to fmsdundata, and the lower -boundsmin and upper -boundsmax bounds (the standard deviations has to be positive). The variable simplex algorithm by Box is used. neldermead.search finally performs the search for the minimum.

  negLL <- function(x = NULL, index = NULL, fmsfundata = NULL){
    mn <- x[1]
    sdv <- x[2]
    out <- -sum(dnorm(fmsfundata$data, mean = mn, sd = sdv, log = TRUE))
  
    return(list(f = out, 
           index = index,
           this = list(costfargument = fmsfundata)))
  }

  set.seed(12345)
  fmsfundata <- structure(
    list(data = rnorm(500,mean = 50,sd = 2)),
    class = 'optimbase.functionargs')

  x0 <- transpose(c(45,3))
  nm <- neldermead()
  nm <- neldermead.set(nm,'numberofvariables',2)
  nm <- neldermead.set(nm,'function',negLL)
  nm <- neldermead.set(nm,'x0',x0)
  nm <- neldermead.set(nm,'costfargument',fmsfundata)
  nm <- neldermead.set(nm,'maxiter',500)
  nm <- neldermead.set(nm,'maxfunevals',1500)
  nm <- neldermead.set(nm,'method','box')
  nm <- neldermead.set(nm,'storehistory',TRUE)
  nm <- neldermead.set(nm,'boundsmin',c(-100, 0))
  nm <- neldermead.set(nm,'boundsmax',c(100, 100))
  nm <- neldermead.search(this = nm)
  summary(nm)
## 
## Number of Estimated Variable(s): 2
## 
## Estimated Variable(s):
##   Initial     Final Lower bound Upper bound
## 1      45 50.164922        -100         100
## 2       3  1.978316           0         100
## 
## Cost Function:
## function(x = NULL, index = NULL, fmsfundata = NULL){
##     mn <- x[1]
##     sdv <- x[2]
##     out <- -sum(dnorm(fmsfundata$data, mean = mn, sd = sdv, log = TRUE))
##   
##     return(list(f = out, 
##            index = index,
##            this = list(costfargument = fmsfundata)))
##   }
## <bytecode: 0x558e2f4bb998>
## 
## Cost Function Argument(s):
## $data
##   [1] 51.17106 51.41893 49.78139 49.09301 51.21177 46.36409 51.26020
##   [8] 49.44763 49.43168 48.16136 49.76750 53.63462 50.74126 51.04043
##  [15] 48.49894 51.63380 48.22728 49.33684 52.24143 50.59745 51.55924
##  [22] 52.91157 48.71134 46.89373 46.80458 53.61020 49.03671 51.24076
##  [29] 51.22425 49.67538 51.62375 54.39367 54.09838 53.26489 50.50854
##  [36] 50.98238 49.35183 46.67590 53.53547 50.05160 52.25702 45.23928
##  [43] 47.87947 51.87428 51.70890 52.92146 47.17380 51.13481 51.16638
##  [50] 47.38640 48.91923 53.89539 50.10718 50.70333 48.65805 50.55591
##  [57] 51.38234 51.64759 54.29013 45.30611 50.29918 47.31494 51.10661
##  [64] 53.17993 48.82624 46.33525 51.77628 53.18698 51.03371 47.40866
##  [71] 50.10923 48.43070 47.90129 54.66102 52.80541 51.88520 51.65252
##  [78] 48.37692 50.95250 52.04252 51.29077 52.08629 49.39126 54.95422
##  [85] 51.94244 53.73420 51.34408 49.38409 51.07305 51.64974 48.07220
##  [92] 48.28983 53.77389 49.21636 48.03873 51.37466 48.98991 54.31544
##  [99] 48.80040 48.61091 50.44785 47.68755 50.84484 47.35049 50.28217
## [106] 48.92790 49.37679 53.11222 49.10393 50.64225 47.53966 47.35188
## [113] 52.52248 52.63846 49.83849 48.98982 49.89569 51.25772 54.36000
## [120] 49.86197 53.08973 52.64290 50.64430 53.06191 49.15752 47.68236
## [127] 46.30926 52.31465 45.75290 47.60794 53.28438 51.76731 51.04975
## [134] 47.63068 55.31158 47.90417 47.97775 51.33784 50.25835 49.15485
## [141] 47.71947 47.41257 48.81060 46.99837 50.03171 51.08034 46.90542
## [148] 51.69931 51.79203 50.27738 46.76134 51.09680 50.39056 48.38700
## [155] 49.78275 49.49811 53.39869 49.31140 50.13554 48.69886 49.02472
## [162] 50.60630 49.51605 49.03653 48.01639 49.43870 51.26603 47.52036
## [169] 53.52863 49.95264 50.39984 52.69439 50.07215 51.64916 46.59466
## [176] 50.96190 54.96710 50.80273 50.43035 46.36858 48.17652 49.90191
## [183] 49.18923 52.26076 51.63093 50.15284 52.90749 50.74824 49.65819
## [190] 48.99557 51.08704 48.98963 51.57359 50.60190 52.62045 51.59687
## [197] 51.70172 49.11286 49.10645 50.02661 47.12771 48.74148 50.48704
## [204] 52.11672 51.66270 50.21042 46.51657 51.29049 50.19421 49.84653
## [211] 51.98390 48.28150 49.43684 54.13249 48.77689 50.63123 51.32059
## [218] 46.55560 45.73075 50.13789 51.73564 45.41991 49.69962 49.46244
## [225] 53.58266 51.34454 49.58140 50.02437 53.06823 50.15458 50.15688
## [232] 48.44148 50.33312 50.53065 51.78156 49.06422 51.51675 48.71653
## [239] 51.25534 50.49666 48.59985 48.86520 49.47721 47.87223 49.78726
## [246] 51.54221 55.49481 49.83213 51.08714 51.50572 48.38265 52.00224
## [253] 50.91211 47.13150 49.46939 51.28354 49.16996 49.08085 48.41501
## [260] 47.68292 51.42178 52.53520 49.71370 48.96994 52.96578 49.67482
## [267] 50.08342 50.96608 47.63975 48.67285 48.73070 48.59607 51.15370
## [274] 45.77384 50.52182 52.29425 50.02959 49.37652 48.08761 50.94683
## [281] 46.97227 50.32856 48.25827 53.18666 51.29320 50.71474 50.20479
## [288] 48.64947 51.94417 51.51174 49.14343 48.57215 49.61923 50.79973
## [295] 48.04431 50.36747 45.69938 48.75407 48.46912 50.92862 51.04456
## [302] 50.01959 49.11895 52.39898 49.76506 50.07642 52.38961 50.68792
## [309] 49.34185 53.34172 48.16388 49.82439 52.64059 53.46157 54.32519
## [316] 49.36854 48.84981 47.18729 54.53572 48.45829 50.76063 51.21027
## [323] 52.03935 50.94989 45.62811 51.86638 50.95087 50.78056 48.54534
## [330] 51.97311 52.84797 50.96946 50.69847 51.72025 50.80922 50.73409
## [337] 46.96160 53.09961 50.99976 50.92175 54.15342 49.38470 51.90474
## [344] 51.06558 49.80890 49.71587 47.63395 51.09064 44.83629 51.55780
## [351] 50.58588 49.82658 47.06728 47.83364 52.11547 49.27936 50.70119
## [358] 50.05652 50.94610 48.16169 49.24835 46.37433 50.57720 49.62075
## [365] 50.03572 51.30086 50.62051 53.33671 51.34523 49.44496 49.70795
## [372] 53.40286 50.94272 51.16340 51.33204 48.44220 52.32665 46.07091
## [379] 51.53834 54.51954 49.04546 49.79484 50.73739 48.92913 51.01320
## [386] 49.69889 51.80849 54.48407 47.60975 49.16295 51.59650 50.99634
## [393] 50.23912 49.26559 50.52466 50.52541 51.28148 50.61418 49.93374
## [400] 47.25050 51.25593 50.00429 50.56876 47.99644 48.76556 51.65639
## [407] 49.83036 49.13056 47.59033 47.95863 48.05754 49.67817 52.40937
## [414] 50.13349 48.23508 49.96016 50.65603 48.13874 47.73717 52.17246
## [421] 51.02848 50.58634 48.87469 50.41429 52.45615 49.13997 52.66382
## [428] 48.47347 52.46464 51.21812 48.35492 49.90470 50.34136 47.54986
## [435] 50.77467 50.98134 54.47825 51.78442 51.59787 53.87089 50.78850
## [442] 51.26504 50.20518 52.02042 52.42936 51.00174 48.18106 45.60558
## [449] 50.87709 47.94117 51.42616 52.05654 48.85873 49.13198 52.01063
## [456] 52.40709 52.99237 49.68680 48.12794 50.70508 47.42762 49.45495
## [463] 50.31383 48.23637 49.75620 48.28950 47.98854 48.17210 45.84353
## [470] 49.71871 52.59928 49.94126 45.52676 51.09006 51.66168 49.00484
## [477] 48.79795 48.97560 49.99014 49.80570 49.95193 49.74522 49.73623
## [484] 50.90748 47.18733 50.42186 51.50467 50.27676 50.31936 50.45687
## [491] 47.67443 51.70932 49.51581 45.96801 49.54462 50.44834 48.52635
## [498] 48.30658 49.77708 49.37680
## 
## attr(,"class")
## [1] "optimbase.functionargs"
## 
## Optimization:
## - Status: "impossibleimprovement"
## - Initial Cost Function Value: 1858.501814
## - Final Cost Function Value: 1050.592365
## - Number of Iterations (max): 137 (500)
## - Number of Function Evaluations (max): 268 (1500)
## 
## Simplex Information:
## 
## - Simplex at Initial Point:
## Dimension: n=2
## Number of vertices: nbve=3
##   Vertex #1/3 : fv=1.858502e+03, x=4.500000e+01 3.000000e+00
##   Vertex #2/3 : fv=1.599340e+03, x=4.600000e+01 3.000000e+00
##   Vertex #3/3 : fv=1.630588e+03, x=4.500000e+01 4.000000e+00
## 
## - Simplex at Optimal Point:
## Dimension: n=2
## Number of vertices: nbve=3
##   Vertex #1/3 : fv=1.050592e+03, x=5.016492e+01 1.978316e+00
##   Vertex #2/3 : fv=1.050592e+03, x=5.016492e+01 1.978316e+00
##   Vertex #3/3 : fv=1.050592e+03, x=5.016492e+01 1.978316e+00

Dependencies of fminsearch

We illustrate in the figures below the network of functions of the neldermead, optimbase, and optimsimplex packages that are called from the fminsearch functions. This large network is broken down in 6 plots, which are shown in the order functions are called. Green boxes represent functions that are not expanded on a given plot but on a previous or later one.

fminsearch function network (1/6; click to view at full scale)

fminsearch function network (2/6; click to view at full scale)

fminsearch function network (3/6; click to view at full scale)

fminsearch function network (4/6; click to view at full scale)

fminsearch function network (5/6; click to view at full scale)

fminsearch function network (6/6; click to view at full scale)

References

1. J.A. Nelder and R. Mead. A Simplex Method for Function Minimization. The Computer Journal. 1965;7:308–13.
2. C.T. Kelley. Iterative Methods for Optimization. Philadelphia, PA: SIAM Frontiers in Applied Mathematics; 1999.
3. R. O’Neill. Algorithm AS47 - Function minimization using a simplex procedure. Applied Statistics. 1971;20:338–45.
4. W. Spendley and G.R. Hext and F.R. Himsworth. Sequential Application of Simplex Designs in Optimisation and Evolutionary Operation. Technometrics. 1962;4:441–61.
5. M.J. Box. A New Method of Constrained Optimization and a Comparison With Other Methods. The Computer Journal. 1965;1:42–52.
6. J.A. Guin. Discussion and correspondence: modification of the complex method of constrained optimization. The Computer Journal. 1968;10:416–7.
7. Z. Michalewicz and D.B. Fogel. How to solve it: modern heuristics. Springer; 2004. p. 231–70.