{"id":2907,"date":"2010-10-05T11:40:52","date_gmt":"2010-10-05T10:40:52","guid":{"rendered":"http:\/\/www.walkingrandomly.com\/?p=2907"},"modified":"2010-10-05T12:44:45","modified_gmt":"2010-10-05T11:44:45","slug":"using-nag-toolbox-to-speed-up-matlabs-fsolve-part-2","status":"publish","type":"post","link":"https:\/\/walkingrandomly.com\/?p=2907","title":{"rendered":"Using NAG Toolbox to speed up MATLAB’s fsolve part 2"},"content":{"rendered":"

Not long after publishing my article, A faster version of MATLAB\u2019s fsolve using the NAG Toolbox for MATLAB<\/a>, I was contacted by Biao Guo of Quantitative Finance Collector<\/a> who asked me what I could do with the following piece of code.<\/p>\n

function MarkowitzWeight = fsolve_markowitz_MATLAB()\r\n\r\nRandStream.setDefaultStream (RandStream('mt19937ar','seed',1));\r\n% simulate 500 stocks return\r\nSimR = randn(1000,500);\r\nr = mean(SimR); % use mean return as expected return\r\ntargetR = 0.02;\r\nrCOV = cov(SimR); % covariance matrix\r\nNumAsset = length(r);\r\n\r\noptions=optimset('Display','off');\r\nstartX = [1\/NumAsset*ones(1, NumAsset), 1, 1];\r\nx = fsolve(@fsolve_markowitz,startX,options, r, targetR, rCOV);\r\nMarkowitzWeight = x(1:NumAsset);\r\nend\r\n\r\nfunction F = fsolve_markowitz(x, r, targetR, rCOV)\r\nNumAsset = length(r);\r\nF = zeros(1,NumAsset+2);\r\nweight = x(1:NumAsset); % for asset weights\r\nlambda = x(NumAsset+1);\r\nmu = x(NumAsset+2);\r\nfor i = 1:NumAsset\r\n    F(i) = weight(i)*sum(rCOV(i,:))-lambda*r(i)-mu;\r\nend\r\nF(NumAsset+1) = sum(weight.*r)-targetR;\r\nF(NumAsset+2) = sum(weight)-1;\r\nend\r\n<\/pre>\n
This is an example from financial mathematics and Biao explains it in more detail in a post over at his blog<\/a>, mathfinance.cn. Now, it isn’t particularly computationally challenging since it only takes 6.5 seconds to run on my 2Ghz dual core laptop.  It is, however, significantly more challenging than the toy problem I discussed in my original fsolve post.  So, let’s see how much we can optimise it using NAG and any other techniques that come to mind.<\/p>\n
Before going on, I’d just like to point out that we intentionally set the seed of the random number generator with<\/p>\n
RandStream.setDefaultStream (RandStream('mt19937ar','seed',1));\r\n<\/pre>\n
to ensure we are always operating with the same data set.  This is purely for benchmarking purposes.<\/p>\n
Optimisation trick 1 – replacing fsolve with NAG’s c05nb<\/h3>\n
The first thing I had to do to Biao’s code in order to apply the NAG toolbox was to split it into two files.  We have the main program, fsolve_markowitz_NAG.m<\/a><\/p>\n
function MarkowitzWeight = fsolve_markowitz_NAG()\r\nglobal r targetR rCOV;\r\n%seed random generator to ensure same results every time for benchmark purposes\r\nRandStream.setDefaultStream (RandStream('mt19937ar','seed',1));\r\n\r\n% simulate 500 stocks return\r\nSimR = randn(1000,500);\r\nr = mean(SimR)'; % use mean return as expected return\r\ntargetR = 0.02;\r\nrCOV = cov(SimR); % covariance matrix\r\nNumAsset = length(r);\r\n\r\nstartX = [1\/NumAsset*ones(1,NumAsset), 1, 1];\r\ntic\r\nx = c05nb('fsolve_obj_NAG',startX);\r\ntoc\r\nMarkowitzWeight = x(1:NumAsset);\r\nend\r\n<\/pre>\n
and the objective function, fsolve_obj_NAG.m<\/a><\/p>\n
function [F,iflag] = fsolve_obj_NAG(n,x,iflag)\r\nglobal r targetR rCOV;\r\nNumAsset = length(r);\r\nF = zeros(1,NumAsset+2);\r\nweight = x(1:NumAsset); % for asset weights\r\nlambda = x(NumAsset+1);\r\nmu = x(NumAsset+2);\r\nfor i = 1:NumAsset\r\n    F(i) = weight(i)*sum(rCOV(i,:))-lambda*r(i)-mu;\r\nend\r\nF(NumAsset+1) = sum(weight.*r)-targetR;\r\nF(NumAsset+2) = sum(weight)-1;\r\nend\r\n<\/pre>\n
I had to put the objective function into its own .m file because because the NAG toolbox currently doesn’t accept function handles (This may change in a future release I am reliably told).  This is a pain but not a show stopper.<\/p>\n
Note also that the function header for the NAG version of the objective function is different from the pure MATLAB version as
\nper my original article<\/a>.<\/p>\n
Transposition problems<\/strong>
\nThe next thing to notice is that I have done the following in the main program, fsolve_markowitz_NAG.m<\/strong><\/p>\n
r = mean(SimR)';\r\n<\/pre>\n
compared to the original<\/p>\n
r = mean(SimR);\r\n<\/pre>\n
This is because the NAG routine, c05nb, does something rather bizarre.  I discovered that if you pass a 1xN matrix to NAGs c05nb then it gets transposed in the objective function to Nx1.  Since Biao relies on this being 1xN when he does<\/p>\n
F(NumAsset+1) = sum(weight.*r)-targetR;\r\n<\/pre>\n
we get an error unless you take the transpose of r in the main program.  I consider this to be a bug in the NAG Toolbox and it is currently with NAG technical support.<\/p>\n
Global variables<\/strong>
\nSadly, there is no easy way to pass extra variables to the objective function when using NAG’s c05nb.  In MATLAB Biao could do this<\/p>\n
x = fsolve(@fsolve_markowitz,startX,options, r, targetR, rCOV);\r\n<\/pre>\n
No such luck with c05nb.  The only way to get around it is to use global variables (Well, a small lie, there is another way, you can use reverse communication in a different NAG function, c05nd, but that’s a bit more complicated and I’ll leave it for another time I think)<\/p>\n
So, It’s pretty clear that the NAG function isn’t as easy to use as MATLAB’s fsolve function but is it any faster? The following is typical.<\/p>\n
>> tic;fsolve_markowitz_NAG;toc\r\nElapsed time is 2.324291 seconds.\r\n<\/pre>\n
So, we have sped up our overall computation time by a factor of 3.  Not bad, but nowhere near the 18 times speed-up that I demonstrated in my original post.<\/p>\n
Optimisation trick 2 – Don’t sum so much<\/h3>\n
In his recent blog post, Biao challenged me to make his code almost 20 times faster and after applying my NAG Toolbox trick I had ‘only’ managed to make it 3 times faster.  So, in order to see what else I might be able to do for him I ran his code through the MATLAB profiler and discovered that it spent the vast majority of its time in the following section of the objective function<\/p>\n
for i = 1:NumAsset\r\n    F(i) = weight(i)*sum(rCOV(i,:))-lambda*r(i)-mu;\r\nend\r\n<\/pre>\n
It seemed to me that the original code was calculating the sum of rCOV too many times.  On a small number of assets (50 say) you’ll probably not notice it but with 500 assets, as we have here, it’s very wasteful.<\/p>\n
So, I changed the above to<\/p>\n
for i = 1:NumAsset\r\n    F(i) = weight(i)*sumrCOV(i)-lambda*r(i)-mu;\r\nend\r\n<\/pre>\n
So,instead of rCOV, the objective function needs to be passed the following in the main program<\/p>\n
sumrCOV=sum(rCOV)\r\n<\/pre>\n
We do this summation once and once only and make a big time saving.\u00a0 This is done in the files fsolve_markowitz_NAGnew.m<\/a> and fsolve_obj_NAGnew.m<\/a><\/p>\n
With this optimisation, along with using NAG’s c05nb we now get the calculation time down to<\/p>\n
>> tic;fsolve_markowitz_NAGnew;toc\r\nElapsed time is 0.203243 seconds.\r\n<\/pre>\n
Compared to the 6.5 seconds of the original code, this represents a speed-up of over 32 times<\/strong>.  More than enough to satisfy Biao’s challenge.<\/p>\n
Checking the answers<\/h3>\n
Let’s make sure that the results agree with each other.  I don’t expect them to be exactly equal due to round off errors etc but they should be extremely close.  I check that the difference between each result is less than 1e-7:<\/p>\n
original_results = fsolve_markowitz_MATLAB;\r\nnagver1_results = fsolve_markowitz_NAGnew;\r\n>> all(abs(original_results - nagver1_results)<1e8)\r\nans =\r\n     1\r\n<\/pre>\n
That’ll do for me!<\/p>\n
Summary<\/h3>\n
The current iteration of the NAG Toolbox for MATLAB (Mark 22) may not be as easy to use as native MATLAB toolbox functions (for this example at least) but it can often provide useful speed-ups and in this example, NAG gave us a factor of 3<\/strong> improvement compared to using fsolve.<\/p>\n
For this particular example, however, the more impressive speed-up came from pushing the code through the MATLAB profiler, identifying the hotspot and then doing something about it.\u00a0 There are almost certainly more optimisations to make but with the code running at less than a quarter of a second on my modest little laptop I think we will leave it there for now.<\/p>\n