Mercurial > hg > octave-lyh

annotate scripts/optimization/__dogleg__.m @ 8306:43795cf108d0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
initial implementation of fsolve remove old fsolve code
author Jaroslav Hajek <highegg@gmail.com>
date Sun, 28 Sep 2008 21:09:35 +0200
parents
children bc982528de11
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
8306
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
1 ## Copyright (C) 2008 Jaroslav Hajek <highegg@gmail.com>
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 ## -*- texinfo -*-
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
20 ## @deftypefn{Function File} {@var{x}} = __dogleg__ (@var{r}, @var{b}, @var{x}, @var{d}, @var{delta})
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
21 ## Solves the double dogleg trust-region problem:
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
22 ## Minimize norm(r*x-b) subject to the constraint norm(d.*x) <= delta,
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
23 ## x being a convex combination of the gauss-newton and scaled gradient.
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
24 ## @end deftypefn
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
25
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
26 ## TODO: error checks
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
27 ## TODO: handle singularity, or leave it up to mldivide?
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
28
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
29 function x = __dogleg__ (r, b, d, delta)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
30 # get Gauss-Newton direction
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
31 x = r \ b;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
32 xn = norm (d .* x);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
33 if (xn > delta)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
34 # GN is too big, get scaled gradient
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
35 s = (r' * b) ./ d;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
36 sn = norm (s);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
37 if (sn > 0)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
38 # normalize and rescale
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
39 s = (s / sn) ./ d;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
40 # get the line minimizer in s direction
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
41 tn = norm (r*s);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
42 snm = (sn / tn) / tn;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
43 if (snm < delta)
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
44 # get the dogleg path minimizer
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
45 bn = norm (b);
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
46 dxn = delta/xn; snmd = snm/delta;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
47 t = (bn/sn) * (bn/xn) * snmd;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
48 t -= dxn * snmd^2 - sqrt ((t-dxn)^2 + (1-dxn^2)*(1-snmd^2));
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
49 alpha = dxn*(1-snmd^2) / t;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
50 else
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
51 alpha = 0;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
52 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
53 else
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
54 snm = 0;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
55 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
56 # form the appropriate convex combination
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
57 x = alpha * x + ((1-alpha) * min (snm, delta)) * s;
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
58 endif
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
59 endfunction
43795cf108d0 initial implementation of fsolve
Jaroslav Hajek <highegg@gmail.com>
parents:
diff changeset
60