Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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}];

position[expr_] := 
  With[{positionData = 
     SortBy[#[[1, 1]] -> #[[All, 2]] & /@ 
         Extract[expr, #, Verbatim] -> # & /@ 
          Position[expr, _, Depth[expr]], First], 
       Min[Length /@ #[[2]]] &] // Dispatch}, 
   Replace[#, positionData] &];

poscoolist = position[coolist];

nnsite[k_, coolist_] := Module[{nncoolist},
   nncoolist = Nearest[Complement[coolist, {p[k]}], p[k]];
    Table[poscoolist[nncoolist[[i]]], {i, 1, Length[nncoolist]}]];

nearestlist = Table[nnsite[i, coolist], {i, 1, num}]; // AbsoluteTiming

Thread[Evaluate@Array[nnlabel, num] = nearestlist];

nnnsite[k_, coolist_] := Module[{nnncoolist},
   nnncoolist = 
      Delete[coolist, Partition[nnlabel[k], 1]], {p[k]}], p[k]];
    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.

share|improve this question
It looks like you're calling Nearest several times. Why not just call it once to get a precomputed NearestFunction, and use that to find the three nearest sites to any given point? The first one will be the site itself, the next two will be the points you want. It will be much more efficient. (I'd post an answer but I don't have access to Mathematica right now.) – Rahul Dec 10 '13 at 0:52
@RahulNarain Oh, thank you. NearestFunction! Let me see it. – matheorem Dec 10 '13 at 0:56
I don't understand what you're asking. – Rahul Dec 10 '13 at 1:55
I see, you want all the points at the nearest distance to $x$ in a single list. And then you want another list containing all the points at the second nearest distance. I suppose you could write a function that calls NearestFunction with larger and larger $n$ until it encounters a point at the third nearest distance... – Rahul Dec 10 '13 at 4:16
Nearest also supports a Nearest[data, x, {n, r}] syntax which returns $n$ nearest points within radius $r$, see: – Yves Klett Dec 10 '13 at 7:11

Here's a simple solution using a single precomputed NearestFunction; much faster than $O(n^2)$. I've written it assuming the sites are in a list ps, rather than embedded inside a function p, because I think this way is easier to generate and manipulate. You may want to modify the code as appropriate.

num = 5000;
ps = RandomInteger[{1, 100}, {num, 2}];

nf = Nearest[Table[ps[[i]] -> i, {i, Length[ps]}]] (* so it returns the index of the site *)

upToNthNearestSites[k_, 0] := {k} (* the "zeroth" nearest neighbour, i.e. itself *)
upToNthNearestSites[k_, n_] := Module[{pk, near, d},
  pk = ps[[k]];
  near = upToNthNearestSites[k, n - 1]; (* get the nearest neighbours up to order n-1 *)
  near = nf[pk, Length[near] + 1]; (* get one more; this is one of the nth nearest *)
  d = N@EuclideanDistance[pk, ps[[Last@near]]]; (* distance to the nth nearest *)
  nf[pk, {Infinity, d}] (* the solution is all the sites up to that distance *)
nthNearestSites[k_, n_] := Module[{pk, near0, near, d},
  pk = ps[[k]];
  near0 = upToNthNearestSites[k, n - 1];
  near = nf[pk, Length[near0] + 1];
  d = N@EuclideanDistance[pk, ps[[Last@near]]];
  near = nf[pk, {Infinity, d}];
  Complement[near, near0] (* same as above except remove neighbours closer than n *)

The nearest neighbours to the $k$th site are given by nthNearestSites[k, 1], the second nearest by nthNearestSites[k, 2], and so on. On my machine, even with a million points, the initial construction of nf takes a little over a second, and after that nthNearestSites[1, 2] takes negligible time.

Edit: I forgot that you want the neighbours of all the sites collected in a big list. Well, then you just do

nearestSites = nthNearestSites[#, 1] & /@ Range[Length[ps]];
nextNearestSites = nthNearestSites[#, 2] & /@ Range[Length[ps]];

On a hundred thousand sites, these take 2.9 and 6.8 seconds on my machine respectively. On a million, they will probably take a couple of minutes.


  1. You could just define nthNearestSites[k_, n_] := Complement[upToNthNearestSites[k, n], upToNthNearestSites[k, n - 1]], but that would end up evaluating the $(n-1)$th neighbours twice (as well as the $(n-2)$th, the $(n-3)$th, and so on). In the implementation above, it makes exactly $2n$ calls to the NearestFunction.

  2. I'm not too happy about having to put the N around EuclideanDistance. Unfortunately, NearestFunction doesn't accept something like $\sqrt2$ as the search radius.

share|improve this answer
Thank you for your answer! I am testing now, and I also found the EuclideanDistance problem, so I use N and add small offset 0.0000001 to it, I found this small offset is necessary, see my answer. – matheorem Dec 10 '13 at 8:57
I see your answer but it doesn't tell me why the offset is necessary. – Rahul Dec 10 '13 at 9:03
Well, you can try this ps={{0, 0}, {Sqrt[2], 0}, {0, Sqrt[2]}, {Sqrt[2], Sqrt[2]}}. Very odd,.but it's true, offset is necessary. – matheorem Dec 10 '13 at 9:49
Ah... Well, another workaround is that if you set ps = N@{{0, 0}, ...} instead then it works. – Rahul Dec 10 '13 at 10:45
Ok, are you considering to improve your answer with N? Though my method is identical to yours, but I make the label a little complex, so I'll accept yours, if there is no better method – matheorem Dec 10 '13 at 11:04

Here is a way to do it. At my first attempt I did not compile the functions and was a bit slower than your reported time. With the compiled function it is faster (at least on a 2Ghz/i7 Macbook Pro) : 5000 takes 23 secs at 10000 is 95 secs.

X := RandomInteger[{0, 100}];
n = 10000; data = Table[{X, X}, {i, 1, n}] ;

tStart = AbsoluteTime[];

"First, compute matrix of square distances (faster than taking the \
square roots) -- compiling is way faster"
dm = Compile[
   {{D, _Real, 2}, n},
    (D[[i, 1]] - D[[j, 1]])^2 + (D[[i, 2]] - D[[j, 2]])^2,
    {i, 1, n}, {j, 1, n}]
distMatrix = dm[data, n]; // AbsoluteTiming

"Now get the first and second closest distances :"
firstDistances = 
   Table[Min[Select[distMatrix[[i]], # != 0 &]], {i, 1, 
     n}]; // AbsoluteTiming
secondDistances = 
     Select[distMatrix[[i]], # != 0 && # != 
         firstDistances[[i]] &]], {i, 1, n}]; // AbsoluteTiming

"Finally find the neighbours :"
firstDistanceIndices = 
   Table[Position[distMatrix[[i]], firstDistances[[i]]], {i, 1, 
     n}]; // AbsoluteTiming
secondDistanceIndices = 
   Table[Position[distMatrix[[i]], secondDistances[[i]]], {i, 1, 
     n}]; // AbsoluteTiming

"Done. Time elapsed :"
AbsoluteTime[] - tStart

Clearly the code is $O(n^2)$ and the performance reflects that : doubling $n$ about quadruples the execution time. $n^2/10^6$ is a reasonable approximation to the actual time, at least for $n\leq10000$.

share|improve this answer
Thank you for your method! I also posted mine :) – matheorem Dec 10 '13 at 8:31

Update Based on Rahul's method. I clear some 'bug's (for example for situations when there are points have the same coordinates), and also did some improvements:

  1. using memorizing trick f[x_]:=f[x]=expr
  2. directly using Nearest[{e1->v1,e2->v2,...},x]

Testing data

num = 200;
ps = N@RandomInteger[{1, 100}, {num, 2}];
pslabeled = Thread[ps -> Range[num]];(*generate a rulelist with everypoint labeled*)

The nearestsitemodule

    nearestsitemodule[pslabeled_] := 
      Module[{ps, labeldispatch, reverselabeldispatch, offset},
       labeldispatch = Dispatch[pslabeled];
       reverselabeldispatch = Dispatch[Reverse /@ pslabeled];
       offset = 0.0000001;
       nf = Nearest[pslabeled]; (*generate nearest function*)


       (*the "zeroth" nearest neighbour,i.e.itself*)
       upToNthNearestSites[coorulelabel_,0] := {Last@coorulelabel}; 

gives the label list of up to nth nearest neighbours relative to coorulelabel
           using memorizing trick f[x_]:=f[x]=expr*)
       upToNthNearestSites[coorulelabel_, n_] := 
        upToNthNearestSites[coorulelabel, n] = Module[{coo, near, d},
          coo = First@coorulelabel;
          near = upToNthNearestSites[coorulelabel, n - 1];(*get the nearest neighbours up to order n-1*)
          near = nf[coo, Length[near] + 1];(*get one more;
          this is one of the nth nearest*)
          d = N[EuclideanDistance[coo, Replace[Last[near], reverselabeldispatch]] + offset];
(*d is distance to the nth nearest, notice a small offset is necessary!!!!*)
          nf[coo, {All, d}](*the solution is all the sites up to that distance*)];

       (*nthNearestSites[coorulelabel,n] gives the label list of nth nearest neighbours to coorulelabel*)
       nthNearestSites[coorulelabel_, n_] := Module[{near0, near, d},
         near0 = upToNthNearestSites[coorulelabel, n - 1];
         near = upToNthNearestSites[coorulelabel, n];
         Complement[near, near0](*get label list of the nth nearest neighbour to  coorulelabel*)];


nthNearestSites[#, 1] & /@ pslabeled; // AbsoluteTiming
nthNearestSites[#, 2] & /@ pslabeled; // AbsoluteTiming

The performance for 5000 points is 0.58s and 0.79s

The performance for 10000 points is 1.32s and 1.27s

Rahul claimed "On a hundred thousand sites, these take 2.9 and 6.8 seconds". I think his computer must be much faster than mine. Because my test shows, both code's running time is comparable, but as I said, I deal with points with the same coordinates.

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.