Mercurial > octave-nkf

comparison scripts/optimization/__doglegm__.m @ 9631:00958d0c4e3c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
split __dogleg__ > __doglegm__
author Jaroslav Hajek <highegg@gmail.com>
date Wed, 09 Sep 2009 14:43:06 +0200
parents
children
comparison
equal deleted inserted replaced
9630:d52e405df4f7 9631:00958d0c4e3c
1 ## Copyright (C) 2008, 2009 Jaroslav Hajek
2 ##
3 ## This file is part of Octave.
4 ##
5 ## Octave is free software; you can redistribute it and/or modify it
6 ## under the terms of the GNU General Public License as published by
7 ## the Free Software Foundation; either version 3 of the License, or (at
8 ## your option) any later version.
9 ##
10 ## Octave is distributed in the hope that it will be useful, but
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 ## General Public License for more details.
14 ##
15 ## You should have received a copy of the GNU General Public License
16 ## along with Octave; see the file COPYING. If not, see
17 ## <http://www.gnu.org/licenses/>.
18
19 ## -*- texinfo -*-
20 ## @deftypefn{Function File} {@var{x}} = __doglegm__ (@var{r}, @var{b}, @var{x}, @var{d}, @var{delta})
21 ## @end deftypefn
22
23 ## Solve the double dogleg trust-region minimization problem:
24 ## Minimize 1/2*norm(r*x)^2 subject to the constraint norm(d.*x) <= delta,
25 ## x being a convex combination of the gauss-newton and scaled gradient.
26
27 ## TODO: error checks
28 ## TODO: handle singularity, or leave it up to mldivide?
29
30 function x = __doglegm__ (r, g, d, delta)
31 ## Get Gauss-Newton direction.
32 b = r' \ g;
33 x = r \ b;
34 xn = norm (d .* x);
35 if (xn > delta)
36 ## GN is too big, get scaled gradient.
37 s = g ./ d;
38 sn = norm (s);
39 if (sn > 0)
40 ## Normalize and rescale.
41 s = (s / sn) ./ d;
42 ## Get the line minimizer in s direction.
43 tn = norm (r*s);
44 snm = (sn / tn) / tn;
45 if (snm < delta)
46 ## Get the dogleg path minimizer.
47 bn = norm (b);
48 dxn = delta/xn; snmd = snm/delta;
49 t = (bn/sn) * (bn/xn) * snmd;
50 t -= dxn * snmd^2 - sqrt ((t-dxn)^2 + (1-dxn^2)*(1-snmd^2));
51 alpha = dxn*(1-snmd^2) / t;
52 else
53 alpha = 0;
54 endif
55 else
56 alpha = delta / xn;
57 snm = 0;
58 endif
59 ## Form the appropriate convex combination.
60 x = alpha * x + ((1-alpha) * min (snm, delta)) * s;
61 endif
62 endfunction
63