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

Clear[nnlabel];
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). nearestlist and nextnearestlist give the result.

Notice to use my code, supply every coordinate p[i] of points.

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.