Implicit expansion
Implicit expansion — the feature that NumPy and others call broadcasting — is the rule by which binary operations on arrays of compatible but not identical shape are extended to arrays of the same shape. The rule has been in MATLAB since R2016b; before that, the function bsxfun (binary singleton expansion function) was required to obtain the same effect. The change made bsxfun largely obsolete and made one of the language’s most-elaborate idioms invisible.
The rule, briefly: two arrays have compatible shape if, for each dimension, either the sizes are equal or one of them is 1. The result has the size of the larger along each dimension. A dimension of size 1 in one operand is implicitly repeated to match the size of the other.
A first example
Subtract the mean of each column from the columns of a matrix:
A = magic(4);
mu = mean(A) % 1×4 row vector
A - mu
% A is 4×4; mu is 1×4. The shapes are compatible:
% dim 1: 4 vs 1 → 4
% dim 2: 4 vs 4 → 4
% The result is 4×4 with mu subtracted from each row.The pre-R2016b idiom — still occasionally seen in legacy code — was:
A - bsxfun(@minus, A, mu) % the old way
A - mu % the new way; identical resultThe rule
Given two arrays A and B with sizes [a1 a2 a3 ...] and [b1 b2 b3 ...], the operation A ⊕ B (for any element-wise operator ⊕) is defined when, for every dimension k, either a_k == b_k, a_k == 1, or b_k == 1. The result has size max(a_k, b_k) along each dimension. The operand with size 1 along a dimension is broadcast (repeated) to the size of the other.
Shapes that do not satisfy the rule raise the error “Arrays have incompatible sizes for this operation”.
The rule applies to every element-wise operator: arithmetic (+, -, .*, ./, .^), comparison (<, <=, >, >=, ==, ~=), and logical (&, |, xor). It applies to most standard-library functions when invoked with two arguments (max(A, B), min(A, B), hypot(A, B), atan2(Y, X)). It applies to user-defined functions via the wrapper bsxfun(@fun, A, B) and, since R2018a, via pagefun for per-page broadcasting on 3-d arrays.
Three patterns
Column-wise normalisation. A 1×N row vector applied to an M×N matrix broadcasts down the rows:
A = rand(100, 5);
A_centred = A - mean(A); % subtract per-column mean
A_normalised = A_centred ./ std(A); % divide by per-column stdRow-wise normalisation. A column vector applied to a matrix broadcasts across the columns:
A = rand(5, 100);
row_means = mean(A, 2); % 5×1
A_centred = A - row_means; % broadcast over columnsOuter product. A column vector and a row vector combine into an outer product:
u = (1:5).'; % 5×1
v = 1:3; % 1×3
u * v % matrix product, 5×3 (the linear-algebra outer product)
u + v % broadcast, 5×3 (an "additive outer product")
u .* v % broadcast, 5×3 (same as u * v in this case)The matrix product u * v of a column by a row is equivalent to the broadcasting u .* v because matrix multiplication of those shapes happens to coincide with element-wise broadcasting — a small coincidence that experienced MATLAB programmers internalise.
Tabulating a function over a grid
A canonical use of broadcasting is producing a 2-d array of values of f(x, y) over a grid:
x = linspace(-2, 2, 200); % 1×200
y = linspace(-2, 2, 200).'; % 200×1
Z = exp(-(x.^2 + y.^2)); % broadcast: 200×200 surface
surf(x, y, Z)The functions meshgrid and ndgrid — venerable companions of MATLAB — produce 2-d arrays of x and y coordinates explicitly; with implicit expansion they are no longer necessary for the simple case, and the column-vs-row form is the clean idiom in modern code.
Higher-dimensional broadcasting
The rule extends to N-d. A 1×1×K array combines with an M×N matrix to produce an M×N×K stack:
A = rand(4, 4); % 4×4
w = reshape([1 2 3], 1, 1, 3); % 1×1×3
A .* w % 4×4×3The construct is the canonical way to apply K different scale factors to a single matrix without an explicit loop. It is the basis of much batched linear algebra (pagemtimes, pagemldivide) and of vectorised numerical methods.
Singletons and the role of permute
To broadcast in a way the rule does not directly support, the operand must be reshaped so that its singleton dimensions line up. permute and reshape are the tools:
A = rand(3, 4); % 3×4
B = rand(4, 5); % 4×5
% We want an outer product over their second dimensions: 3×4 × 4×5 (via broadcasting,
% NOT matrix product). Reshape A to 3×4×1 and B to 1×4×5, then broadcast.
C = A .* reshape(B, 1, 4, 5); % 3×4×5The mental model: each operand is padded with singleton dimensions on the left (or on the right, by reshape) until they have the same number of dimensions, and the rule applies dimension-by-dimension.
Performance
Implicit expansion does not materially copy the array: the implementation walks the smaller operand with a stride of 0 along the broadcast dimension. The performance is the same as the explicit bsxfun form and the same as a fused for-loop in compiled code. The principal performance advice — vectorise, do not loop — applies unchanged.
Pitfalls
The most-cited pitfall is the unintended broadcast: a row vector subtracted from a column vector silently produces an outer-product-shaped array of differences. The reader expecting a scalar or a vector receives a matrix instead and is confused by downstream errors.
u = [1; 2; 3]; % 3×1
v = [10 20 30]; % 1×3
u - v
% 3×3:
% -9 -19 -29
% -8 -18 -28
% -7 -17 -27The remedy is to be deliberate about row/column orientation; the transpose operator .' is cheap and explicit. In R2017a and later, the function-argument-validation block with (1, :) and (:, 1) size requirements is the recommended way to enforce orientation at function boundaries.
A second pitfall is interaction with the matrix product: A * B is the linear-algebra product and does not broadcast; the operator requires size(A, 2) == size(B, 1). The element-wise product A .* B is the broadcasting one. Confusing them is the source of a recurring class of bugs.