Effiziente Bewertung einer Funktion in jeder Zelle eines NumPy-Arrays

Was ist bei einem NumPy- Array A der schnellste / effizienteste Weg, um dieselbe Funktion f auf jede Zelle anzuwenden ?

1. Angenommen, wir weisen A (i, j) das f (A (i, j)) zu .

2. Die Funktion f hat keinen Binärausgang, daher helfen die Maskierungsoperationen nicht.

Ist die "offensichtliche" Doppelschleifeniteration (durch jede Zelle) die optimale Lösung?

Sie können die Funktion einfach vektorisieren und dann jedes Mal, wenn Sie sie benötigen, direkt auf ein Numpy-Array anwenden:

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

Es ist wahrscheinlich besser, einen expliziten Ausgabetyp direkt beim Vektorisieren anzugeben:

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

Eine ähnliche Frage lautet: Zuordnung eines NumPy-Arrays an Ort und Stelle . Wenn Sie ein Ufunc für Ihr f () finden, sollten Sie den Parameter out verwenden.

Wenn Sie mit Zahlen und arbeiten f(A(i,j)) = f(A(j,i)), können Sie scipy.spatial.distance.cdist verwenden , um f als Abstand zwischen A(i)und zu definieren A(j).

Ich glaube, ich habe eine bessere Lösung gefunden. Die Idee, die Funktion in Python Universal Function zu ändern (siehe Dokumentation ), die parallele Berechnungen unter der Haube durchführen kann.

Man kann sein eigenes Customized ufuncin C schreiben , was sicherlich effizienter ist, oder durch Aufrufen np.frompyfuncder eingebauten Factory-Methode. Nach dem Testen ist dies effizienter als 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

Ich habe auch größere Proben getestet und die Verbesserung ist proportional. Einen Vergleich der Leistungen anderer Methoden finden Sie in diesem Beitrag

Wenn das 2d-Array (oder nd-Array) C- oder F-zusammenhängend ist, ist diese Aufgabe, eine Funktion auf ein 2d-Array abzubilden, praktisch dieselbe wie die Aufgabe, eine Funktion auf ein 1d-Array abzubilden - wir nur muss es so sehen, zB via np.ravel(A,'K').

Mögliche Lösung für 1d-Array wurde beispielsweise diskutiert hier .

Wenn jedoch der Speicher des 2d-Arrays nicht zusammenhängend ist, ist die Situation etwas komplizierter, da mögliche Cache-Fehler vermieden werden sollen, wenn die Achse in falscher Reihenfolge behandelt wird.

Numpy verfügt bereits über eine Maschinerie, um Achsen in der bestmöglichen Reihenfolge zu verarbeiten. Eine Möglichkeit, diese Maschine zu benutzen, ist np.vectorize. In der Dokumentation von numpy np.vectorizeheißt es jedoch, dass es "hauptsächlich aus Gründen der Benutzerfreundlichkeit und nicht der Leistung bereitgestellt wird" - eine langsame Python-Funktion bleibt eine langsame Python-Funktion mit dem gesamten damit verbundenen Overhead! Ein weiteres Problem ist der enorme Speicherverbrauch - siehe zum Beispiel diesen SO-Beitrag .

Wenn man eine C-Funktion ausführen möchte, aber die Maschinen von numpy verwenden möchte, ist es eine gute Lösung, numba für die Erstellung von Ufuncs zu verwenden, zum Beispiel:

# 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

Es ist leicht np.vectorizezu schlagen, aber auch, wenn dieselbe Funktion als Multiplikation / Addition von Numpy-Arrays ausgeführt wird, d. H.

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

Im Anhang dieser Antwort finden Sie den Zeitmesscode:

Numbas Version (grün) ist ungefähr 100-mal schneller als die Python-Funktion (dh np.vectorize), was nicht überraschend ist. Es ist aber auch etwa zehnmal schneller als die Numpy-Funktionalität, da die Numbas-Version keine Zwischen-Arrays benötigt und somit den Cache effizienter nutzt.

Der ufunc-Ansatz von numba ist zwar ein guter Kompromiss zwischen Benutzerfreundlichkeit und Leistung, aber immer noch nicht das Beste, was wir tun können. Es gibt jedoch keine Silberkugel oder einen Ansatz, der für eine Aufgabe am besten geeignet ist - man muss verstehen, wo die Grenzen liegen und wie sie gemindert werden können.

Zum Beispiel für transzendentale Funktionen (zB exp, sin, cos) numba bietet keine Vorteile gegenüber der numpy np.exp(es gibt keine temporären Arrays erstellt - die Hauptquelle der Speed-up). Meine Anaconda-Installation verwendet jedoch Intels VML für Vektoren, die größer als 8192 sind - dies ist einfach nicht möglich, wenn der Speicher nicht zusammenhängend ist. Daher ist es möglicherweise besser, die Elemente in einen zusammenhängenden Speicher zu kopieren, um Intels VML verwenden zu können:

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) 

Aus Gründen der Fairness des Vergleichs habe ich die Parallelisierung von VML deaktiviert (siehe Code im Anhang):

Wie man sehen kann, wird der Kopieraufwand nach dem Start von VML mehr als ausgeglichen. Sobald jedoch die Daten für den L3-Cache zu groß werden, ist der Vorteil minimal, da die Aufgabe erneut an die Speicherbandbreite gebunden wird.

Auf der anderen Seite könnte numba auch Intels SVML verwenden, wie in diesem Beitrag erläutert :

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)

und Verwendung von VML mit Parallelisierungsausbeuten:

Die Version von numba hat weniger Overhead, aber für einige Größen schlägt VML SVML trotz des zusätzlichen Kopieraufwands - was keine Überraschung ist, da die ufuncs von numba nicht parallelisiert sind.

Auflistungen:

A. Vergleich der Polynomfunktion:

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. Vergleich von 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)',
    )

Alle obigen Antworten lassen sich gut vergleichen, aber wenn Sie eine benutzerdefinierte Funktion für die Zuordnung verwenden müssen und dies auch tun, müssen numpy.ndarraySie die Form des Arrays beibehalten.

Ich habe nur zwei verglichen, aber es wird die Form von behalten ndarray. Ich habe das Array mit 1 Million Einträgen zum Vergleich verwendet. Hier benutze ich die Quadratfunktion. Ich präsentiere den allgemeinen Fall für ein n-dimensionales Array. Für zweidimensional machen Sie einfach iterfür 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)

Ausgabe

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

Hier können Sie die numpy.fromiterBenutzerquadratfunktion deutlich sehen. Verwenden Sie eine beliebige Funktion Ihrer Wahl. Wenn Ihre Funktion davon abhängt, i, j dass es sich um Array-Indizes handelt, iterieren Sie über die Größe des Arrays for ind in range(arr.size). Verwenden Sie diese Optionnumpy.unravel_index , um i, j, ..anhand Ihres 1D-Index und der Form des Arrays numpy.unravel_index zu ermitteln

Diese Antwort ist inspiriert von meiner Antwort auf eine andere Frage hier