1
$\begingroup$

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.

$\endgroup$
4
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.