{"id":5774,"date":"2015-05-22T09:01:20","date_gmt":"2015-05-22T08:01:20","guid":{"rendered":"http:\/\/www.walkingrandomly.com\/?p=5774"},"modified":"2015-07-02T18:37:43","modified_gmt":"2015-07-02T17:37:43","slug":"matlab-vectorisation-is-a-double-edged-sword","status":"publish","type":"post","link":"https:\/\/walkingrandomly.com\/?p=5774","title":{"rendered":"MATLAB: Vectorisation is a double-edged sword"},"content":{"rendered":"

Update: 2nd July 2015\u00a0<\/strong>The code in github<\/a> has moved on a little since this post was written so I changed the link in the text below to the exact commit that produced the results discussed here<\/a>.<\/p>\n

Imagine that you are a very new MATLAB<\/a> programmer and you have to create an N x N matrix called A where A(i,j) = i+j
\nYour first attempt at a solution might look like this<\/p>\n

N=2000\r\n% Generate a N-by-N matrix where A(i,j) = i + j;\r\nfor ii = 1:N\r\n     for jj = 1:N\r\n         A(ii,jj) = ii + jj;\r\n     end\r\n end<\/pre>\n
On my current machine (Macbook Pro bought in early 2015), the above loop takes 2.03 seconds.<\/strong> You might think that this is a long time for something so simple and complain that MATLAB is slow. The person you complain to points out that you should preallocate your array before assigning to it.<\/p>\n
N=2000\r\n% Generate a N-by-N matrix where A(i,j) = i + j;\r\nA=zeros(N,N);\r\nfor ii = 1:N\r\n     for jj = 1:N\r\n         A(ii,jj) = ii + jj;\r\n     end\r\n end<\/pre>\n
This now takes 0.049 seconds<\/strong> on my machine – a speed up of over 41 times! MATLAB suddenly doesn’t seem so slow after all.<\/p>\n
Word gets around about your problem, however, and seasoned MATLAB-ers see that nested loop, make a funny face twitch and start muttering ‘vectorise, vectorise, vectorise’. Emails soon pile in with vectorised<\/a> solutions<\/p>\n
% Method 1: MESHGRID.\r\n[X, Y] = meshgrid(1:N, 1:N);\r\nA = X + Y;<\/pre>\n
This takes 0.025 seconds<\/strong> on my machine — a healthy speed-up on the loop-with-preallocation solution. You have to understand the meshgrid command, however, in order to understand what’s going on here. It’s still clear (to me at least) what its doing but not as clear as the nice,obvious double loop. Call me old fashioned but I like loops…I understand them.<\/p>\n
% Method 2: Matrix multiplication.\r\nA = (1:N).' * ones(1, N) + ones(N, 1) * (1:N);<\/pre>\n
This one is MUCH harder to read but you don’t worry about it too much because at 0.032 seconds<\/strong> it’s slower than meshgrid.<\/p>\n
% Method 3: REPMAT.\r\nA = repmat(1:N, N, 1) + repmat((1:N).', 1, N);<\/pre>\n
This one appears to be interesting! At 0.009 seconds<\/strong>, it’s the fastest so far – by a healthy amount!<\/p>\n
% Method 4: CUMSUM.\r\nA = cumsum(ones(N)) + cumsum(ones(N), 2);<\/pre>\n
Coming in at 0.052 seconds<\/strong>, this cumsum solution is slower than the preallocated loop.<\/p>\n
% Method 5: BSXFUN.\r\nA = bsxfun(@plus, 1:N, (1:N).');<\/pre>\n
Ahhh, bsxfun<\/a> or ‘The Widow-maker function’ as I sometimes refer to it. Responsible for some of the fastest and most unreadable vectorised MATLAB code I’ve ever written. In this case, it brings execution time down to 0.0045 seconds<\/strong>.<\/p>\n
Whenever I see something that can be vectorised with a repmat, I try to figure out if I can rewrite it as a bsxfun. The result is usually horrible to look at so I tend to keep the original loop commented out above it as an explanation! This particular example isn’t too bad but bsxfun can quickly get hairy.<\/p>\n
Conclusion<\/strong><\/p>\n
Loops in MATLAB aren’t anywhere near as bad as they used to be thanks to advances in JIT compilation<\/a> but it can often pay, speed-wise, to switch to vectorisation. The price you often pay for this speed-up is that vectorised code can become very difficult to read.<\/p>\n
If you’d like the code I ran to get the timings above, it’s on githu<\/a>b (link refers to the exact commit used for this post)\u00a0.\u00a0Here’s the output from the run I referred to in this post.<\/p>\n
Original loop time is 2.025441\r\nPreallocate and loop is 0.048643\r\nMeshgrid time is 0.025277\r\nMatmul time is 0.032069\r\nRepmat time is 0.009030\r\nCumsum time is 0.051966\r\nbsxfun time is 0.004527<\/pre>\n
    \n
  • MATLAB Version: 2015a<\/li>\n
  • Early 2015 Macbook Pro with 16Gb RAM<\/li>\n
  • CPU: 2.8Ghz quad core i7 Haswell<\/li>\n<\/ul>\n
    This post is based on a demonstration given by Mathwork’s Ken Deeley during a recent session at The University of Sheffield.<\/p>\n","protected":false},"excerpt":{"rendered":"
    Update: 2nd July 2015\u00a0The code in github has moved on a little since this post was written so I changed the link in the text below to the exact commit that produced the results discussed here. Imagine that you are a very new MATLAB programmer and you have to create an N x N matrix […]<\/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":true,"jetpack_social_options":{"image_generator_settings":{"template":"highway","default_image_id":0,"font":"","enabled":false},"version":2}},"categories":[53,4,11,7],"tags":[],"class_list":["post-5774","post","type-post","status-publish","format-standard","hentry","category-making-matlab-faster","category-math-software","category-matlab","category-programming"],"jetpack_publicize_connections":[],"jetpack_featured_media_url":"","jetpack_shortlink":"https:\/\/wp.me\/p3swhs-1v8","jetpack_sharing_enabled":true,"_links":{"self":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/5774","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=5774"}],"version-history":[{"count":7,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/5774\/revisions"}],"predecessor-version":[{"id":5805,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/5774\/revisions\/5805"}],"wp:attachment":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=5774"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=5774"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=5774"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}