|
7682
|
1 ## Copyright (C) 2008 Bill Denney |
|
|
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} =} lsqnonneg (@var{c}, @var{d}) |
|
|
21 ## @deftypefnx {Function File} {@var{x} =} lsqnonneg (@var{c}, @var{d}, @var{x0}) |
|
|
22 ## @deftypefnx {Function File} {[@var{x}, @var{resnorm}] =} lsqnonneg (@dots{}) |
|
|
23 ## @deftypefnx {Function File} {[@var{x}, @var{resnorm}, @var{residual}] =} lsqnonneg (@dots{}) |
|
|
24 ## @deftypefnx {Function File} {[@var{x}, @var{resnorm}, @var{residual}, @var{exitflag}] =} lsqnonneg (@dots{}) |
|
|
25 ## @deftypefnx {Function File} {[@var{x}, @var{resnorm}, @var{residual}, @var{exitflag}, @var{output}] =} lsqnonneg (@dots{}) |
|
|
26 ## @deftypefnx {Function File} {[@var{x}, @var{resnorm}, @var{residual}, @var{exitflag}, @var{output}, @var{lambda}] =} lsqnonneg (@dots{}) |
|
|
27 ## Minimize @code{norm (@var{c}*@var{x}-d)} subject to @code{@var{x} >= |
|
|
28 ## 0}. @var{c} and @var{d} must be real. @var{x0} is an optional |
|
|
29 ## initial guess for @var{x}. |
|
|
30 ## |
|
|
31 ## Outputs: |
|
|
32 ## @itemize @bullet |
|
|
33 ## @item resnorm |
|
|
34 ## |
|
|
35 ## The squared 2-norm of the residual: norm(@var{c}*@var{x}-@var{d})^2 |
|
|
36 ## @item residual |
|
|
37 ## |
|
|
38 ## The residual: @var{d}-@var{c}*@var{x} |
|
|
39 ## @item exitflag |
|
|
40 ## |
|
|
41 ## An indicator of convergence. 0 indicates that the iteration count |
|
|
42 ## was exceeded, and therefore convergence was not reached; >0 indicates |
|
|
43 ## that the algorithm converged. (The algorithm is stable and will |
|
|
44 ## converge given enough iterations.) |
|
|
45 ## @item output |
|
|
46 ## |
|
|
47 ## A structure with two fields: |
|
|
48 ## @itemize @bullet |
|
|
49 ## @item "algorithm": The algorithm used ("nnls") |
|
|
50 ## @item "iterations": The number of iterations taken. |
|
|
51 ## @end itemize |
|
|
52 ## @item lambda |
|
|
53 ## |
|
|
54 ## Not implemented. |
|
|
55 ## @end itemize |
|
|
56 ## @seealso{optimset} |
|
|
57 ## @end deftypefn |
|
|
58 |
|
|
59 ## This is implemented from Lawson and Hanson's 1973 algorithm on page |
|
|
60 |
|
|
61 function [x, resnorm, residual, exitflag, output, lambda] = lsqnonneg (c, d, x = [], options = []) |
|
|
62 |
|
|
63 if (isempty (x)) |
|
|
64 ## initial guess is 0s |
|
|
65 x = zeros (columns (c), 1); |
|
|
66 endif |
|
|
67 |
|
|
68 if (isempty (options)) |
|
|
69 ## FIXME: Optimset should be updated |
|
|
70 ## options = optimset (); |
|
|
71 options = struct ("maxiter", 1e5, "tolx", 1e-8); |
|
|
72 endif |
|
|
73 |
|
|
74 ## Initialize the values |
|
|
75 p = []; |
|
|
76 z = 1:numel (x); |
|
|
77 ## compute the gradient |
|
|
78 w = c'*(d - c*x); |
|
|
79 |
|
|
80 iter = 0; |
|
|
81 while (! isempty (z) && any (w(z) > 0) && iter < options.maxiter) |
|
|
82 idx = find (w == max (w)); |
|
|
83 if (numel (idx) > 1) |
|
|
84 warning ("lsqnonneg:nonunique", |
|
|
85 "A non-unique solution may be returned due to equal gradients."); |
|
|
86 idx = idx(1); |
|
|
87 endif |
|
|
88 p(end+1) = z(idx); |
|
|
89 z(idx) = []; |
|
|
90 |
|
|
91 newx = false; |
|
|
92 while (! newx && iter < options.maxiter) |
|
|
93 iter++; |
|
|
94 |
|
|
95 cpos = c; |
|
|
96 cpos(:,z) = 0; |
|
|
97 ## find min norm solution to the positive matrix |
|
|
98 [cpos_q, cpos_r] = qr (cpos, 0); |
|
|
99 xtmp = (cpos_r\cpos_q')*d; |
|
|
100 xtmp(z) = 0; |
|
|
101 if (all (xtmp >= 0)) |
|
|
102 ## complete the inner loop |
|
|
103 newx = true; |
|
|
104 x = xtmp; |
|
|
105 else |
|
|
106 mask = (xtmp < 0); |
|
|
107 alpha = min(x(mask)./(x(mask) - xtmp(mask))); |
|
|
108 x = x + alpha*(xtmp - x); |
|
|
109 idx = find (x == 0); |
|
|
110 p(ismember (p, x)) = []; |
|
|
111 z = [z idx]; |
|
|
112 endif |
|
|
113 endwhile |
|
|
114 w = c'*(d - c*x); |
|
|
115 endwhile |
|
|
116 |
|
|
117 ## generate the additional output arguments |
|
|
118 if (nargout > 1) |
|
|
119 resnorm = norm (C*x-d) ^ 2; |
|
|
120 endif |
|
|
121 if (nargout > 2) |
|
|
122 residual = d-C*x; |
|
|
123 endif |
|
|
124 exitflag = iter; |
|
|
125 if (nargout > 3 && iter >= options.maxiter) |
|
|
126 exitflag = 0; |
|
|
127 endif |
|
|
128 if (nargout > 4) |
|
|
129 output = struct ("algorithm", "nnls", "iterations", iter); |
|
|
130 endif |
|
|
131 if (nargout > 5) |
|
|
132 ## FIXME |
|
|
133 error ("lsqnonneg: lambda is not yet implemented"); |
|
|
134 endif |
|
|
135 |
|
|
136 endfunction |
|
|
137 |
|
|
138 ## Tests |
|
|
139 %!test |
|
|
140 %! C = [1 0;0 1;2 1]; |
|
|
141 %! d = [1;3;-2]; |
|
|
142 %! assert (lsqnonneg (C, d), [0;0.5], 100*eps) |
|
|
143 |
|
|
144 %!test |
|
|
145 %! C = [0.0372 0.2869;0.6861 0.7071;0.6233 0.6245;0.6344 0.6170]; |
|
|
146 %! d = [0.8587;0.1781;0.0747;0.8405]; |
|
|
147 %! xnew = [0;0.6929]; |
|
|
148 %! assert (lsqnonneg (C, d), xnew, 0.0001) |
