octave_packages/optim-1.2.0/vfzero.m

   1 ## Copyright (C) 2008, 2009 VZLU Prague, a.s.
   2 ## Copyright (C) 2010 Olaf Till <olaf.till@uni-jena.de>
   3 ##
   4 ## This program is free software; you can redistribute it and/or modify it under
   5 ## the terms of the GNU General Public License as published by the Free Software
   6 ## Foundation; either version 3 of the License, or (at your option) any later
   7 ## version.
   8 ##
   9 ## This program is distributed in the hope that it will be useful, but WITHOUT
  10 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  11 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  12 ## details.
  13 ##
  14 ## You should have received a copy of the GNU General Public License along with
  15 ## this program; if not, see <http://www.gnu.org/licenses/>.
  16
  17 ## -*- texinfo -*-
  18 ## @deftypefn  {Function File} {} vfzero (@var{fun}, @var{x0})
  19 ## @deftypefnx {Function File} {} vfzero (@var{fun}, @var{x0}, @var{options})
  20 ## @deftypefnx {Function File} {[@var{x}, @var{fval}, @var{info}, @var{output}] =} vfzero (@dots{})
  21 ## A variant of @code{fzero}. Finds a zero of a vector-valued
  22 ## multivariate function where each output element only depends on the
  23 ## input element with the same index (so the Jacobian is diagonal).
  24 ##
  25 ## @var{fun} should be a handle or name of a function returning a column
  26 ## vector.  @var{x0} should be a two-column matrix, each row specifying
  27 ## two points which bracket a zero of the respective output element of
  28 ## @var{fun}.
  29 ##
  30 ## If @var{x0} is a single-column matrix then several nearby and distant
  31 ## values are probed in an attempt to obtain a valid bracketing.  If
  32 ## this is not successful, the function fails. @var{options} is a
  33 ## structure specifying additional options. Currently, @code{vfzero}
  34 ## recognizes these options: @code{"FunValCheck"}, @code{"OutputFcn"},
  35 ## @code{"TolX"}, @code{"MaxIter"}, @code{"MaxFunEvals"}. For a
  36 ## description of these options, see @ref{doc-optimset,,optimset}.
  37 ##
  38 ## On exit, the function returns @var{x}, the approximate zero and
  39 ## @var{fval}, the function value thereof. @var{info} is a column vector
  40 ## of exit flags  that can have these values:
  41 ##
  42 ## @itemize
  43 ## @item 1 The algorithm converged to a solution.
  44 ##
  45 ## @item 0 Maximum number of iterations or function evaluations has been
  46 ## reached.
  47 ##
  48 ## @item -1 The algorithm has been terminated from user output function.
  49 ##
  50 ## @item -5 The algorithm may have converged to a singular point.
  51 ## @end itemize
  52 ##
  53 ## @var{output} is a structure containing runtime information about the
  54 ## @code{fzero} algorithm.  Fields in the structure are:
  55 ##
  56 ## @itemize
  57 ## @item iterations Number of iterations through loop.
  58 ##
  59 ## @item nfev Number of function evaluations.
  60 ##
  61 ## @item bracketx A two-column matrix with the final bracketing of the
  62 ## zero along the x-axis.
  63 ##
  64 ## @item brackety A two-column matrix with the final bracketing of the
  65 ## zero along the y-axis.
  66 ## @end itemize
  67 ## @seealso{optimset, fsolve}
  68 ## @end deftypefn
  69
  70 ## This is essentially the ACM algorithm 748: Enclosing Zeros of
  71 ## Continuous Functions due to Alefeld, Potra and Shi, ACM Transactions
  72 ## on Mathematical Software, Vol. 21, No. 3, September 1995. Although
  73 ## the workflow should be the same, the structure of the algorithm has
  74 ## been transformed non-trivially; instead of the authors' approach of
  75 ## sequentially calling building blocks subprograms we implement here a
  76 ## FSM version using one interior point determination and one bracketing
  77 ## per iteration, thus reducing the number of temporary variables and
  78 ## simplifying the algorithm structure. Further, this approach reduces
  79 ## the need for external functions and error handling. The algorithm has
  80 ## also been slightly modified.
  81
  82 ## Author: Jaroslav Hajek <highegg@gmail.com>
  83
  84 ## PKG_ADD: __all_opts__ ("vfzero");
  85
  86 function [x, fval, info, output] = vfzero (fun, x0, options = struct ())
  87
  88   ## Get default options if requested.
  89   if (nargin == 1 && ischar (fun) && strcmp (fun, 'defaults'))
  90     x = optimset ("MaxIter", Inf, "MaxFunEvals", Inf, "TolX", 1e-8, \
  91     "OutputFcn", [], "FunValCheck", "off");
  92     return;
  93   endif
  94
  95   if (nargin < 2 || nargin > 3)
  96     print_usage ();
  97   endif
  98
  99   if (ischar (fun))
 100     fun = str2func (fun, "global");
 101   endif
 102
 103   ## TODO
 104   ## displev = optimget (options, "Display", "notify");
 105   funvalchk = strcmpi (optimget (options, "FunValCheck", "off"), "on");
 106   outfcn = optimget (options, "OutputFcn");
 107   tolx = optimget (options, "TolX", 1e-8);
 108   maxiter = optimget (options, "MaxIter", Inf);
 109   maxfev = optimget (options, "MaxFunEvals", Inf);
 110   nx = rows (x0);
 111   ## fun may assume a certain length of x, so we will always call it
 112   ## with the full-length x, even if only some elements are needed
 113
 114   persistent mu = 0.5;
 115
 116   if (funvalchk)
 117     ## Replace fun with a guarded version.
 118     fun = @(x) guarded_eval (fun, x);
 119   endif
 120
 121   ## The default exit flag if exceeded number of iterations.
 122   info = zeros (nx, 1);
 123   niter = 0;
 124   nfev = 0;
 125
 126   x = fval = fc = a = fa = b = fb = aa = c = u = fu = NaN (nx, 1);
 127   bracket_ready = false (nx, 1);
 128   eps = eps (class (x0));
 129
 130   ## Prepare...
 131   a = x0(:, 1);
 132   fa = fun (a)(:);
 133   nfev = 1;
 134   if (columns (x0) > 1)
 135     b = x0(:, 2);
 136     fb = fun (b)(:);
 137     nfev += 1;
 138   else
 139     ## Try to get b.
 140     aa(idx = a == 0) = 1;
 141     aa(! idx) = a(! idx);
 142     for tb = [0.9*aa, 1.1*aa, aa-1, aa+1, 0.5*aa 1.5*aa, -aa, 2*aa, -10*aa, 10*aa]
 143       tfb = fun (tb)(:); nfev += 1;
 144       idx = ! bracket_ready & sign (fa) .* sign (tfb) <= 0;
 145       bracket_ready |= idx;
 146       b(idx) = tb(idx);
 147       fb(idx) = tfb(idx);
 148       if (all (bracket_ready))
 149         break;
 150       endif
 151     endfor
 152   endif
 153
 154   tp = a(idx = b < a);
 155   a(idx) = b(idx);
 156   b(idx) = tp;
 157
 158   tp = fa(idx);
 159   fa(idx) = fb(idx);
 160   fb(idx) = tp;
 161
 162   if (! all (sign (fa) .* sign (fb) <= 0))
 163     error ("fzero:bracket", "vfzero: not a valid initial bracketing");
 164   endif
 165
 166   slope0 = (fb - fa) ./ (b - a);
 167
 168   idx = fa == 0;
 169   b(idx) = a(idx);
 170   fb(idx) = fa(idx);
 171
 172   idx = (! idx & fb == 0);
 173   a(idx) = b(idx);
 174   fa(idx) = fb(idx);
 175
 176   itype = ones (nx, 1);
 177
 178   idx = abs (fa) < abs (fb);
 179   u(idx) = a(idx); fu(idx) = fa(idx);
 180   u(! idx) = b(! idx); fu(! idx) = fb(! idx);
 181
 182   d = e = u;
 183   fd = fe = fu;
 184   mba = mu * (b - a);
 185   not_ready = true (nx, 1);
 186   while (niter < maxiter && nfev < maxfev && any (not_ready))
 187
 188     ## itype == 1
 189     type1idx = not_ready & itype == 1;
 190     ## The initial test.
 191     idx = b - a <= 2*(2 * eps * abs (u) + tolx) & type1idx;
 192     x(idx) = u(idx); fval(idx) = fu(idx);
 193     info(idx) = 1;
 194     not_ready(idx) = false;
 195     type1idx &= not_ready;
 196     exclidx = type1idx;
 197     ## Secant step.
 198     idx = type1idx & \
 199         (tidx = abs (fa) <= 1e3*abs (fb) & abs (fb) <= 1e3*abs (fa));
 200     c(idx) = u(idx) - (a(idx) - b(idx)) ./ (fa(idx) - fb(idx)) .* fu(idx);
 201     ## Bisection step.
 202     idx = type1idx & ! tidx;
 203     c(idx) = 0.5*(a(idx) + b(idx));
 204     d(type1idx) = u(type1idx); fd(type1idx) = fu(type1idx);
 205     itype(type1idx) = 5;
 206
 207     ## itype == 2 or 3
 208     type23idx = not_ready & ! exclidx & (itype == 2 | itype == 3);
 209     exclidx |= type23idx;
 210     uidx = cellfun (@ (x) length (unique (x)), \
 211                     num2cell ([fa, fb, fd, fe], 2)) == 4;
 212     oidx = sign (c - a) .* sign (c - b) > 0;
 213     ## Inverse cubic interpolation.
 214     idx = type23idx & (uidx & ! oidx);
 215     q11 = (d(idx) - e(idx)) .* fd(idx) ./ (fe(idx) - fd(idx));
 216     q21 = (b(idx) - d(idx)) .* fb(idx) ./ (fd(idx) - fb(idx));
 217     q31 = (a(idx) - b(idx)) .* fa(idx) ./ (fb(idx) - fa(idx));
 218     d21 = (b(idx) - d(idx)) .* fd(idx) ./ (fd(idx) - fb(idx));
 219     d31 = (a(idx) - b(idx)) .* fb(idx) ./ (fb(idx) - fa(idx));
 220     q22 = (d21 - q11) .* fb(idx) ./ (fe(idx) - fb(idx));
 221     q32 = (d31 - q21) .* fa(idx) ./ (fd(idx) - fa(idx));
 222     d32 = (d31 - q21) .* fd(idx) ./ (fd(idx) - fa(idx));
 223     q33 = (d32 - q22) .* fa(idx) ./ (fe(idx) - fa(idx));
 224     c(idx) = a(idx) + q31 + q32 + q33;
 225     ## Quadratic interpolation + newton.
 226     idx = type23idx & (oidx | ! uidx);
 227     a0 = fa(idx);
 228     a1 = (fb(idx) - fa(idx))./(b(idx) - a(idx));
 229     a2 = ((fd(idx) - fb(idx))./(d(idx) - b(idx)) - a1) ./ (d(idx) - a(idx));
 230     ## Modification 1: this is simpler and does not seem to be worse.
 231     c(idx) = a(idx) - a0./a1;
 232     taidx = a2 != 0;
 233     tidx = idx;
 234     tidx(tidx) = taidx;
 235     c(tidx) = a(tidx)(:) - (a0(taidx)./a1(taidx))(:);
 236     for i = 1:3
 237       tidx &= i <= itype;
 238       taidx = tidx(idx);
 239       pc = a0(taidx)(:) + (a1(taidx)(:) + \
 240                            a2(taidx)(:).*(c(tidx) - b(tidx))(:)) \
 241           .*(c(tidx) - a(tidx))(:);
 242       pdc = a1(taidx)(:) + a2(taidx)(:).*(2*c(tidx) - a(tidx) - b(tidx))(:);
 243       tidx0 = tidx;
 244       tidx0(tidx0, 1) &= (p0idx = pdc == 0);
 245       taidx0 = tidx0(idx);
 246       tidx(tidx, 1) &= ! p0idx;
 247       c(tidx0) = a(tidx0)(:) - (a0(taidx0)./a1(taidx0))(:);
 248       c(tidx) = c(tidx)(:) - (pc(! p0idx)./pdc(! p0idx))(:);
 249     endfor
 250     itype(type23idx) += 1;
 251
 252     ## itype == 4
 253     type4idx = not_ready & ! exclidx & itype == 4;
 254     exclidx |= type4idx;
 255     ## Double secant step.
 256     idx = type4idx;
 257     c(idx) = u(idx) - 2*(b(idx) - a(idx))./(fb(idx) - fa(idx)).*fu(idx);
 258     ## Bisect if too far.
 259     idx = type4idx & abs (c - u) > 0.5*(b - a);
 260     c(idx) = 0.5 * (b(idx) + a(idx));
 261     itype(type4idx) = 5;
 262
 263     ## itype == 5
 264     type5idx = not_ready & ! exclidx & itype == 5;
 265     ## Bisection step.
 266     idx = type5idx;
 267     c(idx) = 0.5 * (b(idx) + a(idx));
 268     itype(type5idx) = 2;
 269
 270     ## Don't let c come too close to a or b.
 271     delta = 2*0.7*(2 * eps * abs (u) + tolx);
 272     nidx = not_ready & ! (idx = b - a <= 2*delta);
 273     idx &= not_ready;
 274     c(idx) = (a(idx) + b(idx))/2;
 275     c(nidx) = max (a(nidx) + delta(nidx), \
 276                    min (b(nidx) - delta(nidx), c(nidx)));
 277
 278     ## Calculate new point.
 279     idx = not_ready;
 280     x(idx, 1) = c(idx, 1);
 281     if (any (idx))
 282       c(! idx) = u(! idx); # to have some working place-holders since
 283                                 # fun() might expect full-length
 284                                 # argument
 285       fval(idx, 1) = fc(idx, 1) = fun (c)(:)(idx, 1);
 286       niter ++; nfev ++;
 287     endif
 288
 289     ## Modification 2: skip inverse cubic interpolation if
 290     ## nonmonotonicity is detected.
 291     nidx = not_ready & ! (idx = sign (fc - fa) .* sign (fc - fb) >= 0);
 292     idx &= not_ready;
 293     ## The new point broke monotonicity.
 294     ## Disable inverse cubic.
 295     fe(idx) = fc(idx);
 296     ##
 297     e(nidx) = d(nidx); fe(nidx) = fd(nidx);
 298
 299     ## Bracketing.
 300     idx1 = not_ready & sign (fa) .* sign (fc) < 0;
 301     idx2 = not_ready & ! idx1 & sign (fb) .* sign (fc) < 0;
 302     idx3 = not_ready & ! (idx1 | idx2) & fc == 0;
 303     d(idx1) = b(idx1); fd(idx1) = fb(idx1);
 304     b(idx1) = c(idx1); fb(idx1) = fc(idx1);
 305     d(idx2) = a(idx2); fd(idx2) = fa(idx2);
 306     a(idx2) = c(idx2); fa(idx2) = fc(idx2);
 307     a(idx3) = b(idx3) = c(idx3); fa(idx3) = fb(idx3) = fc(idx3);
 308     info(idx3) = 1;
 309     not_ready(idx3) = false;
 310     if (any (not_ready & ! (idx1 | idx2 | idx3)))
 311       ## This should never happen.
 312       error ("fzero:bracket", "vfzero: zero point is not bracketed");
 313     endif
 314
 315     ## If there's an output function, use it now.
 316     if (! isempty (outfcn))
 317       optv.funccount = nfev;
 318       optv.fval = fval;
 319       optv.iteration = niter;
 320       idx = not_ready & outfcn (x, optv, "iter");
 321       info(idx) = -1;
 322       not_ready(idx) = false;
 323     endif
 324
 325     nidx = not_ready & ! (idx = abs (fa) < abs (fb));
 326     idx &= not_ready;
 327     u(idx) = a(idx); fu(idx) = fa(idx);
 328     u(nidx) = b(nidx); fu(nidx) = fb(nidx);
 329     idx = not_ready & b - a <= 2*(2 * eps * abs (u) + tolx);
 330     info(idx) = 1;
 331     not_ready(idx) = false;
 332
 333     ## Skip bisection step if successful reduction.
 334     itype(not_ready & itype == 5 & (b - a) <= mba) = 2;
 335     idx = not_ready & itype == 2;
 336     mba(idx) = mu * (b(idx) - a(idx));
 337   endwhile
 338
 339   ## Check solution for a singularity by examining slope
 340   idx = not_ready & info == 1 & (b - a) != 0;
 341   idx(idx, 1) &= \
 342       abs ((fb(idx, 1) - fa(idx, 1))./(b(idx, 1) - a(idx, 1)) \
 343            ./ slope0(idx, 1)) > max (1e6, 0.5/(eps+tolx));
 344   info(idx) = - 5;
 345
 346   output.iterations = niter;
 347   output.funcCount = nfev;
 348   output.bracketx = [a, b];
 349   output.brackety = [fa, fb];
 350
 351 endfunction
 352
 353 ## An assistant function that evaluates a function handle and checks for
 354 ## bad results.
 355 function fx = guarded_eval (fun, x)
 356   fx = fun (x);
 357   if (! isreal (fx))
 358     error ("fzero:notreal", "vfzero: non-real value encountered");
 359   elseif (any (isnan (fx)))
 360     error ("fzero:isnan", "vfzero: NaN value encountered");
 361   endif
 362 endfunction
 363
 364 %!shared opt0
 365 %! opt0 = optimset ("tolx", 0);
 366 %!assert(vfzero(@cos, [0, 3], opt0), pi/2, 10*eps)
 367 %!assert(vfzero(@(x) x^(1/3) - 1e-8, [0,1], opt0), 1e-24, 1e-22*eps)