Mercurial > octave

annotate scripts/optimization/lsqnonneg.m @ 8407:096c22ce2b0b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
style fixes in previous patch
author Jaroslav Hajek <highegg@gmail.com>
date Mon, 15 Dec 2008 13:49:11 +0100
parents 0b7566aea627
children cadc73247d65
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
7682
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
1 ## Copyright (C) 2008 Bill Denney
8406
0b7566aea627 fix & optimize lsqnonneg
Jaroslav Hajek <highegg@gmail.com>
parents: 8304
diff changeset
2 ## Copyright (C) 2008 Jaroslav Hajek
7682
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
3 ##
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
4 ## This file is part of Octave.
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
5 ##
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
6 ## Octave is free software; you can redistribute it and/or modify it
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
7 ## under the terms of the GNU General Public License as published by
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
8 ## the Free Software Foundation; either version 3 of the License, or (at
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
9 ## your option) any later version.
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
10 ##
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
11 ## Octave is distributed in the hope that it will be useful, but
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
14 ## General Public License for more details.
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
15 ##
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
16 ## You should have received a copy of the GNU General Public License
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
17 ## along with Octave; see the file COPYING. If not, see
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
18 ## <http://www.gnu.org/licenses/>.
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
19
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
20 ## -*- texinfo -*-
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
21 ## @deftypefn {Function File} {@var{x} =} lsqnonneg (@var{c}, @var{d})
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
22 ## @deftypefnx {Function File} {@var{x} =} lsqnonneg (@var{c}, @var{d}, @var{x0})
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
23 ## @deftypefnx {Function File} {[@var{x}, @var{resnorm}] =} lsqnonneg (@dots{})
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
24 ## @deftypefnx {Function File} {[@var{x}, @var{resnorm}, @var{residual}] =} lsqnonneg (@dots{})
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
25 ## @deftypefnx {Function File} {[@var{x}, @var{resnorm}, @var{residual}, @var{exitflag}] =} lsqnonneg (@dots{})
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
26 ## @deftypefnx {Function File} {[@var{x}, @var{resnorm}, @var{residual}, @var{exitflag}, @var{output}] =} lsqnonneg (@dots{})
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
27 ## @deftypefnx {Function File} {[@var{x}, @var{resnorm}, @var{residual}, @var{exitflag}, @var{output}, @var{lambda}] =} lsqnonneg (@dots{})
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
28 ## Minimize @code{norm (@var{c}*@var{x}-d)} subject to @code{@var{x} >=
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
29 ## 0}. @var{c} and @var{d} must be real. @var{x0} is an optional
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
30 ## initial guess for @var{x}.
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
31 ##
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
32 ## Outputs:
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
33 ## @itemize @bullet
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
34 ## @item resnorm
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
35 ##
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
36 ## The squared 2-norm of the residual: norm(@var{c}*@var{x}-@var{d})^2
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
37 ## @item residual
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
38 ##
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
39 ## The residual: @var{d}-@var{c}*@var{x}
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
40 ## @item exitflag
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
41 ##
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
42 ## An indicator of convergence. 0 indicates that the iteration count
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
43 ## was exceeded, and therefore convergence was not reached; >0 indicates
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
44 ## that the algorithm converged. (The algorithm is stable and will
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
45 ## converge given enough iterations.)
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
46 ## @item output
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
47 ##
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
48 ## A structure with two fields:
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
49 ## @itemize @bullet
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
50 ## @item "algorithm": The algorithm used ("nnls")
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
51 ## @item "iterations": The number of iterations taken.
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
52 ## @end itemize
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
53 ## @item lambda
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
54 ##
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
55 ## Not implemented.
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
56 ## @end itemize
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
57 ## @seealso{optimset}
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
58 ## @end deftypefn
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
59
7697
0bdfff62cc49 lsqnonneg: use optimset, correctly index Z and P in main loop
bill@denney.ws
parents: 7682
diff changeset
60 ## This is implemented from Lawson and Hanson's 1973 algorithm on page
0bdfff62cc49 lsqnonneg: use optimset, correctly index Z and P in main loop
bill@denney.ws
parents: 7682
diff changeset
61 ## 161 of Solving Least Squares Problems.
7682
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
62
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
63 function [x, resnorm, residual, exitflag, output, lambda] = lsqnonneg (c, d, x = [], options = [])
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
64
7697
0bdfff62cc49 lsqnonneg: use optimset, correctly index Z and P in main loop
bill@denney.ws
parents: 7682
diff changeset
65 ## Lawson-Hanson Step 1 (LH1): initialize the variables.
7682
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
66 if (isempty (x))
7697
0bdfff62cc49 lsqnonneg: use optimset, correctly index Z and P in main loop
bill@denney.ws
parents: 7682
diff changeset
67 ## Initial guess is 0s.
7682
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
68 x = zeros (columns (c), 1);
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
69 endif
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
70
8406
0b7566aea627 fix & optimize lsqnonneg
Jaroslav Hajek <highegg@gmail.com>
parents: 8304
diff changeset
71
8304
eeaee297c0da modify optimset & implement optimget
Jaroslav Hajek <highegg@gmail.com>
parents: 7697
diff changeset
72 MaxIter = optimget (options, "MaxIter", 1e5);
7682
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
73
7697
0bdfff62cc49 lsqnonneg: use optimset, correctly index Z and P in main loop
bill@denney.ws
parents: 7682
diff changeset
74 ## Initialize the values.
8406
0b7566aea627 fix & optimize lsqnonneg
Jaroslav Hajek <highegg@gmail.com>
parents: 8304
diff changeset
75 p = false (1, numel (x));
0b7566aea627 fix & optimize lsqnonneg
Jaroslav Hajek <highegg@gmail.com>
parents: 8304
diff changeset
76 z = ~p;
0b7566aea627 fix & optimize lsqnonneg
Jaroslav Hajek <highegg@gmail.com>
parents: 8304
diff changeset
77 ## If the problem is significantly over-determined, preprocess it using a
0b7566aea627 fix & optimize lsqnonneg
Jaroslav Hajek <highegg@gmail.com>
parents: 8304
diff changeset
78 ## QR factorization first.
0b7566aea627 fix & optimize lsqnonneg
Jaroslav Hajek <highegg@gmail.com>
parents: 8304
diff changeset
79 if (rows (c) >= 1.5 * columns (c))
0b7566aea627 fix & optimize lsqnonneg
Jaroslav Hajek <highegg@gmail.com>
parents: 8304
diff changeset
80 [q, r] = qr (c, 0);
0b7566aea627 fix & optimize lsqnonneg
Jaroslav Hajek <highegg@gmail.com>
parents: 8304
diff changeset
81 d = q'*d;
0b7566aea627 fix & optimize lsqnonneg
Jaroslav Hajek <highegg@gmail.com>
parents: 8304
diff changeset
82 c = r;
0b7566aea627 fix & optimize lsqnonneg
Jaroslav Hajek <highegg@gmail.com>
parents: 8304
diff changeset
83 clear q
0b7566aea627 fix & optimize lsqnonneg
Jaroslav Hajek <highegg@gmail.com>
parents: 8304
diff changeset
84 endif
7697
0bdfff62cc49 lsqnonneg: use optimset, correctly index Z and P in main loop
bill@denney.ws
parents: 7682
diff changeset
85 ## LH2: compute the gradient.
7682
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
86 w = c'*(d - c*x);
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
87
8406
0b7566aea627 fix & optimize lsqnonneg
Jaroslav Hajek <highegg@gmail.com>
parents: 8304
diff changeset
88 xtmp = zeros (columns (c), 1);
0b7566aea627 fix & optimize lsqnonneg
Jaroslav Hajek <highegg@gmail.com>
parents: 8304
diff changeset
89
7682
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
90 iter = 0;
7697
0bdfff62cc49 lsqnonneg: use optimset, correctly index Z and P in main loop
bill@denney.ws
parents: 7682
diff changeset
91 ## LH3: test for completion.
8406
0b7566aea627 fix & optimize lsqnonneg
Jaroslav Hajek <highegg@gmail.com>
parents: 8304
diff changeset
92 while (any (z) && any (w(z) > 0) && iter < MaxIter)
7697
0bdfff62cc49 lsqnonneg: use optimset, correctly index Z and P in main loop
bill@denney.ws
parents: 7682
diff changeset
93 ## LH4: find the maximum gradient.
7682
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
94 idx = find (w == max (w));
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
95 if (numel (idx) > 1)
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
96 warning ("lsqnonneg:nonunique",
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
97 "A non-unique solution may be returned due to equal gradients.");
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
98 idx = idx(1);
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
99 endif
7697
0bdfff62cc49 lsqnonneg: use optimset, correctly index Z and P in main loop
bill@denney.ws
parents: 7682
diff changeset
100 ## LH5: move the index from Z to P.
8406
0b7566aea627 fix & optimize lsqnonneg
Jaroslav Hajek <highegg@gmail.com>
parents: 8304
diff changeset
101 z(idx) = false;
0b7566aea627 fix & optimize lsqnonneg
Jaroslav Hajek <highegg@gmail.com>
parents: 8304
diff changeset
102 p(idx) = true;
7682
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
103
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
104 newx = false;
8304
eeaee297c0da modify optimset & implement optimget
Jaroslav Hajek <highegg@gmail.com>
parents: 7697
diff changeset
105 while (! newx && iter < MaxIter)
7682
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
106 iter++;
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
107
7697
0bdfff62cc49 lsqnonneg: use optimset, correctly index Z and P in main loop
bill@denney.ws
parents: 7682
diff changeset
108 ## LH6: compute the positive matrix and find the min norm solution
0bdfff62cc49 lsqnonneg: use optimset, correctly index Z and P in main loop
bill@denney.ws
parents: 7682
diff changeset
109 ## of the positive problem.
0bdfff62cc49 lsqnonneg: use optimset, correctly index Z and P in main loop
bill@denney.ws
parents: 7682
diff changeset
110 ## Find min norm solution to the positive matrix.
8406
0b7566aea627 fix & optimize lsqnonneg
Jaroslav Hajek <highegg@gmail.com>
parents: 8304
diff changeset
111 xtmp(:) = 0;
0b7566aea627 fix & optimize lsqnonneg
Jaroslav Hajek <highegg@gmail.com>
parents: 8304
diff changeset
112 xtmp(p) = c(:,p) \ d;
7682
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
113 if (all (xtmp >= 0))
7697
0bdfff62cc49 lsqnonneg: use optimset, correctly index Z and P in main loop
bill@denney.ws
parents: 7682
diff changeset
114 ## LH7: tmp solution found, iterate.
7682
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
115 newx = true;
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
116 x = xtmp;
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
117 else
7697
0bdfff62cc49 lsqnonneg: use optimset, correctly index Z and P in main loop
bill@denney.ws
parents: 7682
diff changeset
118 ## LH8, LH9: find the scaling factor and adjust X.
7682
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
119 mask = (xtmp < 0);
7697
0bdfff62cc49 lsqnonneg: use optimset, correctly index Z and P in main loop
bill@denney.ws
parents: 7682
diff changeset
120 alpha = min (x(mask)./(x(mask) - xtmp(mask)));
0bdfff62cc49 lsqnonneg: use optimset, correctly index Z and P in main loop
bill@denney.ws
parents: 7682
diff changeset
121 ## LH10: adjust X.
7682
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
122 x = x + alpha*(xtmp - x);
7697
0bdfff62cc49 lsqnonneg: use optimset, correctly index Z and P in main loop
bill@denney.ws
parents: 7682
diff changeset
123 ## LH11: move from P to Z all X == 0.
8407
096c22ce2b0b style fixes in previous patch
Jaroslav Hajek <highegg@gmail.com>
parents: 8406
diff changeset
124 z |= (x == 0);
096c22ce2b0b style fixes in previous patch
Jaroslav Hajek <highegg@gmail.com>
parents: 8406
diff changeset
125 p = ~z;
7682
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
126 endif
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
127 endwhile
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
128 w = c'*(d - c*x);
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
129 endwhile
7697
0bdfff62cc49 lsqnonneg: use optimset, correctly index Z and P in main loop
bill@denney.ws
parents: 7682
diff changeset
130 ## LH12: complete.
7682
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
131
7697
0bdfff62cc49 lsqnonneg: use optimset, correctly index Z and P in main loop
bill@denney.ws
parents: 7682
diff changeset
132 ## Generate the additional output arguments.
7682
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
133 if (nargout > 1)
8406
0b7566aea627 fix & optimize lsqnonneg
Jaroslav Hajek <highegg@gmail.com>
parents: 8304
diff changeset
134 resnorm = norm (c*x - d) ^ 2;
7682
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
135 endif
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
136 if (nargout > 2)
8406
0b7566aea627 fix & optimize lsqnonneg
Jaroslav Hajek <highegg@gmail.com>
parents: 8304
diff changeset
137 residual = d - c*x;
7682
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
138 endif
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
139 exitflag = iter;
8304
eeaee297c0da modify optimset & implement optimget
Jaroslav Hajek <highegg@gmail.com>
parents: 7697
diff changeset
140 if (nargout > 3 && iter >= MaxIter)
7682
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
141 exitflag = 0;
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
142 endif
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
143 if (nargout > 4)
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
144 output = struct ("algorithm", "nnls", "iterations", iter);
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
145 endif
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
146 if (nargout > 5)
8406
0b7566aea627 fix & optimize lsqnonneg
Jaroslav Hajek <highegg@gmail.com>
parents: 8304
diff changeset
147 lambda = w;
7682
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
148 endif
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
149
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
150 endfunction
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
151
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
152 ## Tests
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
153 %!test
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
154 %! C = [1 0;0 1;2 1];
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
155 %! d = [1;3;-2];
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
156 %! assert (lsqnonneg (C, d), [0;0.5], 100*eps)
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
157
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
158 %!test
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
159 %! C = [0.0372 0.2869;0.6861 0.7071;0.6233 0.6245;0.6344 0.6170];
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
160 %! d = [0.8587;0.1781;0.0747;0.8405];
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
161 %! xnew = [0;0.6929];
f7474c83c782 lsqnonneg: new function
bill@denney.ws
parents:
diff changeset
162 %! assert (lsqnonneg (C, d), xnew, 0.0001)