Pairwise distances of points and classification using GraphDistanceMatrix

Question

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

Answer 1

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:

AbsoluteTiming[
 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.

Answer 2

The first part is fairly easy:

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

Edit

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}
*)

Answer 3

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[
   Flatten[MapIndexed[
     With[{i = First[#2]}, Thread[i -> Select[#1, # > i &]]] &, 
     pos1]]];
g = Graph[res, DirectedEdges -> False];

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

{1.075134, Null}