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