03 · Vectorization vs Loops¶
Verification note
MATLAB was not available in the environment used to write this page. Every code sample below was hand-traced against MATLAB's documented array-operator and broadcasting ("implicit expansion") semantics, and cross-checked against equivalent NumPy behavior where applicable. It was not executed in MATLAB itself.
MATLAB's entire design center is the matrix. That means the fastest, most idiomatic MATLAB code almost never contains an explicit loop over array elements — it expresses the whole computation as one array operation. This module is about recognizing loop patterns that should become vectorized expressions, and knowing when a loop is still the right tool.
Why vectorization matters¶
A for loop over 1:N in MATLAB pays per-iteration overhead: bounds
checks, interpreter dispatch, and (historically) JIT-compilation limits.
A vectorized expression instead calls into MATLAB's underlying compiled
array library once for the whole array. The difference is not cosmetic —
on large arrays it is routinely a 10-100x speedup, and it also tends to
produce shorter, more readable code once you're fluent in the idiom.
n = 1e6;
x = 1:n;
% Loop version
tic;
y_loop = zeros(1, n);
for i = 1:n
y_loop(i) = x(i)^2 + 3*x(i) - 1;
end
t_loop = toc;
% Vectorized version
tic;
y_vec = x.^2 + 3*x - 1;
t_vec = toc;
Both produce identical results element-for-element, but t_vec is
typically one to two orders of magnitude smaller than t_loop. Modern
MATLAB's JIT accelerator narrows this gap for simple loops, but the gap
never fully closes, and vectorized code stays easier to read and to
compose with the rest of the array-oriented toolbox (plotting, reductions,
logical indexing).
The core idiom: element-wise operators¶
Recall from Level 1 that .^, .*, ./ operate element-by-element,
while ^, *, / are matrix operations. Vectorization leans on the
element-wise family constantly:
A loop written as:
is mechanically translated to the vectorized form by replacing x(i) with
x and every *///^ that acts elementwise with its dotted form.
This translation works whenever loop iteration i only reads/writes
index i — no dependency on i-1 or accumulation across iterations.
Logical indexing replaces if inside loops¶
A very common loop pattern is "collect elements meeting a condition":
data = [-3, 7, -1, 9, 0, -5, 4];
positives_loop = [];
for i = 1:length(data)
if data(i) > 0
positives_loop(end+1) = data(i);
end
end
The vectorized equivalent uses a logical mask:
data > 0 produces a logical array the same size as data. Indexing
data with a logical array of matching size selects exactly the true
positions, in order — no manual bookkeeping of an output index.
Assignment works the same way, and is where logical indexing really shines over loops:
The loop equivalent is four lines with an if; the vectorized form is
one line and reads as "wherever data is negative, set it to zero."
find when you need indices, not values¶
Sometimes you need positions, not values — find returns the linear
indices where a condition holds:
idx = find(data > 0); % [2 4 7]
first_idx = find(data > 0, 1); % 2 (first match only)
last_idx = find(data > 0, 1, 'last'); % 7
find(cond, 1) short-circuits after the first match and is the
vectorized replacement for a loop with a break once a condition is hit.
Accumulation patterns: cumsum, cumprod, diff¶
Loops that carry a running total across iterations map onto built-in cumulative functions:
values = [10 20 5 15 30];
running_total_loop = zeros(size(values));
total = 0;
for i = 1:length(values)
total = total + values(i);
running_total_loop(i) = total;
end
% running_total_loop = [10 30 35 50 80]
running_total = cumsum(values); % identical result, one call
diff(values) computes successive differences (values(2:end) -
values(1:end-1)), the vectorized inverse of a loop that computes
period-over-period change:
Vectorizing 2D operations: meshgrid / broadcasting¶
A common beginner loop evaluates a function of two variables over a grid:
xv = -2:0.5:2;
yv = -2:0.5:2;
Z_loop = zeros(length(yv), length(xv));
for i = 1:length(yv)
for j = 1:length(xv)
Z_loop(i, j) = xv(j)^2 + yv(i)^2;
end
end
meshgrid builds coordinate matrices once, then a single vectorized
expression evaluates the whole grid:
Modern MATLAB (R2016b+) also supports implicit expansion
(broadcasting): a row vector and a column vector combine directly without
meshgrid at all —
x = -2:0.5:2; % 1x9 row
y = (-2:0.5:2)'; % 9x1 column
Z2 = x.^2 + y.^2; % 9x9, broadcast automatically
Z2 is numerically identical to Z above; x.^2 broadcasts across rows
and y.^2 broadcasts across columns, and + combines them into a full
9x9 matrix. This is the single biggest vectorization upgrade over older
MATLAB code you'll see in legacy scripts still using meshgrid /
repmat out of habit.
arrayfun and cellfun: vectorizing calls to non-vectorizable functions¶
Sometimes the per-element operation is a function that doesn't support
array inputs natively — often a function with internal branching. arrayfun
applies a function to every element of a numeric array and collects the
results; cellfun does the same for cell arrays:
classify = @(v) merge(v > 0, 'pos', merge(v < 0, 'neg', 'zero'));
% (illustrative only: MATLAB has no `merge`; a real classify would use if/else)
function s = classify_sign(v)
if v > 0
s = "pos";
elseif v < 0
s = "neg";
else
s = "zero";
end
end
labels = arrayfun(@classify_sign, data); % applies elementwise, output as an array
Because classify_sign returns a string scalar (non-uniform across
calls in general), you'd pass 'UniformOutput', false if outputs vary
in type or size:
words = arrayfun(@(v) sprintf('%d units', v), data, 'UniformOutput', false);
% words is a cell array of char row vectors, since sprintf output length varies
arrayfun/cellfun are not automatically faster than an equivalent
for loop — they still call the function once per element under the
hood — but they read cleanly as a vectorized statement of intent, and
they compose well with cellfun(@isempty, c)-style existing-function
checks over cell arrays.
When to keep the loop¶
Vectorization is not a universal replacement. Keep an explicit loop when:
- The iteration is inherently sequential — each output depends on the
previous output, such as a discrete-time simulation
x(k+1) = f(x(k))or a recursive filter. There's no way to computex(k+1)beforex(k)exists, so no array expression replaces it. - The vectorized form would need to materialize a huge intermediate
array — e.g. an
NxNdistance matrix whenNis in the millions — and a loop (possibly chunked) keeps memory bounded. - Readability actually suffers. A three-level broadcasting expression that takes a paragraph to explain is worse code than a four-line loop with a comment, even if it benchmarks faster.
- The body has genuine side effects per iteration — writing files,
updating a GUI, calling external hardware — where "at position
i, do X" is the actual semantics, not a data transformation.
A recursive/sequential example that must stay a loop:
n = 20;
fib = zeros(1, n);
fib(1) = 1; fib(2) = 1;
for i = 3:n
fib(i) = fib(i-1) + fib(i-2); % depends on prior loop outputs
end
There is no vectorized Fibonacci here (short of a closed-form Binet's-formula approximation, which trades exactness for speed and isn't the same computation).
A worked before/after¶
Loop version — compute a weighted moving average of window 3:
x = [4 8 6 5 9 3 7];
n = length(x);
w = [0.2 0.6 0.2];
y_loop = zeros(1, n-2);
for i = 1:n-2
window = x(i:i+2);
y_loop(i) = sum(window .* w);
end
Vectorized version using conv:
conv flips one input internally (convolution vs. correlation), so
fliplr(w) cancels that flip and reproduces the same weighted-sum
semantics as the loop. 'valid' restricts output to positions where the
window fully overlaps x, matching the loop's n-2 output length.
Both give y_loop = y = [0.2*4+0.6*8+0.2*6, ...] = [6.8, 6.6, 6.6, 6.8, 5.6].
How It Actually Works¶
The vectorization speedup isn't folklore — it comes from three concrete, stackable mechanisms:
- Interpreter dispatch overhead. Every statement inside an explicit
forloop pays a fixed per-iteration cost: the interpreter re-enters its execution loop, looks up variable names in the workspace hash table, and re-dispatches the operator based on the runtime types involved. A vectorized call likey = sin(x)pays that dispatch cost once, regardless of whetherxhas 10 or 10 million elements. - JIT eligibility. Since MATLAB's JIT compiler was introduced, some
loops do get compiled to native code — but only when the loop body is
simple enough (no calls to non-JIT-eligible functions, no arrays
changing size or class inside the loop, no
try/catch). A loop that calls a user function or builds up an array withend+1typically falls out of JIT eligibility and pays full interpretation cost every iteration, while an equivalent vectorized expression routes straight into compiled, cache-tuned array-library code regardless of JIT status. - Memory access pattern. A vectorized operation over a whole array
walks memory in one predictable, contiguous (column-major) sweep,
which is friendly to CPU cache prefetching; a loop that grows an array
with
A(end+1) = xforces repeated full reallocation and copy of the growing buffer — genuinelyO(n^2)total work fornappends, versusO(n)for a preallocated vectorized fill.
None of this means loops are "wrong" — a loop with genuine data-dependence between iterations (each step needs the previous step's result, as in many ODE solvers) often cannot be vectorized at all, and a JIT-friendly loop over a preallocated array can be competitive with vectorized code for such cases.
Note: reasoned from MATLAB's documented interpreter/JIT execution model and general cache-locality principles; not benchmarked in a live MATLAB session (no MATLAB installation is available here).
Summary¶
| Loop pattern | Vectorized replacement |
|---|---|
for i, y(i) = f(x(i)) elementwise math |
.^, .*, ./ on whole arrays |
if inside loop collecting matches |
logical indexing: x(x > 0) |
| loop tracking a running total | cumsum, cumprod |
| loop computing successive differences | diff |
| nested loop over a 2D grid | meshgrid + .^/.*, or implicit expansion |
| loop calling a non-vectorized function per element | arrayfun / cellfun |
| loop where output(k) depends on output(k-1) | keep the loop |
The practical workflow: write the loop first if that's the clearest way
to think through the logic, get it correct on a small test case, then ask
"does iteration i only touch index i?" — if yes, vectorize; if the
iteration threads state forward, leave it as a loop.