Computing Fibonacci

Linear Algebra Done Right

Computing the Fibonacci sequence is a common exercise in your typical ‘Intro to Coding’ course. The exercise involves an introduction of concepts in dynamic programming in memoization. While the resulting recursive program is efficient, you can actually compute Fibonacci numbers directly with linear algebra.

This exposition was inspired by Excerise 5.C(16) in the Third Eidtion of Sheldon Axler’s Linear Algebra Done Right textbook.

The Fibonacci sequence

The Fibonacci sequence $F_1,F_2,\dots$ is defined by

$$F_1=1,F_2=1,\text{and}{F}_n=F_{n-2}+F_{n-1}\text{ for }n\ge 3.$$

You may have seen the Fibonacci sequence in the context of the so-called “golden spiral”.

Computing Fibonacci numbers

Notice that to compute a Fibonnaci number $F_n$, we have to compute all Fibonnaci numbers $[3,n-1]$. To perform this computation, we can make a simple recursive program.

def fibonacci(n):
    if n <= 1:
        return n
    return fibonacci(n - 1) + fibonacci(n - 2)

Now let’s perform some tests.

for n in [2,5,10,25,40,41]:
    start_time = time.time()
    print(f"Fibonacci({n}) = {fibonacci(n)}")
    print(f"Execution time: {time.time() - start_time:.5f} seconds")

Fibonacci(2) = 1
Execution time: 0.00002 seconds
Fibonacci(5) = 5
Execution time: 0.00000 seconds
Fibonacci(10) = 55
Execution time: 0.00001 seconds
Fibonacci(25) = 75025
Execution time: 0.00705 seconds
Fibonacci(40) = 102334155
Execution time: 9.58250 seconds
Fibonacci(41) = 165580141
Execution time: 15.58907 seconds

As you can see, computing Fibonacci numbers $F_1,F_5,F_{10},$ and $F_{25}$ are relatively painless. But computing $F_{40}$ $F_{41}$ takes much longer, with the difference between the two execution times being $>5$ seconds. Why?

This is because our program is highly redundant. To compute $F_n$, we compute $F_{n-1}$ and $F_{n-2}$, but in computing $F_{n-1}$ we also compute $F_{n-2}$. Thus, we are computing $F_{n-2}$ twice. The problem only becomes worse as $n$ approaches $0$ in our program.

To remedy this redundancy, we employ dynamic programming. While an exposition of dynamic programming is outside the scope of this notebook, the idea is to create a dictionary, also called a memo, that stores the results of computing $F_j$ for $j<n$ and then to return those results when asked.

def fibonacci_memo(n, memo={}):
    if n in memo:
        return memo[n]
    if n <= 1:
        return n
    memo[n] = fibonacci_memo(n - 1, memo) + fibonacci_memo(n - 2, memo)
    return memo[n]

Now let’s rerun our tests from before.

for n in [2,5,10,25,40,41]:
    start_time = time.time()
    print(f"Fibonacci({n}) = {fibonacci_memo(n)}")
    print(f"Execution time: {time.time() - start_time:.5f} seconds")

Fibonacci(2) = 1
Execution time: 0.00002 seconds
Fibonacci(5) = 5
Execution time: 0.00000 seconds
Fibonacci(10) = 55
Execution time: 0.00000 seconds
Fibonacci(25) = 75025
Execution time: 0.00001 seconds
Fibonacci(40) = 102334155
Execution time: 0.00000 seconds
Fibonacci(41) = 165580141
Execution time: 0.00000 seconds

Cool right? That took barely any time, and all we did was add a dictionary.

But can we do better? Is there a closed form solution to computing Fibonnaci numbers?

Linear maps

Let’s define the linear map $T\in \mathcal{L}({\mathbf{R}}^2)$ by

$$T(x,y)=(y,x+y).$$

Now suppose we had to compute $T^n(0,1)$ for some positive integer $n$. For $n=1,2,3,4,5$, we have $T^1(0,1)=(1,1),T^2(0,1)=(1,2),T^3(0,1)=(2,3),T^4(0,1)=(3,5),T^5(0,1)=(5,8).$ Notice that the $1^{\text{st}}$ coordinate of $T^n(0,1)$ exactly matches $F_n$ and the second coordinate matches $F_{n+1}$. We can even see that $T^6(0,1)$ will produce $8=F_6$. So far, it looks like $T^n(0,1)=(F_n,F_{n+1})$, but can we prove it?

