python - opt - Encontrar la matriz de correlación

scipy optimize minimize in python

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 un código de Python como este:

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

¿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 esté recordando mal exactamente qué es Pearson''s R, que es bastante posible).

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

Por ejemplo, probar 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()

Produce 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 debería hacer exactamente lo que quieres, y es mucho más rápido.

Answer 2

puedes usar el módulo de multiprocesamiento de python, dividir tus filas en 10 grupos, almacenar tus resultados y luego imprimirlos (esto solo aceleraría en una máquina multinúcleo)

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

Por cierto, también deberías convertir tu fragmento en una función y también considerar cómo volver a armar los datos. tener cada subproceso con una lista como esta ... [startcord, stopcord, buff] ... podría funcionar bien

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

Answer 3

Nueva solución

Después de mirar la respuesta de Joe Kington, decidí investigar el código de corrcoef() y me inspiré para hacer 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 bucle genera los coeficientes de Pearson entre la fila i y 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.

Vieja solución

Primero, no recomendaría imprimir las r en 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 su 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 de manera redundante algunas cantidades. 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 está utilizando el valor p (que también devuelve pearsonr() así que vamos a borrar eso también. Usando 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)

Obtuve una aceleración de 4.8x sobre el código scipy directo cuando eliminé el valor p -8.8x si dejo el material p-allí (utilicé 10 columnas con cientos de filas). También verifiqué que da los mismos resultados. Esta 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 tomando 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.