Re: New function: lsqnonlin

[Bill Denney]

John W. Eaton wrote:
> On  1-Apr-2008, Bill Denney wrote:
>
> | Attached is lsqnonlin.
>
> I applied this changeset.  Please remember to add new functions to the
> SOURCES list in the correspoding Makefile.in file.
>
> | I will probably add optimset soon.
>
> There's already a version of optimset in the current sources, though
> it may need some work.
>
> Thanks,
>
> jwe
I actually just noticed a few bugs in the function that I submitted.  
Attached are the changes.  If someone has a source of additional tests, 
I would appreciate it.

Thanks,

Bill

# HG changeset patch
# User address@hidden
# Date 1207194954 14400
# Node ID a4c8a5b154e6a1f8448e7cc2f938ef6cfd158f53
# Parent  f90494008de838634f36f21386d5648daede443f
[mq]: lsqnonneg

diff -r f90494008de8 -r a4c8a5b154e6 scripts/ChangeLog
--- a/scripts/ChangeLog Wed Apr 02 23:24:33 2008 -0400
+++ b/scripts/ChangeLog Wed Apr 02 23:55:54 2008 -0400
@@ -1,3 +1,7 @@ 2008-04-02  John W. Eaton  <address@hidden
+2008-04-02  Bill Denney  <address@hidden>
+
+       * optimization/lsqnonneg.m: fix several bugs
+
 2008-04-02  John W. Eaton  <address@hidden>
 
        * deprecated/Makefile.in (SOURCES): Add spkron.m to the list.
diff -r f90494008de8 -r a4c8a5b154e6 scripts/optimization/lsqnonneg.m
--- a/scripts/optimization/lsqnonneg.m  Wed Apr 02 23:24:33 2008 -0400
+++ b/scripts/optimization/lsqnonneg.m  Wed Apr 02 23:55:54 2008 -0400
@@ -56,42 +56,51 @@
 ## @seealso{optimset}
 ## @end deftypefn
 
-## This is implemented from Lawson and Hanson's 1973 algorithm on page 
+## This is implemented from Lawson and Hanson's 1973 algorithm on page
+## 161 of Solving Least Squares Problems
+## 
(http://books.google.com/books?id=ROw4hU85nz8C&printsec=frontcover&dq=nonnegative+least+squares+algorithm#PPA161,M1)
 
 function [x, resnorm, residual, exitflag, output, lambda] = lsqnonneg (c, d, x 
= [], options = [])
 
+  ## Lawson-Hanson Step 1 (LH1): initialize the variables
   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);
+    options = optimset ();
+    ## FIXME: what are the correct defaults?
+    options.MaxIter = 1e5;
+    options.TolX = 1e-8;
   endif
 
   ## Initialize the values
   p = [];
   z = 1:numel (x);
-  ## compute the gradient
+  ## LH2: compute the gradient
   w = c'*(d - c*x);
 
   iter = 0;
-  while (! isempty (z) && any (w(z) > 0) && iter < options.maxiter)
+  ## LH3: test for completion
+  while (! isempty (z) && any (w(z) > 0) && iter < options.MaxIter)
+    ## LH4: find the maximum gradient
     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) = [];
+    ## LH5: move the index from z to p
+    z(z == idx) = [];
+    p(end+1) = idx;
 
     newx = false;
-    while (! newx && iter < options.maxiter)
+    while (! newx && iter < options.MaxIter)
       iter++;
 
+      ## LH6: compute the positive matrix and find the min norm solution
+      ## of the positive problem
       cpos = c;
       cpos(:,z) = 0;
       ## find min norm solution to the positive matrix
@@ -99,20 +108,24 @@ function [x, resnorm, residual, exitflag
       xtmp = (cpos_r\cpos_q')*d;
       xtmp(z) = 0;
       if (all (xtmp >= 0))
-        ## complete the inner loop
+        ## LH7: tmp solution found, iterate
         newx = true;
         x = xtmp;
       else
+        ## LH8, LH9: find the scaling factor and adjust the 
         mask = (xtmp < 0);
-        alpha = min(x(mask)./(x(mask) - xtmp(mask)));
+        alpha = min (x(mask)./(x(mask) - xtmp(mask)));
+        ## LH10: adjust x
         x = x + alpha*(xtmp - x);
+        ## LH11: move from p to z all x == 0;
         idx = find (x == 0);
-        p(ismember (p, x)) = [];
+        p(ismember (p, idx)) = [];
         z = [z idx];
       endif
     endwhile
     w = c'*(d - c*x);
   endwhile
+  ## LH12: complete
 
   ## generate the additional output arguments
   if (nargout > 1)
@@ -122,7 +135,7 @@ function [x, resnorm, residual, exitflag
     residual = d-C*x;
   endif
   exitflag = iter;
-  if (nargout > 3 && iter >= options.maxiter)
+  if (nargout > 3 && iter >= options.MaxIter)
     exitflag = 0;
   endif
   if (nargout > 4)