Matrix power for sparse matrix in python

linear-algebra, numpy, python, scipy, sparse-matrix

Solution

`**` 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 `dot`s. 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()`.

Problem

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.

Original source