slicing sparse (scipy) matrix

Question

Apprezzerei qualsiasi aiuto, per capire il seguente comportamento quando si affetta una lil_matrix (A) dal pacchetto scipy.sparse.

In realtà, vorrei estrarre una sottomatrice basata su un elenco di indici arbitrari per righe e colonne.

Quando ho usato questa due righe di codice:

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

Tutto andava bene e ho potuto estrarre il sottomatrice destra.

Quando ho provato a fare questo in una sola riga, non è riuscito (La matrice di ritorno era vuoto)

x=A[list 1,list 2] 

Perché è così? Nel complesso, ho usato un comando simile in MATLAB e lì funziona. Quindi, perché non usare il primo, dal momento che funziona? Sembra essere abbastanza che richiede tempo. Dal momento che devo passare attraverso una grande quantità di voci, vorrei accelerarlo usando un singolo comando. Forse uso il tipo di matrice sparsa errata ... Qualche idea?

Answer 1

affettatura avviene con la seguente sintassi:

a[1:4] 

per a = array ([1,2,3,4,5,6,7,8,9]), il risultato è

array([2, 3, 4]) 

Il primo parametro della tupla indica il primo valore da conservare e il secondo parametro indica il primo valore che non deve essere mantenuto.

Se si utilizzano elenchi su entrambi i lati, significa che l'array ha tutte le dimensioni della lunghezza degli elenchi.

Così, con la sintassi, si avrà probabilmente bisogno di qualcosa come questo:

x = A[list1,:,list2] 

a seconda della forma di A.

spero che risulta utile.

Answer 2

Il metodo si sta già utilizzando,

A[list1, :][:, list2] 

sembra essere il modo più veloce per selezionare i valori desiderati da una matrice di pezzi di ricambio. Vedi sotto per un punto di riferimento.

Tuttavia, per rispondere alla tua domanda su come selezionare i valori da righe arbitrari e colonne di Acon un unico indice, si avrebbe bisogno di utilizzare la cosiddetta "advanced indexing":

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

Con l'indicizzazione avanzata , se arr1 e arr2 sono NDarrays, la componente (i,j) di A[arr1, arr2] uguale

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

Così si vorrebbe arr1[i,j] eguagliare list1[i] per tutti j, e arr2[i,j] di eguagliare list2[j] per tutti i.

che può essere organizzato con l'aiuto di broadcasting (vedi sotto) modificando arr1 = np.array(list1)[:,np.newaxis] e arr2 = np.array(list2).

La forma arr1 è (len(list1), 1) mentre la forma di arr2 è (len(list2),) che trasmette a (1, len(list2)) poiché nuovi assi vengono aggiunti automaticamente la partita quando necessario.

Ogni array può essere ulteriormente trasmesso per modellare (len(list1),len(list2)). Questo è esattamente ciò che vogliamo per A[arr1[i,j],arr2[i,j]] per avere senso, poiché vogliamo che (i,j) esegua tutti gli indici possibili per un array di risultati di forma (len(list1),len(list2)).

---

Ecco un caso microbenchmark per una prova che suggerisce che A[list1, :][:, list2] è l'opzione più veloce:

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 

---

Qui è la configurazione che ho usato per il benchmark:

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

per me la soluzione di unutbu funziona bene, ma è lenta.

ho trovato come un veloce alternativa,

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

Si può vedere che Row e Col di ottenere tagliato a parte, ma ognuno convertiti nel formato più veloce sparse, per ottenere indice di questo momento.

Nel mio ambiente di test questo codice è 1000 volte più veloce dell'altro.

Spero di non dire qualcosa di sbagliato o di sbagliare.

Answer 4

L'indicizzazione simultanea come in B[arr1, arr2] funziona ed è più veloce di listener's solution sulla mia macchina. Vedere In [5] nell'esempio Jupyter di seguito. Per confrontarlo con la risposta citata, fare riferimento a In [6]. Inoltre, la mia soluzione non ha bisogno della conversione .tocsc(), rendendo IMO più leggibile.

Si prega di notare che per B[arr1, arr2] al lavoro, arr1 e arr2 deve essere broadcastable array numpy.

A soluzione molto più veloce, tuttavia, utilizza B[list1][:, list2] come pointed out by unutbu. Vedi In [7] qui sotto.

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