Fibonacci Numbers

Introduction

Fibonacci Numbers form a sequence in which each number is the sum of the two preceding ones, defined recursively as: \[\begin{equation*} \begin{aligned} &F(0) = 0\\ &F(1) = 1\\ &F(n) = F(n-1) + F(n-2) \end{aligned} \end{equation*}\]

We would like to compute \(F(n)\) in the most efficient way possible!

Dynamic Programming

You can solve it in O(n) time and O(n) space via dynamic programming.

def solve(n):
  f = [0] * (n+2)
  f[0], f[1] = 0, 1
  for i in range(2, n + 1):
    f[i] = f[i-1] + f[i-2]
  return f[n]

Can also be improved to O(n) time and O(1) space, as you only need to keep track of the last two numbers:

def solve(n):
  f = [0, 1]
  for i in range(2, n + 1):
    f[i % 2] = f[(i+1) % 2] + f[(i) % 2]
  return f[n % 2]

Linear Algebra - O(nlogn)

The second DP approach we implemented above, can be implemented via a matrix and a vector as: \[\begin{bmatrix} 1 & 0 \\ 1 & 1 \\ \end{bmatrix} \begin{bmatrix} F_n\\ F_{n-1} \end{bmatrix} = \begin{bmatrix} F_n + F_{n-1}\\ F_n \end{bmatrix} = \begin{bmatrix} F_{n+1}\\ F_n \end{bmatrix}\]

We can start with vector v = [1, 0] and multiply it with the matrix above n-1 times. The result would be the vector [F(n), F(n-1)]. So, we just read the first element of the vector. Thanks to numpy, we can use numpy array and then the operator @ does the matrix / vector multiplication for us.

Also, since \(M*(M*v)\) is the same as \((M*M)*v\), we can first multiply the matrix n-1 times and then multiply with the vector v. So, the answer would be: \(M^{(n-1)} * v\)

import numpy as np

def solve(n):
  mat = np.array([[1, 1], [1, 0]])
  v = np.array([1, 0])
  m = np.identity(2, dtype=int)
  for i in range(n-1):
    m = m @ mat
  return (m @ v)[0]

So far, it’s still O(n) but we’re gonna use the fast exponention trick to speed it up. The idea is: \(M^{100} = M^{50} * M^{50}\) and \(M^{101} = M^{50} * M^{50} * M\), so we only need to compute \(M^{50}\), for which we need to compute \(M^{25}\) and so on. Hence, we arrive at an algorithm with runtime complexity of O(nlogn):

import math
import numpy as np
def pow(m, p):
  I = np.identity(2)
  if p == 0:
    return I
  half = pow(m, p // 2)
  return half @ half @ (I if p%2==0 else m)

def solve(n):
  mat = np.array([[1, 1], [1, 0]])
  v = np.array([1, 0])
  m = pow(mat, n-1)
  return (m @ v)[0]

Linear Algebra - O(1)

Now, let’s see how we can compute it in O(1) using Matrix Eigen Decomposition.

We know from Linear Algebra that for any Matrix \(M\) with eigen vector \(v_i\) and the corresponding eighen value \(\lambda_i\), we have: \(M * v_i = \lambda_i * v_i\)

So, if we compute the eigen vectors (v1, v2) and eighen values (\(\lambda_1\), \(\lambda_2\)) of our matrix \(M\), we can create a matrix V of eighen vectors as columns: \(V = \begin{bmatrix} v1_0 & v2_0\\ v1_1 & v2_1\\ \end{bmatrix}\)

and a matrix \(\Lambda\) of eigen values on the main diagonal: \[\Lambda = \begin{bmatrix} \lambda_1 & 0\\ 0 & \lambda_2\\ \end{bmatrix}\]

And we have: \[M * V = V * \Lambda\]

Multiply both sides by \(V^{-1}\): \[M * V * V^{-1} = V * \Lambda * V^{-1} \\ M = V * \Lambda * V^{-1}\]

You may ask: “What’s special about this representation?” Let’s take a look at what happens when we multiply M by itself: \[\begin{equation*} \begin{aligned} M * M &= V * \Lambda * V^{-1} * V * \Lambda * V^{-1} \\ &= V * \Lambda * \Lambda * V^{-1} \\ &= V * \Lambda^{2} * V^{-1} \end{aligned} \end{equation*}\]

And in general: \(M^n = V * \Lambda^{n} * V^{-1}\)

So, we reduced the problem of \(M^n\) to \(\Lambda^n\) and that’s easier to compute because \(\Lambda\) is a diagonal matrix and as you may know, to compute the n-th power of a diagonal matrix, you can simply raise every element on the diagonal to the power of n. So: \[\begin{bmatrix} \lambda_1 & 0\\ 0 & \lambda_2\\ \end{bmatrix}^n = \begin{bmatrix} \lambda_1^n & 0\\ 0 & \lambda_2^n\\ \end{bmatrix}\]

Implementation

In order to compute the eighen vector/values and compute the inverse matrix, you can grab a pen & paper and roll up your sleeves. But you can also use WolframAlpha to compute these.

import math
import numpy as np

def solve(n):
  s5 = math.sqrt(5)
  lam1, lam2 = (1 + s5)/2, (1 - s5)/2
  v1, v2 = [lam1, 1], [lam2, 1]
  V = np.array([v1, v2]).transpose()
  V_inv = (1/s5) * np.array([[1, (s5-1)/2], [-1, (s5+1)/2]])
  L = np.array([[math.pow(lam1, n-1), 0], [0, math.pow(lam2, n-1)]])
  v = np.array([1, 0])
  return round((V @ L @ V_inv @ v)[0])

Complexity

Assuming raising a number to the power of n is done in constant time, the function above runs in constant time AKA O(1).

Closed Form Formula

Given what we have from the previous section, we can easily arrive at the closed form formula by computing \(M^{n-1}*v\): \[\begin{equation*} \begin{aligned} M^n * v &= V * \Lambda^{n} * V^{-1} * v\\ &= V * \begin{bmatrix} \lambda_1^{n-1} & 0\\ 0 & \lambda_2^{n-1}\\ \end{bmatrix} * V^{-1} * v \\ &= \begin{bmatrix} \frac{1 + \sqrt{5}}{2} & \frac{1-\sqrt{5}}{2}\\ 1 & 1\\ \end{bmatrix} * \begin{bmatrix} (\frac{1+\sqrt{5}}{2})^{n-1} & 0\\ 0 & (\frac{1-\sqrt{5}}{2})^{n-1}\\ \end{bmatrix} * \begin{bmatrix} \frac{1}{\sqrt{5}} & \frac{\sqrt{5}-1}{2\sqrt{5}}\\ \frac{-1}{\sqrt{5}} & \frac{\sqrt{5}+1}{2\sqrt{5}}\\ \end{bmatrix} * \begin{bmatrix} 1\\ 0\\ \end{bmatrix}\\ &= \frac{1}{\sqrt{5}} \left[ \left( \frac{1+\sqrt{5}}{2} \right ) ^ n - \left( \frac{1-\sqrt{5}}{2} \right ) ^ n \right] \end{aligned} \end{equation*}\]

So, the closed form is: \[\begin{equation*} \begin{aligned} F(n) = \frac{1}{\sqrt{5}} \left[ \left( \frac{1+\sqrt{5}}{2} \right ) ^ n - \left( \frac{1-\sqrt{5}}{2} \right ) ^ n \right] \end{aligned} \end{equation*}\]

Nima Aghdaii