Get rid of nested for loops?

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!

finalplot=zeros(2452,2);
norm=zeros(2452,1);

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

else

      norm(1+round(sqrt((i-a)^2+(j-b)^2)),1)=norm(1+round(sqrt((i-a)^2+(j-b)^2)),1)+1;

      finalplot(1+round(sqrt((i-a)^2+(j-b)^2)),2)=finalplot(1+round(sqrt((i-a)^2+(j-b)^2)),2)+(v.vx(a,b)*v.vx(i,j)+v.vy(a,b)*v.vy(i,j))/(sqrt((v.vx(a,b))^2+(v.vy(a,b))^2)*sqrt((v.vx(i,j))^2+(v.vy(i,j))^2));

end

      end
  end

      end
  end

for i=1:2452
finalplot(i,1)=i;
finalplot(i,2)=finalplot(i,2)/norm(i,1);
end

Answer 1

[Teja Muppirala]

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

Notes:

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

to

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.

Answer 2

[Roger Stafford]

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

 1+round(sqrt((i-a)^2+(j-b)^2))

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)