HN Reader

If your arrays have more than two dimensions, please consider using Xarray [1], which adds dimension naming to NumPy arrays. Broadcasting and alignment then becomes automatic without needing to transpose, add dummy axes, or anything like that. I believe that alone solves most of the complaints in the article.

Compared to NumPy, Xarray is a little thin in certain areas like linear algebra, but since it's very easy to drop back to NumPy from Xarray, what I've done in the past is add little helper functions for any specific NumPy stuff I need that isn't already included, so I only need to understand the NumPy version of the API well enough one time to write that helper function and its tests. (To be clear, though, the majority of NumPy ufuncs are supported out of the box.)

I'll finish by saying, to contrast with the author, I don't dislike NumPy, but I do find its API and data model to be insufficient for truly multidimensional data. For me three dimensions is the threshold where using Xarray pays off.

[1] https://xarray.dev

Hear hear! Some of these complaints have been resolved with numpysane: https://github.com/dkogan/numpysane/ . With numpysane and gnuplotlib, I now find numpy acceptable and use it heavily for everything. But yeah; without these it's unusable.

> Personally, I think np.einsum is one of the rare parts of NumPy that’s actually good.

einsum only being able to do multiplication makes it quite limited. If we leaned into the Einstein notation (e.g. [1]), we could make something that is both quite nice and quite performant.

[1] https://tensora.drhagen.com/

> D = 1/(LM) np.einsum('klm,ln,km->kn', A, B, C)

The first time I came across Einsums was via the Tullio.jl package, and it seemed like magic to me. I believe the equivalent of this would be:

     @tullio D[k, n] = 1/(L*M) * A[k, l, m] * B[l, n] * C[k, m]

which is really close to the mathematical notation.

To my understanding, template strings from PEP 750 will allow for something like:

    D = 1/(L*M) * np.einsum(t'{A}[k,l,m] * {B}[l,n] * {C}[k,m]')

right? If so, that'd be pretty neat to have.

TFA keeps repeating "you can't use loops", but aren't they, like, merely less performant? I understand that there are going to be people out there doing complex algorithms (perhaps part of an ML system) where that performance is crucial and you might as well not be using NumPy in the first place if you skip any opportunities to do things in The Clever NumPy Way. But say I'm just, like, processing a video frame by frame, by using TCNW on each frame and iterating over the time dimension; surely that won't matter?

Also: TIL you can apparently use multi-dimensional NumPy arrays as NumPy array indexers, and they don't just collapse into 1-dimensional iterables. I expected `A[:,i,j,:]` not to work, or to be the same as if `j` were just `(0, 1)`. But instead, it apparently causes transposition with the previous dimension... ?

The answer to all these complaints is simple: use APL. Or rather these days, BQN.

One of the biggest hurdles of numpy is python, which is very clearly not a language designed for numerical analysis.

Comparing it to Julia and Matlab you can see many, many cases where the python design just does not fit.

Implementing array operations basically requires macros and JIT compilation. I have implemented JIT compilation of tensor operations in Scheme using LLVM. Maybe I should redo it in Clojure, but I think, most programmers don't care enough, to leave the Python ecosystem. https://wedesoft.github.io/aiscm/operation.html

Numba is a great option for speeding up (vectorizing) loops and NumPy code, apart from CuPy and JAX. Xarray is also worth trying for tensors beyond 2 dimensions.

NumPy allows a lot of science to happen. Grievance is fine but a little respect is as well.

NumPy is the lingua franca for storing and passing arrays in memory in Python.

Thank you NumPy!

I vaguely recall having a few similar complaints about Perl's PDL back in the day. Makes me wonder why we can't also support something closer to syntactically loop-esque constructs to express the same things without the need for special slicing/indexing notations with their own implicit semantics.

Here's how you do it: focus on the simple case, solve it, then ask Copilot to vectorize the code. Life is too short.

I'm with the author here. Great for simple use cases. Anything beyond that, and I find the syntax inscrutable. It's like using a terse, feature-rich DSL like Regex. Good luck deciphering someone else's numpy code (or your own from a month back) unless you really take the time to commit it to memory.

