Math

Fibonacci: Linear Recurrences and Eigendecomposition in Python

Learn Linear Algebra and its Applications from Scratch

Fibonacci numbers form a spiral. Source: http://mevzuforex.com/wp-content/uploads/2019/05/Fibonacci-800x445.png

Introduction

Today we are going to explore the infamous Fibonacci sequence and use it as an example to explain linear recurrences and eigendecomposition. You should be familiar with what a vector and a matrix is and how we can do matrix multiplication.

The Fibonacci sequence is defined recursively as an = a(n-1) + a(n-2)

We start with a0 = 1 and a1 = 1
a2 = a1 + a0 = 1 + 1 = 2
a3 = a2 + a1 = 2+ 1 = 3 and so on

For every element we just take the sum of the previous two elements.

The Elements up to a13 of the Fibonacci Series computed

Source: https://i.pinimg.com/originals/98/82/d5/9882d569f7e0b5665fe3b2edd5069b06.png

Let us compute the Fibonacci Sequence step by step in Python first

# We start with the initial values a0=1 and a1=1
a0 = 1
a1 = 1
a2 = a1 + a0
a3 = a2 + a1
a4 = a3 + a2
for n,a in enumerate([a0,a1,a2,a3,a4]):
    print("a%s = %i" % (n,a))

Produces the following output:

a0 = 1
a1 = 1
a2 = 2
a3 = 3
a4 = 5

Switching to compute with nested function

1. We define an easy function that just adds two numbers

def fibonacci_simple(an_minus1, an_minus2):
    # just adds two last numbers in sequence a(n-1), a(n-2)
    an = an_minus1 + an_minus2
    return an

2. We compute a4 as an example

a0 = 1
a1 = 1
a2 = fibonacci_simple(a1,a0)
a3 = fibonacci_simple(a2,a1)
a4 = fibonacci_simple(a3,a2)
print("a%s = %i" % (4,a4))

Produces the following output:

a4 = 5

Nested it looks something like this:

This is what recursion does and it is not very efficient as we have four (4) function calls instead of the three (3) we had before and this becomes more and more complicated for larger n in the Fibonacci sequence.

a4 = fibonacci_simple(a3,a2)
a4 = fibonacci_simple(
fibonacci_simple(a2,a1),
fibonacci_simple(a1,a0))
a4 = fibonacci_simple(
fibonacci_simple(
fibonacci_simple(a1,a0),a1),
fibonacci_simple(a1,a0))
print("a%s = %i" % (4,a4))

Produces the same output:

a4 = 5

Computing with a Recursively Defined Function

A visual example of recursion. Source: https://www.smbc-comics.com/comics/1562409923-20190706.png

def fibonacci_recurssion(n):
 if n == 0 or n == 1:
 # after some recurssion steps n is finally 0 or 1 and we know that a0 = 1 and a1 = 1
 return 1
 # if the function does not know what it is, just let it figure out what a(n-1) or a(n-2) is, if it does not know
 # that then we can just try just a(n-3) or a(n-4) until n-x = 0 or n-x = 1, where we know a0 = a1 = 1
 return fibonacci_recurssion(n-1) + fibonacci_recurssion(n-2)
for n in range(0,10):
 print(“a%s = %i” % (n,fibonacci_recurssion(n)))

Produces the following output:

a0 = 1
a1 = 1
a2 = 2
a3 = 3
a4 = 5
a5 = 8
a6 = 13
a7 = 21
a8 = 34
a9 = 55

Computing with a Matrix

Let us try to be even smarter, the Fibonacci sequence seems to be a function that is linear. Thus we can define it as a matrix transformation.

import numpy as np
a0 = 1
a1 = 1
# c1 and c2 define how often we take the last value and the value before that
# The fibonacci sequence is defined as: an = a(n-1) + a(n-2)
# Which is the same as: an = 1*a(n-1) + 1*a(n-2)
# Or: an = c1*a(n-1) + c2*a(n-2) for c1 = 1 and c2 = 1
c1 = 1
c2 = 1
T = np.array([[c1,c2],
 [1,0]])
