原始问题
我将大小为 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 中实现了它(在其他人的帮助下)。MATLAB 和 Julia 是最快的,而且它们有计算列相关性的专用例程。 R 也有专用例程,但速度最慢。
在当前版本的 numpy(来自 Anaconda 的 1.12.1)中,einsum 仍然胜过我使用的专用函数。
所有脚本和时间安排均可在 https://github.com/ikizhvatov/efficient-columnwise-correlation 获得。
我们可以对此进行介绍
np.einsum
;但是,您的里程可能会“有所不同”,具体取决于您的编译以及是否使用 SSE2。额外的 einsum
调用来替换求和操作可能看起来无关紧要,但是 numpy ufunc 直到 numpy 1.8 才使用 SSE2,而 einsum
则使用 SSE2,我们可以避免一些 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
高效的解决方案,但它具有高度可读性并且只有一行代码:
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)