octave_packages/optim-1.2.0/nelder_mead_min.m

   1 ## Copyright (C) 2002-2008 Etienne Grossmann <etienne@egdn.net>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## [x0,v,nev] = nelder_mead_min (f,args,ctl) - Nelder-Mead minimization
  17 ##
  18 ## Minimize 'f' using the Nelder-Mead algorithm. This function is inspired
  19 ## from the that found in the book "Numerical Recipes".
  20 ##
  21 ## ARGUMENTS
  22 ## ---------
  23 ## f     : string : Name of function. Must return a real value
  24 ## args  : list   : Arguments passed to f.
  25 ##      or matrix : f's only argument
  26 ## ctl   : vector : (Optional) Control variables, described below
  27 ##      or struct
  28 ##
  29 ## RETURNED VALUES
  30 ## ---------------
  31 ## x0  : matrix   : Local minimum of f
  32 ## v   : real     : Value of f in x0
  33 ## nev : number   : Number of function evaluations
  34 ##
  35 ## CONTROL VARIABLE : (optional) may be named arguments (i.e. "name",value
  36 ## ------------------ pairs), a struct, or a vector of length <= 6, where
  37 ##                    NaN's are ignored. Default values are written <value>.
  38 ##  OPT.   VECTOR
  39 ##  NAME    POS
  40 ## ftol,f  N/A    : Stopping criterion : stop search when values at simplex
  41 ##                  vertices are all alike, as tested by
  42 ##
  43 ##                   f > (max_i (f_i) - min_i (f_i)) /max(max(|f_i|),1)
  44 ##
  45 ##                  where f_i are the values of f at the vertices.  <10*eps>
  46 ##
  47 ## rtol,r  N/A    : Stop search when biggest radius of simplex, using
  48 ##                  infinity-norm, is small, as tested by :
  49 ##
  50 ##              ctl(2) > Radius                                     <10*eps>
  51 ##
  52 ## vtol,v  N/A    : Stop search when volume of simplex is small, tested by
  53 ##
  54 ##              ctl(2) > Vol
  55 ##
  56 ## crit,c ctl(1)  : Set one stopping criterion, 'ftol' (c=1), 'rtol' (c=2)
  57 ##                  or 'vtol' (c=3) to the value of the 'tol' option.    <1>
  58 ##
  59 ## tol, t ctl(2)  : Threshold in termination test chosen by 'crit'  <10*eps>
  60 ##
  61 ## narg  ctl(3)  : Position of the minimized argument in args            <1>
  62 ## maxev ctl(4)  : Maximum number of function evaluations. This number <inf>
  63 ##                 may be slightly exceeded.
  64 ## isz   ctl(5)  : Size of initial simplex, which is :                   <1>
  65 ##
  66 ##                { x + e_i | i in 0..N }
  67 ##
  68 ##                Where x == args{narg} is the initial value
  69 ##                 e_0    == zeros (size (x)),
  70 ##                 e_i(j) == 0 if j != i and e_i(i) == ctl(5)
  71 ##                 e_i    has same size as x
  72 ##
  73 ##                Set ctl(5) to the distance you expect between the starting
  74 ##                point and the minimum.
  75 ##
  76 ## rst   ctl(6)   : When a minimum is found the algorithm restarts next to
  77 ##                  it until the minimum does not improve anymore. ctl(6) is
  78 ##                  the maximum number of restarts. Set ctl(6) to zero if
  79 ##                  you know the function is well-behaved or if you don't
  80 ##                  mind not getting a true minimum.                     <0>
  81 ##
  82 ## verbose, v     Be more or less verbose (quiet=0)                      <0>
  83
  84 function [x,v,nev] = nelder_mead_min (f, args, varargin)
  85
  86 verbose = 0;
  87
  88                                 # Default control variables
  89 ftol = rtol = 10*eps;           # Stop either by likeness of values or
  90 vtol = nan;                     # radius, but don't care about volume.
  91 crit = 0;                       # Stopping criterion            ctl(1)
  92 tol = 10*eps;                   # Stopping test's threshold     ctl(2)
  93 narg = 1;                       # Position of minimized arg     ctl(3)
  94 maxev = inf;                    # Max num of func evaluations   ctl(4)
  95 isz = 1;                        # Initial size                  ctl(5)
  96 rst = 0;                        # Max # of restarts
  97
  98
  99 if nargin >= 3,                 # Read control arguments
 100   va_arg_cnt = 1;
 101   if nargin > 3,
 102           ctl = struct (varargin{:});
 103   else
 104           ctl = varargin{va_arg_cnt++};
 105   end
 106   if isnumeric (ctl)
 107     if length (ctl)>=1 && !isnan (ctl(1)), crit = ctl(1); end
 108     if length (ctl)>=2 && !isnan (ctl(2)), tol = ctl(2); end
 109     if length (ctl)>=3 && !isnan (ctl(3)), narg = ctl(3); end
 110     if length (ctl)>=4 && !isnan (ctl(4)), maxev = ctl(4); end
 111     if length (ctl)>=5 && !isnan (ctl(5)), isz = ctl(5); end
 112     if length (ctl)>=6 && !isnan (ctl(6)), rst = ctl(6); end
 113   else
 114     if isfield (ctl, "crit") && ! isnan (ctl.crit ), crit  = ctl.crit ; end
 115     if isfield (ctl,  "tol") && ! isnan (ctl.tol  ), tol   = ctl.tol  ; end
 116     if isfield (ctl, "ftol") && ! isnan (ctl.ftol ), ftol  = ctl.ftol ; end
 117     if isfield (ctl, "rtol") && ! isnan (ctl.rtol ), rtol  = ctl.rtol ; end
 118     if isfield (ctl, "vtol") && ! isnan (ctl.vtol ), vtol  = ctl.vtol ; end
 119     if isfield (ctl, "narg") && ! isnan (ctl.narg ), narg  = ctl.narg ; end
 120     if isfield (ctl,"maxev") && ! isnan (ctl.maxev), maxev = ctl.maxev; end
 121     if isfield (ctl,  "isz") && ! isnan (ctl.isz  ), isz   = ctl.isz  ; end
 122     if isfield (ctl,  "rst") && ! isnan (ctl.rst  ), rst   = ctl.rst  ; end
 123     if isfield(ctl,"verbose")&& !isnan(ctl.verbose),verbose=ctl.verbose;end
 124   end
 125 end
 126
 127
 128 if     crit == 1, ftol = tol;
 129 elseif crit == 2, rtol = tol;
 130 elseif crit == 3, vtol = tol;
 131 elseif crit, error ("crit is %i. Should be 1,2 or 3.\n");
 132 end
 133
 134 if iscell (args)
 135   x = args{1};
 136 else                            # Single argument
 137   x = args;
 138   args = {args};
 139 endif
 140
 141 if narg > length (args)         # Check
 142   error ("nelder_mead_min : narg==%i, length (args)==%i\n",
 143          narg, length (args));
 144 end
 145
 146 [R,C] = size (x);
 147 N = R*C;                        # Size of argument
 148 x = x(:);
 149                                 # Initial simplex
 150 u = isz * eye (N+1,N) + ones(N+1,1)*x';
 151
 152 y = zeros (N+1,1);
 153 for i = 1:N+1,
 154   y(i) = feval (f, args{1:narg-1},reshape(u(i,:),R,C),args{narg+1:end});
 155 end ;
 156 nev = N+1;
 157
 158 [ymin,imin] = min(y);
 159 ymin0 = ymin;
 160 ## y
 161 nextprint = 0 ;
 162 v = nan;
 163 while nev <= maxev,
 164
 165   ## ymin, ymax, ymx2 : lowest, highest and 2nd highest function values
 166   ## imin, imax, imx2 : indices of vertices with these values
 167   [ymin,imin] = min(y);  [ymax,imax] = max(y) ;
 168   y(imax) = ymin ;
 169   [ymx2,imx2] = max(y) ;
 170   y(imax) = ymax ;
 171
 172   ## ymin may be > ymin0 after restarting
 173   ## if ymin > ymin0 ,
 174   ## "nelder-mead : Whoa 'downsimplex' Should be renamed 'upsimplex'"
 175   ## keyboard
 176   ## end
 177
 178                                 # Compute stopping criterion
 179   done = 0;
 180   if ! isnan (ftol),
 181      done |= ((max(y)-min(y)) / max(1,max(abs(y))) < ftol);
 182   end
 183   if ! isnan (rtol),
 184      done |= (2*max (max (u) - min (u)) < rtol);
 185   end
 186   if ! isnan (vtol)
 187     done |= (abs (det (u(1:N,:)-ones(N,1)*u(N+1,:)))/factorial(N) < vtol);
 188   end
 189   ## [ 2*max (max (u) - min (u)), abs (det (u(1:N,:)-ones(N,1)*u(N+1,:)))/factorial(N);\
 190   ##  rtol, vtol]
 191
 192                                 # Eventually print some info
 193   if verbose && nev > nextprint && ! done
 194
 195     printf("nev=%-5d   imin=%-3d   ymin=%-8.3g  done=%i\n",\
 196            nev,imin,ymin,done) ;
 197
 198     nextprint = nextprint + 100 ;
 199   end
 200
 201   if done                       # Termination test
 202     if (rst > 0) && (isnan (v) || v > ymin)
 203       rst--;
 204       if verbose
 205         if isnan (v),
 206           printf ("Restarting next to minimum %10.3e\n",ymin);
 207         else
 208           printf ("Restarting next to minimum %10.3e\n",ymin-v);
 209         end
 210       end
 211                                 # Keep best minimum
 212       x = reshape (u(imin,:), R, C) ;
 213       v = ymin ;
 214
 215       jumplen = 10 * max (max (u) - min (u));
 216
 217       u += jumplen * randn (size (u));
 218       for i = 1:N+1, y(i) = \
 219             feval (f, args{1:narg-1},reshape(u(i,:),R,C),args{narg+1:length(args)});
 220       end
 221       nev += N+1;
 222       [ymin,imin] = min(y);  [ymax,imax] = max(y);
 223       y(imax) = ymin;
 224       [ymx2,imx2] = max(y);
 225       y(imax) = ymax ;
 226     else
 227       if isnan (v),
 228         x = reshape (u(imin,:), R, C) ;
 229         v = ymin ;
 230       end
 231       if verbose,
 232         printf("nev=%-5d   imin=%-3d   ymin=%-8.3g  done=%i. Done\n",\
 233                nev,imin,ymin,done) ;
 234       end
 235       return
 236     end
 237
 238   end
 239   ##   [ y' u ]
 240
 241   tra = 0 ;                     # 'trace' debug var contains flags
 242   if verbose > 1, str = sprintf (" %i : %10.3e --",done,ymin); end
 243
 244                                 # Look for a new point
 245   xsum = sum(u) ;               # Consider reflection of worst vertice
 246                                 # around centroid.
 247   ## f1 = (1-(-1))/N = 2/N;
 248   ## f2 = f1 - (-1)  = 2/N + 1 = (N+2)/N
 249   xnew = (2*xsum - (N+2)*u(imax,:)) / N;
 250   ## xnew = (2*xsum - N*u(imax,:)) / N;
 251   ynew = feval (f, args{1:narg-1},reshape(xnew,R,C),args{narg+1:length(args)});
 252   nev++;
 253
 254   if ynew <= ymin ,             # Reflection is good
 255
 256     tra += 1 ;
 257     if verbose > 1
 258       str = [str,sprintf(" %3i : %10.3e good refl >>",nev,ynew-ymin)];
 259     end
 260     y(imax) = ynew; u(imax,:) = xnew ;
 261     ## ymin = ynew;
 262     ## imin = imax;
 263     xsum = sum(u) ;
 264
 265     ## f1 = (1-2)/N = -1/N
 266     ## f2 = f1 - 2  = -1/N - 2 = -(2*N+1)/N
 267     xnew = ( -xsum + (2*N+1)*u(imax,:) ) / N;
 268     ynew = feval (f, args{1:narg-1},reshape(xnew,R,C),args{narg+1:length(args)});
 269     nev++;
 270
 271     if ynew <= ymin ,           # expansion improves
 272       tra += 2 ;
 273       ##      'expanded reflection'
 274       y(imax) = ynew ; u(imax,:) = xnew ;
 275       xsum = sum(u) ;
 276       if verbose > 1
 277         str = [str,sprintf(" %3i : %10.3e expd refl",nev,ynew-ymin)];
 278       end
 279     else
 280       tra += 4 ;
 281       ##      'plain reflection'
 282       ## Updating of y and u has already been done
 283       if verbose > 1
 284         str = [str,sprintf(" %3i : %10.3e plain ref",nev,ynew-ymin)];
 285       end
 286     end
 287                                 # Reflexion is really bad
 288   elseif ynew >= ymax ,
 289
 290     tra += 8 ;
 291     if verbose > 1
 292       str = [str,sprintf(" %3i : %10.3e intermedt >>",nev,ynew-ymin)];
 293     end
 294     ## look for intermediate point
 295                                 # Bring worst point closer to centroid
 296     ## f1 = (1-0.5)/N = 0.5/N
 297     ## f2 = f1 - 0.5  = 0.5*(1 - N)/N
 298     xnew = 0.5*(xsum + (N-1)*u(imax,:)) / N;
 299     ynew = feval (f, args{1:narg-1},reshape(xnew,R,C),args{narg+1:length(args)});
 300     nev++;
 301
 302     if ynew >= ymax ,           # New point is even worse. Contract whole
 303                                 # simplex
 304
 305       nev += N + 1 ;
 306       ## u0 = u;
 307       u = (u + ones(N+1,1)*u(imin,:)) / 2;
 308       ## keyboard
 309
 310       ## Code that doesn't care about value of empty_list_elements_ok
 311       if     imin == 1  , ii = 2:N+1;
 312       elseif imin == N+1, ii = 1:N;
 313       else                ii = [1:imin-1,imin+1:N+1]; end
 314       for i = ii
 315         y(i) = \
 316             ynew = feval (f, args{1:narg-1},reshape(u(i,:),R,C),args{narg+1:length(args)});
 317       end
 318       ##      'contraction'
 319       tra += 16 ;
 320       if verbose > 1
 321         str = [str,sprintf(" %3i contractn",nev)];
 322       end
 323     else                                # Replace highest point
 324       y(imax) = ynew ; u(imax,:) = xnew ;
 325       xsum = sum(u) ;
 326       ##      'intermediate'
 327       tra += 32 ;
 328       if verbose > 1
 329         str = [str,sprintf(" %3i : %10.3e intermedt",nev,ynew-ymin)];
 330       end
 331     end
 332
 333   else                          # Reflexion is neither good nor bad
 334     y(imax) = ynew ; u(imax,:) = xnew ;
 335     xsum = sum(u) ;
 336     ##      'plain reflection (2)'
 337     tra += 64 ;
 338     if verbose > 1
 339       str = [str,sprintf(" %3i : %10.3e keep refl",nev,ynew-ymin)];
 340     end
 341   end
 342   if verbose > 1, printf ("%s\n",str); end
 343 end
 344
 345 if verbose >= 0
 346   printf ("nelder_mead : Too many iterations. Returning\n");
 347 end
 348
 349 if isnan (v) || v > ymin,
 350   x = reshape (u(imin,:), R, C) ;
 351   v = ymin ;
 352 end