Mercurial > hg > octave-lojdl
changeset 7682:f7474c83c782
lsqnonneg: new function
author | bill@denney.ws |
---|---|
date | Tue, 01 Apr 2008 00:29:51 -0400 |
parents | b1c1133641ee |
children | 8136cb19fb7a |
files | scripts/ChangeLog scripts/optimization/Makefile.in scripts/optimization/lsqnonneg.m |
diffstat | 3 files changed, 161 insertions(+), 1 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,8 @@ +2008-04-02 Bill Denney <bill@denney.ws> + + * optimization/lsqnonneg.m: New function. + * optimization/Makefile.in (SOURCES): Add it to the list. + 2008-04-02 David Bateman <dbateman@free.fr> * sparse/spaugment.m: New function
--- a/scripts/optimization/Makefile.in +++ b/scripts/optimization/Makefile.in @@ -32,7 +32,14 @@ INSTALL_PROGRAM = @INSTALL_PROGRAM@ INSTALL_DATA = @INSTALL_DATA@ -SOURCES = __fsolve_defopts__.m glpk.m glpkmex.m optimset.m qp.m sqp.m +SOURCES = \ + __fsolve_defopts__.m \ + glpk.m \ + glpkmex.m \ + lsqnonneg.m \ + optimset.m \ + qp.m \ + sqp.m EXTRAS = glpktest1 glpktest2
new file mode 100644 --- /dev/null +++ b/scripts/optimization/lsqnonneg.m @@ -0,0 +1,148 @@ +## Copyright (C) 2008 Bill Denney +## +## 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 +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{x} =} lsqnonneg (@var{c}, @var{d}) +## @deftypefnx {Function File} {@var{x} =} lsqnonneg (@var{c}, @var{d}, @var{x0}) +## @deftypefnx {Function File} {[@var{x}, @var{resnorm}] =} lsqnonneg (@dots{}) +## @deftypefnx {Function File} {[@var{x}, @var{resnorm}, @var{residual}] =} lsqnonneg (@dots{}) +## @deftypefnx {Function File} {[@var{x}, @var{resnorm}, @var{residual}, @var{exitflag}] =} lsqnonneg (@dots{}) +## @deftypefnx {Function File} {[@var{x}, @var{resnorm}, @var{residual}, @var{exitflag}, @var{output}] =} lsqnonneg (@dots{}) +## @deftypefnx {Function File} {[@var{x}, @var{resnorm}, @var{residual}, @var{exitflag}, @var{output}, @var{lambda}] =} lsqnonneg (@dots{}) +## Minimize @code{norm (@var{c}*@var{x}-d)} subject to @code{@var{x} >= +## 0}. @var{c} and @var{d} must be real. @var{x0} is an optional +## initial guess for @var{x}. +## +## Outputs: +## @itemize @bullet +## @item resnorm +## +## The squared 2-norm of the residual: norm(@var{c}*@var{x}-@var{d})^2 +## @item residual +## +## The residual: @var{d}-@var{c}*@var{x} +## @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} +## @end deftypefn + +## This is implemented from Lawson and Hanson's 1973 algorithm on page + +function [x, resnorm, residual, exitflag, output, lambda] = lsqnonneg (c, d, x = [], options = []) + + if (isempty (x)) + ## initial guess is 0s + x = zeros (columns (c), 1); + endif + + if (isempty (options)) + ## FIXME: Optimset should be updated + ## options = optimset (); + options = struct ("maxiter", 1e5, "tolx", 1e-8); + endif + + ## Initialize the values + p = []; + z = 1:numel (x); + ## compute the gradient + w = c'*(d - c*x); + + iter = 0; + while (! isempty (z) && any (w(z) > 0) && iter < options.maxiter) + idx = find (w == max (w)); + if (numel (idx) > 1) + warning ("lsqnonneg:nonunique", + "A non-unique solution may be returned due to equal gradients."); + idx = idx(1); + endif + p(end+1) = z(idx); + z(idx) = []; + + newx = false; + while (! newx && iter < options.maxiter) + iter++; + + cpos = c; + cpos(:,z) = 0; + ## find min norm solution to the positive matrix + [cpos_q, cpos_r] = qr (cpos, 0); + xtmp = (cpos_r\cpos_q')*d; + xtmp(z) = 0; + if (all (xtmp >= 0)) + ## complete the inner loop + newx = true; + x = xtmp; + else + mask = (xtmp < 0); + alpha = min(x(mask)./(x(mask) - xtmp(mask))); + x = x + alpha*(xtmp - x); + idx = find (x == 0); + p(ismember (p, x)) = []; + z = [z idx]; + endif + endwhile + w = c'*(d - c*x); + endwhile + + ## generate the additional output arguments + if (nargout > 1) + resnorm = norm (C*x-d) ^ 2; + endif + if (nargout > 2) + residual = d-C*x; + endif + exitflag = iter; + if (nargout > 3 && iter >= options.maxiter) + exitflag = 0; + endif + if (nargout > 4) + output = struct ("algorithm", "nnls", "iterations", iter); + endif + if (nargout > 5) + ## FIXME + error ("lsqnonneg: lambda is not yet implemented"); + endif + +endfunction + +## Tests +%!test +%! C = [1 0;0 1;2 1]; +%! d = [1;3;-2]; +%! assert (lsqnonneg (C, d), [0;0.5], 100*eps) + +%!test +%! C = [0.0372 0.2869;0.6861 0.7071;0.6233 0.6245;0.6344 0.6170]; +%! d = [0.8587;0.1781;0.0747;0.8405]; +%! xnew = [0;0.6929]; +%! assert (lsqnonneg (C, d), xnew, 0.0001)