In this post we solve the Fibonacci sequence using linear algebra. The Fibonacci equation is a second-order difference equation which is a particular type of sequence.

## Sequences

A sequence is a (possibly infinite) set of numbers with a defined order.

$a_n = \frac{1}{n}, \textit{ for } n = 0, 1, 2, ...,$

### Difference Equations

A difference equation is a type of sequence. A difference equation expresses the value of the sequence relative to the value of another point in the sequence.

$x_{n+1} = 1.5x_n, \textit{ for } n = 0, 1, 2, ...,$

To evaluate this equation at $n=3$ we need to know the initial value $x_0$, and evaluate $n=1, 2, 3$ in succession. For example, if $x_0=100$:

\begin{align} x_0 & && = 100 \\ x_1 & = 1.5 \times 100 && = 150\\ x_2 & = 1.5 \times 150 && = 225\\ x_3 & = 1.5 \times 225 && = 337.5 \end{align}

In the example above, we evaluate the equation recursively. It would be better to be able to evaluate the equation at $n = 100$ without performing $100$ calculations. Note that we can also write the equation above in this form:

$x_3 = 100 \times 1.5 \times 1.5 \times 1.5 = 337.5$

Solving a difference equation means writing the same equation such that it can be evaluated at a given $n$ without first evaluating $n=1,2,etc$. The solution to the equation above is:

$x_{n} = 1.5^{n} \times x_0, \textit{ for } n = 0, 1, 2, ...,$

### Second Order Difference Equations

The previous example is a first order difference equation because it expresses the value of $x_n$ as a function of one other sample, $x_{n-1}$. A second order difference equation is a function of two other samples, for example

$x_{n+2} = 1.5x_n + 1.4x_{n+1}, \textit{ for } n = 0, 1, 2, ...,$

## Fibonacci Numbers

The easiest way to describe the Fibonacci sequence is with a second order difference equation. Assume $F_0 = 0$ and $F_1 = 1$

$F_{k+2} = F_{k+1} + F_{k}$

The first eight Fibonacci numbers are:

$F_k=0,1,1,2,3,5,8,13,...$

We can solve the sequence with linear algebra. Let us start our sequence with a vector, $u_0=\begin{bmatrix}1\\0\end{bmatrix}$. Find the next vector in our sequence, $u_1$ like this:

$u_{1} = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \begin{bmatrix} 1 \\ 0 \end{bmatrix} = \begin{bmatrix} 1 \\ 1 \end{bmatrix}$

More generally:

$u_{k+1} = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} u_k$

Note that $u_0$ contains $\begin{bmatrix}F_1\\F_0\end{bmatrix}$.

The key to solving the sequence is in the eigenvalues and eigenvectors of $\begin{bmatrix}1&1\\1&0\end{bmatrix}$. We can write $u_3$ as

$u_3 = \begin{bmatrix}1 & 1 \\ 1 & 0\end{bmatrix} \begin{bmatrix}1 & 1 \\ 1 & 0\end{bmatrix} \begin{bmatrix}1 & 1 \\ 1 & 0\end{bmatrix} \begin{bmatrix}1 \\ 0\end{bmatrix} = \begin{bmatrix}3 \\ 2\end{bmatrix} = \begin{bmatrix}F_4 \\ F_3\end{bmatrix}$

import numpy as np
fib = np.matrix([[1,1],[1,0]])
u_0 = np.matrix([[1],[0]])

fib * fib * fib * u_0 # The * operator is dot product for np.matrix (but not np.array)
matrix([[3],
[2]])

### Eigenvector Solution

The general case of the previous example is:

\begin{align} u_k &= \begin{bmatrix}1 & 1 \\ 1 & 0\end{bmatrix}^k\begin{bmatrix}1 \\ 0\end{bmatrix}\\ u_k &= A^ku_0 \end{align}

We can use eigenvectors and eigenvalues to solve the matrix power (see Eigenvalues and Eigenvectors). First, find the eigenvalues of $A$:

\begin{align} det\begin{vmatrix} 1 - \lambda & 1\\ 1 & - \lambda \end{vmatrix} & = 0\\ -\lambda(1-\lambda) - 1 & = 0\\ \lambda{}^2 - \lambda - 1 & = 0 \end{align}

Solve using the quadratic formula $\frac{-b \pm \sqrt{b^2-4ac}}{2a} \implies \frac{1 \pm \sqrt{5}}{2} \approx \frac{1 \pm 2.236}{2}\\$, which gives the eigenvalues for our matrix: \begin{align} \lambda_1 & \approx &1.618 \\ \lambda_2 & \approx &-0.618 \end{align}

Now we plug our constant eigenvalues ($\lambda_1$ and $\lambda_2$) into $(A-\lambda I)\vec{x}=0$ to find the eigen vectors.