We can prove it with induction. We know that $T^1(0,1)=(F_1,F_2)$, so let our induction hypothesis be

$$T^{n-1}(0,1)=(F_{n-1},F_n).$$

Applying $T$ to both sides, we get

$$T^n(0,1)=T(F_{n-1},F_n)=(F_n,F_{n-1}+F_n)=(F_n,F_{n+1}).$$

Using matrices

So given that we know our linear map $T$ computes Fibonnaci numbers, could we simplify the process of iteratively applying $T$ $n$-times? Consider the following matrix $A\in {\mathbf{R}}^{2\times 2}$ defined by

$$A=\left[\begin{array}{cc}0& 1\\ 1& 1\end{array}\right].$$

It’s easy to see that for the vector $x^T=\left[\begin{array}{cc}0& 1\end{array}\right]$, we can write

$$A^{n}x={\left[\begin{array}{cc}0& 1\\ 1& 1\end{array}\right]}^n\left[\begin{array}{c}0\\ 1\end{array}\right]=\left[\begin{array}{c}{F}_n\\ F_{n+1}\end{array}\right].$$

Wouldn’t it be nice if our matrix $A$ were diagonalizable so that

$$A=Q\Lambda Q^{-1}\text{and}{A}^n=Q\Lambda Q^{-1}Q\Lambda Q^{-1}\cdots Q\Lambda Q^{-1}=Q{\Lambda }^n{Q}^{-1}$$

since $Q{Q}^{-1}=I$. Given $\Lambda$ is a diagonal matrix, ${\Lambda }^n$ is simply obtained by raising each diagonal entry to the $n$-th power. This would make the computing of Fibonnaci numbers near trivial. Lucky for us, $A$ is indeed diagonalizable. Let’s start by findind the eigenvalues of $T$.

Eigenvalues

To compute the eigenvalues of $T$, we can write

$$T(x,y)=(\lambda x,\lambda y)\text{and}(y,x+y)=(\lambda x,\lambda y).$$

Combining $y=\lambda x$ and $x+y=\lambda y$, we have

$$x+\lambda x={\lambda }^{2}x$$

which becomes

$$0={\lambda }^2-\lambda -1$$

leaving us to find the roots of the polynomial $f(\lambda )={\lambda }^2-\lambda -1$. Via the tried and true quadratic formula, it follows that

$$\lambda =\frac{1\pm \sqrt{5}}{2}.$$

Hence, the eigenvalues of $T$ are ${\lambda }_1=\frac{1+\sqrt{5}}{2}$ and ${\lambda }_2=\frac{1-\sqrt{5}}{2}$. In particular, our eigenvalue ${\lambda }_1$ is the so-called “golden ratio”.

Eigenvectors

Now let’s find the eigenvectors. For ${\lambda }_1$, we can set $x=1$ and write $T(1,y)=({\lambda }_1,{\lambda }_{1}y)\text{and}(y,1+y)=({\lambda }_1,{\lambda }_{1}y)$ suggesting that $y={\lambda }_1$ and $1+{\lambda }_1=({\lambda }_1{)}^2$. While initially unbelievable, it is indeed the case that

$$1+\frac{1+\sqrt{5}}{2}=(\frac{1+\sqrt{5}}{2}{)}^2,$$ $$\frac{3+\sqrt{5}}{2}=\frac{(1+\sqrt{5}{)}^2}{4},$$ $$\frac{3+\sqrt{5}}{2}=\frac{6+2\sqrt{5}}{4},$$ $$\frac{3+\sqrt{5}}{2}=\frac{3+\sqrt{5}}{2}.$$

Notice that the same result also holds for $1+{\lambda }_2=({\lambda }_2{)}^2$. Hence, the eigenvectors of $T$ are $q_1=(1,{\lambda }_1)\text{and}{q}_2=(1,{\lambda }_2).$

Putting linear algebra to work

At this point, we now have the following algorithm for computing Fibonacci numbers

$$Q{\Lambda }^n{Q}^{-1}x=\frac{1}{{\lambda }_2-{\lambda }_1}\left[\begin{array}{cc}1& 1\\ {\lambda }_1& {\lambda }_2\end{array}\right]\left[\begin{array}{cc}({\lambda }_1{)}^n& 0\\ 0& ({\lambda }_2{)}^n\end{array}\right]\left[\begin{array}{cc}{\lambda }_2& -1\\ -{\lambda }_1& 1\end{array}\right]\left[\begin{array}{c}0\\ 1\end{array}\right]=\left[\begin{array}{c}{F}_n\\ F_{n+1}\end{array}\right]$$

