home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / MATFUN.DI$ / NNLS.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  2.0 KB  |  74 lines
  1. function [x,w] = nnls(E,f,tol)
  2. %NNLS    Non-negative least-squares.
  3. %    X = NNLS(A,b) returns the vector X that solves A*x = b
  4. %    in a least squares sense, subject to x >= 0. 
  5. %
  6. %    A default tolerance of TOL = MAX(SIZE(A)) * NORM(A,1) * EPS
  7. %    is used for deciding when elements of X are less than zero.
  8. %    This can be overridden with X = NNLS(A,b,TOL).
  9. %
  10. %    [X,W] = NNLS(A,b) also returns dual vector W where
  11. %    w(i) < 0 where x(i) = 0 and w(i) = 0 where x(i) > 0.
  12. %
  13. %    See also LSCOV, SLASH.
  14.  
  15. %    L. Shure 5-8-87
  16. %    Revised, 12-15-88,8-31-89 LS.
  17. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  18.  
  19. % Reference:
  20. %  Lawson and Hanson, "Solving Least Squares Problems", Prentice-Hall, 1974.
  21.  
  22. % initialize variables
  23. if nargin < 3
  24.     tol = 10*eps*norm(E,1)*max(size(E));
  25. end
  26. [m,n] = size(E);
  27. P = zeros(1,n);
  28. Z = 1:n;
  29. x = P';
  30. ZZ=Z;
  31. w = E'*(f-E*x);
  32.  
  33. % set up iteration criterion
  34. iter = 0;
  35. itmax = 3*n;
  36.  
  37. % outer loop to put variables into set to hold positive coefficients
  38. while any(Z) & any(w(ZZ) > tol)
  39.     [wt,t] = max(w(ZZ));
  40.     t = ZZ(t);
  41.     P(1,t) = t;
  42.     Z(t) = 0;
  43.     PP = find(P);
  44.     ZZ = find(Z);
  45.     nzz = size(ZZ);
  46.     EP(1:m,PP) = E(:,PP);
  47.     EP(:,ZZ) = zeros(m,nzz(2));
  48.     z = pinv(EP)*f;
  49.     z(ZZ) = zeros(nzz(2),nzz(1));
  50. % inner loop to remove elements from the positive set which no longer belong
  51.     while any((z(PP) <= tol))
  52.         iter = iter + 1;
  53.         if iter > itmax
  54.             error(['Iteration count is exceeded.', ...
  55.                    '  Try raising the tolerance.'])
  56.         end
  57.         QQ = find((z <= tol) & P');
  58.         alpha = min(x(QQ)./(x(QQ) - z(QQ)));
  59.         x = x + alpha*(z - x);
  60.         ij = find(abs(x) < tol & P' ~= 0);
  61.         Z(ij)=ij';
  62.         P(ij)=zeros(1,max(size(ij)));
  63.         PP = find(P);
  64.         ZZ = find(Z);
  65.         nzz = size(ZZ);
  66.         EP(1:m,PP) = E(:,PP);
  67.         EP(:,ZZ) = zeros(m,nzz(2));
  68.         z = pinv(EP)*f;
  69.         z(ZZ) = zeros(nzz(2),nzz(1));
  70.     end
  71.     x = z;
  72.     w = E'*(f-E*x);
  73. end
  74.