6

原始问题

我将大小为 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 的所有列相关联。

  2. 对如何用科学计算常用的不同语言解决相同的问题感到好奇，我在 MATLAB、Julia 和 R 中实现了它（在其他人的帮助下）。MATLAB 和 Julia 是最快的，它们有专门的例程计算列相关性。R 也有专门的例程，但速度最慢。

  3. 在当前版本的 numpy（来自 Anaconda 的 1.12.1）中，einsum 仍然胜过我使用的专用函数。

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

4

1 回答 1

5

我们可以介绍np.einsum一下；但是，您的里程可能会因您的编译以及是否使用 SSE2 而有所不同。替换求和运算的额外einsum调用可能看起来无关紧要，但 numpy ufunc 直到 numpy 1.8 才使用 SSE2 而einsum我们可以避免一些if语句。

在带有 intel mkl blas 的 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

再次在英特尔机器上使用更大的东西：

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
于 2013-10-16T13:25:12.047 回答