Encontrando la matriz de correlación

Question

Tengo una matriz que es bastante grande (alrededor de 50K filas), y quiero imprimir el coeficiente de correlación entre cada fila en la matriz. He escrito el código Python como esto:

for i in xrange(rows): # rows are the number of rows in the matrix. 
    for j in xrange(i, rows): 
     r = scipy.stats.pearsonr(data[i,:], data[j,:]) 
     print r 

Tenga en cuenta que estoy haciendo uso de la función pearsonr disponible en el módulo scipy (http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.pearsonr.html).

Mi pregunta es: ¿hay una manera más rápida de hacer esto? ¿Hay alguna técnica de partición matricial que pueda usar?

Gracias!

Answer 1

puede utilizar el módulo de Python multiproceso, trozo de seguridad de sus filas en 10 juegos, amortiguar sus resultados y luego imprimir la materia hacia fuera (esto sólo acelerarlo en una máquina de múltiples núcleos, aunque)

http://docs.python.org/library/multiprocessing.html

por cierto: también tendría que convertir su fragmento en una función y también considerar cómo hacer el reensamblaje de datos. teniendo cada subproceso haber una lista como esta ... [startcord, stopcord, buff] .. podría funcionar muy bien

def myfunc(thelist): 
    for i in xrange(thelist[0]:thelist[1]): 
    .... 
    thelist[2] = result 

Answer 2

nueva solución

Después de mirar la respuesta de Joe Kington, decidí estudiar el código corrcoef() y se inspiró para realizar la siguiente implementación.

ms = data.mean(axis=1)[(slice(None,None,None),None)] 
datam = data - ms 
datass = np.sqrt(scipy.stats.ss(datam,axis=1)) 
for i in xrange(rows): 
    temp = np.dot(datam[i:],datam[i].T) 
    rs = temp/(datass[i:]*datass[i]) 

Cada ciclo genera los coeficientes de Pearson entre la fila iy las filas i hasta la última fila. Es muy rápido. Es al menos 1.5 veces más rápido que usar corrcoef() solo porque no calcula de forma redundante los coeficientes y algunas otras cosas. También será más rápido y no le dará los problemas de memoria con una matriz de 50,000 filas porque entonces puede elegir almacenar cada conjunto de r o procesarlas antes de generar otro conjunto. Sin almacenar ninguno de los r a largo plazo, pude obtener el código anterior para ejecutar en 50,000 x 10 conjunto de datos generados aleatoriamente en menos de un minuto en mi portátil bastante nuevo.

antigua solución

En primer lugar, yo no recomendaría imprimir el r de la pantalla. Para 100 filas (10 columnas), esta es una diferencia de 19.79 segundos con la impresión frente a 0.301 segundos sin usar el código. Simplemente almacene las "r" y úselas más adelante si lo desea, o haga algún procesamiento con ellas a medida que avance, como buscar algunas de las r más grandes.

En segundo lugar, puede obtener algunos ahorros al no calcular algunas cantidades de forma redundante. El coeficiente de Pearson se calcula en scipy usando algunas cantidades que puede precalcular en lugar de calcular cada vez que se utiliza una fila. Además, no se está utilizando el valor de p (que también es devuelto por pearsonr() así que vamos a rascar eso también con el siguiente código:.

r = np.zeros((rows,rows)) 
ms = data.mean(axis=1) 

datam = np.zeros_like(data) 
for i in xrange(rows): 
    datam[i] = data[i] - ms[i] 
datass = scipy.stats.ss(datam,axis=1) 
for i in xrange(rows): 
    for j in xrange(i,rows): 
     r_num = np.add.reduce(datam[i]*datam[j]) 
     r_den = np.sqrt(datass[i]*datass[j]) 
     r[i,j] = min((r_num/r_den), 1.0) 

consigo una aceleración de alrededor de 4,8 veces por encima del scipy recta código cuando eliminé el valor de p-cosas - 8.8x si dejo las cosas de valor p allí (utilicé 10 columnas con cientos de filas). También verifiqué que da los mismos resultados. Esto no es una gran mejora, pero podría ayudar.

En última instancia, está atascado con el problema de que está calculando (50000) * (50001)/2 = 1,250,025,000 coeficientes de Pearson (si estoy contando correctamente). Eso es mucho. Por cierto, realmente no hay necesidad de calcular el coeficiente de Pearson de cada fila consigo mismo (será igual a 1), pero eso solo le ahorra el cálculo de 50,000 coeficientes de Pearson. Con el código anterior, espero que tome aproximadamente 4 1/4 horas para realizar su cálculo si tiene 10 columnas para sus datos en función de mis resultados en conjuntos de datos más pequeños.

Puede obtener alguna mejora si toma el código anterior en Cython o algo similar. Espero que tengas una mejora de hasta 10 veces con respecto a Scipy si tienes suerte. Además, según lo sugerido por pyInTheSky, puede hacer un multiprocesamiento.

Answer 3

¿Has probado usar numpy.corrcoef? Viendo que no estás usando los valores p, debería hacer exactamente lo que quieras, con el mínimo esfuerzo posible. (A menos que recuerde exactamente qué es Pearson's R, que es bastante posible).

Al comprobar rápidamente los resultados en datos aleatorios, devuelve exactamente lo mismo que el código de @Justin Peel anterior y se ejecuta ~ 100x más rápido .

Por ejemplo, probando cosas con 1000 filas y 10 columnas de datos aleatorios ...:

import numpy as np 
import scipy as sp 
import scipy.stats 

def main(): 
    data = np.random.random((1000, 10)) 
    x = corrcoef_test(data) 
    y = justin_peel_test(data) 
    print 'Maximum difference between the two results:', np.abs((x-y)).max() 
    return data 

def corrcoef_test(data): 
    """Just using numpy's built-in function""" 
    return np.corrcoef(data) 

def justin_peel_test(data): 
    """Justin Peel's suggestion above""" 
    rows = data.shape[0] 

    r = np.zeros((rows,rows)) 
    ms = data.mean(axis=1) 

    datam = np.zeros_like(data) 
    for i in xrange(rows): 
     datam[i] = data[i] - ms[i] 
    datass = sp.stats.ss(datam,axis=1) 
    for i in xrange(rows): 
     for j in xrange(i,rows): 
      r_num = np.add.reduce(datam[i]*datam[j]) 
      r_den = np.sqrt(datass[i]*datass[j]) 
      r[i,j] = min((r_num/r_den), 1.0) 
      r[j,i] = r[i,j] 
    return r 

data = main() 

presenta una diferencia absoluta máxima de ~ 3.3e-16 entre los dos resultados

y tiempos :

In [44]: %timeit corrcoef_test(data) 
10 loops, best of 3: 71.7 ms per loop 

In [45]: %timeit justin_peel_test(data) 
1 loops, best of 3: 6.5 s per loop 

numpy.corrcoef debe hacer precisamente lo que quiere, y es mucho más rápido.