where ${\lambda }_2-{\lambda }_1=\det Q$. Let’s sanity check this algorithm by writing a python function.

def fibonacci_matrix(n):
    sqrt5 = np.sqrt(5)
    lambda1 = (1 + sqrt5) / 2
    lambda2 = (1 - sqrt5) / 2
    
    Q = np.array([[1, 1], [lambda1, lambda2]])
    Lambda_n = np.diag([lambda1**n, lambda2**n])
    Q_inv = np.array([[lambda2, -1], [-lambda1, 1]]) / (lambda2 - lambda1)
    x = np.array([[0], [1]])
    
    result = Q @ Lambda_n @ Q_inv @ x
    return result[0, 0].item()

Like before, let’s rerun those tests.

for n in [2,5,10,25,40,41]:
    start_time = time.time()
    print(f"Fibonacci({n}) = {fibonacci_matrix(n)}")
    print(f"Execution time: {time.time() - start_time:.5f} seconds")

Fibonacci(2) = 0.9999999999999999
Execution time: 0.00020 seconds
Fibonacci(5) = 5.000000000000001
Execution time: 0.00003 seconds
Fibonacci(10) = 55.00000000000002
Execution time: 0.00002 seconds
Fibonacci(25) = 75025.00000000006
Execution time: 0.00001 seconds
Fibonacci(40) = 102334155.00000013
Execution time: 0.00001 seconds
Fibonacci(41) = 165580141.00000024
Execution time: 0.00001 seconds

While fibonacci_matrix is much faster our original fibonacci program, it’s still beat by fibonacci_memo in speed. Also notice that the numbers are not integers as there are small remainders. These remainders occur due to floating-point precision errors in computing with $\sqrt{5}$ in our eigenvalues.

With these issues in mind, we can do better.

From our expression for $Q{\Lambda }^n{Q}^{-1}x$ from before, we can write

$$\frac{1}{{\lambda }_2-{\lambda }_1}\left[\begin{array}{cc}1& 1\\ {\lambda }_1& {\lambda }_2\end{array}\right]\left[\begin{array}{cc}({\lambda }_1{)}^n& 0\\ 0& ({\lambda }_2{)}^n\end{array}\right]\left[\begin{array}{cc}{\lambda }_2& -1\\ -{\lambda }_1& 1\end{array}\right]\left[\begin{array}{c}0\\ 1\end{array}\right]$$ $$\frac{1}{{\lambda }_2-{\lambda }_1}\left[\begin{array}{cc}1& 1\\ {\lambda }_1& {\lambda }_2\end{array}\right]\left[\begin{array}{cc}({\lambda }_1{)}^n& 0\\ 0& ({\lambda }_2{)}^n\end{array}\right]\left[\begin{array}{c}-1\\ 1\end{array}\right]$$ $$\frac{1}{{\lambda }_2-{\lambda }_1}\left[\begin{array}{cc}1& 1\\ {\lambda }_1& {\lambda }_2\end{array}\right]\left[\begin{array}{c}-({\lambda }_1{)}^n\\ ({\lambda }_2{)}^n\end{array}\right]$$ $$\frac{1}{{\lambda }_2-{\lambda }_1}\left[\begin{array}{c}-({\lambda }_1{)}^n+({\lambda }_2{)}^n\\ -({\lambda }_1{)}^{n+1}+({\lambda }_2{)}^{n+1}\end{array}\right]$$

implying that

$$F_n=\frac{-({\lambda }_1{)}^n+({\lambda }_2{)}^n}{{\lambda }_2-{\lambda }_1}.$$

Substituting in the exact values of ${\lambda }_1$ and ${\lambda }_2$, we have

$$F_n=\frac{2}{(1-\sqrt{5})-(1+\sqrt{5})}[-(\frac{1+\sqrt{5}}{2}{)}^n+\left(\frac{1-\sqrt{5}}{2}{)}^n\right]$$

which is equivalently expressed as

$$F_n=\frac{1}{\sqrt{5}}\left[\right(\frac{1+\sqrt{5}}{2}{)}^n-\left(\frac{1-\sqrt{5}}{2}{)}^n\right].$$

Doesn’t it look much nicer? Let’s make a new program.