Something wrong with that font.

The author brings up some fair points. I feel like I had all sorts of small grievances transitioning from Matlab to Numpy. Slicing data still feels worse on Numpy than Matlab or Julia, but this doesn't justify the licensing costs for the Matlab stats/signal processing toolbox.

The issues that are presented in this article mostly related to tensors of rank >2. Numpy is originally just matrices so it is not surprising that it has problems in this domain. A dedicated library like Torch is certainly better. But Torch is difficult in its own ways. IDK, I guess the author's conclusion that numpy is “the worst array language other than all the other array languages” feels correct. Maybe a lack of imagination on my part.

Numpy is mostly just an interface to BLAS/Lapack, but for Python, right? BLAS/Lapack aren’t clever libraries for doing a ton of small operations, they are simple libraries for doing the easy thing (operations on big dense matrices) about as well as the hardware can.

Numpy is what it is. Seems more disappointing that he had trouble with the more flexible libraries like Jax.

Anyway there’s a clear split between the sort of functionality that Numpy specializes in and the sort that Jax does, and they don’t benefit much from stepping on each other’s toes, right?

One of the reasons why I started using Julia was because the NumPy syntax was so difficult. Going from MATLAB to NumPy I felt like I suddenly became a mediocre programmer, spending less time on math and more time on "performance engineering" of just trying to figure out how to use NumPy right. Then when I went to Julia it made sense to vectorize when it felt good and write a loop when it felt good. Because both are fast, focus on what makes the code easiest to read an understand. This blog post encapsulates exactly that experience and feeling.

Also treating things like `np.linalg.solve` as a black box that is the fastest thing in the world and you could never any better so please mangle code to call it correctly... that's just wrong. There's many reasons to build problem specific linear algebra kernels, and that's something that is inaccessible without going deeper. But that's a different story.

For the first example:

    import numpy as np
    A = np.random.random((100, 5, 5))
    b = np.random.random((100, 5, 1))
    x = np.linalg.solve(A, b)

Admittedly the documentation is a little hard to read if you just look at 'type' for b, but then it says "Returned shape is (..., M) if b is shape (M,) and (..., M, K) if b is (..., M, K) where ... is broadcasted between a and b"

So you have to make sure your b vector is actually a matrix and not a vector (if K=1 for your case).

all issues raised are addressed by jax.vmap

I would like to be able to turn off implicit broadcasting entirely. It has caused me so many evil bugs (PyTorch, but the same thing applies to Numpy).

Imagine all of the weirdness of js “truthiness” bugs (with silent unexpected behavior) combined with the inherent difficulty of debugging matrix multiplications in 2, 3, or 4 dimensions. I would rather herd goats in a Lumon basement.

Jax can do it but I’m not willing to migrate everything.

I thought I'd do something smart and inline all the matrix multiplications into the einsums of the vectorized multi-head attention implementation from the article and set optimize="optimal" to make use of the optimal matrix chain multiplication algorithm https://en.wikipedia.org/wiki/Matrix_chain_multiplication to get a nice performance boost.

    def multi_head_attention_golfed(X, W_q, W_k, W_v, W_o, optimize="optimal"):
        scores = np.einsum('si,hij,tm,hmj->hst', X, W_q, X, W_k, optimize=optimize)
        weights = softmax(W_k.shape[-1]**-0.5 * scores, axis=-1)
        projected = np.einsum('hst,ti,hiv,hvd->shd', weights, X, W_v, W_o, optimize=optimize)
        return projected.reshape(X.shape[0], W_v.shape[2])

This is indeed twice as fast as the vectorized implementation, but, disappointingly, the naive implementation with loops is even faster. Here is the code if someone wants to figure out why the performance is like that: https://pastebin.com/raw/peptFyCw

My guess is that einsum could do a better job of considering cache coherency when evaluating the sum.

It’s kind of wild how much work really smart people will do to get python to act like Fortran. This is why R is such a great language IMO. Get your data read and arrays in order in dynamic, scheme-like language, then just switch to Fortran and write actual Fortran like an adult.

