Pourquoi numpy.absolute() est-il si lent ?

Question

Je dois optimiser un script qui fait un usage intensif du calcul de la norme L1 des vecteurs. Comme nous le savons, la norme L1 dans ce cas est juste une somme de valeurs absolues. En chronométrant la vitesse de numpy dans cette tâche, j'ai trouvé quelque chose de bizarre : l'addition de tous les éléments du vecteur est environ 3 fois plus rapide que la prise de la valeur absolue de chaque élément du vecteur. C'est un résultat surprenant, car l'addition est assez complexe par rapport à la prise de la valeur absolue, qui nécessite seulement de mettre à zéro chaque 32ème bit d'un bloc de données (en supposant float32).

Pourquoi l'addition est-elle 3x plus rapide qu'une simple opération par bit ?

import numpy as np

a = np.random.rand(10000000)

%timeit np.sum(a)
13.9 ms ± 87.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

%timeit np.abs(a)
41.2 ms ± 92.3 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

Answer 1

Il y a plusieurs choses à considérer ici. sum renvoie un scalaire abs renvoie un tableau. Donc, même si l'addition de deux nombres et la prise de l'absolu avaient la même vitesse abs serait plus lent car il doit créer le tableau. Et il doit traiter deux fois plus d'éléments (lecture depuis l'entrée + écriture vers la sortie).

Vous ne pouvez donc pas déduire de ces timings quoi que ce soit sur la vitesse de l'addition par rapport aux opérations par bit.

Vous pouvez cependant vérifier s'il est plus rapide d'ajouter quelque chose à chaque valeur d'un tableau que de prendre la valeur absolue de chaque valeur.

%timeit a + 0.1
9 ms ± 155 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit abs(a)
9.98 ms ± 532 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Ou encore, comparez la somme + l'allocation de mémoire par rapport à la prise en compte de la valeur absolue.

%timeit np.full_like(a, 1); np.sum(a)
13.4 ms ± 358 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit abs(a) 
9.64 ms ± 320 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Au cas où vous voudriez accélérer le calcul de la norme, vous pourriez essayer numba (ou Cython, ou écrire vous-même une routine C ou Fortran), de cette façon vous évitez toute allocation de mémoire :

import numba as nb

@nb.njit
def sum_of_abs(arr):
    sum_ = 0
    for item in arr:
        sum_ += abs(item)
    return sum_

sum_of_abs(a)  # one call for the jitter to kick in
%timeit sum_of_abs(a)
# 2.44 ms ± 315 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Answer 2

np.sum renvoie un scalaire. np.abs renvoie un nouveau tableau de la même taille. L'allocation de mémoire pour ce nouveau tableau est ce qui prend le plus de temps ici. Comparez

>>> timeit("np.abs(a)", "import numpy as np; a = np.random.rand(10000000)", number=100)
3.565487278989167
>>> timeit("np.abs(a, out=a)", "import numpy as np; a = np.random.rand(10000000)", number=100)
0.9392949139873963

L'argument out=a indique à NumPy de mettre le résultat dans le même tableau a en écrasant les anciennes données qui s'y trouvent. D'où l'accélération.

Sum est toujours légèrement plus rapide :

>>> timeit("np.sum(a)", "import numpy as np; a = np.random.rand(10000000)", number=100)
0.6874654769926565

mais elle ne nécessite pas autant de écrire accès à la mémoire.

Si vous ne voulez pas écraser a, fournir un autre tableau pour la sortie de abs est une possibilité, en supposant que vous devez prendre de façon répétée abs de tableaux de même type et de même taille.

b = np.empty_like(a)   # done once, outside the loop
np.abs(a, out=b)
np.sum(b)

fonctionne en environ la moitié du temps de np.linalg(a, 1)

Pour référence, np.linalg calcule la norme L1 comme

add.reduce(abs(x), axis=axis, keepdims=keepdims)

ce qui implique l'allocation de mémoire pour le nouveau tableau abs(x) .

---

Idéalement, il y aurait un moyen de calculer la somme (ou le maximum ou le minimum) de toutes les valeurs absolues (ou les résultats d'un autre "ufunc") sans déplacer toute la sortie vers la RAM et ensuite la récupérer pour la somme/max/min. Il y a eu quelques discussions dans le repo de NumPy, plus récemment à ajouter une ufunc max_abs mais il n'a pas encore été mis en œuvre.

El ufunc.reduce est disponible pour les fonctions à deux entrées comme add o logaddexp mais il n'y a pas de addabs fonction ( x, y : x+abs(y) ) à réduire avec.