octave-nkf: scripts/optimization/fsolve.m annotate

annotate scripts/optimization/fsolve.m @ 8648:ff61b53eb294

optimization: use PKG_ADD: comments instead of PKG_ADD file

author	John W. Eaton <jwe@octave.org>
date	Fri, 30 Jan 2009 15:42:51 -0500
parents	06f5dd901f30
children	9c092b111b1d

rev	line source
8590 c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	1 ## Copyright (C) 2008, 2009 VZLU Prague, a.s.
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	2 ##
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	3 ## This file is part of Octave.
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	4 ##
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	5 ## Octave is free software; you can redistribute it and/or modify it
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	6 ## under the terms of the GNU General Public License as published by
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	7 ## the Free Software Foundation; either version 3 of the License, or (at
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	8 ## your option) any later version.
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	9 ##
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	10 ## Octave is distributed in the hope that it will be useful, but
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	13 ## General Public License for more details.
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	14 ##
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	15 ## You should have received a copy of the GNU General Public License
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	16 ## along with Octave; see the file COPYING. If not, see
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	17 ## <http://www.gnu.org/licenses/>.
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	18 ##
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	19 ## Author: Jaroslav Hajek <highegg@gmail.com>
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	20
8466 4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	21 ## -- texinfo --
8513 352d3245d4c1 fsolve.m: doc fix John W. Eaton <jwe@octave.org> parents: 8507 diff changeset	22 ## @deftypefn{Function File} {} fsolve (@var{fcn}, @var{x0}, @var{options})
8466 4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	23 ## @deftypefnx{Function File} {[@var{x}, @var{fvec}, @var{info}, @var{output}, @var{fjac}]} = fsolve (@var{fcn}, @dots{})
8513 352d3245d4c1 fsolve.m: doc fix John W. Eaton <jwe@octave.org> parents: 8507 diff changeset	24 ## Solve a system of nonlinear equations defined by the function @var{fcn}.
8466 4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	25 ## @var{fcn} should accepts a vector (array) defining the unknown variables,
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	26 ## and return a vector of left-hand sides of the equations. Right-hand sides
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	27 ## are defined to be zeros.
8514 39867b4aca52 fsolve.m: additional doc fix John W. Eaton <jwe@octave.org> parents: 8513 diff changeset	28 ## In other words, this function attempts to determine a vector @var{x} such
39867b4aca52 fsolve.m: additional doc fix John W. Eaton <jwe@octave.org> parents: 8513 diff changeset	29 ## that @code{@var{fcn} (@var{x})} gives (approximately) all zeros.
8466 4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	30 ## @var{x0} determines a starting guess. The shape of @var{x0} is preserved
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	31 ## in all calls to @var{fcn}, but otherwise it is treated as a column vector.
8515 ec2715c76039 fzero.m, fsolve.m: additional doc fixes John W. Eaton <jwe@octave.org> parents: 8514 diff changeset	32 ## @var{options} is a structure specifying additional options.
ec2715c76039 fzero.m, fsolve.m: additional doc fixes John W. Eaton <jwe@octave.org> parents: 8514 diff changeset	33 ## Currently, @code{fsolve} recognizes these options:
ec2715c76039 fzero.m, fsolve.m: additional doc fixes John W. Eaton <jwe@octave.org> parents: 8514 diff changeset	34 ## @code{"FunValCheck"}, @code{"OutputFcn"}, @code{"TolX"},
8590 c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	35 ## @code{"TolFun"}, @code{"MaxIter"}, @code{"MaxFunEvals"},
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	36 ## @code{"Jacobian"} and @code{"Updating"}.
8466 4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	37 ##
8513 352d3245d4c1 fsolve.m: doc fix John W. Eaton <jwe@octave.org> parents: 8507 diff changeset	38 ## If @code{"Jacobian"} is @code{"on"}, it specifies that @var{fcn},
352d3245d4c1 fsolve.m: doc fix John W. Eaton <jwe@octave.org> parents: 8507 diff changeset	39 ## called with 2 output arguments, also returns the Jacobian matrix
352d3245d4c1 fsolve.m: doc fix John W. Eaton <jwe@octave.org> parents: 8507 diff changeset	40 ## of right-hand sides at the requested point. @code{"TolX"} specifies
352d3245d4c1 fsolve.m: doc fix John W. Eaton <jwe@octave.org> parents: 8507 diff changeset	41 ## the termination tolerance in the unknown variables, while
8596 8833c0b18eb2 enable default settings queries in optim funcs Jaroslav Hajek <highegg@gmail.com> parents: 8592 diff changeset	42 ## @code{"TolFun"} is a tolerance for equations. Default is @code{1e-7}
8833c0b18eb2 enable default settings queries in optim funcs Jaroslav Hajek <highegg@gmail.com> parents: 8592 diff changeset	43 ## for both @code{"TolX"} and @code{"TolFun"}.
8833c0b18eb2 enable default settings queries in optim funcs Jaroslav Hajek <highegg@gmail.com> parents: 8592 diff changeset	44 ## If @code{"Updating"} is "on", the function will attempt to use Broyden
8590 c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	45 ## updates to update the Jacobian, in order to reduce the amount of jacobian
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	46 ## calculations.
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	47 ## If your user function always calculates the Jacobian (regardless of number
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	48 ## of output arguments), this option provides no advantage and should be set to
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	49 ## false.
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	50 ##
8466 4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	51 ## For description of the other options, see @code{optimset}.
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	52 ##
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	53 ## On return, @var{fval} contains the value of the function @var{fcn}
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	54 ## evaluated at @var{x}, and @var{info} may be one of the following values:
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	55 ##
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	56 ## @table @asis
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	57 ## @item 1
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	58 ## Converged to a solution point. Relative residual error is less than specified
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	59 ## by TolFun.
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	60 ## @item 2
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	61 ## Last relative step size was less that TolX.
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	62 ## @item 3
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	63 ## Last relative decrease in residual was less than TolF.
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	64 ## @item 0
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	65 ## Iteration limit exceeded.
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	66 ## @item -3
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	67 ## The trust region radius became excessively small.
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	68 ## @end table
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	69 ##
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	70 ## Note: If you only have a single nonlinear equation of one variable, using
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	71 ## @code{fzero} is usually a much better idea.
8513 352d3245d4c1 fsolve.m: doc fix John W. Eaton <jwe@octave.org> parents: 8507 diff changeset	72 ## @seealso{fzero, optimset}
8466 4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	73 ## @end deftypefn
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	74
8648 ff61b53eb294 optimization: use PKG_ADD: comments instead of PKG_ADD file John W. Eaton <jwe@octave.org> parents: 8647 diff changeset	75 ## PKG_ADD: __all_opts__ ("fsolve");
ff61b53eb294 optimization: use PKG_ADD: comments instead of PKG_ADD file John W. Eaton <jwe@octave.org> parents: 8647 diff changeset	76
8596 8833c0b18eb2 enable default settings queries in optim funcs Jaroslav Hajek <highegg@gmail.com> parents: 8592 diff changeset	77 function [x, fvec, info, output, fjac] = fsolve (fcn, x0, options = struct ())
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	78
8596 8833c0b18eb2 enable default settings queries in optim funcs Jaroslav Hajek <highegg@gmail.com> parents: 8592 diff changeset	79 ## Get default options if requested.
8833c0b18eb2 enable default settings queries in optim funcs Jaroslav Hajek <highegg@gmail.com> parents: 8592 diff changeset	80 if (nargin == 1 && ischar (fcn) && strcmp (fcn, 'defaults'))
8833c0b18eb2 enable default settings queries in optim funcs Jaroslav Hajek <highegg@gmail.com> parents: 8592 diff changeset	81 x = optimset ("MaxIter", 400, "MaxFunEvals", Inf, \
8647 06f5dd901f30 implement registering of optimization options Jaroslav Hajek <highegg@gmail.com> parents: 8616 diff changeset	82 "Jacobian", "off", "TolX", 1.5e-8, "TolFun", 1.5e-8,
8596 8833c0b18eb2 enable default settings queries in optim funcs Jaroslav Hajek <highegg@gmail.com> parents: 8592 diff changeset	83 "OutputFcn", [], "Updating", "on", "FunValCheck", "off");
8833c0b18eb2 enable default settings queries in optim funcs Jaroslav Hajek <highegg@gmail.com> parents: 8592 diff changeset	84 return;
8833c0b18eb2 enable default settings queries in optim funcs Jaroslav Hajek <highegg@gmail.com> parents: 8592 diff changeset	85 endif
8833c0b18eb2 enable default settings queries in optim funcs Jaroslav Hajek <highegg@gmail.com> parents: 8592 diff changeset	86
8833c0b18eb2 enable default settings queries in optim funcs Jaroslav Hajek <highegg@gmail.com> parents: 8592 diff changeset	87 if (nargin < 2 \|\| nargin > 3 \|\| ! ismatrix (x0))
8833c0b18eb2 enable default settings queries in optim funcs Jaroslav Hajek <highegg@gmail.com> parents: 8592 diff changeset	88 print_usage ();
8833c0b18eb2 enable default settings queries in optim funcs Jaroslav Hajek <highegg@gmail.com> parents: 8592 diff changeset	89 endif
8833c0b18eb2 enable default settings queries in optim funcs Jaroslav Hajek <highegg@gmail.com> parents: 8592 diff changeset	90
8833c0b18eb2 enable default settings queries in optim funcs Jaroslav Hajek <highegg@gmail.com> parents: 8592 diff changeset	91 if (ischar (fcn))
8833c0b18eb2 enable default settings queries in optim funcs Jaroslav Hajek <highegg@gmail.com> parents: 8592 diff changeset	92 fcn = str2func (fcn);
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	93 endif
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	94
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	95 xsiz = size (x0);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	96 n = numel (x0);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	97
8467 77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options John W. Eaton <jwe@octave.org> parents: 8466 diff changeset	98 has_jac = strcmpi (optimget (options, "Jacobian", "off"), "on");
8596 8833c0b18eb2 enable default settings queries in optim funcs Jaroslav Hajek <highegg@gmail.com> parents: 8592 diff changeset	99 maxiter = optimget (options, "MaxIter", 400);
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	100 maxfev = optimget (options, "MaxFunEvals", Inf);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	101 outfcn = optimget (options, "OutputFcn");
8596 8833c0b18eb2 enable default settings queries in optim funcs Jaroslav Hajek <highegg@gmail.com> parents: 8592 diff changeset	102 updating = strcmpi (optimget (options, "Updating", "on"), "on");
8590 c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	103
8467 77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options John W. Eaton <jwe@octave.org> parents: 8466 diff changeset	104 funvalchk = strcmpi (optimget (options, "FunValCheck", "off"), "on");
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	105
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	106 if (funvalchk)
8604 43f831758d42 fix typo in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8596 diff changeset	107 ## Replace fcn with a guarded version.
43f831758d42 fix typo in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8596 diff changeset	108 fcn = @(x) guarded_eval (fcn, x);
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	109 endif
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	110
8466 4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	111 ## These defaults are rather stringent. I think that normally, user
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	112 ## prefers accuracy to performance.
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	113
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	114 macheps = eps (class (x0));
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	115
8613 38482007c834 relax scaling in fsolve and use class-dependent default tolerances Jaroslav Hajek <highegg@gmail.com> parents: 8604 diff changeset	116 tolx = optimget (options, "TolX", sqrt (macheps));
38482007c834 relax scaling in fsolve and use class-dependent default tolerances Jaroslav Hajek <highegg@gmail.com> parents: 8604 diff changeset	117 tolf = optimget (options, "TolFun", sqrt (macheps));
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	118
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	119 factor = 100;
8466 4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	120 ## FIXME: TypicalX corresponds to user scaling (???)
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	121 autodg = true;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	122
8467 77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options John W. Eaton <jwe@octave.org> parents: 8466 diff changeset	123 niter = 1;
77b8d4aa2743 fsolve.m, fzero.m: style fixes; use strcmpi to compare options John W. Eaton <jwe@octave.org> parents: 8466 diff changeset	124 nfev = 0;
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	125
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	126 x = x0(:);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	127 info = 0;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	128
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 diff changeset	129 ## Outer loop.
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	130 while (niter < maxiter && nfev < maxfev && ! info)
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	131
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 diff changeset	132 ## Calculate function value and Jacobian (possibly via FD).
cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 diff changeset	133 ## Handle arbitrary shapes of x and f and remember them.
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	134 if (has_jac)
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	135 [fvec, fjac] = fcn (reshape (x, xsiz));
8592 dacfd030633a handle sparse jacobians in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8590 diff changeset	136 ## If the jacobian is sparse, disable Broyden updating.
dacfd030633a handle sparse jacobians in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8590 diff changeset	137 if (issparse (fjac))
dacfd030633a handle sparse jacobians in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8590 diff changeset	138 updating = false;
dacfd030633a handle sparse jacobians in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8590 diff changeset	139 endif
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	140 nfev ++;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	141 else
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	142 [fvec, fjac] = __fdjac__ (fcn, reshape (x, xsiz));
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	143 nfev += 1 + length (x);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	144 endif
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	145 fsiz = size (fvec);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	146 fvec = fvec(:);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	147 fn = norm (fvec);
8383 a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	148 m = length (fvec);
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	149 n = length (x);
8590 c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	150
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	151 ## For square and overdetermined systems, we update a (pivoted) QR
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	152 ## factorization of the jacobian to avoid solving a full system in each
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	153 ## step. In this case, we pass a triangular matrix to __dogleg__.
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	154 ## Pivoted QR is used for slightly better robustness and invariance
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	155 ## w.r.t. permutations of variables.
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	156 useqr = updating && m >= n && n > 10;
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	157
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	158 if (useqr)
8616 3d75d717cbe0 do not pivot by default in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8615 diff changeset	159 ## FIXME: Currently, pivoting is mostly useless because the \ operator
3d75d717cbe0 do not pivot by default in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8615 diff changeset	160 ## cannot exploit the resulting props of the triangular factor.
3d75d717cbe0 do not pivot by default in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8615 diff changeset	161 ## Unpivoted QR is significantly faster so it doesn't seem right to pivot
3d75d717cbe0 do not pivot by default in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8615 diff changeset	162 ## just to get invariance. Original MINPACK didn't pivot either, at least
3d75d717cbe0 do not pivot by default in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8615 diff changeset	163 ## when qr updating was used.
3d75d717cbe0 do not pivot by default in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8615 diff changeset	164 [q, r] = qr (fjac, 0);
8383 a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	165 endif
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	166
8590 c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	167 ## Get column norms, use them as scaling factors.
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	168 jcn = norm (fjac, 'columns').';
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	169 if (niter == 1)
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	170 if (autodg)
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	171 dg = jcn;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	172 dg(dg == 0) = 1;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	173 endif
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	174 xn = norm (dg .* x);
8590 c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	175 ## FIXME: something better?
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	176 delta = max (factor * xn, 1);
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	177 endif
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	178
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 diff changeset	179 ## Rescale if necessary.
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	180 if (autodg)
8613 38482007c834 relax scaling in fsolve and use class-dependent default tolerances Jaroslav Hajek <highegg@gmail.com> parents: 8604 diff changeset	181 ## FIXME: the original minpack used the following rescaling strategy:
38482007c834 relax scaling in fsolve and use class-dependent default tolerances Jaroslav Hajek <highegg@gmail.com> parents: 8604 diff changeset	182 ## dg = max (dg, jcn);
38482007c834 relax scaling in fsolve and use class-dependent default tolerances Jaroslav Hajek <highegg@gmail.com> parents: 8604 diff changeset	183 ## but it seems not good if we start with a bad guess yielding jacobian
38482007c834 relax scaling in fsolve and use class-dependent default tolerances Jaroslav Hajek <highegg@gmail.com> parents: 8604 diff changeset	184 ## columns with large norms that later decrease, because the corresponding
38482007c834 relax scaling in fsolve and use class-dependent default tolerances Jaroslav Hajek <highegg@gmail.com> parents: 8604 diff changeset	185 ## variable will still be overscaled. So instead, we only give the old
38482007c834 relax scaling in fsolve and use class-dependent default tolerances Jaroslav Hajek <highegg@gmail.com> parents: 8604 diff changeset	186 ## scaling a small momentum, but do not honor it.
38482007c834 relax scaling in fsolve and use class-dependent default tolerances Jaroslav Hajek <highegg@gmail.com> parents: 8604 diff changeset	187
38482007c834 relax scaling in fsolve and use class-dependent default tolerances Jaroslav Hajek <highegg@gmail.com> parents: 8604 diff changeset	188 dg = max (0.1*dg, jcn);
38482007c834 relax scaling in fsolve and use class-dependent default tolerances Jaroslav Hajek <highegg@gmail.com> parents: 8604 diff changeset	189
38482007c834 relax scaling in fsolve and use class-dependent default tolerances Jaroslav Hajek <highegg@gmail.com> parents: 8604 diff changeset	190 ## It also seems that in the case of fast (and inhomogeneously) changing
38482007c834 relax scaling in fsolve and use class-dependent default tolerances Jaroslav Hajek <highegg@gmail.com> parents: 8604 diff changeset	191 ## jacobian, the Broyden updates are of little use, so maybe we could
38482007c834 relax scaling in fsolve and use class-dependent default tolerances Jaroslav Hajek <highegg@gmail.com> parents: 8604 diff changeset	192 ## skip them if a big disproportional change is expected. The question is,
38482007c834 relax scaling in fsolve and use class-dependent default tolerances Jaroslav Hajek <highegg@gmail.com> parents: 8604 diff changeset	193 ## of course, how to define the above terms :)
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	194 endif
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	195
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	196 nfail = 0;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	197 nsuc = 0;
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 diff changeset	198 ## Inner loop.
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	199 while (niter <= maxiter && nfev < maxfev && ! info)
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	200
8590 c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	201 if (useqr)
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	202 tr_mat = r;
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	203 tr_vec = q'*fvec;
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	204 else
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	205 tr_mat = fjac;
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	206 tr_vec = fvec;
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	207 endif
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	208
8590 c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	209 ## Get trust-region model (dogleg) minimizer.
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	210 if (useqr)
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	211 qtf = q'*fvec;
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	212 s = - __dogleg__ (r, qtf, dg, delta);
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	213 w = qtf + r * s;
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	214 s = p * s;
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	215 else
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	216 s = - __dogleg__ (fjac, fvec, dg, delta);
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	217 w = fvec + fjac * s;
8395 88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	218 endif
8590 c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	219
8395 88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	220 sn = norm (dg .* s);
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	221 if (niter == 1)
8395 88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	222 delta = min (delta, sn);
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	223 endif
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	224
8395 88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	225 fvec1 = fcn (reshape (x + s, xsiz)) (:);
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	226 fn1 = norm (fvec1);
8590 c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	227 nfev ++;
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	228
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	229 if (fn1 < fn)
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 diff changeset	230 ## Scaled actual reduction.
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	231 actred = 1 - (fn1/fn)^2;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	232 else
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	233 actred = -1;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	234 endif
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	235
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 diff changeset	236 ## Scaled predicted reduction, and ratio.
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	237 t = norm (w);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	238 if (t < fn)
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	239 prered = 1 - (t/fn)^2;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	240 ratio = actred / prered;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	241 else
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	242 prered = 0;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	243 ratio = 0;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	244 endif
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	245
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 diff changeset	246 ## Update delta.
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	247 if (ratio < 0.1)
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	248 nsuc = 0;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	249 nfail ++;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	250 delta *= 0.5;
8590 c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	251 if (delta <= 1e1machepsxn)
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 diff changeset	252 ## Trust region became uselessly small.
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	253 info = -3;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	254 break;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	255 endif
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	256 else
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	257 nfail = 0;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	258 nsuc ++;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	259 if (abs (1-ratio) <= 0.1)
8395 88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	260 delta = 2*sn;
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	261 elseif (ratio >= 0.5 \|\| nsuc > 1)
8395 88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	262 delta = max (delta, 2*sn);
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	263 endif
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	264 endif
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	265
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	266 if (ratio >= 1e-4)
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 diff changeset	267 ## Successful iteration.
8395 88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	268 x += s;
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	269 xn = norm (dg .* x);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	270 fvec = fvec1;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	271 fn = fn1;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	272 niter ++;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	273 endif
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	274
8590 c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	275 ## FIXME: should outputfcn be only called after a successful iteration?
8615 7baacb6a8a65 fix typo in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8613 diff changeset	276 if (! isempty (outfcn))
8590 c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	277 optimvalues.iter = niter;
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	278 optimvalues.funccount = nfev;
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	279 optimvalues.fval = fn;
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	280 optimvalues.searchdirection = s;
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	281 state = 'iter';
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	282 stop = outfcn (x, optimvalues, state);
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	283 if (stop)
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	284 info = -1;
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	285 break;
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	286 endif
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	287 endif
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	288
8466 4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	289 ## Tests for termination conditions. A mysterious place, anything
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	290 ## can happen if you change something here...
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	291
8466 4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	292 ## The rule of thumb (which I'm not sure M*b is quite following)
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	293 ## is that for a tolerance that depends on scaling, only 0 makes
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	294 ## sense as a default value. But 0 usually means uselessly long
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	295 ## iterations, so we need scaling-independent tolerances wherever
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	296 ## possible.
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	297
8466 4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	298 ## FIXME -- why tolfnxn? If abs (e) ~ abs(x) * eps is a vector
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	299 ## of perturbations of x, then norm (fjace) <= epsn*xn, i.e. by
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	300 ## tolf ~ eps we demand as much accuracy as we can expect.
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	301 if (fn <= tolfnxn)
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	302 info = 1;
8466 4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	303 ## The following tests done only after successful step.
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	304 elseif (actred > 0)
8466 4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	305 ## This one is classic. Note that we use scaled variables again,
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	306 ## but compare to scaled step, so nothing bad.
8395 88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	307 if (sn <= tolx*xn)
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	308 info = 2;
8466 4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	309 ## Again a classic one. It seems weird to use the same tolf
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	310 ## for two different tests, but that's what M*b manual appears
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	311 ## to say.
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	312 elseif (actred < tolf)
8383 a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	313 info = 3;
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	314 endif
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	315 endif
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	316
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 diff changeset	317 ## Criterion for recalculating jacobian.
8590 c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	318 if (! updating \|\| nfail == 2)
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	319 break;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	320 endif
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	321
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 diff changeset	322 ## Compute the scaled Broyden update.
8590 c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	323 if (useqr)
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	324 u = (fvec1 - q*w) / sn;
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	325 v = dg .* ((dg .* s) / sn);
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	326
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	327 ## Update the QR factorization.
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	328 [q, r] = qrupdate (q, r, u, v);
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	329 else
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	330 u = (fvec1 - w);
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	331 v = dg .* ((dg .* s) / sn);
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	332
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	333 ## update the jacobian
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	334 fjac += u * v';
8395 88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	335 endif
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	336 endwhile
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	337 endwhile
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	338
8507 cadc73247d65 style fixes John W. Eaton <jwe@octave.org> parents: 8468 diff changeset	339 ## Restore original shapes.
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	340 x = reshape (x, xsiz);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	341 fvec = reshape (fvec, fsiz);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	342
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	343 output.iterations = niter;
8590 c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	344 output.funcCount = nfev;
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	345
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	346 endfunction
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	347
8466 4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	348 ## An assistant function that evaluates a function handle and checks for
4d008d9f0ccf fsolve.m: style fixes John W. Eaton <jwe@octave.org> parents: 8395 diff changeset	349 ## bad results.
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	350 function fx = guarded_eval (fun, x)
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	351 fx = fun (x);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	352 if (! all (isreal (fx)))
8468 866492035ecf fsolve.m, fzero.m: undo part of previous change John W. Eaton <jwe@octave.org> parents: 8467 diff changeset	353 error ("fsolve:notreal", "fsolve: non-real value encountered");
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	354 elseif (any (isnan (fx)))
8468 866492035ecf fsolve.m, fzero.m: undo part of previous change John W. Eaton <jwe@octave.org> parents: 8467 diff changeset	355 error ("fsolve:isnan", "fsolve: NaN value encountered");
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	356 endif
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	357 endfunction
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	358
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	359 %!function retval = f (p)
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	360 %! x = p(1);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	361 %! y = p(2);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	362 %! z = p(3);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	363 %! retval = zeros (3, 1);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	364 %! retval(1) = sin(x) + y**2 + log(z) - 7;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	365 %! retval(2) = 3x + 2y -z*3 + 1;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	366 %! retval(3) = x + y + z - 5;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	367 %!test
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	368 %! x_opt = [ 0.599054;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	369 %! 2.395931;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	370 %! 2.005014 ];
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	371 %! tol = 1.0e-5;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	372 %! [x, fval, info] = fsolve (@f, [ 0.5; 2.0; 2.5 ]);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	373 %! assert (info > 0);
8383 a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	374 %! assert (norm (x - x_opt, Inf) < tol);
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	375 %! assert (norm (fval) < tol);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	376
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	377 %!function retval = f (p)
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	378 %! x = p(1);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	379 %! y = p(2);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	380 %! z = p(3);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	381 %! w = p(4);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	382 %! retval = zeros (4, 1);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	383 %! retval(1) = 3x + 4y + exp (z + w) - 1.007;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	384 %! retval(2) = 6x - 4y + exp (3*z + w) - 11;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	385 %! retval(3) = x^4 - 4y^2 + 6z - 8*w - 20;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	386 %! retval(4) = x^2 + 2*y^3 + z - w - 4;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	387 %!test
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	388 %! x_opt = [ -0.767297326653401, 0.590671081117440, 1.47190018629642, -1.52719341133957 ];
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	389 %! tol = 1.0e-5;
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	390 %! [x, fval, info] = fsolve (@f, [-1, 1, 2, -1]);
43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	391 %! assert (info > 0);
8383 a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	392 %! assert (norm (x - x_opt, Inf) < tol);
8306 43795cf108d0 initial implementation of fsolve Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	393 %! assert (norm (fval) < tol);
8383 a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	394
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	395 %!function retval = f (p)
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	396 %! x = p(1);
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	397 %! y = p(2);
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	398 %! z = p(3);
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	399 %! retval = zeros (3, 1);
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	400 %! retval(1) = sin(x) + y**2 + log(z) - 7;
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	401 %! retval(2) = 3x + 2y -z*3 + 1;
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	402 %! retval(3) = x + y + z - 5;
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	403 %! retval(4) = xx + y - zlog(z) - 1.36;
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	404 %!test
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	405 %! x_opt = [ 0.599054;
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	406 %! 2.395931;
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	407 %! 2.005014 ];
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	408 %! tol = 1.0e-5;
8590 c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	409 %! [x, fval, info] = fsolve (@f, [ 0.5; 2.0; 2.5 ]);
8383 a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	410 %! assert (info > 0);
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	411 %! assert (norm (x - x_opt, Inf) < tol);
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	412 %! assert (norm (fval) < tol);
a762d9daa700 allow overdetermined systems in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8306 diff changeset	413
8395 88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	414 %!function retval = f (p)
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	415 %! x = p(1);
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	416 %! y = p(2);
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	417 %! z = p(3);
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	418 %! retval = zeros (3, 1);
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	419 %! retval(1) = sin(x) + y**2 + log(z) - 7;
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	420 %! retval(2) = 3x + 2y -z*3 + 1;
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	421 %! retval(3) = x + y + z - 5;
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	422 %!test
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	423 %! x_opt = [ 0.599054;
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	424 %! 2.395931;
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	425 %! 2.005014 ];
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	426 %! tol = 1.0e-5;
8590 c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	427 %! opt = optimset ("Updating", "qrp");
8395 88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	428 %! [x, fval, info] = fsolve (@f, [ 0.5; 2.0; 2.5 ], opt);
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	429 %! assert (info > 0);
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	430 %! assert (norm (x - x_opt, Inf) < tol);
88fd356b0d95 optionally allow pivoting in fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8383 diff changeset	431 %! assert (norm (fval) < tol);
8590 c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	432
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	433 %!test
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	434 %! b0 = 3;
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	435 %! a0 = 0.2;
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	436 %! x = 0:.5:5;
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	437 %! noise = 1e-5 * sin (100*x);
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	438 %! y = exp (-a0*x) + b0 + noise;
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	439 %! c_opt = [a0, b0];
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	440 %! tol = 1e-5;
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	441 %!
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	442 %! [c, fval, info, output] = fsolve (@(c) (exp(-c(1)*x) + c(2) - y), [0, 0]);
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	443 %! assert (info > 0);
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	444 %! assert (norm (c - c_opt, Inf) < tol);
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	445 %! assert (norm (fval) < norm (noise));
c136d313206a improvements to fsolve Jaroslav Hajek <highegg@gmail.com> parents: 8562 diff changeset	446

Mercurial > hg > octave-nkf

annotate scripts/optimization/fsolve.m @ 8648:ff61b53eb294