Multiple dimensions (more than 2) are hard.

I was at a conference the other day, and my friend (a very smart professor) was asking me if it would be possible to move away from Xarray to Pandas or Polars...

Perhaps using Numba or Cython (with loops) might make it fast but less prone to confusion.

Luckily for me, I mostly stay in tabular data (< 3 dimensions).

Try Julia's built-in arrays or LinearAlgebra standard library handles matrix operations. You'll never go back to python (my case).

I think Numpy is representative of the Python ecosystem as a whole. Powerful, but internally complex, bloated, poorly designed and poorly documented.

I'm actually super-interested to see the next post.

TBH, if you would've asked me yesterday if I'm the sort of person who might get sucked in by a cliffhanger story about a numpy replacement, I'm pretty sure I would've been an emphatic no.

But I have, in fact, just tried random things in numpy until something worked...so, you know...tell me more...

Am I the only one who feels like this is a rant where the author is projecting their unfamiliarity with numpy? Most of the examples seem like straw men, and I can't tell if that's because they are or if the author genuinely thinks that these are all significant problems.

For example, in the linalg.solve problem, based on reading the documentation for <60 seconds I had two guesses for what would work, and from experimenting it looks like both work equally. If you don't want to take the time to understand the documentation or experiment to see what works, then just write the for loop.

For indexing problems, how about just don't use weird indexing techniques if you don't understand them? I have literally never needed to use a list of lists to index an array in my years of using numpy. And if you do need to use it for some reason, spend 1 minute to test out what happens. "What shape does B have?" is a question that can be answered with `print(B.shape)`. If the concern is about reading code written by other people, then context should make it clear what's happening, and if it's not clear then that's the fault of the person who wrote the code for not using sensible variable names/comments.

I tried to do something similar with 'first-class' dimension objects in PyTorch https://github.com/pytorch/pytorch/blob/main/functorch/dim/R... . For instance multi-head attention looks like:

    from torchdim import softmax
    def multiheadattention(q, k, v, num_attention_heads, dropout_prob, use_positional_embedding):
        batch, query_sequence, key_sequence, heads, features = dims(5)
        heads.size = num_attention_heads
    
        # binding dimensions, and unflattening the heads from the feature dimension
        q = q[batch, query_sequence, [heads, features]]
        k = k[batch, key_sequence, [heads, features]]
        v = v[batch, key_sequence, [heads, features]]
    
        # einsum-style operators to calculate scores,
        attention_scores = (q*k).sum(features) * (features.size ** -0.5)
    
        # use first-class dim to specify dimension for softmax
        attention_probs = softmax(attention_scores, dim=key_sequence)
    
        # dropout work pointwise, following Rule #1
        attention_probs = torch.nn.functional.dropout(attention_probs, p=dropout_prob)
    
        # another matrix product
        context_layer = (attention_probs*v).sum(key_sequence)
    
        # flatten heads back into features
        return context_layer.order(batch, query_sequence, [heads, features]) 
        

However, my impression trying to get a wider userbase is that while numpy-style APIs maybe are not as good as some better array language, they might not be the bottleneck for getting things done in PyTorch. However, other domains might suffer more, and I am very excited to see a better array language catch on.

Back when our lab transitioned from Matlab to Python I used numpy/scipy quite a bit. I remember heavily leaning on numpy.reshape to get things to work correctly. In some cases I did resort to looping.

Compared to Matlab (and to some extent Julia), my complaints about numpy are summed up in these two paragraphs:

> Some functions have axes arguments. Some have different versions with different names. Some have Conventions. Some have Conventions and axes arguments. And some don’t provide any vectorized version at all.

> But the biggest flaw of NumPy is this: Say you create a function that solves some problem with arrays of some given shape. Now, how do you apply it to particular dimensions of some larger arrays? The answer is: You re-write your function from scratch in a much more complex way. The basic principle of programming is abstraction—solving simple problems and then using the solutions as building blocks for more complex problems. NumPy doesn’t let you do that.

