# 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[[#[]]] - pos[[#[]]]] < 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}], (#[] < #[])&] /.
Thread[Range[Length[v]] -> v];

• 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

## 3 Answers

n = 20000;(*desired:n=20000*)SeedRandom;
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.

• 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*)
SeedRandom;
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}
*)

• 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 reference.wolfram.com/mathematica/tutorial/UsingNearest.html – 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;
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}

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