寻找相关性矩阵

Question

我有一个比较大的矩阵（大约有5万行），我想打印出矩阵中每一行之间的相关系数。我写了这样的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模块的pearsonr函数（http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.pearsonr.html）。

我的问题是：有没有更快的方法来做到这一点？有没有什么矩阵分割的技巧可以使用？

谢谢！

Answer 1

你可以使用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

Answer 2

你有没有试过直接用 numpy.corrcoef 呢？因为你并不需要 p 值，它应该能完全满足你的需求，而且操作起来简单方便。（除非我对皮尔逊相关系数的理解有误，这种可能性是有的。）

我快速检查了一下随机数据的结果，发现它和 @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倍，因为它避免了重复计算系数和其他一些操作。对于一个有50,000行的矩阵，它也会更快，并且不会出现内存问题，因为你可以选择存储每一组r值，或者在生成另一组之前先处理它们。在不长期存储任何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)

当我去掉p值的部分时，我的速度比直接使用scipy的代码快了大约4.8倍，如果保留p值的部分，速度提升为8.8倍（我使用了10列和几百行的数据）。我还检查了结果，发现它们是相同的。这虽然不是一个巨大的提升，但可能会有所帮助。

最终，你面临的问题是你需要计算(50000)*(50001)/2 = 1,250,025,000个皮尔逊系数（如果我没算错的话）。这可真不少。顺便说一下，实际上没有必要计算每一行与自身的皮尔逊系数（它会等于1），这样只会节省计算50,000个皮尔逊系数。根据我在较小数据集上的结果，我预计使用上面的代码计算10列数据大约需要4个小时15分钟。

通过将上述代码转到Cython或类似的工具，你可以获得一些改进。如果运气好的话，可能会比直接使用Scipy快10倍。此外，正如pyInTheSky所建议的，你还可以进行一些多进程处理。