octave_packages/optim-1.2.0/line_min.m

   1 ## Copyright (C) 2000 Ben Sapp <bsapp@lanl.gov>
   2 ## Copyright (C) 2002 Etienne Grossmann <etienne@egdn.net>
   3 ## Copyright (C) 2011 Nir Krakauer nkrakauer@ccny.cuny.edu
   4 ##
   5 ## This program is free software; you can redistribute it and/or modify it under
   6 ## the terms of the GNU General Public License as published by the Free Software
   7 ## Foundation; either version 3 of the License, or (at your option) any later
   8 ## version.
   9 ##
  10 ## This program is distributed in the hope that it will be useful, but WITHOUT
  11 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  12 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  13 ## details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License along with
  16 ## this program; if not, see <http://www.gnu.org/licenses/>.
  17
  18 ## [a,fx,nev] = line_min (f, dx, args, narg, h, nev_max) - Minimize f() along dx
  19 ##
  20 ## INPUT ----------
  21 ## f    : string  : Name of minimized function
  22 ## dx   : matrix  : Direction along which f() is minimized
  23 ## args : cell    : Arguments of f
  24 ## narg : integer : Position of minimized variable in args.  Default=1
  25 ## h    : scalar  : Step size to use for centered finite difference
  26 ## approximation of first and second derivatives. Default=1E-3.
  27 ## nev_max : integer : Maximum number of function evaluations.  Default=30
  28 ##
  29 ## OUTPUT ---------
  30 ## a    : scalar  : Value for which f(x+a*dx) is a minimum (*)
  31 ## fx   : scalar  : Value of f(x+a*dx) at minimum (*)
  32 ## nev  : integer : Number of function evaluations
  33 ##
  34 ## (*) The notation f(x+a*dx) assumes that args == {x}.
  35 ##
  36 ## Reference: David G Luenberger's Linear and Nonlinear Programming
  37
  38 function [a,fx,nev] = line_min (f, dx, args, narg, h, nev_max)
  39   velocity = 1;
  40   acceleration = 1;
  41
  42   if (nargin < 4) narg = 1; endif
  43   if (nargin < 5) h = 0.001; endif
  44   if (nargin < 6) nev_max = 30; endif
  45
  46   nev = 0;
  47   x = args{narg};
  48   a = 0;
  49
  50   min_velocity_change = 0.000001;
  51
  52   while (abs (velocity) > min_velocity_change && nev < nev_max)
  53     fx = feval (f,args{1:narg-1}, x+a*dx, args{narg+1:end});
  54     fxph = feval (f,args{1:narg-1}, x+(a+h)*dx, args{narg+1:end});
  55     fxmh = feval (f,args{1:narg-1}, x+(a-h)*dx, args{narg+1:end});
  56     if (nev == 0)
  57         fx0 = fx;
  58     endif
  59
  60     velocity = (fxph - fxmh)/(2*h);
  61     acceleration = (fxph - 2*fx + fxmh)/(h^2);
  62     if abs(acceleration) <= eps, acceleration = 1; end # Don't do div by zero
  63                                 # Use abs(accel) to avoid problems due to
  64                                 # concave function
  65     a = a - velocity/abs(acceleration);
  66     nev += 3;
  67   endwhile
  68
  69   fx = feval (f, args{1:narg-1}, x+a*dx, args{narg+1:end});
  70   nev++;
  71   if fx >= fx0 # if no improvement, return the starting value
  72         a = 0;
  73         fx = fx0;
  74   endif
  75
  76   if (nev >= nev_max)
  77     disp ("line_min: maximum number of function evaluations reached")
  78   endif
  79
  80 endfunction
  81
  82 ## Rem : Although not clear from the code, the returned a always seems to
  83 ## correspond to (nearly) optimal fx.