Odio NumPy

• El autor explica sus problemas con NumPy y señala varios inconvenientes con ejemplos
• Las operaciones simples con arreglos son fáciles en NumPy, pero cuando aumentan las dimensiones, la complejidad y la confusión crecen rápidamente
• El diseño de NumPy, incluyendo el broadcasting y la indexación avanzada, carece de claridad y abstracción
• Es casi obligatorio escribir código dependiendo de suposiciones y prueba y error, en lugar de especificar los ejes de forma explícita
• Propone ideas para un lenguaje de arreglos mejorado y presentará alternativas concretas en un próximo texto

---

Introducción: una relación de amor y odio con NumPy

• El autor cuenta que ha usado NumPy durante mucho tiempo, pero también ha quedado muy decepcionado por sus limitaciones
• NumPy es una librería esencial e influyente para operaciones con arreglos en Python
• Librerías modernas de machine learning como PyTorch también tienen problemas parecidos a los de NumPy

Lo fácil y lo difícil de NumPy

• Las operaciones simples, como resolver ecuaciones lineales básicas, pueden hacerse con una sintaxis clara y elegante
• Sin embargo, cuando la dimensionalidad del arreglo aumenta o la operación se vuelve más compleja, se vuelve necesario procesar todo en bloque sin usar for loops
• En entornos donde no se pueden usar loops libremente, como en operaciones sobre GPU, se necesitan sintaxis de vectorización extrañas o formas especiales de llamar funciones
• Pero la forma correcta de usar esas funciones es ambigua, y ni siquiera la documentación permite entenderlo con claridad
• De hecho, en el caso de arreglos de alta dimensión, es difícil que alguien tenga plena certeza de cómo usar correctamente la función linalg.solve de numpy

Los problemas de NumPy

• A NumPy le falta una teoría consistente para aplicar operaciones a partes de un arreglo multidimensional o a ejes específicos
• Cuando un arreglo tiene 2 dimensiones o menos, todo resulta claro, pero con 3 dimensiones o más ya no queda claro qué ejes son el objetivo de la operación en cada arreglo
• Obliga a usar métodos complejos para alinear dimensiones, como None, broadcasting o np.tensordot
• Estas formas hacen más probables los errores, reducen la legibilidad del código y aumentan la posibilidad de bugs

Bucles y claridad

• En la práctica, si se permitieran bucles, sería posible escribir código mucho más conciso y claro
• El código con bucles puede parecer menos elegante, pero tiene grandes ventajas en términos de claridad
• En cambio, cuando cambia la dimensionalidad del arreglo, hay que pensar manualmente en transpose o en el orden de los ejes, y la complejidad aumenta

np.einsum: una función excepcionalmente buena

• np.einsum es poderosa porque ofrece un lenguaje específico de dominio flexible que permite nombrar ejes
• einsum deja clara la intención de la operación y generaliza muy bien, por lo que permite implementar explícitamente operaciones complejas sobre ejes
• Sin embargo, el soporte para este tipo de operaciones al estilo de einsum está limitado a solo algunas operaciones; por ejemplo, no puede usarse con linalg.solve

Los problemas del broadcasting

• El broadcasting, uno de los trucos centrales de NumPy, ajusta automáticamente dimensiones cuando no coinciden
• En los casos simples es conveniente, pero en la práctica hace más difícil entender con claridad cuáles son las dimensiones y genera muchos casos propensos a error
• Como el broadcasting es implícito, cada vez que se lee el código hay que verificar de nuevo cómo funciona realmente la operación

La falta de claridad en la indexación

• La indexación avanzada de NumPy hace que sea muy difícil y poco claro predecir la shape de un arreglo
• Como la shape del arreglo resultante cambia según la combinación de índices usada, es difícil anticiparla sin experiencia práctica
• La documentación que explica las reglas de indexación también es larga y compleja, lo que implica una gran inversión de tiempo para dominarla
• Incluso si uno quiere usar solo indexación simple, en ciertas operaciones termina viéndose obligado a usar indexación avanzada

Los límites del diseño de funciones de NumPy

• Muchas funciones de NumPy están optimizadas solo para ciertas shapes de arreglos
• Para arreglos de alta dimensión hay que usar argumentos adicionales como axes, nombres de función distintos o convenciones separadas, y no existe consistencia entre funciones
• Es una estructura que va en contra de los principios de programación basados en abstracción y reutilización
• Incluso si se usa una función que resuelve un problema específico, para reaplicarla a distintos arreglos y ejes a menudo hay que reescribir el código como si fuera algo distinto

Ejemplo real: implementación de self-attention

• Al implementar self-attention en NumPy, usar bucles da claridad, pero forzar la vectorización vuelve el código complejo
• Cuando se necesitan operaciones de alta dimensión, como en multi-head attention, hay que combinar einsum y transformaciones de ejes, y el código se vuelve difícil de entender

