高效縱列相關係數計算

我大小爲n的P列相關，以大小爲n×米的矩陣的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得到了Intel Core時序：

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％會很好。

UPDATE 2017年5月

後相當長的一段時間我回到這個問題，重新運行，延長了任務和測試。

使用einsum，我已經將代碼擴展到P不是行而是矩陣的情況。因此，任務是將O的所有列與P的所有列相關聯。對於如何用通常用於科學計算的不同語言解決同一問題感到好奇，我在其他人的幫助下實施了它（在其他人的幫助下） MATLAB，Julia和R. MATLAB和Julia是最快的，他們有專門的例程來計算列向相關性。 R還有一個專門的例程，但是最慢。
在numpy的當前版本（1.12.1來自Anaconda）中，einsum仍勝過我使用的專用函數。

所有的腳本和計時都在https://github.com/ikizhvatov/efficient-columnwise-correlation。

來源

2013-10-16 ilya

可以改變'O'和'P'嗎？而且這些方法在很大程度上依賴於鏈接到numpy的BLAS庫，如果你還沒有，我會建議使用某種優化的BLAS。 – Daniel

@Ophion，'P'可以改變，'O'不能。我使用從http://www.lfd.uci.edu/~gohlke/pythonlibs/鏈接到英特爾MKL的numpy。 – ilya

目前無法測試，但如果您嘗試使用np.einsum（'ij，ij-> j'，DO，DO）'而不是'np.sum（DO ** 2， 0）'和'np.dot（DP，DP）'而不是'np.sum（DP ** 2）'。 – Jaime

我們可以給這個引入np.einsum;但是，根據您的彙編以及是否使用SSE2，您的milage可能爲vary。額外的einsum調用替換求和操作可能看起來是無關的，但numpy ufuncs不會使用SSE2直到numpy 1.8，而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

與英特爾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

來源

2013-10-16 13:25:12 Daniel

是'tmp = ...; tmp * = ...'應該節省時間嗎？我已經看到其他人也這樣做，我沒有看到'tmp = ... * ...'的優勢。 – askewchan

唯一的原因是前者對我更具可讀性。 – Daniel

高效縱列相關係數計算

回答

相關問題