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)$.

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:

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.

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

And now the Fibonacci sequence computation part is just simply:

~1100ms -> ~50ms!