Non-negative Least Squares (NNLS) algorithm

A little while ago, I looked at this group and the web for an
implementation of the NNLS algorithm, but all I found was a post in
this group from 1997 from somebody also looking for an implementation,
so I had to do it myself.

The problem: Given a matrix A (usually more rows than columns) and
vector f, find vector x such that:
1) All elements of x are non-negative
2) subject to contraint (1), the Euclidian norm (2-norm) of vector
(A.x-f) is minimized.

(Note that the unconstrained problem - find x to minimize (A.x-f) - is
a simple application of QR decomposition.)

The NNLS algorithm is published in chapter 23 of Lawson and Hanson,
"Solving Least Squares Problems" (Prentice-Hall, 1974, republished
SIAM, 1995)

Some preliminary comments on the code:
1) It hasn't been thoroughly tested. 
2) I only have a few months experience programming Mathematica (but
lots of programming experience in general) so there are probably style
issues. I am happy to receive constructive criticism.
3) There is some ugliness where a chunk of code is repeated twice to
be able to have a simple 'while' loop, rather than one that breaks out
in the middle. Arguably breaking the middle of the loop would be
better than repeating the code.
4) I've left in a bunch of debugging statements. 
5) I've cut-and-pasted this from my notebook, which is a bit ugly with
the "\[LeftDoubleBracket]" instead of "[[" etc.
6) I haven't paid much attention to efficiency - e.g. I multiply by
"Inverse[R]"  rather than trying to do a backsubstitution (R is a
triangular matrix.)

====================================================================

