octave_packages/optim-1.2.0/brent_line_min.m

   1 ## Copyright (C) 2009 Levente Torok <TorokLev@gmail.com>
   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 ## -*- texinfo -*-
  17 ## @deftypefn {Function File} {[@var{s},@var{v},@var{n}]} brent_line_min ( @var{f},@var{df},@var{args},@var{ctl} )
  18 ## Line minimization of f along df
  19 ##
  20 ## Finds minimum of f on line @math{ x0 + dx*w | a < w < b } by
  21 ## bracketing. a and b are passed through argument ctl.
  22 ##
  23 ## @subheading Arguments
  24 ## @itemize @bullet
  25 ## @item @var{f}     : string : Name of function. Must return a real value
  26 ## @item @var{args}  : cell   : Arguments passed to f or RxC    : f's only argument. x0 must be at @var{args}@{ @var{ctl}(2) @}
  27 ## @item @var{ctl}   : 5      : (optional) Control variables, described below.
  28 ## @end itemize
  29 ##
  30 ## @subheading Returned values
  31 ## @itemize @bullet
  32 ## @item @var{s}   : 1        : Minimum is at x0 + s*dx
  33 ## @item @var{v}   : 1        : Value of f at x0 + s*dx
  34 ## @item @var{nev} : 1        : Number of function evaluations
  35 ## @end itemize
  36 ##
  37 ## @subheading Control Variables
  38 ## @itemize @bullet
  39 ## @item @var{ctl}(1)       : Upper bound for error on s              Default=sqrt(eps)
  40 ## @item @var{ctl}(2)       : Position of minimized argument in args  Default= 1
  41 ## @item @var{ctl}(3)       : Maximum number of function evaluations  Default= inf
  42 ## @item @var{ctl}(4)       : a                                       Default=-inf
  43 ## @item @var{ctl}(5)       : b                                       Default= inf
  44 ## @end itemize
  45 ##
  46 ## Default values will be used if ctl is not passed or if nan values are
  47 ## given.
  48 ## @end deftypefn
  49
  50 function [s,gs,nev] = brent_line_min( f,dx,args,ctl )
  51
  52 verbose = 0;
  53
  54 seps = sqrt (eps);
  55
  56                                 # Default control variables
  57 tol = 10*eps; # sqrt (eps);
  58 narg = 1;
  59 maxev = inf;
  60 a = -inf;
  61 b = inf;
  62
  63 if nargin >= 4,                 # Read arguments
  64   if !isnan (ctl (1)), tol = ctl(1); end
  65   if length (ctl)>=2 && !isnan (ctl(2)), narg = ctl(2); end
  66   if length (ctl)>=3 && !isnan (ctl(3)), maxev = ctl(3); end
  67   if length (ctl)>=4 && !isnan (ctl(4)), a = ctl(4); end
  68   if length (ctl)>=5 && !isnan (ctl(5)), b = ctl(5); end
  69
  70 end                             # Otherwise, use defaults, def'd above
  71
  72 if a>b, tmp=a; a=b; b=tmp; end
  73
  74 if narg > length (args),
  75   printf ("brent_line_min : narg==%i > length (args)==%i",\
  76     narg, length (args));
  77   keyboard
  78 end
  79
  80
  81 if ! iscell (args),
  82         args = {args};
  83 endif
  84
  85 x = args{ narg };
  86
  87 [R,C] = size (x);
  88 N = R*C;                        # Size of argument
  89
  90 gs0 = gs = feval (f, args);
  91 nev = 1;
  92                                 # Initial value
  93 s = 0;
  94
  95 if all (dx==0), return; end     # Trivial case
  96
  97                                 # If any of the bounds is infinite, find
  98                                 # finite bounds that bracket minimum
  99 if !isfinite (a) || !isfinite (b),
 100   if !isfinite (a) && !isfinite (b),
 101     [a,b, ga,gb, n] = __bracket_min (f, dx, narg, args);
 102   elseif !isfinite (a),
 103     [a,b, ga,gb, n] = __semi_bracket (f, -dx, -b, narg, args);
 104     tmp = a;  a = -b; b = -tmp;
 105     tmp = ga; ga = gb; gb = tmp;
 106   else
 107     [a,b, ga,gb, n] = __semi_bracket (f,  dx, a, narg, args);
 108   end
 109   nev += n;
 110 else
 111   args{narg} = x+a*dx;  ga = feval( f, args );
 112   args{narg} = x+b*dx;  gb = feval( f, args );
 113   nev += 2;
 114 end                             # End of finding bracket for minimum
 115
 116 if a > b,                       # Check assumptions
 117   printf ("brent_line_min : a > b\n");
 118   keyboard
 119 end
 120
 121 s = 0.5*(a+b);
 122 args{narg} = x+ s*dx; gs = feval( f, args );
 123 nev++;
 124
 125 if verbose,
 126   printf ("[a,s,b]=[%.3e,%.3e,%.3e], [ga,gs,gb]=[%.3e,%.3e,%.3e]\n",\
 127           a,s,b,ga,gs,gb);
 128 end
 129
 130 maxerr = 2*tol;
 131
 132 while ( b-a > maxerr ) && nev < maxev,
 133
 134   if gs > ga && gs > gb,        # Check assumptions
 135     printf ("brent_line_min : gs > ga && gs > gb\n");
 136     keyboard
 137   end
 138
 139   if ga == gb && gb == gs, break end
 140
 141                                 # Don't trust poly_2_ex for glued points
 142                                 # (see test_poly_2_ex).
 143
 144   ## if min (b-s, s-a) > 10*seps,
 145
 146                                 # If s is not glued to a or b and does not
 147                                 # look linear
 148   ## mydet = sum (l([2 3 1]).*f([3 1 2])-l([3 1 2]).*f([2 3 1]))
 149   mydet = sum ([s b a].*[gb ga gs] - [b a s].*[gs gb ga]);
 150   if min (b-s, s-a) > 10*seps && abs (mydet) > 10*seps && \
 151         (t = poly_2_ex ([a,s,b], [ga, gs, gb])) < b && t > a,
 152
 153                                 # t has already been set
 154
 155     ## if t>=b || t<=a,
 156     ##  printf ("brent_line_min : t is not in ]a,b[\n");
 157     ##  keyboard
 158     ## end
 159
 160                                 # Otherwise, reduce the biggest of the
 161                                 # intervals, but not too much
 162   elseif s-a > b-s,
 163     t = max (0.5*(a+s), s-100*seps);
 164   else
 165     t = min (0.5*(s+b), s+100*seps);
 166   end
 167
 168   if abs (t-s) < 0.51*maxerr,
 169     #sayif (verbose, "ungluing t and s\n");
 170     t = s + (1 - 2*(s-a > b-s))*0.49*maxerr ;
 171   end
 172
 173   if a > s || s > b,            # Check assumptions
 174     printf ("brent_line_min : a > s || s > b\n");
 175     keyboard
 176   end
 177
 178   xt = args;
 179   args{narg} = x+t*dx;
 180   gt = feval( f, args );
 181   nev++;
 182
 183   if verbose,
 184     printf ("t = %.3e, gt = %.3e\n",t,gt);
 185   end
 186
 187   if t<s,                       # New point is in ]a,s[
 188
 189     if gt > ga + seps,          # Check assumptions
 190       printf ("brent_line_min : gt > ga\n");
 191       keyboard
 192     end
 193
 194     if gt < gs,
 195       b = s; gb = gs;
 196       s = t; gs = gt;
 197     else
 198       a = t; ga = gt;
 199     end
 200   else                          # New point is in ]s,b[
 201     if gt > gb + seps,          # Check assumptions
 202       printf ("brent_line_min : gt > gb\n");
 203       keyboard
 204     end
 205
 206     if gt < gs,
 207       a = s; ga = gs;
 208       s = t; gs = gt;
 209     else
 210       b = t; gb = gt;
 211     end
 212   end
 213
 214   if verbose,
 215     printf ("[a,s,b]=[%.3e,%.3e,%.3e], [ga,gs,gb]=[%.3e,%.3e,%.3e]\n",\
 216             a,s,b,ga,gs,gb);
 217   end
 218   ## keyboard
 219   ## [b-a, maxerr]
 220 end
 221
 222 s2 = 0.5*(a+b);
 223 args{narg} = x + s2*dx; gs2 = feval (f, args );
 224 nev++;
 225
 226 if gs2 < gs,
 227   s = s2; gs = gs2;
 228 end
 229
 230 if gs > gs0,
 231   printf ("brent_line_min : goes uphill by %8.3e\n",gs-gs0);
 232   keyboard
 233 end