{"id":1552,"date":"2009-08-03T12:08:15","date_gmt":"2009-08-03T11:08:15","guid":{"rendered":"http:\/\/www.walkingrandomly.com\/?p=1552"},"modified":"2009-08-03T12:08:15","modified_gmt":"2009-08-03T11:08:15","slug":"a-faster-version-of-matlabs-interp1","status":"publish","type":"post","link":"https:\/\/walkingrandomly.com\/?p=1552","title":{"rendered":"A faster version of MATLAB&#8217;s interp1"},"content":{"rendered":"<p>I often get sent a piece of code written in something like <a href=\"http:\/\/www.mathworks.co.uk\/\">MATLAB<\/a> or <a href=\"http:\/\/www.wolfram.com\/\">Mathematica<\/a> and get asked &#8216;<strong>how can I make this faster?<\/strong>&#8216;.\u00a0 I have no doubt that some of you are thinking that the correct answer is something like &#8216;<strong>completely rewrite it in Fortran or C<\/strong>&#8216; and if you are then I can only assume that you have never been face to face with a harassed, sleep-deprived researcher who has been sleeping in the lab for 7 days straight in order to get some results prepared for an upcoming conference.<\/p>\n<p>Telling such a person that they should throw their painstakingly put-together piece of code away and learn a considerably more difficult language in order to rewrite the whole thing is unhelpful at best and likely to result in the loss of an eye at worst.\u00a0 I have found that a MUCH better course of action is to profile the code, find out where the slow bit is and then do what you can with that.\u00a0 Often you can get very big speedups in return for only a modest amount of work and as a bonus you get to keep the use of both of your eyes.<\/p>\n<p>Recently I had an email from someone who had profiled her MATLAB code and had found that it was spending a lot of time in the <strong><a href=\"http:\/\/www.mathworks.com\/access\/helpdesk\/help\/techdoc\/index.html?\/access\/helpdesk\/help\/techdoc\/ref\/interp1.html\">interp1<\/a> <\/strong>function.\u00a0 Would it be possible to use the <a href=\"http:\/\/www.nag.co.uk\/numeric\/MB\/start.asp\">NAG toolbox for MATLAB<\/a> (which hooks into NAG&#8217;s superfast Fortran library) to get her results faster she wondered?\u00a0 Between the two of us and the NAG technical support team we eventually discovered that the answer was a very definite <strong>yes.<\/strong><\/p>\n<p>Rather than use her actual code, which is rather complicated and domain-specific, let&#8217;s take a look at a highly simplified version of her problem.\u00a0 Let&#8217;s say that you have the following data-set<\/p>\n<pre> x = [0;0.2;0.4;0.6;0.75;0.9;1];\r\n y = [1;1.22140275816017;1.49182469764127;1.822118800390509;2.117000016612675;2.45960311115695; 2.718281828459045];<\/pre>\n<p style=\"text-align: center;\"><img decoding=\"async\" class=\"aligncenter\" src=\"\/images\/matlab\/interp1_data.png\" alt=\"interp1 data\" \/><\/p>\n<p>and you want to interpolate this data on a fine grid between x=0 and x=1 using piecewise cubic Hermite interpolation.  One way of doing this in MATLAB is to use the <strong>interp1<\/strong> function as follows<\/p>\n<pre>tic\r\nt=0:0.00005:1;\r\ny1 = interp1(x,y,t,'pchip');\r\ntoc<\/pre>\n<p style=\"text-align: center;\"><img decoding=\"async\" class=\"aligncenter\" src=\"\/images\/matlab\/interp1_data1.png\" alt=\"interp1 data\" \/><!--<\/p--><\/p>\n<p>The tic and toc statements are there simply to provide timing information and on my system the above code took  0.503 seconds.  We can do exactly the same calculation using the NAG toolbox for MATLAB:<\/p>\n<pre>tic\r\n[d,ifail] = e01be(x,y);\r\n[y2, ifail] = e01bf(x,y,d,t);\r\ntoc<\/pre>\n<p>which took 0.054958 seconds on my system making it around <strong>9 times faster<\/strong> than the pure MATLAB solution &#8211; not bad at all for such a tiny amount of work.\u00a0 Of course y1 and y2 are identical as the following call to norm demonstrates<\/p>\n<pre>norm(y1-y2)\r\nans =\r\n      0<\/pre>\n<p>Of course the real code contained a lot more than just a call to interp1 but using the above technique decreased the run time of the user&#8217;s application from 1 hour 10 minutes down to only 26 minutes.<\/p>\n<p>Thanks to the NAG technical support team for their assistance with this particular problem and to the postgraduate student who came up with the idea in the first place.<\/p>\n","protected":false},"excerpt":{"rendered":"<p>I often get sent a piece of code written in something like MATLAB or Mathematica and get asked &#8216;how can I make this faster?&#8216;.\u00a0 I have no doubt that some of you are thinking that the correct answer is something like &#8216;completely rewrite it in Fortran or C&#8216; and if you are then I can [&hellip;]<\/p>\n","protected":false},"author":3,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"jetpack_post_was_ever_published":false,"footnotes":"","jetpack_publicize_message":"","jetpack_publicize_feature_enabled":true,"jetpack_social_post_already_shared":false,"jetpack_social_options":{"image_generator_settings":{"template":"highway","default_image_id":0,"font":"","enabled":false},"version":2}},"categories":[4,11,28,7],"tags":[],"class_list":["post-1552","post","type-post","status-publish","format-standard","hentry","category-math-software","category-matlab","category-nag-library","category-programming"],"jetpack_publicize_connections":[],"jetpack_featured_media_url":"","jetpack_shortlink":"https:\/\/wp.me\/p3swhs-p2","jetpack_sharing_enabled":true,"_links":{"self":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/1552","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/users\/3"}],"replies":[{"embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=1552"}],"version-history":[{"count":8,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/1552\/revisions"}],"predecessor-version":[{"id":1562,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/1552\/revisions\/1562"}],"wp:attachment":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=1552"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=1552"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=1552"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}