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

У меня есть матрица, которая довольно большая (около 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 ).

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

Благодаря!

3 Solutions collect form web for “Поиск корреляционной матрицы”

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

Посмотрев на ответ Джо Кингтона, я решил изучить код 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() быстрее, чем использование 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 коэффициентов Пирсона (если я правильно рассчитываю). Это много. Кстати, нет необходимости вычислять коэффициент Pearson каждой строки с самим собой (он будет равен 1), но это только избавит вас от вычисления 50 000 коэффициентов Пирсона. С приведенным выше кодом я ожидаю, что для вычисления будет потребоваться около 4 1/4 часа, если у вас есть 10 столбцов для ваших данных, основанных на моих результатах на меньших наборах данных.

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

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

вы можете использовать модуль мультипроцесса 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 
  • Failed scipy.special import "Символ не найден: ___addtf3"
  • Как эффективно вычислять огромное умножение матрицы (функции tfidf) в Python?
  • scipy.misc.imshow RuntimeError ('Не удалось выполнить просмотр изображения')
  • Как ускорить создание матрицы перехода в Numpy?
  • Возможно ли воспроизвести randn () MATLAB с NumPy?
  • python (scipy): изменение размера разреженной матрицы
  • Заполните Pandas SparseDataFrame из SciPy Sparse Matrix
  • Последовательные, перекрывающиеся подмножества массива (NumPy, Python)
  •  
    Interesting Posts for Van-Lav

    py2exe с matplotlib, numpy и pylab

    Расширьте argparse, чтобы записывать имена заданий в тексте справки для выбора опциональных аргументов и определять эти наборы один раз в конце

    Не удалось запустить scons, получив ошибку импорта

    Python: сделайте видео, используя несколько изображений .png

    Отображение обратного отсчета для функции ожидания python

    Добавление двух разреженных матриц `csc` разной формы в python

    Можем ли мы запустить многопроцессорный пул в GAE?

    Python interp1d против UnivariateSpline

    Как я могу сделать Selenium щелчком по переменному числу «следующих» кнопок?

    Почему мой скрипт python случайно убит?

    многострочный отступ python на emacs

    Python metaclasses: Почему не __setattr__ вызывается для атрибутов, установленных во время определения класса?

    создать матрицу из массива элементов под диагональю в numpy

    Преобразование списка в вложенный словарь

    Многопроцессорность Python для параллельных процессов

    Python - лучший язык программирования в мире.