8
$\begingroup$

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];
$\endgroup$
1
  • 1
    $\begingroup$ 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! $\endgroup$
    – dakta
    Jan 28, 2014 at 11:28

3 Answers 3

9
$\begingroup$
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.

$\endgroup$
7
  • $\begingroup$ 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] $\endgroup$ Jan 27, 2014 at 19:58
  • $\begingroup$ @belisarius It seems I'm particularly stupid today. Does it look fine now? $\endgroup$
    – Szabolcs
    Jan 27, 2014 at 20:35
  • $\begingroup$ Yep :) +1 Now :) $\endgroup$ Jan 27, 2014 at 20:46
  • 1
    $\begingroup$ @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. $\endgroup$
    – Szabolcs
    Jan 27, 2014 at 22:03
  • 1
    $\begingroup$ @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? $\endgroup$
    – dakta
    Jan 27, 2014 at 23:51
7
$\begingroup$

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}
*)
$\endgroup$
3
  • $\begingroup$ 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? $\endgroup$
    – dakta
    Jan 27, 2014 at 17:14
  • $\begingroup$ @user11977 Docs are here reference.wolfram.com/mathematica/tutorial/UsingNearest.html $\endgroup$ Jan 27, 2014 at 17:26
  • $\begingroup$ Fair enough. The documentation of NearestFunction does not mention this calling sequence leading me to assume it is undocumented ... Anyway, thank you very much. $\endgroup$
    – dakta
    Jan 27, 2014 at 19:09
4
$\begingroup$

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}

$\endgroup$
1
  • $\begingroup$ Thank you very much! I knew there had to be something better than GraphDistanceMatrix. AdjacencyList escaped me ... $\endgroup$
    – dakta
    Jan 27, 2014 at 21:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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