A faster and more compact way to create a list of distances among all the pairs of points

Hi, could you suggest a faster and more compact way to create a list of distances among all the pairs of points?

My attempt here below:

% Input (x and y coordinates of 6 points)
x = [1 2 2 3 4 5];
y = [1 2 3 7 2 5];
% Plot just to see the 6 points
plot(x,y,'o','MarkerFaceColor','b','markersize',15)
xlim([0 10])
ylim([0 10])

% Calculate the distances among each pair of points
Z = squareform(pdist([x' y']));
% Create a list that includes 3 elements: i-point ID, j-point ID, distance(i,j)
k = 1;
for i = 1 : length(x)-1
    for j = i+1 : length(x)
        list(k,:) = [i j Z(i,j)];
        k = k + 1;
    end
end
list,

list = 15×3

1.0000 2.0000 1.4142 1.0000 3.0000 2.2361 1.0000 4.0000 6.3246 1.0000 5.0000 3.1623 1.0000 6.0000 5.6569 2.0000 3.0000 1.0000 2.0000 4.0000 5.0990 2.0000 5.0000 2.0000 2.0000 6.0000 4.2426 3.0000 4.0000 4.1231

Answer 1

[chicken vector]

N = 1e4;
x = randi(10,N,1);
y = randi(10,N,1);
tic
xIdx = repmat(1 : length(x), length(x), 1);
yIdx = xIdx';
vectorIdx = (1 : size(xIdx, 1))' > (1 : size(xIdx, 2));
xy = [x(:), y(:)];
dist = pdist2(xy, xy);
distPdist = dist(vectorIdx);
list = [xIdx(vectorIdx) , ...
    yIdx(vectorIdx) , ...
    distPdist]

list = 49995000×3

1.0000 2.0000 5.3852 1.0000 3.0000 6.7082 1.0000 4.0000 3.0000 1.0000 5.0000 5.8310 1.0000 6.0000 8.6023 1.0000 7.0000 2.2361 1.0000 8.0000 8.5440 1.0000 9.0000 8.5440 1.0000 10.0000 5.0000 1.0000 11.0000 10.8167

toc

Elapsed time is 2.501384 seconds.

Answer 2

[Image Analyst]

Try pdist2

% Input (x and y coordinates of 6 points)
x = [1 2 2 3 4 5];
y = [1 2 3 7 2 5];
xy = [x(:), y(:)]

xy = 6×2

1 1 2 2 2 3 3 7 4 2 5 5

% Get distances between every (x,y) point and every other (x,y) point:
distances = pdist2(xy, xy)

distances = 6×6

0 1.4142 2.2361 6.3246 3.1623 5.6569 1.4142 0 1.0000 5.0990 2.0000 4.2426 2.2361 1.0000 0 4.1231 2.2361 3.6056 6.3246 5.0990 4.1231 0 5.0990 2.8284 3.1623 2.0000 2.2361 5.0990 0 3.1623 5.6569 4.2426 3.6056 2.8284 3.1623 0

[chicken vector]

pdist2 is definetly the best choice:

N = 5e3;
x = randi(10,1,N);
y = randi(10,1,N);

tic
squareFormZ = squareform(pdist([x' y']));
squareFormTime = toc;

tic;
xy = [x(:), y(:)];
distances = pdist2(xy, xy);
pdistTime = toc;

