Raise The Power Of A Sparse Matrix

I have a sparse matrix of 10001 rows + 10001 columns (with many 0's), I am trying to raise the power of this sparse matrix i.e. A = [[1,1],[1,0]] AS = sparse.csr_matrix(A) AS def

Solution 1:

You can use the power operator: AS ** 10 is the same as AS ^ 10 in normal mathematic notation.

def matrixMul(AS, n):
    if(n <= 1):
        returnASelse:
        returnAS ** n

Solution 2:

In [3]: from scipy import sparse
In [4]: A = [[1,1],[1,0]]
   ...: AS = sparse.csr_matrix(A)
In [5]: AS
Out[5]: 
<2x2 sparse matrix of type'<class 'numpy.int64'>'
    with 3 stored elements in Compressed Sparse Row format>
In [6]: AS.A
Out[6]: 
array([[1, 1],
       [1, 0]])
In [7]: (AS@AS@AS).A
Out[7]: 
array([[3, 2],
       [2, 1]])

@ maps to AS.__matmul__, but np.matmul does not delegate it that way.

In [8]: np.matmul(AS,AS)
Traceback (most recent call last):
  File "<ipython-input-8-37972b025121>", line 1, in <module>
    np.matmul(AS,AS)
ValueError: matmul: Input operand 0 does not have enough dimensions (has 0, gufunc core with signature (n?,k),(k,m?)->(n?,m?) requires 1)

Correction:

In [9]: def matrixMul(AS, n):
   ...:     if(n <= 1):
   ...:         return AS
   ...:     else:
   ...:         return matrixMul(AS, n-1)@ AS
   ...: 
In [10]: matrixMul(AS,3)
Out[10]: 
<2x2 sparse matrix of type'<class 'numpy.int64'>'
    with 4 stored elements in Compressed Sparse Row format>
In [11]: _.A
Out[11]: 
array([[3, 2],
       [2, 1]])
In [12]: matrixMul(AS,10).A
Out[12]: 
array([[89, 55],
       [55, 34]])

As shown in comments, the ** power works:

In [15]: (AS**10).A
Out[15]: 
array([[89, 55],
       [55, 34]])

Sparse matrix is modeled on np.matrix. * is matrix multiplication, and ** is matrix power.

This power times about the same as this semi-smart explicit multiplication:

def foo(AS):
    AS2=AS@ASreturn AS2@AS2@AS2@AS2@AS2
In [33]: timeit foo(AS).A
865 µs ± 378 ns per loop(mean ± std. dev. of 7 runs, 1000 loops each)
In [34]: timeit AS**10767 µs ± 189 ns per loop(mean ± std. dev. of 7 runs, 1000 loops each)