Ich mag NumPy nicht

• Der Autor erläutert seine Kritik an NumPy anhand mehrerer Beispiele.
• Einfache Array-Operationen sind mit NumPy leicht, aber mit steigender Dimensionalität nehmen Komplexität und Verwirrung rapide zu.
• Broadcasting und Advanced Indexing zeigen, dass NumPys Design bei Klarheit und Abstraktion Schwächen hat.
• Statt Achsen explizit anzugeben, muss man beim Schreiben von Code auf Vermutungen und Trial-and-Error setzen.
• Er skizziert Ideen für eine verbesserte Arraysprache und will im nächsten Beitrag konkrete Alternativen vorstellen.

---

Einleitung: Hassliebe zu NumPy

• Der Autor erklärt, dass er NumPy lange verwendet hat, aber von seinen Grenzen vielfach enttäuscht wurde.
• NumPy ist eine unverzichtbare und einflussreiche Bibliothek für Array-Operationen in Python.
• Auch moderne Machine-Learning-Bibliotheken wie PyTorch haben ähnliche Probleme wie NumPy.

Was an NumPy leicht ist – und was nicht

• Einfache Operationen wie das Lösen grundlegender linearer Gleichungen sind mit klarer und eleganter Syntax möglich.
• Sobald Arrays jedoch höherdimensional werden oder Operationen komplexer sind, braucht man Batch-Verarbeitung ohne for-Schleifen.
• In Umgebungen, in denen Schleifen nicht möglich sind, etwa bei GPU-Berechnungen, sind ungewöhnliche Vektorisierungssyntax oder spezielle Funktionsaufrufe nötig.
• Die korrekte Verwendung dieser Funktionen ist jedoch oft unklar und selbst mit der Dokumentation schwer eindeutig zu verstehen.
• Bei der Funktion linalg.solve von NumPy ist es in der Praxis schwer, bei höherdimensionalen Arrays sicher zu wissen, wie sie korrekt eingesetzt werden muss.

Die Probleme von NumPy

• NumPy fehlt eine konsistente Theorie dafür, wie Operationen auf Teile oder bestimmte Achsen mehrdimensionaler Arrays angewendet werden.
• Bei Arrays mit bis zu zwei Dimensionen ist vieles klar, ab drei oder mehr Dimensionen ist jedoch bei jedem Array unklar, auf welche Achsen sich eine Operation bezieht.
• Um Dimensionen explizit anzugleichen, wird man zu komplizierten Verfahren wie None, Broadcasting oder np.tensordot gezwungen.
• Das führt leicht zu Fehlern, verschlechtert die Lesbarkeit des Codes und erhöht die Wahrscheinlichkeit von Bugs.

Schleifen und Klarheit

• Wenn Schleifen tatsächlich erlaubt sind, lässt sich oft kürzerer und klarerer Code schreiben.
• Schleifencode mag weniger elegant wirken, hat aber beim Aspekt der Verständlichkeit große Vorteile.
• Wenn sich dagegen die Array-Dimensionen ändern, muss man sich einzeln mit transpose oder der Reihenfolge der Achsen befassen, was die Komplexität erhöht.

np.einsum: eine ausnahmsweise gute Funktion

• np.einsum ist leistungsstark, weil es eine flexible domänenspezifische Sprache bietet, in der sich Achsen benennen lassen.
• Mit einsum ist die Absicht einer Operation klar, und die Generalisierung funktioniert hervorragend, sodass komplexe Achsenoperationen explizit implementiert werden können.
• Allerdings ist diese Art von Operationsunterstützung nur auf einige Operationen beschränkt; für linalg.solve lässt sie sich zum Beispiel nicht verwenden.

Die Probleme mit Broadcasting

• Broadcasting, der zentrale Trick von NumPy, passt Dimensionen automatisch an, wenn sie nicht übereinstimmen.
• In einfachen Fällen ist das praktisch, in der Praxis macht es aber oft schwer, die Dimensionen klar zu erkennen, und führt zu vielen Fehlersituationen.
• Weil Broadcasting implizit ist, muss man beim Lesen des Codes jedes Mal erneut prüfen, wie die Operation tatsächlich funktioniert.

Die Unklarheit beim Indexing

