Je déteste NumPy

• L’auteur expose ses griefs envers NumPy à l’aide de plusieurs exemples concrets
• Les opérations simples sur les tableaux sont faciles avec NumPy, mais dès que le nombre de dimensions augmente, la complexité et la confusion explosent
• La conception de NumPy, notamment le broadcasting et l’indexation avancée, manque de clarté et d’abstraction
• Écrire du code impose de s’appuyer sur des suppositions et des tâtonnements au lieu de spécifier explicitement les axes
• Il propose des pistes pour un meilleur langage de tableaux et présentera des alternatives concrètes dans un prochain article

---

Introduction : relation amour-haine avec NumPy

• L’auteur explique qu’il utilise NumPy depuis longtemps, mais qu’il a souvent été déçu par ses limites
• NumPy est une bibliothèque essentielle et très influente pour le calcul sur tableaux en Python
• Des bibliothèques modernes de machine learning comme PyTorch souffrent elles aussi de problèmes similaires à ceux de NumPy

Ce qui est facile et ce qui est difficile avec NumPy

• Les opérations simples, comme la résolution d’équations linéaires de base, peuvent être exprimées avec une syntaxe claire et élégante
• Mais quand les tableaux ont davantage de dimensions ou que les opérations deviennent complexes, il faut traiter le tout en bloc, sans boucles for
• Dans les environnements où les boucles ne sont pas possibles, comme le calcul sur GPU, il faut recourir à une syntaxe de vectorisation étrange ou à des appels de fonctions spécifiques
• Or, la bonne manière de les utiliser reste ambiguë, et la documentation seule ne permet pas toujours de lever le doute
• En pratique, avec des tableaux de grande dimension, il est difficile d’avoir la certitude d’utiliser correctement la fonction linalg.solve de numpy

Les problèmes de NumPy

• NumPy manque d’une théorie cohérente pour appliquer des opérations à une partie d’un tableau multidimensionnel ou à certains axes précis
• Jusqu’à 2 dimensions, tout reste clair, mais à partir de 3 dimensions, il devient flou de savoir quels axes de chaque tableau sont concernés par l’opération
• Il faut alors recourir à des techniques complexes comme l’usage de None, le broadcasting ou np.tensordot pour aligner explicitement les dimensions
• Ces approches favorisent les erreurs, réduisent la lisibilité du code et augmentent le risque de bugs

Boucles et clarté

• En réalité, si les boucles sont autorisées, il devient possible d’écrire un code plus concis et plus clair
• Le code avec des boucles peut sembler moins sophistiqué, mais il présente un avantage majeur en matière de clarté
• À l’inverse, dès que la dimension des tableaux change, il faut réfléchir à la main au transpose ou à l’ordre des axes, ce qui augmente la complexité

np.einsum : une fonction exceptionnellement réussie

• np.einsum est puissant, car il fournit un langage spécialisé flexible permettant de nommer les axes
• einsum rend l’intention de l’opération claire et se généralise très bien, ce qui permet d’implémenter explicitement des opérations complexes sur les axes
• Mais ce type de prise en charge, proche de einsum, reste limité à certaines opérations et ne s’applique pas, par exemple, à linalg.solve

Les problèmes du broadcasting

• Le broadcasting, astuce centrale de NumPy, aligne automatiquement les dimensions quand elles ne correspondent pas
• Dans les cas simples, c’est pratique, mais en réalité cela rend souvent plus difficile la compréhension explicite des dimensions et provoque de nombreux cas d’erreur
• Comme le broadcasting est implicite, il faut vérifier à chaque lecture du code comment l’opération va réellement se comporter

Le manque de clarté de l’indexation

• L’indexation avancée de NumPy rend la prédiction de la shape d’un tableau particulièrement difficile et peu claire
• Comme la shape du tableau résultant varie selon les combinaisons d’indexation, il est difficile de la prévoir sans expérience pratique directe
• Même la documentation expliquant les règles d’indexation est longue et complexe, ce qui demande beaucoup de temps pour être assimilé
• Même lorsqu’on veut se limiter à une indexation simple, certaines opérations obligent malgré tout à utiliser l’indexation avancée

