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

2 views (last 30 days)
Musa Güngörürler
Musa Güngörürler on 18 Jan 2022
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;

Answers (3)

James Tursa
James Tursa on 18 Jan 2022
Edited: James Tursa on 18 Jan 2022
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.
  1 Comment
Musa Güngörürler
Musa Güngörürler on 18 Jan 2022
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;
dlmwrite('C:\Users\user\import_ansys.txt',my_parameter) %Exporting the updated parameters (new parameters for Ansys)

Sign in to comment.

Michael Van de Graaff
Michael Van de Graaff on 18 Jan 2022
Edited: Michael Van de Graaff on 18 Jan 2022
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:
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
I found this helpful.
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

Walter Roberson
Walter Roberson on 18 Jan 2022
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);
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.)




Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by