Efektywnie gromadzi kolekcję rzadkich macierzy scipy

Question

Mam zbiór O (N) NxN scipy.sparse.csr_matrix, a każda rzadka macierz ma zestaw N elementów. Chcę dodać wszystkie te macierze, aby uzyskać regularną tablicę NxN. (N jest rzędu 1000). Rozmieszczenie niezerowych elementów w macierzach jest takie, że wynikowa suma z pewnością nie jest rzadka (właściwie nie ma już zerowych elementów).

W tej chwili jestem po prostu robi

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

który działa, ale jest nieco powolny: oczywiście ogromna ilość bezcelowym przetwarzania zer, które się tam dzieje jest absolutnie przerażające.

Czy istnieje lepszy sposób? Nie ma nic oczywistego w docs.

Aktualizacja: zgodnie z sugestią użytkownika 545424, wypróbowałem alternatywny schemat sumowania rzadkich macierzy, a także sumowania rzadkich macierzy na gęstej macierzy. Poniższy kod pokazuje wszystkie podejścia do uruchomienia w porównywalnym czasie (Python 2.6.6 na amd64 Debian/Squeeze na quad-core i7)

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 

i wylogowuje

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

choć można dostać jeden podejście lub inne, aby wyjść z wyprzedzeniem przez współczynnik 2 lub więcej, bawiąc się z parametrami N, S, D ... ale nic takiego jak poprawa wielkości magii, którą byś oczekiwał patrząc na liczbę zerową, powinno to być możliwym do pominięcia.

Answer 1

Myślę, że znalazłem sposób, aby przyspieszyć to o współczynnik ~ 10, jeśli macierze są bardzo rzadkie.

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

Nie możesz po prostu dodać ich przed przekształceniem w gęstą matrycę?

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

Answer 3

To nie jest pełna odpowiedź (i ja również chciałbym zobaczyć bardziej szczegółowe odpowiedzi), ale można uzyskać łatwy czynnik dwóch lub większej poprawy, nie tworząc wyników pośrednich:

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

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

Na moim komputerze (YMMV) powoduje to najszybsze obliczenia. Poprzez umieszczenie podsumowania w postaci () zamiast [], konstruujemy raczej generator niż całą listę przed podsumowaniem.

Answer 4

@ użytkownik545424 opublikował już, co prawdopodobnie będzie najszybszym rozwiązaniem. Coś w tym samym duchu, który jest bardziej czytelny i ma taką samą prędkość ... nonzero() ma wiele przydatnych aplikacji.

def sum_sparse(m): 
     x = np.zeros(m[0].shape,m[0].dtype) 
     for a in m: 
      # old lines 
      #ri = np.repeat(np.arange(a.shape[0]),np.diff(a.indptr)) 
      #x[ri,a.indices] += a.data 
      # new line 
      x[a.nonzero()] += a.data 
     return x