Les limites de la conception des fonctions NumPy

• Beaucoup de fonctions NumPy sont optimisées uniquement pour certaines shapes de tableaux
• Pour les tableaux de grande dimension, il faut ajouter des arguments axes, utiliser des noms de fonctions différents ou suivre des conventions distinctes, sans cohérence d’ensemble entre les fonctions
• Cette structure va à l’encontre des principes de programmation fondés sur l’abstraction et la réutilisation
• Même lorsqu’une fonction résout un problème précis, il faut souvent réécrire totalement le code pour la réappliquer à d’autres tableaux et à d’autres axes

Exemple concret : implémentation du self-attention

• Lorsqu’on écrit une implémentation de self-attention en NumPy, l’usage de boucles rend le code clair, mais imposer la vectorisation le complique fortement
• Pour des opérations de grande dimension comme le multi-head attention, il faut combiner einsum et des transformations d’axes, ce qui rend le code difficile à lire

Conclusion et alternative

• L’auteur affirme que NumPy est « le seul choix devenu incontournable sur le marché, malgré de nombreux défauts par rapport à d’autres langages de tableaux »
• Il annonce avoir créé un prototype de langage de tableaux amélioré pour dépasser plusieurs problèmes de NumPy, comme le broadcasting, le manque de clarté de l’indexation ou l’incohérence des fonctions
• Il présentera plus tard, dans un autre article, des pistes d’amélioration concrètes, sous la forme d’une nouvelle API de langage de tableaux

[youn17]

On dirait une histoire sur les raisons de la création de Julia. Il faut certes se former aux bibliothèques, mais le fait qu’il résolve beaucoup des problèmes de NumPy en fait vraiment une option très séduisante.

[ahwjdekf]

Si on ne maîtrise pas bien la vectorisation avec numpy, les performances s’effondrent. Devoir écrire le code en tenant compte de ça est stressant et difficile.

[domino]

On dirait que pas mal de bibliothèques Python un peu anciennes ont toutes le même genre de problème.

[GN⁺]

Réactions sur Hacker News

• Dans le premier exemple, si on ne regarde que le type de b dans la documentation, c’est difficile à lire, mais comme il y a une explication sur la shape renvoyée, il faut vérifier si le vecteur b est en fait sous forme matricielle, surtout quand K=1
• Dès que les tableaux ont plus de deux dimensions, je recommande d’utiliser Xarray, qui ajoute des noms de dimensions aux tableaux NumPy ; comme le broadcasting et l’alignement sont automatiques sans avoir à ajuster les dimensions ni faire de transpose, cela résout la plupart de ces problèmes ; Xarray est plus faible que NumPy pour l’algèbre linéaire, mais on peut facilement revenir à NumPy, il suffit de créer quelques fonctions utilitaires ; avec Xarray, la productivité augmente fortement dès qu’on manipule des données en 3D ou plus
  • Xarray donne l’impression de combiner les avantages de Pandas et de NumPy ; l’indexation comme da.sel(x=some_x).isel(t=-1).mean(["y", "z"]) est simple, et comme les noms de dimensions sont respectés, le broadcasting est clair ; c’est aussi très solide pour traiter des données géospatiales avec plusieurs CRS ; l’usage avec ArviZ est excellent, ce qui facilite aussi la gestion de dimensions supplémentaires en analyse bayésienne ; on peut regrouper plusieurs tableaux dans un même dataset pour partager des coordonnées communes, ce qui permet d’appliquer facilement quelque chose comme ds.isel(t=-1) à tous les tableaux ayant un axe temporel
  • Grâce à Xarray, mon usage basique de NumPy a beaucoup diminué et ma productivité a nettement augmenté
  • Je me demande s’il existe quelque chose de similaire dans des frameworks comme TensorFlow, Keras ou PyTorch ; je me souviens avoir eu du mal à déboguer ce genre de choses dans le passé
  • Merci pour la présentation, je compte vraiment l’essayer ; je pensais être le seul à trouver pénible une syntaxe comme array[:, :, None], donc ça fait plaisir de voir que je ne suis pas le seul
  • Dans le domaine des biosignaux, NeuroPype permet au-dessus de NumPy de gérer des axes nommés pour des tenseurs n-dimensionnels et de stocker des données par élément sur chaque axe, comme les noms de canaux ou leur position
  • Ça me rappelle l’époque où NumPy dérivait des anciennes bibliothèques Numeric et Numarray ; j’imagine le camp Numarray soutenir sa position pendant vingt ans, obtenir un financement, se renommer Xarray, puis finalement battre NumPy (bien sûr, c’est en grande partie fictif)
