Mercurial > hg > octave-lyh

comparison scripts/optimization/fsolve.m @ 8986:22c8272af34b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
improvements to fsolve & family
author Jaroslav Hajek <highegg@gmail.com>
date Tue, 17 Mar 2009 08:49:08 +0100
parents 20589a8f1a33
children 315828058e0d
comparison
equal deleted inserted replaced
8985:193804a4f82f 8986:22c8272af34b
70 ## @item 0 70 ## @item 0
71 ## Iteration limit exceeded. 71 ## Iteration limit exceeded.
72 ## @item -3 72 ## @item -3
73 ## The trust region radius became excessively small. 73 ## The trust region radius became excessively small.
74 ## @end table 74 ## @end table
75 ## 75 ##
76 ## Note: If you only have a single nonlinear equation of one variable, using 76 ## Note: If you only have a single nonlinear equation of one variable, using
77 ## @code{fzero} is usually a much better idea. 77 ## @code{fzero} is usually a much better idea.
78 ## @seealso{fzero, optimset} 78 ## @seealso{fzero, optimset}
79 ##
80 ## Note about user-supplied jacobians:
81 ## As an inherent property of the algorithm, jacobian is always requested for a
82 ## solution vector whose residual vector is already known, and it is the last
83 ## accepted successful step. Often this will be one of the last two calls, but
84 ## not always. If the savings by reusing intermediate results from residual
85 ## calculation in jacobian calculation are significant, the best strategy is to
86 ## employ OutputFcn: After a vector is evaluated for residuals, if OutputFcn is
87 ## called with that vector, then the intermediate results should be saved for
88 ## future jacobian evaluation, and should be kept until a jacobian evaluation
89 ## is requested or until outputfcn is called with a different vector, in which
90 ## case they should be dropped in favor of this most recent vector. A short
91 ## example how this can be achieved follows:
92 ## @example
93 ## function [fvec, fjac] = my_optim_func (x, optimvalues, state)
94 ## persistent sav = [], sav0 = [];
95 ## if (nargin == 1)
96 ## ## evaluation call
97 ## if (nargout == 1)
98 ## sav0.x = x; # mark saved vector
99 ## ## calculate fvec, save results to sav0.
100 ## elseif (nargout == 2)
101 ## ## calculate fjac using sav.
102 ## endif
103 ## else
104 ## ## outputfcn call.
105 ## if (all (x == sav0.x))
106 ## sav = sav0;
107 ## endif
108 ## ## maybe output iteration status etc.
109 ## endif
110 ## endfunction
111 ##
112 ## ## ....
113 ##
114 ## fsolve (@my_optim_func, x0, optimset ("OutputFcn", @my_optim_func, @dots{}))
115 ## @end example
116 ###
79 ## @end deftypefn 117 ## @end deftypefn
80 118
81 ## PKG_ADD: __all_opts__ ("fsolve"); 119 ## PKG_ADD: __all_opts__ ("fsolve");
82 120
83 function [x, fvec, info, output, fjac] = fsolve (fcn, x0, options = struct ()) 121 function [x, fvec, info, output, fjac] = fsolve (fcn, x0, options = struct ())
127 factor = 100; 165 factor = 100;
128 ## FIXME: TypicalX corresponds to user scaling (???) 166 ## FIXME: TypicalX corresponds to user scaling (???)
129 autodg = true; 167 autodg = true;
130 168
131 niter = 1; 169 niter = 1;
132 nfev = 0; 170 nfev = 1;
133 171
134 x = x0(:); 172 x = x0(:);
135 info = 0; 173 info = 0;
174
175 ## Initial evaluation.
176 fvec = fcn (reshape (x, xsiz));
177 fsiz = size (fvec);
178 fvec = fvec(:);
179 fn = norm (fvec);
180 m = length (fvec);
181 n = length (x);
182
183 if (! isempty (outfcn))
184 optimvalues.iter = niter;
185 optimvalues.funccount = nfev;
186 optimvalues.fval = fn;
187 optimvalues.searchdirection = zeros (n, 1);
188 state = 'init';
189 stop = outfcn (x, optimvalues, state);
190 if (stop)
191 info = -1;
192 break;
193 endif
194 endif
195
196 nsuciter = 0;
136 197
137 ## Outer loop. 198 ## Outer loop.
138 while (niter < maxiter && nfev < maxfev && ! info) 199 while (niter < maxiter && nfev < maxfev && ! info)
139 200
140 ## Calculate function value and Jacobian (possibly via FD). 201 ## Calculate function value and Jacobian (possibly via FD).
143 [fvec, fjac] = fcn (reshape (x, xsiz)); 204 [fvec, fjac] = fcn (reshape (x, xsiz));
144 ## If the jacobian is sparse, disable Broyden updating. 205 ## If the jacobian is sparse, disable Broyden updating.
145 if (issparse (fjac)) 206 if (issparse (fjac))
146 updating = false; 207 updating = false;
147 endif 208 endif
209 fvec = fvec(:);
148 nfev ++; 210 nfev ++;
149 else 211 else
150 [fvec, fjac] = __fdjac__ (fcn, reshape (x, xsiz)); 212 fjac = __fdjac__ (fcn, reshape (x, xsiz), fvec);
151 nfev += 1 + length (x); 213 nfev += length (x);
152 endif 214 endif
153 fsiz = size (fvec);
154 fvec = fvec(:);
155 fn = norm (fvec);
156 m = length (fvec);
157 n = length (x);
158 215
159 ## For square and overdetermined systems, we update a QR 216 ## For square and overdetermined systems, we update a QR
160 ## factorization of the jacobian to avoid solving a full system in each 217 ## factorization of the jacobian to avoid solving a full system in each
161 ## step. In this case, we pass a triangular matrix to __dogleg__. 218 ## step. In this case, we pass a triangular matrix to __dogleg__.
162 useqr = updating && m >= n && n > 10; 219 useqr = updating && m >= n && n > 10;
197 ## jacobian, the Broyden updates are of little use, so maybe we could 254 ## jacobian, the Broyden updates are of little use, so maybe we could
198 ## skip them if a big disproportional change is expected. The question is, 255 ## skip them if a big disproportional change is expected. The question is,
199 ## of course, how to define the above terms :) 256 ## of course, how to define the above terms :)
200 endif 257 endif
201 258
259 lastratio = 0;
202 nfail = 0; 260 nfail = 0;
203 nsuc = 0; 261 nsuc = 0;
204 ## Inner loop. 262 ## Inner loop.
205 while (niter <= maxiter && nfev < maxfev && ! info) 263 while (niter <= maxiter && nfev < maxfev && ! info)
206 264
247 prered = 0; 305 prered = 0;
248 ratio = 0; 306 ratio = 0;
249 endif 307 endif
250 308
251 ## Update delta. 309 ## Update delta.
252 if (ratio < 0.1) 310 if (ratio < min(max(0.1, lastratio), 0.9))
253 nsuc = 0; 311 nsuc = 0;
254 nfail ++; 312 nfail ++;
255 delta *= 0.5; 313 delta *= 0.5;
256 if (delta <= 1e1*macheps*xn) 314 if (delta <= 1e1*macheps*xn)
257 ## Trust region became uselessly small. 315 ## Trust region became uselessly small.
258 info = -3; 316 info = -3;
259 break; 317 break;
260 endif 318 endif
261 else 319 else
320 lastratio = ratio;
262 nfail = 0; 321 nfail = 0;
263 nsuc ++; 322 nsuc ++;
264 if (abs (1-ratio) <= 0.1) 323 if (abs (1-ratio) <= 0.1)
265 delta = 2*sn; 324 delta = 2*sn;
266 elseif (ratio >= 0.5 || nsuc > 1) 325 elseif (ratio >= 0.5 || nsuc > 1)
272 ## Successful iteration. 331 ## Successful iteration.
273 x += s; 332 x += s;
274 xn = norm (dg .* x); 333 xn = norm (dg .* x);
275 fvec = fvec1; 334 fvec = fvec1;
276 fn = fn1; 335 fn = fn1;
336 nsuciter ++;
277 endif 337 endif
278 338
279 niter ++; 339 niter ++;
280 340
281 ## FIXME: should outputfcn be only called after a successful iteration? 341 ## FIXME: should outputfcn be only called after a successful iteration?
319 info = 3; 379 info = 3;
320 endif 380 endif
321 endif 381 endif
322 382
323 ## Criterion for recalculating jacobian. 383 ## Criterion for recalculating jacobian.
324 if (! updating || nfail == 2) 384 if (! updating || nfail == 2 || nsuciter < 2)
325 break; 385 break;
326 endif 386 endif
327 387
328 ## Compute the scaled Broyden update. 388 ## Compute the scaled Broyden update.
329 if (useqr) 389 if (useqr)
345 ## Restore original shapes. 405 ## Restore original shapes.
346 x = reshape (x, xsiz); 406 x = reshape (x, xsiz);
347 fvec = reshape (fvec, fsiz); 407 fvec = reshape (fvec, fsiz);
348 408
349 output.iterations = niter; 409 output.iterations = niter;
410 output.successful = nsuciter;
350 output.funcCount = nfev; 411 output.funcCount = nfev;
351 412
352 endfunction 413 endfunction
353 414
354 ## An assistant function that evaluates a function handle and checks for 415 ## An assistant function that evaluates a function handle and checks for