Vectorizing code for runtime optimization (nested loops)

I am trying to reduce the runtime of the following code:

%Parameters
%Number of Points
T = 10000;
%Array of 3D-Points
v = rand(T,3);
%kth nearest neighbor
k = 101;

%Pre-allocating
d = zeros(1, T-1);
d_xz = zeros(1, T-1);
d_yz = zeros(1, T-1);
d_z = zeros(1, T-1);
e_k = zeros(1, T);

%Initializing counting variables
count_xz = 0;
count_yz = 0;
count_z = 0;

for t = 1:T
    index_aux = 0;
    for z = 1:T
        %Calculates distances from point in time z to points in all other times

        if z == t
            %checks if t = t' which has to be avoided. index_aux helps correcting the index after this case is ignored
            index_aux = 1;

        else

            %v_Xb are auxiliary variables which are used to
            %calculate KNN distance "d(z)" and marginal points
            %distance "d_xz(z), d_yz(z) and d_z(z)".
            
            v_xb = abs(v(t,1)-v(z,1));
            v_yb = abs(v(t,2)-v(z,2));
            v_zb = abs(v(t,3)-v(z,3));

            if v_xb >= v_zb
                d_xz(z-index_aux) = v_xb;
                if v_xb >= v_yb
                    d(z-index_aux) = v_xb;
                end
            else
                d_xz(z-index_aux) = v_zb;
            end

            if v_yb >= v_zb
                d_yz(z-index_aux) = v_yb;
                if v_yb >= v_xb
                    d(z-index_aux) = v_yb;
                end
            else
                d_yz(z-index_aux) = v_zb;
            end

            if v_zb >= v_xb && v_zb >= v_yb
                d(z-index_aux) = v_zb;
            end
            d_z(z-index_aux) = v_zb;
            
            %The following lines do the same calculation but runs slower than the above lines (at least for the values given for T)
            %d(z) = max([abs(v(t,1)-v(z,1)), abs(v(t,2)-v(z,2)), abs(v(t,3)-v(z,3))]);
            %d_xz(z) = max([abs(v(t,1)-v(z,1)), abs(v(t,3)-v(z,3))]);
            %d_yz(z) = max([abs(v(t,2)-v(z,2)), abs(v(t,3)-v(z,3))]);
            %d_z(z) = abs(v(t,3)-v(z,3));
        end
    end
    
    d = sort(d);
    e_k(t) = d(k);
    
    for z = 1:(T-1)
        %Calculating the number of points for 1 given time z.
        %Calculating points for each marginal vector.
        
        if d_xz(z) < e_k(t)
            count_xz = count_xz+1;
        end
        
        if d_yz(z) < e_k(t)
            count_yz = count_yz+1;
        end
        
        if d_z(z) < e_k(t)
            count_z = count_z+1;
        end
    end
   
end

The code basically calculates the kth nearest neighbor distance using Chebyshev metric and then counts how many points made out of marginal vectors are found for each "t" that are less than this distance.

I am trying to figure out how to vectorize the nested for loops and how in general optimize this code. For the distances calculation, an alternative is written as a comment but this takes longer to run. Also I have tried using the MATLAB built in knn algorithm with the same metric, but it also takes longer. Finally for the last part, I also tried using cumsum and other alternatives with other built in MATLAB functions but my primitive way of counting is faster (which is only true when T is of low order of magnitude).

Many thanks in advance!

---------------------------------------------------------------------------------------------------------------------------

EDIT: After re-checking the suggestions made by Jan I ended up with this version which is three times faster than my former code.

%Defining parameters
v = rand(10000, 3);
T = 10000;
k = 101;
count_xz = 0;
count_yz = 0;
count_z  = 0;

%Pre-allocating
d    = zeros(T-1, 1);
d_xz = zeros(T-1, 1);
d_yz = zeros(T-1, 1);
d_z  = zeros(T-1, 1);
e_k  = zeros(k, 1);
v_xb = zeros(T-1, 1);
v_yb = zeros(T-1, 1);
v_zb = zeros(T-1, 1);

