原问题
我将大小为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上找到。
2条答案
按热度按时间ghhaqwfi1#
我们可以将
np.einsum
引入其中;但是,您的里程可能是vary,这取决于您的编译以及它是否使用了SSE 2。额外的einsum
调用来替换求和操作可能看起来无关紧要,但是numpy ufunc直到numpy 1.8才使用SSE 2,而einsum
使用,我们可以避免一些if
语句。在opteron内核上运行这个命令,我得到了一个奇怪的结果,因为我预计
dot
调用会占用大部分时间。运行这个与英特尔系统与英特尔mkl再次我得到一些更合理/预期:
在英特尔的机器上,还有一个更大的东西:
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倍: