Shapes, axes and broadcasting: the wrong-shape bug that runs without complaint

Shape is the first thing to check

m = np.arange(12).reshape(3, 4)
m.shape       # (3, 4): 3 rows, 4 columns

A shape is a tuple with one number per dimension. (3, 4) is three rows of four. (4,) is a flat vector of four. (1, 4) is a row vector — a matrix with one row — and (4, 1) is a column vector. These four shapes hold nearly the same numbers and behave differently in every operation, and a large share of NumPy bugs are one of them standing where another was expected.

reshape reinterprets the same data:

v = np.arange(6)
v.reshape(2, 3)        # 2 rows of 3
v.reshape(3, -1)       # -1 means "work it out": (3, 2)
v.reshape(-1, 1)       # column vector, (6, 1)
v[np.newaxis, :]       # row vector, (1, 6), without reshape

-1 is a placeholder for whichever size makes the element count work; only one -1 is allowed. np.newaxis (or None) in an index inserts a dimension of size 1, which is the usual way to prepare an array for broadcasting.

Axes

Reductions take an axis argument, and the rule is easy to state and easy to get backwards: the axis you name is the one that disappears.

m = np.array([[1, 2, 3],
              [4, 5, 6]])
m.sum()          # 21             everything
m.sum(axis=0)    # [5, 7, 9]      collapse rows → one value per column
m.sum(axis=1)    # [6, 15]        collapse columns → one value per row

axis=0 runs down the rows and leaves the columns; axis=1 runs across the columns and leaves the rows. If you want the mean of each column of a data matrix — the usual thing — that is axis=0. Most people reverse this at least once; the fix is to check the output shape rather than the arithmetic. A (2, 3) matrix reduced on axis=0 has shape (3,); on axis=1, (2,).

keepdims=True keeps the collapsed axis at size 1, so m.mean(axis=1, keepdims=True) has shape (2, 1) — which is exactly what you want to subtract from m in the next section.

Broadcasting: the rule

When two arrays of different shapes meet in an operation, NumPy compares their shapes from the right. Two dimensions are compatible if they are equal, or if one of them is 1. A missing dimension on the left counts as 1. Size-1 dimensions are stretched to match. If any pair is incompatible, it raises.

m = np.ones((3, 4))
m + 10                 # scalar: broadcast to (3, 4)
m + np.ones(4)         # (4,) → (1, 4) → (3, 4): adds to each row
m + np.ones((3, 1))    # (3, 1) → (3, 4): adds to each column
m + np.ones(3)         # (3,) vs (3, 4): compare 3 with 4 → error

The last one is the good case: it raises, with a message naming both shapes. The case that does not raise is the dangerous one.

The bug that runs

scores = np.random.default_rng(0).random(1000)      # (1000,)
weights = load_weights()                            # (1000, 1) — someone saved a column
weighted = scores * weights
weighted.shape                                      # (1000, 1000)

(1000,) against (1000, 1): align from the right, 1000 against 1 — compatible, stretch the 1. Then the missing left dimension of scores counts as 1 against 1000 — compatible, stretch. Result: a million-element outer product, computed without a word of complaint, and the mean of it is a plausible-looking number. Code downstream that takes .mean() or .sum() will run. The wrong answer will be confidently produced.

The check is one line: assert weighted.shape == scores.shape. Cheaper still is the habit of printing shapes as you build a pipeline. Every experienced NumPy user does this; it is not a beginner's crutch.

Where does a (1000, 1) come from? Usually from a reshape(-1, 1) that scikit-learn asked for, a CSV read with one column, or a slice that kept a dimension. Flatten it with .ravel() or .squeeze() before combining it with a flat vector — or be deliberate about which shape is the standard in your code and convert at the boundary.

Centring a matrix: broadcasting used well

The same rule that bites also does real work:

X = np.random.default_rng(0).random((500, 8))     # 500 rows, 8 features
mu = X.mean(axis=0)                               # (8,): one mean per feature
sd = X.std(axis=0)                                # (8,)
Z = (X - mu) / sd                                 # (500, 8) - (8,) → broadcasts row by row

Standardising features is one line, because (8,) broadcasts across all 500 rows. To normalise each row to unit length — the operation under cosine similarity — reduce on axis=1 and keep the dimension:

norms = np.linalg.norm(X, axis=1, keepdims=True)   # (500, 1)
X_unit = X / norms                                 # (500, 8) / (500, 1) → row by row

Without keepdims, norms is (500,), aligned from the right against 8, and it raises — the good failure. With keepdims it is (500, 1) and broadcasts down the rows correctly.

Reading the error

ValueError: operands could not be broadcast together with shapes (500,) (500,8)

That message is a diagnosis in itself: read the two shapes, align from the right, find the pair that is neither equal nor 1. Here 500 sits against 8. The fix is almost always a keepdims=True, a [:, None], or a .ravel(), and knowing which requires knowing which axis you meant.

Try this now

Build the scores and weights example, confirm the million-element result, and fix it two ways: ravel() on weights, and [:, None] on scores with a .squeeze() after. Then standardise a random (100, 5) matrix by column and check Z.mean(axis=0) is near zero and Z.std(axis=0) near one.