{"id":69,"date":"2006-12-12T10:01:16","date_gmt":"2006-12-12T15:01:16","guid":{"rendered":"https:\/\/blogs.mathworks.com\/loren\/?p=69"},"modified":"2017-02-04T14:38:53","modified_gmt":"2017-02-04T19:38:53","slug":"brief-history-of-nonnegative-least-squares-in-matlab","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/loren\/2006\/12\/12\/brief-history-of-nonnegative-least-squares-in-matlab\/","title":{"rendered":"Brief History of Nonnegative Least Squares in MATLAB"},"content":{"rendered":"<div xmlns:mwsh=\"https:\/\/www.mathworks.com\/namespace\/mcode\/v1\/syntaxhighlight.dtd\" class=\"content\">\r\n   <introduction>\r\n      <p>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\r\n         with nonnegativity constraints.  Having been raised properly, I knew immediately where to get a great algorithm\r\n      <\/p>\r\n      <div>\r\n         <ul>\r\n            <li>Lawson and Hanson, \"Solving Least Squares Problems\", Prentice-Hall, 1974, Chapter 23, p. 161.<\/li>\r\n         <\/ul>\r\n      <\/div>\r\n   <\/introduction>\r\n   <h3>Contents<\/h3>\r\n   <div>\r\n      <ul>\r\n         <li><a href=\"#1\">Fortran Code<\/a><\/li>\r\n         <li><a href=\"#3\">A Eureka Moment<\/a><\/li>\r\n         <li><a href=\"#4\">The MATLAB Code<\/a><\/li>\r\n         <li><a href=\"#6\">Lessons Learned<\/a><\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <h3>Fortran Code<a name=\"1\"><\/a><\/h3>\r\n   <p>The Fortran code is supplied in the book and the professor I spoke with did a fairly faithful translation, despite MATLAB\r\n      not yet having a <a href=\"https:\/\/blogs.mathworks.com\/loren\/?p=29\">goto<\/a> statement.  Problem was, it wasn't working correctly.  He asked me to take a look.  I remember being excited, because I love\r\n      <tt>nnls<\/tt> as L-H call the algorithm.  However, pouring over someone else's translated Fortran loses its pizzazz quite quickly.\r\n   <\/p>\r\n   <p>Let me give you some statistics on the Fortran code itself.\r\n   <\/p>\r\n   <div>\r\n      <ul>\r\n         <li>326 lines total<\/li>\r\n         <li>~150 lines of comments<\/li>\r\n         <li>leaving ~175 lines of executable code<\/li>\r\n         <li>main double while-loop, ~125 lines of executable code (lines 83-289)<\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <h3>A Eureka Moment<a name=\"3\"><\/a><\/h3>\r\n   <p>Procrastination can pay off sometimes.  Instead of reading through the M-file code, comparing it to the Fortran line-by-line,\r\n      I opened up the reference for the algorithm.  Here's a glimpse of what I found.\r\n   <\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/69\/nnls.png\"> <\/p>\r\n   <p>Do you see what I see?  I see the problem being talked about in terms of matrices and vectors.  Hmmm, maybe I can implement\r\n      it straight from the algorithm description, bypassing the Fortran implementation.  And so I did.\r\n   <\/p>\r\n   <h3>The MATLAB Code<a name=\"4\"><\/a><\/h3>\r\n   <p>MATLAB shipped with the function <tt>nnls<\/tt>, which ultimately was renamed to <a href=\"https:\/\/www.mathworks.com\/help\/matlab\/ref\/lsqnonneg.html\"><tt>lsqnonneg<\/tt><\/a> and was updated to include some code to standardize it with our other optimization routines.  I had the code written and\r\n      running in less than half a day, including working correctly on the examples sent to me by the professor that were not working\r\n      in his translated code.\r\n   <\/p>\r\n   <p>The code statistics look like this for the M-file.<\/p>\r\n   <div>\r\n      <ul>\r\n         <li>206 lines total<\/li>\r\n         <li>~80 blank and comment lines<\/li>\r\n         <li>~125 lines of executable code, of which ~55 dealing with the infrastructure common to our other optimization routines that\r\n            don't appear in the Fortran\r\n         <\/li>\r\n         <li>~50 lines for the main iteration loops (nested while-loops)<\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <h3>Lessons Learned<a name=\"6\"><\/a><\/h3>\r\n   <div>\r\n      <ul>\r\n         <li>I relearned from this experience that it is worth going back and thinking about code and the algorithm before diving in. \r\n            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\r\n            M-files were based on the same algorithm, I was able to bypass the code written and write a version using the matrix notation\r\n            of the reference.\r\n         <\/li>\r\n         <li>The MATLAB code is shorter, though I am sure some will say that it is not <b>that much<\/b> shorter.\r\n         <\/li>\r\n         <li>The algorithm, when I wrote it, really sunk in for me, since instead of being mired in loops, I watched variables move in\r\n            and out of the active set.  It was very interesting to see the algorithm work from a less nitty-gritty level, yet still capturing\r\n            <b>all<\/b> the details.\r\n         <\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <p>Care to share any of your experiences where thinking about an algorithm differently than the way it was posed help you out?\r\n       Add your thoughts <a href=\"https:\/\/blogs.mathworks.com\/loren\/2006\/12\/12\/brief-history-of-nonnegative-least-squares-in-matlab\/#respond\">here<\/a>.\r\n   <\/p><script language=\"JavaScript\">\r\n<!-- \r\n    function grabCode_cb48cece11b54c01b357893c310f5f5e() {\r\n        \/\/ Remember the title so we can use it in the new page\r\n        title = document.title;\r\n\r\n        \/\/ Break up these strings so that their presence\r\n        \/\/ in the Javascript doesn't mess up the search for\r\n        \/\/ the MATLAB code.\r\n        t1='cb48cece11b54c01b357893c310f5f5e ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' cb48cece11b54c01b357893c310f5f5e';\r\n    \r\n        b=document.getElementsByTagName('body')[0];\r\n        i1=b.innerHTML.indexOf(t1)+t1.length;\r\n        i2=b.innerHTML.indexOf(t2);\r\n \r\n        code_string = b.innerHTML.substring(i1, i2);\r\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\r\n\r\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \r\n        \/\/ in the XML parser.\r\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\r\n        \/\/ doesn't go ahead and substitute the less-than character. \r\n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\r\n\r\n        author = 'Loren Shure';\r\n        copyright = 'Copyright 2006 The MathWorks, Inc.';\r\n\r\n        w = window.open();\r\n        d = w.document;\r\n        d.write('<pre>\\n');\r\n        d.write(code_string);\r\n\r\n        \/\/ Add author and copyright lines at the bottom if specified.\r\n        if ((author.length > 0) || (copyright.length > 0)) {\r\n            d.writeln('');\r\n            d.writeln('%%');\r\n            if (author.length > 0) {\r\n                d.writeln('% _' + author + '_');\r\n            }\r\n            if (copyright.length > 0) {\r\n                d.writeln('% _' + copyright + '_');\r\n            }\r\n        }\r\n\r\n        d.write('<\/pre>\\n');\r\n      \r\n      d.title = title + ' (MATLAB code)';\r\n      d.close();\r\n      }\r\n-->   \r\n      <\/script>\r\n<noscript>\r\n<em>A Javascript-enabled browser is required to use the \"Get the MATLAB code\" link.<\/em>\r\n<\/noscript>\r\n<p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\"><br><a href=\"javascript:grabCode_cb48cece11b54c01b357893c310f5f5e()\"><span style=\"font-size: x-small;        font-style: italic;\">Get \r\n            the MATLAB code<\/span><\/a><br><br>\r\n      Published with MATLAB&reg; 7.3<br><\/p>\r\n<\/div>\r\n<!--\r\ncb48cece11b54c01b357893c310f5f5e ##### SOURCE BEGIN #####\r\n%% Brief History of Nonnegative Least Squares in MATLAB\r\n% In my first year at MathWorks (1987!), a\r\n% professor I know got in touch with me.  He was trying to solve a least\r\n% squares problem with nonnegativity constraints.  Having been raised\r\n% properly, I knew immediately where to get a great algorithm\r\n%\r\n% * Lawson and Hanson, \"Solving Least Squares Problems\", Prentice-Hall,\r\n% 1974, Chapter 23, p. 161.\r\n%% Fortran Code\r\n% The Fortran code is supplied in the book and the professor I spoke with\r\n% did a fairly faithful translation, despite MATLAB not yet having a\r\n% <https:\/\/blogs.mathworks.com\/loren\/?p=29 goto> statement.  Problem was, it\r\n% wasn't working correctly.  He asked me to take a look.  I remember being\r\n% excited, because I love |nnls| as L-H call the algorithm.  However,\r\n% pouring over someone else's translated Fortran loses its pizzazz quite\r\n% quickly.\r\n%%\r\n% Let me give you some statistics on the Fortran code itself.  You can find\r\n% a copy of the code\r\n% <http:\/\/hesperia.gsfc.nasa.gov\/~schmahl\/nnls\/nnls.for here>.\r\n%\r\n% * 326 lines total\r\n% * ~150 lines of comments\r\n% * leaving ~175 lines of executable code\r\n% * main double while-loop, ~125 lines of executable code (lines\r\n% 83-289)\r\n%% A Eureka Moment\r\n% Procrastination can pay off sometimes.  Instead of reading through the\r\n% M-file code, comparing it to the Fortran line-by-line, I opened up the\r\n% reference for the algorithm.  Here's a glimpse of what I found.\r\n%\r\n% <<nnls.png>>\r\n%\r\n% Do you see what I see?  I see the problem being talked about in terms of\r\n% matrices and vectors.  Hmmm, maybe I can implement it straight from the\r\n% algorithm description, bypassing the Fortran implementation.  And so I\r\n% did.  \r\n%% The MATLAB Code\r\n% MATLAB shipped with the function |nnls|, which ultimately was\r\n% renamed to\r\n% <https:\/\/www.mathworks.com\/help\/matlab\/ref\/lsqnonneg.html |lsqnonneg|>\r\n% and was updated to include some code to standardize it with our other\r\n% optimization routines.  I had the code written and running in less than\r\n% half a day, including working correctly on the examples sent to me by the\r\n% professor that were not working in his translated code.\r\n%%\r\n% The code statistics look like this for the M-file.\r\n%\r\n% * 206 lines total\r\n% * ~80 blank and comment lines\r\n% * ~125 lines of executable code, of which ~55 dealing with the\r\n% infrastructure common to our other optimization routines that don't\r\n% appear in the Fortran\r\n% * ~50 lines for the main iteration loops (nested while-loops)\r\n\r\n%% Lessons Learned \r\n% * I relearned from this experience that it is worth going back and thinking\r\n% about code and the algorithm before diving in.  I guess I had to learn\r\n% this one a few times in the course of my progamming life!  While the\r\n% Fortran and both versions of the M-files were based on the same\r\n% algorithm, I was able to bypass the code written and write a version\r\n% using the matrix notation of the reference.\r\n% * The MATLAB code is shorter, though I am sure some will say that it is\r\n% not *that much* shorter.\r\n% * The algorithm, when I wrote it, really sunk in for me, since instead of\r\n% being mired in loops, I watched variables move in and out of the active\r\n% set.  It was very interesting to see the algorithm work from a less\r\n% nitty-gritty level, yet still capturing *all* the details.\r\n%%\r\n% Care to share any of your experiences where thinking about an algorithm\r\n% differently than the way it was posed help you out?  Add your thoughts\r\n% <?p=69#respond here>.\r\n##### SOURCE END ##### cb48cece11b54c01b357893c310f5f5e\r\n-->","protected":false},"excerpt":{"rendered":"<p>\r\n   \r\n      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\r\n         with nonnegativity constraints.  Having been... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/loren\/2006\/12\/12\/brief-history-of-nonnegative-least-squares-in-matlab\/\">read more >><\/a><\/p>","protected":false},"author":39,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[8,12],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/69"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/users\/39"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/comments?post=69"}],"version-history":[{"count":1,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/69\/revisions"}],"predecessor-version":[{"id":2220,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/69\/revisions\/2220"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/media?parent=69"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/categories?post=69"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/tags?post=69"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}