Non-negative Least Squares (NNLS) algorithm
- To: mathgroup at smc.vnet.net
- Subject: [mg43770] Non-negative Least Squares (NNLS) algorithm
- From: mdw at ccu1.auckland.ac.nz (Michael Woodhams)
- Date: Fri, 3 Oct 2003 02:28:54 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
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.