octave_packages/optim-1.2.0/d2_min.m

   1 ## Copyright (C) 2002 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 ## [x,v,nev,h,args] = d2_min(f,d2f,args,ctl,code) - Newton-like minimization
  17 ##
  18 ## Minimize f(x) using 1st and 2nd derivatives. Any function w/ second
  19 ## derivatives can be minimized, as in Newton. f(x) decreases at each
  20 ## iteration, as in Levenberg-Marquardt. This function is inspired from the
  21 ## Levenberg-Marquardt algorithm found in the book "Numerical Recipes".
  22 ##
  23 ## ARGUMENTS :
  24 ## f    : string : Cost function's name
  25 ##
  26 ## d2f  : string : Name of function returning the cost (1x1), its
  27 ##                 differential (1xN) and its second differential or it's
  28 ##                 pseudo-inverse (NxN) (see ctl(5) below) :
  29 ##
  30 ##                 [v,dv,d2v] = d2f (x).
  31 ##
  32 ## args : list   : f and d2f's arguments. By default, minimize the 1st
  33 ##     or matrix : argument.
  34 ##
  35 ## ctl  : vector : Control arguments (see below)
  36 ##      or struct
  37 ##
  38 ## code : string : code will be evaluated after each outer loop that
  39 ##                 produced some (any) improvement. Variables visible from
  40 ##                 "code" include "x", the best parameter found, "v" the
  41 ##                 best value and "args", the list of all arguments. All can
  42 ##                 be modified. This option can be used to re-parameterize
  43 ##                 the argument space during optimization
  44 ##
  45 ## CONTROL VARIABLE ctl : (optional). May be a struct or a vector of length
  46 ## ---------------------- 5 or less where NaNs are ignored. Default values
  47 ##                        are written <value>.
  48 ## FIELD  VECTOR
  49 ## NAME    POS
  50 ##
  51 ## ftol, f N/A    : Stop search when value doesn't improve, as tested by
  52 ##
  53 ##                   f > Deltaf/max(|f(x)|,1)
  54 ##
  55 ##             where Deltaf is the decrease in f observed in the last
  56 ##             iteration.                                     <10*sqrt(eps)>
  57 ##
  58 ## utol, u N/A    : Stop search when updates are small, as tested by
  59 ##
  60 ##                   u > max { dx(i)/max(|x(i)|,1) | i in 1..N }
  61 ##
  62 ##             where  dx is the change in the x that occured in the last
  63 ##             iteration.                                              <NaN>
  64 ##
  65 ## dtol, d N/A    : Stop search when derivative is small, as tested by
  66 ##
  67 ##                   d > norm (dv)                                     <eps>
  68 ##
  69 ## crit, c ctl(1) : Set one stopping criterion, 'ftol' (c=1), 'utol' (c=2)
  70 ##                  or 'dtol' (c=3) to the value of by the 'tol' option. <1>
  71 ##
  72 ## tol, t  ctl(2) : Threshold in termination test chosen by 'crit'  <10*eps>
  73 ##
  74 ## narg, n ctl(3) : Position of the minimized argument in args           <1>
  75 ## maxev,m ctl(4) : Maximum number of function evaluations             <inf>
  76 ## maxout,m       : Maximum number of outer loops                      <inf>
  77 ## id2f, i ctl(5) : 0 if d2f returns the 2nd derivatives, 1 if           <0>
  78 ##                  it returns its pseudo-inverse.
  79 ##
  80 ## verbose, v N/A : Be more or less verbose (quiet=0)                    <0>
  81
  82 function [xbest,vbest,nev,hbest,args] = d2_min (f,d2f,args,ctl,code)
  83
  84 maxout = inf;
  85 maxinner = 30 ;
  86
  87 tcoeff = 0.5 ;                  # Discount on total weight
  88 ncoeff = 0.5 ;                  # Discount on weight of newton
  89 ocoeff = 1.5 ;                  # Factor for outwards searching
  90
  91 report = 0 ;                    # Never report
  92 verbose = 0 ;                   # Be quiet
  93 prudent = 1 ;                   # Check coherence of d2f and f?
  94
  95 niter = 0 ;
  96
  97 crit = 0;                       # Default control variables
  98 ftol = 10 * sqrt (eps);
  99 dtol = eps;
 100 utol = tol = nan;
 101 narg = 1;
 102 maxev = inf;
 103 id2f = 0;
 104
 105 if nargin >= 4                  # Read arguments
 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)), id2f  = ctl(5); end
 112   elseif isstruct (ctl)
 113     if isfield (ctl, "crit")   , crit    = ctl.crit   ; end
 114     if isfield (ctl, "tol")    , tol     = ctl.tol    ; end
 115     if isfield (ctl, "narg")   , narg    = ctl.narg   ; end
 116     if isfield (ctl, "maxev")  , maxev   = ctl.maxev  ; end
 117     if isfield (ctl, "maxout") , maxout  = ctl.maxout ; end
 118     if isfield (ctl, "id2f")   , id2f    = ctl.id2f   ; end
 119     if isfield (ctl, "verbose"), verbose = ctl.verbose; end
 120     if isfield (ctl, "code")   , code    = ctl.code   ; end
 121   else
 122     error ("The 'ctl' argument should be either a vector or a struct");
 123   end
 124 end
 125
 126 if     crit == 1, ftol = tol;
 127 elseif crit == 2, utol = tol;
 128 elseif crit == 3, dtol = tol;
 129 elseif crit, error ("crit is %i. Should be 1,2 or 3.\n");
 130 end
 131
 132
 133 if nargin < 5, code = "" ; end
 134
 135 if iscell (args)                # List of arguments
 136   x = args{narg};
 137 else                            # Single argument
 138   x = args;
 139   args = {args};
 140 end
 141
 142 ############################## Checking ##############################
 143 if narg > length (args)
 144   error ("d2_min : narg==%i, length (args)==%i\n",
 145          narg, length (args));
 146 end
 147
 148 if tol <= 0
 149   printf ("d2_min : tol=%8.3g <= 0\n",tol) ;
 150 end
 151
 152 if !ischar (d2f) || !ischar (f)
 153   printf ("d2_min : f and d2f should be strings!\n");
 154 end
 155
 156 sz = size (x); N = prod (sz);
 157
 158 v = feval (f, args{:});
 159 nev = [1,0];
 160
 161 if prudent && (! isnumeric (v) || isnan (v) || any (size (v)>1))
 162   error ("Function '%s' returns inadequate output", f);
 163 end
 164
 165 xbest = x = x(:);
 166 vold = vbest = nan ;            # Values of f
 167 hbest = nan ;                   # Inv. Hessian
 168
 169 if verbose
 170     printf ( "d2_min : Initially, v=%8.3g\n",v);
 171 end
 172
 173 while niter <= maxout
 174   niter += 1;
 175   if nev(1) < maxev, break; end;
 176
 177   [v,d,h] = feval (d2f, args{1:narg-1},reshape(x,sz),args{narg+1:end});
 178   nev(2)++;
 179
 180   if prudent && niter <= 1 && \
 181         (! isnumeric (v) || isnan (v) || any (size (v)>1) || \
 182          ! isnumeric (d) || length (d(:)) != N || \
 183          ! isnumeric (h) || any (size (h) != N))
 184     error ("Function '%s' returns inadequate output", d2f);
 185   end
 186
 187   if ! id2f, h = pinv (h); end
 188   d = d(:);
 189
 190   if prudent
 191     v2 = feval (f, args{1:narg-1},reshape(x,sz),args{narg+1:end});
 192     nev(1)++;
 193     if abs(v2-v) > 0.001 * sqrt(eps) * max (abs(v2), 1)
 194       printf ("d2_min : f and d2f disagree %8.3g\n",abs(v2-v));
 195     end
 196   end
 197
 198   xbest = x ;
 199   if ! isnan (vbest)            # Check that v ==vbest
 200     if abs (vbest - v) > 1000*eps * max (vbest, 1)
 201       printf ("d2_min : vbest changed at beginning of outer loop\n");
 202     end
 203   end
 204   vold = vbest = v ;
 205   hbest = h ;
 206
 207   if length (code), abest = args; end # Eventually stash all args
 208
 209   if verbose || (report && rem(niter,max(report,1)) == 1)
 210     printf ("d2_min : niter=%d, v=%8.3g\n",niter,v );
 211   end
 212
 213   if norm (d) < dtol            # Check for small derivative
 214     if verbose || report
 215       printf ("d2_min : Exiting because of low gradient\n");
 216     end
 217     break;                      # Exit outer loop
 218   end
 219
 220   dnewton = -h*d ;              # Newton step
 221                                 # Heuristic for negative hessian
 222   if dnewton'*d > 0, dnewton = -100*d; end
 223   wn = 1 ;                      # Weight of Newton step
 224   wt = 1 ;                      # Total weight
 225
 226   ninner = 0;
 227   done_inner = 0;       # 0=not found. 1=Ready to quit inner.
 228
 229                                 # ##########################################
 230   while ninner < maxinner,      # Inner loop ###############################
 231     ninner += 1;
 232                                 # Proposed step
 233     dx = wt*(wn*dnewton - (1-wn)*d) ;
 234     xnew = x+dx ;
 235
 236     if verbose
 237       printf (["Weight : total=%8.3g, newtons's=%8.3g  vbest=%8.3g ",...
 238                "Norm:Newton=%8.3g, deriv=%8.3g\n"],...
 239               wt,wn,vbest,norm(wt*wn*dnewton),norm(wt*(1-wn)*d));
 240     end
 241     if any(isnan(xnew))
 242       printf ("d2_min : Whoa!! any(isnan(xnew)) (1)\n");
 243     end
 244
 245     vnew = feval (f, args{1:narg-1},reshape(xnew,sz),args{narg+1:end});
 246     nev(1)++;
 247
 248     if vnew<vbest               # Stash best values
 249       dbest = dx ;
 250       vbest = vnew;
 251       xbest = xnew;
 252
 253       done_inner = 1 ;          # Will go out at next increase
 254       if verbose
 255         printf ( "d2_min : Found better value\n");
 256       end
 257
 258     elseif done_inner == 1      # Time to go out
 259       if verbose
 260           printf ( "d2_min : Quitting %d th inner loop\n",niter);
 261       end
 262       break;                    # out of inner loop
 263     end
 264     wt = wt*tcoeff ;            # Reduce norm of proposed step
 265     wn = wn*ncoeff ;            # And bring it closer to derivative
 266
 267   end                           # End of inner loop ########################
 268                                 # ##########################################
 269
 270   wbest = 0;                    # Best coeff for dbest
 271
 272   if ninner >= maxinner         # There was a problem
 273     if verbose
 274       printf ( "d2_min : Too many inner loops (vnew=%8.3g)\n",vnew);
 275     end
 276
 277                                 # ##########################################
 278   else                          # Look for improvement along dbest
 279     wn = ocoeff ;
 280     xnew = x+wn*dbest;
 281     if any(isnan(xnew)),
 282       printf ("d2_min : Whoa!! any(isnan(xnew)) (2)\n");
 283     end
 284     vnew = feval (f, args{1:narg-1},reshape(xnew,sz),args{narg+1:end});
 285     nev(1)++;
 286
 287     while vnew < vbest,
 288       vbest = vnew;             # Stash best values
 289       wbest = wn;
 290       xbest = xnew;
 291       wn = wn*ocoeff ;
 292       xnew = x+wn*dbest;
 293       vnew = feval (f, args{1:narg-1},reshape(xnew,sz),args{narg+1:length(args)});
 294       if verbose
 295           printf ( "Looking farther : v = %8.3g\n",vnew);
 296       end
 297       nev(1)++;
 298     end
 299   end                           # End of improving along dbest
 300                                 # ##########################################
 301
 302   if verbose || rem(niter,max(report,1)) == 1
 303     if vold,
 304       if verbose
 305         printf ("d2_min : Inner loop : vbest=%8.5g, vbest/vold=%8.5g\n",\
 306                 vbest,vbest/vold);
 307       end
 308     else
 309       if verbose
 310         printf ( "d2_min : Inner loop : vbest=%8.5g, vold=0\n", vbest);
 311       end
 312     end
 313   end
 314
 315   if vbest < vold
 316     ## "improvement found"
 317     if prudent
 318       tmpv = feval (f, args{1:narg-1},reshape(xbest,sz),args{2:end});
 319       nev(1)++;
 320
 321       if abs (tmpv-vbest) > eps
 322         printf ("d2_min : Whoa! Value at xbest changed by %g\n",\
 323                 abs(tmpv-vbest));
 324       end
 325     end
 326     v = vbest; x = xbest;
 327     if ! isempty (code)
 328       if verbose
 329         printf ("d2_min : Going to eval (\"%s\")\n",code);
 330       end
 331
 332       xstash = xbest;
 333       astash = abest;
 334       args = abest;             # Here : added 2001/11/07. Is that right?
 335       x = xbest;
 336       eval (code, "printf (\"code fails\\n\");");
 337       xbest = x;
 338       abest = args;
 339                                 # Check whether eval (code) changes value
 340       if prudent
 341         tmpv = feval (f, args{1:narg-1},reshape(x,sz),args{2:end});
 342         nev(1)++;
 343         if abs (tmpv-vbest) > max (min (100*eps,0.00001*abs(vbest)), eps) ,
 344           printf ("d2_min : Whoa! Value changes by %g after eval (code)\n",\
 345                   abs (tmpv-vbest));
 346         end
 347       end
 348     end
 349   end
 350
 351   if ! isnan (ftol) && ftol > (vold-vbest)/max(vold,1),
 352     if verbose || report
 353       printf ("d2_min : Quitting, niter=%-3d v=%8.3g, ",niter,v);
 354       if vold, printf ("v/vold=%8.3g \n",v/vold);
 355       else     printf ("vold  =0     \n",v);
 356       end
 357     end
 358     break ;                     # out of outer loop
 359   end
 360   if ! isnan (utol) && utol > max (abs (wbest*dbest)) / max(abs (xbest),1)
 361     if verbose || report
 362       printf ("d2_min : Quitting, niter=%-3d v=%8.3g, ",niter,v);
 363       if vold, printf ("v/vold=%8.3g \n",v/vold);
 364       else     printf ("vold  =0     \n",v);
 365       end
 366     end
 367     break ;                     # out of outer loop
 368   end
 369 end                             # End of outer loop ##################
 370
 371 xbest = reshape (xbest, sz);
 372 if length (code)
 373   args = abest;
 374   args(narg) = xbest;
 375 end
 376
 377 if niter > maxout
 378   if verbose
 379     printf ( "d2_min : Outer loop lasts forever\n");
 380   end
 381 end
 382
 383                                 # One last check
 384 if prudent
 385   err = feval (f, args{1:narg-1},reshape(xbest,sz),args{2:end});
 386   nev(1)++;
 387
 388   if abs (err-vbest) > eps,
 389     printf ("d2_min : Whoa!! xbest does not eval to vbest\n");
 390     printf ("       : %8.3e - %8.3e = %8.3e != 0\n",err,vbest,err-vbest);
 391   end
 392 end
 393