numpy Korrelationskoeffizient: np.dot (A, AT) auf großen Arrays verursacht seg fault

Question

HINWEIS:

Speed is not as important as getting a final result. 
However, some speed up over worst case is required as well. 

Ich habe ein großes Array A:

A.shape=(20000,265) # or possibly larger like 50,000 x 265 

ich brauche die Korrelation zu berechnen, Koeffizienten.

np.corrcoeff # internally casts the results as doubles 

Ich lieh nur ihren Code und schrieb meine eigene cov/Korr nicht in Doppel Gießen, da ich nur 32-Bit wirklich brauchen floats.And ich Graben der conj(), da meine Daten immer real sind.

cov = A.dot(A.T)/n #where A is an array of 32 bit floats 
diag = np.diag(cov) 
corr = cov/np.sqrt(np.mutliply.outer(d,d)) 

ich immer noch über genügend Arbeitsspeicher ausgeführt, und ich bin eine große Speichermaschine, 264GB

ich gesagt habe, dass die schnellen C-Bibliotheken sind wahrscheinlich eine Routine, die den Punkt bricht Produkt in Stücke, und um dies zu optimieren, wird die Anzahl der Elemente auf eine Stärke von 2 aufgefüllt.

Ich brauche nicht wirklich die symmetrische Hälfte der Korrelationskoeffizientenmatrix zu berechnen. Allerdings sehe ich keinen Weg, dies in angemessener Zeit manuell mit Python-Schleifen zu tun.

Kennt jemand einen Weg, um numpy für eine annehmbare Punktproduktroutine zu bitten, die Speicherverbrauch mit Geschwindigkeit ausgleicht ...?

Prost

UPDATE:

Lustig, wie diese Fragen schreiben Sie mir die Sprache für eine bessere Google-Abfrage finden zu helfen, neigt.

diese Gefunden:

http://wiki.scipy.org/PerformanceTips 

nicht sicher, dass ich es folgen .... ja, bitte Antworten zu dieser Lösung kommentieren oder bieten, Ihre eigenen Ideen oder einfach nur allgemeine Kommentare zu dieser Art von Problem.

TIA

EDIT: Ich entschuldige mich, weil mein Array wirklich ist viel größer, als ich dachte. Array-Größe ist eigentlich 151.000 x 265 Ich bin auf einer Maschine mit 264 GB mit mindestens 230 GB frei von Arbeitsspeicher.

Ich bin überrascht, dass die numpy Call zu Blas dgemm und vorsichtig mit C-Order-Arrays nicht hocken.

Answer 1

Sie können versuchen und sehen, ob np.einsum besser funktioniert als Punkt für Ihren Fall:

cov = np.einsum('ij,kj->ik', A, A)/n 

Die internen Abläufe dot sind ein wenig unklar, wie es BLAS optimierte Routinen zu verwenden versucht, die manchmal erfordern Kopien Arrays in Fortran-Reihenfolge, nicht sicher, ob das hier der Fall ist. einsum puffert seine Eingaben und verwendet, wenn möglich, vektorisierte SIMD-Operationen, aber außerhalb davon werden im Grunde die naiven drei verschachtelten Schleifen ausgeführt, um das Matrixprodukt zu berechnen.

Answer 2

Python mit Intel MKL kompiliert wird dies in etwa 30 Sekunden mit 12 GB Speicher laufen:

>>> A = np.random.rand(50000,265).astype(np.float32) 
>>> A.dot(A.T) 
array([[ 86.54410553, 64.25226593, 67.24698639, ..., 68.5118103 , 
     64.57299805, 66.69223785], 
     ..., 
     [ 66.69223785, 62.01016235, 67.35866547, ..., 66.66306305, 
     65.75863647, 86.3017807 ]], dtype=float32) 

Wenn Sie Zugriff auf in Intels MKL herunterladen Python haben nicht anaconda und installieren Sie das Paket beschleunigen, die eine Studie hat Version für 30 Tage oder kostenlos für Akademiker, die eine mkl kompiliert enthält. Verschiedene andere C++ - BLAS-Bibliotheken sollten auch funktionieren - selbst wenn es das Array von C nach F kopiert, sollte es nicht mehr als 30 GB Speicher benötigen.

Das einzige, was ich denken kann, dass Ihre Installation versucht, ist zu versuchen, die gesamte 50.000 x 50.000 x 265-Array im Speicher zu halten, die ganz offen gesagt schrecklich ist. Als Referenz ein float32 50.000 x 50.000 Array ist nur 10 GB, während die zuvor erwähnte Anordnung 2.6TB ist ...

Wenn sie ein gemm Problem, das Sie ein Stück gemm Formel ausprobieren kann:

def chunk_gemm(A, B, csize): 
    out = np.empty((A.shape[0],B.shape[1]), dtype=A.dtype) 
    for i in xrange(0, A.shape[0], csize): 
     iend = i+csize 
     for j in xrange(0, B.shape[1], csize): 
      jend = j+csize 
      out[i:iend, j:jend] = np.dot(A[i:iend], B[:,j:jend]) 
    return out 

Dies wird langsamer , aber wird hoffentlich über Ihre Speicherprobleme hinwegkommen.

Answer 3

UPDATE: Stellt fest, dass das Skalarprodukt ohne Fehler abgeschlossen wurde, aber nach sorgfältiger Prüfung besteht das Ausgangsarray aus Nullen bei 95.000 bis zum Ende der 151.000 Spalten. Das heißt, daß aus ist [:, 94.999] = Nicht-Null, sondern aus [:, 95000] = 0 für alle Zeilen ...

Das ist super ärgerlich ...

Another Blas description

Die Austausch, erwähnt etwas, über das ich auch dachte ... Da Blas Fortran ist, sollte nicht die Reihenfolge der Eingabe F bestellen ...? Wo wie die scpy doc Seite unten, sagt C-Reihenfolge.

Der Versuch F-Befehl verursachte einen Segmentierungsfehler. Also bin ich zurück auf Platz eins.

ORIGINAL POST Ich fand endlich mein Problem, das wie üblich in den Details war.

Ich verwende ein Array von np.float32, die als F-Reihenfolge gespeichert wurden. Ich kann das F meines Wissens nicht kontrollieren, da die Daten aus Bildern mit einer Bildgebungsbibliothek geladen werden.

import scipy 
roi = np.ascontiguousarray(roi)# see roi.flags below 
out = scipy.linalg.blas.sgemm(alpha=1.0, a=roi, b=roi, trans_b=True) 

Diese Stufe 3 Blas-Routine macht den Trick. Mein Problem war zweifach:

roi.flags 
C_CONTIGUOUS : False 
F_CONTIGUOUS : True 
OWNDATA : True 
WRITEABLE : True 
ALIGNED : True 
UPDATEIFCOPY : False 

Und ... ich war blas dgemm nicht sgemm. Das 'd' ist für 'doppelt' und 's' für 'single'. Sehen Sie diese pdf: BLAS summary pdf

ich es früher einmal ausgesehen und war überwältigt ... Ich ging zurück und die Wikipedia-Artikel über blas Routinen lesen Level 3 vs anderen Ebenen zu verstehen: wikipedia article on blas

Jetzt funktioniert es auf A = 150.000 x 265, Durchführung:

Vielen Dank für Ihre Gedanken ... zu wissen, dass es getan werden könnte, war am wichtigsten.