Warum ist numpys einsum schneller als die integrierten Funktionen von numpy?

72

Beginnen wir mit drei Arrays von dtype=np.double. Timings werden auf einer Intel-CPU unter Verwendung von numpy 1.7.1 durchgeführt, das mit iccIntel kompiliert und mit Intel verknüpft ist mkl. Eine AMD-CPU mit der Nummer 1.6.1, die mit gccohne kompiliert mklwurde, wurde ebenfalls verwendet, um die Timings zu überprüfen. Bitte beachten Sie, dass die Timings nahezu linear mit der Systemgröße skalieren und nicht auf den geringen Overhead zurückzuführen sind, der in den Numpy-Funktionsanweisungen anfällt. Dieser ifUnterschied wird in Mikrosekunden und nicht in Millisekunden angezeigt :

arr_1D=np.arange(500,dtype=np.double)
large_arr_1D=np.arange(100000,dtype=np.double)
arr_2D=np.arange(500**2,dtype=np.double).reshape(500,500)
arr_3D=np.arange(500**3,dtype=np.double).reshape(500,500,500)

Schauen wir uns zunächst die np.sumFunktion an:

np.all(np.sum(arr_3D)==np.einsum('ijk->',arr_3D))
True

%timeit np.sum(arr_3D)
10 loops, best of 3: 142 ms per loop

%timeit np.einsum('ijk->', arr_3D)
10 loops, best of 3: 70.2 ms per loop

Befugnisse:

np.allclose(arr_3D*arr_3D*arr_3D,np.einsum('ijk,ijk,ijk->ijk',arr_3D,arr_3D,arr_3D))
True

%timeit arr_3D*arr_3D*arr_3D
1 loops, best of 3: 1.32 s per loop

%timeit np.einsum('ijk,ijk,ijk->ijk', arr_3D, arr_3D, arr_3D)
1 loops, best of 3: 694 ms per loop

Äußeres Produkt:

np.all(np.outer(arr_1D,arr_1D)==np.einsum('i,k->ik',arr_1D,arr_1D))
True

%timeit np.outer(arr_1D, arr_1D)
1000 loops, best of 3: 411 us per loop

%timeit np.einsum('i,k->ik', arr_1D, arr_1D)
1000 loops, best of 3: 245 us per loop

Alle oben genannten sind doppelt so schnell mit np.einsum. Dies sollten Vergleiche zwischen Äpfeln sein, da alles spezifisch ist dtype=np.double. Ich würde die Geschwindigkeit bei einer Operation wie dieser erwarten:

np.allclose(np.sum(arr_2D*arr_3D),np.einsum('ij,oij->',arr_2D,arr_3D))
True

%timeit np.sum(arr_2D*arr_3D)
1 loops, best of 3: 813 ms per loop

%timeit np.einsum('ij,oij->', arr_2D, arr_3D)
10 loops, best of 3: 85.1 ms per loop

Einsum scheint mindestens doppelt so schnell zu sein np.inner, np.outer, np.kronund np.sumunabhängig davon , axesAuswahl. Die np.dot Hauptausnahme ist der Aufruf von DGEMM aus einer BLAS-Bibliothek. Warum ist es also np.einsumschneller als andere gleichwertige Numpy-Funktionen?

Der DGEMM-Fall der Vollständigkeit halber:

np.allclose(np.dot(arr_2D,arr_2D),np.einsum('ij,jk',arr_2D,arr_2D))
True

%timeit np.einsum('ij,jk',arr_2D,arr_2D)
10 loops, best of 3: 56.1 ms per loop

%timeit np.dot(arr_2D,arr_2D)
100 loops, best of 3: 5.17 ms per loop

Die führende Theorie stammt aus dem Kommentar von @sebergs, np.einsumder SSE2 verwenden kann , aber die Ufuncs von numpy werden erst bei numpy 1.8 angezeigt (siehe Änderungsprotokoll ). Ich glaube, dies ist die richtige Antwort, konnte sie aber nicht bestätigen. Einige begrenzte Beweise können gefunden werden, indem der Typ des Eingabearrays geändert und der Geschwindigkeitsunterschied sowie die Tatsache beobachtet werden, dass nicht jeder die gleichen zeitlichen Trends beobachtet.

