Is it possible to process the sparse matrix faster with vectorization instead of for loop?
6 次查看(过去 30 天)
显示 更早的评论
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
0 个评论
回答(3 个)
James Tursa
2022-1-18
编辑:James Tursa
2022-1-18
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.
Michael Van de Graaff
2022-1-18
编辑:Michael Van de Graaff
2022-1-18
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
Walter Roberson
2022-1-18
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.)
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Matrix Indexing 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!