numpy: calcular x.T * x para una gran matriz

Question

En numpy, ¿cuál es la forma más eficiente para calcular x.T * x, donde x es un grande (200.000 x 1000) float32 densa matriz y .T es el operador de transposición?

Para evitar dudas, el resultado es de 1000 x 1000.

edición: En mi pregunta original afirmé que estaba tomando np.dot(x.T, x) horas. Resultó que tenía un poco de NaNs en la matriz, y por alguna razón eso estaba matando por completo el rendimiento de np.dot (¿Alguna idea de por qué?) Esto ahora está resuelto, pero la pregunta original permanece.

Answer 1

Esto puede no ser la respuesta que está buscando, pero una forma de acelerarlo considerablemente es el uso de una GPU en lugar de su CPU. Si tiene una tarjeta de gráficos decentemente poderosa, superará a su CPU en cualquier momento, incluso si su sistema está muy bien ajustado.

Para una buena integración con numpy, puede usar theano (si su tarjeta gráfica está hecha por nvidia). El cálculo en el siguiente código funciona para mí en un par de segundos (aunque tengo un muy potente tarjeta gráfica):

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

que iba a esperar para saber cuánto tiempo >>> numpy.dot(x.T, x) tomó a modo de comparación, pero se aburrió ...

También puede probar PyCuda o PyOpenCL (si no tiene una tarjeta gráfica nvidia), aunque no sé si su soporte numpy es tan sencillo.

Answer 2

hmm, x es aproximadamente 800 Mb, suponiendo que necesita el mismo resultado, ¿está seguro de que tiene suficiente memoria física y no está intercambiando?

Aparte de eso, numpy debería usar una función BLAS, y aunque la biblioteca predeterminada que Numpy usa puede ser relativamente lenta, debería funcionar bien para este tamaño.

edición

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

Answer 3

En primer lugar, asegúrese de utilizar un blas/lapack optimizado, esto puede hacer una gran diferencia (hasta un orden de magnitud). Si utiliza un ATLAS enhebrado, por ejemplo, utilizará todos sus núcleos de manera relativamente eficiente (necesita utilizar un ATLAS reciente, sin embargo, y compilar ATLAS es un PITA).

En cuanto a por qué Nan ralentiza todo: es casi inevitable, el manejo de NaN es mucho más lento que el float "normal" en el nivel de la CPU: http://www.cygnus-software.com/papers/x86andinfinity.html. Depende del modelo de CPU, qué tipo de conjunto de instrucciones está utilizando y, por supuesto, los algoritmos/implementación que está utilizando.