Skip to content
← All articles

When to reach for NumPy, and when it makes things slower

The decision is usually framed as speed and is really about memory layout. 1,000,000 int64: a list costs 8.0 MB for the pointer array alone plus a PyObject each; array.array is 8.18 MB total; np.arange(...).nbytes is exactly 8.0 MB with zero per-element overhead. And below a few hundred elements NumPy is slower.

Ask why NumPy is faster and you will usually be told “it’s C”. That is true and it is not the useful answer, because it does not tell you when NumPy will be slower, which it frequently is.

The useful answer is memory layout.

The numbers, for 1,000,000 int64 values

representation cost
list(range(1_000_000)) 8.0 MB for the pointer array alone, plus a PyObject per element
array.array('q', ...) 8.18 MB total
np.arange(..., dtype=np.int64).nbytes exactly 8.0 MB, zero per-element overhead

A Python list of integers is an array of pointers. Each pointer is 8 bytes and points somewhere else in the heap, to an int object with a refcount, a type pointer and a variable-length digit array — roughly 28 bytes for a small value. So the list costs the pointer array plus the objects, and the objects are scattered.

Scattered is the operative word. Summing that list means following a million pointers to a million cache-unfriendly locations, and unboxing each one. Summing a NumPy array means walking 8 MB of contiguous memory in order, which the prefetcher handles perfectly and which vectorises. The speed difference is a consequence of the layout difference, and that is why it is worth stating in terms of layout: it predicts the cases where NumPy loses.

Where NumPy loses

Below roughly a few hundred elements, NumPy is slower than a list comprehension. Every NumPy operation pays a fixed per-call cost — argument parsing, dtype resolution, broadcasting rules, output allocation — on the order of microseconds. If the actual work is a hundred additions, dispatch is the whole cost. This is not in NumPy’s own documentation, which reasonably focuses on the case it is good at, and it is why “we vectorised our config validation” is a change that made things slower.

Any Python loop over a NumPy array is strictly worse than a list. This one is worth internalising because it looks like the responsible middle ground:

total = 0.0
for x in arr:          # every element is boxed into a fresh np.float64 on access
    total += x

Indexing a NumPy array from Python constructs a new scalar object for the element. You have taken the layout advantage — unboxed contiguous data — and then re-boxed every element, one at a time, with an allocation each. A plain list at least has its objects already made. If you are looping in Python, NumPy is costing you.

The rule that follows: vectorise or do not use it. There is no half-way.

The dependency-free middle ground

Two stdlib options that most people skip past:

array.array — a typed, contiguous, unboxed sequence. array('q', values) for int64, array('d', ...) for float64. It gives you the memory layout without the dependency, supports the buffer protocol, and is a good fit for “I have a million numbers and I do not need linear algebra”.

memoryview and the buffer protocol — a zero-copy view over any buffer-supporting object, so you can slice a large bytes without copying it. PEP 688 (3.12) added collections.abc.Buffer, which finally gives you a type for “supports the buffer protocol”:

from collections.abc import Buffer

def checksum(data: Buffer) -> int: ...

Before 3.12 that parameter had to be annotated bytes | bytearray | memoryview | array | ... and still missed things.

💡A function receives 50 sensor readings and computes a rolling mean. A colleague proposes converting to a NumPy array to "make it fast". What would you say, and does your answer change at 50,000 readings? click to reveal

At 50, no. The conversion itself walks the list, boxes and unboxes 50 values and allocates an array — which is very likely more work than the rolling mean. Then every NumPy call on it pays microseconds of dispatch for microseconds of arithmetic. The change adds a dependency to the function’s signature and makes it slower. If it is slow at 50 elements, the problem is not the container, it is that the function is being called a hundred thousand times.

At 50,000, still not automatically. The question is what happens next. If the array is built and then consumed by vectorised operations — a cumulative sum, a convolution, a dot product, something passed on to another array-aware library — then yes, and the win is large. If the array is built and then iterated in Python, it is worse than the list you started with, because you now box on every access.

And there is a third answer that is often the right one: the algorithm. A rolling mean computed naively is O(n·window); computed from a cumulative sum it is O(n). At 50,000 elements with a 5,000-wide window that is 2.5 × 10⁸ operations against 10⁵ — a factor of 2,500, which no amount of vectorising the wrong algorithm will match. Fix the complexity first, then decide about the container.

The typing reality

numpy.typing.NDArray[np.int64] exists, and it is less useful than it looks. dtype is not tracked through most operations: arr1 + arr2 frequently degrades to Any, and under mypy --strict — which includes --warn-return-any — that Any fires at your function boundary rather than at the operation that produced it.

The practical pattern is a thin typed façade: confine NumPy to a small module, annotate that module’s public functions with concrete return types (list[float], a dataclass, an NDArray[np.float64]), and let the rest of the codebase see only the façade. You get real checking everywhere except the fifty lines where the arrays actually live, instead of Any leaking everywhere.

💡Both a naive and a cumulative-sum rolling mean are correct. What does the cumulative-sum version give up, and how would you decide whether that matters? click to reveal

Numerical accuracy. The naive version sums each window independently, so each result has the error of summing window values. The cumulative version subtracts two prefix sums, and by the end of a long series those prefixes are large — so you are taking the difference of two large nearly-equal numbers, and catastrophic cancellation eats the low-order bits. Over 200,000 values in [0, 1) the prefix reaches ~10⁵, which costs you roughly five decimal digits of the sixteen a float64 has.

How to decide: compare the two on your actual data range and length, and compare the difference against the precision you actually need. For a rolling mean of sensor values reported to three decimal places, losing five digits of sixteen is irrelevant. For a financial cumulative that must reconcile to the cent over ten million rows, it is not.

And there is a middle option worth knowing before you conclude you must accept O(n·window): a rolling sum that adds the entering value and subtracts the leaving one is also O(n), and its error does not grow with the length of the series the way a prefix difference does — though it does drift, since errors accumulate in the running total rather than cancelling. If you need both linear time and bounded error, math.fsum per window or Kahan summation on the rolling total are the honest answers, and both cost something. The point is that “linear or accurate” is a real trade you should make deliberately, not discover in a reconciliation report.