Поиск корреляционной матрицы

Ответ 1

Новое решение

Посмотрев на ответ Джо Кингтона, я решил изучить код corrcoef() и был вдохновлен им сделать следующую реализацию.

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

Каждый цикл через генерирует коэффициенты Пирсона между строками я и строками я до последней строки. Это очень быстро. Он по крайней мере в 1,5 раза быстрее, чем при использовании corrcoef(), потому что он не избыточно вычисляет коэффициенты и несколько других вещей. Он также будет быстрее и не даст вам проблем с памятью с матрицей строк в 50 000, потому что тогда вы можете либо сохранить каждый набор r, либо обработать их до создания другого набора. Не сохраняя какой-либо из долгосрочных обязательств, я смог получить приведенный выше код для запуска на 50 000 х 10 наборов случайно генерируемых данных за минуту на моем довольно новом ноутбуке.

Старое решение

Во-первых, я бы не рекомендовал печатать r на экране. Для 100 строк (10 столбцов) это разница в 19,79 секунды с печатью против 0,301 секунды без использования вашего кода. Просто сохраните r и используйте их позже, если хотите, или выполните некоторую обработку на них, когда вы идете вперед, ища некоторые из самых больших r.

Во-вторых, вы можете получить некоторую экономию, не избыточно вычисляя некоторые количества. Коэффициент Пирсона рассчитывается в scipy с использованием некоторых величин, которые вы можете предварительно вычислять, а не вычислять каждый раз, когда используется строка. Кроме того, вы не используете значение p (которое также возвращается pearsonr(), так что пусть оно тоже царапается). Используя приведенный ниже код:

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)

Я получаю ускорение около 4.8x по прямому scipy-коду, когда я удалял материал p-value - 8.8x, если я оставил там p-значение (я использовал 10 столбцов с сотнями строк). Я также проверил, что он дает те же результаты. Это не очень большое улучшение, но это может помочь.

В конечном счете, вы столкнулись с проблемой, которую вы вычисляете (50000) * (50001)/2 = 1 250 025 000 коэффициентов Пирсона (если я правильно рассчитываю). Это много. Кстати, нет необходимости вычислять каждый коэффициент Пирсона с самим собой (он будет равен 1), но это только избавит вас от вычисления 50 000 коэффициентов Пирсона. С приведенным выше кодом я ожидаю, что для выполнения ваших вычислений потребуется около 4 1/4 часа, если у вас есть 10 столбцов для ваших данных, основанных на моих результатах по более мелким наборам данных.

Вы можете получить некоторое улучшение, взяв вышеуказанный код в Cython или что-то подобное. Я ожидаю, что вы, возможно, получите 10-кратное улучшение по сравнению с прямым Scipy, если вам повезет. Кроме того, как было предложено pyInTheSky, вы можете выполнить некоторую многопроцессорную обработку.

Ответ 2

Вы пробовали использовать numpy.corrcoef? Видя, как вы не используете p-значения, он должен делать именно то, что вы хотите, с наименьшими проблемами. (Если я не помню точно, что такое pearson R, что вполне возможно.)

Просто быстро проверяя результаты на случайных данных, он возвращает точно такую ​​же вещь, как и код @Justin Peel выше и работает ~ 100 раз быстрее.

Например, тестирование вещей с 1000 строк и 10 столбцов случайных данных...:

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 built-in function"""
    return np.corrcoef(data)

def justin_peel_test(data):
    """Justin Peel 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()

Допускается максимальная абсолютная разница ~ 3.3e-16 между двумя результатами

И тайминги:

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 должен делать именно то, что вы хотите, и это намного быстрее.

Ответ 3

вы можете использовать модуль мультипроцесса python, разбивать строки на 10 наборов, буферизировать результаты и затем распечатывать материал (это только ускорит его на многоядерном компьютере)

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

btw: вам также придется превратить ваш фрагмент в функцию, а также подумайте, как выполнить сборку данных. у каждого подпроцесса есть список, подобный этому... [startcord, stopcord, buff].. может работать красиво

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