(* Coded by Michael Woodhams, from algorithm by Lawson and Hanson, *)
(* "Solving Least Squares Problems", 1974 and 1995. *)
bitsToIndices[v_] := Select[Table[i, {i, Length[v]}], v[[#]] == 1 &];
NNLS[A_, f_] := Module[
      {x, zeroed, w, t, Ap, z, q, \[Alpha], i, zeroedSet, positiveSet,
        toBeZeroed, compressedZ, Q, R},
      (* Use delayed evaluation so that these are recalculated on the
fly as \
needed : *)
      zeroedSet := bitsToIndices[zeroed];
      positiveSet := bitsToIndices[1 - zeroed];
      (* Init x to vector of zeros, same length as a row of A *)
      debug["A=", MatrixForm[A]];
      x = 0 A\[LeftDoubleBracket]1\[RightDoubleBracket];
      debug["x=", x];
      (* Init zeroed to vector of ones,
        same length as x *)
      zeroed = 1 - x;
      debug["zeroed=", zeroed];
      w = Transpose[A].(f - A.x);
      debug["w=", w];
      While[zeroedSet != {} 
          && Max[w\[LeftDoubleBracket]zeroedSet\[RightDoubleBracket]]
> 0,
        debug["Outer loop starts."];
        (* The index t of the largest element of w, *)
        (* subject to the constraint t is zeroed *)
        t = 
          Position[w zeroed, Max[w zeroed], 1, 
                1]\[LeftDoubleBracket]1\[RightDoubleBracket]\
\[LeftDoubleBracket]1\[RightDoubleBracket];
        debug["t=", t];
        zeroed\[LeftDoubleBracket]t\[RightDoubleBracket] = 0;
        debug["zeroed=", zeroed];
        (* Ap = the columns of A indexed by positiveSet *)
        Ap = 
          Transpose[
            Transpose[A]\[LeftDoubleBracket]
              positiveSet\[RightDoubleBracket]];
        debug["Ap=", MatrixForm[Ap]];
        (* Minimize (Ap . compressedZ - f) by QR decomp *)
        {Q, R} = QRDecomposition[Ap];
        compressedZ = Inverse[R].Q.f;
        (* 
          Create vector z with 0 in zeroed indices and compressedZ
entries \
elsewhere *)
        z = 0 x;
        z\[LeftDoubleBracket]positiveSet\[RightDoubleBracket] =
compressedZ;
        debug["z=", z];
        While[Min[z] < 0,
          (* There is a wart here : x can have zeros, 
            giving infinities or indeterminates. They don't matter, 
            
            as we ignore those elements (not in postitiveSet) but it
will \
produce warnings. *)
          debug["Inner loop start"];
          (* 
            find smallest x\[LeftDoubleBracket]
                  q\[RightDoubleBracket]/(x\[LeftDoubleBracket]q\
\[RightDoubleBracket] - z\[LeftDoubleBracket]q\[RightDoubleBracket])
*)
          (* such that : q is not zeroed, 
               z\[LeftDoubleBracket]q\[RightDoubleBracket] < 0 *)
          \[Alpha] = Infinity;
          For[q = 1, q <= Length[x], q++,
            
            If[zeroed\[LeftDoubleBracket]q\[RightDoubleBracket] == 0
&&
                z\[LeftDoubleBracket]q\[RightDoubleBracket] < 0,
              \[Alpha] = 
                Min[\[Alpha], 
                  x\[LeftDoubleBracket]q\[RightDoubleBracket]/(x\
\[LeftDoubleBracket]q\[RightDoubleBracket] - 
                        z\[LeftDoubleBracket]q\[RightDoubleBracket])];
              debug["After trying index q=", q, " \[Alpha]=",
\[Alpha]];
              ]; (* if *)
            ]; (* for *)
          debug["\[Alpha]=", \[Alpha]];
          x = x + \[Alpha](z - x);
          debug["x=", x];
          
          toBeZeroed = 
            Select[positiveSet, 
              Abs[x\[LeftDoubleBracket]#\[RightDoubleBracket]] <
10^-13 &];
          debug["toBeZeroed=", toBeZeroed];
          zeroed\[LeftDoubleBracket]toBeZeroed\[RightDoubleBracket] =
1;
          x\[LeftDoubleBracket]toBeZeroed\[RightDoubleBracket] = 0;
          
          (* Duplicated from above *)
          (* Ap = the columns of A indexed by positiveSet *)
          
          Ap = Transpose[
              Transpose[
                  A]\[LeftDoubleBracket]positiveSet\[RightDoubleBracket]];
          debug["Ap=", MatrixForm[Ap]];
          (* Minimize (Ap . compressedZ - f) by QR decomp *)
          {Q, R} = QRDecomposition[Ap];
          compressedZ = Inverse[R].Q.f;
          (* 
            Create vector z with 0 in zeroed indices and compressedZ
entries \
elsewhere *)
          z = 0 x;
          
          z\[LeftDoubleBracket]positiveSet\[RightDoubleBracket] = 
            compressedZ;
          debug["z=", z];
          ]; (* end inner while loop *)
        x = z;
        debug["x=", x];
        w = Transpose[A].(f - A.x);
        debug["w=", w];
        ]; (* end outer while loop *)
      Return[x];
      ]; (* end module *)

====================================================================

And some final comments:

Don't 'reply' to the e-mail address in the header of this post - it is
old and defunct, and not updated because of spammers. I am at
massey.ac.nz, and my account name is m.d.woodhams. (3 year post-doc
appointment, so by the time you read this, that address might also be
out of date.)

I put this code in the public domain - but I'd appreciate it if you
acknowledge my authorship if you use it. I will be writing a
bioinformatics paper about "modified closest tree" algorithm that uses
this, and I might (not decided about this) give a link to Mathematica
code, which will include the above. If this happens, you could cite
that paper.

This software comes with no warrantee etc etc. Your warrantee is that
you have every opportunity to test the code yourself before using it.

If you are reading this long after it was posted, I may have an
improved version by then. Feel free to enquire by e-mail.

0
mdw1 (2)
10/3/2003 6:53:03 AM
comp.soft-sys.math.mathematica 28833 articles. 0 followers. Follow

1 Replies
371 Views

Similar Articles

[PageSpeed] 6
mdw@ccu1.auckland.ac.nz (Michael Woodhams) wrote in message news:<blj6cf$1j6$1@smc.vnet.net>...
> A little while ago, I looked at this group and the web for an
> implementation of the NNLS algorithm, but all I found was a post in
> this group from 1997 from somebody also looking for an implementation,
> so I had to do it myself.
> 
> The problem: Given a matrix A (usually more rows than columns) and
> vector f, find vector x such that:
> 1) All elements of x are non-negative
> 2) subject to contraint (1), the Euclidian norm (2-norm) of vector
> (A.x-f) is minimized.
[...]

Michael,

Congratulations for your NNLS code!
I had a 10-liner naive solution (see below) using FindMinimum. 
It seems to work most of the time, but your solution is more accurate 
and... at least sixty times faster!

NNLSFindMinimum[A_, f_] := 
Module[{nbx = Length[First[A]], xi, x, axf, xinit},
      xi = Array[x,nbx];
      axf = A.xi^2 - f;
      xinit = PseudoInverse[A].f;
      If[And@@( #>=0& /@ xinit),
        	xinit,
        	fm = FindMinimum[Evaluate[axf.axf],
                Evaluate[Sequence@@Transpose[{xi, xinit}]]];
        	xi^2 /. fm[[2]]
      ]
];

(* a small test : *)
a = Table[Random[],{i,10},{j,20}];
f = Table[Random[],{i,10}];
x = NNLSFindMinimum[a,f]//Timing
% // Chop[#, 10.^-5] &
{1.001 Second,{0.0580067,0.0246345,0,0.0227949,0,0,0,0.386312,0,0,0,0,0,0,0,0,
0,0,0.245039,0.215131}}

x = NNLS[a, f] // Timing
{0.016 Second,{0.058114,0.02338271,0,0.0227961,0,0,0,0.386253,0,0,0,0,0,0,0,0,
0,0,0.246019,0.215755}}
---
jcp

0
poujadej (31)
10/4/2003 6:15:38 AM
Reply:
Similar Artilces:

Greedy and non greedy quantifiers
Im just after a bit of clarification with quantifiers in regular expressions. I just want to be sure of the differences between the quantifiers. so for these various regular expressions, [a-z]* - this will match any amount lower case letters [a-z]+ - this will match any amount lower case letters (whats difference between + and * in this case?) [a-z]+? - or \d* - This will match any amount of digits \d*? - This will only match none or one number Please can someone offer some clarification Im still unsure of myself with these expresions as I'm new to ruby, Thanks, Dan ...

resampling / interpolation of non-uniform spaced, non-monotonic data #2
hi there, I've tried to find a solution in the helpfiles on the internet, but i really can't find it. i want to resample or interpolate non-uniform spaced, non-monotonic data to uniform spaced data. (velocity profile) The new datapoints have to be at the old data, no smoothing or so. The problem seems to be that my data is not monotonic as interp and interp1 doesn't do the trick. I've tried some other functions also, but they do't do what i want. some code.. H = [ 5; 8.3; 18; 21; 33; 54.4; 101.3]; S = [ 3; 9; 15; 16; 28; 10; 30;]; x = 0:1:120; y = i...

Using non-dict namespaces in functions
Inspired by the new collections.ChainMap in Python 3.3 http://docs.python.org/dev/library/collections.html#collections.ChainMap I would like to experiment with similar non-dictionary namespaces in Python 3.2. My starting point is these two recipes, adapted for Python 3.2: http://code.activestate.com/recipes/305268/ http://code.activestate.com/recipes/577434/ Or for simplicity, here's a mock version: from collections import Mapping class MockChainMap(Mapping): def __getitem__(self, key): if key == 'a': return 1 elif key == 'b': return 2 ra...

Using ASDF for building non-CL based systems?
Does anyone here have experience in using ASDF to build a non-Common Lisp based system? I have a medium-sized personal project currently using MS Visual Studio that I'd like to port to a more cross platform and flexible build system. I'd prefer to use a build system that is highly flexible and doesn't require me to learn Yet Another Syntax -- in other words I want a Lisp based system. ASDF seems like it would do exactly what I want, if I just knew how to customize it. However, I can't find any examples of using ASDF for non-Lisp systems, and its documentation is extr...

Non linear IVP #2
hi I would like to know how to solve this equation using MATLAB and Newton Raphson method. dh/dt + {3S(h''')+3A/h^4*(h') }^1/n * (h^((n+1)/n) * h') = {-1/2n+1 * (3S(h''')+ (3A/h^4)(h') }^(1-n/n) * h^(2n+1/n) * {3s(h'''')+ 3A/h^4 (h'') } IC: h(x,0)=1+B*sin(qx) Spatial periodic boundary conditions: 0<x<(2*pi/q) n--> variable values (0.5, 0.75, 0.9, 1, 2) B--> varies from (0.1 - 0.8) q--> varies from (0.35 - 1.06) A=A'/(6*pi*rho*h0^3*D) D=rho*(rho*ho^n/m)^(2/n-2) s=(sigma/3ho)*D h=1 ho=5.9*10^-6 m=0.0...

Can we 'twist' a digital picture to be square with the camera?
What 'trick' am I missing to take a photo of my watercolor paintings? Can we 'twist' the digital photo to be square with the camera? I must be doing something wrong. When I edit my digital photographs of my watercolor paintings, I need to crop the photo but I can NEVER get the crop edges to be exactly square with the digital camera photograph of the painting. Therefore I always find myself losing a few millimeters of painting so that the photograph is square. What am I doing wrong? Here is how I take the picture. a) I put my unmatted 30 inch by 22.5 inch watercolor on a whit...

why is operator+=() allowed to be non-member function also ?
I am asking this question for learning purpose only. I will not use operator+=() as non-member. Consider the following program: #include <iostream> #include <cstdlib> using namespace std; class Test { public: Test(int arg = 10); int val; }; Test::Test(int arg) : val(arg) { cout << "one arg ctor called" << endl; } Test operator+=(Test lhs, Test rhs) { cout << "from operator+=" << endl; return Test(lhs.val + rhs.val); } int main() { Test obj; 5 += obj; ...

non-member extractors
#include <fstream> #pragma hdrstop #include <condefs.h> #pragma argsused int main(int argc, char* argv[]) { std::fstream I; I.open("C:\\File.ext"); if (I.is_open()) { unsigned short a; I >> a; return 0; } else return 1; } Is there anything wrong with this? Why does the reading function fail? Fraser. Fraser Ross wrote: > #include <fstream> > #pragma hdrstop > #include <condefs.h> > #pragma argsused > int main(int argc, char* argv[]) > { > std::fstream I; > I.open("C:\\File.ext"); > if (I.is_open...

Automatic gain control / exposure algorithm based on image histogram?
I'm not sure whether this is a good place to post this question... I did some research and it seems that regular digital cameras use hardware circuits and light reflectance meters for this. But I have a need to implement a decent algorithm that would 1. take a picture from a camera, 2. calculate the histogram, 3. based on the histogram adjust the gain / exposure settings on the camera, go to step 1 again and loop until certain criteria are matched. Right now, I have a primitive version in which the exposure is constant and the gain is the only thing that's changed. The algorithm calc...

PS CS3 puzzler
Can anyone shed light on this? 1) In PS, follow the instructions at: http://www.luminous-landscape.com/essays/test-charts.shtml and generate the Granger chart as shown. Oddly, it has a strange distribution given that the gradient layer is linear top-to-bottom. What is creating those diagonals in the chart? 2) For comparison, see; http://farm3.static.flickr.com/2790/4337946696_8ef5e104ff_o.jpg -- gmail originated posts are filtered due to spam. Alan Browne wrote: > Can anyone shed light on this? > > 1) In PS, follow the instructions at: > > htt...

Re: wildcard as empty non terminal #2
> ... your system should see that when it reaches node *.b that there > is a child node,or at least flag that this name has a child. > > Thus the wildcard algorithm should be disabled for that *.b. i agree. -- to unsubscribe send a message to namedroppers-request@ops.ietf.org with the word 'unsubscribe' in a single line as the message text body. archive: <http://ops.ietf.org/lists/namedroppers/> ...

Pointer to non-static member functions
Hi, I'm experimenting with function pointers and found two questions. Let's assume this code: 1 #include <iostream> 2 class A; 3 4 //////////////////////////////////////////// 5 class B 6 { 7 public: 8 void printB( void (A::*)(void) ); 9 }; 10 11 void B::printB( void (A::*func)(void) ) 12 { 13 *func(); // not working 14 } 15 16 //////////////////////////////////////////// 17 class A 18 { 19 public: 20 void printA(void); 21 void invokeB(void); 22 23 private: 24 B myB; 25 }; 26 27 void A::printA(void) 28 { 29 std::c...

non-blocking TCPServer ?
hi, maybe i'm understanding how this is meant to work incorrectly, but i would have assumed TCPServer#accept would return immediately if i had set it to nonblocking. the accept syscall on Linux would not have blocked, according to the manpage.. require 'socket' require 'io/nonblock' s = TCPServer.new('localhost', 9000) s.nonblock = true s.accept $stderr.puts "returned" when i run this in a test script, accept blocks. $ ruby -v ruby 1.8.2 (2004-08-24) [i386-linux] leon Hi, In message "Re: non-blocking TCPServer ?" on Mon, 27 Sep 200...

Highly Non-Stationary Signals for Adaptive Filters
Hello, I encountered a problem with Noise Cancellation. The interfering signal that I am trying to cancel goes on and off really fast (noise of a switch). I have never quantified how long it turns on and off but it it is less than 1 ms or so. I am attempting to cancel its noise by applying an adaptive filter and my application works for stationary signals such as band-limited white Gaussian noise, but it doesn't work for this problem. Assuming that the noise of the interfering signal is the same when it turns on and off, can I turn off the adaptation when the switch is...

Re: [LogoForum] Drawing a square with different styles #2 1987667
> G'day G'day Folks, > > I am plodding along trying to make sense of the structure of Logo. > > Most introductions give an interative approach to forming a square. It > seems like the most reasonable approach since squares have a known > number of sides ie 4. > > to square :side > repeat 4 [forward :side rt 90] > end > > However I have come across the hint that iterations are not the > preferred method in Logo. The honours go to tail end recursion. Put > simply I'm unfamiliar with this technique so I had a look at producing > a squar...

how to reboot a non respond machine
Hi I move mouse and click some places on screen. The screen is then blank and have no respond. I try Ctrl-Alt-Del to reboot, but it did not help. Esc did not help too. I do not want to turn off the machine (unnormally) because it can lead to trouble later.How can I solve this problem? Ngoc Thanks for help "ngoc" <linh@chello.no> wrote in message news:T_WDb.6771$Y06.106384@news4.e.nsc.no... > Hi > I move mouse and click some places on screen. The screen is then blank > and have no respond. I try Ctrl-Alt-Del to reboot, but it did not help. > Esc did not help to...

Getting state of a non-exclusive button (checkbox) widget
I'm building a UI for a program, part of which requires that if a non- exclusive button (checkbox) is unset - this is the default state - then a number of text widgets and their associated labels are hidden or inactive (I've chosen to make them inactive for consistency's sake) - conversely, if the checkbox is set then the widgets should be active and, in the case of the text widgets editable. In pseudocode, the event handler would be something like this... ; Get state information. ; Identify the widget that caused the event. if event.cause == checkbox then { if ch...

Scale factors for Baum Welch Algorithm
Suppose that I have an HMM model of 3 states and 8 symbols. Suppose I have an observation sequence of 10 tokens, corresponding to the following sequence of indices: 2 3 5 0 8 1 2 4 7 6 Now, here is my forward probabilities matrix alpha_ij: alpha_ij: | 0 1 2 3 4 5 6 7 8 9 <--- The matrix index | 2 3 5 0 8 1 2 4 7 6 <---- The observation sequence --------------------------------------------- 0 | 1 | 2 | Of course, each and every slot of the matrix is filled with some double number representing the probabilities. ...

Non-reflective material for IR
Which materials have the non-reflective property to infra-red (IR) light ? Hi I would say holes (as in wheels) and the colour black, as it gets hot in strong sunlight, showing that it absorbs IR (heat). Cheers |-] Dale <Ayem Cipteke> wrote in message news:439a96b3$0$9286$afc38c87@news.optusnet.com.au... > Which materials have the non-reflective property to infra-red (IR) light ? > Ayem Cipteke wrote: > > Which materials have the non-reflective property to infra-red (IR) light ? If you're talking about IR light in the 780-1000 nan...

Algorithm help, any suggestions???
I'm developing a program which contains some fairly tricky problems that I need help solving. I need someone who is very good with algorithms to help figure out some of the problems I'm having. I wouldn't really need any code written, just someone to provide me with the concepts I need to write it myself. I'd be willing to pay or if anyone wanted to do it for fun that would be nice too. If anyone is interested or knows where I should look for this kind of help (either for money or free if possible), please let me know. The problems are too complex that I don't think I...

non
Hi, I have this kind of a problem. I don't know how to implement soft outpu viterbi algorithm for non-binary (over ring) codes. What should the th extrinsic information contain and how to compute and process it??. I went through many papers relating sova and non binary case, but stil I�m stuck. Thank You Filip ...

Algorithms
A while ago, someone in either this ng or a.l.j.p posted a question about converting roman numbers to Arabic numbers. I thought about it for a while and came up with an idea of how it could be done, but I found that my solution was more human oriented than machine oriented. I hope you know what I mean, but if you don't I'll try to explain: If I'm sitting down looking at a problem, I tend to think it through in terms of how I, a human, would do it, like playing chess, or converting XXXVIII to Arabic. Of course, for the most part, many of these problems are just things tha...

Non-stationary, non-linear data... how to analyze? what method?
I am trying to transform some non-stationary data into the frequency domain to see if there are any under lying frequencies in the data. A sample of the data can be seen here: http://www.geocities.com/uebermenchens/00805HR_TN_1.PNG http://www.geocities.com/uebermenchens/0148H10_TN_3.PNG http://www.geocities.com/uebermenchens/1016727-AL_3.PNG I was just wondering what the best way to do this is?? maybe a certian transform? If I am pointed in the right direction i will do further investigation. The data is non-stationary and non-linear...CORRECT?? is there a way to analyze this i...

a-squared HiJackFree 3.0.0.387
a-squared HiJackFree is a detailed system analysis tool which helps advanced users to detect and remove all types of HiJackers, Spyware, Adware, Trojans and Worms. Manage all types of Autoruns on your system, Explorer and Browser plugins (BHOs, Toolbars, etc.). running Processes and their associated modules, control all Services, even those Windows doesn't display, view open ports and the associated listening processes, edit the hosts file, manage installed Layered Service Providers (LSPs) and analyze the system configuration with using our live online analysis. More Info: http:...

Non-blocking socket
Hi, well, working on sockets is not my thing, but I can't get around the fact t= hat IDL "waits" for an answer and this is not the case elsewhere. Let me explain... I open a socket this way: SOCKET, Unit, IP, Port, $ CONNECT_TIMEOUT=3DConnect_Timeout, $ READ_TIMEOUT=3DRead_Timeout, $ WRITE_TIMEOUT=3Dself.Write_Timeout, $ Error=3DRetError, /GET_LUN The IP is a local one and points to another software on the pc (Win 7-64 bi= ts) that acts as socket server. The timeouts are set to 5, 1 and 1 seconds = for connection, read and write, re...