Нахождение корреляционной матрицы

Question

У меня довольно большая матрица (около 50 тыс. Строк), и я хочу напечатать коэффициент корреляции между каждой строкой в ​​матрице.Я написал такой код на Python:

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  

Обратите внимание, что я использую функцию pearsonr, доступную в модуле scipy (http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.pearsonr.html).

. Мой вопрос:Есть ли более быстрый способ сделать это? Есть ли какой-нибудь метод разбиения матрицы, который я могу использовать?

Спасибо!

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

Каждый проходной цикл генерирует коэффициенты Пирсона между строкой i и строкой i до последней строки. Это очень быстро. Это как минимум в 1,5 раза быстрее, чем при использовании только corrcoef(), потому что он не требует избыточного вычисления коэффициентов и некоторых других вещей. Это также будет быстрее и не вызовет проблем с памятью с матрицей в 50 000 строк, потому что тогда вы можете либо сохранить каждый набор r, либо обработать их перед созданием другого набора. Не сохраняя ни одного долгосрочного кода, я смог заставить вышеуказанный код работать на множестве случайно сгенерированных данных размером 50 000 x 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,8 раза по сравнению с простым кодом scipy, когда я удаляю материал с p-значением - 8,8 раза, если я оставляю материал с p-значением (я использовал 10 столбцов с сотнями строк) ). Я также проверил, что это дает те же результаты. Это не очень большое улучшение, но оно может помочь.

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

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

Answer 2

Вы пробовали просто использовать numpy.corrcoef ? Видя, как вы не используете p-значения, он должен делать именно то, что вы хотите, с минимальным суетой, насколько это возможно. (Если я неправильно помню, что такое 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'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()

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

Answer 3

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

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

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

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