octave_packages/m/optimization/lsqnonneg.m

   1 ## Copyright (C) 2008-2012 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} =} lsqnonneg (@var{c}, @var{d})
  23 ## @deftypefnx {Function File} {@var{x} =} lsqnonneg (@var{c}, @var{d}, @var{x0})
  24 ## @deftypefnx {Function File} {[@var{x}, @var{resnorm}] =} lsqnonneg (@dots{})
  25 ## @deftypefnx {Function File} {[@var{x}, @var{resnorm}, @var{residual}] =} lsqnonneg (@dots{})
  26 ## @deftypefnx {Function File} {[@var{x}, @var{resnorm}, @var{residual}, @var{exitflag}] =} lsqnonneg (@dots{})
  27 ## @deftypefnx {Function File} {[@var{x}, @var{resnorm}, @var{residual}, @var{exitflag}, @var{output}] =} lsqnonneg (@dots{})
  28 ## @deftypefnx {Function File} {[@var{x}, @var{resnorm}, @var{residual}, @var{exitflag}, @var{output}, @var{lambda}] =} lsqnonneg (@dots{})
  29 ## Minimize @code{norm (@var{c}*@var{x} - d)} subject to
  30 ## @code{@var{x} >= 0}.  @var{c} and @var{d} must be real.  @var{x0} is an
  31 ## optional initial guess for @var{x}.
  32 ##
  33 ## Outputs:
  34 ## @itemize @bullet
  35 ## @item resnorm
  36 ##
  37 ## The squared 2-norm of the residual: norm(@var{c}*@var{x}-@var{d})^2
  38 ##
  39 ## @item residual
  40 ##
  41 ## The residual: @var{d}-@var{c}*@var{x}
  42 ##
  43 ## @item exitflag
  44 ##
  45 ## An indicator of convergence.  0 indicates that the iteration count
  46 ## was exceeded, and therefore convergence was not reached; >0 indicates
  47 ## that the algorithm converged.  (The algorithm is stable and will
  48 ## converge given enough iterations.)
  49 ##
  50 ## @item output
  51 ##
  52 ## A structure with two fields:
  53 ## @itemize @bullet
  54 ## @item "algorithm": The algorithm used ("nnls")
  55 ##
  56 ## @item "iterations": The number of iterations taken.
  57 ## @end itemize
  58 ##
  59 ## @item lambda
  60 ##
  61 ## Not implemented.
  62 ## @end itemize
  63 ## @seealso{optimset, pqpnonneg}
  64 ## @end deftypefn
  65
  66 ## PKG_ADD: ## Discard result to avoid polluting workspace with ans at startup.
  67 ## PKG_ADD: [~] = __all_opts__ ("lsqnonneg");
  68
  69 ## This is implemented from Lawson and Hanson's 1973 algorithm on page
  70 ## 161 of Solving Least Squares Problems.
  71
  72 function [x, resnorm, residual, exitflag, output, lambda] = lsqnonneg (c, d, x = [], options = struct ())
  73
  74   if (nargin == 1 && ischar (c) && strcmp (c, 'defaults'))
  75     x = optimset ("MaxIter", 1e5);
  76     return
  77   endif
  78
  79   if (! (nargin >= 2 && nargin <= 4 && ismatrix (c) && ismatrix (d) && isstruct (options)))
  80     print_usage ();
  81   endif
  82
  83   ## Lawson-Hanson Step 1 (LH1): initialize the variables.
  84   m = rows (c);
  85   n = columns (c);
  86   if (isempty (x))
  87     ## Initial guess is 0s.
  88     x = zeros (n, 1);
  89   else
  90     ## ensure nonnegative guess.
  91     x = max (x, 0);
  92   endif
  93
  94   useqr = m >= n;
  95   max_iter = optimget (options, "MaxIter", 1e5);
  96
  97   ## Initialize P, according to zero pattern of x.
  98   p = find (x > 0).';
  99   if (useqr)
 100     ## Initialize the QR factorization, economized form.
 101     [q, r] = qr (c(:,p), 0);
 102   endif
 103
 104   iter = 0;
 105
 106   ## LH3: test for completion.
 107   while (iter < max_iter)
 108     while (iter < max_iter)
 109       iter++;
 110
 111       ## LH6: compute the positive matrix and find the min norm solution
 112       ## of the positive problem.
 113       if (useqr)
 114         xtmp = r \ q'*d;
 115       else
 116         xtmp = c(:,p) \ d;
 117       endif
 118       idx = find (xtmp < 0);
 119
 120       if (isempty (idx))
 121         ## LH7: tmp solution found, iterate.
 122         x(:) = 0;
 123         x(p) = xtmp;
 124         break;
 125       else
 126         ## LH8, LH9: find the scaling factor.
 127         pidx = p(idx);
 128         sf = x(pidx)./(x(pidx) - xtmp(idx));
 129         alpha = min (sf);
 130         ## LH10: adjust X.
 131         xx = zeros (n, 1);
 132         xx(p) = xtmp;
 133         x += alpha*(xx - x);
 134         ## LH11: move from P to Z all X == 0.
 135         ## This corresponds to those indices where minimum of sf is attained.
 136         idx = idx (sf == alpha);
 137         p(idx) = [];
 138         if (useqr)
 139           ## update the QR factorization.
 140           [q, r] = qrdelete (q, r, idx);
 141         endif
 142       endif
 143     endwhile
 144
 145     ## compute the gradient.
 146     w = c'*(d - c*x);
 147     w(p) = [];
 148     if (! any (w > 0))
 149       if (useqr)
 150         ## verify the solution achieved using qr updating.
 151         ## in the best case, this should only take a single step.
 152         useqr = false;
 153         continue;
 154       else
 155         ## we're finished.
 156         break;
 157       endif
 158     endif
 159
 160     ## find the maximum gradient.
 161     idx = find (w == max (w));
 162     if (numel (idx) > 1)
 163       warning ("lsqnonneg:nonunique",
 164                "a non-unique solution may be returned due to equal gradients");
 165       idx = idx(1);
 166     endif
 167     ## move the index from Z to P. Keep P sorted.
 168     z = [1:n]; z(p) = [];
 169     zidx = z(idx);
 170     jdx = 1 + lookup (p, zidx);
 171     p = [p(1:jdx-1), zidx, p(jdx:end)];
 172     if (useqr)
 173       ## insert the column into the QR factorization.
 174       [q, r] = qrinsert (q, r, jdx, c(:,zidx));
 175     endif
 176
 177   endwhile
 178   ## LH12: complete.
 179
 180   ## Generate the additional output arguments.
 181   if (nargout > 1)
 182     resnorm = norm (c*x - d) ^ 2;
 183   endif
 184   if (nargout > 2)
 185     residual = d - c*x;
 186   endif
 187   exitflag = iter;
 188   if (nargout > 3 && iter >= max_iter)
 189     exitflag = 0;
 190   endif
 191   if (nargout > 4)
 192     output = struct ("algorithm", "nnls", "iterations", iter);
 193   endif
 194   if (nargout > 5)
 195     lambda = zeros (size (x));
 196     lambda(p) = w;
 197   endif
 198
 199 endfunction
 200
 201 ## Tests
 202 %!test
 203 %! C = [1 0;0 1;2 1];
 204 %! d = [1;3;-2];
 205 %! assert (lsqnonneg (C, d), [0;0.5], 100*eps)
 206
 207 %!test
 208 %! C = [0.0372 0.2869;0.6861 0.7071;0.6233 0.6245;0.6344 0.6170];
 209 %! d = [0.8587;0.1781;0.0747;0.8405];
 210 %! xnew = [0;0.6929];
 211 %! assert (lsqnonneg (C, d), xnew, 0.0001)