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

The challenge

Given $n$ points in $\mathbf{R}^3$ and a cutoff distance, $d$, I would like to determine

  1. all pairs of points that are separated by less than the given cutoff distance, $d$, and
  2. all pairs of points that can be connected by less than 5 segments taken from the list of pairs from above.

While this is straight-forward for small $n$, for large $n$ the simplistic approach given below is slow (searching for the pairs) and memory consuming (GraphDistanceMatrix seems to return a dense matrix). How could this problem be solved efficiently for $n$ of order 30000?

Inefficient sample code:

Generate sample data:

n = 1000; (* desired: n=20000 *)
pos = RandomReal[{0, n^(1/3)-1}, {n, 3}];

Find close points (ad 1):

d = 0.8;
connected = Select[Subsets[Range[n], {2}], (Norm[pos[[#[[1]]]] - pos[[#[[2]]]]] < d)&];

Pairs of points closer than 5 leaps apart (ad 2):

v = Union[Flatten[connected]];
g = Graph[v,UndirectedEdge @@@ connected];
gdm = GraphDistanceMatrix[g, 4];
ind = 
  Select[Position[gdm, Alternatives[1, 2, 3, 4], {2}], (#[[1]] < #[[2]])&] /. 
    Thread[Range[Length[v]] -> v];
share|improve this question
to all: I found it difficult to accept one of the answers. All three contributions were enlightening and helpful. I have eventually accepted szabolcs answer, because I enjoyed being reminded of the direct use of adjacency matrices and the fact that the built-in tools are not always superior. However, I want to point out that belisarius was the first to introduce the use of the NearestFunction, which all other answers have retained (even though it accepts all points at distances <=d, instead of <d), and halmir was the first presenting a complete answer. Thank you all very much! – dakta Jan 28 '14 at 11:28
up vote 9 down vote accepted
n = 20000;(*desired:n=20000*)SeedRandom[42];
pos = RandomReal[{0, n^(1/3) - 1}, {n, 3}];
d = 0.8;
f = Nearest[pos -> Range[n]];
pos1 = f[#, {Infinity, d}] & /@ pos;

Starting from belisarius's pos1, we can efficiently solve task 2 like this:

  1. build the adjacency matrix of the graph as a sparse array:

    am = SparseArray[Flatten[Thread /@ Thread[{pos1, Range[n]}], 1] -> 1];
  2. the $k^\text{th}$ power of the adjacency matrix $A$ will tell us which nodes are connected by precisely $k$ hops. Since this particular adjacency matrix am has all ones on the diagonal, its $k^\text{th}$ power will tell us which nodes are connected by $k$ or fewer hops:

    am5 = Unitize@MatrixPower[am, 5];

All this takes 0.14 seconds on my machine when starting with n=20000 as belisarius did:

 am = SparseArray[Flatten[Thread /@ Thread[{pos1, Range[n]}], 1] -> 1];
 am5 = Unitize@MatrixPower[am, 5];

(* ==> {0.141265, Null} *)

Use ArrayRules to extract the pairs (edges) from a sparse adjacency matrix.

share|improve this answer
the fifth power of the adjacency matrix will tell us which nodes are connected by no more than five steps : AdjacencyGraph@ Unitize@MatrixPower[AdjacencyMatrix[GridGraph[{2, 2}]], 2] – Dr. belisarius Jan 27 '14 at 19:58
@belisarius It seems I'm particularly stupid today. Does it look fine now? – Szabolcs Jan 27 '14 at 20:35
Yep :) +1 Now :) – Dr. belisarius Jan 27 '14 at 20:46
@user11977 Be careful, it's not superfluous. If $A$ is the adjacency matrix, then $A^k$ is the adjacency matrix for connections possible through precisely $k$ hops, not less or more. So we need to sum up all of $A, A^2, \ldots, A^k$ to get the connections possible through $k$ hops or less. – Szabolcs Jan 27 '14 at 22:03
@szabolcs You are certainly right. However, as far as I can tell, "am" from above is the adjacency matrix plus the identity matrix. Thus, the summation is implicit, isn't it? – dakta Jan 27 '14 at 23:51

The first part is fairly easy:

n = 20000;(*desired:n=20000*)
pos = RandomReal[{0, n^(1/3) - 1}, {n, 3}];
d = 0.8;
f = Nearest[pos];
pos1 = f[#, {Infinity, d}] & /@ pos;


The second part, using Combinatorica:

<< Combinatorica`
kNeighborhoods[pos_, d_, dist_] := Module[{f, pos1},
  f = Nearest[pos -> Automatic];
  pos1 = (Rest@f[#, {Infinity, d}]) & /@ pos;
  Neighborhood[pos1, #, dist] & /@ Range@Length@pos
d = 0.8;
dist = 4;
n = 20000;
pos = RandomReal[{0, n^(1/3) - 1}, {n, 3}];
kNeighborhoods[pos, d, dist]; // Timing
 {1.006250, Null}
share|improve this answer
I was not aware of this use of the NearestFunction. Is this gemstone documented? Furthermore, am I right that also points with distance equal to d are returned? – dakta Jan 27 '14 at 17:14
@user11977 Docs are here – Dr. belisarius Jan 27 '14 at 17:26
Fair enough. The documentation of NearestFunction does not mention this calling sequence leading me to assume it is undocumented ... Anyway, thank you very much. – dakta Jan 27 '14 at 19:09

I modified belisarius method to get first part:

n = 20000; SeedRandom[42];
pos = RandomReal[{0, n^(1/3) - 1}, {n, 3}];
d = 0.8;
rule = Thread[pos -> Range[n]];
f = Nearest[rule];
pos1 = f[#, {Infinity, d}] & /@ pos;

second part:

res = Union[
     With[{i = First[#2]}, Thread[i -> Select[#1, # > i &]]] &, 
g = Graph[res, DirectedEdges -> False];

res2 = Union[
    Sort /@ Flatten[
      Thread[{#, AdjacencyList[g, #, 4]}] & /@ VertexList[g], 
      1]]; // AbsoluteTiming

{1.075134, Null}

share|improve this answer
Thank you very much! I knew there had to be something better than GraphDistanceMatrix. AdjacencyList escaped me ... – dakta Jan 27 '14 at 21:01

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.