# HG changeset patch # User Jaroslav Hajek # Date 1252660598 -7200 # Node ID 36d885c4a1ac230c85138e732495b6da97f31928 # Parent da5ba66414a3cd3c55c651221be8efbdd0ff7ff5 implement pqpnonneg diff -r da5ba66414a3 -r 36d885c4a1ac scripts/optimization/Makefile.in --- a/scripts/optimization/Makefile.in Fri Sep 11 09:27:58 2009 +0200 +++ b/scripts/optimization/Makefile.in Fri Sep 11 11:16:38 2009 +0200 @@ -42,6 +42,7 @@ glpk.m \ glpkmex.m \ lsqnonneg.m \ + pqpnonneg.m \ optimset.m \ optimget.m \ __all_opts__.m \ diff -r da5ba66414a3 -r 36d885c4a1ac scripts/optimization/lsqnonneg.m --- a/scripts/optimization/lsqnonneg.m Fri Sep 11 09:27:58 2009 +0200 +++ b/scripts/optimization/lsqnonneg.m Fri Sep 11 11:16:38 2009 +0200 @@ -55,7 +55,7 @@ ## ## Not implemented. ## @end itemize -## @seealso{optimset} +## @seealso{optimset, pqpnonneg} ## @end deftypefn ## PKG_ADD: __all_opts__ ("lsqnonneg"); diff -r da5ba66414a3 -r 36d885c4a1ac scripts/optimization/pqpnonneg.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/optimization/pqpnonneg.m Fri Sep 11 11:16:38 2009 +0200 @@ -0,0 +1,202 @@ +## Copyright (C) 2008 Bill Denney +## Copyright (C) 2008 Jaroslav Hajek +## Copyright (C) 2009 VZLU Prague +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## . + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{x} =} pqpnonneg (@var{c}, @var{d}) +## @deftypefnx {Function File} {@var{x} =} pqpnonneg (@var{c}, @var{d}, @var{x0}) +## @deftypefnx {Function File} {[@var{x}, @var{minval}] =} pqpnonneg (@dots{}) +## @deftypefnx {Function File} {[@var{x}, @var{minval}, @var{exitflag}] =} pqpnonneg (@dots{}) +## @deftypefnx {Function File} {[@var{x}, @var{minval}, @var{exitflag}, @var{output}] =} pqpnonneg (@dots{}) +## @deftypefnx {Function File} {[@var{x}, @var{minval}, @var{exitflag}, @var{output}, @var{lambda}] =} pqpnonneg (@dots{}) +## Minimize @code{1/2*x'*c*x + d'*x} subject to @code{@var{x} >= +## 0}. @var{c} and @var{d} must be real, and @var{c} must be symmetric and positive definite. +## @var{x0} is an optional initial guess for @var{x}. +## +## Outputs: +## @itemize @bullet +## @item minval +## +## The minimum attained model value, 1/2*xmin'*c*xmin + d'*xmin +## @item exitflag +## +## An indicator of convergence. 0 indicates that the iteration count +## was exceeded, and therefore convergence was not reached; >0 indicates +## that the algorithm converged. (The algorithm is stable and will +## converge given enough iterations.) +## @item output +## +## A structure with two fields: +## @itemize @bullet +## @item "algorithm": The algorithm used ("nnls") +## @item "iterations": The number of iterations taken. +## @end itemize +## @item lambda +## +## Not implemented. +## @end itemize +## @seealso{optimset, lsqnonneg, qp} +## @end deftypefn + +## PKG_ADD: __all_opts__ ("pqpnonneg"); + +## This is analogical to the lsqnonneg implementation, which is +## implemented from Lawson and Hanson's 1973 algorithm on page +## 161 of Solving Least Squares Problems. +## It shares the convergence guarantees. + +function [x, minval, exitflag, output, lambda] = pqpnonneg (c, d, x = [], options = struct ()) + + if (nargin == 1 && ischar (c) && strcmp (c, 'defaults')) + x = optimset ("MaxIter", 1e5); + return + endif + + if (! (nargin >= 2 && nargin <= 4 && ismatrix (c) && ismatrix (d) && isstruct (options))) + print_usage (); + endif + + ## Lawson-Hanson Step 1 (LH1): initialize the variables. + m = rows (c); + n = columns (c); + if (m != n) + error ("matrix must be square"); + endif + + if (isempty (x)) + ## Initial guess is 0s. + x = zeros (n, 1); + else + ## ensure nonnegative guess. + x = max (x, 0); + endif + + max_iter = optimget (options, "MaxIter", 1e5); + + ## Initialize P, according to zero pattern of x. + p = find (x > 0).'; + ## Initialize the Cholesky factorization. + r = chol (c(p, p)); + usechol = true; + + iter = 0; + + ## LH3: test for completion. + while (iter < max_iter) + while (iter < max_iter) + iter++; + + ## LH6: compute the positive matrix and find the min norm solution + ## of the positive problem. + if (usechol) + xtmp = -(r \ (r' \ d(p))); + else + xtmp = -(c(p,p) \ d(p)); + endif + idx = find (xtmp < 0); + + if (isempty (idx)) + ## LH7: tmp solution found, iterate. + x(:) = 0; + x(p) = xtmp; + break; + else + ## LH8, LH9: find the scaling factor. + pidx = p(idx); + sf = x(pidx)./(x(pidx) - xtmp(idx)); + alpha = min (sf); + ## LH10: adjust X. + xx = zeros (n, 1); + xx(p) = xtmp; + x += alpha*(xx - x); + ## LH11: move from P to Z all X == 0. + ## This corresponds to those indices where minimum of sf is attained. + idx = idx (sf == alpha); + p(idx) = []; + if (usechol) + ## update the Cholesky factorization. + r = choldelete (r, idx); + endif + endif + endwhile + + ## compute the gradient. + w = -(d + c*x); + w(p) = []; + if (! any (w > 0)) + if (usechol) + ## verify the solution achieved using qr updating. + ## in the best case, this should only take a single step. + usechol = false; + continue; + else + ## we're finished. + break; + endif + endif + + ## find the maximum gradient. + idx = find (w == max (w)); + if (numel (idx) > 1) + warning ("pqpnonneg:nonunique", + "A non-unique solution may be returned due to equal gradients."); + idx = idx(1); + endif + ## move the index from Z to P. Keep P sorted. + z = [1:n]; z(p) = []; + zidx = z(idx); + jdx = 1 + lookup (p, zidx); + p = [p(1:jdx-1), zidx, p(jdx:end)]; + if (usechol) + ## insert the column into the Cholesky factorization. + r = cholinsert (r, jdx, c(p,zidx)); + endif + + endwhile + ## LH12: complete. + + ## Generate the additional output arguments. + if (nargout > 1) + minval = 1/2*(x'*c*x) + d'*x; + endif + exitflag = iter; + if (nargout > 2 && iter >= max_iter) + exitflag = 0; + endif + if (nargout > 3) + output = struct ("algorithm", "nnls-pqp", "iterations", iter); + endif + if (nargout > 4) + lambda = zeros (size (x)); + lambda(p) = w; + endif + +endfunction + +## Tests +%!test +%! C = [5 2;2 2]; +%! d = [3; -1]; +%! assert (pqpnonneg (C, d), [0;0.5], 100*eps) + +## Test equivalence of lsq and pqp +%!test +%! C = rand (20, 10); +%! d = rand (20, 1); +%! assert (pqpnonneg (C'*C, -C'*d), lsqnonneg (C, d), 100*eps)