for t = 1:T
    
    %Ignoring t == index case
    index    = true(1, T);
    index(t) = false;
    
    %Calculating distances between coordinates of all points
    v_xb = abs(v(t, 1) - v(index, 1));
    v_yb = abs(v(t, 2) - v(index, 2));
    v_zb = abs(v(t, 3) - v(index, 3));

    %Calculating the Chebyshev distance between the points and between
    %the projections of the points in the XZ, XY and Z planes
    d_xz = max(v_xb, v_zb);
    d_yz = max(v_yb, v_zb);
    d_z  = v_zb;
    d    = max(d_xz, v_yb);
    
    %Finding Knn-Distance between the points
    e_k = mink(d,k);

    %Counting the number of distances between the points in the different
    %projections that are less than the Knn-Distance
    count_xz = count_xz + sum(d_xz(1:T-1) < e_k(end));
    count_yz = count_yz + sum(d_yz(1:T-1) < e_k(end));
    count_z  = count_z  + sum(d_z(1:T-1)  < e_k(end));
end

Please feel free to suggest any other ways I could optimize this code. All execution time improvements count (since I have to run this code for several data sets independently)

Thank you!

[Jan]

@J. Guerrero: This is your first question. You have formatted the code properly, provided your solution, input data and explained what you have tried so far. The question is nice and clear. Thanks. This is my "question of the week"! :-) +1

Answer 1

[Jan]

v_xb = abs(v(t,1)-v(z,1));

This looks symmetric. Is the output symmteric (I cann test this currently). Then it should be sufficient to run the innerloop over a "triangle" only:

for z = t+1:T

This part can be simplified - perhaps accelerated:

 for z = 1:(T-1)
        %Calculating the number of points for 1 given time z.
        %Calculating points for each marginal vector.
        
        if d_xz(z) < e_k(t)
            count_xz = count_xz+1;
        end
        
        if d_yz(z) < e_k(t)
            count_yz = count_yz+1;
        end
        
        if d_z(z) < e_k(t)
            count_z = count_z+1;
        end
 end

to:

     count_xz = count_xz + sum(d_xz(1:T-1) < e_k(t));
     count_yz = count_yz + sum(d_yz(1:T-1) < e_k(t));
     count_z  = count_z  + sum(d_z(1:T-1)  < e_k(t));

I assume the best help will be to vectorize the inner loop:

     index    = true(1, T);
     index(t) = false;
     v_xb     = abs(v(t, 1) - v(index, 1));
     ...

Unfortunately I do not understand your sorting and comparisons directly (some comments might be a good idea), so I have to run it to find out, how an equivalent vector comparison must look like.

[J. Guerrero]

Thank you for your time and effort, I appreciate it!

Going through your recommendations:

Indeed the code was sped up by around 10% through counting through

     count_xz = count_xz + sum(d_xz(1:T-1) < e_k(t));
     count_yz = count_yz + sum(d_yz(1:T-1) < e_k(t));
     count_z  = count_z  + sum(d_z(1:T-1)  < e_k(t));

which is already a good start. I realized that, since my intention is to find the kth-nearest neighbor distance, I do not need to sort the whole distance matrix "d" and can instead use the mink function as

     e_k = mink(d,k);
     %instead of first sorting and then looking for sorted(k)

I vectorized the inner loop and used the max function instead of 'manually' comparing each terms to find the Chebyshev distances. This cut the runtime by half at least which is amazing!

Please have a look at the latest version of my code, perhaps you can still find optimizations of my vectorization.

Thank you once again for your time and patience!

[Jan]

As far as I can see,

count_xz = count_xz + sum(d_xz(1:T-1) < e_k(end))

can be simplified to

count_xz = count_xz + sum(d_xz < e_k(end))

and for _yz and _z also.

A further tiny improvement:

index = true(1, T);  % Create it once only
for t = 1:T    
    %Ignoring t == index case
    index(t) = false;
    ...
    index(t) = true;  % Reset this index
end

There is no need to pre-allocate arrays, if you do not use them later. So omit:

%Pre-allocating
d    = zeros(T-1, 1);
d_xz = zeros(T-1, 1);
d_yz = zeros(T-1, 1);
d_z  = zeros(T-1, 1);
e_k  = zeros(k, 1);
v_xb = zeros(T-1, 1);
v_yb = zeros(T-1, 1);
v_zb = zeros(T-1, 1);

It does not waste a lot of time, but it is not useful.

[J. Guerrero]

Thank you very much again, Jan! These last comments improved the runtime by around 10%ish. I believe after those modifications there is anything else to do to it. If you come up with anything else please feel free to mention it.