PageRenderTime 95ms CodeModel.GetById 2ms app.highlight 87ms RepoModel.GetById 1ms app.codeStats 0ms

/octave-3.6.2/scripts/optimization/fsolve.m

#
Objective C | 610 lines | 546 code | 64 blank | 0 comment | 81 complexity | 217bc77cb0632a851ae5b2c5e2388a7f MD5 | raw file
  1## Copyright (C) 2008-2012 VZLU Prague, a.s.
  2##
  3## This file is part of Octave.
  4##
  5## Octave is free software; you can redistribute it and/or modify it
  6## under the terms of the GNU General Public License as published by
  7## the Free Software Foundation; either version 3 of the License, or (at
  8## your option) any later version.
  9##
 10## Octave is distributed in the hope that it will be useful, but
 11## WITHOUT ANY WARRANTY; without even the implied warranty of
 12## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 13## General Public License for more details.
 14##
 15## You should have received a copy of the GNU General Public License
 16## along with Octave; see the file COPYING.  If not, see
 17## <http://www.gnu.org/licenses/>.
 18##
 19## Author: Jaroslav Hajek <highegg@gmail.com>
 20
 21## -*- texinfo -*-
 22## @deftypefn  {Function File} {} fsolve (@var{fcn}, @var{x0}, @var{options})
 23## @deftypefnx {Function File} {[@var{x}, @var{fvec}, @var{info}, @var{output}, @var{fjac}] =} fsolve (@var{fcn}, @dots{})
 24## Solve a system of nonlinear equations defined by the function @var{fcn}.
 25## @var{fcn} should accept a vector (array) defining the unknown variables,
 26## and return a vector of left-hand sides of the equations.  Right-hand sides
 27## are defined to be zeros.
 28## In other words, this function attempts to determine a vector @var{x} such
 29## that @code{@var{fcn} (@var{x})} gives (approximately) all zeros.
 30## @var{x0} determines a starting guess.  The shape of @var{x0} is preserved
 31## in all calls to @var{fcn}, but otherwise it is treated as a column vector.
 32## @var{options} is a structure specifying additional options.
 33## Currently, @code{fsolve} recognizes these options:
 34## @code{"FunValCheck"}, @code{"OutputFcn"}, @code{"TolX"},
 35## @code{"TolFun"}, @code{"MaxIter"}, @code{"MaxFunEvals"},
 36## @code{"Jacobian"}, @code{"Updating"}, @code{"ComplexEqn"}
 37## @code{"TypicalX"}, @code{"AutoScaling"} and @code{"FinDiffType"}.
 38##
 39## If @code{"Jacobian"} is @code{"on"}, it specifies that @var{fcn},
 40## called with 2 output arguments, also returns the Jacobian matrix
 41## of right-hand sides at the requested point.  @code{"TolX"} specifies
 42## the termination tolerance in the unknown variables, while
 43## @code{"TolFun"} is a tolerance for equations.  Default is @code{1e-7}
 44## for both @code{"TolX"} and @code{"TolFun"}.
 45##
 46## If @code{"AutoScaling"} is on, the variables will be automatically scaled
 47## according to the column norms of the (estimated) Jacobian.  As a result,
 48## TolF becomes scaling-independent.  By default, this option is off, because
 49## it may sometimes deliver unexpected (though mathematically correct) results.
 50##
 51## If @code{"Updating"} is "on", the function will attempt to use Broyden
 52## updates to update the Jacobian, in order to reduce the amount of Jacobian
 53## calculations.
 54## If your user function always calculates the Jacobian (regardless of number
 55## of output arguments), this option provides no advantage and should be set to
 56## false.
 57##
 58## @code{"ComplexEqn"} is @code{"on"}, @code{fsolve} will attempt to solve
 59## complex equations in complex variables, assuming that the equations possess a
 60## complex derivative (i.e., are holomorphic).  If this is not what you want,
 61## should unpack the real and imaginary parts of the system to get a real
 62## system.
 63##
 64## For description of the other options, see @code{optimset}.
 65##
 66## On return, @var{fval} contains the value of the function @var{fcn}
 67## evaluated at @var{x}, and @var{info} may be one of the following values:
 68##
 69## @table @asis
 70## @item 1
 71## Converged to a solution point.  Relative residual error is less than
 72## specified by TolFun.
 73##
 74## @item 2
 75## Last relative step size was less that TolX.
 76##
 77## @item 3
 78## Last relative decrease in residual was less than TolF.
 79##
 80## @item 0
 81## Iteration limit exceeded.
 82##
 83## @item -3
 84## The trust region radius became excessively small.
 85## @end table
 86##
 87## Note: If you only have a single nonlinear equation of one variable, using
 88## @code{fzero} is usually a much better idea.
 89##
 90## Note about user-supplied Jacobians:
 91## As an inherent property of the algorithm, Jacobian is always requested for a
 92## solution vector whose residual vector is already known, and it is the last
 93## accepted successful step.  Often this will be one of the last two calls, but
 94## not always.  If the savings by reusing intermediate results from residual
 95## calculation in Jacobian calculation are significant, the best strategy is to
 96## employ OutputFcn: After a vector is evaluated for residuals, if OutputFcn is
 97## called with that vector, then the intermediate results should be saved for
 98## future Jacobian evaluation, and should be kept until a Jacobian evaluation
 99## is requested or until outputfcn is called with a different vector, in which
