Skip to content
Vectorization Patterns
step 1/5

Reading — step 1 of 5

Learn

~2 min readVectorization, ODEs, Optimization

Octave (like MATLAB) is much faster with vectorized operations than with explicit loops. Learn the patterns.

Loop vs vectorized

Slow:

N = 1000000;
result = zeros(1, N);
for i = 1:N
    result(i) = i^2;
end

Fast (10-100x):

result = (1:N).^2;

The vectorized form dispatches to optimized BLAS routines.

Logical indexing

Filter without writing a loop:

v = [3, 1, 4, 1, 5, 9, 2, 6];
v(v > 3)            %% [4 5 9 6] — all elements > 3
v(mod(v, 2) == 0)   %% [4 2 6] — evens
v(v >= 3 & v <= 6)  %% [3 4 5 6] — in range

v > 3 produces a logical vector [0 0 1 0 1 1 0 1]. Indexing with logical vectors selects matching positions.

Conditional update

v = -5:5;
v(v < 0) = 0;       %% replace negatives with 0

find for indices

v = [3, 1, 4, 1, 5, 9, 2, 6];
find(v > 3)         %% [3 5 6 8] — positions of matching elements
find(v == 1)         %% [2 4]

Useful when you want positions, not values.

Broadcasting

Operations between scalars and arrays auto-extend:

v = [1 2 3];
v + 10              %% [11 12 13]
v * 2               %% [2 4 6]

Between arrays of compatible sizes:

A = [1 2 3; 4 5 6];     %% 2x3
v = [10 20 30];          %% 1x3
A + v                    %% [11 22 33; 14 25 36] — v added to each row

bsxfun (broadcast function)

For older Octave versions or non-trivial broadcasts:

bsxfun(@plus, A, v)      %% same as A + v
bsxfun(@times, A, v)

meshgrid for grid coordinates

[X, Y] = meshgrid(1:5, 1:3);
%% X = [1 2 3 4 5; 1 2 3 4 5; 1 2 3 4 5]
%% Y = [1 1 1 1 1; 2 2 2 2 2; 3 3 3 3 3]
Z = X.^2 + Y.^2;             %% function on the grid

When loops ARE OK

  • Genuinely sequential algorithms (Newton's method, simulations)
  • Operations that can't vectorize (file I/O per row, complex branching)
  • Small N where loop overhead is irrelevant

Rule of thumb: if you find yourself writing a for loop, ask "can this be a vectorized expression?". 90% of the time the answer is yes.

Performance test

tic
result = zeros(1, 1000000);
for i = 1:1000000
    result(i) = i^2;
end
loop_time = toc;

tic
result = (1:1000000).^2;
vec_time = toc;

printf("loop: %.3f sec\nvec: %.3f sec\nspeedup: %.1fx\n", 
        loop_time, vec_time, loop_time / vec_time);

Typical result: 50-200x speedup.

Discussion

Ask a question, share an insight, or help someone who’s stuck.

Sign in to post a comment or reply.

Loading…