Broadcasting

In one line: NumPy stretches mismatched shapes to fit, without copying, so you can combine a matrix and a vector in one expression.

This is the idea that turns three nested loops into one line. It is also the source of the most confusing error messages in NumPy, so it is worth being precise about the rule.

The rule

Compare shapes right-aligned, axis by axis. Two axes are compatible if they are equal or one of them is 1. Missing leading axes are treated as 1.

        (20000, 500)     expression matrix
             (500,)      per-sample library size
   →    (20000, 500)     ✓  the (500,) is treated as (1, 500) and stretched down

        (20000, 500)
        (20000,   1)     per-gene mean, kept 2-D
   →    (20000, 500)     ✓

        (20000, 500)
          (20000,)       per-gene mean, 1-D
   →      ERROR          ✗  right-aligns as (1, 20000) vs 500 — mismatch

That last case is the one that bites. A (20000,) array right-aligns against the columns, not the rows. To broadcast it down the rows you must make it explicitly 2-D:

expr - gene_mean[:, None]        # (20000, 500) - (20000, 1)  ✓

The [:, None] idiom

v[None, :]     # (n,) → (1, n)   broadcasts across rows
v[:, None]     # (n,) → (n, 1)   broadcasts across columns

Memorise these two. Nearly every broadcasting problem is solved by inserting one of them. np.newaxis is an alias for None if you prefer the explicit name.

No memory is allocated

Broadcasting does not materialise the stretched array. NumPy sets the stride along the broadcast axis to 0, so every step along that axis re-reads the same memory. See Memory Layout and Strides.

a = np.arange(3)[:, None]        # (3, 1)
b = np.broadcast_to(a, (3, 1000))
b.nbytes                          # reports 24000...
b.base.nbytes                     # ...but only 24 bytes are real

np.broadcast_to gives you a read-only view for inspection. np.broadcast_shapes((3,1),(1,4)) tells you the result shape without computing anything — useful for debugging.

Worked example: normalizing an expression matrix

counts = ...                     # (n_genes, n_samples)

# 1. CPM: divide each column by its own library size
lib = counts.sum(axis=0)                 # (n_samples,)
cpm = counts / lib * 1e6                 # (g,s) / (s,) → right-aligns ✓

# 2. log transform
log_cpm = np.log2(cpm + 1)

# 3. centre each gene (row) on its own mean
gene_mean = log_cpm.mean(axis=1)         # (n_genes,)
centred = log_cpm - gene_mean[:, None]   # need the [:, None] ✓

# 4. full z-score per gene
gene_sd = log_cpm.std(axis=1, ddof=1)
z = centred / gene_sd[:, None]

Steps 1 and 3 are the two directions, and the asymmetry between them — one needs [:, None], one does not — is exactly the right-alignment rule in action. Reductions with keepdims=True remove the asymmetry:

gene_mean = log_cpm.mean(axis=1, keepdims=True)   # (n_genes, 1)
centred = log_cpm - gene_mean                      # no [:, None] needed

Use keepdims=True habitually. It makes broadcasting after a reduction just work.

Outer operations

Broadcasting a column against a row gives you every pair:

a = np.array([1, 2, 3])[:, None]      # (3, 1)
b = np.array([10, 20])[None, :]       # (1, 2)
a + b                                  # (3, 2) — every combination

Used for distance matrices, though for real distance work scipy.spatial.distance.pdist is faster and uses less memory:

# pairwise squared distances between n points in d dimensions
diff = X[:, None, :] - X[None, :, :]      # (n, n, d)  ← memory!
d2 = (diff ** 2).sum(axis=-1)              # (n, n)

Note the (n, n, d) intermediate. For n = 10,000 points in 50 dimensions that is 40 GB. Broadcasting is free for the operands but the result is real. This is the main way broadcasting causes out-of-memory errors.

Reading the error message

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

Right-align them:

(20000, 500)
       (20000)     ← 20000 vs 500: incompatible

The fix is almost always a [:, None] or a keepdims=True on whatever reduction produced the second shape.

More bioinformatics examples

# one-hot encode a sequence: (L,1) letters vs (1,4) bases → (L,4)
onehot = (seq_arr[:, None] == BASES[None, :]).astype(np.int8)

# score every position of a sequence against a PWM
# pwm: (4, motif_len); windows: (n_pos, motif_len) of base codes
scores = pwm[windows, np.arange(motif_len)[None, :]].sum(axis=1)

# apply per-sample scaling factors (e.g. DESeq2 size factors)
normalised = counts / size_factors                  # (g, s) / (s,)

# subtract a background profile from every sample
signal = coverage - background[:, None]

Common mistakes

• (n,) where you needed (n,1). The single most common broadcasting bug. Symptom: either a shape error, or — worse — a silently valid result when n happens to equal the other dimension. A square matrix will happily broadcast the wrong way and give you plausible garbage.
• Forgetting keepdims=True and then having to remember which axis needs [:, None].
• Materialising a huge intermediate. X[:, None] - X[None, :] is O(n²d) memory.
• Assuming broadcasting works left-aligned. It is right-aligned. Always.
• In-place ops with broadcasting. a += b requires b to broadcast into a's shape; it cannot grow a.

See also

ufuncs · Axes and Reductions · Vectorization · Memory Layout and Strides · ndarray · Applied - Expression Matrices with NumPy