Évaluation efficace d'une fonction à chaque cellule d'un tableau NumPy

Étant donné un tableau NumPy A , quel est le moyen le plus rapide / le plus efficace d'appliquer la même fonction, f , à chaque cellule?

1. Supposons que l'on attribue à A (i, j) le f (A (i, j)) .

2. La fonction, f , n'a pas de sortie binaire, donc les opérations de masquage n'aideront pas.

L'itération double boucle "évidente" (à travers chaque cellule) est-elle la solution optimale?

Vous pouvez simplement vectoriser la fonction, puis l'appliquer directement à un tableau Numpy à chaque fois que vous en avez besoin:

import numpy as np

def f(x):
    return x * x + 3 * x - 2 if x > 0 else x * 5 + 8

f = np.vectorize(f)  # or use a different name if you want to keep the original f

result_array = f(A)  # if A is your Numpy array

Il est probablement préférable de spécifier un type de sortie explicite directement lors de la vectorisation:

f = np.vectorize(f, otypes=[np.float])

Une question similaire est: mapper un tableau NumPy en place . Si vous pouvez trouver un ufunc pour votre f (), vous devez utiliser le paramètre out.

Si vous travaillez avec des nombres et f(A(i,j)) = f(A(j,i)), vous pouvez utiliser scipy.spatial.distance.cdist définissant f comme une distance entre A(i)et A(j).

Je crois avoir trouvé une meilleure solution. L'idée de changer la fonction en fonction universelle python (voir documentation ), qui permet d'exercer un calcul parallèle sous le capot.

On peut écrire son propre personnalisé ufuncen C, ce qui est sûrement plus efficace, ou en invoquant np.frompyfunc, qui est une méthode d'usine intégrée. Après les tests, c'est plus efficace que np.vectorize:

f = lambda x, y: x * y
f_arr = np.frompyfunc(f, 2, 1)
vf = np.vectorize(f)
arr = np.linspace(0, 1, 10000)

%timeit f_arr(arr, arr) # 307ms
%timeit f_arr(arr, arr) # 450ms

J'ai également testé des échantillons plus grands et l'amélioration est proportionnelle. Pour une comparaison des performances d'autres méthodes, voir cet article

Lorsque le tableau 2d (ou nd-tableau) est contigu en C ou F, alors cette tâche de mappage d'une fonction sur un tableau 2d est pratiquement la même que la tâche de mappage d'une fonction sur un tableau 1d - nous venons de doivent le voir de cette façon, par exemple via np.ravel(A,'K').

Une solution possible pour 1d-array a été discutée par exemple ici .

Cependant, lorsque la mémoire du 2d-array n'est pas contiguë, alors la situation est un peu plus compliquée, car on voudrait éviter d'éventuels échecs de cache si les axes sont traités dans le mauvais ordre.

Numpy dispose déjà d'une machine pour traiter les axes dans le meilleur ordre possible. Une possibilité d'utiliser cette machine est np.vectorize. Cependant, la documentation de numpy np.vectorizeindique qu'elle est "fournie principalement pour la commodité, pas pour les performances" - une fonction python lente reste une fonction python lente avec toute la surcharge associée! Un autre problème est son énorme consommation de mémoire - voir par exemple ce message SO .

Quand on veut avoir une performance d'une fonction C mais utiliser la machinerie de numpy, une bonne solution est d'utiliser numba pour la création d'ufuncs, par exemple:

# runtime generated C-function as ufunc
import numba as nb
@nb.vectorize(target="cpu")
def nb_vf(x):
    return x+2*x*x+4*x*x*x

Il bat facilement np.vectorizemais aussi quand la même fonction serait exécutée que la multiplication / addition de numpy-array, ie

# numpy-functionality
def f(x):
    return x+2*x*x+4*x*x*x

# python-function as ufunc
import numpy as np
vf=np.vectorize(f)
vf.__name__="vf"

Voir l'annexe de cette réponse pour le code de mesure du temps:

La version de Numba (verte) est environ 100 fois plus rapide que la fonction python (ie np.vectorize), ce qui n'est pas surprenant. Mais c'est également environ 10 fois plus rapide que la fonctionnalité numpy, car la version numbas n'a pas besoin de tableaux intermédiaires et utilise donc le cache plus efficacement.

