Eu odeio o NumPy

• O autor explica, com vários exemplos, os problemas do NumPy
• Operações simples com arrays são fáceis com NumPy, mas, à medida que as dimensões aumentam, a complexidade e a confusão crescem rapidamente
• O design do NumPy, incluindo broadcasting e indexação avançada, carece de clareza e abstração
• Em vez de especificar eixos explicitamente, torna-se necessário escrever código baseado em adivinhação e tentativa e erro
• O autor apresenta ideias para uma linguagem de arrays melhorada e pretende mostrar alternativas concretas no próximo texto

---

Introdução: amor e ódio pelo NumPy

• O autor afirma que usa NumPy há muito tempo, mas que ficou muito decepcionado com suas limitações
• NumPy é uma biblioteca essencial e influente para operações com arrays em Python
• Bibliotecas modernas de machine learning, como PyTorch, também têm problemas semelhantes aos do NumPy

O que é fácil e o que é difícil no NumPy

• Operações simples, como resolver equações lineares básicas, podem ser feitas com uma sintaxe clara e elegante
• Porém, quando a dimensionalidade dos arrays aumenta ou as operações ficam mais complexas, torna-se necessário processar tudo em lote sem usar loops for
• Em ambientes onde não se pode usar loops (como operações em GPU), são necessárias sintaxes de vetorização incomuns ou formas especiais de chamar funções
• No entanto, o modo exato de usar essas funções é ambíguo e difícil de entender claramente apenas pela documentação
• Na prática, é difícil para qualquer pessoa ter certeza de como usar corretamente a função linalg.solve do NumPy no caso de arrays de alta dimensão

Problemas do NumPy

• O NumPy carece de uma teoria consistente para aplicar operações a partes de arrays multidimensionais ou a eixos específicos
• Quando a dimensionalidade do array é 2 ou menos, tudo é claro; mas, a partir de 3 dimensões ou mais, fica incerto quais eixos de cada array são o alvo da operação
• Para alinhar dimensões explicitamente, ele força o uso de métodos complexos como None, broadcasting e np.tensordot
• Isso aumenta a chance de erro, reduz a legibilidade do código e eleva a possibilidade de bugs

Loops e clareza

• Na prática, se loops fossem permitidos, seria possível escrever código muito mais conciso e claro
• Código com loops pode parecer menos sofisticado, mas tem uma grande vantagem em termos de clareza
• Em contrapartida, quando a dimensionalidade do array muda, é preciso pensar manualmente em transpose ou na ordem dos eixos, o que aumenta a complexidade

np.einsum: uma função excepcionalmente boa

• np.einsum é poderosa por oferecer uma linguagem específica de domínio flexível que permite nomear os eixos
• O einsum deixa clara a intenção da operação e generaliza muito bem, permitindo implementar explicitamente operações complexas sobre eixos
• No entanto, esse tipo de suporte a operações no estilo do einsum é limitado a algumas operações e, por exemplo, não pode ser usado com linalg.solve

Problemas do broadcasting

• O broadcasting, principal truque do NumPy, é um recurso que ajusta automaticamente as dimensões quando elas não coincidem
• Em casos simples, ele é conveniente, mas na prática torna mais difícil entender claramente quais são as dimensões e gera muitos casos de erro
• Como o broadcasting é implícito, ao ler o código é preciso verificar toda vez como a operação realmente funciona

A falta de clareza na indexação

• A indexação avançada do NumPy torna muito difícil prever o shape de um array, além de ser pouco clara
• Como o shape do array resultante muda conforme a combinação de indexações, é difícil prever o resultado sem experiência prática
• A própria documentação que explica as regras de indexação é longa e complexa, exigindo muito tempo para ser assimilada
• Mesmo tentando usar apenas indexação simples, em certas operações acaba sendo inevitável recorrer à indexação avançada

Limites no design das funções do NumPy

• Muitas funções do NumPy são otimizadas apenas para shapes específicos de array
• Para arrays de alta dimensão, é necessário usar argumentos adicionais como axes, nomes de função separados ou convenções específicas, e não há consistência entre as funções
• Isso vai contra princípios básicos de programação, como abstração e reutilização
• Mesmo ao usar uma função para resolver um problema específico, para reaplicá-la a diferentes arrays e eixos é preciso reescrever o código praticamente do zero

Exemplo real: implementação de self-attention

• Ao escrever uma implementação de self-attention em NumPy, usar loops deixa o código claro, mas forçar a vetorização o torna complexo
• Quando são necessárias operações de alta dimensão, como em attention multi-head, é preciso combinar einsum e transformações de eixos, deixando o código difícil de entender

Conclusão e alternativa

