Aggiunta di una matrice molto ripetitiva a una sparsa in numpy/scipy?

Question

Sto cercando di implementare una funzione in NumPy/Scipy per calcolare Jensen-Shannon divergence tra un singolo vettore (di allenamento) e un gran numero di altri vettori (di osservazione). I vettori di osservazione sono memorizzati in un numero molto grande (500.000x65536) Scipy sparse matrix (una matrice densa non si adatta alla memoria).

Come parte dell'algoritmo ho bisogno di calcolare T + O _i per ogni osservazione vettore O _i, dove T è il vettore di formazione. Non ero in grado di trovare un modo per farlo usando le solite regole di broadcasting di NumPy, dal momento che le matrici sparse non sembrano supportare quelle (se T è lasciato come un array denso, Scipy prova a rendere prima la matrice sparsa densa, che esegue memoria insufficiente, se faccio T in una matrice sparsa, T + O _i non riesce perché le forme sono incoerenti).

Attualmente mi sto prendendo il passo grossolanamente inefficiente delle piastrelle del vettore di formazione in una 500,000x65536 matrice sparsa:

training = sp.csr_matrix(training.astype(np.float32)) 
tindptr = np.arange(0, len(training.indices)*observations.shape[0]+1, len(training.indices), dtype=np.int32) 
tindices = np.tile(training.indices, observations.shape[0]) 
tdata = np.tile(training.data, observations.shape[0]) 
mtraining = sp.csr_matrix((tdata, tindices, tindptr), shape=observations.shape) 

Ma questo richiede una grande quantità di memoria (circa 6 GB), quando è solo la memorizzazione ~ 1500 elementi "reali". È anche piuttosto lento da costruire.

Ho cercato di essere intelligente utilizzando stride_tricks per rendere i membri dell'indptr e dei dati della matrice CSR non utilizzare memoria aggiuntiva sui dati ripetuti.

training = sp.csr_matrix(training) 
mtraining = sp.csr_matrix(observations.shape,dtype=np.int32) 
tdata = training.data 
vdata = np.lib.stride_tricks.as_strided(tdata, (mtraining.shape[0], tdata.size), (0, tdata.itemsize)) 
indices = training.indices 
vindices = np.lib.stride_tricks.as_strided(indices, (mtraining.shape[0], indices.size), (0, indices.itemsize)) 
mtraining.indptr = np.arange(0, len(indices)*mtraining.shape[0]+1, len(indices), dtype=np.int32) 
mtraining.data = vdata 
mtraining.indices = vindices 

