Suppose we have a point set, the points are labeled 1,2,3,....n. And p[i] is the coordinate of the point. Now I want the label of nearest and next nearest point (maybe several) corresponding every point in this point set.

For example:

    p[1] = {0, 0};
    p[2] = {0, 1};
    p[3] = {1, 0};
    p[4] = {0, 5};

and I want a list for nearest neighbour like this:

    {{2, 3}, {1}, {1}, {2}}

it means that p[2] and p[3] is nearest to p[1], p[1] is nearest to p[2], etc. 

And I need a similar next nearest neighbour list.

I write the following code: 

    num = 5000;
    Do[p[i] = RandomInteger[{1, 100}, 2], {i, 1, num}];
    
    coolist = Table[p[i], {i, 1, num}];
    
    Clear[position];
    position[expr_] := 
      With[{positionData = 
         SortBy[#[[1, 1]] -> #[[All, 2]] & /@ 
            GatherBy[
             Extract[expr, #, Verbatim] -> # & /@ 
              Position[expr, _, Depth[expr]], First], 
           Min[Length /@ #[[2]]] &] // Dispatch}, 
       Replace[#, positionData] &];
    
    poscoolist = position[coolist];
    
    Clear[nnsite];
    nnsite[k_, coolist_] := Module[{nncoolist},
       nncoolist = Nearest[Complement[coolist, {p[k]}], p[k]];
       Flatten@
        Table[poscoolist[nncoolist[[i]]], {i, 1, Length[nncoolist]}]];
    
    nearestlist = Table[nnsite[i, coolist], {i, 1, num}]; // AbsoluteTiming
    
    Thread[Evaluate@Array[nnlabel, num] = nearestlist];
    
    Clear[nnnsite];
    nnnsite[k_, coolist_] := Module[{nnncoolist},
       nnncoolist = 
        Nearest[Complement[
          Delete[coolist, Partition[nnlabel[k], 1]], {p[k]}], p[k]];
       Flatten@
        Table[poscoolist[nnncoolist[[i]]], {i, 1, Length[nnncoolist]}]];
    
    nextnearestlist = 
       Table[nnnsite[i, coolist], {i, 1, num}]; // AbsoluteTiming

the function `position` in the above code is provide by Mr.Wizard (see [here][1]). `nearestlist` and `nextnearestlist` give the result.

For 5000 points, it takes 1 minutes. But for 10000 points, it takes 6 minutes. Quite long! 

I feel that my code is so simple, there must be better ways that are more efficient.


  [1]: http://mathematica.stackexchange.com/a/25592/4742