Usually when I write Matlab code, the vectorized version just works, and if there are any changes needed, they're pretty minor and intuitive. With numpy I feel like I have to look up the documentation for every single function, transposing and reshaping the array into whatever shape that particular function expects. It's not very consistent.

I don't use numpy much, but based on the documentation for that function and the fact that you can index by None to add another dimension, it seems like the correct call is:

  y = linalg.solve(A,x[:,:,None])

The documentation given says that A should be (..., M, M), and x should be (..., M, K). So if A is 100x5x5 and x is 100x5, then all you need to do is convert x to 100x5x1.

Is that right? It doesn't seem that bad.

I had to make a HN account just to show this numpy "feature".

Create an array like this np.zeros((10, 20, 30)), with shape (10, 20, 30).

Clearly if you ask for array[0][:, [0,1]].shape you will get an output of shape (20, 2), because you first remove the first dimension, and then : broadcasts across the second dimension, and you get the first two elements of the third dimension.

Also, obviously we know that in numpy array[a][b][c] == array[a,b,c]. Therefore, what do you think array[0, :, [0,1]].shape returns? (20, 2)? Nope. (2, 20).

Why? Nobody knows.

I sort of agree with the author that N>3 dimensional arrays are cumbersome in numpy, that said I think this is partly because we are really not that great thinking in higher dimensions. I'm interested what the authors solution to the problem is, but unlike the author I'm not a big fan of the eigen notation, but maybe I just don't use it often enough. I don't see the issue with a[:,:,None] notation and that's never given me issues, however I agree about the issue with index arrays. I often make something which think should work and then need to go back to the documentation to realise no that's not how it works.

The inconsistency for argument naming is also annoying (even more if we include scipy), e.g. why is it np.fft.fft(x, axis=1) but np.fft.fftshift(x, axes=1)?!

My main issue with numpy is that it's often unclear what operations will be vectorized or how they will be vectorized, and you can't be explicit about it the way you can with Julia's dot syntax.

There are also lots of gotchas related to types returned by various functions and operations.

A particularly egregious example: for a long time, the class for univariate polynomial objects was np.poly1d. It had lots of conveniences for doing usual operations on polynomials

If you multiply a poly1d object P on the right by a complex scalar z0, you get what you probably expect: a poly1d with coefficients scaled by z0.

If however you multiply P on the left by z0, you get back an array of scaled coefficients -- there's a silent type conversion happening.

So

  P*z0 # gives a poly1d
  z0*P # gives an array

I know that this is due to Python associativity rules and laissez-faire approach to datatypes, but it's fairly ugly to debug something like this!

Another fun gotcha with poly1d: if you want to access the leading coefficient for a quadratic, you can do so with either P.coef[0] or P[2]. No one will ever get these confused, right?

These particular examples aren't really fair because the numpy documentation describes poly1d as a part of the "old" polynomial API, advising new code to be written with the `Polynomial` API -- although it doesn't seem it's officially deprecated and no warnings are emitted when you use poly1d.

Anyway, there are similar warts throughout the library. Lots of gotchas having the shape of silent type conversions and inconsistent datatypes returned by the same method depending on its inputs that are downright nightmarish to debug

My main problem with numpy is that the syntax is verbose. I know from the programming language perspective it may not be considered a drawback (might even be a strength). But in practice, the syntax is a pain compared to matlab or Julia. The code for the latter is easier to read, understand, consistent with math notation.

Python in general and its ecosystem are some of the biggest roadblocks to building real-world AI systems, in my opinion.

Python is a scripting language that makes it easy to produce student-quality work. Academia adopted it for its pedagogical utility.

Unfortunately, people made the assumption that, because universities teach it, it's well-suited for production systems.

I'd argue that Python doesn't scale very well for codebases much larger than a homework assignment.

Anyone use xarray? Curious how it compares?

The trouble is that I can never tell when the cause of my issues is "numpy's syntax/documentation is bad" and when it's "I am very bad at thinking about broadcasting and shape arithmetic".

