Scipy Sparse Matrix - Dense Multiplication Vector Effect - Blocks versus Large Matrix

Question

I have some scipy sparse matrices (currently in CSR format) that I need to multiply with a dense 1-dimensional numpy vector. The vector is called G :

 print G.shape, G.dtype (2097152,) complex64 

Each sparse matrix has a shape (16384,2097152) and is very rare. The density is about 4.0e-6. I have a list of 100 such sparse matrices called spmats .

I can easily multiply each matrix by G as follows:

 res = [spmat.dot(G) for spmat in spmats] 

This leads to a list of dense shape vectors (16384,) , as expected.

My application is quite critical, and so I tried an alternative, which should first combine all sparse matrices into one large sparse matrix, and then use only one dot() call as follows:

 import scipy.sparse as sp SPMAT = sp.vstack(spmats, format='csr') RES = SPMAT.dot(G) 

This leads to one long RES vector, which has the form (1638400,) and is a concatenated version of all the result vectors in RES above, as expected. I checked that the results are identical.

I may be completely wrong, but I expected the second case to be faster than the first, since there are far fewer numpy calls, memory allocation, python object creation, python loops, etc. I don't care about the time it takes to concatenate sparse matrices, only the time to calculate the result. According to %timeit , however:

 %timeit res = [spmat.dot(G) for spmat in spmats] 10 loops, best of 3: 91.5 ms per loop %timeit RES = SPMAT.dot(G) 1 loops, best of 3: 389 ms per loop 

I checked that I did not have enough memory in any of the operations, and nothing suspicious was happening. Am I crazy or is it really weird? Does this mean that all sparse products with matrix vectors must run in blocks several rows at a time to make them faster? As I understand it, the sparse time of multiplication of a matrix with a dense vector should be linear in the number of nonzero elements, which in both cases has not changed. What could it be?

I work on the same base Linux machine with 4 GB of RAM using EPD7.3

EDIT:

Here is a small example that reproduces the problem for me:

 import scipy.sparse as sp import numpy as n G = n.random.rand(128**3) + 1.0j*n.random.rand(128**3) spmats = [sp.rand (128**2, 128**3, density = 4e-6, format = 'csr', dtype=float64) for i in range(100)] SPMAT = sp.vstack(spmats, format='csr') %timeit res = [spmat.dot(G) for spmat in spmats] %timeit RES = SPMAT.dot(G) 

I get:

 1 loops, best of 3: 704 ms per loop 1 loops, best of 3: 1.34 s per loop 

The difference in performance in this case is not as large as that of my own sparse matrices, which have some structure (possibly due to caching), but even worse for matrix concatenation.

I tried with both scipy 10.1 and 12.0.

Answer 1

I did not find a reason for the strange behavior mentioned in the question, however I found a way to significantly speed up my calculations, which may be useful to other people.

Because in my particular case, I calculate the product of a sparse matrix float32 and a complex 64 dense vector, I can multiply the real and imaginary components separately. This gives me 4x acceleration.

This takes 2.35 with SPMAT.shape == (16384000, 2097152) :

 RES = SPMAT.dot(G) 

So far it only takes 541 ms:

 RES = n.zeros((SPMAT.shape[0],),dtype=complex64) RES.real = SPMAT.dot(G.real); RES.imag = SPMAT.dot(G.imag) 

And the result is identical. I think it is possible that prealocation n.zeros might be optional, but I don't know how to do this.