Ma questo non funziona perché i panorami strided mtraining.data e mtraining.indices sono la forma sbagliata (e secondo this answer non c'è alcun modo per rendere la forma a destra). Cercare di renderli flat con l'iteratore .flat fallisce perché non sembra abbastanza come un array (ad esempio non ha un membro dtype) e l'uso del metodo flatten() finisce per fare una copia.

C'è un modo per ottenere questo risultato?

Answer 1

L'altra opzione, che non avevo nemmeno preso in considerazione, è implementare la somma nel formato sparse da soli, in modo da poter sfruttare appieno la natura periodica del proprio array. Questo può essere molto facile da fare, se abusare di questo particolare comportamento di matrici sparse di SciPy:

>>> a = sps.csr_matrix([1,2,3,4]) 
>>> a.data 
array([1, 2, 3, 4]) 
>>> a.indices 
array([0, 1, 2, 3]) 
>>> a.indptr 
array([0, 4]) 

>>> b = sps.csr_matrix((np.array([1, 2, 3, 4, 5]), 
...      np.array([0, 1, 2, 3, 0]), 
...      np.array([0, 5])), shape=(1, 4)) 
>>> b 
<1x4 sparse matrix of type '<type 'numpy.int32'>' 
    with 5 stored elements in Compressed Sparse Row format> 
>>> b.todense() 
matrix([[6, 2, 3, 4]]) 

Quindi non hanno nemmeno bisogno di cercare le coincidenze tra il vettore di formazione e di ciascuna delle righe della matrice di osservazione per aggiungerli: basta racchiudere tutti i dati con i puntatori giusti lì, e ciò che deve essere sommato, verrà sommato quando si accede ai dati.

EDIT

Data la lentezza del primo codice, è possibile barattare la memoria per la velocità come segue:

def csr_add_sparse_vec(sps_mat, sps_vec) : 
    """Adds a sparse vector to every row of a sparse matrix""" 
    # No checks done, but both arguments should be sparse matrices in CSR 
    # format, both should have the same number of columns, and the vector 
    # should be a vector and have only one row. 

    rows, cols = sps_mat.shape 
    nnz_vec = len(sps_vec.data) 
    nnz_per_row = np.diff(sps_mat.indptr) 
    longest_row = np.max(nnz_per_row) 

    old_data = np.zeros((rows * longest_row,), dtype=sps_mat.data.dtype) 
    old_cols = np.zeros((rows * longest_row,), dtype=sps_mat.indices.dtype) 

    data_idx = np.arange(longest_row) < nnz_per_row[:, None] 
    data_idx = data_idx.reshape(-1) 
    old_data[data_idx] = sps_mat.data 
    old_cols[data_idx] = sps_mat.indices 
    old_data = old_data.reshape(rows, -1) 
    old_cols = old_cols.reshape(rows, -1) 

    new_data = np.zeros((rows, longest_row + nnz_vec,), 
         dtype=sps_mat.data.dtype) 
    new_data[:, :longest_row] = old_data 
    del old_data 
    new_cols = np.zeros((rows, longest_row + nnz_vec,), 
         dtype=sps_mat.indices.dtype) 
    new_cols[:, :longest_row] = old_cols 
    del old_cols 
    new_data[:, longest_row:] = sps_vec.data 
    new_cols[:, longest_row:] = sps_vec.indices 
    new_data = new_data.reshape(-1) 
    new_cols = new_cols.reshape(-1) 
    new_pointer = np.arange(0, (rows + 1) * (longest_row + nnz_vec), 
          longest_row + nnz_vec) 

    ret = sps.csr_matrix((new_data, new_cols, new_pointer), 
         shape=sps_mat.shape) 
    ret.eliminate_zeros() 

    return ret 

Non è veloce come prima, ma è possibile farlo 10.000 righe in circa 1 s.:

In [2]: a 
Out[2]: 
<10000x65536 sparse matrix of type '<type 'numpy.float64'>' 
    with 15000000 stored elements in Compressed Sparse Row format> 

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

In [4]: csr_add_sparse_vec(a, b) 
Out[4]: 
<10000x65536 sparse matrix of type '<type 'numpy.float64'>' 
    with 30000000 stored elements in Compressed Sparse Row format> 

In [5]: %timeit csr_add_sparse_vec(a, b) 
1 loops, best of 3: 956 ms per loop 

---

EDIT Questo codice è molto, molto lento

def csr_add_sparse_vec(sps_mat, sps_vec) : 
    """Adds a sparse vector to every row of a sparse matrix""" 
    # No checks done, but both arguments should be sparse matrices in CSR 
    # format, both should have the same number of columns, and the vector 
    # should be a vector and have only one row. 

    rows, cols = sps_mat.shape 

    new_data = sps_mat.data 
    new_pointer = sps_mat.indptr.copy() 
    new_cols = sps_mat.indices 

    aux_idx = np.arange(rows + 1) 

    for value, col in itertools.izip(sps_vec.data, sps_vec.indices) : 
     new_data = np.insert(new_data, new_pointer[1:], [value] * rows) 
     new_cols = np.insert(new_cols, new_pointer[1:], [col] * rows) 
     new_pointer += aux_idx 

    return sps.csr_matrix((new_data, new_cols, new_pointer), 
          shape=sps_mat.shape)