Zastępstwo dla radiofonii numpy wykorzystaniem scipy.sparse.csc_matrix

Question

mam w kodzie następujące wyrażenie:

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

gdzie b jest ndarray kształtu (M, N) i x jest ndarray kształtu (M,). Teraz, b jest rzeczywiście rzadki, więc dla wydajności pamięci chciałbym zastąpić w scipy.sparse.csc_matrix lub csr_matrix. Nadawanie w ten sposób nie jest jednak realizowane (nawet jeśli gwarantowane jest dzielenie lub mnożenie w celu zachowania niespójności) (pozycje x są niezerowe) i zwiększa liczbę NotImplementedError. Czy istnieje funkcja sparse Nie jestem świadoma, że ​​zrobiłbym to, co chcę? (dot() sumuje się wzdłuż niewłaściwej osi.)

Answer 1

Jeśli b jest w formacie CSC, a następnie b.data ma niezerowe wpisy b i b.indices ma indeks wiersza każdego z niezerowych wpisów, więc można zrobić podział jako :

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

to hackier niż eleganckie rozwiązanie Warrena, ale prawdopodobnie będzie to również szybszy w większości ustawień:

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 

Możesz przeciąć czas prawie o połowę w powyższym przykładzie, jeśli wykonasz podział w miejscu, tj .:

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

Aby zaimplementować a = (b/x[:, np.newaxis]).sum(axis=1), można użyć a = b.sum(axis=1).A1/x. Atrybut A1 zwraca 1d ndarray, więc wynikiem jest 1d ndarray, a nie matrix. To zwięzłe wyrażenie działa, ponieważ jesteś zarówno skalowanie przez xi zsumowanie wzdłuż osi 1. Na przykład:

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

Ogólniej, jeśli chcesz przeskalować wiersze rozrzedzony matrycy z wektorem x, można pomnóż b po lewej stronie z rzadką macierzą zawierającą 1.0/x na przekątnej. Funkcja scipy.sparse.spdiags może być użyta do utworzenia takiej matrycy. Na przykład:

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