Pourquoi le bouclage bat-il l'indexation ici?

Question

Il ya quelques années, quelqu'un posted sur Recettes d'état actif à des fins de comparaison, trois fonctions python/NumPy; chacun de ceux-ci a accepté les mêmes arguments et a renvoyé le même résultat, une matrice de distance .

Deux d'entre eux ont été extraits de sources publiées; ils sont tous les deux - ou ils me semblent être - un code numérique idiomatique. Les calculs répétitifs requis pour créer une matrice de distance sont pilotés par la syntaxe d'index élégante de numpy. Voici l'un d'eux:

from numpy.matlib import repmat, repeat 

def calcDistanceMatrixFastEuclidean(points): 
    numPoints = len(points) 
    distMat = sqrt(sum((repmat(points, numPoints, 1) - 
      repeat(points, numPoints, axis=0))**2, axis=1)) 
    return distMat.reshape((numPoints,numPoints)) 

Le troisième créé la matrice à distance en utilisant une seule boucle (ce qui est évidemment beaucoup de looping étant donné qu'une matrice de distance de seulement 1000 points 2D, a un million d'entrées). À première vue, cette fonction me semblait être le code que j'écrivais quand j'apprenais NumPy et j'écrivais du code NumPy en écrivant d'abord du code Python, puis en le traduisant ligne par ligne.

Plusieurs mois après la publication de l'état actif, les résultats des tests de performance comparant les trois ont été postés et discutés dans un thread sur la liste de diffusion NumPy.

La fonction avec la boucle en fait significativement surperformé les deux autres:

from numpy import mat, zeros, newaxis 

def calcDistanceMatrixFastEuclidean2(nDimPoints): 
    nDimPoints = array(nDimPoints) 
    n,m = nDimPoints.shape 
    delta = zeros((n,n),'d') 
    for d in xrange(m): 
    data = nDimPoints[:,d] 
    delta += (data - data[:,newaxis])**2 
    return sqrt(delta) 

Un participant dans le fil (Keir Mierle) a offert une raison pour laquelle cela pourrait être vrai:

La raison pour laquelle je suspecte que ce sera plus rapide est qu'il a une meilleure localité, en terminant complètement un calcul sur un ensemble de travail relativement petit avant de passer à la prochaine un. Les liners doivent tirer le tableau MxN potentiellement grand dans le processeur à plusieurs reprises.

Par son propre compte cette affiche, sa remarque est qu'un soupçon, et il ne semble pas qu'il a été discuté plus loin.

Avez-vous d'autres idées sur la façon de comptabiliser ces résultats?

En particulier, y a-t-il une règle utile - en ce qui concerne quand boucler et quand indexer - qui peut être extraite de cet exemple comme guide dans l'écriture de code numpy?

Pour ceux qui ne connaissent pas NumPy, ou qui n'ont pas regardé le code, cette comparaison n'est pas basée sur un cas limite - ce ne serait certainement pas si intéressant pour moi si c'était le cas. Au lieu de cela, cette comparaison implique une fonction qui effectue une tâche commune dans le calcul matriciel (c.-à-d., Création d'un tableau de résultats avec deux antécédents); de plus, chaque fonction est à son tour composée des inserts numpy les plus courants.

Answer 1

TL; DR Le deuxième code ci-dessus ne fait que boucler sur le nombre de dimensions des points (3 fois à travers la boucle for pour les points 3D), donc la boucle n'est pas beaucoup là. L'accélération réelle dans le second code ci-dessus est qu'il exploite mieux la puissance de Numpy pour éviter de créer des matrices supplémentaires lors de la recherche des différences entre les points. Cela réduit la mémoire utilisée et l'effort de calcul.

plus long Explication Je pense que la fonction calcDistanceMatrixFastEuclidean2 vous trompe avec sa boucle peut-être. Il ne fait que boucler sur le nombre de dimensions des points. Pour les points 1D, la boucle ne s'exécute qu'une fois, pour 2D, deux fois, et pour 3D, trois fois. Ce n'est vraiment pas beaucoup en boucle.

Analysons un peu le code pour voir pourquoi l'un est plus rapide que l'autre. calcDistanceMatrixFastEuclidean Je vais appeler fast1 et calcDistanceMatrixFastEuclidean2 sera fast2.

fast1 est basé sur la façon de faire Matlab, comme en témoigne la fonction repmap. La fonction repmap crée dans ce cas un tableau qui est simplement les données d'origine répétées encore et encore. Cependant, si vous regardez le code de la fonction, c'est très inefficace. Il utilise de nombreuses fonctions Numpy (3 reshape s et 2 repeat s) pour ce faire. La fonction repeat est également utilisée pour créer un tableau qui contient les données d'origine avec chaque élément de données répété plusieurs fois. Si nos données d'entrée sont [1,2,3] alors nous soustrayons [1,2,3,1,2,3,1,2,3] de [1,1,1,2,2,2,3,3,3]. Numpy a dû créer beaucoup de matrices supplémentaires entre le code C de Numpy qui aurait pu être évité.

