numpy - Matrix power for sparse matrix in python

Question

I am trying to find out a way to do a matrix power for a sparse matrix M: M^k = M*...*M k times where * is the matrix multiplication (numpy.dot), and not element-wise multiplication.

I know how to do it for a normal matrix:

import numpy as np
import scipy as sp
N=100
k=3
M=(sp.sparse.spdiags(np.ones(N), 0, N, N)-sp.sparse.spdiags(np.ones(N), 2, N, N)).toarray()
np.matrix_power(M,k)

How can I do it for sparse M:

M=(sp.sparse.spdiags(np.ones(N), 0, N, N)-sp.sparse.spdiags(np.ones(N), 2, N, N))

Of course, I can do this by recursive multiplications, but I am wondering if there is a functionality like matrix_power for sparse matrices in scipy. Any help is much much appreciated. Thanks in advance.

Answer 1

** has been implemented for csr_matrix. There is a __pow__ method.

After handling some special cases this __pow__ does:

            tmp = self.__pow__(other//2)
            if (other % 2):
                return self * tmp * tmp
            else:
                return tmp * tmp

For sparse matrix, * is the matrix product (dot for ndarray). So it is doing recursive multiplications.

---

As math noted, np.matrix also implements ** (__pow__) as matrix power. In fact it ends up calling np.linalg.matrix_power.

np.linalg.matrix_power(M, n) is written in Python, so you can easily see what it does.

For n<=3 is just does the repeated dot.

For larger n, it does a binary decomposition to reduce the total number of dots. I assume that means for n=4:

result = np.dot(M,M)
result = np.dot(result,result)

The sparse version isn't as general. It can only handle positive integer powers.

You can't count on numpy functions operating on spare matrices. The ones that do work are the ones that pass the action on to the array's own method. e.g. np.sum(A) calls A.sum().