優化平均外部產品

Question

我目前正在編寫一個簡短的程序來對隨機矩陣特徵值分佈進行一些分析，但我的分析所需的參數選擇導致整個事情變得非常緩慢。基本上我應該循環下面的函數，理想情況下大約5000次，最後收集完整的特徵值列表。

C = np.zeros((N,N)) 
time_series = np.random.normal(mu,sigma, (N + B*(M-1)) ) 

for k in range(int(M)): 
    C += np.outer(time_series[k*B : (N) + k*B], time_series[k*B : (N) + k*B]) 
C = C/M 

eg_v = np.linalg.eigvalsh(C) 

在我需要N = 1000，B約10，M = 100。 然而，利用這種選擇的參數的程序採用類似4-5小時在我的相當-執行筆記本電腦上運行。

拋開硬件限制，我有辦法在代碼方面加速整個事情嗎？

在此先感謝！

Answer 1

您可以使用np.tensordot

這樣一個量化的解決方案代替循環，如下 -

# Get the starting indices for each iteration 
idx = (np.arange(M)*B)[:,None] + np.arange(N) 

# Get the range of indices across all iterations as a 2D array and index 
# time_series with it to give us "time_series[k*B : (N) + k*B]" equivalent 
time_idx = time_series[idx] 

# Use broadcasting to perform summation accumulation 
C = np.tensordot(time_idx,time_idx,axes=([0],[0])) 

的tensordot可以由被替換 -

C = np.zeros((N,N)) 
for k in range(int(M)): 
    C += np.outer(time_series[k*B : (N) + k*B], time_series[k*B : (N) + k*B]) 

可以被替代簡單的點積：

C = time_idx.T.dot(time_idx) 

個

運行測試

功能：

def original_app(time_series,B,N,M): 
    C = np.zeros((N,N)) 
    for k in range(int(M)): 
     C += np.outer(time_series[k*B : (N) + k*B], time_series[k*B : (N) + k*B]) 
    return C 

def vectorized_app(time_series,B,N,M): 
    idx = (np.arange(M)*B)[:,None] + np.arange(N) 
    time_idx = time_series[idx] 
    return np.tensordot(time_idx,time_idx,axes=([0],[0])) 

輸入：

In [115]: # Inputs 
    ...: mu = 1.2 
    ...: sigma = 0.5 
    ...: N = 1000 
    ...: M = 100 
    ...: B = 10 
    ...: time_series = np.random.normal(mu,sigma, (N + B*(M-1)) ) 
    ...: 

時序：

In [116]: out1 = original_app(time_series,B,N,M) 

In [117]: out2 = vectorized_app(time_series,B,N,M) 

In [118]: np.allclose(out1,out2) 
Out[118]: True 

In [119]: %timeit original_app(time_series,B,N,M) 
1 loops, best of 3: 1.56 s per loop 

In [120]: %timeit vectorized_app(time_series,B,N,M) 
10 loops, best of 3: 26.2 ms per loop 

因此，我們看到一個60x加速問題中列出的輸入！