Given a directed graph, how can I pick out pairs of points that have 1 direct path and no other paths with intervening points?

Question

I have a directed graph generated from pairs of points as follows

n = RandomVariate[PoissonDistribution[100]];
pt=RandomReal[{-0.5, 0.5}, {n, 2}].RotationMatrix[45 Degree];
pt = Sort[pt, #1[[2]] < #2[[2]] &];
c = Table[If[(pt[[i, 2]] - pt[[j, 2]])^2 - (pt[[i, 1]] - pt[[j, 1]])^2 > 0, 1, 0], {i, n}, {j, n}];
c = UpperTriangularize[c];
g = AdjacencyGraph[c, VertexLabels -> "Name"];

I want to now obtain a new graph/matrix that only has ones when 2 points have a single direct path between them and no other longer paths (i.e. paths through other intermediary points). I was using the following code for this

l = Table[If[c[[i, j]] == 1 && Length[Flatten[FindPath[g, i, j, n]]] == 2, 1, 0], {i, n}, {j,n}];

but this line of code takes a very long time to run for larger sizes (350 and above). Is there a more efficient (and more importantly faster) way of doing this?

Answer 1

Since c upper triangular, it is a nilpotent matrix. Its degree is less than 20 (by trial and error):

Max@MatrixPower[c, 20]
(* 0 *)

The number of paths between two vertices is then given by

np = Sum[MatrixPower[c, k], {k, 20}];

To get the adjacency matrix of vertex pairs which are connected by a single path only, we can use

1 - Unitize[np - 1]

---

Update: Based on Carl Woll's suggestion, we can also do

result = 1 - Unitize[Total@FixedPointList[#.c &, c] - 1];

Answer 2

Update: TransitiveReductionGraph is much faster:

t1 = First[RepeatedTiming[g1 = TransitiveReductionGraph[g];]]
t2 = First[RepeatedTiming[g2 = Graph@Select[EdgeList[g], 
      Length[FindVertexIndependentPaths[g, ## & @@ #, Infinity]] == 1 &];]]
t3 = First[RepeatedTiming[g3 = AdjacencyGraph[1 -
      Unitize[Sum[MatrixPower[c, k], {k, 20}] - 1]];]]
t4 = First[RepeatedTiming[g4 = prunedGraph[g];]]
t5 = First[RepeatedTiming[g5 = AdjacencyGraph[1 - Unitize[pathCounts[c] - 1]];]]
t0 = First[RepeatedTiming[g0=AdjacencyGraph@Table[If[c[[i, j]] == 1 && 
   Length[Flatten[FindPath[g, i, j, n]]] == 2, 1, 0], {i, n}, {j, n}];]]

Equal @@ (EdgeList/@ {g0,g1, g2, g3, g4, g5})

True

Grid[Prepend[Transpose[{{"FindPath","TransitiveReductionGraph", 
  "FindVertexIndependentPaths", "MatrixPower", "prunedGraph", "pathCounts"}, 
  {t0, t1, t2, t3, t4, t5}}], {"method", "RepeatedTiming"}], Dividers -> All] // TeXForm

$\begin{array}{|c|c|} \hline \text{method} & \text{RepeatedTiming} \\ \hline \text{FindPath} & 2.211 \\ \hline \text{TransitiveReductionGraph} & 0.0023 \\ \hline \text{FindVertexIndependentPaths} & 0.34 \\ \hline \text{MatrixPower} & 0.218 \\ \hline \text{prunedGraph} & 0.128 \\ \hline \text{pathCounts} & 0.0044 \\ \hline \end{array} $

Although all five methods posted so far give the same result, as noted by Szabolcs in a comment, TransitiveReductionGraph has some yet-unfixed bugs. Carl's modification of Szabolc's approach is the fastest among remaining four methods posted so far.

Original answer:

You can use FindVertexIndependentPaths combined with Select:

Select[EdgeList[g], Length[FindVertexIndependentPaths[g, ##& @@ #, ∞]] == 1&] // Length // 
 AbsoluteTiming

{0.474247, 358}

 Select[EdgeList[g], Length[FindVertexIndependentPaths[g, ##& @@ #, ∞]] == 1&] // Short

{1 -> 5, 1 -> 8, 1 -> 9, 1 -> 12, 1 -> 13, 1 -> 14, <<347>>, 101 -> 105, 102 -> 103, 103 -> 104, 103 -> 105, 105 -> 106}

Answer 3

(See addendum below for a refinement of the approach suggested by @Szabolcs)

Another idea is to delete an edge, and check if the graph distance between the vertices is still finite. This won't be as fast as using TransitiveReductionGraph as suggested by @kglr. Here is a function to do this:

prunedGraph[g_] := With[{el = EdgeList[g]},
    Graph @ Pick[el, GraphDistance[EdgeDelete[g, #], Sequence @@ #]& /@ el, Infinity]
]

Using this function on your graph:

p1 = prunedGraph[g]; //AbsoluteTiming
p2 = TransitiveReductionGraph[g]; //AbsoluteTiming

IsomorphicGraphQ[p1, p2]

{0.117671, Null}

{0.002606, Null}

True

Addendum

It is possible to speed up the suggested by @Szabolcs. By using nesting instead of MatrixPower, and making sure to work with packed arrays, we have:

pathCounts[am_?MatrixQ] := With[{s = Developer`ToPackedArray@am},
    FixedPoint[s + # . s &, s, Length[am]]
]

This approach has the advantage that we don't need to know how many terms to Sum, we just keep applying the function until the matrix doesn't change, or until we've done it enough times to know that it won't change. Compare to the approach using MatrixPower:

r1 = pathCounts[c]; //AbsoluteTiming
r2 = Sum[MatrixPower[c, k], {k, 20}]; //AbsoluteTiming

r1 == r2

{0.001114, Null}

{0.230795, Null}

True

Using pathCounts, we can generate your new graph using:

r1 = AdjacencyGraph[1 - Unitize[pathCounts[c] - 1]]; //RepeatedTiming
r2 = TransitiveReductionGraph[g]; //RepeatedTiming

r1 === r2

{0.0021, Null}

{0.0017, Null}

True