Daniel
quelle
Mit welcher BLAS-Bibliothek ist numpy verknüpft? Ist es Multithread?
Ali_m
1
Multithread-MKL BLAS mit AVX.
Daniel
Übrigens tolle Frage und gute Beispiele! Es könnte sich lohnen, dies auf der Mailingliste zu erfragen. Es ist vor ( vor allem in Bezug auf abgedeckt sum), aber ich bin überrascht , dass einsumkonsequent ~ 2x schneller als ist outer, inner, kronusw. Es wäre interessant zu wissen , wo der Unterschied kommen.
Joe Kington
@ JoeKington Ich denke, ich werde es auf die Mailingliste setzen, wenn jemand anderes die ~ 2x Beschleunigung reproduzieren kann. Seltsamerweise zeigt Jamies Antwort dies.
Daniel
etwas verwandt: stackoverflow.com/questions/17527340/… aber in diesem Fall scheint der Grund für Geschwindigkeitsunterschiede die Speicherverwaltung zu sein (wenn Sie anfangen,
Dinge

Antworten:

32

Zunächst einmal gab es in der Vergangenheit viele Diskussionen darüber auf der Numpy-Liste. Siehe zum Beispiel: http://numpy-discussion.10968.n7.nabble.com/poor-performance-of-sum-with-sub-machine-word-integer-types-td41.html http: // numpy- Diskussion.10968.n7.nabble.com/odd-performance-of-sum-td3332.html

Einige davon sind auf die Tatsache zurückzuführen, dass einsumes neu ist und vermutlich versucht, die Cache-Ausrichtung und andere Probleme beim Speicherzugriff besser zu behandeln, während sich viele der älteren Numpy-Funktionen auf eine leicht portierbare Implementierung gegenüber einer stark optimierten konzentrieren. Ich spekuliere dort allerdings nur.


Einiges von dem, was Sie tun, ist jedoch kein Vergleich von Äpfeln zu Äpfeln.

Zusätzlich zu dem, was @Jamie bereits gesagt hat, sumwird ein geeigneterer Akkumulator für Arrays verwendet

Gehen Sie beispielsweise sumvorsichtiger vor, wenn Sie den Typ des Eingangs überprüfen und einen geeigneten Akku verwenden. Betrachten Sie beispielsweise Folgendes:

In [1]: x = 255 * np.ones(100, dtype=np.uint8)

In [2]: x
Out[2]:
array([255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
       255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
       255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
       255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
       255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
       255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
       255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
       255, 255, 255, 255, 255, 255, 255, 255, 255], dtype=uint8)

Beachten Sie, dass das sumrichtig ist:

In [3]: x.sum()
Out[3]: 25500

Während einsumwird das falsche Ergebnis geben:

In [4]: np.einsum('i->', x)
Out[4]: 156

Wenn wir jedoch eine weniger begrenzte verwenden dtype, erhalten wir immer noch das erwartete Ergebnis:

In [5]: y = 255 * np.ones(100)

In [6]: np.einsum('i->', y)
Out[6]: 25500.0
Joe Kington
quelle
Haben Sie einen guten Link, wie sumder Akku ausgewählt wird? Interessanterweise ist die xErweiterung Ihres Arrays auf 1E8Elemente np.einsum('i->',x,dtype=np.uint64)dann nur etwa 10% schneller (15 ms) sum.
Daniel
@Ophion - Die Dokumentation für sumenthält einige Details. Sie können es mit dem dtypekwarg angeben sum. Wenn es nicht angegeben ist und das Array einen Integer-Typ mit weniger Genauigkeit als die "Standard-Plattform-Ganzzahl" hat (normalerweise int64sogar auf 32-Bit-Plattformen, glaube ich), wird standardmäßig die Standard-Ganzzahl verwendet. Siehe: docs.scipy.org/doc/numpy/reference/generated/numpy.sum.html
Joe Kington
Wird auch sumdurch implementiert np.add.reduce, sehen Sie sich hier die Quelle für Reduzierungen ufuncan, wenn Sie an den Details interessiert sind: github.com/numpy/numpy/blob/master/numpy/core/src/umath/…
Joe Kington
1
Wenn ich es richtig verstehe, handelt es sich um Vergleiche von Äpfeln zu Äpfeln, da sich alles speziell darauf beschränkt dtype=np.double?
Daniel
2
Ich glaube schon. Und genau das haben Sie überhaupt getan. Daher ist der Punkt, den ich angesprochen habe, wahrscheinlich doch nicht so relevant!
Joe Kington
22

Nachdem nun numpy 1.8 veröffentlicht wurde, bei dem laut den Dokumenten alle Ufuncs SSE2 verwenden sollten, wollte ich überprüfen, ob Sebergs Kommentar zu SSE2 gültig ist.

Um den Test durchzuführen, wurde eine neue Python 2.7-Installation erstellt - numpy 1.7 und 1.8 wurden mit iccStandardoptionen auf einem AMD-Opteron-Core unter Ubuntu kompiliert .

Dies ist der Testlauf vor und nach dem 1.8-Upgrade:

import numpy as np
import timeit

arr_1D=np.arange(5000,dtype=np.double)
arr_2D=np.arange(500**2,dtype=np.double).reshape(500,500)
arr_3D=np.arange(500**3,dtype=np.double).reshape(500,500,500)

print 'Summation test:'
print timeit.timeit('np.sum(arr_3D)',
                      'import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D',
                      number=5)/5
print timeit.timeit('np.einsum("ijk->", arr_3D)',
                      'import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D',
                      number=5)/5
print '----------------------\n'


print 'Power test:'
print timeit.timeit('arr_3D*arr_3D*arr_3D',
                      'import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D',
                      number=5)/5
print timeit.timeit('np.einsum("ijk,ijk,ijk->ijk", arr_3D, arr_3D, arr_3D)',
                      'import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D',
                      number=5)/5
print '----------------------\n'


print 'Outer test:'
print timeit.timeit('np.outer(arr_1D, arr_1D)',
                      'import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D',
                      number=5)/5
print timeit.timeit('np.einsum("i,k->ik", arr_1D, arr_1D)',
                      'import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D',
                      number=5)/5
print '----------------------\n'


print 'Einsum test:'
print timeit.timeit('np.sum(arr_2D*arr_3D)',
                      'import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D',
                      number=5)/5
print timeit.timeit('np.einsum("ij,oij->", arr_2D, arr_3D)',
                      'import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D',
                      number=5)/5
print '----------------------\n'

Numpy 1.7.1:

Summation test:
0.172988510132
0.0934836149216
----------------------

Power test:
1.93524689674
0.839519000053
----------------------

Outer test:
0.130380821228
0.121401786804
----------------------

Einsum test:
0.979052495956
0.126066613197

Numpy 1.8:

Summation test:
0.116551589966
0.0920487880707
----------------------

Power test:
1.23683619499
0.815982818604
----------------------

Outer test:
0.131808176041
0.127472200394
----------------------

Einsum test:
0.781750011444
0.129271841049

Ich denke, dies ist ziemlich schlüssig, dass SSE eine große Rolle bei den Timing-Unterschieden spielt. Es sollte beachtet werden, dass das Wiederholen dieser Tests die Timings nur um ~ 0,003s sehr stark beeinflusst. Der verbleibende Unterschied sollte in den anderen Antworten auf diese Frage behandelt werden.

Daniel
quelle
4
Fantastisches Follow-up! Dies ist ein weiterer Grund, warum ich einsumhäufiger anfangen muss . Im Übrigen würde ich argumentieren, dass Sie in diesem Fall Ihre eigene Antwort wirklich als richtig markieren sollten.
Joe Kington
19

Ich denke, diese Zeiten erklären, was los ist:

a = np.arange(1000, dtype=np.double)
%timeit np.einsum('i->', a)
100000 loops, best of 3: 3.32 us per loop
%timeit np.sum(a)
100000 loops, best of 3: 6.84 us per loop

a = np.arange(10000, dtype=np.double)
%timeit np.einsum('i->', a)
100000 loops, best of 3: 12.6 us per loop
%timeit np.sum(a)
100000 loops, best of 3: 16.5 us per loop

a = np.arange(100000, dtype=np.double)
%timeit np.einsum('i->', a)
10000 loops, best of 3: 103 us per loop
%timeit np.sum(a)
10000 loops, best of 3: 109 us per loop

So dass Sie im Grunde eine nahezu konstanten 3us Overhead beim Aufruf np.sumüber np.einsum, so dass sie so schnell im Grunde laufen, aber man nimmt ein wenig länger in Gang zu bringen. Warum könnte das sein? Mein Geld ist für Folgendes:

a = np.arange(1000, dtype=object)
%timeit np.einsum('i->', a)
Traceback (most recent call last):
...
TypeError: invalid data type for einsum
%timeit np.sum(a)
10000 loops, best of 3: 20.3 us per loop

Ich bin mir nicht sicher, was genau vor sich geht, aber es scheint, dass np.einsumeinige Überprüfungen übersprungen werden, um typspezifische Funktionen für die Multiplikationen und Additionen zu extrahieren, und direkt mit *und nur +für Standard-C-Typen.


Die mehrdimensionalen Fälle unterscheiden sich nicht:

n = 10; a = np.arange(n**3, dtype=np.double).reshape(n, n, n)
%timeit np.einsum('ijk->', a)
100000 loops, best of 3: 3.79 us per loop
%timeit np.sum(a)
100000 loops, best of 3: 7.33 us per loop

n = 100; a = np.arange(n**3, dtype=np.double).reshape(n, n, n)
%timeit np.einsum('ijk->', a)
1000 loops, best of 3: 1.2 ms per loop
%timeit np.sum(a)
1000 loops, best of 3: 1.23 ms per loop

Also ein größtenteils konstanter Overhead, kein schnelleres Laufen, wenn sie erst einmal dran sind.

Jaime
quelle
1
In der Dokumentation wird außerdem vorgeschlagen , dass einsumauch keine automatische Übertragung durchgeführt wird, und der Benutzer muss die Übertragungsregeln für eine Operation ausdrücken. Es gibt also wahrscheinlich viele Überprüfungen (Typprüfung, Rundfunk usw.), einsumdie übersprungen werden können.
Ely
Seltsamerweise unterscheiden sie sich auf meinem Computer. Bitte sehen Sie sich meine Bearbeitung an.
Daniel
1 oder mehr Dimensionen sind grundsätzlich dasselbe. np.sumAnrufe np.add.reduce, und das wurde überarbeitet 1.7, um mehrere Achsen zu akzeptieren. Die Iteration wird also mit ziemlicher Sicherheit von einem sehr ähnlichen Aufruf wie das C-Äquivalent von np.nditerin beiden Fällen behandelt. Wenn Sie keine Zwischen-Arrays vermeiden, um das Multiplizieren-dann-Addieren von numpy auszuführen, oder wenn Sie eine Multithread-Bibliothek verwenden, sollten Sie abgesehen von der Einrichtung kleine Unterschiede feststellen, was meine Timings zeigen.
Jaime
8
Sie sollten wahrscheinlich eine 2-fache Beschleunigung mit doppelter Genauigkeit (SSE) sehen. Da die Summe naiv ist (möglicherweise nicht auf 1.8+ nicht sicher), während einsum speziell für die Verwendung von SIMD-Anweisungen geschrieben wurde, tun dies die meisten Ufuncs nicht.
Seberg
1
@seberg Du hast es geschafft, beide Prozessoren haben SSE2, also würde man erwarten, dass die Einzelgenauigkeit 4x so schnell ist und es ist. Wenn Sie dies aufschreiben können, werde ich es akzeptieren.
Daniel
1

Ein Update für numpy 1.16.4: Die nativen Funktionen von Numpy sind in fast allen Fällen schneller als einsums. Nur die äußere Variante von einsum und der sum23-Test sind schneller als die Nicht-einsum-Versionen.

Wenn Sie die nativen Funktionen von numpy verwenden können, tun Sie dies.

(Bilder erstellt mit perfplot , einem Projekt von mir.)

Geben Sie hier die Bildbeschreibung ein

Geben Sie hier die Bildbeschreibung ein

Geben Sie hier die Bildbeschreibung ein

Geben Sie hier die Bildbeschreibung ein

Geben Sie hier die Bildbeschreibung ein

Geben Sie hier die Bildbeschreibung ein


Code zur Reproduktion der Diagramme:

import numpy
import perfplot


def setup1(n):
    return numpy.arange(n, dtype=numpy.double)


def setup2(n):
    return numpy.arange(n ** 2, dtype=numpy.double).reshape(n, n)


def setup3(n):
    return numpy.arange(n ** 3, dtype=numpy.double).reshape(n, n, n)

def setup23(n):
    return (
        numpy.arange(n ** 2, dtype=numpy.double).reshape(n, n),
        numpy.arange(n ** 3, dtype=numpy.double).reshape(n, n, n)
    )


def numpy_sum(a):
    return numpy.sum(a)


def einsum_sum(a):
    return numpy.einsum("ijk->", a)


perfplot.save(
    "sum.png",
    setup=setup3,
    kernels=[numpy_sum, einsum_sum],
    n_range=[2 ** k for k in range(10)],
    logx=True,
    logy=True,
    title="sum",
)


def numpy_power(a):
    return a * a * a


def einsum_power(a):
    return numpy.einsum("ijk,ijk,ijk->ijk", a, a, a)


perfplot.save(
    "power.png",
    setup=setup3,
    kernels=[numpy_power, einsum_power],
    n_range=[2 ** k for k in range(9)],
    logx=True,
    logy=True,
)


def numpy_outer(a):
    return numpy.outer(a, a)


def einsum_outer(a):
    return numpy.einsum("i,k->ik", a, a)


perfplot.save(
    "outer.png",
    setup=setup1,
    kernels=[numpy_outer, einsum_outer],
    n_range=[2 ** k for k in range(13)],
    logx=True,
    logy=True,
    title="outer",
)



def dgemm_numpy(a):
    return numpy.dot(a, a)


def dgemm_einsum(a):
    return numpy.einsum("ij,jk", a, a)


def dgemm_einsum_optimize(a):
    return numpy.einsum("ij,jk", a, a, optimize=True)


perfplot.save(
    "dgemm.png",
    setup=setup2,
    kernels=[dgemm_numpy, dgemm_einsum],
    n_range=[2 ** k for k in range(13)],
    logx=True,
    logy=True,
    title="dgemm",
)



def dot_numpy(a):
    return numpy.dot(a, a)


def dot_einsum(a):
    return numpy.einsum("i,i->", a, a)


perfplot.save(
    "dot.png",
    setup=setup1,
    kernels=[dot_numpy, dot_einsum],
    n_range=[2 ** k for k in range(20)],
    logx=True,
    logy=True,
    title="dot",
)

def sum23_numpy(data):
    a, b = data
    return numpy.sum(a * b)


def sum23_einsum(data):
    a, b = data
    return numpy.einsum('ij,oij->', a, b)


perfplot.save(
    "sum23.png",
    setup=setup23,
    kernels=[sum23_numpy, sum23_einsum],
    n_range=[2 ** k for k in range(10)],
    logx=True,
    logy=True,
    title="sum23",
)
Nico Schlömer
quelle
Ein Hinweis auf dem GEMM, wenn Sie numpy.einsum("ij,jk", a, a, optimize=True)die Leistung gleichwertig sind. Es ist etwas seltsam, dass die Latenzzeit geringer ist. Hat sich die Logik dieser Funktionen nach C verschoben? Auch einen Versuch wert np.einsum('i,i->', ...)sowie np.einsum('ij,oij->'einen für mehr Äpfel zu Äpfeln Vergleich.
Daniel
@ Daniel Fügte diese hinzu.
Nico Schlömer