Broadcasting, rule by rule

A batch of 32 examples with 10 features each is an array of shape (32, 10). One bias per feature is an array of shape (10,). You want to add the bias to every row.

There are two obvious ways to do that and both are wrong:

for i in range(32):           # a Python loop, once per row
    x[i] = x[i] + b

x + np.tile(b, (32, 1))      # a real (32, 10) copy of the bias

And the way NumPy is actually written:

x + b                        # (32, 10) + (10,) -> (32, 10)

NumPy reads the two shapes, works out that the (10,) can stand in for every one of the 32 rows, and runs a single compiled loop over the result. There is no Python loop, and there is no tiled copy of the bias anywhere in memory. That is broadcasting.

The same rule is doing the work in almost every line of array code you will write:

x - 0.5                              a scalar against everything
x * scale                            (32, 10) * (10,)   one scale per feature
x - x.mean(axis=0)                   (32, 10) - (10,)   center each column
x / x.sum(axis=1, keepdims=True)     (32, 10) / (32, 1)  normalize each row

None of those is a convenience over a loop you could have written. Vectorized code is written in shapes, and broadcasting is the rule that says which shapes go together.

The whole of broadcasting is three steps, and they are worth memorizing in this order:

• Right-align the two shapes. Write one above the other, ends flush. NumPy never aligns from the left.
• Pad the shorter one on the left with ones until they are the same length. A missing axis is a length-one axis, not an error.
• Compare each pair. They must be equal, or one of them must be 1. The result takes the larger. Any other pair is an error.

Done by hand on (8, 1, 6) and (7, 1):

(8, 1, 6)
   (7, 1)          right-aligned

  8  1  6
  1  7  1          the shorter one padded
  -------
  8  7  6          equal, or one of them is 1, on every axis

So the result is (8, 7, 6), and notice that it is larger than either operand. Broadcasting is not "the bigger shape wins" -- it is decided one axis at a time, and each axis can be won by either side.

Padding on the left rather than the right is the part that has to be learned rather than guessed, and there is a reason for it. The last axis is the innermost one -- features, channels, columns. The leading axes are the outer ones -- batch, time, examples. Padding on the left says "this small array is one item, and it stands for all of them", which is what you almost always mean.

A scalar is the degenerate case. Its shape is (), it pads to all ones, and it matches everything.

The word "stretch" suggests that NumPy makes the small array bigger. It does not, and the previous lesson already gave you the machinery to say what it does instead.

An array is a buffer, a shape, a dtype and a strides tuple. To make a (1, 3) behave like a (4, 3), NumPy writes a new description of the same buffer in which the stride along axis 0 is zero:

b.shape   == (1, 3)      b.strides   == (24, 8)

stretched to (4, 3)      strides     == (0, 8)

move one along axis 1  ->  8 bytes    the next number
move one along axis 0  ->  0 bytes    the same three numbers again

Stepping zero bytes to reach the next row means every row reads the same three values. Nothing was allocated and nothing was copied. You can see the view directly:

np.broadcast_to(b, (4, 3))     shape (4, 3), strides (0, 8)

That view is read-only, and the reason is the stride: a single write would have to land in four places at once, so NumPy refuses rather than pick one.

What this buys you is that the size of the operands does not depend on the size of the result. x + b allocates one new array of the result shape, because the sums are new numbers and have to live somewhere. It does not allocate anything for b.

When a pair is neither equal nor contains a one, you get this, and you will get it often:

ValueError: operands could not be broadcast together with shapes (5,2) (5,)

Read what it does and does not say. It gives you the two shapes as you passed them, unpadded, and it does not tell you which axis stopped it. So you do the alignment yourself:

  5  2
  1  5          the (5,) padded
  ----
  5  ?          2 against 5, and neither is 1

The habit that finds it fastest: start at the last axis and walk left. The first pair that is neither equal nor contains a one is the culprit, and it is usually the last axis, because that is the one that got misaligned.

Two causes account for most of these:

• A rank mismatch. An axis went missing, or an extra one appeared, so the whole alignment is shifted by one and every axis is compared against the wrong partner. Check a.ndim against b.ndim before you check anything else.
• A transpose you did not do, or did twice. (5, 2) against (2, 5) is this, and it is why the numbers in the message often look like each other backwards.

With three or more operands NumPy folds the rule left to right, so the pair that clashed is not necessarily the first two -- an intermediate result can be a shape that appears nowhere in your code.

Two arrays that do not line up usually differ by a single length-one axis, and putting it in the right place is the entire fix. From the previous lesson:

m.shape == (100,)

m[:, np.newaxis]      -> (100, 1)   a column
m[np.newaxis, :]      -> (1, 100)   a row
m.reshape(-1, 1)      -> (100, 1)   the same column

Which one you want is decided by the shape you are meeting. A column repeats across the columns; a row repeats down the rows.

After a reduction, use keepdims. This is where the problem comes up most, because a reduction is exactly the thing that removes the axis you are about to need:

x.shape == (100, 5)

x.mean(axis=1)                  -> (100,)     axis 1 is gone
x.mean(axis=1, keepdims=True)   -> (100, 1)   axis 1 kept, at length one

x - x.mean(axis=1)                 error      (100, 5) against (1, 100)
x - x.mean(axis=1, keepdims=True)  works      (100, 5) against (100, 1)

keepdims=True exists for precisely this: it leaves the result lined up with the array it was computed from, so it can be subtracted from it or divided into it without any further surgery.

The pairwise grid. Deliberately making a column meet a row is how you get every pair at once, with no loop:

a.shape == (n,)

a[:, None] - a[None, :]    (n, 1) against (1, n)  ->  (n, n)

Every axis of both operands is either n or 1, so both stretch and the result is the full grid of differences. Distance matrices, attention scores, kernel and Gram matrices are all this shape, built this way.

Everything so far has failed loudly, which is the good case. The expensive bug is the one that succeeds.

y_true.shape == (100,)
y_pred.shape == (100, 1)

y_true - y_pred    ->   (100, 100)

Run the rule and it is unarguable. Right-aligned, the shapes are (1, 100) over (100, 1). Every pair contains a one, so every pair stretches, and what comes out is a hundred-by-hundred grid of every difference of every pair -- when you wanted a hundred residuals.

What happens next is the problem. ((y_true - y_pred) ** 2).mean() returns a single number, with no error and no warning. The number averages 10,000 terms of which 9,900 compare one example against a different example. It has a plausible magnitude, it goes down when you train, and it is not the loss you wrote down.

The (n, 1) is rarely something you typed. It arrives from a model whose predict returns a column, from a keepdims=True you needed one line earlier and not here, or from a DataFrame column selected as df[["y"]] rather than df["y"].

The defense is cheap, and it is the habit worth taking from this lesson:

assert y_true.shape == y_pred.shape

Where two arrays are meant to be elementwise partners, demand that their shapes are equal rather than merely broadcastable. Broadcasting is a feature you ask for on purpose; between a target and a prediction, nobody ever asked.