Vectorization of index operation for scipy.sparse matrix

The following code is too slow, although everything seems to be vectorized.

from numpy import * from scipy.sparse import * n = 100000; i = xrange(n); j = xrange(n); data = ones(n); A=csr_matrix((data,(i,j))); x = A[i,j] 

The problem is that the indexing operation is implemented as a python function, and calling A[i,j] leads to the following profiled output

  500033 function calls in 8.718 CPU seconds Ordered by: internal time ncalls tottime percall cumtime percall filename:lineno(function) 100000 7.933 0.000 8.156 0.000 csr.py:265(_get_single_element) 1 0.271 0.271 8.705 8.705 csr.py:177(__getitem__) (...) 

Namely, the python function _get_single_element gets called 100,000 times, which is really inefficient. Why is this not implemented in pure C? Does anyone know a way to overcome this limitation and speed up the above code? Should I use a different sparse matrix type?

+6
python scipy indexing sparse-matrix
source share
2 answers

You can use A.diagonal() to get the diagonal faster (0.0009 seconds versus 3.8 seconds on my machine). However, if you want to do double indexing, this is a more complicated question because you do not use slices in the same way as a list of indexes. The _get_single_element function is called 100,000 times because it simply iterates through the iterators (i and j) that you passed to it. The piece will be A [30: 60,10] or something like that.

Also, I would use csr_matrix(eye(n,n)) to make the csr_matrix(eye(n,n)) same matrix you did with iterators, just for simplicity.

Update:

Well, since your question really is being able to quickly access many random elements, I will answer your questions as best as possible.

  • Why is this not implemented in pure C?

The answer is simple: no one turned to him. There is still a lot of work in the field of Scipy sparse matrix modules from what I saw. One part implemented in C is transformations between different sparse matrix formats.

  • Does anyone know a way to overcome this limitation and speed up the above code?

You can try to dive into sparse matrix modules and try to speed them up. I did this and was able to reduce the time to less than one third of the original when you tried your code above for random access using csr matrices. I had to directly access _get_single_element and significantly shorten the code to do this, including removing related checks.

However, it was even faster to use lil_matrix (although it was slower to initialize the matrix), but I had to do list access because the lil matrices are not configured for the type of indexing you are doing. Using list view for csr_matrix still leaves the lil matrix method. Ultimately, the lil matrix is ​​faster for accessing random elements because it is not compressed.

Using lil_matrix in the source form is done in about a fifth of the code above. If I select some related checks and call the lil_matrix _get1 () method directly, I can reduce the time by about 7% from the original time. For clarity, the acceleration is from 3.4-3.8 seconds to 0.261 seconds.

Finally, I tried to make my own function, which accesses lil matrix data directly and avoids function callbacks. The time for this was about 0.136 seconds. This did not help to sort the data, which is another potential optimization (in particular, if you get access to many elements that are on the same lines).

If you want faster than this, you probably have to write your own implementation of a sparse matrix of C code.

  • Should I use a different sparse matrix type?

Well, I suggest the lil matrix if you intend to access many elements, but it all depends on what you need to do. Do you also need to multiply matrices, for example? Just remember that a change between matrices can at least sometimes (under certain circumstances) be quite quick, so do not rule out changes in another matrix format to perform different operations.

If you don’t need to do any algebra operations on your matrix at all, maybe you should just use defaultdict or something like that. The danger with defaultdicts is that whenever an element is requested that it is not in a dict, it sets that element to default and saves it in such a way that it can be problematic.

+7
source share

I think _get_single_element is only called when the default dtype for 'object' is used. Have you tried providing a dtype, e.g. csr_matrix((data, (i,j)), dtype=int32)

0
source share

All Articles