Is it possible to process the sparse matrix faster with vectorization instead of for loop?

Hello,

I need to calculate the Euclidean distances of the 3D point data I have. The distance of each point to all other points must be calculated. If the calculated value is less than 4, it should be defined as zero. I couldn't get any results using the following code with about 200000 points of data. Is it possible to process sparse matrix faster using vectorization instead of for loop?

The Code:

Point_Data = dlmread('E:\pointdata.txt'); % Point Data = xyz coordinate data of 3D points
len1 = length(Point_Data); 
distance = sparse(len1,len1);
for i1=1:len1
    for j1=1:len1
        d = (sqrt((Point_Data(i1,:)^2)-(Point_Data(j1,:)^2)));
        if d < 4
            distance(i1,j1) = d;
        end
    end
end

Answer 1

[James Tursa]

You need to re-evaluate how you are doing things. In the first place, you have the Euclidean calculation wrong. It should be this using your formulation:

d = sqrt(sum((Point_Data(i1,:)-Point_Data(j1,:)).^2));

or maybe just this more simply:

d = norm(Point_Data(i1,:)-Point_Data(j1,:));

And your test should be d>4 instead of d<4.

Secondly, your code iteratively fills up a sparse matrix. At each iteration when you fill an element, the entire data set must be copied into newly allocated memory to make room for the new element. This results in a serious drag on performance. For a large sparse matrix it would not be unusual to see iterative size increasing code to take days or even weeks of computation time. There are better ways of computing the Euclidean distances between sets of points. If you do end up using for loops, the general advice is to save the i,j,v triples off to the side and then build the sparse matrix all at once from these triples.

But most importantly, the size of the result may overwhelm your memory. Before you even try to build such a result, you need to figure out if the result is even something you can handle. In your particular case, what is the expected sparsity of the result? I.e., have you calculated the memory requirements of your expected result, and can you handle it? And what do you intend to do with this result downstream in your code? Maybe there are better ways of getting your problem solved rather than forming this potentially very large matrix explicitly.

[Musa Güngörürler]

Firstly, thank you. The problem I actually want to solve is: Importing data from Ansys software, updating my parameters using Matlab, and exporting to Ansys. Euclidian distances will be used to calculate local averages. E.g; The code works for 100 points, but I don't have enough memory to process about 200000 points of data. Do you have any idea if there is another way to solve it?

my_parameter = dlmread('C:\Users\user\my_par.txt'); % initial values of my parameters (Nx1)

from_ansys = dlmread('C:\Users\user\result.txt'); % Result matrix from Ansys Software (Nx1)

prm_total = sparse(distance*from_ansys); % distance matrix: The matrix I want to create - Euclidian distance matrix(NxN)

for i1=1:len1                                    % calculating mean value for each point - 
    mean_prm = prm_total(i1,1)/dist_total(i1,1); % dist_total(Nx1): Sum of the Euclidean distances of the neighboring points (distance < 4) for each element. 
                                                 % the distance matrix is used to calculate the local average for each point.
    my_value = from_ansys(i1,1)/(2*mean_prm);    

    if (my_value > my_parameter(i1,1))              %Updating the parameters according to the result.
        my_parameter(i1,1) = my_parameter(i1,1) + 1;
    elseif (my_value < my_parameter(i1,1))
        my_parameter(i1,1) = my_parameter(i1,1) - 1;
    end
    
end
dlmwrite('C:\Users\user\import_ansys.txt',my_parameter) %Exporting the updated parameters (new parameters for Ansys)

Answer 2

[Michael Van de Graaff]

So you're not gonna have enough memory, so that's gonna be a big problem. (actually, my computer did 12000 points in about 3 seconds, but nearly crashed for 20000 points, so maybe!)

However, if memory is not a problem, then the following should work:

tic
npoints = 1000; % your final matrix has this number of points squared, so 20000 points gonna hurt
data = (10*rand(npoints,3)); %making some sample points at random locations
permutation = flip(1:length(data(1,:)));
step0 = squeeze(permute(data,permutation)); %squeeze to make later reshape consistent
sze = size(step0);
step0 = reshape(squeeze(step0),[1,sze(1),sze(2)]); %if you only have 2 dimension, you need this to make the next line work 
step1 = data-step0; % this step breaks if 2 dimensions, you need the singleton dimension introduced above
step3 = sqrt(sum(step1.^2,2));

step4 =squeeze(step3); % symmetric matrix of difference, 
threshhold = 4;
sparse_distances = step4.*(step4>threshhold); %this is a symmetric matrix (npoints x npoints)
%  with zeros on diagonal (since every point is 0 away from itself)
%the matrix is symmetric since r_ij = r_ji ofcourse
toc

I found this helpful. https://www.mathworks.com/matlabcentral/answers/432560-how-can-i-compute-all-the-possible-differences-of-the-vectors-in-a-matrix

you should check my code with your loop (using a smaller sample) in case i messed up somehow

this will also work for points in Euclidean N-space, i think.

I also don't have experience with sparse matrices, so you might win there

[Michael Van de Graaff]

lol i skipped step2

Answer 3

[Walter Roberson]

Point_Data = readmatrixread('E:\pointdata.txt'); % Point Data = xyz coordinate data of 3D points

r = 4;
mdl = KDTreeSearcher(Point_Data);
Idx = rangesearch(mdl, Point_Data, r, 'SortIndices', false);

https://www.mathworks.com/help/stats/exhaustivesearcher.rangesearch.html

Idx is an m-by-1 cell array such that cell j (Idx{j}) contains an mj-dimensional vector of indices of the observations in Mdl.X that are within r units to the query observation Y(j,:). If SortIndices is true, then rangesearch arranges the elements of the vectors in ascending order by distance.

[t_r, t_c] = arrayfun(@(row) deal(repmat(row, length(Idx{row}), reshape(Idx{row},[],1)), (1:size(Idx,1)).', 'uniform', 0);
close_r = cell2mat(t_r);
close_c = cell2mat(t_c);

The array indices [close_r(K), close_c(K)] mark the places that are within the given range (and so should be zeroed.)