{"id":1488,"date":"2010-09-21T19:27:27","date_gmt":"2010-09-21T18:27:27","guid":{"rendered":"http:\/\/www.walkingrandomly.com\/?p=1488"},"modified":"2010-09-22T02:16:35","modified_gmt":"2010-09-22T01:16:35","slug":"a-faster-version-of-matlabs-fsolve-using-the-nag-toolbox-for-matlab","status":"publish","type":"post","link":"https:\/\/walkingrandomly.com\/?p=1488","title":{"rendered":"A faster version of MATLAB’s fsolve using the NAG Toolbox for MATLAB"},"content":{"rendered":"

Part of my job at the University of Manchester<\/a> is to help researchers improve their code in systems such as Mathematica<\/a>, MATLAB<\/a>, NAG<\/a> and so on.\u00a0 One of the most common questions I get asked is ‘Can you make this go any faster please?’ <\/em>and it always makes my day if I manage to do something significant such as a speed up of a factor of 10 or more.\u00a0 The users, however, are often happy with a lot less and I once got bought a bottle of (rather good) wine for a mere 25% speed-up (Note:\u00a0 Gifts of wine aren’t mandatory)<\/p>\n

I employ numerous tricks to get these speed-ups including applying vectorisation, using mex files, modifying the algorithm to something more efficient,picking the brains of colleagues and so on.\u00a0 Over the last year or so, I have managed to help several researchers get significant speed-ups in their MATLAB programs by doing one simple thing – swapping the fsolve<\/a> function for an equivalent from the NAG Toolbox for MATLAB<\/a>.<\/p>\n

The fsolve function is part of MATLAB’s optimisation toolbox<\/a> and, according to the documentation, it does the following:<\/p>\n

“fsolve Solves a system of nonlinear equation of the form F(X) = 0 where F and X may be vectors or matrices.”<\/em><\/p>\n

The NAG Toolbox for MATLAB contains a total of 6 different routines that can perform similar calculations to fsolve.\u00a0 These 6 routines are split into 2 sets of 3 as follows<\/p>\n

For when you have function values only:<\/strong>
\nc05nb \u2013 Solution of system of nonlinear equations using function values only (easy-to-use)
\nc05nc \u2013 Solution of system of nonlinear equations using function values only (comprehensive)
\nc05nd \u2013 Solution of system of nonlinear equations using function values only (reverse communication)
\n <\/strong><\/p>\n

For when you have both function values and first derivatives<\/strong>
\nc05pb \u2013 Solution of system of nonlinear equations using first derivatives (easy-to-use)
\nc05pc \u2013 Solution of system of nonlinear equations using first derivatives (comprehensive)
\nc05pd \u2013 Solution of system of nonlinear equations using first derivatives (reverse communication)<\/p>\n

So, the NAG routine you choose depends on whether or not you can supply first derivatives and exactly which options you’d like to exercise.\u00a0 For many problems it will come down to a choice between the two ‘easy to use’ routines – c05nb and c05pb (at least, they are the ones I’ve used most of the time)<\/p>\n

Let’s look at a simple example.\u00a0 You’d like to find a root of the following system of non-linear equations.<\/p>\n

F(1)=exp(-x(1)) + sinh(2*x(2)) + tanh(2*x(3)) – 5.01;
\nF(2)=exp(2*x(1)) + sinh(-x(2) ) + tanh(2*x(3) ) – 5.85;
\nF(3)=exp(2*x(1)) + sinh(2*x(2) ) + tanh(-x(3) ) – 8.88;<\/p>\n

i.e. you want to find a vector X containing three values for which F(X)=0.<\/p>\n

To solve this using MATLAB and the optimisation toolbox you could proceed as follows, first create a .m file called fsolve_obj_MATLAB.m<\/a> (henceforth referred to as the objective function) that contains the following<\/p>\n

function  F=fsolve_obj_MATLAB(x)\r\nF=zeros(1,3);\r\nF(1)=exp(-x(1))  + sinh(2*x(2)) + tanh(2*x(3)) - 5.01;\r\nF(2)=exp(2*x(1)) +  sinh(-x(2) ) + tanh(2*x(3) ) - 5.85;\r\nF(3)=exp(2*x(1)) +  sinh(2*x(2) ) + tanh(-x(3) ) - 8.88;\r\nend<\/pre>\n
Now type the following at the MATLAB  command prompt to actually solve the problem:<\/p>\n
options=optimset('Display','off'); %Prevents fsolve from  displaying too much information\r\nstartX=[0 0 0]; %Our  starting guess for the solution vector\r\nX=fsolve(@fsolve_obj_MATLAB,startX,options)<\/pre>\n
MATLAB finds the following solution (Note that it’s only found one of the solutions, not all of them, which is typical of the type of algorithm used in fsolve)<\/p>\n
X =\r\n0.9001\u00a0\u00a0\u00a0  1.0002\u00a0\u00a0\u00a0 1.0945<\/pre>\n
Since we are not supplying the derivatives of our objective function and we are not using any special options, it turns out to be very easy to switch from using the optimisation toolbox to the NAG toolbox for this particular problem by making use of the routine c05nb.<\/p>\n
First, we need  to change the header of our objective function’s .m file from<\/p>\n
function F=fsolve_obj_MATLAB(x)<\/pre>\n
to<\/p>\n
function  [F,iflag]=fsolve_obj_NAG(n,x,iflag)<\/pre>\n
so our new .m file, fsolve_obj_NAG.m<\/a>, will contain the following<\/p>\n
function  [F,iflag]=fsolve_obj_NAG(n,x,iflag)\r\nF=zeros(1,3);\r\nF(1)=exp(-x(1)) + sinh(2*x(2)) + tanh(2*x(3)) - 5.01;\r\nF(2)=exp(2*x(1)) + sinh(-x(2) ) + tanh(2*x(3) ) - 5.85;\r\nF(3)=exp(2*x(1)) + sinh(2*x(2) ) + tanh(-x(3) ) - 8.88;\r\nend<\/pre>\n
Note that the ONLY change we made to the objective function was the very first line.\u00a0 The  NAG routine c05nb expects to see some extra arguments in the objective function (iflag and n in this example) and it expects to see them in a very particular order but we are not obliged to use them.\u00a0 Using this modified objective function we can proceed as follows.<\/p>\n
startX=[0 0 0]; %Our starting guess\r\nX=c05nb('fsolve_obj_NAG',startX)<\/pre>\n
NAG gives the same result as the optimisation toolbox:<\/p>\n
X =\r\n0.9001\u00a0\u00a0\u00a0 1.0002\u00a0\u00a0\u00a0 1.0945<\/pre>\n
Let’s look at timings.\u00a0 I performed the calculations above on an 800Mhz laptop running Ubuntu Linux 10.04, MATLAB 2009 and Mark 22 of the NAG toolbox and got the following timings  (averaged over 10 runs):<\/p>\n
Optimisation  toolbox: 0.0414 seconds
\nNAG toolbox: 0.0021 seconds<\/p>\n
So, for this particular problem NAG was faster than MATLAB by a factor of 19.7 times on my machine.<\/p>\n
That’s all well and good, but this is a simple problem.  Does this scale to more realistic problems I hear you ask?  Well, the answer is ‘yes’.  I’ve used this technique several times now on production code and it almost always results in some kind of speed-up.  Not always 19.7 times faster, I’ll grant you, but usually enough to make the modifications worthwhile.<\/p>\n
I can’t actually post any of the more complicated examples here because the code in question belongs to other people but I was recently sent a piece of code from a researcher that made heavy use of fsolve and a typical run took 18.19 seconds (that’s for everything, not just the fsolve bit).\u00a0 Simply by swapping a couple of lines of code to make it  use the NAG toolbox rather than the optimisation toolbox I reduced the runtime to 1.08 seconds – a speed-up of almost 17 times.<\/p>\n
Sadly, this technique doesn’t always result in a speed-up and I’m working on figuring out exactly when the benefits disappear.  I guess that the main reason for NAG’s good performance is that it uses highly optimised, compiled Fortran code compared to MATLAB’s interpreted .m code.  One case where NAG didn’t help was for a massively complicated objective function and the majority of the run-time of the code was spent evaluating this function.  In this situation, NAG and MATLAB came out roughly neck and neck.<\/p>\n
In the meantime, if you have some code that you’d like me to try it on then drop me a line.<\/p>\n
If you enjoyed this post then you may also like:<\/strong><\/p>\n