numpy is one idea executed ruthlessly: a typed, contiguous block of memory plus enough metadata to reinterpret it without copying. everything else — broadcasting, views, vectorised speed — falls out of that design. learn the memory model and the api stops being a list of functions and becomes a handful of consequences. 𐃏 all timings and outputs on this page are real (python 3, numpy 2.4, apple m-series).
an ndarray is four things: a buffer (one flat run of bytes), a dtype (how many bytes per element and what they mean), a shape (how to pretend the flat buffer is n-dimensional), and strides (how many bytes to jump to advance one step along each axis). element \((i, j)\) of a 2-d array lives at
\begin{equation} \text{offset}(i, j) = i \cdot s_0 + j \cdot s_1 \quad \text{bytes}, \end{equation}
concretely, for a \(3 \times 4\) array of int32 in C order: one row is \(4 \times 4 = 16\) bytes, so strides are \((16, 4)\):
import numpy as np
a = np.arange(12, dtype=np.int32).reshape(3, 4)
print(a)
print("dtype:", a.dtype, " shape:", a.shape, " strides:", a.strides)
print("itemsize:", a.itemsize, " total bytes:", a.nbytes)
# address arithmetic: element [i, j] lives at offset i*16 + j*4
i, j = 1, 2
offset = i * a.strides[0] + j * a.strides[1]
print(f"a[{i},{j}] = {a[i,j]}, byte offset = {offset}")
# transpose costs nothing: same buffer, strides swapped
print("a.T strides:", a.T.strides, " (no data moved)")
[[ 0 1 2 3]
[ 4 5 6 7]
[ 8 9 10 11]]
dtype: int32 shape: (3, 4) strides: (16, 4)
itemsize: 4 total bytes: 48
a[1,2] = 6, byte offset = 24
a.T strides: (4, 16) (no data moved)
transpose, slicing, reshape (usually), ravel (usually) are all stride tricks — new metadata over the same 48 bytes. that “usually” is the subject of the next section.
the rules, then the gotchas:
a[2:7], a[:, 0], a[::2]) returns a view — new strides, same buffer. writing to it writes to the original.v.base tells the truth: non-None means view.import numpy as np
a = np.arange(10)
v = a[2:7] # basic slice: a VIEW
v[0] = 99 # writes through to a
print("after v[0]=99, a =", a)
print("v.base is a:", v.base is a)
f = a[[2, 3, 4]] # fancy indexing: a COPY
f[0] = -1
print("after f[0]=-1, a =", a, " (unchanged)")
b = a.reshape(2, 5) # reshape of contiguous data: view
print("b.base is a:", b.base is a)
m = np.arange(6).reshape(2, 3)
t = m.T.reshape(6) # transpose is non-contiguous: this line copies
t[0] = 42
print("m unchanged by t[0]=42:", m.ravel())
after v[0]=99, a = [ 0 1 99 3 4 5 6 7 8 9]
v.base is a: True
after f[0]=-1, a = [ 0 1 99 3 4 5 6 7 8 9] (unchanged)
b.base is a: True
m unchanged by t[0]=42: [0 1 2 3 4 5]
the two directions of the gotcha: a view you thought was a copy corrupts your source data; a copy you thought was a view silently drops your writes (the t[0] = 42 line above — no error, no effect on m). when it matters, say what you mean: .copy() explicitly, or np.shares_memory(a, v) to check.
the rules, stated precisely. to combine arrays elementwise, numpy compares shapes right-aligned, axis by axis:
so \((3,1)\) with \((4,)\): pad to \((1,4)\), compare \((3,1)\) vs \((1,4)\) — both axes reconcile, result \((3,4)\). an outer product without a loop:
import numpy as np
# broadcasting: right-align shapes, each axis must match or be 1
col = np.array([[0], [10], [20]]) # shape (3, 1)
row = np.array([1, 2, 3, 4]) # shape (4,) -> (1, 4)
print("col (3,1) + row (4,) ->", (col + row).shape)
print(col + row)
# standardising a data matrix: (5,3) - (3,) works, (5,3) - (5,) does not
X = np.arange(15.0).reshape(5, 3)
mu = X.mean(axis=0) # shape (3,)
print("\nX - mu shape:", (X - mu).shape)
try:
X - X.mean(axis=1) # shape (5,) -- misaligned!
except ValueError as e:
print("X - X.mean(axis=1):", e)
print("fix with keepdims -> shape", (X - X.mean(axis=1, keepdims=True)).shape)
col (3,1) + row (4,) -> (3, 4)
[[ 1 2 3 4]
[11 12 13 14]
[21 22 23 24]]
X - mu shape: (5, 3)
X - X.mean(axis=1): operands could not be broadcast together with shapes (5,3) (5,)
fix with keepdims -> shape (5, 3)
the failed case is the canonical trap: a per-row mean has shape \((5,)\), which right-aligns against the columns axis. keepdims=True preserves the \((5,1)\) shape so the subtraction broadcasts down the correct axis.
a python loop pays interpreter overhead — type dispatch, boxing, reference counting — per element. a numpy ufunc pays it once, then runs a compiled C loop over the buffer (with simd and cache-friendly access). the difference is not subtle:
import numpy as np, timeit
n = 1_000_000
xs = list(range(n))
arr = np.arange(n, dtype=np.int64)
t_loop = timeit.timeit(lambda: sum(x * x for x in xs), number=10) / 10
t_np = timeit.timeit(lambda: np.dot(arr, arr), number=10) / 10
print(f"python loop: {t_loop*1e3:8.2f} ms")
print(f"numpy dot: {t_np*1e3:8.2f} ms")
print(f"speedup: {t_loop/t_np:8.0f}x")
python loop: 30.52 ms
numpy dot: 0.47 ms
speedup: 65x
the working rule: if you wrote for over array elements, you probably meant a ufunc, a reduction, or an einsum. loops over axes (a handful of iterations) are fine; loops over elements (millions) are the bug.
two different machines behind similar syntax:
a[[0,2,3]] gathers rows. with one index array per axis, arrays are paired elementwise — a[[0,1],[2,3]] is \([a_{0,2},\ a_{1,3}]\), not a \(2\times2\) submatrix (that wants np.ix_).a[a > 5] flattens the kept elements. as an assignment target it is the idiomatic in-place filter.import numpy as np
rng = np.random.default_rng(0)
a = rng.integers(0, 10, (4, 5))
print("fancy rows [0,2,3]: shape", a[[0, 2, 3]].shape)
print("fancy pairs a[[0,1],[2,3]] =", a[[0, 1], [2, 3]], " (elementwise pairing!)")
mask = a > 5
print("boolean mask keeps", mask.sum(), "elements ->", a[mask][:6], "...")
a[a > 8] = 8 # in-place clamp via boolean assignment
print("after clamp, max =", a.max())
fancy rows [0,2,3]: shape (3, 5)
fancy pairs a[[0,1],[2,3]] = [5 1] (elementwise pairing!)
boolean mask keeps 10 elements -> [8 6 8 6 9 6] ...
after clamp, max = 8
einsum is index notation as an api: name the axes, repeat a label to multiply along it, omit a label from the output to sum it away.
𐃏
three canonical spellings, each verified against the library primitive:
import numpy as np
rng = np.random.default_rng(42)
A, B = rng.standard_normal((3, 4)), rng.standard_normal((4, 5))
M = rng.standard_normal((4, 4))
batch = rng.standard_normal((10, 3, 4))
# 1. matmul: contract the shared index j
print("matmul ok:", np.allclose(np.einsum("ij,jk->ik", A, B), A @ B))
# 2. trace: repeated index on one operand sums the diagonal
print("trace ok:", np.isclose(np.einsum("ii->", M), np.trace(M)))
# 3. batch matmul: batch index b comes along for the ride
out = np.einsum("bij,jk->bik", batch, B)
print("batch ok:", out.shape == (10, 3, 5) and np.allclose(out, batch @ B))
# linalg essentials
K = np.array([[3.0, 1.0], [1.0, 2.0]])
b = np.array([9.0, 8.0])
x = np.linalg.solve(K, b) # never invert; solve
w, V = np.linalg.eigh(K) # symmetric -> eigh
print("solve:", x, " residual:", np.linalg.norm(K @ x - b))
print("eigenvalues:", np.round(w, 4))
# random Generator API: seeded, local, no global state
g1, g2 = np.random.default_rng(7), np.random.default_rng(7)
print("reproducible:", np.array_equal(g1.normal(size=3), g2.normal(size=3)))
matmul ok: True
trace ok: True
batch ok: True
solve: [2. 3.] residual: 0.0
eigenvalues: [1.382 3.618]
reproducible: True
linalg house rules: solve over inv (faster, numerically saner), eigh over eig for symmetric matrices, lstsq for rectangular systems, and remember @ is matmul. for randomness, the modern api is np.random.default_rng(seed) — a local, explicitly-seeded Generator — not the legacy global-state api (np.random.seed plus module-level draws), whose hidden global state makes reproducibility a whole-program property instead of a local one.
C order (default) stores rows contiguously; fortran order stores columns. the layout is invisible to semantics and very visible to the cache:
import numpy as np, timeit
t = lambda f, n: timeit.timeit(f, number=n) / n
a = np.ones((4000, 4000)) # C order: each row is 32 KB of contiguous memory
print("copy one row a[0,:] : %6.1f us" % (t(lambda: a[0, :].copy(), 100) * 1e6))
print("copy one col a[:,0] : %6.1f us" % (t(lambda: a[:, 0].copy(), 100) * 1e6))
print("ravel order='C' (view): %8.4f ms" % (t(lambda: a.ravel(order="C"), 100) * 1e3))
print("ravel order='F' (copy): %8.2f ms" % (t(lambda: a.ravel(order="F"), 10) * 1e3))
b = np.ones((2000, 2000))
bt = b.T.copy()
print("b + b.T (strided walk): %6.2f ms" % (t(lambda: b + b.T, 20) * 1e3))
print("b + bt (contiguous) : %6.2f ms" % (t(lambda: b + bt, 20) * 1e3))
copy one row a[0,:] : 0.6 us
copy one col a[:,0] : 37.8 us
ravel order='C' (view): 0.0001 ms
ravel order='F' (copy): 74.64 ms
b + b.T (strided walk): 6.48 ms
b + bt (contiguous) : 2.75 ms
reading a column of a C-order array touches one useful element per cache line — sixty times slower than the row it crosses. ravel in the array’s own order is free (a view); against the grain it is a 74 ms full copy. mixing a transposed view into arithmetic (b + b.T) more than doubles the cost of the add. the moral: match your iteration order to your storage order, and when an algorithm insists on column access, pay for one explicit np.asfortranarray up front instead of strided reads forever.
𐃏
np.linalg.solve earning its keep on KKT systemsSource: https://www.w3resource.com/python-exercises/numpy/numpy_100_exercises_with_solutions.php
Create a 1D array with values ranging from 0 to 9.
Convert a 1D array to a 2D array with 2 rows.
Multiply a 5x3 matrix by a 3x2 matrix.
Extract all odd numbers from an array of 1-10.
Replace all odd numbers in an array of 1-10 with -1.
/wiki/Knowledge is a paradox. The more one understand, the more one realises the vastness of his ignorance.
/wiki/ccs/programming/libraries/matplotlib/matplotlib has a reputation for being clunky, and the reputation is earned exactly by people fighting the wrong api. there are two: a stateful pyplot layer that mimics matlab (“draw on whatever was touched last”), and an object-oriented core in which every visible thing is an object you hold a reference to. the second is the real library; the first is sugar for one-liners. write fig, ax = plt.subplots() and stay in object land, and most of the clunk evaporates.
𐃏
all code on this page executed for real; figures were written to /tmp and are described rather than embedded, with true file sizes.
/wiki/ccs/programming/libraries/pandas/pandas is numpy with labels. a Series is a 1-d array married to an index; a DataFrame is a dict of such columns sharing one row index. the single organising idea — the one that explains both the magic and the bugs — is that every operation aligns on labels first and computes second. everything on this page runs against pandas 3.0 (outputs are real).
𐃏
/wiki/ccs/programming/libraries/