def fibonacci_expression(n):
    sqrt5 = np.sqrt(5)
    lambda1 = (1 + sqrt5) / 2
    lambda2 = (1 - sqrt5) / 2
    
    result = (lambda1**n - lambda2**n)/sqrt5
    return result

You know the drill.

for n in [2,5,10,25,40,41]:
    start_time = time.time()
    print(f"Fibonacci({n}) = {fibonacci_expression(n)}")
    print(f"Execution time: {time.time() - start_time:.5f} seconds")

Fibonacci(2) = 1.0
Execution time: 0.00020 seconds
Fibonacci(5) = 5.000000000000001
Execution time: 0.00001 seconds
Fibonacci(10) = 55.000000000000014
Execution time: 0.00001 seconds
Fibonacci(25) = 75025.00000000006
Execution time: 0.00000 seconds
Fibonacci(40) = 102334155.00000013
Execution time: 0.00000 seconds
Fibonacci(41) = 165580141.00000024
Execution time: 0.00000 seconds

As you can see, this new program is much quicker, but we still have those pesky floating-point precision errors. The obvious solution would be to simply round those errors off, but if we are going to resort to rounding, we can actually make an even better program.

Approximating Fibonacci

Let’s take another look at our expression for Fibonacci numbers

$$F_n=\frac{1}{\sqrt{5}}\left[\right(\frac{1+\sqrt{5}}{2}{)}^n-\left(\frac{1-\sqrt{5}}{2}{)}^n\right]=\frac{1}{\sqrt{5}}(\frac{1+\sqrt{5}}{2}{)}^n-\frac{1}{\sqrt{5}}(\frac{1-\sqrt{5}}{2}{)}^n.$$

Do you notice something interesting about the term $\frac{1}{\sqrt{5}}(\frac{1-\sqrt{5}}{2}{)}^n$? Computing the term for $n=1,2,3,4$, we have

$$\frac{1}{\sqrt{5}}(\frac{1-\sqrt{5}}{2}{)}^1=-0.2763\dots ,\frac{1}{\sqrt{5}}(\frac{1-\sqrt{5}}{2}{)}^2=0.1708\dots ,\frac{1}{\sqrt{5}}(\frac{1-\sqrt{5}}{2}{)}^3=-0.1055\dots ,\frac{1}{\sqrt{5}}(\frac{1-\sqrt{5}}{2}{)}^4=0.0652\dots$$

The absolute value of each of these terms is less than $\frac{1}{2}$. So if we are going to round our expression of the Fibonacci number, then we don’t even need to compute the term $\frac{1}{\sqrt{5}}(\frac{1-\sqrt{5}}{2}{)}^n$ since rounding will always make up for what we lost. Hence, we can write

$$F_n=\frac{1}{\sqrt{5}}(\frac{1+\sqrt{5}}{2}{)}^n.$$

Are you skeptical? Let’s make a fifth and final program to compute Fibonacci numbers.

def fibonacci_approximation(n):
    sqrt5 = np.sqrt(5)
    lambda1 = (1 + sqrt5) / 2    
    result = (lambda1**n)/sqrt5
    return round(result)

for n in [2,5,10,25,40,41]:
    start_time = time.time()
    print(f"Fibonacci({n}) = {fibonacci_approximation(n)}")
    print(f"Execution time: {time.time() - start_time:.5f} seconds")

Fibonacci(2) = 1
Execution time: 0.00005 seconds
Fibonacci(5) = 5
Execution time: 0.00001 seconds
Fibonacci(10) = 55
Execution time: 0.00001 seconds
Fibonacci(25) = 75025
Execution time: 0.00001 seconds
Fibonacci(40) = 102334155
Execution time: 0.00000 seconds
Fibonacci(41) = 165580141
Execution time: 0.00000 seconds

Isn’t that just wonderful? We now have a simple closed-form approximation of Fibonacci numbers with an execution time that rivals our memoization method.

The power of linearity

In this notebook, we explored how computing the Fibonacci sequence can be performed with brute-force recursion, memoizaiton, linear maps, matrices, eigenvalues, and rounding. I bet you won’t be surprised, but there are even more methods one can implement to go faster.

From our work, I believe the most striking takeaway is how powerful linearity is. While nonlinear methods have seen widespread adoption in machine learning, the Fibonacci sequence shows us that a great deal can be achieved by finding elegant linear approximations to complex mathematical objects. For more examples of the power of linearity and linear algebra, checkout my solutions to exercises in Sheldon Axler’s Linear Algebra Done Right textbook.