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

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.

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.