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

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

n = RandomVariate[PoissonDistribution];
pt=RandomReal[{-0.5, 0.5}, {n, 2}].RotationMatrix[45 Degree];
pt = Sort[pt, #1[] < #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?

## 3 Answers

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

• You can make this approach much faster by using SparseArrays, and using nesting (e.g., FixedPoint[ c + #.c &, c, Length[c]]) instead. – Carl Woll Aug 23 '17 at 17:18

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}

• I am not comfortable using TransitiveReductionGraph because of this bug which has been present since the introduction of the function, reported multiple times, yet still not fixed in 11.1 ... (Otherwise great idea!) – Szabolcs Aug 23 '17 at 16:26
• Thank you @Szabolcs, good point. – kglr Aug 23 '17 at 16:30
• @kglr can you tell me how FindVertexIndependentPaths is different from FindPath? the documentation isn't very clear on this. Also your code seems to give a different answer from the one I originally had...I'm not sure why – nomaan x Aug 23 '17 at 17:08
• @kglr The edge list output you give is different than what I get. Could you check whether you have the right output? Also, you might want to compare my refinement of Szabolcs' answer as well. – Carl Woll Aug 23 '17 at 18:57
• @CarlWoll, i added pathCounts to the timing comparisons. Re different edge list output, it is probably because i did not set SeedRandom in my original post. In the updated version I used SeedRandom. – kglr Aug 23 '17 at 19:46

(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 = DeveloperToPackedArray@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