Loren on the Art of MATLAB

December 12th, 2006

Brief History of Nonnegative Least Squares in MATLAB

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

By Loren Shure

10:01 UTC | Posted in Readability, Vectorization | Permalink |

You can follow any responses to this entry through the RSS 2.0 feed. You can skip to the end and leave a response. Pinging is currently not allowed.

7 Responses to “Brief History of Nonnegative Least Squares in MATLAB”

  1. nirav pandya replied on March 15th, 2007 at 11:14 UTC :

    hi its very knowlegble but i need to know about simulation in the electrical circuit using the matlab

  2. Eric replied on December 3rd, 2008 at 21:11 UTC :

    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

  3. Loren replied on December 3rd, 2008 at 23:01 UTC :

    Eric-

    I don’t have anything. You can check our documentation and the book I cite above.

    –Loren

  4. Eric replied on April 23rd, 2009 at 18:18 UTC :

    Thanks Loren, I’ll take a look! (I just got back to investigating this problem)

  5. Tom Clark replied on June 26th, 2009 at 13:42 UTC :

    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

  6. Loren replied on June 26th, 2009 at 13:56 UTC :

    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

  7. Loren replied on June 26th, 2009 at 19:33 UTC :

    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

Leave a Reply

Wrap code fragments inside <pre> tags, like this:

<pre class="code">
a = magic(3);
sum(a)
</pre>

If you have a "<" character in your code, either follow it with a space or replace it with "&lt;" (including the semicolon).

About

Loren Shure works on design of the MATLAB language at The MathWorks. She writes here about once a week on MATLAB programming and related topics.

Subscribe


Recent Comments

  • Jun: I totally can not believe it, Loren. You are really helpful. Thank you so much, MATLAB master!
  • Loren: Wow folks- Always lots of interest when there’s a quickie to try out! I will only make 2 general...
  • Loren: Jun- ismember is your friend here: >> [aa,ind] = ismember(Array2,Arra y1) aa = 1 1 1 1 1 1 1 ind = 1 2 1 4 4 3...
  • Dan: I like the first way better than the second way. Combining the arrays into one and running any is nice, although...
  • James Myatt: How about I = (a == 0 | b == 0); a(I) = []; b(I) = [];
  • Tunc: Hello Loren, love your blog because of such inspiring and challenging comments to such ’small’...
  • Pekka Kumpulainen: Here is my tradeoff. I usually want to keep the original variables as they are most probably...
  • Iain: Followup: Of course, to allow NaNs (counting them as non-zero): mask = (a~=0) & (b~=0); The mask says “a...
  • Matt Fig: I would usually go with something like this: y = a&b; x = a(y); y = b(y); But I was surprised to find...
  • kk: c=all([a;b]) a(c) a(b)

Blog Admin

These postings are the author's and don't necessarily represent the opinions of The MathWorks.