Bien que l'approche ufunc de numba soit un bon compromis entre convivialité et performances, ce n'est toujours pas le mieux que nous puissions faire. Pourtant, il n’existe pas de solution miracle ou d’approche idéale pour une tâche quelconque - il faut comprendre quelles sont les limites et comment elles peuvent être atténuées.

Par exemple, pour les fonctions transcendantes (par exemple exp, sin, cos) numba ne fournit pas d'avantages par rapport de numpy np.exp(il n'y a pas de tableaux temporaires créés - la principale source de la vitesse-up). Cependant, mon installation Anaconda utilise le VML d'Intel pour les vecteurs supérieurs à 8192 - il ne peut tout simplement pas le faire si la mémoire n'est pas contiguë. Il serait donc préférable de copier les éléments dans une mémoire contiguë afin de pouvoir utiliser le VML d'Intel:

import numba as nb
@nb.vectorize(target="cpu")
def nb_vexp(x):
    return np.exp(x)

def np_copy_exp(x):
    copy = np.ravel(x, 'K')
    return np.exp(copy).reshape(x.shape) 

Pour l'équité de la comparaison, j'ai désactivé la parallélisation de VML (voir code en annexe):

Comme on peut le voir, une fois que VML entre en jeu, la surcharge de copie est plus que compensée. Pourtant, une fois que les données deviennent trop volumineuses pour le cache L3, l'avantage est minime car la tâche redevient liée à la bande passante mémoire.

D'un autre côté, numba pourrait également utiliser le SVML d'Intel, comme expliqué dans cet article :

from llvmlite import binding
# set before import
binding.set_option('SVML', '-vector-library=SVML')

import numba as nb

@nb.vectorize(target="cpu")
def nb_vexp_svml(x):
    return np.exp(x)

et l'utilisation de VML avec des rendements de parallélisation:

La version de numba a moins de frais généraux, mais pour certaines tailles, VML bat SVML même malgré la surcharge de copie supplémentaire - ce qui n'est pas un peu surprenant car les ufuncs de numba ne sont pas parallélisés.

Annonces:

A. comparaison de la fonction polynomiale:

import perfplot
perfplot.show(
    setup=lambda n: np.random.rand(n,n)[::2,::2],
    n_range=[2**k for k in range(0,12)],
    kernels=[
        f,
        vf, 
        nb_vf
        ],
    logx=True,
    logy=True,
    xlabel='len(x)'
    ) 

B. comparaison de exp:

import perfplot
import numexpr as ne # using ne is the easiest way to set vml_num_threads
ne.set_vml_num_threads(1)
perfplot.show(
    setup=lambda n: np.random.rand(n,n)[::2,::2],
    n_range=[2**k for k in range(0,12)],
    kernels=[
        nb_vexp, 
        np.exp,
        np_copy_exp,
        ],
    logx=True,
    logy=True,
    xlabel='len(x)',
    )

Toutes les réponses ci-dessus se comparent bien, mais si vous devez utiliser une fonction personnalisée pour le mappage, et que vous l'avez fait numpy.ndarray, vous devez conserver la forme du tableau.

Je n'ai comparer que deux, mais il conservera la forme de ndarray. J'ai utilisé le tableau avec 1 million d'entrées à des fins de comparaison. Ici, j'utilise la fonction carrée. Je présente le cas général du tableau à n dimensions. Pour deux dimensions, faites simplement iterpour 2D.

import numpy, time

def A(e):
    return e * e

def timeit():
    y = numpy.arange(1000000)
    now = time.time()
    numpy.array([A(x) for x in y.reshape(-1)]).reshape(y.shape)        
    print(time.time() - now)
    now = time.time()
    numpy.fromiter((A(x) for x in y.reshape(-1)), y.dtype).reshape(y.shape)
    print(time.time() - now)
    now = time.time()
    numpy.square(y)  
    print(time.time() - now)

Production

>>> timeit()
1.162431240081787    # list comprehension and then building numpy array
1.0775556564331055   # from numpy.fromiter
0.002948284149169922 # using inbuilt function

ici, vous pouvez voir clairement la numpy.fromiterfonction carrée de l'utilisateur, utilisez celle de votre choix. Si votre fonction dépend des i, j indices du tableau, itérez sur la taille du tableau comme for ind in range(arr.size), utilisez numpy.unravel_indexpour obtenir en i, j, ..fonction de votre index 1D et de la forme du tableau numpy.unravel_index

Cette réponse est inspirée de ma réponse à une autre question ici