• NumPys Advanced Indexing macht die Vorhersage der Array-Shape sehr schwierig und unklar.
• Je nach Kombination der Indexing-Varianten ändert sich die Shape des Ergebnis-Arrays, sodass Vorhersagen ohne praktische Erfahrung schwerfallen.
• Auch die Dokumentation der Indexing-Regeln ist lang und komplex und kostet viel Zeit beim Erlernen.
• Selbst wenn man nur einfaches Indexing verwenden möchte, zwingen manche Operationen einen letztlich doch zu Advanced Indexing.

Grenzen im Design der NumPy-Funktionen

• Viele NumPy-Funktionen sind nur für bestimmte Array-Shapes optimiert.
• Für höherdimensionale Arrays braucht man zusätzliche axes-Argumente, separate Funktionsnamen oder Konventionen, und das ist von Funktion zu Funktion inkonsistent.
• Diese Struktur läuft den Grundprinzipien der Programmierung entgegen, bei denen Abstraktion und Wiederverwendung zentral sind.
• Selbst wenn man eine Funktion für ein bestimmtes Problem gefunden hat, muss man sie für verschiedene Arrays und Achsen oft mit völlig anderem Code neu schreiben.

Praxisbeispiel: Implementierung von Self-Attention

• Wenn man Self-Attention in NumPy implementiert, ist Code mit Schleifen klarer, während er bei erzwungener Vektorisierung schnell komplex wird.
• Für höherdimensionale Operationen wie Multi-Head-Attention muss man einsum und Achsentransformationen kombinieren, wodurch der Code schwer verständlich wird.

Fazit und Alternativen

• Der Autor erklärt, NumPy sei „in vieler Hinsicht schlechter als andere Array-Sprachen, aber zugleich die einzige Option, die am Markt entsprechend wichtig geworden ist“.
• Um NumPys verschiedene Probleme wie Broadcasting, Unklarheit beim Indexing und Inkonsistenz der Funktionen zu überwinden, kündigt er an, einen Prototyp einer verbesserten Arraysprache gebaut zu haben.
• Konkrete Verbesserungsvorschläge, also eine neue Arraysprachen-API, will er in einem späteren separaten Beitrag vorstellen.

[youn17]

Das klingt so, als würde es erzählen, warum Julia entstanden ist. Man muss zwar die Bibliotheken erst lernen, aber weil es viele Probleme von NumPy löst, wirkt es wie eine wirklich attraktive Option.

[ahwjdekf]

Wenn man vectorization in NumPy nicht richtig nutzt, ist die Performance dahin. Das unter Berücksichtigung solcher Dinge zu schreiben, ist stressig und schwierig.

[domino]

Ich glaube, viele etwas ältere Python-Bibliotheken haben alle ein ähnliches Problem.

[GN⁺]

Hacker-News-Kommentare

• Im ersten Beispiel ist es schwer lesbar, wenn man in der Dokumentation nur auf den Typ von b schaut. Da aber die zurückgegebene shape beschrieben ist, muss man prüfen, ob der Vektor b tatsächlich als Matrixform behandelt wird, insbesondere wenn K=1 ist.
• Wenn ein Array mehr als zwei Dimensionen hat, empfiehlt sich Xarray, das NumPy-Arrays Dimensionsnamen hinzufügt. Broadcasting und Ausrichtung funktionieren dann automatisch, ohne dass man Dimensionen manuell anpassen oder transpose aufrufen muss, wodurch sich die meisten dieser Probleme erledigen. Xarray ist bei linearer Algebra schwächer als NumPy, aber man kann leicht zu NumPy zurückwechseln und braucht dafür nur ein paar Hilfsfunktionen. Mit Xarray steigt die Produktivität beim Arbeiten mit Daten ab drei Dimensionen deutlich.
  • Xarray fühlt sich an wie eine Kombination aus den Stärken von Pandas und NumPy. Indizierung wie da.sel(x=some_x).isel(t=-1).mean(["y", "z"]) ist einfach, und weil Dimensionsnamen berücksichtigt werden, ist auch Broadcasting klarer. Es ist stark bei der Verarbeitung georäumlicher Daten mit mehreren CRSs. Auch mit Arviz lässt es sich hervorragend nutzen, sodass zusätzliche Dimensionen in der Bayes-Analyse leichter handhabbar sind. Mehrere Arrays lassen sich zu einem Dataset bündeln und gemeinsame Koordinaten teilen, sodass man mit etwas wie ds.isel(t=-1) leicht auf alle Arrays mit Zeitachse gleichzeitig arbeiten kann.
  • Dank Xarray brauche ich elementares NumPy viel seltener und bin deutlich produktiver.
  • Ich frage mich, ob es für Frameworks wie Tensorflow, Keras oder Pytorch etwas Ähnliches gibt. Ich erinnere mich, dass ich etwas in dieser Richtung einmal mühsam debuggen musste.
  • Danke für den Hinweis, ich werde es auf jeden Fall ausprobieren. Ich dachte immer, nur ich fände Syntax wie array[:, :, None] unbequem, daher freut es mich, dieselbe Meinung zu lesen.
  • Im Biosignal-Bereich unterstützt NeuroPype auf Basis von NumPy benannte Achsen für n-dimensionale Tensoren sowie die Speicherung von Per-Element-Daten je Achse, etwa Kanalnamen oder Positionen.
  • Das erinnert an die Zeit, als NumPy aus den alten Bibliotheken Numeric und Numarray hervorging. Ich stelle mir vor, wie das Numarray-Lager zwanzig Jahre lang an seiner Position festhält, dann Finanzierung bekommt, sich in Xarray umbenennt und nun doch NumPy besiegt hat. Natürlich ist das größtenteils erfunden.
