How do I calculate the matrix exponential of a sparse matrix?

I’m trying to find the matrix exponential of a sparse matrix:

import numpy as np

b = np.array([[1, 0, 1, 0, 1, 0, 1, 1, 1, 0],
       [1, 0, 0, 0, 1, 1, 0, 1, 1, 0],
       [0, 1, 1, 0, 1, 1, 0, 0, 1, 1],
       [0, 0, 0, 0, 0, 1, 1, 1, 0, 0],
       [1, 1, 0, 0, 0, 0, 1, 0, 0, 0],
       [0, 0, 1, 0, 0, 1, 0, 0, 1, 1],
       [0, 0, 1, 0, 1, 0, 1, 1, 0, 0],
       [1, 0, 0, 0, 1, 1, 0, 0, 1, 1],
       [0, 0, 0, 0, 1, 0, 1, 1, 1, 0],
       [0, 0, 0, 1, 0, 1, 1, 0, 0, 1]])

I can calculate this using scipy.linalg.expm, but it is slow for larger matrices.

from scipy.linalg import expm

S1 = expm(b)

Since this is a sparse matrix, I tried converting b to a scipy.sparse matrix and calling that function on the converted sparse matrix:

import scipy.sparse as sp
import numpy as np

sp_b = sp.csr_matrix(b)
S1 = expm(sp_b);

But I get the following error:

loop of ufunc does not support argument 0 of type csr_matrix which has no callable exp method

How can I calculate the matrix exponential of a sparse matrix?

>Solution :

You need to use scipy.sparse.linalg.expm for your sparse matrix instead of scipy.linalg.expm.

import scipy.sparse as sp
from scipy.sparse.linalg import expm
import numpy as np

b = np.array([[1, 0, 1, 0, 1, 0, 1, 1, 1, 0],
       [1, 0, 0, 0, 1, 1, 0, 1, 1, 0],
       [0, 1, 1, 0, 1, 1, 0, 0, 1, 1],
       [0, 0, 0, 0, 0, 1, 1, 1, 0, 0],
       [1, 1, 0, 0, 0, 0, 1, 0, 0, 0],
       [0, 0, 1, 0, 0, 1, 0, 0, 1, 1],
       [0, 0, 1, 0, 1, 0, 1, 1, 0, 0],
       [1, 0, 0, 0, 1, 1, 0, 0, 1, 1],
       [0, 0, 0, 0, 1, 0, 1, 1, 1, 0],
       [0, 0, 0, 1, 0, 1, 1, 0, 0, 1]])
sp_b = sp.csr_matrix(b)
S1 = expm(sp_b);

Note: As you found, defining your matrix as a CSR matrix gives the warning "SparseEfficiencyWarning: spsolve is more efficient when sparse b is in the CSC matrix format". To get rid of this, you can do as the warning suggests, and define a CSC matrix if that makes sense for your application:

sp_b = sp.csc_matrix(b)

Leave a Reply