sostituto per la trasmissione NumPy utilizzando scipy.sparse.csc_matrix

Question

ho nel mio codice seguente espressione:

a = (b/x[:, np.newaxis]).sum(axis=1) 

dove b è un ndarray di forma (M, N), ed è un x ndarray di forma (M,). Ora, b è in realtà scarso, quindi per l'efficienza della memoria vorrei sostituire in un scipy.sparse.csc_matrix o csr_matrix. Tuttavia, la trasmissione in questo modo non viene implementata (anche se la divisione o la moltiplicazione è garantita per mantenere la scarsità) (le voci di x sono diverse da zero) e genera uno NotImplementedError. Esiste una funzione sparse Non sono a conoscenza di ciò che farebbe ciò che voglio? (dot() sommerebbe lungo l'asse sbagliato.)

Answer 1

Se b è in formato CSC, poi b.data ha le voci non-zero di b, e b.indices ha l'indice di riga di ciascuna delle voci non-zero, in modo da poter fare la vostra divisione come :

b.data /= np.take(x, b.indices) 

E 'hackier di soluzione elegante di Warren, ma probabilmente anche più veloce nella maggior parte delle impostazioni:

b = sps.rand(1000, 1000, density=0.01, format='csc') 
x = np.random.rand(1000) 

def row_divide_col_reduce(b, x): 
    data = b.data.copy()/np.take(x, b.indices) 
    ret = sps.csc_matrix((data, b.indices.copy(), b.indptr.copy()), 
         shape=b.shape) 
    return ret.sum(axis=1) 

def row_divide_col_reduce_bis(b, x): 
    d = sps.spdiags(1.0/x, 0, len(x), len(x)) 
    return (d * b).sum(axis=1) 

In [2]: %timeit row_divide_col_reduce(b, x) 
1000 loops, best of 3: 210 us per loop 

In [3]: %timeit row_divide_col_reduce_bis(b, x) 
1000 loops, best of 3: 697 us per loop 

In [4]: np.allclose(row_divide_col_reduce(b, x), 
    ...:    row_divide_col_reduce_bis(b, x)) 
Out[4]: True 

È possibile ridurre il tempo quasi a metà nell'esempio di cui sopra se si fa la divisione sul posto, cioè .:

def row_divide_col_reduce(b, x): 
    b.data /= np.take(x, b.indices) 
    return b.sum(axis=1) 

In [2]: %timeit row_divide_col_reduce(b, x) 
10000 loops, best of 3: 131 us per loop 

Answer 2

Per implementare a = (b/x[:, np.newaxis]).sum(axis=1), è possibile utilizzare a = b.sum(axis=1).A1/x. L'attributo A1 restituisce l'array 1D n, quindi il risultato è un array 1D nd, non uno matrix. Questa espressione concisa funziona perché si sono entrambi scala da xe sommando lungo l'asse 1. Ad esempio:

In [190]: b 
Out[190]: 
<3x3 sparse matrix of type '<type 'numpy.float64'>' 
     with 5 stored elements in Compressed Sparse Row format> 

In [191]: b.A 
Out[191]: 
array([[ 1., 0., 2.], 
     [ 0., 3., 0.], 
     [ 4., 0., 5.]]) 

In [192]: x 
Out[192]: array([ 2., 3., 4.]) 

In [193]: b.sum(axis=1).A1/x 
Out[193]: array([ 1.5 , 1. , 2.25]) 

Più in generale, se si vuole scalare le righe di una matrice sparsa con un vettore x, si potrebbe moltiplicare b a sinistra con una matrice sparsa contenente 1.0/x sulla diagonale. La funzione scipy.sparse.spdiags può essere utilizzata per creare una tale matrice. Per esempio:

In [71]: from scipy.sparse import csc_matrix, spdiags 

In [72]: b = csc_matrix([[1,0,2],[0,3,0],[4,0,5]], dtype=np.float64) 

In [73]: b.A 
Out[73]: 
array([[ 1., 0., 2.], 
     [ 0., 3., 0.], 
     [ 4., 0., 5.]]) 

In [74]: x = array([2., 3., 4.]) 

In [75]: d = spdiags(1.0/x, 0, len(x), len(x)) 

In [76]: d.A 
Out[76]: 
array([[ 0.5  , 0.  , 0.  ], 
     [ 0.  , 0.33333333, 0.  ], 
     [ 0.  , 0.  , 0.25  ]]) 

In [77]: p = d * b 

In [78]: p.A 
Out[78]: 
array([[ 0.5 , 0. , 1. ], 
     [ 0. , 1. , 0. ], 
     [ 1. , 0. , 1.25]]) 

In [79]: a = p.sum(axis=1) 

In [80]: a 
Out[80]: 
matrix([[ 1.5 ], 
     [ 1. ], 
     [ 2.25]])