• Einer der Gründe, warum ich angefangen habe, Julia zu nutzen, war die viel zu schwierige NumPy-Syntax. Beim Wechsel von MATLAB zu NumPy hatte ich das Gefühl, programmatisch schlechter zu werden, weil ich mehr Zeit mit Performance-Tricks als mit Mathematik verbrachte. In Julia funktionieren sowohl Vektorisierung als auch Schleifen gut, sodass ich mich nur um die Lesbarkeit des Codes kümmern muss. Genau diese Erfahrung und dieses Gefühl kamen in dem Text für mich rüber. Ich halte auch den „Blackbox“-Ansatz für falsch, bei dem Dinge wie np.linalg.solve als maximal schnell gelten und man sie deshalb immer verwenden soll. Es gibt viele Gründe, problemspezifische Kernel selbst zu schreiben.
  • Der Grund ist, dass Julia eine für wissenschaftliches Rechnen entworfene Sprache ist, während NumPy nur eine Bibliothek ist, die auf eine Sprache gesetzt wurde, die nicht dafür entworfen wurde. Hoffentlich setzt sich Julia irgendwann durch und befreit die Leute, die wegen Netzwerkeffekten Python benutzen.
  • MATLAB ist ohne Vektorisierung und mit Schleifen auch ungefähr so langsam wie Python. Die eigentliche Hauptschwäche ist die Langsamkeit von Python. Julia hat klar Vorteile, lässt sich in der Praxis aber nur für stark eingeschränkte Einsatzfälle nutzen. In Python gibt es inzwischen JIT-Hacks, aber sie sind noch immer unvollständig. Eine Alternative zu Python wird dringend gebraucht.
  • Ist MATLAB wirklich so anders? Schleifen bleiben langsam, und am schnellsten ist weiterhin die vollständig optimierte Blackbox wie der \-Operator.
  • Auch moderne Fortran-Versionen sind wie Julia darin, dass sowohl Vektorisierung als auch Schleifen schnell funktionieren, sodass man sich nur auf die Lesbarkeit konzentrieren kann.
• Wenn man die Beschwerden über NumPy im Vergleich zu Matlab und Julia zusammenfasst, dann sind die Achsen-bezogenen Argumente, Benennungen und Arten der Vektorisierung von Funktion zu Funktion inkonsistent. Wenn man eine Funktion auf eine andere Achse anwenden will, muss man den Code oft komplett umschreiben. Abstraktion ist ein Grundprinzip des Programmierens, und NumPy macht genau das schwer. In Matlab läuft vektorisierter Code meist fast unverändert weiter oder die nötigen Anpassungen sind klar, während man in NumPy ständig in die Dokumentation schauen muss und das Angleichen von Typen mit transpose, reshape usw. nicht konsistent ist und dadurch unerquicklich wird.
  • Die Unterstützung von Matlab für Arrays mit mehr als drei Dimensionen ist so schwach, dass die im Artikel genannten Probleme dort eher gar nicht erst auftreten.
  • Für das zweite Problem könnte man jax mit vmap ausprobieren.
  • Wenn man eine bestimmte Funktion für ein 2x2-Array schreibt und sie dann auf Teile eines 3x2x2-Arrays anwenden will, geht das mit Slices und squeeze usw. Dieses Problem ist so vage, dass ich es kaum nachvollziehen kann.
  • Das lässt sich mit reshape lösen.
