Recently, I’ve figured out a way to improve matrix exponentiation-like algorithms for sparse matrices. I’m pretty sure that some of you might find this being elementary linear algebra, but I’m pretty proud of my discovery, so I’m going to write a blog post about it. I've marked some TODOs and ideas in square brackets, in case you want to help me enhance this.
Prerequisites
Cayley-Hamilton theorem, Berlekamp-Massey algorithm, fast linear recurrence solvers.
Terminology
- $$$M$$$ is an $$$n \times n$$$ “sparse” matrix, with $$$m$$$ entries different than $$$0$$$.
- $$$k$$$ is the exponent.
Complexities will be written in terms of $$$n, m, k$$$.
Cayley-Hamilton theorem
One of the most astonishing theorems in my opinion, what the theorem says is that for any $$$n \times n$$$ matrix, there exists some coefficients $$$c_i$$$ such that.
$$$M^n = \sum_{i = 1}^n c_i M^{n - i}$$$
$$$c$$$ is often called the characteristic polynomial of $$$M$$$.
This theorem essentially says that the powers of any $$$n \times n$$$ matrix follow a linear recurrence of order at most $$$n$$$ (to understand why that is, take the above equation and multiply both sides with $$$M^x$$$ for any $$$x \geq 0$$$).
Computing quadratic forms $$$\sum{M^k_{i, j} a_i b_j}$$$
Doing some little algebra, it is easy to see that $$$\sum{Q_{i, j} a_i b_j}$$$ can be written equivalently as $$$a^T Q b$$$, for any matrix $$$Q$$$. This is called a quadratic form.
Many problems involving matrix exponentiation in which you have to output one number essentially boil down to these quadratic forms. For example, computing the $$$k$$$-th fibonacci number essentially asks us to compute $$$a^T M^k b$$$, where $$$a = b = [1, 0]^T$$$, and $$$M$$$ is the fibonacci recurrence matrix.
Note that, due to the Cayley-Hamilton theorem, after doing some algebra it can be seen that $$$s_k = a^T M^k b$$$ follows a linear recurrence with the recurrence coefficients equal to $$$c$$$. Therefore, if we compute $$$s_1, …, s_{2n}$$$, we can use Berlekamp-Massey algorithm to deduce $$$c$$$, and then use a fast linear recurrence algorithm to compute $$$X^k mod c$$$ (polynomial modulo) in $$$O(n^2 \log{k})$$$ or $$$O(n \log{n} \log{k})$$$, and then find out $$$s_k$$$ [there was a tutorial on Codeforces for this, perhaps you can share the link in the comments].
However, to do that, we have to compute $$$s_0, …, s_{2n}$$$. In order to do that, remember that $$$s_i = a^T M^i b$$$. This means that $$$s_i = a^T M (M^{i - 1} b)$$$. If we compute $$$M^i b$$$ for $$$i = 1, 2, …, 2n$$$ (using iterative matrix vector products), we can compute $$$s_1, …, s_{2n}$$$ in time $$$O(m n)$$$ (it is essential to understand that a matrix-vector product $$$Mv$$$ can be computed in complexity $$$O(m)$$$).
Summing up, we can do the following:
- Compute $$$s_0, …, s_{2n}$$$ in $$$O(m n)$$$
- Compute $$$c$$$ in $$$O(n^2)$$$ using Berlekamp-Massey
- Compute $$$X^n mod c$$$ in $$$O(n^2 \log{k})$$$ or $$$O(n \log{n}\log{k})$$$
- Compute $$$s_k$$$ in $$$O(n)$$$
Total complexity is $$$O(mn + n^2 + n^2 \log{k} + n)$$$ or $$$O(mn + n^2 + n \log{n} \log{k} + n)$$$.
Computing matrix-vector products $$$M^k v$$$.
Let’s go even further and compute matrix-vector products. This is the same as computing $$$e^{(j)} M^k v$$$ for all $$$j$$$. Here $$$e^{(j)}$$$ is the $$$j$$$-th element of the canonical base ($$$e^{(j)}_k = [k = j]$$$). The key thing to notice here is that the characteristic polynomial $$$c$$$ is the same for all $$$j$$$, so all elements follow the same recurrence. Also, products $$$M^i v$$$ are computed in the previous procedure, so we can just use them to compute all initial terms of the recurrence!
This essentially means that you can just take any non-zero vector $$$w$$$, compute the initial $$$2n$$$ terms $$$M^i v$$$ that you need, and then do Berlekamp-Massey on $$$w^T M^i v$$$, and you are done! The only thing that changes is the last step, which requires $$$O(n^2)$$$ work now.
Total complexity is $$$O(mn + n^2 + n^2 \log{k} + n^2)$$$ or $$$O(mn + n^2 + n \log{n} \log{k} + n^2)$$$.
Computing $$$M^k$$$?
I’m not sure if this is possible, but I’m eager to hear your ideas!
Example problem
Given a directed graph $$$G$$$ with $$$n$$$ vertices and $$$m$$$ edges, compute the number of paths of length exactly $$$k$$$. Optionally, compute the number of paths of length $$$k$$$ starting in each of the $$$n$$$ vertices. $$$n \leq 1000, m \leq 10000, k \leq 10^9$$$. [any links to this problem?]
[other problems?]
Sample code
Note that the sample code does not compute matrix-vector products in $$$O(m)$$$. While it kind of breaks the purpose of this blog post, it is more a test of correctness than an actual solution to a problem.
