求相关矩阵

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  

请注意，我正在使用scipy模块（http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.pearsonr.html）中的pearsonr函数。

我的问题是：有没有更快的方法？有什么矩阵分割技术可以用吗？

谢谢！

Answer 1

您可以使用python多进程模块，将行分成10组，缓冲结果，然后打印出来（这只会在多核机器上加快速度）

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

顺便说一句：您还必须将代码片段转换为函数，并考虑如何进行数据重组。让每个子进程都有一个这样的列表…[开始代码，停止线，buff]。。可能很管用

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

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

新解决方案

在查看了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行到最后一行之间的皮尔逊系数。它很快。它至少比单独使用corrcoef()快1.5倍，因为它没有多余地计算系数和其他一些东西。它也会更快，并且不会给您50000行矩阵的内存问题，因为然后您可以选择存储每一组r，或者在生成另一组之前处理它们。在不存储任何r的长期代码的情况下，我能够在不到一分钟的时间内将上述代码运行在我相当新的笔记本电脑上的50000 x 10组随机生成的数据上。

旧解决方案

首先，我不建议把r打印到屏幕上。对于100行（10列），打印时的差异为19.79秒，而不使用代码时的差异为0.301秒。如果你愿意的话，只需储存r，以后再使用，或者在寻找一些最大的r的过程中对它们进行一些处理

第二，不必重复计算一些数量，可以节省一些开支。皮尔逊系数是用scipy计算的，它使用一些可以预先计算的量，而不是每次使用行时都计算。另外，您没有使用p值（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.8倍（我使用了10列几百行）。我还检查了它是否给出了相同的结果。这并不是一个很大的进步，但可能会有所帮助。

最终，你会陷入计算（50000）*（50001）/2=1250025000皮尔逊系数的问题（如果我计算正确的话）。太多了。顺便说一下，实际上不需要计算每一行的Pearson系数（它等于1），但这只会节省您计算50000个Pearson系数的时间。使用上面的代码，如果根据我在较小数据集上的结果，您的数据有10列，我预计大约需要4 1/4小时来完成计算。

通过将上述代码转换成Cython或类似的代码，您可以得到一些改进。如果你运气好的话，我想你可能会比直筒子进步10倍。另外，正如pyinthsky所建议的，您可以进行一些多处理。