Broadcasting — The Rule That Runs the World

Course map · DL003 · Broadcasting, allocation and intended sharing

The story

A shape error is helpful: it announces that dimensions disagree. A legal but unintended broadcast is subtler. Predictions(N,) compared with labels(N,1) form an(N,N) all-pairs table. Its mean is a perfectly valid number for the wrong question.

We will derive the indexing rule, predict each original shape exercise, inspect zero-stride sharing, compare pairwise-distance algorithms and implement five attention/normalization patterns. None of this assumes a universal speedup or that a valid shape guarantees the right semantics.

Prerequisites

DL001's shape, dtype, slicing, singleton axes and reductions; DL002's matrix contraction. Here N and K count observations and centers, d counts features. For X:(N,d) and C:(K,d), X[:,None,:]-C[None,:,:] has(N,K,d). That result requires those stated operand shapes; the bare syntax alone cannot determine it.

1 · What is broadcasting, in one sentence?

For compatible elementwise operands, broadcasting reuses values along missing or singleton axes rather than first materializing repeated input copies. The result can still allocate. Think of a reusable stamp applied to many rows: one stamp does not mean the decorated sheet is free.

For X with rows(1,2,3),(4,5,6) and b=(10,20,30), X+b gives rows(11,22,33),(14,25,36). np.tile(b,(2,1)) explicitly creates the repeated b; the broadcasted addition need not. Both results need storage unless an appropriate output buffer is reused. Timing and temporary allocation depend on the expression and implementation.

2 · The three rules

1. Align dimensions from the right; missing leading axes behave as size1.
2. At each position, dimensions must be equal or one must be1.
3. Use the non-singleton size (or the equal size) as the result dimension.

Right alignment is the documented API convention, not a mathematical consequence of C-order memory. Broadcasting works on transposed and Fortran-order arrays too. Feature-last layouts make a(d,) bias convenient; channel-first image APIs demonstrate that “features are universally last” is false.

2.1 Walking through Rule 1 (right-alignment)

Write(3,) as(1,3) beside(2,3). Padding is conceptual metadata, not inserting new numerical values. (2,4) against (2,3,4) pads to (1,2,4), not (2,1,4); an explicit singleton is needed for the latter intent.

2.2 Walking through Rule 2 (compatibility)

For(1,3) and(2,3),1 versus2 is compatible and3 versus3 is equal. For(3,4) and(2,4),3 versus2 fails even though the final dimensions agree. A scalar shape() has no constraints of its own and can broadcast to any array shape, subject to dtype/operation validity.

2.3 Walking through Rule 3 (stretching)

The result of(1,3) and(2,3) is(2,3). Avoid the oversimplified “take the maximum”: dimension0 and1 are compatible and produce 0, not1. Thus(0,3)+(1,3) is an empty(0,3) result, while0 and2 are incompatible. An empty dimension cannot be filled with invented observations.

3 · Worked examples — do these on paper

3.1 Adding a scalar to a matrix

() against(3,4) acts like(1,1) against(3,4), giving(3,4). np.zeros((3,4))+5 is twelve fives. The scalar is reused; a newly returned result still stores twelve entries.

3.2 Row vector + matrix

(3,) pads to(1,3), so(2,3)+(3,) adds the same three offsets to both rows. A one-dimensional vector does not intrinsically have row orientation; right alignment creates this particular effect.

3.3 Column vector + row vector = matrix (the WEIRD one)

(3,1)+(1,4) gives(3,4). For columns1,2,3 and row10,20,30,40, rows are(11,21,31,41),(12,22,32,42),(13,23,33,43). This is an outer sum, not an outer product; replacing+ by* produces the latter.

3.4 The one that FAILS

(3,4)+(2,4) raises ValueError because3 and2 are neither equal nor singleton. Do not blindly reshape until the exception disappears; first ask whether you meant aligned observations or all-pairs combinations.

Original seven predictions, with answers:

1. (5,)+(3,5)→(3,5): the leading1 expands to3.
2. (5,1)+(1,4)→(5,4): each operand is reused along the other's non-singleton axis.
3. (2,3,4)+(4,)→(2,3,4): pad to(1,1,4).
4. (2,3,4)+(3,4)→(2,3,4): pad to(1,3,4).
5. (2,3,4)+(2,1,4)→(2,3,4): the middle1 expands to3.
6. (2,3,4)+(2,4)→error: padded middle2 conflicts with3.
7. (4,3)+(2,3)→error:4 conflicts with2.

import numpy as np
 
