找到相关矩阵

Question

我有一个相当大的矩阵(大约50K行),我想打印矩阵中每行之间的相关系数.我编写了这样的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  

请注意,我正在使用pearsonrscipy模块提供的功能(http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.pearsonr.html).

我的问题是:有更快的方法吗？我可以使用一些矩阵分区技术吗？

谢谢!

Answer 1

新解决方案

在查看了Joe Kington的回答之后,我决定查看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到最后一行之间生成Pearson系数.它非常快.它的速度至少是corrcoef()单独使用的1.5倍,因为它不会冗余地计算系数和其他一些东西.它也会更快,并且不会给你50,000行矩阵的内存问题,因为那时你可以选择存储每组r或者在生成另一组之前处理它们.如果没有存储任何r的长期,我能够在不太新的笔记本电脑上在不到一分钟的时间内运行上面的代码来运行50,000 x 10组随机生成的数据.

旧解决方案

首先,我不建议将r打印到屏幕上.对于100行(10列),这与打印相差19.79秒而不使用代码则为0.301秒.只需存储r并在以后根据需要使用它们,或者在进行时对它们进行一些处理,就像寻找一些最大的r一样.

其次,通过不冗余地计算一些数量,可以节省一些成本.Pearson系数是用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)

当我删除了p值的东西时,我的直接scipy代码加速了大约4.8倍 - 如果我把p值的东西留在那里,我得到8.8x(我使用了10行,有数百行).我还检查过它确实给出了相同的结果.这不是一个真正巨大的改进,但它可能会有所帮助.

最终,你遇到了你正在计算的问题(50000)*(50001)/ 2 = 1,250,025,000 Pearson系数(如果我正确计算).好多啊.顺便说一下,实际上没有必要计算每一行的Pearson系数(它将等于1),但这只能使您免于计算50,000 Pearson系数.使用上面的代码,如果根据我在较小数据集上的结果有10列数据,我预计计算大约需要4 1/4小时.

通过将上面的代码带入Cython或类似的东西,你可以得到一些改进.如果你幸运的话,我希望你可以比直接的Scipy提高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