numpy 高效的列相关系数计算

svujldwt 于 2023-10-19 发布在其他

关注(0)|答案(2)|浏览(178)

原问题

我将大小为n的行P与大小为n×m的矩阵O的每一列相关联。我编写了以下代码：

import numpy as np

def ColumnWiseCorrcoef(O, P):
    n = P.size
    DO = O - (np.sum(O, 0) / np.double(n))
    DP = P - (np.sum(P) / np.double(n))
    return np.dot(DP, DO) / np.sqrt(np.sum(DO ** 2, 0) * np.sum(DP ** 2))

它比简单的方法更有效：

def ColumnWiseCorrcoefNaive(O, P):
    return np.corrcoef(P,O.T)[0,1:O[0].size+1]

以下是我在英特尔内核上使用numpy-1.7.1-MKL获得的时间：

O = np.reshape(np.random.rand(100000), (1000,100))
P = np.random.rand(1000)

%timeit -n 1000 A = ColumnWiseCorrcoef(O, P)
1000 loops, best of 3: 787 us per loop
%timeit -n 1000 B = ColumnWiseCorrcoefNaive(O, P)
1000 loops, best of 3: 2.96 ms per loop

现在的问题是：你能建议一个更快的版本的代码来解决这个问题吗？如果能多挤出20%就好了。

更新2017年5月

过了一段时间，我又回到了这个问题，重新运行并扩展了任务和测试。
1.使用einsum，我已经将代码扩展到P不是行而是矩阵的情况。所以任务是将O的所有列与P的所有列相关联。
1.出于对如何用科学计算中常用的不同语言解决同一问题的好奇，我用MATLAB、Julia和R实现了它（在其他人的帮助下）。MATLAB和Julia是最快的，它们有一个专用的例程来计算列相关性。R也有一个专用的例程，但它是最慢的。
1.在当前版本的numpy（来自Anaconda的1.12.1）中，einsum仍然胜过我使用的专用函数。
所有的脚本和计时都可以在https://github.com/ikizhvatov/efficient-columnwise-correlation上找到。

numpy

来源：https://stackoverflow.com/questions/19401078/efficient-columnwise-correlation-coefficient-calculation

2条答案

按热度按时间

ghhaqwfi1#

我们可以将np.einsum引入其中;但是，您的里程可能是vary，这取决于您的编译以及它是否使用了SSE 2。额外的einsum调用来替换求和操作可能看起来无关紧要，但是numpy ufunc直到numpy 1.8才使用SSE 2，而einsum使用，我们可以避免一些if语句。
在opteron内核上运行这个命令，我得到了一个奇怪的结果，因为我预计dot调用会占用大部分时间。

def newColumnWiseCorrcoef(O, P):
    n = P.size
    DO = O - (np.einsum('ij->j',O) / np.double(n))
    P -= (np.einsum('i->',P) / np.double(n))
    tmp = np.einsum('ij,ij->j',DO,DO)
    tmp *= np.einsum('i,i->',P,P)          #Dot or vdot doesnt really change much.
    return np.dot(P, DO) / np.sqrt(tmp)

O = np.reshape(np.random.rand(100000), (1000,100))
P = np.random.rand(1000)

old = ColumnWiseCorrcoef(O,P)
new = newColumnWiseCorrcoef(O,P)

np.allclose(old,new)
True

%timeit ColumnWiseCorrcoef(O,P)
100 loops, best of 3: 1.52ms per loop

%timeit newColumnWiseCorrcoef(O,P)
1000 loops, best of 3: 518us per loop

运行这个与英特尔系统与英特尔mkl再次我得到一些更合理/预期：

%timeit ColumnWiseCorrcoef(O,P)
1000 loops, best of 3: 524 µs per loop

%timeit newColumnWiseCorrcoef(O,P)
1000 loops, best of 3: 354 µs per loop

在英特尔的机器上，还有一个更大的东西：

O = np.random.rand(1E5,1E3)
P = np.random.rand(1E5)

%timeit ColumnWiseCorrcoef(O,P)
1 loops, best of 3: 1.33 s per loop

%timeit newColumnWiseCorrcoef(O,P)
1 loops, best of 3: 791 ms per loop

赞(0）回复(0）举报 2023-10-19

juzqafwq2#

这是一个效率较低的解决方案，但它是高度可读的，只有一行代码：
np.einsum("ij,i->j", spst.zscore(A), spst.zscore(B)) / A.shape[0]
而且，如果希望B是矩阵而不是向量，并且要计算相应列对的相关性，则很容易修改为：
np.einsum("ij,ij->j", spst.zscore(A), spst.zscore(B)) / A.shape[0]
根据我在2.8 GHz四核Intel Core i5 iMac上的计时测试，它大约慢3倍：

import scipy.stats as spst

def ColumnWiseCorrcoef(O, P):
    n = P.size
    DO = O - (np.sum(O, 0) / np.double(n))
    DP = P - (np.sum(P) / np.double(n))
    return np.dot(DP, DO) / np.sqrt(np.sum(DO ** 2, 0) * np.sum(DP ** 2))

def newColumnWiseCorrcoef(O, P):
    n = P.size
    DO = O - (np.einsum('ij->j',O) / np.double(n))
    P -= (np.einsum('i->',P) / np.double(n))
    tmp = np.einsum('ij,ij->j',DO,DO)
    tmp *= np.einsum('i,i->',P,P)          #Dot or vdot doesnt really change much.
    return np.dot(P, DO) / np.sqrt(tmp)

def OneLineColumnWiseCorrcoef(A, B):
    return np.einsum("ij,i->j", spst.zscore(A), spst.zscore(B)) / A.shape[0]

O = np.reshape(np.random.rand(100000), (1000,100))
P = np.random.rand(1000)

old = ColumnWiseCorrcoef(O,P)
new = OneLineColumnWiseCorrcoef(O, P)
np.allclose(old,new)
True

%timeit ColumnWiseCorrcoef(O,P)
# 325 µs ± 19.8 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

%timeit newColumnWiseCorrcoef(O,P)
# 274 µs ± 713 ns per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

%timeit OneLineColumnWiseCorrcoef(O,P)
# 1.02 ms ± 4.4 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

赞(0）回复(0）举报 2023-10-19

我来回答

numpy 高效的列相关系数计算

2条答案

相关问题

热门标签

最新问答