def broadcast_shape(a,b):
    n=max(len(a),len(b));a=(1,)*(n-len(a))+tuple(a);b=(1,)*(n-len(b))+tuple(b)
    out=[]
    for x,y in zip(a,b):
        if x==y: out.append(x)
        elif x==1: out.append(y)
        elif y==1: out.append(x)
        else: raise ValueError('incompatible dimensions')
    return tuple(out)
cases=[((5,),(3,5),(3,5)),((5,1),(1,4),(5,4)),((2,3,4),(4,),(2,3,4)),
       ((2,3,4),(3,4),(2,3,4)),((2,3,4),(2,1,4),(2,3,4)),
       ((),(3,4),(3,4)),((0,3),(1,3),(0,3))]
for a,b,expected in cases:
    assert broadcast_shape(a,b)==np.broadcast_shapes(a,b)==expected
for a,b in [((2,3,4),(2,4)),((4,3),(2,3)),((0,),(2,))]:
    for fn in (broadcast_shape,np.broadcast_shapes):
        try: fn(a,b)
        except ValueError: pass
        else: raise AssertionError('incompatible shapes accepted')
X=np.array([[1,2,3],[4,5,6]]);b=np.array([10,20,30])
np.testing.assert_array_equal(X+b,[[11,22,33],[14,25,36]])
np.testing.assert_array_equal(np.array([1,2,3])[:,None]+np.array([10,20,30,40]),
                              [[11,21,31,41],[12,22,32,42],[13,23,33,43]])
print('shape algorithm, original exercises and zero-axis controls passed')

4 · Constructing broadcasts on purpose — the [:, None] idiom

4.1 The problem

For five3-D points, define $D_{ij}=\|X_i-X_j\|_2$. A double loop is a readable reference, not inherently ugly. It makes observation identity explicit and is a useful small oracle.

4.2 The broadcasting solution

X[:,None,:] has(5,1,3); X[None,:,:] has(1,5,3). Subtraction produces(5,5,3), with the feature differences at each pair. Squaring and summing the last axis produces squared distances(5,5), then square root gives D. The result is symmetric with zero diagonal up to arithmetic, but these properties alone do not prove every off-diagonal value.

For N points and K centers the temporary has NKd entries. N=K=1000,d=128,float64 requires1,024,000,000 bytes for one differences array, before squares and output. Broadcasting inputs is cheap; this output is not.

4.3 Even shorter, with np.newaxis

np.newaxis is None; choose the spelling that exposes the axis role. keepdims=True performs a related shape-preserving operation after a reduction. Neither mechanism knows that your dimension is “channel,” “token” or “batch.”

5 · How broadcasting makes neural nets legal

X:(32,784) times W:(784,128) produces(32,128). A bias(128,) pads to(1,128) and is shared across observations. A per-example offset would instead need(32,1); one bias per training row has different generalization semantics.

The backward rule follows from reuse: if Yij=Zij+bj and Gij=dL/dYij, then $db_j=\sum_iG_{ij}$. Sum along every replicated leading axis, not arbitrarily average. A loss mean's denominator must appear exactly once.

6 · Broadcasting avoids operand replication; outputs still allocate

np.broadcast_to(row,(1000000,4)) returns a read-only view whose leading stride is0. For an int64 row it reports32,000,000 logical nbytes, while sharing the32-byte row storage. nbytes is not lying: it counts logical element bytes, not unique allocation. Read-only applies to that view; another writable alias can still change the backing row.

NumPy's broadcasting iterator need not literally create a broadcast_to object for every expression. The documented promise is avoiding needless operand copies, not that every operation is free or every broadcasted object in every framework is immutable. In-place operations also cannot generally expand the left-hand output shape. Use out= only with correct shape, casting and aliasing contracts.

row=np.arange(4,dtype=np.int64);wide=np.broadcast_to(row,(1_000_000,4))
assert wide.strides==(0,8) and np.shares_memory(row,wide)
assert wide.nbytes==32_000_000 and row.nbytes==32 and not wide.flags.writeable
try: wide[0,0]=1
except ValueError: pass
else: raise AssertionError('read-only view accepted mutation')
row[0]=7; assert wide[999999,0]==7
result=np.zeros((3,4))+row
assert result.nbytes==3*4*8 and not np.shares_memory(result,row)
a=np.ones((1,3))
try: a+=np.ones((2,3))
except ValueError: pass
else: raise AssertionError('in-place output incorrectly expanded')
print('zero-stride storage, logical bytes, alias mutation and in-place limits passed')

7 · Three silent-shape failure cases

