Sixteen micro-benchmarks that made pure NumPy fast enough

Racing eleven optimisers across six Proben1 datasets, two architectures each, thirty trials per configuration, comes to well over four thousand full training runs. None of that is in this post. It’s the methodology piece that follows. What’s here is the part that happens before any of it: a SpeedTests/ directory of sixteen self-contained timeit scripts, each about thirty lines, each isolating one NumPy idiom against an alternative. No report page cites them. No commit message explains them. They exist purely so that the inner loop, the thing that runs 4,000+ times, was fast before it had to be run that many times.

Last time I mentioned in passing that the bias-as-row question was settled by exactly this kind of script, and that it never prints its answer. That turns out to be true of fifteen of the sixteen. This is the companion post: what all sixteen actually measured, what they measure now on 2026 NumPy, and, since the scripts are short enough to translate, what the same shapes look like in JavaScript, timed live in your browser below.

What the sixteen scripts compare

The 2026 column is a real re-run: uv run --with numpy on this machine, NumPy 2.5.2, Python 3.14.5, Intel i7-14700HX, one run per script, no averaging beyond what timeit itself does with its number= argument.

#	File	Compares	2018	2026 (NumPy 2.5.2)
1	SpeedTest1.py	bias appended as an input row (np.append) vs. added as a bias vector	not recorded	7.01 µs vs 4.50 µs (×1.56)
2	SpeedTest2.py	four ways to get two numbers into a sum: tuple arg, two extra unused positional args, unpack-then-call, instance attributes	not recorded	50.9 / 43.0 / 50.7 / 70.7 ns
3	SpeedTest3.py	tester.sum() vs. tester.dosum(), a cached bound-method reference	not recorded	43.4 vs 49.0 ns (×1.13)
4	SpeedTest4.py	return self.a+self.b vs. self.hold=…; return self.hold	not recorded	43.1 vs 53.1 ns (×1.23)
5	SpeedTest5.py	argument passed each call vs. stored as self.b	not recorded	42.7 vs 43.9 ns (×1.03)
6	SpeedTest6.py	pd.cut(...)*a vs. np.where(...)*a	not recorded	doesn’t run (see below)
7	SpeedTest7.py	b *= a vs. np.multiply(a, b, out=b)	not recorded	163.0 vs 159.4 µs (×1.02, confounded, see below)
8	SpeedTest8.py	a matmul wrapped in c.MatM(a, b) vs. c.MatM() with a, b pre-bound	not recorded	76.71 vs 76.08 µs (×1.01)
9	SpeedTest9.py	four ways to interleave two (100, 200) blocks into a (200, 200) array	not recorded	67.8–84.8 µs vectorised vs. 209.8 µs Python loop (×3.1)
10	SpeedTest10.py	argsort vs. argpartition vs. heapq.nsmallest for a top-1000 pick	not recorded	argsort 266.9 µs vs. argpartition 29.6 µs (×9.0), script has its own bug, see below
11	SpeedTest11.py	five manual np.roll writes vs. the same op via np.vectorize	not recorded	40.6 ms vs. 27.0 ms, single-shot (number=1)
12	SpeedTest12.py	np.where mask-assign vs. np.argmin + fancy 2-D gather (picking the fitter GA parent per weight)	not recorded	27.95 vs. 53.42 ms (×1.91)
13	SpeedTest13.py	b = c (rebind a name) vs. a[:,0] = c.flat (write a column of a 10000×100000 array)	not recorded	57.8 ns vs. 73.38 µs (×1269)
14	SpeedTest14.py	np.dot vs. np.einsum('ij, jk->ik') for a (100,1000)@(1000,300) matmul	not recorded	666.5 µs vs. 11.86 ms (×17.8), see below
15	SpeedTest15.py	a[...,:1000] vs. a[:1000]: two ways to slice a view	not recorded	293.8 vs. 223.9 ns (×1.31, both sub-microsecond)
16	SpeedTest16.py	np.average(a)*10 vs. np.sum(a)*lena	not recorded (only script that prints; output never captured)	50.6 vs. 42.0 µs (×1.21)

A few of these are worth a paragraph on their own.

The bias question, again: it’s about allocation, not arithmetic

# SpeedTests/SpeedTest1.py, lines 10–11 and 23–24
s1 = """\
np.dot(b,np.append(a,[[1]],axis=0))
"""
s2 = """\
np.dot(b,a)+c
"""

np.append doesn’t append in place: there’s no such thing for a NumPy array. It allocates a new, one-element-longer array and copies everything into it, every single forward pass. The bias-as-vector version allocates nothing extra: b and c already exist, np.dot(b,a)+c produces one new output array either way. The ×1.56 in the table is the cost of that copy for a 100-element input vector, small on its own, and exactly the sort of small cost that a training loop pays millions of times.

An honest bug in SpeedTest10.py

# SpeedTests/SpeedTest10.py, lines 27–37
s4 = """\
inds=(a[np.argpartition(a, 1000)[:1000]]).argsort()
b=c[inds]
"""
time4=timeit.timeit(setup=setup,stmt=s3, number=10000)   # <- stmt=s3, not s4

s5 = """\
heapq.nsmallest(1000,np.arange(len(a)),key=a.__getitem__)
"""
time5=timeit.timeit(setup=setup,stmt=s4, number=1000)    # <- stmt=s4, not s5

A copy-paste slip on the stmt= argument means time4 re-times s3 (so s3 gets measured twice) and time5 actually times s4’s code, not the heapq.nsmallest call s5 was written to test. The heap-based selection is never benchmarked at all. It’s a small thing. The surrounding comparison (argpartition beating a full argsort by 9× for a top-1000 pick out of 10,000) still stands, but it’s the kind of bug that a results-only report would never surface, because the report doesn’t quote this file.

