# Loren on the Art of MATLAB

## Timing Code

and Revisiting Fibonacci Festival from 2006

I originally wrote about Fibonacci numbers in 2006. Today, along with Sarah Wait Zaranek, we'd like to discuss some of those same ideas, but focus only on execution time. In R2013b, MATLAB includes a new timing function, timeit, one that I mentioned in some earlier posts when it was only available on the MATLAB File Exchange.

### Contents

#### Using timeit

To use timeit, you simply pass in a function handle that requires no inputs. If your function requires input values, you can pass them in by creating an anonymous function, where the values are captured from the current workspace. You can learn more about function handles in the MATLAB documentation, and find examples and discussions on this blog.

timeit runs the code several times to ensure that any first-run timing effects do not affect the results - and returns the median time. If you wish to time your function in the situation where it returns more than 0 or 1 output values, you can specify that as one of the inputs to timeit.

#### Four Algorithms

In the earlier Fibonacci post, I introduce and times four algorithms. Today we will use the same algorithms, coded as separate programs, in conjunction with timeit. Since we have four algorithms, we will initialize a vector for collecting the timing information.

fibTimes = zeros(1,4);


#### Straightforward for-loop

Let's calculate 102 Fibonacci values (since the first 2 are known already, [1, 1]).

nf = 102;
type fibloop

function fibs = fibloop(nf)
fibs(1) = 1;
fibs(2) = 1;
for n = 3:nf
fibs(n) = fibs(n-1)+fibs(n-2);
end


To use timeit, we simply give it a function handle to the function we wish to time.

fibTimes(1) = timeit(@()fibloop(nf));


#### for-loop with Preallocation

You've all heard it before. One reason the former calculation took so much time is at least alleged to be because of growing the array fibf each time through the loop. So now let's use the same loop, but preallocate the output array so it doesn't need to grow during the calculation.

type fibloopPrealloc

function fibs = fibloopPrealloc(nf)
fibs = ones(1,nf);
for n = 3:nf
fibs(n) = fibs(n-1) + fibs(n-2);
end



Now time it.

fibTimes(2) = timeit(@()fibloopPrealloc(nf));


#### Recursive Function

We can use recursion to calculate successive terms, as I showed in that earlier post. Here's the code.

type fibRecursive

function fibs = fibRecursive(n)
if n == 1,
fibs = 1;   % First element is 1.
return;
elseif n == 2
fibs = [1 1];  % First two elements are 1.
else
% Call fibrec with previous result and compute next one from it.
fm1 = fibRecursive(n-1);
fibs = [fm1 fm1(end)+fm1(end-1)];
end


And now let's use it to calculate the first chunk of Fibonacci numbers.

fibTimes(3) = timeit(@()fibRecursive(nf));


#### Use filter to Calculate the Sequence

I also explained in the earlier Fibonacci post how to use filter to calculate the sequence as well. Here's the code.

type fibFilter
fibTimes(4) = timeit(@()fibFilter(nf));

function fibf = fibFilter(n)
x = [1 zeros(1,n-1)];
a = [1 -1 -1];
b = 1;
fibf = filter(b, a, x);


#### Compare Computation Times

It's time to see how quickly the various methods performed. We have multiplied the times by 1000 to make reading the numbers easier.

fprintf('%s\n\n','Compare Times (in milliseconds)')
meths = ['loop           ',...
'prealloc       ',...
'recursive      ',...
'filter         '];
fprintf('%s\n\n',meths)
fprintf('%2.9f    ', fibTimes*1000), fprintf('\n')

Compare Times (in milliseconds)

loop           prealloc       recursive      filter

0.040766819    0.010003257    0.517234655    0.010549081


#### Do You Do Lots of Timing?

Do you do lots of timing? How do you do it? Profiler? Will timeit help? Let us know here.

Get the MATLAB code

Published with MATLAB® R2013b

### 12 Responses to “Timing Code”

1. Eric replied on :

Hi Loren! Am I correct in thinking that fibRecursive is a decent demonstration of the current function call/resolution overhead? It seems to me that the JIT has improved so much over the years that the overhead of function call resolution(?) has become more prominent than it used to be.

Have you folks considered implementing some sort of optional annotations the user could specify in their function calls to make the resolution more static-y in order to reduce the time required for MATLAB to figure out which function to call? I’m thinking of something like having the ability to add annotations like ‘pwd\myFunc()’, ‘pwd\subFolder\myFunc()’, ‘nested\nestedFunc()’, ‘local\localFunc()’, etc. This would also help prevent variable shadowing functions. You could even have an annotation like ‘fn\anyFunc()’ to tell MATLAB to use the current function precedence resolution rules, except that it can skip the ‘is anyFunc a variable’ check…

2. Eric replied on :

3. Loren Shure replied on :

Eric-

fibRecursive certainly has relatively high function call overhead, but it is also growing an array each time (except for n=1,2).

We also concluded that notation such as you suggested initially would be a poor idea.

–Loren

4. Eric replied on :