Accuracy: equal-shaped prediction/label vectors compare matching examples. Comparing(N,) with(N,1) computes all pairs, whose average agreement is $\sum_c f_{pred}(c)f_{label}(c)$, determined by class frequencies. It need not be10%. A shuffled balanced prediction list can leave this wrong metric unchanged. For one scalar label per example, reshape deliberately to(N,) or squeeze the known label axis, not every singleton dimension.

Squared error: an(N,N) residual matrix compares each prediction with every target. Its mean can decrease while learning the wrong objective. Assert the residual shape before reducing and test N>1.

Channel normalization: for CHW image(3,224,224), per-channel means(3,) align against width and fail. Use(3,1,1) or a keepdims reduction. If all dimensions happen to be3, the wrong operation can silently succeed; compare known channel means numerically. Dividing by a zero mean is a separate domain error; standardization usually subtracts means and scales by a declared nonzero standard deviation policy.

8 · The one weird rule for matmul

For operands of rank at least2, matmul contracts the last two matrix axes and broadcasts leading batch axes. Thus(32,5,3)@(3,7)→(32,5,7), and(2,1,5,3)@(1,4,3,7)→(2,4,5,7). Inner matrix dimensions must match independently of batch compatibility.

One-dimensional operands have special promotion/removal rules: vectors of length d dotted together return a scalar; a matrix times a d-vector returns a vector; scalars are not matmul operands. NumPy dot has different high-rank rules. For Q,K:(B,H,T,d), use K.swapaxes(-2,-1), not K.T, to compute per-head scores.

8b · Five transformer broadcasting patterns

8b.1 · Adding positional bias to attention scores

For scores(B,H,T,T), a bias(H,T,T) is shared over batch. ALiBi uses a head-specific slope times a relative-position penalty added to attention logits, together with the relevant causal mask. For causal self-attention at query i and key j≤i, a common representation is $-m_h(i-j)$; different slopes produce different per-head recency biases. A(T,T) bias would tie that bias across heads. Extrapolation in the paper is empirical and bounded, not unlimited-context capability.

8b.2 · Causal masking

A lower-triangular Boolean(T,T) mask can be reused across B,H. In NumPy, np.where(mask,scores,-np.inf) returns a score-shaped result; ndarray has no masked_fill. A dense Boolean mask costs T² bytes at one byte per element: atT=4096 it is16MiB, not a fixed4KB. A(B,H,T,T) score costs BHT² times its element size. A fused kernel can avoid explicit buffers, but a simple broadcast does not do that automatically. Cached rectangular decoding requires offset-aware causality, not blindly applying a square triangle.

8b.3 · RMSNorm scale

For x:(B,T,d), rms=√(mean(x²,last,keepdims)+ε) has(B,T,1); x/rms then gamma:(d,) preserves(B,T,d). Unlike LayerNorm, RMSNorm does not subtract the mean. Epsilon belongs under the square root in this defined formula. The primary paper motivates rescaling invariance; epsilon and learned scale qualify exact invariance. This is not evidence every model uses identical normalization.

8b.4 · RoPE rotation

Pair adjacent coordinates in an even-dimensional head. For angle θ at position t, map(a,b) to(a cosθ−b sinθ,a sinθ+b cosθ). A(T,d/2) table broadcasts over batch and heads. A “rotate_half” implementation using split halves requires a matching table/layout; mixing adjacent-pair and half-split conventions produces the wrong rotation despite valid shapes. Norm preservation and a2×2 numeric oracle catch it. LongRoPE/YaRN are context-extension methods, not just inserting singleton axes or guaranteed explanations of named10M-context releases.

8b.5 · Grouped-Query Attention repeat_interleave

With32 query heads and8 KV heads, ordinary head broadcasting fails:32 and8 are unequal and neither is1. A simple reference repeats each KV head four times. A grouped formulation reshapes Q into(B,8,4,T,d) and uses K:(B,8,1,T,d), reusing each KV head across its query group. NumPy repeat/PyTorch repeat_interleave materialize; specialized kernels may avoid that repetition. For equal per-head dimensions and dtype, KV storage is4 times smaller than32-head MHA, not automatically4 times less whole-model memory or latency. Quality depends on training/adaptation, not an assumption-free “tiny loss.”

