I am trying to use the following code, which has way too many nested for loops. However, I'm a newbie, and don't have an obvious way to make this code more efficient. Does anyone have any suggestions?
Basically, I have a structural array (v) with four 2441x226 components: vx, vy, x, and y. vx and vy contain the x and y components of a vector field, and x and y label the position of each vector in that field. The code calculates how the correlation between the vectors (finalplot(:,2)) decays with Cartesian distance (finalplot(:,1)).
Thanks for your help!
for b=1:226 for a=1:2441
for j=1:226 for i=1:2441
if isnan(v.vx(a,b)*v.vx(i,j)+v.vy(a,b)*v.vy(i,j)) %this line just gets rid of all the NaN's that are in the vector field
for i=1:2452 finalplot(i,1)=i; finalplot(i,2)=finalplot(i,2)/norm(i,1); end
A few observations:
1. A lot of the operations you are doing can me written more efficiently as matrix operations (dot products, etc.), or operations on entire matrices at once, without needing to loop over elements. For example, normalizing the vectors does not need any explicit loops at all.
2. You are checking every point with every other point. But checking point A with point B is the same as checking point B with point A, so you really only need to do half the calculations.
3. Not pertinent to speeding this up, but "norm" is the name of a very useful MATLAB command. It is generally a good idea to avoid calling your variables the same thing as a built-in MATLAB function.
This gives the same results as your code, but is much faster:
V = [v.vx(:) v.vy(:)]; % Line up the vectors in a big vertival matrix V1 = bsxfun(@rdivide,V,sqrt(v.vx(:).^2+v.vy(:).^2)); %Normalize it P = [v.x(:) v.y(:)]; % Matrix of all the (x,y) points
finalplot=zeros(NumOut,1); loccount=zeros(NumOut,1); %What you called "norm"
notnan = find(~isnan(v.vx(:)) & ~isnan(v.vy(:))); % Find the NaNs
P = P(notnan,:); %Get rid of the NaN locations from the start V1 = V1(notnan,:);
N = numel(notnan); for n = 1:N %Loop only once locations = 1+round(sqrt(sum(bsxfun(@minus,P(n,:),... P(n+1:N,:)).^2,2))); loccount = loccount + accumarray(locations,1,[NumOut 1]); finalplot = finalplot + accumarray(locations,... V1(n+1:N,:)*V1(n,:)',[NumOut 1]); end finalplot = [(1:numel(finalplot))' finalplot./loccount]; finalplot(1,2) = 1; %Fix the first value
1. The code above assumes your x and y data is like NDGRID, not like MESHGRID. Compare
[x,y] = ndgrid(1:5) [x,y] = meshgrid(1:5)
If your data is like MESHGRID, then you'd need to change:
P = [v.x(:) v.y(:)]; % Matrix of all the (x,y) points
P = [reshape(v.x',,1) reshape(v.y',,1)];
to make it match your code exactly.
2. The remaining one FOR loop is highly parallelizable. Every iteration can be done independently. If you have the Parallel Computing Toolbox, or MATLAB Distributed Computing Server, and access to a computer with multiple CPU cores or a computing cluster, you can change the FOR loop into a PARFOR loop and run the calculation in parallel, greatly reducing calculation time.
The nested loops are difficult to avoid for this type of procedure. However, there are other aspects of your computation that can be made more efficient. There is too much repetition of operations, in particular too many square root operations. It is better to process the vector field first so that all vectors are of unit magnitude. I have placed the results of doing so in arrays nx and ny. I put a NaN in vx when needed to inhibit calculations later. Also you should only compute the index
once for each a,b,i,j combination and save it (in t2) rather than four times for each combination.
nx = zeros(2441,226); ny = zeros(2441,226); for b = 1:226 for a = 1:2441 if isnan(v.vx(a,b)) | isnan(v.vy(a,b)) nx(a,b) = NaN; else t1 = sqrt((v.vx(a,b))^2+(v.vy(a,b))^2); nx(a,b) = v.vx(a,b)/t1; ny(a,b) = v.vy(a,b)/t1; end end end finalplot = zeros(2452,2); norm = zeros(2452,1); for b = 1:226 for a = 1:2441 if ~isnan(nx(a,b)) for j = 1:226 for i = 1:2441 if ~isnan(nx(i,j) t2 = 1+round(sqrt((i-a)^2+(j-b)^2)); norm(t2) = norm(t2)+1; finalplot(t2,2) = finalplot(t2,2) + nx(a,b)*nx(i,j)+ny(a,b)*ny(i,j); end end end end end end finalplot(:,1) = 1:2452; finalplot(:,2) = finalplot(:,2)./norm(:,1)