优化平均外部产品

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加速问题中列出的输入！