# the initial values a0 = 1 and a1 = 1
v0 = np.array([a1,a0])
print(“v0 = [a1 a0] = “,v0,”\n”)
# now we take the matrix vector product: [c1*a(n-1) + c2*a(n-2), 1*a1 + 0*a0] = [a2, a1] 
v1 = T @ v0
print(“v1 = [a2 a1] = “,v1)
# now we can apply the transformation as many times as we want
v2 = T @ v1
print(“v2 = [a3 a2] = “,v2)
v3 = T @ v2
print(“v3 = [a4 a3] = “,v3,”\n”)
# We rewrite this which gives us the same result obviously
v3 = T @ T @ T @ v0
print(“v3 = [a4 a3] = “,v3)

Produces the following output:

v0 = [a1 a0] =  [1 1] 

v1 = [a2 a1] =  [2 1]
v2 = [a3 a2] =  [3 2]
v3 = [a4 a3] =  [5 3] 

v3 = [a4 a3] =  [5 3]

Applying the matrix transformation multiple times

Now we got a nice way to compute the Fibonacci numbers with a matrix T, however, we still have 3 calculations for getting to a4, 9 for getting to a10 and so on.

print("T =\n",T,"\n")
print("T^2 =\n",T@T,"\n")
print("T^3 =\n",T@T@T,"\n")
print("A^4 =\n",T@T@T@T,"\n")

Produces the following output:

T =
 [[1 1]
 [1 0]] 

T^2 =
 [[2 1]
 [1 1]] 

T^3 =
 [[3 2]
 [2 1]] 

A^4 =
 [[5 3]
 [3 2]]

It is hard to find any pattern here that could help us speed up the computation of A⁴.

Enter the Diagonal Matrix

Let us consider a diagonal matrix A which has non-zero values only on its diagonal.

A = np.array([[2,0],
              [0,3]])
print("A:\n",A,"\n")
print("A^2:\n",A@A,"\n")
print("A^3:\n",A@A@A,"\n")
print("A^4:\n",A@A@A@A,"\n")

Produces the following output:

A:
 [[2 0]
 [0 3]] 

A^2:
 [[4 0]
 [0 9]] 

A^3:
 [[ 8  0]
 [ 0 27]] 

A^4:
 [[16  0]
 [ 0 81]]

Well, that is interesting, now the numbers on the diagonal just get squared, cubed and so on.

So we can rewrite AAAA = A⁴ as taking every element in A to the power of 4.
Remark: A⁴ means multiplying the matrix 4 times in math notation.
A**4 in Python means doing taking every element in A to the power of 4. (Element-wise operation!)

This returns the exact same result and is only a single computation step

# This return the exact same result and is only a single computation step
print("A**4:\n",A**4,"\n")

Produces the following output:

A**4:
 [[16  0]
 [ 0 81]]

This is an astounding observation. However not yet very helpful, as our matrix T is not a diagonal matrix.

Hunting down the Diagonal Matrix

Basis of a Matrix

Every n*n square matrix has n columns. These n columns represent coordinates. These coordinates show us how far we have to got to the

A two dimensional example:

The two yellow vectors form our standard basis:
[[1, 0]
[0, 1]] = [i, j] = I

We call this the identity matrix I.

The blue vectors have some other coordinates:
[[a, b]
[c, d]] = [i-dot, j-dot] = B

Source: https://images.slideplayer.com/15/4546640/slides/slide_2.jpg

We say the blue matrix is defined in the standard-basis because. This means the vectors in B scale the standard basis I. But what does this mean?

Consider i-dot = [a, c]. a tells us how many steps we have to go in the x direction and c how many steps we have to go in the c direction.

a [1, 0] + c [0, 1] = a * i+ c * j= [a, c]

As we see our steps are defined in the direction of i and j. We go a steps in the direction of i and j steps in the direction of j.

Change of Basis

Would it not be easier to say how many steps we want to go into the direction of i-dot and in the direction of j-dot? This means we would use B=[i-dot, j-dot] as our basis. We know we just go 1 step in the direction of i-dot.

Thus [a, c] = 1* i-dot+ 0* j-dot = [1, 0] in the basis of B.
and [a, c] = 0* i-dot+ 1* j-dot = [0, 1] in the basis of B.

So B in the basis of B = I =
[[1, 0]
[0, 1]]

This is a diagonal matrix on which we can multiply super fast with itself by element-wise exponentiation. After we have done the complicated operation we could just transform back to our standard basis I where B is now transformed to B^n for some n times we multiplied B with itself.

