Solving the Covering Salesman Problem in Mathematica

Question

Say I have a weighted graph

vertexCount = 15;
edgeCount = 50;
ngroups = 4; 

vertexList = Range@vertexCount;
vertexGroups = ms[vertexList, ngroups]
vertexStyles = 
  Flatten@Map[Thread]@Thread[vertexGroups -> discreteColors[ngroups]];

SeedRandom[1];

edgeList = 
  Table[RandomInteger[{1, vertexCount}] <-> 
    RandomInteger[{1, vertexCount}], edgeCount];

    (* make simple graph *) 
edgeList = DeleteDuplicatesBy[Sort][edgeList]; 
edgeList = Select[edgeList, #[[1]] != #[[2]] &];

edgeWeights = 
  Table[RandomInteger[{1, vertexCount}], Length[edgeList]];

(* needs to be run twice *) 
graph = Graph[vertexList, edgeList, EdgeWeight -> edgeWeights, 
  VertexLabels -> Placed["Name", Center], 
  VertexStyle -> MapAt[Opacity[.25, #] &, {All, 2}]@vertexStyles, 
  VertexLabelStyle -> {v_ :> 
     Directive[16, Darker@Darker[v /. vertexStyles]]}, 
  VertexSize -> Large, ImageSize -> Large, ImagePadding -> 5, 
  Epilog -> blobF[graph, discreteColors[ngroups], vertexGroups, .07]]

I want to find the shortest path that transits all the groups - here, the coloured groups).

This is different from the TSP as

1. I'm allowed to pass over the same vertice twice (or more)
2. I'm travelling to groups (rather than all the vertices)

I can't see any option to do this with FindShortestPath. Would be interesting to use with FindGraphCommunities.

Functions required:

fC[pts_, size_ : .04] := Module[{}, CommunityGraphPlot[Graph@{}, {}];
  GraphComputation`GraphCommunitiesPlotDump`generateBlobs[
    Automatic, {pts}, size][[2, 1]]]
(* grouping thanks to @kglr *) 
blobF[g_, cols_, coms_, size_ : .04] := 
 Thread[{cols, EdgeForm[{Gray, Thin}], Opacity[.25], 
   fC[PropertyValue[{g, #}, VertexCoordinates] & /@ #, size] & /@ 
    coms}]
(* give discrete coulors *) 
discreteColors[n_] := 
 With[{partL = Ceiling[Sqrt[n]]}, 
  DeleteCases[
   Flatten[Transpose[
     Partition[
      Table[Lighter[Darker[Hue[c], .1], .25], {c, 0, 1 - 1/n, 1/n}], 
      partL, partL, 1, 0]]], 0]]
(* generate list *) 
ms = PadRight[Partition[#, UpTo[Ceiling[Length[#]/#2]]], #2, {{}}] &;

@kglr gave a solution, but only works for relatively small groups (the example given above will work, but very slow. Anything much bigger will give a memory error). I think this problem is reducible to TSP - but I'm not clear if this is possible using the inbuilt functions.

EDIT - this problem has a name - The Covering Salesman Problem.

Answer 1

After more thought, I managed to come up with a much better and faster solution than the previous one which used Optimization. Here we try to find the cheapest link between groupings until we cover all the groups.

For the @kglr's first example, it returns 8!

For the second example it correctly returns 6:

More importantly, on larger graphs, it performs much better while finding better answers. Here is the comparison of two algorithms on a graph with 300 vertices and 1000 edges:

First@FindMinimumCoveringSpanningTree[graph3, groupings3] // AbsoluteTiming
(* Out: {0.0201544, 158.6} *)

First@FindMinimumCoveringSpanningTreeOptimization[graph3, groupings3, 1] // AbsoluteTiming
(* Out: {0.244467, 297.181} *)

It accepts two arguments:

1. graph: should be directed (I think with some modification, it can support other types)
2. groupings: like {{1,2,3}, {3,4,5}} (groups can overlap)

Output: {total_weight, edges}

Code

ClearAll[FindMinimumCoveringSpanningTree, PositionSmallest2D];

(* Find the position of smallest element in matrix[group1,group2] \
subset *)
PositionSmallest2D[matrix_, group1_, group2_] := 
 Block[{result, temp, x, y},
  result = First@Ordering[Catenate@matrix[[group1, group2]], 1];
  temp = QuotientRemainder[result - 1, Length[group2]] + 1;
  
  x = group1[[First@temp]];
  y = group2[[Last@temp]];
  
  {x, y, matrix[[x, y]]}
  ]


FindMinimumCoveringSpanningTree[graph_Graph, 
  groupings_ : {{__Integer} ..}] := 
 Block[{coveredGroups, groupingsCount = Length[groupings], 
   output = {}, graphDistanceMatrix = GraphDistanceMatrix[graph], 
   range, groupingIndex, result, path},
  
  range = Range[groupingsCount];
  groupingIndex = 
   Association@Catenate[Thread /@ Thread[groupings -> range]];
  
  (* Find the cheapest link between two groups *)
  Block[{smallest = 
     First@TakeSmallestBy[
       Table[PositionSmallest2D[graphDistanceMatrix, groupings[[i]], 
         Catenate@Delete[groupings, i]], {i, groupingsCount}], Last, 
       1]}, path = 
    FindShortestPath[graph, smallest[[1]], smallest[[2]]];
   coveredGroups = DeleteDuplicates@Lookup[groupingIndex, path];
   output = UndirectedEdge @@@ Partition[path, 2, 1];
   ];
  
  (* Find cheapest link to connect more groups until all groups are \
covered *)
  Do[
   result = 
    PositionSmallest2D[graphDistanceMatrix, 
     DeleteDuplicates@VertexList[output], 
     Catenate@groupings[[Delete[range, List /@ coveredGroups]]]];
   path = FindShortestPath[graph, result[[1]], result[[2]]];
   output = Join[output, UndirectedEdge @@@ Partition[path, 2, 1]];
   coveredGroups = 
    DeleteDuplicates@
     Join[coveredGroups, Lookup[groupingIndex, path]];
   If[Length[coveredGroups] == groupingsCount, Break[]];
   , groupingsCount - 1];
  
  {Total@AnnotationValue[{graph, output}, EdgeWeight], output}
  ]

Answer 2

A brute-force approach:

1. Generate all connected subgraphs containing at least one vertex from each vertex group.
2. Pick one subgraph with the smallest total edge-weight.

ClearAll[connectedSubgraphs, shortestPathThroughGroups]

connectedSubgraphs[g_Graph, vgroups_] := Select[ConnectedGraphQ] @
  Map[Subgraph[SimpleGraph @ g, Flatten @ #] &]@
   Tuples@Map[Subsets[#, {1, ∞}] &] @ vgroups

shortestPathThroughGroups[g_Graph, vgroups_] := EdgeList @
  First @ MinimalBy[Total @ PropertyValue[#, EdgeWeight] &]@
    connectedSubgraphs[g, vgroups]

Example 1:

vertexCount = 9;
edgeCount = 30;

vertexList = Range @ vertexCount;
vertexGroups = Partition[vertexList, 3];
vertexStyles = Flatten @ Map[Thread] @ Thread[vertexGroups -> 
   {Blue, Orange, Green}];

SeedRandom[1];

edgeList = Table[RandomInteger[{1, vertexCount}] <-> 
    RandomInteger[{1, vertexCount}], edgeCount];

edgeWeights = Table[RandomInteger[{1, vertexCount}], edgeCount];

graph = Graph[vertexList, edgeList, 
  EdgeWeight -> edgeWeights, 
  VertexLabels -> Placed["Name", Center],
  VertexStyle -> MapAt[Opacity[.25, #] &, {All, 2}] @ vertexStyles, 
  VertexLabelStyle -> {v_ :> Directive[16, Darker@Darker[v /. vertexStyles]]}, 
  VertexSize -> Large,  ImageSize -> Large, ImagePadding -> 5,
  Epilog -> blobF[graph, {Blue, Orange, Green}, vertexGroups, .07]]

shortestPathThroughGroups[graph, vertexGroups]

 {4 <-> 1, 7 <-> 4}

HighlightGraph[SimpleGraph @ graph, 
 Style[#, Red, AbsoluteThickness[15]] & /@ 
  shortestPathThroughGroups[graph, vertexGroups],
 ImageSize -> Large, EdgeLabelStyle -> 16, 
 EdgeLabels -> {Alternatives @@ 
     shortestPathThroughGroups[graph, vertexGroups] -> "EdgeWeight"}]

Example 2:

SeedRandom[2345];
ew = RandomInteger[{1, 10}, 12];

vertexGroups = {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}};

g2 = GridGraph[{3, 3}, 
   EdgeWeight -> ew, 
   VertexLabels -> Placed["Name", Center],
   VertexStyle -> MapAt[Opacity[.25, #] &, {All, 2}]@vertexStyles, 
   VertexLabelStyle -> {v_ :> 
      Directive[16, Darker@Darker[v /. vertexStyles ]]}, 
   VertexSize -> Large,  ImageSize -> Large, 
   Epilog -> blobF[g2, {Blue, Orange, Green}, vertexGroups, .15], 
   ImagePadding -> 40, EdgeLabels -> "EdgeWeight", 
   EdgeLabelStyle -> 16]

shortestPathThroughGroups[g2, vertexGroups]

{1 <-> 4, 4 <-> 5, 5 <-> 6, 6 <-> 9}

HighlightGraph[g2, 
 Style[#, Red, AbsoluteThickness[15]] & /@ 
  shortestPathThroughGroups[g2, vertexGroups]]

Addendum: The blobs used as Epilog above are generated using the function blobF from this answer:

ClearAll[blobF, fC]
fC[pts_, size_ : .04] := Module[{}, CommunityGraphPlot[Graph@{}, {}];
  GraphComputation`GraphCommunitiesPlotDump`generateBlobs[
    Automatic, {pts}, size][[2, 1]]]

blobF[g_, cols_, coms_, size_ : .04] := 
 Thread[{cols, EdgeForm[{Gray, Thin}], Opacity[.25], 
   fC[PropertyValue[{g, #}, VertexCoordinates] & /@ #, size] & /@ coms}]

Answer 3

This solution is obsolete, Use this solution for better and faster results.

Instead of iterating over all possible subgraphs, I use linear optimization to find cheap, separated components covering all the groups that with some post-processing can be connected.

There are some limitations, as I couldn't formulate ConnectedGraphQ in a linear constraint. So the optimization finds separate components that are close to each other but not connected.

Since optimization find seperated paths, it ignores the cost of connecting them. In the post-processing, we'll find the cheapest path to connect these components but this does not guarantee the best solution.

It supports only undirected graphs, but I think with some modification, it can support other types too. You might find some corner cases where the input is not a connected graph.

Arguments:

1. Graph
2. grouping: {{1, 2, 3}, {4, 5, 6}, ...} (groups can overlap)
3. K: Does not have a clear definition but pushes the optimization to find components close to each other (sometimes you can get better results if you play with it).

Outpus: {cost, edges}

Testing the kglr's first example (SimpleGraph applied) clearly shows the limitation as the best answer costs 9 but ours costs 14 (k=1). As for the timing, it took us less than .007 seconds to find which is pretty fast.

On the kglr's second example we can get 6 for k=2, while k=3 gives 9.

We could test it on larger graphs like one with 300 vertices and 1000 edges and it costs us less than a second!

RandomSeed[10];
graph3 = 
  RandomGraph[{300, 1000}, EdgeWeight -> RandomReal[{1, 30}, 1000]];
groupings3 = Partition[RandomSample@Range[300], 10];

Result

AbsoluteTiming@
 First@FindMinimumCoveringSpanningTree[graph3, groupings3, 1]

(* Out: {0.306789, 367.896} *)

Next question is how much good are the answers? I think is highly dependant on the context and you need to test it on your data. If you found it don't forget to comment it below.

Code

I tried to pick informative names and named the main function FindMinimumCoveringSpanningTree! I haven't read papers about this algorithm, if nobody mentioned it before, notify me to change it to my personal name ;)

ClearAll[FindMinimumCoveringSpanningTree, FindEdgesToConnectEdges, 
  CreateSparseArrayFromIndexedList, 
  FindShortestPathFromDistanceMatrix];

(* Find the shortest path connecting one vertex from verticesSet1 to \
another vertex from verticesSet2 *)
FindShortestPathFromDistanceMatrix[distanceMatrix_, verticesSet1_, 
  verticesSet2_] := Block[{temp, minimumPath},
  
  temp = 
   Catenate@
    Outer[{#1, #2, distanceMatrix[[#1, #2]]} &, verticesSet1, 
     verticesSet2];
  
  minimumPath = Quiet@TakeSmallestBy[temp, Last, 1];
  
  If[minimumPath === {} || Head[minimumPath] === TakeSmallestBy, 
   Return@{}];
  
  minimumPath[[1, ;; 2]]
  ]

(* Tries to connect the edges by finding the minimum path,will return \
immediate if it could not connect the first component to the rest *)
FindEdgesToConnectEdges[graph_Graph, edges_] := 
 Block[{distanceMatrix = GraphDistanceMatrix[graph], addedEdges = {}, 
   components, otherComponenets, otherComponentsVertices, 
   mainComponentVertices, minimumPath},
  
  components = ConnectedGraphComponents[edges];
  
  Do[
   mainComponentVertices = VertexList[First@components];
   otherComponentsVertices = 
    DeleteDuplicates@Catenate[VertexList /@ Rest[components]];
   minimumPath = 
    FindShortestPathFromDistanceMatrix[distanceMatrix, 
     mainComponentVertices, otherComponentsVertices];
   
   If[minimumPath =!= {}, 
    minimumPath = 
     UndirectedEdge @@@ 
      Partition[
       FindShortestPath[graph, minimumPath[[1]], minimumPath[[2]]], 2,
        1];
    addedEdges = Join[addedEdges, minimumPath];
    components = ConnectedGraphComponents@Join[edges, addedEdges];, 
    Return[Join[edges, addedEdges]];
    ];
   
   , Length@components - 1
   ];
  
  Join[edges, addedEdges]
  ]

CreateSparseArrayFromIndexedList[lists_ : {{__Integer} ..}, 
  size_List] := 
 SparseArray[
  Catenate@Table[Thread[{i, lists[[i]]}], {i, Length[lists]}] -> 1, 
  size]

FindMinimumCoveringSpanningTree[graph_Graph, 
  groupings_ : {{__Integer} ..}, k_Integer?Positive] := 
 Block[{edgeList = EdgeList[graph], edgeCount = EdgeCount[graph], 
   edgeWeight = AnnotationValue[graph, EdgeWeight], 
   groupingsCount = Length[groupings], edgeListAsList, 
   kLevelEdgeConnection, groupingEdgeIndices, coverKLevel, x, 
   coverGroupings, edgeIndex, vertexEdgeIndexAssociation, 
   findEdgeIndexFromVertex, answer, result, coverKLevelReward},
  
  edgeIndex = PositionIndex@edgeList;
  edgeListAsList = List @@@ edgeList;
  vertexEdgeIndexAssociation = 
   Merge[{GroupBy[edgeList, First], GroupBy[edgeList, Last]}, 
    Catenate@Lookup[edgeIndex, Flatten@Join@#] &];
  
  (* Find all the edge indices based on a vertex *)
  findEdgeIndexFromVertex = 
   EdgeIndex[graph, 
     Cases[edgeList, 
      First[#] \[UndirectedEdge] _ | _ \[UndirectedEdge] First[#] | 
       Last[#] \[UndirectedEdge] _ | _ \[UndirectedEdge] Last[#]]] &;
  
  (* Find K-level edges connected to each edge *)
  kLevelEdgeConnection = 
   Function[x, 
     Nest[DeleteDuplicates@
        Catenate@Lookup[vertexEdgeIndexAssociation, #, Nothing] &, x, 
      k]] /@ (List @@@ edgeList);
  
  (* Find all the edge indices related to each grouping *)
  groupingEdgeIndices = findEdgeIndexFromVertex /@ groupings;
  
  (* Optimization constraint-If an edge is selected,edges to K-
  Level should also be picked *)
  coverKLevel = 
   CreateSparseArrayFromIndexedList[
    kLevelEdgeConnection, {edgeCount, edgeCount}];
  
  (* Optimization constraint-If an edge is selected,statisfy the K-
  th Level requirement *)
  coverKLevelReward = 
   Function[x, 
     Rest@NestList[
       DeleteDuplicates@
         Complement[
          Catenate@
           Lookup[vertexEdgeIndexAssociation, #, Nothing], #] &, x, 
       k]] /@ (List @@@ edgeList);
  
  coverKLevelReward = -Transpose@
       CreateSparseArrayFromIndexedList[#, {edgeCount, edgeCount}] & /@
     Transpose[coverKLevelReward];
  
  (* Optimization constraint-
  selected edges should cover all the groupings *)
  coverGroupings = ConstantArray[0, {groupingsCount, edgeCount}];
  Do[coverGroupings[[i, groupingEdgeIndices[[i]]]] = 1;, {i, 
    groupingsCount}];
  
  answer = 
   LinearOptimization[
    x . edgeWeight, {1 \[VectorGreaterEqual] x \[VectorGreaterEqual] 
      0, (x . coverKLevel) + 
       Total[coverKLevelReward . x] \[VectorLessEqual] 0, 
     coverGroupings . x \[VectorGreaterEqual] 1}, 
    x \[Element] Vectors[edgeCount, Integers], 
    "PrimalMinimizerVector"];
  
  If[FreeQ[answer, Indeterminate], 
   result = 
    FindEdgesToConnectEdges[graph, 
     edgeList[[Catenate@Position[answer, 1]]]];
   
   (* The result may have mutliple edges betwen nodes,
   using FindSpanningTree will fixed,
   subgraph was used to get EdgeWeight *)Block[{tempGraph = result},
    {Total@AnnotationValue[{graph, EdgeList@tempGraph}, EdgeWeight], 
     tempGraph}
    ]
   
   , {0, {}}
   ]
  ]