python - 寻找相关矩阵

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

新解决方案

在看了 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.8 倍(我使用了 10 列和数百行).我还检查了它是否给出了相同的结果。这并不是一个真正巨大的改进，但它可能会有所帮助。

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

将上述代码放入Cython或类似软件中可以得到一些改进。如果幸运的话，我预计你可能会比直接 Scipy 提高 10 倍。此外，正如 pyInTheSky 所建议的，您可以进行一些多处理。