• L’une des raisons pour lesquelles j’ai commencé à utiliser Julia, c’est que la syntaxe de NumPy était trop difficile ; en passant de MATLAB à NumPy, j’avais l’impression de devenir moins bon en programmation et de passer plus de temps à apprendre des astuces de performance qu’à faire des maths ; avec Julia, la vectorisation comme les boucles fonctionnent bien, donc on peut se concentrer uniquement sur la lisibilité du code ; j’ai retrouvé exactement cette expérience et ce ressenti dans l’article ; je ne pense pas qu’une approche « boîte noire » consistant à dire qu’il faut toujours utiliser quelque chose comme np.linalg.solve parce que c’est le plus rapide soit la bonne ; il existe aussi plusieurs bonnes raisons d’écrire soi-même des kernels spécialisés pour un problème donné
  • La raison, c’est que Julia est un langage conçu pour le calcul scientifique, alors que NumPy est une bibliothèque plaquée de force sur un langage qui n’a pas été conçu pour ça ; j’espère qu’un jour Julia l’emportera et libérera, grâce à un effet de réseau inverse, les gens qui utilisent Python
  • MATLAB aussi est aussi lent que Python si on fait tourner des boucles sans vectorisation ; la lenteur de Python est le principal problème, et même si Julia a clairement des qualités, en pratique ses usages restent très limités ; il y a eu des hacks de type JIT dans Python, mais cela reste incomplet ; on a vraiment besoin d’une alternative à Python
  • Est-ce que MATLAB est vraiment différent ? Les boucles y restent lentes, et le plus rapide reste une boîte noire entièrement optimisée comme l’opérateur \
  • Les versions modernes de Fortran aussi permettent à la fois une vectorisation rapide et des boucles rapides, donc on peut se concentrer uniquement sur la lisibilité
• Si on résume les reproches faits à numpy par rapport à Matlab et Julia, c’est que chaque fonction a ses propres paramètres liés aux axes, ses propres conventions de nommage et sa propre manière d’offrir la vectorisation, et que dès qu’on veut appliquer une fonction selon un axe donné, il faut souvent réécrire complètement le code ; l’abstraction est pourtant une base de la programmation, et NumPy la rend difficile ; dans Matlab, le code vectorisé tourne souvent tel quel ou les modifications à faire sont évidentes, alors qu’avec NumPy il faut toujours fouiller la documentation, et les ajustements de type via transpose/reshape manquent de cohérence, ce qui rend l’ensemble ambigu
  • Le support des tableaux en 3D ou plus dans Matlab est tellement faible qu’au contraire les problèmes évoqués dans l’article s’y produisent rarement
  • Pour le deuxième problème, on peut tenter vmap de JAX
  • Si l’on veut écrire une fonction pour un tableau 2x2 puis l’appliquer à une partie d’un tableau 3x2x2, on peut le faire avec des slices et squeeze, etc. ; le problème lui-même est tellement flou que j’ai du mal à le comprendre
  • On peut le traiter avec reshape