However not so fast, this actually does not work yet:

B² ≠ [I²]B = B… I in the basis of B

But we have something that can work.

Similar Matrices

A similar matrix S represents the same linear transformation in a different basis.

So T = [S]b … S in same basis b

Maybe this S could be a diagonal matrix.

Let S = P^(-1)*A*P … where P^(-1) is the inverse of P such that P^(-1)*P = I.

What happens if we square both sides?

S² = P^(-1)*A*P*P^(-1)*A*P = P^(-1)*A*A*P = P^(-1)*A²*P

If we want to know what A² is:

A² = P*S²*P^(-1)

This is the prove that if we find some P and its inverse for which
S = P^(-1)*A*P we can do the computation of A² (and any other power) in another basis.

But is any of those S a diagonal matrix?

The vectors that just get stretched by T: Eigenvalues and Eigenvectors

An eigenvector v is a vector that only gets scaled(stretched) by a linear transformation by the factor λ which we call an eigenvalue.

Tv = λv

Remember that we want to find some similar matrix S that is a diagonal matrix for which:

S = D = P^(-1)*T*P

If we find some eigenvectors in P that just get stretched by A instead of a normal transformation something interesting happens.

If P = [a b c], then T*P = [ λ1*a λ2*b λ3*c].

And by transforming TP into the basis of P.

D = P^(-1)*T*P =

[[ λ1 0 0]
[ 0 λ2 0]
[ 0 0 λ3]]

And now we can potentiate D easily by squaring each element on the diagonal n times. Then we just transform back to A and we are done.

T^n =P*D^n* P^(-1)

So we just have to perform the following steps

1. Find all the (linearly independent) eigenvectors (i.e. P)
2. Find the inverse of P -> P^(-1)
3. Compute D -> D = P^(-1)*T*P
4. Potentiate D -> D^n
5. Transform D^n into the similar matrix D^n -> T^n =P*D^n* P^(-1)

Steps 1 to 3 have to be done only once. So we only have to do steps 4 and 5 for every power n. We can solve Fibonacci rapidly now.

We did it, the magic trick was successful! Expelliarmus

Source: https://qph.fs.quoracdn.net/main-qimg-ee385f708c4e656081d239342511a9a7

Now let’s do it in Code: Diagonalization with numpy

1. Find all the (linearly independent) eigenvectors

eigenvalues, eigenvectors = np.linalg.eig(T)
print("eigenvalues:",eigenvalues, "\n")
print("eigenvectors:\n",eigenvectors)

Produces the following output:

eigenvalues: [ 1.61803399 -0.61803399] 

eigenvectors:
 [[ 0.85065081 -0.52573111]
 [ 0.52573111  0.85065081]]

So we know there are two eigenvalues which scale those two eigenvectors if we multiply T with them.

The first eigenvector v1 get scaled by λ1 like this:

T*v1 = λ1*v1
T*[ 0.85065081 -0.52573111]=
1.61803399*[ 0.85065081 -0.52573111]

Tip: You can find the linearly independent eigenvectors by first finding all real and complex roots (eigenvalues) of the characteristic polynomial det((A-λ)v) = 0 of T through co-factor expansion first and then solving det((A-λ)v) = 0 for each eigenvalue.

As there are a maximum of n independent eigenvectors in T P is just all the eigenvalues as columns. We find the inverse of P where P * P^(-1) = 1 with numpy.

2. Find the inverse of P

P = eigenvectors
P_inv = np.linalg.inv(P)# equals P^(-1)

Tip: You can find the inverse of P by row reduction of [P I] to [I P^(-1)].

3. Compute the Diagonal D

D = P_inv @ T @ P
print("Diagonal = \n",D)

Produces the following output:

Diagonal = 
 [[ 1.61803399e+00  1.11022302e-16]
 [ 2.22044605e-16 -6.18033989e-01]]

Why are the non-diagonal elements not zero? We have some floating point inaccuracy, so let us just extract the diagonal entries we are interested in.

diag_entries = np.diag(D)
print("Diagonal entries: \n",diag_entries)

Produces the following output:

Diagonal entries: 
 [ 1.61803399 -0.61803399]

4. Potentiate D

Let us find some eigenvalues now. We first try to find a9 in the Fibonacci sequence which is the first entry in A⁸[a1,a0].

