Fast Fibonacci Sequence with Matrix Exponentiation (2020-07-12)
In [ ]:
# TITLE: Fast Fibonacci Sequence with Matrix Exponentiation
# COVER: https://i.imgur.com/VM8Y3MM.gif
# DATE: 2020-07-12
# TAGS: algorithms,matrix exponentiation


Matrix exponentiation is a technique that solves a linear recurrence relation in an $O(\log n)$ time, instead of $O(n)$.

A linear recurrence relation is an equation that has the following form:

$$x_n = c_1 x_{n-1} + c_2 x_{n-2} + \ldots + c_k x_{n-k}$$

A well known example of such relation is the Fibonacci sequence:

\begin{align*} F_0 & = 0 \\ F_1 & = 1 \\ F_n & = F_{n-1} + F_{n-2} \end{align*}

So in this post, we're going to use the Fibonacci sequence as an example.

As you can see, to calculate $F_n$, you need to calculate all the previous values in the sequence, which makes its time complexity $O(n)$.

In [1]:
def fib(n):
F_0 = 0
F_1 = 1
for _ in range(n - 1):
F_0, F_1 = F_1, F_0 + F_1
return F_1
%time _ = fib(300_000)

CPU times: user 1.01 s, sys: 1.9 ms, total: 1.02 s
Wall time: 1.02 s


To speed up the above program, all you need to do is the following two things:

• Rewrite the linear recurrence relation into a matrix form (Matrix exponentiation)
• Calculate the exponential part in an $O(\log n)$ time (Divide-and-conquer algorithm)

Matrix Exponentiation¶

Recall the following equation of a linear recurrence relation:

$$x_n = c_1 x_{n-1} + c_2 x_{n-2} + \ldots + c_k x_{n-k}$$

The goal here is to write an equation such that the following can be true: \begin{align*} \begin{bmatrix} x_n \\ x_{n-1} \\ x_{n-2} \\ \vdots \\ x_{n-k} \end{bmatrix} &= C \begin{bmatrix} x_{n-1} \\ x_{n-2} \\ x_{n-3} \\ \vdots \\ x_{n-(k+1)} \end{bmatrix} \end{align*}

Getting the matrix $C$ above is pretty trivial since its relation is linear. A linear recurrence relation can be rewritten with matrices as follows:

\begin{align*} \begin{bmatrix} x_n \\ x_{n-1} \\ x_{n-2} \\ \vdots \\ x_{n-k} \end{bmatrix} &= \begin{bmatrix} c_1 & c_2 & \ldots & c_k \\ 1 & 0 & \ldots & 0 \\ 0 & 1 & \ldots & 0 \\ \vdots & \vdots & \ldots & \vdots \\ 0 & 0 & \ldots & 1 \end{bmatrix} \begin{bmatrix} x_{n-1} \\ x_{n-2} \\ x_{n-3} \\ \vdots \\ x_{n-(k+1)} \end{bmatrix} \end{align*}

Since now we have a matrix that makes our result matrix one $n$ forward, the below equation becomes also true:

\begin{align*} \begin{bmatrix} x_n \\ x_{n-1} \\ x_{n-2} \\ \vdots \\ x_{n-k} \end{bmatrix} &= \begin{bmatrix} c_1 & c_2 & \ldots & c_k \\ 1 & 0 & \ldots & 0 \\ 0 & 1 & \ldots & 0 \\ \vdots & \vdots & \ldots & \vdots \\ 0 & 0 & \ldots & 1 \end{bmatrix}^{\text{ }2} \begin{bmatrix} x_{n-2} \\ x_{n-3} \\ x_{n-4} \\ \vdots \\ x_{n-(k+2)} \end{bmatrix} \\ &= \begin{bmatrix} c_1 & c_2 & \ldots & c_k \\ 1 & 0 & \ldots & 0 \\ 0 & 1 & \ldots & 0 \\ \vdots & \vdots & \ldots & \vdots \\ 0 & 0 & \ldots & 1 \end{bmatrix}^{\text{ }3} \begin{bmatrix} x_{n-3} \\ x_{n-4} \\ x_{n-5} \\ \vdots \\ x_{n-(k+3)} \end{bmatrix} \\ & \vdots \\ &= \begin{bmatrix} c_1 & c_2 & \ldots & c_k \\ 1 & 0 & \ldots & 0 \\ 0 & 1 & \ldots & 0 \\ \vdots & \vdots & \ldots & \vdots \\ 0 & 0 & \ldots & 1 \end{bmatrix}^{\text{ }n} \begin{bmatrix} x_{k} \\ x_{k-1} \\ x_{k-2} \\ \vdots \\ x_{1} \end{bmatrix} \end{align*}

