Mercurial > hg > octave-nkf

comparison scripts/optimization/fsolve.m @ 8466:4d008d9f0ccf

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
fsolve.m: style fixes
author John W. Eaton <jwe@octave.org>
date Mon, 12 Jan 2009 13:04:13 -0500
parents 88fd356b0d95
children 77b8d4aa2743
comparison
equal deleted inserted replaced
8465:cddf85b1524a 8466:4d008d9f0ccf
16 ## along with Octave; see the file COPYING. If not, see 16 ## along with Octave; see the file COPYING. If not, see
17 ## <http://www.gnu.org/licenses/>. 17 ## <http://www.gnu.org/licenses/>.
18 ## 18 ##
19 ## Author: Jaroslav Hajek <highegg@gmail.com> 19 ## Author: Jaroslav Hajek <highegg@gmail.com>
20 20
21 # -*- texinfo -*- 21 ## -*- texinfo -*-
22 # @deftypefn{Function File} {} fsolve(@var{fcn}, @var{x0}, @var{options}) 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{}) 23 ## @deftypefnx{Function File} {[@var{x}, @var{fvec}, @var{info}, @var{output}, @var{fjac}]} = fsolve (@var{fcn}, @dots{})
24 # Solves a system of nonlinear equations defined by the function @var{fcn}. 24 ## Solves a system of nonlinear equations defined by the function @var{fcn}.
25 # @var{fcn} should accepts a vector (array) defining the unknown variables, 25 ## @var{fcn} should accepts a vector (array) defining the unknown variables,
26 # and return a vector of left-hand sides of the equations. Right-hand sides 26 ## and return a vector of left-hand sides of the equations. Right-hand sides
27 # are defined to be zeros. 27 ## are defined to be zeros.
28 # In other words, this function attempts to determine a vector @var{X} such 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. 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 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. 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. Currently, fsolve 32 ## @var{options} is a structure specifying additional options. Currently, fsolve
33 # recognizes these options: FunValCheck, OutputFcn, TolX, TolFun, MaxIter, 33 ## recognizes these options: FunValCheck, OutputFcn, TolX, TolFun, MaxIter,
34 # MaxFunEvals and Jacobian. 34 ## MaxFunEvals and Jacobian.
35 # 35 ##
36 # If Jacobian is 'on', it specifies that @var{fcn}, called with 2 output arguments, 36 ## If Jacobian is 'on', it specifies that @var{fcn}, called with 2 output arguments,
37 # also returns the Jacobian matrix of right-hand sides at the requested point. 37 ## also returns the Jacobian matrix of right-hand sides at the requested point.
38 # TolX specifies the termination tolerance in the unknown variables, while 38 ## TolX specifies the termination tolerance in the unknown variables, while
39 # TolFun is a tolerance for equations. Default is @code{1e1*eps} 39 ## TolFun is a tolerance for equations. Default is @code{1e1*eps}
40 # for TolX and @code{1e2*eps} for TolFun. 40 ## for TolX and @code{1e2*eps} for TolFun.
41 # For description of the other options, see @code{optimset}. 41 ## For description of the other options, see @code{optimset}.
42 # 42 ##
43 # On return, @var{fval} contains the value of the function @var{fcn} 43 ## On return, @var{fval} contains the value of the function @var{fcn}
44 # evaluated at @var{x}, and @var{info} may be one of the following values: 44 ## evaluated at @var{x}, and @var{info} may be one of the following values:
45 # 45 ##
46 # @table @asis 46 ## @table @asis
47 # @item 1 47 ## @item 1
48 # Converged to a solution point. Relative residual error is less than specified 48 ## Converged to a solution point. Relative residual error is less than specified
49 # by TolFun. 49 ## by TolFun.
50 # @item 2 50 ## @item 2
51 # Last relative step size was less that TolX. 51 ## Last relative step size was less that TolX.
52 # @item 3 52 ## @item 3
53 # Last relative decrease in residual was less than TolF. 53 ## Last relative decrease in residual was less than TolF.
54 # @item 0 54 ## @item 0
55 # Iteration limit exceeded. 55 ## Iteration limit exceeded.
56 # @item -3 56 ## @item -3
57 # The trust region radius became excessively small. 57 ## The trust region radius became excessively small.
58 # @end table 58 ## @end table
59 # 59 ##
60 # Note: If you only have a single nonlinear equation of one variable, using 60 ## Note: If you only have a single nonlinear equation of one variable, using
61 # @code{fzero} is usually a much better idea. 61 ## @code{fzero} is usually a much better idea.
62 # @seealso{fzero,optimset} 62 ## @seealso{fzero,optimset}
63 # @end deftypefn 63 ## @end deftypefn
64 64
65 function [x, fvec, info, output, fjac] = fsolve (fcn, x0, options) 65 function [x, fvec, info, output, fjac] = fsolve (fcn, x0, options)
66 66
67 if (nargin < 3) 67 if (nargin < 3)
68 options = struct (); 68 options = struct ();
77 outfcn = optimget (options, "OutputFcn"); 77 outfcn = optimget (options, "OutputFcn");
78 pivoting = optimget (options, "Pivoting", false); 78 pivoting = optimget (options, "Pivoting", false);
79 funvalchk = strcmp (optimget (options, "FunValCheck", "off"), "on"); 79 funvalchk = strcmp (optimget (options, "FunValCheck", "off"), "on");
80 80
81 if (funvalchk) 81 if (funvalchk)
82 # replace fun with a guarded version 82 ## replace fun with a guarded version
83 fun = @(x) guarded_eval (fun, x); 83 fun = @(x) guarded_eval (fun, x);
84 endif 84 endif
85 85
86 # These defaults are rather stringent. I think that normally, user prefers 86 ## These defaults are rather stringent. I think that normally, user
87 # accuracy to performance. 87 ## prefers accuracy to performance.
88 88
89 macheps = eps (class (x0)); 89 macheps = eps (class (x0));
90 90
91 tolx = optimget (options, "TolX", 1e1*macheps); 91 tolx = optimget (options, "TolX", 1e1*macheps);
92 tolf = optimget (options, "TolFun",1e2*macheps); 92 tolf = optimget (options, "TolFun",1e2*macheps);
93 93
94 factor = 100; 94 factor = 100;
95 # FIXME: TypicalX corresponds to user scaling (???) 95 ## FIXME: TypicalX corresponds to user scaling (???)
96 autodg = true; 96 autodg = true;
97 97
98 niter = 1; nfev = 0; 98 niter = 1; nfev = 0;
99 99
100 x = x0(:); 100 x = x0(:);
101 info = 0; 101 info = 0;
102 102
103 # outer loop 103 ## outer loop
104 while (niter < maxiter && nfev < maxfev && ! info) 104 while (niter < maxiter && nfev < maxfev && ! info)
105 105
106 # calc func value and jacobian (possibly via FD) 106 ## calc func value and jacobian (possibly via FD)
107 # handle arbitrary shapes of x and f and remember them 107 ## handle arbitrary shapes of x and f and remember them
108 if (has_jac) 108 if (has_jac)
109 [fvec, fjac] = fcn (reshape (x, xsiz)); 109 [fvec, fjac] = fcn (reshape (x, xsiz));
110 nfev ++; 110 nfev ++;
111 else 111 else
112 [fvec, fjac] = __fdjac__ (fcn, reshape (x, xsiz)); 112 [fvec, fjac] = __fdjac__ (fcn, reshape (x, xsiz));
116 fvec = fvec(:); 116 fvec = fvec(:);
117 fn = norm (fvec); 117 fn = norm (fvec);
118 m = length (fvec); 118 m = length (fvec);
119 n = length (x); 119 n = length (x);
120 if (m < n) 120 if (m < n)
121 error ("fsolve:under", "cannot solve underdetermined systems"); 121 error ("fsolve: cannot solve underdetermined systems");
122 elseif (m > n && niter == 1) 122 elseif (m > n && niter == 1)
123 if (isempty (optimget (options, "TolFun"))) 123 if (isempty (optimget (options, "TolFun")))
124 warning (["An overdetermined system cannot usually be solved exactly. ", ... 124 warning ("an overdetermined system cannot usually be solved exactly; consider specifying the TolFun option");
125 "Consider specifying the TolFun option."]);
126 endif 125 endif
127 endif 126 endif
128 127
129 # get QR factorization of the jacobian 128 ## get QR factorization of the jacobian
130 if (pivoting) 129 if (pivoting)
131 [q, r, p] = qr (fjac); 130 [q, r, p] = qr (fjac);
132 else 131 else
133 [q, r] = qr (fjac); 132 [q, r] = qr (fjac);
134 endif 133 endif
135 134
136 # get column norms, use them as scaling factor 135 ## get column norms, use them as scaling factor
137 jcn = norm (fjac, 'columns').'; 136 jcn = norm (fjac, 'columns').';
138 if (niter == 1) 137 if (niter == 1)
139 if (autodg) 138 if (autodg)
140 dg = jcn; 139 dg = jcn;
141 dg(dg == 0) = 1; 140 dg(dg == 0) = 1;
142 endif 141 endif
143 xn = norm (dg .* x); 142 xn = norm (dg .* x);
144 delta = factor * xn; 143 delta = factor * xn;
145 endif 144 endif
146 145
147 # rescale if necessary 146 ## rescale if necessary
148 if (autodg) 147 if (autodg)
149 dg = max (dg, jcn); 148 dg = max (dg, jcn);
150 endif 149 endif
151 150
152 nfail = 0; 151 nfail = 0;
153 nsuc = 0; 152 nsuc = 0;
154 # inner loop 153 ## inner loop
155 while (niter <= maxiter && nfev < maxfev && ! info) 154 while (niter <= maxiter && nfev < maxfev && ! info)
156 155
157 qtf = q'*fvec; 156 qtf = q'*fvec;
158 157
159 # get TR model (dogleg) minimizer 158 ## get TR model (dogleg) minimizer
160 # in case of an overdetermined system, get lsq solution 159 ## in case of an overdetermined system, get lsq solution
161 s = - __dogleg__ (r(1:n,:), qtf(1:n), dg, delta); 160 s = - __dogleg__ (r(1:n,:), qtf(1:n), dg, delta);
162 if (pivoting) 161 if (pivoting)
163 s = p * s; 162 s = p * s;
164 endif 163 endif
165 sn = norm (dg .* s); 164 sn = norm (dg .* s);
170 169
171 fvec1 = fcn (reshape (x + s, xsiz)) (:); 170 fvec1 = fcn (reshape (x + s, xsiz)) (:);
172 fn1 = norm (fvec1); 171 fn1 = norm (fvec1);
173 172
174 if (fn1 < fn) 173 if (fn1 < fn)
175 # scaled actual reduction 174 ## scaled actual reduction
176 actred = 1 - (fn1/fn)^2; 175 actred = 1 - (fn1/fn)^2;
177 else 176 else
178 actred = -1; 177 actred = -1;
179 endif 178 endif
180 179
181 # scaled predicted reduction, and ratio 180 ## scaled predicted reduction, and ratio
182 if (pivoting) 181 if (pivoting)
183 w = qtf + r*(p'*s); 182 w = qtf + r*(p'*s);
184 else 183 else
185 w = qtf + r*s; 184 w = qtf + r*s;
186 endif 185 endif
191 else 190 else
192 prered = 0; 191 prered = 0;
193 ratio = 0; 192 ratio = 0;
194 endif 193 endif
195 194
196 # update delta 195 ## update delta
197 if (ratio < 0.1) 196 if (ratio < 0.1)
198 nsuc = 0; 197 nsuc = 0;
199 nfail ++; 198 nfail ++;
200 delta *= 0.5; 199 delta *= 0.5;
201 if (delta <= sqrt (macheps)*xn) 200 if (delta <= sqrt (macheps)*xn)
202 # trust region became uselessly small 201 ## trust region became uselessly small
203 info = -3; 202 info = -3;
204 break; 203 break;
205 endif 204 endif
206 else 205 else
207 nfail = 0; 206 nfail = 0;
212 delta = max (delta, 2*sn); 211 delta = max (delta, 2*sn);
213 endif 212 endif
214 endif 213 endif
215 214
216 if (ratio >= 1e-4) 215 if (ratio >= 1e-4)
217 # successful iteration 216 ## successful iteration
218 x += s; 217 x += s;
219 xn = norm (dg .* x); 218 xn = norm (dg .* x);
220 fvec = fvec1; 219 fvec = fvec1;
221 fn = fn1; 220 fn = fn1;
222 niter ++; 221 niter ++;
223 endif 222 endif
224 223
225 # Tests for termination conditions. A mysterious place, anything can 224 ## Tests for termination conditions. A mysterious place, anything
226 # happen if you change something here... 225 ## can happen if you change something here...
227 226
228 # The rule of thumb (which I'm not sure M*b is quite following) is that 227 ## The rule of thumb (which I'm not sure M*b is quite following)
229 # for a tolerance that depends on scaling, only 0 makes sense as a 228 ## is that for a tolerance that depends on scaling, only 0 makes
230 # default value. But 0 usually means uselessly long iterations, 229 ## sense as a default value. But 0 usually means uselessly long
231 # so we need scaling-independent tolerances wherever possible. 230 ## iterations, so we need scaling-independent tolerances wherever
232 231 ## possible.
233 # XXX: why tolf*n*xn? If abs (e) ~ abs(x) * eps is a vector of 232
234 # perturbations of x, then norm (fjac*e) <= eps*n*xn, i.e. by tolf ~ 233 ## FIXME -- why tolf*n*xn? If abs (e) ~ abs(x) * eps is a vector
235 # eps we demand as much accuracy as we can expect. 234 ## of perturbations of x, then norm (fjac*e) <= eps*n*xn, i.e. by
235 ## tolf ~ eps we demand as much accuracy as we can expect.
236 if (fn <= tolf*n*xn) 236 if (fn <= tolf*n*xn)
237 info = 1; 237 info = 1;
238 # The following tests done only after successful step. 238 ## The following tests done only after successful step.
239 elseif (actred > 0) 239 elseif (actred > 0)
240 # This one is classic. Note that we use scaled variables again, but 240 ## This one is classic. Note that we use scaled variables again,
241 # compare to scaled step, so nothing bad. 241 ## but compare to scaled step, so nothing bad.
242 if (sn <= tolx*xn) 242 if (sn <= tolx*xn)
243 info = 2; 243 info = 2;
244 # Again a classic one. It seems weird to use the same tolf for two 244 ## Again a classic one. It seems weird to use the same tolf
245 # different tests, but that's what M*b manual appears to say. 245 ## for two different tests, but that's what M*b manual appears
246 ## to say.
246 elseif (actred < tolf) 247 elseif (actred < tolf)
247 info = 3; 248 info = 3;
248 endif 249 endif
249 endif 250 endif
250 251
251 # criterion for recalculating jacobian 252 ## criterion for recalculating jacobian
252 if (nfail == 2) 253 if (nfail == 2)
253 break; 254 break;
254 endif 255 endif
255 256
256 # compute the scaled Broyden update 257 ## compute the scaled Broyden update
257 u = (fvec1 - q*w) / sn; 258 u = (fvec1 - q*w) / sn;
258 v = dg .* ((dg .* s) / sn); 259 v = dg .* ((dg .* s) / sn);
259 if (pivoting) 260 if (pivoting)
260 v = p'*v; 261 v = p'*v;
261 endif 262 endif
262 263
263 # update the QR factorization 264 ## update the QR factorization
264 [q, r] = qrupdate (q, r, u, v); 265 [q, r] = qrupdate (q, r, u, v);
265 266
266 endwhile 267 endwhile
267 endwhile 268 endwhile
268 269
269 # restore original shapes 270 ## restore original shapes
270 x = reshape (x, xsiz); 271 x = reshape (x, xsiz);
271 fvec = reshape (fvec, fsiz); 272 fvec = reshape (fvec, fsiz);
272 273
273 output.iterations = niter; 274 output.iterations = niter;
274 output.funcCount = niter + 2; 275 output.funcCount = niter + 2;
275 276
276 endfunction 277 endfunction
277 278
278 # an assistant function that evaluates a function handle and checks for bad 279 ## An assistant function that evaluates a function handle and checks for
279 # results. 280 ## bad results.
280 function fx = guarded_eval (fun, x) 281 function fx = guarded_eval (fun, x)
281 fx = fun (x); 282 fx = fun (x);
282 if (! all (isreal (fx))) 283 if (! all (isreal (fx)))
283 error ("fsolve:notreal", "fsolve: non-real value encountered"); 284 error ("fsolve: non-real value encountered");
284 elseif (any (isnan (fx))) 285 elseif (any (isnan (fx)))
285 error ("fsolve:isnan", "fsolve: NaN value encountered"); 286 error ("fsolve: NaN value encountered");
286 endif 287 endif
287 endfunction 288 endfunction
288 289
289 %!function retval = f (p) 290 %!function retval = f (p)
290 %! x = p(1); 291 %! x = p(1);