Efektywna ocena funkcji w każdej komórce tablicy NumPy

Biorąc pod uwagę tablicę NumPy A , jaki jest najszybszy / najbardziej efektywny sposób zastosowania tej samej funkcji f do każdej komórki?

1. Załóżmy, że będziemy przypisać do A (i, j) do f (A (i, j)) .

2. Funkcja f nie ma wyjścia binarnego, więc operacje maskowania nie pomogą.

Czy „oczywista” iteracja podwójnej pętli (przez każdą komórkę) jest optymalnym rozwiązaniem?

Możesz po prostu wektoryzować funkcję, a następnie zastosować ją bezpośrednio do tablicy Numpy za każdym razem, gdy jej potrzebujesz:

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

Prawdopodobnie lepiej jest określić jawny typ wyjścia bezpośrednio podczas wektoryzacji:

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

Podobne pytanie brzmi: mapowanie tablicy NumPy w miejscu . Jeśli możesz znaleźć ufunc dla swojej f (), powinieneś użyć parametru out.

Jeśli pracujesz z liczbami i f(A(i,j)) = f(A(j,i)), możesz użyć scipy.spatial.distance.cdist, definiując f jako odległość między A(i)a A(j).

Uważam, że znalazłem lepsze rozwiązanie. Pomysł, aby zmienić funkcję na uniwersalną funkcję Pythona (patrz dokumentacja ), która może wykonywać równoległe obliczenia pod maską.

Można napisać swój własny dostosowany ufuncw C, co z pewnością jest bardziej wydajne, lub wywołać np.frompyfunc, co jest wbudowaną metodą fabryczną. Po przetestowaniu jest to bardziej wydajne niż 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

Przetestowałem również większe próbki, a poprawa jest proporcjonalna. Aby porównać wyniki innych metod, zobacz ten post

Gdy tablica 2d (lub tablica nd) jest ciągła C lub F, to zadanie mapowania funkcji na tablicę 2d jest praktycznie takie samo jak zadanie mapowania funkcji na tablicę 1d - po prostu trzeba to zobaczyć w ten sposób, np np.ravel(A,'K'). przez .

Możliwe rozwiązanie dla macierzy 1d zostało omówione na przykład tutaj .

Jednak, gdy pamięć tablicy 2d nie jest ciągła, sytuacja jest nieco bardziej skomplikowana, ponieważ chciałoby się uniknąć ewentualnych błędów w pamięci podręcznej, jeśli oś jest obsługiwana w złej kolejności.

Numpy ma już maszynę do obróbki osi w możliwie najlepszej kolejności. Jedną z możliwości wykorzystania tej maszyny jest np.vectorize. Jednak dokumentacja numpy na temat np.vectorizestwierdza, że ​​jest ona "dostarczana głównie dla wygody, a nie dla wydajności" - powolna funkcja Pythona pozostaje wolną funkcją Pythona z całym związanym z nią narzutem! Inną kwestią jest ogromne zużycie pamięci - zobacz na przykład ten post SO .

Kiedy ktoś chce mieć działanie funkcji C, ale użyć maszyny numpy, dobrym rozwiązaniem jest użycie numba do tworzenia ufuncs, na przykład:

# 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

Z łatwością bije, np.vectorizeale także wtedy, gdy ta sama funkcja byłaby wykonywana jako mnożenie / dodawanie tablicy numpy, tj

# 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"

Zobacz dodatek do tej odpowiedzi dla kodu pomiaru czasu:

Wersja Numby (zielona) jest około 100 razy szybsza niż funkcja Pythona (tj. np.vectorize), Co nie jest zaskakujące. Ale jest również około 10 razy szybsza niż funkcja numpy, ponieważ wersja numbas nie wymaga tablic pośrednich, a zatem bardziej wydajnie wykorzystuje pamięć podręczną.

Chociaż podejście ufunc firmy numba jest dobrym kompromisem między użytecznością a wydajnością, nadal nie jest to najlepsze, co możemy zrobić. Nie ma jednak srebrnej kuli ani podejścia, które najlepiej nadaje się do każdego zadania - trzeba zrozumieć, jakie są ograniczenia i jak można je złagodzić.

Na przykład, dla transcendentalnych funkcji (np exp, sin, cos) Numba nie daje żadnych korzyści w porównaniu z numpy użytkownika np.exp(nie ma tymczasowe tablice tworzone - głównym źródłem prędkości-up). Jednak moja instalacja Anacondy wykorzystuje VML Intela dla wektorów większych niż 8192 - po prostu nie może tego zrobić, jeśli pamięć nie jest ciągła. Dlatego lepiej byłoby skopiować elementy do ciągłej pamięci, aby móc używać VML Intela:

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) 

Dla uczciwości porównania wyłączyłem równoległość VML (zobacz kod w załączniku):

Jak widać, po uruchomieniu VML, narzut kopiowania jest więcej niż kompensowany. Jednak gdy dane stają się zbyt duże dla pamięci podręcznej L3, korzyść jest minimalna, ponieważ zadanie ponownie wiąże się z przepustowością pamięci.

Z drugiej strony numba może również używać SVML Intela, jak wyjaśniono w tym poście :

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)

i używając języka VML z wynikami zrównoleglania:

Wersja numba ma mniej narzutów, ale dla niektórych rozmiarów VML bije SVML nawet pomimo dodatkowego narzutu kopiowania - co nie jest zaskoczeniem, ponieważ ufunks numba nie jest zrównoleglony.

Aukcje:

A. porównanie funkcji wielomianu:

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. porównanie 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)',
    )

Wszystkie powyższe odpowiedzi dobrze się porównują, ale jeśli potrzebujesz użyć niestandardowej funkcji do mapowania, a masz numpy.ndarray, i musisz zachować kształt tablicy.

Porównałem tylko dwa, ale zachowa kształt ndarray. Do porównania użyłem tablicy z 1 milionem wpisów. Tutaj używam funkcji kwadratowej. Przedstawię ogólny przypadek dla tablicy n-wymiarowej. W przypadku dwuwymiarowych po prostu wykonaj iter2D.

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)

Wynik

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

tutaj możesz wyraźnie zobaczyć numpy.fromiterfunkcję kwadratu użytkownika, użyj dowolnego wyboru. Jeśli funkcja jest zależna od i, j indeksów tablicy, iteruj po rozmiarze tablicy, takim jak for ind in range(arr.size), użyj, numpy.unravel_indexaby uzyskać i, j, ..na podstawie indeksu 1D i kształtu tablicy numpy.unravel_index

Ta odpowiedź jest inspirowana moją odpowiedzią na inne pytanie tutaj