NumPy: calcolare x.T * x grande matrice

Question

In numpy, qual è il modo più efficace per x.T * x di elaborazione, dove x è un grande (200.000 x 1000) a matrice float32 denso e .T è l'operatore di trasposizione?

Per evitare ogni dubbio, il risultato è 1000 x 1000.

modifica : Nella mia domanda originale ho dichiarato che np.dot(x.T, x) stava prendendo ore. Venne fuori che avevo un po 'sornione NaNs nella matrice, e per qualche motivo che è stato completamente uccidendo le prestazioni di np.dot (eventuali approfondimenti sul perché?) Questo è ora risolto, ma la domanda originale stand.

Answer 1

Questa non può essere la risposta che stai cercando, ma un modo per accelerarlo considerevolmente è quello di utilizzare una GPU al posto della CPU. Se si dispone di una scheda grafica potente decentemente in giro, sarà sorpassare la CPU ogni giorno, anche se il sistema è molto ben sintonizzato.

Per l'integrazione con bel NumPy, si potrebbe usare Theano (se la vostra scheda grafica è fatto da NVIDIA). Il calcolo nei seguenti percorsi di codice per me in paio di secondi (anche se ho una molto potente scheda grafica):

$ THEANO_FLAGS=device=gpu0 python
Python 2.6.5 (r265:79063, Apr 16 2010, 13:57:41) 
[GCC 4.4.3] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import theano
Using gpu device 0: GeForce GTX 480
>>> from theano import tensor as T
>>> import numpy
>>> x = numpy.ones((200000, 1000), dtype=numpy.float32)
>>> m = T.matrix() 
>>> mTm = T.dot(m.T, m)
>>> f = theano.function([m], mTm)
>>> f(x)
array([[ 200000.,  200000.,  200000., ...,  200000.,  200000.,  200000.],
       [ 200000.,  200000.,  200000., ...,  200000.,  200000.,  200000.],
       [ 200000.,  200000.,  200000., ...,  200000.,  200000.,  200000.],
       ..., 
       [ 200000.,  200000.,  200000., ...,  200000.,  200000.,  200000.],
       [ 200000.,  200000.,  200000., ...,  200000.,  200000.,  200000.],
       [ 200000.,  200000.,  200000., ...,  200000.,  200000.,  200000.]], dtype=float32)
>>> r = f(x)
>>> r.shape
(1000, 1000)

I aveva intenzione di aspettare per scoprire quanto tempo >>> numpy.dot(x.T, x) preso come termine di paragone, ma mi sono stufato ...

Si può anche provare PyCuda o PyOpenCL (se non si dispone di una scheda grafica NVIDIA), anche se non so se il loro sostegno NumPy è così semplice.

Answer 2

Per prima cosa, assicurarsi di utilizzare un Blas ottimizzati / LAPACK, questo può fare una differenza enorme (fino a un ordine di grandezza). Se si utilizza un filettata ATLAS, per esempio, userà tutti i nuclei relativamente efficiente (è necessario utilizzare un recente ATLAS, però, e la compilazione di Atlas è una valle di lacrime).

Per quanto riguarda il motivo per cui Nan rallenta tutto fatto: che è praticamente inevitabile, la gestione Nan è molto più lento di galleggiare "normale" a livello di CPU: http://www.cygnus-software.com/papers/x86andinfinity.html . Dipende dal modello di CPU, che tipo di set di istruzioni che si sta utilizzando, e naturalmente l'algoritmi / implementazione che si sta utilizzando.

Answer 3

hmm, x è di circa 800 Mb, ammesso che ha bisogno dello stesso per il risultato, sei sicuro di avere abbastanza memoria fisica e non è lo scambio?

diverso da quello, NumPy dovrebbe utilizzare una funzione BLAS, e anche se la libreria di default che utilizza NumPy possono essere relativamente lenta, dovrebbe funzionare bene per questo formato.

modifica

import numpy as npy
import time

def mm_timing():
  print "   n   Gflops/s"
  print "==============="
  m = 1000
  n = 200000
  a = npy.random.rand(n, m)
  flops = (2 * float(n) - 1) * float(m)**2
  t1 = time.time()
  c = npy.dot(a.T, a)
  t2 = time.time()
  perf = flops / (t2 - t1) / 1.e9
  print "%4i" % n + "     " + "%6.3f" % perf

mm_timing()