$\begin{bmatrix}1-\lambda_1 & 1 \\ 1 & -\lambda_1 \end{bmatrix}\vec{x}_1 = 0 \implies \vec{x}_1 = \begin{bmatrix}\lambda_1 \\ 1\end{bmatrix}$

Note that this solution is not necessarily intuitive. This is not a valid solution for any arbitrary value of $\lambda$: our values, $\lambda_1$ and $\lambda_2$ are valid solutions. We can use the same process to solve for $\lambda_2$, revealing eigenvectors $\begin{bmatrix}\lambda_1 \\1\end{bmatrix}$ and $\begin{bmatrix}\lambda_2 \\1\end{bmatrix}$. Now create the eigenvector matrix, $S$ and eigenvalue matrix, $\Lambda$.

\begin{align} S &= \begin{bmatrix}\lambda_1 & \lambda_2 \\ 1 & 1\end{bmatrix}\\ \Lambda &= \begin{bmatrix}\lambda_1 & 0 \\ 0 & \lambda_2\end{bmatrix} \end{align}

Remeber that $A^k = S^{-1}\Lambda^kS$. To proceed, invert $S$. Remember this shortcut for inverting a $2 \times 2$ matrix.

\begin{align} \begin{bmatrix} a & b \\ c & d \end{bmatrix}^{-1} =&& \frac{1}{ad-cb} &\begin{bmatrix} d & -b \\ -c & a \end{bmatrix} \\ S^{-1} =&& &\begin{bmatrix}\lambda_1 & \lambda_2 \\ 1 & 1\end{bmatrix}^{-1} \\ S^{-1} =&& \frac{1}{\lambda_1-\lambda_2} &\begin{bmatrix}1 & -\lambda_2\\ -1 & \lambda_1 \end{bmatrix}\\ S^{-1} =&& \tfrac{1}{\sqrt{5}} &\begin{bmatrix}1 & -\lambda_2\\ -1 & \lambda_1 \end{bmatrix} \end{align}

Using the components $(A, \Lambda, S, S^{-1},u_0)$:

\begin{align} u_k &= A^ku_0 \\ u_k &= S\Lambda^kS^{-1}u_o \\ u_k &= S\Lambda^k\frac{1}{\sqrt{5}}\begin{bmatrix}1 & -\lambda_2\\ -1 & \lambda_1 \end{bmatrix}u_0 \\ u_k &= \frac{1}{\sqrt{5}}S\Lambda^k\begin{bmatrix}1 & -\lambda_2\\ -1 & \lambda_1 \end{bmatrix}u_0 \\ u_k &= \tfrac{1}{\sqrt{5}} \begin{bmatrix}\lambda_1 & \lambda_2 \\ 1 & 1\end{bmatrix} \begin{bmatrix}\lambda_1 & 0 \\ 0 & \lambda_2\end{bmatrix}^k \begin{bmatrix}1 & -\lambda_2\\ -1 & \lambda_1 \end{bmatrix} u_0 \\ u_k &= \tfrac{1}{\sqrt{5}} \begin{bmatrix}\lambda_1 & \lambda_2 \\ 1 & 1\end{bmatrix} \begin{bmatrix}\lambda_1^k & 0 \\ 0 & \lambda_2^k\end{bmatrix} \begin{bmatrix}1 & -\lambda_2\\ -1 & \lambda_1 \end{bmatrix} \begin{bmatrix}1\\0\end{bmatrix} \\ u_k &= \tfrac{1}{\sqrt{5}} \begin{bmatrix}\lambda_1 & \lambda_2 \\ 1 & 1\end{bmatrix} \begin{bmatrix}\lambda_1^k & 0 \\ 0 & \lambda_2^k\end{bmatrix} \begin{bmatrix}1 \\ -1 \end{bmatrix}\\ u_k &= \tfrac{1}{\sqrt{5}} \begin{bmatrix}\lambda_1 & \lambda_2 \\ 1 & 1\end{bmatrix} \begin{bmatrix}\lambda_1^k \\ -\lambda_2^k\end{bmatrix} \end{align}

## Verify Solution in Python

import numpy as np

lmbd_1 = (1+np.sqrt(5))/2
lmbd_2 = (1-np.sqrt(5))/2

a       = np.matrix([[1,      1],      [1,  0]])
s       = np.matrix([[lmbd_1, lmbd_2], [1,  1]])
s_inv   = np.matrix([[1,     -lmbd_2], [-1, lmbd_1]]) * (1/np.sqrt(5))
eig_val = np.matrix([[lmbd_1, 0],      [0,  lmbd_2]])
u_0     = np.matrix([[1],              [0]])

def eig_k(k):
return eig_val ** k

def u_k(k):
return s * eig_k(k) * s_inv * u_0

def fib(n):
f_n = u_k(n)[1,0]
return int(round(f_n))

[fib(n) for n in range(15)]
[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377]