Mercurial > octave
annotate scripts/optimization/pqpnonneg.m @ 30564:796f54d4ddbf stable
update Octave Project Developers copyright for the new year
In files that have the "Octave Project Developers" copyright notice,
update for 2021.
In all .txi and .texi files except gpl.txi and gpl.texi in the
doc/liboctave and doc/interpreter directories, change the copyright
to "Octave Project Developers", the same as used for other source
files. Update copyright notices for 2022 (not done since 2019). For
gpl.txi and gpl.texi, change the copyright notice to be "Free Software
Foundation, Inc." and leave the date at 2007 only because this file
only contains the text of the GPL, not anything created by the Octave
Project Developers.
Add Paul Thomas to contributors.in.
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Tue, 28 Dec 2021 18:22:40 -0500 |
| parents | 7854d5752dd2 |
| children | e5e35410da90 |
| rev | line source |
|---|---|
|
27923
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
1 ######################################################################## |
|
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
2 ## |
|
30564
796f54d4ddbf
update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents:
29359
diff
changeset
|
3 ## Copyright (C) 2008-2022 The Octave Project Developers |
|
27918
b442ec6dda5c
use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents:
26376
diff
changeset
|
4 ## |
|
27923
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
5 ## See the file COPYRIGHT.md in the top-level directory of this |
|
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
6 ## distribution or <https://octave.org/copyright/>. |
| 9635 | 7 ## |
| 8 ## This file is part of Octave. | |
| 9 ## | |
|
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
10 ## Octave is free software: you can redistribute it and/or modify it |
| 9635 | 11 ## under the terms of the GNU General Public License as published by |
|
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
12 ## the Free Software Foundation, either version 3 of the License, or |
|
22755
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22323
diff
changeset
|
13 ## (at your option) any later version. |
| 9635 | 14 ## |
| 15 ## Octave is distributed in the hope that it will be useful, but | |
| 16 ## WITHOUT ANY WARRANTY; without even the implied warranty of | |
|
22755
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22323
diff
changeset
|
17 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
|
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22323
diff
changeset
|
18 ## GNU General Public License for more details. |
| 9635 | 19 ## |
| 20 ## You should have received a copy of the GNU General Public License | |
| 21 ## along with Octave; see the file COPYING. If not, see | |
|
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
22 ## <https://www.gnu.org/licenses/>. |
|
27923
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
23 ## |
|
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
24 ######################################################################## |
| 9635 | 25 |
| 26 ## -*- texinfo -*- | |
|
20852
516bb87ea72e
2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents:
20735
diff
changeset
|
27 ## @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
|
28 ## @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
|
29 ## @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
|
30 ## @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
|
31 ## @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
|
32 ## @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
|
33 ## @deftypefnx {} {[@var{x}, @var{minval}, @var{exitflag}, @var{output}, @var{lambda}] =} pqpnonneg (@dots{}) |
|
20165
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
34 ## |
| 22764 | 35 ## Minimize @code{1/2*@var{x}'*@var{c}*@var{x} + @var{d}'*@var{x}} subject to |
| 36 ## @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
|
37 ## |
| 22764 | 38 ## @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
|
39 ## positive definite. |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
40 ## |
| 22764 | 41 ## @var{x0} is an optional initial guess for the solution @var{x}. |
| 9635 | 42 ## |
|
22760
c4d80b9d2898
maint: Capitalize variable names appearing in error() messages of m-files.
Rik <rik@octave.org>
parents:
22755
diff
changeset
|
43 ## @var{options} is an options structure to change the behavior of the |
|
28959
5394d688d456
doc: Use @code{} within alternate text for @xref,@pxref macros for better Info display.
Rik <rik@octave.org>
parents:
28942
diff
changeset
|
44 ## algorithm (@pxref{XREFoptimset,,@code{optimset}}). @code{pqpnonneg} |
|
5394d688d456
doc: Use @code{} within alternate text for @xref,@pxref macros for better Info display.
Rik <rik@octave.org>
parents:
28942
diff
changeset
|
45 ## recognizes one option: @qcode{"MaxIter"}. |
| 9635 | 46 ## |
| 47 ## Outputs: | |
|
14366
b76f0740940e
doc: Periodic grammar check of documentation.
Rik <octave@nomad.inbox5.com>
parents:
14363
diff
changeset
|
48 ## |
| 22764 | 49 ## @table @var |
| 50 ## | |
| 51 ## @item x | |
| 52 ## The solution matrix | |
| 53 ## | |
| 9635 | 54 ## @item minval |
| 22764 | 55 ## The minimum attained model value, |
| 56 ## @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
|
57 ## |
| 9635 | 58 ## @item exitflag |
|
20165
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
59 ## 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
|
60 ## 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
|
61 ## 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
|
62 ## enough iterations.) |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10793
diff
changeset
|
63 ## |
| 9635 | 64 ## @item output |
| 65 ## A structure with two fields: | |
|
14366
b76f0740940e
doc: Periodic grammar check of documentation.
Rik <octave@nomad.inbox5.com>
parents:
14363
diff
changeset
|
66 ## |
| 9635 | 67 ## @itemize @bullet |
|
25003
2365c2661b3c
doc: Spellcheck documentation ahead of 4.4 release.
Rik <rik@octave.org>
parents:
24534
diff
changeset
|
68 ## @item @qcode{"algorithm"}: The algorithm used (@nospell{@qcode{"nnls"}}) |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10793
diff
changeset
|
69 ## |
|
17281
bc924baa2c4e
doc: Add new @qcode macro for code samples which are quoted.
Rik <rik@octave.org>
parents:
14366
diff
changeset
|
70 ## @item @qcode{"iterations"}: The number of iterations taken. |
| 9635 | 71 ## @end itemize |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10793
diff
changeset
|
72 ## |
| 9635 | 73 ## @item lambda |
| 22764 | 74 ## @c FIXME: Something is output from the function, but what is it? |
| 75 ## Undocumented output | |
| 76 ## @end table | |
| 77 ## @seealso{lsqnonneg, qp, optimset} | |
| 9635 | 78 ## @end deftypefn |
| 79 | |
|
13027
b9a89ca0fb75
prevent optimization functions from setting ans in workspace at startup
John W. Eaton <jwe@octave.org>
parents:
11588
diff
changeset
|
80 ## 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
|
81 ## PKG_ADD: [~] = __all_opts__ ("pqpnonneg"); |
| 9635 | 82 |
| 22764 | 83 ## This is analogical to the lsqnonneg implementation, which is implemented |
| 84 ## from Lawson and Hanson's 1973 algorithm on page 161 of Solving Least Squares | |
| 85 ## Problems. It shares the convergence guarantees. | |
| 9635 | 86 |
| 22769 | 87 function [x, minval, exitflag, output, lambda] = pqpnonneg (c, d, x0 = [], |
| 88 options = struct ()) | |
| 9635 | 89 |
| 22764 | 90 ## Special case: called to find default optimization options |
| 91 if (nargin == 1 && ischar (c) && strcmp (c, "defaults")) | |
|
26138
804e18e3e320
Reenable query of optimization options (bugs #54952 and #55089).
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
25579
diff
changeset
|
92 x = struct ("MaxIter", 1e5); |
|
17312
088d014a7fe2
Use semicolon after "return" statement in core m-files.
Rik <rik@octave.org>
parents:
17281
diff
changeset
|
93 return; |
| 9635 | 94 endif |
| 95 | |
|
28789
28de41192f3c
Eliminate unneeded verification of nargin, nargout in m-files.
Rik <rik@octave.org>
parents:
27923
diff
changeset
|
96 if (nargin < 2) |
| 9635 | 97 print_usage (); |
| 98 endif | |
| 99 | |
| 22764 | 100 if (! (isnumeric (c) && ismatrix (c)) || ! (isnumeric (d) && ismatrix (d))) |
| 101 error ("pqpnonneg: C and D must be numeric matrices"); | |
| 102 endif | |
| 103 if (! issquare (c)) | |
| 104 error ("pqpnonneg: C must be a square matrix"); | |
| 105 endif | |
| 106 | |
| 107 if (! isstruct (options)) | |
| 108 error ("pqpnonneg: OPTIONS must be a struct"); | |
| 109 endif | |
| 110 | |
| 9635 | 111 ## Lawson-Hanson Step 1 (LH1): initialize the variables. |
| 112 n = columns (c); | |
| 22764 | 113 if (isempty (x0)) |
| 114 ## Initial guess is all zeros. | |
| 9635 | 115 x = zeros (n, 1); |
| 116 else | |
| 117 ## ensure nonnegative guess. | |
| 22764 | 118 x = max (x0, 0); |
| 9635 | 119 endif |
| 120 | |
| 121 max_iter = optimget (options, "MaxIter", 1e5); | |
| 122 | |
| 123 ## Initialize P, according to zero pattern of x. | |
| 124 p = find (x > 0).'; | |
| 125 ## Initialize the Cholesky factorization. | |
| 126 r = chol (c(p, p)); | |
| 127 usechol = true; | |
| 128 | |
| 129 iter = 0; | |
| 130 | |
| 131 ## LH3: test for completion. | |
| 132 while (iter < max_iter) | |
| 133 while (iter < max_iter) | |
|
20735
418ae0cb752f
Replace ++,-- with in-place operators for performance.
Rik <rik@octave.org>
parents:
20165
diff
changeset
|
134 iter += 1; |
| 9635 | 135 |
| 136 ## LH6: compute the positive matrix and find the min norm solution | |
| 137 ## of the positive problem. | |
| 138 if (usechol) | |
| 139 xtmp = -(r \ (r' \ d(p))); | |
| 140 else | |
| 141 xtmp = -(c(p,p) \ d(p)); | |
| 142 endif | |
| 143 idx = find (xtmp < 0); | |
| 144 | |
|
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
145 if (isempty (idx)) |
| 9635 | 146 ## LH7: tmp solution found, iterate. |
| 147 x(:) = 0; | |
| 148 x(p) = xtmp; | |
| 149 break; | |
| 150 else | |
| 151 ## LH8, LH9: find the scaling factor. | |
| 152 pidx = p(idx); | |
| 22764 | 153 sf = x(pidx) ./ (x(pidx) - xtmp(idx)); |
| 9635 | 154 alpha = min (sf); |
| 155 ## LH10: adjust X. | |
| 156 xx = zeros (n, 1); | |
| 157 xx(p) = xtmp; | |
| 158 x += alpha*(xx - x); | |
| 159 ## LH11: move from P to Z all X == 0. | |
| 160 ## 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
|
161 idx = idx(sf == alpha); |
| 9635 | 162 p(idx) = []; |
| 163 if (usechol) | |
| 164 ## update the Cholesky factorization. | |
| 165 r = choldelete (r, idx); | |
| 166 endif | |
| 167 endif | |
| 168 endwhile | |
|
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
169 |
| 9635 | 170 ## compute the gradient. |
| 171 w = -(d + c*x); | |
| 172 w(p) = []; | |
| 173 if (! any (w > 0)) | |
| 174 if (usechol) | |
| 175 ## verify the solution achieved using qr updating. | |
| 22764 | 176 ## In the best case, this should only take a single step. |
| 9635 | 177 usechol = false; |
| 178 continue; | |
| 179 else | |
| 180 ## we're finished. | |
| 181 break; | |
| 182 endif | |
| 183 endif | |
| 184 | |
| 185 ## find the maximum gradient. | |
| 186 idx = find (w == max (w)); | |
| 187 if (numel (idx) > 1) | |
| 188 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
|
189 "a non-unique solution may be returned due to equal gradients"); |
| 9635 | 190 idx = idx(1); |
| 191 endif | |
|
21751
b571fc85953f
maint: Use two spaces after period to indicate sentence break.
Rik <rik@octave.org>
parents:
21178
diff
changeset
|
192 ## move the index from Z to P. Keep P sorted. |
| 9635 | 193 z = [1:n]; z(p) = []; |
| 194 zidx = z(idx); | |
| 195 jdx = 1 + lookup (p, zidx); | |
| 196 p = [p(1:jdx-1), zidx, p(jdx:end)]; | |
| 197 if (usechol) | |
| 198 ## insert the column into the Cholesky factorization. | |
|
10200
7c1b1c084af1
handle cholesky update failure in pqpnonneg
Jaroslav Hajek <highegg@gmail.com>
parents:
9635
diff
changeset
|
199 [r, bad] = cholinsert (r, jdx, c(p,zidx)); |
|
7c1b1c084af1
handle cholesky update failure in pqpnonneg
Jaroslav Hajek <highegg@gmail.com>
parents:
9635
diff
changeset
|
200 if (bad) |
|
7c1b1c084af1
handle cholesky update failure in pqpnonneg
Jaroslav Hajek <highegg@gmail.com>
parents:
9635
diff
changeset
|
201 ## 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
|
202 usechol = false; |
|
7c1b1c084af1
handle cholesky update failure in pqpnonneg
Jaroslav Hajek <highegg@gmail.com>
parents:
9635
diff
changeset
|
203 endif |
| 9635 | 204 endif |
| 205 | |
| 206 endwhile | |
| 207 ## LH12: complete. | |
| 208 | |
| 209 ## Generate the additional output arguments. | |
| 22764 | 210 if (isargout (2)) |
| 9635 | 211 minval = 1/2*(x'*c*x) + d'*x; |
| 212 endif | |
| 22764 | 213 if (isargout (3)) |
| 214 if (iter >= max_iter) | |
| 215 exitflag = 0; | |
| 216 else | |
| 217 exitflag = iter; | |
| 218 endif | |
| 9635 | 219 endif |
| 22764 | 220 if (isargout (4)) |
| 9635 | 221 output = struct ("algorithm", "nnls-pqp", "iterations", iter); |
| 222 endif | |
| 22764 | 223 if (isargout (5)) |
| 9635 | 224 lambda = zeros (size (x)); |
| 225 lambda(p) = w; | |
| 226 endif | |
| 227 | |
| 228 endfunction | |
| 229 | |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
230 |
| 9635 | 231 %!test |
| 232 %! C = [5 2;2 2]; | |
| 233 %! d = [3; -1]; | |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
234 %! assert (pqpnonneg (C, d), [0;0.5], 100*eps); |
| 9635 | 235 |
| 236 ## Test equivalence of lsq and pqp | |
| 237 %!test | |
| 238 %! C = rand (20, 10); | |
| 239 %! 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
|
240 %! 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
|
241 |
|
28942
fc4bb4bd1d5e
maint: Use '##' as lead-in for full-line comments.
Rik <rik@octave.org>
parents:
28896
diff
changeset
|
242 ## Test input validation |
|
28896
90fea9cc9caa
test: Add expected error message <Invalid call> to BIST tests for nargin.
Rik <rik@octave.org>
parents:
28789
diff
changeset
|
243 %!error <Invalid call> pqpnonneg () |
|
90fea9cc9caa
test: Add expected error message <Invalid call> to BIST tests for nargin.
Rik <rik@octave.org>
parents:
28789
diff
changeset
|
244 %!error <Invalid call> pqpnonneg (1) |
| 22764 | 245 %!error <C .* must be numeric matrices> pqpnonneg ({1},2) |
| 246 %!error <C .* must be numeric matrices> pqpnonneg (ones (2,2,2),2) | |
| 247 %!error <D must be numeric matrices> pqpnonneg (1,{2}) | |
| 248 %!error <D must be numeric matrices> pqpnonneg (1, ones (2,2,2)) | |
| 249 %!error <C must be a square matrix> pqpnonneg ([1 2], 2) | |
| 250 %!error <OPTIONS must be a struct> pqpnonneg (1, 2, [], 3) |