# Vectorizing nested for loops

3 views (last 30 days)
Commented: Adam on 3 Apr 2019
How do I vectorize these for loops?
I am looking for the fastest implementation.
Thank you!
My current implementation is
d = zeros(N,N);
for i = 1:N
for j = i:N
d(i,j) = vecnorm( Phi(:,i) - Phi(:,j) );
end
end
d = d + d';

Jos (10584) on 3 Apr 2019
A = magic(5)
d = squareform(pdist(transpose(A))) % transpose to obtain vecnorm between columns
pdist and squareform are part of the Statistics Toolbox
Excellent! This is about an order of magnitude faster than the other answer.

Bjorn Gustavsson on 3 Apr 2019
You find a nice contribution from J Kirk on the file exchange: Distmat
In my historic appreciation efforts there is this Acklamized function for distmat:
function d = distmat1(x, y)
% DISTMAT1 Distance matrix (one for-loop).
%
% D = DISTMAT1(X, Y) returns the distance matrix with all distances
% between the points represented by the rows of X and Y.
%
% DISTMAT1(X) is equivalent to DISTMAT1(X, X), but the former computes the
% distance matrix faster.
%
% Distance is Euclidean.
%
% The calculation is done with one for-loop.
%
% Author: Peter J. A.
% Time-stamp: 2000-07-10 20:20:27
% E-mail: yyyyyyy@xxxx.xx
% URL: http:xxxxxxxxxxxxxxx
error(nargchk(1, 2, nargin));
if nargin == 1 % DISTMAT1(X)
if ndims(x) ~= 2
error('Input must be a matrix.');
end
m = size(x, 1);
d = zeros(m, m); % initialise output matrix
for i = 1:m-1
xi = x(i,:);
diff = xi(ones(m-i,1),:) - x(i+1:m,:); % difference
dist = sqrt(sum(abs(diff).^2, 2)); % distance
d(i+1:m,i) = dist;
d(i,i+1:m) = dist.';
end
else % DISTMAT1(X, Y)
if ndims(x) ~= 2 | ndims(y) ~= 2
error('Input must be two matrices.');
end
[mx, nx] = size(x);
[my, ny] = size(y);
if nx ~= ny
error('Both matrices must have the same number of columns.');
end
m = mx; % number of rows in distance matrix
n = my; % number of columns in distance matrix
p = nx; % dimension of each point
d = zeros(m, n); % initialise output matrix
% The for-loop is applied on the shortest dimension, the other is
% vectorized.
if m < n
idx = ones(1, n);
for i = 1:m
xi = x(i,:);
d(i,:) = sqrt(sum(abs(y - xi(idx,:)).^2, 2)).';
end
else
idx = ones(1, m);
for j = 1:n
yj = y(j,:);
d(:,j) = sqrt(sum(abs(x - yj(idx,:)).^2, 2));
end
end
end
Maybe no longer optimal coding - but the originator of this code was a master at writing good vectorized terse code, back in the olden days before bxfun arrayfun and all the tools making the youth of today .....
HTH

R2018b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!