## Vectorizing nested for loops

on 3 Apr 2019
Latest activity Commented on by Adam

on 3 Apr 2019

### Jos (10584) (view profile)

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';

R2018b

### Jos (10584) (view profile)

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

#### 1 Comment

on 3 Apr 2019
Excellent! This is about an order of magnitude faster than the other answer.

### Bjorn Gustavsson (view profile)

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