Improvement of point cloud distance calculation in 2D

Hello everyone,

I have following question:

Given is an array of n- datapoints, with x and y coordinates for every point.

I know want to calculate the distance from every point to every other point.

For now the most efficient way for me is:

A=rand(5000,2);

Dista = squeeze(sqrt(sum((A - reshape(A',1,2,[])).^2,2)));

With increasing number of points this will take relativ long, and since i have to reapeat it inisde a loop it will in take to long in summary for me.

Theoretical, this could be optimized, since in my code above I calculate the distance from every point to everypoint, which is a waste of approximatly 50% (for an infinite number of points) of the calculations, since 50% of the calculations are done twice since the distcance from point A to point B is equal to B to A.

But how could I get this constraint inside my code, in best case without loops, since I guess, that this will slow it down in total.

All i have to do is to mirror the Dista matrix along the diagonal, but how do i manage that my code only calculates the uper(or lower) triangle of the matrix, in vectorized form.

Many thanks in advance.

Best regards

Marc

[Tommy]

You may be interested:

function [T1, T2, T3, T4] = test
    A = rand(10000,2);
    T1 = timeit(@method1);
    T2 = timeit(@method2);
    T3 = timeit(@method3);
    T4 = timeit(@method4);
    
    function dist = method1
        dist = squeeze(sqrt(sum((A - reshape(A',1,2,[])).^2,2)));
    end
    
    function dist = method2
        dist = sqrt(triu((A(:,1)-A(:,1)')).^2+triu((A(:,2)-A(:,2)')).^2);
        dist = dist+dist';
    end
    
    function dist = method3
        dist = sqrt(((A(:,1)-A(:,1)')).^2+((A(:,2)-A(:,2)')).^2);
    end
    
    function dist = method4
        dist = squareform(pdist(A));
    end
end

Output:

>> [T1, T2, T3, T4] = test

T1 =

    0.6294 % sec


T2 =

    1.0634  % sec


T3 =

    0.3852  % sec


T4 =

    1.0522  % sec

Using triu (or similarly tril) seems a bit pointless as long as every operation (subtraction, raise to 2nd power, addition, sqrt) is applied to every element in the matrix, even if half of those elements are 0. Using two loops or linearizing half of the matrix and then reshaping are two ways I can think of to really avoid half of the triangle, but I'm not sure you'd gain anything.

But maybe the behavior would change with larger matrices and you're thinking higher than 10,000.

[Marc Laub]

I tried using a single parfoor loop with a pdist2 on triu matrizes. It just made it worse...For a small number of points it took like 17 sec, while my method above took like 0.8 sec, for more points I ran out of memory.

Your method 3 saves me like 60% of time for a 10k points, and up to 82%-87% for 40k, 60k points.

But i dont get what you mean with your "linearizing half of the matrix and then reshaping ". Could you tell me what you mean? Mybe it helps for larger matrizes, but i actually dont think so, because in my method the reshape costs like 60-80% of the time, but maybe its worth a try

Thank you very much so far

[Tommy]

Essentially, I was thinking that if you have some symmetric matrix

>> A = rand(3);
>> A = A+A'

A =

    1.5844    0.9952    1.3345
    0.9952    1.6983    1.6917
    1.3345    1.6917    1.4863

then you can get the upper triangle:

>> triu(A,1)

ans =

         0    0.9952    1.3345
         0         0    1.6917
         0         0         0

but you will be dragging around extra 0s, and you can only avoid running your calculations on the 0s if you tell MATLAB exactly which values to work with:

for i = 1:size(A,1)
    for j = i+1:size(A,1)
        disp(A(i,j))
    end
end

%{
0.9952

1.3345

1.6917
%}

or you pull out the relevant values and store them differently before doing vectorized calculations:

>> A(triu(true(size(A)),1))

ans =

    0.9952
    1.3345
    1.6917

or possibly through some other method.

Of course, you are starting with an Nx2 vector, maybe you can bypass the NxN matrix completely. You could probably use nchoosek...

If the third method is good enough, then that also works!

[Marc Laub]

I got rid of the sqrt, since I am only interested in the order of the nearest points and not of the actual distance. This saves a huge amount of time.

Considering linearizing, I am not sure if MatLab is not considering it anayways, when doing operations like A-A'.

Those kind of operations are always symmetrical and I would wonder if a vector operation optimized programm like MatLab wouldnt take advantage of this. Maybe somobody knows this?

A = rand(10000);
B = rand(10000);
tic
C=(A-B').^2;
toc
clear C

tic
C=(A-A').^2;
toc

Second case has a about 20% time improved compared to case 1, so i guess that the symmtrie of the operation is already implemented

[Tommy]

Good point about the sqrt.

Try switching the order, you might find that the second line is always quicker than the first.

[Marc Laub]

I switched the order and now the asymetric one is faster.

Seems like

timeit(@()(A-B').^2)

bypasses this probelm.

timeit(@()(A.*B')

clearly shows that in the multiplacting case Matlabs gets some advantage of the symmetrie, here the symmetric one is alwys faster, independend on the order. In the upper case, there is nearly no advantage, maybe 0.001% and sometims the asymatric one beeing the faster one.

So the question is how to tell Matlab to get advantage of symmetric operations since the symmetrie is not automatically detected