neldermead
packageneldermead
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.
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 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 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, andthis
: 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:
f
,c
,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 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…
The current package takes into account several generic termination criteria. The following termination criteria are enabled by default:
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.
maxiter
element of the neldermead
object: if iterations
\(\ge\) maxiter
, then the status is set to ‘maxiter’ and terminate is set to TRUE.maxfunevals
elements: if funevals
\(\ge\) maxfunevals
, then the status is set to ‘maxfuneval’ and terminate is set to TRUE.tolfunmethod
.
tolxmethod
element.
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.tolsimplexizemethod
element.
optimsimplex
package. This criteria is sensitive to the tolsimplexizeabsolute
and the tolsimplexizerelative
elements.tolssizedeltafvmethod
element.
kelleystagnationflag
element.
boxtermination
element.
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.tolvarianceflag
element.
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.myterminateflag
element.
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.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:
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.
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:
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,
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.
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.
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.
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.
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.
Direct grid search, performed by fmin.gridsearch
, is a functionality added to the original Scilab neldermead
module and constitutes another specialized use of the neldermead
package. This function allows to explore the search space of an optimization problem around the initial point \(x_0\). This optimization problem is defined by an objective function, like for fminsearch
, and not a cost function. fmin.gridsearch
automatically creates a grid of search points selected around the initial point and evaluates the objective function at each point. The boundaries of the grid are set either by a vector of parameter-specific lower and upper limits, or by a vector of factors \(\alpha\) as follows:\([x_{0}/\alpha,x_{0}\times\alpha]\). The number \(npts\) of points evaluated for each parameter (or dimension of the optimization problem) can also be defined. The total number of points in the grid is therefore \(npts^n\). At the end of the search, fmin.gridsearch
returns a table sorted by value of the objective function. The feasibility of the objective function is also determined at each point, as fmin.gridsearch
is a wrapper around optimbase.gridsearch
which assesses the feasibility of a cost function in addition to calculating its value at each particular search point. Because fmin.gridsearch
does not accept constraints, the objective function should always be feasible. Additional information is available in ?fmin.gridsearch
from a R environment.
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.
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
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 (using
neldermead.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
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
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
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
In the following example, we use the Rosenbrock test case introduced as Example 2 to illustrate the direct grid search capacity of neldermead
. We begin by defining the Rosenbrock function, which takes only one input argument and returns the value of the objective function. We request 6 points per dimension of the problem and set the range of search around the standard starting point [-1.2 1.0] by providing limits. fmin.gridsearch
performs the search and return a table sorted by value of the cost function.
rosenbrock <- function(x = NULL){
f <- 100*(x[2]-x[1]^2)^2+(1-x[1])^2
}
x0 <- c(-1.2,1.0)
npts <- 6
xmin <- c(-2,-2)
xmax <- c(2,2)
grid <- fmin.gridsearch(fun = rosenbrock,x0 = x0,xmin = xmin,xmax = xmax,npts = npts,alpha = alpha)
## The grid contains 30 unique combinations.
## Evaluating combination number: 1/30
## Evaluating combination number: 2/30
## Evaluating combination number: 3/30
## Evaluating combination number: 4/30
## Evaluating combination number: 5/30
## Evaluating combination number: 6/30
## Evaluating combination number: 7/30
## Evaluating combination number: 8/30
## Evaluating combination number: 9/30
## Evaluating combination number: 10/30
## Evaluating combination number: 11/30
## Evaluating combination number: 12/30
## Evaluating combination number: 13/30
## Evaluating combination number: 14/30
## Evaluating combination number: 15/30
## Evaluating combination number: 16/30
## Evaluating combination number: 17/30
## Evaluating combination number: 18/30
## Evaluating combination number: 19/30
## Evaluating combination number: 20/30
## Evaluating combination number: 21/30
## Evaluating combination number: 22/30
## Evaluating combination number: 23/30
## Evaluating combination number: 24/30
## Evaluating combination number: 25/30
## Evaluating combination number: 26/30
## Evaluating combination number: 27/30
## Evaluating combination number: 28/30
## Evaluating combination number: 29/30
## Evaluating combination number: 30/30
grid
## x1 x2 f feasible
## 22 1.0 1 0.0 1
## 15 0.0 0 1.0 1
## 20 -1.0 1 4.0 1
## 24 -1.2 1 24.2 1
## 30 -1.2 2 36.2 1
## 16 1.0 0 100.0 1
## 28 1.0 2 100.0 1
## 9 0.0 -1 101.0 1
## 21 0.0 1 101.0 1
## 14 -1.0 0 104.0 1
## 26 -1.0 2 104.0 1
## 18 -1.2 0 212.2 1
## 10 1.0 -1 400.0 1
## 3 0.0 -2 401.0 1
## 27 0.0 2 401.0 1
## 29 2.0 2 401.0 1
## 8 -1.0 -1 404.0 1
## 25 -2.0 2 409.0 1
## 12 -1.2 -1 600.2 1
## 4 1.0 -2 900.0 1
## 23 2.0 1 901.0 1
## 2 -1.0 -2 904.0 1
## 19 -2.0 1 909.0 1
## 6 -1.2 -2 1188.2 1
## 17 2.0 0 1601.0 1
## 13 -2.0 0 1609.0 1
## 11 2.0 -1 2501.0 1
## 7 -2.0 -1 2509.0 1
## 5 2.0 -2 3601.0 1
## 1 -2.0 -2 3609.0 1
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)