D_8 = np.diag(diag_entries**8)
print("Diagonal to the power of 8 :\n",D_8)

Produces the following output:

Diagonal to the power of 8 :
 [[4.69787138e+01 0.00000000e+00]
 [0.00000000e+00 2.12862363e-02]]

5. Transform D⁸ into the similar matrix A⁸

T_8 = P @ D_8 @ P_inv
print("T to the power of 8 : \n",T_8)

Produces the following output:

T to the power of 8 : 
 [[34. 21.]
 [21. 13.]]

Finally!

Source: https://cdn.eventplanner.net/imgs/xnr8784_how-to-build-excitement-for-your-attendees.jpg

v8 = T_8 @ v0
print("v8 containing a9 and a8 : \n",v8)

We computed [a9 a8]:

v8 containing a9 and a8 : 
 [55. 34.]

And this is exactly the same values we got from recursion. We are done! :D

All the Code in one Box

You can change c1, c2 or the initial values a0 and a1 and it will still work.

import numpy as np
def find_D(c1=1,c2=1):
    
    # Define the initial value vector
    v0 = np.array([a1,a0])
# Define the recurrence relation
    T = np.array([[c1,c2],
                  [1,0]])
# Find the linearly independet eigenvectors, P and P inverse
    eigenvalues, eigenvectors = np.linalg.eig(T)
P = eigenvectors
    P_inv = np.linalg.inv(P)
# Find the diagonal
    D = P_inv @ T @ P
    diag_entries = np.diag(D)
    D = np.diag(diag_entries)
return D
# Find the Fibonacci sequence elment you are looking for
def find_Fibonacci(D,n,a0=1,a1=1):
    
    # no further computation necessary
    if n == 0: return a0
    if n == 1: return a1
# Define the initial value vector
    v0 = np.array([a1,a0])
    
    # Take the power of the Diagonal (i.e. do the transformation multiple times)
    D = D**(n-1)
# Find the transformation matrix in the standard basis
    T = P @ D @ P_inv
    
    # Compute the (n-1)th vector containing [Fn Fn-1]
    v = T @ v0
    Fn = int(round(v[0])) # rounding because we expect an integer, not a floating point number with floating point inaccuracy
    return Fn
D = find_D()
for n in range(0,10):
    print("a%s = %i" % (n,find_Fibonacci(D,n)))

And we get exactly what we expected, just a little faster!

a0 = 1
a1 = 1
a2 = 2
a3 = 3
a4 = 5
a5 = 8
a6 = 13
a7 = 21
a8 = 34
a9 = 55

Complexity

