Mercurial > octave
annotate scripts/optimization/lsqnonneg.m @ 8507:cadc73247d65
style fixes
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Tue, 13 Jan 2009 14:08:36 -0500 |
parents | 096c22ce2b0b |
children | a6c1aa6f5915 |
rev | line source |
---|---|
7682 | 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 | 3 ## |
4 ## This file is part of Octave. | |
5 ## | |
6 ## Octave is free software; you can redistribute it and/or modify it | |
7 ## under the terms of the GNU General Public License as published by | |
8 ## the Free Software Foundation; either version 3 of the License, or (at | |
9 ## your option) any later version. | |
10 ## | |
11 ## Octave is distributed in the hope that it will be useful, but | |
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of | |
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
14 ## General Public License for more details. | |
15 ## | |
16 ## You should have received a copy of the GNU General Public License | |
17 ## along with Octave; see the file COPYING. If not, see | |
18 ## <http://www.gnu.org/licenses/>. | |
19 | |
20 ## -*- texinfo -*- | |
21 ## @deftypefn {Function File} {@var{x} =} lsqnonneg (@var{c}, @var{d}) | |
22 ## @deftypefnx {Function File} {@var{x} =} lsqnonneg (@var{c}, @var{d}, @var{x0}) | |
23 ## @deftypefnx {Function File} {[@var{x}, @var{resnorm}] =} lsqnonneg (@dots{}) | |
24 ## @deftypefnx {Function File} {[@var{x}, @var{resnorm}, @var{residual}] =} lsqnonneg (@dots{}) | |
25 ## @deftypefnx {Function File} {[@var{x}, @var{resnorm}, @var{residual}, @var{exitflag}] =} lsqnonneg (@dots{}) | |
26 ## @deftypefnx {Function File} {[@var{x}, @var{resnorm}, @var{residual}, @var{exitflag}, @var{output}] =} lsqnonneg (@dots{}) | |
27 ## @deftypefnx {Function File} {[@var{x}, @var{resnorm}, @var{residual}, @var{exitflag}, @var{output}, @var{lambda}] =} lsqnonneg (@dots{}) | |
28 ## Minimize @code{norm (@var{c}*@var{x}-d)} subject to @code{@var{x} >= | |
29 ## 0}. @var{c} and @var{d} must be real. @var{x0} is an optional | |
30 ## initial guess for @var{x}. | |
31 ## | |
32 ## Outputs: | |
33 ## @itemize @bullet | |
34 ## @item resnorm | |
35 ## | |
36 ## The squared 2-norm of the residual: norm(@var{c}*@var{x}-@var{d})^2 | |
37 ## @item residual | |
38 ## | |
39 ## The residual: @var{d}-@var{c}*@var{x} | |
40 ## @item exitflag | |
41 ## | |
42 ## An indicator of convergence. 0 indicates that the iteration count | |
43 ## was exceeded, and therefore convergence was not reached; >0 indicates | |
44 ## that the algorithm converged. (The algorithm is stable and will | |
45 ## converge given enough iterations.) | |
46 ## @item output | |
47 ## | |
48 ## A structure with two fields: | |
49 ## @itemize @bullet | |
50 ## @item "algorithm": The algorithm used ("nnls") | |
51 ## @item "iterations": The number of iterations taken. | |
52 ## @end itemize | |
53 ## @item lambda | |
54 ## | |
55 ## Not implemented. | |
56 ## @end itemize | |
57 ## @seealso{optimset} | |
58 ## @end deftypefn | |
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 | 62 |
63 function [x, resnorm, residual, exitflag, output, lambda] = lsqnonneg (c, d, x = [], options = []) | |
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 | 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 | 68 x = zeros (columns (c), 1); |
69 endif | |
70 | |
8406
0b7566aea627
fix & optimize lsqnonneg
Jaroslav Hajek <highegg@gmail.com>
parents:
8304
diff
changeset
|
71 |
8507 | 72 max_iter = optimget (options, "MaxIter", 1e5); |
7682 | 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)); |
8507 | 76 z = !p; |
8406
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 | 86 w = c'*(d - c*x); |
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 | 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. |
8507 | 92 while (any (z) && any (w(z) > 0) && iter < max_iter) |
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 | 94 idx = find (w == max (w)); |
95 if (numel (idx) > 1) | |
96 warning ("lsqnonneg:nonunique", | |
97 "A non-unique solution may be returned due to equal gradients."); | |
98 idx = idx(1); | |
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 | 103 |
104 newx = false; | |
8507 | 105 while (! newx && iter < max_iter) |
7682 | 106 iter++; |
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 | 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 | 115 newx = true; |
116 x = xtmp; | |
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 | 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 | 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); |
8507 | 125 p = !z; |
7682 | 126 endif |
127 endwhile | |
128 w = c'*(d - c*x); | |
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 | 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 | 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 | 135 endif |
136 if (nargout > 2) | |
8406
0b7566aea627
fix & optimize lsqnonneg
Jaroslav Hajek <highegg@gmail.com>
parents:
8304
diff
changeset
|
137 residual = d - c*x; |
7682 | 138 endif |
139 exitflag = iter; | |
8507 | 140 if (nargout > 3 && iter >= max_iter) |
7682 | 141 exitflag = 0; |
142 endif | |
143 if (nargout > 4) | |
144 output = struct ("algorithm", "nnls", "iterations", iter); | |
145 endif | |
146 if (nargout > 5) | |
8406
0b7566aea627
fix & optimize lsqnonneg
Jaroslav Hajek <highegg@gmail.com>
parents:
8304
diff
changeset
|
147 lambda = w; |
7682 | 148 endif |
149 | |
150 endfunction | |
151 | |
152 ## Tests | |
153 %!test | |
154 %! C = [1 0;0 1;2 1]; | |
155 %! d = [1;3;-2]; | |
156 %! assert (lsqnonneg (C, d), [0;0.5], 100*eps) | |
157 | |
158 %!test | |
159 %! C = [0.0372 0.2869;0.6861 0.7071;0.6233 0.6245;0.6344 0.6170]; | |
160 %! d = [0.8587;0.1781;0.0747;0.8405]; | |
161 %! xnew = [0;0.6929]; | |
162 %! assert (lsqnonneg (C, d), xnew, 0.0001) |