[Index for uc_optim]
[Return to Master Index]
nelder_mead_min
(uc_optim/nelder_mead_min.m)
Function Synopsis
[x,v,nev] = nelder_mead_min (f, args, ctl)
Help text
[x0,v,nev] = nelder_mead_min (f,args,ctl) - Nelder-Mead minimization of f
ARGUMENTS
f : string : Name of function. Must return a real value
args : list : Arguments passed to f.
or matrix : f's only argument
ctl : vector : (Optional) Control variables, described below
RETURNED VALUES
x0 : matrix : Local minimum of f
v : real : Value of f in x0
nev : number : Number of function evaluations
CONTROL VARIABLES
ctl(1) : 1 or 2 : Select stopping criterion amongst :
ctl(1)==1 : Stopping criterion : Stop search when values at simplex
vertices are all alike, as tested by
ctl(2) > (max_i (f_i) - min_i (f_i)) /max(max(|f_i|),1)
where f_i are the values of f at the vertices.
ctl(1)==2 : Stop search when biggest diameter of simplex, using
infinity-norm, is small, as tested by :
ctl(2) > Diam
ctl(1)==3 : Stop search when volume of simplex is small, as tested by
ctl(2) > Vol
Default=1
ctl(2) : Threshold used in stopping tests. Default=10*eps
ctl(3) : Position of the minimized argument in args Default=1
ctl(4) : Maximum number of function evaluations Default=inf
ctl(5) : Size of initial simplex, which is : Default=1
{ x + e_i | i in 0..N }
Where x == nth (args, narg) is the initial value
e_0 == zeros (size (x)),
e_i(j) == 0 if j != i and e_i(i) == ctl(5)
e_i has same size as x
Set ctl(5) to the distance you expect between the starting
point and the minimum.
ctl(6) : When a minimum is found the algorithm restarts next to it
until the minimum does not improve anymore. ctl(6) is the
maximum number of restarts. Set ctl(6) to zero if you know
the function is well-behaved or if you don't mind not
getting a true minimum. Default=inf
ctl may have length smaller than 6. Default values will be used if ctl is
not passed or if nan values are given.
Listing of function file uc_optim/nelder_mead_min.m
## [x0,v,nev] = nelder_mead_min (f,args,ctl) - Nelder-Mead minimization of f
##
## ARGUMENTS
## f : string : Name of function. Must return a real value
## args : list : Arguments passed to f.
## or matrix : f's only argument
## ctl : vector : (Optional) Control variables, described below
##
## RETURNED VALUES
## x0 : matrix : Local minimum of f
## v : real : Value of f in x0
## nev : number : Number of function evaluations
##
## CONTROL VARIABLES
## ctl(1) : 1 or 2 : Select stopping criterion amongst :
## ctl(1)==1 : Stopping criterion : Stop search when values at simplex
## vertices are all alike, as tested by
##
## ctl(2) > (max_i (f_i) - min_i (f_i)) /max(max(|f_i|),1)
##
## where f_i are the values of f at the vertices.
##
## ctl(1)==2 : Stop search when biggest diameter of simplex, using
## infinity-norm, is small, as tested by :
##
## ctl(2) > Diam
##
## ctl(1)==3 : Stop search when volume of simplex is small, as tested by
##
## ctl(2) > Vol
## Default=1
##
## ctl(2) : Threshold used in stopping tests. Default=10*eps
## ctl(3) : Position of the minimized argument in args Default=1
## ctl(4) : Maximum number of function evaluations Default=inf
## ctl(5) : Size of initial simplex, which is : Default=1
##
## { x + e_i | i in 0..N }
##
## Where x == nth (args, narg) is the initial value
## e_0 == zeros (size (x)),
## e_i(j) == 0 if j != i and e_i(i) == ctl(5)
## e_i has same size as x
##
## Set ctl(5) to the distance you expect between the starting
## point and the minimum.
##
## ctl(6) : When a minimum is found the algorithm restarts next to it
## until the minimum does not improve anymore. ctl(6) is the
## maximum number of restarts. Set ctl(6) to zero if you know
## the function is well-behaved or if you don't mind not
## getting a true minimum. Default=inf
##
## ctl may have length smaller than 6. Default values will be used if ctl is
## not passed or if nan values are given.
##
function [x,v,nev] = nelder_mead_min (f, args, ctl)
## Author : Etienne Grossmann <etienne@isr.ist.utl.pt>
## This software is distributed under the terms of the GPL
verbose = 0;
# Default control variables
crit = 1; # Stopping criterion ctl(1)
tol = 10*eps; # Stopping test's threshold ctl(2)
narg = 1; # Position of minimized arg ctl(3)
maxev = inf; # Max num of func evaluations ctl(4)
isz = 1; # Initial size ctl(5)
rst = inf; # Max # of restarts
if nargin >= 3, # Read arguments
if !isnan (ctl(1)), crit = ctl(1); end
if length (ctl)>=2 && !isnan (ctl(2)), tol = ctl(2); end
if length (ctl)>=3 && !isnan (ctl(3)), narg = ctl(3); end
if length (ctl)>=4 && !isnan (ctl(4)), maxev = ctl(4); end
if length (ctl)>=5 && !isnan (ctl(5)), isz = ctl(5); end
if length (ctl)>=6 && !isnan (ctl(6)), rst = ctl(6); end
end
if is_list (args), # List of arguments
x = nth (args, narg);
else # Single argument
x = args;
args = list (args);
end
if narg > length (args), # Check
error ("nelder_mead_min : narg==%i, length (args)==%i\n",
narg, length (args));
end
[R,C] = size (x);
N = R*C; # Size of argument
x = x(:);
# Initial simplex
u = isz * eye (N+1,N) + ones(N+1,1)*x';
for i = 1:N+1,
y(i) = leval (f, splice (args, narg, 1, list (reshape (u(i,:),R,C))));
end ;
nev = N+1;
ymin0 = min(y) ;
## y
nextprint = 0 ;
v = nan;
while nev <= maxev,
## ymin, ymax, ymx2 : lowest, highest and 2nd highest function values
## imin, imax, imx2 : indices of vertices with these values
[ymin,imin] = min(y) ; [ymax,imax] = max(y) ;
y(imax) = ymin ;
[ymx2,imx2] = max(y) ;
y(imax) = ymax ;
if ymin > ymin0 ,
"nelder-mead : Whoa 'downsimplex' Should be renamed 'upsimplex'"
keyboard
end
# Compute stopping criterion
if crit == 1,
rt = (max(y)-min(y)) / max(1,max(abs(y)));
elseif crit == 2,
rt = max (max (u) - min (u));
else
rt = abs(det(u(1:N,:)-ones(N,1)*u(N+1,:)))^(1/N) ;
## rt = (max(y)-min(y)) / max (max (u) - min(u));
end
# Eventually print some info
if nev > nextprint
if 0 && verbose,
printf("nev=%-5d imin=%-3d ymin=%-8.3g rt=%-8.3g\n",\
nev,imin,ymin,rt) ;
end
nextprint = nextprint + 100 ;
end
if rt < tol , # Termination test
if (rst > 0) && (isnan (v) || v > ymin)
rst--;
if verbose
if isnan (v),
printf ("Restarting next to minimum %10.3e\n",ymin);
else
printf ("Restarting next to minimum %10.3e\n",ymin-v);
end
end
# Keep best minimum
x = reshape (u(imin,:), R, C) ;
v = ymin ;
jumplen = 10 * max (max (u) - min (u));
u += jumplen * randn (size (u));
for i = 1:N+1, y(i) = \
leval (f, splice (args, narg, 1, list (reshape (u(i,:),R,C))));
end
nev += N+1;
[ymin,imin] = min(y) ; [ymax,imax] = max(y) ;
y(imax) = ymin ;
[ymx2,imx2] = max(y) ;
y(imax) = ymax ;
else
if isnan (v),
x = reshape (u(imin,:), R, C) ;
v = ymin ;
end
if verbose,
printf("nev=%-5d imin=%-3d ymin=%-8.3g rt=%-8.3g. Done\n",\
nev,imin,ymin,rt) ;
end
return
end
end
## [ y' u ]
tra = 0 ; # 'trace' debug var contains flags
if verbose, str = sprintf ("%10.3e : %10.3e --",rt,ymin); end
# Look for a new point
xsum = sum(u) ; # Consider reflection of worst vertice
# around centroid.
## f1 = (1-(-1))/N = 2/N;
## f2 = f1 - (-1) = 2/N + 1 = (N+2)/N
xnew = (2*xsum - (N+2)*u(imax,:)) / N;
## xnew = (2*xsum - N*u(imax,:)) / N;
ynew = leval (f, splice (args, narg, 1, list ( reshape (xnew, R,C))));
nev++;
if ynew <= ymin , # Reflection is good
tra += 1 ;
if verbose,
str = [str,sprintf (" %3i : %10.3e good refl >>",nev,ynew-ymin)];
end
y(imax) = ynew; u(imax,:) = xnew ;
## ymin = ynew;
## imin = imax;
xsum = sum(u) ;
## f1 = (1-2)/N = -1/N
## f2 = f1 - 2 = -1/N - 2 = -(2*N+1)/N
xnew = ( -xsum + (2*N+1)*u(imax,:) ) / N;
ynew = leval (f, splice (args, narg, 1, list ( reshape (xnew, R,C))));
nev++;
if ynew <= ymin , # expansion improves
tra += 2 ;
## 'expanded reflection'
y(imax) = ynew ; u(imax,:) = xnew ;
xsum = sum(u) ;
if verbose,
str = [str,sprintf (" %3i : %10.3e expd refl",nev,ynew-ymin)];
end
else
tra += 4 ;
## 'plain reflection'
## Updating of y and u has already been done
if verbose,
str = [str,sprintf (" %3i : %10.3e plain ref",nev,ynew-ymin)];
end
end
# Reflexion is really bad
elseif ynew >= ymax ,
tra += 8 ;
if verbose,
str = [str,sprintf (" %3i : %10.3e intermedt >>",nev,ynew-ymin)];
end
## look for intermediate point
# Bring worst point closer to centroid
## f1 = (1-0.5)/N = 0.5/N
## f2 = f1 - 0.5 = 0.5*(1 - N)/N
xnew = 0.5*(xsum + (N-1)*u(imax,:)) / N;
ynew = leval (f, splice (args, narg, 1, list (reshape (xnew, R,C))));
nev++;
if ynew >= ymax , # New point is even worse. Contract whole
# simplex
nev += N + 1 ;
## u0 = u;
u = (u + ones(N+1,1)*u(imin,:)) / 2;
## keyboard
for i = [1:imin-1,imin+1:N+1],
y(i) = \
leval (f, splice (args, narg, 1, list (reshape (u(i,:),R,C))));
end
## 'contraction'
tra += 16 ;
if verbose,
str = [str,sprintf (" %3i contractn",nev)];
end
else # Replace highest point
y(imax) = ynew ; u(imax,:) = xnew ;
xsum = sum(u) ;
## 'intermediate'
tra += 32 ;
if verbose,
str = [str,sprintf (" %3i : %10.3e intermedt",nev,ynew-ymin)];
end
end
else # Reflexion is neither good nor bad
y(imax) = ynew ; u(imax,:) = xnew ;
xsum = sum(u) ;
## 'plain reflection (2)'
tra += 64 ;
if verbose,
str = [str,sprintf (" %3i : %10.3e keep refl",nev,ynew-ymin)];
end
end
if verbose, printf ("%s\n",str); end
## [tra, nev, ymin, rt]
## [ymin, rt]
## if !isempty (input("","s")), keyboard; end
end
printf ("nelder_mead : Too many iterations. Returning\n");
if isnan (v) || v > ymin,
x = reshape (u(imin,:), R, C) ;
v = ymin ;
end
Produced by oct2html on Sun Feb 11 12:59:56 2001
Cross-Directory links are: OFF