fast2 utilise plus de charges lourdes de Numpy sans créer autant de matrices entre les appels Numpy. fast2 effectue une boucle sur chaque dimension des points, effectue la soustraction et conserve un total cumulé des différences au carré entre chaque dimension. Seulement à la fin est la racine carrée fait. Jusqu'à présent, cela peut ne pas sembler aussi efficace que fast1, mais fast2 évite de faire les choses repmat en utilisant l'indexation de Numpy. Regardons le cas 1D pour plus de simplicité. fast2 crée un tableau 1D des données et les soustrait d'un tableau 2D (N x 1) des données. Cela crée la matrice de différences entre chaque point et tous les autres points sans avoir à utiliser repmat et repeat et ainsi contourner la création de beaucoup de tableaux supplémentaires. C'est là que la différence de vitesse réelle réside à mon avis. fast1 crée beaucoup d'extra entre les matrices (et ils sont créés coûteusement calculés) pour trouver les différences entre les points tandis que fast2 exploite mieux la puissance de Numpy pour les éviter.

D'ailleurs, voici un peu plus rapide version de fast2:

def calcDistanceMatrixFastEuclidean3(nDimPoints): 
    nDimPoints = array(nDimPoints) 
    n,m = nDimPoints.shape 
    data = nDimPoints[:,0] 
    delta = (data - data[:,newaxis])**2 
    for d in xrange(1,m): 
    data = nDimPoints[:,d] 
    delta += (data - data[:,newaxis])**2 
    return sqrt(delta) 

La différence est que nous ne créons plus delta comme une matrice de zéros.

Answer 2

dis pour le plaisir:

dis.dis(calcDistanceMatrixFastEuclidean)

2   0 LOAD_GLOBAL    0 (len) 
       3 LOAD_FAST    0 (points) 
       6 CALL_FUNCTION   1 
       9 STORE_FAST    1 (numPoints) 

    3   12 LOAD_GLOBAL    1 (sqrt) 
      15 LOAD_GLOBAL    2 (sum) 
      18 LOAD_GLOBAL    3 (repmat) 
      21 LOAD_FAST    0 (points) 
      24 LOAD_FAST    1 (numPoints) 
      27 LOAD_CONST    1 (1) 
      30 CALL_FUNCTION   3 

    4   33 LOAD_GLOBAL    4 (repeat) 
      36 LOAD_FAST    0 (points) 
      39 LOAD_FAST    1 (numPoints) 
      42 LOAD_CONST    2 ('axis') 
      45 LOAD_CONST    3 (0) 
      48 CALL_FUNCTION   258 
      51 BINARY_SUBTRACT 
      52 LOAD_CONST    4 (2) 
      55 BINARY_POWER 
      56 LOAD_CONST    2 ('axis') 
      59 LOAD_CONST    1 (1) 
      62 CALL_FUNCTION   257 
      65 CALL_FUNCTION   1 
      68 STORE_FAST    2 (distMat) 

    5   71 LOAD_FAST    2 (distMat) 
      74 LOAD_ATTR    5 (reshape) 
      77 LOAD_FAST    1 (numPoints) 
      80 LOAD_FAST    1 (numPoints) 
      83 BUILD_TUPLE    2 
      86 CALL_FUNCTION   1 
      89 RETURN_VALUE 

dis.dis(calcDistanceMatrixFastEuclidean2)

2   0 LOAD_GLOBAL    0 (array) 
       3 LOAD_FAST    0 (nDimPoints) 
       6 CALL_FUNCTION   1 
       9 STORE_FAST    0 (nDimPoints) 

    3   12 LOAD_FAST    0 (nDimPoints) 
      15 LOAD_ATTR    1 (shape) 
      18 UNPACK_SEQUENCE   2 
      21 STORE_FAST    1 (n) 
      24 STORE_FAST    2 (m) 

    4   27 LOAD_GLOBAL    2 (zeros) 
      30 LOAD_FAST    1 (n) 
      33 LOAD_FAST    1 (n) 
      36 BUILD_TUPLE    2 
      39 LOAD_CONST    1 ('d') 
      42 CALL_FUNCTION   2 
      45 STORE_FAST    3 (delta) 

    5   48 SETUP_LOOP    76 (to 127) 
      51 LOAD_GLOBAL    3 (xrange) 
      54 LOAD_FAST    2 (m) 
      57 CALL_FUNCTION   1 
      60 GET_ITER 
     >> 61 FOR_ITER    62 (to 126) 
      64 STORE_FAST    4 (d) 

    6   67 LOAD_FAST    0 (nDimPoints) 
      70 LOAD_CONST    0 (None) 
      73 LOAD_CONST    0 (None) 
      76 BUILD_SLICE    2 
      79 LOAD_FAST    4 (d) 
      82 BUILD_TUPLE    2 
      85 BINARY_SUBSCR 
      86 STORE_FAST    5 (data) 

    7   89 LOAD_FAST    3 (delta) 
      92 LOAD_FAST    5 (data) 
      95 LOAD_FAST    5 (data) 
      98 LOAD_CONST    0 (None) 
      101 LOAD_CONST    0 (None) 
      104 BUILD_SLICE    2 
      107 LOAD_GLOBAL    4 (newaxis) 
      110 BUILD_TUPLE    2 
      113 BINARY_SUBSCR 
      114 BINARY_SUBTRACT 
      115 LOAD_CONST    2 (2) 
      118 BINARY_POWER 
      119 INPLACE_ADD 
      120 STORE_FAST    3 (delta) 
      123 JUMP_ABSOLUTE   61 
     >> 126 POP_BLOCK 

    8  >> 127 LOAD_GLOBAL    5 (sqrt) 
      130 LOAD_FAST    3 (delta) 
      133 CALL_FUNCTION   1 
      136 RETURN_VALUE 

Je ne suis pas un expert en dis, mais il semble que vous auriez à regarder plus à la f onctions que le premier appelle pour savoir pourquoi ils prennent un certain temps. Il existe également un outil de profil de performance avec Python, cProfile.