• O autor afirma que o NumPy é "a única opção que se tornou tão importante no mercado, apesar de ser pior do que outras linguagens de arrays em muitos aspectos"
• Para superar vários problemas do NumPy (broadcasting, falta de clareza na indexação, inconsistência entre funções etc.), ele antecipa que criou um protótipo de uma linguagem de arrays melhorada
• O plano é apresentar as melhorias concretas (uma nova API de linguagem de arrays) em um texto separado no futuro

[youn17]

Parece a história de por que Julia surgiu. É preciso estudar as bibliotecas, mas parece uma opção realmente atraente por resolver muitos dos problemas do NumPy.

[ahwjdekf]

Se você não souber usar bem a vetorização do numpy, o desempenho vai por água abaixo. Ter que escrever levando isso em consideração é estressante e difícil.

[domino]

Acho que muitas bibliotecas Python mais antigas têm um problema parecido.

[GN⁺]

Comentários no Hacker News

• No primeiro exemplo, olhando apenas o tipo de b e a documentação, fica difícil de entender, mas como há uma explicação sobre o shape retornado, é preciso verificar se o vetor b está de fato em formato de matriz, especialmente no caso de K=1
• Se o array tiver mais de 2 dimensões, recomendo usar o Xarray, que adiciona nomes às dimensões de arrays NumPy; como ele faz broadcasting/alinhamento automaticamente sem precisar ajustar dimensões nem fazer transpose, a maioria desses problemas desaparece, o Xarray é mais fraco que o NumPy no aspecto de álgebra linear, mas é fácil voltar para o NumPy e basta criar algumas funções auxiliares, com Xarray a produtividade aumenta bastante ao lidar com dados de 3 dimensões ou mais
  • O Xarray parece uma mistura das vantagens do Pandas com o NumPy, indexações como da.sel(x=some_x).isel(t=-1).mean(["y", "z"]) são fáceis, o broadcasting também fica claro porque os nomes das dimensões são respeitados, e ele é forte no processamento de dados geoespaciais com vários CRSs, também funciona muito bem com Arviz, facilitando o tratamento de dimensões adicionais em análise bayesiana, além de permitir agrupar vários arrays em um único dataset compartilhando coordenadas em comum, então algo como ds.isel(t=-1) funciona facilmente em todos os arrays com eixo temporal
  • Graças ao Xarray, meu uso mais básico de NumPy diminuiu muito e minha produtividade aumentou bastante
  • Fico curioso se frameworks como Tensorflow, Keras e Pytorch têm algo parecido, lembro de já ter depurado com dificuldade algo como o que foi mencionado antes
  • Obrigado pela apresentação, com certeza vou experimentar, eu achava que só eu me incomodava com sintaxe como array[:, :, None], então foi bom ver alguém com a mesma opinião
  • Na área de biosinais, o NeuroPype oferece em cima do NumPy suporte a eixos nomeados para tensores de n dimensões e também permite armazenar dados por elemento em cada eixo, como nome de canal, posição etc.
  • Isso me faz lembrar da época em que o NumPy estava surgindo a partir das antigas bibliotecas Numeric e Numarray, imagino o grupo do Numarray defendendo sua posição por 20 anos, conseguindo financiamento, mudando o nome para Xarray e finalmente vencendo o NumPy (claro, isso é quase tudo ficção)
• Um dos motivos de eu ter começado a usar Julia foi justamente porque a sintaxe do NumPy era difícil demais, quando passei de MATLAB para NumPy senti que fiquei pior em programação, gastando tempo aprendendo truques de performance em vez de matemática, no Julia tanto a vetorização quanto os loops funcionam bem, então dá para focar só na legibilidade do código, senti exatamente essa experiência e esse sentimento no texto, não acho correto esse tipo de abordagem de “caixa-preta” em que algo como np.linalg.solve é tratado como a forma mais rápida e, por isso, você teria que se adaptar a ela de qualquer jeito, também há vários motivos pelos quais escrever kernels especializados para o problema pode ser melhor
  • A causa é que Julia é uma linguagem projetada para computação científica, enquanto o NumPy é uma biblioteca encaixada à força sobre uma linguagem que não foi feita para isso, espero que um dia Julia vença e liberte quem usa Python por causa do efeito de rede
  • O MATLAB também fica tão lento quanto Python se você rodar loops sem vetorização, a lentidão do Python é o maior problema, Julia claramente tem vantagens, mas na prática ainda serve só para usos bem limitados, surgiram coisas como hacks de JIT no Python, mas continuam incompletas, há uma necessidade real de uma alternativa ao Python
  • Será que o MATLAB é mesmo diferente? Os loops continuam lentos e o mais rápido ainda é uma caixa-preta totalmente otimizada como o operador '\'
  • As versões modernas de Fortran também fazem vetorização e loops rodarem rápido como no Julia, então dá para focar só na legibilidade
