Advantages of scipy.linalg.expm over scipy.sparse.linalg.expm?

user_2789788 0
Viewed 147

I'm having trouble understanding when it would be better to use the sparse version of expm over the regular version.

If I take a dense matrix M,which happens to be sparse, of size 4000 x 4000 matrix and compute scipy.linalg.expm(M), this takes around 25 seconds. If I then do M_sparse = sparse.csc_matrix(M) and compute scipy.sparse.linalg.expm(M_sparse), this runs for at least 10 minutes.

What is the right use case for scipy.sparse.linalg.expm? Is it intended for much larger matrices?

scipy sparse-matrix python
2 Answers

scipy.linalg.expm uses sparse.linalg.expm. Look at its source:

Signature: linalg.expm(A)
Source:   
def expm(A):
    import scipy.sparse.linalg
    return scipy.sparse.linalg.expm(A)

We can look at the sparse source as well, though that's a lot longer, since it's doing the real work. It also handles a number of different cases, including varying degrees of Pade approximation. You may have to read up on that approximation. If I recall correctly it involves an inverse.

When you say your array is 'sparse', what exactly do you mean? Already as scipy.sparse.csr matrix? Or a ndarray with many 0s? What proportion of 0s (or nonzeros)? 50% or 5% sparsity? Sparsity matters a lot when working with scipy.sparse matrices. And (4000,4000) is small by most sparse usage. A dense version will fit in memory, and usually calculate faster.

My guess is that if there's a time difference, it's because of the usage of a sparse matrix versus ndarray, not the calling function.

expm on a (4000,4000) array is fairly slow however it's done.

user_9019250

The reason for common algorithm for sparse and dense is just historical.

The situation will change in version 1.9 via https://github.com/scipy/scipy/pull/15079 and then sparse arrays will use the sparse.linalg version and the dense will go through the linalg.expm. The typical difference is that you can't supply, say, 25000x25000 dense matrix without exceptionally large memory hardware but if there is sparsity in your array and create a sparse array that fits in the memory then you can use the sparse version that doesn't use dense algorithm, in particular, matrix powers, norm estimation and matvec operations.

With this PR, it will accelarate the dense version further especially smaller matrices up to 10x and also will handle stacked numpy arrays.

The plan is to overhaul all matrix functions, solve and eig such that the API is unified and usability is increased.

user_178210620
Related