• Das Verwirrendste an NumPy ist für mich, dass oft unklar ist, welche Operationen vektorisiert funktionieren, und dass man das nicht wie in Julia mit Punkt-Syntax explizit machen kann. Auch rund um Rückgabetypen gibt es viele Fallstricke. Wenn man zum Beispiel ein poly1d-Objekt P von rechts mit z0 multipliziert, erhält man wieder ein poly1d; multipliziert man aber von links in der Form z0*P, wird nur ein Array zurückgegeben, sodass stillschweigend eine Typumwandlung stattfindet. Selbst beim führenden Koeffizienten eines quadratischen Polynoms gibt es zwei Schreibweisen, P.coef[0] und P[2], was leicht verwirrend wird. Offiziell ist poly1d eine „alte“ API und für neuen Code wird die Klasse Polynomial empfohlen, aber tatsächlich gibt es nicht einmal eine Deprecated-Warnung. Solche Typumwandlungen und Inkonsistenzen bei Datentypen sind überall in der Bibliothek versteckte Minenfelder und machen das Debugging zum Albtraum.
• Ich kann die Punkte des Autors nachvollziehen. Beim Wechsel von Matlab zu Numpy gab es für mich viele Unannehmlichkeiten, und auch Data Slicing fühlt sich in Numpy unbequemer an als in Matlab oder Julia. Wenn man allerdings die Lizenzkosten der Matlab-Toolboxes bedenkt, werden die Nachteile von Numpy dadurch ausgeglichen. Die im Text genannten Probleme treten hauptsächlich bei Tensoren mit mehr als zwei Dimensionen auf, und Numpy ist ursprünglich matrixorientiert, also 2D-zentriert, daher ist diese Grenze nicht überraschend. Eine spezialisierte Bibliothek wie Torch ist besser, aber auch nicht gerade einfach. Am Ende wirkt die Aussage plausibel: „NumPy ist ein bisschen schlechter als andere Array-Sprachen, aber wirklich viele Alternativen gibt es auch nicht.“
  • NumPy zielte von Anfang an auf n-dimensionale Arrays ab und stand in der Tradition von numarray, es war also nicht nur auf 2D beschränkt.
• Das größte Problem des Python-Data-Science-Ökosystems ist, dass alles nicht standardisiert ist. Ein gutes Dutzend Bibliotheken verhält sich so unterschiedlich wie vier verschiedene Sprachen, und selbst to_numpy() ist nur so halbwegs vereinheitlicht. Am Ende verbringt man mehr Zeit mit der Umwandlung von Datenformaten als mit dem eigentlichen Lösen des Problems. Auch Julia hat nicht nur Vorteile, aber die Interoperabilität zwischen Bibliotheken für Einheiten, Unsicherheiten usw. funktioniert dort gut, während Python immer viel Boilerplate-Code verlangt.
  • Das Projekt array-api versucht, die API für Array-Manipulationen im gesamten Python-Ökosystem zu standardisieren.
  • R ist wegen seiner vier Klassensysteme eher noch komplizierter.
• Ich frage mich, warum Leute NumPy statt Sage verwenden.
• Manche Probleme lassen sich mit numpysane und gnuplotlib lösen. Seit es diese Kombination gibt, nutze ich NumPy aktiv für alle möglichen Aufgaben; ohne sie wäre es für mich kaum brauchbar.
  • numpysane ist letztlich nur eine Python-Schleife, also keine echte Vektorisierung.
  • Danke für den Hinweis. Ich habe mich wegen solcher Probleme schon oft beschwert und nie daran gedacht, dass es dafür eine einfache übergeordnete Bibliothek geben könnte.
• Für vectorized multi-head attention habe ich alle Matrixmultiplikationen in einsum gepackt und mit optimize="optimal" den Matrix-Chain-Multiplikationsalgorithmus genutzt, um die Performance zu verbessern. Tatsächlich war das gegenüber einer gewöhnlichen vektorisierten Implementierung etwa doppelt so schnell, aber überraschenderweise war eine naive Implementierung mit Schleifen noch schneller. Wer wissen will, warum, sollte sich den Code ansehen. Ich vermute, dass es bei der Cache-Kohärenz innerhalb von einsum noch Raum für Verbesserungen gibt.