octave_packages/optim-1.2.0/powell.m

   1 ## Copyright (C) 2011 Nir Krakauer <nkrakauer@ccny.cuny.edu>
   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{p}, @var{obj_value}, @var{convergence}, @var{iters}, @var{nevs}] = powell (@var{f}, @var{p0}, @var{control})
  18 ## Multidimensional minimization (direction-set method). Implements a direction-set (Powell's) method for multidimensional minimization of a function without calculation of the gradient [1, 2]
  19 ##
  20 ## @subheading Arguments
  21 ##
  22 ## @itemize @bullet
  23 ## @item
  24 ## @var{f}: name of function to minimize (string or handle), which should accept one input variable (see example for how to pass on additional input arguments)
  25 ##
  26 ## @item
  27 ## @var{p0}: An initial value of the function argument to minimize
  28 ##
  29 ## @item
  30 ## @var{options}: an optional structure, which can be generated by optimset, with some or all of the following fields:
  31 ## @itemize @minus
  32 ## @item
  33 ##      MaxIter: maximum iterations  (positive integer, or -1 or Inf for unlimited (default))
  34 ## @item
  35 ##      TolFun: minimum amount by which function value must decrease in each iteration to continue (default is 1E-8)
  36 ## @item
  37 ##      MaxFunEvals: maximum function evaluations  (positive integer, or -1 or Inf for unlimited (default))
  38 ## @item
  39 ##      SearchDirections: an n*n matrix whose columns contain the initial set of (presumably orthogonal) directions to minimize along, where n is the number of elements in the argument to be minimized for; or an n*1 vector of magnitudes for the initial directions (defaults to the set of unit direction vectors)
  40 ## @end itemize
  41 ## @end itemize
  42 ##
  43 ## @subheading Examples
  44 ##
  45 ## @example
  46 ## @group
  47 ## y = @@(x, s) x(1) ^ 2 + x(2) ^ 2 + s;
  48 ## o = optimset('MaxIter', 100, 'TolFun', 1E-10);
  49 ## s = 1;
  50 ## [x_optim, y_min, conv, iters, nevs] = powell(@@(x) y(x, s), [1 0.5], o); %pass y wrapped in an anonymous function so that all other arguments to y, which are held constant, are set
  51 ## %should return something like x_optim = [4E-14 3E-14], y_min = 1, conv = 1, iters = 2, nevs = 24
  52 ## @end group
  53 ##
  54 ## @end example
  55 ##
  56 ## @subheading Returns:
  57 ##
  58 ## @itemize @bullet
  59 ## @item
  60 ## @var{p}: the minimizing value of the function argument
  61 ## @item
  62 ## @var{obj_value}: the value of @var{f}() at @var{p}
  63 ## @item
  64 ## @var{convergence}: 1 if normal convergence, 0 if not
  65 ## @item
  66 ## @var{iters}: number of iterations performed
  67 ## @item
  68 ## @var{nevs}: number of function evaluations
  69 ## @end itemize
  70 ##
  71 ## @subheading References
  72 ##
  73 ## @enumerate
  74 ## @item
  75 ## Powell MJD (1964), An efficient method for finding the minimum of a function of several variables without calculating derivatives, @cite{Computer Journal}, 7 :155-162
  76 ##
  77 ## @item
  78 ## Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (1992). @cite{Numerical Recipes in Fortran: The Art of Scientific Computing} (2nd Ed.). New York: Cambridge University Press (Section 10.5)
  79 ## @end enumerate
  80 ## @end deftypefn
  81
  82 ## PKG_ADD: __all_opts__ ("powell");
  83
  84 function [p, obj_value, convergence, iters, nevs] = powell (f, p0, options);
  85
  86   if (nargin == 1 && ischar (f) && strcmp (f, "defaults"))
  87     p = optimset ("MaxIter", Inf, \
  88                   "TolFun", 1e-8, \
  89                   "MaxFunEvals", Inf, \
  90                   "SearchDirections", []);
  91     return;
  92   endif
  93
  94   ## check number of arguments
  95   if ((nargin < 2) || (nargin > 3))
  96     usage('powell: you must supply 2 or 3 arguments');
  97   endif
  98
  99
 100   ## default or input values
 101
 102   if (nargin < 3)
 103     options = struct ();
 104   endif
 105
 106   xi_set = 0;
 107   xi = optimget (options, 'SearchDirections');
 108   if (! isempty (xi))
 109     if (isvector (xi)) # assume that xi is is n*1 or 1*n
 110       xi = diag (x);
 111     endif
 112     xi_set = 1;
 113   endif
 114
 115
 116   MaxIter = optimget (options, 'MaxIter', Inf);
 117   if (MaxIter < 0) MaxIter = Inf; endif
 118   MaxFunEvals = optimget (options, 'MaxFunEvals', Inf);
 119   TolFun = optimget (options, 'TolFun', 1E-8);
 120
 121   nevs = 0;
 122   iters = 0;
 123   convergence = 0;
 124
 125   p = p0; # initial value of the argument being minimized
 126
 127   try
 128     obj_value = f(p);
 129   catch
 130     error ("function does not exist or cannot be evaluated");
 131   end_try_catch
 132
 133   nevs++;
 134
 135   n = numel (p); # number of dimensions to minimize over
 136
 137   xit = zeros (n, 1);
 138   if (! xi_set)
 139     xi = eye(n);
 140   endif
 141
 142
 143
 144   ## do an iteration
 145   while (iters <= MaxIter && nevs <= MaxFunEvals && ! convergence)
 146     iters++;
 147     pt = p; # best point as iteration begins
 148     fp = obj_value; # value of the objective function as iteration begins
 149     ibig = 0; # will hold direction along which the objective function decreased the most in this iteration
 150     dlt = 0; # will hold decrease in objective function value in this iteration
 151     for i = 1:n
 152       xit = reshape (xi(:, i), size(p));
 153       fptt = obj_value;
 154       [a, obj_value, nev] = line_min (f, xit, {p});
 155       nevs = nevs + nev;
 156       p = p + a*xit;
 157       change = fptt - obj_value;
 158       if (change > dlt)
 159         dlt = change;
 160         ibig = i;
 161       endif
 162     endfor
 163
 164     if ( 2*abs(fp-obj_value) <= TolFun*(abs(fp) + abs(obj_value)) )
 165       convergence = 1;
 166       return
 167     endif
 168
 169     if (iters == MaxIter)
 170       disp ("iteration maximum exceeded");
 171       return
 172     endif
 173
 174     ## attempt parabolic extrapolation
 175     ptt = 2*p - pt;
 176     xit = p - pt;
 177     fptt = f(ptt);
 178     nevs++;
 179     if (fptt < fp) # check whether the extrapolation actually makes the objective function smaller
 180       t = 2 * (fp - 2*obj_value + fptt) * (fp-obj_value-dlt)^2 - dlt * (fp-fptt)^2;
 181       if (t < 0)
 182         p = ptt;
 183         [a, obj_value, nev] = line_min (f, xit, {p});
 184         nevs = nevs + nev;
 185         p = p + a*xit;
 186
 187         ## add the net direction from this iteration to the direction set
 188         xi(:, ibig) = xi(:, n);
 189         xi(:, n) = xit(:);
 190       endif
 191     endif
 192   endwhile
 193