{"id":4016,"date":"2021-03-16T08:15:51","date_gmt":"2021-03-16T12:15:51","guid":{"rendered":"https:\/\/blogs.mathworks.com\/loren\/?p=4016"},"modified":"2021-03-17T08:49:16","modified_gmt":"2021-03-17T12:49:16","slug":"a-prime-case-study-for-making-matlab-code-go-faster","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/loren\/2021\/03\/16\/a-prime-case-study-for-making-matlab-code-go-faster\/","title":{"rendered":"A Prime Case Study for Making MATLAB Code Go Faster"},"content":{"rendered":"
Today I'd like to welcome a guest blogger, <\/span>Mike Croucher<\/span><\/a>, who recently joined MathWorks as a Customer Success Engineer after a long career spent supporting academic computational research.<\/span><\/div>
Table of Contents<\/span><\/div>
Introduction\r\n<\/span><\/a>Don't Fear the Looper\r\n<\/span><\/a>Changing mod to floor \r\n<\/span><\/a>Trying Integer Types for an Integer Algorithm\r\n<\/span><\/a>MATLAB as a Repository of Algorithms\r\n<\/span><\/a>System Information\r\n<\/span><\/a>Your Need for Speed?<\/span><\/a><\/div><\/div>

Introduction<\/span><\/h3>
Most of my career has been dedicated to collaborating with scientists and engineers to improve their computational workflows. There are many things that could be done including encouraging the use of version control or <\/span>unit testing<\/span><\/a> but, by far, the most frequent request has been<\/span> 'Can you make this thing I made go any faster?' <\/span><\/div>
More than twenty years of being asked this simple question has led me to explore many areas of technical computing such as parallel computing, numerical analysis, <\/span>HPC and cloud computing<\/span><\/a>, <\/span>GPU computing<\/span><\/a>, mathematical optimization,<\/span> automatic differentiation<\/span><\/a> and multiple languages to name just a few.<\/span><\/div>
At its core, the art of making a program faster comes down to iterating over two simple steps.<\/span><\/div>
  1. Find out exactly what is slow.<\/span><\/li>
  2. Try something different to see if it's any faster while still giving the right answer<\/span><\/li><\/ol>
    As you might imagine, there's a world of complexity hidden in the above. 'Trying something different' could mean using language features to optimise the algorithm I've written, trying out parallelisation or even switching to a completely different algorithm althogether. I usually add an additional constraint to point 2 which is 'expend as little effort as possible' because I am fundamentally lazy.<\/span><\/div>
    MATLAB's profiler, which <\/span>recently received an interface overhaul<\/span><\/a> to provide flame graph views of code execution time, helps take care of the first point. Exploring the profiler would provide a rich seam of bloggable content as would exploring all of the <\/span>algorithmic performance enhancements<\/span><\/a> we make from release to release.<\/span><\/div>
    \"\"<\/img><\/div>
    Thanks to <\/span>The MATLAB Execution engine<\/span><\/a> (discussed previously on this blog <\/span>here<\/span><\/a> and <\/span>here<\/span><\/a>) and <\/span>various other optimisations<\/span><\/a>, modern MATLAB is faster than ever and, as shown in the chart below, we are constantly making improvements in this area. <\/span><\/div>
    \"\"<\/img><\/div>
    Of course there are always times when we need to optimise our code for speed and I wanted to share an example that was recently discussed internally.<\/span><\/div>

    Don't Fear the Looper<\/span><\/h3>
    Inspired by this post<\/span><\/a>, a colleague wrote the following in MATLAB to count prime numbers.<\/span><\/div>
    type <\/span>myCountPrimes.m<\/span><\/span><\/div>
    % Count how many primes there are up to and including max_n \r\nfunction res = myCountPrimes(max_n)\r\n res = 0;\r\n for i = 2:max_n\r\n if MyIsPrime(i)\r\n res = res+1;\r\n end\r\n end\r\nend\r\n\r\n\r\n% Function to check if n is prime \r\nfunction result = MyIsPrime(n)\r\n if n <= 1\r\n result = false;\r\n return\r\n end\r\n \r\n q = floor(sqrt(n));\r\n for i = 2:q\r\n if mod(n,i) == 0\r\n result = false;\r\n return\r\n end\r\n end\r\n \r\n result = true;\r\nend<\/div><\/div><\/div><\/div>
    numPrimes = myCountPrimes(1000000)<\/span><\/span><\/div>
    numPrimes = 78498<\/div><\/div><\/div><\/div>
    This code is very heavy on loops - exactly the sort of thing I was trained to avoid in my early MATLAB days but thanks to the Execution Engine, it's not too bad. (See <\/span>this recent post for more on the speed of loops<\/span><\/a> in MATLAB)<\/span><\/div>
    f = @() myCountPrimes(1000000);<\/span><\/span><\/div><\/div>
    originalSpeed = timeit(f)<\/span><\/span><\/div>
    originalSpeed = 3.7878<\/div><\/div><\/div><\/div>
    Even so, there are many things we could try to speed up this particular piece of code. We could take advantage of the fact that every iteration of the loop in <\/span>myCountPrimes<\/span> is independent and use a parallel <\/span>parfor<\/span> loop<\/span><\/a>. This is easy to do but might be heavy on computational resources. Once you realise you can leverage <\/span>parfor<\/span>, it becomes tempting to try to make your problem go away by throwing as many CPU cores as you possible at it. My favourite way to do this these days is to fire up a <\/span>Cloud Center cluster<\/span><\/a> and ramp up the core count as high as makes sense. <\/span><\/div>

    Changing <\/span>mod<\/span> to <\/span>floor<\/span> <\/span><\/h3>
    The profiler showed me that a lot of time was spent in the line <\/span><\/div>
    if mod(n,i) == 0<\/span> <\/span><\/div>
    so one approach is to ask <\/span>'Can we rewrite that line to be any faster?'<\/span>. It's not always obvious which alternatives might work so I often simply figure out a different way to do something and give it a try. An alternative way to write the expression <\/span>mod(n,i)==0<\/span> is <\/span>floor(n\/i) == n\/i<\/span> and it turns out to be much faster! I did this in the <\/span>floorMyCountPrimes<\/span> function, checked that it gives the right result, and demonstrated that it is indeed much faster.<\/span><\/div>
    numPrimes = floorMyCountPrimes(1000000)<\/span><\/span><\/div>
    numPrimes = 78498<\/div><\/div><\/div>
    f = @() floorMyCountPrimes(1000000);<\/span><\/span><\/div><\/div>
    floorSpeed = timeit(f)<\/span><\/span><\/div>
    floorSpeed = 0.5296<\/div><\/div><\/div><\/div>

    Trying Integer Types for an Integer Algorithm<\/span><\/h3>
    By default, all numbers in MATLAB are floating point of type <\/span>double<\/span> which is great for a lot of numerical computing but here we are dealing with a purely integer algorithm. It can sometimes pay dividends to make this explicit in the code and this is one of those times! It turns out that the <\/span>mod<\/span> <\/span>function has to do a lot more work when dealing with floating point numbers than is strictly necessary for this particular example. As the programmer, I know that my algorithm is completely integer based so I can switch to using <\/span>uint32<\/span> for some variables. That is, I change the main loop of <\/span>MyIsPrime<\/span> to <\/span><\/div>
    q = floor(sqrt(n));<\/span><\/span><\/div>
    n = uint32(n);<\/span><\/span><\/div>
    for <\/span>i = uint32(2):q<\/span><\/span><\/div>
    <\/span>if <\/span>mod(n,i) == 0<\/span><\/span><\/div>
    result = false;<\/span><\/span><\/div>
    <\/span>return<\/span><\/span><\/div>
    <\/span>end<\/span><\/span><\/div>
    end<\/span><\/span><\/div><\/div>
    This allows MATLAB to use a specialised and faster version of the <\/span>mod<\/span> function since it now knows that we are dealing with integers instead of doubles.<\/span><\/div>
    I did this in the <\/span>uint32MyIsPrime<\/span> function defined at the end of this post and called it via uint32<\/span>MyCountPrimes. <\/span>First let's check that the result agrees with the original function and then time it<\/span><\/div>
    numPrimes = uint32MyCountPrimes(1000000)<\/span><\/span><\/div>
    numPrimes = 78498<\/div><\/div><\/div>
    f = @() uint32MyCountPrimes(1000000);<\/span><\/span><\/div><\/div>
    uint32Speed = timeit(f)<\/span><\/span><\/div>
    uint32Speed = 0.5660<\/div><\/div><\/div><\/div>
    Using either this approach or the floor-based approach, we achieved a significant performance boost for relatively little work and techniques such as <\/span>parfor<\/span> parallelisation are still on the table if we need even more speed.<\/span><\/div>

    MATLAB as a Repository of Algorithms<\/span><\/h3>
    Finally, we could recognise that the fundamental algorithm we are using to count primes is not optimal and use a different one instead. Taking this to the extreme, we could remember that MATLAB is a repository of numerical algorithms as well as a programming language and so find the one that directly solves the problem we want to solve rather than rolling our own. A quick search of the documentation reveals the <\/span>primes<\/span><\/a> function that returns an array of primes. We produce this array and count the elements as follows.<\/span><\/div>
    numPrimes = numel(primes(1000000))<\/span><\/span><\/div>
    numPrimes = 78498<\/div><\/div><\/div><\/div>
    f = @() numel(primes(1000000));<\/span><\/span><\/div><\/div>
    builtinSpeed = timeit(f)<\/span><\/span><\/div>
    builtinSpeed = 0.0029<\/div><\/div><\/div><\/div>
    More than fast enough for my purposes!<\/span><\/div>

    System Information<\/span><\/h3>
    If you run this script on your own machine, you will almost certainly get different results to me. Using the <\/span>version <\/span>function and the excellent <\/span>CPU Info<\/span><\/a> by Ben Tordoff, here are the details of my machine and MATLAB version.<\/span><\/div>
    versionInfo = [<\/span>'MATLAB Release R' <\/span>version(<\/span>'-release'<\/span>)]<\/span><\/span><\/div>
    versionInfo = <\/span>'MATLAB Release R2021a'<\/div><\/div><\/div><\/div><\/div>
    cpuinfo()<\/span><\/span><\/div>
    ans = struct with fields:<\/span><\/span><\/div>
    Name: 'Intel(R) Core(TM) i7-8650U CPU @ 1.90GHz'\r\n Clock: '2112 MHz'\r\n Cache: '1024 KB'\r\n NumProcessors: 4\r\n OSType: 'Windows'\r\n OSVersion: 'Microsoft Windows 10 Pro'\r\n<\/div><\/div><\/div><\/div><\/div><\/div>

    Your Need for Speed?<\/span><\/h3>
    MATLAB is faster than it has ever been and it will continue to get faster in the future. This, combined with improved hardware, will help you to answer larger and more complex questions than you ever have before. However, there will always be occasions where we need even more. What MATLAB speed-up techniques do you find useful and what do you wish was faster in MATLAB? Let us know <\/span>here<\/span><\/a>.<\/span><\/div>
    Helper functions<\/span><\/div>
    function <\/span>result = floorMyIsPrime(n)<\/span><\/span><\/div><\/div>
    % Function to check if n is prime<\/span><\/span><\/div><\/div>
    if <\/span>n <= 1<\/span><\/span><\/div><\/div>
    result = false;<\/span><\/span><\/div><\/div>
    <\/span>return<\/span><\/span><\/div><\/div>
    end<\/span><\/span><\/div><\/div>
    <\/div><\/div>
    q = floor(sqrt(n));<\/span><\/span><\/div><\/div>
    for <\/span>i = 2:q<\/span><\/span><\/div><\/div>
    <\/span>if <\/span>floor(n\/i) == n\/i<\/span><\/span><\/div><\/div>
    result = false;<\/span><\/span><\/div><\/div>
    <\/span>return<\/span><\/span><\/div><\/div>
    <\/span>end<\/span><\/span><\/div><\/div>
    end<\/span><\/span><\/div><\/div>
    <\/div><\/div>
    result = true;<\/span><\/span><\/div><\/div>
    end<\/span><\/span><\/div><\/div>
    <\/div><\/div>
    function <\/span>res = floorMyCountPrimes(max_n)<\/span><\/span><\/div><\/div>
    % Count how many primes there are up to and including max_n<\/span><\/span><\/div><\/div>
    res = 0;<\/span><\/span><\/div><\/div>
    for <\/span>i = 2:max_n<\/span><\/span><\/div><\/div>
    <\/span>if <\/span>floorMyIsPrime(i)<\/span><\/span><\/div><\/div>
    res = res+1;<\/span><\/span><\/div><\/div>
    <\/span>end<\/span><\/span><\/div><\/div>
    end<\/span><\/span><\/div><\/div>
    end<\/span><\/span><\/div><\/div>
    <\/div><\/div>
    function <\/span>result = uint32MyIsPrime(n)<\/span><\/span><\/div><\/div>
    % Function to check if n is prime<\/span><\/span><\/div><\/div>
    if <\/span>n <= 1<\/span><\/span><\/div><\/div>
    result = false;<\/span><\/span><\/div><\/div>
    <\/span>return<\/span><\/span><\/div><\/div>
    end<\/span><\/span><\/div><\/div>
    <\/div><\/div>
    q = floor(sqrt(n));<\/span><\/span><\/div><\/div>
    n = uint32(n);<\/span><\/span><\/div><\/div>
    for <\/span>i = uint32(2):q<\/span><\/span><\/div><\/div>
    <\/span>if <\/span>mod(n,i) == 0<\/span><\/span><\/div><\/div>
    result = false;<\/span><\/span><\/div><\/div>
    <\/span>return<\/span><\/span><\/div><\/div>
    <\/span>end<\/span><\/span><\/div><\/div>
    end<\/span><\/span><\/div><\/div>
    <\/div><\/div>
    result = true;<\/span><\/span><\/div><\/div>
    end<\/span><\/span><\/div><\/div>
    <\/div><\/div>
    function <\/span>res = uint32MyCountPrimes(max_n)<\/span><\/span><\/div><\/div>
    % Count how many primes there are up to and including max_n<\/span><\/span><\/div><\/div>
    res = 0;<\/span><\/span><\/div><\/div>
    for <\/span>i = 2:max_n<\/span><\/span><\/div><\/div>
    <\/span>if <\/span>uint32MyIsPrime(i)<\/span><\/span><\/div><\/div>
    res = res+1;<\/span><\/span><\/div><\/div>
    <\/span>end<\/span><\/span><\/div><\/div>
    end<\/span><\/span><\/div><\/div>
    end<\/span><\/span><\/div><\/div><\/div>
    Copyright 2021 The MathWorks, Inc.<\/span><\/div>\r\n<\/div>