9635
|
1 ## Copyright (C) 2008 Bill Denney |
|
2 ## Copyright (C) 2008 Jaroslav Hajek |
|
3 ## Copyright (C) 2009 VZLU Prague |
|
4 ## |
|
5 ## This file is part of Octave. |
|
6 ## |
|
7 ## Octave is free software; you can redistribute it and/or modify it |
|
8 ## under the terms of the GNU General Public License as published by |
|
9 ## the Free Software Foundation; either version 3 of the License, or (at |
|
10 ## your option) any later version. |
|
11 ## |
|
12 ## Octave is distributed in the hope that it will be useful, but |
|
13 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
14 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
15 ## General Public License for more details. |
|
16 ## |
|
17 ## You should have received a copy of the GNU General Public License |
|
18 ## along with Octave; see the file COPYING. If not, see |
|
19 ## <http://www.gnu.org/licenses/>. |
|
20 |
|
21 ## -*- texinfo -*- |
|
22 ## @deftypefn {Function File} {@var{x} =} pqpnonneg (@var{c}, @var{d}) |
|
23 ## @deftypefnx {Function File} {@var{x} =} pqpnonneg (@var{c}, @var{d}, @var{x0}) |
|
24 ## @deftypefnx {Function File} {[@var{x}, @var{minval}] =} pqpnonneg (@dots{}) |
|
25 ## @deftypefnx {Function File} {[@var{x}, @var{minval}, @var{exitflag}] =} pqpnonneg (@dots{}) |
|
26 ## @deftypefnx {Function File} {[@var{x}, @var{minval}, @var{exitflag}, @var{output}] =} pqpnonneg (@dots{}) |
|
27 ## @deftypefnx {Function File} {[@var{x}, @var{minval}, @var{exitflag}, @var{output}, @var{lambda}] =} pqpnonneg (@dots{}) |
|
28 ## Minimize @code{1/2*x'*c*x + d'*x} subject to @code{@var{x} >= |
|
29 ## 0}. @var{c} and @var{d} must be real, and @var{c} must be symmetric and positive definite. |
|
30 ## @var{x0} is an optional initial guess for @var{x}. |
|
31 ## |
|
32 ## Outputs: |
|
33 ## @itemize @bullet |
|
34 ## @item minval |
|
35 ## |
|
36 ## The minimum attained model value, 1/2*xmin'*c*xmin + d'*xmin |
|
37 ## @item exitflag |
|
38 ## |
|
39 ## An indicator of convergence. 0 indicates that the iteration count |
|
40 ## was exceeded, and therefore convergence was not reached; >0 indicates |
|
41 ## that the algorithm converged. (The algorithm is stable and will |
|
42 ## converge given enough iterations.) |
|
43 ## @item output |
|
44 ## |
|
45 ## A structure with two fields: |
|
46 ## @itemize @bullet |
|
47 ## @item "algorithm": The algorithm used ("nnls") |
|
48 ## @item "iterations": The number of iterations taken. |
|
49 ## @end itemize |
|
50 ## @item lambda |
|
51 ## |
|
52 ## Not implemented. |
|
53 ## @end itemize |
|
54 ## @seealso{optimset, lsqnonneg, qp} |
|
55 ## @end deftypefn |
|
56 |
|
57 ## PKG_ADD: __all_opts__ ("pqpnonneg"); |
|
58 |
|
59 ## This is analogical to the lsqnonneg implementation, which is |
|
60 ## implemented from Lawson and Hanson's 1973 algorithm on page |
|
61 ## 161 of Solving Least Squares Problems. |
|
62 ## It shares the convergence guarantees. |
|
63 |
|
64 function [x, minval, exitflag, output, lambda] = pqpnonneg (c, d, x = [], options = struct ()) |
|
65 |
|
66 if (nargin == 1 && ischar (c) && strcmp (c, 'defaults')) |
|
67 x = optimset ("MaxIter", 1e5); |
|
68 return |
|
69 endif |
|
70 |
|
71 if (! (nargin >= 2 && nargin <= 4 && ismatrix (c) && ismatrix (d) && isstruct (options))) |
|
72 print_usage (); |
|
73 endif |
|
74 |
|
75 ## Lawson-Hanson Step 1 (LH1): initialize the variables. |
|
76 m = rows (c); |
|
77 n = columns (c); |
|
78 if (m != n) |
|
79 error ("matrix must be square"); |
|
80 endif |
|
81 |
|
82 if (isempty (x)) |
|
83 ## Initial guess is 0s. |
|
84 x = zeros (n, 1); |
|
85 else |
|
86 ## ensure nonnegative guess. |
|
87 x = max (x, 0); |
|
88 endif |
|
89 |
|
90 max_iter = optimget (options, "MaxIter", 1e5); |
|
91 |
|
92 ## Initialize P, according to zero pattern of x. |
|
93 p = find (x > 0).'; |
|
94 ## Initialize the Cholesky factorization. |
|
95 r = chol (c(p, p)); |
|
96 usechol = true; |
|
97 |
|
98 iter = 0; |
|
99 |
|
100 ## LH3: test for completion. |
|
101 while (iter < max_iter) |
|
102 while (iter < max_iter) |
|
103 iter++; |
|
104 |
|
105 ## LH6: compute the positive matrix and find the min norm solution |
|
106 ## of the positive problem. |
|
107 if (usechol) |
|
108 xtmp = -(r \ (r' \ d(p))); |
|
109 else |
|
110 xtmp = -(c(p,p) \ d(p)); |
|
111 endif |
|
112 idx = find (xtmp < 0); |
|
113 |
|
114 if (isempty (idx)) |
|
115 ## LH7: tmp solution found, iterate. |
|
116 x(:) = 0; |
|
117 x(p) = xtmp; |
|
118 break; |
|
119 else |
|
120 ## LH8, LH9: find the scaling factor. |
|
121 pidx = p(idx); |
|
122 sf = x(pidx)./(x(pidx) - xtmp(idx)); |
|
123 alpha = min (sf); |
|
124 ## LH10: adjust X. |
|
125 xx = zeros (n, 1); |
|
126 xx(p) = xtmp; |
|
127 x += alpha*(xx - x); |
|
128 ## LH11: move from P to Z all X == 0. |
|
129 ## This corresponds to those indices where minimum of sf is attained. |
|
130 idx = idx (sf == alpha); |
|
131 p(idx) = []; |
|
132 if (usechol) |
|
133 ## update the Cholesky factorization. |
|
134 r = choldelete (r, idx); |
|
135 endif |
|
136 endif |
|
137 endwhile |
|
138 |
|
139 ## compute the gradient. |
|
140 w = -(d + c*x); |
|
141 w(p) = []; |
|
142 if (! any (w > 0)) |
|
143 if (usechol) |
|
144 ## verify the solution achieved using qr updating. |
|
145 ## in the best case, this should only take a single step. |
|
146 usechol = false; |
|
147 continue; |
|
148 else |
|
149 ## we're finished. |
|
150 break; |
|
151 endif |
|
152 endif |
|
153 |
|
154 ## find the maximum gradient. |
|
155 idx = find (w == max (w)); |
|
156 if (numel (idx) > 1) |
|
157 warning ("pqpnonneg:nonunique", |
|
158 "A non-unique solution may be returned due to equal gradients."); |
|
159 idx = idx(1); |
|
160 endif |
|
161 ## move the index from Z to P. Keep P sorted. |
|
162 z = [1:n]; z(p) = []; |
|
163 zidx = z(idx); |
|
164 jdx = 1 + lookup (p, zidx); |
|
165 p = [p(1:jdx-1), zidx, p(jdx:end)]; |
|
166 if (usechol) |
|
167 ## insert the column into the Cholesky factorization. |
|
168 r = cholinsert (r, jdx, c(p,zidx)); |
|
169 endif |
|
170 |
|
171 endwhile |
|
172 ## LH12: complete. |
|
173 |
|
174 ## Generate the additional output arguments. |
|
175 if (nargout > 1) |
|
176 minval = 1/2*(x'*c*x) + d'*x; |
|
177 endif |
|
178 exitflag = iter; |
|
179 if (nargout > 2 && iter >= max_iter) |
|
180 exitflag = 0; |
|
181 endif |
|
182 if (nargout > 3) |
|
183 output = struct ("algorithm", "nnls-pqp", "iterations", iter); |
|
184 endif |
|
185 if (nargout > 4) |
|
186 lambda = zeros (size (x)); |
|
187 lambda(p) = w; |
|
188 endif |
|
189 |
|
190 endfunction |
|
191 |
|
192 ## Tests |
|
193 %!test |
|
194 %! C = [5 2;2 2]; |
|
195 %! d = [3; -1]; |
|
196 %! assert (pqpnonneg (C, d), [0;0.5], 100*eps) |
|
197 |
|
198 ## Test equivalence of lsq and pqp |
|
199 %!test |
|
200 %! C = rand (20, 10); |
|
201 %! d = rand (20, 1); |
|
202 %! assert (pqpnonneg (C'*C, -C'*d), lsqnonneg (C, d), 100*eps) |