Scipy Sparse Matrix - Mnożenie wektorów gęstych - Blocks vs Large Matrix

Question

Mam wiele rzadkich matryc scipy (obecnie w formacie CSR), które muszę pomnożyć za pomocą gęstego, numpy wektora 1D. wektor nazywa G:

print G.shape, G.dtype 
(2097152,) complex64 

Każda macierz rzadką ma kształt (16384,2097152) i jest bardzo rzadki. Gęstość wynosi około 4,0-6. Mam listę 100 tych rzadkich macierzy o nazwie spmats.

mogę łatwo mnożenie każdej matrycy z G tak:

res = [spmat.dot(G) for spmat in spmats] 

Powoduje listy gęstych wektory kształtu (16384,) prawidłowo.

My stosowanie jest dość perfomance krytyczne, a więc próbuje się naprzemiennie, co jest najpierw łączyć wszystkie rzadkich matryc w jednej dużej sparce matrycy, a następnie za pomocą tylko jednego połączenia z dot() tak:

import scipy.sparse as sp 
SPMAT = sp.vstack(spmats, format='csr') 
RES = SPMAT.dot(G) 

Powoduje to, że jeden długi wektor RES ma kształt (1638400,) i jest połączoną wersją wszystkich wektorów wynikowych w powyższej res, zgodnie z oczekiwaniami. Sprawdziłem, czy wyniki są identyczne.

Może jestem całkowicie w błędzie, ale spodziewałem się, że drugi przypadek powinien być szybszy od pierwszego, ponieważ istnieje znacznie mniej numpy, alokacja pamięci, tworzenie obiektów Pythona, pętle Pythona, itp. dbają o czas wymagany do połączenia rzadkich macierzy, a jedynie czas na obliczenie wyniku. Według %timeit jednak:

%timeit res = [spmat.dot(G) for spmat in spmats] 
10 loops, best of 3: 91.5 ms per loop 
%timeit RES = SPMAT.dot(G) 
1 loops, best of 3: 389 ms per loop 

Sprawdziłem, że nie uciekam z pamięci w każdej pracy, i nic podejrzany wydaje się dzieje. Czy jestem szalony, czy to naprawdę dziwne? Czy oznacza to, że wszystkie rzadkie produkty matrycowo-wektorowe powinny być wykonywane w blokach, kilka wierszy na raz, aby były szybsze? O ile rozumiem, rzadki czas mnożenia macierzy gęstym wektorem powinien być liniowy w liczbie niezerowych elementów, który nie zmienia się w dwóch powyższych przypadkach. Co może powodować taką różnicę?

biegnę na pojedynczym rdzeniem maszynie Linux z 4GB RAM, używając EPD7.3

EDIT:

Oto mały przykład, który reprodukuje problem dla mnie:

import scipy.sparse as sp 
import numpy as n 

G = n.random.rand(128**3) + 1.0j*n.random.rand(128**3) 

spmats = [sp.rand (128**2, 128**3, density = 4e-6, format = 'csr', dtype=float64) for i in range(100)] 
SPMAT = sp.vstack(spmats, format='csr') 

%timeit res = [spmat.dot(G) for spmat in spmats] 
%timeit RES = SPMAT.dot(G) 

Otrzymuję:

1 loops, best of 3: 704 ms per loop 
1 loops, best of 3: 1.34 s per loop 

Różnica wydajności w tym przypadku nie jest tak duża jak z moimi własnymi rzadkimi macierzami, które mają pewną strukturę (prawdopodobnie z powodu buforowania), ale jeszcze gorsze jest łączenie macierzy.

Próbowałem z obu scipy 10.1 i 12.0.

Answer 1

Nie znalazłem powodu do dziwnego zachowania wspomnianego w pytaniu, jednak znalazłem sposób, aby znacznie przyspieszyć moje obliczenia, które mogą być przydatne dla innych osób.

Ponieważ w moim konkretnym przypadku używam obliczeń rzadkiej macierzy float32 i gęstego wektora complex64, mogę oddzielnie mnożyć rzeczywiste i urojone elementy. Zapewnia to przyspieszenie 4x.

Dzieje 2.35s z SPMAT.shape == (16384000, 2097152):

RES = SPMAT.dot(G) 

Mimo to zajmuje tylko 541ms:

RES = n.zeros((SPMAT.shape[0],),dtype=complex64) 
RES.real = SPMAT.dot(G.real); RES.imag = SPMAT.dot(G.imag) 

a wynik jest identyczny. Myślę, że być może wstępna alokacja n.zeros może nie być konieczna, ale nie wiem jak inaczej to zrobić.