I don't like NumPy

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

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.

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.

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

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.

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)?!

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.

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.

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.

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.

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.

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.

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.

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).

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.

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... ?

    * Optimized GPU code
    * CPU vectorized code
    * Static CPU unvectorized code
    * Dynamic CPU code

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).

> 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/

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

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.

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...

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).

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.

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".

Anyone use xarray? Curious how it compares?

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.

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.

> 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.

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.

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?

all issues raised are addressed by jax.vmap

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.

I know (and hate) that it's traditional to complain about low-contrast text here, but it really would help me read the python code if it weren't dark grey on black…

(I did submit a note in the "you misspelled" box at the bottom of the page: "you misspelled #6a6a6a as #4a4a4a and the python code is too dark :-)")

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

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

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.

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!

Sighs in F# and Julia

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.

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

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.

    import numpy as np

    def sieve(limit):
        result = np.ones(limit + 2, dtype=bool)
        result[0] = False
        result[1] = False
        for i in range(2, limit + 1):
            if result[i]:
                yield i
                result[::i] = False

    for prime in sieve(100):
        print(prime)

  def sieve_primes(limit):
    nums = np.arange(limit + 1)
    nums[1] = 0
    primes = []
    for i in range(2, limit):
      if nums[i] != 0: # could just be `if nums[i]:` since 0 is false-y
        primes.append(i)
        nums = nums * (np.mod(nums, i) != 0)
    print(primes)
  sieve_primes(1000)

  def sieve_primes(limit):
    nums = np.arange(limit + 1)
    nums[1] = 0
    for i in range(2, limit):
      if nums[i] != 0:
        nums[i+i::i] = 0
    print(np.nonzero(nums))

  result = np.full((limit+1,), True, dtype=bool)
  result[0] = False
  indices = np.arange(len(result))
  for i in range(2, len(result)):
    result[indices[i:] * i] = False

  def primes(n):
      """Input: A natural number n.
         Output: An array of primes, p < n."""
      assert n >= 2
      sieve = np.ones(n / 2, dtype=np.bool)
      for i in range(3, int(n ** 0.5) + 1, 2):
          if sieve[i / 2]:
              sieve[i * i / 2::i] = False
      return np.r_[2, 2 * np.nonzero(sieve)[0][1::] + 1]

extend this to Python itself, and my opinion is revealed.

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

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.

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