SpeedTest6.py doesn’t run any more

# SpeedTests/SpeedTest6.py, lines 9–11
s1 = """\
c=pd.cut(a,[-float('inf'),500-1e-10,500+1e-10,float('inf')],labels=[-1,0,1])
b=c*a
"""

On pandas 3.0.5 (installed alongside NumPy 2.5.2 for this re-run) this raises TypeError: Object with dtype category cannot perform the numpy op multiply before it ever reaches np.where. pd.cut returns a Categorical, and modern pandas no longer lets you multiply a category by an array implicitly: arithmetic on categorical data used to silently coerce and now doesn’t. This is the interesting kind of “gap closed since 2018”: not that pandas got faster, but that it got stricter, and a script that ran fine in 2018 is a TypeError in 2026.

np.dot vs. bare np.einsum: the fix the script never tried

# SpeedTests/SpeedTest14.py, lines 12–18
s1 = """\
c=np.dot(a,b)
"""
s2 = """\
c=np.einsum('ij, jk->ik', a, b)
"""

np.dot on two 2-D float arrays dispatches straight to BLAS. np.einsum with its default arguments does not: it walks a general-purpose contraction path built for subscript expressions BLAS was never meant to evaluate, and for a plain matmul that path is markedly slower: ×17.8 in this re-run. What the script never tried is np.einsum('ij, jk->ik', a, b, optimize=True), which in the same session timed at 0.52 ms, competitive with np.dot’s 0.68 ms. The keyword has existed since long before 2018. The lesson isn’t “einsum is slow”, it’s “the general tool has an opt-in fast path, and a thirty-line speed test that never turns it on will conclude the general tool is slow.”

The real payoff: batchNeuralNetwork

The plan for this post describes the destination as “what let batchNeuralNetwork evaluate 64 GA creatures in a single dot.” Reading the class, that’s not quite what happens, and the correction is worth making, because it’s a nicer story than the approximation.

# Classes/NeuralNetworkClass.py, lines 391–399 (batchNetworkLayer.forwardPropagate)
def forwardPropagate(self, inputMatrix):
    ...
    return 1.0/(1.0+np.exp(-np.einsum('kji, kli->klj', self.WeightMatrixT, inputMatrix)-self.BiasVector))

It’s np.einsum, not np.dot, because a plain 2-D dot has no way to express “K independent matmuls, batched,” and einsum’s subscript notation does: k indexes the 64 GA creatures, j and i are a layer’s output and input dimensions, l is the row index within a batch of training examples. One call replaces what would otherwise be a Python for k in range(64) loop around 64 separate NetworkLayer objects. SpeedTest8.py found that wrapping a big matmul in a method call costs about 1% (noise, because the matmul itself dominates). The 64 creatures in batchNeuralNetwork are the opposite regime: each individual weight matrix is tiny (Classes/GeneticAlgorithmMatrix.py’s population runs on architectures like 9×4×2), so a fixed per-call dispatch cost that’s negligible against a 200×101 matmul stops being negligible when the matmul it’s attached to is nearly free. That’s the actual argument for fusing 64 calls into one: not that einsum beats dot in general (SpeedTest14 says the opposite), but that dispatch overhead matters exactly when the payload per call is small, which population-based training is.

The size of that payload has an exact formula already in the codebase:

# Classes/NeuralNetworkClass.py, lines 278–290
def getSizeOfWeightVector(layerNeurons):
    return np.sum((np.array(layerNeurons[:-1])+1)*layerNeurons[1:])

+1 per layer is the bias: one extra weight per output neuron, stored as its own array, never appended and copied. That’s the whole shape decision from SpeedTest1, generalised: the report picked the layout that doesn’t allocate, and the same layout is what makes batchNeuralNetwork a single flat self.weights array that 64 population members can share a view into.

Try it yourself

With JavaScript enabled, this becomes a table of eleven JS analogues to the SpeedTest scripts above: bias-as-row vs. bias-as-vector, method-call overhead, buffer reuse vs. fresh allocation, cache-friendly vs. naive matmul loop order, and the 64-network fused-vs-looped pass from batchNeuralNetwork, each timed live with a warm-up and a fixed time budget, plus a log-scale bar chart of the ratios. Below it, a calculator using the exact getSizeOfWeightVector formula reports weight count, memory and multiply-adds for any layer sizes, population and batch size you enter.

The benchmark runner calibrates its own batch size per experiment (so performance.now()’s own overhead never dominates a sub-100-ns operation), then takes the median of as many batches as fit in a small time budget, about four seconds total for all eleven rows. Every number in the “your browser” column is computed on the spot; nothing is precomputed or faked. The clearest result across every re-run I’ve done while building this: rebinding a JS variable (b = c) costs on the order of tens of nanoseconds, and writing one row into a column of a large typed array costs microseconds, the same shape as SpeedTest13’s ×1269, for the same underlying reason. Some things about memory don’t change between decades or languages.

What actually changed since 2018

Two things, in opposite directions. np.einsum without optimize=True is still exactly as slow relative to np.dot as the report’s era would have found it: that gap hasn’t closed, it was always avoidable with one keyword NumPy already had. And pd.cut output multiplied against an array, which worked in 2018’s pandas, now raises. Newer isn’t uniformly faster; sometimes it’s stricter about the thing you were getting away with. The sixteen scripts that produced neither number, in either direction, are the ones worth remembering: they’re evidence that someone sat down and checked, even when they never wrote the answer down.