• Ce qu’il y a de plus déroutant dans numpy, c’est qu’il n’est pas clair de savoir quelles opérations sont réellement vectorisées, et qu’on ne peut pas l’indiquer explicitement comme avec la syntaxe par point de Julia ; il y a aussi beaucoup de pièges autour des types de retour ; par exemple, si on multiplie un objet poly1d P par z0 à droite, on obtient un poly1d, mais si on fait la multiplication à gauche sous la forme z0*P, on ne récupère qu’un tableau, avec une conversion de type silencieuse ; même le coefficient dominant d’un polynôme quadratique peut se récupérer de deux façons, P.coef[0] ou P[2], ce qui prête facilement à confusion ; officiellement, poly1d est une API « ancienne » et il est recommandé d’utiliser la classe Polynomial pour le nouveau code, mais en pratique il n’y a même pas d’avertissement de dépréciation ; entre ces conversions de type, les incohérences de types de données et d’autres pièges disséminés dans la bibliothèque, le débogage devient un cauchemar
• Je comprends bien ce que souligne l’auteur ; en passant de Matlab à Numpy, il y avait beaucoup d’aspects peu pratiques, et j’ai aussi trouvé le slicing des données plus pénible dans Numpy que dans Matlab ou Julia ; mais si on prend en compte le coût des licences de toolbox de Matlab, les défauts de Numpy sont compensés ; les problèmes donnés dans l’article apparaissent surtout avec des tenseurs de plus de deux dimensions, et comme Numpy était à l’origine centré sur les matrices 2D, cette limite est assez naturelle ; une bibliothèque spécialisée comme Torch est préférable, même si ce n’est pas simple non plus ; au final, la bonne synthèse ressemble à : « NumPy est un peu moins bon que les autres langages à tableaux, mais il n’y a pas vraiment grand-chose d’autre qu’on puisse utiliser »
  • NumPy visait dès le départ les tableaux à N dimensions et s’inscrivait dans la continuité de numarray, donc il ne s’est pas limité à la 2D
• Le plus gros problème de l’écosystème data science en Python, c’est que rien n’est standard ; une dizaine de bibliothèques se comportent chacune différemment, comme si elles représentaient quatre langages distincts, et la seule chose à peu près unifiée reste to_numpy() ; au final, on passe plus de temps à convertir des formats de données qu’à résoudre le problème ; Julia n’a pas que des avantages, mais l’interopérabilité entre bibliothèques pour les unités, les incertitudes, etc., y fonctionne bien, alors qu’en Python il faut toujours beaucoup de boilerplate
  • Le projet array-api essaie de standardiser les API de manipulation de tableaux dans tout l’écosystème Python
  • R est en réalité encore plus complexe à cause de ses quatre systèmes de classes
• Je me demande pourquoi les gens utilisent numpy plutôt que sage
• Certains problèmes se résolvent avec numpysane et gnuplotlib ; depuis que cette combinaison existe, j’utilise numpy activement pour toutes sortes de travaux ; sans cela, j’aurais trouvé l’outil pratiquement inutilisable
  • numpysane reste au final une boucle Python, ce n’est pas de la vraie vectorisation
  • Merci pour la présentation ; comme je râlais souvent à cause de ce genre de problèmes, je n’avais jamais pensé qu’il puisse exister une simple bibliothèque de plus haut niveau
• Pour de la vectorized multi-head attention, j’ai mis tous les produits matriciels dans einsum avec optimize="optimal" afin d’utiliser l’algorithme de multiplication chaînée de matrices pour gagner en performance ; en pratique, c’était bien environ deux fois plus rapide qu’une implémentation vectorisée classique, mais étonnamment une implémentation naïve avec des boucles est encore plus rapide ; ceux que ça intéresse peuvent regarder le code pour comprendre pourquoi ; je suppose qu’il reste des marges d’amélioration sur la cohérence de cache à l’intérieur de einsum