Mercurial > octave
annotate scripts/optimization/pqpnonneg.m @ 22769:9bd2ea5703e2
lsqnonneg.m: Overhaul function.
* lsqnonneg.m: Use @var entries in docstring. Change to @table rather than
@itemize @bullet. Rename function input x to x0 to match documentation. Use
double quotes rather than single quotes. Expand out input validation to give
more specific error messages. Use isargout() to avoid calculating unnecessary
outputs. Add BIST tests for input validation.
* pqpnonneg.m: Ad XREF to optimset. Wrap function prototype to < 80 chars.
| author | Rik <rik@octave.org> |
|---|---|
| date | Mon, 14 Nov 2016 20:13:26 -0800 |
| parents | d261dcd73dcf |
| children | ef4d915df748 |
| rev | line source |
|---|---|
|
22323
bac0d6f07a3e
maint: Update copyright notices for 2016.
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
1 ## Copyright (C) 2008-2016 Bill Denney |
| 9635 | 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 | |
|
22755
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22323
diff
changeset
|
9 ## the Free Software Foundation; either version 3 of the License, or |
|
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22323
diff
changeset
|
10 ## (at your option) any later version. |
| 9635 | 11 ## |
| 12 ## Octave is distributed in the hope that it will be useful, but | |
| 13 ## WITHOUT ANY WARRANTY; without even the implied warranty of | |
|
22755
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22323
diff
changeset
|
14 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
|
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22323
diff
changeset
|
15 ## GNU General Public License for more details. |
| 9635 | 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 -*- | |
|
20852
516bb87ea72e
2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents:
20735
diff
changeset
|
22 ## @deftypefn {} {@var{x} =} pqpnonneg (@var{c}, @var{d}) |
|
516bb87ea72e
2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents:
20735
diff
changeset
|
23 ## @deftypefnx {} {@var{x} =} pqpnonneg (@var{c}, @var{d}, @var{x0}) |
|
22760
c4d80b9d2898
maint: Capitalize variable names appearing in error() messages of m-files.
Rik <rik@octave.org>
parents:
22755
diff
changeset
|
24 ## @deftypefnx {} {@var{x} =} pqpnonneg (@var{c}, @var{d}, @var{x0}, @var{options}) |
|
20852
516bb87ea72e
2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents:
20735
diff
changeset
|
25 ## @deftypefnx {} {[@var{x}, @var{minval}] =} pqpnonneg (@dots{}) |
|
516bb87ea72e
2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents:
20735
diff
changeset
|
26 ## @deftypefnx {} {[@var{x}, @var{minval}, @var{exitflag}] =} pqpnonneg (@dots{}) |
|
516bb87ea72e
2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents:
20735
diff
changeset
|
27 ## @deftypefnx {} {[@var{x}, @var{minval}, @var{exitflag}, @var{output}] =} pqpnonneg (@dots{}) |
|
516bb87ea72e
2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents:
20735
diff
changeset
|
28 ## @deftypefnx {} {[@var{x}, @var{minval}, @var{exitflag}, @var{output}, @var{lambda}] =} pqpnonneg (@dots{}) |
| 22769 | 29 ## |
| 22764 | 30 ## Minimize @code{1/2*@var{x}'*@var{c}*@var{x} + @var{d}'*@var{x}} subject to |
| 31 ## @code{@var{x} >= 0}. | |
|
20165
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
32 ## |
| 22764 | 33 ## @var{c} and @var{d} must be real matrices, and @var{c} must be symmetric and |
|
20165
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
34 ## positive definite. |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
35 ## |
| 22764 | 36 ## @var{x0} is an optional initial guess for the solution @var{x}. |
| 9635 | 37 ## |
|
22760
c4d80b9d2898
maint: Capitalize variable names appearing in error() messages of m-files.
Rik <rik@octave.org>
parents:
22755
diff
changeset
|
38 ## @var{options} is an options structure to change the behavior of the |
| 22769 | 39 ## algorithm (@pxref{XREFoptimset,,optimset}. @code{pqpnonneg} recognizes |
| 40 ## one option: @qcode{"MaxIter"}. | |
|
22760
c4d80b9d2898
maint: Capitalize variable names appearing in error() messages of m-files.
Rik <rik@octave.org>
parents:
22755
diff
changeset
|
41 ## |
| 9635 | 42 ## Outputs: |
|
14366
b76f0740940e
doc: Periodic grammar check of documentation.
Rik <octave@nomad.inbox5.com>
parents:
14363
diff
changeset
|
43 ## |
| 22764 | 44 ## @table @var |
| 45 ## | |
| 46 ## @item x | |
| 47 ## The solution matrix | |
| 48 ## | |
| 9635 | 49 ## @item minval |
| 22764 | 50 ## The minimum attained model value, |
| 51 ## @code{1/2*@var{xmin}'*@var{c}*@var{xmin} + @var{d}'*@var{xmin}} | |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10793
diff
changeset
|
52 ## |
| 9635 | 53 ## @item exitflag |
|
20165
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
54 ## An indicator of convergence. 0 indicates that the iteration count was |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
55 ## exceeded, and therefore convergence was not reached; >0 indicates that the |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
56 ## algorithm converged. (The algorithm is stable and will converge given |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
57 ## enough iterations.) |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10793
diff
changeset
|
58 ## |
| 9635 | 59 ## @item output |
| 60 ## A structure with two fields: | |
|
14366
b76f0740940e
doc: Periodic grammar check of documentation.
Rik <octave@nomad.inbox5.com>
parents:
14363
diff
changeset
|
61 ## |
| 9635 | 62 ## @itemize @bullet |
|
17281
bc924baa2c4e
doc: Add new @qcode macro for code samples which are quoted.
Rik <rik@octave.org>
parents:
14366
diff
changeset
|
63 ## @item @qcode{"algorithm"}: The algorithm used (@qcode{"nnls"}) |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10793
diff
changeset
|
64 ## |
|
17281
bc924baa2c4e
doc: Add new @qcode macro for code samples which are quoted.
Rik <rik@octave.org>
parents:
14366
diff
changeset
|
65 ## @item @qcode{"iterations"}: The number of iterations taken. |
| 9635 | 66 ## @end itemize |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10793
diff
changeset
|
67 ## |
| 9635 | 68 ## @item lambda |
| 22764 | 69 ## @c FIXME: Something is output from the function, but what is it? |
| 70 ## Undocumented output | |
| 71 ## @end table | |
| 72 ## @seealso{lsqnonneg, qp, optimset} | |
| 9635 | 73 ## @end deftypefn |
| 74 | |
|
13027
b9a89ca0fb75
prevent optimization functions from setting ans in workspace at startup
John W. Eaton <jwe@octave.org>
parents:
11588
diff
changeset
|
75 ## PKG_ADD: ## Discard result to avoid polluting workspace with ans at startup. |
|
b9a89ca0fb75
prevent optimization functions from setting ans in workspace at startup
John W. Eaton <jwe@octave.org>
parents:
11588
diff
changeset
|
76 ## PKG_ADD: [~] = __all_opts__ ("pqpnonneg"); |
| 9635 | 77 |
| 22764 | 78 ## This is analogical to the lsqnonneg implementation, which is implemented |
| 79 ## from Lawson and Hanson's 1973 algorithm on page 161 of Solving Least Squares | |
| 80 ## Problems. It shares the convergence guarantees. | |
| 9635 | 81 |
| 22769 | 82 function [x, minval, exitflag, output, lambda] = pqpnonneg (c, d, x0 = [], |
| 83 options = struct ()) | |
| 9635 | 84 |
| 22764 | 85 ## Special case: called to find default optimization options |
| 86 if (nargin == 1 && ischar (c) && strcmp (c, "defaults")) | |
| 9635 | 87 x = optimset ("MaxIter", 1e5); |
|
17312
088d014a7fe2
Use semicolon after "return" statement in core m-files.
Rik <rik@octave.org>
parents:
17281
diff
changeset
|
88 return; |
| 9635 | 89 endif |
| 90 | |
| 22764 | 91 if (nargin < 2 || nargin > 4) |
| 9635 | 92 print_usage (); |
| 93 endif | |
| 94 | |
| 22764 | 95 if (! (isnumeric (c) && ismatrix (c)) || ! (isnumeric (d) && ismatrix (d))) |
| 96 error ("pqpnonneg: C and D must be numeric matrices"); | |
| 97 endif | |
| 98 if (! issquare (c)) | |
| 99 error ("pqpnonneg: C must be a square matrix"); | |
| 100 endif | |
| 101 | |
| 102 if (! isstruct (options)) | |
| 103 error ("pqpnonneg: OPTIONS must be a struct"); | |
| 104 endif | |
| 105 | |
| 9635 | 106 ## Lawson-Hanson Step 1 (LH1): initialize the variables. |
| 107 n = columns (c); | |
| 22764 | 108 if (isempty (x0)) |
| 109 ## Initial guess is all zeros. | |
| 9635 | 110 x = zeros (n, 1); |
| 111 else | |
| 112 ## ensure nonnegative guess. | |
| 22764 | 113 x = max (x0, 0); |
| 9635 | 114 endif |
| 115 | |
| 116 max_iter = optimget (options, "MaxIter", 1e5); | |
| 117 | |
| 118 ## Initialize P, according to zero pattern of x. | |
| 119 p = find (x > 0).'; | |
| 120 ## Initialize the Cholesky factorization. | |
| 121 r = chol (c(p, p)); | |
| 122 usechol = true; | |
| 123 | |
| 124 iter = 0; | |
| 125 | |
| 126 ## LH3: test for completion. | |
| 127 while (iter < max_iter) | |
| 128 while (iter < max_iter) | |
|
20735
418ae0cb752f
Replace ++,-- with in-place operators for performance.
Rik <rik@octave.org>
parents:
20165
diff
changeset
|
129 iter += 1; |
| 9635 | 130 |
| 131 ## LH6: compute the positive matrix and find the min norm solution | |
| 132 ## of the positive problem. | |
| 133 if (usechol) | |
| 134 xtmp = -(r \ (r' \ d(p))); | |
| 135 else | |
| 136 xtmp = -(c(p,p) \ d(p)); | |
| 137 endif | |
| 138 idx = find (xtmp < 0); | |
| 139 | |
|
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
140 if (isempty (idx)) |
| 9635 | 141 ## LH7: tmp solution found, iterate. |
| 142 x(:) = 0; | |
| 143 x(p) = xtmp; | |
| 144 break; | |
| 145 else | |
| 146 ## LH8, LH9: find the scaling factor. | |
| 147 pidx = p(idx); | |
| 22764 | 148 sf = x(pidx) ./ (x(pidx) - xtmp(idx)); |
| 9635 | 149 alpha = min (sf); |
| 150 ## LH10: adjust X. | |
| 151 xx = zeros (n, 1); | |
| 152 xx(p) = xtmp; | |
| 153 x += alpha*(xx - x); | |
| 154 ## LH11: move from P to Z all X == 0. | |
| 155 ## This corresponds to those indices where minimum of sf is attained. | |
|
19833
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
156 idx = idx(sf == alpha); |
| 9635 | 157 p(idx) = []; |
| 158 if (usechol) | |
| 159 ## update the Cholesky factorization. | |
| 160 r = choldelete (r, idx); | |
| 161 endif | |
| 162 endif | |
| 163 endwhile | |
|
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
164 |
| 9635 | 165 ## compute the gradient. |
| 166 w = -(d + c*x); | |
| 167 w(p) = []; | |
| 168 if (! any (w > 0)) | |
| 169 if (usechol) | |
| 170 ## verify the solution achieved using qr updating. | |
| 22764 | 171 ## In the best case, this should only take a single step. |
| 9635 | 172 usechol = false; |
| 173 continue; | |
| 174 else | |
| 175 ## we're finished. | |
| 176 break; | |
| 177 endif | |
| 178 endif | |
| 179 | |
| 180 ## find the maximum gradient. | |
| 181 idx = find (w == max (w)); | |
| 182 if (numel (idx) > 1) | |
| 183 warning ("pqpnonneg:nonunique", | |
|
11588
d5bd2766c640
style fixes for warning and error messages in script files
John W. Eaton <jwe@octave.org>
parents:
11587
diff
changeset
|
184 "a non-unique solution may be returned due to equal gradients"); |
| 9635 | 185 idx = idx(1); |
| 186 endif | |
|
21751
b571fc85953f
maint: Use two spaces after period to indicate sentence break.
Rik <rik@octave.org>
parents:
21178
diff
changeset
|
187 ## move the index from Z to P. Keep P sorted. |
| 9635 | 188 z = [1:n]; z(p) = []; |
| 189 zidx = z(idx); | |
| 190 jdx = 1 + lookup (p, zidx); | |
| 191 p = [p(1:jdx-1), zidx, p(jdx:end)]; | |
| 192 if (usechol) | |
| 193 ## insert the column into the Cholesky factorization. | |
|
10200
7c1b1c084af1
handle cholesky update failure in pqpnonneg
Jaroslav Hajek <highegg@gmail.com>
parents:
9635
diff
changeset
|
194 [r, bad] = cholinsert (r, jdx, c(p,zidx)); |
|
7c1b1c084af1
handle cholesky update failure in pqpnonneg
Jaroslav Hajek <highegg@gmail.com>
parents:
9635
diff
changeset
|
195 if (bad) |
|
7c1b1c084af1
handle cholesky update failure in pqpnonneg
Jaroslav Hajek <highegg@gmail.com>
parents:
9635
diff
changeset
|
196 ## If the insertion failed, we switch off updates and go on. |
|
7c1b1c084af1
handle cholesky update failure in pqpnonneg
Jaroslav Hajek <highegg@gmail.com>
parents:
9635
diff
changeset
|
197 usechol = false; |
|
7c1b1c084af1
handle cholesky update failure in pqpnonneg
Jaroslav Hajek <highegg@gmail.com>
parents:
9635
diff
changeset
|
198 endif |
| 9635 | 199 endif |
| 200 | |
| 201 endwhile | |
| 202 ## LH12: complete. | |
| 203 | |
| 204 ## Generate the additional output arguments. | |
| 22764 | 205 if (isargout (2)) |
| 9635 | 206 minval = 1/2*(x'*c*x) + d'*x; |
| 207 endif | |
| 22764 | 208 if (isargout (3)) |
| 209 if (iter >= max_iter) | |
| 210 exitflag = 0; | |
| 211 else | |
| 212 exitflag = iter; | |
| 213 endif | |
| 9635 | 214 endif |
| 22764 | 215 if (isargout (4)) |
| 9635 | 216 output = struct ("algorithm", "nnls-pqp", "iterations", iter); |
| 217 endif | |
| 22764 | 218 if (isargout (5)) |
| 9635 | 219 lambda = zeros (size (x)); |
| 220 lambda(p) = w; | |
| 221 endif | |
| 222 | |
| 223 endfunction | |
| 224 | |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
225 |
| 9635 | 226 %!test |
| 227 %! C = [5 2;2 2]; | |
| 228 %! d = [3; -1]; | |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
229 %! assert (pqpnonneg (C, d), [0;0.5], 100*eps); |
| 9635 | 230 |
| 231 ## Test equivalence of lsq and pqp | |
| 232 %!test | |
| 233 %! C = rand (20, 10); | |
| 234 %! d = rand (20, 1); | |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
235 %! assert (pqpnonneg (C'*C, -C'*d), lsqnonneg (C, d), 100*eps); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
236 |
| 22764 | 237 # Test input validation |
| 238 %!error pqpnonneg () | |
| 239 %!error pqpnonneg (1) | |
| 240 %!error pqpnonneg (1,2,3,4,5) | |
| 241 %!error <C .* must be numeric matrices> pqpnonneg ({1},2) | |
| 242 %!error <C .* must be numeric matrices> pqpnonneg (ones (2,2,2),2) | |
| 243 %!error <D must be numeric matrices> pqpnonneg (1,{2}) | |
| 244 %!error <D must be numeric matrices> pqpnonneg (1, ones (2,2,2)) | |
| 245 %!error <C must be a square matrix> pqpnonneg ([1 2], 2) | |
| 246 %!error <OPTIONS must be a struct> pqpnonneg (1, 2, [], 3) | |
| 247 |