performance - 高效的列相关系数计算

原始问题

我正在将大小为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]

以下是在Intel内核上使用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月

一段时间后，我回到了这个问题，重新运行并扩展了任务和测试。

使用einsum，我将代码扩展到了P不是行而是矩阵的情况。因此，任务是将O的所有列与P的所有列相关联。

很好奇如何用科学计算中常用的不同语言解决相同的问题，我在他人，他人的帮助下在MATLAB，Julia和R中实现了该问题。计算列相关的例程。 R也有专用的例程，但是最慢。

在当前版本的numpy(来自Anaconda的1.12.1)中，einsum仍胜过我使用的专用功能。

所有脚本和时间都可以在https://github.com/ikizhvatov/efficient-columnwise-correlation上获得。

最佳答案

我们可以为此引入np.einsum；但是，根据您的编译以及是否使用SSE2，您的里程数可能为vary。用来替换求和操作的额外einsum调用看起来似乎是多余的，但是numpy ufuncs直到numpy 1.8才使用SSE2，而einsum则使用了SSE2，我们可以避免使用一些if语句。

在具有intel mkl bla的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

再次使用具有intel mkl的intel系统运行此命令，我会得到一些更合理/更期望的结果:

%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

在intel机器上再次使用更大的东西:

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

关于performance - 高效的列相关系数计算，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/19401078/