¿Por qué es más lenta la indexación de filas de matrices CSR scipy en comparación con matrices numpy?

No estoy seguro de lo que estoy haciendo mal, pero parece que la indexación de filas de un csr_matrix csr_matrix es aproximadamente 2 veces más lenta en comparación con las matrices numpy (consulte el código a continuación).

¿No debería ser más rápida la indexación de filas de las matrices rsr que las matrices densas porque solo se extraen pocos elementos distintos de cero como en el caso que se muestra a continuación?

¿Hay trucos para hacer que la indexación de filas sea más rápida para matrices CSR scipy?

 import numpy as np import timeit from scipy.sparse import csr_matrix # Generate random matrix A = np.random.rand(5000, 1000) # Make A sparse A[:, 4:] =0 # Create sparse matrix A_sparse = csr_matrix(A) # Create row indexing functions def index_A_dense(): A[4] def index_A_dense_copy(): a = A[4].copy() def index_A_sparse(): A_sparse[4] timeit.timeit(index_A_sparse, number=10000) Out: 0.6668063667304978 timeit.timeit(index_A_dense, number=10000) Out: 0.0029007720424942818 timeit.timeit(index_A_dense_copy, number=10000) Out: 0.00979283193282754 

¡Gracias por adelantado!

La respuesta corta que demuestro a continuación es que la construcción de una nueva matriz dispersa es costosa. Hay una sobrecarga significativa que no depende del número de filas o del número en elementos distintos de cero en una fila en particular.


La representación de datos para matrices dispersas es bastante diferente de la de matriz densa. Las matrices almacenan los datos en un búfer contiguo, y usan de manera eficiente la shape y los strides para iterar sobre los valores seleccionados. Esos valores, más el índice, definen exactamente si estaban en el búfer donde encontrarán los datos. Copiar esos N bytes de una ubicación a otra es una parte relativamente menor de toda la operación.

Una matriz dispersa almacena los datos en varias matrices (u otras estructuras), que contienen los índices y los datos. La selección de una fila requiere buscar los índices relevantes y construir una nueva matriz dispersa con los índices y datos seleccionados. Hay un código comstackdo en el paquete disperso, pero no tanto como el código de bajo nivel que con las matrices numpy.

Para ilustrar haré una matriz pequeña, y no tan densa, por lo que no tenemos muchas filas vacías:

 In [259]: A = (sparse.rand(5,5,.4,'csr')*20).floor() In [260]: A Out[260]: <5x5 sparse matrix of type '' with 10 stored elements in Compressed Sparse Row format> 

El equivalente denso, y una copia de la fila:

 In [262]: Ad=AA In [263]: Ad Out[263]: array([[ 0., 0., 0., 0., 10.], [ 0., 0., 0., 0., 0.], [ 17., 16., 14., 19., 6.], [ 0., 0., 1., 0., 0.], [ 14., 0., 9., 0., 0.]]) In [264]: Ad[4,:] Out[264]: array([ 14., 0., 9., 0., 0.]) In [265]: timeit Ad[4,:].copy() 100000 loops, best of 3: 4.58 µs per loop 

Una fila de matrices:

 In [266]: A[4,:] Out[266]: <1x5 sparse matrix of type '' with 2 stored elements in Compressed Sparse Row format> 

Mire la representación de datos para esta matriz csr , 3 matrices 1d:

 In [267]: A.data Out[267]: array([ 0., 10., 17., 16., 14., 19., 6., 1., 14., 9.]) In [268]: A.indices Out[268]: array([3, 4, 0, 1, 2, 3, 4, 2, 0, 2], dtype=int32) In [269]: A.indptr Out[269]: array([ 0, 2, 2, 7, 8, 10], dtype=int32) 

Así es como se selecciona la fila (pero en el código comstackdo):

 In [270]: A.indices[A.indptr[4]:A.indptr[5]] Out[270]: array([0, 2], dtype=int32) In [271]: A.data[A.indptr[4]:A.indptr[5]] Out[271]: array([ 14., 9.]) 

La ‘fila’ es otra matriz dispersa, con el mismo tipo de matrices de datos:

 In [272]: A[4,:].indptr Out[272]: array([0, 2]) In [273]: A[4,:].indices Out[273]: array([0, 2]) In [274]: timeit A[4,:] 

Sí, el tiempo para la matriz dispersa es lento. No sé cuánto tiempo se dedica a seleccionar los datos, y cuánto se gasta en la construcción de la nueva matriz.

 10000 loops, best of 3: 145 µs per loop In [275]: timeit Ad[4,:].copy() 100000 loops, best of 3: 4.56 µs per loop 

lil formato lil puede ser más fácil de entender, ya que los datos y los índices se almacenan en sublistas, uno por fila.

 In [276]: Al=A.tolil() In [277]: Al.data Out[277]: array([[0.0, 10.0], [], [17.0, 16.0, 14.0, 19.0, 6.0], [1.0], [14.0, 9.0]], dtype=object) In [278]: Al.rows Out[278]: array([[3, 4], [], [0, 1, 2, 3, 4], [2], [0, 2]], dtype=object) In [279]: Al[4,:].data Out[279]: array([[14.0, 9.0]], dtype=object) In [280]: Al[4,:].rows Out[280]: array([[0, 2]], dtype=object) 

Las comparaciones de velocidad como esta tienen algún sentido cuando se trata de código comstackdo apretado, donde los movimientos de bytes de parte de la memoria a otra son importantes consumidores de tiempo. Con la mezcla de Python y el código comstackdo en numpy y scipy no puedes contar las operaciones O(n) .

=============================

Aquí hay una estimación del tiempo que lleva seleccionar una fila de A , y el tiempo que toma devolver una nueva matriz dispersa:

Solo obtén los datos:

 In [292]: %%timeit d1=A.data[A.indptr[4]:A.indptr[5]] i1=A.indices[A.indptr[4]:A.indptr[5]] .....: 100000 loops, best of 3: 4.92 µs per loop 

más el tiempo que lleva hacer una matriz:

 In [293]: %%timeit d1=A.data[A.indptr[4]:A.indptr[5]] i1=A.indices[A.indptr[4]:A.indptr[5]] sparse.csr_matrix((d1,([0,0],i1)),shape=(1,5)) .....: 1000 loops, best of 3: 445 µs per loop 

Prueba una matriz de coo más simple

 In [294]: %%timeit d1=A.data[A.indptr[4]:A.indptr[5]] i1=A.indices[A.indptr[4]:A.indptr[5]] sparse.coo_matrix((d1,([0,0],i1)),shape=(1,5)) .....: 10000 loops, best of 3: 135 µs per loop