Calculating Fibonacci Numbers with Matrices and Linear Algebra

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.

an=1n, for n=0,1,2,..., 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.

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

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

x0=100x1=1.5×100=150x2=1.5×150=225x3=1.5×225=337.5 \begin{aligned} 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{aligned}

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

x3=100×1.5×1.5×1.5=337.5x_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 nn without first evaluating n=1,2,etcn=1,2,etc. The solution to the equation above is:

xn=1.5n×x0, for n=0,1,2,..., 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 xnx_n as a function of one other sample, xn−1x_{n-1}. A second order difference equation is a function of two other samples, for example

xn+2=1.5xn+1.4xn+1, for n=0,1,2,..., 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 F0=0F_0 = 0 and F1=1F_1 = 1

Fk+2=Fk+1+FkF_{k+2} = F_{k+1} + F_{k}

The first eight Fibonacci numbers are:

Fk=0,1,1,2,3,5,8,13,...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, u0=[10]u_0=\begin{bmatrix}1\\0\end{bmatrix}. Find the next vector in our sequence, u1u_1 like this:

u1=[1110][10]=[11] u_{1} = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \begin{bmatrix} 1 \\ 0 \end{bmatrix} = \begin{bmatrix} 1 \\ 1 \end{bmatrix}

More generally:

uk+1=[1110]uk u_{k+1} = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} u_k

Note that u0u_0 contains [F1F0]\begin{bmatrix}F_1\\F_0\end{bmatrix}.

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

u3=[1110][1110][1110][10]=[32]=[F4F3] 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:

uk=[1110]k[10]uk=Aku0 \begin{aligned} u_k &= \begin{bmatrix}1 & 1 \\ 1 & 0\end{bmatrix}^k\begin{bmatrix}1 \\ 0\end{bmatrix}\\ u_k &= A^ku_0 \end{aligned}

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

det∣1−λ11−λ∣=0−λ(1−λ)−1=0λ2−λ−1=0 \begin{aligned} det\begin{vmatrix} 1 - \lambda & 1\\ 1 & - \lambda \end{vmatrix} & = 0\\ -\lambda(1-\lambda) - 1 & = 0\\ \lambda{}^2 - \lambda - 1 & = 0 \end{aligned}

Solve using the quadratic formula −b±b2−4ac2a  ⟹  1±52≈1±2.2362\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:

λ1≈1.618λ2≈−0.618 \begin{aligned} \lambda_1 & \approx &1.618 \\ \lambda_2 & \approx &-0.618 \end{aligned}

Now we plug our constant eigenvalues (λ1\lambda_1 and λ2\lambda_2) into (A−λI)x⃗=0(A-\lambda I)\vec{x}=0 to find the eigenvectors.

[1−λ111−λ1]x⃗1=0  ⟹  x⃗1=[λ11] \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, λ1\lambda_1 and λ2\lambda_2 are valid solutions. We can use the same process to solve for λ2\lambda_2, revealing eigenvectors [λ11]\begin{bmatrix}\lambda_1 \\1\end{bmatrix} and [λ21]\begin{bmatrix}\lambda_2 \\1\end{bmatrix}. Now create the eigenvector matrix, SS and eigenvalue matrix, Λ\Lambda.

S=[λ1λ211]Λ=[λ100λ2] \begin{aligned} S &= \begin{bmatrix}\lambda_1 & \lambda_2 \\ 1 & 1\end{bmatrix}\\ \Lambda &= \begin{bmatrix}\lambda_1 & 0 \\ 0 & \lambda_2\end{bmatrix} \end{aligned}

Remember that Ak=S−1ΛkSA^k = S^{-1}\Lambda^kS. To proceed, invert SS. Remember this shortcut for inverting a 2×22 \times 2 matrix.

[abcd]−1=1ad−cb[d−b−ca]S−1=[λ1λ211]−1S−1=1λ1−λ2[1−λ2−1λ1]S−1=15[1−λ2−1λ1] \begin{aligned} \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{aligned}

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

uk=Aku0uk=SΛkS−1u0uk=SΛk15[1−λ2−1λ1]u0uk=15SΛk[1−λ2−1λ1]u0uk=15[λ1λ211][λ100λ2]k[1−λ2−1λ1]u0uk=15[λ1λ211][λ1k00λ2k][1−λ2−1λ1][10]uk=15[λ1λ211][λ1k00λ2k][1−1]uk=15[λ1λ211][λ1k−λ2k] \begin{aligned} u_k &= A^ku_0 \\ u_k &= S\Lambda^kS^{-1}u_0 \\ 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{aligned}

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]