vectorizing 4 nested for-loops in Matlab

Question

I'm writing a program for school and I have nested for-loops that create a 4-dimensional array (of the distances between two points with coordinates (x,y) and (x',y')) as below:

pos_x=1:20;
pos_y=1:20;
Lx = length(pos_x);             
Ly = length(pos_y);
Lx2 = Lx/2;
Ly2 = Ly/2;

%Distance function, periodic boundary conditions
d_x=abs(repmat(1:Lx,Lx,1)-repmat((1:Lx)',1,Lx));
d_x(d_x>Lx2)=Lx-d_x(d_x>Lx2);
d_y=abs(repmat(1:Ly,Ly,1)-repmat((1:Ly)',1,Ly)); 
d_y(d_y>Ly2)=Ly-d_y(d_y>Ly2);

for l=1:Ly
    for k=1:Lx
        for j=1:Ly
            for i=1:Lx
            distance(l,k,j,i)=sqrt(d_x(k,i).^2+d_y(l,j).^2);
            end
        end
    end
end

d_x and d_y are just 20x20 matrices and Lx=Ly for trial purposes. It's very slow and obviously not a very elegant way of doing it. I tried to vectorize the nested loops and succeeded in getting rid of the two inner loops as:

dx2=zeros(Ly,Lx,Ly,Lx);
dy2=zeros(Ly,Lx,Ly,Lx);
distance=zeros(Ly,Lx,Ly,Lx);

for l=1:Ly
    for k=1:Lx
        dy2(l,k,:,:)=repmat(d_y(l,:),Ly,1);
        dx2(l,k,:,:)=repmat(d_x(k,:)',1,Lx);
    end
end
distance=sqrt(dx2.^2+dy2.^2);

which basically replaces the 4 for-loops above. I've now been trying for 2 days but I couldn't find a way to vectorize all the loops. I wanted to ask:

1. whether it's possible to actually get rid of these 2 loops
2. if so, i'd appreciate any tips and tricks to do so. I have so far tried using repmat again in 4 dimensions, but you can't transpose a 4 dimensional matrix so I tried using permute and repmat together in many different combinations to no avail.

Any advice will be greatly appreciated.

---

thanks for the replies. Sorry for the bad wording, what I basically want is to have a population of oscillators uniformly located on the x-y plane. I want to simulate their coupling and the coupling function is a function of the distance between every oscillator. And every oscillator has an x and a y coordinate, so i need to find the distance between osci(1,1) and osci(1,1),..osci(1,N),osci(2,1),..osci(N,N)... and then the same for osci(1,2) and osci(1,1)...osci(N,N) and so on.. (so basically the distance between all oscillators and all other oscillators plus the self-coupling) if there's an easier way to do it other than using a 4-D array, i'd also definitely like to know it..

Answer 1

If I understand you correctly, you have oscillators all over the place, like this:

Then you want to calculate the distance between oscillator 1 and oscillators 1 through 100, and then between oscillator 2 and oscillators 1 through 100 etc. I believe that this can be represented by a 2D distance matrix, were the first dimension goes from 1 to 100, and the second dimension goes from 1 to 100.

For example

%# create 100 evenly spaced oscillators
[xOscillator,yOscillator] = ndgrid(1:10,1:10); 
oscillatorXY = [xOscillator(:),yOscillator(:)];

%# calculate the euclidean distance between the oscillators
xDistance = abs(bsxfun(@minus,oscillatorXY(:,1),oscillatorXY(:,1)')); %'# abs distance x
xDistance(xDistance>5) = 10-xDistance; %# add periodic boundary conditions
yDistance = abs(bsxfun(@minus,oscillatorXY(:,2),oscillatorXY(:,2)')); %'# abs distance y
yDistance(yDistance>5) = 10-yDistance; %# add periodic boundary conditions

%# and we get the Euclidean distance
euclideanDistance = sqrt(xDistance.^2 + yDistance.^2);