Efficiency of scipy.sparse.linalg.expm_multiply with sparse vs unsparse vectors

Question

From the package scipy.sparse.linalg in Python, calling expm_multiply(X, v) allows you to compute the vector expm(X)@v in a very efficient manner. Moreover, this function works best when X is a sparse datatype (like csr_matrix or csc_matrix). However, I have found that this function works slower when v is a sparse data type, compared with when v is a numpy.ndarray (see example below).

This leads to a few questions. First: why is expm_multiply(X, v) most efficient when v is not a sparse datatype, even though X should be sparse? Second: should I think of this as a general property - in the sense that I should not convert vectors into sparse datatypes in general if efficiency is important?

---

Some example code that shows what happens is below:

import numpy as np
from numpy.random import uniform
from time import process_time
from scipy import linalg
from scipy import sparse
from scipy.sparse import linalg

dim = pow(2, 12)

#create vector where most elements are 0
v = np.zeros(dim)
for i in range(len(v)):
    if i%50 == 0:
        v[i] = uniform(-1, 1)
        
#create matrix where most elements are 0
X = np.zeros([dim,dim])
for i in range(len(v)):
    for j in range(len(v)):
        if i%50 ==0 and (j+1)%50 == 0:
            X[i][j] = uniform(-1, 1)
            
Xsparse = sparse.csr_matrix(X)
vsparse = sparse.csr_matrix(v).T

#compare sparse vs unsparse exp(X)v calculation
start = process_time()
A = sparse.linalg.expm_multiply(Xsparse, v)
print(process_time() - start)

start = process_time()
B = sparse.linalg.expm_multiply(Xsparse, vsparse)
print(process_time() - start)

This will yield an output like:

0.0019309999999990168
0.01379999999999626

So using vsparse takes approximately an order of magnitude longer than using v.

Answer 1

Usually, not always but usually, even if both $A$ and $x$ are sparse, $Ax$ is not. Even when it is, it is denser than $x$. If you consider something like $e^Ax$, which can be rewritten as $(I+A+A^2/2+\dots)x$, there is no guarantee that it will have the same sparsity pattern or the same number of nonzeros as $x$. Which means that either scipy would have to reallocate memory constantly or would have to convert the sparse vector into a dense one. In such small examples, that could account for the order of magnitude difference.