我有一个相当大的矩阵（大约 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 函数。
我的问题是：有更快的方法吗？我可以使用一些矩阵分区技术吗？
谢谢！
【问题讨论】：
新解决方案
看了 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 系数（它将等于 1），但这只会使您免于计算 50,000 个 Pearson 系数。使用上面的代码，根据我在较小数据集上的结果，如果您的数据有 10 列，我预计大约需要 4 1/4 小时来完成计算。
您可以通过将上述代码导入 Cython 或类似的东西来获得一些改进。我希望如果幸运的话，你可能会比直接 Scipy 提高 10 倍。此外，根据 pyInTheSky 的建议，您可以进行一些多处理。
【讨论】：
您是否尝试过仅使用numpy.corrcoef？看看你如何不使用 p 值，它应该完全按照你的意愿做，尽可能少做文章。 （除非我记错了 pearson 的 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 应该做你想做的事，而且速度要快得多。
【讨论】：
corrcoef，但出于某种原因，我记得它变慢了。感觉有点不好意思，因为我相信我的糟糕记忆而不是尝试它。它更快，因为它使用矩阵乘法来消除 python 循环。来自我的 +1。
你可以使用 python 多进程模块，将你的行分成 10 组，缓冲你的结果，然后打印出来（但这只会在多核机器上加速）
http://docs.python.org/library/multiprocessing.html
顺便说一句：您还必须将您的 sn-p 变成一个函数，并考虑如何进行数据重组。让每个子进程都有一个这样的列表 ...[startcord,stopcord,buff] .. 可能会很好地工作
def myfunc(thelist):
for i in xrange(thelist[0]:thelist[1]):
....
thelist[2] = result
【讨论】：