Loren on the Art of MATLAB

Fibonacci and filter 16

Posted by Loren Shure,

Over the years, people have posted many versions of code for calculating Fibonacci numbers. I thought I'd collect a few of the ideas here and start a discussion on the merits of the different versions.

Contents

Fibonacci Numbers

Here's a place to get some history about the man and his work. Over the centuries, people have that these numbers seem to come up in the natural world. Try searching for Fibonacci and sunflower; I just got 31,600 hits. Fibonacci numbers are a sequence of numbers that can be derived via a recurrence relationship. Start the sequence with two |1|s. Then, each subsequent number in the sequence is found by adding the two most recent numbers. Mathematically stated:

Straightforward for-loop

Let's calculate a reasonable chunk of the sequence through the straightforward method suggested by the formula above.

clear
nf = 102;
tic
fibf(1) = 1;
fibf(2) = 1;
for n = 3:nf
    fibf(n) = fibf(n-1)+fibf(n-2);
end
ftimes(1) = toc;  % Collect tic-toc times.

And now let's look at the first few numbers to be sure we're getting the expected answers.

fibf(1:10)
ans =
     1     1     2     3     5     8    13    21    34    55

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. We'll compare all the computation times later on in this post.

tic
fibp = ones(1,nf);
for n = 3:nf
    fibp(n) = fibp(n-1) + fibp(n-2);
end
ftimes(end+1) = toc;

Recursive Function

The mathematical formula above suggests that we could write a Fibonacci sequence algorithm using a recursive function, i.e., one that calls itself. Here's what I came up with. Note that this version grows an array each time. I might have been able to be clever about this. I doubt the code would be as clear, however. Feel free to post a version to prove me wrong!

type fibrec
function f = fibrec(n)
%FIBREC Recursive function for n Fibonacci numbers.

% Minimize the error checking so we don't bog down in it. I have included it
%in comments instead.
% if ~isscalar(n) | ~isreal(n) | n<0 | fix(n)~=n
%    error('ArtBlog:fibrec:MBRealPosInt','N must be a real positive integer')
% end
if n == 1,
    f = 1;   % First element is 1.
    return;
elseif n == 2
    f = [1 1];  % First two elements are 1.
else
    % Call fibrec with previous result and compute next one from it.
    fm1 = fibrec(n-1); 
    f = [fm1 fm1(end)+fm1(end-1)];
end

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

tic
fibr = fibrec(nf);
ftimes(end+1) = toc;

Use filter to Calculate the Sequence

I have one other idea now, to use the function filter. Looking at the help,

help filter
 FILTER One-dimensional digital filter.
    Y = FILTER(B,A,X) filters the data in vector X with the
    filter described by vectors A and B to create the filtered
    data Y.  The filter is a "Direct Form II Transposed"
    implementation of the standard difference equation:
 
    a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)
                          - a(2)*y(n-1) - ... - a(na+1)*y(n-na)
 
    If a(1) is not equal to 1, FILTER normalizes the filter
    coefficients by a(1). 
 
    FILTER always operates along the first non-singleton dimension,
    namely dimension 1 for column vectors and non-trivial matrices,
    and dimension 2 for row vectors.
 
    [Y,Zf] = FILTER(B,A,X,Zi) gives access to initial and final
    conditions, Zi and Zf, of the delays.  Zi is a vector of length
    MAX(LENGTH(A),LENGTH(B))-1, or an array with the leading dimension 
    of size MAX(LENGTH(A),LENGTH(B))-1 and with remaining dimensions 
    matching those of X.
 
    FILTER(B,A,X,[],DIM) or FILTER(B,A,X,Zi,DIM) operates along the
    dimension DIM.
 
    See also FILTER2 and, in the Signal Processing Toolbox, FILTFILT.

    Overloaded functions or methods (ones with the same name in other directories)
       help timeseries/filter.m


    Reference page in Help browser
       doc filter

I see that it could be helpful since prior output values within the filtered array can be used to calculate subsequent ones. I want the latest one to add together the previous two outputs. My input, x is the array [1 1 0 0 ... 0]. And for each output, I need only consider the current input, assuming I pad my original starting value of 1 with zeros. Therefore, for this case, b = 1. The leading coefficient, multiplying the last output, is 1. And then I want to add the two previous outputs, but in the equation, these appear with -a(k), so the a portion of the filter is [1 -1 -1].

tic
x = [1 zeros(1,nf-1)];
a = [1 -1 -1];
b = 1;
fibfil = filter(b, a, x);
ftimes(end+1) = toc;

Let's be sure I got the right answer.

isequal(fibf, fibp, fibr, fibfil)
ans =
     1

isequal Tip

Did you know isequal could take more than two inputs? If you store your items to compare in a cell array, C, you can then do this: passfail = isequal(C{:}):, taking advantage of the comma-separated list notation.)

Compare Computation Times

It's time to see how quickly the various methods performed. I 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    ', ftimes*1000), fprintf('\n')
Compare Times (in milliseconds)

loop           prealloc       recursive      filter         

0.329650836    0.243606380    2.244419333    0.058666674    

filter is Useful for Vectorization

