Accumulo efficiente di una raccolta di matrici scarse spettrali

Question

Ho una raccolta di O (N) NxN scipy.sparse.csr_matrix e ciascuna matrice sparsa è nell'ordine di N elementi impostati. Voglio aggiungere tutte queste matrici insieme per ottenere un normale array numpy NxN. (N è dell'ordine di 1000). La disposizione di elementi diversi da zero all'interno delle matrici è tale che la somma risultante non è certamente sparsa (praticamente nessun elemento a zero rimane in effetti).

Al momento mi sto solo facendo

reduce(lambda x,y: x+y,[m.toarray() for m in my_sparse_matrices]) 

che funziona, ma è un po 'lento: naturalmente l'enorme quantità di elaborazione inutile di zeri che sta andando in là è assolutamente orribile.

C'è un modo migliore? Non c'è niente di ovvio per me nello docs.

Aggiornamento: come suggerito dall'utente545424, ho provato lo schema alternativo di sommare le matrici sparse e anche di sommare matrici sparse su una matrice densa. Il codice seguente mostra tutti gli approcci per l'esecuzione in tempi paragonabili (Python 2.6.6 su amd64 Debian/Squeeze, su un i7 quad-core)

import numpy as np 
import numpy.random 
import scipy 
import scipy.sparse 
import time 

N=768 
S=768 
D=3 

def mkrandomsparse(): 
    m=np.zeros((S,S),dtype=np.float32) 
    r=np.random.random_integers(0,S-1,D*S) 
    c=np.random.random_integers(0,S-1,D*S) 
    for e in zip(r,c): 
     m[e[0],e[1]]=1.0 
    return scipy.sparse.csr_matrix(m) 

M=[mkrandomsparse() for i in xrange(N)] 

def plus_dense(): 
    return reduce(lambda x,y: x+y,[m.toarray() for m in M]) 

def plus_sparse(): 
    return reduce(lambda x,y: x+y,M).toarray() 

def sum_dense(): 
    return sum([m.toarray() for m in M]) 

def sum_sparse(): 
    return sum(M[1:],M[0]).toarray() 

def sum_combo(): # Sum the sparse matrices 'onto' a dense matrix? 
    return sum(M,np.zeros((S,S),dtype=np.float32)) 

def benchmark(fn): 
    t0=time.time() 
    fn() 
    t1=time.time() 
    print "{0:16}: {1:.3f}s".format(fn.__name__,t1-t0) 

for i in xrange(4): 
    benchmark(plus_dense) 
    benchmark(plus_sparse) 
    benchmark(sum_dense) 
    benchmark(sum_sparse) 
    benchmark(sum_combo) 
    print 

e disconnette

plus_dense  : 1.368s 
plus_sparse  : 1.405s 
sum_dense  : 1.368s 
sum_sparse  : 1.406s 
sum_combo  : 1.039s 

anche se è possibile ottenere uno approccio o l'altro per uscire avanti di un fattore di circa 2 facendo confusione con i parametri N, S, D ... ma nulla di simile all'ordine di grandezza che si spera di vedere dal considerare il numero di zero aggiunge che dovrebbe essere possibile saltare.

Answer 1

Penso di aver trovato un modo per accelerarlo di un fattore di ~ 10 se le vostre matrici sono molto rare.

In [1]: from scipy.sparse import csr_matrix 

In [2]: def sum_sparse(m): 
    ...:  x = np.zeros(m[0].shape) 
    ...:  for a in m: 
    ...:   ri = np.repeat(np.arange(a.shape[0]),np.diff(a.indptr)) 
    ...:   x[ri,a.indices] += a.data 
    ...:  return x 
    ...: 

In [6]: m = [np.zeros((100,100)) for i in range(1000)] 

In [7]: for x in m: 
    ...:  x.ravel()[np.random.randint(0,x.size,10)] = 1.0 
    ...:  

     m = [csr_matrix(x) for x in m] 

In [17]: (sum(m[1:],m[0]).todense() == sum_sparse(m)).all() 
Out[17]: True 

In [18]: %timeit sum(m[1:],m[0]).todense() 
10 loops, best of 3: 145 ms per loop 

In [19]: %timeit sum_sparse(m) 
100 loops, best of 3: 18.5 ms per loop 

Answer 2

Non puoi semplicemente aggiungerli prima di convertirli in una matrice densa?

>>> sum(my_sparse_matrices[1:],my_sparse_matrices[0]).todense() 

Answer 3

Questa non è una risposta completa (e anch'io vorrei vedere una risposta più dettagliata), ma è possibile ottenere un fattore semplice di due o più di miglioramento, non creando risultati intermedi:

def sum_dense(): 
    return sum([m.toarray() for m in M]) 

def sum_dense2(): 
    return sum((m.toarray() for m in M)) 

Sulla mia macchina (YMMV), questo risulta essere il calcolo più veloce. Inserendo la sommatoria in un () anziché in un [], prima della sommatoria costruiamo un generatore anziché l'intero elenco.

Answer 4

@ user545424 ha già pubblicato quella che probabilmente sarà la soluzione più veloce. Qualcosa con lo stesso spirito che è più leggibile e ~ stessa velocità ... nonzero() ha tutti i tipi di applicazioni utili.