Vectorisation covariance NumPy pour le tableau 3D

Question

J'ai un tableau numpy 3D de forme (t,n1,n2):

x=np.random.rand(10,2,4) 

je dois calculer un autre tableau 3Dy qui est de forme (t,n1,n1) tel que:

y[0] = np.cov[x[0,:,:]) et ainsi de suite pour toutes les tranches le long du premier axe.

Ainsi, une mise en œuvre loufoque serait -

y=np.zeros((10,2,2)) 
for i in np.arange(x.shape[0]): 
    y[i]=np.cov(x[i,:,:]) 

Est-il possible de vectoriser ce que je peux calculer toutes les matrices de covariance en une seule fois? J'ai essayé de faire:

x1= x.swapaxes(1,2) 
y= np.dot(x,x1) 

Mais cela n'a pas fonctionné.

Answer 1

Piraté dans numpy.cov source code et essayé en utilisant les paramètres par défaut. Il se trouve que, np.cov(x[i,:,:]) serait simplement:

N = x.shape[2] 
m = x[i,:,:] 
m -= np.sum(m, axis=1, keepdims=True)/N 
cov = np.dot(m, m.T) /(N - 1) 

Ainsi, la tâche était de vectoriser cette boucle qui itérer i et traiter toutes les données de x en une seule fois. Pour le même, nous pourrions utiliser broadcasting à la troisième étape. Pour la dernière étape, nous exécutons sum-reduction le long de toutes les tranches du premier axe. Cela pourrait être efficacement mis en œuvre de manière vectorisée avec np.einsum. Ainsi, la mise en œuvre finale est venu à cette -

N = x.shape[2] 
m1 = x - x.sum(2,keepdims=1)/N 
y_out = np.einsum('ijk,ilk->ijl',m1,m1) /(N - 1) 

Runtime Test

In [155]: def original_app(x): 
    ...:  n = x.shape[0] 
    ...:  y = np.zeros((n,2,2)) 
    ...:  for i in np.arange(x.shape[0]): 
    ...:   y[i]=np.cov(x[i,:,:]) 
    ...:  return y 
    ...: 
    ...: def proposed_app(x): 
    ...:  N = x.shape[2] 
    ...:  m1 = x - x.sum(2,keepdims=1)/N 
    ...:  out = np.einsum('ijk,ilk->ijl',m1,m1)/(N - 1) 
    ...:  return out 
    ...: 

In [156]: # Setup inputs 
    ...: n = 10000 
    ...: x = np.random.rand(n,2,4) 
    ...: 

In [157]: np.allclose(original_app(x),proposed_app(x)) 
Out[157]: True # Results verified 

In [158]: %timeit original_app(x) 
1 loops, best of 3: 610 ms per loop 

In [159]: %timeit proposed_app(x) 
100 loops, best of 3: 6.32 ms per loop 

énorme speedup là!