rng=np.random.default_rng(3)
B,H,T,D=2,4,3,4
q,k,v=[rng.normal(size=(B,H,T,D)) for _ in range(3)]
raw=q@k.swapaxes(-2,-1)/np.sqrt(D)
pos=np.arange(T);distance=pos[:,None]-pos[None,:]
slopes=np.array([.5,.25,.125,.0625]);bias=-slopes[:,None,None]*distance[None,:,:]
mask=distance>=0
scores=np.where(mask,raw+bias,-np.inf)
weights=np.exp(scores-scores.max(-1,keepdims=True));weights/=weights.sum(-1,keepdims=True)
out=weights@v
assert out.shape==(B,H,T,D) and np.all(weights[...,~mask]==0)
np.testing.assert_allclose(out[:,:,0],v[:,:,0])
assert not np.allclose(bias[0],bias[1])
x=rng.normal(size=(B,T,D));gamma=np.arange(1,D+1);eps=1e-6
rms=np.sqrt((x*x).mean(-1,keepdims=True)+eps);normalized=x/rms*gamma
assert normalized.shape==x.shape
np.testing.assert_allclose(normalized[0,0],x[0,0]/np.sqrt(sum(x[0,0]**2)/D+eps)*gamma)
angles=pos[:,None]*np.array([1.,.1])[None,:]
c=np.cos(angles)[None,None];s=np.sin(angles)[None,None]
a=q[...,0::2];b=q[...,1::2]
rot=np.empty_like(q);rot[...,0::2]=a*c-b*s;rot[...,1::2]=a*s+b*c
np.testing.assert_allclose((rot*rot).sum(-1),(q*q).sum(-1),atol=1e-12)
manual=np.array([[np.cos(1),-np.sin(1)],[np.sin(1),np.cos(1)]])@q[0,0,1,:2]
np.testing.assert_allclose(rot[0,0,1,:2],manual)
# A small GQA grouping: 4 query heads, 2 KV heads, 2 queries per group.
k2=rng.normal(size=(B,2,T,D));qg=q.reshape(B,2,2,T,D)
grouped=qg@k2[:,:,None].swapaxes(-2,-1)
repeated=q@np.repeat(k2,2,axis=1).swapaxes(-2,-1)
np.testing.assert_allclose(grouped.reshape(B,H,T,T),repeated)
assert not np.shares_memory(k2,np.repeat(k2,2,axis=1))
assert 4096**2==16*1024**2
print('ALiBi, causal mask, RMSNorm, paired rotary and grouped-query oracles passed')

9 · The full mental model — a decision tree

For binary elementwise operations, align and check operand dimensions. Unary exp/log simply act on the existing shape; log also needs a valid numerical domain. For matmul, handle vector promotions, check contraction dimensions, then broadcast batch dimensions. For reductions, identify the axes combined and whether keepdims retains singleton positions. Norm defaults can differ from simple sum defaults, so read the actual API.

The conceptual repeated rows say which values are reused, not how many independent bytes are allocated. For(5,)+(3,5), each of the three rows receives the same five offsets. Multiplying the output dimensions by itemsize estimates logical result bytes; it does not count every temporary or retained parent.

Pairwise-distance tradeoff: broadcast, algebra or blocks

Expanding the square gives $\|x-c\|^2=\|x\|^2+\|c\|^2-2x^Tc$. This replaces an(N,K,d) differences tensor by an(N,K) score calculation, but cancellation can destroy tiny distances between large nearby vectors. Clamping small negative outputs to0 prevents a square-root NaN; it cannot recover lost significant digits or fix arbitrary errors. Centering, wider precision, direct differences or blockwise evaluation can be preferable.

A chunked direct calculation caps the temporary at(chunk,K,d), but storing the full output still costs NK. If only nearest neighbors are needed, maintain a running top-k instead. L1 and infinity distances do not use this same squared-L2 identity. Benchmark equal computations, not a full cosine table against one nearest neighbor.

import timeit
rng=np.random.default_rng(30);A=rng.normal(size=(17,5));C=rng.normal(size=(11,5))
def direct(A,C): return ((A[:,None,:]-C[None,:,:])**2).sum(-1)
def expanded(A,C): return (A*A).sum(1,keepdims=True)+(C*C).sum(1)[None,:]-2*A@C.T
def blocked(A,C,block=4):
    result=np.empty((len(A),len(C)))
    for i in range(0,len(A),block):result[i:i+block]=direct(A[i:i+block],C)
    return result
