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