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.