{"id":191,"date":"2009-07-08T18:19:26","date_gmt":"2009-07-08T18:19:26","guid":{"rendered":"https:\/\/blogs.mathworks.com\/loren\/2009\/07\/08\/a-brief-history-of-polyval\/"},"modified":"2018-01-08T15:22:32","modified_gmt":"2018-01-08T20:22:32","slug":"a-brief-history-of-polyval","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/loren\/2009\/07\/08\/a-brief-history-of-polyval\/","title":{"rendered":"A Brief History of polyval"},"content":{"rendered":"<div xmlns:mwsh=\"https:\/\/www.mathworks.com\/namespace\/mcode\/v1\/syntaxhighlight.dtd\" class=\"content\">\r\n   <introduction>\r\n      <p>When I first started working at MathWorks, there were already a few functions in MATLAB for working with polynomials.  One\r\n         of these, <a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/ref\/polyval.html\"><tt>polyval<\/tt><\/a>, for polynomial evaluation, used one algorithm, based on <a href=\"http:\/\/en.wikipedia.org\/wiki\/Horner_scheme\">Horner's rule<\/a>, known for its efficiency, especially suitable for embedded controller environments. After a customer wondered why the function\r\n         took so long to evaluate a given polynomial for a scalar input value <tt>x<\/tt>, we updated <tt>polyval<\/tt> to use a faster algorithm for the scalar case, and that algorithm uses the function <a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/ref\/filter.html\"><tt>filter<\/tt><\/a>. I thought I'd discuss a few more details here.\r\n      <\/p>\r\n      <p>Note: since I wrote this post (but before I published it) there was a <a>discussion<\/a> on the newsgroup about this same topic.  Be sure to read the reply from <a href=\"https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/authors\/679\"><tt>John D'Errico<\/tt><\/a>.\r\n      <\/p>\r\n   <\/introduction>\r\n   <h3>Contents<\/h3>\r\n   <div>\r\n      <ul>\r\n         <li><a href=\"#1\">Horner's Rule<\/a><\/li>\r\n         <li><a href=\"#4\">The Scalar Case, Using filter<\/a><\/li>\r\n         <li><a href=\"#11\">How Close are the Algorithms?<\/a><\/li>\r\n         <li><a href=\"#13\">Timing<\/a><\/li>\r\n         <li><a href=\"#14\">Multiple Algorithms<\/a><\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <h3>Horner's Rule<a name=\"1\"><\/a><\/h3>\r\n   <p>Let's first look at the code for evaluating a polynomial for multiple values of input <tt>x<\/tt>.  Here's the relevant portion of <tt>polyval<\/tt>, <tt>nc<\/tt> is the length of the vector representing the polynomial (one more than the polynomial degree).\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">dbtype <span style=\"color: #A020F0\">polyval<\/span> <span style=\"color: #A020F0\">62:67<\/span><\/pre><pre style=\"font-style:oblique\">\r\n62    % Use Horner's method for general case where X is an array.\r\n63    y = zeros(siz_x, superiorfloat(x,p));\r\n64    if nc&gt;0, y(:) = p(1); end\r\n65    for i=2:nc\r\n66        y = x .* y + p(i);\r\n67    end\r\n\r\n<\/pre><p>Here's what a general polynomial looks like:<\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/191\/polynomEvaluation_eq28167.png\"> <\/p>\r\n   <p>This equation can be rewritten like this,<\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/191\/polynomEvaluation_eq22573.png\"> <\/p>\r\n   <p>and is the basis for evaluating the polynomial for an array <tt>x<\/tt>.  Notice that the formula basically starts with the coefficient of the highest degree (named <img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/191\/polynomEvaluation_eq53874.png\"> ) multiplied by <tt>x<\/tt>, then adds the next coefficient, multiplies that result by <tt>x<\/tt> again, adds the next coefficient, and so on, until reaching the final coefficient (the constant term, for <img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/191\/polynomEvaluation_eq35173.png\"> ). In this formulation, we never need to explicitly create all the powers of <tt>x<\/tt> to evaluate the polynomial.\r\n   <\/p>\r\n   <h3>The Scalar Case, Using filter<a name=\"4\"><\/a><\/h3>\r\n   <p>Here's the code for the scalar evalation case using the function <a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/ref\/filter.html\"><tt>filter<\/tt><\/a>.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">dbtype <span style=\"color: #A020F0\">polyval<\/span> <span style=\"color: #A020F0\">52:53<\/span><\/pre><pre style=\"font-style:oblique\">\r\n52        y = filter(1,[1 -x],p);\r\n53        y = y(nc);\r\n\r\n<\/pre><p>It may not, at first glance, appear relevant to polynomial evaluation. However, taking a close look at the mathematics reveals\r\n      that it indeed can be useful.  First, let's examine the <tt>filter<\/tt> formula.  In this version, the coefficients for <tt>a<\/tt> have been normalized so <tt>a(1)<\/tt> = 1 (which is true for the inputs to <tt>filter<\/tt> you've just seen).\r\n   <\/p>\r\n   <p>Now, let's try a small example. Suppose we want to evaluate the same third degree polynomial I used in another <a href=\"https:\/\/blogs.mathworks.com\/loren\/category\/solving-equations\/\">set of posts<\/a> recently.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">p = [1 0 1 -1];<\/pre><p>which represents<\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/191\/polynomEvaluation_eq19758.png\"> <\/p>\r\n   <p>When you call <tt>filter<\/tt> like this <tt>y = filter(b,a,x)<\/tt>, the results are\r\n   <\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/191\/polynomEvaluation_eq44672.png\"> <\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/191\/polynomEvaluation_eq56974.png\"> <\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/191\/polynomEvaluation_eq20393.png\"> <\/p>\r\n   <p>Looking again at the <tt>polyval<\/tt> code, use these substitutions for the inputs to <tt>filter<\/tt>: <tt>b = 1<\/tt>, <tt>a = [1 -x]<\/tt>, and <tt>x = p<\/tt>. We evaluate the filter for the values <tt>p<\/tt>, the polynomial coefficients! Let's do a little math to see why <tt>filter<\/tt> works.  Following the formula for the output of the function <tt>filter<\/tt>, and assuming <tt>p(0) = 0<\/tt>, we can find the first output.\r\n   <\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/191\/polynomEvaluation_eq96341.png\"> <\/p>\r\n   <p>And here are the rest, going up to the final coefficient in our example.<\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/191\/polynomEvaluation_eq32626.png\"> <\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/191\/polynomEvaluation_eq54475.png\"> <\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/191\/polynomEvaluation_eq86091.png\"> <\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/191\/polynomEvaluation_eq70838.png\"> <\/p>\r\n   <p>It sure looks like the final answer here is the polynomial <tt>p<\/tt> evaluated at <tt>x<\/tt>.\r\n   <\/p>\r\n   <h3>How Close are the Algorithms?<a name=\"11\"><\/a><\/h3>\r\n   <p>Numerically, the two algorithms end up being essentially the same.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">x = 0:0.001:1;\r\nresultArray = polyval(p,x);\r\nresultScalar = zeros(size(x));\r\n<span style=\"color: #0000FF\">for<\/span> index = 1:numel(x)\r\n    resultScalar(index) = polyval(p,x(index));\r\n<span style=\"color: #0000FF\">end<\/span>\r\nallequal = isequal(resultArray,resultScalar)<\/pre><pre style=\"font-style:oblique\">allequal =\r\n     1\r\n<\/pre><p>Identical!<\/p>\r\n   <h3>Timing<a name=\"13\"><\/a><\/h3>\r\n   <p>I mentioned we made this change for speed.  This particular code is tricky to check timing since, especially if the polynomial\r\n      degree is small, the number of mathematic operations can be swamped by the timing mechanism.  When computers were slower,\r\n      and before we increased performance of for-loops in MATLAB, measuring the performance differences between these code alternatives\r\n      was much easier than it is today.\r\n   <\/p>\r\n   <h3>Multiple Algorithms<a name=\"14\"><\/a><\/h3>\r\n   <p>Do you try to write your code as one size fits all, or do you create code in which you switch algorithms or implementations,\r\n      such as we've seen in <tt>polyval<\/tt>, to account for different conditions such as array size?  I'd love to hear more about your experiences with these ideas <a href=\"https:\/\/blogs.mathworks.com\/loren\/?p=191#respond\">here<\/a>.\r\n   <\/p><script language=\"JavaScript\">\r\n<!--\r\n\r\n    function grabCode_905adc75f19a460bbab59d48708d3271() {\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='905adc75f19a460bbab59d48708d3271 ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 905adc75f19a460bbab59d48708d3271';\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 2009 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-->\r\n<\/script><p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\"><br><a href=\"javascript:grabCode_905adc75f19a460bbab59d48708d3271()\"><span style=\"font-size: x-small;        font-style: italic;\">Get \r\n            the MATLAB code \r\n            <noscript>(requires JavaScript)<\/noscript><\/span><\/a><br><br>\r\n      Published with MATLAB&reg; 7.8<br><\/p>\r\n<\/div>\r\n<!--\r\n905adc75f19a460bbab59d48708d3271 ##### SOURCE BEGIN #####\r\n%% A Brief History of polyval\r\n% When I first started working at MathWorks, there were already a few\r\n% functions in MATLAB for working with polynomials.  One of these,\r\n% <https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/ref\/polyval.html |polyval|>,\r\n% for polynomial evaluation, used one algorithm, based on\r\n% <http:\/\/en.wikipedia.org\/wiki\/Horner_scheme Horner's rule>, known for its\r\n% efficiency, especially suitable for embedded controller environments.  \r\n% After a customer wondered why the function took so long to evaluate a\r\n% given polynomial for a scalar input value |x|, we updated |polyval| to\r\n% use a faster algorithm for the scalar case, and that algorithm uses the\r\n% function\r\n% <https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/ref\/filter.html |filter|>.\r\n% I thought I'd discuss a few more details here.\r\n%\r\n% Note: since I wrote this post (but before I published it) there was a\r\n% <http:\/\/view_thread\/253875 discussion>\r\n% on the newsgroup about this same topic.  Be sure to read the reply from\r\n% <https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/authors\/679 |John D'Errico|>.\r\n%% Horner's Rule\r\n% Let's first look at the code for evaluating a polynomial for multiple\r\n% values of input |x|.  Here's the relevant portion of |polyval|, |nc| is\r\n% the length of the vector representing the polynomial (one more than the\r\n% polynomial degree).\r\ndbtype polyval 62:67\r\n%%\r\n% Here's what a general polynomial looks like:\r\n% \r\n% $$ y = p_{1} x^{n} + p_{2} x^{n-1} + ... + p_{n} x + p_{n+1}$$\r\n%\r\n% This equation can be rewritten like this,\r\n%\r\n%%\r\n%\r\n% $$ y = ((...(p_{1}x + p_{2})x +  ... + p_{n-1})x + p_{n})x + p_{n+1}$$\r\n%\r\n% and is the basis for evaluating the polynomial for an array |x|.  Notice\r\n% that the formula basically starts with the coefficient of the highest\r\n% degree (named $p_{1}$) multiplied by |x|, then adds the next coefficient,\r\n% multiplies that result by |x| again, adds the next coefficient, and so\r\n% on, until reaching the final coefficient (the constant term, for $x^0$).\r\n% In this formulation, we never need to explicitly create all the powers of\r\n% |x| to evaluate the polynomial. \r\n%% The Scalar Case, Using filter\r\n% Here's the code for the scalar evalation case using the function\r\n% <https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/ref\/filter.html |filter|>.\r\ndbtype polyval 52:53\r\n%%\r\n% It may not, at first glance, appear relevant to polynomial evaluation.\r\n% However, taking a close look at the mathematics reveals that it indeed\r\n% can be useful.  First, let's examine the |filter| formula.  In this\r\n% version, the coefficients for |a| have been normalized so |a(1)| = 1\r\n% (which is true for the inputs to |filter| you've just seen).\r\n%\r\n%%\r\n% Now, let's try a small example. Suppose we want to evaluate the same\r\n% third degree polynomial I used in another \r\n% <https:\/\/blogs.mathworks.com\/loren\/category\/solving-equations\/ set of posts> \r\n% recently.  \r\np = [1 0 1 -1];\r\n%%\r\n% which represents \r\n%\r\n% $$p(x) = x^{3} + x -1$$\r\n%\r\n%%\r\n% When you call |filter| like this |y = filter(b,a,x)|, the results are\r\n%%\r\n%\r\n% $$ y_{n} = b_{1} x_{n} $$\r\n%\r\n% $$ + b_{2} x_{n-1} + ... + b_{nb+1} x_{n-nb}$$\r\n%\r\n% $$ - a_{2} y_{n-1} + ... + a_{na+1} y_{n-na}$$\r\n%\r\n%%\r\n% Looking again at the |polyval| code, use these substitutions for the\r\n% inputs to |filter|: |b = 1|, |a = [1 -x]|, and |x = p|.  \r\n% We evaluate the filter for the values |p|, the \r\n% polynomial coefficients! Let's do a little math to see why |filter| \r\n% works.  Following the formula for the output of the function |filter|,\r\n% and assuming |p(0) = 0|, we can find the first output.\r\n%\r\n% $$ y_{1} = p_{1} $$\r\n%\r\n% And here are the rest, going up to the final coefficient in our example.\r\n%\r\n% $$ y_{2} = p_{2} + x p{1} $$ \r\n%\r\n% $$ y_{3} = p_{3} + x (p{2} + x p_{1}) $$\r\n%\r\n% $$ y_{3} = p_{3} + x p{2} + x^{2} p_{1} $$\r\n%\r\n% $$ y_{4} = p_{4} + x p{3} + x^{2} p_{2} +x^{3} p_{1} $$\r\n%\r\n% It sure looks like the final answer here is the polynomial |p| evaluated\r\n% at |x|.\r\n%% How Close are the Algorithms?\r\n% Numerically, the two algorithms end up being essentially the same.  \r\nx = 0:0.001:1;\r\nresultArray = polyval(p,x);\r\nresultScalar = zeros(size(x));\r\nfor index = 1:numel(x)\r\n    resultScalar(index) = polyval(p,x(index));\r\nend\r\nallequal = isequal(resultArray,resultScalar)\r\n%%\r\n% Identical!\r\n%% Timing\r\n% I mentioned we made this change for speed.  This particular code is\r\n% tricky to check timing since, especially if the polynomial degree is\r\n% small, the number of mathematic operations can be swamped by the timing\r\n% mechanism.  When computers were slower, and before we increased\r\n% performance of for-loops in MATLAB, measuring the performance differences\r\n% between these code alternatives was much easier than it is today.\r\n%% Multiple Algorithms\r\n% Do you try to write your code as one size fits all, or do you\r\n% create code in which you switch algorithms or implementations, such as\r\n% we've seen in |polyval|, to account for different conditions such as\r\n% array size?  I'd love to hear more about your experiences with these\r\n% ideas <https:\/\/blogs.mathworks.com\/loren\/?p=191#respond here>.\r\n\r\n\r\n##### SOURCE END ##### 905adc75f19a460bbab59d48708d3271\r\n-->","protected":false},"excerpt":{"rendered":"<p>\r\n   \r\n      When I first started working at MathWorks, there were already a few functions in MATLAB for working with polynomials.  One\r\n         of these, polyval, for polynomial evaluation, used... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/loren\/2009\/07\/08\/a-brief-history-of-polyval\/\">read more >><\/a><\/p>","protected":false},"author":39,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[10,15,12],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/191"}],"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=191"}],"version-history":[{"count":1,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/191\/revisions"}],"predecessor-version":[{"id":2568,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/191\/revisions\/2568"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/media?parent=191"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/categories?post=191"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/tags?post=191"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}