Numpy "vectorisé" produit scalaire rangée par rangée en cours d'exécution plus lente qu'une boucle

Question

Etant donné une matrice A avec la forme (n, k) et un vecteur s de taille n, je souhaite calculer un matrice G avec la forme (k, k) comme suit:

G + = s [i] * A [i] .T * A [i], pour tout i dans { 0, ..., n-1}

I a essayé de mettre en oeuvre que l'utilisation d'une boucle (Méthode 1) et d'une manière vectorisé (Méthode 2), mais la mise en oeuvre boucle est plus rapide pour les grandes valeurs de k (spécialement quand k > 500).

Le code a été écrit comme suit:

import numpy as np 
k = 200 
n = 50000 
A = np.random.randint(0, 1000, (n,k)) # generates random data for the matrix A (n,k) 
G1 = np.zeros((k,k)) # initialize G1 as a (k,k) matrix 
s = np.random.randint(0, 1000, n) * 1.0 # initialize a random vector of size n 

# METHOD 1 
for i in xrange(n): 
    G1 += s[i] * np.dot(np.array([A[i]]).T, np.array([A[i]])) 

# METHOD 2 
G2 = np.dot(A[:,np.newaxis].T, s[:,np.newaxis]*A) 
G2 = np.squeeze(G2) # reduces dimension from (k,1,k) to (k,k) 

Les matrices G1 et G2 sont identiques (ils sont la matrice G), et la seule différence est la façon dont ils ont été calculés. Y a-t-il un moyen plus intelligent et efficace de calculer cela?

Enfin, ce sont les temps que j'ai obtenu avec des tailles aléatoires pour k et n:

Test #: 1 
k,n: (866, 45761) 
Method1: 337.457569838s 
Method2: 386.290487051s 
-------------------- 
Test #: 2 
k,n: (690, 48011) 
Method1: 152.999140978s 
Method2: 226.080267191s 
-------------------- 
Test #: 3 
k,n: (390, 5317) 
Method1: 5.28722500801s 
Method2: 4.86999702454s 
-------------------- 
Test #: 4 
k,n: (222, 5009) 
Method1: 1.73456382751s 
Method2: 0.929286956787s 
-------------------- 
Test #: 5 
k,n: (915, 16561) 
Method1: 101.782826185s 
Method2: 159.167108059s 
-------------------- 
Test #: 6 
k,n: (130, 11283) 
Method1: 1.53138184547s 
Method2: 0.64450097084s 
-------------------- 
Test #: 7 
k,n: (57, 37863) 
Method1: 1.44776391983s 
Method2: 0.494270086288s 
-------------------- 
Test #: 8 
k,n: (110, 34599) 
Method1: 3.51851701736s 
Method2: 1.61688089371s 

Answer 1

Deux versions beaucoup plus améliorée seraient -

(A.T*s).dot(A) 
(A.T).dot(A*s[:,None]) 

Problème (s) avec method2:

Avec method2, nous sommes en train de créer A[:,np.newaxis].T, qui serait de forme (k,1,n), c'est un tableau 3D. Je pense avec un tableau 3D, np.dot va dans une sorte de boucle et n'est pas vraiment vectorisé (le code source pourrait révéler plus d'informations ici).

Pour de telles multiplications de tenseur 3D, il est préférable d'utiliser l'équivalent tenseur: np.tensordot. Ainsi, une version améliorée de method2 devient:

G2 = np.tensordot(A[:,np.newaxis].T, s[:,np.newaxis]*A, axes=((2),(0))) 
G2 = np.squeeze(G2) 

Depuis, nous sommes sum-reducing juste un axe de chacune de ces entrées avec np.tensordot, nous ne avons pas vraiment besoin tensordot ici et simplement np.dot sur la version squeezed-in suffirait. Cela nous ramènera à .

Runtime tests

approches -

def method1(A, s): 
    G1 = np.zeros((k,k)) # initialize G1 as a (k,k) matrix 
    for i in xrange(n): 
     G1 += s[i] * np.dot(np.array([A[i]]).T, np.array([A[i]])) 
    return G1 

def method2(A, s): 
    G2 = np.dot(A[:,np.newaxis].T, s[:,np.newaxis]*A) 
    G2 = np.squeeze(G2) # reduces dimension from (k,1,k) to (k,k) 
    return G2 

def method3(A, s): 
    return (A.T*s).dot(A) 

def method4(A, s): 
    return (A.T).dot(A*s[:,None]) 

def method2_improved(A, s): 
    G2 = np.tensordot(A[:,np.newaxis].T, s[:,np.newaxis]*A, axes=((2),(0))) 
    G2 = np.squeeze(G2) 
    return G2 

minutage et vérification -

In [56]: k = 200 
    ...: n = 5000 
    ...: A = np.random.randint(0, 1000, (n,k)) 
    ...: s = np.random.randint(0, 1000, n) * 1.0 
    ...: 

In [72]: print np.allclose(method1(A, s), method2(A, s)) 
    ...: print np.allclose(method1(A, s), method3(A, s)) 
    ...: print np.allclose(method1(A, s), method4(A, s)) 
    ...: print np.allclose(method1(A, s), method2_improved(A, s)) 
    ...: 
True 
True 
True 
True 

In [73]: %timeit method1(A, s) 
    ...: %timeit method2(A, s) 
    ...: %timeit method3(A, s) 
    ...: %timeit method4(A, s) 
    ...: %timeit method2_improved(A, s) 
    ...: 
1 loops, best of 3: 1.12 s per loop 
1 loops, best of 3: 693 ms per loop 
100 loops, best of 3: 8.12 ms per loop 
100 loops, best of 3: 8.17 ms per loop 
100 loops, best of 3: 8.28 ms per loop