Numpy and torch have opposite conventions as to which end to broadcast from if arrays hve different numbers of dimensions, and as a result I completely forgot both. At this point I have to always None both arrays to the same dimension, which is ugly but actually fairly explicit

In Data for ML, everything has switch from NumPy (Pandas) to Arrow (Polars, DuckDB, Spark, Pandas 2.x, etc). However, Scikit-Learn is still a hold out, so it's Arrow from you data sources all to way to pre-processing pipelines in Scikit-Learn when you have to go back to NumPy. In practice, it now makes more sense to separate feature pipelines in Arrow from training pipelines with Pandas/NumPy and Scikit-Learn.*

*This is ML, not Deep Learning or Transformers.

A dumb thought: technically scipy has a “solve_banded” function that does a banded solve. He could easily recast his problem into a single big diagonal problem, I guess, just with some silly explicit zeros added. I wonder how the performance of that would compare, to iterating over a bunch of tiny solves.

Of course, it would be nice if scipy had a block diagonal solver (maybe it does?). But yea, I mean, of course it would be nice if my problem was always built in functionality of the library, lol.

Maybe a bsr_matrix could work.

Array programming requires true array languages like apl and j.

When I worked as a quant in trading risk, we used MATLAB and I always questioned the reasoning behind it... thinking it was mostly just historical plus the slick IDE, but it also makes this kind of thing dead simple. Vectorized code is just the defaullt.

I feel heard... I'm not the only one finding numpy confusing, esp with more than 3 dimensions.

You could try Elixir’s Nx and related ecosystem.

What I would really like is to type check the shape of a numpy array to prevent dimension array at runtime. You could use third party library like https://github.com/ramonhagenaars/nptyping or https://github.com/beartype/beartype#numpy-arrays but it will not extend to the methode of Numpy.

About optimization, I feel like NumPy is meant to be a de facto standard and reference implementation. It covers all use cases with decent efficiency, not the fastest way possible. There are more limited drop-in replacements that use more CPU parallelism or GPU if NumPy isn't fast enough for your use case. Just wish it were clearer which NumPy build I'm installing, cause apparently `pip3 install numpy` on my Mac gave me something built with the worst flags possible.

About >2 dimensions, I always found this confusing in NumPy but just chalked it up to >2 dim arrays being inherently confusing. Maybe there really is a better way.

Sounds like hammers and nails problem. Your nail gun is called PyTorch.

Sighs in F# and Julia

I am not a big fan of Python data libraries. They’re not cohesive in style across the board. Probably why I found R to be a better “classroom” solution. Julia is nice and so is Mathematica for purely math (and hat tip to Maple).

Why do people use numpy instead of sage?

I would much much prefer the cited "awful" syntax over the proposed loop syntax any day. Don't make me run a little virtual machine in my head to figure out what the end result of a block of code is going to be.

Since this has turned into a session of “my problem with numpy is…” then I’ll add: import time. I have some short scripts whose business logic takes no time, so importing a dependency that imports numpy turns out to take a significant amount of the script’s runtime.

That's one area LLMs hugely helped. You can just ask it and ask to generate tests to verify

"Young man, in numpy you don't understand things. You just get used to them." -- John von Neumann (probably)

I for one like Numpy's broadcast far better than Matlab's. In Numpy it's a logical operation, no unnecessary memory/copying cost.

Last time I checked Matlab (that was surely decade ago) it actually filled memory with copied data.

I skimmed the article and agree with it at a high level, though I haven't faced most of those problems personally.

I have several gripes with NumPy, or more broadly the notion of using Python to call a C/asm library that vectorizes math operations. A lot of people speak of NumPy like the solution to all your high-performance math needs in Python, which I think is disingenuous. The more custom logic you do, the less suitable the tools get. Pure Python numeric code is incredibly slow - like 1/30× compared to Java - and as you find parts that can't be vectorized, you have to drop back down to pure Python.

I would like to give the simple example of the sieve of Eratosthenes:

  def sieve_primeness(limit):
      result = [False] + [True] * limit
      for i in range(2, len(result)):
          if result[i]:
              for j in range(i * i, len(result), i):
                  result[j] = False