oracle=np.array([[sum((float(x)-float(y))**2 for x,y in zip(a,c)) for c in C] for a in A])
np.testing.assert_allclose(direct(A,C),oracle,atol=1e-12)
np.testing.assert_allclose(expanded(A,C),oracle,atol=1e-12)
np.testing.assert_allclose(blocked(A,C),oracle,atol=1e-12)
large=np.array([[1e8,1e8]]);near=np.array([[1e8+1,1e8]])
assert direct(large,near)[0,0]==1
assert expanded(large,near)[0,0]!=1  # cancellation: clamp cannot recover it
pred=np.array([0,1,0,1]);labels=np.array([0,1,0,1])[:,None]
assert (pred==labels).shape==(4,4) and (pred==labels).mean()==.5
assert (pred==labels[:,0]).mean()==1.
img=np.arange(27,dtype=float).reshape(3,3,3)+1
means=img.mean((1,2),keepdims=True)
assert not np.allclose(img/means,img/means[:,0,0])
print('direct/expanded/blocked distances, cancellation and silent shape controls passed')

11 · Retention scaffold

1. Three rules? Right-align, require equal or singleton sizes, then select the non-singleton/equal size. Remember(0,)+(1,)→(0,), so maximum is not the general rule.
2. (5,1)+(1,4)? The first operand reuses each row value across four columns; the second reuses its row over five observations, producing(5,4).
3. Prediction/label trap? (100,) against(100,1) gives(100,100). Assert one scalar target per prediction and select the known target axis; a blind squeeze can remove batch dimensions when N=1.
4. Which pairwise implementation is faster? Bulk kernels often remove Python overhead, but a1000×1000×128 float64 temporary alone costs1.024GB. Measure under a feasible memory budget and consider expanded/blocked algorithms; no100–500× promise follows.
5. When does it raise? Any aligned unequal pair with neither size1 is incompatible. Also an in-place output may be unable to expand even when an out-of-place expression is valid.
6. **Why repeat GQA heads?**32 cannot broadcast against8. Four query heads share each KV head; a grouped/fused implementation can express that sharing without a repeated KV buffer, but quality and hardware performance need separate evidence.

Recall — from memory

Explain the difference between a zero-stride operand view and a newly allocated result, and give(32,5,3)@(3,7)→(32,5,7). Broadcasting does not promise the same speed as a same-shaped contiguous operation: traversal, output size and kernels matter.

Stretch — for one extra hour

Normalize rows of A and B with retained singleton norm axes, reject zero directions, then multiply unit A by unit B.T. Validate a rectangular case, zero row, a noncontiguous view and unit diagonal for a self-comparison. The following bounded timing compares exactly the same cosine table. Scaling to1000×1000 is optional; the code does not allocate the1GB distance temporary.

def pairwise_cosine_similarity(A,B):
    A=np.asarray(A,float);B=np.asarray(B,float)
    assert A.ndim==B.ndim==2 and A.shape[1]==B.shape[1]
    an=np.linalg.norm(A,axis=1,keepdims=True);bn=np.linalg.norm(B,axis=1,keepdims=True)
    if np.any(an==0) or np.any(bn==0):raise ValueError('undefined zero direction')
    return (A/an)@(B/bn).T
A=rng.normal(size=(60,16))[:,::2];B=rng.normal(size=(50,8))
def bulk():return pairwise_cosine_similarity(A,B)
def loop():return np.array([[np.dot(a,b)/(np.linalg.norm(a)*np.linalg.norm(b)) for b in B] for a in A])
np.testing.assert_allclose(bulk(),loop(),atol=1e-12)
np.testing.assert_allclose(np.diag(pairwise_cosine_similarity(A,A)),1,atol=1e-12)
try:pairwise_cosine_similarity(np.zeros((1,8)),B)
except ValueError:pass
else:raise AssertionError('zero vector accepted')
bulk();loop()
print('local seconds bulk/loop',min(timeit.repeat(bulk,number=1,repeat=3)),min(timeit.repeat(loop,number=1,repeat=3)))

In your own words

Say which values are reused and which output entries are newly computed; include the result allocation in your explanation.

Spaced-review pointer

Tomorrow predict all seven shapes; next week reconstruct the GQA grouping and the cancellation counterexample. Re-reading familiar punctuation is not the same as predicting its semantics.

Next-session teaser

DL004 reverses these operations: shared forward operands receive summed backward contributions, with matching parameter shapes.

What to bring back tomorrow — sticky note

Align right, insert axes deliberately, count output bytes, and test unequal dimensions. A batch of one hides many accidental all-pairs bugs.

Primary references and boundary

NumPy2.2 broadcasting, broadcast_to, matmul, ALiBi, RMSNorm, RoFormer and GQA support the contracts. Tests execute these small NumPy mechanisms, not full pretrained architectures or fused GPU kernels. Timing is local and shape-specific.