Re: Non-negative Least Squares (NNLS) algorithm

*To*: mathgroup at smc.vnet.net*Subject*: [mg43782] Re: Non-negative Least Squares (NNLS) algorithm*From*: poujadej at yahoo.fr (Jean-Claude Poujade)*Date*: Sat, 4 Oct 2003 02:05:00 -0400 (EDT)*References*: <blj6cf$1j6$1@smc.vnet.net>*Sender*: owner-wri-mathgroup at wolfram.com

mdw at ccu1.auckland.ac.nz (Michael Woodhams) wrote in message news:<blj6cf$1j6$1 at 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