在Python中的阵列中的每个元素快速计算特征向量

Question

我要计算特征向量数据的阵列（在我的实际情况，我云多边形）

为此我写了这个功能：

import numpy as np 

def eigen(data): 
    eigenvectors = [] 
    eigenvalues = [] 

    for d in data: 
     # compute covariance for each triangle 
     cov = np.cov(d, ddof=0, rowvar=False) 

     # compute eigen vectors 
     vals, vecs = np.linalg.eig(cov) 
     eigenvalues.append(vals) 
     eigenvectors.append(vecs) 

    return np.array(eigenvalues), np.array(eigenvectors) 

的一些测试数据运行此：

import cProfile 
triangles = np.random.random((10**4,3,3,)) # 10k 3D triangles 
cProfile.run('eigen(triangles)') # 550005 function calls in 0.933 seconds 

工作正常，但它会因为迭代循环的速度很慢。有没有更快的方法来计算我所需要的数据而无需迭代数组？如果没有，任何人都可以建议加快速度的方法吗？

Answer 1

本事！

嗯，我攻入covariance func definition，并把在规定输入状态：ddof=0, rowvar=False和事实证明，一切都降低了，只是三线 -

nC = m.shape[1] # m is the 2D input array 
X = m - m.mean(0) 
out = np.dot(X.T, X)/nC 

把它扩大到我们的3D阵列的情况下，我写的向下多圈版本与这三条线被重复用于2D阵列的部分从3D输入阵列，像这样 -

for i,d in enumerate(m): 

    # Using np.cov : 
    org_cov = np.cov(d, ddof=0, rowvar=False) 

    # Using earlier 2D array hacked version : 
    nC = m[i].shape[0] 
    X = m[i] - m[i].mean(0,keepdims=True) 
    hacked_cov = np.dot(X.T, X)/nC 

升压它向上

我们需要加快最后三行。在所有迭代X计算可以用做broadcasting -

diffs = data - data.mean(1,keepdims=True) 

接下来，点积运算的各个版本可以与transpose和np.dot来完成，但transpose可能是这样的多昂贵的事情维数组。存在更好的选择在np.einsum，像这样 -

cov3D = np.einsum('ijk,ijl->ikl',diffs,diffs)/data.shape[1] 

使用它！

综上所述：

for d in data: 
    # compute covariance for each triangle 
    cov = np.cov(d, ddof=0, rowvar=False) 

可以预先计算，如下所示：

diffs = data - data.mean(1,keepdims=True) 
cov3D = np.einsum('ijk,ijl->ikl',diffs,diffs)/data.shape[1] 

这些预计算值可以跨迭代被用来计算像这样的本征矢量 -

for i,d in enumerate(data): 
    # Directly use pre-computed covariances for each triangle 
    vals, vecs = np.linalg.eig(cov3D[i]) 

测试它！

这里有一些运行测试，以评估的前期计算的协方差结果的影响 -

In [148]: def original_app(data): 
    ...:  cov = np.empty(data.shape) 
    ...:  for i,d in enumerate(data):  
    ...:   # compute covariance for each triangle 
    ...:   cov[i] = np.cov(d, ddof=0, rowvar=False) 
    ...:  return cov 
    ...: 
    ...: def vectorized_app(data):    
    ...:  diffs = data - data.mean(1,keepdims=True) 
    ...:  return np.einsum('ijk,ijl->ikl',diffs,diffs)/data.shape[1] 
    ...: 

In [149]: data = np.random.randint(0,10,(1000,3,3)) 

In [150]: np.allclose(original_app(data),vectorized_app(data)) 
Out[150]: True 

In [151]: %timeit original_app(data) 
10 loops, best of 3: 64.4 ms per loop 

In [152]: %timeit vectorized_app(data) 
1000 loops, best of 3: 1.14 ms per loop 

In [153]: data = np.random.randint(0,10,(5000,3,3)) 

In [154]: np.allclose(original_app(data),vectorized_app(data)) 
Out[154]: True 

In [155]: %timeit original_app(data) 
1 loops, best of 3: 324 ms per loop 

In [156]: %timeit vectorized_app(data) 
100 loops, best of 3: 5.67 ms per loop 

Answer 2

我不知道你实际能达到多少加速。

这是一个轻微的修改，可以帮助一点：

%timeit -n 10 values, vectors = \ 
    eigen(triangles) 
10 loops, best of 3: 745 ms per loop 

%timeit values, vectors = \ 
    zip(*(np.linalg.eig(np.cov(d, ddof=0, rowvar=False)) for d in triangles)) 
10 loops, best of 3: 705 ms per loop