Vectorized for loop
2 views (last 30 days)
Show older comments
Good day
Okay, so my problem is the following. I have a massive number of x,y coordinates and need to do some calculations on them. I used for loops, but it takes forever. Below is a simple piece of code that explains what I am trying to do:
x = [ 0.2 0.3 0.2 0.1];
y = [ 0.8 0.4 0.5 0.2];
for i = 1:numel(x)
ind = 0;
for j = i+1:numel(x)-1
ind = ind +1;
d = abs(y(i)-y(j);
c(ind) = d + x(i)+x(j);
end
e = [e c];
end
The idea is to end up with a vector e (that is a function of x and y) that I can pass to another function.
So, my idea is to avoid a double for loop by having long (very long) vectors (sacrifice memory for runtime): thus to calculated d
vectori = [0.8 0.8 0.8 0.4 0.4 0.5];
vectorj = [0.4 0.5 0.2 0.5 0.2 0.2];
d = abs(vectori - vectorj);
The fact that j runs from i+1 to numel(x)-1 makes this also more challenging for me. How can I build vectori and vectorj for this application?
Regards, Barend.
0 Comments
Answers (3)
Robert Cumming
on 17 Nov 2011
before you spend a lot of time recoding - check where the code is slow. to do this use the profiler
profile on
% run you code
profile viewer
Check that you are pre-allocating your arrays (c and e for example are both growing in your loop above - this is very slow)
2 Comments
Robert Cumming
on 17 Nov 2011
preallocation is still valid
check where the time is with the profiler.
Jan
on 17 Nov 2011
As Robert has written already: pre-allocation is crucial. In addition all repeated calculatíon should be moved out of the inner loop.
x = [ 0.2 0.3 0.2 0.1];
y = [ 0.8 0.4 0.5 0.2];
nx = numel(x); % Once only!
c = zeros(1, nx); % Once only!
e = NaN(nx, nx);
for i = 1:nx
ind = 0;
xi = x(i); % Once only!
yi = y(i);
for j = i+1:nx-1
ind = ind + 1;
d = abs(yi - y(j));
c(ind) = d + xi + x(j);
end
e(1:ind, i) = c(1:ind); % Do not let e grow
end
e = e(isfinite(e)); % Remove the NaNs finally
Unfortunately you've posted a simplified computation only. In this example the inner loop can be vectorized easily:
c = abs(yi - y(i+1:nx-1)) + xi + x(i+1:nx-1);
I cannot guess if this will work in the original code also.
0 Comments
Andrei Bobrov
on 17 Nov 2011
e = nonzeros(tril(bsxfun(@plus,x1,x1')+abs(bsxfun(@minus,y1,y1')),-1));
variant with loop
nx = numel(x);
ny = nx - 1;
e = nan(nnz(tril(ones(nx-1),-1)),1);
id = ny;
k = 0;
for i1 = 1:ny
j1 = i1+1:ny;
id = id - 1;
k = k(end) + (1:id);
e(k) = abs(y(i1)-y(j1)) + x(i1)+x(j1);
end
2 Comments
Jan
on 17 Nov 2011
Creating the explicite index vector is slower, because Matlab checks the boundary for each element:
k = k(end) + (1:id); e(k) = ...
Faster:
e(k+1:k+id) = ...; k = k+id;
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!