100## case they should be dropped in favor of this most recent vector.  A short
101## example how this can be achieved follows:
102##
103## @example
104## function [fvec, fjac] = user_func (x, optimvalues, state)
105## persistent sav = [], sav0 = [];
106## if (nargin == 1)
107##   ## evaluation call
108##   if (nargout == 1)
109##     sav0.x = x; # mark saved vector
110##     ## calculate fvec, save results to sav0.
111##   elseif (nargout == 2)
112##     ## calculate fjac using sav.
113##   endif
114## else
115##   ## outputfcn call.
116##   if (all (x == sav0.x))
117##     sav = sav0;
118##   endif
119##   ## maybe output iteration status, etc.
120## endif
121## endfunction
122##
123## ## @dots{}
124##
125## fsolve (@@user_func, x0, optimset ("OutputFcn", @@user_func, @dots{}))
126## @end example
127## @seealso{fzero, optimset}
128## @end deftypefn
129
130## PKG_ADD: ## Discard result to avoid polluting workspace with ans at startup.
131## PKG_ADD: [~] = __all_opts__ ("fsolve");
132
133function [x, fvec, info, output, fjac] = fsolve (fcn, x0, options = struct ())
134
135  ## Get default options if requested.
136  if (nargin == 1 && ischar (fcn) && strcmp (fcn, 'defaults'))
137    x = optimset ("MaxIter", 400, "MaxFunEvals", Inf, \
138    "Jacobian", "off", "TolX", 1e-7, "TolFun", 1e-7,
139    "OutputFcn", [], "Updating", "on", "FunValCheck", "off",
140    "ComplexEqn", "off", "FinDiffType", "central",
141    "TypicalX", [], "AutoScaling", "off");
142    return;
143  endif
144
145  if (nargin < 2 || nargin > 3 || ! ismatrix (x0))
146    print_usage ();
147  endif
148
149  if (ischar (fcn))
150    fcn = str2func (fcn, "global");
151  elseif (iscell (fcn))
152    fcn = @(x) make_fcn_jac (x, fcn{1}, fcn{2});
153  endif
154
155  xsiz = size (x0);
156  n = numel (x0);
157
158  has_jac = strcmpi (optimget (options, "Jacobian", "off"), "on");
159  cdif = strcmpi (optimget (options, "FinDiffType", "central"), "central");
160  maxiter = optimget (options, "MaxIter", 400);
161  maxfev = optimget (options, "MaxFunEvals", Inf);
162  outfcn = optimget (options, "OutputFcn");
163  updating = strcmpi (optimget (options, "Updating", "on"), "on");
164  complexeqn = strcmpi (optimget (options, "ComplexEqn", "off"), "on");
165
166  ## Get scaling matrix using the TypicalX option. If set to "auto", the
167  ## scaling matrix is estimated using the Jacobian.
168  typicalx = optimget (options, "TypicalX");
169  if (isempty (typicalx))
170    typicalx = ones (n, 1);
171  endif
172  autoscale = strcmpi (optimget (options, "AutoScaling", "off"), "on");
173  if (! autoscale)
174    dg = 1 ./ typicalx;
175  endif
176
177  funvalchk = strcmpi (optimget (options, "FunValCheck", "off"), "on");
178
179  if (funvalchk)
180    ## Replace fcn with a guarded version.
181    fcn = @(x) guarded_eval (fcn, x, complexeqn);
182  endif
183
184  ## These defaults are rather stringent. I think that normally, user
185  ## prefers accuracy to performance.
186
187  macheps = eps (class (x0));
188
189  tolx = optimget (options, "TolX", 1e-7);
190  tolf = optimget (options, "TolFun", 1e-7);
191
192  factor = 1;
193
194  niter = 1;
195  nfev = 1;
196
197  x = x0(:);
198  info = 0;
199
200  ## Initial evaluation.
201  ## Handle arbitrary shapes of x and f and remember them.
202  fvec = fcn (reshape (x, xsiz));
203  fsiz = size (fvec);
204  fvec = fvec(:);
205  fn = norm (fvec);
206  m = length (fvec);
207  n = length (x);
208
209  if (! isempty (outfcn))
210    optimvalues.iter = niter;
211    optimvalues.funccount = nfev;
212    optimvalues.fval = fn;
213    optimvalues.searchdirection = zeros (n, 1);
214    state = 'init';
215    stop = outfcn (x, optimvalues, state);
216    if (stop)
217      info = -1;
218      break;
219    endif
220  endif
221
222  nsuciter = 0;
223
224  ## Outer loop.
225  while (niter < maxiter && nfev < maxfev && ! info)
226
227    ## Calculate function value and Jacobian (possibly via FD).
228    if (has_jac)
229      [fvec, fjac] = fcn (reshape (x, xsiz));
230      ## If the Jacobian is sparse, disable Broyden updating.
231      if (issparse (fjac))
232        updating = false;
233      endif
234      fvec = fvec(:);
235      nfev ++;
236    else
237      fjac = __fdjac__ (fcn, reshape (x, xsiz), fvec, typicalx, cdif);
238      nfev += (1 + cdif) * length (x);
239    endif
240
241    ## For square and overdetermined systems, we update a QR
242    ## factorization of the Jacobian to avoid solving a full system in each
243    ## step. In this case, we pass a triangular matrix to __dogleg__.
244    useqr = updating && m >= n && n > 10;
245
246    if (useqr)
247      ## FIXME: Currently, pivoting is mostly useless because the \ operator
248      ## cannot exploit the resulting props of the triangular factor.
249      ## Unpivoted QR is significantly faster so it doesn't seem right to pivot
250      ## just to get invariance. Original MINPACK didn't pivot either, at least
251      ## when qr updating was used.
252      [q, r] = qr (fjac, 0);
253    endif
254
255    if (autoscale)
256      ## Get column norms, use them as scaling factors.
257      jcn = norm (fjac, 'columns').';
258      if (niter == 1)
259        dg = jcn;
260        dg(dg == 0) = 1;
261      else
262        ## Rescale adaptively.
263        ## FIXME: the original minpack used the following rescaling strategy:
264        ##   dg = max (dg, jcn);
265        ## but it seems not good if we start with a bad guess yielding Jacobian
266        ## columns with large norms that later decrease, because the corresponding
267        ## variable will still be overscaled. So instead, we only give the old
268        ## scaling a small momentum, but do not honor it.
269
270        dg = max (0.1*dg, jcn);
271      endif
272    endif
273
274    if (niter == 1)
275      xn = norm (dg .* x);
276      ## FIXME: something better?
277      delta = factor * max (xn, 1);
278    endif
279
280    ## It also seems that in the case of fast (and inhomogeneously) changing
281    ## Jacobian, the Broyden updates are of little use, so maybe we could
282    ## skip them if a big disproportional change is expected. The question is,
283    ## of course, how to define the above terms :)
284
285    lastratio = 0;
286    nfail = 0;
287    nsuc = 0;
288    decfac = 0.5;
289
290    ## Inner loop.
291    while (niter <= maxiter && nfev < maxfev && ! info)
292
293      ## Get trust-region model (dogleg) minimizer.
294      if (useqr)
295        qtf = q'*fvec;
296        s = - __dogleg__ (r, qtf, dg, delta);
297        w = qtf + r * s;
298      else
299        s = - __dogleg__ (fjac, fvec, dg, delta);
300        w = fvec + fjac * s;
301      endif
302
303      sn = norm (dg .* s);
304      if (niter == 1)
305        delta = min (delta, sn);
306      endif
307
308      fvec1 = fcn (reshape (x + s, xsiz)) (:);
309      fn1 = norm (fvec1);
310      nfev ++;
311
312      if (fn1 < fn)
313        ## Scaled actual reduction.
314        actred = 1 - (fn1/fn)^2;
315      else
316        actred = -1;
317      endif
318
319      ## Scaled predicted reduction, and ratio.
320      t = norm (w);
321      if (t < fn)
322        prered = 1 - (t/fn)^2;
323        ratio = actred / prered;
324      else
325        prered = 0;
326        ratio = 0;
327      endif
328
329      ## Update delta.
330      if (ratio < min(max(0.1, 0.8*lastratio), 0.9))
331        nsuc = 0;
332        nfail ++;
333        delta *= decfac;
334        decfac ^= 1.4142;
335        if (delta <= 1e1*macheps*xn)
336          ## Trust region became uselessly small.
337          info = -3;
338          break;
339        endif
340      else
341        lastratio = ratio;
342        decfac = 0.5;
343        nfail = 0;
344        nsuc ++;
345        if (abs (1-ratio) <= 0.1)
346          delta = 1.4142*sn;
347        elseif (ratio >= 0.5 || nsuc > 1)
348          delta = max (delta, 1.4142*sn);
349        endif
350      endif
351
352      if (ratio >= 1e-4)
353        ## Successful iteration.
354        x += s;
355        xn = norm (dg .* x);
356        fvec = fvec1;
357        fn = fn1;
358        nsuciter ++;
359      endif
360
361      niter ++;
362
363      ## FIXME: should outputfcn be only called after a successful iteration?
364      if (! isempty (outfcn))
365        optimvalues.iter = niter;
366        optimvalues.funccount = nfev;
367        optimvalues.fval = fn;
368        optimvalues.searchdirection = s;
369        state = 'iter';
370        stop = outfcn (x, optimvalues, state);
371        if (stop)
372          info = -1;
373          break;
374        endif
375      endif
376
377      ## Tests for termination conditions. A mysterious place, anything
378      ## can happen if you change something here...
379
380      ## The rule of thumb (which I'm not sure M*b is quite following)
381      ## is that for a tolerance that depends on scaling, only 0 makes
382      ## sense as a default value. But 0 usually means uselessly long
383      ## iterations, so we need scaling-independent tolerances wherever
384      ## possible.
385
386      ## FIXME -- why tolf*n*xn? If abs (e) ~ abs(x) * eps is a vector
387      ## of perturbations of x, then norm (fjac*e) <= eps*n*xn, i.e. by
388      ## tolf ~ eps we demand as much accuracy as we can expect.
389      if (fn <= tolf*n*xn)
390        info = 1;
391        ## The following tests done only after successful step.
392      elseif (ratio >= 1e-4)
393        ## This one is classic. Note that we use scaled variables again,
394        ## but compare to scaled step, so nothing bad.
395        if (sn <= tolx*xn)
396          info = 2;
397          ## Again a classic one. It seems weird to use the same tolf
398          ## for two different tests, but that's what M*b manual appears
399          ## to say.
400        elseif (actred < tolf)
401          info = 3;
402        endif
403      endif
404
405      ## Criterion for recalculating Jacobian.
406      if (! updating || nfail == 2 || nsuciter < 2)
407        break;
408      endif
409
410      ## Compute the scaled Broyden update.
411      if (useqr)
412        u = (fvec1 - q*w) / sn;
413        v = dg .* ((dg .* s) / sn);
414
415        ## Update the QR factorization.
416        [q, r] = qrupdate (q, r, u, v);
417      else
418        u = (fvec1 - w);
419        v = dg .* ((dg .* s) / sn);
420
421        ## update the Jacobian
422        fjac += u * v';
423      endif
424    endwhile
425  endwhile
426
427  ## Restore original shapes.
428  x = reshape (x, xsiz);
429  fvec = reshape (fvec, fsiz);
430
431  output.iterations = niter;
432  output.successful = nsuciter;
433  output.funcCount = nfev;
434
435endfunction
436
437## An assistant function that evaluates a function handle and checks for
438## bad results.
439function [fx, jx] = guarded_eval (fun, x, complexeqn)
440  if (nargout > 1)
441    [fx, jx] = fun (x);
442  else
443    fx = fun (x);
444    jx = [];
445  endif
446
447  if (! complexeqn && ! (isreal (fx) && isreal (jx)))
448    error ("fsolve:notreal", "fsolve: non-real value encountered");
449  elseif (complexeqn && ! (isnumeric (fx) && isnumeric(jx)))
450    error ("fsolve:notnum", "fsolve: non-numeric value encountered");
451  elseif (any (isnan (fx(:))))
452    error ("fsolve:isnan", "fsolve: NaN value encountered");
453  endif
454endfunction
455
456function [fx, jx] = make_fcn_jac (x, fcn, fjac)
457  fx = fcn (x);
458  if (nargout == 2)
459    jx = fjac (x);
460  endif
461endfunction
462
463%!function retval = __f (p)
464%!  x = p(1);
465%!  y = p(2);
466%!  z = p(3);
467%!  retval = zeros (3, 1);
468%!  retval(1) = sin(x) + y**2 + log(z) - 7;
469%!  retval(2) = 3*x + 2**y -z**3 + 1;
470%!  retval(3) = x + y + z - 5;
471%!endfunction
472%!test
473%! x_opt = [ 0.599054;
474%! 2.395931;
475%! 2.005014 ];
476%! tol = 1.0e-5;
477%! [x, fval, info] = fsolve (@__f, [ 0.5; 2.0; 2.5 ]);
478%! assert (info > 0);
479%! assert (norm (x - x_opt, Inf) < tol);
480%! assert (norm (fval) < tol);
481
482%!function retval = __f (p)
483%!  x = p(1);
484%!  y = p(2);
485%!  z = p(3);
486%!  w = p(4);
487%!  retval = zeros (4, 1);
488%!  retval(1) = 3*x + 4*y + exp (z + w) - 1.007;
489%!  retval(2) = 6*x - 4*y + exp (3*z + w) - 11;
490%!  retval(3) = x^4 - 4*y^2 + 6*z - 8*w - 20;
491%!  retval(4) = x^2 + 2*y^3 + z - w - 4;
492%!endfunction
493%!test
494%! x_opt = [ -0.767297326653401, 0.590671081117440, 1.47190018629642, -1.52719341133957 ];
495%! tol = 1.0e-5;
496%! [x, fval, info] = fsolve (@__f, [-1, 1, 2, -1]);
497%! assert (info > 0);
498%! assert (norm (x - x_opt, Inf) < tol);
499%! assert (norm (fval) < tol);
500
501%!function retval = __f (p)
502%!  x = p(1);
503%!  y = p(2);
504%!  z = p(3);
505%!  retval = zeros (3, 1);
506%!  retval(1) = sin(x) + y**2 + log(z) - 7;
507%!  retval(2) = 3*x + 2**y -z**3 + 1;
508%!  retval(3) = x + y + z - 5;
509%!  retval(4) = x*x + y - z*log(z) - 1.36;
510%!endfunction
511%!test
512%! x_opt = [ 0.599054;
513%! 2.395931;
514%! 2.005014 ];
515%! tol = 1.0e-5;
516%! [x, fval, info] = fsolve (@__f, [ 0.5; 2.0; 2.5 ]);
517%! assert (info > 0);
518%! assert (norm (x - x_opt, Inf) < tol);
519%! assert (norm (fval) < tol);
520
521%!function retval = __f (p)
522%!  x = p(1);
523%!  y = p(2);
524%!  z = p(3);
525%!  retval = zeros (3, 1);
526%!  retval(1) = sin(x) + y**2 + log(z) - 7;
527%!  retval(2) = 3*x + 2**y -z**3 + 1;
528%!  retval(3) = x + y + z - 5;
529%!endfunction
530%!test
531%! x_opt = [ 0.599054;
532%! 2.395931;
533%! 2.005014 ];
534%! tol = 1.0e-5;
535%! opt = optimset ("Updating", "qrp");
536%! [x, fval, info] = fsolve (@__f, [ 0.5; 2.0; 2.5 ], opt);
537%! assert (info > 0);
538%! assert (norm (x - x_opt, Inf) < tol);
539%! assert (norm (fval) < tol);
540
541%!test
542%! b0 = 3;
543%! a0 = 0.2;
544%! x = 0:.5:5;
545%! noise = 1e-5 * sin (100*x);
546%! y = exp (-a0*x) + b0 + noise;
547%! c_opt = [a0, b0];
548%! tol = 1e-5;
549%!
550%! [c, fval, info, output] =  fsolve (@(c) (exp(-c(1)*x) + c(2) - y), [0, 0]);
551%! assert (info > 0);
552%! assert (norm (c - c_opt, Inf) < tol);
553%! assert (norm (fval) < norm (noise));
554
555
556%!function y = cfun (x)
557%!  y(1) = (1+i)*x(1)^2 - (1-i)*x(2) - 2;
558%!  y(2) = sqrt (x(1)*x(2)) - (1-2i)*x(3) + (3-4i);
559%!  y(3) = x(1) * x(2) - x(3)^2 + (3+2i);
560%!endfunction
561
562%!test
563%! x_opt = [-1+i, 1-i, 2+i];
564%! x = [i, 1, 1+i];
565%!
566%! [x, f, info] = fsolve (@cfun, x, optimset ("ComplexEqn", "on"));
567%! tol = 1e-5;
568%! assert (norm (f) < tol);
569%! assert (norm (x - x_opt, Inf) < tol);
570
571## Solve the double dogleg trust-region least-squares problem:
572## Minimize norm(r*x-b) subject to the constraint norm(d.*x) <= delta,
573## x being a convex combination of the gauss-newton and scaled gradient.
574
575## TODO: error checks
576## TODO: handle singularity, or leave it up to mldivide?
577
578function x = __dogleg__ (r, b, d, delta)
579  ## Get Gauss-Newton direction.
580  x = r \ b;
581  xn = norm (d .* x);
582  if (xn > delta)
583    ## GN is too big, get scaled gradient.
584    s = (r' * b) ./ d;
585    sn = norm (s);
586    if (sn > 0)
587      ## Normalize and rescale.
588      s = (s / sn) ./ d;
589      ## Get the line minimizer in s direction.
590      tn = norm (r*s);
591      snm = (sn / tn) / tn;
592      if (snm < delta)
593        ## Get the dogleg path minimizer.
594        bn = norm (b);
595        dxn = delta/xn; snmd = snm/delta;
596        t = (bn/sn) * (bn/xn) * snmd;
597        t -= dxn * snmd^2 - sqrt ((t-dxn)^2 + (1-dxn^2)*(1-snmd^2));
598        alpha = dxn*(1-snmd^2) / t;
599      else
600        alpha = 0;
601      endif
602    else
603      alpha = delta / xn;
604      snm = 0;
605    endif
606    ## Form the appropriate convex combination.
607    x = alpha * x + ((1-alpha) * min (snm, delta)) * s;
608  endif
609endfunction
610