$LastChangedDate: 16-12-2008 $
optimsimplex
Manage a simplex with arbitrary number of points.
SYNOPSIS
[ newobj , data ] = optimsimplex_new ( coords , fun , data )
this = optimsimplex_destroy (this)
[ this , data ] = optimsimplex_axes ( this , x0 , fun , len , data )
[ this , data ] = optimsimplex_pfeffer ( this , x0 , fun , deltausual , deltazero , data )
[ this , data ] = optimsimplex_randbounds ( this , x0 , fun , boundsmin , boundsmax , nbpoints , data )
[ this , data ] = optimsimplex_spendley ( this , x0 , fun , len , data )
this = optimsimplex_setall ( this , simplex )
this = optimsimplex_setallfv ( this , fv )
this = optimsimplex_setallx ( this , x )
this = optimsimplex_setfv ( this , ive , fv )
this = optimsimplex_setn ( this , n )
this = optimsimplex_setnbve ( this , nbve )
this = optimsimplex_setve ( this , ive , fv , x )
this = optimsimplex_setx ( this , ive , x )
simplex = optimsimplex_getall ( this )
fv = optimsimplex_getallfv ( this )
x = optimsimplex_getallx ( this )
fv = optimsimplex_getfv ( this , ive )
n = optimsimplex_getn ( this )
n = optimsimplex_getnbve ( this )
vertex = optimsimplex_getve ( this , ive )
x = optimsimplex_getx ( this , ive )
sicenter = optimsimplex_center ( this )
optimsimplex_check ( this )
[ this , data ] = optimsimplex_computefv ( this , fun , data )
df = optimsimplex_deltafv ( this )
dfm = optimsimplex_deltafvmax ( this )
m = optimsimplex_dirmat ( this )
sd = optimsimplex_fvmean ( this )
sd = optimsimplex_fvstdev ( this )
[ g , data ] = optimsimplex_gradientfv ( this , fun , method , data )
g = optimsimplex_gradforward ( this )
[ g , data ] = optimsimplex_gradcentered ( this , fun , data )
[ ns , data ] = optimsimplex_oriented ( this , fun , data )
optimsimplex_print ( this )
[ r , data ] = optimsimplex_reflect ( this , fun , data )
[ this , data ] = optimsimplex_shrink ( this , fun , sigma , data )
ssize = optimsimplex_size ( this , method )
this = optimsimplex_sort ( this )
str = optimsimplex_tostring ( this )
cen = optimsimplex_xbar ( this , iexcl )
Purpose
The goal of this component is to provide a building block for
optimization algorithms based on a simplex. The optimsimplex package may
be used in the following optimization methods :
the Spendley et al. simplex method,
the Nelder-Mead method,
the Box algorithm for constrained optimization,
the multi-dimensional search by Virginia Torczon,
etc ...
Design
This toolbox is designed with Oriented Object ideas in mind.
Features
The following is a list of features the Nelder-Mead prototype
algorithm currently provides :
Manage various simplex initializations
initial simplex given by user,
initial simplex computed with a length and along the
coordinate axes,
initial regular simplex computed with Spendley et al.
formula,
initial simplex computed by a small perturbation around the
initial guess point,
initial simplex computed from randomized bounds.
sort the vertices by increasing function values,
compute the standard deviation of the function values in the
simplex,
compute the simplex gradient with forward or centered
differences,
shrink the simplex toward the best vertex,
etc...
Description
This set of commands allows to manage a simplex made of k>=n+1
points in a n-dimensional space. This component is the building block for
a class of direct search optimization methods such as the Nelder-Mead
algorithm or Torczon's Multi-Dimensionnal Search.
A simplex is designed as a collection of k>=n+1 vertices. Each
vertex is made of a point and a function value at that point.
The simplex can be created with various shapes. It can be configured
and quieried at will. The simplex can also be reflected or shrinked. The
simplex gradient can be computed with a order 1 forward formula and with a
order 2 centered formula.
The optimsimplex_new function allows to create a simplex. If vertices
coordinates are given, there are registered in the simplex. If a function
is provided, it is evaluated at each vertex. The optimsimplex_destroy function
destroys the object and frees internal memory. Several functions allow to
create a simplex with special shapes, including axes-by-axes
(optimsimplex_axes), regular (optimsimplex_spendley), randomized bounds simplex with
arbitrary k vertices (optimsimplex_randbounds) and an heuristical small
variation around a given point (optimsimplex_pfeffer).
Functions
The following functions are available.
[ newobj , data ] = optimsimplex_new ( coords , fun , data )
Creates a new simplex object. All input arguments are
optional. If no input argument is provided, returns an empty simplex
object.
newobj
The new object.
coords
optional, matrix of point coordinates in the
simplex.
The coords matrix is expected to be a n x k matrix,
where n is the dimension of the space and k is the number of
vertices in the simplex, with k>= n+1.
fun
optional, the function to compute at vertices.
The function is expected to have the following input and
output arguments :
function y = myfunction (x)
data
optional, user-defined data passed to the
function.
If data is provided, it is passed to the callback
function both as an input and output argument. In that case,
the function must have the following header :
function [ y , data ] = myfunction ( x , data )
The data input parameter may be used if the function
uses some additionnal parameters. It is returned as an output
parameter because the function may modify the data while
computing the function value. This feature may be used, for
example, to count the number of times that the function has
been called.
this = optimsimplex_destroy (this)
Destroy the given object.
this
The current simplex object.
[ this , data ] = optimsimplex_axes ( this , x0 , fun , len , data
)
Configure the current simplex so that it is computed axis by
axis, with the given length.
this
The current simplex object.
x0
the initial point, as a row vector.
fun
optional, the function to compute at vertices.
The function is expected to have the following input and
output arguments :
function y = myfunction (x)
len
optional, the length of the simplex. The default length
is 1.0. If length is a value, that unique length is used in
all directions. If length is a vector with n values, each
length is used with the corresponding direction.
data
optional, user-defined data passed to the
function.
If data is provided, it is passed to the callback
function both as an input and output argument. In that case,
the function must have the following header :
function [ y , data ] = myfunction ( x , data )
The data input parameter may be used if the function
uses some additionnal parameters. It is returned as an output
parameter because the function may modify the data while
computing the function value. This feature may be used, for
example, to count the number of times that the function has
been called.
[ this , data ] = optimsimplex_pfeffer ( this , x0 , fun , deltausual
, deltazero , data )
Configure the current simplex so that it is computed from
Pfeffer's method, i.e. a relative delta for non-zero values and an
absolute delta for zero values.
this
The current simplex object.
x0
the initial point, as a row vector.
fun
optional, the function to compute at vertices.
The function is expected to have the following input and
output arguments :
function y = myfunction (x)
deltausual
optional, the absolute delta for non-zero values. The
default value is 0.05.
deltazero
optional, the absolute delta for zero values. The
default value is 0.0075.
data
optional, user-defined data passed to the
function.
If data is provided, it is passed to the callback
function both as an input and output argument. In that case,
the function must have the following header :
function [ y , data ] = myfunction ( x , data )
The data input parameter may be used if the function
uses some additionnal parameters. It is returned as an output
parameter because the function may modify the data while
computing the function value. This feature may be used, for
example, to count the number of times that the function has
been called.
[ this , data ] = optimsimplex_randbounds ( this , x0 , fun ,
boundsmin , boundsmax , nbpoints , data )
Configure the current simplex so that it is computed by taking
the bounds into account with random scaling. The number of vertices
in the simplex is arbitrary.
this
The current simplex object.
x0
the initial point, as a row vector. It is the first
vertex in the simplex.
fun
optional, the function to compute at vertices.
The function is expected to have the following input and
output arguments :
function y = myfunction (x)
boundsmin
array of minimum bounds
boundsmax
array of maximum bounds
Each component ix =1 , n of the vertex #jve = 2,k is
computed from the formula :
x ( ix , jve ) = boundsmin( ix ) + rand() * (boundsmax( ix ) - boundsmin( ix ))
nbpoints
total number of points in the simplex
data
optional, user-defined data passed to the
function.
If data is provided, it is passed to the callback
function both as an input and output argument. In that case,
the function must have the following header :
function [ y , data ] = myfunction ( x , data )
The data input parameter may be used if the function
uses some additionnal parameters. It is returned as an output
parameter because the function may modify the data while
computing the function value. This feature may be used, for
example, to count the number of times that the function has
been called.
[ this , data ] = optimsimplex_spendley ( this , x0 , fun , len ,
data )
Configure the current simplex so that it is computed from
Spendley's et al. method, i.e. a regular simplex made of k = n+1
vertices.
this
The current simplex object.
x0
the initial point, as a row vector.
fun
optional, the function to compute at vertices.
The function is expected to have the following input and
output arguments :
function y = myfunction (x)
len
optional, the length of the simplex. The default length
is 1.0.
data
optional, user-defined data passed to the
function.
If data is provided, it is passed to the callback
function both as an input and output argument. In that case,
the function must have the following header :
function [ y , data ] = myfunction ( x , data )
The data input parameter may be used if the function
uses some additionnal parameters. It is returned as an output
parameter because the function may modify the data while
computing the function value. This feature may be used, for
example, to count the number of times that the function has
been called.
this = optimsimplex_setall ( this , simplex )
Set all the coordinates and and the function values of all the
vertices.
this
The current simplex object.
simplex
the simplex to set.
The given matrix is expected to be a nbve x n+1 matrix
where n is the dimension of the space, nbve is the number of
vertices and with the following content (where the data is
organized by row with function value first, and x
coordinates)
simplex(k,1) is the function value of the vertex #k,
with k = 1 , nbve
simplex(k,2:n+1) is the coordinates of the vertex
#k, with k = 1 , nbve
this = optimsimplex_setallfv ( this , fv )
Set all the function values of all the vertices. The vertex #k
is expected to be stored in fv(k) with k = 1 , nbve
this
The current simplex object.
fv
the array of function values
this = optimsimplex_setallx ( this , x )
Set all the coordinates of all the vertices. The vertex #k is
expected to be stored in x(1:n,k) with k = 1 , nbve
this
The current simplex object.
x
the coordinates of the vertices.
this = optimsimplex_setfv ( this , ive , fv )
Set the function value at given index and // returns an
updated simplex.
this
The current simplex object.
ive
vertex index
fv
the function value
this = optimsimplex_setn ( this , n )
Set the dimension of the space of the simplex.
this
The current simplex object.
n
the dimension
this = optimsimplex_setnbve ( this , nbve )
Set the number of vertices of the simplex.
this
The current simplex object.
nbve
the number of vertices
this = optimsimplex_setve ( this , ive , fv , x )
Sets the coordinates of the vertex and the function value at
given index in the current simplex.
this
The current simplex object.
ive
the vertex index
fv
the function value
x
the coordinates of the point, as a column vector
this = optimsimplex_setx ( this , ive , x )
Set the coordinates of the vertex at given index, as a column
vector, into the current simplex.
this
The current simplex object.
ive
the vertex index
x
the coordinates of the point, as a column vector
simplex = optimsimplex_getall ( this )
Returns all the coordinates of all the vertices and the
function values in the same matrix.
this
The current simplex object.
simplex
the simplex data.
The simplex matrix has size nbve x n+1, and is organized
by row by row as follows :
simplex(k,1) is the function value of the vertex #k,
with k = 1 , nbve
simplex(k,2:n+1) is the coordinates of the vertex
#k, with k = 1 , nbve
fv = optimsimplex_getallfv ( this )
Returns all the function values of all the vertices, as a
column vector.
this
The current simplex object.
fv
The array of function values. The function value of
vertex #k is stored in fv(k) with k = 1 , nbve.
x = optimsimplex_getallx ( this )
Returns all the coordinates of all the vertices.
this
The current simplex object.
x
the coordinates.
The vertex #k is stored in x(1:n,k) with k = 1 ,
nbve.
fv = optimsimplex_getfv ( this , ive )
Returns the function value at given index
this
The current simplex object.
ive
the vertex index
n = optimsimplex_getn ( this )
Returns the dimension of the space of the simplex
this
The current simplex object.
n = optimsimplex_getnbve ( this )
Returns the number of vertices in the simplex.
this
The current simplex object.
vertex = optimsimplex_getve ( this , ive )
Returns the vertex at given index as a tlist, with fields n, x
and fv
this
The current simplex object.
ive
the vertex index
x = optimsimplex_getx ( this , ive )
Returns the coordinates of the vertex at given index, as a
column vector.
this
The current simplex object.
ive
the vertex index
sicenter = optimsimplex_center ( this )
Returns the center of the given simplex
this
The current simplex object.
optimsimplex_check ( this )
Check the consistency of the internal data. Generates an error
if necessary.
this
The current simplex object.
[ this , data ] = optimsimplex_computefv ( this , fun , data )
Set the values of the function at vertices points.
this
The current simplex object.
fun
optional, the function to compute at vertices.
The function is expected to have the following input and
output arguments :
function y = myfunction (x)
data
optional, user-defined data passed to the
function.
If data is provided, it is passed to the callback
function both as an input and output argument. In that case,
the function must have the following header :
function [ y , data ] = myfunction ( x , data )
The data input parameter may be used if the function
uses some additionnal parameters. It is returned as an output
parameter because the function may modify the data while
computing the function value. This feature may be used, for
example, to count the number of times that the function has
been called.
df = optimsimplex_deltafv ( this )
Returns the vector of difference of function values with
respect to the function value at vertex #1.
this
The current simplex object.
dfm = optimsimplex_deltafvmax ( this )
Returns the difference of function value between the high and
the low vertices. It is expected that the vertex #1 is associated
with the smallest function value and that the vertex #nbve is
associated with the highest function value. Since vertices are
ordered, the high is greater than the low.
this
The current simplex object.
m = optimsimplex_dirmat ( this )
Returns the n x n matrix of simplex directions i.e. the matrix
of differences of vertices coordinates with respect to the vertex
#1.
this
The current simplex object.
sd = optimsimplex_fvmean ( this )
Returns the mean of the function value on the simplex.
this
The current simplex object.
sd = optimsimplex_fvstdev ( this )
Returns the standard deviation of the function value on the
simplex.
this
The current simplex object.
[ g , data ] = optimsimplex_gradientfv ( this , fun , method , data
)
Returns the simplex gradient of the function.
this
The current simplex object.
fun
optional, the function to compute at vertices.
The function is expected to have the following input and
output arguments :
function y = myfunction (x)
method
optional, the method to use to compute the simplex
gradient. Two methods are available : "forward" or "centered".
The forward method uses the current simplex to compute the
simplex gradient. The centered method creates an intermediate
reflected simplex and computes the average.
If not provided, the default method is "forward".
data
optional, user-defined data passed to the
function.
If data is provided, it is passed to the callback
function both as an input and output argument. In that case,
the function must have the following header :
function [ y , data ] = myfunction ( x , data )
The data input parameter may be used if the function
uses some additionnal parameters. It is returned as an output
parameter because the function may modify the data while
computing the function value. This feature may be used, for
example, to count the number of times that the function has
been called.
[ ns , data ] = optimsimplex_oriented ( this , fun , data )
Returns a new oriented simplex, in sorted order. The new
simplex has the same sigma- length of the base simplex, but is
"oriented" depending on the function value. This simplex may be
used, as Kelley suggests, for a restart of Nelder-Mead
algorithm.
this
The current simplex object.
fun
optional, the function to compute at vertices.
The function is expected to have the following input and
output arguments :
function y = myfunction (x)
data
optional, user-defined data passed to the
function.
If data is provided, it is passed to the callback
function both as an input and output argument. In that case,
the function must have the following header :
function [ y , data ] = myfunction ( x , data )
The data input parameter may be used if the function
uses some additionnal parameters. It is returned as an output
parameter because the function may modify the data while
computing the function value. This feature may be used, for
example, to count the number of times that the function has
been called.
optimsimplex_print ( this )
Display the current simplex, with coordinates and function
values.
this
The current simplex object.
[ r , data ] = optimsimplex_reflect ( this , fun , data )
Returns a new simplex by reflexion of current simplex, by
reflection with respect to the first vertex in the simplex. This
move is used in the centered simplex gradient.
this
The current simplex object.
fun
optional, the function to compute at vertices.
The function is expected to have the following input and
output arguments :
function y = myfunction (x)
data
optional, user-defined data passed to the
function.
If data is provided, it is passed to the callback
function both as an input and output argument. In that case,
the function must have the following header :
function [ y , data ] = myfunction ( x , data )
The data input parameter may be used if the function
uses some additionnal parameters. It is returned as an output
parameter because the function may modify the data while
computing the function value. This feature may be used, for
example, to count the number of times that the function has
been called.
[ this , data ] = optimsimplex_shrink ( this , fun , sigma , data
)
Shrink the simplex with given coefficient sigma and returns an
updated simplex. The shrink is performed with respect to the first
point in the simplex.
this
The current simplex object.
fun
optional, the function to compute at vertices.
The function is expected to have the following input and
output arguments :
function y = myfunction (x)
sigma
optional, the shrinkage coefficient. The default value
is 0.5.
data
optional, user-defined data passed to the
function.
If data is provided, it is passed to the callback
function both as an input and output argument. In that case,
the function must have the following header :
function [ y , data ] = myfunction ( x , data )
The data input parameter may be used if the function
uses some additionnal parameters. It is returned as an output
parameter because the function may modify the data while
computing the function value. This feature may be used, for
example, to count the number of times that the function has
been called.
ssize = optimsimplex_size ( this , method )
Returns the size of the simplex.
The diameter is the maximum length of all the vertices. It
requires 2 nested loops over the vertices.
The sigmaminus (resp. sigmamplus) size is the minimum (resp.
maximum) length of the vector from each vertex to the first
vertex. It requires one loop over the vertices.
The "Nash" size is the sum of the norm of the length of the
vector from the given vertex to the first vertex. It requires one
loop over the vertices.
this
The current simplex object.
method
optional, the method to use to compute the size. The
available methods are "Nash", "diameter", "sigmaplus" or
"sigmaminus" (default is "sigmaplus").
this = optimsimplex_sort ( this )
Sorts the simplex with increasing function value order so that
the smallest function value is at vertex #1
this
The current simplex object.
str = optimsimplex_tostring ( this )
Returns the current simplex as a string.
this
The current simplex object.
cen = optimsimplex_xbar ( this , iexcl )
Returns the center of n vertices, by excluding the vertex with
index iexcl.
this
The current simplex object.
iexcl
the index of the vertex to exclude in center
computation. The default value of iexcl is the number of
vertices : in that case, if the simplex is sorted in
increasing function value order, the worst vertex is
excluded.
Example : Creating a simplex with given vertices
coordinates
In the following example, one creates a simplex with known vertices
coordinates. The function values at the vertices are unset.
coords = [
0.0 1.0 0.0
0.0 0.0 1.0
];
s1 = optimsimplex_new ( coords );
computed = optimsimplex_getallx ( s1 );
computed = optimsimplex_getn(s1);
computed = optimsimplex_getnbve (s1);
s1 = optimsimplex_destroy(s1);
Example : Creating a simplex with randomized bounds
In the following example, one creates a simplex with in the 2D
domain [-5 5]^2, with [-1.2 1.0] as the first vertex. One uses the
randomized bounds method to generate a simplex with 5 vertices. The
function takes an additionnal argument myobj, which is counts the number
of times the function is called. After the creation of the simplex, the
value of mydude.nb is 5, which is the expected result because there is one
function call by vertex.
function [ y , myobj ] = mycostf ( x , myobj )
y = rosenbrock(x);
myobj.nb = myobj.nb + 1
endfunction
mydude = tlist(["T_MYSTUFF","nb"]);
mydude.nb = 0;
s1 = optimsimplex_new ();
[ s1 , mydude ] = optimsimplex_randbounds ( s1 , x0 = [-1.2 1.0], fun = mycostf, ...
boundsmin = [-5.0 -5.0] , boundsmax = [5.0 5.0], nbve=5 , data = mydude );
s1 = optimsimplex_destroy ( s1 );
TODO
implement reflection and expansion as in multidimensional
search by Torczon
Initial simplex strategies
In this section, we analyse the various initial simplex which are provided in this
component.
It is known that direct search methods based on simplex designs are very
sensitive to the initial simplex. This is why the current component
provides various ways to create such an initial simplex.
The first historical simplex-based algorithm is the one presented in
"Sequential Application of Simplex Designs in Optimisation and
Evolutionary Operation" by W. Spendley, G. R. Hext and F. R. Himsworth.
The optimsimplex_spendley function creates the regular simplex which is presented in
this paper.
The method used in the optimsimplex_randbounds function is due to M.J. Box
in "A New Method of Constrained Optimization and a Comparison With
Other Methods".
Pfeffer's method is an heuristic which is presented in "Global Optimization Of Lennard-Jones Atomic Clusters" by
Ellen Fan. It is due to L. Pfeffer at Stanford and it is used in fminsearch.
Authors
Michael Baudin, 2008-2009
Bibliography
“Sequential Application of Simplex Designs in Optimisation and
Evolutionary Operation”, Spendley, W. and Hext, G. R. and Himsworth,
F. R., American Statistical Association and American Society for Quality,
1962
"A Simplex Method for Function Minimization", Nelder, J. A. and
Mead, R. The Computer Journal, January, 1965, 308--313
"A New Method of Constrained Optimization and a Comparison With
Other Methods", M. J. Box, The Computer Journal 1965 8(1):42-52, 1965 by
British Computer Society
"Iterative Methods for Optimization", C.T. Kelley, 1999, Chapter 6.,
section 6.2
"Compact Numerical Methods For Computers - Linear Algebra and
Function Minimization", J.C. Nash, 1990, Chapter 14. Direct Search
Methods
"Sequential Application of Simplex Designs in Optimisation and
Evolutionary Operation", W. Spendley, G. R. Hext, F. R. Himsworth,
Technometrics, Vol. 4, No. 4 (Nov., 1962), pp. 441-461, Section 3.1
"A New Method of Constrained Optimization and a Comparison With
Other Methods", M. J. Box, The Computer Journal 1965 8(1):42-52, 1965 by
British Computer Society
“Detection and Remediation of Stagnation in the Nelder--Mead
Algorithm Using a Sufficient Decrease Condition”, SIAM J. on
Optimization, Kelley,, C. T., 1999
" Multi-Directional Search: A Direct Search Algorithm for Parallel
Machines", by E. Boyd, Kenneth W. Kennedy, Richard A. Tapia, Virginia
Joanne Torczon,, Virginia Joanne Torczon, 1989, Phd Thesis, Rice
University
"Grid Restrained Nelder-Mead Algorithm", Árpád Bũrmen, Janez
Puhan, Tadej Tuma, Computational Optimization and Applications, Volume 34
, Issue 3 (July 2006), Pages: 359 - 375
"A convergent variant of the Nelder-Mead algorithm", C. J. Price, I.
D. Coope, D. Byatt, Journal of Optimization Theory and Applications,
Volume 113 , Issue 1 (April 2002), Pages: 5 - 19,
"Global Optimization Of Lennard-Jones Atomic Clusters",
Ellen Fan, Thesis, February 26, 2002, McMaster University