Mercurial > octave-nkf

comparison scripts/optimization/pqpnonneg.m @ 9635:36d885c4a1ac

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
implement pqpnonneg
author Jaroslav Hajek <highegg@gmail.com>
date Fri, 11 Sep 2009 11:16:38 +0200
parents
children 7c1b1c084af1
comparison
equal deleted inserted replaced
9634:da5ba66414a3 9635:36d885c4a1ac
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)