Matlab nearest neighbor / track points

Question

I have a set of n complex numbers that move through the complex plane from time step 1 to nsampl . I want to plot those numbers and their trace over time (y-axis shows imaginary part, x-axis the real part). The numbers are stored in a n x nsampl vector. However in each time step the order of the n points is random. Thus in each time step I pick a point in the last time step, find its nearest neighbor in the current time step and put it at the same position as the current point. Then I repeat that for all other n-1 points and go on to the next time step. This way every point in the previous step is associated with exactly one point in the new step (1:1 relation). My current implementation and an example are given below. However my implementation is terribly slow (takes about 10s for 10 x 4000 complex numbers). As I want to increase both, the set size n and the time frames nsampl this is really important to me. Is there a smarter way to implement this to gain some performance?

Example with n=3 and nsampl=2:

%manually create a test vector X
X=zeros(3,2); % zeros(n,nsampl)
X(:,1)=[1+1i; 2+2i; 3+3i];
X(:,2)=[2.1+2i; 5+5i; 1.1+1.1i]; % <-- this is my vector with complex numbers

%vector sort algorithm
for k=2:nsampl
    Xlast=[real(X(:,k-1)) imag(X(:,k-1))];          % create vector with x/y-coords from last time step
    Xcur=[real(X(:,k)) imag(X(:,k))];               % create vector with x/y-coords from current time step
    for i=1:size(X,1) % loop over all n points
        idx = knnsearch(Xcur(i:end,:),Xlast(i,:));  %find nearest neighbor to Xlast(i,:), but only use the points not already associated, thus Xcur(i:end,:) points
        idx = idx + i - 1;
        Xcur([i idx],:) = Xcur([idx i],:);          %sort nearest neighbor to the same position in the vector as it was in the last time step
    end
    X(:,k) = Xcur(:,1)+1i*Xcur(:,2);                %revert x/y coordinates to a complex number
end

Result:

X(:,2)=[1.1+1.1i; 2.1+2i; 5+5i];

Can anyone help me to speed up this code?

Answer 1

The problem you are tying to solve is combinatorial optimizartion which is solved by the hungarian algorithm (aka munkres). Luckily there is a implementation for matlab available for download. Download the file and put it either on your search path or next to your function. The code to use it is:

for k=2:size(X,2)
   %build up a cost matrix, here cost is the distance between two points.
   pairwise_distances=abs(bsxfun(@minus,X(:,k-1),X(:,k).'));
   %let the algorithm find the optimal pairing
   permutation=munkres(pairwise_distances);
   %apply it
   X(:,k)=X(permutation,k);
end