Deriving a Closed-Form Solution of the Fibonacci Sequence using the Z-Transform

The Fibonacci sequence might be one of the most famous sequences in the field of mathmatics and computer science. Already high school students starting with programming classes compute the first few Fibonacci numbers with their programs using different iterative or recursive approaches. One reason for its popularity might be that the Fibonacci sequence is closely related to many other fields of math and physics, often in very astonishing ways which one might not expect. Usually, the Fibonacci sequence is defined in a recursive manner. Hence, in order to compute the n-th Fibonacci number all previous Fibonacci numbers have to be computed first. In this blog post we will derive an interesting closed-form solution to directly compute any arbitrary Fibonacci number without the necessity to obtain its predecessors first. Interestingly, we will solve this problem with the help of a tool – the so called Z-Transform – which is actually more common in the field of digital signal processing.

The Fibonacci sequence, starting with

can be defined recursively as

A simple recursive function in R, implementing the sequence could look like this:

fibo <- function(n) {
  if(n<0) return (NaN);
  if(n == 0 ) return (0);
  if(n == 1) return (1);
  return (fibo(n-1) + fibo(n-2));
}

sapply(0:19, fibo)

The first 20 elements of the Fibonacci Sequence are thus:
0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181.

Surprisingly (not really, actually, if you think about it) the Fibonacci sequence can also be generated using an IIR (infinite impuls response) filter. Consider the difference equation of an IIR-filter in the form:

The impuls response of this filter is defined as:

where is the Kronecker Delta-Function. So at the time , we give a single impulse into our system and it starts running and computing a sequence of numbers. We get:

Let us now compute the impulse response for the given filter with some code and plot the results. Since the filter coefficients usually have to be passed in a strange way to most functions of the dsp-toolboxes, which requires reading off the coefficients from the transfer function , let us first compute the transfer function (the Z-transform of our filter in time domain):

Now, we can read off the filter coefficients (from the numerator for the forward coefficients and from the denominator for the reverse coefficients). Hence, the forward coefficients are

and the reverse coefficients are

Subsequently, we can compute and plot the impulse response of our system.

library(signal)
library(ggplot2)
b = c(1)
a = c(1,-1,-1)
h <- impz(filt=b,a=a, n=10)$x # compute impulse response for our Filter

# Plot the impulse response
df <- data.frame(x=h$t, y=h$x)
p <- ggplot(data = df, aes(x=x, y=y)) +
  geom_segment( aes(x=x, xend=x, y=0, yend=y), color="darkgrey", size=1) +
  geom_point(size=4, shape=21, fill="orange") +
  geom_text(aes(label=y),hjust=0, vjust=-.7) +
  theme_light(base_size = 18)+
  scale_x_continuous(breaks=c(0:10))+
  xlab("n")+
  ylab("h[n]")
p
alt and id

Figure 1: The impulse response (first 10 values) of our second-order IIR filter.

Not surprisingly, the impulse response corresponds exactly to the Fibonacci sequence.

In Eq. \eqref{eq:z-transform}, we already computed the transfer function of our “Fibonacci”-filter. Let us now see, if we can obtain another representation of the impulse response in time domain, which is no longer recursive and which depicts a closed-form description of the Fibonacci numbers. Such a representation would really be astonishing. We will do the following:

  1. Compute the partical fraction decomposition of our transfer function in Eq. \eqref{eq:z-transform}
  2. Look at the Z-transform of an certain type of infinite series
  3. Use the insights from 2. to transform our transfer function back into time domain and removing the recursive structure of the impulse response
  4. Take the new impulse response to generate arbitrary Fibonacci numbers immediately

1. Partial Fraction Decomposition of H(z)

Initially, we have to find the poles of our system hence, the roots in the denominator:

The 2 (non-complex) roots can be trivially found to be:

Now we can write our transfer function as:

In order to find the constants and , we can multiply with the denominator of each term and evaluate the expression at the poles we found previously (done examplarily here for the constant ):

Since the term with the constant vanishes due to the vanishing numerator, this leads to

Similarly, we obtain the constant :

This let’s us finally write our previous transfer function as:

2. Determing the inverse Z-Transform of a special Transfer Function

In order to transfer back our transfer function in Eq. \eqref{eq:partialfrac} back into the time domain, we have to find a suitable inverse transformation for an expression of the type

This is not a trivial task. However, if we remember the geometric series, we might be able to continue. We notice that the infinite sum

converges to

if .

If we look at Eq. \eqref{eq:ztransinfinite}, we can see that we have something which is very similar to Eq. \eqref{eq:geom}. By identifying

we can write

Note that this step assumes that or .

For this sum that we found, the inverse Z-transform can now be easily determined:

where is the unit-step function which turns on at This is necessary, since the sum in Eq. \eqref{eq:sumZ} starts with and is not defined (or is zero) for negative indexes.

We can summarize the findings of this section in basically one equation:

3. Inverse Transform of the previous Partial Fraction Decomposition

No we are ready, to transform Eq. \eqref{eq:partialfrac} back into the time domain by using the relation in Eq. \eqref{eq:finalrelation}. Let’s do it. Remember:

This leads us to:

which – when looking at it closer – is a truely remarkable result. There is no recursive formulation in the impulse response any longer. This means that we can compute any arbitrary Fibonacci number directly using this closed-form solution.

Another interesting observation from above equation is, that the term

is the so called golden ratio The golden ration itself has many interesting properties, which one should take a look at. In this case, we will just use it to simplify our expression above.

Using

and the additional relation

we can re-write Eq. \eqref{eq:backintime} in the following way:

which is – as I find – a really nice formulation of the Fibonacci sequence. (There are also some issues with above representation, but lets forget about them at the moment).

4. Compute arbitrary Fibonacci Numbers using the Closed-Form Solution

A simple R-function implementing this closed form solution could look like this:

fibo <- function(n) {
  phi = (1+sqrt(5))/2
  (phi^(n) - (-phi)^(-n))/(2*phi-1)
}

Try it out! For exampe, I get the following results in the following for the following cases:

> sapply(0:19, fibo)
[1] 0  1  1  2  3  5  8  13  21  34  55  89  144  233  377  610  987 1597 2584 4181
> fibo(57)
[1] 365435296162

Markus Thill

Markus Thill
I studied computer engineering (B.Sc.) and Automation & IT (M.Eng.). Generally, I am interested in machine learning (ML) approaches (in the broadest sense), but particularly in the fields of time series analysis, anomaly detection, Reinforcement Learning (e.g. for board games), Deep Learning (DL) and incremental (on-line) learning procedures.

Derivation of a Weighted Recursive Linear Least Squares Estimator

\( \let\vec\mathbf \def\myT{\mathsf{T}} \def\mydelta{\boldsymbol{\delta}} \def\matr#1{\mathbf #1}\)In this post we derive an increme...… Continue reading

Conway's Game of Life

Published on April 13, 2019