Conclusión y alternativa

• El autor afirma que NumPy es "la única opción que se volvió importante en el mercado, aunque tiene muchos defectos frente a otros lenguajes de arreglos"
• Adelanta que creó un prototipo de un lenguaje de arreglos mejorado para superar varios problemas de NumPy, como el broadcasting, la falta de claridad en la indexación y la inconsistencia de las funciones
• Planea presentar más adelante una propuesta concreta de mejora, es decir, una nueva API de lenguaje de arreglos, en un texto aparte

[youn17]

Parece la historia de por qué nació Julia. Hay que estudiar las bibliotecas, pero me parece una opción realmente atractiva considerando que resuelve muchos de los problemas de NumPy.

[ahwjdekf]

Si no sabes usar bien la vectorización de numpy, el rendimiento se va al piso. Tener que escribirlo teniendo eso en cuenta es estresante y difícil.

[domino]

Creo que muchas bibliotecas de Python algo antiguas tienen un problema parecido.

[GN⁺]

Opiniones en Hacker News

• En el primer ejemplo, si solo miras el tipo de b y vas a la documentación, es difícil de leer, pero como sí hay una explicación sobre el shape que se devuelve, hay que confirmar si el vector b realmente tiene forma de matriz, especialmente cuando K=1
• Si el arreglo tiene más de 2 dimensiones, recomiendan usar Xarray, que agrega nombres de dimensión a los arreglos de Numpy; como el broadcasting y la alineación se hacen automáticamente sin tener que ajustar dimensiones ni hacer transpose, la mayoría de estos problemas se resuelven así. Xarray es más débil que NumPy en álgebra lineal, pero es fácil volver a NumPy y basta con hacer funciones auxiliares. Con Xarray, la productividad sube mucho al trabajar con datos de 3 dimensiones o más
  • Xarray se siente como una mezcla de las ventajas de Pandas y NumPy; facilita indexación como da.sel(x=some_x).isel(t=-1).mean(["y", "z"]), el broadcasting es claro porque se respetan los nombres de dimensión, y tiene fortalezas para procesar datos geoespaciales con varios CRS. También funciona excelente con Arviz, así que manejar dimensiones adicionales en análisis bayesiano es fácil. Además, se pueden agrupar varios arreglos en un solo dataset y compartir coordenadas comunes, por lo que es sencillo aplicar algo como ds.isel(t=-1) a todos los arreglos que tengan eje temporal
  • Gracias a Xarray, el uso básico de NumPy se reduce bastante y la productividad mejora mucho
  • Me pregunto si frameworks como Tensorflow, Keras o Pytorch tienen algo parecido; recuerdo haber tenido que depurar con mucha dificultad algo relacionado con esto hace tiempo
  • Gracias por la recomendación, pienso probarlo sí o sí. Creía que solo a mí me incomodaba una sintaxis como array[:, :, None], así que me alegra ver a alguien que piensa lo mismo
  • En el campo de biosignales, NeuroPype sobre NumPy permite ejes con nombre para tensores de n dimensiones y guardar datos por elemento en cada eje, como nombres de canales o posiciones
  • Esto me recuerda la época en que NumPy surgía a partir de las antiguas librerías Numeric y Numarray. Me imagino al bando de Numarray sosteniendo su postura durante 20 años, consiguiendo financiamiento, cambiando de nombre a Xarray y por fin derrotando a NumPy (claro, casi todo esto es ficción)
• Una de las razones por las que empecé a usar Julia fue que la sintaxis de NumPy era demasiado difícil. Al pasar de MATLAB a NumPy, sentí que programaba peor y que pasaba más tiempo aprendiendo trucos de rendimiento que haciendo matemáticas. En Julia, tanto la vectorización como los loops funcionan bien, así que uno puede concentrarse solo en la legibilidad del código. Sentí exactamente esa experiencia y esa emoción en el artículo. Tampoco creo que sea correcto ese enfoque de “caja negra” de asumir que algo como np.linalg.solve es lo máximo en velocidad y por eso hay que adaptarse a la fuerza; también hay varias razones por las que escribir kernels específicos para el problema puede ser mejor
  • La causa es que Julia es un lenguaje diseñado para cómputo científico, mientras que NumPy es una librería forzada encima de un lenguaje que no fue diseñado para eso. Ojalá algún día Julia gane y la gente que usa Python por efecto de red quede liberada
  • MATLAB también es tan lento como Python si corres loops sin vectorizar. La lentitud de Python es el problema principal. Julia claramente tiene ventajas, pero en la práctica solo puede usarse en casos bastante limitados. Han aparecido hacks tipo JIT para Python, pero siguen siendo incompletos. Hace mucha falta una alternativa a Python
  • ¿De verdad MATLAB es tan distinto? Los loops siguen siendo lentos, y lo más rápido sigue siendo una caja negra completamente optimizada como el operador \
  • Las versiones modernas de Fortran, igual que Julia, hacen que tanto la vectorización como los loops funcionen rápido, así que uno puede concentrarse solo en la legibilidad