Thanks Loren, good point about the array growing.

Are there any thoughts at TMW about ways to reduce the existing function call overhead? I’m sure it must be optimized as much as possible using the current language specification, so that’s why I was wondering about ways that the user could choose to explicitly reduce the dynamicness-related overhead for function calls that they know are in a critical section / hot loop. It just seems to me that you can only optimize things so far, at some point the precedence rules means that there has to be X lookups in the dispatch process unless the user can choose to give guidance to the engine, or am I wrong here? Thanks!

I have a folder on my Matlab path littered with quick scripts comparing the speed of two or more implementations of something I am doing in a bigger demo. These tend to be a bit sporadic though and I always worry I might not be taking something into account that will bias one method or other.

Generally these just take the form of a tic toc wrapping up a for loop that runs the function a sufficient number of times for me to trust the results (even then I run the script multiple times and often swap the order in which it compares the two methods in case I accidentally left something behind from the first method that would bias the second).

timeit might well be useful for this if it wraps up all that in a more predictable and reliable manner.

I also use the profiler which is pleasingly easy to use compared with when I try to profile C++ code. I have only looked at this a little though and used it to pin down my most serious bottlenecks in algorithms that are massively iterated.

6. Knut replied on :

Quite astonishing that a for-loop is (slightly) faster than filter().

Assuming that it is the recursion that slows down filter (not being able to use SIMD), I tested FIR-filtering with [1 1]. Turns out that the loop is slightly faster than filter() for that case as well.

So is the JIT stuff really smart, or is filter() less efficient than one might expect?

Turns out that both conv() and filter([1 1], 1) are significantly faster than doing the same in a loop, but only for large vectors (40000 elements). And once you throw in the recursion, writing out the loop is just as fast as filter, even for long vectors.

7. Loren Shure replied on :

Eric-

Right now, I have told you all I know about things to consider. Among these is use of function handles. But there isn’t one answer that will work best in all cases.

Thanks for sharing your work methods. timeit should really help you not have to run your timing in loops – that’s part of what it does for you to shake out the startup costs.

Knut-

As you can see, which method is best does depend on size.

Thanks for all the good thoughts, everyone.

–Loren

8. Eric replied on :

Thanks Loren, I hadn’t thought about the use of function handles in this light, I’ll have to do some testing :)

9. Matt Tearle replied on :

One thing I used to do a lot when timing code was compare not just two algorithms, but how they scaled with a size parameter, such as the effect Knut noticed. So timeit could make that kind of testing code neater:

sz = [20,50,100,200,500,1000,2000];
n = length(sz);
tm = zeros(n,1);
for k = 1:n
tm(k) = timeit(@()fibFilter(sz(k)));
end
plot(sz,tm,’o-’)

Anonymous function handles FTW!

10. Yair Altman replied on :

For large values of nf, a recursive implementation that uses memoization (http://en.wikipedia.org/wiki/Memoization) significantly speeds up the performance. On my system this runs fastest of all for large-enough nf (of course it’s limited by Matlab’s max recursion depth):

function value = fibonacci_memoization(n)
persistent cachedData
if isempty(cachedData)
cachedData(1) = 1;
cachedData(2) = 1;
end
try
value = cachedData(n);
catch
%value = fibonacci(n-2) + fibonacci(n-1);
if length(cachedData) < n-1
% populate the cache for n-1,n-2; discard result
fibonacci_memoization(n-1);
end
value = cachedData(n-2) + cachedData(n-1);
cachedData(n) = value; % update the cache
end
end % fibonacci_memoization

Readers might also be interested in comparing John D’Errico implementation of the Fibonacci function (http://www.mathworks.com/matlabcentral/fileexchange/34766-the-fibonacci-sequence). The source code may seem complex, but the resulting performance gain is well worth the extra complexity. I believe that readers who will read this utility’s source code and understand its underlying logic will gain insight into several performance tricks that could be very useful in general.

11. Yair Altman replied on :

…and of course, the best performance is achieved using the mathematical identity of fibonacci’s sequence values:

function values = fibonacci(n)
s5 = sqrt(5); % can be cached for further optimization
Psi = (1+s5)/2; % can be cached for further optimization
value = round((Psi.^n – (-Psi).^-n) / s5);
end

12. Loren Shure replied on :

Thanks, Yair.

I also discussed memoization quite a while ago (2006):

http://blogs.mathworks.com/loren/2006/02/08/use-nested-functions-to-memoize-costly-functions/

It can be a big win sometimes.

–Loren

 Name (required) E-mail (required, will not be published) Website (optional) Spam protection (required): What is 2 + 10 ?

Wrap code fragments inside <pre> tags, like this:

<pre class="code">
a = magic(3);
sum(a)
</pre>


If you have a "<" character in your code, either follow it with a space or replace it with "&lt;" (including the semicolon).

Loren Shure works on design of the MATLAB language at MathWorks. She writes here about once a week on MATLAB programming and related topics.

These postings are the author's and don't necessarily represent the opinions of MathWorks.