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.
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.
my experience was positive in context of pymc and arviz (as a way to represent posterior and posterior predictive distributions). API is sane and quite easy to get around.
I had to write a normal codebase in R exactly one time and it was one of the weirdest and most unpleasant coding experiences I've had. After that, I decided R is tidyverse and a handful of stats libraries. If you need more, reach for another tool.
I never understood the appeal of tidyverse. I have a colleague who, like you, thinks that R is tidyverse. He also has the nasty habit of starting Excel formulas with a '+' sign.
The main attraction of tidyverse is that it's easy to copy-paste code for common cases. If there's no ready recipe for something, it's usually just not done, or some "technical person" is called to do it.
R is used mostly as a fancy calculator, which is fine in itself, but it makes the comparisons to general purpose languages like Python quite apples-to-oranges.
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... ?
I really disagree. That's not the only point of NumPy. A lot of people use it like Matlab, to answer questions with minimal coding time, not minimal runtime.
I mean sure, the fact that it is performant means tons of functionality is built on it that is hard to find elsewhere.
But the point is still that the main purpose in building it was to be performant. To be accelerated. Even if that's not why you're personally using it.
I mean, I use my M4 Mac's Spotlight to do simple arithmetic. That's not the main point in building the M4 chip though.
As best I can see, you weren't originally talking about the reason for creating NumPy. Instead you seemed to be talking about the reason for using NumPy.
I used it once purely because I figured out that "turn a sequence of per-tile bitmap data into tiles, then produce the image corresponding to those tiles in order" is equivalent to swapping the two inner dimensions of a four-dimensional array. (And, therefore, so is the opposite operation.) The task was extremely non-performance-critical.
Of course I wasn't happy about bringing in a massive dependency just to simplify a few lines of code. Hopefully one day I'll have a slimmer alternative, perhaps one that isn't particularly concerned with optimization.
But once you need to use the GPU you need to go to another framework anyway (e.g. jax, tensorflow, arrayfire, numba...). AFAIK many of those can parallise loops using their jit functionality (in fact, e.g. numbas jit for a long time could not deal with numpy broadcasing, so you had to write out your loops). So you're not really running into a problem?
It's a slippery slope. Sometimes a python loop outside some numpy logic is fine but it's insane how much perf you can leave on the table if you overdo it.
It's not just python adding interpretor overhead, you also risk creating a lot of temporary arrays i.e. costly mallocs and memcopies.
You can draw out a sort of performance hierachy, from fastest to slowest:
* Optimized GPU code
* CPU vectorized code
* Static CPU unvectorized code
* Dynamic CPU code
where the last one refers to the fact that a language like Python, in order to add two numbers together in its native, pure-Python mode, does a lot of boxing, unboxing, resolving of class types and checking for overrides, etc.
Each of those is at least an order of magnitude slower than the next one up the hierarchy, and most of them appreciably more than one. You're probably closer to think of them as more like 1.5 orders of magnitude as a sort of back-of-the-envelope understanding.
Using NumPy incorrectly can accidentally take you from the top one, all the way to the bottom one, in one fell swoop. That can be a big deal, real quick. Or real slow, as the case may be.
In more complicated scenarios, it matters how much computation is going how far down that hierarchy. If by "processing a video frame by frame" you mean something like "I wrote a for loop on the frames but all the math is still in NumPy", you've taken "iterating on frames" from the top to the bottom, but who cares, Python can iterate on even a million things plenty quickly, especially with everything else that is going on. If, by constrast, you mean that at some point you're iterating over each pixel in pure Python, you just fell all the way down that hierarchy for each pixel and you're in bigger trouble.
In my opinionated opinion, the trouble isn't so much that it's possible to fall down that stack. That is arguably a feature, after all; surely we should have the capability of doing that sort of thing if we want. The problem is how easy it is to do without realizing it, just by using Python in what looks like perfectly sensible ways. If you aren't a systems engineer it can be hard to tell you've fallen, and even if you are honestly the docs don't make it particularly easy to figure out.
As a simple example : once upon a time a
We needed to generate a sort of heat map.
Doing it in pure python takes a few seconds at the desired size (few thousand cells where each cell needs a small formula). Dropping to numpy braucht that downs to hundreds of milliseconds. Pushing it to pure c got us to tens of milliseconds.
Yeah one of other beauties of numpy is you can pass data to/from native shared libraries compiled from C code with little overhead. This was more klidgy in Matlab last I checked
though what does "static cpu" vs "dynamic cpu" mean? it's one thing to be pointer chasing and missing the cache like OCaml can, it's another to be running a full interpreter loop to add two numbers like python does
That's what it means, basically. I draw a distinction between static code like C++ or Rust may generate and code like what Python may generate.
There is a middle ground of languages that box everything, but lack the rich complexity of Python or Ruby, such as Erlang, and I believe, O'Caml if you aren't semi-carefully programming for performance, that fits fairly cleanly into the middle ground between them. However compared to the uptake of static languages on the one side and the full dynamic scripting languages on the other, these are relatively speaking less common and don't get their own separate "in between" tier in my head, as it would end up being a .75-ish tier that would break the pattern. That is not to say they are bad or uninteresting, there's plenty of interesting languages there, they just aren't as popular.
This is great aspect of it, but that doesn’t diminish that it feels so tedious to work with compared to Julia and Matlab. Some of that is just from trying to shoehorn Python as a scientific computing language, but I think it’s time to revisit whether Python should have first party support for vectorization, arrays in memory, and broadcasting as Julia and Matlab have.
I've never understood the "so tedious" argument of Python vs Matlab. Sure, having to import numpy and use np.array to create numpy arrays is a bit typing, but other than that I don't see major differences.
Given what inconsistent mess Matlab is aside from matrix manipulation, I take NumPy any day.
Where are these people who use Matlab properly? All I ever see in other people’s Matlab code is for loops. Their idea of speeding up their code is using parfor.
Obligatory licensing rant: you have to pay extra to use parfor, as I found out one day when the license server at work told me I couldn’t run somebody else’s code because we were out of licenses for that iteration construct. This is just Mathworks taking advantage of ignorant people who don’t understand what Matlab is to inflict more misery with their insane licensing scheme. It’s not just that it’s expensive, it’s that it’s weirdly unpredictable. I mean, come on, the code was slow and chock full of errors as it was. To have it randomly fail because the license server was out of parallel toolbox licenses is just insulting.
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.
Have you heard of JIT libraries like numba (https://github.com/numba/numba)? It doesn't work for all python code, but can be helpful for the type of function you gave as an example. There's no need to rewrite anything, just add a decorator to the function. I don't really know how performance compares to C, for example.
I'm not a numpy expert, I just use it on occasion but this works and eliminates the explicit inner loop, but not the outer loop. It collects the list of prime numbers while removing multiples of the prime from the numpy array of numbers:
import numpy as np
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)
This takes advantage of the fact that True and False are treated as 1 and 0 in Python for the multiplication.
EDIT: The other solutions people have shown are closer to the sieve itself than my solution. I was having trouble with the slicing notation, now that I've seen how that should be done I'd redo the above as:
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))
However, the first version is faster by my testing so I'd end up just going back to it if I really cared about performance.
> Please tell me how to vectorize this Python code. Go on, I'm waiting.
You can't entirely get rid of the loops, but at least the inner one can be vectorized, I think. That would reduce the "pythonic runtime" from squared to linear.
Maybe something like this?
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
This is just out of my head, so may contain bugs. Also, the indices array may need to be clamped, not sure right now if numpy accepts out-of-bounds indices.
Ironically enough, I wrote a Sieve in NumPy 12 years ago or so in grad school messing around with Project Euler. I haven't touched it since and have not regularly used either NumPy or even Python in quite some time, but I still have the solutions in a private Github repo. The code is:
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]
It still relies on looping and I have no idea if you can fully vectorize everything, but it does at least get quite a bit of speedup thanks to the ability to index one array using another, so I can avoid the inner loop you have to use in pure Python.
I have absolutely no clue what that final line is doing and I'm not sure I understood it 12 years ago, either.
sieve[i] is nonzero precisely when 2i+1 is an odd prime or i=0 (because it's initialized to all-1 and nothing ever zeros the first element).
So np.nonzero(sieve) is all the indices for which i=0 or 2i+1 is an odd prime.
Except that the behaviour of np.nonzero is a bit weird: it returns a tuple of arrays, one for each axis. (So, e.g., if you feed it a 2d array it gives you (a,b) such that the nonzero elements are at positions (a[0],b[0]), (a[1],b[1]), etc.) That's kinda reasonable behaviour, but it means that for the simple case of a 1-d array you get a 1-element tuple.
So np.nonzero(sieve)[0] is actually all the indices i for which i=0 or 2i+1 is an odd prime, and now you know why I didn't write "all the indices i for which ..." above.
The case i=0 would correspond to 1 being prime, which it isn't considered to be. So we throw that one out.
So np.nonzero(sieve)[0][1::] is all the indices i for which 2i+1 is an odd prime.
And so 2*np.nonzero(sieve)[0][1::]+1 is all the 2i+1 for which 2i+1 is an odd prime; that is, all the odd primes.
I don't think I've encountered np.r_ before, but it's pretty obvious what it's doing: prepending 2 to that array, so now it's "2 followed by all the odd primes" or, more simply, all the primes.
(Obviously, in the above phrases like "all the primes" mean "all the primes up to the limit defined by the size of the array".)
> Please tell me how to vectorize this Python code. Go on, I'm waiting.
Here's a start:
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)
As you said, you can't vectorize the outer loop because each iteration depends on the result of previous iterations, so I think you can't do too much better than that with this algorithm. (There are a few easy algorithm optimizations, but I think that's kind of orthogonal to your point, and anyway you can do them in any language.)
I would expect there's still a significant performance gap between that and, say, a simple C implementation, but that shouldn't be surprising, and personally I haven't encountered these supposed droves of people claiming that NumPy fully solves the performance gap. From what I've seen there was a general consensus that vectorizing with NumPy can significantly tighten the gap but can rarely fully close it.
> As you said, you can't vectorize the outer loop because each iteration depends on the result of previous iterations
I don't see why not. That check for whether result[i] is prime is just a performance optimization. If you want to do a different performance optimization, you can do that without hurting anything.
Note that while it's true that multiples of prime numbers are composite, it's also true that multiples of composite numbers are composite.
One big problem is that Python's semantics make it difficult to implement the language efficiently. When I translate the sieve code into safe Common Lisp, I get a speedup of 15. After adding a few declarations, making the code as unsafe as C, the speedup factor is 30. Those improvements are possible because the compiler can assume many things about the functions and operators used in the code. For instance, that the Lisp equivalents of range() and len() have not been redefined.
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.
Yeah, in case anyone else has the misfortune of having to work with multi-dimensional data in MATLAB, I'd recommend the Tensor Toolbox, Tensorlab, or the N-way Toolbox.
There are certainly some aspects of it that are inelegant, in the interests of backwards compatibility, but otherwise I don't know what you are talking about. Matlab supports >2d arrays just fine, and has for at least 20 years.
> Now, how do you apply it to particular dimensions of some larger arrays?
What exactly is the complaint here? If I write a function on a 2x2 array, there's no way to apply it to tiles in a 3x2x2 array? But there is - take a slice and squeeze. What's the issue?
Edit: if the issue is some weird version of "why can't I apply my function to an arbitrary set of indices" (like the 0,1 in a 2x2x2 array) I'm at a loss for words (because that's a completely different shape - indeed it's the entire array).
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.
true, a nice jit compiler solves a lot of the problems mentioned in the article. These days i often use jax.jit for the gpu support and numpy like syntax with the added benefit of fast loop constructs.
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
There is at least one project that managed to take the Iverson ghost[1] from Numpy and actually implement a complete array language[2] on top of it, so there is no need to stray far away from Python's ecosystem.
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.
R sucks for any real programming tasks that you don't hardcode every single thing. It sucks at loading modules from your own project. It sucks to find the path of the script you're executing.
Bacially everything is stateful in R. You call standard library functions to install third party libraries ffs. And that operation can invoke your C compiler.
Putting R and a repo in a docker container to run it in a pipeline where nothing is hardcoded (unlike our datascience guy's workspace) was the worst nightmare we had to deal with.
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.
Is there anything similar to this for something like Tensorflow, Keras or Pytorch? I haven't used them super recently, but in the past I needed to do all of the things you just described in painful to debug ways.
Same - I’ve been using einops and jaxtyping together pretty extensively recently and it helps a lot for reading/writing multidimensional array code. Also array_api_compat, the API coverage isn’t perfect but it’s pretty satisfying to write code that works for both PyTorch and numpy arrays
Thanks for sharing this library. I will give it a try.
For a while I had a feeling that I was perhaps a little crazy for seeming to be only person to really dislike the use of things like ‘array[:, :, None]’ and so forth.
along those lines, for biosignal processing, NeuroPype[0] also builds on numpy and implements named axes for n-dimensional tensors, with the ability to store per-element data (i.e. channel names, positions, etc.) for each axis
Xarray is great. It marries the best of Pandas with Numpy.
Indexing like `da.sel(x=some_x).isel(t=-1).mean(["y", "z"])` makes code so easy to write and understand.
Broadcasting is never ambiguous because dimension names are respected.
It's very good for geospatial data, allowing you to work in multiple CRSs with the same underlying data.
We also use it a lot for Bayesian modeling via Arviz [1], since it makes the extra dimensions you get from sampling your posterior easy to handle.
Finally, you can wrap many arrays into datasets, with common coordinates shared across the arrays. This allows you to select `ds.isel(t=-1)` across every array that has a time dimension.
Life goes full circle sometimes. I remember that Numpy roughly came out of the amalgamation of the Numeric and Numarray libraries; I want to imagine that the Numarray people kept fighting these past 20 years to prove they were the right solution, at some point took some money from Elon Musk and renamed to Xarray [0], and finally started beating Numpy.
I understand the Julia part, but does F# have any particularly good facilities or libraries for this kind of (fast) array manipulation? It's been getting a good amount of love on HN lately, but I don't know if any specific strengths of it (beyond the general sense of "functional language with access to .NET ecosystem").
I was more curious why they chose F# among the many many languages out there, and this kind of list unfortunately is not very useful for that - similar lists exist in pretty much every awesome-$language page out there. They don't give us any idea of how robust or featureful the ecosystem and support is, which are the things that would differentiate "this too supports it, kinda" (which is every major language out there) and "this is really well-supported and ergonomic".
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.
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.
The syntax of actual array languages can be beautifully concise and expressive. You can express mathematical formulas in ways that makes sense when you read a single line, and once you get used to the notation and some originally non-intuitive quirks (from a programming background perspective), you can write in very few lines what otherwise would take you several lines of rather ugly, less readable code.
In my view, python+numpy is not actually an array language. Numpy as a library adds vectorization operations to python in order to help with speed. This is different. It does not (intend to) bring the advantages that array language syntax has, even if it was a bit more consistent.
I mean, the article and discussion are about numpy's syntax around vectorised code and problems people have with that. Many comments make comparisons with matlab, and I am pointing that a language allowing you do use arrays and write vectorised code, and a language being an array language are not the same thing. Writing vectorised code in a language were everything is based on arrays is in general more natural. The variable definitions are simpler (you do not specify what is array and what is not because everything is an array), and operations tend to work more consistently. Eg in matlab operations are by default done column-wise, because that's how the language is designed and works internally. So functions acting on 2+D arrays act by default column-wise, it does not depend on other context, and they are designed so in order to be faster, not merely consistent to the user. Consistency comes from how arrays are internally represented in memory and the need to have fast code, not just as an arbitrary design choice at the highest level.
Most developers do not touch array languages but I guess most developpers don't in general (need to) vectorise code this way and avoid loops (because they work in other problem domains, use lower level languages etc). If anything, not all problems can be vectorised anyway (or at least elegantly). But if one writes vectorised code, doing that in an array language makes more sense.
You don't have to wrap the lists in np.array if you use NumPy functions (or if one of the arguments already is a NumPy array, which usually is the case):
from numpy import *
A = [[1, 2], [3, 4]]
x = [[5], [6]]
y = dot(A, x)
That's nice, but only works, as I understand it, if you use numpy-only functions, which means that you should not use those who denote also base-pythonic, eg +,* etc operations, because then they are interpreted differently. Eg `A + x` gives
[[1, 2], [3, 4], [5], [6]]
instead of
array([[ 6, 7],[ 9, 10]])
You have to keep track of the context there to know what you can do and what not I guess, which is not ideal.
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.
"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."
If you’re going to give praise, include the GNU Fortran developers, who have really set the standard for open source production quality Fortran compiler development.
The only times I set foot in the Fortran world is when some Julia-related development [1] [2] happens, and then I peruse some surrounding links and related projects - so I have very little familiarity with GNU Fortran myself.
> 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.
Matlab is about as slow without readily vectorized operations as Python.
Slowness of Python is a huge pain point, and Julia has a clear advantage here. Sadly Julia is practically just unusable beyond quite niche purposes.
Python does now have quite serviceable jit hacks that let one escape vectorization tricks, but they are still hacks and performant Python alternative would be very welcome. Sadly there aren't any.
One of the things i suspect will happen very soon is that all languages will become as fast as each other and your reason to use one over the other for performance reasons might not exist. I know this sounds optimistic and wishful thinking but hear me out here;
Modern AIs are literal translation engines. That's their origin. The context that lets them do what they do was originally there to allow good translations. It doesn't matter if the translation is programming language output, or actual language output. I can today ask an AI to rewrite a chunk of Matlab code into Rust and it works! There's even code generators that will utilize the GPU where sensible.
We're really not that far off that we can write code in any language and transparently behind the scenes have it actually running on a more performant backend when needed. Ideally you keep the Matlab/Python frontend and the translation will be on the intermediate layers in a two way fashion so step through/debug works as expected.
Based on playing with the above steps manually with good results we're almost at the stage of just needing some glue to make it all work. Write in Matlab/Python, and run as fast as any LLVM backed language.
Sounds like transpiling Typescript to JavaScript, except that you translate to a vastly different language, thereby making you, the programmer, largely unable to reason about the performance characteristics of your code ("Can the transpiler optimize this Python loop or nah?"), and you also throw in the indeterminism of LLMs. Oooof, I'm not sure I'd want that.
When you're using Python the type of performance characteristics you're usually thinking of are "can i call an external library that runs this on the GPU" or "is this log(n) or n^2 complexity?" and such language translations shouldn't change those parts.
As in the language is regularly benchmarked at 10-100x slower than other languages. It's not a languages for caring about subtle performance characteristics outside "am i calling the external library correctly or not?" (as in this article) or "is this algorithm exponential time?".
I don’t understand this. Julia is a general-purpose language that I increasingly use to replace bash scripts. I’ve also happily used it for web scraping tasks¹, interactive web² projects³, generation of illustrations and animations, and, yes, scientific calculations. Aside from the language itself, which is fun to program in, its brilliant package system massively facilitates composing libraries into a solution.
The reason is quite simple. Julia is a language designed for scientific computation. Numpy is a library frankenstein-grafted onto a language that isn't really designed for scientific computation.
We can only hope that Julia somehow wins and those of forced to work in python because of network effects can be freed.
To be general purpose with a very large package ecosystem so that you can get just about anything started pretty quickly. It is relatively easy to do things that aren't performance critical with python, which is great if you want an MVP to grow off of or if you're messing around with something and want to make a little flask server for it, or maybe run some image recognition, or a little anything. It's really just insanely flexible and puts a lot of the cognitive load in libraries so you can get straight to doing stuff, again, often at the cost of performance.
Julia really needs to generalize and clarify its deployment story before it could possibly take off. It was built with a promise of generality but its tethered packaging is an albatross.
Is MATLAB materially different? Loops are slow (how slow depends on the version), and the fastest thing in the world is the utter perfection of the black box called '\'.
The last time I seriously used Matlab, over 10 years ago, they had implemented JIT compilation and often straightforward loops were a lot faster than trying to vectorize. And definitely less error-prone.
Iirc, ‘\’ was just shorthand for solving a system of equations, which you could mentally translate into the appropriate LAPACK function without loss of precision. You did have to be a little more careful about making extra copies than in Fortran or C (or even Python). But nothing was magic.
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.
Broadcasting is a useful and powerful feature. It's precisely specified and easily learned. However, very few type systems in real world languages are powerful enough to express the concept of broadcasting in the type system.
221 comments
[ 3.0 ms ] story [ 256 ms ] threadIt isn’t true vectorization.
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.
https://github.com/pydata/xarray/issues/3232
I've done only very simple ML stuff using R and it never feels like exactly the right tool.
R is used mostly as a fancy calculator, which is fine in itself, but it makes the comparisons to general purpose languages like Python quite apples-to-oranges.
Doing the lower level stuff that NumPy is used for is really mindfuck territory in R.
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... ?
Even code you write now, you may need to GPU accelerate later, as your simulations grow.
Falling back on loops is against the entire reason of using NumPy in the first place.
But the point is still that the main purpose in building it was to be performant. To be accelerated. Even if that's not why you're personally using it.
I mean, I use my M4 Mac's Spotlight to do simple arithmetic. That's not the main point in building the M4 chip though.
Of course I wasn't happy about bringing in a massive dependency just to simplify a few lines of code. Hopefully one day I'll have a slimmer alternative, perhaps one that isn't particularly concerned with optimization.
(And that's before you even consider GPUs.)
It's not just python adding interpretor overhead, you also risk creating a lot of temporary arrays i.e. costly mallocs and memcopies.
Each of those is at least an order of magnitude slower than the next one up the hierarchy, and most of them appreciably more than one. You're probably closer to think of them as more like 1.5 orders of magnitude as a sort of back-of-the-envelope understanding.
Using NumPy incorrectly can accidentally take you from the top one, all the way to the bottom one, in one fell swoop. That can be a big deal, real quick. Or real slow, as the case may be.
In more complicated scenarios, it matters how much computation is going how far down that hierarchy. If by "processing a video frame by frame" you mean something like "I wrote a for loop on the frames but all the math is still in NumPy", you've taken "iterating on frames" from the top to the bottom, but who cares, Python can iterate on even a million things plenty quickly, especially with everything else that is going on. If, by constrast, you mean that at some point you're iterating over each pixel in pure Python, you just fell all the way down that hierarchy for each pixel and you're in bigger trouble.
In my opinionated opinion, the trouble isn't so much that it's possible to fall down that stack. That is arguably a feature, after all; surely we should have the capability of doing that sort of thing if we want. The problem is how easy it is to do without realizing it, just by using Python in what looks like perfectly sensible ways. If you aren't a systems engineer it can be hard to tell you've fallen, and even if you are honestly the docs don't make it particularly easy to figure out.
though what does "static cpu" vs "dynamic cpu" mean? it's one thing to be pointer chasing and missing the cache like OCaml can, it's another to be running a full interpreter loop to add two numbers like python does
There is a middle ground of languages that box everything, but lack the rich complexity of Python or Ruby, such as Erlang, and I believe, O'Caml if you aren't semi-carefully programming for performance, that fits fairly cleanly into the middle ground between them. However compared to the uptake of static languages on the one side and the full dynamic scripting languages on the other, these are relatively speaking less common and don't get their own separate "in between" tier in my head, as it would end up being a .75-ish tier that would break the pattern. That is not to say they are bad or uninteresting, there's plenty of interesting languages there, they just aren't as popular.
It could be a 12 hour run vs. 12000000 hour run.
NumPy is the lingua franca for storing and passing arrays in memory in Python.
Thank you NumPy!
Given what inconsistent mess Matlab is aside from matrix manipulation, I take NumPy any day.
Obligatory licensing rant: you have to pay extra to use parfor, as I found out one day when the license server at work told me I couldn’t run somebody else’s code because we were out of licenses for that iteration construct. This is just Mathworks taking advantage of ignorant people who don’t understand what Matlab is to inflict more misery with their insane licensing scheme. It’s not just that it’s expensive, it’s that it’s weirdly unpredictable. I mean, come on, the code was slow and chock full of errors as it was. To have it randomly fail because the license server was out of parallel toolbox licenses is just insulting.
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:
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
This takes advantage of the fact that True and False are treated as 1 and 0 in Python for the multiplication.EDIT: The other solutions people have shown are closer to the sieve itself than my solution. I was having trouble with the slicing notation, now that I've seen how that should be done I'd redo the above as:
However, the first version is faster by my testing so I'd end up just going back to it if I really cared about performance.You can't entirely get rid of the loops, but at least the inner one can be vectorized, I think. That would reduce the "pythonic runtime" from squared to linear.
Maybe something like this?
This is just out of my head, so may contain bugs. Also, the indices array may need to be clamped, not sure right now if numpy accepts out-of-bounds indices.I have absolutely no clue what that final line is doing and I'm not sure I understood it 12 years ago, either.
So np.nonzero(sieve) is all the indices for which i=0 or 2i+1 is an odd prime.
Except that the behaviour of np.nonzero is a bit weird: it returns a tuple of arrays, one for each axis. (So, e.g., if you feed it a 2d array it gives you (a,b) such that the nonzero elements are at positions (a[0],b[0]), (a[1],b[1]), etc.) That's kinda reasonable behaviour, but it means that for the simple case of a 1-d array you get a 1-element tuple.
So np.nonzero(sieve)[0] is actually all the indices i for which i=0 or 2i+1 is an odd prime, and now you know why I didn't write "all the indices i for which ..." above.
The case i=0 would correspond to 1 being prime, which it isn't considered to be. So we throw that one out.
So np.nonzero(sieve)[0][1::] is all the indices i for which 2i+1 is an odd prime.
And so 2*np.nonzero(sieve)[0][1::]+1 is all the 2i+1 for which 2i+1 is an odd prime; that is, all the odd primes.
I don't think I've encountered np.r_ before, but it's pretty obvious what it's doing: prepending 2 to that array, so now it's "2 followed by all the odd primes" or, more simply, all the primes.
(Obviously, in the above phrases like "all the primes" mean "all the primes up to the limit defined by the size of the array".)
Here's a start:
As you said, you can't vectorize the outer loop because each iteration depends on the result of previous iterations, so I think you can't do too much better than that with this algorithm. (There are a few easy algorithm optimizations, but I think that's kind of orthogonal to your point, and anyway you can do them in any language.)I would expect there's still a significant performance gap between that and, say, a simple C implementation, but that shouldn't be surprising, and personally I haven't encountered these supposed droves of people claiming that NumPy fully solves the performance gap. From what I've seen there was a general consensus that vectorizing with NumPy can significantly tighten the gap but can rarely fully close it.
I don't see why not. That check for whether result[i] is prime is just a performance optimization. If you want to do a different performance optimization, you can do that without hurting anything.
Note that while it's true that multiples of prime numbers are composite, it's also true that multiples of composite numbers are composite.
> 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.
What exactly is the complaint here? If I write a function on a 2x2 array, there's no way to apply it to tiles in a 3x2x2 array? But there is - take a slice and squeeze. What's the issue?
Edit: if the issue is some weird version of "why can't I apply my function to an arbitrary set of indices" (like the 0,1 in a 2x2x2 array) I'm at a loss for words (because that's a completely different shape - indeed it's the entire array).
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
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 types returned as you described are exactly what I’d expect and idiomatic for Python. This isn’t a problem with numpy.
Previous discussion: https://news.ycombinator.com/item?id=42293723
[1] https://analyzethedatanotthedrivel.org/2018/03/31/numpy-anot...
[2] https://github.com/briangu/klongpy
But yeah, I learned Fortran to use with R lol. And it is nice. Such easy interop.
Bacially everything is stateful in R. You call standard library functions to install third party libraries ffs. And that operation can invoke your C compiler.
Putting R and a repo in a docker container to run it in a pipeline where nothing is hardcoded (unlike our datascience guy's workspace) was the worst nightmare we had to deal with.
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
The docs say that it's a prototype feature, and I think it has been that way for a few years now, so no idea how production-ready it is.
https://docs.pytorch.org/docs/stable/named_tensor.html
https://docs.kidger.site/jaxtyping/
https://data-apis.org/array-api-compat/
For a while I had a feeling that I was perhaps a little crazy for seeming to be only person to really dislike the use of things like ‘array[:, :, None]’ and so forth.
[0] https://www.neuropype.io/docs/
Indexing like `da.sel(x=some_x).isel(t=-1).mean(["y", "z"])` makes code so easy to write and understand.
Broadcasting is never ambiguous because dimension names are respected.
It's very good for geospatial data, allowing you to work in multiple CRSs with the same underlying data.
We also use it a lot for Bayesian modeling via Arviz [1], since it makes the extra dimensions you get from sampling your posterior easy to handle.
Finally, you can wrap many arrays into datasets, with common coordinates shared across the arrays. This allows you to select `ds.isel(t=-1)` across every array that has a time dimension.
[1] https://www.arviz.org/en/latest/
[0] most of the above is fiction
Also worth noting the Array API standard exists now. This is generally also trying to straighten out the sharp edges.
In my view, python+numpy is not actually an array language. Numpy as a library adds vectorization operations to python in order to help with speed. This is different. It does not (intend to) bring the advantages that array language syntax has, even if it was a bit more consistent.
Most developers do not touch array languages but I guess most developpers don't in general (need to) vectorise code this way and avoid loops (because they work in other problem domains, use lower level languages etc). If anything, not all problems can be vectorised anyway (or at least elegantly). But if one writes vectorised code, doing that in an array language makes more sense.
A = [1, 2; 3, 4];
x = [5; 6];
y = A * x;
with this uglier version:
import numpy as np
A = np.array([[1, 2], [3, 4]])
x = np.array([[5], [6]])
y = A @ x
(Perhaps you too are making a joke and it's gone over my head, in which case my apologies.)
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?
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.
This is true of modern Fortran also.
[1] https://lfortran.org/
[2] "in recent years" because that's when I became aware of this stuff, not to say there wasn't effort before then
The only times I set foot in the Fortran world is when some Julia-related development [1] [2] happens, and then I peruse some surrounding links and related projects - so I have very little familiarity with GNU Fortran myself.
[1] https://fortran-lang.discourse.group/t/update-on-prima-a-jul... [2] https://discourse.julialang.org/t/ann-julia-back-end-for-lfo...
Matlab is about as slow without readily vectorized operations as Python.
Slowness of Python is a huge pain point, and Julia has a clear advantage here. Sadly Julia is practically just unusable beyond quite niche purposes.
Python does now have quite serviceable jit hacks that let one escape vectorization tricks, but they are still hacks and performant Python alternative would be very welcome. Sadly there aren't any.
Modern AIs are literal translation engines. That's their origin. The context that lets them do what they do was originally there to allow good translations. It doesn't matter if the translation is programming language output, or actual language output. I can today ask an AI to rewrite a chunk of Matlab code into Rust and it works! There's even code generators that will utilize the GPU where sensible.
We're really not that far off that we can write code in any language and transparently behind the scenes have it actually running on a more performant backend when needed. Ideally you keep the Matlab/Python frontend and the translation will be on the intermediate layers in a two way fashion so step through/debug works as expected.
Based on playing with the above steps manually with good results we're almost at the stage of just needing some glue to make it all work. Write in Matlab/Python, and run as fast as any LLVM backed language.
As in the language is regularly benchmarked at 10-100x slower than other languages. It's not a languages for caring about subtle performance characteristics outside "am i calling the external library correctly or not?" (as in this article) or "is this algorithm exponential time?".
It'll work for toy examples. For nontrivial code the generated transpilation will flame out in a different way on every run.
Why? Just the ecosystem?
It is pretty REPL focused, but what is the thing you want that you don't have? Like better debugger support in VSCode or something?
https://news.ycombinator.com/item?id=44013455
I don’t understand this. Julia is a general-purpose language that I increasingly use to replace bash scripts. I’ve also happily used it for web scraping tasks¹, interactive web² projects³, generation of illustrations and animations, and, yes, scientific calculations. Aside from the language itself, which is fun to program in, its brilliant package system massively facilitates composing libraries into a solution.
1 https://lee-phillips.org/amazonJuliaBookRanks
2 https://www.admin-magazine.com/Archive/2025/85
3 https://lee-phillips.org/pluckit
We can only hope that Julia somehow wins and those of forced to work in python because of network effects can be freed.
Iirc, ‘\’ was just shorthand for solving a system of equations, which you could mentally translate into the appropriate LAPACK function without loss of precision. You did have to be a little more careful about making extra copies than in Fortran or C (or even Python). But nothing was 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.
Maybe it being explicit gets in the way and become inelegant once you get some experience.