I've shown one example here in which filter was useful for vectorizing code and make it perform well. Another one I know about is code from John D'Errico for calculating a moving window standard deviation. There are probably others. Do you know of other good uses for filter, besides straightforward applications for filtering and convolving signals? Let me know.


Published with MATLAB® 7.2

16 CommentsOldest to Newest

Hi,
Other examples of running/moving stats. was provided in CSSM by legendary Peter Acklam:

with the title “Running (moving) sum, mean, variance and standard deviation (Was: Vectorized version of moving or running standard deviation)” on Friday 12 Mar 2004 10:16.

However, for my own applications, where I may find occasional missing data (NaNs), FILTER return all NaNs and I use a for loop (Ouch :-)

/Lars

Here’s my non-array-growing version of fibrec. I only had to change the two lines in the ‘else’ clause.

Steve

function f = fibrec(n)
%FIBREC Recursive function for n Fibonacci numbers.

% Minimize the error checking so we don’t bog down in it. I have included it
%in comments instead.
% if ~isscalar(n) | ~isreal(n) | n

Well that didn’t work. I’ll try again but only with the two lines I changed in the ‘else’ clause:

f = [fibrec(n-1) 0];
f(n) = f(n-2) + f(n-1);

Steve

If you work with Fibonacci numbers (or something similar where recursion might come handy) and have enough memory it might be cunning to do something like this (speeding up things by a factor of much):

function f = fibrec(n)
%FIBREC Recursive function for n Fibonacci numbers.

% Minimize the error checking so we don’t bog down in it. I have included it
%in comments instead.
% if ~isscalar(n) | ~isreal(n) | nlength(F)

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

end

f = F(1:n);

Somethings obviously got lost…

Before ‘if n==1′ there should be two lines with:

persistent F

if n>length(F)

To make above work…

Your filter code tests fine, but the text seems inconsistent with the code. You write:
> My input, x is the array [1 1 0 0 ... 0]. And for each
> output, I need only consider the current input, assuming
> I pad my original starting values of [1 1] with zeros.

In the code, x has only a single 1 followed by zeros. The second one is unnecessary because the filter picks it up from y[1], no?

Michael,

You are correct and I have fixed the text in the blog to reflect that. Thanks for the correction!

–Loren

If you don’t want all n numbers, the following can give you a single one: Note, this gives f0 = 1, f = 1, f2 = 2, etc…
Also, at around n = 70 (near bitmax) this starts deviating from the results returned by addition based implementations

function fn = fibonnacci(n)
alpha1 = (1 + sqrt(5)) / 2;
alpha2 = (1 – sqrt(5)) / 2;

c1 = (1 – alpha2) / (alpha1 – alpha2);
c2 = (alpha1 – 1) / (alpha1 – alpha2);

fn = round(c1 * alpha1^n + c2 * alpha2^n);

% in fact the alpha2 term is really small, so you don’t
% need it, except for n = 0;
fn = round(c1 * alpha1^n);

On my PC, I get the following timings for nf=102 :
prealloc :0.015, filter:0.053ms.
For nf=10000, prealloc:3.12,filter:4.61 ms
That is, the prealloc is faster. I have checked this several times. Can this depend on memory, PC internals and so on?
I found another method described in Roman Maeder’s Mathematica book.
Let m=[1 1;1 0];
The matrix, m^(n-1) gives [f(n) f(n-1);f(n-1) f(n-2)]
where f(n) is fibonacci number for n, f(n-1) for (n-1) and so on.
fib=m^(n-1); fib(1,1) gives the fibonacci number for n. As expected,this computation is fast.
It would be interesting to see how the code can be vectorised with the above method to get the entire series upto n. Any volunteers?

My prealloc results might have been higher because of the way I ran the code. Did you run the preallocation in a script, function, or directly from the command line? That can make a huge difference sometimes, especially command-line vs. file.

What’s wrong with using the closed form expression? Even if you want to calculate the whole series of numbers it will be much faster. For example,

tic;
a=(1+sqrt(5))/2;
nf=102;
b=1-a;
F=(a.^[1:nf]-b.^[1:nf])/sqrt(5);
toc;

gives me report

Elapsed time is 0.016000 seconds

Eric

Correction to the previous post.

In fact, the 2nd time I run the code I get

Elapsed time is 0.000000 seconds

I believe the 1st time I ended up with 16ms because of the Hard Disk (or memory) read, while the 2nd time the program was running from cache.

Eric

hey loren-
i am interested to hear what you or steve think about the idea of matlab supporting gpgpu (using the stream processing hardware of modern video cards to perform vector/matrix operations up to 50 times faster than on the cpu) from matlab. obvious example related to your article: filtering/fft. wired has a great story on recent exciting results in this field:
http://www.wired.com/news/technology/computers/0,72090-0.html?tw=wn_index_16

(i’m not sending to steve because most of the image processing blog is over my head, so i don’t keep up with it — but it seems like that community would also be very interested.)
-erik

hi Loren,
I just wanted to kmow how matlab calls the recursive function when it is repeated twice as in fibonacci series???

POONKODI-

It just calls it twice. The first call completes (recursion and all) before the second one gets called. This version of recursion is quite wasteful since it doesn’t store and therefore use values that have been previously calculated.

–Loren

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