In a linear recurrence relation, the necessary initial values ${x_1, \ldots, x_k}$ are given. So if we know the value of $C^n$, we no longer require the previous values!

So if we plug our Fibonacci sequence example with the above, we get:

\begin{align*} \begin{bmatrix} F_n \\ F_{n-1} \end{bmatrix} &= \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \begin{bmatrix} F_{n-1} \\ F_{n-2} \end{bmatrix} \\ &= \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}^{\text{ }n} \begin{bmatrix} 0 \\ 1 \end{bmatrix} \end{align*}

But as you can see, we still have to compute $n$ times to calculate $C^n$. That's why we have to talk about a way to calculate an exponential in a logarithmic time by using a divide-and-conquer algorithm.

Fast Exponentiation¶

The idea is very simple. Since matrix multiplication has associativity, we can divide the calculation as follows:

$$a^b = \begin{cases} 1 & \text{if b = 0} \\ a^{b/2} a^{b/2} & \text{if b is an even number} \\ a^{\lfloor b/2 \rfloor} a^{\lfloor b/2 \rfloor} a & \text{if b is an odd number} \end{cases}$$

With the above, instead of doing the multiplication $n$ times, you can do it only $\log_{2}{n}$ times, which makes it the time complexity of $O(\log n)$.

Here's the Python version of the above equation:

In [2]:
def pow(a, b):
if b == 0:
c = 1
else:
v = pow(a, b // 2)
c = v * v
if b % 2 == 1:
c *= a
return c

pow(2, 10)

Out[2]:
1024

Fast Fibonacci Sequence¶

Now that we have all the tools we need in our hands, let's write a program that calculates the $n$th Fibonacci sequence number in an $O(\log n)$ time.

Let's first define functions to calculate matrices in the equation. Of course, if you're not in a competitive programming setup, you can just use numpy.

In [3]:
def matmul21(A, B):
return [
[A[0][0] * B[0][0] + A[0][1] * B[1][0]],
[A[1][0] * B[0][0] + A[1][1] * B[1][0]]
]

def matmul22(A, B):
return [
[A[0][0] * B[0][0] + A[0][1] * B[1][0], A[0][0] * B[0][1] + A[0][1] * B[1][1]],
[A[1][0] * B[0][0] + A[1][1] * B[1][0], A[1][0] * B[0][1] + A[1][1] * B[1][1]]
]


Let's write a power function that uses a divide-and-conquer algorithm:

In [4]:
def matpow22(a, b):
if b == 0:
c = [[1, 0], [0, 1]] # 2x2 Identity matrix
else:
v = matpow22(a, b // 2)
c = matmul22(v, v)
if b % 2 == 1:
c = matmul22(c, a)
return c


And now the Fibonacci sequence computation part is just simply:

In [5]:
C         = [[1, 1],
[1, 0]]
F_initial = [[0],
[1]]

def fib_fast(n):
return matmul21(matpow22(C, n), F_initial)[0][0]

%time _ = fib_fast(300_000)

CPU times: user 45.5 ms, sys: 1.04 ms, total: 46.6 ms
Wall time: 45.5 ms


~1100ms -> ~50ms!