Edit: Someone suggested that his iterative approach is faster than my approach. However that is not true and we can prove it numerically. (Mathematicians would think that is joke) If we want to lookup the 10th number of the Fibonacci sequence it is just as fast as looking up the 1000th number practically. This means we have constant lookup time 0(1) that does not depend on n.
(O(log(n) actually because exponentiation takes log(n) time, but this is a minor difference.)

In the code we compare calculating the Fibonacci numbers 100k times for either small n ranging from index 0 to 100 or big n ranging from 0 to 1000.

t1 = time.time()
for i in range(1000):
    for n in range(100):
        a = find_Fibonacci(D,n)
time_for_small_n = time.time()-t1
t1 = time.time()
for i in range(100):
    for n in range(1000):
        a = find_Fibonacci(D,n)
time_for_big_n = time.time()-t1
print("time_for_small_n = %fs ~ time_for_big_n = %fs" % (time_for_small_n, time_for_big_n))

Produces the following output:

time_for_small_n = 3.126421s ~ time_for_big_n = 2.989327s

Takes roughly the same time.

Other approaches

Recursive approach — Big O(2^n)

This scales terrible and it already takes us 117s to calculate the first 40 Fibonacci numbers.
(We calculated the first 1000 100 times in 3 seconds before!)

t1 = time.time()
for n in range(1,40):
    a = fibonacci_recurssion(n)
print("Time: %fs" % time.time()-t1)
Output: Time: 117.35020399093628s

Also if n becomes bigger it becomes impossible to calculate it that way because the recursive approach scales with 2^n. We only had a few billion years to calculate Fibonacci yet.

Iterative approach — Big O(n)

Just calculating a0, then a1, then a2 until we are at an is much easier and more scaleable. However if we want to lookup the nth Fibonacci number it takes n steps. We had constant lookup times.

We define a function for this:

def fibonacci_iterative(n):
    if n == 0:
        return 1
    
    computation_steps = n - 1
    
    Fn_previous = 1
    Fn_before_previous = 1
    
    for _ in range(computation_steps):
        temp = Fn_previous
        Fn_previous = Fn_previous + Fn_before_previous
        Fn_before_previous = temp
    
    Fn_current = Fn_previous
    return Fn_current

And test it on the same task as before, finding Fibonacci numbers 100k times for the first 100 and 1000 Fibonacci numbers. This is actually a super fast approach for small n’s and also feasible for finding big n. It just takes 10 times longer to find the 1000th then the 100th Fibonacci number.

t1 = time.time()
for i in range(1000):
    for n in range(100):
        a = fibonacci_iterative(n)
time_for_small_n = time.time()-t1
t1 = time.time()
for i in range(100):
    for n in range(1000):
        a = fibonacci_iterative(n)
time_for_big_n = time.time()-t1
print("time_for_small_n = %fs ~ time_for_big_n = %fs" % (time_for_small_n, time_for_big_n))

Produces the following output:

time_for_small_n = 0.798730s ~ time_for_big_n = 6.187797s

As we see it is really fast for small n and even faster than our sophisticated code. However, for big n it takes longer.

Reducing some of the unnecessary flexibility of our approach like setting initial values and computing T for a given

Speeding up O(log(n)) even more

We know that the nth Fibonacci number is a linear combination of the eigenvalues to the power n such that Fn = a*λ1^n + b*λ2^n for some a and some b.

We know this to be true because this is the same thing our diagonal does, first taking the power of the eigenvalues and then linearly transforming them back into the standard basis.

We know that F0= a*(λ1)⁰ + b*(λ2)⁰ = a + b = 1 and F1= a*λ1 + b*λ2 = 1 such that:

a*λ1 + b*λ2 = a + b = 1

If we can find a and b we just have to compute Fn = a*λ1^n + b*λ2^n to find the nth Fibonacci number Fn and that makes it even easier. We already know what the two eigenvalues λ1 and λ2 are.

We use numpy to solve the linear system

Source: Pen and Paper

This finds the solution [a b] for the linear system on the left

eigenvalue1, eigenvalue2 = eigenvalues
linear_system = np.array([[eigenvalue1,eigenvalue2],
                          [1,1]])
a,b = np.linalg.solve(np.array([eigenvalues,[1,1]]),[1,1])

Then we define a function that computes Fn

def find_Fibonacci_fast(n,a=a,b=b,eigenvalue1=eigenvalue1,eigenvalue2=eigenvalue2):
    Fn = a*eigenvalue1**n + b*eigenvalue2**n
    Fn = int(round(Fn))
    return Fn

Which works very well

for n in range(0,10):
    print("a%s = %i" % (n,find_Fibonacci_fast(n)))

Produces the following output:

a0 = 1
a1 = 1
a2 = 2
a3 = 3
a4 = 5
a5 = 8
a6 = 13
a7 = 21
a8 = 34
a9 = 55

And performs not only with close to constant lookup times 0(log(n)), but extremely fast just like we expected (3 times faster than the matrix operations)

t1 = time.time()
for i in range(1000):
 for n in range(2,102):
 a = find_Fibonacci_fast(n)
time_for_small_n = time.time()-t1
t1 = time.time()
for i in range(100):
 for n in range(2,1002):
 a = find_Fibonacci_fast(n)
time_for_big_n = time.time()-t1
print(“time_for_small_n = %fs ~ time_for_big_n = %fs” % (time_for_small_n, time_for_big_n))

Produces the following output:

time_for_small_n = 0.995449s ~ time_for_big_n = 0.967372s

If you made it that far and understand most of it congratulations!
You should now have a solid computational basic grasp of linear algebra.

Wow!

Read some of my other (way shorter 😜) articles:

Ways to Change a Dollar: Counting in Python

Why Mathematicians hate Statistics and Machine Learning

The Shortest Introduction To Deep Learning You Will Find On The Web

Build a Scalable Search Engine for Fuzzy String Matching

Don’t forget to follow me to not miss out on any new posts on AI, machine learning, math and entrepreneurship! Christoph Ostertag

Christoph Ostertag - Medium