53. NumPy: arrays, broadcasting, and why vectorized code is fast

53.1 The ndarray

import numpy as np

a = np.array([1.0, 2.0, 3.0])
M = np.arange(12, dtype=np.float32).reshape(3, 4)
print(a.dtype, a.shape, M.shape, M.ndim, M.size, M.nbytes)   # float64 (3,) (3, 4) 2 12 48
print(M)
print(M[1], M[:, 2], M[1:, ::2])          # row, column, sub-block (views, not copies)
print(np.zeros((2, 3)), np.ones(3), np.eye(2), np.linspace(0, 1, 5))
rng = np.random.default_rng(seed=0)       # the modern RNG API: always seed it
print(rng.normal(size=(2, 2)))

An ndarray is a block of raw numbers of one dtype, stored contiguously, plus a header: shape, dtype, and strides (how many bytes to jump to reach the next element along each axis).

 M = np.arange(12, dtype=np.float32).reshape(3, 4)

 memory:  [ 0  1  2  3  4  5  6  7  8  9 10 11 ]   48 contiguous bytes
 shape:   (3, 4)
 strides: (16, 4)     next row = +16 bytes, next column = +4 bytes

 M.T  -> same memory, shape (4, 3), strides (4, 16)   transpose is free: no copy

Compare with a Python list of floats: an array of pointers to separate 24-byte float objects scattered around the heap (Ch 26.3). NumPy's layout is what C, BLAS, and GPUs need, so whole-array operations run as one C (or SIMD, or GPU) loop.

53.2 Vectorization: no Python loops over numbers

import numpy as np

prices = np.array([120.0, 80.5, 99.9, 15.0])
qty = np.array([2, 1, 3, 10])

revenue = prices * qty                 # element-wise, in C
print(revenue, revenue.sum(), revenue.mean())
print(prices[prices > 90])             # boolean mask: filter without a loop
print(np.where(qty > 2, "bulk", "single"))
print(np.sqrt(prices).round(2), np.log1p(qty))

The rule from Ch 25.4: every loop over individual numbers in Python is a bug waiting for a big input. Express the computation as array operations and let NumPy run the loop.

53.3 Broadcasting

When shapes differ, NumPy stretches size-1 dimensions to match, without copying:

import numpy as np

X = np.array([[1.0, 200.0], [2.0, 300.0], [3.0, 400.0]])   # 3 samples, 2 features
mu = X.mean(axis=0)                     # shape (2,): per-feature mean
sigma = X.std(axis=0)
Z = (X - mu) / sigma                    # (3,2) - (2,) broadcasts across rows: standardization
print(Z.round(3))

# Pairwise distances between 4 points in 2-D, no loops:
P = np.random.default_rng(1).random((4, 2))
D = np.sqrt(((P[:, None, :] - P[None, :, :]) ** 2).sum(axis=-1))   # (4,1,2) - (1,4,2) -> (4,4,2)
print(D.shape, np.allclose(D, D.T))

Broadcasting rules: compare shapes from the right; two dimensions are compatible if equal or if one of them is 1. None (or np.newaxis) adds a size-1 axis. Most "shape mismatch" errors in ML code are broadcasting errors; print .shape everywhere until the shapes are second nature.

53.4 Linear algebra: the language of ML

import numpy as np

rng = np.random.default_rng(0)
X = rng.normal(size=(100, 3))                 # 100 samples, 3 features
true_w = np.array([2.0, -1.0, 0.5])
y = X @ true_w + 0.1 * rng.normal(size=100)   # @ is matrix multiplication

# Least squares: the closed-form linear regression from the Maths guide.
w, *_ = np.linalg.lstsq(X, y, rcond=None)
print(w.round(2))                             # close to [2, -1, 0.5]

# The same thing by gradient descent: the loop every neural network runs.
w = np.zeros(3)
for step in range(500):
    grad = 2 / len(X) * X.T @ (X @ w - y)     # gradient of mean squared error
    w -= 0.1 * grad
print(w.round(2))

That second loop is AI Ch 10 in five lines: predict, measure the error, compute the gradient, step downhill. Chapter 57 computes the gradient automatically instead of by hand.

53.5 Views, copies, and dtype traps

import numpy as np

a = np.arange(6)
v = a[1:4]          # a VIEW: shares memory
v[0] = 99
print(a)            # [ 0 99  2  3  4  5]  <- a changed

c = a[[1, 2, 3]]    # fancy indexing makes a COPY
c[0] = -1
print(a[1])         # still 99

x = np.array([200, 100], dtype=np.uint8)
print(x + x)        # [144 200]: silent overflow wraps around in fixed-size dtypes
print(np.array([1, 2]) / 2, np.array([1, 2]) // 2)   # float division vs integer division

Chapter 2's "names versus objects" applies again, one level down: slices are views. Fixed-size dtypes (uint8, int32, float16) overflow and lose precision silently, unlike Python's int. Image pipelines that add uint8 pixels and mixed-precision training (float16/bfloat16) both hit this.