numpy matrix multiplication zu triangular/sparse storage?

Question

Ich arbeite mit einem sehr großen Sparse Matrix Multiplikation (Matmul) Problem. Als ein Beispiel lassen Sie uns sagen:

• A ist eine binäre (75 x 200.000) Matrix. Es ist spärlich, also benutze ich csc für den Speicher. Ich brauche den folgenden matmul Betrieb zu tun:

• B = A.transpose() * A

• Der Ausgang wird eine spärliche und symmetrische Matrix der Größe 200Kx200K sein.

Leider B wird Weg zu groß sein, im RAM zu speichern (oder "in Kern") auf meinem Laptop. Auf der anderen Seite bin ich glücklich, da es einige Eigenschaften für B gibt, die dieses Problem lösen sollten.

Da B entlang der Diagonalen und spärlich symmetrisch ist, könnte ich eine Dreiecksmatrix (oben/unten) verwenden, um die Ergebnisse der Matmul-Operation zu speichern, und ein Sparse-Matrix-Speicherformat könnte die Größe weiter reduzieren.

Meine Frage ist ... kann numpy oder scipy gesagt werden, vor der Zeit, wie die Ausgabe Speicheranforderungen aussehen werden, so dass ich eine Speicherlösung mit numpy auswählen und die "Matrix ist zu groß" vermeiden Laufzeitfehler nach mehreren Minuten (Stunden) der Berechnung?

Mit anderen Worten, können die Speicheranforderungen für die Matrix multipliziert werden, indem der Inhalt der zwei Eingangsmatrizen unter Verwendung eines approximativen Zählalgorithmus analysiert wird?

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

Wenn nicht, ich bin auf der Suche in eine Brute-Force-Lösung. Etwas Einbeziehung Karte/reduzieren, out-of-Core-Speicher oder eine matmul Unterteilung Lösung (Strassen-Algorithmus) aus den folgenden Web-Links:

Ein paar Map/Reduce Problem Unterteilung Lösungen

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

A out-of-Kern (PyTables) Aufbewahrungslösung

• Very large matrices using Python and NumPy

A matmul Unterteilung Lösung:

• 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/

Vielen Dank im Voraus für jede Empfehlungen, Kommentare oder Anleitung!

Answer 1

Da Sie nach dem Produkt einer Matrix mit ihren Transponierten sind, wird der Wert bei [m, n] im Grunde geht das Punktprodukt der Spalten m und n in Ihrer ursprünglichen Matrix sein.

Ich werde die folgende Matrix als Spielzeug Beispiel

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]]) 

Es ist die Form (3, 12) und 7 Nicht-Null-Einträge verwenden. Das Produkt seiner Transponierung ist natürlich von der Form (12, 12) und hat 16 Einträge ungleich null, davon 6 in der Diagonalen, so dass nur 11 Elemente gespeichert werden müssen.

Sie können eine gute Vorstellung davon, was die Größe der Ausgabematrix in eine von zwei Arten sein wird:

CSR FORMAT

Wenn die ursprüngliche Matrix C Nicht-Null-Spalten , Ihre neue Matrix hat höchstens C**2 Einträge ungleich Null, von denen C in der Diagonalen sind, und es wird sichergestellt, dass sie nicht null sind, und von den restlichen Einträgen müssen Sie nur die Hälfte behalten, also höchstens (C**2 + C)/2 Null Elemente. Natürlich werden viele von ihnen auch Null sein, also ist dies wahrscheinlich eine grobe Überschätzung.

Wenn die Matrix in csr Format gespeichert ist, dann ist das indices Attribut des entsprechenden scipy Objekts hat eine Reihe mit den Spaltenindizes aller Nicht-Null-Elemente, so dass man leicht die obige Schätzung als berechnen:

>>> 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 

so gibt es 6 Spalten mit nicht-Null-Einträge, und so würde die Schätzung für höchstens 36 nicht-Null-Einträge sein, weit mehr als die reale 16.

CSC FORMAT

Wenn wir anstelle von Spaltenindizes von Elementen ungleich Null Zeilenindizes haben, können wir tatsächlich eine bessere Schätzung vornehmen. Damit das Skalarprodukt zweier Spalten nicht null ist, müssen sie ein Element ungleich null in derselben Zeile enthalten. Wenn R Nicht-Null-Elemente in einer bestimmten Zeile vorhanden sind, tragen sie R**2 Nicht-Null-Elemente zum Produkt bei. Wenn Sie dies für alle Zeilen summieren, müssen Sie einige Elemente mehr als einmal zählen, also ist dies auch eine obere Grenze.

Der Reihenindizes der Nicht-Null-Elemente Ihrer Matrix sind in dem indices Attribute einer sparse csc-Matrix, so kann dieser Schätzwert wie folgt berechnet werden:

>>> 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 

Dies ist stopfen Nähe des realen 16! Und es eigentlich kein Zufall ist, dass diese Schätzung ist eigentlich das gleiche wie:

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

In jedem Fall sollte der folgende Code ermöglicht es Ihnen, zu sehen, dass die Schätzung ziemlich gut:

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