Sparse-Matrix-Vektor-Produkt von Hand ausrechnen

Wir wollen nun einige erste Schritte machen, um zu sehen, wie man mit speziellen Datenstrukturen für dünnbesetzte Matrizen arbeiten kann.

Das Matrizen aus numpy dicht besetzt sind, speichern diese jeden Eintrag der Matrix, also auch die die Einträge welche 0 sind. In einem Matrix-Vektorprodukt, liefern diese Einträge aber keinen Beitrag zum Ergebnis, daher möchten wir nun versuchen diese zu vermeiden.

Die scipy Bibliothek stellt im Paket sparse verschiedene Algorithmen und Datenstrukturen zur Verfügung um mit Sparse-Matrizen zu arbeiten:

import numpy
from scipy import sparse

Ein Sparse-Matrizen Datenformat soll einen leichten Zugriff auf alle Nicht-Null-Einträge erlauben. Wir wollen dies exemplarisch für ein Matrix-Vektor-Produkt nutzen.

Die folgende Version ist zwar langsam, da wir alles in python schreiben, aber die zeigt die Konzeptionelle Idee:

def spmv(M,v):
   x = numpy.zeros(v.shape)

Wir iterieren über alle Zeilen:

for r in range(M.shape[0]) :
    row = M[r]

und holen uns die Liste der Spaltenindizes aller Nicht-Null-Einträge in Zeile r:

cols = row.nonzero()[1]

dann iterieren wir über alle Nicht-Null-Einträge und berechnen den Beitrag zur Ergbenis:

    for c in cols:
        x[r] = x[r] + M[r,c] * v[c]
return x

Jetzt wollen wir unsere Implementierung testen und mit der Implementierung von scipy vergleichen

Wir legen eine Matrix im CSR Format an:

B = sparse.csr_matrix(
   [[0,3,1,0],
    [0,2,0,0],
    [1,0,0,0],
    [0,0,0,1]],dtype="double")

Und einen Vektor als numpy.array:

v = numpy.array(range(4))

Wir können beides ausgeben und das Matrix-Vektor-Produkt berechnen:

print 'B =',B.todense()
print 'v =',v

zunächst unsere Implementierung:

print 'by hand: B * v =',spmv(B,v)

und dann die scipy Implementierung:

print 'scipy:   B * v =',B*numpy.array(v)