This code is scalar, and porting it to a language like C/C++/Rust/Java gives decent performance straight out of the box. Performance in Python is about 1/30× the speed, which is not acceptable.

People pretend that you can hand-wave the performance problem away by using NumPy. Please tell me how to vectorize this Python code. Go on, I'm waiting.

You can't vectorize the `if result[i]` because that controls whether the inner for-loop executes; so it must execute in pure Python. For the inner loop, you can vectorize that by creating a huge temporary array and then AND it with the result array, but that is extremely memory-inefficient compared to flipping bits of the result array in place, and probably messes up the cache too.

Alternatively, you can run the code in PyPy, but that usually gives a speed-up of maybe 3×, which is nowhere near enough to recover the 30× speed penalty.

Another example is that NumPy is not going to help you implement your own bigint library, because that also requires a lot of custom logic that executes between loop iterations. Meanwhile, I've implemented bigint in pure Java with acceptable performance and without issues.

Again, my point is that anytime you venture into custom numerical/logical code that is not the simple `C = A + B`, you enter a world of pain when it comes to Python performance or vectorization.

The main problem (from my point of view) of python data science ecosystem is a complete lack of standardization on anything.

You have ten different libraries that try to behave like 4 other languages and the only point of standardization is that there is usually something like .to_numpy() function. This means that most of the time I was not solving any specific problem related to data processing, but just converting data from a format one library would understand to something another library would understand.

In Julia (a language with it's own problems, of course) things mostly just work. The library for calculating uncertainties interacts well with a library handling units and all this works fine with the piloting library, or libraries solving differential equations etc. In python, this required quite a lot of boilerplate.

It seems the main complaint is about confusing shapes / dimensions. Xarray has already be mentioned, but this is a broader concept, called named dimensions, sometimes also named tensors, named axes, labeled tensors, etc, which has often been proposed before, and many implementations exists.

https://nlp.seas.harvard.edu/NamedTensor

https://namedtensor.github.io/

https://docs.pytorch.org/docs/stable/named_tensor.html

https://github.com/JuliaArrays/AxisArrays.jl

https://github.com/ofnote/tsalib

https://github.com/google-deepmind/tensor_annotations

https://github.com/google-deepmind/penzai

https://github.com/tensorflow/mesh/

https://github.com/facebookresearch/torchdim

https://xarray.dev/

https://returnn.readthedocs.io/en/latest/getting_started/dat... (that's my own development; unfortunately the doc is very much outdated on that)

I have an unpopular opinion: I don't like numpy.einsum because it is too different from the rest of numpy. You label your axes with letters but none of the other regular numpy functions do that. I usually avoid using numpy.einsum in favor of using a combination of indexing notation with numpy.newaxis (None), broadcasting, and numpy.swapaxes.

And I learned from someone more senior than me that you should instead label your variables with single-letter axes names. This way, the reader reads regular non-einsum operations and they still have the axes information in their mental model. And when you use numpy.einsum these axes labeling information become duplicated.

jax's vmap and jitting were attempts to improve things like what the author describes about solve. Pytorch had attempts to integrate this. Not sure what came out of it

I completely agree with you mate.

I have to agree with this as someone coming from a strongly-typed background (F#). PyTorch and NumPy are very flexible and powerful, but their API’s are insanely underspecified because every function seems to take multiple combinations of vaguely-typed objects. The library just figures out the right thing to do at runtime using broadcasting or other magic.

This kind of “clever” API seems to be both a benefit and curse of the Python ecosystem in general. It makes getting started very easy, but also makes mastery maddeningly difficult.

My issue with it is how easy it is to allocate all over the place if you forget to use inplace operations. It's even worse with cupy - rather than applying a series of operations to some data to produce some other data, you end up producing a set of data for each operation. Yes, there are workarounds, but they aren't as ergonomic (cupy.fuse() almost does the right thing, cleanly, but is a step you have to remember to use, and doesn't really work for anything that requires multiple shapes of array).

Coming from the excellent tidyverse or even data.table in R, Numpy always has felt like twenty steps back into mindfuck territory.