• Resumiendo las quejas sobre numpy frente a Matlab o Julia: en cada función cambian los argumentos relacionados con ejes, los nombres y la forma en que se ofrece la vectorización, y si quieres aplicar una función sobre cierto eje tienes que reescribir completamente el código. La base de programar es la abstracción, y NumPy lo dificulta. En Matlab, el código vectorizado casi siempre funciona tal cual o las modificaciones son claras, pero en NumPy siempre hay que ir a buscar en la documentación, y el ajuste de tipos con transpose/reshape y demás no es consistente, así que todo se vuelve ambiguo
  • El soporte de Matlab para arreglos de 3 dimensiones o más es tan débil que, de hecho, los problemas mencionados en el artículo no aparecen mucho
  • Para el segundo problema, podría valer la pena probar vmap de jax
  • Si quieres escribir una función para un arreglo 2x2 y luego aplicarla a parte de un arreglo 3x2x2, eso se puede hacer con slices, squeeze y similares. El problema en sí me parece tan ambiguo que ni lo entiendo bien
  • Se puede resolver con reshape
• Lo más confuso de numpy es que no queda claro qué operaciones funcionan de forma vectorizada, y no se puede indicar explícitamente como en Julia con la sintaxis de punto. También hay muchas trampas con los tipos de retorno. Por ejemplo, si multiplicas un objeto poly1d P por z0 desde la derecha, sale un poly1d, pero si haces la multiplicación desde la izquierda en forma z0*P, solo devuelve un arreglo y la conversión de tipo ocurre silenciosamente. Incluso el coeficiente líder de una cuadrática puede obtenerse tanto con P.coef[0] como con P[2], así que es fácil confundirse. Oficialmente poly1d es una API “antigua” y para código nuevo recomiendan la clase Polynomial, pero en la práctica ni siquiera hay advertencia de deprecated. Como estas conversiones de tipo y las inconsistencias de datatype están regadas por toda la librería, depurarla es una pesadilla
• Coincido con lo que señala el autor. Al pasar de Matlab a Numpy hubo muchas incomodidades, y también siento que el slicing de datos en Numpy es más incómodo que en Matlab o Julia. Pero considerando el costo de licencias de los toolbox de Matlab, los defectos de Numpy quedan compensados. Los problemas planteados en el artículo aparecen principalmente con tensores de más de 2 dimensiones, y como Numpy originalmente estaba basado en matrices (2D), esa limitación es esperable. Una librería especializada como Torch es mejor, pero tampoco es fácil. Al final, la sensación correcta es algo como: "NumPy es un poco peor que otros lenguajes de arrays, pero tampoco hay muchas otras opciones utilizables"
  • NumPy, desde el principio, apuntaba a arreglos de N dimensiones como continuación de numarray, así que no se quedó solo en 2D
• El mayor problema del ecosistema de data science en Python es que nada es estándar. Una decena de librerías se comportan de maneras tan distintas como si fueran cuatro lenguajes diferentes, y lo único más o menos unificado es algo como to_numpy(). Al final se gasta más tiempo convirtiendo formatos de datos que resolviendo el problema. Julia tampoco tiene puras ventajas, pero la interoperabilidad entre librerías para cosas como unidades e incertidumbre funciona bien, mientras que en Python siempre hace falta mucho código boilerplate
  • El proyecto array-api está intentando estandarizar la API de manipulación de arreglos en todo el ecosistema de Python
  • R, de hecho, es todavía más complejo por sus cuatro sistemas de clases
• Me pregunto por qué la gente usa numpy en vez de sage
• Algunos problemas se resuelven usando numpysane y gnuplotlib. Desde que existe esa combinación, uso numpy de forma activa para todo tipo de trabajo; sin eso sería prácticamente imposible de usar
  • Al final, numpysane no deja de ser un loop en Python; no es vectorización real
  • Gracias por la recomendación. Como yo también me quejaba seguido por estos problemas, no se me había ocurrido que existiera una librería superior sencilla
• Para hacer vectorized multi-head attention, metí todas las multiplicaciones de matrices en einsum y usé el algoritmo de multiplicación en cadena de matrices con optimize="optimal" para mejorar el rendimiento. En efecto, llegó a ser unas 2 veces más rápido que una implementación vectorizada general, pero sorprendentemente una implementación ingenua con loops fue aún más rápida. Quien tenga curiosidad por la razón puede revisar el código. Sospecho que todavía hay margen para mejorar la cache coherency dentro de einsum