Broadcasting

Core Principle: Broadcasting automatically expands trailing and unit dimensions of mismatched ndarray operands without memory replication, enabling direct elementwise operations between multidimensional arrays.

Broadcasting replaces multi-level nested loops with concise array expressions by aligning dimension strides in compiled C routines.

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 primary source of unexpected behavior. 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.

24

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:

array([[11, 21],
       [12, 22],
       [13, 23]])

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

Exercise

Write normalise(matrix, vector), which divides the matrix by the vector along whichever axis the shapes allow. If the vector's length matches the number of columns, reshape it to (1, n) and use mode cols; otherwise if it matches the number of rows, reshape it to (n, 1) and use mode rows. When both apply, because the matrix is square, choose cols. Return ["incompatible"] when neither length matches and ["zero divisor"] when any entry of the vector is zero.

Otherwise the first line is mode=<cols|rows> shape=<r>x<c>, where the shape is the result of broadcasting the matrix shape against the reshaped divisor rather than assumed. Follow it with one line per row of the scaled matrix, values rounded to three decimals and joined by commas, and close with col_max=<list>, the largest scaled value in each column rounded the same way.