Ersatz für numpy Rundfunk mit scipy.sparse.csc_matrix

Question

ich den folgenden Ausdruck in meinem Code haben:

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

wo b eine ndarray der Form ist (M, N) und x ist ein ndarray der Form (M,). Nun, b ist eigentlich spärlich, so für Speichereffizienz möchte ich in einer scipy.sparse.csc_matrix oder csr_matrix ersetzen. Das Senden auf diese Weise wird jedoch nicht implementiert (obwohl Division oder Multiplikation garantiert Sparsity beibehalten) (die Einträge x sind nicht Null) und löst eine NotImplementedError. Gibt es eine sparse Funktion, die mir nicht bewusst ist, würde das tun, was ich will? (dot() würde entlang der falschen Achse summieren.)

Answer 1

Wenn b in CSC-Format ist, dann hat b.data die Nicht-Null-Einträge von b und b.indices hat den Zeilenindex jedes der Nicht-Null-Einträge, so dass Sie Ihre Abteilung tun können, als :

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

es ist hackier als elegante Lösung Warren, aber es wird in den meisten Einstellungen wahrscheinlich auch schneller sein:

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 

Sie können die Zeit fast die Hälfte in dem obigen Beispiel schneiden, wenn Sie die Aufteilung in-place tun, d.h .:

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

Zum Implementieren a = (b/x[:, np.newaxis]).sum(axis=1) können Sie a = b.sum(axis=1).A1/x verwenden. Das A1 Attribut gibt das 1D ndarray zurück, also ist das Ergebnis ein 1D ndarray, kein matrix. Dieser kurze Ausdruck funktioniert, weil Sie beide Skalierung von xund Summieren entlang der Achse sind 1. Zum Beispiel:

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]) 

Allgemeiner gesagt, wenn Sie die Zeilen einer Sparse Matrix mit einem Vektor x, skalieren möchten könnten Sie multiplizieren Sie b auf der linken Seite mit einer dünnen Matrix, die 1.0/x auf der Diagonale enthält. Mit der Funktion scipy.sparse.spdiags kann eine solche Matrix erstellt werden. Zum Beispiel:

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]])