In my first year at MathWorks (1987!), a professor I know got in touch with me. He was trying to solve a least squares problem with nonnegativity constraints. Having been raised properly, I knew immediately where to get a great algorithm
- Lawson and Hanson, "Solving Least Squares Problems", Prentice-Hall, 1974, Chapter 23, p. 161.
Contents
Fortran Code
The Fortran code is supplied in the book and the professor I spoke with did a fairly faithful translation, despite MATLAB not yet having a goto statement. Problem was, it wasn't working correctly. He asked me to take a look. I remember being excited, because I love nnls as L-H call the algorithm. However, pouring over someone else's translated Fortran loses its pizzazz quite quickly.
Let me give you some statistics on the Fortran code itself. You can find a copy of the code here.
- 326 lines total
- ~150 lines of comments
- leaving ~175 lines of executable code
- main double while-loop, ~125 lines of executable code (lines 83-289)
A Eureka Moment
Procrastination can pay off sometimes. Instead of reading through the M-file code, comparing it to the Fortran line-by-line, I opened up the reference for the algorithm. Here's a glimpse of what I found.
Do you see what I see? I see the problem being talked about in terms of matrices and vectors. Hmmm, maybe I can implement it straight from the algorithm description, bypassing the Fortran implementation. And so I did.
The MATLAB Code
MATLAB shipped with the function nnls, which ultimately was renamed to lsqnonneg and was updated to include some code to standardize it with our other optimization routines. I had the code written and running in less than half a day, including working correctly on the examples sent to me by the professor that were not working in his translated code.
The code statistics look like this for the M-file.
- 206 lines total
- ~80 blank and comment lines
- ~125 lines of executable code, of which ~55 dealing with the infrastructure common to our other optimization routines that don't appear in the Fortran
- ~50 lines for the main iteration loops (nested while-loops)
Lessons Learned
- I relearned from this experience that it is worth going back and thinking about code and the algorithm before diving in. I guess I had to learn this one a few times in the course of my progamming life! While the Fortran and both versions of the M-files were based on the same algorithm, I was able to bypass the code written and write a version using the matrix notation of the reference.
- The MATLAB code is shorter, though I am sure some will say that it is not that much shorter.
- The algorithm, when I wrote it, really sunk in for me, since instead of being mired in loops, I watched variables move in and out of the active set. It was very interesting to see the algorithm work from a less nitty-gritty level, yet still capturing all the details.
Care to share any of your experiences where thinking about an algorithm differently than the way it was posed help you out? Add your thoughts here.
Get
the MATLAB code
Published with MATLAB® 7.3
hi its very knowlegble but i need to know about simulation in the electrical circuit using the matlab
Hi Loren,
Do you have several test cases for the nnls problem that you could email to me? (or do you know where I could get them?)
Thank you,
-Eric
Eric-
I don’t have anything. You can check our documentation and the book I cite above.
–Loren
Thanks Loren, I’ll take a look! (I just got back to investigating this problem)
Loren,
Do you happen to know if there is an alternative version of lsqnonneg which will handle large, sparse matrices?
The current version generates a full matrix at line 159, so is causing me memory issues.
I’ll work on one if there’s not an algorithm floating around… in that event, I’ll send it to you once I finish (IF I finish!).
Kind regards
Tom Clark
Tom-
I’m not aware of one. We’d love to have your contribution on MATLAB Central File Exchange if/when you have one!
–Loren
Just saw Marcelo’s note on the newsgroup:
lsqnonneg was upgraded in version R2008b; it now uses sparse linear algebra if the
matrix is sparse (it uses “\”, backslash).
R2008b Release notes:
http://www.mathworks.com/access/helpdesk/help/techdoc/rn/brqyzsl-1.html
(See “Upgrades to lsqnonneg”.)
-Marcelo
I have a problem about non-neglative least square. I am so exited to see the comand (lsqnonneg) on the website.
Thanks for your contribution sincerely!
Hi Loren,
Thanks for your contribution on the lsqnonneg. It does help me a lot. I am currently looking into some problem invovled with minimizing the standard non-negative least squares and an extra constraint, e.g. the energy constraint to smooth the spikes provided by the nnls and the misfit derived from the modified minimization equation would be a bit larger than the standard nnls. Is there any function built-in matlab could do this?
Thank you,
Jing
Hi Jing,
I am assuming that you don’t have Optimization Toolbox. If you have there are several options to tackle your problem.
In MATLAB I can think of some workarounds that may help you.
(1) If your energy constraint is linear (this is not likely but may be you can linearize it) you can append it to the pair C,d in lsqnonneg(C,d). Note that at the end lsqnonneg will reduce NORM(d-C*X): your linear constraint may be satisfied at this point (no guarantees, though). Another note: if you have a linear inequality, don’t forget to add a slack variable to make your constraint an equality.
(2) If you cannot (don’t want to) linearize your energy constraint, you can use fminsearch by creating your own objective function and adding your energy constraint with a penalty parameter in front, i.e. something like
min NORM(d-C*X) + aPenParameter * g(X)
where g(X) = 0.
I hope this helps.
Best Regards,
Derya
Hi Loren,
The code TSNNLS claims to be substantially faster than the matlab implementation (as of the 2009 writing presumably) with comparable solutions.
http://www.jasoncantarella.com/webpage/index.php?title=Tsnnls
Have you compared the most recent matlab lsqnonneg with the speed and results of this code, or any others, or is the matlab implementation developed without reference to other codes? If these comparisons have been done, do you know where we can find a report about them, and is it possible to run these codes through a mex wrapper from within matlab?
Thanks!
Hi Alex,
lsqnonneg is not benchmarked against other code. I’m unaware of any results comparisons with TSNNLS or other software, although I’m sure they exist somewhere.
You might want to look for comparisons of accuracy as well, if it is a concern. From the short description, it says that it uses a Cholesky factorization, which leads me to believe that the normal equations (including
) are explicitly formed. This can be bad if A is ill-conditioned.
I’m also sure you can make a MEX wrapper to call this code as well.
Regards,
+Steve