tic;
X = repmat(x, length(x), 1);
Y = repmat(y, length(y), 1);
deltaX = tril(x' - X, -1);
deltaY = tril(y' - Y, -1);
vectorZ = sqrt(deltaX(:).^2 + deltaY(:).^2);
vectorTime = toc;

fprintf("Time with squareform:        %.3f\n", squareFormTime)

Time with squareform:        0.448

fprintf("Time with pdist2:            %.3f\n", pdistTime)

Time with pdist2:            0.104

fprintf("Time with vectorisation:     %.3f\n", vectorTime)

Time with vectorisation:     0.290

[Sim]

Thanks a lot both @Image Analyst and @chicken vector!

I tried your code with the three variants and the test.dat set of points (here attached)

clear all; close all; clc
a = load('test.txt');
x = a(:,1);
y = a(:,2);

tic
squareFormZ = squareform(pdist([x' y']));
squareFormTime = toc;
tic;
xy = [x(:), y(:)];
distances = pdist2(xy, xy);
pdistTime = toc;
tic;
X = repmat(x, length(x), 1);
Y = repmat(y, length(y), 1);
deltaX = tril(x' - X, -1);
deltaY = tril(y' - Y, -1);
vectorZ = sqrt(deltaX(:).^2 + deltaY(:).^2);
vectorTime = toc;
fprintf("Time with squareform:        %.3f\n", squareFormTime)
fprintf("Time with pdist2:            %.3f\n", pdistTime)
fprintf("Time with vectorisation:     %.3f\n", vectorTime)

and I got the following results:

Time with squareform:        0.015
Time with pdist2:            0.012
Time with vectorisation:     70.236

I do not understand why, with this dataset, I got such a long computation time by using vectorisation..... :-)

[Sim]

More concisely:

clear all; close all; clc
a = load('test.txt');
tic
squareFormZ = squareform(pdist([a(:,1)' a(:,2)']));
squareFormTime = toc;
tic;
distances = pdist2(a,a);
pdistTime = toc;
fprintf("Time with squareform:        %.3f\n", squareFormTime)
fprintf("Time with pdist2:            %.3f\n", pdistTime)

I got this:

Time with squareform:        0.015
Time with pdist2:            0.017

[chicken vector]

Because it was designed for row vectors as per your question.

It is possible to generalise for both by doing instead:

X = repmat(x(:), 1, length(x))';
Y = repmat(y(:), 1, length(y))';

Anyway, when you have a built-in function in Matlab, it will almost always outperform hardcoded ones so just go with pdist2.

If you need to create a list with the indeces, see the method below to avoid loops, since that was the real time-consuming part of your code, rather than the distance method itself.

[Sim]

thanks a lot @chicken vector! However, the vectorisation method to create the list takes ages on my machine, and using the attached test.txt dataset....

clear all; close all; clc
a = load('test.txt');
x = a(:,1);
y = a(:,2);

tic;
xIdx = repmat(1 : length(x), length(x), 1);
yIdx = xIdx';
X = repmat(x, length(x), 1);
Y = repmat(y, length(y), 1);
deltaX = tril(x' - X, -1);
deltaY = tril(y' - Y, -1);
deltaX = deltaX((1 : size(deltaX, 1))' > (1 : size(deltaX, 2)));
deltaY = deltaY((1 : size(deltaY, 1))' > (1 : size(deltaY, 2)));
list = [xIdx((1 : size(xIdx, 1))' > (1 : size(xIdx, 2))) , ...
        yIdx((1 : size(yIdx, 1))' > (1 : size(yIdx', 2))) , ...
        sqrt(deltaX.^2 + deltaY.^2)]
vectorTime = toc;
fprintf("Time with vectorisation (just for creating the list):            %.3f\n", vectorTime)

I had to interrupt the code since it was taking ages... I am quite confused now....

[chicken vector]

Same reason as before, it was designed for row vectors.

I'd like to remark that built-in function is almost always the best choice.

As you can see is less prone to errors and also 30 % faster.

N = 1e4;
x = randi(10,N,1);
y = randi(10,N,1);

tic

xIdx = repmat(1 : length(x), length(x), 1);
yIdx = xIdx';
vectorIdx = (1 : size(xIdx, 1))' > (1 : size(xIdx, 2));

xy = [x(:), y(:)];
dist = pdist2(xy, xy);

list1 = [xIdx(vectorIdx) , ...
         yIdx(vectorIdx) , ...
         dist(vectorIdx)];
toc

Elapsed time is 2.443436 seconds.

But if you really want to use hardcoded method try this:

tic

xIdx = repmat(1 : length(x), length(x), 1);
yIdx = xIdx';
vectorIdx = (1 : size(xIdx, 1))' > (1 : size(xIdx, 2));

deltaX = tril(x(:) - repmat(x(:)', length(x), 1), -1);
deltaY = tril(y(:) - repmat(y(:)', length(y), 1), -1);

distVec = sqrt(deltaX(vectorIdx).^2 + deltaY(vectorIdx).^2);

list2 = [xIdx(vectorIdx) , ...
         yIdx(vectorIdx) , ...
         distVec];
toc

Elapsed time is 3.172369 seconds.

Quick check to see if everything works fine:

all(list1 == list2)

ans = 1×3 logical array
   1   1   1

[Sim]

Thanks a lot @chicken vector!!

Exactly this code you wrote works greatly for me!

If you post the following code as an Answer, I will accept it!

N = 1e4;
x = randi(10,N,1);
y = randi(10,N,1);
tic
xIdx = repmat(1 : length(x), length(x), 1);
yIdx = xIdx';
vectorIdx = (1 : size(xIdx, 1))' > (1 : size(xIdx, 2));
xy = [x(:), y(:)];
dist = pdist2(xy, xy);
distPdist = dist(vectorIdx);
list = [xIdx(vectorIdx) , ...
        yIdx(vectorIdx) , ...
        distPdist];
toc

[chicken vector]

You can accept this one, don't worry.

Thank you anyway!

[Image Analyst]

@Sim you can click the "Vote" icon for his answer below. Votes give the same 2 reputation points as Answers do.

[chicken vector]

Actually, according to this, accepted answer give +4 and votes +2.

Is this old?

[Image Analyst]

No, I must be thinking of the old way. Anyway, you can post a "final" fixed up program for a new answer and he can accept that.

[Sim]

Thanks a lot both @Image Analyst and @chicken vector!!

If @chicken vector you want to re-post an Answer as @Image Analyst suggested I will accept it :-) Meanwhile, obviously, I will upvote both :-)

Answer 3

[chicken vector]

You can build the indeces without for loop:

N = 5e2;
x = randi(10,1,N);
y = randi(10,1,N);

% Loop method:
tic;
k = 1;
for i = 1 : length(x)-1
    for j = i+1 : length(x)
        loopList(k,:) = [i j];
        k = k + 1;
    end
end
loopTime = toc;

% Vectorised method:
tic;
xIdx = repmat(1 : length(x), length(x), 1);
yIdx = xIdx';
vectorList = [xIdx((1 : size(xIdx, 1))' > (1 : size(xIdx, 2))) , ...
              yIdx((1 : size(yIdx, 1))' > (1 : size(yIdx', 2)))];
vectorTime = toc;

fprintf("Time with for loop:        %.3f seconds\n", loopTime)

Time with for loop:        0.988 seconds

fprintf("Time with vectorisation:   %.3f seconds\n", vectorTime)

Time with vectorisation:   0.009 seconds

You can also increase the speed for computing the distance with the following:

% Squareform method:
tic
squareFormZ = squareform(pdist([x' y']));
squareFormTime = toc;

% Vectorised method:
tic;
X = repmat(x, length(x), 1);
Y = repmat(y, length(y), 1);
deltaX = tril(x' - X, -1);
deltaY = tril(y' - Y, -1);
vectorZ = sqrt(deltaX(:).^2 + deltaY(:).^2);
vectorTime = toc;

fprintf("Time with squareform:        %.3f seconds\n", squareFormTime)

Time with squareform:        0.072 seconds

fprintf("Time with vectorisation:     %.3f seconds\n", vectorTime)

Time with vectorisation:     0.009 seconds

You can build your original list with the following wrapped up:

% Data:
x = [1 2 2 3 4 5];
y = [1 2 3 7 2 5];

% Initialise indeces:
xIdx = repmat(1 : length(x), length(x), 1);
yIdx = xIdx';

% Initialise elements distribution:
X = repmat(x, length(x), 1);
Y = repmat(y, length(y), 1);

% Compute distances:
deltaX = tril(x' - X, -1);
deltaY = tril(y' - Y, -1);

% Re-arrange to vector:
deltaX = deltaX((1 : size(deltaX, 1))' > (1 : size(deltaX, 2)));
deltaY = deltaY((1 : size(deltaY, 1))' > (1 : size(deltaY, 2)));

% Build lsit:
list = [xIdx((1 : size(xIdx, 1))' > (1 : size(xIdx, 2))) , ...
        yIdx((1 : size(yIdx, 1))' > (1 : size(yIdx', 2))) , ...
        sqrt(deltaX.^2 + deltaY.^2)]

list = 15×3

1.0000 2.0000 1.4142 1.0000 3.0000 2.2361 1.0000 4.0000 6.3246 1.0000 5.0000 3.1623 1.0000 6.0000 5.6569 2.0000 3.0000 1.0000 2.0000 4.0000 5.0990 2.0000 5.0000 2.0000 2.0000 6.0000 4.2426 3.0000 4.0000 4.1231