slicing sparse (scipy) matrix

Question

Agradecería cualquier ayuda, para comprender el siguiente comportamiento al cortar una lil_matrix (A) del paquete scipy.sparse.

En realidad, me gustaría extraer una submatriz basada en una lista de índices arbitrarios para filas y columnas.

Cuando utiliza estas dos líneas de código:

x1 = A[list 1,:] 
x2 = x1[:,list 2] 

Todo fue bien y pude extraer la submatriz derecha.

Cuando traté de hacer esto en una sola línea, que fracasaron (La matriz de regresar estaba vacía)

x=A[list 1,list 2] 

¿Por qué es esto así? En general, he usado un comando similar en matlab y funciona. Entonces, ¿por qué no usar el primero, ya que funciona? Parece ser bastante lento. Como tengo que pasar por una gran cantidad de entradas, me gustaría acelerarlo con un solo comando. Tal vez use el tipo de matriz dispersa incorrecta ... ¿Alguna idea?

Answer 1

slicing sucede con esta sintaxis:

a[1:4] 

para a = array ([1,2,3,4,5,6,7,8,9]), el resultado es

array([2, 3, 4]) 

El primer parámetro de la tupla indica el primer valor que se conservará, y el segundo parámetro indica que no se debe retener el primer valor.

Si utiliza listas en ambos lados, significa que su matriz tiene tantas dimensiones como la longitud de las listas.

Así, con su sintaxis, es probable que necesite algo como esto:

x = A[list1,:,list2] 

dependiendo de la forma de A.

espero que le hayan ayudado.

Answer 2

El método que ya se está utilizando,

A[list1, :][:, list2] 

parece ser la forma más rápida para seleccionar los valores deseados de una matriz de repuestos. Vea a continuación un punto de referencia.

Sin embargo, para responder a su pregunta sobre cómo seleccionar los valores de las filas arbitrarias y columnas de Acon un solo índice, que tendría que utilizar los llamados "advanced indexing":

A[np.array(list1)[:,np.newaxis], np.array(list2)] 

Con la indexación avanzada Si arr1 y arr2 son NDarrays, el componente de (i,j)A[arr1, arr2] es igual

A[arr1[i,j], arr2[i,j]] 

Así que usted quiere arr1[i,j] a ser igual para todos list1[i]j, y arr2[i,j] a ser igual para todos list2[j]i.

que se pueden organizar con la ayuda de broadcasting (véase más adelante) mediante el establecimiento de arr1 = np.array(list1)[:,np.newaxis], y arr2 = np.array(list2).

La forma de arr1 es (len(list1), 1) mientras que la forma de arr2 es (len(list2),) que emite a (1, len(list2)) ya que los nuevos ejes se añaden a la izquierda automáticamente cuando sea necesario.

Cada matriz se puede retransmitir a la forma (len(list1),len(list2)). Esto es exactamente lo que queremos para A[arr1[i,j],arr2[i,j]] para tener sentido, ya que queremos (i,j) para ejecutar todos los índices posibles para un conjunto de resultados de la forma (len(list1),len(list2)).

---

Este es un caso de microanálisis por una prueba de lo que sugiere que A[list1, :][:, list2] es la opción más rápida:

In [32]: %timeit orig(A, list1, list2) 
10 loops, best of 3: 110 ms per loop 

In [34]: %timeit using_listener(A, list1, list2) 
1 loop, best of 3: 1.29 s per loop 

In [33]: %timeit using_advanced_indexing(A, list1, list2) 
1 loop, best of 3: 1.8 s per loop 

---

Esta es la configuración que utilicé para el punto de referencia:

import numpy as np 
import scipy.sparse as sparse 
import random 
random.seed(1) 

def setup(N): 
    A = sparse.rand(N, N, .1, format='lil') 
    list1 = np.random.choice(N, size=N//10, replace=False).tolist() 
    list2 = np.random.choice(N, size=N//20, replace=False).tolist() 
    return A, list1, list2 

def orig(A, list1, list2): 
    return A[list1, :][:, list2] 

def using_advanced_indexing(A, list1, list2): 
    B = A.tocsc() # or `.tocsr()` 
    B = B[np.array(list1)[:, np.newaxis], np.array(list2)] 
    return B 

def using_listener(A, list1, list2): 
    """https://stackoverflow.com/a/26592783/190597 (listener)""" 
    B = A.tocsr()[list1, :].tocsc()[:, list2] 
    return B 

N = 10000 
A, list1, list2 = setup(N) 
B = orig(A, list1, list2) 
C = using_advanced_indexing(A, list1, list2) 
D = using_listener(A, list1, list2) 
assert np.allclose(B.toarray(), C.toarray()) 
assert np.allclose(B.toarray(), D.toarray()) 

Answer 3

para mí la solución de unutbu funciona bien, pero es lenta.

he encontrado como una alternativa rápida,

A = B.tocsr()[np.array(list1),:].tocsc()[:,np.array(list2)] 

Se puede ver que la fila de col y de cortarse por separado, pero cada uno de ellos convertido en el formato más rápido escasa, para obtener el índice este momento.

En mi entorno de prueba, este código es 1000 veces más rápido que el otro.

Espero no decir algo incorrecto o cometer un error.

Answer 4

La indexación simultánea como en B[arr1, arr2] funciona y es más rápido que listener's solution en mi máquina. Consulte In [5] en el ejemplo de Jupyter a continuación. Para compararlo con la respuesta mencionada, consulte In [6]. Además, mi solución no necesita la conversión .tocsc(), por lo que es más fácil de leer IMO.

Tenga en cuenta que para que B[arr1, arr2] funcione, arr1 y arr2 deben ser broadcastable matrices numpy.

A solución mucho más rápida, sin embargo, está utilizando B[list1][:, list2] como pointed out by unutbu. Ver In [7] a continuación.

In [1]: from scipy import sparse 
     : import numpy as np 
     : 
     : 

In [2]: B = sparse.rand(1000, 1000, .1, format='lil') 
     : list1=[1,4,6,8] 
     : list2=[2,4] 
     : 
     : 

In [3]: arr1 = np.array(list1)[:, None] # make arr1 a (n x 1)-array 
     : arr1 
     : 
     : 
Out[3]: 
array([[1], 
     [4], 
     [6], 
     [8]]) 

In [4]: arr2 = np.array(list2)[None, :] # make arr2 a (1 x m)-array 
     : arr2 
     : 
     : 
Out[4]: array([[2, 4]]) 

In [5]: %timeit A = B.tocsr()[arr1, arr2] 
100 loops, best of 3: 13.1 ms per loop 

In [6]: %timeit A = B.tocsr()[np.array(list1),:].tocsc()[:,np.array(list2)] 
100 loops, best of 3: 14.6 ms per loop 

In [7]: %timeit B[list1][:, list2] 
1000 loops, best of 3: 205 µs per loop