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

Answer by Jos (10584)
on 3 Apr 2019

Accepted Answer

A = magic(5)

d = squareform(pdist(transpose(A))) % transpose to obtain vecnorm between columns

pdist and squareform are part of the Statistics Toolbox

Answer by 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.

%

% See also DISTMAT0, DISTMAT2, DISTMAT3.

% 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

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.