• Resumindo as reclamações sobre o NumPy em comparação com Matlab e Julia: cada função trata de forma diferente os argumentos relacionados a eixos, a nomenclatura e a forma de oferecer vetorização, e para aplicar uma função em determinado eixo muitas vezes é preciso reescrever completamente o código, sendo que a base da programação é abstração, algo que o NumPy dificulta, no Matlab o código vetorizado em geral roda quase do mesmo jeito ou a adaptação é clara, mas no NumPy você vive vasculhando a documentação, e ajustes de tipo com transpose/reshape etc. não são consistentes, o que torna tudo ambíguo
  • O suporte do Matlab a arrays com 3 dimensões ou mais é tão fraco que, ironicamente, os problemas citados no texto quase não aparecem lá
  • Para o segundo problema, talvez valha tentar o vmap do jax
  • Se a ideia é escrever uma função para um array 2x2 e depois aplicá-la a partes de um array 3x2x2, isso dá para fazer com slice, squeeze etc., a formulação do problema em si é tão ambígua que fica difícil até entender
  • Dá para resolver com reshape
• O ponto mais confuso do NumPy é que não fica claro quais operações funcionam de forma vetorizada, e não há como explicitar isso com uma sintaxe de ponto como no Julia, também há várias armadilhas ligadas ao tipo de retorno, por exemplo, se você multiplicar um objeto poly1d P pela direita com z0, o resultado é um poly1d, mas se fizer pela esquerda no formato z0*P, ele retorna apenas um array e ocorre uma conversão de tipo silenciosa, o coeficiente líder de uma quadrática também pode ser acessado tanto por P.coef[0] quanto por P[2], o que facilita confusão, oficialmente poly1d é uma API “antiga” e para código novo a recomendação é usar a classe Polynomial, mas na prática não há nem aviso de deprecated, por causa dessas conversões de tipo e inconsistências de datatype, a biblioteca está cheia de minas terrestres e o debug vira um pesadelo
• Concordo com os pontos levantados pelo autor, ao migrar de Matlab para Numpy senti muito incômodo, e também acho slicing de dados no Numpy mais desconfortável do que em Matlab/Julia, mas considerando o custo das licenças de toolbox do Matlab, os defeitos do Numpy acabam sendo compensados, os problemas mostrados no texto aparecem principalmente com tensores acima de 2 dimensões, e como o Numpy originalmente se baseia em matrizes (2D), essa limitação é compreensível, bibliotecas especializadas como Torch são melhores, mas também não são fáceis, no fim a sensação é: “NumPy é um pouco pior do que outras linguagens de array, mas também não há muitas alternativas utilizáveis”
  • O Numpy foi pensado desde o início para arrays N-dimensionais, como continuação do numarray, então não ficou restrito a 2D
• O maior problema do ecossistema de data science em Python é que nada é padronizado, uma dúzia de bibliotecas se comporta de formas tão diferentes quanto quatro linguagens distintas, e a única coisa mais ou menos unificada é algo como to_numpy(), no fim gasta-se mais tempo convertendo formatos de dados do que resolvendo o problema, Julia também não tem só vantagens, mas a integração entre bibliotecas para unidades, incerteza etc. funciona bem, enquanto em Python sempre é preciso escrever muito boilerplate
  • O projeto array-api está tentando padronizar a API de manipulação de arrays em todo o ecossistema Python
  • O R, por outro lado, é ainda mais complexo por causa dos seus quatro sistemas de classes
• Fico me perguntando por que as pessoas usam numpy em vez de sage
• Alguns desses problemas podem ser resolvidos usando numpysane e gnuplotlib, desde que conheci essa combinação passei a usar NumPy ativamente para todo tipo de trabalho, sem isso seria praticamente impossível de usar
  • numpysane no fim das contas ainda é loop em Python, então não é vetorização de verdade
  • Obrigado pela indicação, eu já reclamava dessas coisas de vez em quando e nunca tinha pensado que poderia existir uma biblioteca de mais alto nível tão simples
• Para fazer vectorized multi-head attention, tentei colocar todas as multiplicações de matriz em einsum e usar optimize="optimal" com o algoritmo de multiplicação em cadeia de matrizes para melhorar a performance, de fato ficou cerca de 2x mais rápido que uma implementação vetorizada comum, mas surpreendentemente uma implementação ingênua com loops ainda é mais rápida, quem quiser saber o motivo pode olhar o código, suspeito que ainda haja espaço para melhorar a coerência de cache dentro do einsum