mnożenie macierzy w rozmiarze trójkątnym/rzadkim?

Question

Pracuję z bardzo dużym problemem mnożenia macierzy rzadkich (matmul). Jako przykład załóżmy:

• A jest macierzą binarną (75 x 200 000). Jest rzadki, więc używam csc do przechowywania. Trzeba wykonać następujące matmul operacji:

• B = A.transpose() * a

• Wyjście będzie rzadki i symetrycznych macierzy o rozmiarze 200Kx200K.

Niestety, B będzie sposób do dużych przechowywać w pamięci RAM (lub „w rdzeniu”) na moim laptopie. Z drugiej strony, mam szczęście, ponieważ istnieją pewne właściwości B, które powinny rozwiązać ten problem.

Ponieważ B będzie symetryczne wzdłuż przekątnej i rzadkie, mógłbym użyć trójkątnej matrycy (górna/dolna) do przechowywania wyników operacji matmul, a rzadki format przechowywania matrycy może dodatkowo zmniejszyć rozmiar.

Moje pytanie brzmi ... potrafię zmylić lub scipy powiedzieć z wyprzedzeniem, jakie będą wymagania dotyczące pamięci wyjściowej, aby móc wybrać rozwiązanie pamięci masowej za pomocą numpy i uniknąć "matrycy jest zbyt duża" błąd runtime po kilku minutach (godzinach) obliczeń?

Czy można zwiększyć wymagania dotyczące pamięci masowej macierzy poprzez analizę zawartości dwóch macierzy wejściowych przy użyciu algorytmu obliczania przybliżonego?

• https://en.wikipedia.org/wiki/Approximate_counting_algorithm

Jeśli nie, szukam do brutalnej siły rozwiązania. Coś z udziałem Mapa/zmniejszenia, out-of-core przechowywania lub roztworem matmul podziałowych (algorytm strassens) z poniższych linków:

Kilka Mapa/Reduce rozwiązań podziałowe problemowych

• http://www.norstad.org/matrix-multiply/index.html
• http://bpgergo.blogspot.com/2011/08/matrix-multiplication-in-python.html

poza wewnątrzrdzeniowego roztwór (PyTables) składowanie

• Very large matrices using Python and NumPy

matmul rozwiązanie Dzielnica:

• https://en.wikipedia.org/wiki/Strassen_algorithm
• http://facultyfp.salisbury.edu/taanastasio/COSC490/Fall03/Lectures/FoxMM/example.pdf
• http://eli.thegreenplace.net/2012/01/16/python-parallelizing-cpu-bound-tasks-with-multiprocessing/

Dzięki z góry za wszelkie zalecenia, komentarze lub wskazówek!

Answer 1

Ponieważ jesteście po produkcie macierzy z transpozycją, wartość [m, n] jest w zasadzie produktem iloczynowym kolumn m i n w oryginalnej macierzy.

mam zamiar użyć następującego macierz jako przykład zabawek

a = np.array([[0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1], 
       [0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0], 
       [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1]]) 
>>> np.dot(a.T, a) 
array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 
     [0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0], 
     [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1], 
     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 
     [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1], 
     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 
     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 
     [0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0], 
     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 
     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 
     [0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0], 
     [0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 2]]) 

Jest kształtu (3, 12) i ma 7 niezerowe wpisy. Produkt jego transpozycji ma oczywiście kształt (12, 12) i ma 16 wpisów niezerowych, z których 6 ma przekątną, więc wymaga tylko przechowywania 11 elementów.

Można uzyskać dobry pomysł, co wielkość swojej macierzy wyjściowej będzie w jeden z dwóch sposobów:

FORMAT CSR

Jeśli oryginalna matryca ma C niezerowe kolumny , twoja nowa matryca będzie miała co najwyżej C**2 niezerowych pozycji, z których C ma przekątną i jest pewna, że ​​nie ma zera, a pozostałych wpisów wystarczy zachować połowę, więc jest to najwyżej (C**2 + C)/2 zero elementów. Oczywiście wiele z nich będzie również zerowych, więc jest to prawdopodobnie zawyżona wartość.

Jeżeli matryca jest przechowywana w formie csr, wówczas cecha odpowiedniego scipy obiektu indices ma tablicę wskaźników kolumny wszystkich nie zerowych elementów, dzięki czemu można łatwo obliczyć wyżej oszacowania jako:

>>> a_csr = scipy.sparse.csr_matrix(a) 
>>> a_csr.indices 
array([ 2, 11, 1, 7, 10, 4, 11]) 
>>> np.unique(a_csr.indices).shape[0] 
6 

więc istnieje 6 kolumn niezerowych wpisów, a więc oszacowanie byłoby co najwyżej 36 niezerowych wpisów, o wiele bardziej niż realnej 16.

CSC FORMAT

Jeśli zamiast indeksów kolumn niezerowych elementów mamy indeksy wierszy, możemy faktycznie lepiej oszacować. Aby iloczyn punktowy dwóch kolumn był niezerowy, muszą one mieć niezerowy element w tym samym wierszu. Jeśli w danym wierszu występują niezerowe elementy R, przyczynią się one do umieszczenia niezerowych elementów R**2 w produkcie. Kiedy sumujesz to dla wszystkich wierszy, jesteś zobowiązany liczyć kilka elementów więcej niż jeden raz, więc jest to również górna granica.

indeksów rzędzie niezerowych elementów swojej matrycy znajdują się w atrybucie rozrzedzony matrycy csc indices, więc ta ocena może być obliczona w następujący sposób:

>>> a_csc = scipy.sparse.csc_matrix(a) 
>>> a_csc.indices 
array([1, 0, 2, 1, 1, 0, 2]) 
>>> rows, where = np.unique(a_csc.indices, return_inverse=True) 
>>> where = np.bincount(where) 
>>> rows 
array([0, 1, 2]) 
>>> where 
array([2, 3, 2]) 
>>> np.sum(where**2) 
17 

To jest cholernie blisko prawdziwy 16! I nie jest to rzeczywiście przypadek, że szacunek ten jest właściwie taka sama, jak:

>>> np.sum(np.dot(a.T,a),axis=None) 
17 

w każdym przypadku, następujący kod powinien pozwalają zobaczyć, że ocena jest bardzo dobra:

def estimate(a) : 
    a_csc = scipy.sparse.csc_matrix(a) 
    _, where = np.unique(a_csc.indices, return_inverse=True) 
    where = np.bincount(where) 
    return np.sum(where**2) 

def test(shape=(10,1000), count=100) : 
    a = np.zeros(np.prod(shape), dtype=int) 
    a[np.random.randint(np.prod(shape), size=count)] = 1 
    print 'a non-zero = {0}'.format(np.sum(a)) 
    a = a.reshape(shape) 
    print 'a.T * a non-zero = {0}'.format(np.flatnonzero(np.dot(a.T, 
                   a)).shape[0]) 
    print 'csc estimate = {0}'.format(estimate(a)) 

>>> test(count=100) 
a non-zero = 100 
a.T * a non-zero = 1065 
csc estimate = 1072 
>>> test(count=200) 
a non-zero = 199 
a.T * a non-zero = 4056 
csc estimate = 4079 
>>> test(count=50) 
a non-zero = 50 
a.T * a non-zero = 293 
csc estimate = 294