moltiplicazione della matrice numpy in memoria triangolare/sparsa?

Question

Sto lavorando con un problema di moltiplicazione (matmul) di matrice sparse molto grande. Ad esempio, diciamo:

• A è una matrice binaria (75 x 200.000). È scarso, quindi sto usando csc per l'archiviazione. Devo eseguire la seguente operazione matmul:

• B = A.transpose() * Un

• L'uscita sarà una matrice sparsa e simmetrica di dimensione 200Kx200K.

Purtroppo, B sta per essere modo alla grande per memorizzare nella RAM (o "in core") sul mio portatile. D'altra parte, sono fortunato perché ci sono alcune proprietà di B che dovrebbero risolvere questo problema.

Poiché B sta per essere simmetrico lungo la diagonale e sparse, potrei usare una matrice triangolare (superiore/inferiore) per memorizzare i risultati dell'operazione di matmul e un formato di archiviazione di matrice sparsa potrebbe ridurre ulteriormente le dimensioni.

La mia domanda è ... si può dire a numpy o scipy, prima del tempo, quali saranno i requisiti di archiviazione di output in modo che io possa selezionare una soluzione di archiviazione usando numpy ed evitare che la "matrice sia troppo grande" errore di runtime dopo alcuni minuti (ore) di calcolo?

In altre parole, i requisiti di memoria per il moltiplicatore di matrice possono essere approssimati analizzando il contenuto delle due matrici di input utilizzando un algoritmo di conteggio approssimativo?

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

In caso contrario, sto cercando una soluzione di forza bruta. Qualcosa che coinvolge Map/Reduce, out-of-core di archiviazione, o una soluzione matmul suddivisione (Algoritmo di Strassen) dai seguenti link web:

Un paio Mappa/Ridurre soluzioni di suddivisione problema

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

Una soluzione (PyTables) stoccaggio out-of-core

• Very large matrices using Python and NumPy

Una soluzione suddivisione matmul:

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

Grazie in anticipo per qualsiasi suggerimenti, commenti, o guida!

Answer 1

Dato che si utilizza il prodotto di una matrice con la sua trasposizione, il valore di [m, n] sarà fondamentalmente il prodotto di punti delle colonne m e n nella matrice originale.

ho intenzione di utilizzare la seguente matrice come esempio giocattolo

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

È forma (3, 12) e dispone di 7 voci non nulli. Il prodotto della sua trasposizione con esso è naturalmente forma (12, 12) e ha 16 voci diverse da zero, 6 delle quali nella diagonale, quindi richiede solo la memorizzazione di 11 elementi.

è possibile ottenere una buona idea di ciò che la dimensione del vostro matrice di uscita sta per essere in uno dei due modi:

CSR FORMATO

Se la matrice originale ha C colonne non-zero , la tua nuova matrice avrà al massimo C**2 voci diverse da zero, di cui C sono nella diagonale, e sono assicurati di non essere zero, e delle voci rimanenti devi solo tenere la metà, quindi è al massimo (C**2 + C)/2 non- zero elementi. Naturalmente, molti di questi saranno anche zero, quindi probabilmente è una sovrastima lorda.

Se la matrice viene memorizzato in formato csr, l'attributo indices del corrispondente scipy oggetto è un array con gli indici colonna di tutti gli elementi non nulli, così si può facilmente calcolare la stima precedente come:

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

quindi ci sono 6 colonne con ingressi non-zero, e quindi la stima sarebbero per al massimo 36 voci diverse da zero, modo più di quello reale 16.

CSC FORMATO

Se invece di indici di colonna di elementi diversi da zero abbiamo indici di riga, possiamo effettivamente fare una stima migliore. Perché il prodotto punto di due colonne sia diverso da zero, deve avere un elemento diverso da zero nella stessa riga. Se ci sono R elementi diversi da zero in una determinata riga, essi contribuiranno con elementi non zero a R**2. Quando sommi questo per tutte le righe, sei tenuto a contare alcuni elementi più di una volta, quindi questo è anche un limite superiore.

Gli indici schiera di elementi non nulli del tuo matrice sono nell'attributo indices di una matrice CSC sparsa, quindi questa stima può essere calcolato come segue:

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

Questo è davvero vicino al reale 16! E non è in realtà un caso che questa stima è in realtà lo stesso di:

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

In ogni caso, il seguente codice dovrebbe consentire di vedere che la stima è piuttosto buona:

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